Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Efficient high-order accurate unstructured finite-volume algorithms for viscous and inviscid compressible… Michalak, Christopher 2009

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

Item Metadata

Download

Media
24-ubc_2009_spring_michalak_christopher.pdf [ 3.38MB ]
Metadata
JSON: 24-1.0067108.json
JSON-LD: 24-1.0067108-ld.json
RDF/XML (Pretty): 24-1.0067108-rdf.xml
RDF/JSON: 24-1.0067108-rdf.json
Turtle: 24-1.0067108-turtle.txt
N-Triples: 24-1.0067108-rdf-ntriples.txt
Original Record: 24-1.0067108-source.json
Full Text
24-1.0067108-fulltext.txt
Citation
24-1.0067108.ris

Full Text

Efficient High-Order Accurate Unstructured Finite-Volume Algorithms for Viscous and Inviscid Compressible Flows by Christopher Michalak  B.A.Sc., The University of British Columbia, 2002  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April, 2009 c Christopher Michalak 2009  Abstract High-order accurate methods have the potential to dramatically reduce the computational time needed for aerodynamics simulations. This thesis studies the discretization and efficient convergence to steady state of the high-order accurate finite-volume method applied to the simplified problem of inviscid and laminar viscous two-dimensional flow equations. Each of the three manuscript chapters addresses a specific problem or limitation previously experienced with these schemes. The first manuscript addresses the absence of a method to maintain monotonicity of the solution at discontinuities while maintaining high-order accuracy in smooth regions. To resolve this, a slope limiter is carefully developed which meets these requirements while also maintaining the good convergence properties and computational efficiency of the least-squares reconstruction scheme. The second manuscript addresses the relatively poor convergence properties of Newton-GMRES methods applied to high-order accurate schemes. The globalization of the Newton method is improved through the use of an adaptive local timestep and of a line search algorithm. The poor convergence of the linear solver is improved through the efficient assembly of the exact flux Jacobian for use in preconditioning and to eliminate the additional residual evaluations needed by a matrix-free method. The third manuscript extends the discretization method to the viscous fluxes and boundary conditions. The discretization is demonstrated to achieve the expected order of accuracy. The fourth-order scheme is also shown to be more computationally efficient than the second-order scheme at achieving grid-converged values of drag for two-dimensional laminar flow over an airfoil.  ii  Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Contents  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  List of Tables  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . .  ix  Co-Authorship Statement  . . . . . . . . . . . . . . . . . . . . . .  x  1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  1.1  1.2  Methods in Computational Aerodynamics . . . . . . . . . . .  2  1.1.1  Modeling . . . . . . . . . . . . . . . . . . . . . . . . .  2  1.1.2  Discretization  4  1.1.3  Numerical Solution  Overview of the Solver  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  6  . . . . . . . . . . . . . . . . . . . . .  7  1.2.1  Governing Equations  . . . . . . . . . . . . . . . . . .  8  1.2.2  Vertex Centered Finite-Volume Discretization  . . . .  9  1.2.3  Reconstruction . . . . . . . . . . . . . . . . . . . . . .  10  1.2.4  Interior and Boundary Flux Evaluation . . . . . . . .  14  1.2.5  Flux Integration . . . . . . . . . . . . . . . . . . . . .  14  1.3  Objectives and Outline  . . . . . . . . . . . . . . . . . . . . .  16  1.4  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . .  19  iii  Contents 2 Monotonicity Enforcement Using a Limiter  . . . . . . . . .  24  . . . . . . . . . . . . . . . . . . . . . . . . . . .  24  2.1  Introduction  2.2  High-Order Accurate Solution Reconstruction  2.3  Second-Order Limiting  2.4  2.5  . . . . . . . .  25  . . . . . . . . . . . . . . . . . . . . .  27  High-Order Limiting . . . . . . . . . . . . . . . . . . . . . . .  30  2.4.1  Monotonicity . . . . . . . . . . . . . . . . . . . . . . .  30  2.4.2  Accuracy . . . . . . . . . . . . . . . . . . . . . . . . .  32  2.4.3  Uniform Regions and Smooth Extrema  . . . . . . . .  34  2.4.4  Boundary Treatment  . . . . . . . . . . . . . . . . . .  35  2.4.5  Complete Algorithm . . . . . . . . . . . . . . . . . . .  39  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  40  2.5.1  Ringleb’s Flow . . . . . . . . . . . . . . . . . . . . . .  41  2.5.2  Transonic Flow Over an Airfoil  . . . . . . . . . . . .  48  Results  2.6  Conclusion  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  56  2.7  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . .  59  3 Efficient Convergence to Steady State . . . . . . . . . . . . .  62  3.1  Introduction  3.2  High-Order Accurate Finite-Volume Method  3.3  3.4  3.5  . . . . . . . . . . . . . . . . . . . . . . . . . . .  62  . . . . . . . . .  64  3.2.1  Governing Equations  . . . . . . . . . . . . . . . . . .  64  3.2.2  Finite-Volume Formulation . . . . . . . . . . . . . . .  64  3.2.3  Solution Reconstruction . . . . . . . . . . . . . . . . .  65  3.2.4  Numerical Flux Calculation and Integration  . . . . .  68  Globalized Newton Method . . . . . . . . . . . . . . . . . . .  69  3.3.1  Pseudo-Timestepping . . . . . . . . . . . . . . . . . .  69  3.3.2  Line Search . . . . . . . . . . . . . . . . . . . . . . . .  70  3.3.3  Selecting the Timestep  . . . . . . . . . . . . . . . . .  71  Preconditioned GMRES . . . . . . . . . . . . . . . . . . . . .  74  3.4.1  Matrix-Free Method . . . . . . . . . . . . . . . . . . .  75  3.4.2  Forming the High-Order Jacobian Matrix . . . . . . .  75  3.4.3  Preconditioning  . . . . . . . . . . . . . . . . . . . . .  78  3.4.4  Startup . . . . . . . . . . . . . . . . . . . . . . . . . .  79  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  80  Results  iv  Contents 3.5.1  Test Cases  . . . . . . . . . . . . . . . . . . . . . . . .  3.5.2  Globalization Strategy  3.5.3  Solving the Linear System  . . . . . . . . . . . . . . . . .  80 81  . . . . . . . . . . . . . . .  82  3.6  Conclusion  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  95  3.7  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . .  98  4 Extension to Viscous Flows  . . . . . . . . . . . . . . . . . . . 101  4.1  Introduction  . . . . . . . . . . . . . . . . . . . . . . . . . . . 101  4.2  Governing Equations  4.3  Solution Reconstruction . . . . . . . . . . . . . . . . . . . . . 103  . . . . . . . . . . . . . . . . . . . . . . 102  4.3.1  Reconstruction for Interior Control Volumes  . . . . . 103  4.3.2  Boundary Condition Enforcement . . . . . . . . . . . 107  4.4  Flux Integration . . . . . . . . . . . . . . . . . . . . . . . . . 111  4.5  Solution Convergence  4.6  Results  . . . . . . . . . . . . . . . . . . . . . . 112  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113  4.6.1  Cylindrical Couette Flow . . . . . . . . . . . . . . . . 113  4.6.2  Flow Over an Airfoil  . . . . . . . . . . . . . . . . . . 116  4.7  Conclusion  4.8  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 124  5 Conclusion 5.1  5.2  5.3  . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127  Significance of Results . . . . . . . . . . . . . . . . . . . . . . 127 5.1.1  Limiter . . . . . . . . . . . . . . . . . . . . . . . . . . 128  5.1.2  Globalized Matrix-Explicit Newton-GMRES  5.1.3  Viscous Discretization . . . . . . . . . . . . . . . . . . 131  . . . . . 130  Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.2.1  Turbulent Flows . . . . . . . . . . . . . . . . . . . . . 132  5.2.2  Three Dimensional Flows . . . . . . . . . . . . . . . . 135  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . 137  v  List of Tables 2.1  Convergence order of error in pressure for Ringleb’s flow . . .  48  2.2  Computational time for transonic airfoil test case . . . . . . .  54  3.1  Relative computational time of flux and Jacobian evaluation without limiter for subsonic case . . . . . . . . . . . . . . . .  3.2  83  Relative computational time of flux and Jacobian evaluation with limiter for transonic case . . . . . . . . . . . . . . . . . .  83  3.3  Average number of inner GMRES iterations per Newton iteration 87  3.4  Memory use and CPU time for subsonic case . . . . . . . . .  93  3.5  Memory use and CPU time for transonic case . . . . . . . . .  94  4.1  Order of error norms in velocity for cylindrical Couette flow . 114  4.2  Computational time for airfoil case . . . . . . . . . . . . . . . 120  4.3  Comparison of fine grid results for airfoil case . . . . . . . . . 120  vi  List of Figures 1.1  Median dual control volume . . . . . . . . . . . . . . . . . . .  10  1.2  Example reconstruction stencil . . . . . . . . . . . . . . . . .  12  1.3  Gauss integration for a control volume at a curved boundary  15  2.1  Venkatakrishnan’s smooth approximation to min(1, y) . . . .  29  2.2  min(1, y) compared to min(1, y) . . . . . . . . . . . . . . . . .  33  2.3  Transition function used to disable the limiter in nearly uniform regions . . . . . . . . . . . . . . . . . . . . . . . . . . . .  35  2.4  Two coarsest meshes used for Ringleb’s flow test case  . . . .  41  2.5  Mach number contours from Ringleb’s flow . . . . . . . . . .  43  2.6  Limiter value for pressure for Ringleb’s flow . . . . . . . . . .  45  2.7  Entropy for Ringleb’s flow . . . . . . . . . . . . . . . . . . . .  46  2.8  Error convergence for pressure in Ringleb’s flow . . . . . . . .  47  2.9  Mesh for NACA 0012 airfoil case . . . . . . . . . . . . . . . .  49  2.10 Entropy for transonic flow over NACA 0012 airfoil . . . . . .  51  2.11 Total pressure for transonic flow over NACA 0012 airfoil . . .  52  2.12 Surface pressure profiles for transonic flow over a NACA0012 airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  53  2.13 Convergence history for transonic airfoil case . . . . . . . . .  54  2.14 Surface pressure sensitivity to limiter tuning parameter . . .  55  2.15 Total pressure sensitivity to limiter tuning parameter . . . . .  57  3.1  Mesh for subsonic multi-element airfoil case . . . . . . . . . .  83  3.2  Pressure contours for subsonic flow over a multi-element airfoil 84  3.3  Mesh for transonic NACA 0012 airfoil case  . . . . . . . . . .  85  vii  List of Figures 3.4  Surface pressure coefficient for transonic flow over a NACA 0012 airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . .  87  3.5  Convergence plots for subsonic case using line search method  88  3.6  Convergence plots for transonic case using line search method  89  3.7  Convergence plots for subsonic case using different timestepping methods . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.8  Convergence plots for transonic case using different timestepping methods . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.9  90 91  Convergence for matrix-explicit method compared to matrixfree method . . . . . . . . . . . . . . . . . . . . . . . . . . . .  96  4.1  Four coarsest grids used in cylindrical Couette flow case . . . 115  4.2  Convergence of error in velocity for cylindrical Couette flow . 116  4.3  Effect of change in boundary reconstruction weight . . . . . . 117  4.4  Three coarsest grids used for airfoil case . . . . . . . . . . . . 118  4.5  Mach number contours for airfoil case . . . . . . . . . . . . . 119  4.6  Computational time for airfoil case . . . . . . . . . . . . . . . 119  4.7  Viscous drag coefficient grid convergence . . . . . . . . . . . . 121  4.8  Pressure drag coefficient grid convergence . . . . . . . . . . . 122  4.9  Grid convergence of the separation point on the upper airfoil surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122  5.1  Effect of curved boundary on anisotropic control volume . . . 135  viii  Acknowledgements First and foremost I would like to thank all of the faculty that have guided me in these years as a graduate student. My supervisor, Dr Carl Ollivier-Gooch deserves a special thank you. Carl is a wonderful person, who also happens to be a great research adviser. I would also like to acknowledge all of my student colleagues which I have had a pleasure of working with during this journey. Amongst these I was especially lucky to have a friend in the lab who was along for the whole ride. Serge, thanks for not letting me forget French! I am also grateful to all my friends with whom I’ve shared my life with outside of work, usually in the mountains. The perspective on life you provide helps me make sense of it all. A very special thank you is due to Krystil with whom I have shared my life for the last few years. Krystil, you complete me. These acknowledgements could not possibly be complete without mentioning my parents. Each in your own way, you have made me who I am, and have always supported me. Thank you.  ix  Co-Authorship Statement The research ideas and methods explored in the three co-authored manuscripts of this thesis are the fruits of a close working relationship between Dr Carl Ollivier-Gooch and Christopher Michalak. The implementation of the methods, the data analysis, and the manuscript preparation were done by Michalak with invaluable guidance from Ollivier-Gooch throughout the process.  x  Chapter 1  Introduction For almost all cases of practical interest, the equations governing fluid flow cannot be solved analytically. Therefore engineers and scientists needing to predict flow properties must depend on physical experiments or on the numerical solution of the flow equations. The latter method, commonly known as computational fluid dynamics (CFD) has grown to be an important and broad field of study. Nonetheless, although immensely useful, these methods have yet to reach a level of accuracy sufficient to completely replace the need for costly and time consuming experimental work. Although many aspects of the work presented in this thesis could be used for other applications, the immediate focus is on aerodynamics. To date, the most common goal of computational aerodynamics is the prediction of net aerodynamic forces of lift and drag to help guide the design process. Early successes in performing two-dimensional viscous simulations over airfoils and three-dimensional inviscid simulations over aircraft were documented in a review paper by Jameson in the mid 1980s [19]. The results available at the time were sufficiently accurate to gain a qualitative understanding of the flow but highly inaccurate for the quantitative prediction of drag and lift. Since then, advances in modeling, numerical methods, and computational power have greatly improved the fidelity of the results. Nonetheless the accuracy of simulations involving three-dimensional turbulent flow around aircraft configurations is still inadequate for many practical engineering purposes. In recent years a series of drag prediction workshops has been held to evaluate the state of the art in this domain. Studies of the combined results [18] of these workshops have demonstrated that both the numerical methods and the turbulence models used are responsible for these inaccuracies. This work represents a modest contribution towards improving the accuracy of these 1  1.1. Methods in Computational Aerodynamics types of simulations. However the outlook for computational aerodynamics, and therefore the need for efficient numerical methods, does not end with the accurate prediction of lift and drag. As methods improve and computational power increases, CFD applications will expand to the automatic optimization of aerodynamic configuration, the study of fluid structure interaction, and the study of flight dynamics.  1.1  Methods in Computational Aerodynamics  There are three essential components to the numerical simulation of fluid flow: 1. Modeling. The partial differential equations describing the physics of the problem must be defined. 2. Discretization. These partial differential equations must be discretized thereby transforming them into a large system of algebraic equations. 3. Numerical Solution. The value of the discrete flow properties satisfying the system of algebraic equations must be found. All three of these aspects remain active fields of research. This thesis introduces new contributions in the areas of discretization and numerical solution which together can be called the numerical method. In this section all three aspects will be introduced with an emphasis on methods most relevant to the current work.  1.1.1  Modeling  For engineering purposes, the Navier-Stokes equations accurately describe the physics of air flow. However the turbulent nature of aerodynamic flows results in large disparities between time and length scales of small eddies and macroscopic flow features. Since, in aggregate, these turbulent eddies have an important effect on the flow they cannot simply be ignored. Currently available computational power is grossly inadequate to solve the discretized  2  1.1. Methods in Computational Aerodynamics Navier-Stokes equations sufficiently accurately to capture the effects of turbulence for a typical aircraft configuration. Using an assumption that computational power increases by a factor of 5 every five years, Spalart estimated [35] that the direct simulation of turbulence for aircraft will not be possible until the year 2080. For these reasons the Reynolds averaged Navier-Stokes (RANS) equations are preferred. To close this system of equations a turbulence model is needed to estimate the effects of the fluctuation in the velocity field. The development of these models, which themselves are typically systems of partial differential equation, is an important field of research. The overall accuracy of the solution will depend on the fidelity of the turbulence model and on how accurately the RANS and turbulence model equations are solved by the numerical method. Large eddy simulation (LES) is a hybrid method in which effects of smaller eddies are simulated with empirical models while larger eddies are simulated directly. This approach can yield much more accurate results than the RANS method, but requires much greater accuracy. In 2000, Spalart estimated [35] that the grid size needed to perform LES simulations using existing numerical methods would be a factor of 30, 000 larger than that needed to perform RANS simulations. The cost of LES simulations is further exacerbated by the lack of a steady-state solution necessitating averaging of a time-accurate evolution of the solution. For these reasons Spalart estimates that the computational resources to achieve this using existing numerical approaches will not be available until 2045. Since the benefit of high-order accurate methods increases with the desired solution resolution, the prospects for these methods is even more promising for LES simulations than it is for RANS simulations. Although the numerical methods introduced in this work are ultimately intended to be applied in conjunction with a turbulence model, they are only demonstrated for inviscid and laminar flows. However some aspects of the extension of our methods to the turbulent flows will be addressed in the concluding chapter.  3  1.1. Methods in Computational Aerodynamics  1.1.2  Discretization  The purpose of discretization is to define a set of algebraic equations that approximate the partial differential equations of the model. A number of broad approaches exist for this. The most commonly used approach in computational aerodynamics, and the one used in this work, is the finitevolume method. In this method the spatial domain considered is tessellated into non-overlapping control volumes. Algebraic equations are then used to relate the control volume average states of the flow properties. 1.1.2.1  Choice of grids  Structured grids are attractive because they simplify many aspects of the solver. However generating structured grids around complex geometries is challenging. Even when a structured grid can be fitted, it often leads to grid spacing which is poorly adapted to the expected local flow gradients. Some methods, such as multi-block grids and overlapping grids, can be used to help alleviate some of these limitations. However structured grids remain difficult to generate for non-trivial geometries. For these reasons unstructured grids have gained popularity in computational aerodynamics and are used in the present work. They can be fitted to geometries of any complexity and the local grid density can readily be controlled. Unlike with structured grids, the process of generating the mesh can be almost completely automated. However unstructured grids also lead to greater complexity of the solver since the mesh connectivity must be stored and since the discretization cannot take advantage of a structured topology. The first solutions of the Euler equations on unstructured grids were carried out using the finite-element method in the 1980s [1, 22]. Shortly thereafter work using the finite-volume method was carried out by Jameson and Mavriplis [21]. 1.1.2.2  Order of accuracy  For any consistent numerical method, the error resulting from the discretization of the partial differential equations should decrease when the grid is 4  1.1. Methods in Computational Aerodynamics refined. Specifically, a discretization will produce an error of O (hp ) where h is the mesh length scale and p is the order of accuracy. The disadvantage of methods achieving a higher order of accuracy is that they are more complicated to implement and are more computationally expensive on a given grid. However to achieve the same level of accuracy a high-order accurate method can use a coarser grid than a lower-order accurate method. Therefore, from a computational cost perspective, the choice of discretization accuracy is not an obvious one. Due to their relative simplicity, second-order accurate schemes are by far the most commonly used methods for computational aerodynamics. The question of whether high-order accurate schemes could achieve desired levels of accuracy more efficiently has not been adequately addressed. A third-order accurate method for structured grids has been studied by De Rango and Zingg [12]. In this study two-dimensional turbulent flow over airfoils was considered. The high-order accurate method was shown to obtain adequately accurate solutions on coarser grids than the second-order scheme. Furthermore the overall computational cost of attaining this level of accuracy was shown to be lower for the third-order scheme than the second-order scheme due to the coarser mesh requirement. On unstructured grids a nominally third-order accurate finite-volume methods on unstructured grids for the Euler equations was introduced by Barth and Frederickson [3–5]. This work introduced k-exact reconstruction, which was further studied by other researchers [13, 16, 38]. These latter works used a third-order k-exact reconstruction for the advective fluxes and Green-Gauss gradient to compute viscous fluxes. Although a qualitative assessment of accuracy for practical cases was made in these works, the effective nominal accuracy of the scheme has not been demonstrated. Furthermore the superiority of the method, in terms of computational effort to attain a level of accuracy relative to second-order schemes, has also not been demonstrated. Considerable work using the Euler equation has also previously been carried out by others in the present author’s group [26–28]. These works successfully demonstrated that the nominal level of accuracy for schemes up to fourth-order can be obtained for shockless flows. Work in the 5  1.1. Methods in Computational Aerodynamics group has also been previously carried out demonstrating the effectiveness of the k-exact reconstruction in obtaining high-order accurate solutions of the advection-diffusion equation [29]. The discontinuous Galerkin finite-element method has also been successfully applied to the high-order accurate solution of compressible flow [7, 8, 15]. As with the high-order finite-volume method, the discontinuous Galerkin method has not yet matured enough to be used in production codes. The computational efficiency, in particular, has not been shown to exceed that of second-order finite-volume methods in common use. A disadvantage of this method is that, unlike in the finite-volume method, the number of degrees of freedom increases with the order of accuracy on a given mesh.  1.1.3  Numerical Solution  Once a spatial discretization has been selected, the remaining step is to address the time evolution of the system. For aerodynamics problems modelled using RANS, often only the steady-state solution is of interest. In this case, the time dependent terms could, in theory, be eliminated from the system of equations. If this is done, the remaining problem is to find the root of the equations defining the space discretization. In practice this root finding problem is most easily solved by reintroducing a time-like term to the equations. Therefore, whether a time-dependent or steady-state solution is sought, some form of time integration is needed. Explicit time integration methods, such as the Runge-Kutta schemes, are relatively simple to implement since they only require a means of evaluating the residual of the space discretized equations. However their stability is limited by the Courant–Friedrichs–Lewy (CFL) condition. This limits the size of the timesteps that can be taken relative to the grid spacing. For typical grids used in computational aerodynamics this results in the need for an inordinate number of residual evaluations. Therefore the computational effort needed by these schemes is prohibitively high. For steady-state problems, the efficiency of explicit time-advance methods can be greatly improved through the use of multigrid methods [11, 20, 23, 24, 30].  6  1.2. Overview of the Solver Alternatively, implicit time integration methods can be used. Unlike explicit methods, the timestep is not restricted by the CFL condition. Implicit methods can therefore be substantially more efficient, particularly when only a steady-state solution is needed. At the limit when the timestep is infinitely large, the approach reduces to the Newton method. Typically the solution strategy is to use an implicit timestep until the solution is sufficiently converged to transition to a full-Newton method. Alternatively, the timestep can be ramped up gradually and the method approaches the Newton method asymptotically. However implicit methods are considerably more complex to implement since the solution of a large linear system involving the Jacobian of the residual is needed. Early approaches for solving the linear systems arising from the implicit method were based on direct solution techniques [2, 36, 37]. However the memory and processing requirements of direct methods quickly become intractable. Therefore attention quickly turned to solving these linear systems approximately using iterative methods. In particular Generalized Minimal Residual (GMRES) [33], part of the Krylov subspace family of methods, has been shown to be very effective. On structured grids Pueyo and Zingg [31] demonstrated an efficient Newton-GMRES method for aerodynamic flow computations. A number of researchers have successfully applied NewtonGMRES on unstructured grids [6, 9, 26, 39].  1.2  Overview of the Solver  This section provides a brief overview of the high-order accurate finitevolume solver used in this work. Focus will be placed on the components of the solver that are not thoroughly addressed in the manuscript chapters. Specifically, the focus here will be on the spatial discretization of the inviscid flow equations. The extension to the viscous case is addressed in Chapter 4 and strategies for solving the discrete equations for steady state problems are addressed in Chapter 3.  7  1.2. Overview of the Solver  1.2.1  Governing Equations  The two-dimensional Euler equations for a control volume Vi can be written in integral form as d dt   ˆ  ρ  ˛ U dA +  F dS = 0  Vi  Ωi      (1.1)   ρun    ρuu + P n n x  F =  ρvun + P ny      ρu    U =   ρv     (E + P )un  E         (1.2)  where U is the solution vector in conservative variables, F is the flux vector normal to the control volume boundary, n ˆ x and n ˆ y are the unit normal components, and un = uˆ nx + vˆ ny . The conserved variables consist of the density (ρ), two components of momentum (ρu and ρv), and energy (E). Air is modeled as an ideal gas with the equation of state P = ρRT where P is pressure, R is the specific gas constant for air, and T is the temperature. The following thermodynamic relationships are also used γ=  Cp , Cv  e = Cv T,  Cp =  γR , γ−1  Cv =  1 et = e + (u2 + v 2 ), 2  R , γ−1  R = Cp − Cv  E = ρet ,  h = et +  P ρ  where γ is the ratio of specific heats for air, Cp and Cv are the specific heats, e is the specific internal energy, et is the specific total energy, E is the total energy, and ht is the specific total enthalpy. An important property for compressible flows is the Mach number M which is the ratio of the flow speed to the speed of sound a √ M  =  u2 + v 2 a 8  1.2. Overview of the Solver a =  γRT  To reduce round-off error and to improve the conditioning of the discretization, the flow properties are normalized. For this, the following reference quantities are used ρref = ρ∞ ,  Tref = T∞ ,  uref = vref = a∞ ,  Pref = γP∞  Rref = R∞ = const  where [·]ref is the reference value, and [·]∞ is the far-field value for external aerodynamics problems or some other representative value for internal flows. The governing equations can then be expressed in terms of the non-dimensionalized quantities ρn =  ρ  ,  Tn =  ρref u v un = , vn = , uref vref  T , Tref  Pn =  Rn =  P Pref  R =1 Rref  The important non-dimensionalized constitutive relations become Pn =  ρn Tn , γ  an =  Tn ,  En =  1 Pn + ρn u2n + vn2 , γ−1 2  hn =  En + Pn ρn  For the remainder of this thesis we will refer only to the normalized values and the subscript n is therefore dropped.  1.2.2  Vertex Centered Finite-Volume Discretization  The grids used in the present work are generated using a Delaunay refinement algorithm [10, 34]. In two dimensions these algorithms generate meshes composed of triangular cells. The solver used in the present work uses the median dual control volume formulation shown in Figure 1.1. Control volumes are associated with each vertex and are defined by connecting cell centroids with face centroids of the primal mesh. The dependent variables, which will be solved for, represent the control volume averaged values of the  9  1.2. Overview of the Solver  Primal Mesh Median Dual  Figure 1.1: Median dual control volume flow properties Ui =  1.2.3  1 Ai  ˆ U dA Vi  Reconstruction  Reconstruction is a central component to the high-order accurate finitevolume method. Some of the details of reconstruction are therefore repeated to various extent in the manuscript chapters of this thesis. However, for completeness, the fundamental details of reconstruction for the inviscid scheme are summarized here. The extension to the viscous case is described in Chapter 4. To obtain a high-order accurate numerical approximation to the residual we reconstruct a polynomial representation of the solution within each control volume. As is common, the reconstruction is carried out in the primitive variables U = [ρ, u, v, P ]T rather than the conserved variables. As in previous work [29], the high-order accurate reconstruction is obtained by formulating and solving a least-squares problem. To obtain an accurate estimate of the solution at flux quadrature points, we represent the solution within each control volume by the Taylor series  10  1.2. Overview of the Solver expansion UiR (x − xi , y − yi ) =  U |i +  ∂U ∂U (x − xi ) + (y − yi ) + ∂x i ∂y i  ∂2U ∂x2  (x − xi )2 ∂2U (x − xi ) (y − yi ) + + 2 ∂x ∂y i  i  ∂2U ∂y 2  i  (y − yi )2 + ··· 2  where UiR is the value of the reconstructed solution and  (1.3) ∂ k+l Ui ∂xk ∂y l  are its  derivatives at the reference point (xi , yi ) of control volume i. The reference point typically used for vertex-centered schemes is the vertex location. The order of accuracy that the numerical method will achieve is, amongst other factors, limited by the truncation error of this reconstruction. Truncating this polynomial to degree k will result in a scheme which is at most (k + 1)-order accurate. In the present work we consider reconstructions of up to cubic degree resulting in schemes up to fourth-order accurate. The coefficients of the polynomial are computed so that the mean value of the solution in the control volume is conserved and the reconstruction approximates nearby control volume averages. The neighboring control volumes used for reconstruction are selected based on their topological distance as shown in Figure 1.2. Entire topological layers are added until the number of neighbors equals or exceeds a set minimum. In the present work we use a minimum of 4 neighbors for linear reconstruction, 10 neighbors for quadratic reconstruction, and 16 neighbors for cubic reconstruction. Conservation of the mean within a control volume is satisfied when ˆ 1 ∂U ∂U Ui = UiR dA ≡ U |i + xi + y + Ai Vi ∂x i ∂y i i ∂2U ∂x2 ∂2U ∂y 2  i  i  x2 i ∂2U + xy + 2 ∂x ∂y i i y2i + ··· 2  (1.4)  11  1.2. Overview of the Solver 3  3  3 2  2  3  3  3 R  2  3  2 3 2  2 3  3  3  3  3  Figure 1.2: Example reconstruction stencil for second- and third-order scheme. In this example the fourth-order stencil would be identical to the third-order stencil. where xn y m  i  1 ≡ Ai  ˆ (x − xi )n (y − yi )m dA Vi  As we shall see, this mean constraint can be eliminated analytically from the least-squares system. The terms used to construct the reconstruction least-squares problem have the form 1 Aj  ˆ UiR (x − xi ) dA = Vj  U |i + +  ∂2U ∂x2  ∂U ∂U xij + yij ∂x i ∂y i  i  x2 ij ∂2U ∂2U xy ij + + 2 ∂x ∂y i ∂y 2  (1.5)  i  y 2 ij + ··· 2  The geometric terms in this equation are of the general form xn y m ij  ≡  1 Aj m  =  ˆ ((x − xj ) + (xj − xi ))n · ((y − yj ) + (yj − yi ))m dA Vj n  n! m! (xj − xi )k · (yj − yi )l · xn−k y m−l j l! (m − l)! k! (n − k)! l=0 k=0  The resulting least-squares problem to be solved for each solution variable in 12  1.2. Overview of the Solver each control volume takes the form    xi  yi  x2 i  xy i  y2i  ···  wi1 xi1  wi1 yi1  wi1 x2 i1  wi1 xy i1  wi1 y 2 i1  ···  1    w  i1   wi2    wi3   .  ..   wiN    wi2 xi2  wi2 yi2  wi3 xi3 .. .  wi3 yi3 .. .  wiN xiN  wiN yiN  Ui  wi2  x2  wi3  x2  i2  wi2 xy i2  .. .  wi3 xy i3 .. .  wiN x2 iN  wiN xy iN  i3  wi2  y2  i2  wi3  y2  i3  ···  .. .  ··· .. .  wiN y 2 iN  ···                   U    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x ∂y 1 ∂2U 2 ∂y 2                  .. .  i       w U  i1 1      wi2 U 2    =   wi3 U 3      ..   .    (1.6)  wiN U N  where the weights can be used to emphasize geometrically nearby data wij =  1 |xj − xi |n  (1.7)  where typically n ∈ [0, 2]. For isotropic or lightly anisotropic grids (such as those used for laminar flows) the choice of n does not impact the accuracy of the scheme but does affect its conditioning and stability. In the present work we use n = 1 for the inviscid cases and n = 0 for viscous cases. The first equation in the linear system is the mean constraint, which can be removed by Gauss elimination, leaving an unconstrained least-squares problem of reduced size. The efficient solution of these least-squares systems is discussed in Chapter 3. Once the polynomial coefficients are known, the reconstructed flow properties at any point in the control volume can easily be found using Equation 1.3.  13  1.2. Overview of the Solver  1.2.4  Interior and Boundary Flux Evaluation  To maintain the numerical stability of the system, careful upwinding of the Euler system of equations is needed. The discretization of the flux follows the pioneering work of Godunov [17] who suggested solving what amounts to a shock tube–or Riemann–problem at the interface of control volumes. Since solving this Riemann problem exactly is computationally costly, a number of approximate Riemann solvers have been developed. In the present work, Roe’s approximate Riemann solver [32] is used. The reconstructed solution from each of two control volumes incident on the Gauss point is used to define the two states of the Riemann solver. The inviscid boundary conditions are enforced by using a special flux function at boundary Gauss points. In general this flux is simply taken as the analytic flux of Equation 1.2 computed using some information from the single reconstructed value from the interior and some information from the boundary condition. Specifically, for walls, we use the reconstructed value of ρ and P and a zero velocity normal to the wall. For inflows and outflows the boundary flux is constructed based on an eigenvalue analysis. Specifically, for supersonic inflows all values are taken from the boundary condition. For subsonic inflows the reconstructed value of pressure from the interior is used along with the total pressure, total temperature, and flow direction from the upstream condition. Analogously, for supersonic outflows the flux is based entirely on reconstructed values from the interior while for subsonic outflows pressure is defined by the boundary condition while the remaining flow properties are obtained from the reconstruction.  1.2.5  Flux Integration  As part of the numerical method, the contour integral in Equation 1.1 needs to be discretized. This is obtained through the use of Gauss quadrature. To achieve the desired order of accuracy, a single Gauss point per median dual edge is used for first- and second-order schemes and two Gauss points per edge are used for third- and fourth-order schemes. This represents, respectively, two and four points per edge of the primal mesh. 14  1.2. Overview of the Solver  Primal Mesh Median Dual Control−Volume Curved Boundary Gauss Point Normal  Figure 1.3: Approximate location and normal of Gauss integration points for a third- or fourth-order scheme on a boundary control volume To achieve third- and high-order accurate schemes it is well known that Gauss points must be correctly located on curved boundaries [27, 29]. Recently it has also been shown that for vertex-centered schemes the curvature of the boundary must be accounted for when determining the Gauss point locations and the normals to the interior median dual edge incident on the boundary [14]. For this reason, the Gauss point locations, normals and weights are all computed from the curved boundary information as shown in Figure 1.3. In the present work we are concerned with steady-state flows. In these problems, the time derivative term in Equation 1.1 is eliminated and the remaining task is to find the root of the system of equations defining the discrete flux integrals. The novel approach taken in this work, and how it compares to more commonly used methods, is discussed in Chapter 3.  15  1.3. Objectives and Outline  1.3  Objectives and Outline  The objective of this thesis is to demonstrate that carefully implemented high-order accurate unstructured finite-volume methods for computational aerodynamics can be more computationally efficient than the second-order methods commonly used. This broad objective is attained by considering a number of specific, previously unresolved issues. Although ultimately this study is aimed at resolving issues with the simulation of three-dimensional turbulent flows, the current study is carried out in two-dimensions. Furthermore, only inviscid and laminar viscous flows are considered. We use these simplifications since certain aspects of implementing the high-order method for turbulent three-dimensional flows have not yet been addressed. These challenges will be discussed in the concluding chapter. The core chapters of this thesis are composed of three manuscripts submitted for journal publication. Further motivation and relevant background information is therefore addressed in each individual chapter. Many aerodynamic flows of engineering interest occur in the transonic or supersonic regimes. Chapter 2 addresses the previously unresolved issue of how to eliminate unwanted oscillations near solution discontinuities such as shocks while maintaining the accuracy of the high-order method in smooth regions. The chapter studies the extension of limiters used for second-order methods to the high-order case. The requirements for accuracy and efficient convergence are discussed. A new limiting procedure is proposed. Ringleb’s flow problem is used to demonstrate that nearly nominal orders of accuracy for schemes up to fourth-order can be achieved in smooth regions using the new limiter. The results for the fourth-order accurate solution of transonic flow demonstrates good convergence properties and significant qualitative improvement of the solution relative to the second-order method. It is also demonstrated that the new limiter can reduce the dissipation of second-order schemes with minimal sacrifices in convergence properties relative to existing approaches. Although the matrix-free Newton-GMRES method commonly used to converge the solution of the second-order scheme to steady-state has been 16  1.3. Objectives and Outline shown to be relatively effective for the third-order scheme, it has also been shown to be less effective for the fourth-order scheme [25, 26]. Chapter 3 successfully addresses the poor convergence to steady sate of the fourth-order scheme through the use of novel Newton globalization techniques and the use of the exact flux Jacobian with a matrix-explicit GMRES method. We show that robustness of the globalization can be improved by supplementing the pseudo-timestepping method commonly used with a line search method. The number of timesteps required for convergence can be reduced by using a timestep that scales with the local residual. We also show that it is possible to form the high-order Jacobian explicitly at a reasonable computational cost. This is demonstrated for cases using both limited and unlimited reconstruction. This exact Jacobian can be used for preconditioning and directly in the GMRES method. The benefits resulting from improvements in preconditioning and the elimination of residual evaluations in the inner iterations of the matrix-free GMRES method are substantial. Computational results focus on second- and fourth-order accurate schemes with some results for the third-order scheme. Overall computational cost for the matrix-explicit method is lower than the matrix-free method for all cases. The fourth-order matrix-explicit scheme is a factor of 1.6 to 3 faster than the matrix-free scheme while requiring about 50-100% more memory. In Chapter 4, the methods successfully applied for the high-order solution of inviscid flow are extended to laminar viscous flow. A means of enforcing the viscous boundary conditions by introducing additional rows in the reconstruction matrices is devised. The solution gradients needed for the viscous fluxes are obtained from the same reconstruction polynomials used for the inviscid flux. Rapid convergence of steady-state problems is made possible by computing the exact flux Jacobian and the use of a Newton-GMRES algorithm. A grid convergence study of cylindrical Couette flow is used to demonstrate the order of accuracy of the scheme. Laminar flow over an airfoil on a sequence of anisotropic grids is also used to qualitatively demonstrate the superior accuracy of the fourth-order scheme. The conclusion in Chapter 5 discusses the combined significance of the methods and results presented in the manuscript chapters. A number of ex17  1.3. Objectives and Outline tensions to the present work are suggested. In particular the extensions of the methods to three dimensions and turbulent flows are discussed. The expected difference in the computational efficiency trade-offs for three-dimensional methods relative to their two-dimensional counterparts are addressed. Specific issues of high-order methods applied to the highly anisotropic grids needed to solve high Reynolds number flows are discussed. The need for a stable high-order accurate discretization of turbulence model is also addressed.  18  1.4. Bibliography  1.4  Bibliography  [1] F. Angrand and A. Dervieux. Some explicit triangular finite element schemes for the Euler equations. International Journal for Numerical Methods in Fluids, 4:749–764, Aug. 1984. [2] H. E. Bailey and R. M. Beam. Newton’s method applied to finite difference approximations for steady state compressible Navier-Stokes equations. Journal of Computational Physics, 93:108–127, 1991. [3] T. J. Barth. Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations. In Unstructured Grid Methods for Advection-Dominated Flows, pages 6–1 – 6–61. AGARD, Neuilly sur Seine, France, 1992. AGARD-R-787. [4] T. J. Barth. Recent developments in high order k-exact reconstruction on unstructured meshes. AIAA paper 93-0668, Jan. 1993. [5] T. J. Barth and P. O. Frederickson. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA paper 90-0013, Jan. 1990. [6] T. J. Barth and S. W. Linton. An unstructured mesh Newton solver for compressible fluid flow and its parallel implementation. AIAA paper 95-0221, Jan. 1995. [7] F. Bassi, A. Crivellini, S. Rebay, and M. Savini. Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ω turbulence model equations. Computers and Fluids, 34(4-5):507 – 540, 2005. Residual Distribution Schemes, Discontinuous Galerkin Schemes and Adaptation. [8] F. Bassi and S. Rebay. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. Journal of Computational Physics, 131:267–279, 1997.  19  1.4. Bibliography [9] M. Blanco and D. W. Zingg. Fast Newton-Krylov method for unstructured grids. American Institute of Aeronautics and Astronautics Journal, 36(4):607–612, 1998. [10] C. Boivin and C. F. Ollivier-Gooch. Guaranteed-quality triangular mesh generation for domains with curved boundaries. International Journal for Numerical Methods in Engineering, 55(10):1185–1213, Dec. 2002. [11] A. Brandt. Guide to multigrid development. In W. Hackbusch and U. Trottenberg, editors, Multigrid Methods, volume 960 of Lecture Notes in Mathematics. Springer-Verlag, 1982. [12] S. De Rango and D. W. Zingg. Higher-order spatial discretization for turbulent aerodynamic computations. American Institute of Aeronautics and Astronautics Journal, 39(7):1296–1304, July 2001. [13] M. Delanaye and J. A. Essers. Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids. American Institute of Aeronautics and Astronautics Journal, 35(4):631–639, Apr. 1997. [14] B. Diskin and J. Thomas. Accuracy analysis for mixed-element finitevolume discretization schemes. Technical Report 2007-08, National Institute of Aerospace, 2007. [15] K. J. Fidkowski, T. A. Oliver, J. Lu, and D. L. Darmofal. p-multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations. Journal of Computational Physics, 207(1):92 – 113, 2005. [16] P. Geuzaine, M. Delanaye, and J.-A. Essers. Computation of high Reynolds number flows with an implicit quadratic reconstruction scheme on unstructured grids. In Proceedings of the Thirteenth AIAA Computational Fluid Dynamics Conference, pages 610–619. American Institute of Aeronautics and Astronautics, 1997.  20  1.4. Bibliography [17] S. K. Godunov. A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics. Matematicheskie Sborniki, 47:357–393, 1959. [18] M. J. Hemsch. Statistical analysis of CFD solutions from the Drag Prediction Workshop. In Fortieth AIAA Aerospace Sciences Meeting, Jan. 2002. AIAA paper 2002-842. [19] A. Jameson. Successes and challenges in computational aerodynamics. In Proceedings of the Eighth AIAA Computational Fluid Dynamics Conference, June 1987. AIAA paper 87-1184. [20] A. Jameson. Artificial diffusion, upwind biasing, limiters and their effect on accuracy and multigrid convergence in transonic and hypersonic flows. AIAA paper 93-3359-CP, July 1993. [21] A. Jameson and D. Mavriplis. Finite volume solution of the twodimensional Euler equations on a regular triangular mesh. American Institute of Aeronautics and Astronautics Journal, 24(4):611–618, 1986. [22] R. Loehner, K. Morgan, J. Peraire, and O. C. Zienkiewicz. Finite element methods for high speed flows. In Proceedings of the Eighth AIAA Computational Fluid Dynamics Conference, pages 403–410, 1985. AIAA paper 1985-1531. [23] D. Mavriplis. Accurate multigrid solution of the Euler equations on unstructured and adaptive meshes. American Institute of Aeronautics and Astronautics Journal, 28(2):213–221, Feb. 1990. [24] D. Mavriplis. Directional agglomeration multigrid techniques for high reynolds number viscous flow solvers. Thirty-sixth AIAA Aerospace Sciences Meeting, Jan. 1998. [25] A. Nejat and C. Ollivier-Gooch. Effect of discretization order on preconditioning and convergence of a high-order unstructured NewtonGMRES solver for the Euler equations. Journal of Computational Physics, 227(4):2366–2386, 2008. 21  1.4. Bibliography [26] A. Nejat and C. Ollivier-Gooch. A high-order accurate unstructured finite volume Newton-Krylov algorithm for inviscid compressible flows. Journal of Computational Physics, 227(4):2592–2609, 2008. [27] C. Ollivier-Gooch, A. Nejat, and C. Michalak. On obtaining and verifying high-order finite-volume solutions to the Euler equations on unstructured meshes. American Institute of Aeronautics and Astronautics Journal, Submitted, 2008. [28] C. F. Ollivier-Gooch. High-order ENO schemes for unstructured meshes based on least-squares reconstruction. AIAA paper 97-0540, Jan. 1997. [29] C. F. Ollivier-Gooch and M. Van Altena. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. Journal of Computational Physics, 181(2):729–752, 2002. [30] N. Pierce, M. Giles, A. Jameson, and L. Martinelli. Accelerating threedimensional Navier-Stokes calculations. In Proceedings of the Thirteenth AIAA Computational Fluid Dynamics Conference, pages 676–698. American Institute of Aeronautics and Astronautics, July 1997. AIAA 97-1953-CP. [31] A. Pueyo and D. W. Zingg. Efficient Newton-Krylov solver for aerodynamic computations. American Institute of Aeronautics and Astronautics Journal, 36(11):1991–1997, 1998. [32] P. L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43:357–372, 1981. [33] Y. Saad and M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal of Scientific and Statistical Computing, 7(3):856–869, July 1986. [34] J. R. Shewchuk. Delaunay refinement algorithms for triangular mesh generation. Computational Geometry: Theory and Applications, 22(1– 3):21–74, May 2002. 22  1.4. Bibliography [35] P. R. Spalart. Strategies for turbulence modelling and simulations. International Journal of Heat and Fluid Flow, 21(3):252 – 263, 2000. [36] V. Venkatakrishnan. Newton solution of inviscid and viscous problems. Twenty-sixth AIAA Aerospace Sciences Meeting, 1988. AIAA Paper 88-0143. [37] V. Venkatakrishnan and T. J. Barth. Application of direct solvers to unstructured meshes for the Euler and Navier-Stokes equations using upwind schemes. Twenty-seventh AIAA Aerospace Sciences Meeting, 1989. AIAA Paper 89-0364. [38] D. Vigneron, J.-M. Vaassen, and J.-A. Essers. An implicit finite volume method for the solution of 3D low Mach number viscous flows using a local preconditioning technique. Journal of Computational and Applied Mathematics, 215(2):610 – 617, 2008. [39] P. Wong and D. W. Zingg. Three-dimensional aerodynamic computations on unstructured grids using a Newton-Krylov approach. Computers and Fluids, 37:107–120, 2008.  23  Chapter 2  Monotonicity Enforcement Using a Limiter† 2.1  Introduction  High-order discretizations have been shown to reduce computational effort on structured grids [6, 22]. High-order finite-volume methods on unstructured grids, although well known [3, 4, 7], have not yet effectively been applied to large scale aerodynamics problems. An outstanding issue with these methods is how to deal with discontinuities, such as shocks in the flow, while maintaining good accuracy and convergence. One means of dealing with discontinuities is to use the classic MUSCL [20] scheme with the addition of a slope limiter. For the second-order case, Barth and Jespersen [5] demonstrated the use of limited reconstruction for the solution of the Euler equations. Efficient convergence to steady state was achieved by Venkatakrishnan [21] by modifying the limiter to be differentiable. Third-order accurate schemes using k-exact reconstruction have been demonstrated first by Barth and Frederickson [4] and subsequently by others [2, 7, 9, 13, 14, for example]. Transonic and supersonic solutions have been computed in some of these works by various extensions of second-order limiters. However, the work of Barth [2] presents a limiting approach which causes difficulties in steady-state convergence, while other works [7, 13] present approaches that do not strictly enforce monotonicity and therefore allow some undesirable oscillations to occur. Furthermore none of these works †  A version of this chapter will be submitted for publication. C. Michalak and C. Ollivier-Gooch. Accuracy preserving limiter for the high-order accurate solution of the Euler equations.  24  2.2. High-Order Accurate Solution Reconstruction formally demonstrate that high-order accuracy is maintained in smooth regions of the flow. An alternative to MUSCL for obtaining high-order accurate solutions is the essentially non-oscillatory (ENO) scheme [1, 10, 18, for instance]. These methods avoid the need for slope limiters by selecting a smooth flux stencil at each iteration. Due to the inherent non-differentiability of this process, convergence of the solution to steady state is not possible. Weighted ENO (WENO) schemes [8, 11, 12, for instance] were introduced in part to resolve this issue. However these schemes do not converge to steady state as efficiently as MUSCL schemes. The computational cost per residual evaluation is also much higher than for reconstruction based solvers. A hybrid between WENO and MUSCL schemes named Quasi-ENO [15] has similar limitations. This method’s reconstruction step is much more expensive than traditional MUSCL since the reconstruction least-squares matrix changes at each iteration. For these reasons ENO-like schemes have not seen widespread application to aerodynamics problems. The present work formulates the requirements and presents a candidate for a limiter that achieves fourth-order accurate solutions in smooth regions while maintaining monotonicity and good convergence properties. An overview of the high-order MUSCL scheme is given in Section 2.2. Second-order limiters are reviewed in Section 2.3. Our extension of these methods to high-order schemes is presented in Section 2.4. Finally, we present results in Section 2.5 that demonstrate that our method achieves nominal order of accuracy in smooth regions of the flow while effectively eliminating oscillations at shocks and maintaining good convergence properties.  2.2  High-Order Accurate Solution Reconstruction  The third- and fourth-order accurate reconstruction procedure we use here is documented by Ollivier-Gooch and Van Altena [16] and is briefly reviewed in this section. Only the equations that are needed for the discussion of limiters  25  2.2. High-Order Accurate Solution Reconstruction are presented. In the finite-volume method, the domain is tessellated into non-overlapping control volumes. Each control volume Vi has a geometric reference point xi . While in principle any point can be chosen as the reference point, the usual choices (which we recommend) are the cell centroid for cell-centered control volumes and the vertex for vertex-centered control volumes. For any smooth function U (x) and its control volume averaged values U i , the k-exact least-squares reconstruction will use a compact stencil in the neighborhood of control volume i to compute an expansion UiR (x − xi ) that conserves the mean in control volume i and reconstructs exactly polynomials of degree ≤ k (equivalently, UiR (x − xi ) − U (x) = O ∆xk+1 ) . Conservation of the mean requires that the average of the reconstructed function UiR and the original function U over control volume i be the same: Ui ≡  1 Vi  ˆ UiR (x − xi ) dA = Vi  1 Vi  ˆ U (x) dA  (2.1)  Vi  The expansion UiR (x − xi ) can be written as:  UiR (x − xi ) =  U |xi + + +  ∂2U ∂x2 ∂2U ∂y 2  ∂U ∂x  xi  xi  (x − xi ) + xi  ∂U ∂y  (x − xi )2 ∂2U + 2 ∂x∂y  (y − yi ) xi  ((x − xi )(y − yi )) (2.2) xi  (y − yi )2 + ··· 2  Taking the control volume average of this expansion over control volume i and equating it to the mean value gives Ui =  U |xi + +  ∂2U ∂x2  ∂U ∂x  xi  xi + xi  ∂U ∂y  x2 ∂2U + 2 ∂x∂y  yi  (2.3)  xi  xy + xi  ∂2U ∂y 2  xi  y2 + ··· 2  26  2.3. Second-Order Limiting where xn y m i ≡  1 Ai  ˆ (x − xi )n (y − yi )m dA.  (2.4)  Vi  are control volume moments. This condition, which must be satisfied exactly, is combined with the reconstruction goal of approximating nearby control volume averages to obtain a constrained least-squares problem for the solution of the Taylor series expansion coefficients. Since the resulting least-squares matrix depends only on geometric terms, its pseudoinverse may be found in a preprocessing step. Therefore the reconstruction step at each flux evaluation is reduced to a matrix-vector product and the exact flux Jacobian can be computed as described in the manuscript reproduced in Chapter 3. Ollivier-Gooch [15] presents a modification to the reconstruction procedure resulting in a quasi-ENO scheme. This scheme eliminates the requirement for a limiter by varying the weights of the rows in the least-squares matrix at each iterations based on a measure of smoothness. However, since the pseudoinverse can no longer be precomputed, this scheme is computationally expensive.  2.3  Second-Order Limiting  To avoid introducing oscillation in the solution process, no new local extrema must be formed during reconstruction. Barth and Jespersen [5] introduced the first limiter for unstructured grids. The scheme consists of finding a limiter value Φi for each primitive flow variable in each control volume that will limit the gradient in the piecewise-linear reconstruction of the solution. In the second-order reconstruction case, if the reference location xi is taken to be the control volume centroid, the point-wise value U |xi is equal to the control volume average U i . This leads to the limited reconstruction of the form UiR (x − xi , Φi ) = U i + Φi  Ui · (x − xi ), Φ ∈ [0, 1]  The goal is to find the largest Φi which prevents the formation of local extrema at the flux integration Gauss points. The following procedure is used by Barth and Jespersen : 27  2.3. Second-Order Limiting 1. Find the largest negative (δUimin = min(U j −U i )) and positive (δUimax = max(U j − U i )) difference between the solution in the immediate neighbors j and the current control volume i. 2. Compute the unconstrained reconstructed value at each Gauss point (Uik = UiR (xk − xi )). 3. Compute a maximum allowable value of Φik for each Gauss point k.  Φik =   δU max  min 1, U i−U ,   i ik  min      δUi ik −U i  if Uik − U i > 0  min 1, U  , if Uik − U i < 0  1,  if Uik − U i = 0  4. Select Φi = min(Φik ) . 5. Compute the limited reconstruction UiR (x − xi , Φi ) at flux Gauss integration points. Clearly, steps 1, 3, and 4 introduce non-differentiability in the computation of the reconstructed function. Consequently, the second-order flux is also nondifferentiable. This has severe adverse effect on the convergence properties of the solver. This is particularly evident for implicit schemes, but even explicit time advance schemes are unable to obtain more than two or three orders of magnitude in residual reduction. In practice, the non-differentiability of step 3 causes the greatest degradation in convergence performance. For this reason, Venkatakrishnan [21] introduces a smooth alternative to step 3 of the Barth-Jespersen procedure by replacing the function min(1, y) with φ(y) =  y 2 + 2y y2 + y + 2  (2.5)  The effect of this modification can be seen in Figure 2.1. This function is differentiable and is entirely contained within the monotonicity bounds of the original limiter. Furthermore by satisfying the condition φ(2) = 1 this limiter can be shown to preserve second-order accuracy in regions where 28  2.3. Second-Order Limiting  min(1,y) Venkatakrishnan  1.2 1 0.8 0.6 0.4 0.2 0 0  0.5  1  1.5  2 y  2.5  3  3.5  4  Figure 2.1: Venkatakrishnan’s smooth approximation to min(1, y) no extrema exist for perfectly uniform meshes. Specifically for smooth solutions on a uniform mesh, the value of the Gauss point is expected to be located at the midpoint connecting two control volume centroids. Therefore for any smooth function  δUimax Uik −U i  is expected to be 2 ± O (xi − xik ), and  φ(y) is therefore expected to be 1 ± O (xi − xik ). Hence, modifying the gradient in this manner for second-order solutions introduces an error in the reconstruction which is on the order of truncation for smooth flows on uniform grids. However, for general unstructured grids the Gauss point can be located at a distance O (xi − xik ) away from the midpoint between centroids. Therefore, as will be seen in the results, this limiter leads to accuracy loss relative to the unlimited scheme. A further modification introduced by Venkatakrishnan is a method to avoid applying the limiter in regions of nearly uniform flow and smooth extrema. In these regions we expect the solution to vary such that Uik − U i = O ∆x2 where ∆x is the characteristic length of the control volume i. Therefore if the effect of the limiter can be eliminated when Uik − U i ≤ (K∆x)1.5 , where K is a tunable parameter, accuracy in smooth extrema regions can be improved without compromising monotonicity enforcement in 29  2.4. High-Order Limiting other regions. Venkatakrishnan uses a method inspired by van Albada [19] to make the switch between limited and unlimited regions of flow using a differentiable function so that convergence properties are not adversely affected. For the case Uik − U i > 0 the limiter becomes φik =  (∆2+ + 2 )∆− + 2∆2− ∆+ 1 ∆− ∆2+ + 2∆2− + ∆− ∆+ + 2  In this equation ∆− = Uik − U i , ∆+ = δUimax and  (2.6) 2  = (K∆x)3 . In  addition to improving accuracy this modification is known to be critical to achieving good convergence. When the solution is perturbed in nearly uniform regions or near smooth extrema the results of steps 1 and 4 of the limiting procedure can be expected to change more frequently than in nonuniform regions. Therefore the non-differentiability of these steps is a much greater hindrance to convergence in uniform regions than in non-uniform regions. By effectively disabling the limiter in uniform regions, the addition of the term in Equation 2.6 greatly improves convergence. The choice of the parameter K is a compromise. Large values of K are favorable to accuracy in smooth regions and good convergence. However, since for any K > 0 the limiter no longer strictly enforces monotonicity, large values of K can lead to significant overshoots near discontinuities in the solution.  2.4 2.4.1  High-Order Limiting Monotonicity  The first challenge of extending the limiting procedure to third- and fourthorder accurate schemes is to express the monotonicity requirement including the high-order reconstruction terms. In second-order schemes the assumption that the solution at the reference point is equal to the control volume average is often made. For cell-centered schemes, where the reference point is the centroid, this assumption is correct. For vertex-centered scheme, the reference point is usually chosen as the vertex location, therefore the assumption is 30  2.4. High-Order Limiting not strictly correct. However, for third- and higher-order schemes the control volume average solution is in general not equal to the centroidal value of the reconstruction. Therefore, when devising a limiter for these schemes, making the distinction between control volume average values and reference-point values becomes critical for maintaining high-order accuracy. With this in mind, the reconstruction polynomial in Equation 2.2 can be rewritten in terms of the control volume average with the help of equation 2.3 to yield:  UiR (x − xi ) = U i +   +  ∂U ∂x  ∂2U ∂x2  ∂2U + ∂y 2  ((x − xi ) − xi ) + xi  ( xi  ∂U ∂y  ((y − yi ) − y i ) xi  (x − xi )2 x2 ∂2U − )+ 2 2 ∂x∂y  ((x − xi )(y − yi ) − xy) xi    xi  (y − yi )2 y 2 ( − ) + · · · 2 2  (2.7)  This can be interpreted as meaning that the reconstructed solution at any point is the control volume average plus second and high-order contributions from the reconstruction : UiR (x − xi ) = U i + S(x − xi ) + H(x − xi ) Therefore, analogous to the second-order case, no new extrema will be formed if δUimin ≤ S(x − xi ) + H(x − xi ) ≤ δUimax As in the work of Barth[2], the limited form of the high-order accurate reconstruction can be expressed as UiR (x − xi , Φi ) = U i + Φi S(x − xi ) + Hi (x − xi )  (2.8)  Given this formulation, the same limiting procedure used in Section 2.3 can be applied to the high-order reconstruction. Some previous works [7, 13] have suggested a formulation where the 31  2.4. High-Order Limiting limiter value multiplies only the second-order terms while the high-order terms are “switched off” when discontinuities are detected. This formulation has the form of UiR (x − xi , Φi , σi ) = U i + (Φi (1 − σi ) + σi ) S(x − xi ) + σi H(x − xi ) where σi , the discontinuity detector, is zero near discontinuities and one in smooth regions of the flow. However, this approach may violate the monotonicity requirement as the high-order terms may have contributed to reducing the overshoot in the unlimited reconstruction used in determining the value of Φi . Therefore, the value of Φi computed may be insufficient to reduce the slope such that overshoots occur when the high-order reconstruction terms are disabled. The additional parameters introduced in a smooth switching function for σi can be considered a further disadvantage.  2.4.2  Accuracy  As previously mentioned, on uniform grids a limiter for the second-order scheme maintains nominal accuracy as long as |φ − 1| ≤ O (∆x) since this results in an error that is on the order of the quadratic term in the Taylor expansion. However, when a third- or fourth-order scheme is used, the limiter must satisfy |φ − 1| ≤ O ∆x2 or |φ − 1| ≤ O ∆x3 , respectively, for the effect of the limiter in smooth regions to be on the order of truncation error. For this reason, Venkatakrishnan’s limiter will not provide sufficient accuracy even in smooth regions without any local extrema. While the BarthJespersen limiter does satisfy these conditions, the lack of differentiability will make achieving a steady-state solution difficult. Therefore, we seek a new approximation for min(1, y) used in step 3 of the limiting procedure, which we will call min(1, y). Like Venkatakrishnan’s function, given in Equation 2.5, we require that it be differentiable at all points and that it be entirely contained under the function min(1, y). However, unlike Equation 2.5, we also require this new limiting function to have a value of exactly 1 for a range of values y ≥ yt where 1 < yt < 2 represents a threshold value. For  32  2.4. High-Order Limiting  min(1,y) P(y) with yt=1.5 P(y) with yt=1.75  1.2 1 0.8 0.6 0.4 0.2 0 0  0.5  1  1.5  2 yt  2.5  3  3.5  4  Figure 2.2: min(1, y) with yt = 1.5 and yt = 1.75 compared to min(1, y) this function, we propose the form min(1, y) =  P (y) 1  y < yt y ≥ yt  where P (y) is a polynomial satisfying P |0 = 0 dP dy 0 = 1  P |yt = 1 dP dy yt = 0  P (y) ≤ min(1, y),  y ∈ [0, yt ]  The resulting polynomials for yt = 1.5, 1.75 are plotted in Figure 2.2. This function allows the preservation of high-order accuracy on uniform grids by satisfying min(1, y) − 1 ≤ O ∆x3 . Additionally, this function is also effective in maintaining high-order accuracy in regions of mild mesh nonuniformity. The degree of non-uniformity that can be accommodated is dictated by the choice of the threshold value yt . Smaller values of yt are less likely to unduly activate the limiter on non-uniform meshes but result in a limiter that approaches non-differentiability. Therefore the the choice of yt 33  2.4. High-Order Limiting is a compromise between maintaining good accuracy on non-uniform grids and maintaining good convergence properties. For the results presented in this work, we use yt = 1.5 which yields the following cubic polynomial P (y) = −  2.4.3  4 3 y +y 27  Uniform Regions and Smooth Extrema  For the reasons already mentioned in Section 2.3, it is desirable to eliminate the effect of the limiter in regions of uniform flow or near smooth extrema. To maintain high-order accuracy it is essential to permit monotonicity to be violated near smooth extrema. Furthermore the limiter value changes rapidly and is non-differentiable in uniform regions which prevents good convergence of the solver. Therefore, analogous to the second-order limiter of Venkatakrishnan, we wish to eliminate the effect of the limiter when the local solution variation is O ∆x2 or smaller. Specifically, we propose to disable the limiter when 3  δU ≡ (δUimax − δUimin ) < (K∆x) 2  where K is a tunable parameter. However using a simple switch would introduce a non-differentiable step in the residual evaluation which would cause convergence problems. Therefore, to maintain differentiability, the following procedure is proposed: Φi = σi + (1 − σi )Φi  (2.9)  where Φi is the limiter value as calculated in step 4 of the procedure in Section 2.3 and σi is the following function:  σi =          1 s  δU 2 −(K∆x)3 (K∆x)3  0  δU 2 ≤ (K∆x)3 (K∆x)3 < δU 2 < 2(K∆x)3 δU 2  ≥  (2.10)  2(K∆x)3  where the transition function s is defined by 34  2.4. High-Order Limiting  1  s(y)  0.8  0.6  0.4  0.2  0 0  0.2  0.4  0.6  0.8  1  y  Figure 2.3: Transition function s(y) = 2y 3 − 3y 2 + 1 used to smoothly disable the limiter in nearly uniform regions  s(y) = 2y 3 − 3y 2 + 1  (2.11)  and is plotted in Figure 2.3. The limited reconstruction is then computed for each Gauss point by evaluating UiR (x − xi , Φi ). Although this two stage limiting procedure is somewhat more computationally expensive than Venkatakrishnan’s limiter in the general case, some “short circuiting” is possible in uniform regions of flow. Since σi depends only on neighboring control volume averages, unlike Φi which also depends on an evaluation of the unconstrained reconstruction at each Gauss point, it is relatively inexpensive to compute. When σi evaluates to 1, computational effort can be saved by not computing Φi since it does not affect the value of the final limiter Φi .  2.4.4  Boundary Treatment  Maintaining high-order accuracy near domain boundaries represents a special challenge. Local extrema are expected to exist on the boundary for smooth 35  2.4. High-Order Limiting flows. Unlike smooth extrema in the interior of the domain, these points are not characterized by a zero first derivative of flow property with respect to space. Therefore, the method used in Subsection 2.4.3 will not be effective in disabling the limiter in these regions. As a first measure, in our implementation of both Venkatakrishnan and the new limiter we elect to only iterate over interior face Gauss points in Step 3 of the limiting procedure. Therefore an extremum forming on a boundary Gauss point will not cause the limiter to activate. In our experience this does not cause any oscillatory issues at shocks. However, the Gauss points on interior faces of boundary control volumes will often also have reconstructed values of smooth solutions that form an extremum relative to the control volume averages of the solution in the reconstruction stencil. In some cases even Gauss points of control volumes adjacent to boundary control volumes will exhibit this behavior. In our experience eliminating all of these Gauss points from Step 3 of the limiting procedure causes unacceptable oscillations in the solution near shocks. Therefore, to maintain high-order accuracy near boundaries while maintaining solution monotonicity another approach is needed. Thus far the proposed limiting procedure is entirely physics agnostic. However we have been unable to find a satisfactory means of dealing with the issue of accuracy and monotonicity of boundary and near-boundary control volumes in this manner. Therefore we propose a method specific to the Euler equations. For external and many internal flows the far-field boundary usually has nearly uniform flow conditions. For these boundaries the limiter is effectively disabled by the method applied in Subsection 2.4.3. Therefore we will focus our attention on the wall boundary. Two flow conditions can be present at a wall boundary condition: the flow can be tangential to the wall, or it can be a stagnation point. We propose methods for maintaining high-order accuracy in both cases.  36  2.4. High-Order Limiting 2.4.4.1  Tangential Flow  For every wall boundary control volume we will consider a “ghost” control volume which is a mirror image of the boundary control volume about the boundary. It will have an approximate control volume average value of the solution assigned to it consistent with a shockless flow. These values will be used to expand the stencil used to determine monotonicity outside of the flow domain. This will effectively make the reconstructed Gauss point values no longer extrema for smooth flows. To preserve high-order accuracy, these ghost control volume values will not be used in the least-squares reconstruction process. They will only be used when determining the values of δUmin and δUmax in the limiting procedure. For any control volume that includes the boundary control volume as its first neighbor the respective mirror control volume will also be included in determining δUmin and δUmax . To determine an appropriate solution value for the mirror control volume we begin by noting that any well resolved curved boundary locally resembles a circular arc. We approximate the radius of curvature from the boundary Gauss point normals of the curved boundary edges. Since this information is needed by the high-order boundary flux integration scheme, it is readily available. Next, we make the assumption of steady momentum in the direction normal to the streamline ∂P ρV 2 =− ∂n R where n is the direction normal to the wall, ρ is density, V is velocity and R the radius. Since the ghost value is only used in computing the limiter, a first-order approximation is sufficient. We can therefore approximate the ghost value of pressure by 2  P gi = P i − 2d ·  ρi V i R  (2.12)  where P i , ρi , and V i are the control volume average pressure, density and velocity respectively, and d is the distance of the control volume centroid  37  2.4. High-Order Limiting from the wall in the convex direction (d is negative for concave boundaries). The ghost values of Mach number and density can be obtained by considering the isentropic transformation from the control volume state to the ghost value state with a pressure of P gi  ρgi = ρi  P gi Pi   2  M gi =  1 γ  2  P ti γ−1 P gi  γ−1 γ    − 1  where M gi is the Mach number in the ghost control volume and P ti is the total pressure as calculated using the boundary control volume average flow properties. This, together with the assumption that the flow direction remains tangential to the surface, fully establishes the state of the ghost control volume. 2.4.4.2  Stagnation Point  Additional steps need to be taken to prevent the application of the limiter at stagnation points. Although affecting a very small fraction of control volumes, obtaining high-order accuracy near stagnation points can be critical to global accuracy. In practical aerodynamics problems, the flow discontinuity which requires the proper application of the limiter is the shock. Since supersonic flow is required to produce a shock, and stagnation points are necessarily in subsonic regions, it is possible to simply disable the limiter in the latter areas. Specifically, we propose that the reconstruction of a control volume can only be adversely affected by a shock if at least one of the control volumes in its reconstruction stencil contains supersonic flow. Practically, we wish to smoothly disable the limiter as the highest control volume average Mach number in the reconstruction stencil is reduced. For this purpose we can reuse the approach of Subsection 2.4.3. Specifically, the limiter value is once  38  2.4. High-Order Limiting again modified such that Φi = σi + (1 − σi )Φi ,  (2.13)  where Φi is the limiter value after applying the uniform-flow fix in Equation 2.9 and σi =          1 s  Mi,max −M1 M2 −M1  0  Mi,max ≤ M1 M1 < Mi,max < M2  (2.14)  Mi,max ≥ M2  where Mi,max is the maximum Mach number of the control volume averages of the reconstruction stencil of control volume i, and s is the function in Equation 2.11. The parameters M1 and M2 define the Mach numbers at which the effect of the limiter is fully disabled and fully enabled respectively. In the present work we use M1 = 0.8 and M2 = 0.85. The new limited reconstruction used by the flux integration scheme becomes UiR (x − xi , Φi ). For typical aerodynamics flows there is no harm in applying this step to all control volumes, even if they are not near boundaries. Doing so can reduce computational effort since Φi does not need to be evaluated if σi evaluates to 1.  2.4.5  Complete Algorithm  The complete algorithm is, in the general case, more complex and costly than the second-order limiting procedure presented in Section 2.3. However, the various additions can, in some cases, be used to “short circuit” the evaluation of the limiter, therefore reducing computational effort. The complete algorithm for applying the limiter to the reconstruction for each flow property of each control volume is 1. Find the maximum Mach number of the reconstruction neighbors control volume averages of the solution and evaluate σi using Equation 2.14. If σi = 1 then Φi = 1 and the algorithm jumps to step 9. 2. Find the largest negative (δUimin = min(U j −U i )) and positive (δUimax = max(U j − U j )) difference between the solution in the immediate neigh39  2.5. Results bors j and the current control volume i. If the control volume or any of its immediate neighbors are adjacent to a wall boundary include their ghost values using Equation 2.12 and isentropic relations. 3. Compute σi using Equation 2.10. If σi = 1 then Φi = Φi = 1 and the algorithm jumps to step 9. 4. Compute the unconstrained reconstructed value at each Gauss point (Uik = UiR (xk − xi )). 5. Compute a maximum allowable value of φik for each Gauss point k.  φik =   δU max  min 1, U i−U ,   i ik  min      δUi ik −U i  if Uik − U i > 0  min 1, U  , if Uik − U i < 0  1,  if Uik − U i = 0  6. Select Φi = min(φik ) . 7. Compute Φi using Equation 2.9. 8. Compute Φi using Equation 2.13. 9. Compute the limited reconstruction UiR (x − xi , Φi ) at flux Gauss integration points using Equation 2.8.  2.5  Results  The presented results were obtained using a Newton-GMRES [17] vertexcentered finite-volume solver. The solution process consists of two stages. In the preiteration stage the linear system resulting from a local timestepping is solved at each iteration. The Jacobian from the first-order accurate scheme is used on the left-hand side and the full-order accurate flux is used on the right-hand side. At each Newton iteration, the linear system is approximately solved using incomplete-lower-upper factorization (ILU) preconditioned GMRES. During the second stage, the left-hand side is  40  2.5. Results  Figure 2.4: Two coarsest meshes used for Ringleb’s flow test case replaced with the full-order accurate Jacobian as described in the manuscript reproduced in Chapter 3. In addition to results using the new limiter, results are also presented for the high-order scheme using the procedure in Section 2.4.1 but with Venkatakrishnan’s limiting function. For Venkatakrishnan’s function we use a tuning parameter of K = 2 and for the new limiter we use K = 1.  2.5.1  Ringleb’s Flow  We begin by considering Ringleb’s flow which is transonic but shockless and has a known exact solution. This will enable us to quantify the negative effects of Venkatakrishnan’s limiter and the new limiter on the accuracy of the solution in smooth regions of flow. We consider four meshes consisting of 369, 1426, 5467 and 20690 control volumes. The two coarsest meshes are shown in Figure 2.4. The exact solution to Ringleb’s flow is given by the streamlines defined  41  2.5. Results by: x =  11 2ρ  y = ±  1 kρq  1 2 − 2 2 q k 1−  + q k  J 2  2  where k is constant along streamlines and q =  V  1 1 1 1+c 1 + 3 + 5 − ln c 3c 5c 2 1−c γ−1 2 c = 1− q 2 ρ = c2/(γ−1)  J  =  The computational domain is constructed using solid walls at the streamlines for k = 0.55, 1.2 and placing the inlet and outlet at q = 0.5. This results in a geometry that is symmetric about the horizontal axis. The flow is subsonic at the inlet and outlet but is supersonic near the inside wall of the throat. Using this domain, rather than just the upper half, is a more stringent test of the accuracy of the computational scheme since errors produced at the throat are allowed to propagate to the lower half of the domain. Since this test case contains no stagnation points, the Mach number dependent deactivation of the limiter should not be necessary. For this reason, and to make the test as stringent as possible, we do not apply steps 1 and 8 of the high-order limiting procedure as presented in Subsection 2.4.5. As an initial qualitative assessment, we compare the Mach number contours generated on the 1426 control volume mesh using Venkatakrishnan’s limiter and the new limiter in Figure 2.5. The additional dissipation of the Venkatakrishnan limited solutions is evident in the reduced Mach number at the inner wall near the outflow. Using the same criterion we can visually detect a slight superiority of the fourth-order scheme over the second-order scheme when using the new limiter. The Venkatakrishnan limited procedure, on the other hand, does not benefit from the use of the high-order accurate 42  2.5. Results  (a) 2nd-order with Venkatakrishnan limiter  (b) 2nd-order with new limiter  (c) 4th-order with Venkatakrishnan limiter  (d) 4th-order with new limiter  Figure 2.5: Mach number contours from numerical solution of Ringleb’s flow on 1426 control volume mesh 43  2.5. Results reconstruction. Next we examine the value of the limiter value Φ for the pressure component at steady state for the fourth-order solution on the 5467 control volume mesh using the two limiting schemes in Figure 2.6. In addition to Venkatakrishnan’s limiter and the new limiter, we also consider the new limiter without the special boundary ghost value treatment introduced in steps 3 and 7 of the high-order limiting procedure. The new limiter successfully avoids limiting in almost all control volumes. All the control volumes that are limited have values Φ ≥ 0.99. Without the boundary ghost treatment significant limiting occurs near the inner wall boundary. Venkatakrishnan’s limiter, on the other hand, applies some limiting to virtually all control volumes. Next, in Figure 2.7, entropy is plotted for the second- and fourth-order schemes on the second-finest grid. The new limiter causes no distinguishable additional entropy production relative to the unlimited scheme for the second-order method, and a very slight increase for the fourth-order method. On the other hand, Venkatakrishnan’s limiter increases entropy by approximately an order of magnitude for the second-order scheme and two to four orders for the fourth-order scheme. Once again we note that when applying Venkatakrishnan’s limiter there is no apparent benefit to using the fourth-order scheme over the second-order scheme. The L2 norm of error of the converged solution compared to the known exact solution is used to generate Figure 2.8. The grid-convergence orders of the error norms using a regression fit of the results from all grids are given in Table 2.1. We begin by noting that the unlimited schemes attain their nominal orders of accuracy in L1 and L2 norms and one order accuracy less than nominal for the L∞ norm. Although the error of second-, third-, and fourth-order schemes using Venkatakrishnan’s limiter is lower than that of the first-order scheme, the error does not converge with nominal accuracy. In fact the grid convergence of the error when using Venkatakrishnan’s limiter is only first-order. The new limiter, on the other hand, allows for the secondand third-order schemes to perform virtually identically to the unlimited case and has only minimal adverse effect on the accuracy of the fourth-order scheme. The results using the new limiter without the boundary ghost value 44  2.5. Results  (a) Venkatakrishnan’s limiter  (b) New limiter without boundary ghost values  (c) New limiter  Figure 2.6: Limiter value for pressure for the fourth-order converged solution on the 5467 control volume mesh. Only values Φ = 1 are plotted  45  2.5. Results  (a) 2nd-order no limiter  (b) 4th-order no limiter  (c) 2nd-order Venkatakrishnan limiter  (d) 4th-order Venkatakrishnan limiter  (e) 2nd-order new limiter  (f) 4th-order new limiter  Figure 2.7: Difference in dimensionless entropy from the freestream value for Ringleb’s flow 46  2.5. Results 0.1  L2 Norm Error in P  0.01  0.001  0.0001  1e-05  1e-06 200  1st-Order 2nd-Order, No Limiter 2nd-Order, Venkatakrishnan Limiter 2nd-Order, New Limiter w/o Ghosts 2nd-Order, New Limiter 400  800  1600 3200 Num CVs  6400  12800  25600  6400  12800  25600  6400  12800  25600  (a) Second-order 0.01  L2 Norm Error in P  0.001  0.0001  1e-05  1e-06 200  2nd-Order, No Limiter 3rd-Order, No Limiter 3rd-Order, Venkatakrishnan Limiter 3rd-Order, New Limiter w/o Ghosts 3rd-Order, New Limiter 400  800  1600 3200 Num CVs  (b) Third-order 0.001  L2 Norm Error in P  0.0001  1e-05  1e-06  1e-07 200  2nd-Order, No Limiter 4th-Order, No Limiter 4th-Order, Venkatakrishnan Limiter 4th-Order, New Limiter w/o Ghosts 4th-Order, New Limiter 400  800  1600 3200 Num CVs  (c) Fourth-order  Figure 2.8: Error convergence for pressure in Ringleb’s flow 47  2.5. Results Nominal Order 1st 2nd 2nd 2nd 2nd 3rd 3rd 3rd 3rd 4th 4th 4th 4th  Limiter  L1 Norm  L2 Norm  L∞ Norm  None None Venkatakrishnan New w/o Ghosts New None Venkatakrishnan New w/o Ghosts New None Venkatakrishnan New w/o Ghosts New  1.24 2.22 1.10 1.71 2.22 3.18 1.14 1.83 3.14 4.07 0.62 1.07 3.24  1.24 1.96 1.07 1.23 1.96 3.12 1.10 1.43 3.08 3.74 0.62 0.74 3.11  0.99 1.24 0.47 0.47 1.24 2.70 0.51 0.59 2.40 3.06 0.12 0.23 2.08  Table 2.1: Convergence order of norms of error in pressure for Ringleb’s flow computed using regression fit of all mesh results. treatment outperform Venkatakrishnan’s limiter only slightly. This indicates that the application of the limiter in even an isolated region has severe implications for global accuracy.  2.5.2  Transonic Flow Over an Airfoil  Next, we present results for transonic flow over a NACA 0012 airfoil at Mach 0.8 and an angle of attack of 1.25 degrees. The computational mesh consists of 4656 control volumes and is shown in Figure 2.9. We will consider second- and fourth-order schemes using Venkatakrishnan’s limiter, the new limiter, and the new limiter without the stagnation region fix which disables the limiter at low Mach number in steps 1 and 8 of the high-order limiting procedure.  48  2.5. Results  Figure 2.9: Mesh consisting of 4656 control volumes used for the NACA 0012 test case  49  2.5. Results 2.5.2.1  Accuracy  We begin by assessing the quality of the solution upstream of the shock. Entropy near the leading edge of the airfoil is plotted in Figure 2.10. The fourth-order solution with Venkatakrishnan’s limiter once again fails to outperform its second-order counterpart. The new limiter produces approximately an order of magnitude less entropy than Venkatakrishnan’s limiter for the second-order scheme. The new limiter applied to the fourth-order scheme results in even less entropy production. Disabling the stagnation region fix results in a modest increase in entropy production. The quality of the solutions can also be compared by examining the stagnation pressure ratio along the upper surface of the airfoil shown in Figure 2.11. The decrease in total pressure across the shock is comparable for all schemes. However, the schemes limited with the Venkatakrishnan limiter result in substantial stagnation pressure loss upstream of the shock. For the second-order scheme the stagnation region fix has little effect while for the fourth-order scheme applying this step in the limiting procedure results in a further improvement in the conservation of total pressure. 2.5.2.2  Shock Monotonicity  The performance of the new limiter in controlling oscillations and overshoots in the solution near the strong shock on the upper surface of the airfoil is demonstrated in the pressure plot of Figure 2.12. The new limiter and Venkatakrishnan’s limiter are both effective in eliminating any substantial oscillations near the shock. The pressure computed on the upper surface of the airfoil using the different limiters is virtually indistinguishable. Once again the lower dissipation of the new limiter is demonstrated by the sharper profile of the weak shock on the lower surface of the airfoil. 2.5.2.3  Convergence  Next, we consider the residual convergence properties of the new limiting scheme coupled with our Newton-GMRES solver. Figure 2.13 shows the  50  2.5. Results  (a) Second-order Venkatakrishnan limiter (b) Fourth-order Venkatakrishnan limiter  (c) Second-order new limiter without stag- (d) Fourth-order new limiter without stagnation fix nation fix  (e) Second-order new limiter  (f) Fourth-order new limiter  Figure 2.10: Difference in dimensionless entropy from the freestream value for Mach 0.8, α = 1.25 flow around a NACA 0012 airfoil  51  2.5. Results  0.06  2nd-order Venkatakrishnan 2nd-order New Limiter w/o Stagnation Fix 2nd-order New Limiter  0.05  1-(Pt/Pt_inf)  0.04 0.03 0.02 0.01 0 -0.01 0  0.06  0.1  0.2  0.3  0.4  0.5 0.6 x/c (a) Second-order  0.7  0.8  0.9  1  0.9  1  4th-order Venkatakrishnan 4th-order New Limiter w/o Stagnation Fix 4th-order New Limiter  0.05  1-(Pt/Pt_inf)  0.04 0.03 0.02 0.01 0 -0.01 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  x/c (b) Fourth-order  Figure 2.11: Decrease in total pressure along the upper surface of the NACA 0012 airfoil at Mach 0.8 α = 1.25  52  2.5. Results  -1.5 -1  Cp  -0.5 0 0.5 1 2nd-order Venkatakrishnan 2nd-order New Limiter  1.5 0  0.2  0.4  0.6  0.8  1  x/c (a) Second-order  -1.5 -1  Cp  -0.5 0 0.5 1 4th-order Venkatakrishnan 4th-order New Limiter  1.5 0  0.2  0.4  0.6  0.8  1  x/c (b) Fourth-order  Figure 2.12: Surface pressure profiles for transonic flow over a NACA0012 airfoil  53  2.5. Results  2nd-order Venkatakrishnan 2nd-order New 4th-order Venkatakrishnan 4th-order New  10000  Residual  100 1 0.01 0.0001 1e-06 1e-08 0  20  40  60 Iterations  80  100  120  Figure 2.13: Convergence history for transonic airfoil case Order 2nd 2nd 4th 4th  Limiter Venkatakrishnan New Venkatakrishnan New  Computational Time (s) 37 40 84 99  Table 2.2: Computational time for transonic airfoil test case convergence of the scheme with the new limiter relative to Venkatakrishnan’s limiter for second and fourth order accurate schemes. We have also provided the computational time required for convergence in Table 2.2. The scheme with the new limiter exhibits a slightly poorer convergence rate than Venkatakrishnan’s limiter. This is likely due to the lower dissipation of the scheme which results in poorer conditioning of the nonlinear system of equations. Similarly, we find that the fourth-order scheme requires slightly more iterations to converge than the second-order scheme.  54  2.5. Results  -1.2 -1.15 -1.1  Cp  -1.05 -1 -0.95 -0.9 K=0.5 K=2 K=8  -0.85 -0.8  0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 x/c (a) Venkatakrishnan limiter -1.2 -1.15 -1.1  Cp  -1.05 -1 -0.95 -0.9 K=0.25 K=1 K=4  -0.85 -0.8  0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 x/c (b) New limiter  Figure 2.14: Upper surface pressure profile upstream of the shock for transonic flow over a NACA0012 airfoil using fourth-order scheme  55  2.6. Conclusion 2.5.2.4  Sensitivity to Tuning Parameter  Venkatakrishnan’s limiter and the new limiter both require a tuning parameter K to detect regions of uniform flow. Ideally the solutions produced by the schemes should be relatively insensitive to the choice of this parameter. Therefore we consider the effect of increasing and decreasing the parameter by a factor of 4 from the original choice. For Venkatakrishnan’s scheme we consider K = 0.5, 2, 8 and for the new limiter we consider K = 0.25, 1, 4. In Figure 2.14 we examine the behaviour of pressure on the upper surface of the airfoil just upstream of the shock using the fourth-order scheme. The solution using Venkatakrishnan’s limiter exhibits a strong sensitivity to the value of the parameter. Specifically, K = 0.5 results in a lower peak pressure ahead of the shock, while K = 8 results in significant oscillations. The result using the new limiter shows less sensitivity. Specifically, the results using K = 0.25 are virtually indistinguishable from those with K = 1 while the results with K = 4 result in some oscillation. To further examine the effect of this parameter on the dissipation of the scheme, the stagnation pressure along the upper surface of the airfoil is shown in Figure 2.15. For Venkatakrishnan’s scheme the value of the tuning parameter has little effect on the total pressure loss near the leading edge of airfoil. However the stagnation pressure loss across the shock exhibits a strong dependence on this parameter. The new limiter exhibits almost no sensitivity to the parameter ahead of the shock. Downstream of the shock the results using parameters K = 0.25 and K = 1 are indistinguishable while the stagnation pressure with K = 4 is higher.  2.6  Conclusion  A new unstructured grid limiter broadly based on the framework introduced by the Barth and Jespersen limiter has been developed specifically for high-order schemes. Like the Venkatakrishnan limiter, it is designed to be differentiable to enable good convergence to steady state by implicit schemes. To achieve this and to maintain high-order accuracy, a new function defining  56  2.6. Conclusion  0.06  K=0.5 K=2 K=8  0.05  1-(Pt/Pt_inf)  0.04 0.03 0.02 0.01 0 -0.01 0  0.1  0.2  0.3  0.4  0.5 0.6 0.7 x/c (a) Venkatakrishnan limiter  0.06  0.8  0.9  1  0.9  1  K=0.25 K=1 K=4  0.05  1-(Pt/Pt_inf)  0.04 0.03 0.02 0.01 0 -0.01 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  x/c (b) New limiter  Figure 2.15: Sensitivity of total pressure along the upper surface of the NACA 0012 airfoil to the choice of limiter tuning parameter  57  2.6. Conclusion the relation between the limiter value and the unlimited reconstruction at Gauss points has been constructed. To preserve accuracy at smooth extrema, a smooth switch disables the limiter in regions of nearly uniform flow. To avoid limiting regions of smooth flow near boundaries, the range defining monotonicity is carefully expanded based on a constant momentum extrapolation of conditions beyond the boundary. Application of the limiter is avoided at stagnation points by smoothly disabling it in regions of low Mach number. The results indicate that the new limiter is effective in maintaining high-order accuracy in smooth flows while effectively suppressing oscillations near shocks. Additionally, the new scheme is also relatively insensitive to tuning parameters which should lead to robustness and reproducibility of solutions. Although the greatest gains in accuracy are seen when using a highorder scheme, the new limiter is also effective in reducing the dissipation and increasing the accuracy of second-order schemes relative to the currently most widely used limiter. Since only minor modifications to the reconstruction procedure are needed, the addition of the new method to existing solvers should be relatively straightforward. Although the present work presents only inviscid results in two dimensions, the extension to viscous flows and three dimensions should be relatively simple. Specifically, for viscous flows the stagnation region fix should be sufficient in eliminating the effect of the limiter near walls making the wall ghost value fix unnecessary. For three-dimensional inviscid flows the wall ghost value fix will need to account for streamline curvature even though streamlines will not be coordinate aligned.  58  2.7. Bibliography  2.7  Bibliography  [1] R. Abgrall. On essentially non-oscillatory schemes on unstructured meshes: Analysis and implementation.  Journal of Computational  Physics, 114(1):45–58, 1994. [2] T. J. Barth. Recent developments in high order k-exact reconstruction on unstructured meshes. AIAA paper 93-0668, Jan. 1993. [3] T. J. Barth. Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations. In Lecture Series 1994-05, Rhode-Saint-Gen`ese, Belgium, Mar. 1994. von Karman Institute for Fluid Dynamics. [4] T. J. Barth and P. O. Frederickson. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA paper 90-0013, Jan. 1990. [5] T. J. Barth and D. C. Jespersen. The design and application of upwind schemes on unstructured meshes. AIAA paper 89-0366, Jan. 1989. [6] S. De Rango and D. W. Zingg. Higher-order spatial discretization for turbulent aerodynamic computations. American Institute of Aeronautics and Astronautics Journal, 39(7):1296–1304, July 2001. [7] M. Delanaye and J. A. Essers. Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids. American Institute of Aeronautics and Astronautics Journal, 35(4):631–639, Apr. 1997. [8] O. Friedrich. Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids. Journal of Computational Physics, 144(1):194–212, July 1998. [9] P. Geuzaine, M. Delanaye, and J.-A. Essers. Computation of high Reynolds number flows with an implicit quadratic reconstruction scheme  59  2.7. Bibliography on unstructured grids. In Proceedings of the Thirteenth AIAA Computational Fluid Dynamics Conference, pages 610–619. American Institute of Aeronautics and Astronautics, 1997. [10] A. G. Godfrey, C. R. Mitchell, and R. W. Walters. Practical aspects of spatially high-order accurate methods. American Institute of Aeronautics and Astronautics Journal, 31(9):1634–1642, Sept. 1993. [11] C. Q. Hu and C. W. Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150(1):97–127, Mar. 1999. [12] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted ENO schemes. Journal of Computational Physics, 126(1):202–228, June 1996. [13] A. Nejat and C. Ollivier-Gooch. A high-order accurate unstructured finite volume Newton-Krylov algorithm for inviscid compressible flows. Journal of Computational Physics, 227(4):2592–2609, 2008. [14] C. F. Ollivier-Gooch. High-order ENO schemes for unstructured meshes based on least-squares reconstruction. AIAA paper 97-0540, Jan. 1997. [15] C. F. Ollivier-Gooch. Quasi-ENO schemes for unstructured meshes based on unlimited data-dependent least-squares reconstruction. Journal of Computational Physics, 133(1):6–17, 1997. [16] C. F. Ollivier-Gooch and M. Van Altena. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. Journal of Computational Physics, 181(2):729–752, 2002. [17] Y. Saad and M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal of Scientific and Statistical Computing, 7(3):856–869, July 1986. [18] C.-W. Shu, T. A. Zang, G. Erlebacher, D. Whitaker, and S. Osher. Highorder ENO schemes applied to two- and three-dimensional compressible flow. Applied Numerical Mathematics, 9(1):45–71, 1992. 60  2.7. Bibliography [19] G. D. van Albada, B. van Leer, and W. W. Roberts, Jr. A comparative study of computational methods in cosmic gas dynamics. Astronomy and Astrophysics, 108:76–84, Apr. 1982. [20] B. van Leer. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. Journal of Computational Physics, 32:101–136, 1979. [21] V. Venkatakrishnan. On the accuracy of limiters and convergence to steady-state solutions. AIAA paper 93-0880, Jan. 1993. [22] D. Zingg, S. De Rango, M. Nemec, and T. Pulliam. Comparison of several spatial discretizations for the Navier-Stokes equations. Journal of Computational Physics, 160:683–704, 2000.  61  Chapter 3  Efficient Convergence to Steady State† 3.1  Introduction  High-order discretizations on structured grids have been shown [6, 24] to reduce computational effort for a given level of solution accuracy. Although high-order finite-volume methods on unstructured grids are relatively well known [1, 2, 7], means to efficiently converge the solution to steady state have not been adequately studied. For this reason, the advantages of high-order methods have not been convincingly demonstrated. Although Newton-GMRES [21] methods provide rapid convergence, they are known to be less robust than explicit schemes. This lack of robustness is directly related to the choice of globalization technique. For compressible flows, the technique of choice has been the addition of a timestep that is increased as the solution converges. Successful convergence often depends on the careful choice of parameters controlling the growth of this timestep. High-order schemes, due to their lower diffusivity, are particularly sensitive to these robustness issues. In this work, we demonstrate how the timestep approach can be improved by the addition of a line search method. We also demonstrate that scaling the local timestep based on the local residual can dramatically improve performance and robustness. The use of the GMRES method to approximately solve the linear systems arising from implicit schemes is an increasingly popular approach for solution †  A version of this chapter will be submitted for publication. C. Michalak and C. Ollivier-Gooch. Globalized matrix-explicit Newton-GMRES for the high-order accurate solution of the Euler equations.  62  3.1. Introduction convergence for second-order unstructured flow solvers. Previous work has focused on using the matrix-free [3] approach. With this method the full flux Jacobian does not need to be explicitly constructed. Instead, the matrix-vector products required by GMRES are computed using Frechet derivatives. However, a preconditioning matrix is still required for good GMRES convergence; the typical choice here is to compute the Jacobian of the first-order scheme (though perhaps using the reconstructed second-order data) and use an incomplete LU factorization of this matrix to precondition the linear system. Our group’s previous work in developing an efficient high-order accurate flow solver extended this approach to high-order residual calculations, including careful study of preconditioning methods [15, 16]. While those studies were successful in dramatically reducing CPU time requirements for high-order schemes, especially for third-order accuracy, there are also clear indications that the fourth-order scheme in particular is poorly preconditioned by the first-order Jacobian. The research described in this paper explores the possibility of forming and exploiting the full order Jacobian. On the one hand, we expect that preconditioning will be much better if we use the full order Jacobian. Also, having the Jacobian on hand will eliminate all the flux evaluations done within GMRES inner iterations in matrix-free schemes; this savings in CPU time should help offset the cost of computing the Jacobian. The major drawback to forming the full order Jacobian is the memory requirements. Section 3.2 briefly reviews the high-order discretization used in the present work, focusing on the reconstruction step. In Section 3.3 globalization techniques for the Newton method are discussed. In this section we discuss how the line search method can be effectively combined with the timestepping method. We also present our novel method of deriving a local timestep from the residual. In Section 3.4 we discuss the matrix-free and matrix-explicit GMRES methods and their preconditioning. Some detail on how to efficiently form the exact Jacobian needed by the matrix-explicit method is provided. Numerical results for a multi-element subsonic test case and a single element transonic test case are provided in Section 3.5. Concluding remarks are found in Section 3.6. 63  3.2. High-Order Accurate Finite-Volume Method  3.2  High-Order Accurate Finite-Volume Method  The high-order accurate finite-volume discretization relies on a piecewise polynomial reconstruction in each control volume to provide an accurate approximation to the solution at Gauss integration points. A series of Riemann problems for the Euler equation are solved at these points using the Roe scheme [20]. This results in a discrete flux integral for each control volume which must be reduced to zero at steady state. In the present work, the control volumes are formed using a vertex-centered median dual of an unstructured triangular grid.  3.2.1  Governing Equations  The non-dimensionalized two-dimensional Euler equations for a control volume Vi can be written in integral form as d dt   ˛  ˆ  ρ  U dA + Ωi         ρu    U =   ρv     E  F dS = 0  Vi  ρun    ρuu + P n n x  F =  ρvun + P ny   (E + P )un  (3.1)         (3.2)  where U is the solution vector in conservative variables, and F is the flux vector normal to the control volume boundary. The total energy per unit volume, E, and the pressure, P , are related by the ideal gas equation of state P = ρRT  3.2.2  (3.3)  Finite-Volume Formulation  The governing equations are discretized in space using the finite-volume method on an unstructured grid. The control volume average of the solution  64  3.2. High-Order Accurate Finite-Volume Method in each control volume is defined by 1 Ui = Ai  ˆ U dA Vi  which can be used to rewrite Equation 3.1 as ˛ dU i 1 =− F dS ≡ −Ri U dt Ai Ω i  (3.4)  Since in the present work we are concerned only with steady-state solutions the time derivative term can be dropped. Therefore the system of equations to be solved is simply R(U ) = 0 where R is the vector function representing the flux integral, or residual, in all control volumes.  3.2.3  Solution Reconstruction  To obtain a high-order accurate numerical approximation to the residual we reconstruct a polynomial representation of the solution within each control volume. As in previous work [17], the high-order accurate reconstruction is obtained by formulating and solving a least-squares problem. To obtain an accurate estimate of the solution at flux quadrature points, we represent the solution within each control volume by the Taylor series expansion UiR (x − xi , y − yi ) =  U |i +  ∂U ∂U (x − xi ) + (y − yi ) + ∂x i ∂y i  ∂2U ∂x2  (x − xi )2 ∂2U + (x − xi ) (y − yi ) + 2 ∂x ∂y i  i  ∂2U ∂y 2  i  (y − yi )2 + ··· 2  65  3.2. High-Order Accurate Finite-Volume Method where Ui is the value of the reconstructed solution and  ∂ k+l Ui ∂xk ∂y l  are its deriva-  tives at the reference point (xi , yi ) of control volume i. The coefficients of the polynomial are computed so that the mean value of the solution in the control volume is conserved and the reconstruction approximates nearby control volume averages. Conservation of the mean within a control volume is satisfied when ˆ 1 ∂U ∂U Ui = UiR dA ≡ U |i + xi + y + Ai Vi ∂x i ∂y i i ∂2U ∂x2 ∂2U ∂y 2 where xn y m  i  1 ≡ Ai  i  i  x2 i ∂2U + xy + 2 ∂x ∂y i i y2i + ··· 2  (3.5)  ˆ (x − xi )n (y − yi )m dA Vi  As we shall see, this mean constraint can be eliminated analytically from the least-squares system. The terms used to form the reconstruction least-squares problem have the form 1 Aj  ˆ UiR (x − xi ) dA = Vj  U |i + +  ∂2U ∂x2  ∂U ∂U yij xij + ∂x i ∂y i  i  x2 ij ∂2U ∂2U + xy ij + 2 ∂x ∂y i ∂y 2  (3.6)  i  y 2 ij + ··· 2  The geometric terms in this equation are of the general form xn y m ij  ≡  1 Aj m  =  ˆ ((x − xj ) + (xj − xi ))n · ((y − yj ) + (yj − yi ))m dA Vj n  n! m! (xj − xi )k · (yj − yi )l · xn−k y m−l j l! (m − l)! k! (n − k)! l=0 k=0  The resulting least-squares problem to be solved for each solution variable in  66  3.2. High-Order Accurate Finite-Volume Method each control volume takes the for    xi  yi  x2 i  xy i  y2i  ···  wi1 xi1  wi1 yi1  wi1 x2 i1  wi1 xy i1  wi1 y 2 i1  ···  1    w  i1   wi2    wi3   .  ..   wiN    wi2 xi2  wi2 yi2  wi3 xi3 .. .  wi3 yi3 .. .  wiN xiN  wiN yiN  Ui  wi2  x2  wi3  x2  i2  wi2 xy i2  .. .  wi3 xy i3 .. .  wiN x2 iN  wiN xy iN  i3  wi2  y2  i2  wi3  y2  i3  ···  .. .  ··· .. .  wiN y 2 iN  ···                   U    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x ∂y 1 ∂2U 2 ∂y 2                  .. .  i       w U  i1 1      wi2 U 2    =   wi3 U 3      ..   .    (3.7)  wiN U N  where the weights can be used to emphasize geometrically nearby data: wij =  1 |xj − xi |n  (3.8)  where typically n ∈ [0, 2]. The first equation in this linear system is the mean constraint, which can be removed by Gauss elimination, leaving an unconstrained least-squares problem. In previous work [14], the least-squares problem was solved at each flux evaluation using QR factorization [9]. However since the matrix contains only geometric terms, is identical for each solution variable in a given control volume and does not change between iterations, substantial savings in computational time can be achieved by precomputing and storing the pseudoinverse of the reconstruction matrix for each control volume. To obtain the pseudoinverse in a numerically stable manner, the present work uses the singular value decomposition (SVD) method [9]. Given a SVD of a  67  3.2. High-Order Accurate Finite-Volume Method reconstruction matrix A, the pseudoinverse A† can easily be obtained A = U ΣV T A† = V Σ † U T  (3.9)  where the diagonal entries of Σ are the singular values of A, the columns of U and V are the left- and right-singular vectors, Σ† is a diagonal matrix containing the reciprocal values of Σ. Once this precomputation has been performed, the reconstruction coefficients in each control volume can be obtained by performing a matrix-vector multiplication; the number of columns in this matrix equals the number of control volumes in the reconstruction stencil, while the number of rows equals the number of required reconstruction coefficients. From these coefficients the reconstructed values of the solution can easily be computed at the flux quadrature points. When solving flows with shocks, a slope limiter is applied to the reconstruction polynomial. We refer the reader to the manuscript reproduced in Chapter 2 for details of the high-order limiting algorithm used for the transonic results of the present work.  3.2.4  Numerical Flux Calculation and Integration  The control volume boundaries are formed by connecting the cell centroids with the face centroids surrounding vertices of the triangular grid. The contour integral in Equation 3.4 is discretized using Gauss quadrature over the dual edges. One Gauss point per edge is used for the second-order scheme. Two Gauss points are used for the third- and fourth-order schemes. The procedure to compute the numerical flux F is: 1. Translate the control volume average solution from conserved variables (density, momentum, energy) to primitive variables (density, velocity, and pressure). 2. Compute the reconstruction coefficients using the precomputed pseudoinverses of the reconstruction matrices and, if necessary, the limiter. 68  3.3. Globalized Newton Method 3. At each Gauss point, compute the left and right solution states using the reconstruction coefficients from both adjacent control volumes. 4. At each Gauss point, approximately solve the Riemann problem using the reconstructed solution from both adjacent control volumes. In the present work, the Roe scheme [20] is used.  3.3  Globalized Newton Method  Once a procedure to evaluate the discrete numerical residual has been established the remaining problem is to find the root of the system of equations R(U ) = 0  (3.10)  When a good starting guess for U is known, the Newton method can be used. This results in an iterative method of the form ∂R n δU = −R U ∂U n+1 n U = U + δU  (3.11)  However, in practical problems a sufficiently good initial solution for U is not known. Therefore the Newton method must be modified to improve its robustness. Two methods are commonly used for this in computational fluid dynamics. For aerodynamics flows most authors use a form of pseudotimestepping [10, 12]. For incompressible flow a line search algorithm is sometimes used [18, 22]. In the present work we modify the Newton method by applying both a pseudo-timestep and a line search algorithm.  3.3.1  Pseudo-Timestepping  In pseudo-timestepping, instead of directly solving Equation 3.10 using the Newton method, we consider the series of root-finding problems defined by the implicit timestep U  n+1  n  = U − diag(∆t) R(U  n+1  ) 69  3.3. Globalized Newton Method n  where U is the solution from the previous timestep, U  n+1  is the solution at  the next timestep, and ∆t is a vector of local timesteps. This results in a root-finding problem at each timestep which we will define as R U  n+1  ≡R U  n+1  +  I n+1 n U −U =0 ∆t  (3.12)  The addition of a timestep component has two related benefits. First, n  the root of the new residual R is closer to the to previous iterate U which makes it a better candidate for the Newton method. Second, since it mimics the time derivative in the original time-dependent Euler equations, it causes the transient solution to evolve in a physically feasible way. This prevents physically unfeasible states, such as negative pressure or density, which result in an undefined residual, from occurring. The most common approach is to approximately solve R at each timestep using a single Newton step and to increase ∆t as the solution progresses towards steady state. In this case the scheme reduces to the backward-Euler time advance scheme ∂R I δU + ∆t ∂U U  3.3.2  = −R U  n+1  n  (3.13)  n  = U + δU  Line Search  To add robustness, instead of solving Equation 3.12 using a pure Newton method we apply a single iteration of a line search algorithm [8]. The resulting solution update is given by U  n+1  n  = U + αδU  where δU is computed using a Newton iteration and α is a positive scalar step length. To select the value of α we consider an optimization problem related to Equation 3.12. Specifically, the objective function z (x) =  1 R (x) 2  2 2  70  3.3. Globalized Newton Method is required to satisfy a sufficient decrease condition given by z U  n+1  ≤z U  n  + c1 α∇z T δU  (3.14)  where c1 is taken to be 10−4 in the present work. The product of the gradient of the objective function with the Newton update vector is simply ∇z T δU = RT  n  ∂R ∂U  δU n  and is therefore easily obtained. Satisfying this inequality, known at the Armijo condition, is necessary for the good convergence of the algorithm. If the initial guess α = 1 does not satisfy this condition a quadratic line search algorithm [8] is used to find a value α ∈ [0, 1] which does. Although the details of this algorithm are beyond the scope of this article, it should be noted that in the case α = 1 additional evaluations of the objective function, and therefore of the flux, are needed. Since the ultimate goal of the scheme is to reduce the residual R, we are not particularly concerned with the accuracy of the solutions to the intermediate problem R U  n+1  = 0. Therefore only a single iteration of  the line search algorithm is applied per timestep. As long as the residual is sufficiently well behaved for the full Newton step to satisfy the Armijo condition the result of this algorithm is the same as the backward-Euler time advance scheme with only a minimal additional expense of verifying that the sufficient decrease condition is satisfied. On the other hand, when the line search determines that backtracking is necessary the additional flux evaluations involved are a reasonable cost for the improved stability of the resulting modified step.  3.3.3 3.3.3.1  Selecting the Timestep Global timestepping  In a nonlinear sense, the timestep serves to improve the convexity of the modified residual in Equation 3.12 so that the Newton method can be applied. 71  3.3. Globalized Newton Method In a linear sense it modifies the left-hand side of Equation 3.13 by making it more diagonally dominant. These modifications are most important when the current solution is far from steady state and therefore the linearized residual at the current state does not accurately reflect the behavior of the residual at the steady state solution. Therefore, since the linearization is expected to improve as the solution converges, it is common to use the norm of the residual to control the size of the timestep. The most common algorithm, named successive evolution-relaxation (SER) [13], is defined as  τ  n+1  n  = c2 · τ ·  R(U  n−1  )  2  n  R(U )  (3.15)  2  where τ is a scalar timestep. For global timestepping we set ∆ti = τ for all control volumes. Although this method has been shown to be effective, it requires careful tuning of the initial timestep and the parameter c2 to prevent the solution process from diverging. Conversely, if c2 or the initial timestep are taken to be too small the number of iterations needed to converge may be unnecessarily large. The robustness of this approach can be improved by reducing the timestep when a solution update would result in unphysical flow properties [11]. Specifically, when we consider SER in the present work, we reduce the timestep by a factor of two if the solution update would result in the reconstructed values of density or pressure becoming negative at any of the Gauss integration points. The selection of α satisfying the condition in Equation 3.14 is a measure of how “well behaved”, in the sense of convexity, the modified residual of Equation 3.12 is. Since the use of a sufficiently small timestep has the effect of improving the behaviour of the modified residual function, α can be used to determine the quality of the previous timestep. Therefore we propose to use the result of the line search to help control the timestep such that  τ  n+1  n  = c2 · τ ·  R(U  n−1  )  2  n  R(U )  (3.16)  2  72  3.3. Globalized Newton Method where c2 (α) =  c2 ,  if α = 1  max (α, 0.5) , if α = 1  This modification essentially detects when the timestep has been increased too rapidly and decreases the size of the subsequent timestep. Since this feedback loop prevents the timestep from growing too quickly it is possible to be more aggressive with the choice of c2 than with the standard SER method. 3.3.3.2  Local Timestepping  Since we do not seek a time-accurate solution, it is not necessary to use a scalar timestep. The concept of local timestepping is rooted in the stability of explicit schemes. These schemes become unstable for CFL numbers larger than one. Although this condition is not binding for implicit schemes, the CFL number is nonetheless commonly used to select a local timestep. Specifically, it is common to derive a local timestep from the global timestep such that ∆tn+1 = τ n+1 · i  speedni ∆xi  (3.17)  where speedni is the largest characteristic speed in the control volume, ∆xi a characteristic size, and τ n+1 is calculated using Equation 3.15. Alternatively, a global timestepping scheme can be obtained by selecting the smallest timestep, computed using Equation 3.17, amongst all control volumes. In the present work we propose an alternative means of computing a local timestep. Since the residual was an effective guide for selecting the global timestep in Equation 3.15, we propose to use the local residual to guide the selection of the local timestep. Specifically, ∆tn+1 = ζin+1 τ n+1 · i   ζin+1  =  speedni ∆xi   max min   R U  n  Ri U   2  n    , 102  , 10−2   (3.18)  2  73  3.4. Preconditioned GMRES where Ri U and  τ n+1  n  is the residual in control volume i in conserved variables,  is evaluated using Equation 3.16. In other words, if the flux in  the control volume was small in the previous timestep we expect the local linearization of the residual to be a good representation and a large timestep can be taken. The intended result is that small timesteps are taken in regions of high non-linearity while larger timesteps are taken elsewhere. The goal is to give the scheme some multigrid-like properties such that error waves can be rapidly propagated away from areas where the solution is changing rapidly. For example, it is expected that this scheme will result in small timesteps near an evolving shock while allowing larger timesteps to be taken elsewhere in the domain.  3.4  Preconditioned GMRES  Solving the sparse linear system in Equation 3.13 exactly is prohibitively expensive in computational time and memory. Therefore, as is common, we solve the system approximately using the GMRES [21] method. Solving a linear system with GMRES only requires a means of computing the product of the left-hand-side matrix with an arbitrary vector. Therefore it is common to avoid forming the exact Jacobian matrix and instead approximate the matrixvector product using Frechet derivative. However, since the convergence properties of GMRES are highly sensitive to the conditioning of the matrix, the system must be preconditioned using an approximate factorization of a simplified Jacobian matrix. This matrix is usually chosen to be the Jacobian of the first-order scheme or an approximation to it. In the present work we consider both this “matrix-free” method and a method which explicitly forms the exact high-order Jacobian.  74  3.4. Preconditioned GMRES  3.4.1  Matrix-Free Method  In the matrix-free method the matrix-vector products are approximated using a Frechet derivative R(U + h · a) − R(U ) ∂R a= h ∂U where h is a small value. This scheme requires one flux evaluation per inner GMRES iteration plus one flux evaluation per outer iteration.  3.4.2  Forming the High-Order Jacobian Matrix  Alternatively, the exact Jacobian of the high-order scheme can be formed explicitly. The explicit Jacobian can then be used as an improved preconditioner and to avoid the use of Frechet derivatives and their associated residual evaluations. 3.4.2.1  Without Limiters  The analytic Jacobian can be represented explicitly as ∂R ∂FluxInt ∂FluxInt ∂Flux ∂RecSol ∂RecCoef ∂PVars ≡ = ∂CVars ∂Flux ∂RecSol ∂RecCoef ∂PVars ∂CVars ∂U where FluxInt is the flux integral, Flux are the numerical fluxes, RecSol are the reconstructed solutions at Gauss points, RecCoef are the reconstruction coefficients, PVars are the control volume averages of the primitive variables, and CVars are the control volume averages of the conserved variables. To compute the Jacobian, the following procedure is used at each timestep: 1.  ∂ PVars ∂ CVars  is computed for each control volume and stored. This is the  standard Jacobian for the change of variables from (ρ ρu ρv E)T to (ρ u v P )T . 2. For each Gauss point, do the following for each of the two adjacent control volumes:  75  3.4. Preconditioned GMRES (a)  ∂ RecSol ∂ PVars  =  ∂ RecSol ∂ RecCoef ∂ RecCoef ∂ PVars  is computed. The  ∂ RecCoef ∂ PVars  term is simply the pseudoinverse of the reconstruction matrix RecSol term is a geoprecomputed in Equation 3.9, while the ∂∂RecCoef metric term that depends on the location of the Gauss quadrature point. Although the entire ∂∂RecSol PVars term is purely geometric and could be precomputed and stored, this would dramatically increase the memory requirements since, unlike ∂ RecCoef , this ∂ PVars term is unique for each Gauss point of every control volume. Since the computational efficiency gains would be small, this increase in memory is not worthwhile. ∂ Flux (b) ∂ RecSol , the Jacobian of the Roe flux, is computed. Due to its complexity, the code required to compute this term was automatically generated using the ADIC [4] automatic differentiation tool. Flux = ∂ Flux ∂ RecSol is computed efficiently (c) The product ∂∂PVars ∂ RecSol ∂ PVars by taking advantage of the sparsity of the reconstruction terms that is due to the lack of coupling between solution variables. ∂ Flux couples all solution variables in the reconstruction stencil. ∂ PVars Flux = ∂ Flux ∂ PVars is computed. Since ∂ Flux (d) The product ∂∂CVars ∂ PVars ∂ CVars ∂ PVars couples all solution variables in the reconstruction stencil, this step is computationally intensive. ∂ FluxInt ∂ Flux (e) ∂∂FluxInt CVars = ∂ Flux ∂ CVars is the contribution to the flux integral due to one side of this Gauss point. This component is computed by using the appropriate Gauss integration weight. The result is added to the total flux Jacobian. The sparse analytic Jacobian is found once for every Newton iteration and used to produce the matrix-vector products needed by the GMRES solver. These products require fewer operations to compute than the flux evaluation needed by the matrix-free method. They are also easy to optimize to fully benefit from caching and vector operations that may be available on the processor. 76  3.4. Preconditioned GMRES 3.4.2.2  With Limiters  For transonic and supersonic flows, a limiter must be used to eliminate overshoots in the solution. We will describe a means of including the effect of our limiter, described in Chapter 2, which has been shown to maintain high-order accuracy in smooth regions. The procedure using another limiter, such as the Venkatakrishnan limiter [23], would be the same. The limiting procedure consists in computing the maximum limiting factor Φ ∈ [0, 1] which eliminates overshoots. The limiter is then applied as:   1 Ψx Ψy Ψx2 Ψxy i Ψy 2 i · · ·         RecCoef =         . .  .  Φ    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x ∂y 1 ∂2U 2 ∂y 2                    Φ .. .     ···    ···   Φ  .. .  U   ···     ···    Φ  .. .    .. .  Φ .. .     ···    .. .  .. .  unlmt (3.19)  where unlmt represents the coefficients found by using unlimited reconstruction and Ψ ≡ Φ − 1. The Jacobian is therefore computed exactly as for the case without limiters except for the the ∂ RecCoef term which is computed as: ∂ PVars ∂RecCoef unlmt ∂RecCoef ∂RecCoef ∂Φ ∂RecCoef = + ∂PVars ∂Φ ∂PVars ∂RecCoef unlmt ∂PVars where RecCoef unlmt are the coefficients of the unlimited reconstruction. To ∂Φ balance computational time and memory use, the ∂ PVars term is computed and stored once for each flow variable in each control volume in the same way PVars was for the unlimited case. The computation of the remaining as ∂∂CVars terms and the assembly of ∂ RecCoef is performed “on the fly” while iterating ∂ PVars over Gauss points: •  ∂ RecCoef ∂Φ  term is easily found based on the linear relationship in  Equation 3.19. 77  3.4. Preconditioned GMRES ∂ RecCoef ∂ RecCoef  •  term is also evident from Equation 3.19.  unlmt  ∂ RecCoef  •  unlmt is the pseudoinverse of the reconstruction matrix pre-  ∂ PVars  computed in Equation 3.9. The limiter value Φ computed using our scheme or that of Venkatakrishnan, requires the control volume average value, the neighboring maximum and minimum control volume average values, and the unlimited reconstructed value of the solution at the Gauss points. During the limiting procedure, the minimum value of Φ from all the Gauss points is selected. For the purposes of the Jacobian, the choice of which Gauss point produced the smallest value of Φ is “frozen”. Similarly, the control volume which contains the minimum ∂Φ or maximum control volume average is also preselected. The ∂ PVars term can therefore be decomposed as:  ∂Φ ∂PVars  =  ∂Φ ∂CVAve ∂Φ ∂NeighMin + ∂CVave ∂PVars ∂NeighMin ∂PVars ∂Φ ∂NeighMax ∂Φ ∂RecVal + + ∂NeighMax ∂PVars ∂RecVal ∂PVars  where CVave, NeighMin, and NeighMax are simply select components of PVar and RecVal is the unlimited reconstructed value of the solution at the ∂Φ Gauss point producing the smallest value of Φ. The components ∂ CVave , ∂Φ ∂Φ ∂Φ , ∂ NeighMax , and ∂ RecVal are obtained by computing the ∂ NeighMin derivative of the limiter function. The component ∂∂RecVal PVars is one row ∂ RecCoef ∂ RecSol unlmt and is found as discussed in the previous of ∂ RecCoef ∂ PVars unlmt section.  3.4.3  Preconditioning  Since the rate of convergence of the GMRES method is strongly dependent on the condition number of the matrix, preconditioning is used to alter the spectrum and hence accelerate the convergence rate of the iterative technique. Left preconditioning is applied by modifying the linear system to be solved  78  3.4. Preconditioned GMRES such that Ax = b becomes M −1 Ax = M −1 b where M −1 is an approximate inverse of the preconditioning matrix M ≈ A. A common approach [3, 19] is to use the flux Jacobian of the first-order scheme for M and to use ILU decomposition to form the approximate inverse. This approach has the advantages of being easier to compute and requiring less memory than using the full-order accurate Jacobian. To form the preconditioning matrix M , a procedure similar to that presented in Subsection 3.4.2 is used except that the terms  ∂F lux ∂RecSol ∂RecSol ∂RecCoef  are eliminated.  This method is used for the matrix-free results presented in the present work. Since the high-order Jacobian already needs to be computed in the matrixexplicit method, its ILU decomposition can easily be used as a preconditioner. The increase in memory use can be partially mitigated by using a lower level of fill; as we will show, our results demonstrate that even with low levels of fill, the matrix-explicit method is much better conditioned than the matrix-free method.  3.4.4  Startup  Since Newton-like quadratic convergence cannot be expected during initial timestepping, computing the full-order Jacobian or using a matrix-free method is unnecessarily expensive at that stage. Instead, the first-order Jacobian is used on the left-hand side and the full-order residual is used on the right-hand side. In the present work, we use timestepping with the firstorder Jacobian for the first 10 and 50 iterations for subsonic and transonic cases respectively. We then switch to using either matrix-free GMRES or full-order matrix-explicit GMRES.  79  3.5. Results  3.5 3.5.1  Results Test Cases  Two two-dimensional inviscid flows will be used to evaluate the performance of the proposed scheme. First, subsonic flow at Mach 0.2 around a multielement airfoil with zero angle of attack is considered. The mesh consists of 4913 vertices and is shown in Figure 3.1. The pressure from the steady-state solution using the second- and fourth-order schemes is shown in Figure 3.2; the lower dissipation of the fourth-order solution is evident from the larger extent of peak high and low pressure areas. Second, transonic flow around a NACA 0012 airfoil at Mach 0.8 and an angle of attack of 1.25 degrees is also considered. To prevent oscillations in the solution the limiter described in Chapter 2 is used. The mesh consists of 4656 vertices and is shown in Figure 3.3. The pressure along the upper and lower airfoil surfaces is plotted in Figure 3.4. The benefit of the high-order scheme is evident from the sharper shock on the lower surface. To study the effects of the various methods and parameters discussed in Sections 3.3 and 3.4 we will use what we propose as the best method as a base scheme and compare it with variations in each aspect of the solution strategy. Therefore, unless otherwise noted, the globalization scheme used is a line search method with a local residual-dependent timestep (c2 = 1.1). The initial timestep is based on a CFL number of 0.1. The startup phase consists of 10 iterations for the subsonic case and 50 iterations for the transonic case. During this phase, the Jacobian of the residual is approximated by that of the first-order scheme for both preconditioning and the matrix-vector products needed in the GMRES method. Preconditioning during the startup phase is always achieved using ILU with zero levels of fill (ILU(0)). During the second phase, the full-order exact Jacobian is formed and used for preconditioning and matrix-explicit GMRES. Preconditioning is achieved using ILU(0). The linear system at each iteration is solved to a relative residual drop of 10−3 .  80  3.5. Results  3.5.2 3.5.2.1  Globalization Strategy Line Search  To study the effects of adding the line search step to each timestep, we compare it to a more standard backward-Euler approach. Specifically, we compare the line search algorithm described in Subsection 3.3.2 using the timestep evolution in Equation 3.16 with the backward-Euler approach of Equation 3.13 using the standard SER timestep control of Equation 3.15. In both cases the local timestep is derived from these values using the local residual as in Equation 3.18. The results are presented in the form of residual convergence plots in Figures 3.5 and 3.6. For the backward-Euler method without line search we use values of c2 of 1.0 and 1.1. Since the line search method is expected to be more robust to larger values of c2 we use values of 1.1 and 1.2. The results indicate that there is little difference in the number of iterations required for convergence for the subsonic flow. However, for transonic flow, the method without line search can fail to converge. Specifically it can enter a cycle where the solution alternates between two states. The line search algorithm provides the necessary robustness to avoid this pitfall. The results also illustrate that the robustness of the line search remains unaffected by the choice of parameter c2 . 3.5.2.2  Local Timestep  Next, we consider three methods of deriving the local timesteps ∆ti from the global timestep τ computed using our modified SER method in Equation 3.16. The simplest method is to use a global timestep. The second method is to compute a local timestep using the local CFL condition for each control volume. Finally, the local residual based method of Equation 3.18 is also considered. The resulting convergence plots are shown in Figures 3.7 and 3.8. In the subsonic case, the use of the local residual based timestep results in more rapid convergence than with either CFL based local timestepping or global timestepping. For transonic flow, this trend is even more pronounced, as the non-residual-based methods fail to converge within 350 iterations.  81  3.5. Results Since it is impossible to study all permutations of parameters these results should only be viewed as a representative sample. With the careful selection of startup CFL number and other hand tuning of the startup procedure, convergence is often possible with the methods that failed to converge in these examples. However the gains in robustness and performance demonstrated from the use of the line search method and of residual based local timestep are representative of our general experience.  3.5.3 3.5.3.1  Solving the Linear System Cost of Jacobian Evaluation  Having determined how to minimize the number of linear systems to solve to converge the nonlinear residual, we now turn our attention to the efficiency of the linear solver itself. We begin by comparing the relative CPU time needed to compute a flux integral, a first-order Jacobian, and a high-order Jacobian. The results are given in Tables 3.1 and 3.2. The large increase in computational time needed for flux and Jacobian evaluations for the third-order scheme relative to the second-order scheme can in large part be attributed to the doubling of the number of Gauss points. Since the fourthorder scheme uses the same number of Gauss points as the third-order scheme, a smaller increase in flux and Jacobian evaluation times is observed. For flux evaluations, this increase is purely due to the cost of the reconstruction procedure. For the Jacobian, the increased cost of matrix-matrix products needed for its assembly is also important. 3.5.3.2  Quality of Preconditioning  Next, we compare the relative effectiveness of the ILU decomposition of the first-order and high-order Jacobians. The average number of inner GMRES iterations needed per Newton iteration to obtain a relative residual drop of 10−3 is shown in Table 3.3. The results also include the case where the full LU decomposition of the first-order Jacobian is used to precondition the high-order matrix-free scheme. Since doing the full LU decomposition is prohibitively expensive in both time and memory, this is included solely to 82  3.5. Results  Figure 3.1: Multi-element airfoil mesh consisting of 4913 vertices Order 1 2 3 4  Flux seconds 0.0116 0.0194 0.0399 0.0568  1st Order Jacobian seconds relative to flux 0.217 18.7 0.217 11.2 0.217 5.4 0.217 3.8  Full Order Jacobian seconds relative to flux 0.346 1.122 1.219  17.8 28.1 21.5  Table 3.1: Relative computational time of flux and Jacobian evaluation without limiter for subsonic case Order 1 2 3 4  Flux seconds 0.0109 0.0356 0.0733 0.1137  1st Order Jacobian seconds relative to flux 0.205 18.8 0.205 5.8 0.205 2.8 0.205 1.8  Full Order Jacobian seconds relative to flux 0.375 1.623 2.241  10.5 22.1 19.7  Table 3.2: Relative computational time of flux and Jacobian evaluation with limiter for transonic case 83  3.5. Results  (a) Second-order scheme  (b) Fourth-order scheme  Figure 3.2: Pressure contours for subsonic flow over a multi-element airfoil  84  3.5. Results  Figure 3.3: NACA 0012 airfoil mesh consisting of 4656 vertices  85  3.5. Results demonstrate the upper bounds on the quality of preconditioning that can be achieved when using the first-order Jacobian. The results indicate that the first-order Jacobian is a reasonably good preconditioner for the second-order scheme if ILU with enough levels of fill is used. However, this preconditioner does a poor job for the third- and fourth-order schemes. Even with a full LU decomposition, the number of GMRES iterations required to solve the fourth-order case remains high. Since the matrix-free method must be used when the high-order Jacobian is not available, this large number of GMRES iterations results in a large number of residual evaluations. The costs of these residual evaluations exceed the relative additional cost of computing the high-order Jacobian. Using the full-order Jacobian for the preconditioner, the convergence properties of the high-order schemes is comparable to that of the second-order scheme. In all cases the transonic test case requires fewer GMRES iterations per Newton iteration indicating that the linear systems are better conditioned. Due to the high nonlinearity of the shock, the adaptive timestep method yields smaller timesteps for the transonic case than the subsonic case. This likely explains why the linear system is better conditioned in the transonic case. However, since these smaller timesteps also result in a larger number of Newton iterations, the total number of GMRES iterations for the transonic case exceeds that of the subsonic case. 3.5.3.3  Overall Memory and Computational Cost  The major components contributing to the memory requirement for both matrix-explicit and matrix-free methods are: • The pseudoinverse of the reconstruction matrix. To avoid solving the least-squares problem at each flux evaluation, these matrices need to be precomputed and stored for each control volume. • The Jacobian. For the matrix-free scheme, this will always be the first-order Jacobian. For the matrix-explicit scheme, the high-order Jacobian will inevitably have more fill and require more storage. • ILU decomposition of the Jacobian. The memory required for this 86  3.5. Results  -1.5 -1  Cp  -0.5 0 0.5 1 2nd-order 4th-order  1.5 0  0.2  0.4  0.6  0.8  1  x/c  Figure 3.4: Surface pressure coefficient for transonic flow over a NACA 0012 airfoil  Order 2 3 4  First-Order Jacobian (MF) ILU(1) ILU(4) LU 85.5 49.7 29.2 81.3 55.0 38.8 152.4 105.4 95.7  High-Order Jacobian (ME) ILU(0) ILU(1) 54.2 33.7 34.0 23.4 32.6 23.0  (a) Subsonic case  Order 2 3 4  First-Order Jacobian (MF) ILU(1) ILU(4) LU 41.9 28.6 26.8 40.9 32.9 29.5 58.9 47.5 44.9  High-Order Jacobian (ME) ILU(0) ILU(1) 25.2 13.3 10.9 7.1 20.6 11.9  (b) Transonic case  Table 3.3: Average number of inner GMRES iterations per Newton iteration  87  3.5. Results  No Linesearch c=1.0 No Linesearch c=1.1 Linesearch c=1.1 Linesearch c=1.2  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  2  4  6  8 10 Iterations  12  14  16  18  (a) Second-order scheme for subsonic case No Linesearch c=1.0 No Linesearch c=1.1 Linesearch c=1.1 Linesearch c=1.2  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  5  10  15 Iterations  20  25  30  (b) Fourth-order scheme for subsonic case  Figure 3.5: Convergence plots comparing the line search method with a backward-Euler method for subsonic case  88  3.5. Results  No Linesearch c=1.0 No Linesearch c=1.1 Linesearch c=1.1 Linesearch c=1.2  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  10  20  30  40  50  60  70  Iterations  (a) Second-order scheme for transonic case No Linesearch c=1.0 No Linesearch c=1.1 Linesearch c=1.1 Linesearch c=1.2  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  50  100  150  200  250  Iterations  (b) Fourth-order scheme for transonic case  Figure 3.6: Convergence plots comparing the line search method with a backward-Euler method for transonic case  89  3.5. Results  Residual Based Local TS CFL Based Local TS Global TS  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  10  20  30 Iterations  40  50  60  (a) Second-order scheme for subsonic case Residual Based Local TS CFL Based Local TS Global TS  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  10  20  30  40 Iterations  50  60  70  80  (b) Fourth-order scheme for subsonic case  Figure 3.7: Convergence plots comparing global timestepping, CFL-based local timestepping and residual-based local timestepping for subsonic case  90  3.5. Results  Residual Based Local TS CFL Based Local TS Global TS  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  50  100  150  200 Iterations  250  300  350  400  (a) Second-order scheme for transonic case Residual Based Local TS CFL Based Local TS Global TS  10000  L2 Norm Residual  100  1  0.01  0.0001  1e-06  0  50  100  150  200 Iterations  250  300  350  400  (b) Fourth-order scheme for transonic case  Figure 3.8: Convergence plots comparing global timestepping, CFL-based local timestepping and residual-based local timestepping for transonic case  91  3.5. Results depends not only on the fill of the Jacobian but also on the additional fill due to the decomposition which increases with n for ILU(n). • Krylov subspace. The maximum number of inner GMRES iterations required to solve a Newton iteration and the number of solution unknowns per control volume determines the memory requirement of the Krylov solver. The breakdown of memory requirement along with the total CPU time is shown in Tables 3.4 and 3.5. The additional memory required by the matrixexplicit scheme is due to the increased fill of the Jacobian and resulting preconditioning matrix. However, this is partially offset by the lower fill ratio of the ILU decomposition and the reduced memory use of the Krylov solver. For the subsonic case, the matrix-explicit method preconditioned with ILU(0) requires less than 50% more memory than the matrix-free method with ILU(4) for second- and fourth-order schemes while for the third-order scheme the difference is about 100% more memory. For the transonic case, the matrix-explicit second-order method requires 31% more memory than the matrix-free method, while for the third- and fourth-order scheme, the difference is closer to 100%.  92  Scheme  2 2 2 2 3 3 3 3 4 4 4 4  ME ME MF MF ME ME MF MF ME ME MF MF  Precon -ditioner ILU(0) ILU(1) ILU(1) ILU(4) ILU(0) ILU(1) ILU(1) ILU(4) ILU(0) ILU(1) ILU(1) ILU(4)  Recon -struction 11.7 11.7 11.7 11.7 88 88 88 88 179.2 179.2 179.2 179.2  Jacobian 299.2 299.2 109 109 585 585 109 109 615.2 615.2 109 109  ILU Matrix 299.2 480.3 144.1 300.4 585 1011.7 144.1 300.4 615.2 1087.6 144.1 300.4  Krylov Subsp. 252 156 488 252 172 124 476 316 164 108 808 572  Total (floats) 862.1 947.2 752.9 673.2 1430 1808.6 817.2 813.5 1573.7 1990 1240.4 1160.7  Total (bytes) 6897 7577 6023 5385 11440 14469 6537 6508 12589 15920 9923 9285  Time (seconds) 10.2 9.5 21.8 15.2 16.7 15.7 35 26.2 19.7 20.9 92.1 66.5  Time (res eval) 526 490 1124 784 419 393 877 657 347 368 1621 1171  Table 3.4: Memory in number of floating point numbers per control volume along with run time in seconds and equivalent residual evaluations for subsonic case  3.5. Results  Order  93  Scheme  2 2 2 2 3 3 3 3 4 4 4 4  ME ME MF MF ME ME MF MF ME ME MF MF  Precon -ditioner ILU(0) ILU(1) ILU(1) ILU(4) ILU(0) ILU(1) ILU(1) ILU(4) ILU(0) ILU(1) ILU(1) ILU(4)  Recon -struction 11.7 11.7 11.7 11.7 88 88 88 88 179.2 179.2 179.2 179.2  Jacobian 297.2 297.2 108.7 108.7 580.2 580.2 108.7 108.7 616.7 616.7 108.7 108.7  ILU Matrix 297.2 470.2 142.6 291.7 580.2 986 142.6 291.7 616.7 1073.4 142.6 291.7  Krylov Subsp. 120 72 216 144 116 52 268 220 100 56 316 256  Total (floats) 726.1 851.1 479 556.1 1364.3 1706.2 607.3 708.4 1512.6 1925.3 746.5 835.6  Total (bytes) 5809 6809 3832 4449 10915 13649 4858 5668 12101 15403 5972 6685  Time (seconds) 27.5 27.1 47.2 40.2 184.5 156.4 179.1 155.9 53.5 54.8 101.6 88.1  Time (res eval) 772 761 1326 1129 2517 2134 2443 2127 471 482 894 775  Table 3.5: Memory use in number of floating point numbers per control volume along with run time in seconds and equivalent residual evaluations for transonic case  3.5. Results  Order  94  3.6. Conclusion Next we consider the overall computational time required to converge to the solution. The results presented in Tables 3.4 and 3.5 present both the CPU time in seconds on a single core of an Intel Core 2 processor as well as a multiple of the cost of a residual evaluation. Despite the fact that no attempt to hand tune the parameters of the globalization strategy was made, the runtime in terms of equivalent residual evaluations is competitive with that reported elsewhere. In particular, Blanco and Zingg report [5] that the Newton-Krylov algorithm applied to a second-order accurate matrixdissipation scheme required the equivalent of 660 residual evaluation to converge for the same transonic case. However due to the differences in the discretization scheme and the mesh used, too much emphasis should not be placed on such direct comparisons. The results indicate that in terms of the memory to computational time trade-off the matrix-explicit method with ILU(0) preconditioning and the matrix-free method with ILU(4) preconditioning represent, in most cases, the best trade-offs. Figure 3.9 shows the residual as a function of time for these two schemes for the second- and fourth-order schemes. On these plots the startup procedure is shown with a line while the main stage is shown with a line and a point marker at each iteration. In all cases the matrix-explicit method outperforms the matrix-free method. The difference is most substantial for the subsonic case where the startup phase represents a smaller portion of the total time and where the conditioning of the linear system is poorer. In this case, the matrix-explicit method provides a 38% and 70% reduction in CPU time for the second- and fourth-order schemes respectively. For the transonic case, the matrix-explicit method provides a 32% and 39% reduction in CPU time for second- and fourth-order schemes respectively.  3.6  Conclusion  We have shown that the most commonly used globalization technique may fail to converge, particularly for difficult cases such as the fourth-order solution of transonic flow. This was successfully remedied by the addition of a line search 95  3.6. Conclusion  MF 2nd-order ILU(4) ME 2nd-order ILU(0) MF 4th-order ILU(4) ME 4th-order ILU(0)  10000  100  L2 Norm Residual  1  0.01  0.0001  1e-06  1e-08  1e-10 0  10  20  30 40 Computational Time (seconds)  50  60  70  (a) Subsonic test case MF 2nd-order ILU(4) ME 2nd-order ILU(0) MF 4th-order ILU(4) ME 4th-order ILU(0)  10000  100  L2 Norm Residual  1  0.01  0.0001  1e-06  1e-08  1e-10 0  10  20  30 40 50 60 Computational Time (seconds)  70  80  90  (b) Transonic test case  Figure 3.9: Residual versus computational time for matrix-explicit method and matrix-free method  96  3.6. Conclusion method and the use of the result of the line search to adapt the timestep for the next iteration. A new method of computing the local-timestep from the local residual was shown to dramatically improve robustness and performance of both the second- and fourth-order schemes relative to the use of a global timestep or a local timestep based on the CFL number. A method for efficiently computing the exact Jacobian of the high-order scheme was devised. The incomplete factorization of this Jacobian was shown to be a more effective preconditioner than the factorization of the first-order Jacobian. The difference is particularly important for the fourthorder scheme. In addition the use of this exact Jacobian makes the use of Frechet derivatives and their associated residual evaluations unnecessary in the GMRES method. These factors lead to a reduction in computational time needed for convergence relative to the matrix-free scheme. The difference is particularly important for the fourth-order scheme where a 70% and 39% reduction of computational time is achieved for the subsonic and transonic cases respectively. These reductions in computational time come at the cost of increased memory consumption of 36% and 81% for subsonic and transonic cases respectively. Although these proposed methods help reduce the computational time of both the second- and fourth-order schemes, they have a greater effect on the latter. This helps reduce the gap between the computational effort needed for the fourth-order scheme and the second-order scheme on the same grid thereby making the high-order schemes more attractive. With the use of the proposed globalization techniques and the matrix-explicit method, the fourth-order scheme is only a factor of 1.9 more expensive than the second-order scheme. When compared to the matrix-free second-order scheme, the matrix-explicit fourth-order scheme requires only a factor of 1.1 to 1.3 more computational time and about three times as much memory. Therefore the fourth-order matrix-explicit method is less expensive in time and memory than the second-order matrix-free scheme applied to a mesh with one additional 4 : 1 refinement step.  97  3.7. Bibliography  3.7  Bibliography  [1] T. J. Barth. Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations. In Lecture Series 1994-05, Rhode-Saint-Gen`ese, Belgium, Mar. 1994. von Karman Institute for Fluid Dynamics. [2] T. J. Barth and P. O. Frederickson. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA paper 90-0013, Jan. 1990. [3] T. J. Barth and S. W. Linton. An unstructured mesh Newton solver for compressible fluid flow and its parallel implementation. AIAA paper 95-0221, Jan. 1995. [4] C. Bischof, L. Roh, and A. Mauer. ADIC — An extensible automatic differentiation tool for ANSI-C. Preprint ANL/MCS-P626-1196, Argonne National Laboratory, 1996. [5] M. Blanco and D. W. Zingg. A Newton-Krylov algorithm with a looselycoupled turbulence model for aerodynamic flows. AIAA Paper 2006-0691. Presented at the 44th AIAA Aerospace Sciences Meeting, Jan. 2006. [6] S. De Rango and D. W. Zingg. Higher-order spatial discretization for turbulent aerodynamic computations. American Institute of Aeronautics and Astronautics Journal, 39(7):1296–1304, July 2001. [7] M. Delanaye and J. A. Essers. Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids. American Institute of Aeronautics and Astronautics Journal, 35(4):631–639, Apr. 1997. [8] J. E. Dennis and R. B. Schnabel. Numerical Methods for Unconstrained Optimization. Prentice-Hall, 1983. [9] G. H. Golub and C. F. V. Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, 1983. 98  3.7. Bibliography [10] W. Gropp, D. Keyes, L. C. McInnes, and M. D. Tidriri. Globalized Newton-Krylov-Schwarz Algorithms and Software for Parallel Implicit CFD. International Journal of High Performance Computing Applications, 14(2):102–136, 2000. [11] H. Jiang and P. A. Forsyth. Robust linear and nonlinear strategies for solution of the transonic Euler equations. Computers and Fluids, 24(7):753 – 770, 1995. [12] C. T. Kelley and D. E. Keyes. Convergence analysis of pseudo-transient continuation. SIAM Journal on Numerical Analysis, 35:508–523, 1998. [13] W. A. Mulder and B. van Leer. Experiments with implicit upwind methods for the Euler equations. Journal of Computational Physics, 59:232–246, June 1985. [14] A. Nejat and C. Ollivier-Gooch. A high-order accurate unstructured GMRES algorithm for inviscid compressible flows. In Proceedings of the Seventeenth AIAA Computational Fluid Dynamics Conference. American Institute of Aeronautics and Astronautics, June 2005. [15] A. Nejat and C. Ollivier-Gooch. Effect of discretization order on preconditioning and convergence of a high-order unstructured NewtonGMRES solver for the Euler equations. Journal of Computational Physics, 227(4):2366–2386, 2008. [16] A. Nejat and C. Ollivier-Gooch. A high-order accurate unstructured finite volume Newton-Krylov algorithm for inviscid compressible flows. Journal of Computational Physics, 227(4):2592–2609, 2008. [17] C. F. Ollivier-Gooch and M. Van Altena. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. Journal of Computational Physics, 181(2):729–752, 2002. [18] R. P. Pawlowski, J. N. Shadid, J. P. Simonis, and H. F. Walker. Globalization techniques for Newton–Krylov methods and applications to the  99  3.7. Bibliography fully coupled solution of the Navier–Stokes equations. SIAM Review, 48(4):700–721, 2006. [19] A. Pueyo and D. W. Zingg. Improvements to a Newton-Krylov solver for aerodynamic flows. Thirty-sixth AIAA Aerospace Sciences Meeting, Jan. 1998. AIAA Paper 98-0619. [20] P. L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43:357–372, 1981. [21] Y. Saad and M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal of Scientific and Statistical Computing, 7(3):856–869, July 1986. [22] R. S. Tuminaro, H. F. Walker, and J. N. Shadid. On backtracking failure in Newton-GMRES methods with a demonstration for the Navier-Stokes equations. Journal of Computational Physics, 180(2):549 – 558, 2002. [23] V. Venkatakrishnan. On the accuracy of limiters and convergence to steady-state solutions. AIAA paper 93-0880, Jan. 1993. [24] D. Zingg, S. De Rango, M. Nemec, and T. Pulliam. Comparison of several spatial discretizations for the Navier-Stokes equations. Journal of Computational Physics, 160:683–704, 2000.  100  Chapter 4  Extension to Viscous Flows† 4.1  Introduction  Unstructured grid finite-volume methods have seen wide acceptance in recent times due to their efficient handling of complex geometries. Though third-order schemes for the advective fluxes have been known [2] for almost 20 years, they have seen little use. The number of details which must be properly implemented, such as the correct treatment of curved boundaries, is likely partially to blame. Furthermore, the application of convergence acceleration techniques is also more complicated than for the second-order schemes. Despite this we have recently had success in demonstrating the accuracy of third- and fourth-order schemes for the Euler equations [15]. Just as importantly, we have had success in devising convergence acceleration techniques for these schemes [11, 13, 14]. The result is that we have successfully demonstrated that third- and fourth-order schemes can be more computationally efficient at solving the Euler equations for a given level of accuracy than second-order schemes. We now aim to extend this work to the Navier-Stokes equations. Though other researchers [3, 5] have previously used third-order discretizations for the advective terms, they used traditional second-order discretizations for the viscous terms. The present work is an extension of work previously carried out [16] demonstrating the fourth-order accurate solution of the advectiondiffusion equation. Although the simulation of compressible turbulent flows for aerodynamics applications is our ultimate goal, the present work is †  A version of this chapter will be submitted for publication. C. Michalak and C. OllivierGooch. A high-order accurate unstructured finite-volume method for the Navier-Stokes equations.  101  4.2. Governing Equations restricted to subsonic laminar compressible flows. In Section 4.3 we demonstrate how k-exact reconstruction including boundary conditions can be efficiently carried out. The evaluation of advective and diffusive fluxes using this reconstruction is addressed Section 4.4. The matrix-explicit Newton-Krylov method used to converge the solution to steady state is presented in 4.5. Finally, in Section 4.6, we present numerical experiments demonstrating the accuracy of the method and its computational efficiency.  4.2  Governing Equations  The two-dimensional compressible Navier-Stokes equations for a control volume Vi can be written in integral form as d dt    ρ  =  G =    ρun    ρuu + P n n x    ρvun + P ny   (E + P )un  E  (4.1)  Ωi  Vi    F  (F − G) dS = 0  U dA +       ρu    U =   ρv     ˛  ˆ           0           τxx n ˆ x + τyx n ˆy         τxy n ˆ x + τyy n ˆy (uτxx + vτxy − qx ) n ˆ x + (uτyx + vτyy − qy ) n ˆy  where U is the solution vector in conservative variables, F and G are the inviscid and viscous flux vectors normal to the control volume boundary, n ˆx and n ˆ y are the unit normal components, and un = uˆ nx + vˆ ny . The strain tensor and heat flux vector are τ  = µ  2ux  uy + vx  vx + uy  2vy  2 − µ 3  ux + vy  0  0  ux + vy  102  4.3. Solution Reconstruction q =  1 µ ∇T (γ − 1) P r  where µ is the viscosity and T is the temperature.  4.3  Solution Reconstruction  The basis of the present work is k-exact least-squares reconstruction [1] which has been demonstrated to achieve high-order accuracy for the solution of the Euler equation [15]. This approach has previously been extended to include the diffusive fluxes of the advection-diffusion equation [16]. Our current approach is similar to this previous work, but includes some additional efficiency considerations, particularly with regards to the handling of boundary condition enforcement. Though the flux is evaluated in the conserved variables (ρ, ρu, ρv, E), as is common we reconstruct the primitive variables (ρ, u, v, P ).  4.3.1  Reconstruction for Interior Control Volumes  In the finite-volume method, the domain is tessellated into non-overlapping control volumes. Although our approach is equally applicable to vertexand cell-centered control volume discretizations, only results using vertexcentered grids will be presented herein. Each control volume Vi has a geometric reference point xi . This reference point is taken as the cell centroid for cell-centered control volumes and the vertex for vertex-centered control volumes. For any smooth function U (x) and its control volume averaged values U i , the k-exact least-squares reconstruction uses a compact stencil in the neighborhood of control volume i to compute an expansion UiR (x−xi ) that conserves the mean in control volume i and reconstructs exactly polynomials of degree ≤ k (equivalently, UiR (x − xi ) − U (x) = O ∆xk+1 ) . Conservation of the mean requires that the average of the reconstructed function Ri and the original function U over control volume i be the same: 1 Vi  ˆ UiR (x − xi ) dA = Vi  1 Vi  ˆ U (x) dA ≡ U i .  (4.2)  Vi  103  4.3. Solution Reconstruction The expansion UiR (x − xi ) can be written as: UiR (x − xi ) =  U |xi + +  ∂2U ∂x2  ∂2U + ∂y 2  ∂U ∂x  xi  xi  (x − xi ) + xi  ∂U ∂y  (x − xi )2 ∂2U + 2 ∂x∂y  (y − yi ) xi  ((x − xi )(y − yi )) (4.3) xi  (y − yi )2 + ··· 2  Taking the control volume average of this expansion over control volume i and equating it to the mean value gives Ui =  U |xi + +  ∂2U ∂x2  where xn y m  i  ∂U ∂x  xi  xi + xi  ∂U ∂y  ∂2U x2 + 2 ∂x∂y  1 ≡ Ai  yi  (4.4)  xi  xy + xi  ∂2U ∂y 2  xi  y2 + ··· 2  ˆ (x − xi )n (y − yi )m dA.  (4.5)  Vi  are control volume moments. kth-order accuracy requires that we compute the (k − 1)th derivatives by minimizing the error in predicting the mean value of the reconstructed function for control volumes in the stencil of neighboring control volumes {Vj }i . In other words, we aim to minimizing the difference between the control volume average U j and the average of the reconstruction Ri over control volume j. The minimum number of control volumes in the reconstruction stencil is equal to the number of derivative terms to be approximated. In practice, to improve conditioning, we increase this minimum to 4, 10, and 16 for second-, third-, and fourth-order schemes respectively. The stencil is constructed by first considering the first-neighbors of the control volume. If additional neighbors are needed they are added in order of topological proximity. All neighbors at given topological distance are added at once. The resulting stencil is therefore topologically centered but may exceed the minimum specified size. The mean value, for a single control volume Vj , of 104  4.3. Solution Reconstruction the reconstructed function UiR is ˆ  1 Aj  ∗−4ex Vj  ∂U + ∂x ∂U + ∂y ∂2U + ∂x2  xi  1 Aj  xi  1 Aj  xi  ∂2U ∂x∂y  +  ∂2U ∂y 2  +  UiR (x − xi )dA = U |xi  xi  xi  (4.6)  ˆ (x − xi )dA Vj  ˆ  (y − yi )dA Vj  1 2Aj 1 Aj 1 2Aj  ˆ (x − xi )2 dA  ˆ  Vj  (x − xi )(y − yi )dA Vj  ˆ  (y − yi )2 dA + · · · Vj  To avoid computing moments of each control volume in {Vj }i about xi , replace x − xi and y − yi with (x − xj ) + (xj − xi ) and (y − yj ) + (yj − yi ), respectively. Expanding and integrating, we obtain 1 Aj  ˆ ∗ −4ex Vj  + + + +  UiR (x − xi ) = U |xi +  ∂U ∂y  (xj + (xj − xi )) xi  y j + (yj − yi ) xi  ∂2U ∂x2  xi  x2 j + 2xj (xj − xi ) + (xj − xi )2 2  ∂2U ∂x∂y ∂2U ∂y 2  ∂U ∂x  xy j + xj (yj − yi ) + (xj − xi )y j +(xj − xi )(yj − yi )) xi  xi  y 2 j + 2y j (yj − yi ) + (yj − yi )2 + ··· 2  (4.7)  This equation is written for every control volume in the stencil. These equations combined with the mean constraint form a constrained leastsquares system of the form:  105  4.3. Solution Reconstruction     x ¯  y¯  x2  xy  y2  ···  wi1 xi1  wi1 yi1  wi1 x2 i1  wi1 xy i1  wi1 y 2 i1  ···  1    w  i1   w  i2   wi3   .  ..   wiN    wi2 xi2  wi2 yi2  wi2  x2 x2  wi2 xy i2  i2  wi2  y2  i2  y2  i3  wi3 xi3 .. .  wi3 yi3 .. .  wi3  wi3 xy i3 .. .  wi3  .. .  wiN xiN  wiN yiN  wiN x2 iN  wiN xy iN  wiN y 2 iN  i3  .. .  ··· ··· .. . ···  U    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2                                   .. .  i    Ui    w U  1 i1   w U 2 i2  =  wi3 U 3   ..  .   wiN U N               (4.8)  where the first row is the mean constraint, and the geometric terms are xn y m  ij  ≡  1 Aj m  =  ˆ ((x − xj ) + (xj − xi ))n ((y − yj ) + (yj − yi ))m dA Vj n  m! n! (xj − xi )k (yj − yi )l xn−k y m−l j l! (m − l)! k! (n − k)! l=0 k=0  and w, the weights, can be chosen to emphasize geometrically nearby data wij =  1 . |xj − xi |n  (4.9)  where n is typically chosen to be 0, 1, or 2. Although the use of n = 0 has been shown to lead to some accuracy issues for highly anisotropic grids [9], in our experience it is the best choice for successful convergence of the scheme and is therefore used in the present work. Since the matrix in Equation 4.8 depends purely on geometric terms it is possible to reduce the reconstruction step at each flux evaluation to a matrix-vector multiplication by precomputing a pseudoinverse during a 106  4.3. Solution Reconstruction preprocessing stage. To obtain this form, we first analytically eliminate the mean constraint from Equation 4.8. Next, we compute the pseudoinverse of the resulting unconstrained least-squares problem. For this we use singularvalue decomposition (see, for instance, [6] for details). The mean constraint can then be reintroduced to yield a matrix, which we store for each control volume, that can be used to compute the reconstruction coefficients using a single matrix-vector multiplication operation of the form                  U      ∂U   ∂x  ∂U   ∂y 2 1∂ U   2 ∂x2  ∂2U  ∂x∂y   1 ∂2U  2 ∂y 2  .   ..  =  Ui    U  1   U 2  A  U3   ..  .   UN         .       (4.10)  i  Since the matrix A is based on purely geometric terms, the same pseudoinverse matrix can be used to reconstruct all the flow variables in the control volume which helps minimize memory use. Furthermore, this reconstruction form helps in constructing the exact flux Jacobian which was shown in Chapter 3to be helpful in obtaining rapid convergence of high-order accurate schemes.  4.3.2  Boundary Condition Enforcement  In general, boundary conditions can be enforced either through special flux functions or by modifying the reconstruction step. In practice, we use special flux functions for boundary conditions of the inviscid flux and for maintaining the adiabatic condition at walls. We modify the reconstruction step to maintain the no-slip condition at the wall by enforcing Dirichlet boundary conditions for the components of velocity u and v. However, our approach to boundary condition enforcement is general enough to be applied to any Dirichlet or Neumann boundary condition. In previous work with the advection-diffusion equation [16] boundary conditions were enforced by introducing additional constraints on the least-squares problem in Equation 107  4.3. Solution Reconstruction 4.8. This resulted in the boundary conditions being satisfied exactly at the Gauss integration points. In the present work we modify this approach by instead augmenting the least-squares problem with additional reconstruction goals. While this does not ensure that the boundary conditions are met exactly, it does guarantee that they are satisfied within the truncation error of the scheme and it simplifies the reconstruction process. Suppose that the solution must satisfy a Dirichlet boundary condition such that along part of the boundary ∂Ω1 the reconstruction approximately satisfies R(x − xi ) = f1 (x). Rather than attempting to satisfy this condition everywhere along the boundary, we consider only the solution at each Gauss integration point xg on the boundary. In terms of the reconstruction polynomial coefficients, this reconstruction goal can be expressed as f1 (xg ) =  U |xi + +  ∂2U ∂x2  ∂2U + ∂y 2  ∂U ∂x  xi  xi  (xg − xi ) + xi  ∂U ∂y  (xg − xi )2 ∂2U + 2 ∂x∂y  (yg − yi ) xi  (xg − xi )(yg − yi ) xi  (yg − yi )2 + ··· 2  (4.11)  The boundary condition can be incorporated into Equation 4.8 by adding a row to the matrix for each boundary Gauss point giving                    1  x ¯  y¯  x2  xy  y2  ···  wia  wia δxia wia δyia  wia δx2ia  wia δxia δyia  2 wia δyia  ···  wib  wib δxib  wib δx2ib  wib δxib δyib  2 wib δyib  ···  wi1 wi2 wi3 .. . wiN  wi1 xi1 wi2 xi2 wi3 xi3 .. . wiN xiN  wib δyib wi1 yi1 wi2 yi2 wi3 yi3 .. . wiN yiN  wi1  x2  wi2  x2  wi3  x2  i1 i2 i3  .. . wiN  x2  iN  wi1 xy i1 wi2 xy i2 wi3 xy i3 .. . wiN xy iN  wi1  y2  i1  ···  wi2  y2  i2  ···  wi3  y2  i3  .. . wiN  y2  iN  ··· .. .                       U    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2                  ···  108  .. .  i  4.3. Solution Reconstruction     Ui    wia f1 (xa )    wib f1 (xb )    wi1 U 1 =  w U  2 i2   w U 3 i3   ..  .                     (4.12)  wiN U N where the weights for the Gauss point rows can be computed as wig = cb  1 |xg − xi |n  (4.13)  and δxig ≡ xg − xi δyig ≡ yg − yi where n is the same value used in Equation 4.9 and cb determines how strongly the boundary condition is enforced. Except where otherwise noted, we use cb = 6. Neumann boundary conditions can be applied in an analogous way. If the reconstruction polynomial is to approximate the boundary condition  ∂UiR (x−xi ) ∂n  = f2 (x) we simply add reconstruction goals at the  boundary Gauss points of the form f2 (xg ) = ∇UiR (xg − xi ) · n ˆ   ∂U ∂2U = nx  + ∂x ∂x2  xi  ∂2U (xg − xi ) + ∂x∂y    ∂U ∂2U ny  + ∂y ∂x∂y  xi  ∂2U (xg − xi ) + ∂y 2    (yg − yi ) + · · · + xi    (yg − yi ) + · · · xi  where n ˆ is the unit boundary normal. As for the interior case, the pseudoinverse of the constrained least-squares 109  4.3. Solution Reconstruction problem of Equation 4.12 can be computed to reduce the reconstruction to a matrix-vector product of the form   U                    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2                  .. .    =  f1 (xa )    f1 (xb )    Ui    U1   U  2   U3   ..  .   Ab Ai  i           ,          UN  where the component Ab of the pseudoinverse has as many columns as there are boundary Gauss points. Since the value of f1 (x) does not change at each iteration, the boundary reconstruction step can be further optimized by precomputing a vector v b = Ab  f1 (xa ) f1 (xb )  .  The reconstruction can then be found by    U                    ∂U ∂x ∂U ∂y 1 ∂2U 2 ∂x2 ∂2U ∂x∂y 1 ∂2U 2 ∂y 2                  .. .    Ui    U  1   U 2  = Ai   U3   ..  .   UN          + vb .       (4.14)  i  Therefore the computational cost of reconstruction at the boundaries is almost identical to that of the interior. However since the different flow variables do not need to satisfy the same boundary conditions, the pseudoinverse matrix  110  4.4. Flux Integration Ai and vector vb will need to be stored independently for each flow variable. Since only a small fraction of control volumes are adjacent to the boundary, this is not a major issue. Boundary conditions which linearly couple flow variables or their derivatives can also be implemented. This is most easily accomplished by considering the constrained least-squares problem which simultaneously solves the reconstruction polynomial for all flow variables at a control volume. However this causes the memory and computational cost of boundary control volume reconstruction to increase quadratically with the number of flow variables. Therefore we elect to only simultaneously solve the reconstruction for the flow variables which are coupled by boundary conditions and independently solve those that are not.  4.4  Flux Integration  Once the reconstruction polynomial coefficients are computed as in Section 4.3, the solution and its gradients can easily be obtained at any Gauss integration point along the control volume boundary. In the present work, we use a vertex-centered centroidal-dual control volume which is obtained by connecting mesh face centroids with cell centroids. This results in two control volume edges per face in the primal mesh. Gauss integration is carried out on these dual edges using one Gauss point per control volume edge for the second-order scheme and two Gauss points for the third-, and fourth-order schemes. For the third-, and fourth-order schemes, it is well known that Gauss points must be located on the curved boundary rather than on the straight edge connecting vertices. Recently, it has been shown that for vertex-centered schemes the boundary face centroid used for the construction of the dual mesh must also be correctly located on the boundary, even for second-order schemes [4]. For this reason, we use curved boundary information for determining Gauss point locations, weights, and normals for all orders of accuracy. The inviscid component of the flux is computed by evaluating the solution at the Gauss point using the reconstruction from both incident control 111  4.5. Solution Convergence volumes as input to Roe’s approximate Riemann solver [19]. The viscous flux requires not only the solution, but also the solution gradients. For this we use the average gradient obtained from the reconstruction of the two incident control volumes. For a p-order accurate reconstruction scheme, fluxes that depend only on the solution will be p-order accurate. Since the solution gradient is only (p − 1)-order accurate, the accuracy of the viscous flux is an order of accuracy lower than that of the inviscid flux. However, previous work [16] has analytically demonstrated that the even-order discretization of the Laplace equation on equilateral triangular grids benefits from error cancellation leading to an overall scheme that is p-order accurate. Numerical results in the same work confirmed the p-order accuracy of the overall scheme for the advection-diffusion equation on an isotropic unstructured grid. At the onset of this work it was not known whether this could also be expected of the Navier-Stokes equations on an unstructured grid.  4.5  Solution Convergence  The solution is converged to steady state using a Newton-GMRES [20] algorithm. The exact full-order Jacobian is used for forming the preconditioning matrix and for performing the matrix-vector products needed by the GMRES algorithm. The preconditioner used is ILU with one level of fill for the second-order scheme and zero levels of fill for the third- and fourth-order schemes. As in our work with the Euler equation, reproduced as Chapter 3 in this thesis, the globalization technique used is a combination of a local timestepping scheme with a line search algorithm. The local timestep is a function of a global timestep and the relative size of the local residual. The global timestep is increased as the solution converges using the successive evolution-relaxation (SER) [12] scheme. Convergence is facilitated by scaling the local timestep by the size of the local residual from the previous iteration. The line search algorithm is used to detect when an exceedingly large timestep is taken which may adversely affect the nonlinear convergence of the scheme. When this occurs, the subsequent timestep size is reduced. In the manuscript presented in Chapter 3, it was shown that the exact 112  4.6. Results high-order accurate flux Jacobian for the Euler equation could be found using the chain rule ∂FluxInt ∂FluxInt ∂Flux ∂RecSol ∂RecCoef ∂PVars = ∂CVars ∂Flux ∂RecSol ∂RecCoef ∂PVars ∂CVars where FluxInt is the flux integral, Flux are the numerical fluxes, RecSol are the reconstructed solutions at Gauss points, RecCoef are the reconstruction coefficients, PVars are the control volume averages of the primitive variables, and CVars are the control volume averages of the conserved variables. In that work, a method to efficiently compute these terms was also presented. For viscous equations the flux depends not only on the solution but also its gradient. To accommodate this, we simply need to replace RecSol, the vector of reconstructed solutions at Gauss points, with a larger vector RecSolGrad which contains the solution and its two gradient terms. The term ∂ RecSolGrad ∂ RecCoef becomes the sensitivity of the reconstructed solution and solution gradients ∂ Flux to the reconstruction coefficients, and the term ∂ RecSolGrad becomes the Jacobian of the inviscid and viscous fluxes with respect to the reconstructed solution and its gradients. With these changes, the method to efficiently assemble the Euler flux Jacobian described in Chapter 3 can readily be used to assemble the Navier-Stokes flux Jacobian. The sparsity pattern of the assembled Jacobian matrix does not change through the addition of the viscous terms.  4.6 4.6.1  Results Cylindrical Couette Flow  To numerically evaluate the order of accuracy of the scheme we consider flow around concentric cylinders. The outer cylinder with a radius of 2 is fixed and adiabatic while the inner cylinder has a radius of 1 and is rotated clockwise with a tangential Mach number of 0.5. Since isothermal wall conditions for the inner cylinder would be difficult to enforce with the current choice of reconstruction variables, we obtain an equivalent result by augmenting the  113  4.6. Results Reconstruction Order 2 3 4  L1 1.60 2.22 3.54  L2 1.59 2.25 3.57  L∞ 1.37 2.38 3.26  Table 4.1: Order of error norms in velocity computed using a regression fit through data from all grids fixed velocity wall condition with an isobaric condition. To avoid Taylor instabilities a Reynolds number of 10 based on the inner cylinder radius is used. Viscosity is considered constant for this case which allows comparison with the analytic exact solution [8]. In particular, in this case the velocity profile reduces to that of the incompressible flow with ro  vθ = Ωi ri rro ri  − −  r ro ri ro  where vθ is the tangential flow velocity, Ωi is the angular velocity of the inner wall, r is the radius, and ri and ro and the radius of the inner and outer cylinders. The numerical solution is carried out on a series of four meshes consisting of 219, 840, 3256, and 12812 vertices shown in Figure 4.1. Since the grids are unstructured and not aligned with the expected flow gradients and since the discretization is carried out in Cartesian coordinates, the error in the numerical solution should not have any unexpected benefits from the simplicity of the solution expressed in radial coordinates. The L2 norms of error in the velocity terms for the second-, third-, and fourth-order schemes is shown in Figure 4.2 and the convergence order of the error norms is shown in Table 4.1. The results indicate that the effective order of accuracy of the schemes falls between the nominal order of the flux calculation of the solution and gradient terms. To demonstrate that the enforcement of boundary conditions using a least-squares approach is not detrimental to the overall accuracy of the scheme, we consider the effects of varying cb in Equation 4.13. The L2  114  4.6. Results  (a) 219 vertex mesh  (b) 840 vertex mesh  (c) 3256 vertex mesh  (d) 12812 vertex mesh  Figure 4.1: Four coarsest grids used in cylindrical Couette flow case  115  4.6. Results  0.01  Second-Order Third-Order Fourth-Order  L2 Norm Error  0.001  0.0001  1e-05  1e-06 200  400  800  1600 3200 Num CVs  6400 12800  Figure 4.2: Convergence of error in velocity for cylindrical Couette flow norm of error in the velocity terms sampled at 10, 000 points along the inner cylinder boundary along with the norms of error in the interior of the domain is plotted in Figure 4.3. The results indicate that increasing the value of cb only serves to decrease the error in the constrained terms at the boundary without affecting the overall accuracy of the solution in the interior.  4.6.2  Flow Over an Airfoil  Next we consider subsonic laminar compressible flow around an airfoil. Specifically, flow around a NACA 0012 airfoil at zero angle of attack at a Mach number of 0.5 and Reynolds number of 5000 is simulated on a series of grids. These grids were generated by adaptive refinement using the second-order scheme [17]. A total of six grids is considered with the coarsest containing 1345 vertices and the finest containing 37399 vertices. The grids are unstructured, anisotropic, and are not self-similar. The coarsest grid is shown in Figure 4.4 and a sample flow solution is given in Figure 4.5. The airfoil surface is considered adiabatic and Sutherland’s law is used to compute the viscosity.  116  4.6. Results  0.01  Overall cb = 1 Bdry cb = 1 Overall cb = 6 Bdry cb = 6 Overall cb = 100 Bdry cb = 100  0.001  L2 Norm Error  0.0001 1e-05 1e-06 1e-07 1e-08 1e-09 200  400  800  1600  3200  6400 12800  Num CVs  Figure 4.3: Convergence of error in velocity for the fourth-order scheme at the boundary and in the interior using different least-squares weights for the boundary condition terms The convergence process is facilitated by using the converged solution from the previous coarser mesh as the initial solution on the next mesh. The solution is transferred by Gauss integration on the fine grid of the values obtained from the solution reconstruction on the coarse grid. Using this grid-sequencing method, the CPU time needed for a solution scales nearly linearly with the number of control volumes. The time required for a solution at each grid level is shown in Figure 4.6. These results along with the cumulative time required for the grid sequencing are shown in Table 4.2. In general, the third- and fourth-order schemes are a factor of 2.5 and 2.8, respectively, more costly than the second-order scheme. The majority of this additional cost can be attributed to the doubling in the number of Gauss points required in the high-order schemes. From the converged flow solution, the pressure and viscous drag coefficients are computed and the flow separation point is found. Table 4.3 demonstrates that the results are in general agreement with other published values. Next we consider the convergence of viscous drag with mesh refine117  4.6. Results  (a) 1345 vertex grid  (b) 2681 vertex grid  (c) 5363 vertex grid  Figure 4.4: Three coarsest grids used for airfoil case  118  4.6. Results  Figure 4.5: Mach number contours as computed on the 10541 vertex mesh using the fourth-order scheme  CPU Time (s)  1000  Second-Order Third-Order Fourth-Order  100  10  1 1000  10000  100000  Num CVs  Figure 4.6: Computational time required for convergence at each grid level  119  4.6. Results  Number of Control Volumes 1345 2681 5363 10541 20016 37399  Incremental 2nd 3rd 1.6 4.8 2.8 8.1 7.6 17.3 13.9 31.9 32.4 83.6 94.2 190.0  Time 4th 5.3 8.8 18.8 37.1 82.3 227.4  Cumulative Time 2nd 3rd 4th 1.6 4.8 5.3 4.4 12.9 14.1 10.4 25.4 27.6 21.5 49.2 55.9 46.3 115.5 119.4 126.6 273.6 309.7  Table 4.2: Computational time in seconds required for convergence of airfoil test case using grid sequencing  Method  Grid Size  Second-Order Third-Order Fourth-Order ARC2D [7] Mavriplis [10] Radespiel [18]  37399 37399 37399 320 x 128 320 x 64 512 x128  Viscous Drag Coefficient 0.0311 0.0323 0.0321 0.0321 0.0332 0.0330  Pressure Drag Coefficient 0.0237 0.0224 0.0225 0.0221 0.0229 0.0224  Separation Point 0.768 0.778 0.789 0.824 0.814 0.814  Table 4.3: Comparison of fine grid results for airfoil case  120  4.6. Results  0.038  Second-Order Third-Order Fourth-Order  Viscous Drag Coefficient  0.037 0.036 0.035 0.034 0.033 0.032 0.031 0  10000  20000 Num CVs  30000  40000  Figure 4.7: Viscous drag coefficient grid convergence ment in Figure 4.7 and of pressure drag in Figure 4.8. We note that the fourth-order scheme exhibits excellent grid convergence of drag coefficients. In fact, both drag coefficients are within 4 counts of the converged values beginning at the third-coarsest mesh (5363 vertices) and within 2 counts of the converged values beginning at the fourth-coarsest mesh (10541 vertices). The third-order scheme also exhibits good viscous drag convergence but the pressure drag varies by 8 counts between the second-finest and finest grids. Like the third-order scheme, the second-order scheme does not exhibit grid-convergence of drag. Furthermore, the fine grid solution differs significantly from that of the third- and fourth-order schemes. Next we consider the separation point in Figure 4.9. The fourth-order scheme exhibits convergence of the separation point to within 1% of chord length beginning at the fourth-coarsest mesh while the third- and second-order schemes continue to see fluctuations of up to 5%. Overall, the grid convergence study demonstrates that the fourth-order scheme is remarkably superior to the second- and third-order schemes. The fourth-order solution on the third-coarsest mesh is arguably superior to the second-order solution on the finest mesh and requires only  1 5  the computa121  4.6. Results  0.029  Second-Order Third-Order Fourth-Order  Pressure Drag Coefficient  0.028 0.027 0.026 0.025 0.024 0.023 0.022 0  10000  20000 Num CVs  30000  40000  Figure 4.8: Pressure drag coefficient grid convergence  1  Second-Order Third-Order Fourth-Order  Separation Point (x/c)  0.95 0.9 0.85 0.8 0.75 0.7 0  10000  20000  30000  40000  Num CVs  Figure 4.9: Grid convergence of the separation point on the upper airfoil surface  122  4.7. Conclusion tional time.  4.7  Conclusion  The scheme based on a fourth-order reconstruction was shown to obtain an order of accuracy of approximately 3.5 under a mesh refinement study of the cylindrical Couette flow. The subsonic laminar airfoil test case demonstrated that this accuracy translates to rapid grid convergence of lift, drag and separation points relative to the second-order scheme. When combined with a matrix-explicit Newton-GMRES convergence acceleration technique, the fourth-order scheme was shown to be only a factor of 2.8 more computationally expensive than a second-order scheme on the same grid. The results indicate that the fourth-order scheme can be less computationally expensive for obtaining a level of accuracy typically desired for aerodynamics simulations. Future work will examine the extension of the present scheme to turbulent flows. One significant challenge will lay in the stability of the high-order discretization of the turbulence model. Furthermore the need for curved boundary treatment for the high-order scheme combined with the highlyanisotropic grids required for efficient discretizations of turbulent flows can lead to malformed control volumes unless the interior faces are also curved.  123  4.8. Bibliography  4.8  Bibliography  [1] T. J. Barth. Recent developments in high order k-exact reconstruction on unstructured meshes. AIAA paper 93-0668, Jan. 1993. [2] T. J. Barth and P. O. Frederickson. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. AIAA paper 90-0013, Jan. 1990. [3] M. Delanaye and J. A. Essers. Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids. American Institute of Aeronautics and Astronautics Journal, 35(4):631–639, Apr. 1997. [4] B. Diskin and J. Thomas. Accuracy analysis for mixed-element finitevolume discretization schemes. Technical Report 2007-08, National Institute of Aerospace, 2007. [5] P. Geuzaine, M. Delanaye, and J.-A. Essers. Computation of high Reynolds number flows with an implicit quadratic reconstruction scheme on unstructured grids. In Proceedings of the Thirteenth AIAA Computational Fluid Dynamics Conference, pages 610–619. American Institute of Aeronautics and Astronautics, 1997. [6] G. H. Golub and C. F. V. Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, 1983. [7] C. F. Gooch. Solution of the Navier-Stokes Equations on Locally-Refined Cartesian Meshes Using State-Vector Splitting. PhD thesis, Stanford University, 1993. [8] C. Illingworth. Some solutions of the equations of flow of a viscous compressible fluid. Proceedings of the Cambridge Philosophical Society, 46:469–478, 1950. [9] D. J. Mavriplis. Revisiting the least-squares procedure for gradient reconstruction on unstructured meshes. In Proceedings of the 16th AIAA Computational Fluid Dynamics Conference, 2003. 124  4.8. Bibliography [10] D. J. Mavriplis, A. Jameson, and L. Martinelli. Multigrid solution of the Navier-Stokes equations on triangular meshes. ICASE Report No. 89-11, NASA Langley Research Center, 1989. [11] C. Michalak and C. Ollivier-Gooch. Matrix-explicit GMRES for a higher-order accurate inviscid compressible flow solver. In Proceedings of the Eighteenth AIAA Computational Fluid Dynamics Conference, 2007. [12] W. A. Mulder and B. van Leer. Experiments with implicit upwind methods for the Euler equations. Journal of Computational Physics, 59:232–246, June 1985. [13] A. Nejat and C. Ollivier-Gooch. Effect of discretization order on preconditioning and convergence of a high-order unstructured NewtonGMRES solver for the Euler equations. Journal of Computational Physics, 227(4):2366–2386, 2008. [14] A. Nejat and C. Ollivier-Gooch. A high-order accurate unstructured finite volume Newton-Krylov algorithm for inviscid compressible flows. Journal of Computational Physics, 227(4):2592–2609, 2008. [15] C. Ollivier-Gooch, A. Nejat, and C. Michalak. On obtaining high-order finite-volume solutions to the Euler equations on unstructured meshes. In Proceedings of the Eighteenth AIAA Computational Fluid Dynamics Conference, 2007. AIAA-2007-4464. [16] C. F. Ollivier-Gooch and M. Van Altena. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. Journal of Computational Physics, 181(2):729–752, 2002. [17] D. Pagnutti and C. Ollivier-Gooch. A generalized framework for high order anisotropic mesh adaptation. Computers and Structures, To appear, 2009. [18] R. Radespiel. A cell-vertex multigrid method for the Navier-Stokes equations. NASA TM-101557, 1989.  125  4.8. Bibliography [19] P. L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43:357–372, 1981. [20] Y. Saad and M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal of Scientific and Statistical Computing, 7(3):856–869, July 1986.  126  Chapter 5  Conclusion In the previous chapters a number of specific improvements to the state of the art of the unstructured high-order finite-volume method were presented. In Chapter 2 a means of eliminating oscillations of the numerical solution near shocks while maintaining the high-order accuracy of the scheme and computational efficiency of least-squares reconstruction was devised. The explicit formation of the high-order Jacobian combined with novel globalization techniques was shown in Chapter 3 to provide significant improvement to the Newton-Krylov convergence strategy. In Chapter 4 a method to discretize diffusive fluxes of the Navier-Stokes equations and a means of enforcing the associated boundary conditions was presented and shown to attain the expected order of accuracy. In this chapter the combined significance of these works is discussed, how these results affect the relative competitiveness of the high-order finite-volume method compared with other high-order methods is addressed, and future developments are suggested.  5.1  Significance of Results  The intended combined significance of the manuscript chapters was to resolve specific limitations of the high-order least-squares reconstruction based finitevolume scheme so as to pave the way for future research aimed at exploring the feasibility and competitiveness of this method for aerodynamics simulations. The common theme of the work is the demonstration that for the simplified problem of two-dimensional inviscid and laminar flow the method can be accurate and efficient. While this does not prove that the method will be efficient or even accurate for the more difficult case of practical interest, three-dimensional turbulent flow, it provides a critical stepping stone in that 127  5.1. Significance of Results direction.  5.1.1  Limiter  The manuscript in Chapter 2 introduces a slope limiter for the least-squares reconstruction which does not inhibit high-order accuracy in smooth regions while maintaining monotonicity at shocks. Just as importantly, this is achieved in a computationally efficient manner and without degrading the steady-state convergence properties of the scheme. Previous researchers [3, 7, 17] have used a nominally high-order least-squares reconstruction together with a limiter to resolve flows with shocks. However none of these works have demonstrated or made explicit claims that the high-order accuracy of the method was maintained in smooth regions when such a limiter is used. The broad family of essentially non-oscillatory (ENO) scheme [1, 10, 20, for instance] including the Weighted ENO (WENO) schemes [9, 11, 13, for instance] provides a viable solution to the accuracy and monotonicity challenge. However these schemes do not converge to steady state as efficiently as least-squares reconstruction schemes. Furthermore, the reconstruction step is much more costly than that of the least-squares reconstruction. A hybrid approach, named Quasi-ENO [19], has also been developed. In this approach, the reconstruction is carried out in a least-squares sense, but the weights assigned to the reconstruction goals are modified at each iteration to emphasize smooth data. Although found to be effective, this method is computationally expensive since the pseudoinverse of the least-squares matrix cannot be precomputed and stored. The cost savings enjoyed from this precomputation by regular least-squares reconstruction schemes increases dramatically as the order of accuracy, and therefore the size of the leastsquares matrices, increases. For these reasons the limiter approach present in this work is more attractive than ENO-type schemes for steady-state solutions with relatively simple shock formations typically seen in transonic aerodynamics. The possible application of the limiter to unsteady problems with more complex shock formations for which ENO schemes have primarily been developed has yet to be studied.  128  5.1. Significance of Results The significance of the development of the new limiter may extend beyond high-order accurate schemes. The results in Chapter 2 demonstrated that even for second-order schemes, the new limiter provided significantly better accuracy than what is considered the de facto standard limiter [22] for unstructured finite-volume methods. The degradation of the accuracy of the solution with the use of limiters is a well-known problem. In fact, it is common practice to avoid using a limiter for cases where only weak shocks are expected. Such shocks occur, for example, in flow over a wing designed to operate efficiently in transonic regimes. The choice to avoid the use of limiters is made due to the improved convergence properties and accuracy in smooth regions of the unlimited scheme. However, this choice comes at the cost of introducing oscillations at the shock. The new limiter applied to the second-order scheme was shown to maintain the same accuracy as the unlimited scheme, therefore making the trade-off in choosing a good quality solution away from the shock at the expense of the quality of the solution at the shock unnecessary. Since the new limiter fits in essentially the same framework as other limiters, existing finite-volume codes which incorporate Venkatakrishnan’s [22] or Barth-Jespersen’s [3] limiter can easily be extended to accommodate it. The extension of the limiter to turbulent viscous flows is expected to be straightforward. The only significant difference for the viscous case has to do with the no-slip wall boundary condition. The method of maintaining accuracy at stagnation points, which was developed and tested on the inviscid case, is expected to be effective at allowing the limiter to maintain accuracy near no-slip walls. In fact, the application of the limiter to viscous flows should be simpler than for inviscid flows since the special treatment of tangential flow at walls will no longer be necessary. Likewise, the extension to three-dimensional flows is not expected to be problematic. The two-dimensional assumption is only used in the development of the special treatment of tangential flow at the wall. Therefore, for inviscid three-dimensional flows some adaptation of this method will need to take place. However for viscous three-dimensional flows the method as presented in this work should be directly applicable. 129  5.1. Significance of Results  5.1.2  Globalized Matrix-Explicit Newton-GMRES  The manuscript in Chapter 3 introduces improvements to two aspects of the efficient steady-state solution process developed for the high-order finitevolume scheme. The first aspect deals with improving the robustness and efficiency of the globalization of the Newton method. The second aspect deals with the efficiency of the solution of the linear systems arising from the globalized Newton method. As part of the improvement in the globalization, a method of combining timestepping with the line search method was proposed to improve robustness. Although both timestepping and line search methods are commonly used, this is to the author’s knowledge the first study combining both methods in a CFD application. The lack of robustness of the Newton-GMRES method, particularly for flows with shocks, has prevented its rapid adoption in production codes. For high-order schemes, where the numerical method is less dissipative, the problem is exacerbated. Although far from being a full solution to this issue, the combined line search and timestepping method provides an additional tool in the arsenal that researchers and developers can use to improve the robustness of the solution process. Another improvement to the globalization method presented in this work relates to the choice of local timestep. Selecting a timestep that depends on the relative local residual at the previous iterate was shown to be significantly more effective than using a CFL based local timestep or a global timestep. For transonic cases, in particular, the convergence of the scheme was drastically improved. The idea behind the use of the residual based timestep is simple. Regions where the residual, or flux integral, is large are changing rapidly and therefore require the use of small timesteps for the correct evolution to take place. Regions where small residuals prevail, on the other hand, benefit from large timesteps so that errors can be propagated to the boundaries quickly. An analogy to multigrid would be that the goal of the residual based local timestep scheme is to smooth high frequency error in regions where the residual at the previous iterate was large and to smooth low frequency error in regions where the residual was small. Further research is warranted to  130  5.1. Significance of Results determine if this approach is effective in reducing or eliminating the grid size dependence on the number of nonlinear iterations needed for convergence. However this is certainly not a complete replacement for multigrid since the stiffness of the resulting linear system maintains its grid size dependence. Both improvements to the globalization technique are not tied to the finitevolume method and could potentially be used with other discretizations. Chapter 3 also demonstrates that the exact Jacobian for the high-order finite-volume method can be effectively used to improve the performance of the preconditioned GMRES linear solver. The benefit is most significant for the fourth-order scheme for which the first-order Jacobian had previously been shown to provide poor preconditioning [17]. The high-order Jacobian not only improves conditioning but also eliminates the need for residual evaluations resulting from the matrix-free GMRES method previously employed. Perhaps the most significant revelation from this work is the demonstration that the exact Jacobian, which was previously thought to be prohibitively expensive to compute, can be obtained at reasonable cost in both computational time and memory. The significance of this extends beyond its demonstrated application for efficient ILU preconditioned GMRES. For example, the availability of the explicit Jacobian introduces the possibility of using an algebraic multigrid method to solve the linear system. Another potential use of the exact Jacobian is for aerodynamic shape optimization. Adjoint-based methods have been proposed for this application [2, 12, 14]. Currently, the discrete adjoint method is favored over the continuous adjoint method since it does not require that the flux Jacobian be explicitly known. However, in light of the fact that the exact flux Jacobian was found to be relatively inexpensive to compute, the continuous adjoint methods may prove to be more efficient. This research question has begun to be addressed by another graduate student in the author’s group.  5.1.3  Viscous Discretization  The extension of the least-squares based unstructured high-order finitevolume discretization to the full Navier-Stokes equations was carried out  131  5.2. Outlook in Chapter 4. Since, unlike the inviscid fluxes of the Euler equation, the viscous fluxes require the solution gradient in addition to the solution value, this represents an important extension to the discretization which warranted study. In addition to demonstrating that the expected order of accuracy could be achieved, it was also shown that such schemes can readily be converged through the use of the Newton-Krylov method with the help of the exact Jacobian, which can be obtained in a method similar to that described for the inviscid case in Chapter 3. This is, to the author’s knowledge, the first high-order unstructured leastsquares reconstruction based finite-volume discretization to demonstrate high-order accuracy for viscous flows. However, the discontinuous Galerkin method has also been successfully used for this [4, 5, 8]. A comparative study of the accuracy and efficiency of the high-order finite-volume method with the discontinuous Galerkin method would certainly be warranted. However, since both methods are in their infancy, such a study would provide only preliminary information as to which scheme could fare better on more complex problems which are ultimately of interest.  5.2 5.2.1  Outlook Turbulent Flows  The Reynolds averaged Navier-Stokes (RANS) equations provide a computationally feasible way of solving typical aerodynamics flows. For this application, the most common closure models are the one-equation SpalartAllmaras [21] and two-equation k − ω and k − models [23]. These closure equations, consisting of convective, diffusive, and source terms, are difficult to solve due to their numerical stiffness. Their numerical solution must be carefully designed to maintain the positivity of turbulence variables. For this reason it is common for codes which use a second-order discretization for the RANS equations to use a first-order discretization for some or all of the terms of the turbulence model. While good results have been achieved with this method, it is anticipated that a first-order discretization will be a limiting  132  5.2. Outlook factor on the accuracy of the overall scheme if high-order accurate methods are used for the RANS equations. Therefore a stable high-order accurate discretization of the turbulence model needs to be developed. Although this question has not been explored for high-order accurate finite-volume schemes, it has been to be addressed in the context of the discontinuous Galerkin method [4] which may provide valuable insight. The use of Newton-Krylov methods for the solution of turbulent flows will also require careful study. For second-order methods, the question has been addressed with some success. Due to numerical stiffness, some evidence suggests that a fully coupled approach to the solution of the turbulence and RANS equations may not provide the best convergence properties. For this reason a loosely-coupled approach has been suggested [6]. This approach is also attractive due to its reduced memory requirement resulting from the reduced block size of the Jacobian matrix. In the fourth-order case, where the use of the full-order Jacobian was shown to be important for good convergence, this savings in memory would be especially important. Another complication is that the turbulence model equations, unlike the Navier-Stokes equations, contain source terms. The high-order Gauss integration rules for these terms on arbitrary control volumes, including those at curved boundaries, have been developed [18]. However, in their current state, these rules make no guarantee that the Gauss points will be inside the control volume itself or even within the domain boundary. For simple applications, such as initializing a solution, transferring a solution to another grid, or evaluating the error relative to an analytic solution, this is not problematic. In fact, this integration method is already being successfully used as part of the pre- and post-processing of the solution in our code. However if applied to the integration of turbulence source terms, problems may arise near the domain boundary. For control volumes near a no-slip wall some Gauss points will be located outside the domain on the other side of the wall. The evaluation of turbulence model source terms at these points is ill defined since they depend on the distance to the nearest wall and the velocity magnitude. Therefore some special methods may need to be developed for high-order integration of turbulence source terms near 133  5.2. Outlook wall boundaries. The added complications of high Reynolds number flows are not limited to the issues directly related to the turbulence model. Since the boundary layer in these flows is very thin relative to the macroscopic flow properties, it is essential to use highly anisotropic grids to achieve computational efficiency. For aerodynamics applications, it is not unusual to use grids with cell aspect ratios of 10, 000 : 1 in the boundary layer. This is expected to cause problems with the stability, accuracy, and convergence properties of the scheme, all of which have also been issues with second-order schemes. For example, a relatively recent paper [16] addresses the accuracy of second-order least-squares reconstruction on anisotropic grids. Although the findings indicate that the use of distance weighing of the leastsquares terms is essential for accuracy, this weighing is not typically used due to stability concerns. In practice, structured or quasi-structured grids are almost always used in the boundary layer. The resulting grid alignment with flow properties alleviates many of the accuracy and stability issues of reconstruction based finite-volume schemes. Whether these accuracy and stability issues will be exacerbated with the use of a high-order accurate reconstruction needs to be studied. The difficulty of converging solutions to steady state using Newton-Krylov methods on anisotropic grids results in part from poor conditioning due to the strong directional coupling of flow properties. This problem will be aggravated by high-order methods since they have higher inherent numerical stiffness than second-order schemes. For second-order finite-volume methods, this has been addressed through the use of directional coarsening multigrid or the use of directional implicit smoothing [15]. Similar approaches will undoubtedly need to be considered for high-order schemes. A final challenge of highly anisotropic grids, which is unique to highorder schemes, arises due to the treatment of curved boundaries. To obtain high-order accuracy it is essential to account for the curved boundary when locating flux integration points and when computing the control volume moments used for reconstruction. However, for highly-anisotropic grids, this may lead to ill-defined control volumes where the boundary face bulges past 134  5.2. Outlook  (a) straight boundary  (b) curved boundary  Figure 5.1: Effect of curved boundary treatment on a quadrilateral highly anisotropic boundary control volume the interior face as seen in Figure 5.1. The solution to this problem may be to also curve the interior faces of anisotropic control volumes near the boundary. For this, a method to compute the appropriate location of the interior Gauss points used for flux integration and moment calculation needs to be devised.  5.2.2  Three Dimensional Flows  The fundamental aspects of the high-order discretization do not change substantially from the two-dimensional case to the three-dimensional case. Indeed, the discretization has already been implemented and shown to be accurate for simple problems such as the Poisson equation. However, the challenge of implementing curved boundary support has not yet been addressed. The Gauss integration rules developed for curved boundary control volumes [18] will need to be extended to the three-dimensional case. Furthermore, for turbulent flows, the schemes devised to accommodate highly anisotropic curved boundary control volumes will need to be extended to three dimensions. In terms of computationally efficiency, there is good news and bad news for the extension of the method to three dimensions. The good news is that the benefit of high-order accuracy, in terms of accuracy versus number of control volumes, should increase. For example, let’s consider the results in Chapter 4 where a fourth-order solution of the airfoil problem on a mesh with 5, 363 control volumes was more accurate than that of a second-order scheme on a mesh with 37, 399 control volumes. To achieve similar mesh spacing in three dimensions the number of control volumes needs to be increased by a 135  5.2. Outlook power of 32 . Therefore, the equivalent mesh size in three dimensions would be 390, 000 control volumes for the fourth-order scheme and 7, 200, 000 control volumes for the second-order scheme. While the two-dimensional problem required roughly 7 times fewer control volumes for the fourth-order solution, the three-dimensional problem may require 18 times fewer control volumes. The bad news is that the computational cost of the high-order scheme relative to the second-order scheme on a per control volume basis also increases. In all likelihood the median dual vertex centered control volume, which we prefer in two dimensions due to conditioning issues, will need to be abandoned due to the large number of faces of such a tessellation in three dimensions. Even if a cell-centered approach is taken, which results in only four faces per control volume, the relative increase in cost of moving to a high-order accurate flux integration in three dimensions exceeds that of two dimensions. In two dimensions a fourth-order scheme requires twice as many Gauss integration points as a second-order scheme, while in three dimensions this factor increases to four. Furthermore, in terms of memory requirement, the additional cost of the matrix-explicit method relative to the matrix-free method will also increase. This is due to the increase in Jacobian block size from the additional velocity term, as well as to the increased flux dependency stencil which determines the matrix sparsity. Due to this plethora of differences in the three-dimensional case, it is difficult to speculate on the exact potential of high-order methods for these flows.  136  5.3. Bibliography  5.3  Bibliography  [1] R. Abgrall. On essentially non-oscillatory schemes on unstructured meshes: Analysis and implementation.  Journal of Computational  Physics, 114(1):45–58, 1994. [2] W. Anderson and V. Venkatakrishnan. Aerodynamic design optimization on unstructured grids with a continuous adjoint formulation. Computers and Fluids, 28(4–5):443–480, 1999. [3] T. J. Barth and D. C. Jespersen. The design and application of upwind schemes on unstructured meshes. AIAA paper 89-0366, Jan. 1989. [4] F. Bassi, A. Crivellini, S. Rebay, and M. Savini. Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ω turbulence model equations. Computers and Fluids, 34(4-5):507 – 540, 2005. Residual Distribution Schemes, Discontinuous Galerkin Schemes and Adaptation. [5] F. Bassi and S. Rebay. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. Journal of Computational Physics, 131:267–279, 1997. [6] M. Blanco and D. W. Zingg. A Newton-Krylov algorithm with a looselycoupled turbulence model for aerodynamic flows. AIAA Paper 2006-0691. Presented at the 44th AIAA Aerospace Sciences Meeting, Jan. 2006. [7] M. Delanaye. Polynomial Reconstruction Finite Volume Schemes for the Compressible Euler and Navier-Stokes Equations on Unstructured Adaptative Grids. PhD thesis, Universite de Liege, Faculte des Sciences Appliquees, 1998. [8] K. J. Fidkowski, T. A. Oliver, J. Lu, and D. L. Darmofal. p-multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations. Journal of Computational Physics, 207(1):92 – 113, 2005.  137  5.3. Bibliography [9] O. Friedrich. Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids. Journal of Computational Physics, 144(1):194–212, July 1998. [10] A. G. Godfrey, C. R. Mitchell, and R. W. Walters. Practical aspects of spatially high-order accurate methods. American Institute of Aeronautics and Astronautics Journal, 31(9):1634–1642, Sept. 1993. [11] C. Q. Hu and C. W. Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150(1):97–127, Mar. 1999. [12] A. Jameson, L. Martinelli, and N. A. Pierce. Optimum aerodynamic design using the Navier-Stokes equations. Theoretical and Computational Fluid Dynamics, 10(1–4):213–237, 1998. [13] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted ENO schemes. Journal of Computational Physics, 126(1):202–228, June 1996. [14] D. Mavriplis. A discrete adjoint-based approach for optimization problems on three-dimensional unstructured meshes. AIAA Paper 2006-0050. Presented at the 44th AIAA Aerospace Sciences Meeting, Jan. 2006. [15] D. J. Mavriplis. On convergence acceleration techniques for unstructured meshes. Technical Report ICASE No. 98-44, Institute for Computer Applications in Science and Engineering (ICASE), NASA Langley Research Center, NASA Langley Research Center, Hampton VA 23681-2199, 1998. [16] D. J. Mavriplis. Revisiting the least-squares procedure for gradient reconstruction on unstructured meshes. In Proceedings of the 16th AIAA Computational Fluid Dynamics Conference, 2003. [17] A. Nejat and C. Ollivier-Gooch. Effect of discretization order on preconditioning and convergence of a high-order unstructured NewtonGMRES solver for the Euler equations. Journal of Computational Physics, 227(4):2366–2386, 2008.  138  5.3. Bibliography [18] C. Ollivier-Gooch. On k-exact integration over arbitrary polygons and polyhedra. ACM Trans. Math. Softw., In preparation, 2009. [19] C. F. Ollivier-Gooch. Quasi-ENO schemes for unstructured meshes based on unlimited data-dependent least-squares reconstruction. Journal of Computational Physics, 133(1):6–17, 1997. [20] C.-W. Shu, T. A. Zang, G. Erlebacher, D. Whitaker, and S. Osher. Highorder ENO schemes applied to two- and three-dimensional compressible flow. Applied Numerical Mathematics, 9(1):45–71, 1992. [21] P. Spalart and S. R. Allmaras. A one-equation turbulence model for aerodynamic flows. AIAA paper 92-0439, Jan. 1992. [22] V. Venkatakrishnan. On the accuracy of limiters and convergence to steady-state solutions. AIAA paper 93-0880, Jan. 1993. [23] D. C. Wilcox. Turbulence Modeling for CFD. DCW Industries, second edition, 2004.  139  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0067108/manifest

Comment

Related Items