UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An adaptive higher-order unstructured finite volume solver for turbulent compressible flows Jalali, Alireza 2017

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

Item Metadata


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

Full Text

An Adaptive Higher-OrderUnstructured Finite VolumeSolver for TurbulentCompressible FlowsbyAlireza JalaliB.Sc., Mechanical Engineering, University of Tehran, 2010M.ASc., Mechanical Engineering, The University of British Columbia, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)January 2017c© Alireza Jalali 2017AbstractThe design of aircraft depends increasingly on the use of Computational Fluid Dy-namics (CFD) in which numerical methods are employed to obtain approximate solu-tions for fluid flows. One route to improve the numerical accuracy of CFD simulationsis higher-order discretization methods. Moreover, finite volume discretizations are themethod of choice in commercial CFD solvers and also in computational aerodynamicsbecause of intrinsic conservative and shock-capturing properties. Considering thatnearly all practical flows with aerodynamic applications are classified as turbulent, wedevelop a higher-order finite volume solver for the Reynolds Averaged Navier-Stokes(RANS) solution of turbulent compressible flows on unstructured meshes.Higher-order flow solvers must account for boundary curvature. Since turbulentflow simulations require anisotropic cells in shear layers, we use an elasticity analogyto project the boundary curvature into the interior faces and prevent faces from in-tersecting near curved boundaries. Furthermore, we improve the accuracy of solutionreconstruction and output quantities on highly anisotropic cells with curvature us-ing a local curvilinear coordinate system. A robust turbulence model for higher-orderdiscretizations is fully coupled to the system of RANS equations and an efficient solu-tion strategy is adopted for the convergence to the steady-state solution. We presentour higher-order results for simple and complicated configurations in two dimensions.These results are verified by comparison against well-established numerical and ex-perimental values in the literature. Our results show the advantages of higher-ordermethods in obtaining a more accurate solution with fewer degrees of freedom andalso fast and efficient convergence to the steady-state solutions.Moreover, we propose an hp-adaptation algorithm for the unstructured finite vol-ume solver based on residual-based and adjoint-based error indicators. In this ap-proach, we enhance the local accuracy of the discretization via h-refinement or p-enrichment based on the smoothness of the solution. Mesh refinement is carried outby local cell division and introducing non-conforming interfaces in the mesh whileiiAbstractorder enrichment is obtained by local increase of the polynomial order in the recon-struction process. Our results show that this strategy leads to accuracy and efficiencyimprovements for several types of compressible flow problems.iiiPrefaceThe research ideas and methods explored and presented in this thesis are the out-come of the research done by Alireza Jalali during his PhD program in the depart-ment of Mechanical Engineering. The implementation of methods, data analysis, andmanuscript preparation were done by the author, Alireza Jalali, under the supervi-sion of Professor Carl Ollivier-Gooch throughout the process. The following are thejournal articles published or in progress based on the contents of this thesis:• A. Jalali, M. Sharbatdar, C. Olliovier-Gooch. “Accuracy analysis of unstruc-tured finite volume discretization schemes for diffusive fluxes”, Computers &Fluids, 101 (2014) 220-232.This paper was co-authored with another PhD candidate, M. Sharbatdar, underthe supervision of C. Ollivier-Gooch. Some results of this paper are available inChapter 4 of this thesis. The code development, data analysis and manuscriptpreparation were done by A. Jalali and an eigen analysis tool developed by M.Sharbatdar for her PhD thesis was used to interpret data.• A. Jalali, C. Ollivier-Gooch. “Higher-order Unstructured Finite Volume RANSSolution of Turbulent Compressible Flows”, Computers & Fluids, 143 (2017)32-47.The results of this paper are also available in Chapters 3 and 4 of this the-sis.• A. Jalali, C. Ollivier-Gooch. “An hp-adaptive Unstructured Finite VolumeSolver for Compressible Flows”, Under review, 2016.A version of chapter 5 of this thesis has been submitted for publication.• M. Sharbatdar, A. Jalali, C. Ollivier-Gooch, “Smoothed Truncation Error inivPrefaceFunctional Error Estimation and Correction using Adjoint Methods in an Un-structured Finite Volume Solver”, Computers & Fluids, 140 (2016) 406-421.This paper is an excerpt from M. Sharbatdar’s PhD thesis about functionalcorrection using adjoint methods. The substantial part of this work includingthe code development for continuous adjoint and solution interpolation tech-niques, and also data analysis and manuscript preparation were done by M.Sharbatdar. A. Jalali provided the discrete adjoint solver which was alreadydeveloped for hp-adaptation method described in Chapter 5 of this thesis.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Turbulence Modeling in CFD . . . . . . . . . . . . . . . . . . . . . . 61.2 Higher-order Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 Solution Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . 152.3 Flux Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3.1 Convective Fluxes . . . . . . . . . . . . . . . . . . . . . . . . 212.3.2 Viscous Fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . 242.4 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.4.1 Flux Integral . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.4.2 Source Term . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.4.3 Moments of Area . . . . . . . . . . . . . . . . . . . . . . . . . 29viTable of Contents2.5 Time Advance Schemes . . . . . . . . . . . . . . . . . . . . . . . . . 302.5.1 Jacobian Matrix . . . . . . . . . . . . . . . . . . . . . . . . . 322.6 Extension to Turbulent Flows . . . . . . . . . . . . . . . . . . . . . . 333 Anisotropic Mesh Treatment . . . . . . . . . . . . . . . . . . . . . . . 373.1 Interior Curving Strategy . . . . . . . . . . . . . . . . . . . . . . . . 373.2 Curved Mesh Support . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2.1 Flux Integral . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2.2 Source Term . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2.3 Moments of Area . . . . . . . . . . . . . . . . . . . . . . . . . 453.2.4 Distance Function . . . . . . . . . . . . . . . . . . . . . . . . 453.3 Reconstruction on Anisotropic Meshes . . . . . . . . . . . . . . . . . 473.3.1 Least-Squares Conditioning . . . . . . . . . . . . . . . . . . . 493.3.2 Curvilinear Coordinates . . . . . . . . . . . . . . . . . . . . . 523.4 Numerical Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.4.1 Straight Meshes . . . . . . . . . . . . . . . . . . . . . . . . . 563.4.2 Curved Meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 613.4.3 General Meshes . . . . . . . . . . . . . . . . . . . . . . . . . 684 RANS Simulation of Turbulent Flows . . . . . . . . . . . . . . . . . 734.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.2 Flux Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.2.1 Convective Fluxes . . . . . . . . . . . . . . . . . . . . . . . . 774.2.2 Viscous Fluxes . . . . . . . . . . . . . . . . . . . . . . . . . . 784.3 Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.4.1 High Reynolds Number Flow Over a Flat Plate . . . . . . . . 844.4.2 Subsonic Flow Over a NACA 0012 Airfoil . . . . . . . . . . . 904.4.3 Transonic Flow Over a RAE 2822 Airfoil . . . . . . . . . . . 934.4.4 High-lift Multi-element Airfoil . . . . . . . . . . . . . . . . . 1005 Adaptive Mesh Refinement and Order Enrichment . . . . . . . . . 1055.1 Error Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1085.2 Adaptation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105.2.1 h-refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . 110viiTable of Contents5.2.2 p-enrichment . . . . . . . . . . . . . . . . . . . . . . . . . . . 1145.2.3 hp-refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . 1155.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1165.3.1 Inviscid Subsonic Flow . . . . . . . . . . . . . . . . . . . . . . 1175.3.2 Inviscid Transonic Flow . . . . . . . . . . . . . . . . . . . . . 1185.3.3 Laminar Subsonic Flow . . . . . . . . . . . . . . . . . . . . . 1275.3.4 Turbulent Subsonic Flow . . . . . . . . . . . . . . . . . . . . 1326 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1386.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1386.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1396.3 Recommended Future Work . . . . . . . . . . . . . . . . . . . . . . . 141Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144AppendicesAppendix A: Cubic Interpolation Functions . . . . . . . . . . . . . . . 159Appendix B: Non-dimensional Equations . . . . . . . . . . . . . . . . 163viiiList of Tables2.1 Number of quadrature points required for source term integration oncell-centered meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.2 Number of quadrature points required for computing moments by in-tegration around the control volumes . . . . . . . . . . . . . . . . . . 302.3 Breakdown of memory requirement for an inviscid subsonic flow problem 343.1 Error in calculating the moments of area for a quarter-annulus . . . . 463.2 Maximum condition number for reconstruction along principal axes onstraight meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.3 Maximum condition number for reconstruction along principal(Cartesian)axes on curved meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 683.4 Maximum condition number for reconstruction along local curvilinearaxes on curved meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 694.1 Friction drag coefficient of turbulent flat plate for different orders andmesh sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.2 Convergence properties of turbulent RANS solver for high Reynoldsnumber flat plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.3 Convergence of drag and lift coefficients with mesh refinement for sub-sonic flow around NACA 0012 . . . . . . . . . . . . . . . . . . . . . 924.4 Convergence properties of turbulent RANS solver for subsonic flowover a NACA 0012 airfoil . . . . . . . . . . . . . . . . . . . . . . . . . 924.5 Convergence properties of turbulent RANS solver for transonic flowover a RAE 2822 airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . 100ixList of Figures1.1 General steps in a CFD simulation . . . . . . . . . . . . . . . . . . . 21.2 Example of mesh types over a 2D airfoil . . . . . . . . . . . . . . . . 42.1 Control volume illustration in 2D . . . . . . . . . . . . . . . . . . . . 132.2 Cell-centered stencil . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3 Solution and gradient reconstruction from two sides for flux integration 212.4 Propagation of a linear wave in positive direction . . . . . . . . . . . 222.5 Schematic illustration of Gauss quadrature points on straight faces . . 262.6 Schematic illustration of higher than second-order quadrature pointson curved boundary faces . . . . . . . . . . . . . . . . . . . . . . . . 272.7 Schematic illustration of quadrature points for fourth-order integrationof source terms over 2D cells . . . . . . . . . . . . . . . . . . . . . . . 282.8 Combination of unstructured finite volume solver elements and re-quired pieces for extension to turbulent flows . . . . . . . . . . . . . . 353.1 Curved boundary intersection with interior faces for anisotropic meshes 383.2 Illustration of cubic reference elements . . . . . . . . . . . . . . . . . 413.3 Schematic representation of boundary displacement for curving theinterior faces of a linear mesh . . . . . . . . . . . . . . . . . . . . . . 423.4 Representation of cubic cells obtained by interior curving of a linearmesh around the geometry of NACA 0012 (dashed lines: linear mesh,solid lines: cubic mesh) . . . . . . . . . . . . . . . . . . . . . . . . . . 423.5 Illustration of mapping from a reference line segment into a generalcubic face . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.6 Illustration of mapping from a reference triangle into an arbitrary cubictriangular cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.7 Illustration of distance function calculation for an arbitrary point . . 47xList of Figures3.8 Anisotropic triangular meshes with uniform stencil: (a) aligned withthe Cartesian coordinate system; (b) rotated 15◦ in the counter clock-wise direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.9 Condition number of reconstruction matrix for uniform stencil meshesaligned with the Cartesian coordinate system . . . . . . . . . . . . . . 513.10 Condition number of reconstruction matrix for uniform stencil meshesaligned rotated 15◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.11 Illustration of tangential-normal coordinate construction . . . . . . . 543.12 Anisotropic manufactured solution for a straight mesh (y-axis has beenscaled) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.13 Accuracy of LS reconstruction for anisotropic straight meshes . . . . 593.14 Reconstructed vs. exact values at different distances from the wall(x = 0.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.15 Anisotropic boundary layer type solution on a circular arc . . . . . . 623.16 Solution reconstruction error using principal coordinate system . . . 633.17 Weighted LS reconstructed values in the principal coordinate systemagainst exact values at different distances from the wall (θ = −1◦) . 643.18 Reconstruction error (relative for normal derivatives) using curvilineart− n coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . 663.19 Reconstructed values in the curvilinear coordinate system against ex-act values at different distances from the wall (θ = −1◦) . . . . . . . 673.20 Unstructured triangular mesh over NACA 0012 . . . . . . . . . . . . 693.21 Separate regions for NACA 0012 meshes . . . . . . . . . . . . . . . . 703.22 Manufactured solution for reconstruction accuracy test on NACA 0012 713.23 Error norms of reconstructed solution on anisotropic cells of NACA 0012 724.1 Anisotropic solution and derivative on a stretched grid with 20×20 cells 804.2 Comparison of eigenvalue spectra for the discretization of viscous fluxeson anisotropic meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . 814.3 Turbulent flat plate test case mesh and geometry . . . . . . . . . . . 854.4 Distribution of wall friction factor and turbulent viscosity for flat plate 874.5 Convergence history of finite volume solver for high-Reynolds turbulentflat plate for n.DoF = 32,640 . . . . . . . . . . . . . . . . . . . . . . 89xiList of Figures4.6 Example of unstructured hybrid meshes used for subsonic flow overNACA 0012 airfoil with 25, 088 degrees of freedom . . . . . . . . . . 904.7 Distribution of scaled turbulence working variable for subsonic flowover NACA 0012 on a mixed-element mesh (n.DoF = 25, 088) . . . . 934.8 Distribution of surface pressure and friction coefficients for second- andfourth-order and comparison with FUN3D . . . . . . . . . . . . . . . 944.9 Convergence history of finite volume solver for subsonic flow overNACA 0012 for n.DoF = 100,352 . . . . . . . . . . . . . . . . . . . . 954.10 Third-order solution of transonic flow over RAE 2822 on a mixed-element mesh (n.DoF = 35, 840) . . . . . . . . . . . . . . . . . . . . . 974.11 Distribution of surface pressure and friction coefficients for differentorders and comparison with experiment . . . . . . . . . . . . . . . . . 984.12 Convergence of drag and lift coefficients with mesh refinement for tran-sonic flow around RAE 2822 . . . . . . . . . . . . . . . . . . . . . . . 994.13 Convergence history of finite volume solver for transonic flow over RAE2822 for n.DoF = 35,840 . . . . . . . . . . . . . . . . . . . . . . . . . 1014.14 Mixed element mesh illustration for the high-lift three-element config-uration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.15 Fourth-order solution of turbulence working variable over multi-elementairfoil on a mixed-element mesh (n.DoF = 45, 802) . . . . . . . . . . 1024.16 Distribution of surface pressure coefficient over multi-element airfoil . 1034.17 Convergence history of finite volume solver for flow over multi-elementairfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045.1 Schematic illustration of h-refinement pattern for 2D cells . . . . . . 1115.2 Schematic illustration of the first rule for the refinement of 2D cells . 1125.3 Schematic illustration of the second rule for the refinement of 2D cells 1135.4 Different number of degrees of freedom for linear elasticity problem onhalf-length non-conforming face . . . . . . . . . . . . . . . . . . . . . 1145.5 Schematic illustration of p-enrichment for an arbitrary cell from second-order to third-order . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1155.6 Convergence of drag coefficient for inviscid subsonic test case (Ma∞ =0.5 and α = 2◦) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119xiiList of Figures5.7 Mesh resolution and discretization order for inviscid subsonic test case(Ma∞ = 0.5 and α = 2◦) . . . . . . . . . . . . . . . . . . . . . . . . . 1205.8 Contours of Mach number for inviscid subsonic test case (Ma∞ = 0.5and α = 2◦) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1215.9 Convergence of drag coefficient for inviscid transonic test case (Ma∞ =0.8 and α = 1.25◦) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1245.10 Final hp-adapted mesh and order for inviscid transonic test case (Ma∞ =0.8 and α = 1.25◦) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1255.11 Contours of Mach number for inviscid transonic test case (Ma∞ = 0.8and α = 1.25◦) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1265.12 Comparison of pressure profiles near the upper surface strong shockbetween adjoint-based hp-adapted meshes for inviscid transonic testcase at y = 0.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1275.13 Comparison of pressure profiles near the upper surface strong shockbetween adjoint-based and residual-based hp-adaptation at y = 0.3 . 1285.14 Convergence of drag coefficient for laminar subsonic test case (Ma∞ =0.5, Re = 5000 and α = 1◦) . . . . . . . . . . . . . . . . . . . . . . . 1305.15 Mesh resolution and discretization order for laminar subsonic test case(Ma∞ = 0.5, Re = 5000 and α = 1◦) . . . . . . . . . . . . . . . . . . . 1315.16 Contours of Mach number for laminar subsonic test case (Ma∞ = 0.5,Re = 5000 and α = 1◦) . . . . . . . . . . . . . . . . . . . . . . . . . . 1325.17 The coarsest quadrilateral mesh for turbulent subsonic flow . . . . . 1335.18 Convergence of lift coefficient for turbulent subsonic test case (Ma∞ =0.15, Re = 6× 106 and α = 10◦) . . . . . . . . . . . . . . . . . . . . 1355.19 Illustration of final hp-adapted mesh for turbulent subsonic flow . . . 1365.20 Contours of turbulence working variable for turbulent subsonic testcase (Ma∞ = 0.15, Re = 6× 106 and α = 10◦) . . . . . . . . . . . . . 137A.1 Reference line segment for face integration . . . . . . . . . . . . . . . 159A.2 Reference cubic triangular element . . . . . . . . . . . . . . . . . . . 160A.3 Reference cubic quadrilateral element . . . . . . . . . . . . . . . . . 161xiiiNomenclatureRoman SymbolsA Areaa Wave speed, Sound speedAR Aspect ratioc Chord lengthCD Drag coefficientCf friction coefficientCL Lift coefficientCp Pressure coefficientcp Specific heat at constant pressureD Destruction, Dissipation matrixd Linear elasticity coefficients, distance functionE Young’s modulusEt Total energyF Flux vectorh Mesh sizeht Total enthalpyI Identity matrixxivNomenclatureJ Jacobian of transformationK Stiffness matrixMa Mach numberP Pressure, Productionp Order of polynomialPr Prandtl numberPrTTurbulent Prandtl numberR Residual operatorRe Reynolds numberS Source term, Vorticity strength, Smoothness indicators Control surface length, Stretching factorT Temperaturet Time, Reconstruction weight exponent, Curve parameterst, n tangential-normal coordinatesU Vector of conserved variables, Left singular vectoru x-velocityV Right singular vectorv y-velocityW Vector of primitive variablesw Reconstruction weightx, y Cartesian coordinatesx′, y′ Principal coordinatesxvNomenclatureZ Adjoint solutionJ FunctionalK Smoothness indicator constantnˆ Unit normal vector~r Reference location connecting vector~V Velocity vectorGreek Symbolsα Angle of attackγ˙ Shear rate tensorǫ Control volume error contributionηl Linear solver toleranceγ Ratio of specific heatsκ Condition number, von Karman constantΛ Diagonalized matrix of eigenvaluesλ Eigenvalues of flux Jacobianµ Dynamic viscosityµT Turbulent viscosityν Poisson’s ratio, Kinematic viscosityω Under-relaxation factorφ Arbitrary primitive variable, Finite element basis functionsρ DensityΣ Diagonal matrix of singular values reciprocalxviNomenclatureσ Stress tensor components, Spalart-Allmaras diffusion constantτ Shear stress tensorθ Azimuthal locationν˜ Turbulence model working variable~δ Displacement vectorξ, η Reference element coordinatesxviiAcknowledgementsFirst and foremost, I would like to acknowledge my research supervisor, Dr. CarlOllivier-Gooch, whose constant support, patience and guidance has been invaluablethroughout the course of this research. I enjoyed a lot during the years I worked withyou as a graduate student and I learned a lot not just in CFD.I would also like to thank my supervisory committee members, Kendal Bushe,BrianWetton and Chen Greif. Your comments and criticisms in the annual committeemeetings and the time and efforts for reading the manuscript have surely led to abetter presentation of this thesis.Also, I am very grateful for the financial support of the University of BritishColumbia through the Four Year Doctoral Fellowship and also the Discovery Grantprovided by the Canadian Natural Sciences and Engineering Research Council. More-over, I would like to thank my dear colleagues at the Advanced Numerical SimulationLab: Reza, Gary, Shayan and Logan. You guys were always available to bounce ideasaround and generously offered me your opinion whenever needed.On a personal side, I would also like to express my sincere appreciation to mywonderful family, especially my dedicated mother, Parvin, and my great father, Mo-hammadhassan, for their love, encouragement and dedication all through my life.Finally, my deepest appreciation goes for my lovely spouse, Mahkame, for her end-less love and support in my life. You have truly been my best companion throughoutthis long journey from the very beginning to this end.xviiiChapter 1IntroductionThe behavior of fluid flows is governed by a set of partial differential equations(PDEs), which relate flow field variables to their derivatives in time and space. Thesolution of these equations enables us to predict the dynamics of the fluid by knowingthe values of flow quantities, such as pressure, temperature and velocity at desiredlocations and times. These PDEs, which are typically obtained by some complicatedphysical models, have an analytic solution only in a limited number of cases derivedby non-trivial mathematical techniques [1]. As a result, we need to use alternativeoptions such as experimental or numerical methods in most of the cases with complexgeometries and/or physical models.Experimental fluid dynamics, which dates back to ancient times, has been usedover the centuries to give insight into flow patterns and measure field quantities. Inrecent years, a considerable amount of efforts have been put to improve the accuracyand efficiency of experimental techniques. However, these methods still suffer frommany restrictions related to their accuracy and reproducibility. There are some errorsthat stem from the disturbances induced by inserting probes into the flow field andreduce the fidelity of experimental data. The measurements are often restricted toone quantity at sampled points every time and thus the measurement of the wholeflow field requires a large amount of time and resources. Nevertheless, there are somecases where once a model is built, a lot of data can be measured. The experimentsprovide prediction only for a laboratory-scale model of simple prototypes whereasmany real life cases (e.g., such as reentry of space vehicles) cannot be investigated inlaboratories since similarity conditions are not achieved.With the recent progress in computational resources and numerical algorithms,Computational Fluid Dynamics (CFD) has emerged as a reliable tool to examinecomplex flows and realistic operating conditions. In this method, which has shown aremarkable ability in quantitative prediction of flow phenomena with high resolutionin time and space, numerical techniques are employed to provide the solution to the1Chapter 1. IntroductionFigure 1.1: General steps in a CFD simulationPDEs describing flow behavior. CFD is also popular due to its lower cost comparedto experiments and its great advantage in computing many flow variables which arenot accessible in an experiment. The ultimate goal of CFD is to find an accurateand reliable solution to fluid dynamic problems in a short amount of time and withminimum computational resources.The numerical simulation of flow fields by means of CFD requires three essentialelements [2]: physical modeling of the flow; domain decomposition, which is referredto as mesh generation; and the numerical approximation of the governing equationsarising from the physical modeling using a robust, efficient and accurate solver. Thesesteps are illustrated in Figure 1.1.The first element includes the mathematical modeling of the flow physics. In2Chapter 1. Introductionthis phase of the process, the system of governing equations is specified along withthe geometry of the domain to be decomposed later and the appropriate physicalboundary conditions. This process is fairly straightforward for simple flow scenarios(e.g., subsonic inviscid or viscous laminar flows) with trivial geometries. However,it exhibits complications for more complex situations such as those with multiphase,turbulence and hydromagnetic field effects for which an exact mathematical descrip-tion of the physical problem is impossible or infeasible, and so physical modelingis required. Therefore, this modeling process leads to a known source of errors inCFD simulations, referred to as physical modeling error. For instance, the turbu-lence models used to predict the behavior of turbulent flows are tuned based on onlya limited number of physical cases and thus fail to provide a general and accurateestimation of turbulence properties in many other applications. Nevertheless, theyseem to be the only feasible approach [3, 4] since the Direct Numerical Simulation(DNS) of turbulent flows will not be practical for the next couple of decades due tocomputing limitations.The next step is to create a mesh on which the solution is approximated. In thisprocess, the physical domain is tessellated with shapes that are recognizable by thesolver. Considering that the governing equations are discretized on the mesh, anappropriate domain decomposition is essential to compute an accurate solution. As aresult, high quality mesh generation is considered as one of the most crucial elementsin CFD.There are basically two different types of meshes: structured and unstructured. Ina structured mesh, all cells and vertices have the same topology and their connectivityis specified by their indices. For example, cell (i, j) is always topologically to the rightof cell (i− 1, j). Such implicit information leads to faster processing by the solverbut restricts the topology of the mesh which in turn makes the automation of gridgeneration challenging for complex geometries. On the other hand, the vertices inunstructured meshes are topologically different and thus their connectivity must bedeclared. These meshes are typically formed as a collection of polygons in 2D andpolyhedrons in 3D. Despite the fact that the processing time of the solver is longer forunstructured meshes, they exhibit a higher flexibility and thus are better candidatesfor arbitrary geometries. The use of unstructured meshes is becoming more popularin modern CFD applications as they promise to be more capable and successful forcomplex aerodynamic problems [5]. Figure 1.2 illustrates the two types of the meshes3Chapter 1. Introduction(a) Structured (b) UnstructuredFigure 1.2: Example of mesh types over a 2D airfoilaround a 2D airfoil.After physical modeling and mesh generation are complete, the infinitely accurategoverning PDEs must be turned into a finite approximation of flow quantities overthe mesh. This process is known as discretization whose error can be expressedas O (hp+1) for problems with smooth solution, where h is a characteristic lengthfor the mesh and p + 1 is the asymptotic order of accuracy. Clearly, the accuracyof a numerical solution can be increased via h-refinement (finer mesh) and/or p-enrichment (higher-order accurate discretization). In the CFD community, thosediscretization schemes where p > 1 (higher than second-order accurate), are typicallyknown as higher-order schemes. The discretization is typically performed by one ofthe following methods: finite difference, finite element and finite volume.The finite difference method is the most traditional method and has historicalimportance. In this method, the point wise quantities are approximated using adifference relation obtained by the Taylor series expansion of the solution on a struc-tured mesh. Finite difference methods can be easily programmed and their extensionto higher-order is straightforward; however, the corresponding solution does not nec-essarily conserve mass, momentum and energy. Also, the method is impractical forunstructured meshes. Due to recent interest in simulating complex geometries, whichis most easily accomplished by the use of unstructured meshes, finite difference is less4Chapter 1. Introductionfrequently used in CFD codes nowadays.The finite element method originated from the structural analysis of solids but isalso applicable to fluid flows. In this method, the solution is represented by local basisfunctions over the elements of the mesh. The degree of the polynomial basis functionis chosen based upon the desired order of solution accuracy. While the solutionrepresented by the polynomial does not completely satisfy the governing equations,the polynomial coefficients are chosen so as to minimize the residuals weighted by atest function. The extension of the finite element method to higher-order solutionscan be easily carried out using a basis function with higher degree. Also, it can be usedfor elements with arbitrary shape and thus is well suited for unstructured meshes.Although the conservation of mass, momentum and energy may be achieved in thefinite element formulation, it is not trivial, particularly in flows with discontinuities.The finite volume method has been designed based on the control volume analysisused for thermodynamics or fluid dynamics systems. The finite volume method dis-cretizes the governing equations in conservative form over cells in the mesh to yieldthe conservation of mass, momentum and energy across the boundaries of control vol-umes. The conservation property enhances the capability of the finite volume methodin capturing discontinuous phenomena such as shock waves. In this technique, con-trol volume averages are used to find the reconstructed piecewise polynomial overeach finite volume. The volume boundary fluxes are computed based on the recon-structed polynomial, and the flux integrals are used to update the control volumeaverages. The finite volume technique, which shows remarkable flexibility for com-plex geometries and unstructured meshes, can be extended to higher-order accuracy[6, 7], although with more difficulty compared to the finite difference or finite elementmethods.In this thesis, we focus on finite volume methods for compressible flows encoun-tered in computational aerodynamics. In particular, we are interested in higher-ordermethods on unstructured meshes for turbulent flows in 2D. In what follows, a briefreview of turbulence modeling approaches employed in CFD simulations is given.Moreover, the state of the art in higher-order CFD methods used in computationalaerodynamics is reviewed. This chapter continues with a clear description of objec-tives set for this thesis and ends with the outline of the following chapters.51.1. Turbulence Modeling in CFD1.1 Turbulence Modeling in CFDTurbulence is frequently observed in nearly all areas of fluid mechanics ranging fromturbomachinery and combustion engines to the design of cars and aircraft. A greatdeal of effort has been put into experimental and theoretical studies to improve theunderstanding of the turbulent motions. However, a universal turbulence theorywhich is able to predict the features of turbulent flows in different applications is notyet available.Numerical simulation of turbulent flows has practical importance in science andengineering. In the context of CFD, there exist three approaches to simulate turbulentflows. These approaches are the Direct Numerical Simulation (DNS), the Large EddySimulation (LES), and the Reynolds-Averaged Navier-Stokes equations (RANS).In the DNS approach, the deterministic Navier-Stokes equations are solved with-out the incorporation of any turbulence model. As a result, the whole range of spatialand temporal scales of the turbulent motions must be resolved to obtain an accurateprediction of the flow field. The smallest turbulence scales, which are referred to asthe Kolmogorov micro-scales, must be captured by the computational length scale(mesh size) and time step. The calculations show that the smallest spatial and tem-poral scales scale with the Reynolds number of turbulent flows as Re−3/4 and Re−1/2,respectively [8]. This induces huge memory and computational power requirementsfor realistic flow problems at relatively high Reynolds numbers (Re > 106). As anexample, the computation of the turbulent flow around an aircraft in one second ofthe flight time would require several thousand years and 1016 grid points using asupercomputer with 1012 Flops [9]. Therefore, the application of the DNS approachhas been limited to simple geometries and flows with low Reynolds number [10, 11].In the DNS simulations, higher-order approximations are typically employed to over-come the numerical dissipation of the lower-order methods [12, 13]. In addition,explicit time integration is mostly used due to the large memory requirements andthe need to resolve small time scales.In LES, the computational cost of DNS is reduced by ignoring the smallest turbu-lence scales. The main idea behind LES is to extract the large scale energy containingcomponents by the convolution of the dependent variables of the Navier-Stokes equa-tions with a predefined low-pass filter, which can be effectively viewed as a spatialand temporal averaging. Therefore, the information about the small scales, which is61.1. Turbulence Modeling in CFDremoved from the numerical solution, must be modeled using a subgrid scale model[8]. In other words, the numerical solution accounts for the geometry dependent largescale eddies of the flow, while the universal smaller scale eddies are captured implic-itly. This separation of scales is possible away from walls where the turbulence isnot in equilibrium but cannot be applied close to solid boundaries where turbulencemanifests in the form of coherent structures [14]. Therefore, an accurate modelingof wall effects is a dominant challenge in LES. Even though LES is computationallyless expensive compared to DNS, it still needs a considerable amount of resources forlarge scale engineering problems with high Reynolds number. In recent years andwith significant progress in computational power, LES is reaching a level of maturityand is becoming a reliable tool for engineering and industrial computations [15]. Theoverall success of LES is dependent on the accuracy of discretization method and theperformance of the simulation in a time- and cost-effective manner [8].The most common approach for the modeling of turbulent flows with engineeringapplications is RANS. In this approach, which is extensively used for steady flowproblems, the physical quantities are decomposed via the Reynolds decompositioninto time-averaged and fluctuating components. The substitution of the decomposedquantities into the Navier-Stokes equations yields the governing equations for themean flow variables. This process leads to the appearance of non-linear Reynoldsstresses which are dependent on the fluctuating velocity components. Additionalmodeling is required for such terms to provide the closure for the system of RANSequations.Most of the turbulence models used in engineering practice employ the Boussinesqassumption in which the Reynolds stresses in the Reynolds-averaged momentum andenergy equations are assumed to be proportional to the mean strain tensor. In thisapproximation, the constant of proportionality is isotropic and is called the turbulentviscosity. The turbulent viscosity can be related to the mean flow quantities, as inalgebraic (zero-equation) models [16, 17] or be obtained by one or more auxiliaryfield equations for the turbulent kinetic energy and time scales [18, 19, 20].On the other hand, the Reynolds stress transport models, which do not use theBoussinesq assumption, are considered as the highest level of available closure forthe system of RANS equations [21]. These models are superior to the turbulentviscosity models because they remove the assumption that the Reynolds stressesrespond immediately to changes in the mean strain rate. Moreover, they take into71.2. Higher-order Methodsaccount the anisotropy of turbulence. However, this comes with the price of solvingadditional PDEs for each component of the Reynolds stress tensor and thus they aremore expensive compared to the turbulent viscosity models in terms of CPU time andmemory usage. In addition, numerical stability problems arise due to the absence ofthe turbulent viscosity and strong coupling [22].Most of the numerical investigations in the field of computational aerodynamicshave been based on steady low-order RANS simulations during the last few decades.This is due to the fact that the second-order RANS simulations became robust and af-fordable on small CPU clusters. Most commercial CFD solvers use RANS turbulencemodels for different applications. However, it is well understood that the RANS mod-els are unable to provide good predictions for flows with separation and/or vortexdominated flows [23]. Considering that such problems are encountered in aerody-namic problems (e.g., landing and take-off with high lift configurations), more pow-erful modeling tools are required. With the recent advances in hardware technology,higher level CFD approaches such as higher-order methods [24], hybrid RANS/LES[25, 26] and LES [27] are being examined in the field of computational aerodynamics.1.2 Higher-order MethodsHigher-order (higher than second-order) discretizations of the fluid dynamic equa-tions have recently received substantial attention due to their potential advantagesin obtaining more accurate solutions with lower cost. With higher-order accuratemethods, the computational cost increases on a fixed mesh; however, a coarser meshcan be used to save time and memory and increase accuracy as well. Several familiesof higher-order methods have been developed over the years in the CFD communityto solve different types of problems. In this section, a short review of previous workusing higher-order methods for the solution of compressible flows is given.Many higher-order methods used in computational aerodynamics have been suc-cessfully developed for computations on structured meshes [28, 29]. Although higher-order unstructured methods are advantageous for complex geometries, it has beenshown that the finite difference method on structured meshes are superior in termsof computational cost and efficiency for boundary layer simulations [30]. Consideringthat structured single-block generation of high quality meshes is challenging for realis-tic problems, higher-order multi-block finite difference solvers have become attractive81.2. Higher-order Methodsin recent years [31]. In these methods, higher-order finite difference approximationsare employed on each block and some type of interface conditions is applied to transferinformation between adjacent blocks [32, 33, 34]. Some other higher-order schemessuch as the Essentially Non-Oscillatory finite volume methods (ENO/WENO) havebeen designed for the solution of hyperbolic equations on multi-dimensional struc-tured meshes. In these methods, which rely on the conservation of flow quantities ineach control volume, the one-dimensional ENO/WENO reconstruction operators areapplied along the coordinate lines to provide a higher-order approximation of numer-ical fluxes across the interfaces. The reconstruction operators use adaptive stencilsso that higher-order accuracy and non-oscillatory properties near discontinuities areobtained [35, 36]. For structured meshes, it has been shown that, for practical levelsof accuracy, using a higher-order accurate method can be more efficient both in termsof solution time and memory usage [37, 6].Research in high-order unstructured solvers is motivated by the accuracy and ef-ficiency advantages demonstrated in the application of these schemes on structuredmeshes with the flexibility of unstructured meshes in adaptivity and for complexgeometries. As described, the finite element approach can be easily implementedon the unstructured meshes and offers a simple path to higher-order discretiza-tions. However, the standard, continuous Galerkin method is inappropriate for useon convection-dominated problems due to its stability issues and shortcomings inshock capturing. In recent years, a class of the finite element methods known asthe Discontinuous Galerkin (DG) schemes have become very popular particularlyin the context of compressible flow simulations. In these methods, the basis func-tions are chosen so that jumps in the solution are allowed across element interfaces.The discontinuous approximation space provides stability for convection problems.The work by Cockburn et al. [38, 39], where DG methods were employed to solvenon-linear time dependent hyperbolic equations, led to the rapid development ofhigher-order DG solvers. The compactness and flexibility of DG methods for higher-order discretizations on unstructured meshes make them attractive even though thenumber of degrees of freedom grows rapidly with polynomial degree, making theseschemes very expensive per element. These methods have been extensively used forthe higher-order solution of compressible Euler [40, 41] and Navier-Stokes [42, 43]equations. Also, recent studies have demonstrated the ability of DG in computingRANS turbulent flows with shock waves [44, 45].91.3. ObjectivesConsidering that the finite volume methods are traditionally the most popularmethod in computational aerodynamics because of intrinsic conservative and shock-capturing properties, the use of higher-order finite volume on unstructured meshes isquite attractive. The first step in the extension of finite volume methods to higher-order discretizations on unstructured meshes dates back to the quadratic reconstruc-tion operator developed by Barth and Frederickson [46] for the Euler Equations.Third- and fourth-order discretizations of these methods have been successfully ap-plied to the Euler and Navier-Stokes equations [47, 7, 48, 49]. These methods haveshown their ability for efficient computations of flows with/without discontinuitieswith promising accuracy on irregular meshes. From the mathematical point of view,finite volume schemes can be written as a Petrov-Galerkin variant of the discontinu-ous Galerkin method regardless of the type of reconstruction operator [50]. The nextmajor step towards the general adoption of higher-order finite volume methods foruse in CFD solvers is turbulence modeling which is the subject of this thesis.1.3 ObjectivesThe higher-order unstructured finite volume solver previously developed at the Ad-vanced Numerical Simulation Lab at UBC is capable of solving inviscid and viscouslaminar compressible flow problems on unstructured meshes. Considering that com-mon types of flow in aerodynamic applications are turbulent, the accuracy that canbe achieved in a flow simulation strongly depends on the prediction of the character-istics of the turbulent flow field. Important physical phenomena, such as boundarylayer separation and shock-boundary layer interaction, can be predicted with highlyaccurate simulation of turbulent flows. The first objective of this thesis is to extendthe capabilities of the current flow solver so that it can compute higher-order RANSsolutions on unstructured meshes over two-dimensional aerodynamic configurations.Having accomplished the first objective, a higher-order unstructured finite vol-ume solver for a wide range of aerodynamic problems from inviscid subsonic to high-Reynolds viscous turbulent will be in hand. However, higher-order accuracy is ad-vantageous only in smooth regions of the solution where there is no discontinuityin the solution. In aerodynamic applications, several sources of discontinuities suchas shocks and contact discontinuities exist which negate the benefit of using higher-order discretizations. Therefore, we can enhance the efficiency (and also robustness)101.4. Outlineof the solver using a so-called hp-adaptation technique. hp-adaptation combines gridrefinement (i.e., h-refinement) and order enrichment (i.e., p-enrichment) into a singleadaptive method where the level and location of refinement are determined based onan a posteriori error estimation. The second objective of this thesis is to develop anhp-adaptive algorithm for our compressible flow solver and apply that on differentclasses of problems which are of particular interest in computational aerodynamics.1.4 OutlineThis thesis is organized as follows.Chapter 2 provides an overview of the finite volume solver that has been usedas the existing infrastructure for the fulfillment of our objectives. Several aspectsof the solver including solution reconstruction, flux functions, integration rules andtime advancement schemes are discussed in detail. Also, the elements required forextension to RANS simulation of compressible aerodynamic flows, which is our firstobjective, are introduced briefly. The details of each of these elements is the subjectof the next two chapters.Chapter 3 explains the treatment of anisotropic meshes considering that cells withvery high aspect ratio are necessary for turbulent flow simulations. These treatmentsinclude curving the interior faces of a linear mesh for higher-order computations,cubic mesh pre-processing algorithms and accurate and well-conditioned solution re-construction on highly anisotropic meshes.In addition, we require a robust and accurate RANS turbulence model for higher-order computations and an effective solution strategy for obtaining the steady-statesolution of the flow field. These elements are described in Chapter 4 along with ournumerical results for a variety of 2D turbulent aerodynamic problems.Chapter 5 is devoted to the hp-adaptive algorithm proposed to accomplish oursecond objective. This chapter describes in detail the error estimation methods usedto drive adaptation, and the mesh refinement and order enrichment techniques em-ployed. Also, the results of our adaptive solver, including its efficiency and accuracyadvantages, are presented for different flow conditions ranging from inviscid to tur-bulent flows.Finally, Chapter 6 summarizes the research in this thesis, provides conclusionsbased on the results, and proposes some future work.11Chapter 2BackgroundThis chapter provides background material relevant for the current research. Thiswork is based on the higher-order unstructured finite volume solver previously devel-oped at the Advanced Numerical Simulation Lab (ANSLib) for inviscid and viscouslaminar flows. The background includes an overview of relevant aspects of this CFDsolver and those elements that are necessary for its extension to turbulent flows.Finite volume methods use the fully conservative control volume form of thegoverning equations that describe fluid flows. These methods compute control volumeaverages of the unknowns over a set of finite volumes that cover the domain. ANSLibis capable of supporting two types of control volumes that are widely used to tessellatea two-dimensional domain: (1) cell-centered type where the cells of the mesh form thecontrol volumes with their centroid as reference locations and (2) vertex-centered typein which the mesh vertices are chosen as the control volumes reference points and theperimeter of each control volume (median-dual) is defined as the lines surroundinga vertex connecting the midpoint of edges to the centroid of the cells. Figure 2.1illustrates the definition of both types of control volumes on a triangular mesh. Inthree-dimensions, ANSLib only supports cell-centered data structure as the formationof vertex-centered control volumes requires a large number of quadrature points percell.To discretize the flow equations using the finite volume method, the governingequations should be recast in fully conservative form as:∂U∂t+∇ · ~F = 0 (2.1)in which U denotes the conserved solution vector and ~F is the flux vector.Integrating Equation 2.1 over an arbitrary control volume in 2D and using thedivergence theorem gives the finite volume formulation of the governing equations inthe form of an area and a surface integral, Equation 2.2, where dA is the infinitesimalarea and nˆ and ds represent the outward unit normal vector and infinitesimal length12Chapter 2. Background(a) Cell-center control volume (b) Vertex-center control volumeFigure 2.1: Control volume illustration in 2Dalong the surface, respectively.ˆ ˆCV∂~U∂tdA+˛CS~F · nˆ ds = 0 (2.2)Assuming the discretized physical domain does not change in time, U can be broughtout from the integral in Equation 2.2, as the average solution vector of the controlvolume:ACVidU idt= −˛CSi~F · nˆ ds (2.3)The left hand side of Equation 2.3 is the time derivative of the average solutionvector in the ith control volume. The right hand side is called the flux integral orresidual of the control volume, which is dependent on the solution vector in generaland represents the spatial discretization of the same control volume:ACVidU idt= −R(U i)(2.4)The residual is evaluated by numerical integration along the interfaces of a controlvolume using sufficient quadrature points.132.1. Governing Equations2.1 Governing EquationsConsidering that the focus of this research is on two-dimensional compressible flowproblems, we turn our attention to 2D descriptions. The conservative form of thecompressible Navier-Stokes equations describing the conservation of mass, momentumand total energy in two dimensions is given as:∂U∂t+∇ ·(~Fc (U)− ~Fv (U,∇U))= 0 (2.5)where U is the conserved solution vector, ~Fc is the convective flux, which depends onsolution only and ~Fv is the viscous flux, which is dependent on both the solution andgradient. For a compressible viscous laminar flow, the solution and flux vectors are:U =ρρuρvEt , Fxc =ρuρu2 + Pρuvu (Et + P ) , Fyc =ρvρuvρv2 + Pv (Et + P )F xv =0τxxτxyuτxx + vτxy + cp(µPr)∂T∂x , Fyv =0τyxτyyuτyx + vτyy + cp(µPr)∂T∂y (2.6)where ρ is the fluid density, ~V = (u, v) are the Cartesian velocity components, P isthe fluid pressure, Et is the total energy, cp is the specific heat at constant pressure,T is the fluid temperature, Pr is the Prandtl number and τij is the viscous stresstensor. Assuming a Newtonian fluid, the viscous stress term becomes:τij = 2µγ˙ij (2.7)γ˙ij =12(∂Vi∂xj+∂Vj∂xi)− 13∂Vk∂xkδij142.2. Solution Reconstructionin which µ is the fluid dynamic viscosity. The fluid pressure can be related to thetotal energy by introducing the ideal gas equation of state given as:P = (γ − 1)[Et − 12ρ(u2 + v2)](2.8)Note that the conserved solution vector is different from primitive variables, W =(ρ u v P )T . We assume that the working fluid is air, γ = 1.4, Pr = 0.72. The viscousflux is zero in the case of an inviscid flow governed by Euler’s equation.2.2 Solution ReconstructionThe spatial accuracy of the finite volume solution depends on the accuracy of the fluxintegral. A high-order accurate flux integral in a control volume requires an accuratenumerical flux and accurate integration. An accurate numerical flux is obtained by anaccurate approximation of the unknown variable in the control volume. Assume thatthe unknown variable represented by φ is one of the primitive variables of the flowfield such as density, pressure, or velocity components. The flow solver approximatesthe unknown variables of the flow field in the control volume by reconstructing apiecewise polynomial about the control volume’s reference point (xi, yi).φRi (x, y) = φ|i +∂φ∂x∣∣∣∣∣i(x− xi) + ∂φ∂y∣∣∣∣∣i(y − yi) + ∂2φ∂x2∣∣∣∣∣i(x− xi)22+∂2φ∂x∂y∣∣∣∣∣i(x− xi) (y − yi) + ∂2φ∂y2∣∣∣∣∣i(y − yi)22+ ... (2.9)In Equation 2.9, φi is the value of the reconstructed solution and∂n+mφi∂xn∂ymare itsderivatives at the reference point of control volume i. These values are the coefficientsof the Taylor polynomial. The degree of the reconstructed polynomial determinesthe order of accuracy of the solution. For example, a second-order accurate solutionapproximation can be achieved by knowing the gradient of the solution at the controlvolume reference point and reconstructing a linear polynomial in the control volume.φRi (x, y) = φ|i +∂φ∂x∣∣∣∣∣i(x− xi) + ∂φ∂y∣∣∣∣∣i(y − yi) +O(∆x2,∆y2)(2.10)152.2. Solution ReconstructionHigher order solutions require the values of higher derivatives at the reference pointand a higher order polynomial.For compressible flow simulations, the hyperbolic system of equations may leadto discontinuities in the solution. Therefore, the reconstruction operator must ensurethe desired order of accuracy while capturing discontinuities such as shocks. Essen-tially non-oscillatory (ENO) finite volume schemes were first developed for structuredmeshes [36] and then extended to irregular unstructured meshes by Abgrall [51] andSonar [52]. In this family of schemes, including weighted ENO (WENO) schemes, thereconstruction stencil is selected based on the solution smoothness and thus stencilscontaining discontinuities are avoided. In this way, the reconstructed solution fallswithin the expected order of accuracy and sharp changes in the solution are capturedautomatically. Nevertheless, these higher-order finite volume schemes suffer fromdifficulties in selecting the appropriate stencil for multi-dimensional problems witha large number of variables. For these problems, the selection of a different stencilfor each flow variable results in complexities and computational costs that limit theirwidespread use in general applications.The construction of a k-exact reconstruction operator is another candidate toreconstruct a solution polynomial based on the control volume data such that thetruncation error of the solution in each control volume remains at the desired orderof accuracy. This method dates back to Barth and Frederickson’s work designing aquadratic reconstruction operator to estimate the advective flux of the Euler equa-tions [46]. Also, Barth incorporated the reconstruction operator into a upwind finitevolume scheme to solve a range of advection-diffusion equations [53]. Later, Ollivier-Gooch and Van Altena [54] described a new k-exact reconstruction procedure for thehigher-order approximation of gradients within a cell which are essential for diffusivefluxes. They demonstrated the high-order accuracy of the reconstruction operatorfor an advection-diffusion problem. In this approach, a fixed stencil is used and aleast-squares system is formed that gives the reconstruction coefficients at the refer-ence point of each control volume. For hyperbolic system of equations, the k-exactreconstruction procedure can be combined with a higher-order limiting strategy tocapture flow discontinuities [55].Our flow solver uses k-exact reconstruction where the coefficients are computedby minimizing the error in predicting the mean value of nearby control volumes [54].162.2. Solution ReconstructionThe conservation of the mean within a control volume requires that1AiˆViφRi dA = φ¯i (2.11)By expanding the left-hand side of Equation 2.11 term by term, one can easily showthat1AiˆViφRi dA = φ|i+∂φ∂x∣∣∣∣∣ixi+∂φ∂y∣∣∣∣∣iyi+∂2φ∂x2∣∣∣∣∣ix2i2+∂2φ∂x∂y∣∣∣∣∣ixyi+∂2φ∂y2∣∣∣∣∣iy2i2+... (2.12)wherexnymi =1AiˆVi(x− xi)n (y − yi)mdA. (2.13)In addition, the error of the mean value of the reconstructed solution for control vol-umes in the stencil {Vj}i should be minimized. In other words, the difference betweenthe actual control volume average φ¯j and the average of φRi over control volume jis minimized. The mean value for a single control volume Vj of the reconstructedsolution φRi is1AjˆVjφRi dA = φ|i +∂φ∂x∣∣∣∣∣i{1AjˆVj(x− xi) dA}+∂φ∂y∣∣∣∣∣i{1AjˆVj(y − yi) dA}+∂2φ∂x2∣∣∣∣∣i{12AjˆVj(x− xi)2 dA}+∂2φ∂x∂y∣∣∣∣∣i{1AjˆVj(x− xi) (y − yi) dA}+∂2φ∂y2∣∣∣∣∣i{12AjˆVj(y − yi)2 dA}+ ... (2.14)To avoid computing moments of each control volume in {Vj}i about the referencepoint of control volume i, replace x − xi and y − yi with (x− xj) + (xj − xi) and172.2. Solution Reconstruction(y − yj) + (yj − yi), respectively. Using Equation 2.13, we obtain1AjˆVjφRi dA = φi +∂φ∂x∣∣∣∣∣i(xj + (xj − xi)) + ∂φ∂y∣∣∣∣∣i(yj + (yj − yi))+∂2φ∂x2∣∣∣∣∣i(x2j + 2xj (xj − xi) + (xj − xi)22)+∂2φ∂x∂y∣∣∣∣∣i(xyj + xj (yj − yi) + (xj − xi) yj + (xj − xi) (yj − yi))+∂2φ∂y2∣∣∣∣∣iy2j + 2yj (yj − yi) + (yj − yi)22+ ... (2.15)The geometric terms in this equation are just dependent on the mesh and can becomputed once and stored. Using the parallel axis theorem, it is possible to showthat these terms have the general form ofx̂nymij ≡1AjˆVj((x− xj) + (xj − xi))n · ((y − yj) + (yj − yi))m dA (2.16)=n∑l=0m∑k=0n!l! (n− l)!m!k! (m− k)! (xj − xi)k · (yj − yi)l · xn−kym−ljHence, Equation. 2.15 can be re-written asφ¯j = φ|i +∂φ∂x∣∣∣∣∣ix̂ij +∂φ∂y∣∣∣∣∣iŷij (2.17)+∂2φ∂x2∣∣∣∣∣ix̂2ij2+∂2φ∂x ∂y∣∣∣∣∣ix̂yij +∂2φ∂y2∣∣∣∣∣iŷ2 ij2+ · · ·This equation is written for every control volume within the stencil of control volumei. The minimum number of neighboring control volumes in the stencil is equal to thenumber of reconstruction coefficients. The control volumes are chosen [54] based ontheir topological proximity to the reconstruction control volume as all neighbors ata given level are added one by one until the number of control volumes in the stencilreaches the number requested for each order of accuracy. The requested number istypically 50%more than the minimum to result in more robust solution reconstructionin the presence of non-smooth and/or vigorously oscillatory data. Figure 2.2 givesthis choice of stencil for an interior control volume labeled R for a cell-centered case;the numeric labels show the order of accuracy at which a control volume is added to182.2. Solution Reconstructionthe stencil. The resulting least-squares problem isR 2232 333444444 4334Figure 2.2: Cell-centered stencil1 xi yi x2i xyi y2i · · ·wi1 wi1x̂i1 wi1ŷi1 wi1x̂2i1 wi1x̂yi1 wi1ŷ2i1 · · ·wi2 wi2x̂i2 wi2ŷi2 wi2x̂2i2 wi2x̂yi2 wi2ŷ2i2 · · ·wi3 wi3x̂i3 wi3ŷi3 wi3x̂2i3 wi3x̂yi3 wi3ŷ2i3 · · ·................... . .wiN wiN x̂iN wiN ŷiN wiN x̂2iN wiN x̂yiN wiN ŷ2iN · · ·φ∂φ∂x∂φ∂y12∂2φ∂x2∂2φ∂x∂y12∂2φ∂y2...i=φ¯iwi1φ¯1wi2φ¯2wi3φ¯3...wiN φ¯N(2.18)where N is the number of nearby control volumes in the stencil and the line separatesthe conservation of mean constraint from equations to be solved by least-squares.In this problem, geometric weights wij can be set to emphasize the importance ofgeometrically nearby data:wij=1|~rj − ~ri|t(2.19)Throughout this thesis, the unweighted LS refers to the case where t = 0 whilethe weighted LS implies t = 1. Note that increasing the value of t in Equation 2.19,192.3. Flux Functionshighlights the importance of geometrically closer cells in the reconstructed polynomialregardless of the order of accuracy. The least-squares system of Equation 2.18 is astandard least-squares problem with the mean constraint. Gaussian elimination isapplied for the constraint, replacing the sub-diagonal entries in the first column withzeros, and the remaining unconstrained (reduced) least-squares problem is solved forevery control volume by the singular value decomposition (SVD) method [56]. In thecase of viscous flow simulations, the Dirichlet boundary conditions pertaining to theno-slip condition at the walls are also added to the least-squares system. Althoughthis can be done by enforcing those conditions as extra hard constraints [54], we justleave them for least-squares minimization with a much larger weight.Having computed the coefficients of the piecewise polynomial, the flux vectorsmust be evaluated at the quadrature points along the interfaces of a control volumeto obtain the numerical flux integral. Considering that the flux vectors depend onthe solution and/or gradient of the solution, the Taylor series expansions must beused to yield such values in an arbitrary quadrature point within the control volume,(xg, yg):φg = φ|i +∂φ∂x∣∣∣∣∣i(xg − xi) + ∂φ∂y∣∣∣∣∣i(yg − yi) + ∂2φ∂x2∣∣∣∣∣i(xg − xi)22+∂2φ∂x∂y∣∣∣∣∣i(xg − xi) (yg − yi) + ∂2φ∂y2∣∣∣∣∣i(yg − yi)22+ ... (2.20)∂φ∂x∣∣∣∣∣g=∂φ∂x∣∣∣∣∣i+∂2φ∂x2∣∣∣∣∣i(xg − xi) + ∂2φ∂x∂y∣∣∣∣∣i(yg − yi) + ... (2.21)∂φ∂y∣∣∣∣∣g=∂φ∂y∣∣∣∣∣i+∂2φ∂x∂y∣∣∣∣∣i(xg − xi) + ∂2φ∂y2∣∣∣∣∣i(yg − yi) + ... (2.22)2.3 Flux FunctionsThe solutions and gradients reconstructed at the two sides of an arbitrary quadraturepoint using Equations 2.20 to 2.22 are not necessarily equal (Figure 2.3). However,the conservation property of finite volume method requires that the flux leaving acontrol volume must enter its neighbor. A numerical flux function is required to202.3. Flux Functionsconstruct a unique flux vector from the two different ones obtained from each side.Flux functions used for the computation of numerical fluxes are divided into twomajor categories: convective (inviscid) and diffusive (viscous) flux functions. Thesetwo groups are characterized based on the physical properties of convection anddiffusion. The transport of properties with the velocity of flow field is known asconvection whereas the mixing or mass transport without requiring bulk motion isreferred to as diffusion. Convective fluxes only depend on the solution while viscousfluxes are also dependent on the gradient of the solution as shown in Equation 2.6.The flux functions used for each category should respect the physical characteristicsof the fluid flow to provide a robust and accurate discretization of the governingequations. Here, we explain the flux functions used in ANSLib for each type for theEuler/Navier-Stokes equations.Figure 2.3: Solution and gradient reconstruction from two sides for flux integration2.3.1 Convective FluxesThe convective (inviscid) part of the governing equations for compressible flow prob-lems is a classic example of a system of non-linear hyperbolic partial differentialequations (PDEs). To have an idea about the appropriate treatment of convectivefluxes, it is helpful to start from the first-order linear wave equation known as the212.3. Flux Functionssimplest hyperbolic PDE:∂φ∂t+ a∂φ∂x= 0 (2.23)This equation governs the propagation of waves traveling at a wave speed a. Forpositive values of a, the wave propagates strictly from left to right along the x axis asshown in Figure 2.4. According to the physical characteristics of the model equation,it is evident that the information in the field is propagating in the wave direction (leftto right). Therefore, the solution at point i is influenced by the solution at point i−1and the solution at point i+1 will not physically affect point i. As a result, in finitedifference methods, it is reasonable to replace ∂φ∂xby φi−φi−1∆xrather than forward orcentral differencing. Such a treatment is called upwind discretization where the in-formation is borrowed from upstream in the direction of “wind”. It is well understoodthat using central or forward differencing can cause non-physical and/or oscillatorybehavior and often leads to instability in the solution as these schemes do not followthe physical characteristics of the governing equation. Upwind discretization in thedirection where the information travels is ideal for other hyperbolic equations.Figure 2.4: Propagation of a linear wave in positive directionThe Euler equations, which form the convective part of the governing equations,are a more complex example of hyperbolic equations due to the non-linearity andvector form. In this case, the information travels along the characteristic lines whose222.3. Flux Functionsslopes are the eigenvalues of the normal flux Jacobian, ∂Fnc∂U:F nc = Fxc nˆx + Fyc nˆy =ρunˆx + ρvnˆy(ρu2 + P ) nˆx + ρuvnˆyρuvnˆx + (ρv2 + P ) nˆyu (Et + P ) nˆx + v (Et + P ) nˆy∂F nc∂U=0 nˆx nˆy 0(γ − 1) qnˆx − uVn Vn − (γ − 2)unˆx unˆy − (γ − 1) vnˆx (γ − 1) nˆx(γ − 1) qnˆy − vVn vnˆx − (γ − 1)unˆy Vn − (γ − 2) vnˆy (γ − 1) nˆy((γ − 1) q − ht)Vn htnˆx − (γ − 1) uVn htnˆy − (γ − 1) vVn γVnq =12(u2 + v2), Vn = u nˆx + v nˆy, ht = h + q (2.24)It is possible to show that the eigenvalues of this matrix are :λi =[Vn Vn Vn + a Vn − a](2.25)where a is the sound speed. These characteristics lines provide the right direction forupwind discretization. Several classes of upwind methods such as flux vector splitting[57, 58], flux difference splitting [59] and advection upstream splitting methods [60, 61]have been proposed to discretize the hyperbolic system of compressible flows.In our solver, we use Roe’s approximate Riemann solver [62] for inviscid fluxdiscretizations. The Roe flux function, which is classified as one of the flux differencesplitting methods, includes the average of the two convective flux vectors computedfrom each side minus a dissipation term which splits (upwinds) the difference of thetwo fluxes:Fc (UL, UR) =12(Fc (UL) + Fc (UR))− 12∣∣∣A˜∣∣∣ (UR − UL) (2.26)The details of the flux difference splitting is very complicated and beyond the scopeof this section. In this formulation,∣∣∣A˜∣∣∣ is the flux Jacobian matrix in diagonalizedform as: ∣∣∣A˜∣∣∣ = X˜−1 ∣∣∣Λ˜∣∣∣ X˜, ∣∣∣Λ˜∣∣∣ = Diag (∣∣∣λ˜i∣∣∣) (2.27)232.3. Flux FunctionsNote that (˜) means evaluations at Roe’s average state defined as:ρ˜ =√ρLρRu˜ =√ρLuL +√ρRuR√ρL +√ρRv˜ =√ρLvL +√ρRvR√ρL +√ρRh˜t =√ρLht,L +√ρRht,R√ρL +√ρR(2.28)2.3.2 Viscous FluxesAs described earlier, the evaluation of diffusive fluxes requires estimates of the so-lution gradient at the quadrature points of a face. The diffusive (viscous) part ofthe governing equations by itself can be considered as a system of elliptic PDEs. Itis well-known that central differencing is an optimal choice for the discretization ofelliptic problems such as Poisson’s equation. Even though for a structured scheme,a central differencing based on the cell averages of neighboring control volumes isusually computed to the desired order [63], this approach is problematic for unstruc-tured grids, because the cell centroids are often far from lying on the perpendicularbisector of the face between them. This drawback decreases the order of accuracy ofthe computed gradient and therefore affects the total accuracy of solvers.Instead, the face gradient can be determined by averaging the two gradients re-constructed from each side of a quadrature point using Equations 2.21 and 2.22. Themost obvious way is arithmetic averaging in which equal weights are set for the twocell gradients; however, the weight of averaging can be tuned by some geometricalconsiderations within the cell. Volume weighted averaging (area weighted in 2D) orlinear interpolation based on the distance of the opposite cell centroid to the facemidpoint can be used as alternatives [64]. In ANSLib, we calculate the interface gra-dients of the primitive variables as the arithmetic average of the two reconstructedgradients:(∇φ)F =12((∇φ)R + (∇φ)L) (2.29)Likewise, the values of those primitive variables involved in the calculation of viscous242.4. Integrationfluxes are obtained by arithmetic averaging:φF =12(φR + φL) (2.30)Using the averaged gradients and values of Equations 2.29 and 2.30, the uniquecomponents of viscous flux vectors in Equation 2.6 can be calculated easily. Notethat this strategy is robust and efficient for low-Reynolds laminar problems wherenearly isotropic meshes are used. In those cases where highly anisotropic meshes arerequired, such a scheme leads to instability as will be shown later.2.4 IntegrationNumerical integration is involved in several parts such as residual evaluation, sourceterm integration and control volume moment calculation. In all of these elements,numerical integration must be accurate enough not to degrade the overall accuracyof the discretization. A summary of quadrature rules and integration method usedfor each one is given in this section.2.4.1 Flux IntegralTo compute the flux integral for each control volume, numerical fluxes should beintegrated over control volume faces. The accuracy of flux integration should beequal or higher than the accuracy of solution reconstruction. In other words, for ap+ 1-order accurate solution, we should be able to integrate a polynomial of degreep accurately. Gauss quadrature gives the capability of evaluating a definite integralwith the integrand evaluated at only a few points [65]. For instance, one quadraturepoint located at the mid-face is sufficient for a linear reconstruction, i.e., second-ordersolution. For higher than second-order, more quadrature points per edge are required.Figure 2.5 shows the quadrature points for the flux integration schematically. It ispossible to verify that these quadrature points are sufficient to obtain higher-orderaccurate flux integrals [66]. Comprehensive information regarding the locations andweights of the Gauss quadrature integration points has been discussed in the work ofStroud and Secrest [65].The control volume flux integral in Equation 2.3 is approximated as the summa-tion of flux integrals over the faces. The flux is integrated along each interior face only252.4. Integration(a) 2nd order/1 Gauss point per face (b) 3rd and 4th order/2 Gauss points perfaceFigure 2.5: Schematic illustration of Gauss quadrature points on straight facesonce and then added to the control volume with outward unit normal vector and sub-tracted from the other one. In this way, we can decrease the cost of flux integrationby a factor of two, as well as ensuring local conservation to machine precision.For those problems with curved boundaries, a piecewise linear representation ofthe boundary is not adequate for higher-order discretizations. This comes from thefact that the line segments between boundary vertices are separated from the actualcurve by a distance that isO (h2) for face length h. In general, the discrete approxima-tion of the boundary must converge to the real boundary shape with the same order ofaccuracy as the order of discretization. For flux integration along curved boundaries,special care must be taken that the locations, unit normals and integration weightsreflect the shape of the boundary appropriately. For general curved boundaries suchas aerodynamic configurations in two-dimensions, we represent the geometry by apiecewise cubic spline which is sufficient up to fourth-order discretization (the dis-tance from the actual curve is O (h4) ). Flux integration along such curves uses theboundary representation directly in setting up the information of Gauss integrationpoints [66]. It should be noted that the locations of the quadrature points are deter-mined based on fractions of arc length, which are calculated iteratively for a cubiccurve. Having the parametric representation of the curved geometry, the unit normalvectors can be obtained. Finally, the weights of integration are assigned based on arclength. Figure 2.6 shows the correct (points Gb1 and Gb2) and wrong (points G′b1 andG′b2) locations of Gauss points along a curved boundary face with their corresponding262.4. Integrationnormals. It is worth mentioning that two quadrature points along curved boundariesare used for flux integration of third- and fourth-order methods.(a) Correct quadrature points (b) Wrong quadrature pointsFigure 2.6: Schematic illustration of higher than second-order quadrature points oncurved boundary faces2.4.2 Source TermEven though the Euler and Navier-Stokes equations do not have a source term, weneed the capability of accurate source term integration for some model problems (e.g.,Poisson’s equation) and also testing manufactured solutions. In the case of havinga source term, S (U,∇U), the residual evaluation includes both the flux and sourceterm integrals on control volumes:R(U i)=˛CS~F · nˆ ds−ˆCVS (U,∇U) dA (2.31)In contrast with the flux integral that is evaluated along the interfaces of a controlvolume, source term integral requires accurate numerical integration over controlvolumes. In this way, the source term integral can be approximated as:ˆCVS (U,∇U) dA ≈Nq∑q=1wq S (xq, yq) ACV (2.32)where Nq represents the number of quadrature points and wq is the correspondingweight of integration for quadrature point q. Table 2.1 shows the number of quadra-272.4. Integrationture points required for each order of accuracy for different types of cell-centeredelements in 2D and 3D. It is possible to verify the nominal accuracy of source termintegration numerically using the reported number of quadrature points. Note thatthe locations of quadrature points for quads and hexes are obtained from the tensorproduct of Gauss points in 1D on a mapped reference element. The details about thelocation of quadrature points and weights of integration for triangles and tetrahedracan be found in Ref. [66]. The quadrature rule used for source term integrationin 2D is exactly valid for flux integration in 3D where the faces of control volumesare triangles and quadrilaterals. Figure 2.7 shows the location of these points in 2Dschematically for a fourth-order scheme.2nd order 3rd order 4th order2Dtriangle 1 3 4quadrilateral 1 4 43Dtetrahedron 1 4 5hex 1 8 8Table 2.1: Number of quadrature points required for source term integration on cell-centered meshes Figure 2.7: Schematic illustration of quadrature points for fourth-order integrationof source terms over 2D cellsFor accurate integration over those boundary cells that are adjacent to a curved282.4. Integrationboundary, we can use the same approach employing ideas from isoparametric finiteelements. However, we choose instead to use an approach for integration over ir-regularly shaped control volumes which guarantees exact integration for polynomialsof degree p (and thus is p + 1-order accurate). Such a strategy has the advantageof reducing the number of quadrature points for vertex-centered meshes in 2D as itremoves the necessity of breaking a vertex-centered control volume into constituentparts. Furthermore, it only requires the calculation of moments of area, which mustbe computed for reconstruction anyway. More information about the implementationof this method for arbitrary control volumes is available in Ref. [66].2.4.3 Moments of AreaThe computation of moments of area, which is necessary for reconstruction, requiresintegration over the surface of each control volumes. To have a high fidelity solutionreconstruction, we need to calculate these moments exactly. It is possible to calculatethe moments using the same rule as source term integration:xnymi =1AiNq∑q=1wq (xq − xi)n (yq − yi)mAi (2.33)However, the integration rule for arbitrarily shaped control volumes such as vertex-centered ones or those adjacent to curved boundaries needs the computation of mo-ments beforehand as described earlier. As a result, we compute the moment by usingthe Gauss’s theorem to convert them into integrals around each control volume:xnymi =1(n + 1)Ai˛CS(x− xi)n+1 (y − yi)m nˆx ds (2.34)where nˆx is the x-component of the outward unit normal vector on the control volumeboundary. In this way, we will have a unified framework for the calculation of mo-ments for all types of control volumes as opposed to computing them using Equation2.33 for straight cell-centered control volumes and using Equation 2.34 for irregu-larly shaped ones. The integral of Equation 2.34 is evaluated by using a Gaussianquadrature of appropriate order along each segment of the surface of each controlvolume as summarized in Table 2.2. The values given in this table are for straightcontrol volume boundaries where the normal vector is constant. For control volumes292.5. Time Advance Schemeswith curved boundaries, three Gauss points are used along the boundary face for allmoments to account for variations in the direction of the unit normal vector.Moment Required for order > No. of Gauss pointsArea 1 1x¯, y¯ 2 2x2, xy ,y2 3 2x3, x2y ,xy2,y3 4 3Table 2.2: Number of quadrature points required for computing moments by integra-tion around the control volumes2.5 Time Advance SchemesFor time-dependent solutions, Equation 2.4 must be integrated in time to yield thesolution at different time levels. This time integration can be done either explicitly orimplicitly. In explicit time advance schemes, the spatial discretization represented bythe residual vector is performed based on the solution at the current time level. Theeasiest explicit time integration scheme is known as explicit Euler where the knownsolution data at the current time level n is used to approximate the time derivativeof the solution vector at control volume i:AiU¯n+1i − U¯ni∆t= −R(U¯ni)(2.35)The solution in the next time level n + 1 is obtained directly by solving the aboveequation. Explicit Euler is only first-order in time and will degrade the overall ac-curacy of time-dependent solutions if combined with higher than first-order spatialdiscretizations. In general, the same order of time integration as the order of spatialdiscretization is required to obtain a high-order time-dependent solution. In ANSLib,we use multi-stage explicit Runge-Kutta schemes with appropriate order for time in-tegration. Equation 2.36 gives a two-stage scheme which is second-order accurate intime and can be used for spatial discretizations up to second-order:AiU¯n+1i − U¯∗i∆t= −R(U¯∗i)(2.36)302.5. Time Advance SchemesAiU¯∗i − U¯ni∆t/2= −R(U¯ni)In this formulation, the intermediate state, U¯∗i , is obtained based on the solutionat current time level n; the solution in the next time level requires another residualevaluation at the intermediate time level. Explicit time advance schemes are easy toprogram and exhibit robustness in many cases. However, they suffer from time-steprestrictions enforced by the Courant-Friedrichs-Lewy (CFL) stability constraint thatresult in taking very small time steps, particularly for higher-order discretizationsand fine meshes. In many engineering applications in which we seek the steady-state solution of the system, taking small time steps prohibitively slows down theconvergence rate.On the other hand, implicit time advance schemes use the next time level solutionfor spatial discretization. Implicit Euler, which is again first-order in time, can bewritten as:AiU¯n+1i − U¯ni∆t= −R(U¯n+1i)(2.37)= −[R(U¯ni)+∂R∂U(U¯n+1i − U¯ni)+O((U¯n+1i − U¯ni)2)]or (I∆t/Ai+∂R∂U)δUi = −R(U¯ni), U¯n+1i = U¯ni + δUi (2.38)where ∂R∂U¯is a large sparse matrix called the flux Jacobian representing the sensitivityof the spatial discretization to the control volume averages exist in the stencil and Iis the identity matrix. The solution at the next time level is obtained by solving thelinear system of Equation 2.38 and computing the update vector, δUi.For steady-state problems, R(U¯i)= 0, we do not care about the temporal ac-curacy of the solution. In this case, one may want to solve the non-linear system ofequations by the direct application of Newton’s methods for steady-state problemsas:∂R∂UδUi = −R(U¯ni), U¯n+1i = U¯ni + δUi (2.39)However, Newton’s method will diverge if the initial guess is too far from the realsolution. As a result, the linear system is augmented by a damping term whichmimics the time derivative in the original time-dependent equations and prevents312.5. Time Advance Schemesthe evolution of non-physical solution at each iteration. This is called the implicitpseudo time-stepping method, which has considerably less strict time-step restrictionscompared to explicit methods and offers a significantly faster convergence to thesteady-state solution.Nevertheless, we must solve the large sparse linear system of Equation 2.38. Thedirect solution of the linear system is computationally intensive, so we use iterativemethods for solving the linear system. In particular, we use the Generalized Min-imum Residual (GMRES) method, which is a preferable choice for non-symmetriclinear systems. GMRES, which approximates the solution by a vector with minimalresidual in the Krylov subspace [67], has shown good performance for complex CFDproblems [68, 69, 47]. The convergence properties of GMRES are highly sensitive tothe conditioning of the matrix. As a result, a good approximation of the left-handside matrix and a robust preconditioning strategy are required. Stationary methodssuch as Gauss Seidel and SOR are easy to implement and they are effective in damp-ing high frequency errors. However, they often have restrictive stability conditionsreducing the benefits of Newton’s method. Alternatively, Incomplete Lower-Upperfactorization methods, ILU(k), have proven to be a robust strategy for GMRES pre-conditioning [70, 7]. The number p shows the fill-level which determines the memoryusage and the accuracy of ILU decomposition; using larger fill-level often leads to amore accurate factorization increasing the quality of preconditioning. However, thereis a restriction in increasing the fill-level in practice due to memory limitations.2.5.1 Jacobian MatrixIn addition, we need to obtain the global Jacobian matrix from our spatial discretiza-tion. This section describes how we calculate this matrix explicitly [49]. As statedbefore, the global Jacobian matrix represents the sensitivity of the residual vectorto the control volume averages in the conservative form. For now, we assume thatthe residual vector only consists of the flux integral. In the case of having a sourceterm that is dependent on the solution and/or gradient, the same procedure can begeneralized. Using the chain rule, the global Jacobian can be expanded as:∂R∂U¯=∂FluxInt∂CVar=∂FluxInt∂Flux∂Flux∂RecSol∂RecSol∂RecCoef∂RecCoef∂PVar∂PVar∂CVar(2.40)322.6. Extension to Turbulent Flowsin which FluxInt is the flux integral, Flux is the numerical normal flux obtained bythe application flux functions, RecSol are the reconstructed solutions and derivativesat Gauss points, RecCoef are the coefficients of reconstructed polynomial obtainedby solving the least-squares system, PVar and CVar are the average of primitiveand conserved variables over control volumes. To compute the Jacobian, each term isfound for each Gauss point along the interfaces of a control volume using the followingprocedure:1. ∂PVar∂CVaris computed for each problem based on the relation between primitiveand conserved variables.2. The∂RecCoef∂PVarterm can be found by using the pseudoinverse of the reconstructionmatrix in the least-squares system of Equation 2.18. The pseudoinverse canbe obtained by the SVD method once and stored since it is only dependent ongeometric terms and does not change between iterations.3. ∂RecSol∂RecCoefterm, which is also a geometric term and only dependent on the locationof the Gauss point, is found using Equations 2.20 to 2.22.4. Flux is defined as the dot product of numerical fluxes and unit normal vectors,~F · nˆ, at each Gauss point. Therefore, ∂Flux∂RecSolis found based on the choices offlux function for inviscid and viscous fluxes described in Section 2.3.5. ∂FluxInt∂Fluxis easily computed as the weight of integration for the correspondingGauss point.Table 2.3 gives the breakdown of memory requirement for a simple subsonic inviscidflow around a multi-element airfoil with zero angle of attack [49]. The values in thetable show the memory requirement in terms of the number of floating point num-bers per control volume for the reconstruction matrix pseudoinverse, global Jacobianmatrix, ILU factorized matrix and Krylov subspace.2.6 Extension to Turbulent FlowsAll the elements described so far lead to a higher-order unstructured finite volumesolver for compressible inviscid and viscous laminar flows which is efficient in com-puting steady-state solutions. Figure 2.8 shows how the major elements of the solver332.6. Extension to Turbulent FlowsOrder Preconditioner Recon. Jacobian ILU matrix Krylov subsp. Total2 ILU(0) 11.7 299.2 299.2 252 862.12 ILU(1) 11.7 299.2 480.3 156 947.23 ILU(0) 88 585 585 172 14303 ILU(1) 88 585 1011.7 124 1808.64 ILU(0) 179.2 615.2 615.2 164 1573.74 ILU(1) 179.2 615.2 1087.6 108 1990Table 2.3: Breakdown of memory requirement for an inviscid subsonic flow problemare combined together to obtain the solution. In the pre-processing stage, an un-structured mesh of arbitrary type (e.g., cell-centered) and a desired physics, whichincludes all the information about flux functions and source term, are used to re-construct an appropriate initial solution field to a desired order. Note that all thosereconstruction elements (e.g., pseudoinverse and geometric terms) that are indepen-dent of solution are created in this stage and stored. The reconstructed solution isused to compute the residual vector and global Jacobian matrix using an appropri-ate quadrature rule. The residual and Jacobian are exported to the implicit pseudotime-stepper to update the solution vector. The new solution vector is brought backto recompute the reconstructed solution and thus residual (fluxes and source term)for the next iteration. This process continues until we converge to the steady-statesolution which satisfies the discrete residual operator to some tolerance.As stated earlier, one objective of this thesis is to extend the capabilities of thecurrent flow solver so that we can compute subsonic and transonic turbulent flows onmixed-element meshes over aerodynamic configurations in two dimensions. For thispurpose, we need several additional pieces highlighted by the gray boxes in Figure2.8.For turbulent flow simulations, it is common to have highly anisotropic cells inthe mesh as the change of properties in one direction is considerably larger thanthe other one in turbulent boundary layer and wake regions. So the first requiredelement is the ability to curve the interior faces of a mesh. As discussed, a high-orderscheme needs a high-order representation of curved surfaces. For isotropic cells orcells with small aspect ratio which are preferable for the Euler and laminar Navier-Stokes equations, we can use the high-order representation in setting up the Gaussintegration points of individual boundary cells with the method described in Section342.6. Extension to Turbulent FlowsFigure 2.8: Combination of unstructured finite volume solver elements and requiredpieces for extension to turbulent flows2.4. For turbulent flows with high aspect ratio cells, on the other hand, this strategyresults in the intersection of boundary face with interior faces and thus non-physicalmeshes as will be shown later. Therefore, we need a mechanism to propagate theboundary curvature to interior faces. The details of this mechanism is described inChapter 3.Having curved the interior faces, we should pre-process the mesh. Consideringthat the median-duals for curved vertex-centered control volumes are difficult evento define uniquely, we choose to add the capability of handling curved meshes (with352.6. Extension to Turbulent Flowscurved interior faces) for cell-centered control volumes. Although there is no differ-ence in the connectivity information between curved and straight cell-centered cases,some other geometric properties (such as moments, quadrature points, etc.) arecomputed differently, as discussed in Chapter 3.It is also well-known that solution reconstruction on high aspect ratio meshessuffers from two known issues: poor conditioning of the least-squares system andpoor accuracy for high-aspect ratio meshes in the presence of curvature. We willshow later that both of these problems originate from the natural choice of Cartesiancoordinate system for reconstruction. To resolve the issues, we need the capabilityof performing accurate solution reconstruction in other coordinate systems such aslocal principal and curvilinear systems. For this purpose, we must be able to findthe mapping coefficients from Cartesian coordinates and also compute the momentsof area in the appropriate coordinate system. The details of such a treatment areexplained in Chapter 3.In addition, we need a RANS turbulence model for the simulation of turbulentflows. In computational aerodynamics, it is quite standard to use either the Spalart-Allmaras one equation turbulence model [19] or some variant of the k − ω two-equation turbulence model [71, 72]. Both of these models have non-linear sourceterms dependent on the solution and gradient of the flow field, which is not the casefor the other governing equations solved by our solver. Therefore, we need to enhancethe ability of Jacobian calculations to include the sensitivity of the source term withrespect to solution using the strategy described in Section 2.5. In this thesis, we use avariant of the Spalart-Allmaras (SA) model designed for higher-order discretizationsand fully couple the evolution equation of turbulence working variable with RANSequations as will be discussed in Chapter 4.The last element that needs to be modified is the implicit solver. For turbulentflow simulations, highly anisotropic cells induce significant stiffness into the discreteequations and have the potential to hamper the solution procedure. This also be-comes worse for higher-order discretizations with smaller numerical dissipation. Inaddition, there is significant difference in the order of magnitude of flow field vari-ables. Consequently, special attention should be taken to the solution strategy forthese problems. This will be done by a robust selection of time-step, appropriateunder-relaxation of the solution update and proper scaling of governing equations aswill be described in Chapter 4.36Chapter 3Anisotropic Mesh TreatmentIn the context of turbulent flows for aerodynamic applications, it is quite commonto have cells with high aspect ratio to resolve the flow properties in one directionas opposed to the direction where quantities change slowly. A typical example isturbulent boundary layer where the gradient of velocity components is quite highacross the boundary layer whereas the gradients are small along the direction of flow.Obtaining sufficient resolution in such regions using isotropic cells leads to additionalunnecessary work including huge computational cost and memory usage due to thelarge number of degrees of freedom.As mentioned earlier, the presence of anisotropic cells in the mesh requires specialtreatment. An effective interior curving strategy, pre-computing geometrical proper-ties of a mesh with curved faces and finally accurate and well-conditioned solutionreconstruction on anisotropic cells are among those elements that need to be consid-ered. This chapter explains each of these in detail.3.1 Interior Curving StrategyThe effect of curvature must be taken into account in higher-order solvers by a moreaccurate representation of boundary faces. When the mesh is comprised of isotropictriangles, it is typically possible to just deform the boundary of the elements whichare in contact with the curved boundary. However, this strategy fails when the cellsare highly anisotropic as is common in boundary layer regions. Figure 3.1 showsthe layers of anisotropic triangles over a curved boundary represented by the dashedred line. It can be seen that even for such a moderate aspect ratio, curving theboundary faces causes intersection with the interior faces and this will be more severefor triangles/quadrilaterals with higher aspect ratio. In this situation, it is necessaryto propagate the mesh deformation into the domain interior to prevent faces fromintersecting near curved boundaries.373.1. Interior Curving Strategy---- Curved BoundaryFigure 3.1: Curved boundary intersection with interior faces for anisotropic meshesTo curve the interior faces of a mesh, Luo [73] suggested a Bezier-based approachto assign curvature to interior faces based on the curvature of the boundary. Also,Sherwin et al. [74] proposed a type of hybrid meshing using quad elements close to thecurved boundary and a curvature based refinement procedure. More recently, Perssonand Peraire [75] devised a method for generating well-shaped curved unstructuredmeshes using a non-linear elasticity analogy. They used the high-order finite elementmethod to solve a non-linear elasticity problem with the boundary curvature actingas an input deformation; the equilibrium displacement from the elasticity problemdefines the curvature of internal faces.In the current thesis, a modified linear elasticity method is used to project theboundary curvature into the interior edges. In this method, which has been used byWang et al. [76], the geometry of the domain to be meshed is considered as an elasticsolid that obeys the isotropic linear elasticity relations. For the initial linear mesh inwhich the Cartesian coordinates are denoted by (X, Y ), the elasticity equation canbe recast in the following form:∂∂X(d11∂δX∂X+ d12∂δY∂Y)+∂∂Y(d33(∂δX∂Y+∂δY∂X))= 0∂∂X(d33(∂δX∂Y+∂δY∂X))+∂∂Y(d21∂δX∂X+ d22∂δY∂Y)= 0 (3.1)δ = (δX , δY ) represents the nodal displacement vector in the Cartesian coordinate383.1. Interior Curving Strategydirection and the coefficients, d, are defined as follows:d11 = d22 =E (1− ν)(1 + ν) (1− 2ν)d12 = d21 =Eν(1 + ν) (1− 2ν) (3.2)d33 =E2 (1 + ν)in which E is assumed to be a constant throughout the domain and thus canceledand ν denotes Poisson’s ratio, which is set to be 0.25. Equation 3.1 can be re-writtenin the following form:∇ · S = 0 (3.3)where S represents the stress tensor.S = σXX σXYσY X σY Y (3.4)For the finite element discretization, Equation 3.3 is multiplied by an arbitrary testfunction, z (X, Y ), and integrated over the domain A. Integration by parts gives theweak formulation of the system of equations:ˆAS (δX , δY ) · ∇z (X, Y ) dA = 0 (3.5)where the boundary term is not present because the test functions vanish on bound-aries due to essential Dirichlet boundary conditions. Following the standard finiteelement technique, the discrete displacement vector is represented at the nodes usingnodal basis functions, φ (X, Y ):~δ (x, y) =∑k~δkφk (X, Y ) (3.6)Using the Galerkin finite element method, the test function is set equal to the nodalbasis function at each node and thus the system of Equation 3.5 can be written asKδ = f where the stiffness matrix entries are 2× 2 matrices:393.1. Interior Curving StrategyKij =´A(d11∂φi∂X∂φj∂X+ d33∂φi∂Y∂φj∂Y)dA´A(d12∂φi∂X∂φj∂Y+ d33∂φi∂Y∂φj∂X)dA´A(d33∂φi∂X∂φj∂Y+ d21∂φi∂Y∂φj∂X)dA´A(d33∂φi∂X∂φj∂X+ d22∂φi∂Y∂φj∂Y)dA (3.7)To compute the stiffness matrix, nodal basis functions are found for each element upto the desired order of accuracy. For this purpose, Lagrange cubic (10-node) elementsfor triangles and serendipity cubic (12-node) elements for quadrilaterals are usedto represent the boundary geometry up to fourth-order accuracy. The serendipitycubic (12-node) element is chosen instead of a tensor-product cubic (16-node) elementto decrease the cost of solving the linear elasticity problem. Figure 3.2 shows thereference element for both cases in the reference space of (ξ, η). The nodal basisfunctions for both the cubic reference elements are given in Appendix A. Note that thederivatives of basis functions in the physical space,(∂φ∂X, ∂φ∂Y), are needed in Equation3.7 rather than those in the reference space,(∂φ∂ξ, ∂φ∂η). By knowing the mapping fromthe reference element to any of the physical cells, (ξ, η)→ (X, Y ), which is linear fortriangles and bilinear for quadrilaterals, these derivatives can be computed at anyarbitrary point by using:∂φ∂X∂φ∂Y =∂X∂ξ∂Y∂ξ∂X∂η∂Y∂η−1 ∂φ∂ξ∂φ∂η (3.8)The integrals of Equation 3.7 are computed using sixth-order quadrature rules [65]by evaluating Equation 3.8 in quadrature points of the reference elements, (ξq, ηq).The computed elemental stiffness matrices are assembled into a global matrix to formthe left-hand side of the system of equations.On the boundary of the domain, the displacement can be obtained for the twomiddle nodes of each boundary face by representing the geometry by a higher-ordercurve (Figure 3.3). The location of these nodes in the linear mesh is found by linearinterpolation between the two end points of the corresponding face:Xk = Xa + kXb −Xa3, Yk = Ya + kYb − Ya3for k = 1, 2 (3.9)403.1. Interior Curving Strategy(a) Triangle (b) QuadrilateralFigure 3.2: Illustration of cubic reference elementsUsing piecewise cubic splines for the geometry ((xB (t) , yB (t)), where t = [0, 1] foreach curved boundary face), the two middle nodes are assumed to be displaced tothe locations obtained by interpolation in the parametric space of the same boundaryface, (xB,i (t) , yB,i (t)):xk = xB,i(t =k3), yk = yB,i(t =k3)for k = 1, 2 (3.10)Therefore, the boundary displacement which is imposed as a Dirichlet boundarycondition is found as δi = (xi −Xi, yi − Yi). Solving the system of equations arisingfrom Equation 3.1 gives the displacement, and thus the location of nodes in the curvedmesh. The linear system is solved via GMRES method preconditioned by ILU(1).Note that the Cartesian coordinates are represented by (x, y) in the higher-ordercurved mesh and by (X, Y ) in the initial linear mesh. Figure 3.4 shows higher-orderelements near the leading edge of the NACA 0012 geometry obtained by curvingthe initial linear mesh shown by dashed lines. As seen, all the interior faces arecurved in this way to prevent the intersection of boundary and interior. Far from thecurved surface, the deformations are considerably smaller although the cells are stillrepresented by cubic faces.413.2. Curved Mesh SupportFigure 3.3: Schematic representation of boundary displacement for curving the inte-rior faces of a linear meshFigure 3.4: Representation of cubic cells obtained by interior curving of a linearmesh around the geometry of NACA 0012 (dashed lines: linear mesh, solid lines:cubic mesh)3.2 Curved Mesh SupportAfter interior curving, the cells of the mesh are not triangles/quadrilaterals withstraight faces. Therefore, special care should be taken to compute those geometri-cal features required for a finite volume solver. These include flux integral quadra-ture points on the interfaces along with corresponding normal vector and integrationweight, and quadrature points and weights required for source term integration andmoments of area, which are essential for solution reconstruction.3.2.1 Flux IntegralA face of a higher-order mesh is comprised of four nodes coming directly from solvingthe linear elasticity equations at the nodes of each cubic element. Therefore, a cubicmapping from a reference straight line segment into each curved face (Figure 3.5) can423.2. Curved Mesh Supportbe created as:x (ξ) =4∑k=1xkφk (ξ) , y (ξ) =4∑k=1ykφk (ξ) (3.11)where (xi, yi) are the coordinates of the nodes and φi (ξ) are the Lagrangian cubicinterpolation functions in 1D (see Appendix A). Having this, we can also map fluxintegration along curved faces into flux integration along the straight reference linesegment as:˛CS~F · nˆ ds =˛CS~F · nˆ√dx2 + dy2 =˛CS~F · nˆ√√√√(dxdξ)2+(dydξ)2dξ (3.12)To evaluate the integral of Equation 3.12, we can use appropriate Gauss quadraturerules for the reference line segment. However, the flux vector must be evaluated atthe quadrature points mapped to the physical space whose locations are obtainedfrom Equation 3.11. The unit normal vectors at the quadrature points are foundas (−dy/dξ, dx/dξ). Due to the changes in the direction of unit normal vector, weuse a quadrature rule that is able to integrate a polynomial of degree 2p+ 1 exactlyon the reference line segment for a (p+ 1)-accurate discretization. Our numericalexperiments showed that this number of quadrature points is always sufficient.Figure 3.5: Illustration of mapping from a reference line segment into a general cubicface3.2.2 Source TermTo compute the source term integral over a curved cell, we use the idea of isoparamet-ric finite elements in which element geometry and displacement vector componentsare represented by the same type of basis functions. In other words, the Cartesiancoordinates inside each cubic triangle/quadrilateral can be interpolated by the same433.2. Curved Mesh Supportcubic basis functions used to solve the linear elasticity problem:x (ξ, η) =N∑k=1xkφk (ξ, η) , y (ξ) =N∑k=1ykφk (ξ, η) (3.13)Note that the number of nodes, N , is 10 for a cubic triangle and 12 for a cu-bic serendipity quadrilateral. Figure 3.6 shows such a mapping for cubic trianglesschematically. To calculate the integral of a source term over a curved cell, we canuse coordinate transformation to the reference element space as:ˆCVS (U,∇U) dA =ˆCVS (U,∇U) dxdy =ˆCVS (U,∇U)∣∣∣∣∣∂ (x, y)∂ (ξ, η)∣∣∣∣∣ dξdη (3.14)in which∣∣∣∂(x,y)∂(ξ,η)∣∣∣ = | J | denotes the determinant for the Jacobian of the coordinatetransformation:| J | =∣∣∣∣∣∂x∂ξ ∂y∂η − ∂y∂ξ ∂x∂η∣∣∣∣∣ (3.15)Therefore, the integral of Equation 3.14 can be evaluated using a quadrature rulethat is able to integrate a polynomial of degree 2p exactly over the reference element.Note that the source term, S (U,∇U), must be evaluated at the quadrature pointsmapped to the physical space whose coordinates are found using Equation 3.13.Figure 3.6: Illustration of mapping from a reference triangle into an arbitrary cubictriangular cell443.2. Curved Mesh Support3.2.3 Moments of AreaAs described in Chapter 2, the least-squares solution reconstruction requires themoments of each control volume about its reference point. For this purpose, we canuse the same strategy as source term integral since the moments of area are definedas volume integrals:xnymi =1AiˆCVi(x− xi)n (y − yi)m| J | dξdη (3.16)where the (xi, yi) is the coordinate of the reference point of a control volume. Thispoint is obtained by mapping point(13, 13)for triangles and point (0, 0) for quadrilat-erals from the reference element to the physical space using Equation 3.13. To obtainsufficiently accurate moments, we use a quadrature rule that is able to integrate amonomial of degree 2p + 1 exactly on the reference element for a (p+ 1)-accuratereconstruction.To assess the accuracy of curved boundary representation and moment calcula-tions, the moments of area of a quarter-annulus, whose inner and outer radii are 7/8and 9/8, are found about the origin and compared with the exact values. For thispurpose, the moments of each element are computed first and then translated to theorigin by the use of the parallel axis theorem. The error involved in computing themoments has been tabulated in Table 3.1 for two different N values which representthe number of divisions in both the radial and azimuthal directions. The ratio oferror observed for different moments verifies the fourth-order accuracy of the curvedboundary representation and also sufficient accuracy of numerical integration up tofourth order.3.2.4 Distance FunctionAnother geometrical property that must be computed for anisotropic curved meshesis the distance of the cell reference location from the closest curved boundary. Thisquantity is needed for accurate solution reconstruction on high aspect ratio meshesnear curved geometries as will be shown in the next section. In addition, the sourceterm of the Spalart-Allmaras turbulence model is dependent on the distance fromthe closest wall. As a result, we need to pre-compute and store the distance from theclosest wall boundary for all of the quadrature points used to compute the integral453.2. Curved Mesh SupportMoment|E|N |E|20|E|40N = 20 N = 40A 7.37× 10−8 4.61× 10−9 15.97x 5.71× 10−8 3.66× 10−9 15.60y 5.71× 10−8 3.66× 10−9 15.60x2 1.04× 10−7 6.24× 10−9 16.61xy 4.80× 10−8 3.38× 10−9 14.19y2 1.03× 10−7 6.23× 10−9 16.48x3 1.34× 10−7 8.05× 10−9 16.69x2y 5.67× 10−8 3.70× 10−9 15.31xy2 5.62× 10−8 3.70× 10−9 15.19y3 1.33× 10−7 8.05× 10−9 16.49Table 3.1: Error in calculating the moments of area for a quarter-annulusof the source term. In this section, the method by which the distance from a higher-order curve is calculated for an arbitrary point (reference location and quadraturepoints) is described.As discussed, we represent curved geometries with a piecewise cubic spline. Notethat the cubic polynomial representing boundary face i (one piece of the spline) isdifferent from other boundary faces:xB,i (t) = α3,it3 + α2,it2 + α1,it+ α0,i (3.17)yB,i (t) = β3,it3 + β2,it2 + β1,it+ β0,iTo find the closest point to an arbitrary point of (x0, y0) on the curved geometry,we need to identify the boundary face on which the closest point exists (Figure 3.7).For this purpose, we compute the unit normal vectors at the two end points of eachboundary face (nˆ1 and nˆ2). It is worth-mentioning that the first derivative and thusunit normal vector of a cubic spline are continuous at the points. We also find thevectors that connect each of the two end points to point (x0, y0). Computing thecross product of this vector, ~ri , and corresponding unit normal vector, nˆi, the closestpoint lies on the boundary face where the two cross products have different signs,i.e., ((~r1 × nˆ1) · (~r2 × nˆ2)) < 0. Note that the cross product of two vectors in two-dimensions is a scalar. Having found the appropriate boundary faces, the distance463.3. Reconstruction on Anisotropic Meshessquared of each point on that face from (x0, y0) is:D2 =(α3,it3 + α2,it2 + α1,it+ α0,i − x0)2+(β3,it3 + β2,it2 + β1,it+ β0,i − y0)2(3.18)The closest point is obtained as the one with minimum distance squared whose cor-responding parameter satisfies dD2dt= 0. This minimization problem is solved forparameter tc via Newton’s method to machine zero tolerance. The correspondingparameter is substituted in Equation 3.17 to find the Cartesian coordinates of theclosest point on the curved geometry. If we do not find any boundary face where thecross product changes sign, the closest point on the surfaces will be one of the cornerpoints of the cubic spline (e.g., trailing edge). The distances from all of the cornerpoints are measured and compared to find the closest point on the wall.Figure 3.7: Illustration of distance function calculation for an arbitrary point3.3 Reconstruction on Anisotropic MeshesA known issue for k-exact solution reconstruction on highly anisotropic meshes is thepoor conditioning of the least-squares (LS) system. This is exacerbated for higher-order computations as the presence of geometric terms with higher exponents resultsin a considerable difference in the magnitude of matrix columns. Furthermore, aero-dynamic configurations typically consist of curved surfaces on which a flow boundarylayer develops. Previous studies have demonstrated that least-squares reconstruction473.3. Reconstruction on Anisotropic Mesheson high aspect ratio meshes with finite curvature suffers from poor accuracy, evenwith straight-sided cells. Mavriplis [77] examined the accuracy of the least-squarestechnique for second-order discretization on unstructured meshes. He used a non-linear function of distance to the wall as a test function to evaluate the accuracy ofthe reconstruction on highly stretched meshes in the presence of surface curvature.The results showed that unweighted least-squares reconstruction fails to estimate thesolution gradient; however, the accuracy can be recovered using an inverse distanceweighting for vertex-centered discretizations. For a cell-centered method on triangu-lar meshes, none of the control volumes within the stencil may be sufficiently closeto the cell center under consideration. Therefore, inverse-distance weighting does nothelp and least-squares reconstruction exhibits poor accuracy. To improve the accu-racy of reconstructed gradient on the curved meshes with high aspect ratio, Diskinet al. [78] proposed a least-squares minimization in a mapped domain using a localcurvilinear coordinate system aligned with the wall. The mapped domain is con-structed using the distance function, which is the nearest distance to the boundary,to estimate the solution and gradient at the face mid-points. Although this methodproduces more accurate gradient for second-order schemes, it does not satisfy con-servation of mean within each control volume nor can it be used for higher-ordermethods. Petrovskaya introduced the concept of numerically distant points to ex-plain the poor accuracy of a weighted LS method on stretched meshes [79]. Shealso studied quadratic least-squares reconstruction on triangular cells with high as-pect ratio and concluded that the weighting of distant stencil points does not leadto more accurate approximation because of numerically distant points in the stencil[80]. More recently, Petrovskaya [81] designed a new approach that allows one tomeasure the distance between points in the data space instead of the physical spacefor second-order accuracy and devised a new weighting function that de-emphasizesthe points which are far from the original point in the data space.In this and the next sections, we will investigate the accuracy of high-order so-lution reconstruction on highly anisotropic meshes with or without curvature forcell-centered discretization along with the conditioning of the least-squares problems.We will show that the condition number of the LS system grows rapidly with meshrefinement and increasing aspect ratio even on meshes without curvature. We willsuggest some modification for solving the LS system that improves the condition-ing considerably. In addition, our numerical results reinforce the notion that the483.3. Reconstruction on Anisotropic Meshescell-centered least-squares reconstruction in its traditional form degrades solution ac-curacy on high aspect ratio meshes with finite curvature in the sense that the errorin the reconstructed coefficients are relatively large compared to the exact values anddo not attain the expected order of accuracy. While this has already been proven forsecond-order methods on triangles, we will demonstrate that the issue remains forhigher-order LS reconstruction and adding a weight function does not improve theaccuracy. Consequently, a new least-squares reconstruction framework is proposedin which the solution and derivatives are approximated in a local tangential-normalcoordinate system obtained by solving an auxiliary LS problem. We will show thatthe k-exact reconstruction in the new coordinate system improves the accuracy andconditioning significantly for high aspect ratio cells with finite curvature and satisfiesthe expected order of error reduction with mesh refinement.3.3.1 Least-Squares ConditioningThe reduced least-squares system of Equation 2.18 can be written in the form ofAx = b (3.19)where A is the reconstruction matrix after one step of the Gaussian elimination, xis the reconstruction coefficients vector and b is the weighted/unweighted differenceof control volume averages of the solution within the stencil. Given an SVD of thereconstruction matrix, the pseudoinverse A† can be obtainedA† = VΣ†UT (3.20)in which the columns of U and V are the left and right singular vectors of A andΣ† is a diagonal matrix containing the reciprocal of singular values of A. Since thereconstruction matrix only includes geometric weights, it does not change betweeniterations [49]. Therefore, its pseudoinverse is precomputed and stored to yield thereconstruction coefficients by a matrix-vector multiplicationx = A†b (3.21)To have an accurate estimate of the reconstruction coefficients, the least-squares493.3. Reconstruction on Anisotropic Meshesproblem should be well-conditioned in the sense that a small error in b (the controlvolume averages in this case) does not cause a large error in x. This property is typi-cally investigated by the condition number of the reconstruction matrix A (equivalentto the condition number of A†) which is defined as the ratio of the maximum to theminimum singular value. The condition number bounds the relative inaccuracy in theestimate of the unknown components. A least-squares problem with a large conditionnumber is said to be ill-conditioned.-0.5 0.5-0.20.2(a)(b)hFigure 3.8: Anisotropic triangular meshes with uniform stencil: (a) aligned with theCartesian coordinate system; (b) rotated 15◦ in the counter clock-wise directionAs a preliminary test, we consider the least-squares reconstruction on anisotropictriangles obtained from similar splitting of quadrilaterals as shown in Figure 3.8.We employ the unweighted LS method to investigate the effect of aspect ratio andmesh size on the condition number. Note that in such a mesh, the stencil is uniformthroughout all interior triangles and thus all control volumes have the same conditionnumber. The columns of the reconstruction matrix of Equation 2.18 contain themoments of area that are scaled by both the size and aspect ratio of the triangularcells. Increasing the aspect ratio of the cells and/or decreasing the mesh spacingresults in a huge difference in the order of magnitude of matrix columns that inturn leads to large condition numbers. This difference is asymptotically larger forhigher order reconstruction whose matrix includes higher moments of area. Thisbrief discussion suggests the idea of matrix scaling with the largest entry of eachcolumn [82] which also was used by Ivan and Groth to improve the conditioningof a central ENO reconstruction for viscous flows [83]. In this way, all the entriesof the reconstruction matrix are in the range of [−1, 1] and the conditioning of theleast-squares problem improves. Also, the entry with the maximum absolute value503.3. Reconstruction on Anisotropic Meshesis stored to recover the right reconstruction coefficients later. Figure 3.9 shows howthe condition number of the unweighted reconstruction matrix changes with aspectratio and mesh spacing with/without column scaling for the mesh that is alignedwith the Cartesian coordinate system. Our results show that the condition numberof matrix A is proportional to ARk−1/hk−2 for a k-th order reconstruction. However,the column scaling makes the conditioning independent of the mesh size and aspectratio.ARCondition Number101 102 103101103105107109101110131015 unscaled-2nd orderunscaled-3rd orderunscaled-4th orderscaled-2nd orderscaled-3rd orderscaled-4th orderSlope = 2.0Slope = 3.0Slope = 1.0(a) Effect of AR, h = 0.1hCondition Number0.04 0.08 0.1210110310510710910111013Slope = 1.0Slope = 2.0(b) Effect of h, AR = 1000Figure 3.9: Condition number of reconstruction matrix for uniform stencil meshesaligned with the Cartesian coordinate systemNow, we turn our attention to the mesh being rotated 15◦ in Figure 3.8. For thiscase, the condition number has been plotted in Figure 3.10 for the same variationsof aspect ratios and mesh sizes. As seen, the condition number of the scaled recon-struction matrix still grows with the aspect ratio in a way analogous to the unscaledsystem although the values are smaller by several orders of magnitude. This is incontrast with the case where the cells were aligned with the Cartesian coordinatesystem and column scaling led to independence from both the aspect ratio and meshsize. This can be explained by the fact that the Cartesian coordinate system is al-most aligned with the principal axes1 as the aspect ratio increases in that mesh. Sothe combination of column scaling and principal coordinate system makes the mesh1Principal axes are the two directions around which the second moments of area are minimumand maximum and obtained via Mohr’s circle in 2D513.3. Reconstruction on Anisotropic Meshesisotropic and resolves the conditioning issue for high-order reconstruction on highaspect ratio meshes. It is worth mentioning that this happens just for meshes wherethe principal axes are almost parallel for the cells within the stencil (Figure 3.8).Nevertheless, for more general cases such as anisotropic meshes over curved surfaceswhere the orientation of principal axes changes across the stencil, we will show thatthe combination of column scaling and principal coordinate system is still helpful inreducing the condition number.ARCondition Number101 102 103101103105107109101110131015 unscaled-2nd orderunscaled-3rd orderunscaled-4th orderscaled-2nd orderscaled-3rd orderscaled-4th orderSlope = 2.0Slope = 3.0Slope = 1.0(a) Effect of AR, h = 0.1hCondition Number0.04 0.08 0.1210110310510710910111013Slope = 1.0Slope = 2.0(b) Effect of h, AR = 1000Figure 3.10: Condition number of reconstruction matrix for uniform stencil meshesaligned rotated 15◦3.3.2 Curvilinear CoordinatesAs mentioned earlier, the simulation of high Reynolds number turbulent flows requiressufficiently accurate polynomial approximation on high aspect ratio meshes with finitecurvature. It is well understood that the cell-centered k-exact LS reconstructionsuffers from poor accuracy on highly anisotropic meshes over curved surfaces forsecond-order. Our numerical results in Section 3.4 yield the same conclusion forhigher than second-order reconstruction even after weighting the least-squares systemby geometrically close data. This problem is alleviated for those cells that are farfrom the wall and thus have smaller curvature.523.3. Reconstruction on Anisotropic MeshesTo improve the accuracy of the cell-centered reconstructed gradient on curvedmeshes with high aspect ratio, the local mapping proposed by Diskin et al. [78] isquite inspiring. However, this method is only applicable to second-order accurate re-construction as it only gives the point-wise first-order face gradient in the curvilinearcoordinate system constructed at the middle of each face. This also takes the advan-tage of the fact that the solution values at the cell centroids are equal to the controlvolume averages for second-order and the reconstruction matrix is only comprisedof ∆x and ∆y. For higher-order methods, this scheme cannot be incorporated withthe moments of area and real control volume averages that differ from the point-wisevalues at the centroid. Therefore, a higher-order polynomial cannot be reconstructedwithin a cell to give the solution and the gradient vector.Instead, we construct a higher-order local curvilinear coordinate system at thereference point of each control volume. For this purpose, we define a mapping fromthe physical space into a tangential-normal coordinate system for cells with highcurvature near the walls:t = a1 (x− xi) + a2 (y − yi) + a3 (x− xi)2 + a4 (x− xi) (y − yi) + ... (3.22)n = b1 (x− xi) + b2 (y − yi) + b3 (x− xi)2 + b4 (x− xi) (y − yi) + ...A cubic mapping from (x, y) to (t, n) is sufficient for reconstruction up to fourth-order. The mapping is obtained by the distance function and constructed tangentialdirection at the cell’s reference points. The difference in distance from the walldetermines the normal coordinate, n(i)j = Dj − Di, while the tangent coordinateis obtained by the projection of the vector connecting two reference points on theconstructed tangential direction, t(i)j = ~rij · tˆi. This tangential direction is definedas the perpendicular direction to the normal to the wall direction as seen in Figure3.11. The two sides of Equation 3.22 are evaluated for the same handful of controlvolumes used in the reconstruction stencil of a particular control volume. To find themapping coefficients, an auxiliary least-squares system is solved to give the values ofai and bi in Equation 3.22.It should be noted that the auxiliary least-squares problem will suffer from thesame conditioning issue explained before since the moments of area in the reconstruc-tion matrix are replaced by some other geometric terms that have the same orderof magnitude and are scaled in an unsatisfactory manner with mesh refinement and533.3. Reconstruction on Anisotropic MeshesFigure 3.11: Illustration of tangential-normal coordinate constructionincreasing aspect ratio. Consequently, the auxiliary LS problem is also brought tothe principal coordinate system denoted by (x′, y′) and its origin is fixed on the cell’sreference point (Figure 3.11). The new mapping is defined from (x′, y′) to (t, n) ast = a′1x′ + a′2y′ + a′3x′2 + a′4x′y′ + ... (3.23)n = b′1x′ + b′2y′ + b′3x′2 + b′4x′y′ + ...and column scaling is applied to improve the conditioning.Having computed the mapping coefficients, it is possible to find the coordinatesof any arbitrary point with respect to the local tangential-normal system of a certaincontrol volume by applying the curvilinear mapping of Equation 3.23. For solutionreconstruction in the new curvilinear coordinates, we require the moments of controlvolume i about its reference point:tnnmi =1AiˆVitnnmdA (3.24)The integration is performed over the curved cells using the method described inSection 3.2.3. For this purpose, the location of each quadrature point in the physi-cal space (xq, yq) is obtained by Equation 3.13 with its corresponding Jacobian oftransformation Then, the appropriate rotation transformation is applied to map(xq − xi, yq − yi) to(x′q, y′q)and finally the non-linear mapping of Equation 3.23 isemployed to find (tq, nq). This procedure can also be extended to compute the mo-543.4. Numerical Testsments of control volume j about the reference point of control volume i, i.e., t̂nnmij .For these moments, the quadrature points of control volume j must be combined withthe curvilinear mapping coefficients of control volume i. Contrary to the reconstruc-tion in the Cartesian coordinate system, these moments need to be calculated one byone since the parallel axis theorem cannot be used for non-parallel local curvilinearcoordinates.The final point about the curvilinear coordinate system is the transformation ofsolution (or derivatives) to an arbitrary location in the physical space, (xp, yp). Asdescribed for the quadrature points, it is possible to find the curvilinear coordinatesof the desired point, (tp, np). Therefore, the solution can be found as:up (t, n) = u|i +∂u∂t∣∣∣∣∣itp +∂u∂n∣∣∣∣∣inp +∂2u∂t2∣∣∣∣∣it2p2+∂2u∂t∂n∣∣∣∣∣itpnp +∂2u∂n2∣∣∣∣∣in2p2+ ... (3.25)Similarly, the gradient at any location can be computed in the local coordinate sys-tem, (∂u/∂t, ∂u/∂n)p. One can transform the gradient to the Cartesian coordinatesystem by applying the change of variables and the chain rule:∂u∂x∣∣∣∣∣p=∂u∂t∣∣∣∣∣p.∂t∂x′∣∣∣∣∣p.∂x′∂x∣∣∣∣∣p+∂u∂t∣∣∣∣∣p.∂t∂y′∣∣∣∣∣p.∂y′∂x∣∣∣∣∣p(3.26)+∂u∂n∣∣∣∣∣p.∂n∂x′∣∣∣∣∣p.∂x′∂x∣∣∣∣∣p+∂u∂n∣∣∣∣∣p.∂n∂y′∣∣∣∣∣p.∂y′∂x∣∣∣∣∣p∂u∂y∣∣∣∣∣p=∂u∂t∣∣∣∣∣p.∂t∂x′∣∣∣∣∣p.∂x′∂y∣∣∣∣∣p+∂u∂t∣∣∣∣∣p.∂t∂y′∣∣∣∣∣p.∂y′∂y∣∣∣∣∣p+∂u∂n∣∣∣∣∣p.∂n∂x′∣∣∣∣∣p.∂x′∂y∣∣∣∣∣p+∂u∂n∣∣∣∣∣p.∂n∂y′∣∣∣∣∣p.∂y′∂y∣∣∣∣∣p3.4 Numerical TestsIn this section, we present the results for k-exact least-squares reconstruction onanisotropic meshes. First, we show the accuracy and conditioning results correspond-ing to unweighted and weighted LS reconstructions of an empirical flat plate boundary553.4. Numerical Testslayer profile on straight meshes where the variation of principal axis orientation issmall between adjacent triangles. Then, we consider anisotropic triangular meshesover curved surfaces. For this purpose, we generate anisotropic triangular meshesover a circular arc and reconstruct a boundary layer type function. The reconstruc-tion procedure both in its traditional form but along principal axes and on the newcurvilinear coordinate system is performed and the results are compared. Then, weextend the analysis to more general meshes by testing reconstruction scenarios onthe anisotropic boundary layer region over the NACA 0012 airfoil.3.4.1 Straight MeshesWe start the analysis with anisotropic straight meshes where the alignment changesslowly when moving from one triangle to another. These meshes are typically used foranisotropic flows over non-curved surfaces on which a boundary layer develops (e.g.,flat plate) and/or wake regions behind single objects. Anisotropic meshes over curvedsurfaces can fall in this category provided that the radius of curvature is considerablylarger than the cell spacing.To investigate the accuracy of reconstructed derivatives, the method of manufac-tured solution is employed. Considering that anisotropic meshes are used to capturesolution anisotropy, an anisotropic function must be manufactured for reconstructiontests. One of the obvious candidates for such a function is the turbulent boundarylayer velocity profile. In this paper, we make use of Reichardt’s empirical boundarylayer profile [84] :u ={2.5 ln(1 + 0.4 y+)+ 7.8(1− exp(−y+11)− y+11exp(−0.33 y+))}u∗ (3.27)where y+ is the non-dimensional distance from the wall and u∗ is the friction velocity:y+ =yu∗νu∗ =√τwρ(3.28)In Equation 3.28, τw is the wall shear stress which can be related to the friction factorCf,x =τw12ρU2∞(3.29)563.4. Numerical Testsand approximated by the one-seventh power law:Cf,x =0.027(Rex)1/7, Rex =U∞xν(3.30)Using Equations 3.28 to 3.30, it is possible to simplify u∗ and y+ :u∗ =√√√√ 0.0135(ReL.x¯)1/7, y+ = y¯.ReL.√√√√ 0.0135(ReL.x¯)1/7(3.31)where x¯ and y¯ are the non-dimensional coordinates based on the length of plate, L,for a plate whose leading edge is at the origin. We assume that L = 1 and takeReL = 107 to provide highly anisotropic behavior for u¯ = u/U∞ near the wall. Notethat u∗ and y+ are both singular at x = 0 and thus the region close to the leadingedge must be avoided. For our reconstruction tests, we assume that the leading edgeis at (−1, 0) and only consider the zone between x = −0.5 and x = 0.5. Therefore, x¯is replaced with x+1 in Equation 3.31. The grid is generated with (N+1)×(2N + 1)nodes which are uniformly distributed along the plate but stretched out from the wallwith a factor of s. The thickness of the first layer is set such that one control volumeexists in the viscous sublayer (y+ < 5) for the coarsest mesh; this is shrunk by afactor of two at each level of refinement. The quadrilaterals formed by these nodesare randomly divided into two triangles and nodes are perturbed in both directionswith a factor of rhy where r is a random number in [−0.25, 0.25] and hy is localvertical spacing. We consider three different mesh sizes in which N = 8, 16, 32 ands = 1.6, 1.25, 1.12, respectively. It is worth mentioning that for all these meshes, themaximum aspect ratio is about 10, 000. This value is obtained by considering thenumber of divisions in the horizontal direction and the necessity of having one cell inthe viscous sublayer region for the coarsest mesh. Figure 3.12 shows the anisotropicsolution on a randomly triangulated mesh where N = 16 both in a far view and veryclose to the wall.To assess the accuracy of reconstruction, the L2-norm of error in the reconstructedsolution and first derivatives at the reference point of each triangle is considered forboth unweighted and weighted least-squares systems. In all these cases, the re-construction has been performed along the principal axes of each control volume toimprove the conditioning of the least-squares problem. In this framework, the deriva-573.4. Numerical TestsX-0.4 -0.2 0 0.2 0.400.020.040.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1-0.4 -0.2 0 0.2 0.402E-054E-05Figure 3.12: Anisotropic manufactured solution for a straight mesh (y-axis has beenscaled)tives have been brought to the Cartesian coordinate system and compared against theexact derivatives obtained from the manufactured solution. Figure 3.13 shows asymp-totic error convergence for second-, third- and fourth-order solution reconstruction. Inall cases, the order of error in the reconstructed solution matches the expected orderof reconstruction. Likewise, the reconstructed derivatives are one order less accuratethan the solution, as expected. Note that the y−derivatives are larger due to theanisotropic property of the test function and thus their relative errors are reported tobe comparable with the x−derivatives. As seen in this figure, the difference in solu-tion error is less noticeable between weighted (red) and unweighted (black) for second-and fourth-order but is about one order of magnitude for third-order reconstruction.Moreover, the unweighted LS yields more accurate derivatives particularly in thex−direction for all orders of accuracy. These results suggest that the unweighted LSprovides more accurate reconstructed values for anisotropic straight meshes.In addition to error norms, the local accuracy of reconstruction is important.Figure 3.14 compares the solutions and y-derivatives obtained by second- and fourth-order unweighted LS reconstruction with the exact values at different distances fromthe wall at x = 0.1. As expected, the reconstructed values match well with theexact profiles even at small distances in which the triangles are highly skewed. Also,583.4. Numerical TestshL 20.05 0.1 0.1510-410-310-210-1100101u (unweighted LS)u/ x (unweighted LS)u/ y (unweighted LS)u (weighted LS)u/ x (weighted LS)u/ y (weighted LS) Second-orderhL 20.05 0.1 0.1510-510-410-310-210-11003.351.952.093.302.26(b) Third-orderhL 20.05 0.1 0.1510-610-510-410-310-210-11003.723.182.942.89(c) Fourth-orderFigure 3.13: Accuracy of LS reconstruction for anisotropic straight meshes593.4. Numerical Teststhe error is reduced by increasing the order of accuracy. This can be more clearlynoticed for the y-derivatives that are one order less accurate and more scattered inmagnitude.u ( x = 0.1 , y)y0 0.2 0.4 0.6 0.8 110-610-510-410-310-210-1Exact2nd-order4th-order(a) Solutionu/ yy0 5000 10000 1500010-610-510-410-310-210-1Exact2nd-order4th-order(b) y-derivativeFigure 3.14: Reconstructed vs. exact values at different distances from the wall(x = 0.1)To address the conditioning of the LS reconstruction for straight meshes, weconsider the maximum condition number, κ∞ (A), corresponding to the least-squaressystems. This is given by Table 3.2 for different conditions. In this table, ARmaxrefers to the maximum aspect ratio which occurs in the first layer. Increasing ReLin the solution reduces the thickness of viscous sub-layer and thus the first interiorlayer which increases the maximum aspect ratio in the mesh. Clearly, column scalingimproves the conditioning of the LS system in all cases. Also, the unweighted systemafter column scaling exhibits remarkably better conditioning for higher than second-order as it leads to considerably smaller condition numbers for similar meshes andthe condition number remains almost constant when changing the maximum aspectratio and/or mesh size. On the other hand, the weighted system even after columnscaling produces condition numbers proportional to aspect ratio for third- and fourth-order reconstructions. This result implies more accurate reconstruction at extremelyanisotropic meshes using unweighted LS combined with column scaling.603.4. Numerical TestsARmax h Col. Scaleκ∞ (A)k = 2 k = 3 k = 4unweighted LS3.15× 103 0.0625 yes 2.21 13.78 107.129.65× 103 0.125 yes 2.00 16.69 107.469.65× 103 0.0625 yes 2.20 17.35 106.779.65× 103 0.0625 no 1.43× 104 1.44× 1010 6.19× 10159.65× 103 0.03125 yes 2.39 15.55 104.712.67× 104 0.0625 yes 2.21 15.53 95.522.67× 104 0.0625 no 4.07× 104 1.26× 1011 1.39× 1017weighted LS3.15× 103 0.0625 yes 2.74 1.41× 104 6.48× 1043.15× 103 0.03125 yes 2.76 1.43× 104 6.83× 1049.65× 103 0.0625 yes 2.78 4.28× 104 2.07× 1059.65× 103 0.0625 no 8.10× 103 5.96× 109 2.67× 10152.67× 104 0.03125 yes 2.76 1.27× 105 6.20× 1052.67× 104 0.03125 no 2.48× 104 1.17× 1011 3.26× 1017Table 3.2: Maximum condition number for reconstruction along principal axes onstraight meshes3.4.2 Curved MeshesNow, we focus on meshes varying anisotropically over curved surfaces. For thesemeshes, the alignment varies considerably between adjacent triangles. This is usuallyencountered in high-Reynolds turbulent flow simulations over curved walls on whicha boundary layer exists and/or wake regions of upstream objects in multi-elementconfigurations.As the first step for testing different reconstruction scenarios, we generate anisotropiccells over a circular arc and reconstruct a boundary layer type function. For this pur-pose, we use the boundary layer profile of Equation 3.27 and replace y with distancefrom the wall and x with the horizontal distance from the leading edge. Even thoughthis function does not give the velocity profile over a cylinder with stream-wise pres-sure gradient, it has similar characteristics and exhibits strong normal gradients, andvanishing stream-wise gradients. To prevent singularities given by Equation 3.27613.4. Numerical Testsat the cylinder’s leading edge, we consider a 40◦ circular arc whose initial edge issufficiently away from the leading edge, θ = −π/2 (Figure 3.15). Again, we use(N + 1) × (2N + 1) nodes with the aforementioned values for N and diagonalizequadrilaterals randomly (in the case of triangular meshes) along with a node pertur-bation based on local radial spacing. Note that the nodes are placed with uniformspacing along the arc but are stretched towards the wall with the stretching factorsalready mentioned for straight meshes. We still ensure the presence of at least onecontrol volume in the region where y+ based on distance from the wall is less than5 and set ReL = 107. Figure 3.15 shows the anisotropic solution and mesh whereARmax ≃ 6320, N = 16 and s = 1.25.-0.4 0.40.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1Axis of symmetryFigure 3.15: Anisotropic boundary layer type solution on a circular arcHaving manufactured an appropriate solution, we examine the accuracy of higher-order reconstruction in its traditional form described in Chapter 2 but using theprincipal coordinate system to make the conditioning better. For the results shownin this section, the boundary faces on the circular arc are curved using the methoddescribed in Section 3.1 which is accurate up to fourth-order. Figure 3.16 shows theL2 norm of solution reconstruction error for the manufactured function of Figure 3.15and the three mesh sizes described earlier for both the triangular and quadrilateralmeshes.For triangles, it is seen that the asymptotic order of error in the unweightedLS is not consistent with the order of reconstruction in the range of mesh sizesconsidered here. Even though there is a chance of achieving the order expected inreconstruction with more levels of refinement, this negates the advantage of higher-623.4. Numerical Testsorder computations. Weighting the least-squares system corrects the asymptoticorder of error but at the same time is associated with a noticeable increase in thevalue of error. Comparing the error values with the solution magnitude illustrated inFigure 3.15 implies a large relative error which is not acceptable.For quadrilateral meshes, the unweighted LS delivers the right asymptotic orderof convergence. However, increasing the order of accuracy degrades the accuracyof reconstructed solution at the cell reference points. Likewise, adding the distanceweight to the LS system increases the magnitudes of reconstruction error. Therefore,it is accurate to say that traditional higher-order least-squares reconstruction suffersfrom accuracy issues for anisotropic meshes over curved surface regardless of theweighting function and cell type.hL 2  (u e - u)0.02 0.04 0.06 0.08 0.110-410-310-210-1100101102 2nd order - unweighted LS2nd order - weighted LS3rd order - unweighted LS3rd order - weighted LS4th order - unweighted LS4th order - weighted LS0.781.722.881.874.402.91(a) Triangular mesheshL 2  (u e - u)0.02 0.04 0.06 0.08 0.110-410-310-210-11001011021.832.943.443.651.921.26(b) Quadrilateral meshesFigure 3.16: Solution reconstruction error using principal coordinate systemTo shed more light on this issue, the reconstructed solution and normal derivativeson a fixed position across the arc, θ = −1◦, have been plotted against the exact valuesin Figure 3.17 for the second triangular mesh. These are computed at differentdistances from the wall for a weighted LS reconstruction (Figure 3.16). At smalldistances from the wall where high aspect ratio triangles exist, the second-orderapproximation of solution and corresponding normal derivatives are highly off theexact solution in some locations. This can be explained by the lack of sufficientlyclose control volumes to the cell center under consideration as Mavriplis showed in his633.4. Numerical Testswork [77] for second-order reconstruction on triangular meshes. For fourth-order, thestencil becomes larger and includes geometrically close control volumes that improvethe estimate of normal derivatives; however, the reconstructed values are notably offthe exact curve for the weighted LS reconstruction.u ( r,  = -1 deg )Distance from Wall0 0.2 0.4 0.6 0.8 110-610-510-410-310-210-1Exact2nd-order4th-order(a) Solutionu/ nDistance from Wall0 5000 10000 1500010-610-510-410-310-210-1Exact2nd-order4th-order(b) Normal derivativeFigure 3.17: Weighted LS reconstructed values in the principal coordinate systemagainst exact values at different distances from the wall (θ = −1◦)On the other hand, the asymptotic order of reconstruction error using the newcurvilinear coordinate system is consistent with the order of reconstruction and theerror values are smaller by several orders of magnitude compared to the principalcoordinates. Figure 3.18 shows the reconstruction error for second- to fourth-orderreconstruction in the tangential-normal coordinate system for triangular and quadri-lateral cells. Also, the first derivative in the normal direction is one order less accuratethan the solution, as expected. For both types of cells, the unweighted LS producesmore accurate reconstruction coefficients, particularly for quadrilateral cells wherethe difference in the error of reconstructed solution is more than one order of mag-nitude in some cases. For normal derivatives, the differences between weighted andunweighted LS become smaller by increasing the order of reconstruction althoughthe unweighted LS outperforms in all of the cases. In general, the unweighted LSreconstruction in the curvilinear coordinate system solves the accuracy issues relatedto high aspect ratio meshes in the presence of curvature.Again, the solution and normal derivatives reconstructed by unweighted LS onthe curvilinear coordinate system at θ = −1◦ are plotted against the exact values643.4. Numerical TestshL 20.02 0.04 0.06 0.08 0.110-410-310-210-1100101u (unweighted LS)u/ n (unweighted LS)u (weighted LS)u/ n (weighted LS)1.920.97(a) Second-order/triangleshL 20.02 0.04 0.06 0.08 0.110-410-310-210-1100101u (unweighted LS)u/ n (unweighted LS)u (weighted LS)u/ n (weighted LS)1.921.511.600.87(b) Second-order/quadshL 20.02 0.04 0.06 0.08 0.110-510-410-310-210-11003.271.882.57(c) Third-order/triangleshL 20.02 0.04 0.06 0.08 0.110-510-410-310-210-11002.891.852.61(d) Third-order/quads653.4. Numerical TestshL 20.02 0.04 0.06 0.08 0.110-610-510-410-310-210-11003.842.91(e) Fourth-order/triangleshL 20.02 0.04 0.06 0.08 0.110-610-510-410-310-210-11002.663.863.44(f) Fourth-order/quadsFigure 3.18: Reconstruction error (relative for normal derivatives) using curvilineart− n coordinate systemon the second triangular mesh (Figure 3.19). In contrast to the results obtained byreconstruction on principal coordinates, both the solution and normal derivatives arereasonably close to the exact values and become more accurate by increasing theorder of accuracy.Now, we turn our attention to the conditioning of the least-squares reconstructionin each coordinate system. Table 3.3 gives the maximum condition number of LS sys-tem for reconstruction along principal or Cartesian coordinates for triangular meshes.Note that “no” in principal coordinates column implies the traditional reconstructionin the Cartesian coordinate system. Clearly, the use of principal coordinates and col-umn scaling highly enhance the conditioning as they reduce κ∞ (A) particularly forhigher orders by several orders of magnitude. Even though the combination of thesetwo does not make the condition number independent of mesh properties in this case,it is highly recommended even for the auxiliary LS problem having the same features.Having this in mind, one can compare the maximum condition numbers for thescaled system in principal coordinates for different ARmax and h. Second-order recon-struction is always well-conditioned regardless of weight function and mesh properties.For higher-order, the unweighted LS gives the condition number being scaled linearly663.4. Numerical Testsu ( r,  = -1 deg )Distance from Wall0 0.2 0.4 0.6 0.8 110-610-510-410-310-210-1Exact2nd-order4th-order(a) Solutionu/ nDistance from Wall0 5000 10000 1500010-610-510-410-310-210-1Exact2nd-order4th-order(b) Normal derivativeFigure 3.19: Reconstructed values in the curvilinear coordinate system against exactvalues at different distances from the wall (θ = −1◦)with mesh spacing at a constant aspect ratio whereas the changes are smaller for theweighted system. Likewise, the changes with ARmax at a constant h are more se-vere for the unweighted LS. In general, weighting gives smaller condition numbers forhigher-order on curved meshes; however, it does not guarantee a perfect conditioningparticularly for fourth-order reconstruction on high aspect ratio meshes.On the other hand, the unweighted LS along the curvilinear coordinate system isalways well-conditioned as shown by Table 3.4. Based on the conditioning results re-ported so far, we know that column scaling is essential for higher-order reconstructionon anisotropic meshes. Therefore, we just consider the maximum condition numberfor the scaled LS system. The conditioning results for the curvilinear coordinatesare completely similar to those obtained for principal axes on straight meshes. Thecondition numbers become independent of mesh properties for all orders of accuracyfor the unweighted system. On the other hand, they are scaled with ARmax butremain independent of h for the weighted case. As a result, the unweighted LS onthe curvilinear coordinates outperforms other schemes for anisotropic curved meshesin terms of both accuracy and conditioning. Note that the same type of conclusioncan be drawn for the conditioning of solution reconstruction on quadrilateral meshesalthough we skip their results for brevity.673.4. Numerical TestsARmax h Col. Scale Princ. Coord.κ∞ (A)k = 2 k = 3 k = 4unweighted LS2.06× 103 0.0436 yes yes 2.43 254.94 3.46× 1032.06× 103 0.0436 yes no 54.53 5.46× 103 2.81× 1062.06× 103 0.0218 yes yes 2.43 135.15 1.88× 1032.06× 103 0.0218 no yes 193.84 6.14× 106 4.57× 10116.32× 103 0.0436 yes yes 4.98 842.41 1.05× 1046.32× 103 0.0436 yes no 79.49 5.61× 103 8.36× 1066.32× 103 0.0436 no no 158.90 8.02× 105 6.84× 10101.75× 104 0.0218 yes yes 2.44 1.15× 103 1.56× 1041.75× 104 0.0218 yes no 117.13 2.55× 104 1.56× 108weighted LS2.06× 103 0.0436 yes yes 2.21 30.65 561.582.06× 103 0.0218 yes yes 2.22 42.88 501.296.32× 103 0.0872 yes yes 15.50 53.99 6.64× 1036.32× 103 0.0436 yes yes 17.35 162.56 3.81× 1036.32× 103 0.0436 yes no 55.85 3.53× 103 2.25× 1066.32× 103 0.0436 no yes 119.23 2.42× 106 4.35× 10106.32× 103 0.0218 yes yes 13.70 399.98 2.74× 1031.75× 104 0.0436 yes yes 2.22 48.93 4.77× 1031.75× 104 0.0436 yes no 53.69 3.00× 103 6.52× 106Table 3.3: Maximum condition number for reconstruction along principal(Cartesian)axes on curved meshes3.4.3 General MeshesAs a more general case for anisotropic meshes over curved surfaces, we use the struc-tured grids generated for turbulent RANS simulations by NASA Langley ResearchCenter [85]. For experiments on triangular meshes, we also randomly triangulatethe quadrilaterals (Figure 3.20). We curve the interior faces of the mesh to fourth-order using a cubic spline representation of the airfoil surface. We only consider theunweighted LS because of the conditioning issues discussed earlier.The meshes typically used for viscous flow simulations are comprised of differentregions: anisotropic cells with curvature very close to solid walls (R1), anisotropic683.4. Numerical TestsARmax h column scalingκ∞ (A)k = 2 k = 3 k = 4unweighted LS2.06× 103 0.0436 yes 1.74 13.99 101.812.06× 103 0.0218 yes 1.74 13.86 97.486.32× 103 0.0872 yes 2.06 13.58 91.496.32× 103 0.0436 yes 2.23 15.61 97.176.32× 103 0.0218 yes 2.31 17.47 105.361.75× 104 0.0218 yes 1.73 13.84 97.43weighted LS2.06× 103 0.0436 yes 2.21 4.85× 103 4.17× 1042.06× 103 0.0218 yes 2.22 5.29× 103 4.50× 1046.32× 103 0.0436 yes 2.20 2.86× 104 1.51× 1051.75× 104 0.0436 yes 2.22 4.07× 104 3.53× 105Table 3.4: Maximum condition number for reconstruction along local curvilinear axeson curved meshes-400 -200 0 200 400-400-20002004000.5 1 1.5 2 2.5-1-0.500.510 0.5 1-0.500.5Figure 3.20: Unstructured triangular mesh over NACA 0012693.4. Numerical Testsstraight cells often seen in wake regions (R2) and isotropic cells sufficiently far fromthe walls (R3). Since the calculation of mapping coefficients and moments of areaare computationally more expensive for the curvilinear coordinate system, we needto isolate the regions that require this treatment. We already know that higher-orderleast-squares reconstruction even in Cartesian coordinate system provides accurateresults along with good conditioning for isotropic cells. In addition, our previousresults for anisotropic straight meshes suggest a principal coordinate system andcolumn scaling for accurate and well-conditioned LS reconstruction. As a results, onlycells with high aspect ratio and significant curvature demand curvilinear coordinates.Such cells can be recognized as those that have a large aspect ratio (AR > 10 here)and are fairly well aligned with a wall. The latter can be understood by the anglebetween the line that connects the reference location to the closest point on thewall and the principal axis around which the second moment of area is minimum(γ < 10◦). Note that the cells in R2 only satisfy the first criterion. Figure 3.21 showsthe separation of regions for NACA 0012 meshes using the mentioned procedure.0.2 0.3 0.4 0.5 0.600.10.21 1.04-0.0200.020.02 0.04 0.06R1 R2 R3Figure 3.21: Separate regions for NACA 0012 meshesFor reconstruction tests, we just consider R1 colored by blue because that is theregion where traditional LS suffers from accuracy issues. Note that the cells in R2are identical to straight meshes we considered in Section 3.4.1. To prove that theimprovements in LS reconstruction on curvilinear coordinate system is not limitedonly to simple geometries, we perform our accuracy tests on the cells of R1 forthe unstructured meshes around the NACA 0012. For this purpose, an anisotropicsolution is manufactured for this region and some layers of neighbors in other regions.We re-employ the boundary layer profile of Equation 3.27 for the cells with a referencelocation in −0.1 < y < 0.1 and x > −0.15 and initialize other control volume averageswith zero. In Equation 3.31, y¯ is replaced with distance from the wall for the cells703.4. Numerical Testsnear the wall with a reference location in 0 < x < 1 and x¯ is replaced with x+0.2 toprevent singularities near the leading edge:uei (x, y) ={2.5 ln(1 + 0.4 y+)+ 7.8(1− exp(−y+11)− y+11exp(−0.33 y+))}u∗u∗ =√√√√ 0.0135(ReL. (x+ 0.2))1/7, y+ = d.ReL.u∗ (3.32)Also, ReL = 6× 106 as this is the Reynolds number for which numerical results havebeen reported in Ref. [85] and thus matches the anisotropy in the mesh. Figure3.22 illustrates the manufactured solution at different locations near the wall on thecoarsest mesh.0.13 0.14 0.15 0.16 0.17 0.2 0.4 0.6 0.8 1(a)0 0.0025 0.0050.00250.0050.00750.010.0125(b)0.98 1 1.02-0.004-0.00200.0020.004(c)Figure 3.22: Manufactured solution for reconstruction accuracy test on NACA 0012The accuracy test is performed on a sequence of quadrilateral and triangularmeshes. For this purpose, the average of the absolute value of the reconstructionerror, E¯h , is calculated over each control volume as:∣∣∣E¯h∣∣∣i=1ACViˆ ˆCVi∣∣∣uRi (x, y)− uei (x, y)∣∣∣ dA (3.33)713.4. Numerical Tests(Num. CV)1/2L 2 (|Eh|)100 200 300 400 50010-610-510-410-310-22nd - Principal3rd - Principal4th - Principal2nd - Curvilinear3rd- Curvilinear4th - Curvilinear2:13:14:1(a) Quadrilateral cells(Num. CV)1/2L 2 (|Eh|)200 400 600 80010-710-610-510-410-310-22nd - Principal3rd - Principal4th - Principal2nd - Curvilinear3rd- Curvilinear4th - Curvilinear3:14:12:1(b) Triangular cellsFigure 3.23: Error norms of reconstructed solution on anisotropic cells of NACA 0012where ACVi is the surface are for control volume i, uRi (x, y) is the reconstructed solu-tion profile in the control volume and uei (x, y) is the exact solution profile obtainedby Equation 3.32 at each point of the control volume. The integral of Equation 3.33is computed by employing a sixth-order quadrature rule and evaluating the differ-ence between exact and reconstructed values at the quadrature points of each controlvolume.One can compare the error norms obtained for unweighted LS reconstructionusing principal and curvilinear coordinate systems. Using the local curvilinear co-ordinate system leads to a significant improvement in the accuracy of reconstructedsolutions as seen in Figure 3.23. The error values become considerably smaller forboth the quad and triangular anisotropic cells and the nominal order of reconstruc-tion is obtained by mesh refinement. On the other hand, higher than second-orderreconstruction yields slightly better than second-order reconstructed solutions usingprincipal coordinates. It is clear that the approximation of the flux vector and sourceterm with low-order reconstructed solutions cannot lead to a higher-order accuratesolution at the end. Although this improvement was shown only for the reconstruc-tion of a particular manufactured solution, its advantage in evaluating more accuratesolutions and outputs will be shown later in the next chapter.72Chapter 4RANS Simulation of TurbulentFlowsWith the success of Discontinuous Galerkin (DG) methods in obtaining more accu-rate solutions in computational aerodynamics, more realistic flow conditions such asturbulent flows governed by the Reynolds Averaged Navier-Stokes (RANS) equationshave been tackled in recent years. The earliest success in this area was the implicithigher-order RANS solver of Bassi et al. [86] based on the k − ω turbulence modelin which the solutions of turbulent flows on a flat plate and turbine blades wereinvestigated. Later, Nguyen et al. [87] used the one equation turbulence model ofSpalart-Allmaras (SA) in their higher-order DG solver for the solution of turbulentflat plate and 2D airfoil test cases. In recent years, a great deal of effort has beendevoted to increase the efficiency and robustness of higher-order DG RANS solvers[88, 89, 90] and also to expand the capability of solving turbulent flow problems withmore complex geometries [91]. In addition to these, some other discretization meth-ods such as streamline upwind/Petrov Galerkin (SUPG) [92] or residual distribution(RD) [93] schemes have been used to deliver a higher-order solution of RANS equa-tions. These methods are less expensive as the solution is assumed to be continuousacross the interfaces and thus the number of degrees of freedom grows less rapidlywhen increasing the order of polynomial.On the other hand, the higher-order finite volume methods, which have showntheir ability for efficient computations of inviscid and viscous laminar flow withpromising accuracy on irregular meshes, have not been used for the solution of tur-bulent compressible flows. In this chapter, we describe the remaining challengesencountered in the extension of higher-order unstructured finite volume methods toRANS simulations and also the treatments required to tackle each of these issues.This chapter presents the development of an implicit higher-order unstructured finitevolume solver for turbulent aerodynamic flows.734.1. Governing EquationsOne of the challenges in the application of higher-order methods to RANS sim-ulations with the Spalart-Allmaras model is the abrupt change in the slope of theturbulence working variable at the edge of boundary layer. This behavior resultsin negative values of the turbulence working variable and ultimately causes solverfailure. Several modifications [88, 94] have been proposed for this model to make itbenign for higher-order discretizations across the slope discontinuity. In the presentwork, we use the SA-neg model developed by Allmaras et al. [95]. We also fully cou-ple the turbulence model equation with the mean flow equations and the RANS-SAsystem is considered as a complete system of equations.Another important task in a higher-order flow solver is robust and efficient con-vergence to the steady-state solution. For a turbulent flow simulation, anisotropiccells induce significant stiffness into the discrete equations and hamper solution con-vergence. Consequently, special attention should be paid to the solution strategy forthese problems. In this work, we adopt a solution strategy originally developed for ahigher-order DG RANS solver [90], with minor modifications.An outline of this chapter is as follows. In Section 4.1, the governing equations in-cluding the negative variant of the SA model are described. We present the numericalflux functions used for the fully-coupled RANS-SA system in Section 2.3. Section 4.3briefly reviews the solution strategy used for the steady-state solution of turbulentflows. Four numerical examples, including a turbulent flat plate, subsonic and tran-sonic airfoils, and finally a multi-element configuration are presented in Section 4.4to highlight the ability of higher-order methods to obtain a more accurate solutionon coarser meshes and also efficient convergence to the steady-state solution.4.1 Governing EquationsIn this work, the RANS equations are coupled to the one equation turbulence modelof Spalart and Allmaras (SA model) [19]. The original SA model admits only non-negative values of the working variable. Recent studies in the context of higher-orderDG methods have found difficulties in the robust application of the Spalart-Allmarasmodel due to non-smooth solution behavior of this turbulence model at the edge of theflow boundary layer [89]. This behavior results in negative values of the SA workingvariable and ultimately causes solver failure. These negative values are generatedby Gibbs phenomena that stem from employing high-order approximations across a744.1. Governing Equationsdiscontinuity in slope. Alternatively, we use the negative variant of the SA modelproposed by Allmaras et al. [95]:∂ρν˜∂t+∇ ·(ρν˜ ~V)−∇ ·((µ+ µ′fnρν˜)∇ν˜σ)= ρ (P −D) + µ′σcb2ρ∇ν˜ · ∇ν˜−1σ(ν + µ′ν˜)∇ρ · ∇ν˜ (4.1)where µ′ = 1000 is a scaling factor used to improve the convergence rate of theimplicit Newton-Krylov solver by making the turbulence working variable comparablein magnitude to other physical quantities [96]. In Equation 4.1, P and D denote theproduction and destruction of eddy viscosity which depend on the sign of the eddyviscosity and are defined as:P =cb1 (1− ft2) S˜ν˜ ν˜ ≥ 0cb1 (1− ct3)Sν˜ ν˜ < 0 , D =µ′(cw1fw − cb1κ2 ft2) (ν˜d)2ν˜ ≥ 0−µ′cw1(ν˜d)2ν˜ < 0(4.2)In Equation 4.2, S˜ is the modified vorticity, which must always remain positive:S˜ =S + S¯ S¯ ≥ −cv2SS +S(c2v2S+cv3S¯)(cv3−2cv2)S−S¯S¯ ≤ −cv2S(4.3)S =∣∣∣∣∣∂u∂y − ∂v∂x∣∣∣∣∣ (4.4)S¯ =µ′ν˜fv2κ2d2and d is the distance to the closest wall. The functions fv1 and fv2 are :fv1 =χ3χ3 + c3v1, fv2 = 1− χ1 + χfv1 , χ =µ′ν˜ν(4.5)and functions fn and ft2 are defined as:fn =cn1 + χ3cn1 − χ3 , ft2 = ct3 exp(−ct4χ2)(4.6)754.1. Governing EquationsThe destruction term coefficients are:r = min(µ′ν˜S˜κ2d2, 10)g = r + cw2(r6 − r)(4.7)fw = g[1 + c6w3g6 + c6w3]The constant values used in this model are as follows: σ = 0.66, κ = 0.41, cv1 = 7.1,cv2 = 0.7, cv3 = 0.9, cb1 = 0.1355, cb2 = 0.622, cw1 =cb1κ2+ 1+cb2σ, cw2 = 0.3, cw3 = 2.0,ct3 = 1.2, ct4 = 0.5, cn1 = 16.The conservative form of the compressible Reynolds Averaged Navier-Stokes (RANS)equations combined with the SA-neg turbulence model equation can be re-arrangedas:∂U∂t+∇ ·(~Fc (U)− ~Fv (U,∇U))= S (U,∇U) (4.8)The solution, flux and source term vectors for the coupled system of RANS-SA intwo-dimensions are given as:U =ρρuρvEtρν˜, F xc =ρuρu2 + Pρuvu (Et + P )ρuν˜, F yc =ρvρuvρv2 + Pv (Et + P )ρvν˜F xv =0τxxτxyuτxx + vτxy + cp(µPr+ µTPrT)∂T∂x1σ(µ+ µ′fnρν˜)∂ν˜∂x, F yv =0τyxτyyuτyx + vτyy + cp(µPr+ µTPrT)∂T∂y1σ(µ+ µ′fnρν˜)∂ν˜∂yS =0000ρ (P −D) + µ′σcb2ρ∇ν˜ · ∇ν˜ − 1σ (ν + µ′ν˜)∇ρ · ∇ν˜(4.9)764.2. Flux Functionswhere PrTis the turbulent Prandtl number, µT is the turbulent eddy viscosity andthe other variables are the same as those described for the Navier-Stokes equationsin Chapter 2. We assume that the working fluid is air with PrT= 0.9; the turbulenteddy viscosity is given as:µT =µ′ρν˜fv1 ν˜ ≥ 00 ν˜ < 0(4.10)Also, the total viscous stress tensor, τij , including the Boussinesq approximatedReynolds stresses is found as:τij = 2 (µ+ µT ) γ˙ij (4.11)4.2 Flux FunctionsMany turbulent flow solvers discretize the turbulence model equation in a loosely-coupled manner, which treats the convection term as though it were a scalar transportequation evolving with a prescribed velocity field. However, this treatment neglectsthe fact that the velocity field, which convects the turbulence model quantities, isheavily influenced by the turbulence model solution. Recent work by Burgess andMavriplis [45] demonstrated the importance of full coupling of the RANS equationswith the one equation SA model. In this section, we describe the numerical fluxfunctions used for the discretization of convective and viscous fluxes.4.2.1 Convective FluxesThe RANS equations closed with the SA turbulence model results in a total of fiveequations in two spatial dimensions whose flux function must be re-derived. To usea flux-difference splitting method, we need the eigenvalues of the Jacobian of thenormal convective flux for the system of RANS equations:F nc =ρunˆx + ρvnˆy(ρu2 + P ) nˆx + ρuvnˆyρuvnˆx + (ρv2 + P ) nˆyu (Et + P ) nˆx + v (Et + P ) nˆyρuν˜nˆx + ρvν˜nˆy(4.12)774.2. Flux FunctionsIt is possible to show that the eigenvalues of ∂Fnc∂Uare:λ1 = ~V · nˆ− a λ2 = λ3 = λ4 = ~V · nˆ λ5 = ~V · nˆ + a (4.13)in which a is the sound speed. We use the numerical flux function of Roe and Pike[97] that was derived by Burgess and Mavriplis [45] for the RANS-SA system:Fc (UL, UR) =12(Fc (UL) + Fc (UR)−D) (4.14)where D is the dissipative component of the numerical flux, which is given as:D =∣∣∣λ˜2∣∣∣ (ρL − ρR) + δ1∣∣∣λ˜2∣∣∣ (ρuL − ρuR) + δ1u˜+ δ2nˆx∣∣∣λ˜2∣∣∣ (ρvL − ρvR) + δ1v˜ + δ2nˆy∣∣∣λ˜2∣∣∣ (Et,L −Et,R) + δ1H˜ + δ2 (u˜nˆx + v˜nˆy)∣∣∣λ˜2∣∣∣ (ρν˜L − ρν˜R) + δ1 ˜˜ν(4.15)δ1 =− ∣∣∣λ˜2∣∣∣+∣∣∣λ˜1∣∣∣+ ∣∣∣λ˜5∣∣∣2 ∆Pa˜2+∣∣∣λ˜5∣∣∣− ∣∣∣λ˜1∣∣∣2ρ˜a˜(nˆx∆u+ nˆy∆v)δ2 =− ∣∣∣λ˜2∣∣∣+∣∣∣λ˜1∣∣∣+ ∣∣∣λ˜5∣∣∣2 ρ˜ (nˆx∆u+ nˆy∆v) +∣∣∣λ˜5∣∣∣− ∣∣∣λ˜1∣∣∣2∆Pa˜Note that in Equation 4.15, ∆ () = ()L − ()R and X˜ refers to quantity X evaluatedat the Roe state [59]. In particular, ˜˜ν is the turbulence model working variable atRoe’s state. This flux function is computationally less expensive compared to theimplementation of Roe’s flux function available in our solver as it does not requireany matrix-matrix or matrix-vector multiplication.4.2.2 Viscous FluxesAs described in Chapter 2, the gradient at the face quadrature points can be calcu-lated by the average of the two reconstructed gradients and that value is employedto calculate the diffusive fluxes. However, this strategy fails for highly anisotropicmeshes. To shed light on this issue, we consider the Poisson equation, ∇2φ = S, asa model of viscous discretization along with the method of manufactured solutions784.2. Flux Functionson anisotropic meshes. For this purpose, we create irregular anisotropic meshes ona rectangular domain (x, y) ∈ [0, 1] × [0, 0.05] using the procedure that has beendescribed in detail by Diskin et al [98]. Figure 4.1 depicts the anisotropic mesh testcase with 20 × 20 cells. For generating the test case, the first step is stretching aregular rectangular grid with (N + 1)× (N + 1) nodes toward the bottom line y = 0using a stretching factor s and the maximum aspect ratio ARmax. The y-coordinatesof the horizontal grid lines in the top half of the domain are defined as:yk = yk−1 +skARmax · N (4.16)in whichy0 = 0 , k = 1, 2, ..., N (4.17)Then, irregularities are introduced by random shifts on interior nodes in vertical andhorizontal directions and finally each perturbed quadrilateral is randomly triangu-lated with one of the two choices for its diagonal. In Figure 4.1 where N = 20, thestretching factor and the maximum aspect ratio are set as s = 1.14 and ARmax = 100to yield triangles whose aspect ratios vary approximately between 10 and 100. Thismesh consists of anisotropic triangles which are long in the x-direction and thin inthe y-direction.In this study, these meshes with N = 10, 20, 40, 80 are used where the stretchingfactors are s = 1.28, 1.14, 1.065, 1.0335, respectively. The s values are chosen so asto keep the the length scale reduction in the y-direction the same as the x−direction.Since the anisotropic meshes are used for anisotropic solutions, the function beingused for test purposes is manufactured to comply with the cells’ aspect ratio. Theanisotropic mesh of Figure 4.1 is stretched toward the horizontal line y = 0 usingEquation 4.16 which is strongly dependent on the y-coordinates. The variations inthe x-direction are introduced by an exponential function to yield isotropic behaviorin this direction and homogeneous boundary conditions:φ (x, y) =exp (x (x− 1) y (y − 0.05))− 1sARmax+N(1− 1s)y(4.18)In this part, Equation 4.18 is used as the test function with geometric valuescorresponding to the coarsest mesh (N = 10 and s = 1.28). Figure 4.1 depicts the794.2. Flux Functionsmanufactured solution on one of the constructed anisotropic mesh along with its y-derivative which is too large near the bottom line with highly stretched cells and isreduced smoothly in the vertical direction compliant with cells’ aspect ratio.0 0.2 0.4 0.6 0.8 0.0003 0.0007 0.0010 0.0013 0.0017 0.0020 0.0023 0.0027 0.0030(a) Solution0 0.2 0.4 0.6 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 / y(b) y-derivativeFigure 4.1: Anisotropic solution and derivative on a stretched grid with 20× 20 cellsFigure 4.2 shows the eigenvalues of the global Jacobian matrix for a second-orderdiscretization of Poisson’s equation on this mesh. It is worth-mentioning that theglobal Jacobian matrix of a linear problem is independent of the solution. Notethat pure averaging of the two reconstructed gradients for viscous flux discretizationleads to the appearance of eigenvalues in the right-half plane that has potential tomake the solution procedure unstable. However, a jump term which comes from the804.3. Solution Methoddiscontinuous solution at the Gauss points can be added to the face gradient [99]:(∇φ)F =12((∇φ)R + (∇φ)L) +α|~rij .nˆ| (φL − φR) nˆ (4.19)where α is an arbitrary value (we use α = 1 ), ~rij is the vector that connects thereference points of the two control volumes that share the face and nˆ is the outwardunit normal vector at the Gauss point. The jump term has the role of stabilizing bydamping high-frequency errors. This can also be seen in Figure 4.2 as adding thisterm pushes all the eigenvalues to the left-half plane [100]. We generalize this ideato the complicated case of the non-linear RANS equations on anisotropic meshes anduse Equation 4.19 to compute the gradients involved in the viscous fluxes of Equation4.9.Re( )Img()-440000 0 440000-400000-2000000200000400000AveragingAveraging + JumpFigure 4.2: Comparison of eigenvalue spectra for the discretization of viscous fluxeson anisotropic meshes4.3 Solution MethodIn the context of RANS simulations with finite volume solvers, it is common to useexplicit Runge-Kutta time-stepping even for the steady-state solution of the flow field.However, these methods slow down the convergence significantly unless combinedwith multigrid solution techniques [101]. Alternatively, we use the implicit pseudotime-stepping method described in Section 2.5 to accelerate the convergence of ourflow solver: (I∆t/A+∂R∂U¯)δU = −R(U¯n)(4.20)814.3. Solution MethodFor the sake of brevity, we re-write this equation as:L δUi = −R(U¯n)(4.21)The convergence to the steady-state solution can be posed as an optimizationproblem where we search for the minimum of an objective function. We define ourobjective function as:f(U˜)=12∣∣∣Rt (U˜)∣∣∣2 = 12RTt(U˜)Rt(U˜)(4.22)where the unsteady residual is defined as:Rt(U˜)=I∆t/A(U˜ − U¯n)+R(U˜)(4.23)In other words, we seek a trial solution, U˜ , that minimizes the norm of the unsteadyresidual vector. Multiplying both sides of Equation 4.21 by the transpose of itsleft-hand side givesδUT LT L δU = −δUT LT R(U¯ni)= −δUT ∂f∂U˜∣∣∣∣∣U¯n> 0 (4.24)where the inequality comes from the fact that the left-hand side is the dot productof a non-zero vector with itself. Considering thatδUT∂f∂U˜∣∣∣∣∣U¯n< 0, (4.25)δU is a descent direction for the scalar function f(U˜). As a result, we can incorporatea line search algorithm with the implicit solver to enhance the robustness. For thispurpose, an under-relaxation factor, ωn , is used at each iteration to update thesolution:U¯n+1 = U¯n + ωnδU (4.26)The under-relaxation factor can be found in a way to guarantee sufficient decrease inthe norm of the unsteady residual vector such that∥∥∥Rt (U˜)∥∥∥2<∥∥∥R (U¯n)∥∥∥2. However,824.3. Solution MethodCeze and Fidkowski showed that relaxing the sufficient decrease condition to∥∥∥Rt (U˜)∥∥∥2< κLS∥∥∥R (U¯n)∥∥∥2(4.27)with κLS > 1 is helpful in accelerating the convergence of a fully coupled RANS solver[90]. In our work, we use κLS = 1.2 as the relaxation parameter for the sufficientnon-increase condition. We also need to choose the under-relaxation factor so as toensure a physical solution at the next iteration. Note that the physicality of thesolution state is verified by checking the positivity of density and pressure at allquadrature points. Therefore, we start with the initial guess of ωn = 1 and checkif the new solution state is physical and satisfies the relaxed decrease condition. Ifnot, we halve the under-relaxation factor every time until the both conditions aresatisfied.For time accurate solutions, the time step in Equation 4.20 must be the same forall control volumes. However, local time-stepping with a global CFL number can beused when we seek the steady-state solution to accelerate the convergence process.The control volume wise time steps at each iteration are obtained as:∆tni =CFLn · hiλmax, i(4.28)where λmax, i is the eigenvalue of the convective flux Jacobian with the largest mag-nitude obtained from Equation 4.13, and hi is a characteristic size for the controlvolume. In the first stages of the computations, the initial solution does not satisfythe boundary conditions and thus strong transients take place in the solution. It isimportant to keep the global CFL number small at these stages to make the solutionfollow a physical path. As we proceed, the CFL number must be increased to makethe convergence to the steady-state solution faster. Therefore, an evolution strategyis required for the CFL from its initial value to a large value such that Equation 4.20becomes Newton’s method and the solution approaches the steady state. We use theexponential progression with under-relaxation strategy proposed by Ceze and Fid-kowski [96] in which the size of CFL growth is tuned by the under-relaxation factor.This strategy can be summarized as:834.4. ResultsCFLn+1 =β · CFLn ωn = 1CFLn ωmin < ωn < 1κ · CFLn ωn < ωmin(4.29)In our work, we set the parameters to β = 1.5, κ = 0.1 and ωmin = 0.01. Notethat when we search for the appropriate under-relaxation factor which satisfies boththe relaxed sufficient decrease and physicality conditions, we terminate the searchprocess for ωn < ωmin and cut the global CFL number according to Equation ResultsIn this section, we present our results for four fully turbulent cases: subsonic flow overa flat plate, subsonic flow over the NACA 0012 airfoil, transonic flow over the RAE2822 airfoil and flow around a high-lift multi-element configuration [102]. Our resultsinclude verification of the turbulence model discretization, accuracy advantages ob-tained by higher-order discretizations and the convergence behavior of each case. Forall these test cases, the no-slip and zero turbulent viscosity boundary conditions onthe walls are applied by adding an extra constraint to the least-squares system forsolution reconstruction for each boundary Gauss point. In addition, the adiabaticwall condition is applied weakly by zeroing the heat flux components in Equation 4.9for wall boundary faces. For inflow, we specify the values of total pressure, velocitycomponents and turbulent viscosity while the static pressure is the only pre-specifiedcondition at the outflow; other quantities are obtained from solution reconstruction.No point vortex correction is done for the far-field boundaries. Also, the dynamicviscosity is obtained by Sutherland’s law.4.4.1 High Reynolds Number Flow Over a Flat PlateFor this test case, we consider the 2D flat plate verification case from the NASATurbulence Modeling Resource (TMR) website [85]. Figure 4.3 shows the problemgeometry and boundary conditions for this problem with Re = 5×106 andMa∞ = 0.2.The range of the computational domain in the x-direction is [−0.33, 2] with theleading edge of the flat plate at x = 0. The size of the domain in the y-direction is 1.At y = 0, the symmetry boundary condition is applied for −0.33 ≤ x ≤ 0 weakly by844.4. Resultsforcing the normal derivatives to be zero in the flux vectors and the no slip adiabaticboundary condition is imposed for 0 ≤ x ≤ 2. Note that this causes a singularityin the solution at the leading edge. The free-stream value of the turbulent workingvariable is ν˜/ν∞ = 3.0.XY-0.5 0 0.5 1 1.5 200.51farfield Reimann BCoutflowsymmetryinflowadiabatic solid wallFigure 4.3: Turbulent flat plate test case mesh and geometryFor this test case, we generated a sequence of nested grids with quadrilateral cells.The coarsest mesh, which has 2040 cells has 61 nodes in the x- and 35 nodes in they-direction. The grid is highly stretched near the wall such that the height of thefirst cell is 9.2 × 10−6. Considering that our higher-order discretization employs alarger stencil for the Taylor series expansion, having anisotropic cells near singularitiescauses instability in the solution process. Therefore, we cluster the cells near theleading edge to produce nearly isotropic cells (AR = O(1)) in that region. Finermeshes are obtained by uniform refinement of the coarsest mesh.Figure 4.4 shows the skin friction along the wall and turbulent viscosity at x = 0.97near the wall for second- to fourth-order discretizations on a mesh with 32, 640 controlvolumes. These profiles are converging to the grid-converged values obtained bythe FUN3D solver from NASA Langley [85], which uses the SA-Neg model, on amuch finer mesh. As expected, higher-order discretizations lead to more accurateestimates of the flow field on the same mesh. Closer views of the plots show thatthird- and fourth-order results are closer to the expected values. In addition, thefourth-order results for 32, 640 control volumes are comparable with second-orderresults on the next finer mesh (n.DoF = 130, 560). Note that in the context offinite volume methods, the number of degrees of freedom (n.DoF) per equation isdefined as the number of control volumes in the mesh independent of the order of854.4. Resultsaccuracy. These figures verify the correctness of our RANS flow solver and also showthe accuracy advantages obtained by employing higher-order discretizations.We also examine the convergence of friction drag coefficient on the sequence ofsystematically refined meshes for different orders. Table 4.1 shows the values of CD ondifferent meshes and error in the drag coefficient, defined as CD,Err = |CD − CD,Ref |.Note that the reference value is obtained by fourth-order discretization on the finestmesh with 130, 560 control volumes. As expected, the error decreases by increasingthe order of accuracy on the same mesh. Also, the drag values converge faster usinga higher-order discretization so that the error associated with fourth-order on thesecond finest mesh is smaller than the error of the second-order scheme on the finestmesh.n.DoFCD CD,Err2nd 3rd 4th 2nd 3rd 4th2040 0.003617 0.002163 0.002329 7.54× 10−4 7.00× 10−4 5.34× 10−48160 0.002979 0.002803 0.002849 1.16× 10−4 6.0× 10−5 1.4× 10−532640 0.002873 0.002862 0.002862 1.0× 10−5 2× 10−6 1× 10−6130560 0.002865 0.002864 0.002863 2× 10−6 1× 10−6 0Table 4.1: Friction drag coefficient of turbulent flat plate for different orders andmesh sizesTo investigate the iterative convergence properties of our flow solver, relevantparameters have been tabulated in Table 4.2 for each order and three mesh sizes. Thistable gives the number of linear and non-linear iterations, total CPU time on a singlecore of an i7-4790 (3.60 GHz) CPU, and total number of work units until convergence(CPU time divided by the time of a single residual evaluation). In addition, thetime for residual evaluation, Jacobian matrix formation and linear system solutionaveraged over the non-linear iterations is listed. Note that the Jacobian matrix iscalculated to full order to enhance the convergence speed of the solver. The initialsolution is the uniform free stream condition everywhere and the linear system solveris GMRES preconditioned with ILU(2). The size of the Krylov subspace is 100 forall the cases and the relative residual of the linear solver is ηl = 10−3. Quotientminimum degree (QMD) re-ordering is applied at the time of preconditioning. Theconvergence process terminates when the L2-norm of the residual vector drops below10−8.An interesting observation is that the number of non-linear iterations is constant864.4. ResultsxC f0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20.0020.00250.0030.00350.0040.00450.0050.00550.006 2nd - n.DoF = 130,5602nd - n.DoF = 32,6403rd - n.DoF = 32,6404th - n.DoF = 32,640FUN3D - n.DoF = 208,896(a) Wall friction coefficient(b) Turbulent viscosity at x = 0.97Figure 4.4: Distribution of wall friction factor and turbulent viscosity for flat plate874.4. Resultsindependent of the order on the same mesh and the number of linear iterationsdecreases with increasing order. With mesh refinement, the number of non-lineariterations slightly increases although considerably more linear iterations are required.As expected, the total time of computation increases with mesh and order refinement.It should be noted that the time of residual evaluation and Jacobian matrix formationare fully scalable as they increase by a factor of 4 in each level of mesh refinement.However, the time spent to solve the linear system does not linearly increase withthe number of degrees of freedom. This is also expected since mesh refinement leadsto a larger linear system with larger condition number while the size of the Krylovsubspace is fixed.Figure 4.5 shows the convergence history for the flat plate test case on a meshwith 32, 640 cells. The drop in the L2-norm of the residual vector in terms of non-linear iteration count is very similar for different orders. However, the difference canbe found by considering the numbers of work unit. For all meshes, the number ofwork units is almost the same between second- and third-order and is about twice aslarge for fourth-order. This can be explained by the extra time required to form theJacobian matrix with more non-zeros and factorize it for preconditioning in the caseof fourth-order.Order n.DoFNo. of CPU Work Average Time (Sec.)iterations TimeUnitRes. Jac. LSLinear Non-linear (Sec.) Eval. Form. Solver2nd2, 040 291 28 6.08 1330 0.004 0.112 0.0998, 160 696 29 27.28 1690 0.016 0.411 0.51032, 640 2701 31 207.03 2840 0.073 1.644 4.9543rd2, 040 155 28 16.68 1770 0.009 0.274 0.3118, 160 318 29 74.64 2030 0.036 1.071 1.46632, 640 668 31 366.17 2430 0.150 4.327 7.3234th2, 040 98 28 31.70 1810 0.017 0.347 0.7728, 160 269 29 164.82 3630 0.045 1.379 4.25632, 640 576 31 865.84 4850 0.178 5.653 22.078Table 4.2: Convergence properties of turbulent RANS solver for high Reynolds num-ber flat plate884.4. ResultsIterationsL 2  Norm Res.5 10 15 20 25 30 3510-910-710-510-310-11011032nd-order3rd-order4th-order(a) Non-linear iterationsWork UnitL 2  Norm Res.102 103 10410-910-710-510-310-1101103 2nd-order3rd-order4th-order(b) Work unitsFigure 4.5: Convergence history of finite volume solver for high-Reynolds turbulentflat plate for n.DoF = 32,640894.4. Results4.4.2 Subsonic Flow Over a NACA 0012 AirfoilThe second test case considered is the subsonic flow over a NACA 0012 airfoil atangle of attack α = 10◦; the Mach number based on the free stream condition isMa∞ = 0.15 and the Reynolds number based on the airfoil chord is Rec = 6 × 106.This is also among the verification test cases of the TMR, which documents resultsobtained by established codes for different families of structured meshes. We use themeshes with the smallest spacing near the trailing edge (meshes of Family II [85])where the solution is singular. We also convert the structured quad meshes into hybridmeshes comprised of a few quad layers near the wall and triangles (obtained by split-ting quads) in other regions to have truly unstructured meshes. Computations areperformed on a sequence of four meshes with 6, 272, 25, 088, 100, 352, 401, 408 controlvolumes. Figure 4.6 shows an example of these meshes. No-slip adiabatic boundaryconditions are imposed on the airfoil surface and far-field conditions (inflow/outflowbased on the direction of velocity vector) are applied on the outer boundary located500c away from the airfoil. Again, the free-stream value of the turbulent workingvariable is ν˜/ν∞ = 3.0.Figure 4.6: Example of unstructured hybrid meshes used for subsonic flow over NACA0012 airfoil with 25, 088 degrees of freedomFigure 4.7 shows the contours of turbulence working variable obtained by second-and fourth-order discretizations on a hybrid mesh with 25, 088 degrees of freedom.The fourth-order method resolves the wake region more accurately and provides asmoother distribution of scaled turbulent viscosity.Figure 4.8 shows the distribution of pressure and friction coefficients on the wallfor second- and fourth-order methods using a mesh with 100, 352 cells. For both904.4. Resultsmethods, the values are consistent with the distributions obtained by FUN3D on asuper fine mesh (with about 1M cells); fourth-order gives more accurate values nearthe leading and trailing edges as expected. Note that these results are obtained byemploying a curvilinear coordinate system for highly anisotropic cells near the wallwith the procedure explained in Section 3.3.Table 4.3 gives the values of pressure drag, viscous drag and lift coefficients ondifferent meshes for all discretization orders using curvilinear and Cartesian coordi-nate systems. This table highlights the importance of accurate reconstruction on thefinal solution obtained by the solver. Comparing the values of coefficients with thereference values reported by TMR in the limit of mesh refinement reveals the factthat reconstruction in a curvilinear coordinate system leads to grid converged valuesconsiderably faster. As shown in Table 4.3a, coefficients of lift and drag convergeto the reference values with mesh and order refinement. Also, the convergence rateincreases for higher-order discretizations. Note that the second- and third-order re-sults on the finest mesh and fourth-order result on the second finest mesh are veryclose to the reference values. On the other hand, using a Cartesian coordinate systemgives poor output values, particularly for second- and fourth-order, even on very finemeshes (Table 4.3b). Although the values of drag and lift coefficients may convergeto the reference values on super fine meshes, this is in contrast with the motivation ofhigher-order methods, which is supposed to produce more accurate results on coarsermeshes.Similarly, parameters related to the iterative convergence of our flow solver islisted in Table 4.4. The linear solver settings are the same as those reported forthe flat plate test case and initial solution is set to uniform free stream values forall cases. The number of work units changes more slowly with mesh refinement forthird-order than for second- and fourth-order. This comes from the larger number oflinear iterations required for these cases and is related to the structure of Jacobianmatrix. For all the cases, the solver converges in a reasonable number of non-lineariterations and work units although the linear solver is not perfectly scalable.Figure 4.9 shows the convergence histories in terms of non-linear iteration countand number of work units for a mesh with 100, 352 control volumes.914.4. Resultsn.DoF CD,p CD,v CLReference 0.00606 0.00620 1.09102nd order6, 272 0.02748 0.00998 1.046625, 088 0.01095 0.00704 1.0778100, 352 0.00685 0.00623 1.0881401, 408 0.00614 0.00617 1.09083rd order6, 272 0.00019 0.00220 0.943225, 088 0.00380 0.00556 1.0817100, 352 0.00587 0.00615 1.0877401, 408 0.00607 0.00619 1.09044th order6, 272 0.00090 0.06200 0.581325, 088 0.00544 0.00582 1.0951100, 352 0.00615 0.00623 1.0905(a) Curvilinear coordinatesn.DoF CD,p CD,v CLReference 0.00606 0.00620 1.09102nd order6, 272 0.02221 0.00084 1.138425, 088 0.00443 0.00103 1.1765100, 352 0.00204 0.00207 1.1625401, 408 0.00321 0.00417 1.13533rd order6, 272 0.00443 0.00818 0.816625, 088 0.00740 0.00789 1.0199100, 352 0.00644 0.00662 1.0776401, 408 0.00612 0.00619 1.08994th order6, 272 0.03832 0.00477 0.685725, 088 0.01084 0.00862 1.0257100, 352 0.00749 0.00695 1.0724(b) Cartesian coordinatesTable 4.3: Convergence of drag and lift coefficients with mesh refinement for subsonicflow around NACA 0012Order n.DoFNo. of CPU Work Average Time (Sec.)iterations TimeUnitRes. Jac. LSLinear Non-linear (Sec.) Eval. Form. Solver2nd6, 272 837 36 52.72 2140 0.025 0.61 0.8225, 088 2521 43 301.72 3040 0.10 2.37 4.50100, 352 7281 46 2126.52 5260 0.40 9.88 35.643rd6, 272 662 50 88.10 2260 0.039 0.92 0.7725, 088 2210 56 486.05 3270 0.15 3.83 4.39100, 352 6638 47 2432.71 4010 0.60 14.70 36.044th6, 272 389 55 203.19 4000 0.05 1.31 2.7225, 088 1269 46 777.48 3810 0.20 5.29 11.12100, 352 11673 76 8306.67 10220 0.81 21.25 85.30Table 4.4: Convergence properties of turbulent RANS solver for subsonic flow over aNACA 0012 airfoil924.4. Results(a) Second-order solution(b) Fourth-order solutionFigure 4.7: Distribution of scaled turbulence working variable for subsonic flow overNACA 0012 on a mixed-element mesh (n.DoF = 25, 088)4.4.3 Transonic Flow Over a RAE 2822 AirfoilIn this test case, we consider a transonic, turbulent flow over a RAE 2822 airfoil.The free stream Mach number is Ma∞ = 0.73, the Reynolds number based on thechord is Rec = 6.5 × 106 and the angle of attack is α = 2.79◦. For this problem,which is the same as case 9 in the experimental investigation of Cook et al. [103],a shock wave appears on the upper surface of the airfoil and interacts with theturbulent boundary layer, forming a small recirculation behind the shock. As aresult, the convergence to the steady-state solution is more challenging, particularlyfor a higher-order discretization. This test case can examine the robustness of oursolver for non-trivial flow situations with the presence of shock and singularities.For this test case, we use a family of C-meshes consisting of triangles and quadswith three levels of refinement (four levels for second-order). The far-field boundaryis placed 20c away from the airfoil surface and adiabatic no-slip conditions are appliedon the wall. The flow is assumed to be fully turbulent with free stream turbulentviscosity of ν˜/ν∞ = 3.0.934.4. ResultsxCp0 0.2 0.4 0.6 0.8 1-6-5-4-3-2-101(a) Pressure coefficientxCf0 0.2 0.4 0.6 0.8 1-0.0100.010.020.032nd order (n.DoF = 100,352)4th order (n.DoF = 100,352)FUN3D (n.DoF = 917,504)(b) friction coefficientFigure 4.8: Distribution of surface pressure and friction coefficients for second- andfourth-order and comparison with FUN3D944.4. ResultsIterationsL 2  Norm Res.20 40 60 8010-910-710-510-310-11011032nd-order3rd-order4th-order(a) Non-linear iterationWork UnitL 2  Norm Res.102 103 10410-910-710-510-310-11011032nd-order3rd-order4th-order(b) Work unitsFigure 4.9: Convergence history of finite volume solver for subsonic flow over NACA0012 for n.DoF = 100,352954.4. ResultsFigure 4.10 shows the contours of pressure and scaled turbulence working variablefor the third-order solution on a mesh with 35, 840 control volumes. Note that theshock wave and small recirculation bubble behind the shock have been capturedautomatically without any change in the flux functions or using slope limiters. Inaddition, the turbulence working variable predicted by the SA-neg model is largestin the wake region and negative at the edge of boundary layer.The distributions of pressure and friction coefficients obtained by different ordersof discretizations are shown in Figure 4.11. There is a good match, including thelocation of the shock on the upper surface, between our computed pressure coefficientsand experimental data. However, the difference between second- and higher-ordervalues is less clear.Figure 4.12 shows the convergence of drag and lift coefficients for this test casewith mesh refinement. Similar to the other cases, the advantage of higher-ordermethods is observed as they converge faster to grid converged values. In addition,the error associated with fourth-order discretization is smaller on coarser meshescompared to second- and third-order discretizations.As described earlier, the initial condition for the previous test cases were the uni-form free stream states. However, the presence of shock wave in this case makes theconvergence of our iterative solver slower as a considerable number of iterations arerequired to locate the shock properly. On the other hand, starting from a good ini-tial guess, which gives the location of shock approximately, can help the convergenceand even robustness in some cases. As is common in the higher-order community, alower-order solution can be used as the initial state for the higher-order discretiza-tions. Figure 4.13 shows the convergence history of our solver for different conditionson a mesh with 35, 840 degrees of freedom. Note that starting directly from freestream values leads to about the same number of non-linear iterations for second-and third-order and non-convergence for fourth-order. Alternatively, starting fromthe converged second-order solution results in fast convergence in a small number ofiterations and work units. The iterative convergence properties tabulated for higher-order discretizations in Table 4.5 are for the conditions started from the second-ordersolution on the same mesh. Note that for the second-order cases, more iterations andwork units are necessary compared to the previous cases, which illustrates the chal-lenge in computing the solution of this case. To obtain a higher-order solution on thesame mesh, only a few more work units are required which are perfectly affordable964.4. Results(a) Pressure(b) Turbulence working variableFigure 4.10: Third-order solution of transonic flow over RAE 2822 on a mixed-elementmesh (n.DoF = 35, 840)974.4. ResultsxCp0 0.2 0.4 0.6 0.8-1-0.500.512nd order3rd order4th orderExperiment(a) Pressure coefficientxCf0 0.2 0.4 0.6 0.800.0020.0040.0060.0080.01(b) Friction coefficientFigure 4.11: Distribution of surface pressure and friction coefficients for differentorders and comparison with experiment984.4. Results1/(n.DoF)1/2CD0.005 0.01 0.0150.0140.01450.0150.01550.0160.01650.0170.01750.018(a) Drag coefficient1/(n.DoF)1/2CL0.005 0.01 0.0150.70.750.80.850.92nd-order3rd-order4th-order(b) Lift CoefficientFigure 4.12: Convergence of drag and lift coefficients with mesh refinement for tran-sonic flow around RAE 2822994.4. Resultsfor grid converged mesh resolutions.Order n.DoFNo. of CPU Work Average Time (Sec.)iterations TimeUnitRes. Jac. LSLinear Non-linear (Sec.) Eval. Form. Solver2nd8, 960 3000 60 123.24 3460 0.035 0.87 1.0735, 840 13853 125 1429.55 9460 0.15 3.41 7.32143, 360 27105 146 9126.80 14340 0.63 13.58 45.983rd8, 960 1164 44 104.77 1700 0.06 1.34 0.8835, 840 2864 35 415.20 1810 0.23 5.25 6.29143, 360 12167 59 4706.12 4990 0.94 21.24 55.654th8, 960 628 32 130.68 1560 0.084 1.91 2.0635, 840 1831 36 696.62 2220 0.31 7.55 11.35143, 360 6022 55 5840.08 4560 1.28 30.76 71.38Table 4.5: Convergence properties of turbulent RANS solver for transonic flow overa RAE 2822 airfoil4.4.4 High-lift Multi-element AirfoilAs our final test case, we consider turbulent flow over the three-element airfoil MDA30P30N configuration. The free stream Mach number is Ma∞ = 0.2, the Reynoldsnumber is Re = 9 × 106 and the angle of attack is α = 16◦. The geometry of themulti-element airfoil, which consists of a leading edge slat, a main center element anda trailing edge flap, has a number of sharp corners. This is considered a hard testcase in 2D because of the complex geometry, various flow structures and high angleof attack [104, 105].For this test case, we only demonstrate the capability of our flow solver for onesingle mesh with 45, 802 mixed cells, Figure 4.14. It is known that higher-ordersolution reconstruction on anisotropic meshes near singularities leads to solver failureeven for simple geometries [106]. As a result, the mesh consists of isotropic cells nearthe sharp corners. Other cells that are sufficiently far from the sharp corners andstill close to the walls have high aspect ratio. The new curvilinear coordinate is onlyconstructed for these cells with high aspect ratio; isotropic cells still use the Cartesiancoordinate system. Figure 4.14c shows the choice of reconstruction coordinate forcontrol volumes near the slat.Figure 4.15 shows the contours of the turbulence working variable for a fourth-1004.4. ResultsIterationsL 2  Norm Res.50 100 15010-910-710-510-310-1101103(a) Non-linear iterationWork UnitL 2  Norm Res.102 103 10410-910-710-510-310-1101103direct 2nddirect 3rd2nd to 3rddirect 4th2nd to 4th(b) Work unitsFigure 4.13: Convergence history of finite volume solver for transonic flow over RAE2822 for n.DoF = 35,8401014.4. Results(a)(b) (c)Figure 4.14: Mixed element mesh illustration for the high-lift three-element configu-rationFigure 4.15: Fourth-order solution of turbulence working variable over multi-elementairfoil on a mixed-element mesh (n.DoF = 45, 802)1024.4. Resultsorder discretization. Note that the slat wake continues all the way over the mainelement and flap and joins the wakes of the two other elements. Also, Figure 4.16compares the distribution of pressure coefficient obtained by different orders of dis-cretizations with experimental data available in Ref. [92]. As seen, there is a goodmatch between the computed pressure coefficients and experimental data.xC p0 0.2 0.4 0.6 0.8 1 1.2-12-8-402nd order3rd order4th orderExperimentFigure 4.16: Distribution of surface pressure coefficient over multi-element airfoilFinally, the convergence history of this test case is shown in Figure 4.17. Notethat order sequencing has been used to accelerate the convergence to the steady-state solution. For this test case, we solve the linear system of Equation 4.20 withthe tighter tolerance of ηl = 10−6 to obtain a more accurate solution update at eachnon-linear iteration as suggested by Ceze and Fidkowski [90]. The solver exhibitsremarkable robustness for this complex test case as the solutions of all orders ofaccuracy converge in a small number of non-linear iterations.1034.4. ResultsIterationsL 2  Norm Res.50 100 15010-910-710-510-310-11011032nd-order3rd-order4th-orderFigure 4.17: Convergence history of finite volume solver for flow over multi-elementairfoil104Chapter 5Adaptive Mesh Refinement andOrder EnrichmentAs noted, higher-order discretizations of fluid dynamic equations have received a greatdeal of attention due to their potential advantages in obtaining more accurate solu-tions with less cost. However, higher-order accuracy is only obtained in the smoothregions of the solution where there is no discontinuity in the solution or gradient. Inaerodynamic applications, several sources of discontinuities such as shocks, contactdiscontinuities, and the turbulence working variable at the edge of boundary layer(e.g., in SA model), can deteriorate the order of accuracy in non-smooth regions andalso cause solver failure in some cases. This motivates the idea of simultaneous meshrefinement (h-refinement) and order enrichment (p-enrichment) in higher-order com-pressible flows solvers. In other words, the order of solution approximation can beincreased in those parts of the domain where the solution is smooth whereas the meshresolution is enhanced in non-smooth regions in which a lower-order discretization isemployed.For compressible flow simulations, the combination of such a strategy with output-based adaptation (known as hp-adaptation) has been used by several DiscontinuousGalerkin (DG) solvers where the number of degrees of freedom increases rapidly withthe order of discretization. In this way, it is possible to optimally place degrees offreedom within a problem and achieve required accuracy with minimal costs [107, 108,44]. In addition, it has the advantage of solver robustness by employing a first-orderdiscretization of flow field variables near a shock wave without using slope limiters orartificial dissipation [109, 110]. hp-adaptive methods have been successfully appliedto the other variants of Galerkin-based schemes such as hybridized-DG [111] andPetrov-Galerkin [112, 113] and also the correction procedure via reconstruction (CPR)method [114].In contrast to compact schemes such as DG and its variants, higher-order finite105Chapter 5. Adaptive Mesh Refinement and Order Enrichmentvolume discretizations extend the stencil to obtain a more accurate estimate of thesolution. In such schemes, there is no coupling between the number of degrees offreedom and order of discretization as the number of control volumes remains con-stant on a mesh with a fixed number of elements regardless of the order of accuracy.Nevertheless, a higher-order polynomial for solution approximation requires a largernumber of derivatives to reconstruct at each iteration [54]. Also, in the case of usingan implicit time advance scheme for the convergence to the steady-state solution,higher-order methods lead to a larger number of non-zero entries in the Jacobian ma-trix. This increases the memory usage and computational cost due to a larger storagerequirement and also solving a denser linear system [49, 7] which must be precondi-tioned by a factorization method (e.g., Incomplete LU). As a result, hp-adaptationmethodology, which has originally been designed for compact discretization schemes,has certain advantages for being used in higher-order finite volume flow solvers. Inthis way, the extra number of derivatives are only reconstructed in those regionswhere needed and the Jacobian matrix becomes sparser with a smaller bandwidth.Furthermore, this strategy is capable of automatic limiting for flows with inherentdiscontinuities as we can start from a lower-order solution everywhere and increasethe order of polynomial in smooth regions while the mesh is refined near disconti-nuities. This can be viewed as an alternative to typical limiting approaches wherethe slope limiter designed to be active only near discontinuities may deteriorate so-lution accuracy even in smooth regions [55] and/or hamper the convergence to thesteady-state due to non-differentiability [115].This chapter discusses the development of an hp-adaptive unstructured finite vol-ume solver and its application to two-dimensional compressible flow problems rangingfrom inviscid flows governed by the Euler equations to viscous turbulent flows gov-erned by the RANS equations and the Spalart-Allmaras turbulence model. In par-ticular, we adapt the mesh resolution and discretization order based on an estimateof error in the computed solution at each level of refinement.The quality of the adaptive schemes relies on the accuracy of error estimates. Theconventional method of adaptation is based on feature-based criteria that highlightthe distinctive features such as shock waves and boundary layers in the flow field [116,117, 118]. This approach, which is simple and effective in some CFD applications,requires trial and error to determine the appropriate flow features and thus fails toprovide a general and robust error estimate [119, 120]. Alternatively, the adaptation106Chapter 5. Adaptive Mesh Refinement and Order Enrichmentprocedure can be carried out using residual-based approaches in which the truncationerror of the fluid flow quantities constructs the adaptation criteria [121, 122, 123].Moreover, output-based adaptation criteria based on solutions to the so-called adjoint(dual) problem, which is derived for an output of interest, has become very maturein recent years [124]. In this type of error estimates, the solution of the adjointproblem is multiplied by the local contribution of the truncation error estimate toprovide information about the interaction of the error in different components of thesolution [125] and subsequently those locations that require more resolution (typicallyh-refinement) for a more accurate estimate of the output quantity [126, 127]. It hasbeen shown that the inclusion of the adjoint solution improves the effectiveness ofthe adaptation procedure over the traditional residual-based approaches [120]. Thesetechniques have been used to perform adaptive mesh refinement (only h-refinement)in second-order unstructured finite volume methods for inviscid [128] and viscouslaminar flows [120]. In addition, they have been used in the context of DG methodsfor inviscid [41], viscous laminar [129] and turbulent RANS [130, 89, 88] simulations.In our work, we use residual-based and adjoint-based error estimation methods forhp-adaptation in our unstructured finite volume solver. In the former approach, weemploy a higher-order residual operator to estimate the truncation error of a lower-order discretization scheme. The magnitude of the estimated truncation error is usedas a local error indicator for h- or p-refinement. For the adjoint-based hp-adaptation,we compute the discrete adjoint solution obtained by a single linear system solve ateach refinement level and multiply the adjoint solution by the higher-order estimateof the truncation error to find a different error indicator. The adjoint problem isformed by evaluating one order higher operators based on the lower-order solution.In either approach, a certain fraction of control volumes contribute most to the totalerror are flagged for refinement. The decision for h-refinement versus p-enrichmentis based on local smoothness of the primal problem. It is worth mentioning thatwe allow for non-conforming interfaces (i.e., hanging nodes) in the mesh once we dothe h-refinement to be able to handle triangles and quadrilaterals in the same way,including meshes containing both.1075.1. Error Estimation5.1 Error EstimationAs mentioned, a reliable error indicator is required for any adaptation procedure.For residual-based adaptation, a local estimate of the error in flow field quantities isused whereas such an estimate is weighted by the solution of an adjoint problem foroutput-based adaptations. In our work, we use a higher-order operator to obtain anestimate of the truncation error.Consider the following continuous non-linear problem for which U is the exactsolution:R (U) = 0 (5.1)A lower-order discrete approximation of the exact solution, Up−1, satisfies the lower-order discrete non-linear residual as:Rp−1 (Up−1) = 0 (5.2)The exact truncation error for the lower-order discrete problem is defined as theamount by which the discrete lower-order solution does not satisfy the continuousPDE. Such an error property can be estimated by applying a higher-order discreteoperator, Rp, on the lower-order solution [125, 131, 132]:R (Up−1) ≈ Rp (Up−1) (5.3)Therefore, the local contribution of control volume i to the total error can be obtainedby the norm of the estimated truncation error in the corresponding control volume:ǫi = |Rp (Up−1)|i (5.4)For adjoint-based adaptation, the error involved in the computation of an out-put functional based on a lower-order solution, Jp−1 (Up−1), is estimated. For thispurpose, a higher-order estimate of the output functional based on a higher-order so-lution can be obtained by expanding a Taylor series expansion about the lower ordersolution projected to the higher-order space, Up−1→p. Note that in the context of fi-nite volume methods, the lower-order solution containing the control volume averagescan be identically mapped into the higher-order space (Up−1→p = Up−1). Therefore,1085.1. Error Estimationthe higher-order functional can be expanded as:Jp (Up) = Jp (Up−1) +(∂Jp∂Up)Up−1(Up − Up−1) + ... (5.5)where Jp (Up−1) and(∂Jp∂Up)Up−1are the higher-order functional and its sensitivity withrespect to the higher-order solution both evaluated at the lower-order solution state.Considering that a higher-order solution is not in hand, a Taylor series expansion canbe used for the higher-order residual to eliminate (Up − Up−1) in Equation 5.5:Rp (Up) = Rp (Up−1) +(∂Rp∂Up)Up−1(Up − Up−1) + ... (5.6)The higher-order solution satisfies its corresponding discrete operator, Rp (Up) = 0.So Equation 5.6 can be re-arranged to solve for the higher-order unknown solutionas:(Up − Up−1) ≈ −(∂Rp∂Up)−1Up−1Rp (Up−1) (5.7)Substituting Equation 5.7 into Equation 5.5 yields the following expression for theestimate of the error in the functional:Jp (Up)−Jp (Up−1) ≈ −(∂Jp∂Up)Up−1(∂Rp∂Up)−1Up−1Rp (Up−1) (5.8)Next, the higher-order discrete adjoint solution is defined as the variable Zp suchthat: (∂Rp∂Up)TUp−1ZP =(∂Jp∂Up)TUp−1(5.9)Note that both sides of Equation 5.9, which show the sensitivity of the higher-orderresidual and functional to the higher-order solution, are evaluated at the availablelower-order solution. The functional error can be re-written in terms of the discreteadjoint solution as:Jp (Up)− Jp (Up−1) ≈ −ZTp Rp (Up−1) (5.10)For adaptation purposes, we require a local error indicator. Therefore, the magnitudeof the contribution to the functional error from a particular control volume i can be1095.2. Adaptation Methodsapproximated as:ǫi =∣∣∣ZTp Rp (Up−1)∣∣∣i (5.11)5.2 Adaptation MethodsIn our hp-adaptive process, we use both error indicators of Equation 5.4 and 5.11 tocompare the effectiveness of the two estimates. Following Ceze and Fidkowski [108],we flag a certain percentage of control volumes, fadapt, with the largest error indicatorsfor refinement at each step of the adaptive procedure. For this purpose, a sorted listof the control volumes based on the value of their corresponding error indicator fromthe highest to the lowest is created. A loop over the list is executed and the controlvolumes are flagged for p-enrichment or h-refinement procedure where appropriateuntil the pre-specified target (NCV × fadapt) is reached. This ensures that only thecontrol volumes with the highest error magnitude are refined for highly non-uniformerror distributions and thus those parts of the computational domain which havenegligible contribution to the error in a functional output of interest are excludedfrom refinement.As described, accuracy enhancement in each control volume can be attained bydecreasing the mesh size (h) or increasing the order of accuracy (p). In what follows,we discuss the details of each method in isolation and then how they are combinedto give the hybrid hp-adaptive scheme for our unstructured finite volume solver.5.2.1 h-refinementTo be able to handle triangular and quadrilateral cells in the same way, we use a non-conforming mesh refinement framework where hanging nodes are allowed regardlessof the type of elements. In such a refinement, a triangle or quadrilateral is dividedinto four sub-cells. A triangle is refined using mid-point subdivision so that a nodeis placed at the mid-point of the faces whereas an additional node is required at thecenter of the cell for a quadrilateral. Figure 5.1 shows the refinement pattern forboth types of the cells in a four-to-one basis recognizing that the sub-cells createdthrough refinement inherit the approximation order from the original cell.The existence of hanging nodes complicates the computation of flux integrals andalso selection of reconstruction stencils. Any face connecting two vertices is identified1105.2. Adaptation Methods(a) Triangular cells(b) Quadrilateral cellsFigure 5.1: Schematic illustration of h-refinement pattern for 2D cellsas a unique face on which integration quadrature points are placed. In addition, theindices of the cells at the two sides of such a face is stored to retrieve the indices of thefirst layer of neighbors in the process of stencil selection. In this approach, a trianglewith one hanging node is represented by four faces (or similarly a quadrilateral withone hanging node is represented by five faces).To provide a smooth distribution of cell size throughout the mesh, we set up twostandard rules for the refinement. First, a cell is automatically flagged for refinementif all of its faces are split due to the refinement of the neighbors. This rule has theadvantage of limiting the number of hanging nodes to two and three for trianglesand quadrilaterals, respectively in addition to removing potential un-refined holesfrom the mesh. Figure 5.2 shows the implementation of this rule for triangular andquadrilateral cells. Note that the red cells are refined due to large error indicator ornon-smooth primal while the blue cell is refined just because all of its neighbors havebeen refined.The second rule enforces that all non-conforming faces have a 2:1 face length ratioby controlling the subdivision of half length non-conforming faces. If a fine cell with1115.2. Adaptation Methods(a) Triangular cells(b) Quadrilateral cellsFigure 5.2: Schematic illustration of the first rule for the refinement of 2D cellsa half length face needs to be refined, the coarse cell on the other side containingthe half length face must be divided as well to balance the face length ratio betweenadjacent cells. Figure 5.3 depicts the implementation of the second rule for triangularand quadrilateral cells. In this case, the red cells are those need to be refined due toa large error indicator or non-smooth primal while the blue cells are refined to keepthe 2:1 face length ratio.As described earlier, higher-order simulation on highly anisotropic meshes re-quires curving the interior faces of the mesh to prevent intersection with higher-orderboundary representations. When considering non-conforming meshes, special atten-tion must be paid to how to curve the cells that have hanging nodes. One tedious wayis taking the linear mesh, applying the necessary refinement and re-curving the mesh1125.2. Adaptation Methods(a) Triangular cells(b) Quadrilateral cellsFigure 5.3: Schematic illustration of the second rule for the refinement of 2D cellsby solving the linear elasticity equation at each level of refinement. Considering thatwe use the continuous Galerkin finite element method to solve the elasticity equation,this approach necessitates the enforcement of a constraint over non-conforming inter-faces to ensure that the faces defined from both sides conform to the same mapping.This comes from the fact that the refined faces on a non-conforming interface containdifferent number of degrees of freedom from each side (Figure 5.4). Alternatively,we can assume that the initial mesh represents the geometry sufficiently accuratelyand no further geometry information is added throughout the refinement. In this ap-proach, we only curve the initial conforming mesh and use the mapping informationto create higher-order non-conforming meshes at the next levels. The mapping fromeach linear cell of the initial mesh to its corresponding cubic cell is computed andstored. At each level of refinement, the top parent for each sub-cell is found. Thecubic mapping corresponding to the top parent is used to create the cubic sub-cells1135.2. Adaptation Methodswithin the cubic top parent cell without the loss of geometry due to non-conformingmapping. Figure 5.4: Different number of degrees of freedom for linear elasticity problem onhalf-length non-conforming face5.2.2 p-enrichmentThe implementation of p-enrichment is easier as it does not lead to any geometricalcomplexity and remains the same for triangular and quadrilateral cells. Recall thatthe primitive variables of the flow field are approximated throughout a control volumeby a Taylor series expansion about the reference point as:φRi (x, y) = φ|i +∂φ∂x∣∣∣∣∣i(x− xi) + ∂φ∂y∣∣∣∣∣i(y − yi) + ∂2φ∂x2∣∣∣∣∣i(x− xi)22+∂2φ∂x∂y∣∣∣∣∣i(x− xi) (y − yi) + ∂2φ∂y2∣∣∣∣∣i(y − yi)22+ ... (5.12)To enhance the accuracy by p-enrichment, the order of this polynomial is incrementedby one and a larger least-squares system with more control volumes in the stencil issolved to find the derivatives of the primitive variables at the reference point. Figure5.5 shows the p-enrichment of a cell from p = 1 (second-order) to p = 2 (third-order)schematically.We also apply some rules over the order of accuracy which mimic the rules usedin h-refinement for uniform distribution of cell size and lead to uniform distributionof order of accuracy throughout the mesh. First, if the order of accuracy in allneighboring cells is larger, we increment the order within the cell by one. Second, wemodify the order of polynomial in neighboring control volumes not to end up with ajump of more than one order between adjacent cells.1145.2. Adaptation MethodsFigure 5.5: Schematic illustration of p-enrichment for an arbitrary cell from second-order to third-orderThe number of integration quadrature points is determined based on the maxi-mum order (max (p+, p−)) of the two cells on each side of a face to ensure that thehigher-order reconstructed polynomial is integrated sufficiently accurately along thefaces.5.2.3 hp-refinementThe final step in the adaptation procedure is the decision between h-refinement andp-enrichment. As mentioned, this is typically decided based on the smoothness of theprimal solution [111, 110, 109] so that the mesh size is reduced near discontinuitiesand order is increased in the smooth regions of the primal solution. In our work,the smoothness is determined by an inter-element jump indicator designed for shockdetection in a DG solver [133]. This indicator has also been successfully used forthe decision between mesh refinement and order enrichment in other hp-adaptive DGsolvers [109, 110] and can be easily fit in our unstructured finite volume solver. Thevalue of the jump for an arbitrary quantity of the flow field, φ, in a cell (cell i) isgiven as:Si =1|∂Ωi|˛∂Ωi∣∣∣∣∣ φ+ − φ−12(φ+ + φ−)∣∣∣∣∣ ds (5.13)where φ+ and φ− are the reconstructed solutions at the two sides of a quadraturepoint, |∂Ωi| is the perimeter of the cell and ds is the infinitesimal length along thefaces of the cell. For smooth solutions, the jump in the reconstructed solution andthus Si are expected to be small whereas the indicator should return a value of O (1)1155.3. Numerical Resultsnear a discontinuity such as a shock wave since the jump and average flow propertiesare of the same order of magnitude. Therefore, the choice between whether to applyh-refinement or p-enrichment is made by:Si >1Kφ, h-refinementSi <1Kφ, p-enrichment(5.14)where Kφ is a constant required to be large enough to capture discontinuities properly.In this work, we apply the jump indicator only on pressure for non-turbulent flowsand on pressure and turbulence working variable of SA model for turbulent flows.The proper values of the jump indicator constants are found by several numericalexperiments. In our solver, the values of of 200 ≤ KP ≤ 400 for pressure and10≤ Kν˜ ≤ 20 for turbulence working variable seem to be effective choices for ourproblems. However, the best choices of constants in these ranges can be obtained byconsidering the flow regimes governed by important non-dimensional groups (Machand Reynolds number).5.3 Numerical ResultsThe hp-adaptation method proposed is evaluated for four different flow conditionsover the NACA 0012 geometry: inviscid subsonic and transonic flows governed bythe Euler equations, laminar subsonic flow governed by the Navier-Stokes equationsand turbulent subsonic flow governed by the RANS equations and SA-negative tur-bulence model. In each case, the efficiency of the residual-based and adjoint-basedhp-adaptive schemes are compared with second- and higher-order uniform refine-ments in terms of the number of degrees of freedom per equation and also CPU timerequired to obtain a grid-converged value of an output of interest. The presentednumber of degrees of freedom is the number of control volumes in each level and thewall clock time for the various test cases is based on simulations on a single core of ani7-4790 (3.60 GHz) CPU. In the adjoint-based adaptation method, the time includesthe solution of both the primal and adjoint whereas only the primal solve time isincluded for the uniform refinement and residual-based adaptations.Considering that our solver supports up to fifth-order discretizations (p = 4),once a control volume with the polynomial degree of pmax = 3 is flagged for extra1165.3. Numerical Resultsrefinement in smooth regions of the solution, we perform h-refinement instead ofp-enrichment since we need one order higher discrete operator for error estimation(Equations 5.4 and 5.11). However, this strategy fails for highly anisotropic mesheswith curvature and leads to solver failure. For such problems, the highly anisotropiccells close to curved walls (those using the curvilinear coordinate system as describedin Chapter 3) are taken out of the error indicator list once they reach p = 3. Instead,some other control volumes with lower discretization order and/or other regions ofthe mesh are flagged for p-enrichment or h-refinement as appropriate. In this way,we can still reduce the error without breaking the structure of anisotropic cells nearthe wall.5.3.1 Inviscid Subsonic FlowAs our first test case, we consider the inviscid subsonic flow over the NACA 0012with chord length of c = 1 and zero thickness at the trailing edge. The flow ischaracterized by the free stream Mach number of Ma∞ = 0.5 and angle of attackof α = 2◦ . For residual-based adaptation, we use the L2-norm of the truncationerror obtained by the application of one order higher residual operator on all fourcomponents of the solution for each control volume. For adjoint-based adaptation,we consider the target functional of the pressure drag coefficient given by:J (U) =ˆ∂Ωψ · (P nˆ) ds (5.15)where nˆ is the outward unit normal vector and ψ = 1C∞(cos α, sin α)Talong the wallboundaries and 0 everywhere else. C∞ is a normalized reference value defined as:C∞ =12γMa∞P∞c (5.16)In the hypothetical case of infinitely far outer boundary, the drag force for this isen-tropic inviscid flow must converge to zero. However, we place the outer boundary100c away from the wall boundaries and thus the pressure drag coefficient convergesto a finite small value. The adaptation process starts with the second-order solutionon a baseline coarse mesh with 2, 776 triangular cells. At each level of the adap-tive procedure, fadapt = 15% of control volumes with the largest error indicators areflagged for p-enrichment or h-refinement. Considering that the solution of this prob-1175.3. Numerical Resultslem is smooth almost everywhere in the computational domain, we choose the upperbound of the pressure smoothness indicator, KP = 200, for the decision betweenp-enrichment and h-refinement.The efficiency of the hp-adaptive refinement schemes is shown in Figure 5.6 incomparison with uniform refinement of different discretization order. A referencevalue for this problem is obtained by a fourth-order solution on a mesh with 197, 891cells. The adjoint-based adaptation outperforms the other refinement strategies asit converges to the grid converged drag coefficient after a few cycles with a smallernumber of degree of freedom. The change in the values of the functional becomessmall after 6 cycles of adaptations and remains almost constant in the last three con-secutive cycles. The residual-based adaptation converges to the same value but aftermore refinement cycles and with a larger number of control volumes. In terms of theCPU time, the adjoint-based method takes slightly longer due to the extra adjointsolver. For this problem with a smooth solution, the fourth-order uniform refinementproduces reasonably accurate output values with sufficient number of degrees of free-dom which takes longer in term of CPU time. However, the second- and third-orderuniform refinement procedures are far from the grid converged value except for thefinest mesh with about 80, 000 control volumes.Figure 5.7 illustrates the baseline mesh and also the order of accuracy and meshresolution after 6 cycles of the adjoint-based hp-adaptation for this test case. Notethat the region near the leading edge in which the flow experiences high gradientstowards the stagnation point is refined the most considering that the resolution ofthe initial mesh is not sufficient there. It is also worth-mentioning that the order ofaccuracy far from the wall remains equal to 2 which leads to less numerical complexitycompared to uniformly higher-order discretizations .In addition, the contours of the Mach number on the two meshes with the pre-scribed order of accuracy for each control volume is shown in Figure 5.8. As expected,the hp-adapted mesh produces a smoother distribution of the Mach number near theairfoil.5.3.2 Inviscid Transonic FlowAs our second test case, we consider the transonic inviscid flow around the NACA0012 with free stream Mach number of Ma∞ = 0.8 and angle of attack of α = 1.25◦.1185.3. Numerical ResultsDegrees of FreedomC D103 104 10500.511.522.533.544.55 hp-adjointhp-residualuniform-2nduniform-3rduniform-4thReference10 - 4(a) Drag coefficient evolution with respect to degrees of freedomCPU Time (sec.)C D100 101 102 10300.511.522.533.544.55hp-adjointhp-residualuniform-2nduniform-3rduniform-4th10 - 4(b) Drag coefficient evolution with respect to CPU timeFigure 5.6: Convergence of drag coefficient for inviscid subsonic test case (Ma∞ = 0.5and α = 2◦)1195.3. Numerical Results(a) Baseline mesh with 2, 776 triangular cell (purely second-order)(b) adjoint-based hp-adapted after 6 adaptation cycles (second- to fourth-order)Figure 5.7: Mesh resolution and discretization order for inviscid subsonic test case(Ma∞ = 0.5 and α = 2◦)1205.3. Numerical Results(a) Initial mesh with purely second-order discretization(b) adjoint-based hp-adapted after 6 levels of adaptationFigure 5.8: Contours of Mach number for inviscid subsonic test case (Ma∞ = 0.5 andα = 2◦)1215.3. Numerical ResultsFor this test case, a strong shock and a weak shock appear on the upper and lowersurfaces of the airfoil, respectively. Therefore, this is a good example to demonstratethe capability of the proposed hp-adaptation method for shock capturing. Similar tothe inviscid subsonic case, we use the L2-norm of the truncation error for residual-based adaptation and pressure drag coefficient as the functional of interest for theadjoint-based method. The other adaptation parameters are chosen as fadapt = 15%and KP = 300. As described, the hp-adaptation framework can be an alternative toslope limiting procedure for compressible problems with shocks and so can be usedto enhance the robustness of an unstructured finite volume solver. This suggests weallow for first-order approximations near the shock as traditional slope limiters reducethe order of accuracy to one close to singularities. Although it is possible to startfrom a first-order solution everywhere and increasing the order in smooth regions ofthe flow field, this approach leads to the over-refinement of the initial mesh due tothe the non-smooth behavior of the first-order solution. This is also repeated in thenext cycles since the order is not increased in all of the h-refined cells (with order =1). Instead, we start from a purely second-order solution on the same baseline meshas the subsonic test case but we reduce the order of discretization by one (providedthat the order is larger than 1) whenever a cell is flagged for h-refinement due tonon-smooth primal solution. In this way, we can attain the first-order approximationnear the singularities after the first cycle.For uniform refinement studies, we only employ the second- and third-order dis-cretizations with the slope limiter of Venkatakrishnan [134] and avoid the fourth-orderdiscretization since it suffers from stability issues for such a problem with a strongshock. Figure 5.9 compares the efficiency of different refinement scenarios for thiscase. A reference drag coefficient is obtained by a second-order solution of this prob-lem on an unstructured mesh with 308, 470 cells. It is evident that the third-orderdiscretization is not converging to a fixed value as the slope limiter degrades theaccuracy of the solution. Furthermore, the tunable parameter in the slope limiteris adjusted for each mesh to get the solution to converge and this exacerbates theasymptotic convergence to a single value. Therefore, increasing the order of accuracyfor problems with locally non-smooth solutions does not lead to more accurate an-swers as expected. On the other hand, the adjoint-based hp-adaptive method reachesthe grid converged value after a small number of adaptation cycles so that the changein the drag coefficient becomes less than 0.2% in the last three consecutive cycles.1225.3. Numerical ResultsNote that such a good estimate of the drag coefficient is found on an hp-adapted meshwith about 20, 000 control volumes. In term of the CPU time, the adjoint-based adap-tation obviously outperforms uniform refinement considering the fact that the dragpredicted by the second-order scheme on the finest mesh is not still sufficiently closeto the reference value and takes about 5 times longer. The residual-based adaptationis also approaching the reference value but with a larger number of refinement cycles,degrees of freedom and CPU time compared the adjoint-based.The difference between the drag coefficients predicted by the two adaptation errorindicators can be explained by the refinement pattern obtained from each one. Figure5.10 compares the mesh resolution and local order of accuracy in the last cycle of theadjoint-based and residual-based adaptations. The residual-based indicator mainlyfocuses on the strong shock of the upper surface with minimal refinement of theleading and trailing edges. This behavior was previously reported by Woopen et al.[135] in residual-based mesh adaptation for their hybridized discontinuous Galerkin(HDG) method. In their study, the drag coefficients obtained on the residual-basedadapted meshes differ by a constant amount from the values of the adjoint-basedadapted meshes after a few refinement cycles. Conversely, the adjoint-based errorindicator exhibits a more accurate refinement pattern as the the area close to the weakshock on the lower surface is refined properly in addition to the sufficient refinementof the leading and trailing edges. This implies that weighting the local residuals withthe discrete adjoint solution provides a more accurate error indicator that capturesall important features of the flow field during the refinement process. It should benoted that the first-order approximation is maintained close to the upper surfacestrong shock as magnified in Figure 5.10a which is ideal for regions with singularsolution. Such a strategy is advantageous as it results in the efficiency of the solverin obtaining a grid converged output value with small number of degrees of freedomand also its robustness by employing a low-order solution in the singular regions ofthe flow field.Figure 5.11 shows the contours of Mach number close to the airfoil for the initialand final adjoint-based adapted meshes of this test case. As expected, the two shockssharpen as we proceed and increase the mesh resolution in their vicinity. Also, theuse of first-order approximation in these regions prevents any overshoots in the finalsolution. In the smooth regions of the flow, the solution becomes smoother as weprovide extra resolution by p-enrichment.1235.3. Numerical ResultsDegrees of FreedomC D103 104 1050.0220.0240.0260.028 hp-adjointhp-residualuniform-2nduniform-3rdReference(a) Drag coefficient evolution with respect to degrees of freedomCPU Time (Sec.)C D100 101 102 1030.0220.0240.0260.028 hp-adjointhp-residualuniform-2nduniform-3rd(b) Drag coefficient evolution with respect to CPU timeFigure 5.9: Convergence of drag coefficient for inviscid transonic test case (Ma∞ = 0.8and α = 1.25◦)1245.3. Numerical Results(a) Adjoint-based adaptation(b) Residual-based adaptationFigure 5.10: Final hp-adapted mesh and order for inviscid transonic test case (Ma∞ =0.8 and α = 1.25◦)1255.3. Numerical Results(a) adjoint-based hp-adapted after 1 adaptation cycle (first adaptedmesh)(b) adjoint-based hp-adapted after 9 adaptation cycles (final adaptedmesh)Figure 5.11: Contours of Mach number for inviscid transonic test case (Ma∞ = 0.8and α = 1.25◦) 1265.3. Numerical ResultsxPressure0.6 0.65 1Level 3Level 5Level 7Level 9Figure 5.12: Comparison of pressure profiles near the upper surface strong shockbetween adjoint-based hp-adapted meshes for inviscid transonic test case at y = 0.3Figure 5.12 compares the pressure profiles near the strong shock at y = 0.3 be-tween different levels of adjoint-based hp-adaptation. On the initial mesh, the shockspans a wide distance with several jumps in the magnitude of pressure. However,it gradually becomes thinner and sharper such that the shock computed on the lastlevel is around 50 times thinner than the first one. This plot reveals the advantageof hp-adaptation in shock capturing in flows with singularities.Also, Figure 5.13 compares the pressure profile at y = 0.3 between adjoint-basedand residual-based methods after the final cycle of adaptation. Note that the residual-based adaptation refines the region near the strong shock more than the adjoint-basedas expected. This can be understood by the larger number of piecewise constantpressure values (first-order solutions) in the vicinity of the shock.5.3.3 Laminar Subsonic FlowNow we turn our attention to viscous laminar subsonic flow around the NACA 0012airfoil with a free stream Mach number ofMa∞ = 0.5, Reynolds number of Re = 50001275.3. Numerical ResultsxPressure0.64 0.65 0.660. 5.13: Comparison of pressure profiles near the upper surface strong shockbetween adjoint-based and residual-based hp-adaptation at y = 0.3with constant dynamic viscosity and angle of attack of α = 1◦. In this flow, athin boundary layer appears on the surface of the airfoil and is followed by a wakedownstream of the trailing edge. Similarly, the L2-norm of the truncation error isused for residual-based adaptation, whereas the drag coefficient obtained as the totalof the pressure and viscous drag coefficients is employed for the adjoint-based method.The total drag functional is defined as:J (U,∇U) =ˆ∂Ωψ · (P nˆ− τnˆ) ds (5.17)in which nˆ is the outward unit normal vector, τ is the viscous stress tensor andψ = 1C∞(cos α, sin α)Talong the wall boundaries and 0 everywhere else. We startthe adaptations with a second-order solution on an unstructured mesh with 10, 797triangular cells. The adaptation parameters for this problem with smooth solutionare set to fadapt = 10% and KP = 200.Figure 5.14 compares the efficiency of the hp-adaptation methods with uniformrefinement of second-, third- and fourth-order discretizations. For this purpose, the1285.3. Numerical Resultsconvergence of the total drag coefficient versus the number of degrees of freedomand CPU time is plotted. For this problem, a reference drag coefficient is found bya fourth-order discretization on a mesh with 308, 470 cells. The advantage of hp-adaptation over uniform refinement is very clear for this test case with thin shearlayers: the hp-adapted meshes reach grid convergence with about as many degreesof freedom as the first uniform refinement. At this resolution, the second-order dis-cretization on the uniformly refined mesh has not yet begun to improve the results.The higher than second-order schemes (particularly fourth-order) converge faster asexpected for this problem with smooth solution; however, the CPU time required toget close to the reference drag coefficient is more than 6 times larger compared withhp-adaptation. Similar to the previous test cases, the adjoint-based adaptation ismore successful than the residual-based. Even though the residual-based convergesto the same output value, it adds a number of unnecessary degrees of freedom in eachcycle that do not improve functional accuracy. Note that the increase in the numberof degrees of freedom for this problem typically originates from those fourth-ordercontrol volumes which need more resolution and thus are refined via h-refinement.On the other hand, the adjoint-based adaptation indicator optimally increases theresolution to be able to deliver a very accurate functional for which the adaptationis performed. It is worth mentioning that the time required for solving the discreteadjoint problem is a small portion of the total CPU time. Interestingly, the outputvalue becomes almost constant in the last three cycles of the adjoint-based adaptationwhich highlights its efficiency.Figure 5.15 depicts the mesh resolution and local order for the control volumesof the initial mesh and also adjoint-based hp-adapted mesh after 7 cycles. As re-quired, the order of accuracy increases in the boundary layer and wake regions viap-enrichment. Since some of the control volumes in these regions still contribute tothe error the most even after the p-enrichment, we apply h-refinement so that we canreduce the error in the next cycles effectively and converge to the reference outputvalue as fast as possible. To provide such a resolution with uniform refinement, werequire considerably more cells in the mesh, most of which are far from the airfoiland do not contribute to the total error. Note that the mesh size and order do notchange sufficiently far from the viscous dominated regions in the hp-adapted mesh.Finally, the contours of Mach number on the initial and final hp-adapted meshesare shown in Figure 5.16. The longer wake region and smoother velocity distribution1295.3. Numerical ResultsDegrees of FreedomC D, tot104 105 1060.0540.0560.0580.060.0620.064hp-adjointhp-residualuniform-2nduniform-3rduniform-4thReference(a) Drag coefficient evolution with respect to degrees of freedomCPU Time (Sec.)C D, tot101 102 1030.0540.0560.0580.060.0620.064hp-adjointhp-residualuniform-2nduniform-3rduniform-4th(b) Drag coefficient evolution with respect to CPU timeFigure 5.14: Convergence of drag coefficient for laminar subsonic test case (Ma∞ =0.5, Re = 5000 and α = 1◦)1305.3. Numerical Results(a) Baseline mesh with 10, 797 triangular cells(b) adjoint-based hp-adapted after 7 adaptation cycles (second- to fourth-order)Figure 5.15: Mesh resolution and discretization order for laminar subsonic test case(Ma∞ = 0.5, Re = 5000 and α = 1◦)1315.3. Numerical Results(a) Initial mesh with purely second-order discretization(b) adjoint-based hp-adapted final meshFigure 5.16: Contours of Mach number for laminar subsonic test case (Ma∞ = 0.5,Re = 5000 and α = 1◦)near the airfoil on the final mesh are clearly visible.5.3.4 Turbulent Subsonic FlowAs our final test case, we compute the turbulent subsonic flow around the NACA0012 airfoil described in Chapter 4. This flow is characterized by a free stream Machnumber of Ma∞ = 0.15, Reynolds number of Re = 6 × 106 and angle of attack ofα = 10◦. As demonstrated by the previous test cases, the adjoint-based method issuperior to the residual-based adaptation. Therefore, we only consider the adjoint-based error indicator for hp-adaptation. The lift coefficient is used as the output1325.3. Numerical Resultsfunctional of interest and is defined as:J (U,∇U) =ˆ∂Ωψ · (P nˆ− τnˆ) ds (5.18)with ψ = 1C∞(− sin α, cos α)Talong the wall boundaries and 0 everywhere else. Thecomputations are performed on the second family of structured quadrilateral meshesavailable on the NASA turbulence modeling resource [85]. For simplicity, we assumethat the dynamic viscosity is constant.For uniform refinement, we consider second- to fourth-order discretization on threenested meshes with 14, 336, 57, 344 and 229, 376 control volumes. The adjoint-basedhp-adaptation starts with the second-order solution on the coarsest mesh (Figure5.17). As described earlier, those highly anisotropic cells near the wall for which thecurvilinear coordinate is employed come out of consideration for further refinementonce they reach p = 3. The other adaptation parameters are set to KP = 200,Kν˜ = 20 and fadapt = 10%.Figure 5.17: The coarsest quadrilateral mesh for turbulent subsonic flowFigure 5.18 compares the efficiency of the adjoint-based hp-adaptation with uni-form refinement of different discretization orders. The reference lift coefficient is thevalue reported by the NASA turbulence modeling resource for these meshes. Again,adaptive refinement is considerably more successful than uniform refinement as itrequires 5 times smaller number of control volumes to obtain the grid converged out-put. This is expected since adaptive refinement places sufficient resolution in theboundary layer and wake regions without unnecessary refinement of other regions.1335.3. Numerical ResultsSuch a resolution can only be obtained on super fine meshes in the case of uniformrefinement. More importantly, the CPU time on the final adapted mesh is 6, 9 and24 times smaller than that of second-, third- and fourth-order discretizations on thefinest mesh, respectively. This is also expected since the extra derivatives required forhigher-order approximations are only computed in a small fraction of control volumesand also the Jacobian matrix has fewer entries.Figure 5.19 shows the final adapted mesh and also order of accuracy for thecontrol volumes close to the airfoil. Note that the adaptive refinement focuses onthe boundary layer, the wake behind the trailing edge and the stagnation line at theleading edge. It is also worth mentioning that the order of accuracy at the edge ofthe boundary layer and wake regions where the derivative of the turbulence workingvariable is discontinuous does not increase even though the mesh is refined.Figure 5.20 compares the turbulence working variable between the initial and finalhp-adapted solutions. The viscous dominated parts of the flow, including the wallboundary layer and wake are captured properly. Moreover, the region at the edgeof viscous regions where ν˜ is negative becomes thinner as we refine. The smoothcontours of the turbulence working variable and the peak value in the wake regionshow the improvement in the solution of the hp-adapted mesh compared to the initialmesh.1345.3. Numerical ResultsDegrees of FreedomC L104 105 1061. Lift coefficient evolution with respect to degrees of freedomCPU Time (Sec.)C L102 103 104 1051. Lift coefficient evolution with respect to CPU timeFigure 5.18: Convergence of lift coefficient for turbulent subsonic test case (Ma∞ =0.15, Re = 6× 106 and α = 10◦)1355.3. Numerical Results(a) Mesh(b) Local order of accuracyFigure 5.19: Illustration of final hp-adapted mesh for turbulent subsonic flow1365.3. Numerical Results(a) Initial mesh with purely second-order discretization(b) adjoint-based hp-adapted final meshFigure 5.20: Contours of turbulence working variable for turbulent subsonic test case(Ma∞ = 0.15, Re = 6× 106 and α = 10◦)137Chapter 6Concluding Remarks6.1 SummaryThis thesis described RANS simulation of turbulent compressible aerodynamic flowsusing a higher-order unstructured finite volume discretization of governing equations.In addition an hp-adaptation framework was developed to improve the efficiency androbustness of the flow solver for the two dimensional problems that are of interest incomputational aerodynamics.The solution of turbulent flows necessitates the use of highly anisotropic meshessince the gradients of flow properties normal to the shear layers are considerablylarger than the gradients along them. Considering that all higher-order accuratesolvers must account for boundary curvature, we curved the interior faces for highlyanisotropic meshes to prevent faces from intersecting near curved boundaries. Forthis purpose, we used an elasticity analogy to project the boundary curvature intothe interior assuming that the domain to be meshed behaves as an elastic solid. Thegeometrical features of the higher-order mesh were calculated using the coordinatetransformation from the curved cells to the original straight cells.Moreover, the higher-order cell-centered finite volume solution reconstruction pro-cedure was revisited on highly anisotropic meshes over curved surfaces to address twoknown issues: poor conditioning of the least-squares system and poor accuracy of re-constructed values on anisotropic meshes over curved surfaces. For highly anisotropicmeshes, the least-squares system suffers from ill-conditioning and this issue is exac-erbated for higher-order computations as the presence of the higher moments of arearesults in a considerable difference in the order of magnitude of matrix columns. Wealso investigated the accuracy and conditioning of the k-exact reconstruction on awide range of unstructured anisotropic meshes from straight triangles to more gen-eral cases used for turbulent simulations over aerodynamic configurations. For thispurpose, we manufactured anisotropic functions that mimic the characteristics of real1386.2. Conclusionsanisotropic solutions in aerodynamic problems.For turbulent flow simulations, the negative variant of the Spalart-Allmaras tur-bulence model, developed for higher-order discretizations, was fully coupled to thesystem of RANS equations. The convective fluxes of the coupled system were ob-tained by the Roe-Pike flux function and the gradients used in the viscous fluxes werecalculated by the average of the two reconstructed gradients plus a jump term hav-ing a stabilization property. An implicit time stepping method combined with a linesearch algorithm was employed to solve the discrete system of equations efficiently.At each non-linear iteration, the GMRES method preconditioned with ILU(k) wasused to solve the linear system. We also presented our numerical results for simpleand complicated turbulent flow problems in 2D.Finally, an hp-adaptation method was proposed to enhance the capabilities of thesolver and improve its efficiency. Considering that the higher-order approximationsare only advantageous in smooth regions of the solution, we employed either h- orp-refinement based on the smoothness of the solution. This procedure was carriedout for those cells that contribute to the total error the most using residual-based andadjoint-based error indicators. For mesh refinement, we allowed for non-conforminginterfaces in the mesh to be able to handle triangles and quadrilaterals in the same.We presented the result of our adaptive solver for different compressible flow problemsincluding inviscid, laminar and turbulent cases.6.2 ConclusionsTo conclude this thesis, a short review of the highlights observed in each step ispresented.For solution reconstruction on anisotropic meshes, our results showed that columnscaling using the maximum entry of each column and performing the reconstructionprocedure along principal axes alleviates the poor conditioning of the least-squaressystems significantly. As a result, the condition number of the scaled matrix in theprincipal coordinate system can be used as the measure of conditioning.Furthermore, our reconstruction results demonstrated that the unweighted sys-tem produces more accurate reconstruction coefficients in the case of straight mesheswhere the orientation of principal axes of control volumes change only slightly. More-over, the condition number of the scaled LS system is small and independent of mesh1396.2. Conclusionssize and aspect ratio for the unweighted system. For curved meshes, the recon-structed values using the unweighted LS method in its traditional form do not givethe expected order of accuracy. Although adding a weight function to the LS systemresolves this issue, it leads to large relative errors for the reconstructed solution andpoor estimate of normal derivatives. Also, the condition number of the scaled LSsystem grows with aspect ratio and mesh refinement in both cases particularly forfourth-order reconstruction. Instead, we proposed a new baseline where a curvilinearcoordinate system aligned with the wall is constructed and used for reconstruction.The accuracy results in the proposed coordinates have small values of error and matchthe expected order of accuracy with or without geometric weights. Nevertheless, theunweighted LS errors are noticeably smaller and thus outperform the weighted LSmethod. Likewise, the unweighted LS problem is well-conditioned as the conditionnumbers are small and do not change with mesh size and aspect ratio. Finally, weextended our analysis to general meshes such as those used for turbulent simulationsover NACA 0012 airfoils. These meshes are comprised of anisotropic regions nearthe curved walls and in wake region and also isotropic triangles sufficiently far fromthese regions. We presented a method for separating the mesh into different regionswhere different coordinate systems are required for accurate reconstruction. More-over, we demonstrated the advantage of the new curvilinear coordinate system inapproximating anisotropic solutions in the boundary layer region of these meshes.In terms of turbulent flow simulations, numerical results were presented for severaltest cases: the turbulent flow over the flat plate, the subsonic flow over the NACA0012 airfoil, the transonic flow over the RAE 2822 airfoil and the turbulent flow overthe three-element airfoil. The solutions of these case were verified by comparisonagainst the solutions of well-established codes or experimental data. For the flatplate and single airfoil test cases, the accuracy advantage of higher-order discretiza-tions was shown by obtaining a more accurate solution on a coarser mesh. It was alsoshown that solution reconstruction in the newly developed curvilinear coordinates re-sults in more accurate values of drag and lift coefficients for the curved geometries.In addition, the convergence to the steady-state solution from a free stream uniformsolution was shown to be fast and efficient for the flat plate and subsonic airfoil testcases for all order of accuracy (second- to fourth-order). However, for the more chal-lenging test case of transonic airfoil in which a shock wave interacts with turbulentboundary layer, starting from a lower-order solution helps the convergence signifi-1406.3. Recommended Future Workcantly. We also used order sequencing to accelerate the convergence for the difficultcase of multi-element airfoil. The results of this test case were presented to show therobustness of our solver for complex geometries with several flow structures.In the previous chapter, the results of the hp-adaptation algorithm developed foran unstructured finite volume solver in this thesis were presented for the subsonicinviscid, transonic inviscid, subsonic laminar and subsonic turbulent flows aroundthe NACA 0012 airfoil. In each case, uniform refinement of second- and higher-ordermethods was compared against residual-based and adjoint-based adaptations in termsof the number of degrees of freedom and CPU time required to reach a referencevalue. In all of the cases, it was shown that the adjoint-based hp-adaptive methodoutperforms the other refinement procedures and is more successful in capturing theflow features that affect the accuracy of an output functional. Moreover, it was shownthat the hp-adaptive method is capable of shock capturing and automatic limiting offlows with discontinuities (such as the transonic inviscid case) by mesh refinement andlowering the order of accuracy to one in those regions, simultaneously. For such flowswhere increasing the order of accuracy in non-smooth regions of the solution does notlead to more accurate solutions or sometimes causes solver failure, such a strategyis superior. Also, the advantage of the adaptive methods were sufficiently clear inthe case of viscous laminar and turbulent cases where thin shear layers are presentin the solution. For these cases, the adaptive methods place enough resolution inthe viscous dominated regions and thus provide a good estimate of functional valuesas opposed to the uniform refinements where many control volumes are needed toresolve those regions accurately.6.3 Recommended Future WorkThis thesis demonstrated an hp-adaptive unstructured finite volume solver for RANSsimulations over 2D configurations with aerodynamic applications. There are a num-ber of areas to which this research could be extended.1. The linear elasticity approach used to curve the interior faces of meshes withhighly anisotropic cells was successful for the cases discussed in the thesis. How-ever, for complex configurations with high curvature in several regions, a morepowerful curving strategy is required. There are some other analogy approaches1416.3. Recommended Future Worksuch as those that use non-linear elasticity equations [75] or incorporate ther-mal stress terms into the state equation [136] to enhance the quality of thehigher-order mesh.2. Although the non-linear convergence of our solver is quite good in most cases,the CPU time spent by the linear system solver does not scale linearly withthe number of degrees of freedom. Such behavior could be obtained by usinga more powerful preconditioner such as hp-multigrid with agglomeration anddirectional coarsening in the boundary layer region where highly anisotropicmeshes exist [137].3. Higher-order simulation of 2D turbulent flows is challenging as described andrequires several modifications. Considering that the ultimate goal in the nu-merical simulation of flows with engineering applications is the accurate andefficient computation of 3D turbulent flows, this work could be extended to 3Dproblems later by generalizing the ideas explained in this thesis for 2D problems.However, this requires several considerations:• The curvilinear coordinate system developed in this thesis must be ex-tended to 3D. This can be achieved by introducing the unit normal direc-tion, nˆi, in the same way as 2D by finding the unit vector that connectsthe closest point on the wall to the reference location of the target controlvolume. For the two other axes (we call them tˆi and bˆi), we can use thetwo principal directions of the polygon obtained by the intersection of thecross plane (normal to nˆi) and the target control volume. The first axis,tˆi, is defined as the one that forms the minimum angle with the projectionof the free stream velocity vector on the cross plane and the second axisis defined so that bˆi = tˆi × nˆi.• For turbulence modeling in 3D, the negative variant of SA model could stillbe used. However, the current status of computation power and resourcesallows one uses LES or hybrid RANS/LES for those flows encountered incomputational aerodynamics.• Also, a very powerful solution strategy particularly for the solution of thelinear system is a need. For 3D problems, the size of the linear system isvery large and thus better preconditioning techniques are required. Note1426.3. Recommended Future Workthat the solution of 3D turbulent problems requires parallel processingwith many nodes and thus a scalable linear solver cannot be avoided.4. In the current setting of our hp-adaptive algorithm, second-order discretizationis used in those regions that do not affect the functional accuracy. It mightbe advantageous to add a de-refinement strategy to the adaptive methods todecrease the number of cells in such regions and make the solver even lessexpensive.5. Finally, the hp-adaptive technique proposed here must be extended to 3D prob-lems since obtaining grid-converged output values is not affordable by uniformrefinement in 3D.143Bibliography[1] C.R. Illingworth. Some solutions of the equations of flow of a viscous com-pressible fluid. Proceedings of the Cambridge Philosophical Society, 46:469–478,1950.[2] Doug Pagnutti. Anisotropic adaptation: Metrics and meshes. Master’s thesis,The University of British Columbia, Department of Mechanical Engineering,2008.[3] F.R. Bailey. High-end computing challenges in aerospace design and engineer-ing. In Third International Conference on Computational Fluid Dynamics (In-vited Lecture), 2004.[4] P.R. Spalart. Topics in detached-eddy simulation. In Third International Con-ference on Computational Fluid Dynamics (Invited Lecture), 2004.[5] V. Venkatakrishnan. A perspective on unstructured grid flow solvers. Technicalreport, ICASE 95-3 NASA, 1995.[6] S. De Rango and D. W. Zingg. Higher-order spatial discretization for turbulentaerodynamic computations. American Institute of Aeronautics and Astronau-tics Journal, 39(7):1296–1304, July 2001.[7] Amir Nejat and Carl Ollivier-Gooch. A high-order accurate unstructured finitevolume Newton-Krylov algorithm for inviscid compressible flows. Journal ofComputational Physics, 227(4):2592–2609, 2008.[8] Stephen B. Pope. Turbulent Flows. Cambridge University Press, 2000.[9] Parviz Moin and John Kim. Tackling turbulence with supercomputers. Scien-tific American, 276(1):46–52, 1997.144Bibliography[10] Sergio Pirozzoli and Francesco Grasso. Direct numerical simulation of impingingshock wave/turbulent boundary layer interaction at M = 2.25. Physics ofFluids, 18(6):065113, 2006.[11] Philipp Schlatter and Ramis Örlü. Assessment of direct numerical simulationdata of turbulent boundary layers. Journal of Fluid Mechanics, 659:116, 2010.[12] Xuliang Liu, Shuhai Zhang, Hanxin Zhang, and Chi-Wang Shu. A new class ofcentral compact schemes with spectral-like resolution I: Linear schemes. Journalof Computational Physics, 248:235–256, 2013.[13] J.-B. Chapelier, M. de la Llave Plata, F. Renac, and E. Lamballais. Evaluationof a high-order discontinuous Galerkin method for the DNS of turbulent flows.Computers & Fluids, 95:210–226, 2014.[14] Nicholas J. Georgiadis, Donald P. Rizzetta, and Christer Fureby. Large-EddySimulation: Current capabilities, recommended practices, and future research.AIAA Journal, 48(8):1772–1784, 2010.[15] Roland Bouffanais. Advances and challenges of applied large-eddy simulation.Computers & Fluids, 39:735–738, 2009.[16] Barrett S. Baldwin and Harvard Lomax. Thin-layer approximation and alge-braic model for separated turbulent flows. In Proceedings of Sixteenth AIAAAerospace Science Meeting, 1978. AIAA Paper 78-257.[17] T. Cebeci. Analysis of turbulent boundary layers. Elsevier, 2012.[18] Barrett S. Baldwin and Timothy J. Barth. A one-equation turbulence transportmodel for high Reynolds number wall-bounded flows. NASA TM-102847, 1990.[19] P.R. Spalart and S. R. Allmaras. A one-equation turbulence model for aerody-namic flows. In Proceedings of the Thirtieth AIAA Aerospace Sciences Meetingand Exhibit, January 1992. AIAA paper 92-0439.[20] W. P. Jones and B. E. Launder. The prediction of laminarization with a two-equation model of turbulence. International Journal of Heat and Mass Transfer,15(2):301–314, 1972.145Bibliography[21] B. E. Launder, G. Jr. Reece, and W. Rodi. Progress in the development of aReynolds-stress turbulence closure. Journal of Fluid Mechanics, 68(3):537–566,1975.[22] B. Andersson, R. Andersson, L. Håkansson, M. Mortensen, R. Sudiyo, andB. van Wachem. Computational Fluid Dynamics for Engineers. CambridgeUniversity Press, 2011.[23] Pierre Sagaut and Sébastien Deck. Large eddy simulation for aerodynamics:status and perspectives. Philosophical Transactions of the Royal Society ofLondon A: Mathematical, Physical and Engineering Sciences, 367(1899):2849–2860, 2009.[24] Z.J. Wang, Krzysztof Fidkowski, Rémi Abgrall, Francesco Bassi, Doru Caraeni,Andrew Cary, Herman Deconinck, Ralf Hartmann, Koen Hillewaert, H.T.Huynh, Norbert Kroll, Georg May, Per-Olof Persson, Bram van Leer, andMiguel Visbal. High-order CFD methods: Current status and perspective.International Journal for Numerical Methods in Fluids, 72(8):811–845, 2013.[25] M. Terracol and S. Deck. Numerical investigation of the flow around a three-element high-lift airfoil using two zonal hybrid RANS/LES methods: ZDES andNLDE. In Progress in Hybrid RANS-LES Modelling, pages 345–355. Springer,2012.[26] Shia-Hui Peng, Bastian Nebenführ, and Lars Davidson. Lessons learned fromhybrid RANS-LES computations of a three-element airfoil flow. In Proceedingsof Twenty-First AIAA Computational Fluid Dynamics Conference, 2013. AIAAPaper 2013-2841.[27] Hui Zhu, Song Fu, Lei Shi, and Z. J. Wang. Implicit large-eddy simulation forthe high-order flux reconstruction method. AIAA Journal, pages 2721–2733,2016.[28] S. Lele. Compact finite difference schemes with spectral-like resolution. Journalof Computational Physics, 103:16–42, 1992.146Bibliography[29] M. Visbal and D. Gaitonde. On the use of higher-order finite-difference schemeson curvilinear and deforming meshes. Journal of Computational Physics,181(1):155–185, 2002.[30] Xiaogang Deng, Meiliang Mao, Guohua Tu, Hanxin Zhang, and Yifeng Zhang.High-order and high accurate CFD methods and their applications for complexgrid problems. Communications in Computational Physics, 11(4):1081–1102,2012.[31] D. A. Brown and D. W. Zingg. Performance of a Newton-Krylov-Schur algo-rithm for solving steady turbulent flows. AIAA Journal, pages 1–14, 2016.[32] M. M. Rai. A relaxation approach to patched grid calculations with the Eulerequations. Journal of Computational Physics, 66:99–131, 1986.[33] A. Lerat and Z. N. Wu. Stable conservative multidomain treatments for implicitEuler equations. Journal of Computational Physics, 123:45–64, 1996.[34] J. E. Hicken and D. W. Zingg. Parallel Newton-Krylov solver for the Eulerequations discretized using simultaneous-approximation terms. AIA Journal,46(11):2773–2786, 2008.[35] Ami Harten, Stanley Osher, Björn Engquist, and Sukumar R. Chakravarthy.Some results on uniformly high-order accurate essentially non-oscillatoryschemes. Applied Numerical Mathematics, 2:347–377, 1986.[36] Ami Harten, Björn Enquist, Stanley Osher, and Sukumar R. Chakravarthy.Uniformly high order accurate essentially non-oscillatory schemes, III. Journalof Computational Physics, 71(2):231–303, August 1987.[37] D. Zingg, S. De Rango, M. Nemec, and T. Pulliam. Comparison of several spa-tial discretizations for the Navier-Stokes equations. Journal of ComputationalPhysics, 160:683–704, 2000.[38] B. Cockburn, S. Hou, and Shu C. W. TVB Runge-Kutta local projectiondiscontinuous Galerkin finite element method for conservation laws IV: Themultidimensional case. Mathematics of Computation, 54(190):545–581, 1998.147Bibliography[39] B. Cockburn and C. W. Shu. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis,35(6):2440–2463, 1998.[40] F. Bassi and S. Rebay. High-order accurate discontinuous finite element solutionof the 2D Euler equations. Journal of computational physics, 138(2):251–285,1997.[41] Ralf Hartmann and Paul Houston. Adaptive discontinuous Galerkin finite ele-ment methods for the compressible Euler equations. Journal of ComputationalPhysics, 183(2):508–532, 2002.[42] F. Bassi and S. Rebay. A high-order accurate discontinuous finite elementmethod for the numerical solution of the compressible Navier-Stokes equations.Journal of Computational Physics, 131:267–279, 1997.[43] Krzysztof J. Fidkowski, Todd A. Oliver, James Lu, and David L. Darmo-fal. p-multigrid solution of high-order discontinuous Galerkin discretizations ofthe compressible Navier-Stokes equations. Journal of Computational Physics,207(1):92 – 113, 2005.[44] Nicholas K. Burgess and Dimitri J. Mavriplis. An hp-adaptive discontinuousGalerkin solver for aerodynamic flows on mixed-element meshes. In Proceedingsof the Forty-Ninth AIAA Aerospace Sciences Meeting, 2011.[45] N. K. Burgess and D. J. Mavriplis. High-order discontinuous galerkin methodsfor turbulent high-lift flows. In Proceeding of Seventh International Conferenceon Computational Fluid Dynamics (ICCFD7), 2012.[46] Timothy J. Barth and Paul O. Frederickson. Higher order solution of the Eulerequations on unstructured grids using quadratic reconstruction. AIAA paper90-0013, January 1990.[47] Amir Nejat. A Higher-Order Accurate Unstructured Finite Volume Newton-Krylov Algorithm for Inviscid Compressible Flows. PhD thesis, The Universityof British Columbia, Department of Mechanical Engineering, 2007.148Bibliography[48] Christopher Michalak. Efficient high-order accurate unstructured finite-volumealgorithms for viscous and inviscid compressible flows. PhD thesis, The Uni-versity of British Columbia, Department of Mechanical Engineering, 2009.[49] Christopher Michalak and Carl Ollivier-Gooch. Globalized matrix-explicitNewton-GMRES for the high-order accurate solution of the Euler equations.Computers and Fluids, 39:1156–1167, 2010.[50] Timothy J. Barth and Mats G. Larson. A posteriori error estimates for higherorder Godunov finite volume methods on unstructured meshes. Finite Volumesfor Complex Applications III, London, 2002.[51] Rémi Abgrall. On essentially non-oscillatory schemes on unstructured meshes:Analysis and implementation. Journal of Computational Physics, 114(1):45–58,1994.[52] Thomas Sonar. On the construction of essentially non-oscillatory finite volumeapproximations to hyperbolic conservation laws on general triangulations: poly-nomial recovery, accuracy and stencil selection. Computer Methods in AppliedMechanics and Engineering, 140:157–181, 1997.[53] Timothy J. Barth. Recent developments in high order k-exact reconstructionon unstructured meshes. AIAA paper 93-0668, January 1993.[54] Carl F. Ollivier-Gooch and Michael Van Altena. A high-order accurate unstruc-tured mesh finite-volume scheme for the advection-diffusion equation. Journalof Computational Physics, 181(2):729–752, sep 2002.[55] Christopher Michalak and Carl Ollivier-Gooch. Accuracy preserving limiter forthe high-order accurate solution of the Euler equations. Journal of Computa-tional Physics, 228(23):8693–8711, 2009.[56] Gene H. Golub and Charles F. Van Loan. Matrix Computations. The JohnsHopkins University Press, Baltimore, 1983.[57] J.L. Steger and R.F. Warming. Flux vector splitting of the inviscid gasdy-namics equations with applications to finite difference methods. Journal ofComputational Physics, 40(2):263–293, 1981.149Bibliography[58] Bram van Leer. Flux vector splitting for the Euler equations. ICASE ReportNo. 82-30, NASA Langley Research Center, 1982.[59] P. L. Roe. Approximate Riemann solvers, parameter vectors, and differenceschemes. Journal of Computational Physics, 43:357–372, 1981.[60] Meng-Sing Liou. A sequel to AUSM: AUSM+. Journal of ComputationalPhysics, 129(2):364–382, December 1996.[61] Meng-Sing Liou. A further development of the AUSM+ scheme towards robustand accurate solutions for all speeds. In Proceedings of the Sixteenth AIAAComputational Fluid Dynamics Conference. American Institute of Aeronauticsand Astronautics, June 2003. AIAA paper 2003-4116.[62] P. L. Roe. Characteristic-based schemes for the Euler equations. In AnnualReview of Fluid Mechanics, volume 18, pages 337–365. Annuals Reviews, Inc.,1986.[63] Charles Hirsch. Numerical computation of internal and external flows: Com-putational methods for inviscid and viscous flows. Wiley, 1990.[64] Andrew W. Cary, Andrew J. Dorgan, and Mori Mani. Towards accurate flowpredictions using unstructured meshes. In Proceedings of the Nineteenth AIAAComputational Fluid Dynamics Conference, 2009. AIAA paper 2009-3650.[65] A. H. Stroud and Don Secrest. Gaussian Quadrature Formulas. Prentice-Hall,Englewood Cliffs, N.J., 1966.[66] Carl Ollivier-Gooch, Amir Nejat, and Christopher Michalak. On obtaining andverifying high-order finite-volume solutions to the Euler equations on unstruc-tured meshes. American Institute of Aeronautics and Astronautics Journal,47(9):2105–2120, 2009.[67] Youcef Saad and Martin H. Schultz. GMRES: A generalized minimal residualalgorithm for solving nonsymmetric linear systems. SIAM Journal of Scientificand Statistical Computing, 7(3):856–869, July 1986.[68] D. A. Knoll, P. R. McHugh, and D. E. Keyes. Newton-Krylov methods for low-Mach-number compressible combustion. AIAA Journal, 34(5):961–967, 1996.150Bibliography[69] D. A. Knoll and D. E. Keyes. Jacobian-free Newton-Krylov methods: a surveyof approaches and applications. Journal of Computational Physics, 193(2):357–397, 2004.[70] Amir Nejat and Carl Ollivier-Gooch. Effect of discretization order on precon-ditioning and convergence of a high-order unstructured Newton-GMRES solverfor the Euler equations. Journal of Computational Physics, 227(4):2366–2386,2008.[71] David C. Wilcox. Comparison of two-equation turbulence models for boundarylayers with pressure gradient. AIA Journal, 31(8):1414–1421, 1993.[72] David C. Wilcox. Simulation of transition with a two-equation turbulencemodel. AIAA Journal, 32(2):247–255, 1994.[73] Xiaojuan Luo. An Automatic Adaptive Directional Variable p-version Methodin 3D Curved Domains. PhD thesis, Rensselaer Polytechnic Institute, 2005.[74] Spencer J. Sherwin and Joaquim Peiro. Mesh generation in curvilinear domainsusing high-order elements. International Journal for Numerical Methods inEngineering, 52(2-3):207–223, 2002.[75] Per Olaf Persson and Jaime Peraire. Curved mesh generation and mesh refine-ment using Lagrangian solid mechanics. In Proceedings of the Forty-SeventhAIAA Aerospace Sciences Meeting, 2009. AIAA 2009-949.[76] Li Wang, W. Kyle Anderson, J. Taylor Erwin, and Sagar Kapadia. Solutionsof high-order methods for three-dimensional compressible viscous flows. InProceedings of the Twentieth AIAA Computational Fluid Dynamics Conference,2012.[77] Dimitri J. Mavriplis. Revisiting the least-squares procedure for gradient re-construction on unstructured meshes. In Proceedings of the Sixteenth AIAAComputational Fluid Dynamics Conference, 2003.[78] Boris Diskin, James L. Thomas, Eric J. Nielsen, Hiroaki Nishikawa, and Jef-frey A. White. Comparison of node-centered and cell-centered unstructuredfinite-volume discretizations: Viscous fluxes. American Institute of Aeronau-tics and Astronautics Journal, 48(7):1326–1338, July 2010.151Bibliography[79] Natalia B. Petrovskaya. Discontinuous weighted least-squares approximationon irregular grids. Computer Modeling in Engineering and Sciences, 32:69–84,2008.[80] Natalia B. Petrovskaya. A quadratic least-squares solution reconstruction in aboundary layer region. International Journal for Numerical Methods in Biomed-ical Engineering, 26:1721–1735, 2010.[81] Natalia B. Petrovskaya. Data dependent weights in discontinuous weightedleast-squares approximation with anisotropic support. Calcolo, 48:127–143,2011.[82] Alireza Jalali and Carl F. Ollivier-Gooch. Higher-order finite volume solutionreconstruction on highly anisotropic meshes. In Proceedings of the Twenty-First AIAA Computational Fluid Dynamics Conference. American Institute ofAeronautics and Astronautics, 2013.[83] Lucian Ivan and Clinton Groth. High-order solution-adaptive central essentiallynon-oscillatory (CENO) method for viscous flows. Journal of ComputationalPhysics, 257:830–862, 2014.[84] H. Reichardt. Vollständige darstellung der turbulentengeschwindigkeitsverteilung in glatten leitungen. ZAMM - Journal of Ap-plied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematikund Mechanik, 31:208–219, 1951.[85] Christopher Rumsey. Turbulence Modeling Resource, NASA Langley ResearchCenter, http://turbmodels.larc.nasa.gov/, 2014.[86] Francesco Bassi, Andrea Crivellini, Stefano Rebay, and Marco Savini. Dis-continuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ωturbulence model equations. Computers and Fluids, 34(4-5):507 – 540, 2005.Residual Distribution Schemes, Discontinuous Galerkin Schemes and Adapta-tion.[87] Ngoc Nguyen, Per-Olof Persson, and Jaime Peraire. RANS solutions usinghigh order discontinuous Galerkin methods. In Proceedings of Forty-fifth AIAAAerospace Sciences Meeting and Exhibit, 2007. AIAA Paper 2007-914.152Bibliography[88] Todd A. Oliver. A High-Order, Adaptive, Discontinuous Galerkin Finite Ele-ment Method for the Reynolds-Averaged Navier-Stokes Equations. PhD thesis,Massachusetts Institute of Technology, 2008.[89] T. Oliver and D. Darmofal. An unsteady adaptation algorithm for discontin-uous Galerkin discretizations of the RANS equations. In Proceedings of theEighteenth AIAA CFD Conference, 2007. AIAA Paper 2007-3940.[90] M. Ceze and K. Fidkowski. Constrained pseudo-transient continuation. Inter-national Journal for Numerical Methods in Engineering, 102:1683–1703, 2015.[91] R. Hartmann. Higher-order and adaptive discontinuous Galerkin methods withshock-capturing applied to transonic turbulent delta wing flow. InternationalJournal for Numerical Methods in Fluids, 72:883–894, 2013.[92] Li Wang, W. Kyle Anderson, J. Taylor Erwin, and Sagar Kapadia. Discon-tinuous Galerkin and Petrov Galerkin methods for compressible viscous flows.Computers & Fluids, 100:13–29, 2014.[93] D. De Santis. High-order linear and non-linear residual distribution schemesfor turbulent compressible flows. Computer Methods in Applied Mechanics andEngineering, 285:1–31, 2015.[94] D. Moro, N. C. Nguyen, and J. Peraire. Navier-Stokes solution using hybridiz-able discontinuous Galerkin methods. In Proceedings of the Twentieth AIAAComputational fluid Dynamics Conference, 2011. AiAA Paper 2011-3407.[95] S. R. Allmaras, F. T. Johnson, and P. Spalart. Modifications and clarificationsfor the implementation of the Spalart-Allmaras turbulence model. In Proceed-ings of the Seventh International Conference on Computational Fluid Dynamics(ICCFD7), 2012.[96] M. Ceze and K. Fidkowski. Pseudo-transient continuation, solution updatemethods, and CFL strategies for DG discretizations of the RANS-SA equations.In Proceedings of the Twenty-First AIAA Fluid Dynamics Conference, 2013.AIAA paper 2013-2686.[97] F. Eleuterio Toro. Riemann Solvers and Numerical Methods for Fluid Dynam-ics. Springer-Verlag, New York, NY, 1999.153Bibliography[98] Boris Diskin, James L. Thomas, Eric J. Nielsen, Hiroaki Nishikawa, and Jef-frey A. White. Comparison of node-centered and cell-centered unstructuredfinite-volume discretizations. Part I: Viscous fluxes. In Proceedings of the Forty-Seventh AIAA Aerospace Sciences Meeting, 2009. AIAA paper 2009-597.[99] Hiroaki Nishikawa. Beyond interface gradient: A general principle for con-structing diffusion scheme. In Proceedings of the Fortieth AIAA Fluid DynamicsConference and Exhibit, 2010. AIAA paper 2010-5093.[100] Alireza Jalali, Mahkame Sharbatdar, and Carl Ollivier-Gooch. Accuracy anal-ysis of unstructured finite volume discretization schemes for diffusive fluxes.Computers and Fluids, 101:220–232, 2014.[101] R.C. Swanson and C.-C. Rossow. An efficient solver for the RANS equationsand a one-equation turbulence model. Computers & Fluids, 42:13–25, 2011.[102] Alireza Jalali and Carl Ollivier-Gooch. Higher-order unstructured finite volumeRANS solution of turbulent compressible flows. Computer & Fluids, 143:32–47,2017.[103] P.H. Cook, M.A. McDonald, and M.C.P. Firmin. Aerofoil RAE 2822 - pres-sure distributions, and boundary layer and wake measurements. In Experimen-tal Data Base for Computer Program Assessment, AGARD Report AR 138.AGARD, 1979.[104] First International Workshop on High-Order CFD Methods.http://people.ku.edu/ z651w035/hiocfd.html, 2012. Sponsored by theAIAA Fluid Dynamics Technical Committee, the US Air Force Office ofScientific Research, and the Deutsches Zentrum fur Luft- und Raumfahrt.[105] Third International Workshop on High-Order CFD Methods.https://www.grc.nasa.gov/hiocfd/. Sponsored by NASA, AIAA, DLR,and Army Research Office (ARO).[106] Alireza Jalali and Carl Ollivier-Gooch. Higher-order unstructured finite volumemethods for turbulent aerodynamic flows. In Proceedings of the Twenty-SecondAIAA Computational Fluid Dynamics Conference, number AIAA Paper 2015-2284, 2015.154Bibliography[107] Stefano Giani and Paul Houston. High–order hp–adaptive discontinuousGalerkin finite element methods for compressible fluid flows. In ADIGMA-AEuropean Initiative on the Development of Adaptive Higher-Order VariationalMethods for Aerospace Applications, pages 399–411. Springer, 2010.[108] Marco Ceze and Krzysztof J. Fidkowski. Anisotropic hp-adaptation frameworkfor functional prediction. AIAA Journal, 51(2):492–509, 2013.[109] Li Wang and Dimitri J Mavriplis. Adjoint-based hp adaptive discontinuousGalerkin methods for the 2d compressible Euler equations. Journal of Compu-tational Physics, 228(20):7643–7661, 2009.[110] Nicholas K. Burgess and Dimitri J. Mavriplis. hp-adaptive discontinuousGalerkin methods for the Navier-Stokes equations. AIAA Journal, 50(12):2682–2694, 2012.[111] Michael Woopen, Arvind Balan, Georg May, and Jochen Schütz. A comparisonof hybridized and standard DG methods for target-based hp-adaptive simula-tion of compressible flow. Computers & Fluids, 98:3–16, 2014.[112] W. Kyle Anderson, Behzad R. Ahrabi, and James C. Newman. Finite ele-ment solutions for turbulent flow over the NACA 0012 airfoil. AIA Journal,54(9):2688–2704, 2016.[113] Behzad R. Ahrabi and W. Kyle Anderson. An adjoint-based hp-adaptivePetrov–Galerkin method for turbulent flows. In Proceedings of Twenty-SecondAIAA Computational Fluid Dynamics Conference, 2015. AIAA Paper 2015-2603.[114] Lei Shi and Z. J. Wang. Adjoint-based error estimation and mesh adaptation forthe correction procedure via reconstruction method. Journal of ComputationalPhysics, 295:261–284, 2015.[115] Michael J. Aftosmis, Datta Gaitonde, and T. Sean Tavares. On the accuracy,stability, and monotonicity of various reconstruction algorithms for unstruc-tured meshes. AIAA paper 94-0415, January 1994.[116] M. J. Berger and P. Colella. Local adaptive mesh refinement for shock hydro-dynamics. Journal of Computational Physics, 82(1):64–84, 1989.155Bibliography[117] Mahkame Sharbatdar and Carl Ollivier-Gooch. Anisotropic mesh adaptation:Recovering quasi-structured meshes. In Proceedings of the Fifty-First AIAAAerospace Sciences Meeting, 2013.[118] Doug Pagnutti and Carl Ollivier-Gooch. A generalized framework for high orderanisotropic mesh adaptation. Computers and Structures, 87:670–679, 2009.[119] X. Zhang, M.-G. Vallet, J. Dompierre, P. Labbe, D. Pelletier, J.-Y. Trepanier,R. Camarero, J. Lassaline, L. Manzano, and D. Zingg. Mesh adaptation usingdifferent error indicators for the Euler equations. In Proceedings of the FifteenthAIAA Computational Fluid Dynamics Conference, 2001. AIAA Paper 2001-2549.[120] D. A. Venditti and D. L. Darmofal. Grid adaptation for functional outputs: Ap-plication to two-dimensional viscous flows. Journal of Computational Physics,187:22–46, 2003.[121] M. Ainsworth and J. T. Oden. A unified approach to a posteriori error esti-mation using element residual methods. Numerische Mathematik, 65(1):23–50,1993.[122] Haiyang Gao and Z.J. Wang. A residual-based procedure for hp-adaptation on2D hybrid meshes. In Proceedings of the Forty-ninth AIAA Aerospace SciencesMeeting, 2011. AIAA Paper 2011-492.[123] N. Ganesh, N. V. Shende, and N. Balakrishnan. A residual estimator basedadaptation strategy for compressible flows. In Computational Fluid Dynamics2006, pages 383–388. Springer, 2009.[124] Krzysztof J. Fidkowski and David L. Darmofal. Review of output-based er-ror estimation and mesh adaptation in computational fluid dynamics. AIAAJournal, 49(4):673–694, 2011.[125] Niles Pierce and Michael Giles. Adjoint and defect error bounding and correc-tion for functional estimates. Journal of Computational Physics, 200(2):769–794, 2004.[126] Michael Andrew Park. Adjoint-based, three-dimensional error prediction andgrid adaptation. AIAA Journal, 42(9):1854–1862, 2004.156Bibliography[127] D. A. Venditti and D. L. Darmofal. Adjoint error estimation and grid adapta-tion for functional outputs: Application to quasi-one-dimensional flow. Journalof Computational Physics, 164(1):204–227, October 2000.[128] D. A. Venditti and D. L. Darmofal. Grid adaptation for functional outputs: Ap-plication to two-dimensional inviscid flows. Journal of Computational Physics,175(1):40–69, February 2002.[129] Ralf Hartmann. Adaptive discontinuous Galerkin methods with shock-capturing for the compressible Navier–Stokes equations. International Journalfor Numerical Methods in Fluids, 51(9–10):1131–1156, 2006.[130] R. Hartmann., J. Held, and T. Leicht. Adjoint-based error estimation andadaptive mesh refinement for the RANS and k–ω turbulence model equations.Journal of Computational Physics, 230(11):4268–4284, 2011.[131] Mahkame Sharbatdar and Carl Ollivier-Gooch. Comparison of truncation errorestimators for defect correction and output error estimation on unstructuredmesh finite volume method. In Proceedings of the Fifty-Third AIAA AerospaceSciences Meeting, 2015.[132] Mahkame Sharbatdar, Alireza Jalali, and Carl Ollivier-Gooch. Smoothed trun-cation error in functional error estimation and correction using adjoint meth-ods in an unstructured finite volume solver. Computers & Fluids, 140:406–421,2016.[133] L. Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon, and J.E. Flaherty.Shock detection and limiting with discontinuous Galerkin methods for hyper-bolic conservation laws. Applied Numerical Mathematics, 48:323–338, 2004.[134] V. Venkatakrishnan. Convergence to steady-state solutions of the Euler equa-tions on unstructured grids with limiters. Journal of Computational Physics,118:120–130, 1995.[135] M. Woopen, G. May, and J. Schütz. Adjoint-based error estimation and meshadaptation for hybridized discontinuous Galerkin methods. International Jour-nal for Numerical Methods in Fluids, 76(11):811–834, 2014.157[136] D. Moxey, D. Ekelschot, U. Keskin, S. J. Sherwin, and J. Peiro. A thermo-elastic analogy for high-order curvilinear meshing with control of mesh validityand quality. Procedia Engineering, 82:127–135, 2014.[137] Nicholas K. Burgess. An Adaptive Discontinuous Galerkin Solver for Aerody-namic Flows. PhD thesis, University of Wyoming, December 2011.[138] Roland W. Lewis, Perumal Nithiarasu, and Kankanhalli N. Seetharamu. Fun-damentals of the Finite Element Method for Heat and Fluid Flow. Wiley, 2004.[139] Eugenio Oñate. Structural Analysis with the Finite Element Method. LinearStatics. Springer, 2013.158Appendix A: Cubic InterpolationFunctionsAs discussed in Chapter 3, cubic interpolation functions on 1D/2D reference elementsare required for flux integration over curved faces, curving the interior faces of ananisotropic mesh and source term integration over cubic cells. Here we present theinterpolation functions for each type.Figure A.1: Reference line segment for face integrationFor a reference line segment (Figure A.1), interpolation functions are:φ1 (ξ) =12(1− ξ) (1− 3ξ) (2− 3ξ)φ2 (ξ) =92(1− ξ) ξ (2− 3ξ)φ3 (ξ) =92(1− ξ) ξ (3ξ − 1)φ4 (ξ) =12ξ (3ξ − 1) (3ξ − 2)For a reference triangular element (Figure A.2), we use cubic Lagrangian inter-polation functions. To simplify the formulation of these functions, three dependentvariables are defined as [138]:L1 = 1− ξ − η , L2 = ξ , L3 = η159Appendix A: Cubic Interpolation FunctionsFigure A.2: Reference cubic triangular elementThe interpolation functions are given as:φ1 (ξ, η) =12L1 (3L1 − 1) (3L1 − 2)φ2 (ξ, η) =92L1L2 (3L1 − 1)φ3 (ξ, η) =92L1L2 (3L2 − 1)φ4 (ξ, η) =12L2 (3L2 − 1) (3L2 − 2)φ5 (ξ, η) =92L2L3 (3L2 − 1)φ6 (ξ, η) =92L2L3 (3L3 − 1)φ7(ξ,η) =12L3 (3L3 − 1) (3L3 − 2)φ8 (ξ, η) =92L3L1 (3L3 − 1)φ9 (ξ, η) =92L3L1 (3L1 − 1)φ10 (ξ, η) = 27L1L2L3For a reference quadrilateral element (Figure A.3), we use cubic serendipity in-160Appendix A: Cubic Interpolation FunctionsFigure A.3: Reference cubic quadrilateral elementterpolation functions. These functions are given as [139]:φ1 (ξ, η) =132(1− ξ) (1− η)(−10 + 9ξ2 + 9η2)φ2 (ξ, η) =932(1− 3ξ)(1− ξ2)(1− η)φ3 (ξ, η) =932(1 + 3ξ)(1− ξ2)(1− η)φ4 (ξ, η) =132(1 + ξ) (1− η)(−10 + 9ξ2 + 9η2)φ5 (ξ, η) =932(1 + ξ)(1− η2)(1− 3η)φ6 (ξ, η) =932(1 + ξ)(1− η2)(1 + 3η)φ7 (ξ, η) =132(1 + ξ) (1 + η)(−10 + 9ξ2 + 9η2)φ8 (ξ, η) =932(1 + 3ξ)(1− ξ2)(1 + η)φ9 (ξ, η) =932(1− 3ξ)(1− ξ2)(1 + η)φ10 (ξ, η) =132(1− ξ) (1 + η)(−10 + 9ξ2 + 9η2)161Appendix A: Cubic Interpolation Functionsφ11 (ξ, η) =932(1− ξ)(1− η2)(1 + 3η)φ12 (ξ, η) =932(1− ξ)(1− η2)(1− 3η)162Appendix B: Non-dimensionalEquationsIn most CFD solvers, a non-dimensional set of equations are discretized to keep themagnitude of flow quantities close to each other. In this Appendix, we describe thenon-dimensional parameters and groups used in our solver.To non-dimensionalize the system of RANS equations given in Chapter 4, we needto define reference variables to scale the flow field quantities. For this purpose, weset the speed of sound, density, temperature and viscosity at the free stream as thereference variables:ρref = ρ∞, uref = a∞, vref = a∞, Tref = T∞, µref = µ∞ (B.1)The reference variables for the other flow quantities can be obtained using the vari-ables of Equation B.1 and a characteristic length, L, as:Pref = ρ∞a2∞, ν˜ref =µ∞ρ∞, µT,ref = µ∞, tref =La∞, xref = L, yref = L (B.2)We also define our non-dimensional groups as:Ma∞ =V∞a∞, Re=ρV∞Lµ∞(B.3)A non-dimensional variable for a quantity of interest is defined as X∗ = X/Xref.However, we choose to show the non-dimensional variables without an asterisk tosimplify our notation. Using the reference values and non-dimensional groups, wecan re-write the conservation of mass, momentum and energy in the non-dimensionalform as:∂ρ∂t+∂ (ρu)∂x+∂ (ρv)∂y= 0 (B.4)163Appendix B: Non-dimensional Equations∂ (ρu)∂t+∂ (ρu2)∂x+∂ (ρuv)∂y= −∂P∂x+Ma∞Re(∂τxx∂x+∂τyx∂y)(B.5)∂ (ρv)∂t+∂ (ρuv)∂x+∂ (ρv2)∂y= −∂P∂y+Ma∞Re(∂τxy∂x+∂τyy∂y)(B.6)∂Et∂t+∂ (uEt)∂x+∂ (vEt)∂y= −(∂ (uP )∂x+∂ (vP )∂y)(B.7)+Ma∞Re(∂ (uτxx + vτxy)∂x+∂ (yτyx + vτyy)∂y+γγ − 1(µPr+µTPrT)(∂2T∂x¯2+∂2T∂y2))Similarly, the governing equation of the SA-negative turbulence model must benon-dimensionalized using the same reference variables:∂ (ρν˜)∂t+∂ (ρuν˜)∂x+∂ (ρvν˜)∂y=Ma∞Re · σ(∂∂x((µ+ µ′fnρν˜)∂ν˜∂x)+∂∂y((µ+ µ′fnρν˜)∂ν˜∂y)+ cb2 µ′ ρ(∂ν˜∂x)2+(∂ν˜∂y)2−(µρ+ µ′ν˜)(∂ρ∂x∂ν˜∂x+∂ρ∂y∂ν˜∂y))+ ρ (P −D) (B.8)The non-dimensional production and destruction terms are defined as:P =cb1 (1− ft2) S˜ν˜ ν˜ ≥ 0cb1 (1− ct3)Sν˜ ν˜ < 0 , D =Ma∞Reµ′(cw1fw − cb1κ2 ft2) (ν˜d)2ν˜ ≥ 0−Ma∞Reµ′cw1(ν˜d)2ν˜ < 0(B.9)whereS˜ =S + S¯ S¯ ≥ −cv2SS +S(c2v2S+cv3S¯)(cv3−2cv2)S−S¯S¯ ≤ −cv2S(B.10)164Appendix B: Non-dimensional EquationsS =∣∣∣∣∣∂u∂y − ∂v∂x∣∣∣∣∣ (B.11)S¯ =Ma∞Reµ′ν˜fv2κ2d2and the wall function contributing to the destruction term is obtained as:r = min(Ma∞Reµ′ν˜S˜κ2d2, 10)g = r + cw2(r6 − r)(B.12)fw = g[1 + c6w3g6 + c6w3]The other constants and functions are the same as the dimensional form of the SA-negative model.165


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items