A MULTI-GRID METHOD FOR COMPUTATION OF FILM COOLING By Jian Ming Zhou B. Sc. (Mathematics) Fudan University Shanghai, China A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES INSTITUTE OF APPLIED MATHEMATICS DEPARTMENT OF MECHANICAL ENGINEERING DEPARTMENT OF MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 1990 © Jian Ming Zhou, 1990 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 The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract This thesis presents a multi-grid scheme applied to the solution of transport equations in turbulent flow associated with heat transfer. The multi-grid scheme is then applied to flow which occurs in the film cooling of turbine blades. The governing equations are discretized on a staggered grid with the hybrid differencing scheme. The momentum and continuity equations are solved by a nonlinear full multi-grid scheme with the S I M P L E algorithm as a relaxation smoother. The turbulence k — e equations and the thermal energy equation are solved on each grid without multi-grid correction. Observation shows that the multi-grid scheme has a faster convergence rate in solving the Navier-Stokes equations and that the rate is not sensitive to the number of mesh points or the Reynolds number. A significant acceleration of convergence is also produced for the k — t and the thermal energy equations, even though the multi-grid correction is not applied to these equations. The multi-grid method provides a stable and efficient means for local mesh refinement with only little additional computational and.memory costs. Driven cavity flows at high Reynolds numbers are computed on a number of fine meshes for both the multi-grid scheme and the local mesh-refinement scheme. Two-dimensional film cooling flow is studied using multi-grid processing and significant im-provements in the results are obtained. The non-uniformity of the flow at the slot exit and its influence on the film cooling are investigated with the fine grid resolution. A near-wall turbulence model is used. Fi lm cooling results are presented for slot injection with different mass flow ratios. ii Table of Contents Abs t r ac t i i L i s t of Tables iv L i s t of F i g u r e s V Lis t of Symbols v i Acknowledgements x 1 In t roduct ion 1 1.1 Numerical Computation of F i lm Cooling 1 1.1.1 Background 1 1.1.2 Motivation of the Present Work 3 1.2 Multi-Grid Methods 4 1.2.1 Background 4 1.2.2 Multi-Grid Method in Computational Fluid Dynamics 5 1.2.3 Local Mesh-Refinement Techniques 6 1.3 Purpose and Scope of Present Work 8 2 Mathema t i ca l M o d e l 9 2.1 Background 9 2.2 Governing Equations 10 2.3 Near-Wall Turbulence Models 11 iii 2.4 Boundary Conditions 15 3 F in i t e Difference A p p r o x i m a t i o n 18 3.1 Finite Volume Formulation 18 3.2 Ellipticity of the Discrete System 22 3.2.1 Ellipticity of the Continuous System 22 3.2.2 Ellipticity of the Discrete System 23 3.3 False Diffusion 25 4 M u l t i - G r i d Me thods 27 4.1 Basic Principle of Multi-Grid Method 27 4.1.1 Error Smoothing and Coarse Grid Correction 28 4.1.2 Full Approximation Scheme (FAS) 29 4.1.3 Full Multi-Grid Scheme 30 4.2 Important Aspects of Multi-Grid Method 32 4.2.1 Smoothing Operator 32 4.2.2 Grid Coarsening 35 4.2.3 Restriction and Prolongation 36 4.2.4 Boundary Condition Enforcement 38 4.2.5 Calculation of Thermal Energy and Turbulence 39 4.3 Performance of Multi-Grid Methods 40 4.3.1 Cost of Multi-Grid Method 41 4.3.2 Convergence of Multi-Grid Method 41 4.3.3 Computational Complexity 43 4.3.4 Closing Remarks - 45 5 L o c a l Mesh-Refinement Techniques 46 i v 5.1 Basic Principle of F A C Scheme 46 5.1.1 Construction of Composite Grid 47 5.1.2 Two-level Nonlinear F A C Scheme 49 5.1.3 Multi-grid F A C Scheme 50 5.2 Interface Treatments 52 5.2.1 Interpolation of the Slave Grid Points 53 5.2.2 Evaluation of the Residuals at Coarse Grid Interface Points . . . . 53 5.3 Self-Adaptivity of F A C 54 5.3.1 Truncation error estimate 54 5.3.2 Absolute Error Estimate 55 5.4 Performance of F A C Algorithm 56 5.4.1 Complexity of F A C Scheme 56 5.4.2 Convergence of F A C Scheme 57 6 Numerical Results 60 6.1 A Model Test - Driven cavity flow 60 6.2 Numerical Computation of F i lm Cooling Flow 70 6.2.1 Computational Grid Structure 70 6.2.2 Performance of Multi-Grid Method 70 6.2.3 Non-uniformity of the Slot Injection 76 6.2.4 F i lm Cooling Performance 80 6.2.5 Comparison of the Near-Wall Turbulence Treatments 85 7 Conclusions and Recommendations 87 Bibliography 89 v List of Tables 2.1 k - e Model Constants 12 4.1 Theoretical Smoothing Factors of S I M P L E Algorithm 43 6.1 Number of Iterations and C P U Times Required to Reduce Residual Norm by 10" 4 62 6.2 Comparison of Multi-Grid Performance with T E A C H Code 63 6.3 Asymptote Convergence Rates of F M G , F A C and T E A C H codes . . . . 63 6.4 Convergence Rate of Multi-Grid Method for F i lm Cooling Computation 75 vi List of Figures 1.1 Illustration of F i lm Cooling 2 2.1 Computational Domain of Fi lm Cooling 15 3.1 Relative Location of Staggered Grid and Calculated Variables 19 3.2 Typical Control Volume 19 4.1 A Coarse-Grid Cell Divided into Fine-Grid Cells 37 5.1 Construction of a Composite Grid 48 5.2 Flow Diagram of Multi-Grid FAC Processing 51 5.3 Composite grid interface 52 6.1 Locally Refined Gird 61 6.2 (/-velocity Profile along Vertical Centerline, Re = 100 64 6.3 [/-velocity Profile along Vertical Centerline, Re = 1,000 . . . 65 6.4 [/-velocity Profile along Vertical Centerline, Re = 10,000 65 6.5 Streamline Pattern for Cavity Flow, Re = 100 66 6.6 Streamline Pattern for Cavity Flow, Re = 1,000 67 6.7 Streamline Pattern for Cavity Flow, Re = 10,000 68 6.8 Computational Grid for Fi lm Cooling 71 6.9 Performance of Multi-Grid Method ((/—momentum Equation) 72 6.10 Performance of Multi-Grid Method (V—momentum Equation) . ' 72 6.11 Performance of Multi-Grid Method (Mass Continuity Equation) 73 vii 6.12 Performance of Multi-Grid Method (Thermal Energy Equation) 73 6.13 Performance of Multi-Grid Method (fc-equation) 74 6.14 Performance of Multi-Grid Method (e-equation) 74 6.15 Velocity Variation across the Slot ( M = 0.2) 77 6.16 Angle Variation across the Slot ( M = 0.2) 77 6.17 Flow Pattern at the Interface of Slot Injection (Af = 0.2,Refined grid(HB)) 78 6.18 Flow Pattern at the Interface of Slot Injection (M = 0.2,Coarse grid(SK)) 78 6.19 Flow Pattern at the Interface of Slot Injection (M = 0.2,Coarse grid(HB)) 79 6.20 Temperature Distribution across Slot Exit ( M = 0.2) 79 6.21 F i lm Cooling Effectiveness (Af = 0.2) 80 6.22 Longitudinal Velocity Distribution Downstream of the Slot (M = 0.2) . 81 6.23 Turbulence Intensity Distribution Downstream of the Slot ( M = 0.2) . . 82 6.24 Temperature Distribution Downstream of the Slot ( M = 0.2) 83 6.25 F i lm Cooling Effectiveness at Different Mass Ratios 84 6.26 F i lm Cooling Effectiveness for Different Near-Wall Turbulence Treatments 86 6.27 Flow Pattern Downstream of Slot Injection (Af = 0.2, Model 2) 86 viii List of Symbols aw Coefficient in finite difference equation. C Convective flux coefficient. Ci, C2, Cfj. Turbulence model constants. Upwinding and central difference for convective flux. D Diffusive flux coefficient. -Dinlet Height of main stream. E Integration constant in logarithmic law of the wall. eh Absolute error (= qh — q) on grid £lh. F Total flux coefficient, also source term of P D E . G Generation of turbulence kinetic energy. k Turbulence kinetic energy. TH 1H Interpolation operator from grid ClH to fl / , . TH 1h Restriction operator from grid QH to QH. L Differential operator. Lh Difference operator. L Fourier transform of L. L Determinant of L. Lh Fourier transform of Lh. ix lm Turbulence mixing length scale. M Mass ratio of slot flow. p Static pressure. P Pressure(= p+ \pk). Pe Cell Peclet number. Pri,Pri Laminar and Turbulent Prandtl numbers. Qh Exact solution of difference equation. q Exact solution of differential equation. qh Numerical solution of difference equation. Re Reynolds numbers (= pUL/n). ruv,rp Under-relaxation factor of momentum and continuity equations. s Dimension of slot. 5* Coefficient of linearized source term. S<f, Source term in transport equation. T Temperature. Taw Adiabatic wall temperature. Ts Coolant temperature. Too Free stream temperature. x u,v Velocity in x, and y directions. Uo,V0 Frozen velocity in x, and y directions. u'/U Turbulence intensity. ur Friction velocity (=\JTw/p). Cartesian coordinates. Difference operator (Sxcf>(x) = (<f>(x + h) — (f>(x — h))/2h) Ax, Ay Computational cell dimensions. y+ Wall coordinate (= t/C/T/i/). 7 Heat diffusity. e Turbulence dissipation rate. V Cooling effectiveness. 6 Fourier transform coordinates (^1,^2)-K von Karmann constant. Laminar and turbulence viscosity. P Smoothing factor of relaxation. V Kinematic viscosity (fi/p)-p Fluid density. o~k,o-e Constants in k — t equations. rw Wall shear stress. T Defect correction. <!> General variable. r Diffusion coefficient for <f> variable. xi Sub/superscripts fi Computational domain. C , 7, F Coarse, interface and fine grid points eff Effective value. h, H, h Fine, coarse and composite grids. n,s,e,w Cell faces. P, N, 5, E, W Grid points. s Slot injection value. w Wall value. oo Free stream value. xii Acknowledgements I would express gratitude to my supervisor, Professor M . Salcudean. Her guidance and amiable nature have made this project an enjoyable experience and her financial assis-tance has allowed to participate in this project. I would express my appreciation to Professor I.S. Gartshore and Professor U . Ascher for their helpful suggestions, comments and discussions. Thanks are due to all members of the C F D group of the Department of Mechanical Engineering, especially to Z. Abdullah and D. Sidilkovar and to P. Carter of the Institute of Applied Mathematics. Their suggestions and discussions have been very helpful. I wish also to thank Professor A.I . Livesey for proof reading my thesis. xiii Chapter 1 Introduction 1.1 Numerical Computation of Film Cooling 1.1.1 Background Fi lm cooling is a technique that has long been used to protect turbine blades from the high temperature of the surrounding gas stream. This cooling technique allows higher temperatures to be used in the gas turbine engine and hence higher efficiencies can be obtained [1,2]. In the film cooling process, cool gas is injected into the hot stream from a series of holes placed strategically over the blade surface. The hot and cool gases merge over the turbine blade, however under certain conditions, separation of the injected flow may occur resulting in the loss of cooling to the blade. The domain used for the numerical simulation includes the plenum from which injec-tion occurs, the slot, and the external region extending into the free stream. A simplified two-dimensional slot injection case is shown in Figure [1.1] and is useful mainly as a guide to elucidate the very complicated flows which occur in the actual three-dimensional case. The film cooling removes, reduces and/or redistributes thermal strains in the blades. A n important quantity allowing the evaluation of the cooling process is the film cooling effectiveness defined as J- aw •Loo / 1 i \ V= rp _rp , ^ (1-1) s J- oo where Taw is the adiabatic wall temperature, is the free stream temperature, and Ts is the cooling fluid temperature. Some of the more important parameters affecting the Chapter 1. Introduction 2 Figure 1.1: Illustration of F i lm Cooling cooling are the film cooling injection geometry, the mass flow ratio M = (pU)s/(pU)00, and Reynolds number[3]. Since 1960, film cooling research has been directed towards improved protection of components in gas turbines. Numerical simulation can provide information about fluid flow and heat transfer phenomena and, in conjunction with experimental research, can contribute to the study of the parameters which affect the film cooling process. Accurate prediction of the thermal boundary conditions surrounding the blades is important in the prediction of the thermal stress in a turbine blade. Experimental and numerical observations indicate that the flow near the injection holes is complex, with high gradients of velocity, turbulent intensity and temperature. Non-uniformities of the injected flow have been shown to have significant effects on the prediction of wall shear stress and heat transfer downstream of the injection. The objective of the present work is Chapter 1. Introduction 3 to predict film cooling effectiveness with improved accuracy and computational efficiency. 1.1.2 Motivation of the Present Work In the past decade, much progress has been made in computational fluid dynamics and heat transfer and developments in advanced algorithms have enabled analysis and nu-merical solution of highly complex flows and heat transfer. However, because of such factors as non-linearity and coupling between the equations, complexity of the flow field and large Reynolds numbers, difficulties in direct applications of present numerical codes, e.g., T E A C H , etc., to film cooling computation still exist. The use of a coarse grid can yield a converged solution but such a solution may often be inaccurate in practice. In-accuracies may result in the appearance of false diffusion, and inaccurate prediction of heat transfer near solid walls. Refining the computational grid, although usually leading to improved representation, is accompanied by deterioration in the convergence rate of the solution and unreasonable use of computer memory. Previous studies [4] have shown that, in order to obtain accurate numerical prediction of the film cooling process, resolution near the slot and turbine wall is needed because the flow in these regions has considerable influence on the heat transfer. It has been shown that finer resolution in some regions of the domain provides accurate prediction of the local features (such as non-uniformity or flow separation), and at the same time increases the overall accuracy. These observations present the motivation for the development of a technique which can provide mesh refinement as well as acceleration of the convergence rate of the solution. Chapter 1. Introduction 4 1.2 Multi-Grid Methods 1.2.1 Background Proper approximation to a physical problem and adequate numerical resolution require a grid which is fine compared to the phenomenon of interest. The number of discrete variables and algebraic relations between the variables is almost always very large. Since the partial differential equations indicate the local relation of physical conservation laws, the resulting discretization system will also be local and each algebraic relation will involve only a small number of neighboring variables. Usually, the iteration method for driving such a discrete system toward a solution on some given fine grid is a local process, and is inefficient because it affects the non-local features of the solution very slowly. The need for repeating the iteration sweeps therefore increases when the grid is refined. The multi-grid method efficiently eliminates most of the work related to the repe-tition of iterations. In multi-grid processing, discretization and solution processes are intermixed with a sequence of grids of geometrically decreasing mesh sizes. The coop-erative solution processes on these grids involve relaxation sweeps over each of them, coarse-grid-to-fine-grid interpolation of corrections, and fine-to-coarse transfer of residu-als. Typically, each grid supplies corrections to the equations of the coarser grid and to the approximate solution of the next finer grid. The combination of the relaxations on each level produces local and global information for efficient solution of the equations. The overall work required for the multi-grid process is only a few 'work units', where a work unit is defined as the amount of computational work needed for the solution of the algebraic equations on the finest grid. Multi-grid methods rapidly solve the algebraic system of equations and the solution converges at a rate not sensitive to the number of grid points. Chapter 1. Introduction 5 1 . 2 . 2 Multi-Grid Method in Computational Fluid Dynamics The potential of the multi-grid method in the field of computational fluid dynamics has been adequately exposed and demonstrated by many researchers in the literature since the later 1970's. However, application of this method to the solution of a coupled nonlinear system of Navier-Stokes equations still poses several questions which have been studied by various investigators. As a pioneer in this area, Brandt [5] introduced a Multi-Level Adaptive Technique (MLAT) with a nonlinear Fast Approximation Scheme (FAS) and subsequently devel-oped a Distributive Gauss-Seidel (DGS) relaxation method as a smoother for solving the Navier-Stokes equation on a staggered grid [6]. Fuchs [7] examined the DGS scheme as applied to recirculating flows with Dirichlet and parabolic boundary conditions, and observed that the multi-grid method has a fast rate of convergence with both types of boundary conditions and is not sensitive to the number of mesh points. Later he pro-posed a new Pressure Gradient Averaging (PGA) scheme [8] for updating the pressure gradients in the DGS scheme and found that the new P G A scheme has an even better convergence rate for those high Reynolds numbers for which the DGS scheme is unstable. Since the choice of the relaxation smoother for the Navier-Stokes equations is problem-dependent and since the treatment of the coupling between pressure and velocity fields is crucial to rapid convergence of the algorithm, researchers have used different approaches. Ghia et al [9] have used a Coupled Strongly Implicit Procedure to simultaneously relax the stream-function and vorticity equation. Their fine-mesh solutions for cavity flow at high Reynolds number, using the multi-grid method, provided flow structures which compared well with the results of other researchers and experiments. Vanka [10] devel-oped the Symmetrical Coupled Gauss-Seidel (SCGS) relaxation method for the primitive variable solution. In the SCGS scheme, four velocities and one pressure corresponding i I Chapter 1. Introduction 6 to one finite-difference node are simultaneously updated by the inversion of a small ma-trix equation. The simultaneous relaxation showed a good smoothing rate for strongly coupled equations. Another widely used algorithm in a number of diverse engineering situations is the S I M P L E algorithm of Patankar and Spalding [11] in which the variables are located on a staggered grid and a decoupled solution strategy is used for solving the non-linearly coupled Navier-Stokes equations. Recently, Sivaloganathan and Shaw [12], in studying the smoothing properties of the scheme, introduced a multi-grid approach with a continuity-preserving procedure for grid coarsening. However, in their calculation of the cavity flow, the grid seemed insufficiently fine and the efficiency of convergence of the multi-grid method for very fine grids was not clear. Since the choices of the multi-grid components are problem dependent and very im-portant to the efficiency of the scheme itself, studies for the Navier-Stokes equations are ongoing, especially for practical problems. The application to turbulent flow and heat transfer is still a new topic. Special difficulties arise related to the modelling of turbulence. 1.2.3 Local Mesh-Refinement Techniques If computational power were free and unlimited, the differential system on a global uni-form grid could be resolved with the smallest desired mesh size. But most applications require concession to coarser global scales with the associated loss of important physics and its impact. The development of a variety of mesh-refinement techniques is motivated by need for local resolution of special phenomena. Numerical discretization and solution processes must take account of local resolution without a dramatic increase in the over-all computational cost. Classical local mesh-refinement methods are developed to avoid difficulties in the differencing structure and the solution algorithms of the conventional approaches (such as coordinate transformation and finite element method) [13]. These Chapter 1. Introduction 7 methods use global (coarse) grid approximation as boundary values for a solution of local (fine) grid. Since the communication between levels is only through the emerging solu-tion, one implication is that either the approximation must accept coarse-grid accuracy at grid interfaces (thereby limiting fine-grid accuracy) or some extra cycling between these levels is necessary. Increased resolution usually requires increased interlevel cycling. With the development of multi-grid methods, multi-grid local mesh-refinement tech-niques appeared as good strategies, different mesh sizes being allocated to different re-gions of the domain. The multi-grid method provides a means of local mesh refinement so that local resolution of the solution is increased. In his M L A T scheme, Brandt [5,14] presented a technique for local mesh refinement similar to the conventional multi-grid algorithm but with the additional characteristic that the fine grids need not be global. Local mesh-refinement, by extending uniform grids, is relatively easy and inexpensively implemented. With the modification of an exchange-rate algorithm (A-FMG) and with the optional-mesh criteria, the method shows that singular problems are solved as effi-ciently as regular problems. In 1985, Fuchs [15] reported the application of the local mesh refinement technique to some separated flows and examined the stability and efficiency of the method. After Brandt's work, many different types of multi-grid local mesh-refinement meth-ods, such as the Local Defect Correction (LDC) method and the Domain Decomposition method [16] appeared in the literature [19,18]. In 1984, McCormick [20] introduced the Fast Adaptive Composite Grid (FAC) method which has the distinctive feature that the communication on the internal boundary between the coarse region and the refined re-gion is through the equations themselves rather than the solution. F A C uses global and local uniform grids to define the composite grid problem. Both the discrete system and the iterations on the composite grid are conservative. In the F A C iteration processing, the coarse grid is used not only for fast solution on the fine grid but also for transmitting Chapter 1. Introduction 8 fine grid accuracy throughout the domain. 1.3 Purpose and Scope of Present Work This work establishes a multi-grid algorithm with local mesh-refinement and investigates the applicability of the algorithm for the solution of the transport equations in turbulent flows associated with heat transfer. A nonlinear multi-grid algorithm is applied to a finite difference system on a staggered grid and the S I M P L E scheme is used as a relaxation smoother for solving the momentum and continuity equation. The local mesh-refinement is done through the application of the F A C method. As a model test for the algorithm, driven cavity flows at high Reynolds number on very fine meshes are calculated. In the application to film cooling, a two-dimensional geometry is chosen in which the slot orifice is normal to the main flow. The flow at the slot exit and its influence on the effectiveness of film cooling is studied by refining the computational grid near the slot and solid wall to eliminate the false diffusion and also to provide better flow information near the wall for precise prediction of the heat transfer. The multi-grid method is used to speed the convergence of the iterations in solving the difference system of momentum and continuity equations and to provide more accurate flow fields for calculation of the turbulence and thermal energy equations. A new treatment of the near wall turbulence computation is introduced. The influence of the mass ratio of the slot flow on film cooling effectiveness is studied. Because of the lack of the experimental data, the numerical results are compared with results from other C F D codes. Chapter 2 Mathematical Model This chapter presents the background for the numerical solution of the differential equa-tions for turbulent flow and heat transfer. The general governing equations and the turbulence k — e model are described. The standard turbulence wall function is imple-mented by combining the algebraic turbulence model with the k — e model. Finally, the boundary conditions for the two-dimensional film cooling flow in the vertical slot injection case are described. 2.1 Background The fundamental differential equations of fluid dynamics and heat transfer represent the laws of conservation of mass, momentum and energy. Despite all the recent advances in computer technology, direct numerical solution of these equations is not possible for practical flows because of the complex nature of turbulence. For this reason, the statis-tically regular nature of turbulence flow is studied through the time-averaging procedure first proposed by Reynolds over a hundred years ago. However, because the equations are non-linear, the averaging procedure produces extra unknown terms (such as the turbu-lent or Reynolds stresses) representing the transport of the mean momentum and heat by the turbulent motion, and the equations no longer constitute a closed system. Turbulent models simulate the transport terms in the mean flow equations and close the system so that the effect of turbulence on the mean flow behavior can be simulated. 9 Chapter 2. Mathematical Model 10 2.2 Governing Equations * For a two-dimensional, incompressible, and statistically steady flow, the flow field and temperature distribution are governed by the continuity equation, the Reynolds-averaged Navier-Stokes equations with Boussinesq eddy-viscosity formulation and thermal energy equations: (for details see [22]) continuity equation: | ( ^ > + &V) - 0, (2.1) x-momentum equation: d . TTTT, < 9 / r r T / N dP d ,n dU, dr ,dU dV Yx{PUU) + -(PUV) = - - + r x ( 2 M - ) + + - ) ] , (2.2) y-momentum equation: d , rrT,N d , T,T„ dP d t n dV, dr ,dU dV„ Yx{PUV) + tfpVV) = + ^ e / f - ) + _ [ * , , ( _ + - ) ] , (2.3) thermal energy equation: f^UT) + fapVT) - + » (2.4) with P = P+\pk (2.5) where x and y are the Cartesian coordinates; U and V are the velocity components in the x and y directions, respectively; /), p, and T are the fluid density, pressure, temperature. The effective viscosity and heat diffusivity are obtained from PeSS = Pt + Pi (2-6) Chapter 2. Mathematical Model 11 and where P r ; is the laminar Prandtl number; Prt is the turbulent Prandtl number; Hi is the laminar viscosity, and fit is the turbulent eddy viscosity obtained from the value of the turbulent kinetic energy k and its dissipation rate e according to: k2 fit = pC^— (2.8) The distributions of k and c are obtained with the Launder and Spalding [23] high Reynolds number k — t turbulence model which can be written as: fc-transport equation: UfUk) + UpVk) = U^ld±) + f ( ^ ) + G- pe, (2.9) ox oy ox (7k ox oy Ck oy e-transport equation: >«>+" *<^£> + £C£> + " <4 <"°> where the production rate of the turbulence kinetic energy can be defined as ^ , n r,3{/., ,dV.j. ,dU dV.7, . * * \ G = "'<2l<a7> + <8»>1+ <& + *>> ' ' The standard Launder and Spalding constants in equations (2.9-2.10) are listed in Table 2.1. 2.3 Near-Wall Turbulence Models The original k — e turbulence model neglects viscous effects which are important in the vicinity of solid boundaries. In this region, the local turbulence Reynolds number is Chapter 2. Mathematical Model 12 Table 2.1: k - t Model Constants C2 0.09 1.44 1.92 1.0 1.3 low enough for molecular viscosity to affect the processes of generation, destruction and transport of turbulence. Although the thickness of this viscosity-affected zone is usually much smaller than the overall width of the flow, its effects can extend over the whole flow [24] with the result that in the numerical analysis where an accurate prediction of near-wall turbulence is required, e.g., heat transfer analysis of high Prandtl number fluids, some errors may appear in the predictions. The majority of existing techniques for the treatment of the near wall region are focused on the wall function formulation. The standard wall function formulation of Launder and Spalding [23] assumes that the flow in the wall region behaves locally as one-dimensional Couette flow and conditions similar to those of an equilibrium turbulent boundary layer prevail. The computational nodes immediately adjacent to the solid wall, are located in the fully turbulent region sufficiently far away from the wall for viscous effects to be negligible, i.e., at y + >30. This treatment, which makes no attempt to resolve the flow within the laminar viscous sublayer, requires fewer mesh cells and so it is less likely to be affected by the computational stiffness associated with highly stretched grids. However, heat transfer effects may not be modelled accurately by a wall function. To allow for local variations of the turbulent kinetic energy k and its dissipation rate e to evaluate the generation and destruction of k near the wall, Chieng and Launder [25] developed a two-layer near-wall model which Amano [24] later improved to a three-layer model. However, the performance of these multi-layer models is not always satisfactory and the problem concerning the near-wall turbulence modelling is Chapter 2. Mathematical Model 13 still open. One important factor is that the coarse grid limits the numerical accuracy in this region. In a number of flow situations, e.g. when adverse pressure gradients are strong, the law of the wall is invalid and the conservation equations with low Reynolds number terms may have to be solved to the wall. The k — t model can be modified by the addition of model functions in the k, e equations and turbulent viscosity formulations. These model functions are meant to meet the limiting behavior of the near wall turbulence [28]. Low Reynolds number turbulence models are established by new formulations of the turbu-lence equations in the near wall region [29]. However in various practical problems, there is little experimental data to guide the development, and the reliability and applicability of the models are difficult to estimate [30]. The present work uses the standard turbulence wall function for coarse grid solution and, when the grid is refined, uses the algebraic viscosity formulation to model the turbulent viscosity fit for computational nodes close to the viscous-sublayer. The reasons are that the appropriate modelling of fit from the viscous sublayer to the inertial layer is crucial for accurate prediction of the velocity profiles, and normal wall functions do not perform well in this region. If a grid point p next to the wall is within the inertial layer, i.e., 30 < < 150, then the standard wall function is applied. The computational procedure for this cell is modified under local equilibrium conditions in this layer. One important modification is for the wall shear stress which evaluates the near wall diffusion terms in the momentum equations (2.2,2.3) and in the generation term of the k equation (2.9): Another modification is for e at the wall-adjacent node which is evaluated from local Chapter 2. Mathematical Model 14 equilibrium conditions and is used for the boundary condition of the e equation e = C3J4-^- (2.13) and for the mean dissipation rate used in evaluation of the dissipation term — pe in the k-equation (2.9): £3/4^3/2 e = - * ln(JS?y+). (2.14) The turbulent viscosity at point p is calculated by the turbulent eddy-viscosity formula-tion (2.8). If point p is within y+ < 30 , an algebraic model for the turbulent viscosity based on a Van Driest approach to the inner region with a damped law of the wall is applied. The formula is: where the mixing length scale 7 m = « y ( l - e x p [ - ^ ] ) (2.16) with K = 0.41 and A* = 25. The following alternatives for the modification of wall shear stress are considered. 1. Calculate the wall shear stress on the point q which is away from the wall, i.e., 30 < y* < 150, assuming that the Reynolds shear stress is nearly constant in the near-wall region. 2. Calculate the wall shear stress on the point q with the eddy-viscosity concept, i.e., rw = itt— (2.17) dy Chapter 2. Mathematical Model 15 H Umain s o LL C "co 2 B Symmetric Boundary I Adiatic Boundary / / / / / / 10s 30s 20s 3s « •o c 3 o m 3 o t E Uslot Secondary Flow Inlet Figure 2.1: Computational Domain of Fi lm Cooling Comparison between these two alternatives do not show significant difference. The tur-bulence energy dissipation e on the boundary and the mean dissipation rate for the k equation are expressed as follows: e = C2J4^-. C?J4k3'2 (2.18) (2.19) This treatment is not accurate; however no better choice is available presently. 2.4 Boundary Conditions A typical computational domain for the film cooling flow is shown in Figure (2.1). Four types of boundaries and the corresponding boundary conditions are: Chapter 2. Mathematical Model 16 1. Main inflow (AB): values of all the variables, except pressure, are specified. For the present flow configuration, uniform profiles for U, T, k, and e are imposed; V is set to zero. The turbulent kinetic energy, fcoo, is determined by: ^ = 0.00015-C& (2.20) and it also can be determined in the later work from the experimentally measured boundary layer turbulence intensity; c has to be estimated from k and a character-istic length scale: too = k%2/l, (2.21) where / is assumed as the length of entrance Diniet-2. Slot inflow (DE): similar to the main inflow, uniform V velocity is imposed, U = 0; the inlet turbulent kinetic energy is imposed as *s = 1.5(^)2t/ 8 2 (2.22) and the inlet dissipation rate of energy is imposed as e8 = (2.23) s where / is the turbulence length scale that is assumed to be equal to s. The value of turbulence intensity (jj)a is specified as 0.05. 3. Outflow (GH): the outflow boundary can be located sufficiently far downstream from the merging and the recirculating flow region to ensure that the flow in the upstream region is not affected by downstream conditions. Zero stream-wise gra-dient condition | ^ = 0 is then applied to all variables at this location. In addition, since the S I M P L E algorithm needs mass conservation in the computational domain, the mass correction is imposed for the U velocity at this point in order to help the convergence. Chapter 2. Mathematical Model 17 4. Symmetry Axis (AH): the normal velocity, V , is set to zero along the axis of symmetry and a zero cross-stream gradient condition, | ^ = 0, is imposed for all variables. This condition is different from that of the experimental case because the experiments are carried out in a wind tunnel with a wall. No influence between them is assumed if boundary condition is set far away from the solid wall. Numerical experiments have shown that this assumption is valid. 5. Solid Walls (BI, ID, J E and JG): zero mass and heat fluxes are imposed at all solid boundaries to satisfy the no-slip and adiabatic boundary condition. As described in the previous section, a combination of the standard turbulence wall function of Launder and Spalding [23] and the algebraic viscosity formulation of Cecibi and Smith [30] is applied for determination of turbulence variables. The standard wall function is used at the walls of the coolant orifice. Chapter 3 Finite Difference Approximation This chapter describes the finite volume formulation for the general transport equations. For the multi-grid approach, it establishes the ellipticity of the discrete system of the continuity and the Navier-Stokes equations. This chapter also discusses the false diffusion resulting from the difference approximation. 3.1 Finite Volume Formulation The discretized system of the governing equations is formed with a staggered grid with variables located as shown in Figure 3.1. Due to the staggering of the mesh, three different types of control volumes are required for momentum equations, continuity equation and scalars (such as T, k and e) in the interior. Modifications of these control volumes near the boundaries are straightforward. For numerical purposes, the governing equations described in the previous chapter can be represented by the general transport equation ^ W ) + | W ) = | [ r ( | ) ] + | l r ( | ) ] + 5 „ ( 3 , ) where <f> can be replaced for different equations, e.g., <f> = U for U—momentum equation; T is a general diffusivity coefficient, and a general source term. The finite volume form of equation (3.1) is obtained by integrating the respective control volumes (Figure 3.2 shows a typical control volume for any variables), i.e., JUJSW*) + |W) - £ir<!)J - |[r(|)] - (3.2) 18 Chapter 3. Finite Difference Approximation $ I * 1 ! Figure 3.1: Relative Location of Staggered Grid and Calculated Variables w W N x Figure 3.2: Typical Control Volume Chapter 3. Finite Difference Approximation 20 and with Gauss' divergence theorem, the volume integrals can be transformed into surface integrals: £(PU<f> - T^)Xcdy - f\PU4> - r^)Xwdy (3.3) £> v* - - />v< - (3-4) dy,Vn * JXwK" r dy -JJ S^dxdy = 0, (3.5) cv where Fe,...,Fs represent the sum of convective and diffusive fluxes across the faces e,. . . , 5, e.g., Fe = Ce + De with Ce = r(pU<t>)Xcdy (3.6) Jys and B-/>ri>-* <3J> The convective and diffusive flux terms (Ce,De, etc.) are approximated by the hybrid differencing scheme which is the combination of two types of finite difference approxima-tions: central difference and upwind difference. The central differencing approximation assumes, e.g., that <f> varies linearly between the grid points E and p (see Figure [3.2]), i.e., C e = PeUe{^±^)Ay (3.8) and D. = r . ( * = ^ i £ ) A , . (3.9) LAX and is used at low cell Peclet number, i.e., Pt = \pUA/V\ is smaller than 2. The upwind differencing applied at higher Peclet numbers assumes that the upwind value of <f> at face Chapter 3. Finite Difference Approximation 21 e prevail: Ce = peUe<t>PAy for Ue > 0 (3.10) Ce = peUe(f>EAy for Ue < 0 (3.11) The source term S<j, is linearized as follows : JJ S+dxdy = SUP + S£t, (3.12) where Sp and Sfj are derived using central differencing approximations. Substituting the source and flux terms into equation (3.5) yields the general finite difference equation: (aP - Sp)<t>p = aN<pN + a5</>5 + aE4>E + a\v<f>w + (3.13) with op = a,- (i = N,s,E,w) i and aN = max(|C*n/2|,Z)n) - Cn/2 (3.14) as = max(|C,/2|, D.) + Cs/2 (3.15) aE = max(|C e /2 | , De) - Ce/2 (3.16) aw = maxdC^/21, Dw) + Cw/2 (3.17) The modification of the discretization near the boundaries is clear and straightfor-ward. The algebraic equations (3.13) are solved by an iterative procedure, which will be described in next chapter. Chapter 3. Finite Difference Approximation 22 3.2 Ellipticity of the Discrete System The smoothing properties of an iterative method depend upon the ellipticity of the dis-crete representation of the partial differential equations of the problem. The system of continuity and momentum equations (2.1) (2.2) (2.3) can be written as operator form: ft dxdy d_ dx u 0 u 9 2 ' dxdy d_ dy V — 0 a_ dx d dy 0 p 0 where (3.18) (3.19) 3.2.1 Ellipticity of the Continuous System The non-linear system (3.18) is linearized by freezing p and \i at p0 and fi0 respectively resulting from the previous iteration. Similarly, the contribution of the velocities U and V to the non-linear term is at UQ and Vo respectively resulting from the previous iteration. The linearized system may then be written as LQ = F, (3.20) where L = Wdxdy -fJ>o dxdy d_ dx d_ dy d_ dx C0 ~ ^ 0 ( 2 g j 2 + QyS) — Co- Vo(2§y2 + §12) §1 0 (3.21) Chapter 3. Finite Difference Approximation 23 with the linearized convective operator Co = po(Uo-§^ + ^ 0 ^ ) 5 Q = (U,V,P)T, and F = (0,0,0) T . The Fourier transform L of L is given by Co + p0{2Q\ + 0j) nQ0x02 \0X L= po0x02 Co + iio{20l + 0*) i92 , (3-22) i6>! i0 2 0 where C 0 = P o K ^ c A + K A ) is the Fourier transform of C 0 from into (0\,02)-Investigation of the ellipticity of the system considers the determinant L of L: L = detl = (01 + 0l)[po(O2 + 01) + PoWi + VQ02)}. (3.23) The system (3.20) is defined to be elliptic if L is non zero for all real 0 = (0X, 02) not equal to (0,0) [31]. Clearly this is the case for all non-zero fi0 so that the linearized problem is elliptic and the non-linear problem (3.18) is also considered to be elliptic. For the high Reynolds number, ellipticity is almost lost in the 02) plane perpen-dicular to a local streamline. 3.2.2 Ellipticity of the Discrete System The resulting discrete equations on a grid £lh C D of mesh size h are: LhQh = Fn^ (3.24) where (3.25) < -p6hJhy 6* -fi6hx6hy , ahv 6hy % 0 Qh = (uh,vh,ph), and uh, vh, ph, which approximate U, V , P, are defined on tth. The component operators of Lh are defined by: a^<f>p = aP(j>p — aE<t>E — aw<t>w ~ UN4>N — as<t>s (3.26) Chapter 3. Finite Difference Approximation 24 and Shx(f)p = <i>e - <j>M h (3.27) Similar expressions can be written for a£, 8y. Definitions of the other elements can be found in [11]. According to the definition of Brandt [31], the discrete Fourier transform Lh(9) of Lh from (x,y) into (9x,92) is: Lh{9) = a^O 4 / i 0 5 i 5 2 / / i 2 2\Si/h Afi0sis2/h2 a^9 2\s2/h 213x1 h 2\s2/h 0 (3.28) w here ahu(9) = M m a x ( ^ , l\p0U0\) + X m a X < X ' ^ ° y o l ) - i ^ ( f / 0 s i n 5 1 + V 0 s in5 2 ) , n (3.29) °HM = xmax(?F' \LPOVOL) + !T m a x ( ?f' ^ ot/ol) - i ^ ( J 7 0 s i n ^ + V 0 s in« 2 ) , a (3.30) S l = sin(5j/2), 5 2 = sin(0 2/2). (3.31) The determinant Z' l(5) of the discrete Fourier transform Lh(0) from (x, y) into (9x,92) is: J>(0) = 4^(0) + sl~ahv(e) - %w\s\lh2\lh\ (3.32) Thus Im(/>(0)) = -4p 0(s? + ^ ( t / o s i n A + VJ, sin 0 2 ) / / i 3 , (3.33) Chapter 3. Finite Difference Approximation 25 which is zero along lines V0/U0 = — s i n ^ / sin#2. Since Re(~ahu(9)) > 4^ 0(25 2 + s22)/h2, Mrt(9)) > ^ o(2s22 + s2)/h2, (3.34) (3.35) it follows that Re(Lh(9)) > ^ 0[4(.2 + sl)/h2)\ (3.36) Clearly \Lh(9)\ is small in the region of (91,92) = (0,0) or in cases of larger Reynolds numbers, therefore L has a good ft-ellipticity measure and the discretization is discretely elliptic in the sense of Brandt. 3.3 False Diffusion In the construction of the convective fluxes Ce,Cw,Ca, and Cn (see Figure [3.2]), the central difference scheme, for example, is used for small Peclet number flow and is second order accurate. However for large Peclect numbers, the central difference scheme is numerically unstable and the upwinding difference scheme is applied instead {pU)w{(j)w + <f>p) (3.37) (3.38) (pU)w<f>w for Uw > 0 (3.39) (pV)Js for V.>0 (3.40) Chapter 3. Finite Difference Approximation 26 The upwind difference scheme can be interpreted as a correction of the central difference scheme by adding an artificial diffusive term: -l-[&x(PU<j>x)x + Ay(pV<t>y)y}, (3.41) i.e., CZ = CZ - \{pU)w{<t>P - <j>w) = (pU)w<i>w for Uw > 0 (3.42) CV. = C°s ~ \(pV).(<f>p - <£s) = (pVUs for Vs> 0 (3.43) This artificial term (3.41), called "false diffusion" has a form similar to a diffusion term with the effective diffusion coefficients T x = l-AxPU and Ty = ^AypV. (3.44) The coefficients show that false diffusion is a multi-dimensional phenomenon, which oc-curs when the flow is oblique to the grid lines or when there is a non-zero gradient of the dependent variable in the direction normal to the flow. The false diffusion artificially increases the physical diffusion coefficients; and results in the smearing of the gradients in the flow field. The artificial diffusion can become comparable to, or even larger than physical diffusion, and can lead to significant errors. It is noteworthy that skewness and strong gradients are often prevailing conditions in recirculating flows. One basic approach that reduces the false diffusion to an acceptable level uses small mesh A x and Ay in critical regions of the flow. With the multi-grid technique, reduc-tion of the false diffusion becomes possible in practical problems. Other approaches to eliminate false diffusion include the alignment of the grid with the flow and higher order difference schemes [32,33]. Chapter 4 Multi-Grid Methods This chapter describes the principle of multi-grid method. It presents the important com-ponents in the design of the present multi-grid scheme, such as the relaxation smoother, interpolation and restriction, grid coarsening, boundary conditions, and the decoupled solution of the turbulence and energy equations. This chapter also evaluates the perfor-mance of the multi-grid method. 4.1 Basic Principle of Multi-Grid Method A boundary-value problem for the differential equation of a mathematical model Lq = F (4.1) is discretized on a given grid O'1 by a difference scheme LhQh = ph (4 _2) with certain boundary conditions. These difference equations are then solved by a cer-tain iteration performed in a certain ordering. The truncation error of the difference approximation is Th = Lq- LhQh (4.3) and the residual error after several iteration sweeps is rh = Fh- Lhqh, (4.4) 27 Chapter 4. Multi-Grid Methods 28 where q is the current solution approximation of Qh. The accuracy of the numerical solution is generally determined by these two errors. Computation of the numerical solution within a prescribed accuracy requires a dif-ference scheme, a computational grid and an algebraic method for solving the difference system so that both the truncation error and the algebraic error can be within the accu-racy limit. Higher order difference schemes and refined mesh are essential for the elimination of truncation errors of the difference system. A n iteration method is usually chosen as the method for solving the difference system. Iterating sweeps reduce the residual error and drive the difference system toward a solution. The multi-grid method is a numerical strategy for enhancing the convergence rate of the iteration procedure. 4.1.1 Error Smoothing and Coarse Grid Correction Numerical experiments show that the convergence of relaxation is fast at the beginning but very slow after a few sweeps because the relaxation appears to reduce efficiently the non-smooth error components whose wave-lengths are smaller than or comparable to the grid size. When the error is smooth, the convergence is slow. The problem in the relaxation is the elimination of low frequency error components. However, a wave-length which is long relative to a fine mesh is short relative to a coarser mesh and the smoothing error can be significantly reduced on a coarser grid. The multi-grid method is based on this observation. The method employs a hierarchy of grids and, as the convergence of the solution on the fine grid becomes slow, it switches to a coarser grid where the longer-wavelength errors are more effectively damped. The solution on the finer grid then corrected appropriately reflects the removal of these long-wavelength components of the error. For example, if the equation(4.2) is assumed to be linear, after several relaxation Chapter 4. Multi-Grid Methods 29 sweeps, the residuals of the fine grid equations are given by (4.4). The error eh = Qh — qh will satisfy Lhvh = rh. (4.5) Since eh is smooth, it can be approximated by a coarse grid function vH satisfying the equation LHvH = rH, (4.6) where LH is a grid H difference approximation to the differential operator and r " = / * V \ (4.7) where If? is a fine-to-coarse restriction operator. After a certain approximate solution vH of (4.5) is obtained, it can be used to correct the fine grid solution in the following way qh*-qh + IhHv\ (4.8) where Ijj is a coarse-to-fine prolongation operator. The process of calculating rh, transferring it to the coarse grid, solving the coarse grid equations for vH and interpolating the result and adding it to the fine grid approximation is called "coarse grid correction". 4.1.2 Full Approximation Scheme (FAS) The construction of coarse grid equations for the error eh in the fine grid approximation qh to the solution of (4.2) is called the Correction Scheme (CS). In nonlinear problems, and in the local refinement method, each of the variables should be correctly represented in all the terms for coarse grids. The Full Approximation Scheme (FAS) Brandt[5] uses Chapter 4. Multi-Grid Methods 30 coarse grid equations in terms of a different function. Instead of vH in equation (4.6), the exact variable q11 qH = I^qh + vH (4.9) is used. The coarse grid function which approximates the full solution on the coarse grid satisfies LHQH = pH^ ( 4 1 Q ) where FH = I»rh-rLH(I?qh). (4.11) After this equation is approximately solved, qH is used to correct the fine grid approxi-mation in the following way q h - q h + IH(qH-IhHqh). (4.12) A n estimate of the truncation error is an important by-product of the FAS scheme and is useful in defining natural stopping criteria and grid adaptation criteria. It can be shown that FH - I»Fh = L"(4V) - I?Lhqh, (4.13) which is exactly the Vth approximations to the £lH truncation errors. 4.1.3 Full Multi-Grid Scheme If a sequence of grids h\ > h2> ... > where h,M = h and hk = 2 / i* + 1 is assumed, the k grid equation for the FAS scheme is written as LkQk = Fk^ (4.14) Chapter 4. Multi-Grid Methods 31 where Fk = l£+1(Fk+1 - Lk+1qk+1) + Lk(I^+1qk+1) for k < M (4.15) and pM _ ph. (4.16) The basic multi-grid scheme starts with relaxation on the finest grid, and proceeds to coarser grids to reduce smooth error components. However the coarse grid solution is a good initial approximation for the fine grid solution. Therefore the solution should be first started on a coarse grid. The solution should proceed to finer grids only when suffi-ciently good approximations to the solution of the original equation have been achieved. This procedure yields the so-called Full Multi-Grid (FMG) cycling scheme. The basic procedure of the two-grid (k, k + 1) F M G scheme is: • Step 1. Coarse-Grid Pre-Relaxation: Perform Nx relaxation sweeps on LkQk = Fk starting with initial guess uk and obtain qk. • Step 2. Prolongation: Interpolate qk into fifc+1 and obtain 9*+^: q k ^ t e = Ik+1 qk. • Step 3. Fine-Grid Relaxation: Perform N2 relaxation sweeps on Lk+1Qh+1 = Fk+1 with initial guess q^1 and obtain qk+1. • Step 4. Restriction: Restrict the residual r f c + 1 = Fk+1 - Lk+1qk+1 into tlk and obtain rk: rk = rk+1. • Step 5. Coarse-Grid Correction: Compute the correction vk from LkQk = rk starting with initial guess vk = 0. Chapter 4. Multi-Grid Methods 32 • Step 6. Prolongation: Interpolate vk into Slk+1 and obtain q ^ e c t = vk. • Step 7. Fine-Grid Post-Relaxation: Perform N3 relaxation sweeps on Q,k+1 with o*<iLct a n d obtain qk*w-The basic aim of F M G scheme is to guarantee a good initial guess for the finest grid fiM before any processing is done. It follows that F M G scheme will often be used without a starting guess or, more precisely, with uM = 0 only. 4.2 Important Aspects of M u l t i - G r i d M e t h o d The development of a good multi-grid code needs prior knowledge of the obtainable efficiency followed by implementation. The first step is the construction of a relaxation scheme with a high smoothing rate. This step has to be followed by careful design of grid coarsening, inter-grid transfers within the interior and with the boundary. 4.2.1 Smoothing Operator The choice of an efficient relaxation scheme is problem dependent and an important aspect in the multi-grid technique. In the past, a number of smoothers have been used on a variety of problems. A S I M P L E pressure-correction scheme in which the momentum and continuity equations are relaxed in a coupled manner is used as the smoother in this work. The reason of this choice is not only its wide use in engineering applications but also its good smoothing behavior. In S IMPLE, an approximate pressure field is first used to solve the momentum equa-tions. Second, a pressure-correction equation is derived by the combination of the finite-differenced continuity equation and approximate forms of the momentum equations. The Chapter 4. Multi-Grid Methods 33 implementation of the SIMPLE algorithm in the multi-grid procedure is presented in the following. The system (3.24) can be written as LhQh = Fh^ (4>17) where the source terms which arise during the multi-grid process are included in Fh. If qh = (uh,vh,ph) is the approximate solution to (4.17) and qh = (uh,vh,ph) = qh + Sq = qh + (Su,Sv,Sp) is the approximation obtained from qh after a pressure-correction iteration. Then substitute qh into equation (4.17): (ahuuh - fi6hjy + 6hxPh - t) + (ahJu - ^hx8hy8v + 6hJp) = 0 (4.18) (ahvvh - p,6hx6y + 8hxph - fhv) + (ahv8v - p.6hx6hy8u + 6hy6P) = 0 (4.19) (6hxuk + 6hyvh - ft) + (ShxSu + 6hySv) = 0, (4.20) assuming that the current approximation qh already satisfies the momentum equation, the correction Sq satisfies ahu8u - »8hx6hy8v + 6hx8p = 0 (4.21) ahvSv - nShxShySu + ShySp = 0 (4.22) 8hxSu + ShySv = -Shxuh + Shyvh - / ; (4.23) The S I M P L E pressure-correction method neglects the mixed derivative terms and diagonalizes the operator a£ and aj to give dhuSu = -ShxSp (4.24) dhvSv = -ShySp (4.25) Chapter 4. Multi-Grid Methods 34 (ShJu + Shy6v) = -(6hxuh + 6hyvh - /J), (4.26) where dhu = ahu and dhv = <z£ . Hence the corrections 6q satisfy 6u = -(d*)-*6*6p (4.27) 6v = -{dhv)-l8hy8p (4.28) [M)- 1 ** + lhy{dhv)-'8hy))8p = 8hxuh + 8hyvh - fhp. (4.29) The algorithm proceeds by applying a few Line-by-Line Gauss-Seidel sweeps to the 'Poisson' equation (4.29) for 6p and hence obtain 6u, 8v from equations (4.27) and (4.28). Under-relaxation, which has the effect of a false time dependent term in the equations and has stabilizing characteristics, is implemented by updating the solution qh by qh: uh = u h + ruvSu, (4.30) vh = vh + ruv6v, (4.31) ph = ph + rp6p, (4.32) where ruv,rp are under-relaxation parameters. This procedure, applied after a few relaxation sweeps of the momentum equations uses a relaxation parameter r m o n and aims to satisfy the continuity equation assuming that the momentum equations are satisfied. The pressure-correction algorithm used as a smoother employs an under-relaxation parameter r m o m of approximately 0.5 in relaxing the momentum equations. Empirically this has been found to be most effective at low Reynolds numbers. A gradual reduction of ?"mom a ^ higher Reynolds numbers tends to improve the smoothing properties of the algorithm. At high Reynolds number, convergence can be accelerated by varying the Chapter 4. Multi-Grid Methods 35 relaxation parameter on the grids in relation to mesh Reynolds number (which is related to the amount of artificial viscosity added on each grid). As yet, little theoretical analysis has been carried out on optimizing relaxation parameters. The system arising from equation (4.29) is easily shown to be singular; and a unique solution for the pressure correction Sp exists if and only if the sum of all the right-hand-side terms is equal to zero. Although this requirement is satisfied on the finest grid for which the source terms are zero, in a multi-grid context this may not be true on coarser grids. However, the condition may be enforced on coarser grids by a slight modification of the discretization near the boundaries. 4.2.2 Grid Coarsening The method developed in the present work satisfies the continuity when grid coarsening is applied. Each coarse grid continuity control volumes is composed of four fine grid con-tinuity control volumes (see Figure [4.1.a]). It can be shown that defining the restricted coarse grid velocities as the mean of their two nearest neighboring fine grid velocities satisfies continuity: (S?uH + S?vH)\I}J = + * y = 0- (4-33) Continuity is automatically satisfied on the coarser grid. Although this grid coarsening results in no coincident mesh points on any of the hierarchy of grids, it does result in compatible momentum control volumes on all grids. For the boundary layer problem, such coarsening is inadvisable. Brandt [34] recom-mends a semi-coarsening technique using grid H which is coarser than h in only one coordinate (see Figure [4.1.b]). Semi-coarsening involves relatively little work on the coarser grids since the grid is coarsened along one of the coordinates. A combination of block relaxation and semi coarsening may be the best choice. However this technique is Chapter 4. Multi-Grid Methods 36 C.j+1) P H d . J ) uh(U) P h ( i , j ) 0+1 .i) A v H < I J 1 (a) General Coarsening U H(U) „h „ P (U) „h 4 v h (M) A V H ( I , J ) (b) Semi-Coarsening Figure 4.1: A Coarse-Grid Cell Divided into Fine-Grid Cells Chapter 4. Multi-Grid Methods 37 not pursued because of time limitations. 4.2.3 Restriction and Prolongation Restriction is used for the transfer of fine grid values to a coarse grid; prolongation is used for the extrapolation of coarse grid correction to a fine grid. The two operators were denoted earlier by If? and respectively. For the staggered grid, because of the different locations of variables on coarse and fine grids, the restriction is carried out by averaging nearby values (See Figure [4.1.a]). « " j = \\<i + (4-34) = \[vl + t£ij], (4-35) and P?,J = \[Pi,J + rfj+i + J + Awl- (4-36) The prolongation relations are derived using a bi-linear interpolation. For each coarse grid node, four fine grid values are derived. For U-velocity, the four values are: 4j = \^uf,J + u?,J-^ (4-37) «t;+i = \luU+i + 3«£,] , (4.38) < i j = + + u?+i,J-i + 3«f+i,j], (4-39) uJVij+i = + 3 ^ + u f + 1 > J + 1 + 3 t i f + l i J ] . (4.40) The u-velocity can be prolongated by equivalent relations obtained by similar expressions. The pressure prolongation have different weightings because of their center cell locations. Chapter 4. Multi-Grid Methods 38 For each coarse grid cell Pi = ^[9p?j + 3 * - i + Spl^j + p f - L j . J , (4.41) itj+i = jftfa + 3rfrJ+1 + 3 p ? _ M + p ? _ M + 1 ] , (4.42) p!Vi.i = ^ [ 9 P f J + 3 P m . J + 3 A - i + (4-43) p t i , i + i = i t 9 * + 3 A I . J + 3 A + i + pf+iv+ii. ( 4 - 4 4 ) The above relations are slightly modified near the boundaries to avoid the use of the boundary pressures. The prolongation operator for the interpolation of the coarse grid solution to the next finer grid may or may not be the same as the prolongation operator for the corrections. For simplicity, the same operators are used for both correction and solution. More accu-rate cubic prolongation operators can further accelerate convergence but are inconvenient to use in practical flows with complex geometrical shapes [10]., 4.2.4 Boundary Condition Enforcement The treatment of global constraints is discussed by Brandt [5]. For example, in the case of discrete non-linear constraint of the form ch(qh) = 0 (4.45) on a grid-function qh, the following constraint should be imposed on the coarser grid cH(qH) = dM, (4.46) where dH = J j f t - c V ) ) + c"(/*V)). (4.47) Chapter 4. Multi-Grid Methods 39 This procedure for the treatment of constraints is consistent with the FAS procedure. In the S I M P L E scheme two main global enforcements are made: inflow/outflow boundary and pressure enforcement. For the problems with inflow/outflow boundary, there is a constraint on the outflow velocities imposed by the overall continuity constraint In the case of inflow-outflow problems, the discrete compatibility condition does not in general hold since the flow is developing. Even in this case the pressure-correction method used as an iterative solver generally converges. The multi-grid method will almost certainly fail to converge. Since the pressure field in the incompressible Navier-Stokes equations is determined only up to an additive constant, the pressure must also be constrained uniquely. This constraint can be achieved by a pointwise condition However imposing this constraint on all grids does not affect the smoothing properties since this constraint results only in a lowest-frequency error shift. 4.2.5 Calculation of Thermal Energy and Turbulence In the present code, the k — e equations and thermal energy equations ar€ solved every multi-grid cycle but only on the finest grid without the coarse grid corrections. Most difficulties in the convergence of the solution for the turbulent flow with heat transfer are (4.48) In the discrete system, the following compatibility condition must be satisfied ^T(6XU + 6yv)h2 = ^ ( U i n f l o w - " o u t f l o w ) / * + X ) ( U i n f l o w ~ V o u t f l o w ) ^ = 0. p(so,yo) = o (4.49) Chapter 4. Multi-Grid Methods 40 due to the nonlinearly coupled system of the Navier-Stokes equations and the continuity equation. The turbulence equations are coupled with the momentum equations through the effective viscosity terms. The momentum equations are not coupled with the thermal energy equation when no buoyancy is considered. Therefore, when the flow field has been determined, no further convergence problems exist in solving these weakly coupled linear equations. Another reason for not using coarse grid correction for the turbulence equations is that a coarse grid turbulence wall function, consistent with the fine grid flow field, is not easily defined and prohibits the simultaneous multi-grid processing of both the flow and k — t equations. In addition, because the turbulent wall shear stress cannot be adequately resolved on the coarser grid, neglect of its coarse correction on the coarser grid at worst creates a high frequency error that can be eliminated on the fine grid very quickly. 4.3 Performance of Multi-Grid Methods The multi-grid method can be shown theoretically and numerically as an optimal iterative method for solving an increasingly large number of problems. Typically, the cost of each multi-grid cycle does not increase much as the number of unknowns increases. The convergence rate of the multi-grid scheme is not sensitive to the number of mesh points. This property of the convergence rate usually means that multi-gridding can improve the algebraic accuracy effectively. The optimal iterative performance depends on the particular application and multi-grid construction. In this section, the performance and the computational complexity are studied based on assumptions on the accuracy and the stability of the discretization and the grid-independent performance of the basic multi-grid cycle. Chapter 4. Multi-Grid Methods 41 4.3.1 Cost of Multi-Grid Method The multi-grid method requires a marginal increase in storage and computation work over fine grid relaxation. On each grid £lk, a function uk requires about nd locations for a d dimensional problem (d>l). If Q,k is coarsened in powers of two, the total storage requirement is about nd(l + 2-d + 2 - 2 d + •••)= r ^ Z I . (4.50) In particular, for a two-dimensional problem the requirement is less than 4/3 of the fine grid problem alone. The storage cost of the multi-grid algorithm decreases relatively as the dimension of the problem increases. The arithmetic operation of one sweep of relaxation on $7M is defined as one work unit (WU) and the intergrid transfers and other noniteration operations are disregarded. The cost of performing one basic multi-grid cycling with one sweep on each grid is about WU(1 + 2-d + 2~2d + •••) = r ^ p 7 , (4.51) and the cost of one complete full multi-grid cycling is about wu , n_ u , WU i _ 2 ~ d { l + 2 + 2 + • ' - ) = ( 4 - 5 2 ) A n F M G cycle costs about 7/2 of the fine grid simple iteration for a two-dimensional case. 4.3.2 Convergence of Multi-Grid Method The significant advantage of the multi-grid method is that it converges at a rate which is not strongly dependent of the mesh size h. The reason is that the convergence factor for most relaxation schemes is independent of h for the highly oscillatory modes. Through Chapter 4. Multi-Grid Methods 42 the use of coarse grids, the overall convergence factor for a good multi-grid scheme is usually close to the convergence factor of relaxation restricted to the oscillatory modes. A more precise statement about multi-grid convergence is difficult to make. The analysis of any multi-grid algorithm and application can be complex since it must account for many factors, e.g., relaxation and its ordering, prolongation, restriction, and operator coefficients. Brandt [5] presented a Local Fourier Mode Analysis method to predict the smoothing rate of relaxation. Since the reduction of high frequency error components is a local process dependent principally on the local difference nodes, the analysis of this reduction does not need to consider distant boundaries, or varying difference nodes. If the algebraic error is assumed to be a combination of Fourier components exp(\9-x/h), the smoothing factor /2 of the relaxation scheme can be computed through the amplification factor fi(0) of each Fourier component after each relaxation sweep as follows: u = max */2<|*|<T This smoothing factor is an accurate measure of the ability of the relaxation method to damp high frequency error components. This simple quantitative prediction provides the obtainable multi-grid efficiency, an approximation to the asymptotic convergence factor /I" (where v is the number of the sweeps performed on the finest grid per multi-grid cycle). Therefore the smoothing factor may be used to design the interior relaxation smoother, to choose between different relaxation methods, and to judge the performance of the multi-grid scheme. With the two-grid method where the coarse grid problem is solved accurately by per-forming a large number of S I M P L E iterations, a prediction of the smoothing capabilities of the S I M P L E algorithm and hence the efficiency of the proposed multi-grid method can be obtained. The two-grid method has been found efficient with a smoothing rate ranging from 0.5 to 0.9 for Re = 0 to 10,000 (for details, see [35]). Table [4.1] shows the Chapter 4. Multi-Grid Methods Table 4.1: Theoretical Smoothing Factors of S I M P L E Algorithm 43 Reynolds Z^ min /^ max 100 0.630 0.664 1000 0.720 0.935 10000 0.839 0.996 bound of the theoretically predicted smoothing rate of the S I M P L E algorithm. In addi-tion, the smoothing rate for a given Reynolds number does not deteriorate significantly with grid refinement. Thus the indications are that, if the coarse grid problem can be solved accurately by the multi-grid method, it will result in a method that exhibits grid-independent convergence properties over a wide range of Reynolds numbers. However the S I M P L E scheme does not seem to be very satisfactory at large Reynolds numbers, due to the loss of the ellipticity of the difference system (3.24). 4.3.3 Computational Complexity The purpose of a typical calculation is to obtain the approximations qh of discrete system (4.2) which agree with the exact solution u of the continuous problem (4.1) with a specified tolerance e, and with a condition such as \\q—qh\\ < t. The condition \\q—qh\\ < e can be satisfied if it is guaranteed that | |£ f c | | + | |c f c | |<«, (4.53) Chapter 4. Multi-Grid Methods 44 where the global error of the discretization process is Eh = q — Qh and the algebraic error of the iteration process is eh = Qh — qh. It implies that II? -qh\\<\\q-Qh\\ + \\Qh ~ qh\\ = \\Eh\\ + ||e1| < e (4.54) One way to ensure (4.53) is to require that \\Eh\\ < e/2 and ||e f c|| < e/2 be satisfied independently. The first condition determines the grid spacing on the finest grid. If the global error of the dicretization is assumed to be bounded in norm in the form \\Eh\\<chp c G R , (4.55) The choice must be h<h* = (Z)i/ P , (4.56) The second condition ||e^|| < e/2 determines how well our approximations uh must converge to the exact discrete solution Qh, if relaxation or multi-grid cycles have been performed until the condition He^H < e/2 is met on grid Vth, where h < h*. If (4.56) replaces the value of e above, a convergence criterion that directly reflects the computa-tional objective is achieved: \\Qh-qh\\<chp. (4.57) Two problems are associated with this criterion. First, the direct use of (4.57) requires an estimate of c. Second, to satisfy (4.57) with the initial guess | |e^0j| | fixed, it can be expected that the number of multi-grid cycles required , u, must satisfy v = O(-logh). It is supposed that the convergence rate of the multi-grid scheme, 7, is independent of h. With the full multi-grid scheme, a better initial guess, which has already reached the level of coarse grid truncation can be achieved. Hence l|efo)ll = l l« f c-9 ( % l l = ll«*-«(»)ll. (4-58) Chapter 4. Multi-Grid Methods 45 <\\Qh-QH\\ + \\QH-<&)\\ < II?-^ 11 + 119-^ 11 +11^ -95)11' < chp + cHp + cHp = c(l + 2rp)hp. This will require v V-cycles, where 7" = 1 + 1 2 r P or v = 0(1). The computational cost of the F M G method is therefore 0(Nd) which is of optimal order. In practice, this optimal order may not be reached because of the efficiency of the smoother, the order of the prolongation and restriction, and coarsening. However, some corroboration of these results can be provided by numerical experiments. 4.3.4 Closing Remarks In the computation of film cooling there are several limitations in the design of an efficient multi-grid scheme. One is that the coarsest grid is still not coarse enough, since the geometry of the slot is too small and at least two grid points and at most sixteen grid points (in our present study) have to be arranged in the slot for the coarsest and the finest grid respectively. That means that only three levels of grid are allowed in the multi-grid cycle. Since this limitation prohibits reduction of lower frequency errors, the multi-grid scheme cannot perform in an optimum manner. The second limitation is that the grid on each level cannot be designed as a uniform grid because of the boundary layer and small geometry of the slot. The use of non-uniform grid introduces truncation errors in the difference discretization and also results in unnecessary refinement in some regions as a result of the design of the non-uniform grid itself. A local mesh refinement technique would be useful for this reason. Chapter 5 Local Mesh-Refinement Techniques Since in practice, fine grid resolution and accuracy are usually needed only in a local subdomain, the standard multi-grid scheme (which acts as a global grid solver) is in-efficient. This chapter presents a multi-grid local mesh-refinement method, F A C , first presented by McCormick in 1984. The fundamental treatment of the F A C scheme is described. For simplicity, the two-level version of the scheme is studied then extended to the multi-level. The interface treatment of F A C is applied to the present problem. Finally, the performance of the method is discussed. 5.1 Basic Principle of F A C Scheme The multi-grid local mesh-refinement method can be used if a certain accuracy is re-quired for correct resolution of the system (4.1) in most of the global domain fi and a substantially greater accuracy is needed in a given local region f2/, ft; C ft (see Figure [5.1.a]). The domain can be divided into three partitions: ft = ftcUft/UftF, (5.1) where coarse grid partition ftc = ft\ftj (ftj = ft; U (9ft), interface grid partition ft/ = dtti\d£l and fine grid partition ftF = ft;. If it is assumed that a resolution with a refined mesh size h is needed for the local phenomena in ft/, then the fine and coarse grids can be constructed globally as (see Figure [5.1.b]: v ft* = ft£ u ft) u ftF, (5.2) 46 Chapter 5. Local Mesh-Refinement Techniques 47 and nH = n% u n f u n £ . (5.3) The difference systems on both coarse and fine grids by a block decomposition of the operators and grid functions are: LHQH = pH ( TH TH ^cc -° CI 0 0 TH TH TH L'lC ^11 L'lF TH TH ( Qf QF ) F» on nH (5.4) hr\h LhQ rh rh ^CC ^ CI rh rh rh L'lC ^11 \ e Lh FI e \ (Qh) — F h / K Q F J on nh (5.5) The main goal is to perform fine grid computation just on the local domain f2 F instead of the global domain Clh to avoid unnecessary relaxation, intergrid transfer and residual calculations outside of flp-UJlj. The Fast Adaptive Composite Method (FAC) developed by McCormick [40] provides a practical local grid-refinement treatment. 5.1.1 Construction of Composite Grid In the F A C process, a composite grid (see Figure [5.1.c]) is used instead of the global fine grid where = coarse grid points away from the refinement region; f i j = Oj 7 = coarse grid interface points and (5.6) Chapter 5. Local Mesh-Refinement Techniques n ft, an, (a) Domain ft with Local Region ft, • • • • • • • • • • • • • D fl • • • • • • • •c • • o o 0 o • • • • • • • o • • • • • o • o • • • • it • F • • o • • • • • • • o • • • • • • r • • O o nf o i • r o • • • o • (b) Coarse Grid ftH • n b • • • • O Clh o o • \ i j O • • • • • • • I • • • p • • • • • • (c) Fine Grid QH (d) Composite Grid ft^-Figure 5.1: Construction of a Composite Grid Chapter 5. Local Mesh-Refinement Techniques 49 ft nhF fine grid patch. The difference systems on the composite grids is: / rk LkQh = Fn CC ^CI rk rk i-iir. i-Ji\ \ J1C 0 h FI IF rh rh J-IVT LJFF ) Or Ff 1-1 h on Vtk (5.7) The composite grid operator agrees with the coarse grid operator on Q,Q (i.e., LQC = Lcc) a n d w i t h the fine grid operator on (i.e., LpF = LhFF). LFF is similar to LQC except that there are more grid points and therefore more entries in the LFF block. The difference between LH and L- is in blocks L ; F and Ljt which are the important parts of the composite operator L-. 5.1.2 Two-level Nonlinear FAC Scheme The F A C scheme performs as a natural extension of multi-grid method to composite grid equations. It uses the coarse grid H to approximate the smooth global component of qh and fine grid h to approximate the oscillatory local components. The two-level FAC scheme is processed on the composite grid ft- with the locally refined grid Qf. A F A C cycle that allows for non-linear operators and is similar to the general multi-grid scheme discussed in section 4.1.3. can be written as: • Step 1. Coarse grid relaxation Relax on the coarse grid ft": qH «- LHQH = lJ[{F- - l£q*) + LH(l£q±), where If? is the restriction operator from the composite grid to the coarse grid. • Step 2. Coarse grid correction Prolong correction and add to composite grid ft-: q- = q-+ 7^(g H ' - Jjfg-), where Ijj is the prolongation operator from the coarse grid to the composite grid. Chapter 5. Local Mesh-Refinement Techniques 50 • Step 3 . Local refined grid relaxation Relax on the local refined grid ftF: qp <— LpQp = Fp — Lpj(qj). Step 1 augments the right hand side with the residual to correct the most recent ap-proximation to the composite grid solution. Step 3 solves the system on ftF using the previous approximation for so that the information is transferred through the internal boundaries. It can be seen that, in the F A C processing, coarse grids are used for fast convergence of the fine grid solution and transmitting fine grid accuracy throughout the domain. The F A C scheme acts as an iterative mesh-refinement method. The classical mesh-refinement methods start by solving the coarse grid equation, then solving the fine grid equations using the coarse grid solution to define Dirichlet boundary conditions at the in-terface. The problem is that the results of this two-grid process usually do not achieve an accuracy commensurate with the added resolution. The F A C scheme uses the composite grid equations and, with the coarse and fine grid processing steps, to obtain accuracy more effectively. 5.1.3 Multi-grid FAC Scheme The multi-grid F A C algorithm can be obtained from the two-grid F A C algorithm. If a family of locally nested grids ft*1, ft*1,..., ft*M covering the nested region ftl5 ft2, • • •, (in the definition, ft*' fl ft*,+1 C ft*' and ft, C = 1 , 2 , . . . , M - 1 ) and if h = ht signifies coarsest level with successive mesh sizes differing by a power of 2 , the resulting composite grids are: ft*± = nh>uft*-1 u...unh i (t = I , 2 , . . . , M ) . ( 5 . 8 ) A general V—cycle F A C algorithm on level i is given in the flow diagram (Figure [5 .2 ] ) which also shows the structure of the present code. The relaxation sweeps from ftF_1 Chapter 5. Local Mesh-Refinement Techniques 51 ^ e ^ i n ^ k=i Relax on n1 times Flux Correction on ftj* Restriction A4*-1 : Ci>u p>-. k>1 I k-k-1 Z3Z Relax on coarsest grid n2 T times Flux Approximation on fif k - 1 Prolongation J. Relax on n£+1 n3 times k< i T Flux Approximation on ny*+2 k-k+1 Figure 5.2: Flow Diagram of Multi-Grid F A C Processing Chapter 5. Local Mesh-Refinement Techniques 52 A n Slave Point Interface Point Fine Grid Point —0 Coarse Grid Point Figure 5.3: Composite grid interface to f l F are referred to as 'coarse grid corrections'. The residuals on the finer grid points are restricted to the coarser grid, and interface grid points are corrected through the composite equations. The relaxations from Clp1 to f l F - 1 provide the approximation for the refined grid. The interface grid points are interpolated from the coarse grid points. 5.2 Interface Treatments A most important part of F A C implementation concerns the treatment of the grid inter-face because it controls the transition between the two grids. The treatment is guided by the composite grid equations although the equations themselves need not be used explicitly in the algorithm. The only significant impact of non-uniformity on the compu-tation is the interpolation of the slave point values (see Figure [5.3]) and the evaluation of certain residuals at coarse grid interface points. Chapter 5. Local Mesh-Refinement Techniques 53 5.2.1 Interpolation of the Slave Grid Points On the staggered composite grid, the slave points values are usually obtained by linear interpolation from the coarse grid (although higher order interpolation is preferred). These points are used as the boundary conditions for the local computation (see Figure [5.3]). Since the S I M P L E pressure correction scheme is used as the direct smoother on the nested grids and since a unique solution of the pressure correction equations exists if and only if the mass is balanced in the computational domain, the mass must be conserved on the locally refined region. The mass balance on the slave points should be corrected before entering the next locally refined grid to eliminate the error within the interpolation. In the present scheme, the mass balance along the boundary of the local region is calculated and the correction is distributed to each point on the boundary so that the mass is balanced. 5.2.2 Evaluation of the Residuals at Coarse Grid Interface Points After linear interpolation and mass correction, some errors introduced at the slave points will also affect the equation. One important step is to correct these slave point values by the correction equation with the residuals being evaluated at the coarse grid interface points with the slave point values. For example the mass conservation at point P on the grid interface in Figure [5.3] can be written as 4 = pH[u» - a? + - (vl + vhJ/2] (5.9) because the velocity V at points sw and se on the bottom side of the control volume are interpolated from the coarse grid. The bottom side flux correction is essential to avoid mass unbalance in the processing of the local refinement. In this scheme,-the calculation of the fluxes on the coarse grid is referred to as "flux approximation" and that on the local Chapter 5. Local Mesh-Refinement Techniques 54 refined grid as "flux correction". The same correction can be made for the momentum flux, heat flux and other physical variables. 5.3 Self-Adaptivity of FAC With self-adaptive schemes, FAC can be organized so that the decision to refine further a particular grid can be made after that grid has been processed. Two self-adaptive schemes can be used in the F A C algorithm: truncation error estimate (r-extrapolation) and absolute error estimate (Richardson extrapolation). 5.3.1 Truncation error estimate In FAS nonlinear multi-grid scheme, the fine-to-coarse defect correction TJ? (4.13) r" = LH(I»qh) - I?(Lkqh), (5-10) serves as an approximation to the local truncation error TH on the coarse grid H which is defined by r» = LH(IHq)-IH(Lq), (5.11) where q is the true differential solution and qH is the approximation solution to q. It is clear that at convergence T " « T f c + Tf, (5.12) where rh is the fine-grid local truncation error and the sign « means equality up to a higher order in h. If the local approximation order at the point x is p, i.e. rh(x) = hvC(l\x) + h^C(2\x) + . . . , (5.13) Chapter 5. Local Mesh-Refinement Techniques 55 where C^\x) and C^2\x) is independent of h, then T H ( X ) = HPC^(x)-rH^1C^(x) + ... = m'h'CMix) + m'+1h'+1CW{x) + ..., (5.14) where m is mesh ratio between grid H and / i ( m > 1). Hence r f (x) » (m p - l )C ( 1 ) (x ) f t p (5.15) and T H ( X ) « m p ( m p - l ) _ 1 r f (x). (5.16) This operation which is called "T—extrapolation" allows the determination of the relative local truncation error rf? in an inexpensive way. This local truncation error of grid H relative to grid h can be used for the local truncation error estimate and then be used as Full Multi-Grid (FMG) stopping criteria and in grid adaptation criteria. 5.3.2 Absolute Error Estimate Obtaining the absolute error at a given point for each of the two grids is similar to obtaining the truncation error estimate. If the errors at point x on the coarser grid H and finer grid h are given respectively by (m p - l)-*7f, (5.17) eH = q _ q n = HrFW(x) + Hp+1F(2\x) + ... (5.18) and eh = q-qh = hpF^{x) + hp+1F^2\x) + ..., (5.19) Chapter 5. Local Mesh-Refinement Techniques 56 where F^(x) and F^2\x) are general undetermined coefficients, subtracting (5.18) from (5.19) gives eH - eh = qh - qH = (mp - l)F{1\x)hp + (mp+1 - l J F ^ x ) * ' * 1 + . . . . (5.20) Comparison equation (5.19) and (5.20) gives ^ + 0 ( A ' + 1 ) . (5.21) so that mP- 1 v ; is (k + l)st order approximation to the error and determine points at which refinement is needed. 5.4 Performance of FAC Algorithm As an iterative solver of the composite grid equations, F A C cycles have essentially the same computational complexity as the solver used for global multi-grid equations. Main-tenance of the linear independence of the iteration on the number of gridpoints requires improvement of the scheme and theoretical analysis. This section presents the complexity of the F A C algorithm and the convergence improvement of the scheme. 5.4.1 Complexity of FAC Scheme If a full F A C algorithm is considered and if for each composite grid ft- = UJ t = 1ft' 1* the discretization satisfies the relation \\Eh-w<ci:hi k=i Chapter 5. Local Mesh-Refinement Techniques 57 where E- = q — Q- is the global error in the composite grid solution, p > 0 is the order of approximation, c is a constant independent of h. To obtain accuracy to the level of truncation, the convergence criterion for the com-posite grid f i - must satisfy where e- = Q-— q^ is the algebraic error in the i/th F A C iteration as an approximation to UK As with the full multi-grid scheme, an initial guess q^ is attempted. A grid signif-icantly coarser than $7- is required (so that computation is relatively inexpensive), but one that delivers reasonable discretization accuracy. If the approximation g ^ , which is the next coarser composite solution in vih iteration, satisfies the convergence criterion (5.4.1 Relatively to its level I I ^ I I <ci>r,, jt=i h H * then the initial guess g ^ = g^for grid h satisfies kfo)H = H ^ - 9 f o , l l = WQ-- (5-23) < k - ^ l l + l k - Q - | l + l l ^ - - ? f t l l <(l + 2r')cJ>Z, k=i where r = hk\' Hk(k = 1,...,/) is composite grid mesh ratio. In general, the convergence rate of F A C scheme fiFAC < 7 < 1 is sufficient to perform v = log(l + 2r p)/(—log7) cycles of F A C on each composite grid. 5.4.2 Convergence of FAC Scheme Multi-grid local mesh-refinement methods assume that the convergence rates are not sensitive to the number of grid points on each refinement level. Direct solvers produce Chapter 5. Local Mesh-Refinement Techniques 58 approximations accurate to the level of truncation error at a cost proportional to the number of composite grid points (see [14,41]). Numerical experiments can be used to study the performance of F A C methods, as theoretical information on the convergence behavior of the F A C algorithm is scarce, especially for practical problems. When the full multi-grid algorithm is used for solving a regular problem with all grids covering the entire domain, the work invested on the coarse level is usually less than that on the finest one despite the fact that the coarser grids are relaxed many more times. The explanation is that the coarser the grid , the fewer its points. As a result, total work is proportional to the number of grid points in the domain. On the other hand, when local grid refinement is used with finer levels covering much smaller subdomains, the number of points on the coarser levels is not necessarily smaller in comparison. More work may be invested on the coarse grid than that on the locally finest grid. Therefore the work for the coarse grid correction cannot be ignored compared to the work on the finest grid. For this reason the regular F A C scheme may not perform as efficiently as the global multi-grid methods. In the present work, no theoretical determination of the ideal convergence rate of the F A C scheme for the Navier-Stokes equations has been carried out. Only numerical exper-iments are presented. The following suggestions are made to optimize the performance of present F A C code: (1) Optimization of the mesh sizes of local regions: Brandt introduced the A-FMG algorithm [14] which uses the so-called mesh optimization equation for optimizing the mesh size h by minimizing the solution error under fixed invested work. (2) Increased order of interpolation: Since the interpolation introduces a higher order of local error than that of the local truncation error, higher order interpolation can be used. (3) Choice of coarsening: The iteration on the local region can be accelerated by the Chapter 5. Local Mesh-Refinement Techniques coarse grid on the local region itself. Chapter 6 Numerical Results This chapter evaluates the efficiency of the multi-grid method in the calculation of the Navier-Stokes equations for recirculating flows. The performance of the multi-grid scheme for driven cavity flow, which is well known as model test of numerical schemes, is illus-trated. A local grid refinement applies the F A C method to a region where greater res-olution is needed. This chapter also gives the results of film cooling computations with the multi-grid method. The influence of refining the computational grid on the predicted flow pattern and the film cooling results are investigated. Coarse grid and refined grid results are compared. The effects of the turbulence wall treatment on the computation is discussed. 6.1 A Model Test - Driven cavity flow The numerical solution of the flow in a driven cavity is studied to validate the efficiency of numerical algorithms for the Navier-Stokes equations. A viscous and incompressible fluid is contained in a square cavity. The problem, namely the determination of the flow induced by the steady tangential shearing motion of the top wall, characterizes the elliptic and non-linear nature of many engineering flows. Numerical studies on the flow structure at high Reynolds numbers (Re > 5000) have been published [9,10]. Several finite difference grids have been considered to test the present code for driven cavity flow at different Reynolds numbers. Calculations have been made for Reynolds numbers of 100, 1,000, and 10,000 with grids consisting of 34,66, and 130 finite difference 60 Chapter 6. Numerical Resnlts 61 Figure 6.1: Locally Refined Gird nodes. The boundary conditions are of the Dirichlet type for velocities. No boundary conditions were necessary for the pressure (the system is solved by referencing the pres-sure to the value at one comer node). The under-relaxation factor for the momentum equations ranges from 0.5 to 0.1 depending on the Reynolds numbers of the flow. At higher Re lower under-relaxation factors must be used. The convergence criterion was based on the maximum residual in the three equations, i.e., # c o n v = max(r^ , r^ . , r f j ) . The momentum residuals have been normalized by pU^ and the mass error is normalized by pUo, where Uo is the top wall velocity and p is the density. The convergence criterion Rconv is set to 10~4 and, when the fine grid (or the composite grid) residuals decrease below this value, the calculations are terminated. The number of relaxation sweeps on each grid is fixed in the present scheme. In multi-grid processing, an estimate of the local truncation errors is automatically obtained. The truncation errors are large near the region of the moving boundary. Local Chapter 6. Numerical Results 62 Table 6.1: 10" 4 Number of Iterations an d C P U Times Required to Reduce Residual Norm by Re Finest Grid (Work unit / C P U time) 18x18 34x34 66x66 130x130 100 20.3 / 0.8s 19.9 / 2.75s 23 / 13.7s 30 / 1.72m 1000 41.1 / 1.5s 59 / 8.4s 75 / 45s 67.3 / 3.73m 1 0 0 0 0 60.8 / 2.2s 80.9 / 12.3s 107/ 1.02m 129 / 7.02m mesh refinement is applied in this region (Figure [6.1]). The multi-grid local mesh-refinement method improves the efficiency of the numerical solution on the composite grid with only little added computational cost over the coarse grid solution. The calculation is performed on an IBM3090 machine with a F O R T R A N compiler with opt(3) optimization. Tables [6.1] - [6.3] summarize the convergence characteristics of the multi-grid method and the local refinement method together with the C P U times and compare their per-formance against T E A C H code (a widely-used single grid code based on the S I M P L E algorithm) [36]. As indicated in Table [6.1], for low Re, the number of work units re-quired for the solution do not increase significantly with the grid size. The efficiency of the smoother is lower at higher Re, due to the loss of the ellipticity of the system. Table [6.2] shows that, for the whole range of Re, the multi-grid scheme is significantly faster than the standard T E A C H code. Table [6.3] compares the convergence rates of the M G Chapter 6. Numerical Results Table 6.2: Comparison of Multi-Grid Performance with T E A C H Code. 63 T E A C H Code Multi-Grid Scheme Re ( Work units / C P U times) 100 506/4.14m 23/13.7s 1000 >1000 />8.9m 75 / 45s 10000 >1000 / >11.43m 107 / 1.02m Table 6.3: Asymptote Convergence Rates of F M G , F A C and T E A C H codes F M G TEACH FE (18x18) (66x66) ( 66x66 ) 100 0.60 0.73 0.98 0.77 1000 0.71 0.76 0.988 0.82 10000 0.84 0.87 0.999 0.83 Chapter 6. Numerical Results 64 U / Uref Figure 6.2: i7-velocity Profile along Vertical Centerline, Re = 100 scheme, the F A C scheme, and the T E A C H code, where the convergence rate is defined as the reduction factor of the residual error after each iteration. The efficiency of the local refinement method is lower than that of non-linear global multi-grid methods. Figures [6.2], [6.3] and [6.4] show the (7-velocity distribution along the central line of the cavity on the coarse, locally refined, and globally refined (co-extensive) grids. The first advantage of the locally refined grid is that it not only resolves the local region more accurately near the plate in which the finer grid is used but this accuracy also reaches out to the rest of the region, to which only the coarser grid is applied. A second advantage is that computational time and memory space are saved. The computation cost for the local mesh refinement processing is about 1/4 of the co-extensive multi-grid processing. Figures [6.5], [6.6] and [6.7] show the flow pattern contours at different Re on the different grids. The finest grid 130 x 130 clearly shows the secondary vortices of the large Re (= 10,000) number cavity flow. The F A C scheme is applied to refine the coarse grid (34 x 34) locally near the moving boundary. The results are compared with the Chapter 6. Numerical Results 65 Figure 6.4: /^-velocity Profile along Vertical Centerline, Re = 10,000 Chapter 6. Numerical Results (c) Global grid (34x34) (d) Locally refined grid of (c) Figure 6.5: Streamline Pattern for Cavity Flow, Re = 100 Chapter 6. Numerical Results (c) Global grid (34x34) ( d ) Locally refined grid of (c) Figure 6.6: Streamline Pattern for Cavity Flow, Re = 1,000 Chapter 6. Numerical Results (a) Global grid (130x130) (c) Global grid (34x34) Figure 6.7: Streamline Pattern (b) Global grid (66x66) (d) Locally refined grid of (c) for Cavity Flow, Re = 10,000 Chapter 6. Numerical Results 69 coarse grid and the co-extensive grid (66 x 66). The secondary vortices can be found on the left-top corner in the numerical results but cannot be seen on the figure because of plotting difficulties of the composite grid results. Chapter 6. Numerical Results 70 6.2 Numerical Computation of Film Cooling Flow 6.2.1 Computational Grid Structure The selection of the numerical grid affects not only the accuracy but also the convergence of the numerical solutions. The present work uses a non-uniform grid for film cooling computation because of the existence of the small slot and boundary layers. Both x— and y—direction grids are refined locally around the slot and in the separation region. This grid arrangement is aimed at the reduction of the truncation errors of the discretization as well as at a more precise prediction of flow separation. Figure [6.8] shows the typical coarsest grid and finest grid distribution used for the film cooling flow domain. This distribution structure is obtained from previous computations [3,4]. Only six or eight grid points were located within the slot in the previous computation and the point next to the wall was located at 1/8.S. In the present computation, the grid is refined, 16 points are located within the slot, so that the next wall point is located at l/16s. 6.2.2 Performance of Multi-Grid Method A fixed multi-grid computation is performed on the non-uniform grids. The convergence criterion is based on the values of absolute residual errors of the governing equations on the finest grid and is set to 1 0 - 4 . Figures [6.9]-[6.14] show the convergence behavior of the multi-grid processing and the T E A C H code applied to each of the governing equations. Four-level grids are used in the multi-grid cycle and the T E A C H code is run only on the finest grid. The multi-grid processing produced much faster convergence (more than 10 times) than the T E A C H code. Table [6.4] shows the convergence rate for solving each governing equation for a Reynolds number of 1,800 (Re = pUs\ots/fi). The convergence rate obtained for the Navier-Stokes equations for laminar flows with simple geometries, such as driven cavity Chapter 6. Numerical Results 71 (a) Coarsest Grid (N^ot — 2) • • • • • • • • • • • i i i m i i i i m i i i t i i •••liiiiiiiiiiiiiiiiiiiiiiint l l l l l l l l l i l l i i l l l i i l l l i i i m i n I I H I I I I I I I I I I I I I I I H I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I M I I l l l ••••••••••lllllllllllllllllllll ••••••••••••IIIIUII1IIHI1I1I • • • • • • • • • • • • • • i i m i i i i m i m i • • • • • • • • • • • • • i i i i i i i i i n i i i m i • • • • • • • • • • • • • • i i i m i i i i i i i i i i i • • • • • • • • • • • • • • • • I IMIIIII t l l l l l •••••••••••••••••••llllllllllil • • • • • • • • • • • • • • • • • • • • I I M l l J l l l l • • • • • • • • • • • • • • • I f l l f m i l l M I I I I • •••••••••••••••••••tiimiiiii ••••^••••••••••••••••••••lltlll •••••••••••••••tiiaiiiiiiNiiii • • • • • • • • • • • • • • • • • I tll«lllllllll • ••••^•••••••••••••••untHiiii • • • • • • • • • • • • • • • • • • • • • • t i m i m mm ii- • lllfilir : ".'iiiimin MUM' ' iniiiiiiiiii HIHll 1 ii: 1111111)11 Itl.'ll II 1 • r inn Jin (Iti'lM <i imiiifii HHU • ••\umm iiiiini t::i IIIINIII llli'lM -• •'•i I1I1IIMI min ii • • n iiimm mum TUllll l l l l l l II:M! . MMIIIIIKIII IJ II 11 ' • •• i-mnmi IIIIIIU' ' ' "tilHIJlltl •III • 1 llllltllll HII'I H • < 1 lllllltll HHtl ! 1 - • 1 Miimm - iiuuiiu < Ml lllllltll • "'! "!!!"" IIIIIIIUI I IUIIHIB IIIIIIIUI i m i i i i a a IIIIIIIIN u m i a a a a IIIIIIIUI i i a i i m i i I I I I I I I U I n i m i a a a i i i i i i a a a a _ _ l (b) Finest Grid (Nslol = 16) Figure 6.8: Computational Grid for Film Cooling (a Chapter 6. Numerical Results Single Grid • • • • Hultl-Grld e.ei — 1 CO 68E-3 * i.aeE-4 * 1.86E-5 Hork Unit Figure 6.9: Performance of Multi-Grid Method ([/-momentum Equation) _ 1.B8E-3 * in D "O VI dl QL i.esE-4 1.88E-5 Work Unit Figure 6.10: Performance of Multi-Grid Method (V-momentum Equation) Chapter 6. Numerical Results Figure 6.11: Performance of Multi-Grid Method (Mass Continuity Equation) Figure 6.12: Performance of Multi-Grid Method (Thermal Energy Equation) Chapter 6. Numerical Results 74 Work Unit Figure 6.14: Performance of Multi-Grid Method (e-equation) Chapter 6. Numerical Results 75 Table 6.4: Convergence Rate of Multi-Grid Method for F i lm Cooling Computation Equation U V M T K E Convergence Rate 0.91 0.93 0.89 0.84 0.85 0.87 flow, is not attained in this case. One reason is that the coarsest grid is not coarse enough. Another reason is the extra-coupling between the momentum equations and turbulence equations through the near-wall turbulence function and effective viscosity values. However, the code is efficient enough for practical purposes. Figures [6.12]-[6.14] show the convergence of the relaxation for the turbulence equations and the energy equation without coarse grid correction although each equation is calculated on each level. A significant amount of acceleration in the convergence of the solution of energy and turbulence equations is produced by the multi-gridding of the momentum and continuity equations. The flow and energy equations are not coupled, and the flow and turbulence equations are weakly coupled. Therefore, once the flow field reaches convergence, no further difficulties exist in the calculation of energy and turbulence equations; hence no introduction of multi-gridding to the energy and turbulence equations is necessary. In the calculation of turbulence, the accuracy of the numerical method is very im-portant to the assessment of a given turbulence model. In the classical T E A C H code, because of the slow convergence of the Navier-Stokes equations, the algebraic error can-not be lowered to a reasonable level and convergence difficulties are experienced with the turbulence e equation. These difficulties are not found when the multi-grid method Chapter 6. Numerical Results 76 is used. Multi-grid processing offers a good approach to assess the turbulence mod-elling performance and thus is of critical importance in the computation of film cooling effectiveness. 6.2.3 Non-uniformity of the Slot Injection Figures [6.15]-[6.19] show the predicted variation of magnitude and direction of the ve-locity at the interface with coarse and refined grids. The hybrid difference scheme(HB), which is first order accurate, is used on both the coarse grid and the refined grid. The skewed difference scheme(SK), which reduces the false diffusion [27], is used only on the coarse grid. The flow pattern at the slot exit shows more non-uniformity with higher resolution. The flow distribution at the corner of the slot is crucial to the prediction of the separation downstream of the slot. Grid refinement, therefore, can provide a more accurate flow field for the computation of heat transfer along the wall surface. The vector field shows that the flow in certain places is inclined to the grid mesh. False diffusion of the difference scheme can reach an unacceptable level in this case. False diffusion results in smearing the gradients of flow and temperature fields. Figure [6.20] shows smearing of the temperature gradient by false diffusion. Even with the low order difference scheme, grid refinement reduces false diffusion errors. er 6. Numerical Results Chapter 6. Numerical Results 78 ^ ^ ^ ^ ^ ^ ^ *,0*ffff ^ *tftttt1\ t A i 4 4 4 it tttttM Figure 6.17: Flow Pattern at the Interface of Slot Injection ( M = 0.2,Refined grid(HB)) S S / / Figure 6.18: Flow Pattern at the Interface of Slot Injection ( M = 0.2,Coarse grid(SK)) Chapter 6. Numerical Results 79 / S / z0 ^ . ^ / S / / ^, ^ * * s / f / Figure 6.19: Flow Pattern at the Interface of Slot Injection ( M = 0.2,Coarse grid(HB)) CO o o X/8-0 0.0 0.2 T / Tref e 0.0 0.8 0.0 0.6 0.0 0.6 0.0 0.6 P - l 1__ 1 _J 1 _ 1 1 X/8-1/J , x/t-2/ 0.4 0.6 I 0.6 X/S 1.0 1.2 CO Refined Grid(HB) • Coarse grid(SK) • Coarse grid(HB) Figure 6.20: Temperature Distribution across Slot Exit (M = 0.2) Chapter 6. Numerical Results 80 2-i i i i i i i i i < 0.0 S.6 6.0 T.6 10.0 12.6 16.0 17.6 20.0 22.6 26.0 X/8 Figure 6.21: F i lm Cooling Effectiveness (M = 0.2) 6.2.4 Film Cooling Performance The film cooling effectiveness is computed for a mass ratio of M = 0.2. Figure [6.21] shows the comparison of film cooling effectiveness computed with the coarse grid and with the refined grid. On the coarse grid, the cool slot-injection is diffused numerically and the effectiveness at the exit corner does not reach 1, because of the false diffusion of the hybrid scheme near the slot. On the refined grid the false diffusion is greatly reduced and the effectiveness reaches 1 which coincides with the experimental observations. Figures [6.22] - [6.24] show the longitudinal velocity, turbulence intensity and temper-ature distribution downstream of the slot. The reference velocity and temperature in the figures are the main flow velocity and temperature respectively. Although the prediction of the effectiveness on coarse grid with skewed scheme reaches 1 at the leading edge of the slot, the coarse grid solution is unable to predict the separation downstream of the slot. The separation is found on the refined grid solution where a recirculation region is located between 1 < x/s < 1.62. The difference of the mean velocity between the coarse Chapter 6. Numerical Results U / Uref CO o o.o o.e o.o o.e o.o o.e 1 1 ' _ ' L _ 0.0 0.6 0.0 0.6 s * 1 I i s i i i i i 1 i I i I 1 X/S=0| X/8=1 | X/8=2 i 1 i m % % 1 1 X/8=3 g X/8=4 g % % 1 i i c i 4 f < < i f * i • * 1 1 — • • i i # e CO 0 0 0.5 1.0 1.6 2.0 2.6 3.0 3.6 4 .0 4.6 6.0 X/S U / Uref ° ° ° . e ° ° o.e o.o o.e o.o o.e o.o o.e o.o oe ' ' 1 -m-i 1 » J 1 . I i _ CO o to d o e . X/8=5^ X/8=6£ X/S= X/8=tU X/8=&f < , 9 X/8-1|J e CO I * 1 ¥ 1 f — 1 9 1 9 w | 6.0 6.6 6.0 6.6 7.0 7.6 8.0 6.6 8.0 8.6 10.0 10.6 11.0 X/S • Refined Grid(HB) • Coarse grid(SK) + Coarse grid(HB) Figure 6.22: Longitudinal Velocity Distribution Downstream of the Slot (Af = 0.2 Chapter 6. Numerical Results uV Uref 0 0.0 0.1 0.0 0.1 0.0 0.1 0.0 0.1 oi-r-» L 1 L - 1 L_ L _ CO e o e. 0.0 0.1 0.0 0.1 I , I L— I i i i i £ x/8=5 £ x/8=6 £ x/s=7 £ x/s=8 £ x/s=9 i X/8=10 • K % % % < I 5 I < i { 1 i e d o CO n r — * 1 f 1 f I f -6.0 6.6 6.0 6.6 7.0 7.6 8.0 8.6 9.0 9.5 10.0 10.6 11.0 X/S o 0.0 0.1 0.0 0.1 uV Uref 0.0 0.1 0.0 0.1 i i 0.0 0.1 1 1 w ~ 8 I I i i to $ I t i i i i i I * < c* I < B * B | X/8=0 X/8=1 m9 X/8=2 B « X/8=3 X/8=4 CO e • • • i i I i •# I V to e»- % V \ < o T < , % • • * * 0.0 0.6 1.0 1.6 2.0 2.6 3.0 3.5 4.0 4.6 6 o o CO X/S • Refined Grid(HB) • Coarse grid(SK) 4- Coarse grid(HB) Figure 6.23: Turbulence Intensity Distribution Downstream of the Slot (M = 0.2) Chapter 6. Numerical Results CO o o b T / Tref o.o o.e o.o o.e 0.0 0.6 0.0 0.6 0.0 0.6 8 8 8 8 8 8 8 8 i 8 8 s i 1 8 8 8 8 X/8=0 j X/8=1 J s 8 8 X/8=2 J i 8 8 X/8=3 J $ 8 8 X/8=4 | 8 8 8 i 8 8 i i < 8 : 8 4 / 8 • tv 8 . ° ? , If , 4 -o tO u> © o 6 0.0 0.6 1.0 1.6 2.0 2.6 3.0 3.6 4.0 4.6 6.0 X/S T / Tref _ 0.0 0.6 0.0 0.6 0.0 0.6 0.0 0.6 0.0 0.6 0.0 0.6 0 *H " 1 w-> ' r - » ' w* 1 r - 1 1 m-rd CO o d X/8= X/8=7g i X/8=8g < X/8= X/8=5> 1f / X/8=1{ i i © _ | _ | , . ., d 6.0 6.6 8.0 6.5 7.0 7.6 8.0 8.6 0.0 8.5 10.0 10.6 11.0 CO X/S • Refined Grid(HB) O Coarse grid(SK) + Coarse grid(HB) Figure 6.24: Temperature Distribution Downstream of the Slot (M = 0.2) Chapter 6. Numerical Results 84 M-0.4 M-o.e T I I I I I I I I I T 0.0 2.6 6.0 7.6 10.0 12.6 16.0 17.6 SO.O 22.6 26.0 X/S Figure 6.25: F i lm Cooling Effectiveness at Different Mass Ratios grid solution and fine grid solution is small in the recirculating region, but larger turbu-lence kinetic energy is found near the wall with the refined grid. The underprediction of the turbulence kinetic energy of the coarse grid solution results in underestimation of mixing and heat transfer. Therefore, smaller effectiveness is predicted with the refined grid downstream of the slot, i.e., 1 < x/s < 3.5. Figure [6.25] shows the effectiveness at different mass flow ratios. The cooling effec-tiveness increases almost in proportion to the mass flow ratio at large x/s. Within the near-slot region (1 < x/s < 2), the cooling effectiveness at Af = 0.4 and Af = 0.6 are lower than that at Af = 0.2 because of the stronger separations which occur for Af = 0.4 and Af = 0.6. Figure [6.25] also shows no significant difference between Af = 0.4 and Af = 0.6 downstream of the slot x/s < 5. This means that larger mass flow ratios do not necessarily provide better overall cooling performance. Because, in practical film cool-ing processing, coolant orifices are located on the blade surface close to each other, the performance near each injection hole is very crucial to the overall cooling performance. Chapter 6. Numerical Results 85 6.2.5 Comparison of the Near-Wall Turbulence Treatments The turbulence model near the wall has some but not significant influence on the com-putation. Figure [6.26] shows the film cooling effectiveness for the two turbulence mod-els used. Model 1 uses the shear stress determined on the grid distant from the wall (30 < y+ < 150). It is assumed that the standard wall function can predict better turbu-lence shear stress and that within the inertial layer the shear stress changes little. Model 2 uses the algebraic model to predict the shear stress obtained from the velocity gradient on the wall and the turbulent viscosity. Model 2 produces a separation point closer to the leading edge of the slot (x/s = 1.23) (see Figure [6.27]) than Model 1 (x/s = 1.37) (see Figure [6.17]) and the recirculation region of Model 2 is larger (1.23 < x/s < 4.4) than that of Model 1 (1.37 < x/s < 1.62). These differences are caused by the shear stress correction on the wall. Model 1 uses the shear stress based on a coarse grid which does not resolve the recirculation zone. This causes an error in the wall shear stress calcula-tion. Model 2 is based on the fine grid which resolves the recirculation zone. Therefore the correct wall shear stress is used in Model 2. In Figure [6.26], the effectiveness of Model 2 near the slot (1.2 < x/s < 4.5) is smaller than that of Model 1 because the stronger recirculation found in Model 2 increases mixing and heat transfer. At this stage, no further comparison can be made because the assumptions are still controversial and the experimental work in support of the modelling has not been yet carried out. Chapter 6. Numerical Results 86 to ° CO IU z > 2 p ° o t. LEGEND ° MODEL 1 MODEL 2 ° T i i i i i i i i i 0.0 8.6 6.0 7.6 10.0 t2.S 18.0 17.6 80.0 82.6 U . 0 X/S Figure 6.26: Fi lm Cooling Effectiveness for Different Near-Wall Turbulence Treatments Figure 6.27: Flow Pattern Downstream of Slot Injection (Af = 0.2, Model 2) Chapter 7 Conclusions and Recommendations The multi-grid method solves efficiently the non-linear strongly-coupled system of Navier-Stokes equations and provides fine resolution in flow computation. The performance of the multi-grid scheme tested on the driven cavity flow shows that the convergence rate is not sensitive to the number of grid points and flow Reynolds number. The local mesh-refinement method applied to cavity flow is a stable and efficient technique to increase the local resolution. The results obtained for cavity flow agree with those reported in the literature. The present computations of film cooling show the multi-grid method to be almost ten times faster than the T E A C H code. Accurate flow fields are provided for the calculation of the turbulence and energy equations. A strong non-uniformity of the flow at the slot exit is found. As this non- uniformity affects the film cooling process, it needs to be accounted for. The near-wall turbulence treatment is critical to the prediction of heat transfer near the wall because the cooling effectiveness needs correct prediction of turbulence kinetic energy. The multi-grid method provides refinement near the wall and makes the accurate prediction of turbulence kinetic energy possible. However the present turbulence modelling near the wall is still not fully satisfactory. The multi-grid method also provides increased accuracy and convergence. This allows for the assessment of the performance of different turbulence models. The multi-grid method contributes significantly to the reduction of false diffusion. The following further work is recommended: 87 Chapter 7. Conclusions and Recommendations 88 (1) Higher order schemes should be introduced in combination with the multi-grid scheme in order to further improve accuracy. (2) Higher order interpolation should be introduced in the local grid- refinement scheme at the interface between the locally refined region and the rest of the domain in order to improve accuracy and increase convergence. (3) A coarse grid correction of the turbulence equations and the heat transfer equation should be introduced to increase further computational efficiency. (4) Further work should be carried out in turbulence modelling for film cooling. The high accuracy of the multi-grid method allows for a better assessment of the performance of the different turbulence models. Low Reynolds turbulence models should also be considered. (5) The film cooling computation should be extended to two-dimensional inclined cooling orifices using non-orthogonal numerical grid generation. (6) A three-dimensional multi-grid code should be developed in order to investigate the film cooling on the blade by using a more realistic geometry. (7) Experiments should be undertaken to validate the numerical code and to improve the turbulence modelling. Bibliography [1] R . J . Goldstein, "Fi lm Cooling", Adv. Heat Transfer, 7, pp. 321-377, 1971. [2] E . R . G . Eckert, "Analysis of film cooling and full coverage film cooling of gas turbine blades", Jour, of Eng. for Gas Turbines and Power, vol. 106, pp. 207-213, 1984. [3] I. Gartshore, M . Salcudean and P. Martys, "Fi lm Cooling Effectiveness", Proc. CSME 1990, Toronto, Ontario, 1990. [4] D. Sinitsin, "A numerical and experimental study of flow and heat transfer from a flush, inclined film cooling slot ", M.A.Sc. Thesis, University of British Columbia, Vancouver, 1988. [5] A . Brandt, "Multi-level adaptive solutions to boundary-value problems", Math. Comp., Vol.31, No.138 (1977), pp.333-390. [6] A . Brandt, "Multi-level adaptive computations in fluid dynamics", AIAA J., Vol.18, No.10 (1980), pp.1165-1172. [7] L . Fuchs and H.-S. Zhao, "Solution of three-dimensional viscous incompressible flows by a multi-grid method", Int. J. Num. Meth. Fluids, vol.4, (1984), pp.539-555. [8] L. Fuchs, "Multi-grid scheme for incompressible flows", Notes on Numerical Fluid Mechanics, vol.104, (1984), pp.38-51. [9] U . Ghia, K . N . Ghia and C.T. Shin, "High-Re solution for incompressible flow us-ing the Navier-Stokes Equation and a multigrid method", J. Comp. Phys., vol.48, (1982), pp. 387-411. [10] S.P. Vankar, "Block-implicit multigrid solution of Navier-Stokes equations in prim-itive variables," Journal of Computational Physics, vol. 65, pp.138-158, 1986. [11] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980. [12] S. Sivaloganathan and G.J . Shaw, " A multi-grid method for recirculating flows ", Int. J. Numer. Methods. Fluids, Vol.8, (1988), pp.417-440. [13] Kotake, S. and Kodama, T., "Coarse-fine mesh method for heat and mass transfer in locally complex flows ", Int. J. Num. Meth. Eng., vol. 25, 1988, pp. 347-355. 89 Bibliography 90 [14] D. Bai and A . Brandt, "Local mesh refinement multilevel techniques", SIAM J. Sci. Stat. Comput., Vol.18, No.2 (1987), pp.109-134. [15] L. Fuchs, " A local mesh-refinement technique for incompressible flows", Computer & Fluids, vol.14, N o . l , (1986), pp.69-81. [16] W. Hackbusch, "Local defect correction method and domain decomposition tech-niques", in Defect Correction Methods: Theory and Applications, K . Bohmer and H.J . Stetter, eds., Springer-Verlag, Wien, 1984, pp.89-113. [17] W. Hackbusch and U . Trottenberg, eds., Multigrid Methods, Lecture Notes in Math-ematics 960, Springer-Verlag, New York, 1982. [18] M . J . Berger, "Data structures for adaptive mesh refinement," Adaptive Computa-tional Methods for Partial Differential Equations, I. Babuska, J . Chandra and J .E. Flaherty, eds., S I A M , Philadelphia. [19] M . J . Berger, and J . Oliger, "Adaptive mesh refinement for hyperbolic partial differ-ential equations," J. Comp. Phys.,53 (1984), pp.484-512. [20] S. McCormick, "Fast adaptive composite grid methods", in Defect Correction Meth-ods: Theory and Applications, K . Bohmer and H.J . Stetter, eds., Springer-Verlag, Wien, 1984, pp.115-121. [21] J . M . Zhou and M . Salcudean, "A multi-grid local mesh refinement method for com-putation of recirculating flows", Proc. 11th Annual Meeting of Canadian Applied Mathematics Society, Halifax, Nova Scotia, 1990. [22] W. Rodi, Turbulence Models and Their Application in Hydraulics- a state of art review, International Association for Hydraulics Research, Delft, the Netherlands, 1984. [23] B . E . Launder and D.B. Spalding, "The numerical computation of turbulence flow", Comp. Meths. Appl. Mech. Engng., vol. 3, pp. 263-289, 1974. [24] R.S. Amano, "A study of turbulence flow downstream of an abrupt pipe expansion", AIAA Journal, Vol.21, No.10, pp.1400-1405, 1983. [25] C .C. Chieng and B . E . Launder, "On the calculation of turbulent heat transport down-stream from an abrupt pipe expansion", Numer. Heat Transfer, vol.5, pp. 493-496, 1982. [26] N . Djilali, I. Gartshore and M . Salcudean, "Calculation of convective heat transfer in recirculating turbulent flow using various near-wall turbulence models", Numer. Heat Transfer, vol. 16, pp. 189-212, 1989. Bibliography 91 N . Djilali, " A n investigation of two-dimensional flow separation with reattachment", Ph.D. Thesis, University of British Columbia, Vancouver, 1987. Y . Nagano and M . Tagawa, "An improved k — emodel for boundary layer flows", Journal Fluids Engineering, Vol.113, pp.33-39, 1990. W . C . Reynolds, "Commutation of turbulence flows", Ann. Rev. of Fluids Mech, vol. 8, pp. 183-208, 1976. P. Bradshaw, T. Cebeci and J .H. Whitelaw, Engineering Calculation Methods for Turbulence, Academic Press, 1981. A . Brandt and N . Dinar, "Multigrid solutions to elliptic flow problems", in Numerical Methods for PDEs, S.V. Parter, eds., Academic Press, New York, 1979, pp. 53-147. G.D. Raithby, "Skew upwind differencing schemes for problems involving fluid flow", Comp. Meths. Appl. Mech. Engng, vol.9, pp. 153-164, 1976. D. Sidilkover, "Numerical solution to steady-state problems with discontinuities", Ph.D. Thesis, the Weizmann Institute of Science, Rehovot, Israel, 1989. A . Brandt, "Multigrid Techniques: 1984 Guide, with applications to fluid dynamics", Weizmann Institute of Science, Rehovot, Israel, 1984. G.J . Shaw and S. Sivaloganathan, "On the smoothing properties of the S I M P L E pressure-correction algorithm ", Int. J. Numer. Methods. Fluids, Vol.7, (1988), pp.441-461. A . D . Gosman, " T E A C H - T , a computational program for the calculation of two-dimensional turbulent recirculating flows", Mech. Eng. Dept., Imperial College, Lon-don, Endland. W. Briggs,4 Multigrid Tutorial, S I A M , 1977. W. Briggs and S. McCormick, "Introduction", in Multigrid Methods, S. McCormick, ed., S I A M , 1987, pp.1-30. J .W. Yokota, "Multigrid calculation of 3-d turbulent viscous flows", Proc. First Canadian Symposium of Aerodynamics, Ottawa, Ontario, pp. 29.1-29.16, 1989. S. McCormick and J . Thomas,"The fast adaptive composite grid (FAC) method for Elliptic Equations", Math. Comp., Vol.46, pp.439-456. 1986, S. McCormick, Multilevel Adaptive Methods for Partial Differential Equations, S I A M , 1989.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A multi-grid method for computation of film cooling
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A multi-grid method for computation of film cooling Zhou, Jian Ming 1990
pdf
Page Metadata
Item Metadata
Title | A multi-grid method for computation of film cooling |
Creator |
Zhou, Jian Ming |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | This thesis presents a multi-grid scheme applied to the solution of transport equations in turbulent flow associated with heat transfer. The multi-grid scheme is then applied to flow which occurs in the film cooling of turbine blades. The governing equations are discretized on a staggered grid with the hybrid differencing scheme. The momentum and continuity equations are solved by a nonlinear full multi-grid scheme with the SIMPLE algorithm as a relaxation smoother. The turbulence k — Є equations and the thermal energy equation are solved on each grid without multi-grid correction. Observation shows that the multi-grid scheme has a faster convergence rate in solving the Navier-Stokes equations and that the rate is not sensitive to the number of mesh points or the Reynolds number. A significant acceleration of convergence is also produced for the k — Є and the thermal energy equations, even though the multi-grid correction is not applied to these equations. The multi-grid method provides a stable and efficient means for local mesh refinement with only little additional computational and.memory costs. Driven cavity flows at high Reynolds numbers are computed on a number of fine meshes for both the multi-grid scheme and the local mesh-refinement scheme. Two-dimensional film cooling flow is studied using multi-grid processing and significant improvements in the results are obtained. The non-uniformity of the flow at the slot exit and its influence on the film cooling are investigated with the fine grid resolution. A near-wall turbulence model is used. Film cooling results are presented for slot injection with different mass flow ratios. |
Subject |
Turbines -- Blades -- Cooling Differential equations, Partial -- Numerical solutions Numerical grid generation (Numerical analysis) |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-21 |
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. |
IsShownAt | 10.14288/1.0080376 |
URI | http://hdl.handle.net/2429/29414 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1990_A6_7 Z45_5.pdf [ 5.02MB ]
- Metadata
- JSON: 831-1.0080376.json
- JSON-LD: 831-1.0080376-ld.json
- RDF/XML (Pretty): 831-1.0080376-rdf.xml
- RDF/JSON: 831-1.0080376-rdf.json
- Turtle: 831-1.0080376-turtle.txt
- N-Triples: 831-1.0080376-rdf-ntriples.txt
- Original Record: 831-1.0080376-source.json
- Full Text
- 831-1.0080376-fulltext.txt
- Citation
- 831-1.0080376.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-0080376/manifest