"Applied Science, Faculty of"@en . "Mechanical Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Hoshyari, Shayan"@en . "2017-08-28T22:21:51Z"@en . "2017"@en . "Master of Applied Science - MASc"@en . "University of British Columbia"@en . "High-order accurate numerical discretization methods are attractive for their potential to significantly reduce the computational costs compared to the traditional second-order methods. Among the various unstructured higher-order discretization schemes, the k-exact reconstruction finite volume method is of interest for its straightforward mathematical formulation, and its compatibility with the current lower-order industrial solvers. However, current three-dimensional finite volume solvers are limited to the solution of inviscid and laminar viscous flow problems. Since three-dimensional turbulent flows appear in many industrial applications, the current thesis takes the first step towards the development of a three-dimensional higher-order finite volume solver for the solution of both inviscid and viscous turbulent steady-state flow problems. The k-exact finite volume formulation of the governing equations is rederived in a dimension-independent manner, where the negative Spalart-Allmaras turbulence model is employed. This one-equation model is reasonably accurate for many flow conditions, and its simplicity makes it a good starting point for the development of numerical algorithms. Then, the three-dimensional mesh preprocessing steps for a finite volume simulation are presented, including higher-order accurate numerical quadrature, and capturing the boundary curvature in highly anisotropic meshes. Also, the issues of k-exact reconstruction in handling highly anisotropic meshes are reviewed and addressed. Since three-dimensional problems can require much more memory than their two-dimensional counter-parts, solution methods that work in two dimensions might not be feasible in three dimensions anymore. As an attempt to overcome this issue, a practical and parallel scalable method for the solution of the discretized system of nonlinear equations is presented. Finally, the solution of four three-dimensional test problems are studied: Poisson\u00E2\u0080\u0099s equation in a cubic domain, inviscid flow over a sphere, turbulent flow over a flat plate, and turbulent flow over an extruded NACA 0012 airfoil. The solution is verified, and the resource consumption of the flow solver is measured. The results demonstrate the benefit and practicality of using higher-order methods for obtaining a certain level of accuracy."@en . "https://circle.library.ubc.ca/rest/handle/2429/62846?expand=metadata"@en . "A Higher-Order Unstructured Finite Volume Solver forThree-Dimensional Compressible FlowsbyShayan HoshyariB.Sc., Sharif University of Technology, 2015a thesis submitted in partial fulfillmentof the requirements for the degree ofMASTER OF APPLIED SCIENCEinthe faculty of graduate and postdoctoral studies(Mechanical Engineering)The University of British Columbia(Vancouver)August 2017\u00C2\u00A9 Shayan Hoshyari, 2017AbstractHigh-order accurate numerical discretization methods are attractive for their potential to significantlyreduce the computational costs compared to the traditional second-order methods. Among the variousunstructured higher-order discretization schemes, the k-exact reconstruction finite volume method is ofinterest for its straightforward mathematical formulation, and its compatibility with the current lower-order industrial solvers. However, current three-dimensional finite volume solvers are limited to thesolution of inviscid and laminar viscous flow problems. Since three-dimensional turbulent flows appearin many industrial applications, the current thesis takes the first step towards the development of a three-dimensional higher-order finite volume solver for the solution of both inviscid and viscous turbulentsteady-state flow problems.The k-exact finite volume formulation of the governing equations is rederived in a dimension-independentmanner, where the negative Spalart-Allmaras turbulence model is employed. This one-equation modelis reasonably accurate for many flow conditions, and its simplicity makes it a good starting point for thedevelopment of numerical algorithms. Then, the three-dimensional mesh preprocessing steps for a finitevolume simulation are presented, including higher-order accurate numerical quadrature, and capturingthe boundary curvature in highly anisotropic meshes. Also, the issues of k-exact reconstruction inhandling highly anisotropic meshes are reviewed and addressed.Since three-dimensional problems can require much more memory than their two-dimensional counter-parts, solution methods that work in two dimensions might not be feasible in three dimensions anymore.As an attempt to overcome this issue, a practical and parallel scalable method for the solution of thediscretized system of nonlinear equations is presented.Finally, the solution of four three-dimensional test problems are studied: Poisson\u00E2\u0080\u0099s equation in a cubicdomain, inviscid flow over a sphere, turbulent flow over a flat plate, and turbulent flow over an extrudedNACA0012 airfoil. The solution is verified, and the resource consumption of the flow solver is measured.The results demonstrate the benefit and practicality of using higher-order methods for obtaining a certainlevel of accuracy.iiLay SummaryThe finite volume method is a popular numerical scheme for the solution of aerodynamic flow problems.Although conventional finite volume schemes aremostly second-order accurate, there has been a growinginterest in higher-order accurate numerical methods, since they can be considerably more efficient interms of computational resources for achieving a certain level of accuracy.Higher-order finite volumemethods have been successfully developed for awide range of two-dimensionalflow problems. Nevertheless, numerical methods that work in two dimensions might not be feasible inthree dimensions anymore, since three-dimensional problems can require much more memory than theirtwo-dimensional counterparts. The current thesis aims at identifying, and resolving such shortcomingsto construct a working three-dimensional finite volume solver with an emphasis on turbulent flows. Sub-sequently, three-dimensional benchmark flow problems are solved to verify and assess the performanceof the developed numerical method.iiiPrefaceAll the work presented in this thesis is an intellectual product of a working relationship between ShayanHoshyari and Dr. Carl Ollivier-Gooch. The implementation of the methods, the data analysis, and themanuscript preparations were done by Shayan Hoshyari with invaluable guidance from Carl Ollivier-Gooch throughout the process.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 The Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 K-exact Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3 Studied Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3.1 Poisson\u00E2\u0080\u0099s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3.2 Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3.3 Extension to Turbulent Flows . . . . . . . . . . . . . . . . . . . . . . . . . . 112.4 Numerical Flux Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4.1 Inviscid Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4.2 Viscous Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 Three-Dimensional Mesh Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17v3.1 Element Mapping and Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.2 Creating Curved Anisotropic Meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3 Modified Basis Functions for Highly Anisotropic Meshes . . . . . . . . . . . . . . . . 244 Solving the Discretized System of Equations . . . . . . . . . . . . . . . . . . . . . . . . . 284.1 Pseudo Transient Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.2 Linear Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.3 Preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.3.1 Point Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.3.2 Block Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.3.3 ILU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.4 Improved Preconditioning Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 344.4.1 Inner GMRES Iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.4.2 Lines of Strong Coupling Between Unknowns . . . . . . . . . . . . . . . . . . 354.5 Numerical Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365 Three-Dimensional Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.1 Poisson\u00E2\u0080\u0099s Equation in a Cubic Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 415.2 Inviscid Flow Around a Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425.3 Turbulent Flow Over a Flat Plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485.4 Turbulent Flow Over an Extruded Airfoil . . . . . . . . . . . . . . . . . . . . . . . . . 516 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60A ANSLib Command-Line Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66A.1 Two-Dimensional Turbulent Flow over a NACA 0012 Airfoil from Section 4.5 . . . . . 66A.2 Poisson\u00E2\u0080\u0099s Equation from Section 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 69A.3 Inviscid Flow Over a Sphere from Section 5.2 . . . . . . . . . . . . . . . . . . . . . . 69A.4 Turbulent Flow Over a Flat Plate from Section 5.3 . . . . . . . . . . . . . . . . . . . . 70A.5 Turbulent Flow Over an Extruded NACA 0012 Airfoil from Section 5.4 . . . . . . . . 71B Sample Script for Running a Parallel Job on Grex . . . . . . . . . . . . . . . . . . . . . 73viList of TablesTable 2.1 Dimensionless empirical parameters used in the negative S-A model . . . . . . . . 13Table 3.1 Performance of the linear elasticity solver . . . . . . . . . . . . . . . . . . . . . . . 24Table 4.1 Preconditioning methods considered . . . . . . . . . . . . . . . . . . . . . . . . . 38Table 4.2 Performance comparison for different preconditioning schemes . . . . . . . . . . . 40Table 4.3 Computed drag and lift coefficients for the NACA 0012 airfoil . . . . . . . . . . . . 40Table 5.1 Number of iterations and the resource consumption of the flow solver for the sphereproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47Table 5.2 Computed value and convergence order of the drag coefficient and the skin frictioncoefficient at the point x = (0.97, 0, 0.5) for the flat plate problem . . . . . . . . . . 51Table 5.3 Number of iterations and the resource consumption of the flow solver for the flat plateproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52Table 5.4 Computed drag and lift coefficients for the extruded NACA 0012 airfoil problem . . 54Table 5.5 Resource consumption of the flow solver for the extruded NACA 0012 problem . . . 56viiList of FiguresFigure 2.1 Illustration of terminology in finite volume discretization . . . . . . . . . . . . . . 6Figure 2.2 Reconstruction stencils for a control volume . . . . . . . . . . . . . . . . . . . . . 9Figure 3.1 Control volumes for cell-centered and cell-vertex methods . . . . . . . . . . . . . 18Figure 3.2 First order reference Lagrange elements . . . . . . . . . . . . . . . . . . . . . . . 19Figure 3.3 Geometry of the curved mesh generation test case . . . . . . . . . . . . . . . . . . 23Figure 3.4 Curved mesh displacement magnitude . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 3.5 Residual per iteration for the linear elasticity solver . . . . . . . . . . . . . . . . . 25Figure 3.6 Construction of approximate wall coordinates for a two-dimensional control volume 26Figure 4.1 Schematic of Block Jacobi decomposition . . . . . . . . . . . . . . . . . . . . . . 33Figure 4.2 Mesh and lines of strong unknown coupling for the geometry of a NACA 0012 airfoil 38Figure 4.3 Comparison of residual histories for the NACA 0012 problem obtained using differ-ent preconditioning algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 4.4 Residual history for the NACA 0012 problem on the finest mesh of 400K controlvolumes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 5.1 Exact solution of Poisson\u00E2\u0080\u0099s problem . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 5.2 Different meshes for the solution of Poisson\u00E2\u0080\u0099s problem . . . . . . . . . . . . . . . 43Figure 5.3 Discretization error versus mesh size for Poisson\u00E2\u0080\u0099s problem . . . . . . . . . . . . . 44Figure 5.4 Symmetric cut of the coarsest mesh near the sphere surface . . . . . . . . . . . . . 44Figure 5.5 Relative entropy norm versus mesh size for the sphere problem . . . . . . . . . . . 45Figure 5.6 Computed Mach contours on the x3 = 0 symmetry plane for the sphere problem . . 46Figure 5.7 Norm of the residual vector per PTC iteration for the sphere problem . . . . . . . . 47Figure 5.8 The parallel speedup of the solver for the sphere problem . . . . . . . . . . . . . . 48Figure 5.9 Coarsest mesh for the flat plate problem. . . . . . . . . . . . . . . . . . . . . . . . 49Figure 5.10 Distribution of the turbulence working variable on the plane x3 = 0.5 for the flatplate problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 5.11 Distribution of the nondimensional eddy viscosity on the line (x1 = 0.97)\u00E2\u0088\u00A7(x3 = 0.5)for the flat plate problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 5.12 Norm of the residual vector per PTC iteration for the flat plate problem . . . . . . 51viiiFigure 5.13 The parallel speedup of the solver for the flat plate problem . . . . . . . . . . . . . 52Figure 5.14 Meshes for the extruded NACA 0012 problem . . . . . . . . . . . . . . . . . . . . 53Figure 5.15 Distribution of the turbulenceworking variable for the extrudedNACA0012 problemon the x3 = 0 plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54Figure 5.16 Distribution of surface pressure coefficient on the intersection of the extruded NACA0012 airfoil and the x3 = 0.5 plane. . . . . . . . . . . . . . . . . . . . . . . . . . 55Figure 5.17 Norm of the residual vector per PTC iteration for the extruded airfoil problem . . . 55ixList of SymbolsOrdersk Order of reconstructionl Order of Lagrangian polynomialsp Order of finite element basis functionsq Order of accuracy for a quadrature ruleNumber of EntitiesNunk Number of unknownsNDOF Number of degrees of freedomNdim Number of dimensionsNrec Number of reconstruction basis functionsNlag Number of Lagrangian basis functionsNqua Number of quadrature pointsNCV Number of control volumesNumerical Discretizationx Cartesian coordinatesL Differential operatorR Set of rational numbersu Solution vectorxh Mesh representative length scaleTh Set of all control volumesuh Discrete solution vectorUh Vector of degrees of freedomeh Discretization errorR Residual vector\u00E2\u0084\u00A6 Domain\u00E2\u0088\u0082\u00E2\u0084\u00A6 Boundary of domainn Surface unit normal\u00CF\u0084 A control volume\u00E2\u0088\u0082\u00CF\u0084 Boundary of a control volume\u00E2\u0084\u00A6\u00CF\u0084 Volume or area of a control volume\u00CF\u0088 A Finite element or Lagrangian basis function\u00CF\u0086 A Cartesian coordinates monomial\u00CF\u0095\u00E2\u0088\u0097 A Principal coordinates monomial\u00CF\u0095+ A Curvilinear coordinates monomialx\u00CF\u0084 Reference location of a control volumew\u00CF\u0084 Principal coordinates for a control volume\u00CE\u00BE\u00CF\u0084 Curvilinear coordinates for a control volumet\u00CF\u0084 Tangent to wall coordinate for a control volumed\u00CF\u0084 Normal to wall coordinate for a control volumeF Inviscid flux matrixQ Viscous flux matrixS Source term vectorF Numerical inviscid flux functionQ Numerical viscous flux functionxiB Boundary condition constraint functionB Boundary conditionsMinNeigh Minimum number of reconstruction stencil membersBdryQuad Set of all boundary quadrature points for a control volumeStencil Reconstruction stencil for a control volumeSolution of System of Equationsp ILU fill levelx Vector of unknownsb Right hand side vector for a linear systemr Residual vector for a linear systemA Left hand side matrix for a linear systemA\u00E2\u0088\u0097 The matrix used for constructing the preconditionerP Preconditioner matrixL\u00CB\u009C, U\u00CB\u009C Lower and upper matrices resulting from ILU factorizationCFL The Courant-Friedrichs-Lewy numberFlow Related VariablesRe Reynolds numberPr Prandtl numberPrT Turbulent Prandtl numberMa Mach numberDiff Diffusion term in the Spalart-Allmaras turbulence modelProd Production term in the Spalart-Allmaras turbulence modelDest Destruction term in the Spalart-Allmaras turbulence modelTrip Trip term in the Spalart-Allmaras turbulence model\u00CF\u0084 Viscous stress tensorxiiv Velocity vector\u00CF\u0081 DensityE Total energyP PressureT Temperature\u00C2\u00B5 Viscosity\u00C2\u00B5T Turbulent viscosity\u00CE\u00B3 Specific heat ratio\u00CE\u00BD\u00CB\u009C Spalart-Allmaras working variablet Timed Distance from wallOverscripts\u00CB\u009C(.) Roe average state\u00CB\u0086(.) Reference element statexiiiAcknowledgmentsFirst and foremost, I would like to thank my supervisor, Dr. Carl Ollivier-Gooch, for his professionalmentorship, continuous support, and endless patience. It has been a pleasure for me to work under hissupervision for the past two years.I am thankful for the financial support of the Natural Sciences and Engineering Research Council ofCanada and the University of British Columbia through the Canada Graduate Scholarships-Master\u00E2\u0080\u0099s(CGS-M) program and the Gartshore Scholarship, respectively.I would like to express my gratitude to my friends and fellow labmates for their valuable help andsuggestions. Special thanks to Alireza for patiently helping me with my initial transition into theresearch environment at ANSLab, and Gary for proofreading this thesis.I wish to thank my dear friend Mohammad. I will never forget your help and advice both inside andoutside of work. Also, thank you for your technical suggestions and for proofreading a major portion ofthis thesis.Last but not least, I would like to thank my parents, Iraj and Mahboubeh, and my sister, Leiana, for theirunconditional love, support, and caring.xivChapter 1Introduction1.1 MotivationThe advent of computational fluid dynamics (CFD) has considerably improved the design and manufac-turing processes in many industries such as aerospace, turbomachinery, oil and gas, and bioengineering.CFD finds solutions to engineering problems by numerically solving the governing and/or modelingpartial differential equations (PDEs). CFD simulations offer a cheap alternative to experimental setupsin many, though not all, situations while not suffering from the limitation of analytical methods inhandling complex geometries. Nevertheless, CFD methods are not always considered a rival to the othermethods of analysis. Computational methods are sometimes used for tuning parameters or selectingmethods of measurement in experimental setups. For example, Yoo et al. [80] employ CFD simulationsto estimate the right scale for building a nuclear fuel cask experimental model. Conversely, many PDEsthat are solved in CFD are derived from experimental data, such as the turbulence model of Spalart andAllmaras [68].A majority of CFD methods are categorized as mesh-based discretization schemes. A mesh-baseddiscretization scheme seeks to approximate the solution to a PDELu(x) = 0 inside a domain\u00E2\u0084\u00A6 \u00E2\u008A\u0082 RNdim ,where L is a differential operator, u(x) \u00E2\u0088\u0088 RNunk is the unknown exact solution, and no time dependencehas been assumed for simplicity. In the first step, the domain is subdivided (meshed) into a set ofnon-overlapping sub-volumes Th, followed by defining a discrete function in the form of uh(x;Uh).Here, the term discrete means that this function will be uniquely defined if a finite number of degreesof freedom, Uh \u00E2\u0088\u0088 RNDOF , are all specified. Furthermore, the subscript h, which is the representativelength scale of the subdivision, emphasizes the dependency of the solution on the mesh. Finally, thegoal is to find Uh, such that uh closely approximates the exact solution u as the mesh gets smaller. Therigorous definition of \u00E2\u0080\u009Cto closely approximate\u00E2\u0080\u009D is that the discretization error, eh = uh \u00E2\u0088\u0092 u, must havean asymptotic behavior such that: \u00E2\u0080\u0096eh \u00E2\u0080\u0096 = O(hp). The discretization scheme is then said to be pth-orderaccurate.1Although conventional discretization schemes for aerodynamical flows are only second-order accurate,there has been a growing interest in higher-order accurate methods. Higher-order methods can beconsiderably more efficient compared to second-order methods in achieving accurate solutions, sincethe former results in more accurate results on coarser meshes [21, 48]. Even small improvements inaccuracy can be of considerable importance. For example, Vassberg et al. [74] show that for a long-rangejet-aircraft that delivers a payload between distant city pairs, a 1% error in predicting the drag can reducethe carried payload by 7%. With the airlines operating on profit margins of only a few percent, such aloss can make the service completely unprofitable.Regarding the adoption of higher-order methods, there has always been the legitimate concern ofmodeling errors, i.e., the discrepancies between the exact solution of the governing equations and thetrue physical quantities. For a problem where the modeling errors are dominant, only a limited level ofreduction in the discretization errors is of interest. Even in such a case, preliminary results have shownthat an hp-adaptive strategy can achieve a desired level of accuracy faster than conventional second-order methods [33, 38]. An hp-adaptive method increases the local order of accuracy p and decreasesthe mesh size h at only certain regions of the domain. In the context of design and optimization,thousands of simulations might be required for optimizing a subject geometry [42]. Thus, even asmall runtime improvement for a single simulation can result in a considerable reduction of the overallresource consumption. Moreover, modeling errors are not always the dominant mode. In recent work byMavriplis [45], where he studied the numerical solution of the flow around a wing-body configurationfrom the AIAA drag prediction workshop (DPW), the dominant errors were found to be those ofdiscretization.For structured meshes, highly-accurate finite difference methods have long been developed [39, 75], andare known to have superior properties in terms of computational cost and efficiency [19]. Nonetheless,generation of structured meshes for complex geometries is a challenging task, and requires extensiveamounts of human input. Therefore, grids for complex geometries are often obtained by the fairlyautomatic and convenient alternative of unstructured mesh generation techniques.For unstructuredmeshes, various approaches have been devised for achieving higher-order accuracy, mostnotably: continuous [5] and discontinuous [29]Galerkin finite elementmethods, the correction procedurevia reconstruction formulations of the discontinuous Galerkin and spectral volume methods [31], andfinite volume schemes [32]. (See the work of Andren et al. [8] for a detailed comparison of numericalresults). The use of finite volume methods is partly motivated by their straightforward mathematicalformulation compared to the complex structure of the other mentioned methods. Furthermore, most ofthe current industrial CFD codes are based on the finite volume method. While implementing otherapproaches would require the development of completely new commercial codes, higher-order finitevolume methods can be integrated into the current industrial solvers. Finite volume methods are alsoattractive because of the smaller number of degrees of freedom that they require for the same meshcompared to finite element methods.2In two dimensions, unstructured high-order finite volume methods have been successfully applied to arange of aerodynamic problems: the Euler equations [27, 46], laminar Navier-Stokes equations [33, 40],and turbulent Reynolds Averaged Navier-Stokes (RANS) equations [32]. Although taking the effectsof turbulence into consideration is necessary for correctly capturing many aerodynamic flows, three-dimensional results are scarce and limited to the solution of Euler and laminar Navier-Stokes equations[25, 40]. In the long run, the ANSLab team at UBC is interested in developing a three-dimensionalhigher-order finite volume solver for the solution of both inviscid and viscous turbulent steady-state flowproblems, while providing approximate error bounds on target unknowns such as lift and drag. The firststep towards this goal will be taken in this thesis: solution of well-known three-dimensional benchmarkflow problems.1.2 ObjectivesThe ultimate goal of this thesis is the generalization of our current 2-D flow solver, ANSLib, for thesolution of 3-D benchmark problems involving inviscid and viscous turbulent flows. The pursuit of thisgoal is divided into the following steps:\u00E2\u0080\u00A2 Derive the finite volume formulation of the governing equations in a dimension-independentmanner. Of course, the choice of the turbulence model in this part considerably affects thesolution. In this thesis, the objective is to employ the negative Spalart-Allmaras turbulence model[6]. This one-equation model is reasonably accurate for many flow conditions, and its simplicitymakes it a good starting point for the development of numerical algorithms. The negative versionof this model is chosen because it permits negative values for the model working variable, whichcan be present when higher-order methods are employed.\u00E2\u0080\u00A2 Identify and undertake the preprocessing steps involved in handling three-dimensional grids re-quired for turbulent simulations.\u00E2\u0080\u00A2 Design a practical and parallel scalable method for the solution of the discretized system ofnonlinear equations. Since three-dimensional problems can require much more memory thantheir two-dimensional counterparts, solution methods that work in two dimensions might not befeasible in three dimensions anymore.\u00E2\u0080\u00A2 Verify the performance of the solver, and the accuracy of the solutions obtained.1.3 Thesis OutlineThis thesis is organized in the following manner:Chapter 2 provides an overview of our in-house flow solver ANSLib, in which the algorithms of thisthesis have been implemented. The key concepts of the solver, i.e., k-exact reconstruction and finite3volume discretization are discussed in detail. Also, governing and modeling equations of interestare introduced. Finally, the flux functions and the boundary conditions are revisited in a dimension-independent notation.Chapter 3 discusses the preprocessing steps of three-dimensional meshes for a finite volume simulation.First, numerical integration is addressed, where Gauss quadrature rules are employed in conjunctionwith reference element mappings to construct quadrature information for the control volumes and theirfaces. Then, the capturing of boundary curvature in highly anisotropic meshes is addressed, which isa necessity for higher-order solution methods. Finally, the issues of k-exact reconstruction in handlinghighly anisotropic meshes are reviewed and addressed.In Chapter 4, the solution of the discretized system of nonlinear equations is discussed. First, the pseudotransient continuation method is revisited, which transforms the solution of the nonlinear system ofequations into the solution of a series of linear systems. Then, a memory lean, yet effective, method forthe solution of the corresponding linear systems is proposed, and compared to the previously availablelinear solution methods in ANSLib.Chapter 5 presents the solution of four three-dimensional test problems. To test the correct imple-mentation of the mesh preprocessing algorithms, Poisson\u00E2\u0080\u0099s equation is solved in a simple geometry,where the discretization error is explicitly evaluated and its asymptotic behavior is verified. Then, thesubsonic inviscid flow around a sphere is studied, where the solution accuracy is verified by measuringthe deviation of the entropy from the inflow conditions throughout the domain. Finally, the problemsof viscous turbulent flow over a flat plate and an extruded NACA 0012 airfoil are solved, and the solu-tions are verified against the reference data provided in the NASA turbulence modeling (TMR) website[60].In the end, Chapter 6 summarizes the research of the thesis, provides conclusions, and proposes possiblefuture work.In addition, the command-line options that were provided to the ANSLib executable for running eachtest case are provided in Appendix A. A sample script for running a parallel job on the WestGrid Grexcluster [3] is also given in Appendix B, as the three-dimensional cases of this thesis were all run on thiscluster.4Chapter 2BackgroundThis chapter presents the fundamentals upon which this research has been founded. As the algorithmsinvolved in this thesis have been implemented as parts of ANSLib, it is essential to present the workingmechanism of this solver. Thus, this chapter starts with a description of the finite volume method andk-exact reconstruction, which are the numerical solution schemes in ANSLib, and then moves on to themodel equations of interest.2.1 The Finite Volume MethodThe finite volumemethod is applied to equations which can be expressed in the conservative form:\u00E2\u0088\u0082u\u00E2\u0088\u0082t+ \u00E2\u0088\u0087 \u00C2\u00B7 (F(u) \u00E2\u0088\u0092 Q(u,\u00E2\u0088\u0087u)) = S(u,\u00E2\u0088\u0087u), (2.1)where F \u00E2\u0088\u0088 RNunk\u00C3\u0097Ndim is the inviscid fluxmatrix, Q \u00E2\u0088\u0088 RNunk\u00C3\u0097Ndim is the viscous fluxmatrix, and S \u00E2\u0088\u0088 RNunkis the source term vector. In this method, the degrees of freedom are the same as the average value of thediscrete solution inside every control volume. This fundamental constraint is known as the conservationof the mean, and can be written as:1\u00E2\u0084\u00A6\u00CF\u0084\u00E2\u0088\u00AB\u00CF\u0084uh(x)d\u00E2\u0084\u00A6 = Uh,\u00CF\u0084 \u00CF\u0084 \u00E2\u0088\u0088 Th, (2.2)where \u00E2\u0084\u00A6\u00CF\u0084 and Uh,\u00CF\u0084 \u00E2\u0088\u0088 RNunk represent the volume and local DOF vector for control volume \u00CF\u0084, respec-tively. Establishing a relation between the discrete solution uh, and the degrees of freedom Uh, whensatisfying the conservation of the mean constraint, is called reconstruction in the finite volume frame-work. To accomplish this, ANSLib uses the k-exact reconstruction scheme, which will be introduced inSection 2.2.The finite volume method uses the divergence theorem to discretize Equation (2.1). Consider a control5n(uh+,\u00E2\u0088\u0087uh+)(uh-,\u00E2\u0088\u0087uh-)(uh,\u00E2\u0088\u0087uh)\u00E2\u0088\u0082\u00CF\u0084\\u00E2\u0088\u0082\u00CE\u00A9\u00E2\u0088\u0082\u00CF\u0084\u00E2\u0088\u00A9\u00E2\u0088\u0082\u00CE\u00A9Figure 2.1: Illustration of terminology in finite volume discretizationvolume \u00CF\u0084 \u00E2\u0088\u0088 Th as depicted in Figure 2.1. Integrating Equation (2.1) inside the control volume, and usingthe divergence theorem gives:dUh,\u00CF\u0084dt+1\u00E2\u0084\u00A6\u00CF\u0084\u00E2\u0088\u00AB\u00E2\u0088\u0082\u00CF\u0084\\u00E2\u0088\u0082\u00E2\u0084\u00A6(FI (u+h, u\u00E2\u0088\u0092h) \u00E2\u0088\u0092QI (u+h,\u00E2\u0088\u0087u+h, u\u00E2\u0088\u0092h,\u00E2\u0088\u0087u\u00E2\u0088\u0092h)) dS+1\u00E2\u0084\u00A6\u00CF\u0084\u00E2\u0088\u00AB\u00E2\u0088\u0082\u00CF\u0084\u00E2\u0088\u00A9\u00E2\u0088\u0082\u00E2\u0084\u00A6(FB(uh,B) \u00E2\u0088\u0092QB(uh,\u00E2\u0088\u0087uh,B)) dS \u00E2\u0088\u0092 1\u00E2\u0084\u00A6\u00CF\u0084\u00E2\u0088\u00AB\u00CF\u0084S(uh,\u00E2\u0088\u0087uh)d\u00E2\u0084\u00A6 = 0, (2.3)where B represents the boundary conditions, in the case where the control volume has faces lying on theboundary. Although the discrete solution uh and its gradient \u00E2\u0088\u0087uh are continuous inside every controlvolume, they can be discontinuous on the control volume boundaries \u00E2\u0088\u0082\u00CF\u0084. These discontinuous values areshown using the (.)+ and (.)\u00E2\u0088\u0092 notations. FI andQI represent the interior numerical flux functions, whileFB and QB represent their boundary counterparts. Numerical flux functions are designed to mimic theproduct of the original flux matrices and the normal vector while taking into account the discontinuityof the discrete solution and the boundary conditions. The flux functions used in this thesis will beintroduced in Section 2.4.Looking back at Equation (2.3), the following system of ODEs for control volume averages can bederived:dUhdt+ R(Uh) = 0, (2.4)which can be solved with a variety of time advance schemes when a time-accurate solution is of interest.Butcher [15] presents a detailed study of these methods. Although only the steady state solution is ofinterest in this thesis, the unsteady terms can be used to improve the robustness of the solver with respectto the initial solution guess. This method, known as pseudo transient continuation (PTC) [35], will beintroduced in Chapter 4.62.2 K-exact ReconstructionThe design of reconstruction schemes in the finite volume framework started with Van Leer\u00E2\u0080\u0099s MUSCLscheme [72, 73]. His scheme took advantage of the uniform pattern of control volumes in structuredmeshes andwas not directly applicable to unstructured grids. Later on, the k-exact reconstruction [12] andthe WENO/ENO [66] family of methods were designed, and did not suffer from MUSCL\u00E2\u0080\u0099s shortcomingin handling unstructured meshes. The latter, however, has poor steady state convergence properties [59].Targeted to solve aerodynamic steady state problems, ANSLib uses the k-exact reconstruction methodin its finite volume formulation.Suppose a polynomial function of order smaller or equal to k, v(x), is integrated in every control volumeto find the control volume averagesVh. If these average values are fed to a k-exact reconstruction schemeto construct the function vh(x), the identity vh(x) = v(x) must hold, which is where the name k-exactcomes from. A k-exact scheme is also called (k + 1)th-order accurate, since it results in a nominaldiscretization error of order O(h(k+1)) [13]. For simplicity, let us consider this method when there isonly one unknown variable, i.e., the vector u reduces to a scalar u. Generalization to multiple unknownvariables will then be straightforward. In this case, the solution in every control volume is defined as thesuperposition of a set of basis functions in the form:uh(x;Uh,B)|x\u00E2\u0088\u0088\u00CF\u0084 = uh,\u00CF\u0084(x;Uh,B) =Nrec\u00E2\u0088\u0091i=1ai\u00CF\u0084(Uh,B)\u00CF\u0086i\u00CF\u0084(x) \u00CF\u0084 \u00E2\u0088\u0088 Th . (2.5)Where \u00CF\u0086i\u00CF\u0084(x) and ai\u00CF\u0084 represent the ith basis function and reconstruction coefficient for control volume \u00CF\u0084,respectively. Most commonly, the basis functions will be chosen as monomials in Cartesian coordinateswith origin at the reference point of each control volume:{\u00CF\u0086i\u00CF\u0084(x)\u000C\u000C i = 1 . . . Nrec } = { 1a!b!c! (x1 \u00E2\u0088\u0092 x\u00CF\u00841)a(x2 \u00E2\u0088\u0092 x\u00CF\u00842)b(x3 \u00E2\u0088\u0092 x\u00CF\u00843)c \u000C\u000C\u000C\u000C a + b + c \u00E2\u0089\u00A4 k } , (2.6)where x\u00CF\u00841, x\u00CF\u00842, and x\u00CF\u00843 are the coordinates of the control volume\u00E2\u0080\u0099s reference location, which is usuallychosen as the centroid of volume. This particular choice of basis functions has the advantage that thediscrete solution uh will resemble the Taylor expansion of the exact solution u, which is particularlyuseful in theoretical analysis [54]. However, there are situations, e.g., highly anisotropic meshes, inwhich it would be more beneficial to use other basis functions [32]. Such a case will be introduced inChapter 3.The discrete solution must satisfy the conservation of the mean constraint, given in Equation (2.2).Moreover, for every control volume \u00CF\u0084, a specific set of its neighbors are chosen as its reconstructionstencil Stencil(\u00CF\u0084). The k-exact reconstruction requires uh,\u00CF\u0084 to predict the average values of the membersof Stencil(\u00CF\u0084) closely. Furthermore, if the control volume is located on the boundary (\u00E2\u0088\u0082\u00CF\u0084 \u00E2\u0088\u00A9 \u00E2\u0088\u0082\u00E2\u0084\u00A6 , \u00E2\u0088\u0085),enforcing uh,\u00CF\u0084 or \u00E2\u0088\u0087uh,\u00CF\u0084 \u00C2\u00B7n to have certain values at boundary quadrature points, can improve convergencein the presence of certain boundary conditions [54]. Thus, the reconstruction coefficients can be found7by solving the following constrained minimization problem:minimizea1\u00CF\u0084 ...aNrec\u00CF\u0084\u00E2\u0088\u0091\u00CF\u0083\u00E2\u0088\u0088Stencil(\u00CF\u0084)(1\u00E2\u0084\u00A6\u00CF\u0083\u00E2\u0088\u00AB\u00CF\u0083uh,\u00CF\u0084(x)d\u00E2\u0084\u00A6 \u00E2\u0088\u0092Uh,\u00CF\u0083)2+\u00E2\u0088\u0091q\u00E2\u0088\u0088BdryQuad(\u00CF\u0084)(B(uh,\u00CF\u0084(q),\u00E2\u0088\u0087uh,\u00CF\u0084(q),B))2subject to1\u00E2\u0084\u00A6\u00CF\u0084\u00E2\u0088\u00AB\u00CF\u0084uh(x)d\u00E2\u0084\u00A6 = Uh,\u00CF\u0084,(2.7)where BdryQuad(\u00CF\u0084) is the set of boundary quadrature points for control volume \u00CF\u0084, and B is the boundarycondition constraint function. The number of control volumes in Stencil(\u00CF\u0084) must be greater thanthe number of reconstruction coefficients Nrec, so that the minimization problem does not becomeundetermined. To construct Stencil(\u00CF\u0084), all the neighbors at a given topological distance from \u00CF\u0084 areadded to the stencil until the number of stencil members gets bigger than a given value MinNeigh(k).For a well behaved problem on a mesh with high regularity, choosing MinNeigh(k) = Nrec can resultin an acceptable solution. However, a value of MinNeigh(k) = 1.5Nrec is chosen to cope with meshirregularities and oscillatory solution behavior. Figure 2.2 shows a control volume and its reconstructionstencil for different k values.By introducing the integral of each basis function of every control volume inside itself and its stencilmembers:I i\u00CF\u0084\u00CF\u0083 =\u00E2\u0088\u00AB\u00CF\u0083\u00CF\u0086i\u00CF\u0084(x)d\u00E2\u0084\u00A6 \u00CF\u0083 \u00E2\u0088\u0088 Stencil(\u00CF\u0084) \u00E2\u0088\u00AA { \u00CF\u0084 } , (2.8)the constrainedminimization problem in Equations (2.7), can bewritten in a compact matrix form:\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0I1\u00CF\u0084\u00CF\u0084 . . . INrec\u00CF\u0084\u00CF\u0084I1\u00CF\u0084\u00CF\u00831 . . . INrec\u00CF\u0084\u00CF\u00831.... . ....I1\u00CF\u0084\u00CF\u0083NS(\u00CF\u0084) . . . INrec\u00CF\u0084\u00CF\u0083NS(\u00CF\u0084)\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0a1\u00CF\u0084...aNrec\u00CF\u0084\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0Uh,\u00CF\u0084Uh,\u00CF\u00831...Uh,\u00CF\u0083NS(\u00CF\u0084)\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB, (2.9)where the boundary constraints have been neglected for the sake of simplicity. The symbols \u00CF\u00831, \u00CF\u00832,. . . , \u00CF\u0083NS(\u00CF\u0084) represent the members of Stencil(\u00CF\u0084). The first row of the matrix is a constraint and mustbe exactly satisfied, while the other rows correspond to equations that have to be minimized in a leastsquares fashion. As the left hand side matrix in Equation (2.9) is only dependent on geometric terms,and does not include the average solution values Uh, the solution can be written as:\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0a1\u00CF\u0084...aNrec\u00CF\u0084\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB = A\u00E2\u0080\u00A0\u00CF\u0084\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0Uh,\u00CF\u0084Uh,\u00CF\u00831...Uh,\u00CF\u0083NS(\u00CF\u0084)\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB, (2.10)where the matrix A\u00E2\u0080\u00A0\u00CF\u0084 is the pseudo-inverse of the left hand side matrix. A change of variables inspired bythe Gaussian Elimination method is used to convert Equation (2.9) into an unconstrained optimizationproblem [54]. Then the unconstrained problem is solved by singular value decomposition (SVD) [24]8\u00CF\u0084Figure 2.2: Reconstruction stencils for a control volume \u00CF\u0084: A set of values MinNeigh(1) = 3,MinNeigh(2) = 9, and MinNeigh(3) = 18 have been used, which have resulted in reconstruc-tion stencils { blue }, { blue, magenta }, and { blue, magenta, cyan } for the 1-exact, 2-exact,and 3-exact reconstruction schemes, respectively.to yield the pseudo-inverse matrix. This process has to be evaluated for every control volume, only as apreprocessing step, which prevents it from becoming a bottleneck in our computations.In the case of equations with discontinuous solutions, shock capturing numerical schemes such as slopelimiters [76] or artificial diffusion [69] are required in conjunction with the k-exact reconstructionmethod to ensure convergence to the correct solution. For more recent implementation of these methodsin conjunction with higher-order schemes see [11, 46, 53]. In this thesis, however, the emphasis is onmodel problems without discontinuities. Thus, no shock capturing methods have been used.2.3 Studied EquationsThis section introduces the equations of interest. Namely, Poisson\u00E2\u0080\u0099s, Euler, laminar and Reynoldsaveraged Navier-Stokes, and the Spalart-Allmaras turbulence closure equations.2.3.1 Poisson\u00E2\u0080\u0099s EquationThe Poisson\u00E2\u0080\u0099s equation is a simple yet powerful tool to verify the correct implementation of manyalgorithms in ANSLib. This equation has a scalar unknown u and a solution-independent source term fin the form of:\u00E2\u0088\u0092 \u00E2\u0088\u0087 \u00C2\u00B7 (\u00E2\u0088\u0087u) = f (x). (2.11)9To simplify the notation, we consider the gradient operator on a scalar function as a row vector. Forexample, \u00E2\u0088\u0087u = [ \u00E2\u0088\u0082u\u00E2\u0088\u0082x1 , \u00E2\u0088\u0082u\u00E2\u0088\u0082x2 , \u00E2\u0088\u0082u\u00E2\u0088\u0082x3 ]. Thus the Poisson equation can be recovered from the steady state formof the conservative Equation (2.1) by replacing F = \u00E2\u0088\u0087u, Q = 0, and S = f (x).2.3.2 Navier-Stokes EquationsThe Navier-Stokes and the continuity equations are widely used for the simulation of fluid flow. In thecompressible form of these equations, the vector of unknowns is u = [\u00CF\u0081, \u00CF\u0081vT , E]T , where \u00CF\u0081 is the fluiddensity, v = [v1, v2, v3]T is the velocity vector, and E is the total energy. When combined with theideal gas internal energy and state equations, the nondimensionalized Navier-Stokes equations can beidentified by a zero source term vector and flux matrices:F =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0\u00CF\u0081vT\u00CF\u0081vvT + PI(E + P)vT\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB Q =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B00MaRe \u00CF\u0084(E + P)\u00CF\u0084v + 1\u00CE\u00B3\u00E2\u0088\u00921( \u00C2\u00B5Pr) \u00E2\u0088\u0087T\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB , (2.12)where Ma, Re, Pr, and \u00CE\u00B3 represent the Mach number, the Reynolds number, the Prandtl number, and thespecific heat ratio, respectively. T is the temperature, (.)T denotes matrix transpose, P is the pressure,\u00CF\u0084 is the viscous stress matrix, \u00C2\u00B5 is the dimensionless fluid viscosity, and I is the identity matrix. Sincethe emphasis is on air as the working fluid, the values \u00CE\u00B3 = 1.4, Pr = 0.72 are used, and \u00C2\u00B5 is found fromSutherland\u00E2\u0080\u0099s law. The pressure is related to the dependent variables via the formula:P = (\u00CE\u00B3 \u00E2\u0088\u0092 1)(E \u00E2\u0088\u0092 12\u00CF\u0081(v \u00C2\u00B7 v)). (2.13)Similarly, temperature is related to pressure and density in the form:T =\u00CE\u00B3P\u00CF\u0081. (2.14)For Newtonian compressible fluids, the viscous stress tensor is related to the velocity as:\u00CF\u0084 = 2\u00C2\u00B5(12(\u00E2\u0088\u0087v + (\u00E2\u0088\u0087v)T)\u00E2\u0088\u0092 13trace(\u00E2\u0088\u0087v)I), (2.15)where \u00E2\u0088\u0087v represents the gradient of the velocity vector, such that (\u00E2\u0088\u0087v)i j = \u00E2\u0088\u0082vi\u00E2\u0088\u0082x j .If the viscous terms are neglected, i.e., either Q, or equivalently \u00C2\u00B5 is set to zero, the Euler equationswould be recovered. The Euler equations are used in this work to asses the capability of the solver inhandling reasonably big three-dimensional problems.102.3.3 Extension to Turbulent FlowsIn this thesis, the Reynolds averaged Navier Stokes (RANS) equations are used for modeling turbulentflows, and are coupled with the negative Spalart-Allmaras (negative S-A) turbulence model [6]. Thismodel includes a number of dimensionless empirical constants, listed in Table 2.1, and a few empiricalfunctions, shown by the symbol f and an appropriate subscript, which will be introduced shortly. Thenondimensionalized flux matrices for the RANS + negative S-A equations are defined as:F =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0\u00CF\u0081vT\u00CF\u0081vvT + PI(E + P)vT\u00CE\u00BD\u00CB\u009C\u00CF\u0081vT\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BBQ =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B00MaRe \u00CF\u0084(E + P)\u00CF\u0084v + 1\u00CE\u00B3\u00E2\u0088\u00921(\u00C2\u00B5Pr +\u00C2\u00B5TPrT)\u00E2\u0088\u0087T\u00E2\u0088\u0092 MaRe\u00CF\u0083 (\u00C2\u00B5 + \u00C2\u00B5T )\u00E2\u0088\u0087\u00CE\u00BD\u00CB\u009C\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB, (2.16)where \u00CE\u00BD\u00CB\u009C is the negative S-A working variable, \u00C2\u00B5T is the turbulent viscosity, and PrT is the turbulentPrandtl number, which has a value of 0.9 for air. The source term of the nondimensionalized RANS +negative S-A is defined as:S =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0000Diff +\u00CF\u0081(Prod \u00E2\u0088\u0092Dest+Trip)\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB, (2.17)whereDiff , Prod, Dest, and Trip represent diffusion, production, destruction, and trip terms, respectively.The viscous stress tensor can be found via the modified equation:\u00CF\u0084 = 2(\u00C2\u00B5 + \u00C2\u00B5T )(12(\u00E2\u0088\u0087v + (\u00E2\u0088\u0087v)T)\u00E2\u0088\u0092 13trace(\u00E2\u0088\u0087v)I). (2.18)The turbulent viscosity is expressed as:\u00C2\u00B5T =\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3\u00C2\u00B5\u00E2\u0080\u00B2 fv1\u00CF\u0081\u00CE\u00BD\u00CB\u009C \u00CE\u00BD\u00CB\u009C \u00E2\u0089\u00A5 00 \u00CE\u00BD\u00CB\u009C < 0, (2.19)where \u00C2\u00B5\u00E2\u0080\u00B2 is the reference value that is used to nondimensionalize the turbulence working variable. In thiswork, we have used \u00C2\u00B5\u00E2\u0080\u00B2 = 1000 to make \u00CE\u00BD\u00CB\u009C comparable in size to the other nondimensionalized variables,which enhances the performance of the numerical solver [17]. The production term in Equation (2.17)is defined as:Prod =\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3cb1(1 \u00E2\u0088\u0092 ft2)S\u00CB\u009C\u00CE\u00BD\u00CB\u009C \u00CE\u00BD\u00CB\u009C \u00E2\u0089\u00A5 0cb1(1 \u00E2\u0088\u0092 ct3)S \u00CE\u00BD\u00CB\u009C < 0. (2.20)11The destruction term is found from the equation:Dest =\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3\u00C2\u00B5\u00E2\u0080\u00B2MaRe(cw1 fw \u00E2\u0088\u0092 cb1\u00CE\u00BA2 ft2) (\u00CE\u00BD\u00CB\u009Cd)2\u00CE\u00BD\u00CB\u009C \u00E2\u0089\u00A5 0\u00E2\u0088\u0092\u00C2\u00B5\u00E2\u0080\u00B2MaRe cw1(\u00CE\u00BD\u00CB\u009Cd)2\u00CE\u00BD\u00CB\u009C < 0, (2.21)where d is the minimum distance to wall boundaries. The diffusion term is given as:Diff =MaRe\u00CF\u0083(\u00C2\u00B5\u00E2\u0080\u00B2cb2\u00CF\u0081\u00E2\u0088\u0087\u00CE\u00BD\u00CB\u009C \u00C2\u00B7 \u00E2\u0088\u0087\u00CE\u00BD\u00CB\u009C \u00E2\u0088\u0092 \u00C2\u00B5\u00CF\u0081(1 + \u00CF\u0087 fn) \u00E2\u0088\u0087\u00CF\u0081 \u00C2\u00B7 \u00E2\u0088\u0087\u00CE\u00BD\u00CB\u009C), (2.22)where \u00CF\u0087 is:\u00CF\u0087 =\u00C2\u00B5\u00E2\u0080\u00B2\u00CF\u0081\u00CE\u00BD\u00CB\u009C\u00C2\u00B5. (2.23)The trip term inEquation (2.17)models flows that include the laminar to turbulent transition phenomenon.Since the emphasis in this thesis is on fully turbulent flows, the trip term is neglected and set to zero.The vorticity S, and its modified forms, S\u00CB\u009C, S\u00C2\u00AF, are found from the equations:S =\u00E2\u0088\u009A\u00E2\u0088\u009A12\u00E2\u0088\u0091i\u00E2\u0088\u0091j(\u00E2\u0088\u0082vi\u00E2\u0088\u0082xj\u00E2\u0088\u0092 \u00E2\u0088\u0082vj\u00E2\u0088\u0082xi)2S\u00C2\u00AF =\u00C2\u00B5\u00E2\u0080\u00B2MaRe\u00CE\u00BD\u00CB\u009C\u00CE\u00BA2d2fv2S\u00CB\u009C =\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3S + S\u00C2\u00AF S\u00C2\u00AF \u00E2\u0089\u00A5 \u00E2\u0088\u0092cv2SS +c2v2S+cv3S\u00C2\u00AF(cv3\u00E2\u0088\u00922cv2)S\u00E2\u0088\u0092S\u00C2\u00AF S\u00C2\u00AF < \u00E2\u0088\u0092cv2S.(2.24)The empirical function fn ensures the positivity of \u00C2\u00B5T , and is defined as:fn =\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B31 \u00CE\u00BD\u00CB\u009C \u00E2\u0089\u00A5 0cn1+\u00CF\u00873cn1\u00E2\u0088\u0092\u00CF\u00873 \u00CE\u00BD\u00CB\u009C < 0. (2.25)The functions fv1, fv2, fv3 are given as:fv1 =\u00CF\u00873\u00CF\u00873 + c3v1fv2 = 1 \u00E2\u0088\u0092 \u00CF\u00871 + \u00CF\u0087 fv1 ft2 = ct3 exp(\u00E2\u0088\u0092ct4\u00CF\u00872), (2.26)and the function fw is given as:r =\u00C2\u00B5\u00E2\u0080\u00B2MaRe\u00CE\u00BD\u00CB\u009CS\u00CB\u009C\u00CE\u00BA2d2g = r + cw2(r6 \u00E2\u0088\u0092 r) fw = g(1 + c6w3g6 + c6w3). (2.27)Having defined the equations of interest, let us now introduce the numerical flux functions.12Table 2.1: Dimensionless empirical parameters used in the negative S-A modelName Value Name Value Name Valuecb1 0.1355 cb2 0.622 cw1 cb1\u00CE\u00BA2 +1+cb2\u00CF\u0083cw2 0.3 cw3 2.0 ct3 1.2ct4 0.5 cn1 16 cv1 7.1cv2 0.7 cv3 0.9 \u00CE\u00BA 0.41\u00CF\u0083 0.662.4 Numerical Flux FunctionsAs discussed earlier, the finite volume method relies on flux functions that not only take into accountthe discontinuity of the solution along internal faces, but also correctly capture the influence of theboundary conditions. In this section, we will introduce the numerical flux functions used in ANSLibfor the solution of the more general RANS + negative S-A equations. Flux functions for the simplerEuler and Poisson equations can simply be derived by omitting the relevant terms from the more generalcase.2.4.1 Inviscid FluxThe inviscid fluxes represent propagation of information via finite speed waves, and are convective in na-ture. Most numerical flux functions seek a solution to the approximate one-dimensional equation:\u00E2\u0088\u0082u\u00E2\u0088\u0082t+\u00E2\u0088\u0082F(u)n\u00E2\u0088\u0082s= 0, (2.28)and then use this solution to find the numerical flux. In Equation (2.28), F is the inviscid flux matrix,n is the face normal vector, and s is the unit of length. The artificial viscosity method [79] adds extranonlinear dissipative terms to Equation (2.28) to give it an elliptic nature. On the other hand, theGodunov method [23] is based on the exact solution of Equation (2.28), with piecewise constant initialconditions corresponding to the left and right states (the Riemann Problem). Due to the expensive cost ofexactly solving the Riemann problem, many researchers have developed approximate solutions such asthe Rusanov [61], HLL family [28, 70], and Roe [58] methods. In this thesis, a computationally efficientformulation [57] of the approximate Riemann solver of Roe has been used due to its effectiveness andsimplicity. This formulation has the form ofFI (u+h, u\u00E2\u0088\u0092h) =12(F(u+h)n + F(u\u00E2\u0088\u0092h)n \u00E2\u0088\u0092 D(u+h, u\u00E2\u0088\u0092h)), (2.29)13where D represents the diffusion vector given by the formula:D(u+h, u\u00E2\u0088\u0092h) =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0|\u00CE\u00BB\u00CB\u009C2 |\u00E2\u0088\u0086\u00CF\u0081 + \u00CE\u00B41|\u00CE\u00BB\u00CB\u009C2 |\u00E2\u0088\u0086(\u00CF\u0081v) + \u00CE\u00B41v\u00CB\u009C + \u00CE\u00B42n|\u00CE\u00BB\u00CB\u009C2 |\u00E2\u0088\u0086E + \u00CE\u00B41H\u00CB\u009C + \u00CE\u00B42v\u00CB\u009C \u00C2\u00B7 n|\u00CE\u00BB\u00CB\u009C2 |\u00E2\u0088\u0086(\u00CF\u0081\u00CE\u00BD\u00CB\u009C) + \u00CE\u00B41 \u00CB\u009C\u00CB\u009C\u00CE\u00BD\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB. (2.30)Here, \u00E2\u0088\u0086(.) = (.)+h\u00E2\u0088\u0092 (.)\u00E2\u0088\u0092h, and the variables \u00CE\u00BB1, . . . , \u00CE\u00BB6 represent the six eigenvalues of the Jacobian matrix\u00E2\u0088\u0082F(u)n\u00E2\u0088\u0082u , expressed as:\u00CE\u00BB1 = v \u00C2\u00B7 n \u00E2\u0088\u0092 c \u00CE\u00BB2 = . . . = \u00CE\u00BB5 = v \u00C2\u00B7 n \u00CE\u00BB6 = v \u00C2\u00B7 n + c. (2.31)The variable c is the speed of sound, defined as:c =\u00E2\u0088\u009A\u00CE\u00B3P\u00CF\u0081. (2.32)In Equation (2.30) the symbol \u00CB\u009C(.) represents the Roe average state, evaluated from the equation:\u00CB\u009C(.) =\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3\u00E2\u0088\u009A(.)\u00E2\u0088\u0092h(.)+h(.) = \u00CF\u0081\u00E2\u0088\u009A\u00CF\u0081\u00E2\u0088\u0092h(.)\u00E2\u0088\u0092h+\u00E2\u0088\u009A\u00CF\u0081+h(.)+h\u00E2\u0088\u009A\u00CF\u0081\u00E2\u0088\u0092h+\u00E2\u0088\u009A\u00CF\u0081+h(.) = v, \u00CE\u00BD\u00CB\u009C,H, (2.33)where H = P+E\u00CF\u0081 is the enthalpy. Note that v\u00CB\u009C and \u00CB\u009C\u00CB\u009C\u00CE\u00BD are the velocity and the S-A working variable at theRoe average state, respectively. The variables \u00CE\u00B41 and \u00CE\u00B42 in Equation (2.30) are given as:\u00CE\u00B41 =\u00E2\u0088\u0086Pc\u00CB\u009C2(\u00E2\u0088\u0092|\u00CE\u00BB\u00CB\u009C2 | + |\u00CE\u00BB\u00CB\u009C1 | + |\u00CE\u00BB\u00CB\u009C6 |2)+\u00CF\u0081\u00CB\u009C2c\u00CB\u009C2\u00E2\u0088\u0086(v \u00C2\u00B7 n) ( |\u00CE\u00BB\u00CB\u009C6 | \u00E2\u0088\u0092 |\u00CE\u00BB\u00CB\u009C1 |)\u00CE\u00B42 =\u00E2\u0088\u0086P2c\u00CB\u009C2( |\u00CE\u00BB\u00CB\u009C6 | \u00E2\u0088\u0092 |\u00CE\u00BB\u00CB\u009C1 |) + \u00CF\u0081\u00CB\u009C\u00E2\u0088\u0086(v \u00C2\u00B7 n) (\u00E2\u0088\u0092|\u00CE\u00BB\u00CB\u009C2 | + |\u00CE\u00BB\u00CB\u009C1 | + |\u00CE\u00BB\u00CB\u009C6 |2 ) . (2.34)ANSLib evaluates the the inviscid boundary fluxes using the method of characteristics [30]. For walland symmetry boundary conditions, the mass flux must be zero. Thus, the inviscid flux is evaluatedas:FB(uh,B) = [0, PhnT , 0]T B is symmetry or wall. (2.35)On farfield boundaries, however, the boundary flux formulation changes based on the sign of the normalvelocity, v \u00C2\u00B7 n. In this work, we consider subsonic inflow and outflow boundary conditions which areidentified by the inequalities 0 \u00E2\u0089\u00A4 v \u00C2\u00B7 n < c, and \u00E2\u0088\u0092c < v \u00C2\u00B7 n \u00E2\u0089\u00A4 0, respectively. In either case, the boundaryflux is evaluated in terms of an intermediate state, u\u00E2\u0088\u0097h, defined as a function of both the boundaryconditions and the discrete solution:FB(uh,B) = F(u\u00E2\u0088\u0097h)n u\u00E2\u0088\u0097h = (uh,B) B is inflow or outflow. (2.36)14For subsonic inflow, five values of farfield turbulence working variable \u00CE\u00BD\u00CB\u009Cfar, total temperature Tt , totalpressure Pt , side slip angle \u00CF\u0088, and angle of attack \u00CE\u00B1 must be specified. Subsequently, the intermediatestate u\u00E2\u0088\u0097hcan be found as:P\u00E2\u0088\u0097h = Ph T\u00E2\u0088\u0097h = Tt(P\u00E2\u0088\u0097hPt) \u00CE\u00B3\u00E2\u0088\u00921\u00CE\u00B3\u00CF\u0081\u00E2\u0088\u0097h = \u00CE\u00B3P\u00E2\u0088\u0097hT\u00E2\u0088\u0097h\u00CE\u00BD\u00CB\u009C\u00E2\u0088\u0097h = \u00CE\u00BD\u00CB\u009Cfar\u00E2\u0080\u0096v\u00E2\u0088\u0097h \u00E2\u0080\u00962 =\u00E2\u0088\u009A2\u00CE\u00B3 \u00E2\u0088\u0092 1(TtT\u00E2\u0088\u0097h\u00E2\u0088\u0092 1)v\u00E2\u0088\u0097h1 = \u00E2\u0080\u0096v\u00E2\u0088\u0097h \u00E2\u0080\u00962 cos\u00CE\u00B1 cos\u00CF\u0088 v\u00E2\u0088\u0097h2 = \u00E2\u0080\u0096v\u00E2\u0088\u0097h \u00E2\u0080\u00962 sin\u00CE\u00B1 cos\u00CF\u0088 v\u00E2\u0088\u0097h3 = \u00E2\u0080\u0096v\u00E2\u0088\u0097h \u00E2\u0080\u00962 sin\u00CF\u0088.(2.37)For subsonic outflow, only the back-pressure value Pb must be specified. The intermediate state willthen be defined as:\u00CF\u0081\u00E2\u0088\u0097h = \u00CF\u0081h v\u00E2\u0088\u0097h = vh P\u00E2\u0088\u0097h = Pb \u00CE\u00BD\u00CB\u009C\u00E2\u0088\u0097h = \u00CE\u00BD\u00CB\u009Ch . (2.38)In this thesis, the same values are chosen for \u00CE\u00BD\u00CB\u009Cfar, Tt , Pt , and Pb as the work of Jalali and Ollivier-Gooch[32]:\u00CE\u00BD\u00CB\u009Cfar =3\u00C2\u00B5\u00E2\u0080\u00B2Tt = 1 +\u00CE\u00B3 \u00E2\u0088\u0092 12Ma2 Pt =1\u00CE\u00B3T\u00CE\u00B3\u00CE\u00B3\u00E2\u0088\u00921t Pb =1\u00CE\u00B3. (2.39)2.4.2 Viscous FluxThe viscous flux functions are evaluated in a different manner compared to their inviscid counterparts.When evaluating internal flux values, an intermediate state u\u00E2\u0088\u0097his used in the form:QI (u+h,\u00E2\u0088\u0087u+h, u\u00E2\u0088\u0092h,\u00E2\u0088\u0087u\u00E2\u0088\u0092h) = Q(u\u00E2\u0088\u0097h,\u00E2\u0088\u0087u\u00E2\u0088\u0097h)n. (2.40)Finding the intermediate solution as the numerical average of the left and right states,u\u00E2\u0088\u0097h =12(u+h + u\u00E2\u0088\u0092h), (2.41)results in a sufficiently accurate solution [51]. When evaluating the intermediate solution gradient,however, simple averaging may lead to spurious solutions and instabilities. Nishikawa [50] suggestedthe following modified formula, in line with the interior penalty formulation [9] used in the DG frame-work:\u00E2\u0088\u0087u\u00E2\u0088\u0097h =12(\u00E2\u0088\u0087u+h + \u00E2\u0088\u0087u\u00E2\u0088\u0092h ) + \u00CE\u00B7 ( u+h \u00E2\u0088\u0092 u\u00E2\u0088\u0092h\u00E2\u0080\u0096x\u00CF\u0084+ \u00E2\u0088\u0092 x\u00CF\u0084\u00E2\u0088\u0092\u00E2\u0080\u00962)n, (2.42)where x\u00CF\u0084\u00E2\u0088\u0092 and x\u00CF\u0084+ are the reference locations of the adjacent control volumes, respectively, and \u00CE\u00B7 is aheuristic damping factor, known as the jump term. Jalali et al. [34] have numerically tested differentvalues of \u00CE\u00B7 in high-order finite volume simulations, and have recommended a value of \u00CE\u00B7 = 1, which isalso used in this thesis.15The boundary viscous fluxes are evaluated in the same manner as reference [32], by simply replacingthe interior state into the viscous flux matrix,QB(uh,\u00E2\u0088\u0087uh,B) = Q(uh,\u00E2\u0088\u0087uh)n. The boundary conditionsare then enforced through a combination of soft and hard constraints. For adiabatic walls, the velocity,heat flux, and turbulence working variable have to be zero, which are expressed as:Hard constraints: vh = 0 \u00CE\u00BD\u00CB\u009Ch = 0Soft constraint: \u00E2\u0088\u0087Th \u00C2\u00B7 n = 0B is adiabatic wall, (2.43)where the hard constraints are applied to the k-exact reconstruction process through the boundarycondition constraint function B. The soft constraint, however, is applied by simply replacing the valueof \u00E2\u0088\u0087Th \u00C2\u00B7 n with zero when evaluating the boundary flux. On symmetry boundaries, the heat flux, normalderivative of the turbulence working variable, and the tangential viscous force have to be zero, which areenforced via the following soft constraints:Soft constraints: \u00E2\u0088\u0087\u00CE\u00BD\u00CB\u009Ch \u00C2\u00B7 n = 0 \u00E2\u0088\u0087Th \u00C2\u00B7 n = 0((nT\u00CF\u0084hn)I \u00E2\u0088\u0092 \u00CF\u0084h)n = 0 B is symmetry. (2.44)Finally, no constraints are applied when evaluating the viscous flux function on the farfield bound-aries.16Chapter 3Three-Dimensional Mesh ProcessingThe finite volume method requires the subdivision of the domain into a set of control volumes, whichare obtained from a mesh. There are two different approaches for constructing the control volumes:cell-centered, and vertex-centered (cell-vertex). In the cell-centered approach, every cell of the meshis considered as a control volume, as shown in Figure 3.1a. Conversely, the vertex-centered approachassociates a control volume to each vertex of the mesh, and is dependent on the definition of the controlvolume faces. Figure 3.1b shows an example of vertex-centered control volumes, which are created byconnecting the barycenter of each triangle to the midpoints of its edges. This method for constructingthe vertex-centered control volumes is known as the median-dual approach.The only data that the solver requires about the geometry of the control volumes is their adjacencyand quadrature information. Thus, relevant data structures and algorithms have to be implementedin a numerical solver package to offer access to such information, with reasonable time and memorycost. Although the cell-centered and cell-vertex methods are both straightforward to implement in twodimensions, three-dimensional implementation of the latter is more difficult to code, and requires morequadrature points for integration at a given order of accuracy. As a result, the cell-centered method hasbeen chosen for the purpose of this thesis, and is implemented in ANSLib. This chapter discusses thepreprocessing steps required for a three-dimensional cell-centered finite volume solver.3.1 Element Mapping and QuadratureIn terms of solution accuracy, hexahedral and quadrilateral meshes are known to have superior propertiesin three and two dimensions, respectively, compared to tetrahedral and triangular meshes. Althoughstructured mesh generation techniques can create such grids, they can be expensive in terms of timeresources, since applying these algorithms to complex geometries requires a considerable amount ofhuman input. Conversely, unstructured mesh generation techniques are rather automatic, but are notable to handle complex geometries without resorting to robust algorithms that only create simplex cells,17(a) (b)Figure 3.1: Schematic of control volumes: (a) cell-centered, (b) median-dual cell-vertex.such as the Delaunay mesh generation methods [18]. Therefore, simulation grids are usually composedof a mixture of different cell types. The purpose of this section is to introduce the connectivity andquadrature formulas for the control volumes and their faces in a computational mesh. Hereinafter, boththe control volumes and their faces will be referred to as elements.The process starts with mapping each element E to its reference version E\u00CB\u0086 , which resides in a nondimen-sional space. The mapping is a polynomial of degree l, that must guarantee geometric continuity acrossneighboring elements. Thus, it is usually constructed from Lagrangian basis functions (polynomials).Using the \u00CB\u0086(.) notation for the reference space, the mapping can be written as:x(x\u00CB\u0086) =Nlag\u00E2\u0088\u0091i=1yi\u00CF\u0088\u00CB\u0086i(x\u00CB\u0086) x\u00CB\u0086 \u00E2\u0088\u0088 E\u00CB\u0086, (3.1)where \u00CF\u0088\u00CB\u0086i and yi are the Lagrange basis functions and nodal locations, respectively, and Nlag is the numberof Lagrange points. Hexahedra, prisms and tetrahedra are the most commonly used three-dimensionalelements in both finite element and finite volume simulations. Pyramids, on the other hand, are notvery popular, and are used for making transitions between different element types in a mixed mesh.Quadrilateral and triangular elements are also used to represent the control volume faces. Figure 3.2shows the first order (l = 1) versions of these elements in the reference space. A well defined mesh mustprovide the location of every node. Furthermore, it must store the local to global node numbering forevery element, along with its type and interpolation order.For a reference element, a qth-order accurate quadrature rule, consists of Nqua point locations z\u00CB\u0086i, andweights w\u00CB\u0086i, such that for every polynomial P\u00CB\u0086(x\u00CB\u0086) of degree less than or equal to q \u00E2\u0088\u0092 1:\u00E2\u0088\u00ABE\u00CB\u0086P\u00CB\u0086(x\u00CB\u0086)dE\u00CB\u0086 =Nqua\u00E2\u0088\u0091i=1P\u00CB\u0086(z\u00CB\u0086i)w\u00CB\u0086i . (3.2)These quadrature rules are tabulated for different orders and types of elements in the literature, such asthe work of Solin et al. [67]. By a change of variables, quadrature rules can be constructed for elements18x\u00CB\u00861x\u00CB\u00862x\u00CB\u0086312 3456 78 y\u00CB\u00861 = (\u00E2\u0088\u00921, \u00E2\u0088\u00921, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00861 = (1 \u00E2\u0088\u0092 x\u00CB\u00861)(1 \u00E2\u0088\u0092 x\u00CB\u00862)(1 \u00E2\u0088\u0092 x\u00CB\u00863)/8y\u00CB\u00862 = (+1, \u00E2\u0088\u00921, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00862 = (1 + x\u00CB\u00861)(1 \u00E2\u0088\u0092 x\u00CB\u00862)(1 \u00E2\u0088\u0092 x\u00CB\u00863)/8y\u00CB\u00863 = (+1, +1, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00863 = (1 + x\u00CB\u00861)(1 + x\u00CB\u00862)(1 \u00E2\u0088\u0092 x\u00CB\u00863)/8y\u00CB\u00864 = (\u00E2\u0088\u00921, +1, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00864 = (1 \u00E2\u0088\u0092 x\u00CB\u00861)(1 + x\u00CB\u00862)(1 \u00E2\u0088\u0092 x\u00CB\u00863)/8y\u00CB\u00865 = (\u00E2\u0088\u00921, \u00E2\u0088\u00921, +1) \u00CF\u0088\u00CB\u00865 = (1 \u00E2\u0088\u0092 x\u00CB\u00861)(1 \u00E2\u0088\u0092 x\u00CB\u00862)(1 + x\u00CB\u00863)/8y\u00CB\u00866 = (+1, \u00E2\u0088\u00921, +1) \u00CF\u0088\u00CB\u00866 = (1 + x\u00CB\u00861)(1 \u00E2\u0088\u0092 x\u00CB\u00862)(1 + x\u00CB\u00863)/8y\u00CB\u00867 = (+1, +1, +1) \u00CF\u0088\u00CB\u00867 = (1 + x\u00CB\u00861)(1 + x\u00CB\u00862)(1 + x\u00CB\u00863)/8y\u00CB\u00868 = (\u00E2\u0088\u00921, +1, +1) \u00CF\u0088\u00CB\u00868 = (1 \u00E2\u0088\u0092 x\u00CB\u00861)(1 + x\u00CB\u00862)(1 + x\u00CB\u00863)/8123456y\u00CB\u00861 = (0, 0, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00861 = (1 \u00E2\u0088\u0092 x\u00CB\u00863)(1 \u00E2\u0088\u0092 x\u00CB\u00862 \u00E2\u0088\u0092 x\u00CB\u00861)/2y\u00CB\u00862 = (1, 0, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00862 = (1 \u00E2\u0088\u0092 x\u00CB\u00863)x\u00CB\u00861/2y\u00CB\u00863 = (0, 1, \u00E2\u0088\u00921) \u00CF\u0088\u00CB\u00863 = (1 \u00E2\u0088\u0092 x\u00CB\u00863)x\u00CB\u00862/2y\u00CB\u00864 = (0, 0, +1) \u00CF\u0088\u00CB\u00864 = (1 + x\u00CB\u00863)(1 \u00E2\u0088\u0092 x\u00CB\u00862 \u00E2\u0088\u0092 x\u00CB\u00861)/2y\u00CB\u00865 = (1, 0, +1) \u00CF\u0088\u00CB\u00865 = (1 + x\u00CB\u00863)x\u00CB\u00861/2y\u00CB\u00866 = (0, 1, +1) \u00CF\u0088\u00CB\u00866 = (1 + x\u00CB\u00863)x\u00CB\u00862/212 345y\u00CB\u00861 = (\u00E2\u0088\u00921, \u00E2\u0088\u00921, 0) \u00CF\u0088\u00CB\u00861 = (x\u00CB\u00863 + x\u00CB\u00861 \u00E2\u0088\u0092 1)(x\u00CB\u00863 + x\u00CB\u00862 \u00E2\u0088\u0092 1)/4((1 \u00E2\u0088\u0092 x\u00CB\u00863) + \u000F )y\u00CB\u00862 = (+1, \u00E2\u0088\u00921, 0) \u00CF\u0088\u00CB\u00862 = (x\u00CB\u00863 \u00E2\u0088\u0092 x\u00CB\u00861 \u00E2\u0088\u0092 1)(x\u00CB\u00863 + x\u00CB\u00862 \u00E2\u0088\u0092 1)/4((1 \u00E2\u0088\u0092 x\u00CB\u00863) + \u000F )y\u00CB\u00863 = (+1, \u00E2\u0088\u00921, 0) \u00CF\u0088\u00CB\u00863 = (x\u00CB\u00863 \u00E2\u0088\u0092 x\u00CB\u00861 \u00E2\u0088\u0092 1)(x\u00CB\u00863 \u00E2\u0088\u0092 x\u00CB\u00862 \u00E2\u0088\u0092 1)/4((1 \u00E2\u0088\u0092 x\u00CB\u00863) + \u000F )y\u00CB\u00864 = (+1, \u00E2\u0088\u00921, 0) \u00CF\u0088\u00CB\u00864 = (x\u00CB\u00863 + x\u00CB\u00861 \u00E2\u0088\u0092 1)(x\u00CB\u00863 \u00E2\u0088\u0092 x\u00CB\u00862 \u00E2\u0088\u0092 1)/4((1 \u00E2\u0088\u0092 x\u00CB\u00863) + \u000F )y\u00CB\u00865 = (0, 0, 1) \u00CF\u0088\u00CB\u00865 = x\u00CB\u00863 \u000F = 10\u00E2\u0088\u0092201234y\u00CB\u00861 = (0, 0, 0) \u00CF\u0088\u00CB\u00861 = 1 \u00E2\u0088\u0092 x\u00CB\u00861 \u00E2\u0088\u0092 x\u00CB\u00862 \u00E2\u0088\u0092 x\u00CB\u00863y\u00CB\u00862 = (1, 0, 0) \u00CF\u0088\u00CB\u00862 = x\u00CB\u00861y\u00CB\u00863 = (0, 1, 0) \u00CF\u0088\u00CB\u00863 = x\u00CB\u00862y\u00CB\u00864 = (0, 0, 1) \u00CF\u0088\u00CB\u00864 = x\u00CB\u00863x\u00CB\u00861x\u00CB\u00862x\u00CB\u0086312 34 y\u00CB\u00861 = (\u00E2\u0088\u00921, \u00E2\u0088\u00921, 0) \u00CF\u0088\u00CB\u00861 = (1 \u00E2\u0088\u0092 x\u00CB\u00861)(1 \u00E2\u0088\u0092 x\u00CB\u00862)/4y\u00CB\u00862 = (+1, \u00E2\u0088\u00921, 0) \u00CF\u0088\u00CB\u00862 = (1 + x\u00CB\u00861)(1 \u00E2\u0088\u0092 x\u00CB\u00862)/4y\u00CB\u00863 = (+1, +1, 0) \u00CF\u0088\u00CB\u00863 = (1 + x\u00CB\u00861)(1 + x\u00CB\u00862)/4y\u00CB\u00864 = (\u00E2\u0088\u00921, +1, 0) \u00CF\u0088\u00CB\u00864 = (1 \u00E2\u0088\u0092 x\u00CB\u00861)(1 + x\u00CB\u00862)/4123 y\u00CB\u00861 = (0, 0, 0) \u00CF\u0088\u00CB\u00861 = 1 \u00E2\u0088\u0092 x\u00CB\u00861 \u00E2\u0088\u0092 x\u00CB\u00862y\u00CB\u00862 = (1, 0, 0) \u00CF\u0088\u00CB\u00862 = x\u00CB\u00861y\u00CB\u00863 = (0, 1, 0) \u00CF\u0088\u00CB\u00863 = x\u00CB\u00862Figure 3.2: First order reference Lagrange elements and polynomials. From top to bottom: hexa-hedron, prism, pyramid, tetrahedron, quadrilateral, and triangle.19in the physical space. For a three-dimensional element, we can write:\u00E2\u0088\u00ABEf (x)dE =\u00E2\u0088\u00ABE\u00CB\u0086f (x(x\u00CB\u0086)) det(J(x\u00CB\u0086))dE\u00CB\u0086 \u00E2\u0089\u0088Nqua\u00E2\u0088\u0091i=1f (x(z\u00CB\u0086i)) det(J(z\u00CB\u0086i))w\u00CB\u0086i . (3.3)Here, J = \u00E2\u0088\u0082x\u00E2\u0088\u0082x\u00CB\u0086 represents the Jacobian of the transformation given in Equation (3.1). A similar statementcan be made for two-dimensional elements, giving way to the following procedure for finding quadraturerules in the physical space:2D Element: wi = \u00E2\u0080\u0096\u00E2\u0088\u0082x\u00CB\u00861x \u00C3\u0097 \u00E2\u0088\u0082x\u00CB\u00862x\u00E2\u0080\u00962w\u00CB\u0086i zi = x(z\u00CB\u0086i)3D Element: wi = det(J(z\u00CB\u0086i))w\u00CB\u0086i zi = x(z\u00CB\u0086i). (3.4)It is noteworthy to mention that the higher the interpolation order of an element, l, the higher the valuesof q are required to evaluate an integral up to a certain order of accuracy. This is due to the presenceof the term det(J) in Equation (3.3) that prevents a quadrature rule to maintain its nominal order ofaccuracy in the physical space. Furthermore, for integrals on the faces that involve the normal vector,the following formula can be used:n = \u00E2\u0088\u0082x\u00CB\u00861x \u00C3\u0097 \u00E2\u0088\u0082x\u00CB\u00862x\u00E2\u0080\u0096\u00E2\u0088\u0082x\u00CB\u00861x \u00C3\u0097 \u00E2\u0088\u0082x\u00CB\u00862x\u00E2\u0080\u00962. (3.5)In this thesis we have used the Lagrange polynomials and reference element quadrature rules that areprovided by the libMesh finite element library [36].3.2 Creating Curved Anisotropic MeshesHigh-order numerical discretization schemes must correctly account for the curvature of domain bound-aries to maintain their nominal order of accuracy [14]. As conventional mesh generators create onlymeshes with planar faces, modification schemes are required before such meshes can be used for ahigh-order numerical simulation. This modification can be as simple as just replacing the boundaryfaces with higher-order representations, as in isotropic meshes. Successful resolution of quantities ofinterest in turbulent flows, however, requires highly anisotropic meshes that have big spacing in the walltangent direction compared to small spacing in the wall normal direction. Replacing only the boundaryfaces in this case would create self intersecting invalid domain subdivisions. Thus, the curvature of theboundary faces must somehow be propagated throughout the domain to prevent mesh tangling.For structuredmeshes, a simple remedy is to agglomerate a finemesh and use the extra points to define thehigher order representation of the cells and their faces. Although this method is effective, its limitationin only working with structured meshes confines its use to verification and benchmark problems. A moregeneral approach is to use a solid mechanics analogy, and model the faceted mesh as an elastic solid.Then, the boundary of the solid can be deformed to match its curved representation, leading to an internaldeformation which is consequently used to define the curvature of the inner faces and cells. Hartmann20and Leicht [29] gave a detailed summary of this family of methods. They also introduced modificationsthat greatly reduce the computational cost, while increasing the quality of the final mesh produced.Another approach was presented by Toulorge et al. [71], which is based on minimizing a novel nonlinearenergy function. Although this approach is more complicated than using a linear elasticity analogy, theauthors argue that it has better properties in terms of robustness and efficiency. Jalali and Ollivier-Gooch[32] have developed a simplified two-dimensional linear elasticity solver, which is currently used inANSLib to generate high-order curved meshes for two-dimensional turbulent flow simulations. In thissection, their approach will be extended to work with three-dimensional unstructured meshes.A linear elastic solid is governed by the Navier equations, which are derived from the conservation oflinear momentum, and expressed as [37]:\u00E2\u0088\u0092 \u00E2\u0088\u0087 \u00C2\u00B7 \u00CF\u0083 = 0, (3.6)where \u00CF\u0083 is the stress tensor, and is defined as:\u00CF\u0083 =\u00C2\u00B52(\u00E2\u0088\u0087u + \u00E2\u0088\u0087uT)\u00E2\u0088\u0092 \u00CE\u00BB3trace(\u00E2\u0088\u0087u)I, (3.7)where u = [u1, u2, u3]T is the displacement vector, and the variables \u00C2\u00B5 and \u00CE\u00BB are the Lame\u00E2\u0080\u0099s coefficients.The elastic solid properties are usually reported in terms of the Young\u00E2\u0080\u0099s modulus, E , and Poisson\u00E2\u0080\u0099s ratio,\u00CE\u00BD, which are related to the Lame\u00E2\u0080\u0099s coefficients as:\u00C2\u00B5 =E2(1 + \u00CE\u00BD) \u00CE\u00BB =\u00CE\u00BDE(1 + \u00CE\u00BD)(1 \u00E2\u0088\u0092 2\u00CE\u00BD) . (3.8)On the boundaries, either the displacement, the normal force, or a linear relation between them has tobe specified, which results in Dirichlet, Neumann, or Robin boundary conditions, respectively. Theseconditions can all be expressed through a single equation:Q(x)u + \u00CF\u0083n = f(x) x \u00E2\u0088\u0088 \u00E2\u0088\u0082\u00E2\u0084\u00A6, (3.9)where f is the normal force, and Q is a generalized spring constant. In the curved mesh generationprocess, only the Dirichlet boundary conditions are of interest, which can be recovered by settingQ = 1\u000F I and f(x) = 1\u000F ub(x) in Equation (3.9). Here, ub(x) is the prescribed boundary displacement, and\u000F is a small number such as 10\u00E2\u0088\u009210.The solution method for the Navier equations and their boundary conditions will now be presented. Toease the notation, the Einstein summation convention is used for indices i, j, k, and l, which are reservedonly for the geometric dimensions. In the first solution step, Equations (3.6) and (3.7) are combined toyield:\u00E2\u0088\u0092 \u00E2\u0088\u0082\u00E2\u0088\u0082xk(Ki jkl\u00E2\u0088\u0082u j\u00E2\u0088\u0082xl)= 0 i = 1 \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 Ndim. (3.10)21where Ki jkl are the components of the fourth-order stiffness tensor:Ki jkl = \u00CE\u00BB\u00CE\u00B4ik\u00CE\u00B4jl + \u00C2\u00B5(\u00CE\u00B4i j\u00CE\u00B4kl + \u00CE\u00B4il\u00CE\u00B4jk). (3.11)The next step involves the use of the continuous Galerkin finite element method to solve the governingequations. The numerical solution uh(x) will be defined as the linear superposition of a set of basisfunctions,uh(x) =NDOF/Ndim\u00E2\u0088\u0091r=1Uh,r\u00CF\u0088r (x), (3.12)where \u00CF\u0088r is the rth finite element basis function. Note that the finite element basis functions arenot necessarily the same as Lagrange polynomials that are used to represent the curved geometry inEquation (3.1). Moreover, Uh,r is the rth sub-component of the vector of degrees of freedom, Uh, andis expressed as:Uh,r = [(Uh,r )1, (Uh,r )2, (Uh,r )3]T . (3.13)Subsequently, u is replaced with uh in Equation (3.10), and the result is multiplied by every basis function\u00CF\u0088r , which when integrated over the whole domain gives:\u00E2\u0088\u0092\u00E2\u0088\u00AB\u00E2\u0084\u00A6\u00E2\u0088\u0082\u00E2\u0088\u0082xk(Ki jkl\u00E2\u0088\u0082uh, j\u00E2\u0088\u0082xl)\u00CF\u0088rd\u00E2\u0084\u00A6 = 0 i = 1 \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 Ndim r = 1 \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 NDOF.Then, using the integration by parts theorem and Equation (3.7), we arrive at the discretized weakform:NDOF\u00E2\u0088\u0091s=1(\u00E2\u0088\u00AB\u00E2\u0084\u00A6\u00E2\u0088\u0082\u00CF\u0088r\u00E2\u0088\u0082xkKi jkl\u00E2\u0088\u0082\u00CF\u0088s\u00E2\u0088\u0082xld\u00E2\u0084\u00A6 +\u00E2\u0088\u00AB\u00E2\u0088\u0082\u00E2\u0084\u00A6\u00CF\u0088rQi j\u00CF\u0088sdS)(Uh,s)j =\u00E2\u0088\u00AB\u00E2\u0088\u0082\u00E2\u0084\u00A6fi\u00CF\u0088rdS, (3.14)for every i \u00E2\u0089\u00A4 Ndim and r \u00E2\u0089\u00A4 NDOF, which is a linear system of equations for finding Uh.The described solution scheme is implemented using the libMesh finite element library [36] in thisthesis. The hierarchical polynomials [22] of order p are chosen as the basis functions, where p is aninput parameter. LibMesh supports quadrature rules of arbitrary orders of accuracy for assembling thediscretized system of equations. Therefore, a quadrature scheme accurate up to order (2p + 1) is used,which ensures exact evaluation of all integrals, so that we do not have to worry about the effect ofquadrature error on the final solution.To illustrate the performance of the implemented method, a model problem is presented, with a simplegeometry resembling a three-dimensional airfoil. This geometry along with its anisotropic mesh areshown in Figure 3.3. The mesh is constituted of 1288 hexahedra. The curved part of the boundary isanalytically parameterized as:222.51.200 00.83XYZFigure 3.3: Geometry of the curved mesh generation test casey(u, v) =[u + (R0 \u00E2\u0088\u0092 bu) sin(v), cu, (R0 \u00E2\u0088\u0092 bu) cos(v) \u00E2\u0088\u0092 H0(1 \u00E2\u0088\u0092 buR0)]T0 < u < 1.25 \u00E2\u0088\u0092 pi6< v n or \u00CE\u00B3m/\u00CE\u00B2 < rtolReturn x(m) = PVm\u00CE\u00B1mfunction Arnoldi(Vm,w)Returns the vector vm+1, which is the unit normal component of w to the space spanned by theprevious vi (i \u00E2\u0089\u00A4 m) vectors.Returns the column vector hm, where (hm)i is the length of the projection of w over the vi vectorif i \u00E2\u0089\u00A4 m + 1, and is zero otherwise.function DenseSolve(\u00CE\u00B2,Hm)Returns \u00CE\u00B1m = argmin\u00CE\u00B1 \u00E2\u0080\u0096b \u00E2\u0088\u0092 APVm\u00CE\u00B1\u00E2\u0080\u0096 and \u00CE\u00B3m = min\u00CE\u00B1 \u00E2\u0080\u0096b \u00E2\u0088\u0092 APVm\u00CE\u00B1\u00E2\u0080\u0096.The objective of this function is achieved by solving the equivalent problemargmin\u00CE\u00B1 \u00E2\u0080\u0096\u00CE\u00B2e1 \u00E2\u0088\u0092 Hm\u00CE\u00B1\u00E2\u0080\u00962 using a rotational change of variables. Solving this problem includesthe solution of a dense m\u00C3\u0097m linear system, and yields both \u00CE\u00B3m and \u00CE\u00B1m. Also, e1 = [1, 0, . . . , 0]T .where D is the control volume diagonal part of A, L is the strict lower part, and U is the strict upperpart. Then, the matrix-by-vector product Pv of the PGS preconditioner is defined as:z(0) = 0z(1) = (D + L)\u00E2\u0088\u00921(v \u00E2\u0088\u0092 A\u00E2\u0088\u0097z(0))...z(n) = (D + L)\u00E2\u0088\u00921(v \u00E2\u0088\u0092 A\u00E2\u0088\u0097z(n\u00E2\u0088\u00921))Pv = z(n),(4.12)where n is the number of PGS iterations. Although PGS is a strong preconditioner for solving linearsystems arising from discretization on isotropicmeshes, it performs poorly when dealingwith anisotropicmeshes [44]. Therefore, PGS is only used for the solution of Poisson\u00E2\u0080\u0099s and Euler equations in this thesis,for which an isotropic mesh correctly captures the solution.32UL DD(a) (b)Figure 4.1: Schematic of Block Jacobi decomposition: (a) The mesh, and the partition belongingto each processor (b) The lower, upper, and block diagonal matrices.4.3.2 Block JacobiThe Block Jacobi (BJ) preconditioner also uses a decomposition of matrix A\u00E2\u0088\u0097, with the difference thatthe diagonal blocks of matrix D must be composed of all the degrees of freedom belonging to the sameprocessor. Figure 4.1a shows a partitioned mesh, and Figure 4.1b shows the corresponding structureand decomposition of the LHS matrix. The matrix-by-vector product for this preconditioner is denotedas:z(0) = 0Dz(1) = v \u00E2\u0088\u0092 A\u00E2\u0088\u0097z(0)...Dz(n) = v \u00E2\u0088\u0092 A\u00E2\u0088\u0097z(n\u00E2\u0088\u00921)Pv = z(n),(4.13)where n is the number of BJ iterations. Note that the block diagonal matrix D cannot be invertedexactly anymore, and each intermediate vector z(k) must be solved for by using an inner iterative solver.Nevertheless, since every block of D only belongs to a single processor, the inner iterative solver canbe serial. Thus, BJ is a means of parallelizing serial preconditioners, and is used for this purpose inANSLib.4.3.3 ILUThe incomplete lower-upper factorization with fill level p (ILUp) produces an approximation of the LUfactorization of A\u00E2\u0088\u0097, such that the approximate matrices have a much smaller nonzero structure than theexact LU factorization. The factored lower and upper triangular matrices L\u00CB\u009C and U\u00CB\u009C are computed byperforming Gaussian elimination on A\u00E2\u0088\u0097, but ignoring certain matrix entries. Then, the preconditionermatrix is set to P = (L\u00CB\u009CU\u00CB\u009C)\u00E2\u0088\u00921. A larger fill level results in L\u00CB\u009C and U\u00CB\u009C matrices with bigger nonzero structure,33but is likely to construct a more effective preconditioner [63]. Nevertheless, the ILUmethod is inherentlya serial algorithm, so it is implemented in conjunction with the BJ method for parallel simulations.Since the ILU preconditioner is based on matrix factoring, its performance is dramatically affected bythe ordering of DOFs in the matrix A\u00E2\u0088\u0097. Quotient minimum degree (QMD) and reverse Cuthill-McKee(RCM) are among the reordering algorithms that are offered in PETSc [10], and are used in ANSLib. Theformer reduces the fill of the factored matrices, while the latter reduces the fill of matrix A\u00E2\u0088\u0097 itself.In the course of development of ANSLib, Nejat and Ollivier-Gooch [49] used an ILU preconditionerfactored from the LHS matrix of the 0-exact scheme to solve the Euler equations. They showedthe good performance of their preconditioning method for 1- and 2-exact finite volume schemes, butconcluded that the 3-exact scheme requires a more effective preconditioner. Later, Michalak andOllivier-Gooch [47] demonstrated that incomplete factorization of the higher-order LHS matrix resultsin faster convergence compared to factoring the 0-exact LHS matrix for inviscid problems. They furtherconjectured that factoring the LHS matrix to the full order using a fill level of 3 has a feasible memoryconsumption, and can be a practical preconditioner for three-dimensional flow problems. Hereinafter,the ILU preconditioner of fill level p factored from the 0-exact and full-order LHS matrices will bedenoted as LO-ILUp and HO-ILUp, respectively.Viscous flow problems are more challenging than their inviscid counterparts. Jalali and Ollivier-Gooch[32] observed that a fill level of three or larger is required forHO-ILUpreconditioning of viscous turbulentflow problems. The HO-ILU3 method is a practical preconditioner for two-dimensional problems, butits memory cost soars drastically in three-dimensions, hindering its implementation for such problems.To mitigate this issue, a new ILU based algorithm will be introduced in the next section that not only hasa less memory consumption than HO-ILU3, but also does not suffer from poor performance of LO-ILUwhen used with the 3-exact reconstruction scheme.4.4 Improved Preconditioning AlgorithmsIn this section, an effective preconditioning algorithm is proposed based on inner GMRES iterations. Fur-thermore, the lines of strong unknown coupling are presented for reordering of the LHS matrix. As willbe shown in Section 4.5, this reordering method can improve the speed of the solver considerably.4.4.1 Inner GMRES IterationsSome researchers have observed that the lower-order LHS matrix can construct a more effective precon-ditioner compared to its higher-order counterpart because the structure of the lower-order LHS matrixonly includes the immediate neighbors of every control volume [41, 56, 78]. Nevertheless, as will beshown in Section 4.5, the LO-ILU method can behave poorly for high-order reconstruction schemesapplied to viscous flow problems. Our conjecture is that using the exact inverse of the lower-order LHS34matrix as the preconditioner: P = (A\u00E2\u0088\u0097)\u00E2\u0088\u00921, rather than the product of the ILU matrices: P = (L\u00CB\u009CU\u00CB\u009C)\u00E2\u0088\u00921, canmitigate the mentioned issue. In this case, the matrix-by-vector product z = Pv can be found by solvingthe system:A\u00E2\u0088\u0097z = v. (4.14)Nevertheless, the linear system of Equation (4.14) can be as large as the original linear system Ax = b,and its solution cannot be carried out exactly. Instead of seeking the exact solution, the proposal is to usefurther inner ILU preconditioned GMRES iterations to find an approximate solution for Equation (4.14).The resulting preconditioner will be referred to as GMRES-LO-ILU in this thesis, and its performancewill be examined in Section 4.5.When an inner GMRES solver is used as the preconditioner, P will not be a linear operator anymore, andits effect changes from iteration to iteration. Therefore, the equation x(m) = PVm\u00CE\u00B1m cannot be used anylonger. The solution to this issue is using a Flexible GMRES (FGMRES) outer solver [62]. The resultingcombination, i.e., a GMRES-LO-ILU preconditioned FGMRES solver, is shown in Algorithm 2. Unlikenormal GMRES, the z vectors are found by applying the inner GMRES solver to the v vectors, and areexplicitly kept track of as the columns of the Z matrix. Moreover, the final approximate solution is writtenas x(m) = Zm\u00CE\u00B1m while the ILU factored matrices are still supplied to the algorithm for preconditioningthe inner GMRES solver. Also, tolerance values and limits for the maximum number of iterations haveto be specified for both the inner and outer solvers which are identified by the subscripts (.)i and (.)o,respectively.Algorithm 2 GMRES-LO-ILU preconditioned FGMRES with no restart and a zero initial guess.function FGMRES(A, A\u00E2\u0088\u0097, L\u00CB\u009C, U\u00CB\u009C, b, no, rtolo, ni, rtoli)\u00CE\u00B2 = \u00E2\u0080\u0096b\u00E2\u0080\u00962v1 = b/\u00CE\u00B2V1 = [v1], H0 = [], Z0 = [],m = 0repeatm = m + 1zm =GMRES(A\u00E2\u0088\u0097, (L\u00CB\u009CU\u00CB\u009C)\u00E2\u0088\u00921, vm, ni, rtoli)w = Azmhm, vm+1 = Arnoldi(Vm,w)Vm+1 = [Vm, vm+1], Hm = [Hm\u00E2\u0088\u00921, hm], Zm = [Zm\u00E2\u0088\u00921, zm]\u00CE\u00B1m, \u00CE\u00B3m =DenseSolve(\u00CE\u00B2,Hm)until m > n or \u00CE\u00B3m/\u00CE\u00B2 < rtolReturn x(m) = Zm\u00CE\u00B1m4.4.2 Lines of Strong Coupling Between UnknownsForming non-overlapping lines that contain control volumes with strongly coupled degrees of freedomis beneficial in constructing effective preconditioners. Mavriplis [43] introduced the concept of lines inanisotropic unstructured meshes to enhance the performance of multigrid solvers via implicit smoothing.35His approximate algorithm for finding such lineswas purely geometric-based, and only considered controlvolume aspect ratios. Thus, his method only captured coupling through diffusion. Okusanya et al. [52]later developed an algorithm that considered both the advection and diffusion phenomena. Fidkowskiet al. [21] proved that such an algorithm results in a set of unique lines, and used the lines for nonlinear p-multigrid solution of Euler equations. Diosady and Darmofal [20] showed that reordering the unknowns,such that members of a line have consecutive numbers, is an effective reordering strategy for the ILUpreconditioner. In this thesis, the line reordering algorithm of [20] is implemented to improve the speedof the linear solver.In the first step of the line creation algorithm, a weight is assigned to every face inside the mesh.This weight is derived from the Jacobian of the 0-exact discretization of the advection-diffusion equa-tion,\u00E2\u0088\u0087 \u00C2\u00B7 (vu \u00E2\u0088\u0092 \u00C2\u00B5L\u00E2\u0088\u0087u) = 0, (4.15)where u is the unknown variable, v is the velocity taken from the initial conditions, and \u00C2\u00B5L is an inputparameter that controls the sensitivity of the lines to mesh anisotropy. The weight of a face f is thendefined as:W( f ) = max(\u00E2\u0088\u0082R\u00CF\u0083\u00E2\u0088\u0082u\u00CF\u0084,\u00E2\u0088\u0082R\u00CF\u0084\u00E2\u0088\u0082u\u00CF\u0083), (4.16)where W is the weight of the face, \u00CF\u0083 and \u00CF\u0084 are the adjacent control volumes of the face, and R theresidual vector. If a face is located on the boundary, Equation (4.16) will not directly be applicable to it.Thus, a mirrored ghost control volume is placed adjacent to every boundary face. After finding the faceweights, calling Algorithm 3 will construct the lines, where F(\u00CF\u0084) is the set of faces of control volume\u00CF\u0084, and \u00CF\u0083 = C(\u00CF\u0084, f ) is the control volume which has the face f in common with control volume \u00CF\u0084. Thisalgorithm ensures that two adjacent control volumes are connected to each other if and only if they havetheir first or second highly weighted face lying between them [21].To illustrate the performance of the algorithm, an anisotropic mesh for the geometry of a NACA 0012airfoil is considered, as shown in Figure 4.2. The velocity has been taken from the 2-exact solution ofthe flow resulting from a Mach number of 0.15, a Reynolds number of 6 \u00C3\u0097 106, and an angle of attackof 10\u00E2\u0097\u00A6. The parameter \u00C2\u00B5L is set to the heuristic value of 10\u00E2\u0088\u00925. As desired, the lines follow the directionof mesh anisotropy near wall boundaries, while their pattern changes, and follows the flow direction inother parts of the geometry.4.5 Numerical ComparisonsIn this section, the performance of the proposed preconditioning algorithm is studied by considering theturbulent flow around a NACA 0012 airfoil, with \u00CE\u00B1 = 10\u00E2\u0097\u00A6, Ma = 0.15, and Re = 6 \u00C3\u0097 106. Three nestedmeshes with approximately 25K, 100K, and 400K control volumes (NCV = 25K,100K, and 400K) areconsidered, where the coarsest mesh and its corresponding lines of strong unknown coupling are shownin Figure 4.2. The farfield boundaries are located almost 500 chords away from the airfoil, and the chord36Algorithm 3 Line creationfunction CreateLines . Calling this function creates all lines.Unmark all control volumes.while there a exists an unmarked control volume \u00CF\u0084 doMakePath(\u00CF\u0084) . Forward path creationMakePath(\u00CF\u0084) . Backward path creationfunctionMakePath(\u00CF\u0084) . Helper function.repeatFor control volume \u00CF\u0084, pick the face f \u00E2\u0088\u0088 F(\u00CF\u0084)with the highest weight, such that control volume\u00CF\u0083 = C(\u00CF\u0084, f ) is not part of the current line.if f is a boundary face thenTerminateelse if \u00CF\u0083 is marked thenTerminateelse if f is not the first or second ranked face of \u00CF\u0083 in weight thenTerminateMark \u00CF\u0083.Add \u00CF\u0083 to the current line.\u00CF\u0084 \u00E2\u0086\u0090 \u00CF\u0083until Terminationlength has a nondimensional value of one. This problem is taken from the NASA Turbulence ModelingResource (TMR) website [60], and has been previously studied by Jalali and Ollivier-Gooch [32], wherethey used a HO-ILU3 preconditioner with QMD reordering. Only the highest-order discretizationscheme of ANSLib (3-exact) is considered here, since it results in the stiffest linear systems that are mostdifficult to solve. Also, the initial conditions are taken from the steady-state solution of the lower-order2-exact scheme.The ILU based preconditioning methods tested are shown in Table 4.1. In all cases, the outer GMRESsolver stops when the linear residual of Equation (4.3) is reduced by a factor of 10\u00E2\u0088\u00923, or a maximumnumber of 500 iterations is performed. The maximum Krylov subspace size is 100 for the outer GMRESsolver, and the initial CFL number is set to 10\u00E2\u0088\u00922. Furthermore, the simulation ends when the norm ofthe residual vector \u00E2\u0080\u0096R(Uh)\u00E2\u0080\u00962 is reduced by a factor of 10\u00E2\u0088\u00928. The command-line options supplied to theANSLib executable for each case are listed in Appendix A. Cases A-E are compared to each other onthe 25K and 100Kmeshes, while case E\u00E2\u0088\u0097 is used to solve the problem on the finest mesh. Note that othercombinations of preconditioner and reordering schemes that are not considered in Table 4.1, did notresult in convergence even for the coarse mesh. Moreover, an extensive search for finding the optimumnumber of inner iterations ni has not been carried out, and the large number of inner iterations for caseE\u00E2\u0088\u0097 is chosen as a safe measure to ensure convergence on the finest mesh.The residual history of the solver for each preconditioning algorithm is shown in Figure 4.3, where thewalltime is obtained from a single core of an Intel i7-4790 (3.60 GHz) CPU. Our proposed preconditioning37Block ColorFigure 4.2: Mesh and lines of strong unknown coupling for the geometry of a NACA 0012 airfoilTable 4.1: Preconditioning methods consideredCase Preconditioning Reordering rtoli, niname method algorithmA HO-ILU3 QMD \u00E2\u0080\u0094B LO-ILU0 RCM \u00E2\u0080\u0094C LO-ILU0 lines \u00E2\u0080\u0094D GMRES-LO-ILU0 RCM 10\u00E2\u0088\u00923, 10E GMRES-LO-ILU0 lines 10\u00E2\u0088\u00923, 10E\u00E2\u0088\u0097 GMRES-LO-ILU0 lines 10\u00E2\u0088\u00923, 60380.0 0.2 0.4 0.6 0.8 1.0Wall time (s)0.00.20.40.60.81.0|R(Uh)|20 200 400 600 80010\u00E2\u0088\u0092610\u00E2\u0088\u0092410\u00E2\u0088\u00922100102104NCV=25K0 1000 2000 3000 4000 5000NCV=100KCase ACase BCase CCase DCase EFigure 4.3: Comparison of residual histories for the NACA 0012 problem obtained using differentpreconditioning algorithms05\u00C3\u00971041\u00C3\u0097105Wall time (s)10-610-410-2100102\u00E2\u0080\u0096R(Uh)\u00E2\u0080\u00962case E\u00E2\u0088\u0097Figure 4.4: Residual history for the NACA 0012 problem on the finest mesh of 400K controlvolumesalgorithm (case E) outperforms all the other methods on both mesh sizes, and takes only half the timeof the HO-ILU3 algorithm (case A) to find the steady-state solution. Furthermore, the line reorderingalgorithm shows its prominence on the 100K control volumemesh, where the RCM reordering algorithmfails. The effectiveness of the new preconditioning algorithm is further demonstrated by solving theproblem on a 400K control volume mesh. The residual history for this case is shown in Figure 4.4.Table 4.2 shows the number of PTC iterations (N-PTC), the number of outer GMRES iterations (N-GMRES), totalmemory consumption, CPU time spent on the linear solver (LST), and total computationaltime (TST). The proposed preconditioning scheme (case E) speeds up the linear solve process by a factorof three, and results in a twofold reduction in memory consumption, compared to the HO-ILU3 method.The shortcomings of the LO-ILU methods (cases B and C) are evident because of their large number ofPTC iterations and failure in convergence, on the coarse and medium meshes, respectively. Also, the39Table 4.2: Performance comparison for different preconditioning schemes. Cases with entriesmarked by \u00E2\u0080\u009C\u00E2\u0088\u0092\u00E2\u0080\u009D did not converge to the steady state solution.Preconditioner N-PTC N-GMRES Memory(GB) LST(s) TST(s)NCV = 25KA 37 956 1.4 416 635B 44 10, 893 0.7 264 572C 51 13, 692 0.7 328 705D 34 1, 856 0.7 134 379E 34 1, 787 0.7 118 361NCV = 100KA 39 4, 640 5.9 3, 122 4, 068B \u00E2\u0088\u0092 \u00E2\u0088\u0092 2.8 \u00E2\u0088\u0092 \u00E2\u0088\u0092C \u00E2\u0088\u0092 \u00E2\u0088\u0092 2.8 \u00E2\u0088\u0092 \u00E2\u0088\u0092D \u00E2\u0088\u0092 \u00E2\u0088\u0092 3.2 \u00E2\u0088\u0092 \u00E2\u0088\u0092E 36 4, 396 3.2 1, 308 2, 348NCV = 400KE\u00E2\u0088\u0097 38 6, 869 13.1 87, 312 91, 502Table 4.3: Computed drag and lift coefficients for the NACA 0012 airfoilNCV CDp CDv CL25K 0.00545 0.00583 1.0951100K 0.00615 0.00623 1.0905400K 0.00609 0.00619 1.0922NASA TMR 0.00607 0.00621 1.0910memory consumption scales linearly with respect to the problem size for the proposed preconditioningscheme, whereas the linear solve time seems to be increasing at a much more rapid rate. Although theproposed preconditioning scheme has a huge linear solve time for NCV = 400K, it still outperforms theHO-ILU3 method, which cannot handle this mesh size at all.The viscous drag (CDv), the pressure drag (CDp), and the lift (CL) coefficients computed from thenumerical solution along with their reference values are listed in Table 4.3. The reference values aretaken from the NASA TMR website, and are obtained by the FUN3D [2] solver running on a veryfine mesh (approximately 2M quadrilaterals). As expected, the change in the coefficients between thefine and medium meshes is much smaller than that of the medium and coarse meshes. Furthermore,the gap between the coefficients and their reference values gets smaller with mesh refinement, with theexception ofCL . Solution singularities, different boundary condition implementations between ANSLiband FUN3D, or inexact evaluation of the distance from wall function may partly be accountable for thenon-ideal behavior of the lift coefficient. Nevertheless, the purpose of this section was the solution ofthe discretized system of equations, which according to Figure 4.4 has been performed correctly.40Chapter 5Three-Dimensional ResultsThis chapter presents the k-exact finite volume solution of four three-dimensional problems. Theproblems are studied in order of difficulty, and the solutions are verified to assess the accuracy andperformance of the solver. First, Poisson\u00E2\u0080\u0099s equation is solved in a cubic domain, where the source termis obtained from a manufactured solution [64]. Then, the Euler equations are solved to simulate thesubsonic inviscid flow around a unit sphere. Subsequently, the solution of Navier-Stokes + negative S-Aequations are presented to simulate the subsonic viscous turbulent flow over a flat plate and the extrudedNACA 0012 airfoil. All the command-line options provided to the solver for running each case areprovided in Appendix A.All the flow simulations were performed using the Grex cluster from the WestGrid computing facilities[3]. Each node of the cluster consists of two 6-core Intel Xeon X5650 2.66GHz processors, and has aminimum of 48GB memory. All compute nodes are connected by a non-blocking Infiniband 4X QDRnetwork. The wall time and memory cost of all the flow problems are reported using 8 processors. Inaddition, a strong parallel scaling test is conducted for the flat plate and the sphere test cases, i.e., thesame problem is solved by using different numbers of processors ranging from 1 to 10 while the resultingspeedup is recorded.5.1 Poisson\u00E2\u0080\u0099s Equation in a Cubic DomainFor this test case, Poisson\u00E2\u0080\u0099s equation is considered inside the domain \u00E2\u0084\u00A6 = [0, 1]3. The manufacturedexact solutionu = sinh(sin(x1)) sinh(sin(x2)) sinh(sin(x3)) (5.1)is used to find the source term, while homogeneous Dirichlet conditions are imposed over all theboundaries. Slice plots of the exact solution are shown in Figure 5.1. The solution has a maximum valueof (sinh(1))3 at the point (12, 12, 12 ), and decays to zero on the boundaries.41Figure 5.1: Exact solution of Poisson\u00E2\u0080\u0099s problemA study of the solution accuracy is performed on four families of unstructured grids. The first gridfamily is only composed of hexahedra while the other families are composed of only prisms, mixtureof pyramids and tetrahedra, and only tetrahedra, respectively. The hexahedral grids are constructed byperturbing the vertices of a uniform N \u00C3\u0097 N \u00C3\u0097 N structured mesh. Other grid families are subsequentlycreated by decomposing each hexahedron into the desired element types. Figure 5.2 shows all the gridsfor N = 5.The problem is solved on grid sizes: N = 10, 20, 40, and reconstruction orders: k = 1, 2, 3. Figure 5.3shows the L2 norm of the discretization error \u00E2\u0080\u0096eh \u00E2\u0080\u00962 as a function of the mesh length scale h = 1N andthe number of degrees of freedom. Optimal accuracy order of k + 1 is attained for the k = 1 and 3reconstruction schemes as expected, whereas the k = 2 scheme only achieves a second-order accuratesolution. Nevertheless, the non-optimal behavior of k = 2 is expected for Poisson\u00E2\u0080\u0099s problem, and waspreviously addressed and theoretically explained by Ollivier-Gooch and Van Altena [54].5.2 Inviscid Flow Around a SphereThe second test problem is the inviscid flow around a sphere with unit radius, and centered at the point(0, 0, 0). The Mach number is equal to Ma = 0.38, and the free-stream flow is in the x1 direction. Thisproblem is chosen as a three-dimensional extension of the inviscid flow around a circle, the higher-ordersolution of which was studied by Bassi and Rebay [14]. Three prismatic grids (64K, 322K, and 1Mcontrol volumes) are used, where the farfield is located at a distance of 100 dimensionless units from thesphere surface. A symmetric cut of the coarsest mesh near the sphere surface is shown in Figure 5.4.The faces and control volumes adjacent to the wall boundaries are curved using only second-orderLagrange polynomials, as quadratic functions are sufficient to represent a spherical surface. To ensurethat numerical quadrature does not introduce additional error in the solution, a quadrature scheme42(a) (b)(c) (d)Figure 5.2: Different meshes for the solution of Poisson\u00E2\u0080\u0099s problem: (a) hexahedra, (b) prisms, (c)mixture of pyramids and tetrahedra, (d) tetrahedra.accurate up to order k + 2 is employed.In all cases, the initial solution is the uniform free-stream condition everywhere, and the initial CFLnumber is set to 10. The GMRES linear solver stops when the linear residual of Equation (4.3) is 10\u00E2\u0088\u00923times its initial value, or a maximum of 500 iterations is performed. Furthermore, the GMRES solverhas a maximum Krylov subspace size of 100, and is preconditioned using the BJ preconditioner withno = 4 iterations. The inner solver for the BJ preconditioner is PGS with ni = 4 iterations. Althoughthe LHS matrix of the linear system is calculated to full-order, its preconditioner is constructed from theLHS matrix of the 0-exact reconstruction scheme.As the flow does not contain any discontinuity, the entropy must be constant throughout the domain.Therefore, the L2 norm of the entropy relative to the free-stream conditions, \u00E2\u0080\u0096S \u00E2\u0088\u0092 S\u00E2\u0088\u009E\u00E2\u0080\u00962, has an exactvalue of zero, and is used to study the solution accuracy. The convergence of \u00E2\u0080\u0096S \u00E2\u0088\u0092 S\u00E2\u0088\u009E\u00E2\u0080\u00962 with meshrefinement is shown in Figure 5.5, where the mesh length scale is defined as h = (NCV)\u00E2\u0088\u00921/3. The k = 1and 3 reconstruction schemes converge even faster than their expected ratios of 2 and 4, respectively.4312 3 4h0/h10\u00E2\u0088\u0092610\u00E2\u0088\u0092510\u00E2\u0088\u0092410\u00E2\u0088\u0092310\u00E2\u0088\u0092210\u00E2\u0088\u00921|eh|2k=112 3 4h0/hk=212 3 4h0/hk=3Hex Prism Pyramid+Tet Tet O(hk+1)(a)103104105106NDOF10\u00E2\u0088\u0092610\u00E2\u0088\u0092510\u00E2\u0088\u0092410\u00E2\u0088\u0092310\u00E2\u0088\u0092210\u00E2\u0088\u00921|eh|2k=1103104105106NDOFk=2103104105106NDOFk=3(b)Figure 5.3: Discretization error versus mesh size for Poisson\u00E2\u0080\u0099s problem: (a) Error versus meshlength scale, (b) Error versus number of degrees of freedom.Figure 5.4: Symmetric cut of the coarsest mesh near the sphere surface441 2 3 4 5 6h0/h10-710-610-510-410-3\u00E2\u0080\u0096Sh\u00E2\u0088\u0092S\u00E2\u0088\u009E\u00E2\u0080\u00962ave. slope=2.2ave. slope=2.4ave. slope=4.8k=1k=2k=3Figure 5.5: Relative entropy norm versus mesh size for the sphere problemThe convergence ratio for the k = 2 reconstruction scheme, however, is half an order smaller than itsnominal value. The slight deviation of the convergence orders from their theoretical values can partlybe attributed to the fact that the meshes employed are not nested.The effect of the reconstruction order k on solution accuracy is also evident from the qualitativeMach contours on the x3 = 0 symmetry plane, shown in Figure 5.6. Not surprisingly, only the k = 3reconstruction scheme has been able to capture the symmetry of the flow on the finestmesh. Furthermore,the artificialwake behind the sphere seems smaller for the k = 2 solution compared to that of k = 1.The norm of the residual vector per PTC iteration is shown in Figure 5.7. Convergence is independent ofthe reconstruction order k, but the number of iterations increases as the mesh is refined. It is noteworthyto mention that the increase in the initial norm corresponding to mesh refinement is due to the definitionof the residual vector R(Uh) in Equations (2.3) and (2.4). Each entry of the residual vector is divided bythe size of its corresponding control volume. Therefore, when the mesh is refined, the control volumesizes become smaller, and the residual vector entries will increase.The number of iterations and the resource consumption of the flow solver are listed in Table 5.1. Thetime spent on linear solves makes up less than twenty percent of the computational time, suggesting thatthe Jacobian evaluation is the major computational bottleneck for this problem. Either increasing thereconstruction order, or the number of control volumes results in stiffer linear systems, which requiremore GMRES iterations to be solved. Nevertheless, both the solution time and memory consumptionincrease linearly with the number of control volumes for all k. Also, the advantage of using a high-orderreconstruction scheme can be clearly seen: According to Figure 5.5, the k = 1 scheme needs at leastone level of uniform mesh refinement on the finest grid to achieve an entropy norm of 10\u00E2\u0088\u00927. While sucha mesh refinement would increase the resource consumption by a minimum factor of eight, the k = 3scheme achieves the same level of accuracy by consuming only two times the same resources.To test the parallel scalability of the flow solver, the sphere problem was solved on different numbers of45NCV = 64K NCV = 1Mk=1k=2k=3Figure 5.6: Computed Mach contours on the x3 = 0 symmetry plane for the sphere problem460.0 0.2 0.4 0.6 0.8 1.0PTC iteration0.00.20.40.60.81.00 10 2010-1210-1010-81-610-410-2100102\u00E2\u0080\u0096R(Uh)\u00E2\u0080\u00962k=10 10 20k=20 10 2k=3NCV = 64KNCV = 322KNCV = 1MFigure 5.7: Norm of the residual vector per PTC iteration for the sphere problemTable 5.1: Number of iterations and the resource consumption of the flow solver for the sphereproblemk N-PTC N-GMRES Memory(GB) LST(s) TST(s)NCV = 64K1 15 182 2.94 24 1072 15 181 3.91 20 1123 15 255 6.00 31 320NCV = 322K1 16 207 13.24 134 5792 16 212 14.05 132 6203 16 300 28.39 271 1, 879NCV = 1M1 17 275 39.48 486 1, 9692 17 277 45.52 536 2, 1503 17 385 85.78 793 5, 904processors ranging from 1 to 10 with NCV = 322K and k = 3. The parallel speedup of the total solutionprocess, the Jacobian evaluation algorithm, and the linear solver scheme were computed separately, andare shown in Figure 5.8. Not surprisingly, the linear solver algorithm is less scalable than that of theJacobian integrator mainly because the performance of the BJ preconditioner decreases with an increasein the number of processors. The sudden dip of the parallel speedup for 9 processors may be due toa less favorable mesh partitioning, or communication problems on the cluster. Regardless, the overallsolution algorithm scales very well in parallel, as the linear solver only constitutes a small part of theoverall process.471 2 3 4 5 6 7 8 9 10# of processors12345678910Speed-upIdealTotalJac. int.Lin. solveFigure 5.8: The parallel speedup of the solver for the sphere problem: NCV = 322K and k = 3.5.3 Turbulent Flow Over a Flat PlateFor this test case, we consider a three-dimensional extension of the flat plate verification case fromthe NASA TMR website. The computational domain is \u00E2\u0084\u00A6 = [\u00E2\u0088\u00920.33, 2] \u00C3\u0097 [0, 1] \u00C3\u0097 [0, 1], and thenondimensional parameters have the values: Re = 5 \u00C3\u0097 106, Ma = 0.2, \u00CE\u00B1 = 0, and \u00CF\u0088 = 0. The flat plateis located on the (0 \u00E2\u0089\u00A4 x1 \u00E2\u0089\u00A4 2) \u00E2\u0088\u00A7 (x2 = 0) boundary, where adiabatic solid wall conditions are imposed.Symmetry boundary conditions are imposed on the (x3 = 0), (x3 = 1), and (\u00E2\u0088\u00920.33 \u00E2\u0089\u00A4 x1 \u00E2\u0089\u00A4 0) \u00E2\u0088\u00A7 (x2 = 0)boundaries, while other boundaries are considered as farfield. As symmetry boundary conditions areimposed on all the boundaries normal to the x3-axis, the exact solution of this problem must be constantin the x3 direction. Therefore, the solutions obtained in this section can be compared to the two-dimensional results from the NASA TMR website, which are obtained by the CFL3D solver [1] on a544 \u00C3\u0097 384 grid.The two-dimensional meshes provided by the NASA TMR website are extruded in the x3 direction toconstruct a series of nested three-dimensional grids with dimensions: 60 \u00C3\u0097 34 \u00C3\u0097 7, 120 \u00C3\u0097 68 \u00C3\u0097 14, and240 \u00C3\u0097 136 \u00C3\u0097 28. The coarsest mesh is depicted in Figure 5.9. An anisotropic mesh is necessary tocorrectly capture the flow pattern in the boundary layer and at the plate leading edge. On each mesh, theproblem is solved for k = 1, then for k = 2, and finally for k = 3. The free-stream conditions are usedas the initial conditions for k = 1, while the initial solution for each k \u00E2\u0089\u00A5 2 is taken from the convergedsolution of the (k \u00E2\u0088\u0092 1)-exact scheme. The GMRES-LO-ILU2 method is used as the preconditioner,with 10 inner GMRES iterations, and the RCM reordering algorithm. Moreover, the ILU factorizationis performed only on the diagonal block of each processor. Other solver parameters are the same asSection 5.2.To schematically demonstrate the correctness of the solution, the distribution of the turbulence workingvariable obtained from k = 3 on the finest mesh is plotted on the x3 = 0.5 plane, and shown inFigure 5.10a. The results for a 544\u00C3\u0097384 grid provided by theNASATMRare also shown in Figure 5.10b.48XYZFigure 5.9: Coarsest mesh for the flat plate problem.(a) 3-exact: min : 7 \u00C3\u0097 10\u00E2\u0088\u00923, max : 379 (b) NASA TMR: min : 0, max : 378Figure 5.10: Distribution of the turbulence working variable on the plane x3 = 0.5 for the flat plateproblemThere is a close agreement between the reference solution of NASA TMR and the solution of this thesis,although a few undershoots of small magnitude are present in the latter. The next quantity of interestis the distribution of the nondimensional eddy viscosity \u00C2\u00B5T\u00C2\u00B5 on the line (x1 = 0.97) \u00E2\u0088\u00A7 (x3 = 0.5),which is depicted in Figure 5.11. Note how the solution discontinuity at control volume boundaries isquite evident for k = 1, while higher-order reconstruction schemes result in much smoother solutions.The reference solution from the NASA TMR website is also shown in Figure 5.11. An acceptableagreement is observable between our computed results and the reference values. Also, it is evident thata higher-order reconstruction scheme leads to a more accurate estimate of the eddy viscosity.To further verify the numerical solution, the convergence of the drag coefficient CD and the skinfriction coefficient Cf at the point x = (0.97, 0, 0.5) are studied with mesh refinement. Table 5.2 showsthe computed values on different meshes. The computed values converge faster for a higher-order490.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0x20.000.010.020.0360\u00C3\u009734\u00C3\u00977NASA TMRk=1k=2k=30.000.010.020.03120\u00C3\u009768\u00C3\u0097140 50 100 150 200 250\u00C2\u00B5T/\u00C2\u00B50.000.010 020.03240\u00C3\u0097136\u00C3\u009728Figure 5.11: Distribution of the nondimensional eddy viscosity on the line (x1 = 0.97)\u00E2\u0088\u00A7(x3 = 0.5)for the flat plate problemreconstruction scheme, such that only k = 3 offers an accurate solution on the medium mesh. Thereference values from the NASA TMR library are also shown in the same table. As desired, there is agood agreement between the computed values using the method of this thesis and the reference valuesof the NASA TMR website. The convergence orders of CD and Cf are computed using the procedure ofCelik et al. [16], and are also shown in Table 5.2. All the reconstruction schemes have attained optimalconvergence rates. Nevertheless, care must be taken when interpreting the estimated convergence ratesbecause the procedure of Celik et al. is most reliable when the error has already achieved an asymptoticbehavior on the coarsest mesh, which might not necessarily be the case here.The norm of the residual vector per PTC iteration is shown in Figure 5.12. Similar to Section 5.2,convergence is not affected by the reconstruction order k, but is slightly degraded as the number ofdegrees of freedom increases. The number of iterations and the resource consumption of the flow solverare listed in Table 5.3. The majority of the total time is spent on linear solves, showing that the linearsystems arising from viscous turbulent flows are much stiffer compared to their inviscid counterparts.50Table 5.2: Computed value and convergence order of the drag coefficient and the skin frictioncoefficient at the point x = (0.97, 0, 0.5) for the flat plate problemCD CfNASA TMR 0.00286 0.00271Meshk 1 2 3 1 2 360 \u00C3\u0097 34 \u00C3\u0097 7 0.00396 0.00233 0.00233 0.00350 0.00228 0.00222120 \u00C3\u0097 68 \u00C3\u0097 14 0.00301 0.00281 0.00285 0.00283 0.00268 0.00271240 \u00C3\u0097 136 \u00C3\u0097 28 0.00287 0.00286 0.00286 0.00274 0.00273 0.00273Convergence order 2.8 3.3 5.4 3 3 5.10 10 20 30 40 50 60 70 80 90PTC iteration10-910-710-510-310-1101103105\u00E2\u0080\u0096R(Uh)\u00E2\u0080\u00962k=1 k=2 k=360\u00C3\u009734\u00C3\u00977120\u00C3\u009768\u00C3\u009714240\u00C3\u0097136\u00C3\u009728Figure 5.12: Norm of the residual vector per PTC iteration for the flat plate problemAs desired, memory consumption increases linearly with mesh refinement, whereas the linear solve timedoes not scale as nicely as for the inviscid sphere problem. Most importantly, the benefit of higher-ordermethods can once again be observed: While the k = 3 solution on the medium mesh offers the samelevel of accuracy as the k = 1 solution on the fine mesh, the computational time of the former is smallerthan that of the latter by a factor of four.A strong parallel scaling test is conducted for the k = 3 scheme on the 120\u00C3\u009768\u00C3\u009714 mesh. The resultingspeedup for different parts of the solver is shown in Figure 5.13, which is quite similar to the resultsof the sphere problem, and suggests that the preconditioning algorithm is performing well for viscousflow problems. Nevertheless, the net scaling of the solver degrades marginally compared to the sphereproblem, as linear solves take up a bigger portion of the total solution time.5.4 Turbulent Flow Over an Extruded AirfoilFor the last test case, we consider a three-dimensional extension of the flow around the NACA 0012airfoil. The computational domain of Section 4.5 is extruded in the x3 direction with an extrusion length51Table 5.3: Number of iterations and the resource consumption of the flow solver for the flat plateproblemk N-PTC N-GMRES Memory(GB) LST(s) TST(s)60 \u00C3\u0097 34 \u00C3\u0097 7 mesh1 26 844 0.42 31 552 26 1, 009 1.35 42 1243 26 1, 071 2.10 57 202120 \u00C3\u0097 68 \u00C3\u0097 14 mesh1 28 1, 436 5.24 510 7132 29 1, 864 8.30 742 1, 4893 29 2, 041 14.63 825 2, 124240 \u00C3\u0097 136 \u00C3\u0097 28 mesh1 29 2, 492 38.77 6, 222 7, 8842 27 3, 305 60.70 12, 353 18, 1373 27 2, 906 121.81 23, 141 35, 4121 2 3 4 5 6 7 8 9 10# of processors12345678910Speed-upIdealTotalJac. int.Lin. solveFigure 5.13: The parallel speedup of the solver for the flat plate problem: 120 \u00C3\u0097 68 \u00C3\u0097 14 mesh andk = 3.of one nondimensional unit. Symmetry boundary conditions are imposed on the two boundaries normalto the x3-axis, while other boundaries retain their conditions from Section 4.5. The nondimensionalparameters are: Re = 6 \u00C3\u0097 106, Ma = 0.15, \u00CE\u00B1 = 10, and \u00CF\u0088 = 0. As the exact solution must be constant inthe x3 direction, the solutions are compared to the reference values provided by NASA TMR, which areobtained by the FUN3D [2] solver on a 7169 \u00C3\u0097 2049 grid.Two meshes are employed: a hexahedral mesh with NCV = 100K, and a mixed prismatic-hexahedralmesh with NCV = 176K. In both cases, quadrature points are obtained from the tensor product of one-dimensional quadrature formulas and the quadrature points of the curved two-dimensional meshes. Aportion of the meshes is depicted in Figure 5.14. In both cases, there are 7 layers of extruded controlvolumes in the x3 direction. A value of MinNeigh(3) = 30 performed acceptably for the hexahedralmesh, whereas the same value resulted in an ill-conditioned version of Equation (2.7) for the mixed52XYZ(a) Hexahedral meshXYZ(b) Mixed prismatic-hexahedral meshFigure 5.14: Meshes for the extruded NACA 0012 problemmesh. Therefore, MinNeigh(3) = 60 was used in the latter case as a safe measure to ensure a well-posed minimization problem for Equation (2.7). Other solver parameters, including the order rampingprocedure for solution initialization, are the same as Section 5.3.The computed contours of turbulence working variable on the x3 = 0 plane are shown in Figure 5.15for the k = 3 solutions. The mixed-mesh resolves the wake region more accurately, while providing asmoother distribution of the turbulence working variable. The superior performance of the mixed meshcan be explained by its bigger number of degrees of freedom, and the presence of extra flow alignedfaces that are missing in the hexahedral mesh.The distribution of the surface pressure coefficient on the intersection of the airfoil and the x3 = 0.5plane is computed using the k = 1 and 3 solutions, and depicted in Figure 5.16. Generally, the computeddistributions are consistent with the reference values obtained by FUN3D. Closer views of the plots areshown for four distinct points to better compare the results. As expected, the k = 3 scheme predicts themost accurate values, particularly on the mixed mesh. For example, the FUN3D results show that thepressure coefficient on the lower surface becomes bigger than that of the upper surface, near the trailingedge of the airfoil. This phenomenon is captured by the k = 3 solutions, but not observed by those ofk = 1.53(a) Hexahedral mesh(b) Mixed prismatic-hexahedral meshFigure 5.15: Distribution of the turbulence working variable for the extruded NACA 0012 problemon the x3 = 0 planeTable 5.4: Computed drag and lift coefficients for the extruded NACA 0012 airfoil problemk CDp CDv CLNASA TMR\u00E2\u0088\u0092 0.00607 0.00621 1.0910Hex mesh1 0.01703 0.00582 1.06192 0.01702 0.00497 1.05073 0.00301 0.00472 1.0417Mixed mesh1 0.01129 0.00574 1.07352 0.00365 0.00565 1.07763 0.00550 0.00536 1.0869The computed and the reference lift and drag coefficients are shown in Table 5.4. It seems that thehexahedral mesh is too coarse, as none of the schemes can obtain accurate enough solutions. For themixed mesh, the k = 3 scheme gives satisfactory results with less than 10% error for all the coefficients.The solution of the other schemes, however, is quite off, particularly for the pressure drag coefficientCDp. Note that the reference results are obtained using a 7169 \u00C3\u0097 2049 mesh, which is more than athousand times finer than the meshes employed in this thesis.The norm of the residual vector per PTC iteration is shown in Figure 5.17. Convergence is marginallyaffected by either the mesh type or the reconstruction order k. The latter is a desirable result of usingthe steady-state solution of a k-exact scheme as the initial guess for the (k + 1)-exact solution.540.0 0.2 0.4 0.6 0.8 1.0x\u00E2\u0088\u00926\u00E2\u0088\u00925\u00E2\u0088\u00924\u00E2\u0088\u00923\u00E2\u0088\u00922\u00E2\u0088\u00921012Cp(a)(b)(c)(d)FUN3Dk=3, hexk=3, mixedk=1, hexk=1, mixed(a) (b)(c) (d)Figure 5.16: Distribution of surface pressure coefficient on the intersection of the extruded NACA0012 airfoil and the x3 = 0.5 plane.0 20 40 60 80 100PTC iteration10-710-510-310-1101103105\u00E2\u0080\u0096R(Uh)\u00E2\u0080\u00962k=1 k=2 k=3HexMixedFigure 5.17: Norm of the residual vector per PTC iteration for the extruded airfoil problem55Table 5.5: Resource consumption of the flow solver for the extruded NACA 0012 problemk N-PTC N-GMRES Memory(GB) LST(s) TST(s)Hex mesh, NCV = 100K1 33 1, 154 4.77 317 7442 31 1, 788 6.82 730 2, 0973 31 2, 415 12.23 1, 057 3, 215Mixed mesh, NCV = 176K1 34 1, 132 8.70 458 1, 1642 32 1, 769 10.87 800 2, 4273 31 2, 185 26.47 1, 311 4, 666Finally, parameters related to the iterative convergence of the solver are listed in Table 5.5. The linearsolve time makes up around 30% of the total computation time. The total computation time is stronglydependent on k, but does not differ considerably between the two meshes for the same reconstructionorder. Moreover, the memory consumption seems to be linearly increasing with both k and NCV.56Chapter 6ConclusionsA higher-order solution framework was presented for steady-state three-dimensional compressible flows.The highlights of this solution framework include: a k-exact finite volume discretization, a PTC solutionalgorithm, a memory lean preconditioner based on inner GMRES iterations, and an ILU reorderingalgorithm based on lines of strong coupling between them.Higher-order accuracy was achieved by constructing a piecewise continuous representation of the controlvolume average values. This piecewise continuous representation, also referred to as the discrete solution,was defined as the superposition of a set of basis functions. Monomials of Cartesian, principal, orcurvilinear coordinates were used as the basis functions while the coefficient of each monomial in thesuperposition was found by solving a constrained minimization problem. In the finite volume method,the convective fluxes were evaluated by the Roe flux function. The viscous fluxes were obtained byadding a stabilizing damping term to the averaged gradients of the two adjacent states. Specifically, theflux functions were presented for the RANS equations fully coupled with the negative S-A turbulencemodel.The curvature of the domain boundaries must be correctly accounted for in high-order numericaldiscretization schemes, yet simply replacing the boundary faces of highly anisotropic meshes withhigher-order representations can result in invalid intersecting elements. Thus, a three-dimensional finiteelement elasticity solver was developed to propagate the curvature of the boundary throughout thedomain, and to prevent mesh tangling. The developed solver was tested by solving a model problem thatimitated the geometry of a three-dimensional airfoil. Also, the principal coordinates and the curvilinearcoordinates basis functions were presented to prevent the k-exact reconstruction scheme from behavingpoorly on curved anisotropic meshes.The PTCmethodwas revisited for the solution of the discretized nonlinear equations, where each iterationrequired the solution of a linear system. The challenges involved in the solution of these linear systemswere addressed for three-dimensional problems, and a novel method was introduced to mitigate theseissues. In this method, the FGMRES linear solver was preconditioned using inner GMRES iterations,57which were subsequently preconditioned by the ILU method. The L\u00CB\u009C and U\u00CB\u009C matrices were factoredfrom the LHS matrix of the 0-exact discretization scheme, and the lines of strong unknown couplingwere employed to reorder the unknowns. By considering the 2-D fully turbulent flow around the NACA0012 airfoil, it was shown that the proposed preconditioning algorithm is more efficient both in termsof solution time and memory consumption, compared to the other ILU based preconditioning methodsconsidered.The correct implementation of element quadrature and connectivity information in the solver was verifiedby solving Poisson\u00E2\u0080\u0099s equation inside a cubic domain. Subsequently, inviscid and viscous turbulent testproblems were considered, where the solution was verified against reference values either obtainedanalytically or provided by the NASA TMR website. In all cases, a satisfactory agreement was observedbetween the solutions of this thesis and the reference values. Moreover, accuracy tests were performedwith mesh refinement, and optimal convergence order was attained for reconstruction orders k = 1 and 3.The k = 2 scheme, however, was slightly inferior in performance compared to its theoretical convergenceorder.Timing results demonstrated the benefit and practicality of using higher-order methods for obtaininga certain level of accuracy. A second-order discretization scheme required a finer grid and morecomputational time to attain the same level of accuracy compared to a higher-order discretizationscheme. Furthermore, strong parallel scaling tests were performed, and excellent scaling was observedfor the Jacobian integration algorithm. The linear solver algorithm also demonstrated satisfactory parallelscaling results, although the performance in this area can still be improved. It was also observed thatthe required number of PTC iterations for convergence is independent of the reconstruction order whileonly slightly affected by the grid size.Much work remains to be done to extend the solution method of this thesis to solve industrially practicalproblems:1. A theoretical study can be conducted to investigate the inferior performance of the k = 2 dis-cretization scheme. The theoretical approach for unstructured mesh-based discretization methodsis usually through the functional analysis framework. A possible starting point can be the workof Barth and Larson [13] where they prove theoretical convergence orders for the k-exact finitevolume method.2. The boundary condition implementations of this work, though functional, might not be the bestchoice. An investigation of other possible boundary condition implementations has the potentialto improve convergence properties and solution quality [29].3. The proposed preconditioning algorithm does not require the explicit form of the full-order LHSmatrix. Thus, a matrix-free method can be implemented, where the matrix-by-vector productof the full-order LHS matrix is evaluated by finite difference Fr\u00C3\u00A9chet derivatives while only the0-exact LHS matrix is assembled. This strategy can reduce the memory cost considerably [55],58while possibly increasing the computation time. Sacrificing computation time for less memoryusage might be unavoidable for very large problems.4. In this thesis, it was possible to directly use the 2-D curvilinear coordinate basis functions for 3-Dproblems, as the solutions of the test problems considered were constant in the x3 direction. Apractical problem, on the other hand, requires the definition of truly 3-D curvilinear coordinatebasis functions. It should be possible to generalize the method of Section 3.3 by providing adefinition for the third wall coordinate. Alternatively, a big hexahedral meta-element can beconstructed by extruding a portion of the nearest wall surface to include all the members of atarget reconstruction stencil. The reference coordinates of the constructed meta-element can thenbe used as the curvilinear coordinates.5. The original 0-exact discretization scheme does not take the gradient dependent source terms intoaccount. Modifying this scheme to somehow capture the influence of the gradient dependentsource terms can greatly enhance the performance of the preconditioning algorithm [78].6. Employing a line Gauss-Seidel preconditioner instead of the ILU method, in conjunction with theproposed inner GMRES preconditioner has the potential to improve the parallel scalability of thelinear solver [4].7. As already mentioned, real world problems contain solution discontinuities, and require shockcapturing methods to be used in conjunction with the k-exact reconstruction scheme. Shockcapturing is even more important for higher-order methods, as they are more prone to developingunstable and spurious solutions.8. Haider et al. [26] recently introduced a multi-level reconstruction scheme that uses compactreconstruction stencils, and eliminates the need to directly work with large stencils. Their high-order finite volume method is therefore better suited for parallel simulations. Although their workis limited to the linear advection equation, its application in the context of turbulent compressibleflows is worth investigating.9. An hp-adaptive strategy can be much more efficient in achieving a prescribed level of accuracy,compared to uniform mesh refinement, or using high reconstruction orders for all the controlvolumes [33, 38]. Desirably, hanging nodes resulting from adaptively refined meshes fit naturallyinto the k-exact finite volume framework. Nevertheless, an effective error estimator must bedesigned that selects the candidate regions for refinement. In the case of anisotropic meshrefinement, the error estimator must also provide an estimate of the direction in which the errorchanges most abruptly.59Bibliography[1] Computational Fluids Laboratory 3-D, 2011. URL https://cfl3d.larc.nasa.gov/.[2] Fully Unstructured Navier\u00E2\u0080\u0093Stokes in 3-D, 2011. URL http://fun3d.larc.nasa.gov/.[3] WestGrid computing resources, 2017. URL https://www.westgrid.ca/.[4] B. R. Ahrabi and D. J. Mavriplis. Scalable solution strategies for stabilized finite-element flowsolvers on unstructured meshes. In Proceedings of the Fifty-Fifth AIAA Aerospace SciencesMeeting, page 517, 2017.[5] B. R. Ahrabi, W. K. Anderson, and J. C. Newman III. An adjoint-based hp-adaptive stabilizedfinite-element method with shock capturing for turbulent flows. Computer Methods in AppliedMechanics and Engineering, 318:1030\u00E2\u0080\u00931065, 2017.[6] S. R. Allmaras, F. T. Johnson, and P. Spalart. Modifications and clarifications for theimplementation of the Spalart-Allmaras turbulence model. In Proceedings of the SeventhInternational Conference on CFD, 2012.[7] P. R. Amestoy, I. S. Duff, and J.-Y. L\u00E2\u0080\u0099Excellent. MUMPS multifrontal massively parallel solverversion 2.0. Technical report, CERFACS: Centre Europ\u00C3\u0083\u00C4\u00BEen de Recherche et de FormationAvanc\u00C3\u0083\u00C4\u00BEe en Calcul Scientifique, 1998. URL http://mumps.enseeiht.fr/.[8] J. Andren, H. Gao, M. Yano, D. Darmofal, C. Ollivier Gooch, and Z. Wang. A comparison ofhigher-order methods on a set of canonical aerodynamics applications. In Proceedings of theTwentieth AIAA Computational Fluid Dynamics Conference, page 3230, 2011.[9] D. N. Arnold. An interior penalty finite element method with discontinuous elements. SIAMJournal on Numerical Analysis, 19(4):742\u00E2\u0080\u0093760, 1982.[10] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout,W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, andH. Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.6, Argonne NationalLaboratory, 2015. URL http://www.mcs.anl.gov/petsc.[11] G. E. Barter and D. L. Darmofal. Shock capturing with PDE-based artificial viscosity forDGFEM: Part I. Formulation. Journal of Computational Physics, 229(5):1810\u00E2\u0080\u00931827, 2010.[12] T. J. Barth and P. O. Frederickson. Higher order solution of the Euler equations on unstructuredgrids using quadratic reconstruction. In Twenty-eighth AIAA Aerospace Sciences Meeting, Jan.1990. AIAA paper 90-0013.60[13] T. J. Barth and M. G. Larson. A posteriori error estimates for higher order Godunov finite volumemethods on unstructured meshes. Finite Volumes for Complex Applications III, London, 2002.[14] F. Bassi and S. Rebay. High-order accurate discontinuous finite element solution of the 2D Eulerequations. Journal of Computational Physics, 138(2):251\u00E2\u0080\u0093285, 1997.[15] J. C. Butcher. Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, 2016.[16] I. B. Celik, U. Ghia, and P. J. Roache. Procedure for estimation and reporting of uncertainty dueto discretization in CFD applications. Journal of Fluids Engineering-Transactions of the ASME,130(7), 2008.[17] M. Ceze and K. J. Fidkowski. Constrained pseudo-transient continuation. International Journalfor Numerical Methods in Engineering, 102(11):1683\u00E2\u0080\u00931703, 2015.[18] S.-W. Cheng, T. K. Dey, and J. Shewchuk. Delaunay Mesh Generation. CRC Press, 2012.[19] X. Deng, M. Mao, G. Tu, H. Zhang, and Y. Zhang. High-order and high accurate CFD methodsand their applications for complex grid problems. Communications in Computational Physics, 11(04):1081\u00E2\u0080\u00931102, 2012.[20] L. T. Diosady and D. L. Darmofal. Preconditioning methods for discontinuous Galerkin solutionsof the Navier\u00E2\u0080\u0093Stokes equations. Journal of Computational Physics, 228(11):3917\u00E2\u0080\u00933935, 2009.[21] K. J. Fidkowski, T. A. Oliver, J. Lu, and D. L. Darmofal. p-Multigrid solution of high-orderdiscontinuous Galerkin discretizations of the compressible Navier\u00E2\u0080\u0093Stokes equations. Journal ofComputational Physics, 207(1):92\u00E2\u0080\u0093113, 2005.[22] J. E. Flaherty. Course notes on finite element analysis. URL http://www.cs.rpi.edu/~flaherje/pdf.Accessed 05/12/17.[23] S. K. Godunov. A finite difference method for the numerical computation of discontinuoussolutions of the equations of fluid dynamics. Matematicheskie Sborniki, 47:357\u00E2\u0080\u0093393, 1959.[24] G. H. Golub and C. F. V. Loan. Matrix Computations. The Johns Hopkins University Press,Baltimore, 1983.[25] F. Haider, J.-P. Croisille, and B. Courbet. Stability analysis of the cell centered finite-volumeMuscl method on unstructured grids. Numerische Mathematik, 113(4):555\u00E2\u0080\u0093600, 2009.[26] F. Haider, P. Brenner, B. Courbet, and J.-P. Croisille. Efficient implementation of high orderreconstruction in finite volume methods. In Finite Volumes for Complex Applications VIProblems & Perspectives, pages 553\u00E2\u0080\u0093560. Springer, 2011.[27] F. Haider, P. Brenner, B. Courbet, and J.-P. Croisille. Parallel implementation of k-exact finitevolume reconstruction on unstructured grids. In High-Order Nonlinear Numerical Schemes forEvolutionary PDEs, pages 59\u00E2\u0080\u009375. Springer, 2014.[28] A. Harten, P. D. Lax, and B. Van Leer. On upstream differencing and Godunov-type schemes forhyperbolic conservation laws. In Upwind and High-Resolution Schemes, pages 53\u00E2\u0080\u009379. Springer,1997.61[29] R. Hartmann and T. Leicht. Generation of unstructured curvilinear grids and high-orderdiscontinuous Galerkin discretization applied to a 3D high-lift configuration. InternationalJournal for Numerical Methods in Fluids, 2016.[30] C. Hirsch. Numerical Computation of Internal and External Flows: Computational Methods forInviscid and Viscous Flows. Wiley, 1990.[31] H. Huynh, Z. J. Wang, and P. E. Vincent. High-order methods for computational fluid dynamics:a brief review of compact differential formulations on unstructured grids. Computers and Fluids,98:209\u00E2\u0080\u0093220, 2014.[32] A. Jalali and C. Ollivier-Gooch. Higher-order unstructured finite volume RANS solution ofturbulent compressible flows. Computers and Fluids, 143:32\u00E2\u0080\u009347, 2017.[33] A. Jalali and C. F. Ollivier Gooch. An hp-adaptive ustructured finite volume solver forcompressible aerodynamic flows. In Proceedings of the Fifty-Fifth AIAA Aerospace SciencesMeeting, page 0082, 2017.[34] A. Jalali, M. Sharbatdar, and C. Ollivier-Gooch. Accuracy assessment of finite volumediscretization schemes on unstructured meshes. Computers and Fluids, 101:220\u00E2\u0080\u0093232, 2014.doi:10.1016/j.compfluid.2014.06.008.[35] C. T. Kelley and D. E. Keyes. Convergence analysis of pseudo-transient continuation. SIAMJournal on Numerical Analysis, 35(2):508\u00E2\u0080\u0093523, 1998.[36] B. S. Kirk, J. W. Peterson, R. H. Stogner, and G. F. Carey. libMesh: a C++ library for paralleladaptive mesh refinement/coarsening simulations. Engineering with Computers, 22(3-4):237\u00E2\u0080\u0093254, 2006.[37] W. M. Lai, D. H. Rubin, D. Rubin, and E. Krempl. Introduction to Continuum Mechanics.Butterworth-Heinemann, 2009.[38] T. Leicht and R. Hartmann. Anisotropic mesh refinement for discontinuous Galerkin methods intwo-dimensional aerodynamic flow simulations. International Journal for Numerical Methods inFluids, 56(11):2111\u00E2\u0080\u00932138, 2008.[39] S. K. Lele. Compact finite difference schemes with spectral-like resolution. Journal ofComputational Physics, 103(1):16\u00E2\u0080\u009342, 1992.[40] W. Li. Efficient Implementation of High-Order Accurate Numerical Methods on UnstructuredGrids. Springer, 2014.[41] H. Luo, J. D. Baum, and R. L\u00C3\u00B6hner. A fast, matrix-free implicit method for compressible flows onunstructured grids. In Sixteenth International Conference on Numerical Methods in FluidDynamics, pages 73\u00E2\u0080\u009378. Springer, 1998.[42] D. Mavriplis, D. Darmofal, D. Keyes, and M. Turner. Petaflops opportunities for the NASAfundamental aeronautics program. In Proceedings of the Eighteenth AIAA Computational FluidDynamics Conference, page 4084, 2007.[43] D. J. Mavriplis. Multigrid strategies for viscous flow solvers on anisotropic unstructured meshes.Journal of Computational Physics, 145(1):141\u00E2\u0080\u0093165, 1998.62[44] D. J. Mavriplis. Directional agglomeration multigrid techniques for high-Reynolds-numberviscous flows. AIAA Journal, 37(10):1222\u00E2\u0080\u00931230, 1999.[45] D. J. Mavriplis. Grid resolution study of a Drag Prediction Workshop configuration using theNSU3D unstructured mesh solver. In Proceedings of the Twenty-Third AIAA AppliedAerodynamics Conference, June 2005. AIAA paper 2005-4729.[46] C. Michalak and C. Ollivier-Gooch. Accuracy preserving limiter for the high-order accuratesolution of the Euler equations. Journal of Computational Physics, 228(23):8693\u00E2\u0080\u00938711, 2009.doi:10.1016/j.jcp.2009.08.021.[47] C. Michalak and C. Ollivier-Gooch. Globalized matrix-explicit Newton-GMRES for thehigh-order accurate solution of the Euler equations. Computers and Fluids, 39:1156\u00E2\u0080\u00931167, 2010.doi:10.1016/j.compfluid.2010.02.008.[48] C. R. Nastase and D. J. Mavriplis. High-order discontinuous Galerkin methods using anhp-multigrid approach. Journal of Computational Physics, 213(1):330\u00E2\u0080\u0093357, 2006.[49] A. Nejat and C. Ollivier-Gooch. A high-order accurate unstructured finite volume Newton-Krylovalgorithm for inviscid compressible flows. Journal of Computational Physics, 227(4):2582\u00E2\u0080\u00932609,2008.[50] H. Nishikawa. Beyond interface gradient: A general principle for constructing diffusion scheme.In Proceedings of the Fortieth AIAA Fluid Dynamics Conference, 2010. AIAA paper 2010-5093.[51] H. Nishikawa. Two ways to extend diffusion schemes to Navier-Stokes schemes: Gradientformula or upwind flux. In Proceedings of the Twentieth AIAA Computational Fluid DynamicsConference, page 3044, 2011.[52] T. Okusanya, D. Darmofal, and J. Peraire. Algebraic multigrid for stabilized finite elementdiscretizations of the Navier\u00E2\u0080\u0093Stokes equations. Computer Methods in Applied Mechanics andEngineering, 193(33):3667\u00E2\u0080\u00933686, 2004.[53] C. F. Ollivier-Gooch. Quasi-ENO schemes for unstructured meshes based on unlimiteddata-dependent least-squares reconstruction. Journal of Computational Physics, 133(1):6\u00E2\u0080\u009317,1997.[54] C. F. Ollivier-Gooch and M. Van Altena. A high-order accurate unstructured mesh finite-volumescheme for the advection-diffusion equation. Journal of Computational Physics, 181(2):729\u00E2\u0080\u0093752,sep 2002. doi:10.1006/jcph.2002.7159.[55] P.-O. Persson and J. Peraire. Newton-GMRES preconditioning for discontinuous Galerkindiscretizations of the Navier\u00E2\u0080\u0093Stokes equations. SIAM Journal on Scientific Computing, 30(6):2709\u00E2\u0080\u00932733, 2008.[56] A. Pueyo and D. W. Zingg. Efficient Newton-Krylov solver for aerodynamic computations. AIAAJournal, 36(11):1991\u00E2\u0080\u00931997, 1998.[57] P. Roe and J. Pike. Efficient construction and utilisation of approximate Riemann solutions. InProceedings of the Sixth International Symposium on Computing Methods in Applied Sciencesand Engineering, VI, pages 499\u00E2\u0080\u0093518. North-Holland Publishing Co., 1985.63[58] P. L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal ofComputational Physics, 43:357\u00E2\u0080\u0093372, 1981. doi:10.1016/0021-9991(81)90128-5. URLhttp://www.sciencedirect.com/science/article/pii/0021999181901285.[59] A. Rogerson and E. Meiburg. A numerical study of the convergence properties of ENO schemes.Journal of Scientific Computing, 5(2):151\u00E2\u0080\u0093167, 1990.[60] C. Rumsey. Turbulence Modeling Resource, 2014. URL http://turbmodels.larc.nasa.gov/.[61] V. Rusanov. Calculation of Interaction of Non-steady Shock Waves with Obstacles. Technicaltranslation. National Research Council of Canada, 1962.[62] Y. Saad. A flexible inner-outer preconditioned GMRES algorithm. SIAM Journal on ScientificComputing, 14(2):461\u00E2\u0080\u0093469, 1993.[63] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2003.[64] K. Salari and P. Knupp. Code verification by the method of manufactured solutions. Technicalreport, Sandia National Labs., Albuquerque, NM (US); Sandia National Labs., Livermore, CA(US), 2000.[65] K. Shahbazi, D. J. Mavriplis, and N. K. Burgess. Multigrid algorithms for high-orderdiscontinuous Galerkin discretizations of the compressible Navier\u00E2\u0080\u0093Stokes equations. Journal ofComputational Physics, 228(21):7917\u00E2\u0080\u00937940, 2009.[66] C.-W. Shu. High order ENO and WENO schemes for computational fluid dynamics. InHigh-order Methods for Computational Physics, pages 439\u00E2\u0080\u0093582. Springer, 1999.[67] P. Solin, K. Segeth, and I. Dolezel. Higher-Order Finite Element Methods. CRC Press, 2003.[68] P. Spalart and S. R. Allmaras. A one-equation turbulence model for aerodynamic flows. InThirtieth AIAA Aerospace Sciences Meeting, Jan. 1992. AIAA paper 92-0439.[69] P. K. Sweby. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAMJournal on Numerical Analysis, 21(5):995\u00E2\u0080\u00931011, Oct. 1984.doi:http://dx.doi.org/10.1137/0721062.[70] E. Toro and A. Chakraborty. The development of a Riemann solver for the steady supersonicEuler equations. The Aeronautical Journal, 98(979):325\u00E2\u0080\u0093339, 1994.[71] T. Toulorge, C. Geuzaine, J.-F. Remacle, and J. Lambrechts. Robust untangling of curvilinearmeshes. Journal of Computational Physics, 254:8\u00E2\u0080\u009326, 2013.[72] B. Van Leer. Towards the ultimate conservative difference scheme. IV. A new approach tonumerical convection. Journal of Computational Physics, 23(3):276\u00E2\u0080\u0093299, 1977.[73] B. Van Leer. Towards the ultimate conservative difference scheme. V. A second-order sequel toGodunov\u00E2\u0080\u0099s method. Journal of Computational Physics, 135(2):229\u00E2\u0080\u0093248, 1997.[74] J. Vassberg, M. DeHaan, and T. Sclafani. Grid generation requirements for accurate dragpredictions based on OVERFLOW calculations. In Proceedings of the Sixteenth AIAAComputational Fluid Dynamics Conference, page 4124, 2003.64[75] M. R. Visbal and D. V. Gaitonde. On the use of higher-order finite-difference schemes oncurvilinear and deforming meshes. Journal of Computational Physics, 181(1):155\u00E2\u0080\u0093185, 2002.[76] J. Von Neumann and R. D. Richtmyer. A method for the numerical calculation of hydrodynamicshocks. Journal of Applied Physics, 21(3):232\u00E2\u0080\u0093237, 1950.[77] M. Wallraff, R. Hartmann, T. Leicht, N. Kroll, C. Hirsch, F. Bassi, C. Johnston, and K. Hillewaert.Multigrid solver algorithms for DG methods and applications to aerodynamic flows. IDIHOM:Industrialization of High-Order Methods-A Top-Down Approach, 128:153\u00E2\u0080\u0093178, 2015.[78] P. Wong and D. W. Zingg. Three-dimensional aerodynamic computations on unstructured gridsusing a Newton\u00E2\u0080\u0093Krylov approach. Computers and Fluids, 37(2):107\u00E2\u0080\u0093120, 2008.[79] P. Woodward and P. Colella. The numerical simulation of two-dimensional fluid flow with strongshocks. Journal of Computational Physics, 54(1):115\u00E2\u0080\u0093173, 1984.[80] S. H. Yoo, H. C. No, H. M. Kim, and E. H. Lee. CFD-assisted scaling methodology andthermal-hydraulic experiment for a single spent fuel assembly. Nuclear Engineering and Design,240(12):4008\u00E2\u0080\u00934020, 2010.65Appendix AANSLib Command-Line OptionsThe command line interface provided by the ANSLib library consists of an executable file Solver, whichcan be executed through the Linux shell, with the input options provided in a separate text file:$ ./Solver -opt < address of options file >ANSLib achieves parallelism using the Message Passing Interface (MPI). Thus, parallel jobs can beexecuted via the MPI launcher :$ mpiexec -np < # of processors > ./Solver -opt < address of options file >In the remainder, the input options are presented for the test problems considered in this thesis. Notethat comments are distinguished by the character #, in the options file.A.1 Two-Dimensional Turbulent Flow over a NACA 0012 Airfoil fromSection 4.5For this problem, several cases of preconditioners were considered. The following options are sharedbetween all the cases:Common Options# --------------------------- MESH-d 2 # Dimensions-mesh_type c # Cell-centered mesh-f < Location of .mesh file without extension >-B < Location of .bdry file without extension >-fec < Location of .fec file without extension ># --------------------------- HISTORY AND SOLUTION FILE NAME-sol < Location of initial solution input file >-sol_out < Location of final solution output file >66-ita_history_name < Location of residual history output file ># --------------------------- ACCURACY-a < Order of accuracy of the scheme: a=k+1 ># --------------------------- PHYSICS-physics RoeTurbSA2D # 2-D RANS + negative S-A-reynolds 6e6 # Reynolds number-mach 0.15 # Mach number-angle 10 # Angle of attack# --------------------------- NON-LINEAR SOLVER OPTIONS-exnut 1 # Ignore \tilde{nu} entries in the line-search-C 0.01 # Initial CFL number-no_distance_weight # No distance weights in the k-exact# reconstruction least squares problem-pseudolts_fixed # No additional CFL ramping-max_iter 80 # Maximum number of PTC iterations-jacobian_type recanalytic # Evaluate the Jacobian exactly-ita_target_residual 1e-5 # Target norm for the non-linear residualTo implement different preconditioning strategies, distinct sets of options have to be provided to the flowsolver in each case:Extra Options for Case A# --------------------------- LINEAR SOLVER OPTIONS# --------------------------- KSP-ksp_max_it 500 # Maximum # of GMRES iterations-ksp_rtol 1e-3 # rtol for the outer KSP-ksp_type gmres # KSP solver type# --------------------------- PC-pc_type ilu # Preconditioner type: ILU-pc_factor_levels 3 # ILU fill level: 3-pc_factor_mat_ordering_type qmd # ILU reordering type: QMDExtra Options for Case B# --------------------------- LINEAR SOLVER OPTIONS# --------------------------- KSP-ksp_max_it 500 # Maximum # of GMRES iterations-ksp_rtol 1e-3 # rtol for the outer KSP-ksp_type gmres # KSP solver type# --------------------------- PC-pre_order 1 # Order of scheme for A*-pc_type ilu # Preconditioner type: ILU-pc_factor_levels 0 # ILU fill level: 0-pc_factor_mat_ordering_type rcm # ILU reordering type: RCM67Extra Options for Case C# --------------------------- LINEAR SOLVER OPTIONS# --------------------------- KSP-ksp_max_it 500 # Maximum # of GMRES iterations-ksp_rtol 1e-3 # rtol for the outer KSP-ksp_type gmres # KSP solver type# --------------------------- PC-mu 1e-5 # \mu_L for line construction-pre_order 1 # Order of scheme for A*-pc_type ilu # Preconditioner type: ILU-pc_factor_levels 0 # ILU fill level: 0-pc_factor_mat_ordering_type lines # ILU reordering type: LINESExtra Options for Case D# --------------------------- LINEAR SOLVER OPTIONS# --------------------------- KSP-ksp_max_it 500 # Maximum # of GMRES iterations-ksp_rtol 1e-3 # rtol for the outer KSP-ksp_type fgmres # KSP solver type: FGMRES# --------------------------- PC-pre_order 1 # Order of scheme for A*-pc_type ksp # Preconditioner type: inner-KSP-ksp_ksp_type gmres # inner-KSP type: GMRES-ksp_ksp_max_it 10 # inner-KSP n_i: 10-ksp_pc_type ilu # inner-KSP preconditioner: ILU-ksp_pc_factor_levels 0 # ILU fill level: 0-ksp_pc_factor_mat_ordering_type rcm # ILU reordering: RCMExtra Options for Case E# --------------------------- LINEAR SOLVER OPTIONS# --------------------------- KSP-ksp_max_it 500 # Maximum # of GMRES iterations-ksp_rtol 1e-3 # rtol for the outer KSP-ksp_type fgmres # KSP solver type: FGMRES# --------------------------- PC-mu 1e-5 # \mu_L-pre_order 1 # Order of scheme for A*-pc_type ksp # Preconditioner type: inner-KSP-ksp_ksp_type gmres # inner-KSP type: GMRES-ksp_ksp_max_it 10 # inner-KSP n_i: 10-ksp_pc_type ilu # inner-KSP preconditioner: ILU-ksp_pc_factor_levels 0 # ILU fill level: 0-ksp_pc_factor_mat_ordering_type lines # ILU reordering: LINES68A.2 Poisson\u00E2\u0080\u0099s Equation from Section 5.1# --------------------------- MESH-d 3 # Dimensions-mesh_type c # Cell-centered mesh-f < Location of .vmesh file without extension ># --------------------------- HISTORY AND SOLUTION FILE NAME-sol_out < Location of final solution output file >-ita_history_name < Location of residual history output file >#---------------------------- PHYSICS-physics Poisson # Solve Poisson\u00E2\u0080\u0099s equation-poisson_problem_data sinhsin # Name of the manufactured solution# --------------------------- ACCURACY-a < Order of accuracy of the scheme: a=k+1 ># --------------------------- NON-LINEAR SOLVER OPTIONS-C 1e10 # Initial CFL number-jacobian_type recanalytic # Evaluate the Jacobian exactly-no_distance_weight # No distance weights in the k-exact# reconstruction least squares problem-max_iter 3 # Maximum number of PTC iterations-ita_target_residual 1e-10 # Target norm for the residual vector R# --------------------------- LINEAR SOLVER OPTIONS-ksp_type gmres # KSP type: GMRES-ksp_rtol 1e-12 # GMRES rtol-pre_order 1 # Order of scheme for A*-pc_type sor # Preconditioner type: PGS-pc_sor_its 10 # Number of PGS iterations: 10A.3 Inviscid Flow Over a Sphere from Section 5.2# --------------------------- MESH-d 3 # Dimensions-mesh_type c # Cell-centered mesh-f < Location of .vmesh file without extension ># --------------------------- HISTORY AND SOLUTION FILE NAME-sol_out < Location of final solution output file >-ita_history_name < Location of residual history output file >----------------------------- PHYSICS-physics Euler3D # Solve 3-D Euler equations-mach 0.38 # Mach number-angle 0 # Angle of attack# --------------------------- ACCURACY-a < Order of accuracy of the scheme: a=k+1 >69-mcell3d_sphere_hack 1 # Curve the boundary of the sphere-mcell3d_extraq_face 1 # Use quadrature rule of order k+1 for faces-mcell3d_extraq_cell 1 # Use quadrature rule of order k+1 for cells# --------------------------- NON-LINEAR SOLVER OPTIONS-C 10 # Initial CFL number-jacobian_type recanalytic # Evaluate the Jacobian exactly-pseudolts_fixed # No additional CFL ramping-max_iter 25 # Maximum number of PTC iterations-ita_target_residual 1e-6 # Target norm for the residual vector R# --------------------------- LINEAR SOLVER OPTIONS-ksp_type fgmres # KSP type: FGMRES-ksp_rtol 1e-3 # GMRES rtol-pre_order 1 # Order of scheme for A*-pc_type sor # Preconditioner type: PGS parallelized by BJ-pc_sor_its 4 # Number of local PGS iterations-pc_sor_lits 4 # Number of BJ iterationsA.4 Turbulent Flow Over a Flat Plate from Section 5.3# --------------------------- MESH-d 3 # Dimensions-mesh_type c # Cell-centered mesh-f < Location of .vmesh file without extension ># --------------------------- HISTORY AND SOLUTION FILE NAME-sol < Location of initial solution input file >-sol_out < Location of final solution output file >-ita_history_name < Location of residual history output file >----------------------------- PHYSICS-physics TurbSA3D # Solve 3-D RANS + negative SA-mach 0.2 # Mach number-reynolds 5e6 # Reynolds number-angle 0 # Angle of attack-turbsa3d_problem_data flat_plate # Specify post-processing operations# --------------------------- ACCURACY-a < Order of accuracy of the scheme: a=k+1 ># --------------------------- NON-LINEAR SOLVER OPTIONS-C 0.1 # Initial CFL number-jacobian_type recanalytic # Evaluate the Jacobian exactly-pseudolts_fixed # No additional CFL ramping-max_iter 25 # Maximum number of PTC iterations-ita_target_residual 1e-8 # Target norm for the residual vector R-no_distance_weight # No weight for the k-exact# reconstruction least squares problem# --------------------------- LINEAR SOLVER OPTIONS-ksp_type fgmres # KSP type: FGMRES70-ksp_rtol 1e-3 # GMRES rtol-pre_order 1 # Order of scheme for A*-pc_type ksp # Preconditioner type: inner-KSP-ksp_ksp_max_it 10 # Inner-KSP n_i: 10-ksp_ksp_type gmres # Inner-KSP type: GMRES-ksp_pc_type bjacobi # Inner-KSP preconditioner: BJ# BJ number of iterations defaults to 1-ksp_sub_pc_type ilu # BJ inner solver: ILU-ksp_sub_pc_factor_levels 2 # ILU fill level: 2-ksp_sub_pc_factor_mat_ordering_type rcm # ILU reordering type: RCMA.5 Turbulent Flow Over an Extruded NACA 0012 Airfoil fromSection 5.4# --------------------------- MESH-d 3 # Dimensions-mesh_type e # Cell-centered extruded mesh-mext_nlayer 7 # Number of extruded layers-mext_btag 2 # Boundary tag for the extruded walls (symmetry)-mext_length 1 # Extrusion length-f < Location of .mesh file without extension >-B < Location of .bdry file without extension >-fec < Location of .fec file without extension ># --------------------------- HISTORY AND SOLUTION FILE NAME-sol < Location of initial solution input file >-sol_out < Location of final solution output file >-ita_history_name < Location of residual history output file >----------------------------- PHYSICS-physics TurbSA3D # Solve 3-D RANS + negative SA-mach 0.15 # Mach number-reynolds 6e6 # Reynolds number-angle 10 # Angle of attack-turbsa3d_problem_data naca0012 # Specify post-processing operations# --------------------------- ACCURACY-a < Order of accuracy of the scheme: a=k+1 ># --------------------------- NON-LINEAR SOLVER OPTIONS-C 0.1 # Initial CFL number-jacobian_type recanalytic # Evaluate the Jacobian exactly-pseudolts_fixed # No additional CFL ramping-max_iter 100 # Maximum number of PTC iterations-ita_target_residual 1e-5 # Target norm for the residual vector R-no_distance_weight # No weight for the k-exact# reconstruction least squares problem# --------------------------- LINEAR SOLVER OPTIONS-ksp_type fgmres # KSP type: FGMRES71-ksp_rtol 1e-3 # GMRES rtol-pre_order 1 # Order of scheme for A*-pc_type ksp # Preconditioner type: inner-KSP-ksp_ksp_max_it 10 # Inner-KSP n_i: 10-ksp_ksp_type gmres # Inner-KSP type: GMRES-ksp_pc_type bjacobi # Inner-KSP preconditioner: BJ# BJ number of iterations defaults to 1-ksp_sub_pc_type ilu # BJ inner solver: ILU-ksp_sub_pc_factor_levels 2 # ILU fill level: 2-ksp_sub_pc_factor_mat_ordering_type rcm # ILU reordering type: RCM72Appendix BSample Script for Running a Parallel Jobon GrexThe computational work on WestGrid systems is carried out through a non-interactive batch job system.Rather than directly executing the solver in the interactive mode, a request has to be made in the form ofa job script, which specifies the commands to be executed as well as the required memory, computationtime, number of nodes, and number of processors. The following shows a sample job script that executesthe solver using 8 processors, and redirects the output to a text file:Sample Job Script#!/bin/bash#PBS -S /bin/bash ## This is a request file for job submission#PBS -l nodes=8:ppn=1 ## Number of nodes and processors per node#PBS -l walltime=00:40:00 ## Maximum requested wall time#PBS -l mem=10gb ## Maximum requested total memory#PBS -M ## Email the owner after the job is finished#PBS -m abe ## Send the email at any occasion## Run an external script that correctly sets the environment## variables and loads Grex\u00E2\u0080\u0099s modules.## For more information on modules see Grex\u00E2\u0080\u0099s website.. ~/scripts/setenv.sh## On many WestGrid systems a variable PBS_NP is automatically## assigned the number of cores requested of the batch system.## On systems where $PBS_NP is not available , one could use:CORES=$(/bin/awk \u00E2\u0080\u0099END {print NR}\u00E2\u0080\u0099 $PBS_NODEFILE)echo \"Running on $CORES cores\"## Announce start timeecho \"Starting run at: $(date)\"73## Recommended MPI options from Petsc for running a parallel job.NPMAX=${CORES}MPI_BINDING=\u00E2\u0080\u0099--bysocket --bind-to-socket --report-bindings\u00E2\u0080\u0099MPIEXEC=mpiexec## Call the ANSLib executable , provide the correct option file,## and finally redirect the output to a file, so that it can## be inspected later.${MPIEXEC} ${MPI_BINDING} -np ${CORES} \~/code/ANSLib/hooshi/apps/solver/Solver \-opt options/h2_p2_facet.opt \|& tee pp_data/h2_p2_facet.out## Announce that the job is finished.echo \"pbs/h2_p2_facet.pbs finished at: $(date)\"After a script has been created, it should be submitted for execution:$ qsub
The batch job systemwill place the script in the execution queue, where the execution priority is evaluatedbased on the the requested amount of resources.74"@en . "Thesis/Dissertation"@en . "2017-11"@en . "10.14288/1.0354983"@en . "eng"@en . "Mechanical Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-ShareAlike 4.0 International"@* . "http://creativecommons.org/licenses/by-nc-sa/4.0/"@* . "Graduate"@en . "A higher-order unstructured finite volume solver for three-dimensional compressible flows"@en . "Text"@en . "http://hdl.handle.net/2429/62846"@en .