Improving Finite-Volume DiffusiveFluxes Through Better ReconstructionbyChandan Balachandra SejekanB.Tech., Mechanical Engineering, National Institute of Technology Karnataka,Surathkal (India), 2013A 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 2016c© Chandan Balachandra Sejekan 2016AbstractThe overarching goal of CFD is to compute solutions with low numerical error. For finite-volume schemes, this error originates as error in the flux integral. For diffusion problemson unstructured meshes, the diffusive flux (computed from reconstructed gradients) is oneorder less accurate than the reconstructed solution. Worse, the gradient errors are notsmooth, and so no error cancellation accompanies the flux integration, reducing the fluxintegral to zero order for second-order schemes. Our aim is to compute the gradient andflux more accurately at the cell boundaries and hence obtain a better flux integral fora slight increase in computational cost. We propose a novel reconstruction method andflux discretization to improve diffusive flux accuracy on cell-centred, isotropic unstructuredmeshes. Our approach uses a modified least-squares system to reconstruct the solutionto second-order accuracy in the H1 norm instead of the prevalent L2 norm, thus ensuringsecond-order accurate gradients. Either circumcentres or containment centres are chosen asthe control-volume reference points based on a criteria to facilitate calculation of second-order gradients at flux quadrature points using a linear interpolation scheme along witha high-accuracy jump term to enhance stability of the system. Numerical results show asignificant improvement in the order of accuracy of the computed diffusive flux as well asthe flux integral. When applied to a channel flow advection-diffusion problem, the schemeresulted in an increased order of accuracy for the flux integral along with gains in solutionaccuracy by a factor of two. The characteristics of the new scheme were studied throughstability, truncation error and cost analysis. The increase in computational costs weremodest and affordable. The behaviour of the scheme was also tested by implementing avariation of it within the ANSYS Fluent discretization framework.iiPrefaceThe research ideas and methods discussed in this thesis are the fruits of a close workingrelationship between Dr. Carl Ollivier-Gooch and Chandan Balachandra Sejekan. Theimplementation of methods, the data analysis and the manuscript preparation were doneby Chandan Balachandra Sejekan with invaluable guidance from Dr. Carl Ollivier-Goochthroughout the process.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication - The Master Voyage . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 The Finite-Volume Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1.1 Mesh Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1.2 Solution Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . 142.1.3 Flux Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.1.3.1 Diffusive Flux Evaluation . . . . . . . . . . . . . . . . . . 222.1.3.2 Face Gradients . . . . . . . . . . . . . . . . . . . . . . . . 23ivTable of Contents2.1.3.3 Jump Terms . . . . . . . . . . . . . . . . . . . . . . . . . . 242.1.4 Flux Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.1.5 Time Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.2 Flux Error Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273 H1 Finite-Volume Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.1 H1 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.2 H1 Flux Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2.1 Theory in One Dimension . . . . . . . . . . . . . . . . . . . . . . . 313.2.2 Extension to Two Dimensions . . . . . . . . . . . . . . . . . . . . . 323.2.2.1 Circumcentres . . . . . . . . . . . . . . . . . . . . . . . . . 333.2.2.2 Containment Centres . . . . . . . . . . . . . . . . . . . . . 353.2.2.3 High-Accuracy Jump Term . . . . . . . . . . . . . . . . . . 364 Analytical Tests and Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 384.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.2 Analytical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.2.1 Linear Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 394.2.2 Quadratic Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 434.2.3 H1 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454.2.3.1 Regular Stencil H1 Reconstruction . . . . . . . . . . . . . 454.2.3.2 Extended Stencil H1 Reconstruction . . . . . . . . . . . . 475 Numerical Tests and Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 495.1 Preliminary Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 505.1.1 One-dimensional H1 Flux Discretization . . . . . . . . . . . . . . . . 505.1.2 Two-dimensional H1 Reconstruction Tests . . . . . . . . . . . . . . 515.2 2D Poisson Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.3 2D Advection-Diffusion Test . . . . . . . . . . . . . . . . . . . . . . . . . . 576 Characterization Tests and Behavioural Analysis . . . . . . . . . . . . . . 616.1 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61vTable of Contents6.2 Truncation Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 656.3 Cost Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 696.4 H1 Finite-Volume Scheme in ANSYS Fluent Discretization . . . . . . . . . 736.4.1 Fluent Diffusive Flux Discretization . . . . . . . . . . . . . . . . . . 736.4.1.1 Reconstruction Gradients . . . . . . . . . . . . . . . . . . 736.4.1.2 Diffusion Gradients . . . . . . . . . . . . . . . . . . . . . . 746.4.1.3 Diffusive flux formulation . . . . . . . . . . . . . . . . . . 756.4.1.4 Boundary Fluxes . . . . . . . . . . . . . . . . . . . . . . . 766.4.2 H1 Scheme in Fluent Diffusive Flux . . . . . . . . . . . . . . . . . . 766.4.3 Comparative Numerical Test Results . . . . . . . . . . . . . . . . . 776.4.3.1 2D Poisson Test . . . . . . . . . . . . . . . . . . . . . . . . 786.4.3.2 2D Advection-Diffusion Test . . . . . . . . . . . . . . . . . 807 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 857.3 Recommended Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . 87Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89viList of Tables2.1 Flux Evaluation Strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.2 Gauss-point locations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.1 Stencils and truncation error for Laplacian on a cell-centred mesh . . . . . . 405.1 Orders of accuracy for Poisson model problem . . . . . . . . . . . . . . . . . 535.2 Orders of accuracy for advection-diffusion problem . . . . . . . . . . . . . . 606.1 Wall times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69viiList of Figures1.1 Stages of CFD evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Mesh generation approaches over a two-dimensional airfoil . . . . . . . . . . 31.3 Order of reconstruction accuracy of unstructured least-squares reconstructionschemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.4 Solution error comparison of second and third order least-squares reconstruction 82.1 Control volume illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 High-order reconstruction stencil for cell centred mesh type . . . . . . . . . 172.3 Mesh geometrical features . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4 Gauss quadrature points schematic . . . . . . . . . . . . . . . . . . . . . . . 263.1 1D Non-uniform Grid Stencil . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2 Various control-volume reference centres . . . . . . . . . . . . . . . . . . . . 333.3 Circumcentre cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.4 The obtuse containment centre . . . . . . . . . . . . . . . . . . . . . . . . . 364.1 Analytical Mesh coordinates and cell labels . . . . . . . . . . . . . . . . . . 394.2 Second-order least-squares on uniform equilateral triangle mesh (decoupling) 424.3 Third-order least-squares on uniform equilateral triangle mesh . . . . . . . . 444.4 Regular H1 reconstruction on uniform equilateral triangle mesh (decoupling) 464.5 Extended H1 reconstruction on uniform equilateral triangle mesh (decoupling) 485.1 Isotropic unstructured meshes used . . . . . . . . . . . . . . . . . . . . . . . 505.2 Preliminary Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515.3 Manufactured solutions for Poisson equation on an unstructured mesh . . . 52viiiList of Figures5.4 Case (i) Sinusoidal Poisson 2D test results . . . . . . . . . . . . . . . . . . . 545.5 Case (ii) Gaussian Bump Poisson 2D test results . . . . . . . . . . . . . . . 555.6 Case (iii) Plain Laplacian Poisson 2D test results . . . . . . . . . . . . . . . 565.7 Solution error comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.8 Advection-diffusion solution (α = 0.01) . . . . . . . . . . . . . . . . . . . . . 585.9 Advection-Diffusion 2D test results . . . . . . . . . . . . . . . . . . . . . . . 596.1 Eigenvalue plots - centroid vs. circumcentres . . . . . . . . . . . . . . . . . 626.2 H1 reconstruction eigenvalue plots . . . . . . . . . . . . . . . . . . . . . . . 636.3 H1 reconstruction on equilateral mesh . . . . . . . . . . . . . . . . . . . . . 646.4 Advention-diffusion eigenvalue plots . . . . . . . . . . . . . . . . . . . . . . 656.5 Truncation error convergence rate . . . . . . . . . . . . . . . . . . . . . . . . 676.6 Truncation error magnitude plots . . . . . . . . . . . . . . . . . . . . . 686.7 Jacobian calculation comparisons . . . . . . . . . . . . . . . . . . . . . . . . 706.8 Comparison of total wall clock time required for convergence . . . . . . . . 716.9 Flux error (steady state) vs total wall clock time required for convergence . 716.10 Time profile analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 726.11 Illustration of Fluent diffusive flux terms . . . . . . . . . . . . . . . . . . . . 756.12 Gradient stencils for diffusive fluxes about face f . . . . . . . . . . . . . . . 766.13 Poisson solution convergence rates in Fluent software (Fluent-code) and Flu-ent discretization implementation in ANSLib (ANSLib-Fluent) . . . . . . . 786.14 Poisson test results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.15 Advection-Diffusion solution convergence rates in Fluent software and Fluentdiscretization implementation in ANSLib . . . . . . . . . . . . . . . . . . . 816.16 Advection-Diffusion Solution convergence rates . . . . . . . . . . . . . . . . 82ixList of SymbolsRoman Symbolseˆb unit-normal vector connecting the cell centroid and the boundary face centroideˆs unit-normal vector connecting the two adjacent cell centroidsnˆ unit normal vectorFD diffusive fluxPInv pseudo-inverseS source term~A area vector of the face~F flux vector~r distance/position vector~U solution vectorA areads infinitesimal lengthh signed normal distance to the Gauss pointt timeu velocityx, y cartesian coordinatesxList of SymbolsGreek Symbolsα jump-term coefficient, Diffusion coefficient constant with a value of 10−9η normal direction coordinateΩ size of control volume (area in 2D and volume in 3D)φ control-volume-averaged value of the variableφ unknown variable, exact solutionφ˜ reconstructed value of the variableSuperscripts+ left side evaluation- right side evaluationm, n polynomial exponentsSubscripts2ndLS second-order least-squares reconstruction3rdLS third-order least-squares reconstructionCV control volumeg Gauss-point evaluationH1 proposed H1 reconstructioni original control volumej neighbouring control volumeL left CV reference point valueR right CV reference point valuexiAcknowledgementsI owe my gratitude to all those people who have made this thesis possible and because ofwhom my graduate experience has been one that I will cherish forever.It is difficult to overstate my appreciation of my research supervisor Dr. Carl Ollivier-Gooch, whose constant support, patience and guidance has been invaluable throughout thecourse of this research. I count myself lucky to have him as a supervisor and mentor. I havelearned much from him not just in academics and research but also valuable life lessons inthe virtues of patience, compassion and humility.I would also like to acknowledge my colleagues at ANSLab research group Gary, Varunand Alireza for bearing with me and clarifying even the most trivial of doubts. Withouttheir support, my progress would have been much slower. My thanks also go to Reza andMahkame for proof-reading this document.I am highly grateful to my friends for all the wonderful trips and fond memories. Theyare the spice that made my journey less humdrum.Last but most importantly, I wish to thank my beloved family back in India for theirunceasing, unwavering love and encouragement. They have been instrumental in mouldingme into who I am today. It is to them that I dedicate my thesis.xiiThe Master VoyageRainy welcome to a city known anew.Long arduous journey cast asideWith aspirations and zeal renew.Stepping out boldly astride,Across uncharted lands, uncharted life.Friends from the past reckon.Friends from the future beckon.Friends galore together striveThrough highs and lows, a hive,Over uncharted lands, uncharted life.Homes of fun and laughter, incompleteWithout fond family far far away.Constant craving to carve out a tunnelThrough the Earth’s infernal core,Through uncharted lands, uncharted life.Travel and learning hand in hand.New sights, hills, and mountains,Knowledge to explore and surmount.Aspiring to grow steadfast,Learning uncharted lands, uncharted life.Innumerable steps this voyage take.Countless crowds cater this eager soul,A morsel here, a feast there.With all their unwavering support, I purportConquering uncharted lands, uncharted life.- Chandan B. Sejekan23rd March, 2016xiiiChapter 1IntroductionComputational Fluid Dynamics (CFD) has grown tremendously in the past decades to playa vital role in the study and design of fluid flows. Thanks to the accessibility of affordablehigh-performance computing hardware, CFD has become a mainstay in industries fromcivil construction to aerospace. It thus forms the third front in fluid analysis along withtheoretical and experimental analysis. In fact CFD is preferred in the design fields due toits significant advantages of quick turn-around times, greater level of detail and flexibilityas well as significantly lower costs of designs.Most CFD analysis can be broken down to comprise of three basic stages of evalua-tion [57]; pre-processing, flow solver run and post-processing. We will briefly touch on eachof these stages so as to better identify where the current research resides in the whole schemeof things. Figure 1.1 presents the same information in a schematic flow chart.Pre-processing: The pre-processing stage prepares all of the necessary inputs that arerequired by the flow solver stage, in a format that can be understood by it. This stageessentially comprises of physical modelling and mesh generation. In the former, we math-ematically define or model the physics of the flow problem that needs to be solved. Thisentails identification of the physical and chemical phenomena to be modelled and definingthe geometry of the computational domain as well as its boundary conditions. We are essen-tially writing mathematical governing equations, typically in the form of partial differentialequations (PDEs), that define the flow physics or approximate them. The error that arisesdue to this approximation of the flow physics is known as physical modelling error. Thisform of error varies greatly depending on the flow being analyzed.The boundary conditions are equally important and should not be neglected. Theydirectly influence the solution obtained at every stage. There are two common classifications1Chapter 1. IntroductionFigure 1.1: Stages of CFD evaluationof boundary conditions: Dirichlet and Neumann. Dirichlet boundary conditions directlyspecify a value for the given variable at the boundary, for instance isothermal wall orconstant velocity inlet conditions. Neumann conditions, on the other hand, specify a normalgradient for that variable at the boundary like the adiabatic wall or outlet conditions. Thethird and rare type, the Robin boundary condition is usually an expression which is acombination of the Dirichlet and Neumann boundary types. Its implementations can befound in convective heat transfer models.Once the geometry of the computational domain has been defined, we can generate amesh for it, which essentially involves decomposing the domain into finite non-overlappingsub-domains. These sub-domains are commonly known as cells or control volumes. Theflow solver discretizes the governing equations on these individual sub-domains. Thus thefiner these sub-domains are the more accurate the solution will be, while also increasingthe computational costs. The mesh generation can follow two approaches, structured orunstructured. A structured mesh, as the name suggests, has a regular topology, whichmakes identification of each cell predictable and implicit. For instance the cell (i, j, k)would always have (i ± 1, j ± 1, k ± 1) as its neighbours. This implicit cell connectivitymakes it easier for the flow solver to process structured meshes. However, meshing ofcomplex geometries is virtually impossible to automate for all possible geometries and thusdemands a significant chunk of total CFD analysis time. In contrast unstructured meshes2Chapter 1. Introduction(a) Structured meshing (b) Unstructured meshingFigure 1.2: Mesh generation approaches over a two-dimensional airfoildo not enforce a regular topology. This makes it much easier to mesh any geometry as anarbitrary combination of cells, commonly using triangles in two dimensions and tetrahedrain three dimensions. This assures flexibility and automation in meshing. However, now thecell connectivity needs to be explicitly defined which makes writing efficient unstructuredflow solvers all the more challenging. Fortunately with advances with solver techniques, useof unstructured meshes has become common and proves to be more successful for complexproblems [56]. These two approaches to mesh generation have been illustrated in Figure 1.2.Flow Solve: With the mathematical model equations describing the flow field finalized,we need to discretize these PDEs within each cell, to enable their numerical approximationto the exact solution. A wide variety of discretization schemes and methods have beendevised by researchers over the years. Broadly they fall under three main families; finite-difference, finite-element (and spectral element) and finite-volume schemes.The common outline binding all three families is,1. Discretize (approximate) the governing equations, in space and time if necessary, intoa system of algebraic equations,2. Solve these algebraic equations by some iterative method.The approximation in the first step leads to the error known as the truncation error, whichis the difference between the discrete algebraic equations and the continuous PDEs. The3Chapter 1. Introductionlatter step results in the final numerical solution and an associated error known as thediscretization error which can be defined as the difference between the exact solution of thegoverning equations and the numerical solution obtained from the scheme. As might beexpected, the discretization error is related to the truncation error by the error transportequation, as shown by Roy [49]. Roy also highlighted the importance of truncation error forerror estimation and also its applications in mesh adaptation to reduce the overall error [48].The most common approach to reducing the discretization error to improve the numericalaccuracy is either to use a finer mesh (h-refinement) or to increase the order of accuracy(p-refinement).The finite-difference approach has historical value being the very first numerical tech-nique to be developed for differential equations. It relies on using Taylor series expansionsto approximate the point-wise values of the solution variables. It is simple and mathemat-ically intuitive. However, conservation of variables needs to be explicitly handled; also itsextension to unstructured meshes is quite a challenge. Excellent overviews of the finite-difference process can be found in the books [2, 27] with numerous works in the literatureon advanced topics such as [45, 24].The finite-element techniques finds its origins in the structural analysis of solids and waslater applied to fluids. It encompasses methods such as the continuous Galerkin, streamline-upwind/Petrov–Galerkin (SUPG)[9, 16], discontinuous Galerkin [15, 10], spectral volume,spectral difference [23, 26] and more recently flux reconstruction [17, 58] methods. In thesemethods the solution is approximated by polynomial basis functions of the required degreewithin each cell, and the approximate solution is allowed to be discontinuous at the cellinterfaces. Solution between neighbouring cells is coupled via the use of a common numericalflux applied at these cell interfaces. It is easily extended to unstructured meshes. However,the conservation of properties is still not guaranteed, and is particularly a challenge in flowswith discontinuities such as shock waves.The finite-volume scheme, which is at the heart of our research, takes the control-volumeanalysis approach commonly used in thermofluids systems. Here we integrate the conser-vative form of the governing equations over each control volume. This step distinguishesthe finite-volume method from the others since the resulting numerical scheme is inher-4Chapter 1. Introductionently conservative. This conservative property makes it easy to understand, intuitive andattractive for flows with discontinuities. Once again there are several schemes that fallunder this umbrella such as the Green-Gauss, least-squares reconstruction (also known ask-exact reconstruction), essentially non-oscillatory schemes (ENO) [14, 1] and its weightedmodifications (WENO) [25, 12]. A detailed description of the most commonly implementedfinite-volume scheme, the least-squares technique, will follow in Chapter 2.Post-processing: The final stage is the presentation and analysis of the numerical so-lution obtained. This generally involves visualization in the form of vector plots, surfaceplots, contour plots, and the like. There are innumerable plotting tools at one’s disposalamong which we found gnuplot, matplotlib (python) and Paraview to be useful for us.Software Frameworks UsedThe Advanced Numerical Simulations Library (ANSLib) [41], a CFD tool developed andmaintained at the Advanced Numerical Simulations Laboratory (ANSLab) in the Universityof British Columbia, Vancouver, was used as the primary test bed in this research. It is anunstructured finite-volume flow solver that is capable of handling two and three dimensionalmeshes with accuracies ranging from first to fourth order for a variety of flow physicsproblems.All the meshes used in this research were generated using the meshing software, Gener-ation and Refinement of Unstructured Mixed-element Meshes in Parallel (GRUMMP) [40]also developed at ANSLab, UBC. It is a general purpose automatic mesh generation softwarefor unstructured meshes with mixed element types in two and three dimensions. The soft-ware can produce high-quality meshes that meet user-defined mesh density requirements,using elements appropriate for the geometry and physics of a particular problem. It alsosupports mesh refinement, optimization, coarsening and local remeshing.51.1. Motivation1.1 MotivationFirst and second order finite-volume schemes (commonly known as low order) are widelyutilized in the industrial and commercial Computational Fluid Dynamics (CFD) codes fortheir robustness and algorithm simplicity. However they fail to satiate the ever increas-ing demands of accuracy expected from CFD algorithms. High-order methods (having aspatial accuracy of third-order or higher) can help meet this demand for accuracy. Vari-ous researchers at the First (2012), Second (2013) and Third (2015) International Work-shops on High-Order CFD Methods suggest better performance and accuracy comparedto second-order methods with lower computational costs for the same target solution er-ror [5, 51, 59, 18]. High-order schemes thus show great promise. Yet industry level adoptionof high-order accurate schemes has been slow. This may be attributed to the increased al-gorithmic complexity and less robust behaviour of high-order schemes. There has beenconscientious effort by the high-order community over the years to address this issue andbridge the current chasm that exists between the technology being driven in academia andthe current algorithms adopted in practical applications.High-order finite-volume methods, presented in section 2.1, are in a unique positionto answer this calling since they are inherently an extension of the second-order finite-volume methods currently utilized in the industry. This significant advantage should maketheir wide-spread adoption easier. Unfortunately they do demand a substantial increase incomputational resources per control volume and also a greater algorithmic complexity. Thishas proved to be a major deterrent in their adoption. We hope to mitigate this situationwith the development of an intermediate finite-volume scheme that combines the best ofboth worlds, while being an easy upgrade from the existing second-order scheme.It is well known among the CFD community that a p-order-accurate solution reconstruc-tion would produce a p-order-accurate convective flux that depends only on the solution.However diffusive fluxes, since they also depend on the gradient, turn out to be of only(p-1)-order-accuracy. This can be visualized in Figure 1.3a, for which we set up a second-order finite-volume least-squares reconstruction for an unstructured mesh (details followin Chapter 2) and plotted the respective orders of accuracy of the reconstructed solution,61.1. Motivation-5-4-3-2-102 3 4 5log 10( L 2 error )log10( Number of Control Volumes )1st order2nd orderL2 GradXL2 GradYL2 SolnL2 Flux Integral(a) Second-Order-7-6-5-4-3-2-102 3 4 5log 10( L 2 error )log10( Number of Control Volumes )1st order2nd order3rd orderL2 GradXL2 GradYL2 SolnL2 Flux Integral(b) Third-orderFigure 1.3: Order of reconstruction accuracy of unstructured least-squares reconstructionschemes.the gradient and the flux integral. We can observe that although the solution is indeedsecond-order as anticipated, the gradient reconstruction is one order less accurate while theflux integral is even worse at zero-order. A similar trend can be seen in the convergencerate of a third-order scheme (Figure 1.3b), with loss of an order of accuracy as we movefrom solution to gradient to flux-integral. We aim to address this shortcoming of the fluxwith our research.Any improvements in the accuracies of the flux and gradients would be beneficial invarious applications where the physical quantities rely on the gradient evaluation such asviscous stress and heat fluxes in the Navier-Stokes equations as well as the gradient-basedsource terms that are used in many turbulence flow models such as the Spalart-Almarasturbulence model.It is also reasonable to expect that any improvement in the flux would ultimately trans-late into improvements in the flux integral as well. This in turn should cascade into improve-ments in the truncation error, which is known to be identical to the flux integral obtainedfrom smooth initial data. This would entail improvements in all of the areas where trunca-tion error estimation plays a vital role namely mesh adaptation [48], error estimation[49],71.2. Objectives(a) Second-order solution error (b) Third-order solution errorFigure 1.4: Solution error comparison of second and third order least-squares reconstructiondefect correction as well as the adjoint problem [13]. Section 6.2 will further touch on thistopic.Moreover, as Figure 1.4 shows, the error characteristics of the solution also vary greatlywith order of accuracy. The second-order and the third-order least-squares schemes showa striking difference in the smoothness of the solution discretization and truncation errors,with the latter providing remarkably smoother errors. As can be observed, the smoothersolution error also reflects the shape of the underlying solution, which in this case waspi8 sin (pix) sin (piy). A smoother solution error enables better application of various errorestimation algorithms [50].1.2 ObjectivesThe primary goal of this research is to develop a finite-volume framework that would improvethe diffusive fluxes to an order of accuracy on par with the solution order of accuracy.In particular we will focus on obtaining second-order accurate diffusive fluxes at the cellinterfaces especially at points of interest known as Gauss points at which the flux valuesare integrated using Gauss quadrature to compute the flux integral.81.3. OutlineWe also aim to achieve this at a modest computational cost which in this context shouldbe between the cost required for a full second-order finite-volume least-squares scheme andthird-order finite-volume scheme.Essentially as an general guideline, we strive to strike a balance between the best ofsecond and third order finite-volume schemes.Once developed we will look at evaluating and analyzing the characteristics of the newnumerical scheme for important qualities of robustness (or in other words, stability) as wellas the truncation error generated.1.3 OutlineWe begin by laying down a detailed foundation of the general finite-volume flow solvercurrently implemented in ANSLib in Chapter 2. We will focus on key aspects of the solverincluding the solution reconstruction and the flux discretization, which are most relevantto this research, while briefly discussing the rest namely the flux integration and the timeadvance techniques of the solver.Chapter 3 provides details of the proposed new finite-volume scheme which we havenamed the H1 finite-volume scheme, so called because our solution reconstruction is secondorder accurate in the H1 norm.Analytical tests were conducted, the results of which are presented in Chapter 4. We jux-tapose the new scheme with the second and third order finite-volume least-squares schemesand draw inferences. The analytical tests are only practical for a uniform equilateral trian-gular mesh, while the numerical tests are more versatile and give a better picture.Chapter 5 highlights the results of the numerical tests and once again the H1 finite-volume scheme is compared with the other two finite-volume schemes. We explore thecharacteristics and behaviour of the new scheme in Chapter 6, focusing on the traits ofstability, truncation error and computational cost. Finally, we conclude with Chapter 7,wherein we summarize, provide some concluding remarks and also list possible avenues forfuture research and development.9Chapter 2BackgroundThe relevant CFD practices that are essential for the understanding of the current researchwill be expounded in this chapter. The Advanced Numerical Simulations Library (ANSLib),the unstructured finite-volume flow-solver developed at the Advanced Numerical Simula-tions Lab at UBC, has been the primary test bed for this research. We begin with someinsight into the workings of a finite-volume flow-solver with emphasis on the solution re-construction and flux evaluation followed by an overview of flux integral evaluation andtime advance techniques. We will then move on to an understanding of error reduction insubsequent sections of this chapter.2.1 The Finite-Volume SolverFinite-volume methods have long been recognized as one of the most robust and reliablediscretization techniques available in CFD. As mentioned earlier, these methods inherentlyutilize the conservative form of the equations that govern fluid mechanics and as suchrequire the physical domain to be divided into a number of control volumes. There arevarious approaches to forming the individual control volumes based on a mesh but the mostcommonly implemented are cell-centred and vertex-centred.Cell-centred control volumes are formed directly from the underlying mesh grid withthe cells of the mesh defined as the control volumes. In this approach the cell centroid iscommonly chosen as the cell reference point for the cell, serving as the local origin for areconstruction polynomial described later in this chapter. The vertex-centred approach, onthe other hand, directly utilizes the mesh vertices as the cell reference points and buildscontrol volumes around them, typically using the median-dual as the cell perimeter. Themedian-dual is defined as the lines connecting the cell centroids and edge midpoints.102.1. The Finite-Volume Solver(a) Cell centered control volume (b) Vertex centered control volumeControl Volume Centroid Median DualPrimal MeshFigure 2.1: Control volume illustrationFigure 2.1 provides a graphical illustration of the two common types of control volumes.Both forms have their inherent advantages and disadvantages. An inquisitive reader isdirected to the work of Diskin et al. [11] for an in depth comparison of the two. In ourresearch we have primarily focused on the two-dimensional cell-centred unstructured meshdomain and shall restrict our discussion to the same.The conservative form of the fluid governing equations can be cast into a generic form [2]which can be written as∂~U∂t+∇ · ~F = ~S (2.1)where ~U denotes the solution vector, ~F the flux dyad and ~S the source term. In generalour flow solver ANSLib works with the dimensionless form of the governing equations andhence in this thesis all the governing equations are cast in the dimensionless form.Taking the integral of Equation 2.1 over an arbitrary control volume CVi and utilizingthe divergence theorem we obtain a 2D finite-volume formulation of the governing equationsas ¨CVi∂~U∂tdA+˛∂CVi~F · nˆds =¨CVi~SdA (2.2)where nˆ is the outward unit normal vector and ds is the infinitesimal for the line integralalong the control volume edges. Put simply, this equation says that the rate of change ofthe unknown variable within a control volume should be equal to the net flux that crosses112.1. The Finite-Volume Solverall its boundaries. This concept lies at the very heart of the finite-volume method.Equation 2.2 can be further simplified upon assuming the discretized physical domainto be constant. The solution vector ~U can then be extracted as its control volume average.Thus we havedUidt+1Ai˛∂CVi~F · nˆds = 1Ai¨CVi~SdA (2.3)where Ai is the area of the control volume. The first term on the left-hand side is thetime derivative of the averaged solution vector U¯ and the second term is known as the fluxintegral or the residual, representing the spatial discretization of the control volume. Theright hand side is the source term integral which is usually added on to the flux integralbefore advancing the solution vector in time.For the purpose of our research, the steady-state two-dimensional Poisson equation∇2φ =(∂2φ∂x2+∂2φ∂y2)= S (2.4)is used as the model diffusion equation, where φ is the unknown scalar variable and S isagain the source term. Here and throughout this thesis we will represent the unknownvariable by φ. The Poisson equation is a model for the steady-state diffusion equationand hence it appears in various forms within the fluid physics given by the Navier-Stokesequations. An unsteady form of the steady-state diffusion equation can be written as∂φ∂t−(∂2φ∂x2+∂2φ∂y2)= S (2.5)giving us the solution and flux vectors as~U = φ ; ~F =−∂φ∂x−∂φ∂y = −∇φ (2.6)Applying a similar procedure as above we obtain the finite-volume formulation of the122.1. The Finite-Volume SolverPoisson equationdφidt− 1Ai˛∂CVi∇φ · nˆds = 1Ai¨CViSdA (2.7)This form can now be implemented and solved for by a finite-volume flow solver.The overall finite-volume numerical scheme can be summarized in the following iterationstages [3]:1. Mesh Pre-processing: The physical grid discretization is pre-processed to ensureconsistency, define the control volumes as required and precompute a few geometricterms.2. Solution Reconstruction: In the most simplified terms, the process of obtainingthe piecewise polynomial coefficients for use in evaluating the fluxes in Equation 2.3from the control volume averages for each cell, is known as solution reconstruction.3. Flux evaluation: The numerical flux at individual Gauss-points are evaluated usingthe reconstructed solution.4. Flux integration: The flux contributions in each cell are integrated by performinga high-order accurate flux quadrature along each edge, generally by Gaussian quadra-ture.5. Time evolution: The solution is evolved in time using an adequate time steppingscheme such as Implicit Euler, Runge-Kutta etc., resulting once again in a new col-lection of cell averages.6. Repeat stages 2 through to 5 until a steady-state or a time-bound solution is obtained.Our current research presented in this thesis has its roots in the second and third stagesof the finite-volume iteration laid out above. Hence, in this background discussion we willfocus primarily on the same with a brief introduction to the other steps to provide a moreholistic understanding.132.1. The Finite-Volume Solver2.1.1 Mesh Pre-processingAlthough not an integral component of the finite-volume numerical scheme, the mesh pre-processing stage is nonetheless essential from the computational point of view since it is inthis stage that we analyze and process the physical domain discretization, more commonlyknown as the mesh.The various control volumes are formed in this stage and also their respective geometricterms such as their areas, their cell reference point locations etc. are evaluated. Thegeometric moments are also precomputed and stored to be used later during the solutionreconstruction stage.This pre-processing is only performed once at the very beginning in the flow solver andas such does not demand much computational time.2.1.2 Solution ReconstructionSolution reconstruction is the primary, as well as the most complex, component of thefinite-volume method [3]. Reconstruction provides an accurate estimate of the solution andgradients of the unknown variable, which is crucial since the flux integral in turn is derivedfrom those very estimates and hence reconstruction ultimately limits the spatial accuracyof the finite-volume solution.Numerous approaches to perform solution reconstruction have been developed and im-plemented over the years [8]. The second-order Green-Gauss method laid out by Barth [7]had been widely accepted by the CFD community for its simplicity and robustness. Barth [6,3, 4] also presented a high-order (loosely defined as third-order or more) approach in theform of k-exact reconstruction (also known as minimum energy or least-squares reconstruc-tion) capable of (k+1)th-order accuracy.The latter has been shown to be more accurate and stable with mean constraints andutilizing weighted reconstructions. The research of Ollivier-Gooch and Van Altena [43]extended this technique with boundary constraints in their finite-volume scheme for theadvection-diffusion equation. Their least-squares reconstruction scheme is described in thissection. We will, in Chapter 3, propose novel modifications to this technique for our research.142.1. The Finite-Volume SolverWe begin by fitting a piecewise polynomial representation of the solution within eachcontrol volume using a Taylor series expansion about the cell reference point (xi, yi).φ˜i(x, y) = φ|i +∂φ∂x∣∣∣∣i(x− xi) + ∂φ∂y∣∣∣∣i(y − yi)+∂2φ∂x2∣∣∣∣i(x− xi)22+∂2φ∂y2∣∣∣∣i(y − yi)22+∂2φ∂x∂y∣∣∣∣i(x− xi)(y − yi) + . . . (2.8)where φ˜i is the reconstructed solution at any point (x, y) and ∂k+lφi∂xk∂ylare its derivatives atthe reference point (xi, yi) of the control volume i. The solution and its derivatives arecollectively known as the coefficients of the Taylor polynomial.The eventual order of accuracy that the numerical scheme will be able to achieve isdictated by the number of derivatives evaluated in this Taylor series polynomial. Thusthe higher the degree of this polynomial the greater will be the order of accuracy of thesolution and truncating the polynomial to a degree k will result in at most (k+1)th-orderaccuracy of the scheme. For instance, truncating Equation 2.8 to just the first derivativesin x and y and thus obtaining a linear polynomial, will result in a second-order accuratereconstruction.φ˜i(x, y) = φ|i +∂φ∂x∣∣∣∣i(x− xi) + ∂φ∂y∣∣∣∣i(y − yi) +O(∆x2,∆y2) (2.9)The coefficients of the Taylor polynomial are evaluated while conserving the mean valueof the solution in the given control volume and also such that the reconstruction closelyapproximates the nearby neighbours’ control volume averages. This is a task best performedby casting the system as a least-squares problem. Hence this scheme is known as least-squares reconstruction [38, 39].Conservation of the mean is an essential requirement of the finite-volume method. Toensure this, a mean constraint equation is included in the system before solving. For the152.1. The Finite-Volume Solverreconstructed control volume this is given byφi =1Ai¨Viφ˜i(x, y)dA (2.10)Upon substituting for the reconstructed value from Equation 2.8 we note thatφi = φ|i +∂φ∂x∣∣∣∣ixi +∂φ∂y∣∣∣∣iyi +∂2φ∂x2∣∣∣∣ix2i2+∂2φ∂x∂y∣∣∣∣ixyi +∂2φ∂y2∣∣∣∣iy2i2+ ... (2.11)where the geometric moments have been represented asxmyni =1Ai¨Vi(x− xi)m(y − yi)ndA (2.12)In addition to conserving the mean of the reconstructed control volume, we also need tominimize the error in reconstructing the mean of the neighbouring control volumes. Thisempowers us with more equations that we can solve to obtain up to the kth derivative ofEquation 2.8 as demanded by the reconstruction accuracy we seek. It also ensures that thereconstruction is relevant globally across the physical domain. A specific set of neighboursis chosen for each control volume which forms its reconstruction stencil {Vj}i. The numberof neighbours and the neighbours themselves that make up this reconstruction stencil arechosen based on the required reconstruction accuracy and their topological distance fromthe reconstructed cell, as illustrated in Figure 2.2.Selecting control volumes that are physically close to the given cell is known as com-pact support [32]. Compact support has been shown to improve stability by emphasizingdata that should heavily influence the reconstruction. Neighbours are thus added onto thereconstruction stencil until an over-determined least-squares system is created wherein thenumber of equations is more than the unknown coefficients for the desired order of accuracy.Thus for second-order accuracy just three neighbours are sufficient while third and fourthorder accuracy demand at least six and ten neighbours respectively.In practice however the stencils are precomputed by appending an entire layer of neigh-162.1. The Finite-Volume SolverTNSNFNFirst Layer Neighbours (Second order reconstruction)Second Layer Neighbours (Third order reconstruction)Third Layer Neighbours (Fourth order reconstruction)SNTNTNTNSNTNTNTNSN SNTNSNTNFNFNFNSNTNCViFigure 2.2: High-order reconstruction stencil for cell centred mesh typebours to the stencil for each increase in order of accuracy such that the requirements arealways met with an excess. When implementing for the cell centred mesh type, layers ofneighbours are formed based on whom they share an edge with. The first neighbours arethose that share a common edge with the given control volume, second neighbours in turnshare a common edge with the first neighbours and so on. This generally results in stencilsizes of at least 3, 9 and 18 neighbours for second, third and fourth order accuracy. Doing somakes the system more robust and less susceptible to noise and non-smooth or oscillatorydata.Thus for every neighbouring control volume Vj in the reconstruction stencil we need tosatisfy as best we can a mean constraint for its control volume average. Mathematicallythis can be written asφj ≈1Aj¨Vjφ˜i(x, y)dA (2.13)As before, expanding the right hand side using the Taylor series expansion of Equa-172.1. The Finite-Volume Solvertion 2.8 we have1Aj¨Vjφ˜i(x, y)dA = φ|i +∂φ∂x∣∣∣∣i{1AjˆVj(x− xi)dA}+∂φ∂y∣∣∣∣i{1AjˆVj(y − yi)dA}+∂2φ∂x2∣∣∣∣i{12AjˆVj(x− xi)2dA}+∂2φ∂y2∣∣∣∣i{12AjˆVj(y − yi)2dA}+∂2φ∂x∂y∣∣∣∣i{1AjˆVj(x− xi)(y − yi)dA}+ ... (2.14)To simplify the moment calculation of every neighbouring control volume about thereference point of control volume i, we replace (x − xi) with (x − xj) + (xj − xi) and(y− yi) with (y− yj) + (yj − yi). Doing so effectively eliminates the otherwise complicatedmoment calculation by utilizing the known moments of each control volume about theirown reference points. These geometrical moments of every cell can easily be computed inadvance and stored during the initial stage of mesh pre-processing.1Aj¨Vjφ˜i(x, y)dA = φ|i +∂φ∂x∣∣∣∣i(xj + (xj − xi)) + ∂φ∂y∣∣∣∣i(yj + (yj − yi))+∂2φ∂x2∣∣∣∣i(x2j + 2xj(xj − xi) + (xj − xi)22)+∂2φ∂x∂y∣∣∣∣i(xyj + xj(yj − yi) + yj(xj − xi) + (xj − xi)(yj − yi))+∂2φ∂y2∣∣∣∣i(y2j + 2yj(yj − yi) + (yj − yi)22)+ ... (2.15)These mesh dependent geometric moments can be generalized in the formx̂mynij =1Aj¨Vj((x− xj) + (xj − xi))m · ((y − yj) + (yj − yi))n dA=n∑l=0m∑k=0n!l!(n− l)!m!k!(m− k)! (xj − xi)k · (yj − yi)l · xm−kyn−lj (2.16)Thus the error minimization equation for the control volume averages of every neighbour182.1. The Finite-Volume Solverwithin the stencil of the control volume i, can be written asφj ≈ φ|i +∂φ∂x∣∣∣∣ix̂ij +∂φ∂y∣∣∣∣iŷij +∂2φ∂x2∣∣∣∣ix̂2ij2+∂2φ∂x∂y∣∣∣∣ix̂yij +∂2φ∂y2∣∣∣∣iŷ2ij2+ ... (2.17)Hence the final least-squares system is formed for every control volume by collatingtogether all these neighbouring error minimization equations. This system is also formed foreach unknown variable in all the control volumes. The least-squares system, correspondingto the general form of Ax = β, can be written as1 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∂φ∂y12∂2φ∂x2∂2φ∂x∂y12∂2φ∂y2...i=φiwi1φ1wi2φ2wi3φ3...wiNφN(2.18)where the very first line is the mean constraint of the given control volume and the weightswij are used to emphasize the effect of geometrically nearby control volumes. These weightsare defined aswij =1|~rj − ~ri|n (2.19)where typically n ∈ [0, 2] and ~r is the position vector of the control-volume reference point.The choice of weighting scheme has been known to affect the conditioning and stability ofthe system [30, 29]. The inverse-distance weighting (n = 1) is used as a default for mostinviscid cases [30, 37].In ANSLib, our unstructured flow-solver, the mean constraint is eliminated by Gausselimination resulting in an unconstrained least-squares system which is then solved for192.1. The Finite-Volume Solverevery control volume using the Singular-Value-Decomposition (SVD) method. Moreover ifwe observe the least-squares system that is formed for every control volume, we notice thatthe leading matrix A on the left-hand side is purely geometrical and mesh dependent. Thuswe can easily speed up the reconstruction process by precomputing and storing the pseudo-inverse of this matrix. It can then be multiplied with the modified control-volume averages,the column matrix β, at each time step to obtain the new reconstruction coefficients.φ∂φ∂x∂φ∂y12∂2φ∂x2∂2φ∂x∂y12∂2φ∂y2...i︸ ︷︷ ︸Rec−Coeff=1 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 · · ·−1︸ ︷︷ ︸Pseudo−inverseφiwi1φ1wi2φ2wi3φ3...wiNφN︸ ︷︷ ︸CV−avgs(2.20)We will represent the pseudo-inverse in the following general form where the columnscorrespond to the neighbours and the rows correspond to the reconstruction coefficients:PInvLS =A11 A12 A13 · · · A1NA21 A22 A23 · · · A2NA31 A32 A33 · · · A3N......... . . ....(2.21)Another very crucial aspect that should be addressed before any computational schemeis complete are its boundary treatment methods. The research presented here does not deal202.1. The Finite-Volume Solverwith the boundary control volumes and hence prior knowledge of the same is not essentialto its understanding. A curious reader is directed to a detailed assessment of the same byPuneria [47, 46].2.1.3 Flux EvaluationAs depicted by Equation 2.2, to find the rate of change of the unknown variable we need toevaluate the flux integral that tells us the net flux that crosses the boundaries of the givencontrol volume. To do so we need the value of the numerical flux at the Gauss points on eachedge of the cell so that a Gauss-quadrature integration can be performed using them. Thus,an accurate flux evaluation plays a vital role in the finite-volume computational scheme.The evaluation of the numerical flux greatly depends on the physics of the problembeing solved and the corresponding discretization scheme being used. There are broadlytwo main categories of fluxes. Convective fluxes describe the direct transport of somephysical quantity across the control volume boundary and depend only on the solution.Diffusive fluxes describe processes such as mixing and heat conduction and depend on boththe solution and its gradient. This research focus being the Poisson diffusion physics overa cell-centred finite-volume approach, we will in this section only discuss flux evaluationstrategies for the same. A detailed explanation of convective fluxes can be found in thework of Jalali [19].When using structured meshes, it is common practice to directly use a central differenc-ing scheme on the cell averages from neighbouring control volumes to obtain a face gradient.Due to the inherent smoothness of a structured mesh, doing so rewards them with errorcancellations that prevent any loss of accuracy. Unfortunately the same does not hold truefor unstructured mesh solvers. The centroids which are used as cell reference points forunstructured cell-centred meshes are inconveniently not located on the perpendicular bisec-tors of the shared edge. This setback causes much harm to the accuracy of the face gradientand eventually reflects in the overall accuracy of the solver. This very concern has beenaddressed in the current research and a resolution will be presented in Chapter 3 with re-sults in subsequent chapters. For now, we will discuss the current practices in unstructuredmesh diffusive flux evaluation.212.1. The Finite-Volume SolverLRgbaFigure 2.3: Mesh geometrical features2.1.3.1 Diffusive Flux EvaluationThe diffusive flux for a face is given byFD =ˆab(∇φ)g · nˆds (2.22)Thus in the case of diffusive fluxes, a unique value of gradients needs to be evaluated atall the face Gauss points. The easiest approach would be to take advantage of the gradientcoefficients obtained from the reconstruction process explained in the previous section. Thusthe gradients at the Gauss points, approximated using a Taylor series expansion from saythe left control volume can be written as∇˜φg+ =∂φ∂x∣∣∣L+ ∂2φ∂x2∣∣∣L(xg − xL) + ∂2φ∂x∂y∣∣∣L(yg − yL) + · · ·∂φ∂y∣∣∣L+ ∂2φ∂y2∣∣∣L(yg − yL) + ∂2φ∂x∂y∣∣∣L(xg − xL) + · · · (2.23)Figure 2.3 illustrates the various geometric terms utilized in this equation. A similarexpression can be written from the perspective of the right control volume. But it should222.1. The Finite-Volume Solverbe noted that we require a unique value since the conservation property of the finite-volumemethod tells us that, across a shared edge, the out-flux from one cell should be equal to thein-flux into the neighbouring cell, that is,ˆab(∇φ)g+ · nˆds =ˆab(∇φ)g− · nˆds (2.24)This may seem straightforward, since we do have gradients at the cell reference pointsfrom the reconstruction and we just have to extrapolate the same using the Taylor seriesexpansion, but there is a catch. The two gradient values at the shared Gauss point obtainedfrom the reconstruction in the left and right cells are bound to be different. There is adiscontinuity that arises due to the inherent reconstruction errors and non-smooth data inthe reconstruction. Jalali et al. [21] have performed a comparative assessment of differentapproaches to resolve the same, including various averaging methods as well as additionof finite-difference term and jump terms for the calculation of face gradients. His findingshave been presented in the following sections.2.1.3.2 Face GradientsThe easiest way to resolve this discontinuity across the cell face would be an arithmeticaveraging of the two gradient evaluations from the left and the right cell:∇φg = 12(∇˜φg+ + ∇˜φg−)Other variations incorporate some geometrical weighted averaging scheme such as linearinterpolation or volume weighting. Linear interpolation works similar to a regular interpo-lation by taking into consideration the distance of the opposing cell reference points fromthe shared Gauss point, that is,∇φg = 1|~rLG|+ |~rRG|(|~rRG|∇˜φg+ + |~rLG|∇˜φg−)(2.25)232.1. The Finite-Volume Solverwhile the volume weighted averaging emphasizes each gradients with their correspondingvolumes (or areas in case of two-dimensions).∇φg = 1AL +AR(AL∇˜φg+ +AR∇˜φg−)(2.26)Unfortunately using just the averaged value for the flux is insufficient for a robustdiscretization scheme making the scheme susceptible to high-frequency errors [21]. The useof finite-difference terms and jump terms have been proposed in the literature to mitigatethis effect by introducing some damping into the system. Among these the use of the jumpterm has been shown to be most effective.2.1.3.3 Jump TermsThe jump term was first developed by Nishikawa [35] and was shown to be a very effectivedamping agent for high-frequency errors. It accounts for the discontinuity in the solutionacross the control volume interface and is added directly to the flux in the face normaldirection using the following formulation(∇φ)jump =α|~rLR · nˆg|(φ˜g+ − φ˜g−)(2.27)where α is a constant known as the jump term coefficient, ~rLR is the distance vector joiningthe left and right cell reference points, nˆg is the face normal unit-vector at the Gauss-pointg and φ˜g+ , φ˜g− are the solution values at the Gauss-point obtained using Taylor seriesexpansion from the left and right control volumes respectively. The jump term coefficientgreatly affects the efficacy of the jump term. Jalali [21] suggests a value of 4/3 as optimal forthe jump coefficient. Jump terms with optimal jump coefficients have been demonstratedto considerably increase the solution accuracy in the case of diffusion schemes [21, 60]. Aswill be later discussed our implementation of the jump term improves upon this currentdesign. Our findings indicate a significant improvement in the accuracy as well as stabilityof our scheme upon the addition of a jump term.The various flux evaluation strategies have been summarized in Table 2.1 with Figure 2.3242.1. The Finite-Volume SolverStrategy FormulaArithmetic averaging 12(∇˜φg+ + ∇˜φg−)Linear interpolation 1|~rLg |+|~rRg |(|~rRg|∇˜φg+ + |~rLg|∇˜φg−)Volume-weighted 1AL+AR(AL∇˜φg+ +AR∇˜φg−)Jump Term α|~rLR·nˆg |(φ˜g+ − φ˜g−)Table 2.1: Flux Evaluation Strategiesillustrating the individual geometric terms.2.1.4 Flux IntegrationA Gauss quadrature integration rule is a common choice for numerical integration sinceit requires the integrands to be evaluated at only a few specific points known as Gaussquadrature points. The greater the accuracy required, the greater is the number of Gausspoints chosen per edge. For instance, one Gauss point located at the edge-mid-side wouldsuffice for first and second-order accuracy, while third and fourth order accurate integrationwould demand two Gauss points per edge of the cell. The exact location of the Gauss pointsfor up to 4th order is given by Table 2.2. Figure 2.4 illustrates the same schematically.The accuracy of the flux integration has to be set equal to or higher than the solutionreconstruction accuracy lest all the reconstruction accuracy achieved be lost trivially.Once we have the numerical flux values at the Gauss points from the previous stage,a flux integration can be performed over that edge of the control volume using this Gaussquadrature rule. The total flux integral in the Equation 2.3 for the control volume is thenOrder of Flux Integration Location (in terms of edge vertices)Second-order ~rg =(~rB+~rA2)Third and Fourth order ~rg =(~rB+~rA2)±(~rB−~rA2√3)Table 2.2: Gauss-point locations252.1. The Finite-Volume Solver(a) Second-order (1 Gauss-point per edge) (b) Third or fourth order (2 Gauss-point per edge)Figure 2.4: Gauss quadrature points schematicsimply a summation of the individual integrals along each edge of the cell. A more in-depth analysis and explanation of the flux integration process can be found in the works ofOllivier-Gooch and Van Altena [43] and Van Altena [55].2.1.5 Time EvolutionIn the present work we are concerned with the solution to the Poisson equation which isa steady-state model equation. In such cases the time-derivative in Equation 2.3 is maderedundant and thus may be eliminated. The task then is to find the solution to the systemof equations defining the discrete flux integrals, which can be solved by any matrix solutiontechnique to obtain the steady-state Poisson solution. How well and how quickly we areable achieve that steady-state solution is described by the term, convergence of solution.Although this is essentially straight forward it greatly lacks the robustness that is oftendemanded from all flow-solvers. To circumvent this drawback ANSLib implements a tech-nique known as Pseudo-Timestepping wherein we take (implicit) time steps towards thesteady state solution instead of solving for it at one go. This imparts the desired robust-ness in the steady-state time evolution without sacrificing solution accuracy. Hence in ourresearch we use this Pseudo-Timestepping technique. Moreover, by doing we we can hopeto apply similar accuracy results we see in the model cases even in the future applications262.2. Flux Error Reductionof Euler and Navier-Stokes equations which are not necessarily steady state problems.This and other time evolution algorithms that are used within ANSLib have been verywell documented by Nejat [33] and Michalak [30, 31].2.2 Flux Error ReductionFew other researchers have attempted to address similar motivations as those described inSection 1.1.Katz and Sankaran [22] proposed a novel approach to obtain third-order accurate so-lution to inviscid problems by means of a flux correction term. Unlike our motivation,their ultimate goal was to obtain higher accuracy for the solution rather than the gradi-ents. Their flux correction term is aimed at specifically improving the truncation errorthrough error-cancellations which they show translates to improvements in the accuracy ofthe solution. It is also unique since it can be directly added to any existing second-orderaccurate vertex-centred (also known as node-centred) scheme. However, as they mention,this flux correction term is unfortunately limited in its applications to vertex-centred linearGalerkin schemes. Moreover, the computation of this term is itself not trivial and requiresat least a second-order accurate estimate of the gradient of the flux. This flux gradientis computed using a quadratic least squares approximation. The additional computationaltime required per iteration for this corrected scheme was reported to be around 50% morethan the regular vertex-centred linear Galerkin schemes.This flux correction scheme was also extended to unsteady viscous flows by Pincock andKatz [44]. As can be anticipated, the correction terms required for the viscous fluxes aresignificantly more computationally intensive, with a 40-70% increase over a second-orderfinite-volume method. Unfortunately, comparisons to a third-order finite-volume methodwere not shown, which makes it tough to compare our schemes. Here, for the viscousterms, they propose a new cubic polynomial “Jacobian-averaged” finite-element mappingto compute the gradients within each control-volume, which contributes to the significantcosts. They report some instabilities with regards to certain problems and also express theirwork in progress to enable the use of larger time steps.272.2. Flux Error ReductionAnother equally interesting approach was developed by Nishikawa for the diffusion equa-tion [34] and later extended to encompass advection as well [36]. The key concept he pro-poses is the use of a first-order hyperbolic system that is equivalent to the the regularsecond-order diffusion equation, albeit only in the steady state. A secondary yet remark-able feature is that his new scheme gives second-order accuracy for both the solution aswell as the gradients. Although developed for a second-order residual-distribution scheme,he suggests that it could be developed in a similar manner for finite-difference and finite-volume schemes as well. At its current state of development it only works with an explicittime advancing strategy and an implicit time advance would ruin the solution accuracy.However, they do report an advantage of being able to use much larger explicit time stepsthan other regular Galerkin schemes.28Chapter 3H1 Finite-Volume SchemeIn this chapter we present the theory behind our proposed H1 finite-volume scheme. Wewill describe our scheme in one dimension before examining the two-dimensional scheme indetail.Our overall objective can be broadly put forth as striking a balance between the bestfeatures of third-order accurate least-squares with nearly the speed and simplicity of second-order methods. An overview of the scheme that we propose is as follows1. Compute a reconstruction that is second order for both the solution and the gradients.2. Utilize this reconstruction to evaluate a second-order diffusive flux at the Gauss points.3. Use single point Gauss-quadrature (located on the face mid-side) to obtain the fluxintegral for each control volume.4. Update the solution states in every control volume.The first step is in essence to obtain the gradients to second-order accuracy at the control-volume reference points with the least computational effort while the second step is a newflux discretization we propose. Our modifications do not apply to the last two steps, theregular flux-integral evaluation and time-stepping strategies, which were discussed in briefin section 2.1.3.1 H1 ReconstructionThe first step in our new scheme is to obtain the gradients to second-order accuracy at thecontrol volume reference points. A regular second-order reconstruction scheme computesjust the solution and its first-derivative (gradient) coefficients and as such requires just three293.1. H1 Reconstructionface neighbours from the first layer in Figure 2.2 to complete its reconstruction stencil. Itthus boasts a relatively cheap pseudo-inverse that has the general form ofPInv2ndLS =A11 A12 A13A21 A22 A23A31 A32 A33(3.1)In contrast a third-order reconstruction typically uses nine neighbours to enable theestimation of the three additional second-derivative (quadratic) coefficients. Hence its least-squares system and pseudo-inverse are more complex and computationally intensive. Thisresults in a third-order accurate approximation to the solution and a second-order accurateapproximation to the solution gradient.PInv3rdLS =A11 A12 A13 A14 A15 · · · A19A21 A22 A23 A24 A25 · · · A29A31 A32 A33 A34 A35 · · · A39A41 A42 A43 A44 A45 · · · A49A51 A52 A53 A54 A55 · · · A59A61 A62 A63 A64 A65 · · · A69(3.2)With H1 reconstruction we do not need to compute the three additional second-derivativecoefficients as will be evidenced in the next subsection. Thus this enables us to eliminatethe last three rows corresponding to the quadratic coefficients. This saves us considerableresources and computational time when repeatedly performing the matrix-vector multipli-cation between the pseudo-inverse and the control-volume averages to obtain the solutionreconstruction coefficients at each time step. It also provides a substantial gain in efficiency303.2. H1 Flux Discretizationi+1/2i-1 i i+1 i+2hi-1 hi hi+1 hi+2Figure 3.1: 1D Non-uniform Grid Stencilsince we get a relatively large jump in accuracy without the laborious computation that thequadratic terms would demand. Note that the three pseudo-inverse rows that we are leftwith are exactly the same as those from a third-order reconstruction.PInvH1 =A11 A12 A13 A14 A15 · · · A19A21 A22 A23 A24 A25 · · · A29A31 A32 A33 A34 A35 · · · A39(3.3)3.2 H1 Flux DiscretizationOnce we have obtained the necessary second-order gradients at the control-volume referencepoints we must translate those to the Gauss points while maintaining our accuracy. Herewe discuss our new approach to the flux discretization at the Gauss point which countersthe loss of error cancellations due to the irregularity of the mesh.3.2.1 Theory in One DimensionThe quickest way to test the waters with a proof of concept is to take a look at what can beaccomplished on a one-dimensional grid. In our case that turns out to be relatively straight-forward. We will expand on the same with reference to the one-dimensional non-uniformgrid stencil shown in Figure 3.1; this is the one-dimensional analog to an unstructured mesh.Consider the Taylor series expansion of the first derivative about the cell centres in terms313.2. H1 Flux Discretizationof the derivative at the cell interface. Note that the reconstructed first derivatives(d˜φdx)are a second-order approximation to the exact derivatives(dφdx)at the cell reference point.d˜φdx∣∣∣∣∣i+O(h2) =dφdx∣∣∣∣i=dφdx∣∣∣∣i+1/2− hi2d2φdx2∣∣∣∣i+1/2+O(h2) (3.4)d˜φdx∣∣∣∣∣i+1+O(h2) =dφdx∣∣∣∣i+1=dφdx∣∣∣∣i+1/2+hi+12d2φdx2∣∣∣∣i+1/2+O(h2) (3.5)Now normally for a least-squares reconstruction we would take a direct average of theabove two Equations 3.4 and 3.5, to obtain the derivative at the cell interfacedφdx∣∣∣∣i+1/2=12{d˜φdx∣∣∣∣∣i+d˜φdx∣∣∣∣∣i+1}−(hi+1 − hi4)d2φdx2∣∣∣∣i+1/2+O(h2) (3.6)For structured (uniform) meshes the second term on the right hand side of Equation 3.6would cancel out upon taking an arithmetic average; unfortunately the same cannot be saidof unstructured meshes. Thus structured meshes enjoy a second-order accurate derivativewhile unstructured mesh counterparts have to contend with first-order derivatives.Instead of using a simple arithmetic average, if we were to eliminate each of the first-order terms on the right hand side of Equations 3.4 and 3.5 with a linear interpolation, thenwe can expect a gain in accuracy even for unstructured (non-uniform) grids.dφdx∣∣∣∣i+1/2=1(hi+1 + hi){hi+1.d˜φdx∣∣∣∣∣i+ hi.d˜φdx∣∣∣∣∣i+1}+O(h2) (3.7)This enables us to gain back second-order accuracy even in unstructured meshes withoutthe use of the higher derivatives. Thus the higher derivatives obtained from the reconstruc-tion are made redundant and need not be computed. In the next subsection we look atextending this technique to two dimensions.3.2.2 Extension to Two DimensionsIn two dimensions we essentially only require the component of the gradient in the nor-mal direction to compute the diffusive flux at any given control-volume interface. The323.2. H1 Flux DiscretizationCircumcenterContainment centerCentroidCircumcircleContainment circleFigure 3.2: Various control-volume reference centrestraditional use of centroids as the control-volume reference point has proven to be a stableand robust choice. However, obtaining the second-order normal gradient components usingcentroids necessitates the use of the second derivatives. Another control-volume referencepoint candidate that has been rarely looked into is the circumcentre. A circumcentre isdefined as the point where the perpendicular bisectors of the faces of a triangle intersect.Obtaining the second-order normal gradient components using them should be easier, sincethey always lie on the face normal through the edge mid-side. Unfortunately the use ofcircumcentres as control volume reference points leads to instabilities, which will be furtherdiscussed in section 6.1. The use of containment centres instead provides better stabilitywithout hampering the accuracy gains. Containment centres are defined as the centres ofthe smallest circle that can enclose a given triangle. The trio of control-volume referencepoints — namely the centroid, the circumcentre and the containment-centre — are illus-trated in Figure 3.2. Our flux discretization strategy using a combinations of circumcentresor containment centres is presented below.3.2.2.1 CircumcentresCircumcentres are the prime candidate for control-volume reference points for this scheme,since they always lie on the perpendicular bisector of the faces. However with unstructuredmeshes the triangles are not guaranteed to be acute. Hence there will definitely arise sce-narios where the circumcentre falls outside the cell or worse coincides with the circumcentre333.2. H1 Flux Discretization× Circumcentres • Gauss-quadrature pointshRhLLeft Cell Right CellGLRητ(a) The Good circumcentre casehRhLLeft CellRight CellGLRητ(b) The Bad circumcentre casehRhLLeft CellRight CellGLRητ(c) The Ugly circumcentre caseFigure 3.3: Circumcentre casesof a neighbour. These situations fall into the following three categories. Treating hL andhR as signed normal distance to the Gauss point, each positive into its own cell, enablesthe identification and handling of these scenarios.The Good: This is the best case scenario — and the most common for isotropic meshes— where the cells on either side of the face are acute and hence the circumcentres arewell behaved as illustrated by Figure 3.3a. The circumcentres lie on either side of theGauss point, so we can directly implement a linear interpolation formula similar tothat in one dimension given by Equation 3.7. We note that the same formula can alsobe used in cases where both the triangles are obtuse in which case the signed distances,being both negative, would just negate each other. The second-order gradient at the343.2. H1 Flux DiscretizationGauss point is given by∂φ∂η∣∣∣∣g=1(hL + hR){hR∂˜φ∂η∣∣∣∣∣L+ hL∂˜φ∂η∣∣∣∣∣R}+O(h2) (3.8)The Bad: Here we face an obtuse triangle on one side of the face as illustrated by Fig-ure 3.3b. Thus both the circumcentres are on the same side of the Gauss point andone of those signed distances is negative. We still obtain a second-order gradient atthe Gauss point using Equation 3.8 thanks to the linear extrapolation formula formedimplicitly by the use of signed distances.The Ugly: This is the last and worst scenario, which occurs when the two circumcen-tres coincide. In such cases we revert to arithmetic averaging of cell gradients, asinterpolation is impossible.∂φ∂η∣∣∣∣g=12{∂˜φ∂η∣∣∣∣∣L+∂˜φ∂η∣∣∣∣∣R}+O(h) (3.9)However, these cases were later identified as the root cause of instabilities to beexplained in Section 6.1, and hence were eliminated with the use of containmentcentres.3.2.2.2 Containment CentresContainment centres are defined as the centres of the smallest circle that can enclose agiven triangle. They are identical to the circumcentre of an acute triangle but have theadded advantage of not spilling out of the triangle when it is obtuse. In that scenario thecontainment centre is the mid-side of the longest edge of the triangle, which enables us todirectly utilize its value for our face gradient. This should theoretically be ideal for thatspecific Gauss point but hurt accuracy for the flux at the remaining Gauss points on theother faces.∂φ∂η∣∣∣∣g={∂˜φ∂η∣∣∣∣∣R}+O(h2) (3.10)Two approaches to the implementation of containment centres can be thought of. First353.2. H1 Flux Discretization(hR= 0)hL- Cell containment-center- Gauss quadrature pointLeft CellRight CellGLRητFigure 3.4: The obtuse containment centrewould be a blanket implementation where containment centres are implemented through-out the mesh thus eliminating the bad and ugly cases. The other approach would usecontainment centres to counter the ugly circumcentre cases, which were the root cause ofthe stability problem. The latter proved much more accurate and hence results presentedin this thesis use this approach.3.2.2.3 High-Accuracy Jump TermAs discussed in subsection 2.1.3.3, the jump term that accounts for the difference in thereconstructed solution in the two cells at the Gauss point plays a significant role in enhancingthe stability of a diffusion scheme and making it more robust. Hence we wished to takeadvantage of its damping effects. However a regular jump term would adversely affectthe accuracy of our discretization since it would only be first-order accurate, given thatthe solution is second-order accurate. Hence we have formulated a new jump term withhigher-order accuracy that would enhance stability while serving our accuracy needs.The first requirement to obtain the high-accuracy jump term is an estimate of thesecond-derivative at the Gauss point in the normal direction. We obtain this estimate usinga central-difference formulation between the second-order gradients at the control-volume363.2. H1 Flux Discretizationreference points.∂2φ∂η2∣∣∣∣L≈∂2φ∂η2∣∣∣∣mid={∂˜φ∂η∣∣∣L− ∂˜φ∂η∣∣∣R}hLR + +O(h2)We then use this to estimate more accurate left and right solution states at the Gausspoint of the face.φ˜g+ ≈(φ˜L +∂φ∂x∣∣∣∣LhLx +∂φ∂y∣∣∣∣LhLy)+∂2φ∂η2∣∣∣∣Lh2L2where the term in brackets on the right-hand side is evaluated implicitly from the gradientreconstruction.We can then use these solution states in the jump term formulation to enhance itsaccuracy.(∇φ)jump =αhEdgeLength(φ˜g+ − φ˜g−)The constant , with a value of 10−9, as well as the face edge-length being used asthe length base, are precautions in place to avoid scenarios where the circumcentres orcontainment centres coincide and we would be dividing by zero.37Chapter 4Analytical Tests and AnalysisAs we discussed in the previous chapter, a p-order-accurate solution reconstruction willproduce (p-1)-order-accurate gradients which in turn leads to (p-1)-order diffusive flux eval-uation. Ollivier-Gooch and Van Altena [43] were able to explain this loss in accuracy byanalytically solving the least-squares system on uniform triangular meshes. They analyzedthe flux integral, and its associated truncation error, for the Laplacian operator obtained us-ing linear (second-order), quadratic (third-order) and cubic (fourth-order) reconstructions.The Laplace equation being a simplified form of the Poisson equation without a sourceterm, can be formulated as:∇2φ = 0. (4.1)Given the uniform and symmetric nature of the mesh, computations become simplifiedenough to be handled by symbolic mathematical software tools such as Maple, but are stilltoo laborious to be done by hand.Here we extend the analysis presented by [43], to our proposed H1 reconstruction schemeand compare it with the second and third-order least-square schemes.4.1 MethodologyTo begin, a coordinate and labelling system is created for a uniform, equilateral triangularmesh as illustrated by Figure 4.1. Each cell-centred control-volume triangle has edges oflength h. The index coordinates (i, j), used for identifying the cells, do not run parallel tothe Cartesian coordinates (x, y), with the j axis running from the bottom left to top rightcorner. Each parallelogram with sides parallel to the i and j axes is labelled with a unique(i, j) pair and the respective triangles are distinguished as being either up (∆) or down (∇)384.2. Analytical ResultsPQ QQRRRRR RSSSS SSTTTUWUVUWU V UWUVFigure 4.1: Analytical Mesh coordinates and cell labelspointing. Triangles have also been classified to belong to a particular family based on theircoefficients for the flux integral. The control volumes from the same family share the samecoefficients due to the symmetry of the mesh. We will perform our reconstruction aboutthe central P cell.The control volume moments are computed by analytical integration. We can thenprepare the finite-volume framework that was described in chapter 2 to analytically solvefor the solution, diffusive flux and the flux integral. The next section details our results andanalysis.4.2 Analytical ResultsThe results of our analysis have been summarized in Table 4.1. Further discussion on eachof the cases follow.4.2.1 Linear ReconstructionBeginning with the linear (second-order) reconstruction, as already discussed the first threeneighbours (Q’s) are sufficient to satisfy the least-squares constraint requirements. All394.2. Analytical ResultsReconstruction typeLinear Regular H1 Extended H1 QuadraticRecon. stencil(weight) Q(3/h2) Q (3/h2),R(1/h2) Q (3/h2),R (1/h2),T(3/4h2) Q (3/h2),R (1/h2)Leading T.E. term h216∇2∇2φh√312 β (φ)+h28 ∇2∇2φh√312 β (φ)+73h2564 ∇2∇2φh√345 β (φ)+7h2144∇2∇2φLaplacianstencil multiplier23h213h21267h249h2Family P (1) -6 -18 -504 -15Family Q (3) 0 -2 -56 3Family R (6) 1 3 69 1Family S (6) - 2 36 -Family T (3) - -2 -16 -Family W (3) - - 30 -β (φ) ≡ 3 ∂3φ∂x2∂y− ∂3φ∂y3Table 4.1: Stencils and truncation error for Laplacian on a cell-centred mesh404.2. Analytical Resultsthe geometric weights were set to inverse-distance-squared (distance between cell-referencelocations). The second-order least-squares solution for an arbitrary control volume (i,j,∆)is φ∂φ∂x∂φ∂y=φ(i,j,∆)φ(i,j,∇)−φ(i−1,j,∇)hφ(i,j,∇)+φ(i−1,j,∇)−2φ(i,j−1,∇)h√3. (4.2)An analogous result is obtained for down-pointing triangles. Proceeding to compute thediffusive flux and the single-Gauss-point flux integral as discussed in Section 2.1, we obtainthe following expression in terms of control-volume averages for the computed Laplacian inthe control volume (0, 0,∆)∇˜2φ0,0 =(23h2)(−6φP + ΣφR) (4.3)where ΣφR is the sum of the control-volume averaged solution for all the cells labelled R(six in total) in Figure 4.1.We can observe that this expression only involves the up-pointing triangle families, Pand R, and as such the solution in the down-pointing triangles are “decoupled” from itsup-pointing neighbours. This is analogous to the odd-even decoupling commonly observedin structured meshes. The result of this decoupling can be visualized in Figure 4.2, whenwe run our numerical solver ANSLib (described in the next chapter) on this uniform equi-lateral triangle mesh for the Poisson test case. The decoupling clearly affects the solutionconvergence since the neighbours are not “communicating” with each other. In particularthe down-pointing triangles do not have access to the boundary conditions without the helpof the up-pointing triangles.The remedy to this is to break the symmetry that brings about the decoupling or tointroduce some form of external coupling. The former is solved by utilizing unstructuredmeshes while the latter can be achieved with the help of a jump term described in Sec-tion 2.1.3.3. The use of the scheme on unstructured meshes still admits the possibility ofsome decoupling effects, but presence of a jump term should completely negate decoupling414.2. Analytical ResultsFigure 4.2: Second-order least-squares on uniform equilateral triangle mesh (decoupling)effects since it directly incorporates the two neighbouring solutions in its formulation. Thesolution thus obtained is essentially the same as Figure 4.3.We can expand the control-volume averages in Equation 4.3 in terms of derivatives of theunderlying smooth solution φ at the origin by using Equation 2.17. This gives us the errorin the computed Laplacian ∇˜2φ0,0 from the control-volume average of the exact Laplacian∇2φ0,0 as∇˜2φ0,0 −∇2φ0,0 =h216(∂4φ∂x4+ 2∂4φ∂x2∂y2+∂4φ∂y4)+O(h4). (4.4)The presence of an order-h2 (O(h2)) term on the right-hand side shows that this schemefor calculating the Laplacian is second-order accurate on structured cell-centred meshes.Moreover, by substituting the control-volume average expansions into the gradient for-mulations in Equation 4.2 it can be easily shown that the first derivatives (gradients) areonly first-order accurate. This supports our numerical findings presented in Section 1.1,which was the original motivation for this work.424.2. Analytical Results4.2.2 Quadratic ReconstructionProceeding similarly for the (third-order) quadratic reconstruction, all second-layer neigh-bours (R’s) are required to make the least-squares system over-constrained. We applyinverse-distance-squared weighting to these new neighbours. A look at the solution of thethird-order least-squares problem in an arbitrary control volume (i,j,∆), shows the need forthe additional layer of neighbours.φ∂φ∂x∂φ∂y∂2φ∂x2∂2φ∂x∂y∂2φ∂y2=φi,j,∆φi+1,j,∆−φi−1,j,∆−φi,j−1,∆+φi+1,j−1,∆+3φi,j,∇−3φi−1,j,∇6h√3(2φi,j+1,∆+2φi−1,j+1,∆−φi,j−1,∆−φi+1,j−1,∆+3φi,j,∇+3φi−1,j,∇−φi+1,j,∆−φi−1,j,∆−6φi,j−1,∇)18h−30φi,j,∆+7φi+1,j,∆+7φi−1,j,∆+9φi,j,∇+9φi−1,j,∇−2φi,j+1,∆−2φi−1,j+1,∆+φi,j−1,∆+φi+1,j−1,∆9h2√3(2φi,j+1,∆−2φi−1,j+1,∆+3φi,j−1,∆−3φi+1,j−1,∆+3φi,j,∇−3φi−1,j,∇−φi+1,j,∆+φi−1,j,∆)9h2φi,j−1,∆−φi+1,j,∆−φi−1,j,∆−10φi,j,∆+2φi,j+1,∆+2φi−1,j+1,∆+φi+1,j−1,∆+4φi,j−1,∇+φi,j,∇+φi−1,j,∇3h2(4.5)Once again we proceed to find the diffusive fluxes, which now incorporate the second-derivative terms, and the flux integral using a two-point Gauss quadrature. This gives usthe following expression in terms of control-volume averages for the computed Laplacian inthe control volume (0, 0,∆)∇˜2φ0,0 =(49h2)(−15φP + 3 ΣφQ + ΣφR) (4.6)where, for example, ΣφQ is the sum of the control-volume averaged solution for all the cellslabelled Q (three in total) in Figure 4.1.Equation 4.6 incorporates both the up as well as the down-pointing triangles, unlikethe second-order version seen in Equation 4.3. Thus we can expect this third-order least-squares scheme to not be decoupled. This is verified by the successful convergence of ournumerical simulation of this reconstruction scheme on an equilateral triangle mesh, shown434.2. Analytical ResultsFigure 4.3: Third-order least-squares on uniform equilateral triangle meshin Figure 4.3.Upon substituting the respective expansions of the control-volume averages from Equa-tion 2.17 into Equation 4.6, we can find its exact error and its order of accuracy. The errorin the Laplacian for third-order least-squares is given by∇˜2φ0,0−∇2φ0,0 =h√345(3∂3φ∂x2∂y− ∂3φ∂y3)+7h2144(∂4φ∂x4+ 2∂4φ∂x2∂y2+∂4φ∂y4)+O(h3). (4.7)This is a surprising result, since the Laplacian is only first-order accurate here due tothe presence of an O(h) term on the right-hand side. Moreover, even if we were to ignorethis first-order term, we still have to contend with the second-order accurate term which isonly slightly smaller than its counter-part in the linear reconstruction. As explained in [43]this first-order term is an unfortunate combination of third-derivative terms that cannot becomputed independently for any single family of control volumes. This can be primarilyattributed to the reconstruction’s symmetric stencils being unable to distinguish between thenon-zero x2y and y3 moments in the Taylor series expansion of the control-volume averageson uniform triangular meshes. However, fortunately this does not extend to unstructured444.2. Analytical Resultsmeshes where, as we have seen in Figure 1.3b, the third-order (quadratic) least-squaresreconstruction performs much better, thanks to the absence of this fatal symmetry.4.2.3 H1 ReconstructionFor our proposed H1 reconstruction we present two cases, one with a regular quadraticstencil (Q +R families) and the other with an extended stencil which encompasses all thevertex neighbours of the given control volume (Q+R+T families). This provides us withinsights into how our scheme is affected by extending the reconstruction stencil.4.2.3.1 Regular Stencil H1 ReconstructionThe reconstruction stencil used for our proposed H1 reconstruction is the same as the oneused for the quadratic reconstruction. Thus the reconstruction solution we obtain is alsosimilar save that we do not obtain the second-degree derivatives, as given by:φ∂φ∂x∂φ∂y=φi,j,∆φi+1,j,∆−φi−1,j,∆−φi,j−1,∆+φi+1,j−1,∆+3φi,j,∇−3φi−1,j,∇6h2(φi,j+1,∆+φi−1,j+1,∆)−φi,j−1,∆−φi+1,j−1,∆+3(φi,j,∇+φi−1,j,∇)−φi+1,j,∆−φi−1,j,∆−6φi,j−1,∇6h√3(4.8)Using a two-point Gauss quadrature, we obtain the following expression for the com-puted Laplacian in the control volume (0, 0,∆)∇˜2φ0,0 =(13h2)(−18φP − 2 ΣφQ + 3ΣφR + 2 ΣφS − 2 ΣφT ) (4.9)in terms of the control-volume averages and∇˜2φ0,0−∇2φ0,0 =h√312(3∂3φ∂x2∂y− ∂3φ∂y3)+h28(∂4φ∂x4+ 2∂4φ∂x2∂y2+∂4φ∂y4)+O(h3) (4.10)as the error in the computed Laplacian.454.2. Analytical ResultsFigure 4.4: Regular H1 reconstruction on uniform equilateral triangle mesh (decoupling)We can instantly recognize the similarities with the parent quadratic reconstruction,from which H1 reconstruction is derived, especially with a similar first-order term in the errorexpression for Laplacian, albeit larger by almost a constant factor of 4. The second-orderterm seems larger than even its counter-part in linear reconstruction. However, fortunatelythe scheme fares better when used with realistic unstructured meshes, as we will see inChapter 5.Another anomaly we observed was an unforeseen pattern of solution decoupling, as seenin Figure 4.4, when numerical simulations were conducted on the uniform equilateral mesh.However this is not supported by our understanding of Equation 4.9, which shows no signsof decoupling as was evident in Equation 4.3 for the linear reconstruction scheme. Onceagain, this decoupling can be countered with the aid of a jump term, just like in the linearreconstruction case.464.2. Analytical Results4.2.3.2 Extended Stencil H1 ReconstructionFor the extended stencil H1 Reconstruction, we include all the vertex neighbours of thegiven control volume, namely the Q, R and T family of cells seen in Figure 4.1. We includejust three additional neighbours (T family) instead of the entire third-layer of neighbourswhich includes the S family of neighbours as well. The reconstruction solution we get isgiven by:φ∂φ∂x∂φ∂y=φi,j,∆45(φ(i+1,j−1,∇) − φ(i−1,j−1,∇))+ 44(φ(i+1,j−1,∆) − φ(i,j−1,∆))+ 64(φ(i+1,j,∆) − φ(i−1,j,∆))+252(φ(i,j,∇) − φ(i−1,j,∇))+ 20(φ(i,j+1,∆) − φ(i−1,j+1,∆))534h183(φ(i−1,j,∇) + φ(i,j,∇))− 15(φ(i+1,j−1,∇) + φ(i−1,j−1,∇))− 41(φ(i+1,j,∆) + φ(i−1,j,∆))−61(φ(i,j−1,∆) + φ(i,j−1,∆))+ 102(φ(i−1,j+1,∆) + φ(i,j+1,∆))+ 30φ(i−1,j+1,∇) − 366φ(i,j−1,∇)534h√3(4.11)Once again using a two-point Gauss quadrature, we obtain the following expression forthe computed Laplacian in the control volume (0, 0,∆)∇˜2φ0,0 =(1267h2)(−504φP − 56 ΣφQ + 69ΣφR + 36 ΣφS − 16 ΣφT + 30 ΣφU) (4.12)in terms of the control-volume averages and∇˜2φ0,0 −∇2φ0,0 =h√312(3∂3φ∂x2∂y− ∂3φ∂y3)+287h22136(∂4φ∂x4+ 2∂4φ∂x2∂y2+∂4φ∂y4)+O(h3)(4.13)as the error in the computed Laplacian.Here too we find the same decoupling quirk as seen with the regular case, as observedin Figure 4.5.It is apparent that an extension of the stencil does not have much impact on the accuracyof the system, however we have observed a slight increase in stability of the system withthe extension of the stencil during numerical simulation using our flow-solver ANSLib.However the additional stability was negligible compared to the effects of other measuressuch as the use of a jump term. Moreover the added benefits do not justify the more474.2. Analytical ResultsFigure 4.5: Extended H1 reconstruction on uniform equilateral triangle mesh (decoupling)intensive computations required due to the addition of these additional vertex neighbours(and consequently, additional rows in the left-hand side of every control-volume’s leastsquares system).Henceforth, we will stick to the regular H1 reconstruction for all of our numerical testsand simulations that follow in the next chapter.48Chapter 5Numerical Tests and AnalysisThe analytic tests in the previous chapter were designed for meshes with regular topologyto provide insights into the behaviour of the discretization schemes on a structured mesh.Similar analysis on an unstructured mesh would prove highly impractical for every controlvolume would have its own unique qualities and reconstruction stencil. Here instead we com-pare and analyze the various accuracy results from our proposed H1 finite-volume schemeand the conventional second and third order least-squares reconstruction schemes. All thenumerical experiments were conducted using the Advanced Numerical Simulations Library(ANSLib), the unstructured numerical simulation code developed in Advanced NumericalSimulations Laboratory (ANSLab), UBC.The precise and consistent evaluation of the numerical errors generated is essentialto analyze the accuracy of any given scheme. The discretization error, which forms themajority of the numerical error, along with its primary source, the truncation error, canbe estimated by a technique known as the method of manufactured solutions (MMS). Thistechnique has not only been implemented to estimate the numerical error but also used fororder of convergence studies and code verification [54].In MMS, an exact solution that satisfies the given governing equations is known (ormanipulated using a source term) and preset. Hence errors can be computed by evalu-ating and comparing the final discrete solution obtained through the simulation with thisexact solution. The same is done for the corresponding flux and flux integrals. The L2-norm of these errors are then evaluated for a series of five consistently-refined cell-centredunstructured meshes and analyzed to obtain the asymptotic order of convergence of thescheme. The mesh refinement was performed by reducing the length scale by a factor oftwo to form the next finer mesh. The triangular meshes are regular isotropic and non-nested495.1. Preliminary Studies(a) Square mesh with 780 triangles (b) Square mesh with 3132 triangles(c) Channel mesh with 2162 trianglesFigure 5.1: Isotropic unstructured meshes usedmeshes and were generated by Delaunay refinement. Figure 5.1 shows a few of the isotropicunstructured meshes used.We begin with some preliminary studies, which verify the foundations of our new H1finite-volume scheme. Complete accuracy analysis was then performed for two separatemodel problems, namely the Poisson equation and the advection-diffusion equation.5.1 Preliminary Studies5.1.1 One-dimensional H1 Flux DiscretizationPreliminary studies were undertaken to verify the advantages of H1 flux discretization in one-dimension that was described in subsection 3.2.1. A Matlab code was written to comparethe accuracy achieved with the use of the simple arithmetic averaging as well as the linearinterpolation given by Equations 3.6 and 3.7 respectively, assuming second-order gradientsat the cell reference points. The tests were conducted on a series of randomly perturbed505.1. Preliminary Studies1e+01 1e+02 1e+031e−061e−051e−041e−031e−021e−01Number of Control VolumesL 2 error in Flux 1st order2nd orderL2 GradientH1 Gradient(a) One-dimensional H1 flux discretization-6-5-4-3-2-102 3 4 5log 10( L 2 error )log10( Number of Control Volumes )1st order2nd orderH1 GradXH1 GradYH1 SolnL2 GradXL2 GradYL2 Soln(b) Order of reconstruction accuracy in two dimensionsFigure 5.2: Preliminary Studiesone-dimensional grids with the cell mid-side as the cell reference point. A test function ofe−x2 was used as the manufactured solution with Dirichlet boundary conditions.Figure 5.2a shows the results of the preliminary tests in one-dimension, which verifiesour theory. The linear interpolation of H1 reconstruction does indeed provide us withsecond-order gradients at the cell interface in contrast to the first-order gradients obtainedwith simple arithmetic averaging of second-order least-squares scheme.5.1.2 Two-dimensional H1 Reconstruction TestsThe consistency of the H1 reconstruction in two-dimensions, detailed in section 3.1, wasalso verified. A manufactured solution test function of pi8 sin (pix) sin (piy) was used to as-sess the accuracy of the reconstruction, with homogeneous Dirichlet boundary conditionenforcement.Figure 5.2b compares the reconstruction order of accuracy of our new H1 reconstructionwith that of second-order least-squares reconstruction for the solution and the two gradientcomponents.As is clear, we have been successful in increasing the order of accuracy of the recon-structed gradients to second order using H1 reconstruction while the same are first-order515.2. 2D Poisson Test(a) Case (i) Sinusoidal:pi8sin (pix) sin (piy)(b) Case (ii) Gaussian Bump:e−(x2+y2)/0.2(c) Case (iii) Plain Laplacian:cos(pix) sinh(piy)sinh(pi)Figure 5.3: Manufactured solutions for Poisson equation on an unstructured meshaccurate in second-order (linear) least-squares reconstruction.5.2 2D Poisson TestThe Poisson equation serves as the ideal scalar model of diffusion problems. Three uniquemanufactured solution test functions were used to assess the accuracy of the reconstruction,namely:Case (i) Sinusoidal: pi8 sin (pix) sin (piy), with Dirichlet boundary conditionsCase (ii) Gaussian Bump: exp(− (x2 + y2) /0.2), with Dirichlet boundary conditionsCase (iii) Plain Laplacian: cos (pix) sinh (piy) / sinh (pi), with Dirichlet boundaries on thetop and bottom walls and Neumann boundaries on the side wallsFigure 5.3 shows the resultant manufactured solutions on a unit square mesh consisting of780 control volumes, with the origin located at the bottom left corner.The orders of accuracy of the Poisson solution, flux and flux integral have been plottedand compared in Figures 5.4, 5.5 and 5.6. The tests were run on 5 successively refinedmeshes containing 176, 780, 3132, 12766 and 50715 isotropic triangles.The diffusive flux at the Gauss points using the smooth initial data is second-orderaccurate as seen in Figures 5.4a, 5.5a and 5.6a.525.2. 2D Poisson TestQuantity Second-orderleast squaresH1reconstructionThird-orderleast squaresCase(i)Case(ii)Case(iii)Case(i)Case(ii)Case(iii)Case(i)Case(ii)Case(iii)Flux(initial data)1.14 1.07 1.21 1.82 2.01 1.73 2.06 2.10 2.22Flux integral(initial data)0.01 0.11 0.14 0.62 0.56 0.47 1.36 1.20 1.10Flux(converged)1.22 1.38 1.41 1.53 1.58 1.32 2.05 2.11 2.20Solution 2.17 2.03 2.31 2.02 2.08 2.05 2.05 2.01 3.00Table 5.1: Orders of accuracy for Poisson model problemWith the flux-integral order obtained at this stage, as seen in Figures 5.4b, 5.5b and 5.6b,we can observe that H1 reconstruction is indeed first order accurate but loses accuracy onthe finest mesh. None the less it still shows much better accuracy compared to a second-order least-squares method. Note that this flux integral was still computed by using just asingle Gauss quadrature point.When the solution is fully converged, the flux shows a trend similar to the flux integral,losing its order as it progresses to finer meshes; this is illustrated by Figures 5.4c, 5.5c and 5.6c.As for the converged solution itself, seen in Figures 5.4d, 5.5d and 5.6d, all three schemesare comparable and essentially produce second-order accurate solutions.The orders of accuracy achieved for the solution, diffusive flux and the flux integral witheach scheme and for all cases, are tabulated in Table 5.1. The values listed are the averagedrates of convergence across all the five refinement levels.Another interesting feature of the H1 reconstruction can be found when we consider itssolution error. When compared to the other least-squares schemes, seen in Figure 5.7 for thetest case (i), it is able to overcome the rough noise plaguing the second-order least-squaresdiscretization error and is closer to a third-order least-squares in terms of its smoothness.Overall, the differences in the error distribution for the three schemes are quite striking,considering that the error norms are nearly identical.535.2. 2D Poisson Test-6-5-4-3-2-12 3 4 5log 10( L 2 error in Flux )log10( Number of Control Volumes )1st order2nd orderH1 recon2nd Least-squares3rd Least-squares(a) Flux convergence rate at Gauss points (initialdata)-3-2-102 3 4 5log 10( L 2 error in Flux Integral )log10( Number of Control Volumes )1st orderH1 recon2nd Least-squares3rd Least-squares(b) Flux integral convergence rate (initial data)-6-5-4-3-2-12 3 4 5log 10( L 2 error in Flux )log10( Number of Control Volumes )1st order2nd orderH1 recon2nd Least-squares3rd Least-squares(c) Flux convergence rate at Gauss points (convergedsolution)-6-5-4-3-22 3 4 5log 10( L 2 error in Solution )log10( Number of Control Volumes )2nd orderH1 recon2nd Least-squares3rd Least-squares(d) Solution convergence rateFigure 5.4: Case (i) Sinusoidal Poisson 2D test results545.2. 2D Poisson Test-6-5-4-3-2-12 3 4 5log 10( L 2 error in Flux )log10( Number of Control Volumes )1st order2nd orderH1 recon2nd Least-squares3rd Least-squares(a) Flux convergence rate at Gauss points (initialdata)-3-2-1012 3 4 5log 10( L 2 error in Flux Integral )log10( Number of Control Volumes )1st orderH1 recon2nd Least-squares3rd Least-squares(b) Flux integral convergence rate (initial data)-6-5-4-3-2-12 3 4 5log 10( L 2 error in Flux )log10( Number of Control Volumes )1st order2nd orderH1 recon2nd Least-squares3rd Least-squares(c) Flux convergence rate at Gauss points (convergedsolution)-6-5-4-3-22 3 4 5log 10( L 2 error in Solution )log10( Number of Control Volumes )2nd orderH1 recon2nd Least-squares3rd Least-squares(d) Solution convergence rateFigure 5.5: Case (ii) Gaussian Bump Poisson 2D test results555.2. 2D Poisson Test-6-5-4-3-2-12 3 4 5log 10( L 2 error in Flux )log10( Number of Control Volumes )1st order2nd orderH1 recon2nd Least-squares3rd Least-squares(a) Flux convergence rate at Gauss points (initialdata)-2-1012 3 4 5log 10( L 2 error in Flux Integral )log10( Number of Control Volumes )1st orderH1 recon2nd Least-squares3rd Least-squares(b) Flux integral convergence rate (initial data)-6-5-4-3-2-12 3 4 5log 10( L 2 error in Flux )log10( Number of Control Volumes )1st order2nd orderH1 recon2nd Least-squares3rd Least-squares(c) Flux convergence rate at Gauss points (convergedsolution)-7-6-5-4-3-22 3 4 5log 10( L 2 error in Solution )log10( Number of Control Volumes )2nd orderH1 recon2nd Least-squares3rd Least-squares(d) Solution convergence rateFigure 5.6: Case (iii) Plain Laplacian Poisson 2D test results565.3. 2D Advection-Diffusion Test(a) Second-order least-squares (b) H1 reconstruction (c) Third order least-squaresFigure 5.7: Solution error comparison5.3 2D Advection-Diffusion TestAs a further demonstration of H1 reconstruction, we solve the advection-diffusion equationu · ∇φ = α∇2φwhere α is the diffusion coefficient and u is the velocity vector. The advection-diffusionequation is a model equation which is closer to real world flows, serving as a preliminarystep towards future applications to the Navier-Stokes equations. The inherent stability ofthe upwind discretization of the convective terms is an overall benefit to the stability ofthis problem using H1 reconstruction. We use the same H1 reconstruction to compute theconvective terms as we use for the diffusive terms. The tests were conducted on a series ofunstructured meshes in a rectangular channel of height (H) 1 and length (L) 3, initializedwith smooth exact solution data. The top and bottom wall were enforced with homogeneousisothermal wall boundary conditions. The outlet boundary on the right was maintained atzero flux across the boundary, while at the inlet a velocity profile of sin (piy) was specified.The exact solution for this case can be given by(eR1xR2eR2L − eR2xR1eR1LR2eR2L −R1eR1L)sin (piy) (5.1)575.3. 2D Advection-Diffusion TestFigure 5.8: Advection-diffusion solution (α = 0.01)whereR1 =H2α+√H24α2+ pi2; R2 =H2α−√H24α2+ pi2. (5.2)The final solution is shown in Figure 5.8 with the coefficient of diffusion (α) set equalto 0.01 on an unstructured mesh consisting of 2162 control volumes. This value of 0.01for the coefficient of diffusion causes the flow to be advection dominated while diffusiondominance is achieved with a value of at least 1.The orders of accuracy of the advection-diffusion case, flux and flux integral have beenplotted and compared in Figure 5.9. The tests were run on 5 successively refined meshescontaining 138, 530, 2162, 8655 and 35156 isotropic triangles.The flux integral obtained using the initial data has been shown in Figure 5.9a. Itshows good overall improvement with respect to a second-order least squares scheme. Thediffusion dominant case reflects better order of accuracy while even the advection dominantcase shows accuracy gains by a constant factor of two.The solution convergence in Figure 5.9b, shows that using H1 discretization for thediffusive fluxes provides significant gains in solution accuracy by at least a factor of two,even with an advection dominant flow. In both the advection and the diffusion dominantcases, H1 reconstruction performs at a balance between second and third order least-squaresschemes.Table 5.2 lists the order of accuracy achieved in both advection and diffusion dominantcases for each scheme. Once again the rates of convergence have been averaged over the585.3. 2D Advection-Diffusion Test-2-102 3 4 5log 10( L 2 error in Flux Integral )log10( Number of Control Volumes )1st orderH1 recon2nd Least-squares3rd Least-squares-4-3-2-102 3 4 5log 10( L 2 error in Flux Integral )log10( Number of Control Volumes )1st orderH1 recon2nd Least-squares3rd Least-squaresDiffusion Dominant (α = 1) Advection Dominant (α = 0.01)(a) Flux-integral convergence rates (initial data)-6-5-4-3-22 3 4 5log 10( L 2 error in Solution )log10( Number of Control Volumes )2nd orderH1 recon2nd Least-squares3rd Least-squares-6-5-4-3-2-12 3 4 5log 10( L 2 error in Solution )log10( Number of Control Volumes )2nd orderH1 recon2nd Least-squares3rd Least-squaresDiffusion Dominant (α = 1) Advection Dominant (α = 0.01)(b) Solution convergence ratesFigure 5.9: Advection-Diffusion 2D test results595.3. 2D Advection-Diffusion Testfive refinement stages.Quantity Second-orderLeast SquaresH1ReconstructionThird-orderLeast SquaresDiffusionDominantFlux integral(initial data)0.06 0.51 1.17Solution 1.83 2.01 3.01AdvectionDominantFlux integral(initial data)1.01 1.15 1.63Solution 2.08 2.31 2.50Table 5.2: Orders of accuracy for advection-diffusion problemOverall apart from the accuracy gains, the presence of the upwind convective flux alsogreatly enhances the stability of the H1 discretization, almost playing the part of the jumpterm itself. A more in-depth discussion about stability of the scheme will be presented inSection 6.1. The coefficient of diffusion has also been observed to affect the added stabilitybrought in by the upwind discretization, being largest in the advection dominant situation.60Chapter 6Characterization Tests andBehavioural AnalysisApart from pure numerical accuracy, there are other desirable qualities that one wouldexpect from a numerical scheme. The foremost among these is the issue of stability androbustness which we will look at first. We will also analyze the truncation error that resultsfrom our scheme followed by cost comparison studies. A further opportunity to study thebehaviour of H1 reconstruction arose when the ANSYS Fluent software development teamexpressed their interest in our new scheme. They kindly offered to provide an insight intothe flux discretization scheme used in the Fluent software. This enabled us to performdetailed comparative studies with their discretization scheme as well as possible avenues ofimplementing the H1 finite-volume scheme into their discretization. Our findings from thisstudy have been presented in the last section.6.1 Stability AnalysisIn this section, we analyze the stability of the various schemes discussed. Stability analysisis crucial since it provides valuable insights into the robustness of the scheme on the varietyof complex meshes that are encountered in real world applications.In 2012, Ollivier-Gooch and Roy [42] proposed an eigenvalue analysis tool for unstruc-tured finite-volume schemes. This tool, in addition to its original purpose of analyzing theasymptotic order of the truncation error, can also be used to study the stability of a schemethrough the analysis of the eigenvalues of its Jacobian matrix.The location of the eigenvalues of a space discretization scheme is indicative of thestability of the discretization scheme [27]. Specifically, the stability of the system of ODE’s616.1. Stability Analysisobtained from the method of lines is guaranteed if and only if all its eigenvalues havenegative real components.−16000 −14000 −12000 −10000 −8000 −6000 −4000 −2000 0 2000Re(λ)−600−400−2000200400600Im(λ)(a) Second-order least squares reconstruction with centroids−120000 −100000 −80000 −60000 −40000 −20000 0 20000 40000 60000Re(λ)−10000−50000500010000Im(λ)(b) Second-order least squares reconstruction with circumcentresFigure 6.1: Eigenvalue plots - centroid vs. circumcentresFor instance, the eigenvalue spectrum of a regular second-order least squares recon-struction scheme with respect to its Jacobian matrix, as shown in Figure 6.1a, has allwell-behaved eigenvalues with negative real components and consequently, this is a stable626.1. Stability Analysisproblem. Although the eigenvalues of the Poisson equation should be pure real, that is notthe case here. This can be attributed to the unstructured nature of the underlying meshthat our Jacobian matrix is based upon and as a result it is non-symmetric. Thus the lackof symmetry contributes to the imaginary components of the eigenvalues. This and mostof the other tests that follow were run on a cell-centred unstructured mesh for the Poissonmodel problem, unless specified otherwise.−50000 −40000 −30000 −20000 −10000 0 10000 20000Re(λ)−1500−1000−500050010001500Im(λ)(a) H1 reconstruction with circumcentres−70000 −60000 −50000 −40000 −30000 −20000 −10000 0 10000Re(λ)−800−600−400−2000200400600800Im(λ)(b) H1 reconstruction with containment centres, high accuracy jump termFigure 6.2: H1 reconstruction eigenvalue plots636.1. Stability Analysis−80000 −70000 −60000 −50000 −40000 −30000 −20000 −10000 0 10000Re(λ)−0.010−0.0050.0000.0050.010Im(λ)Figure 6.3: H1 reconstruction on equilateral meshIn contrast Figure 6.1b illustrates the eigenvalue spectrum of the same second-orderleast squares reconstruction but with circumcentres as the control-volume reference pointsinstead of the traditional centroids. We can instantly observe a major degradation of theeigenvalue spectrum, with seven of the eigenvalues lying on the positive real half of the plane(the unstable half), which provides visual proof of the instability of the scheme involvingcircumcentres. Unfortunately H1 reconstruction on a cell-centred unstructured mesh withcircumcentres is also unstable as shown by the positive eigenvalues in Figure 6.2a.We were successful in enhancing the stability of the H1 reconstruction scheme withthe use of containment centres instead of circumcentres as control-volume reference pointsand the inclusion of a high-accuracy jump term. The resulting eigenvalues are plottedin Figure 6.2b. As can be seen, they are all well behaved and lie on the stable half ofthe complex plane. Likewise the eigenvalue spectrum of H1 reconstruction is stable on anequilateral (structured) mesh, as seen in Figure 6.3.Figure 6.4a shows the eigenvalue spectrum of the stable second-order least squaresscheme on the advection-diffusion equation with a coefficient of diffusion set to 0.01. Thestabilizing influence of the upwind discretization of the convective terms is clear when welook at the well behaved eigenvalues of H1 reconstruction with only circumcentres andwithout any jump term, in Figure 6.4b.646.2. Truncation Error Analysis−250 −200 −150 −100 −50 0Re(λ)−80−60−40−20020406080Im(λ)(a) Second-order least squares with centroids (advection-diffusion)−700 −600 −500 −400 −300 −200 −100 0 100Re(λ)−40−2002040Im(λ)(b) H1 reconstruction with circumcentres (advection-diffusion)Figure 6.4: Advention-diffusion eigenvalue plots6.2 Truncation Error AnalysisThe truncation error, as discussed earlier, is the error in approximating the continuousPDEs with the finite discretized equations. It has recently emerged as another measureof accuracy of the scheme and its estimates find varied applications such as improving thesolution through defect correction, finding the error in an engineering output quantity via656.2. Truncation Error Analysisadjoint problems [13], mesh adaptation [48], as well as estimation of discretization errors viathe error transport equation [49]. Ollivier-Gooch and Van Altena [43] performed truncationerror analysis by analyzing the Taylor series of the Laplacian on structured triangularmeshes. Their work was extended to arbitrary meshes and any linear PDE by Jalali andOllivier-Gooch [20]. The result of their analysis is a set of Taylor series coefficients for thetruncation error obtained for every control volume, which when combined with the localsolution derivatives gives the truncation error itself. The same technique is herein used toevaluate and analyze how the H1 finite-volume scheme fares in comparison to the other tworegular least-squares schemes. The truncation error magnitudes were found for the Poissonmodel problem on an unstructured mesh in a square domain.The truncation error is known to be identical to the flux integral or the residual obtainedfrom smooth initial data. Running the exact solution data through the discretization schemeonce provides the finite discrete estimate of the flux integral. This estimate can then becompared to the exact known flux integral to obtain the exact truncation error for thatscheme. Hence the asymptotic convergence rate of truncation error is the same as that forthe flux integral with exact solution data. The flux integral convergence plot for Poissonhas been presented again in Figure 6.5 for convenience. As with the flux integral, theasymptotic convergence rate of truncation error in least-squares schemes is one order lowerthan the gradient and hence two orders lower than the solution. Thus as can be seen inFigure 6.5, the truncation error for second-order least squares is zero order, while that ofthird-order least squares is first order. The H1 finite-volume scheme as we have alreadyseen in section 5.2 performs on par with third-order least squares asymptotically.Figure 6.6 juxtaposes the truncation errors from each of the three schemes being ana-lyzed. The first row shows the truncation error, while the second row plots the magnitudeof truncation error against the real part of the corresponding eigenvalues. A logarithmicscale has been used for both axes in order to enhance the visualization; however taking intoconsideration that the real part of the eigenvalues Re(λ) is negative (for our stable scheme),−Re(λ) is plotted instead. Plotting in such a manner is much more insightful than justlooking at the convergence rates since it provides an explanation for the behaviour of theconvergence [21]. It provides insights into the distribution of the error with respect to the666.2. Truncation Error Analysis-3-2-102 3 4 5log 10( L 2 norm of Truncation Error )log10( Number of Control Volumes )1st orderH1 recon2nd Least-squares3rd Least-squaresFigure 6.5: Truncation error convergence ratehigh and low frequency modes of the eigenvalues. The leftmost eigenvalues correspond tothe high frequency modes while the rightmost ones are the smoother low frequency errors.It is a general trend for the high frequency errors to dominate the eigendecomposition of thetruncation error as evidenced by the truncation error magnitudes of least-squares schemesin Figure 6.6(a) and (c). For instance, if we take a careful look at the plot for the third orderscheme we can observe that the majority of the eigenvalues are on the high-frequency bandto the left and those dominate the truncation error as well. Those high-frequency modesthat occupy the left half of the spectrum correspond to the boundary regions, while theright-most low frequency modes correspond to interior regions. The low frequency eigen-values have considerably smaller truncation error magnitude. Not surprisingly, errors arelarger for the boundary regions.However, H1 reconstruction seems to defy this trend to some extent, with the high as wellas the low frequency modes being of almost similar magnitudes. Yet in overall magnitudeit is still only comparable to second-order least squares scheme. The smoother modesunfortunately seem to have suffered in H1 reconstruction and have higher error. There alsoexist certain outlier high-error eigenvalues which we believe result from the not-so-accurate676.2.TruncationErrorAnalysis12345log10(−Re(λ))-4-3-2-1012log10(Mag(TE))012345log10(−Re(λ))-4-3-2-10112345log10(−Re(λ))-7-6-5-4-3-2-101(a) Second-order least squares (b) H1 finite volume (c) Third-order least squaresFigure 6.6: Truncation error magnitude plots686.3. Cost ComparisonsNumberof CellsConvergence Wall Time perIteration (seconds)Jacobian Assembly Time(seconds)2ndLS H1 3rdLS 2ndLS H1 3rdLS176 1.08×10−2 1.30×10−2 2.17×10−2 3.86×10−3 5.21×10−3 1.09×10−23132 1.45×10−1 1.84×10−1 3.40×10−1 6.74×10−2 9.19×10−2 1.88×10−150715 2.52 3.23 5.66 1.09 1.50 3.01Table 6.1: Wall timescontainment centres.Moreover, if we look at the truncation error visualizations, on the top row of Figure 6.6,then the H1 reconstruction seems to more closely resemble the third order discretization interms of the underlying smoothness but the extreme values are almost four times that ofthird-order least squares reconstruction. Meanwhile, the second-order scheme extremes arealmost a similar factor larger than H1 reconstruction.6.3 Cost ComparisonsNext, we compare the computational cost incurred with the H1 finite-volume scheme againstthe conventional least-squares schemes on the Poisson model problem.Let us begin by analyzing the Jacobian calculations and its memory requirements. Asindicated by Figure 6.7, although the memory required to store the Jacobian matrix is thesame, the total Jacobian matrix assembly time of the new scheme is much lower than a fullthird-order least-squares scheme. The Jacobian memory allocation process accounts onlyfor the size of the reconstruction stencil used and since H1 reconstruction utilizes the samenine cell stencil as the third-order scheme, we would expect their memory requirements tobe identical. However, the H1 scheme uses just one Gauss quadrature point while the third-order least-squares scheme uses two. This translates to savings in the Jacobian assemblytime, making it twice as fast as third-order and only slightly slower than the second-orderJacobian assembly time. This relation is clear when we look at the Jacobian assembly timeslisted in Table 6.1.Moreover, our modifications detailed in section 3.1 are reflected in the convergence wall696.3. Cost Comparisons-3-2-1012 3 4 5log 10( Jacobian Assembly Wall Time ) in secondslog10( Number of Control Volumes )H1 recon2nd Least-squares3rd Least-squares(a) Jacobian calculation wall times34562 3 4 5log 10( Number of non-zeros )log10( Number of Control Volumes )H1 recon2nd Least-squares3rd Least-squares(b) Jacobian memory usageFigure 6.7: Jacobian calculation comparisonstime per iteration of the new scheme which is also lower than a third-order least-squaresscheme as seen in Figure 6.8b. Table 6.1 also lists the convergence wall time per iterationfor each scheme. Unfortunately the H1 scheme currently takes a few additional iterations(pseudo-time steps) to reach full convergence which causes the total convergence time to belarger.All the more insightful is the Figure 6.9 where we can see that for any given accuracy ofthe flux accuracy that we may require (say 10−3) then the quickest way to achieve that is byusing the third-order least-squares scheme followed closely by our H1 scheme and then, byfurther off, the second-order least-squares scheme. We believe that by rectifying the causeof the few additional iterations that H1 scheme requires, we can manage to improve thisplot, by becoming comparable to, if not better than the third-order least-squares scheme.Figure 6.9 also portrays the reason why high-order methods are generally more preferableto low-order schemes.The impact of the savings in Jacobian assembly time is evident when we take a look atthe time profile analysis. The Jacobian evaluation clearly dominates the convergence walltime as seen in Figure 6.10. The pie chart diameters for each scheme is proportional to itsconvergence wall time per iteration. The other prominent processes that make up the rest706.3. Cost Comparisonsof the solver run time include the function evaluation and the KSP solve. The former isthe flux integration process, which includes the evaluation of the fluxes. The KSP solve islinear solver part of our time-advance algorithm, solving the linear system that is obtainedat each time step. The Generalized Minimal Residual (GMRES) algorithm is the Krylovsubspace iterative linear solver that is used in ANSLib [32].-2-10122 3 4 5log 10( Convergence Wall Time ) in secondslog10( Number of Control Volumes )H1 recon2nd Least-squares3rd Least-squares(a) Total convergence time-2-1012 3 4 5log 10( Convergence Wall Time ) in secondslog10( Number of Control Volumes )H1 recon2nd Least-squares3rd Least-squares(b) Convergence time per iterationFigure 6.8: Comparison of total wall clock time required for convergence-5-4-3-2-1-2 -1 0 1 2log 10( L 2 error in Flux )log10( Convergence Wall Time ) in secondsH1 recon2nd Least-squares3rd Least-squaresFigure 6.9: Flux error (steady state) vs total wall clock time required for convergence716.3.CostComparisons26%56%10%8%2nd-order Least-squares26%61%10%3%H1 Reconstruction22%65%8%5%3rd-order Least-squaresFunction Evaluation Jacobian Evaluation KSP(Linear) Solve OthersFigure 6.10: Time profile analysis726.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretization6.4 Behavioural Analysis of the H1 Finite-Volume Schemein the ANSYS Fluent DiscretizationThe interest expressed by the numerics software development team of ANSYS Fluent wasan invaluable opportunity to further study the behaviour of the H1 finite-volume schemewhen implemented within a different flux discretization framework. It also provided anopportunity to baseline our ANSLib numerical simulation software with a well establishedcommercial software like Fluent. We begin with a brief description of the flux discretiza-tion scheme used in the Fluent software, followed by our approach to implement the H1finite-volume scheme within their discretization framework. We will then look at how ourimplementation of the Fluent diffusion discretization scheme within ANSLib matches withthe results from the Fluent software itself. Once satisfied with our baseline calibrations, wewill look at how our modifications fare in comparison with other schemes.6.4.1 Fluent Diffusive Flux DiscretizationThe Fluent flux discretization is based on the scheme of Mathur and Murthy [28], withadditional improvements over the years. Foremost among the modifications is the inclusionof a second-order least-squares solution reconstruction which is very similar to the one usedin ANSLib solver.A central feature of the diffusive flux scheme implemented in Fluent is the “two-pass”gradient computation, which evaluates two distinct forms of gradients. These two gradientcomputations are discussed below followed by the final diffusive flux formulation and adiscussion about how it is handled at the boundaries.6.4.1.1 Reconstruction GradientsThe first pass performs a second-order least-squares solution reconstruction, which providesbase gradients that are only first-order accurate, as anticipated from our discussions inChapter 4. This reconstruction uses the first layer of face-neighbours to form the recon-struction stencil just like the one used in ANSLib second-order least-squares reconstruction.These gradients are only used to reconstruct the solution value at the face mid-point (Gauss736.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretizationpoint) for pass two. Thus at the Gauss point we haveφ˜g = φL +∇φL · ~rLGwhere ~rLG is the distance from the cell centroid to the Gauss point location.6.4.1.2 Diffusion GradientsThe second pass computes a distinct second gradient by applying the divergence theoremover the control volume. ¨CVi∇φdV =˛∂CViφ · ~A (6.1)which upon assuming the gradient ∇φ to be constant over the control volume (first-orderaccurate) can be simplified as∇φCViΩ =˛∂CViφ · ~A (6.2)where Ω is the size of the control volume (area in 2D and volume in 3D).In this divergence formulation, Fluent uses the averaged reconstructed solution valuesat the face. Thus the diffusion gradient for each control volume is formulated as∇φ = 1Ω∑f(φ¯g ~A)(6.3)where ~A is the face area-vector andφ¯g =(φ˜g+ + φ˜g−2). (6.4)This is a generalized formulation for both two and three dimensions, hence face-area vectorsand volumes have been used.Such a formulation extends the stencil that ultimately contributes to the formation ofthe gradients to include both the first and second layers of face-neighbours.746.4. H1 Finite-Volume Scheme in ANSYS Fluent DiscretizationL R- Cell Centroid- Gauss quadrature pointhLRGFigure 6.11: Illustration of Fluent diffusive flux terms6.4.1.3 Diffusive flux formulationThe Fluent diffusive flux algorithm expresses the diffusion term at the face-centreFDf = α(∇φ · ~A)(6.5)where α is the diffusion coefficient, in terms of the solution values (φ) at the control-volumecentroids of the left and right cell. The final diffusive flux formulation is given byFDf = α(φL − φRhLR) ~A · ~A~A · eˆs+ α(∇φ · −→A −∇φ · eˆs~A · ~A~A · eˆs)(6.6)where the averaged diffusion gradient of the two adjacent control volumes is used∇φ =(∇φL +∇φR2)(6.7)and eˆs is the unit-normal vector along the line connecting the two control-volume centroidsas depicted in Figure 6.11.The extended stencil about the face f that results as a consequence of the diffusiongradient can be visualized in Figure 6.12a. In contrast a regular second-order least-squaresflux discretization that was discussed in Section 2.1.3 would result in a first-layer face-neighbour stencil as illustrated in Figure 6.12b.756.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretizationf(a) Fluent schemef(b) Second-order least-squares schemeFigure 6.12: Gradient stencils for diffusive fluxes about face f6.4.1.4 Boundary FluxesThe boundary face centroids are considered as regular cell centroids at the boundaries andhence the diffusive flux takes a similar form as Equation 6.6 yieldingFDf = α(φL − φbhLb) ~A · ~A~A · eˆb+ α(∇φL · −→A −∇φL · eˆb~A · ~A~A · eˆb)(6.8)where φb is the solution value at the boundary face and eˆb is the unit-normal vector con-necting the cell centroid and the boundary face centroid.6.4.2 H1 Scheme in Fluent Diffusive FluxThere are various design preferences and constraints to contend with when implementingchanges within an existing algorithm or scheme. Keeping in mind these constraints, thefollowing choices were made with regards to how we implement our scheme within theirs.1. We retain the use of the modified H1 reconstruction which provides us with secondorder gradients since that was our motivation to begin with.2. We use centroids instead of circumcentres/containment-centres since centroids are theaccepted norm.766.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretization3. We replace just their formulation of the diffusion gradients discussed in subsection 6.4.1.2with the H1 flux discretization as∇φH1 =1(hL + hR){hR.∇φL + hL.∇φR}+O(h2) (6.9)where ∇φL are the gradients obtained from the H1 reconstruction. The remainingformulation of the Fluent diffusive flux is maintained. Hence we haveFDf = α(φL − φRhLR) −→A · −→A−→A · eˆs+ α(∇φH1 ·−→A −∇φH1 · eˆs−→A · −→A−→A · eˆs)(6.10)6.4.3 Comparative Numerical Test ResultsSince we are dealing with two entirely different software packages that have very few thingsin common, to make any comparisons meaningful, we set up calibration tests. That involvedgetting the two codes, Fluent and ANSLib, to work on the same test cases as well as on thesame meshes. The latter, as it turned out, was the tougher bottle-neck for which specificmesh converters had to be written. As for the test cases, the Poisson and advection-diffusionmodel problems were deemed sufficient and the problem set up was maintained exactly thesame as those mentioned in Sections 5.2 and 5.3.We then implemented the Fluent diffusion discretization into ANSLib and performedconvergence studies, while keeping the reconstruction parameters such as inverse-distanceweighting consistent in both pieces of software. Once satisfied with our calibrations, we rancomparison tests to assess the accuracy of various schemes including the H1 scheme withinthe Fluent discretization discussed in the previous section. The following finite-volumeschemes were compared:1. H1 finite-volume scheme (ANSLib-H1),2. Second-order least-squares scheme with jump term (ANSLib-A2L2-wJ),3. Second-order least-squares scheme without jump term (ANSLib-A2L2-nJ),4. Fluent diffusive flux discretization in ANSLib (ANSLib-Fluent), and776.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretization-6-5-4-3-22 3 4 5log 10( L 2 error in solution )log10( Number of Control Volumes )2nd order ref lineFluent-codeANSLib-FluentFigure 6.13: Poisson solution convergence rates in Fluent software (Fluent-code) and Fluentdiscretization implementation in ANSLib (ANSLib-Fluent)5. H1 scheme in Fluent diffusive flux (ANSLib-H1Fluent)The results of these calibration tests as well as the full comparison tests are presented inthe following subsections.6.4.3.1 2D Poisson TestThe Poisson test setup is similar to the one used in Section 5.2, wherein the test functionpi8 sin (pix) sin (piy) and associated source term was used for the method of manufacturedsolutions along with homogeneous Dirichlet boundary conditions and smooth initial solutiondata. Figure 5.3 shows the resultant manufactured solution on a unit square mesh consistingof 780 triangles.Figure 6.13 shows how the Fluent discretization implementation in ANSLib fares com-pared to the actual Fluent software. The Poisson solution accuracy convergence rates havebeen plotted and as can be seen they overlap each other. Thus we can be satisfied withthe comparative studies that follow, which only plot the ANSLib-Fluent discretization for786.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretization-4-3-2-12 3 4 5log 10( L 2 error in flux )log10( Number of Control Volumes )1st order ref lineANSLib-H1ANSLib-A2L2-wJANSLib-A2L2-nJANSLib-FluentANSLib-H1Fluent(a) Flux convergence rate at Gauss points (initialdata)-3-2-102 3 4 5log 10( L 2 error in FI )log10( Number of Control Volumes )1st order ref lineANSLib-H1ANSLib-A2L2-wJANSLib-A2L2-nJANSLib-FluentANSLib-H1Fluent(b) Flux integral convergence rate (initial data)-4-3-2-12 3 4 5log 10( L 2 error in flux )log10( Number of Control Volumes )1st order ref lineANSLib-H1ANSLib-A2L2-wJANSLib-A2L2-nJANSLib-FluentANSLib-H1Fluent(c) Flux convergence rate at Gauss points (convergedsolution)-6-5-4-3-22 3 4 5log 10( L 2 error in solution )log10( Number of Control Volumes )2nd order ref lineANSLib-H1ANSLib-A2L2-wJANSLib-A2L2-nJANSLib-FluentANSLib-H1Fluent(d) Solution convergence rateFigure 6.14: Poisson test results796.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretizationsimplicity.The asymptotic convergence of the diffusive flux obtained using smooth initial data hasbeen shown in Figure 6.14a. We observe that the full H1 reconstruction flux, with second-order accuracy, does better than the rest of the schemes, which are first-order accurate.However our H1 scheme in the Fluent discretization does not fare well, because cell centroidsare used instead of circumcentres. Moreover, the additional terms from Equation 6.10present in the H1-Fluent flux, apart from the H1 gradients, may be adversely affecting theflux accuracy. A similar trend can be seen in the flux integral accuracy in Figure 6.14b.Following that result, we can observe the flux accuracy at the end of complete solutionconvergence in Figure 6.14c. All the schemes here have significantly smaller error than whenusing the exact initial solution (Figure 6.14a). The full H1 reconstruction flux is no longersecond order in the presence of the solution error, though is is still the most accurate. Incontrast, the H1-Fluent discretization performs much better relative to other schemes.All the more interesting is the plot of the converged solution accuracies seen in Fig-ure 6.14d where we can observe that the rankings have almost been reversed save for thesecond-order least-squares with jump term which does exceptionally well.6.4.3.2 2D Advection-Diffusion TestOnce again we set up the same advection-diffusion case that was described in Section 5.3,with an unstructured rectangular-channel mesh.The ANSLib-Fluent implementation does compare well enough with the results fromthe Fluent software as we can see in Figure 6.15. However, variations arise in the advectiondominated flow (α = 0.01) which indicates that the advection flow solution for the two donot quite match up. The advection solution and its corresponding flux are heavily dependenton the solution reconstruction, which is identical in the interior. Hence we can infer that thesecond-order least-squares reconstruction used in the two software packages vary slightlyat the boundaries. Nonetheless, the two results were deemed to be close enough to moveahead with our comparisons.The full comparison test plots for the convergence of the advection-diffusion solutionhave been presented in Figure 6.16. H1 finite-volume is slightly more accurate than the other806.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretization-6-5-4-3-22 3 4 5log 10( L 2 error in solution )log10( Number of Control Volumes )2nd order ref lineFluent-codeANSLib-FluentDiffusion Dominant (α = 1)-6-5-4-3-2-12 3 4 5log 10( L 2 error in solution )log10( Number of Control Volumes )2nd order ref lineFluent-codeANSLib-FluentAdvection Dominant (α = 0.01)Figure 6.15: Advection-Diffusion solution convergence rates in Fluent software and Fluentdiscretization implementation in ANSLib816.4. H1 Finite-Volume Scheme in ANSYS Fluent Discretization-6-5-4-3-22 3 4 5log 10( L 2 error in solution )log10( Number of Control Volumes )2nd order ref lineANSLib-H1ANSLib-A2L2-wJANSLib-A2L2-nJANSLib-FluentANSLib-H1Fluent-6-5-4-3-2-12 3 4 5log 10( L 2 error in solution )log10( Number of Control Volumes )2nd order ref lineANSLib-H1ANSLib-A2L2-wJANSLib-A2L2-nJANSLib-FluentANSLib-H1FluentDiffusion Dominant (α = 1) Advection Dominant (α = 0.01)Figure 6.16: Advection-Diffusion Solution convergence ratesschemes, while the H1-Fluent implementation does well in the advection dominated case.The convective flux seems to have a positive influence on the H1-Fluent implementation.The second-order least-squares without jump term scheme does not perform as well as theothers in both the cases.The presence of the upwind convective flux in the advection-diffusion problem can beseen as a positive influence on the stability and accuracy of the H1 finite-volume scheme aswe saw in Sections 6.1 and 5.3. However in the case of the H1-Fluent scheme, once againthe influence of the unmodified terms in Equation 6.10 seems to have significant effects onthe solution.82Chapter 7Summary and ConclusionsAny implementation of the finite-volume gradients has always been expected to be an orderless accurate than the solution itself. This applies to the diffusive fluxes as well as anyother quantity derived from the gradients. The research described in this thesis was aimedat addressing that concern and we have presented a solution. In the following sections wewill summarize, highlight the major conclusions and also present our visions for the futureexpansions of the research described in this thesis.7.1 SummaryA novel finite-volume scheme to evaluate the diffusive fluxes to second-order accuracy wasdeveloped and presented in this thesis. It has the potential to improve their accuracy,raising it by an entire order even on unstructured meshes and thus counter the loss of anorder in comparison to the solution accuracy. A more accurate diffusive flux should cascadeinto improvements in the flux integral and this has been shown to be true with our results.A finite-volume formulation was utilized as the foundation to develop the new scheme sinceit is currently widely implemented for practical application in the industry. Although moreaccurate high-order schemes have been around for a while, they have not gained traction inthe industry due to the higher computational complexity involved and lack of robustness.Hence the industry standard is only second-order accuracy. This scheme was developed tobridge that gap, among other benefits including a better truncation error estimation.The framework of a general least-squares finite-volume scheme was described in detail.The proposed H1 reconstruction procedure was then described in relation to a second andthird-order least-squares scheme, to provide second-order gradients and without evaluat-ing the second-degree derivatives to save on overall computational resources utilized. The837.1. Summarycircumcentres were used as the primary control-volume reference points, with containment-centres being implemented only when adjacent triangles had nearly coincident circumcen-tres. Second-order diffusive fluxes are then computed from these gradient values withoutthe use of second-degree derivatives using a linear interpolation flux discretization. Theflux integral is obtained from these flux estimates using a single-point Gauss quadrature.We also introduce a high-accuracy jump term to stabilize the scheme yet not degrade itsasymptotic accuracy.Rigorous tests were conducted to verify the proposed scheme as well as to ascertainits characteristics and behaviour. Both analytical and numerical approaches were used asverification procedures.The analytical approach gave us preliminary insights into the accuracy of the schemeon a uniform equilateral triangle grid. On this grid, the computations can be handledanalytically with the aid of symbolic algebra software tools. We analyzed and tweakedthe reconstruction stencil used in our scheme as well as studying the resulting error in theLaplacian (truncation error) while comparing with the second and third-order least-squaresschemes. We were also able to study the decoupling phenomenon that can been seen insecond-order least-squares as well as H1 reconstruction schemes.The numerical approach on the other hand is not restricted by the number of controlvolumes and their topology. Hence it is ideal to study the scheme on unstructured meshes.Preliminary tests were conducted to verify the foundations of our scheme’s linear interpo-lation flux discretization in one dimension as well as the H1 reconstruction accuracy in twodimensions. The Poisson model equation was used as the primary test case to isolate andanalyze the diffusive flux. The Method of Manufactured Solutions (MMS), with a set ofthree exact solutions for Poisson’s equation, was used to evaluate the discretization errorand study the convergence characteristics of the solution, the flux and the flux integral.The expected convergence was verified and compared with those from the least-squaresschemes. The effect of the two alternative implementations of control-volume referencepoints was studied. The implementation and verification of H1 finite-volume scheme in theadvection-diffusion model problem was also studied.Once satisfied with the accuracy and implementation of the scheme, we analyzed its847.2. Conclusionsvarious characteristics and behaviour, concentrating on its stability, truncation error andresource usage. The stability analysis was visualized by studying the eigenvalues spectrumalong the lines of the eigenanalysis strategy proposed by Ollivier-Gooch and Roy [42]. Wewere able to gain much insight into the effect that choice of control-volume reference pointsand the use of a jump term had on the overall stability of different schemes. A truncationerror analysis of our scheme was also performed using methods proposed by Jalali andOllivier-Gooch [20]. We compared the truncation error of our scheme with that of secondand third-order least-squares schemes. Finally an in-depth cost analysis was conductedto verify the advantages gained over other least-squares schemes and possible avenues ofimprovements.The behaviour of H1 finite-volume scheme within the Fluent diffusive flux discretizationscheme was also explored. The Fluent diffusive discretization was implemented withinANSLib and baseline tests were performed to ensure the implementation worked as well asthe actual Fluent software itself. Several schemes and their variations were then comparedon the chosen model problems of Poisson and advection-diffusion to study and compare thebehaviour of each.7.2 ConclusionsAlthough the research and development of a new finite-volume scheme is challenging andinevitably fraught with unforeseen complications, we are happy to be able to present onewhich does fulfil the expectations.By means of the analytical tests we were able to verify the reconstruction order of accu-racy for linear (second-order) least-squares reconstruction on cell-centred uniform equilat-eral triangle meshes. However we do observe inconsistencies in the quadratic reconstructionwhich has a leading-order-error term which is first-order resulting from degeneracy in recon-struction on uniform triangular mesh. H1 reconstruction being a derivative of the quadraticreconstruction, seems to inherit this inconsistency along with an unforeseen solution decou-pling unique to H1 reconstruction.Fortunately the numerical tests, conducted with cell-centred triangular unstructured857.2. Conclusionsmeshes on Poisson and advection-diffusion problems, show more promising results. Anoverview of the results convey that the diffusive fluxes do indeed have an improved order ofaccuracy while the solution itself is more accurate than second-order least-squares schemeby a constant factor, albeit a small factor. The flux integral at zero iterations, which hasmuch smaller error, shows much improvement in order for the case where we replace onlythe “ugly circumcentres” with containment centres. In addition, the solution error is alsomuch smoother and lies between that of second and third order least-squares schemes interms of smoothness.The complication of stability was by far the toughest hurdle to overcome with very fewoptions available that would not adversely sacrifice the gains in accuracy of the diffusiveflux. The use of circumcentres as the control-volume reference points turned out to be ahighly unstable choice. The solution presented itself in the form of containment centres andmore importantly the high-accuracy jump term.The truncation error again showed some improvement compared to the second-orderleast-squares reconstruction, although not as much as we would have wished for. Thehigh and low frequency modes of the eigenvalues appear to share similar magnitude in thetruncation error for H1 reconstruction.As for computational resource usage, as anticipated, H1 reconstruction is cheaper thana full third-order least-squares reconstruction. The Jacobian evaluation time is twice asfast as the third-order scheme and almost on par with the second-order scheme. Howeverthe memory requirement is identical to the third-order scheme. This is expected sincethe memory requirement is essentially only dependent on the number of neighbours in thereconstruction stencil and H1 reconstruction uses the exact same stencil as the third-orderscheme. These results also translate into the total wall convergence time per iteration sincethe Jacobian evaluation time is the largest single part of the computational time. Howeverat the moment H1 reconstruction takes a couple more iterations to converge than either ofthe least-squares schemes and this adversely affects the overall time taken for convergence.The study of the Fluent diffusive flux discretization was a very insightful experience.Every discretization scheme has its unique advantages and disadvantages, so comparingthese idiosyncrasies and their causes is overall always beneficial. For instance, the extended867.3. Recommended Future Workstencil that results from the use of the distinct Fluent diffusion gradients and its inferredeffect on the stability of the scheme was enlightening. The implementation of the Fluentdiscretization into ANSLib was completely successful and the matching results reinforcedthe similarities of other underlying algorithms. However combining the H1 scheme with theFluent discretization provided mixed success that was not entirely as we anticipated. Theseemingly inverse relation between the flux accuracy and solution accuracy convergence forthe Poisson scheme was the highlight of those results and its cause will be probed furtherin the near future.In conclusion, looking at the whole picture, the proposed new H1 finite-volume schememeets our research objectives of striking a balance between the best of second and thirdorder finite-volume schemes.7.3 Recommended Future WorkOur work with the H1 finite-volume scheme paves the way for further analysis, extensionsand applications.1. Work can be geared towards improving the convergence characteristics as well asimproving the robustness of the scheme.2. The current research has only focused on the interior regions of the mesh. Implementa-tion in the boundary regions is another important step. The boundary discretizationis the most challenging part of any numerical scheme and thus the impact of theboundaries on the solution and its gradients cannot be neglected. A possible route tothe same could be through the effective usage of ghost cells.3. We could also look at possible implementations for vertex-centred meshes. Vertex-centred mesh discretizations have the advantage of being slightly more robust thantheir cell-centred counterparts due to more numerous first neighbours, as well ascompact stencils for third-order (and hence for H1 reconstruction).4. The present research was limited to second-order solution and gradients. The ex-tension to higher order would be more challenging yet potentially lucrative. The877.3. Recommended Future Workanticipated challenge would be dealing with the two Gauss point values required forhigher order accuracy.5. Extension from two dimensions to three dimensions is another avenue to look at.Using circumcentres of tetrahedrons should be straightforward. However handlingissues of stability in three dimensions could be a significant challenge.6. Future explorations could include extension from simple model problems to the morecomplex real fluid flows defined by the Navier-Stokes equations, possibly in combi-nation with existing viscous flux discretization schemes, and study of their combinedcharacteristics.7. Another area of future work is application to problems with gradient-based sourceterms, a class which includes the production term for many RANS turbulence mod-els including the Spalart-Allmaras turbulence models [53, 52]. The requirements forimproved accuracy in source term evaluation and its impact on solution accuracy willbe of primary interest.88Bibliography[1] Rémi Abgrall. On essentially non-oscillatory schemes on unstructured meshes: Analysisand implementation. Journal of Computational Physics, 114(1):45–58, 1994.[2] John D. Anderson. Computational Fluid Dynamics: The Basics with Applications.McGraw-Hill Education, 1995.[3] Timothy J. Barth. Aspects of unstructured grids and finite-volume solvers for the Eulerand Navier-Stokes equations. In Unstructured Grid Methods for Advection-DominatedFlows, pages 6–1 – 6–61. AGARD, Neuilly sur Seine, France, 1992. AGARD-R-787.[4] Timothy J. Barth. Recent developments in high order k-exact reconstruction on un-structured meshes. AIAA paper 93-0668, January 1993.[5] Timothy J. Barth and Herman Deconinck. High-order methods for computationalphysics, volume 9. Springer Science & Business Media, 1999.[6] Timothy J. Barth and Paul O. Frederickson. Higher order solution of the Euler equa-tions on unstructured grids using quadratic reconstruction. AIAA paper 90-0013, Jan-uary 1990.[7] Timothy J. Barth and Dennis C. Jespersen. The design and application of upwindschemes on unstructured meshes. AIAA paper 89-0366, January 1989.[8] Timothy J. Barth and Mario Ohlberger. Finite Volume Methods: Foundation andAnalysis. John Wiley & Sons, Ltd, 2004.[9] Alexander N. Brooks and Thomas J. R. Hughes. Streamline upwind/Petrov-Galerkinformulations for convection dominated flows with particular emphasis on the incom-89Bibliographypressible Navier-Stokes equations. Computer methods in applied mechanics and engi-neering, 32(1):199–259, 1982.[10] Bernardo Cockburn and Chi-Wang Shu. Runge–Kutta discontinuous Galerkin methodsfor convection-dominated problems. Journal of scientific computing, 16(3):173–261,2001.[11] Boris Diskin, James L. Thomas, Eric J. Nielsen, Hiroaki Nishikawa, and Jeffrey A.White. Comparison of node-centered and cell-centered unstructured finite-volume dis-cretizations: Viscous fluxes. American Institute of Aeronautics and Astronautics Jour-nal, 48(7):1326–1338, July 2010.[12] Olivier Friedrich. Weighted essentially non-oscillatory schemes for the interpolation ofmean values on unstructured grids. Journal of Computational Physics, 144(1):194–212,July 1998.[13] Michael B. Giles and Niles A. Pierce. An introduction to the adjoint approach todesign. Flow, Turbulence, and Combustion, 65(3–4):393–415, 2000.[14] Ami Harten, Björn Enquist, Stanley Osher, and Sukumar R. Chakravarthy. Uniformlyhigh order accurate essentially non-oscillatory schemes, III. Journal of ComputationalPhysics, 71(2):231–303, August 1987.[15] Jan S. Hesthaven and Tim Warburton. Nodal discontinuous Galerkin methods: algo-rithms, analysis, and applications. Springer Science & Business Media, 2007.[16] Thomas J. R. Hughes, Michel Mallet, and Mizukami Akira. A new finite elementformulation for Computational Fluid Dynamics: II. beyond SUPG. Computer Methodsin Applied Mechanics and Engineering, 54(3):341–355, 1986.[17] H.T. Huynh. A flux reconstruction approach to high-order schemes including discon-tinuous Galerkin methods. AIAA paper, 4079:2007, 2007.[18] H.T. Huynh, Z. J. Wang, and P.E. Vincent. High-order methods for computationalfluid dynamics: A brief review of compact differential formulations on unstructuredgrids. Computers & Fluids, 98:209–220, 2014.90Bibliography[19] Alireza Jalali. Truncation error analysis of unstructured finite volume discretizationschemes. Master’s thesis, The University of British Columbia, Department of Mechan-ical Engineering, 2012.[20] Alireza Jalali and Carl Ollivier-Gooch. Accuracy assessment of finite volume discretiza-tions of diffusive fluxes on unstructured meshes. In Proceedings of the Fiftieth AIAAAerospace Sciences Meeting. American Institute of Aeronautics and Astronautics, 2012.AIAA Paper 2012-0608.[21] Alireza Jalali, Mahkame Sharbatdar, and Carl Ollivier-Gooch. Accuracy analysis ofunstructured finite volume discretization schemes for diffusive fluxes. Computers &Fluids, 101:220–232, 2014.[22] Aaron Katz and Venkateswaran Sankaran. An efficient correction method to obtain aformally third-order accurate flow solver for node-centered unstructured grids. Journalof Scientific Computing, 51(2):375–393, 2012.[23] David A. Kopriva and John H. Kolias. A conservative staggered-grid Chebyshev mul-tidomain method for compressible flows. Journal of Computational Physics, 125(1):244– 261, 1996.[24] Sanjiva K. Lele. Compact finite difference schemes with spectral-like resolution. Journalof computational physics, 103(1):16–42, 1992.[25] Xu-Dong Liu, Stanley Osher, and Tony Chan. Weighted essentially non-oscillatoryschemes. Journal of Computational Physics, 115:200–212, 1994.[26] Yen Liu, Marcel Vinokur, and ZJ Wang. Spectral difference method for unstructuredgrids i: basic formulation. Journal of Computational Physics, 216(2):780–801, 2006.[27] Harvard Lomax, Thomas H. Pulliam, and David W. Zingg. Fundamentals of Compu-tational Fluid Dynamics. Springer Science & Business Media, 2001.[28] S. R. Mathur and J. Y. Murthy. A pressure-based method for unstructured meshes.Numerical Heat Transfer, 31(2):195–215, 1997.91Bibliography[29] Dimitri J. Mavriplis. Revisiting the least-squares procedure for gradient reconstructionon unstructured meshes. In Proceedings of the Sixteenth AIAA Computational FluidDynamics Conference, 2003.[30] Christopher Michalak. Efficient high-order accurate unstructured finite-volume algo-rithms for viscous and inviscid compressible flows. PhD thesis, The University ofBritish Columbia, Department of Mechanical Engineering, 2009.[31] Christopher Michalak and Carl Ollivier-Gooch. Globalized matrix-explicit Newton-GMRES for the high-order accurate solution of the Euler equations. Computers andFluids, 39:1156–1167, 2010.[32] Amir Nejat. A Higher-Order Accurate Unstructured Finite Volume Newton-KrylovAlgorithm for Inviscid Compressible Flows. PhD thesis, The University of BritishColumbia, Department of Mechanical Engineering, 2007.[33] Amir Nejat and Carl Ollivier-Gooch. A high-order accurate unstructured finite volumeNewton-Krylov algorithm for inviscid compressible flows. Journal of ComputationalPhysics, 227(4):2592–2609, 2008.[34] Hiroaki Nishikawa. A first-order system approach for diffusion equation. I: Second-order residual-distribution schemes. Journal of Computational Physics, 227(1):315–352, 2007.[35] Hiroaki Nishikawa. Beyond interface gradient: A general principle for constructingdiffusion schemes. In Proceedings of the Fortieth AIAA Fluid Dynamics Conference,2010.[36] Hiroaki Nishikawa. A first-order system approach for diffusion equation. II: Unificationof advection and diffusion. Journal of Computational Physics, 229(11):3989–4016, 2010.[37] Carl Ollivier-Gooch, Amir Nejat, and Christopher Michalak. On obtaining and verify-ing high-order finite-volume solutions to the Euler equations on unstructured meshes.American Institute of Aeronautics and Astronautics Journal, 47(9):2105–2120, 2009.92Bibliography[38] Carl F. Ollivier-Gooch. High-order ENO schemes for unstructured meshes based onleast-squares reconstruction. AIAA paper 97-0540, January 1997.[39] Carl F. Ollivier-Gooch. Quasi-ENO schemes for unstructured meshes based on unlim-ited data-dependent least-squares reconstruction. Journal of Computational Physics,133(1):6–17, 1997.[40] Carl F. Ollivier-Gooch. An unstructured mesh improvement toolkit with applicationto mesh improvement, generation and (de-)refinement. AIAA 98-0218, January 1998.[41] Carl F. Ollivier-Gooch. A toolkit for numerical simulation of PDE’s: I. Fundamen-tals of generic finite-volume simulation. Computer Methods in Applied Mechanics andEngineering, 192:1147–1175, February 2003.[42] Carl F. Ollivier-Gooch and Christopher Roy. Reducing truncation error on unstruc-tured meshes by vertex movement. In Proceedings of the Forty-Second AIAA FluidDynamics Conference, 2012.[43] Carl F. Ollivier-Gooch and Michael Van Altena. A high-order accurate unstructuredmesh finite-volume scheme for the advection-diffusion equation. Journal of Computa-tional Physics, 181(2):729–752, 2002.[44] B. Pincock and A. Katz. High-order flux correction for viscous flows on arbitraryunstructured grids. Journal of Scientific Computing, 61(2):454–476, 2014.[45] Thomas H Pulliam and Joseph L Steger. Implicit finite-difference simulations of three-dimensional compressible flow. AIAA Journal, 18(2):159–167, 1980.[46] Varun P. Puneria. Truncation error analysis for diffusion schemes in boundary regions ofunstructured meshes. Master’s thesis, The University of British Columbia, Departmentof Mechanical Engineering, 2015.[47] Varun P. Puneria and Carl Ollivier-Gooch. Truncation error analysis for diffusionschemes in boundary regions of unstructured meshes. In Proceedings of the Twenty-Second Annual Conference of the Computational Fluid Dynamics Society of Canada.Société canadienne de CFD / CFD Society of Canada, 2014.93Bibliography[48] Christopher Roy. Strategies for driving mesh adaptation in CFD. In Proceedings of theForty-Seventh AIAA Aerospace Sciences Meeting. American Institute of Aeronauticsand Astronautics, 2009.[49] Christopher Roy. Review of discretization error estimators in scientific computing. InProceedings of the Forty-Eighth AIAA Aerospace Sciences Meeting. American Instituteof Aeronautics and Astronautics, 2010. AIAA Paper 2010-0126.[50] Mahkame Sharbatdar and Carl Ollivier-Gooch. Smooth estimate of the truncation errorfor unstructured mesh finite volume methods. In Proceedings of the Sixth EuropeanConference on Computational Fluid Dynamics, 2014.[51] Chi-Wang Shu. High-order finite difference and finite volume WENO schemes anddiscontinuous galerkin methods for CFD. International Journal Computational FluidDynamics, 17:107–118, 2003.[52] P. R. Spalart. Strategies for turbulence modelling and simulations. InternationalJournal of Heat and Fluid Flow, 21(3):252 – 263, 2000.[53] P. R. Spalart and S. R. Allmaras. A one-equation turbulence model for aerodynamicflows. In Proceedings of the Thirtieth AIAA Aerospace Sciences Meeting, January 1992.[54] James L. Thomas, Boris Diskin, and Christopher Rumsey. Towards verification ofunstructured-grid solvers. In Forty-sixth AIAA Aerospace Sciences Meeting, 2008.AIAA 2008-666.[55] Michael Van Altena. High-order finite-volume discretisations for solving a modifiedadvection-diffusion problem on unstructured triangular meshes. Master’s thesis, TheUniversity of British Columbia, Department of Mechanical Engineering, October 1999.[56] V. Venkatakrishnan. A perspective on unstructured grid flow solvers. Technical report,ICASE 95-3 NASA, 1995.[57] Henk Kaarle Versteeg and Weeratunge Malalasekera. An introduction to computationalfluid dynamics: the finite volume method. Pearson Education, 2007.94Bibliography[58] P. Vincent, P. Castonguay, and A. Jameson. A new class of high-order energy stableflux reconstruction schemes. Journal of Scientific Computing, 47(1):50–72, 2011.[59] Z. J. Wang, Krzysztof Fidkowski, Rémi Abgrall, Francesco Bassi, Doru Caraeni, An-drew Cary, Herman Deconinck, Ralf Hartmann, Koen Hillewaert, HT Huynh, et al.High-order CFD methods: current status and perspective. International Journal forNumerical Methods in Fluids, 72(8):811–845, 2013.[60] Gary Kai-Kin Yan, Varun P. Puneria, Alireza Jalali, and Carl Ollivier-Gooch. Trun-cation and discretization error for diffusion schemes on unstructured meshes. In Pro-ceedings of the Fifty-Second AIAA Aerospace Sciences Meeting. American Institute ofAeronautics and Astronautics, 2014.95
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Improving finite-volume diffusive fluxes through better...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Improving finite-volume diffusive fluxes through better reconstruction Sejekan, Chandan Balachandra 2016
pdf
Page Metadata
Item Metadata
Title | Improving finite-volume diffusive fluxes through better reconstruction |
Creator |
Sejekan, Chandan Balachandra |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | The overarching goal of CFD is to compute solutions with low numerical error. For finite-volume schemes, this error originates as error in the flux integral. For diffusion problems on unstructured meshes, the diffusive flux (computed from reconstructed gradients) is one order less accurate than the reconstructed solution. Worse, the gradient errors are not smooth, and so no error cancellation accompanies the flux integration, reducing the flux integral to zero order for second-order schemes. Our aim is to compute the gradient and flux more accurately at the cell boundaries and hence obtain a better flux integral for a slight increase in computational cost. We propose a novel reconstruction method and flux discretization to improve diffusive flux accuracy on cell-centred, isotropic unstructured meshes. Our approach uses a modified least-squares system to reconstruct the solution to second-order accuracy in the H₁ norm instead of the prevalent L₂ norm, thus ensuring second-order accurate gradients. Either circumcentres or containment centres are chosen as the control-volume reference points based on a criteria to facilitate calculation of second-order gradients at flux quadrature points using a linear interpolation scheme along with a high-accuracy jump term to enhance stability of the system. Numerical results show a significant improvement in the order of accuracy of the computed diffusive flux as well as the flux integral. When applied to a channel flow advection-diffusion problem, the scheme resulted in an increased order of accuracy for the flux integral along with gains in solution accuracy by a factor of two. The characteristics of the new scheme were studied through stability, truncation error and cost analysis. The increase in computational costs were modest and affordable. The behaviour of the scheme was also tested by implementing a variation of it within the ANSYS Fluent discretization framework. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-04-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0300049 |
URI | http://hdl.handle.net/2429/57719 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2016-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2016_may_sejekan_chandan.pdf [ 6.09MB ]
- Metadata
- JSON: 24-1.0300049.json
- JSON-LD: 24-1.0300049-ld.json
- RDF/XML (Pretty): 24-1.0300049-rdf.xml
- RDF/JSON: 24-1.0300049-rdf.json
- Turtle: 24-1.0300049-turtle.txt
- N-Triples: 24-1.0300049-rdf-ntriples.txt
- Original Record: 24-1.0300049-source.json
- Full Text
- 24-1.0300049-fulltext.txt
- Citation
- 24-1.0300049.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0300049/manifest