Heat Transfer Analysis for Preliminary Design of Gas Turbine Combustion Chamber Liners
by
Nikhil Sharma
A thesis submitted in conformity with the requirements for the degree of Masters of Applied Sciences Graduate Department of Aerospace Engineering University of Toronto
c Copyright 2015 by Nikhil Sharma
Abstract Heat Transfer Analysis for Preliminary Design of Gas Turbine Combustion Chamber Liners Nikhil Sharma Masters of Applied Sciences Graduate Department of Aerospace Engineering University of Toronto 2015 The objective of this thesis was to assess numerical techniques that can be utilized to predict gas turbine combustion chamber liner temperature in preliminary design. There are three main aspects of liner temperature prediction that were explored: (1) hot gas temperature prediction; (2) radiation modelling; and (3) cooling technology modelling. Reactor networks, zonal method along with simple one/two-dimensional models were picked for the these three factors. Preliminary tests for zonal method very low run time with results showing accurate trends. Reactor Networks provided valid trends for combustor outlet temperatures; further validation would be required to assess its capability to predict local temperature. One/two-dimensional models for cooling technologies were tested potential benefits over empirical correlations were discussed. Additional validation would be required for all the sub modules to be integrated into a larger Preliminary Multidisciplinary Design Optimization (PMDO) tool for gas turbine combustion chambers.
ii
Dedication To my family and friends
iii
Acknowledgements I am thankful to my supervisor Professor Sam Sampath for giving me the opportunity to work on this project and for his guidance and support throughout its execution. I am also ¨ thankful to Professor Clinton Groth, Professor Adam Steinberg, Professor Omer G¨ ulder, Professor Gottleib, Sri Sreekanth (PWC) and Haley Ozem (PWC) for their comments and input in the project. I am thankful to my colleagues at UTIAS, and to my friends and family for their constant moral support. Nikhil Sharma University of Toronto Institute for Aerospace Studies September 30, 2015
iv
Contents List of Figures
viii
List of Symbols
xii
1 Introduction
1
1.1
Gas Turbine Combustor Design . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
Combustor Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3
Thesis Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.4
Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2 Background 2.1
6
Cooling Technology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.1
Film Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.2
Double Walled Cooling . . . . . . . . . . . . . . . . . . . . . . . .
10
2.1.3
Effusion and Transpiration . . . . . . . . . . . . . . . . . . . . . .
10
2.1.4
Thermal Barrier Coatings . . . . . . . . . . . . . . . . . . . . . .
12
2.1.5
Liner Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2
Combustor Heat Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.3
Heat Transfer Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.4
Liner Durability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
3 Literature Review and Modelling Strategy
20
4 Radiation Modelling
24
4.1
Zonal Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
4.2
Radiative Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.3
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
4.3.1
Direct Exchange Area Calculation . . . . . . . . . . . . . . . . . .
28
4.3.2
Grid Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
v
4.4
Results and Discussion 4.4.1 Case 1 . . . . . 4.4.2 Case 2 . . . . . 4.4.3 Case 3 . . . . .
. . . .
31 32 34 37
. . . .
40 40 41 42 43
6 Cooling Technology Modelling 6.1 Film Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Double Walled Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Effusion and Transpiration Cooling . . . . . . . . . . . . . . . . . . . . .
47 47 49 53
7 Conclusions 7.1 Thesis Accomplishments . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57 57 58
Bibliography
63
5 Combustion Modelling 5.1 Background . . . . . . 5.2 Cantera Model . . . . 5.3 Implementation . . . . 5.4 Results and Discussion
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
vi
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
List of Figures 1.1
Pressure ratio and outlet temperature trend for large gas turbine combustors
2
1.2
Inlet temperature and aspect ratio trend for gas turbine combustors . . .
3
1.3
Size comparison of an older and a modern combustor . . . . . . . . . . .
4
2.1
Combustor configurations . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2
Axial, reverse flow combustor layout and general airflow in combustor primary zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.3
Traditional cooling technologies . . . . . . . . . . . . . . . . . . . . . . .
9
2.4
Tiled cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.5
Transpiration and effusion cooling . . . . . . . . . . . . . . . . . . . . . .
11
2.6
Combustor wall temperature effect through variation in inlet pressure, inlet temperature, and air mass flow rate . . . . . . . . . . . . . . . . . .
14
2.7
One dimensional heat transfer on combustor liners . . . . . . . . . . . . .
15
4.1
Surface to surface interaction . . . . . . . . . . . . . . . . . . . . . . . .
25
4.2
Approximation used to calculate transmissivity between major grid elements 29
4.3
Effect of integration step size on transmissivity calculation . . . . . . . .
30
4.4
Effect of grid density on calculation accuracy of gas to gas and surface to surface direct exchange area . . . . . . . . . . . . . . . . . . . . . . . . .
31
4.5
Geometry layout for direction vector used in direct exchange area verification 32
4.6
Comparison of discretization schemes that can be used for a cylindrical geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
4.7
Cylindrical geometry to represent a combustor [46]. . . . . . . . . . . . .
33
4.8
Effect of zonal parameters on the net radiative heat flux . . . . . . . . .
34
4.9
Axial temperature, soot profiles and normalized radial profiles used to predict . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
4.10 Effect of zonal parameters on the net radiative heat flux . . . . . . . . .
36
4.11 Effect of radial soot profile variation and constant radial temperature and soot profile on net radiative heat flux on the wall . . . . . . . . . . . . .
37
vii
4.12 Measured vs predicted results for incident radiative heat flux . . . . . . . 4.13 A course discretization grid of primary and secondary zone of a reverse flow combustor for zonal method implementation . . . . . . . . . . . . . 4.14 Combustor cross section schematic with eight sectors . . . . . . . . . . .
38
5.1 5.2 5.3 5.4
Engine A and Engine B profile schematics for comparison . . . . . . . . . Engine A combustor outlet temperature; measured and predicted results Engine A reactor network . . . . . . . . . . . . . . . . . . . . . . . . . . Effect of mechanisms and reactor network setup on combustor outlet temperature for Engine B . . . . . . . . . . . . . . . . . . . . . . . . . . . . Network configurations for Engine B . . . . . . . . . . . . . . . . . . . . Effect of reactor network setup on combustor N Ox emission index (EIN Ox )
43 44 44
Schematic representing Counter Flow Film Cooling (CFFC) . . . . . . . Temperature profiles predicted for panels of various lengths for (a) CFFC and (b) PFFC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Effect of number of nodes used in conduction solver on the overall temperature profile for CFFC . . . . . . . . . . . . . . . . . . . . . . . . . . Predicted results from the effusion model for the wall temperature and the effect of blowing angle on cooling effectiveness . . . . . . . . . . . . . . . Parametric study for effusion cooling. . . . . . . . . . . . . . . . . . . . . Predicted results from transpiration model for wall temperature distribution with and without a ceramic coating. . . . . . . . . . . . . . . . . . .
49
5.5 5.6 6.1 6.2 6.3 6.5 6.6 6.7
viii
39 39
45 46 46
51 52 55 56 56
List of Symbols Variables A Cp D Eb F G gg H h i J j L Lmb k M m ˙ m Nu P Pr Q q Ru Re S St
Area Heat capacity (constant pressure) Diameter Black body radiation emission View factor Incident radiation on a gaseous element Gas to gas direct exchange area Incident radiation on a surface element Convective heat transfer coefficient Enthalpy Radiosity Chilton-Colburn j factor Luminosity factor Mean beam length Thermal conductivity Mass flux ratio Mass flow rate mass Nusselt number Pressure Prandtl number Heat transfer rate Heat flux rate Universal gas constant Reynolds number Path length used in Zonal method Stanton number ix
(m2 ) (J/(m ∗ K)) (m) (W/m2 ) (W ) (m2 ) (W ) (W/(m2 ∗ K)) (J) (W )
(m) (W/(m ∗ K)) (Kg/s) (Kg) (P a) (W ) (W/(m2 )) (J/(mol ∗ K)) (m)
ss sg s T t tr V x Y z
Surface to surface direct exchange area Surface to gas direct exchange area Cooling slot height Temperature Lip thickness Residence time Volume Axial distance on liner wall Mass fraction Normal distance through liner wall
Greek α κ ρ σ τ η ν µ ω˙
Absorptivity Emissivity Absorption Coefficient Density Stefan Boltzmann constant Transmissivity Cooling film effectiveness Velocity Dynamic Viscosity Chemical production rate
Subscripts a an c ca h1 h2 hg i in l out p
Air Combustor annulus Coolant Combustor Casing Effusion hole exit Effusion hole entrance Hot gas i species Inlet condition Liner Outlet condition Potential core
(m2 ) (m2 ) (m) (K) (m) (s) (m3 ) (m) (m)
(Kg/(m3 )) (W/(m2 K 4 ))
(m/s) (Kg/(m ∗ s)) (mol/s)
x
s w wad w1 w2
Slot Wall Adiabatic wall condition Hot side wall Cool side wall
Acronyms AFR CFD CFFC CO DEA DOM DTM FAR LCF MDO MW NOx PFFC PFR PMDO PSR TBC UHC WSR
Air Fuel Ratio Computational Fluid Dynamics Counter Flow Film Cooling Carbon Monoxide Direct Exchange Area Discrete Ordinates Method Discrete Transfer Method Fuel Air Ratio Low Cycle Fatigue Multidisciplinary Design Molar Weight Nitrogen Oxides Parallel Flow Film Cooling Plug Flow Reactor Preliminary Multidisciplinary Design Perfectly Stirred Reactor Thermal Barrier Coating Unburned Hydrocarbons Well Stirred Reactor
xi
xii
Chapter 1 Introduction 1.1
Gas Turbine Combustor Design
Gas turbine engine design involves multiple disciplines, is highly complex and is iterative in nature. Gas turbine design involves three main stages which are: conceptual design, preliminary design, and detailed design stage. Success of an engine can be defined by the conceptual design since decisions made early on in the project have a significant impact on the project. Traditionally, an engineering group utilizes empirical correlations or low fidelity one dimensional models to complete the conceptual design phase, followed by specialized advanced design groups responsible for bringing the concept to production. This results in end-loaded designs where information and effort is added towards the end of the design cycle; changes made at the end of the design cycle are time consuming and expensive [54]. In order to make the design process more efficient modern practices involve more information being added early on in the conceptual stage itself. Using empirical correlations in the conceptual stage can be disadvantageous as it may lead to designs that are conservative in nature and not optimal or competitive. In order to add information in the conceptual design stage, advanced tools are required that can either provide higher accuracy or greater exploration of the design space. However, since the design is in flux in the conceptual design stage the design optimization and decision making have to have a quick turnaround time. Major work has been done in aerospace industry to utilize automated computer frameworks that can incorporate design tools from various disciplines and provide an optimized solution based on the object. Such Multi-Disciplinary Optimization Tools (MDO) are discussed by Panchenko et al. [54]. Multiple engineering disciplines govern the design of each turbo-machinery components in a gas turbine engine. These disciplines include but are not limited to: aerodynamics, dynamics, structures, heat transfer, materials, manufacturing, and acoustics. 1
2
Chapter 1. Introduction 60 Normalized Outlet Temperature
2.00
Pressure Ratio
50 40 30 20 10 0 1940
1950
1960
1970 1980 Year
1990
2000
2010
(a)
1.75
1.50
1.25
1.00 1940
1950
1960
1970 1980 Year
1990
2000
2010
(b)
Figure 1.1: (a) Pressure ratio and (b) outlet temperature trend for large gas turbine combustors [20] The combustion chamber of the engine might not have any moving components but the design and modelling is nonetheless complicated due to the combustion process. Furthermore, there are five major conflicting design requirements that are critical for the combustion chamber; these are: (1) emission control; (2) liner durability; (3) flame stability; (4) altitude relight capability; and (5) exit temperature profile requirement. To prevent combustion liner from reaching critical temperature cooling air is usually added along the liner wall as insulation. This cooling air quenches the combustion products, prevents complete combustion resulting in higher carbon m onoxide (CO) emissions and smoke emissions. Higher cooling flows also reduce the air temperature near the wall, increasing the difference in peak and average exit temperature and making the exit temperature profile worse. Hence, significant effort is required to make the combustor more durable while meeting all other operational requirements [39].
1.2
Combustor Cooling
The problem of durability (and consequently of cooling) is even greater for small gas turbine combustion chambers. Combustor volume scales with total flow level due to engine starting requirements, whereas, combustor cooling requirement would scale with total surface area (to maintain the same temperature the cooling air used per area must remain the same). If the combustor size is scaled, keeping the shape proportionately the same, then the combustor volume would scale to the cube of characteristic length and the surface area to the square of the characteristic length. This leads to a higher surface
3
1.2. Combustor Cooling 6 Combustor Length To Injector Spacing Ratio (L/P)
Normalized Inlet Temperature
1.6 1.5 All Subsonic Non-Recuperated Engines
1.4
Up to 5000 Hp 7000 lb Thrust Engines
1.3 1.2 1.1 1.0
1960
1970 Year
(a)
1980
1990
5 All Engines
4 3 2
Up to 5000 Hp 7000 lb Thrust Engines
1 0
1960
1970
1980
1990
Year
(b)
Figure 1.2: (a) Inlet temperature and (b) aspect ratio trend for gas turbine combustors [9]
area to volume ratio for smaller combustor. The surface area to volume ratio for large combustors may range from 5 to 10 f t−1 meanwhile it may range from 1 to 5 f t−1 for small combustors [39]. This in turn requires that a larger portion of total air has to be used to maintain certain liner temperature. Although for smaller gas turbine engines the combustor operating point temperature is lower as compared to larger engines, this might not be enough to offset higher percentage of cooling that is required [9]. Furthermore, lower cooling flows would require smaller orifices or cooling holes in the liner. However, the tolerances of machining do not scale down with size. Therefore, the tolerance levels would cause greater temperature variation in a smaller combustor than a larger combustor [9]. Higher customer requirements for engine durability and power output also require the combustor engineer to further optimize the cooling flow. Gas turbine manufacturers have been increasing pressure ratio for the past 50 years to increase the power output of the engine. Increase in pressure ratio is accompanied by an increase in the combustor inlet air temperature (usually termed T3 ) which lowers the heat sink capability of the air and increases the combustion temperature. Both leading to higher heat load on the combustor liner. The general industry trends for pressure ratio and combustor outlet temperatures (usually termed T4 ) are shown in Figs. 1.1(a) and 1.1(b) respectively [20]. The data presented in Fig. 1.1 is for large gas turbines but the overall trend is followed by smaller engines as well. Manufacturers have also reduced the combustor size to lower the emission output and to shorten the overall engine. The general industry trends for combustor inlet temperature and combustor size in shown in Figs. 1.2(a) and 1.2(b).
4
Chapter 1. Introduction
Introduced in 1970's
Modern combustor
Figure 1.3: Size comparison of an older and a modern combustor [39] The data presented in Fig. 1.2 is for General Electric combustors but these trends are representative of the entire industry. The drastic change in combustor size can also be seen by an overlay of an older combustor introduced in 1970s and a modern combustor shown in Fig. 1.3. Furthermore, typical combustor life requirements before 1990s would be about 5000 mission hours, whereas, for recent applications the requirement for mission hours per overhaul is in the range of 6000 to 12,000 hours [7, 39].
1.3
Thesis Objective
In light of higher customer requirements and shorter design cycles, use of numerical tools for analysis and MDO becomes essential in the design of gas turbine engine. Therefore, the objective of this project is to explore heat transfer analysis tools that can predict gas turbine liner temperature. A liner temperature prediction tool can be sub-divided into three modules: (1) hot gas radiation prediction; (2) hot gas temperature prediction; and (3) cooling technology modelling. Integration of these modules in larger PMDO tool could be possible, hence, factors effecting this integration would be taken into account in this thesis.
1.4
Thesis Outline
Following this introductory chapter, Chapter 2 of the thesis describes general gas turbine combustor geometrical features and various layouts used in the industry. It also describes common cooling techniques used to protect the liner, factors effecting the liner temper-
1.4. Thesis Outline
5
ature, and a general strategy used to predict the liner temperature. The chapter ends with a brief discussion of durability issues with aviation gas turbine combustors. Chapter 3 focuses on previous work that has been done to model a combustor using high and low fidelity tools and the various levels of success they have achieved. Most of the work discussed focuses on liner temperature prediction or combustor life prediction, however, emission prediction techniques are also discussed as these strategies can be used in a tool to predict liner temperature. Liner temperature prediction is dependent on modelling of (1) combustion (for hot gas temperature prediction); (2) radiation; and (3) cooling technology; Chapter 3 outlines the reasoning behind specific strategies that were selected for each of these sub modules in the context of the literature. Details of these models are discussed in Chapters 5, 4, and 6, respectively. Chapter 7 discuss the conclusions drawn based on the results and future work that can be done to improve the predictions by these proposed sub modules.
Chapter 2 Background Overall combustor layout is largely dependent on the engine type and performance requirements. In terms of the cross section of the combustor three main configuration types have been used in the industry as shown in Fig. 2.1. ‘Can’ combustors consist of cylindrical liners concentrically placed inside a cylindrical casing. They are relatively inexpensive to develop and maintain; however, they generally weigh more than annular counterparts and are seldom used in modern aviation engines [20]. An annular combustor consists of an annular liner located inside an annular casing. These combustors tend to weigh less and are aerodynamically clear resulting in less pressure losses [20]. Tuboannular configuration consists of combustor cans with an annular outer casing, having some advantages from both configurations discussed above. The general industry trend for aircraft engines after 1960s has been to use annular combustors. In terms of heat transfer there are two major differences in cylindrical and annular liners. First, a cylindrical liner allows for simple cylindrical correlations to be used for convective heat transfer calculations. Second, in cylindrical liners radiation in a sector is contained and does not propagate to adjacent sectors; this is not true for annular combustors. Details of combustor heat transfer and a common technique to model is discussed in Sections 2.2 and 2.3 respectively. There are two main axial combustor layouts used in the industry, straight through and reverse flow combustors. These configuration layouts are shown in Fig. 2.2. Advantages of straight through annular combustor are that it provides a compact combustor geometry and low surface area to volume ratio which is beneficial for cooling requirements; disadvantages include longer engine length and higher sensitivity to diffuser operation. On smaller engines higher shaft speed requires close coupling of the compressor and the turbines due to whirling issues [20]. Reverse flow combustors allow for shorter shaft lengths, are compact for engines with centrifugal compressors and reduce overall engine 6
2.1. Cooling Technology
7
Figure 2.1: Combustor configurations [20] length. Hence, this configuration is heavily used in smaller aviation engines. In terms of heat transfer however, reverse flow configuration requires larger cooling flow due to larger surface area when compared to straight-through configuration.
2.1
Cooling Technology
There are various configurations that are used in the industry to cool a gas turbine combustor liner. Here a brief description is provided for three main categories of these configuration, these are (1) film cooling; (2) double walled cooling; and (3) effusion/ transpiration cooling. Film cooling has been traditionally employed on combustors of various sizes due to their manufacturing simplicity. Double walled cooling requiring additional tiles, hence, are more suitable for larger combustors. Effusion cooling requires large number of holes to be drilled on the liner, and has been made possible due to modern manufacturing techniques that lower the cost of such an operation.
2.1.1
Film Cooling
Film cooling refers to cooling concepts that employ a cold film of air to insulate the liner from hot gases. The film gets destroyed downstream of the slot due to turbulence, therefore, usually the slots are provided in 40-80 mm intervals. Some commonly used film cooling devices are wiggle strips, stacked rings, splash-cooling, and machined rings [20].
8
Chapter 2. Background
(a)
(b)
(c)
Figure 2.2: (a) Axial combustor layout, (b) reverse flow combustor layout and (b) general airflow in combustor primary zone [20]
Wiggle strip configuration utilizes the total pressure available in the cooling air (rather than only the static pressure in the case of splash cooling). The combustor liner consists of separate sections with annular height difference between them. A corrugated metal strip is used to connect the various sections, as can be seen in Fig. 2.3(b). Wiggle strip configuration can be very sensitive to small manufacturing variances in the strip’s thickness [20]. Stacked ring shown in Fig. 2.3(c), also uses the total coolant pressure available to form a film. In this configuration, consecutive liner sections which vary in annular height are joined via a metal plate. The metal plate contains holes that are drilled with high precision. The plenum to the aft of the strip is required for the individual jets to coalesce into a sheet. Although being less structurally sound than wiggle strip configuration, the dimensional accuracy of the holes is higher in stacked ring configuration, therefore, providing less variation in cooling air flow rate [20]. Splash cooling only utilizes the static pressure drop across the liner to form the cooling film; the configuration is shown in Fig. 2.3(a). Holes on the combustor liner bleed air and a lip redirects the cooling air in the required direction. Typical lip length is about 3-4 times the slot height which is on the order of 1.5-3 mm [20].
9
2.1. Cooling Technology
(a)
(c)
(b)
(d)
Figure 2.3: Traditional cooling technologies are: (a) splash cooling ring; (b) wiggle strip; (c) stacked ring; and (d) machined ring [20]
The machined ring configuration avoids the braze joint required in stacked ring and therefore, the entire step section is machined from solid metal. This configuration is shown in Fig. 2.3(d). Splash cooling (Louvre cooling) or any other cooling configuration requiring a lip is most likely to crack at the tip of the lip. This is due to high temperature gradients at this location; the film provided by previous slot deteriorates at this point while incoming cooling air is almost at combustor inlet temperature. To avoid this drawback, the holes can be drilled with smaller diameter, resulting in thinner cooling jets which require less distance to coalesce into a sheet. A configuration that utilizes a number of small holes without any lip is the Z ring configuration [20].
10
Chapter 2. Background
HotTile Coolantfilm CoolantDirection
CoolingFins Slot
CoolAir Tile HotGas (a)
(b)
Figure 2.4: Tiled cooling: (a) Impingement with tiles and (b) Counter Flow Film Cooling (CFFC)
2.1.2
Double Walled Cooling
Double walled cooling refers to a class of cooling schemes where cooling air flows through a passage between two surfaces before exiting and forming a film on the surface of the liner. Various configurations are possible depending on the direction of the cooling air in relation to the direction of the film that is formed on the top surface. In Counter Flow Film Cooling (CFFC), the coolant air flows in the opposite direction of the film as shown in Fig. 2.4(b), whereas in Parallel Flow Film Cooling (PFFC) it flows in the same direction. If impingement is used to further enhance the heat transfer to the coolant, then the configuration shown in Fig. 2.4(a) is achieved. These configurations allow for the coolant to act as a heat sink and pick up heat before it forms a layer on top, hence, improving the overall cooling effectiveness and reducing the total amount of cooling air used as compared to traditional cooling louvres. To further enhance the convective heat transfer to the coolant, various finned configurations can be placed in the passage between the two panels as shown in Fig. 2.4(b).
2.1.3
Effusion and Transpiration
With Effusion cooling large number of small discrete holes are drilled directly on the liner as shown in Fig. 2.5(b). Furthermore, these holes are drilled at an angle to provide a twofold advantage in terms of cooling: (1) angled holes provide higher surface area for heat transfer to occur between the liner and the coolant and (2) shallow hole angles result in jets that are less likely to penetrate the hot gas and hence, are better at forming a film downstream from the hole. Advancements in laser drilling have made effusion
11
2.1. Cooling Technology
CoolantDirection
(a)
(b)
Figure 2.5: (a)Transpiration and (b) effusion cooling [20] cooling a viable technology for aircraft combustors. Hole diameters can be upwards of 0.4 mm whereas the lowest angle that is attainable is around 20◦ [20]. Drawbacks of effusion cooling include the increase weight of the combustor, this is mainly from increased thickness that is required to provide buckling strength and higher cost associated with drilling of high number of holes. The cooling effectiveness of this configuration can be improved with holes drilled with diffuser shaped expansion at the exit portion, however, this would further increase the cost of manufacturing. Several geometrical factors effect the cooling effectiveness of effusion cooling; these are briefly discussed here. Cooling effectiveness increases with increased hole size, mainly due to the decrease in coolant velocity [17]. Lower coolant velocity would lead to lower turbulence and hence, decrease the amount of mixing of hot gases with the coolant. Effectiveness can also be increased by using a diffuser shaped cooling hole outlet. This would increase the lateral spreading of the coolant and prevent coolant jet blow-off [17, 53]. The result of increasing span-wise or stream-wise distance between the holes is that the overall effectiveness goes down [17]. Higher span-wise pitch prevents the coolant jets to coalesce and higher stream-wise pitch increases the area that a single film has to cool. Finally, increasing the blowing/ inclination angle results in lowering the overall effectiveness. With lower inclination angle the coolant passage in the liner increases in length resulting in higher convective cooling and it also results in better attachment of coolant film to the hot surface [17]. However, with increase in blowing ratio the effect of inclination angle is decreased [1]. Transpiration cooling is defined as a process whereby a fluid transpires through a porous medium and hence, the temperature of the medium and the coolant is equal at the exit. Large internal area of a porous wall would facilitate large amount of heat removal required for this configuration. Emerging jets would coa-
12
Chapter 2. Background
lesce to form a protective film over the surface similar to effusion cooling configuration. Effusion cooling would approach the behaviour of transpiration cooling as the hole sizes reduce. Transpiration cooling ideally enables a uniform temperature to be maintained on the entire surface where the control of the temperature relies on cooling air flow rate. This method would be very efficient in terms of cooling air requirements, however, extremely small hole sizes or air passages in the porous wall increase the chances of blockage due to external debris and oxidation. Due to these issues with porous materials manufacturers have opted to use multi-laminated sheets that provide quasi-transpiration configuration. High manufacturing cost and lack of mechanical strength make this configuration impractical for application in current combustors.
2.1.4
Thermal Barrier Coatings
Thermal Barrier Coatings (TBC) are ceramic materials that can be applied to the inside of the combustor liner. TBC material is usually low in emissivity and thermal conductivity, thus, reflecting large portion of the incident radiation and forming a thermal insulator between the liner and hot combustion gases. A typical TBC has a metallic base coat and one or two layers of ceramic coatings. Oxidation resistant base coat can also be applied to prevent oxidation/ corrosion damage of the liner. Overall the coating thickness can be around 0.4-0.5 mm which can give temperature reduction on the order of 40-70 K [20].
2.1.5
Liner Materials
New super-alloys are being developed for modern combustor liners, however, metals inherently have limitations for extremely high temperature applications. Monolithic ceramic materials on the other hand, can operate on much higher temperature and provide substantial weight benefits. Some disadvantages of ceramic materials are that they are brittle, prone to attack by hot combustion gases and expensive to manufacture. One method of avoiding catastrophic failure mode due to brittle nature of monolithic ceramics is to introduce particles or whiskers that deflect and arrest cracks.
2.2
Combustor Heat Transfer
Common gas turbine combustor features and the complex airflow that can be expected in the combustion chamber is shown in Fig. 2.2. In modern gas turbine combustors gases can have temperatures peaking over 2100 ◦ C, while nickel and cobalt based alloys that are used for combustor liner material cannot operate above 1373 ◦ C [20]. As a result,
2.2. Combustor Heat Transfer
13
various methods are employed to prevent the transfer of heat to the liner and effective removal of heat from the liner as discussed in Section 2.1. In general, the liner is heated by radiation and convection from hot gases inside it and is cooled by radiation to outer casing and convective heat transfer to annulus air. Thermal gradients produced by hot gases are an order of magnitude larger than the thermal gradients in the liner due to conduction, therefore, a common practice in conducting heat transfer calculation is to ignore liner conduction in the axial direction. Internal radiation forms a major portion of total heat transferred from combustion gases to the liner. Areas where cooling films form an effective barrier between the combustion gases and the liner, radiation is the only means by which heat can be transferred. Radiation from gas turbine fuels consists of luminous and non-luminous component. Nonluminous radiation emanates from heteropolar gases such as carbon dioxide and water vapour, while luminous radiation generally depends on solid particles in the flame. The production of soot particles increases at elevated pressures resulting in high amount of heat transfer through luminous radiation [20]. Outside the liner, majority of the heat is transferred through convection to the annulus air. Due to convective transfer, the temperature of air outside the liner progressively gets higher as it traverses the combustion chamber. Radiation from the liner to the engine casing is less prominent as compared to external convection and it depends on surface temperatures and material properties of the liner and the engine casing. Internal convective heat transfer in gas turbine combustion chambers is extremely difficult to predict since the hot gases are going under rapid physical and chemical changes. Various local factors such as boundary layer formation or destruction can have significant impact on convective heat transfer. Any mathematical model used to predict internal convective heat transfer must take into account the combustor aerodynamics, which is heavily influenced by geometric features such as inlet swirler, diffuser and cooling louvres. The liner material in the primary zone (specifically the dome) is subjected to mixture that has high concentration of fuel as it is in close proximity to the fuel nozzle. In the case where the fuel sprayed by the nozzle impinges on the liner, local hot spot can be created. Therefore, the dome section of the primary zone requires special attention by the designer. Fig. 2.6 shows experimental results obtained for liner temperature with variation in three parameters that change with the operating conditions: (1) pressure; (2) inlet air temperature; and (3) air mass flow rate [21]. These results were part of a paper written by Lefebvre and Herbert [21] which employed a methodology that relied on empirical correlations (this is discussed in detail in Section 2.3). The parameters shown were varied
14 1.6
TL Normalized Liner Temperature ( Tmeas,min )
TL Normalized Liner Temperature ( Tmeas,min )
Chapter 2. Background
1.4
1.2
1.0
0.8 0.80
Measured Predicted 1.12
1.44 1.76 2.08 3 ) Inlet Pressure, ( PP3min
2.40
1.6
1.4
1.2
1.0
0.8 0.8
Measured Predicted 1.0
1.2 1.4 1.6 3 ) Inlet Temperature, ( TT3min
2.0
(b) TL Normalized Liner Temperature ( Tmeas,min )
(a)
1.8
1.6
Measured Predicted
1.4
1.2
1.0
0.8 0.5
1.5 2.5 3 ) Inlet Air Mass Flow Rate, ( WW3min
3.5
(c)
Figure 2.6: Combustor wall temperature effect through variation in: (a) inlet pressure; (b) inlet temperature; and (c) air mass flow rate [21] individually on a test combustor in a lab, whereas, during a flight these parameters would change simultaneously. Fig. 2.6(a) shows the effect of pressure on liner temperature. An increase in pressure increases the emissivity in the flame and in turn increases the radiation heat transfer to the liner, it also suppresses chemical dissociation, increasing the total heat release. Due to both of these reasons liner temperature increases, following an exponential law, whereby at high pressures an increase in pressure would not result in similar wall temperature rise as compared to low pressures [21]. Fig. 2.6(b) shows the effect of inlet air temperature on liner temperature. An increase in inlet temperature increases the flame temperature and total heat transferred to the liner, while it also decreases the effectiveness of cooling air and reducing the amount of heat removed from the liner. Due to these reasons the wall liner temperature increases
2.3. Heat Transfer Modelling
15
Figure 2.7: One dimensional heat transfer on combustor liners [20] with an increase in inlet temperature [21]. Fig. 2.6(c) shows the effect of air mass flow rate on liner temperature. External convection removes more heat from the liner when compared to external radiation, on the contrary, internal convection and internal radiation both transport heat to the liner from hot gases. As a result, increase in mass flow increases effectiveness of convective heat transfer, but overall it increases the amount of heat removed from the liner lowering the temperature [21].
2.3
Heat Transfer Modelling
A number of researchers [21, 36, 16, 35, 37, 45, 38, 3, 48] have used the empirical strategy discussed in this section (or a variation that involves improved correlations). This methodology also becomes a reference to compare the new proposed methodology that would be discussed in later sections. Details of the methodology were taken from a paper by Herbert and Lefebvre [21], which was the oldest source the author could find on the issue. If it is assumed that heat flux in constant in the axial and circumferential direction of the liner, then the balance of the heat flux in the normal direction is given by qr1 + qc1 = qr2 + qc2 = qk1−2 .
(2.1)
This balance is also shown schematically in Fig. 2.7. Here qc1 and qc2 refer to internal and external convective heat flux respectively, while qr1 and qr2 refer to internal and external radiative heat flux respectively. Eq. (2.1) shows this by equating the radiation and convective heat flux inside the liner to the radiation and convective heat flux outside the line and to the conduction through the liner. Convection heat flux is predicted
16
Chapter 2. Background
through Newton’s law of cooling as described by qc1 = hhg (Thg − Tw1 ).
(2.2)
Hot gas temperature (Thg ) can be obtained from combustion rise charts assuming near 100% combustion efficiency in the gas turbine combustion chamber. Lefebvre assumed flow inside the combustor to be highly turbulent and to be similar to flow inside a straight circular duct [21]. Stanton number relation for this assumption is given by St = 0.0283 Re−0.2 .
(2.3)
Using relationship between Stanton number, Reynolds number and Prandtl number the following is obtained: khg qc1 = 0.02 0.2 Dl
m˙hg Al µhg
0.8 (Thg − Tw1 ) ,
(2.4)
where Dl is the liner diameter and Al is the liner cross sectional area. To adjust this methodology for cooled liner the concept of film effectiveness was introduced; film effectiveness is defined as η=
Thg − Twad . Thg − Tc
(2.5)
The combsutor liner temperature (Tw1 ) in Eq. (2.4) is replaced with the adiabatic temperature (Twad ) obtained from the cooling effectiveness correlations. Hence, the problem of convective heat flux prediction is broken into two components: (1) of estimating the heat transfer coefficient; and (2) of estimating adiabatic wall temperature. The heat transfer coefficient can be obtained from Nusselt number correlation given below: x 0.7 Nu = 0.069 Res . s
(2.6)
Eq. (2.6) is valid for a blowing ratio (defined as density, velocity ratio of coolant gas to ρUa hot gas m = ρU ) ranging from 0.5 to 1.3. Conceptually, the film issuing from a slot g is assumed to have three main regions: (1) Potential core; (2) Transition zone; and (3) Fully turbulent region [20]. In case of a cooling louvre, important variables that the film effectiveness depends on are slot thickness, slot depth, distance downstream of the slot, blowing ratio and the slot Reynolds number. For the calculation of external convection (liner to casing) similar analysis is done as
17
2.3. Heat Transfer Modelling
described above without considering cooling geometry. However, hydraulic mean diameter is used instead of cross sectional diameter to take into account the annulus area. Further assumption is made that the temperature of the cooling gas does not change axially and is the same as inlet temperature (T3 ) to give kc qc2 = 0.02 0.2 Dan
m˙ c Aan µc
0.8 (Tw2 − T3 ) ,
(2.7)
where all the gas dependent variables in the equation correspond to the cooling gas in the annulus. The Stefan Boltzmann law is used along with emissivity and absorptivity adjustments to predict radiative heat flux from hot gases to the combustor liner through the following: 4 4 qr1 = σ(hg Thg − αhg Tw1 )
(2.8)
Empirical correlation was employed to relate the hot gas absorptivity and emissivity g 1.5 ( αgg = ( TTw1 ) ), while a correction factor (0.5(1 + w )) is added to correct for non black emission from the liner surface resulting in 1.5 2.5 2.5 qr1 = 0.5σ(1 + w )hg Thg (Thg − Tw1 )
(2.9)
Original equation describing hot gas emissivity for non luminous gases by Reeves was adjusted to take luminosity into account to form hg = 1 − e(−290P
−1.5 L(FAR Lmb )0.5 Thg )
,
(2.10)
here luminosity factor (L) which depends on fuel hydrocarbon content is defined as L = 7.53(C/H − 5.5)0.84
(2.11)
and mean beam length can be obtained from Lmb = 3.4
Volume Surface Area
(2.12)
For the calculation of external radiation (liner to the casing) view factor analysis for two surface enclosure was employed [47] to obtain 4 σ(Tw2 − T34 ) qr2 Aw = 1−w . c + Aw F1. wc + 1− w Aw c Ac
(2.13)
18
Chapter 2. Background
Using a view factor of unity for long annular space, estimated emissivity values for steel, and diameter ratio the following equation was obtained 4 − T34 ) qr2 = 0.6σ(Tw2
(2.14)
Since this technique does not require large input data it can be useful in preliminary design. Good results for liner temperature prediction using this technique were shown in Fig. 2.6 and discussed in the previous section; although the operating conditions for these results were not as extreme as in modern combustors.
2.4
Liner Durability
The thrust and power settings of an aircraft engine varies during aircraft flight depending on the aircraft mission, which heavily influences the internal loads that the engine is put under. The gas turbine combustor is put under high dynamic loads due to combustion processes. Some desirable properties for the materials used are high strength, creep resistance, oxidation resistance, low thermal expansion, high thermal conductivity and ease of welding. Generally, super-alloy metals are used with high contents of nickel and chromium such as Hastelloy X and Inconel 617 [42]. The major failure modes of the combustor are creep, Low Cycle Fatigue (LCF) and oxidation which either directly result in surface cracking or exacerbate the issue by weakening the structure. The interaction between each of these failure modes can be complex which makes life prediction or failure modelling very difficult [42]. Creep can be defined as permanent deformation of metal at elevated temperature for extended periods of time. This can happen during aircraft cruise phase and with stress levels well below the yield strength of the material. Creep can become an issue for hot gas turbine components when homologous temperature (ratio of actual to melting temperature) exceeds 0.5. During cruise phase of an aircraft hot components can be well above this temperature. Two methods commonly used in the industry to predict creep are Larson-Miller parameter and Theta-projection method [42]. Low cycle fatigue is associated with a low number of cycles to failure (between 10,000 and 100,000 cycles). Each flight of an aircraft can be considered a single LCF cycle; for a combustor, start up and shut-down cause temperature changes resulting in thermal stresses. LCF can be correlated with the inelastic strain range, however, inelastic strain range can be hard to determine, therefore, LCF is also correlated with total strain range [42].
2.4. Liner Durability
19
Creep and LCF interaction is specially important where the damage by both the mechanisms is almost equal. Synergistic damage process caused by cyclic thermal and mechanical loading is called Thermo-mechanical fatigue. Life time under such loadings can be very different from what is obtained with isothermal LCF tests conducted at maximum temperature where creep damage would be the maximum. Simple analytical methods apply a linear accumulation law to life due to LCF damage in absence of creep and life due to creep damage in absence of LCF to predict overall life [42]. Due to the oxygen content and high temperature of combustion gases, hot sections of the engine are susceptible to oxidation and corrosion. Contaminants in the fuel and air such as Sulphur and Vanadium can cause corrosion attack on high surface materials. Alloys with higher content of chromium and aluminium tend to be more resistant to oxidation and can be used as surface coatings to further protect the combustor liner fromn oxidation[42, 30].
Chapter 3 Literature Review and Modelling Strategy Hot gas temperature prediction can be considered a preliminary step for combsutor liner temperature prediction. The designer may have this information from previous designs or can estimate this using combustion temperature rise charts (assuming certain combustion efficiency as shown in Section 2.3. This technique was also used by Odgers and Kretschmer [16], however, instead of using a single efficiency factor throughout the combustor they divided the combustor in three main zones (primary, secondary and dilution zones) and varied the efficiency and temperature in these zones. Odgers et al. [37] went over the limitation of the liner temperature calculation approach outlined in Section 2.3 for smaller combustors and Gosselin et al. [38] improved the correlations for better results with small combustors. However, it is difficult to extrapolate the results form these correlations and they might not be applicable to new combustor designs. Another method to predict temperature, along with flow properties is to divide the combustor in smaller sub flows, and model these with empirical pressure-drop/ flow rate equations. The PhD thesis by Stuttaford [45] used such a method for modelling gas turbine combustors for preliminary design analysis. The network of sub flows also employed momentum and recirculation effects through empirical factors. Terms relating heat flow between nodes were modelled as sources in the flow equations and the conjugate heat transfer method utilizing equations discussed in Section 2.3 was used for liner temperature prediction. Results of the study showed agreement with one-dimensional industry code used for combustor design and also with experimental data. A drawback of using such a technique is the extensive use of empirical correlations for flow modelling. As discussed earlier accurate data might not be available for newer combustors. A methodology that can relate hot gas residence time (which is related to combustor 20
21 volume and mass flow rate of fuel and air), temperature and pressure to hot gas temperature is reactor network analysis. Reactor networks make the assumption that the turbulence is infinitely high and therefore, the temperature and species concentrations are limited by the chemical reaction rates. The equations for mass, energy and species concentration reduce to non linear algebraic equations which can be solved numerically [49]. The flow is usually represented through interconnectedness of multiple reactors, where each reactor solves the equations aforementioned. Hammond and Mellor [14], and Rizk et al. [32] used reactor networks for emission prediction, while Swithenbank et al. [13] and Sturgess and Shouse [44] used reactor networks to predict flame stability (blow off limits). More recently, Marchand [27] and Lanewala [18] used reactor networks to predict CO and NOx emissions respectively; these projects used CFD flow information that was available to direct the layout of the networks. Although liner temperature was not verified in these models, they demonstrate the ability of reactor networks to predict species concentration (CO and NOx ) with accuracy. Bradshaw [3], in his study focused on impact of manufacturing variability on combustor liner durability. Simplified models were used to link combustor life, liner temperature variability and effects of manufacturing variability; probabilistic analysis was then applied on these models to asses combustor life. A simple reactor network with one reactor to model the entire combustor was used to obtain bulk gas temperature based on various air fuel ratios and inlet temperatures. Heat transfer analysis done in this study was similar to what has been discussed in Section 2.3. Although none of the literature discussed above compared liner temperature with predicted results from reactor networks, the technique has significant advantages for preliminary design temperature prediction. Reactor networks do not require large number of inputs, are scalable (as the networks can be made as complex as required) and can potentially model emissions and blow out. Emission prediction and blow out are not the focus of the current project, however, these capabilities might be used in a larger PMDO tool in the future. Chapter 5 goes over the details of the implementation of reactor networks for the current project. The importance of radiation modelling for liner temperature prediction is highlighted by Lefebvre [20]. Viskanta and Menguc show that radiation flux can be 30-50% of the total heat flux on the liner [26]. Carvalho and Coelho [6] conducted high fidelity simulation of a representative gas turbine combustor sector; radiative exchange was calculated through Discrete Transfer Method (DTM). Percentage of radiative flux to total heat flux for 5, 15, and 25 bar operating points was calculated to be 48%, 43%, and 44% respectively. Although this work was based on numerical simulations, it highlighted the
22
Chapter 3. Literature Review and Modelling Strategy
importance of modelling radiation for combustor heat transfer. The results also revealed higher radiative flux in the primary zone of the combustor. Lefebvre [19], in his review on radiation in gas turbine combustors, discusses the effect of spray atomization in axial radiative flux variation in a combustor. Radiation in gas turbine combustors is found to be dominated by luminous radiation, which is radiation due to soot particles. Atomizers control the spray characteristics and soot forming regions. Since, majority of the fuel is concentrated near the atomizer (in the primary zone), radiation is higher in the primary zone and drops towards the end of the combustor. Axial variation in radiation heat flux can have an impact on the final cooling configuration picked by the combustor designer. Empirical correlations that produce axially average radiation flux might not provide enough detail required to select cooling schemes at various axial positions of the combustor. There are a number of papers in literature that discuss implementation of high fidelity tools for accurate radiation modelling. Details of these implementation, are not discussed, since high fidelity modelling would be out of the scope of this thesis. However, Stuttaford [46] in his paper, assesses the viability of using Discrete Ordinates Method (DOM) for preliminary combustor design. Stuttaford shows validation of the method through multiple test cases where comparison was done between measured thermal paint data for a combustor and predicted results from DOM. However, the computational time was not reported for these cases. On the other hand, the zonal method is a widely used method to solve radiative transfer for practical engineering problems and the computational time required for this method is usually smaller than the time required by the alternative methods [50]. Disadvantages of the zonal method include: (1) high computational time if it is coupled with high fidelity finite difference solvers that contain fine grids; (2) high computational time if the geometry is complex; and (3) inability to treat non-gray, temperature dependent radiative properties [50]. For preliminary design, an assumption was made that impact of minute geometrical detail is not required for radiation prediction. Also, that hot gases within the combustor are gray and the combustor liner is a gray emitter as well (this is true if the liner is covered with soot deposit). Under these assumptions, the zonal method is a viable tool for preliminary design and the details are discussed in Chapter 4. Although, development of a complete PMDO optimization tool is outside the scope of this thesis, various strategies to take heat transfer into account in such tools were explored. Tietz and Behnredt [48] developed a software tool for preliminary design of a gas turbine combustor. The objective of the tool was to optimize the amount of cooling air used for a given metal temperature, minimize NOx emissions and maintain
23 combustor volume for stable combustion. Reactor modelling was used to predict hot gas temperature in the combustor using a simple network while a more dense reactor network was proposed to predict UHC and CO emissions in the future. In order to calculate heat transfer through the combustor liner, a methodology similar to the one described in Section 2.3 was used. The tool focused on modelling effusion cooling with further option of expansion to other configurations in the future. Major advantage of the tool is the detailed modelling of effusion cooling as it takes into account conduction through the liner surface (normal to hole face) and radially outward from the hole, thus providing a two-dimensional analysis. Since heat transfer modelling of the cooling configuration and emission modelling both depend on local fuel to air ratio and air flow split, the optimization was done on specific modules (contour design, cooling air and emissions) and iterated over one by one until the solution converged. Pegemanyfar and Pfitzner [40] developed a software tool for preliminary design of gas turbine combustors as well, however, the emphasis was not on optimization but on automated design process. Operating requirements such as flame stability (relight and weak extinction limit) determined the volume of the combustion chamber and the local fuel to air ratio. After the initial layout of the combustor was picked, the designer would be able to pick a cooling configuration and iterate the design until local air to fuel targets were met. Results of preliminary design were used to develop a parametric CAD geometry. The CAD model was converted to a CFD grid using automatic grid generation software, and a final CFD simulation was run. From the two tools discussed above, it can be inferred that any optimization tool would require an iteration between local temperature prediction, fuel and air mass flow rate adjustment and cooling performance prediction. It is possible to achieve this by using the Reactor network and zonal method. The details of modelling various cooling strategies is discussed in Chapter 6
Chapter 4 Radiation Modelling 4.1
Zonal Method
Zonal method was initially formulated by Hottel and Cohen in 1958 to be used for gas furnaces but it lends itself well to preliminary design analysis for gas turbine combustion chambers. In this method a radiative enclosure is divided in several isothermal surface and volume zones. Radiative exchange between any two zones is dependent on radiative exchange areas. An energy balance on each zone leads to a system of equations equalling the number of zones and can be solved for unknown temperatures or heat fluxes. The equations presented here are taken by a textbook on radiation by Modest [31]. The formulation would only be briefly discussed and the reader can refer to the book for further details. Modest [31] discusses equations for gray, absorbing, emitting, nonscattering medium with constant absorption coefficient, enclosed within a gray, diffusely emitting and reflecting surface. The most significant step in zonal method is the formulation of Direct Exchange Areas (DEA). DEA can be described as the fraction of radiative heat flux arriving at receiving element over the total flux leaving the emitting element; given by Qi→j = si sj Ji .
(4.1)
Here the transfer is taking place between two surface components as shown in Fig. 4.1. Note that the total flux leaving the element is its radiosity (sum of emitted black body radiation and reflection of all incident radiation), defined by Jj = j Ebj + ρj Hj . 24
(4.2)
25
4.1. Zonal Method dAj cos θj nj
Aj , Tj dAj θj
S
ni
ni
θi
dωj−i
dAi dAi
Ai , Ti (a)
(b)
Figure 4.1: Surface to surface interaction
Considering the geometrical location of the two surface, radiation transfer from surface to surface can be described by dQi→j
Ji dAj cos θj = (dAi cos θi ) τij . π S2
(4.3)
The first term on the right side of the equation is radiation that leaves the first element (per solid angle), the second term is the projected area of the first element in the direction of the second element, the third term is the solid angle subtended by the two elements and the last term is transmissivity (i.e. radiation transferred after attenuation due to absorption on path length S). The Beer-Lambert law defines transmissivity for the case with no scattering as τ = e−
R
κ(S) dS
.
(4.4)
Scattering behaviour of particulate matter is governed by Mie theory; in cases where the particles are small relative to the irradiated wavelength such as soot particles in gas turbine combustors, scattering can be neglected [19]. For the special case where absorption coefficient is constant over the path length transmissivity reduces to τ = e−κS
(4.5)
Comparing the definition of a DEA and the formulation in Eq. (4.3), we get the
26
Chapter 4. Radiation Modelling
following expression for surface to surface DEA: Z Z si sj =
τij Ai
Aj
cos θi cos θj dAj dAi . πS 2
(4.6)
For interaction of gas to gas elements, DEA definition, radiative heat transfer and DEA formulation is given respectively by the following: Qi→j = gi gj Ebi
(4.7)
dAj dQi→j = 4κi Ebi dVi τij κj dSj 4πS 2 Z Z κi κj gi gj = τij 2 dVj dVi . πS Vi Vj
(4.8) (4.9)
For interaction of gas to surface elements, DEA definition, radiative heat transfer and DEA formulation is given respectively by the following: Qi→j = gi sj Ebi
(4.10)
dAj cos θj dQi→j = 4κi Ebi dVi τij 4πS 2 Z Z cos θj τij gi sj = κi dAj dVi . πS 2 Vi Aj
(4.11) (4.12)
Note that with these definitions all DEA’s follow the rule of reciprocity given by si sj = sj si ,
gi sj = sj gi ,
gi gj = gj gi .
(4.13)
The total radiation reaching a surface is the sum of radiation from all other surfaces (given by the surface to surface DEA and multiplied by emitting surface’s radiosity) and radiation from all the volumes (given by gas to surface DEA and multiplied by the black body radiation for the gas volume). This formulation is given by Ai Hsi =
N X
si sj Jj +
j=1
K X
si gk Ebgk ,
i = 1, 2, . . . , N.
(4.14)
k=1
Similarly, summing the radiation on a gaseous element results in κVi Gi =
N X j=1
gi sj Jj +
K X k=1
gi gk Ebgk ,
i = 1, 2, . . . , K.
(4.15)
27
4.2. Radiative Properties
Radiosity term in Eqs. (4.14) and (4.15) can be eliminated with Eq. (4.2) and the final form of equations is obtained where the unknowns are surface incident radiation and gas incident radiation. The resultant is a linear system of equations with number of equations (N + K) equalling the number of unknowns (N surface radiation and K gaseous radiation).
4.2
Radiative Properties
For gas turbine combustion application radiation can be divided in two major modes, luminous and non luminous. Non luminous radiation is from water vapour and carbon dioxide gases present in combustion mixture. Since these are molecular gases their radiation is in specific spectral bands. Luminous radiation is emitted by soot and at higher pressures can contribute much more than non luminous radiation [20]. Accurate prediction of radiative heat flux requires accurate knowledge of radiative properties. This can be a challenge since radiative properties depend on species concentration, temperature, pressure, wavelength and path length [31]. For zonal method the only radiative property required for interacting gas is the absorption coefficient. Extensive review of methods to predict absorption coefficient is out of the scope of this thesis, however, general techniques applicable for preliminary design are discussed below. One strategy to predict average absorption coefficient is to use a relevant correlation based on local temperature. Correlation shown below was produced by Gibb and Joyner based on data obtained from combustor testing: Thg
κhg = 0.32 + 0.28e− 1135 .
(4.16)
A drawback of this correlation can be that it does not scale directly with pressure. Although the temperature would rise with higher pressure resulting in lower absorption coefficient (for temperatures ranging from 700 to 2100 K the absorption coefficient from this correlation varies from 0.47 to 0.36 m−1 ), the full impact of pressure might not be taken into account. Another strategy that is employed is computational radiation codes and can be employed with zonal method is to obtain average gray absorption coefficient from total gas emissivity. Equation shown below describes this relationship: κhg = −
1 ln(1 − hg ). Lmb
(4.17)
28
Chapter 4. Radiation Modelling
The Kayakol et al. [33], Viskanta and Menguc [50], and Lallemant et al. [34] all mention this with Kayakol using it with Discrete Ordinates Method (DOM) to obtain absorption coefficient for each cell. In absence of detailed information an average mean beam length for an arbitrary shape (in this case the combustor’s shape)can be obtained by the Eq. (2.12). Using this formulation with non homogenous gray medium is not mathematically sound, however, due to prohibitive computational effort required for detailed formulation this formulation is used in several high fidelity codes [34]. If water and carbon dioxide concentrations are assumed constant throughout the combustor, total emissivity correlations can be used to obtain an average absorption coefficient. To predict total emissivity for a luminous flames , the following equation can be used [19]: tot = gas + soot − gas soot .
4.3 4.3.1
(4.18)
Implementation Direct Exchange Area Calculation
As mentioned earlier calculating DEA’s is the most significant step in zonal method; the rest of the calculations involve simple linear algebra and would not be discussed in detail here. A strategy that was adopted to estimate DEA was to conduct direct numerical integration; this technique is described by M´echi’s [41] in detail. A given geometry is divided into surface and gaseous elements. The elements are considered isothermal zones where properties such as emissivity, temperature and absorption coefficient are assumed to be constant. To calculate DEA with direct integration each element is further subdivided into a finer grid. If a surface element is subdivided in a grid of ‘k’ by ‘k’ sub elements then Eq. (4.6) is transformed into si sj =
k X k X m=1 n=1
τij
cos θi cos θj dAj dAi . πS 2
(4.19)
With this definition of DEA, finer resolution of discretization should result in higher accuracy of total DEA calculated. In the formulation discussed by M´echi [41], transmissivity is calculated for non homogeneous medium between each sub element which can be computationally very expensive, therefore, for this project an average absorption coefficient was calculated between given two elements. It was assumed that the average coefficient did not vary with the location on each
29
4.3. Implementation
B4
B4
B3
B3 B2
B2
B1
B1 A3 A1
(a)
A4
A3
A2
A1
A4 A2
(b)
Figure 4.2: Approximation used to calculate transmissivity between major grid elements element, and it was calculated on path length joining the center of the two elements as shown in Fig. 4.2(b). Fig. 4.2(a) shows the original case where absorption coefficient could vary from one sub element to the other. To calculate the total attenuation due to non homogeneous absorption coefficient the integral in Eq. (4.4) is calculated through direct numerical summation as well. The path length ray S is traced from the emitting to the receiving element and further divided into smaller elements. At each element the distance from the current node’s center point and each adjacent node’s center point is calculated. The smallest distance indicates the node that current step lies in. Absorption coefficient for that node is added to a running total which is averaged based on ray divisions in the end. This method of calculating inhomogeneous absorption coefficient requires that information regarding adjacent element nodes be available at each step of the calculation. This method allows for accurate the integration of transmissivity term given that the ray is subdivided into fine enough elements. Effect of a coarse and fine path length step size is shown in Fig. 4.3. A finer step size result in more information being added in the formulation (violet elements being added ) Ideally, the ray subdivisions should be an order of magnitude smaller than the smallest edges of the volume or surface zones, enabling representative traversal of these zones and better weighting for the average coefficient. A major disadvantage of using direct numerical integration for absorption coefficient estimation is that as the size of isothermal zones becomes smaller and smaller step size is required, computational time can become high. To verify the implementation for DEA described above, results were compared to the ones reported by Yuen et al. [52] for normalized exchange factors. Yuen et al. [52]
30
Chapter 4. Radiation Modelling
B
B
A
A
Figure 4.3: Effect of integration step size on transmissivity calculation. Reduced step size adds information from new cells (shown in purple). provided normalized exchange factors for various combination of direction vectors and optical depth (κD). Comparison of calculated and reported results for gas to gas and surface to surface interactions are shown in Figs. 4.4(a) and 4.4(b) respectively. Fig. 4.5 shows the position used for direction vector corresponding to (x=5, y=5, z=5). Fig. 4.5 showed expected effect of number of sub elements on accuracy of DEA calculation; it also served as a verification test for the current implementation for DEA calculation (with the exception of transmissivity calculation).
4.3.2
Grid Generation
Implementation of a rectangular prism grid is the most straightforward. A small rectangular prism is created to represent an isothermal zone; this is then copied and shifted in three major axis for the corresponding required lengths to form the final geometry. For a cylindrical geometry, small volumetric divisions are created for each radial position. This arc is then copied and rotated circumferentially and then copied and shifted axially to form a cylinder. If the three major geometrical lengths (axial, radial, and circumferential for a cylinder) are evenly divided for refinement, the method described above creates a cylindrical geometry that is extremely fine in the center, however, may still be coarse on the surface of the cylinder. In order to refine a cylinder more evenly a specific scheme was adopted whereby the cross sectional area, and by extension the volume (since axially it is evenly divided) of each sub-volume is kept constant. This requires that for each radial position the number of circumferential divisions varies and it increases as the radial position moves from the inside of the cylinder to the outside. For preliminary studies an initial division
31
4.4. Results and Discussion 6
12
Dir (1,1,5), κD =0.1
Dir (1,1,5), κD =0.1 Dir (3,3,3), κD =0.1
10
Dir (3,3,3), κD =0.1
5
Dir (3,3,5), κD =0.1
Dir (3,3,5), κD =0.1 Dir (1,1,5), κD =2.0
Dir (1,1,5), κD =2.0
6
Error (%)
Error (%)
Dir (3,3,5), κD =2.0
2
2
1 101 102 Grid Density (#)
103
(a)
Dir (3,3,5), κD =2.0
3
4
0 0 10
Dir (3,3,3), κD =2.0
4
Dir (3,3,3), κD =2.0
8
0 0 10
101 Grid Density (#)
102
(b)
Figure 4.4: Effect of grid density on calculation accuracy of (a) gas to gas and (b) surface to surface direct exchange area of three was picked which meant the circumferential divisions varied from 3,5,7,9 for radial positions 1,2,3,4 respectively. Fig. 4.6 shows the comparison of a regular circular grid and the new scheme, both with equal refinement. A note should be made here about the axisymmetric nature of the problem. Axisymmetric assumption allows for significant time reduction in the computation of all radiative matrices. For example, in the case of an axisymmetric cylinder with 3 sub-grid model shown in Fig. 4.6, radiative interaction between all elements within sub grid A are equal to radiative interaction between all elements within B and C. Similarly only radiative interaction between AB is equal to BC and CA. It should be noted that all execution times reported here take this implementation into account. The direct exchange area is assumed to be the sum of exchange areas from all the fine subdivision of the emitting element to all the fine subdivision of the receiving element. The principal of reciprocity is held for DEA which means that the DEA1→2 is the same as DEA2→1 . Hence, to minimize computation time the DEA is calculated only once for each pair of zones. Direct integration technique requires that if an accurate DEA is required then the second sub grid level should be finer.
4.4
Results and Discussion
Implementation of zonal method as described in the previous section was tested against data from three sources. In all cases the radiative heat flux over a cylindrical geometry
32
Chapter 4. Radiation Modelling
4 3 2 1 0 3
4 3 2 1 0 3
3 0 (a)
2
X1 ax 2 is
1
ax is
1
0
Y
X1 ax 2 is
ax is
2
Y
0
Z axis
5
Z axis
5
3 0 (b)
Figure 4.5: Geometry layout for direction vector used in direct exchange area verification. Direction (X5,Y5,Z5) is shown here with axisymmetric input data was calculated; the details of the first geometry are given in Fig. 4.7.
4.4.1
Case 1
Stuttaford [46] reported Monte Carlo results for net radiative heat flux along the axial length of the cylindrical combustor (Case 1). Since the input data is axisymmetric the resultant radiative heat flux is also axisymmetric. There are three main parameters in the current implementation of zonal method: (1) major grid panels which determine the number of isothermal zones to represent input data and hence, overall accuracy of the solution; (2) path interval length which determines the accuracy of transmissivity factor calculation from one panel to the other; and (3) the grid density (divisions) of each panel which determine the accuracy of the DEA integral calculated. These parameters were varied and solved for the geometry shown in Fig. 4.7; the results are shown in Fig. 4.8. For the results obtained in Fig. 4.8(a), the path length was kept constant at 0.005 m and 2 axial, radial and circumferential divisions per panel; the grid number shown in the figure corresponds to axial and radial divisions respectively. The execution time varied from 36 seconds for the (10,2) grid to 3 hours for the (20,8) grid. For the results obtained in Fig. 4.8(b), the major grid was kept constant with 20 axial and 4 radial divisions and
33
4.4. Results and Discussion
A
B
A
B
C
C
Figure 4.6: Comparison of discretization schemes that can be used for a cylindrical geometry. On the left circumferential divisions vary with radial position, on the right they are constant 5m 2.5 m
Twall = 500 K, = 0.8
1m 0.5 m
Tgas = 1700 K κ = 0.6 m−1
Tgas = 1100 K κ = 0.05 m−1
Figure 4.7: Cylindrical geometry to represent a combustor [46]. each panel is resolved with 2 points in each direction. The execution time varied from 7 minutes for the path length of 0.01 m to 54 minutes for the path length of 0.001 m. For the results obtained in Fig. 4.8(c), the major grid was kept constant with 10 axial and 2 divisions with the path length of 0.05 m. The execution time varied from 1 minute for 4 divisions to 6 minutes for 6 divisions per panel. Overall, the predicted solutions follow the general trend of Monte Carlo solution, however, results under predicts radiative flux by 20% at peak conditions. Overall impact of parameter variation on the final results is minor; out of the three parameters, path interval length had the most effect, followed by major grid divisions and grid density for each panel. While the offset of the predicted and measured results cannot fully be explained and might be a limitation of the zonal method, the small variation in results with variation in input parameters can be explained due to the nature of the input data. In the geometry under consideration, there are only two distinct zones where the temperatures and absorption coefficient vary while the temperature and emissivity of the
34
Chapter 4. Radiation Modelling
Net Radiative Heat Flux % ( qpeak,q
100
Zonal (20,4) Zonal (20,8)
80 60 40 20 0 0.0
0.2
0.4 0.6 Axial Location ( lol )
0.8
Monte Carlo Zonal ∆S =0.01 Zonal ∆S =0.005 Zonal ∆S =0.001
MC
Monte Carlo Case 1 Zonal (10,2)
MC
Net Radiative Heat Flux % ( qpeak,q
)
120
)
120
100 80 60 40 20 0 0.0
1.0
0.2
0.4 0.6 Axial Location ( lol )
(a)
0.8
1.0
(b) )
120 Net Radiative Heat Flux % ( qpeak,q
MC
Monte Carlo Case 1 Zonal 4 divs Zonal 5 divs Zonal 6 divs
100 80 60 40 20 0 0.0
0.2
0.4 0.6 Axial Location ( lol )
0.8
1.0
(c)
Figure 4.8: Effect of (a) grid density; (b) path length divisions; and (c) divisions per panel on net radiative heat flux on the wall wall do not vary. Therefore, discretization of the geometry does not have a significant impact on the final solution.
4.4.2
Case 2
Meng¨ u¸c et al. [26] presented numerical solution using P3 approximation technique for radiative heat flux on a cylindrical combustor (Case 2). The input data presented in the paper was used as a second test case and the predicted profile from P3 approximation solution has been presented in subsequent results for this case as reference. The axial temperature and soot profiles that were used for all the conditions are shown in Fig. 4.9(a). Normalized radial profiles that were used for the four subsequent cases are shown in Fig. 4.9(b). The axial profiles were based on experimental data that was gath-
35
4.4. Results and Discussion 2400 2000
3.0 Case 2.1
T (r)−Tw Tm −Tw
20
2.5
Case 2.2
Cs (r) Cs,avg
2.0
Case 2.3
Cs (r) Cs,avg
Case 2.4
T (r)−Tw Cs (r) Tm −Tw & Cs,avg
15
1200
10
Cs
T(K)
g m3
1600
25
1.5 1.0
800 400 0
5
1
2
3
4 z/ro
(a)
5
6
0 7
0.5 0.0 0.0
0.2
0.4
0.6
0.8
1.0
r/ro
(b)
Figure 4.9: (a) Axial temperature and soot profiles and (b) normalized radial profiles used by Meng¨ u¸c et al. [26]
ered previously, while the radial profiles were selected by the authors to represent various conditions that might occur in an actual combustor. A correlation was used to obtain average absorption coefficient based on temperature and local soot concentration. Fig. 4.10 shows results for Case 2.1 where a constant radial soot profile but varying radial temperature profile were used. The three major parameters for zonal method code were varied similar to Case 1. For the results obtained in Fig. 4.10(a), the path length was kept constant at 0.005 m and 2 axial, radial and circumferential divisions per panel. The execution time varied from 29 seconds for (7,3) grid to 1 hour for (22,10) grid. For the results obtained in Fig. 4.10(b), the major grid was kept constant with 15 axial and 7 radial divisions and each panel is resolved with 3 points in each direction. The execution time varied from 27 minutes for path length of 0.01 m to 49 minutes for path length of 0.001 m. For the results obtained in Fig. 4.10(c), the major grid was kept constant with 15 axial and 7 divisions with path length of 0.005 m. The execution time varied from 9 minutes for 2 divisions to 6 minutes for 5 divisions per panel. Similar to Case 1 results, variation of path length interval has the greatest effect on the net radiative heat flux. This effect is greater than Case 1 since there is much more radial and axial variation in absorption coefficient throughout the combustor. The effect of major grid divisions and grid density on net radiative heat flux is extremely low. If the major grid divisions are high in number, which is the case for Fig. 4.10(c), then changing grid density only has a minor effect on DEA accuracy and hence on the final radiative heat flux. Grid density has an effect on DEA accuracy when there is large variation in calculation
36
Chapter 4. Radiation Modelling 120
Zonal (7,3) Zonal (8,4)
96
Zonal (22,10)
72 48 24 0 0
1
2
3 4 5 Axial Location ( rl0 )
6
P3 Approx. Case 2.1
3
Net Radiative Heat Flux % ( qpeak,q P )
P3 Approx. Case 2.1
3
Net Radiative Heat Flux % ( qpeak,q P )
120
Zonal ∆S =0.01 Zonal ∆S =0.005 Zonal ∆S =0.001
96 72 48 24 0 0
7
1
2
3 4 5 Axial Location ( rl0 )
(a)
6
7
(b) P3 Approx. Case 2.1
3
Net Radiative Heat Flux % ( qpeak,q P )
120
Zonal Zonal Zonal Zonal
96
2 3 4 5
divs divs divs divs
72 48 24 0 0
1
2
3 4 5 Axial Location ( rl0 )
6
7
(c)
Figure 4.10: Effect of (a) grid density; (b) path length divisions; and (c) divisions per panel on net radiative heat flux on the wall of relative angles (i.e. θi and θj in Eqs. (4.6) and (4.12)) or path length over the elements (S in Eqs. (4.6), (4.9), and (4.12)). For small panels that are relatively far away change in angle and path length for points within the panel is negligible. Furthermore, for cases where large amount of absorbing medium separates the panels the DEA integral is dominated by transmissivity term (τ ). As a result, variation in grid density has almost no effect, while variation in major grid panels has very low effect on the net radiative heat flux on the wall. Effect of varying the soot profile on net radiative heat flux is shown in Fig. 4.11. The trend of axial variation in net radiative flux matches what was obtained with P3 approximation, however, the magnitude of the flux varies significantly. As the soot profile was changed for a higher peak in the center (Case 2.3), Meng¨ u¸c et al. reported an
37
3
120
P3 Approx. Case 2.2 P3 Approx. Case 2.3 Zonal Case 2.2 Zonal Case 2.3
96 72 48 24 0 0
Net Radiative Heat Flux % ( Case2.1 qqpeak, P )
3
Net Radiative Heat Flux % ( Case2.2 qqpeak, P )
4.4. Results and Discussion
1
2
3 4 5 Axial Location ( rl0 )
(a)
6
7
450
P3 Approx. Case 2.1 P3 Approx. Case 2.4 Zonal Case 2.1 Zonal Case 2.4
360 270 180 90 0 0
1
2
3 4 5 Axial Location ( rl0 )
6
7
(b)
Figure 4.11: Effect of (a) radial soot profile variation and (b) constant radial temperature and soot profile on net radiative heat flux on the wall increase in overall radiative flux and the reason for this change is explained to be lower absorption of radiation near the wall due to lower soot. Results from the zonal method show the opposite trend where the radiative heat flux decreased as the amount of soot near the wall is lowered; shown in Fig. 4.11(a). For the final Case 2.4 the radial soot and temperature profile were kept at a constant average value; results are shown in Fig. 4.11(b). Results obtained by the authors from P3 approximation show an increase in net radiative heat flux as the temperature near the wall is increased from Case 2.1 to Case 2.4 (80% increase in peak flux) and that trend is matched by results obtained from zonal method although the increase in peak flux is much higher at 300%. By comparing results from the two methods for all four cases it can be seen that the current implementation of zonal method is much more sensitive to absorption coefficient near the wall.
4.4.3
Case 3
Kayakol et al. [33] compared measured results for incident radiative heat flux on a representative gas turbine combustor with predicted results from DOM simulation. Combustor length to radius ratio was approximately 8.1 therefore, it is comparable to the geometry used by Meng¨ u¸c et al. Other details regarding the geometry of the combustor and input parameters are described in detail by Kayakol et al. [33]. Absorption coefficient for thirteen radial and fourteen axial positions were measured, therefore, the major grid used
38 Incident Radiative Heat Flux % ( qpeak,q meas )
Chapter 4. Radiation Modelling 120 96 72 48
Measured Zonal Method
24 0 0.0
0.2
0.4 0.6 Axial Location (l/l0)
0.8
1.0
Figure 4.12: Measured vs predicted results for incident radiative heat flux for zonal method matched that configuration with 3 divisions per panel (each direction) and path length interval equal to 0.005 m. Since the major grid was extremely dense the case finished execution in 5 hours. Figure 4.12 shows a comparison of measured and predicted result from zonal method. Results show a trend that matches measured results, with heat flux under predicted in the beginning of the combustor but the location and magnitude of the peak matching measured results. So far the discussion has focused on radiation in cylindrical enclosures, however, majority of the implementation remains the same for implementation on a geometry representative of a gas turbine combustor. An example geometry to model a gas turbine combustor is shown in Figure 4.13, where the primary and secondary zones for one sector are modelled ignoring rest of the combustor geometry. Curvature of the combustor requires that surface and gaseous elements that cannot receive any radiation from other zones (due to the combustor core being in the way) be recognized for the computation, these elements were ignored in the shown geometry. Since the combustor is the hottest in the primary zone, this assumption might still produce results that are representative of actual radiation flux. Another geometry feature that is being ignored in Figure 4.13 is various adjacent combustor sectors that would increase the radiative flux on the current sector. If it is assumed that all nozzles in the combustor produce the same pattern and by extension the same temperature profile, the radiative flux would be axisymmetric in nature, however,
39
4.4. Results and Discussion
A 0.2
B H
Z axis
0.1 0.0
− 0.1
X
− 0.2 − 0.1 0.0 0.1 0.20.0
C G
is ax
0.1
0.2 Y axis
0.3
− 0.2 0.4
Figure 4.13: A course discretization grid of primary and secondary zone of a reverse flow combustor for zonal method implementation
D F E Figure 4.14: Combustor cross section schematic with eight sectors
the interaction of various sectors would need to be taken into account. Figure 4.14 shows a combustor cross section with eight sectors. Application of axisymmetric condition was discussed in the context of cylindrical geometry with three major sub grids; this can be further extended to any number of subgrids. Therefore, for eight sectors the computation time would be eight times the time required for only a single sector. This is the worst case scenario as in reality the combustor core blocks a significant number of sectors. In case this strategy is computationally expensive, a scaling factor could be applied to the results of a single sector as the trend for incident radiative flux in axial direction should not vary significantly from single sector results.
Chapter 5 Combustion Modelling 5.1
Background
There are three main characteristic times that describe combustion within the combustor, these are, chemical, evaporation and mixing characteristic times. In simulations employing CFD an assumption that is usually made is that the mixing or evaporation rate is the limiting factor while the chemical rate is infinitely fast. In reactor modelling the assumption is made that chemical rate is the limiting factor while high turbulence causes mixing to be instantaneous. Conceptual reactors make broad assumptions about the flow and if used individually can seldom model real flows, but in combination might be able to model complex flows. Several theoretical reactors are described by Turns [49] but in the context of combustion chambers that operate on almost constant pressure two are the most relevant which are: (1) Perfectly Stirred Reactor (PSR); and (2) Plug Flow Reactor (PFR). A PSR is an ideal reactor in which perfect, infinitely fast mixing is achieved inside a control volume [49]. Mass conservation is given by dmi,cv =m ˙ 000 V + m ˙ i,in − m ˙ out . dt
(5.1)
For a non compressible system the left hand side of the equation reduces to zero. Mass generation can be given by each reactants chemical production rate changing the conservation equation to ωM ˙ Wi V + m(Y ˙ i,in − Yi,out ) = 0, (5.2) where Y defines the mass fraction of each species. Energy balance can be done based on 40
41
5.2. Cantera Model internal enthalpy of incoming reactants and outgoing products to give Q = m(i ˙ out − iin ).
(5.3)
To close the system of equations equation of state is used to relate density, pressure and temperature as follows: P M Wmix (5.4) ρ= Ru Tr Finally, residence time of the reactor can be defined as tR =
ρV m ˙
(5.5)
Regions in the combustor where the turbulence is high, for example the primary zone, might approach the idealisation of a PSR reactor. A PFR is an ideal reactor which is one dimensional in nature (compared to zero dimensional for PSR) for which velocity, concentration, temperature and other properties describing the state might vary axially. It also has the following attributes: (1) steady state, steady flow, (2) no mixing in axial direction, (3) uniform properties in perpendicular direction to the flow, (4) ideal frictionless flow and (5) ideal gas behaviour [49]. Regions in the combustor where the flow is one dimensional and has relatively low turbulence, for example the dilution zone, can be modelled with a PFR. These conceptual reactors might be linked in a network with output from one being fed as the input for the other to resemble a complex combustor flow. Ideally, the closer the reactor would be to actual flow features, the better the accuracy of species production and temperature data.
5.2
Cantera Model
For the current project an open source software, CANTERA [5] was used to solve chemical equilibrium equations and model reactor network. CANTERA numerically solves differential and algebraic equations discussed earlier; the software has the capability to model various forms of reactors, namely, a general reactor, ideal gas reactor, constant pressure reactor, ideal gas constant pressure reactor, and flow reactor. Differences in various models lie in various assumptions being made to solve the mass, momentum (in the case of flow reactor i.e. a PFR) species conservation and energy equations being solved. To match the equations discussed in previous section ideal gas constant pressure reactor and flow reactor were selected for the current study. CANTERA inherently solves the
42
Chapter 5. Combustion Modelling
time dependent versions of equations discussed previously, as a result, to obtain steady state results all cantera simulations were run until temperature had converged.
5.3
Implementation
Since the reactor networks would eventually be part of a larger combustor design system, the creation of the network was generalized. The generalization would allow for a new reactor network to be created based on the combustor geometry. The user can specify the interconnectedness of the reactors via an input text file. Flow splits from one reactor to every other reactor in the network are stored in a matrix structure (this is similar to an adjacency matrix for a graph which would store a 1 or 0 but instead stores the split value). Details of the amount of fuel and air to each reactor are stored and read as a table in another input text file. Flow splits are initially defined as percentage split of the input flow; this was done to make the creation of input files easier. To get the actual mass flow rate from the percentage is trivial if the network does not contain any recirculation zones or any part of the flow coming back to the input. Since recirculation flow are to be expected an iterative algorithm was used to find the absolute flow rate. First, each reactor’s output is calculated based on the flow split and the input while ignoring any recirculation. In second iteration previous results are used to calculate new outputs while taking recirculation into account. This is continued until the mass flow rates in the entire network converge. Recirculation also poses a problem when it comes to solving the networks in CANTERA directly. In cases where flow from downstream reactors is connected to upstream reactors the CANTERA numerical solver fails to converge. If the evaporation and heat loss are being solved iteratively it implies that the entire network would be solved a large number of times. Since CANTERA inherently solves for time dependent system, a highly stiff system would be computationally costly to solve in an iterative scheme. In order to achieve convergence faster each reactor was ignited and solved separately, and then the network solved iteratively similar to the approach used when calculating mass flow rates from flow splits. CANTERA utilizes mechanism files that contain list of elementary reactions along with reaction rate data and thermodynamic data to predict chemical compositions. Marchand [27] and Lanewala [18] had success in estimating emissions with CANTERA and utilizing GRI3.0 mechanism. GRI3.0 mechanism is designed for detailed combustion analysis of propane and methane,
43
5.4. Results and Discussion
(a)
(b)
Figure 5.1: Engine A and Engine B profile schematics for comparison contains 53 species and 325 reactions [8]. Since propane has higher heat value than jet fuel the reactor network calculations were adjusted by decreasing the total amount of fuel entering each reactor in proportion to the fraction of higher heating value of the two fuels [27, 18] . This method ensures that the total energy entering the system is similar to what it would be with the jet fuel. Jet fuel mechanism by Westbrook was tested by Marchand but not used for the final calculation as it was found significantly more computationally expensive (it contains 1421 species and 7851 reactions) and not suitable for preliminary design calculations [27]. In order to model temperature more accurately two additional reduced mechanisms were tested; these were Kollrack’s jet fuel mechanism [11, 28] and University of California San Diego’s (UCSD) JP 10 mechanism [29]. To ignite each reactor in the reactor network, a short Gaussian pulse of hot gas is initially injected into each reactor. The pulse consists of combustion products of the specified fuel with air at 1000K. This ensures that no concentration of species introduced is somewhat comparable to what is expected at the end from the reactor.
5.4
Results and Discussion
Figs. 5.1(a) and 5.1(b) shows profiles for Engine A and Engine B that were used as test cases discussed in this section. Fig. 5.2 shows results that were obtained for Engine A using a simple linear reactor network of five reactors with GRI 3.0 mechanism and gaseous propane as fuel. The reactor network configuration shown in Fig. 5.3 was run for idle, approach, climb and take off power settings. Operating conditions and the combustor outlet measured temperature were obtained from a technical report by Gratton et al.
44
Chapter 5. Combustion Modelling 2.5 Combustor Temperature Rise,
T4 T3
Measured Predicted
2.0 1.5 1.0 0.5 0.0
Idle
Approach
Climb
Take Off
Figure 5.2: Engine A combustor outlet temperature; measured and predicted results
1
2
3
4
5
Figure 5.3: Engine A reactor network [24] . Since the simulation was run for propane it is surprising that the results are within 6% of experimental data. This might be due to higher residence time and lower fuel to air ratio that is used in relatively older combustors. If the residence time is higher then the temperature prediction would be controlled by elementary reactions that take place later in the combustion process, which are similar for hydrocarbon fuels. In Figure 5.4 combustor outlet temperature (T4 ) is plotted against air to fuel ratio (AFR) for Engine B. As compared to Engine A, Engine B is of a modern design with the overall residence time being much lower and the fuel to air ratio being higher. Operating conditions and combustor outlet measured temperature for engine B was obtained from a book by Schumann [43]. Fig. 5.4(a) shows the effect of various reduced mechanisms on outlet temperature. Out of the three mechanisms tested GRI 3.0 and UCSD mechanisms follow the overall trend with UCSD being closer to experimental results even though it models JP 10 fuel and not Jet A while Kollrack mechanism severely under-predicts T4 . Fig. 5.4(b) shows the effect of network configuration on temperature prediction. Configuration 1 shown in the figure consists of 3 linear reactors similar to the linear reactor configuration of Engine A in Fig. 5.3. Configuration 2 was a network that attempted to improve the modelling of recirculation zone by sending a percentage output of reactor 2 to reactor 3 and then to 1 as shown in Fig. 5.5(a). The new model did not have any impact on the final outlet temperature prediction. However, network configuration does
45
5.4. Results and Discussion 2.50
2.50
2.25
T4 T3
Combustor Temperature Rise,
Combustor Temperature Rise,
T4 T3
Measured GRI 3.0 Kollrack UCSD
2.00
1.75
1.50 35
45
55
65
2.25
2.00
1.75
1.50 35
GRI 3.0 Config 1 GRI 3.0 Config 2 UCSD Config 1 UCSD Config 2
45
55
AFR
(a)
65
AFR
(b)
Figure 5.4: Effect of (a) mechanisms and (b) reactor network setup on combustor outlet temperature for Engine B
have an impact on the local temperature in distinct combustor zones. Local temperature is also an important factor for radiation prediction and for emission prediction (especially NOx ). The effect of network configuration on NOx prediction is shown in Figure 5.6. Here the third configuration of reactor networks tested is shown in Fig. 5.5(b). Detailed analysis on emission prediction with network reactors was conducted by Marchand and Lanewala and is not the focus of current project, however, since heat transfer prediction and emission prediction would be part of a single PMDO design tool it is important to note that the current tool gives the user the ability to predict emissions while the configuration can be set based on a strategy outlined in aforementioned projects.
46
Chapter 5. Combustion Modelling
1
2
4
3
1
2
4
3 (a)
(b)
Figure 5.5: Network configurations for Engine B: (a) configuration 2 and (b) configuration 3
Ox EIN Ox % ( EIN OEIN x, Config
1 max
)
1.0 0.8
Config 1 Config 2 Config 3
0.6 0.4 0.2 0.0 1.8
2.0 2.2 2.4 Combustor Temperature Rise, TT43
2.6
Figure 5.6: Effect of reactor network setup on combustor N Ox emission index (EIN Ox )
Chapter 6 Cooling Technology Modelling 6.1
Film Cooling
Film cooling regardless of the geometry being used has mainly been modelled using correlations. The definition of cooling film effectiveness (Eq. (2.5)) and the consequent calculations for liner temperature are described in detail in Section 2.3. Correlations for louvres fall in two broad categories as described by Lefebvre [20]. There are correlations based on turbulent boundary-layer model and wall-jet model; Lefebvre’s correlations for these models are: x 0.3 M µ −0.15 c η = 0.6 Re (6.1) Ms µhg and
η = 1.28
µc µhg
0.15 −0.2 x −0.2 t . s s
(6.2)
Boundary-layer models are derived based on idealized turbulent boundary layer downstream of the slot and might not represent conditions close to the slot; Eq. (6.1) was based on measured skin friction coefficient to avoid this drawback. Lefebvre used the following correlation for machined ring cooling slot: η = 1.0 − 0.12SN0.65 ,
(6.3)
and the following correlation for stacked ring cooling slots: η = 1.0 − 0.094SN0.65 , 47
(6.4)
48 where
Chapter 6. Cooling Technology Modelling
x − xp SN = Ms
µc Re µhg
−0.15
Ao . Aeff
Li and Mongia [22] employed mass transfer analogy to measure film cooling effectiveness in machined ring liners. Computational results were gathered for certain design of experiment configurations and on the basis of experimental and numerical results correlation with the functional form shown below was suggested: η = e−aζ
b
(6.5)
where ζ is defined as x = ζ= Ms
µinf µ2
Rex Re2
(6.6)
and a and b are functions of blowing ratio (m), slot lip thickness to slot height ratio (s/t), coolant injection angle, lip taper angle and starting edge angle on nugget-slot exit. Lefebvre’s Eqs. (6.2) and (6.3) were also compared with the collected data. Equation (6.2) was found not be a good match while equation (6.3) matched data to within 15% accuracy. However, the new correlation provides predictive capabilities with greater number of geometrical parameters. Some authors have found film mass flow rate and hot stream turbulence level to be of greater importance than slot geometry; based on this Juhasz and Marek [15] provided a correlation of the form given by η=
1 1 + Cm Mxs
(6.7)
and for the case where the cold and hot air differ significantly in composition the correlation included heat capacity terms determined using η=
1 C
1 + Cm Mxs Cpph
.
(6.8)
s
The proposed correlation was able to predict test results within 20%. Other correlations for various geometries take into account different geometric factors; However, blowing ratio, local Reynolds number, and downstream distance are common correlating parameters. Various empirical correlations have been derived for different flow conditions (Nusselt number formulas) and cooling geometry (film effectiveness formulas). These parameters are the key in order to optimize for minimum required flux. Various conditions under which these equations are valid must be taken into account in order to
49
6.2. Double Walled Cooling
L Tc,2
3
Tw,2
2
4 Tw,1
1 Tc,i
Figure 6.1: Schematic representing Counter Flow Film Cooling (CFFC) arrive at a realizable and valid solution.
6.2
Double Walled Cooling
The analytical model to solve for Counter Flow Film Cooling (CFFC) and Parallel Flow Film Cooling (PFFC) was obtained from the NASA technical report by Colladay [4]. Fig. 6.1 shows the geometry setup used for the CFFC solver. The main CFFC tile of length L is shown in striped pattern while the bottom surface shown in gray does not take part in heat transfer. Heat gained through convection by the coolant as it passes by the underside of the tile (labelled ‘1’ in the figure) is given by dQc = hcf (Tw2 − Tc ) dAeff .
(6.9)
A Coulbourn j factor correlation appropriate for fin geometry was used to obtain heat transfer coefficient for the cool side. This correlation was taken from a paper by London [23] and is given by m ˙c hcf = Cp,c j Pr−2/3 . (6.10) Ap The total area for the transfer is also increased and reflected by the effective area term (Aeff ) in Eq. (6.9). Energy transfer to the coolant is also limited by the heat capacity of the coolant given by the following equation: dQc = m˙ c Cp,c dTc .
(6.11)
Similarly, on the hot side of the tile (labelled ‘3’) the convective heat transfer is described
50
Chapter 6. Cooling Technology Modelling
by dQhg = hhg (Taw − Tw1 ) dA.
(6.12)
The value for hot side heat transfer coefficient is given by the following Nusselt number correlation: hhg x = 0.0296 Rex0.8 Pr1/3 (6.13) Nu = k To obtain the adiabatic wall temperature (Taw ) in Eq. (6.12), the following correlation for film effectiveness was used: x −8/10 η = 16.9 Ms
(6.14)
Note that film effectiveness was defined in equation Eq. (2.5), here Tc,2 is taken as the coolant temperature (Tc in the original definition). This equation was obtained from [12]. A two-dimensional conduction code was written for the tile itself using finite difference technique. The heat equation for two-dimensional, steady state, with no generation and constant thermal conductivity that was used for the code is given by ∂ 2T ∂ 2T + = 0. ∂x2 ∂y 2
(6.15)
The finite difference scheme was taken from textbook by Bergman et al. [47]. The boundary condition for the finite difference code were heat flux for the top and the bottom sides of the tile (labelled ‘3’ and ‘1’ respectively) while the other edges were assumed to be at constant temperature, in this case being Tc,2 . As the resulting system of equations is coupled it was solved for iteratively; the program logic that was used can be found in the original report and will be described here briefly by the following steps: (1) take an initial guess for coolant temperature profile; (2) calculate cold side heat transfer coefficient and fin effectiveness; (3) calculate film effectiveness and hot side heat transfer coefficient; (4) run the two-dimensional conduction code for the tile; and (5) repeat steps (1) - (4) until the solution converges. To obtain a stable solution the coolant temperature profile was calculated until it had converged before the conduction code was run, this also ensured that less time was spent on conduction code which relatively takes longer to execute. Fig. 6.2 shows the results for temperature profiles that can be expected on the top surfaces of tiles used in CFFC and PFFC. The coolant mass flow rate was adjusted by Colladay [4] to achieve a maximum surface temperature of 1255 K. In the results presented here maximum temperatures for all the lengths except 2.5 and 5 cm are higher
51
6.2. Double Walled Cooling 1.6
1.7
Temperature( TTw2 ) c,i
1.5 1.4 1.3 1.2 1.1 0.00
1.6 Temperature( TTw2 ) c,i
25 cm 20 cm 15 cm 10 cm 5 cm 2.5 cm
1.5 1.4 1.3
25 cm 20 cm 15 cm 10 cm 5 cm 2.5 cm
1.2 1.1
0.05
0.10 0.15 Distance(m)
(a)
0.20
0.25
1.0 0.00
0.05
0.10 0.15 Distance(m)
0.20
0.25
(b)
Figure 6.2: Temperature profiles predicted for panels of various lengths for (a) CFFC and (b) PFFC. than 1255. This discrepancy is due to ambiguity in certain input parameters that could not be clarified. However, overall the trends shown here match what was presented in the original report. As discussed earlier determining and minimizing the coolant required for a combustor is an important task for the designers. Colladay found that for tiles with longer length (larger area) the total coolant required to maintain a certain temperature was larger, however, the coolant required per area was lower [4]. Therefore, longer tiles were more efficient in terms of coolant requirements. Increasing the length of the tiles increases the time that the tile has to transfer energy to the coolant before it exits on top to form the film layer. This also results in the temperature of the tile near the slot being higher for longer section lengths. A note should be made here that in his analysis Colladay predicted the required pressure for the required mass flow rate, which increased for larger tiles [4]. In reality this maximum pressure drop across the combustor cannot be changed by the combustor designer, however, this analysis suggests that the longest tile that can sustain the mass flow rate according to a given pressure should be picked. Effect of wall thickness and operating pressure on overall heat transfer was also reported by Colladay. A thin tile wall was reported to have higher temperature on its coolant side which allow for higher heat flux transfer to the coolant, resulting in more efficient cooling. Higher operating pressure had the effect of increasing coolant density and hence, the total coolant mass flow rate while the heat flux increase to the tile was Tc,2 −Tc,i ) not as high. As a result, the overall convective effectiveness (defined by η = Tw,2max −Tc,i
52
Chapter 6. Cooling Technology Modelling 1.8
1.6
1.4
1.2
1.0 0.00
0.05
0.10 0.15 Distance(m)
0.20
Tw2, peak Tc,i )
100 300 500 1000 2000 3000
Temperature(
Temperature(
Tw2, peak Tc,i )
1.8
1.6
1.4
1.2
1.0 2 10
0.25
(a)
104
1.8
1.45 1.40
Tw2, peak Tc,i )
1.50
Tpeak Tavg
Temperature(
2 4 8 12 20 100
1.55 Tw2, peak Tc,i )
103 Number of Nodes(#)
(b)
1.60
Temperature(
Tpeak Tavg
1.6
1.4
1.2
1.35 1.30 0.000
0.005
0.010 0.015 0.020 Distance (m)
0.025
0.030
(c)
1.0 0 10
101 Number of Nodes (#)
102
(d)
Figure 6.3: Effect of number of axial nodes on a 25 cm section length CFFC’s temperature profile is shown in (a) while the change in peak and average temperature is shown in (b). Effect of number of normal nodes on a 2.5 cm section length CFFC’s temperature profile is shown in (c) while the change in peak and average temperature is shown in (d).
decreased with an increase in pressure [4]. In order to determine the density of the grid used for the conduction solver a number of cases were run varying axial and normal nodes; the results for these cases are shown in Fig. 6.3. Axial nodes were varied from 100 to 3000 for a CFFC panel with length of 25 cm. The peak of the temperature profile decreased with an increase in axial nodes with not much improvement for higher than 1000 nodes as shown in Fig. 6.3(a) and Fig. 6.3(b). Since computational cost was not a major issue, highest density of axial nodes was selected for all the analysis shown here (120 nodes per cm). Number of normal nodes did not have a significant impact on the temperature profile.
53
6.3. Effusion and Transpiration Cooling Hot Gas, Tm Twm
Coolant Direction Text
Tent Twc Coolant Gas, Tc (a)
(b)
Figure 6.4: Schematic representing effusion geometry
6.3
Effusion and Transpiration Cooling
Fig. 6.4 shows a schematic representation for the geometry that was solved for by Martiny et al. [25]. The analytical model consists of one cooling hole (from an array of holes that might be used to cover a large area shown in Fig. 6.4(a)). It is assumed that property variation over the stream-wise and pitch-wise direction is negligible and due to symmetry in both directions all planes can be considered adiabatic. Therefore, under this assumption the heat transfer is only in the normal direction. For an infinitesimal wall and air element shown in Fig. 6.4(b), equation representing the heat balance is given by dQw dQa + = 0, (6.16) dz dz w a where dQ represents heat lost by the wall and dQ represents heat gained by the air. dz dz Heat lost by the wall is due to conduction and given by
dQw d2 Tw = −kw A 2 , dz dz
(6.17)
while the heat gained by the air is controlled by the total heat capacity and convective
54
Chapter 6. Cooling Technology Modelling
heat transfer. These relationships are given by the following equations: dQa dTa = m˙ a Cp ∗ ∗ dz dz dTa hh πd (Tw − Ta ) = ∗ dz m˙ a Cp
(6.18) (6.19)
This is a system of coupled equation and can be solved analytically to give equations for temperature profile for the wall and for temperature profile for the air. The solution to the resulting two coupled ODE’s defining the cooling in the effusion hole are m˙ a Cp Tw = sin α kw A
C 1 β1 z ∗ C 2 β2 z ∗ e + e β1 β2
and ∗
∗
Ta = C1 eβ1 z + C2 eβ2 z −
+ C4
γ1 C3 , sin αγ2
(6.20)
(6.21)
where, hh πd , m˙ a Cp hh πd , γ2 = kw A γ1 =
β1,2
γ1 =− ± 2
r
γ12 + sin αγ2 . 4
The four constants of the two differential equations can be found with equations that describe the boundary conditions of the system. The total enthalpy rise of the coolant has to equal to the convective heat flux due to hot gases; this relationship is given by the following equation: hhg A(Thg − Tw1 ) − m˙ a Cp (Th1 − Tc ) = 0 (6.22) Convective heat flux of the coolant on the cold side is equal to the enthalpy rise in the coolant before it enters the hole; this relationship is given by the following equation: hc A(Tw2 − Th2 ) − m˙ a Cp (Th2 − Tc ) = 0
(6.23)
The convective heat flux of the coolant equals the conductive heat flux on the cold side; this relationship is given by the following equation: dTw hc A(Tw2 − Tent ) − kw A =0 dz z=0
(6.24)
55
c Non Dimensional Wall Temperature, Θw ( TTmw −T −Tc )
6.3. Effusion and Transpiration Cooling 0.85
0.290
Cooling Effectiveness, η
Predicted
0.285 0.280 0.275 0.270 0.265 0.0
0.80 0.75 0.70 0.65
Predicted 0.2
0.4 0.6 Normal Location (z/t)
(a)
0.8
1.0
0.60 0
20
40 60 Incident Angle (◦)
80
100
(b)
Figure 6.5: Predicted results from the effusion model for (a) wall temperature distribution and (b) effect of blowing angle on cooling effectiveness Eq. (6.19) can be evaluated for z = 0 to get a relation between Twc and Tent . With these four equations the formulation of the analytical model is complete. Martiny et al. presented numerous trends with the analytical model and these were all matched to verify the implementation. Two of these predicted results are shown in Fig. 6.5. Fig. 6.5(a) shows the trend for non dimensional wall temperature; one of the benefits of the model is that it is able to capture the non linear behaviour of this profile due to coupling of conduction and convection. Fig. 6.5(b) shows the effect of blowing angle on cooling effectiveness; this matches the expected trend discussed in Section 2.1 where with shallower angles the length of the hole increases allowing for greater heat transfer to take place from the wall to the coolant and increasing the cooling effectiveness. To further demonstrate the capabilities of this model for a designer the formulation of the problem was changed to produce the required mass flow rate per area for a specified hot side temperature. Implementation of this model as a module gives the designer the opportunity to optimize the cooling scheme geometry for a given objective function (in this case mass flow rate per area). DAKOTA [2] is one such open source software that is available for optimization and parametric analysis (amongst other features). Parametric study was done on the effusion cooling model using data presented in the original paper, the required hot side temperature for the cooling surface was set at 600K. Further assumption was made the the pitch in stream-wise and pitch-wise direction are the same. With this assumption there are three main geometrical parameters that can be varied with this model: (1) hole diameter; (2) blowing angle; and (3) liner thickness. During combustor design phase liner thickness is most likely to be fixed so it was not
56 2.50
8.000
Diameter (mm)
2.04 1.58
2.0 00
7.000
0
0
4.0
00
3.0 1.12
5.0 0
6.00
00
0.66
1.000 0.20 10
30
50 Incident angle (◦)
70
90
Figure 6.6: Parametric study for effusion cooling.
˜ Non Dimensional Wall Temperature, Θ, Θ
Chapter 6. Cooling Technology Modelling 0.7 0.6
Θs Θf
0.5
˜s Θs , Θ ˜f Θf , Θ
0.4 0.3 0.2 0.1 0.0 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Normal Location, X/m
Figure 6.7: Predicted results from transpiration model for wall temperature distribution with and without a ceramic coating.
considered in this analysis. The contour plot from the parametric study based on variation of hole diameter and blowing angle is shown in Fig. 6.6. The effect of variation in blowing angle is more pronounced for larger hole diameters than smaller ones. The accuracy of the predictions cannot be evaluated without experimental data relevant to actual gas turbine combustor conditions, however, this demonstrates that models such as the one presented here can provide more physical insight than a semi-empirical correlation while allowing for fast execution time, hence being relevant in preliminary design environment. Furthermore, the effect of manufacturing imperfections on cooling performance can be tested through coupling of cooling modules with software such as DAKOTA that allow uncertainty analysis. Wolfersdof [51] compared effect of various cold side boundary conditions on an analytical model for transpiration cooling. One-dimensional analytical model for transpiration cooling presented in the paper and a number of other studies [10] is very similar to the one presented for effusion cooling, hence, the details would not be presented here. Coupled relationship between conduction of the wall and convection of the air remains the same, however, the convection relationship was modelled using volumetric heat transfer coefficient rather than a heat transfer coefficient for a discrete hole. Results for nondimensional wall and coolant temperature can be seen in Fig. 6.7 along with a case where a ceramic coating is added.
Chapter 7 Conclusions 7.1
Thesis Accomplishments
The objective of this thesis, as stated in Chapter 1, was to assess heat transfer analysis tools that can be utilized to predict gas turbine combustion chamber liner temperature in the context of preliminary design. Liner temperature prediction calculation was subdivided into three steps (or modules): (1) hot gas radiation prediction; (2) hot gas temperature prediction; and (3) cooling technology modelling. Although currently the modules operate independently, their integration in larger PMDO tool could be possible, hence, factors effecting this integration were also taken into account. In this thesis, a zonal method was assessed as a radiation prediction tool for gas turbine combustors. Three cases were tested for cylindrical geometry with given radiative properties. Trends produced by the zonal method for variation in net radiative flux on the wall were similar to the ones obtained through higher fidelity numerical techniques in all cases. The zonal method under-predicted radiative flux for the cases where the absorption coefficient and temperature near the wall were low (case 1, 2.1, 2.2 and 2.3) and over-predicted for the case where these properties were high near the wall (case 2.4). The method also demonstrated to be insensitive to various program parameters that were varied. This can be useful as accurate results can be achieved with low execution time. Reactor network analysis was assessed as a combustor hot gas temperature prediction tool. Comparison was made between predicted and measured combustor outlet temperatures for two combustors, with accurate predictions for the combustor with higher residence time and valid trends for the combustor with lower residence time. The models were able to predict the general trend of outlet temperature variation with air to fuel ratio and the drop in local temperature with axial distance. Hence, reactor networks can be used as a tool to assess the difference in average temperature in local sections 57
58
Chapter 7. Conclusions
of the combustor such as primary and secondary zones based on the amount of cooling flow that is provided. Furthermore, a programming environment was created to give the designer the ability to rapidly create reactor networks with various structures; this can be a valuable tool when assessing the configuration of a network for new combustors. Various methods that are available in literature for cooling technology modelling were discussed in this study. Models for effusion/ transpiration cooling have the ability to provide temperature profiles for the wall and coolant air in normal direction while the model for CFFC/PFFC has the ability to predict axial and normal temperature variation on a cooling tile. Both of the aforementioned models make broad assumptions, however, they can provide more physical insight as compared to empirical correlations. It is also possible for the designer to asses the impact of geometrical parameters on cooling performance which might not be possible with correlations due to complicated geometry (for example CFFC) or such detailed correlations being unavailable in literature.
7.2
Future Work
Although preliminary tests were conducted on each sub module as discussed in the main body of the thesis, comparison should be made between relevant combustor geometry and operating conditions for further validation. Integration of the three sub modules would also be a future objective following this validation. The combustion modelling can also be improved through implementation of evaporation models and detailed network configurations that allow for emission prediction and flame stability. This would allow the temperature prediction tool to be integrated with a higher level preliminary design tool. The zonal method should be implemented on a full annular combustor geometry to assess the effect of adjacent sectors on the radiation prediction. Third party software can be used for generation of parametric geometry which can be linked with numerical solver for zonal method, reducing the programming overhead required for each new geometry. The models for cooling technology lack the capability to predict mass flow rate based on pressure drop and this relationship should be added for each model for greater accuracy and proper integration.
Bibliography [1] A. Picchi L. Tarchi A. Andreini, B. Facchini and F. Turrini. Experimental and Theoretical Investigation of Thermal Effectiveness in Multiperforated Plates for Combustor Liner Effusion Cooling. Journal of Turbomachinery, 136(9):091003–1:13, 2014. [2] W. J. Bohnhoff K. R. Dalbey M. S. Ebeida J. P. Eddy M. S. Eldred P. D. Hough K. T. Hu J. D. Jakeman L. P. Swiler B. M. Adams, L. E. Bauman and D. M. Vigil. Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: Version 5. User’s Manual SAND2010-2183, Sandia National Laboratories, apr 2013. [3] S. Bradshaw and I. Waitz. Impact of Manufacturing Variability on Combustor Liner Durability. Journal for Engineering for Gas Turbines and Power, 131(3):032503:1– 12, 2009. [4] R. S. Colladay. Analysis and comparison of wall cooling schemes for advanced gas turbine applications. Technical Report NASA-TN-D-6633, NASA, OH, USA, 1972. [5] H. K. Moffat D. G. Goodwin and R. L. Speth. Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes. http:// www.cantera.org, 2013. Version 2.1.0. [6] M. da G. Carvalho and P. J. Coelho. Heat Transfer in Gas Turbine Combustors. Journal of Thermophysics and Heat Transfer, 3(2):123–131, apr 1989. [7] R. I. McCormick S. H. Moustapha P. Sampath E. Hosking, D. P. Kenny and A. A. Smailys. The pw100 engine: 20 years of gas turbine technology evolution. In Design Principles and Methods for Aircraft Gas Turbine Engines, Toulouse, France, may 1998. RTO AVT. [8] M. Frenklach N. W. Moriarty B. Eiteneer M. Goldenberg C. T. Bowman R. K. Hanson S. Song W. C. Gardiner Jr. V. V. Lissianski G. P. Smith, D. M. Golden and Z. Qin. Gri-mech. http://www.me.berkeley.edu/gri_mech/, 2013. Version 3.0. 59
60
Bibliography
[9] S. Howell. Combustor technology for small aircraft engines. In Technology Requirements for Small Gas Turbines. NATO AGARD, oct 1993. [10] R. S. Colladay J. B. Esgar and A. Kaufman. An analysis of capabilities and limitations of film air cooling methods for turbine engines. Technical Report NASA-TND-5992, NASA, OH, USA, 1971. [11] W. H. Heiser J. D. Mattingly and D. T. Pratt. Aircraft Engine Design. AIAA, 2002. [12] R. C. Birkebak J. P. Hartnett and E. R. G. Eckert. Velocity distributions, temperature distributions, effectiveness and heat transfer for air injected through a tangential slot into a turbulent boundary layer. Journal of Heat Transfer, 83(3):293–305, 1961. [13] M. W. Vincent J. Swithenbank, I. Poll and D. D. Wright. Combustion design fundamentals. Symposium (International) on Combustion, 14(1):627 – 638, 1973. [14] D. C. Hammond JR and A. M. Mellor. Analytical Predictions of Emissions from and Within an Allison J-33 Combustor. Combustion Science and Technology, 6(5):279– 286, 1973. [15] A. J. Juhasz and C. J. Marek. Gas turbine combustor liner film cooling for slot geometries in presence of high free stream turbulence. Technical Report NASA-TND-6360, NASA, OH, USA, 1971. [16] D. Kretschmer and J. Odgers. A simple method for the prediction of wall temperatures in gas turbines. In ASME 1978 International Gas Turbine Conference and Products Show. ASME. [17] R. Krewinkel. A review of gas turbine effusion cooling studies. International Journal of Heat and Mass Transfer, 66:706 – 722, 2013. [18] H. Lanewala. Multi-dimensional nitric oxide emissions predictor for preliminary gas turbine combustor design optimization. Master’s thesis, University of Toronto Institute for Aerospace Studies, 2013. [19] A. H. Lefebvre. Flame radiation in gas turbine combustion chambers. International Journal of Heat and Mass Transfer, 27(9):1493 – 1510, 1984. [20] A. H. Lefebvre and D. R. Ballal. Gas Turbine Combustion Alternative Fuels and Emissions Third Edition. CRC Press, FL, USA, 2010.
Bibliography
61
[21] A. H. Lefebvre and M. V. Herbert. Heat-Transfer Processes in Gas-Turbine Combustion Chambers. In Proceedings of the institution of Mechanical Engineers, volume 174, page 463:478. Institution of Mechanical Engineers, 1960. [22] S. C. Li and H. C. Mongia. An improved method for correlation of film-cooling effectiveness of gas turbine combustor liners. In 37th Joint Propulsion Conference and Exhibit, Joint Propulsion Conferences. AIAA, AIAA, 2001. [23] A. L. London and R. K. Shah. Journal of Engineering for Gas Turbines and Power, 90(3):218–228, 1968. [24] I. Critchley M. Gratton and P. Sampath. Alternate fuels combustion research. Technical report, Aero Propulsion Laboratories, Air Force Wright Aeronautical Laboratories, Wright Patterson AFB, OH, 1984. [25] A. Schulz M. Martiny and S. Wittig. Mathematical Model Describing the Coupled Heat Transfer in Effusion Cooled Combusor Walls. In International Gas Turbine & Aeroengine Congress & Exhibition, volume 3, FL, USA, jun 1997. ASME. [26] W. G. Cumming III M. P. Meng¨ uc¸ and R. Viskanta. Radiative transfer in a gasturbine combustor. Journal of Propulsion, 2(3):241–247, 1986. [27] M. Marchand. Multi-dimensional carbon monoxide emissions predictor for preliminary gas turbine combustor design optimization. Master’s thesis, University of Toronto Institute for Aerospace Studies, 2013. [28] D. Mavris. Enhanced emission prediction modeling and analysis for conceptual design. Technical report, Georgia Institute of Technology, 2010. [29] Mechanical and University of California at San Diego Aerospace Engineering (Combustion Research). The San Diego Mechanism, jun 2015. [30] A. M. Mellor. Design of modern gas turbine combustor. Academic Press, 1990. [31] M. F. Modest. Radiative heat transfer. Academic press, San Diego, CA, 2013. [32] A. W. Marshall N. K. Rizk, J. S. Chin and M. K. Razdan. Predictions of NOx Formation Under Combined Droplet and Partially Premixed Reaction of Diffusion Flame Combustors. In ASME 1999 International Gas Turbine and Aeroengine Congress and Exhibition, volume 2, IN, USA, jun. ASME.
62
Bibliography
¨ L. G¨ [33] I. Campbell N. Kayakol, N. Sel¸cuk and O. ulder. Performance of discrete ordinates method in a gas turbine combustor simulator. Experimental Thermal and Fluid Science, 21(13):134 – 141, 2000. [34] A. Sayre N. Lallemant and R. Weber. Evaluation of emissivity correlations for H2O CO2 N2/air mixtures and coupling with solution methods of the radiative transfer equation. Progress in Energy and Combustion Science, 22(6):543 – 574, 1996. [35] Y. S. H. Najjar and R. M. Droubi. Prediction of liner temperature in gas turbine combustors. Fuel, 66(8):1156 – 1160, 1987. [36] C. T. Norgren. Comparison of Primary-Zone Combustor Liner Wall Temperatures with Calculated Predictions. Technical Report NASA-TM-X-2711, NASA, OH, USA, 1973. [37] J. Odgers and D. Kretschmer. Problems associated with wall temperature measurement and prediction in small combustors. In International Gas Turbine and Aeroengine Congress and Exposition, FL, USA, jun. ASME. [38] A. de Champlain P. Gosselin and D. Kretschmer. Prediction of Wall Heat Transfer for a Gas Turbine Combustor. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 213(169):180, 1999. [39] T. C. J. Hu P. Sampath and H. Ozem. Combustion technology challenges for small aviation gas turbines. In 23rd Congress of International Council of the Aeronautical Sciences, volume 5.8. International Council of the Aeronautical Sciences, sep 2002. [40] N. Pegemanyfar and M. Pfitzner. Development of a combustion chamber design methodology and automation of the design process. In Proceedings of the 25th International Congress of the Aeronautical Sciences, 2006. [41] K. Guedri K. Halouani R. M´echi, H. Farhat and R. Said. Extension of the zonal method to inhomogeneous non-grey semi-transparent medium. Energy, 35(1):1 – 15, 2010. [42] H. I. H. Saravanamuttoo, G. F. C. Rogers, and H. Cohen. Gas turbine theory. Pearson Education, 2001. [43] U. Schumann, editor. Impact of NOx Emissions from Aircraft Upon the Atmosphere at Flight Altitudes 8-15km (AERONOx). EL-DLR, August 1995.
Bibliography
63
[44] G. J. Sturgess and D. T. Shouse. A hybrid model for calculating lean blow-outs in practical combustors. In 32nd AIAA/ASME/SAE Joint Propulsion Conference, Joint Propulsion Conferences. AIAA, jul 1996. [45] P. J. Stuttaford. Preliminary Gas Turbine Combustor Design Using Network Approach. PhD thesis, Cranfield University, 1997. [46] P. J. Stuttaford and P. A. Rubini. Assessment of a radiative heat transfer model for gas turbine combustor preliminary design. Journal of propulsion and power, 14(1):66–73, 1998. [47] F. P. Incropera T. L. Bergman and A. S. Lavine. Fundamentals of heat and mass transfer. John Wiley & Sons, 2011. [48] S. Tietz and T. Behrendt. Development and application of a pre-design tool for aero-engine combustors. CEAS Aeronautical Journal, 2(1-4):111–123, 2011. [49] S. R. Turns. An introduction to combustion concepts and applications. McGraw-Hill, New York, 1996. [50] R. Viskanta and M. P. Meng¨ u¸c. Radiation heat transfer in combustion systems. Progress in Energy and Combustion Science, 13(2):97 – 160, 1987. [51] J. von Wolfersdorf. Effect of coolant side heat transfer on transpiration cooling. Heat and Mass Transfer, 41(4):327–337, 2005. [52] and E. E. Takara W. W. Yuen. The zonal method: A practical solution method for radiative transfer in nonisothermal inhomogeneous media. Annual review of heat transfer, 8(8), 1997. [53] A. Wadia. Advanced Combustor Liner Cooling Technology for Gas Turbines. Defence Science Journal, 38(4):363–380, 1988. [54] S. Mah K. Patel M. J. Dowhan Y. Panchenko, H. Moustapha and D. Hall. Preliminary multi-disciplinary optimization in turbomachinery design. In reduction of Military Vehicle Acquisition Time and Cost through Advanced Modelling and Virtual Simulation, Paris, France, apr 2002. RTO AVT.
