D E V E L O P M E N T A N D VERIFICATION O F T H E G E N E R I C T H R E E DIMENSIONAL FINITE V O L U M E S O L V E R By M. K. Harsha Perera MPhil., The University of Southampton, UK., 2002 M.Eng., The University of Southampton, UK., 1999 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (The Department of Mechanical Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 2003 © M. K. Harsha Perera, 2003 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of mgcu^Mi mi. B^xnxtMB£Ri N (h The University of British Columbia Vancouver, Canada Date O-o //( /g>3 DE-6 (2/88) 11 A b s t r a c t The generic finite volume solver, ANSLib, has been extended to three dimensions and used to verify the accurate computation of three-dimensional advection-diffusion and Poisson problems. A simple cubic domain has been selected as the domain of interest. Over this domain a steady state solution is computed for each model problem with 2nd, 3rd and 4th order accurate schemes. Gauss quadrature is used to evaluate the flux integral to 2nd, 3rd and 4th order accuracy. The flux scheme makes use of the centred scheme to model the viscous term and a simple upwind scheme to model the advective fluxes. In both cases the existing k-exact least-square reconstruction code is used to obtain the solution at each Gauss integration point to 2nd, 3rd and 4th order accuracy. Newton-Krylov implicit time advance is used to obtain the steady state solution, implemented using the GMRES algorithm. Prior to subjecting ANSLib code to physical model problems, the code is first tested for the correct enforcement of boundary constraints, accurate implementation of the solution reconstruction, and flux integral evaluation in three dimensions. Relevant changes are made to the code to obtain the desired behaviour. Analytical solutions for the advection-diffusion and Poisson equations are derived to satisfy Dirichlet and Neumann boundary conditions; these are used to compute the resulting discretized error in numerically converged solutions. Accuracy assessment is completed through a comparison of error norms for different mesh sizes for each test problem. Localized errors on the domain for each problem at steady state are plotted. The efficiency of the code to obtain a numerical converged solution is reported. iii T a b l e o f C o n t e n t s A B S T R A C T I I L I S T O F F I G U R E S V L I S T O F T A B L E S V I I A C K N O W L E D G E M E N T V I I I 1 I N T R O D U C T I O N 1 1.1 U N S T R U C T U R E D H I G H O R D E R M E T H O D S 1 1.2 OBJECTIVE A N D SCOPE 2 1.3 THESIS STRUCTURE 3 2 G E N E R I C S O L V E R 4 2.1 G E N E R I C SOLVER D E V E L O P M E N T 4 2.2 OPERATION OF THE A N S L I B SOLVER 5 2.3 C O N C L U S I O N 5 3 T H R E E - D I M E N S I O N A L F L O W S O L V E R 8 3.1 SPATIAL DISCRETIZATION 8 3.1.1 Higher order methods for the convective terms 9 3.1.2 Higher order methods for the viscous terms 11 3.2 SOLUTION UPDATE TECHNIQUES 12 3.2.1 Time stepping 13 3.3 CONCLUSION 13 4 H I G H E R O R D E R R E C O N S T R U C T I O N A N D V E R I F I C A T I O N 14 4.1 K - E X A C T RECONSTRUCTION 14 4.2 VERIFICATION TEST PROGRAMS 16 4.2.1 Boundary Condition treatment 17 4.3 C O N C L U S I O N 19 5 M O D E L V E R I F I C A T I O N 2 0 5.1 F L U X INTEGRAL TEST 2 0 5.1.1 Advection-diffiision 21 5.1.2 Poisson 23 5.2 S T E A D Y STATE SOLUTION 2 6 5.2.1 Advection-diffiision results 2 6 5.2.2 Poisson results 3 4 5.3 C O N C L U S I O N 41 6 C O N C L U S I O N S A N D F U T U R E W O R K 4 2 7 B I B L I O G R A P H Y 43 A A P P E N D I X A 4 5 A. 1 G A U S S QUADRATURE FOR A TRIANGLE 45 A .2 G A U S S QUADRATURE FOR A TETRAHEDRON 45 A.3 P A R A L L E L AXIS THEOREM 4 6 A .4 E X A C T INTEGRATION F O R M U L A OVER TETRAHEDRON 4 6 iv A . 5 C A L C U L A T I N G ORDER OF A C C U R A C Y F R O M ERROR NORMS 47 A.6 A N A L Y T I C A L SOLUTION METHODS 48 A.6.1 Advection-diffusion 48 A.6.2 Poisson 49 V L i s t o f F i g u r e s FIGURE 2-1 F L O W CHART SHOWING THE INTERACTION OF VARIOUS OBJECTS IN THE A N S L I B CODE. T H E BOXES WITH B O L D B O U N D A R Y A N D G R E Y BOXES SHOW SOME OF T H E CLASSES. T H E DASHED LINES SHOW THE INTERACTION OF THE D R I V E R P R O G R A M WITH T H E CLASSES, A N D T H E B L A C K LINES SHOW SOME INFORMATION PASSING BETWEEN CLASSES. T H E B O X E S WITH B O L D BOUNDARIES A R E THE ONES AFFECTED B Y THE THIRD DIMENSION 7 FIGURE 3-1 U N S T R U C T U R E D MESH OVER A CUBIC D O M A I N 8 FIGURE 4-1 O R D E R OF A C C U R A C Y OF THE RECONSTRUCTION OF FUNCTION IN E Q . (4 .4A) FOR T H E INTERIOR CONTROL V O L U M E S 17 FIGURE 4-2 T H E B O U N D A R Y FACE DIRECTIONS A N D L A B E L S FOR THE CUBIC D O M A I N 18 FIGURE 5-1 T H E L 2 ERROR N O R M FOR THE ADVECTION-DIFFUSION F L U X TEST FOR VARIOUS M E S H SIZES 23 FIGURE 5-2 T H E L 2 ERROR N O R M FOR THE POISSON F L U X TEST FOR VARIOUS M E S H SIZES 2 4 FIGURE 5-3 T H E L 2 ERROR N O R M FOR THE ADVECTION-DIFFUSION P R O B L E M A T STEADY STATE FOR VARIOUS M E S H SIZES 2 6 FIGURE 5-4 T H E SUB FIGURES (A)-(D) SHOWS THE SOLUTIONS TO T H E ADVECTION-DIFFUSION P R O B L E M A L O N G THE GEOMETRY A T X=0, 0 .3 ,0 .6 A N D 1 2 7 FIGURE 5-5 E R R O R BETWEEN THE E X A C T A N D 2 N D ORDER A C C U R A T E SOLUTION, AFTER N U M E R I C A L C O N V E R G E N C E OF T H E ADVECTION-DIFFUSION P R O B L E M 2 9 FIGURE 5-6 E R R O R BETWEEN THE E X A C T A N D 3 R D ORDER A C C U R A T E SOLUTION, A F T E R N U M E R I C A L C O N V E R G E N C E OF THE ADVECTION-DIFFUSION P R O B L E M 3 0 FIGURE 5-7 E R R O R BETWEEN THE E X A C T A N D 4 T H ORDER A C C U R A T E SOLUTION, AFTER N U M E R I C A L C O N V E R G E N C E OF THE ADVECTION-DIFFUSION P R O B L E M 3 0 FIGURE 5-8 T H E COMPUTATIONAL TIME T A K E N FOR THE STEADY-STATE SOLUTION OF T H E ADVECTION-DIFFUSION P R O B L E M FOR VARIOUS MESHES 31 FIGURE 5-9 T H E L 2 RESIDUAL N O R M FOR THE ADVECTION-DIFFUSION P R O B L E M FOR T H E 2 N D ORDER A C C U R A T E METHOD 3 2 FIGURE 5-10 T H E L 2 RESIDUAL N O RM FOR THE ADVECTION-DIFFUSION P R O B L E M FOR THE 3 R D ORDER A C C U R A T E METHOD 33 FIGURE 5-11 T H E L 2 RESIDUAL N O R M FOR THE ADVECTION-DIFFUSION P R O B L E M FOR T H E 4 T H ORDER A C C U R A T E METHOD 33 FIGURE 5-12 T H E L 2 RESIDUAL N O R M FOR THE ADVECTION-DIFFUSION P R O B L E M FOR T H E 2 N D ORDER A C C U R A T E METHOD WITH K R Y L O V SUBSPACE VECTORS SIZES OF 3 0 A N D 60 3 4 FIGURE 5-13 T H E L 2 ERROR N O R M FOR THE POISSON P R O B L E M A T STEADY STATE FOR VARIOUS M E S H SIZES 3 5 FIGURE 5-14 T H E SUB FIGURES (A)-(C) SHOW THE SOLUTIONS TO THE POISSON P R O B L E M A L O N G T H E GEOMETRY A T X=0, 0.3, A N D 0.6 35 FIGURE 5-15 E R R O R BETWEEN THE E X A C T A N D 2 N D ORDER A C C U R A T E SOLUTION, A F T E R N U M E R I C A L C O N V E R G E N C E FOR THE POISSON P R O B L E M 3 7 V I FIGURE 5-16 E R R O R BETWEEN THE E X A C T A N D 3 R D ORDER A C C U R A T E SOLUTION, A F T E R N U M E R I C A L C O N V E R G E N C E FOR THE POISSON P R O B L E M 38 FIGURE 5-17 E R R O R BETWEEN THE E X A C T A N D 4 T H ORDER A C C U R A T E SOLUTION, AFTER N U M E R I C A L C O N V E R G E N C E FOR THE POISSON P R O B L E M 38 FIGURE 5-18 T H E COMPUTATIONAL TIME T A K E N FOR THE STEADY-STATE SOLUTION OF T H E POISSON P R O B L E M FOR VARIOUS MESHES. . . . . 3 9 FIGURE 5-19 T H E L 2 RESIDUAL NORM FOR THE POISSON P R O B L E M FOR THE 2 N D ORDER A C C U R A T E METHOD 4 0 FIGURE 5-20 T H E L 2 RESIDUAL N O R M FOR THE POISSON P R O B L E M FOR THE 3 R D ORDER A C C U R A T E METHOD 4 0 FIGURE 5-21 T H E L 2 RESIDUAL N O R M FOR THE POISSON P R O B L E M FOR THE 4 T H ORDER A C C U R A T E METHOD 41 FIGURE A - 7 - 1 T E T R A H E D R O N CONTROL V O L U M E 4 6 V l l L i s t o f T a b l e s T A B L E 4-1 N U M B E R OF CELLS USED IN T H E THREE DIMENSIONAL LEAST-SQUARE RECONSTRUCTION STENCIL 15 T A B L E 4-2 A L L DIRICHLET CONSTRAINED B O U N D A R Y CONDITION TEST 19 T A B L E 4-3 D I R I C H L E T - N E U M A N N CONSTRAINED B O U N D A R Y CONDITION TEST 19 T A B L E 5-1 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR A D V E C T I O N -DIFFUSION F L U X TEXT FOR 2 N D ORDER A C C U R A T E METHOD 2 2 T A B L E 5-2 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR A D V E C T I O N -DIFFUSION F L U X TEXT FOR 3 R D ORDER A C C U R A T E METHOD 2 2 T A B L E 5-3 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR A D V E C T I O N -DIFFUSION F L U X TEXT FOR 4™ ORDER A C C U R A T E METHOD 2 2 T A B L E 5-4 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR POISSON F L U X TEST FOR THE 2 N D ORDER A C C U R A T E METHOD 25 T A B L E 5-5 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR POISSON F L U X TEST FOR THE 3 R D ORDER A C C U R A T E METHOD 25 T A B L E 5-6 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR POISSON F L U X TEST FOR THE 4™ ORDER A C C U R A T E METHOD 25 T A B L E 5-7 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR A D V E C T I O N -DIFFUSION P R O B L E M FOR THE 2 N D ORDER A C C U R A T E METHOD 2 8 T A B L E 5-8 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR A D V E C T I O N -DIFFUSION P R O B L E M FOR THE 3 R D ORDER A C C U R A T E METHOD 2 8 T A B L E 5-9 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR A D V E C T I O N -DIFFUSION P R O B L E M FOR THE 4™ ORDER A C C U R A T E METHOD 2 8 T A B L E 5-10 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR POISSON P R O B L E M FOR THE 2 N D ORDER A C C U R A T E METHOD 3 6 T A B L E 5-11 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR POISSON P R O B L E M FOR T H E 3 R D ORDER A C C U R A T E METHOD 3 6 T A B L E 5-12 T H E C A L C U L A T E D ERROR NORMS A N D ORDER OF A C C U R A C Y FOR POISSON P R O B L E M FOR T H E 4™ ORDER A C C U R A T E METHOD 3 6 T A B L E A - l G A U S S POINTS A N D WEIGHT FOR A TRIANGLE 45 T A B L E A - 2 G A U S S POINTS A N D WEIGHT FOR A TETRAHEDRON 45 T A B L E A - 3 S A M P L E T A B L E OF L 2 ERROR NORMS FOR E A C H M E S H 4 7 V l l l A c k n o w l e d g e m e n t I would like to say thanks to my supervisors, Dr. Carl Ollivier-Gooch and Dr. Kendal Bushe, for their great guidance, support and patience throughout my years of study in U B C . I am grateful for my student colleagues in the research group for their assistance on countless occasions. 1 1 I n t r o d u c t i o n Computational Fluid Dynamics (CFD) is used mainly as a design tool. It computes the flow patterns for a given set of initial and boundary conditions. The computed flow field provides a good measure of the governing fluid and thermal forces. The combination of a versatile CFD program and an experienced user is necessary for the generation and interpretation of computational results, for efficient application of CFD. The CFD code must be verified to prove that the algorithm performs to expectations, and validated to prove that computed solutions are accurate, with respect to experimental data. Before each new application, the user must carry out a verification and validation study. However, a code with prior verification on similar problems can provide insight into the code's ability to solve problems. Code verification forms an important part of CFD application and development. In nearly every thermo-fluid problem, three-dimensionality is relevant. A generic finite volume solver, called ANSLib (Advanced Numerical Simulation Library), has been extended to three dimensions to solve physically realistic problems. ANSLib has generic high order solution methods, which are currently working well in two dimensions, in both structured and unstructured discretized domains. High order accurate methods provide an efficient means for obtaining an accurate solution on a coarse grid while maintaining monotonicity. In this thesis, the process of ANSLib program development, to include three dimensional functionality is revealed through a verification study. Sample problems, including the advection-diffusion problem and the Poisson problem are solved using the three-dimensional ANSLib solver, and some characteristics of the generic algorithms are shown. 2.1 Unstructured High Order Methods Computational methods are constantly required to deliver more accurate solutions over complex configurations. Unstructured mesh methods have emerged as a means of meeting this requirement. Discretizing the complex geometry with an unstructured mesh method has replaced the traditional block-structured mesh techniques, which require more user interaction to produce a good quality mesh. A n unstructured mesh enables straightforward adaptive meshing techniques to be applied to its elements, which are often triangular in two-dimensions and tetrahedral in three-dimensions. 2 Most commercial software packages (e.g. Fluent 5™, CFX-4.2™) make the costly step of using a finer mesh necessary to provide an accurate solution, since in these programs, the order of accuracy is second order at best. Considering that a smooth and a fast second order scheme that require a larger mesh, a more accurate solution method can be expected to achieve the same level of accuracy on a coarser mesh, and, possibly, at the same time be computationally inexpensive. Higher-order solution reconstruction methods are successfully applied to compressible flow problems, and provide a low cost option over low order, fine mesh computations. Originally, these schemes had been applied to accurately resolve discontinuities1. However, the principles upon which these methods are built are equally applied to incompressible flows, where it is still necessary to resolve steep solution gradients in shear layers with a minimal amount of numerical smearing. The advective fluxes are computed with an upwind scheme on unstructured meshes. In general, high order upwind schemes are oscillatory and require limiters to damp them down. They need approximate Riemann solvers to accurately evaluate fluxes at mesh cell faces. Such methods are computationally very expensive. The k-exact reconstruction method, used in a modified way in ANSLib has the properties of non-oscillatory and has good convergence, thus producing a high order solution. For diffusion fluxes, the same k-exact reconstruction scheme can be used to obtain a higher order accurate solution. Gauss quadrature is used to evaluate the flux integral over the cell faces, as well as to evaluate the source term. In ANSLib, an explicit multistage scheme and an implicit Newton-Krylov scheme are available to solve the discretized equations. Implicit time advance will be used in three dimensions since the Newton-Krylov scheme, which uses the Generalized Minimal Residual (GMRES) 2 , has the added advantage of being able to accelerate the rate of convergence. 2.2 Objective and scope The ANSLib solver, verified for two-dimensional problems5, is developed and verified in three dimensions. The verification process for the three-dimensional case is conducted in two stages. The first stage is verification of the performance of the generic algorithms in three dimensions. The algorithms are tested for high-order solution reconstruction, mesh pre-processing, Gauss integration, and flux integration. The second stage is to verify that the model equations, solved to steady state, show the desired order of accuracy. The global and local numerical errors are 3 reported, as well as the residual convergence histories and computational efficiency for each order of accuracy. 2.3 Thesis structure The operation of the ANSLib solver and the requirements for three-dimensional functionality is explained in Chapter Two. Some generic solver attributes are discussed. In Chapter Three, the discretization and iterative schemes, found in the ANSLib solver, are compared with existing numerical schemes, found in the literature. The discretization techniques discussed are high-order methods with an emphasis on unstructured grid techniques. The upwinding used in most unstructured solvers requires either explicit or implicit solution methods; implicit methods, incorporating solution acceleration, have proven to be the most efficient. The existing k-exact type reconstruction technique in the ANSLib solver is explained in Chapter Four. After making the necessary adjustments to the three-dimensional mesh class, the implementation of the reconstruction scheme is checked in three dimensions. Tests are conducted to show that the given reconstruction algorithm function in three dimensions satisfactorily, as well as how the reconstruction scheme reconstructs the boundary solutions. In Chapter Five, flux integration and steady state verification tests are conducted on the advection-diffusion and Poisson model problems, where the exact analytical solutions are used to compare the accuracy of the computed solution. The model verification results include performance of the 2nd, 3rd and 4 t h order schemes with respect to overall accuracy, convergence, and computational time. In Chapter Six, the thesis is concluded with a discussion of the overall status of the code and future plans. 4 2 G e n e r i c S o l v e r Generic solvers are able to perform numerical computations on many different physical problems governed by different Partial Differential Equations, given that they can be expressed in some generic form using the same numerical procedure. Most generic solvers are built using C/C++ Object-Oriented Programming (OOP) techniques, which allows programmers to define shared procedures on common mesh classes (structured and unstructured meshes), and include linear system solvers, to function independently of other classes. Programmers can also control volume integration, surface integration, flux accumulation, and error norm evaluation for all mesh elements. There is a common post-processing class as well. ANSLib contains this integration of object classes and more, where common tasks and common data management are done to make it an effective generic solver. 2.1 Generic solver development The design of a generic solver is an evolving art, undertaken by many researchers. Programs using the finite element (FE) method were employed in the past as generic solvers3. These programs required the user to have in-depth knowledge of the internal finite-element algorithms when tackling new physical problems or when introducing new discretization techniques. Jinn-Liang Liu et al. 4 used OOP to implement both the FE and finite volume (FV) methods for solving various mathematical problems using adaptive implementation, which had not been previously addressed. ANSLib, on the other hand, allows the maximum flexibility in the numerical schemes with the least dependence on the physics5. It provides a robust and general numerical support while being independent of the underlying mesh topology. The ANSLib code provides an environment that can be used to represent many different mathematical models found in CFD and Computational Solid Mechanics (CSM). The C S M models that were originally solved using FE techniques can now be solved using the same F V algorithms found in CFD. Developments undertaken in the two dimensional ANSLib infrastructure include automatic mesh generation on complex geometries and coupling C S M with CFD models on a single global domain6. Generic solvers need robust solver algorithms, since they must be able to solve many differently conditioned systems of equations on irregular meshes. Hence, the solver algorithms require a thorough verification and validation before being applied to realistic problems. 5 2.2 Operation of the ANSLib solver The ANSLib generic solver software is operated by the Driver program, which reads user input and selects the instructions to execute. For example, the user provides the following command information: name of a model program (e.g. Poisson), three-dimensional, unstructured, steady state, 3 r d order, Krylov, 30 Krylov subspaces, 50 iterations, residual limit (e.g. l.Oe'1 3), and restart. With this information, the Driver program initializes the necessary classes and program variables, then coordinates the program's execution. Figure 2-1 shows a flow diagram of this process. The dashed lines represent the Driver program's execution, from left to right. The physics object is created by the Driver program, which is linked to the problem definition functions in a model program, e.g. Poisson. The model program's various functions define the flux discretization of various terms in the model problem, for interior and boundary control volumes, provide a platform specifying the boundary conditions, initialise solutions, and define the source term. The program also has optional functions for the definition of exact solution and exact solution divergence, which is used to test the flux integral and to determine the discretization error of a numerically converged solution. When the Driver program executes the chosen iterative method (Multistage, e.g. Runge-Kutta or Newton-Krylov), the iterative method needs the flux class. The purpose of the flux class is mainly to accumulate fluxes, using the physics class. The physics class is linked to functions in model programs such as Poisson. In this way, the information regarding the discretized model gets into the iterative solver. The black arrows in Figure 2-1 show only a few interactions between classes. The FieldQuant class deals with issues related to the storage of the solution. The Recon class deals with reconstruction of the solution, and the Integrate program has a set of functions, which can supply information for solution integration over control volumes, and computes error norms. When the iteration count or residual limit is reached, control returns to the Driver program. The final step is to execute the appropriate action in the output class, where restart files are written and the solution is formatted into files appropriate for visualisation and analysis. 2.3 Conclusion It is clear from Figure 2-1 that the ANSLib solver uses object oriented programming techniques, and that solver actions are governed by what each class decides to do. This important object-oriented programming property allows ANSLib to expand, and facilitate independent 6 maintenance and testing of classes with respect to certain other classes. The original code has no prior implementation of any three dimensional flow models, nor has any code segment been verified in three dimensions. Adding three dimensionality to the original generic two dimensional ANSLib code implementation requires modification, assessment and testing of only the boxes with bold boundary. Driver Physics Solution Update i i Problem Setup \iSititj≥- • Krylov L i Interior Flux j-•j Boundary Flux \-*\ Source | •f Boundary Type i H Initial Solution ± Held Quant Exact Solution Exact Divergence • t Output Mesh Recon Integrate I —•) Number of unknowns 2-1 Flow chart showing the interaction of various objects in the ANSLib code. The boxes with bold boundary and grey boxes show some of the classes. The dashed lines show the interaction of the Driver program with the classes, and the black solid lines show some information passing between classes. The boxes with bold boundary are the ones affected by the third dimension. -.^i'V . "• " '•'V. ;:-v-"'~:-8 3 T h r e e - D i m e n s i o n a l F l o w S o l v e r An unstructured tetrahedral mesh is used to discretize the domain. In each tetrahedral control volume, the model partial differential equation is integrated and transformed into a surface integral form. The flux integral is solved using Gauss-quadrature and the fluxes are evaluated at each gauss point that lies on the triangular face. Existing k-exact reconstruction is used in ANSLib to reconstruct the solution at each Gauss point to high order accuracy, as discussed in section 3.1, in relation to other high order schemes used to evaluate the fluxes. The reconstructed high order solution is used to evaluate the fluxes to high order accuracy at each Gauss integration point. In section 3.2, solution advance schemes are explained, with emphasis on implicit time advance. 3.1 Spatial Discretization The tetrahedral unstructured mesh, over the cubic geometry, was created using the latest version of G R U M M P 7 software. The cell-centre is used as the reference point to compute the control volume averaged variables. A sample mesh is shown in Figure 3-1, containing 2112 tetrahedral control volumes. Figure 3-1 Unstructured mesh over a cubic domain For a given mesh, the general three dimensional model PDE can be written in the following form: 9 dU(x,y,z) | dF(x,y,z) t 9G(x,y,z) | dH(x,y,z) S(x,y,z) (3-1) where U is the solution vector, S is the source term, and F, G, and H are the flux terms. The finite volume form can be written as: where Vj is the cell size for the i control volume, A is the cell face area, and n is the outward pointing normal. The flux integral is approximated by the Gauss-quadrature: where wj is the Gauss weight, and U(XJ) is the solution evaluated at the j Gauss location. (Table of Gauss locations and weights for a triangle and for a tetrahedron is given in Table A - l and Table A-2 respectively, in Appendix A.), the flux vector is evaluated with a simple upwind scheme: where U + ( X J ) is the U reconstructed on the outside of C V , and U"(XJ) is U reconstructed on the inside of C V . The flux direction is chosen according to the direction of the solution velocity vector. The convective and viscous flux terms on the right hand side of Eq. (3.4) are evaluated at a given Gauss point by a higher order method. 3.1.1 Higher order methods for the convective terms A general low order upwind method usually has excessive numerical dissipation and exhibits strong smearing of the solution, leading to poor accuracy and resolution. Higher order methods, in upwind schemes, are constructed by reconstructing approximations to higher-order derivatives of the solution variables on either side of the control volume face. They can be used to extrapolate more accurate values to either side of the control volume face. The resulting interface values are discontinuous, but the jump in discontinuities decreases with the increasing order of o accuracy of the scheme. Third order upwind schemes, such as the Q U I C K scheme, applied in structured discretization, improves the accuracy of the computation by increasing the order of the +— d F i + G j + Hk n d A - S = 0 (3.2) (3.3) (3-4) 10 leading order truncation error, however, the solution displays unphysical oscillations in regions of steep gradients, causing numerical instability. Thus, all high-order schemes need to reduce or prevent spurious oscillations. Applying a limiter, based on the solution gradient, improves the monotonicity of the upwind schemes, resulting in a much more robust scheme. The Flux Corrected Transport (FCT) 9 and Total Variation Diminishing (TVD) 1 7 methods are other high order methods that began with structured meshes. In both schemes, a limited amount of anti-diffusive flux is added to attain monotonicity. In addition, these schemes employ non-linear limiters on all regions of the flow field, which may adversely affect the convergence. Van Leer 1 0 devised the Monotonic Upstream-Centred Scheme for Conservation Laws (MUSCL) to achieve a higher order accurate solution. M U S C L schemes rely on a piecewise-polynomial reconstruction procedure using nonlinear limiters to enforce the monotonicity. A Riemann solver is often employed to calculate the flux at the interface for given right and left interface states on solution domains with discontinuities. However, an exact solution to a Riemann problem is too costly, therefore the approximate Riemann solver from Osher11 and Roe 1 2 is used for this purpose. M U S C L schemes were applied initially to structured meshes, and then later to unstructured meshes. Another scheme to improve the flux evaluation is the fluctuation splitting scheme13, which provides an alternative to Riemann solvers incorporating upwinding, however they have second order accuracy only at steady state. Most of the schemes discussed, do not guarantee in the absence of oscillations in the multi-dimensional case14. However, Barth and Jespersen15 created a multidimensional limiter for unstructured meshes that satisfied the monotonicity principle in multidimensions. In this scheme, the reconstructed distribution in the control volume is bounded everywhere by the values of the neighbours. However, their scheme suffered from convergence difficulties and a method to solve the problem came at the expense of monotonicity16. There are many high-order, non-linear, stable "shock-capturing" schemes developed for the compressible flows 1 7. The philosophy behind these schemes is to accurately resolve the rapid transition regions or shocks in a stable and globally correct fashion, yet at the same time have a high resolution over the smooth part of the flow. Thus, when fully resolving the flow is impossible or too costly (e.g. combustion flows), ENO (Essentially Non-Oscillatory) or WENO (Weighted ENO) schemes on a coarse grid can be used to obtain at least partial information about the flow18. 11 E and Shu successfully applied high order ENO schemes to two dimensional incompressible flows in order to resolve shear roll-ups. Both WENO and ENO schemes use the concept of adaptive stencils in the reconstruction procedure, based on the local smoothness of the numerical solution, to automatically achieve high order accuracy and nonoscillatory properties near discontinuities. ENO uses one out of many candidate stencils whereas WENO uses the convex combination of all candidate stencils, each assigned with a nonlinear weight depending on the local smoothness of the numerical solution based on that stencil. Another higher order scheme, initially applied to capture discontinuities, is Barth and Fredrickson's1 9 k-exact reconstruction. It involves the reconstruction of variables, satisfying the properties of conservation of the mean, k-exactness (able to reconstruct exact polynomials of degree < k), and compact support. Combining the properties of k-exact reconstruction and data dependent ENO schemes, Ollivier-Gooch 2 0, has constructed the DD-L2-ENO scheme. It employs a least square technique to estimate the higher order derivatives and provide a better approximation on distorted meshes. The k-exact reconstruction scheme has the following properties: • conservation of the mean value of the function in each C V • uniformly accurate with no overshoots • good steady state convergence • computationally efficient • easily implemented on arbitrary meshes The higher-order k-exact reconstruction, used in the three-dimensional ANSLib solver, is slightly different from the D D - L 2 - E N 0 2 0 scheme in the sense that the data dependent weight function is not calculated, and the latter is especially needed to evaluate discontinous solutions. Details on this reconstruction scheme and its verification are given in chapter 4. 3.1.2 Higher order methods for the viscous terms Viscous terms are diffusive in nature, and are usually centrally differenced. Higher order methods for diffusive fluxes have not received much attention. Ollivier-Gooch and Van Altena 2 1 showed that k-exact type reconstruction could be used to compute viscous fluxes to higher than second order accuracy. This scheme uses the same k-exact type reconstruction discussed for advective fluxes and it has been implemented in ANSLib. 12 3.2 Solution update techniques The Eq. (3.2) is discretized to obtain the system of coupled ordinary differential equations ^ - + — R(U) = 0 (3.5) dt V where R(U) is the spatial discretized terms, and V is the cell volume. The Eq. (3.5) can be solved by either explicit or implicit methods. Explicit methods include the family of linear multi-step methods, such as the Runge-Kutta schemes. They are stable for a range of Courant numbers where the maximum time step is proportional to local cell size. Thus, because of the time step restriction, convergence to a steady state is slow. A n implicit scheme can be obtained by evaluating the spatially discretized term at the new time level, n+1. After linearizing about the current time level, Eq. (3.6) is obtained. —+ ^ | A U = -R ( U n ) (3.6) At a u j 3R where -^ j is the Jacobian matrix, A U is the change in solution vector which the system must solve for, and At is the time step. The updated solution is U n + 1 = U n + A U . If the Jacobian is exactly linearized, a large sparse matrix must be inverted, which is computationally expensive. Therefore, the Jacobian is often replaced with a first order discretization. The system of equations can then be solved, using iterative methods, but the convergence is slow. In order to increase the convergence rate pre-conditioned GMRES methods are used2 2. The Incomplete Lower Upper (ILU) factorization could be used as the pre-conditioner, to create a more efficient rate of convergence. However, such a scheme still consumes considerable memory, especially when many dependent variables are present. A solution is to use the matrix free method, which evaluates the matrix vector product by finite-difference techniques: 3R A U = = R ( U + £ A U ) - R ( U ) dU e 13 where e is chosen as the square root of machine zero. An algorithm with pre-conditioned, matrix-free GMPvES would be most efficient. The three-dimensional ANSLib code uses this matrix free implicit Newton-Krylov solver, but without pre-conditioning. 3.2.1 Time stepping The time step given in Eq. (3.6) is calculated by Eq. (3.8). A t ^ - p ^ (3.8) dC • n • dA j max.i i 8A (X, j - x B ,)-n = y _ R - J ' xC (3.9) | X L , i X R , j | where C m a x is the fastest wave entering the control volume, determined from the flux function. For advective fluxes it is the advective velocity, and for the diffusive fluxes the C m a x term will contain a pseudo wavespeed component, C p s eudo,i for the control volume, i. An estimate of an equivalent pseudo-wavespeed is given by Eq. (3.9). In Eq. (3.9), the i and j are control volume indices on either side of the cell face, and the constant C is related to the diffusion coefficient. The denominator of Eq. (3.8), which shows an integration of wavespeeds around control volume, i, is evaluated using Gauss-quadrature. For steady-state solutions, local value of the time step for each control volume is used. 3.3 Conclusion The spatial discretization is performed on an unstructured mesh, using an upwind scheme over a cubical domain. The solutions at quadrature points on the cell faces are evaluated to high order of accuracy, using k-exact reconstruction. The advective flux at each Gauss point is evaluated using simple upwind method, where the left or right solution states are reconstructed to high order accuracy. The diffusive fluxes are evaluated using a central method. The flux integral and source term are computed using Gauss-quadrature. The reconstruction order of accuracy matches the order of accuracy of the Gauss integration. For compressible flows, ENO and WENO schemes provided the most robust and non-oscillatory accurate solution. The principles used in these schemes can be applied for incompressible problems as well. The resulting discretized equations are solved to the steady state by implicit time advancement, using a matrix-free Newton-Krylov solver. 14 4 H i g h e r O r d e r R e c o n s t r u c t i o n a n d V e r i f i c a t i o n The k-exact reconstruction scheme, discussed in section 3.1, has been implemented in the ANSLib code. However, the implementation in three dimensions has not been verified. In section 4.1, the k-exact reconstruction is briefly described in stages, and techniques to verify each stage, are given with results. For more detailed information on k-exact reconstruction and DD-L2-ENO, the reader is referred to the publications of Bath and Frederickson19 and Ollivier-Gooch20 respectively. 4.1 k-exact reconstruction The reconstruction polynomial is a solution expansion, <bjR(x-Xi), and is given by Eq. (4.1). • f <*-*.)=4 + | | + + a 2u dx2 a 2u +-3y2 a 2u dz2 , \ du, (x-x i) + — I dy ( y - y , ) + s l (x - X i ) 2 a 2u 2 dxdy (y-yi) 2 , a 2u dy3z (z- Z i ) 2 9 2u — + dzdx (z- Z j ) (x-XiXy-Yi) i (y-yiKz-zO (z-z j)(x-x i) + H.O.T (4.1) where is the value of the reconstruction solution. The reconstruction solution and its derivatives are evaluated at the reference point (xj,yi,Zj). The coefficients of the expansion are chosen to conserve the mean value in the control volume and to minimise error in representing smooth solutions, as shown in Eq. (4.2). 15 i i J_f ( x_X i). d V + ^ X\(y-yi).dV + ^ - - f ( Z - Z i ) d V + + a 2 u dx2 a 2 u +• 3y 2 a 2 u , 2 V M dxdy a 2 u dz2 ^ j ( y - y i ) 2 - d v + a y 3 z M c x - x ^ y - y ^ - d V i i i M c y - y i X z - z ^ - d v (4.2) , v , , • V: J I ' 1 J - J ( z - z i ) ( x - x i ) d V + H.O.T i V i . where V; is the cell volume. The volume averaged monomials, the moments, can be expressed using Green's theorem to convert it to a boundary integral, and they can be integrated exactly, using Gaussian quadrature with points and weights over triangles as shown in Table A - l . The expansion must be carried out through the k t h derivative, then the reconstruction can be said to be k-exact. The mean value of solution given by the left hand term of Eq. (4.2) and the moments are known. The least-squares method attempts to minimise the error in predicting the mean value of the reconstruction polynomial for control volumes, within a fixed stencil. In this way the derivatives at Xj are chosen20. A l l geometrical terms are pre-processed, and computed only once. The reconstruction stencil is compact, i.e. only cells that lie within a small neighbourhood near the cell centre are used. The required number of cells for the least-square reconstruction for each order in three-dimensions are shown in Table 4-1, and they are the upper limit used to solve the least-square reconstruction matrix. This number exceeds the number of derivatives that must be solved for in order to obtain the k-exact reconstruction. Table 4-1 Number of cells used in the three dimensional least-square reconstruction stencil Order Required number of cells in stencil 2 4 3 15 4 30 The least-square problem is solved using Householder transformation in order to reduce the coefficient matrix to upper triangular form. With this method, solutions that are more accurate can be obtained, especially for ill-conditioned matrices, because normal equation methods will have conditioning problems. The back substitution is used to yield the required derivatives. 16 These derivatives are then used to compute solutions at the Gauss points, up to k order accuracy. 4.2 Verification test programs A l l test programs are conducted on a cubical domain with an unstructured mesh, with 1166 control volumes. In the reconstruction scheme, integration over surface triangles and tetrahedral volumes is carried out. Using the IntegTest program, the Gauss weights and Gauss points were checked with up to 4 t h order accuracy, as well as confirming the integration itself. The results from the IntegTest program confirms that mesh class is able to supply appropriate Gauss weights and Gauss point locations to perform the integration. The test functions are monomials, therefore the numerical integration using Gauss quadrature is exact. The computed solution is then compared with the exact monomial integrals. Integration is checked for each order of accuracy (2 n d , 3 r d , and 4 t h) at a tolerance limit of l.Oe"13. A l l monomial integrations are passed. In the MeshTest program verifies that mesh pre-processing done correctly. The exact volume integration of monomials is compared with the computed monomials. The number of monomials required for the 3 r d order accuracy is 10, and for the 4 t h order accuracy is 20. The monomials are computed as a pre-processing step in the mesh class. The exact volume integral of monomials is found from formulas given in section A.4, and for higher order monomials Maple™ software is used to compute the exact integral, using the formula given in Eq. (A.7). Gauss quadrature is used to compute integrals of monomials to high order accuracy. The tolerance limit between the exact and the computed monomials is l.Oe"13. For each order of accuracy the respective number of monomials pass the test. In the MomentTest program the moments of control volumes about the cell centre (the reference location) are checked. The moments are computed for the entire domain and the results are compared with the generalization of the parallel axis theorem, given in section A.3. The program was run once again for each order of accuracy. The largest error is smaller than ±1 .Oe"14. After checking the components of the reconstruction in three-dimensions, the order of accuracy of the whole reconstruction routine is checked up to 4 t h order accuracy by the ReconTest program. A set of test functions, corresponding to monomials of order 3 or lower, is integrated to 17 give control volume averages over an entire mesh. This data is then reconstructed for each control volume and checked by evaluating the monomial at a test point near the C V reference point on each cell. The gradients of monomials are tested in a similar way. The comparison between the reconstructed value and the exact value gave a maximum error smaller than 1 .Oe" . A sinusoidal test function, given by Eq. (4.4a) below, is used to test the interior control volume reconstruction. Such a function cannot be represented exactly using monomials. The L2 error norm between the reconstructed and exact (Eq. (4.4a)) is plotted against number of control volumes as shown in Figure 4-1. The last three points were used to draw a least-square fit line, and the gradient of this line for 2 n d , 3 r d , and 4 t h order methods gave computed order of accuracy values of 1.7, 2.9, and 3.8 respectively. The reconstruction scheme performs well to expected order of accuracy. Figure 4-1 Order of accuracy of the reconstruction of function in Eq. (4.4a) for the interior control volumes 4.2.1 Boundary Condition treatment For the boundary cells the solutions must be constrained to enforce the Dirichlet and Neumann boundary conditions. The boundary constraints are enforced at the boundary Gauss integration 18 points. On the adjacent control volume, the constraints are added to the least-squares reconstruction. The test program ConstraintTest computes the reconstructed values with and without the boundary constraints for each order of accuracy. The tolerance limit, once again, is l.Oe"13. In all cases the unconstrained reconstruction is inaccurate, as expected. Three solutions used to initialize the domain are given by Eq. (4.4a-c). The use of three solutions also tests the Field Quant class in handling more than one solution in the three dimensional runs. T(x,y, z) = sin(7rx)sin(7ry)sin(7iz) T(x, y, z) = sin(27tx) sin(27ty) sin(27tz) T(x, y, z) = sin(47rx) sin(47ty) sin(47iz) (4.4a) (4.4b) (4.4c) Figure 4-2 The boundary face directions and labels for the cubic domain The error in constrained reconstruction is evaluated. Table 4-2 and Table 4-3 shows the number of errors reported for the two test cases when the tolerance limit is l.Oe"13. The first test case has Dirichlet conditions on all six boundary faces, and the second test case has one Neumann condition on the 'East' face and Dirichlet boundary conditions on all other faces (boundary faces can be identified in Figure 4-2). It can be seen that in Table 4-2 the high-order pre-processing (column one) results in more boundary Gauss points at which constraints must be satisfied. This results in having more constraints than degrees of freedom, and makes it impossible to satisfy all the constraints 19 simultaneously. Therefore when all the six boundary conditions on a cube are Dirichlet then 2 n and 4 t h order accurate reconstruction with the same order accurate moments leads to no errors in boundary reconstruction. However the third order scheme produces some errors. When one boundary is Neumann at 'West' face and the rest are Dirichlet in Table 4-3, then once again the 2 n d and 4 t h orders perform equally well for all three solutions. Table 4-2 All Dirichlet constrained boundary condition test Order of mesh Moments Order of reconstruction Dirichlet wall Dirichlet 'East' face Dirichlet 'West' face 2 2 0 0 0 3 2 46 18 16 3 10 2 7 4 2 92 36 32 3 0 0 0 4 0 0 0 Table 4-3 Dirichlet-Neumann constrained boundary condition test Order of mesh Moments Order of reconstruction Dirichlet wall Dirichlet 'East' face Neumann 'West' face 2 2 0 0 0 3 2 39 18 110 3 10 2 0 4 2 82 36 172 3 0 0 56 4 0 0 0 4.3 Conclusion The k-exact reconstruction implementation is discussed, together with various means of testing it. It appears that Gauss-quadrature can provide a very accurate monomial integral and flux integral. The least-squares technique worked well in three dimensions and the boundary constrained reconstruction achieved acceptable level of accuracy, apart from the 3 r d order accurate scheme. 20 5 M o d e l V e r i f i c a t i o n The numerical algorithm used to solve the given problem consists of a spatial discretization and an iterative solver. It is essential that the numerical algorithm used to solve the flow equations be reliable, accurate and efficient. The efficiency is measured by the computational effort required to solve a problem to a given numerical convergence level. Therefore, the convergence rate of the iterative solver measures the computational effort. The accuracy of the numerically converged solution depends on spatial discretization, which introduces flux evaluation errors and reconstruction errors. In this chapter, the accuracy of the discretized solution and computational efficiency are assessed. In section 5.1, the order of accuracy in evaluating the flux integral is calculated. In section 5.2, numerically converged steady state solution of the advection-diffusion and Poisson problems are used to assess both the accuracy of the method used and the efficiency of the three-dimensional solver. The accuracy of each order method is computed by calculating the error in computation with respect to the exact analytical solutions. The rate of convergence of the Newton-Krylov iterative scheme is graphed to establish the computational effort required to achieve the converged solution. A flux test is carried out over a control volume for the model problems to verify that the flux integration is done accurately. The error in surface flux evaluation is computed by comparing it with the volume integral of flux divergence (as shown in Eq. (5.1)) evaluated to 6 t h order accuracy. where the F, G, and H are the flux terms. For an advection-diffusion and a Poisson problem, the exact analytical solutions are used to compute the flux divergence. The initial solution is set to the exact solution. The error in the flux integration test is used to calculate the global error norms, the L I , L2 and Linf (L°°), as defined in Eq. (A.9) to Eq. ( A . l 1) in Appendix A . Out of all the error norms, the L2 norm is chosen to measure the error and used in computing the order of accuracy of the method. The error is proportional to the size of control volume to the power k, where k is the order of accuracy. Since the size of a control volume is proportional to the cube root of the number of cells, the error can be expressed in terms of the number of cells occupying 5.1 Flux integral test (5.1) 21 the domain. Therefore, the order of accuracy can be determined by evaluating the slope of a log-log plot of the L2 error norms versus the number of control volumes, which then is multiplied by three. (The formula is given in Eq. (A. 12) in Appendix A). 5.1.1 Advection-diffusion Consider an advection-diffusion problem over a three-dimensional cube. Assuming that the velocity field is constant, and the dependent variable is the scalar temperature T with a diffusion coefficient a. The advection-diffusion equation is solved over the domain with Dirichlet and Neumann boundary conditions, as shown below. 8T 3T 3T u—— + v — + w—- = a| dx dy dz rd2T d2T d2T^ dx2 3y2 dz2 (5.2) 0 < x < a, 0<y<b, 0<z<c T(0,y,z) = s i n f ^ S m | V b , Tx(a,y,z) = 0 T(x,0,z) = 0 T(x,b,z) = 0 T(x,y,0) = 0 T(x,y,c) = 0 V c ; These equations are solved to obtain the exact solution given by Eq. (A. 18), and the flux divergence is given by Eq. (A.22). The domain size is such that 'a', 'b', and 'c ' are all one, and the'm' and 'n ' values in the sine function of the first boundary condition equation are set to one. The diffusion coefficient for the flux test is set to 0.001. The constant velocity is just constant x-component of the velocity, and it is also of magnitude one. The error norms were calculated for each order of accuracy method and are shown in Table 5-1 to Table 5-3. The plot of the L2 error norms are shown in Figure 5-1. The last three points for each order of accuracy method are used to draw a least-square fit linear curve. The last three points are expected to lie close to the grid independent region. The slope of curves represents the accuracy of the advection-diffusion flux for each order of accuracy method used. The linear equation and the R 2 value of the fitted line are shown in Figure 5-1. From the gradients of the curves it is clear that the 4 t h order has a steeper gradient than the 3 r d order, which, in turn, has a 22 steeper gradient than the 2 order method. The second order method behaves like 3.26 order, 3 r order method behaves like 4.07 order, and 4 t h order method behaves like 4.59 order. The flux integration seems to be working very well. Table 5-1 The calculated error norms and order of accuracy for advection-diffusion flux text for 2nd order accurate method Number of control volumes LI L2 Linf Order 280 6.4287E-04 8.1713E-04 2.6734E-03 1166 1.0349E-04 1.7014E-04 1.4166E-03 3.299924 2112 4.5679E-05 6.5627E-05 2.6730E-04 4.810875 4409 2.0060E-05 3.8282E-05 3.0421E-04 2.197023 6986 1.1313E-05 2.3429E-05 2.8402E-04 3.200484 11404 5.4223E-06 9.1606E-06 8.6344E-05 5.748642 16523 3.4483E-06 5.7102E-06 5.0388E-05 3.824223 23025 2.4524E-06 4.2775E-06 3.2843E-05 2.611747 Table 5-2 The calculated error norms and order of accuracy for advection-diffusion flux text for 3 rd order accurate method Number of control volumes LI L2 Linf Order 280 2.2534E-04 3.0047E-04 7.5830E-04 1166 2.3964E-05 3.2082E-05 1.1908E-04 4.704496 2112 9.5510E-06 1.3709E-05 8.6815E-05 4.293762 4409 3.2290E-06 5.3131E-06 3.4079E-05 3.863536 6986 1.5513E-06 2.5029E-06 2.1604E-05 4.906396 11404 6.4325E-07 9.4960E-07 5.9050E-06 5.932936 16523 3.5250E-07 5.2859E-07 3.9751E-06 4.739888 23025 2.3233E-07 3.6758E-07 4.1182E-06 3.284215 Table 5-3 The calculated error norms and order of accuracy for advection-diffusion flux text for 4 order accurate method Number of control volumes LI L2 Linf Order 280 6.6011E-05 1.3239E-04 7.3711E-04 1166 4.5653E-06 7.1024E-06 4.1358E-05 6.151951 2112 1.8449E-06 3.6429E-06 3.1269E-05 3.371608 4409 4.3568E-07 9.3191E-07 8.8577E-06 5.556865 6986 2.0369E-07 5.1491E-07 7.0392E-06 3.86675 11404 8.1128E-08 2.3523E-07 3.1393E-06 4.796075 16523 4.0211E-08 1.1217E-07 1.6559E-06 5.991274 23025 2.4670E-08 8.0851E-08 1.6337E-06 2.960333 23 Figure 5-1 The L2 error norm for the advection-diffusion flux test for various mesh sizes 5.1.2 Poisson The Poisson problem over a three-dimensional cube is considered. The dependent variable is the scalar temperature, T. The advection-diffusion equation is solved over the domain solely with Dirichlet boundary conditions as shown below. The values of 'a', 'b', 'c ' , 'm' , and 'n ' are all one, as explained in the previous section. V 2 T = f(x,y,z) (5.7) 0 < x < a, 0<y<b, 0<z<c T(0,y,z) = sin T(a,y,z) = 0 T(x,0,z) = 0 T(x,b,z) = 0 T(x,y,0) = 0 T(x,y,c) = 0 n f ) i n nrcz 24 where the source term is f(x,y,z) = sin nr. sin 'rrm H Ml (5.8) The 'r' value in the source term is set to one. The final solution is given by Eq. (A.39), and the flux divergence used in flux test is given by Eq. (A.41). Figure 5-2 The L2 error norm for the Poisson flux test for various mesh sizes The error norms were calculated for each order method and they are shown in Table 5-1 to Table 5-3. The plot of the L2 error norms is shown in Figure 5-2. The last three points for each order of accuracy method are used to draw a least-square fit linear curve. The linear equation and the R 2 value of the fitted line are shown. Unlike the advection-diffusion case, the slopes of the 3 r d and 4 t h order methods are almost the same and the data points are much close together. It appears that the 4 t h order flux integral does not create any significant reduction in error for the Poisson problem. The 2 n d method behaves like 2.99 order, the 3 r d order method behaves like 2.94, and the 4 t h order method behaves like 2.94 order. 25 Table 5-4 The calculated error norms and order of accuracy for Poisson flux test for the 2nd order accurate method Number of control volumes LI L2 Linf Order 280 4.55E-03 8.01E-03 3.88E-02 1166 1.35E-03 2.24E-03 1.23E-02 2.680487 2112 6.41E-04 1.00E-03 5.27E-03 4.06857 4409 3.55E-04 6.08E-04 4.55E-03 2.032273 6986 2.17E-04 3.93E-04 3.54E-03 2.84362 11404 1.28E-04 2.43E-04 2.03E-03 2.935679 16523 8.75E-05 1.65E-04 1.79E-03 3.149341 23025 6.59E-05 1.21E-04 1.45E-03 2.803133 Table 5-5 The calculated error norms and order of accuracy for Poisson flux test for the 3 rd order accurate method Number of control volumes LI L2 Linf Order 280 2.96E-03 4.55E-03 1.50E-02 1166 7.17E-04 1.02E-03 4.06E-03 3.146163 2112 3.25E-04 4.31E-04 2.40E-03 4.343603 4409 1.64E-04 2.30E-04 1.58E-03 2.56188 6986 1.01E-04 1.61E-04 1.87E-03 2.315774 11404 5.37E-05 7.06E-05 3.86E-04 5.047794 16523 3.72E-05 4.91E-05 2.04E-04 2.933207 23025 2.68E-05 3.55E-05 3.13E-04 2.950528 Table 5-6 The calculated error norms and order of accuracy for Poisson flux test for the 4 order accurate method Number of control volumes LI L2 Linf Order 280 2.23E-03 3.17E-03 1.06E-02 1166 5.99E-04 8.58E-04 2.77E-03 2.7474 2112 2.72E-04 3.51E-04 1.00E-03 4.516881 4409 1.32E-04 1.75E-04 6.07E-04 2.843967 6986 8.15E-05 1.08E-04 5.51E-04 3.14577 11404 4.68E-05 6.18E-05 2.67E-04 3.410046 16523 3.30E-05 4.40E-05 1.63E-04 2.749631 23025 2.35E-05 3.10E-05 1.24E-04 3.15219 26 5.2 Steady state solution The advection-diffusion and Poisson problems are again selected for the steady state verification. In all the computational runs, the residual error norms were less than l.Oe" . The residual here defines the difference in flux integral from one iteration to the next. The initial solution value chosen was equal to (1-x). The error is evaluated using the difference between the exact solution and the converged numerical solution. The diffusion coefficient for the advection-diffusion problem was 0.001 and the velocity field is only in the x-direction with a magnitude of one. In all cases, where appropriate, local time stepping was chosen without any under-relaxation. 5.2.1 Advection-diffusion results The cross-sectional solution, for the advection-diffusion problem at various positions along the x-direction of the geometry, is shown in Figure 5-4. The solution is very smooth and decreases in magnitude along the x-direction. -6 J 1 1 1 1 ! 1 0 1 2 3 4 5 6 Log 10 of number of control volumes Figure 5-3 The L2 error norm for the advection-diffusion problem at steady state for various mesh sizes 27 x \ / / \ \ \ \\\\ J -0.000714 O.I« 0.286 0.571 0.714 0.857 1 ' (a) -0000714 0.142 0 286 0. 428 0.571 0 714 0.857 (c) HI/ WW / I I ! ( ( V \ \ \ ))•!} -0.000714 0.^ 42 (b) 0.571 0.714 0.857 1 T -oWm 0.142 0.285 0.428 0.571 0.714 0.857 (d) Figure 5-4 The sub figures (a)-(d) shows the solutions to the advection-diffusion problem along the geometry at x=0, 0.3, 0.6 and 1. 28 Table 5-7 The calculated error norms and order of accuracy for advection-diffusion problem for the 2' order accurate method Number of control volumes LI L2 Linf Order 280 0.038915 0.053234 0.194567 1166 0.014205 0.021685 0.20923 1.888644 2112 0.008886 0.013002 0.17607 2.583372 4409 0.005105 0.008113 0.090074 1.922037 6986 0.003943 0.006394 0.084893 1.552309 11404 0.00258 0.004086 0.047984 2.742009 23025 0.001632 0.002662 0.038703 1.828511 Table 5-8 The calculated error norms and order of accuracy for advection-diffusion problem for the 3 rd order accurate method Number of control volumes LI L2 Linf Order 280' 0.010911 0.014771 0.056298 1166 0.003214 0.004309 0.019803 2.59104 2112 0.001607 0.00218 0.019316 3.44124 4409 0.000775 0.001086 0.014465 2.839161 6986 0.000502 0.000705 0.007428 2.816198 11404 0.000286 0.000402 0.003326 3.445605 23025 0.000144 0.000208 0.002108 2.816537 Table 5-9 The calculated error norms and order of accuracy for advection-diffusion problem for the 4 order accurate method * Number of control volumes LI L2 Linf Order 280 0.004062 0.006167 0.034796 1166 0.000674 0.000915 0.006093 4.013034 2112 0.000352 0.000509 0.004252 2.96263 4409 0.000136 0.000203 0.001681 3.753562 6986 8.26E-05 0.000132 0.001884 2.812106 11404 4.32E-05 7.60E-05 0.001929 3.361907 23025 1.85E-05 3.65E-05 0.000835 3.130634 29 Table 5-7 to Table 5-9 show the computed orders of accuracy of the 2 n d , 3 r d and 4 t h order methods, and Figure 5-3 summarizes the L2 error norm. From these results it is clear that 2 n d and 3 r d order methods behave better than expected, while the 4 t h order method behaves only slightly better than a 3 r d order scheme. The solution error after numerical convergence is graphically presented in Figure 5-5 to Figure 5-7 for 2 n d to 4 t h order accurate solution. This is the local error between the computed and the exact. The darker grey triangles show the largest local error in the domain. For the 2 n d order method, high magnitude errors occur frequently at the block corners. The 3 r d order accurate solution produces the significant magnitude of error, as can be seen in Figure 5-6, on the outlet face. Similarly 4 t h order accurate solution produces the largest error on this outlet face, and compared with the 3 r d order method, the 4 t h order method has errors of high magnitude on this outlet face. On the outlet face, Neumann boundary conditions are enforced. It appears that the reason 4 t h order method does not produce 4 t h order accurate solution is that the error introduced at the outlet boundary does not decrease sufficiently. Hence, there may be a boundary condition evaluation or a boundary condition enforcement problem that requires further investigation. T 0 1 0.00425 0.00851 0.0128 0.017 0.0213 0.0255 0 0 2 9 8 Figure 5-5 Error between the exact and 2" order accurate solution, after numerical convergence of the advection-diffusion problem 30 Figure 5-6 Error between the exact and 3 r order accurate solution, after numerical convergence of the advection-diffusion problem Figure 5-7 Error between the exact and 4 order accurate solution, after numerical convergence of the advection-diffusion problem 31 The computational efficiency of the methods implemented by the ANSLib code is assessed by plotting the L2 error norm versus the computational time taken to achieve the numerical converged solution, as shown in Figure 5-8. The data points represent various meshes (defined in Table 5-7) of increasing number of control volumes as the computational time increases. For a given level of accuracy corresponding to an L2-norm value, it is clear that a 4 t h order method achieves a numerically converged solution on a much coarser mesh than either 3 r d or 2 n d order methods. A 4 t h order method progressively consumes larger amounts of time on finer meshes. However, for a given level of accuracy the 4 t h order method takes less time than either of the lower order methods. For the advection-diffusion problem there exists an exponential reduction in accuracy for increasing computational time. This implies that the lower order solution methods would need very fine meshes and require long computational times in order to reach a given level of high order accuracy. However, the low order methods may reach asymptotic convergence faster than the high order methods. Figure 5-8 The computational time taken for the steady-state solution of the advection-diffusion problem for various meshes 32 The rate of convergence for each order method for various mesh sizes is shown in Figure 5-9 to Figure 5-11. For all methods up to 23025 control volumes, the number of iterations needed for convergence is less than 15. From Figure 5-11 and Figure 5-8 it is apparent that the reason 4th order method takes longer is that the cost per iteration is high. In the 4th order method the number of monomials that should be evaluated is almost twice as many as are needed for the 3 r d order method leading to a larger least-squares matrix. In addition, there are 16 Gauss points per control volume on which the fluxes must be calculated, which add to the computational time. The 3rd and 4th order methods converge with a slightly smaller number of iterations, compared to the 2n d order method. For the finest mesh all schemes show deviation from the usual quicker rate of convergence. The solution update become smaller after the residual decrease beyond 1 .Oe"6. It seems some error term gets build-up up to this point and starts to affect the solution update. The matrix vector product of the Jacobian and the solution update vectors in GMRES algorithm may be badly conditioned. Residual L2 Norm of AdvDiff, a=2 Figure 5-9 The L2 residual norm for the advection-diffusion problem for the 2" order accurate method 33 Residual L2 Norm of AdvDiff, a=3 Iterations Figure 5-10 The L2 residual norm for the advection-diffusion problem for the 3 rd order accurate method Residual L2 Norm of AdvDiff, a=4 Iterations Figure 5-11 The L2 residual norm for the advection-diffusion problem for the 4th order accurate method 34 2nd Order on 27021 Cvs, with different krylov subspaces, m=30, m=60 Number of iterations Figure 5-12 The L2 residual norm for the advection-diffusion problem for the 2nd order accurate method with Krylov subspace vectors sizes of 30 and 60 For the finest mesh with 27021 control volumes advection-diffusion problem was run again with 30 and 60 Krylov subspaces. The convergence history is shown in Figure 5-12. It can be seen that a larger Krylov subspace is needed for the finer mesh to improve the convergence. However, larger Krylov subspace vector means larger computational time per iteration. This may lead to inefficient computation. In this case, however, there is a drastic reduction is number of iterations, over 20, and the computational time measured showed 1 hour reduction in computation time for the case with 60 Krylov subspace vector. Therefore, larger Krylov subspace leads to an improvement in efficiency for finer meshes. 5.2.2 Poisson results The solution at three different cross-sections along the geometry in the x-axis is shown in Figure 5-14. The solution at the 'East' face, where x=0 (refer to Figure 4-2), quickly reduces in magnitude leaving the 'West' face, where x=1.0, with only a small magnitude of solution. 35 Figure 5-13 The L2 error norm for the Poisson problem at steady state for various mesh sizes The order of accuracy of the 2 n d , 3 r d and 4 t h order methods are shown in Table 5-10 to Table 5-12, and Figure 5-13 summarizes the order of accuracy of these methods. The table makes clear that the gradient used to compute the order of accuracy fluctuates a little, but it is reduced in magnitude as the mesh becomes finer. In consideration of the last three points, least-square trend lines are drawn to show the order of accuracy of methods. Once again, the last three points are assumed to be in the asymptotic convergence region, where an expected order of accuracy is reached. It can be seen that the 2 n d order method behaves like 2.03 order accuracy, the 3 r d order method behaves like 3.21 order accuracy, and the 4 t h order scheme behaves like 3.54 order accuracy. 35 r / f \\\ n I Z PC T . -0:00377 0 14 o i » 3 0 426 ' 'osT (a) 0.713 0.857 1 -0W37W ' 0.2041 0308 0 412 ' O.Slt 0 618 '0.723 0.*27 BAH (C) T -o'ooWo/l ' 'oio4 o^SoT 0.41* oiii 'OJM am oizt o&ti (b) Figure 5-14 The sub figures (a)-(c) show the solutions to the Poisson problem along the geometry at x=0, 0.3, and 0.6. 36 Table 5-10 The calculated error norms and order of accuracy for Poisson problem for the 2" order accurate method Number of control volumes LI L2 Linf Order 280 0.014209 0.026896 0.217003 1166 0.006177 0.013538 0.153936 1.443619 2112 0.003935 0.008263 0.098601 2.493532 4409 0.002783 0.006265 0.106814 1.128384 6986 0.002304 0.005034 0.058556 1.42595 11404 0.001514 0.003249 0.048518 2.680852 23025 0.001014 0.002217 0.038729 1.63136 Table 5-11 The calculated error norms and order of accuracy for Poisson problem for the 3 rd order accurate method Number of control volumes LI L2 Linf Order 280 0.005436 0.009923 0.063499 1166 0.00158 0.003033 0.037175 2.492926 2112 0.000823 0.001513 0.022905 3.511748 4409 0.000533 0.001043 0.01619 1.514797 6986 0.000364 0.000757 0.018307 2.090403 11404 0.000215 0.000388 0.009702 4.089946 23025 0.000116 0.000208 0.005364 2.657524 Table 5-12 The calculated error norms and order of accuracy for Poisson problem for the 4 order accurate method Number of control volumes LI L2 Linf Order 280 0.002108 0.003856 0.030735 1166 0.000327 0.000706 0.010416 3.570896 2112 0.000164 0.000355 0.008115 3.464576 4409 9.12E-05 0.000195 0.004366 2.441283 6986 4.47E-05 0.000115 0.003947 3.439482 11404 2.24E-05 5.12E-05 0.001121 4.969616 23025 1.15E-05 2.75E-05 0.001792 2.65582 37 The error between the numerically converged solution and the exact solution are graphically presented in Figure 5-15 to Figure 5-17 for 2 n d to 4 t h order accurate solutions. The Figure 5-15 show local error is clustered around the 'East' face, and the largest magnitude of error is darker grey colour. In all figures, the maximum error is clustered to the corners of the 'East' face. Unlike the advection-diffusion problem the 2 n d , 3 r d and 4 t h order accurate method produces errors at the 'East' face with the error decreasing in magnitude and occurrence. Considering that the same boundary condition is applied to the 'East' face in advection-diffusion case, but the decay in solution here is larger for the Poisson's case, the solution near 'East' face boundary in Poisson's case would have large gradients in the x-direction. This may lay some strain on the accuracy of the boundary constraints, while at the 'West' face the solution and its gradient are small therefore the solution error is small. The reason why the 4 t h order method has large errors at the corners (Figure 5-17) leading to 3 r d order behaviour may be due to some incompatible boundary conditions at this block corner, where two differently valued boundary conditions are constrained in the least-square problem. u 0 1 0.00315 0.00631 00*0946 0.0126 0.0158 0.0189 0.0221 Figure 5-15 Error between the exact and 2nd order accurate solution, after numerical convergence for the Poisson problem 38 0.000641 0.00128 0.00192 0.00256 0.0032 0 0 0 3 8 4 0.00449 Figure 5-16 Error between the exact and 3 r d order accurate solution, after numerical convergence for the Poisson problem Mr I f ft* 0.000263 0.000527 0.00079 * 0 . 0 0 1 0 5 0.00132 0.00158 0.00184 Figure 5-17 Error between the exact and 4th order accurate solution, after numerical convergence for the Poisson problem 39 The computational efficiency of the methods used in solving the Poisson problem is shown in Figure 5-18. The plot points correspond to the meshes shown in Table 5-10. Once again, in comparison to the lower order methods the 4 t h order method achieves higher order solution on a coarser grid with lesser computational time for a given level of error norm. Figure 5-18 The computational time taken for the steady-state solution of the Poisson problem for various meshes The rate of convergence of the solution is shown in Figure 5-19 to Figure 5-21 for 2 n d , 3 r d , and 4 t h order methods. Unlike the advection-diffusion problem the number of iterations needed to reach solution convergence increases with the number of control volumes. The 4 t h order method requires more iterations than the 3 r d order method, which in turn require more iterations than the 2 n d order method. This means the rate of convergence decreases with order of accuracy. The 4 t h order method taking more iterations may be attributed to the 4 t h order Jacobian matrix having more complicated structure than the 3 r d order Jacobian matrix. However, in all cases the number of iterations needed for converged solution is small. Figure 5-20 The L2 residual norm for the Poisson problem for the 3rd order accurate method 41 Iterations Figure 5-21 The L2 residual norm for the Poisson problem for the 4th order accurate method 5.3 Conclusion The flux test results show that the surface flux evaluations performed using the Flux class and the Integrate program are good enough to carry out a steady state analysis. For a given error norm value the 4 t h order method achieves convergence on a coarse mesh, while finer meshes are required for the 3 r d and 2 n d order methods. Also, the time taken to achieve a given level of accuracy for the 4 t h order method is considerably smaller than for the 2 n d order method, while the 3 r d order method lies somewhere in between. In the advection-diffusion problem for 3 r d and 4 t h order methods showed the largest local error localized at the 'West' face or the outlet, while for the Poisson problem it was at the 'East' face or the inlet. In both cases for 2 n d order method, the errors cluster at the four corners. The number of overall iterations needed for numerical convergence is small for both the advection-diffusion problem and the Poisson problem. For finer meshes, larger Krylov subspace vector is necessary. The 4 t h order method consumes the largest amount of time per iteration owing to a large number of Taylor expansion terms in the least-squares matrix, and the presence of a large number of Gauss points on which the fluxes are evaluated. 42 6 C o n c l u s i o n s a n d F u t u r e W o r k The program classes in ANSLib code use abstract data types, which play a central role in object-oriented programming in allowing what is being done to be separate from how it is done. ANSLib generic solver separates the implementation of physics from numerics, so that the end user will not need to modify the code when new physics are introduced. The flow chart in Figure 2-1 showed that flux integration and time advance are separate entities, and the dimensional changes to the program became independent as well. This allowed code development and verification to be carried out in three dimensions easily and quickly. Various test cases are run in order to verify the solution reconstruction, flux integration and iterative time advancement to steady state in three dimensions. The flux integration testing is a crude approximation to test the exact flux integral, but it serves as a good means of determining inconsistencies in flux integration, using Gauss quadrature. It becomes apparent that Gauss quadrature is a very accurate means of integrating monomial functions, and to resolve 2n d and 4th order Dirichlet and Neumann boundary conditions. The order of accuracy of the 4th order does come very close to being 4th order accurate in steady state analysis. Also it is clear that the 4th order accurate method is the best to obtain a numerically converged solution for a given level of error norm, since the mesh sizes required and computational time used in achieving the results, are superior to 2nd or 3rd order methods. The final results, from flux integration and steady state runs, allows the three dimensional code to be fully verified for solving linear problems. The code verification must be done to solve a Navier-Stokes problem, then it may be proved that problems with sufficient complexity can be solved. For example, the exact analytical unsteady three dimensional Navier-Stokes solutions given by Ethier and Steinman23 present a good benchmark problem to solve. In order to tackle non-linear problems it is the solver classes and the definition of fluxes in the physics description program that are expected to be modified. The Newton-Krylov solver may needs to be upgraded with preconditioning, since preliminary runs on a three-dimensional Navier-Stokes problem showed poor numerical convergence. Improvements in the area of accuracy may be made with the three dimensional data-dependent reconstruction, where the data-dependent weight no longer serve as a means of determining discontinuities in the solution domain, but instead, resolves sharp gradients in incompressible flows. Accuracy may also be improved by creating a layer of small size control volumes at all boundaries. 43 7 B i b l i o g r a p h y [I] Harten, A. , "Uniformly High Order Accurate Essentially Non-oscillatory Schemes, III", Journal of Computational Physics, 131:3-47 (1997) [2] Saad, Y . , Schultz, M . , "GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems", S IAM Journal on Scientific and Statistical Computing, 7(3): 856-869 (1986). [3] Zimmermann, T., Dubois-Pelerin, Y . , Bomme, P., "Object-oriented finite element programming: I. Governing Principles", Computer Methods in applied Mechanics and Engineering, 98:291-303 (1992). [4] Liu, J-L., Lin, I-J., Shih, M-Z. , Chen, R-R., Hsieh, M - C , "Object-oriented programming of adaptive finite element and finite volume methods", Applied Numerical Mathematics, 21:439-467 (1996). [5] Ollivier-Gooch, C , " A toolkit for numerical simulation of PDEs: I. Fundamentals of generic finite-volume simulation", Computer Methods in Applied Mechanics and Engineering, 192:1147-1175 (2003). [6] Boivin, C , "Infrastructure for Solving Generic Multiphysics problems", PhD thesis, Department of Mechanical Engineering, University of British Columbia, Vancouver, Canada. [7] G R U M M P : Two dimensional and three dimensional unstructured grid generator. Available from website <http://tetra.mech.ubc.ca/GRUMMP/index.html> [8] Boris, J.P., and Book, D.L., "Flux corrected transport: I, SHASTA, a fluid transport algorithm that works", Journal of Computational Physics, 11:38-69 (1973) [9] Van Leer, B. , "Towards the ultimate conservative difference scheme V . A second order sequel to Godunov's method", Journal of Computatioinal Physics, 32:101-136 (1979) [10] Osher, S., "Riemann solvers, the entropy condition and difference approximations", S IAM Journal of Numerical Analysis, 21(2):217-235 (1984) [II] Roe, P.L., ""Approximate Reimann solvers, parameter vectors and difference schemes", Journal of Computational Physics, 43(2):357-372 (1981) [12] Deconinck, H. , Paillere, H. , Struijs, R., Roe, P.L., "Multidimensional upwind schemes based on fluctuation-splitting of conservation laws", Computational Mechanics, 11:323-340 (1993) [13] Venkatakrishnan, V . , " A Perspective on Unstructured Grid Flow Solvers", A I A A Journal, 34(3):533 (1996) [14] Barth, T.J., Jespersen, D.C., "The design and application of upwind schemes on unstructured meshes", A I A A Paper 89-0366 (1989) 44 [15] Venkatakrishnan, V. , "On the accuracy of limiters and convergence to steady state solutions", A I A A Paper 93-0880 (1993) [16] Harten, A. , "High resolution schemes for hyperbolic conservation laws", Journal of Computational Physics, 49:357 (1983) [17] E, W., Shu, C-W., " A Numerical Resolution Study of High Order Essentially Non-oscillatory Schemes Applied to Incompressible Flow", Journal of Computational Physics, 110:39-46(1994) [18] Barth, T.J., Fredrickson, P.O., "Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction", A I A A Paper 90-0013 (1990) [19] Ollivier-Gooch, C.F., "Quasi-ENO schemes for unstructured meshes based on unlimited data-dependent least-squares reconstruction", Journal of Computational Physics, 133:6-17 (1997) [20] Ollivier-Gooch, C , Altena, M . V . , " A High-Order-Accurate Unstructured Mesh Finite-Volume Scheme for the Advection-Diffusion Equation", Journal of Computational Physics, 181: 729-752 (2002) [21] Venkatakrishnan, V . , Mavriplis, D.J., "Implicit solvers for unstructured meshes", Journal of computational physics, 105(1):83-91 (1993) [22] Ethier, C. R., Steinman, D. A. , "Exact fully 3D Navier-Stokes solutions for benchmarking", International Journal for Numerical Methods in Fluids, 19(5):369-375 (1994) A A p p e n d i x A A.1 Gauss quadrature for a triangle Table A-l Gauss points and weight for a triangle n Order Points & & % Wj 1 Linear (1 or 2 ) 1 1/3 1/3 1/3 1 1 1/2 1/2 0 1/3 3 Quadratic (3) 2 0 1/2 1/2 1/3 3 1/2 0 1/2 1/3 1 1/3 1/3 1/3 -9/16 4 Cubic (4) 2 0.6 0.2 0.2 9/16 3 0.2 0.6 0.2 9/16 4 0.2 0.2 0.6 9/16 1 A1 B1 B1 C1 2 B1 A1 B1 C1 6 Quintic (6) 3 B1 B1 A1 C1 4 A2 B2 B2 C2 5 B2 A2 B2 C2 6 B2 B2 A2 C2 A l 0.8168475730 CI 0.1099517437 BI 0.0915762135 C 2 0.2233815897 A 2 0.1081030182 B 2 0.4459484909 A.2 Gauss quadrature for a tetrahedron Table A-2 Gauss points and weight for a tetrahedron n Order Points $2 $ 3 $4 Wj 1 Linear 1 1/4 1/4 1/4 1/4 1 1 A1 A2 A2 A2 1/4 3 Quadratic 2 A2 A1 A2 A2 1/4 3 A2 A2 A1 A2 1/4 4 A2 A2 A2 A1 1/4 1 1/4 1/4 1/4 1/4 -4/5 4 Cubic 2 1/3 1/6 1/6 1/6 9/20 3 1/6 1/3 1/6 1/6 9/20 4 1/6 1/6 1/3 1/6 9/20 5 1/6 1/6 1/6 1/3 9/20 5+3V5 20 5 - V 5 20 46 A.3 Parallel axis theorem The total moment of x is the sum of moment about the centroid plus the contribution from points around the centroid, for each control volume, n. Some formulae for x and xy moments are shown here. x = TrZ V nk+jx .dv ] (A.1) »T n=l x y = 7 7 - X v n [ j x y - d v + x o j y - d v + y c H v + x c y c ] (A.2) • * T n=l A.4 Exact integration formula over tetrahedron Figure A-7-1 Tetrahedron control volume The centroid of a tetrahedron is located at point (0,0,0), hence the sum of the global coordinates equals zero. x l + x2 + x3 + x4 = 0 y l + y2 + y3 + y4 = 0 (A.3) z l + z2 + z3 + z4 = 0 Then the integration formulas for x, x 2 , and xy are as follows Jx.dV = 0 (A.4) 4 7 f x 2 . d V = ( x l 2 + x 2 2 + x 3 2 + x 4 2 ) — ! 2 0 ( A . 5 ) J x 2 . d V = ( x l y l + x 2 y 2 + x 3 y 3 + x 4 y 4 ) V _ 2 0 ( A . 6 ) w h e r e V is the v o l u m e . I n g e n e r a l the i n t e g r a l o f p o l y n o m i a l i n t e r m s o f the v o l u m e c o o r d i n a t e s £ i , § 3 , a n d is fe^&.dV = ( 6 V ) ( a + b + C + d + 3)! A.5 Calculating order of accuracy from error norms T h e e r r o r is d e f i n e d as the d i f f e r e n c e b e t w e e n the e x a c t a n d the c o m p u t e d s o l u t i o n . E i = Uexact ( X i . y i > Z i ) - U c 0 m p ( X i . Y i » Z i ) T h e e r r o r n o r m s are d e f i n e d as L I = L 2 = Z E . 2 N ; L i n f = L ° o = m a x E j i I f the L 2 - n o r m is g i v e n f o r t w o m e s h e s as s h o w n i n T a b l e A - 3 . Table A-3 Sample table of L2 error norms for each mesh Mesh Control Volumes L2-norm Mi L2i M2 L22 ( A . 7 ) ( A . 8 ) ( A . 9 ) ( A . 10) ( A . 11) T h e n the o r d e r o f a c c u r a c y is c a l c u l a t e d b y the f o l l o w i n g e q u a t i o n . 48 log Order of accuracy, = 3-L2, log v M 2 y (A.12) A.6 Analytical solution methods A.6.1 Advection-diffusion 9T dT 3T u—— + v— + w—— = a| dx dy dz ^32T 3 2T 3 2 T A v d x 2 + d y 2 + d z 2 y (A.13) 0 < x < a, 0<y<b, 0<z<c T(0,y,z) = sin| T,(a,y,z) = 0 T(x,0,z) = 0 T(x,b,z) = 0 T(x,y,0) = 0 T(x,y,c) = 0 sin ^ c ; (A. 14) By the method of separation of variables assume the solution to be in the form T(x,y,z) = Q(x).R(y).S(z) (A.15) Substituting Eq. (A.15) to Eq. (A.13), and then dividing by Eq. (A.15) leads to an equation, which is separated into left and right side of the equality sign, based on the property that one side must be independent of a variable, while the other side is dependent on this variable. Since one side depends on the variable which could be x, y, or z, and the other side does not, each side must be equal to a constant. Using these constants (X, [i) and assuming that velocity components are constants and u=U, and v=w=0, the following three equations are formed. Q ' - - Q ' + | i Q = 0 a R * + ( X - H ) R = O (A.16) (A.17) S*-XS = 0 (A. 18) 49 These equations are solved to obtain the exact solution T(x, y,2) = [C, (eY|X - eY2X)+ eY2" ]sm^^-sin V c ; (A. 19) where c = (A.20) Jl,2 —" 2a 2 + + ,2a, I b J I c J (A.21) The divergence of fluxes for advection-diffusion problem is given as V(uT - aVT) = (UC, (7,3, - y2a2)+Uy2a2 )A + 2a7t2 (C, (a, - a2) + a 2 )A (A.22) where a is the diffusion coefficient, and A = sin(7ry) sin(Ttz) a,=eY'x a 2 = e** •Y2e Y2 Y,eYl - Y2eY2 (A.23) (A.24) (A.25) (A.26) 11,2 — " 2a (A.21) A.6.2 Poisson The Poisson problem over a three dimensional rectangular domain is described as follows. V 2 T = f(x,y,z) 0 < x < a, (A.28) 0<y<b, 0<z<c 50 T(0,y,z) = sin T(a,y,z) = 0 T(x,0,z) = 0 T(x,b,z) = 0 T(x,y,0) = 0 T(x,y,c) = 0 'rrmy^ f sin cr V riTtz (A.29) where the source term f(x,y,z) = sin sin sm V D J (A.30) Assume the solution can be written in the form T(x, y, z) = U(x, y, z) + V(x, y, z) (A.31) where U(x,y,z) is a particular solution to the Poisson equation and V(x,y,z) is the solution of the associated homogeneous equation, i.e. V 2 U = f(x,y,z) V 2 V = 0 Let U(x,y,z) = A sin M sin sin V a J I b J I c J A = --1 rrc 2 ( 2 + + l a ; \ I c J (A.32) (A.33) , then substituting to Eq. (4.12) gives constant A as (A.34) Since U(x,y,z) is obtained the solution of the Dirichlet problem must be solved V 2 V = 0 on Domain V = - U + T(x,y,z) on Boundary (A.35) 51 V(0,y,z) = sin V(a,y,z) = 0 V(x,0,z) = 0 V(x,b,z) = 0 V(x,y,0) = 0 V(x,b,c) = 0 ^mTcy^ sin ^n7tz^ ^ c ; (A.36) By the method of separation of variables assume the solution to be in the form V(x,y,z) = Q(x).R(y).S(z) (A.37) The method of separation of variable, as was done in previous section, leads to the equation, cosh(y) V(x ,y ,z ) = -sinh sinh(v) This gives the final solution 'yx^ + cosh sin 'nmy^ .sin I a J U J. I b J I c J (A.38) T(x,y,z) = where A sin 71x^1 cosh(y) a ) sinh(y) sinh + cosh fyxY sin .sin (xnz\ I a J I a J. I b J I c J (A.39) Y = . nm + n7t (A.40) The divergence of Poisson problem is simply made to zero by reversing the sign of the source term. ( V ( V T ) - f ) = 0 (A.41)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Development and verification of the generic three dimensional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Development and verification of the generic three dimensional finite volume solver Perera, M. K. Harsha 2003
pdf
Page Metadata
Item Metadata
Title | Development and verification of the generic three dimensional finite volume solver |
Creator |
Perera, M. K. Harsha |
Date Issued | 2003 |
Description | The generic finite volume solver, ANSLib, has been extended to three dimensions and used to verify the accurate computation of three-dimensional advection-diffusion and Poisson problems. A simple cubic domain has been selected as the domain of interest. Over this domain a steady state solution is computed for each model problem with 2nd, 3rd and 4th order accurate schemes. Gauss quadrature is used to evaluate the flux integral to 2nd, 3rd and 4th order accuracy. The flux scheme makes use of the centred scheme to model the viscous term and a simple upwind scheme to model the advective fluxes. In both cases the existing k-exact least-square reconstruction code is used to obtain the solution at each Gauss integration point to 2nd, 3rd and 4th order accuracy. Newton-Krylov implicit time advance is used to obtain the steady state solution, implemented using the GMRES algorithm. Prior to subjecting ANSLib code to physical model problems, the code is first tested for the correct enforcement of boundary constraints, accurate implementation of the solution reconstruction, and flux integral evaluation in three dimensions. Relevant changes are made to the code to obtain the desired behaviour. Analytical solutions for the advection-diffusion and Poisson equations are derived to satisfy Dirichlet and Neumann boundary conditions; these are used to compute the resulting discretized error in numerically converged solutions. Accuracy assessment is completed through a comparison of error norms for different mesh sizes for each test problem. Localized errors on the domain for each problem at steady state are plotted. The efficiency of the code to obtain a numerical converged solution is reported. |
Extent | 6672715 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-11-17 |
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.0091244 |
URI | http://hdl.handle.net/2429/15077 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2004-0006.pdf [ 6.36MB ]
- Metadata
- JSON: 831-1.0091244.json
- JSON-LD: 831-1.0091244-ld.json
- RDF/XML (Pretty): 831-1.0091244-rdf.xml
- RDF/JSON: 831-1.0091244-rdf.json
- Turtle: 831-1.0091244-turtle.txt
- N-Triples: 831-1.0091244-rdf-ntriples.txt
- Original Record: 831-1.0091244-source.json
- Full Text
- 831-1.0091244-fulltext.txt
- Citation
- 831-1.0091244.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-0091244/manifest