HIGH-ORDER FINITE-VOLUME DISCRETISATIONS SOLVING A MODIFIED ADVECTION-DIFFUSION ON UNSTRUCTURED TRIANGULAR PROBLEM MESHES by MICHAEL VAN ALTENA B.Math., The University of Waterloo, 1997 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Mechanical Engineering) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A October 1999 ©Michael Van Altena, 1999 FOR In presenting degree this at the thesis in partial fulfilment of University of British Columbia, I agree freely available for reference copying of this department publication or of and study, this his or her representatives. A1eeWo.ro Ga\ bf\c{\Aeer"ir\g The University of British Columbia Vancouver, Canada Date DE-6 (2/88) Oc-f \x rm that the may be It thesis for financial gain shall not permission. Department of requirements i further agree thesis for scholarly purposes by the is that an advanced Library shall make it permission for extensive granted by the understood be for head that allowed without of my copying or my written Abstract A fourth-order accurate diffusive flux calculation scheme has been developed which can be incorporated into any existing reconstruction-based unstructured mesh inviscidflowsolver. The flux scheme makes use of a Laplace problem to model the viscous terms and of the existing least-squares reconstruction code used in Dr. Carl F. Ollivier-Gooch's inviscid flow solver. Boundary conditions on the solution are enforced during the reconstruction, which requires boundary constraints to be included in the least-squares matrix. The Taylor Constraints method for generating boundary constraints is developed and then discarded due to its drawbacks. The Gauss Constraints method is then developed which is simpler and avoids the drawbacks of the Taylor Constraints method. The flux integration is validated to demonstrate fourth-order accuracy on square, circle and cardioid meshes with a combination of Dirichlet and Neumann boundary conditions. The flux integration scheme is shown to be fourth-order accurate for a modified advection-diffusion problem of the same form as the Navier-Stokes' momentum equations. The diffusive flux calculation scheme can later be combined with an existing inviscid flow solver to solve the viscous compressible Navier-Stokes' equations. ii Contents Abstract ii List of Tables v List of Figures vi Acknowledgements ix 1 Introduction 1.1 1 Objectives and Scope 3 2 Flow Solvers 5 2.1 Finite-Volume Methods 5 2.2 Inviscid Flow Solver 11 2.3 Viscous Modelling 12 2.4 Hybrid Scheme 13 3 High-Order Flux Integrals 17 3.1 Gradient Approximation 18 3.2 Flux Integration: Cell-Centred Control Volumes 18 3.3 Flux Integration: Vertex-Centred Control Volumes 21 4 Boundary Conditions 4.1 24 Taylor Constraints 25 4.1.1 Dirichlet Boundary Conditions 25 4.1.2 Neumann Boundary Conditions 30 4.1.3 Dirichlet & Neumann Combined Boundary Conditions 34 4.1.4 Drawbacks 38 iii CONTENTS iv 4.2 Gauss Constraints . 39 4.2.1 Dirichlet Boundary Conditions 40 4.2.2 Neumann Boundary Conditions 40 4.2.3 Dirichlet & Neumann Combined Boundary Conditions 42 4.2.4 Constraint Comparison 44 4.3 Constrained Least-Squares Reconstruction Method 46 5 Validation 49 5.1 Interior Control Volume F l u x Integral Validation 51 5.2 Dirichlet Boundary Conditions on a Square Mesh 53 5.2.1 Accuracy Analysis — Taylor Constraints 53 5.2.2 Accuracy Analysis — Gauss Constraints 55 5.3 Neumann Boundary Conditions on a Square Mesh 5.4 Combined Dirichlet & Neumann Boundary Conditions on a Square Mesh 5.5 Singular Point on a Square Mesh 62 5.6 Curved Mesh Boundaries 63 5.6.1 Circle Mesh Boundary 65 5.6.2 Cardioid Mesh Boundary 67 5.7 Summary . . 6 Application 58 . . 60 70 71 6.1 Modified Advection-Diffusion Problem 71 6.2 A p p l y i n g the Finite-Volume Scheme 72 6.3 Validation 73 7 Summary and Conclusions 76 7.1 Issues During Development 78 7.2 Computer Resource Usage 80 7.3 Direction for Future Work 84 Bibliography 87 Appendix A 88 List of Tables 4.1 Constraint Comparison for Control Volumes with Taylor Constraints 4.2 Constraint Comparison for Control Volumes with Gauss Constraints 45 5.1 Error Norms of the Interior Control Volumes 52 5.2 Error Norms with Dirichlet Taylor Constraints 54 5.3 Error Norms with Dirichlet Gauss Constraints 56 5.4 Error Norms with Neumann Gauss Constraints 59 5.5 Error Norms with Gauss Constraints 61 5.6 Error Norms with Gauss Constraints 64 5.7 Error Norms on a Circle Mesh 66 5.8 Error Norms on a Cardioid Mesh 69 6.1 Error Norms for the Modified Advection-Diffusion Problem 75 7.1 Mesh Information for Square Meshes 80 7.2 Memory Data for Square Meshes 81 v .... 38 List o f Figures 2.1 Cell-Centred Control Volume 7 2.2 Vertex-Centred Control Volume 7 2.3 High-order reconstruction with vertex-centred control volumes 8 2.4 Hybrid Scheme 15 3.1 Cell-Centred Flux Calculation 19 3.2 Cell Gradients 20 3.3 Cell Average Gradients 20 3.4 Vertex-Centred Gauss Points 23 4.1 Cell-centred control volume's Gauss point locations 45 4.2 Vertex-centred control volume's Gauss point locations 45 4.3 Example of the least-squares matrix 48 4.4 Example of the constrained least-squares matrix 48 5.1 Error Norm Plot Legend 51 5.2 Initial temperature function for interior control volumes 52 5.3 Exact flux for interior control volumes 52 5.4 L error norm plots for interior cell-centred control volumes 52 5.5 L2 error norm plots for interior vertex-centred control volumes 52 5.6 Flux error for Cell scheme using Dirichlet Taylor constraints 54 5.7 Flux error for Vertex scheme using Dirichlet Taylor constraints 54 5.8 L error norm plots for Cell scheme using Dirichlet Taylor constraints . . . . 54 5.9 L error norm plots for Vertex scheme using Dirichlet Taylor constraints . . . 54 2 2 2 5.10 Flux error for Cell scheme using Dirichlet Gauss constraints 56 5.11 Flux error for Vertex scheme using Dirichlet Gauss constraints 56 5.12 L error norm plots for Cell scheme using Dirichlet Gauss constraints . . . . 56 2 vi LIST OF FIGURES vii 5.13 L error norm plots for Vertex scheme using Dirichlet Gauss constraints . . . 56 5.14 Initial function for Vertex scheme with Neumann boundary conditions . . . . 58 5.15 Flux error for Cell scheme using Neumann Gauss constraints 59 5.16 Flux error for Vertex scheme using Neumann Gauss constraints 59 5.17 L error norm plots for Cell scheme using Neumann Gauss constraints . . . . 59 5.18 L error norm plots for Vertex scheme using Neumann Gauss constraints . . 59 5.19 Initial function for Vertex scheme for combined boundary conditions 60 5.20 Flux error for Cell scheme using combined Gauss constraints 61 5.21 Flux error for Vertex scheme using combined Gauss constraints 61 5.22 L error norm plots for Cell scheme using combined Gauss constraints . . . . 61 5.23 L error norm plots for Vertex scheme using combined Gauss constraints . . 61 5.24 Initial function of singular point test case on cell-centred control volumes . . 63 5.25 Flux error for singular test case using Cell scheme with Gauss constraints . . 64 5.26 Flux error for singular test case using Vertex scheme with Gauss constraints 64 5.27 L error norm plots for Cell scheme using Gauss constraints 64 5.28 L error norm plots for Vertex scheme using Gauss constraints 64 5.29 Initial function for cell-centred control volumes on a circle mesh 65 5.30 Flux error for Cell scheme on a circle mesh 66 5.31 Flux error for Vertex scheme on a circle mesh 66 2 2 2 2 2 2 2 5.32 L error norm plots for Cell scheme on a circle mesh 2 . 66 5.33 L error norm plots for Vertex scheme on a circle mesh 66 5.34 Initial function for Vertex scheme on a cardioid mesh 68 5.35 Flux error for Cell scheme on a cardioid mesh 69 5.36 Flux error for Vertex scheme on a cardioid mesh 69 5.37 L error norm plots for Cell scheme on a cardioid mesh 69 5.38 L error norm plots for Vertex scheme on a cardioid mesh 69 6.1 Initial function for advection-diffusion problem 74 6.2 Flux error for Cell scheme for advection-diffusion problem 75 6.3 Flux error for Vertex scheme for advection-diffusion problem 75 6.4 L error norm plots for Cell scheme for advection-diffusion problem 6.5 L error norm plots for Vertex scheme for advection-diffusion problem . . . . 75 7.1 Plot of number of control volumes vs. amount of memory 81 2 2 2 2 2 75 LIST OF FIGURES viii 7.2 Plot of L error norm vs. CPU time for Cell scheme 83 7.3 Plot of L error norm vs. CPU time for Vertex scheme 83 2 2 A.l squareO.mesh 88 A.2 squarel.mesh 89 A.3 square2.mesh 89 A.4 circleO.mesh 90 A.5 circlel.mesh A.6 cardioidO.mesh . . . 90 91 Acknowledgements I would like to thank my research supervisor, Dr. Carl F. Ollivier-Gooch, for his constant support and help throughout the course of this thesis. Being the guinea pig graduate student was a pleasure, and his future students will be lucky to have him. I would also like to thank Jody Kirchner; her friendship throughout my graduate school experience I will never forget. ix Chapter 1 Introduction The ultimate goal of the computational fluid dynamics (CFD) field is to obtain accurate solutions to engineering problems with a rapid turnaround time. This is done by finding an approximate solution to the equations which govern the dynamics of fluids, known as the Navier-Stokes' equations, and their prescribed boundary conditions. A fluids problem is solved numerically by first decomposing the computational domain into a mesh. The Navier-Stokes' equations are then discretised and solved over the mesh. The accuracy of the solution and the efficiency with which it is obtained are very important. The computer program should compute an accurate solution in the least amount of time possible. The geometry of a problem can be as simple as a flat plate, or as complex as a space shuttle. This domain must be decomposed into a mesh for the flow solver. In general there are two types of meshes, structured and unstructured. In a structured mesh the identity of a cell's neighbours can be determined from the identity of the cell. The mesh cells are almost always quadrilaterals and hexahedra. Structured meshes are difficult and labour intensive to 1 CHAPTER 1. INTRODUCTION 2 generate for complex geometries, due to difficulties in meshing around holes, sharp corners, etc. Typically, different regions of a complex domain are meshed separately and then the component meshes are sewn together. Meshing can be greatly simplified by using unstructured meshes. An unstructured mesh differs from a structured mesh in that the neighbours for any particular cell cannot be computed, but must be stored in a data structure in the computer's memory. An unstructured mesh can be constructed with a variety of polygons in two dimensions (usually triangles and/or quadrilaterals) or polyhedra in three dimensions (usually tetrahedra and/or hexahedra). The topological flexibility of unstructured meshes makes it easier to mesh complex geometries consistently. This thesis uses two dimensional triangular meshes exclusively, although the results should generalize easily to other types of meshes. With the mesh generated, the governing equations need to be discretised to compute the solution over the mesh. One common approach to discretisation in CFD is the finite-volume method, which performs fluids control volume analysis by calculating thefluxesthrough the many small control volumes in a mesh. The development of such a method is at the heart of this thesis. The two most important issues in developing a numerical flow solver are efficiency and accuracy. The problem is that efficiency and accuracy are inversely related. Most flow solvers compute solutions in the control volumes to second-order accuracy, as this can be done quickly. To improve the accuracy of the solution, the mesh is made finer, increasing the number of control volumes while decreasing their size. However the larger number of control volumes also means a larger number of calculations, which require more CPU time. CHAPTER 1. INTRODUCTION 3 Another way to improve accuracy is to use the original coarse mesh and compute solutions in the control volumes more accurately. This increases the computation time per control volume, but uses fewer control volumes than mesh refinement. Which method is more efficient? The trade-off between computational accuracy, mesh refinement and efficiency for viscous flow solvers on unstructured meshes must be quantified. To study this trade-off, a method must first be developed to solve the Navier-Stokes' equations to high-order accuracy on unstructured meshes. This work takes the first step towards building a high-order accurate flow solver by developing a high-order accurate flux integration scheme for diffusive fluxes. 1.1 Objectives and Scope The main objective of this research is to develop and validate a fourth-order accurate diffusive flux calculation scheme. The flux scheme will make use of a Laplace problem to model the viscous terms and of the existing reconstruction code used in Dr. Carl F . Ollivier-Gooch's inviscid flow solver. This diffusive flux calculation scheme can later be combined with the existing inviscid flow solver to solve the complete viscous compressible Navier-Stokes' equations. The existing inviscid flow solver is discussed in Chapter 2. Next, the viscous terms in the Navier-Stokes' equations are examined and a simple model problem is proposed. Then the commonly-used second-order Hybrid scheme for diffusive flux calculation is explained. This scheme cannot efficiently be extended to higher-order accuracy, motivating the development CHAPTER 1. INTRODUCTION 4 of new finite-volume methods. A new high-order finite-volume scheme for diffusive flux calculations is discussed i n Chapter 3. T h i s reconstruction-based scheme uses the same reconstruction data calculated for the inviscid terms. T h e description includes a detailed example of the flux integral calculation. Chapter 4 describes how the boundary conditions are enforced during reconstruction. Two methods are discussed for deriving the equations to constrain the reconstructed solution near the boundary. The method of imposing these constraints is then explained. The newly developed finite-volume method is validated in Chapter 5. Methodology for measuring the accuracy of the method is described. The method is tested for the interior control volumes on a simple square mesh. Next, several test cases are shown using square, circle and cardioid meshes to prove the scheme works using either cell-centred or vertexcentred control volumes at second-, third- and fourth-order accuracy. A comparison is made between the Taylor Constraints and Gauss Constraints methods for creating the boundary constraint equations. A t the end of the chapter an example of a problem w i t h a singularity is discussed. Finally, i n Chapter 6, the control volume finite-volume method is applied to a modified advection-diffusion problem. T h i s scalar model problem has the same structure as the Navier-Stokes' momentum equations, to prove the scheme can be used to evaluate the viscous terms i n the Navier-Stokes' equations. The last chapter summarises the work done and presents conclusions about the developed schemes. Chapter 2 Flow Solvers A method is to be developed to extend a high-order accurate inviscidfinite-volumeflow solver to a high-order accurate viscous flow solver. This chapter summarises how finitevolume methods work in general before going into the inner workings of the existing inviscid flow solver. Next, a model problem for the viscous terms is discussed, since the viscous terms themselves are too complex to study easily. Finally, the commonly used Hybrid finite-volume method is explored for computing high-order fluxes. 2.1 Finite-Volume Methods Finite-volume methods are based on the fully conservative control volume form of the governing equations being solved. These methods compute control volume averages of the unknowns, rather than computing values at the mesh vertices. Both time advance and steady statefinite-volumemethods compute the net flux of the 5 CHAPTER 2. FLOW SOLVERS 6 unknown variables through the boundary of each control volume. Thefluxintegral can be used for both time advance and relaxation schemes. One method to accurately advance the solution in time is to use thefluxintegral with the fourth-order accurate Runge-Kutta scheme to advance the mean values in the control volumes. For steady-state problems the flux integral can be used to drive a successive over-relaxation process. The accuracy of the solution hinges on the accuracy of the flux integral. A high-order accuratefluxintegral in a control volume requires an accuratefluxand accurate integration. An accuratefluxrequires an accurate approximation of the solution in the control volume. One way to achieve this is to reconstruct a polynomial approximation to the solution in the control volume. Before describing a reconstruction method, two types of control volumes a finite-volume scheme can use to break up an unstructured mesh are described. Thefinite-volumeschemes described in this thesis can use either cell-centred control volumes or vertex-centred control volumes. Cell-centred control volumes are defined as the cells of the mesh. The perimeter of vertex-centred control volumes are defined as lines surrounding a vertex, connecting the midpoint of the triangle edges to the centroid of the cells. Examples of the two control volumes are shown in Figures 2.1 and 2.2. Thefinite-volumeschemes use a local reference point in the control volume as the origin for the reconstruction polynomial. These local reference points are the centroid of the triangle and the encompassed vertex, for the cellcentred and vertex-centred control volumes respectively. Theflowsolver approximates the inviscid terms' unknown variables in the control volumes with the cubic Taylor polynomial centered about the control volume's reference point CHAPTER 2. FLOW 7 SOLVERS Figure 2.1: Cell-Centred Control Volume T(x,y) dx l c dT Figure 2.2: Vertex-Centred Control Volume dT 3 c c (x-x ) (y-y ) dx dy c 2 c dy + c 2 (y - Vcf , 0»T + (x - x f c dx 3 (x - x ) (y - y y aT dxdy 3 | c dx (y-y ) 2 2 c 2 b*T dT xy 2 c 2 {x - x ) (y - y ) + (x - 2 (y - y ) + c 2 dxdy dT , dT (x - x ) + — dy c + + (2.1) dy 3 Here and throughout this thesis, the unknown variable will be represented by the temperature T. The coefficients of the Taylor polynomial are determined by reconstructing the solution using a least-squares reconstruction method [17, 18]. This reconstruction method is used in this thesis. The coefficients of the Taylor polynomial approximation are chosen to predict the nearby control volume averages well. The Taylor polynomial for the control volume where reconstruction is being done is averaged analytically over near control volumes, and the CHAPTER 2. FLOW SOLVERS # Reconstruction C V • Added for second order ^ Added for third order # Added for fourth order Figure 2.3: High-order reconstruction with vertex-centred control volumes on a partial mesh. The line along the bottom represent a boundary of the mesh. The neighbouring control volumes used during reconstruction are highlighted for different orders of accuracy for two control volumes. results are compared with the real averages. The equations are weighted according to how close the nearby control volume is to the one being reconstructed. The closer neighbours are more important (i.e. weighted heavily) than the control volumes that are further away. The number of neighbouring control volumes used to create the system of equations depends on the order of accuracy desired. Figure 2.3 shows the vertex centred-control volumes used for each order of accuracy. As the order of accuracy increases, the neighbouring control volumes CHAPTER 2. FLOW SOLVERS 9 of the lower-orders are also taken into account. For example, for a fourth-order accurate reconstruction, all the neighbours of the reconstructed control volume, at the bottom of the mesh shown in Figure 2.3, for second-, third- and fourth-order accuracy are taken into account. By design, this process creates a system of equations with more equations than there are unknowns for the desired order of accuracy. The system of equations is over-determined, which is why a least-squares algorithm is used to solve for the unknown coefficients. A mean constraint equation is also added to the system before solving, to ensure the mean value in the reconstructed control volume is conserved — a requirement of finitevolume schemes. The mean constraint is: dT dx ^xxx 6 & T dx + 3 dT &T /** d T Iyy d T + L;xy + + dy 2 dx dxdy 2 dy T dT dT I dT xxy lyyy U ± 2 dxdy + 6 dy 2 dx dy + xyy 2 2 + Iy 2 2 2 3 3 3 L u 2 ± + 2 3 (2.2) where T is the unknown value and the L's are the control volume's moments, for example: xxy l = jjx ydA 2 (2.3) cv where A is the area of the control volume. The least-squares matrix contains the mean constraint equation, but this alone does not ensure the mean value in the control volume will be conserved. For the mean constraint to be satisfied exactly, the first column of the matrix is reduced using Gaussian elimination with the mean constraint equation. This ensures the mean constraint is enforced during reconstruction. The least-squares algorithm is then used on the column-reduced matrix to solve for the unknown coefficients. The Taylor polynomial with the reconstructed coefficients gives an CHAPTER 2. FLOW SOLVERS 10 approximation to the solution in the control volume. This process is repeated for each control volume in the mesh. Flux integrals are formed from the Taylor polynomial approximation of the solutions in each control volume. The flux integrals are evaluated using Gauss quadrature [21]. Gauss quadrature gives the capability of evaluating a definite integral to high-order accuracy with the integrand evaluated at only a few points. A n alternate approach to reconstruction is the essentially non-oscillatory (ENO) scheme, which ensures that all control volumes with a smooth neighbourhood have a uniformly high-order accurate reconstruction. Some E N O schemes search for the smoothest stencil for reconstruction in each control volume [1, 6, 22]. These schemes can have difficulty with steady-state calculations. To overcome these problems weighted E N O (WENO) schemes were developed [7, 10]. Instead of selecting the smoothest stencil, they use a weighted sum of all the stencils. E N O reconstruction is compatible with the method developed in this thesis for evaluating viscous terms. In summary, the general steps for finite-volume methods are as follows: 1. Reconstruct a high-order approximation to the solution, based on the current mean values in local control volumes. 2. Compute the high-order flux integral through the boundaries of the control volumes. 3. Use a time advance or relaxation scheme to advance the current solution with the flux integral. 4. Repeat all steps until desired time or steady-state solution is reached. 11 CHAPTER 2. FLOW SOLVERS 2.2 Inviscid Flow Solver Dr. Carl F. Ollivier-Gooch has developed a fourth-order accuratefinite-volumeflow solver for the inviscid terms of the compressible Navier-Stokes' equations: ^ + V-(pu) = 0 ^ ^ ^ where + V-(pwg>w) + V-(u(E + P)) (u <g> i t ) . i F B (2.4) = -VP + F (2.5) = (2.6) B 0 = UiUj = net vector body force for the control volume In inviscid compressibleflow,the onlyfluxnormal to a boundary wall is thefluxdue to pressure. The pressure is obtained by extrapolation and used explicitly in thefluxcalculation at the boundary wall. This is all that is needed for boundary conditions at a wall in inviscid compressibleflow.Inflow and outflow boundary conditions are based on characteristic theory. Roe's scheme [20] is used to compute thefluxat the control volume boundary to match the compressible inviscid physics The netfluxesare then computed, and an appropriate time advance scheme is used to update the solution, completing a single time step. In addition to upwind methods, such as the one used in the existing inviscid flow solver [3, 4, 5, 17, 18], there are other successful methods based on centred differences [12, 13, 15, 16]. Centred difference codes do not perform reconstructions, so the method for viscous discretisation presented here will not fit into that context. CHAPTER 2. FLOW SOLVERS 2.3 12 Viscous Modelling A high-order viscousflowsolver is to be developed by extending the methods in the existing inviscidflowsolver to include the viscous terms of the governing equations. The full viscous compressible Navier-Stokes' equations can be written as given in equations 2.7 to 2.9. The boxed terms in equations 2.8 and 2.9 represent the viscous terms to be evaluated. ^ + V-(pu) = 0 (2.7) d(pu) + V-(pw<g>u) = - V P + |V -TI+FB dt dE i 1 i 1 -—- + V • (u (E + P)) = V • (r • u) + V • (KVT) dt 1 where 1 1 (2.8) (2.9) 1 (u <8> w) • = UiUj i nj = 2ti ^Dij - ^5ijD ^j mT1 D = 5 ( V . (V«) ) T + o r D, 3= 5 ( ^ + ^ Fi3 = net vector body force for the control volume The viscous terms themselves are too complex to study easily since it is a coupled system involving two unknowns, so a model equation is used. As a model problem the Laplace equation for temperature is used with no forcing term: V T =0 2 (2.10) This equation is integrated over each control volume in the mesh and the Divergence Theorem is applied to obtain afluxintegral: / d(CV) VT-nds = 0 (2.11) CHAPTER 2. FLOW SOLVERS 13 where h is the unit normal vector to the shared control volume edge. To solve the homogeneous Laplace equation to high-order accuracy, the flux integral must be calculated to high-order accuracy. Because the integrand contains the temperature gradient, VT must be known to high-order accuracy. The core of this thesis describes how to calculate the flux integral J VT • n ds and 8(CV) its boundary conditions to high-order accuracy. The types of boundary conditions discussed are Dirichlet alone, Neumann alone and linear combinations of Dirichlet and Neumann. 2.4 Hybrid Scheme Before discussing the new high-order method, this section looks at one of the most common finite-volume methods for discretising the Laplace equation. Such a method is the Hybrid scheme, so named because the control volumes are vertex-centred and the gradients are calculated over the cells of the mesh. This scheme calculates thefluxthrough vertex-centred control volumes. Many researchers use essentially this scheme, often for computing viscous terms of the Navier-Stokes' equations, using eitherfinite-elementorfinite-volumemethods. The finitevolume community seems to use second-order accuracy exclusively [2, 13, 14, 15, 16], and evenfinite-elementresearchers do not typically exploit high-order accuracy [8, 9, 11, 19]. For solving the Laplace equation, as shown in the previous section, the temperature gradient is needed to calculate thefluxthrough the control volume boundary. The Hybrid scheme approximates the gradient of the solution on the triangles of the mesh and uses it CHAPTER 2. FLOW SOLVERS 14 to calculate thefluxthrough the control volume edges inside the triangle. The gradient is approximated by applying Green's Theorem: /(!£ " f ) CV d A = (^* / + (2.12) Qdv) 9(CV) By letting both Q and P be equal to temperature T, equation 2.12 can be broken up into: Sf h iA= iy x CV d(CV) -f% f iA= CV Tix (213) d(CV) The temperature gradient is assumed constant throughout the control volume, which causes a second-order error. This assumption allows the integral equations in 2.13 to be written as: — -— I lh ~A J Tdy 8{CV) dT dJi I f = -A } d(CV) T d X ( 2 ' 1 4 ) These equations provide the needed temperature gradient for thefluxcalculation. The vertices and the area of the triangle are known. In order to evaluate the contour integrals, the Hybrid scheme assumes the control volume average value applies at the respective vertex, meaning the temperature located at vertex 0 is equal to the mean temperature T . This 0 assumption also causes a second-order error, since the vertex is not located at the centre of the control volume. Using the assumption, the temperatures at each of the triangle vertices CHAPTER 2. FLOW SOLVERS 15 are known. With the equations in 2.14, the gradient can be evaluated as: yi yi = 1/ ^ = UJTdy dx = + J Tdy + J Tdy yi i [(ti - t ) ( 0 f yo Vl yi yi - yo) + {T - t ) (y - ) + ( t - t ) (y - y )] 2 yi 2 Vl 0 2 0 (2.15) 2 yo y\ yi - [ ( t i - t ) (x - x ) + ( t - fx) (x - xi) + ( t - t ) (x - x )\ (2.16) 0 x 0 2 2 0 2 0 2 With the temperature gradient known for the triangle, the flux through the edges between the control volumes located inside the triangle is calculated, using the flux integral. Because the gradient is constant, the flux integral is simply the dot product of the gradient and the unit normal to the shared edge, multiplied by the length of the shared control volume edge. Figure 2.4: Hybrid Scheme Consider the situation in Figure 2.4. The gradient in the shaded triangle is approximated to second-order accuracy using equations 2.15 and 2.16. This gradient is then used in CHAPTER 2. FLOW SOLVERS 16 evaluating the flux integral through the shared edge of the vertex-centred control volumes inside the triangle. Thefluxesare then subtracted from the source control volume, and added to the destination control volume, as determined by the direction of the normal to the edge used in the flux integral calculation. In Figure 2.4 the flux calculated across the upper edge is subtracted from the current totalfluxfor control volume 0 and added to the current totalfluxfor control volume 2. Similarly the totalfluxesfor control volumes along the remaining two edges are adjusted. Once this process has been repeated for each cell in the mesh, the flux integral for each control volume is known. Thefluxesare then used with an appropriate time advance scheme to update the mean values in their corresponding control volumes. Improvement of the accuracy of this scheme requires high-order gradients, which implies a reconstruction over the triangles. This means if the high-order extension of the Hybrid method were applied to solving the Navier-Stokes' equations, two reconstructions would need to be performed: one over the vertex-centred control volumes for the inviscid terms and one over triangles for the viscous terms. For this reason attempting to improve the accuracy of the Hybrid scheme efficiently would be difficult. It was decided to develop other high-order schemes, which would share the reconstructed solution from the inviscid terms and thus be more efficient. Chapter 3 High-Order Accurate Gradient Approximation and Flux Integration This chapter introduces two reconstruction-based schemes that use the Taylor polynomial coefficients, found during reconstruction, to calculate the viscous fluxes through the control volume boundaries. A high-order accurate approximation to the flux integral J" VT • ft ds d(CV) requires a high-order accurate gradient. The existing reconstruction scheme in the inviscid flow solverfindsa cubic Taylor polynomial approximation to the solution expanded about the reference point of the control volumes. T(x,y) dx l c dT , (x -x ) c dr — oy + (x - x ) (y - y ) + c dxdy dT c (x-x ) (y-y ) z c 2 dT 3 {y - 3 r yf c dy dT 3 dxdy 2 c + 2 (y-y y 2 + dx dy dy dT Xf dx c , z c (x - 2 2 2 dT {y - y ) + ,2 , + &T c dx (x - x ) (y c (x - x y 3 y) + z c + (3.1) c 6 17 CHAPTER 3. HIGH-ORDER FLUX INTEGRALS 18 This equation is the basis for obtaining an approximation to the gradient of the solution. This chapter describes how the gradient is approximated and how the flux integral is calculated for both cell-centred and vertex-centred control volumes. 3.1 Gradient Approximation The solution gradient is approximated in the control volume by taking the gradient of the high-order Taylor polynomial. The gradient of the cubic Taylor polynomial (equation 3.1) gives the gradient of the solution within each control volume: dT 2 VT(*,y) « ( + f dT (x - x ) + c dx 2 dxdy dT s dx dy 2 dy {x - x ) (y - y ) + c dT c + (x - x ) + c dxdy 2 (v dxdy - (x - x ) [y - y ) + c c (x - x y c dx 3 + yf c 2 dT (X - X f c {y - yc) + dy 2 dT 3 3 dxdy {y - Vc) + 2 2 dT 3.2 3 2 3 dT dT dT dy 3 (y - Vcf dx dy 2 + (3.2) Flux Integration: Cell-Centred Control Volumes In this section, an example will be used to describe the method by which the flux integral is calculated for cell-centred control volumes. It is assumed that fourth-order accuracy is desired. Consider the two interior control volumes shown in Figure 3.1. The previous section described how the gradient is approximated in each control volume. The gradients used CHAPTER 3. HIGH-ORDER FLUX INTEGRALS 19 Figure 3.1: Cell-Centred Flux Calculation along edge ab in the flux integrals § VT • n ds and d(CVO) § VT - fi ds must be identical. 9(CV1) This requirement comes from global conservation: the flux leaving a control volume must enter its neighbour. However, along the shared edge between two control volumes there are two gradients, which are not in general the same, as can be seen in the one dimensional example in Figure 3.2. The curves represent the cubic Taylor polynomial approximations and the straight lines are the gradients at the shared edge for each control volume. The flux through this edge is calculated by taking the average of the two gradients, weighted by the control volumes' areas: VT = W o + W , A + A\ 0 If this is done along the shared edge in the Figure 3.2 example, the result would be as shown in Figure 3.3. With the gradient along the edge known, and the unit normal to the shared edge given by the geometry of the mesh, the flux integral can be evaluated. The flux integral is evaluated to fourth-order accuracy using Gauss quadrature. Ob- 20 CHAPTER 3. HIGH-ORDER FLUX INTEGRALS 0 1 Figure 3.2: Cell Gradients Figure 3.3: Cell Average Gradients taining fourth-order accuracy implies that a cubic polynomial must be integrated exactly. To do this using Gauss quadrature requires two Gauss points for every straight edge to be integrated. For the shared edge between the cell-centred control volumes in Figure 3.1 the integral is broken up into two equal length segments, each with a specified Gauss point, normal and weight. The integral can then be evaluated as: Flux = (VT (x ) • n ) W + (VT (x ) • n ) W gl where x i and x g g2 gl gl g2 g2 (3.4) g2 are the Gauss point locations along the shared edge and W \ and W g g2 are the corresponding Gauss weights. For fourth-order accuracy, the Gauss points are given by equations 3.5 and 3.6. The Gauss weights are simply half the length of the edge. For second- and third-order accuracy only one Gauss point is needed, located at the midpoint of the edge with a corresponding Gauss weight equal to the length of the edge. The normal is the unit normal of the edge, CHAPTER 3. HIGH-ORDER FLUX 21 INTEGRALS given in equation 3.7. Xb x Xgl flgl + xa t , Kb 2 gl b Xb 2~~ g2 = (3.5) 2\/3 x + x,a h "^a X (3-6) a 2y/3 (Va - Vb, X a — X - X) a b (3.7) X b Thefluxis then subtracted from thefluxintegral for control volume 0 and added to the flux integral for control volume 1 (following the direction of the normal). This process is repeated for each internal edge of the mesh. For the boundary edges thefluxintegration process is the same except the gradient does not need to be averaged, and thefluxesonly adjust the boundary control volumes' flux integral. The method for boundary condition enforcement is discussed in Chapter 4. 3.3 Flux Integration: Vertex-Centred Control Volumes The only difference when applying this scheme to vertex-centred control volumes is the number and location of the Gauss points. Figure 3.4 shows an example of a shared edge between two control volumes with the normals and Gauss points on the edge for second-, third- and fourth-order accuracy. The formulae for the Gauss points, weights and normals, CHAPTER 3. HIGH-ORDER FLUX 22 INTEGRALS for the different orders of accuracy are given here: \X - X \ a L \Xb - x \ 2 c ^ 2 -order nd b _|_ ^2 Xg+Xh x„ Li + L \x c (.Vc Vat X a rd X) c w„ g 3 -order 2 x\ a n b+Xc X Xg + X xi b 2 g X g2 W2 9 (ife-yaiiCa-a?t) flgl (Ve -Vb,Xb- fl X) c g2 4 -order th x\ g 2 SCa + 2\/3 + X b 2 aj + 6 x 2 ^ x ! c 2 g3 #6 + <C x± 2v/3 ^6 - x C c g 2x73 = *V fl gi - a?c W^2 = = = = %2 = - ™<?4 = Li (3.21) 2 ^2 (3.22) 2 (l/6 -Va,X -Xb) a In {Vc -Vb,x b (3.23) (3.24) CHAPTER 3. HIGH-ORDER FLUX INTEGRALS Figure 3.4: Vertex-Centred Gauss Points 23 Chapter 4 Boundary Conditions For the scheme described in Chapter 3 to be completely high-order accurate, the high-order reconstruction in boundary control volumes must be consistent with the boundary condition. That is, the Taylor approximation of the solution in the control volume must be constrained to satisfy the boundary condition. These boundary constraints are in addition to the mean value constraint for each control volume. Most of this chapter deals with the derivation of the boundary constraints. Section 4.1 describes Taylor Constraints, which combine the Taylor polynomial approximation of the solution in the control volume with the Taylor polynomial approximation of the boundary condition to create appropriate constraints for Dirichlet boundary conditions, Neumann boundary conditions and a linear combination of Dirichlet and Neumann boundary conditions. Unfortunately there are a number of drawbacks to the Taylor Constraints approach. Section 4.2 describes a different and simpler way of creating the boundary constraints called Gauss Constraints. Whereas the Taylor Constraints method enforces the boundary con24 CHAPTER 4. BOUNDARY CONDITIONS 25 dition all along the boundary edge of the control volume, the Gauss Constraints method enforces it only at the Gauss points where the solution gradient is needed for the Gauss quadrature flux integral calculation. This method is also derived for Dirichlet boundary conditions, Neumann boundary conditions and a linear combination of Dirichlet and Neumann boundary conditions. Finally, the last section describes how the constraints are added to the least-squares matrix, turning it into a constrained least-squares problem. 4.1 Taylor Constraints The Taylor Constraints method was the first attempt at enforcing the boundary conditions in the reconstruction over the control volume. The constraints are created using a Taylor polynomial to approximate the boundary condition. This Taylor polynomial is matched term by term to the Taylor polynomial approximation of the solution in the control volume evaluated along the boundary edge, giving a set of constraints on the unknown coefficients of the Taylor polynomial approximation of the solution. The following derivations are for an arbitrary straight boundary edge, where the boundary condition is a function of the variable s along the edge. 4.1.1 Dirichlet Boundary Conditions The Dirichlet boundary condition is given by T (s) = f (s), where / is C . The boundary 3 function / is replaced with a cubic Taylor polynomial about the midpoint of the edge s : m /OO « fL+% (s - s ) + m (s - s y ds m 2 af A + ds 3 {s - s ) m (4.1) CHAPTER 4. BOUNDARY CONDITIONS 26 Now a similar approximation is needed for the solution function, so the terms in the two equations can be matched to form the boundary constraints. The Taylor expansion of the solution must be in terms of Taylor coefficients at the reconstruction reference point of the control volume, because these are the terms solved for during reconstruction. The fourth-order accurate Taylor polynomial approximation of the solution about the reference point of the control volume is given again in equation 4.2. The idea is to rewrite this equation, evaluated along the boundary, into a Taylor polynomial about the midpoint of the boundary edge, while keeping the desired derivatives at the reference point of the control volume. OT T{x,y) (x - X ) + c 0T dT dy c dr 2 dxdy dT dx dy 3 (x - x ) (y - y ) + c 2 dy 2 (x-x ) (y-y ) c c xy c ,.a»T dx + (x - x f c 3 dT dxdy (x - x ) (y - y f 3 | 2 crT dy 2 c 2 (x - 2 (y-y ) 2 c dT dx 2 (y - y ) + c 2 (y - y f c + + (4.2) c 3 To evaluate this polynomial along the boundary, x and y in equation 4.2 are replaced by: x{s) = \m X s) y m L + (s - s )D y(s) = y\ + (s - Sm) ( - = = x vL + (s m s ) m D (4.3) (4.4) where (x ,y ) and (xb,y ) are the endpoints of the boundary edge in Cartesian coordinates a a b CHAPTER 4. BOUNDARY CONDITIONS 27 and s and s& are the endpoints in the boundary coordinate. After considerable manipulation a CHAPTER 4. BOUNDARY CONDITIONS the solution on the boundary can be written as: T(s) (x dT dx (x -x ) dT dx (X -X f 2 c (Vm ~ Vc) + , dT + dxdy (x - x ) (y dT (x X) dx dy 2 m dT dy x ) -\- m 2 c m 2 3 m e - y) + m c 3 c | 3 m fym c yc) 2 (X - X ) {Vm m ~ Vcf e &T dy + dxdy 2 {ym - Vcf &T dy 2 + + 3 dT dT D + D + dx dy x y dT dx dT (D (x - x ) + D (y -y )) + dxdy (X ~ X ) dT 2 2 y 2 dT dy 2 m e D (y y m ~yc)+g^ m m e 2 3 2 x C + dT dx dy Dx {x - x ) (y - y ) + Dy ^ &T dxdy D (x - x ) (y - y ) + D ^ 3 m 2 y 2 dT dy 3 D, 3 c m (ym m c - yy m c . Vm Vc x |+ ^ |+ (s - s ) + c dT dxdy dT 2 2 X m Dl H n 2 dx ^ Xm c x 2 ; + D d T dT D (x - x ) + (2D D (x - x ) + D (y - y )) + dx dx dy d T (2D D (y - y ) + D (x - x )) + dxdy 3 3 2 2 m 3 c x 2 y m 3 2 x 2 dT dy y m c (s - 3 D (y - y ) y 3 m y sf m c m c + dT D + SD D + dx dy &T dx 3 3 2 x 3 y 2 d T dxdy c 3 d T\ 3 SD D V + 2 2 X d y 3 D 3 (s - s f m 6 e X m c CHAPTER 4. BOUNDARY CONDITIONS 29 Now the boundary constraints can be determined by matching up the terms in equation 4.1 with the terms in equation 4.5. The full Dirichlet Taylor constraints for any straight boundary are: dT f\r , \x dx l c dT _ dT . dT x ) -\dy m dT 2 2 (x - x ) (y m dxdy c dT {x -x f dx 6 3 m dT (x 3 - y) + m dT t c c 2/c) dx dy c m - yf d^T c (y + + c (jjm X) m dxdy ds - yy m {x + - yy m (4.6) c dy 2 dT (y 2 - x ) (y m x) m c 3 c 3 (x 2 c 3 dT dy r, dx dT dT 2 2 dxdy (D (x Y dT c D (y y +D - x) m X (y -y )) m m - y) + c + -Q-£ c 3 D. dx 3 (x dT •+ x) m c ^ ^ J+ m 3 D (x m - x ) (y m - y) + D D (x y m - x ) (y m - y) + D ^ m - yf x dx dy 2 dT c c y 3 dxdy 2 dT 3 dy df ds 2 r»2 dx D 2 . dr 3 x (4.7) dT — 2 dT 2 2D D x y+ 3 A + y - y )) + (2D D m - y ) + D [x dT - x )) + — c m 2 x 3 3 c c 3 dxdy dx m 3 - x ) + D\ (y 3 3 Dl (x - x ) + dx m y dT ds dT dxdy (2D D (x x 2 dT + Vc x c H dx dy d*f ^ Vm c 2 dT 2 2 (y D, 3 c (y y dT 2 c y 3 H D 3 x dx dy c dT 3 SD D 2 x 2 m y + 3 n r>2 , dxdy 2 D (y 2 y m D - y) c (4.8) (4.9) The constraints shown above are greatly simplified if the boundary edge is parallel to either the x- or y-axis. For example, if the boundary edge were parallel to the £-axis, yt — y a CHAPTER 4. BOUNDARY CONDITIONS 30 would be zero, meaning that D = 0, and the resulting boundary constraints are: y f\m (x = X ) -\- m dT c (x -x ) 2 , 2 m c dx 2 dT [x -x f 2 3 m dT dT , 2 (x m c aT 3 - y) + m c X) (jjm c {x 3 m - x ) (y c m - yf <PT c dxdy dy yc) dT (l/m - Vcf dT 2 x - Vc? + + (4.10) 3 D + dx m dy 2 dT (y 2 2 dT df_ ds - x )(y + dx dy dx 3 dT 2 dxdy c (Vm ~ Vc) + dy 2 dx dxdy | 2 dT Dx {Vm - Vc) + 3 dx 3 D, dxdy 2 2 dT 2 ds dx d*f dT ds dx 2 3 4.1.2 2 x 2 dT df D (x dx dy 3 D . x - x ) (y c - y) + m c (l/m - Vcf dT (4.11) dT 3 3 H m Dl (x - x ) + m dx 3 c dx dy 2 D (y 2 X m - y) e (4.12) 3 3 (4.13) Dl Neumann Boundary Conditions The Neumann boundary condition = g( ), where g is C , is dealt with in much the s 3 same manner as the Dirichlet boundary condition. However, a little extra work is required in dealing with the normal derivative. Applying the definition of the directional derivative, with the knowledge that x = x (s), the Neumann boundary condition is rewritten as: dT ^ ^ dT „ + - n , = s( S ) (4.14) The boundary function g is written using the quadratic Taylor polynomial about the CHAPTER 4. BOUNDARY 31 CONDITIONS midpoint of the boundary edge: d + 9 {») dg (s - 2 A ds (s -s ) + m sf (4.15) m ds 2 This is done so the terms can be compared with the terms in the Taylor polynomial approximation to the left-hand side of equation 4.14, which is also a quadratic function. The spatial derivatives are approximated with a Taylor polynomial. The Taylor polynomial is obtained by taking the spatial partial derivatives of equation 4.2. The derivative approximations are: dTbdvy (x, y) ^ dx dT dT 2 dT + dx c9 T dx 2 c bdry (x, y) ^ dT (x - x ) (y - y ) + c dT c dT + dT (x - x ) + c dxdy dxdy 2 (x - x ) (y - y ) + c c xy c + 3 (y - y f (4.16) c dxdy 2 dT dT 3 (y - y ) + dy c 2 dT 3 3 (x - dx c 2 2 dy {y - y ) + 3 dx dy dT 3 dxdy 3 2 dT 2 (x - x ) + (y - dy 3 yf c dx dy 2 (x - xy c + (4.17) To convert this to a polynomial containing (s — s ) terms, the substitutions shown in m equations 4.3 and 4.4 are used in equations 4.16 and 4.17. After considerable manipulation CHAPTER 4. BOUNDARY 32 CONDITIONS the solution on the boundary can be written as: dT{x,y) dx dT dT dT 2 + dx 2 {x x ) -\- m dx 2 c d*T dx dT D + x dx 2 - y)+ c m (Vm - VcY dxdy 2 + dT 2 2 c m 2 dT - x ) (y (x dx dy 3 ("Um ~ Vc) + dxdy 3 Dy + dxdy dx 3 aT 3 &T {{x m dx dy 2 dT - x ) D + (y c dx dy dy c x dT . H 2 D x dx dy 2 2 D x D y dxdy + dT dT 2 + (Xm dxdy (x -x f dx dy 2 m 2 dT c 2 dxdy dT (s - y (s - s ) + m sy (4.18) m 2 dT 3 - Vc) D m 2 Dl 2 (y dxdy 3 3 n 3 dT - y) D) + m dT 3 dT(x,y) y D + x X ) ~\~ dy c (y - y ) + m 2 c dT dT 3 3 + (x ' dxdy 2 &T dy - x ) (y c (y m - y) + m - yy c c + &T Dy + 2 m dx dy 2 {XjYl ^c) dT ^X 3 3 dxdy 2 ((x m ar 3 dx dy - x ) D + (y r.2 D A 2 . x c y m - y) D) + c x &T dT 3 dy 3 D. (y m (s - - y) D c s) y (s - s ) + m 2 m (4.19) dxdy 2 Now equations 4.18 and 4.19 are substituted back into 4.14. The boundary constraints are then determined by matching up the terms with those in 4.15. The full Neumann Taylor CHAPTER 4. BOUNDARY 33 CONDITIONS constraints for any straight boundary are: dT 9\ m = dT + dx dT x ) -\- m dx 2 c (Vm ~ Vc) + dxdy x) m c + dx 3 dT 3 (x - 2 dT x ) (y m dx dy c - y) + m dT c 2 dT 2 x)+ (x + dxdy dT fl dxdy 2 dy c m dy dT 3 (x - x ) (y m dxdy 2 dT c 2 - y) + m c dT D + x 2 dT dT + (x x) m c dx dy 2 (Vm - VcY + (4.20) dy 3 dT 2 dx x 3 {Vm ~ Vc) + 2 3 3 D„ dxdy ' dx y 3 3 dT 3 ((x m dx dy 2 dT - x ) D + (y c y m dT 2 c x {Vm ~ Vc) D dxdy y 2 h + x 3 y dT - y) D) + dT 2 dxdy {Xm dx dy + 2 3 x ) D -|c dT x 3 ((x m dxdy 2 &9_ dT r»2 dx dx D H 3 2 {x 3 2 (x 3 dj_ dx dT dT 2 3 x ; dT 3 dx dy 2 - x ) D + (y c y m - y) D) + c x dT {y m ~ Vc) D hy y (4.21) 3 dT 3 dx dy crT Dl + dxdy Dl hx + 2 n dT 3 2 2 2D D x y + dy 3 Dl (4.22) As for the Dirichlet case, the constraints shown above are greatly simplified if the boundary edge is parallel to either the x- or y-axis. For example, if the boundary edge were parallel to the x-axis, y — y would be zero, meaning that D = 0, and the resulting boundary conb a y CHAPTER 4. BOUNDARY 34 CONDITIONS straints are as shown in equations 4.23 to 4.25. dT 9\m = dT 2 {x + dx dx X) m 2 dT c 3 - x ) (y m dx dy 2 dT c dT m - y) + c 2 dy c (x x) + m m dT c dxdy - x ) (y c dT dT 2 dT D ° + d x (y dy dT - y) + c dg dx 2 2 4.1.3 dT dx dT (x {Vm - Vcf dy 3 dT x) m c dx dy + n,. (4.23) 2 3 {x x) D + m 3 c x dx dy 2 (y dT m - Vc) D x |h + x 3 D + x Dl (x m dx dy 2 3 3 - Vc) + m 2 2 dxdy + dx 3 3 2 dx m c (Vm ~ Vcf 3 (x X) m dxdy 3 2 (x 3 2 2 + dxdy dT dx (Vm ~ Vc) + dT 3 (x dr 3 dT -jdxdy 2 aT 3 h x + dx dy 2 Xc) Dx "h dxdy (y m - Vc) D.. h y (4.24) (4.25) Dl Dirichlet & Neumann Combined Boundary Conditions The Taylor constraint equations for the Dirichlet and Neumann boundary conditions can be coupled together (combining appropriate terms) for a general boundary condition of the form aT(s) + B^^- = h(s), where a and d are constants and h is C , to create new 3 boundary constraint equations for the least-squares reconstruction. The resulting boundary CHAPTER 4. BOUNDARY CONDITIONS constraint equations are: \m h = dT ( ) \c + ( ( m ~ X ) + Bfl ) a T a X a {x x) m fi(x - a e m m r 6 x + dx 2 + fi (y - Vc) n + P (x c x m dT + fi (Vm ~ Vc) fly x) c c m - x) h) c y cPT + dxdy 2 {Vm ~ Vcf {x dT x) h m (a (x - x ) (y -y ) a y 2 + c m dT + (a {Vm ~ Vc) + fin ) -Q-; + dx x c — + a "" {x fi m - x) c 2 + dy 2 „ \ d T n + 1 3a; x 3 \2 a — — ^ ^ (x Oi m m — — — X ) (Vm ivm-vcf a V c + P(x - Vc) c x , o {Vm ~ Vc) + , (vm-y )\ + fi R x ) (ym - Vc) n m fi c n, y 3 fi^ m , / \r n + fi(x - x ) (y 0 x \ d*T dy + m c m ^ h y \ ~ - y) n c y CHAPTER 4. BOUNDARY dh ds dT (aD ) dr + (a (x - x ) D + + {<*Dy) dy dx x 36 CONDITIONS m c d?T dx 3D h ) x x x + 2 dT 2 (a ((x - x ) D + (y - y ) D ) + 0D h m c y m c x y dT (a (y - y ) D + BDyhy) dy + x 3D h ) x + dxdy y 2 m a c (x y x) m + 2 <9 T ' dx 3 D + 3 (x - x ) D h, c x m c x x 3 2 a [x - x ) (y - y ) D + a^ ^ Xm m c m c D + X x y dT dx dy 3 3 ((x - x ) D + (y - y ) D ) n + 3 (x - x ) D h m c y m c x x m a (x - x ) (y - y ) D + a ^ c m c x y + 2 ^ D + Vm m c Vc y x dT 3 3 {y - y ) D h m c (y -y ) y + 3 {{x - x ) D + (y - y ) D ) h x m c dh ds m c = 2 c x + dxdy y 2 3 D + 3{y - n y dT (-Dl) dx y dT dxdy 2 (2aD D ) x y y 3 dT + (-D ) dy 2 + (4.27) n m 2 2 m , . \ dT Vc) D n dy 2 a y 2 + 2 y dT dx 2 3 (a (x - x ) Dl + 3D h ) 2 m c x (2a (x - x ) D D m c x x + a(y - y + 3 y ) D + 23D D h 2 m c X x y + 3D n ) 2 x x &T + y d T 3 (2a (y - y ) D D m c x + a(x - y x ) Dl + 3D h 2 m c dy (ot (y - Vc) D\ + 3D h ) d T (-Dl) ^ ^ d x W y dx y x + 23D D h ) x y + -^-^ y (4.28) 3 d*h ds 3 2 = m y 3 + D y dT 3 aT 3 3 + ^ D ^ d x W 2 + (4.29) ^ d y - The constraints shown above are greatly simplified if the boundary edge is parallel to either the x- or y-axis. For example, if the boundary edge were parallel to the x-axis, y — y b a CHAPTER 4. BOUNDARY CONDITIONS 37 would be zero, meaning that D — 0, and the resulting boundary constraints are: y h\ m = dT dx (a) T\ + (a (x - x ) + fin ) m a (x x) m c c x -\-fi{x c m &T dx x) n m dT + (a (y - y ) + fifty) dy x + c + 2 dT (a (x - x ) (y -y ) +fi(y - y ) h + (3 (x - x ) h ) dxdy 2 m c m a m o m c 3 (x - x f „ \ d T Vc) + fi(x - x ) (y - y ) h + fi 2 ~ } dx dy 3 m (x - x ) (y - y ) c m (ym-yc? oi 6 x a 1 x m (aD ) dT dx (x + c hp m m c m m c m c m (4.30) 3 c x x m c x x &T dT dx + (a (y - y ) D ) dxdy c x m n n x dT dx m c x y x x dT dx dy m c x y dx dy 2 + (4.31) 2 2 c + 3 + (a {x - x ) Dl + fiD h ) m &T + fi(x - x ) D h x dT dxdy 2 2 x + 3 c c m 2 &T dx x (y - Vcf , . \ a — D + fi(y - y ) D n 2 m + 2 3 D + fi(x - X ) D h c dT dxdy 3 y 2 c + 2 \ dT n, dy m m 2 Uy c a (x - x ) (y - y ) D +fi(y -y ) D h 2 x x c x dh ds c Jy - y f fi -; n + fi{x - x ) (y - y ) h , (y -y )\ 0 c + (a (x - x ) D + fiD h ) x) m + y %c) - \ d T + n dx r 2 = c + y (%m %c) , a ip^rn +' fi2 6 m m t (%m ^c) {.Vm dh ds x 2 c OJ- a c , ,Ad T + fi(y - y ) n I -Q-2 (ym-y ) 2 d c x d*T dx 3 + 3 (a (y - y ) D\ + fiD h ) 2 m ds 3 e d^T dx 3 x y 2 (4.32) (4.33) CHAPTER 4.1.4 4. BOUNDARY 38 CONDITIONS Drawbacks There are several drawbacks with the Taylor Constraints method for boundary condition enforcement. • The total number of constraints can be excessive. Table 4.1 compares the orders of accuracy, the total number of constraints (assuming Dirichlet boundary conditions and including the mean constraint) for control volumes with one and two boundary edges and the number of unknowns being solved in the reconstruction. Several of the values in the tables are stressed in bold to indicate there are a greater number of constraints than unknown variables. The least-squares problem in these cases is over-constrained. There are situations, emphasized in italics, where there are the same number of constraints, or one less, than there are unknowns, which makes the reconstruction less reliable, since there is only one degree of freedom for the leastsquares algorithm to calculate the unknown value accurately. A t second-order this is not a problem, since the second-order least-squares matrix will never degenerate with two constraints and one integral constant. Table 4.1: Constraint Comparison for Control Volumes with Taylor Constraints Order of Accuracy 2 3 4 Cell Boundary Edges One Two 3 4 5 Vertex Boundary Edges Two 5 7 5 7 9 9 Unknowns 3 6 10 Another limitation is that the boundary condition equation must be C (to assume it 3 CHAPTER 4. BOUNDARY CONDITIONS 39 can be approximated with a cubic Taylor polynomial for fourth-order accuracy). This is the case for any scheme to be fourth-order accurate. • The derived Taylor constraints were for arbitrary straight boundary edges only. While it is possible in principle to derive the Taylor constraints for curved boundary edges, by replacing equations 4.3 and 4.4 by approximating x (s) with a cubic Taylor polynomial in s, the method for doing so is quite complex and, considering the existing problems with the scheme, the derivation was not done. Because of these problems, a different approach was developed for creating the boundary constraints to add to the least-squares reconstruction problem. 4.2 Gauss Constraints The Gauss Constraints method was the second attempt at generating boundary constraints. This method enforces the boundary condition at the Gauss points used in the flux integral, instead of enforcing it along the entire boundary edge of the control volume. The reconstruction polynomial is evaluated at the Gauss points and is set equal to the boundary condition at that point. This gives one constraint on the coefficients of the reconstruction polynomial per Gauss point. These boundary constraints may be applied to any boundary edge, straight or curved, where the Gauss point locations and normals are known. The Gauss Constraints method is simpler to implement than the Taylor Constraints method, and does not suffer from the same drawbacks. CHAPTER 4.2.1 4. BOUNDARY 40 CONDITIONS Dirichlet Boundary Conditions The solution in the boundary control volume T (x,y) is approximated with a cubic Taylor polynomial about the reference point of the control volume. This polynomial is then evaluated at each Gauss point along the boundary edge of the control volume and equated to the known boundary temperature / (xbdry,Vbdry)- The result is one boundary constraint per Gauss point along the boundary edge: f(x ,y ) g 9 = T| + c dT dx , dT dx {Xg-X f &T {Xg~X f 2 _ 3T X) + [Xg ^ c , c + 2 2 dx g dxdy c g {x - x ) (y - y ) dT dy {Vg ~ Vcf 2 d'T c g c g c 2 ar 3 (x - x ) (y - y ) g c g c dxdy dT 3 dy 2 (Vg ~ Vcf + + 3 2 4.2.2 (x - x ) (y - y ) + + dx dy 3 g 2 dT , C (y ~ Vc) + (4.34) Neumann Boundary Conditions The Gauss constraint method for the Neumann boundary condition = g (s) is almost as simple as for the Dirichlet boundary condition. The Gauss constraint forces the normal derivative at the Gauss point on the boundary to the specified condition. Applying the definition of the directional derivative, with the knowledge that x = x(s), the Neumann boundary condition is rewritten as: dT dT n + —h dx dy x y = g [x (s), y (s)) (4.35) The spatial partial derivatives of T can be determined by taking the partial derivatives of the Taylor polynomial approximation of T. This gives the following Taylor polynomial CHAPTER 4. BOUNDARY 41 CONDITIONS approximation to the partial derivatives: dT dx dT dT dT 2 + dx dx c 2 aT [x - x ) (y - y ) + c 2 dT dy dT dT dxdy c 2 + dy dT c 3 c dxdy 2 c + 2 d*T (y - Vc) (4.36) dxdy 2 dT dT (y - y ) + dx dy dy dT xy c 2 c 2 c (x - 3 3 (x - x ) (y - y ) + xy 3 c 2 (x - x ) + (x - 3 (y - y ) + dxdy 3 dx dy dT dx 2 (x - x ) + + (y - y ) 2 (4.37) c dy 3 Substituting equations 4.36, 4.37 into equation 4.35, the full Gauss constraint can be written as follows, for any given Gauss point 9(x ,y ) 9 9 dT dT (x ,y ): g dT 2 2 (x + dx dx x) + g 2 dT g c (y ~ Vc) + dxdy g dT 3 3 dx dy 2 (x - x ) [y - y ) + g c dT dT g c 2 dxdy dy x ) "f- g c X) c dx (y - Vcf 9 dxdy 2 dT 2 (x (Xg 3 dy 2 ar aT 3 (y - y) + g c (9 X c) X dx dy 2 3 dT 3 dxdy 2 (x - x ) (y - y ) + g c g c dy 3 (y 9 - yet (4.38) CHAPTER 4. BOUNDARY CONDITIONS 42 Collecting like terms, the general Neumann Gauss constraint becomes: 9(x ,y ) 9 9 = (%) dT dT (fl ) + dx y dT 2 + {n (x - x )) x dy g c dT + dx 2 2 (n (y - y ) + n (x - x )) x g c y (x -x ) \ n c dT dx 2 g g 2 + {% (y - yc)) dxdy 9 (x -x ) \ dT dx dy 2 h (x - x ) (y - y ) + n c n ^—z~x - (y -y ) \ 2 9 c 2 + + 3 g dy 3 c x x dT g c g 3 c y 2 + h (x - x ) (y - y ) y g c g c SrT dxdy 2 + + d*T dy 3 4.2.3 (4.39) Dirichlet & Neumann Combined Boundary Conditions The Gauss constraint equations for the Dirichlet and Neumann boundary conditions can be coupled together (combining appropriate terms) for a general boundary condition of the form aT (s)+0^^^ = h (s), where a and 3 are constants, to create new boundary constraint equations for the least-squares reconstruction. The resulting boundary constraint equation CHAPTER 4. BOUNDARY CONDITIONS 43 for the Gauss point (x ,y ), along the boundary, is: g h{x ,y ) g g g 71 + ^ = a (Xg dx l c dT dT X ) ~f~ &T {Vg ~ Vc) + dy c dT 2 (x - x ) (y - y ) + g dxdy &T dx dy c g (x - x f g c (y - y ) c g c 9 , dT + dx dy (x - x ) (y - 3 g c g 2 dT 3 (Vg ~ Vc) X) c yf c + + + 6 dy 3 + 3 dT dxdy 2 C (Xg 3 c 2 + X) dx (y -y f 2 (Xg 2 dT (xg - x )) dx 2 P M dx- + (fix y dy~ + ) {n c 2 + dT dT 2 (fl (V ~ Vc) + fly X 9 . (x -x ) \ 2 c ffr dx 2 g (Xg ~ X )) c U x 3 (Xg-X ) \ 2 h (x - x ) (y - y ) + n g c (Vg-Vcf & n. (y -y ) \ 2 n, 9 c g c + (hy (Vg ~ Vc)) dxdy c y c FT dy 3 FT dx dy + FT dxdy + 2 (Xg ~ X ) (Vg ~ Vc) + fly + + . x 2 2 (4.40) CHAPTER 4. BOUNDARY CONDITIONS 44 Collecting like terms results in: h(x ,y ) g g = (a)T\ +(a(x -x ) a c g (Xg X) dT + 0n ) c x dx g dT c y x c + dx 2 (a (x - x ) (y - y ) + 0 (y - y ) n + 0 {x - x ) n ) g c (Vg ~ Vcf a- c g + P(y - (x -X f , c g - x g dT Vc) n c y 2 dxdy + Xg c g dT S 3 y) - + dx " 2 x y\2 (y 2 \ c + dy y ( -X )\ 6 (x c dT 2 g g + 2 + 0 (xg - x ) h C g dT dy + (a {y - y ) + 0h ) , c . (x - x f c g „ \ dT 3 dx dy 2 Jx g x )(y c 9 f yc + p ( l U ^ ^ + 0 {Xg _ X c ) { y g _ y c ) n y &T dxdy 2 y 4.2.4 dy 3 + + (4.41) Constraint Comparison It was previously shown there are too many Taylor constraints, but what about the number of Gauss constraints? Table 4.2 shows there are enough degrees of freedom in each case to be able to perform the least-squares reconstruction reliably. At second-order accuracy, there are still situations where there is only one or no degrees of freedom for the least-squares reconstruction. Fortunately this is not a problem, since the second-order least-squares matrix will never degenerate with two constraints and one integral constant. Figure 4.1 shows, on a partial mesh, the two ways cell-centred control volumes can appear in the corners. The upper left shows an example of a control volume with two CHAPTER 4. BOUNDARY CONDITIONS 45 Table 4.2: Constraint Comparison for Control Volumes with Gauss Constraints Order of Accuracy 2 3 4 Cell Boundary Edges One Two 2 2 3 3 3 5 Figure 4.1: Cell-centred control volume's Gauss point locations Vertex Boundary Edges Two Unknowns 3 3 5 3 6 10 Figure 4.2: Vertex-centred control volume's Gauss point locations different boundary edges. This gives it a total of four Gauss points. However the upper right control volume shows how it is possible to have only two Gauss points, the same number as any regular cell-centred boundary control volume. Figure 4.2 shows, on the same partial mesh, the two different ways vertex-centred control volumes appear in corners. In each case the control volume has two different boundary edges. For regular vertex-centred control volumes, this is also the case, as shown by the hollow circle Gauss points. So no matter the location of a vertex-centred control volume, the number of Gauss points remains constant. CHAPTER 4. BOUNDARY CONDITIONS 46 4.3 Constrained Least-Squares Reconstruction Method This section describes how the constraints are applied to the least-squares reconstruction problem and the method by which the problem is solved. The first step is to add the previously derived constraints to the top of the least-squares matrix. With the constraints in place the reconstruction is a constrained least-squares problem. An example of this matrix is shown in Figure 4.3. The constraint equations must be satisfied exactly, not just in a least-squares sense. To accomplish this, Gaussian elimination is performed on as many rows as the total number of constraints m (including the mean constraint). The method is described here for each constraint i = [l..m]. 1. Determine the largest constraint coefficient in row i and divide that row through by it, including the right-hand side. 2. Exchange the column with this largest constraint coefficient with column i (column pivoting). This has the effect of exchanging corresponding rows in the solution vector. 3. The new column i, now containing a 1 in row i, is then reduced using Gaussian Elimination from row i + 1 onwards. The result is a partially reduced least-squares matrix. An example of this is shown in Figure 4.4. The least-squares algorithm, by use of Householder transformations, is then applied to the matrix, which computes the unknown Taylor polynomial coefficients for the approximation to the solution for the control volume, with proper constraint enforcement. CHAPTER 4. BOUNDARY CONDITIONS 47 The least-squares algorithm by use of Normal equations can also be used with a little more work. The entire constrained matrix, shown in Figure 4.4, is not used. Instead, the least-squares algorithm by use of Normal equations solves the matrix given by y\^ to y ^Q n for the corresponding unknowns. The solved unknowns are then used to solve the constraint equations for the remaining unknown coefficients. CHAPTER 4. BOUNDARY c l,l CONDITIONS Cl,3 °l,2 Cl,4 Cl,5 • • C^io \ / T\ u e 2,l C2,2 02,3 C ,4 C2,5 • •• C ,io C3,l C3,2 C3,3 C3,4 C3,5 • •• C io C 2 2 dx \c 3) dT dy c dT dx 2 2 Xl,l X2,l Xl,2 Xi,3 %2,2 ^2,3 ^1,4 Zl,5 • •• ^1,10 Z ,4 Z ,5 • •• ^2,10 2 2 dT dxdy c 2 2 '• 8PT 9y \ Xn,l %n,2 n,3 x Xn,4 x 3 Xn,5 J Figure 4.3: Example of the least-squares matrix. The c's represent the constraint equations and the x's represent the rest of the least-squares matrix 1 di, 0 1 rf ,3 0 0 1 4 rfi,3 2 rfi,5 4,io 4,5 4,10 d ,2 4,5 4,10 di, ^2,2 2 \ aT Qx 2 1 3 dT dy D 3 er I 0 0 0 2/1,2 1/1,5 0 0 0 2/2,2 • • • - ax \c ••• l/uo dT dxdy Y t 2 1/2,5 ••• 1/2,10 Y 2 • • dT By 3 \ 0 0 0 2/„,2 lfo,B ••• l/n.10 3 \y.) Figure 4.4: Example of the constrained least-squares matrix, where columns 2 and 4 were exchanged by column pivoting. The d's represent the altered constraint equations and the y's represent the rest of the altered least-squares matrix Chapter 5 Validation Thefinite-volumeflux scheme, developed in the previous chapters for two types of control volumes, was validated in several stages. The accuracy of the fluxes wasfirstvalidated for the interior control volumes, then over an entire simple square mesh with Dirichlet boundary conditions, for both boundary constraint methods. The Gauss Constraints method was shown to be superior and verification continued with more complex problems, including combined Dirichlet and Neumann boundary conditions and meshes with curved boundaries. Accuracy was assessed by examining the behaviour of the norms of the error in the flux integral with mesh refinement. Thefluxerror is calculated by subtracting the numerically calculatedfluxintegral from the exactfluxintegral. The exactfluxintegral is calculated by evaluating J V TdA 2 using Gauss quadrature over control volumes, where T is the known cv initial solution. The Gauss quadrature for the exactfluxintegral is performed to sixth-order accuracy to ensure that errors in this integral do not contaminate the results. The error of each test case is examined using the L error norms, where k = 1,2,..., oo. fc 49 CHAPTER 5. VALIDATION 50 The L i norm describes the error globally — all errors are weighted equally. As k increases, the error norm includes more local errors, until the error is strictly a local error at infinity. The L error norm, weighted by control volume area, is calculated using the formula shown k in equation 5.1, where N is the number of control volumes, Aj is the control volume area and Ei is thefluxerror, calculated as the exactfluxintegral minus the numericalfluxintegral, for the i th control volume. The error norms used during validation are the L is a measure of both global and local error) and the L ^ . The 1 ; L (which 2 error norm is simply the highest error magnitude in the mesh. ( Lfe Error Norm EAE^ 8 =1 (5.1) N E A The order of accuracy was determined by plotting the number of control volumes versus the error norm on a log-log plot. The data points are obtained by running afluxcalculation program using the same initial solution and boundary conditions over several identical geometries, meshed with a variety of coarsenesses. Examples of the coarsest meshes can be found in Appendix A. Ideally the data points on the log-log plot give a straight line whose slope represents negative one half of the order of accuracy obtained. The reason for this is explained as follows: the area of a control volume in the mesh is proportional to 1/7V, where N is the number of control volumes in the mesh. So the length scale AL is proportional to IN ^ (N~~^\^ l/y/N\ Suppose that the error is 0(Ax ). reduces to n Then the L error norm is y fe = N~%. Taking the log of both sides results in ' " ^ j 1 N , which = — f • In practice the slope of the line of best fit through the data points is used to determine the order of accuracy. CHAPTER 5. VALIDATION 51 For the rest of this thesis, the legend for the error norm plots is shown in Figure 5.1. It is also important to note that the flux error plots shown are for coarse meshes, since the finer meshes do not print well on paper due to the smallness of the control volumes. 2 -Order _ _ _ 3 -Order 4*-Order nd rd Figure 5.1: Error Norm Plot Legend 5.1 Interior Control Volume Flux Integral Validation The interior control volume fluxes were validated before attempts at computing the boundary conditions were made. The interior control volumefluxeswere obtained by running the flux integral calculation program with an initial temperature distribution of T (x, y) = sin (TVX) sin (iry) on square meshes with a domain and range of [—1,1] with all boundary conditions set to zero. Figure 5.2 shows a plot of the initial temperature distribution with cell-centred control volumes. Figure 5.3 shows a plot of the exactfluxwith vertex-centred control volumes. Table 5.1 contains the calculated error norms of the error data collected during testing of the cell-centred and vertex-centred schemes for the interior control volumes. The error norms were calculated using only the interior control volumes of the meshes. Figure 5.4 shows the L error norm plot for the interior cell-centred control volumes and Figure 5.5 2 shows the same for vertex-centred control volumes. The plots clearly show the slopes of the CHAPTERS. 1.0043e+000 8.9856e-001 7.9284e-001 6.8713e-001 5.8142e-001 4.7571e-001 3.6999e-001 2.6428e-001 1.5857e-001 VALIDATION mm. i 5 • • • • • • 3.8960e-001 3.4859e-001 3.0758e-001 26657e-001 2.2556e-001 1.8455c-001 1.4354e-001 1.0253e-001 6.1517e-002 20506e-002 -2.0506e-002 -6.1517e-002 -1.0253e-001 -1.4354e-001 • -1.8455e-001 • -2.2556e-001 • -2.6657e-001 • -3.0758e-001 -3.4859e-001 -3.8960e-001 :j!« i l l M 5.2a%e-002 -5.2856e-002 -1.5857e-001 -2.6428e-001 3.6999c- (XH ^.7571e-001 -5.8142e-001 -6.8713e-001 -7.9284e-001 -8.9856e-001 -1.0043e4O0O Figure 5.2: Initial temperature of sin(7rx) sin(7ry) on cell-centred control volumes l e +02 le+03 le+04 le+05 Figure 5.3: sin(7nr) sin(7ry) control volumes le+01 le+02 Number of Control Volumes Exact flux for on vertex-centred le+03 le+04 Figure 5.4: L error norm plots for interior cell-centred control volumes 2 Figure 5.5: L error norm plots for interior vertex-centred control volumes 2 Table 5.1: Error Norms of the Interior Control Volumes Error Norm U U Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 1.94 1.94 1.71 th 2.94 2.94 2.67 le+05 Number of Control Volumes 3.89 3.89 3.67 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.42 2.37 2.02 th 3.46 3.41 2.95 4.29 4.27 3.91 CHAPTERS. VALIDATION 53 lines increasing with order of accuracy. The slopes are almost exactly half the desired orders of accuracy, indicating the correctness of the scheme with both control volumes. 5.2 Dirichlet Boundary Conditions on a Square Mesh With the interior control volumes proven to be high-order accurate, the boundary flux calculation becomes the next target to validate high-order accuracy. The first test case run which included boundary conditions used the same mesh and test function as for the interior control volume testing, with compatible Dirichlet boundary conditions. Figure 5.2 shows a plot of the function using cell-centred control volumes and Figure 5.3 shows a plot of the exact flux used to compare the numerical results for vertex-centred control volumes. 5.2.1 Accuracy Analysis — Taylor Constraints This section analyses the results of running the flux calculation program on a square mesh using the Taylor Constraints method for generating the boundary constraints for the reconstruction. For cell-centred control volumes Figure 5.6 shows the error between the numerical flux, calculated to fourth-order accuracy, and the exact flux. The flux errors in some of the boundary control volumes are significantly larger than the interior error, known to be highorder accurate. The L error norms for each of the three orders of accuracy for the cell-centred scheme 2 are plotted in Figure 5.8. The scheme is shown to be second-order accurate by the solid CHAPTERS. VALIDATION 54 - 1.9707e-001 - 1.7633e-001 • 1.5558e-001 - 1.3484e-001 - 1.1409e~001 - 9.3349e-O02 • 7.2605e-O02 5.1861e-002 3.1116e-002 1.0372e-O02 -1.0372e-002 -3.1116e-002 -5.1861e-002 • -7.2605e-002 • -9.3349e-002 • -1.14O9e-001 • -1.3484e-001 • -1.5558eO01 • -1.7633e-001 • -1.9707e-001 Figure 5.6: Flux error for the cellcentred control volume scheme using Dirichlet Taylor constraints Figure 5.7: Flux error for the vertexcentred control volume scheme using Dirichlet Taylor constraints le+OO le-01 le-02 g le-03 p le-04 I 5S K \ s le-05 \ Ie-06 le-07 le-03 1 le+02 lc+03 le+04 le+01 le-t-02 Number of Control Volumes le+03 le+04 Number of Control Volumes Figure 5.8: L error norm plots for the cell-centred control volume scheme using Dirichlet Taylor constraints 2 Figure 5.9: L error norm plots for the vertex-centred control volume scheme using Dirichlet Taylor constraints 2 Table 5.2: Error Norms with Dirichlet Taylor Constraints Error Norm Li L Loo 2 Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.16 2.16 1.88 th 3.31 3.09 2.00 4.14 3.88 3.09 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.20 1.76 1.08 th 2.23 1.56 0.84 2.12 1.51 0.89 CHAPTERS. VALIDATION 55 straight line with a slope of -1. However the scheme does not perform as well at higher orders of accuracy, as indicated by the slightly crooked lines for third- and fourth-orders of accuracy. The high orders of accuracy do improve the solution, but not consistently, most likely because of the over-constraining problem described in the Taylor constraint method for generating boundary control volumes. For the vertex-centred control volumes, the story is even worse. Figure 5.7 shows the error between the numerical flux, calculated to fourth-order accuracy, and the exact flux. Thefigureshows that the flux error of the boundary control volumes greatly outweighs that of the high-order accurate interior control volumes. The L error norms for each of the three orders of accuracy for the vertex-centred scheme 2 are plotted in Figure 5.9. The scheme is shown to be second-order accurate, regardless of desired order of accuracy, most likely because of the over-constraining problem described in the Taylor constraint method for generating boundary constraints. Table 5.2 shows the different error norms for each order of accuracy for each control volume type. The orders of accuracy are computed from the slopes of the lines of best fit to the data points. The vertex-centred control volumes do not achieve high orders of accuracy with the Taylor constraints. 5.2.2 Accuracy Analysis — Gauss Constraints This section analyses the results of running the flux calculation program on a square mesh using the Gauss Constraints method for generating the boundary constraints to add to the CHAPTER 5. VALIDATION 50 - 7.3003e-003 - 6.5318e-003 - 5.7634efl03 - 4.9949e-003 - 4.2265eO03 - 3.4580e-003 • 26896e-003 - 1.9211e-003 1.1527e-003 3.8422e-004 -3.8422e-004 -1.1527eO03 -1.92n&-003 • -2.6896e-0a3 lav -til • -3.4580e-003 •4.2265e-003 • -4.9949iM»3 • -5.7634e-003 • -6.531SWXB • -7.3003e-003 Figure 5.10: Flux error for the cellcentred control volume scheme using Dirichlet Gauss constraints Figure 5.11: Flux error for the vertexcentred control volume scheme using Dirichlet Gauss constraints Figure 5.12: L error norm plots for the cell-centred control volume scheme using Dirichlet Gauss constraints Figure 5.13: L error norm plots for the vertex-centred control volume scheme using Dirichlet Gauss Constraints 2 2 Table 5.3: Error Norms with Dirichlet Gauss Constraints Error Norm Li L Loo 2 Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order 1.99 3.24 4.15 1.99 3.26 4.17 1.78 2.78 3.91 nd rd th Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order 2.59 3.70 4.60 3.62 2.53 4.51 2.17 3.01 3.85 nd rd th CHAPTER 5. VALIDATION 57 reconstruction. Figure 5.10 shows the error between the numericalflux,calculated to fourth-order accuracy, and the exact flux. As with the Taylor boundary constraints, thefigureshows that some of thefluxerrors in the boundary control volumes outweigh the interior error, known to be high-order accurate. However, unlike using Taylor constraints, the boundary error diminishes as rapidly as the interior error when the mesh is refined. The L error norms for 2 each of the three orders of accuracy for the cell-centred scheme are plotted in Figure 5.12. The scheme is shown to reach each of the orders of accuracy successfully. For the vertex-centred control volumes, the story is the same. Figure 5.11 shows the error between the numericalflux,calculated to fourth-order accuracy, and the exact flux. The figure shows that some of the flux errors of the boundary control volumes outweigh that of the high-order accurate interior control volumes. Again this error decreases rapidly as the mesh is refined. The L error norms for each of the three orders of accuracy for the 2 vertex-centred scheme are plotted in Figure 5.13. The scheme is again shown to reach each of the orders of accuracy successfully. The reconstruction is no longer over-constrained and is able to complete the calculations to the proper orders of accuracy. Table 5.3 shows the different error norms for each order of accuracy for each control volume type. The orders of accuracy are computed from the slopes of the lines of bestfitto the data points. Each error norm meets the target order of accuracy, although the LQO error norm is less well-behaved than the others, because of its sensitivity to local error, which is greatest along the boundary. This proves the Gauss Constraints method for generating the constraint equations works to high-order accuracy. CHAPTERS. VALIDATION 58 Since the Gauss constraints proved superior to the Taylor constraints, the Taylor constraints were discarded. The rest of the validation deals only with cases using Gauss constraints. 5.3 Neumann Boundary Conditions on a Square Mesh The function for testing the Neumann boundary conditions was T (x, y) = sin(7nr) sinh(ny). The boundary conditions were set up such that the left and right vertical walls used §^ = ±7rcos(7ra;) sinh(7ry) and the upper and lower walls used | ^ = ± 7 r sin(TTX) cosh(iry) (the alternating signs depend on the direction of the inward boundary normal). A plot of the function is shown in Figure 5.14 for vertex-centred control volumes. Since the Laplacian of the initial function is zero, the exact flux is also zero everywhere in the mesh. 1.0586el()01 9.4713e+000 8.3S71e+0O0 7.2428o f (XX) 6.1285c >000 5.0142e+000 3.9(XXkMlXX) 2.7857e+000 1.6714c l (XX) 5.5714c-001 -5.5714c-001 -1.6714c I (XX) -2.7857m (XX) 3.9000et()00 5.0142c+(XX) 6.12&5c l (XX) -7.2428e+000 -8.3571 et (TOO 9.4713c+000 1.0586c 1001 Figure 5.14: Initial function of sin(7nr) sinh(7rj/) for vertex-centred control volumes for Neumann boundary conditions Figures 5.15 and 5.16 show the flux integral error for the Cell scheme and the Vertex scheme respectively. In both cases the error along the lower and upper boundaries is rather CHAPTERS. VALIDATION 59 3.5672e-002 3.1917e-002 2.8162e-002 2.4407e-002 2.0652e-002 1.6897e-002 1.3142e-002 9.3873e-003 5.6324e-003 1.8775e-003 -1.8775e-003 -5.6324e-003 -9.3873e-003 -1.3142e-002 -1.6897e-002 -2.0652e-002 -2.4407e-002 -2.8162e-002 -3.1917e-002 -3.5672e-002 • - 5.3924e-002 48248e-002 4.2572e-002 3.6895e-002 3.1219e-002 Z5543e-002 1.9867e-002 1.4191e^»2 8.5143e-003 28381e-003 -2.8381e-003 -8.5143e-003 - -1.4191e-002 • -1.9867e-002 • -2.5543e-()02 • -3.1219e-002 • -3.6895e^XX2 • -4.2572e-002 • -4.8248e-002 • -5.3924e-002 Figure 5.15: Flux error for the cellcentred control volume scheme using Neumann Gauss constraints Figure 5.16: Flux error for the vertexcentred control volume scheme using Neumann Gauss constraints Figure 5.17: L error norm plots for the cell-centred control volume scheme using Neumann Gauss constraints Figure 5.18: L error norm plots for the vertex-centred control volume scheme using Neumann Gauss constraints 2 2 Table 5.4: Error Norms with Neumann Gauss Constraints Error Norm u L 2 LQO Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 1.93 1.92 1.73 th 3.05 3.04 2.69 4.46 4.27 3.57 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.31 2.00 1.13 th 3.48 3.31 2.83 4.38 4.12 3.39 CHAPTER 5. VALIDATION 60 large. This remains so as the mesh is refined. This is because the gradient along these walls is very large for the boundary control volumes, a value of over 36 at the peaks where the boundary error is largest. Looking at the L error norm plots in Figures 5.17 and 5.18, the 2 slopes indicate the orders of accuracy are obtained. This is confirmed by the data collected in Table 5.4. The L ^ error norms represent a measure of the local error, which is sensitive, due to the gradients at the boundary control volumes. High orders of accuracy are still achieved, even though the approximation of the very large gradients at the boundaries are locally inaccurate. 5.4 Combined Dirichlet & Neumann Boundary Conditions on a Square Mesh - 8.4135e-001 7.5278e-001 6.6422e-001 5.7566e-001 4.8710e-001 3.9853e-001 3.0997e-001 2.2141e-001 1.3284e-001 4.4281e-002 ^.4281e-002 -1.3284e-001 -2.2141e-001 -3.0997e-O01 • -3.9853e-001 • -4.8710e-001 • -5.7566e-001 • -6.6422e-001 • -7.5278e-001 • -8.4135e-001 Figure 5.19: Initial function for sin(|x) sinh(^y) vertex-centred control volumes for combined boundary conditions In an effort to see if the size of the gradients really is the problem with the Neumann CHAPTERS. VALIDATION 61 - 22762e-003 20366e-003 1.7970e-003 1.5574eO03 1.317&MM3 1.0782e-003 8.3860e-004 5.9900c-0O4 3.5940e4X)4 1.1980e-004 -1.1980e-004 -3.5940e-004 • -5.9900e-004 --8.3860e^X)4 • -1.0782e-003 • -1.3178e-003 |, .i, j a r . F i i - - , , £z.A&VS.,u.'...uwmz..„v 1 • -1.5574e-003 • -1.7970e-003 • -2.0366e-003 • -2.2762e-0a3 Figure 5.20: Flux error for the cellcentred control volume scheme using combined Gauss constraints Figure 5.21: Flux error for the vertexcentred control volume scheme using combined Gauss constraints le-01 le-01 8 le-06 V le-05 ' V ie-06 • • : .; . \ \ .• - • » s lc-07 lc-(W I ; - : le+WZ le+03 le+04 le+05 S. •. -x. lc-07 le-08 1 le+01 le+02 Number of Control Volumes Ie+03 le+04 Figure 5.22: L error norm plots for the cell-centred control volume scheme using combined Gauss constraints 2 Figure 5.23: L error norm plots for the vertex-centred control volume scheme using combined Gauss constraints 2 Table 5.5: Error Norms with Gauss Constraints Error Norm Li L 2 LQO Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.01 2.04 1.92 th 3.06 3.08 2.80 le+05 Number of Control Volumes 4.26 4.27 3.72 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.38 2.08 1.12 th 3.48 3.40 2.88 4.49 4.39 3.62 CHAPTERS. VALIDATION 62 boundary conditions, the initial problem was adjusted slightly to T (x, y) = sin(f a;) sinh(f y). This time Neumann boundary conditions were used on the upper and left walls, and the lower and right walls used a linear combination of Dirichlet and Neumann boundary conditions of the form \T + A plot of the function is shown in Figure 5.19 for vertex-centred control volumes. Since the Laplacian of the initial function is still zero, the exact flux is also zero everywhere in the mesh. Figures 5.20 and 5.21 show the flux integral error for the Cell scheme and the Vertex scheme respectively. In both cases the error along the lower and upper boundaries has decreased, indicating that the large gradients were indeed the problem in the flux integral calculation for the case in Section 5.3. Looking at the L error norm plots in Figures 5.22 2 and 5.23, the slopes indicate the orders of accuracy are obtained. This is confirmed by the data collected in Table 5.5. 5.5 Singular Point on a Square Mesh Previous cases had smooth solutions. What if that is not the case? A test case with | ^ = 0 on the upper and left walls, ^ = 1 on the right wall and T = 0 on the bottom proved to be interesting. The exact solution to this problem was computed analytically to be L.2 2 •lfk}\—~—• A plot of the function is shown in Figure 5.24. T(x,y) — 2~2 fc=odd k n s m h \ . ~ T ) Since the Laplacian of the initial solution is zero, then the exact flux integral is also zero. Figures 5.25 and 5.26 show there is a problem with the flux integral calculation in the control volumes around the bottom right corner of the mesh. The reason for this is that CHAPTERS. 63 VALIDATION 1.5973e+000 1.4292e+000 1.2611e+000 1.0929e+000 9.2478e-001 7.5664e-001 5.8850e-001 4.2035e-001 2.5221e-001 84071e-002 -8.4071e-002 -2.5221e-001 -4.2035e-001 -5.8850e-001 -7.5664e-001 -9.2478e-001 -1.0929et000 -1.261 let 000 -1.4292e4000 -1.5973et000 Figure 5.24: Initial function of singular point test case on cell-centred control volumes the corner vertex is a singular point, at which f £ ->• oo. As with the first Neumann test case in Section 5.3, there will be a large gradient (at least in the y direction) which causes problems with the flux integral calculation. As the mesh is refined, the error in the corner decreases, as indicated in the L error norm plots in Figures 5.27 and 5.28. However the 2 error at that location remains much greater than elsewhere in the mesh, resulting in, at best, second-order accuracy as shown in Table 5.6. This is a recurring problem whenever the gradients get too large for the flux schemes to handle. 5.6 Curved Mesh Boundaries With the Cell and Vertex flux integration schemes shown to work on a simple square mesh with boundaries parallel to the Cartesian axes, the next test involves curved mesh boundaries. Curved mesh boundaries do not affect the flux integral calculations of the interior control volumes: only the boundary control volumes are affected. Two different curved CHAPTER 5. VALIDATION 64 1.8220-002 1.6302e-002 1.4384e-002 1.2466e-002 1.0548e-002 8.6303e-003 6.7125e-003 4.7946e-003 28768cM»3 9.5893e-004 -9.5893e-004 -2.8768e-003 ^.7946e-003 -6.7125e-003 -8.6303e-003 -1.0548e-O02 -1.2466e-002 -1.4384e-002 -1.6302e-002 I -1.8220e-002 • • 7.9483e-003 7.1116e-003 6.2749e-003 5.4383e-003 4.6016e-003 3.7650e-003 2.9283e-003 2.0916e-003 1.2550e-0a3 4.1833e-004 ^.1833e-004 -1.2550e-003 • -2.0916e-003 •-2.9283e-003 • -3.7650e-003 • -4.6016e-003 • -5.4383e-003 • -6.2749e-003 • -7.1116e-(M3 • -7.9483e-003 Figure 5.25: Flux error for singular test case using the cell-centred control volume scheme with Gauss constraints Figure 5.26: Flux error for singular test case using the vertex-centred control volume scheme with Gauss constraints Figure 5.27: L error norm plots for the cell-centred control volume scheme using Gauss constraints Figure 5.28: L error norm plots for the vertex-centred control volume scheme using Gauss constraints 2 2 Table 5.6: Error Norms with Gauss Constraints Error Norm Li L Loo 2 Cell-Centred Control Volumes 2 -0rder 3 -Order 4 -Order nd rd 1.97 1.66 0.70 th 2.46 1.60 0.59 2.32 1.56 0.56 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.21 1.76 0.76 th 2.46 1.51 0.40 2.24 1.30 0.14 CHAPTERS. VALIDATION 65 boundaries were examined: circle boundaries and cardioid boundaries. 5.6.1 Circle Mesh Boundary- Figure 5.29: Initial function of cos (jx) sinh ( | y ) for cell-centred control volumes on a circle mesh The circle mesh boundary was defined as a unit circle centred on the origin. Examples of the two coarsest meshes are given in Figures A.4 and A.5 in Appendix A . The test function used was T (x, y) = cos (jx) sinh (f y). A plot of the initial test function is given in Figure 5.29 for cell-centred control volumes. The boundary conditions used were nn 4 8T 7T on cos (jr3?) + xsin (jr^j sinh (jryj ~ Vcos (jr%) for c ° s h (jry) > (5-2) 9 = [0, ?r) for the upper half of the circle, and the lower half is divided into two equal parts with the lower left quarter set to the exact Dirichlet condition T = cos (jri?) sinh (J^yj , for 9 = 3n 7T ' 2 (5.3) CHAPTERS. VALIDATION 66 • • • • 6.2934e-005 5.6310e-005 4.9685e-005 4.3060e-005 3.6436e-O05 29811e-005 23186e-005 1.6562e-005 9.9370e-006 3.3123e-O06 -3.3123e-006 -9.9370e-006 -1.6562e-005 • -2.3186e-005 • -2.9811e-005 • -3.6436e-005 •-4.3060e-005 • -4.9685e-005 • -5.6310e-005 • -6.2934e-005 Figure 5.30: Flux error for the cellcentred control volume scheme on a circle mesh Figure 5.31: Flux error for the vertexcentred control volume scheme on a circle mesh Number of Control Volumes Number of Control Volumes Figure 5.32: L error norm plots for the cell-centred control volume scheme on a circle mesh 2 Figure 5.33: L error norm plots for the vertex-centred control volume on a circle mesh 2 Table 5.7: Error Norms on a Circle Mesh Error Norm In L 2 Loo Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.00 2.00 1.81 th 3.14 3.20 2.87 4.60 4.56 3.96 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.18 1.71 1.09 th 3.51 3.54 2.98 4.43 4.12 3.35 CHAPTER 5. VALIDATION 67 and the lower right quarter set to the exact Neumann condition dT TT ^rsin (^z) sinh ^yj dn - ^/cos Q z ) cosh (^y) , "3?T _ \ for (9= y,27rj (5.4) T h e Laplacian of the initial function is zero, so the exact flux is also zero everywhere in the mesh. Figures 5.30 and 5.31 show the flux integral error for the Cell and Vertex schemes respectively on the coarsest mesh. T h e error for the interior control volumes is smaller than the error in the boundary control volumes. A s the mesh is refined, the error in the boundary control volumes decreases, more quickly than the interior error. Figures 5.32 and 5.33 and Table 5.7 show the orders of accuracy for the L i and L type. The L 2 2 error norms for each control volume error norms for the vertex-scheme at third-order accuracy are approaching fourth-order. However, the L ^ error norms for second- and fourth-order accuracy are somewhat lower than desired. T h e reason for this is currently unknown, but may be typical since these results were also obtained running different test cases. 5.6,2 Cardioid Mesh Boundary T h e more complex cardioid mesh boundary was defined with the polar function r = l+cos.9. T h e coarsest mesh is given in Figure A . 6 in Appendix A . W h a t makes this curved boundary mesh interesting is the cusp located at the origin. T h e mesh is refined as it approaches the cusp to accurately compute the flux integral. T(x,y) = sin ( f x ) sinh ( f y ) . T h e test function used over this mesh was A plot of the initial test function is given in Figure 5.34 for CHAPTERS. VALIDATION 68 — 8.7526e-001 — 7.83l3e-001 — 6.9099e-001 5.9886C-001 5.0673e-001 4.146Ge-001 3.2246e-001 2J033e-001 1.3820<H001 4.6066C-002 ^.6066e-002 -13820e-001 —-2J033e-001 — -3.2246e-u01 — -4.1*60e-001 5.0673e^»l —-5.9886e-001 6.9099<H)m 7.8313e^)01 8.7526e-001 Figure 5.34: Initial function of sin (fa:) sinh (fy) for vertex-centred control volumes on a cardioid mesh vertex-centred control volumes. The boundary conditions used were 4 T + dT ^ [ = Sin il ) ~ X X S i U (l )] X " S i D h y S i n (l*) ° C S h (l ) y ( 5 ' 5 ) for the upper half of the cardioid, ~4 ~L \ T+d = ~ xsin {[ s i n (i ) x (l )} x sinh ( i y ) " y s i n ii ) x c o s h for the lower half. The Laplacian of the initial function is zero, so the exact flux is also zero everywhere in the mesh. Figures 5.35 and 5.36 show the flux integral error for the Cell and Vertex schemes respectively on the coarsest mesh. The error is greatest in the boundary control volumes where the initial solution is largest. As the mesh is refined, the error in the boundary control volumes decreases. Figures 5.37 and 5.38 and Table 5.8 confirm the orders of accuracy for most of the error norms for each control volume type. However, the error norms are once again somewhat low at most orders of accuracy, for both control volume types. The L error norm 2 at fourth-order accuracy for the vertex-centred control volume is also too low. Again, the } CHAPTERS. VALIDATION 69 - 4.418&V005 - 3.9536e4»5 - 3.4885e-O05 - 3.0234e-005 - 2.5582e-005 - 2.0931e-005 - 1.6280e-005 - 1.1628e-005 6.9770e-006 2.3257C-006 -2.3257e-006 -6.9770e-006 -1.1628e-O05 —1.6280e-005 -2.0931e-005 —2.5582e-005 --3.0234<~005 --3.4885e-005 -3.9536e-005 -Jt.418«e-005 Figure 5.35: Flux error for the cell centred control volume scheme on a cardioid mesh Figure 5.36: Flux error for the vertexcentred control volume scheme on a cardioid mesh | Z le-08 le-08 le+03 te+04 1 le+02 le+03 Number of Control Volumes le+04 Number of Control Volumes Figure 5.37: L error norm plots for the cell-centred control volume scheme on a cardioid mesh 2 Figure 5.38: L error norm plots for the vertex-centred control volume scheme on a cardioid mesh 2 Table 5.8: Error Norms on a Cardioid Mesh Error Norm Li L 2 Loo Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 1.98 1.97 1.80 th 2.97 2.95 2.37 4.26 4.08 3.21 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.26 1.76 0.98 th 3.29 3.16 2.86 4.01 3.58 2.39 CHAPTER 5. VALIDATION 70 reason for this is unknown. 5.7 Summary Thefinite-volumefluxscheme was successfully validated for the interior control volumes. When Dirichlet boundary conditions were introduced, the Taylor Constraints method gave erratic results while the Gauss Constraints method provided clean, high-order accurate flux integrals. The Gauss Constraints method did provide problems at boundaries where the gradients are high — leading to second-order accurate results regardless of the desired order of accuracy. Verification continued as the scheme was executed over meshes with curved boundaries. For both the circle mesh and the more complex cardioid mesh, the scheme produced high-order accuratefluxintegrals. Chapter 6 Application The previous chapters have demonstrated a high-orderfinite-volumemethod for the Laplace equation. What remains to be shown is that the scheme can handle the viscous terms of the Navier-Stokes' equations to high-order accuracy. This chapter describes a modified advection-difFusion problem. It is a scalar model problem which has the same structure as the Navier-Stokes' momentum equations. 6.1 Modified Advection-Diffusion Problem Recall from Chapter 2 the Navier-Stokes' equations: dp dt d(pu) dt -l-V-(pu) = 0 + V • (pu <g> «) dE — + V-(u(E (6.1) VP + V •r +FB V • (r • « ) + V • ( K V T ) + P)) 71 (6.2) (6.3) CHAPTER 6. APPLICATION where (tt <g> u)^ F 72 = UiUj . o . - ^ g H - g D - B = net vector body force for the control volume i(v U + (V ) U ) T The viscous terms of the Navier-Stokes' equations (the boxed terms) are of the form V • r, where r for three dimensions is: 2du _ I (9v 3 dx 3 \dy T = 2fi 2 \dy dx) 2 dv 3dy dx) , 9iA 1 (du , 9"") 2 \dy ^~ dx) 1 (du , 9tu 2 \dz ^~ dx I (du 3 \ dx j_ / dv _ i _ cte? \ 2 \^9z "T" 9 ; dw\ ' dz ) i 1 (9v 4. <tU.\ 1 (du fete < dw\ 2 \dz ^~ dx) 2 \dz ^ (6.4) y l§w _ I (du , dv\ 3 dz 3 \dx ^~ dy) dy) The scalar model problem used to validate the claim that the viscous terms can be evaluated is: dT dT dT , 2 . dT 2 n n dx (6.5) where a, b, \i\ and [i are constants. 2 6.2 Applying the Finite-Volume Scheme The flux integral for this model problem is derived by using the same method as for the Laplace equation. Equation 6.5 is integrated over the control volume. The Divergence Theorem is then applied to obtain the flux integral shown in equation 6.6, where n is the unit vector normal to the edge of the control volume. The jfjf- term cannot be handled CHAPTER 6. APPLICATION 73 uniquely, and is broken up into two terms with the Divergence Theorem. Flux Integral = dT dT dT dT\ dx dy dy dx) (6.6) The flux integral is evaluated using the same technique as for the Laplace equation. A least-squares reconstruction is performed, including the Gauss constraints specified by the boundary conditions to give a Taylor polynomial approximation to T in each control volume. From this the value of T and V T can be determined at the Gauss points. Along control volume boundaries, the area-weighted averages of T and V T are used to calculate the flux. The flux integral for each edge segment of the control volume is calculated using Gauss quadrature and added to the control volumes' total flux integral. 6.3 Validation This section demonstrates the modified advection-diffusion flux integral can be computed to high-order accuracy on a square mesh. The numerical flux integral is compared to the exact flux integral, calculated to sixth-order accuracy by integrating equation 6.5 over the control volumes. The error norms are calculated and plotted to demonstrate that each order of accuracy can be obtained. The test function used to validate the modified advection-diffusion problem was T (x, y) = sin (irx) cos (iry). A picture of the initial function is shown in Figure 6.1 for cell-centred control volumes. The boundary conditions used were Dirichlet, T = sin (nx) cos (ny), along the top and bottom boundaries and Neumann, | ^ = sin (nx) cos (ny), along the left and right boundaries. The constants a, b, n\ and /J,2 are all equal to one. CHAPTER 6. APPLICATION 1.0043e+000 8.9856e-001 7.9284e-001 6.8713e-001 5.8142e-001 4.7571e-001 3.6999e-001 2.6428c-001 1.5857e-001 5.2856e-002 -5.2856e-002 4.5857e-001 -2.6428e-001 -3.6999e-001 -4.7571e-001 -5.8t42e-001 -6.8713e-001 -7.9284c-001 -8.9856e-001 -1.0043ef000 Figure 6.1: Initial function of sin (nx) cos (ny) for the modified advection-diffusion problem for cellcentred control volumes Figures 6.2 and 6.3 show the flux integral error for the cell- and vertex-centred control volumes respectively. In each case, the error in the boundary control volumes is comparable to that in the interior control volumes. The L error norm plots in Figures 6.4 and 6.5 show 2 the correct relationship between the error norm and the number of control volumes for each control volume type. A l l the orders of accuracy are validated by the three error norms shown in Table 6.1. CHAPTER 6. APPLICATION • 9.1174e-003 • 8.1577e-003 • 7.1979e-003 • 6.2382e-003 • 5.2785e-003 • 4.3188e~003 - 3.a590e-003 - 2.3993e-003 1.4396e-003 4.7986e-004 4.7986e-004 -1.4396e-003 -2.3993e-003 • -3.3590c-003 • -4.3188e-003 • -5.2785e-003 • -6.2382e-003 • -7.1979e-003 • -8.1577e-003 • -9.1174e-003 , Figure 6.2: Flux error for the cellcentred control volume scheme for the modified advection-diffusion problem le+03 Figure 6.3: Flux error for the vertexcentred control volume scheme for the modified advection-diffusion problem le+02 le+04 Number of Control Volumes le+03 le+04 Number of Control Volumes Figure 6.4: L error norm plots for the cell-centred control volume scheme for the modified advection-diffusion problem 2 Figure 6.5: L error norm plots for the vertex-centred control volume scheme for the modified advection-diffusion problem 2 Table 6.1: Error Norms for the Modified Advection-Diffusion Problem Error Norm Li L Loo 2 Cell-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.16 2.16 1.98 th 3.21 3.25 2.90 4.14 4.11 3.76 Vertex-Centred Control Volumes 2 -Order 3 -Order 4 -Order nd rd 2.58 2.52 2.11 tn 3.67 3.63 3.17 4.57 4.51 4.03 Chapter 7 Summary and Conclusions A fourth-order accurate diffusive flux calculation scheme has been developed which can be incorporated into any existing reconstruction-based unstructured mesh inviscid flow solver. Commonly used finite-volume methods were examined, concluding that a new method is needed to easily calculate high-order accurate flux integrals, using the same least-squares reconstruction method as with the inviscid scheme in order to avoid performing multiple time-consuming reconstructions. Using the Laplace equation as a model of the viscous terms of the Navier-Stokes' equations, a scheme was developed to calculate flux integrals to high-order accuracy by using data from a high-order accurate least-squares reconstruction. The boundary conditions are enforced during the reconstruction, which requires boundary constraints to be included in the least-squares matrix. Initially the Taylor Constraints method, which approximates the boundary condition with a Taylor polynomial and matches it term by term to the Taylor polynomial approximation of the solution evaluated along the boundary edge, was developed to derive the boundary constraints for straight boundaries. 76 CHAPTER 7. SUMMARY AND CONCLUSIONS 77 However, there were drawbacks: too many constraints made the constrained least-squares reconstruction matrix over-constrained, the boundary condition had to be C , and extend3 ing the scheme to curved boundaries was complex. The Taylor Constraints method was discarded and the simpler Gauss Constraints method was developed that eliminated the drawbacks. The Gauss Constraints method simply equates the reconstruction polynomial at the boundary Gauss points to the boundary condition at that point. The flux integration was validated by demonstrating high-order accuracy on square, circle and cardioid meshes with a combination of Dirichlet and Neumann boundary conditions. During testing two cases were revealed which cause problems with the flux integration scheme: initial solutions which have large gradients at the boundaries produce large errors locally at those boundaries, but still maintain high-order accuracy, and singular points along the boundary cause the high-order scheme to fail. The flux integration scheme was shown to be high-order accurate for a modified advection-diffusion problem of the same form as the Navier-Stokes' momentum equations. This demonstrates that the developedfluxintegral calculation scheme can be used to evaluate the viscous terms of the Navier-Stokes' equations. The rest of this chapter describes the important issues that arose during development of the code. Conclusions are then made about the scheme, describing the limitations and comparing different aspects of the two different control volume types, including memory usage and CPU time. CHAPTER 7. SUMMARY AND CONCLUSIONS 7.1 78 Issues During Development There are several areas in the code which must be implemented precisely correctly in order to achieve and verify high-order accuracy. Each of the items listed below was at some point a significant problem. They have been listed here in the hopes that others may avoid them. 1. The boundary vertices of the generated meshes must lie precisely on the geometry's boundary. The mesh file must also contain enough significant digits for the boundary vertices. Without either of these, attempts at verifying that the flux integral is highorder accurate will fail. The Gauss points on the boundary are computed from the boundary vertices. If the boundary's vertices are not accurate enough, the Gauss points will be in the incorrect place. This causes a small error in the high-order solution, but causes big problems during verification, as the computations of the initial and exact solutions are also affected. 2. Likewise, the location of the Gauss points must be computed accurately so they lie on the straight or curved boundary of the geometry. 3. Gauss quadrature for integration of the moments along the curved mesh boundaries must be performed with respect to the correct variable. The accuracy of the reconstruction, and hence flux integral calculation, absolutely depends on this at high orders of accuracy. In the code used for this thesis, the moments are integrated analytically in x and numerically in y. It is important to compute the Gauss points and weights with respect to y and not along the arc parameter s, as done for the flux integration. CHAPTER 7. SUMMARY AND CONCLUSIONS 79 During validation, the mean values of the initial solution in the boundary control volumes of the mesh need to be calculated in the same way as the moments. In either case, if the integration is done with respect to s, only second-order accuracy can be obtained. 4. During validation, the initial function and its exact flux function (the Laplace or advection-diffusion operator applied to the initial solution) must be correct with respect to each other. If they are not, the computed error norms will not give the correct order of accuracy and the plots of the flux integral error will be erratic. This was detected when the interior control volumes gave large flux integration errors when they were already validated to be correct. As part of the testing code, the analytic integration with respect to a; of both the initial function and its flux function must also be done correctly, otherwise similar problems will occur. 5. Another technicality in validating the code is to ensure the specified boundary conditions match the initial solution. If this is not done, the flux integration error will be large where the boundary condition is incorrect, possibly leading to the false conclusion that the boundary flux integration code is incorrect. Also the boundary conditions cannot have any discontinuities along the boundary, otherwise results such as the validation case in Section 5.5 will occur. CHAPTER 7.2 7. SUMMARY AND CONCLUSIONS 80 Computer Resource Usage It is difficult to make conclusions about which control volume is preferable. In addition to accuracy, memory usage and computation time need to be weighed carefully. When memory usage is considered, the cell-centred control volumes are the better choice. Memory usage was computed for two of the sample square meshes used during testing; size statistics for each mesh are listed in Table 7.1. The finer mesh is about 66 times larger than the" coarser mesh. Table 7.1: Mesh Information for Square Meshes Mesh Name squarel.mesh square4.mesh # Cells # Edges # Boundary Edges # Vertices 550 37726 857 56841 64 504 308 19116 There are two main data structures which consume memory: the coefficients of the reconstructed solution and the mesh data structure. Table 7.2 shows the approximate amount of memory in kilobytes used for each order of accuracy by each mesh. It is important to note that these can be considered an upper bound for memory usage, as the data structures used in the thesis code are by no means optimal. The reconstruction data structures increases in size as the order of accuracy increases, since the number of terms, and hence coefficients, of the Taylor polynomial approximation increases. In the mesh data structure the number of moments increases with order of accuracy, as does the number of neighbouring control volumes for each individual control volume. The cell-centred control volumes use less memory than the vertex-centred control vol- CHAPTER 7. SUMMARY AND CONCLUSIONS 81 Table 7.2: Memory Data for Square Meshes (kilobytes) Mesh 1 4 Order 2 3 4 2 3 4 Cell-Centred C V ' s Recon Data Mesh Data Total 13.2 26.4 44.0 905.4 1810.8 3018.1 94.8 123.8 164.5 5955.5 7916.7 10663.7 Vertex-Centred C V ' s Recon Data Mesh Data Total 108.0 150.2 208.5 6860.9 9727.5 13681.8 7.4 14.8 24.6 458.8 917.6 1529.3 141.3 163.0 250.9 8791.2 10248.7 15570.6 148.7 177.8 275.5 9250.0 11166.3 17099.9 umes. And since there are about twice as many cells as vertices, cell-centred control volumes are much less memory intensive per control volume. The vertex-centred control volumes use significantly more memory for fourth-order accuracy than for third-order accuracy. This is because the interior Gauss points for the cell-centred case are computed during the flux integral calculation, while the interior Gauss information for the vertex-centred case is computed during preprocessing and stored in the mesh data structure. The Vertex scheme could use less memory by computing the Gauss information on the fly like the Cell scheme, but this would increase the C P U time. / .../... / Cell 2 -Order Cell 3 -Order Cell 4 -Order Vertex 2 -Order Vertex 3 -Order Vertex 4 -Order nd «< . ld ,h ,h / ,h th 0 ^ le+05 2e+05 3e+05 Number of Control Volumes Figure 7.1: Plot of number of control volumes versus the amount of memory used for each control volume scheme. CHAPTER 7. SUMMARY AND CONCLUSIONS 82 Figure 7.1 shows the relationship between the order of accuracy and amount of memory used by each control volume scheme. The memory increase between meshes is approximately linearly proportional to the number of control volumes, about 65 times for cell-centred control volumes and 62 times for vertex-centred control volumes. Overall memory use may or may not be an issue, depending on the amount of memory in the computer and the size of the finest mesh for the problem's geometry. Although the Vertex scheme uses more memory, it executes more quickly than the Cell scheme. Figures 7.3 and 7.2 show the relationship between desired orders of accuracy and the C P U time to compute the flux integrals for the two different schemes. The times were recorded for the sinusoidal test case on the finest square mesh used in Section 5.2.2. A Pentium 233 MHz personal computer, running the Windows98 operating system, with 40 MBytes of R A M was used to perform the calculations. As with memory usage, the times can be considered an upper bound, as efficiency was not taken into consideration during implementation. However, implementation is similar between the two schemes, and thus a basic comparison can be made. In each figure the lines cross at approximately the same point. After this point, fourthorder accurate methods reach lower error norms more quickly than second- and third-order accurate methods. This strongly suggests that solving problems over coarse meshes to highorder accuracy is more efficient than solving to second-order accuracy over fine meshes. As the order of accuracy increases, the C P U time to compute the flux integral over cell-centred control volumes essentially doubles. The vertex-centred scheme's times increase CHAPTER 7. SUMMARY AND CONCLUSIONS 83 2 -Order ___ 3 -0rder - - 4 -Order nd ld ,h le-03 lc-04 le-05 le-06 i 0.1 CPU Time (seconds) Figure 7.2: Plot of L error norm versus the amount of C P U time used for cellcentred control volumes. 2 I CPU Time (seconds) Figure 7.3: Plot of L error norm versus the amount of C P U time used for vertex-centred control volumes. 2 linearly (a factor of one) with order of accuracy. The main explanation for the time difference between control volumes types is that there are fewer vertex-centred control volumes than cells in any mesh — roughly half as many. Because the most time consuming part of the calculation is the reconstruction over each control volume, having fewer control volumes to reconstruct greatly decreases the amount of C P U time needed. Because the Cell scheme is working with twice as many control volumes as the Vertex scheme, it should be more accurate on a given mesh. When the C P U time difference is a factor of two or less, the cell scheme is more efficient on a per control volume basis. One last consideration is the possibility of extending the two-dimensional schemes to three dimensions. The details for cell-centred control volumes in three dimensions are far easier than for vertex-centred control volumes because the flux integration for oddly-shaped control volumes is quite complex. CHAPTER 7. SUMMARY AND CONCLUSIONS 7.3 84 Direction for Future Work The flux integral of a modified advection-diffusion equation in a form similar to the NavierStokes' momentum equations was calculated to high-order accuracy, indicating the developed diffusive flux scheme can be incorporated into Dr. Carl F. Ollivier-Gooch's existing inviscid compressible flow solver to evaluate the viscous terms and make it a viscous compressible flow solver. This could be done by performing the constrained least-squares reconstruction and then using the approximate control volume solutions to calculate the flux integrals of each of the terms (inviscid and viscous) in the Navier-Stokes' equations in a manner similar to the modified advection-diffusion example. The diffusive flux scheme has been developed and validated only for scalar quantities. This can be a problem for boundary conditions which are a function of more than one unknown (e.g. no-slip wall). One way the scheme can be expanded for coupled boundary conditions is to include all of the unknown coefficients of both Taylor polynomial approximations of the variables in the unknown vector of the least-squares problem for the reconstruction. This essentially doubles the size of the constrained least-squares matrix, as well as the total number of constraints, as all the unknown coefficients are solved simultaneously. This coupled reconstruction problem will be very sparse, which can be taken advantage of to help the efficiency of the reconstruction. It was shown in Section 5.5 that the scheme does not work well if there are singular points in the solution. The error in the control volumes around the singular point are larger CHAPTER 7. SUMMARY AND CONCLUSIONS 85 than elsewhere in the solution. This problem might be solved using a different coordinate system (e.g. polar coordinates) around the problem area. The boundary constraint scheme is also currently limited to linear combinations of Dirichlet and Neumann boundary conditions, although other boundary conditions, such as radiation or acceleration, can be derived by similar methods. Extending the diffusive finite-volume scheme to three dimensions is, in principle, straightforward. The flux integral is now an integral through the faces of the control volumes, as opposed to the edges. Gauss quadrature can still be performed with Gauss points located on the surface of the face and the Gauss Constraints method can still be used to add constraints to the reconstruction matrix. However implementation must be done very carefully in order to avoid the problems described in Section 7.1. Bibliography [1] R. Abgrall. On Essentially Non-Oscillatory Schemes on Unstructured Meshes - Analysis and Implementation. Journal of Computational Physics, 114(l):45-58, September 1994. [2] W . Kyle Anderson, Russ D . Rausch, and Daryl L . Bonhaus. Implicit/multigrid algorithms for incompressible turbulent flows on unstructured grids. In 12th AIAA Computational Fluid Dynamics Conference, pages 1067-1081. American Institute of Aeronautics and Astronautics, June 1995. A I A A 95-1740-CP. [3] Timothy 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. A G A R D , 1992. AGARD-R-787. [4] Timothy J . Barth. Recent Developments in High Order &-Exact Reconstruction on Unstructured Meshes. A I A A paper 93-0668, January 1993. [5] Timothy J . Barth and Paul O. Frederickson. Higher Order Solution of the Euler Equations on Unstructured Grids Using Quadratic Reconstruction. A I A A paper 90-0013, January 1990. [6] L . J. Durlofsky, B . Enquist, and S. Oscher. Triangle Based Adaptive Stencils for the Solution of Hyperbolic Conservation Laws. Journal of Computational Physics, 98:64-73, January 1992. [7] O. Friedrich. Weighted Essentially Non-Oscillatory Schemes for the Interpolation of Mean Values on Unstructured Grids. Journal of Computational Physics, 144(1): 194212, July 1998. [8] R. Lohner H . Luo, J. D . Baum and J . Cabello. Adaptive Edge-Based Finite-Element Schemes for the Euler and Navier-Stokes Equations on Unstructured Grids. A I A A paper 93-0336, January 1993. [9] I. Harari and T. J . R. Hughes. Stabilized Finite-Element Methods for Steady AdvectionDiffusion with Production. Computer Methods in Applied Mechanics and Engineering, 115(1-2):165-191, May 1994. [10] C . Q. Hu and C . W . Shu. Weighted Essentially Non-Oscillatory Schemes on Triangular Meshes. Journal of Computational Physics, 150(1):97-127, March 1999. 86 BIBLIOGRAPHY 87 T. J . R. Hughes, L . R Franca, and G . M . Hulbert. A New Finite Element Formulation for Computational Fluid-Dynamics: The Galerkin Least-Squares Method for Advection-Diffusive Equations. Computer Methods in Applied Mechanics and Engineering, 73(2):173-189, May 1989. Antony Jameson. Artificial Diffusion, Upwind Biasing, Limiters and their Effect on Accuracy and Multigrid Convergence in Transonic and Hypersonic Flows. A I A A paper 93-3359, 1993. Dimitri J . Mavriplis. Accurate Multigrid Solution of the Euler Equations on Unstructured and Adaptive Meshes. AIAA Journal, 28(2):213-221, February 1990. Dimitri J . Mavriplis. Turbulent Flow Calculations Using Unstructured and Adaptive Meshes. International Journal for Numerical Methods in Fluids, 13(9), November 1991. Dimitri J. Mavriplis. A Three Dimensional Multigrid Reynolds-Averaged Navier-Stokes Solver for Unstructured Meshes. I C A S E Report No. 94-29, N A S A Langley Research Center, 1994. Dimitri J . Mavriplis, Antony Jameson, and L . Martinelli. Multigrid Solution of the Navier-Stokes Equations on Triangular Meshes. I C A S E Report No. 89-11, N A S A Langley Research Center, 1989. Carl F . Ollivier-Gooch. High-Order E N O Schemes for Unstructured Meshes Based on Least-Squares Reconstruction. A I A A paper 97-0540, January 1997. Carl 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. D. Pellentier, L . Ignat, and F . Ilinca. Adaptive Finite Element Method for Conjugate Heat Transfer. Numerical heat Transfer Part A - Applications, 32(3): 267-287, August 1997. P. L . Roe. Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes. Journal of Computational Physics, 43:357-372, 1981. A . H . Stroud and D . Secrest. Gaussian Quadrature Formulas. Prentice-Hall, Englewood Cliffs, N . J . , 1966. M . Zamecnik and J . A . Hernandez. Essentially Non-Oscillatory Schemes for the Euler Equations on Unstructured Meshes. Proceedings of the Institution of Mechanical Engineers Part G - Journal of Aerospace Engineering, 213:1-12, 1999. Appendix A These are some of the meshes used in validating the flux integral calculations. They are the meshes coarse enough to be viewable when printed: Figure A . l : squareO.mesh — 145 cells, 89 vertices 88 APPENDIX A. Figure A.3: square2.mesh — 2280 cells, 1204 vertices APPENDIX A. Figure A.5: circlel.mesh — 1617 cells, 854 vertices APPENDIX A. Figure A.6: cardioidO.mesh — 1673 cells, 879 vertices
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- High-order finite-volume discretisations for solving...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
High-order finite-volume discretisations for solving a modified advection-diffusion problem on unstructured… Van Altena, Michael 1999
pdf
Page Metadata
Item Metadata
Title | High-order finite-volume discretisations for solving a modified advection-diffusion problem on unstructured triangular meshes |
Creator |
Van Altena, Michael |
Date Issued | 1999 |
Description | A fourth-order accurate diffusive flux calculation scheme has been developed which can be incorporated into any existing reconstruction-based unstructured mesh inviscid flow solver. The flux scheme makes use of a Laplace problem to model the viscous terms and of the existing least-squares reconstruction code used in Dr. Carl F. Ollivier-Gooch's inviscid flow solver. Boundary conditions on the solution are enforced during the reconstruction, which requires boundary constraints to be included in the least-squares matrix. The Taylor Constraints method for generating boundary constraints is developed and then discarded due to its drawbacks. The Gauss Constraints method is then developed which is simpler and avoids the drawbacks of the Taylor Constraints method. The flux integration is validated to demonstrate fourth-order accuracy on square, circle and cardioid meshes with a combination of Dirichlet and Neumann boundary conditions. The flux integration scheme is shown to be fourth-order accurate for a modified advection-diffusion problem of the same form as the Navier-Stokes' momentum equations. The diffusive flux calculation scheme can later be combined with an existing inviscid flow solver to solve the viscous compressible Navier-Stokes' equations. |
Extent | 10503020 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-06-29 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080945 |
URI | http://hdl.handle.net/2429/9800 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1999-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1999-0633.pdf [ 10.02MB ]
- Metadata
- JSON: 831-1.0080945.json
- JSON-LD: 831-1.0080945-ld.json
- RDF/XML (Pretty): 831-1.0080945-rdf.xml
- RDF/JSON: 831-1.0080945-rdf.json
- Turtle: 831-1.0080945-turtle.txt
- N-Triples: 831-1.0080945-rdf-ntriples.txt
- Original Record: 831-1.0080945-source.json
- Full Text
- 831-1.0080945-fulltext.txt
- Citation
- 831-1.0080945.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080945/manifest