Truncation Error Analysis for DiffusionSchemes in Boundary Regions ofUnstructured MeshesbyVarun Prakash PuneriaB.Tech., Aerospace Engineering, Indian Institute of Technology Kharagpur (India), 2010A 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)April 2015c© Varun Prakash Puneria 2015AbstractAccuracy of numerical solution is of paramount importance for any CFD simula-tion. The error in satisfying the continuous partial differential equations by theirdiscrete form results in truncation error and it has a direct influence on discretiza-tion errors. Discretization error, which is the difference in the numerical solutionand the exact solution to a CFD problem, is generally the largest source of numeri-cal errors. Understanding the relationship between the discretization and truncationerror is crucial for reducing numerical errors. Studies have been carried out to un-derstand the truncation-discretization error relationship in the interior regions of acomputational domain but fewer for the boundary regions. The effect of differentsolution reconstruction methods, face gradient averaging schemes and boundarycondition implementation methods on boundary truncation error in specific andoverall truncation and solution error are the subject of research for this thesis. Toachieve the goals laid out, the error has to be quantified first and then tests per-formed to compare the schemes. The Poisson’s equation is chosen as the modeldiffusion equation. The truncation error coefficients, analogous to the analyticalcoefficients of the spatial derivatives in Taylor series expansion of truncation error,are quantified using error metrics. The solution error calculation is made possi-ble by a careful selection of exact solutions and their appropriate source terms forPoisson’s equation. Analytic tests are performed on a family of topologically reg-ular meshes to test and verify the theoretical implementation of different schemesand to eliminate schemes performing poorly from consideration for numerical tests.The numerical tests are performed on unstructured triangular, mixed and pure quadmeshes to extend the accuracy assessment for general meshes. The results obtainedfrom both the tests are utilized to arrive at schemes where the overall truncation er-ror and discretization error can be minimized simultaneously.iiPrefaceThe research ideas and methods discussed in this thesis are the fruits of a closeworking relationship between Dr. Carl Ollivier-Gooch and Varun Prakash Puneria.The implementation of methods, the data analysis and the manuscript preparationwere done by Varun Prakash Puneria with invaluable guidance from Carl Ollivier-Gooch throughout the process.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1 Finite Volume Approach . . . . . . . . . . . . . . . . . . . . . . 152.2 Solution Reconstruction . . . . . . . . . . . . . . . . . . . . . . 162.2.1 Green-Gauss Method . . . . . . . . . . . . . . . . . . . 172.2.2 Least-Squares Method . . . . . . . . . . . . . . . . . . . 182.2.2.1 Framework for Least-Squares Method . . . . . 192.2.2.2 Weighted Least-Squares Method . . . . . . . . 222.2.2.3 Advantages of Least-Squares Method . . . . . 232.2.3 Reconstruction Stencils . . . . . . . . . . . . . . . . . . 242.3 Flux Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.3.1 Numerical Flux Calculation . . . . . . . . . . . . . . . . 262.3.2 Flux Integration . . . . . . . . . . . . . . . . . . . . . . 282.4 Boundary Treatment Methods . . . . . . . . . . . . . . . . . . . 292.4.1 Use of Hard Constraints . . . . . . . . . . . . . . . . . . 29ivTable of Contents2.4.1.1 Strong Boundary Condition . . . . . . . . . . . 302.4.1.2 4-Unknown Model . . . . . . . . . . . . . . . 302.4.2 Increased Stencil Size . . . . . . . . . . . . . . . . . . . 312.4.3 Use of Single and Two Jump Terms . . . . . . . . . . . 322.4.4 Ghost Cell Methods . . . . . . . . . . . . . . . . . . . . 333 Methods of Error Analysis . . . . . . . . . . . . . . . . . . . . . . . 353.1 Analytic Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.2 Numerical Tests . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2.1 Truncation Error . . . . . . . . . . . . . . . . . . . . . . 433.2.2 Discretization Error . . . . . . . . . . . . . . . . . . . . 454 Triangular Mesh Results . . . . . . . . . . . . . . . . . . . . . . . . 474.1 Analytic Test Analysis . . . . . . . . . . . . . . . . . . . . . . . 474.1.1 Reconstruction Accuracy . . . . . . . . . . . . . . . . . 484.1.2 Flux Integral Accuracy . . . . . . . . . . . . . . . . . . . 534.1.3 Boundary Treatment Methods . . . . . . . . . . . . . . . 554.2 Numerical Test Analysis . . . . . . . . . . . . . . . . . . . . . . 594.2.1 Isotropic Boundary Layer Mesh . . . . . . . . . . . . . . 594.2.2 Stretched Boundary Layer Mesh . . . . . . . . . . . . . 755 Mixed Mesh Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.1 Analytic Test Analysis . . . . . . . . . . . . . . . . . . . . . . . 815.1.1 Reconstruction Accuracy . . . . . . . . . . . . . . . . . 825.1.2 Flux Integral Accuracy . . . . . . . . . . . . . . . . . . . 835.1.3 Boundary Treatment Methods . . . . . . . . . . . . . . . 875.2 Numerical Test Analysis . . . . . . . . . . . . . . . . . . . . . . 895.2.1 Isotropic Boundary Layer Mesh . . . . . . . . . . . . . . 915.2.2 Stretched Boundary Layer Mesh . . . . . . . . . . . . . 1136 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 1206.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1206.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1216.3 Recommended Future Work . . . . . . . . . . . . . . . . . . . . 124Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126vList of Tables2.1 Face gradient calculation . . . . . . . . . . . . . . . . . . . . . . 273.1 Coefficients for boundary numerical flux integral . . . . . . . . . 414.1 Gradient error coefficients for boundary control volumes . . . . . 514.2 Gradient error coefficients for interior control volumes . . . . . . 514.3 Gradient error coefficients using corrected Green-Gauss method . 534.4 Truncation error order comparison for flux integral . . . . . . . . 544.5 Manufactured solutions used for numerical tests . . . . . . . . . . 614.6 Optimal jump coefficients for minimum truncation and discretiza-tion error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.7 Optimal jump coefficients for minimum truncation and discretiza-tion error, using case 1 exact solution . . . . . . . . . . . . . . . . 805.1 Gradient error coefficients for boundary control volumes . . . . . 845.2 Gradient error coefficients for interior control volumes . . . . . . 845.3 Gradient error coefficients for interface control volumes . . . . . . 855.4 Truncation error order comparison for flux integral . . . . . . . . 865.5 Optimal jump coefficients for minimum truncation and discretiza-tion error, using case 1 exact solution . . . . . . . . . . . . . . . . 106viList of Figures1.1 Flow chart for standard procedures in CFD . . . . . . . . . . . . . 21.2 Description of mesh types . . . . . . . . . . . . . . . . . . . . . 41.3 One dimensional grid for convection-diffusion equation . . . . . . 92.1 Control volume illustration . . . . . . . . . . . . . . . . . . . . . 162.2 Gradient reconstruction for mixed cell interface using Green-Gaussmethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3 Stencil choices for reconstruction . . . . . . . . . . . . . . . . . . 252.4 Linear interpolation method for numerical flux . . . . . . . . . . . 272.5 Numerical flux calculation for Green-Gauss method . . . . . . . . 282.6 Illustration of Gauss quadrature points for flux integration . . . . . 292.7 Illustration of stencil increment . . . . . . . . . . . . . . . . . . . 322.8 Illustration for jump term coefficients . . . . . . . . . . . . . . . 332.9 Illustration of ghost cell method . . . . . . . . . . . . . . . . . . 343.1 Uniform mesh regular stencil . . . . . . . . . . . . . . . . . . . . 363.2 Uniform mesh boundary stencil . . . . . . . . . . . . . . . . . . . 393.3 Uniform mixed cell mesh regular stencil . . . . . . . . . . . . . . 404.1 Family of distorted meshes . . . . . . . . . . . . . . . . . . . . . 494.2 Gradient error metrics for different schemes . . . . . . . . . . . . 564.3 Comparison of boundary treatment methods . . . . . . . . . . . . 584.4 Isotropic triangular mesh cases . . . . . . . . . . . . . . . . . . . 594.5 Truncation error distribution for Sq30 mesh using least-squaresmethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.6 Exact solutions used for discretization error analysis . . . . . . . . 624.7 Convergence plots for discretization errors . . . . . . . . . . . . . 634.8 Plot for truncation error using single jump term method . . . . . . 644.9 Case 1: Solution error plot using single jump term method . . . . 654.10 Case 2: Solution error plot using single jump term method . . . . 664.11 Case 3: Solution error plot using single jump term method . . . . 664.12 Comparison of truncation error for different schemes . . . . . . . 67viiList of Figures4.13 Comparison of schemes for discretization error . . . . . . . . . . 684.14 Pareto plot for discretization and truncation errors correspondingto case-1 exact solution, Sq30 mesh . . . . . . . . . . . . . . . . 714.15 Pareto case error deviation from minimum possible errors . . . . 714.16 Comparison of truncation errors for different minimization cases,Sq30 mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.17 Comparison of discretization errors for different minimization cases,Sq30 mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.18 Truncation error contour plot variation with two jump terms, Sq30mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.19 Case 1: Discretization error contour plot using two jump termmethod, Sq30 mesh . . . . . . . . . . . . . . . . . . . . . . . . . 744.20 Case 2: Discretization error contour plot using two jump termmethod, Sq60 mesh . . . . . . . . . . . . . . . . . . . . . . . . . 744.21 Case 3: Discretization error contour plot using two jump termmethod, Sq120 mesh . . . . . . . . . . . . . . . . . . . . . . . . 754.22 Stretched boundary layered triangular mesh . . . . . . . . . . . . 764.23 Truncation error distribution for Sq30 stretched mesh using least-squares method . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.24 Truncation error variation using single jump term for triangularmeshes with stretched boundary layers . . . . . . . . . . . . . . . 784.25 Comparison of different schemes with regards to truncation anddiscretization error . . . . . . . . . . . . . . . . . . . . . . . . . 794.26 Pareto case error deviation from minimum possible errors . . . . . 805.1 Gradient error metrics for different schemes . . . . . . . . . . . . 875.2 Comparison of boundary treatment methods . . . . . . . . . . . . 895.3 Comparison of single and two jump term methods on randomlystretched meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 905.4 Comparison of interior and interface truncation errors . . . . . . . 905.5 Isotropic meshes for numerical analysis . . . . . . . . . . . . . . 925.6 Truncation error distribution for Sq30 mesh with two isotropic quadlayers, using LSM . . . . . . . . . . . . . . . . . . . . . . . . . 935.7 Truncation error distribution for Sq30 mesh with five isotropic quadlayers, using LSM . . . . . . . . . . . . . . . . . . . . . . . . . 935.8 Convergence of solution using least-squares and Green-Gauss meth-ods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.9 Interior truncation error using single jump term . . . . . . . . . . 965.10 Boundary truncation error using single jump term . . . . . . . . . 97viiiList of Figures5.11 Plot for solution error on square meshes with five boundary quadlayers, using the single jump term method . . . . . . . . . . . . . 995.12 Plot for solution error on pure quad square meshes, using the singlejump term method . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.13 Comparison for truncation error on square meshes with two quadlayers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1015.14 Comparison for truncation error on square meshes with five quadlayers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.15 Comparison for truncation error on pure quad meshes . . . . . . . 1035.16 Comparison of solution using least-squares and Green-Gauss meth-ods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.17 Pareto plot for discretization and boundary truncation error corre-sponding to case-1 exact solution . . . . . . . . . . . . . . . . . . 1075.18 Pareto case error deviation from minimum possible errors . . . . . 1095.19 Variation of boundary truncation error in two-jump term space . . 1105.20 Variation of discretization error in two-jump term space for Sq30mesh with five Quad layers . . . . . . . . . . . . . . . . . . . . . 1115.21 Variation of discretization error in two-jump term space for a 32x32pure quad mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.22 Stretched quad layer mixed meshes . . . . . . . . . . . . . . . . . 1135.23 Truncation error distribution for Sq30 mesh with five stretchedquad layers using LSM . . . . . . . . . . . . . . . . . . . . . . . 1145.24 Boundary truncation error variation using single jump term method 1155.25 Comparison of boundary truncation error for different schemes onmixed mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1175.26 Comparison of solution error for different schemes on mixed mesh 1185.27 Pareto case error deviation from minimum possible errors . . . . . 119ixList of SymbolsRoman Symbolsnˆ unit normal vectore edge vectorL continuous linear differential operatorfB continuous or discrete boundary condition functionA Areaa curved mesh aspect ratioCV control volumeEX, EY gradient truncation error metricE flux integral truncation error metricg boundary gauss pointh characteristic mesh length scalek mesh scaling factorm row number for analytic meshesN number of stencil membersR residual, curved mesh radius ratios length of control volume edges, mesh stretching factorV control volume area for integrationw weights for least-squares methodx, y Cartesian coordinatesxList of Symbolsr position vectorGreek Symbolsα interior jump term coefficientαB boundary jump term coefficientε discretization errorλ flux integral truncation error coefficientsν diffusion coefficientφ control volume average of the variableΦ unknown variables vectorφ solution variable, exact solution valueρ random perturbationsτ truncation errorφ˜ discrete form of the variableξ , η cell gradient truncation error coefficientsSuperscriptsB boundary condition valueL left sidep, q polynomial exponentR right side, reconstructed functionSubscriptsbdry boundaryF faceh discretized form of the continuous equation or variableint interiorxiList of Symbolsi reference control volumex,y Cartesian coordinatesj neighboring control volumexiiAcknowledgementsFirst and foremost, I would like to acknowledge my research supervisor, Dr. CarlOllivier-Gooch, whose constant support, patience and guidance has been invalu-able throughout the course of this research. Carl is a wonderful person and a greatresearch adviser.I am also grateful to my colleagues Alireza Jalali and Gary Yan at ANSLabresearch group for the suggestions and discussions which have been instrumentalfor my research. I would also like to thank David Zuniga from ANSLab researchgroup for his invaluable time and the contributions to the mesh solver GRUMMPwhich made the test cases possible for the research.Finally, my deepest appreciation is for my family, especially my mother “Durga”and father “Om”, and friends for their constant support and for their endless loveand encouragement throughout my life.xiiiChapter 1IntroductionKnowledge of fluid dynamics was not uncommon among the ancient civilizations.Design of arrows, spears, boats, irrigation systems, etc. demonstrate their knowl-edge in fluid flows. Notable contributions were made by Archimedes who formu-lated the laws of buoyancy and applied them to floating and submerged bodies.Leonardo da Vinci derived the equation of conservation of mass in one dimen-sional steady flow and his notes contained accurate descriptions of waves, jets,hydraulic pumps, eddy formation and both low and high drag designs. Sir IsaacNewton noted the effects of viscosity in diminishing the velocity of running waterand his second law of motion was instrumental to the development of equationsof fluid motion by Leonhard Euler. Euler’s idea to express the motion of fluids inpartial differential equations was a major breakthrough, although it was limited inthe sense that it did not account for the viscous forces.The individual efforts of George Stokes and Claude Navier resulted in a com-plete set of governing equations for compressible viscous flows. However, theequations were too complicated to develop analytical solutions for arbitrary flows.This led to simplified forms of the equations with analytical solutions that appliedto special cases of fluid flows, which were laminar or inviscid in nature. Thesespecial cases of flows had limited applications since most flows of engineeringinterests are viscous and turbulent.The need for experiments has been indispensable to the advancement of under-standing fluid dynamics. With limited analytical solutions, experimental fluid dy-namics provided the much needed understanding of arbitrary fluid flows for com-plicated geometries. However, the costs that are necessary for both the experi-mental facility and the physical model to adequately simulate the prototype flowfield restricts their application to modern day problems. Nevertheless, the data ob-tained from experiments is important for validating the mathematical solutions ofthe governing equations.Up until twentieth century, theoretical and experimental fluid dynamics werethe two approaches for the study of fluids. From the 1960s onward, the advent ofcomputers and development of numerical algorithms has led to a third approachcalled Computational Fluid Dynamics (CFD). It has been an equal partner to theanalytical and experimental approaches. Unlike experimental fluid dynamics, the1Chapter 1. IntroductionFigure 1.1: Flow chart for standard procedures in CFDflow conditions and geometry of the problem can be easily varied to optimize thedesign process. Impractical situations like the re-entry of space vehicles can besimulated through CFD whereas analytical and experimental approaches are inad-equate in present capacities.The procedure of solving a fluid flow problem using CFD is broadly catego-rized into three parts [44, 41] : 1) Pre-processing which includes preparation of thecomputational domain and its discretization by meshing, 2) Numerical solver in-corporates the mathematical expression of the governing laws along with initial andboundary conditions of the fluid flow problem to compute the numerical solution,and 3) Post-processing which is the analysis and visualization of the data generatedfrom the solver. Figure 1.1 illustrates the schematic flow chart for procedures inCFD.2Chapter 1. IntroductionPre-processing constitutes a major part of the efforts realized for any CFDproject. The extraction of a computational domain after simplification of the phys-ical geometry forms the first step for any CFD project. The computational domainshould then be decomposed into smaller non-overlapping subdomains through theprocess of mesh generation. The accuracy of a numerical simulation is governedby the number of cells used to divide the domain. The larger the number of cellsthe better the resolution of the flow and the better the accuracy of the solution.There are two types of mesh based on their topology, namely the structuredand unstructured mesh types. In the structured mesh type, each of the grid el-ement is identified by the indices (i, j,k) which corresponds to their position inspace. The regular connectivity of the mesh implicitly gives the address of theneighboring control volumes which makes it easier for performing calculations inthe solver. The evaluation of fluxes, gradients and treatment of boundary condi-tions is greatly simplified by using structured grids. Structured grids commonlyemploy quadrilateral elements for two dimensions and hexahedral elements forthree dimensions. Unstructured meshes lack the regular connectivity as observedin structured meshes. This means that the solver needs to process the mesh anddetermine the connectivity data with regards to the vertices, edges and neighboringcells. In two dimensions, the unstructured mesh is composed of triangular ele-ments whereas for three dimensional meshes tetrahedral, prismatic and polyhedralelements can be used.Compared to unstructured meshes, generation of structured meshes for sim-pler geometries is easier. The associated memory requirements, numerical errorand computational time for structured meshes are less compared to unstructuredmeshes. However, for complex geometries, generating a structured mesh can bevery time consuming or infeasible. Unstructured meshes offer the greatest flex-ibility in the treatment of complex geometries. The main advantage of unstruc-tured meshes is the fact that the mesh generation can be performed automaticallywith minimal inputs, independent of the complexity of the domain. Unstructuredmeshes comprised of different element types are referred to as mixed grids whereasgrids which have combined structured-unstructured grids are termed as hybridgrids. Figure 1.2 presents the different types of mesh that can be constructed arounda 2D airfoil.A popular method to accurately resolve the boundary layers for viscous flowsis to use quadrilateral elements in two dimensions and prismatic or hexahedralelements in three dimensions near the solid walls [4]. Usually, viscous flow regionsrequire high aspect ratio grid cells and the ways to achieving this condition is bygenerating 90 degree triangles or rectangles in shear layers [7, 38, 42]. In threedimensions, this can be achieved by using prismatic or hexahedral cells near theboundary.3Chapter 1. Introduction(a) Structured(b) Unstructured(c) MixedFigure 1.2: Description of mesh types4Chapter 1. IntroductionThe development of flow simulation software begins with the mathematicalmodeling of physical phenomena. It is here that the decision to omit or considerdifferent aspects of the flow to be modeled are made. The extent to which thephysical phenomena is modeled will determine the extent to which the numericalsimulation will represent the actual flow scenarios. For instance, whether the flowis compressible or incompressible, inviscid or viscous, laminar or turbulent willdetermine the set of mathematical equations, typically partial differential equa-tions (PDEs), to be solved. Even if we chose to solve the flow of certain nature,say incompressible turbulent flows, the extent to which the flow features need tobe resolved will determine the type of turbulence models needed. For example,the Spalart-Allmaras [39] turbulence model is a good RANS model for aerody-namic flows but computer intensive models like Large Eddy Simulation (LES) orDirect Numerical Simulation (DNS) are essential for complex or separated turbu-lent flows. The requirements to model the physical phenomena are not limited toturbulence but may also include multiphase interactions, chemical reactions andelectromagnetic phenomena. The tremendous requirements on computational re-sources to solve complex phenomena like turbulence or chemical reactions neces-sitates approximation by physical modeling and the errors associated with it. Thestudy of these errors, called physical modeling errors, is not within the scope ofthis research.The initial values of flow variables must be specified at all solution points inthe domain. The initial conditions determine the state of the fluid at time t = 0 orat the first step of an iterative scheme. A good initial guess implies that conver-gence to a steady-state solution will be faster. Boundary conditions along with thecomputational domain give the information on the physical problem that is beingsolved for. For instance, if the computational domain is that of an airfoil with anappropriate mesh, the boundary conditions will dictate if we are solving a subsonicor a supersonic flow. Some of the common boundary conditions are the wall, in-let and outlet, symmetry, periodic, etc. Many common boundary conditions areclassified as Dirichlet or Neumann boundary conditions. If the variable is specifiedalong the boundary, it is known as the Dirichlet type whereas if the normal gradientof the variable is specified, it is called the Neumann type boundary condition. Alinear combination of the Dirichlet and Neumann types can also be specified andis known as the Robin type boundary condition.After the governing equations are chosen and initial and boundary conditionsare set, the PDEs need to be discretized so that they can be solved on a computer.There are three popular methods of discretization for the PDEs : finite differencemethods, finite element methods and finite volume methods.The finite difference method was the first method applied to calculate numeri-cal solution of differential equations. The solution variables Φ are stored as point5Chapter 1. Introductionvalues at the grid points. Using Taylor series expansions, the derivatives of thevariables Φ are approximated by their finite difference equations in terms of theirpoint wise solution values. These finite difference approximations are substitutedfor the derivatives of the variables Φ in the governing equations at every grid point,resulting in a set of algebraic equations in terms of point wise solution variablesat the grid points. Solving these algebraic expressions for the entire domain givesus the solution data. The finite difference method is very simple, effective andeasy to obtain higher order schemes on regular grids. However, the conservation ofmass, momentum and energy is not enforced unless special care is taken and thisapproach is impractical for unstructured grids.The finite element method was initially developed for structural stress analysisand is also applied for fluid dynamics. This method uses simple piecewise ba-sis functions, which are valid on elements and describe the local variations of theunknown variables Φ. The choice of basis functions determines the accuracy ofthe scheme. This solution approximation is substituted into the governing equa-tions; and since it does not satisfy the equations exactly, a residual is obtained. Theresiduals are minimized in some sense by multiplying them with weighing func-tions and integrating. The resulting set of algebraic equations are solved for theunknown coefficients of the approximating function. This method is easily extend-able to unstructured meshes and higher order schemes. However, the conservationof mass, momentum and energy, though possible, is not trivial in flows with dis-continuities.The finite volume method uses integration of the governing equations of fluidflow over each control volume of the domain. Applying Gauss’s theorem leadsto an equation relating the amount of a conserved quantity in a control volume tothe flux of that quantity through the control volume boundary. The fluxes, in turn,are approximated in terms of control volume averages, which results in a systemof algebraic equations. The solution field is obtained by solving this system ofequations. It is the easiest method to understand and the approximated terms havephysical meanings which makes it a popular method of discretization for CFD.Mass, momentum and energy are directly conserved using this method and it canbe easily implemented on structured and unstructured grids.Irrespective of the type of discretization used, the approximation leads to anerror due to the difference in the discrete and continuous equations called the trun-cation error. The error in satisfying the exact solution to the continuous PDEsby the solution to the discretized equations is called the discretization error. Dis-cretization error is the largest source of numerical errors (others being round-offerror, iterative convergence error, etc). Truncation error acts as a local source forthe transport of discretization error [35] and it is usually estimated by replacing allthe nodal values in the discrete approximation by a Taylor series expansion about a61.1. Motivationsingle point (see page 9). The convergence of truncation errors is generally a goodindicator of convergence of discretization errors on structured meshes; however,the degradation of convergence for truncation errors doesn’t necessarily imply adegradation in discretization errors for unstructured meshes. Nevertheless, manymethods use the truncation error for mesh adaptation to reduce the discretizationerrors [35, 46]. By reducing the truncation errors, it is possible to reduce the dis-cretization errors.The finite volume method of discretization is used for the research presentedin this thesis. A good CFD solver incorporates discretization techniques suitablefor the treatment of the key physical phenomena, convection and diffusion, as wellas any source terms associated to give a convergent, consistent and stable solution.The underlying physical phenomena is usually complex and non-linear, hence thesolver should have a stable time marching or iterative scheme for computing the so-lution. The CFD solver code ANSLib [33] developed by the Advanced NumericalSimulations Laboratory (ANSLab) has all the characteristics of being a good solverand is chosen for the numerical experiments conducted in this research. It uses anedge based formulation to handle mixed grids with a robust numerical scheme.ANSLib allows the users to specify physics of their problem by providing codesnippets to compute interior and boundary fluxes, source terms, initial conditionsand by specifying constraints on the solution at the boundaries.The final stage of CFD is post-processing. The visualization of vector plots,domain geometry, 2D and 3D surface plots, solution and error contours come underpost-processing. There are many tools available in the market for post-processing.We have chiefly employed Paraview and MATLAB for post-processing along withcommon plotting tools.1.1 MotivationThe numerical simulation of fluid flows is affected by three forms of errors, namely:the geometry modeling, physical modeling and numerical errors. Geometry model-ing errors are related to the transformation from the physical domain to a computa-tional domain by either omitting or simplification of complex parts of the geometry.They are usually unavoidable and can help allocate resources to capture importantfeatures of the flow. Omitting rivets in the structures or modification of a pointedtrailing edge are examples of geometry modeling.Physical modeling errors pertain to the effective modeling of physical phenom-ena, for instance turbulence or real gas effects, by mathematical models. Oftenmodeling is required for turbulence, boundary conditions, multi-phase models, etc.The inadequacy in predicting the physical phenomena accurately by such models71.1. Motivationleads to physical modeling errors. Both geometry and physical modeling errors areoutside the scope of this thesis, which focuses on numerical errors.Current experiments have shown that the numerical errors are as large as thephysical modeling errors and it is significantly important to reduce them. The firstthree AIAA Drag Prediction Workshops (DPW) [23, 21, 26] used similar winggeometries to compare the computational data with the experimental data. The firsttwo workshops used a DLR-F4 and DLR-F6 wing geometry respectively and theresults revealed significant scatter in terms of drag prediction and agreement withthe experimental data. The scatter was attributed to regions of turbulent separatedflow in the test cases and a wing body fairing was designed and added to the DLR-F6 wing body configuration in the third DPW workshop, so as to suppress flowseparation in the wing-body juncture. The results for the third workshop and theoverview of the three workshops [28] revealed that the scatter was not just a resultof physical modeling errors but also numerical errors. The first three workshopsalso established that the use of different solvers, turbulence models, grid types orgrid constructions all contribute to the variation in the results.The fifth AIAA Drag Prediction Workshop [24] focused on reducing grid-related errors by performing a grid refinement study on a common grid sequencederived from a multiblock structured grid. The study had six different levels of gridrefinement ranging from 0.64×106 cells to 136×106 cells, with structured oversetand hexahedral, prismatic, tetrahedral and hybrid unstructured grid formats testedwith different turbulence models. The study concluded that there was no clearbreak-outs with regards to grid type or turbulence model with discretization errorsand turbulence modeling still being the significant contributors to the scatter in re-sults. Overall, different grid types, grid constructions, and discretization schemeshad different results even when the solution was converged.Discretization errors, which are a major source of numerical errors, are oftenmost difficult to calculate for general flows. Discretization error arises out of theinterplay between the chosen discretization scheme, the mesh quality (cell size,shape, anisotropy, etc.), the mesh resolution and the solution and its derivatives.Roy [36] has discussed the issues related to discretization error estimators andhas emphasized that for any method to give a reliable estimate, it is necessary toachieve the asymptotic range of solution convergence. As difficult as it is to esti-mate discretization error, we utilize the understanding that it is transported throughthe domain in much the same fashion as the solution in the underlying mathemat-ical model [14, 35] and is locally generated according to the truncation error. Roy[36] showed that for a linear differential operator L , the truncation error τ can beused to estimate the discretization error ε , asL (ε) =−τ (1.1)81.1. MotivationFigure 1.3: One dimensional grid for convection-diffusion equationknown as the error transport equation. Banks et al. [1] extended the error transportmethod to estimate discretization error for finite volume discretizations of nonlin-ear hyperbolic equations and demonstrated the approach to provide useful errorestimates. The scope of this research is limited to linear problems, however. Equa-tion 1.1 shows that the discretization error evolves the same way as the solutiondoes, with the truncation error being the driving source. For some cases, the care-ful refinement of the mesh in regions of high truncation error can lead to a reductionin discretization error as shown in [46]. Mesh adaptation based on truncation er-ror is observed to give the lowest discretization errors [46, 35] and motivates us tounderstand the relation between the two types of errors.In this thesis, we concern ourselves with a general approach to reduce trunca-tion error irrespective of the solution type and do not attempt at either estimationof truncation error or mesh adaptation. Consider a one dimensional convection-diffusion equation discretized using the finite volume approach on a one dimen-sional grid [17]. Most flows fall into the category of convection and diffusion andhence it serves as an ideal model to explain the idea behind our approach.L (φ) = ∂φ∂x −ν∂ 2φ∂x2 = 0 (1.2)A finite volume discretization of the equations requires the averaging over thecontrol volume i. The advection term when averaged over the control volume takesthe form,∂φ∂x =1∆xˆ xi+ 12xi− 12(∂φ∂x)dx =φi+ 12 −φi− 12∆x(1.3)The solution values are stored at the nodes for the uniform mesh shown in theFigure 1.3. To calculate the averaged gradient in the Equation 1.3, one would needthe values of the solution variable on the control volume boundaries. For a general91.1. Motivationunstructured mesh, we would need solution reconstruction to better approximatethe required data. However, for purposes of illustration here, we use a first or-der upwind approximation to estimate the solution variables at the control volumeboundaries. Thus we have,φi+ 12 ≈ φ i, φi− 12 ≈ φ i−1 (1.4)and the averaged gradient would be∂φ∂x =φ i−φ i−1∆x(1.5)To determine the accuracy of the scheme, we use Taylor series expansion tocalculate the control volume averages in terms of the solution and derivatives at thereference point of control volume i.φ i =1∆xˆ xi+ 12xi− 12(φi +∂φ∂x∣∣∣∣i(x− xi)+12∂ 2φ∂x2∣∣∣∣i(x− xi)2 +16∂ 3φ∂x3∣∣∣∣i(x− xi)3 . . .)= φi +124∂ 2φ∂x2∣∣∣∣i(∆x)2 +O(∆x)4 (1.6)φ i−1 =1∆xˆ xi− 12xi− 32(φi +∂φ∂x∣∣∣∣i(x− xi)+12∂ 2φ∂x2∣∣∣∣i(x− xi)2 +16∂ 3φ∂x3∣∣∣∣i(x− xi)3 . . .)= φi−∂φ∂x∣∣∣∣i(∆x)+1324∂ 2φ∂x2∣∣∣∣i(∆x)2−524∂ 3φ∂x3∣∣∣∣i(∆x)3 +O(∆x)4 (1.7)Substituting the above expansions into the Equation 1.5, the resultant gradientturns out to be first order accurate with a leading error term of−12∂ 2φ∂x2∣∣∣i(∆x) whichis the truncation error for the gradient approximation.∂φ∂x =∂φ∂x∣∣∣∣i−12∂ 2φ∂x2∣∣∣∣i(∆x)+524∂ 3φ∂x3∣∣∣∣i(∆x)2 +O(∆x)3 (1.8)Applying a similar procedure for the diffusion term in Equation 1.2, the cen-tered difference results in a second order accurate approximation:∂ 2φ∂x2 =∂ 2φ∂x2∣∣∣∣i+18∂ 4φ∂x4∣∣∣∣i(∆x)2 +O(∆x)4 (1.9)101.1. MotivationThe discrete approximation to the convection-diffusion equation can finally bewritten as:Lh(φ) =∂φ∂x −ν∂ 2φ∂x2 =(∂φ∂x∣∣∣∣i−ν ∂2φ∂x2∣∣∣∣i)−12∂ 2φ∂x2∣∣∣∣i(∆x)+ (1.10)524∂ 3φ∂x3∣∣∣∣i(∆x)2−ν8∂ 4φ∂x4∣∣∣∣i(∆x)2 + . . .The operatorLh is the discrete operator whereas on the right hand side, we canidentify the continuous PDE for convection-diffusion. The entire equation can bere-written in the form derived in the works of Roy [35] as shown in the Equation1.11. A careful substitution of the discrete and exact solutions to the PDE will leadus to the derivation shown in the Equation 1.1τ (φ) =Lh(φ)−L (φ) (1.11)The truncation error in Equation 1.10 highlights that the discretization is for-mally first order and that the leading error term adds extra diffusion into the system.The even order derivatives have a tendency to smear the sharp gradients of the so-lution and are said to be diffusive in nature whereas odd order derivatives inducean oscillatory behavior in the smooth regions of the solution and are said to bedispersive in nature.Truncation error is dependent on the derivatives of the solution and their asso-ciated coefficients. Schemes which are independent of the solution derivatives buteffective in reducing the truncation error coefficients will in general reduce trunca-tion error and in turn discretization error. Previous work has shown that the trun-cation error coefficients are influenced by the selection of reconstruction methods,face gradient averaging schemes, jump terms, etc. [17, 30, 45] in the numericalschemes. This observation lays the foundation for this thesis where the goal isto analyze the effects of different discretization schemes and boundary conditionimplementations in bringing out the desired results in solution accuracy.Jalali and Ollivier-Gooch [18, 19] have performed analytical and numericalexperiments on cell-centered finite volume discretizations of diffusion and con-vection equations. Their analysis involved the computation of the coefficients ofTaylor series expansion for the truncation error and have demonstrated the methodsto quantify both the truncation and discretization errors. Yan et al. [45] performedtruncation and discretization error analysis for diffusion problems based on vertex-centered discretization. In solving Poisson’s equation using the finite volume ap-proach, Diskin et al. [10] compared cell-centered and vertex-centered schemes and111.2. Objectivesconcluded that there is little difference in the accuracy of the two approaches at anequivalent number of degrees of freedom.The previous works on diffusion and convection were restricted to the interiordomain of the mesh and the influence of discretization schemes on the boundarycontrol volumes is largely unclear. Moreover, it is not known which implementa-tions of the boundary conditions are more accurate. The effect of different solutionreconstruction methods, face gradient calculations, boundary condition implemen-tation methods and modifications to meshes on the boundary truncation error andoverall discretization error needs to be studied. Given that the discretization errorand truncation error are related to each other, the understanding of the relationshipbetween the two errors is not very clear, either.Thus, our primary objective for this research is to explore the ways in whichtruncation error in the boundary regions of a mesh can be reduced while the overallgoal is to reduce both the truncation and discretization errors globally. A detailedobjective of this thesis is presented in the next section.1.2 ObjectivesIn this thesis, we define our problem statement as the truncation error analysis fordiffusion schemes in the boundary regions of unstructured meshes and attempt toreduce interior and boundary truncation error as well as the solution error. Pois-son’s equation is used as the model diffusion equation and vertex-centered controlvolumes are defined for the finite volume discretization. We study the influence ofdifferent discretization schemes and boundary condition implementation methodson truncation and discretization error and also try to decipher the relationship be-tween the truncation error and the solution accuracy. The scope of this research islimited to isotropic and stretched 2-D meshes, both triangular and mixed.Among the type of discretization schemes, we chiefly inspect the differencesin the use of popular solution reconstruction methods which are the least-squaresand Green-Gauss methods. Green-Gauss method has been traditionally popularfor viscous flows [9, 13, 27]. It has also been noted to suffer from decouplingissues when quadrilateral meshes are used. Thomas et al. [40] and Diskin andThomas [9] have shown remedies to counter this issue associated with Green-Gaussmethod. The least-squares method is one of the preferred methods for solutionreconstruction in the literature [8, 31, 20, 25] and was found to be comparableor even better than the Green-Gauss method for accurate gradient calculation. Itpaves the way for higher order reconstruction whereas the Green-Gauss method isrestricted to second order. Apart from gradient calculation at the center of controlvolumes, the calculation of face gradients on control volume boundaries influence121.2. Objectivesthe accuracy significantly [6, 17] and is a concern addressed in this thesis.The use of ghost cells as a method of boundary condition implementation iswell known and a brief description of the method can be found in Blazek [4] andLeveque [22]. Cary et al. [6] used ghost cells for boundary condition implementa-tion in Boeing’s BCFD solver in a cell-centered approach and found an improveddrag prediction for a NACA 0012 airfoil section using Euler’s equations. How-ever, the impact of ghost cells in reducing the truncation error is unknown and westudy this boundary condition implementation on vertex-centered approach with asecond order reconstruction in the ghost cells.An alternate method of applying the boundary condition is to use the boundarycondition information for solution reconstruction in the boundary control volumesand thereby calculating the flux at the boundary edges. This approach is presentedin the works of Ollivier-Gooch and Van Altena [34], where the boundary conditionsat the Gauss points for flux integration are used as hard constraints, which alongwith information from neighboring control volumes, is used to reconstruct gradi-ents in a least-squares sense. This approach and its variants are explored to deviseschemes efficient in reducing truncation error and improve solution accuracy.Finally, the use of jump terms [30] to introduce a damping term based on thejump in solution values at control volume interfaces has been proven to give accu-rate solution for diffusion problems. Jalali [17] and Yan et al. [45] have found jumpterms to be effective in modifying truncation error coefficients for interior controlvolumes. The current research explores the potential of this method to devise newboundary treatment methods and its impact on boundary truncation error as well.The evaluation of different discretization schemes and boundary treatment meth-ods is conducted in two stages. In the first stage, a qualitative and quantitativecomparison is performed analytically for different test cases for meshes which aregeneral in nature. These tests are performed to ascertain the asymptotic behav-ior of the schemes developed. They also present a platform to verify the schemeson a variety of systematically perturbed meshes so that the implementation can becorrected, compared and schemes performing below expectations be eliminated.In the second stage, the numerical experiments are performed with schemes fil-tered from analytical tests to re-affirm the observations found in analytical results.The numerical tests provide a robust testing mechanism where the schemes aretested with different boundary conditions, solution types and meshes to assess theapplications and limitations of schemes developed. Isotropic and stretched meshsequences are generated for pure triangular and mixed meshes along with isotropicpure quadrilateral meshes. Each sequence has four systematically refined meshesto study the asymptotic convergence of schemes as well.Thus, results for a variety of different treatments to the boundary cells, solutionreconstruction and boundary stencil are presented in this work with emphasis laid131.3. Outlineon reducing the truncation error in the boundary regions. At the same time, effortsare made to realize the best scenarios where the discretization error and overalltruncation error, in interior and boundary regions are minimized.1.3 OutlineChapter 2 lays the theoretical foundation for the research presented in this thesis.An overview of the finite volume method used for solving the diffusion equationby constructing vertex-centered control volumes is given in this chapter. The keyaspects of the solver which include solution reconstruction, face gradient calcu-lation and flux integration, are discussed in detail. The differences in the interiorand boundary control volumes are highlighted with respect to the reconstructionstencils and other relevant factors. The chapter discusses in detail the schemesdeveloped or considered to test for reducing the truncation errors which includemodification to the stencil, use of boundary constraints and the use of jump terms.The methodologies for the calculation of truncation error, both analytically andnumerically, are presented in Chapter 3. The Method of Manufactured Solutions isalso discussed for the calculation of discretization error.Chapter 4 presents the analytical and numerical results for pure triangularmeshes. The analytic section presents the challenges faced in the implementationof the methods and the remedies applied to standard methods for obtaining betterresults. The schemes performing better in the analytic section are taken forward forthe implementation in the numerical solver. Chapter 5 presents the analytical andnumerical results for the mixed and pure quadrilateral meshes and the challengeswith the implementation for mixed meshes are highlighted in the analytic and nu-merical sections. The results for analytic and numerical sections are compared andcontrasted with the results for triangular meshes. The differences and similaritiesin isotropic and stretched boundary layered meshes are also presented in both thechapters. Overall, the schemes performing best for reducing the truncation anddiscretization errors are explored to seek guidelines for achieving simultaneousreduction in different forms of errors for all the types of mesh tested in this thesis.The summary of approaches taken in this thesis are presented in the Chapter 6along with concluding remarks drawn from the results. The scope for future workis also laid out in this chapter.14Chapter 2BackgroundA background of the CFD practices used in the current research is presented in de-tail in this chapter. First a brief overview of the relevant aspects of CFD for solvinga diffusion problem using the finite volume discretization is given. The formulationfor solution reconstruction using methods like least-squares and Green-Gauss areexplained. Next, numerical flux calculation and rules for flux integration are given.Finally, the rationale and methodology for different boundary treatment methodsexperimented with in this thesis are presented.2.1 Finite Volume ApproachThe finite volume method is a robust method for solving fluid flow problems as itdirectly utilizes the conservation laws of fluids. The first step of the finite volumeapproach is to divide the physical domain into arbitrary polyhedral (3D) or polyg-onal (2D) grid elements and then to discretize the integral form of the governingequations over each control volume. There are two basic approaches to definingthe control volumes with respect to the grid. They are the (i) cell-centered ap-proach and (ii) vertex-centered approach. In cell-centered approach, the controlvolumes are identical to the primal grid cells and the solution values are stored atthe centroids of the primal cells. On the other hand in vertex-centered approach,the solution values are stored at the vertices and the control volumes are the meshdual. In our current research we limit ourselves to two dimensions although theextension to three dimensions for the schemes presented in this thesis is prettystraightforward. We use the vertex-centered approach with vertices as the refer-ence points and the control volumes defined by the median dual. The median dualis constructed by connecting the mid points of the faces and the centroids of theprimal cells sharing the reference vertex. Figure 2.1 gives an illustration of thedifferent approaches discussed here.In current research, the Poisson equation,∇2φ = S, is used as a model equationfor diffusion. The differential form of Poisson’s equation is integrated over thecontrol volume and the application of the divergence theorem gives us the finitevolume formulation for Poisson’s equation, taking the form as shown:152.2. Solution Reconstruction(a) Cell-Centered control volume (b) Median Dual control volumeFigure 2.1: Control volume illustration˛∂V∇φ .nˆds =˛VSdA (2.1)where nˆ is the outward normal vector and ds is the infinitesimal for the integrationalong the edges of the control volume. The solution for the Poisson’s equation isa steady state solution where the left hand side of Equation 2.1 is the flux integralalong the perimeter of the control volume. The elliptic nature of Poisson’s equa-tion implies that the solution in the interior is strongly affected by the boundaryconditions.2.2 Solution ReconstructionThe simplest method of reconstruction is to assume a constant solution value inthe control volumes, which is basically a first order reconstruction. However, thescheme is too diffusive and leads to excessive growth of shear layers [4]. Thus,a second or higher order reconstruction is needed. The Green-Gauss method ispopularly employed for second order accurate flux calculations in the literature[9, 13, 27, 15]. Higher order accurate methods first appeared in literature as earlyas the 1990s, when Barth and Fredrickson [3] presented a method for k-exact re-construction and use of Gaussian quadrature rules for high order flux evaluations.The works of Ollivier-Gooch and Van Altena [34] detail the inclusion of mean andboundary constraints in the framework of developing higher order schemes for theadvection-diffusion equation.162.2. Solution ReconstructionThe evaluation of solution and gradients to the desired level of accuracy is nec-essary for accurate flux integration and ultimately a steady solution to the Poisson’sequation. The following sections discuss the traditionally popular Green-Gaussand least-squares methods to obtain a second order accurate solution reconstruc-tion. The Green-Gauss method presented in this work assumes a constant gradientin the primal grid cells and the gradient calculation takes an algebraically simplerform for triangular mesh.Next, the framework for the least-squares method is presented where the con-ditions to be satisfied for a consistent solution reconstruction of desired accuracyare laid out. This is followed by the weighted second order least-squares formu-lation where the unknown derivatives of solution variables are evaluated followedby the advantages of the least-squares method over the Green-Gauss method. Atthe end of this section we also highlight the subtle differences in the reconstructionstencils for interior and boundary control volumes which we anticipate to affect theaccuracy of gradient calculation and in turn truncation and discretization errors.2.2.1 Green-Gauss MethodThe gradients required for numerical flux evaluation are directly calculated usingthe Green-Gauss method. The gradients are calculated using the Green-Gauss for-mula [11, 43, 15] over the primal cells even though the FVM formulation is vertexbased. Let us consider a mixed element grid as shown in the Figure 2.2. For fullytriangular grids, the formulation is equivalent to a Galerkin finite element schemewith a linear basis function [2]. The Green-Gauss formula isˆV∇φdA =˛∂Vφ nˆds (2.2)The gradients at the dual face BC contained in the triangular cell is calculated as∇φBC = ∇φ021 (2.3)=12A021[(φ 0 +φ 2)(y0− y2)+(φ 2 +φ 1)(y2− y1)+(φ 1 +φ 0)(y1− y0)−(φ 0 +φ 2)(x0− x2)− (φ 2 +φ 1)(x2− x1)− (φ 1 +φ 0)(x1− x0)]where ∇φ021 is the gradient evaluated on the primal cell identified by the subscripts.For the dual face AB contained in the quad cell, the gradient is evaluated using acombination of the Green-Gauss formula and an edge derivative correction to avoiddecoupling and give improved gradients [4, 11].∇φAB = ∇φ0432 +[φ 2−φ 0|r2− r0|−∇φ0432 · e02]e02 (2.4)172.2. Solution ReconstructionFigure 2.2: Gradient reconstruction for mixed cell interface using Green-Gaussmethodri is the position vector for the node i whereas e02 is the unit vector in the directionof the edge [0,2] given bye02 =r2− r0|r2− r0|(2.5)Note that for both the dual faces AB and BC, the edge derivative is recovered.Mathematically, this means∇φBC · e02 = ∇φAB · e02 =φ 2−φ 0|r2− r0|(2.6)The above gradient calculation is found to be make the Green-Gauss gradientsless susceptible to the odd-even decoupling [16]. The scheme presented in thissection has been shown by Thomas et al. [40] and Diskin and Thomas [9] topossess second order accuracy for general isotropic mixed meshes with regards todiscretization error.Our experiments revealed that the truncation error was more sensitive to thefact that the vertex doesn’t coincide with centroid of the control volume. Thisdiscrepancy is observed to give O(h−1) order for truncation error where h is thecharacteristic mesh size and was remedied by using centroid locations instead ofnode locations for the gradient calculation.2.2.2 Least-Squares MethodThe first step in using the least-squares method is to formulate the conditions nec-essary for a consistent reconstruction scheme. The reconstruction design crite-182.2. Solution Reconstructionria [3, 32] requires the method to conserve the mean of control volumes, havea compact support and exactly reconstruct polynomials of required degree. Theframework for least-squares is built in the following sub-section by satisfying thereconstruction design criteria followed by the second order weighted formulationto obtain the unknown derivatives.2.2.2.1 Framework for Least-Squares MethodHigher order accurate reconstruction and an equally accurate flux integration methodare essential for an accurate finite volume scheme. In this sub-section, we limit thediscussion to the development of a general reconstruction framework to obtain thedesired level of accuracy. Our goal is to first express the solution as a continuousTaylor series polynomial about the reference point (xi,yi)φRi (x,y) = φ |i +∂φ∂x∣∣∣∣i(x− xi)+∂φ∂y∣∣∣∣i(y− yi)+ (2.7)12∂ 2φ∂x2∣∣∣∣i(x− xi)2 +∂ 2φ∂x∂y∣∣∣∣i(x− xi)(y− yi)+12∂ 2φ∂y2∣∣∣∣i(y− yi)2 + ...In Equation 2.7, the solution φRi is the reconstructed solution at any point (x,y)about the reference control volume i provided we know the solution values and thederivatives ∂k+lφi∂ kx∂ ly at that point. The derivatives and the solution values are the coef-ficients of the Taylor polynomial and the degree of this polynomial determines theorder of reconstruction. Second order accuracy for instance can be obtained if weknow the gradients in x and y, as can be seen in Equation 2.8. Such a reconstructionis linearly varying in the control volume. The implications can be seen in a cell-centered scheme where the pointwise solution value at the centroid naturally rep-resents the control volume average at second order accuracy. In a vertex-centeredapproach however, the control volume average differs from the pointwise solutionvalue at the vertex with a first order error.φRi (x,y) = φ |i +∂φ∂x∣∣∣∣i(x− xi)+∂φ∂y∣∣∣∣i(y− yi)+O(δx2,δy2) (2.8)The coefficients of the Taylor polynomial are chosen such that the error insatisfying the mean solution value in the stencil control volumes is minimized.Enforcement of the mean in the reference control volume requires that192.2. Solution Reconstructionφ i =1AiˆViφRi (x,y)dA (2.9)Substituting Equation 2.7 in the mean constraint requirement, one can see that,φ i =1AiˆViφRi (x,y)dA = φ |i +∂φ∂x∣∣∣∣ixi +∂φ∂y∣∣∣∣iyi + (2.10)12∂ 2φ∂x2∣∣∣∣ix2i +∂ 2φ∂x∂y∣∣∣∣ixyi +12∂ 2φ∂y2∣∣∣∣iy2i + ...where the geometric moments for the reference control volumes are defined asxpyqi =1AiˆVi(x− xi)p(y− yi)qdA (2.11)By observation, the above integral can be easily calculated by using Greenstheorem to convert into a line integral along the boundary of control volume. Theintegration can now be evaluated using Gaussian quadrature rules for desired accu-racyxpyqi =1Ai(p+1)ˆ∂Vi(x− xi)p+1(y− yi)qdy (2.12)For k exact reconstruction, the Taylor series expansion of φRi must be carriedout through until the kth derivatives. The derivatives are computed by seeking tominimize the error in predicting the mean value of the solution for control volumesin the stencil {Vj}i. The mean value for neighboring control volumes is calculatedas1A jˆVjφRi (x,y)dA = φ |i +∂φ∂x∣∣∣∣i1A jˆVj(x− xi)dA+ (2.13)∂φ∂y∣∣∣∣i1A jˆVj(y− yi)dA+12∂ 2φ∂x2∣∣∣∣i1A jˆVj(x− xi)2dA+∂ 2φ∂x∂y∣∣∣∣i1A jˆVj(x− xi)(y− yi)dA+12∂ 2φ∂y2∣∣∣∣i1A jˆVj(y− yi)2dA+ ...The calculation of moments for every neighboring control volume in the stencil{Vj}i about reference control volume is simplified by replacing (x− xi) and (y−202.2. Solution Reconstructionyi) with (x− x j)+ (x j− xi) and (y− y j)+ (y j− yi) respectively. Expanding andintegrating Equation 2.131A jˆVjφRi (x,y)dA = φ |i +∂φ∂x∣∣∣∣i(x j +(x j− xi))+∂φ∂y∣∣∣∣i(y j +(y j− yi))+ (2.14)12∂ 2φ∂x2∣∣∣∣i(x2 j +2x j(x j− xi)+(x j− xi)2)+∂ 2φ∂x∂y∣∣∣∣i(xy j + x j(y j− yi)+ y j(x j− xi)+(x j− xi)(y j− yi))+12∂ 2φ∂y2∣∣∣∣i(y2 j +2y j(y j− yi)+(y j− yi)2)+ ...The geometric terms from the Equation 2.14 can be written in a general formasx̂pyqi j =1A jˆVj((x− x j)+(x j− xi))p((y− y j)+(y j− yi))qdA (2.15)=q∑l=0p∑k=0(ql)(pk)(x j− xi)k(y j− yi)l xp−kyq−l jUsing the hat notation for the geometric terms, Equation 2.14 is re-written asφ j = φ |i +∂φ∂x∣∣∣∣ix̂i j +∂φ∂y∣∣∣∣iŷi j + (2.16)12∂ 2φ∂x2∣∣∣∣ix̂2i j +∂ 2φ∂x∂y∣∣∣∣ix̂yi j +12∂ 2φ∂y2∣∣∣∣iŷ2i j + ...The mean constraint for every control volume in the the stencil {Vj}i is ex-pressed in the form given by Equation 2.16. The choice of the stencil should becompact and limited to local neighbors so that only valid data is used for recon-struction. Neighbors which are physically far from the vertex are not consideredin the stencil. Using only a small topological neighborhood of a vertex ensuresa compact physical support. The stencil for second order reconstruction includesthe control volumes adjoining the reference control volume whereas for third andfourth order accuracy, second and third neighbors may also be included. Stencildetermination is done as a pre-processing step and for control volumes which aretoo close to the boundary, more neighbors can be added to the stencil to have suf-ficient support for achieving desired accuracy of reconstruction. Once the stencil212.2. Solution Reconstructionis determined, the information from each neighbor is emphasized on the basis ofproximity by using inverse distance weighting wi j. The inverse distance weighingis used as it is observed to improve conditioning and recover accuracy in gradientcalculations for isotropic meshes [25]. The resulting system of equations to besolved for the Taylor coefficients takes the form1 xi yi 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 · · ·φ∂φ∂x∂φ∂y∂ 2φ∂x2∂ 2φ∂x∂y∂ 2φ∂y2...i=φ iwi1φ 1wi2φ 2wi3φ 3...wiNφN(2.17)where the line separator is used to emphasize the mean constraint applied to thesystem of equations to be solved using least-squares andwi j =1|r j− ri|(2.18)is the geometric weight used. In the solver, the mean constraint is eliminated ana-lytically and the unexhausted equations are solved in a least-squares sense for everycontrol volume by the singular value decomposition (SVD) method.2.2.2.2 Weighted Least-Squares MethodThe second order weighted least-squares system is derived from the equation (2.17)which gives a comprehensive formulation to obtain higher order reconstruction forany control volume. The resulting system of equations with a mean constraint tobe solved is simplified as below:1 xi yiwi1 wi1x̂i1 wi1ŷi1wi2 wi2x̂i2 wi2ŷi2wi3 wi3x̂i3 wi3ŷi3.........wiN wiN x̂iN wiN ŷiNφ∂φ∂x∂φ∂yi=φ iwi1φ 1wi2φ 2wi3φ 3...wiNφN(2.19)The mean constraint is applied using Gauss elimination and the set of equationstakes a compact form as shown222.2. Solution Reconstructionwi1(x̂i1− xi) wi1(ŷi1− yi)wi2(x̂i2− xi) wi2(ŷi2− yi)wi3(x̂i3− xi) wi3(ŷi3− yi)......wiN(x̂iN− xi) wiN(ŷiN− yi)(∂φ∂x∂φ∂y)i=wi1(φ 1−φ i)wi2(φ 2−φ i)wi3(φ 3−φ i)...wiN(φN−φ i)(2.20)For second order reconstruction, the natural choice for choosing the stencilwould be all the vertex neighbors connected by edges to the reference vertex, whichin general vary from 5-7 for a triangular mesh and 4 for a quadrilateral mesh. Thedetails of stencil selection, differences in stencil choice for triangular and mixedmeshes, and interior and boundary control volumes are presented in the Section2.2.3.2.2.2.3 Advantages of Least-Squares MethodUsing the least-squares method of reconstruction offers advantages and simplifi-cations when compared to the Green-Gauss method. The reconstruction is gridtransparent in the sense that it does not rely on the mesh type; whether it is aquadrilateral or a triangular mesh. The set of equations solved are related to a sten-cil which identifies relevant neighboring points for the reconstruction. Althoughthe choice of stencil can be arbitrary, the most obvious choice is the set of con-trol volumes adjacent or physically closer to the reference point. The least-squaresmethod can be extended to calculate higher order derivatives and successively ahigher order solution reconstruction which can be used to calculate both the ad-vective and diffusive fluxes. On the other hand, the Green-Gauss method has beenused predominantly for viscous fluxes along with additional schemes for calcu-lating inviscid fluxes. The extension of the Green-Gauss method beyond secondorder accuracy is rarely explored. Thus the least-squares method has the benefitof providing with a single reconstruction that can be used for calculating both theinviscid and diffusive fluxes.Further advantages are discussed in later sections where we show methods toapply boundary constraint and the use of jump terms (discussed in section 2.3.1)which are proven to be crucial in our objective of reducing truncation and dis-cretization errors. The boundary constraints, applied at the Gauss quadrature pointsare such that the reconstructed solution satisfies the boundary condition exactly atthe quadrature points. The Green-Gauss and least-squares methods were found tohave comparable accuracies in reconstruction for isotropic meshes [25], [17] butfor anisotropic meshes, the gradients were found to be inaccurate. The reader is232.2. Solution Reconstructiondirected to the works of Jalali and Ollivier-Gooch [20] for details about accurate re-construction on unstructured anisotropic meshes for second, third and fourth orderaccurate schemes using weighted and unweighted least-squares method.2.2.3 Reconstruction StencilsBarth and Frederickson [3] and Ollivier-Gooch [31] highlight the importance ofchoosing a compact support implying that the reconstruction should use data froma stencil whose members are physically and topologically closer. The size of sucha stencil is defined by the number of derivatives to be evaluated, usually larger bya certain factor for least squares problem. Including additional neighbors allowsleeway for ignoring non-smooth data as well. For second order accuracy, the firstvertex neighbors are usually sufficient, whereas for third and fourth order accuracy,the second or third neighbors may be included. Using second neighbors is suffi-cient for reconstruction up to fourth order for interior control volumes but not forboundary control volumes. Quite often, an additional layer of neighbors are addedfor boundary control volumes to provide the necessary support.The nearly symmetric position of the stencil neighbors pave the way for errorcalculation on all the sides of a control volume in the interior. In an equilateraltriangular mesh or a quad mesh with rectangular elements, the cancellation is exactwhich yields higher orders of truncation error. Perturbed meshes exhibit lower or-ders of truncation error while the boundary control volumes with stencil membersdominant on one side have higher error magnitudes. Lack of symmetry can also beobserved in the interface layers of mixed meshes where the truncation error is alsoof lower order.The choice of stencil shown in the Figure 2.3 for triangular and mixed mesheshighlight the members needed for different levels of accuracy. The red markedneighbors for second order reconstruction are the first neighbors which share theboundary edges of the control volume. Hence we have four neighbors for a trian-gular boundary element and three for a quad boundary element. Also notice thatthe interface layer control volumes where the triangular and quad layer meet haveabout five neighbors. The stencil choice for third and fourth order reconstructionfor interior control volumes coincide whereas for the boundary control volumes,the third neighbors are required for fourth order accuracy. The use of 3, 8 and14 neighbors for second, third and fourth order accuracy, (to calculate 2, 5 and 9derivatives) seem to be sufficient [31].242.2.SolutionReconstruction(a) Stencil Selection for pure triangular mesh (b) Stencil selection for mixed meshFigure 2.3: Stencil choices for reconstruction252.3. Flux Evaluation2.3 Flux Evaluation2.3.1 Numerical Flux CalculationEvaluation of numerical flux is a two part problem: firstly the solution reconstruc-tion at the edge Gauss points must be performed and secondly ensuring that theflux across the edge is same for the control volumes sharing the edge. Previoussections have shown methods of solution reconstruction using least-squares andGreen-Gauss methods. This section explains the methods for calculating numeri-cal flux for the two reconstruction methods such that flux across the edge is unique.In the least-squares reconstruction method, the gradients at the edge are not thesame. Hence to calculate the flux, the gradients need to be averaged. Jalali [17]has performed multiple tests for assessing different averaging methods, additionof finite difference terms and jump terms for the calculation of face gradients toevaluate diffusive fluxes. The analysis has favored the use of linear and arithmeticaveraging for gradient calculation and found them to be comparably less erroneousthan volume weighted averaging. The results also found linear interpolation to bemarginally better than arithmetic averaging for cell-centered schemes. The inclu-sion of jump terms also reduces the errors.To evaluate the numerical flux, we perform the solution reconstruction. Oncethe Taylor polynomial for reconstruction is evaluated, as elaborated in Section 2.2,the gradient for the control volume can be approximated as∇φRi (x,y) =∂φ∂x∣∣∣i+ ∂2φ∂x2∣∣∣i(x− xi)+∂ 2φ∂x∂y∣∣∣i(y− yi)+ · · ·∂φ∂y∣∣∣i+ ∂2φ∂y2∣∣∣i(y− yi)+∂ 2φ∂x∂y∣∣∣i(x− xi)+ · · · (2.21)Equation 2.21 is a general form for evaluating the gradients. For second orderaccurate reconstruction where the solution varies linearly in the control volume, thegradients are constant over the control volume. The need for averaging of gradientsarises only when there is a discontinuity in gradients across the control volumeboundaries, as is the case for least-squares reconstruction. In our research, we useboth linear interpolation and arithmetic averaging to calculate the face gradients atthe control volume boundaries. The averaged gradient at the flux integration Gausspoint g (see Figure 2.4) is shown in Table 2.1 for arithmetic and linear interpolationmethods.A damping term based on the jump in solutions across the edge was proposedby Nishikawa [30] and was found to be effective in high frequency error damping.This has been demonstrated to give more accurate solutions for diffusion problems262.3. Flux EvaluationAveraging Method FormulaArithmetic Averaging (AA) 12∇φi +12∇φ jLinear Interpolation (LI) |r jF ||riF |+|r jF |∇φi + |riF ||riF |+|r jF |∇φ jTable 2.1: Face gradient calculationFigure 2.4: Linear interpolation method for numerical fluxon isotropic and anisotropic meshes for cell-centered [17] and vertex-centered [45]schemes. In some cases, the lack of a damping term leads to unstable or incon-sistent results. This term can be efficiently and systematically integrated with anadvection scheme to solve advection-diffusion problems. The damping term, alsoknown as jump term, is added to the averaged flux in the face normal direction as∇φJump =α|ri j · nˆF |(φR−φL)nˆF (2.22)where ri j is the vector joining the reference centers of the two control volumes andnˆF is the face normal at the Gauss point g. The left and right control volume solu-tions at the Gauss point are φL and φR. The jump coefficient α plays an importantrole in controlling the damping effect of the jump term. Certain values of α arefound to increase the order of accuracy for truncation errors on uniform meshes,and in general reduce the errors in truncation and discretization significantly.The calculation of numerical flux has been focused on using the least-squaresgradient reconstruction so far. However, the Green-Gauss method presents a rathersimple approach to calculate the numerical flux. Since the gradient in the primalcells is constant, the gradient across the dual edges of control volumes is the sameon both sides and hence the numerical flux is the dot product of the gradient andunit normal to the edge. For instance, the numerical flux for each of the edges272.3. Flux EvaluationFigure 2.5: Numerical flux calculation for Green-Gauss methodof control volumes 0,1,2 in the primal triangular cell shown in the Figure 2.5 issimply the dot product of gradient of the cell calculated by the formula given inEquation 2.3 with the unit normal for corresponding dual edges in the primal cell.2.3.2 Flux IntegrationDiscretized equations of the form given in Equation 2.1 are to be evaluated for ev-ery control volume in the mesh. The accuracy of solution hinges on the accuracyof the flux integral which in turn depends on the numerical flux accuracy and in-tegration methods. Gauss quadrature rules of integration are used to calculate fluxintegral of desired accuracy. The quadrature rules calculate high order accurateflux integrals by evaluating the integrand at the Gauss points and by multiplyingthem with appropriate weights.To obtain a second order accurate flux integral over a single edge, the numericalflux is evaluated at the mid point of the edge and multiplied by the edge length. Ahigher order flux integral needs more than one Gauss point. Figure 2.6 gives anillustration of the Gauss quadrature points. The flux integral for a control volumeis a direct summation of the integration on each edge. Detailed information aboutthe Gauss points and the appropriate weights for integration can be found in theworks of Ollivier-Gooch and Van Altena [34] and Van Altena [43]Once the flux integral is calculated for every control volume, we obtain a sys-tem of linear equations which can be solved by any matrix solution technique toobtain the steady state Poisson solution.282.4. Boundary Treatment MethodsFigure 2.6: Illustration of Gauss quadrature points for flux integration2.4 Boundary Treatment MethodsThe boundary control volumes differ from the interior ones in two major aspects;the stencil and the application of boundary conditions. They are also the regionsof a mesh where the error is maximum compared to the interior and hence needto be treated differently. Thus more resources need to be spent for an effectiveimplementation of boundary conditions for an accurate solution. We present dif-ferent methods of boundary treatments applied in our work to assess their effect onimproving accuracy in truncation and discretization error. The Poisson problem issolved using a homogeneous or non-homogeneous Dirichlet boundary condition.The treatment focuses on four different approaches to reduce the error: 1) im-proving gradient calculation by inclusion of the boundary constraints 2) includingmore stencil members to increase the support for least squares method 3) applyingboundary condition by using ghost cells and 4) by using jump terms.2.4.1 Use of Hard ConstraintsEquation 2.19 for least-squares reconstruction includes the mean constraint. For aboundary control volume with four members in the stencil the equation takes theform292.4. Boundary Treatment Methods1 xi yiwi1 wi1x̂i1 wi1ŷi1wi2 wi2x̂i2 wi2ŷi2wi3 wi3x̂i3 wi3ŷi3wi4 wi4x̂i4 wi4ŷi4φ∂φ∂x∂φ∂yi=φ iwi1φ 1wi2φ 2wi3φ 3wi4φ 4(2.23)The gradients calculated by using these set of equations with just the mean con-straint weakly satisfy the boundary condition. This method of boundary treatmentis termed as ’Weak Boundary Condition (WBC)’ and has a jump in the recon-structed boundary solution and the given boundary condition at the Gauss point.The truncation errors are generally higher and best suited for using along with ajump term.2.4.1.1 Strong Boundary ConditionOllivier-Gooch and Van Altena [34] proposed the method of using constraints atthe boundary Gauss points to satisfy the boundary condition exactly. In this way,the flux integral accurately reflects the true boundary condition. The constraint isincluded in the reconstruction matrix and is termed as ’Gauss Constraint’. We havetwo Gauss points on the boundary per control volume, one per edge on either sideof the vertex for calculating second order accurate flux integral. The reconstruc-tion matrix with the constraints at Gauss points G1(xg1,yg1) and G2(xg2,yg2) andboundary condition given by the function fB(x,y) is1 xi yi1 xg1− xi yg1− yi1 xg2− xi yg2− yiwi1 wi1x̂i1 wi1ŷi1wi2 wi2x̂i2 wi2ŷi2wi3 wi3x̂i3 wi3ŷi3wi4 wi4x̂i4 wi4ŷi4φ∂φ∂x∂φ∂yi=φ ifB(xg1,yg1)fB(xg2,yg2)wi1φ 1wi2φ 2wi3φ 3wi4φ 4(2.24)The constraints are directly used for calculating the unknown derivatives andsolution variables for this scheme. This implementation is referred to as ’StrongBoundary Condition (SBC)’ where the boundary condition and the control volumeaverage are sufficient for solution reconstruction in boundary control volumes.2.4.1.2 4-Unknown ModelPreliminary tests revealed the error in the normal derivative (say y direction) to behigher than the tangential derivative. One remedy to fix the error in the normal302.4. Boundary Treatment Methodsderivative is to include a higher order derivative in the normal direction to increasethe accuracy. This method is termed as the 4Unk model for the four unknowns thatare solved in the reconstruction matrix with the mean and boundary constraintsapplied. The set of equations, for a boundary control volume with four neighbors,takes the form:1 xi yi y2i1 xg1− xi yg1− yi (yg1− yi)21 xg2− xi yg2− yi (yg2− yi)2wi1 wi1x̂i1 wi1ŷi1 wi1ŷ2i1wi2 wi2x̂i2 wi2ŷi2 wi2ŷ2i2wi3 wi3x̂i3 wi3ŷi3 wi3ŷ2i3wi4 wi4x̂i4 wi4ŷi4 wi4ŷ2i4φ∂φ∂x∂φ∂y∂ 2φ∂y2i=φ ifB(xg1,yg1)fB(xg2,yg2)wi1φ 1wi2φ 2wi3φ 3wi4φ 4(2.25)The constraints are eliminated using Gauss elimination with pivoting. The resultantsystem of equations are:wi1Ŷ 2i1wi2Ŷ 2i1wi3Ŷ 2i1wi4Ŷ 2i4(∂ 2φ∂y2)i=wi1Φi1wi2Φi2wi3Φi3wi4Φi4 (2.26)where the calculation of the terms Ŷ 2i j and Φi j is straightforward. Equation 2.26 isfirst solved using the least squares method for second derivative in y and the othergradients are calculated later from the constraint equations.2.4.2 Increased Stencil SizeA boundary control volume is at a natural disadvantage because it has fewer neigh-bors when compared to a well surrounded interior CV. This directly implies that thereconstruction matrix for a boundary CV carries less information from its neigh-borhood than its interior counterparts. Experiments were conducted to assess ifthe boundary CVs can have better gradients and in turn better truncation errorby increasing the number of neighbors in the reconstruction stencil. A triangu-lar boundary CV has four neighbors or fewer and a quad boundary CV has threeor fewer compared to their interior counterparts, which have on an average six andfour neighbors respectively. The stencil was increased to add more neighbors. TheFigure 2.7 illustrates an experiment of increasing the stencil size for triangular andquad boundary CV.312.4. Boundary Treatment Methods(a) Triangular boundary CV(b) Quadrilateral boundary CVFigure 2.7: Illustration of stencil increment2.4.3 Use of Single and Two Jump TermsThe use of jump terms has been discussed in Section 2.3.1 where the importance ofα in reducing truncation and discretization errors has been highlighted. Howeverfor boundary control volumes, the edges on the boundary are different for the cal-culation of numerical flux when compared to the edges which are shared by othercontrol volumes. For interior edges the gradient and solution for two control vol-umes are available to calculate the flux. At the boundary however, the gradient andsolution from the reference control volume is the only information available apartfrom the boundary condition.One approach would be to use the jump term for interior edges of a boundarycontrol volume with a SBC implementation for applying boundary condition. Thisapproach is coined as the ’Single Jump Term Method’. In an alternate approach,a jump term instead of SBC at the boundary can also be used. This approach hasbeen shown to significantly reduce the truncation error in Chapters 4 and 5 com-pared to the case where there was no jump term used. For this approach, a newjump coefficient is introduced to account for the jump in satisfying the boundarycondition by solution reconstruction. Thus we will have two jump coefficients, aninterior jump coefficient α and a boundary jump coefficient αB, and hence the im-plementation is called as a ’Two Jump Term Method’. The jump term for boundaryedges takes the form322.4. Boundary Treatment MethodsFigure 2.8: Illustration for jump term coefficients∇φJumpB =αB|riB · nˆF |(φB−φR)nˆF (2.27)where φB and φR are the boundary condition and reconstructed solution at theGauss point B. Experiments were conducted to assess the impact of two jump termapproach versus a single jump term implementation at the boundary, in improvingthe accuracy of truncation and discretization errors. Figure 2.8 illustrates the useof jump terms for interior and boundary dual faces.2.4.4 Ghost Cell MethodsAn intuitive way of applying the boundary condition is to extend the domain be-yond the boundaries to form ghost cells. Ghost cells are attractive in the sensethat they reduce the need for a special treatment of boundary control volumes.Once ghost cell control volume averages are determined using the boundary condi-tion, they can be directly included in the framework of least-squares reconstruction.They have been employed to apply boundary conditions for complex CFD prob-lems to improve accuracy in solution [6] and also for flux-limiting schemes [29]. Abrief description of the ghost cell implementation for different boundary conditionsis provided in the books of Blazek [4] and Leveque [22].Ghost cells are virtual cells outside the physical domain but are associated withthe geometric quantities like volume, area, moments etc, and the physical quantitieslike temperature, pressure and velocity for CFD problems, and in our case thesolution values. For the vertex-centered method, two layers of ghost cells are addedoutside the domain. The solution values are extrapolated from the neighboring non-ghost cell control volumes such that the boundary condition is satisfied by default.332.4. Boundary Treatment MethodsFigure 2.9: Illustration of ghost cell methodAs a result, the reconstruction for boundary control volumes has a well surroundedstencil which implicitly applies the boundary condition. The flux is calculated onthe physical control volume edges at the boundary. Figure 2.9 gives an illustrationof the Gauss points on the boundary and the reflection of cells used to create ghostcells.34Chapter 3Methods of Error AnalysisThe steps to obtaining a finite volume solution to a Poisson problem with differ-ent boundary treatments have been laid out in the previous chapter. The study oftheir effect on reducing truncation error and discretization error require that boththe errors are evaluated consistently and precisely. An interesting way to evaluatetruncation and discretization errors for any PDE is the Method of ManufacturedSolutions (MMS). These were successfully employed for assessing the magnitudeof error, order of convergence and verification of the code [5, 37, 40, 9]. Using themethod of manufactured solutions, the exact solutions are known and the calcula-tion of discretization error is thus straightforward.Truncation error measures the accuracy of the discrete approximation to thedifferential equation. The truncation error analysis approach in this thesis is similarto the work of Jalali [17], where the method is applicable to any arbitrary solution.The Taylor polynomial coefficients associated with the solution and its derivativesare calculated for both the exact and numerical flux integrals and the truncationerror is computed by subtracting the exact flux integral from the numerical cellresidual.The modeled Poisson’s equation is linear in nature and since we are only in-terested in the spatial accuracy of the discretization schemes, its steady state formalone is considered. Truncation error analysis is performed on both the linear al-gebra system and the numerical simulation code ANSLib whereas the calculationof discretization error is only possible on ANSLib. Both the codes are grid trans-parent and are capable of reading mixed meshes which are used in the course ofthis research. ANSLib uses an edge based data structure which allows it to handlestructured, unstructured and hybrid meshes seamlessly.The linear algebra system uses a regular topology uniform mesh which hasbeen perturbed, scaled, sheared, stretched or curved to evaluate the truncation error.The algebra system uses an analytic approach in calculating the truncation errorand is not easily generalizable to include a general unstructured mesh which hasirregular connectivity of five, six or seven incident edges on a vertex. The type ofmesh used for analytic tests is shown in Figure 3.1, which is primarily employed fortruncation error analysis of the pure triangular mesh. A different mesh (see Figure3.3) with quad layers in the bottom half and triangular mesh in the top half is used353.1. Analytic TestsFigure 3.1: Uniform mesh regular stencilfor truncation error analysis of mixed meshes. The regions of concern for this meshare the triangular-quad interface layer control volume, and the boundary quad cells.The mesh is chosen so as to include control volumes for boundary and interiorquad duals, and an interface layer control volume to perform the analysis required.The numerical approach used for truncation error analysis on ANSLib can easilycalculate the truncation error for any arbitrary stencil in a general unstructuredmesh regardless of the partial differential equations for linear problems.3.1 Analytic TestsThe analytical approach for truncation error is explained using the Figures 3.1and 3.3 for the interior, interface and boundary layer control volumes. The meshis composed of equilateral triangles and quadrilaterals which have side length h.Truncation error analysis for interior control volume is performed on CVs whosefirst and second neighbors are not boundary control volumes. However, the inter-face layer where a triangular cell and quadrilateral cell adjoin is treated differentlyfrom an interior control volume. The control volumes on the boundary and theirimmediate neighbors are considered crucial for the purpose of reducing inaccuracyin boundary condition implementation and reducing truncation and discretizationerrors. The methodology itself reveals the differences in order of accuracy not justbecause of the boundary condition but also due to non-uniformity of the control363.1. Analytic Testsvolume shape.The Poisson’s equation is the model equation for diffusion problems. The com-putation of diffusive fluxes is a step by step procedure where the control volumeaverages are first conserved and the property is used for the calculation of least-squares reconstruction. Alternatively, Green-Gauss method is used for evaluatingthe gradients. Finally the gradients are averaged on the dual edges for the flux cal-culation. The benefit of analytic tests lies in a quick and easy evaluation of orderof accuracy for gradient, solution and flux integral calculation.φ i = φ |i +∂φ∂x∣∣∣∣ixi +∂φ∂y∣∣∣∣iyi +O(h2) (3.1)Consider Equation 3.1 which gives the control volume average for the vertexi. For cell-centered schemes, the first moments xi and yi are both zero and bydefault the average is represented with second order accuracy at the centroid ofthe cell. However, for a vertex-centered scheme, the first moments are zero for auniform mesh but not when the mesh is perturbed randomly. This difference plays acrucial role for analytic test accuracy. The least-squares method of reconstructionfor an interior and boundary CV are presented here. As the stencil is of regularconnectivity, the interior and boundary control volumes have six and four neighborsrespectively.wi2(x̂i2− xi) wi2(ŷi2− yi)wi3(x̂i3− xi) wi3(ŷi3− yi)wi4(x̂i4− xi) wi4(ŷi4− yi)wi5(x̂i5− xi) wi5(ŷi5− yi)wi6(x̂i6− xi) wi6(ŷi6− yi)wi7(x̂i7− xi) wi7(ŷi7− yi)∂φ∂x∂φ∂yi=wi2(φ 2−φ i)wi3(φ 3−φ i)wi4(φ 4−φ i)wi5(φ 5−φ i)wi6(φ 6−φ i)wi7(φ 7−φ i)(3.2)The gradient of the control volume i = 1 as depicted in Figure 3.1 is calculatedusing Equation 3.2. The uniform nature of the mesh can be seen in the calculationof the moments and weights and this results in perfect error calculation for accurateresults. The gradient calculated is expressed in terms of the neighboring controlvolume averages as∇˜φ 1 =∂˜ φ∂x∂˜ φ∂y1=− 16h(2φ 7 +φ 2 +φ 6−2φ 4−φ 3−φ 5)√36h (φ 2 +φ 3−φ 5−φ 6) (3.3)373.1. Analytic TestsWhen the neighboring control volume averages are expressed as a Taylor seriesexpansion about the reference control volume i = 1, the order of accuracy for thegradients can be easily assessed. The results reveal that the gradients are secondorder accurate for the reference control volume as shown in Equation 3.4. Thisorder of accuracy is one more than the expected and observed order of accuracyfor unstructured meshes in general.∇˜φ 1 =∂˜ φ∂x∂˜ φ∂y1=23h2144(∂ 3φ∂x3 +∂ 3φ∂x∂y2)+ ∂φ∂x23h2144(∂ 3φ∂y3 +∂ 3φ∂x2∂y)+ ∂φ∂y (3.4)The boundary control volume has four neighbors as shown in Figure 3.2 andthe neighbors differ in the moment calculation and also the type of boundary con-dition applied. The least-squares system for the reference control volume i = 1in Figure 3.2, using weak boundary condition enforcement is given in Equation3.5. Evidently, there is less information from the neighbors for calculating the gra-dients. The gradient when expressed in terms of the neighboring control volumeaverages iswi2(x̂i2− xi) wi2(ŷi2− yi)wi3(x̂i3− xi) wi3(ŷi3− yi)wi4(x̂i4− xi) wi4(ŷi4− yi)wi7(x̂i7− xi) wi7(ŷi7− yi)∂φ∂x∂φ∂yi=wi2(φ 2−φ i)wi3(φ 3−φ i)wi4(φ 4−φ i)wi7(φ 7−φ i) (3.5)∇˜φ 1 =∂˜ φ∂x∂˜ φ∂y1=− 15h(2φ 7 +φ 2−2φ 4−φ 3)9√320h(φ 2 +φ 3−2φ 1) (3.6)As expected, the accuracy for boundary control volumes for gradient calcula-tion is O(h). The trend observed for the mixed cell mesh is similar for both theinterior and boundary control volumes in terms of the order of accuracies for gra-dient calculation.After evaluating the gradients using the least-squares method, there is a needfor averaging the gradients using Table 2.1 for face gradient calculation. An exam-ple is provided for the face gradient calculation at the interface of control volumesi = 1 and j = 7 using a simple arithmetic averaging scheme for the mesh shown inFigure 3.2.383.1. Analytic TestsFigure 3.2: Uniform mesh boundary stencil∇˜φ 17 =∂˜ φ∂x∂˜ φ∂y17=110h(−2φ 7−2φ 18−φ 19 +φ 3 +2φ 4 +2φ 1)−9√340h (2φ 7 +2φ 1−φ 19−φ 3−2φ 2) (3.7)The gradient at the interface of the control volumes i = 1 and j = 7 lacksthe symmetry which is usually observed for the interior control volumes. Theface gradients and the control volume gradients are all observed to be first orderaccurate. Important to note is that the averaging of gradient is not required onthe boundary dual edges as the boundary Gauss point lies completely within thereference control volume. Once the face gradients are calculated for every Gausspoint on the control volume boundary, the flux for each dual edge is calculated bythe dot product of the normal unit vector and the face gradients.∇˜2φi =1ACVi˛∂Vi∇φF · nˆds (3.8)The flux integral calculation for a dual edge is explained in Section 2.3.2 wherethe flux is multiplied by the appropriate weight which is the edge length for secondorder accuracy. Once the flux integral for the control volume is computed, weobtain an expression for the numerical flux integral in terms of the control volumeaverages. The numerical flux integral obtained using least-squares reconstructionis a linear combination of 19 control volumes for control volume i = 1 in the Figure3.1. Using Green-Gauss reconstruction gives a much more compact stencil for thenumerical flux integral.∇˜2φ i−LSM =−19h2(6φ 1 +7∑p=2φ p−19∑q=8φ q) (3.9)393.1. Analytic TestsFigure 3.3: Uniform mixed cell mesh regular stencil∇˜2φ i−GG =−23h2(6φ 1−7∑p=2φ p) (3.10)A strong decoupling is observed when the numerical flux integral is calculatedfor the quad interior cells using least-squares reconstruction. A discussion of thedecoupling and its effects and method to prevent decoupling can be found in thework of Haselbacher and Blazek [15]. The stencil for the flux integral is again a lin-ear combination of 5 control volumes and doesn’t involve the immediate neighborsof the reference control volume i = 25 shown in Figure 3.3.∇˜2φ i =14h2(φ 23 +φ 27 +φ 35 +φ 15−4φ 25) (3.11)Decoupling can lead to spurious solution modes and large truncation errors,403.1. Analytic Testsboth of which are undesirable. By implementing the edge derivative correction forface-gradient calculation on the dual edges (see Section 2.2.1, Equation 2.4), thenumerical flux integral is recalculated using gradients calculated by least-squaresmethod. The resultant stencil observed is a compact one involving the referenceand the adjacent control volumes as shown in Figure 3.3. The numerical flux inte-gral of an interior control volume in the quad region of the mixed mesh, calculatedusing least-squares method of reconstruction and edge correction for face gradi-ents, is∇˜2φ i−EC =1h2(φ 24 +φ 26 +φ 30 +φ 20−4φ 25) (3.12)Now that the remedies to prevent instabilities in the interior regions of the meshare described, we proceed onto the evaluation of the flux integral for a boundarycontrol volume. The numerical flux integral for the boundary control volume, asshown in Figure 3.2, is evaluated using least-squares reconstruction. A compar-ison of the stencil for interior and boundary flux integral can be drawn from theexpression for the numerical flux integral for boundary given in equation below∇˜2φ i =1h2((a1φ 1 +a2(φ 2 +φ 3)+a4(φ 4 +φ 7)+ (3.13)a8(φ 8 +φ 10)+a9φ 9 +a11(φ 11 +φ 19)+a12(φ 12 +φ 18))a1 = 1.103 a2 =−1.189a9 = 0.319 a4 =−0.104 a8 = 0.271a11 = 0.178 a12 = 0.133Table 3.1: Coefficients for boundary numerical flux integralThe coefficients a′is for Equation 3.13 are given in the Table 3.1. Even for aperfectly uniform mesh, the coefficients and a one-sided stencil hint at the lack oferror cancellation required for higher order truncation error.Having presented the stencil sizes and behavior of the schemes, we shall nowdiscuss the evaluation of the analytic flux integrals, both exact and numerical. Theevaluation of the error in the numerical flux integral requires the control volumeaverages in the stencil to be expressed as a Taylor series expansion about the ref-erence control volume, using the Equation 2.16. For the purpose of demonstration,we calculate the truncation error only for the interior control volume i = 1 in Figure3.1. The control volume averages are now replaced in Equation 3.9 by their Taylor413.1. Analytic Testsseries expansions about the reference control volume, which after some simplifica-tion yields∇˜2φ 1 =∂ 2φ∂x2∣∣∣∣1+∂ 2φ∂y2∣∣∣∣1+41h2144(∂ 4φ∂x4∣∣∣∣1+2∂ 4φ∂x2y2∣∣∣∣1+∂ 4φ∂y4∣∣∣∣1)(3.14)To evaluate the truncation error, the exact flux integral has to be calculated tocompare with the numerical one. The Laplacian operator is performed over thesolution, which is expressed in a Taylor series expansion about its reference pointas given in the Equation 2.7, and averaged over the corresponding control volumeto obtain the exact flux integral. The Laplacian of the solution is given as∇2φ i =∂ 2φ∂x2∣∣∣∣i+∂ 2φ∂y2∣∣∣∣i+∂ 3φ∂x3∣∣∣∣i(x− xi)+∂ 3φ∂x2∂y∣∣∣∣i(y− yi)+ (3.15)∂ 3φ∂x∂y2∣∣∣∣i(x− xi)+∂ 3φ∂y3∣∣∣∣i(y− yi)+∂ 4φ∂x4∣∣∣∣i(x− xi)22+∂ 4φ∂x3y∣∣∣∣i(x− xi)(y− yi)+∂ 4φ∂x2y2∣∣∣∣i((x− xi)22+(y− yi)22)+∂ 4φ∂xy3∣∣∣∣i(x− xi)(y− yi)+∂ 4φ∂y4∣∣∣∣i(y− yi)22+ · · ·The next step is to take the area average over the control volume to obtain thegeneral form of exact flux integral for arbitrary mesh.∇2φ i =∂ 2φ∂x2∣∣∣∣i+∂ 2φ∂y2∣∣∣∣i+(∂ 3φ∂x3∣∣∣∣i+∂ 3φ∂x∂y2∣∣∣∣i)xi + (3.16)(∂ 3φ∂x2∂y∣∣∣∣i+∂ 3φ∂y3∣∣∣∣i)yi +∂ 4φ∂x4∣∣∣∣ix2i2+∂ 4φ∂x3y∣∣∣∣ixyi +∂ 4φ∂x2y2∣∣∣∣i(x2i2+y2i2)+∂ 4φ∂xy3∣∣∣∣ixyi +∂ 4φ∂y4∣∣∣∣iy2i2+ · · ·Equation 3.16 gives the exact flux integral, irrespective of whether it is a quador a triangular control volume, interior or boundary. The moments for an interiorcontrol volume for a uniform triangular mesh are substituted in Equation 3.16 toobtain the exact flux integral.423.2. Numerical Tests∇2φ 1 =∂ 2φ∂x2∣∣∣∣1+∂ 2φ∂y2∣∣∣∣1+ (3.17)5h2144(∂ 4φ∂x4∣∣∣∣i+2∂ 4φ∂x2y2∣∣∣∣i+∂ 4φ∂y4∣∣∣∣i)+ · · ·The truncation error for interior control volumes is given in Equation 3.18which is observed to be second order accurate. If the vertex locations are per-turbed randomly, this accuracy is impossible to reproduce and the truncation errorstops converging with decreasing mesh size and has the order O(1).∇˜2φ 1−∇2φ 1 =h24(∂ 4φ∂x4∣∣∣∣i+2∂ 4φ∂x2y2∣∣∣∣i+∂ 4φ∂y4∣∣∣∣i)+O(h4)(3.18)The truncation error depends on the type of reconstruction used, face gradientcalculation, use of jump terms and edge correction and also the type of boundarycondition implementation. The influence of different schemes and methods ontruncation error and discretization error are studied in detail in Chapters 4 and 5.3.2 Numerical TestsThe analytic tests were designed for meshes with regular topology to give an ideaabout the behavior of the discretization schemes and boundary treatment methodson a topologically regular mesh. The analysis is best performed by the computerlinear algebra system developed which replaces the laborious hand calculationsrequired otherwise. A similar analysis for truncation error on general unstruc-tured meshes is considered impractical as it would entail mathematically intenseoperations on every control volume. An alternative approach is explained in thefollowing section to calculate truncation error based on the methods described byJalali [17] for a similar analysis of cell-centered schemes in the interior regions.A brief discussion on calculation of the discretization error using the Method ofManufactured Solutions is also presented in Section 3.2.2.3.2.1 Truncation ErrorThe numerical approach for calculation of truncation error is applicable for generalunstructured meshes and is independent of the nature of partial differential equa-tions solved. The approach exploits the fact that the discrete flux integral can beexpressed as a linear combination of the control volume averages for an arbitrary433.2. Numerical Testsunstructured mesh stencil. The discrete flux integral in a general form is expressedas∇˜2φ i =∑j∂Ri∂φ jφ j (3.19)where the summation is over the control volumes in the stencil {Vj}i of the ref-erence control volume i. ∂Ri∂φ jis one row of the global flux Jacobian. ANSLib iscapable of calculating the global flux Jacobian explicitly, for schemes up to fourthorder accuracy. Detailed information on the calculation of the explicit flux Jaco-bian can be found in the works of Michalak [29]. The procedure pertaining to thePoisson problem is described here. The analytic Jacobian is explicitly representedas∂R∂φ=∂FluxInt∂Var=∂FluxInt∂Flux∂Flux∂RecSol∂RecSol∂RecCoef∂RecCoef∂Var (3.20)in which FluxInt is the discrete flux integral, Flux is the numerical flux calcu-lated at the Gauss point, RecSol are the reconstructed solutions at the Gauss points,RecCoef are the coefficients of the reconstructed polynomial obtained by the least-squares method and Var is the control volume average of the unknown variable.To compute the Jacobian, the following procedure is used for each Gauss point foreach of the two adjacent control volumes:1. ∂RecSol∂Var =∂RecSol∂RecCoef∂RecCoef∂Var is first computed. The term∂RecCoef∂Var is simply thepseudo-inverse of the reconstruction matrix precomputed in the least-squaressystem of Equation 2.17 using SVD. The numerical flux is always a func-tion of the solution and/or solution gradient at the Gauss points and requiresreconstruction of the same variables at the Gauss points. Thus RecSol is de-fined as(φ ∂φ∂x∂φ∂y)T. The term ∂RecSol∂RecCoef is essentially a geometric termthat only depends on the location of the Gauss point, as RecSol is expressedas a Taylor polynomial about the reference point.2. In the Poisson problem, the Flux is defined as the dot product of the facegradients and the unit normal vector at the Gauss points. Thus ∂Flux∂RecSol de-pends on the choice of discretization and use of jump terms or edge derivativecorrection to calculate the flux. So if only arithmetic averaging is used for443.2. Numerical Testscalculating the face gradient, the term evaluates as∂Flux∂RecSolL=∂Flux∂RecSolR=(0 −12 nˆx −12 nˆy)(3.21)3. ∂FluxInt∂Flux depends on the Gauss quadrature formula. It is calculated by usingthe appropriate Gauss integration weight and the length of the correspondingedge.4. After calculating all the terms, the ∂FluxInt∂Var is obtained by combining them.Once the explicit Jacobian ∂R∂φ is calculated using the above procedure1, the steps toexpress the flux integral in terms of the control volume averages φ j is performed.Then, the control volume averages for the neighboring control volumes are re-placed by their Taylor series expansions about the reference control volume interms of its solution and derivatives, using Equation 2.16. In calculating the sum,the terms corresponding to the Taylor coefficients are collected so as to attain aform seen in numerical flux integral given by the Equation 3.14. The exact fluxintegral is calculated again by computing the average of the flux integral over theentire control volume as given in Equation 3.16. The difference of the numeri-cal flux integral from the exact flux integral simply gives the truncation error foreach control volume. Subject to different boundary treatments, jump terms or dis-cretization methods, the truncation error is evaluated over the family of meshes toperform convergence and accuracy analysis for various schemes.3.2.2 Discretization ErrorDiscretization error is usually difficult to measure for general flow scenarios, es-pecially when the exact solution is unknown. Roy [36] has reviewed differentdiscretization error estimators and noted that these estimations are only reliablewhen the numerical solution is in asymptotic range, the verification of which re-quires at least three systematically refined meshes. It has also been observed [12]that for regular grids, the convergence of truncation error is an accurate indicatorof convergence for discretization errors provided the discrete boundary conditionsare adequate. On irregular grids, though, the truncation error convergence is nota good predictor for discretization error convergence. Design order discretizationerror convergence has been observed in cases where the truncation error converges1An alternative way to calculate the Jacobian is the “finite difference Jacobian” where the vari-ables concerned are perturbed individually and the change in the flux integral function is measuredto calculate the Jacobian. This method is slower than explicit Jacobian method and is used for trun-cation error estimation using the Green-Gauss reconstruction scheme453.2. Numerical Testsat a lower than expected order or doesn’t converge at all. Thus its necessary to es-timate the discretization error accurately and study its relationship with truncationerror.Given the simple nature of Poisson’s equation, we use the Method of Manu-factured Solutions with three exact solutions and the associated source terms tosolve the problem numerically. The discrete solution obtained is then comparedwith the exact solution to evaluate the discretization error. To verify the conver-gence order of discretization error, we use four systematically refined meshes foreach sequence of mesh type. The relationship between the truncation error and dis-cretization error is explored and the effect of different discretization schemes andboundary treatment methods are analyzed for reducing discretization error as well.An attempt is also made to realize the conditions under which the truncation errorand discretization error are minimized simultaneously. The effects of minimizingtruncation error and its influence on discretization error are also studied.46Chapter 4Triangular Mesh ResultsTo achieve the least truncation error at the boundary with reduced discretizationerror, we proceed in a step-by-step manner. The schemes presented in Section 2.4are subjected to analytic tests described in Section 3.1 and those that perform wellare implemented in the unstructured solver. Numerical tests are then performedto confirm the behavior observed in analytic tests and to generalize the results forisotropic unstructured meshes. Using analytic tests, we exploit the algebraic na-ture of results to estimate the order of convergence for truncation error for variousdiscretization schemes and boundary treatment methods, and the order of error ingradient calculation for different reconstruction methods. This gives a predictionfor the behavior to be expected in the unstructured solver, paving the way to verifyand solve the problem accurately.The analytic tests for triangular meshes are performed on the meshes shownin Figures 3.1 and 3.2 for interior and boundary control volumes. Using the sametopology, the mesh is perturbed in many systematic and random ways to test fordifferent cases. The schemes performing poorly are eliminated from considerationfor numerical tests where the mesh is fully unstructured. The following sectionsdetail the findings of analytic tests followed by numerical tests performed on purelytriangular meshes.4.1 Analytic Test AnalysisAnalytic tests provide an excellent analysis method to verify if the reconstruc-tion, face-gradient calculation and flux integration behave the way as expected. Ashighlighted earlier, they are used to troubleshoot errors in implementation and alsopredict the behavior of schemes on general unstructured meshes. To ensure thatanalytic results are applicable for general scenarios, they are tested on five topo-logically regular meshes, which are created by distorting the vertices of a uniformtriangular mesh in a specific manner and are shown in Figure 4.1. The referenceposition for CV1 is preferably but not necessarily kept at the origin. The family ofdistorted meshes are described here:1. Sheared meshes are distortions of a uniform triangular mesh where the y-474.1. Analytic Test Analysiscoordinates remain unchanged but the x-coordinates are sheared. The x-coordinates in a horizontal line are all moved by a constant factor which isa function of their distance from the reference x-axis. If each horizontal lineis numbered with the reference x-axis as 0th row and increasing upwards,the new x-coordinates are then obtained by adding a constant amount bmhwhere b > 0 is a constant shearing factor and m is the row number. The newx-coordinates are such that the slanted lines remain straight and parallel.2. Scaled meshes have the x-coordinates of a uniform mesh unperturbed but they-coordinates are scaled by a constant factor of k > 0.3. Stretched meshes also have the x-coordinates the same as that of uniformmesh. However, the y-coordinates appear stretched as one moves away fromthe reference x-axis. A stretching factor s is defined to determine the heightof each row of vertices. If the horizontal rows are numbered, the height ofthe triangle between rows m+ 1 and m would be given by√32 hs−m, where0 < s < 1.4. Curved meshes are such that the mesh vertices are all transformed to a radialcoordinate system. The vertices with the same y-coordinate are mapped ontoa circular arc with the same radial coordinate while the vertices with thesame x-coordinate are mapped on the radial line with the same polar angle θ .These meshes are characterized by the constant ratio (R) of radius betweenthe consecutive layers of mesh and by the ratio a = hr0 for the mesh.5. Random meshes have both the coordinates of the vertices perturbed ran-domly, by ρxh and ρyh in the x and y directions respectively. The ran-dom variations ρx,ρy ∈ [−0.1,0.1] such that the resultant mesh has non-overlapping cells after the operation.4.1.1 Reconstruction AccuracyA first order gradient is expected for a second order solution reconstruction usingthe least-squares and the Green-Gauss method. As discussed earlier, the gradientsfor a control volume calculated by the least-squares method involve solving thesystem of equations required to minimize the error in satisfying the mean in neigh-boring control volumes. The Green-Gauss method calculates the gradient for theprimal cell instead of the control volume itself. Both the methods calculate thegradients as a linear combination of control volume averages in the stencil of thereference CV.484.1. Analytic Test Analysis(a) Perfect mesh (b) Sheared mesh (c) Scaled mesh(d) Curved mesh (e) Stretched mesh (f) Random meshFigure 4.1: Family of distorted meshesWe have already discussed the gradient accuracy for interior control volumesfor a perfect mesh which has second order accuracy for reconstruction. We nowobserve the gradient accuracy for control volumes in the stencil of the referenceboundary control volume i = 1 in the Figure 3.2, using the least-squares recon-struction. The evaluation of the gradients is explained in Section 3.1 on page 38and the same results are presented here:∇˜φ 1 =∂˜ φ∂x∂˜ φ∂y1=− 15h(2φ 7 +φ 2−2φ 4−φ 3)9√320h(φ 2 +φ 3−2φ 1) (4.1)=∂φ∂x∣∣∣1+ 11√354∂ 2φ∂x∂y∣∣∣1h+O(h2)∂φ∂y∣∣∣1+ 9√380(∂ 2φ∂x2∣∣∣1+3 ∂2φ∂y2∣∣∣1)h+O(h2)494.1. Analytic Test AnalysisEquation 4.1 gives the gradient calculated for a boundary control volume of theuniform mesh. The control volume averages in the above gradient calculation areexpressed as a Taylor series expansion about the reference control volume i = 1using the Equation 2.16, to determine its accuracy. The Taylor series expansion ofthe solution variable about the control volume reference point is given as:φi(x,y) = φ |i +∂φ∂x∣∣∣∣i(x− xi)+∂φ∂y∣∣∣∣i(y− yi)+ (4.2)12∂ 2φ∂x2∣∣∣∣i(x− xi)2 +∂ 2φ∂x∂y∣∣∣∣i(x− xi)(y− yi)+12∂ 2φ∂y2∣∣∣∣i(y− yi)2 + ...To calculate the exact gradient, the Del operator is performed over the aboveequation. The resultant is averaged over the control volume to calculate the exactgradients as shown:∇φ i =∂φ∂x∂φ∂yi=∂φ∂x∣∣∣i+ ∂2φ∂x2∣∣∣ixi +∂ 2φ∂x∂y∣∣∣iyi + · · ·∂φ∂y∣∣∣i+ ∂2φ∂y2∣∣∣iyi +∂ 2φ∂x∂y∣∣∣ixi + · · · (4.3)The difference in the exact and discrete gradient calculation gives the cell gra-dient error and one can verify the order of accuracy for reconstruction from thiserror. The general form of this error is:∇˜φ 1−∇φ 1 =(ξx ∂φ∂x∣∣∣1+ξy ∂φ∂y∣∣∣1)+(∂ 2φ∂x2∣∣∣1ξxx + ∂2φ∂x∂y∣∣∣1ξxy + ∂2φ∂y2∣∣∣1ξyy)h+O(h2)(ηx ∂φ∂x∣∣∣1+ηy ∂φ∂y∣∣∣1)+(∂ 2φ∂x2∣∣∣1ηxx + ∂2φ∂x∂y∣∣∣1ηxy + ∂2φ∂y2∣∣∣1ηyy)h+O(h2)(4.4)504.1. Analytic Test AnalysisMesh ξx ξy ξxx ξxy ξyyLeast-SquaresUniform 0 0 0 0.1283 0Random 0 0 0.10779 0.14341 0.00208Green-GaussUniform 0 0 0 0.64150 0Random −0.07744 −0.00816 −0.46849 −0.00608 −0.00191(a) Error coefficients of ∂φ∂xMesh ηx ηy ηxx ηxy ηyyLeast-SquaresUniform 0 0 0.19486 0 0.36004Random 0 0 0.18030 −0.00084 0.37131Green-GaussUniform 0 −0.25926 0.14434 0 0.20849Random −0.10161 −0.24883 −0.09825 −0.35668 0.21764(b) Error coefficients of ∂φ∂yTable 4.1: Gradient error coefficients for boundary control volumesMesh ξx ξy ξxx ξxy ξyyLeast-SquaresUniform 0 0 0 0 0Random 0 0 0.09227 0.070601 −0.00856Green-GaussUniform 0 0 0 0.86603 0Random 0.16051 −0.01405 −0.01254 1.04424 0.02339(a) Error coefficients of ∂φ∂xMesh ηx ηy ηxx ηxy ηyyLeast-SquaresUniform 0 0 0 0 0Random 0 0 0.03151 −0.01332 0.12935Green-GaussUniform 0 0 0.14434 0 0.43301Random 0.00902 0.04706 0.12362 0.02548 0.47531(b) Error coefficients of ∂φ∂yTable 4.2: Gradient error coefficients for interior control volumes514.1. Analytic Test AnalysisTables 4.1 and 4.2 give the error coefficients ξ and η for gradients as in Equa-tion 4.4 for a random mesh with fixed perturbations. The least-squares method ofreconstruction performs the way as expected with the error coefficients for the zeroorder error being zero (ξx = ξy = 0 and ηx = ηy = 0). The Green-Gauss methodfollows this trend only for perfect and non-randomly perturbed meshes, in the in-terior control volumes. The Green-Gauss method as discussed in the Section 2.2.1yields a zero-order gradient for general random meshes and is thus undesirable.The error metrics for interior control volumes are given for comparison andone can notice that the gradient for the least-squares method is one order better foruniform interior control volumes than expected and has a minimum of first orderaccuracy for all other cases, for interior and boundary control volumes. It is alsoimportant to note that the error coefficients for gradients in the normal direction yto the boundary is twice the error in the tangential direction x to the boundary. Thisresult served as a motivation to improve the gradient in the normal direction and isthe reason for the implementation of the method described in Section 2.4.1.2An analysis of the results also reveal that the Green-Gauss method is underperforming for random meshes in the interior and in general for boundary con-trol volumes with the gradient being zero order. Experiments revealed that sucha discrepancy arises because the control volume average is used for the gradientcalculation but the reference locations used are the vertex points. This has theimplication of having one order less accurate gradients as the vertex locations rep-resent the control volume average with first order accuracy and not second orderaccuracy; the centroids of the control volumes do so by default. So instead of usingthe vertex locations, the centroid locations where used for calculating the gradientsand upon verification, the gradient accuracy was found to be consistent with the ex-pected order. The flux integration however, is performed on the same median dualfor the corrected Green-Gauss method. To implement this correction, one woulduse the centroid locations of the corresponding control volumes as (x j,y j) insteadof the vertex locations in the Equation 2.3 to calculate the gradients. The resultsfor corrected Green-Gauss method for interior and boundary control volumes is asshown in the Table 4.3.The error coefficients for both the Green-Gauss methods behave identically forinterior uniform mesh as it was expected. Clearly, the gradient using the referencelocation correction for Green-Gauss method is first order accurate, as can be seenfrom the Table 4.3. The gradient error coefficients are also observed to be betterfor the least-squares method when compared to the Green-Gauss method, with orwithout correction.524.1. Analytic Test AnalysisMesh ξx ξy ξxx ξxy ξyyInterior control volumeUniform 0 0 0 0.86603 0Random 0 0 −0.00683 0.90451 0.025258Boundary control volumeUniform 0 0 0 0.64150 0Random 0 0 −0.00508 0.65342 0.01655(a) Error coefficients of ∂φ∂xMesh ηx ηy ηxx ηxy ηyyInterior control volumeUniform 0 0 0.14434 0 0.43301Random 0 0 0.11828 0.01766 0.45525Boundary control volumeUniform 0 0 0.19486 0 0.36004Random 0 0 0.16576 0.01851 0.37194(b) Error coefficients of ∂φ∂yTable 4.3: Gradient error coefficients using corrected Green-Gauss method4.1.2 Flux Integral AccuracyThe truncation error for the Poisson’s equation is calculated as described in Section3.1 where the exact flux integral is subtracted from the computed residual for thecontrol volume. The Table 4.4 gives the results for order of truncation error withdifferent schemes and reconstruction methods, for a variety of distorted meshesfor interior and boundary control volumes. The table shows the zero order trunca-tion error obtained for boundary control volumes using the least-squares methodand the corrected Green-Gauss (GGC) method. The flux integral for the least-squares method (LSM) is calculated using both the arithmetic averaging (AA) andlinear interpolation (LI) and the order of convergence was found to be identical.The Green-Gauss method without correction (GG), however, has negative order ofconvergence for truncation errors with curved and random meshes. The order oftruncation error for LSM and corrected Green-Gauss method is in agreement withthe results of Diskin and Thomas [9]. Diskin and Thomas [12] also highlight thatlower order of convergence or lack of convergence doesn’t necessarily indicate thatthe solution doesn’t converge. In fact, the solution is found to converge for all thereconstruction methods for isotropic meshes, as shown in the numerical results in534.1. Analytic Test AnalysisSchemeBoundary MeshUniform Sheared Scaled Stretched Curved RandomLSM O(1) O(1) O(1) O(1) O(1) O(1)GG O(1) O(1) O(1) O(1) O(1/h) O(1/h)GGC O(1) O(1) O(1) O(1) O(1) O(1)(a) Boundary control volumesSchemeInterior MeshUniform Sheared Scaled Stretched Curved RandomLSM O(h2) O(h2) O(h2) O(1) O(1) O(1)GGC O(h2) O(h2) O(h2) O(1) O(1) O(1)(b) Interior control volumesTable 4.4: Truncation error order comparison for flux integralSection 4.2.∇˜2φi =1ACVi˛∂Vi∇φF · nˆds (4.5)The order of convergence for truncation error can be shown to depend on theaccuracy of the gradient provided the flux integral quadrature formula is of ade-quate accuracy. By careful inspection of Equation 4.5, we observe that the surfacearea ACVi and s are proportional to h2 and h, respectively. Hence if the accuracyfor the gradient is first order, the flux integral has a zero order accuracy, in the ab-sence of error cancellation. However, in case of the normal Green-Gauss method,the gradient is zero order and thus we observe a negative order of convergence fortruncation errors. There are cases where the zero order terms of the flux integralcancel out resulting in a higher order truncation error. This is most commonly ob-served for interior control volumes when the uniform, scaled and sheared mesheswere used. Such a cancellation was never observed for boundary control volumesbecause they lack a well surrounded stencil like the interior control volumes.The general form of truncation error, based on the expected order of accuracyfrom Table 4.4, can be written as:544.1. Analytic Test Analysis∇˜2φ i−∇2φ i =(λxx∂ 2φ∂x2∣∣∣∣i+λxy∂ 2φ∂x∂y∣∣∣∣i+λyy∂ 2φ∂y2∣∣∣∣i)+O(h) (4.6)The truncation error can be calculated if the exact solution and the error coef-ficients are known. However, the problem of minimizing will then be dependenton the nature of problem solved. The error coefficients are the property of thestencil {Vj}i specifically and generally depend on the mesh. Thus the truncationerror can be minimized indirectly by minimizing the error coefficients so that theresults will be applicable to any solution. A truncation error metric is used, similarto the one used by Jalali [17] to compare different schemes in their effectiveness toreduce truncation error. The error metric is the L2 norm of the error coefficients,as shown:Ei =(λ 2xx +λ 2xy +λ 2yy) 12 (4.7)4.1.3 Boundary Treatment MethodsThe error coefficients change with the change in schemes used for flux calcula-tion, boundary treatment methods and especially with the reconstruction method.For triangular meshes, there was less scope to experiment and improve the resultswith the Green-Gauss method compared to the least-squares method. Using least-squares reconstruction, we can either try to increase the accuracy of the gradients,increase the stencil size or implement the boundary conditions differently. Methodslike SBC and 4-Unknown models are aimed at improving the accuracy of gradi-ents. Increasing the stencil size is also expected to improve the gradient accuracywhereas ghost cells and use of jump terms are expected to implement boundarycondition effectively.Reflecting on Section 4.1.1 and the gradient error in Equation 4.4, error metricsare defined for the gradients to compare the accuracy of different methods.EXi =(ξ 2xx +ξ 2xy +ξ 2yy) 12 (4.8)EYi =(η2xx +η2xy +η2yy) 12 (4.9)Figure 4.2 presents the results for the use of different schemes to improve thegradient at the boundary. For weak boundary condition enforcement, the error inthe normal gradient is almost twice the error in the tangential gradient. The er-rors in the gradients are comparable for both the SBC model and the 4-Unknownmodels. Stencil size was increased to provide sufficient support for reconstruc-tion of gradients. However, the increase in stencil size seems to affect the gradient554.1. Analytic Test AnalysisFigure 4.2: Gradient error metrics for different schemesaccuracy adversely as shown in the Figure 4.2 where error metrics for gradientsin x (SC-GradX) and y (SC-GradY) direction are higher for all schemes. This isprobably because the data from distant control volumes is less relevant for recon-struction.A jump term is added to the flux in the face normal direction. An interestingresult is observed when a jump term is used for interior control volume on a uni-form mesh. The truncation error is second order accurate when the least-squaresmethod is used. By using a single jump term, the truncation error takes the analyticform as shown:∇˜2φ i−∇2φ i =3h216(43−α)(∂ 4φ∂x4∣∣∣∣i+2∂ 4φ∂x2∂y2∣∣∣∣i+∂ 4φ∂y4∣∣∣∣i)+O(h4) (4.10)The jump coefficient modifies the error coefficients and in this case, it can increasethe order of accuracy if the jump coefficient is chosen to be α = 43 . The resultanttruncation error is observed to be of the order of O(h4) by the use of single jumpterm for uniform interior control volumes. A gain in order of accuracy was impos-sible to achieve for a random mesh; the truncation error was however found to bereduced significantly. The truncation error for a boundary control volume, usingsingle jump term is given as564.1. Analytic Test Analysis∇˜2φ i−∇2φ i =((λxx−Cxxα)∂ 2φ∂x2∣∣∣∣i+(λyy−Cyyα)∂ 2φ∂y2∣∣∣∣i)+O(h) (4.11)where Cxx,Cyy are coefficients for α in the truncation error for flux integral. Oneintuitive way to achieve a higher order truncation error is to introduce another con-trol variable, like the jump term α in the flux calculation. Using a boundary jumpterm is one way to introduce the required control parameter αB through which ahigher order truncation error can be achieved. The truncation error with two jumpterms takes the form∇˜2φ i−∇2φ i =((λxx−Cxxα−BxxαB)∂ 2φ∂x2∣∣∣∣i+(λyy−Cyyα−ByyαB)∂ 2φ∂y2∣∣∣∣i)+O(h) (4.12)Ei− jump =((λxx−Cxxα−BxxαB)2 +(λyy−Cyyα−ByyαB)2+(λyy−Cyyα−ByyαB)2) 12 (4.13)where Bxx,Byyare the corresponding coefficients for αB in the truncation error forflux integral. By a careful selection of the jump coefficients it is possible to havea higher order truncation error for uniform boundary control volumes. However,the stability of the solution is a concern for the jump coefficients obtained for thiscase and hence, the discussion is left for future work. The corresponding errormetric when jump terms are used for a regular boundary control volume is givenin Equation 4.13 which predicts the behavior of the error metric to be parabolicwith the variations in jump coefficients. The jump coefficient αB appears only forboundary control volumes and hence its effect is local.The performance of the boundary treatment methods are compared for a ran-dom mesh using the error metric for truncation error in the flux integral. The resultsare presented in the Figure 4.3. Except for the corrected Green-Gauss method, allother schemes employ the least-squares method of reconstruction. The error metricis presented on the primary axis whereas the secondary axis gives the percentageratio for all the methods in comparison to weak boundary condition (WBC).The Green-Gauss method follows the trend observed in the gradients wherethe error coefficients for the gradients were higher than the least-squares method.In this plot, the Green-Gauss method has errors almost twice the error compared574.1. Analytic Test AnalysisFigure 4.3: Comparison of boundary treatment methodsto weak boundary condition implementation. Using ghost cells or increasing thestencil size hasn’t performed better than the weak boundary condition either. Evenwith an increased stencil, the errors in 4-Unknown model and SBC were higherthan the normal stencil. Although the gradients were of comparable accuracy forthe 4-Unknown and SBC methods, the use of the second derivative to get an accu-rate gradient in y at the Gauss points doesn’t improve the truncation error much.Hence, the SBC method performs better than its counterpart. Using a single andtwo jump term for boundary control volumes result in at least 75% reduction intruncation error. Two jump term method fares better than the single jump termmethod.In view of the results for different boundary treatment methods, the use of ghostcells and the 4-Unknown model is dropped for numerical tests as the results showlittle promise compared to other schemes. The corrected Green-Gauss method ispresented only for comparison whereas the use of SBC and jump terms is furtherexplored, with and without the increase of stencil size. The numerical results fordifferent reconstruction methods, gradient averaging methods and boundary treat-ment methods are presented in the next section.584.2. Numerical Test Analysis(a) Sq-15 with 462 triangles and 264 vertices (b) Sq-120 with 28670 triangles and 14592 ver-ticesFigure 4.4: Isotropic triangular mesh cases4.2 Numerical Test AnalysisA number of different methods were tested in the previous section where someperformed as expected and some underperformed. Schemes successful in reducingtruncation error are now considered for the numerical tests where the truncationerror is calculated for a sequence of unstructured meshes. Discretization error isalso calculated by using the method of manufactured solutions, taking advantage ofexact solutions satisfying Poisson’s equation to compare. Different manufacturedsolutions are tested to ensure that the results are applicable for a broader range ofproblems.4.2.1 Isotropic Boundary Layer MeshA sequence of unstructured meshes were generated for a unit square domain, thecoarsest and finest mesh of which are presented in the Figure 4.4. The mesh sizeincreases such that the characteristic length halves in the sequence and the meshis named as Sq-# where the number is approximately equal to the square root ofthe number of vertices. The triangular elements are of high quality and the sizechanges gradually across the domain. The mesh is isotropic, both in the interiorand near the boundary, with an aspect ratio of about one near the boundary.We shall discuss the convergence of the schemes to assess the order of accuracy594.2. Numerical Test Analysisfor truncation and discretization error followed by confirmation of results observedin analytic tests with an emphasis on the flux integral and solution accuracy. Tocompare the truncation error, we use the error metric calculated at the control vol-ume reference point, as the one used for analytic tests discussed in Section 4.1.2.An L2 norm of error metrics for the control volumes of the mesh is used to com-pare the effectiveness of the schemes in achieving the set goals. However, thereare two truncation error L2 norms defined for a mesh, the interior and boundary.The interior error norm considers only the interior control volumes such that noneof the vertex neighbors lie on the boundary. It by default excludes errors for theboundary control volumes and their first neighbors which is accounted for, by theboundary error norm. The discretization error calculation is straightforward wherethe control volume average of the discrete solution is compared with the exact con-trol volume average and the L2 norms for the overall mesh is used as a measure forcomparison.Figure 4.5 gives the truncation error distribution for a Sq30 mesh using least-squares method and strong boundary condition implementation. The figure vali-dates the necessity to have a separate error norm for the boundary control volumesand their first neighbors as the magnitude of the error is higher for these controlvolumes. Using the strong boundary condition, the error on the boundary Gausspoints is driven to zero by explicitly satisfying the boundary condition through re-construction. This is the reason as to why the source term for solution error, whichis the truncation error, has higher magnitude in the boundaries. The negative mag-nitude of truncation error metric near the boundary is common among differentmeshes and can be intuitively deduced through the relation given in the Equation1.1.In the previous section, the schemes which were effective in reducing trunca-tion error were WBC, SBC, single and two jump term methods where the increasedstencil size and Green-Gauss methods were used for just comparison. However,the use of WBC for boundary condition implementation alone doesn’t convergeand needs the use of a jump term for convergence.Three manufactured solutions are used for conducting the numerical tests. Thefirst two apply homogeneous Dirichlet boundary condition whereas the third ap-plies a non-homogenous Dirichlet boundary condition. The manufactured solu-tions used are shown in the Table 4.5 and presented in the Figure 4.6. The sourceterm for each of the exact solution is manufactured by applying the Laplacian op-erator over the exact solution, and is fed to the numerical solver ANSLib.The exact solution U =−sin(pix)sin(piy) is used to conduct preliminary numer-ical tests and accumulate initial data for rigorous testing with other manufacturedsolutions. The convergence plots for discretization error is presented in the Figure4.7 to ascertain that all the schemes have consistent order of accuracy. The least-604.2. Numerical Test Analysis(a) Error in second derivative of x (b) Error in second derivative of yFigure 4.5: Truncation error distribution for Sq30 mesh using least-squares methodCase Boundary Condition Exact SolutionU1 Dirichlet Homogenous −sin(pix)sin(piy)2 Dirichlet Homogenous 20(x−1)(y−1)log(x+1)log(y+1)3 Dirichlet Non-Homogenous sin(pi(x+ y))Table 4.5: Manufactured solutions used for numerical testssquares method results are presented for both linear interpolation and arithmeticaveraging methods of face gradient calculation. Both the methods are observed toconverge with second order accuracy and the difference in discretization error isalso minor. A comparison with the Green-Gauss method with and without correc-tion reveal that the latter exhibits a first order convergence when the difference indiscrete and exact control volume averages are used as a measure of error. Hencefor our calculations, we use only the corrected version of Green-Gauss method(GGC) for comparison and the error is always calculated by comparing discreteand exact control volume averages of the solution. The corrected Green-Gaussmethod appears to fare slightly better than the least-squares method with SBCwhen it comes to having accurate solution. Its important to note that the correctedGreen-Gauss method applies the boundary flux in a way as to enforce the bound-ary condition value at the vertices as the control volume average. The increase instencil neighbor size also leads to lack of convergence for least-squares method. Inaccordance with the convergence results obtained, the two jump term method withWBC implementation, the single jump term method, the strong boundary condi-614.2. Numerical Test Analysis(a) Case 1: Sinusoidal solution with homogeneous boundary(b) Case 2: Logarithmic solution with homogeneous boundary(c) Case 3: Sinusoidal solution with non homogeneous boundaryFigure 4.6: Exact solutions used for discretization error analysis624.2. Numerical Test AnalysisFigure 4.7: Convergence plots for discretization errorstion, and the corrected Green-Gauss method are the only methods used for testingnumerically. Experiments with increased stencil are few and shown to confirm theobservations found in analytical results. Unless specified, the face gradients areaveraged by linear-interpolation.Figure 4.8 presents the results for the use of a single jump term with a strongboundary condition applied at the boundary. The error metric corresponding toα = 0 is equivalent to using only the SBC implementation for boundary. Observ-ing both the figures, we can deduce that the use of jump terms effectively reducesthe truncation error for interior and boundary control volumes. The truncation er-ror with jump terms is lower than the case with SBC, confirming the better perfor-mance of jump terms over SBC. The jump coefficient corresponding to minimuminterior truncation error is within 10% of the value α = 1.2 for the sinusoidal ex-act solution under consideration. The optimal jump coefficient to minimize theboundary truncation error is within the 10% of the value α = 1.75. This differ-ence in optimal jump coefficients in interior and boundary regions along with theintuition of using multiple variables to reduce truncation error, as evident from theEquation 4.12, is the motivation for experimenting with two dedicated jump coef-ficients, each for the boundary and interior edges of the control volumes. Since thetruncation error metrics are by definition independent of the solution used, it wasalso verified that they are not affected by changing the exact solutions tested.Referring the Equation 4.13, in the absence of a boundary jump term, the error634.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 4.8: Plot for truncation error using single jump term method644.2. Numerical Test AnalysisFigure 4.9: Case 1: Solution error plot using single jump term methodmetric for each control volume is parabolic with respect to the jump coefficientα . This behavior is observed globally for the mesh as evident from the parabolicnature of the plots for boundary and interior truncation error. Using the single jumpterm method, the discretization error for the sinusoidal exact solution is presentedin the Figure 4.9. The spike near the jump coefficient α = 0 is related to thesteep drop in error when a jump term is used. It can also be observed that thesolution error minimum has a consistent behavior and that the minimum occurs inthe vicinity of the optimal jump coefficient α = 1.2, consistent with the truncationerror results.Since truncation error metrics are independent of the problem solved for, thetrend for reducing truncation errors with the variation in jump term is the same forexact solution cases 2 and 3. Figures 4.10 and 4.11 show the results for discretiza-tion error with a single jump term for cases 2 and 3 respectively. The spike nearα = 0 shows that SBC has higher errors and the reduction obtained by using singlejump term is about two orders of magnitude. The jump coefficient correspondingto the minimum for both the cases are within 20% of the optimal jump coefficientα = 1.2 for case 1. Nevertheless, the use of single jump term method is effectivefor the test cases with regards to truncation and discretization error. It also givesconsistent results across the mesh sequence and exact solution cases, as expectedfrom analytic tests.654.2. Numerical Test AnalysisFigure 4.10: Case 2: Solution error plot using single jump term methodFigure 4.11: Case 3: Solution error plot using single jump term method664.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 4.12: Comparison of truncation error for different schemes674.2. Numerical Test AnalysisFigure 4.13: Comparison of schemes for discretization errorThe two jump term implementation is now tested and the results for truncationerror is presented in comparison to all the other schemes in the Figure 4.12. Trun-cation errors for WBC seem to be smaller than SBC but the former case doesn’tyield a converged solution. The Green-Gauss method is shown for the purposeof comparison and as such, it doesn’t give the least truncation error for the bound-ary. When compared with the use of a single jump term method, the two jump termimplementation has half the error magnitude. Thus the two jump term implementa-tion has the least boundary truncation error followed by single jump term method,GGC and then SBC. Interior truncation error is not influenced by the boundaryjump term and has the lowest truncation error using the single jump term methodfollowed by SBC and GGC methods. The order of convergence for truncation erroris observed to be zeroth.Figure 4.13 gives an insight as to how the schemes compare when solution erroris concerned. All the schemes have a consistent second order rate of convergenceas evident from the plot. Using least-squares method with SBC alone has the worstdiscretization error when compared with the corrected Green-Gauss method andthe use of jump terms. The use of single and two jump term methods does not showa greater disparity in solution error since the effect of the boundary jump term islocal and limited to the boundary control volumes. Hence we can say with certainsurety that the use of two jump terms can potentially reduce the discretization andtruncation error. The important factor though is whether all the three measures684.2. Numerical Test AnalysisMeshInterior TE Boundary TE Discretization ErrorError α Error α αB Error α αBSq15 0.0554 1.09 0.225 1.54 1.28 1.27e-03 1.21 2.24Sq30 0.0780 1.18 0.192 1.6 1.34 2.90e-04 1.24 1.25Sq60 0.0772 1.18 0.200 1.57 1.34 7.49e-05 1.21 1.28Sq120 0.0746 1.18 0.199 1.57 1.34 1.81e-05 1.18 2.96Table 4.6: Optimal jump coefficients for minimum truncation and discretizationerrorof error have a minimum at same jump coefficients or not. Table 4.6 gives thecorresponding jump coefficients associated with the error minima for interior andboundary truncation error and discretization error. We can notice that the minimafor interior truncation error and discretization error occur in the close vicinity ofα = 1.2 but the boundary truncation error minimum occurs in the neighborhood of{α = 1.57,αB = 1.34}. The discretization error corresponding to α = 1.57 is fourtimes the minimum error achievable and in most cases holds the priority.Figures 4.16 and 4.17 explain the difference in error distribution for a Sq30mesh corresponding to cases with minimum for boundary truncation error and dis-cretization error and encourage us to find a trade-off for reducing both the errormetrics simultaneously. The variation in boundary truncation error is gradual withthe change in jump coefficients. On the other hand, the discretization error changesrapidly with change in jump coefficients and in general is four times the minimumdiscretization error possible when we use jump coefficients corresponding to min-imum boundary truncation error case. The boundary truncation error is also foundto be at least two times higher than the minimum possible error when discretizationerror is minimized.To select the jump coefficients where both the boundary and discretization er-ror are minimized, we use a Pareto plot for each of the mesh and exact solutioncombination. Each point on a Pareto plot corresponds to the normalized values ofboundary truncation and discretization error for the discrete jump coefficient set onthe two jump term space. The boundary truncation error and discretization error arenormalized by their corresponding minimum error obtained using the single jumpterm method. The Pareto plot for a Sq30 mesh for case 1 exact solution is shownin the Figure 4.14, where the points closer to the origin are considered suitablefor minimizing both discretization and boundary truncation errors simultaneously.The point with jump coefficients α = 1.24 and αB = 1.11 seems to have discretiza-tion error and interior truncation error within 1% of their minimum possible valueswhereas boundary truncation error was within 8% of the minimum, thus achieving694.2. Numerical Test Analysissimultaneous reduction of both the types of errors. The error distribution corre-sponding to these jump coefficients is also presented in the Figures 4.16 and 4.17as Optimal error cases.The interior truncation error hardly varies with boundary jump coefficient asclearly visible in the Figure 4.18a. Figure 4.18b reveals the behavior of boundarytruncation error with the jump coefficients and we can see that there is a region ofgradual change near the bottom of the surface which spans a significant area in theα and αB space. This implies that boundary truncation error can be lowered overa wider range of jump coefficients. The variation of discretization error in Figure4.19 also has a flatter surface near the line α = 1.2 and has a range of valuesfor jump coefficient αB where it stays in the vicinity of the global minimum fordiscretization error. This behavior for all the three measures of error is found to bethe similar across the mesh sequence.The error distribution for case-2 and case-3 are presented in the Figures 4.20and 4.21 respectively. We can identify regions where the variation in discretizationerror is gradual and it is the overlap of this region with the region for minimumboundary truncation error that we seek. We find that such regions of minimum er-rors match for different meshes and different cases of exact solution. This result ispresented in the Figure 4.15 where the deviation of boundary truncation error anddiscretization error from their minimum possible error in the two-jump coefficientspace is given for the jump coefficient set {α = 1.24 , αB = 1.11}. Irrespectiveof where the minimum for either of the errors occur, the deviation from their min-imum is less than 30% for discretization error and 10% for boundary truncationerror. The interior truncation error deviates by less than 10% for the jump coeffi-cient set chosen.704.2. Numerical Test AnalysisFigure 4.14: Pareto plot for discretization and truncation errors corresponding tocase-1 exact solution, Sq30 meshFigure 4.15: Pareto case error deviation from minimum possible errors714.2.NumericalTestAnalysis(a) Minimum boundary TE (b) Minimum discretization error (c) Optimal error caseFigure 4.16: Comparison of truncation errors for different minimization cases, Sq30 mesh(a) Minimum boundary TE (b) Minimum discretization error (c) Optimal error caseFigure 4.17: Comparison of discretization errors for different minimization cases, Sq30 mesh724.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 4.18: Truncation error contour plot variation with two jump terms, Sq30mesh734.2. Numerical Test AnalysisFigure 4.19: Case 1: Discretization error contour plot using two jump term method,Sq30 meshFigure 4.20: Case 2: Discretization error contour plot using two jump term method,Sq60 mesh744.2. Numerical Test AnalysisFigure 4.21: Case 3: Discretization error contour plot using two jump term method,Sq120 mesh4.2.2 Stretched Boundary Layer MeshAs mentioned in the Chapter 1, the use of stretched grids in the boundary regionsis common for viscous flows. The viscous regions require high aspect ratio gridcells and the ways to achieving this condition is by either generating 90 degree tri-angles or rectangles near the boundary [7, 38, 42]. The behavior of discretizationschemes and boundary treatment methods on stretched grids is crucial to under-stand, for them to be applicable to wider scenarios of flow problems. To generatestretched triangular boundary layered mesh, first a quad layered boundary mesh isconstructed around the isotropic triangular square mesh sequence described in theprevious section. These quad layers form the new boundary enclosing the trian-gular mesh as they grow from an aspect ratio of two at the boundary to an aspectratio of one in the interior. This transition takes place at the growth rate of 1.2 andfive quad layers are required ( see Figure 5.22 in Section 5.2.2). To generate ourstretched boundary layer with triangular elements, these quad elements are ran-domly split and they acquire the same aspect ratio and growth rates as that of thestretched quad layered mesh. Note that the stretched quad layered mesh is used fornumerical tests in Section 5.2.2 and complements the results in this section. Figure4.22 shows the intermediate meshes in the sequence with a growth layer visiblein the boundary. The mesh nomenclature for stretched meshes is Sq#-Str where754.2. Numerical Test Analysis(a) Sq30-Str mesh with 1730 vertices (b) Sq60-Str mesh with 5117 verticesFigure 4.22: Stretched boundary layered triangular meshthe # corresponds to the original triangular mesh which is modified to include thestretched boundary layer.The truncation error variation using least-squares reconstruction with strongboundary condition is presented in the Figure 4.23. There is a remarkable similarityin the distribution of errors when compared to the isotropic triangular mesh; how-ever, the magnitude of errors are twice as large as they were for isotropic meshes.Truncation error is maximum in the boundary regions, as expected.The use of single jump term method has shown considerable deviation with re-gards to truncation errors. Figure 4.24 gives the variation of interior and boundarytruncation error when single jump coefficient is used. The optimal jump coefficientfor interior truncation error is α = 1.0 whereas that for the boundary truncationerror is α = 2.0. These jump coefficients are very distinct from each other and dif-ferent from the jump coefficients observed for the isotropic triangular mesh by atleast 20%. Using single jump term method, the discretization errors behave similarto the isotropic mesh cases such that the deviation in jump coefficients is less than10% compared to the isotropic mesh cases. It is worthwhile to note that the interiortruncation errors are higher for stretched meshes compared to isotropic meshes.The comparison of different schemes for discretization and boundary trunca-tion errors are presented in the Figure 4.25. The results for the two jump termscheme reveal that the truncation errors are lowest for this scheme. They also fol-low the trend observed for isotropic meshes and are five times smaller than the casewhere simple SBC scheme is used for boundary treatment. Although the Green-764.2. Numerical Test Analysis(a) Error in second derivative of x (b) Error in second derivative of yFigure 4.23: Truncation error distribution for Sq30 stretched mesh using least-squares methodGauss method has lower error compared to a single jump term method, the imple-mentation of boundary condition enforces a flux such that control volume averagematches the boundary solution. Hence, the magnitude of boundary truncation er-ror is because of the flux enforcement rather than the Green-Gauss reconstructionmethod.Figure 4.25b shows the comparison and convergence of different schemes fordiscretization error. The convergence is second order accurate and the results aresimilar to the isotropic mesh cases with the exception of the least-squares methodwith SBC when linear interpolation is used for face gradient calculation. The meshSq120-Str does not yield a converged solution when linear interpolation is usedfor face gradient calculation. However, the use of arithmetic averaging yields aconverged solution. Using jump coefficients also yielded converged solution withthe least-squares method when linear interpolation is used and the instability canbe attributed to the nature of boundary layers. The discretization errors are lowestfor single and two jump term schemes.Table 4.7 presents the optimal jump coefficients for interior and boundary trun-cation error as well as for discretization errors. The jump coefficients for boundarytruncation errors are vastly different from the results for isotropic meshes. How-ever, they show a consistent pattern across the mesh sequence. Applying the sameprocedure as for isotropic meshes, a Pareto plot is used to identify points in thetwo jump coefficient space where the truncation and discretization errors can besimultaneously reduced. The jump term coefficient set of {α = 1.24,αB = 0.84}774.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 4.24: Truncation error variation using single jump term for triangularmeshes with stretched boundary layers784.2. Numerical Test Analysis(a) Boundary truncation error(b) Discretization errorFigure 4.25: Comparison of different schemes with regards to truncation and dis-cretization error794.2. Numerical Test AnalysisMeshInterior TE Boundary TE Discretization ErrorError α Error α αB Error α αBSq15_Str 0.0968 1 0.4413 0.61 0.57 7.0E-04 1.21 3Sq30_Str 0.1023 1 0.4465 0.64 0.6 2.4E-04 1.21 0.75Sq60_Str 0.0993 1 0.4648 0.55 0.54 6.8E-05 1.18 1.65Sq120_Str 0.0941 1.03 0.4290 0.73 0.66 1.8E-05 1.18 0.42Table 4.7: Optimal jump coefficients for minimum truncation and discretizationerror, using case 1 exact solutionFigure 4.26: Pareto case error deviation from minimum possible errorsis observed to reduce truncation and discretization error consistently, across differ-ent meshes and different exact solution cases. The deviation of discretization andboundary truncation errors for the proposed jump coefficient set from the minimumpossible errors for corresponding meshes and exact solution cases is presented inthe Figure 4.26. The figure reveals that deviation of errors is less than 10% forboundary truncation error and less than 65% for discretization errors when the pro-posed jump coefficients are used. The interior truncation errors are well within10% of their minimum values. Given that the boundary truncation error drops byfive times compared to the use of SBC alone and that the discretization error dropsby 90% compared to SBC scheme, these deviations seem like reasonable trade-offs.80Chapter 5Mixed Mesh ResultsA triangular mesh near a boundary is sometimes inefficient to capture the bound-ary information for cases where the mesh is required to be aligned, for example inthe direction of flow near a wall. A quad boundary layer can however be alignedwith the flow quantities and is expected to have better accuracy compared to a tri-angular mesh. A purely triangular mesh can be modified to have a quad boundarylayer which is an effective and time saving alternative to having a pure quad meshfor complicated geometries, hence generating a mixed mesh. In this chapter, weexplore the use of mixed mesh with quad cells at the boundary and triangular cellsin the interior with the objective to reduce boundary truncation error while simul-taneously achieving accurate solution for a simple diffusion problem. We limitourselves to studying the effect of quad boundary layers for a diffusion problemleaving the effects on advection for future study.We proceed as before by using a linear algebra system to analytically evaluatethe truncation error in the flux integral for boundary, interior and also the inter-face control volumes. We employ the schemes used for a purely triangular meshon a mixed cell mesh and the ones with better reduction in truncation errors aretaken forward to implement in the numerical solver ANSLib. As with the triangu-lar mesh results, schemes performing better for the triangular mesh are expected tohave similar performance for mixed cells as well. Figure 3.3 shows the mesh usedfor analytic tests. The square mesh sequence used for triangular tests are modi-fied to have two and five quad boundary layers for numerical tests. The two quadboundary layer mesh is used to study the effect of boundary condition on the inter-face layer whereas for a five quad layered boundary mesh, the interface layer is inthe interior. The following sections present the results for analytic and numericaltests.5.1 Analytic Test AnalysisThis section presents a comprehensive analysis of the schemes for mixed cells.Control volume i = 25 represents an interior control volume in Figure 3.3 whereasi = 40 is a boundary control volume. The interface control volume is indexed as815.1. Analytic Test Analysisi = 10. The uniform mixed mesh is perturbed in the same way as the triangularmesh was perturbed, as described in Section 4.1 and summarized in the Figure 4.1for the five different distortions performed on the mesh. The following sectionsscrutinize gradient accuracy and flux integral accuracy followed by the effect ofboundary treatment methods5.1.1 Reconstruction AccuracyThe reconstruction accuracy for quad cells in the interior, boundary and interfacecontrol volumes is expected to be first order for both the corrected Green-Gaussmethod and the least-squares method. The gradients calculated from either meth-ods take the form of a linear combination of control volume averages in the stencilof i. Following the exact same procedure as for the triangular mesh, the control vol-ume averages in the gradient calculations are expressed as a Taylor series expansionabout the reference control volume to compare the accuracies. The gradient formfor a typical uniform boundary and interface control volume, using least-squaresmethod are given here.∇˜φ 40 =∂˜ φ∂x∂˜ φ∂y40=12h(φ 41−φ 39)43h(φ 35−φ 40) (5.1)=∂φ∂x∣∣∣40+ 14∂ 2φ∂x∂y∣∣∣40h+O(h2)∂φ∂y∣∣∣40+ 23∂ 2φ∂y2∣∣∣40h+O(h2)∇˜φ 10 =∂˜ φ∂x∂˜ φ∂y10=1h(−0.4φ 9−0.2φ 5 +0.4φ 11 +0.2φ 6)1h(0.352(φ 5 +φ 6)−0.323φ 10−0.381φ 15) (5.2)=∂φ∂x∣∣∣10+0.149 ∂2φ∂x∂y∣∣∣10h+O(h2)∂φ∂y∣∣∣10+(0.084 ∂2φ∂x2∣∣∣10+0.069 ∂2φ∂y2∣∣∣10)h+O(h2)Equation 5.1 gives the discrete gradient for a boundary control volume calcu-lated using least-squares for the perfect mesh whereas Equation 5.2 gives the gra-dient for an interface control volume. The discrete gradient is then subtracted from825.1. Analytic Test Analysisthe exact gradient for an arbitrary control volume, which is given by the Equation4.3. The general form of the error in gradient is the same as in Section 4.1.1 andpresented here for first order gradients.∇˜φ i−∇φ i =(∂ 2φ∂x2∣∣∣iξxx + ∂2φ∂x∂y∣∣∣iξxy + ∂2φ∂y2∣∣∣iξyy)h+O(h2)(∂ 2φ∂x2∣∣∣iηxx + ∂2φ∂x∂y∣∣∣iηxy + ∂2φ∂y2∣∣∣iηyy)h+O(h2) (5.3)The Tables 5.1, 5.2 and 5.3 give the gradient error coefficients for boundary,interior and interface control volumes for least-squares and corrected Green-Gaussmethod. The Green-Gauss method without correction has the same behavior onquad cells as the triangular cells, giving a zero-order gradient for all the three typesof control volumes presented. Except for uniform interior quad cells, the orderof gradient is zero for uniform or random mesh for all types of control volumes.The gradient errors for least-squares method are again higher in the normal direc-tion compared to the tangential direction. The error for corrected Green-Gaussmethod is higher than least-squares method in general, a confirmation of the re-sults observed for the triangular test mesh. The error coefficients for the gradientsof interface control volumes are comparable to interior control volumes. Hence wecan expect the truncation error for interface control volumes to be similar to thatof an interior one. It also proves that the least-squares and Green-Gauss methodare effective in calculating the gradients at the interface layer without producingsignificant errors.5.1.2 Flux Integral AccuracyTable 5.4 gives the order of error in the flux integral which is calculated usingthe method described in Section 3.1 for boundary, interior and interface controlvolumes. The results for the least-squares method (LSM) are identical for botharithmetic averaging and linear interpolation methods. The order of truncation forcorrected Green-Gauss method is consistent with the expected zero order accuracyfor general cases. Without correction, the Green-Gauss method gives one orderless accurate flux integral for random and curved meshes, consistent with the ob-servations for a triangular mesh. Since the order of flux integrals for a randommesh is zero, the error metric used for comparing the boundary treatment methodsis the same as the one used for triangular mesh.835.1. Analytic Test AnalysisMesh ξxx ξxy ξyyLeast-SquaresUniform 0 0 0Random −.12120 −0.01514 −0.00620Green-Gauss CorrectedUniform −0.5000 0.37500 0Random −.46759 .39745 −0.00631(a) Error coefficients of ∂φ∂xMesh ηxx ηxy ηyyLeast-SquaresUniform 0 0 0.41667Random −0.04690 0.01531 0.43451Green-Gauss CorrectedUniform 0 −0.50 0.41667Random 0.00732 −0.47419 0.43092(b) Error coefficients of ∂φ∂yTable 5.1: Gradient error coefficients for boundary control volumesMesh ξxx ξxy ξyyLeast-SquaresUniform 0 0 0Random 0.21196 0.00479 −0.00307Green-Gauss CorrectedUniform −0.5 0.5 0Random −0.50646 0.50200 0.00234(a) Error coefficients of ∂φ∂xMesh ηxx ηxy ηyyLeast-SquaresUniform 0 0 0Random −0.00128 −0.02412 −0.30029Green-Gauss CorrectedUniform 0 −0.5 0.5Random 0.00128 −0.50026 0.50071(b) Error coefficients of ∂φ∂yTable 5.2: Gradient error coefficients for interior control volumes845.1. Analytic Test AnalysisMesh ξxx ξxy ξyyLeast-SquaresUniform 0 0.17916 0Random 0.01349 0.11940 −0.00992Green-Gauss CorrectedUniform −0.5 0 0Random −0.50446 0.03520 0.02191(a) Error coefficients of ∂φ∂xMesh ηxx ηxy ηyyLeast-SquaresUniform 0.08411 0 0.09926Random 0.03802 0.04928 −0.19762Green-Gauss CorrectedUniform −0.14369 −0.5 0.44424Random −0.14255 −0.42280 0.46126(b) Error coefficients of ∂φ∂yTable 5.3: Gradient error coefficients for interface control volumesThe order of truncation errors are in agreement with the results of Diskin andThomas [9]. It should also be noted that a negative order of convergence for trun-cation error doesn’t necessarily mean a lack of convergence in solution. However,the quad elements have a truncation error which is decoupled and results in solu-tion divergence. Using jump terms has been observed to alleviate this problem,which means using the reconstruction methods alone will not be able to achieve aconverged solution. We also implement the gradient correction described by Hasel-bacher and Blazek [15] and Diskin and Thomas [9] to obtain strongly coupled sten-cils for both least-squares and Green-Gauss method. This method is presented inthe Section 2.2.1 for the Green-Gauss method and the same correction is appliedfor the least-squares method as well.Using a jump coefficient of α = 1 gives the same exact results for flux integralas it gives with gradient correction for least-squares method on a uniform mesh.However, this behavior deviates for random meshes. The gradient correction re-places the gradient along the edge direction of the two control volumes with anedge derivative whereas a jump term introduces a damping effect to reduce er-rors. For instance, referring to Section 2.2.1, and defining the directions eˆ alongthe edge of [0,2] and nˆ normal to the edge in Figure 2.2 and replacing the gradi-ent on the edge AB with face gradient ∇φ02, the gradient correction used for both855.1. Analytic Test AnalysisSchemeBoundary MeshUniform Sheared Scaled Stretched Curved RandomLSM O(1) O(1) O(1) O(1) O(1) O(1)GG O(1) O(1) O(1) O(1) O(1/h) O(1/h)GGC O(1) O(1) O(1) O(1) O(1) O(1)(a) Boundary control volumesSchemeInterior MeshUniform Sheared Scaled Stretched Curved RandomLSM O(h2) O(h2) O(h2) O(1) O(1) O(1)GGC O(h2) O(h2) O(h2) O(1) O(1) O(1)(b) Interior control volumesSchemeInterface LayerUniform Sheared Scaled Stretched Curved RandomLSM O(1) O(1) O(1) O(1) O(1) O(1)GGC O(1) O(1) O(1) O(1) O(1) O(1)(c) Interface control volumesTable 5.4: Truncation error order comparison for flux integralleast-squares and Green-Gauss methods to evaluate flux on edge AB is given as:∇φAB = ∇φ02 +[φ 2−φ 0|r2− r0|− (∇φ02 · eˆ)]eˆ (5.4)Without the gradient correction, ∇φ02 would be the gradient on the edge ABcalculated either through gradient averaging scheme for least-squares or the gra-dient for the primal cell using the Green-Gauss method. However, to obtain acompact stencil, we can use the correction for the gradient for every edge. To il-lustrate the effects better, we decompose the gradient ∇φ02 into the edge and edgenormal direction as shown in Equation 5.5 and substitute in the above equation tosee the effect of the gradient correction.∇φ02 = (∇φ02 · eˆ) eˆ+(∇φ02 · nˆ) nˆ (5.5)∇φAB = (∇φ02 · nˆ) nˆ+[φ 2−φ 0|r2− r0|]eˆ (5.6)Thus the gradient along the edge direction is replaced by an edge-derivativeand the normal component of the gradient is unchanged, which is shown to result865.1. Analytic Test AnalysisFigure 5.1: Gradient error metrics for different schemesin a strong coupling. The flux integral without gradient correction for an interiorcontrol volume i = 25 in Figure 3.3 is comprised of second neighbors, as shown inthe equation below:∇˜2φ i =14h2(φ 23 +φ 27 +φ 35 +φ 15−4φ 25) (5.7)whereas the flux integral with gradient correction is comprised of immediateneighbors, thus resulting in a compact stencil as shown:∇˜2φ i−GC =1h2(φ 24 +φ 26 +φ 30 +φ 20−4φ 25) (5.8)A similar effect is also seen for triangular mesh, though the solution is stableeven without the gradient correction.5.1.3 Boundary Treatment MethodsThe schemes tested for triangular mesh analysis are applied for quad meshes toverify their consistency. The schemes intended to improve gradient accuracy arethe SBC, 4-Unknown and increased stencil size methods. These schemes are testedand compared in reference to the WBC method. The effect of stencil change isshown for all three methods in Figure 5.1 where using increased stencil has highererror metrics compared to a stencil choice of first neighbors. Also, the error in thegradients is minimum for 4-Unknown scheme for the normal gradient followed bySBC and then WBC.875.1. Analytic Test AnalysisJump terms can also achieve improved truncation error as observed for theuniform triangular mesh. For a uniform quadrilateral mesh, using jump terms yielda truncation error of the form given in the Equation 5.9. However, there is no gainin accuracy for an interface control volume by using a single jump term methodand the order of error remains zero. This is not the case at the boundary, however,where the quadrilateral control volumes allow for a cancellation in error in thetangential direction, leading to a truncation error of the form given in Equation5.10.∇˜2φ int −∇2φ int =h24(43−α)(∂ 4φ∂x4∣∣∣∣+∂ 4φ∂y4∣∣∣∣)+O(h4) (5.9)∇˜2φ bdry−∇2φ bdry =7(1−α)15∂ 2φ∂y2∣∣∣∣+O(h) (5.10)The optimal jump coefficient for a uniform interior quadrilateral control vol-ume is α = 43 , the same value observed in uniform triangular mesh results and byNishikawa [30], to yield higher order truncation error. The optimal coefficient forthe boundary control volume, however, is α = 1. The search for an optimal jumpcoefficient for random meshes using single or two jump terms to gain higher or-der truncation error was not fruitful. However, the observed jump coefficients arecontrasted with the optimal jump coefficients to observe for any similarity.The different boundary treatment methods are now compared for their perfor-mance in reducing truncation error and Figure 5.2 presents the results for a fixedrandom perturbation of the uniform mesh with weak boundary condition as a ref-erence scheme. SBC, single jump term and two jump term method perform muchbetter than increasing the stencil size or 4-Unknown model. Using gradient correc-tion with SBC have truncation errors comparable to the case where just the SBCscheme is used. Hence we continue with the schemes which are performing wellfor the boundary truncation error which include single and two jump term methods.We can also explore the influence of using single jump term method and twojump term method when the quad layers are stretched and randomized simultane-ously. This analysis is performed in two ways: one by stretching a random meshand the other by increasing the randomness of a stretched mesh. The random mesh,generated by perturbing coordinates of a uniform mesh by a maximum of 10% ofthe edge length, is stretched increasingly from stretch ratios of 0.75 to 1.5. The ran-domness of a stretched mesh with a stretch ratio of 1.2, is increased by perturbingcoordinates proportional to a certain percentage of their edge length. The results885.2. Numerical Test AnalysisFigure 5.2: Comparison of boundary treatment methodspresented in the Figure 5.3 reveal that the two jump term scheme is consistentlybetter at reducing boundary truncation error compared to other methods.The truncation error for interior and interface control volumes are also com-pared in Figure 5.4 to see if the effect of having mixed cells produces abrupt in-crease in error in the interior domain. The errors for interface control volumescompared to interior control volumes are higher for the least-squares method andthe single jump term method whereas they are lower for the corrected Green-Gaussmethod. But when contrasted with the boundary truncation error in Figure 5.2 forthe reference scheme of WBC, the errors for both the control volume types aresignificantly lower. Hence it can be concluded that the reconstruction schemes areeffective in calculating the gradients efficiently for the interface layers and thusthey can be omitted from the regions of concern for lowering truncation errors.5.2 Numerical Test AnalysisThe results observed for analytic tests on a mixed mesh have been consistent withthe results observed for a triangular mesh, especially when it comes to provingthat using a second jump coefficient gives an extra control variable to reduce theerrors, both in truncation and discretization. It has been also pointed out that fora mixed mesh, we encounter the problem of decoupling which results in lack ofsolution convergence. We use the exact solutions given in Table 4.5 which satisfythe Poisson’s equation with a suitable source term for numerical tests. The useof exact solutions simplifies discretization error calculation; the truncation error895.2. Numerical Test Analysis(a) Increasing stretch ratio for random mesh(b) Increasing randomness of a stretched meshFigure 5.3: Comparison of single and two jump term methods on randomlystretched meshesFigure 5.4: Comparison of interior and interface truncation errors905.2. Numerical Test Analysismetrics are calculated as explained in Section 3.2.1 and are independent of theexact solution used but depend on the mesh.5.2.1 Isotropic Boundary Layer MeshWe continue to use a unit square domain for our numerical tests; however the tri-angular mesh is modified to include quad boundary layers. Each of the triangularsquare meshes used for numerical tests in the previous chapter now has two and fiveisotropic quad layers added at the boundary, resulting in two different sequences ofmeshes. The mesh is scaled such that the domain is a unit square and has the sametriangular mesh in the interior with the quad layers grown near the boundary. Themesh is of high quality with gradual change in size and is presented in the Figure5.5. The square mesh with five isotropic quad layers is used primarily to analyzethe schemes to reduce truncation error whereas the square meshes with two quadlayers is used to observe the impact of having an interface layer interact with theboundary control volumes. The interface layer in two quad layered meshes haveflux integrals which include the boundary control volume information directly. Fi-nally, we also study a sequence of pure quad meshes which have characteristiclength scales comparable with the isotropic pure triangular mesh. They are labeledas Quad-#x# where the number is equal to the square root of the total number ofprimal cells in the mesh. The nomenclature for the mixed mesh is straightforward.Truncation error distributions for square meshes with two and five quad bound-ary layers is presented in the Figures 5.6 and 5.7 respectively. The error distri-butions observed for second derivatives in the x and y direction reveal that theboundary control volumes and their first neighbors are the regions of maximumerror. Hence we continue with our approach of using boundary error metric toaccount for the truncation error in boundary control volumes and their first neigh-bors whereas an interior truncation error metric measures the error for the interiorcontrol volumes. Some conclusions can be readily made from the figures, whichare that the interface layer blends well with the interior for a square mesh withfive quad layers. The fact that error distributions for a square mesh with two quadlayers are similar to that of a five quad layered mesh with no spikes in the interfacelayer confirms that a transition from triangular to quad mesh is smooth.915.2. Numerical Test Analysis(a) Sq30 mesh with five isotropic quad layers (b) Sq60 mesh with five isotropic quad layers(c) Sq30 mesh with two isotropic quad layers (d) A 32x32 pure quad meshFigure 5.5: Isotropic meshes for numerical analysis925.2. Numerical Test Analysis(a) Error in second derivative of x (b) Error in second derivative of yFigure 5.6: Truncation error distribution for Sq30 mesh with two isotropic quadlayers, using LSM(a) Error in second derivative of x (b) Error in second derivative of yFigure 5.7: Truncation error distribution for Sq30 mesh with five isotropic quadlayers, using LSMThe schemes found analytically to be effective in reducing truncation error fora mixed mesh are the least-squares method with strong boundary condition (SBC)and the single and two jump term methods. The two jump term method uses a weakboundary condition whereas the single jump term method uses a strong boundarycondition; both employing least-squares gradient reconstruction. Since the useof least-squares method with SBC alone doesn’t yield a converged solution, the935.2. Numerical Test Analysismethod is used only for comparison of truncation errors. However, we focus on theuse of least-squares method with gradient correction and strong boundary condi-tion (LSM-GC) to compare with other schemes for reduction of error in truncationand discretization. The corrected Green-Gauss method is also presented as an al-ternative method of solution reconstruction to compare with least-squares methodin general.The convergence of the solution is first verified for all three mesh sequences inFigure 5.8 and for this purpose the least-squares method with gradient correction(LSM-GC) is used to avoid decoupling. We compare the face-gradient averagingmethods which include the linear interpolation (LI) and arithmetic averaging (AA)for least-squares method and find that the differences are negligible for discretiza-tion errors and the order of convergence is indeed second order with the gradientcorrection as well. The Green-Gauss method with correction (GGC) and withoutcorrection (GG) are also compared for solution convergence. As expected, the nor-mal Green-Gauss method shows a one order lower accuracy in convergence of theerrors in control volume averages; this is well remedied by the corrected Green-Gauss method for gradient calculations, showing a second order convergence. Thegradient corrected least-squares method has better accuracy when compared to theGreen-Gauss method. The errors for pure quad meshes are comparable but lowerthan the mixed meshes.Figures 5.9 and 5.10 give the variation of truncation errors using the singlejump term method for the mesh sequence with two and five quad layers along withpure quad meshes. It can be observed that the jump coefficient required to min-imize the boundary truncation error is starkly different from the jump coefficientrequired to minimize interior truncation error. The optimal jump coefficient forinterior truncation error for the mixed meshes is in the neighborhood of α = 1whereas for the pure quad meshes, it varies little from the optimal coefficient ofα = 0.9 across the sequence. This is different from the results observed for the in-terior truncation error obtained for a pure triangular mesh where the optimal jumpcoefficient was around α = 1.2.The optimal jump coefficients for boundary truncation error metrics is aroundthe value of α = 1.75 for both the mixed and pure quad meshes, and it is strikinglycloser to the value observed for a pure triangular mesh. However, the parabolicnature of the curve for truncation error metrics, both interior and boundary, withrespect to the jump coefficients confirms the results predicted by analytic tests andthose observed for triangular meshes. It also implies that the flatter minimum forthe curve can accommodate a broader range of jump coefficients to achieve reducedtruncation error.The variation of discretization error for mixed meshes with five quad layers ispresented in Figure 5.11. The optimal jump coefficients for a sequence of meshes945.2. Numerical Test Analysis(a) Square meshes with five isotropic boundary quad layers(b) Square meshes with two isotropic boundary quad layers(c) Pure Quad mesh sequenceFigure 5.8: Convergence of solution using least-squares and Green-Gauss methods955.2. Numerical Test Analysis(a) Square meshes with two quad layers(b) Square meshes with five quad layers(c) Pure quad meshesFigure 5.9: Interior truncation error using single jump term965.2. Numerical Test Analysis(a) Square meshes with two quad layers(b) Square meshes with five quad layers(c) Pure quad meshesFigure 5.10: Boundary truncation error using single jump term975.2. Numerical Test Analysisfor a particular case of exact solution is consistent; however, there are variationsas the exact solution changes. The optimal jump coefficient varies from the coef-ficients α = 0.9 for logarithmic solutions to α = 1.2 for sinusoidal solution. Theoptimal jump coefficient for non-homogeneous sinusoidal exact solution case isα = 1. These results were found to be similar for the mixed mesh sequence withtwo quad layers as well.The variation of discretization error for the pure quad mesh sequence is given inFigure 5.12. The optimal jump coefficients are again consistent across the sequencefor a particular exact solution but are dissimilar across the cases. The optimal co-efficients in comparison with the mixed mesh are very similar for cases 2 and 3 butfor case 1, they resemble the optimal jump coefficient for a uniform quadrilateralmesh. The optimal jump coefficient of α = 43 was observed to yield a fourth ordertruncation error for uniform mesh, as shown by the Equation 5.9 and is also thepoint where case 1 yields minimum discretization error.The conclusions for mixed and pure quad meshes are harder to draw as thereis a subtle variation in the jump coefficients with the change in exact solutions.However, it should be noted that there is a certain degree of consistency for theoptimality of jump coefficients when the mesh is changed in a sequence. Thereduction in error by the use of jump coefficients is certainly an effect seen for allthe cases of exact solutions and mesh sequences. The drop in error for the exactsolutions tested is at least twenty times compared to no-jump term scenarios, acrossthe mesh sequences for triangular, mixed and pure quad meshes.So far, we have discussed the results for a single jump term method. The othermethods under consideration are the use of least-squares method with gradientcorrection and two jump term method. The results for interior and boundary trun-cation error using different schemes are compared in the Figures 5.13, 5.14 and5.15 for mixed two and five quad layered, and pure quad meshes respectively. Inte-rior truncation error is maximum for the case of Green-Gauss method followed byleast-squares method without gradient correction (LSM) and with gradient correc-tion (LSM-GC). The use of jump term has considerably lower truncation error asexpected. The over all truncation error in the interior is lesser for the quad elementsfor all the schemes when compared to the mixed meshes.Near the boundary, the use of jump terms and gradient corrected least-squaresmethod reduces the truncation error drastically when contrasted with the Green-Gauss method or the normal least-squares method. Consistent with the observa-tions for triangular mesh, the use of two jump terms produces the least truncationerror near the boundary followed by the gradient corrected least-squares methodused to reduce decoupling in the flux integral. The use of gradient correction issimilar but not the same as the use of jump terms and hence the lower truncationerror when compared to other methods. The magnitudes and behavior for all the985.2. Numerical Test Analysis(a) Case 1: Sinusoidal solution, homogeneous boundary condition(b) Case 2: Logarithmic solution, homogeneous boundary condition(c) Case 3: Sinusoidal solution, non-homogeneous boundary conditionFigure 5.11: Plot for solution error on square meshes with five boundary quadlayers, using the single jump term method995.2. Numerical Test Analysis(a) Case 1: Sinusoidal solution, homogeneous boundary condition(b) Case 2: Logarithmic solution, homogeneous boundary condition(c) Case 3: Sinusoidal solution, non-homogeneous boundary conditionFigure 5.12: Plot for solution error on pure quad square meshes, using the singlejump term method1005.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 5.13: Comparison for truncation error on square meshes with two quadlayers1015.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 5.14: Comparison for truncation error on square meshes with five quadlayers1025.2. Numerical Test Analysis(a) Interior truncation error(b) Boundary truncation errorFigure 5.15: Comparison for truncation error on pure quad meshes1035.2. Numerical Test Analysisthree mesh sequences is identical since the boundary layer is basically composedof quadrilateral elements.A summary of the comparison of schemes with respect to the discretizationerror is presented in the Figure 5.16. In case of discretization error, the use of jumpterms acts as a damping term which alleviates the decoupling issues otherwise ob-served when the simple least-squares method is used. We observed in triangularmesh results that the use of the least-squares method alone had errors significantlyhigher than the Green-Gauss method and the jump term methods. For mixed andpure quad meshes, we lack such a comparison as the least-squares fails to convergeby itself. This also explains why the use of least-squares method with gradientcorrection has errors comparable to those observed with the jump term methods.The discretization error continues to be less dependent on the boundary jump co-efficient compared to interior jump coefficient and hence the errors are very closefor both single and two jump term methods.The discussion is now focused on exploring the jump coefficients for achievingsimultaneous reduction in truncation and discretization errors. Table 5.5 presentsthe results for the the optimal jump coefficients associated with each of the inte-rior truncation error, boundary truncation error and discretization error. The widedisparity in the jump coefficients required to minimize different aspects of errorscan be observed for different mesh sequences. The results for the mixed squaremesh sequence with two and five quad layers are very similar with regards to allthe three measures of errors. Only the truncation error is minimized in a narrowband of jump coefficients for mixed and pure quad meshes.The question now remains to find the jump coefficients that simultaneously re-duce truncation errors for interior and boundary as well as discretization errors. Weaim to keep the optimal scenario errors within 20% of the minimum possible errorsfor interior and boundary truncation error and the same is anticipated for discretiza-tion error. The discretization errors and boundary truncation errors correspondingto each discrete jump coefficient set in the two jump term space are normalized bytheir respective minimum errors obtained using the single jump term method andare plotted as a Pareto plot. Figure 5.17 shows the Pareto plots for case 1 exactsolution on a mixed mesh with five quad layers and a pure quad mesh. The pointsof interest are the points in the bottom-left corners which are closer to both the axisand are the points of minimum for boundary truncation errors and discretizationerrors.This analysis, similar to the isotropic triangular mesh case, is performed forall the exact solution cases and different mesh sequences. A particular point onthe jump term space is explored for the mixed mesh and pure quad mesh sequenceand the deviation from the minimum errors possible are calculated. Figure 5.18shows the results for deviation of the errors from minimum possible for the mesh1045.2. Numerical Test Analysis(a) Square meshes with two isotropic boundary quad layers(b) Square meshes with five isotropic boundary quad layers(c) Pure quad mesh sequenceFigure 5.16: Comparison of solution using least-squares and Green-Gauss methods1055.2.NumericalTestAnalysisMeshInterior TE Boundary TE Discretization ErrorError α Error α αB Error α αBSq15_ISO_2Quad 0.0893 1 0.2029 1.87 1.32 8.8E-04 1.18 3Sq30_ISO_2Quad 0.0992 1 0.2019 1.84 1.32 2.6E-04 1.18 3Sq60_ISO_2Quad 0.0967 1 0.2041 1.87 1.32 7.1E-05 1.18 2.97Sq120_ISO_2Quad 0.0931 1.03 0.2017 1.87 1.32 1.8E-05 1.18 2.94(a) Square mesh with two quad layersMeshInterior TE Boundary TE Discretization ErrorError α Error α αB Error α αBSq15_ISO_5Quad 0.0820 0.94 0.1788 1.81 1.29 5.0E-04 1.18 3Sq30_ISO_5Quad 0.0916 0.97 0.1830 1.78 1.29 1.9E-04 1.18 3Sq60_ISO_5Quad 0.0923 0.97 0.1797 1.78 1.29 6.0E-05 1.15 3Sq120_ISO_5Quad 0.0910 1 0.1822 1.78 1.29 1.7E-05 1.18 2.94(b) Square mesh with five quad layersMeshInterior TE Boundary TE Discretization ErrorError α Error α αB Error α αBQuad16x16 0.0541 0.88 0.1716 1.84 1.32 3.0E-04 1.36 1.47Quad32x32 0.0609 0.88 0.1681 1.84 1.32 7.3E-05 1.33 1.71Quad64x64 0.0594 0.88 0.1694 1.78 1.29 1.8E-05 1.33 1.74Quad128x128 0.0595 0.88 0.1642 1.84 1.32 4.7E-06 1.33 2.25(c) Pure quad meshesTable 5.5: Optimal jump coefficients for minimum truncation and discretization error, using case 1 exact solution1065.2. Numerical Test Analysis(a) Sq30 mesh with five quad layers(b) Quad 32x32 meshFigure 5.17: Pareto plot for discretization and boundary truncation error corre-sponding to case-1 exact solution1075.2. Numerical Test Analysissequences and for each of the exact solution cases. The gain in accuracy for bound-ary truncation error using two jump term method is about eight times the errorwhen simple least-squares method is used. The gain in solution accuracy is about20 times when jump terms are used compared to the case of none or almost zerojump coefficients. Thus using Pareto plot, we choose points so that boundary trun-cation errors deviate by a maximum of 20% from their global minimum while thediscretization error are kept as close to the global minimum as possible.The results for the mixed mesh sequence with five quad layers is presented forthe jump coefficients {α = 1.24,αB = 1.11} chosen from the Pareto plot whereasfor the pure quad mesh sequence we used the jump coefficients {α = 1.36,αB = 1.11}.The solution errors for mixed meshes deviate by a maximum of 35% from the min-imum possible errors for all exact solution cases. For pure quad meshes, the devi-ation from the minimum possible error is about 20% for cases 1 and 3 but for case2, the deviation is about 100%. Despite the variations, the drop in discretizationerror is at least ten times or more the error for reference case of almost zero jumpcoefficient, in all the scenarios. The deviation for boundary truncation error formixed and pure quad meshes is about 20% from the minimum possible errors. Theinterior truncation error deviates by about 10% for the mixed mesh sequence oftwo and five quad layers, for the jump coefficients used for Pareto case. The resultsfor mixed meshes with two quad layers follow the results observed for the mixedmeshes with five quad layers. The interior truncation error for the pure quad meshsequence is on the higher side of 20-25% for the jump coefficient of α = 1.36The variation of boundary truncation error using two jump term method fora five quad layered mixed mesh and a pure quad mesh is shown in the Figure5.19. The contours are similar since the boundary has quad layers for both themeshes. The contours for interior truncation error are parabolic and independentof the variation in the boundary jump coefficient. The variation of discretizationerrors are also given for the five quad layered mixed mesh and pure quad mesh inthe Figures 5.20 and 5.21 respectively. The region for lower errors are concentratedalong lines of constant alpha for different cases of exact solutions. There is slightdependence on the boundary jump coefficients for discretization errors.1085.2. Numerical Test Analysis(a) Square mesh sequence with five quad layers(b) Pure quad mesh sequenceFigure 5.18: Pareto case error deviation from minimum possible errors1095.2. Numerical Test Analysis(a) Sq30 mesh with five quad layers(b) Quad 32x32 meshFigure 5.19: Variation of boundary truncation error in two-jump term space1105.2. Numerical Test Analysis(a) Exact Solution Case 1(b) Exact Solution Case 2(c) Exact Solution Case 3Figure 5.20: Variation of discretization error in two-jump term space for Sq30mesh with five Quad layers1115.2. Numerical Test Analysis(a) Exact Solution Case 1(b) Exact Solution Case 2(c) Exact Solution Case 3Figure 5.21: Variation of discretization error in two-jump term space for a 32x32pure quad mesh1125.2. Numerical Test Analysis5.2.2 Stretched Boundary Layer Mesh(a) Sq30 mesh with two stretched quad layers (b) Sq30 mesh with five stretched quad layersFigure 5.22: Stretched quad layer mixed meshesThe meshes tested in the previous section had an isotropic quad boundary layer.Sometimes, a growth in the boundary layer is desired to allow for high aspect ratioscloser to the boundary without increasing the cell count drastically. This sectiondiscusses the results for such a sequence of meshes where the boundary layer hashigh aspect ratio cells near the wall which increase in height as they grow towardsthe interior. Two and five quad layers are grown around the pure triangular meshsequence with a growth rate of 1.2 to complement the isotropic mixed meshes thathave been previously tested. The aspect ratio at the interface layer is maintainedas one whereas for five quad layered mesh, the aspect ratio is two for the quadsat the boundary. The aspect ratio is 1.25 for the boundary quads in the two quadlayered mesh. The meshes are shown in the Figure 5.22 where the growth in thequad layers is moderate but significant enough to differ from the isotropic meshsequence. The two quad layered mesh doesn’t show a greater deviation from theisotropic mesh but the effect of stretched mesh can be better studied with the meshsequence of five stretched quad layers.The distribution of truncation error metrics for a Sq30 mesh with five stretchedquad layers is presented in Figure 5.23. The results are very similar to the case of asquare mesh with isotropic quad layers when least-squares reconstruction is used.The solution convergence was also found to be second order when the least-squaresmethod with gradient correction or the corrected Green-Gauss method was used.1135.2. Numerical Test Analysis(a) Error in second derivative of x (b) Error in second derivative of yFigure 5.23: Truncation error distribution for Sq30 mesh with five stretched quadlayers using LSMThe normal Green-Gauss method had the same behavior as before, with first orderconvergence. This trend was the same for both mesh sequences with two and fivestretched quad layers.The variation of interior truncation error with a single jump term was almostidentical to the cases where the boundary layer wasn’t stretched. However, aninspection of the boundary truncation error reveals that the optimal jump coeffi-cients for minimum truncation error differs from the isotropic quad layer cases.Figure 5.24 reveals the differences in the results for stretched quad layers. Theoptimal jump coefficient of α = 1.75 for two quad layered meshes is similar tothe case of isotropic boundary layered meshes but the optimal jump coefficient forfive stretched quad layered meshes is α = 1.98. This optimal jump coefficient forfive quad layered stretched meshes is the same as the optimal jump coefficient ofα = 2 for stretched triangular meshes which were in fact derived from the samefive quad layered stretched mixed meshes. It is also noted that the truncation errorsare marginally higher for stretched meshes.We compared the boundary truncation error for jump term schemes, gradientcorrected least-squares reconstruction and the corrected Green-Gauss method asshown in Figure 5.25. Contrary to our expectations from results for isotropic meshwhere the two jump term method had the lowest boundary truncation errors, the useof a single jump term appears to have the lowest truncation error for a stretchedmesh with five quad layers. This also contradicts the prediction in the analytictests where randomly stretched meshes (see Figure 5.3) showed two jump termscheme to be superior to single jump term method in reducing boundary truncation1145.2. Numerical Test Analysis(a) Mesh sequence with two stretched quad layers(b) Mesh sequence with five stretched quad layersFigure 5.24: Boundary truncation error variation using single jump term method1155.2. Numerical Test Analysiserror. The triangular stretched meshes were consistent with the results observedfor isotropic meshes and also had lowest errors for two jump term scheme than anyother method. This anomaly in results warrants further inspection especially forstretched quad meshes.Also presented in the Figure 5.25 are the results for the least-squares methodwith SBC which is given for the sake of comparison even though the solutiondoesn’t converge for this scheme. The Green-Gauss method shows a considerablyhigher truncation error whereas least-squares method with gradient correction hasreasonably good accuracy, just after the two jump term method. For a five quadlayered sequence, the single jump term method has the lowest error whereas fortwo quad layered meshes, the errors for the two jump term method are the lowest.The mixed stretched quad layered meshes have results similar to isotropicmixed meshes when it comes to discretization error and interior truncation errors.The interior truncation errors behave the same way as they did for isotropic quadlayered meshes whereas the results for discretization error are given in the Fig-ure 5.26. The two jump term and single jump term method are again the best inreducing the discretization error followed by the gradient correction used for least-squares method.Since the single jump term method has lower boundary truncation error com-pared to two jump term scheme, the simultaneous reduction of errors is difficultto be limited to a minimum deviation of 20%. Applying the same procedure usedfor isotropic boundary layer mesh, jump coefficients are identified in the two jumpterm space to reduce truncation and discretization errors. The deviation of errors,corresponding to jump coefficient set of {α = 1.24,αB = 0.84}, from the mini-mum possible errors for boundary truncation and discretization error are presentedin the Figure 5.27. The discretization errors deviate by a maximum of 100% fromtheir minimum possible values for all exact solution cases, whereas the boundarytruncation error deviates by about 35% from their minimum values achieved us-ing the single jump term method. An alternate approach was to use a single jumpterm scheme and the best errors were realized at the interior jump coefficient ofα = 1.4. The boundary truncation error deviated by 35% but discretization errorsdeviate by at least 100-400%. Hence, the two jump term scheme serves as thebest scheme in reducing truncation errors and discretization errors simultaneously,across different mesh types. However, prior to making any assumptions about thebest approach suited when the boundary has stretched quad layers, it warrants fur-ther inspection of results and causes for the behavior of two jump term scheme.Hence, this discussion is left for future studies.1165.2. Numerical Test Analysis(a) Mesh sequence with two stretched quad layers(b) Mesh sequence with five stretched quad layersFigure 5.25: Comparison of boundary truncation error for different schemes onmixed mesh1175.2. Numerical Test Analysis(a) Mesh sequence with two stretched quad layers(b) Mesh sequence with five stretched quad layersFigure 5.26: Comparison of solution error for different schemes on mixed mesh1185.2. Numerical Test AnalysisFigure 5.27: Pareto case error deviation from minimum possible errors119Chapter 6Summary and Conclusions6.1 SummaryThis thesis presents techniques aimed at reducing the truncation error at the bound-ary specifically, and the discretization and overall truncation error in general. Afinite volume formulation was used to solve the governing equations with a vertex-centered control volume approach for triangular and mixed meshes. The compar-ison of truncation error was made possible by defining an error metric. This met-ric was based on the magnitude of coefficients associated with each of the spatialderivatives of the Taylor series expansion of the truncation error for an arbitrarysolution. The discretization error was calculated by using the Method of Manu-factured Solutions, where a set of three exact solutions were used throughout thethesis for a robust comparison.The techniques have been discussed for a diffusion problem, using Poisson’sequation as the model. A sequence of triangular, quad, and mixed meshes withquad boundary layers and triangular interiors were considered for analyzing theeffects of different schemes developed to reduce truncation and discretization er-rors. The second order least-squares and the Green-Gauss solution reconstructionmethods were described in detail. The subtleties of using the vertex as a referencepoint instead of the centroid for the control volume were also highlighted for theGreen-Gauss methods.The schemes proposed for reducing error were evaluated systematically usinganalytical and numerical approaches. The analytical approach gives an algebraicrepresentation for the behavior of schemes with respect to truncation error. It is adirect extension of the Taylor series analysis used for structured meshes. Regulartopological meshes, both triangular and mixed, were used for evaluating differentschemes for asymptotic order of convergence and magnitude of errors. These testsyield the general form of the truncation error for different scenarios and serve asthe first screening process to eliminate schemes which perform poorly for system-atically perturbed meshes. The tests were performed on a family of systematicallyperturbed meshes to make the results applicable for general scenarios.The use of analytical tests is restricted to regular topological meshes with fixedconnectivity and is extended to more general and unstructured meshes through1206.2. Conclusionsnumerical testing. The flux integral for Poisson’s problem is expressed as a linearcombination of control volume averages. Evaluation of the explicit Jacobian isinstrumental for determining the coefficients of stencil control volume averagesin the flux integral. Numerical tests compared the truncation error for meshesusing different discretization and boundary treatment methods. These schemeswere compared for their performance with respect to truncation and discretizationerrors and trade-offs were sought where both the types of errors can be minimizedsimultaneously.The discretization schemes tested differ in their choice of solution reconstruc-tion methods, methods of face gradient calculation and inclusion of jump terms.The use of boundary constraints and inclusion of a second derivative in the normaldirection as an unknown are the means of including boundary condition throughreconstruction. Using ghost cells to apply the boundary condition, increasing sten-cil size to offer support for reconstruction of solution and the use of single and twojump term were the other boundary condition implementation methods tested.The triangular mesh sequence provides a framework for testing the schemes inthe most common mesh format and mixed meshes explore the inclusion of quadlayers in the boundary. The interface layer for a mixed mesh is a region of in-terest as well. Both isotropic and stretched boundary layers were tested for puretriangular, quad and mixed meshes.6.2 ConclusionsThe objective of attaining simultaneous reduction of truncation and discretizationerror for general isotropic and stretched unstructured meshes is accomplished sat-isfactorily. A wide range of schemes were applied to implement the boundary con-dition either through improving the solution reconstruction, increasing stencil size,using ghost cells or jump terms as alternate approaches. The following discussioncompares and contrasts the analytical and numerical results for the triangular, quadand mixed meshes.Analytical tests showed the gradients to have lower error when the least-squaresmethod was used in comparison to the Green-Gauss method and also demonstratedthat augmenting stencil size was not fruitful beyond a point. It was also shownthat strong boundary condition (SBC) and 4-Unknown models (4Unk) have bet-ter gradients at the boundary compared to weak boundary condition (WBC) forthe least-squares method. It was noted that using vertices as reference points forcalculating gradients by the Green-Gauss method yielded zero order gradients andfirst order solution convergence whereas using centroids as reference points pro-duced the much desired first order gradients and second order convergence for the1216.2. Conclusionssolution.Analytic tests determined that the truncation error for triangular, quad andmixed meshes is at best zeroth order for randomly perturbed meshes in both theinterior and boundary regions as well as for interface layers on mixed meshes.The potential of jump terms to increase the effective order of accuracy for trunca-tion errors on uniform meshes was also highlighted. Only the uniform mesh quadboundary layers showed an improvement in accuracy of truncation error by oneorder for positive jump coefficients; this was not possible for uniform triangularmeshes. The interface layer for uniform meshes, although in the interior region,has zero order truncation error accuracy for any discretization scheme. Secondorder truncation error accuracy was observed only for interior control volumes onuniform, scaled, and sheared meshes across different element types.The general form of truncation error revealed decoupled flux integrals for quadelements which was remedied by a gradient correction scheme. Both the least-squares and the Green-Gauss methods were prone to decoupling and particularlyfor least-squares, this led to divergence in numerical tests in the absence of a gra-dient correction.Implementation of different boundary treatment methods, including jump termmethods and gradient correction schemes, have shown second order convergencefor discretization error, across different mesh types and exact solution cases. Usingmixed cell meshes, the interface layer was shown to be similar to interior regionswith regards to truncation and discretization error. The errors were marginallyhigher for the interface layer when compared to interior regions with regards totruncation or discretization errors but were much less than the boundary regions.The face gradient calculation for the least-squares method using arithmetic av-eraging and linear interpolation showed little difference for discretization or trun-cation error. The Green-Gauss method had reasonable solution accuracy for trian-gular, quad and mixed meshes while the truncation errors were higher compared toother schemes for all element types.The numerical results for all the mesh sequences showed significant reduc-tion in discretization and interior truncation errors using the single and two jumpterm methods. The Green-Gauss and least-squares methods with SBC and gradientcorrection when needed were the next best schemes to reduce discretization errors.The boundary truncation error was best minimized when the two jump term methodwas used for implementing the boundary condition on isotropic meshes. Interiortruncation error was independent whereas discretization error was less sensitive tothe changes in the value of boundary jump coefficient. Using the single jump termmethod, the optimal jump coefficients for minimizing the interior and boundarytruncation errors had distinct values; discretization error showed slight variation inoptimal jump coefficients when the exact solutions were changed, for all the mesh1226.2. Conclusionssequences.Using two jump term method, the set of optimal jump coefficients for reducingeach of the interior truncation error, the boundary truncation error and the dis-cretization error are different. However, a Pareto plot enabled us to find a point inthe two jump term space such that all the three errors are minimized simultaneouslywithin a certain range of their minimum possible values. On isotropic triangularand mixed meshes, a two jump coefficient set of {α = 1.24,αB = 1.11} was iden-tified. For this jump coefficient set, the overall truncation and discretization errorhad a deviation of about 20-30% from their respective minimum possible valuesacross the sequence and for different exact solution cases tested. The jump set of{α = 1.36,αB = 1.11} was also identified for quad meshes but the variation fordiscretization error was about 100% whereas the truncation error was allowed todeviate by about 20% from their respective minimums.Thus using two jump term methods, the simultaneous reduction of truncationerror and discretization error for an isotropic mesh was made possible in this thesis.When mixed five stretched quad layered meshes were used, the advantage of usingtwo jump terms was lost as the single jump term scheme had the lowest errors. Thiswas in contrast to the results observed for stretched triangular boundary layeredmeshes, where two jump term scheme consistently resulted in lowest boundarytruncation errors. The analytical results for stretched quad layers also pointed atthe superior performance of two jump term scheme over single jump term method.A Pareto plot analysis for stretched five quad layered meshes was performedand the deviation from their minimum possible errors for discretization and bound-ary truncation errors using the jump coefficient set {α = 1.24,αB = 0.84} wereabout 100% and 35% respectively. This scenario using two jump term method wasbetter than using the single jump term method to reduce errors simultaneously. Asimilar analysis for the stretched triangular boundary layered meshes, using thejump coefficient set of {α = 1.24,αB = 0.84} showed a maximum deviation of65% for discretization errors and about 10% for overall truncation error.This conclusion proves that the two jump term scheme has successfully re-duced truncation and discretization errors for any mesh type or exact solution typetested in this thesis. However, the anomaly in the performance of the two jumpterm scheme for stretched quad layered boundary meshes warrants further researchbefore any reliable suggestions can be made for the success in this category ofmeshes.For the isotropic and stretched mesh families, the truncation errors near theboundary with quad layers were in general lower than in the case of triangularmeshes, though the gain was not phenomenal. The interior truncation error for apure quad domain is only slightly better than that for the pure triangular meshes.This implies that there are no significant gains with regards to diffusion when an1236.3. Recommended Future Workunstructured mesh is used with either triangular or quad boundary layers. Usingstretched meshes only increased the truncation error near the boundary.In conclusion, it is noted that using the two jump term scheme has reducedtruncation errors by at least five times for isotropic triangular and mixed meshesas well as quad and stretched triangular meshes when compared to the use of SBCalone. Using the single or two jump term methods has reduced boundary truncationerrors by about 2-3 times for stretched five quad layered meshes compared to SBCmethod. The discretization errors decrease at least by a factor of 10-20 times forall the mesh types using either of the jump term schemes for the exact solutionstested. The interior truncation error was approximately halved by using jump termscompared to the no jump term scenario.6.3 Recommended Future WorkThe techniques and results presented in this thesis have answered most of the ques-tions that were posed at the onset of this research. The techniques developed in thiswork can be used to reduce different forms of errors simultaneously for isotropictriangular, quad and mixed meshes. However, the scope of research was limitedby certain factors and new avenues were identified where the research can be takenforward in order to answer questions that arose during the course of this thesis. Thefollowing are considerable ways the research can be extended1. The control volume averages for ghost cells were extrapolated from the inte-rior regions using a simple second order scheme. Effective ways of extrapo-lating information to the ghost cells with higher order accuracy is anticipatedto give better truncation errors than those observed in this work. It would beinteresting to see if this suggestion can lead to truncation errors better thanthe two jump term method.2. The boundary treatment methods were tested mostly on isotropic meshes.Performance of these methods on anisotropic meshes is a missing piece tothe puzzle in understanding the behavior of truncation and discretization er-rors.3. The analysis for the boundary truncation error can be easily extended to cell-centered control volume methods. The inclusion of boundary constraintsand jump terms for a cell-centered framework is already available in theliterature. The assessment of different schemes on cell centered methods isa considerable study that can be pursued.1246.3. Recommended Future Work4. The research presented was limited to second order methods and the exten-sion to higher order methods is crucial in understanding the behavior of theseschemes. The extension would entail the calculation of higher order flux Ja-cobians, higher order moments and higher order reconstruction methods.5. Extension of the current methodology from 2D to 3D versions will be chal-lenging yet necessary if the effect of the schemes on real-life problems areto be understood. 3D reconstruction and 3D Jacobian calculation are reason-ably straightforward, however the storage requirements could be an issue ofconcern.6. Ultimately, the analysis can be extended from a simple diffusion problem tomore general fluid flow phenomena governed by the PDEs such as the Euleror Navier-Stokes equations.125Bibliography[1] J.W. Banks, J.A.F. Hittinger, J.M. Connors, and C.S. Woodward. Numericalerror estimation for nonlinear hyperbolic pdes via nonlinear error transport.Comp. Meth. Appl. Mech. Eng., 213-216:1–15, March 2012.[2] Timothy J. Barth. Numerical aspects of computing high reynolds num-ber flows on unstructured meshes. In In Proceedings of tweny-ninth AIAAAerospace Science Meeting, January 1991. AIAA paper 91-0721.[3] Timothy J. Barth and Paul O. Frederickson. Higher order solution of the Eulerequations on unstructured grids using quadratic reconstruction. AIAA paper90-0013, January 1990.[4] J. Blazek. Computational Fluid Dynamics: Principles and Applications. El-sevier Academic Press, 2001.[5] Ryan Bond, Patrick Knupp, and Curtis Ober. A manufactured solution forverifying CFD boundary conditions, part II. 2005. 43rd AIAA AerospaceSciences Meeting and Exhibit. January.[6] Andrew Cary, Andrew Dorgan, and Mori Mani. Towards accurate flow pre-dictions using unstructured meshes. In Fluid Dynamics and Co-located Con-ferences. American Institute of Aeronautics and Astronautics, June 2009.AIAA paper 2009-3650.[7] William Coirier and Philip Jorgenson. A mixed volume grid approach for theeuler and navier-stokes equations. 1996. 34th Aerospace Sciences Meetingand Exhibit. January.[8] AE Dahoe and RS Cant. On least-squares gradient reconstruction and itsapplication in conjunction with a Rosenbrock method. 2004.[9] B. Diskin and J.L. Thomas. Accuracy analysis for mixed-element finite-volume discretization schemes. Technical Report 2007-08, National Instituteof Aerospace, 2007.126Bibliography[10] B. Diskin, J.L. Thomas, E.J. Nielsen, H. Nishikawa, and J.A. White. Compar-ison of node-centered and cell-centered unstructured finite-volume discretiza-tions. Part I: Viscous fluxes. In Proc. 47th Aero. Sci. Meeting, 2009. AIAApaper 2009-597.[11] Boris Diskin and James Thomas. Accuracy of gradient reconstruction ongrids with high aspect ratio. Technical report, National Institute of Aerospace,August 2007. NIA Report 2007-08, National Institute of Aerospace.[12] Boris Diskin and James Thomas. Notes on accuracy of finite-volume dis-cretization schemes on irregular grids. Applied numerical mathematics,60:224–226, 2010. Rebuttal of Svard 2008.[13] Peter Eliasson. EDGE: A Navier-Stokes solver for unstructured grids. Tech-nical report, FOI-R-0298-SE, 2001.[14] Joel H Ferziger and Milovan Peric´. Computational methods for fluid dynam-ics, volume 3. Springer Berlin, 2002.[15] Andreas Haselbacher and Jiri Blazek. Accurate and efficient discretizationof navier-stokes equations on mixed grids. AIAA journal, 38(11):2094–2102,2000.[16] Andreas Haselbacher, James J. McGuirk, and Gary J. Page. Finite volumediscretization aspects for viscous flows on mixed unstructured grids. AIAA J.,Vol. 37(No. 2):177–184, February 1999.[17] Alireza Jalali. Truncation error analysis of unstructured finite volume dis-cretization schemes. Master’s thesis, University of British Columbia, May2012.[18] Alireza Jalali and Carl Ollivier-Gooch. Accuracy assessment of finite volumediscretizations of diffusive fluxes on unstructured meshes. In Aerospace Sci-ences Meetings. American Institute of Aeronautics and Astronautics, January2012.[19] Alireza Jalali and Carl Ollivier-Gooch. Accuracy assessment of finite volumediscretizations of convective fluxes on unstructured meshes. In AerospaceSciences Meetings. American Institute of Aeronautics and Astronautics, Jan-uary 2013.[20] Alireza Jalali and Carl Ollivier-Gooch. Higher-order finite volume solutionreconstruction on highly anisotropic meshes. June 2013. 21st AIAA Compu-tational Fluid Dynamics Conference. June.127Bibliography[21] Kelly R. Laflin, Steven M. Klausmeyer, Tom Zickuhr, John C. Vassberg,Richard A. Wahls, Joseph H. Morrison, Olaf P. Brodersen, Mark E. Rakowitz,Edward N. Tinoco, and Jean-Luc Godard. Data Summary from Second AIAAComputational Fluid Dynamics Drag Prediction Workshop. 2005. Journal ofAircraft 2005 42:5 , 1165-1178.[22] Randall J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cam-bridge University Press, 2004.[23] David W. Levy, Thomas Zickuhr, John Vassberg, Shreekant Agrawal,Richard A. Wahls, Shahyar Pirzadeh, and Michael J. Hemsch. Data Sum-mary from the First AIAA Computational Fluid Dynamics Drag PredictionWorkshop. 2003. Journal of Aircraft 2003 40:5 , 875-882.[24] D.W. Levy, Kelly R. Laflin, Edward N. Tinoco, John C. Vassberg, Mori Mani,Ben Rider, Christopher L. Rumsey, Richard A. Wahls, Joseph H. Morrison,Olaf P. Brodersen, Simone Crippa, Dimitri J. Mavriplis, and Mitsuhiro Mu-rayama. Summary of Data from the Fifth AIAA CFD Drag Prediction Work-shop. In Proceedings of the Fifty-First AIAA Aerospace Sciences Meeting,2013.[25] Dimitri Mavriplis. Revisiting the least-squares procedure for gradient recon-struction on unstructured meshes. In Fluid Dynamics and Co-located Con-ferences. American Institute of Aeronautics and Astronautics, June 2003.[26] Dimitri Mavriplis. Results from the 3rd Drag Prediction Workshop using theNSU3D unstructured mesh solver. In Aerospace Sciences Meetings. Amer-ican Institute of Aeronautics and Astronautics, January 2007. 45th AIAAAerospace Sciences Meeting and Exhibit.[27] Dimitri J. Mavriplis, Antony Jameson, and L. Martinelli. Multigrid solutionof the Navier-Stokes equations on triangular meshes. ICASE Report No. 89-11, NASA Langley Research Center, 1989.[28] Dimitri J. Mavriplis, John C. Vassberg, Edward N. Tinoco, Mori Mani, Olaf P.Brodersen, Bernhard Eisfeld, Richard A. Wahls, Joseph H. Morrison, TomZickuhr, and David Levy. Grid quality and resolution issues from the DragPrediction Workshop series. In 46th Aero. Sci. Meet., 2008.[29] Christopher Michalak. Efficient high-order accurate unstructured finite-volume algorithms for viscous and inviscid compressible flows. PhD thesis,U. Brit. Col., 2009.128Bibliography[30] Hiroaki Nishikawa. Beyond interface gradient: A general principle for con-structing diffusion scheme. In Proceedings of the Fortieth AIAA Fluid Dy-namics Conference and Exhibit, 2010. AIAA paper 2010-5093.[31] Carl F. Ollivier-Gooch. High-order ENO schemes for unstructured meshesbased on least-squares reconstruction. AIAA paper 97-0540, January 1997.[32] Carl F. Ollivier-Gooch. Quasi-ENO schemes for unstructured meshes basedon unlimited data-dependent least-squares reconstruction. J. Comp. Phys.,133(1):6–17, 1997.[33] Carl F. Ollivier-Gooch. A toolkit for numerical simulation of PDE’s: I. Fun-damentals of generic finite-volume simulation. Comp. Meth. Appl. Mech.Eng., 192:1147–1175, February 2003.[34] Carl F. Ollivier-Gooch and Michael Van Altena. A high-order accurate un-structured mesh finite-volume scheme for the advection-diffusion equation.J. Comp. Phys., 181(2):729–752, 2002.[35] Christopher Roy. Strategies for driving mesh adaptation in CFD. In Proc.47th Aero. Sci. Meeting. Amer. Inst. Aero. Astro., 2009.[36] Christopher Roy. Review of discretization error estimators in scientific com-puting. In Proc. 48th Aero. Sci. Meeting, 2010.[37] Kambiz Salari and Patrick Knupp. Code verification by the method of man-ufactured solutions. Technical report, Sandia National Laboratories, June2000. Sandia Report.[38] D. Sharov and K. Nakahashi. Hybrid prismatic/tetrahedral grid generationfor viscous flow applications. 1998. AIAA Journal 1998 36:2 , 157-162.[39] P.R. Spalart and S. R. Allmaras. A one-equation turbulence model for aero-dynamic flows. In Proceedings of the Thirtieth AIAA Aerospace SciencesMeeting and Exhibit, January 1992. AIAA paper 92-0439.[40] James L. Thomas, Boris Diskin, and Christopher Rumsey. Towards verifi-cation of unstructured-grid solvers. In 46th Aero. Sci. Meet., 2008. AIAA2008-666.[41] Jiyuan Tu, Guan Heng Yeoh, and Chaoqun Liu. Computational fluid dynam-ics: a practical approach. Butterworth-Heinemann, 2007.129Bibliography[42] Lars Tysell. Hybrid grid generation for viscous flow computations aroundcomplex geometries. Technical report, 2010. Department of Mechanics,Royal Institute of Technology.[43] Michael Van Altena. High-order finite-volume discretisations for solvinga modified advection-diffusion problem on unstructured triangular meshes.Master’s thesis, U. Brit. Col., October 1999.[44] H.K. Versteeg and W. Malalasekera. An Introduction to Computational FluidDynamics: The Finite Volume Method. Pearson Education Limited, 2007.[45] Gary Yan, Varun Prakash Puneria, Alireza Jalali, and Carl Ollivier-Gooch.Truncation and discretization error for diffusion schemes on unstructuredmeshes. In AIAA SciTech. American Institute of Aeronautics and Astronau-tics, January 2014.[46] X.D. Zhang, J.-Y. Trépanier, and R. Camarero. A posterior error estimationfor finite-volume solutions of hyperbolic conservation laws. Comp. Meth.Appl. Mech. Eng., 185:1–19, 2000.130