Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical prediction of film cooling of turbine blades using multiblock curvilinear grids He, Pingfan 1995

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

Item Metadata

Download

Media
831-ubc_1995-05974x.pdf [ 3.3MB ]
Metadata
JSON: 831-1.0080014.json
JSON-LD: 831-1.0080014-ld.json
RDF/XML (Pretty): 831-1.0080014-rdf.xml
RDF/JSON: 831-1.0080014-rdf.json
Turtle: 831-1.0080014-turtle.txt
N-Triples: 831-1.0080014-rdf-ntriples.txt
Original Record: 831-1.0080014-source.json
Full Text
831-1.0080014-fulltext.txt
Citation
831-1.0080014.ris

Full Text

NUMERICAL PREDICTION OF FILM COOLING OF TURBINEBLADE.S USING MULTIBLOCK CURVILINEAR GRIDSByPingfan HeB.Sc., ZhengZhou University, M.Sc., IAPCM, ChinaA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFPH.DinTHE FACULTY OF GRADUATE STUDIESINSTITUTE OF APPLIED MATHEMATICSDEPARTMENT OF MECHANICAL ENGINEERINGDEPARTMENT OF MATHEMATICSWe accept this thesis as conformingto the required standardTHE ITY OF BRITISH COLUMBIA1995© Pingfan He, 1995In presenting this thesis in partial fulfilment of the requirements for an advanced degree atthe University of British Columbia, I agree that the Library shall make it freely availablefor reference and study. I further agree that permission for extensive copying of thisthesis for scholarly purposes may be granted by the head of my department or by hisor her representatives. It is understood that copying or publication of this thesis forfinancial gain shall not be allowed without my written permission.Institute of Applied MathematicsDepartment of MathematicsDepartment of Mechanical EngineeringThe University of British Columbia2075 Wesbrook PlaceVancouver, CanadaV6T 1W5Date: 19 9AbstractIn the present thesis, a computational capability is developed for the prediction of filmcooling of turbine blades. This includes the development of two comprehensive numericalcodes and the associated methods. To address the difficulties associated with complexconfigurations of turbine blades and convergence problems, several numerical techniquesare used, including curvilinear coordinate-based calculations, multigrid acceleration, domain segmentation and grid generation. Investigation and validation of these methodsare carried out and novel techniques are proposed to achieve an efficient numerical solver.Two computer codes have been developed. One is a 3D curvilinear coordinate-basedCFD code, called CMGFD, which can be used to calculate laminar/turbulent flows inarbitrary geometries using non-structured curvilinear grids. The methods developed herehave been implemented into the codes. To support the application of the CMGFD code,a multigrid elliptic grid generation code, MBEGG, was developed which can be usedto generate multi-block curvilinear grids for the CMGFD code. A multigrid methodis used to solve the three elliptic grid generation equations thus providing an efficientgrid generator. The developed computational codes are applied to study film cooling ofan experimental turbine blade model. The computational domain follows the physicalgeometry which includes a curved blade surface and a number of injection holes. A blockstructured curvilinear grid is generated by the MBEGG code which exactly represents theinclined, round film-holes and the curved blade surface. The computational results showthat the developed numerical tool has the potential to accurately model the complexcooling process in actual blade geometries.:11AbstractTable of ContentsiiList of TablesList of FiguresNomenclatureAcknowledgment2 A Computational Method Using Curvilinear Grids2.1 Introduction2.2 Numerical Formulations2.2.1 Physical Geometric Quantities2.2.2 Calculation of the Geometric Quantities2.2.3 Discretizationviiviiixivxvii113589141 Introduction1.1 Background1.2 Motivation for the Present Study1.3 Previous Numerical Studies1.4 Remarks on Film Cooling Studies1.5 Objective and Scope1.6 Contributions of the Thesis1818222325261112.2.4 Discretization of the Nonorthogonal Fluxes 302.2.5 Formulations of the Discretization Coefficients 332.3 Calculation of The Flow Field 342.3.1 Variable Arrangement and The Choice of Velocity Unknowns . 342.3.2 Discretization of Momentum Equations 362.3.3 Treatment of Pressure Gradient and Associated Boundary 402.3.4 Discretization of Continuity Equations 412.4 Solution Procedure 422.4.1 Coupling of the Velocity and Pressure 422.4.2 Symmetric Coupled Gauss-Seidel Iteration Method 432.5 Validation of the Method 442.5.1 Cavity flow with inclined side walls 442.5.2 Laminar flow through a tube with a constriction 462.5.3 Developing laminar flow in a pipe 482.5.4 Flow in a pipe with a smooth 90 degree bend 492.6 Closure 543 Computation of Turbulent Flows Using Curvilinear Grids 553.1 Introduction 553.2 Governing Equations 573.3 Discretization of the k — e Equations . . 583.4 Wall Boundaries 623.5 Application to Turbulent Flows 643.5.1 Developing Turbulent Pipe Flow 643.5.2 Developing Turbulent Flow in a 90° Curved Square Duct 693.5.3 Developing Turbulent Flow in a 90° Curved Pipe 71iv3.6 Closure 794 Multigrid Acceleration 804.1 Introduction 804.2 The Multigrid Algorithm 844.3 Consistent Formulation of Governing Equations 864.4 Consistent Calculation of the Grid Geometric Quantities 904.5 Multigrid Acceleration of Turbulent Flows 914.6 The Treatment of Wall Boundary on the Coarse Grids 924.7 The Treatment of k — e Quantities and Source Terms 954.8 Restriction and Prolongation 974.9 Test Calculations 1004.9.1 Laminar Flow in a Skewed Cubic Cavity 1014.9.2 Laminar Flow in a Duct with Obstruction 1054.9.3 Laminar Flow in a Strongly Curved 90° Pipe 1094.9.4 Developing Turbulent Flow in A Strongly Curved 90° Pipe . . 1134.10Closure 1155 Domain Segmentation and the CMGFD Code 1165.1 Introduction 1165.2 Domain Segmentation Technique 1175.3 Communication Between Blocks 1195.4 The CMGFD Code 1245.5 Test Computations 1285.5.1 Computations of a L-Shaped Duct . . 1285.5.2 Calculation of A Hydrocyclone . . . . 1355.5.3 Laminar Flows Over A Turbine Blade. 137v5.6 Closure. 1436 Numerical Grid Generation6.1 Introduction6.2 Elliptic Grid Generation6.3 Control Functions6.4 Multigrid Solver6.5 Multi-Block Grid Generation6.6 The Code Structure6.7 Application of the Code6.8 Closure1441441461471501521541541557 Prediction of Film Cooling Over7.1 Problem Description7.2 Computational Grid7.3 Solution Model7.4 Boundary Conditions7.5 Computational Procedure .7.6 Results and Discussion7.7 Closure8 Conclusions And Recommendations8.1 Conclusions8.2 Recommendations for Future WorkBibliography 191A Turbine Blade 159160162162165167169184185185189viList of Tables2.1 Number of iterations necessary for convergence 462.2 Comparison with experimental results and numbers of iterations for calculation 484.1 Convergence characteristics of all calculations (the superscript * indicatesthat the values are estimated) 105viiList of Figures1.1 Illustration of a discrete-hole cooling arrangement for a turbine blade. . . 22.1 Illustration of the physical geometric quantities for a control cell 242.2 A general control volume. Nodes to the north, south, east, west, top andbottom are represented by N, 5, E, W, T and B respectively. Lower casedenotes control volume surfaces 282.3 Illustration of discretization of cross-derivatives 312.4 Illustration of the three types of physical velocity components in curvilinear grids, OA1 and 0A2: physical contravariant velocity components,OB1 and 0B2: physical tangential velocity components, and OC1 and0C2: physical covariant velocity components 352.5 Velocity decomposition in a locally fixed coordinate system 392.6 A scalar control volume with 7 unknowns 422.7 Streamlines for flow in cavities with inclined side wall, (a) /3 = 60°, (b)/3 = 45°, (c) /3 = 30° 452.8 Flow through tube with constriction: (a) geometry and grid, (b) velocityfield, (c) streamlines 472.9 (a) Illustration of grid and velocity, (b) development of velocity profile, (c)velocity contour in the cross-pipe plane at x = 1OD, (d) comparison withanalytic solution at x = 1OD, (e) Comparison for flow in the developingregion, where = 502.10 Grid and flow geometry for a pipe with a smooth 90° bend 51viii2.11 Velocity vectors on the symmetry plane 512.12 Contours of mean velocity U/UB in the bend 522.13 Comparison with the experimental results of Eanyet et al (1982) 533.1 (a) The non-orthogonal rectangular type grid. (b) Predicted velocity fieldat the symmetry plane 663.2 Predicted velocity profiles for Re = i0, (a) developing velocity profileat the pipe center, (b) developed velocity profile at cross-section wherez/D 60 673.3 Predicted velocity profiles for Re = 3.0 x iO, (a) developing velocityprofile at the pipe center, (b) developed velocity profile at cross-sectionwhere z/D = 60 683.4 Illustration of computational grid for the curved duct 703.5 Flow field in the middle vertical plane, Re = 40, 000, (a) flow field in thewhole duct, (b) flow field in the curved section 723.6 Secondary flow field and primary flow velocity contours of U/Ub at twodifferent cross-sections 733.7 Comparisons of velocity profiles with the measurements at four differentlocations, (a) 3 = 30°, (b) /3 = 60°, (c) 3 = 75°, (d) downstream xd = 0.25D 743.8 Illustration of computational grid for the curved pipe 753.9 Flow field in the middle vertical plane, Re = 40, 000 763.10 Secondary flow field and primary flow streamline velocity contours at twodifferent cross-sections 773.11 Comparisons of velocity profiles for different cross-sections, (a) upstream0.5D, (b) 3 = 30°, (c) /3 = 60°, (d) 3 = 75° 784.1 Grids in cylindrical coordinates: (a) 16 x 16 grid; (b) 2 x 2 grid 87ix4.2 Curvilinear grids: (a) 16 x 16 grid; (b) coarsest 2 x 2 grid; (c) 2 x 2 gridrestricted from the fine grid 16 x 16 884.3 Illustration of using inconsistent defect equations at interior scalar nodes934.4 Illustration of locations of UE3 — residuals on fine and coarse grids forrestriction 994.5 Physical geometry of the skewed cube 1014.6 Convergence history of MG calculations using various grids 1024.7 Comparison of convergence behaviors of multigrid and single grid calculations 1034.8 Illustration of the grid 65 x 17 x 17 used in the study 1074.9 Illustration of the coarsest grid 5 x 3 x 3 in a xy— plane 1074.10 Convergence history for computation of laminar flow in a duct using finestgrid 129 x 33 x 33 1084.11 Finest grid in the MG calculation of curved pipe flow, 33 x 33 x 97. . 1104.12 Coarsest grid 3 x 2 x 7: (a) grid restricted from the finest grid; (b) gridobtained from the coarsest grid nodes 1114.13 Convergence histories for computations using different grid levels 1124.14 Convergence history for the computation of a developing turbulent flow ina strongly curved pipe 1145.1 Illustration of multi-block curvilinear grids with discontinuous grid lineslopes across the interface 1215.2 Illustration of the semi-matched multi-block curvilinear grids 1255.3 Illustration of the adopted two-block non-orthogonal grid at the symmetryplane for the b-shaped duct 129x5.4 The two-block Cartesian grid at the symmetry plane for the b-shaped duct.1305.5 Predicted velocity field at the symmetry plane by the non-orthogonal grid. 1315.6 Predicted velocity field at the symmetry plane by the orthogonal grid. . . 1325.7 Comparison of the u-velocity profiles at the interface shown in Figure 5.3. 1335.8 Comparison of maximum residuals by multigrid and single-grid calculations for the non-orthogonal grid 1345.9 Illustration of a conventional hydrocyclone 1365.10 (a) 3-block mesh used in the present study, (b) predicted velocity field. . 1375.11 Computational domain for a turbine blade model 1385.12 A 4-block curvilinear grid generated for the domain shown in Figure 5.11. 1395.13 The velocity field at the xz-plane through the middle of an injection hole. 1405.14 The velocity at a surface distanced from the turbine wall by two cells . 1415.15 The streamline-velocity contour at the cross-plane of an injection hole; a)one cell size distance from the turbine surface; b) two cell size distance fromthe turbine surface; c) three cell size distance from the turbine surface; d)four cell size distance from the turbine surface 1426.1 Illustration of the face numbers in a 3D block 1486.2 The mesh generated for a finer grid with 32 x 32 cells 1566.3 Convergence behavior, (a) multigrid calculation with the finest grid of64 x 64. (b) single-grid calculation with 32 x 32 grid 1566.4 A two-dimensional domain to be considered 1577.1 Illustration of the turbine blade model used in the experimental study ofSalcudean et al (1994) 161xi7.2 Illustration of the computational domain, including a main flow region,• two half injection-hole ducts in the first row and one injection hole ductin the second row 1617.3 The grid generated for the present calculation containing 101 x 17 x 37 gridnodes in the main flow region and 9 x 9 x 17 grid nodes in one injectionhole 1637.4 Computational grid at the curved blade surface and the injection hole ducts 1647.5 M, = 0.52; Solutions at the streamwise plane (xz-plane) through thecenterline of the 150 coolant orifice: (a) local velocity field around thecurved blade surface, (b) pressure distribution in the whole xz-plane, (c)local velocity field (near the orifice exit) and (d) local pressure, (e) localtemperature, (f) local turbulence kinetic energy 1707.6 M = 0.97; Solutions at the streamwise plane through the centerline ofthe 150 coolant orifice: (a) local velocity field around the curved bladesurface, (b) pressure distribution in the whole xz-plane, (c) local velocityfield (near the orifice exit) and (d) local pressure, (e) local temperature,(f) local turbulence kinetic energy 1717.7 M = 0.52; Solutions at the streamwise plane through the centerline ofthe 440 coolant orifice: (a) local velocity field around the curved bladesurface, (b) pressure distribution in the whole xz-plane, (c) local velocityfield (near the orifice exit) and (d) local pressure, (e) local temperature,(f) local turbulence kinetic energy 173xii7.8 M = 0.97; Solutions at the streamwise plane through the centerline ofthe 440 coolant orifice: (a) local velocity field around the curved bladesurface, (b) pressure distribution in the whole xz-plane, (c) local velocityfield (near the orifice exit) and (d) local pressure, (e) local temperature,(f) local turbulence kinetic energy 1747.9 M = 0.52; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a 15°, (b) a = 22°, (c) a = 28° and (d) a = 36°.1767.10 M = 0.97; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a = 15°, (b) a = 22°, (c) a = 28° and (d) a = 36°.1777.11 M = 0.52; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a = 440, (b) a = 55°, (c) a = 60° and (d) a = 80°.1787.12 M = 0.97; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a = 440, (b) a 55°, (c) a = 60° and (d) a = 80°.1797.13 Contours of predicted FCE for M = 0.52 at the blade surface 1807.14 Contours of measured FCE for M = 0.52 at the blade surface 1817.15 M = 0.52; Comparison of the spanwise averaged film cooling effectivenesswith the measured result 1827.16 Contours of predicted FCE for M = 0.97 at the blade surface 1827.17 Contours of measured FCE for M = 0.97 at the blade surface 1837.18 M = 0.97; Comparison of the spanwise averaged film cooling effectivenesswith the measured result 183xliiNomenclatureC,, C1, C2 The Ic — 6 turbulence model constants.F2 Volume flow rates across cell faces.G Generation rate of turbulent kinetic energy.J Generalized flow flux.M Mass flow ratio, M = 000UooM15 Mass flow ratio of the first row coolant orifices.M44 Mass flow ratio of the second coolant orifices.Overall mass flow ratio averaged from all coolant holes.F2, i=1,2,3 Control functions in grid generation equations.S’ Surface area vectors.T Temperature.U2 Contravariant velocities.tangential velocity components.TwuT Friction velocity, UT = —eV Volume of a control cells.X Streamwise distance measured from stagnation line.Z Spanwise distance measured from one edge of a first row hole.a Covariant vectors.a1 Contravariant vectors.a23 Coefficients of the grid generation equations.a& Coefficients of the discretized governing equations.xive Tangential vectors.e’ Normalized surface vectors.d Injection orifice diameter.g22 Surface area metric tensor.k Turbulent kinetic energy.p Mean static pressure.u Velocity vector.xi Cartetian coordinates.x, y, z Cartetian coordinates.Dimensionless distance from the solid wall, y eYUTa2 Non-orthogonal angles.Dissipation rate of turbulent kinetic energy.TooTawFilm cooling effectiveness,=Von Karman constant, c = 0.41.Dynamic viscosity.Curvilinear coordinates.Fluid Density.Turbulence model constants.Wall shear stress.Non-dimensional temperature,= —General dependent variable.1’ General diffusivity coefficient.Subscriptaw Adiabatic wall.c Coolant injection.xvco Main stream.e, w, ... Center of control volume surfaces.E, W, ... Control volume center.15 Referring to row 1.44 Referring to row 2.Superscriptc Coarser grid.f Finer grid.oriented variables.xviAcknowledgmentI wish to express a deep sense of gratitude to my supervisors Dr. M. Salcudean andI. S. Gartshore for their encouragement, guidance and support throughout the course ofthis research. Working with them has been an enjoyable experience which I will cherishthroughout my life.I would like to thank my co-supervisor Professor B. Seymour in the Institute ofApplied Mathematics, for the numerous pieces of advice he has give me. Thanks are alsodue to Professor J. Varah and Professor M. Quick.I would also like to thank all my colleagues of the CFD group in the Department ofMechanical Engineering, especially Dr. Z. Abdullah and Dr. P. Nowak. Their suggestionsand discussions have been in valuable. Thanks to Mr. M. Findlay for his proofreadingof this thesis.Finally, I wish to thank my wife, Hui, for her patience and support.xviiChapter 1Introduction1.1 BackgroundOne of the most effective ways to improve gas turbine engine efficiency is to increasethe working fluid temperature. Consequently, new turbine designs are continually being pushed to operate at higher temperatures. The modern gas turbine engines operateunder working temperatures of 1400— 1500°C which are beyond the allowable metal temperatures. Therefore, to maintain reasonable engine life and safety, critical components,such as the turbine blades, have to be protected.Film cooling is a widely used technique to protect turbine blades from the high temperature of the surrounding gas stream. This technique involves injection of coolant atone or more discrete locations along a surface exposed to a high temperature environmentin order to protect the surface in the downstream region. A discrete-hole cooling arrangement for a turbine blade is illustrated in Figure 1.1. In this blade there are a numberof discrete, small injection holes on the surface from which coolant is injected into theboundary layer on the blade surface to form a thin layer of coolant which protects theblade surface from the high temperature environment. The cooling process occurs in acomplex 3D configuration which includes not only a curved airfoil but also a number ofsmall injection holes. The aim of film cooling design is to minimize the coolant injectedwhile maximizing the thermal protection for the blades. A large number of experimentaland computational studies have been devoted to the understanding of the film cooling1Chapter 1. Introduction 2Injection tubesHot Main FlowFigure Li: Illustration of a discrete-hole cooling arrangement for a turbine blade.Chapter 1. Introduction 3process in order to achieve higher film cooling effectiveness. A brief literature review isprovided in the following paragraphs.1.2 Motivation for the Present StudyFilm cooling has been a subject of research for many years. There are a large numberof research papers on film cooling of turbine blades in the open literature. Goldstein(1971) gave a detailed review of film cooling research and summarized the experimentalresults on film cooling effectiveness with normal, tangential or nearly tangential injectionup to 1971. There are two major techniques in film cooling studies; one is experimentalinvestigation and the other is numerical prediction. Most early investigations of filmcooling were carried out using experimental techniques. Current engine design is highlyempirical and relies on a large experimental database. Recent works are more concernedwith discrete-hole film cooling on a flat plate or a curved blade surface. There havebeen numerous experimental investigations of discrete hole cooling. However, a reviewof experimental studies is not the subject of the present thesis. Instead, only a selectionof some experimental works associated with the physical model of the present numericalstudy are mentioned.In the film cooling of a turbine blade, the stagnation region near the leading edge is aparticularly critical region as it suffers the most intense exposure to the hot free streamand therefore requires particular attention. A number of studies have been carried out forfilm cooling at the leading edge. Ou & Han (1991, 1992) and Mehendale & Han (1992)carried out a systematic investigation of the effects of high mainstream turbulence andthe effects of injection orifice geometry on the leading edge film cooling. They measuredthe film cooling effectiveness over a blunt body with a semi-cylindrical leading edge andChapter 1. Introduction 4a fiat afterbody. Here, the film cooling effectiveness is defined as follows,TjoTwT—T(1.1)where Taw is the adiabatic wall temperature, T is the mainstream temperature andT is the coolant temperature. Two rows of orifices at +150 and +400 were used withspanwise blowing and with spacing-to-diameter ratios of 4 and 3. The effectiveness wasfound to decrease with increasing mainstream turbulence; however, this effect decreasedwith increasing blowing rate. Their study also indicates that the spacing-to-diameterratios have a significant influence on the film cooling effectiveness.Recently, a systematic experimental investigation of film cooling effectiveness nearthe leading edge of a turbine blade was carried out by Gartshore et al (1993) and Salcudean et al (1994) at the University of British Columbia. This study was carried out tofulfill a research collaboration between Pratt & Whitney Canada and the University ofBritish Columbia. A large turbine blade model was used in this experimental study. Thismodel has a semi-circular leading edge with four rows of laterally-inclined film coolingorifices positioned symmetrically about the stagnation line. Both air and CO2 were usedas coolant. A heat-mass transfer analogy was used to measure the cooling effectiveness.Their experimental results indicate that the film cooling effectiveness changes substantially with respect to the overall mass flow ratio. Here, the mass flow ratio is defined aseU/ where Pc and p are the coolant density and the mainstream flow densityrespectively. The best effectiveness was obtained in a range of mass flow ratios 0.5 — 0.6.A continuous decrease in the film cooling effectiveness was found when the mass flowratio was further increased.The above experimental studies indicate that the film cooling effectiveness is influenced by a number of factors, such as, spacing-to-diameter ratios, mass flow ratios andChapter 1. Introduction 5free stream turbulence. Therefore, in order to achieve an optimal design of the film cooling system, a parametric analysis has to be carried out to determine the film coolingperformance with various blade geometries and flow parameters. For this purpose, experimental measurements are too expensive and time-consuming to obtain the detailedflow and heat transfer data. In addition, measurement is also limited by the availabilityand accuracy of the measuring instruments for many situations. An alternative methodin film cooling research is numerical prediction.Recent advances in high speed digital computers and computational fluid dynamicsmake it possible to model the flow field and heat transfer phenomena in film cooling ofturbine blades. Numerical simulation can provide detailed information about flow fieldsand constitutes a fast and inexpensive tool for understanding the film cooling process.The present study is motivated by the above need and is devoted to the numericalsimulation of film cooling of turbine blades. It is our aim to develop a computationalcapability for the prediction of film cooling of a turbine blade with a realistic geometry.1.3 Previous Numerical StudiesNumerical simulation of the film cooling of turbine blades is not new. There are a numberof numerical studies in the literature and the first research paper can be traced back toKacker et at (1969). However, most of the early numerical studies were confined to twodimensional flows with wall jet configurations. Among the three-dimensional studies,Bergeles et at (1978) used a partially parabolic numerical scheme to predict the meanvelocity and temperature of laminar flow for a single row of inclined holes and for asurface with multiple rows of holes. Later, they used the same method in conjunctionwith an anisotropic turbulence model to calculate a single jet injected into a crossflow at900 and 300. Good agreement with measured results for an inclined injection hole wasChapter 1. Introduction 6obtained over most of the domain with the exception of the region in the vicinity of the jetorifice where the agreement was poor especially for high mass flow ratios. This is partlycaused by the assumption of a uniform velocity profile at the jet exit plane. Demuren etal (1986) carried out a systematic study for a single row of holes in a flat plate inclinedat different angles using the locally-elliptic procedure. As in the treatment of Bergeleset al (1978), uniform velocities at the jet exit plane were assumed in order to excludethe coolant orifice ducts from the computational domain. The predicted film coolingeffectiveness was in poor agreement with the measured data in some cases, especially forhigh mass flow ratios. Sathyamurthy and Patankar (1990) investigated discrete-hole filmcooling over a flat plate with lateral injection. A laterally periodic parabolic procedurewas used in their calculation. Their results indicate that lateral injection can operate athigh blowing rates and can achieve better film coverage than streamwise injection.Computations of film cooling by solving elliptic Navier-Stokes equations have alsobeen reported. White (1981) solved the fully elliptic transport equations for a singleinjection hole on a flat plate. The novel feature of the method is its use of separatecomputational domains in the injection pipe and cross-stream regions. Comparison ofjet exit profiles with measurements showed good agreement. Jubran (1989) used a commercial CFD code, PHOENICS, to predict the film cooling effectiveness and the velocityfield for two rows of holes inclined in the streamwise and spanwise directions over a flatplate. More recently, Zhou et al (1993) used a multigrid and a segmentation techniquefor the prediction of film cooling for the case of discrete holes with lateral injection. Withthis technique, they were able to use a very fine grid and the solution also extended intothe injection hole regions for the vertical injection hole case. The computational resultsshowed good agreement with the experiments.More recently, computations have been reported for more complex geometries usingbody-fitted coordinates to approach real blade configurations. This includes the use ofChapter 1. Introduction 7some commercial CFD codes. Leylek and Zerkle (1993) conducted computations for astreamwise inclined orifice on a fiat plate. Their solutions include flow in the plenum,coolant tube and cross-stream regions. A curvilinear grid was used to represent thecircular injection tube which formed an angle with the blade wall in the streamwisedirection. They reported important information on the flow in the coolant hole whichcontains counter-rotating vortices and local jetting effects. The results showed that fora hole length of about 4 times the hole diameter, there is a large separation region inthe film cooling tube which leads to redistribution of the coolant flow at the hole exit.However, the influence of wall curvature was not considered since the computations didnot include the curved blade surface.Garg and Gaugler (1993, 1994) used general body-fitted coordinates to model thecomplex flow domain around a real turbine airfoil. The fully three-dimensional NavierStokes equations were solved together with closure proposed by the Baldwin and Lomax(1978) algebraic turbulence model. The computational domain included the whole curvedblade surface but not the coolant ducts. Turbulent (1/7 power law) profiles were specifiedat the duct exits, in conformity with observations of Leylek and Zerkle (1993). Experimentally measured wall temperatures were specified as boundary conditions at the airfoilsurface and wall heat fluxes were calculated. Fair agreement with the experimental results was reported. Their study provided important guidelines for numerical simulationsof film cooling for real turbine blades.Vogel (1994) carried out calculations of turbine flows with film cooling using multiblock body-fitted coordinates. Two cases are reported in Vogel’s paper, the first dealingwith two-dimensional film cooling with a slot configuration. A curvilinear grid was usedto represent the curved leading edge. The second calculation deals with film coolingfrom a single orifice on the suction side of the blade. Since the calculation deals onlywith a single streamwise inclined injection orifice, the author was able to use a symmetryChapter 1. Introduction 8condition through the middle of the injection orifice. The computational domain includedthe suction side of a blade surface and half of the streamwise injection orifice. The paperfocussed on the detailed flow field and did not describe the temperature field.1.4 Remarks on Film Cooling StudiesActual measurements provide the most reliable information about a physical process.However, it should be mentioned that there are serious difficulties with measurementsin many situations and the measuring instruments are not completely free from errors.In addition, experimental investigations are usually expensive and time-consuming inobtaining the complete flow data. Numerical simulation can furnish detailed informationabout fluid flow and heat transfer phenomena and provide a fast and inexpensive tool forunderstanding the film cooling process. However, numerical simulation has not reacheda stage where it can be directly applied in engine design. One of the major problems isthat there is lack of appropriate methods to model the cooling process for the complexgeometry of a turbine blade. Therefore the majority of the previous calculations werelimited to simplified geometries.The complex configuration of an actual turbine blade presents a number of difficultiesfor numerical simulation. First, the flow region includes a curved blade surface and anumber of circular injection holes for which an ordinary coordinate system, such asCartesian or cylindrical coordinates can not be applied. Secondly, the flow occurs indifferent flow regions: a main flow region around the blade surface and a number of subregions in the small coolant injection holes for which ordinary one-segment methods areineffective. Thirdly, the diameters of the injection holes are very small compared withthe characteristic dimension of the main flow region. When a fairly dense grid is used inthe injection holes, the number of grid nodes needed overall will be too large to be usedChapter 1. Introduction 9with the present computer capacities without using excessively nonuniform grids. Theuse of a very non-uniform and dense grid near the cooling orifices can cause additionalconvergence difficulty. Therefore, a solution method should be developed to address boththe convergence and geometrical difficulties.The use of simplified geometries, such as a fiat plate with discrete wall jets, can provideuseful flow information about the complex cooling process. However, the curvature of areal blade is not accounted for in such a model; consequently the film cooling might besignificantly miscalculated. Due to the shape of turbine passages, the flow through themexperiences very strong streamline curvature. This curvature can significantly influencethe cooling process. Furthermore, as reported by Leylek and Zerkie (1993), the velocityprofile at the orifice exit is very complex and is neither uniform nor in alignment with thecoolant tube. Our study also showed that the velocity, temperature and k — c profiles ata coolant-hole exit plane strongly depend on the hole location and the mass flow ratio.Thus, an assumption of constant or power-law solution profiles may be very inaccurate.Therefore, in order to obtain accurate flow and heat transfer data, the simulation shouldinclude both the curved blade surface and the coolant ducts.1.5 Objective and ScopeThe present work has two main objectives. (1). One objective is to address the needdiscussed above to develop an efficient numerical tool suitable for modeling the filmcooling process over realistic physical geometries. (2). The other objective is to use thisnumerical tool to study film cooling of the TJBC experimental turbine blade for whichconsiderable measured data has been obtained. This provides not only a validation ofthe experimental results but also a detailed flow and temperature field.Chapter 1. Introduction 10The first objective involves both code and method developments. To address the difficulties associated with complex geometries and convergence problems, several numericaltechniques are used including grid generation, curvilinear coordinate-based calculation,multigrid method and domain segmentation. Two computer codes are developed. Oneis a 3D curvilinear coordinate-based CFD code, called CMGFD, which can be used tocalculate laminar/turbulent flows in arbitrary geometries using non-structured curvilinear grids. The other is a multigrid elliptic grid generation code, called MBEGG, whichcan be used to generate multi-block curvilinear grids for the CMGFD code.The CMOFD code is developed based on an existing CFD code, MGFD, developed byNowak (1992) in our research group. The MGFD code is a comprehensive 3D numericalcode with multigrid and domain segmentation techniques. However, its application isstill limited to rectangular type geometries since it is based on the Cartesian coordinatesystem. To address the present need for the prediction of film cooling of turbine blades,the CMGFD code is developed using block-structured curvilinear grids. All the techniques used in the MGFD code, such as the multigrid method and domain segmentationfeature, are generalized to the CMGPD code. It should be mentioned that due to thecomprehensive nature of the MGFD code, substantial work was required to implementthe curvilinear grid-based calculation. In addition to changing the basic grid system,other significant changes, such as inclusion of the multigrid method, are also made.The CMGFD code involves several numerical techniques. These techniques are implemented in the code in four steps. In each step, considerable effort is involved in developingefficient methods suitable for the code. A detailed investigation of the proposed methodand validation of the code is carried out. The development of these methods comprisesone of the important tasks in the present study and is described in the following paragraphs.In the first step, the calculation of laminar flow in complex geometries using generalChapter 1. Introduction 11curvilinear grids is investigated. This step is necessary to develop an efficient, curvilinear coordinate-based finite-volume method while avoiding the additional complexityof turbulence modeling and domain segmentation treatments. Finite-volume methodsin conjunction with curvilinear grids have been become popular in recent years for flowsimulation in complex geometries. However, problems still exist in terms of the accuracy and efficiency of such computations. In the present study, a finite volume methodusing general curvilinear grids is proposed to overcome some of these difficulties. Thecoordinate invariant conservation equations and the physical geometric quantities of thecontrol cells are used directly to formulate the numerical scheme without reference tothe commonly-used covariant and contravariant vectors. The definitions of the physicalgeometrical quantities and the tangential velocity unknowns are introduced. A new difference scheme for the non-orthogonal diffusion and pressure terms is proposed. This schemecontributes to the main diagonal terms in the resulting coefficient matrix and allows animplicit treatment of the non-orthogonal quantities without increasing the number ofcomputational molecules. A coupled equation solver is used in place of the commonly-used pressure-correction equation associated with grid non-orthogonality. The developedmethod is implemented in the CMGFD code. Several two- and three-dimensional laminar flows are computed and compared with other numerical, experimental and analyticalresults to validate the solution method and code. The main work in this step is presentedin chapter 2.In the second step, the method developed in the first step is generalized to turbulentflows. It is well known that turbulent flows are more difficult to handle. This step isintended to increase our knowledge in dealing with turbulent flows in the environment ofgeneral curvilinear grids and achieve an efficient numerical method for turbulent flows.The mathematical model adopted for the present study is introduced. The k — c twoequation model together with the ‘wall function’ is used to simulate turbulent flows.Chapter 1. Introduction 12The discretization of the k — e equations and the formulation of the ‘wall function’ arepresented. The treatment of the source terms in the k — e equations is discussed indetail. The performance of the method is investigated through several three-dimensionalturbulent flows to validate the method as well as the CMGFD code. This is presented inchapter 3.The third step is concentrated on the development of an efficient multigrid methodfor flow simulation in general curvilinear grids. The multigrid method is an efficient iterative solution procedure which exhibits convergence rates insensitive to grid refinement.It has been widely used in the field of computational fluid dynamics. However, therehave been relatively few studies carried out using general curvilinear grids. There areadditional difficulties for multigrid acceleration in curvilinear coordinate systems whichrequire careful consideration.In the present study, a multigrid method for calculating laminar/turbulent flows usinggeneral curvilinear grids is developed. First, the multigrid calculation of laminar flowsis studied with emphasis on the influence of complex flow boundaries. It is found thatthe discrete governing equations on different grid levels can become inconsistent in somecurvilinear grids due to complex flow boundaries, thus reducing the efficiency of themultigrid algorithm. A novel treatment is proposed to solve this problem. Secondly, thedifficulties associated with the multigrid acceleration of turbulent flows are discussed andseveral techniques which help the performance of the multigrid method are introduced.Particularly, it is found that the formulation of the coarse-grid defect equation using the‘wall function’ can cause inconsistency of the governing equations between fine and coarsegrids. This problem is discussed in detail and a novel approach is proposed to allow thesuccessful implementation of the multigrid method for turbulent flows. Work in this stepis presented in chapter 4.In the last step, the method developed in the three previous steps is generalized toChapter 1. Introduction 13non-structured curvilinear grids by using a domain segmentation technique. This technique divides the domain of interest into different sub-domains. Solutions are obtainedby iteratively applying the solver described in the previous steps to each sub-domain.The implementation of the domain segmentation strategy in the method using generalcurvilinear grids is studied. Particular attention is given to the communication betweenneighbouring sub-domains. The performance of the multi-block method is investigatedthrough several computational examples.The methods described above have been implemented in the CMGFD code. With theabove developments the code can be used to deal with arbitrary complex geometries ifproper block-structured curvilinear grids are generated. However, generating appropriategrids is often a difficult task. The existing difficulty in generating appropriate gridsis one of the reasons that few numerical studies are carried out over complex threedimensional domains. In the present study an elliptic grid generation code, MBEGG,is developed to support the application of the CMGFD code. Even though the code isinitially developed to address the need for the prediction of film cooling of a turbine bladeit can be used to generate multi-block curvilinear grids for various complex flows. Anelliptic grid generation method is used in the code, which solves three nonlinear ellipticgrid generation equations. The source terms in the equations are problem-dependentfunctions and can be used to control the grid characteristics, such as orthogonality andgrid stretching. In the present study, methods to achieve certain grid quantities, suchas skewness, through the control functions are implemented in the code. The multigridmethod is an ideal solver for the elliptic grid generation equations due to the fully ellipticnature of the equations. In the present study, a multigrid method is developed to solvethese equations. A domain segmentation technique is adopted in the code to generateblock-structured curvilinear grids. Block-structured curvilinear grids are generated forflows over a cooled turbine blade using the code.Chapter 1. Introduction 14As mentioned earlier, the main objective of the present thesis is to develop an efficientnumerical tool for the prediction of fluid flows in complex geometries with emphasis oncalculation of film cooling problems. To this end, the developed computational code isapplied to study film cooling of a UBC experimental turbine blade model. This modelhas a semi-circular leading edge with four rows of film cooling orifices positioned symmetrically about the stagnation line. It uses a lateral injection and the cooling orificesare inclined spanwise 300 to the turbine blade surface. The computational domain istaken directly from the physical geometry without simplification and includes not onlythe curved blade surface but also the entrance ducts of the circular coolant holes. Thephysical domain is segmented into a number of sub-domains and separate curvilineargrids are generated for different flow regions. Grids are generated by the MBEGG codewhich represents exactly the curved blade surface as well as the circular injection orifices.ec UComputations have been carried out for a low overall mass flow ratio (M = ) ofe:nftyU0.52 and a high mass flow ratio of 0.97, where e is density, U is velocity and the subscriptc and co indicate the coolant and hot flow respectively. The computational results arecompared with the available experimental results.1.6 Contributions of the ThesisIn summary, the present thesis contributes to both the engineering application and theoryof CFD in the following areas:• A numerical method for computing fluid flow in complex three-dimensional geometries using curvilinear grids is developed. Unlike previous numerical studies,the original coordinate-invariant conservation equations and the physical geometricquantities of the control cells are used directly to formulate the numerical schemeChapter 1. Introduction 15without reference to the coordinate derivatives of transformation. The physical geometric quantities, including volumes, surface areas and surface normal directionsof control cells, are introduced to formulate the governing equations in place of thecommonly-used covariant and contravariant vectors. A calculation and interpolation procedure for the physical geometric quantities is developed. This methodallows the use of significantly non-smooth grids.• A new scheme for handling the non-orthogonal terms in the momentum and pressure equations is developed. This scheme contributes to the main diagonal termsin the resulting coefficient matrix and allows an implicit treatment of the non-orthogonal quantities without increasing the number of computational molecules.This treatment of the non-orthogonal terms enhances computational stability andconvergence rate.• The physical tangential velocity components resulting from the velocity expansionin the unit tangent vector basis are proposed as dependent variables in the momentum equations. The discrete momentum equations using the tangential velocitycomponents are obtained by an algebraic manipulation which avoids the tensor expression of Christoffel symbols. A coupled solution procedure is used in place of thecomplicated pressure-correction equation associated with grid non-orthogonality.• Numerical simulation of turbulent flows in complex geometries is studied. The k —two-equation model together with the ‘wall function’ treatment is used as the turbulence closure. Discretization of the k — e equations is presented with particularattention to the linearization of the source terms and calculation of the turbulenceenergy generation rate. Methods are introduced to enhance the computational stability and facilitate the calculation of the energy generation rate. The formulationof the ‘wall function’ in general curvilinear grids is developed.Chapter 1. Introduction 16A multigrid procedure for the computation of three-dimensional laminar/turbulentflows in general curvilinear grids is developed. It is shown that the discrete governing equations can become inconsistent between fine and coarse grids in certaincurvilinear grids, as well as when the k — 6 model with the cwall function’ is usedfor turbulent flows. Novel procedures have been proposed to solve these problems.Computational results using fine grids of up to approximately one million grid nodesshowed that the developed multigrid algorithm is very efficient and a reduction ofup to 99% in CPU time was achieved.• A domain segmentation strategy in conjunction with general curvilinear grids is developed. This approach allows the treatment of arbitrary three-dimensional geometries using non-structured curvilinear grids. Particular attention has been given forthe communication between different sub-domains. The performance of the multi-block method is investigated through several computational examples which showthat the developed method has the capability to deal with complex flow domainsusing multi-block structured curvilinear grids and allows the use of significantlydiscontinuous grid slopes at the interface.• A block-structured curvilinear coordinate-based finite-volume code, CMGFD, isdeveloped based on an existing Cartesian coordinate-based CFD code. The newcode has the capability to deal with arbitrary geometries by using separate curvilinear grids in different flow regions. There is no restriction for curvilinear gridsto be used and no limitation for the number of segments to be used in a flow region. An efficient multigrid method is implemented in the code which promises fastconvergence for both laminar and turbulent flows.Chapter 1. Introduction 17• A multi-block elliptic grid generation code, MBEGG, is developed for generatingblock-structured curvilinear grids in complex geometries. The code uses an elliptic grid generation method which solves three elliptic partial differential equations.A multigrid method is developed to solve these equations which provides an efficient grid generator. Techniques to control grid quantities, such as skewness, areintroduced. A multi-block grid generation strategy is adopted to allow the generation of block-structured curvilinear grids in arbitrary geometries. This code hasthe capability to generate block-structured curvilinear grids in arbitrary 2D or 3Ddomains.• Computations are carried out for film cooling of a UBC turbine blade model. Block-structured grids are generated by the developed MBEGG code which exactly represent the inclined, round film-holes and the curved blade surface. Computationsover the cooled turbine blade model are carried out for overall mass flow ratios of0.52 and 0.97. The detailed flow field is obtained which improves our understandingof film cooling. The structure of vortices is investigated numerically. This providesinformation on the mixing of coolant with hot flow. The predicted results arecompared with the experimental results of Salcudean et at (1994) and show goodagreement. The present computations show that the developed numerical tool hasthe potential to model the complex cooling process in actual blade geometries.Chapter 2A Computational Method Using Curvilinear GridsIn this chapter, a finite-volume method for laminar flow in complex 3D geometries using general curvilinear grids is presented. A new second-order accurate difference schemefor the cross-derivative terms is proposed to describe the non-orthogonal components,allowing parts of those terms to be treated implicitly without increasing the numberof computational molecules. The coordinate-invariant conservation equations and thephysical geometric quantities of control cells are used directly to formulate the numerical scheme without reference to the coordinate derivatives of transformation. A coupledequation solver is used in place of the complicated pressure-correction equation associated with grid non-orthogonality. Several two- and three-dimensional laminar flows arecomputed and compared with other available numerical, experimental and analytical results to validate the solution method. Part of the work presented in this chapter hasbeen reported previously in He and Salcudean (1994).2.1 IntroductionFor problems with complex geometries, finite-element methods appear to be the naturalchoice due to their intrinsic geometric flexibility. Finite-volume/difference methods, however, are well established in computational fluid dynamics and an alternative approachwould be to use these methods with an appropriate body-fitted coordinate system. Suchmethods can be developed based on well-established solution algorithms and numericalcodes for the regular geometries.18Chapter 2. A Computational Method Using Curvilinear Grids 19Curvilinear coordinates, with orthogonal and non-orthogonal grids, have been usedextensively for fluid flow in complex geometries and a number of research papers exist in the open literature; for examples, see Maliska and Raithby (1984) and Shyy etal (1985). The governing equations are considerably simpler in orthogonal coordinates;however these methods have serious geometric limitations. Orthogonal grids are difficultto generate, especially for three-dimensional domains. For complex three-dimensional geometries, non-orthogonal grids are often necessary because they provide greater flexibilityin the distribution of the grid points.The use of a curvilinear coordinate system can entail a number of difficulties whichhave to be addressed. First, most grid generation techniques provide discrete grid pointsrather than the analytic functions of transformation. In such cases, the coordinate deriva• .. axayaztives of transformation (commonly called covariant base vectors), aj= (—, —, -), areusually computed approximately and other metric quantities, such as contravariant basevectors and the Jacobian determinant, are then evaluated from these coordinate derivatives. The definition and calculation of these derivatives, however, becomes ambiguouswhen a non-smooth grid is used. The procedure used to approximate these quantitiesis critical and may lead to significant numerical errors or even unrealistic solutions, asdemonstrated by Segal et al (1992) and Lee et al (1992) for example.Secondly, significant deterioration of the convergence rate of some of the availableiterative solution procedures can occur with a non-orthogonal grid, as demonstrated byPeric (1990), especially when the grid is highly non-orthogonal. This is partially due tothe explicit treatment of the large non-orthogonal diffusion terms, which can be describeda 128u 0 138uanalytically by cross-derivative terms such as —(Pg-a-—) and --—(Pg --—) in the82momentum equations, where I’ is the diffusion coefficient and g23 is the controvariantmetric tensor.The choice of dependent variables in the momentum equations also requires carefulChapter 2. A Computational Method Using Curvilinear Grids 20consideration. When Cartesian velocity components are used as dependent variablesthe conservation of momentum is considered along fixed directions everywhere in thefield and no curvature terms appear in the momentum equations. The applicability andperformance of the scheme, however, depend on the orientation of the computationalgrid relative to the reference Cartesian coordinate system. In the staggered approach,zero convective fluxes across the control volume faces may occur when the grid turnsby 900 and one Cartesian velocity component is stored on each cell face. Grid-orientedvelocities can be used as dependent variables in order to eliminate this difficulty. The useof grid-oriented velocity components, however, leads to grid-sensitive curvature terms inthe momentum equation. Such terms depend on the second derivative of grid coordinates,and thus are difficult to discretize in a conservative manner, possibly leading to severeinaccuracies as noted by Segal et al (1992).Most of the work reported in the literature uses a variant of the SIMPLE algorithm ofPatankar (1980) to couple the pressure and velocity fields. The extension of this techniqueto non-orthogonal grids results in a complex pressure correction equation, especially forthree-dimensional cases (see Peric (1990)).In the present study a method for the computation of fluid flow in complex geometries is proposed that attempts to solve some of the problems discussed above. Efforts aremade to achieve greater flexibility in grid design using significantly non-orthogonal andnon-smooth grids. To solve the problem related to the non-orthogonal grid, a new secondorder accurate numerical scheme is proposed to describe the cross-derivative terms. Thisscheme allows implicit treatment of part of the non-orthogonal terms, without increasing the size of the computational molecules. To avoid inaccurate discretization for anon-smooth grid, the physical geometric quantities, i.e., volumes, surface areas and surface normal vectors, are calculated directly and used to formulate the numerical schemes.Chapter 2. A Computational Method Using Curvilinear Grids 21These geometric quantities are calculated in such a way as to satisfy the geometric conservation laws as described by Vinokur (1989). When the divergence theorem is used overa control cell, the conservation equation can be expressed accurately for an arbitrarygrid using these geometric quantities. Such an approach yields an overall conservativeapproximation for any grid and provides a better physical understanding of the resultingformulation.The proposed method uses the physical tangential velocity components as dependentvariables in the momentum equations. These variables are the contravariant-type velocityunknowns and are volume flow rates across cell faces with appropriate normalization.Similar to other grid-oriented velocity unknowns, the choice of such dependent variablesgives rise to additional curvature terms. These terms can be expressed by Christoffelsymbols using tensor notation and discretized directly as done by Segal et al (1992), forinstance. In the present study, a different approach is followed which avoids the explicitdiscretization of the second-order coordinate derivatives. A coupled equation solver,combined with a staggered grid approach, is used to solve the momentum and continuityequations directly, eliminating the need for the pressure correction equation.The proposed method is described by first presenting the derivation of the discretizedgeneral governing equations and the numerical scheme for the non-orthogonal terms.This scheme is then proven to be second-order accurate. Next, the treatment of themomentum equations and curvature terms using the tangential velocity components asthe dependent variables is described, followed by a description of the overall coupledsolution procedure in the curvilinear coordinate system. Finally, four computationalexamples are presented to demonstrate the capabilities of the method.Chapter 2. A Computational Method Using Curvilinear Grids 222.2 Numerical FormulationsOnly laminar flows in single domains are considered in this chapter in order to avoidthe complexity of turbulence modeling and multi-domain treatment. The generalizationto those more complex cases will be presented in the following chapters. The flow isassumed steady, incompressible and Newtonian. Under such assumptions, the laminarflow is governed by the following Navier-Stokes equations,eu7u—9(IVu)=—7p, (2.1)(2.2)where g is the flow density, is the dynamic viscosity, p is pressure and u is the velocityvector. Equations (2.1) and (2.2) express the momentum and mass conservations, respectively, for steady incompressible flows. It should be noted that these equations arewritten in a coordinate-free form independent of coordinate systems.The governing equations (2.1-2.2) can be transformed to a general curvilinear coordinate system with the aid of tensor calculus. The transformed equations can then bediscretized in the transformed rectangular domain. In the present study, a direct discretization method is developed. Instead of using the transformed governing equations,the coordinate-free governing equations are used directly to formulate the numericalschemes in a general curvilinear grid. All the numerical schemes are formulated using the physical geometric quantities without using the commonly-used covariant andcontravariant vectors. This method is more desirable for non-smooth grids since thecoordinate derivatives, which are not well defined in irregular grids, are avoided.This section describes the discretization of a general scalar governing equation, whichforms the basis for the discretization of the momentum equations. The geometric quantities used in our discretization are described first. The divergence theorem is then usedChapter 2. A Computational Method Using Curvilinear Grids 23to integrate the governing equation. Finally, the new discretization method for the non-orthogonal terms is described and the second order accuracy is proven.2.2.1 Physical Geometric QuantitiesFor curvilinear coordinates = (r), where r = r(xi, x2, x), is a vector and X: are thecomponents in a Cartesian coordinate system, the commonly-used covariant vector andcontravariant vector in the literature are defined as aj = and a = 7j respectively.These vectors and associated metric tensors are required in the transformed governingequations and subsequent numerical formulation. For computations using curvilineargrids, the covariant vectors are usually calculated approximately from the discrete gridpoints of a given grid. Other quantities required in the numerical formulation suchas the contravariant metric tensor and Jacobian determinant are then evaluated fromthese vectors. In the solution procedure, interpolations are necessary to obtain thesequantities in various locations, especially for a staggered grid system. The procedureof calculation and interpolation of these quantities is critical and may violate certaingeometric conservation laws (see Vinokur (1989)), leading to significant numerical errors.In the present study, instead of using coordinate derivatives aj, the physical geometricquantities, which include cell surface areas, cell volumes and surface normal directions,are used directly to formulate the numerical scheme. A uniform grid with mesh size= 1 is assumed in the transformed computation domain whenever the notations ofcoordinate transformation are used. The geometric quantities used in the present studyare illustrated in Figure 2.1. The unit tangent vectors, denoted by e2 (i= 1, 2, 3), arecalculated at the centers of the control volume surfaces and are locally parallel to thecoordinate lines . These tangent vectors correspond to the normalized covariant basevectors aj in the literature. The surface area vectors, denoted by S (i = 1, 2, 3), aredefined at the same point as e: and are normal to the control volume surfaces, with24Chapter 2. A Computational Method Using Curvilinear Gridsr2Figure 2.1: Illustration of the physical geometric quantities for a control cell.Chapter 2. A Computational Method Using Curvilinear Grids 25magnitude S2 equal to the corresponding surface area. The volume of the control cellis denoted by V.e2, S and V are the basic grid quantities and are calculated directly using discrete gridpoints. For the convenience of formulation, two additional quantities are defined fromthe above geometric quantities. The non-orthogonal angles, denoted by a2 are defined asthe angles between the surface area vector, S2, and the tangential vector e2. These anglesare a measure of the degree of grid non-orthogonality; for an orthogonal grid these anglesare zero. The surface area vectors, S are rescaled and denoted as: e =SS’cosc2.2.2 Calculation of the Geometric QuantitiesIn general, the geometric quantities have to be evaluated from discrete grid points. Theunit tangent vectors e2 are calculated by second order accurate averaging (using thecenter points of two neighboring control cells).The faces of a control cell are generally surfaces rather than planes, as illustratedin Figure 2.2. In fact, it is not always possible to fit a plane through four points. Thesurface area vectors and volume of a control cell can be approximated directly by variousformulas.To approximate the surface area vectors and volume of a cell, the following formulasfor a triangle and a pyramid are used. A properly oriented surface area vector for atriangular face with vertices r1, r2 and r3, is given byS = (r2 — ri) x (r3 — ri)/2, (2.3)and the volume for a pyramid with vertices r1,r2,r3 and r4 is given by(r2 — ri) x (r3 — ri) (r4 — ri) = S . (r4 — ri). (2.4)Using equations (2.3) and (2.4) we can use various averaging processes to calculate thegeometric quantities for an arbitrary cell having straight line edges. Each polygonalChapter 2. A Computational Method Using Curvilinear Grids 26face is divided into plane triangular faces, and the total volume treated as a sum oftetrahedra. In order to obtain an accurate discretization, the geometric quantities arecalculated according to the geometric conservation laws described by Vinokur (1989) andthe following formulas for surface area vectors and volume are used:S, = (r6 — r,) x (r3 —r4)/2 (2.5)V + S + S) (r8 — r,), (2.6)where the subscripts w, s, b indicate locations, west, south, bottom, respectively (seeFigure 2.1), of the vectors. From the two basic vectors, S and e, the quantities coscqand & can be calculated.Noting that r2 —r1, r3 — r1 and r4 — r, are the approximations of covariant vectors a1 =a2= 82and a3=respectively, we can see S1 a2 x a3, and V (ai x a2) afrom the formulas (1) and (2). Since the volume element and contravariant vectors& satisfy formulas = (ai x a2) a3 and a’ = a2xa3 (see Thompson et al (1985)), wehave, 1 and V The volume element is the inverse of the Jacobiandeterminant.The geometric quantities at the positions described in section 2.2.1 are calculatedfrom known discrete points. For other locations, the values for & and V are computedby a second order accurate interpolation and the unit tangent vectors are calculated frome2 by orthogonal relations, & e3 = 0 (i j), instead of calculation from e2 by averaging.2.2.3 DiscretizationThe general transport equation for a dependent variable is written as follows:V (2.7)where u, e and I’, are the velocity, density and diffusion coefficient respectively. The firstterm represents transport of variable by convection and the second term representsChapter 2. A Computational Method Using Curvilinear Grids 27transport by molecular diffusion. The last term, C, is a source term which generallydepends on q. This equation is the natural expression of the transport of variable q andtakes different forms in various coordinate systems. The discretized equations were obtained by first integrating equation (2.7) over a control volume in the physical space, andthen deriving an algebraic approximation to this integral equation. Geometric quantitiesassociated with the control volume are used directly to evaluate the fluxes on the controlcell surfaces.Let J = puqf — I’ denote the total flux. Equation (2.7) can be then written asJ+C—0 (2.8)Equation (2.8) is integrated over a general control volume in the physical space, SV, asshown in Figure 2.2, and the divergence theorem is applied,fJfv.JdvjJ.dS (2.9)J S’e — J S’jw + J S2n — J . S28 + J S — J . S3Next consider transport of variable q along the coordinate direction:JS = qfu.S—I’4.S (210)=whereF=u•S, (i=1,2,3)are the volume flow rates across cell faces. By the definition of gradient,cq5= (jcdS)/VApplying Gauss divergence theorem, we obtain(4’)p = [S’(e—+ S2(q5—q) + S3(q — b)j/V (2.11)Chapter 2. A Computational Method Using Curvilinear Grids 28w.-I II /I I•NFigure 2.2: A general control volume. Nodes to the north, south, east, west, top andbottom are represented by N, S, E, W, T and B respectively. Lower case denotes controlvolume surfaces.Chapter 2. A Computational Method Using Curvilinear Grids 29w ere çSe— , çbn—an—cbb are approximations o ----, ----, and -s-—, respectively,vci uc2 uç3since A2 = 1 in the transformed space. Hence, the gradient of the dependent variable q5can be expressed as:= + S2- + (2.12)Substituting equation (2.12) into (2.10) we obtain the following expression for the transport of variable 4:3. = — pguif- — (2.13)wheresi.si=, (i,j = 1,2,3)is the surface area metric tensor. When the grid is orthogonal, g3 = 0 for i j, therefore, the last term in (2.13) is the result of grid non-orthogonality. The total transportof variable can be decomposed into an orthogonal component and a non-orthogonalcomponent:J•S=J+J, i=1,2,3 (2.14)whereJ = — 1’g-and12—VngSubstituting equation (2.14) into (2.9), the conservation equation over the control cellcan be written as:(J + Jçr)e — (J, +J1)+ (J + Jr)In — (J, + Jr)s (2.15)Chapter 2. A Computational Method Using Curvilinear Grids 30The orthogonal components, J, have the same form as for the Cartesian coordinatesystem. Therefore, schemes such as the hybrid scheme, or the power-law scheme ofPatankar (1980) for regular geometries can be applied to these terms. The non-orthogonalcomponents, however, cause some numerical difficulties. Ordinary discretization of thenon-orthogonal terms, involving only the corner points, results in a 19-diagonal coefficientmatrix which is not unconditionally diagonally dominant. This may result in numericalinstability and unrealistic solutions. Direct solution methods are impractical for this19-diagonal coefficient matrix as demonstrated by Braaten and Shyy (1986). Hence, thecommon practice is to treat the non-orthogonal components of the equation explicitlyby combining them into the source term. This may cause serious deterioration in theconvergence rate if these explicit components become large. It is therefore necessary toinvestigate other treatments. In the present study, a new numerical scheme is proposedto address this difficulty and is presented in the following section.2.2.4 Discretization of the Nonorthogonal FluxesReferring to Figure 2.3, consider the approximation of one of the non-orthogonal terms,g12p, at the mid-point of the east surface. When the central difference scheme is used,there is no contribution to the main diagonal terms of the resulting coefficient matrix.Also, four additional corner points are involved, with at least two negative coefficients inthe resulting discretization for the cross-derivative term, _(g12), making an implicittreatment difficult. The following alternative approximation was introduced to solve thisproblem:If g12e 0, thene(NEE+q’P Suc2Chapter 2. A Computational Method Using Curvilinear Grids 31NW N NE___fl hW w P e ES hSW S hSEFigure 2.3: Illustration of discretization of cross-derivatives.Chapter 2. A Computational Method Using Curvilinear Grids 32and if g12 > 0, then= (N - P + E - SE)/2h (2.17)Different numerical schemes are chosen depending on the sign of the coefficient g’2,based on the preference for a positive contribution to the main diagonal term. The ideais similar to the upwind scheme for the convection terms in which different expressionsare employed for different signs of the velocity. In contrast to the one-side difference, thisscheme is symmetrical around point e and is second-order accurate. This can be shownthrough the following order analysis.With reference to Figure 2.3, the following is obtained using Taylor expansions:NE 75E HEh + jEh2+ 0(h3),andq5p— bs = — ph2 + 0(h3).Combining the above expressions yields,kNE — bE + 4p — q5s = + )h+(E — 4)h2 + 0(h3). (2.18)Since82 +2je + 88h2 + 0(h3), (2.19)4lEP 8621eh+0) (2.20)we can write equation (2.16) as= — E + P — s)/2h + 0(h2). (2.21)A similar estimate exists for approximation (2.17):= — + — sE)/2h + 0(h2). (2.22)Chapter 2. A Computational Method Using Curvilinear Grids 33The above estimates show that the schemes (2.16-2.17) presented here are second-orderaccurate. The other derivatives involved in the expressions for the non-orthogonal components can be treated similarly. This approach provides a second order approximationfor the non-orthogonal term Also, using this method, the number of cornerpoints used is reduced by half and the main diagonal term of the resulting coefficientmatrix, ap, is augmented by the following amount:2P(g’ + jg’3 + g23).2.2.5 Formulations of the Discretization CoefficientsAfter treating the non-orthogonal diffusion as described above, the general governingequation can be discretized following the methods developed for standard Cartesian coordinates. In the present study, the Power Law profile of Patankar (1980) is used for theorthogonal fluxes which ensures the central difference operator at low Peclet number andgradually changes to upwind differencing at high Peclet numbers.A general algebraic equation is obtained by substituting the discretized forms of theorthogonal and non-orthogonal fluxes into equation (2.15). The resulting equation canbe written asapqp = aEbE + awqw + aNqN + ascbs + aTqT + aBçbB + ab + b (2.23)where the aE, aw, etc. denote the combined convection-diffusion coefficients, includingthe non-orthogonal terms. The summation index “nc” represents the corner points fromthe discretization of non-orthogonal fluxes (such as c’wN, and q5ws) and b includes onlythe source term C.Chapter 2. A Computational Method Using Curvilinear Grids 342.3 Calculation of The Flow Field2.3.1 Variable Arrangement and The Choice of Velocity UnknownsA staggered grid arrangement is adopted in which the scalar quantities, such as pressureare located in the geometric center of each control cell and the tangential velocitycomponents are located in the middle of the cell surfaces. There are several optionsfor choosing velocity unknowns for the method using non-orthogonal grids. The use ofcurvilinear velocity unknowns is preferred here since the applicability of such a methodis not limited by grid-orientation, which is a requirement for later calculations.Different sets of velocity unknowns may be chosen. A discussion and review of thevarious possibilities is given by Rodi et al (1989). In the present study, the physicaltangential velocity components are used as the dependent variables for the momentumequation. These variables, denoted as U, as defined by Borisenko and Tarpov (1968),are the resulting coefficients from the velocity expansion in the unit tangent basis vectorsu = U’e1 + Ue2e2 + Ue3. (2.24)From the above equation and the orthogonal relation of {ej and {et}, the followingexpression for U can be obtained:UE = u e’ (2.25)eFrom the above formula, it can be seen that the physical tangential velocity componentsare similar to the physical contravariant velocity components u & which were usedby Demirdzic et al (1987) as the primary velocity unknowns and differ from them onlyby a factor of e: e. Unlike the non-physical contravariant velocity components, thesephysical velocity unknowns have the same order of magnitude as the Cartesian velocitycomponents which are independent of the mesh size. As pointed out by Demirdzic etChapter 2. A Computational Method Using Curvilinear Grids 35Figure 2.4: Illustration of the three types of physical velocity components in curvilinear grids, OA1 and 0A2: physical contravariant velocity components, OB1 and 0B2:physical tangential velocity components, and OC1 and 0C2: physical covariant velocitycomponents.al (1987), the use of physical velocity unknowns is preferred since there are drawbacksinevitably associated with the non-physical velocities. There are three different typesof physical curvilinear velocities which are illustrated in Figure 2.4. In this figure, theline segments OA1 and 0A2 represent the physical contravariant velocity componentsand are the velocity projections on the unit surface-normal vectors. The line segmentsOB1 and 0B2 represent the physical tangential velocity components and are the resultsof velocity expansion in the tangential vector basis. Segments OC1 and 0C2 representthe physical covariant velocity components and are the velocity projections on the unittangential vectors e1 and e2. The advantage of using the physical tangential velocitycomponents is that if the velocity is tangent to one of the grid lines, the correspondingphysical tangential velocity component will become identical with the velocity whileother tangential velocity components vanish. From Figure 2.4, it can be seen that onlye2A2 UAlChapter 2. A Computational Method Using Curvilinear Grids 36the physical tangential velocity components have such a property, which reduces theeffects of false diffusion when one of the grid lines is aligned with the flow streamline.From equation (2.25) it can be seen that the physical tangential velocity componentsare the volume fluxes, u S, normalized by appropriate geometric quantities. The volumefluxes m S were used as velocity unknowns by Rosenfeld et at (1988), and are equivalentto the quantities As reported by Segal et at (1992), the use of quantities asvelocity unknowns provides a more accurate solution than the use of the contravariantvelocity projections U.The tangential velocity components are contravariant unknowns which are equivalentto the quantities . Unlike the contravariant velocity projections U or volumeScosafluxes, U’ are the physical components and their physical lengths are of the same order ofmagnitude as the velocity. This facilitates the implementation of the boundary conditionin programming.Similar to the use of volume fluxes, the use of the tangential velocity componentssatisfies the velocity-recovery requirement as described by Segal et at (1992). In fact,from the calculations of geometric quantities described in section 2.2.2 and the formula(2.24) and (2.25), it is easy to see that the transformation u —* U —* v gives exactlyv = u if u is a constant vector. Computational tests also show that the use of thetangential velocity components as velocity unknowns provides satisfactory performance.2.3.2 Discretization of Momentum EquationsAs with other types of grid-oriented velocity components, the use of tangential velocitycomponents as dependent variables gives rise to additional curvature terms which canbe expressed as Christoffel symbols using tensor notation. Since these curvature termsinvolve second order derivatives of the grid coordinates, they are difficult to discretizeaccurately when grids become non-smooth. There are two basic approaches for theChapter 2. A Computational Method Using Curvilinear Grids 37discretization of momentum equations for the use of curvilinear velocity unknowns. Onemethod is based on the transformed momentum equations for the curvilinear velocityunknowns. Discretization is carried out in a rectangular domain in the transformedspace. The other is based on the discrete momentum equations written in a locallyfixed coordinate system. The present study follows the second approach and derives thediscretized equations based on algebraic equations for Cartesian velocity components.This approach was first used by Karki and Patankar (1988) for calculations of two-dimensional flows using the physical covariant velocity unknowns as the primary variablesin the momentum equations. It seems more desirable for non-smooth grids since theexplicit discretization of second order derivatives of the grid coordinate is avoided.In order to discretize the momentum equations, auxiliary discretizations for the Cartesian velocity components are considered. Suppose that all three Cartesian velocity components are located at the U’ position. The governing equations for Cartesian velocitycomponents can be viewed as special cases of the general governing equation (2.7). Thediscretization can be obtained according to the method outlined in the last section byintegrating the governing equations over the control cell for Ut’. Since the Cartesian yelocity components are assumed to share the same control cell as for Us’, their discretizedequations will have identical coefficients. Therefore, the discretization can be written,using vector notation, as:a’u = a’ue + au + a$L’ufl + a’u3 + + a’ub + + be’, (2.26)where the lower case is used for the subscripts instead of the upper case to indicate thatthe positions of UC1 are at cell faces instead of at cell centers. The source term be’,contains only the pressure gradient term, which is discretized by applying the divergencetheorem, as discussed later. The above discretization for the momentum equation iscoordinate-invariant and independent of the choice of dependent variables.Chapter 2. A Computational Method Using Curvilinear Grids 38To find the discretization using the tangential velocity components as dependentvariables, the velocity expansion in the local tangential vector basis at point “p”,is considered. Taking the inner product of vector e and vector equation (2.26) andapplying formula (2.25), the following equation is obtained:= aj,(UE) + + (be’)’, (2.27)where(U,)’ = e. u, (U$j’ = e• u, (be’)’ = e be’, (2.28)the index “nb” represents the six nearest neighbors of the node p, namely, e (east), w(west), n (north), s (south), t (top) and b (bottom), and the index “nc” represents thecorner points, as previously defined. The primed velocities, (Us)’ and (Us)’, are velocityprojections of neighboring velocities over the vector, e’, at point p. The vector e’ changesfrom point to point in the flow field, and therefore, the velocities (U,)’ generally differfrom the actual neighboring variables, U?. Since (U,)’ are not dependent variables, theymust be replaced by U,. This can be done by rearranging equation (2.27) as follows:a’U’ = aU + a$U,! + b1,c + bE1 (2.29)where= > a((U,!)’ — U?) + a((U,)’ — U,). (2.30)Equation (2.30) represents the curvature terms when the tangential velocity componentsare used as velocity unknowns. Using equations (2.25) and (2.28), equation (2.30) canbe rewritten as= a(eb — e) . u,1, + — e) u. (2.31)The above equation shows that the curvature terms for U’ are produced by the change ofvector e1 from point “p” to the neighboring points. Figure 2.5 illustrates the formationof curvature source terms.Chapter 2. A Computational Method Using Curvilinear Grids 39Figure 2.5: Velocity decomposition in a locally fixed coordinate system.For an actual computation, the velocity u at neighboring points is calculated fromequation (2.24), and the curvature term b1,c is then calculated from equation (2.31). Fora staggered grid the velocity unknowns U2 and are not defined at the UE1 positionTherefore, some interpolation is needed for the calculation of velocities Unb and u,. Inthe present study, full weighting is used for the interpolation.The discretized equations for the other two velocity unknowns can be derived in asimilar manner. The resulting discretization for the momentum equations, with part ofthe pressure difference term written explicitly, can be expressed as follows:= + + b1,c + A1(pp—pw) + b’ (2.32)a2U$2 = + a$U, + b2,c + A2(pp—PS) + b2 (2.33)= aU + a$L’U7 + b3,c + A3(pp—PB) + bE3 (2.34)a 2,e 7 Ui,eUe(a)Ui,e(b) alewhere the superscripts indicate the different coefficients associated with different velocityChapter 2. A Computational Method Using Curvilinear Grids 40components and the second summation on the right hand side includes partial non-orthogonal flux terms. The third term represents the curvature body force and the lastterm contains only pressure terms.2.3.3 Treatment of Pressure Gradient and Associated BoundaryIn equations (2.32-2.34), the pressure term is concealed in the source term. However,the pressure is also an unknown in a flow field and must be found together with thevelocity variables. The pressure term, b”, in these equations can be obtained from theabove derivation of the discrete momentum equations. The p can be expressed usingthe physical geometric quantities as follows,v(S1+S2+S3f). (2.35)The pressure source term, bi1’, in equations (2.32-2.34) can be written explicitly as=+ (2.36)whereg”=S•S2, (i,j=l,2,3)is the surface area metric tensor. The first term in equation (2.36) is the orthogonalcomponent, and the last term in equation (2.36) is the non-orthogonal component whichvanishes for an orthogonal grid. The orthogonal component is discretized naturally usingthe central difference scheme. For the non-orthogonal component, the scheme proposedfor the non-orthogonal diffusion terms is applied. Referring to Figure 2.3, the discretization of g12f at point ‘e’ takes the following form.12 — f g’2(pNE —PE +pp —ps)!2 ifg12 <0guc2.g’2(pN — pp + PE — psE)/2 otherwiseChapter 2. A Computational Method Using Curvilinear Grids 41With such a discretization, a quantity g’2e(pE — pp)/2 is generated which appearsas a direct pressure force to drive the velocity component Ut’. This quantity, which isequivalent to an orthogonal component g12/2, can be treated in a similar manner asthe orthogonal term. Such a treatment of non-orthogonal pressure terms will increase theimplicit portion in an iteration procedure and contribute to the main diagonal terms in theresulting pressure equation, thereby enhancing the numerical stability and convergencespeed.The discretization of the pressure gradient requires pressure values outside the flowdomain. This can be seen from the discretization of non-orthogonal pressure termspresented in the above formulation. For example, the discretization of cross-derivativeg120 in the Ue1 — momentum equation next to the north boundary requires the pressurevalues beyond this boundary. In the present study, a linear extrapolation is used tocalculate the pressure values outside the boundary using the pressures inside the domain.2.3.4 Discretization of Continuity EquationsIntegrating the continuity equation7 •gu = 0 (2.37)over a scalar control volume, as shown in Figure 2.2 and applying the divergence theorem,we obtain(u — (gu S’) + (eu S2)7, — (eu• S2)8 + (eu. S3) — (eu. = 0. (2.38)Substituting equations (2.24) into the above equation yields the following discretizedequation:(iU’)e — (ciU’) + (c2U2) — (C2Ue2)5+ (c3UE3)t — (c3U)b = 0, (2.39)Chapter 2. A Computational Method Using Curvilinear Grids 42where,Figure 2.6: A scalar control volume with 7 unknowns.Ct = eSlcoscxWith this formulation, the mass conservation equation is expressed exactly using thetangential velocity components, without any extra source terms.2.4 Solution Procedure2.4.1 Coupling of the Velocity and PressureThe evaluation of the pressure field has always been a difficult issue in the primitive-variable approach for incompressible flow since the pressure is indirectly involved in thecontinuity equation. This difficulty becomes more acute when a curvilinear coordinate,/IChapter 2. A Computational Method Using Curvilinear Grids 43system is used since expressions for the pressure gradient and the continuity equationbecome much more complicated. In this section, a block-implicit pressure-velocity coupled solution procedure is presented which essentially updates the velocities and pressuresimultaneously. This method was first used by Vanka (1986a) to solve two-dimensional,steady, incompressible flow using a Cartesian coordinate system and was observed toprovide good convergence rates.2.4.2 Symmetric Coupled Gauss-Seidel Iteration MethodThe pressure and velocities for a typical control volume are shown in Figure 2.6. Thepressure, located at the center, and the velocities, located at the control volume surfaces,are treated as the unknowns, and the pressure and velocities at all the other pointsare treated explicitly. The momentum equations (2.32-2.34) provide the six algebraicequations for the six unknown velocities. These equations can be simplified by treatingonly two variables implicitly for each equation and are written as follows:a’U’ — App = R1, (2.40)a’eU’ + A’pp = 1?2, (2.41)a28U — A2pp= R3, (2.42)a2U12 + A2pp = R4, (2.43)3 e3 ,iE3 —b — ‘b PP — .LL5,+ App = R6, (2.45)where the R (i = 1, 2, ...6) include all the other terms left in the equations, calculatedusing the currently-available values of the velocity and pressure at the neighboring points.Chapter 2. A Computational Method Using Curvilinear Grids 44This approach is intended to provide a simple algebraic system which can be solvedefficiently. The continuity equation (2.39) provides an additional algebraic equation, i.e.,Ci,eUj’ — + C2,U,2 —c2,3U2 +c3,tU! —c3,bU = 0. (2.46)These equations are arranged in a block structure as follows:4’ 0 0 0 0 0 —A U5j R1o a’6 0 0 0 0 A!’ U’o 0 0 0 0 — A!2 U2 R3o 0 0 a2 0 0 A2 X U2 = R4o 0 0 0 0 —A Ubo 0 0 0 0 R6Cl, Ci,w C2,n —C2,s C,t C3,b 0 PP 0The above block of equations is solved analytically using Gaussian elimination. Theprocedure is repeatedly applied in the checkerboard order to all cells in the flow field.After sweeping the whole field, the coefficients a, etc., are recalculated and the entireprocedure is repeated until the residues become sufficiently small.2.5 Validation of the MethodThe proposed method has been implemented in the CMGFD code. This code was writtenfor general three-dimensional geometries, and can treat both two- and three-dimensionalproblems. A number of examples are reported in this chapter. They are computed usingthe CMGFD code to validate the code and the proposed method using non-orthogonalgrids.2.5.1 Cavity flow with inclined side wallsC.nCD IszI-i CD0CD 0 C, CD C, CD CD0C)0(7’Chapter 2. A Computational Method Using Curvilinear Grids 46Grid 20 x 20 40 x 40 80 x 80Inclination /3 600 45° 30° 60° 45° 30° 60° 45° 30°# of Iterat. 75 85 110 141 160 190 310 396 450Table 2.1: Number of iterations necessary for convergence.The flow in a cavity with the top wall moving at a constant velocity, U.n,, is often used asa case study for testing the efficiency of a solution algorithm. Here, the side walls of thecavity are taken to be inclined, forming an angle, 3, with the horizontal plane. Threecases are computed with inclinations of /3 = 60°, 45° and 30°, at a Reynolds numberRe = 100, where the Reynolds number is defined as Re =.Calculations wereperformed for various grid densities and Table 2.1 lists the convergence rates for eachtest calculation. The convergence criterion was a maximum residual of less than i0 foreach equation. The momentum residuals were normalized by and the mass residualby oU. Unlike the results reported by Peri (1990), the present method allows a widerange of under-relaxation factors with fast convergence for significantly non-orthogonalgrids. Computational results using the finest grid, 80 x 80, are plotted in Figure 2.7.They show a strong main vortex, driven by the lid movement, and a sequence of weakervortices in the sharp corner between the bottom plate and the upstream side wall. Thissecond vortex system is the main difference between the square cavity and cavities withinclined walls. These results are very similar to the previous numerical results of Perié(1990), which used a central difference scheme for non-orthogonal fluxes.2.5.2 Laminar flow through a tube with a constrictionAs the second example, laminar flow through a tube with an axisymmetric constrictionwas considered. Such flows were studied experimentally by Young et al (1973). TheCDCDoi-. CDo t-.÷ocCDC4Zcox-‘ 0C0CD CDjCD0 0CD0p,CD CD CDCD CD—0CDI—coCD—.—.nOq0CDCDCD CDCDCD0 I-iD’ICDCDcnI-—CD-00ACDCDA0 —0CDCD0‘-1 CDO 0 c3 )-1 CDCD—.I-i—-CDC, I-1 CD I000bb,o0-oco0C)c3-0 0 >< C’) Cl) CD 0 a CDHWU jJ,llllll11r)C).O0).4(Vo---r’)t,rn11q1mpChapter 2. A Computational Method Using Curvilinear Grids 4880 x 20 Separation point Reattachment pointRe # of Iterations Experiment Prediction Experiment Prediction50 210 0.33 0.32 2.28 2.27100 241 0.34 0.33 4.19 4.10Table 2.2: Comparison with experimental results and numbers of iterations for calculation.The boundary conditions imposed were no-slip on the wall, zero stream-wise gradientat the outlet, fully developed parabolic flow at the inlet, and symmetry conditions at theaxis. An 80 x 20 grid in x-r coordinates was used for the calculations (see Figure 2.8a).The density of the grid points was higher near the wall than at the centre of the tube,and the grid was stretched in the axial direction, with more grid points in the constrictedregion. Computations were carried out at Re=50 and 100 where the Reynolds numberwas defined as2UR0Re=iiand U is the mean velocity of the inflow. The predicted flow field and streamlines forRe = 50 are plotted in Figures 2.8b & 2.8c respectively. In these figures, X8 and X,. areseparation and reattachment points respectively. The flow pattern for Re 100 is verysimilar, except for the size of the recirculation zone.The predicted separation and reattachment lengths, together with experimental results from Young et al (1973), are presented in Table 2.2. Excellent agreement is shownbetween the present computations and experimental results.2.5.3 Developing laminar flow in a pipeFor the third computational example, laminar flow through a pipe was considered. Thisflow can be considered to be two-dimensional because of axial symmetry. However, theChapter 2. A Computational Method Using Curvilinear Grids 49use of a non-orthogonal grid in the present calculation makes the problem appear fullythree-dimensional (see Figure 2.9a). This problem was chosen because the availableanalytical solution makes it possible to test the accuracy of the proposed method. Furthermore, such a grid will be used as part of the mesh for subsequent computations ofthe film.cooling process for turbine blades.The 16 x 16 x 50 mesh illustrated in Figure 2.9a was generated by solving an ellipticgrid generation system of equations.The three-dimensional, incompressible Navier-Stokes equations were solved, with theno-slip condition at the wall, uniform velocity profile (w = 1) at the inlet, and zerostream-wise gradient condition at the outlet. The pipe length was chosen as twelve timesthe diameter (D) for the computations and the Reynolds number was ReD = 100. Usingthe same convergence criterion as in the previous problem, a solution was obtained after480 iterations. Figure 2.9b shows the velocity vectors in the symmetry plane. Figure 2.9cshows the calculated velocity contours in a cross-plane located ten diameters downstream(x/D = 10). Figure 2.9d compares the calculated solution in the well-developed region(x/D = 10) with the analytical solution u(r) = 2.0(1 — i’)2. The result shows goodagreement between the analytic and the numerical solutions. Figure 2.9e compares themaximum velocity in the developing region with the results reported by Shah and London(1978) which also shows a good agreement.2.5.4 Flow in a pipe with a smooth 90 degree bendFor this example, laminar flow in a pipe with a smooth 90° bend was considered. This isa strongly three-dimensional flow, where the main flow in the streamwise direction alongthe pipe is influenced by strong secondary flows in the pipe cross section, arising fromthe centrifugal forces due to the bend curvature. This flow was studied experimentallyby Enayet et at (1982). The same geometry used in the experimental study was adoptedChapter 2. A Computational Method Using Curvilinear Grids 50Level w7 1.956 1.835 1.504 1.173 0.842 0.511 0.18Figure 2.9: (a) Illustration of grid and velocity, (b) development of velocity profile, (c)velocity contour in the cross-pipe plane at x = 1OD, (d) comparison with analytic solutionat x = 1OD, (e) Comparison for flow in the developing region, where - =(c)1,A. *i1ii 1iittIl•IiiitflflRtflh1P1fl1ti1i1111m11iI 11111lliiiiII:iiiu1i11II(a) (b)2.01.51.0:: . .-0.50 -0.25 0.00 0.25 0.50(d)2.00U/Ubcompanal0.025 0.050c)0.075 xfDRe-4.CDI-iCDUllll1UllllBIIWUllS1-4.‘I.qIllhluhlIuI(I)_hu1lluhbullllll1Mh1Ui011111111111111gUllllffluhb’‘UllilMilli’C, -4.0CD‘J1llMlli’CD1blllli’_0‘iI•lli’00CD-4.CD330 oCD(0CD0 CDH-HI.Chapter 2. A Computational Method Using Curvilinear Grids 52Figure 2.12: Contours of mean velocity U/UB in the bend.for the present computation. The geometry and grid are illustrated in Figure 2.10. Dueto the symmetry about the x-z plane, the computational domain has been limited to thetop half of the pipe only. A 32 x 16 x 64 mesh was used in the computations, which werecarried out at Re = pUbD/u = 500, corresponding to the experiment, where Ub is thebulk velocity at the entrance and D is the pipe diameter. The convergence criterion wasthe same as that used in the first example, with convergence satisfied after 620 iterations.Figure 2.11 shows the calculated velocity vectors in the plane of symmetry. Contoursof streamwise velocity in four cross-stream locations are plotted in Figure 2.12. Thepredicted velocity at these same four locations in the plane of symmetry, together withthe results measured by Eanyet et at (1982), are plotted in Figure 2.13, where reasonableagreement can be seen.Inside a) 3=3O Inside b) 13=60Inside c) 13=75 Inside d) h=DChapter 2. A Computational Method Using Curvilinear Grids 53u.2OO71a) f3=300.5 xJD 0.7—a------ comp-----s---- exp1.5w1.00.51.5w1.00.50.2 1:0b) f=6O1.5w1.00.5c) J3=750.2 0.5 0.7d) h=D10Figure 2.13: Comparison with the experimental results of Eanyet et al (1982).Chapter 2. A Computational Method Using Curvilinear Grids 542.6 ClosureA numerical method for computing laminar flows in complex three dimensional domains ispresented. Computational examples show that the proposed method can be used to solveflow problems in complex three-dimensional geometries using significantly non-orthogonaland moderately non-smooth grids. A new scheme for handling the non-orthogonal fluxesis proposed which is useful when the grid is significantly non-orthogonal. Efforts havebeen made to allow for the use of moderately non-smooth grids by directly employing thegeometric quantities of control cells and avoiding explicit discretization of the derivativesof the grid coordinates.Chapter 3Computation of Turbulent Flows Using Curvilinear GridsIn this chapter, the method developed for laminar flows in the last chapter is generalized to solve turbulent flows in complex three-dimensional geometries. The Reynolds-averaged Navier-Stokes equations together with the k — e turbulence model are solved tosimulate turbulent flows. The cwall function’ is used to take into account the near-wallregion where the k — 6 model is not valid. The discretization of the Ic — e equationsand the treatment of the source terms are discussed. Special attention is given to thetreatment of wall boundaries. The performance of the numerical method is investigatedthrough several three-dimensional turbulent flows. Comparisons with available experimental results are made which show that the present solution method for turbulent flowsis reliable.3.1 IntroductionIn this chapter, the ideas used in chapter 2 are generalized to solve turbulent flows incomplex three-dimensional geometries. It is well known that turbulent flow is moredifficult to handle than laminar flow. Turbulence modeling is one of the difficult areasin CFD. However, most engineering flows, such as film cooling of a turbine blade, areturbulent. This chapter presents work intended to increase our knowledge of computingturbulent flows in the environment of general curvilinear grids and to develop an efficientsolver for turbulent flows.55Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 56The 3D, incompressible, time-averaged Navier-Stokes equations are solved. Turbulence closure is attained by the use of the standard k — e model with the ‘wall function’treatment described by Launder and Spalding (1974). The k — e model is the most widelyused turbulence model in the CFD community. Its capability has been demonstrated bymany authors. However, the k — c model also has some shortcomings for turbulent flowswith curvature, recirculation, etc. due to anisotropic and non-equilibrium turbulence. Itis our aim to show how the k — e model performs with our solution method for curvilineargrids.The k — e turbulence model was derived for fully turbulent flows and is not validin the near-wall region in which the influence of laminar viscosity is important. In thepresent study, the ‘wall function’ described by Launder and Spalding (1974) and othersis used to link the equations in the fully turbulent flow and in the near-wall region. Themethod evaluates the turbulence quantities near the wall by assuming different velocityprofiles (log-law or linear) in the near-wall region. The advantage of such a method isthat there is no need for excessively fine grids to take into account the steep gradientsin the near-wall region. In this chapter, the detailed formulations of the ‘wall function’using the grid geometric quantities and the physical tangential velocity components arepresented. The velocity boundaries on the wall are implemented by treating the wallshear stress as an auxiliary force in the momentum equations.For turbulence calculations, the k — e equations do not generate additional complexityfor discretization. They can be viewed as special cases of a general scalar equation anddiscretized according to the method described in previous chapter. However, the numerical treatment of the source terms in these equations is important and deserves specialattention. In this chapter, the treatment of these source terms is reported in somewhatmore detail. In non-orthogonal curvilinear grids, the turbulence energy generation ratebecomes difficult to calculate from the tangential velocity unknowns. To overcome thisChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 57difficulty and facilitate the calculation of the curvature source terms in the momentumequations, the Cartesian velocity components are used as auxiliary dependent variables.The turbulence energy generation rate can be easily calculated using the Cartesian velocity components.The performance of the developed method is investigated through several three-dimensional turbulent flows, including benchmark computations for developing turbulent pipe flow, and developing turbulent flows in a square duct and a tube, each witha 900 bend. Computational results are compared with available experimental results tovalidate the method.3.2 Governing EquationsIn this section, the governing equations and turbulence model used in the present studyare described. The full, Reynolds-averaged Navier-Stokes equations together with thestandard k — e equations are solved for turbulent flows. Using the eddy viscosity concept and the closure of the k — c two-equation model of Launder and Spalding (1974),the averaged governing equations for steady, incompressible flows can be written in thefollowing coordinate-free form:u.vu—v.(jieiivu)—vp, (3.1)vu = 0, (3.2)v(euk—vk)=G—ee, (3.3)2[It 6 6v .(eue — — v c) = C1G — 0C--, (3.4)fl, lbwhere k is the turbulence kinetic energy, c is the turbulence energy disspation rate, eis the flow density, G is the turbulence energy generation rate, C1, C2, k and a6 areChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 58empirical constants which, following Launder and Spalding (1974), are taken asC1 = 1.44, C2 = 1.92, k = 1.0, o = 1.3, C,1 = 0.09.Equations (3.1) and (3.2) represent the momentum and mass conservations, respectively,for steady, incompressible flows, and are the well-known Reynolds-averaged Navier-Stokesequations. Equations (3.3) and (3.4) are the transport equations for turbulence kineticenergy and its dissipation rate respectively. In these equations, u is the mean velocityvector which is time-averaged from the instantaneous velocity, u is a tensor, and neffis the effective viscosity which is given by:Peff = ILt + i’i,where tj and i-it are the laminar and the turbulent viscosities respectively. The turbulentviscosity is evaluated from the relationpt = pC,k2/c, (3.5)where C,4 was found empirically to be approximately constant at high Reynolds numbers.3.3 Discretization of the k — e EquationsThe momentum equations and the continuity equation are same as that for laminarflows presented in Chapter 2 except that the viscosity ueff in equation (3.1) includes theadditional turbulent viscosity. Therefore, the velocity components and pressure can betreated in the same way as in the last chapter. The discretization of equations (3.1-3.2)can be obtained from the equations derived in the previous chapter. Therefore, only thediscretization of the k — € equations requires further consideration.The k — 6 equations (3.3-3.4) can be treated together as a general scalar equation asfollows:(3.6)Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 59where the dependent variable, q, is a scalar which designates the turbulence energy k ordissipation rate c. The corresponding diffusion coefficient and source terms in the aboveequation areF=, S=G—e, (3.7)and1’ = , S = CG — C2T. (3.8)06for the k and 6 equations respectively.The discretization of the scalar equation (3.6) has been presented in the last chapter.Therefore, the discretization of the k — e equations can be obtained without additionaldifficulties if the source terms can be solved appropriately. There is not much informationin the literature regarding the efficient treatment of these source terms. It was found,however, that the treatment of these source terms is critical in terms of computationalstability and efficiency. In the present study, the numerical treatment of these terms isreported.The source terms in the k — 6 equations include the turbulence energy generationrate which has to be evaluated from the velocity. The evaluation of the turbulenceenergy generation rate G becomes complex for computations using curvilinear velocityunknowns. In a Cartesian coordinate system, the generation rate G can be expressedusing the derivatives of mean velocity components:aUi 0U2 c9UC= + —)----, (3.9)where (U1,U2,U3) are the Cartesian mean velocity components. In a curvilinear coordinate system, the generation rate is usually calculated using the transformational formula.The expression of the turbulence generation rate based on the physical tangential velocitycomponents can be obtained with the aid of tensor calculus. However, the formulationChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 60includes the second-order derivatives of grid coordinates and becomes complex and difficult to calculate. To save computational efforts, Cartesian velocity components areused as auxiliary variables in the solution procedure in addition to the primary velocityvariables. Such an approach allows efficient calculation of the turbulence generation ratewithout using the transformed formulas for the physical tangential velocity components.In addition, the curvature source terms in the momentum equations can be efficiently andaccurately evaluated. The disadvantage is the requirement of additional computer storage for these auxiliary variables. This disadvantage can be remedied using a subdomainsolution procedure as introduced in a later chapter.The generation rate can be evaluated directly from the auxiliary Cartesian velocitycomponents as follows. First, the gradients of Cartesian velocity components VU arecalculated from the formula for presented in equation (2.12). From the gradientVU, the derivatives 0U2/6x1 can be obtained easily and the energy generation rate Gcan then be calculated from equation (3.9). Such an approach allows efficient evaluationof the energy generation rate.After calculating the turbulence energy generation rate, the source terms can betreated in the following manner. For the turbulence energy equation, integrating thesource term Sk over a control volume SV givesS =—oe) r (G — 0e)V (3.10)where the subscript p indicates that the value is associated with the control-volume centerand V, is the volume of the control cell with center ‘p’. There are two approaches whichcan be used to deal with the source term. One is to move all the quantities in equation(3.10) to the right-hand side in the discrete equation and treat them as source terms inthe solution procedure. The other is to split the source term into two different portionsby appropriate linearization and to treat one of them implicitly. For the first approach,Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 61however, the source term may become negative. As discussed by Patankar (1980) for analways-positive variable, the source term in the discretization equation should be positivein order to avoid unrealistic solutions. To achieve a positive source term, the dissipationrate c is replaced by C,Lek2/uL using equation (3.5) and the negative term —Ge can thenbe linearized with respect to the dependent variable k. Therefore, the total source termcan be split as follows:(3.11)where(Y 2c’1Li,& C” IrYTI\)p, cV)pThe term Sr., positively contributes to the main diagonal term ap in the discretization ofthe k-equation, which enhances the numerical stability. The source term of the c-equationcan also be linearized in a similar manner according to the always-positive rule describedby Patankar (1980).After calculating the source terms, the final discrete equations for k and c can bewritten as the following seven-point algebraic equations:4kp = 4kE + 4kw + 4kN + agks + 4kT + aB + bc, (3.13)aep = aECE + aew + aNEN + acs + aTET + a + b, (3.14)where the subscripts indicate the seven scalar positions, namely the center point ‘P’ andthe six neighbors; the superscripts indicate coefficients associated with k and e.For a non-orthogonal grid, the source terms in the above two equations contain thenon-orthogonal diffusion quantities which may become negative which allows the possibility that the physical positive variables could acquire an erroneous negative value. Tosolve this problem, an artificial damping quantity is added to both sides of the discreteequations. Taking the k—equation as an example, a positive quantity 3k is added toChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 62both the source term and the main diagonal term at the same time, where /3 is a positiveconstant and is chosen depending on the problem to be solved. Therefore, the sourceterm is augmented by an amount I3kp while the coefficient of the main diagonal term,ap, is increased to ap + /3. This artificial damping technique is different from the ordinary under-relaxation in the literature. The former increases both the main diagonaland source terms and the latter reduces the change of the dependent variables. Such apractice is found to be very useful for turbulent calculations using non-orthogonal grids.3.4 Wall BoundariesThe k — 6 model is valid only for fully turbulent regions and does not take into accountthe laminar viscosity effect which is important in the near-wall region. To take this effectinto account without using an excessive number of grid points near the wall, the ‘wallfunction’ is adopted in the present study.The ‘wall function’ has been used in the literature for turbulent flow calculations incurvilinear coordinates. Demirdzic et al (1987) presented a calculation procedure for two-dimensional turbulent flows based on the physical contravariant velocity unknowns. Aformulation of the ‘wall function’ for two-dimensional curvilinear grids was also presentedin this paper. Agouzoul et al(1990) presented a calculation for a three-dimensional turbulent flow using a non-staggered and non-orthogonal grid based on the Cartesian velocityunknowns. To avoid the complication associated with the formulation of momentumequations near the wall for three-dimensional flows, an ‘equivalent viscosity’ approachwas proposed in the paper to take the variations of the velocity in the near-wall regioninto account. In the present study, the velocity boundary at the wall is implementedby taking the wall shear stress as an auxiliary force in the momentum equations andappropriately altering the flow fluxes on the control-volume surfaces near the wall. TheChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 63formulation in a general curvilinear coordinate system based on the physical tangentialvelocity unknowns and the grid physical geometric quantities is described below.We consider the formulation at a scalar node Cp next to the wall. Following Launderand Spalding (1974), two different velocity profiles near the wall are assumed dependingon the values of the dimensionless distance y+,= lm(Ey), if y > 11.36, (3.15)= y, if y < 11.36 (3.16)where the dimensionless distance y+ is defined as:=ILand the friction velocity is defined as:UT=(3.17)Here, y, is the normal distance of the node ‘p’ from the wall, U, is the velocity componentparallel to the wall at the node ‘p’, and i and E are empirical constants which are takenas i = 0.41, E = 9.97 following Launder and Spalding (1974).The normal distance, yr,, is evaluated directly from a given grid. To calculate U,, thevelocity vector u which is calculated from formula (2.24) is decomposed into the normalvector, u, and the tangential vector, ui,, at the wall,u = (ew . up)e’, u = u, — u (3.18)where et° is the unit surface normal vector at the wall. Here the magnitude of thetangential velocity vector u is U.After U, and y, are obtained, the only unknown in equation (3.15) or (3.16) is thefriction velocity UT. Therefore, UT can be solved from these equations iteratively usingChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 64Newton’s method. The values of k and c at the first node from the wall are then specifiedby the following algebraic expressions:= (3.19)= T (3.20)pvInstead of using the ‘no-slip’ wall condition, the wall boundary for velocity is implemented by appropriately modifying the flux transport terms at the cell surfaces adjointto the boundary and then treating the wall shear stress as an auxiliary force in the momentum equations. The force due to the wall shear stress is in the opposite directionto the tangential velocity u and can be expressed as: —rAu/U, where A is the wallsurface area of the control cell which can be obtained from the surface area vector, andthe wall shear stress can be obtained from the friction velocity using equation (3.17).This term is treated as a force in the momentum equations and is used to modify thesource terms in the discretized equations. Taking the bottom wall as an example, thesource terms due to the wall shear stress for the Ue1 and UC2 momentum equations areejS3/U and eSj/U, respectively.3.5 Application to Turbulent FlowsThe method described in this chapter has been implemented in the CMGFD. In thissection, the code is applied to several three-dimensional problems to investigate theperformance of the developed method and the k — c model.3.5.1 Developing Turbulent Pipe FlowIn the first example, the turbulent developing flow in a smooth pipe was computed. Sucha flow has been studied extensively by many authors. A review paper by Klein (1981)Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 65summarized most of the experimental results for the turbulent developing pipe flow. Aspointed out by Klein, developing pipe flow is very complex and, unlike earlier studieswhich believed that turbulent developing pipe flow was solely determined by the growthof wall-boundary layers, the velocity profile peak in the pipe centerline may reach amaximum and then decrease again.Since turbulent developing pipe flow has been extensively studied and many reliableresults are available, this problem was chosen to test the proposed method. As in thelast chapter, a non-orthogonal rectangular type grid (see Figure 3.la) is used to replacethe ordinary cylindrical coordinate grid. Such non-orthogonal grids are very useful forsolving many internal flows.Based on the experimental results for this flow, the pipe length L was chosen as80 times the diameter D from the inlet, within which fully developed flow is generallyachieved. Flows with two different Reynolds numbers, namely Re = eUbD/ = 10 and3.0 x iO, are calculated, where Ub is the bulk velocity. At the inlet, a uniform velocity,a constant turbulence kinetic energy, and a constant dissipation rate are prescribed. kand e are evaluated from the following equations:k = l.5U2(u/ ), e = C,1c/l (3.21)where u/U is the turbulence intensity and 1 is the turbulence length scale. Here, the turbulence intensity u/U was taken as 0.02 and 0.04 for Re = iO and 3.0 x iO respectively,and the turbulence length scale was taken as one-tenth of the pipe diameter.For flow with Re = 1O, the computed velocity profile at the centerline, together withexperimental results from Richman et at (1973) are plotted in Figure 3.2a. A velocitypeak was found at about 29 pipe diameters from the entrance which differs from themeasurements. However, the overall agreement is reasonable if one notices that the largevelocity scale exaggerates the difference which is under 6%. Further, this phenomenonChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 66I [ttItlItIIIIutiltttItIIIIttitItIttIIIItt It’IlIItIIItIti°IItIIIIlIIttiI11ttt1tt11it25 tIIItItllhIItiIIttItIIIItIi111111111111111111II1IIII I0-0.50 0.00 0.50(b)Figure 3.1: (a) The non-orthogonal rectangular type grid. (b) Predicted velocity field atthe symmetry plane.(a)0.700 10 20 30 40 D 50(a)Figure 3.2: Predicted velocity profiles for Re = iOn, (a) developing velocity profile at thepipe center, (b) developed velocity profile at cross-section where z/D = 60.Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 671.201.00’0.90Present calculations0.80Measurements, J. W. Richman etc.60 70 80-0.50 -0.25 -0.00 yj 0.25(b)0.50E Measurements, J. W. Richman etc.Present calculationsFigure 3.3: Predicted velocity profiles for Re = 3.0 x iO, (a) developing velocity profileat the pipe center, (b) developed velocity profile at cross-section where z/D = 60.Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 681.201.10D1.000.900.800.700 10 20 30 40Z1D 50(a)60 70 8(b)Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 69was confirmed by other experimental results as discussed by Klein (1981). The flowbecomes fully developed after 50 pipe diameters. The profile of the predicted fully-developed velocity is plotted in Figure 3.2b with the experimental data. This figureshows that the fully-developed velocity profile is in good agreement with the experimentalresults.Similar computational results for the higher Reynolds number 3.0 x iO together withthe experimental measurements of Richman et al (1973) are plotted in Figure 3.3. Thevelocity peak occurs at about 32 pipe diameters; a little later than the case for the lowerReynolds number. Good agreement between the predicted and experimental results isalso seen here at large z/d.From both computations, accurate solutions were obtained for the flow in the well-developed region. This is because the k — e model and the ‘wall function’ describe thewell-developed turbulent pipe flows very well. Thus the method can be very accurate ifa suitable turbulence model is used.3.5.2 Developing Turbulent Flow in a 90° Curved Square DuctAs the second computation, turbulent flow in a 90°-bend square duct was simulated.Experimental measurements were reported by Taylor et al (1982). The geometry whichwas used in the experiment was also used for the computation, which consisted of a 90°curved square duct with a mean bend radius of 92 mm with an upstream straight duct0.3 m long and a downstream straight duct 0.2 m long. The dimensions of the crosssection are 40x 40 mm2. The Reynolds number defined by the bulk velocity, Ub, and thecross-sectional dimension is 40, 000.Exploiting the advantage of flow symmetry, only half of the physical domain was computed to save computational time. The configuration and grid used for the computationare illustrated in Figure 3.4. A grid with 21 x 11 x 81 nodes was generated, which hadChapter 3. Computation of Turbulent Flows Using Curvilinear GridsEEla,cD70upstream—3.3D=48 mmsymmetry(a) cross-section (b) vertical cross-sectionFigure 3.4: Illustration of computational grid for the curved duct.Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 71a higher density of grid points in the curved section and less grid points near the outletand inlet.As in the last example, a uniform velocity and constant k — € values were prescribed atthe inlet. Computational results are plotted for different cross-sections. The flow field inthe middle vertical plane is plotted in Figure 3.5. The maximum velocity occurs initiallynear the inner side wall in the upper portion of the curved section and gradually shifts tothe outside wall due to the radial pressure difference arising from the bend curvature. Asin the experimental results, strong secondary flows are found in the present computationwhich influence the whole flow pattern. The secondary current is developed graduallyand becomes strongest at about 9 = 60°, where the maximum secondary velocity reachesapproximately thirty percent of the bulk velocity. The secondary flow fields at twodifferent duct cross-sections together with the streamline velocity contours are plotted inFigure 3.6.The predicted velocity profiles at four different locations are compared in Figure 3.7with the measured results of Taylor et at (1982).The results agree well at the upstream portion of the curved section and deterioratewith increasing downstream distance. However, the overall agreement is still reasonable.The discrepancy in the downstream section is believed to be due to limitations of the k — €model for flows with streamline curvature and strong secondary flows. A higher-orderscheme and denser grid spacing may also improve the computational results.3.5.3 Developing Turbulent Flow in a 90° Curved PipeThe last example simulates developing turbulent flow in a 900 curved pipe. Measurementsof such a flow was reported by Enayet et al (1982). Laminar flow with the same geometrywas studied in the last chapter where a good agreement between the experimental andcomputational results was obtained. In the present study, the turbulent flow is computedChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 72(a)Figure 3.5:whole duct,Flow field in the middle vertical plane, Re = 40, 000, (a) flow field in the(b) flow field in the curved section.(b)CDiCDH/0Cl)IJ!J1/1’noJJ1I1//_\I1III”i1liii’’.’I— 0oVllliw\\0oQq(t I— oI-. CboCo0o(Il0 0CD C’)c. CD V\\\N\\Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids1.251.000.750.500.25-0.25 -0.00 0.25(b)Figure 3.7: Comparisons of velocity profiles with thelocations, (a) ,8 = 30°, (b) 3 = 60°, (c) ,8 = 750, (d) downstream xd = 0.25D.measurements at four different74S MeasurementCalculation(a) -o 0. 0-0.00(c)-0.00(d)Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 75Figure 3.8: Illustration of computational grid for the curved pipe.to test the present method. The geometry is similar to that in the previous exampleexcept that the cross-section is taken to be a circular disk. The diameter of the cross-section is 48 mm and the radius of the pipe curvature is equal to 2.8 times the pipediameter. The bend was fitted with a straight pipe 240 mm long upstream and 480 mmlong downstream.The grid and geometry are illustrated in Figure 3.8. A non-uniform and non-orthogonalgrid was used which was generated using an elliptical grid generation technique, to bedescribed in chapter 6. As was the case for the previous example, only half of the flowdomain was used for the computation.The computation was carried out for a Reynolds number of 43, 000; the same as thatused in the experiments of Enayet et al (1992). The results are plotted in a similarmanner as for the previous example. The velocity vectors at the middle vertical planezxGrid at Cross SectionChapter 3. Computation of Turbulent Flows Using Curvilinear Grids 76U/Ub=1 .0-Figure 3.9: Flow field in the middle vertical plane, Re = 40, 000.CfcjDCMCDCDnnoF-” 0P 0 CD ci, hi E ci, F-I 0 CM F-i CD ci, CD CD 0 0 0 0 ii 0 CM ci,ID Ci) a CD C)-o CCl) a m 2)-DoCD Cl) a CDD Cl) a CDa-Do Ca-o 0Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 78Figure 3.11: Comparisons of velocity profiles for different cross-sections, (a) upstreamO.5D, (b) 3 = 30°, (c) j3 = 60°, (d) 3 = 75°.(a) (b)(C) (d))Chapter 3. Computation of Turbulent Flows Using Curvilinear Grids 79are plotted in Figure 3.9. Unlike the pattern of laminar flow in the same geometry (seeFigure 2.11), the velocity maximum in the curved, upstream section occurs near theinner side and not near the outside wall as observed in laminar flow. As in the previousexample, the pressure-driven secondary flow extends over the entire flow region in thecurved and downstream sections. The secondary flow fields together with the streamlinevelocity contours at two different locations are plotted in Figure 3.10 for reference. Thepredicted velocity at different locations with the measured results of Enayet (1982) et alare plotted in Figure 3.11. The overall agreement is good, even though there are againsome discrepancies in the downstream portion of the curved section.3.6 ClosureNumerical simulation of turbulent flows in complex geometries using general curvilinearcoordinates is studied. The curvilinear coordinate-based method developed in the lastchapter is generalized to turbulent flows using the k — e two equation model with the ‘wallfunction’ treatment. The discretization of the k — e equations and the treatment of the‘wall function’ in curvilinear grids are presented. Methods enhancing the computationalstability are introduced.Several three-dimensional turbulent flows are computed using the developed solutionmethod and detailed results are presented. Computational tests show that the methodis efficient and reliable. Accurate solutions can be obtained for flows without curvatureand recirculation. The computational results for the last two examples show certaindiscrepancies with the experimentairesults. This is attributed to the limitations of the k—e turbulence model and the ‘wall function’ treatment for flows with strong secondary flowsand streamline curvature. However, the overall agreement observed for both examples isstill acceptable.Chapter 4Multigrid AccelerationIn this chapter, the development and performance of a multigrid procedure for thecomputation of three-dimensional laminar/turbulent flows using general curvilinear gridsare presented. Brandt’s (1981) Full Approximation Storage (FAS) in conjunction with theFull Multigrid cycling (FMG) is adopted. The numerical method for laminar/turbulentflows described in the last two chapters is used as the smoother. The present studyshows that the discrete governing equations can become inconsistent between fine andcoarse grids in curvilinear grids, as well as when the k — c model with the ‘wall function’is used for turbulent flows. Novel procedures are proposed to solve these problems.Several laminar/turbulent flows are calculated using the multigrid method to illustrateits performance.4.1 IntroductionThe convergence rate of traditional single-grid iterative solution methods seriously deteriorates as grids are refined. By using a Fourier analysis, Brandt (1980) showed thatiterative solution methods are only efficient in smoothing those error components whosewavelengths are comparable to the mesh size. The multigrid method is an efficientiterative solution procedure which exhibits convergence rates insensitive to grid refinement. The method has been widely used in the field of computational fluid dynamics.More recently, a few studies have been reported on the use of multigrid acceleration in80Chapter 4. Multigrid Acceleration 81curvilinear coordinates. Those works differ mostly in the choice of the basic numerical algorithms (coupled or decoupled), grid layout (staggered or collocated grids), andprimary velocity unknowns on non-orthogonal grids (Cartesian, contravariant velocitycomponents etc.). Rayner (1991) and Shyy et at (1993a,b) conducted multigrid calculations for two-dimensional flows in complex geometries based on decoupled sequentialsolution procedures. Those methods were based on the staggered grid system and thepressure-correction scheme. Rayner (1991) used the contravariant velocity componentsas the primary velocity unknowns, while Shyy et at (1993a) adopted a combination strategy in which both the Cartesian and contravariant velocity components were used tofacilitate the formulation of the mass and momentum equations. In addition to thedecoupled solution procedure, Joshi et at (1991) and Oosterlee et at (1992) reportedmultigrid algorithms for flows in two-dimensional complex geometries based on a coupled solution solver. These methods were based on the staggered grid system and thecontravariant-type velocity unknowns. In contrast to the above two-dimensional workswhich are based on the staggered grid system, Smith et at (1993) presented a multigridprocedure for three-dimensional flows on non-orthogonal collocated grids. The algorithmused by the authors solves the momentum equations for the Cartesian velocities as thedependent variables and stores all variables at the control volume centers.The implementation of multigrid algorithms for turbulent flows faces many additionaldifficulties, even in a Cartesian coordinate system. The k — e two-equation model with the‘wall function’ treatment is widely used for simulation of turbulent flows. The multigridacceleration for turbulent flows using such a two-equation turbulence model has beeninvestigated by a few authors. Based on experience with the difficulties of multigridcycling associated with the k—c equations and the ‘wall function’, Vanka (1987) developeda multigrid method which did not include the k — e equations within the multigrid cycle,but instead solved them on the fine grid and transferred only the turbulent viscosity. ThisChapter 4. Multigrid Acceleration 82method was also adopted by Yokota (1990) who reported good convergence rates. Eventhough this method can speed-up the convergence rate for turbulent flows, the multigridefficiency may be limited by the slow convergence of the k —6 equations on very fine grids.A more effective solution procedure can be expected if the set of governing equations asa whole are solved together using the multigrid cycle.Rubini et al (1992) presented a multigrid calculation for turbulent, variable densityflow over a three-dimensional, backward-facing step using the k — c model with ‘wallfunction’ treatment. The authors used a staggered grid system and the SIMPLE algorithm as the smoother. The convergence of the k — c equations was accelerated togetherwith the momentum and continuity equations using the multigrid strategy. The authorsfound that additional linearization of the source terms was necessary to prevent negative values of k and e. In addition, the negative values of k and e can arise due tonegative correction. This problem was circumvented by simply not updating values inlocations that would result in a negative value. Their results indicate a dramatic reduction in the computer time required to obtain a solution by a factor of over 25 for thefinest grid considered. Shyy et al (1993b) presented a multigrid solution procedure fortwo-dimensional turbulent flows in conjunction with the k — e two-equation turbulencemodel in general curvilinear grids. The method used a staggered grid arrangement andthe Cartesian velocity components as the primary variables in the momentum equations.They reported some special treatments which facilitate the successful implementation ofthe multigrid algorithm. Particularly, they noted the inconsistency of the ‘wall function’among different grid levels which may decrease the efficiency of the multigrid method.In the present study, the development of a multigrid procedure for computation ofthree-dimensional laminar/turbulent flows in general curvilinear grids is presented. Problems associated with multigrid calculation in curvilinear grids are analyzed. Novel techniques are proposed to solve these problems.Chapter 4. Multigrid Acceleration 83It was reported by Shyy et al (1993b) that the geometrical complexities can imposean upper limit on the number of effective grid levels useful for the MG method. Thepresent study shows that this limitation is caused by the inconsistency between the fine-and coarse-grid governing equations. When the coarse-grid defect governing equationsare directly formulated using the coarse grid nodes only, they will degrade relative to thefine-grid equations with the degradation of the spatial resolution of the coarse grid. Whenthe degraded coarse-grid equations are solved, the resulting correction will be inadequateto correct the fine-grid solutions thus reducing the efficiency of the multigrid algorithm.In the present study, this problem is analyzed and a new method is proposed to solve it.As mentioned in the previous paragraph, when the multigrid acceleration is appliedto all the governing equations, including the k — e equations, it is found that the formulation of the coarse-grid defect equation using the ‘wall function’ may cause inconsistencybetween the governing equations of fine and coarse grids. This problem is discussed inthis chapter and a novel practice is presented to allow the successful implementation ofthe multigrid method for turbulent flows using the k — e equations. Several techniqueswhich improve the multigrid solution efficiency for turbulent flows are introduced.Consistent restriction of residuals and solutions between coarse and fine grids is important for the nonlinear multigrid algorithm. It is impossible to preserve both the massand momentum conservations in the course of the transfer of solutions. However, massconservation can be maintained between grids by deliberately choosing the interpolatingformulation. In general, mass conservation between grids is usually not maintained inrestriction for methods using curvilinear grids. In the present study, a mass-preservingrestriction for the physical tangential velocity components is developed.In order to validate the multigrid performance, computations were carried out for several laminar/turbulent flows with significant non-orthogonal and strongly curved grids.Chapter 4. Multigrid Acceleration 84Comparisons were made between the multigrid and single-grid methods through measurement of the number of finest grid iterations, the equivalent work units and CPU time.The computational results show that the developed multigrid method is very efficient.4.2 The Multigrid AlgorithmThe idea behind using multigrid methods is to remove the low-frequency errors by solvingthe governing equations on various auxiliary coarse grids. For a multigrid method, it iscritical to use consistent governing equations between the coarse and fine grids. Asdiscussed later, this criterion may be violated in certain curvilinear grids. Therefore,additional care should be taken in the development of a multigrid method for generalcurvilinear grids.In the present study, the Full Approximation Storage (FAS) of Brandt (1981) inconjunction with the Full Multi-Grid (FMG) cycling is adopted. The solution procedureis initiated on the coarsest grid and subsequent initial solutions on finer grids are obtainedby interpolating converged solutions. On each finer grid and on the finest grid, the FASsolution procedure is applied to obtain a converged solution.To illustrate the problem to be solved, the FAS procedure is briefly described in thissection using a two-grid (c, f) system. The discrete governing equations on a fine gridare written as an operator equation as follows:L(q) = 0, (4.1)where L is the non-linear operator designating the discrete governing equations, q is thesolution vector which includes all dependent variables, and the superscript ‘f’ indicatesthat the solution and operator are associated with the fine grid. Suppose is anapproximate solution of equation (4.1). Unless ‘i is the exact solution, there will be aChapter 4. Multigrid Acceleration 85residual r such that= (4.2)Let c be a correction of the approximation solution ‘çb so that z/i + satisfies equation(4.1), namely,Lb1 + c) = 0. (4.3)Subtracting equation (4.3) from equation (4.2), the following equation is obtained:L(b + c1) — L(b) = —r’. (4.4)If the correction c can be found from the above equation, it is obvious that the exactsolution çif can be obtained. The idea of the multigrid method is to solve for the correction c, rather than solve for the solution on a coarse grid. In order to do so, the aboveequation is restricted to the coarse grid ‘c’,L’(Ib + cc) — Lc(I/) = —Ir, (4.5)where the superscript ‘c’ indicates that the quantities are associated with the coarse grid,while I is the restriction operator. For linear problems, the approximation solution Iz/.cancels out in equation (4.5) and a correction equation can be obtained. For non-linearproblems, the full approximation equationLC(q5C) = Lc(Ic?.,bf ) — (4.6)is solved in the FAS procedure and the correction cc is then obtained from the solutiondifference= c,jc —The auxiliary equation (4.6) containing the restricted fine-grid residual (or defect) isfrequently referred to as the defect equation in the literature. After the correction isChapter 4. Multigrid Acceleration 86solved from equation (4.6), it is prolongated to the fine grid to adjust the old solutionnamely, = ‘/‘f + Ibc — Ibf) where is the prolongation operator. If theadjusted solution q does not satisfy equation (4.1), the above procedure can be reapplied iteratively until the residue is below the required criterion. The above procedureforms a two-level FAS procedure.4.3 Consistent Formulation of Governing EquationsFrom the above multigrid procedure, it can be seen that the coarse-grid defect equation(4.6) used for the calculation of the correction is formulated from the fine-grid equation(4.1). Therefore, the defect equation should be a faithful representation of the fine-gridequation on the coarse grid. Three components, the residual r, approximation solutionbf, and solution operator L1 have to be transferred to the coarse grid to formulate thedefect equation. Attention has been given to the accurate restriction of the solution andthe residual in the literature. However, the restriction of solution operators also requirescareful consideration. The defect equation (4.6) is similar to the solution equation (4.1)except for additional source terms which include the residual restricted from the fine grid.Therefore, the coarse-grid operator LC is usually formulated by discretizing the governingequation on the coarse grid in the same manner as that used for the solution operator L.Such a practice is usually adequate for computations in Cartesian coordinates. However,this may result in inconsistent operators between coarse and fine grids for computationsusing curvilinear grids.Shyy et al (1993b) reported that geometrical complexities can impose an upper limiton the number of effective grid levels for the multigrid method. Their computationsdemonstrated that the efficiency of a multigrid method is not improved by using highergrid levels if the spatial resolution on coarser grids becomes too degraded. The presentChapter 4. Multigrid Acceleration 87(a) r1 (b)Figure 4.1: Grids in cylindrical coordinates: (a) 16 x 16 grid; (b) 2 x 2 grid.study indicates that this limitation is caused by the direct formulation of the defectoperator LC on the coarse grid which is no longer consistent with the fine-grid operator Vdue to the degradation of the coarse grid. On a curvilinear grid, a coarse grid operator LCdepends on the grid geometric quantities associated with the coarse grid. In general, theanalytic function of the coordinate transformation is unknown and the physical geometricquantities have to be calculated from the discrete grid nodes. Due to the deteriorationof the spatial resolution on the coarse grids, however, the physical geometric quantitiescalculated from the coarse-grid nodes can become inconsistent with those on the fine gridin the sense that they are not defined by the same coordinate functions.In order to illustrate this problem, a two-dimensional flow domain, namely, a curvedchannel is considered. For this particular problem, an analytic coordinate transformationexists:= rcos6,y = rsin&, r1 r r2,O , (4.7)which transforms the curved channel domain in the x— y plane to a rectangular domainin the r — 8 plane. The corresponding curvilinear grids are shown in Figures 4.1a & 4.lbfor a fine grid 16 x 16 and a coarse grid 2 x 2, respectively. From Figure 4.1, it can be seenjrChapter 4. Multigrid Acceleration 88Figure 4.2: Curvilinear grids: (a) 16 x 16 grid; (b) coarsest 2 x 2 grid; (c) 2 x 2 gridrestricted from the fine grid 16 x 16.that even the 2 x 2 coarse grid still fully represents the curved channel geometry. This isbecause the grid lines on the coarse grid are defined by the analytical function (4.7), whichis independent of the number of grid nodes. When the physical geometric quantities forfine and coarse grids are calculated from the coordinate transformation function (4.7),they are consistent since they are defined by the same coordinate system. In general,however, the physical geometric quantities have to be calculated from the discrete gridnodes When the physical geometric quantities are calculated from the discrete grid nodesusing the linear profile assumption, such as that used in chapter 1, the grid lines used forthe calculation are actually the line segments produced by connecting the neighbouringgrid nodes, as shown in Figures 4.2(a) &4.2(b) for the 16 x 16 fine grid and 2 x 2 coarsegrid, respectively. From these figures, it can be seen that the fine grid in Figure 4.2ais a reasonable approximation to the corresponding grid in the cylindrical coordinatesshown in Figure 4.la. Thus, the physical geometric quantities calculated from the fine-grid nodes are approximately equal to those defined by the coordinate transformation(4.7). However, the coarse grid in Figure 4.2b can no longer give a faithful representationof the corresponding grid in the cylindrical coordinate system shown in Figure 4.lb, andthus the physical geometric quantities calculated from the coarse-grid nodes are different(a) (b) (c)Chapter 4. Multigrid Acceleration 89from those calculated from the coordinate transformation (4.7). Therefore, the physicalgeometric quantities on the coarse and fine grids are inconsistent since they are essentiallyequivalent to those defined by two different coordinate systems.It should be noted that the grid lines on a coarse grid always coincide with thoseon the fine grid in an analytical coordinate system, such as the cylindrical coordinatesystem, no matter how coarse the grid is. Similarly, this argument should also be valid forgeneral curvilinear grids, namely, the grid lines on a coarse grid used for the computationof the physical geometric quantities should coincide with those on the fine grid in orderto be defined by the same coordinate transformation. For the curved channel geometry,the grid lines in a 2 x 2 grid defined by the equivalent transformation of the 16 x 16 finegrid in Figure 4.2a should be the restriction of the fine-grid lines as shown in Figure 4.2c,rather than those shown in Figure 4.2b. However, the physical geometric quantitiessuch as the surface areas (lengthes in 2D) and cell volumes defined in Figure 4.2c cannot be obtained from the coarse grid nodes alone. Calculation of the grid geometricquantities from the coarse-grid nodes is equivalent to using the low-resolution grid linesas shown in Figure 4.2b rather than the curvilinear grid lines as shown in Figure 4.2c.For instance, the length (area in 3D) of the cell-surface ‘s’ calculated from the grid nodesin Figure 4.2b is equal to the distance of the two neighbouring nodes AB rather than thelength of the corresponding curved segment AB in Figure 4.2c. Therefore, the physicalgeometric quantities calculated from the coarse-grid nodes only will be inconsistent withthose on the fine grid. When the defect equations on the coarse grid are formulated usingthe coarse-grid geometric quantities, an inconsistency will arise between the governingequations on the fine and coarse grids. The correction calculated from the inconsistentcoarse-grid defect equations will be inadequate to correct the fine-grid solution, thus,reducing the multigrid efficiency.Chapter 4. Multigrid Acceleration 904.4 Consistent Calculation of the Grid Geometric QuantitiesFrom the above discussion, we know that in order to obtain consistent governing equationsbetween fine and coarse grids, consistent physical geometric quantities have to be used.In the present study, a calculation procedure is developed to achieve this goal. Theidea of the method is to use the same grid resolution in the formulation of governingequations for all the grid levels. Instead of using the coarse-grid nodes to discretize theequations, the finest grid is used to formulate the discrete governing equations on allthe grid levels. All the coarse-grid geometric quantities required to formulate the defectgoverning equations are calculated from the finest grid instead of from the coarse-gridnodes. First, the grid geometric quantities on the finest grid, namely, the surface areas,cell-volumes and surface normal vectors, are calculated using the method presented inchapter 2. The geometric quantities in the subsequent coarser grids are calculated in thefollowing manner. The surface areas for a coarse grid are calculated by summing theareas of the relevant cell surfaces of the finest grid,A, (4.8)finestand the unit surface normal vector e2 on a coarse grid is calculated as the area-mean ofthe relevant unit surface normal vectors & in the finest grid,e = (4.9)C finestwhere the subscripts ‘c’ and ‘f’ indicate that the quantities are associated with the coarseand the finest grid respectively, and the summation is over all the cell surfaces of thefinest grid coincident with the coarse-grid surface. The cell volume of a coarse grid iscalculated as the summation of the relevant cell volumes of the finest grid,> Vf, (4.10)finestChapter 4. Multigrid Acceleration 91where the summation is over all the finest-grid cells coincident with the coarse-grid cell.Such a calculation procedure of physical geometric quantities is equivalent to computations using the coordinate transformation designated by the finest grid through all thegrid levels. For instance, the length of the cell-surface ‘s’ for the 2 x 2 coarse grid shownin Figure 4.2 will be the length of the curved segment AB in Figure 4.2c rather than theline segment AB in Figure 4.2b. With the present method, the coarse grids are actuallynot involved in the numerical formulations. Therefore, the deterioration of the geometricrepresentation of coarse grids does not influence the multigrid efficiency.4.5 Multigrid Acceleration of Turbulent FlowsThe problem discussed in the last section occurs for both laminar and turbulent flows. Inthis section, the problems particularly associated with multigrid acceleration of turbulentflows are discussed. As mentioned earlier, the k — c equations are solved together with themomentum and continuity equations in the present study using the multigrid procedure.The multigrid acceleration of the k — 6 equations faces additional difficulties which haveto be carefully addressed. First, the use of the ‘wall function’ to link the flow in thenear wall region to the main flow region may not be suitable for all grid levels. The‘wall function’ is valid only within a certain region near the wall. The calculation ofk and e near the wall relies on the assumption of the local equilibrium of productionof and generation of turbulence energy. In a six-level multigrid solution procedure, forinstance, the distance from the first scalar node to the wall on the coarsest grid willbe 2 = 32 times that of the finest grid. When the first scalar node of the finest gridis adequately located in the equilibrium region, the first grid node of the coarsest gridmay be too far from the wall to obtain a reasonable solution profile. Secondly, the useof the ‘wall function’ in different grid levels will cause inconsistent governing equationsChapter 4. Multigrid Acceleration 92between grids as will be demonstrated later. Thirdly, the strong nonlinearity and thedominance of the source terms in the k and € equations make these equations difficult tohandle in the multigrid cycling. The k — e equations have a strong nonlinear couplingthrough the production and dissipation in the source terms. These source terms areoften dominant in the discrete governing equations. Furthermore, unlike the velocitiesand pressure, certain k and c values updated by iteration or by adding corrections willcause complete failure of the algorithm. These problems are discussed in the followingsections and novel treatments are introduced to solve them.4.6 The Treatment of Wall Boundary on the Coarse GridsIn this section, the problem associated with the inconsistent governing equations causedby the use of the ‘wall function’ in different grid levels is considered. To illustrate thisproblem, a two-level multigrid procedure for a two-dimensional turbulent flow is considered. Figure 4.3 shows a two-level multi-grid near a wall boundary, where the big andsmall circles indicate the scalar nodes next to the wall on which the ‘wall function’ isapplied.. The k — e equations are applied to the cells marked with an ‘x’. If the ‘wallfunction’ is used to formulate the coarse-grid defect equation, the restricted residual irin equation (4.6) at the scalar node ‘E’, for instance, will be the average of the fine-gridresidual R at the nodes ‘A’, ‘B’, ‘C’ and ‘D’. However, the residuals at nodes ‘A’ and‘B’ are calculated from the k — 6 equations which are inconsistent with the operator LCat node ‘E’ formulated by the ‘wall function’.In equation (4.6), the restricted residual I includes the residuals of the k—€ equationson the fine grid, but the solution operator Lc has to be formulated using the ‘wall function’. Therefore, it is unlikely that a reasonable correction can be obtained by applyingthe ‘wall function’ in certain areas while the residuals in the fine grid are calculated fromChapter 4. Multigrid Acceleration 93x x x x xx x x x x xA B—)0 0 0 0 0 0C DF / F / F / F / F F / ft / FFigure 4.3: Illustration of using inconsistent defect equations at interior scalar nodes ‘x’.the k — c equations there. This will cause an inconsistency of the governing equationsbetween the fine and coarse grid levels.The deterioration of speed-up from the single-grid to the multigrid procedure forturbulent flow calculations using the k — 6 two-equation model with ‘wall function’ wasalso observed by Shyy et at (1993b). To solve this problem, they used a special gridcoarsening system in which the first grid lines away from the wall are retained duringthe course of grid restriction. With this practice, the treatment of wall boundary forturbulent flows is adequately handled at the coarser grid levels, resulting in an increasedspeed-up from the single-grid to the multigrid procedure. However, this non-standardgrid restriction will cause additional complexity in programming. Furthermore, the muchincreased mesh expansion ratios among the adjacent cells on the coarser grids is not wellhandled.In the present study, a new method is proposed to solve this problem. The idea isbased on the argument that consistent governing equations through all grid levels haveto be used. To avoid the inconsistency between the use of k — e equations and the use ofChapter 4. Multigrid Acceleration 94the ‘wall function’ treatment on different grid levels, the ‘wall function’ is only appliedon the finest grid. For the wall boundary on the coarser grids, a Dirichiet boundarycondition for the k — e equations is used in place of the ‘wall function’ treatment. Insteadof solving k and c from the balance between the local production and dissipation of Icin the equilibrium region near the wall, the restricted Ic and 6 from the finer grid areprescribed at the scalar nodes next to the wall as the boundary condition. A similartreatment to the above has been developed simultaneously and independently by Sun(1994).For the momentum equations, the wall boundary on the coarse grids can be treatedin the same way as for the Ic — e equations, that is, the velocity components next to thewall can be simply prescribed using the restricted velocity components from the fine grid.Such a treatment is observed to provide better performance than the standard treatment.However, this treatment will result in no correction for the tangential velocity componentson the cells next to the wall which limits the potential speed-up of the convergence rateof the multigrid method. In the present study, an alternative approach is used whichsolves the momentum equations in all the velocity positions without using the prescribedvalues from the fine grid. The tangential velocity components next to the wall are stronglydependent on the wall shear stress. If the wall shear stress can be determined, a governingequation for the tangential velocity components next to the wall can be formulated. Inthe ‘wall function’ treatment, the wall shear stress ‘r is calculated from the frictionvelocity UT which is defined by the ‘wall function’. Since the ‘wall function’ will notbe used on the coarser grids in the proposed method, the wall shear stress can only beobtained from the finest grid. Fortunately, it is possible to calculate the wall shear stresson a coarse grid from that on the finest grid. It is known that the wall shear stress isindependent of the grid level. Therefore, the correct wall shear stress on the coarser gridcan be calculated from values on the finest grid. In the present study, the wall shearChapter 4. Multigrid Acceleration 95stress on a coarser grid is taken as the area-weighted average of the corresponding wallshear stresses on the finest grid,1T,, = (4.11)C finestwhere the scripts ‘c’ and cf indicate the association between the coarse and finest gridsrespectively. The summation is over all the finest-grid wall surfaces coincident with thecoarse-grid surface. After the wall shear stress is obtained, the momentum equations forthe tangential velocity components next to the wall can be formulated by consideringthe wall shear stress as an auxiliary force in the momentum equation. The componentsresulting from the vector expansion of the wall shear stress r in the unit tangential vectorbasis {e} are used to modify the source terms in the discrete momentum equations (2.34-2.36). With such a treatment, the governing equations on the coarse grids can be solvedto provide accurate correction.4.7 The Treatment of k — e Quantities and Source TermsIn this section, several special treatments facilitating the multigrid calculation of turbulent flows are introduced. In chapter 3, a treatment of the source terms of the k —equations was presented for calculation in general curvilinear grids. A key feature of thetreatment is to assure positive source terms in both the discrete k and e equations. Thispractice is useful in avoiding negative values for turbulence energy and dissipation sinceit is known that any negative k and c acquired during the iteration will lead to completefailure of the solution method.In the coarse-grid defect equation, the source terms contain additional quantitiesfrom the restricted residuals and solutions which may become negative. Thus the valuesof k and c possibly acquire an erroneous negative value. A similar treatment to thatreported by Rubini et al (1992) is used to solve this problem. However, different fromChapter 4. Multigrid Acceleration 96their practice in which only the resultant error of the restricted solution was considered,the source term as a whole is considered. Suppose the coarse-grid defect equations on acoarser grid H for k and e can be written asa = > aq + Sc, (4.12)where S is the source term containing quantities from various sources, such as gridnon-orthogonality and restricted residual. If the source term 5C becomes negative, thefollowing discrete equation is used for the iterative solution solver:Sc(a—-,—)q = aflbcbflb. (4.13)‘1p nbThis practice will eliminate any negative source term in the k — e equations.In a multigrid method, the values of k and e may also become negative after addinga correction. Rubini et al (1992) and Shyy et at (1993b) adopted the same strategy bywhich they simply do not update any locations that would result in a negative value. Inthe present study, this strategy is modified by including an additional limitation. Sinceboth k and e are involved in the denominators of some terms in the governing equations,they should not be too small. For instance, if k becomes zero after adding the correction,the source term in e equation (4.12) will become infinite. This will also lead to the failureof the multigrid algorithm. Therefore, the correction scheme is modified for k and C inthe present study. Letq = + J!( -then for k and eJ q ifq>1O4&( otherwisewhere is the old approximation solution for Ic and e, qf — i/ is the correction on thecoarse grid. The present practice implies that if the Ic and e at some locations becomeChapter 4. Multigrid Acceleration 97unreasonably small by adding a correction, then the values at these locations will not beupdated.Due to the strong nonlinearity and the dominance of the source terms in both the kand c equations, under-relaxation of both the solution and source terms is adopted.4.8 Restriction and ProlongationIn addition to the consistent formulation of governing equations on different grid levels,the consistent restriction of solutions and residuals is also important for the multigridalgorithm. In a nonlinear problem, the discrete solution operator Lc on a coarse griddepends on the restricted velocity unknowns from the corresponding fine grid. Therefore, the restricted solution should be a good representation of the fine-grid solution.It is impossible to conserve both mass and momentum in the course of the restrictionof solutions. However, mass conservation can be maintained between grids by deliberately choosing the interpolating procedure. The mass preserving restriction has beenemphasized and adopted in several previous works, such as Sivaloganathan (1988).Mass conservation requires the mass flux through a control-volume surface on a coarsegrid to be equal to the total mass flux through the relevant surfaces of the fine grid,namely,.Flux = (u . ezA)c = >(eu . eA)f, (4.14)finewhere the subscripts ‘c’ and ‘f’ indicate the association with the coarse and fine grids,respectively, and the summation is over all the control-volume surfaces on the fine gridcoincident with the coarse-grid surface. Using equation (2.25) for the physical tangentialvelocity components, equation (4.14) can be rewritten as(eUiAiei. e). =(0UAe2 . ejf. (4.15)fineChapter 4. Multigrid Acceleration 98From the above equation, the physical tangential velocity component on the coarsegrid can be defined as follows:= ( e .eA) (eUCiAiei e) (4.16)2 CfjaFor the cell-centered variables and residuals the coarse-grid values are taken as thevolume-weighted mean of the eight corresponding fine-grid values. Taking the pressurep as an example, we have:Pc = >(V.p)f, (4.17)Cfinewhere the summation is over all the fine-grid control cells enclosed in the correspondingcontrol cell on the coarse grid. In order to preserve the high-frequency contents of the fine-grid residuals, an interpolation of higher order than that used for velocities is employed.The residuals of momentum equations on a coarse grid are taken as the area-weightedmean of the relevant twelve values on the fine grid. Taking the residual of the U —momentum equation as an example, the following formula is taken with reference toFigure 4.4:[Arf,f,kf_l + Arf_l,Jf,kf_l + Arf,f_l,kf_l+ ArIf_l,f_l,kf_l + 2Arf,f,kI + 2Ar’If,f_l,kf + 2Arf_l,f,kf+ 2Ar’f_l,f_l,kf + Arf,f_l,kf+l + ArIf_l,f_l,kf+l+ Ar f,jf,k+i + Ar1 if—1,jf,kf +1, (4.18)where the summation is over all the coefficients of the residuals. A bilinear interpolation is used as the prolongation operator. In all formulations for restriction and prolongation, the area-weighted average procedure is used for all quantities defined on thecontrol-volume surface, such as velocity unknowns, surface normal vectors, and momentum residuals. The volume-weighted average procedure is adopted for all quantitiesdefined at the control-volume center, such as pressure and mass residual.Chapter 4. Multigrid Acceleration 99kf+1Figure 4.4: Illustration of locations of — residuals on fine and coarse grids for restriction.11kf=2kc-1Chapter 4. Multigrid Acceleration 1004.9 Test CalculationsThe developed multigrid algorithm in general curvilinear grids has been implemented inthe CMGFD code. As shown in chapters 2 and 3, the code has been validated througha number of model problems for both laminar and turbulent flows and has been provento be a reliable solver. Since the objective of this chapter is to develop an efficientmultigrid algorithm, comparisons with experimental or other computational results arenot made. Instead, the detailed convergence history for various grid sizes and multigridlevels is presented. In order to provide a uniform measure of computational effort, a workunit, WU, is defined as being equivalent to the computational effort for one traditionalsingle-grid iteration on the finest grid. For instance, for three-dimensional problems,eight iterations on the next coarse grid are equivalent to one iteration on a fine grid. Theconvergence criterion is based on the absolute maximum residual of the four equations,i.e.,R = max2,3,k(R , Rv, Rn’, R”, Rk, R)where the maximum is taken over all computational cells.. A solution is considered to beconverged in the present study when the maximum residual is below 10—6. The reason forusing such a tight convergence criterion is that certain calculations require lower residualsthan others in order to show some flow phenomena. For instance, the third vortices inthe skewed cavity contain very small velocities which can not be seen without obtaininglower order residuals (see Demirdzic et al (1992)).The computations are started at the coarsest grid and the converged solution is prolongated to the next finer grid where the FAS procedure is applied. This procedure iscontinued until the finest grid is reached. A fixed V-cycle FAS procedure with two presmooth and two post-smooth iterations is used through all the computations. On thecoarsest correction level, ten iterations are applied to solve the defect equation.Chapter 4. Multigrid Acceleration 101Figure 4.5: Physical geometry of the skewed cube.4.9.1 Laminar Flow in a Skewed Cubic CavityIn the first application, the multigrid algorithm is applied to the laminar flow in a skewedcubic cavity. This example tests the performance of the developed multigrid method fornon-orthogonal grids. Flow in a cavity (or cube in 3D) has been considered as a standardtest case for numerical algorithms due to the characteristic nonlinear elliptic nature itshares with many flows of practical interest. In the present study, the laminar flow in askewed cubic cavity is considered. The counterpart of this problem in two dimensions,i.e., the flow in a skewed driven cavity, was proposed recently as a bench-mark problemby Demirdzic (1992) for testing numerical methods on non-orthogonal grids and has beenstudied by a number of authors. Figure 4.5 shows the problem currently being considered.The side wall is inclined at an angle /3 = 45° with the top wall moving at a velocity U.x= 1Chapter 4. Multigrid Acceleration10110°1 011 021 o-1 o1 o1 061 o5 10 15 20 25 30 35 40 45102(3)(2)(1)(2)(3)(4)17x17x17 grid33x33x33 grid65x65x65 grid97x97x97 grid c,)0.Cl)cD•(1)Work unitFigure 4.6: Convergence history of MG calculations using various grids.Chapter 4. Multigrid Acceleration 103Figure 4.7: Comparison of convergence behaviors of multigrid and single grid calculations.Chapter 4. Multigrid Acceleration 104All the walls are squares of dimension D x D where D = 1 in Figure 4.5. The Reynoldsnumber, defined as pUD/, was taken to be 100 in the present study.Multigrid calculations have been carried out for grids consisting of (1) 16 x 16 x 16,(2) 32 x 32 x 32, (3) 64 x 64 x 64, and (4) 96 x 96 x 96 cells using 4-level, 5-level, 6-leveland 6-level grids respectively. In these multigrid calculations the coarsest grid consistedof only 2 x 2 x 2 cells for the first three grids and 3 x 3 x 3 cells for the last grid on whichthe converged solution can be obtained with just a few iterations. All calculations wereinitiated from zero velocity and pressure fields. The convergence histories of the multigridcalculations for the four different grids are plotted in Figure 4.6 with reference to the workunit (WU). It can be seen that, unlike the performance of single-grid methods, the rateof convergence is independent of the grids and convergence is achieved for all the gridsin 40 work units. Vanka (1986b,1991) reported multigrid calculations for laminar flowsin a cubic cavity. His calculations showed that the residuals decreased by about threeorders of magnitude in roughly twenty finest-grid iterations. Considering the differencebetween the work unit and fine grid iteration, and different efforts used on one single-griditeration, the performance of the present multigrid calculation using a non-orthogonalgrid is found to be comparable with the results reported by Vanka (1991) for flow in acube where the Cartesian coordinates were used.In order to compare the performance of the multigrid and single-grid methods, computations with the traditional single-grid were also made. The convergence behaviorof the multigrid and the traditional single-grid methods are plotted in Figures 4.7 (a),(b), (c) and (d) with reference to the equivalent work units for grids of 17 x 17 x 17,33 x 33 x 33, 65 x 65 x 65 and 97 x 97 x 97, respectively. The convergence characteristics of these computations are also summarized in Table 4.1. CPU times for anIRIX(SGI) machine are reported. Single-grid computations for the fine grids 65 x 65 x 65and 97 x 97 x 97 were terminated after the first 2000 and 1000 iterations respectively asChapter 4. Multigrid Acceleration 105Finest MG MG CPU SG CPU CPU time Finest grid MG SGGrid Levels seconds seconds speed-up iterations WU WU17 x 17 x 17 4 47.2 213.9 4.5 18 34.91 24233 x 33 x 33 5 390.9 7356.2 18.8 20 35.71 109265 x 65 x 65 6 4843.5 3.8 x 105* 78* 21 37.46 4300*97 x 97 x 97 6 16495.2 1.6 x 106* 100* 24 41.46 6000*Table 4.1: Convergence characteristics of all calculations (the superscript * indicates thatthe values are estimated).too much computational time is required to bring the residuals below 106 as specifiedfor the multigrid method. However the number of iterations required for convergencecan be estimated from the logarithmic variation after a certain number of iterations asobserved in Figure 4.7. The estimated speed-up ratio was 78 and 100 for the fine grids65 x 65 x 65 and 97 x 97 x 97 respectively. From these results, it can be seen that themultigrid algorithm produced a substantial improvement in the convergence rate. Thespeed-up ratio was increased when the grid was refined. It should be noted that it isimpractical to seek a converged solution using the traditional single-grid method in theIRIX (SGI) machine for the fine grid (4) which would require months to finish. However,a solution can be obtained overnight on the same computer when the multigrid methodis used. This result clearly shows the necessity of using multigrid techniques for flowcalculations when very fine grids are required.4.9.2 Laminar Flow in a Duct with ObstructionThe previous example requires the use of significantly non-orthogonal grids. However,the problem does not involve other difficulties associated with grid curvature and thedegradation of spatial resolution on coarse grids as the grid is uniform and the gridgeometric quantities are constant everywhere. In the following example, the laminar flowChapter 4. Multigrid Acceleration 106in a duct with an arc on the bottom wall is computed for which a curved grid is required.The geometry shown in Figure 4.8 consists of a duct with a square cross-section modifiedby a circular arc on the bottom wall near the inlet. The counterpart of this problem intwo dimensions, i.e., flow in a planar channel with a circular arc on the bottom wall, wasrecently computed by Shyy et al(1993b) using a multigrid algorithm. Here, the maximumheight of the arc is taken as 30% of the inlet dimension which is higher than that usedby Shyy et al (1993b). Multigrid calculations were carried out for a fine grid consistingof 129 x 33 x 33 nodes using different grid levels in order to observe the efficiency of themultigrid method using higher grid levels. Figure 4.8 shows a coarser grid 65 x 17 x 17used in the present study where a non-uniform grid was adopted downstream. Figure 4.9shows the coarsest grid, 9 x 3 x 3, used in the present study, where Figure 4.9a wasobtained by the restriction of grid lines on the finest grid and Figure 4.9b was obtaineddirectly from the coarsest grid nodes. It can be seen that the grid shown in Figure 4.9bis not a faithful representation of the flow domain. When the defect equation (4.6) isdirectly formulated on the coarsest grid, it is equivalent to using the grid geometricquantities in Figure 4.9b. Therefore, the formulated governing equation is inconsistentwith that on the fine grid due to the degradation of the geometric representation of thecoarse grid. However, the physical geometric quantities calculated from the finest grid bythe proposed method are equivalent to those in the curvilinear grid shown in Figure 4.9a.Thus, the representation of the flow domain does not deteriorate on the coarse grids.Figure 4.10 shows the convergence histories of multigrid (MG) calculations using different grid levels for the grid of 129 x 33 x 33. From this plot, it can be seen that thecomputational efficiency is improved when a higher grid level is used even though thecoarsest grid was seriously degraded. Therefore, the effective grid levels useful for theMG method are not limited by the geometric complexity, unlike the results reported byChapter 4. Multigrid AccelerationFigure 4.8: Illustration of the grid 65 x 17 x 17 used in the study.107(a)(b)Figure 4.9: Illustration of the coarsest grid 5 x 3 x 3 in a xy— plane.Chapter 4. Multigrid Acceleration 10810210110010.110.2111106600Figure 4.10: Convergence history for computation of laminar flow in a duct using finestgrid 129 x 33 x 33.100 200 300 400 500Chapter 4. Multigrid Acceleration 109Shyy et at (1993b). The higher-level MG computations generated a substantial improvement in the convergence rate and a 5-level MG method resulted in a speed-up of theconvergence rate which was two times that of the 3-level MG calculation. This clearlyshows the advantage of the proposed method for multigrid calculations when the coarsegrid becomes degraded in spatial resolution.4.9.3 Laminar Flow in a Strongly Curved 90° PipeIn the last example, flow in a strongly curved pipe with a 900 bend was computed.Such flows have been studied experimentally and numerically by a number of authors(for example Smith et at (1993)) and contain important characteristics, including strongsecondary flow and strong grid curvature. The adopted geometry is similar to that usedin chapter 2 but with a more strongly curved bend of 90° with the radius of the bendcurvature equal to 1.5 times the pipe diameter (see Figure 4.11). The bend was fittedwith a straight pipe of two diameters length in the upstream direction and a sufficientlylong pipe downstream.Figure 4.11 shows the finest grid of 33 x 33 x 97 used for the present calculation. Anon-uniform grid in the downstream direction was adopted in order to place the outletboundary far enough into the uniform flow region. Computations were carried out usingthree different multigrid levels. For the 5-level multigrid calculation, the coarsest gridconsisted of 2 x 1 x 6 cells (as shown in Figure 4.12). As in the last example, Figure 4.12bwas obtained from the finest grid by restriction of the grid lines, and Figure 4.12a wasobtained directly from the the coarsest grid nodes. It can be seen that the geometricrepresentation of the coarsest grid in Figure 4.12a is seriously degraded.The convergence behavior for computations using 4 different grid levels is shown inFigure 4.13. Again, the present results show that the number of effective grid levels usefulfor the multigrid method is not limited by the geometrical complexity. The accelerationChapter 4. Multigrid Acceleration 110DownstreamEci,C’)c.DVertical planeCross planeFigure 4.11: Finest grid in the MG calculation of curved pipe flow, 33 x 33 x 97.Chapter 4. Multigrid Acceleration :iii(a) (b)Figure 4.12: Coarsest grid 3 x 2 x 7: (a) grid restricted from the finest grid; (b) gridobtained from the coarsest grid nodes.Cross plane Vertical plane Cross plane Vertical planeChapter 4. Multigrid Acceleration 11210210110°1 011 021 011 01 o-500 1000 1500 2000Figure 4.13: Convergence histories for computations using different grid levels.Chapter 4. Multigrid Acceleration 113ratio of the multigrid method increases when the number of multigrid levels is increased.With the 5-level multigrid calculation, the convergence rate accelerates more than 40times compared to the traditional single-grid method. However, when the number ofmultigrid levels is reduced to 4, the speed-up ratio is reduced to about 28. This resultshows the necessity of using higher grid levels for an efficient multigrid method. Thistype of flow was also studied by Smith et al (1993) using a multigrid algorithm and areduction of the overall CPU time up to 20% was observed compared with the single-gridmethod. The speed-up ratio of the multigrid method obtained using the 5-level grid inthe present study is found to be much higher than that reported by Smith et al (1993).4.9.4 Developing Turbulent Flow in A Strongly Curved 90° PipeThis example considers the developing turbulent flow in a strongly curved pipe. Thesame geometry as that used for laminar flow in the previous example was adopted.Computations were carried out for a Reynolds number of i0 with grids consisting of17 x 17 x 49 and 33 x 33 x 129 nodes using four and five-level multigrid cycles respectively.Figure 4.14 presents a typical plot of the convergence behavior for computations usingthe multigrid and single grid methods for the 17 x 17 x 49 grid. The fast convergenceof the multigrid method can be seen. A reduction of 85% in CPU time was achieved.The converged solution for the fine grid 33 x 33 x 129 was not carried out using thesingle-grid method due to the long computer time required, but the estimated speed-upratio of the multigrid acceleration based on the first 500 iterations of the single gridcalculation is more than 20 (or 95% reduction). This result is quite encouraging formultigrid calculations of turbulent flows.Chapter 4. Multigrid Acceleration 1141 o1021011001011 021111 0Figure 4.14: Convergence history for the computation of a developing turbulent flow ina strongly curved pipe.250 500Chapter 4. Multigrid Acceleration 1154.10 ClosureAn efficient multigrid procedure for three-dimensional laminar/turbulent flows in generalcurvilinear grids is presented based on the curvilinear coordinate-based method developedin the previous two chapters. It is shown that in certain cases, the ordinary formulation ofthe coarse-grid equations may result in inconsistent governing equations between fine andcoarse grids which reduces multigrid efficiency. Special procedures need to be taken forthe formulation of the coarse-grid equations in curvilinear grids and the treatment of the‘wall function’ in turbulent flows. Computational tests show that the multigrid algorithmdeveloped here efficiently reduces overall CPU time for both laminar and turbulent flows.On a very fine grid, the CPU time for the multigrid method was found to be about onehundred times less than that for a single-grid method. This clearly shows the necessity ofmultigrid method for computations using dense grids. Unlike the observations reportedin the literature, the number of effective grid levels useful for the multigrid method isnot limited by the geometrical complexities with the present method, even for seriouslydegraded coarse grids.Chapter 5Domain Segmentation and the CMGFD CodeIn this chapter, the methods developed in the previous three chapters are generalizedto non-structured curvilinear grids by using a domain segmentation technique (multiblock method). This technique divides the domain of interest into different sub-domains.Solutions are obtained by iteratively applying a solver to each sub-domain. The domainsegmentation strategy adopted here is described in this chapter. Particular attentionis given to the communication among neighbouring sub-domains. The performance ofthe multi-block method is investigated through several computational examples whichshow that the developed method has the capability to deal with complex configurationsusing multi-block curvilinear grids. A curvilinear, coordinate-based CFD code, calledCMGFD, is developed based on the methods described in this and previous chapters,and an existing CFD code called MGFD.5.1 IntroductionWith the methods developed in the previous chapters, laminar/turbulent flows in arbitrary single-block geometries can be solved with appropriate curvilinear grids. However,only structured, curvilinear grids can be used with those methods. It would be difficultto treat multiple flow regions using single block grids, such as that in film cooling of aturbine blade.One approach to deal with multiple flow regions with finite volume methods is theuse of a domain decomposition technique. This technique has become popular in dealing116Chapter 5. Domain Segmentation and the CMGFD Code 117with complex geometries in recent years. Leschziner and Dimitriadis (1989) presented acalculation procedure for flows in a domain consisting of a main duct and a side branchwhich may intersect at any angle. They used a multi-block, non-orthogonal grid torepresent domains in two conjunction ducts. At the duct interface, the computationalsub-domains overlap one set of finite volumes; those pertaining to the velocities normalto the interface. Perng et al (1991) and Hinatsu et al (1991) developed a multigriddomain-splitting technique for simulating incompressible flows in complex geometries.In this method, the momentum equations are solved independently on each subdomainwhile the pressure field is computed simultaneously on the entire flow field by a multigridmethod. However, they only applied the method to several two-dimensional flows inrelatively simple geometries.In this chapter, a domain segmentation technique in conjunction with the methodsdeveloped in the previous three chapters are developed to allow the use of block-structured(multi-block) curvilinear grids. Particular attention is given to the data transfers betweenneighbouring blocks. Several computational examples are carried out to validate themethod. The methods described in this and the previous three chapters are implementedin the curvilinear, coordinate-based CFD code CMGFD. The development of this codeis summarized.5.2 Domain Segmentation TechniqueThe idea of the domain segmentation technique is not complicated. The domain of interest is segmented into an arbitrary number of sub-domains, called computation blocks, andeach block is covered by its own Cartesian or curvilinear grid. A solver is then repeatlyapplied to each block to obtain a converged solution in the whole domain. However,the development of a generally-applicable 3D multi-block CFD code is not trivial. A 3DChapter 5. Domain Segmentation and the CMGFD Code 118multi-block and multigrid CFD code, called MGFD, was developed by Nowak (1992).The code is written generally which allows an arbitrary number of blocks to be used. Itcan be used to solve flows in arbitrary geometries which may be decomposed into a number of rectangular sub-domains. In the present study, the domain segmentation techniqueof Nowak (1992) is generalized to curvilinear grids by using the methods developed in theprevious chapters. This practice will allow the use of non-structured, curvilinear gridsto deal with arbitrary geometries with multiple flow regions. The developed multi-blockmethod is described in the following paragraphs.First, the domain of interest is segmented into a number of subdomains, which arecalled computational blocks. The segmentation of a given domain can be arbitrary anddepends on the configuration. Each block overlaps with its neighbouring blocks by twosets of cells to facilitate the communication among blocks.In each block, the single-block solution method developed in chapters 2 and 3 is applied to solve the governing equations. First, from either an initial-guess solution oran intermediate solution obtained from the last iteration, the auxiliary Cartesian velocity components and the turbulence viscosity are calculated. Next, the coefficients andsource terms in all the discrete governing equations are calculated using the formulations presented in chapter 2 and 3. The calculation of the coefficients and source termsrequires the geometric quantities and dependent variables at neighbouring control cells.This procedure needs special treatment at block boundaries. For boundaries of the blockwhich fall on flow boundaries, such as a wall or an inlet boundary, the actual boundaryconditions are applied. For boundaries which fall in a neighbouring block, the dependentvariables and the physical geometric quantities are obtained from the neighbouring block.The details will be discussed in the next section. After the discrete governing equationsare obtained, the alternative line-coupled iterative procedure described in chapter 3 canbe used to update the dependent variables.Chapter 5. Domain Segmentation and the CMGFD Code 119The above procedure is applied to all the blocks sequentially. After sweeping allthe blocks once, the turbulence viscosity and auxiliary Cartesian velocity componentsare recalculated for each block and all the coefficients and source terms are recalculatedbased on the updated turbulence viscosity and velocity unknowns. The entire procedureis cycled through all the blocks until the residuals become sufficiently small.Since the unknowns in only one block are solved at one time, the storage required forthe quantities in the discrete equations is much less than that required when the wholeflow domain is solved together. The number of cells contained in each block is usuallysmall and is independent of the number of cells used in the whole flow domain, so thecomputer storage required can be reduced substantially. More importantly, this solutionprocedure is able to deal with flows in multiple flow regions.5.3 Communication Between BlocksEfficient communication between the blocks has to be established in order to obtainaccurate solutions and fast convergence. This is done by appropriately transferring databetween adjacent blocks in each iteration. There are different methods to transfer databetween neighbouring blocks. In the present study, the method used by Perng (1991)and Nowak (1992) is followed. In this method, each computational block is extendedinto its neighbouring blocks by two sets of cells. The discrete governing equations on theextended cells are formulated using the physical geometrical quantities and dependentvariables obtained from the neighbouring blocks. The discrete governing equations onthe extended domains are solved in the same way as those on ordinary cells by using aline-coupled solver.In the CMGFD code, three types of multi-block curvilinear grids are used whichrequire different treatments: (1) the grid lines in adjacent blocks match at the interfaceChapter 5. Domain Segmentation and the CMGFD Code 120with complete continuity, (2) the grid lines match at the interface but with discontinuousgrid line slopes, and (3) the grid lines do not match at interface. For the first type ofgrid, an interior boundary condition is used which ensures that both mass and momentumconservation are satisfied at the interface. This interior boundary condition in each blockis implemented by using the overlapping cells in each of the neighbouring blocks. Forcurvilinear grids, the formulation of governing equations on the extended cells in a blockrequires the grid physical geometric quantities in its neighbouring blocks. The transferof the grid geometric quantities from one block to another is carried out at the beginningof the solution procedure. The unknowns on the extended boundaries are prescribed byusing the intermediate solution of the associated neighbouring blocks. When the coupledGauss-Seidel solution solver has swept over all the computational blocks, the unknowns inthe whole flow domain, including all the interfaces, are updated. The converged solutionsatisfies both mass and momentum conservation near the interfaces. Therefore, such anapproach will not produce any additional inaccuracy near the interface.For the second type of grid, since the grid lines change directions abruptly across theinterfaces (see in Figure 5.1), the derivatives of the grid coordinates are not continuousacross the interface. Therefore, the ordinary numerical formulations using derivativesof grid coordinates may produce significant numerical errors. With the numerical formulations described in chapter 2, the coordinate derivatives are avoided. However, interpolations are required to obtain the physical geometrical quantities and the velocitycomponents at the positions where they are not defined. These interpolations near theinterface are important as they may cause significant numerical error. This is discussedas follows.For convenience of discussion, the scheme is presented for the two-dimensional case butit can be generalized to three-dimensional problems without additional difficulties. Figure5.1 shows a two-block grid with discontinuous grid line slopes at the block interface. DueChapter 5. Domain Segmentation and the CMGFD Code 121Figure 5.1: Illustration of multi-block curvilinear grids with discontinuous grid line slopesacross the interface.11i-i 1-112 I 1+1/2Chapter 5. Domain Segmentation and the CMGFD Code 122to the change of grid line slopes, the surface normal vectors e1 (or contravariant vectora1) change direction abruptly across the interface. From equation (2.25) we know thatthe velocity component UE is grid oriented and depends on the surface normal vectore1. Therefore, the change of grid line direction can produce a large jump of UE acrossthe interface even with constant a velocity field, and the ordinary interpolation becomesinadequate.The numerical formulation of the momentum equations needs two (three for 3D)velocity components at each velocity position. Since only one velocity component isdefined at a velocity position (control volume surface) in a staggered grid arrangement,other components have to be obtained by interpolation. In particular, interpolation isneeded to obtain the velocity UE at each U’l position on the interface. As an example,we consider the interpolation procedure at the point (i, j+1/2) on the interface. Here,the UE_ velocity positions are marked with the (i-1/2, j), (i+1/2, j), etc., and the U’—velocity positions are marked with (i, j-1/2), (i, j+1/2), etc. On the point (i, j+1/2), onlyU” is defined and Ue has to be obtained from its neighbouring values by interpolation.One natural interpolation scheme is to take the average of the four neighbouring values,i.e.,UE(i,j+)= 2,j 2,3+1+2,1+1) (5.1)As discussed in the previous paragraph, the values of U(i — 1/2,j) and U(i + 1/2,j)in segment 1 can be significantly different from the values at the neighbouring pointsUE(i— 1/2,j + 1) and UC(i + 1/2,j + 1) in segment 2 due to the sudden change of thesurface normal vectors e’ across the interface. Averaging these sometimes significantlylarge different values unavoidablely produces certain numerical errors. Furthermore, thevelocity unknowns in scheme (5.1) have inconsistent base vectors. As discussed in chapter2, each velocity component U pertains to a base vector. The surface normal vector e’ isChapter 5. Domain Segmentation and the CMCFD Code 123used in the present thesis as the base vector to define U. A base vector for the averagedvelocity unknown UE(i,j + 1/2) has to be properly defined by interpolation. However,a proper interpolation scheme to define the surface normal vector e1(i,j + 1/2) at theinterface can not be justified in view of the suddenly changed surface normal vectors e1across the interface.In the present study, an interpolation scheme is proposed to avoid these uncertaintiesof interpolation near the interface. Instead of the full average scheme (5.1), a one-sideinterpolation is used, namely, the interpolation at the interface is only taken from oneside of the interface. For example, the velocity component U(i, j + 1/2) at the interfaceis averaged only from one segment, either from segment 1,1 Ui — • + U1i+ 1U(i,j + ) = 2’1 2 2’s) (5.2)or from segment 2,U(i,j+ ) = U(i— ,j+1)+U(i+ ,j+1) (53)If the interpolation is taken in segment 1, for instance, using equation (5.2), the surfacenormal vector e’ will be calculated in the same way by averaging the two surface normalvectors in segment 1 without using surface normal vectors from segment 2. Even thoughthe scheme (5.2) has only first order of accuracy, the base vectors (surface normal vectors) for the three velocity components U(i,j + 1/2), U(i — 1/2,j) and U(i + 1/2,j)are consistent since there is no significant difference among them. This avoids the inconsistency of the base vectors (the surface normal vectors) associated with the velocitycomponents used in the full average scheme (5.1).The choice of scheme (5.2) or (5.3) can follow the idea of the upwind differencescheme. If UE(i,j + 1/2) > 0, the fluid flows from segment ito segment 2. The velocitycomponents U at point (i,j+1/i) are more likely to be influenced by those in segmentChapter 5. Domain Segmentation and the CMGFD Code 1241. Therefore, the interpolation scheme (5.2) based on segment 1 is used. Otherwise thescheme (5.3) based on segment 2 is adopted.For the third type of grid, since the grid lines do not meet at the interface theextended cells can not be obtained directly from their neighbouring blocks. To reducethe inaccuracy caused by interpolation, only semi-matched grids at the interface areconsidered, namely grid lines in one block which meet with a coarser or finer grid ofanother block (the same concepts as those used in a multigrid method). An illustrationof this type of grid is given in Figure 5.2. In this figure, the grid in Segment 1 (bottom)is denser than that in Segment 2. For such grids, the extended cells for a block can beconstructed using a finer or coarser grid from its neighbouring blocks as illustrated inFigure 5.2. From Figure 5.2 it can be seen that the extended cells of Segment 1 are just thefiner grid of Segment 2. The grid physical geometric quantities and the unknowns in theadditional cells are obtained from the corresponding values in the neighbouring blocks byrestriction or prolongation with the same strategy as that used in the multigrid methoddescribed in chapter 4 (see Section 4.4). Such an assumption reduces the computationaleffort significantly for non-matched curvilinear grids.5.4 The CMGFD CodeAs mentioned earlier, the CMGFD code was developed based on an existing CFD code,MGFD, at UBC. The code was developed in four steps. In each step, detailed investigation of the method and validation of the code were carried out in order to achievean efficient and accurate code. The methods described in chapters 2-5 have been implemented in the code. This is discussed in the following paragraphs.Unlike Cartesian coordinate-based calculations, the computational grid for a complexChapter 5. Domain Segmentation and the CMGFD Code 12511Segment 2-Extended cellsU —irom Segment 1InterfaceExtended cellsfrom Segment 2Segment 1Figure 5.2: Illustration of the semi-matched multi-block curvilinear grids.Chapter 5. Domain Segmentation and the CMGFD Code 126geometry is usually not simple to generate. For generality, the code assumes that a curvilinear grid is available for a given problem and leaves the grid generation task for variousgrid generation techniques, such as one which is described in chapter 6. Therefore, theprogram GEO developed for generating various non-uniform Cartesian grids in MGFDcode actually will not be used in the CMGFD code.For a given grid, the Cartesian coordinates of all the grid nodes are loaded into theCMGFD code. A subroutine package has been developed to deal with the curvilineargrids and arrange the basic block structures. The grid geometric quantities at cells in eachblock are calculated from the given grid nodes. First, the surface normal vectors, surfaceareas and volumes are calculated using the grid nodes. The unit tangent vectors arethen calculated from the surface normal vectors by the orthogonal relations. Secondly,the physical geometric quantities are extended two sets of cells outside the block tofacilitate the formulation of the governing equations. For a boundary coinciding with aflow boundary, the physical geometric quantities outside the boundary are taken to beequal to the nearest values. For a block boundary falling in a neighbouring block, thegeometric quantities are directly taken from this neighbouring block.There are two different approaches for storing the geometric quantities. One approachis calculating only one set of basic geometric quantities, for instance the surface areavectors, and storing them in the computer memory. All other quantities used in theformulation of governing equations are calculated from these basic geometric quantitiesin each iteration. The advantage of such an approach is that substantial computermemory can be saved. However, recalculation of these quantities in each iteration requiressignificant computer time. In the present study, an alternative approach is followed andall the geometric quantities are calculated at the beginning and fixed throughout thesolution procedure.As a by-product, the grid is checked using the geometric quantities, un-usable gridsChapter 5. Domain Segmentation and the CMGFD Code 1.27such as the grids with zero or negative volume cells, are reported to the user. This oftenhappens due to the wrong output of the grid generation code or inconsistent geometricinformation loaded to the CMGFD code. In addition, there are a number of physicalcharacteristics that affect the quality of the grid, namely orthogonality, smoothness, cellexpansion ratios, etc. All these quantities are checked to identify low quality grids forpossible improvement.The computed geometric quantities are used to formulate the governing equations.This was done by properly modifying the MGFD code. All the code involved with gridinformation (such as cell dimensions dx, dy, dz) are modified throughout the code. Theformulations for the coefficients and source terms are substantially changed from thosein the MGFD code. The source terms in the discrete momentum and continuity equations are particularly complex in the present formulation as they include non-orthogonalpressure, non-orthogonal diffusion, and grid curvature terms. Their calculation requiresvarious interpolations of velocity, pressure and geometric quantities. It is found that theinterpolation procedures are important in terms of computational accuracy and efficiency.An interpolating procedure was developed through testing a number of problems.The multigrid method proposed in chapter 4 is implemented in the CMGFD code byusing a new subroutine package. This package contains three subroutines namely, ‘FMG’,‘FAS’ and ‘GEOM’. The FMG and FAS subroutines execute the standard FMG-FASprocedure. They are written generally which allows easy implementation of multigridmethods to other single-grid CFD codes. This package will also be used in the multigridelliptic grid generation code, MBEGG (see the next chapter). The GEOM subroutine iswritten to calculate the physical geometric quantities on all the coarse grids using themethod developed in chapter 4.So far, the code is limited to 3D steady, incompressible laminar/turbulent flows withscalar transfer. For turbulent flows, the standard k — e model with the wall functionChapter 5. Domain Segmentation and the CMGFD Code 128treatment is used for the turbulence closure. The CMGFD code has the capability todeal with arbitrary 2D/3D geometries by using separate curvilinear grids in different flowregions. There is no restriction for curvilinear grids to be used and no limitation for thenumber of segments in each flow region. An efficient multigrid method is implementedin the code and promises fast convergence for both laminar and turbulent flows.5.5 Test ComputationsThe validation of the CMGFD code has been carried out throughout previous chapters forvarious geometries using single-block curvilinear grids. In this section, several calculationsusing multi-block curvilinear grids are presented.5.5.1 Computations of a L-Shaped DuctIn the first example, laminar flow in a b-shaped duct is calculated. This example allowsone to validate the solution method using multi-block curvilinear grids with discontinuousgrid slopes across block interfaces. In this study, two computational grids are adopted,one is a two-segment non-orthogonal grid with discontinuous grid slopes and the otheris a two-block Cartesian type grid. Figure 5.3 shows the two-block non-orthogonal gridused for the present study. It can be seen that the grid lines change direction significantly (by 900) at the interface. For comparison, computation is also carried out usinga two-block Cartesian grid as shown in Figure 5.4. Figure 5.5 shows the predicted flowfield in the symmetry plane of the duct using the two-block non-orthogonal grid. Tworecirculation zones are found in the two corner regions. Figure 5.6 shows the predictedflow field using the orthogonal grid at the same plane. From Figures 5.5 and 5.6 the similarity of results predicted by the non-orthogonal and Cartesian grids can be seen. ForChapter 5. Domain Segmentation and the CMGFD Code4-Cci)EC)ci)(I)Two-segment grid, segmenti: 33x1 7x65, segment2: 33x1 7x65129Figure 5.3: Illustration of the adopted two-block non-orthogonal grid at the symmetryplane for the L-shaped duct.2InterfaceCD0C 00Cl)0_HHCD i-i CDI-.0 CDCD CoN 0CD •CD F’.)CDCD i-I CD Co CD 0Chapter 5. Domain Segmentation and the CMCFD Code 131Calculation using two-segment curvilinear gridFigure 5.5: Predicted velocity field at the symmetry plane by the non-orthogonal grid.Chapter 5. Domain Segmentation and the CMCFD Code 132Calculation using two-segment cartesian coordinatesFigure 5.6: Predicted velocity field at the symmetry plane by the orthogonal grid.Chapter 5. Domain Segmentation and the CMGFD Code 1331.00.0 0.2 0.4 0.6 0.8Figure 5.7: Comparison of the u-velocity profiles at the interface shown in Figure 5.3.Chapter 5. Domain Segmentation and the CMGFD Code 134Figure 5.8: Comparison of maximum residuals by multigrid and single-grid calculationsfor the non-orthogonal grid.10110010110211 o110.6500 1000 1500 2000Chapter 5. Domain Segmentation and the CMGFD Code 135a quantitative comparison, the predicted U-velocity (U is the Cartesian velocity component in the outflow direction) profiles at the centerline of the interface from the twogrids are compared in Figure 5.7. Good agreement can be seen. These results show thatthe developed method allows the use of multiblock grids with significantly discontinuousgrid line slopes. Computations are carried out for the non-orthogonal grids with bothmultigrid and single grid methods. The convergence behaviour of multigrid and singlegrid calculations are shown in Figure 5.8. It can be seen that the multigrid methodconverged much faster than the single grid method. This shows the advantage of themultigrid method for computations using multi-block curvilinear grids.5.5.2 Calculation of A HydrocycloneHydrocyclones have many applications in industry. Figure 5.9 illustrates a conventionalcone-cylindrical hydrocyclone. The cylindrical part is closed at the top by a cover,through. which the vortex finder protrudes some distance into the cyclone body. Theunderfiow leaves through the opening in the apex of the cone. The liquid enters thehydrocyclone through a tangential inlet (feed opening) and leaves through the vortexfinder and underfiow orifice.The developed CMGFD code is applied to simulate fluid flow in a hydrocyclone. Except for the region in and just around the tangential inlet duct, the motion of the fluidwithin the hydrocyclone has axial symmetry. Therefore, the flow can be treated approximately as two-dimensional, axi-symmetrical, while the angular velocity component istreated as a scalar quantity according to a discussion with Abdullah (1993). Figure 5.lOashows the non-orthogonal grid used for the present calculation. The grid contains threeblocks and the third block grid is generated by solving a two-dimensional grid generationsystem. The predicted velocity field in a vertical plane is shown in Figure 5.lOb. Fromthis figure, the separation of fluid can be seen. Most of the incoming fluid moves in anChapter 5. Domain Segmentation and the CMGFD CodeI=O.4*DDO=O.34*D136D4.Vortex finder.I DOD=O.075mz=L=5*DTangential vel. inlet, feed openingz=2.9*Dxunderfiow orifice (apex opening)Figure 5.9: Illustration of a conventional hydrocyclone.Chapter 5. Domain Segmentation and the CMGFD Code 137(b)Figure 5.10: (a) 3-block mesh used in the present study, (b) predicted velocity field.outer helical flow into the outer portion of the inverted cone where it begins to feed acrosstowards the center. Some of the downward flow leaves through the underfiow orifice inthe apex of the cone while the rest reverses its vertical direction and goes up via theinner helical flow and out through the vortex finder.5.5.3 Laminar Flows Over A Turbine BladeThe last example in this chapter is a computation for laminar flows over a turbine blade.Such a flow has no direct application to turbine engines since the flow in film coolingof turbine blades is actually turbulent. Nevertheless, the present calculation is useful inproviding guidelines for calculation of flows in complex configurations using multi-blockcurvilinear grids as well as the prediction of film cooling of an actual turbine blade.The computational domain is chosen according to the UBC experimental turbineblade model of Salcudean et al (1993). This model uses four rows of circular injectionblock3: 9x1 1 (a)Chapter 5. Domain Segmentation and the CMGFD Code 138Figure 5.11: Computational domain for a turbine blade model.tubes located at a semi-circular nose coupled to a flat after-body. Exploiting the advantage of its symmetric property about the stagnation line, only half of the flow domain istaken as the computational domain. Figure 5.11 shows the physical domain selected forthe present numerical study where three injection tubes are adopted in a row of circularcooling holes in order to simulate flows with multi-injection holes. The physical domainis segmented into four subregions, a main flow region and three injection coolant regions.In order to obtain a continuous grid at the interface, a rectangular type grid is used forthe injection holes instead of cylindrical grids. A multi-block curvilinear grid has beengenerated using the MBEGG code described in the next chapter. This grid is shown inFigure 5.12.Computations were carried out for a Reynolds number of 1000 and mass flow ratio of0.5, where the Reynolds number is defined with respect to the radius of the semi-circularMain gaCoolant holesChapter 5. Domain Segmentation and the CMGFD Code 139Figure 5.12: A 4-block curvilinear grid generated for the domain shown in Figure 5.11.CCd0U-,if0ci)C.)PNx.4-.I0C)ci)>>4-..00ci)>Ct’ci)0C0.4-C.)ci)CICt’ci)CCl)0.4-.C)ci)>>.4-00ci)>Ct’C)0-J.a,I0U,ifC00if9CU,if900U)000009U,090U,if9a,C00a,0a,-da,0a,-da,00a,a,Hka,bOoo00U)0U)Cif00C00U)CU)0U)C’J00ffd000Chapter 5. Domain Segmentation and the CMGFD Code 14].xFigure 5.14: The velocity at a surface distanced from the turbine wall by two cellsCD i-iCDi—Q”•CD-“O CDOacCDc3I— 0o0.CDCDjCDQ_•p,CD÷CD0(00CDI—hC0o 0I-I-s. 0O-I-:D-CD000) 0Chapter 5. Domain Segmentation and the CMGFD Code 143body and main-stream velocity. Figure 5.13 shows the calculated velocity vectors on axz-plane through the middle of an injection hole, from which it can be seen that thevelocity is very non-uniform at the coolant hole exits. The flow in the injection holeexits is affected by the main stream velocity and strong pressure gradient and, thus,changed rapidly with respect to magnitude, direction, etc. The velocities at a surfacedistanced from the turbine wall by two cells are plotted in Figure 5.14. Figure 5.15shows the development of the non-uniform profile near the exit of injection holes. Thestreamline velocity contour is plotted at the cross-planes of an injection hole at fourdifferent locations. A strongly non-uniform velocity profile near the injection hole canbe clearly seen, which shows that the use of a uniform velocity profile as the boundarycondition at the injection exit, used in many previous studies, could be very inaccurate.5.6 ClosureA domain segmentation technique in conjunction with the curvilinear coordinate-basedmethod developed in the previous three chapters is described. Communication betweendifferent blocks is emphasized. The developed multi-block method allows the treatmentof arbitrary complex configurations using non-structured curvilinear grids. A multiblock and multigrid code is developed based on an existing CFD code. The multi-blocktechnique in conjunction with the methods developed in the previous three chapters isimplemented in the code. Computational examples show that the developed method andcode can be used to solve flows in complex geometries using block-structured curvilineargrids with significant discontinuities in grid line slopes.Chapter 6Numerical Grid GenerationIn this chapter, the development of a 3D elliptic grid generation code is presented. Thecode is developed to address the need to generate appropriate grids for analyzing the filmcooling of turbine blades. An elliptic grid generation method is used in the code, whichsolves three nonlinear, elliptic, grid generation equations. A multigrid method is developed to solve these equations which provides an efficient grid generator. Techniques tocontrol grid quantities are introduced. A multi-block grid generation strategy is adoptedto allow the generation of block-structured curvilinear grids in arbitrary geometries.6.1 IntroductionIn order to simulate the film cooling process over a turbine blade, suitable grids haveto be generated to represent the actual blade geometry. One of the important tasksin the present thesis is to develop an efficient grid generator to generate suitable gridsover turbine blades. This includes the development of a multigrid and multi-block gridgeneration code, called MBEGG.Grid generation has been investigated extensively in recent years and a large numberof research papers have been published on the generation of body-fitted grids for a widevariety of geometries in two as well as three dimensions. A review paper by Thompson(1982) summarized most of the work done before 1981. The grid generation methodsused in the literature can be classified into three categories: (1) conformal mapping,(2) algebraic grid generation methods and, (3) elliptic grid generation methods. Among144Chapter 6. Numerical Grid Generation 145these methods, the third type of method has the best potential for the generation ofgrids in complex three-dimensional configurations. In the present study, the elliptic gridgeneration method proposed by Thompson (1985) is adopted. The method is based onthree elliptic partial differential equations. When proper grid nodes are specified at theboundary of a domain, a smooth grid can be generated by solving these equations.In the grid generation equations, the inhomogeneous terms are problem-dependentfunctions (which are commonly called control functions). The choice of an appropriatescheme for the evaluation of the control functions is one of the important considerationsin an elliptic grid generation method. Methods to achieve certain grid quantities throughthe control functions, such as stretching and orthogonality, have been implemented inthe code. When these grid quantities are initially specified by users, the code can choosethe appropriate source terms to obtain the desired grids automatically.In an elliptic grid generation method, a suitable solution method has to be developedto solve the three nonlinear, elliptic, differential equations. The multigrid method isan ideal solver for the grid generation equations due to the fully elliptic nature of theseequations. In the present study, a multigrid method is developed to solve these equations.It is found that the method produces a much faster convergence rate than the traditionalsingle grid method.For complex configurations, it becomes necessary to partition the physical domaininto a number of subregions. Grids are generated in the different subregions. A numberof papers have been reported in the literature for generating multi-block (composite)grids in various complex configurations. A detailed discussion can be found in the reviewpaper of Thompson et at (1982). In the MBEGG code, a multi-block method is used toallow the treatment of multi-region geometries. The development of the MBEGG codetogether with the methods adopted is described in the following sections.Chapter 6. Numerical Grid Generation 1466.2 Elliptic Grid GenerationThe elliptic grid generation method proposed by Thompson (1985) is based on the following well-known Possion equations:\7k = Pk(xl,x2,3), Ic = 1,2,3, (6.1)where (x1,x, x3) are the Cartesian coordinates in a physical domain and (i=1, 2, 3)are the curvilinear coordinates in a regular domain. Equation (6.1) has the same form asthe heat conduction equation, where P2(xi, x, x) are heat source terms. The equationsestablish a correspondence between the Cartesian coordinates (x1,x2,x3) in a complexphysical domain and the curvilinear coordinates 2, in a cuboid domain. However,the equations in this form, yielding (i,2,3) for given (x1,23)would be difficult touse, especially when specifying the boundaries. Therefore, the following transformedPoisson equations are usually used in the actual grid generation procedure:O2Xk ä2Xk ö2Xk____a11 2 + a22 2 + a33 2 +2a12ö + 2a138 + 2a381 2 163 23+ + + = 0, (6.2)where the coefficients a23 are given by the following expression:a23 = bm:bmjm1and where b23 is the co-factor of the (i, j) element in the transformation matrixM = ()3x3with J=det(M) as the Jacobian determinant. In the above equations, the P2 are problemdependent functions which can be used to control the grid characteristics, such as orthogonality and grid stretching. Thus, an important part of the work in the elliptic gridChapter 6. Numerical Grid Generation 147generation method is to develop an appropriate scheme for the evaluation of the controlfunctions. When a grid is specified at the boundary of a domain, the smoothest possiblegrid can be obtained by solving the above grid generation equations.6.3 Control FunctionsTo obtain the desired grids, proper control functions have to be constructed. In theMBEGG code, two types of control functions are calculated from the boundary gridnodes which may be specified by the user. The first type of control function controls thegrid stretching, and the second controls the grid orthogonality near the boundaries. Suchcontrol functions are implemented in the code to allow the users to generate the requiredgrids. The two types of control functions are calculated from the boundary grid nodesprior to the solution procedure. The calculation of the two types of control functions isdescribed in the following paragraphs.For the convenience of our discussion, the surfaces in a three-dimensional block arenumbered. This is illustrated in Figure 6.1. Face 1, for example, is the boundary surfacewhere is fixed to zero.The stretching functions used here are based on the method of Thomas and Middlecoff(1980). These functions have the following form,= —(xxE + + + y + z.), i = 1,2,3 (6.3)where the superscript C3 indicates that the functions control stretching and the subscriptsindicate differentiation. With these stretching functions, the stretching of the grid atboundaries can be transmitted into the interior.In the MBEGG code, the F’ functions are calculated on all the boundaries whereis not a constant. For instance, P1’ is calculated on faces 3) to 6) but not on face 1) or 2).This is because the P function is defined by the deratives with respect to which canChapter 6. Numerical Grid Generation 148Figure 6.1: Illustration of the face numbers in a 3D block.not be calculated from the boundary nodes on surfaces 1) and 2). This function can becalculated on faces 3) to 6) by using certain difference approximations. In the MBEGGcode, a second order difference scheme is used to discretize all the derivatives involvedin the F functions. After the F functions are calculated, they are then interpolatedinto the interior of the grid using certain interpolation methods. These control functionsare calculated only once at the beginning of the solution process. This type of controlfunction is very useful for most grid generation problems.The control functions used to control grid orthogonality near the boundaries usedin the MBEGG code are based on those proposed by Steger and Sorenson (1979). Thecontrol functions are assumed to have the following form:(6.4)where the index n is the surface number from 1) to 6) (see Figure 6.1). The first termon the right-hand side of equation (6.4) is the stretching function described in the lastparagraph, and the second term is the function to control the orthogonality near the3Chapter 6. Numerical Grid Generation 149boundary. Each R is a function which can be used to control the cell height andorthogonality on face ‘n’. The R are determined as follows. Taking R (i=l,2,3) asexamples, the functions are assumed to have the following form,R(1,23)= (6.5)where at are positive constants. To define appropriate Q2, equations (6.2) are appliedon face 1) to obtain three linear equations for Qi, Q2 and Q. These equations canbe solved to obtain Q (i=l,2,3) if all the derivatives involved in equation (6.2) canbe calculated from grid nodes on face 1). This means that all possible first and secondpartial derivatives of x with respect to must be determined. Derivatives involving onlyor can be found by simply differencing the given boundary grid nodes. However,derivatives involving can not be calculated in an obvious manner from the boundarygrid nodes. These derivatives are found from the constraints imposed upon a line ofvarying which intersects face 1) as follows: (a) it must be normal to the tangent vectoralong the 2 direction, (b) it must be normal to the tangent vector along the direction,(c) the length along that line from the surface to the next node must be controlled. Thesegeometric constraints can be expressed in the following equations:OXi Oxi Ox2 Ox2 Ox3 Ox3 —+ + — 0, 6.6uci ‘c uç ‘ç uciOx1 Ox1 Ox2 Ox2 Ox3 Ox3 —+ .-+ — U,uci uç3 uci uç3 uci 1Jç3=2, (6.8)where S is the desired distance to the first point from face 1). From these equations,the first derivatives with respect to i can be solved. The only derivatives remaining tobe determined are the second partial derivatives with respect to . Those can be foundby differencing the existing solution in the iteration process. After all the derivativesChapter 6. Numerical Grid Generation 150involved in equation (6.2) are calculated, these equations are treated as a 3 x 3 linearalgebraic system for the unknowns Qi, Q2 and Q. The orthogonal control functions Rcan then be found using equation (6.5).With the control functions expressed in equation (6.4), boundary stretching can beconserved in the interior and orthogonal grids at boundaries can be obtained. These twotypes of control functions can be used separately or in combination.6.4 Multigrid SolverThe numerical solution of equation (6.2) can be obtained by iterative methods, such asGauss-Seidel iteration. Although iterative solution procedures have been widely used,they become computationally expensive for finer grids and for higher accuracy. Thus,there is a need to devise faster solution procedures. The multigrid method is an idealsolver for the grid generation system, due to the fully elliptic nature of these equations.In order to obtain an efficient grid generator, a multigrid method is developed to solvethese equations. The Full Approximation Scheme (FAS) of Brandt (1980) in conjunctionwith Full Multigrid Cycling (FMG) is adopted. The basic idea of the FMG scheme is toprovide a good initial solution for the finer grids by prolongating from the solutions oncoarse grids. This avoids the need for an algebraic grid generation procedure for initialsolutions, which is often necessary for single-grid methods in order to avoid singularities.A multigrid method includes four major components: discretization, smoothing (solution solver), restriction, and prolongation, which are summarized in the following paragraphs.Discretization of equations (6.2) is obtained using the second-order, central differencescheme for all derivatives involved. The resulting discrete equations are solved on acuboid domain with a Dirichlet boundary condition. The discrete governing equationsChapter 6. Numerical Grid Generation 151are then written as seven-point formulas as follows,apbp = aEbE + awqw + aNcbN + asqs + arbT + aBcbB + C (6.9)where the variable q expresses the three unknown xk (k=1,2,3) and the subscripts W,E, etc. indicate the six neighbouring positions corresponding to node P.Since the coefficients a23 in equation (6.2) are functions of the coordinate derivatives,the singular case, i.e. J=O, and the anisotropic case where the coefficients a11 ,a22 and a33have different orders of magnitude, e.g., a11 << a22, may occur. The point Gauss-Seideltype relaxation is not efficient for the anisotropic case, such as a11 << a22, since it is onlysmoothed in the 2 direction and not in the i direction. This can be remedied partiallyby simultaneously solving for those unknowns along a line parallel to the 2 direction.To prevent all the possible anisotropic cases, the alternating line-coupled Gauss-Seideliterative procedure is used as the smoother.Singular cases may also happen during iteration (for example, zero initial guess solution will result in a11 = a22 = a33 = 0). For such cases ap 0, and the Gauss-Seideltype relaxation fails as a smoother. To avoid this case, an artificial damping quantityis introduced. That is, a positive quantity 13qp is added to both sides of equation (6.9),where 3 is a positive constant and is chosen depending on the problem to be solved.Therefore, the source term C is augmented by an amount of/3çbp while the coefficient apof the main diagonal term is increased to ap + /3.Restriction is a process by which contents of a fine grid are transferred to a coarsergrid after a few sweeps of relaxation. Two quantities, namely, the approximate solution(x1,x2,x3) and the residual are transferred to the coarse grid by two different types ofoperators. Pure injection is used for the restriction of the solution. In order to preservethe high frequency contents of the fine grid residuals, full weighting is applied in thetransfer operator for residue restriction.Chapter 6. Numerical Grid Generation 152Prolongation is a process by which the corrections and solutions (x1,x2,x3) in thecoarse grid are interpolated onto the fine grid. Its order depends on the order of thedifferencing scheme and for a second order difference scheme, it is sufficient to applylinear interpolation for prolongation to avoid large magnification of the high frequencycomponents. In the program, bilinear interpolation is used.6.5 Multi-Block Grid GenerationAlthough in principle it is possible to establish a correspondence between any physicalregion and a single rectangular type domain by solving the elliptic grid generation equations, the resulting grid is likely to be too skewed and irregular to be usable for certaincomplex 3D configurations, such as that for the film cooling of a turbine blade. For complex configurations it becomes necessary to partition the physical domain into a numberof subregions. Grids are generated in different subregions.In the MBEGG code, a multi-block grid generation method is used to generate gridsfor complex geometries. In this method, the entire complex physical domain is segmentedinto different sub-domains called grid blocks. Grids are generated in each block andpatched together to produce a multi-block grid.The topology to define the block structure in a multi-block method can be complicatedfor complex configurations. In the MBEGG code, the multiblock topology is definedin the following manner. First, a data file is created to provide the basic geometricinformation, such as the number of blocks, the neighbouring blocks, the grid dimensions ofeach block, and block types. Secondly, a Cartesian coordinate system is chosen referringto the given configuration to define the coordinates of each point in the domain.A subroutine is used to define the block structure from the information provided bythe user. Each block is defined by six boundary surfaces. Each surface is defined byChapter 6. Numerical Grid Generation 153four edges. The shape of each edge is specified by the user as either a line or a curvesegment. For line segments, oniy the coordinates of the two end points are required.For a curved segment, a function or certain constraints have to be given. A rule for thegrid distribution in each edge is specified by the user, such as uniform distribution ornonuniform distribution with a given cell-dimension ratio. From this information, thegrid nodes in an edge are calculated. After all the edges have been defined, the gridfor each boundary surface can be generated. There are two types of boundary surfaces;one is a plane surface and the other is a curved surface. For a plane surface, the two-dimensional, elliptic, grid generation equations are solved to generate the smoothest gridon the surface. The grid generation for a curved surface is somewhat more involved. Inthe present study, an indirect technique is used. First, the curved surface is mapped intoa 2D plane area according to certain rules. A grid is then generated for the plane areaby solving the 2D, elliptic, grid generation equations. The grid for the curved surface isobtained by mapping the 2D plane grid back onto the curved surface.With the specified grid coordinates on the six surfaces of each block, the inner grid ineach block can be obtained by solving equation (6.2) while the grid coordinates of the sixsurfaces are used as boundary conditions. For multi-block grids, the blocks share commonboundary surfaces called interfaces. The grids near the interfaces have to be controlledin order to establish good relations between adjacent blocks. The grid lines in adjacentblocks might be made to meet at the interface with complete continuity, with continuousgrid line but with discontinuous slope, or perhaps not to meet at all. Matched grid linesacross the interface can be achieved by simply specifying the same grid coordinates fortwo adjacent blocks at the interface. This type of grid is adopted in all our calculations.For a multi-block grid, smooth grid lines across the interfaces are preferred. This canbe achieved by specifying both the same grid coordinates and same grid line slopes atthe interface. Unfortunately, the elliptic grid generation system only permits one typeChapter 6. Numerical Grid Generation 154of boundary condition to be specified in order to have a solution. In the present study,generating grids with smooth grid lines across the interfaces was not sought. Instead,efforts are made in the CMGFD code to allow for the use of discontinuous grid line slopesacross the interface.6.6 The Code StructureThis code includes a 2D, elliptic, grid generation package and a 3D, elliptic, grid generation package which can be used to generate grids in arbitrary 2D or 3D single-blockdomains. The 2D grid generation code solves a two-dimensional, elliptic, grid generation system and can be used to generate grids on the boundary surfaces of each blockwhich serve as boundary conditions for three-dimensional grid generation. The 3D gridgeneration code solves the three-dimensional, elliptic, grid generation system and can beused to generate grids in hexahedral domains of arbitrary shape. The automatic choiceof control functions is implemented in both packages. Designed grids can be obtained byspecifying the control functions. Both the 2D and 3D code packages are written usinga full multigrid and full approximation scheme. The code package described in chapter5 for the multigrid method is modified for the present use. The MBEGG code is constructed by using the two single-block generation codes and some subroutines for domaindecomposition and linkage.6.7 Application of the CodeThe code has been applied to generate grids in various geometries, such as curved pipesand turbine blades, to support the application of the CMGFD code. The generation ofthose grids, such as the multi-block grids over a turbine blade, have been shown in thechapters where the problems are solved. In this section, only the generation of a grid inChapter 6. Numerical Grid Generation 155a 2D domain is presented to demonstrate the performance of the multigrid method.Figure 6.2 shows the domain to be considered. This domain is related to the configuration of a turbine blade to be studied in the next chapter. Here, the circular edgerepresents the blade surface in the cross-blade section. Such a domain will be used as aboundary surface in a multi-block grid over the turbine blade.Computation is carried out using a 6-level multigrid with the coarsest grid containingonly 2 x 2 cells. Convergence was considered to be achieved when the maximum residuesfor all the three governing equations were less than a prescribed tolerance, which for thesetests was set at 5 x iO. The converged solution is plotted for a finer grid (32 x 32) inFigure 6.3.It was found that the full multigrid method which first started on the coarsest gridhad significant contributions to numerical stability for the grid generation. Because of thestrong nonlinearity of the grid generation equation, an inadequate initial guess will causenumerical instability. Computations using the single-grid method showed that it wasdifficult to converge for fairly dense grids. An algebraic grid generation procedure wasused in the single-grid method to achieve a good initial solution. With heavy relaxationa solution was obtained for a grid of 32 x 32 cells. However, convergence could not beobtained for grid dimensions beyond 32 x 32. The convergence behavior for the multigridand the single-grid calculations are compared in Figure 6.4. It can be seen that themultigrid method produces a much faster convergence rate.6.8 ClosureA multigrid and multi-block grid generation code was developed. An elliptic grid generation method is used, which generates grids by solving three, elliptic, differential equations.Methods to achieve certain grid quantities through the control functions are implemented—.—.CDo..:i.CDp-.o_CDni-iI—’CD0CDCDcz___I-i—.i÷9d I-i i-A.0 CD CD j) I-i 0 xI-’.I-i CD H CD CD CD CD I-i p., CD p., CD I-i I-i i-A.C C r) C 0 C.) 0 0I.-. 0..p.,0) C CI.Chapter 6. Numerical Grid Generation 1570aSymmetry LineFigure 6.4: A two-dimensional domain to be considered.Chapter 6. Numerical Grid Generation 158in the code. A multigrid method was developed to solve these equations. The developedmultigrid algorithm dramatically reduces the computational effort required to solve thesystem of elliptic grid generation equations. Fast convergence has been achieved with azero initial guess. A domain segmentation technique was developed to generate blockstructured grids. The code can be used to generate grids in arbitrary 2D or 3D domains.The combination of the MBEGG and CMGFD codes allows for flow simulations in arbitrary geometries.Chapter 7Prediction of Film Cooling Over A Turbine BladeExtensive experimental studies of leading edge film cooling have been carried out atUBC by Gartshore et al(1993) and Salcudean et al (1993,1994a) based on a large turbineblade model. In this chapter, computations are carried out based on such a blade modeland are compared with the experimental data. This model has a semi-circular leadingedge with four rows of film cooling orifices positioned symmetrically about the stagnationline. Lateral injection is used and the cooling orifices are inclined at 300 to the bladesurface in the spanwise direction.The computational methods and codes described in previous chapters are used toanalyze the flow and related heat transfer in this complex geometry. The computationaldomain follows the blade geometry and includes not only the curved blade surface butalso the inclined circular coolant orifices. The computational domain is segmented intoa number of sub-domains and separate curvilinear grids are employed for different flowregions. Grids are generated using the methods and codes described in the previouschapter. A block-structured, non-orthogonal grid is used to exactly represent the curvedblade surface as well as the circular injection orifices. Computations over the cooledturbine blade model are carried out for overall mass flow ratios of 0.52 and 0.97. Therelative mass flow ratios from each orifice are specified to match experimental values.Density ratios of coolant to free stream were taken to be unity (constant density). Thepredicted results are compared with the experimental results of Salcudean et al (1994).159Chapter 7. Prediction of Film Cooling Over A Turbine Blade 1607.1 Problem DescriptionThe film cooling model under consideration is that of an experimental turbine blade modelat TJBC, as shown in Figure 7.1. This experiment was designed to study the film coolingeffectiveness near the leading edge of turbine blades. The model has a semi-circularleading edge diameter of 127 mm, spanwise top and bottom lengths of 1.2 m (betweenfences), and a chordwise length of 2.3 m. Coolant was injected through four rows ofcircular orifices which are inclined at 300 in the spanwise direction and perpendicularto the blade surface in the streamwise direction on the cylindrical leading edge. Thesefour rows of orifices are located symmetrically at +150 and +440 with respect to thestagnation line. Each injection tube has a diameter d of 12.7 mm and each row hasorifices whose spanwise spacing was 4d.Exploiting the advantage of the flow symmetry about the stagnation line and penodicity in the spanwise direction, only half of the physical region with one period in thespanwise direction was taken as the computational domain. Figure 7.2 shows the domainused in the present study, where d is the diameter of the injection orifice. The computational domain includes a main flow region, two half injection orifices in the first row,and one injection orifice in the second row. The main flow region extends lOd upstreamfrom the leading edge, 40d downstream (in the positive x-direction) from the end of thesemi-cylinder, and 30d in the z-axis vertical direction from the leading edge. Since thecircular injection tubes are inclined at 300 to the spanwise direction, the orifices becomeellipses with a semi-major axis of 2d in the y-direction, as shown in Figure 7.2. Thelength of the injection hole is 4d.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 161Figure 7.1: Illustration of the turbine blade model used in the experimental study ofSalcudean et al (1994).Figure 7.2: Illustration of the computational domain, including a main flow region, twohalf injectionho1e ducts in the first row and one injection hole duct in the second row.55dxChapter 7. Prediction of Film Cooling Over A Turbine Blade 1627.2 Computational GridThe grid was generated using the MBEGG code described earlier. The computationaldomain is segmented into four subdomains, a “hot” flow region over the blade, a coolanthole in the second row and two half coolant holes in the first row. In order to obtain acontinuous grid at the interface, a rectangular type grid is used for the injection holesinstead of cylindrical grids. The generated grid is shown in Figure 7.3. The fine grid usedfor the present study contains 101 x 17 x 37 nodes in the main flow region and 9 x 9 x 19nodes for each injection orifice duct. Figure 7.4 shows the grid on the blade surface andin the injection orifice regions. The grids in the orifice ducts match with the grids in themain flow region at the interface.7.3 Solution ModelSince the comparable film cooling experiments were done in a low, speed incompressibleflow with isothermal conditions and a heat/mass transfer analogy, the computations tocompare with the experiments are made with a steady, incompressible, constant densitycode. The 3D, incompressible, Reynolds-averaged, Navier-Stokes equations together withthe energy equation are solved. Turbulence closure is attained by the use of the standardk — model with the ‘wall function’ treatment described by Launder and Spalding (1974).The shortcomings of the k — model are well known for the prediction of film coolingproblems. These are mainly due to the anisotropic and non-equilibrium turbulence causedby flow curvature, film jet spreading, flow complexity, and unsteadiness near the injectionorifice exits, as well as the heat transfer to the wall in the near-wall, viscosity-affectedsublayer. However, the present calculations are for an adiabatic wall and, therefore, donot involve heat transfer to the wall. It is useful to determine how well the standard k —model can predict the film cooling over an actual geometry with a curved blade surfaceChapter 7. Prediction of Film Cooling Over A Turbine Blade 163Figure 7.3: The grid generated for the present calculation containing 101 x 17 x 37 gridnodes in the main flow region and 9 x 9 x 17 grid nodes in one injection hole.oq CD 0 0 CD I-i CD 0’CD CD CD CD 0 0 0 CD 0CD 0 0 0xI.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 165and the injection orifice ducts.In addition to the Reynolds-averaged Navier-Stokes equations and the k — e equationsdescribed in chapter 3, the following thermal energy equation is solved:•(uT— (- + -) T) = 0, (7.1)where T is temperature, Pr and Pr are the Prandtl number and turbulent Prandtlnumber respectively.7.4 Boundary ConditionsFive types of boundary conditions; inlet, outlet, wall, no-flux, and periodic, are used.The treatment of boundary conditions on each side can be described as follows:• Mainstream:The inlet plane of the mainstream (the west plane of the main flow region) islocated lOd upstream from the leading edge where all the dependent variablesexcept pressure are prescribed to match the experimental values. The mean velocityis taken as U = 10 rn/s following the experiments. The turbulence energy andturbulence dissipation were evaluated using the following formulas:k = 1.5U2(u/U), e = C,k/l (7.2)where u/U is the turbulence intensity and 1 is the turbulence length scale. Theturbulence intensity is assigned a value of u/U = 0.5% and the turbulence dissipation is calculated based on a length scale equal to the height of the wind tunnel1 = 1.6 m. Since the turbulence intensity is very small, these assumptions shouldnot affect the results significantly. The temperature is assigned a non-dimensionalvalue of 1.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 166• Coolant Inlets:The velocities at coolant duct inlets are specified from the mass flow ratios anddischarge rates measured in the experimental study of Salcudean et at (1994a).The k and c values at the duct inlets are evaluated from equation (7.2). Theturbulence intensity u/U is specified as 5% and the lenth scale 1 is taken to be thediameter d of the injection orifice. The temperature is assigned a non-dimensionalvalue of 0.• Adiabatic Wall:The wall at the turbine blade surface (from the leading edge line to the outlet)is assumed to be adiabatic and impermeable with zero tangential velocity. Thestandard ‘wall function’ described by Launder and Spalding (1974) is used. Thefront section before the stagnation line at the bottom of the main flow region istreated as a zero flux and free-slip condition due to the flow symmetry about thestagnation line.• Periodic Condition:The periodic (cyclic) condition is used on the boundaries in the spanwise direction(i.e., the south and north planes). Such a boundary condition assumes that thereare an infinite number of orifices in the spanwise direction.• No-Flux Condition:The top boundary is located at a distance 30d in the z- direction from the symmetryplanes, where the impermeable, free-slip condition is imposed.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 167• Outlet Condition:The outflow boundary is located at a distance 40d downstream from the semi-cylindrical end. The zero-gradient condition in the streamwise direction is imposed for all dependent variables at the outlet boundary. It was sufficiently fardownstream to ensure that the flow in the upstream region was not affected bydownstream conditions.7.5 Computational ProcedureThe CMGFD code described in Chapter 5 was used to carry out computations of flowand film cooling effectiveness for mass flow ratios of M = 0.52 and M = 0.97, whereM is the overall mass flow ratio averaged from all injection orifices which is definedas M = -. Here U2nfty is the free-stream velocity and U is the average coolantinjection velocity. Due to the strong pressure gradient near the leading edge, the massflow from the first row (located at +15° of the semi-circle from the stagnation line) isless than that in the second row (located at +44°). Following the notation of Salcudeanet al (1994a), the mass flow ratio for the first and the second rows are denoted as M15and M44, respectively. The mass flow ratios M15 and M44 were obtained by Salcudean etat (1994a) by measuring the discharge flow rate and are given by the following relations.(1) if = 0.52, then M15 = 0.3 and M44 = 0.75.(2) if M = 0.97, then M15 = 0.86 and M44 = 1.08.These values were used to calculate the turbulence kinetic energy and dissipation rateat the coolant duct inlets. For both mass flow ratios, the coolant duct Reynolds numberRe is equal to 4200 as in the experimental study, where Re pUd is based on the overallChapter 7. Prediction of Film Cooling Over A Turbine Blade 168average coolant velocity U, injection orifice diameter d, and the viscosity and density ofthe air at STP.A two-level multigrid was used for the computations of the mass flow ratio M = 0.52.No special effort was made for the initialization of the solution. On the coarse grid, allthe dependent variables except Ic and 6 were set to zero, and Ic and e were set equal tothe values in the free stream. The converged solution was prolongated to the fine gridwhere the FAS procedure was applied.It was observed that heavy under-relaxation was necessary to obtain a convergedsolution. This is mainly due to the first row of injection holes where the jet exits arelocated in the high pressure region. Convergence was not obtained until a and ct, werereduced to 0.5 and 0.3 respectively, where c is the under-relaxation factor for momentum equations, and a, is the under-relaxation factor for pressure and the Ic — e equations.Furthermore, under-relaxation of the turbulence production rate, the turbulence viscosity, and near wall turbulence quantities was applied to facilitate convergence of thecomputations. For the mass flow ratio M = 0.97 convergence was even more difficultto obtain. To overcome this difficulty, the converged solution for M = 0.52 was usedas the initial solution for this calculation. This approach was found to be very useful.The maximum residuals of all the discrete governing equations dropped by 4 orders ofmagnitude in about 400 iterations for M = 0.52 and 480 iterations for M 0.97.Further iterations were carried out to achieve a residual drop of 8 orders of magnitude.Grid independence studies have shown that, for the number of computational nodes equalto half of that used in the final calculations, the film cooling effectiveness was loweredas much as 20% while decreasing the number of nodes by 30% lowered effectiveness by5%. One can conclude that the results predicted here for the finest grid are approachinggrid independence. Full grid independence can not be achieved with the standard walltreatment in turbulent flows.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 1697.6 Results and DiscussionThe computational results are presented in this section. It should be mentioned that thesolutions are not presented in the complete computational domain in most of the figuresfor reasons of clarity. The flow region near the curved blade surface and the injectionorifice exits contains the most rapid changes of dependent variables, while the region farfrom the curved blade surface tends to be more uniform. Therefore, most of the figuresshow the region near the blade surface and injection orifice exits.Figure 7.5 presents the calculated solutions in the xz-plane and in a plane throughthe centerline of the injection orifice (that is, the north or south plane) in the first row for= 0.52. It should be noted that the orifice length projected on the xz-plane is onlyhalf of the length of the orifice as viewed in the figures due to the inclination angle inthe spanwise direction. Here, Figure 7.5a shows a local velocity field around the curvedblade surface. Figure 7.5b shows the pressure distribution in the xz-plane for the entirecomputational domain. The non-dimensional pressure is expressed as —p---, where PrefPrefis the pressure at the center of the stagnation line. Figures 7.5c & 7.5d show close-upviews of the velocity and pressure fields, respectively, near the orifice exit and curvedblade surface. It can be seen that the orifice exit is located in a high pressure region.The coolant in the front region of the orifice is not blown out directly into the main flowdue to the high pressure near the leading edge which causes some of the hot fluid toflow into the coolant orifice. The velocity in this region is nearly parallel to the bladesurface. Figure 7.5e shows the distribution of non-dimensional temperature , where° = It can be seen that the temperature is higher than the coolant temperaturein a small interior region near the coolant tube exit. Non-dimensional temperature inthis case is equal to non-dimensional concentration, the number 0 being assigned to thepure coolant and the number 1 being assigned to the free stream fluid. Figure 7.5f showsChapter 7. Prediction of Film Cooling Over A Turbine Blade 170Figure 7.5: M = 0.52; Solutions at the streamwise plane (xz-plane) through the centerline of the 150 coolant orifice: (a) local velocity field around the curved blade surface,(b) pressure distribution in the whole xz-plane, (c) local velocity field (near the orificeexit) and (d) local pressure, (e) local temperature, (f) local turbulence kinetic energy.zYx(a)Chapter 7. Prediction of Film Cooling Over A Turbine Blade 171Figure 7.6: M = 0.97; Solutions at the streamwise plane through the centerline of the150 coolant orifice: (a) local velocity field around the curved blade surface, (b) pressuredistribution in the whole xz-plane, (c) local velocity field (near the orifice exit) and (d)local pressure, (e) local temperature, (f) local turbulence kinetic energy.(d)Chapter 7. Prediction of Film Cooling Over A Turbine Blade 172the distribution of the turbulence kinetic energy normalized as It can be seen thatthe turbulence is mainly produced in the regiondownstream of the orifice and is close tothe blade wall.Figure 7.6 presents the results for M = 0.97 in the same manner. Unlike the previousresult for M 0.52, the coolant in the upstream region of the orifice is blown directlyinto the hot flow field, although the velocityprofile at the exit is still strongly nonuniform. From Figure 7.6e it can be seen thatthere is little dilution of the coolant flowin the coolant duct. From both the distribution of non-dimensional temperature andturbulence kinetic energy, it can be seen that the coolant penetrates farther into the hotair for M = 0.97 than for M 0.52. Thisis one of the reasons why a higher massflow ratio may not produce higher film coolingeffectiveness. Figures 7.7 and 7.8 presentsimilar results in the xz-plane through the coolant tube of the second row for M = 0.52and M = 0.97 respectively. Unlike the first row, the second row injection is locatedin a lower pressure region. The lowest pressure occurs just downstream of the injectionorifice.The turbulence kinetic energy is larger at M =0.97 than that at M = 0.52. Itis also larger for the second row than for thefirst row. It should be noted that theturbulence kinetic energy at the entrance of thecoolant orifice which is specified as theinlet boundary condition is not convected to the interior of the orifice. This is because thesource terms are dominant in the k — c equationsand the convection terms are small. Theinlet k and c values have little influence on the interior turbulence and the turbulence ismainly produced by the generation term (contained in the source terms of k—c equations)and ‘wall function’. Therefore, the use of estimated, uniform k — e values at the inletboundaries is likely adequate.Comparison of Figures 7.5-7.8 shows that, as expected, the coolant penetrates furtherChapter 7. Prediction of Film Cooling Over A Turbine Blade 173L(a)Figure 7.7: M = 0.52; Solutions at the streamwise plane through the centerline of the440 coolant orifice: (a) local velocity field around the curved blade surface, (b) pressuredistribution in the whole xz-plane, (c) local velocity field (near the orifice exit) and (d)local pressure, (e) local temperature, (f) local turbulence kinetic energy.-UChapter 7. Prediction of Film Cooling Over A Turbine Blade(a)174Figure 7.8: M = 0.97; Solutions at the streamwise plane through the centerline of the440 coolant orifice: (a) local velocity field around the curved blade surface, (b) pressuredistribution in the whole xz-plane, (c) local velocity field (near the orifice exit) and (d)local pressure, (e) local temperature, (f) local turbulence kinetic energy.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 175into the hot fluid with increasing mass flow ratio. The velocity at the orifice exit isvery complex and depends on the orifice location and mass flow ratio; an assumptionof a uniform velocity profile at the exit plane can be very inaccurate. The pressuredistribution shows that the highest pressure occurs upstream of the first row and thelowest pressure occurs at the downstream edge of the second row. This produces a morenon-uniform distribution of the flow at the coolant tube exit in the first row than in thesecond row. This effect is more noticeable at low mass flow ratios than at haigh massflow ratios.Figure 7.9 shows the development of vortices for M = 0.52 in planes which areperpendicular to the blade surface and which are between the first and the second rows.Here, angle “a” indicates the angular position of the plane on the semi-circular bladesurface. The figure is shown with a longer domain in the y-direction for reasons ofclarity by extending the computational domain using the periodic boundary condition.Figure 7.9a shows the velocity field in the plane through the centerline of the first row oforifices. The velocity is generally directed downward except in the region near the bladesurface. This is not surprising since the plane through the first row forms a sharp angle(15°) with the free stream direction. No vortices are observed in this plane. The highpressure at the exit of the orifices causes the field to flow in both directions and the effectof the lateral injection can not be observed.Figure 7.9b shows the velocity field in the cross plane at a = 22°, which is approximately at x/d = 0.6 downstream from the centerline of the first row of orifices. Figures7.9c and 7.9d show the velocity vectors in the cross planes at a = 28° (x/d = 1.2) anda = 36° (x/d = 1.8) respectively, which clearly show the streamwise vortices. Figure 7.10presents velocity vectors for M = 0.97 in the same manner. Unlike the previous case,the lateral injection is obvious at a = 15° (x/d = 0) and in the positive y — axis directionclose to the wall. No vortices can be observed in this plane, but they can be seen along_________________________________(d)2.01.5N1.00.512-1.0Figure 7.9: M = 0.52; Velocity fields at cross-stream planes which are perpendicular tothe blade surface: (a) a = 15°, (b) a = 22°, (c) a = 28° and (d) a = 36°.176Chapter 7. Prediction of Film Cooling Over A Turbine Blade(a) V=O.5—(b)2.01.5N1.0I I\\ 1////0.53-1 0 1 3 4Y 5Position of the injection holes of the first row(c) V=0.5V=O.52.01.5N1.00.5V=0.5Ill/,,,,,j, jjjjlfI jj’’’’Ij IllIl/I’’’’’’ jljjjll illliii’’’ lilt lIlt ‘t’’’tliii’’’ jill jill lllI’’’’ Ill II/I/jjj tIll 1111111 Ill II/ ///tt t I I I I Ill lI’I’’’I j I I I I/1’’’’’’’ 1111111 ‘‘‘‘‘‘‘‘jI ‘111111111sIll,‘‘S551 I.\f_\\t___,,‘\_%\5l -\_____,I s___- .I - - - ..,,ttlllItl js,i.st itlltllII,t,jjjhlSit jjjjjStttj,,,.._...j tIll55 I Iill_,.ttt lt,.,,,__t,,I jIllli_lilt li.,,’’’’’s _\tsst--- _\_I . -.2.011/1’’’’’’ 1III/// /‘‘‘‘‘‘‘ ,I/,I/I, ‘IiIIIJ// ‘‘III,I///,,,,,,lljIIIJ/// ‘‘‘‘‘‘Ill I////,,,,,II,lI I/////,,,,iilj I////,,,,,iiIjI 1I////,,,,Illj I////,,,,llIIj ///“‘‘‘‘‘ljj 1///‘‘‘IlIjj ////‘‘‘‘‘htlIJj// — — ///‘ — — I I‘i i //: : iIllll III_—1 -‘‘II 1/\\:Figure 7.10: M = 0.97; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a = 15°, (b) a = 22°, (c) a = 28° and (d) a = 36°.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 177(a) V=0.5//// //1.5N1,-’——’!0.5-1 0 1 3 4v 5Position of the injection holes of the first row(c) V=0.5 (d)1.5N1.00.5V=O.52.0,\\\\\\tI,15 .,..\\\\\lt.N /\ \ \ I I II -.3 4Y 5Chapter 7. Prediction of Film Cooling Over A Turbine Blade 178Figure 7.11: M = 0.52; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a = 44°, (b) a = 550, (c) a = 60° and (d) a = 80°.the other three locations. It can be seen that the vortices formed above the orifices aredisplaced spanwise to the right side as the plane moves downstream.Figures 7.11 and 7.12 show the velocity fields in cross planes downstream of thesecond row of coolant orifices. The flow pattern differs significantly from those presentedin the previous figures. Vortices are now formed closer to the orifice centerline and arestronger. It is worthwhile mentioning that the formation of vortices is quite different(a) V=O.5 — (b) V=O.5 —iIlllii 1i1lllII1Ull llIlU. I1iHW2.0 lirirul iTt2.0 [iVII I1i\ I III HIlI1i\ l1Iii\\\ 1l1ii\1.5 I I I I II 1 5 1 1 I I I \ \ 1ii;\\ H IU\\i IT Ill‘\I IN 11liiiIIl\ lrlllIT\ Nii,iiiI lII1lI IiIIIi\\ IIli\1.0 t 1.0 III I Iii ii//,,‘l,-\\ \\0.: 0.5V 1 2 3 6y 7Position of the injection holes of the 2nd row(c) V=O.5 (d) V=O.5 —2.0 Iil\U - Iili\ 2.0 IT lii\\\\ Iljii1l1i\ I\\\ \\\\IIllliiir\\\\\IIIlliiiri\iU\\ l1\ \\\\THIIIIII\\\\\\IIIIIIIII\\ ll\\\ Ii\\ \\\\ITlIlIll\\\\\\\IIIIIIII\\15 I 11 III \ 1.5 \\ \ i I liii \\ \ iI 1 1 I \ 1 I I I I I \\ \ i\ I I 11 1 \\\ \ I ... \\\ \ I I I i i..N II I II I I I ii N \\ I Iii ,,. .‘..\\\ I IiiI! I I I \ \\\ I I I I I \ \\ I I i , , . - I , i1.0\ ti/I 1!I\\\ fI,/I1IIs\\ 1.0\i//I 1i...\\\ ////, \ 1,,,,I//I,! Iii,, /,,,_ ‘I05ujfr’W o5//Za0•p12 356v7 p 1 2 3 5 6v 7from that reported by Vogel (1994) for streamwise injection. In the latter, the largeChapter 7. Prediction of Film Cooling Over A Turbine Blade\\\\\ ill\\\\\\\\\\\i II\\ \ \ II II \\\ \\\ \ \ ii7\\ \ \ li \ \\\\\ \ 1 1 \\\\\I I IIII I\\\\\\\\ illLi \\\\\\\\ iiiN I!II\\\\ TiI,\\\\ 1!!,s’.\Ii,,II/,.. fri I/i/—5’ 5——.5—’I /5-C . •— s. -5ss.-.- 5, 5, -)1 2 3 5 6v 7(d) V=O.5 —Figure 7.12: M = 0.97; Velocity fields at cross-stream planes which are perpendicularto the blade surface: (a) a = 44°, (b) a = 550, (c) a = 60° and (d) a = 80°.179(a)2.0V=O.51.5N1.00.5(b) V=O.5II U\ I HNJ1 UHH hi10,f/ S\\\05Wll/’- :(c) V=O.5 —62.01.5N1.0-5\\\ \ \ \ I \ \ \ ‘ \\\\ \ I I \ \ \ S-_S-S\\\\ II III s ‘‘.\\\\ I 111‘‘‘S-S__..-..\\\\I___‘s\\\\II . _._ssss\\\\l_....—..\\\ I I ..., _\\\ I I.-\\\\I,.., .-\\\II,..Chapter 7. Prediction of Film Cooling Over A Turbine Blade 180Figure 7.13: Contours of predicted FCE for M = 0.52 at the blade surface.vortices tend to lift the coolant from the surface, enhance mixing and thereby decreasingcooling effectiveness. For the present, lateral injection, the vortices are smaller, formfurther downstream, and further from the blade surface than with streamwise injection.The predicted film cooling effectiveness on the blade surface for M = 0.52 is shownin Figure 7.13. It can be seen that the coolant from the first row of orifices flows throughthe surface located between the cooling orifices of the second row, providing good coverage of the blade surface. This broad coverage, with flow from the first row covering thespan between second row of orifices was also obtained by Salcudean et al (1994a). Themeasured film cooling effectiveness is shown in Figure 7.14 for comparison. The spanwiseaverage effectiveness is shown in Figure 7.15 for M = 0.52 together with the experimental curve. Here the x-axis is the curved distance x/d from the stagnation line andthe y-axis is the spanwise averaged film cooling effectiveness. The agreement betweenthe two results is good.Chapter 7. Prediction of Film Cooling Over A Turbine Blade 181B40N2CFigure 7.14: Contours of measured FCE for M = 0.52 at the blade surface.The film cooling effectiveness on the blade surface for M = 0.97 is shown inFigure 7.16. For comparison, the experimental contours of film cooling effectiveness areshown in Figure 7.17 for a mass flow ratio of M, = 0.97. There are some discrepanciesbetween the predicted and measured results. It can be seen that the coolant from thefirst row of orifices is blown toward the orifices of the second row, leaving a significantregion of the blade surface with poor coolant coverage. This observation is in agreementwith experiments (Salcudean et al 1994a). The predicted spanwise averaged film coolingeffectiveness together with the experimental curve is shown in Figure 7.18. The agreement is fair but not as good as for M = 0.52. At high mass flow ratios the coolantpenetrates deeper into the outer flow field and produces a much more complex flow involving anisotropic film jet spreading and non-equilibrium turbulence. Therefore, thetreatment of turbulence using the standard k — c model does not represent the flow asaccurately.X/dChapter 7. Prediction of Film Cooling Over A Turbine Blade 1820.8jJFigure 7.15: M = 0.52; Comparison of the spanwise averaged film cooling effectivenesswith the measured result.3 4 5xj 6 7 8 9 10Figure 7.16: Contours of predicted FCE for M = 0.97 at the blade surface.Chapter 7. Prediction of Film Cooling Over A Turbine Blade-o\NB432CCFigure 7.17: Contours of measured FCE for M = 0.97 at the blade surface.183Figure 7.18: M = 0.97; Comparison of the spanwise averaged ifim cooling effectivenesswith the measured result.0 2X/d2 4 6 8 100.70.6LI0.5Eci) 0.1Chapter 7. Prediction of Film Cooling Over A Turbine Blade 1847.7 ClosureComputations of film cooling are carried out for a UBC experimental blade model using the methods described in this thesis. From the present calculations, the followingobservations are made:• The flow field shows significant differences between the first row and the secondrow because of the pressure variation on the blade surface. Flow at the exit of theorifices is very non-uniform, particularly at the front row.• The interaction between the main flow and coolant shows a different vortex formation in the present spanwise injection from that observed in streamwise injection.The vortices are weaker, form further downstream, and are further from the bladesurface. Therefore, the lifting of the coolant from the surface which occurs forstreamwise injection is not as significant for spanwise injection.• Examination of the film cooling effectiveness shows good coverage for M = 0.52,where the coolant from the first row flows between the orifices of the second row.For M = 0.97 the coolant from the first row of orifices blows toward the orificesof the second row and therefore the coverage is poor. A comparison of film coolingeffectiveness with experimental observations shows good agreement for M = 0.52and fair agreement for M = 0.97.• The present study shows that the computational methods and codes developed inthis thesis have the potential to model the complex cooling process in an actualturbine blade geometry.Further improvements of the computational results may be possible with improvements of the turbulence and near wall modeling, and extension of the computationaldomain to include a part of the plenum.Chapter 8Conclusions And Recommendations8.1 ConclusionsA computational capability has been developed for the prediction of film cooling of turbine blades. This includes the development of two comprehensive numerical codes andthe associated methods. The codes involve several numerical methods, including curvilinear coordinate-based calculations, multigrid acceleration, domain segmentation, andgrid generation. A detailed investigation of these methods is carried out and novel techniques are proposed to improve the methods. Conclusions on theory development andapplication are given in the following.• A computational method using general curvilinear grids is proposed. The methoduses a novel discretization to overcome the difficulties associated with the use of non-smooth grids. The numerical scheme is formulated directly using the coordinate-invariant governing equations and the physical geometric quantities which includethe cell-surface areas, the surface normal vectors, and the cell volumes, insteadof the commonly-used covariant and contravariant vectors. This method avoidsthe coordinate derivatives of transformation which are not well defined for nonsmooth grids and allows for the use of significantly non-smooth grids. A newscheme for handling the non-orthogonal terms in the momentum equations is proposed. This scheme contributes to the main diagonal terms in the resulting coefficient matrix and allows for an implicit treatment of the non-orthogonal quantities,185Chapter 8. Conclusions And Recommendations 186without increasing the number of computational molecules. This treatment of thenon-orthogonal terms enhances computational stability and convergence rate. Thephysical tangential velocity components, resulting from the velocity expansion inthe unit tangent vector basis, are proposed as dependent variables in the momentumequations. A coupled solution procedure is used in place of the pressure-correctionequation associated with grid non-orthogonality. A number of flow problems havebeen studied using the proposed method which include flows in a skewed cavity,flows in a curved pipe and ducts, and flows through a circular cylinder.• Calculations of turbulent flow in curvilinear coordinates have been investigated.The curvilinear coordinate-based method developed for laminar flows in complexgeometries has been generalized to turbulent flows. The k— c two-equation modeltogether with the ‘wall function’ treatment is used as the turbulence closure. Discretization of the k — equations is presented with particular attention to the linearization of the source terms and calculation of the turbulence energy generationrate. Methods are introduced to facilitate the calculation of the energy generationrate and enhance the computational stability. A formulation for the ‘wall function’in general curvilinear grids is presented. The performance of the developed methodis studied through several three-dimensional, turbulent flows.• Multigrid acceleration of three-dimensional, laminar flows in general curvilineargrids is studied. It was found that the discrete governing equations on differentgrid levels can become inconsistent in some curvilinear grids due to complex flowboundaries, thus reducing the efficiency of the multigrid algorithm. A novel technique has been proposed to solve this problem. Computations using grids withboth significant non-orthogonality and strong curvature showed that the developedChapter 8. Conclusions And Recommendations 187multigrid algorithm is very efficient. Unlike the phenomenon reported in the literature, the number of effective grid levels useful for the multigrid method wasnot limited by the geometrical complexities of the computational domain with thepresent method, even for seriously degraded coarse grids.• Several techniques which help the performance of the multigrid method in turbulentflows are introduced. Particularly, it is found that the formulation of the coarse-griddefect equation using the cwall function’ can cause inconsistency of the governingequations between fine and coarse grids. This problem is discussed in detail and anovel practice is proposed to allow the successful implementation of the multigridmethod in turbulent flows. A similar approach was also developed independentlyby Sun (1994).• Calculations using block-structured (multi-block) curvilinear grids are investigated.The method developed for single-domain computations in general curvilinear gridsis generalized to multi-domains by using a domain segmentation technique. Particular attention has been given to the communication between different domains.The performance of the multi-block method is investigated through several computational examples which show that the developed method works well and allowsfor the use of significantly discontinuous grid slopes at interfaces.• One of the contributions of the present thesis is that a block-structured, curvilinearcoordinate-based, finite-volume code, CMGFD, has been developed based on anexisting Cartesian coordinate-based CFD code. The methods described in thisthesis, such as the curvilinear coordinate-based method, have been implementedin the code in several steps, from simple cases to complex ones. At each step,detailed investigations of the method and validations of the code are performedin order to achieve an efficient and accurate code. An efficient multigrid methodChapter 8. Conclusions And Recommendations 188is implemented in the code and promises fast convergence for both laminar andturbulent flows. This code has the capability to deal with arbitrary geometries byusing separate curvilinear grids in different flow regions. There is no restriction forcurvilinear grids to be used and no limitation on the number of segments whichcan be used in the flow region.• A multi-block and multigrid grid generation code, MBEGG, has been developedto support the application of the CMGFD code. The code is based on an ellipticgrid generation method which solves three elliptic, partial differential equations.Methods to achieve certain grid quantities, such as stretching and orthogonality,through the use of control functions have been implemented in the code. Whenthese grid quantities are initially specified by users, the code can choose the appropriate source terms to obtain the desired grid. A multigrid method is developedto solve the three elliptic equations which provides an efficient grid generator. Thecapability of the code is demonstrated through the application to a number of problems. It can be used to generate block-structured curvilinear grids in arbitrary 2Dor 3D domains.• Computations are carried out for film cooling of a turbine blade model at UBC.In these calculations, the flow and the associated heat transfer are solved usinga block-structured, curvilinear grid generated by the MBEGG code which exactlyrepresents the inclined, round film-holes and the curved blade surface. Detailedflow field and heat transfer data are obtained which improves the understanding ofthe complex flow. The predicted vortex structures show how the coolant is liftedand mixed with hot flows. Comparisons with experiments show good agreement atlow mass flow ratios, however, the agreement deteriorates at higher mass flow ratiosdue to deficiencies of the k—turbulence model with ‘wall function’ treatment.Chapter 8. Conclusions And Recommendations 189Further improvement is expected by using more sophisticated turbulence models.8.2 Recommendations for Future WorkSome recommendations for future work are made as follows:• The turbulence modeling is a major source of numerical errors for many engineeringflow calculations. The present study shows that the k — e model with the ‘wallfunction’ treatment is inadequate for predicting film cooling of turbine blades withhigh mass flow ratios. Further studies are needed in this area, possibly using non-isotropic, higher-order closure turbulence models or multiple-time-scale models.The near-wall turbulence treatment to represent the increased mixing occurring athigh mass flow ratios should be further investigated.• At the time of this study the CMGFD code is limited to incompressible flows withconstant density. However, many engineering flows involve variable density, such asthe film cooling problem with cool air or CO2 as coolant. Therefore, it is worthwhileto generalize the CMGFD code to include variable density flows.• The film cooling efficiency of a turbine blade is influenced by many factors, such asthe mass flow ratio, and the shape and location of the injection orifices. Therefore,in order to achieve an optimal design for film cooling, a parametric study shouldbe carried out. This is possible with the numerical codes developed here.• The transfer of data between subdomains is critical for multi-domain calculations.It is often the source of numerical error and slow convergence, especially for nonmatched grids at interfaces. Some methods use one set of overlapping cells whileother methods use two sets of overlapping cells to facilitate the data transfer between neighbouring subdomains. However, there is a lack of detailed informationChapter 8. Conclusions And Recommendations 190about the efficiency of these methods. Further investigations would be useful inorder to achieve an efficient and accurate solver.• The numerical scheme adopted in this thesis is based on the power-law profile ofPatankar (1980) for convection-diffusion which is a low-order scheme. Most of theapplications of higher-order schemes are limited to Cartesian coordinates. Thedevelopment and investigation of higher-order schemes in curvilinear grids wouldbe very worthwhile.• Grid generation is one of the major tasks in numerical simulation for flows incomplex geometries. The developed grid generation code, MBEGG, can be usedfor generating block-structured, curvilinear grids in complex geometries. However,there are several aspects which require further consideration. In particular, thecode has no capability for generating grids with continuous grid slopes across theinterfaces for multi-domains. One possible method is to iteratively adjust the sourceterms near the interface to achieve smooth grids.• Although the CMGFD and MBEGG codes are able to simulate laminar/turbulentflows in arbitrary 2D/3D geometries, the application of these codes is problem/userdependent. It would be very useful to develop these codes into a user-friendlypackage.Bibliography[1] Z. Abdullah, 1993, “Calculations of a hydro-cyclone”, Private Communication.[2] M. Agouzoul, M. Reggio and R Camarero, 1990, “Calculation of turbulent flows ina hydraulic turbine draft tube”, J. Fluids Eng., Vol. 112, pp.257-63.[3] B.S. Baldwin and H. Lomax, 1978, “Thin-layer approximation and algebraic modelfor separated turbulent flows,” AIAA paper, pp.78-257.[4] G.Bergeles, A.D.Gosman and B.E.Launder, 1978, “The turbulent jet in a crossstream at low injection rates: a three dimensional numerical treatment”, NumericalHeat Transfer, Vol.1, pp.217-42[5] M.E. Braaten and W. Shyy, 1986, “Comparison of iterative and direct solutionmethods for viscous flow calculations in body-fitted coordinates”, Tnt. J. Nume.Meth. in Fluids, Vol. 6, pp. 325-349.[6] M. E. Braaten and S. V. Patankar, 1989, “A block-corrected subdomain solutionprocedure for recirculating flow calculation”, Nume. Heat Transfer, Part B, Vol.15, pp.1-20.[7j A. Brandt, 1981, Guide to Multigrid Development, Lecture Notes in Mathematics,960, Springer-Verlag, New York.[8] A.I. Borisenko and I. E. Tarapov, 1968, Vector and Tensor Analysis with Applications, Prentice-Hall, Inc.[9] L. S. Caretto, R. M. Curr, and D. B. Spalding, 1972 “Two numerical methodsfor three-dimensional boundary layers”, Comp. Meth. Appi. Mech. Eng., Vol.1,pp.39-57.[10] L. Davidson and P. fledberg, 1989, “Mathematical derivation of a finite volume formulation for laminar flow in complex geometries”, Tnt. J. Numer. Methods Fluids,Vol.9, pp.531-40[11] I. Demirdzic, A. D. Gosman, R. I. Issa and M. Peric, 1987, “A calculation procedurefor turbulent flow in complex geometries”, Compu. & Fluids, Vol.15, pp.251-73.[12] I. Demirdzic, Z. Lilek and M. Peric, 1992, “Fluid flow and heat transfer test problems for non-orthogonal grids: bench-mark solutions”, Tnt. J. Numer. Methods inFluids, Vol.15, pp.32954.191Bibliography 192[13] A.D.Demuren, W.Rodi and B.Schonung, 1986, “Systematic study of film coolingwith a three dimensional calculation procedure”, Journal of Turbomachinery, Vol.108, 124-130.[14] Q. V. Dihn, R. Glowinsi, and J. Periaux, Solving Elliptic Problems by DomainDecomposition Methods with Applications, in 0. Birkhoff and A. Schoenstadt (eds.),Elliptic Problem Solvers II, Academic Press, New York, 1984.[15] M.M. Eanyet, M.M. Gibson, A.M.K.P. Taylor and M. Yianneskis, 1982, “Laser-doppler measurements of laminar and turbulent flow in a pipe bend’, Tnt. J. Heatand Fluid Flow, Vol.3 pp.213-9[16] M. Faghri, E.M. Sparrow, and S.T. Prata, 1984, “Finite difference solutions ofconvection diffusion problems in irregular domains using nonorthogonal coordinatetransformation’, Numer. Heat Transfer, Vol.7, pp.183-209[17] P. F. Galpin, J. P. Van Doormaal, and G. D. Raithby, 1985, “Solution of theincompressible mass and momentum equations by application of a coupled linesolver”, mt. J. Numer, Methods Fluids, Vol.15, pp.615-2.[18] V.K. Garg, and R.E. Gaugler, 1993, “Heat transfer in film-cooled turbine blades”,ASME paper 93-GT-81.[19] V.K. Garg and RE. Gaugler, 1994, “Prediction of film cooling on gas turbineairfoils”, ASME paper 94-GT-16.[20] I. Gartshore, M. Salcudean, Y.Barnea, K.Zhang, and Aghdasi, F., 1993, “Someeffects of coolant density on film cooling effectiveness”, ASME 93-GT-76.[21] R.J.Goldstein, Advances in Heat Transfer, Vol.7, pp.32l-3’19, 1971.[22] P. He and M. Salcudean, “A numerical method for 3D viscous incompressible flowsusing non-orthogonal grids”, J. of Numer. Method in Fluids, Vol.18, pp.449-6.[23] P. He, M. Salcudean and LS. Gartshore, 1994, “Multigrid calculation of laminarand turbulent flows using general curvilinear grids”, Proc. of CFD94 of Canada,Toronto.[24] P. He, M. Salcudean and LS. Gartshore, “ Computations of film cooling for theleading edge region of a turbine blade model”, Preprint, (1994).[25] P. He, M. Salcudean and I. S. Gartshore, 1993, “A nonlinear multigrid algorithm forgrid generation”, Technical Report, The Department of Mechanical Engineering,UBC.Bibliography 193[26] P. He, CMGFD: A 3D Curvilinear Coordinate-Based Multi-Block Multigrid CFDCode, Technique, CFD group, Dept. of Mech. Eng., The University of British Colubia.[27] M. Hinatsu and J. H. Ferziger, 1991, “Numerical computation of unsteady incompressible flow in complex geometry using a composite multigrid technique”, mt. J.Numer. Methods Fluids, Vol.13, pp.971-97.[28] D. S. Joshi and S. P. Vanka, 1991, “Multigrid calculation procedure for internalflows in complex geometries”, Numer. Heat Transfer, B, Vol.20, pp.61-80.[29] B.A. Jubran, 1989, “Correlation and prediction of film cooling from two rows ofholes”, J. Turbomachinery, Vol.111, pp.502-9.[30] S.C.Kacker, B.R.Par and J.H.Whitelaw, 1969, “The prediction of wall jets flowswith particular reference to film cooling”, Progress in Heat Transfer, Vol.2, pp.163-186.[31] K.C. Karki and S.V. Patankar, 1988, “Calculation procedure for viscous incompressible flows in complex geometries”, Numer. Heat Transfer, Vol.14, pp.295-307.[32] A. Klein, 1981, “REVIEW: Turbulent developing pipe flow”, J. Fluids Eng.,Vol.103, pp.242-9.[33] B. E. Launder and D. B. Spalding, 1974, “The numerical computation of turbulentflows”, Comp. Meth. Appi. Mech. Eng., Vol. 3, pp.269-75.[34] D. Lee and J.J. Chiu, 1992, “Covariant velocity-based calculation procedure withnonstaggered grids for computation of pulsatile flows”, Numer. Heat Transfer, PartB, Vol.21, pp.269-86.[35] M. A. Leschziner and K. P. Dimitriadis, 1989, “Computation of three-dimensionalturbulent flow in non-orthogonal junctions by a branch-coupling method”, Comp.Fluids, Vol.17, pp.371-96.[36] J.H. Leylek and R.D. Zerkle, 1993, “Discrete-jet film cooling: A comparison ofcomputational results with experiments”, ASME 93-GT-207.[37] C.R.Maliska and G.D.Raithby, 1984, “A method for computing three dimensionalflows using non-orthogonal boundary-fitted coordinates”, Tnt. J. Numer. MethodsFluids, Vol.4, pp.51837[38] A.B. Mehendale and J.C. Han, 1992, “Influence of high mainstream turbulence onleading edge film cooling heat transfer” ASME Journal Of Turbomachinery, Vol.114, pp.707-15.Bibliography 194[39] M.C. Melaaen, 1992, “Calculation of Fluid Flows with Staggered and NonstaggeredCurvilinear Nonorthogonal Grids- Theory”, Numer. Heat Transfer, Part A, Vol. 21,pp.1-20.[40] W.J. Mick and R.E. Mayle, 1988, “Stagnation film cooling and heat transfer including its effect within the hole pattern,” ASME Journal of Turbomachinery, Vol.110, pp.66-72.[41] P. Nowak, 1992, “A multigrid and multi-block method”, Technique Report, TheUniversity of British Columbia.[42] C. W. Oosterlee and P. Wesseling, 1992, “A multigrid method for an invariantformulation of the incompressible navier-stokes equations in general co-ordinates”,Comm. Applied Nume. Methods, Vol.8, pp.72134[43] S. Ou and J.C. Han, 1991, “Influence of mainstream turbulence on leading edge filmcooling heat transfer through two rows of inclined film slots”, ASME 91-GT-254.[44] 5. Ou and J.C. Han, 1992, “Influence of mainstream turbulence on leading edgefilm cooling heat transfer through two rows of inclined film slots”, Transactions ofthe ASME, Vol. 114, pp.724-33.[45] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington,D.C., 1980.[46] C.Y. Perng and R.L. Street, 1991, “A coupled multigrid-domain-splitting technique for simulating incompressible flows in geometrically complex domains”, Int.J. Numer. Methods Fluids, Vol.13, pp.269-86.[47] M. Perié, 1990, “Analysis of pressure-velocity coupling on nonorthogonal grids”,Numer. Heat Transfer, Part B, Vol.17, pp.63-82.[48] D. Rayner, 1991, “Multigrid flow solutions in complex two-dimensional geometries”,Tnt. J. Numer. Methods in Fluids, Vol.13, pp.507-18.[49] J.W. Richman and R.S. Azad, 1973 “Developing turbulent flow in smooth pipes”,Api. Sci. Res., Vol.28, pp.419-41[50] W. Rodi, S. Majumdar, and B. Schonung, 1989, “Finite volume methods for twodimensional incompressible flows with complex boundaries”, Comp. Meth. in Appl.Mech. and Eng., Vol.75, pp.369-92.[51] M. Rosenfeld, D. Kwak and M. Vinokur, 1988, “A solution method for the unsteadyincompressible Navier-Stokes equations in generalized coordinate systems”, AIAAPaper AIAA-88-0718.Bibliography 195[52] p. A. Rubini, H. A. Becker and E. W. Grandmaison, A. Pollard, A. Sobiesiakand C. Thurgood, 1992, “Multigrid acceleration of three-dimensional, turbulent,variable-density flows”, Numer. Heat Transfer, Part B, Vol.22, pp.163-77.[53] M. Salcudean, I. Gartshore, K.Zhang, and Y.McLean, 1993, “An experimentalstudy of film cooling effectiveness near the leading edge of a turbine blade”, J. ofTurbomachinery, Vol. 116, pp.71-6.[54] M. Salcudean, I. Gartshore, K.Zhang, and Y.Barnea, 1994a, “Leading edge filmcooling of a turbine blade model through single and double row injection: effectsof coolant density”, ASME 94-GT-2.[55] M. Salcudean, P. He, P. Nowak, I.S. Gartshore and Z. Abdullah, 1994b, CMGFD:A Curvilinear Coordinate-based, Multigrid and Multi-Block CFD Code, The University of British Columbia.[56] P. Sathyamurthy and S.V. Patankar, 1990, “Prediction of film cooling with lateralinjection”, Heat Transfer in Turbulent Flow, pp. 61-70.[57] A. Segal, P. Wesseling, J. Van Kan, C.W. Oosterlee and K. Kassels, 1992, “Invariantdiscretization of the incompressible Navier Stokes equations in boundary fittedcoordinates”, Tnt. J. Numer. Meth. Fluids, Vol.15, pp.411-26.[58] R.K. Shah and A.L. London, Advances in Heat Transfer, Supplement 1, AcademicPress, Inc., 1978.[59] W. Shyy, S.S. Tong and S.M. Correa, 1985, “Numerical recirculating flow calculations using a body-fitted coordinate system”, Numer. Heat Transfer, Vol.8,pp.99-113.[60] W. Shyy and M. Braaten, 1996, “Three-dimensional analysis of the flow in a curvedhydraulic turbine draft tube”, Tnt. J. Numer. Methods Fluids, Vol.6, pp.861-82[61] W. Shyy and C. Sun, 1993a, “Development of a pressure-correction/staggered-gridbased multigrid solver for incompressible recirculating flows”, Comp. Fluids, Vol.22,pp.51-76.[62] W. Shyy and C. Sun, 1993b, “Multigrid computation for turbulent recirculatingflows in complex geometries”, Numer. Heat Transfer, B, Vol.23, pp.79-98.[63] S. Sivaloganathan and 0. J. Shaw, 1988, “A multigrid method for recirculatingflows”, Tnt. J. Nemer. Method in Fluids, Vol.8, pp.417-40Bibliography 196[64] K. M. Smith, W. K. Cope and S. P. Vanka, 1993, “A Multigrid procedure for three-dimensional flows on nono-orthogonal collocated grids”, Tnt. J. Numer. Methods inFluids, Vol.17, pp.887-904.[65] Y. Sun, 1994, “The wall function treatment for multigrid calculations”, PrivateCommunication.[66] J.L. Steger and R.L. Sorenson, 1979, “Automatic mesh-point clustering near aboundary in grid generation with elliptic partial differential equations”. J. Comp.Phys., Vol. 33, pp.405-10.[67] A. M. K. P. Taylor, J. H. Whitelaw and M. Yianneskis, 1982, “Curved ductswith strong secondary motion: velocity measurements of developing laminar andturbulent flow”, J. Fluids Eng., Vol.104, pp.350-9.[68] J.F.Thompson, Z.TJ.A.Warsi and C.W.Mastin, 1982, “Boundary- fitted coordinatesystems for numerical solution of partial differential equations- A Review”, J. ofComp. Phys., Vol.47, pp.1-08.[69] J.F.Thompson, 1984, “Grid generation techniques in computational fluid dynamics”, AIAA Journal, Vol.22, pp. 1505-1521.[70] J.F. Thompson, Z.U.A. Warsi and C.W.Mastin, Numerical Grid Generation, Foundations and Applications, Elsevier, New York, 1985.[71] P.D. Thomas and J.L. Middlecoff, 1979, “Direct control of the grid point distribution in meshes generated by elliptic equations”, AIAA Journal, Vol.18, pp.652-6.[72] S.P. Vanka, 1986a, “Block-implicit multigrid solution of Navier- Stokes equationsin primitive variables”, J. Compu. Phys., Vol.65, pp.138-5.[73] S.P. Vanka, 1986b, “A calculation procedure for three-dimensional steady recirculating flows using multigrid methods”, Comp. Meth. Applied Mech. Eng., Vol.55,pp.321-338.[74] S. P. Vanka, 1987, “Block-implicit computation of viscous internal flows-recentresults”, ATAA paper No. ATAA-87-0058, New York.[75] S.P. Vanka, 1991, “Fast numerical computation of viscous flow in a cube”, Numer.Heat Transfer, Part B, Vol.20, pp.255-61.[76] M. Vinokur, 1989, “An analysis of Finite-difference and Finite-Volume Formulations of conservation laws”, J. Compu. Phys., Vol.81.Bibliography 197[77] D.T. Vogel, 1994, “Navier-stokes calculation of turbine flows with film cooling”,ICAS 94-2.5.3.[78] A.J. White, 1981, “The prediction of the flow and heat transfer in the vicinity of ajet in crossflow”, ASME 80-WA/HT-26.[79] D.F. Young and F.Y. Tsai, 1973, “Flow characteristics in models of arterialstenoses-I. steady flow”, J. Biomech., Vol.6, pp. 395-410.[80] J. W. Yokota, 1990, “Diagonally inverted lower-upper factored implicit multigrid scheme for the three-dimensional navier-stokes equations”, AIAA J., Vol. 29,pp.1642-1649[81] J. Zhou, M. Salcudean and 1.5. Gartshore, 1993, “Prediction of film cooling bydiscrete-hole injection”, ASME 93-GT-75.

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items