"Applied Science, Faculty of"@en . "Mechanical Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Georgallis, Mike"@en . "2009-11-27T23:21:26Z"@en . "2004"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "Rotary kilns have wide use in industry from the calcination of limestone to\r\ncement manufacturing to calcining of petroleum coke etc. These machines have survived\r\nand have been continuously improved (fuel efficiency, automation) for over a century.\r\nModelling has aided the design and operation of rotary kilns over the years. In the\r\npresent study, a three-dimensional steady-state model to predict the flow and heat transfer\r\nin a rotary lime kiln is presented. All important phenomena are included for the pre-heat\r\nand calcination zones. The model is based on a global solution of three submodels for the\r\nhot flow, the bed and the rotating wall/refractories. Information exchange between the\r\nmodels results in a fully coupled 3-D solution of a rotary lime kiln.\r\nThe hot flow model uses block-structured body-fitted coordinates with domain\r\nsegmentation. This results in a capability of simulating flows in complex three-dimensional\r\ngeometries with local refinement that capture details of the burner geometry\r\nand hood/kiln interactions. The CFD solver was developed by Nowak within the\r\nDepartment of Mechanical Engineering at the University of British Columbia. The\r\nmethod is based on the finite volume method. The model includes buoyancy effects,\r\nturbulence (using the standard \u00CE\u00BA-\u00CE\u00B5 model), evolution and combustion of gaseous\r\nspecies (through Magnussen's model), and radiation through the ray-tracing technique\r\n(Discrete Ordinate Method). Combustion of natural gas or oil may be modelled.\r\nThe bed model solves for a three-dimensional energy and species transport\r\nbalance in the non-Newtonian flow field. The field is modelled as two distinct Newtonian\r\nregions: an active layer adjacent to the hot flow and a plug flow region below. The bed\r\nmodel is based on the finite element method and was developed as part of this thesis. The\r\ncalcium carbonate reaction is modelled and assumed to be dominated by heat transfer.\r\nReaction proceeds as a function of temperature and a bulk particle size diameter. An\r\neffective thermal conductivity (modified by a mass diffusion factor) is used for the\r\ngranular bed heat transfer in the active layer. The bed flow field is supplied via a\r\nvorticity-stream function formulation and satisfies continuity equation. A 3-D\r\ntemperature field results along with distributions for CaCO3 and CaO. Carbon dioxide is\r\nreleased to the hot flow. A three-dimensional wall model, also developed as part of this thesis and based\r\non the finite element method, includes varying conductivity within the metal and\r\nrefractories. Rotation is included through a convective plug flow term in the energy\r\nequation (V = \u00CE\u00A9xR). Heat is lost to the ambient through free convection and radiation.\r\nThe overall model is validated using UBC's pilot kiln trials (5m laboratory kiln).\r\nThe solutions indicate that the present model may be used to identify problems of kiln\r\noperation or design. Prediction of a wealth of information is possible for various (primary\r\nair/fuel) input conditions, kiln geometries and burner designs. When the flow field,\r\ntemperature distributions and species concentrations are visualized and understood for a\r\ngiven kiln operating condition, optimization and diagnosis of operational problems can be\r\ngreatly enhanced through the use of the present computational tool. Analyzing complex\r\naerodynamic flows (secondary flows) of hood/kiln interactions through flame\r\nvisualization (and possible impingement) also has the potential of increasing refractory\r\nlife."@en . "https://circle.library.ubc.ca/rest/handle/2429/15899?expand=metadata"@en . "14272182 bytes"@en . "application/pdf"@en . "M A T H E M A T I C A L MODELLING OF L I M E KILNS By Mike Georgallis B.Eng., Concordia University, 1986 M.Eng., Concordia University, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES Department of Mechanical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 2004 \u00C2\u00A9 Mike Georgallis, 2004 Abstract Rotary kilns have wide use in industry from the calcination of limestone to cement manufacturing to calcining of petroleum coke etc. These machines have survived and have been continuously improved (fuel efficiency, automation) for over a century. Modelling has aided the design and operation of rotary kilns over the years. In the present study, a three-dimensional steady-state model to predict the flow and heat transfer in a rotary lime kiln is presented. All important phenomena are included for the pre-heat and calcination zones. The model is based on a global solution of three submodels for the hot flow, the bed and the rotating wall/refractories. Information exchange between the models results in a fully coupled 3-D solution of a rotary lime kiln. The hot flow model uses block-structured body-fitted coordinates with domain segmentation. This results in a capability of simulating flows in complex three-dimensional geometries with local refinement that capture details of the burner geometry and hood/kiln interactions. The CFD solver was developed by Nowak within the Department of Mechanical Engineering at the University of British Columbia. The method is based on the finite volume method. The model includes buoyancy effects, turbulence (using the standard K-e model), evolution and combustion of gaseous species (through Magnussen's model), and radiation through the ray-tracing technique (Discrete Ordinate Method). Combustion of natural gas or oil may be modelled. The bed model solves for a three-dimensional energy and species transport balance in the non-Newtonian flow field. The field is modelled as two distinct Newtonian regions: an active layer adjacent to the hot flow and a plug flow region below. The bed model is based on the finite element method and was developed as part of this thesis. The calcium carbonate reaction is modelled and assumed to be dominated by heat transfer. Reaction proceeds as a function of temperature and a bulk particle size diameter. An effective thermal conductivity (modified by a mass diffusion factor) is used for the granular bed heat transfer in the active layer. The bed flow field is supplied via a vorticity-stream function formulation and satisfies continuity equation. A 3-D temperature field results along with distributions for C a C 0 3 and CaO. Carbon dioxide is released to the hot flow. ii A three-dimensional wall model, also developed as part of this thesis and based on the finite element method, includes varying conductivity within the metal and refractories. Rotation is included through a convective plug flow term in the energy equation (V = QxR). Heat is lost to the ambient through free convection and radiation. The overall model is validated using UBC's pilot kiln trials (5m laboratory kiln). The solutions indicate that the present model may be used to identify problems of kiln operation or design. Prediction of a wealth of information is possible for various (primary air/fuel) input conditions, kiln geometries and burner designs. When the flow field, temperature distributions and species concentrations are visualized and understood for a given kiln operating condition, optimization and diagnosis of operational problems can be greatly enhanced through the use of the present computational tool. Analyzing complex aerodynamic flows (secondary flows) of hood/kiln interactions through flame visualization (and possible impingement) also has the potential of increasing refractory life. iii Table of Contents Abstract ii List of Figures viii List of Tables xii Nomenclature xiii Acknowledgements xviii 1.0 INTRODUCTION 1 1.1 Aerodynamics 6 1.2 Heat Transfer 11 1.3 Bed Transport Processes 15 1.4 Multi-Dimensional Modelling 17 1.5 Summary 21 2.0 SCOPE AND OBJECTIVES 23 3.0 PHYSICAL MODEL AND GOVERNING EQUATIONS 25 3.1 Hot Flow Model 26 3.2 Bed Model 32 3.2.1 Bed Species Transport Model 33 3.2.2 Bed Thermal Model 38 3.2.3 Bed Flow Field Model 40 3.2.4 Bed Reaction Rate Model 41 3.2.5 Summary of Governing Equations for Bed Model 44 3.3 Wall/Refractories Model 45 iv 4.0 COMPUTATIONAL PROCEDURE 47 4.1 Hot Flow Model 48 4.1.1 Inlet Boundaries 50 4.1.2 Outlet Boundary 51 4.1.3 Wall Boundaries 51 4.2 Bed Model 53 4.2.1 Bed Flow Field 56 4.2.1.1 Boundary Conditions 5 7 4.2.1.2 Discretization 58 4.2.1.3 Three Dimensional Velocity Field 62 4.2.2 Bed Species and Thermal Model 64 4.2.2.1 Boundary Conditions 66 4.2.2.2 Discretization 67 4.2.2.3 Contour Integrals for Mass and Thermal Fluxes 74 4.2.2.4 Granular Temperature Distribution 75 4.2.2.5 Calcination Temperature Distribution 79 4.3 Wall/Refractories Model 81 4.3.1 Boundary Conditions 83 4.3.1.1 Free Convection Heat Transfer Coefficient 85 4.3.1.2 Total Heat Loss to Surroundings 89 4.3.2 Discretization 90 4.4 Coupling and Information Exchange 93 5.0 GRID GENERATION 98 5.1 UBC Pilot Kiln 98 5.2 Grid Generation for the Hot Flow/Wall Coupling 100 5.3 Grid Generation For The Fully Coupled Run 103 v 6.0 CODE STRUCTURE 111 6.1 Overall Structure of Model and Convergence of Hot Flow/Wall Coupling 111 6.2 Overall Structure of Model and Convergence of Fully Coupled Run 115 7.0 INVESTIGATION OF AERODYNAMICS IN THE HOT END OF THE KILN 118 7.1 Aerodynamics, Turbulence and Combustion 118 7.2 Prediction Capability for the Isothermal Confined Jet 121 7.3 Error Involved in Using the Standard K \u00E2\u0080\u0094 \u00C2\u00A3 Model 128 7.4 Summary Remarks on Modelling of Lime Kilns 137 8.0 RESULTS AND DISCUSSION 141 8.1 Hot Flow/Wall Coupled Results 141 8.1.1 A Transient Solution of the Hot Flow 148 8.2 Fully Coupled Results 153 8.2.1 Comparison to Available Experimental Data 155 8.2.2 Hot Flow Results 159 8.2.3 Bed Results 161 8.2.3.1 Effect of Artificial Viscosity 168 8.2.4 Wall Results 171 8.2.5 Dominant Thermal Diffusion Directions and Possible Simplifications to the Models 173 8.3 Summary 174 9.0 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 175 BIBLIOGRAPHY 179 vi APPENDICES A 20-Node, 3-D Isoparametric Element Description 186 B Calculation of Surface Integrals 188 C Granular Temperature Distribution in Bed Active Layer 190 D Calcination Temperature Distribution for the Bed 192 E Schematic Diagram of the UBC Pilot Kiln 193 F Overall Heat Balance Model 194 Vll List of Figures 1.1 Kraft Process 1 1.2 The Rotary Kiln (Taken from Peray and Waddell [1]) 2 1.3 Fire Hood (Taken from Kramm, D.J. [2]) 4 1.4 Recirculation Patterns in the UBC Pilot Kiln (scale distorted) (A-Flow separation point; B-Flow reattachment point) 8 1.5 Heat Transfer Paths 12 3.1 Coupling of Models 25 3.2 Evaluation of Solids Volume (and mass) Fraction 34 3.3 Partially Decomposed Limestone Particle 35 3.4 Decomposition Pressure of Calcium Carbonate [59] 44 4.1 Surface Grid for the UBC pilot kiln 49 4.2 Boundary Conditions for Bed Active Layer Flow Field Model 57 4.3 Scalar Stream Function Distribution in Bed Cross-Section 62 4.4 Velocity Field in Bed 63 4.5 Boundary Conditions for Bed Model 66 4.6 Finite Element Matrix Assembly for K=KPLANE 73 4.7 Transfer ofC02 to Hot Flow Model 75 4.8 Experimentally Observed Granular Temperature Distribution at the Bed Surface 76 4.9 Imposed Bed Surface Granular Temperature Distribution 77 4.10 Depth Distribution of Granular Temperature in Active Layer 78 4.11 Cross-Sectional Granular Bodenstein Number Distribution 79 4.12 Cross-Sectional Calcination Temperature Distribution 80 4..13 Wall Zones for Conductivity Calculations 82 4.14 Wall Model Boundary Conditions 84 4.15 Illustration of Angle 6 for Inclined Plane of Length x 85 4.16 Comparison of Heat Transfer Correlations 88 viii 4.17 Grid Cross-Section with Hot Flow, Bed and Wall 94 4.18 Data Transfer Between Models 94 4.19 Transfer of Temperatures and Heat Fluxes between the Finite Element Description and the Control Volume Discretization 96 5.1 UBC Pilot Kiln 99 5.2 Firing Hood of UBC Pilot Kiln 99 5.3 Grid in a Cross-Section for Hot Flow/Wall Coupling 100 5.4 Grid in Central Vertical Plane for Hot Flow/Wall Coupling 101 5.5 Hot Flow Grid for Hot Flow/Wall Coupling 102 5.6 Wall Grid for Hot Flow/Wall Coupling 102 5.7 Grid Cross-Section with Hot Flow, Bed and Wall 103 5.8 Bed Geometry Description 106 5.9 Bed Grid for Fully Coupled Case 107 5.10 Cross-Section of Grid at Burner/Fire Hood Region 108 5.11 Hot Flow Grid for the Fully Coupled Case 109 6.1 Global Loop for Hot Flow/Wall Coupling 112 6.2 Global Convergence for Hot Flow /Wall Coupling 114 6.3 Global Loop for Fully Coupled Run 116 6.4 Global Convergence of Fully Coupled Run 117 7.1 Water Model Flow Diagrams of Moles et al. [6] (top half) Compared to CFD Calculations (bottom half) (Secondary to Primary Diameter Ratio=9.1) 120 7.2 Measurement for Spreading of Jet Lj and Duct Boundary Layer L B L 122 7.3 Spreading of Jet and Duct Boundary Layer (Exp. from Razinsky et al. [77]) 123 7.4 Center Line Velocity Decay (Exp. from Razinsky et al. [77]) 123 7.5 Center Line Velocity Decay with Pope's Correction [81 ] (Exp. from Yule et al. [80]) 126 7.6 Center Line Velocity Decay with Pope's Correction [81 ] (Exp. from Razinsky et al. [77]) 127 ix 7.7 Mean Velocity Distribution at Axial Position of 1.875 Inlet Diameters for the Conical Duct of Binder et al. [82] 129 7.8 Turbulent Stress Profile at Axial Position of 1.875 Inlet Diameters for the Conical Duct of Binder et al. [82] 130 7.9 Intermediate Grid of 51x20 Volumes Used ifor the Conical Duct Calculation of Binder et al. [82] 131 7.10 Mean velocity (left) and turbulent stress (right) solution with Ci=1.30 Experimental data from Binder et al. [82] 132 7.11 Mean velocity (left) and turbulent stress (right) solution with Ci=l .44 Experimental data from Binder et al. [82] 133 7.12 Mean velocity (left) and turbulent stress (right) solution with Ci=l .55 Experimental data from Binder et al. [82] 134 7.13 3-D Grid Used for Generic Industrial Kiln 136 7.14 Temperature (K) in Central Vertical Plane of Generic Industrial Kiln For Values of Ci of 1.30, 1.44 and 1.55 (Axes are in meters) 137 7.15 Temperature (K) in Central Vertical Plane of Generic Industrial Kiln With and Without Buoyancy (Axes are in meters) 138 7.16 Temperature (K) in Cross Sectional Plane at Axial Position of 6.6 m for Generic Industrial Kiln (Axes are in meters) 139 8.1 Fuel (CH4) Mass Fraction (length scale distorted) 143 8.2 Axial Velocity Field (m/s) (length scale distorted) 144 8.3 Wall Hot Spots in Pilot Kiln (Inner Wall Temperatures in K) 145 8.4 Axial Temperature Profiles in Pilot Kiln (Exp. Run#6 Alyaser [61]) 145 8.5 Hot Flow Temperature Contours in Flame Zone (K) 146 8.6 Inner Wall Temperatures of Pipes, Burner and Fire Hood (K) 147 8.7 Fuel Mass Fraction (left side) and Temperature (K) (right side) on Central Vertical Plane (Vertical axis is amplified 3X) 149 8.8 Fuel Mass Fraction (left side) and Temperature (K) (right side) at 1.78 Kiln Diameters Downstream of Burner Exit 151 8.9 Axial Temperature Profiles in Pilot Kiln (Exp. Run T21 Barr [62]) 156 8.10 Bed Heat Transfer Rates (Exp. Run T21 Barr [62]) 158 8.11 Heat Loss Through Shell (Exp. Run T21 Barr [62]) 158 8.12 Hot Flow Temperature Contours in Fully Coupled Results (K) 160 x 8.13 Inner Wall Temperatures (K) in the Fire Hood Region for Fully Coupled Results 161 8.14 Bed Surface Temperature Distribution (K) 162 8.15 Cross-Section Bed Temperature Profiles (K) at Axial Positions of (a) 0.5 m (b) 2.5 m and (c) 4.5 m 163 8.16 Solids Volume Fraction of CaC03 in Bed 164 8.17 Cross-Section Solids Volume Fraction of CaCC\u00C2\u00BB3 at Axial Positions of (a) 0.0 m(b) 0.5 m and (c) 1.3 m 165 8.18 Surface Contour for the Reaction Rate at a Value of 0.0005 per sec 166 8.19 Reaction Rate Contours (per second) at Axial Positions of (a) 0.5 m and (b) 1.5 m 167 8.20 Peclet Numbers in Active Layer 169 8.21 Peclet Numbers in Plug Flow Region 170 8.22 Inner Wall Temperature (K) 171 8.23 Cross-Section Wall Temperature Distribution (K) at Axial Position of 2.5 m 172 8.24 Inner Wall Temperature as a Function of Circumferential ( 6 ) Direction 173 A. 1 20-Node, 3-D Element Used for the Bed and Wall Models 186 B. l Mapping of 3-D surface 188 C l Depth Distribution for Granular Temperature, S 191 E. 1 Simplified Schematic of the UBC Pilot Kiln (not to scale) 193 E.2 Simplified Schematic of the Fire Hood/Burner Arrangement in the UBC Pilot Kiln (not to scale) 193 x i L i s t o f T a b l e s 3.1 Components of Equation 3.1 28 3.2 Constants in the turbulence model [46] 29 4.1 Constants Used in Bed Model (for case T21 of Barr's experiments [62]) 65 4.2 Reference Variables for the Non-Dimensional Bed Equations 67 4.3 Fluxes Evaluated in the Post-Processing Phase (Dimensional Form) 74 4.4 Constants in Wall Model 82 4.5 Wall Conductivity Curve Fit Data for Alyaser's Experiments [61] 83 4.6 Wall Conductivity Curve Fit Data for Barr's Experiments [62] 83 4.7 Correlation Constants for Equation 4.46 (-90 < 6 < 90) 86 4.8 Reference Variables for the Non-Dimensional Wall Equation 90 5.1 Discretization in Each Segment for the Hot Flow/Wall Coupling Grid 101 5.2 Input Data for Bed Geometry 105 5.3 Discretization in Each Segment of the Fully Coupled Hot Flow Grid 109 6.1 Computing Times for Hot Flow/Wall Coupling (Based on a 3 GHz Pentium IV Processor) 112 6.2 Computing Times for Fully Coupled Run (Based on a 3 GHz Pentium IV Processor) 116 8.1 Pilot Kiln Trial Run #6 from Alyaser [61] 142 8.2 Model Heat Transfer Rate Imbalances for Hot FlowAVall Coupling 142 8.3 Comparison of Heat Transfer Rates to Experimental Data of Alyaser [61] 142 8.4 Pilot Kiln Trial Run T21 from Barr [62] 154 8.5 Model Heat Rate Imbalances for Fully Coupled Results 154 8.6 Comparison of Heat Rates to Overall Heat Balance Model of Appendix G 155 xii Nomenclature OR Rosseland mean absorption coefficient (1/m) absorption coefficient (1/m) A area (m2) c, constant in K - e turbulence model C2 constant in K - e turbulence model Q constant in Pope's modification to K-e model specific heat (J/kgK) constant in K-e turbulence model d diameter of particles (m) artificial diffusion coefficient (m Is) D8 granular diffusion coefficient (m2/s) Dm,D\u00E2\u0080\u009E inner kiln diameter (m) EP coefficient of restitution (-) E energy (J) f solids mass fraction of CaCCb, fractional fill for bed (-) g solids mass fraction of CaO (-) g gravitational constant vector (m/s ) ga(v) radial distribution function (-) G generation term (equation 3.1) (kg/ms3) h enthalpy (J/kg) effective heat transfer coefficient (W/m K) h spectral radiative intensity (W/m3rad) black body radiative intensity (W/m rad) K conductivity (W/mK) xiii Ka artificial conductivity (W/mK) Kg granular thermal conductivity (W/mK) Kw wall conductivity (W/mK) K' bulk bed conductivity (W/mK) /, turbulent length scale (m) L kiln length (m) m2 fraction of original mass of CaC03 converted to CaO (-) w3 fraction of original mass of CaCC\u00C2\u00BB3 converted to C O 2 (-) m entrain entrainment rate (kg/s) m, mass fraction of species i (-) m p primary jet mass flow rate (kg/s) ms secondary mass flow rate (kg/s) MW; molecular weight of species i (kg/kgmole) n unit normal vector (-) P pressure (N/m2) q excess volumetric flow rate of primary discharge (m3/s) qfc heat loss to surroundings by free convection (W/m2) qrad heat loss to surroundings by radiation (W/m2) qr radiative heat flux vector (W/m ) Q total volumetric flow rate (m3/s) QCOI rate of heat lost with C02 (W) QCOND r a t e \u00C2\u00B0f heat transfer by granular and bulk diffusion (W) QP rate of heat to a particle (W) QR radiative source term for equation 3.1 (W/m3) QR rate of heat absorbed by reaction (W) rt radius at particle reaction interface (m) xiv rp primary jet radius (m) ki ln jet radius, limestone particle radius (m) R reaction rate (1/s) RIN inner ki ln radius (m) R' excess discharge ratio (q/Q)(-) R universal gas constant (J/(kgmole)(K)) R radius vector (m) s ray direction (m) non-dimensional rate of strain tensor (-) SS source terms (equation 3.1) T temperature (K) Tc calcination temperature (K) Tf f i lm temperature (K) Tsh outer wall shell temperature (K) T oo ambient temperature (K) U velocity in coordinate direction x (m/s) UP average primary velocity (m/s) Us average secondary velocity (m/s) u,um mean inlet velocity (m/s) inlet center line velocity (m/s) V velocity in coordinate direction y (m/s) V velocity vector (m/s) w velocity in coordinate direction z (m/s) XV Greek: aBED angle of repose (deg.) diffusion coefficients (equation 3.1) S' boundary layer displacement thickness (m) \u00C2\u00A3 turbulent energy dissipation rate (m2/s3) es emissivity (-) 6 angular position (deg, rad) K turbulent kinetic energy (m2/s2) X wavelength (m) (1 molecular dynamic viscosity (kg/ms) jue effective viscosity (//,+//) (kg/ms) fl, turbulent viscosity (kg/ms) v solids fraction (-) \u00C2\u00A3 unknowns (equation 3.1) E granular temperature (m2/s2) p density (kg/m3) Pcaco, effective density of CaC03 (kg/m3) pCa0 effective density of CaO (kg/m3) pw wall density (kg/m3) a Stefan-Boltzmann constant (W/m2K4) 3, into NaOH. In the calcium loop of the recovery cycle, lime kilns are used to convert lime mud back to quick lime for reuse in the causticizing process. The rotary lime kiln considered in this work is effectively a direct-contact counter-flow heat exchanger. The endothermic reaction that takes place in the lime uses energy from the hot burnt gases that flow above the mud. The rotary kiln as shown in Figure 1.2, is essentially a long cylinder typically 3 m in diameter and 90 m in length. It rotates about its axis at roughly 1 rpm and is sloped at 2 degrees. Lime mud is fed from the elevated cold end and moves down the kiln (into the far field of fig. 1.2) due to rotation and gravity. The retention time for the mud in the kiln is two to three hours. The hot end is typically maintained at about 1200 C by burning fuel (natural gas or oil). Figure 1.2: The Rotary Kiln (Taken from Peray and Waddell [1]) 3 There are three basic zones in the operation of the kiln, in the following order starting from the cold end; drying, preheat and calcination. In the initial drying zone, water is evaporated from the mud. Heat transfer rates are enhanced in this zone by the use of steel chains attached to the interior of the kiln shell. Lime is then heated until it reaches a temperature of about 8 1 0 C . That allows for calcination to take place. In the burning zone the metal shell is lined with refractory brick. This protects the metal shell from the hot gases. Refractory failure is considered the most critical upset in kiln operation primarily due to the expense of replacing and re-installing the lining [ 1 ] . The life of the refractory lining in the hot section can vary from 1 or 2 months to more than a year. Although this wide range can be explained by mechanical factors such as refractory type and proper installation, kiln operation is also considered to be an important contributor to life expectancy. Flame impingement on the wall for example can cause early deterioration of the lining. The fire hood in the hot end is shown in Figure 1 . 3 (no integral tube coolers are used in the kiln shown here). The refractories are visible on the interior of the rotating cylinder. The lime product (in the form of spherical pellets) drops from the cylinder to the chute below at the bottom of the fire hood and moves on to the cooler. Primary air and fuel are supplied through the centrally placed burner. It is important to note here that many problems associated with the kiln operation, such as poor product quality, low fuel efficiency and short refractory life can be related to the design and operation of the burner. Secondary air enters through numerous openings. Typically the inspection door in the back of the fire hood and the access door on the side are at least partially open to allow for adequate secondary air. Secondary air may enter through the bearings at the 4 interface between the hood and the rotating kiln. In the absence of integral tube coolers, air also tends to come up through the chute and into the kiln. Figure 1.3: Fire Hood (Taken from Kramm, D.J. [2]) Physical phenomena present in the kiln include; three dimensional turbulent gas flow with combustion, three dimensional non-Newtonian flow in the lime mud (bed motion), radiative, convective and conductive heat transfer within and between the various components of the kiln, namely the hot gas, the bed, the refractories, the metal shell and the surroundings. In addition there is heat transfer to the bed from the rotating wall. This is commonly called heat regeneration. The modelling in the present work includes all the important phenomena in the preheat and calcination zones. The remainder of this chapter surveys the literature on various operational aspects of the lime kiln, namely, the aerodynamics, heat transfer and bed transport processes. Multi-dimensional modelling efforts to date are then summarized and motivation for current work is given. Chapter 2 outlines the scope and objectives of this work. Chapters 3 and 4 describe the mathematical models and computational procedures used. Chapter 5 is dedicated to the topic of grid generation followed by the overall code structure in Chapter 6. Chapter 7 investigates the accuracy of modelling assumptions made in predicting the flow patterns in the hot end. Chapter 8 presents the coupled results and final concluding remarks are given in Chapter 9. 6 1.1 Aerodynamics Understanding and predicting aerodynamic phenomena in rotary kilns can have many pay-offs for the industry through gains in knowledge for design and operation of these furnaces. Efforts to date in the prediction of the turbulent diffusion flame in the rotary ki ln have ranged from an estimation of the flame length alone [3] to detailed three-dimensional flow studies [4]. A n empirical flame length model for example was developed by Ruhland [3] for flames in the cement rotary ki ln. This empirical model was based on experiments which used two liquid streams (primary and secondary) whose interaction represented the mixing in the burner region. Assuming that \"burnt is mixed\", the well mixed region also represents the flame itself, so that the effective length could be visualized experimentally. The primary jet was represented by a mixture of caustic soda, water and thymolphthalein for the visual indicator. The secondary stream included a dilute solution of hydrochloric acid. Although a liquid representation was used, equivalence between the observed \"flame length\" and the flame length in a real ki ln was confirmed by measurements of the flame in a full scale rotary ki ln. The change in volume of substances during combustion was accounted for by using a volumetric flow ratio between entrance conditions and end-of-flame conditions (using adiabatic flame temperature information). Although useful, the correlation predicts the length of flames in kilns for a non-swirling jet of fuel from a smooth tubular nozzle. Additional aerodynamic complications such as hood/kiln interactions and detailed distortions of flames from buoyancy would require further empirical parametric investigation. 7 Both the geometry and Reynolds number are claimed to be similar between the model liquid flame and the rotary kiln flame. As described below, this makes Ruhland's correlation suspect. The primary core flow (fuel plus air) entrains secondary air as it would in a free jet. Ricou and Spalding [5] measured this 'entrainment appetite' of a confined jet by surrounding the turbulent jet with a porous-walled cylindrical chamber. Air was injected through the wall until no pressure gradients were detected. The measured injected flow rate was then assumed to be equal to that which would be entrained by a free jet. The primary jet momentum (consisting of fuel and primary air) is typically at least one order of magnitude higher than the surrounding secondary air momentum entering the kiln hot end. The resulting primary to secondary momentum ratio may or may not result in a recirculation zone further down the kiln. The size, strength or existence of this recirculation eddy has important implications for rotary kiln flame stability, flame length and emissions. Figure 1.4 shows a hot flow calculation for the UBC pilot kiln (including buoyancy). This was calculated as part of the present research and is given here for illustration. Further details for this flow field (code used, boundary conditions, geometry, etc.) will be given later in this thesis. The axial velocity field reveals recirculation in the burner, the fire hood and a recirculation eddy attached to the bottom wall further down the kiln. As the core jet expands into the confined furnace it entrains secondary air. The momentum of the core flow (relative to the secondary stream) is high enough to reduce the momentum of the secondary stream (and cause adverse pressure gradients) along the wall which result in reverse flow at the wall. The entrainment capacity of the core flow is 8 said to exceed the secondary supply. The degree of recirculation therefore depends on the entrainment capacity of the core flow and more specifically on the momentum ratio between primary and secondary streams. Axial Velocity Field (m/s) (scale distorted) 0.2 If c o T> 0 a> b -0.4 -0.4 -0.2 0 0.2 0.4 0.S 0.8 1 1.2 1.4 1.6 Axial Direction (m) Figure 1.4: Recirculation Patterns in the UBC Pilot Kiln (scale distorted) (A-Flow separation point; B-Flow reattachment point) Reynolds number typically plays a secondary role. Flow exiting the burner is always turbulent. The Reynolds number of the turbulent diffusion flame (i.e. the core flow) based on the mean primary flow velocity and the burner diameter, is typically high enough (>105) in rotary kilns to be Reynolds number independent [6]. Experiments such as Ruhland's could have been carried out at lower Reynolds number without a great loss of similarity provided they were kept above 104. The confined jet problem described above has been identified and studied since the early 1950's. Several attempts have been made to obtain suitable similarity parameters for modelling this problem. Core Air Entrainment 1 \u00E2\u0080\u00A2' \u00E2\u0080\u00A2 i \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2' i \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 i \u00E2\u0080\u00A2 \u00E2\u0080\u00A2' \u00E2\u0080\u00A2 i ' > \u00E2\u0080\u00A2' i \u00E2\u0080\u00A2''' i ' ' \u00E2\u0080\u00A2' i ' \u00E2\u0080\u00A2 \u00E2\u0080\u00A2' i '' \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 i \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 Thring and Newby [7] assumed that until the core flow reaches the wall (reattachment point of the recirculation eddy in fig. 1.4) it behaves as a free jet. Using Hinze's formula for the entrainment rate m entrain [8], x ntemrmn = m p ( 0 . 2 1) (1.1) rp and the universal jet spreading angle, x =A.5rs (1.2) they proposed a similarity parameter based on the separation to reattachment distance ratio. Here rp and mp are the radius and mass flow rate of the primary jet. rs is the kiln radius and x the axial distance from the nozzle exit to the point where the core flow spreads to the wall. The mass flow ratio at which the secondary flow exactly satisfies the entrainment requirement (ms = wen / r am) as the core flow spreads to the wall, could be evaluated by equating x to x. In this unique situation (77z = 1 below) no recirculation would be observed. A similarity parameter Th was proposed as follows, Th f . . \ m,, + ms V m P J rp _ Separation..distance ( 1 3 ) r Re attachment, .dis tan ce The reattachment distance is the axial distance between burner exit and point B in fig. 1.4. The separation distance is the axial distance between burner exit and point A in fig. 1.4. The entrainment requirements are satisfied for Th>\. Recirculation is present for Th < 1. This parameter generally does not apply to typical rotary kilns for which rp jrs > 0.05. Under such conditions there are excessive flow distortions in the nozzle region [6]. 10 By simplifying the equations of motion for a constant density fluid and introducing experimental observations, Craya and Curtet [9, 10] defined a dimensionless parameter m representing the thrust term (momentum flux and pressure force). Parameter m is a function of the excess discharge ratio R , which is the ratio of the excess volumetric flow rate of the primary discharge q to the total volumetric flow rate Q. Their analysis relates m to R as follows: m = R- 1.5R +0.579 - (1.4a) ( V O 2 where R =q/Q, (1.4b) q=nr2p(up-us) (1.4c) and Q = 7t(r, -S')2us +q (1.4d) Here up and us are the average velocities for the primary and secondary flows. The boundary layer displacement thickness 8* is usually negligible. The Craya-Curtet parameter m is also identified in the literature as a Curtet number [11]. Ct = -L (1.5) Recirculation occurs for a Craya-Curtet parameter of approximately m>1.5 (Ct<0.8). Moles et al. [6] have surveyed 53 rotary kilns and found that they operate in the range of medium to low recirculation intensity (0.4 C02 + 2H20 (-50 010 kJ/kg CH4) (3.4) The number that follows in brackets is the enthalpy of reaction. The source term Sm in the species equations is evaluated using Magnussen's Eddy-Dissipation Model [47]. The combustion rates are assumed to be limited by the mixing of turbulent eddies and are predicted using a turbulent time scale tc/e in conjunction with reaction 3.4. That is, Smi=AmiP- (3.5) K 30 The coefficient A is set at a constant value of 4.0 in the entire field for the turbulent diffusion flame calculation [47]. The rate may be limited by the available oxygen. For the incompressible flow under consideration the density of the gas mixture is a function of temperature T and mass fraction. (3.6) where MWt is the molecular weight of each species and R is the universal gas constant. For this calculation, pressure P is constant at 101.3 kPa. The energy equation solves for enthalpy h. JANAF thermochemical tables [51] are used to relate enthalpy to temperature. The source term Sh in the energy equation includes combustion and radiation heat transfer rates. The combustion portion of the energy source term is obtained by taking the product of the consumption rate of the fuel with the enthalpy of reaction for equation 3.4. The pilot kiln modelled here has a relatively 'clean' hot flow environment due primarily to low dust levels rising from the bed. The burning of natural gas also results in a negligible amount of carbon (soot) particles being produced. This 'clean' environment permits radiation calculations to be carried out for a non-luminous hydrocarbon flame where only absorption and emission may be considered from the chief gaseous radiative constituents (CO2 and H2O). 31 The full equation of radiative transfer [52] is solved (without scattering) using the Discrete Ordinate Method [48-50]. ^T = ax^{s)-axh{s) (3.7) ds The unknown variable i\ is the spectral radiative intensity. The LHS of equation (3.7) represents the change of radiative intensity along a ray (direction s) for a given wavelength A. The prime indicates a single direction. One such equation is solved for each spectral interval of the gaseous absorption bands. The first term on the RHS is the gain in intensity by emission. Here, i'u denotes the black body intensity \u00C2\u00A3 iihdX = \u00E2\u0080\u0094 . The second term on the RHS is the loss in intensity by absorption. The absorption coefficient ax is generally a function of wavelength, temperature, pressure and composition. A total pressure of 1 atm is again assumed. The gray-band model [53] with 5 intervals for each species is used to integrate the spectral intensity. When calculating the radiative portion of the energy equation source term, the divergence of the radiative flux crossing a given area as a result of intensities from all directions is considered (- V \u00E2\u0080\u00A2 qr) where, v \u00E2\u0080\u00A2 = L^[ 4 ^ \u00E2\u0080\u00A2 f <3-8) The spectral intensity is first integrated over all solid angles co and then over the 5 spectral bands for CO2 and H2O. 32 3.2 Bed Model The bed model developed in the present study solves for a three-dimensional energy and species transport balance in the non-Newtonian flow field. The non-Newtonian field is approximated and modelled as two distinct Newtonian regions: an active layer adjacent to the hot flow and a plug flow region below [4, 36]. The calcium carbonate reaction is modelled and is assumed to be dominated by heat transfer [54]. The reaction proceeds as a function of temperature and a bulk particle size diameter. Within the reaction rate model, diffusion effects of CO2 within the bed are accounted for by imposing a partial pressure distribution as a function of bed depth. An effective thermal conductivity (modified by a mass diffusion factor [55, 56]) is used for the granular bed heat transfer in the active layer. The bed flow field is supplied via a vorticity/stream function formulation and satisfies the volumetric continuity equation. Assumptions made in deriving the governing equations for the bed are summarized as follows: 1) The bed is under steady-state conditions. 2) The non-Newtonian flow field is approximated as two distinct Newtonian regions. 3) The diameter of a limestone particle is assumed to remain constant through the calcination process. As a consequence the volumetric flow rate of the solids is assumed constant. 4) A constant vorticity value (equal to 2 times the rotation rate of the wall) is assumed in deriving a simplified velocity field. The axial component of 33 velocity is assumed to be constant in the active layer. Forward motion is assumed to be zero in plug flow region. 5) Constant effective densities for each of the solid species are assumed. These account for voids and fines in the granular matrix. 6) CO2 gas is assumed to be instantaneously transported to the hot flow but not in the development of the reaction rate model. 7) A constant specific heat is assumed for the solids. 8) Granular motion is assumed to follow a diffusive type behaviour (Fick's law for mass diffusion and Fourier's law for heat conduction). 9) Reaction in the bed is assumed to be dominated by conductive heat transfer. Diffusion of C02 within the bed (for the reaction rate model only) is assumed to follow an imposed partial pressure distribution. 10) Radiation heat transfer within the bed is ignored. The model developed here results in a 3-D temperature field along with distributions for CaC03 and CaO. Carbon dioxide is released to the hot flow. 3.2.1 Bed Species Transport Model Due to heating from the hot flow (and wall) the bed undergoes an endothermic calcination reaction according to, CaC03 (s)=> CaO (s) + C02 (g) (+1800 kJ/kg CaC03) (3.9) 34 The solids (CaCCb and CaO) are present in the bed. CO2 gas generally escapes to the hot flow region. Assuming an initial mass M0 occupying a volume V enters the bed inlet, the change in bed composition due to the above reaction is illustrated in Figure 3.2. \u00C2\u00AE CaO (6) CaC03 (s) Figure 3.2: Evaluation of Solids Volume (and mass) Fraction Three assumptions are made for the species occupying volume V (rectangle in the figure). First, an effective density for each species is considered to account for voids and fines in the granular matrix. Volume V, in figure 3.2 is occupied completely by CaC03 provided that the effective density pCaC0} is used. The value of pCaCOi used here is 1680 kg/m3. Its true density is 2800 kg/m3 at 25 C. Similarly the effective density in volume V2 is pCa0 (940 kg/m3). Second, CO2 gas is assumed to be instantaneously transported to the hot flow. No transport equation for the gas phase in the bed is considered. (It is important to note that for the reaction rate model [section 3.2.4] however a CO2 distribution as a function of bed depth is imposed). The total volume V is occupied by the solids only. That is, V = V,+V2 (3.10) 35 The third and final assumption is that the volume occupied by the solids does itself not change. This is a reasonable assumption considering calcination takes place as a topochemical or shrinking core reaction [54]. Figure 3.3 shows a partially decomposed limestone particle. Heat Reaction Interlace \u00E2\u0080\u00A2wer Figure 3.3: Partially Decomposed Limestone Particle As the reaction proceeds inward, a layer of porous CaO is maintained with outer diameter d. Although in practice parts of the particle may break away, the assumption made here is that the volumetric flow rate of the solids is constant, even though the density is changing to satisfy mass conservation and the chemical reaction. In developing the transport equations for the solids, variable

)M0m3. The effective density p of the solids may be derived from equation 3.10. The result is, P=PcaCOi+Q-)PcaO (3 \u00E2\u0080\u00A21 2) The solids volume fraction of CaC03 is in fact , and for CaO, (1 - (/>). The solids mass fraction for CaC03 and CaO are given below as / and g respectively. / = (3-13) (p + (\- 0 ) is the rate of change of mass for CaC03 then from mass conservation of this species, \u00C2\u00AB*P = - r + *.V.[D,V#] (3.15) Equation 3.15 may now be used in an Eulerian frame of reference. With the use of Fick's law of mass diffusion, diffusive effects are added on the RHS of the equation for the granular flow behaviour in the bed. Dg is the granular diffusion coefficient. This will be described in greater detail below. On a per volume basis equation 3.15 becomes, D(j) - V - W = - J I (3.16) Dt where R = r/(V/?CoCOj) is the reaction rate (in units of 1/sec). Details on the form of R will be given in section 3.2.4. Since steady-state conditions are assumed, equation 3.16 may be recast as, 3 7 VV)M0m2CpT (3.22) where T is the temperature. A constant specific heat Cp is used for the solids. The first two terms on the RHS of equation 3.21 are both negative (sinks of energy in V ) and are modelled according to reaction rate R . QR=-MohCaC0R (3.23) Qco2 =-M0m,RCpcoT (3.24) 39 Enthalpy of reaction hCaC0^ for calcination is used in equation 3.23. Specific heat of C O 2 Cpco in equation 3.24 is a function of temperature. This ensures proper communication when enthalpies are supplied to the hot flow model. The final term on the RHS of equation 3.21 is for conductive heat transfer. Noting that the final form of the energy equation developed here will be used in an Eulerian frame of reference, Fourier's law of heat conduction gives QC0ND in the following form. QCOND=W[KVT] (3.25) Conductivity K is the sum of a bulk conductivity value K in the bed and a granular conductivity K for the active layer. K = K + Kg (3.26) Kg has also been developed in Ffsiau and Hunt's kinetic theory analysis [55]. Kg=~^ (3.27) Substituting equations 3.22 through 3.25 into 3.21 yields (on a per volume basis), r v ^ r l - v . ^ r l - ^ ^ + ^ ^ j ) (3.28) Steady-state conditions are assumed as in the Eulerian species transport equation. Equation 3.12 was used for the density in the derivation of the convective term. 40 3.2.3 Bed Flow Field Model Since the volume occupied by the solids (V in figure 3.2) does not change, a divergence free velocity field is sought for. VV = 0 (3.29) To obtain a solution for the 3-components of velocity (u, v and w) in the 3-directions (x, y and z) the momentum equations would ideally be used along with equation 3.29. Any solution would be an estimate of the non-Newtonian flow field. The approach taken here is to assume 1 velocity component and the vorticity so that the volume continuity equation (3.29) is satisfied. This is described in more detail below. The implicit assumption made is that an approximation to the velocity field for the convective effects in the species and thermal transport equations (3.17 and 3.28) would be a good approximation to the actual complex flow distribution in the bed. The velocity field in the axial direction (w-component) is imposed. A constant value is supplied in the active layer corresponding to the solids residence time. In the plug flow region no forward motion is allowed since particle motion follows the wall. With a constant w-component imposed, equation 3.29 becomes, fHr-\u00C2\u00B0 <3'0> ox dy Rather than imposing another velocity component, equation 3.30 is transformed to a vorticity/stream function formulation. Defining a scalar stream function *F such that, dl1 u = -\u00E2\u0080\u0094 (3.31) dy 41 3 ^ and, v = \u00E2\u0080\u0094 (3.32) ox the volumetric continuity equation (3.30) is satisfied identically. With no velocity gradients in the axial direction the only component of vorticity is coz. dv du coz= \u00E2\u0080\u0094 -\u00E2\u0080\u0094 (3.33) ox ay Substituting equations 3.31 and 3.32 into 3.33, Vorticity is imposed as a constant in the entire field (equal to 2 times the rotation rate of the wall). Due to the highly diffusive nature of the active layer vorticity would normally vary. The boundary conditions for *F are also taken as a constant at the periphery of a bed cross-section. Further details will be given in Chapter 4. 3.2.4 Bed Reaction Rate Model The topochemical or shrinking core reaction of a limestone particle has been described in section 3.2.1. It is assumed that the overall heat of reaction is controlled by heat conduction. A Fourier conduction model based on a shrinking core reaction is developed here for the consumption rate of CaC03. The reaction interface as shown in figure 3.3 is between CaO and CaC03. The decomposed CaCC\u00C2\u00BB3 shell has a low effective thermal conductivity and acts as a thermal barrier [54]. For calcination to occur, heat 42 from the surroundings must be transferred to the reaction interface through this barrier where it would then be absorbed by the endothermic reaction. The rate of heat transfer to an individual spherical particle Q (>0) is, Qp=KCa04m-dT_ dr (3.35) where KCaQ is the conductivity of the porous calcined layer of figure 3.3. After integration from radius rt to rs (from reaction interface to surface), 4xKCa0(Ts -Tt) r. - r rsri (3.36) where 7]. is the core temperature and Ts the particle surface temperature. Volume fraction

0). This heat must be absorbed by the reaction as given by equation 3.23. That is, QpN-MohCaCOR = 0 Substituting equations 3.36 through 3.38 into 3.39, (3.39) 43 R = \ \"CaCO, PcaCO, J\ Ks Moo \T,-Tt (3.40) The first term in the brackets is a constant. The second term indicates that the reaction rate is proportional to the temperature differential above the core (or calcination) temperature and inversely proportional to r2 (or surface area) of the particle. The final term in the brackets represents the reduced resistance a particle would have to incoming heat at high values of

-V [DgV]=-R V \u00E2\u0080\u00A2 v\pCPT]- V \u00E2\u0080\u00A2 [KVT] = -pCaCO) R(hCaCOi + m3CpcQ2 T) respectively, where the density of the solids is given by, P = PcaCO, +Q--)PcaO and the reaction rate R is given by, (3.17) (3.28) (3.12) R = 3vK, CaO nCaCOy PcaCO (3.41) 3.3 Wall/Refractories Model The wall model includes varying conductivity Kw within the metal and refractories using manufacturer-supplied curve fits. Rotation is included through a convective plug flow term in the energy equation. v\pwVwCpJ-KwVT]=0 (3.42) where, Vw=ClxR (3.43) 46 Equation 3.42 is solved in 3-D for temperature with the velocity vector for the wall defined by equation 3.43. A constant density pw and specific heat C are imposed. Heat is lost to the ambient through free convection and radiation. Free convection is modelled as a series of individual surface inclined planes according to Fujii and Imura [60]. This will be described in Chapter 4. Nusselt number correlations are calibrated to full cylinder correlations. 4.0 COMPUTATIONAL PROCEDURE 47 Governing equations for the three models used in this study (hot flow, bed and wall) have been described in Chapter 3. Details of the computational procedures used in solving these equations (computational methods, boundary conditions, discretization, etc.) are given in this chapter. The hot flow model uses a CFD solver based on the control volume method (CVM) [45]. Discretization and boundary conditions used are described in section 4.1. Unlike the hot flow model, the bed and wall models are developed here as part of the current research. The models are based on the finite element method (FEM). The development of the set of algebraic equations and methods used to solve them are described in sections 4.2 and 4.3. Coupling and information exchange between models is described in section 4.4. Two coupled runs are introduced that model experimental data from the UBC pilot kiln trials. Alyaser's experiments [61] include the hot flow with a non-rotating wall. Run #6 is modelled here with a burner load of 77 kW. Run T21 from Barr's experiments [62] is also modelled with a burner load of 88 kW. This includes the hot flow, bed (with limestone feed) and rotating wall. Details of boundary conditions imposed, etc. for the aerodynamic studies presented in chapter 7 will be covered in chapter 7. 48 4.1 Hot Flow Model The hot flow model [45] used to solve equation 3.1 (with components given in Table 3.1) including combustion and radiation, transforms the physical space into a computational domain for discretization. For better numerical efficiency, segments of the flow regions (inlet pipes, burner, fire hood, kiln) are divided into blocks and discretization is performed independently on each block. This domain segmentation strategy results in a capability of simulating flows in complex three-dimensional geometries with local refinement that capture details of the burner geometry and hood/kiln interactions. Conservation equations are discretized by the 1st order hybrid scheme [63]. This scheme is conservative and positive (off-diagonal coefficients of the discrete equations are non-negative). This prevents the oscillation behaviour of converged solutions. Convergence is easier to achieve than for higher order schemes, however significant false diffusion occurs especially when the flow is inclined to the grid. Second order schemes (higher accuracy and producing finer details of the solution) are available in this model but generally suffer from deteriorating convergence when applied to complex flow fields and geometries. The hot flow code also allows the division of the flow field boundary into patches. Figure 4.1 shows the computational domain boundary for the UBC pilot kiln modelling. Details of the development of this grid will be described in chapter 5. The figure shows the burner, fire hood and individual jets for fuel and air inlets. Eight primary and 8 49 secondary air jets are modelled. Boundary conditions imposed on various patches are described below. Laminar dynamic viscosity used for these calculations was 2.0(10s) kg/ms. This is the approximate value for air at about 300 K, but because of the high effective viscosity generated by turbulence, and calculated in the model, this laminar value is of little or no consequence here. Prandtl and Schmidt numbers were assumed constant in the flow field and both set to 1.0. Figure 4.1: Surface Grid for the UBC pilot kiln 50 4.1.1 Inlet Boundaries The total mass flow rate through each individual jet inlet is imposed as a uniform flow. Free stream jet inlet values for the turbulent kinetic energy K\u00E2\u0080\u009E and its dissipation rate f\u00E2\u0080\u009E are imposed by, Jr\u00E2\u0080\u009E =1.5(72 fi]1 (4..) TI and E_=c{* \u00E2\u0080\u0094 (4.2) where \u00C2\u00A3/ is the mean inlet velocity, = is the turbulence intensity (set to 5%) and /, the turbulent length scale (set to 10% of a characteristic inlet geometric length such as the diameter of the inlet pipe). Since inlet turbulence intensities are very small relative to the amount generated in the field, these assumptions should not affect the results significantly. The fuel input (natural gas) is assumed to be pure methane. Air inlets are assumed to have volume fractions of 0.21 and 0.79 for 02 and N2 respectively. All inlet gases enter at room temperature (300 K) and emissivity for radiation calculations is set to 1.0 (black aperture). 51 4.1.2 Outlet Boundary The gas outlet boundary (not shown in figure 4.1) has an imposed zero gradient condition for all dependent variables. It is assumed that the outlet boundary is sufficiently far downstream to ensure that the flow in the upstream region is not affected by the downstream conditions. All radiative energy coming to the outlet boundary from the flow field is reflected back into the domain. This is a reasonable assumption since in the experiments most of the exit plane is a wall. The relatively small amount of radiative energy (relative to the hot end of the kiln) reaching the exit plane would in part be reflected back into the kiln. 4.1.3 Wall Boundaries When using the K-e model with the standard wall function approach [46] a velocity profile normal to the wall is imposed (universal law of the wall) and the constant wall shear stress between the first interior computational node and the boundary is calculated. This is treated as a force in the momentum equations and is used to modify the source terms in the discretized equations. A Reynolds analogy is applied (similar velocity and temperature profiles) to calculate the constant turbulent heat flux through the wall. Values for K and e at the first interior node are also functions of the shear stress. A zero-gradient is imposed normal to the wall for all species considered. The bed surface is also modelled as a wall with the exception that a flow rate is imposed for C 0 2 52 through each control volume boundary surface area. These values are calculated by the bed model and updated for the hot flow model at every global iteration (a single run for the hot flow, bed and wall models). Most of the boundary surface area interacts with either the wall or bed geometry (at each global iteration). The fire hood however does not. In this case a one-dimensional heat conduction problem is solved (1-D refractories) for each control volume boundary surface area. In the UBC pilot kiln a 0.1524 m (6 inch) insulating wall with a conductivity of 2.5 W/mK was used in the fire hood [62]. In the numerical model this constant thermal resistance is used for the 1-D conduction calculation. An outside wall emissivity of 1.0 is assumed and heat is lost to the ambient surroundings by radiation only. On the inside of all thermally radiating walls, emissivity was set to 0.8. This is an estimate for the refractory (or refractory with lime deposit) at the experimental temperature operating range. 53 4.2 Bed Model The FEM is used to obtain a numerical solution to equations 3.30, 3.17 and 3.28 (volume, species and energy conservation respectively). The method makes use of a weighted residual formulation to arrive at a system of matrix equations. In general the differential equation [ D.E. ] is weighted by Wi as follows. JJJ^ . [D.E.\kdydz = 0 (4.3) The most common weighted residual formulation has been the Galerkin method in which the weighting Wi and the interpolation functions (variables within the differential equation) are from the same class of functions. A three-dimensional 20 node iso-parametric element (quadratic distribution for both the dependent variables and the geometry) is used in the discretization of the flow and thermal fields. Element description is provided in Appendix A. A 27-point integration is performed (3 sampling points in each direction) for the volume integrals. Gauss-Legendre quadrature integration, integrates the fifth order polynomials exactly. Matrix inversion is done by a banded Guassian elimination solver. Normally a matrix inversion would be carried out for the discretized equations that result for the entire flow field. Here matrix inversion is applied to one kiln cross section at a time (plane by plane relaxation) and is possible (as will be described below) due to a Newton's linearization of the governing equations. 54 There is generally no reason why a FEM cannot be used in conjunction with a CVM approach for the hot flow model. The interface geometries (surface areas between model domains) are identical, as will be shown in chapter 5. One of the advantages of the FEM is the ease of modelling complex geometries. This is due to a local transformation from the physical domain to computational space (each element corresponds to a locally transformed computational grid). In a CVM discretization, a global transformation must first be obtained for the entire physical domain. The FEM also has the advantage of applying boundary conditions to complex geometries with relative ease. A zero gradient, for example, in the normal direction of a curved boundary would be imposed by simply ignoring a contour integral in the formulation. A Galerkin finite element discretization results in an equivalent central differencing scheme for the differential operators. As in the CVM, upwind schemes are available. These account for highly convective flows as would be the need in a CVM calculation. Three basic techniques have been used to achieve the upwind effect in finite elements [64]: In the first, artificial diffusion is added to the physical diffusion and discretization proceeds through a conventional Galerkin finite element approach. This approach has been termed first order accurate but is highly diffusive. That is, the effective Reynolds number is reduced unless a very fine grid is used. This approach however provides the necessary dissipation for numerical stability. More recently [65] this approach has been 55 extended to a higher order scheme in which a fourth-order dissipation is recast as the difference of two Laplacian operators. In the second approach [66] the numerical quadrature rule for the convection term is modified to achieve the upwind effect. Finally in the third approach the weighting function for a given node is modified to weigh an element more heavily upwind of a node than an element downstream of the node. Such formulations are consistent Petrov-Galerkin weighted residual methods [67, 68]. The most widely used of these methods is the Streamline Upwind/Petrov-Galerkin (SU/PG) formulation [64] where an artificial diffusion operator is constructed to act only in the flow direction, minimizing any crosswind diffusion. (In multidimensional problems solutions often exhibit excessive diffusion perpendicular to the flow direction. The hybrid scheme used in the hot flow model where the flow direction may not be aligned with the grid is a case in point.) All three approaches described above suffer from either false diffusion or oscillations in the converged solution when applied to general problems that may include source terms, time dependence or with the generalization to multi-dimensions. Because of the complexity of the current problem and the major difficulties to obtain coupled solutions, the simplest approach is adopted in which artificial diffusion is simply added to the physical diffusion as required (for numerical stability). Inclusion of a higher order upwind scheme was considered risky at this point of model development. 56 4.2.1 Bed Flow Field As described in section 3.2.3 the velocity field in the axial direction (w-component) is imposed. The value is based on the ratio of volumetric flow rate of solid material to the bed inlet active layer area. (In the plug flow region no forward motion is allowed.) In the active layer this value for experimental case T21 [62] is w=-0.0038 m/s (positive z-direction follows the general hot flow direction). With no velocity gradients in the axial direction and the assumption that vorticity is constant in the entire field, the equation to be solved is (from equation 3.34), where co is the rotation rate of the wall. Equation 4.4 is solved numerically for the active layer only. In the plug flow region the velocity field will be described below. 57 4.2.1.1 Boundary Conditions Boundary conditions are shown in Figure 4.2. At the periphery of the bed a stream function value of 1.0 is used. y Wall Rotation = oo Figure 4.2: Boundary Conditions for Bed Active Layer Flow Field Model In the plug flow region the velocity field is assumed in cylindrical coordinates as ue = cor (ur = uz = 0). The single vorticity component in the flow field gives, 1 r d{ru6) dur dr dd = 20) (4.5) as expected. A scalar stream function may also be defined in cylindrical coordinates such dV dV that continuity is satisfied identically. That is ug = and - rur = . Upon dr dd substitution, equation 4.5 reduces to an ordinary differential equation, \_d_ r dr dr 2co (4.6) 58 whose solution (with *P = 1 at the wall) is, \u00C2\u00A5 = | ( r 2 - t f 2 ) + l (4.7) Equation 4.7 not only describes the stream function distribution in the plug flow region, but also supplies the boundary condition at the active layer/plug flow interface for the active layer calculation as shown in figure 4.2. 4.2.1.2 D i s c r e t i z a t i o n Equation 4.4 is solved using the FEM. The procedure for obtaining a system of matrix equations using a Galerkin weighted residual formulation is summarized here and repeated for the bed heat transfer and wall models. STEP 1: Non-dimensionalize governing equations (not done for equation 4.4). The primary reason this step is carried out prior to discretization is to be able to readily extract relative ratios of terms in the governing equations through non-dimensional numbers such as the Peclet and Bodenstein numbers. It also aids the numerical computations since most variables are now of order 1. 59 STEP 2: A Newton's linearization is applied to the dependent variables in the governing equation(s). For the stream function in equation 4.4, \u00C2\u00BBF\"+1 =VN +#F (4.8) where N represents values at the current iteration level and N +1 at the new iteration level. #F are the unknowns. STEP 3: Apply the Galerkin Method of Weighted Residuals (GMWR). This method consists of minimizing the residuals of the system of equations over the solution domain. Here only the diffusive terms are integrated by parts producing contour integrals for the Neumann boundary conditions. STEP 4: Apply a finite element discretization using the element description of Appendix A. Equation 4.4 is already linear but the standard procedure described above is followed. A Newton's linearization given by equation 4.8 gives, a 2 #p d2^ dx2 + dy2 ^ = -K (4-9) where R% is the original equation at iteration level N . dx2 + dy2 * * = V i - V i \u00E2\u0080\u0094 2 \" (4-10) 60 Convergence for the discretized system of equations can be measured directly by R% . An L x norm is used here as defined by, DOF. (4.11) where DOF is the total number of cell-vertex unknowns (element nodes) or degrees of freedom in entire domain. Applying the GMWR to equation 4.9 gives, JJJ*. d2dV d26V dx2 + dy2 dxdydz = - jjJN,R^ dxdydz (4.12) Notation e in the triple integrals indicates that integration is performed over an element. Note that volume integrals are applied to the 2-D problem considered. The 3-D element description of Appendix A was used here. The problem was solved for one element in the z-direction with zero gradients on the upstream and downstream boundaries. This was done since the subroutines required were already written for the bed thermal and wall models (3-D problems). Also note that the weight function Wi of equation 4.3 is replaced by the shape function TV. that describes the behaviour of the variable in the element. Integrating second order derivative terms of equation 4.12 by parts, JJJ dN, dN. d&\u00C2\u00A5 ox dx dy dy pxdydz = JJJ oW,. dVN | dN, dVN ^ dx dx dy dy dxdydz - j^N^V \u00E2\u0080\u00A2 ndA (4.13) 61 The zero gradient boundary conditions for the upstream and downstream boundaries are imposed by simply ignoring the contour integral in equation 4.13. ( = 0 is of course dz implicit in this equation for the 2-D problem.) Finally, substituting the finite element interpolation polynomials into equation 4.13 for both spatial discretization and unknowns #F etc. yields, NDOF NDOF Z Z JJJ 1=1 7 = 1 dN, dNj dNi dNj dx dx dy dy NDOF = Z JJI dN, d*\u00C2\u00A5N dN, dVN A, \u00E2\u0080\u009E ^ J . J , (4.14) Equation 4.14 represents the local influence matrix (of a given element). NDOF (number of degrees of freedom) is 20 for the 20-node element description used. This 20X20 local influence matrix is repeated for all elements in the domain. An element-by-element matrix assembly gives, NELEM NDOF NDOF z Z J J J ^ + 1=1 7 = 1 'dN, dN,. dN. dN, dx dx dy dy NELEM = Z e=l NDOF Z JJJ '\u00E2\u0080\u00A2=1 dN, dVN dN, dVN Ar \u00E2\u0080\u009E t J J J , \ d t ^ + \dj-dy- + N ' 2 a ) ^ (4.15) and results in a system of equations for all unknowns (CrU~L (4.26a) K P * 8 = P c \u00C2\u00B0 C \u00C2\u00B0 ^ p U ~ L (4.26b) Pe. = Pc\"co^U~L (4.26c) 1 1 1 1 \u00E2\u0080\u0094 = \u00E2\u0080\u0094 r + + (4.26d) Pe Pe Pe Pe Da3 is Damkohler's Third, which represents the ratio of heat absorbed/released by reaction to the heat transported by convection. 69 Da, = flCaC03RL CpTJJm Da, = m,C\u00E2\u0080\u009E TRL 3 Pco2 C T U P OO CK (4.27) (4.28) Equation 4.27 represents the heat absorbed by the reaction of CaC03 and equation 4.28 the heat released from the control volume with C02 in the process. In non-dimensional form the density is updated (from equation 3.12) by, p = m2 + m3

(4.30) TNH=TN+dT (4.31) with S(p and ST being the unknowns. Substituting equation 4.30 and 4.31 into 4.22 and 4.23 the linearized species transport equation becomes, Bo (4.32) where, R\u00C2\u00BB =V-V(t>N - V Bo VN + Da, (4.33) and the linearized energy equation becomes, VV\pNbT}-V-Pe \u00E2\u0080\u00A2VST \u00E2\u0080\u0094 \u00E2\u0080\u0094RT (4.34) where, R? =V-V\pNTN]-V \u00E2\u0080\u0094VTN Pe + Da, +Da, 3 \u00E2\u0080\u009E J r (4.35) 70 and R* are once again the original equations at iteration level N. Convergence for the discretized system of equations is measured through L x norms as indicated by equation 4.11. Density and source terms Dat, Da3g, Daico are evaluated at iteration level N. Applying the GMWR to equation 4.32 and 4.34 and integrating diffusive terms by parts gives, ffl N,V \u00E2\u0080\u00A2V{fy) + \u00E2\u0080\u0094 {VN,)-(Vfy) dxdydz JJJ Bo NtV \u00E2\u0080\u00A2 VN + \u00E2\u0080\u0094 (V7V,.)\u00E2\u0080\u00A2 (v>*)+ /V,Dfl, \lxdydz Bo (4.36) and, JJJ \u00E2\u0080\u00A2'7(pNSr)+j-{7Ni)-{Vdr) dxdydz N,V \u00E2\u0080\u00A2 v(pNTN)+ yJ^N, )\u00E2\u0080\u00A2 (VTN)+ N,Da3it + NtDa. co2 vixdydz + ^\u00E2\u0080\u0094\u00C2\u00A9 \u00E2\u0080\u00A2 -(b - o \u00E2\u0080\u0094 ? f t f 1 K P L A N E I K P L A N E - 1 K P L * N E * 1 Figure 4.6: Finite Element Matrix Assembly for K=KPLANE A complete set of elements in a given cross-sectional plane (KPLANE) is assembled by looking at the influence on the nodes of interest. This requires assembly from information found in KPLANE-1 through KPLANE+1. Any information for nodes outside KPLANE is ignored for the influence matrix. On the right hand side all information is assembled since this represents the original equation at 0 0<0 9<86\u00C2\u00B0 laminar: c=0.55 a =0.20 6=1.00 turbulent: Gr'x Prcos2 6 > 1010 c=0.17 a =0.25 6=2.00 6 > -86\u00C2\u00B0 laminar: c=0.55 o=0.20 6=1.00 turbulent: G^*Pr>1010 c=0.17 a=0.25 b =0.00 0>86\u00C2\u00B0 always laminar: c=0.55 o=0.20 b =0.00 0 < -86\u00C2\u00B0 always turbulent: c=0.17 o=0.25 b =0.00 Table 4.7: Correlation Constants for Equation 4.46 (-90 <0 <90) 87 Although these correlations were derived for constant heat flux conditions there is some evidence to indicate that the relations may also apply to constant temperature surfaces [73]. Since qfc in equation 4.47 is a function of the heat transfer coefficient (see equation 4.45) an evaluation of hfc for a given outer wall temperature must be obtained by iteration. For a given outer wall element surface, 3-4 iterations are typically needed for the heat transfer coefficient to be within 0.1 W/m2K. It was also found that anything more than 8 planes (elements) in the wall rotation direction had little change in the bulk/average heat transfer coefficient data (average of 8 planes). As a check, bulk heat transfer coefficient data was compared to Churchill and Chu's horizontal cylinder correlation [74] for laminar and turbulent free convection. K Nu^2 =0.60 + 0.387 5 c c 4 O a> o O 12 Inciried plants \u00E2\u0080\u00A2round a cyllndsr at constant temperature hfWAn'K) O d- i .oo \u00E2\u0080\u00A2 d-1.23 1 1 1 ' \u00E2\u0080\u00A2 1 1 1 1 1 2 3 4 5 B Heat Transfer Coefficient from Cylinder Figure 4.16: Comparison of Heat Transfer Correlations A new coefficient d in equation 4.46 accounts for the difference between this data. The new correlation for an inclined plane becomes, Nu = dc(Gr' Pr cos\" d)\" (4.50) 89 where d was found to have a value of 1.23. This new parameter may be viewed as accounting for the interaction between inclined planes on an average basis. It should be pointed out however that it is not uncommon to have an experimental data scatter of \u00C2\u00B1 20% in these empirical relations [73]. This extra step taken in comparing to a cylinder correlation simply adds confidence to the model. 4.3.1.2 Total Heat Loss to Surroundings On the outer wall heat is lost to the surroundings by the sum of free convection and radiation. With the evaluation of the free convection heat transfer coefficient the total heat flux q is the sum of equations 4.44 and 4.45. q = \u00C2\u00A3sCj(Tsl-T:)+hfc(Tsh-T\u00E2\u0080\u009E) (4.51) For the iterative computations carried out in the wall model this is recast as, q = heff{T-Tj (4.52a) where, heff = hfc + \u00C2\u00A3SCT(T2 +T2 \T + T\u00E2\u0080\u009E) (4.52b) The equation is linearized in the computations by lagging temperature T in heff by one iteration. 90 4.3.2 Discretization Equation 4.41 is discretized by following the four steps given in section 4.2.1.2. Non-dimensionalization of the governing equation is first carried out using the reference values of Table 4.8. VARIABLE NON-DIMENSIONALIZED BY Independent variables x,y,z Kiln inner diameter Din Velocities u,v,w Wall rotation rate times inner diameter Temperature 7 Ambient temperature (300 K) T\u00E2\u0080\u009E Conductivity Kw K\u00E2\u0080\u009E = 1 W/mK K\u00E2\u0080\u009E Heat transfer coefficient Ratio of K\u00E2\u0080\u009E to inner diameter KjDin Heat flux q K\u00E2\u0080\u009ETjDin Table 4.8: Reference Variables for the Non-Dimensional Wall Equation The non-dimensional equation becomes, ?P + V \u00E2\u0080\u00A2 (vj)-J-V \u00E2\u0080\u00A2 {KWVT) = 0 (4.53) dt Pew For clarity notation to indicate non-dimensional values is dropped. In the diffusive term Pew is a 'rotational' Peclet number representing the ratio of heat convection in the wall rotation direction to heat conduction. p C (OD2 p \u00C2\u00BB ( 4 5 4 ) Heat loss to the surroundings is also non-dimensionalized. Equation 4.52 becomes, q = heff{T-\) (4.55a) 91 where, heff=h/c+St(T2+llT + \) (4.55b) Once again for clarity, notation to indicate non-dimensional values is dropped. St is the Stefan number which represents the ratio of heat radiation (to surroundings) to heat conduction. (4.56) Using equation 4.31, a Newton's linearization to equation 4.53 gives, afar) at Pe.., (4.57) where, R (4.58) The original equation at iteration level N (R\") may again be used to measure convergence of the discretized system of equations through the L\u00E2\u0080\u009E norm of equation 4.11. Applying the GMWR to equation 4.57 and integrating diffusive terms by parts gives, JJJ a( are the coefficients in the manufacturer supplied curve fits for conductivity in the hot zone (Kw> - aw+bwT). Similarly aw^ and bWi are the coefficients for the insulating castable zone. fr in equations 6.2 and 6.3 eventually goes from 0 to 1 during the global iterative process. Global convergence for the hot flow/wall coupled run is shown in Figure 6.2 along with the changing values for / and f r . Relaxation factors for temperature and heat flux data transfer between models were both set to 0.02. Global Convergence Number - IG Figure 6.2: Global Convergence for Hot Flow/Wall Coupling 115 Both the change in temperature (and change in heat flux) between global iterations \"IG\" were monitored through the following norms. \"AT ERROR\"=X|r,/C-TiJC_\ (6.4) I \"Aq ERROR\"=Xk,/c-?,/C-.| <6-5) The difference between temperature (and heat flux) values between global iterations, were summed over all interface boundary points (areas). At convergence, found acceptable in this run as a 3 order of magnitude drop in these measurements, the difference in local temperature of approximately 1 degree K is present between global iterations and a change in local heat flux of approximately 0.5 kW/m2 (corresponding to change in heat rate of approximately 0.008 kW, based on an average cell interface area). 6.2 Overall Structure of Model and Convergence of Fully Coupled Run The hot flow, bed and wall models are coupled to obtain a solution for Barr's T21 experimental case for comparison [62]. Overall relaxation factors for global data transfer (temperature and heat flux) were set to 0.1. No relaxation was used for C02 transfer to the hot flow boundary from the bed surface. A flowchart of this coupling is given in Figure 6.3. An initial guess for temperature at all boundaries is once again the starting point. A run of the hot flow sub-model supplies heat flux boundary information to the wall and bed submodels. A run of the bed and then wall model results in one complete global iteration. 116 Assume Temperatures Run Hot Flow Model Run Wall Model Relax T\u00E2\u0080\u009E .T\u00E2\u0080\u009E Figure 6.3: Global Loop for Fully Coupled Run As in the hot flow/wall coupling case, all models are run to convergence in the first global iteration. After which, 100 local iterations were made in the hot flow model, 6 in the bed model and 8 in the wall model. Fifty-four such global iterations were required for a converged coupled solution. Computing times for this run are given in Table 6.2. NUMBER OF ITERATIONS CPU TIME Is' Global Iteration: Hot Flow Model Run Bed Model Run Wall Model Run 1 global iteration 20000 local iterations 20 local iterations 40 local iterations 3728 min 3600 min 83 min 45 min Global Iterations # 2 to #54 Each Hot Flow Model Run Each Bed Model Run Each Wall Model Run 53 global iterations 100 local iterations 6 local iterations 8 local iterations 2756 min total 18 min 25 min 9 min Table 6.2: Computing Times for Fully Coupled Run (Based on a 3 GHz Pentium IV Processor) 117 Global convergence for the fully coupled run is shown in Figure 6.4. Change in temperature (heat flux and C 0 2 mass flow rate) boundary data between global iterations (five parameters plotted in figure 6.4) were made through equations similar to 6.4 and 6.5. 7 6 5 \u00C2\u00BB 4 22 \u00E2\u0080\u0094\u00E2\u0080\u00A2\u00E2\u0080\u0094'\u00E2\u0080\u0094 i\u00E2\u0080\u0094 I\u00E2\u0080\u0094'\u00E2\u0080\u0094\u00E2\u0080\u00A2\u00E2\u0080\u0094'\u00E2\u0080\u0094 i \u00E2\u0080\u0094 I \u00E2\u0080\u0094 \u00E2\u0080\u00A2 \u00E2\u0080\u0094 i\u00E2\u0080\u0094 i \u00E2\u0080\u0094'\u00E2\u0080\u0094 I\u00E2\u0080\u0094 i\u00E2\u0080\u0094\u00E2\u0080\u00A2\u00E2\u0080\u0094\u00E2\u0080\u00A2\u00E2\u0080\u0094 \u00E2\u0080\u00A2 \u00E2\u0080\u0094 I \u00E2\u0080\u0094 \u00E2\u0080\u00A2 \u00E2\u0080\u0094'\u00E2\u0080\u0094'\u00E2\u0080\u0094'\u00E2\u0080\u0094 I 0 10 20 30 40 50 X/Rs Figure 7.4: Center Line Velocity Decay (Exp. from Razinsky et al. [77]) 124 Yule et al. [78] have carried out numerical computations of axisymmetric jet flows with the standard K - \u00C2\u00A3 model [46] and the algebraic stress model [79] including four dissipation source related modifications that have been proposed for the free jet. The dissipation source term is the term in the dissipation rate equation (equation 3.1) that includes the coefficient C\. Further details may be found in [78]. Computations were compared with measurements [80] at a secondary to primary diameter ratio of 22.3. One of the recommendations for general use for jet flows (with and without confinement) was the vortex stretching dissipation source correction of Pope [81]. Two important advantages in using Pope's modification for free or confined jets are that reference is only made to local events in a flow independent manner and that a physical explanation is provided to justify the modification to the K - \u00C2\u00A3 model. An additional term is included in the dissipation rate equation as follows, p\u00C2\u00A3V-^-V\u00C2\u00A3 C . G g ClP\u00C2\u00A32 C4pe2 A* K K (7.1) where in Cartesian tensor notation, X = 0),jO)jksk, (7.2) 1 K CO,, = 1 2 \u00C2\u00A3 du, duj (7.3) 1 K s\" =Tl du, du. L + 1 ydXj dx, j (7.4) The modification refers to a non-dimensional measure of vortex stretching x that is a function of the rotation tensor co,j and the rate of strain tensor stj. The terms are non-125 dimensionalized with the turbulent quantities computed in the model. As a result, in addition to a production and dissipation term in the dissipation rate equation there is now an additional term for a measure of the vortex stretching mechanism. Pope assumes that this new source for dissipation is a linear function of this measure. The coefficient C 4 has been optimized by Pope [81] to be 0.79 using experimental data of Rodi [79]. The mechanism is described below. It is argued that \"the vorticity of the large turbulent motions tends to be aligned with the vorticity in the mean flow\". That is if the mean vorticity is being stretched (in the direction of the vorticity vector) then so is the turbulent vorticity. The mechanism is clarified by considering the largest scales of turbulence. Since the rate of dissipation of turbulence energy is governed by the rate of scale reduction, only the breaking up of the large (energy containing) motions is considered. Smaller scales are only important in that eventually turbulence energy is dissipated at the smallest (Kolmogoroff) scales. What is suggested here is that the stretching of turbulent eddies or vortex tubes by the mean flow plays a large role on the rate of scale reduction. If the mean motion has a positive normal strain in the direction of the turbulent vorticity vector then the vortex will stretch. Conservation of angular momentum implies that this vortex will rotate faster and decrease in width. This vortex stretching mechanism then will on average increase the rate of scale reduction and therefore the dissipation rate. Physically, vortex stretching promotes the formation of smaller vortices. If the mean flow has a negative normal strain the vortex will rotate slower and increase in width. This would decrease the dissipation rate. 126 For an axisymmetric jet in cylindrical coordinates (with no swirl) % reduces to, K K\u00C2\u00A3J dr dx (7.5) where ux and ur are the velocity components in the axial and radial directions respectively. Equation 7.5 reveals that vortex stretching in an axisymmetric jet occurs along the circumferential vorticity component coxr = cog. This mechanism is not present in the 2-D plane jet where % - 0 \u00E2\u0080\u00A2 Computations of Yule et al. [78] were duplicated here for both the standard K-e model and with the correction of Pope. Results are given for a Ct number of 1.7 in Figure 7.5. The improvement of the center line velocity decay with Pope's correction is evident. 1.2 1.0 0.8 Ucl/Uclo 0.6 0.4 0.2 h 0.0 CENTER LINE VELOCITY Diameter Ratio = 22.3 Exp. Data [Yule+Damou, 1991] Std. k-e model Pope's correction [Ce 4=0.79] 1 * * \u00E2\u0080\u00A2 * * * * * \u00E2\u0080\u00A2 * 1 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 ' * ' ' 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 ' \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 * \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 * 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 ' \u00E2\u0080\u00A2 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 X/Rs Figure 7.5: Center Line Velocity Decay with Pope's Correction [81] (Exp. from Yule et al. [80]) 127 When applied to the experimental case of Razinsky et al. [77] however Pope's correction does not improve predictions when compared to the standard K-e model results. Figure 7.6 shows that at a diameter ratio of 6.0 Pope's correction has a negative effect at axial positions less than 14 secondary radii. 1.2 Ucl/Uclo CENTER LINE VELOCITY Diameter Ratio=6.0 o Exp.Data [Razinsky+Brighton] Std. k-\u00E2\u0082\u00AC model Pope's correction [Ce 4=0.79] Figure 7.6: Center Line Velocity Decay with Pope's Correction [81] (Exp. from Razinsky et al. [77]) The confinement level or secondary to primary diameter ratio may play an important role in the turbulent characteristics. It appears that the lower the confinement ratio the more risky it is to use Pope's correction. 128 7.3 Error Involved in Using the Standard K - \u00C2\u00A3 Model Section 7.2 considered an investigation of the prediction capability for the isothermal confined jet. The standard K - \u00C2\u00A3 model does a 'reasonable' job at predicting the mean and turbulent profdes along the length of a kiln but there is room for improvement. Non-isotropic modelling considerations have been investigated in the past for this class of turbulent flows [78, 79, 81] and further development is still required. The larger question for hot flow kiln modelling however (that would also include combustion and radiation modelling) is what would the error be in using the known limitations from cold flow predictions using the standard K \u00E2\u0080\u0094 e turbulence model. The range of error involved in using the standard K - \u00C2\u00A3 model may be revealed by adjusting the coefficient of the dissipation source term d in the dissipation rate equation. Several solutions may be obtained for different values of Ci with the aim of more accurately predicting the flow field. This range of possible C\ values may then be used to analyze the effect on hot flow modelling. Experimental data of Binder et al. [82] for the isothermal confined jet has been modelled using axisymmetric conditions. The secondary to primary diameter ratio is 10.0 with a Curtet number of 0.59. The duct wall has a 2.5 degree divergence. This experimental case has also been modelled by Zhu et al. [83]. As has been done with all modelling considered so far, a grid independence study was considered to ensure numerical accuracy. Both first and second order [84] upwind schemes were used for 3 different grids. Figure 7.7 gives the solution for the mean velocity distribution (for all six cases) at an axial position of 1.875 inlet secondary 129 diameters. Mean velocity U is non-dimensionalized with the inlet mean value. The radius is non-dimensionalized with the inlet duct radius. The results indicate that with the exception of the 26x10 volumes case the mean velocity distribution is well represented whether first or second order discretization is used. 26X10 -1st order 26X10-2nd order 51X20-1st order 51X20-2nd order 102X40-1st order 102X40-2nd order Figure 7.7: Mean Velocity Distribution at Axial Position of 1.875 Inlet Diameters for the Conical Duct of Binder et al. [82] The turbulent stress profile is also considered in Figure 7.8 which is non-dimensionalized with the square of the inlet center line velocity U 0 . u v is evaluated by, r du dv^ u V = dr dz (7.6) where u and v are the mean velocities in the axial (z) and radial (r) directions and turbulent viscosity ju, is evaluated from equation 3.2. Results of figure 7.8 once again 130 indicate that at least the intermediate grid of 51x20 volumes is required for an accurate solution. The solution given for the second order discretization used on this grid is nearly identical to the second order discretization solution on the finest grid of 102x40. As a result the intermediate grid of 51x20 volumes is used with second order discretization for further calculations. This grid is shown in Figure 7.9. Care was taken to properly grid the region around the shear layer, as indicated in the insert. i Figure 7.8: Turbulent Stress Profile at Axial Position of 1.875 Inlet Diameters for the Conical Duct of Binder et al. [82] 131 Center line Figure 7.9: Intermediate Grid of 51x20 Volumes Used for the Conical Duct Calculation of Binder et al. [82] Three solutions are given for values of Ci of 1.30, 1.44 and 1.55 in Figures 7.10, 7.11 and 7.12 respectively. Both the mean velocity and turbulent stress profdes are given at axial locations of 0.625, 1.250, 1.875 and 2.5 inlet diameters. Figure 7.11 with a value of Ci of 1.44 represents the standard K - \u00C2\u00A3 model solution. In the solution of figure 7.10 the source to dissipation in the dissipation rate equation (equation 3.1) was reduced. As a result the larger turbulent kinetic energy (when compared to the standard K - \u00C2\u00A3 model solution) results in an increased core flow spread rate. The higher radial momentum transfer quickly reduces center line velocities. The opposite effect appears in figure 7.12. 132 (a) Z/Do=0.625 (b)Z/Do= 1.250 (c) Z/Do= 1.875 10 15 lOOOuv/Uo' (d) Z/Do=2.500 Figure 7.10: Mean velocity (left) and turbulent stress (right) solution with Ci=1.30 Experimental data from Binder et al. [82] 133 (a) Z/Do=0.625 (b)Z/Do= 1.250 (c)Z/Do= 1.875 lOOOuv/Uo' 10 15 lOOOuv/Uo2 (d) Z/Do=2.500 Figure 7.11: Mean velocity (left) and turbulent stress (right) solution with Ci=1.44 Experimental data from Binder et al. [82] 134 (a) Z/Do=0.625 (b)Z/Do= 1.250 (c)Z/Do= 1.875 10 15 lOOOuvflJo* (d) Z/Do=2.500 Figure 7.12: Mean velocity (left) and turbulent stress (right) solution with Ci=1.55 Experimental data from Binder et al. [82J 135 The mean velocity distribution appears to be more accurately represented with a Ci value of 1.30 in the early stages of jet spreading (Z/Do=0.625) as shown in figure 7.10a. By axial position Z/D0=2.5 the mean velocity distribution is more accurately modelled with a value of Ci of 1.55 as shown in figure 7.12d. Similarly, the turbulence stress is modelled accurately with the standard C i value of 1.44 at an axial position of Z/Do=0.625 as indicated in figure 7.1 la. By axial position Z/D0=2.5 the turbulence stress would require a value of C i between 1.44 and 1.55. The important point to be made here is that these three solutions represent the upper and lower limits on any simple turbulence model improvement through changes in constant C i . This range for C i (1.30 - 0 . 2 i i i i I i i i i I i i i 11 i i i 11 i i i i I i i i i I i i i i I i i i 11 i i i 11 i i i i I i i i i > i -0.2 0 0.2 0.4 O.B 0.8 1 1.2 1.4 1.6 1.8 2 Axial Direction (rn) Figure 8.1: Fuel (CH4) Mass Fraction (length scale distorted) The axial velocity component shown in the contour plot of Figure 8.2 identifies various details of the calculated flow, including the Craya-Curtet type recirculation pattern downstream of the flame. Reverse flow in the burner indicates possible fouling problems. This was indeed observed during the experimental trials on the UBC pilot kiln. Recirculation in the completely sealed hood is also evident. 144 -0.4 -0.2 0 0.2 0.4 O.B 0.8 1 1.2 1.4 1.6 Axial Direction (rn) Figure 8.2: Axial Velocity Field (m/s) (length scale distorted) Two hot spots are indicated in the inner wall temperatures of Figure 8.3. The first is the hot spot due to flame impingement. Downstream of the flame zone a second hot spot is visible. This is present due to a refractory change to lower conductivity. As conduction drops, heat flux through the wall drops but inner wall temperatures increase. (At 1000 K conductivity drops from 1.685 W/mK to 0.275 W/mK.) Such hot spots can severely reduce the life expectancy of the refractory lining. In addition, these temperature profiles have implications for product quality (the possibility of overburned lime from rapid calcination at high temperatures) and the possibility of ring formation in the hot end. It is evident from these observations that proper design of the refractory lining requires a detailed knowledge of the wall thermal profile. 145 Figure 8.3: Wall Hot Spots in Pilot Kiln (Inner Wall Temperatures in K) A comparison with the experimental data is shown in Figure 8.4. The numerical data (lines) capture the radial thermal gradients of the buoyancy-affected field as well as the expected difference between upper and lower inner wall temperatures. 8 % 1 2 3 4 5 Axial Position (m) Figure 8.4: Axial Temperature Profiles in Pilot Kiln (Exp. Run#6 Alyaser [61]) 146 The largest deviations between the experimental data and the computed results appear to be in the flame zone of the hot flow results. During the global iterative process these 'lines' oscillated above and below the experimental data. This effect is investigated further in section 8.1.1. A closer look at the temperature field in the flame zone is shown in Figure 8.5. The effects of the buoyant flame are visible in both the central vertical plane and the two cross-sectional planes (located at 0.88 and 1.78 kiln diameters downstream of the burner exit plane). The hottest region in the flow field is at the periphery of the core flow exiting the burner. This is expected since this is where most of the combustion takes place. Within the core flow a fuel rich zone is typically outside the flammability limits for combustion. Similarly, outside the core, fuel lean conditions exist. Y Figure 8.5: Hot Flow Temperature Contours in Flame Zone (K) 147 Inner wall temperatures of the burner/fire hood region are shown in Figure 8.6. Due to the buoyant flame the temperature at the fire hood TDC (top dead center) position is about 30 K hotter than that at BDC. The much hotter burner surface and cooler inlet pipes are also evident. Y Figure 8.6: Inner Wall Temperatures of Pipes, Burner and Fire Hood (K) 148 8.1.1 A Transient Solution of the Hot Flow In the flame zone of figure 8.4 the predicted hot flow temperature curves oscillated above and below the experimental data during the global iterative process. This indicated that there might be an unstable flame present under these conditions. A transient solution of the hot flow (with a fixed wall temperature) was then calculated. The time step (0.1 sec) was chosen such that information traveled one control volume for each time step taken. The computations were carried out for a total real time of 10 seconds. Figures 8.7 and 8.8 show the transient solution at 2-second intervals on the central vertical plane and the cross-sectional plane (located at 1.78 kiln diameters from the burner exit) respectively. Results for fuel mass fraction (CH4) and temperature are given. An animation of these results was made using the 0.1-second interval. Results show that although the flame is located above the center line of the kiln during most of this time period there is a great deal of unsteadiness in its behaviour. The flame repeatedly impinges onto the upper refractory, as well as being situated below the center line of the kiln during brief time periods. In the fuel mass fraction plots of figure 8.7 the flame is shown to impinge the upper refractory leaving fuel to burn in this region as the core jet moves away from the wall. The thermal plots show that this effect can have a 400 K difference at the upper inner wall boundary over this time period. An analysis of the animation video shows that parts of the upper wall see as much as a 600 K oscillation in temperature. Figure 8.8 shows that oscillations are not only present in the vertical direction. Temperature contours closely follow the flame boundary. 149 Figure 8.7: Fuel Mass Fraction (left side) and Temperature ( K ) (right side) on Central Vertical Plane (Vertical axis is amplified 3X) 150 (d) Time = 6.0 seconds (e) Time = 8.0 seconds (f) Time = 10.0 seconds Figure 8.7 (cont'd): Fuel Mass Fraction (left side) and Temperature (K) (right side) on Central Vertical Plane (Vertical axis is amplified 3X) 15 1 (a) Time = 0.0 seconds (b) Time = 2.0 seconds (c) Time = 4.0 seconds Figure 8.8: Fuel Mass Fraction (left side) and Temperature (K) (right side) at 1.78 Kiln Diameters Downstream of Burner Exit (d) Time = 6.0 seconds (e) Time = 8.0 seconds (f) Time = 10.0 seconds Figure 8.8 (cont'd): Fuel Mass Fraction (left side) and Temperature (K) (right side) at 1.78 Kiln Diameters Downstream of Burner Exit 153 This instability in the flame is not necessarily limited to hot flow conditions. It is well known that cold flow confined jets at low Ct numbers display an increased unsteadiness in the flow field [6]. 8.2 Fully Coupled Results Coupled model results are described for the thermal measurements of Barr [62] on the UBC 0.41 m internal diameter pilot kiln. These measurements include gas flow temperatures 10 cm off the wall and 2.5 cm off the bed. These locations are shown in the insert of Figure 8.9. Measurements were also made at four radial locations on the refractory (8 locations along the length of the kiln). Bed temperature was measured at 12 axial locations. Since the thermocouples alternate between the freeboard gas and bed, the temperature was calculated through a lumped capacity response equation using the transient signal. The bed temperature may be taken as the bulk (cross-sectional average) temperature. Run T21 of Barr's experiments was modelled with a burner load of 88 kW and a primary air ratio of 0.4. Conditions for this run are given in Table 8.4. The Curtet number was calculated as was done in section 8.1. Compared to Alyaser's run #6 it is 12 % higher for this run. This indicates that there may be less oscillatory behaviour in the flame zone. Indeed, no instabilities were observed in these hot flow results. It is important to note that the reaction rate used in this model was calibrated to produce the experimentally observed degree of calcination at the bed exit. The control 154 parameter for the bed depth distribution of the calcination temperature in the reaction rate model, as described in section 4.2.2.5, is adjusted to reflect the degree of calcination at the bed exit. It is hoped that the value of this constant would predict the degree of calcination for other industrial kilns but further checks should be carried out first. Natural Gas Flow 5.5 SCFM Burner Load 88 kW Total Primary Air Flow 36.0 SCFM Total Secondary Air Flow 54.0 SCFM Excess Air 72% Ct number 0.73 Limestone Feed Rate 55Kg/hr Wall Rotation 1.5 rpm Table 8.4: Pilot Kiln Trial Run T21 from Barr [62] Heat transfer rate imbalance data in the coupled model are given in Table 8.5. An overall heat transfer rate imbalance of less than 1% of the burner load is present in the converged results. The mass flow rate imbalance in the hot flow model is 1.3% of the total input flow rate. In the bed model the mass flow rate imbalance is 4.5% of the bed input feed rate. This relatively larger error along with the 1.4 kW in table 8.5, which is 5.4% of the heat absorbed by the calcination reaction, are primarily caused by the diffusive flux calculations at the bed boundary as was explained in chapter 4. It is important to note that minimizing or eliminating false diffusion in the bed calculations would reduce these errors. Model Imbalance Percent of Burner Load Hot Flow 1.2 kW 1.4% Bed 1.4 kW 1.6% Wall 0.8 kW 0.9 % Table 8.5: Model Heat Rate Imbalances for Fully Coupled Results 155 An overall heat balance model was developed and selected values compared to the CFD results. Details of the heat balance model are given in Appendix F in the form of a spreadsheet. Input to this model includes inlet and outlet temperatures, mass flow rates, the 2 kW loss of heat through the fire hood region and the degree of calcination. This boundary information was taken from the CFD data. The aim here was to compare the 3-D computed heat transfer rate data with a global control volume calculation. A comparison of the calculated heat rates from this model to the CFD results is given in Table 8.6. The satisfactory comparison gives confidence in the computed 3-D fields. Heat Balance Model 3-D CFD Result Difference Heat carried with products of combustion 54.7 kW 55.1 kW 0.4 kW Burner Thermal Load 88.1 kW 89.0 kW 0.9 kW Heat Absorbed by Calcination Reaction 24.9 kW 25.7 kW 0.8 kW Heat loss through shell 15.8 kW 16.8 kW 1.0 kW Table 8.6: Comparison of Heat Rates to Overall Heat Balance Model of Appendix G 8.2.1 Comparison to Available Experimental Data Axial temperature profiles in the pilot kiln are compared to Barr's trial T21 in Figure 8.9. The numerical results (lines) are generally less than 100 K cooler than the experimental data. Due to the high value of excess air, the flame is only about 0.5 m in length (based on 1% fuel mass fraction contour results). The flame is buoyant and in this region the gas 156 temperatures 10 cm off wall are as much as 500 K higher than that at 2.5 cm off bed. Since the measurements started at about 0.5 m this difference in temperature was not captured in the experiment. Aside from the first 0.5 m the two hot flow curves are very similar over most of the kiln length. Differences are once again present in the cold end of the kiln with a maximum difference at the kiln exit of about 25 K. The gas temperatures 2.5 cm off the bed are cooler in this region (compared to the 10 cm off the wall data) likely due to the heat absorption by the cold bed feed. The high heat rates in the cold end result in an increase in slope (in absolute terms) for both hot flow curves. 2 3 4 Axial Position (m) Figure 8.9: Axial Temperature Profiles in Pilot Kiln (Exp. Run T21 Barr [62]) The overall heat transfer rate from the hot flow zone to the wall is 40.5 kW, 78% of which is by radiation. Similarly, the overall heat transfer rate from the hot flow zone to the bed surface is 17.2 kW, 81% of which is by radiation. 157 Bulk bed temperatures and inner wall temperatures are also very similar to one another in the middle portion of the kiln. These curves indicate that the kiln may be divided into three distinct zones. In the cold end a sharp increase in bed and wall temperature is observed. The high rate of heat input to the bed quickly diminishes at about 4.2 m. From this location to about 1.5 m the bed and inner wall temperatures are indistinguishable. In the final quarter length of the bed (hot end) wall temperatures are higher than the bulk bed temperatures. Bed temperatures are nearly constant in this region. The sudden change in slope of the bed temperature (and to a lesser extent the wall temperature) is due to the onset of the calcination zone. Due to the endothermic reaction, high rates of heat input to the bed are also present in this region. Figure 8.10 compares bed heat transfer rates to experimentally determined data. Both the heat input from the covered wall to the bed and the net heat input indicate that high rates are only present in the initial heating of the bed material and in the calcination zone. The middle portion (half of the kiln length) is highly inefficient in transferring heat to the bed. In the calcination zone the predicted heat transfer rates are considerably higher than those reported in the experimental data. A possible reason for the larger differences present in the hot end will be given in section 8.2.3.1. A comparison between the calculated and experimentally observed heat loss through the shell is also made in Figure 8.11. Once again the largest deviation between the predicted values and the experimental data occurs in the calcination zone but in general the comparison is favorable. 35 Axial Position (m) Figure 8.10: Bed Heat Transfer Rates (Exp. Run T21 Barr [62]) Figure 8.11: Heat Loss Through Shell (Exp. Run T21 Barr [62]) 159 8.2.2 Hot Flow Results Hot flow temperature contours are shown in Figure 8.12. The vertical plane in this view is a slanted plane as indicated in the insert. The two cross-sectional planes shown are once again at 0.88 and 1.78 kiln diameters downstream from the burner exit. Bed surface temperature results in this region are also included. Buoyancy effects are not as large here when compared to the hot flow/wall coupled results. To quantify this, an estimate of the buoyancy force to the inertial force (Grashoff number divided by Reynolds number squared, GWRe2, effectively proportional to the inverse of a Froude number) was made for both cases. The value of Gr /Re2 for the hot flow/wall coupled case is twice that of the hot flow/bed/wall fully coupled case described here. These calculations were based on the kiln diameter, a bulk wall temperature, an adiabatic flame temperature and a mean core flow velocity. The insert in figure 8.12 gives a closer look at the calculated temperatures at the 0.88 kiln diameter cross-section plane. In this view the flame is deflected vertically 12% of the kiln radius (based on fuel mass fraction data). It is important to note that in this run no fluctuations were observed in the flame zone. Inner wall temperatures at the fire hood region are shown in Figure 8.13. Temperatures at the top of the fire hood are once again about 30 K hotter than those at BDC. 161 Figure 8.13: Inner Wall Temperatures (K) in the Fire Hood Region for Fully Coupled Results 8.2.3 Bed Results The temperature distribution at the bed surface is shown in Figure 8.14. As indicated by the bulk bed temperature profile of figure 8.9 the hot end is at a roughly uniform temperature of about 1090 K. 162 Y Figure 8 .14 : Bed Surface Temperature Distribution (K) Temperature profiles at three bed cross-sections are given in Figure 8.15. Axial positions of 0.5, 2.5 and 4.5 m are representative of the calcination zone, the low and high heat transfer rate zones to the bed, respectively. In all three results the highest thermal gradients (and highest temperatures) occur at the initial wall heat regeneration position (lower left portion of wall). As will be seen below a large portion of the calcination reaction takes place in this region. The temperature range in each of the three cross-sections is different. In the calcination zone plot a 40 K range is present. In the low heat transfer zone plot only a 7 K range is observed. In the cold end of figure 8.15c the temperature range increases once again to 30 K. These numbers are proportional to the local heat transfer rates. In a large portion of the plug flow region there are small circumferential (re-direction) thermal gradients due to the highly convective nature of the flow field. That is not at all the case in the active layer as observed in figure 8.15a. (a) Axial Position 0.5 m (c) Axial Position 4.5 m Figure 8.15: Cross-Section Bed Temperature Profiles (K) at Axial Positions of (a) 0.5 m (b) 2.5 m and (c) 4.5 m 164 The solids volume fraction of CaCC\u00C2\u00BB3 in the bed is given in Figure 8.16. The insert in the figure gives the averaged cross-sectional value along the length of the kiln. The calcination zone occupies a short section in the hot end that makes up only about 25% of the total kiln length. Bed flow direction 8 o c / l M i 3 o a* in - 1 3 5 03 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 D 2 S 4 Aid a) CXra atari (m) s Figure 8.16: Solids Volume Fraction of CaC0 3 in Bed The cross-section solids volume fraction of CaC03 at three axial positions is shown in Figure 8.17. The highly convective nature of the plug flow region is once again apparent in all three plots. It is important to note that the range of volume fractions in each cross-section is less than 0.01. At the axial position of 0.0 m (exit plane) an averaged value of 0.095 indicates the overall degree of calcination for the process. 165 (b) Axial Position 0.5 m (c) Axial Position 1.3 m Figure 8.17: Cross-Section Solids Volume Fraction of CaCC>3 at Axial Positions of (a) 0.0 m (b) 0.5 m and (c) 1.3 m 166 A surface contour plot for the reaction rate is shown in Figure 8.18. A value of 0.0005 per second is plotted as calculated by equation 3.41. A perspective view is given looking towards the hot end from within the bed. The plot reveals in detail where calcination takes place. Two distinct zones are identified. First, calcination takes place primarily in the active layer where high heat transfer to the bed is present from the hot flow. A second zone is present in the initial wall heat regeneration (left wall) region. The reaction rate data in this region (sum of data along all elements attached to the left wall) reveals that 28% of the calcination occurs due to direct wall heat regenerative action. Figure 8.18: Surface Contour for the Reaction Rate at a Value of 0.0005 per sec f 167 A detailed plot of the reaction rate at two cross-sections is given in Figure 8.19. At axial position 0.5 m calcination does occur in part in the plug flow region. Figure 8.19b is however more typical of the reaction zone. That is, the majority of the reaction occurs in the active layer. (b) Axial Position 1.5 m Figure 8.19: Reaction Rate Contours (per second) at Axial Positions of (a) 0.5 m and (b) 1.5 m 168 The current reaction rate model used in these calculations is realistic, as can be shown by comparison with experimental data of Murthy et al [89]. In Murthy's experiments, cylindrically shaped pellets of three different diameters were cooked at a temperature of 1113 K. Plots giving the degree of calcination versus time were given for the three diameters of 8.95 mm, 11.95 mm and 17.0 mm. To check that the reaction rate model used in this thesis is of the same order of magnitude as the data provided by Murthy et al., the time for complete calcination of the pellets was considered. This gives bulk values of reaction rates for the experiments. Assuming a partial pressure of 0.15 and therefore a calcination temperature of 1053.5 K the maximum difference between the experimental bulk reaction rates and the current model is 40% in the 3 sets of data. Values used for the conductivity of CaO etc. in equation 3.41 were those used in the present CFD calculation. It should be emphasized that this very rough comparison is made simply to show that the reaction rate model used here, whose constants were chosen to provide the observed degree of calcination, agrees reasonably well with other available data. 8.2.3.1 Effect of Artificial Viscosity As mentioned in chapter 4 artificial viscosity was used to stabilize the numerical scheme used for the bed calculation. Figures 8.20 and 8.21 give the artificial Peclet number distributions used in the active layer and plug flow region respectively. This value is constant in each respective region at any cross-section. Artificial 169 viscosity/conductivity in the plug flow region is generally one order of magnitude lower than that in the active layer. As shown in figure 8.20 the artificial Peclet number is about one order of magnitude lower (artificial conductivity is one order of magnitude higher) than the minimum granular Peclet number in the cold end. In the hot end this difference is about 2.5 orders of magnitude. 6 r fij I a- h Q_ n r i i i i i i i i i i i i i i i \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 i i i i \u00E2\u0080\u00A2 i \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 1 2 3 4 5 Axial Po8ltlon (m) Figure 8.20: Peclet Numbers in Active Layer In the plug flow region a similar situation occurs as shown in figure 8.21. The artificial Peclet number is about one order of magnitude lower than the bulk bed Peclet number in the cold end and about 2 orders lower in the hot end. 170 LU o w 1 4 n. 3 c 3 z 5 Q_ Bed Bulk Peclet Number - PE ' i i i i I i i Artificial Peclet Number -PE . 1 1 1 1 1 1 1 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 1 1 \u00E2\u0080\u00A2 1 \u00E2\u0080\u00A2 1 2 3 4 Axial Position (m) Figure 8.21: Peclet Numbers in Plug Flow Region Artificial viscosity was also used for the species equation. In the active layer an artificial Bodenstein value of 500 was used for the entire active layer region. This is about one order of magnitude lower than the minimum granular Bodenstein number. In the plug flow region no viscosity effects should be present for the species equation but an artificial Bodenstein value of 104 was required to stabilize the scheme. Artificial viscosity generally affects the solution in two ways. First, the thermal and volume fraction fields tend to be highly diffusive. The range of temperatures for example would be lower in the present results than if no artificial diffusion was used. Granular effects in the active layer are somewhat smeared in these solutions. Second, heat fluxes generally increase to or from the bed. In the bed heat transfer rates plot of 171 figure 8.10 for example the heat transfer rate to the bed was much higher than the experimental data in the hot end. Figures 8.20 and 8.21 may explain this situation since artificial viscosity is highest in the hot end. 8.2.4 Wall Results Inner wall temperatures are shown in Figure 8.22. Due to a high rotational component of heat transfer, temperature gradients in the circumferential direction are minimal compared to the hot flow/wall coupled results where no wall rotation was present. Figure 8.22: Inner Wall Temperature (K) This is observed throughout the wall thickness in Figure 8.23. Circumferential gradients are only present in what has been termed the active wall layer [13], which is roughly the inner 5% of the wall region. A closer look at the inner wall circumferential Y 172 temperature distribution is given in Figure 8.24. The inner wall temperature increases in the hot flow region and then releases heat to the bed periodically. Figure 8.23: Cross-Section Wall Temperature Distribution (K) at Axial Position of 2.5 m 173 1000 exposed wall covered wall 900 -P 5 8 0 0 f Axial direction 4.5 m CO CD CL IB - r ^ 600 -CD C c 500 h Axial direction 5.2 rn 400 Circumferential direction Figure 8.24: Inner Wall Temperature as a Function of Circumferential (0) Direction 8.2.5 Dominant Thermal Diffusion Directions and Possible Simplifications to the Models Components of the temperature gradient term, VT, were analyzed for the wall model in cylindrical coordinates. This analysis supports above observations by showing that the circumferential direction component is 1 to 2 orders of magnitude smaller than the radial component throughout the wall domain. In addition, the axial component is also about 2 orders of magnitude smaller than the radial component. This indicates that it may be sufficient to solve the wall as has been done in the past [19, 20, 39] by using a one-dimensional approach with the equation, (8.1) 174 A similar analysis was carried out using the bed data. Both temperature and volume fraction gradients were found to have no dominant direction. A definite conclusion however cannot be made here since artificial viscosity has large effects on these gradients in both the thermal and species fields. 8.3 S u m m a r y Fully 3-D coupling of the physical phenomena and hardware interactions has allowed a detailed investigation of the processes present in a rotary lime kiln. A summary with concluding remarks follows in chapter 9. The bed model is available from the Department of Mechanical Engineering at the University of British Columbia. 175 9.0 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS The hot flow, bed and wall numerical models for a lime kiln have been coupled and used to model selected experimental conditions of the UBC pilot kiln trials [61, 62]. The physical processes taking place in the rotary kiln were investigated using the wealth of information available from the computed results. The solutions indicate that the present model may be used to identify problems of kiln operation or design. Prediction of a wide range of information is possible for various (primary air/fuel) input conditions, kiln geometries and burner designs. Analysis of the results leads to an improved understanding of the process. When the flow field, temperature distributions and species concentrations are visualized and understood for a given kiln operating condition, optimization and diagnosis of operational problems can be greatly enhanced through the use of the present computational tool. Analyzing complex aerodynamic flows (secondary flows) of hood/kiln interactions through flame visualization (and possible impingement) also has the potential of increasing refractory life. A brief summary, concluding remarks and recommendations for future work are given below: Summary: \u00E2\u0080\u00A2 A hot flow solution has been formulated using an existing CFD solver [45] that captures details of the kiln, burner and fire hood interactions. 176 The accuracy of the hot flow model was investigated by considering the aerodynamics in the hot end of a rotary lime kiln. Flame impingement on the refractories (caused by the buoyancy effect) was observed in the hot flow/wall coupled results. This is important because resulting wall hot spots can severely reduce the life expectancy of the refractory lining. The oscillatory behaviour of the flame was observed and shown to produce large temperature oscillations near the walls of the kiln. A bed model has been developed that includes bed motion, species and thermal transport as well as the calcium carbonate chemistry specific to lime kilns. The three dimensional non-Newtonian field was approximated and modelled as two distinct Newtonian regions: an active layer adjacent to the hot flow and a plug flow region below. The bed flow field was supplied via a vorticity/stream function formulation. Granular motion was modelled and assumed to follow diffusive type behaviour. Granular diffusivity and thermal conductivity are functions of a locally imposed granular temperature. A bed reaction rate model was developed and assumed to be dominated by heat transfer. Diffusion effects of CO2 were accounted for by imposing a partial pressure distribution as a function of bed depth. Although diffusion effects of CO2 have been approximated, the bed reaction rate model has been calibrated to an experimentally observed overall degree of calcination. The model has also been shown to agree reasonably well with other experimental observations and may be useful in future predictions of calcination in lime kilns. A wall model has been developed that includes rotation and losses to the ambient. 177 \u00E2\u0080\u00A2 A plane-by-plane relaxation technique has been adapted in the finite element numerical procedure for the solutions of the bed and wall models. This was made possible by applying a Newton's linearization to the governing equations. In this approach computer memory requirements are drastically reduced. Concluding remarks: \u00E2\u0080\u00A2 An investigation of the prediction capability for the isothermal confined jet showed that the standard K-e model does a 'reasonable' job at predicting core flow spread rates, mean velocity and turbulent stress distributions. Although there is room for improvement in the turbulence model, there is still no universally accurate turbulence simulation method for this class of problem. The errors involved however in using the standard K-e model in lime kiln modelling with combustion and radiation has been shown to be acceptable. The possible error in the position of a wall hot spot for example may be limited to within one kiln diameter. \u00E2\u0080\u00A2 Buoyancy effects have been shown to be important for both the generic industrial kiln considered in chapter 7 and for the coupled solutions of chapter 8. \u00E2\u0080\u00A2 For the coupled solution of chapter 8, roughly 80 % of the heat transfer from the hot flow to the bed and wall occurred through radiation. 178 \u00E2\u0080\u00A2 Two calcining zones are identified from the bed reaction calculations. Calcination takes place primarily in the active layer where high heat transfer to the bed is present from the hot flow. A second zone is present in the initial wall heat regeneration region. For the fully coupled case described in chapter 8, 28 % of the calcination occurs due to direct wall heat regenerative action. \u00E2\u0080\u00A2 Analysis of the 3-D results of the rotating wall indicates that it may be sufficient to solve the wall by using a one-dimensional approach. \u00E2\u0080\u00A2 The fully coupled solution (hot flow, bed and wall) suggests that the middle portion of the UBC pilot kiln (about half of the kiln length) is highly inefficient in transferring heat to the bed. High heat transfer rates are only present in the initial heating of the bed material and in the endothermic calcining region. Recommendations: o A higher order upwind scheme for the finite element computations in the bed should be included in the future. This would prevent highly diffusive solutions and the generally increased heat fluxes to the bed. o A one-dimensional approach using equation 8.1 for the wall model would reduce complexity and computer time. o Investigate the possibility of extending the model to the chain zone. 179 BIBLIOGRAPHY 1. PERAY, K.E., WADDELL, J.J., \"The Rotary Cement Kiln\", Chemical Publishing Co., Inc., New York (1972). 2. KRAMM, D.J., \"Update on Lime Sludge Kilns in the Pulp Mill Environment\", Paper Trade Journal, 52-58, May (1979). 3. RUHLAND, W., \"Investigation in the Cement Rotary Kiln\", Journal of the Institute of Fuel, 40:69-75 (1967). 4. BUI, R.T., SIMARD, G., CHARETTE, A., KOCAEFE, Y., PERRON, J., \"Mathematical Modelling of the Rotary Coke Calcining Kiln\", Can. J. Chem. Eng., 73:534-544 (1995). 5. RICOU, F.P., SPALDING, D.B., \"Measurements of Entrainment by Axisymmetrical Turbulent Jets\", J. Fluid Mech., 11:21-32 (1961). 6. MOLES, F.D., WATSON, D., LAIN, P.B., \"The Aerodynamics of the Rotary Cement Kiln\", Journal of the Institute of Fuel, Paper 6, 353-362, Dec. (1973). 7. THRING, M.W., NEWBY, M.P., \"Combustion Length of Enclosed Turbulent Jet Flames\", 4th Int. Symp. on Comb., 789-796 (1953). 8. HINZE, J.O., \"Turbulence\", McGraw-Hill, Inc., second edition (1975). 9. CRAY A, A., CURTET, R., \"On the Spreading of a Confined Jet\", C.R., Acad. Sci., Paris, 241:621-622(1955). 10. CURTET, R., \"Confined Jets and Recirculation Phenomena with Cold Air\", Combustion and Flame, 2:383-411 (1958). 11. BECKER, H.A., HOTTEL, H.C., WILLIAMS, G.C., \"Mixing and Flow in Ducted Turbulent Jets\", 9th Int. Symp. on Comb., 7-20 (1963). 12. GOROG, J.P., BRIMACOMBE, J.K., ADAMS, T.N., \"Radiative Heat Transfer in Rotary Kilns\", Met. Trans. B, 12B:55-70 (1981). 13. BARR, P.V., BRIMACOMBE, J.K., WATKINSON, A.P., \"A Heat Transfer Model for the Rotary Kiln: Part I. Pilot Kiln Trials\", Met. Trans. B, 20B:391-402 (1989). 14. PEARCE, K.W., \"A Heat Transfer Model for Rotary Kilns\", Journal of the Institute of Fuel, Paper 7, 363-371, Dec. (1973). 180 15. SUNAVALA, P.D., \"A Mathematical Analysis of the Temperature Distribution of Gas and Charge in Dry Process Rotary Cement Kilns\", Journal of the Institute of Fuel, 23-32, Mar. (1977). 16. CROSS, M., YOUNG, R.W., \"Mathematical Model of Rotary Kilns Used in the Production of Iron Ore Pellets\", Ironmaking and Steelmaking, 3:129-137 (1976). 17. WATKINSON, A.P., BRIMACOMBE, J.K., \"Heat Transfer in a Direct Fired Rotary Kiln: II. Heat Flow Results and Their Interpretation\", Met. Trans. B, 9B:209-219 (1978). 18. BRIMACOMBE, J.K., WATKINSON, A.P., \"Heat Transfer in a Direct Fired Rotary Kiln-I. Pilot Plant and Experimentation\", Met. Trans. B, 9B:201-208 (1978). 19. BARR, P.V., BRIMACOMBE, J.K., WATKINSON, A.P., \"A Heat Transfer Model for the Rotary Kiln: Part II. Development of the Cross-Section Model\", Met. Trans. B, 20B:403-419(1989). 20. GOROG, J.P., ADAMS, T.N., BRIMACOMBE, J.K., \"Regenerative Heat Transfer in Rotary Kilns\", Met. Trans. B, 13B: 153-163 (1982). 21. GOROG, J.P., ADAMS, T.N., BRIMACOMBE, J.K., \"Heat Transfer from Flames in a Rotary Kiln\", Met. Trans. B, 14B:411-424 (1983). 22. SOOD, R.R., CLARK, R., STOKES, D.M., \"Computerized Mathematical Model to Determine the Optimum Design and Operating Parameters for Rotary Kilns for Petroleum Coke Calcination\", Light Met., Proc. Tech. Sess. AIME 101st Annual Meeting, Denver, CO. 151-161 (1972). 23. LI, K.W., FRIDAY, J.R., \"Simulation of Coke Calciners\", Carbon, 12:225-231 (1974). 24. BROOKS, D.G., \"Mathematical Simulation of a Rotary Coke Calciner\", Light Met., Proc. Tech. Sess. TMS 118th Annual Meeting, Las Vegas, NV., 461-469 (1989). 25. SILCOX, G.D., PERSHING, D.W., \"The Effects of Rotary Kiln Operating Conditions and Design on Burden Heating Rates as Determined by a Mathematical Model of Rotary Kiln Heat Transfer\", J. Air Waste Manage. Assoc., 40:337-344 (1990). 26. BUI, R.T., PERRON, J., READ, M., \"Model-Based Optimization of the Operation of the Coke Calcining Kiln\", Carbon, 31 (7): 1139-1147 (1993). 27. GHOSHDASTIDAR, P.S., ANANDAN UNNI, V.K., \"Heat Transfer in the Non-Reacting Zone of a Cement Rotary Kiln\", J. Engin. Industry, 118:169-172 (1996). 28. HENEIN, H., BRIMACOMBE, J.K., WATKINSON, A.P., \"Experimental Study of Transverse Bed Motion in Rotary Kilns\", Met. Trans. B, 14B:191-205 (1983). 181 29. RUTGERS, R., \"Longitudinal Mixing of Granular Material Flowing Through a Rotating Cylinder- Part I. Descriptive and Theoretical\", Chem. Eng. Science, 20:1079-1087 (1965). 30. MU, J., PERLMUTTER, D.D., \"The Mixing of Granular Solids in a Rotary Cylinder\", AIChE Journal, 26(6):928-934 (1980). 31. CARLEY-MACAULY, K.W., DONALD, M.B., \"The Mixing of Solids in Tumbling Mixers-I\", Chem. Eng. Science, 17:493-506 (1962). 32. CARLEY-MACAULY, K.W., DONALD, M.B., \"The Mixing of Solids in Tumbling Mixers-II\", Chem. Eng. Science, 19:191-199 (1964). 33. TSCHENG, S.H., WATKINSON, A.P., \"Convective Heat Transfer in Rotary Kilns\", Can. J. Chem. Eng. 57:433-443 (1979). 34. JENKINS, B.G., MOLES, F.D., \"Modelling of Heat Transfer from a Large Enclosed Flame in a Rotary Kiln\", Trans. Insti. Chem. Eng., 59:17-25 (1981). 35. HOTTEL, H.C., SAROFIM, A.F., \"Radiative Transfer\", McGraw-Hill, Inc. (1967). 36. BOATENG, A.A., BARR, P.V., \"A Thermal Model for the Rotary Kiln Including Heat Transfer Within the Bed\", Int. J. Heat Mass Transfer, 39(10):2131-2147 (1996). 37. BOATENG, A.A., BARR, P.V., \"Modelling of Particle Mixing and Segregation in the Transverse Plane of a Rotary Kiln\", Chem. Eng. Science, 51(17):4167-4181 (1996). 38. AL YASER, A.H., \"Fluid Flow and Combustion in Rotary Kilns\", PhD Thesis, Faculty of Graduate Studies, Department of Metals and Materials Engineering, University of British Columbia (1998). 39. MASTORAKOS, E., MASSIAS, A., TSAKIROGLOU, CD., GOUSSIS, D.A., BURGANOS, V.N., PAYATAKES, A.C, \"CFD Predictions for Cement Kilns Including Flame Modelling, Heat Transfer and Clinker Chemistry\", Appl. Math. Model, 23(1):55-76(1999). 40. CFDS, FLOW-3D Users Manual, AEA Harwell, UK. 41. BUI, R.T., SIMARD, G., CHARETTE, A., KOCAEFE, Y., PERRON, J., \"A Computer Model of the Rotary Coke Calciner\", in Proc. Int. Symp. Comp. Software in Extractive Metall. CW. Bales and G.A. Irons, eds., CIM Proceedings, Montreal, 237-248 (1993). 42. KOCAEFE, Y.S., SIMARD, G., BUI, R.T., CHARETTE, A., POTOCNIK, V., PERRON, J., \"Analyzing the Heat Transfer in a Coke Calcining Kiln\", Light Metals 1992, Euel R. Cutshall, ed., The Minerals, Metals & Materials Society, 627-632 (1991). 182 43. BUI, R.T., SIMARD, G., KOCAEFE, Y., CHARETTE, A., LACROIX, M., JAIN, S., PERRON, J., PROULX, A. and BARR, P., \"3D Simulation of the Thermal Performance of a Coke Calcining Kiln\", in Proc. Int. Symp. on Extraction Refining and Fabrication of Light Metals, M. Sahoo and P. Pinfold, eds., CIM Proceedings, Montreal, 24:367-380 (1991). 44. SPALDING, D.B., \"The PHOENICS Reference Manual\", Report No. CHAM TR/200, CHAM Ltd., London, U.K. (1991). 45. HE, P., SALCUDEAN, M., GARTSHORE, I.S., NOWAK, P., \"Multigrid Calculation of the Fluid Flows in Complex 3D Geometries using Curvilinear Grids\", Computers and Fluids, 25(4):395-419 (1996). 46. LAUNDER, B.E., SPALDING, D.B., \"The Numerical Computation of Turbulent Flows\", Comp. Meth. App. Mech. Eng., 3:269-289 (1974). 47. MAGNUSSEN, B.F., HJERTAGER, B.H., \"On Mathematical Modeling of Turbulent Combustion with Special Emphasis on Soot Formation and Combustion\", 16th Int. Symp. on Comb., Cambridge, MA., 719-729 (1976). 48. CHANDRASEKHAR, S., \"Radiative Transfer\", Dover, New York (1960). 49. LATHROP, K.D., \"Use of Discrete Ordinates Methods for Solution of Photon Transport Problems\", Nucl. Sci. Eng., 24:381-388 (1966). 50. FIVELAND, W.A., \"Three-Dimensional Radiative Heat Transfer Solutions by the Discrete-Ordinates Method\", J. Thermophys. Heat Transfer, 2(4):309-316 (1988). 51. STULL, D.R., \"JANAF Thermochemical Tables\", The Joint Army-Navy-Air Force-ARPA-NASA Thermochemical Working Group, Aug. (1965). 52. SIEGEL, R., HOWELL, J.R., \"Thermal Radiation Heat Transfer\", McGraw-Hill, Inc. second edition (1981). 53. FIVELAND, W.A., JAMALUDDIN, A.S., \"Three-Dimensional Spectral Radiative Heat Transfer Solutions by the Discrete-Ordinates Method\", J. Thermophys. Heat Transfer, 5(3):335-339 (1991). 54. THEMELIS, N.J., \"Transport and Chemical Rate Phenomena\", Gordon and Breach Science Publishers, Switzerland (1995). 55. HSIAU, S.S., HUNT, M.L., \"Kinetic Theory Analysis of Flow-Induced Particle Diffusion and Thermal Conduction in Granular Material Flows\", J. Heat Transfer, 115:541-548 (1993). 183 56. SAVAGE, S.B., \"Disorder Diffusion and Structure Formation in Granular Flows\", in Disorder and Granular Media, D. Bideau and A. Hansen, eds., Elsevier Science Publishers, Amsterdam (1993). 57. CARNAHAN, N.F., STARLING, K.E., \"Equations of State for Nonattracting Rigid Spheres\", J. Chem. Phys., 51(2):635-636 (1969). 58. LUN, C.K.K., SAVAGE, S.B., \"The Effects of an Impact Velocity Dependent Coefficient of Restitution on Stressed Developed by Sheared Granular Materials\", Acta Mechanica, 63:15-44 (1986). 59. KELLOGG, H.H., \"Elements of Pyrometallurgy, Lecture Notes, Columbia University (1989). 60. FUJII, T., IMURA, H., \"Natural-Convection Heat Transfer from a Plate with Arbitrary Inclination\", Int. J. Heat Mass Transfer, 15:755-767 (1972). 61. ALYASER, A.H., \"Fluid Flow and Combustion in Rotary Kilns\", PhD Thesis, Faculty of Graduate Studies, Department of Metals and Materials Engineering, University of British Columbia (1998). 62. BARR, P.V., \"Heat Transfer Processes in Rotary Kilns\", PhD Thesis, Faculty of Graduate Studies, Department of Chemical Engineering, University of British Columbia (1986). 63. PATANKAR, S.V., \"Numerical Heat Transfer and Fluid Flow\", Hemisphere Publishing Corp. (1980). 64. BROOKS, A.N., HUGHES, T.J.R., \"Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations\", Comp. Meth. Appl. Mech. Eng.,32:199-259 (1982). 65. BARUZZI, G.S., HABASHI, W.G., GUEVREMONT, J.G., HAFEZ, M.M., \"A Second Order Finite Element Method for the Solution of the Transonic Euler and Navier-Stokes Equations\", Int. J. Num. Methods Fluids, 20:671-693 (1995). 66. HUGHES, T.J.R., \"A Simple Scheme for Developing 'Upwind' Finite Elements\", Int. J. Num. Methods Engin., 12:1359-1365 (1978). 67. CHRISTIE, I., GRIFFITHS, D.F., MITCHELL, A.R., ZIENKIEWICZ, O.C., \"Finite Element Methods for Second Order Differential Equations with Significant First Derivatives\", Int. J. Num. Methods Engin., 10:1389-1396 (1976). 68. HEINRICH, J.C., HUYAKORN, P.S., ZIENKIEWICZ, O.C., MITCHELL, A.R., \"An 'Upwind' Finite Element Scheme for Two-Dimensional Convective Transport Equation\", Int. J. Num. Methods Engin., 11:131-143 (1977). 184 69. ZHANG, Y., CAMPBELL, C.S., \"The Interface Between Fluid-Like and Solid-Like Behaviour in Two-Dimensional Granular Flows\", J. FluidMech., 237:541-568 (1992). 70. BOATENG, A.A., \"Rotary Kiln Transport Phenomena: Study of the Bed Motion and Heat Transfer, PhD Thesis, Faculty of Graduate Studies, University of British Columbia (1993). 71. BARR, P. V., Unpublished Research, University of British Columbia (2003). 72. VLIET, G.C., ROSS, D.C, \"Turbulent Natural Convection on Upward and Downward Facing Inclined Constant Heat Flux Surfaces\", J. Heat Transfer, 549-555, Nov. (1975). 73. HOLMAN, J.P., \"Heat Transfer\", McGraw-Hill, Inc., fifth edition (1981). 74. CHUCHILL, S.W., CHU, H.H.S., \"Correlating Equations for Laminar and Turbulent Free Convection from a Horizontal Cylinder\", Int. J. Heat Mass Transfer, 18:1049-1053 (1975). 75. WINSLOW, A.M., \"Numerical Solution of the Quasilinear Poisson Equation in a Nonuniform Triangle Mesh\", J. Comp. Phys., 2:149-172 (1967). 76. THOMPSON, J.F., THAMES, F.C, MASTIN, C.W., \"Automatic Numerical Generation of Body-Fitted Curvilinear Coordinate System for Field Containing Any Number of Arbitrary Two-Dimensional Bodies\", J. Comp. Phys., 15:299-319 (1974). 77. RAZINSKY, E., BRIGHTON, J.A., \"Confined Jet Mixing for Nonseparating Conditions\", J. Basic Engr., 333-349, Sept. (1971). 78. YULE, A.J., DAMOU, M., KOSTOPOULOS, D., \"Modeling Confined Jet Flow\", Int. J. Heat and Fluid Flow, 14(1):10-17 (1993). 79. RODI, W., \"The Prediction of Free Turbulent Boundary Layer by Use of Two-Equation Models of Turbulence\", PhD Thesis, Imperial College of Science and Technology, London (1972). 80. YULE, A.J., DAMOU, M., \"Investigations of Ducted Jets\", Exp. Thermal and Fluid Science, 4:469-490 (1991). 81. POPE, S.B., \"An Explanation of the Turbulent Round-Jet/Plane-Jet Anomaly\", AIAA Journal, 16(3):279-281 (1978). 82. BINDER, G., KIAN, K., \"Confined Jets in a Diverging Duct\", Proceedings, Turbulent Shear Flows 4, Karlsruhe, 7.18-7.23 (1983). 185 83. ZHU, J., SHIH, T.-H., \"A Numerical Study of Confined Turbulent Jets\", J. Fluids Eng., 116:702-706(1994). 84. ZIJLEMA, M., SEGAL, A., WESSELING, P., \"Finite Volume Computation of Incompressible Turbulent Flows in General Co-ordinates on Staggered Grids\", Int. J. Num. Methods Fluids, 20:621-640 (1995). 85. DEISSLER, R.G., \"Diffusion Approximation for Thermal Radiation in Gases with Jump Boundary Condition\", J. Heat Transfer, 86(2):240-246 (1964). 86. ROSSELAND, S., \"Theoretical Astrophysics; Atomic Theory and the Analysis of Stellar Atmospheres and Envelopes\", Clarendon Press, Oxford (1936). 87. LAUNDER, B.E., PRIDDIN, C.H., SHARMA, B.L, \"The Calculation of Turbulent Boundary Layers on Spinning Curved Surfaces\", J. Fluids Eng., 99:231-239 (1977). 88. HE, P., SALCUDEAN, M., GARTSHORE, I.S., \"A Numerical Simulation of Hydrocyclones\", Trans. IchemE, 77:429-441, Part A (1999). 89. MURTHY, M.S., HARISH, B.R., RAJANANDAM, K.S., AJOY PA VAN KUMAR, K.Y., \"Investigation on the Kinetics of Thermal Decomposition of Calcium Carbonate\", Chem. Eng. Sci., 49(13):2198-2204 (1994). 90. KAPLAN, W., \"Advanced Calculus\", Addison-Wesley Publishing Company, Inc. (1957). 186 APPENDIX A 20-Node, 3-D Isoparametric Element Description The 20-node element used in the discretization for the bed and wall models is shown in Figure A. 1. The quadratic description is used for both the unknowns and the geometry. (a) Element in Cartesian Coordinates (b) Element in Local Coordinates Figure A.l: 20-Node, 3-D Element Used for the Bed and Wall Models Shape functions Nt for nodes / = 1 to / = 20 are as follows, Corner nodes, / = 1-8 Ntf,V,C)=L-(\ + tfi)(\+m)(\ + \u00C2\u00A3 C M i + W, + \u00C2\u00AB*, -2) (A.l.a) On sides, / = 9, 12, 17, 19 N^n^^^l-^Xl + nnJl + CCi) (A.l.b) 187 On sides, / = 10, 11, 18, 20 Ntt,nX) = \(l-Tl2l\ + &)(l + tfi) (A.l.c) On sides, i = 13, 14, 15, 16 N,(\u00E2\u0082\u00AC,T,,C) = -\(l-C2Xl + &l)Ll + 'ml) (A.l.d) Description of unknown variables

,(\u00C2\u00A3\u00C2\u00BB7,\u00C2\u00A3>, (A.2.a) I=I 20 x^,V,C) = ZNM,TJ,Ck (A.2.b) 20 AZ,T,\u00C2\u00A3) = 1\u00C2\u00A3Nl{4,n,c:)yl (A.2.c) 20 i=l where ., x(, y., and z(. are values at nodes /'. 188 APPENDIX B Calculation of Surface Integrals Parametric equations x = x(^,n), y = y{^,n), and z = z(\u00C2\u00A3, 77) are taken as a mapping of an area in (^,TJ) coordinates onto a 3-D surface as shown in Figure B.l [90]. .U..j.... ................. (...\u00E2\u0080\u00A2-\u00E2\u0080\u00A2\u00E2\u0080\u00A2. i : : Figure B.l: Mapping of 3-D surface Tangent vectors Px (on 77= constant curve) and P2 on (\u00C2\u00A3 = constant curve) are described by, - dx ~ dy ~ dz r - dx ~ dy ~ dz r dn dn dn (B.l.a) (B.l.b) The cross-product P\~xP2 dictates the normal direction h to the surface and surface area A is evaluated through, ^ = JJjy3 x P2fedn = jjdA (B.2) 189 Care must be taken so that the unit normal vector h = away from the element surface. = -(5*5) PxxP2 , always points For a known heat flux q = q the heat rate, Q, through the surface is, Q=\\qhdA P{xP2 = j ^ x P ^ d n (B.3) If q is given by a convective flux (q = pCpTV ) the heat rate is evaluated by, = jjpCpT[v-(PlxP2)\i&ri (B.4) Finally, when q has the form of a diffusive flux (q = -KVT), Q = -jJKVT(PlxP2)dZdT1 dTd{y,z) + KdTd{z,x) , vBTB{x,y) + K-dx d{\u00C2\u00A3,r}) By B{\u00C2\u00A3,n) Bz B^n) (B .5) where etc. are the jacobians of transformation. For example, 1 0 0 Bx By Bz a7 a? al 3JC 3y 3z 3/7 877 9/7 (B.6) 190 APPENDIX C Granular Temperature Distribution in Bed Active Layer The bed surface granular temperature distribution is assumed to take the form, 2^=2,+2 2 ( 4- 4\u00C2\u00B0) where S, and E2 are the parabolic functions shown in figure 4.9. Besides the angular extent of the bed < 0 < 02) as shown in figure 4.10, three input variables are required for the parabolic description of equation 4.40. Namely, the two peak values of granular temperature, !Em and SOTj and the position of 6\ (of figure 4.9). S, and E2 are evaluated from. for d > e:, for 9 < e[, From the geometry, {e-djd-e2) ( C 1 ) E2 _ (e-0\\d-e2) En,2 (6mi -0\Xen,2 ~02) (C.2.a) :2=0 (C.2.b) em=6-^ ( c . 3 ) em=e;+^A (c.4) 9t is evaluated from the \"fractional range for s 2 \" orr&8 as follows, 191 (C.5) With a calculated free surface distribution, the depth distribution for E (in active layer region) is assumed to have an exponential form as follows, urface e A -1 e\u00C2\u00B0GT -1 (C.6) where A is the local thickness of the active layer at the angular position 6 as shown in figure 4.10. The granular temperature distribution in the active layer is a function of the radial position 8 and a constant -79632.231 kJ/Kgmol CH4 species C02 02 N2 H20 moles heat form. h at T h at 298 sensi. heat 1 1.4436 12.947936 2 -393520 0 0 -241820 34373.81 25900.39 25010.13 29612.98 9354.61 8678.468 8619.556 9898.194 25019.2 17221.92 16390.58 19714.78 kJ/Kgmol CH4 -368500.8 24861.568 212224.13 -444210.43 Sum of Products => -575625.54 kJ/Kgmol CH4 Net Heat -495993.31 kJ/Kgmol CH4 Net Heat -30999.582 kJ/Kg CH4 Heat Rate -54 484865 kW Heat Absorbed in Bed Reactants at T_b,in species CaC03 Products at T_b,out species C02 CaO CaC03 300.0 moles 1 1089.0 moles 0.905 0905 0.095 heat form. -1207000 heat form. h at T h at 298 sensi. heat kJ/Kgmol CaC03 34500 34270 230 -1206770 Sum of Reactants => -1206770 kJ/Kgmol CaC03 (C02 at TJrf.out) 835.0 h at T h at 298 sensi. heat -393520 34373.81 -633559 \u00E2\u0080\u00A21207000 70131.6 125235 9354.61 19191.2 34270 25019.2 50940.4 90965 Sum of Products => kJ/Kgmol CaC03 -333493 23 -527269.83 -106023.33 -966786.38 kJ/Kgmol CaC03 Net Heat 239983.615 kJ/Kgmol CaC03 Net Heat 2399.83615 kJ/Kg CaC03 Heat Rate 36.6641634 kW Overall Balance Heat Rate Loss (other than wall) 2.00 Heat Released in H.F. -54.484865 Heat Absorbed in Bed 36.6641634 HEAT RATE LOST = HEAT RATE SOURCE BY REACTIONS ==> Q w,l 15.82 Burner Thermal Load -802310 kJ/Kgmole CH4 -50144.375 kJ/KgCH4 Bed Thermal Load and Q b,r Q_hf,r 88.13375 (complete combustion assumed) Bed Thermal Load (max) 179921 kJ/Kgmole CaC03 1799.21 kJ/Kg CaC03 <== Enthalpy of reaction used in Bed Model 27.487931 kW Q b,r 24.876577 kW (equals to Bed Thermal Load if f=1) Hot Flow Input Heat Rate Reactants at T_hf,in 288.7 species moles h at T kJ/Kgmol CH4 CH4 1 9703.901 9703.901 02 3.4436 8402.684 28935.48 N2 12.947936 8348.818 108100 Sum of Reactants => 146739.3 kJ/Kgmol CH4 9171.209 kJ/KgCH4 Q_hf,in=> 16.11932 kW Hot Flow Output Heat Rate Products at T_hf,out 835.0 species moles h at T kJ/Kgmol CH4 C02 1 34373.81 34373.81 02 1.4436 25900.39 37389.8 N2 12.947936 25010.13 323829.6 H20 2 29612.98 59225.95 Sum of Products => 454819.2 kJ/Kgmol CH4 28426.2 kJ/Kg CH4 Q_hf,out =: 49.96188 kW Bed Input Heat Rate Reactants at T_b,in 300.0 species moles CaC03 1 Bed Output Heat Rate Products at T_b,out 1089.0 species moles CaO 0.905 CaC03 0.095 hatT 34500 Sum of Reactants => kJ/Kgmol CaC03 34500 34500 kJ/Kgmol CaC03 345 kJ/Kg CaC03 Q_b,in => 5.270833 kW h at T kJ/Kgmol CaC03 70131.6 63469.1 125235 11897.33 Sum of Products => 75366.42 kJ/Kgmol CaC03 753.6642 kJ/Kg CaC03 Q_b,out=> 11.51431 kW Heat Rate of C02 loss from Bed Products at T_top,C02 1089.0 species moles C02 0.905 Hot Flow Output Heat Rate of Bed C02 Products at T_hf,out 835.0 species moles C02 0.905 kJ/Kgmol CaC03 43195.93 hatT 47730.31 Sum of C02 Prod. => Q_top,C02 =====> hatT 34373.81 Sum of C02 Prod. => 43195.93 kJ/Kgmol CaC03 431.9593 kJ/Kg CaC03 31108.3 31108.3 kJ/Kgmol CaC03 311.083 kJ/KgCaC03 6.599378 kW kJ/Kgmol CaC03 Q_b.C02 =====> 4.752656 kW "@en . "Thesis/Dissertation"@en . "2004-05"@en . "10.14288/1.0103853"@en . "eng"@en . "Mechanical Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en . "Graduate"@en . "Mathematical modelling of lime kilns"@en . "Text"@en . "http://hdl.handle.net/2429/15899"@en .