Numerical Stability Analysesand Improvement ofCell-centered Finite VolumeMethods on UnstructuredMeshesbyReza ZangenehB.Sc., The University of Tehran, 2012M.Sc., The University of British Columbia, 2014A 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)September 2019© Reza Zangeneh, 2019The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled:Numerical stability analyses and improvement of cell-centered finite volume methods on unstructured meshessubmitted by Reza Zangeneh in partial fulfillment of the requirements forthe degree of Doctor of Philosophyin Mechanical EngineeringExamining Committee:Carl Ollivier-Gooch, Mechanical EngineeringSupervisor Mauricio Ponga de la Torre, Mechanical EngineeringSupervisory Committee Member Michael J. Ward, MathematicsUniversity ExaminerRajeev Jaiman, Mechanical EngineeringUniversity ExaminerJean-Pierre CroisilleExternal ExaminerAdditional Supervisory Committee Members:Chris Hill, ANSYSSupervisory Committee MemberMichael Friedlander, Computer ScienceSupervisory Committee MemberiiAbstractThe purpose of this thesis is to develop a framework in which one can detect and au-tomatically improve the numerical stability of cell-centered finite-volume calculations onunstructured meshes through optimization schemes that modify the mesh, the solutionreconstruction, or the boundary conditions. In this process, eigenanalysis and the gra-dients of the eigenvalues with respect to different parameters are used to ensure energystability of the system, consequently resulting in convergence.First, gradients of eigenvalues with respect to the local changes in the mesh are calcu-lated to find directions and magnitudes of mesh movements that will make the Jacobianof a semi-discrete system of equations negative semi-definite. These mesh movement vec-tors are used to modify the mesh locally. The numerical results show that the proposedmethods are able to locate the problematic parts of the mesh responsible for instabilitiesfor several physical problems and to correct the instability.Secondly, I develop a mathematical method, introduced by Haider et al., to mea-sure the stability impact of the reconstruction for high order and nonlinear problems,regardless of the solution. Second order and third order accurate advection and Burgersand Euler problems are used to present detailed practical results and discussion aroundthe use of the local reconstruction map for stability analysis. This method shows thatincreasing the stencil size will lead to more stable problems. An empirical study is per-formed which sheds light on connections between the mesh properties and the stabilityof the reconstruction. I also propose a systematic approach to optimize both the shapeand the size of the reconstruction stencil for better numerical stability through eigenvalueanalysis. In this approach, one can directly optimize the solution reconstruction stenciliiifor every control volume to obtain better numerical stability and convergence propertiesfor steady state problems.The rightmost eigenpairs of the spatially discretized system of equations are used toobtain an optimized boundary condition which will ensure the energy stability of thesystem. Lastly, the sensitivity of the rightmost eigenvalues to the solution is measuredto investigate the effect of using surrogate solutions for the purpose of linearizing thesemi-discretized Jacobian.ivLay SummaryComputer simulations are used extensively in scientific modeling or engineering design,where a mathematical equation to a physical problem is solved numerically. Significantcomputational cost and effort by a domain expert can be required to find a solution thatsatisfies a physical condition or a numerical threshold. Often, this process is iterativeand pursued by trial and error. In this work, a systematic approach is designed for thefirst time which can detect the unfavorable artifacts that may cause problems during thesimulation and automatically improve them to obtain a valid solution.vPrefaceThe research ideas and methods explored in this thesis are the fruits of a close workingrelationship between Dr. Carl Ollivier-Gooch and Reza Zangeneh. The implementationof the methods, the data analysis, and the manuscript preparation were done by theauthor, Reza Zangeneh, under the supervision of Dr. Carl Ollivier-Gooch throughoutthe process.Chapter 2 on modifying the mesh to obtain numerical stability has been published inComputers & Fluids (see [1]); preliminary versions of this work have also been presentedat the 2017 AIAA Aerospace Sciences Meeting, [2], and 2018 AIAA Aerospace SciencesMeeting, [3]. Chapter 3, on solution reconstruction optimization has been accepted forthe journal of Computational Physics (see [4]). Part of Chapter 3 has also been presentedat the 2017 AIAA Aerospace Sciences Meeting (see [5]). Finally, Chapters 4 and 5 havebeen submitted; this article is currently under review.• R. Zangeneh, C. Ollivier-Gooch, Mesh optimization to improve the stability offinite-volume methods on unstructured meshes, Computers & Fluids 156 (2017)590–601.• R. Zangeneh, C. F. Ollivier-Gooch, Stability analysis and improvement of the solu-tion reconstruction for cell-centered finite volume methods on unstructured meshes,Journal of Computational Physics in press (2019).• R. Zangeneh, C. F. Ollivier-Gooch, Boundary condition optimization to improvethe stability of finite-volume methods on unstructured meshes, submitted (2019).viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii1 Introduction and Background . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Flow Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.1.1 Higher-order Methods . . . . . . . . . . . . . . . . . . . . . . . . 71.1.2 Solution Reconstruction . . . . . . . . . . . . . . . . . . . . . . . 71.1.3 Flux Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.1.4 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.1.5 Analytical Flux Jacobian . . . . . . . . . . . . . . . . . . . . . . 121.1.6 Time Advance Schemes . . . . . . . . . . . . . . . . . . . . . . . 131.2 Numerical Stability Methods . . . . . . . . . . . . . . . . . . . . . . . . 17vii1.2.1 Von Neumann Stability Analysis . . . . . . . . . . . . . . . . . . 191.2.2 Matrix Stability Analysis: . . . . . . . . . . . . . . . . . . . . . . 211.2.3 Entropy Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . 241.2.4 L−Stable and A−Stable Methods . . . . . . . . . . . . . . . . . 261.3 Steady State Numerical Stability in this Thesis . . . . . . . . . . . . . . 271.3.1 Energy Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . 281.3.2 Eigenvalue Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 301.3.2.1 Eigenvalue Problem . . . . . . . . . . . . . . . . . . . . 301.3.2.2 Eigenvalue Calculation Algorithms . . . . . . . . . . . . 311.4 Surrogate Solution: Potential Flow . . . . . . . . . . . . . . . . . . . . . 321.5 Objectives and Novelty . . . . . . . . . . . . . . . . . . . . . . . . . . . 381.6 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 402 Mesh Optimization to Improve Stability . . . . . . . . . . . . . . . . . 422.1 Stability Improvement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 422.1.1 Eigensystems and Perturbations . . . . . . . . . . . . . . . . . . 432.1.2 Eigenvalue Derivatives With Respect To Mesh Movement . . . . 462.1.3 Optimization Scheme . . . . . . . . . . . . . . . . . . . . . . . . 492.1.4 Vertices to Move . . . . . . . . . . . . . . . . . . . . . . . . . . . 522.1.5 Allowing Boundary Movement . . . . . . . . . . . . . . . . . . . 542.2 A Priori Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 542.2.1 Lower Order Solution . . . . . . . . . . . . . . . . . . . . . . . . 582.2.1.1 Mesh Improvement . . . . . . . . . . . . . . . . . . . . 582.2.1.2 With Boundary Vertex Movement . . . . . . . . . . . . 632.2.2 Surrogate Solution: Panel Methods . . . . . . . . . . . . . . . . . 652.2.3 Surrogate Solution: One-equation Potential Flow . . . . . . . . . 652.2.4 Surrogate Solution: Two-equation Potential Flow . . . . . . . . . 722.3 A posteriori Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . 79viii2.3.1 Stability Analysis with Half Converged Solution . . . . . . . . . . 792.4 Computational Cost and Time Study . . . . . . . . . . . . . . . . . . . 843 Solution Reconstruction Optimization . . . . . . . . . . . . . . . . . . . 883.1 Local Reconstruction Map . . . . . . . . . . . . . . . . . . . . . . . . . . 903.1.1 1st Order Advection . . . . . . . . . . . . . . . . . . . . . . . . . 903.1.2 2nd Order Advection . . . . . . . . . . . . . . . . . . . . . . . . . 923.1.3 3rd Order Advection . . . . . . . . . . . . . . . . . . . . . . . . . 973.1.4 Inviscid Burgers’ Problem . . . . . . . . . . . . . . . . . . . . . . 1003.1.5 Euler and Navier-Stokes Problems . . . . . . . . . . . . . . . . . 1033.1.6 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 1053.2 Ordering the Stencil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1123.2.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 1143.3 Reconstruction Stencil Optimization by Eigenanalysis . . . . . . . . . . 1183.3.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 1253.3.1.1 Reconstruction Stencil Size . . . . . . . . . . . . . . . . 1263.3.1.2 Reconstruction Stencil Direction . . . . . . . . . . . . . 1324 Boundary Condition Optimization . . . . . . . . . . . . . . . . . . . . . 1394.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1414.1.1 Farfield Boundary Conditions . . . . . . . . . . . . . . . . . . . . 1414.1.1.1 Characteristics BC with vortex correction . . . . . . . . 1424.1.1.2 Pressure BC . . . . . . . . . . . . . . . . . . . . . . . . 1444.1.1.3 Dirichlet BC . . . . . . . . . . . . . . . . . . . . . . . . 1444.1.1.4 Radial Velocity BC . . . . . . . . . . . . . . . . . . . . 1444.1.1.5 Non-reflecting outflow BC . . . . . . . . . . . . . . . . . 1454.1.2 Wall Boundary Condition . . . . . . . . . . . . . . . . . . . . . . 1464.1.2.1 Rigid Wall . . . . . . . . . . . . . . . . . . . . . . . . . 146ix4.1.2.2 Characteristic Rigid Wall . . . . . . . . . . . . . . . . . 1464.1.2.3 Soft Wall . . . . . . . . . . . . . . . . . . . . . . . . . . 1474.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1484.2.1 Boundary Condition Optimization . . . . . . . . . . . . . . . . . 1484.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1504.3.1 Boundary Condition Optimization . . . . . . . . . . . . . . . . . 1504.3.1.1 Farfield . . . . . . . . . . . . . . . . . . . . . . . . . . . 1524.3.1.2 Wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1565 Eigenvalue and Solution Sensitivity . . . . . . . . . . . . . . . . . . . . . 1675.1 Sensitivity of Eigenvalues to Solution . . . . . . . . . . . . . . . . . . . . 1675.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1695.2.1 Sensitivity of Eigenvalues to Solution . . . . . . . . . . . . . . . . 1696 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1796.1 Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1796.2 Surrogate Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1806.3 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1806.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1816.5 Solution Sensitivity and Eigenanalysis . . . . . . . . . . . . . . . . . . . 1816.6 Recommendations for Future Work . . . . . . . . . . . . . . . . . . . . . 182Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198xList of Tables3.1 Local reconstruction map for a 2D, 2nd order advection . . . . . . . . . . 1063.2 Local reconstruction map for a 3D, 2nd order advection . . . . . . . . . . 1063.3 Local reconstruction map for a 3D, 3rd order advection . . . . . . . . . . 1103.4 2D, 2nd order inviscid Burgers’ problem with 1382 cells. . . . . . . . . . . 1123.5 Local reconstruction map for a 2D, 2nd order Euler . . . . . . . . . . . . 1153.6 Local reconstruction map for a 2D, 2nd order laminar . . . . . . . . . . . 1184.1 Wall boundary conditions and the right most eigenvalues . . . . . . . . . 1625.1 Jacobian and solution comparison of the high order and lower order schemeswhen optimizing the stencil size . . . . . . . . . . . . . . . . . . . . . . . 1775.2 Jacobian and solution comparison of the high order and lower order schemeswhen optimizing the stencil direction . . . . . . . . . . . . . . . . . . . . 1775.3 Jacobian and solution comparison of the steady state solution and the halfconverged solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1785.4 Jacobian and solution comparison of the high order and lower order schemeswhen optimizing the boundary condition types . . . . . . . . . . . . . . . 178xiList of Figures1.1 Vertex Centered vs Cell Centered control volumes . . . . . . . . . . . . . 41.2 Stability regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241.3 Panel method for a NACA 0012 airfoil . . . . . . . . . . . . . . . . . . . 352.1 Sensitivity map of the gradient of the Jacobian matrix with respect toperturbation size for an inviscid Burgers’ problem . . . . . . . . . . . . . 502.2 How to choose vertices for perturbation in finding the derivative of aneigenvalue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 532.3 Largest eigenvalue gradient components of the rightmost eigenvalue forNACA 0012, inviscid, Ma=0.8, α=2◦ . . . . . . . . . . . . . . . . . . . . 552.4 Boundary vertex movement . . . . . . . . . . . . . . . . . . . . . . . . . 562.5 Before and after the mesh movement for a 2nd order advection problem in a3D rectangular channel with five hundred control volumes. Time steppingis done by Crank–Nicolson. . . . . . . . . . . . . . . . . . . . . . . . . . . 592.6 Before and after mesh movement for an inviscid Burgers’ problem on arectangular channel with two thousand control volumes. Time stepping isdone by Crank–Nicolson. . . . . . . . . . . . . . . . . . . . . . . . . . . . 602.7 The snippet of mesh modification to stabilize the problem for an inviscidBurgers’ problem in a rectangular domain with two thousand control vol-umes. The solid lines are the original mesh and the dashed red lines arethe mesh after optimization. . . . . . . . . . . . . . . . . . . . . . . . . . 61xii2.8 Before and after mesh optimization for an Euler problem with Mach num-ber 0.6, and angle of attack of 5◦ around NACA 0012 airfoil with fivehundred control volumes. Time stepping is done by Crank–Nicolson. . . . 612.9 Mesh modification to stabilize the problem for an Euler problem withMach number 0.6, and angle of attack of 5◦ around NACA 0012 airfoilwith five hundred control volumes. The solid lines are the original meshand the dashed red lines are the mesh after optimization. . . . . . . . . . 622.10 Before and after mesh optimization for an Euler problem with Mach num-ber 0.6, and angle of attack of 5◦ around NACA 0012 airfoil with onehundred thousand control volumes. Time stepping is done by implicitEuler time advance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 632.11 Before and after mesh optimization for an Euler problem with Mach num-ber 0.4, and angle of attack of 1◦ around multi element airfoil with twentyfive thousand control volumes. Time stepping is done by implicit Eulertime advance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 642.12 Optimized meshes for NACA 0012, Inviscid, Ma=0.8, α=2◦ . . . . . . . . 662.13 Eigenvalue optimization for NACA 0012, Inviscid, Ma=0.8, α=2◦; everydot represents an iteration of the mesh optimization . . . . . . . . . . . . 672.14 Convergence rate for NACA 0012, Inviscid, Ma=0.8, α=2◦, Time steppingis done by backward nonlinear solver. . . . . . . . . . . . . . . . . . . . . 682.15 Comparison of second order Euler and Panel solutions for a NACA 0012airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 692.16 Residual history for both unstable and stabilized second order Euler problem 702.17 Original and optimized meshes for NACA 0012 airfoil with 10000 controlvolume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 712.18 Eigenvalues before and after mesh perturbation . . . . . . . . . . . . . . 72xiii2.19 Perturbed Mesh, NACA 0012, Inviscid, Mach number of 0.4, α = 5◦. Theone-equation potential model has been used in this case. . . . . . . . . . 732.20 Convergence history, NACA 0012, inviscid, Mach number of 0.4, α = 5degrees. The one-equation potential model has been used in this case. . . 732.21 The wake cut for NACA 0012 airfoil, Inviscid, Mach 0f 0.72, and angle ofattack of α = 5ř, and with 6200 control volumes . . . . . . . . . . . . . . 742.22 Original and optimized meshes for NACA 0012 airfoil, Inviscid, Mach 0.72,and α = 5 degrees, with 6200 control volumes. The two-equation potentialmodel has been used in this case. . . . . . . . . . . . . . . . . . . . . . . 762.23 L2 norm of residual for original and optimized meshes for NACA 0012airfoil, Inviscid, Mach 0.72, and α = 2ř with 6200 control volumes. Thetwo-equation potential model has been used in this case. . . . . . . . . . 762.24 Potential and Euler solution comparison . . . . . . . . . . . . . . . . . . 772.25 Residual history of Euler problem initialized with potential flow solution 782.26 Convergence history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 812.27 Mesh before and after stabilization of a third order Euler problem . . . . 822.28 Mesh before and after stabilization of a third order Laminar problem . . 822.29 Pressure coefficients around 0012 airfoil for two unstable and stabilizedproblems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 832.30 The computational time . . . . . . . . . . . . . . . . . . . . . . . . . . . 873.1 Norms of the local reconstruction map for second order 2D advection prob-lem with different stencil sizes on a two dimensional channel with 1500control volumes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1073.2 Norms of the local reconstruction map for second order 3D advection prob-lem with different stencil sizes on a three dimensional channel with 650control volumes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1083.3 Local reconstruction map norms for a 3D, 3rd order advection problem . 109xiv3.4 First 25 eigenvalues for a 3D mesh 3rd order advection problem for 2 dif-ferent stencil sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1103.5 2D, 2nd order inviscid Burgers’ problem on an unstructured mesh of a twodimensional channel with with 1382 cells. The highlighted control volumeshave a larger stencil size of 9. . . . . . . . . . . . . . . . . . . . . . . . . 1133.6 Convergence history for transonic Euler problem with Mach 0.8 and angleof attack of 2◦ around a NACA 0012 airfoil. . . . . . . . . . . . . . . . . 1153.7 Convergence history for subsonic Euler problem with Mach 0.5 and angleof attack of 2◦ around a NACA 0012 airfoil. . . . . . . . . . . . . . . . . 1163.8 Local reconstruction map norms for a 2D, 2nd order Euler problem . . . . 1173.9 Local reconstruction map norms for a 2D, 2nd order laminar problem . . 1193.10 Determining the upwind stencil for a control volume, α, where β is acontrol volume in the stencil of the α. −→U is the local velocity at controlvolume α. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1213.11 Upwind stencil versus regular stencil for three control volume samples inan unstructured mesh with 1500 control volumes around a NACA 0012airfoil for an Euler scheme with M = 0.72, α = 3◦. . . . . . . . . . . . . . 1213.12 Local Jacobian and reconstruction calculation . . . . . . . . . . . . . . . 1253.13 Second order transonic Euler stencil size optimization . . . . . . . . . . . 1273.14 Second order supersonic Euler stencil size optimization . . . . . . . . . . 1283.15 Residual history before and after stencil size optimization for second orderEuler problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1293.16 Third order subsonic laminar stencil size optimization . . . . . . . . . . . 1303.17 A second order supersonic Euler wind tunnel with ramp for M = 2 and3500 control volumes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131xv3.18 Stencil size optimization; Second order supersonic Euler problem on an un-structured mesh with 3500 control volumes and with inflow Mach numberof 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1323.19 Second order transonic Euler stencil direction optimization . . . . . . . . 1333.20 Second order supersonic Euler stencil direction optimization . . . . . . . 1343.21 Residual history before and after stencil direction optimization for secondorder Euler problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1353.22 Second order subsonic laminar stencil direction optimization . . . . . . . 1373.23 Residual history before and after stencil direction optimization for secondorder laminar problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1384.1 Characteristic lines for a subsonic two-dimensional flow at an outlet forflow in the x-direction. The y direction is normal to the page. . . . . . . 1424.2 Different types of eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . 1514.3 Derivative of the eigenvalue with vertex movement for different eigenmodetypes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1534.4 Boundary condition configurations: far-field boundary condition optimiza-tion for an Euler problem with Mach 0.55, and an angle of attack of 2◦ ona NACA 0012 airfoil with 2000 control volumes . . . . . . . . . . . . . . 1544.5 History of convergence: far-field boundary condition optimization for anEuler problem with Mach 0.55, and an angle of attack of 2◦ on a NACA0012 airfoil with 2000 control volumes . . . . . . . . . . . . . . . . . . . . 1554.6 Boundary condition optimization between Rigid Walls (RW), and Char-acteristic Walls (CW), for two different flow conditions around a NACA0012 airfoil with 5000 control volumes. The leading edge of the airfoil isdisplayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156xvi4.7 Boundary condition optimization between Rigid Walls (RW), and Char-acteristic Walls (CW), for two different flow conditions around a NACA0012 airfoil with 5000 control volumes. The central part of the airfoil isdisplayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1574.8 Boundary condition optimization between Rigid Walls (RW), and Char-acteristic Walls (CW), for two different flow conditions around a NACA0012 airfoil with 5000 control volumes. The trailing edge of the airfoil isdisplayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1574.9 Boundary condition optimization between Rigid Walls (RW), and SoftWalls (SW), for two different flow conditions around a NACA 0012 airfoilwith 5000 control volumes. The leading edge of the airfoil is displayed. . 1584.10 Boundary condition optimization between Rigid Walls (RW), and SoftWalls (SW), for two different flow conditions around a NACA 0012 airfoilwith 5000 control volumes. The central part of the airfoil is displayed. . . 1594.11 Boundary condition optimization between Rigid Walls (RW), and SoftWalls (SW), for two different flow conditions around a NACA 0012 airfoilwith 5000 control volumes. The trailing edge of the airfoil is displayed. . 1594.12 Boundary condition optimization between Characteristic Walls (CW), andSoft Walls (SW), for two different flow conditions around a NACA 0012airfoil with 5000 control volumes. The leading edge of the airfoil is displayed.1604.13 Boundary condition optimization between Characteristic Walls (CW), andSoft Walls (SW), for two different flow conditions around a NACA 0012airfoil with 5000 control volumes. The central part of the airfoil is displayed.1604.14 Boundary condition optimization between Characteristic Walls (CW), andSoft Walls (SW), for two different flow conditions around a NACA 0012airfoil with 5000 control volumes. The trailing edge of the airfoil is displayed.161xvii4.15 L2 norm of flux residual for different wall boundary condition types, M =0.81, and α = 1.25◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1624.16 L2 norm of flux residual for different wall boundary condition types, M =0.88, and α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1634.17 Pressure coefficient along the NACA 0012 airfoil with M = 0.81, andα = 1.25◦with different wall boundary conditions . . . . . . . . . . . . . . 1644.18 Pressure coefficient along the NACA 0012 airfoil with M = 0.88, andα = 0◦with different wall boundary conditions. In this case, Rigid walland Characteristic wall boundary conditions did not converge, hence theyare absent in this analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . 1654.19 Pressure distribution of Euler problem around NACA 0012 airfoil with5000 control volumes and with M = 0.81, and α = 1.25◦ . . . . . . . . . 1665.1 The residual history for an Euler problem with Mach 0.8 and an angle ofattack of α = 1.25◦ on NACA 0012 airfoil with 4200 control volumes. . . 1705.2 Sensitivity of the rightmost eigenvalue with respect to the density, ρ, atthree different stages of convergence (see Figure 5.1). . . . . . . . . . . . 1715.3 Sensitivity of the rightmost eigenvalue with respect to the axial velocity,u, at three different stages of convergence (see Figure 5.1). . . . . . . . . 1725.4 Sensitivity of the rightmost eigenvalue with respect to the vertical velocity,v, at three different stages of convergence (see Figure 5.1). . . . . . . . . 1735.5 Sensitivity of the rightmost eigenvalue with respect to the pressure, p, atthree different stages of convergence (see Figure 5.1). . . . . . . . . . . . 174xviiiNomenclatureθα the second order gradient coefficients of the control volume α in its solution re-constructioncp Specific heat at constant pressurec∞ free stream sound speedCl the lift coefficienteˆr the radial unit vectorEn amplitud of nth error harmonicEt Total energyFa addvective fluxFd diffusive flux~F flux vectorG Amplification factorI identity matrixJ˜ The perturbed Jacobian matrixJ The Jacobian matrixJ global Jacobian of residualxixJ (1) The perturbation matrix in general or derivative Jacobian matrix specificallyJ∗ Conjugate transpose of JJT Transpose of JL Left eigenvectorL matrix of left eigenvectorsM∞ free stream Mach numbernαβ the unit normal vector between control volume α and βnˆ unit normal vectorP Pressure, ProductionPtot total pressurePr Prandtl number< flux integral or residualR˜ The right eigenvector of the perturbed Jacobian matrixR The right eigenvectorR matrix of right eigenvectorsS source termT Temperaturet timeTtot total temperatureu x-velocityxxUR reconstructed solution~U discretized solution vector~u solution vectorv y-velocityx, y Cartesian coordinates~x Cartesian vector of locationα the control volume of interestβ face neighbor of control volume αδU change in the solution vector or solution update vector The infinitesimal perturbationγ˙ Shear rate tensorγ Ratio of specific heatsγ control volumes in the reconstruction stencilLαγ distance between the cell reference of the control volume α and γκ jump term coefficientλ Eigenvalueλ∗ Eigenvalue of conjugate transpose of Jλ˜ The eigenvalue of the perturbed Jacobian matrixµ Dynamic viscosityΩ size of the control volumexxi∂Ω control volume peripheralϕ the potential variable<{.} Real part of the complex numberρ Densityρ∞ free stream densityρ˜ upwind isentropic densityΨ the adjoint variableσα the first order gradient coefficients of the control volume α in its solution recon-structionτ Shear stress tensor4tα local time step for control volume α% pseudo-spectral radiusξ The mesh vertex perturbation(.)+ the positive part of an inner productΓαβ the boundary area between control volume α and βxxiiAcknowledgmentsFirst and foremost, I would like to express my sincere appreciation to my wonderful fam-ily, especially my mother, Halimeh, and my father, Hossein, for their love and encourage-ment throughout my life. I would also like to acknowledge my research supervisor, Dr.Carl Ollivier-Gooch, whose constant support, patience, and guidance has been invaluablethroughout the course of this research as well as the research committee, especially Dr.Chris Hill. Above all, I am grateful beyond what words can describe for my wonderfulwife and partner in life, Kimia, who helped me with proofreading, endured this processwith me, encouraged me, and always believed in me.I would also like to thank the financial support that made this research possible fromNatural Sciences and Engineering Research Council of Canada, the Four Year Fellowshipfrom the University of British Columbia and ANSYS Canada.xxiiiChapter 1Introduction and BackgroundMany physical phenomena in fluid flows are mathematically modeled with partial dif-ferential equations (PDE). These PDEs describe properties of the flow such as velocity,pressure, density and temperature at any location in space and at any point in time.Thus, such PDEs, which are called the governing equations of specific physical problems,are of paramount importance in the field of fluid dynamics. Only rarely do these PDEshave an analytical solution. The use of Computational Fluid Dynamics (CFD) is a robustalternative to solve these PDEs for practical and complex fluid flows. These numericalmethods first decompose the computational domain, which is known as mesh genera-tion. Secondly, they numerically model the continuous PDE into a finite representationof flow on the elements of the mesh. This is known as discretization which is normallyimplemented by finite difference, finite element, or finite volume methods.Mesh generation is broadly categorized into two different types: structured meshesand unstructured meshes. In a structured mesh, all the vertices and mesh elementshave a fixed and predefined topology. This means that a vertex can traverse to nearbyvertices in a structured manner through specific vertex or cell indexing. For example,cell (i+ 1, j) is always to the right of the cell (i, j) in a two dimensional space. This fixedtopology makes the flow solver, which numerically solves the discretized PDE, muchfaster and easier to implement. However, this fixed topology renders the automationof mesh generation extremely difficult for general problems with complex domains. Onthe other hand, unstructured meshes do not require a fixed topology and are thereforemore suitable for automatic mesh generation. Since there is no fixed topology to the1unstructured meshes, the element and the mesh connectivity information have to beprovided to the flow solver.As the capabilities of CFD tools have grown, so too have the size and complexityof the problems to which industrial users apply CFD. Even for expert users — whounderstand how to generate meshes and choose flow solver options to get better solutionsfor routine problems in their area of expertise — new large scale problems are challenging.Trial and error is often required to identify and resolve important flow features andfind a stable solution. Historically the main tool for increasing the solution accuracyhas been either employing grid refinement or high order schemes. However, for realworld application problems, the baseline simulation often pushes the limits of availablecomputing resources. Such studies are often prohibitive due to instability issues. Thisproblem is particularly challenging since commercial CFD tools typically handle complexproblem geometries using unstructured meshes, for which accuracy and stability issuesare not as well understood as for structured meshes.Therefore, with this rapid development of high order numerical methods comes theneed for stability analysis. However, these studies are not yet sufficient to fully under-stand and predict unstable features of finite volume methods on general unstructuredmeshes; hence the need to remedy them. The stability of numerical discretization meth-ods depends not only on the methodology but also on the mesh: bad features in a meshfar from any flow features of interest can still have a deleterious effect on the stability orconvergence of a CFD solver. A thorough understanding of this connection can provideguidance in the design of numerical methods or mesh generation that would improvesolver performance and robustness. Energy stability analysis has often been used in thefinite element (e.g., [6]), finite volume (e.g., [7, 8]), and finite differences (e.g., [9]) com-munities to find some upper bound and thresholds of stability. These methods hingemathematically on eigenanalysis. Entropy stability methods (e.g. see [10, 11, 12] andthe references therein) utilize entropy variables to devise schemes for nonlinear partial2differential equations for which the mathematical entropy is non-decreasing and whichare therefore stable. Energy and entropy stability methods are very powerful for analysisof stability, but to the best of our knowledge have never been applied to modify mesh ordiscretization scheme features as a controlling feedback tool to stabilize an unstable case.This work extends energy analysis to study and improve the mathematical stability of thesemi-discrete system of equations arising from unstructured mesh space discretizations.In the end, I provide a practical tool to predict instability and automatically stabilize adiscretized linear or nonlinear PDE on a general unstructured mesh.This thesis is built upon a higher order unstructured finite volume framework (ANSLib),which has been previously developed in the Advanced Numerical Simulation Laboratory.In finite volume methods, the volume integrals containing divergence terms are trans-formed to surface integrals using the divergence theorem. The net change of the propertyflux across the surfaces of the control volume plus any generation of that property re-flects how the property changes with time. Moreover, the flux leaving a control volumeis equal to the flux entering the adjacent control volume. Thus, finite volume methodsdiscretize and solve partial differential equations (PDEs) in a fully conserved form. Thisis a desirable property of finite volume methods as they always satisfy the conservation ofmass, momentum and energy, which is not easily obtained exactly by other discretizationschemes when applied to complex flows.Generally there are two different ways that control volumes for a tesselated physicaldomain can be formulated (see Figure 1.1): vertex centered, or cell centered. In the formercase, the control volumes’ reference points are the vertex locations and the boundaries ofthe control volumes are formed by connecting the midpoint of edges incident on the vertex(median-dual). On the other hand, the cell-centered control volumes are formed by thecells where the edges of the control volume are the edges of the tesselation and the controlvolume reference points are located at the centroid of the cells. ANSLib supports bothvertex centered and cell centered control volumes for two dimensional flows. However, it3Figure 1.1: Vertex Centered vs Cell Centered control volumesonly supports cell centered control volumes for three dimensional flows. In this thesis, thefocus is only on cell centered finite volume methods. In this chapter, I further describethe spatial discretization, flow scheme, boundary condition implementation and timeadvanced schemes in ANSLib as well as other supporting materials for the rest of thethesis.The rest of this chapter is organized as: first, in Section 1.1 I briefly describe thephysical flows and governing equations such as Euler and Navier-Stokes. Then, thenumerical schemes used to discretize these equations in both space and time are explainedwhich include the solution reconstruction, analytical flux Jacobian, and time advanceschemes. Secondly, in Section 1.2, the general and common numerical stability methodsare briefly explained. In Section 1.3, I define the stability method that I am using in thisthesis. In the next chapters, I will use surrogate solutions for the purpose of investigatingthe numerical stability of finite volume discretization of Euler schemes in a priori. Onesuch surrogate solution is the potential solution; which I briefly explain the differentmethods of calculating in Section 1.4. Section 1.5 points out the specific objectives andnovelties of this work, while Section 1.6 describes a general outline of the entire thesis.41.1 Flow SchemeThe partial differential equations solved in compressible CFD can be reduced to thevector equation:∂u(−→x , t)∂t+∇ · F (−→x , u, t) = S(−→x , t) (1.1)In Eq. 1.1, u(−→x , t) is the exact solution at the location −→x and time t. F is the fluxfunction, which is composed of the convective flux and the viscous flux terms. S is thesource term. For example, for a two dimensional Newtonian compressible laminar flow,the solution and the flux vectors are:U =ρρuρvEt, (1.2)F x =ρuρu2 + P − τxxρuv − τxyu (Et + P )−(uτxx + vτxy + cp(µPr)∂T∂x), (1.3)F y =ρvρuv − τyxρv2 + P − τyyv (Et + P )−(uτyx + vτyy + cp(µPr)∂T∂y)(1.4)where pressure, P , and the stress tensor, τ , are calculated using:5τij = 2µγ˙ij (1.5)γ˙ij =12(∂u∂x+ ∂v∂y)− 13(∇ ·−→u )δij (1.6)P = (γ − 1)[Et − 12ρ(u2 + v2)](1.7)where ρ is the fluid density, −→u (−→x , t) = (u, v) are the velocity components in the Cartesiansystem, Et is the total energy, cp is the specific heat at constant pressure, T is the fluidtemperature, Pr is the Prandtl number and µ is the fluid dynamic viscosity. A gasproperty equation such as the ideal gas law is used as the last equation to close the setof equations of Eq. 1.2-1.7. Equations and parameters above are non-dimensionalized bythe free stream flow properties.Integrating Eq. 1.1 over a control volume and applying the divergence theorem willresult in:dUidt= 1|Ωi|ˆΩidudtdV = − 1|Ωi|˛∂ΩiFdA+ 1|Ωi|ˆΩiSdV = <(Ui) (1.8)where Ωi denotes the size of the ith control volume and Ui is the control volume, CV,average of the solution. The flux integral is calculated by numerical quadrature rules.Because the compressible Navier-Stokes equations are homogeneous of degree one, onecan write the spatially discretized Eq. 1.8 as:dUdt= ∂<∂UU = JU (1.9)where < is the flux integral vector, or residual function, and J = ∂<∂Uis the global Jacobianof the residual.61.1.1 Higher-order MethodsIn this thesis we investigate the numerical stability of both the second order and higherorder finite volume methods. Methods that have an order of accuracy more than two arecalled higher order methods. With higher order discretization, the errors in satisfyingthe partial differential equations decrease faster than their lower order counter-parts, andgenerally reach a given accuracy threshold in less time and memory (e.g., see [13, 14]).However, this reduction in overall time comes with an increased hardship in discretiza-tion and sometimes stability. The early works on developing finite volume methods onunstructured meshes goes back to the quadratic reconstruction methods developed byBarth and Frederickson [15] for the Euler equations. Since then, many other researchers(e.g., see [16, 17, 13] and the references therein) have developed higher order finite volumemethods for compressible flows and have shown the merits of using these higher ordermethods, as opposed to their lower order counterparts, for convergence to the solutionof the PDE.1.1.2 Solution ReconstructionThe spatial accuracy of the finite volume method depends on the accuracy of the numeri-cal fluxes and the accuracy of flux integrations. Therefore, for a high order discretization,there is a need for a more accurate approximation to the solution and its gradients atGauss integration points and within the control volumes rather than the control volumeaverage Uα. To obtain this, the solution inside each control volume is approximated bypiecewise polynomials. These polynomials are obtained by expanding the solution in aTaylor series within each control volume, shown here in two dimensions:7URα (x, y) = Uα +∂U∂x∣∣∣∣∣α(x− xα) + ∂U∂y∣∣∣∣∣α(y − yα) + ∂2U∂x2∣∣∣∣∣α(x− xα)22 +∂2U∂x∂y∣∣∣∣∣α(x− xα) (y − yα) + ∂2U∂y2∣∣∣∣∣α(y − yα)22 + · · · (1.10)where URα (x, y) is the reconstructed solution at location (x, y) while the coefficients ofthe polynomial, Uα and ∂n+mUα∂xn∂ymare the solution and its derivatives at the cell referencepoint for control volume α. The degree of the polynomial determines the order of theaccuracy of the solution reconstruction.In compressible flow simulations, the solution discretization leads to discontinuitiessuch as shock waves. To ensure the desired order of accuracy, essentially non-oscillatory(ENO) or weighed essentially non-oscillated (WENO) reconstruction schemes were intro-duced (see [18, 19, 20] for more details). These schemes select a set of control volumesin the stencil of the reconstructions that satisfy a smoothness criterion. In practice, thismeans that these schemes choose control volumes for the reconstruction only from oneside of a discontinuity. However, there are some difficulties and complexities for selectingthe best stencil for each of control volumes in the flow domain. In particular the sten-cil can change from one time step to the next, making the steady-state condition verydifficult to achieve.Another approach for reconstruction, which is the choice of reconstruction in thisthesis and in ANSLib, is the the k-exact solution reconstruction. This method, whichwas first implemented in a finite volume upwind scheme by Barth [21], was later extendedby Ollivier-Gooch and Van Altena [22] to describe a new k-exact solution reconstructionmethod for the higher-order approximation of solution and its derivatives within a controlvolume. These k-exact reconstruction methods are used to obtain high-order accuracyfor finite volume methods (e.g., see [22, 23, 24]).In k-exact solution reconstruction, the solution in each control volume in the re-8construction stencil of cell α is approximated by Eq. (1.10) (using the same solutiongradients, ∂U∂x∣∣∣α, ∂2U∂x2∣∣∣α, . . .) and averaged over each CVγ in the stencil. In this way, onecan approximate the mean value of the solution in CVγ by the reconstruction functionin CVα (Eq. 1.10):1|Ωγ|¨ΩγURα (~x− ~xα) dV = Uα +∂U∂x∣∣∣∣∣αx̂αγ + · · ·+ ∂n+mU∂xn∂ym∣∣∣∣∣αx̂nymαγ (1.11)where x̂nymαγ terms are the relative moments of the CVγ from CVα calculated as sug-gested by Ollivier-Gooch and Van Altena [22]:x̂nymαγ =m∑l=0n∑k=0m!l!(m−l)!n!k!(n−k)! (xγ − xα)k (yγ − yα)l xn−kym−lγ (1.12)and xn−kym−lγ terms are the local moments of the CVγ :xn−kym−lγ =1|Ωγ|¨Ωγ(x− xγ)n−k (y − yγ)m−l dV (1.13)A collection of Eq. 1.11 over all the nearby control volumes in the stencil of thecontrol volume α forms the reconstruction matrix equation (Eq. (1.14)) from which thereconstruction coefficients ( ∂U∂x∣∣∣α, ∂2U∂x2∣∣∣α, . . .) can be calculated by least-squares.1 x̂α ŷα · · · x̂nymαα · · ·ωαγ1 ωαγ1x̂αγ1 ωαγ1 ŷαγ1 · · · ωαγ1x̂nymαγ1...... ... ... . . . ... ...ωαγN ωαγN x̂αγN ωαγN ŷαγN · · · ωαγN x̂nymαγN · · ·︸ ︷︷ ︸AU∂U∂x∂U∂y...∂n+mU∂xn∂ym...︸ ︷︷ ︸σ=UαUγ1...UγN(1.14)9The first equation is the conservation of the mean constraint for control volume α whichensures that the reconstruction function calculates the control volume averaged solutionin the CVα exactly. Gauss elimination is applied to remove this constraint from theremaining equations, which are solved in a least-squares sense. The weights, ω, areinversely related to the geometrical distances of the control volume and its neighbors togive more emphasis to the ones closer:ωαγi =1Lαγ =1‖−→x α −−→x γ‖ (1.15)where Lαγ is the distance between the cell reference of the control volume α and γ. See[22] for more details.1.1.3 Flux FunctionsAfter the solution reconstruction is calculated for every control volume, and since con-tinuity of the solution across the control volume boundary is not imposed, the solutionsand the gradient of solutions at both sides of every arbitrary point at the control volumeboundaries are different. If the flux, which is leaving a control volume and entering itsneighbor, is calculated twice using data from the two sides of the control volume bound-ary, they will be different. However, due to the conservation property these fluxes mustbe equal. This in turn means that special flux functions have to be computed usingthe solutions at both sides of a control volume surface to unify the flux values. In thisthesis, I use Roe’s flux splitting scheme [25] to calculate the convective flux (Eq. 1.16).The diffusive flux is computed by averaging the flux at both sides of the control volumeboundary and adding a jump term (Eq. 1.17):10Fa =F (Uα) + F (Uβ)2 +12∣∣∣A˜∣∣∣ (Uα − Uβ) (1.16)Fd =F (Uα) + F (Uβ)2 + κUα − Uβ∣∣∣−→L αβ · nˆ∣∣∣ (1.17)with∣∣∣A˜∣∣∣ = Γ−1 |Λ|Γ where Γ is is the matrix of the eigenvectors of the flux Jacobian,∂F∂U, and |Λ| is a diagonal matrix with the absolute values of the eigenvalues of the fluxJacobian as the diagonal entries. In Eq. 1.17, −→L αβ is the vector connecting the referencepoints of the two adjacent control volumes, nˆ is the normal vector of the control volumeboundary, and κ is the jump coefficient. It has been shown (see [26] for more details),that the jump term improves convergence when the diffusive fluxes are dominant.1.1.4 IntegrationOnce numerical fluxes are evaluated at control volume boundaries, they will be integratedaround the control volumes. Moreover, source terms and moments of area are integratedover each control volume. The order of accuracy of these integrations have to be equalto, or higher than the order of accuracy of the solution reconstruction; otherwise thetotal accuracy of the solution will be compromised. Gauss quadrature integration [27]evaluates these integrals over a specific set of points to satisfy a given order of accuracy.For example, for a second order accuracy, evaluating the integral at only one quadraturepoint at the mid-point of the control volume surface is required. A comprehensive list ofrequired points and their weights can be found in [27].The flux integral for each control volume is calculated by iterating over the controlvolume faces only once, and adding the value of the flux to the control volume into whichthe face normal vector points, and subtracting the flux from the adjacent control volume.111.1.5 Analytical Flux JacobianAdditionally, I calculate the global Jacobian matrix, which is the derivative of the fluxintegral with respect to the solution variables. In other words, it shows the sensitivityof the flux integrals at the control volume boundaries with respect to the changes of thesolution. Throughout this thesis the Jacobian matrix is extensively analyzed for stabilityand convergence properties. The global Jacobian can be calculated indirectly by Frechetderivatives or directly by applying chain differentiation to the flux integral. In this paper,I use the second approach, introduced by [17], namely:Jn×n =∂<∂U= ∂FluxInt∂CVar =∂F luxInt∂F lux︸ ︷︷ ︸Qn×g∂F lux∂RSol︸ ︷︷ ︸Fg×g∂RSol∂RCoef︸ ︷︷ ︸Rg×c∂RCoef∂PV ar︸ ︷︷ ︸Cc×n∂PV ar∂CV ar︸ ︷︷ ︸Vn×n(1.18)where n is the number of degrees of freedom (number of control volumes multiplied bythe number of unknowns) in Eq. 1.8, g is the overall number of integration points and cis the number of reconstruction coefficients. In Eq. 1.18:• First, for flow problems, the conserved variables (such as density, momentum, andenergy) for each control volume are transformed to the primitive variables (such asdensity, velocity, pressure) that are reconstructed within each control volume. Inour solver, we can apply the solution reconstruction scheme on either the primitivevariables or the conserved variables. However, we have observed that applying thesolution reconstruction on the primitive variables results in better convergence insome cases. Dirichlet boundary conditions such as constant pressure or velocity canalso be applied as constraints on the solution reconstruction, Eq. (1.14). Addition-ally, it is more physical to use solution limiters on pressure instead of energy. Forthese reasons, we have opted to reconstruct the primitive solution variables suchas pressure, velocity and density at the Gauss quadrature points. However, we use12conserved forms of numerical fluxes for calculating the flux integral at control vol-umes. Thus, the Jacobian V in Eq. (1.18) is used to transform the reconstructedprimitive solution at the Gauss quadrature points to the reconstructed conservedsolution. The Jacobian V of the transformation is calculated by ∂PVar∂CVar .• Second, for both the adjacent control volumes at each Gauss quadrature point:– The derivative of the reconstruction coefficients — the derivatives in Eq. (1.10)— with respect to the control volume average solution of the primitive vari-ables is calculated: ∂RCoef∂PV ar. This matrix C = A† = AT(AAT)−1is the pseudo-inverse1 of the reconstruction matrix, Eq. (1.14).– Next, the derivative of the reconstructed solution with respect to the recon-struction coefficients, R = ∂RSol∂RCoef, is calculated. This matrix, evaluated at theGauss quadrature points, contains the geometric terms from the Taylor seriesexpansion of Eq. (1.10) which are functions of the distance between the Gaussquadrature points and the reference location of the control volume.– The Jacobian of the numerical flux function with respect to the reconstructedsolution, F = ∂Flux∂RSol, is calculated along with the flux using automatic differ-entiation.– The last part of the chain rule is the derivative of the flux integral with respectto the numerical flux at each side, Q = ∂FluxInt∂F lux, which is simply the appropriateGauss quadrature weight.1.1.6 Time Advance SchemesFor a time accurate solution, the residual equation, Eq. 1.8 has to be discretized andintegrated in time. Time discretization schemes generally fall into two broad categories1However, in practice we use singular value decomposition to calculate the pseudo-inverse of thereconstruction matrix, A.13of explicit and implicit schemes. In the explicit schemes as shown in Eq. 1.19, the fluxintegral is evaluated at the previous time step, n, to measure the solution at the nexttime step. In this way all the time derivative equations are decoupled and thus theimplementation of these schemes are much easier. However, these schemes are limited inthe size of the time step, 4t, by the Courant-Friedrichs-Lewy (CFL) stability constraintand its analogs for diffusion-dominated problems. This constraint makes these schemesvery slow to converge to steady state. The simplest formulation of explicit schemes,shown in Eq. 1.19, is called Explicit Euler, and is only first order accurate in time.Un+1 − Un4t = <(Un) (1.19)Conversely, implicit time advance schemes evaluate the flux integrals at time stepn + 1, at which the solution is not yet known. Eq. 1.20 shows the Implicit Euler timeadvance scheme, which is first order in time:Un+1 − Un4t = <(Un+1) (1.20)This scheme is not bounded by CFL constraints. As such, much larger time steps canbe taken, which in turn reduces the number of iterations for steady state convergence. Tosolve Eq. 1.20, the nonlinear residual is linearized by Taylor series expansion as shownin Eq. 1.21.Un+1 − Un4t = <(Un) + ∂<∂UδU +O(δU2)(1.21)(I4t −∂<∂U)δUn = <(Un) (1.22)Un+1 = Un + δUn (1.23)In Eq. (1.22), the time step ∆t is calculated from the global CFL. I is the identity14matrix. The resulting linear system (Eq. 1.22) is solved to compute the solution updates,δU , where ∂<∂Uis the sparse Jacobian matrix. As a result, implicit time advance schemesbind all the time discretization equations together, which in turn make the implemen-tation of these schemes much harder due to the fact that all the equations have to besolved simultaneously.However, for steady state solutions, which are the focus of this thesis, there is no needfor time accurate solutions. As such, the non-linear residual equation <(U) = 0 could inprinciple be solved directly using Newton’s method:∂R∂UδU = −< (Un) (1.24)Us = U0 + δU, (1.25)where the solution at the steady state, Us, can be reached by only a few solution updates,with convergence dependent only on the quality of the linearization. However, Newton’smethod is strongly dependent on the initial solution; an initial solution too far from theconverged solution can cause the method to diverge. Pseudo-transient continuation offersa middle ground between these two, in which an implicit time stepping scheme is usedin combination with local time stepping, under-relaxation, and rapidly increasing CFLnumber.(I4tnα− ∂R∂U)δU = <(Un) (1.26)Un+1 = Un + ωnδUn, (1.27)where the solution is updated partially using an under-relaxation factor, ωn, at eachpseudo time-step to ensure that the residual after each iteration does not increase andthe solution is physical in the sense of physically realizable not physically accurate. This is15achieved by calculating ωn using a line search algorithm which in turn is an optimizationproblem in the space of the solution update vector. This line search algorithm finds ωnto approximately minimize the residual at the next time step. Calculating the linearsystem of Eq. 1.26 at each pseudo time-step directly is computationally expensive. Forthis reason, an iterative method such as GMRES (see [16, 28] for more details) is used tosolve this linear system and to calculate the solution update at each pseudo time-step.Additionally, time accuracy of the solution is not sought in the transient solutionsof a pseudo transient continuation method, and thus bigger time steps and local timestepping can be used to converge to the steady state solution. Hence, the local time step,4tnα, for every control volume, which is calculated from the global CFL, a local lengthscale, hα, and the maximum local wave velocity, cmax,α, is:4tnα =CFLn hαcmax,αFor early time steps, the initial solution does not conform to the boundary conditions.Thus, strong transient errors start propagating through the flow field if the initial timesteps are large. For this reason, the very first time steps are kept small while the latertime steps have to rapidly increase to expedite the convergence rate. This is controlledby an exponential growth of CFL and some stopping criteria which was proposed by [29].Moreover, for visualization reasons that will be explained more in Chapter 2, some-times we use a simpler time-stepping method without any line search algorithm or localtime stepping: the Crank-Nicolson implicit time stepping scheme:Un+1 − Un4t =12[<(Un) + <(Un+1)](1.28)where Eq. 1.28 is linearized around the current solution state, Un:Un+1 − Un4t =12[<(Un) + <(Un) + ∂<∂UδU +O(δU2)](1.29)16By discarding the second order error terms, Eq. 1.29 is reduced to:(I4t −12∂<∂U)δUn = <(Un)which is the equation to be used when calculating the solution update vector using theCrank-Nicolson method.1.2 Numerical Stability MethodsThere are multiple ways to determine the stability of the numerical schemes. Section 1.3will focus on energy stability, which is the method used in this thesis. In this section, Idescribe other common methods in the literature. One of the early works on the stabilityof numerical finite difference schemes [30], denotes a method as stable if the numericalscheme does not allow the errors to increase indefinitely. This means that for any timestep, the errors have to bounded and not amplified. If for a discrete solution, Un, at timestep n, there is an exact solution U , then the discretized solution satisfies the followingequality:U = U + (1.30)where is the discretization error which is the difference between the discrete solutionand the exact solution. Hence, this definition of stability requires that discretizationerror be bounded at all times for all the control volumes. In other words, this stabilitycondition can be re-written as:‖ni ‖ ≤ K, (1.31)for every control volume i, and time step n, where K is a constant number which isindependent of n. However, this condition does not necessarily ensure that the error17will remain small enough for the solution to be physically valid. For this reason, Laxand Richtmyer [31] proposed a more general numerical stability condition whereby thenumerical stability is based on the evolution of the solution in time instead of the errorto the solution. In this stability definition, no components of the initial smooth solutionshould be amplified beyond bounds. Mathematically put, for any time step, n, thefollowing condition should hold:‖Un‖ ≤ K∥∥∥U0∥∥∥ (1.32)where U0 is the initial solution and K is a constant positive real number. For a linearsystem of equations which results in a linear discretization operator, the discrete solution,Eq. 1.30, evolves in time according to a linear operator, L:Un+1 = LUn (1.33)where operator L is a function of the time and space discretization parameters. If theinitial solution at time t = 0, is denoted as U0, then the solution after the first time stepis:U1 = LU0 (1.34)Similarly, the solution at the second and the nth time steps are respectively:U2 = LU1 = L2U0 (1.35)Un = LnU0 (1.36)For the solution at time step n, Un, to be bounded, given that the initial solution issmooth and bounded, the multiplications of the linear operator has to be bounded by18some arbitrary positive real number K:‖Ln‖ ≤ K (1.37)Eq. 1.37 holds for any arbitrary value of n and a constant value of K, if and only ifthe norm of the linear operator is always smaller than one:‖L‖ ≤ 1 (1.38)Eq. 1.38 has inspired many other works on numerical stability, some of which areexplained in the subsequent sections. In the next sections, I briefly describe the analysis ofstability through Von Neumann stability analysis, matrix stability, and spectral stability.1.2.1 Von Neumann Stability AnalysisOne of the earliest and most commonly used techniques for stability analysis of discretenumerical schemes is the Von Neumann stability analysis. This algorithm, which wasdeveloped during World War 2 by Von Neumann but kept confidential until it firstappeared in [32], is restricted to infinite or periodic domains. Hence, the effects of theboundaries have to be neglected if one chooses to use this approach, and the errors ofthe discretization in the interior domain is modeled by Fourier series terms.This scheme works by inserting Eq. 1.30 into the linear system of equations, Eq. 1.33,to obtain:Un+1 + n+1 = LUn + Ln (1.39)Since, the exact solution, U , satisfies the linear discrete equation, Un+1 = LUn, Eq.1.39 becomes:n+1 = Ln (1.40)19which signifies that the discretization errors evolve in time in the same way that thediscrete solution does. Hence, to have a bounded solution, the error has to be bounded,which means that as time goes on the discretization error should not increase. Now ifone models the errors by Fourier series and harmonic terms, the error at any time stepfor any control volume is comprised of multiple harmonic terms:n =∑jEnj eiθj (1.41)where Enj and θj are the amplitude and frequency of the jth harmonic term. Providedthat this method is used for a linear scheme, every harmonic term of the error shouldalso satisfy the discretization equation, Eq. 1.40:En+1j = LEnj (1.42)Using the stability condition of Eq. 1.38, this is only achieved if for every subsequenttime step and for every harmonic, the amplitude of the error harmonic, Enj , does notincrease:Gj =∥∥∥∥∥En+1jEnj∥∥∥∥∥ ≤ 1 (1.43)where Gj is the amplification factor for mode j, and is a function of the time step,mesh size and physical properties such as wave speed. A common approach in usingthe amplification factor Gj for all modes j is to calculate the stability curves of Gjin the complex plane and ensure that for all such curves the curve is within the unitcircle, indicating stability. This method usually leads to conditional stability with somerestrictions on the Courant-Friedrichs-Lewy (CFL) parameter. The CFL condition, whichwas first introduced in [33], requires that the numerical scheme determining Un+1 fromthe previous time steps has to have all the physical information. For example, for a onedimensional convection problem, this means that for any time interval ∆t, the distance20which is traveled by a wave should not exceed the length of a control volume.The Von Neumann stability analysis is limited to linear problems. Thus, in the caseof non-linear systems, the stability information may not be fully extended to the generalnon-linear problem as one has to investigate stability locally by freezing non-linearityand non-constant coefficients.1.2.2 Matrix Stability Analysis:The Von Neumann stability analysis does not take into account the boundary effects andit is only suitable for either the interior of the domain or periodic boundary conditionswhich are rarely practical cases. While resembling the previous stability analysis insome respects, the matrix stability method briefly described in this section considers theentire domain, including boundaries. As the name suggests, this method models thenumerical scheme as a matrix equation. If there is a linear operator such as ∂u∂u= L (u)for the solution vector u, after spatial discretization and using the method of lines, thislinear partial differential equation is transformed into an ordinary system of differentialequations:dUdt= JU + S (1.44)where as explained earlier, J is the Jacobian matrix, U is the discrete solution vector,and the source term S is the non-homogeneous part of the spatial discretization. If forthe N × N matrix J , where N is the degrees of freedom (number of control volumesmultiplied by the number of components of the solution) of the problem, there is acomplete eigensystem2 with N eigenvalues and N eigenvectors, the homogeneous solutionto the Eq. 1.44 can be written as a linear combination of the eigenvectors:2Explained in Section 1.3.221U =N∑j=1Uj (t)Rj (1.45)In Eq. 1.45, Rj is the jth right eigenvector of the matrix J where the jth right eigensystemwith eigenvalue λj is:JRj = λjRj (1.46)Similarly, the source term S can be written as:S =N∑j=1SjRj (1.47)Now, substituting Eq. 1.45 and 1.47 into Eq. 1.44 and using the the right eigensystem,one reaches the following ODE for every component of the solution:ddt N∑j=1Uj (t)Rj = N∑j=1Uj (t) JRj +N∑j=1SjRj (1.48)ddtUj (t)Rj = Uj (t)(JRj)+ SjRj (1.49)ddtUj (t) = Uj (t)λj + Sj (1.50)Additionally, we know that the ordinary differential Eq. 1.44 has a homogeneous, Uh,and particular, Up, solution of:Uh = Cjeλjt (1.51)Up = −Sjλj(1.52)where Cj is a constant based on the initial solution, U0j , and can be calculated as Cj =U0j+Sjλj. Hence, the solution of the Eq. 1.44 is:22U =N∑j=1[U0j eλjt + Sjλj(eλjt − 1)]Rj (1.53)Now, for the solution, U , to be bounded and to meet the condition of Eq. 1.32 at anytime, when the initial solution is smooth and bounded, the real parts of the eigenvalueshave to be negative. This clearly states that the stability of the Eq. 1.44 is fully deter-mined by the Jacobian matrix of the spatially discrete system and its eigensystem; hencethe homogeneous part only. Therefore, if the space discretization generates a positiveeigenvalue, there is no guarantee that any time stepping method will provide a stablesolution. Additionally, when one is looking for the steady-state solution using time step-ping methods, they can infer that the convergence rate to the steady state solution isdirectly related to the real part of the eigenvalue, R (λj). Transient solution componentswith eigenvalues of larger negative real parts decay much faster than the ones with eigen-values of smaller negative real parts. Therefore, the convergence rate is dominated by thenegative eigenvalues closer to zero (the right most eigenvalues). Since the homogeneouspart of the discretization is responsible for stability if one uses the homogeneous part ofthe discrete solution, Un = ∑Nj=1 Unj Rj, in the stability condition of Eq. 1.33:N∑j=1Un+1j Rj = LN∑j=1Unj Rjand using the the stability condition of Eq. 1.36 on every solution component, Unj , onecan identify the amplification factor as a matrix format:‖GL (λ,∆t)‖ ≤ 1However, we know from linear algebra that the pseudo-spectral radius of a matrix isalways equal or smaller than the vector induced norm (such as L2 norm) of the matrix:% (GL (λ,∆t)) ≤ ‖GL (λ,∆t)‖ ≤ 1 (1.54)23(a) Fully-discrete scheme (b) Semi-discrete schemeFigure 1.2: Stability regionswhere the pseudo-spectral radius of a matrix is defined as the maximum of the modulusof the eigenvalues of the matrix:% (GL) = max{∥∥∥λ0GL∥∥∥ , · · · , ∥∥∥λNGL∥∥∥} (1.55)Eq. 1.54 brings us to the pseudo-spectral stability condition which states that thenumerical scheme is stable if and only if the pseudo-spectral radius of the scheme is lessthan one. In other words, all the modulus of all the eigenvalues of the fully discretesystem have to be smaller than or equal to one. Figures 1.2a and 1.2b show the stabilityregion in the fully discrete scheme and the spatially discrete scheme respectively.It is also shown in the literature (see [34, 35]) that the matrix stability of dynamicalsystems of equations is equivalent the classical Von Neumann stability analysis on aperiodic Cartesian gird. Moreover, the extension of matrix stability to nonlinear problemsis the entropy stability analysis explained briefly in the next subsection. For more aboutmatrix stability, see [36, 37].1.2.3 Entropy StabilityIt is shown in the literature that the solution to the system of conservation laws, suchas Eq. 1.1, can become non-smooth even if initialized by a smooth solution (e.g. see24[38]). For this reason, the weak form solution (the integral form) may not be unique. Anadditional criterion such as entropy conditions often has to be met to ensure uniquenessof the solution. For example, for a one dimensional form of Eq. 1.1 with no source term,the conservation law is:∂u(x, t)∂t+ ∂F (x, t)∂x= 0 (1.56)Now, if one assumes that for the system 1.1, there exists an integrability equalitysuch that:∂q(u)∂x=(∂η(u)∂x)T∂F (x, t)∂x(1.57)where q is the entropy flux function and η is the entropy function, while χ = ∂η(u)∂xis theentropy variable3. If the Eq. 1.56 is left-multiplied by χT , the so-called conservation lawfor entropy in the smooth regions in found:∂η(u)∂t+ ∂q(u)∂x= 0 (1.58)However, in the presence of discontinuities (shocks), entropy is dissipated which resultsin the following entropy inequality:∂η(u)∂t+ ∂q(u)∂x≤ 0 (1.59)Hence, if one integrates the entropy inequality in space, an upper bound for the totalentropy is found:∂∂tˆη(u)dx ≤ 0 (1.60)which states that the total entropy is non-decreasing with respect to time. In other3This is not the thermodynamic entropy but mostly an analogy.25words, given the first initial conditions, u0, the inequality of Eq. 1.60 results in an apriori upper bound for the solution of the system 1.56:ˆη(u)dx ≤ˆη(u0)dxAs mentioned earlier, this is similar to the L2 stability. For further information, see[39, 12, 11] and the references therein.1.2.4 L−Stable and A−Stable MethodsA time integration method for system of equations of dUdt= JU , is A-stable if for anypositive time step and negative semi-definite Jacobian matrix of J , the amplificationfactor of the system, Eq. 1.43, is not greater than one. In this way, if the semi-discretesystem of equations of the numerical method are stable, using the method of lines, thenthe full discrete system of equations are also stable. An L-stability is an extension of theA-stability where the numerical method, in addition to being A-stable, has to have anasymptotic amplification factor of zero. This means that when the time step is indefinitelylarge, the amplification factor is zero (see [40] for more details).For example, the trapezoidal method is A-stable but not L-stable, as its stability re-gion contains the entire left half of the complex plane. This time integration scheme, doesnot damp any mode of the solution that does not need to be decayed. The amplificationfactor of the trapezoidal time integration approaches one for some components of thesolution with large eigenvalue modules (i.e. when the eigenvalues of the spatially discretesystem are positive). This borderline stability causes these stiff modes of the solutionto decay slowly which in turn are shown as spurious oscillations in the solution of thenumerical method. In other words, any eigenvalues of the spatial discretization which arepositive will remain unstable, while any negative eigenvalues of the spatial discretizationwill remain stable with this time discretization scheme. Hence, this method is used as a26sanity check in this thesis to ensure that a method is stable only because our stability al-gorithm has made it stable and not because our sophisticated time stepping algorithm ishiding the instability. The implicit Euler scheme and its derivatives are L-stable schemes.1.3 Steady State Numerical Stability in this ThesisAs explained earlier, a numerical scheme that does not allow the errors or the non-physical modes to grow without any bound is called stable. In this thesis for the purpose ofstability, I implemented method of lines (see [34]) which uses only the spatially discretizedsystem of equations, Eq. 1.8, to determine stability of finite volume methods at steadystate regardless of the choice of time stepping method. In this way, the complexity oftime stepping methods for determining the stability of the numerical methods at steadystate can be avoided. Additionally, using the Lax equivalence theorem (see [31, 41]),I use the convergence as a proxy to the stability of the numerical methods. The Laxtheorem, which originally states that for a linear finite difference method, consistencyand convergence are necessary and sufficient conditions for the stability of the method,has been extended to the nonlinear finite volume method (see [42, 43, 44, 45]) for whichconsistency is not a necessary condition anymore. Hence, this fundamental theorem ofnumerical analysis states that convergence and stability of the numerical methods areequivalent. A numerical scheme is said to be consistent if the discretized equationstend to the corresponding continuous equations for sufficiently small space and timediscretization, where this condition applies to any discretization scheme of accuracy ofmore than zero. A common method to show the consistency of a discretization schemeis to use the Taylor series expansions of the discrete solution and substitute them in thediscretized set of equation to arrive at the continuous set of equations. The remainingterms are said to be truncation errors. For a consistent continuous discretization scheme,the truncation error goes to zero when ∆x and ∆t approach zero. Likewise, if the solution27to the discretized equations tend to the exact solution of the corresponding continuousequations for sufficiently small space and time discretization, the numerical scheme isconvergent. The amount to which the discrete solution does not match the exact solutionof the continuous equations is called the discretization error.The core stability analysis of this thesis is based on energy stability which herein isused in two different ways: I seek to directly minimize the growth of quadratic energyfor different discretization schemes. Additionally, I use eigenvalue analysis of the Jaco-bian matrix of semi-discrete scheme as a sufficient condition for energy stability. Thesemethods are described in the following subsections.1.3.1 Energy StabilityIn energy stability analysis (e.g., see [46, 47]), an energy4 is defined based on the solution,and then stability is proved by showing that this energy is not increasing. The (linear orlinearized) PDE is discretized in space, producing a coupled set of ordinary differentialequations, such as Eq. (1.9). The energy associated with the vector of CV solutionaverages as functions of time ~Uα(t) is defined by:E (t) = 12∑αUTα (t)Uα(t) (1.61)The rate of change of energy is then:dEdt=∑αUTαdUαdt=∑αUTα JUα (1.62)For stability, the energy should not be increasing(dEdt≤ 0).Moreover, energy stability can also be expressed using the Lyapunov theorem [48]which proves that the two following statements for a symmetric matrix are equivalent:1. All eigenvalues λ of J have non-positive real parts, < (λ) ≤ 0.4This is not the thermodynamic energy but merely an analogy282. There exists a positive definite matrix M such that matrix MJ + JTM is negativesemi-definite.Following the energy stability analysis and the second expression of the extension ofLyapunov theorem (see [48] for more details), for symmetric Jacobian matrices, it issufficient to find a positive definite matrix M such that the matrix MJ + JTM (with Jbeing the Jacobian of the semi-discrete systems of equations) is negative semi-definite.However, it is difficult to find such a matrix M for an arbitrary system of equations.Instead, Haider et al. [7] suggested finding a strategy to decrease the growth of energy,at the very least, to gain better stability properties. The growth of energy was rewrittenas:ddt(U,MU) =∑α(UTαMJUα + (JUα)TMUα)= 2∑α(UTαMJUα)(1.63)where a positive diagonal matrixM with elementsMαα = |να| is used. να is the size of thecontrol volume α. This makes the calculations invariant under mesh scaling. However,Eq. (1.63) does not hold for general cases especially when the Jacobian matrix is notsymmetric or when the partial differential equation is not homogeneous of degree one.For these reasons, I use a weighted form of the growth of the quadratic energy instead,Eq. (1.62), which is:ddt(U,MU) = dEMdt= 2∑α(UTαMdUαdt)(1.64)The quadratic energy equation, defined by Eq. (1.64), does not require a partialdifferential equation to be homogeneous of degree one or the numerical scheme to havea symmetric Jacobian matrix. This quadratic energy equation is further developed inChapter (3) for different flow fields and numerical schemes.291.3.2 Eigenvalue AnalysisAs stated earlier, for stability, the energy should not be increasing(dEdt≤ 0). In general,this is possible if and only if the matrix J is negative semi-definite. As by definition, amatrix J is negative semi-definite if and only if for every non-zero vector U such that:UTJU ≤ 0 (1.65)where Eq. (1.65) is equivalent to the matrix expression in Eq. (1.62). This in turnconfirms that in order for the energy to not increase the Jacobian matrix has to benegative semi-definite.Moreover, if the matrix J is symmetric, the semi-discrete system is stable if and onlyif the eigenvalues of the global Jacobian J have non-positive real parts. However forgeneral asymmetric Jacobian matrices, the system may be stable with some eigenvalueshaving positive real parts but it is stable if all eigenvalues have non-positive real parts(see [49] for more details).For general Jacobian matrices, I satisfy the sufficient condition for energy stability bymoving the positive rightmost eigenvalues of the Jacobian matrix to the negative (left)half of the complex plane. In this way, I may be unnecessarily changing the rightmostpositive eigenvalues of an already stable problem to guarantee energy stability. Yet, bydoing this extra work I am able to automatically reach energy stability for every problem.This energy stability is used in all the subsequent chapters.1.3.2.1 Eigenvalue ProblemEigenvalues and eigenvectors are numbers and vectors respectively that will specify thecharacteristics of a square matrix Jn×n. The right and left eigenvectors can be sought bysolving the right, Eq. 1.66, and the left, Eq. 1.67, eigenvalue problem, respectively:30J(−→b )R(−→b ) = ΛR(−→b ) (1.66)L(−→b )J(−→b ) = L(−→b )Λ (1.67)where Λ is a diagonal matrix of eigenvalues, λ, of the global Jacobian J , and the collectionof right and left eigenvectors of J are in columns and rows of the R and L matricesrespectively. In general, the right and the left eigenvectors are not equal (R 6= L) sincethe Jacobian matrix, J , is not symmetric. The vector −→b represents all the independentvariables in the Jacobian matrix J and its eigenpairs. Taking the transpose of the lefteigenvector problem will result in Eq.JTLT = L(−→b )Λ (1.68)Comparing Eq. 1.66 and 1.68 shows that the left eigenvector of matrix J is the righteigenvector of the transposed matrix, JT , while the eigenvalues are the same. In thisthesis, I use this fact to calculate both the right and the left eigenvectors of the Jacobianmatrix using software that computes the right eigenvectors of its input matrix.1.3.2.2 Eigenvalue Calculation AlgorithmsCalculating all the eigenvalues of Eq. 1.66 and 1.67 directly (e.g., using the LQ andQR factorization algorithms [50]) is infeasible or impossible as the process is intensivelymemory expensive (it is cubic with respect to the degrees of freedom). Instead, I useiterative methods such as Krylov subspace method (see, for instance, [51] for a summary)for parts of the spectrum that are of interest. These iterative methods converge fasterfor eigenvalues that are well separated and degrade or diverge for closely packed eigen-values. For these reasons, matrix transformations such as the shift-invert and Cayleytransformations (see [52] for transformations of this type) are used to map the lumped31eigenvalues to well-separated regions. To solve the eigenvalue problem in my calculations,I use the Scalable Library for Eigenvalue Problem Computations package [53], SLEPc,which has support for iterative solvers as well as spectral transformations. Also, SLEPcis consistent with the linear algebra package incorporated in ANSLib, PETSc [54], whichmakes it particularly easy to integrate. Unfortunately, the iterative solvers in SLEPcdo not support the calculation of the left eigenvectors. To rectify this problem I cal-culate the right eigenvectors of the transposed Jacobian matrix. More information onhow to calculate the rightmost eigenvalues of a sparse eigensystem with tightly clusteredeigenvalues is provided in the Appendix (see Section(6.6)).1.4 Surrogate Solution: Potential FlowAs will be explained in Section 2.2, in order to investigate the numerical stability at thesteady-state conditions, one needs to have the steady state solution. However this is notpossible if there is an instability problem. In this thesis, I intend to predict the stabilitybefore attempting to solve the principal problem. To this end, I use surrogate solutionsto approximate a more complicated physical problem at the steady-state conditions.Hence, in this thesis I am using the potential solution as an approximation of the inviscidcompressible flow solution at the steady state.In recent years, there has been a rising interest in using potential flows as a mediumfidelity analysis tool, as they are significantly cheaper in terms of computational timeand required resources. Hence, to obtain a fast and appropriate estimate of the steadystate solution to Euler problems, I use the solution of either the linear potential equationor the full potential equation.For flows with small disturbances over obstacles, the full potential equation can bereduced to the linear Prandtl-Glauert equation. Thus, for flow with small disturbancesaround a solid body, the velocity flow can be measured as:32−→V = U∞iˆ+∇ϕ (1.69)where U∞ is the free stream velocity assumed in the x direction (although if not, a rotationof coordinates is first used) and ϕ is the potential field. Then the Euler equations forinviscid, steady flow will be:ρ∞∇ · ∇ϕ+−→U ∞ · ∇ (ρ− ρ∞) = 0 (1.70)ρ∞−→U ∞ · ∇ (∇ϕ) +∇ (p− p∞) = 0 (1.71)where p∞ and ρ∞ are the free stream pressure and density respectively. By the assump-tion of the isentropic flow, pressure and density can be calculated as dp = c2dρ ≈ c2∞dρwhere c∞ is the free stream sound velocity. Eliminating the p, ρ from 1.70 and 1.71 leadsto the single equation:∇ ∇ϕ−−→U ∞c2∞(−→U ∞ · ∇ϕ) = 0 (1.72)where the free stream Mach is M∞ = U∞c∞ . Since the free stream velocity is assumed tobe in the x direction, Eq. 1.72 is reduced to:(1−M2∞)ϕxx + ϕyy = 0 (1.73)where Equation 1.73 is solved by the well-known panel methods. With a change ofvariables such as:x´ =√1−M2∞x, y´ = y (1.74)Eq. 1.73 transforms to a linear Poisson equation:33ϕx´x´ + ϕy´y´ = 0 (1.75)where the solution to this equation can be superimposed by different potential solutions ofthe source, doublet, and vortex singularities. In this method the surface of the boundaryis divided into segments (see Figure 1.3) where each segment has a source, doublet, orvortex singularity attached to them. This way, the problem of finding the solution forthe entire domain is replaced by finding the solution singularities on boundary segments.Additionally, the partial differential equation is simplified to a linear integral equationsuch as:ϕ = U∞x´ cosα + U∞y´ sinα +N∑j=1w [qj (s)2pi ln r −γ2piθ]ds (1.76)when the singularities are of type source and vortex. In Eq. 1.76, the first two terms onthe right hand side of the potential solution are the free stream components with α beingthe angle of attack. The first and the second integrands in Eq. 1.76 are the potentialsolutions from the source and vortex singularities respectively, where q (s) is the sourcestrength, γ is the vortex singularity strength, s is the panel length, and θ = arctan(y´x´)and r =√x´2 + y´2 are the polar coordinates. There are different ways to implement apotential solution by either using different singularities, different strength of singularitiesor different geometries of singularities (such as curved panels). In this implementation,I use straight panels, with different source strength for each panel while vortex strengthis the same throughout. In this formulation, the source terms are means of satisfyingthe flow tangency condition on the surface, while the constant vortex determines thecirculation around the airfoil by satisfying the Kutta condition. The Kutta conditionrequires that the flow of air from the top part of the trailing edge must leave the surfacesmoothly with the flow from the bottom part of the trailing edge. In practice, thiscondition is normally met by enforcing the same pressure at the upper and lower parts34Figure 1.3: Panel method for a NACA 0012 airfoilof the trailing edge. In the case of potential flows, this condition is met by ensuringthat the tangential velocities from the upper and the lower parts are equal (the normalvelocity is already satisfied by the conservation of mass). This potential field has Nsource strength terms and one vortex strength term to be found which makes a problemof N + 1 unknowns. Figure 1.3 shows the panel segments on a NACA 0012 airfoil. Thepoints shown in Figure 1.3 are the mid-sections of the panels as well as the trailing edge.Hence, to solve the N+1 unknowns of the potential field, I use these N mid-section pointsby satisfying the tangential velocity on the panel at each of these points. Additionally Iuse the Kutta condition as explained earlier for the last equation:Vi · ni = 0 for every i, i = 1, · · · , NV1 · n1 = −VN · nN Kutta condition(1.77)where Vi is the velocity at the ith panel, which is calculated by taking the gradients ofthe potential field V = ∇ϕ, and ni is the normal vector.As explained, the boundary surface is divided into panels formed by the boundarydiscretization. Through the integration of source and doublet singularities on these panelnetworks, the potential solution anywhere in the domain can be calculated. Hess andSmith [55] developed the first three dimensional panel method for incompressible flows.Later, others improved and extended the panel methods to be applicable for compressible35flows ranging from subsonic to supersonic flows (e.g. see [56, 57, 58, 59, 60]). Advancedalgorithms of these families have been coded for subsonic and supersonic flows in PANAIR [56, 61] and for transonic flows in TRAN AIR [62] by Boeing.However a more accurate solution can be obtained by solving the full potential equa-tion. Historically (e.g., see [63, 64, 65, 66]), the full potential problem which describesconservation of mass and enthalpy with constant entropy was reduced to one equationin the conservation form. In this implementation, the density is first solved using theinviscid, irrotational, and energy equation (Bernoulli):∇(∂ϕ∂t+ 12∥∥∥−→∇ϕ∥∥∥2 + pρ+ gy)= 0 (1.78)where by assuming negligible gravity effects and steady-state conditions and also usingthe free stream conditions, Eq. 1.78 is reduced to:12∥∥∥−→∇ϕ∥∥∥2 + pρ= 12U2∞ +p∞ρ∞(1.79)If one eliminates pressure, p, from Eq. 1.79 using the isentropic relation, dp = c2dρ ≈c2∞dρ, and defining free stream Mach as M∞ = U∞c∞ , Eq. 1.79 is solved for density as:ρ =[1 + γ − 12 M2∞(1−∥∥∥−→∇ϕ∥∥∥2)] 1γ−1 (1.80)Here γ is the specific heat ratio,M∞ is the free stream Mach number, and the velocityfield can be determined by −→U = −→∇ϕ. Eq. 1.80 is scaled with non-dimensionalized freestream density, ρ∞, and sound speed, c∞, of one. Finally this equation can be used inthe conservation of mass to lead to the one-equation model for potential flow:−→∇ · (ρ−→∇ϕ) = 0 (1.81)However, for transonic flows and the treatment of hyperbolic regions the conserved36potential equation (Eq. 1.81) must be modified for stability in the presence of shocksand discontinuities. This is done by biasing the value of isentropic density in the upwinddirection:ρ˜ = ρ− h˜−→∇ϕ∥∥∥−→∇ϕ∥∥∥ ·−→∇ρ (1.82)h˜ = µh√M2 − 1 (1.83)and using the upwind biased model:−→∇ · (ρ˜−→∇ϕ) = 0 (1.84)where h is the local length scale in the mesh, µ is a hyper parameter in the range of0.8 − 0.9, and M is the local Mach number. Equation 1.84 is then solved by finitevolume methods after applying the divergence theorem. The Kutta and wake boundaryconditions are also satisfied by imposing a zero density at the trailing edge jump, andmass flux continuity along the wake sheet respectively. Eq. 1.84 is used in a finite volumediscretization with flux of ρ˜−→∇ϕ and a source term of zero.For better stability properties, Parrinello and Mantegazza [67, 68] introduced a newtwo-equation model for full potential flow by defining ρ and ϕ independently to satisfyconservation of mass and enthalpy respectively. This new formulation allows settingup the numerical solution more easily. Galbraith et. al. [69] further improved andimplemented this model for finite element methods. The two equations to be solved forρ and ϕ are:37ρt +∇. (ρ˜∇ϕ) = 0 (1.85)ϕt +12(|∇ϕ|2 −M2∞)+ M−2∞γ − 1k ( ρ˜ρ∞)γ−1− 1 = 0 (1.86)where the upwinded density is obtained using Eq. 1.82. For supersonic flows, along withthe two equations 1.85, and 1.86 a transport equation has to be solved for entropy (k):kt +∇ϕ ∇k = 0The continuity of pressure and enthalpy at wake sheets is imposed by enforcing Eq. 1.87and 1.88 above and below the wake sheets:(kργ)lower = (kργ)upper (1.87)(ϕt +12 |∇ϕ|2)lower=(ϕt +12 |∇ϕ|2)upper(1.88)Similarly, both Eq. 1.85 and 1.86 are solved using a cell-centered finite volume dis-cretization with boundary conditions, with Eq. 1.87 and 1.88, enforced weakly as bound-ary flux functions. Section 2.2 describes how these three potential models are used tocalculate a surrogate solution for Euler problems to investigate the numerical stability ofthe finite volume discretizations.1.5 Objectives and NoveltyAll previous methods of numerical stability focus on the mathematical proof of stability orreshaping space and time discretizations to meet certain stability thresholds. However,in this thesis I design a stabilization scheme that can detect the instability issues in38different aspects of finite volume problems and automatically stabilize these originallyunstable problems in practice. Additionally, a vast amount of engineering time andcomputational resources are normally needed in either a research or engineering work tofind a convergent solution in a trial and error manner. This thesis aims to address thisspecific issue of iterative trial and error by providing a systematic way whereby unstableproblems can be stabilized in practice.Thus, the first objective and novelty of this thesis is to use the right most eigenvaluesof the semi-discrete Jacobian and their gradient with respect to mesh vertex locationsto stabilize the originally unstable cell-centered finite volume methods on unstructuredmeshes. Often times, the instabilities arise from the unfavorable characteristics in themesh and so far the only approach to stabilization has been to regenerate the meshto meet certain characteristics. For the first time, I show in this thesis that many ofthese meshes are transformable by smoothing the location of mesh vertices through anoptimization scheme with the goal of stabilization. This stabilization is achieved eitherafter instability has been observed in practice, during the convergence, or even beforeattempting to solve the initial problem through using a surrogate solution.Secondly, I investigate the stability properties of the solution reconstruction whichis a pivotal part of the k-exact finite volume discretization of high order scheme onunstructured meshes. This investigation has two parts. In the first part, the noveltylies in extending both a mathematical and experimental stability assessment techniqueto nonlinear and high order problems from an existing work in order to show the effectof the solution reconstruction stencil on stability. In the second part, I introduce anew automated optimization scheme that can select the optimal solution reconstructionstencil for every control volume in the problem using eigenanalysis.Often times the hardship of convergence stems from unfavorable boundary conditionsand thus one needs to modify these boundary conditions to arrive at a stable solution.Hence, the last objective and novelty of this work is to design a scheme that can automat-39ically find the best hybrid boundary condition for every boundary face in the problem.1.6 OutlineThe rest of this thesis is organized as follows.Chapter 2 introduces the use of eigenvalue analysis in automatically optimizing thenumerical stability of a finite volume scheme by moving mesh vertices locally. Thischapter details the fundamental ideas of this thesis, including the calculation of thegradients of eigenvalues with respect to mesh vertex location through the gradients ofthe global Jacobian matrix, the effects of these gradients in projecting the eigenvalues ofa perturbed mesh, and the use of surrogate solutions in determining the stability a priorto solving the principal problem.Chapter 3 first shows how the solution reconstruction affects stability by using aquadratic form of energy stability analysis. Two different methods are then used toautomatically optimize the stencil of the solution reconstruction in order to stabilize anoriginally unstable problem. The first method, which is purely geometrical, uses a proxy,which is called the solution local reconstruction map, to find the best reconstructionstencil for every control volume in the problem. The second method uses gradients of theright most eigenvalues of the semi-discrete Jacobian with respect to the size and directionof the solution reconstruction stencil and an optimization scheme to stabilize unstablediscretizations.Chapter 4 shows how a hybrid boundary condition at every boundary control volumecan be obtained by optimizing different types of boundary conditions for the purpose ofnumerical stability, as simply optimizing the mesh vertices of the solution reconstructionmay not always lead to stabilization.Chapter 5 depicts the sensitivities of the right most eigenvalues of the semi-discreteJacobian on perturbations in the solution. This chapter is rather brief and is mainly to40confirm the applicability of the methods used in previous chapters. This is due to the factthat throughout Chapters 2, 3, and 4 I frequently use different perturbations on differentelements of the mesh, discretization, or solution to project the numerical stability of thefinite volume scheme.Finally, Chapter 6 concludes the main contributions of this thesis, and presents op-portunities for the future work.41Chapter 2Mesh Optimization to ImproveStabilityAs previously explained in Chapter 1 and more specifically in Section 1.3, to define thestability of a spatially discretized problem, the rightmost eigenvalues of the Jacobian ofthis semi-discretized problem are sought. The rightmost eigenvalues are the ones in theeigensystem with the largest positive real parts. This thesis and specifically this chap-ter hinges on using these rightmost eigenvalues, eigenvectors and their correspondinggradients to determine and to improve the numerical stability of the partial differentialequations in fluid dynamics spatially discretized using finite volume methods. This chap-ter first derives a methodology to automatically improve the numerical stability of thesesemi-discretized problems in Section 2.1. Second, this method is expanded to improvethe numerical stability both a priori (prior to solving the physical flow) and a posteriori(after solving the physical flow) in Sections 2.2 and 2.3 respectively. Finally, a study oncomputational cost of this stability method is provided in Section 2.4.2.1 Stability ImprovementFrom energy stability analysis in Section 1.3.1 and our definition of linear numericalstability in Section 1.3, a semi-discrete system of equations, Eq. 1.9, is stable if and onlyif the Jacobian matrix J = ∂R∂uis negative semi-definite; we assume the transient growthis negligible. To obtain the semi-discrete Jacobian matrix, J , we use an approximate42solution such as a lower order solution or a separate surrogate at the steady state tolinearize the Jacobian matrix for nonlinear problems. This chapter shows how theseapproximate solutions can be used in different scenarios to help with improving thestability of the second order problem.From the eigen-decomposition of the semi-discrete systems (e.g., see [70, 34]) for real-istic finite volume discretization one can indeed calculate and predict the stability char-acteristics for an implicit or explicit solver. This idea motivated us to use eigenanalysisto determine stability and for the first time use gradients of eigenvalues to predict howspectral stability will change with changes in the mesh. The innovation in this chapter isthe use of small perturbations in the location of mesh vertices to improve the stability ofsemi-discretized systems. To this end, I am using small movements in the mesh to guidethe unstable eigenvalues of these semi-discrete systems to the stable regions.To better understand how small perturbations can drive eigenvalues and to explain themethodology, perturbations of eigenvalues are briefly explained in Section 2.1.1. Section2.1.2 uses the perturbation theorem to derive derivatives of eigenvalues with respect tothe location of mesh vertices. In Section 2.1.3, an optimization scheme is introducedwhich uses the gradients of the eigenvalues to find the best mesh movement vectors forstabilizing an unstable problem by moving mesh vertices. Sections 2.1.4 and 2.1.5 arefocused on applying this stability scheme using a small subset of mesh vertices includingvertices located on boundaries.2.1.1 Eigensystems and PerturbationsBased on the perturbation theory, some calculations are described here to explain howsmall perturbations in an underlying system of equations can in turn change the eigen-values of the system. For an n× n (Jacobian) matrix J , where λ, R are respectively aneigenvalue, which is not repeated, and eigenvector of JR = λR, a family matrices can bedefined such as:43J˜ () = J + J (1) (2.1)where is infinitesimally small and∥∥∥J (1)∥∥∥ ≤ 1 is the perturbation matrix. Similarly theeigensystem of the perturbed matrix J˜ is defined as:J˜ () R˜ () = λ˜ () R˜ () (2.2)where its right eigenvector R˜ and eigenvalue λ˜ are analytical functions of the right eigen-vector R and eigenvalue λ of the unperturbed matrix J :λ˜ () = λ+ λ(1) + 122λ(2) + · · · (2.3)R˜ () = R + R(1) + 122R(2) + · · · (2.4)Now by substituting Eq. 2.1, 2.3 and 2.4 into Eq. 2.2, one can show that the followingequations should hold:O(1) JR = λRO() JR(1) + J (1)R = λR(1) + λ(1)RO (2) 12JR(2) + J (1)R(1) = 12λR(2) + λ(1)R(1) + 12λ(2)R(2.5)Eq. 2.5 is rewritten to:(J − λI)R(1) =(λ(1)I − J (1))R (2.6)The adjoint of matrix J is denoted as J∗. Since matrix J is a real non-symmetricmatrix, its transpose and adjoint matrices are the same J∗ = JT . Likewise, J , J∗ andJT all have the same eigenspace (i.e., λ = λ∗ = λT ). The right eigenproblem for JT isdefined as:44JTLT = λLT (2.7)where taking a transpose of Eq. 2.7 yields:LJ = λL (2.8)Eq. 2.8 in turn is the left eigenproblem of the matrix J , where L is the left eigenvector.Now multiplying Eq. 2.6 by the left eigenvector yields:L (J − λI)R(1) = L(λ(1)I − J (1))R (2.9)where by using Eq. 2.8 it can be proved that the left hand side equals to zero:L (J − λI)R(1) = LJR(1) − LλR(1)= λLR(1) − LλR(1)= 0Thus the right hand side of Eq. 2.9 which is used to measure the first order perturbationsof the eigenvalue is rewritten to:λ(1) = LJ(1)RL ·R (2.10)From the linear algebra it can be proved that the left and right eigenvectors of an eigen-value are not orthogonal. Hence, using Eq. 2.1 and 2.10, an upper bound for theperturbations of the eigenvalue can be found:∥∥∥λ(1)∥∥∥ =∥∥∥LJ (1)R∥∥∥‖L ·R‖ ≤‖L‖∥∥∥J (1)∥∥∥ ‖R‖‖L ·R‖ ≤‖L‖ ‖R‖‖L ·R‖ = κ (λ) (2.11)45where κ (λ) is called the condition number of the eigenvalue. The next section shows howgradients of eigenvalues with respect to the location of mesh vertices can be calculatedby using small local changes in the mesh vertex locations ξ as infinitesimal perturbations.2.1.2 Eigenvalue Derivatives With Respect To MeshMovementThe key and the main innovation in this chapter is the ability to predict changes in theeigenvalues with respect to changes in the mesh. Following on the previous section andusing infinitesimal perturbations on the matrix to calculate changes in the eigenvalue, Iintroduce mesh vertex movement as the driver of the perturbations of the Jacobian matrixand in turn eigenvalues. This means that by moving the vertices, the Jacobian matrixchanges, which then changes the eigensystem. However, different methods can be used toapproximate the gradients of the eigenvalues such as directly applying finite differences onthe eigenvalues or by reverse automatic differentiation. The former is an excruciatinglyexpensive process where the eigenvalues of the system have to be recalculated for eachsmall perturbation. However, as we will see in Section 2.4, solving the eigensystem ofthe large and sparse Jacobian matrices even once accounts for most of the computationalcost of the process. This makes any algorithm that relies on the calculation of theeigensystem multiple times slower and impractical. The other alternative for measuringthe gradients of the eigenvalues is using reverse automatic differentiation. Automaticdifferentiation, which is used extensively in computer science, measures the gradients ofan output following the natural flow of its calculations using chain rule differentiation.This simply means following the chain rule to measure the gradients of a dependentvariable by multiplying the gradients of the input and intermediary variables. One simpleway to achieve this is to replace each primitive operation in the original program by itsdifferential counterpart. To accommodate both of these operations in the same program,real variables are augmented by the well-known dual numbers. This transforms a program46that evaluates an expression to an overloaded program (in programming languages thatsupport overloading) that can also evaluate the derivative of the expression at the sametime. The forward-mode automatic differentiation comes with a big disadvantage whichis re-evaluation of the entire process for every output variable. This would be verycostly if we wanted to calculate the gradient of a large complicated expression of manyvariables, which happens in this case as well. On the other hand, one can use reverseautomatic differentiation where the input-output roles of variables are inverted. Thismeans that instead of using the chain rule to start from an input variable and measuringthe gradient of the output variable, the operation is reversed to drive the gradients of theinput variables from the output variable. This operation can be very helpful in scenarioswhere only one output variable which depends on many input variables is to be sought,as all the gradients of all the input variables can be solved in one go. These gradientscan be reconstructed to find the gradients of the output variable. However, one problemis that the derivative calculations can no longer be easily blended with the evaluations ofthe original expression, since all the derivative calculations appear to be going in reverseof the original program. Another main challenge with automatic differentiation is thatit is quite intrusive to programs and in some cases it needs to be added to the programfrom the ground-up. This makes it very expensive to use in existing applications. Forthese reasons, in this thesis I use the perturbation theory and direct application of finitedifferences on the Jacobian matrix to derive the gradients of the eigenvalues with respectto changes in the mesh.If the matrix eigenvalue problem of interest for the Jacobian matrix of the semi-discretized partial differential equation is:J(−→u ,−→ξ )Ri(−→ξ ) = λiRi(−→ξ ) (2.12)where Ri is the ith right eigenvector associated with the ith eigenvalue λi, and the JacobianJ(−→u ,−→ξ ) is a function of both solution and the mesh, then the eigenvalue derivatives with47respect to the mesh coordinate vector ξ are obtained in summary as follows:ddξ(JRi = λiRi) (2.13)Li(dJdξRi + JdRidξ= dλidξLi + λidRidξ) (2.14)(LiJ − Liλi) ∂Ri∂ξ= Li∂λi∂ξRi − Li∂J∂ξRi (2.15)where multiplication of Eq. 2.14 by the left (row) eigenvector Li leads to Eq. 2.15where the left hand side is zero by definition of the left eigenvectors as shown in theprevious section. Moreover, the eigenvectors are normalized so that Li · Ri = 1. Hencethe gradients of the eigenvalues with respect to the mesh vertices can be calculated using:dλidξ= LidJdξRi (2.16)To measure the gradients of the eigenvalues in Eq. 2.16, one needs to first evaluate thederivatives of the Jacobian matrix. As explained previously in Section 1.1.5, the globalJacobian matrix itself is calculated directly using Eq. 1.18. This analytical flux Jacobianis a function of the numerical solution, mesh properties, flux function, and discretizationmethod. For a specific flux function and discretization scheme, the derivative of theJacobian matrix with respect to the mesh entities in Eq. 2.16 is written as:dJ(−→u ,−→ξ )dξ= ∂J∂u∂u∂ξ+ ∂J∂ξ(2.17)where ∂u∂ξcan be found from the derivative of the residual (Eq. 1.8) at steady state:R(−→u ,−→ξ ) = 0d<dξ= J ∂u∂ξ+ ∂<∂ξ= 0(2.18)Substituting for ∂u∂ξfrom Eq. 2.18 in Eq. 2.17 results in:48dJ(−→u ,−→ξ )dξ= −∂J∂u(J−1) ∂<∂ξ+ ∂J∂ξ(2.19)However in the scope of this thesis, I assume that the first term is negligible comparedto the second term in Eq. 2.19 (this assumption is tested in Chapter 5) since I am usinga close approximation of the final solution to approximate the Jacobian of the unstableproblem at the steady state. Herein, I am using three different solutions to approximatethe Jacobian of an unstable problem: first, a lower order steady state solution for thehigher order problem; second, a surrogate solution of a low fidelity problem such aspotential flows; third, a half converged solution of the same order problem. Additionally,it is noteworthy to mention that for linear problems such as advection, where the Jacobianis independent of the solution, the first term in Eq. 2.19 is exactly zero. In either case,I approximate the derivative of the Jacobian matrix J with respect to the mesh entitiesusing partial finite differences:∂J∂ξ=J(−→ξ + δ−→ξ)− J(−→ξ)∥∥∥δ−→ξ ∥∥∥ (2.20)Hence, using equation 2.16, I am able to predict changes induced in eigenvalues bymesh movement. Figure 2.1 shows that for a good range of δξ, the gradient of theJacobian more or less does not change for a specific mesh location. The horizontal axisin Figure 2.1 is the size of δξ in Eq. 2.20 where the length scale is the length of theshortest edge incident on each vertex. This analysis validates the use of finite differencesto calculate the gradients of the Jacobian matrix.2.1.3 Optimization SchemeUsing energy stability results along with the knowledge of eigenvalue derivatives upon anymesh perturbation, I can tune the mesh movement so that the real parts of eigenvalues(especially the unstable ones) decrease. However, this is not an easy task, as naively49Finite difference size(% of length scale)L 2 norm of gradients10-8 10-6 10-4 10-2 10010001500200010-9 10-8 10-7 10-6 10-52073.22073.222073.242073.262073.282073.3Figure 2.1: Sensitivity map of the gradient of the Jacobian matrix with respect to per-turbation size for an inviscid Burgers’ problem50changing the mesh to improve one eigenvalue may lead to destabilizing other (stable)eigenvalues.One intuitive way to change the mesh is to consider all (right-most) eigenvalues sep-arately. Surely, the easiest route to stabilizing a single eigenvalue regardless of the othereigenvalues is to change the mesh in exactly the opposite direction of the gradient of theeigenvalue (steepest descent method) which means that the following inequality shouldhold:<{λnew} ≈ <{λorig}+4−→ξ · ∂λ∂ξ≤ 0 (2.21)where <{.} is the real part of the eigenvalue. This results in a mesh movement vectorwith the direction and size of:4ξ = −|k|<{∂λ∂ξ}With k > <{λ (ξorig)}(∂λ∂ξ)2 (2.22)A complication arises when there are multiple unstable or nearly unstable eigenvaluesbecause their movement vectors could partly or completely contradict each other. Oneway to solve this is to take a weighted average of these multiple vectors with weightsproportional to how positive (unstable) the corresponding eigenvalues are to obtain asingle mesh movement vector. Another more sophisticated approach is to reformulatethe problem to stabilize all the unstable eigenvalues collectively. To do so, I negate thereal part of the right-most eigenvalues directly. The single movement vector, −→d = 4ξ,should satisfy all the linear inequalities (Eq. 2.23) required to stabilize the problem:<{λj}+−→d · ∂λj∂ξ≤ 0 1 ≤ j ≤M (2.23)where M is the number of unstable eigenvalues. The optimization variables are theentities of the mesh movement vector, the linear optimization problem is defined as:51minM∑jsj where sj =(<{λj}+−→d · ∂λj∂ξ)(2.24)where sj are the negative of the slack variables (positivity of each inequality), subjectto linear constraints sj ≤ 0. Since we have a linear optimization problem, the optimumsolution to the summation of the objective functions is equivalent to the solution of themulti-objective minimization of each eigenvalue. In other words, instead of minimizingthe slack variable for each eigenvalue separately, I minimize the summation of the slackvariables. The upper bounds for the optimization variables are based on the local lengthscale to avoid any non-conformality or irregularity in the mesh after the modification.This empirical criterion limits all vertex displacements to be less than 10% of the lengthof the shortest incident edge. However, in a rare instance of creating a negative volume,the update to the mesh is rolled back to smaller scaled displacement vectors at vertices.This second test ensures the positivity of all cell sizes.2.1.4 Vertices to MoveThe key to this stability analysis is to approximate the gradients of the eigenvalues withrespect to the mesh vertices ∂λj∂ξi. However, calculation of the gradients of the eigenvaluesfor all degrees of freedom in the mesh is highly time consuming and expensive. Yet,there is no need to calculate this for the entire mesh as generally only part of the meshis responsible for instabilities. In practice, moving most parts of the mesh results innegligible changes in the least stable eigenvalues, which in turn implies that these specificeigenvalues are independent of the mesh in that region.The right eigenvector is a mode of the solution. Therefore if a Jacobian matrix tendsto have an unstable solution, the largest components in the right eigenvectors of theunstable modes will specify the parts of the mesh where things have gone wrong astheir effects in these parts are magnified. The gradients and hence movement vectors52CVEigenvector Component0 50 100 150 200-1-0.8-0.6-0.4-0.200.20.40.60.811st = 5.0262st = 1.5773rd = 0.2297(a) Span of the right eigenvectors for an in-viscid Burgers’ problem(b) Vertices on the control volumes of theJacobian fill of the control volume withthe largest eigenvector componentFigure 2.2: How to choose vertices for perturbation in finding the derivative of an eigen-valuecalculated in this way do not include all degrees of freedom, yet they represent the partsof the mesh responsible for degeneration of the stability. Moreover, there is no need forthe exact calculation of the gradient of the eigenvalues, as any approximate one is able toguide the mesh modification in the right direction for better stability properties. Henceto approximate ∂λj∂ξione needs to1. List all the right eigenvector components (e.g., see Figure 2.2a)2. Select the largest components of the eigenvector3. List all CVs corresponding to these components as well as the ones in their Jacobianfill (e.g., see Figure 2.2b) as these ones also affect the value of the flux integral inthe selected control volume.4. Perturb vertices located on these CVsThe steps mentioned above make the algorithm much faster, as the gradients are53only calculated for a limited number of vertices. Also, only certain rows of the Jacobianmatrix need to be calculated in this way which is explained in detail at the end of thischapter. The next sub-section explains what happens if the vertex which is selected forthe perturbation and optimization lies on a boundary.2.1.5 Allowing Boundary MovementAs explained in Subsection 2.1.4, the mesh movement vector only includes a small subsetof control volumes found from the rightmost eigenvectors. However, there are manyinstances that the dominant control volumes reside at the boundary which result in higheigenvalue gradient adjacent to the wall (e.g. Figure 2.3). Figure 2.3 shows that thecomponents of the gradient of the rightmost eigenvalue are largest near the leading edgeof a NACA 0012 airfoil for an inviscid Euler problem with M = 0.8 and α = 2◦ while therest of the physical domain has significantly smaller gradients. Hence, moving only theinterior vertices may not be enough or be the fastest route to stabilize the problem. Forthese reasons one needs to move the boundary vertices along the boundary curve as well.This constraint can be either imposed in the optimization process or as a post processingstep. The latter is easier to implement and it does not stiffen the optimization process.In this way the movement of vertices at the boundaries is done in two steps. First theoriginal vertex is moved by the unconstrained perturbation vector. Second, this interimpoint is projected to the closest point on the boundary curve. The final movement vectoris from the original point to the new projected point on the curve. Figure 2.4 shows thesesteps schematically.2.2 A Priori Stability AnalysisThe key to this stability analysis and improvement is calculating the gradients of theeigenvalues using Eq. 2.16 and 2.20 for a specific Jacobian matrix and eigensystem.54Gradients5.0E+035.0E+025.0E+015.0E+005.0E-015.0E-02Figure 2.3: Largest eigenvalue gradient components of the rightmost eigenvalue for NACA0012, inviscid, Ma=0.8, α=2◦55Figure 2.4: Boundary vertex movement56However when a problem has numerical stability issues, there is no steady state solutionavailable. This means that one has to predict and improve the stability of a numericalmethod before any attempt at solving it. The lack of solution in these cases is the mainbarrier for using this stabilization algorithm. To rectify this issue, surrogate solutions areused in the absence of the steady state solution to approximate the Jacobian matrix andaccordingly the eigenspace of the final steady state. This surrogate solution can be anysolution that can approximate the final steady state solution of the numerical method.This family of solutions as explained in this section includes a lower order solution fromthe same physical problem, or using a low fidelity solution of another physical problemwhich can roughly mimic the steady state solution. This method of using a surrogatesolution is called a priori stability analysis, as the stability of the problem is sought andimproved before any attempt at solving the primary expensive problem.One thing to be wary of when analyzing the results of this stabilization method andjudging whether it is working, is the use of advanced time stepping schemes. To evolve tothe steady state solution from the spatially discretized Eq. 1.8, one can normally use theimplicit Euler time advance or Newton solvers (see [17, 71, 72] and references therein) forfaster convergence. As explained in Section 1.1.6, in ANSLab the main steady state solveris such a scheme that invokes local time-stepping, line search, and other sophisticatedmethods to improve numerical stability. However these techniques for accelarating thesteady state convergence are empirically shown to mask the instability for the semi-discrete system due to slightly positive eigenvalues. Moreover, the stability region ofthe implicit time advance scheme often includes much of the right side of the spectralmap of the semi-discrete Jacobian which in turn makes these schemes unsuitable for thepurpose of this chapter by hiding the positive eigenvalues of the semi-discrete Jacobian.For these reasons, when instabilities are being masked in this thesis, the time advancescheme is switched to Crank–Nicolson time integration which is an A-stable scheme andensures that a positive eigenvalue leads to growth of the solution, hence instability. This57change in time advance scheme shows that the present methodology for stabilization isworking. However, if not stated otherwise, the main time advance scheme explained inSection 1.1.6 is being used.2.2.1 Lower Order SolutionAs a first approach to stabilizing finite volume methods on unstructured meshes, a lowerorder solution is used to approximate the Jacobian of a higher order unstable problem.For example, when the spatial discretization is second order and the Jacobian matrix ofsecond order flux residual is to be sought, a first order steady state solution to the samephysical problem on the same mesh is used. To show the applicability of the stabilizationalgorithm, I present the stability results for three different physical problems: the linearadvection problem in a 3D rectangular channel, the inviscid Burgers problem in a 2Drectangular domain, and the Euler problem with the Roe flux scheme for compressibleflows around NACA 0012, and a multi-element airfoil with varying resolutions.2.2.1.1 Mesh ImprovementAs a preliminary test case, a 3D second order advection problem for a rectangular channelwith five hundred tetrahedral control volumes has been stabilized by pushing its singleunstable eigenvalue to the left half of the complex plane (see Figure 2.5a). In this problem,only four vertices with maximum displacement of 5% of the local length scale have beenmoved. The stability can also be observed from the convergence history (as is seen fromFigure 2.5b).Note that the Jacobian in the linear advection problem is only a function of meshcoordinates and constant wave speeds and is completely independent of the solution. Thisin turn implies that changing the mesh in a direction predicted by the gradients of theeigenvalues is indeed a proper approach to gain stability. However, for more complicatednonlinear problems, more care must be taken as the Jacobian is also a function of the58ReIm-5 0 5-10-50510BeforeAfter(a) Spectral mapIterationL 2 norm of solution updates5 10 15 20 25 30 35 40 45 5010-510-410-310-210-1100101102103UnstableStable(b) Convergence historyFigure 2.5: Before and after the mesh movement for a 2nd order advection problem ina 3D rectangular channel with five hundred control volumes. Time stepping is done byCrank–Nicolson.solution.Another complication arises when there are multiple unstable or nearly unstable eigen-values; thereby there are multiple candidate vectors which could partly or completelycontradict each other. To mitigate this problem, an optimization process as explained inSubsection 2.1.3 is used to obtain a single movement vector that can collectively improveall the unstable eigenvalues. Secondly, I chose to put more emphasis on the rightmosteigenvalues, so as to disregard an eigenvalue in calculating the movement vector, if itwas contradicting (i.e., has a negative dot product with) the movement vector calculatedsolely from the right-most ones.The first trial for non-linear problems is done using the inviscid Burgers’ problemwhere the Jacobian of the semi-discrete system is no longer independent of the solution.To linearize the second order Jacobian, a first order solution is used to approximate thesecond order Jacobian, which specifies the unstable eigenvalues as well as parts of themesh responsible for these instabilities. Figure 2.7 shows how modifying the mesh locallyhas changed the unstable eigenvalues in Figure 2.6a to the left half of the complex plane59ReIm-80 -60 -40 -20 0-60-40-20204060Original MeshOptimized Mesh(a) Spectral mapIterationsL 2 Norm of Residual20 40 60 80 10010-510-310-1101103105Original MeshOptimized Mesh(b) Convergence historyFigure 2.6: Before and after mesh movement for an inviscid Burgers’ problem on arectangular channel with two thousand control volumes. Time stepping is done by Crank–Nicolson.while the rest of the eigenvalues are more or less the same. The stability effect of themesh movement can also be seen from Figure 2.6b where the slightly changed mesh isconverging as opposed to the original mesh.Likewise in a more complicated problem, an originally unstable subsonic flow on acoarse mesh with nearly five hundred control volumes around the NACA 0012 airfoilwith Mach number 0.6, and an angle of attack of 5◦ is stabilized. In this problem whichhas two unstable eigenvalues (Figure 2.8a), local changes around the leading edge of theairfoil (Figure 2.9) based on the corresponding eigenvectors are enough to stabilize theproblem.It is noteworthy to mention that the slow convergence in the convergence plots ofFigures 2.6b, and 2.8b has two causes: first, the rightmost eigenvalue of the perturbedmesh for both Burgers’ and Euler problem is quite close to zero which makes the problemill-conditioned. Secondly, the presence of conjugate pairs of eigenvalues in the rightmostregions of the spectral map causes decaying spurious oscillations. The implicit timeadvance scheme with local time stepping and a line search algorithm converges much60Figure 2.7: The snippet of mesh modification to stabilize the problem for an inviscidBurgers’ problem in a rectangular domain with two thousand control volumes. The solidlines are the original mesh and the dashed red lines are the mesh after optimization.ReIm-80 -60 -40 -20 0 20-200-150-100-50050100150200Original MeshOptimized Mesh-0.1 -0.05 0.05 0.1-101(a) Spectral MapIterationsL 2 Norm of Residual20 40 60 80 10010-510-310-1101103105 Original MeshOptimized Mesh(b) Convergence HistoryFigure 2.8: Before and after mesh optimization for an Euler problem with Mach number0.6, and angle of attack of 5◦ around NACA 0012 airfoil with five hundred control volumes.Time stepping is done by Crank–Nicolson.61Figure 2.9: Mesh modification to stabilize the problem for an Euler problem with Machnumber 0.6, and angle of attack of 5◦ around NACA 0012 airfoil with five hundred controlvolumes. The solid lines are the original mesh and the dashed red lines are the meshafter optimization.62IterationsL 2 Nrom of Residual100 101 10210-710-510-310-1101Original Mesh, max=0.42Optimized Mesh, max=-0.0102(a) Convergence history (b) Before and after mesh optimization, the solidlines are the original mesh and the dashed red linesare the mesh after optimization.Figure 2.10: Before and after mesh optimization for an Euler problem with Mach number0.6, and angle of attack of 5◦ around NACA 0012 airfoil with one hundred thousandcontrol volumes. Time stepping is done by implicit Euler time advance.faster in general.To further validate the applicability of our stabilizing algorithm, two numerical ex-periments for Euler problem on finer meshes for NACA 0012 airfoil and a multi-elementairfoil are conducted. Figure 2.10 shows that even for a fine mesh with nearly one hundredthousand control volumes, optimizing only a small number of vertices (28 with maximumdisplacement of 10% of the local length scale) in problematic parts of the mesh is suffi-cient to stabilize the problem. Likewise, Figure 2.11 illustrates that the same algorithmis applicable to more complicated geometries with sharp edges as well. In this case, 46vertices have been moved with a maximum displacement of 6% of the local length scaleto stabilize the problem.2.2.1.2 With Boundary Vertex MovementAs explained in Subsection 2.1.5, for problems for which the mesh related instability issuesare higher near the boundaries (see Figure 2.3), meaning that the eigenvalue gradients63IterationsL 2 Nrom of Residual100 101 10210-710-510-310-1101Original Mesh, max=3.085Optimized Mesh, max=-0.022(a) Convergence history (b) Before and after mesh optimization, the solidlines are the original mesh and the dashed red linesare the mesh after optimization.Figure 2.11: Before and after mesh optimization for an Euler problem with Mach number0.4, and angle of attack of 1◦ around multi element airfoil with twenty five thousandcontrol volumes. Time stepping is done by implicit Euler time advance.are much larger nearby the boundaries, allowing the boundary vertices to move is ofparamount importance. For example, for a transonic inviscid Euler problem with Roeflux and Mach number 0.8 and angle of attack of 2◦ which is originally unstable, onecan change the mesh in three different ways (snippets of the updated meshes are shownin Figure 2.12): moving all vertices, moving only the interior vertices, and moving onlythe boundary vertices. Figure 2.13 shows the optimization process of the rightmosteigenvalue where each point represents the real part of the corresponding eigenvalue afteran iteration of the mesh perturbation. At each step, the Jacobian of the modified mesh iscalculated by the solution of the initial mesh. The spikes in the curves illustrate the useof the new converged solutions, calculated on the new modified mesh, in approximatingthe Jacobian on the same mesh. Eventually, the curves in Figure 2.13 fall below the zeroline which suggest that the Jacobian of the final mesh has negative rightmost eigenvalueby either solution. As is seen from Figure 2.13, not being able to move the boundaryvertices inhibits or slows down the optimization process, since the important boundary64vertices are not amongst the subset of vertices to be optimized. Moreover, it can alsochange the convergence rate of the final stable problem (Figure 2.2.1.2). The reason thatthe residual of the unstable case in Figure 2.2.1.2 does not increase in the last iterationsis that the implicit time advance scheme keeps shrinking down the time step size tocompensate for the unfavorable residual in the previous step.2.2.2 Surrogate Solution: Panel MethodsIn the first attempt to use a low fidelity solution as the surrogate solution for the mainphysical problem of interest, I investigated a prior stability analysis of subsonic flowswith an estimated solution of the Euler problem using panel methods. Panel methods,explained in Section 1.4, are linearized potential flows, where the potential solution any-where in the domain is approximated by solving a much smaller system of equationsfor source and doublet singularities at the boundary. As is seen from Figures 2.15a and2.15b, even though the panel solution does not completely agree with the Euler solution,it is still be sufficient to detect the mesh related instabilities for subsonic flows. Figures2.2.2 and 2.17 show that the mesh optimization based on the panel solution was able tostabilize an originally unstable problem. In this subsonic problem with Mach number0.55 and angle of attack of 3 degrees, the rightmost eigenvalue of the semi discretizedJacobian, approximated by panel solution, was λ0unstable = 0.285. This predicted positiverightmost eigenvalue means that the corresponding Euler problem would be unstable,which in fact it is (see Figure 2.2.2). Yet, by optimizing the mesh according to thispredicted eigenvalue and its gradients with respect to the mesh vertices, the problem wasstabilized.2.2.3 Surrogate Solution: One-equation Potential FlowIn another attempt to a priori stabilization, I am using a medium fidelity solution ob-tained by solving the one-equation full potential flow (see Section 1.4) to approximate65Perturbed MeshOriginal Mesh(a) No ConstraintsPerturbed MeshOriginal Mesh(b) Only Interior VerticesOriginal MeshPerturbed Mesh(c) Only Boundary VerticesFigure 2.12: Optimized meshes for NACA 0012, Inviscid, Ma=0.8, α=2◦66Stabilizing IterationsRightmost Eigenvalue50 100 15000.511.5No ConstraintOnly Bdry VertsOnly Inter Verts(a) Entire optimization historyStabilizing IterationsRightmost Eigenvalue5 10 15 20 25 30 35 4000.511.5No ConstraintOnly Bdry VertsOnly Inter Verts(b) Close-upFigure 2.13: Eigenvalue optimization for NACA 0012, Inviscid, Ma=0.8, α=2◦; every dotrepresents an iteration of the mesh optimization67IterationsL 2 Norm of Residual10 20 30 40 50 60 7010-810-710-610-510-410-310-210-1100101Unperturbed Mesh ( max=1.822)All Verts Free To Move ( max=-0.028)Only Bdry Verts Perturbed ( max=-0.030)Only Inter Verts Perturbed ( max=-0.005)Figure 2.14: Convergence rate for NACA 0012, Inviscid, Ma=0.8, α=2◦, Time steppingis done by backward nonlinear solver.68X/LMach0 0.2 0.4 0.6 0.80.20.40.60.81O(h2) Euler, upper surfaceO(h2) Euler, lower surfacePanel, upper surfacePanel, lower surface(a) Mach numberX/LCp0 0.2 0.4 0.6 0.8 1-1.5-1-0.500.51O(h2) Euler, upper surfaceO(h2) Euler, lower surfacePanel, upper surfacePanel, lower surface(b) Pressure coefficientFigure 2.15: Comparison of second order Euler and Panel solutions for a NACA 0012airfoil with 10000 control volumes and free stream Mach number of 0.55 and angle ofattack of 3◦. The second order Euler solution has been obtained after stabilization.69IterationsL 2 Norm of Residual20 40 60 80 10010-710-510-310-1101103 Optimized MeshOriginal MeshFigure 2.16: Residual history for both unstable and stabilized second order Euler problemon a NACA 0012 airfoil with 10000 control volumes and free stream Mach number of0.55 and angle of attack of 3◦. The panel method for the linearized potential equationhas been used in this case.the Jacobian of the Euler problem. Figure 2.18 shows that using both the one-equationfull potential model and a lower order solution in the second order Jacobian which wasused in Section 2.2.1 can predict the instability of the second order Euler problem. Thepotential solution (or in general any other estimate of the solution), despite being dis-crepant with the lower order solution in calculating the rightmost eigenvalue, can produceperturbation vectors to stabilize the problem without solving the main problem. Figure2.19 shows the perturbed mesh, and it is evident that fixing a part of the mesh locallycan push the unstable eigenvalue to the left half of complex plane. Moreover, as pre-dicted by the eigenvalue analysis of the second order Jacobian with potential solution,Figure 2.20 shows that the second order problem is in fact unstable on the original mesh.Yet, by applying the mesh perturbations, I was able to stabilize the unstable problemfor the same second order operator before solving the Euler problem. Consequently onecan conclude that in this algorithm, any rough estimates of the solution is sufficient topredict and fix the stability a priori.70Figure 2.17: Original and optimized meshes for NACA 0012 airfoil with 10000 controlvolumes and free stream Mach number of 0.55 and angle of attack of 3◦. The panelmethod for the linearized potential equation has been used in this case. The black linesare the original mesh and the red lines are the optimized mesh.71ReIm-4 -2 0 2 4-1-0.500.51original: O(2nd) Jac, O(1st) solutionoriginal: O(2nd) Jac, potential solutionperturbed: O(2nd) Jac, O(2nd) solution-0.04 -0.03 -0.02-0-000.0.Figure 2.18: Eigenvalues before and after mesh perturbation, NACA 0012, Inviscid, Machnumber of 0.4, α = 5◦. The potential solution and first order solution are used to predictthe eigenvalues of the second order discretized Jacobian. The one-equation potentialmodel has been used in this case.2.2.4 Surrogate Solution: Two-equation Potential FlowTo use the two equation potential model, described in Section 1.4, for cell-centered finitevolume methods I needed to devise a new way to satisfy the wake boundary conditions. In72Original meshPerturbed meshFigure 2.19: Perturbed Mesh, NACA 0012, Inviscid, Mach number of 0.4, α = 5◦. Theone-equation potential model has been used in this case.IterationsL 2 Norm of the residual5 10 15 2010-810-710-610-510-410-310-210-1100101Original MeshPerturbed MeshFigure 2.20: Convergence history, NACA 0012, inviscid, Mach number of 0.4, α = 5degrees. The one-equation potential model has been used in this case.73Figure 2.21: The wake cut for NACA 0012 airfoil, Inviscid, Mach 0f 0.72, and angle ofattack of α = 5°, and with 6200 control volumesboth finite element and the vertex-centered finite volume methods, boundary conditionson the wake (Eq. 1.87 and 1.88) are imposed on the vertices on the wake. However, incell-centered finite volume methods, finding vertices on the wake to apply wake boundaryconditions is not always possible. Hence, in my implementation of boundary conditionsfor the wake, I try to apply these conditions on the control volume average values inthe flux calculations. Therefore, I first find a surface in the mesh starting from thetrailing edge along the free stream continuing on to the far field boundary to create adiscontinuity within the potential flow (e.g., see Figure 2.21). This surface is createdfrom the trailing edge of the airfoil by following along the edges that are more alignedto the direction of the free stream and connecting vertices at the end of these edge toreach to the farfield. Then, control volumes on both sides of this surface are taggedas the wake control volumes. Afterwards the tangential velocities and densities withinthese control volumes are averaged (Eq. 2.25) to weakly enforce the continuity of thetangential velocities and densities. It is noteworthy to mention that the continuity of thenormal velocity is automatically obtained by conservation of the mass.(ρ, V T)lower=(ρ, V T)upper=(ρu + ρl2 ,V Tu + V Tl2)(2.25)Finally in this last attempt to a priori stability improvement, I was also able to showthe same stability effects on a slightly transonic flow with Mach number of 0.72, and an74angle of attack of 5 degrees. In this problem, a two-equation model of potential flow (asexplained in Section 1.4) was solved to obtain a surrogate solution for the second orderEuler problem which was originally unstable with one unstable eigenvalue on a mesh of6200 control volumes around NACA 0012 airfoil. Figures 2.22 and 2.23 show the meshand the residual history of the steady state solutions for unstable and stabilized problemsrespectively. By comparing the solutions (Figures 2.24a and 2.24b) of the second orderEuler problem (after stabilization) and the potential flow which is used as its surrogatesolution to stabilize the Euler problem, it is observed that even a solution which is notcompletely accurate can be used to approximate the steady state discretized Jacobianof an unstable problem. Hence, the mesh can be optimized by doing eigenanalysis onthis approximate Jacobian before attempting to solve the unstable problem. Moreover,another advantage of having a cheap surrogate solution is to use it to initialize the mainproblem (in this case the Euler problem). In this way, many of the main features of thesteady state solution have already been resolved, and since the residuals are much smallerat start-up, the time advance scheme can be used much more aggressively. In an extremecase, a direct Newton solve could be applied to this solution to get to the steady state inone or two steps. Figure 2.25 showcases the advantage of using the surrogate (potential)solution as an initial solution over using the free stream or far field conditions. It isalso note worthy to mention that in this case both the first and the second order Eulerproblems were unstable. Hence, initializing the second order Euler problem with the firstorder (which itself was unstable) could not expedite the convergence. The convergenceplots in Figure 2.25 were obtained on the original mesh in this case and without meshoptimization.75Original Mesh (Black)Optimized Mesh (Red)Figure 2.22: Original and optimized meshes for NACA 0012 airfoil, Inviscid, Mach 0.72,and α = 5 degrees, with 6200 control volumes. The two-equation potential model hasbeen used in this case.IterationsL 2 Norm of Residual20 40 60 80 10010-710-510-310-1101103Optimized MeshOriginal MeshFigure 2.23: L2 norm of residual for original and optimized meshes for NACA 0012 airfoil,Inviscid, Mach 0.72, and α = 2° with 6200 control volumes. The two-equation potentialmodel has been used in this case.76X/LCp0 0.2 0.4 0.6 0.8-1.5-1-0.500.51O(h2) Euler, upper surfaceO(h2) Euler, lower surfacePotential, upper surfacePotential, lower surface(a) Pressure coefficientX/LMach0 0.2 0.4 0.6 0.80.60.811.21.4O(h2) Euler, upper surfaceO(h2) Euler, lower surfacePotential, upper surfacePotential, lower surface(b) Mach numberFigure 2.24: Potential and Euler solutions on a mesh of 6200 control volumes for a NACA0012 airfoil with Mach number of 0.72, and an angle of attack of 5◦. The Euler solutionis obtained after it has been stabilized. The two-equation potential model has been usedin this case.77IterationsL 2 Norm of Residual20 40 60 80 10010-710-510-310-1101103PotentialO(h2) Euler initialized by potentialO(h2) Euler initialized by free streamO(h2) Euler initialized by O(h1) EulerO(h1) Euler initialized by free streamFigure 2.25: Residual history of Euler problem initialized with potential flow solutionand far field free stream solution for NACA 0012 airfoil, Mach number 0.72, and angle ofattack of α = 5°, and with 6200 control volumes. Both O (h1) and O (h2) Euler problemswere unstable in this case. The two-equation potential model has been used in this case.782.3 A posteriori Stability AnalysisOne concern that might arise from the aforementioned calculations is that predictingand improving the stability of the problem of interest by first solving another problem(a lower order problem or another physical problem) is not practical. To alleviate thisconcern and also show the flexibility of this stabilization algorithm, I use half convergedsolutions to approximate the Jacobian at the steady state when the convergence stallsor slows down. This use of the half-converged solutions for stabilizing the numericalmethod is called a posteriori stability analysis as one has to first encounter the instabilityduring the run time and then invoke this method to resolve the instability issues. In thisway, there is no need for first solving another physical or numerical problem to helpstabilizing the primary unstable numerical method. This indeed reduces the user’s effortfor implementing this stabilization scheme.2.3.1 Stability Analysis with Half Converged SolutionThis section describes one of the most prevalent scenarios that happen when numeri-cally solving partial differential equations which is the case when the problem does notfully converge or convergence plateaus to a half converged solution. To improve uponthe steady state stability of these cases, these partially converged solutions are used toapproximate the Jacobian matrix. This Jacobian is in turn used in Eq. 2.20 with smallmesh perturbations to approximate the derivative of the Jacobian. However, limitationsof this approach are discussed in Chapter 5.1 where the sensitivity of the eigenvalues tothe surrogate solutions (in this case half-converged solution) are investigated.Two cases are investigated in this section: a third order inviscid Euler problem withMach number of 0.68 and angle of attack of 3.5 degrees, and a third order laminar flowwith Mach number of 0.68 and angle of attack of 3.5 degrees and Reynolds number of5000. Both cases are solved on two different unstructured meshes of nearly 5000 control79volumes. Originally, these problems both converge only partially as seen in Figures2.26a and 2.26b. The inviscid and the laminar problems on the original meshes haverespectively positive rightmost eigenvalues of λ0 = 0.0041 and λ0 = 0.0028 while after thestabilization, their right most eigenvalues are improved to λ0 = −0.016 and λ0 = −0.021respectively. To stabilize these problems, the partially converged solutions are used toapproximate the Jacobian of their spatially discretized problems. The eigenvalues ofthe approximate Jacobian and their derivatives are used to find mesh movement vectorsrequired to project these eigenvalues to the stable region. From Figures 2.26a and 2.26b,it can be seen that after the stabilization, both problems were able to break the plateauand converge rapidly to steady state solutions.Figures 2.27 and 2.28 show how locally moving a limited number of vertices basedon the optimization scheme described in Section 2.1.3 suffices to stabilize the previouslyunstable problems.As mentioned earlier, to approximate the Jacobian matrix, the partially convergedsolutions where the lines bifurcate in Figure 2.26 are used. To see how closely thesehalf-converged solutions approximate the steady-state solution, I examine the coefficientof pressure around the airfoil. From Figure 2.29, it can be observed that even though thehalf converged solutions are not satisfying the discretized partial differential equationswhere the residual of the flux integrals are still high, they are more or less capturingthe flow around the airfoils. It can be deduced that for many of the engineering caseswhere the convergence slows down or stops before fully converging, many of the physicalproperties of the solution are already converged, thus, the solution at these half convergedstates are valuable for troubleshooting the stability of the rest of the path to the steadystate.80IterationsL 2 Norm of Residual20 40 60 80 10010-710-510-310-1101103OriginalStabilized(a) Third order Euler problem with M = 0.68, α =3.5◦ around a 0012 airfoil on an unstructured meshwith 5000 control volumesIterationsL 2 Norm of Residual20 40 60 80 10010-810-610-410-2100102OriginalStabilized(b) Compressible third order laminar problem withM = 0.68, α = 3.5◦ and Re = 5000 around a 0012airfoil on an unstructured mesh with 5000 control vol-umesFigure 2.26: Convergence history for two unstable and stabilized problems. The solutionwhere the line bifurcates is used as the surrogate solution for stabilization.81Figure 2.27: Mesh before and after stabilization of a third order Euler problem withM = 0.68, α = 3.5◦ around a 0012 airfoil on an unstructured mesh with 5000 controlvolumes. Red lines are the changes in the mesh after optimization.Figure 2.28: Mesh before and after stabilization of a third order Laminar problem withM = 0.68, α = 3.5◦ and Re = 5000 around a 0012 airfoil on an unstructured mesh with5000 control volumes. Red lines are the changes in the mesh after optimization.82X/LC p0.2 0.4 0.6 0.8-1.5-1-0.500.51Half ConvergedStabilized Steady State(a) Third order Euler problem with M = 0.68, α =3.5◦ around a 0012 airfoil on an unstructured meshwith 5000 control volumesX/LC p0.2 0.4 0.6 0.8-1.5-1-0.500.51Half ConvergedStabilized Steady State(b) Compressible third order laminar problem withM = 0.68, α = 3.5◦ and Re = 5000 around a 0012airfoil on an unstructured mesh with 5000 controlvolumesFigure 2.29: Pressure coefficients around 0012 airfoil for two unstable and stabilizedproblems.832.4 Computational Cost and Time StudyThe problem of stabilizing unstable simulations by mesh vertex displacement seems com-putationally expensive at first glance. In fact, the algorithm can be expensive dependingon how long calculating the eigenvalues takes. This step has proved to be a bottleneckof this algorithm where the cost of the eigenvalue solver consumes at least 95% of theoverall time. However, the overall cost and algorithm time can be drastically reduced bytaking advantage of the locality of the changes, sparsity of the matrix, and optimizingthe iterative eigensolver. These steps include:First, at every step of the optimization, only a very small subset of vertices are moved(by selecting the largest components of the eigenvector; see Subsection 2.1.4). This in turnmeans that the gradients of the eigenvalues for only these vertices need to be calculated.The number of such vertices does not change as the problem size grows. From Eq.2.16, the cost of this step is dominated by the growth of right and left eigenvectors asthe problem size increases while the number of non-zero rows in ∂J∂ξremains the same.Hence, only limited vector-vector products are required.Second, in calculating ∂J∂ξfor a vertex only a small number of rows of the Jacobianmatrix are recalculated and are taken into account. This is because moving every vertexonly affects the flux residual for control volumes that either have this vertex directly orhave it in the reconstruction stencil of control volumes that appear in the flux calculationof this CV. This information is also readily available from the sparsity pattern of theJacobian matrix. Non-zero entities in each column of the Jacobian matrix indicate thecontrol volumes which appear in the flux residual of the control volume corresponding tothat column, hence no extra work is required. Most importantly, in this way the numberof flux residual reevaluation is completely independent of the size of the problem andis only a function of how large the reconstruction stencil is. This means that this costremains roughly constant as the size of the mesh grows for a required order of accuracy.It is also a small fraction of one Newton iteration of the flow solver.84Third, the most expensive part of the algorithm is calculating the right-most eigenval-ues of the large and sparse Jacobian matrix. However, in our experience, most unstableproblems only have a very small and limited number of positive eigenvalues. In this way,the iterative eigensolver only needs to calculate a small portion of the eigen-spectrumwhere multiple techniques such as filtering unnecessary regions of the spectrum (e.g.,the whole left half plane), spectral transformations (where eigenvalues are mapped toanother space) and similarity transformations (where the matrix is balanced to havebetter conditioning), and subspace initialization and deflations are available to expeditethe process. For these reasons, we are using SLEPc eigensolver which puts all of thesetechniques at our disposal. According to the SLEPc documentation [53], the overall costof an iterative eigensolver with a matrix size of n and a number of requested eigenvaluesm is governed by:• storing a basis of the subspace that is O(m) (recommended to be 2m) vectors oflength n,• orthogonalization of the basis vectors with a cost of roughly O(n ·m2),• constructing a projected eigen-problem of size O(m),• solving a dense matrix of size O(m3)It is clear that the overall cost is mostly dominated by the number of the requestedeigenvalues. Hence, solving for a small number of the rightmost eigenvalues is not anoverly expensive process and can scale linearly with the size of the Jacobian matrix.Yet, this step takes up almost all of the computational time of the algorithm. To betterillustrate the computational cost and scalability of the algorithm, the detailed time ofthe algorithm to stabilize varying mesh sizes of NACA 0012 airfoil for an Euler problemwith Mach 0.6 and an angle of attack of 5◦ is profiled in Figures 2.30a, and 2.30b. InFigures 2.30a, and 2.30b the overall cost of the flow solver to converge for the stabilized85problem and its first iteration are also included as baseline times for better comparisonof the computational costs.In each iteration of the stabilizing algorithm, the eigenproblem is first solved to obtainthe rightmost eigenvalue and its corresponding right eigenvector where the timings forthis step is presented in Figures 2.30a, and 2.30b. To calculate the left eigenvector ofthe rightmost eigenvalue, the eigenproblem is solved for the transposed Jacobian matrixas well, where the timings for this step are only reflected in the overall time in Figures2.30a, and 2.30b. Fluctuations for the overall time in Figure 2.30a occur mainly since eachmesh is different and therefore different number of iterations are required for each case.Nonetheless for all the cases, the cost of eigenvalue calculation is reduced considerablyafter the first iteration of the algorithm due to subspace initialization of the eigensolverfrom previous iterations.86# CVTime(s)103 104 10510-310-210-1100101102Overall Flow SolverGradients of EigenvaluesSolving for Displacement VectorEigenvalue CalculationOveral Mesh Optimization Algorithm(a) The overall computational time of the stabilizingalgorithm# CVTime(s)103 104 10510-310-210-11001011021st iteration of Flow SolverGradients of EigenvaluesSolving for Displacement VectorEigenvalue CalculationMesh Optimization Algorithm(b) Computational cost of first iteration of the stabi-lizing algorithmFigure 2.30: The computational time and cost of stabilizing algorithm for NACA 0012,Inviscid, Mach 0.6, α=5◦87Chapter 3Solution ReconstructionOptimizationHigh order discretization methods for partial differential equations are shown to improvethe asymptotic rate of convergence of the discrete solution to the continuous solution bydecreasing the error at a faster rate as the mesh is refined (e.g., see [73, 74, 23, 75, 76, 13]),especially for turbulent flows [77, 78] where important physical phenomena in the floware to be resolved. As seen in Chapter 1, this high order discretization relies on thek-exact reconstruction of the solution in each control volume to obtain a more accuratesolution representation (e.g., see [15, 22, 79]) at the flux quadrature points. In the ongoingsearch for better solution reconstruction methods, there have been numerous attempts(e.g., see [80, 81, 82, 83]) to address the accuracy of the solutions or solution smoothnesssuch as ENO/WENO schemes (e.g., see [84, 85, 86, 87, 88, 89, 90, 91]). These schemeshave been successful in solving higher order problems in the presence of discontinuitiesand have improved convergence to steady state in some cases. However, there has beenlittle work to rigorously investigate the impact of the solution reconstruction process onthe stability of the discretization. It is also shown in the literature (e.g., see [92]) thatmuch of the difference in robustness for unstructured grids compared to their structuredcounterparts is due to the reconstruction schemes. Additionally, there is little or no workon schemes that can automatically improve stability and convergence through optimizingthe solution reconstruction. In this chapter, I first expand on the reconstruction mapmethodology introduced by Haider et al. [7] to generalize their theoretical as well as88empirical results beyond second order accurate linear advection. I investigate the effectsof the least-squares reconstruction stencil on the stability of the semi-discretized systemof equations5, seeking the optimum choice of the reconstruction stencil that can improvethe stability and convergence rate without any loss in accuracy.Additionally, I employ the gradients of the eigenvalues of the Jacobian of the semi-discretized finite volume methods with respect to the solution reconstruction parametersto systematically optimize the reconstruction stencil for every control volume in themesh. In this way, the rightmost eigenvalues (the most unstable ones) of the semi-discreteJacobian of the flux integral are stabilized by using the gradients of the eigenvalues withrespect to the reconstruction stencil parameters.The novelty of this chapter is twofold: first, the local reconstruction map analysisis extended beyond second order accurate advection problem. Secondly, a procedure isdescribed to automatically choose the optimum reconstruction stencil for each controlvolume which reduces growth of the quadratic energy of the system, and therefore hasbetter stability properties (i.e. negative semi-definite systems). To the best of my knowl-edge, this is the first scheme of its kind. Additionally, in the long run, I hope that thisresearch paves the way for more research to identify trends in local mesh features as wellas reconstruction stencils that correlate with the stability of the reconstruction.The layout for this chapter is as follows. In Section 3.1, the local reconstruction mapstability analysis is further developed for high order advection as well as the non-linearinviscid Burgers’ and Euler problems. Section 3.2 presents an empirical method, sup-ported with numerical results, to stabilize the PDE discretization by merely choosingnearby control volumes for the reconstruction stencil based on purely geometric proper-ties. Section 3.3 introduces an optimization method which systematically determines thesize and shape of the reconstruction stencil for every control volume in the mesh basedon the gradients of the rightmost eigenvalues of the semi-discretized Jacobian.5That is, the coupled system of ordinary differential equations in time that arises from discretizingthe PDE’s in space.893.1 Local Reconstruction MapOne approach to studying stability finite volume numerical methods which use solutionreconstruction for high order methods is to use the local reconstruction map. The localreconstruction map, which was first introduced by Haider et al. [7], is a set of purelygeometric calculations which ascribes a matrix norm value to every control volume inthe mesh. Haider et al. [7] showed that for a second order discretization of a linearadvection problem, the maximum norm of this local reconstruction map adversely affectsthe numerical stability by contributing to the growth of the quadratic energy. To developthe local reconstruction map analysis for high order and nonlinear schemes in the nextsection, I first briefly describe the terminology and annotations which were developed forthe first and the second order linear advection schemes respectively.3.1.1 1st Order AdvectionFor a linear advection equation with constant velocity c:∂tu(x, t) + c · ∇u(x, t) = 0, (3.1)the flux, f , at the face between cells α, β is:fαβ = (c · nαβ)+Uα + (c · nαβ)−Uβ (3.2)where control volumes β are the first neighbors (the ones with a common face) of thecontrol volume α and nαβ is the normal vector to face αβ. In Eq. 3.2, (c · nαβ)± is thepositive/negative part of the inner product of the wave velocity and the normal vectorto face αβ. With this definition, the semi-discrete form of Eq. (3.1) reduces to:dU(t)dt α= − 1|να|∑βˆΓαβ(fαβ) dΓ (3.3)90Using the quadratic energy equation, defined in sub-section (1.3.1), the energy for thefirst order advection scheme is calculated as follows:ddt(U,MU) = 2∑α(UTαMdUαdt)= 2∑α−UTα∑βˆΓαβ(fαβ) dΓ= −2∑α∑β(UTα [(c · nαβ)+Uα + (c · nαβ)−Uβ])(3.4)To simplify Eq. (3.4), I use two changes of variables and generalizations. As statedearlier, β only iterates over the face neighbors of the control volume α, while α iteratesover all the control volumes in the problem. However, β can be extended to iterate overall the control volumes by redefining the normal vector between two adjacent controlvolumes, nαβ so that:nαβ =⊥αβ if β is a face neighbor of α0 Otherwise(3.5)Additionally, for an upwind scheme the positive part of the dot product of the wavevelocity and the normal vector to the face is identical to the negative of the negative ofthis product:(c · nαβ)+ = −(c · nαβ)− (3.6)Using Eq. (3.6) in Eq. (3.4), and factoring the multiplier, 2, into a summation ofterms, transform Eq. (3.4) to:ddt(U,MU) = −∑α∑β(c ·nαβ)+(|Uα|2 + |Uα|2)+∑α∑β(c ·nαβ)+(UTα Uβ + UTα Uβ)(3.7)91Now using the change of variables as explained in Eq. (3.5) into Eq. (3.7) andextending β to iterate over all the control volumes in the latter, the quadratic energy forthe first order advection scheme is:ddt(U,MU) = −∑α∑β(c · nαβ)+(|Uα|2 + |Uβ|2 − UTα Uβ − UαUTβ)= −∑α∑β[(c · nαβ)+ (Uα − Uβ)2]≤ 0 (3.8)Eq. (3.8) shows that the rate of the growth of energy for the first order advectionscheme is always negative and hence this scheme is always energy stable. However, thisstability does not hold for higher order problems where the solution variables within thecontrol volumes are reconstructed to have more accurate approximation of the flux atthe flux integration points.3.1.2 2nd Order AdvectionFollowing the solution reconstruction explained in Section 1.1.2, the solution anywherewithin the control volume including the Gauss quadrature (GQ) points can be calculatedusing the solution at the cell reference location and the reconstructed solution gradients;the solution within the control volume for a second order scheme is:wα(~x, t) = Uα + σα · (~x− ~xα) (3.9)where the σα terms are the gradient coefficients calculated from the solution reconstruc-tion. The pseudo-inverse, Cαγ, of the reconstruction matrix as explained in Section 1.1.2and 1.1.5 is calculated by using a singular value decomposition. For every control vol-ume, α, the gradient coefficients are then reconstructed using the pseudo-inverse of thereconstruction matrix and the control volume averages in the reconstruction stencil of92control volume α:σα =∑γCαγ (Uγ − Uα) , (3.10)where γ iterates over all the cells in the reconstruction stencil of the cell α. Similar tothe first order, the upwind second order flux function between control volumes α and βfor linear advection problem is defined as:fαβ = (c · nαβ)+wα + (c · nαβ)−wβ, (3.11)where wα and wβ are the linear reconstructed solutions at the Gauss quadrature pointsof control volumes α and β respectively:wα = Uα + σα · (~x− ~xα)wβ = Uβ + σβ · (~x− ~xβ) (3.12)By using the quadratic energy equation, Eq. (1.64), the flux function, Eq. (3.11),and the reconstructed solutions at the Gauss quadrature points, the energy equation forthe second order advection scheme is calculated as follows:93ddt(U,MU) = 2∑α(UTαMdUαdt)= −2∑αUTα∑βˆΓαβ((c · nαβ)+wα + (c · nαβ)−wβ) dΓ= −2∑αUTα∑β(c · nαβ)+Uα + σαˆΓαβ(~x− ~xα) dΓ+−2∑αUTα∑β(c · nαβ)−Uβ + σβˆΓαβ(~x− ~xβ) dΓ (3.13)The gradients of the solution reconstruction, σ, in Eq. (3.13) are replaced by Eq.(3.10). The geometric integration terms at the interface of the control volumes α and βare computed by using the Gauss quadrature rules.kαβ =ˆΓαβ(~x− ~xα) = xαβ − xα (3.14)In this way, kαβ is the distance between the cell reference and the location of theGQ point, xαβ − xα. Using Eq. (3.14) and (3.10) into Eq. (3.13), further simplifies thequadratic energy equation of the second order scheme to:94ddt(U,MU) = −2∑αUTα∑β(c · nαβ)+[Uα +∑γCαγ (Uγ − Uα) kαβ]+−2∑αUTα∑β(c · nαβ)−[Uβ +∑γCβγ (Uγ − Uβ) kβα](3.15)= −2∑α∑β(UTα [(c · nαβ)+Uα + (c · nαβ)−Uβ])+−2∑α∑β(UTα[(c · nαβ)+∑γCαγ (Uγ − Uα) kαβ])+−2∑α∑β(UTα[(c · nαβ)−∑γCβγ (Uγ − Uβ) kβα])(3.16)The first term in Eq. (3.16) is the contribution of the first order terms to the growthof energy which is similar to Eq. (3.4), hence it is always negative. However, the secondterm of Eq. (3.16), which is the growth of the quadratic energy from the second orderterms, is further expanded using the transformations defined in Eq. (3.6) and (3.5):−2∑α∑β(c · nαβ)+{UTα kαβ∑γCαγUγ − UTα kαβ∑γCαγUα +−UTβ kαβ∑γCαγUγ + UTβ kαβ∑γCαγUα}=−2∑α∑β(c · nαβ)+{(UTα − UTβ)kαβ∑γCαγUγ −(UTα − UTβ)kαβ∑γCαγUα}=−2∑α∑β(c · nαβ)+[(UTα − UTβ)kαβ∑γCαγ (Uγ − Uα)](3.17)Now assembling the first order and the second order terms results in the quadraticenergy of the second order advection scheme:95ddt(U,MU) = −∑α∑β[(c · nαβ)+ (Uα − Uβ)2]+−2∑α∑β(c · nαβ)+[(UTα − UTβ)kαβ∑γCαγ (Uγ − Uα)], (3.18)where the first term in Eq. (3.18) is always negative while there is no guarantee for thesecond term which is calculated from the second order terms.Using the quadratic energy (Eq. (1.63)), Haider et al. showed that for a second ordersolution of a scalar linear advection problem, the local reconstruction map for a cell α isdefined as:Rα =∑β∑γkαβCαγ (3.19)where as mentioned earlier, β iterates over the first neighbors (the ones with a commonface) of the cell α and γ iterates over all the cells in the reconstruction stencil of the cellα. In the matrix form, the local reconstruction map for each cell is a lα × nα matrixRα = KαCα, where the Kα is an lα × d matrix, (lα is the number of first neighbors, nαis the number of cells in the stencil of cell α and d is the dimension).KTα = [kαβ1 , kαβ2 , · · · , kαβlα ]T (3.20)Hence, the local reconstruction map norm obtained by multiplication of the pseudo-inverse of the reconstruction matrix by the distance of the Gauss quadrature points andcell reference points shows the sensitivity of Gauss quadrature point solution values tocontrol volume averages for cells in the reconstruction stencil after the linear reconstruc-tion.Haider et al. showed that for a second order advection problem, the quadratic energyis a linear function of the local reconstruction map and its growth can only originate96from the second order terms that include the reconstruction parts. The matrix normused here is the row maximum of the sum of squares of entries in each row, ‖A‖Lα =supx∈lα√ ∑y∈nα|axy|2 for consistency with Haider’s results. In the subsequent subsections theLocal Reconstruction Map stability analysis is extended to high order and nonlinearproblems for the first time.3.1.3 3rd Order AdvectionUsing the same MUSCL6 flux scheme shown in Eq. (3.11) in the advection scheme inEq. (3.1) at the face between cells α, β and a quadratic solution reconstruction at cell α(Eq. (3.9)), the reconstructed solutions at the Gauss quadrature points between cell αand β are:wα = Uα + σα · (~x− ~xα) + 12 (~x− ~xα)T θα · (~x− ~xα)wβ = Uβ + σβ · (~x− ~xβ) + 12 (~x− ~xβ)T θβ · (~x− ~xβ)(3.21)In Eq. (3.21), σ and θ terms are the first order and the second order gradient co-efficients of the solution reconstruction respectively. For example, θα and θβ are thecoefficients of the second order gradients of the solution reconstruction within controlvolumes α and β respectively. Likewise, the pseudo-inverse of the solution reconstruc-tion matrices C(1)αγ and C(2)αγ respectively for the first and second order gradients are definedas:σα =∑γC(1)αγ (Uγ − Uα) (3.22)θα =∑γC(2)αγ (Uγ − Uα) (3.23)6MUSCL stands for Monotonic Upwind Scheme for Conservation Laws. It is a total variation dimin-ishing scheme where reconstructed left and right states are obtained from cell-averaged states at previoustime step. These reconstructed states can then be used to calculate fluxes at the cell boundaries by aRiemann Solver.97Multiplying Eq. (3.3) by Uα and substituting the resulting equations into Eq. (1.64)results in:ddt(U,MU) = dEMdt= 2∑αUTα |να|∑β− 1|να|∑βˆΓαβ(fαβ) dΓ (3.24)Now using the first order and second order gradient coefficients, Eq. (3.22) and (3.23),and substituting Eq. (3.11) and (3.21) for the flux function, fαβ, and the reconstructedsolution respectively into Eq. (3.24) leads to:ddt(U,MU) = −2∑α∑β{(c · nαβ)+UTα Uα + (c · nαβ)−UTα Uβ}+−2∑α∑β{(c · nαβ)+UTα kαβ∑γC(1)αγ (Uγ − Uα)}+−2∑α∑β{(c · nαβ)−UTα kαβ∑γC(1)βγ (Uγ − Uβ)}+−2∑α∑β{(c · nαβ)+UTα k(2)αβ∑γC(2)αγ (Uγ − Uα)}+−2∑α∑β{(c · nαβ)−UTα k(2)αβ∑γC(2)βγ (Uγ − Uβ)}(3.25)The generalization of Eq. (3.5) and the change of variables of Eq. (3.6) helps recastthe first term on the right hand side of Eq. (3.25) which contains only the first orderterms:−2∑α∑β{(c · nαβ)+UTα Uα + (c · nαβ)−UTα Uβ}=−∑α∑β{(c · nαβ)+UTα Uα + (c · nαβ)+UTβ Uβ − 2(c · nαβ)+UTα Uβ}=−∑α∑β(c · nαβ)+|Uα − Uβ|2 (3.26)98Similarly using Eq. (3.5) and (3.6), the third terms on the right hand side of Eq.(3.25) are reformulated to:−2∑α∑β(c · nαβ)−{UTα kαβ∑γC(1)αγUγ − UTα kαβ∑γC(1)βγ Uβ}=2∑α∑β(c · nαβ)+{UTβ kαβ∑γC(1)βγ Uγ − UTβ kαβ∑γC(1)αγUα}(3.27)Now combining the second term on the right hand side of Eq. (3.25) with Eq. (3.27),which represent the second order terms in the discretization, results in:2∑α∑β(c · nαβ)+kαβ(UTα − UTβ){∑γC(1)αγUγ −∑γC(1)αγUα}=2∑α∑β(c · nαβ)+(UTα − UTβ)kαβ∑γC(1)αγ (Uγ − Uα) (3.28)Likewise, the third order terms in the third order spatial discretization of the advectionschemes which are the the fourth and the fifth terms on the right hand side of Eq. (3.25)are transformed to yield:2∑α∑β(c · nαβ)+(UTα − UTβ)k(2)αβ∑γC(2)αγ (Uγ − Uα) (3.29)Thus, the third order quadratic energy for the advection problem, analogous to Eq.(1.64), is obtained by combining Eq. (3.26), (3.28) and (3.29) for its first order, secondorder and third order terms respectively:ddt(U,MU) = −∑α∑β(c · nαβ)+|δUαβ|2 +2∑α,β(c · nαβ)+δUTαβ∑γ[C(1)αγ kαβ + C(2)αγ k(2)αβ ]δUγα (3.30)99where δUαβ is the solution difference, (Uα − Uβ) between cell α and a nearby cell, β, inits reconstruction stencil. The local reconstruction map coefficients are now:rαβ = C(1)αγ kαβ + C(2)αγ k(2)αβ (3.31)wherek(2)αβ =ˆΓαβ12(~x− ~xα)T (~x− ~xα)dx (3.32)Eq. (3.32) is the second order equivalent of the Eq. (3.20) where the geometricterms are the distances from the cell reference point of control volume α to the controlvolume boundary with the control volume β. These terms appear when the Taylor seriesexpansion of the solution reconstruction, Eq. (1.10), are integrated within the controlvolumes. However, these terms are not generally measured analytically. Instead, Gaussquadrature rules are used to calculate these terms at the Gauss quadrature points at thecontrol volume boundaries. The first term in Eq. (3.30) is always negative. Therefore,if one can make the local reconstruction map (Eq. (3.31)) smaller in the second term inEq. (3.30), the rate of growth of energy would be smaller if not negative.3.1.4 Inviscid Burgers’ ProblemFor investigating nonlinear problems, I start with the inviscid Burgers’ equation:∂tu(x, t) + u(x, t) · ∇u(x, t) = 0 (3.33)where after the upwind spatial discretization, the semi-discrete set of ODE’s are:|να| dU(t)dt α= −∑ ˆΓαβ((wα · nαβ)+wα + (wβ · nαβ)−wβ) dΓ (3.34)100where w is again the reconstructed solution at the Gauss quadrature point. Then thequadratic energy for the inviscid Burgers scheme using Eq. (1.64) is defined as:ddt(U,MU) = −2∑α∑βUTα [(wα · nαβ)+wα + (wβ · nαβ)−wβ] (3.35)First orderFor the first order inviscid Burgers’ problem there are no reconstruction terms; thereforeusing the generalizations described in Eq. (3.6) and (3.5), the first order terms in Eq.(3.35) are reduced to:ddt(U,MU) = −2∑α∑βUTα [(Uα · nαβ)+Uα + (Uβ · nαβ)−Uβ]= −2∑α∑βUTα [(Uα · nαβ)+Uα − (Uβ · nαβ)+Uβ]= −∑α∑β[(Uα · nαβ)+(UTα Uα + UαUTα)−(Uβ · nαβ)+(UTα Uβ + UTβ Uα)]}(3.36)Now there are two scenarios:(Uα · nαβ)+ ≥ (Uβ · nαβ)+(Uα · nαβ)+ ≤ (Uβ · nαβ)+(3.37)In either case there is a maximum such as U+α ·nαβ, which leads to the maximum growthof energy. In this way the growth of the quadratic energy for the first order Burgers’scheme, Eq. (3.36), is transformed to:ddt(U,MU) ≤ −∑α∑β[(U+α · nαβ)+(U2α + U2β)−(U+α · nαβ)+(UTα Uβ + UTβ Uα)]101= −∑α∑β(U+α · nαβ)+ (Uα − Uβ)2 ≤ 0 (3.38)Hence, the first order problem is always energy stable.Second orderIn the second order problem, with linear reconstruction, the quadratic energy is composedof two parts: the first part comes from the first order terms (the first two terms on theright hand side of Eq. (3.39)) which are always negative. Therefore they do not have anycontribution in the growth of energy. On the contrary, the second part coming from thereconstructed terms (the last three lines on the right hand side of Eq. (3.39)) can havesome unfavorable unstable effects. Using the two generalizations introduced in Eq. (3.5)and (3.6) and the first order energy equation, Eq. (3.38), the overall quadratic energyfor the second order Burgers’ scheme is similarly calculated by:ddt(U,MU) =−∑α∑β(U+α · nαβ)+(δUαβ)2−∑α∑β(kαβ∑γCαγ (δUγα) · nαβ)+(δUαβ)2−∑α∑β(U+α · nαβ)+(δUTαβ)kαβ∑γCαγ︸ ︷︷ ︸rαβ(δUγα)−∑α∑β(kαβ∑γCαγ (δUγα) · nαβ)+×(δUTαβ)kαβ∑γCαγ︸ ︷︷ ︸rαβ(δUγα)(3.39)102Any reconstruction scheme which is the minimizer of the last two terms in Eq. (3.39)is the optimal one. Unfortunately there is no optimal solution to the minimization ofthese two terms. However, it has been empirically observed that any increase in the sizeof the reconstruction stencil decreases the maximum norm of the reconstruction map ineach cell (see [7] for more details). This in turn leads to a lower growth of energy whichmeans more stable solutions.3.1.5 Euler and Navier-Stokes ProblemsTo solve the Euler and Navier-Stokes problems I use Roe’s scheme (see [25]) to ap-proximate the fluxes at the cell interfaces. In this flux differencing scheme, the flux iscomposed of an advective operator and a diffusive operator which amounts to splittingthe flux differences at the cell interfaces into left-moving and right-moving components.With the same notation as previous sections, the flux at the control volume boundary ofCVα and nearby CVβ is:Fαβ =12 (Fα + Fβ)−12∣∣∣Aαβ (U˜)∣∣∣ (Uβ − Uα) (3.40)where U˜ is defined by Roe’s scheme. After the spatial discretization, the semi-discreteset of ODE’s are:|να| dU(t)dt α= −∑βˆΓαβ12[Fα + Fβ −∣∣∣A˜αβ∣∣∣ (Uβ − Uα)] dΓ (3.41)By using the third condition of Roe’s scheme, and decomposing the A˜ matrix to the leftand right moving components I have:A˜αβ (Uβ − Uα) = Fβ − FαA˜αβ = A˜+αβ + A˜−αβ,∣∣∣A˜αβ∣∣∣ = A˜+αβ − A˜−αβ , (3.42)Now by substituting Eq. (3.42) into Eq. (3.41), the spatial discretization is reduced to:103|να| dU(t)dt α= −∑βˆΓαβFαdΓ︸ ︷︷ ︸=0− 12∑βˆΓαβ(A˜αβ −∣∣∣A˜αβ∣∣∣) (Uβ − Uα) dΓ= −∑βˆΓαβA˜−αβ (Uβ − Uα) dΓ (3.43)Then by applying the quadratic energy I have:ddt(U,MU) = −2∑α∑βUTα A˜−αβ (Uβ − Uα) (3.44)where α iterates over all CVs, and β iterates over the CVs sharing a face with α. However,one can extend β to iterate over all the CVs as well by defining A˜αβ as zero when the twocontrol volumes α, β do not have any face in common. In this way, one can interchangeα, and β. Hence, the quadratic energy for the first order problem without any datareconstruction is:ddt(U,MU) ==∑α∑β(UTα A˜−αβUα + UTβ A˜−αβUβ − UTα A˜−αβUβ − UTβ A˜−αβUα)=∑α∑β(Uβ − Uα)T A˜−αβ (Uβ − Uα) ≤ 0 (3.45)Since A˜−αβ by definition is negative semi-definite, Eq. (3.45) is always less than zerowhich proves the energy stability for the first order Euler problem. However, a piece-wiselinear reconstruction of the solution at Gauss quadrature points (Eq. (3.9)) changes thequadratic energy to be a function of the local reconstruction map as well:104ddt(U,MU) =∑α∑β(Uβ − Uα)T A˜−αβ (Uβ − Uα)+∑α∑β(Uα − Uβ)T A˜−αβkαβ∑γCαγ︸ ︷︷ ︸rαβ(Uγ − Uα) (3.46)In Eq. (3.46), γ iterates over all the CVs in the reconstruction stencil of CVα, and rαβis the piece-wise linear local reconstruction map as defined earlier. From Eq. (3.46), thegrowth in the quadratic energy can only come from the second term which is a functionof the local reconstruction map. For example, if only the first neighbor control volumesare used in the reconstruction, γ, and β can be used interchangeably which in turn meansthe solution part of the second part in Eq. (3.46)(∑α∑β(Uα − Uβ)T A˜−αβ (Uγ=β − Uα))isalways positive. Hence, the smaller the value of the local reconstruction map rαβ is, theless positive the growth of the quadratic energy is going to be.3.1.6 Numerical ResultsThe results of Haider et al. showed that for the 2nd order advection problem the morecells there are in the stencil, the smaller the norm of the local reconstruction map is.The numerical experiments confirm this (see Table 3.1). Moreover, Table 3.1 backs thisobservation by showing that the rightmost eigenvalue, λmax, gets smaller (more negative)as the stencil size increases. In this chapter, the right most eigenvalue of the Jacobianmatrix of the semi-discretized scheme is used as a proxy for the numerical stability of thesteady state solution. There seems to be an apparent threshold of stability too, meaningthat when the maximum norm of the local reconstruction map, Rmax, is larger than one,the problem is likely to be unstable, whereas when the norm is smaller than one it ismore likely to be stable. Moreover, Figure (3.1) shows the same mesh that is used forthe stability analysis shown in Table (3.1) and how the norm of the local reconstruction105Stencil Size Rmax Rave λmax #+3 1.28 0.46 1.47 14 1.11 0.41 1.05 15 0.63 0.32 -2.70 06 0.45 0.27 -4.06 0Table 3.1: Stencil size and stability analysis for a second order 2D advection problem onan unstructured mesh for a channel. Rmax and Rave are the maximum and the averageof the norm of the local reconstruction map respectively. λmax shows the rightmosteigenvalue of the semi-discretized Jacobian while #+ stands for the number of positiveeigenvalues. For a triangular mesh in two dimensions, the first layer of the stencil consistsof three control volume neighbors with the exception of boundary control volumes.Stencil Size Rmax Rave λmax #+4 4.40 0.87 3.51 25 1.68 0.63 1.08 16 1.17 0.53 0.12 17 0.95 0.45 -0.28 0Table 3.2: Stencil size and stability analysis for a second order 3D advection problem onan unstructured mesh for a channel. Rmax and Rave are the maximum and the averageof the norm of the local reconstruction map respectively. λmax shows the rightmosteigenvalue of the semi-discretized Jacobian while #+ stands for the number of positiveeigenvalues. For a tetrahedral mesh in three dimensions, the first layer of the stencilconsists of four control volume neighbors with the exception of boundary control volumes.map of every control volume decreases by increasing the stencil size. Furthermore, theseresults carry over to second order 3D advection problems as seen in Table (3.2). Fora tetrahedral mesh in three dimensions, the first layer of the stencil consists of fourcontrol volume neighbors with the exception of the boundary control volumes. Table(3.2) shows that by increasing the stencil size from only the first neighbors to the secondneighbor control volumes, both the maximum norm of the local reconstruction map andthe rightmost eigenvalue of the semi-discretized problem decrease. This in turn signals amore stable problem. Likewise, Figure (3.2) demonstrates how the local reconstructionnorm of all the boundary control volumes decreases by increasing the stencil size for asecond order 3D advection problem.Unlike second order advection, where there is a clear threshold of stability, for high106(a) Stencil size 3(b) Stencil size 4(c) Stencil size 5(d) Stencil size 6Figure 3.1: Norms of the local reconstruction map for second order 2D advection problemwith different stencil sizes on a two dimensional channel with 1500 control volumes.107(a) Stencil size 4(b) Stencil size 5(c) Stencil size 6(d) Stencil size 7Figure 3.2: Norms of the local reconstruction map for second order 3D advection problemwith different stencil sizes on a three dimensional channel with 650 control volumes.108(a) Stencil size 11(b) Stencil size 12(c) Stencil size 15(d) Stencil size 24Figure 3.3: Norms of the local reconstruction map for third order 3D advection problemwith different stencil sizes on a three dimensional channel with 2400 control volumes.Since the local reconstruction map norms vary over a large range of values, these plotsuse a smaller range than the maximum range of values to better show the differences indata.109St. Sizecoarse mesh (653 CVs) fine mesh (2399 CVs)Rmax Rave λmax #+ Rmax Rave λmax #+11 33.07 2.06 8.25 4 44.17 3.42 5.93 412 18.53 3.24 1.78 2 20.50 2.71 -6.77 015 6.01 1.89 -4.70 0 6.492 1.58 -9.69 024 2.17 1.04 -6.34 0 2.73 0.88 -9.82 0Table 3.3: Stability analysis using the local reconstruction map and the rightmost eigen-values of the semi-discrete Jacobian for a 3D, 3rd order advection problem on an unstruc-tured mesh.ReIm-10 -5 0 5-15-10-505101524 Neighbors11 NeighborsFigure 3.4: First 25 eigenvalues for a 3D mesh 3rd order advection problem for 2 differentstencil sizes110order schemes there is no apparent bound for the value of the local reconstruction mapthat signals stability. However, as shown in Table 3.3, as the number of neighbors in thereconstruction stencil increases, the value of the reconstruction map norm decreases; thistrend is accompanied by a decrease in the number of positive eigenvalues of the semi-discrete Jacobian as well as the positivity of the rightmost eigenvalue. Likewise, Figure(3.3) demonstrates the decreasing trend of the local reconstruction map norms reportedin Table (3.3) for boundary control volumes. These results support the argument thatincreasing the stencil size always tends to improve stability (see Figure 3.4). However,in this first stage of the stability analysis, there are no mechanisms for deciding howto increase the stencil size for each control volume in the mesh. This is in fact themain drawback of the regular methods for increasing the stencil size of the solutionreconstruction where the control volumes are only added to stencil based on the order thatthey have been found either in the neighborhood of a control volume or randomly. In thesemethods in a pre-processing step, multiple layers of vertex neighbors or face neighborsof control volumes in the mesh are found where vertex neighbors are the neighbors of acontrol volume connected through a vertex and face neighbors are the ones that share aface with the control volume. In this thesis, the face neighbors and the face neighborsof these face neighbors and so on are found for each control volume in a pre-processingstep. However, in the next parts of this chapter, I present multiple strategies using thenorm of the local reconstruction map, geometrical properties of the control volumes,and using eigenvalue analysis to automatically select the control volumes in a solutionreconstruction stencil to obtain better numerical stability properties.From Table 3.4, it can be seen that increasing the stencil size has local effects onreconstruction. This means that one can increase the stencil size only in those cells (seeFigure (3.5)) that have large reconstruction map norm values and which are more likely tomake the solution unstable. By increasing the stencil size for only five percent of controlvolumes for an inviscid Burgers’ problem, both Rmax and λmax decrease. This will save us111SmallStencil Size% Cells withLarge StencilLargeStencil SizeRmax Rave λmax #+3 0 9 1.12 0.44 5.026 23 1 9 0.97 0.43 5.026 23 5 9 0.67 0.41 -0.067 03 10 9 0.53 0.39 -8.63 03 100 9 0.42 0.26 -9.43 0Table 3.4: 2D, 2nd order inviscid Burgers’ problem with 1382 cells.both memory and computation cost compared with using large stencils everywhere whileimproving the robustness of the method compared with using small stencils everywhere.3.2 Ordering the StencilAny correlation between the value of the local reconstruction map and the mesh proper-ties will enable us to form the reconstruction stencil wisely and know in advance whichstencil configuration will lead to more stable solutions. In this effort, it was observed em-pirically that the mesh property correlates best with the value of the local reconstructionmap is the so-called relative moments x̂nymαγ (Eq. (1.12)) in the reconstruction matrixof control volume α (Eq. (1.14)) and the distance between CVα and the neighbor CVγto calculate the mean value of the reconstruction function for a control volume γ in thestencil. This experimental ordering stems from the calculation of the local reconstruc-tion map. In the local reconstruction map of control volume α, Eq. (3.19), there aretwo geometric terms: kαβ, which is the distance between the cell reference of α and thelocation of the Gauss quadrature points; this value is constant when changing the stencilsize for the same order of accuracy, as the number and the location of Gauss quadraturepoints for a control volume is only dependent on the order of accuracy. However, thesecond term, Cαγ, which is the pseudo inverse of the reconstruction matrix of controlvolume α and its neighbors γ, Eq. (1.14), is strongly dependent on the control volumes112(a) 1% of control volumes with an extended stencil size(b) 5% of control volumes with an extended stencil size(c) 10% of control volumes with an extended stencil sizeFigure 3.5: 2D, 2nd order inviscid Burgers’ problem on an unstructured mesh of a twodimensional channel with with 1382 cells. The highlighted control volumes have a largerstencil size of 9.113in the solution reconstruction stencil. In fact, by observing Eq. (1.14) and (1.15), onecan deduce that the contribution of each control volume γ in the stencil of α, comes fromtwo terms in the rows of matrix Eq. (1.14): the distance between the control volume αand γ; and the relative moments of control volume γ from control volume α as calculatedin Eq. (1.12). For these reasons, the effects of these two terms on the local reconstruc-tion map were studied. It was observed that if one reorders the control volumes in theneighborhood of CVα in a descending order with a specific criterion and choose the firstones in this list for the reconstruction stencil, one can decrease the local reconstructionmap norm for a given stencil size as well as improve the stability and convergence rate.By this criterion CVγ has priority over CVγ+1 if:LL+min{∣∣∣x̂lyk∣∣∣}Ll+kα,γ+1≥LL+min{∣∣∣x̂lyk∣∣∣}Ll+kα,γ, l ⊆ [1n] , k ⊆ [1m] (3.47)where L is the distance between the reference location of CV α, and γ, as is shown inEq. (1.15), and L is the minimum local length scale at the cell α.3.2.1 Numerical ResultsTable 3.5 shows how implementing this criterion in ordering the control volumes in thereconstruction stencil can influence the stability as well as the value of the local recon-struction map for two subsonic and transonic Euler problems around a NACA 0012 airfoilwith Mach numbers of 0.5 and 0.8 respectively, and an angle of attack of 2◦. From Table3.5, it can also be observed that even for the Euler problem the value of the reconstruc-tion map norm as well as the stability of the problem improves as the number of thecontrol volumes in the reconstruction stencil increases. Moreover, Figures 3.6 and 3.7further justify the use of the suggested ordering of control volumes in the reconstructionstencil.114Problem StencilSizeUnordered Ordered by Eq 3.47Rmax Rave λmax Rmax Rave λmaxSubsonic3 0.93 0.46 8.386 0.59 0.14 -0.0334 0.63 0.43 5.527 0.60 0.15 -0.0335 0.50 0.34 -0.032 0.36 0.13 -0.033Transonic3 0.93 0.46 11.219 0.59 0.14 -0.0324 0.63 0.43 7.594 0.60 0.15 -0.0315 0.50 0.34 -0.031 0.36 0.13 -0.032Table 3.5: Stability analysis using the local reconstruction map norm and the right-most eigenvalues of the semi-discrete Jacobian for a 2D, 2nd order Euler problem on anunstructured mesh around NACA 0012 airfoil with 1500 control volumes.IterationsL 2 Norm of Residuals5 10 15 20 2510-710-510-310-1101103105Stencil 3, Re-orderedStencil 3, Un-orderedStencil 4, Re-orderedStencil 4, Un-orderedStencil 5, Re-orderedStencil 5, Un-orderedFigure 3.6: Convergence history for transonic Euler problem with Mach 0.8 and angle ofattack of 2◦ around a NACA 0012 airfoil.115IterationsL 2 Norm of Residuals5 10 15 20 2510-710-510-310-1101103105Stencil 3, Re-orderedStencil 3, Un-orderedStencil 4, Re-orderedStencil 4, Un-orderedStencil 5, Re-orderedStencil 5, Un-orderedFigure 3.7: Convergence history for subsonic Euler problem with Mach 0.5 and angle ofattack of 2◦ around a NACA 0012 airfoil.116(a) Unordered(b) Ordered by Eq. (3.47)Figure 3.8: The local reconstruction map of a 2D, 2nd order Euler problem on an un-structured mesh around NACA 0012 airfoil with 1500 control volumes with and withoutordering of the solution reconstruction stencils by Eq. (3.47) with 5 control volumes inthe stencil.117Problem StencilSizeUnordered Ordered by Eq 3.47Rmax Rave λmax Rmax Rave λmaxSubsonic3 0.65 0.41 1.620 0.63 0.40 0.0274 0.63 0.40 0.027 0.60 0.34 -0.0285 0.54 0.32 -0.032 0.53 0.29 -0.081Transonic3 0.65 0.41 2.072 0.63 0.40 0.1554 0.63 0.40 0.112 0.60 0.34 -0.0295 0.54 0.32 -0.008 0.53 0.29 -0.042Table 3.6: Stability analysis using local Recon-map and the rightmost eigenvalues of thesemi-discrete Jacobian for a 2D, 2nd order laminar problem (angle of attack of 5◦, andReynolds of 5000) on an unstructured mesh around NACA 0012 airfoil with 5000 controlvolumes for both subsonic (Mach of 0.3) and transonic flow conditions (Mach of 0.72).Similar trends also carry over to the laminar flow where two examples are used toshowcase this. Table (3.6) shows that by ordering the solution reconstruction stencil usingthe criterion in Eq. (3.47), I was able to stabilize two initially unstable subsonic laminarflows. Table (3.6) shows that for the same stencil sizes, the right most eigenvalue of eachof the problems becomes less positive by ordering the stencils while the maximum valuesof the norm of the local reconstruction map drops as well. The increase of the stencilsize in both Table (3.5) and (3.6) from three to five is based on the general practice ofincreasing the stencil size described in Table (3.1), (3.2) and (3.3), where control volumesare added to the stencil in the order they are found in the mesh for the unordered cases,while for the ordered cases, the stencils are extended based on Eq. (3.47). Additionally,Figures (3.8) and (3.9) show how the local reconstruction map for each control volume inthe two meshes used in Tables (3.5) and (3.6), decreases globally by ordering the stencils.3.3 Reconstruction Stencil Optimization byEigenanalysisIn the previous sections, the effects of the reconstruction stencil on numerical stabilityfor different physical problems by means of reconstruction maps was showed. It was118(a) Unordered(b) Ordered by Eq. (3.47)Figure 3.9: The local reconstruction map of a 2D, 2nd order laminar Navier-Stokesproblem on an unstructured mesh around NACA 0012 airfoil with 5000 control volumeswith and without ordering of the solution reconstruction stencils by Eq. (3.47) with 5control volumes in the stencil.119shown that the shape and size of the stencil can move a discretization towards stabilityor otherwise by contributing to the growth of energy. However, the reconstruction mapanalysis is independent of the solution. In other words, the effects of the reconstruction onthe solution are not observed; hence any secondary stability effects of the reconstructionthrough the solution variables are not studied. In Chapter (2), it was shown that gradientsof the rightmost eigenvalues of the global Jacobian matrix, Eq. (1.18), can be used tooptimize the mesh vertices in a finite volume discretization on unstructured meshes togain better stability and convergence results. In this section, the same notion with asimilar optimization approach is used to obtain a stabilized reconstruction stencil sizeand direction. Instead of simply increasing the stencil size for all of the control volumesbased on the reconstruction map, the goal is to systematically determine the stencil sizeand direction of the reconstruction stencil for each control volume based on the gradientsof the rightmost eigenvalues of semi-discretized Jacobian to stabilize the final steady stateproblem.As shown previously, a reconstruction stencil can include multiple layers of the controlvolume neighbors to construct the reconstruction matrix and to approximate the solutiongradients within the control volume (see 1.1.2). The first neighbors of a control volume areits face neighboring CVs; the second neighbors are the face neighbors of its first neighbors,and so on. However, this will result in a more or less central reconstruction stencil. Tosystematically determine an upwind, downwind, or central reconstruction stencil, I firstcategorize all the neighbors of a control volume as upwind or downwind using the cosineof the angle between the local velocity and the line connecting the control volume to itsneighbor (see Figure 3.10). For upwind reconstruction, control volumes with negativecosines are selected to be part of the reconstruction stencil. In each category, neighborswhich are closer to the control volume reconstruction have higher priority to be part ofthe reconstruction stencil. For central, the cosine is discarded and the reconstructionstencil is constructed based only on the layers of the neighbors. Figure (3.11) better120» »¼ à L ® »¼ ® »¼ ¼ Figure 3.10: Determining the upwind stencil for a control volume, α, where β is a controlvolume in the stencil of the α. −→U is the local velocity at control volume α.(a) Upwind stencil (b) Regular (central) stencilFigure 3.11: Upwind stencil versus regular stencil for three control volume samples in anunstructured mesh with 1500 control volumes around a NACA 0012 airfoil for an Eulerscheme with M = 0.72, α = 3◦.illustrates the difference between the regular and the upwind stencil selection for threecontrol volumes in the mesh. In the rest of this section, I describe how optimization candetermine whether to use an upwind or a central reconstruction stencil for each controlvolume in the problem.The idea of this work hinges on the gradient of the eigenvalues of the semi-discretizedJacobian matrix, Eq. (1.9), with respect to the parameters of interest which are specifiedweights for each control volume in this work (see Section 2.1.2 for more details on howto calculate the gradients of eigenvalues with respect to the parameters of interest).Gradients of eigenvalues for large generalized eigenvalue problems can be approximated121by measuring the perturbations of the eigenmodes with respect to perturbations in thesystem [93, 1]. For each control volume α, two weights of wα, and 1−wα are attributedrespectively to the smaller stencil (the first layers of the neighbors) and the larger stencil(the first, second, and farther layers of the neighbors). Similar weights are applied toevery control volume in the problem to optimize between upwind and centered stencils.The gradients of the eigenvalues with respect to these weights are then calculated byusing Eq. (3.48).∂λ∂wα= L ∂J∂wαR (3.48)where L and R are the left and right (Eq. (1.67) and (1.66)) eigenvectors of the Jacobianmatrix respectively, and λ is the rightmost eigenvalue. To calculate the derivative of theJacobian matrix with respect to the weights, I look into the chain rule differentiation inEq. (1.18) where only the pseudo-inverse of the reconstruction matrix C is dependent onthe reconstruction stencil. Therefore, I can define C as a linear combination of differentstencil sizes for each control volume α, Cα =∑iwα,i·Cα,i with weights defined as wi.When optimizing between two choices of stencil, for each control volume Cα is measuredby:Cα = wαCi1 + (1− wα)Ci2 (3.49)where Ci1 and Ci1 are the pseudo-inverse of the reconstruction matrices for each stencilchoice.Using chain rule differentiation and the global Jacobian differentiation, I calculatethe Jacobian Ji with weight wα and its derivative for each control volume:122Jα,i = IFRCα,iV (3.50)Jα =∑iwα,iJα,i (3.51)∂J∂wα= IFR∂Cα∂wαV (3.52)where all the terms except C in the Jacobian chain are independent of the weight wi. InEq. (3.52), the gradient of the pseudo-inverse of the reconstruction matrix, when thereare two stencil choices, is calculated by taking the derivative of Eq. (3.49) with respectto the weights for each control volume:∂Cα∂wα= Cii − Ci2 (3.53)Replacing Eq. (3.53) in Eq. (3.52), transforms Eq. (3.52) to:∂J∂wα= Jwα − J1−wα (3.54)Next, to move the rightmost eigenvalue to the left half of the complex plane, one cancalculate the gradients of the rightmost eigenvalues to update the weights wi in backpropagation through gradient descent:λn+1 = λn +−−→4w ∂λ∂w≤ 0 (3.55)−−→4w = − |K| ∂λ∂w(3.56)K = λn(∂λ∂w)2 (3.57)wopt = w +−−→4w (3.58)where K specifies the update size of the weights while the ∂λ∂wspecifies the direction in123which the gradient descent updates the weights. Additionally, this can be achieved moresystematically by defining an optimization problem such that the weights are minimizedwhile the projected rightmost eigenvalue is negative. In this problem smaller weightsare preferable as they are associated with the smaller stencil size. Smaller stencils arepreferred as they produce a more sparse Jacobian matrix; hence, a faster linear solve ofthe Jacobian matrix. To this end, I minimize the L2 norm of the weights update vector−→S = −−→4w:min{N∑i(4wi)2}(3.59)which is constrained by the following linear inequalities to ensure that the new rightmosteigenvalue is in the left half of the complex plane:λn+1j +−→S · ∂λj∂w≤ 0 1 ≤ j ≤M (3.60)where N is the number of control volumes and M is the number of unstable eigenvalues.This quadratic optimization problem is solved using the Interior Point Method. Aftercalculating the weights wi in either of the above methods and mapping them between 0to 1, the smaller or larger stencil size for each control volume is chosen based on whichend of the spectrum the weights are leaning more towards. For example, weights ofmore than 0.5 signal a shift towards the larger stencil. However, this threshold can becustomized by the user to favor one stencil size over the other. Similarly, the directionof the reconstruction stencil for any control volume in the mesh can be optimized toyield the best stability and convergence properties by determining a central or upwindreconstruction stencil using the same optimization scheme described here.Since the optimization is done for every control volume, the process is of the order ofO (N2) where N is the total number of control volumes in the problem. This problemarises when for every control volume, Eq. (3.48) recalculates the whole Jacobian and124» ¼Ú ¼Û ¼Ü Figure 3.12: Local Jacobian and reconstruction calculation. The control volume wherethe reconstruction stencil optimization is happening is shown with α, while β1, β2, andβ3 are the face neighboring control volumes.reconstruction matrix. To expedite the process, I need to recalculate both the Jacobianand the reconstructions locally. This means that for every control volume, a limitednumber of rows of the Jacobian matrix needs to be updated. Once the shape or sizeof the reconstruction stencil for a control volume CVα changes, only the pseudo-inverseof the reconstruction matrix for this control volume changes which in turn changes thesolution and the gradients of the solution at the Gauss quadrature points on the edgesof the control volume. Hence, only the flux integral and the rows of the Jacobian matrixfor the control volume itself and its face neighbor control volumes are changed; thusthey need to be updated. For example in Figure 3.12, if the size or the shape of thereconstruction changes for control volume α, the reconstruction changes for CVα onlyand the flux integral or Jacobian rows change only for control volumes α, β1, β2, and β3.3.3.1 Numerical ResultsIn this section, the numerical results for optimizations of both the reconstruction stencilsize and its direction are presented for the Euler and laminar Navier-Stokes problemsunder different flow conditions. To evolve to the steady state solution from the initialconditions on the spatially discretized Eq. (1.8) I use Crank-Nicolson time integration.Crank-Nicolson, which is an A-stable scheme, ensures that positive eigenvalues of the125semi-discrete Jacobian lead to the growth of the solution over time; hence instability. Inother words, any unstable mode in the spatially discrete equations will remain unstableand likewise all stable modes will remain stable. This will affirm that any instability ob-served when marching to the steady state solution originates from the spatial discretiza-tion rather than time discretization. In this section, in all the mesh figures depictedwith red and blue control volumes, red represents the control volumes with an increasedstencil size or an upwind direction while blue represents the control volumes which havea regular stencil size or a symmetric direction.3.3.1.1 Reconstruction Stencil SizeAs first test cases, reconstruction stencil size has been optimized for all the controlvolumes in the mesh for the second order transonic and supersonic Euler problems.Figures 3.13a and 3.14a show reconstruction stencil sizes selected for each control volumeafter optimization for a transonic Euler problem with free stream Mach number of 0.72and an angle of attack of 3°, and a supersonic Euler problem with free stream Machnumber of 1.21 and an angle of attack of 3° respectively. Likewise, Figures 3.13b and3.14b present the gradients of the rightmost eigenvalue with respect to the reconstructionstencil weights for each control volume. There is a strong tendency in both of the picturesto use larger stencil sizes in upwind of the shocks. As shown in Figure 3.15 both problemsare stabilized after the optimization. In both of these problems, the second order semi-discrete Jacobian is linearized by the first order steady state solution to the problem forpurpose of optimization.Similarly, the reconstruction stencil size is optimized for a O (h3) subsonic laminarNavier-Stokes problem with Mach number of 0.3, angle of attack of 2°, and Reynoldsnumber of 5000. Figures 3.16a and 3.16b show a strong tendency towards larger stencilsizes in the boundary layer and wake. In this problem, a second order steady statesolution was used to linearize the third order Jacobian for purpose of optimization.126(a) Reconstruction stencil size map. Red control volumeshave two layers of control volumes in their stencil whilethe blue ones have only the first layer of control volumeneighbors in their stencil.1e-111e-101e-91e-81e-71e-61e-50.00010.0010.011.000e-122.718e-01Gradient(b) Gradients of Eigenvalues with respect to stencil sizeFigure 3.13: Stencil size optimization; Second order transonic Euler on an unstructuredmesh with 1500 control volumes and with free stream Mach number of 0.72 and an angleof attack of 3°.127(a) Reconstruction stencil size map. Red control volumeshave two layers of control volumes in their stencil whilethe blue ones have only the first layer of control volumeneighbors in their stencil.1e-91e-81e-71e-61e-50.00010.0010.010.11.000e-109.930e-01Gradient(b) Gradients of Eigenvalues with respect to stencil sizeFigure 3.14: Stencil size optimization; Second order supersonic Euler on an unstructuredmesh with 1500 control volumes and with free stream Mach number of 1.21 and an angleof attack of 3°.128IterationsL 2 Norm of the residual20 4010-810-610-410-2100102Transonic, optimized stencilTransonic, fist layer stencilSupersonic, optimized stencilSupersonic, first layer stencilFigure 3.15: Residual history before and after stencil size optimization; Second ordertransonic (free stream Mach number of 0.72 and an angle of attack of 3°) and supersonic(free stream Mach number of 1.21, an angle of attack of 3°) Euler on an unstructuredmesh with 1500 control volumes.129(a) Reconstruction stencil size map. Red control volumeshave three layers of control volumes in their stencil whilethe blue ones have only the first and the second layer ofcontrol volume neighbors in their stencil.1e-111e-101e-91e-81e-71e-61e-50.00010.0010.010.11.000e-121.000e+00Gradient(b) Gradients of Eigenvalues with respect to stencil sizeFigure 3.16: Stencil size optimization; Third order subsonic laminar flow on an unstruc-tured mesh with 2400 control volumes and with free stream Mach number of 0.3, angleof attack of 2°, and Reynolds number of 5000.130(a) Pressure field(b) Axial velocity non-dimensionalized by the free stream sound speedFigure 3.17: A second order supersonic Euler wind tunnel with ramp for M = 2 and3500 control volumesAnother common problem in two dimensional inviscid cases is the multiple Mach re-flection problem (e.g., see [94]). This problem has different variants. In one scenario,the multiple Mach reflection happens when the supersonic incoming flow in a tunnelapproaches a sudden barrier such as an inclined ramp as is shown in Figure 3.17. Figure3.17 shows the pressure field as well as the axial velocity in the wind tunnel. The secondorder finite volume solver for this Euler problem with Mach number 2 is originally unsta-ble on an unstructured mesh with 3500 control volumes. The three rightmost eigenvaluesof this problem are λ0 = 0.125, λ1,2 = −0.016 ± i1.482 while after the optimization therightmost eigenvalues are λ0 = −0.005, λ1,2 = −0.015± i1.478 respectively. A first ordersteady state solution to this problem was used to approximate the semi-discrete Jaco-bian of the second order discretization for purpose of optimization. Figure 3.18 showsthe control volumes, in the mesh highlighted in red, where the optimization has foundit necessary to have a larger stencil size (first neighbors and second neighbors) for thescheme to be stable, while the rest of the control volumes in blue have just one layer ofcontrol volume neighbors in their stencils.131Figure 3.18: Stencil size optimization; Second order supersonic Euler problem on anunstructured mesh with 3500 control volumes and with inflow Mach number of 2.3.3.1.2 Reconstruction Stencil DirectionHere, I stabilize finite volume discretizations on unstructured meshes by optimizing re-construction stencil directions for each control volume in the mesh. The optimizationprocess chooses central or upwind reconstruction stencil for each control volume. Thefirst two examples, shown in Figures 3.19 and 3.20, cover the optimization for second or-der transonic and supersonic flows respectively with free stream Mach number of 0.72 onan unstructured mesh with 1500 CVs and Mach number of 1.21 on an unstructured meshwith 5000 control volumes around an 0012 airfoil. Both cases have an angle of attack of3°. Figures 3.19a and 3.20a show the locations in the mesh where upwind discretizationhas been preferred over central discretization for transonic and supersonic cases respec-tively. The sensitivity of the numerical stability to these discretizations are shown in thecontours of the gradients of the rightmost eigenvalues with respect to the reconstructionstencil choices in Figures 3.19b and 3.20b. Figure 3.21 shows the convergence historyfor both transonic and supersonic problems with and without reconstruction stencil op-timization.The same stencil direction optimization is also applicable to other physical problemsand flow conditions. Thus, I applied the stencil direction optimization between centraland upwind discretization to a second order subsonic laminar flow with 5000 controlvolumes and with a free stream Mach number of 0.42, an angle of attack of 1°, anda Reynolds number of 5000. In Figure 3.22a, I show the stencil direction choice afteroptimization. Figure 3.22b shows control volumes in the mesh where the sensitivity of the132(a) Reconstruction stencil direction map. Red control vol-umes have an upwind stencil while the blue ones have a sym-metric stencil.1e-61e-50.00010.0010.010.113.693e-073.321e+00Gradient(b) Gradients of Eigenvalues with respect to stencil direc-tionFigure 3.19: Stencil direction optimization; Second order transonic Euler on an unstruc-tured mesh with 1500 control volumes and a free stream Mach number of 0.72 and anangle of attack of 3°.133(a) Reconstruction stencil direction map. Red control vol-umes have an upwind stencil while the blue ones have a sym-metric stencil.1e-61e-50.00010.0010.010.11101001.0e-073.0e+03Gradient(b) Gradients of Eigenvalues with respect to stencil direc-tionFigure 3.20: Stencil direction optimization; Second order supersonic Euler on an unstruc-tured mesh with 5000 control volumes and a free stream Mach number of 1.21 and anangle of attack of 3°.134IterationsL 2 Norm of the residual20 40 60 80 10010-810-610-410-2100102Transonic, optimized stencilTransonic, original stencilSupersonic, optimized stencilSupersonic, original stencilFigure 3.21: Residual history before and after stencil direction optimization; Second ordertransonic (free stream Mach number of 0.72 and angle of attack of 3°) and supersonic(free stream Mach number of 1.21, angle of attack of 3°) Euler on an unstructured meshwith 1500 control volumes.135stability (i.e. the gradient of the rightmost eigenvalue) to the choice of the reconstructionstencil direction (upwind or central) is higher. In this problem, the half-converged secondorder solution was used to linearize the second order Jacobian. In Figure 3.23 wherethe two residual lines detach is the iteration at which I applied the optimization tothe half-converged problem. It also shows how the problem quickly converges after theoptimization step unlike the original discretization, which continues to oscillate with nosign of convergence. Similar to Table (5.1), I present the differences of the Jacobianmatrices and their eigensystem of the schemes used in this section in Table (5.2).136(a) Reconstruction stencil direction map. Red control vol-umes have an upwind stencil while the blue ones have asymmetric stencil.1e-71e-61e-50.00010.0010.010.11.000e-081.000e+00Gradient(b) Gradients of Eigenvalues with respect to stencil direc-tionFigure 3.22: Stencil direction optimization; second order subsonic laminar flow on anunstructured mesh with 5000 control volumes and with free stream Mach number of0.42, angle of attack of 1°, and Reynolds number of 5000.137IterationsL 2 Norm of Residual20 40 60 80 100 120 14010-710-510-310-1101Original StencilOptimized Stencil after 100 IterationsFigure 3.23: Residual history before and after stencil direction optimization; second ordersubsonic laminar flow on an unstructured mesh with 5000 control volumes and with freestream Mach number of 0.42, angle of attack of 1°, and Reynolds number of 5000.138Chapter 4Boundary Condition OptimizationIn the previous chapters, I showed how the mesh and the solution reconstruction canaffect the numerical stability and how stability can be improved by optimizing the meshvertices or the solution reconstruction stencil respectively. However, often numericalinstabilities arise from the boundary discretization or boundary conditions as the interiorof the computational domain is normally well-posed. To this end, we aimed to investigateboundary conditions in terms of stability in this chapter.Boundary conditions play an important role in defining and closing the mathematicalproblems of computational fluid dynamics. The proper implementation of the boundaryconditions for different physical problems and discretizations is an active field of research(e.g., see [95, 96, 97, 98, 99, 100]). Likewise, in computational aerodynamics problems,the boundary conditions specified to approximate the physical solution at the boundariesof the computational domain are of paramount importance and have been the subjectof many studies. Thomas and Salas [101] showed that the numerical accuracy of thesteady state solution can be greatly affected by careless implementation of boundaryconditions at the truncated farfield boundaries. Spurious non-physical waves originatingat or reflected from the farfield boundaries and the solid surfaces can contaminate theinterior solution, in particular the aerodynamic forces. In fact, they show that the errorsand residuals of conserved solutions of the Euler and Navier-Stokes equations undergo thesame dissipation, dispersion, and reflection as the physical waves in the domain. Mazaheri[102] showed that there are mainly three multiplier factors in the attenuation of thesewaves: numerical dissipation in the interior domain, reflection of waves from the exterior139of the domain and the farfield, and the reflection of waves from solid surfaces within thedomain. To minimize the reflection of the waves at solid surfaces, wall treatments suchas soft walls [102] were introduced. Likewise, there have been numerous efforts to findschemes that would control the dissipation within the domain such as local time-stepping(e.g., see [18, 29, 103]), preconditioning (e.g., see [104, 16, 105, 106]), and multigridmethods (e.g., see [107, 108, 109, 110, 111]). Moreover, there have been many differentefforts to derive boundary conditions at the farfield which asymptotically approximatethe farfield free stream conditions, and thus reduce the reflection of the error waves at thefarfield boundary (e.g., see a review on these boundary conditions by Tsynkov [112] andthe references therein). One of the first methods developed to approximate the analyticalsolution at the farfield corrects the constant free stream conditions with a vortex flow,which is dependent on the airfoil lift (see [101]). This method is successful in reducing thesensitivity of the forces on the airfoil with respect to the location of the farfield boundary,albeit with the added cost of introducing non-local features (i.e. lift coefficient) to theboundary conditions. In another method for reducing non-physical oscillations, non-reflecting boundary conditions are imposed at the farfield. These boundary conditions,which were first introduced by Engquist and Majda [113] and further improved by [114,115], direct the oscillations out of the physical domain. Recently, there have been morework on these types of boundary conditions. In another attempt, Allmaras et al. [116]derived a boundary condition for the farfield which is both local and independent of theforces on the airfoil, relying on their modified variables and radial velocity. They showedthat by applying this boundary condition, the farfield boundary can be brought muchcloser to the airfoil.One of the shortcomings in all research papers in addressing the use of boundaryconditions is the lack of an applied approach to systematically choose the most stableboundary conditions amongst the set of all available boundary conditions. Thus, thenovelty of this chapter is that instead of introducing a specific type of boundary condition140which may not be always stable in practice, a framework is devised which can find themost stable boundary condition for a given problem. This new framework uses thegradients of the rightmost eigenvalues of the semi-discretized Jacobian with respect tothe different types of the boundary conditions to optimize boundary conditions throughgradient descent.The layout for this chapter is as follows. Section 4.1 describes different types ofcommonly-used farfield boundary conditions as well as wall boundary conditions for 2Dexternal aerodynamics. Section 4.2 introduces a systematic approach to optimize differ-ent types of boundary conditions to improve the numerical stability by calculating thederivatives of the rightmost eigenvalues of the semi-discretized Jacobian with respect toboundary condition parameters. Finally, Section (4.3) showcases the numerical stabilityresults for the aforementioned optimization approach.4.1 Boundary Conditions4.1.1 Farfield Boundary ConditionsIn this sub-section, I review some of the common boundary conditions used for externalcompressible subsonic flows. If all variables were known at the boundary, evaluation ofEq. (1.8) would be seamless. However, this is not usually the case. To determine all thestates at boundaries for the discretized equations, numerical and physical conditions needto be set. As each characteristic direction carries a piece of information into or out of thediscretized domain, only these transported variables can be set freely at the boundary;these are the physical boundary conditions. The rest of the variables needed by the setof discretized equations, namely the numerical boundary conditions, can be determinedby the discretization scheme itself. For example, Figure (4.1) shows the characteristiclines for a subsonic outflow in a two-dimensional flow.All numerical and physical boundary conditions in this thesis are enforced weakly141Figure 4.1: Characteristic lines for a subsonic two-dimensional flow at an outlet for flowin the x-direction. The y direction is normal to the page.through proper flux functions. The following shows different types of boundary conditionsused for a two-dimensional subsonic flow.4.1.1.1 Characteristics BC with vortex correctionThe boundary conditions can be fixed by specifying the propagated variables directlyfrom the characteristic information or similarly by the Riemann invariants. In Eq. (4.1)and (4.2) the primitive variables that are directly transferred with the characteristicvariables are determined both at in-flow and out-flow. The rest of the required variablesare determined from the interior discretization.• In-flow:Right-moving acoustic wave: (p− p∞)− ρc∞ (Vn − V∞n) = 0 (4.1)Entropy: (p− p∞)− c2∞ (ρ− ρ∞) = 0Tangential velocity: (Vt − V∞t) = 0142• Out-flow:Left-moving acoustic wave: (p− p∞) + ρc∞ (Vn − V∞n) = 0 (4.2)where variables subscripted with∞ represent farfield variables and c is the sound velocity.In external flows (and similarly internal flows) disturbances at the farfield (or inlet andoutlet) can be transferred to the interior of the flow. To minimize these effects, thephysical boundaries are normally placed far from the domain of interest.However, by better approximating the solution at the farfield, the farfield boundarycan be taken closer to the flow region of interest (e.g. airfoil) without affecting the flowsolution. For the inviscid, compressible flow condition, the farfield flow perturbations canbe modeled by the linear potential equation. The farfield uniform free stream variablesin Eq. (4.1) and (4.2) can be corrected by asymptotic perturbations to the lowest orderfrom the thin airfoil theory and the linear potential equation (see [101] for more details):V∞x = M∞c∞ cos (α) +z sin (θ)V∞y = M∞c∞ sin (α) − z cos (θ)(4.3)where M∞ is the free stream Mach number, α is the angle of attack, and c∞ is the freestream sound velocity. The parameter z is:z = ClLM∞c∞4pi√1−M2∞r11−M2∞ sin2 (θ − α)(4.4)The perturbations in Eq. (4.4) are governed by the circulation of a vortex centered atthe airfoil quarter-chord. The circulation itself is dependent on the lift coefficient, Cl.The lift coefficient is measured around any arbitrary surface enclosing the airfoil. In Eq.(4.4), parameters L, r, θ are respectively the chord length, the radius measured from thecenter of the vortex to the boundary location, and the angle between the chord and theline connecting the center of the vortex to the boundary location.1434.1.1.2 Pressure BCAnother commonly used boundary condition at the farfield is specifying the pressure. Aspressure is constant across a viscous wake, this does not disturb the transport of viscousperturbations. Thus, imposing a constant pressure downstream of the flow is a commonpractice. In this boundary condition, we specify the free stream back pressure at outflow,and we use total pressure, total temperature, and angle of attack at the free stream todetermine the three physical boundary condition at the inflow using isentropic relations.• In-flow:Ptot, Ttot, α (4.5)• Out-flowp = p∞ (4.6)4.1.1.3 Dirichlet BCIn this type of boundary condition, all primitive (or conservative) variables are directlydetermined from free stream variables at farfield. The boundary fluxes are then calculatedusing these fixed values. This problem is over-specified at both inflow and outflow.However, it has been shown that the numerical errors tend to be localized.• In-flow:ρ, u, v, p (4.7)• Out-flowρ, u, v, p (4.8)4.1.1.4 Radial Velocity BCTo bring the far field boundary closer to the airfoil without affecting the lift and dragcoefficients, Allmaras et al. [116] suggested setting properties at the far field that are tofirst order independent of the airfoil lift:144• In-flow:Total Enthalpy:[γγ − 1pρ+ V22]=[γγ − 1pρ+ V22]∞(4.9)Entropy: pργ= p∞ργ∞Radial Velocity: Vr = ~V∞ eˆr• Out-flowVr = ~V∞ eˆr (4.10)In Eq. (4.9) and (4.10), eˆr is the radial unit vector calculated from the vortex center,situated at a quarter of the chord, to the boundary location.4.1.1.5 Non-reflecting outflow BCAs explained earlier, physical information crosses the boundaries with characteristic vari-ables. This in turn causes reflection of the perturbation waves at the boundaries. Forexample, specifying only the back pressure at outflow does not satisfy the perturbationsof the characteristic waves ∆W = (p− p∞)− ρc (Vn − V∞n); hence, pressure waves withmagnitudes of ρc (Vn − V∞n) will be reflected back into the domain. For steady-statecomputations with constant pressure boundary conditions, Rudy et al. [114] proposed amethod which damps out the pressure waves leaving the computational domain in thetransient and thus satisfying constant pressure at steady-state:• In-flow:Ptot, Ttot, α (4.11)• Out-flow∂p∂t− ρc∂Vn∂t+ ω (p− p∞) = 0 (4.12)where ω decreases from roughly 0.7− 0.8 in subsonic to 0.1− 0.3 in transonic flows.1454.1.2 Wall Boundary Condition4.1.2.1 Rigid WallAt a solid surface there is no flow into or out of the wall and the normal velocity at thewall is set to zero:un|w = 0 (4.13)Hence, at a solid wall only the pressure needs to be set. For a first order boundarycondition, this pressure at the wall (pw) is equal to the pressure in the next control volume(pβ):pw = pβ (4.14)where w denotes the conditions at the wall, while β defines the conditions at the controlvolume next to the wall.4.1.2.2 Characteristic Rigid WallIn this boundary condition, the wall is still rigid and there is no flow at the wall. Thus,the normal velocity at the wall is set to the normal velocity of the wall. However, thewall pressure is determined using the characteristic waves incident on the wall.Entropy, along the particle path: (pw − pβ)− c2β (ρw − ρβ) = 0Acoustic wave on the wall: (pw − pβ) + (ρc)β (uw − uβ)n = 0 (4.15)From the acoustic wave incident on the wall (Eq. (4.15)) and the zero normal velocity(un|w = 0) at the wall, the wall pressure and in turn the pressure components of the fluxare calculated:146pw = pβ + (ρc)β (uβ)n (4.16)4.1.2.3 Soft WallIf we assume that the wall is not fully rigid and has some elasticity, and will come to restat steady state, we can measure and model a normal velocity at the wall. This normalvelocity will decay to zero as elasticity of the wall becomes zero when approaching thesteady state condition. Mazaheri [102] showed that the normal velocity and the pres-sure at the wall can be modeled using a linear spring-damping analogy with a dampingcoefficient of ( 1µ), stiffness of ( 1τµ) and a force of ( 1ρcdpdt):dudt+ uτ= µρcdpdt(4.17)where µ and τ are free parameters that are specifically designed to minimize the transientand steady state ratios of the strength of the reflected waves over incident waves on thewall. Mazaheri [102] has shown that high values of µ increase wave reflection in thetransient but decrease the steady state wave reflection, while very small values of µ havethe opposite effect. It is experimentally found that µ = 0.5 is the best compromise andresults in better convergence. Likewise, the parameter τ is also bounded not to be toolarge or too small. If τ is much smaller than a (pseudo) time step, truncation errorsincrease. If τ is much greater than the round trip time for a wave, as is determined bythe slowest wave, the next returning wave from farfield will not see the proper boundarycondition. Again, it is empirically found that τ = 100∆t gives the best convergence. Thestiffness equation (Eq. (4.17)) along with the incident characteristic waves on the wallare solved to determine the flow conditions at the wall:147uwn − uβn4t +uβnτ= µ(ρc)βpw − pβ∆t (4.18)(pw − pβ)− c2β (ρw − ρβ) = 0 (4.19)(pw − pβ) + (ρc)β (uw − uβ)n = 0 (4.20)(uw − uβ)t = 0 (4.21)where uwn , uwt , pn, and ρw are then used in the flux functions to determine the boundaryflux at the wall.4.2 Methodology4.2.1 Boundary Condition OptimizationIn this work, I optimize boundary conditions to improve the numerical stability of thesemi-discrete system of equations arising from finite volume discretization on unstruc-tured meshes. Herein, I am looking into the stability effects of different boundary condi-tion implementations. Throughout this work, I apply boundary conditions weakly usingflux calculations. The premise of this work is to use the semi-discretized Jacobian of theflux integral, and its right most eigenvalues and eigenvectors to calculate the gradients ofeigenvalues of the Jacobian matrix with respect to the scheme design parameters b (theboundary flux weights).∂∂b(JRi = λiRi) (4.22)Li(∂J∂bRi + J∂Ri∂b= ∂λi∂bRi + λi∂Ri∂b) (4.23)(LiJ − Liλi) ∂Ri∂b= Li∂λi∂bRi − Li∂J∂bRi (4.24)148We take the derivative of the right eigensystem with respect to the variable of interestin Eq. (4.22). While keeping the solution U fixed, we multiply from the left by the lefteigenvector in Eq. (4.23). The left hand side of Eq. (4.24) contains the left eigenproblem,leaving us with Eq. (4.25). Since the eigenvectors are normalized so that Ri Li = 1, wetake the gradients of the eigenvalues:∂λi∂b= Li∂J∂bRi (4.25)where parameter b is the scheme design variable, which we are taking the derivativewith respect to. One key step in evaluating Eq. (4.25) is approximating the gradient ofthe Jacobian matrix with respect to the variables of interest, ∂J∂b. I chose to use finitedifferencing on the Jacobian matrix; however this gradient can also be calculated usingreverse automatic differentiation.Moreover, I am aiming to systematically select and optimize the boundary conditionsto stabilize the finite volume methods. To find the most stable boundary conditionconfiguration, I first linearly superimpose the different boundary fluxes for each boundarycontrol volume. Second, I optimize these superposition parameters based on the gradientof the right-most eigenvalues of the Jacobian matrix based on these parameters. For asimple case in which optimization is between only two types of boundary conditions, Icombine the two different boundary fluxes with weights bcvi and 1−bcvi for every boundarycontrol volume:Flux = bcviF1 + (1− bcvi)F2 (4.26)∂J∂bcvi= J1 − J2 (4.27)Now by using Equations (4.25) and (4.27), where L and R are respectively the leftand right eigenvectors of the right most eigenvalue, we can measure the gradient of the149eigenvalue with respect to the boundary flux weights b. Hence, by expanding the Taylorseries of the new eigenvalue of the Jacobian matrix J , we can find the best boundaryflux configuration for each boundary control volume so that the new projected eigenvalueλn+1 is on the left hand side of the complex plane of eigenvalues:λn+1 = λn + ∂λ∂b−→4b ≤ 0 (4.28)−→4b = − |K| ∂λ∂b(4.29)K = λn(∂λ∂b)2 (4.30)bopt = b+−→4b (4.31)In this way, the best boundary condition configuration at each boundary controlvolume is attained by optimizing the boundary flux weights b.4.3 Results and Discussion4.3.1 Boundary Condition OptimizationIn Chapter (2), while working to improve the stability of finite volume methods on un-structured meshes by optimizing the location of the mesh vertices, it was realized thatnot all the unstable eigenmodes of all the problems were fixable by merely modifyingthe mesh vertices. In the process of stabilizing the eigenvalues by modifying the meshlocally, described in Chapter (2), I came across different types of eigenvectors: spikyeigenvectors and flat-distributed eigenvectors. I noticed that the former category alwaysresults in mesh perturbation vectors that eventually stabilize the corresponding eigen-values. However, the latter category of eigenvectors always produce eigenvalue gradientswhich are noticeably smaller; hence the resulting perturbation vectors are larger than150Degree of FreedomEigenvector Components0 500 1000 1500-1-0.8-0.6-0.4-0.200.20.40.60.81Physics dominatedMesh dominatedFigure 4.2: Different types of eigenvectorsthe local length scale, which inevitably invalidates the mesh. Moreover, applying anyreasonably scaled variant of these vectors does not change the corresponding eigenvaluesignificantly. These findings suggest that this eigenvalue and its corresponding eigenvec-tor are nearly independent of the mesh and are dominated by the physics of the problem.Figure 4.2 shows examples of these two different types of eigenvectors. Moreover, a goodempirical threshold (Eq. 4.32) has been observed to easily distinguish these two types ofeigenvectors.Mesh Dominant:max{∂λj∂ξi}Length Scale > 10−3 (4.32)151where length scale is the minimum length of the incident edges at each vertex and ξi isthe displacement at vertex i. Since boundary conditions always play an important rolein the simulation of physical problem, I investigate whether any of these eigenvalues arerelated to the boundary conditions. To do so, I apply different boundary conditions to seehow sensitive these different eigenvalues and eigenvectors are to the boundary conditions.To better showcase the difference between different eigenmodes, I show here twoeigenvectors for an inviscid Euler problem with Mach 0.8 and an angle of attach of 2◦(Figure (4.3)). I measured the gradient of eigenvalues with respect to the mesh vertexlocations for different boundary conditions at the farfield. Figure (4.3a) clearly showsthat modifying the vertex locations for different boundary conditions result in the samevalues of eigenvalue gradients. In other words, the gradients in this case are independentof the types of boundary conditions, but the eigenvalue is strongly influenced by changesin the mesh (notice the high values of gradients in Figure (4.3a)). On the other hand,changing the boundary conditions at the farfield has noticeably changed the gradientsof another eigenvalue which is more dependent on the physics of the problem (Figure(4.3b)). Additionally, in Figure (4.3b), all the gradients are orders of magnitude smallerthan the ones in Figure (4.3a). This in turn signals that this specific eigenmode is notsensitive to the changes in the mesh vertices. Another important point to notice in thisanalysis is that when stability is mesh related, the larger gradients of the eigenvalues withrespect to the mesh modification are associated with the larger eigenvector components(see Figure (4.3a)). However, this is not the case in Figure (4.3b).4.3.1.1 FarfieldTo investigate the sensitivity of the eigenvalues with respect to the boundary conditionsand to validate the boundary condition stabilization in sub-section (4.2.1), I stabilizeor improve the convergence rates for an inviscid Euler problem. Figure (4.4) shows theboundary condition types at the far field and their best configuration for an Euler problem152Eigenvector Component(d/d) / (Lentgh Scale)10-5 10-4 10-3 10-2 10-1 10010-510-410-310-210-1100101102103104Pressure, 0 = 1.822Characteristics with Vortex, 0 = 1.891Radial Velocity, 0 = 2.052Nonreflecting pressure, 0 = 1.822(a) Boundary condition independent eigenmodeEigenvector Component(d/d) / (Lentgh Scale)10-5 10-4 10-3 10-2 10-1 10010-510-410-310-2Pressure, 1 = -0.0303Characteristics with Vortex, 1 = -0.044Radial Velocity, 1 = 0.1079Nonreflecting pressure, 1 = -0.0303(b) Boundary condition dependent eigenmodeFigure 4.3: Derivative of the eigenvalue with vertex movement for two eigenmodes ofEuler problem with Mach 0.8, and an angle of attack of 2◦ on a NACA 0012 airfoil with2000 control volumes. Vertical axes are the average of the gradients of eigenvalues withrespect to the mesh vertices for all the vertices of each control volume. The gradientsof eigenvalues with respect to mesh vertices are calculated at each vertex. However, toproject them on the eigenvector components (i.e. each unknown at each control volume)the gradients are averaged for each CV. The horizontal axes are the average of the valuesof all components of the eigenvector corresponding to each control volume. Each controlvolume has the same number of entities in the eigenvector as the number of unknowns.153Figure 4.4: Boundary condition configurations: far-field boundary condition optimizationfor an Euler problem with Mach 0.55, and an angle of attack of 2◦ on a NACA 0012 airfoilwith 2000 control volumeson NACA 0012 airfoil with Mach 0.55, and an angle of attack of 2◦. In this problem,four different far field boundary conditions of vortex characteristics, radial velocity [116],pressure, and Dirichlet were used. It can be observed from Figure (4.5) that theseoptimized boundary condition configurations have been able to stabilize an originallyunstable problem or in some cases improve the convergence rate. In all the cases, thebundled boundary conditions are showing better convergence characteristics. After theoptimization process, the boundary conditions are linear averages of each individualboundary condition, where the weights of these linear combinations are the results of theoptimization.154IterationsL 2 Norm of Residual20 40 60 80 10010-810-610-410-2100102PressureRadialVortex CharacteristicsDirichletPressure-RadialVortex C.-PressureVortex C.-RadialPressure-DirichletDirichlet-RadialVortex C.-DirichletFigure 4.5: History of convergence: far-field boundary condition optimization for anEuler problem with Mach 0.55, and an angle of attack of 2◦ on a NACA 0012 airfoil with2000 control volumes155(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.6: Boundary condition optimization between Rigid Walls (RW), and Charac-teristic Walls (CW), for two different flow conditions around a NACA 0012 airfoil with5000 control volumes. The leading edge of the airfoil is displayed.4.3.1.2 WallI implemented three different types of wall boundary conditions namely, Rigid wall,Characteristic wall, and Soft wall, to investigate their stability properties. Using theoptimization scheme described in Section (4.2.1), I apply pairwise optimization of thesewall boundary conditions to find a hybrid boundary condition which shows superiorconvergence and stability properties. The boundary control volumes in Figures (4.6),(4.7), (4.8), (4.9), (4.10), (4.11), (4.12), (4.13) and (4.14) are colored based on the weightsof the boundary flux, where a weight of 0 or 1 means only one type of boundary conditionfor that control volume is used, while a weight between zero and one shows a hybrid type.I also use two flow conditions for a second order Euler problem where one is showing moreinstability issues. These analyses help us better understand how the optimized hybridboundary conditions can improve either stability or convergence for an already stableproblem.For example, Figures (4.6), (4.7) and (4.8) shows the optimized boundary conditionconfiguration between two types of Rigid wall and Characteristic wall for two secondorder Euler problems around a NACA 0012 airfoil with 5000 control volumes. For this156(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.7: Boundary condition optimization between Rigid Walls (RW), and Charac-teristic Walls (CW), for two different flow conditions around a NACA 0012 airfoil with5000 control volumes. The central part of the airfoil is displayed.(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.8: Boundary condition optimization between Rigid Walls (RW), and Charac-teristic Walls (CW), for two different flow conditions around a NACA 0012 airfoil with5000 control volumes. The trailing edge of the airfoil is displayed.157(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.9: Boundary condition optimization between Rigid Walls (RW), and Soft Walls(SW), for two different flow conditions around a NACA 0012 airfoil with 5000 controlvolumes. The leading edge of the airfoil is displayed.optimization problem a lower order Euler solution is used to approximate the steady-state solution and to linearize the Jacobian matrix used in optimization procedure, Eq.(4.28). The initial value of the boundary flux weights are set to 0.5 to reflect an equalcontribution from both boundary types in the beginning. In Figures (4.6a), (4.7a) and(4.8a), the boundary condition configuration for leading, central, and trailing edge of theairfoil respectively for a Mach of 0.81 and angle of attack of 1.25◦, most of the controlvolumes are impartial to the type of boundary condition, while two of the control volumesare strongly leaning towards the Characteristic wall. On the contrary, in Figures (4.6b),(4.7b) and (4.8b) for a flow conditions of Mach of 0.88 and angle of attack of 0◦, almostall the control volumes are strongly biased towards the Characteristic wall conditions.Likewise, the same optimization procedure is applied to get the most stable hybridboundary conditions for Rigid walls, Soft walls (see Figure (4.9), (4.10) and (4.11)) andCharacteristic walls, Soft walls (see Figure (4.12), (4.13) and (4.14)).Figures (4.15) and (4.16) show the convergence history of all the different types ofwall boundary conditions for the two flow conditions studied in this paper. In these twofigures I show the flux residuals approaching the steady-state solution for all single types158(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.10: Boundary condition optimization between Rigid Walls (RW), and Soft Walls(SW), for two different flow conditions around a NACA 0012 airfoil with 5000 controlvolumes. The central part of the airfoil is displayed.(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.11: Boundary condition optimization between Rigid Walls (RW), and Soft Walls(SW), for two different flow conditions around a NACA 0012 airfoil with 5000 controlvolumes. The trailing edge of the airfoil is displayed.159(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.12: Boundary condition optimization between Characteristic Walls (CW), andSoft Walls (SW), for two different flow conditions around a NACA 0012 airfoil with 5000control volumes. The leading edge of the airfoil is displayed.(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.13: Boundary condition optimization between Characteristic Walls (CW), andSoft Walls (SW), for two different flow conditions around a NACA 0012 airfoil with 5000control volumes. The central part of the airfoil is displayed.160(a) M = 0.81, and α = 1.25◦ (b) M = 0.88, and α = 0◦Figure 4.14: Boundary condition optimization between Characteristic Walls (CW), andSoft Walls (SW), for two different flow conditions around a NACA 0012 airfoil with 5000control volumes. The trailing edge of the airfoil is displayed.of boundary conditions and the pairwise optimized ones. Additionally, the convergencehistory of the average of these pairwise boundary type conditions are included as abenchmarking reference. The average here means that both of the boundary conditiontypes are equally contributing to the boundary flux functions (they both have the sameweight). This average is a benchmark as it is the initial condition in optimizing thehybrid boundary condition types. Figure (4.16) shows that this boundary conditionoptimization is able to stabilize an unstable problem. Another important point to noticeis that the softwall boundary condition, which damps out the wave reflections at the wall,has the best convergence properties amongst the three individual boundary conditions.Moreover in both Figures (4.15) and (4.16), the optimized hybrid boundary conditionsare yielding the best convergence rates. These convergence histories are consistent withtheir rightmost eigenvalues (See Table (4.1))To affirm that the steady state solution does not change when using hybrid boundaryconditions, the pressure coefficient around the airfoil for both of the flow problems aresketched in Figures (4.17) and (4.18). It can be seen that the optimized hybrid bound-ary conditions are converging to the same pressure coefficients that the individual wall161IterationsL 2 Norm of Residual20 4010-810-610-410-2100102104 RWCWSWRW-CW (Average)RW-CW (Optimized)RW-SW (Average)RW-SW (Optimized)CW-SW (Average)CW-SW (Optimized)Figure 4.15: L2 norm of flux residual for NACA 0012 with 5000 control volumes withM = 0.81, and α = 1.25◦ and with different wall boundary conditions. Average of theboundary conditions means that both of the boundary conditions used, have the sameweights and are contributing equally. The average is used as a benchmark point.Boundary Conditions λ0M = 0.81, and α = 1.25◦ M = 0.88, and α = 0◦Rigid Wall (RW) -0.0060 0.0280Characteristic Wall (CW) -0.0121 0.0020Soft Wall (SW) -0.0138 -0.0136RW-CW Optimized -0.0130 -0.0141RW-SW Optimized -0.0138 -0.0141CW-SW Optimized -0.0144 -0.0088Table 4.1: Wall boundary conditions and the right most eigenvalues162IterationsL 2 Norm of Residual20 40 60 80 10010-810-610-410-2100102104 RWCWSWRW-CW (Average)RW-CW (Optimized)RW-SW (Average)RW-SW (Optimized)CW-SW (Average)CW-SW (Optimized)Figure 4.16: L2 norm of flux residual for NACA 0012 with 5000 control volumes withM = 0.88, and α = 0◦ and with different wall boundary conditions. Average of theboundary conditions means that both of the boundary conditions used, have the sameweights and are contributing equally. The average is used as a benchmark point.163X/LC p0 0.2 0.4 0.6 0.8-1-0.500.51RWCWSWRW-CW (Optimized)RW-SW (Optimized)CW-SW (Optimized)Figure 4.17: Pressure coefficient along the NACA 0012 airfoil with M = 0.81, andα = 1.25◦with different wall boundary conditionsboundary conditions (i.e., Rigid wall, Characteristic Wall, and Soft wall) are convergingto.Additionally, for all the optimized hybrid wall boundary conditions the pressure fieldsare shown in Figure (4.19). We can see that even though the mesh is not refined aroundthe shock discontinuity, all the hybrid boundary conditions are capturing the same pres-sure field and the shock line around the airfoil. This further confirms that the use of thehybrid boundary conditions does not change the steady-state solution.164X/LC p0 0.2 0.4 0.6 0.8 1-1-0.500.51SWRW-CW (Optimized)RW-SW (Optimized)CW-SW (Optimized)Figure 4.18: Pressure coefficient along the NACA 0012 airfoil with M = 0.88, andα = 0◦with different wall boundary conditions. In this case, Rigid wall and Characteristicwall boundary conditions did not converge, hence they are absent in this analysis.165(a) Optimized boundary with Rigid wall and Characteristic wall(b) Optimized boundary with Rigid wall and Soft wall(c) Optimized boundary with Rigid wall and Soft wallFigure 4.19: Pressure distribution of Euler problem around NACA 0012 airfoil with 5000control volumes and with M = 0.81, and α = 1.25◦166Chapter 5Eigenvalue and Solution Sensitivity5.1 Sensitivity of Eigenvalues to SolutionAs I have typically used surrogate solutions or lower order solutions to approximatethe semi-discrete Jacobian for eigenvalue analysis, it is of paramount importance tounderstand the effect of these approximate solutions on the rightmost eigenvalues of theJacobian matrix. One common way to compute sensitivity is to use adjoint methods (e.g.,see the work done by Giles and Pierce [117]). In an adjoint problem setting, we specifythe output as the rightmost eigenvalues, λ, and the sensitivity parameter in theory canbe any measurable variable. However, since we calculate the sensitivity with respectto different parameters of the discrete solution, such as mesh vertex location, solutionreconstruction, or the type of boundary condition as used in the previous chapters, weuse the same variables, b, for the sensitivity parameter:dλdb= −ΨT ∂<∂b+ ∂λ∂b(5.1)[∂<∂U]TΨ = JTΨ =[∂λ∂U]T(5.2)where in the functional output sensitivity equation, Eq. 5.1, and the adjoint linearproblem, Eq. 5.2, Ψ is the adjoint variable and < is the residual of the primal PDE.Now, the sensitivity of the eigenvalue with respect to the solution can be calculated by:167∂λ∂U=(∂<∂bT)−1 (dλdb− ∂λ∂b)T(5.3)where the derivatives can be calculated by doing finite differences on the mesh andsolution:dλdb= y(J (U (b+ db) , b+ db)− J (U (b) , b)‖db‖)x (5.4)∂λ∂b= y(J (U (b) , b+ db)− J (U (b) , b)‖db‖)x (5.5)∂<∂b= < (U (b) , b+ db)−< (U (b) , b)‖db‖ (5.6)However, calculating the sensitivity of the eigenvalues with respect to the solution usingEq. 5.3 requires calculating multiple derivatives for each control volume in the problem.Instead, I use Eq. 4.25 and the Frechet derivatives to calculate these sensitivity values.In this way, the sensitivity of eigenvalues to the solution can be calculated by substitutingparameters b in Eq. 4.25 with control volume solution averages:∂λ∂U= y ∂J∂Ux (5.7)where the right hand side of Eq. 5.7 is calculated by multiplying the left eigenvector bythe following Frechet derivative:∂J∂Ueˆkx =J (U + εeˆkx)− J (U)ε(5.8)In Eq. 5.7 and 5.8, U is the vector of the control volume solution averages and ε is aninfinitesimal number for differencing. Eq 5.8 is calculated row by row for each controlvolume to break up an otherwise third rank tensor. However, in evaluation of Eq. 5.8,there is no need to recalculate the whole Jacobian matrix, J , for all control volumes in168the mesh every time, as this would result in a O(N2) computational time where N is thenumber of control volumes. Instead, when calculating each row of the right hand side ofEq. 5.7, we determine the specific rows in the Jacobian matrix which need to be updatedafter the solution within a control volume has been perturbed. In other words, once thesolution within control volume α changes, the flux integral at all nearby control volumeswhich have control volume α in their reconstruction stencil change. As the discretizationis constant and does not change when the solution changes, the sparsity of the Jacobianmatrix is known a priori. This means that for identifying the control volumes for whichthe flux integral changes after a solution perturbation at control volume α, we look atnon-zero entries of the column of the Jacobian matrix associated with control volume α.In this way, the evaluation of Eq. 5.8 is of order of O(N).5.2 Results and Discussion5.2.1 Sensitivity of Eigenvalues to SolutionTo present the effect of the solution on the eigenvalues, I developed a methodology tocalculate the sensitivity of each of the eigenvalues to the solution in Section 5.1. Here, Ishowcase this sensitivity for a second order inviscid transonic Euler problem with Machnumber 0.8 and an angle of attack of α = 1.25◦. This problem, with 4200 control volumes,is initialized with farfield free stream flow conditions. It is evolved to steady-state usingpseudo-time stepping with an implicit time advance scheme. Figure 5.1 shows the L2norm of residual for this problem and the three iterations at which I have measured thesensitivity of the eigenvalues with respect to the solution.Figures 5.2, 5.3, 5.4, and 5.5 are showing the gradients of the rightmost eigenvaluewith respect to the four primitive variables, ρ, u, v, p respectively at all control volumes.In addition to the three iterations, as shown in Figure 5.1, at which the sensitivity ismeasured, I am also showing the sensitivity of the rightmost eigenvalue of the same169IterationsL 2 Norm of Residual5 10 15 20 25 30 35 4010-710-510-310-1101103Figure 5.1: The residual history for an Euler problem with Mach 0.8 and an angle ofattack of α = 1.25◦ on NACA 0012 airfoil with 4200 control volumes.1701e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(a) First order solution, λ0 = −0.03011e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(b) At 10th iteration, λ0 = −0.03121e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(c) At 35th iteration, λ0 = −0.03181e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(d) At 42nd iteration, λ0 = −0.0318Figure 5.2: Sensitivity of the rightmost eigenvalue with respect to the density, ρ, at threedifferent stages of convergence (see Figure 5.1).second-order discretization evaluated using the first-order solution as a linearization.This sheds light on the previous Chapters (2), (3) and (4), where I used a lower-ordersteady-state solution to linearize a higher-order discretization (e.g., using a first-ordersteady-state solution as a surrogate solution to the second order-solution for linearizingthe Jacobian).Comparing the sensitivity maps between all four different solutions and for each ofthe primitive variables reveals that both the half converged solution and the lower ordersolution have similar features to the sensitivity of the fully converged solution at nearfield. The rightmost eigenvalues of all the four solutions are relatively close, while theeigenvalues at iteration 35 (with a flux integral residual of order of 10−1) and the steady-1711.0e-071.0e-061.0e-051.0e-041.0e-031.0e-081.0e-02Gradient(a) First order solution, λ0 = −0.03011.0e-071.0e-061.0e-051.0e-041.0e-031.0e-081.0e-02Gradient(b) At 10th iteration, λ0 = −0.03121.0e-071.0e-061.0e-051.0e-041.0e-031.0e-081.0e-02Gradient(c) At 35th iteration, λ0 = −0.03181.0e-071.0e-061.0e-051.0e-041.0e-031.0e-081.0e-02Gradient(d) At 42nd iteration, λ0 = −0.0318Figure 5.3: Sensitivity of the rightmost eigenvalue with respect to the axial velocity, u,at three different stages of convergence (see Figure 5.1).1721e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(a) First order solution, λ0 = −0.03011e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(b) At 10th iteration, λ0 = −0.03121e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(c) At 35th iteration, λ0 = −0.03181e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(d) At 42nd iteration, λ0 = −0.0318Figure 5.4: Sensitivity of the rightmost eigenvalue with respect to the vertical velocity,v, at three different stages of convergence (see Figure 5.1).1731e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(a) First order solution, λ0 = −0.03011e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(b) At 10th iteration, λ0 = −0.03121e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(c) At 35th iteration, λ0 = −0.03181e-71e-61e-50.00010.0011.0e-081.0e-02Gradient(d) At 42nd iteration, λ0 = −0.0318Figure 5.5: Sensitivity of the rightmost eigenvalue with respect to the pressure, p, atthree different stages of convergence (see Figure 5.1).174state solution are exactly identical. However, sensitivity maps at the first iterations ofconvergence (Figures 5.2b, 5.3b, 5.4b, and 5.5b) show high gradients, generally furtheraway from the airfoil and falling short of capturing the important features adjacent tothe airfoil. This is indeed aligned with the intuition that the solution near the airfoil hasnot yet had the chance to evolve from the free stream conditions. The study of eigenvaluesensitivity with respect to the solution further confirms the idea of using either a lower-order or half-converged solution as a surrogate solution for linearizing the Jacobian inour stability algorithms.As stated earlier, when a numerical scheme has instability issues, the problem gen-erally can not converge to any steady state solution. Without any solution, one cannot calculate the Jacobian of the semi-discretized scheme and without measuring theJacobian matrix, one can not calculate the eigenvalues of this Jacobian matrix to pre-dict and improve numerical stability. As stated earlier, to break this deadlock, we use asurrogate solution (such as the steady state solution of a lower order solution), which isas close as possible to the steady state solution of the primary problem, to approximatethe Jacobian matrix of the primary problem. The accuracy of estimating the Jacobianmatrix in this stability analysis is fundamental to the success of this method. Thus, Ipresent (as shown in Table (5.1), (5.2), (5.3) and (5.4)) a comparison of the Jacobianmatrices when measured by their stabilized steady state solution, versus when they areestimated by their lower order steady state solution counterparts. I calculate the normof the difference of the Jacobian matrix measured by either solutions, ∆J = Jh− Jl. Weuse the L1 (the maximum absolute column sum), and the L∞ (the maximum absoluterow sum) induced matrix norms as well as the Frobenius norm, LF , which are defined asfollows:175L1(Jn×n) = ‖J‖1 = max1≤j≤nn∑i=1|Jij| , (5.9)L∞(Jn×n) = ‖J‖∞ = max1≤i≤nn∑j=1|Jij| , (5.10)LF (Jn×n) = ‖J‖F = n∑j=1n∑i=1|Jij|2 12 , (5.11)where Jij is the Jacobian value at the ith row and jth column. The norms of the differencesare then normalized by the norm of the Jacobian of the primary problem:(∆J)L =‖∆J‖‖Jh‖ (5.12)Additionally, I compare the eigensystem of the primary Jacobian with the eigensystemof the approximated Jacobian matrix. I compare the right most eigenvalues of bothmatrices as well as their associated right eigenvectors. To compare the eigenvectors, Iuse the cosine of the angle between them:cos (Rhmax , Rlmax) =Rhmax Rlmax|Rhmax | |Rlmax|(5.13)Table (5.1) shows these three matrix norms as well as the right most eigenvalues of theprimary Jacobian, λhmx , and the right most eigenvalue of the Jacobian when estimated bya lower order steady state solution, λlmax . In Eq. (5.13) and Table (5.1), λhmx and Rhmaxare respectively the rightmost eigenvalue and its corresponding right eigenvector of theJacobian matrix of the original discretization when the stabilized high order steady statesolution is used. It is evident that the differences in the Jacobian matrices measured byboth solution states are small enough to result in very similar eigensystems.Likewise, in optimization of the direction of the solution reconstruction stencil, it canbe observed that the differences in the Jacobian matrices are again small enough to result176Scheme (∆J)L∞ (∆J)LF (∆J)L1 λhmx λlmax cos (Rhmax , Rlmax )O(h2) Euler, M = 0.72, α = 3◦,1500 CV, (Figure (3.13))0.85% 0.57% 0.78% 0.013 0.013 0.987O(h2) Euler, M = 1.21, α = 3◦,1500 CV, (Figure (3.14))1.72% 0.97% 1.63% 0.033 0.036 0.942O(h2) Euler, M = 2, 3500 CV,(Figure (3.18))4.18% 2.39% 4.02% 0.138 0.125 0.902O(h3) Laminar, M = 0.3, α = 2◦,Re = 5000, 2400 CV, (Figure(3.16))0.61% 0.37% 0.58% 0.002 0.002 0.988Table 5.1: Comparison of the Jacobian matrix of the high order scheme when stabilizingthe stencil size using the steady state high order solution versus when using the lower or-der steady state solution to approximate the Jacobian of the high order scheme. (∆J)L∞ ,(∆J)LF , (∆J)L1 are the percentage ratios of the norms of the changes of the Jacobianmatrix to the Jacobian matrix of the primary problem.Scheme (∆J)L∞ (∆J)LF (∆J)L1 λhmx λlmax cos (Rhmax , Rlmax )O(h2) Euler, M = 1.21, α = 3◦,5000 CV, (Figure (3.20))0.82% 0.47% 0.77% 0.087 0.088 0.989O(h2) Laminar, M = 0.42,α = 1◦, Re = 5000, 5000 CV,(Figure (3.16))71.12% 0.61% 1.04% 0.019 0.018 0.969Table 5.2: Comparison of the Jacobian matrix of a higher order scheme when stabilizingthe stencil direction using the primary steady state solution (after stabilization) versuswhen using the lower order steady state solution to approximate the Jacobian of theprimary problem. (∆J)L∞ , (∆J)LF , (∆J)L1 are the percentage ratios of the norms ofthe changes of the Jacobian matrix to the Jacobian matrix of the primary problem.in a similar eigensystem.7In this problem, a half converged solution is used instead of a lower order solution. Thus, λlmax and Rlmax are referring to the halfconverged solution.177Surrogate Solutions (∆J)L∞ (∆J)LF (∆J)L1 λsmax cos (Xpmax , Xsmax )O(h1) Euler steady state solution(Figure (5.2a))1.88% 0.97% 1.79% -0.0301 0.927O(h2) solution at 10th iteration(Figure (5.2b))0.87% 0.48% 0.83% -0.0312 0.972O(h2) solution at 35th iteration(Figure (5.2c))0.08% 0.05% 0.07% -0.0318 0.998Table 5.3: Comparison of the Jacobian matrix of the primary problem of O (h2) Eulerscheme with free stream M = 0.8 and an angle of attack of α = 1.25◦ and 4200 controlvolumes (Figure (5.2d)) when using a surrogate solution (a lower order solution and twohalf-converged solutions of the same primary problem) to approximate the Jacobian ofthe problem. (∆J)L∞ , (∆J)LF , (∆J)L1 are the percentage ratios of the norms of thechanges of the Jacobian matrix to the Jacobian matrix of the primary problem. Therightmost eigenvalue of the primary problem is λpmx = −0.0318.Scheme (∆J)L∞ (∆J)LF (∆J)L1 λpmx λsmax cos (Xpmax , Xsmax )O(h2) Euler, M = 0.8, α = 2◦,2000 CV, (Figure (4.3))1.05% 0.67% 0.98% 1.822 1.735 0.971O(h2) Euler, M = 0.55, α = 2◦,2000 CV, (Figure (4.4))0.62% 0.35% 0.57% -0.0115 -0.012 0.989O(h2) Euler, M = 0.81,α = 1.25◦, 5000 CV, (Figure(4.6a))1.38% 0.76% 1.35% -0.0060 -0.0082 0.932O(h2) Euler, M = 0.88, α = 0◦,5000 CV, (Figure (4.6b))1.21% 0.66% 1.18% 0.0280 0.0220 0.958Table 5.4: Comparison of the Jacobian matrix of the high order scheme when optimizingthe boundary condition types using the steady state high order solution versus whenusing the lower order steady state solution to approximate the Jacobian of the high orderscheme. (∆J)L∞ , (∆J)LF , (∆J)L1 are the percentage ratios of the norms of the changesof the Jacobian matrix to the Jacobian matrix of the primary problem.178Chapter 6Concluding RemarksIn this work, I studied numerical stability of the semi-discrete system of equations arisingfrom finite volume discretizations of the PDE’s governing computational fluid dynamicsproblems at steady state conditions and proposed a new family of approaches to au-tomatically stabilize these schemes in practice. Various contributions of this work arepresented in the following sections:6.1 MeshThe proposed approach, which to the best of my knowledge is the first of its kind, movesthe mesh vertices locally so that the new mesh is more suitable for the PDE of interest.In my quest to improve stability, I exploit the gradients of eigenvalues to determine inwhich direction and by how much the mesh vertices should be moved so that the Jacobianof the semi-discretized set of equations has more stable eigenvalues. The less positivethese eigenvalues are, the more stable the semi-discrete system of equations are. Theresults for linear advection and the nonlinear inviscid Burgers’ and Euler problems showthe applicability and merits of this algorithm. I have also demonstrated that being ableto move vertices on the boundary can significantly improve convergence of the eigenvalueoptimization process and increase the convergence rate on the resulting optimized mesh.Overall, I showed that changing the mesh locally for a subset of control volumes (for lessthan 2% of control volumes in many cases), can stabilize the initially unstable problemswith minimal efforts and computational costs.1796.2 Surrogate SolutionsIn this work, an a priori stability analysis of Euler problems using a surrogate solutionwas introduced. I have demonstrated that an estimate of the steady state solution,such as potential flows for Euler problems or lower order steady state solution of theprimary problem, can be sufficient for diagnosing and improving the numerical stabilityof the steady state solution. Three different solutions of potential flow for differentflow conditions were used to show that the stability can be predicted by fast and cheapsurrogate solutions. Indeed, the key outcome of this work is that if there is a steadystate instability arising from the mesh, mesh vertex locations can be modified based onan eigenanalysis of the semi-discrete Jacobian, which is approximated using data froma surrogate solution. This study in turn paves the way for systematically diagnosingand improving numerical stability of more complicated flows such as the Navier-Stokesequations in advance by incorporating integral boundary layer methods with potentialflow solutions.6.3 ReconstructionI have been able to mathematically and numerically show that there is a correlationbetween the number of cells in the reconstruction stencil and the stability of the semi-discrete Jacobian. For second order advection there is a clear threshold of stability basedon the value of the norm of the reconstruction map, first identified by Haider and hiscolleagues: stable problems usually have a local reconstruction map norms smaller thanone. However, this threshold does not hold for higher order and nonlinear problems.Yet even for these problems, I have shown that a larger stencil size always leads to asmaller norm of the reconstruction map as well as a more stable problem. Additionally,the dependency of the local reconstruction map on the mesh was investigated and anempirical correlation was found between this value and the minimum of relative moments180of the control volumes used in the reconstruction. It was shown that prioritizing thecontrol volumes in the stencil based on this criterion can decrease the value of the normof the reconstruction map and increase the stability, which is the end goal. Secondly,I proposed a more systematic approach for designing the direction and the size of thereconstruction stencil to improve numerical stability. Through a gradient descent andback-propagation algorithm, the weights of the reconstruction stencil are updated towardsone size or direction of the stencil to minimize (or negate) the rightmost eigenvalue. Inthis way the reconstruction scheme will yield an energy stable scheme and will also bestable in practice.6.4 Boundary ConditionsIn this work, I first proposed and studied a new approach to optimize boundary conditionsfor stabilization of finite volume discretization on unstructured meshes. This method su-perimposes different boundary conditions for aerodynamic problems at far-field and solidwalls. It automatically optimizes this mixed boundary condition configuration throughgradients of the right most eigenvalues of the semi-discretized Jacobian. In this way, aboundary condition configuration is designed that produces a stable finite volume methodin terms of steady state convergence without affecting the solution at the steady-state.6.5 Solution Sensitivity and EigenanalysisAdditionally, I have presented a novel methodology to investigate the sensitivities of theeigenvalues with respect to the solution at any point when converging to the steady-state solution. The results show that both the lower-order and half-converged solutionscapture enough variation in the solution to be able to approximate the semi-discretizedJacobian of a higher-order solution at steady state.1816.6 Recommendations for Future WorkIt was observed in Section (2.4) that computing the eigensystem overwhelmingly com-prises majority of the computational time of the stabilization algorithm. Hence, if onecan find a better and faster eigenvalue solver, the entire algorithm can run up to twentytimes faster. Additionally, it is seen that the current eigensolver even after implementingmultiple techniques to improve its performance is extremely slow for much larger prac-tical problems in three dimensions. Calculating more than a handful of eigenvalues isalmost impossible. Therefore, if this research is to be continued to turbulent and 3Dflows, definitely the main concern should be a fast and scalable eigensolver.Additionally, there are two other challenges for extension of the stability analyses to3D problems. First, one has to be able to handle the edge cases such as mesh movementin sliver cells and if any negative volumes are to be created in the process. One mayalso consider applying edge/face swapping in these cases to recover some of the cells withnegative volume. Second, one needs to re-think the optimization scheme. It has been seenin practice that in the presence of several unstable eigenvalues, the optimization problemhas a hard time finding the optimum point as many of the modes will be competingagainst each other. This can get harder as there are more optimization parameters (onemore degree of freedom for each vertex) and many more vertices to worry about.Another recommendation for the future work is to extend the use of medium fidelitysolutions as surrogate solutions to predicting and fixing instabilities for laminar andturbulent flows. This requires correcting the potential solution and coupling it with theintegral boundary layer flow.In Chapter (5), sensitivity of the eigenvalues with respect to the solution was pre-sented. However, this work can be greatly augmented by using Adjoint methods and Eq.(5.1), (5.2) and (5.3) to calculate the total derivative of the eigenvalues of the Jacobianof the semi-discrete system with respect to the parameters of interest (e.g., mesh vertexlocations). In this thesis, I was only able to use the partial gradients of the eigenvalues.182This means that some pieces of information might have been missed. For example, whenthe mesh changes locally, it will cause a change in the solution which in turn will have asecondary effect on the eigenvalues. This was not captured in the analysis used in thisthesis as a fixed solution was assumed when taking the derivatives of the eigenvalues withrespect to mesh vertex coordinates.Most importantly, in the eigenanalysis (Section (1.3.2)) and the stability analysis(Section (1.3.1)) the necessary and sufficient conditions are only true at the same time ifthe Jacobian matrix is symmetric which is not the case in general as explained previously.However, in this thesis only the sufficient condition was used in practice to develop thestability analysis schemes.In the eigenvalue analysis, the sufficient condition is used which states that if theproblem is energy stable, all the eigenvalues are non-positive. In this way, in stabilizingthe numerical schemes I might end up doing more work than is necessary for stabilityby pushing all the eigenvalues to the left half of the complex plane. However, by doingso the energy stability is guaranteed by the sufficient conditions. For example, I mightunnecessarily use this scheme to change the rightmost positive eigenvalues of an alreadystable problem to guarantee energy stability. Yet, by doing this extra work I reached toenergy stability automatically for every problem.Nevertheless, negative-definiteness is defined for general matrices as well. A necessaryand sufficient condition for a general real square matrix J to be negative definite is thatits symmetric part is negative definite. As stated earlier, a general real square matrix Jis negative definite if for all non zero real vectors of x, the following equation holds:xTJx < 0 (6.1)Hence, when the symmetric part of the matrix J is defined as Js = 12(J + JT), thematrix M is said to be negative definite if and only if:183xTJsx < 0 (6.2)Since matrix Js is symmetric, it is negative definite if and only if all its eigenvaluesare negative. Consequently, a general real square matrix M (the Jacobian matrix inthis case) is negative definite (or semi-definite) if and only if all the eigenvalues of itssymmetric part are negative (or non-positive). This can be the focus of a body of researchto extend the energy stability analysis by using the negative definiteness of the symmetricpart of the matrix in the future. Moreover, commercial eigensolvers are generally morerobust and convergent with symmetric matrices which can greatly reduce the cost of thestabilization method in this thesis.184Bibliography[1] R. Zangeneh, C. Ollivier-Gooch, Mesh optimization to improve the stability offinite-volume methods on unstructured meshes, Computers & Fluids 156 (2017)590–601.[2] R. Zangeneh, C. Ollivier-Gooch, A posteriori stability analysis and improvementfor finite volume methods on unstructured meshes., in: 55th AIAA AerospaceSciences Meeting, 2017.[3] R. Zangeneh, C. Ollivier-Gooch, A priori stability analysis of finite volume methodson unstructured meshes, in: 56th AIAA Aerospace Sciences Meeting, 2018.[4] R. Zangeneh, C. F. Ollivier-Gooch, Stability analysis and improvement of thesolution reconstruction for cell-centered finite volume methods on unstructuredmeshes, Journal of Computational Physics in press (2019).[5] R. Zangeneh, C. Ollivier-Gooch, Reconstruction map stability analysis for cellcentered finite volume methods on unstructured meshes., in: 55th AIAA AerospaceSciences Meeting, 2017.[6] T. J. R. Hughes, W. K. Liu, Implicit-explicit finite elements in transient analysis:Stability theory, Journal of Applied Mechanics 45 (1977) 371–374.[7] F. Haider, J.-P. Croisille, B. Courbet, Stability analysis of the cell centered finite-volume MUSCL method on unstructured grids, Numerische Mathematik 113 (2009)555–600.185[8] P. H. Lauritzen, A stability analysis of finite-volume advection schemes permittinglong time steps, American Meteorological Society 135 (2007) 2658–2673.[9] X. Merle, F. Alizard, J. C. Robinet, Finite difference methods for viscous incom-pressible global stability analysis, Computers & Fluids 39 (2010) 911–925.[10] C. Berthon, Stability of the MUSCL schemes for the Euler equations, Communi-cations in Mathematical Sciences 3 (2005) 133–157.[11] D. Ray, P. Chandrashekar, Entropy stable schemes for compressible Euler equa-tions, International Journal of Numerical Analysis and Modeling 4 (2013) 335–352.[12] M. Lukacova-Medvidova, E. Tadmor, On the entropy stability of Roe-type finitevolume methods, in: Proceedings of Symposia in Applied Mathematics, volume 67,2009.[13] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Decon-inck, R. Hartmann, K. Hillewaert, H. Huynh, N. Kroll, G. May, P.-O. Persson,B. van Leer, M. Visbal, High-order CFD methods: Current status and perspective,International Journal for Numerical Methods in Engineering 72 (2013) 811–845.[14] D. W. Zingg, Comparison of high-accuracy finite-difference methods for linear wavepropagation, SIAM Journal on Scientific Computing 22 (2000) 476–502.[15] T. Barth, P. Frederickson, Higher order solution of the Euler equations on unstruc-tured grids using quadratic reconstruction, in: 28th Aerospace Sciences Meeting,1990.[16] A. Nejat, C. Ollivier-Gooch, Effect of discretization order on preconditioning andconvergence of a high-order unstructured Newton-GMRES solver for the Eulerequations, Journal of Computational Physics 227 (2008) 2366–2386.186[17] C. Michalak, C. Ollivier-Gooch, Globalized matrix-explicit Newton-GMRES forthe high-order accurate solution of the Euler equations, Computers & Fluids 39(2010) 1156–1167.[18] M. Dumbser, A. Hidalgo, O. Zanotti, ADER-WENO finite volume schemes withspace-time adaptive mesh refinement, Journal of Computational Physics 268 (2013)359–387.[19] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order ac-curate essentially non-oscillatory schemes, Upwind and high-resolution schemes(1987) 218–290.[20] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calcula-tions: a quasi conservative approach, Journal of Computational Physics 125 (1996)150–160.[21] T. J. Barth, Recent developments in high order k-exact reconstruction on unstruc-tured meshes, AIAA paper 93-0668, 1993.[22] C. F. Ollivier-Gooch, M. Van Altena, A high-order accurate unstructured meshfinite-volume scheme for the advection-diffusion equation, Journal of Computa-tional Physics 181 (2002) 729–752.[23] E. Lee, H. T. Ahn, H. Luo, Cell-centered high-order hyperbolic finite volumemethod for diffusion equation on unstructured grids, Journal of ComputationalPhysics 355 (2018) 464–491.[24] G. Hu, N. Yi, An adaptive finite volume solver for steady Euler equations with non-oscillatory k-exact reconstruction, Journal of Computational Physics 312 (2016)235–251.187[25] P. L. Roe, Approximate Riemann solvers, parameter vectors, and differenceschemes, Journal of Computional Physics 43 (1981) 357–372.[26] H. Nishikawa, Beyond interface gradient: A general principle for constructing diffu-sion schemes, in: In Proceedings of the Fortieth AIAA Fluid Dynamics Conferenceand Exhibit, 2010, pp. 5082–5093.[27] A. H. Stroud, D. Secrest, Gaussian Quadrature Formulas, Prentice-Hall, EnglewoodCliffs, N.J., 1966.[28] M. Blanco, D. W. Zingg, Fast Newton-Krylov method for unstructured grids.,American Institute of Aeronautics and Astronautics Journal 35 (1998) 607–612.[29] M. Ceze, K. J. Fidkowski, Constrained pseudo-transient continuation., Interna-tional Journal for Numerical Methods in Engineering 102 (2015) 1683–1703.[30] G. O. George, M. A. Hyman, S. Kaplan, A study of the numerical solution of partialdifferential equations, Journal of Mathematics and Physics 29 (1950) 223–251.[31] P. D. Lax, R. D. Richtmyer, Survey of the stability of linear finite differenceequations, Communications on Pure and Applied Mathematics 9 (1956) 267–293.[32] J. Crank, P. Nicolson., A practical method for numerical evaluation of solutions ofpartial differential equations of the heat-conduction type, Mathematical Proceed-ings of the Cambridge Philosophical Society 43 (1947) 50–67.[33] R. Courant, K. Friedrichs, H. Lewy, Uber die partiellen differenzengleichungen dermathematischen physik, Mathematische annalen 100 (1928) 32–74.[34] S. C. Reddy, L. N. Trefethen, Stability of the method of lines, Numerische Math-ematik 62 (1992) 235–267.[35] B. Gustafsson, The Godunov-Ryabenkii condition: The beginning of a new stabilitytheory, in: Godunov Methods, Springer, 2001, pp. 425–443.188[36] C. Hirsch, Numerical computation of internal and external flows: The fundamentalsof computational fluid dynamics, Elsevier, 2007.[37] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time dependent problems and differencemethods, John Wiley & Sons, 1995.[38] C. M. Dafermos, Hyperbolic conservation laws in continuum physics, Springer-Verlag, 2010.[39] T. Chen, C.-W. Shu, Entropy stable high order discontinuous Galerkin methodswith suitable quadrature rules for hyperbolic conservation laws, Journal of Com-putational Physics 345 (2017) 427–461.[40] C. W. Gear, Numerical solution of ordinary differential equations: is there anythingleft to do?, SIAM review 23 (1981) 10–24.[41] R. D. Richtmyer, E. H. Dill, Difference methods for initial-value problems, PhysicsToday 50 (1958).[42] B. Despres, Lax theorem and finite volume schemes, Mathematics of Computation247 (2004) 1203–1234.[43] F. Boyer, Convergence analysis of the upwind finite volume scheme for generaltransport problems, Finite Volumes for Complex Applications VI Problems &Perspectives (2011) 155–163.[44] M. Ben-Artzi, J. Li, Consistency and convergence of finite volume approximationsto nonlinear hyperbolic balance laws, arXiv preprint arXiv:1902.09047 (2019).[45] M. Ndjinga, Weak convergence of nonlinear finite volume schemes for linear hyper-bolic systems, Finite Volumes for Complex Applications VII-Methods and Theo-retical Aspects (2014) 411–419.189[46] O. O’Reilly, T. Lundquist, E. M. Dunham, J. Nordstrom, Energy stable andhigh-order accurate finite difference methods on staggered grids, Journal of Com-putational Physics 346 (2017) 572–589.[47] B. C. Vermeire, P. E. Vincent, On the properties of energy stable flux reconstructionschemes for implicit large eddy simulation, Journal of Computational Physics 327(2016) 368–388.[48] R. Horn, C. Johnson, Topics in Matrix Analysis, Cambridge Universal Press, 1991.[49] C. Y. Zhang, S. Luo, A. Huang, J. Lu, The eigenvalue distribution of block di-agonally dominant matrices and block hâ matrices, Electronic Journal of LinearAlgebra 20 (2010) 621–639.[50] G. H. Golub, C. F. V. Loan, Matrix Computations, The Johns Hopkins UniversityPress, Baltimore, 1983.[51] Y. Saad, Numerical Methods for Large Eigenvalue Problems, second edition ed.,Manchester University Press, 1992.[52] K. Meerbergen, D. Roose, Matrix transformation for computing rightmost eigen-value of large sparse non-symmetric eigenvalue problems, IMA Journal of NumericalAnalysis 16 (1996) 297–346.[53] C. Campos, J. E. Roman, E. Romero, A. Tomas, SLEPc homepage,http://http://www.grycap.upv.es/slepc/, 2011.[54] B. Smith, S. Balay, J. Brown, L. Dalcin, D. Karpeev, M. Knepley, H. Zhang, PETSchome page, http://www.mcs.anl.gov/petsc, 2011.[55] J. L. Hess, A. M. Smith, Calculation of non-lifting potential flow about arbitrarythree-dimensional bodies, Technical Report ES-40622, Douglas Aircraft Co LongBeach CA, 1962.190[56] R. L. Carmichael, L. L. Erickson, PAN AIR-A higher order panel method for pre-dicting subsonic or supersonic linear potential flows about arbitrary configurations,in: AIAA 14th Fluid and Plasma Dynamics Conference, 1981.[57] F. T. Johnson, R. M. James, J. E. Bussoletti, A. C. Woo, A transonic rectangulargrid embedded panel method, in: (AIAA/ASME) 3rd Joint Thermophysics, Fluids,Plasma and Heat Transfer Conference, 1982.[58] L. L. Erickson, S. M. Strande, A theoretical basis for extending surface-panelingmethods to transonic flow, AIAA Journal 23 (1985) 1860–1867.[59] B. Oskam, Transonic panel method for the full potential equation applied to mul-ticomponent airfoils., AIAA Journal 23 (1985) 1327–1334.[60] F. E. Ehlers, F. T. Johnson, P. E. Rubbert, A higher order panel method forlinearized supersonic flow., in: AIAA 9th Fluid and Plasma Dynamics Conference,1976.[61] P. K. Baruah, J. E. Bussoletti, D. T. Chiang, W. A. Massena, F. D. Nelson, D. J.Furdon, K. Tsurusaki, PAN AIR: A computer program for predicting subsonic orsupersonic linear potential flows about arbitrary configurations using a higher orderpanel method. Volume 4: Maintenance document (version 1.1), Technical Report,NASA-CR-3254, 1981.[62] F. T. Johnson, S. S. Samant, M. B. Bieterman, R. G. Melvin, D. P. Young, J. E.Bussoletti, C. L. Hilmes, TranAir: A full-potential, solution-adaptive, rectangulargrid code for predicting subsonic, transonic, and supersonic flows about arbitraryconfigurations. Theory document, Technical Report, NASA-CR-4348, 1992.[63] A. Jameson, D. A. Caughey, A finite volume method for transonic potential flowcalculations., AIAA Journal (1977).191[64] G. Volpe, A. Jameson, Transonic potential flow calculations by two artificial densitymethods, AIAA Journal 26 (1988) 425–429.[65] X. C. Cai, M. Paraschivoiu, M. Sarkis, An explicit multi-model compressible flowformulation based on the full potential equation and the Euler equations on 3Dunstructured meshes., in: Eleventh International Conference of Domain Decompo-sition Methods, 1999.[66] V. Shankar, V. Ide, H. Gorski, S. Osher, A fast, time-accurate, unsteady fullpotential scheme, AIAA Journal 25 (1987) 230–238.[67] A. Parrinello, P. Mantegazza, Independent two-fields solution for full-potentialunsteady transonic flows, AIAA Journal 48 (2010).[68] A. Parrinello, P. Mantegazza, Improvements and extensions to a full-potentialformulation based on independent fields, AIAA Journal 50 (2012) 571–580.[69] M. C. Galbraith, S. R. Allmaras, R. Haimes, Full potential revisited: A mediumfidelity aerodynamic analysis tool, in: 55th AIAA Aerospace Sciences Meeting,2017.[70] S. C. Reddy, L. N. Trefethen, LAX-stability of fully discrete spectral methods viastability regions and pseudo-eigenvalues, Computer Methods In Applied Mechanicsand Engineering 80 (1990) 147–164.[71] A. Crivellini, F. Bassi, An implicit matrix-free discontinuous Galerkin solver forviscous and turbulent aerodynamic simulations, Computers & Fluids 50 (2011)81–93.[72] H. Nishikawa, P. L. Roe, Third-order active-flux scheme for advection diffusion:Hyperbolic diffusion, boundary condition, and Newton solver, Computers & Fluids125 (2016) 71–81.192[73] D. Zingg, S. De Rango, M. Nemec, T. Pulliam, Comparison of several spatialdiscretizations for the Navier-Stokes equations, Journal of Computational Physics160 (2000) 683–704.[74] C. Huang, L. L. Chen, A new adaptively central-upwind sixth-order WENO scheme,Journal of Computational Physics 357 (2018) 1–15.[75] C. Huang, F. Xiao, T. Arbogast, Fifth order multi-moment WENO schemes forhyperbolic conservation laws, Journal of Computational Physics 64 (2015) 477–507.[76] J. Zhu, C.-W. Shu, A new type of multi-resolution WENO schemes with increas-ingly higher order of accuracy, Journal of Computational Physics 375 (2018) 659–683.[77] S. De Rango, D. W. Zingg, Higher-order spatial discretization for turbulent aero-dynamic computations, AIAA Journal 39 (2001) 1296–1304.[78] A. Jalali, C. F. Ollivier-Gooch, Higher-order unstructured finite volume methodsfor turbulent aerodynamic flows, in: 22nd AIAA Computational Fluid DynamicsConference, AIAA Aviation, 2015, pp. 2284–2297.[79] Q. Wang, Y. X. Ren, J. Pan, W. Li, Compact high order finite volume methodon unstructured grids III: Variational reconstruction, Journal of ComputationalPhysics 337 (2017) 1–26.[80] R. K. Agarwal, D. W. Halt, A compact high-order unstructured grids method forthe solution of Euler equations, International Journal for Numerical Methods inFluids 31 (1999) 121–147.[81] K. H. Kim, C. Kim, Accurate, efficient and monotonic numerical methods formulti-dimensional compressible flows, Journal of Computational Physics 208 (2005)570–615.193[82] C. L. Touze, A. Murrone, H. Guillard, Multislope MUSCL method for generalunstructured meshes, Journal of Computational Physics (2014) 44–99.[83] Z. Sun, S. Inaba, F. Xiao, Boundary Variation Diminishing (BVD) reconstruction:A new approach to improve Godunov schemes, Journal of Computational Physics322 (2016) 309–325.[84] O. Friedrich, Weighted essentially non-oscillatory schemes for the interpolation ofmean values on unstructured grids, Journal of Computational Physics 144 (1998)194–212.[85] Z. Du, J. Li, A Hermite WENO reconstruction for fourth order temporal accurateschemes based on the GRP solver for hyperbolic conservation laws, Journal ofComputational Physics 355 (2018) 385–396.[86] H. Luo, Y. Xia, S. Li, R. Nourgaliev, C. Cai, A Hermite WENO reconstruction-based discontinuous Galerkin method for the Euler equations on tetrahedral grids,Journal of Computational Physics 231 (2012) 5489–5503.[87] H. Luo, Y. Xia, S. Spiegel, R. Nourgaliev, Z. Jiang, A reconstructed discontinuousGalerkin method based on a hierarchical WENO reconstruction for compressibleflows on tetrahedral grids, Journal of Computational Physics 236 (2013) 477–492.[88] R. Kumar, P. Chandrashekar, Simple smoothness indicator and multi-level adap-tive order WENO scheme for hyperbolic conservation laws, Journal of Computa-tional Physics 375 (2018) 1059–1090.[89] B. S. van Lith, J. H. T. T. Boonkkamp, W. L. IJzerman, Embedded WENO: Adesign strategy to improve existing WENO schemes, Journal of ComputationalPhysics 330 (2017) 529–549.194[90] H. Jackson, On the eigenvalues of the ADER-WENO galerkin predictor, Journalof Computational Physics 333 (2017) 409–413.[91] J. Zhu, J. Qiu, A new third order finite volume weighted essentially non-oscillatoryscheme on tetrahedral meshes, Journal of Computational Physics 349 (2017) 220–232.[92] F. Haider, J.-P. Croisille, B. Courbet, Stability of the MUSCL method on generalunstrucutred grids for applications to compressible flows, in: 5th ICCFD, 2008.[93] S. Gorgizadeh, T. Flisgen, U. van Rienen, Eigenmode computation of cavitieswith perturbed geometry using matrix perturbation methods applied on generalizedeigenvalue problems, Journal of Computational Physics 364 (2018) 347–364.[94] C. Hu, C.-W. Shu, Weighted essentially non-oscillatory schemes on triangularmeshes, Journal of Computational Physics 150 (1999) 97–127.[95] D. Zhang, Q. Shangguan, Y. Wang, An easy-to-use boundary condition in dissipa-tive particle dynamics system., Computers & Fluids 166 (2018) 117–122.[96] N. Odier, M. Sanjose, L. Gicquel, T. Poinsot, S. Moreau, F. Duchaine, A charac-teristic inlet boundary condition for compressible, turbulent, multispecies turbo-machinery flows., Computers & Fluids 178 (2019) 41–55.[97] M. Darwish, L. Mangani, F. Moukalled, Implicit boundary conditions for coupledsolvers, Computers & Fluids 168 (2018) 54–66.[98] F. Haider, B. Courbet, J.-P. Croisille, A high-order interpolation for the finitevolume method: The coupled least squares reconstruction, Computers & Fluids176 (2018) 20–39.[99] E. Constant, J. Favier, M. M, P. Meliga, E. Serre, An immersed boundary methodin OpenFOAM: Verification and validation, Computers & Fluids 157 (2017) 55–72.195[100] C. Brehm, O. Browne, N. Ashton, Towards a viscous wall model for immersedboundary methods, in: 56th AIAA Aerospace Sciences Meeting, 2018.[101] J. L. Thomas, M. D. Salas, Far-field boundary conditions for trransonic liftingsolutions to the Euler equations, AIAA Journal 24 (1986) 1074–1080.[102] K. Mazaheri, Numerical Wave Propagation and Steady-State Solutions, Ph.D. the-sis, University of Michigan, 1992.[103] A. Jalaali, C. Ollivier-Gooch, An hp-adaptive unstructured finite volume solverfor compressible flows, International Journal for Numerical Methods in Fluids 85(2017) 563–582.[104] L. T. Diosady, D. L. Darmofal, Preconditioning methods for discontinuous Galerkinsolutions of the Navier-Stokes equations, Journal of Computational Physics 228(2009) 3917–3935.[105] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM Journalon Scientific Computing 14 (1993) 461–469.[106] P. Wong, D. W. Zingg, Three-dimensional aerodynamic computations on unstruc-tured grids using a Newton-Krylov approach, Computers and Fluids 37 (2008)107–120.[107] K. J. Fidkowski, T. A. Oliver, J. Lu, D. L. Darmofal, p-Multigrid solution ofhigh-order discontinuous Galerkin discretizations of the compressible Navier-Stokesequations, Journal of Computational Physics 207 (2005) 92–113.[108] D. J. Mavriplis, Directional agglomeration multigrid techniques for high-Reynolds-number viscous flows, AIAA Journal 37 (1999) 1222–1230.196[109] T. Okusanya, D. Darmofal, J. Peraire, Algebraic multigrid for stabilized finite ele-ment discretizations of the Navier-Stokes equations, Computer Methods in AppliedMechanics and Engineering 193 (2004) 3667–3686.[110] K. Shahbazi, D. J. Mavriplis, N. K. Burgess, Multigrid algorithms for high-orderdiscontinuous Galerkin discretizations of the compressible Navier-Stokes equations,Journal of Computational Physics 228 (2009) 7917–7940.[111] M. Wallraff, R. Hartmann, T. Leicht, N. Kroll, C. Hirsch, F. Bassi, C. Johnston,K. Hillewaert, Multigrid solver algorithms for DGmethods and applications to aero-dynamic flows, Industrialization of High-Order Methods-A Top-Down Approach128 (2015) 53–178.[112] S. V. Tsynkov, Numerical solution of problems on unbounded domains. A review.,Applied Numerical Mathematics 27 (1998) 465–532.[113] B. Engquist, A. Majda, Absorbing boundary conditions for numerical simulationof waves, in: Proceedings of the National Academy of Sciences, volume 74, 1977,pp. 1765–1766.[114] D. H. Rudy, J. C. Strikwerda, A nonreflecting outflow boundary condition forsubsonic Navier-Stokes calculations, Journal of Computational Physics 36 (1980)55–70.[115] M. B. Giles, Nonreflecting boundary conditions for Euler equation calculations,AIAA Journal 28 (1990) 2050–2058.[116] S. R. Allmaras, V. Venkatakrishnan, F. T. Johnson, Far-field boundary conditionsfor 2-D airfoils, in: 17th AIAA Computational Fluid Dynamics Conference, 2005.[117] M. B. Giles, N. A. Pierce, An introduction to the adjoint approach to design, Flow,Turbulence and Combustion 65 (2000) 393–415.197AppendixAs explained in Section (1.3.2.2), iterative methods such as Krylov’s subspace methodsare used in this thesis to calculate a small subset of the eigenvalues of the Jacobianmatrix. To have these iterative methods to converge for the rightmost eigenvalues, Iused various numerical techniques such as spectral transformations and spectrum slicing,and subspace initialization.The first technique and the most effective one used in this thesis is spectral transfor-mations. Spectral transformations are operations that map the eigenvalues of a matrixto another space while the eigenvectors are kept intact. These transformations are usedmainly for three different purposes: first, to compute the eigenvalues in a specific part ofthe spectrum of the eigensystem; second, to expedite the convergence especially when theeigenvalues are closely located near each other; and third, when the matrix is singular. Ihave used two spectral transformations, the Shift-and-Invert and Cayley transformation,specifically for the first two reasons as it was seen in practice that the rightmost eigenval-ues of the Jacobian of the semi-discrete system of finite volume methods on unstructuredmeshes are normally lumped near zero. The techniques used to calculate the rightmosteigenvalues are described in the subsequent sections. However, I will briefly explain howthe Krylov subspace methods are used in this thesis to iteratively find the rightmosteigenvalues.The Jacobian matrices arising from finite volume discretization of compressible flowson unstructured meshes are very large and sparse. To facilitate a cost-efficient method ofcomputation for eigenvalues, and to solve the eigenproblems of Eq. (1.66) and (1.67), thecomputational domain is generally projected to a low dimensional subspace. Assuming198that the orthogonal basis of this subspace is Vm = {v1, v2, · · · , vm}, the eigenproblem ofJR = λR, where λ is the eigenvalue and R is the right eigenvector, is transformed to:V TmJVmVTmR = V TmRλ (6.3)where by considering the new eigenvalues and eigenvectors to be λ˜ = λ and R˜ = V TmRrespectively, the eingeproblem in the lower dimension is:BR˜ = λ˜R˜ (6.4)where B = V TmJVm. This is known as the Rayleigh-Ritz method. In order to solve theeigensystem through this method, first a subspace Vm is found; second, the matrix Bis calculated, and finally the eigen system of Eq. (6.4) is solved in a lower dimensionalsubspace of size m. With this implementation of the Vm subspace, one is only capableof finding m eigenvalues at most. Given this use of subspace projection, most of theiterative eigensolvers differ only in their choice of the subspace. It is also known thatthese projections work best if they mimic a subspace of the original matrix J . For thisreason, the Krylov subspace, which is shown in Eq. (6.5) with x0 as an initial vector, ismuch more effective than a random initizalization.κm ={x0, Jx0, J2x0, · · · , Jm−1x0}(6.5)Using the Krylov subspace, Eq. (6.5), the projected eigensolver converges ratherquickly to the projected eigensystem, provided that the initial vector x0 is as much aspossible in the same direction as the eigenvectors to be sought. This rarely occurs inpractice and thus the Krylov methods iteratively solve this projected eigenproblem byusing larger subspaces in κm. This is called the m-step Arnoldi method. However, thiscan lead to the growth of memory footprint and CPU time. Thus, another method whichaugments the Arnoldi method is to stop the iterations after m steps and restart the entire199process with a new calculated initial vector x0, which is calculated from the previouslycalculated eigenvectors. This method, called Krylov-Schur, is the method of choice inthis thesis. The next sections describe the numerical techniques used to expedite theconvergence of this method.Cayley TransformationThis transformation, unlike the Rayleigh-Ritz transformation described earlier, maps theeigensystem to a new space where eigenvalues are different, while the eigenvectors remainthe same. As explained earlier, one of the issues with the eigensystems arising from thefinite volume discretizations used in this thesis was that the eigenvalues are located quitenear to each other. Iterative eigensolvers such as Krylov notoriously diverge when lookingfor eigenvalues in this packed space. To remedy this, I first use a spectral transformationcalled the Cayley transformation. In this transformation, the eigenproblem of JR = λRis transformed to:(J − σI)−1 (J + υI) R˜ = λ˜R˜ (6.6)where σ is called the shift parameter and υ is called the anti-shift parameter. Thetransformed eigensystem is:R˜ = R (6.7)λ˜ = (λ+ υ)λ− σ (6.8)As mentioned earlier the eigenvector does not change; the eigenvalues are mapped toanother space where the eigenvalues closer to σ will be far from zero in a sparse space andthus easier to find. However, eigenvalues closer to −υ are moved near the origin and thusharder to find. Additionally, this transformation has the most interesting property, where200the eigenvalues to the left of the center-line of σ, −υ (i.e. Re (λ) ≤ σ−υ2 ) are transformedto an area within a unit circle in the new space (harder to find), while the eigenvaluesto the right of this center-line are moved to the outer side of the unit circle with centerof origin in the new space. For this reason, I first use the Cayley transformation withequal shifts and anti-shifts to find the positive rightmost eigenvalues. This means thatby choosing σ = υ, I can more easily find eigenvalues that are positive. However, if noeigenvalues were found using this transformation, I resort to a second transformationwhich is called Shift-and-Invert.Shift-and-Invert TransformationIn this transformation, the eigenproblem of JR = λR is transformed to:(J − σI)−1 R˜ = λ˜R˜ (6.9)where σ, called the shift parameter, simply shifts the origin from zero to the σ value. Theinversion transformation maps the eigenvalues that are closer to σ far from the origin toa sparse region in the mapped space. The transformed eigensystem is:R˜ = R (6.10)λ˜ = 1λ− σ (6.11)This transformation is extremely useful when eigenvalues are close to a target valueare of interest. For this reason, I use this transformation to find the eigenvalues that areclose to zero. To use this transformation effectively, I use a technique called spectrumslicing.201Spectrum slicingSince there is no prior knowledge of the shifts when applying the Shift-and-Invert trans-formation above, the Shift-and-Invert transformation is done in multi-stages. In practice,it has been observed that the eigenvalues of interest from the Jacobian matrices beingstudied, the rightmost eigenvalues, almost always lie in a small range (for example be-tween −0.1 or −0.05 to 1 or 5) close to zero. For this reason, I apply the Shift-and-Inverttransformation in multiple eigensolves for smaller subsets within this range. This helpsin finding the eigenvalues lying in these smaller range better as any eigenvalue in thissmall range will be transformed to a sparse region where the eigensolver is capable offinding it.Subspace InitializationAs explained earlier, if the initial vector x0 in the Krylov subspace is far from the eigen-vector of the eigensystem to be solved, all the variants of the Krylov methods will havehard time converging to a solution. Additionally, it was mentioned in Section (1.3.2.1)that to solve for the left eigenvectors of the Jacobian matrix which are required in thestability process, the right eigenvectors of the transposed Jacobian is solved. This meansthat at every level, the eigensystem has to be solved twice. However, the second eigen-system can be solved much more efficiently and quickly if one initializes the initial vectorx0 with the subspaces of the first eigensystem. Thus, after solving the first eigensystem,the initial vector x0 of the second eigensystem is initialized by the subspace of the firsteigensystem. This has proved to make the second eigensolver drastically faster and cancut the number of iterations required in half.202Command Line InterfaceIn this section, the command line interface and its different options are described. Thecommand line interface provided by the application is consist of an executable which canbe used through a Linux shell. Different options are provided to the executable to rundifferent case scenarios.To execute the eigenanalysis stabilization, use the following command line.$ . / EigenSo lver −f <mesh f i l e > <f low and phys i c s opt ions><Eigenana ly s i s opt ions> −ana lys i s_type <analys i s_types><Mesh , Reconstruct ion , Boundary opt imiza t i on opt ions>To execute the local reconstruction map analysis, use the following command line# ./ reconMap −mesh <mesh f i l e > −d <dimensions>−a <order o f accuracy> −mesh_type <mesh type> <opt ions>All the other options such as the ones used for reordering of the reconstruction, meshoptimization, solution reconstruction and boundary condition optimization are describedin the following sections.Local Reconstruction Map• To specify the number of control volumes in the stencil, use−ReconS <number o f c on t r o l volumes>• To re-order the control volumes in the stencil, use−Reorder_Stenc i lEigenanalysisFirstly, to use a direct eigensolver such as QR or an iterative eigensolver such as Krylov,a macro in setPetscOptions() need. to be changed.203de f i n e LEFT_VECS• Afterwards, the type of the eigensolver can be specified by setting a program optionas follows−eps_type <lapack / kry lovschur /gd/ jd ( d e f au l t i s k ry lovschur )>• The number of the rightmost eigenvalues can be specified by:−eps_nev <# of e igenva lue s>• The number of the largest modes of the eigenvalues to use−nm <# of components>• You can also change the linear solver in the eigensystem−stT <opt ions>j : d i r e c t l i n e a r s o l v e r with ksp p r e cond i t i on e rd : d i r e c t l i n e a r s o l v e r with ILU pr e cond i t i on e ri : GMRES s o l v e r with ILU pr e cond i t i on e rb : GMRES s o l v e r with block Jacobi p r e c ond i t i on e r<de fau l t >: d i r e c t l i n e a r s o l v e r• Shift and antishift are automatically applied based a range of values that are foundto work best for the inviscid schemes. The are hard-coded in the SolveEigenSys-temSubspace function.• All of the eigensystem configurations are applied in the folllowing subroutines:– setEigenSolver– solveeigenSystem– solveTransposeEigenSystem204– solveEigenSystemSubspace– setSTSpace– setRGSpace– setInitialSpaceForTrans– setInitialSpaceMesh Optimization• Starting a mesh optimization analysis−ana lys i s_type mesh_opt• Specifying the surrogate solution−su r rogate_so l <opt ions>lower : i t uses a lower order s o l u t i o npo t e n t i a l : i t uses a f u l l p o t e n t i a l s o l u t i o npanel : i t uses a panel s o l u t i o n• Allowing boundary condition movement−B < f i l e . bdry>Solution Reconstruction• Optimizing the solution reconstruction directions−ana lys i s_type s t enc i l_d i r e c t i on_opt −wide_stenc i l_r−ReconS <s t en c i l_ s i z e >• Optimizing the solution reconstruction size−ana lys i s_type s t enc i l_s i z e_opt −ReconS <s t en c i l_ s i z e >205Boundary Condition Optimization• Optimizing the farfield boundary conditions−ana lys i s_type far_bdry_opt −Bdry_Cond <bdry_options>bdry_options :1 ( d e f au l t ) : Pres sure2 : Radial v e l o c i t y3 : D i r i c h l e t boundary cond i t i on s4 : NonRef lect ing out f low boundary cond i t i on5 : Cha r a c t e r i s t i c s with vortex c o r r e c t i o n s===============Hybrids======================12 : Pressure−Radial13 : Pressure−D i r i c h l e t15 : Pressure−Cha r a c t e r i s t i c s with vortex23 : Radial−D i r i c h l e t25 : Radial−Cha r a c t e r i s t i c s with vortex35 : D i r i c h l e t−Cha r a c t e r i s t i c s with vortex• Optimizing the solid surfaces boundary conditions−ana lys i s_type wall_bdry_opt −Wall_Bdry_Cond <bdry_options>bdry_options :1 ( d e f au l t ) : Rigid wa l l2 : Cha r a c t e r i s t i c s wa l l5 : So f t wa l l===============Hybrids======================12 : Rigid−Cha r a c t e r i s t i c s15 : Rigid−So f t25 : Cha r a c t e r i s t i c s−So f t206Gradients of Eigenvalues w.r.t. Solution• Calculating the gradients of the eigenvalues with respect to the solution.−ana lys i s_type e igen_sol_sens207
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical stability analyses and improvement of cell-centered...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical stability analyses and improvement of cell-centered finite volume methods on unstructured meshes Zangeneh, Reza 2019
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Numerical stability analyses and improvement of cell-centered finite volume methods on unstructured meshes |
Creator |
Zangeneh, Reza |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | The purpose of this thesis is to develop a framework in which one can detect and automatically improve the numerical stability of cell-centered finite-volume calculations on unstructured meshes through optimization schemes that modify the mesh, the solution reconstruction, or the boundary conditions. In this process, eigenanalysis and the gradients of the eigenvalues with respect to different parameters are used to ensure energy stability of the system, consequently resulting in convergence. First, gradients of eigenvalues with respect to the local changes in the mesh are calculated to find directions and magnitudes of mesh movements that will make the Jacobian of a semi-discrete system of equations negative semi-definite. These mesh movement vectors are used to modify the mesh locally. The numerical results show that the proposed methods are able to locate the problematic parts of the mesh responsible for instabilities for several physical problems and to correct the instability. Secondly, I develop a mathematical method, introduced by Haider et al., to measure the stability impact of the reconstruction for high order and nonlinear problems, regardless of the solution. Second order and third order accurate advection and Burgers and Euler problems are used to present detailed practical results and discussion around the use of the local reconstruction map for stability analysis. This method shows that increasing the stencil size will lead to more stable problems. An empirical study is performed which sheds light on connections between the mesh properties and the stability of the reconstruction. I also propose a systematic approach to optimize both the shape and the size of the reconstruction stencil for better numerical stability through eigenvalue analysis. In this approach, one can directly optimize the solution reconstruction stencil for every control volume to obtain better numerical stability and convergence properties for steady state problems. The rightmost eigenpairs of the spatially discretized system of equations are used to obtain an optimized boundary condition which will ensure the energy stability of the system. Lastly, the sensitivity of the rightmost eigenvalues to the solution is measured to investigate the effect of using surrogate solutions for the purpose of linearizing the semi-discretized Jacobian. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-10-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0383237 |
URI | http://hdl.handle.net/2429/71814 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_november_zangeneh_reza.pdf [ 18.99MB ]
- Metadata
- JSON: 24-1.0383237.json
- JSON-LD: 24-1.0383237-ld.json
- RDF/XML (Pretty): 24-1.0383237-rdf.xml
- RDF/JSON: 24-1.0383237-rdf.json
- Turtle: 24-1.0383237-turtle.txt
- N-Triples: 24-1.0383237-rdf-ntriples.txt
- Original Record: 24-1.0383237-source.json
- Full Text
- 24-1.0383237-fulltext.txt
- Citation
- 24-1.0383237.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0383237/manifest