M A T H E M A T I C A L M O D E L L I N G O F 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 © 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 threedimensional 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 released to the hot flow. ii CaC03 and CaO. Carbon dioxide is 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.3 3.2.4 Bed Reaction Rate Model 41 3.2.5 44 Summary of Governing Equations for Bed Model 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 56 Bed Flow Field 4.2.1.1 Boundary Conditions 57 4.2.1.2 Discretization 58 4.2.1.3 Three Dimensional Velocity Field 4.2.2 62 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 4.3 5.0 4.2.2.4 Granular Temperature Distribution 75 4.2.2.5 Calcination Temperature Distribution 79 Wall/Refractories Model 81 4.3.1 83 4.3.2 4.4 74 Boundary Conditions 4.3.1.1 Free Convection Heat Transfer Coefficient 85 4.3.1.2 Total Heat Loss to Surroundings 89 Discretization 90 Coupling and Information Exchange 93 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 6.1 6.2 7.0 8.0 Overall Structure of Model and Convergence of Hot Flow/Wall Coupling 111 Overall Structure of Model and Convergence of Fully Coupled Run 115 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 — £ Model 128 7.4 Summary Remarks on Modelling of Lime Kilns 137 RESULTS AND DISCUSSION 141 8.1 Hot Flow/Wall Coupled Results 141 8.1.1 148 8.2 8.3 9.0 111 A Transient Solution of the Hot Flow 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 Summary 174 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS BIBLIOGRAPHY 175 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 1.4 Fire Hood (TakenfromKramm, D.J. [2]) Recirculation Patterns in the UBC Pilot Kiln (scale distorted) (A-Flow separation point; B-Flow reattachment point) 4 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 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 ofC0 to Hot Flow Model 4.8 Experimentally Observed Granular Temperature Distribution 2 at the Bed Surface 8 49 75 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 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 7.3 Measurement for Spreading of Jet Lj and Duct Boundary Layer L B L Spreading of Jet and Duct Boundary Layer (Exp. from Razinsky et al. [77]) 7.4 Center Line Velocity Decay (Exp. from Razinsky et al. [77]) 7.5 Center Line Velocity Decay with Pope's Correction [81 ] (Exp. from Yule et al. [80]) Center Line Velocity Decay with Pope's Correction [81 ] (Exp. from Razinsky et al. [77]) 7.6 ix 103 122 123 123 126 127 7.7 Mean Velocity Distribution at Axial Position of 1.875 Inlet Diameters for the Conical Duct of Binder et al. [82] 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 datafromBinder 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 datafromBinder 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) Temperature (K) in Central Vertical Plane of Generic Industrial Kiln With and Without Buoyancy (Axes are in meters) Temperature (K) in Cross Sectional Plane at Axial Position of 6.6 m for Generic Industrial Kiln (Axes are in meters) 7.8 7.15 7.16 129 137 138 139 8.1 Fuel (CH ) 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) 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 8.8 4 x 149 8.13 Inner Wall Temperatures (K) in the Fire Hood Region for Fully Coupled Results 8.14 Bed Surface Temperature Distribution (K) 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 161 162 163 8.16 Solids Volume Fraction of CaC0 in Bed 8.17 Cross-Section Solids Volume Fraction of CaCC»3 at Axial Positions of 3 (a) 0.0 m(b) 0.5 m and (c) 1.3 m 164 165 8.18 Surface Contour for the Reaction Rate at a Value of 0.0005 per sec 8.19 Reaction Rate Contours (per second) at Axial Positions of (a) 0.5 m and (b) 1.5 m 166 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 Inner Wall Temperature as a Function of 172 Circumferential ( 6 ) Direction 173 8.24 A. 1 20-Node, 3-D Element Used for the Bed and Wall Models B. l Mapping of 3-D surface 188 Cl 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 xi 186 List o f Tables 3.1 3.2 Components of Equation 3.1 Constants in the turbulence model [46] 4.1 Constants Used in Bed Model (for case T21 of Barr's experiments [62]) 4.2 Reference Variables for the Non-Dimensional Bed Equations 4.3 Fluxes Evaluated in the Post-Processing Phase (Dimensional Form) 4.4 Constants in Wall Model 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 Computing Times for Fully Coupled Run (Based on a 3 GHz Pentium IV Processor) 116 6.2 28 29 65 67 74 82 8.1 Pilot Kiln Trial Run #6 from Alyaser [61] 8.2 Model Heat Transfer Rate Imbalances for Hot FlowAVall Coupling 8.3 Comparison of Heat Transfer Rates to Experimental Data of Alyaser [61] 142 8.4 Pilot Kiln Trial Run T21 from Barr [62] 8.5 Model Heat Rate Imbalances for Fully Coupled Results 8.6 Comparison of Heat Rates to Overall Heat Balance Model of Appendix G 155 xii 142 142 154 154 Nomenclature OR Rosseland mean absorption coefficient (1/m) absorption coefficient (1/m) A area (m) c, constant in K - e turbulence model C 2 2 Q constant in K - e turbulence model constant in Pope's modification to K-e model specific heat (J/kgK) constant in K-e turbulence model diameter of particles (m) d artificial diffusion coefficient (m Is) granular diffusion coefficient (m/s) 2 D 8 D ,D„ inner kiln diameter (m) E P coefficient of restitution (-) E energy (J) f solids mass fraction of CaCCb,fractionalfill for bed (-) g solids massfractionof CaO (-) g gravitational constant vector (m/s ) m g (v) radial distribution function (-) G generation term (equation 3.1) (kg/ms) h enthalpy (J/kg) a 3 effective heat transfer coefficient (W/m K) h spectral radiative intensity (W/mrad) 3 black body radiative intensity (W/m rad) K conductivity (W/mK) xiii K a artificial conductivity (W/mK) K g granular thermal conductivity (W/mK) K wall conductivity (W/mK) K' bulk bed conductivity (W/mK) /, turbulent length scale (m) L kiln length (m) w m fraction of original mass of CaC03 converted to CaO (-) w fraction of original mass of CaCC»3 converted to C O 2 2 3 (-) m entrain entrainment rate (kg/s) m, massfractionof species i (-) m p primary jet massflowrate (kg/s) m s secondary mass flow rate (kg/s) MW; molecular weight of species i (kg/kgmole) n unit normal vector (-) P pressure (N/m) q excess volumetric flow rate of primary discharge (m/s) q heat loss to surroundings byfreeconvection (W/m) 2 3 2 fc q heat loss to surroundings by radiation (W/m) q radiative heat flux vector (W/m ) Q total volumetric flow rate (m/s) Q rate of heat lost with C0 (W) 2 rad r 3 COI QCOND 2 r a t e °f heat transfer by granular and bulk diffusion (W) Q rate of heat to a particle (W) Q R radiative source term for equation 3.1 (W/m) Q R rate of heat absorbed by reaction (W) r radius at particle reaction interface (m) P t 3 xiv primary jet radius (m) p r kiln jet radius, limestone particle radius (m) R reaction rate (1/s) IN inner kiln radius (m) R' excess discharge ratio R universal gas constant (J/(kgmole)(K)) R radius vector (m) s ray direction (m) R (q/Q)(-) non-dimensional rate of strain tensor (-) S S source terms (equation 3.1) T temperature (K) T calcination temperature (K) c film temperature (K) f T T outer wall shell temperature (K) T ambient temperature (K) sh oo U velocity in coordinate direction x (m/s) P average primary velocity (m/s) s average secondary velocity (m/s) U U u,u m 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: a BED angle of repose (deg.) diffusion coefficients (equation 3.1) S' boundary layer displacement thickness (m) £ turbulent energy dissipation rate (m /s ) 2 e emissivity (-) 6 angular position (deg, rad) K turbulent kinetic energy (m /s ) X wavelength (m) (1 molecular dynamic viscosity (kg/ms) ju effective viscosity (//,+//) (kg/ms) fl, turbulent viscosity (kg/ms) v solids fraction (-) £ unknowns (equation 3.1) E granular temperature (m /s ) p density (kg/m) Pcaco, effective density of CaC03 (kg/m) p Ca0 effective density of CaO (kg/m) p w wall density (kg/m) a Stefan-Boltzmann constant (W/m K ) s 2 e 2 3 2 2 3 3 3 3 2 4 <J GT granular temperature distribution factor (-) <J T calcination temperature distribution factor (-) o~ constant in K - £ turbulence model <T constant in K - £ turbulence model (p solids volume fraction of CaC03 (-) X non-dimensional measure of vortex stretching (-) *F stream function (m/s) K £ 2 xvi wall rotation rate (1/s) CO vorticity component in x-y plane (1/s) co z co non-dimensional rotation tensor (-) Q. rotation rate vector (rad/s) ij Non-dimensional numbers: artificial Bodenstein number granular Bodenstein number ° B s Curtet number (l/Vw) Ct Da Damkohler's First Da, Damkohler's Third Gr Grashof number GK modified Grashof number m Craya-Curtet parameter Nu Nusselt number Pe artificial Peclet number * granular Peclet number x a P g (GrNu) Pe 'rotational' Peclet number for wall Pe bulk bed Peclet number Pr Pradtl number Pr, turbulent Pradtl number (constant in K - £ model) Re Reynolds number w Sm turbulent Schmidt number (constant in K - £ model) St Stefan number Th Thring-Newby parameter t xvii Acknowledgements I consider myself very fortunate to have had the opportunity to work with some of my disciplines greats. First and foremost my supervisors, Dr. Martha Salcudean and Dr. Ian Gartshore, who would call me Michael when they were pleased with my progress and Mike when they were less than enthusiastic. The influence they have had in my life will always be valued. My sincere appreciation for Dr. Philip Hill who was there for me, especially in the early stages of the research. I would like to thank Dr. Peter Barr who helped me decipher through the experimental data used in this thesis for validation purposes and for his valuable contributions for the development of the bed reaction rate model. I would also like to express my gratitude to Dr. Paul Nowak, the developer of the hot flow model used in this thesis, who always had time for discussions no matter how busy he was. Finally, I would like to thank my family for their loving and continued support. xvui 1 1.0 I N T R O D U C T I O N In a Kraft process of pulp and paper production, white liquor containing active chemicals sodium hydroxide, N a O H , and sodium sulfide, Na2S, are recovered in the sodium loop for re-use, as shown in Figure 1.1, due to environmental and economic reasons. The wood pulp fibers go through a bleaching process. The pulp then flows through the headbox and eventually dried to the final product. Figure 1.1: Kraft Process 2 In the sodium loop, black liquor is concentrated and burned in the recovery boiler. The resulting smelt is dissolved and reacted with quick lime, CaO, to convert sodium carbonate, Na2CC>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 farfieldoffig.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 varyfrom1 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) dropsfromthe cylinder to the chute below at the bottom of thefirehood 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 thefirehood 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 o f these furnaces. Efforts to date in the prediction of the turbulent diffusion flame in the rotary k i l n have ranged from an estimation of the flame length alone [3] to detailed threedimensional flow studies [4]. A n empirical flame length model for example was developed by Ruhland [3] for flames in the cement rotary kiln. 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 kiln was confirmed by measurements o f the flame in a full scale rotary kiln. The change in volume of substances during combustion was accounted for by using a volumetric flow ratio between entrance conditions and endof-flame conditions (using adiabatic flame temperature information). Although useful, the correlation predicts the length o f flames in kilns for a nonswirling jet o f 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 liquidflameand the rotary kilnflame.As described below, this makes Ruhland's correlation suspect. The primary coreflow(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 afreejet. 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 thisflowfield(code used, boundary conditions, geometry, etc.) will be given later in this thesis. The axial velocityfieldreveals recirculation in the burner, thefirehood 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 coreflow(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> a> b 0 Core Air Entrainment -0.4 1 -0.4 •' • i • * • • * • • •' i • • • • i • •' • i ' > •' i •''' i ' ' •' i ' • •' i'' • • i • • -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 (>10) in rotary kilns to be Reynolds number independent [6]. Experiments such 5 as Ruhland's could have been carried out at lower Reynolds number without a great loss of similarity provided they were kept above 10. 4 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. 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 (0.2 p 1) (1.1) p r and the universal jet spreading angle, x =A.5r s (1.2) they proposed a similarity parameter based on the separation to reattachment distance ratio. Here r and m are the radius and mass flow rate of the primary jet. r is the kiln p p s 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 secondaryflowexactly satisfies the entrainment requirement (m = w s en/ra m) 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,, + m s V m P J r p r _ Separation..distance Re attachment, .dis tan ce (13) 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 r jr > 0.05. Under such conditions there are excessiveflowdistortions in the nozzle p s 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 (VO - (1.4a) 2 where R =q/Q, (1.4b) q=nr (u -u ) (1.4c) 2 p p s and Q = 7t(r, -S') u 2 s +q (1.4d) Here u and u are the average velocities for the primary and secondary flows. The p s 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<Ct<0.8). They used cold flow (air and water) modelling techniques with operating conditions to cover the above range and have shown that similarity parameters such as Ct are useful indicators of recirculation levels. 11 The apparent need to work in this range emphasizes the important role recirculation plays in the rotary kiln. If Ct is too high (little or no recirculation) there may not be enough mixing for complete combustion. This leads to excessive emissions, problems with fuel efficiency etc. If Ct is too low (very large recirculation region) the recirculating gases lower the oxygen concentration around the flame and may result in precisely the same conditions [6]. At very low Ct values the flow patterns within the burning zone are also affected significantly through oscillations in the core jet flow [6]. Numerical simulations over a range of Ct values and with constant density have been made as part of the present thesis and will be described in Chapter 7. 1.2 Heat Transfer By the 1970's efforts were underway in the modelling of heat transfer for rotary kilns. Heat transfer paths in a rotary kiln are illustrated in Figure 1.5. The flame/hot flow region transmits energy by radiation and convection to the exposed wall and bed surfaces. Radiation heat exchange also occurs between these exposed surfaces. (In most industrial rotary kilns radiative heat transfer is significant due to the high hot flow temperatures). Heat regeneration to the bed via the rotating wall is an integral part of the process in the rotary kiln. Within the bed, heat is both conducted and convected due to the motion of the particles. Relative motion between particles (granular diffusion) adds to this heat transfer. 12 A certain degree of radiation is also present within the bed. Finally heat is lost to the ambient at the outside of the kiln. Figure 1.5: Heat Transfer Paths Many assumptions (including the omission of selective heat transfer paths) were made in the early efforts of heat transfer modelling. One-dimensional modelling, where the dimension of interest is the axial direction, typically ignores combustion aerodynamics by prescribing the flame length. In addition, the length to diameter ratio in rotary kilns is large enough to consider only local heat transfer exchanges in such models. It has been shown for example that radiative exchange is limited to a finite distance upstream and downstream from a given cross-section [12]. The geometry of the kiln as mentioned above, along with gas absorption effects makes this assumption possible. The magnitude of axial heat transfer in the rotary kiln (including wall and bed) has been assumed to be negligible relative to radial transfer in these models. For these reasons dividing the kiln into isothermal slices became the standard in heat transfer modelling [13]. A cross-sectional model using local conditions 13 such as temperatures and gas concentrations can then be applied repeatedly over the length of the kiln. One of the earliest studies was that of Pearce [14] who considered heat transfer between the gas and bed surface (by radiation and convection) and heat transfer within the bed (by conduction and radiation). Pearce pointed out that improving such models requires validation of heat transfer coefficients through experiments. Sunavala [15] went as far as to use a constant overall heat transfer coefficient of 350 W/m K and applied it to the theory of a counter-flow heat exchanger. Simplifications 2 such as uniform mass flow rates along the length of the kiln for both the bed and hot flow, allowed expressions to be derived for estimating the length of the different heat transfer zones in a cement kiln. Cross and Young [16] developed a model of a rotary kiln used in the processing of iron ore pellets. A gamma distribution for the heat released by the flame was used along with constant cross-sectional gas and bed temperatures. Although bed reactions were ignored, Cross and Young were thefirstto include regenerative heat transfer via the rotating wall. Rates of heat flow from the gas to the wall and the bed have been analyzed by Watkinson and Brimacombe [17] from temperature measurements in the UBC pilot rotary kiln [18]. They found that the gas to bed heat flux is limited by mixing in the solids at low feed rates (and corresponding rotational speeds). This was attributed to the slumping action of the bed where limited particle mixing is present at the surface as the kiln rotates, until a critical level is reached where particles suddenly slump down to a lower level. At higher feed rates, a rolling motion ensures proper mixing at the surface 14 and heat transfer is limited by the gas side heat transfer. As a result, higher heat fluxes were observed. After showing that radiative exchange is limited in the domain of a rotary kiln, Gorog et al. [12] developed a local radiation model (for a given transverse cross-section) that considers only a short distance upstream and downstream from a given axial position. Barr et al. [19] carried such radiative heat transfer calculations approximately 3 kiln inside diameters in each axial direction. In their work, Barr et al. developed a unified model for the heat transfer and successfully compared results with detailed measurements on the UBC pilot kiln [13] in which 23 experimental trials were made using limestone, Ottawa sand and petroleum coke as the feed. The model was used to explain how suppressed bed temperatures (relative to the inside wall) in both the endothermic calcination zone and the initial bed inlet zone produced an experimentally observed increase in the net heat transfer to the bed. Gorog et al. [20] in their series of papers also looked at regenerative heat transfer by modelling the rotating wall. They showed that regeneration is substantial in the cold end of a rotary kiln. The wall model was simplified by considering two regions, the first being a thin inner wall 'active layer', where a cyclic temperature change occurs as the wall rotates. In the remaining wall region no cyclical temperature variation takes place. Conduction in both regions was considered only in the radial direction. Heat transfer at the part of the wall surface covered by bed material was calculated by considering an average heat transfer coefficient and assumed that the contacting bed material remained at a constant temperature. 15 Solving the same problem for the wall, Barr et al. [19] avoided the problem of specifying a boundary condition at the interface between the bed and the wall by extending the conduction model into the plug flow region of the bed. The boundary condition then was simply the bulk temperature of the bed. Heat transfer in the flame zone was also considered by Gorog et al. [21] by assuming the flame to be cylindrical in shape. The length of the flame and secondary air entrainment were once again characterized by empirical equations. Three types of fuel were considered and various general conclusions were made on the effects of secondary air temperatures, increases in primary air and oxygen enrichment. Many more one-dimensional models have appeared in the literature for rotary kilns [22-27]. A series of equations (typically 14 ordinary differential equations) representing conservation of mass, energy and species averaged over the cross-section were solved using appropriate numerical methods. Correlations and algebraic equations were always used in these 1-D modelling efforts. 1.3 Bed Transport Processes Granular flow in the rotary kiln has numerous detailed heat transfer paths. Besides conduction of heat within a particle there is particle-to-particle conduction through contact, particle-to-particle radiation and gas-to-particle convection through voids. Details of granular flow behaviour however are beyond the scope of this thesis. 16 Six modes of bed motion have been identified [28] including slipping, slumping, rolling, cascading and centrifuging but the rolling mode is most common in rotary kiln operations. Rutgers [29] identified particle motion in a rolling bed. In the cross-section a plug flow region [30] exists where the particles rotate as a rigid body with the same angular rate as the kiln wall. Particles within this region have no axial velocity. This flow feeds the bed active layer at the surface. Particle mixing is high in this thin layer, which is typically about 10% of the total bed depth [28]. Mixing primarily occurs in the crosssectional plane. Macauly and Donald [31, 32] found that mixing in this plane is at least two orders of magnitude higher than in the axial direction. The high mixing in the bed also promotes segregation, where smaller particles tend to drop through the matrix of larger particles. As a result fine material tends to concentrate within a central core. This effect is not considered in this thesis. The motion of the bed particles in the cross-sectional plane plays an important role in bed heat transfer. When mixing is high, energy will be distributed thoroughly within the bed. The effective conductivity of the bed will increase. In the limit the bed will be isothermal at any given cross-section. One other important heat transfer mechanism is the high heat flux to the bed (relative to the wall) from the hot flow due to motion of particles at the bed surface. Tscheng and Watkinson [33] measured the convective heat transfer coefficient to the bed (using hot air gas temperatures of 350-590 K) and found it to be up to ten times the value of the coefficient to the exposed wall. Barr et al. [19] derived the convective heat transfer coefficients from the UBC pilot kiln trials by subtracting the calculated radiative 17 contribution from the net heat fluxes and found convection to the bed surface to be about two times that to the exposed wall. Watkinson and Brimacombe [17] postulated this to be associated with motion of particles on the surface of the bed. A rolling bed, as indicated earlier, ensures high mixing and higher heat transfer from the hot flow. Also, a larger effective surface area of the bed is exposed to the gas than would be the case for a simple plane. (Close packed spheres in a planar array for the bed for example would enhance the area by about 80%.) This alone has the potential of increasing both radiative and convective heat flux to the bed. 1.4 Multi-Dimensional Modelling Modelling has aided the design and operation of rotary kilns over the years. As discussed in Section 1.2 many one-dimensional models have appeared in the literature. These models can be efficient tools for the optimization and operation of the kiln by studying the effects of key control variables such as kiln slope, speed and production rate. They have in fact been the basis for the design and optimization of rotary kilns. The effort however has not included a detailed understanding of the phenomena associated with kiln operation but is rather based on the experience of the designer or operator. Such models have been generated for the overall kiln but they are not always suitable for detailed study. They do not, for example provide the tools to optimize the combustion process and burner design. 18 The critical assumption that must be made in any one-dimensional model is that uniform conditions exist across the cross-section in the freeboard gas, the walls/refractories and the bed. The bed for example is assumed to be well mixed and isothermal in any given cross-section. Flame positioning cannot be predicted by a onedimensional model. Although these models have been successfully used in industry, they provide a limited amount of information. Jenkins and Moles [34] were one of the first to develop an axisymmetric model for the hot flow in a cement kiln. They used the zone method for radiation of Hottel and Sarofim [35]. The flow pattern however was specified using the entrainment rate of Ricou and Spalding [5]. Measured concentration profiles from a cement kiln were also used as input. Generation rates of species such as CO2 were used to evaluate heat release in the flame. Resulting axisymmetric temperature contours compared favorably with measured data. The rotating wall and bed were not included in this study. Evidence (such as a non-uniform product) has suggested that large temperature gradients exist near and within the bed. As a result, a number of researchers have begun the quest for a more encompassing modelling effort. Boateng and Barr [36] have coupled a conventional one-dimensional plug flow model with a two-dimensional representation of the bed's cross-section [37]. This improves the simulation of conditions within the bed. The simulation looks at the non-uniformities within the bed without accounting for the complex phenomena occurring in the hot gas. Alyaser [38] has modelled the hot gas flow in the absence of a bed for axisymmetric conditions. His model was validated using thermal measurements from UBC's pilot kiln. This effort demonstrates how a model may be used to capture flame 19 phenomena for rotary kilns (adjustment of primary air ratio, primary jet momentum, extracting flame length etc.). Mastorakos et al. [39] have also coupled a 2-D axisymmetric solution of the gaseous phase with effectively one-dimensional models for the rotating wall and bed. This work was carried out for a cement kiln. Included were the bed clinker chemistry and an additional 'attached clinker' conductive layer on the wall. The radiation intensity in the gas was calculated with a Monte-Carlo method. (This was included as a radiation module in the CFD code FLOW-3D [40]). A constant absorption factor was used for the gas in the whole domain. Effects of local gas composition and temperature on the absorption factor were ignored. Bed temperature and composition were solved by considering only axial gradients and neglecting conduction. Despite all the uncertainties involved in the models the authors claim that quantities such as clinker exit composition, shell temperature and exhaust gas composition were predicted to acceptable accuracy. The development of CFD and the need to model the coupled effects of various rotary kiln hardware (as well as the coupled effects of associated physics such as aerodynamics, heat transfer, etc.) has led to bolder modelling efforts for related furnaces. Bui et al. [4, 41-43] were the first to construct a three-dimensional model for a rotary petroleum coke calciner. A coupled simulation of the physical phenomena that occur in the kiln was undertaken using the CFD code PHOENICS [44]. Four submodels (gas, bed, refractory and radiation) were coupled using four grid domains. The model takes into account the major phenomena of interest including the gas flow, all modes of heat transfer, combustion, bed motion and thermal effects of the refractories. 20 In a coke calcining kiln, volatiles evolve from the bed (hydrogen, methane and tar) at about mid-length. Evolution of these volatiles from the bed results in subsequent energy self-sufficiency for the process. In their model, Bui et al. assume that the volatiles burn in the hot flow only. Unlike the lime kiln, heat from the burner is primarily required for the start-up phase of the operation. Combustion from the burner is not included in the work of Bui et al. Primary and secondary air is provided in the burner end of the kiln. However, 'tertiary' air is also injected into the calcining zone (mid-length) through nozzles fixed on the rotating wall. These nozzles plough through the bed at each rotation and are modelled by Bui et al. using a pseudo steady-state approach. Although in their initial work the combustion pattern was imposed empirically, Bui et al. have recently published results [4] which include calculations of species transport and combustion in the hot flow. A three-dimensional calculation is made for the granular bed assuming an active and a plug flow layer. Each layer is assumed to behave as a Newtonian fluid. Mass, momentum and energy equations are solved in the bed. Three viscosity values (one for each of the active layer, plug flow and axial direction) are used in order to match experimentally observed velocities. An assumption is made that the bed's cross-sectional velocity profile does not change along the kiln length. Thermal conductivity of the bed is assumed constant throughout at 0.1 W/mK. Mass evolution of the three volatiles from the bed to the hot flow is imposed from experimental data in the form of temperature dependent bell-shaped curves. No reactions are explicitly calculated in the bed model. 21 Results such as cross-sectional average coke bed temperatures versus axial distance were compared by Bui et al. to an industrial kiln. Computed temperature profiles along the length of the kiln were within 2 0 0 C of experimental data. The combustion zone, tertiary air use and bed chemistry makes the coke calcining kiln considerably different from the rotary lime kiln considered in the current research. The work presented in this thesis, in part, concentrates on the lime kiln calcination process within the bed and the resulting C0 transfer to the hot flow. Fire hood/kiln 2 interaction capability is illustrated and flame zone details including flame instability issues are investigated. A detailed tool results for the design and optimization of the lime kiln operation. 1.5 Summary There are numerous economical and environmental reasons for the development of a coupled detailed model of fluid flow, combustion, heat transfer and other transport phenomena associated with the rotary lime kiln. Kiln fuel costs for example are typically a major production cost in a mill. The kiln is the major user of purchased fuel. Combustion efficiency must therefore be maintained at a maximum. Pollutant emissions from the combustion in lime kilns are of concern for environmental protection. As mentioned earlier, a badly operating kiln can destroy parts of the kiln itself such as the refractory lining. In an effort to examine many of these issues in a unified manner a model must account for the interacting hardware (kiln, burner, fire hood, etc.) in lime 22 kilns. In addition the model must represent the coupled effects of aerodynamics, heat transfer and the non-uniform cross-sectional conditions. N o such model currently exists for the lime kiln. In the present study, fully coupled three-dimensional modelling is applied to the rotary lime kiln. Care is taken when modelling the geometry: The capability of capturing details of the burner and hood/kiln interactions is illustrated. Buoyancy effects are included in the calculations, enhancing the predictive capability considerably. The granular material is modelled through a 3-D energy and species transport in the bed including the calcium carbonate chemistry specific to lime kilns. The conductive rotating wall is also modelled in 3 - D . The resulting model can provide insights into the physical processes taking place in the rotary lime k i l n and can be used as a validation tool for less complete models. It can also be used to optimize kiln design and operations. 23 2.0 SCOPE AND OBJECTIVES The need for a unified rotary lime kiln model has been described in Chapter 1. The overall objective of this thesis is to model the rotary lime kiln in order to provide a tool for operation and design optimization. As a result, the prime focus of the current research is to couple the interactions of the hardware (kiln, burner, fire hood, bed, wall) along with all the physical phenomena present (aerodynamics, heat transfer, fluid flow, radiation, combustion, bed motion, bed reaction, wall rotation, losses to ambient) in the preheat and calcination zones. All models are solved in 3-D. Resulting coupled solutions may be used to identify problems of kiln operation or design. Objectives of the current research are to, 1) Consider the accuracy of the existing model for cold and hot flow conditions using experimental data where available and investigate the aerodynamics in a rotary lime kiln through, - core flow spread rates - center line velocity decay - cross sectional velocity profiles - turbulent intensities and stresses, and - recirculation patterns for a range of Ct numbers. 2) Investigate the accuracy of the standard K-e model and consider modifications for improvement. 3) Formulate a hot flow solution using an existing CFD solver that captures details of kiln, burner and fire hood interactions. 24 4) Develop a bed model that includes bed motion, species and thermal transport as well as the calcium carbonate chemistry specific to lime kilns. 5) Develop a wall model including rotation and losses to the ambient. 6) Couple the hot flow, bed and wall models and critically investigate the physical processes taking place in the rotary lime kiln. 7) Evaluate the effect of flame buoyancy and the possible impingement on the refractories. 8) Investigate the transient response of a flame at low Ct numbers. 9) Use available measurements from the UBC pilot kiln to validate the model. Development and structure of the model (objectives 3 to 5) is addressed in chapters 3 to 6. Objectives 1 and 2 (accuracy and aerodynamics) are covered in chapter 7. Finally, the coupled model results and discussions (objectives 6 through 9) are addressed in chapter 8. Conclusions and recommendations for future work are presented in chapter 25 3.0 PHYSICAL MODEL AND GOVERNING EQUATIONS A unified model for the lime kiln must include modelling o f the hot flow (with kiln, fire hood and burner interactions), the bed transport processes (with calcination chemistry) and the rotating wall (with losses to the ambient). Descriptions of each o f these models and the assumptions made in creating them are outlined i n this chapter. The three models (hot flow, bed and wall) are coupled by information exchange between them (temperature, heat flux and flow rate of CO2) as shown i n Figure 3.1. The hot flow model calculation uses temperatures from the bed and wall as boundary conditions. The wall model calculation uses heat fluxes as boundary conditions from the bed and hot flow. This coupling w i l l be described further in Chapter 4. Equations solved in each domain are referred to i n the figure through equation numbers given i n this chapter. HOT FLOW MODEL: (continuity, momentum, turbulent kinetic energy, turbulent dissipation rate, species,energy) BED MODEL: Equations 3.31, 3.19,330 (continuity, species, energy) WALL MODEL: Equation 3.44 (energy) Figure 3.1: Coupling of Models 26 Only the preheat and calcination zones are included in the present research and all the important phenomena in these zones are modelled. The chain zone with its evaporation of water is not considered in this work. The UBC pilot kiln runs modelled here [13] had dry feed (zero moisture content). 3.1 Hot Flow Model The CFD solver used in this study to obtain hot flow solutions was developed by Nowak within the Department of Mechanical Engineering at the University of British Columbia [45]. The approach is based on the finite volume method. The model includes buoyancy effects, turbulence (using the standard K-e model [46]), evolution and combustion of gaseous species (through Magnussen's model [47]), and radiation through the ray-tracing technique (Discrete Ordinate Method [48-50]). Combustion of natural gas or oil may be modelled. Two major assumptions are made in the hot flow model. First, although a transient solution of the hot flow alone will be given in Chapter 8, the fully coupled kiln model assumes steady-state conditions. The periodic thermal response, when following any given point of the rotating wall, is modelled by adding a convective term in the conduction equation. Density in the governing equations is calculated by considering constant pressure. Density gradients in the flow field are allowed but only due to temperature and species gradients. Flow compressibility is ignored. This is an accurate assumption for low Mach number flows. The only region in the .domain where it may be 27 inappropriate to model the lime kiln in this fashion is near the burner exit where Mach numbers as high as 0.6 may be reached. Nevertheless the flow throughout the domain is treated as incompressible. Physics in the turbulent, chemically reacting flow is governed by the laws of conservation of mass, momentum (Newton's second law) and energy (first law of thermodynamics). Reynolds averaging of variables is used in dealing with the turbulent field. The steady-state conservation equations for the flow under consideration may be written as, v.[pff-r V0 = s e where the unknowns £, the diffusion coefficients (3.i) { and the source terms of the 11 partial differential equations to be solved are given in Table 3.1. Mass conservation for example is composed of a convective term only ( V • [pV] = 0). The three components of velocity for V are u, v and w in each of the coordinate directions x, y and z respectively. The source term in the momentum equations is the sum of pressure and gravitational forces (buoyancy effects) as well as 3 additional terms arising from the turbulent and laminar stresses. Notation T in the last of these terms refers to the momentum equation under consideration (i = 1,2,3). The diffusion process is driven by turbulence. 28 EOUATIONfS) £ CONTINUITY l 0 0 MOMENTUM V Me - VP+pg -\v{p ) -\v{fi.V • v)+ V • L (3) K 3 TURBULENT KINETIC ENERGY K DISSIPATION RATE SPECIES (4) ENERGY BL £ m , h 3 ^ dx J j G-p£ C G£ K X C p£ K 2 2 Me Sm, Me Pr, Table 3.1: Components of Equation 3.1 The much smaller molecular dynamic viscosity p is added to the turbulent viscosity p, (p -p, + p). As mentioned above, turbulence is simulated through the standard twoe equation model of Launder and Spalding together with the wall function approximation [46]. Differential transport equations are solved for the turbulent kinetic energy K and the turbulent energy dissipation rate £. The eddy viscosity concept is used to represent the turbulent stresses that result from Reynolds averaging and the transported turbulent quantities are related to the turbulent viscosity by, 29 (3.2) M,=C p— fl £ Modelling constants are the standard values, as summarized in Table 3.2. The generation term. G is given by, + G=ju, + \dz j C 0.09 1.44 dy dx + w du dv dw^\ + + dx dz J dz dy 2 1.92 1.0 1.22 Sm t Pr, 1.0 1.0 (3.3) Table 3.2: Constants in the turbulence model [46] Four species transport equations (for mass fractions m through m ) are solved for 0 , i 4 2 CFL;, CO2, and H2O. N makes up the remaining composition. Generation and 2 consumption of species proceed according to the following reaction, CFL + 20 => C0 + 2H 0 2 2 2 (-50 010 kJ/kg CH4) (3.4) The number that follows in brackets is the enthalpy of reaction. The source term S in m 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, S =Am mi (3.5) iP 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 MW is the molecular weight of each species and R is the universal gas constant. t 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 S in the energy equation h 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 = x^{ )- xh{s) ds a s (3.7) a 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. Thefirstterm on the RHS is the gain in intensity by emission. Here, i' u £ i dX = — ih denotes the black body intensity . The second term on the RHS is the loss in intensity by absorption. The absorption coefficient a is generally a function of wavelength, temperature, x 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 • q ) where, r v • L^[ ^ • f = 4 <3-8) The spectral intensity isfirstintegrated 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 nonNewtonian 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 flowfieldis 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 C0 within the bed (for the reaction rate model only) is assumed 2 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 heatingfromthe hot flow (and wall) the bed undergoes an endothermic calcination reaction according to, CaC0 (s)=> CaO (s) + C0 (g) 3 2 (+1800 kJ/kg CaC0 ) 3 (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 M occupying a volume V enters the bed inlet, 0 the change in bed composition due to the above reaction is illustrated in Figure 3.2. ® CaO (6) CaC0 (s) 3 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 p is used. The value of p CaC0} CaCOi used here is 1680 kg/m. Its true density is 2800 kg/m at 25 C. Similarly the effective density in 3 volume V is p 3 2 (940 kg/m). Second, CO2 gas is assumed to be instantaneously 3 Ca0 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,+V 2 (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 •wer 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 <p is defined as the fraction of the original mass of CaC03 that remains as CaC03 as volume V is transported through the bed. In addition m is defined as the fraction of the original mass 2 of CaC03 converted to CaO and m as the fraction of the original mass of CaC0 2 3 converted to CO2 so that, m +m -\ 2 3 (3.11) 36 At any place in the bed, the mass of CaCCh and CaO in volume V is (pM and 0 (\-(p)M m 0 2 respectively. Similarly the mass of CO2 that has escaped is (\-<f>)M m . 0 3 The effective density p of the solids may be derivedfromequation 3.10. The result is, P=<t>PcaCO +Q-</>)PcaO i ( • 3 1 ) 2 The solids volumefractionof CaC03 is in fact <f>, and for CaO, (1 - (/>). The solids mass fraction for CaC03 and CaO are given below as / and g respectively. / = (3-13) (p + (\-<p)m 2 <p + (l-<p)m 2 During this reaction, if r ( r > 0 ) is the rate of change of mass for CaC03 then from mass conservation of this species, «*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. D is the granular diffusion coefficient. This will g be described in greater detail below. On a per volume basis equation 3.15 becomes, D(j) Dt where R = r/(V/? -V-W=-JI (3.16) ) is the reaction rate (in units of 1/sec). Details on the form of R CoCOj will be given in section 3.2.4. Since steady-state conditions are assumed, equation 3.16 may be recast as, 3 7 VV<P-V[D (3.17) V<p]=-R The procedure for obtaining the velocity field V will be described in section 3.2.3. The diffusion coefficient D follows from Hsiau and Hunt's kinetic theory G analysis for granular materials [55]. (3.18) 8(l + e , ) v £ » The granular diffusivity D (and granular thermal conductivity K which is used G G later) was found to increase with the square root of the granular temperature E [55]. Granular temperature quantifies the kinetic energy of the flow due to granular motion. It is defined as the ensemble average of the square of the fluctuating velocities in the bed due to its bulk motion. This kinetic energy is generated by shearing or vibration of the flow and is then diffused into the bulk of the material. At this point it is important to note that in the current research the varying granular temperature in the bed is imposed a priori. The distribution in the 3-Dfieldwill be given in Chapter 4. A coefficient of restitution e between particles (of diameter d) accounts for the inelastic collisions. p Solids fraction v is the ratio of bulkflowdensity (which includes voids etc.) to the particle density. In the context of the bed granular diffusivity this is assumed constant (v = 0.6). A radial-distribution function g (v) (a correction factor) is used in equation 3.18. a Carnahan and Starling [57] proposed the following empirical form. (2-v) go<y) = 2(1 -v) Lun and Savage [58] suggest, (3.19) 38 r g (y) \-2.5v 0 (3.20) vJ V where v* is the maximum shearable solid fraction for spherical particles. The value of g (v) used in the current research was taken to be the average of these two forms. a 3.2.2 Bed Thermal Model Starting once again with a Lagrangian approach and volume V in figure 3.2 an energy equation may be developed as follows: Rate of change of Energy in V = Rate of heat absorbed + by reaction Rate of heat lost + with C0 Rate of heat transfer by granular and bulk diffusion 2 (3.21a) DE — = Q +Qco +Qco R 2 (3.21b) ND Radiation heat transfer within the bed is ignored. Energy E is the sum of energies of CaC03 and CaO present in the volume. That is, E=E CaCOj +E =<pM C T + (l-<t>)M m C T Ca0 0 p 0 2 p (3.22) where T is the temperature. A constant specific heat C is used for the solids. The first p 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 . Q =-M h R R o (3.23) CaC0 Qco =-M m,RC T 2 0 pco (3.24) 39 Enthalpy of reaction h ^ for calcination is used in equation 3.23. Specific heat of C O 2 CaC0 C pco 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 Q C0ND in the following form. Q =W[KVT] COND (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 +K g (3.26) K has also been developed in Ffsiau and Hunt's kinetic theory analysis [55]. g K =~^ g (3.27) Substituting equations 3.22 through 3.25 into 3.21 yields (on a per volume basis), rv^rl-v.^rl-^^+^^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 divergencefreevelocityfieldis sought for. VV = 0 (3.29) To obtain a solution for the 3-components of velocity (u, v and w) in the 3directions (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 flowfield.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 velocityfieldfor 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 velocityfieldin 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-° ox dy <'> 30 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, dl 1 u = -— dy (3.31) 41 and, 3^ v=— (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 co . z dv du ox ay co = — -— z (3.33) 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»3 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, dT_ Q =K 4mp where K CaQ Ca0 (3.35) dr is the conductivity of the porous calcined layer of figure 3.3. After integration from radius r to r (from reaction interface to surface), t s 4xK (T Ca0 s r. - -T ) (3.36) t r si r r where 7]. is the core temperature and T the particle surface temperature. s Volume fraction <p defined earlier is related to the radii geometrically as follows, (3.37) The number of spherical particles in a control volume may be given by, (3.38) N Volume V is divided by the volume per particle. The solids fraction v represents the reduction in the number of particles present in the control volume due to fines etc. The total rate of heat supplied to the control volume is Q N (>0). This heat must be absorbed p by the reaction as given by equation 3.23. That is, Q N-M h p o R CaCO =0 Substituting equations 3.36 through 3.38 into 3.39, (3.39) 43 Moo R= \T,-T (3.40) t \ "CaCO, PcaCO, J\ s K Thefirstterm 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 r (or surface area) of the particle. The final 2 term in the brackets represents the reduced resistance a particle would have to incoming heat at high values of <p (i.e. a thin CaO layer). Although this term remains around 1.0 for as high a value of <p as 0.75, numerically it gets very large when <p is close to 1 and causes difficulties (unrealistically high reaction rates). This is a consequence of ignoring the heat up period of a particle in the current analysis. For a workable simplification to this model the third term is ignored. Thefinalform of the reaction rate is, (3.41) where T is the computed temperature of the energy equation (replacing the surface temperature). This would also be the experimentally measured temperature. T is the c calcination temperature (the core temperature) and <p is a switching function used when a there is no more CaC03 present in the bed (a = 10^). The calcination temperature T is allowed to vary through the bed depth c according to an estimated local partial pressure of CO2. Figure 3.4 shows the effect of temperature on the decomposed pressure of CaC03 to CO2. At about 1169 K this pressure is greater than 1 atmosphere. 44 P . I bar . 0.987 mm ?r • i • i I t f300 1250 • • • I i • i i 1200 1150 I i i . • I • • • • I • • • • I • . • . I • 1100 1050 1000 950 800 . •• I Temperature (K) Figure 3.4: Decomposition Pressure of Calcium Carbonate [59] Calcination of limestone may either take place by reaching temperatures above 1169 K in a CO2 environment or by locally attaining lower partial pressures of CO2. In the rotary kiln lower partial pressures are possible near the surface of the bed due to diffusion of CO2 to the hot flow. In the present model the distribution of T in the bed is c imposed a priori. More detail will be given in Chapter 4. Calcination is also possible by reducing the atmospheric pressure by means of a vacuum pump but this is not usually done in practice. 3.2.5 Summary of Governing Equations for Bed Model Once the stream function is evaluated from, 3 ¥ 2 d¥ 2y (3.34) 45 The velocityfield(with an imposed w-component) is calculated by, and, dV u — —dy (3.31) v= (3.32) dx This mass conserved divergencefreefieldis used to solve the volumefractionof CaCCh and temperaturefieldfrom, (3.17) V V</>-V [D V<t>]=-R g V • v\pC T]- V • [KVT] = -p P CaCO) R(h CaCOi +mC 3 pcQ2 T) (3.28) respectively, where the density of the solids is given by, P = </>PcaCO, (3.12) +Q--<t>)PcaO and the reaction rate R is given by, 3vK,CaO R= CaCOy 3.3 (3.41) PcaCO n Wall/Refractories Model The wall model includes varying conductivity K w within the metal and refractories using manufacturer-supplied curve fits. Rotation is included through a convective plug flow term in the energy equation. v\p V C J-K VT]=0 w where, w p w V =ClxR w (3.42) (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 p and specific heat C w 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. 47 4.0 COMPUTATIONAL PROCEDURE 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 datafromthe 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 T21fromBarr'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 flowfieldboundary 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,firehood 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(10) kg/ms. This s 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„ and its dissipation rate f„ are imposed by, Jr„ =1.5(7 fi] 1 2 (4..) TI and (4.2) _=c{* — E where £/ 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 0 and N respectively. All inlet gases enter 2 2 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 W as follows. i JJJ^. [D.E.\kdydz = 0 (4.3) The most common weighted residual formulation has been the Galerkin method in which the weighting W and the interpolation functions (variables within the differential i 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 thefifthorder 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 (wcomponent) 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 entirefield,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 u = cor (u = u = 0). The single vorticity component in the flowfieldgives, e r z 1 d{ru ) du 6 r r dr dd (4.5) = 20) as expected. A scalar stream function may also be defined in cylindrical coordinates such that continuity is satisfied identically. That is u = dV g and - ru = dV r dr . Upon 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, ¥ =|(r -tf )+l 2 2 (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/plugflowinterface for the active layer calculation as shown in figure 4.2. 4.2.1.2 Discretization 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, »F" +1 +#F =V N (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 #p 2 dx 2 d^ 2 +^ + dy 2 = -K (4-9) where R% is the original equation at iteration level N . * * = Vdxi - + Vdyi — 2 " 2 + 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, d dV dx 2 JJJ*. 2 d 6V dxdydz = - jjJN,R^ dxdydz dy 2 + 2 (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 W of equation 4.3 is replaced i 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 oW,. dV N = JJJ dx dx dN, ox dx dN, dV ^ dy dy dN. d&¥ pxdydz dy dy N | dxdydz - j^N^V • 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 dN, dNj Z Z JJJdx 1=1 i dx 7=1 NDOF dN dNj dN, d*¥ dN, dV N = Z JJI dy dy ,„ ^ . , N A 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-byelement matrix assembly gives, NELEM NDOF NDOF z Z J JdxJ dx ^ 1=1 NELEM 7=1 NDOF = Z Z JJJ e=l 'dN, dN,. '•=1 dN, dV \ d t ^ N dN. dN, + dy dy dN, dV \dj-dy- N + Ar + N „ ' t J 2 J a J ) , ^ (4.15) and results in a system of equations for all unknowns (<W. 's) in the domain. Since the original equation being solved is linear only one matrix inversion operation is required. 62 4.2.1.3 Three Dimensional Velocity Field The solution to the discrete system of equations (4.15) with the boundary conditions given in figure 4.2 is shown in Figure 4.3. Y Figure 4.3: Scalar Stream Function Distribution in Bed Cross-Section This includes the plug flowfieldanalytical distribution (equation 4.7). The grid used to obtain this solution will be described in Chapter 5. The velocityfieldis calculated from the stream function distribution (and applied to the entire 3-Dfieldas follows, u= dv (4.16) dy dv__ ^rdN N v= and i dx ~( dx (4.17) w = 0.0038 m/s for active layer 0 m/s for walls and plug flow (4.18) 63 The velocity field solution at any bed cross-section is given in Figure 4.4a. A 3-D particle trace plot is given in Figure 4.4b. Note that motion in the axial direction is only present in the active layer. (a) Velocity Vector Plot at a Bed Cross-Section (b) Computed Particle Trace in Bed Figure 4.4: Velocity Field in Bed As mentioned earlier the 20-node element is also used for spatial discretization. Any boundary would be represented by a quadratic function (in a given direction). Although the post-processing plots given above show linear segments, they are in fact quadratic in the computations. The wall geometry for example is described exactly. 64 4.2.2 Bed Species and Thermal Model Governing equations for the species and thermal models in the bed are given by equations 3.17 and 3.28 respectively. For numerical stability, artificial diffusion is added to the physical diffusion for both equations as discussed in section 4.2. The equations become, VV(p-V-l(D +D r7<p\=-R g VV\pC T]-V[{K+K p +K )VT]=-p R(h g a (4.19) a CaCO CaCO) + m C T) 3 pco (4.20) Artificial diffusion D and K in these equations are chosen to be as small as a a possible but high enough so that no oscillations in the converged solutions are present. Their distribution is generally a function of the axial coordinate. Details will be given in Chapter 8. The velocityfieldV for the solids is supplied by the bedflowfieldmodel and is fixed in above equations. Granular diffusivity D and thermal conductivity K g g require an imposed granular temperature distribution in the current model. This description will be given in section 4.2.2.4. Reaction rate R requires an imposed calcination temperature distribution (based on an estimated partial pressure of C0 in the bed). This description will be given in 2 section 4.2.2.5. Specific heat of C0 , C ^ , in equation 4.20 varies with temperature so that 2 pc proper communication with enthalpy calculations in the hot flow model is ensured. The distribution used within the bed model for temperature T (K) is, 195 9773C10) =1005.59 + 0.1997957;^ ', 5 C pcoi J/kgK (4.21) Remaining constants required in solving equations 4.19 and 4.20 are summarized in Table 4.1. BED GENERAL: Bulk density ofCaCO, True density of CaCO, Bulk density of CaO Bulk specific heat of bed 1680.0 kg/m * 2800.0 kg/m 940.8 kg/m 1150.0 J/kgK 3 3 3 REACTION RATE: Heat of reaction (CaCO,) Diameter of granules Conductivity of CaO layer 1800.0 kJ/kg 0.0025 m 0.5 W/mK DIFFUSION: Coefficient of restitution Bulk conductivity of bed, K 0.6 0.7 W/mK GRANULAR TEMPERATURE: Peak value for G.T. function 1, E„ Peak value for G.T. function 2, £ Fractional range for G.T. function 2 G. T. depth distribution factor (o ) GT 6.0(10") m /s 24.0(10")m/s 0.3 5.0 4 4 2 2 2 2 CALCINATION TEMPERATURE: CT. at bed surface CT. at bed bottom CT. depth distribution factor (<J ) T 1083.0 1169.0 4.0 Table 4.1: Constants Used in Bed Model (for case T21 of Barr's experiments [62]) * Henein et al. [28] 66 4.2.2.1 Boundary Conditions In addition to an inlet and outlet boundary the bed models require boundary conditions at the hot flow/bed interface and the wall/bed interface as shown in Figure 4.5. Wall/Bed Interface Temperature boundary condition supplied by wall model Figure 4.5: Boundary Conditions for Bed Model For the species transport calculation an inlet value (cold end of kiln) of <p = 1 is used (where <p has been defined previously as the solids volume fraction of CaC0 ). At all 3 other faces (outlet, hot flow/bed interface and wall/bed interface) a zero normal gradient, r^- - 0, is imposed. an For the thermal transport calculation an inlet bed temperature of 300 K is used. At dT the outlet of the bed (hot end of kiln) a zero normal gradient, — = 0, is imposed. A heat dn flux distribution is supplied by the hot flow model at the hot flow/bed interface and is 67 updated at every global iteration. At the wall/bed interface a temperature distribution is updated by the wall model. 4.2.2.2 Discretization The four steps described in section 4.2.1.2 for discretizing the governing equations are followed here. Step 1 requires the non-dimensionalization of equations 4.19 and 4.20. Reference variables used are summarized in Table 4.2. VARIABLE Velocities Independent Variables Temperature NON-DIMENSIONALIZED BY Absolute value of the u,v,w Constant axial component of velocity in active layer x,y,z Length of kiln Density T L Bed surface calcination temperature Bulk density ofCaC0 P Too PcaCO 3 } Table 4.2: Reference Variables for the Non-Dimensional Bed Equations The non-dimensional versions of equation 4.19 and 4.20 are, (4.22) = -Da, V - V [ p T ] - V - ' 1 K Pe 1 1 Pe g Pe ) a W (4.23) = -Da, - Da-, J co 2 68 For clarity, notation to indicate non-dimensional values is dropped. In the diffusive term of equation 4.22 the Bodenstein number, Bo , is the ratio of convective mass transfer to diffusive mass transfer. Bo = ^ (4.24a) g U L Bo (4.24b) — = —— + — Bo Bo (4.24c) Bo g a Da is Damkohler's First, which describes the ratio of the reaction rate to the flow rate. x RT Da. = — 1 (4.25) U In the diffusive term of equation 4.23 the Peclet number, Pe, is the ratio of heat convection to heat conduction. Pe = " °> r ~ K Pc P c C U 4.26a) L ( * = ° °^ ~ P c C p U (4.26b) L 8 Pe. = " ^ ~ Pc 1 co 1 U — =— r + Pe Da 3 Pe 4.26c) L 1 ( 1 + Pe Pe (4.26d) is Damkohler's Third, which represents the ratio of heat absorbed/released by reaction to the heat transported by convection. 69 Da, = CaC0 fl (4.27) C TJJ p Da, RL 3 = m m,C„ TRL (4.28) Pco 3 2 C T U P OO CK Equation 4.27 represents the heat absorbed by the reaction of CaC0 and equation 3 4.28 the heat released from the control volume with C0 in the process. 2 In non-dimensional form the density is updated (from equation 3.12) by, (4.29) p = m + m <p 2 3 A Newton's linearization is now applied to equations 4.22 and 4.23 as follows, <p N+i (4.30) =<p +5<t> N T =T +dT NH (4.31) N 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, (4.32) Bo where, R» =V-V(t> N -V (4.33) V<t> + Da, N Bo and the linearized energy equation becomes, VV\p bT}-VN (4.34) •VST — —R Pe T where, R? =V-V\p T ]-V N N —VT Pe N (4.35) + Da, +Da, 3 „ J r 70 and R* are once again the original equations at iteration level N. Convergence for the discretized system of equations is measured through L norms as indicated by x equation 4.11. Density and source terms Da , Da , t 3g are evaluated at iteration Da ico level N. Applying the GMWR to equation 4.32 and 4.34 and integrating diffusive terms by parts gives, ffl N,V •V{fy) + — {VN,)-(Vfy) dxdydz Bo JJJ N V • V<t> N t +— Bo (V7V,.)• ( v > * ) + /V,Dfl, \lxdydz (4.36) and, JJJ •'7(p Sr)+j-{7N )-{Vdr) N i N,V • v(p T )+ N N yJ^N, )• (VT )+ N,Da N 3it dxdydz + N Da. t vixdydz co 2 + <H—N.VTndA (4.37) Contour integrals description and their evaluation will be given below (and in Appendix B). Substituting the element interpolation polynomials (of Appendix A) and summing over all elements gives the discrete system of equations with unknowns and STj at all nodal points. For equation 4.36, 5<pj 71 NELEM NDOF e=\ NELEM NDOF (VNj) dxdydz 2 j=\ 2 JJJN,V • VNj +—(VN,)• Bo 1=1 NDOF - 2 2 -JJJ z ifl^w-^ NELEM NDOF eA=l 1=1 -2 5o I i (4.38) and for equation 4.37, NELEM e=\ NELEM =2 NDOF NDOF Z Z {j^^V-V^N^j-iVN^iVN^dydz^ NDOF J. e=l (4.39) NELEM *2 eA=l The contour integrals in equations 4.38 and 4.39 are applied to the boundaries of the domain. A 9-point integration is performed (3X3) using a Gauss-Legendre quadrature procedure. Once again this results in exact solutions. Note that the contour integrals are summed over boundary elements 'eb\ These surface integrals are related to the mass flux J and heat flux q through the boundary. Further details will be given in section 4.2.2.3. In equation 4.38 the contour integral term is subtracted (rather than added) to the system of equations. This term would normally simply be dropped to achieve the zero normal gradient boundary condition, —— = 0. During the development work it was on realized that a zero gradient boundary condition was difficult to achieve with the simpler 72 approach of ignoring the contour integral in equation 4.38. In parts of the calcination zone reaction takes place (as will be shown in Chapter 8) primarily near and at the boundary (hot wall and bed surface conditions). This is in contradiction to a normal zero gradient condition unless a veryfinegrid is used to resolve the details near the boundary. In the current work the contour integral calculation is performed and the evaluated diffusive flux is subtracted at any given boundary element surface area. That is, if a condition exists where the diffusive flux is non-zero, it is subtracted from the appropriate lines in the matrix of discrete equations. This resulted in a 57% drop of the overall diffusive mass flux through the boundaries of the domain where mass flux should be zero. It did not go to zero exactly as will be described in Chapter 8. A separate matrix assembly for S<p. 's and dTj 's is performed at each local (bed model) iteration. A matrix inversion is then carried out for these two matrices through a banded guassian elimination solver. Every line (describing a discrete equation for an element node) of these matrices is influenced by numerous elements surrounding that node. This influence for all nodes requires that a complete matrix needs to be assembled for all nodes in the computed field. However, since these matrices involve influence on the change in the solution and not on the solution itself the matrices may effectively be solved in sections. In the current work information for one cross-sectional element plane is assembled at a time. Figure 4.6 illustrates this procedure. For clarity a 2-D grid (perhaps a top view of the bed) is given. The k-direction is the bed/kiln axial direction. 73 23 N O O E S O F INTEREST -o — <b fd> ^—©•- (b ? f 1 KPLANE-1 t f KPLANE I KPL*NE*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 <p and T . In the N N example given in figure 4.6 a 23X23 influence matrix would require assembly. In this manner computer memory requirements are drastically reduced yet convergence is not expected to suffer dramatically since gradients in the bed crosssectional plane are presumed to be larger than in the axial direction. It should be emphasized that the resulting solution at convergence is completely 3-D. 74 4.2.2.3 Contour Integrals for Mass and Thermal Fluxes Equations 4.38 and 4.39 require the evaluation of contour integrals at the boundaries of the domain. The detailed numerics on how these integrals are evaluated are given in Appendix B. Their form depends on whether for example the heat flux q= VT is known or needs to be calculated. Integrals are also evaluated in the prePe processing phase for areas on the boundary elements as well as in the post-processing phase when checking for overall balances. Table 4.3 gives various flow rate integrals (in dimensional form) for the post-processing phase. Convective flow rate ofCaCO, (kg/s) Diffusive flow rate ofCaCO, (kg/s) Convective flow rate of energy (W) Diffusive flow rate of energy (W) CaC0 rate of increase due to reaction (kg/s) Heat rate of increase due to reaction (W) C0 rate of increase due to reaction (kg/s) 3 2 j^Pcaco^ndA <Q-DV<pndA j!jpC TVndA p <Q-KVTndA Pc co dydz JJJ- Rdx a JJJ- PcaCO, ( CaC R h 0i + 3 o m C Pc 2 T^dydz \ll iPcaCO, y m Rdxd dz Table 4.3: Fluxes Evaluated in the Post-Processing Phase (Dimensional Form) The last integral in the table evaluates the rate of increase of C0 in an element 2 due to reaction. These results are then summed as shown in Figure 4.7 so that C0 boundary data are supplied to the hot flow model. 2 75 Figure 4.7: Transfer of C0 to Hot Flow Model 2 J-index element contributions of C 0 flow rate are summed at every I-index location and 2 released at the surface of the bed. 4.2.2.4 Granular Temperature Distribution Granular diffusivity D and granular thermal conductivity K (equation 3.18 and g g 3.27 respectively) are functions of a local granular temperature, S. As mentioned in section 3.2.1 granular temperature quantifies the kinetic energy of the flow due to granular motion. Its value in the plug flow region is zero since particles simply follow the wall rotation direction. 76 In general granular temperature is large in regions where the material exhibits fluid-like behaviour (as in the active layer) and low in regions where the behaviour is solid-like (plug flow) [69]. Experimental evidence [70] has shown that in the flow direction of the active layer (as shown in Figure 4.8) the granular temperature reaches a peak (starting from zero) then dissipates back to zero before re-entering the plug flow region. Plug Flow Region Figure 4.8: Experimentally Observed Granular Temperature Distribution at the Bed Surface In the current research the surface granular temperature is imposed by adding two parabolas as shown in Figure 4.9. They are both functions of angular position 6 in the bed as indicated in Figure 4.10. 77 Angular Position Figure 4.9: Imposed Bed Surface Granular Temperature Distribution The two peak values 5 , 5 M| mz and the position of #, are inputs to the model. Surface distribution is then evaluated by, 2^=5,+S 2 (4.40) Depth distribution is also imposed as a function of this surface value through an exponential decay (down to the plug flow interface) as shown in Figure 4.10. 78 Figure 4.10: Depth Distribution of Granular Temperature in Active Layer In addition to the three input parameters discussed above, one more is required for the exponential distribution. Input data are only estimates at this point but are chosen such that the mean surface granular bed conductivity in the active layer is about 5 times that of the bulk bed conductivity K [71]. Further details are given in Appendix C. Resulting cross-sectional distribution for the granular Bodenstein number is shown in Figures 4.11. A similar distribution exists for the granular Peclet number. 79 Y BODENSTEIN NUMBER Figure 4.11: Cross-Sectional Granular Bodenstein Number Distribution 4.2.2.5 Calcination Temperature Distribution As discussed in section 3.2.4 the bed reaction rate model requires a calcination temperature distribution T . This is allowed to vary through the bed depth according to c an estimated local partial pressure of C0 . Calcination temperature at the surface and 2 bottom of the bed is set to 1083 K and 1169 K respectively. These values correspond to a partial pressure of 0.25 at the surface and 1.0 at the bed bottom. It is important to note that a partial pressure of 0.25 at the surface corresponds to roughly twice that present in the hot flow results of the calcination zone. The bed bottom value of 1.0 is also an 80 estimate at this point. However, a control parameter for the distribution of the calcination temperature and therefore the reaction rate is described below. The variation (from bed surface to bed bottom) is imposed from an exponential distribution that is a function of bed depth. One input parameter is required for this distribution. Further details are given in Appendix D. Input to this model is calibrated to the experimental value for the degree of calcination at the bed exit. The resulting cross-sectional distribution of the calcination temperature is given in Figure 4.12. CALCINATION T E M P E R A T U R E (K) Figure 4.12: Cross-Sectional Calcination Temperature Distribution 81 4.3 Wall/Refractories Model The FEM is used to obtain a numerical solution to equation 3.42. Wall rotation is included through a convective plug flow term. In the fully coupled modelling (Barr's experiments [62] with hot flow, bed and rotating wall) the wall temperature solutions fluctuated considerably during the global iterative process. This occurs primarily due to sensitivity of the wall model to purely heat flux boundary conditions accepted from the hot flow and bed models. To aid the iterative process, a time dependent term was added to equation 3.42. The numerical time step used was such that information traveled one element along the wall rotation direction at every global iteration. It is worth noting that Bui et al. [4] resorted to the same time stepping procedure in their wall modelling. The time dependent term effectively acts as an energy accumulation term reducing the sensitivity of input heat fluxes to the temperature distribution in the wall. At the end of the global iterative process, steady state is achieved since the wall goes through approximately 3 to 4 revolutions. With the time dependent term, equation 3.42 becomes, (4.41) The velocity field has one component in cylindrical coordinates (V =cord). Constants g used in the wall model are summarized in Table 4.4. 82 1334 kg/m Density, p 1350 J/kgK Specific heat, C Rotation rate, co 0.1571 rad/sec (1.5 rpm) Numerical time 2.5 sec Step, At 3 w pw Table 4.4: Constants in Wall Model In modelling Monem's experiments [61] the time dependent term was not used and co was set to zero (no wall rotation). A 3-D conduction equation was solved. V-[K VT] W = (4.42) 0 The wall model includes varying conductivity K W within the metal and refractories using manufacturer-supplied curve fits as follows, (4.43) K =a +bJ w w where a and b are constants. Figure 4.13 illustrates the refractory and steel shell zones w w in the wall. Dimensions will be given in chapter 5 with the grid generation description. (?) Carbon Stool Shell ® Refractory Castable Insulating Castable Kiln Inner Radius Kin Outer Radius Kiln Center Line Figure 4.13: Wall Zones for Conductivity Calculations 83 Tables 4.5 and 4.6 give data for these curvefitsfor the hot flow/wall coupling [61] and the fully coupled run [62] respectively. ZONE TYPE K (W/mK) (W/mK ) 1 1 2 3 Refractory Castable (1.8m in hot end) Insulating Castable Carbon steel shell 1.24 4.45(10") 0.16 40.0 1.15C10") 0.0 4 4 Table 4.5: Wall Conductivity Curve Fit Data for Alyaser's Experiments [61] ZONE TYPE K (W/mK) (W/mK ) 2 1 and 2 Refractory 3 Steel shell 0.2475 40.0 1.45(10^) 0.0 Table 4.6: Wall Conductivity Curve Fit Data for Barr's Experiments [62] 4.3.1 Boundary Conditions Boundary conditions on the inside wall are heat flux distributions supplied by the hot flow and bed models as shown in Figure 4.14. 84 Figure 4.14: Wall Model Boundary Conditions On the outer wall, heat is lost to the surroundings by free convection and radiation. On the kiln hot end and cold end faces a zero normal gradient for temperature, —— = 0, is imposed. an Heat loss to the surroundings by radiation, q , rad is calculated from, 1 =£s°(Tl-T:) rad (4.44) where e is the emissivity at the outer surface (set to 0.8) and o the Stefan-Boltzmann s constant. T is the outer wall shell temperature and 7^ the ambient temperature (300 K). sh Heat loss by free convection, q , is calculated from, fc o =h {T -Tj /c /c sh (4.45) 85 where h is the heat transfer coefficient. The procedure to obtain the free convection heat fc transfer coefficient is described below. 4.3.1.1 Free Convection Heat Transfer Coefficient Free convection heat transfer from individual element surface areas (outer kiln curved element surfaces) is modelled as free convection from inclined planes. The angle 6, which the plane makes with the vertical, is shown in Figure 4.15. Positive Angles Negative Angles Figure 4.15: Illustration of Angle 6 for Inclined Plane of Length x Positive angles indicate that the hot surface faces downward and negative angles have the surface facing upward. 86 Experimental data have been correlated [60, 72] to, (4.46) Nu = C(G/ ?T COS* 0) a X fc h X where Nu is the Nusselt number (Nu =——). a, b and c are constants corresponding k to specific physical conditions. Pr is the Prandtl number (set to 0.8) and Gr ' is a x modified Grashof number (Gr* = GrNu) given by, Gr gq T kv (4-47) fc 2 f where g is the gravitational constant. Physical properties such as conductivity k and T +T kinematic viscosity v for air are evaluated at the film temperature, T (=— —). r Third order polynomial expressions were curve fit to allow for the variation of k and v for air. A summary of correlation constants a, b and c in equation 4.46 is given in Table 4.7. 9>0 0<0 6 > -86° laminar: 9<86° laminar: c=0.55 c=0.55 a =0.20 o=0.20 6=1.00 6=1.00 turbulent: Gr' Prcos 6 > 10 turbulent: G^*Pr>10 c=0.17 c=0.17 a =0.25 a=0.25 6=2.00 b =0.00 0>86° 0 < -86° 2 10 10 x always laminar: c=0.55 o=0.20 b =0.00 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 q in equation 4.47 is a function of the heat transfer coefficient (see fc equation 4.45) an evaluation of h for a given outer wall temperature must be obtained fc 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/mK. It was also found that anything more 2 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 turbulentfreeconvection. K Gr Pr Nu^ =0.60 + 0.387<! 2 16/ 1+ 0.599 Pr (4.48) '16 This correlation applies to both constant heat flux and isothermal surfaces. The Grashof number Gr is, Gr = s(Tsh-TM T v< f out (4.49) 88 where D is the outer kiln diameter. The shell temperature T andfilmtemperature T nm sh f are based on a bulk value (mean temperature of planes around a cylinder offixedaxial length) for this comparison. Figure 4.16 compares results between a bulk heat transfer coefficient from 12 inclined planes around a cylinder with the heat transfer coefficientfroma horizontal cylinder. For each data point shown, calculations were carried out such that the temperature was constant around the cylinder (by varying the heat flux to individual planes). to E 7 CL TJ Q> 5 c 12 Inciried plants •round a cyllndsr at constant temperature hfWAn'K) c 4 O a> o O O • 1 1 2 1 1 ' • 3 1 1 1 d-i.oo d-1.23 1 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 ± 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. (4.51) q = Cj(T l-T:) h (T -T„) s £s + fc sh For the iterative computations carried out in the wall model this is recast as, (4.52a) q = h {T-Tj eff where, h = h + £ CT(T +T \T + T„) eff fc S 2 2 (4.52b) The equation is linearized in the computations by lagging temperature T in h by eff 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 Independent variables x,y,z Velocities u,v,w Temperature 7 Conductivity K Heat transfer coefficient Heat flux q w NON-DIMENSIONALIZED BY Kiln inner diameter D Wall rotation rate times inner diameter Ambient temperature (300 K) T„ K„ = 1 W/mK K„ Ratio of K„ to inner diameter KjD in in K„TjD in Table 4.8: Reference Variables for the Non-Dimensional Wall Equation The non-dimensional equation becomes, ?P + V • (vj)-J-V dt Pe • {K VT) = 0 (4.53) W w For clarity notation to indicate non-dimensional values is dropped. In the diffusive term Pe is a 'rotational' Peclet number representing the ratio of heat convection in the wall w rotation direction to heat conduction. p C (OD » 2 p ( 4 5 4 ) Heat loss to the surroundings is also non-dimensionalized. Equation 4.52 becomes, q = h {T-\) eff (4.55a) 91 where, h =h +St(T +llT 2 eff /c (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 (4.57) Pe.., where, (4.58) R The original equation at iteration level N (R") may again be used to measure convergence of the discretized system of equations through the L„ norm of equation 4.11. Applying the GMWR to equation 4.57 and integrating diffusive terms by parts gives, a(<sr) JJJ dxdydz dt ?)T N Pe — + N,Pe V • VT + K (VN,)- [VT )dxdydz N •JJJ N t w + jJN^yTndA w w N W (4.59) 92 The contour integral is evaluated at a boundary as shown in Appendix B depending on whether the heat flux (q =-KjVT) is known. For a zero gradient boundary condition the contour integral is ignored. For a known heat flux (inner wall boundary) the contour integral is evaluated as given in Appendix B by, jJN,K VT w • hdA = -j^N^dA (4.60) where q is the magnitude of the heat flux vector. For an outer wall boundary condition the heat flux is given by equation 4.55a so that equation 4.59 becomes, JJJ Pe,N, | j + PeJ/,V w V<ST + K.{VN,)• (V<5T)dxdydz + = -JJJ e _ PeN, 1 OLD 1 At + PewNlVw-VT» +K (VN )-(VT ) dxdydz N W I -§N dA (4.61) iq where q = h (T — \). In this general form h = 0 on the left hand side of equation 4.61 eff eff if q is known at the boundary. In equation 4.61 the time dependent term has been discretized by a backward finite difference approximation. At is the non-dimensional time step and T 0LD the value of the temperature at the previous time step (or global iteration). Substituting the element interpolation polynomials (of Appendix A) and summing over all elements gives the discrete system of equations with unknowns dTj at all nodal points. 93 NELEM NDOF e=\ NELEM NDOF Z Z \ NDOF NDOF r + z Z Z# 1=1 7=1 eb=\ Z i=i w l w J w l J 1 w^k v JJJ (=1 NDOF J ffN NDOF NELEM ^N,N +Pe N V -VN +K (VN ).(™ ) dxdydz^STj At J _ T OLD 1 ^ + Pe Ny w VT +K {VN, N w W ){VT ) dxdydz w,M r N , (4.62) Solution to this system of equations is obtained as described in the bed model (section 4.2.2.2) including the use of the element plane-by-plane relaxation technique. 4.4 C o u p l i n g and Information Exchange The grid in the cross-section with hot flow, bed and wall is shown in Figure 4.17. Details in the development of this grid will be given in Chapter 5. Information exchange and directions of transfer are shown in Figure 4.18. Interface temperature boundary conditions for the kiln are always used in the hot flow model. Heat flux boundary conditions are used for both the inner and outer surfaces in the wall model. The bed model accepts heat flux values from the hot flow side and temperatures on the wall/bed interface. Running the three models once amounts to one global iteration. Figure 4.17: Grid Cross-Section with Hot Flow, Bed and Wall Figure 4.18: Data Transfer Between Models 95 In the hot flow cross-section the domain is divided into 28x28 computational cells (volumes). The 56x56 grid shown in Figure 4.17 is used for geometry description. That is, for every volume surface a parabola describes the geometry in a given direction. The 20-node element description in the bed and wall models also describes the geometry by parabolas. Kiln geometry shown in Figure 4.17 is not only exact but also identical between models. In the bed cross-section the domain is divided into 14x7 computational cells (elements). Once again Figure 4.17 shows a 28x14 grid (double the elements in each direction). Temperature exchange at the bed/hot flow interface is illustrated in Figure 4.19. An element surface consisting of 8 nodal points communicates with 4 surface volumes. (There are also twice as many control volume surfaces in the axial direction compared to the number of element faces in the bed model.) Temperatures at the center of each volume face are evaluated in the bed model through linear interpolation using a 4-node element description. That is, only the four corner nodes (infigure4.19) are used for temperature exchange data. 96 Figure 4.19: Transfer of Temperatures and Heat Fluxes between the Finite Element Description and the Control Volume Discretization The temperature for example at the center of control volume surface 4 (with area A ) in figure 4.19 is evaluated by, 4 T =r(^/7) = r(0.5-0.5) = X^,(0.5-0.5)r,. 4 (4.63) where £ and rj are the computational element coordinates (shown in figure 4.19) and N ( are the linear element shape functions. Heat flux data exchanges from the hot flow model to the bed/hot flow interface are evaluated by accounting for the four areas present. In general, the heat flux boundary condition, q , for the bed model is calculated by, BC EM, (4.64) 97 where q 's are the heat fluxes calculated by the hot flow model corresponding to surface { areas A 's. t Temperature and heat flux exchanges at the hot flow/wall and bed/wall interfaces are evaluated in the same manner except that the volume to element interfaces are not always 2 to 1 as indicated in figure 4.17. The number of control volume faces in the hot flow model for example, communicating with one wall element in the rotational direction is 7. Further details of coupling and information exchange as they relate to the computational grid and overall code structure will be given in chapters 5 and 6 respectively. 98 5.0 GRID GENERATION This chapter is dedicated to the topic of grid generation. Specifically, the two sets of grids generated for the hot flow/wall coupling and the fully coupled run (hot flow, bed and wall) for which experimental data exists on the UBC pilot kiln [61, 62]. A brief description of the pilot kiln isfirstgiven in section 5.1. The grid generated for the hot flow/wall coupling uses an algebraic grid generation approach (cylindrical coordinates). This is covered in section 5.2. In the fully coupled case of section 5.3 a Laplace type solution is obtained in transforming the computational coordinates to the physical domain. 5.1 UBC Pilot Kiln The pilot kiln from which experimental data was obtained and used here for comparisons to the modelling, wasfirstdescribed by Brimacombe and Watkinson [18]. The kiln has been refurbished and re-instrumented for various trials [61, 62]. The 5 m pilot kiln with an inner diameter of 0.410 m is shown in Figure 5.1. It is instrumented with suction thermocouples along its length (for gas temperature measurements) in addition to shell/refractory thermocouples. Thermocouples also access the bed material, but since they alternate between the hot gas and bed during wall rotation, a lumped capacitance response equation was used to approximate bed temperatures [62]. 99 Figure 5.1: UBC Pilot Kiln Details of the firing hood is shown in Figure 5.2. The burner is a simplified version of an industrial burner. Eight primary air nozzles surround the central fuel port within the burner. In addition, 8 air ports outside the burner supply secondary air. Figure 5.2: Firing Hood of UBC Pilot Kiln 100 5.2 Grid Generation for the Hot Flow/Wall Coupling A n algebraic grid generation approach was used in the generally cylindrical domain for the hot flow/wall coupling. Dimensions of the geometry for both the hot flow domain and wall are given in Appendix E. The outer burner square geometry was accounted for, as shown in Figure 5.3. Figure 5.3: Grid in a Cross-Section for Hot Flow/Wall Coupling This grid in the cross-section was maintained for the length of the kiln (along the symmetry or z-axis). Figure 5.4 shows the grid on the center-line vertical plane. Twenty segments were used to describe the geometry (domain segmentation strategy). 101 tl Figure 5.4: Grid in Central Vertical Plane for Hot Flow/Wall Coupling The number of volumes used in the r, 6 and z directions for each segment are given in Table 5.1. SEGMENT # GEOMETRY 1 2-9 10 11-18 19 20 Fuel nozzle Primary air nozzles Burner Secondary air nozzles Firing hood Kiln r 2 3 8 3 11 24 9 z 24 3 24 3 24 24 3 3 10 3 4 65 Table 5.1: Discretization in Each Segment for the Hot Flow/Wall Coupling Grid The resulting 3-D grid (without the wall) is shown in Figure 5.5. The kiln (segment 20) is surrounded by the wall. 102 Y Figure 5.5: Hot Flow Grid for Hot Flow/Wall Coupling The 4x12x31 (r,0,z) 3-D wall grid is shown in Figure 5.6. A total of 1488 elements were used with 7236 nodal unknowns (20 degrees of freedom in each element). Y Figure 5.6: Wall Grid for Hot Flow/Wall Coupling 103 There are half as many elements used in both the 6 and z directions compared to the number of control volumes used in the hot flow grid of segment 20. That is, 4 control volume surface areas interact with each finite element surface area in the wall grid. 5.3 Grid Generation For The Fully Coupled Run In the fully coupled run (hot flow, bed and rotating wall) Laplace type solutions are obtained for the grid cross-sections of the hot flow, the bed active layer and the bed plug flow region as shown in Figure 5.7. These grids are surrounded by the wall grid, obtained through an algebraic approach, as was done in section 5.2. Hot Flow Bed Active Layer Bed Plug Flow Wall Figure 5.7: Grid Cross-Section with Hot Flow, Bed and Wall 104 The differential equation method [75, 76] is used to generate these grids where the mapping between the physical space and computational space is typically controlled by a Laplace or Poisson equation. Coordinate points (x,y) on the boundary of the physical domain are specified and the interior points are determined by solving, a£ a£ 2 2 dx dy dn dn 2 2 dx 2 (5.1) 2 2 + (5.2) dy 2 where (^,T]) are the computational space coordinates and V(^,rj) and Q(£,rj) control the grid spacing in the physical domain (prescribed source terms). The transformation to computational space is done by interchanging the roles of the independent and dependent variables. The system of equations to solve in computational space is, dx 2 „ -2a, a ^ - 2 a d? 3 JC dx 2 2 d@n 2 ^ d@n ' dn = RHS, (5.3) = RHS (5.4) 2 + a ^ J dn 2 y dy_ where, (5.5) [drj. dx dx a, = 2 a = dy dy (5.6) + ——— d^dn d?dn dx (5.7) 3 RHS r dg drj (5.8) 105 RHS = -J' 3£ y and, J= dx dy (5.9) dn dy dx (5.10) 9£ dn 9£ dn Equations 5.3 and 5.4 are discretized by second order central finite difference equations and solved on a uniformly-spaced rectangular grid in computational space. The resulting system of equations are solved for, through a Gauss-Seidel iteration. Although manual input of and Q(^,n) was attempted, grid control was found to be easier by setting these to zero everywhere and controlling the boundary point positioning (through expansion factors). Further details will be given below. All boundary curves in the geometry of figure 5.7 were described analytically a priori, so that a minimum of input data is required to produce the grids. Input data to define the bed geometry for example, is given in Table 5.2. Definitions of these input variables are given, with the aid of Figure 5.8. INPUT VARIABLE m R f DESCRIPTION Internal kiln radius Fractional fill VALUE COMMENT 0.205 m 0.12 Bed cross-sectional area divided by 7tR . For h 2 IN &BED fal f J mg L Angle of repose Fraction of active layer thickness Fraction for minimum Active layer gap Ratio of plug flow wall distances Table 5.2: Input Data for Bed Geometry 30 deg 0.14 0.5 0.5 BED Angular bed position to the vertical Kl = fal BED h K =L K, g f =^w g ( See g- fi 5 -) 8 (see fig. 5.8) (see fig. 5.8) 106 F i g u r e 5.8: B e d G e o m e t r y Description The number of boundary nodes to be used along with expansion factors are also supplied along the curved boundary segments. A boundary segment here is any curve in figure 5.8 between ' X ' symbols. An analytical function is prescribed for each such segment that is a function of the input variables of Table 5.2. The number of nodes and expansion factor used for each segment are chosen with specific grid topology in mind. For example, a finer grid (compared to other regions in the bed) is required at the lower left wall region of figure 5.8 for the expected high temperature gradients due to wall heat regeneration. A finer grid in the active layer with a smooth transition to the plug flow grid region is also accounted for. 107 The solutions obtained above are used to describe the 20-node 3-D elements used in the bed. The resulting 3-D grid for the bed is shown in Figure 5.9. Figure 5.9: Bed Grid for Fully Coupled Case The total number of elements in I, J and K directions is 14x7x31. A total of 3038 elements with 14504 nodal unknowns was used. The hot flow cross-section grid shown in figure 5.7 is obtained as mentioned earlier by a separate solution to equations 5.3 and 5.4. The boundary data coordinates are uniformly distributed along the wall and are prescribed by the bed information at the hot flow/bed interface. At the burner and fire hood region the grid must account for other geometries, including nozzle inlet pipes and the square outer burner surface. A forth solution to equations 5.3 and 5.4 is obtained for this region. The results are given in Figure 5.10. 108 Y Figure 5.10: Cross-Section of Grid at Burner/Fire Hood Region Coordinates within the computational domain for various pipe geometry, etc. are treated as known points in the discretized field. In the kiln axial direction a linear interpolation is done between the grid of figure 5.10 to the region where the bed starts (hot flow grid offigure5.7). The resulting 3-D hot flow grid is shown in Figure 5.11. The 'BED INTERFACE' interacts with the bed grid of figure 5.9. 109 Y \ s BED INTERFACE Figure 5.11: Hot Flow Grid for the Fully Coupled Case The number of control volumes used in various segments is given in Table 5.3. SEGMENT # GEOMETRY 1 2-9 10 11-18 19-22 23 Fuel nozzle Primary air nozzles Burner Secondary air nozzles Firing hood Kiln I J K 4 2 to 4 20 2 12 to 32 28 4 2 to 4 20 2 12 to 56 28 3 3 10 3 4 65 Table 5.3: Discretization in Each Segment of the Fully Coupled Hot Flow Grid 110 Finally, the wall grid (also shown in cross-section in figure 5.7) is obtained by the algebraic grid generation approach described in section 5.2. Discretization in the r, 0, z directions is 4x16x31 for a total of 1984 elements with 9648 nodal unknowns. Ill 6.0 CODE STRUCTURE Overall code structure and convergence criteria for the two coupled runs (hot flow/wall coupling and the fully coupled run introduced in previous chapters) are described in this chapter. 6.1 Overall Structure of Model and Convergence of Hot Flow/Wall Coupling The hot flow and wall models are coupled to obtain a solution for Alyaser's experimental run #6 for comparison [61]. Besides relaxation factors within each code, overall relaxation factors are required for global data transfer. A flowchart for this coupling is given in Figure 6.1. A guessed inner wall temperature T (a constant value) is used for the first global iteration. After the first run w of the hot flow model a heat flux distribution is available for the wall model calculation. The wall model run produces new temperatures and the procedure continues until convergence is obtained. 112 AssumeT, Run Hot R o w Model Figure 6.1: Global Loop for Hot Flow/Wall Coupling Forty such global iterations were required to achieve a converged coupled solution. Computing times for this run are given in Table 6.1. I ' Global Iteration: s Hot Flow Model Run Wall Model Run Global Iterations # 2 to #40 Each Hot Flow Model Run Each Wall Model Run NUMBER OF ITERATIONS 1 global iteration 13000 local iterations 50 local iterations 39 global iterations 100 local iterations 5 local iterations Table 6.1: Computing Times for Hot Flow/Wall Coupling (Based on a 3 GHz Pentium IV Processor) CPU TIME 1160 min 1138 min 22 min 429 min total 9 min 2 min 113 In the first global iteration each sub-model (hot flow and wall) is run to convergence. With realistic solutions in both submodels (based on the initial guess) the remaining global iteration steps do not necessarily require complete convergence for each run of the hot flow and bed models. In the hot flow sub-model 100 local iterations were made (using the successive line-by-line relaxation method) and 5 local iterations were used (for the plane-by-plane relaxation technique) in the wall model. Beside the two global relaxation factors for temperature and heat flux, two additional relaxation strategies had to be implemented in the wall model (no time term was used in the wall conduction equation here). In the first 10 to 15 global iterations excessively large heat fluxes (both positive and negative) were present as input to the wall model. These were minimized as follows, q'i=f,q, (6-D where q are the heat flux values evaluated by the hot flow model and q the reduced t { values supplied as the boundary condition to the wall model. At each wall surface area, q is multiplied by f t q (0<f <1). This was eventually set to 1.0 as heat flux value q variation diminished. Two different refractories were used in the UBC pilot kiln for this run (see figure 4.13 and table 4.5). To avoid convergence difficulties during global iterations, the high conductivity in the hot zone (refractory castable) was initially set to the conductivity of the remaining kiln refractory (insulating castable). This conductivity was gradually adjusted to the hot zone values as follows, 114 'w =a +f (a -a ) (6.2) a l ^ Here a and b w> = b Wl W r 2 Wi Wi + f r K - b J (6.3) are the coefficients in the manufacturer supplied curve fits for conductivity in the hot zone (K - a +b T). w> w w Similarly a ^ and b w Wi are the coefficients for the insulating castable zone. f in equations 6.2 and 6.3 eventually goes r 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 . Relaxation factors for temperature and r 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,/ -T _\ C (6.4) iJC I "Aq ERROR"=Xk,/c-?,/ -.| <-) 6 5 C 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/m (corresponding to 2 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 C0 transfer to the 2 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 submodel 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 „ . T „ 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. I ' Global Iteration: s Hot Flow Model Run Bed Model Run Wall Model Run Global Iterations # 2 to #54 Each Hot Flow Model Run Each Bed Model Run Each Wall Model Run NUMBER OF ITERATIONS 1 global iteration CPU TIME 3728 min 20000 local iterations 20 local iterations 40 local iterations 3600 min 83 min 45 min 53 global iterations 2756 min total 100 local iterations 6 local iterations 8 local iterations 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 mass flow rate) boundary data between global iterations 2 (five parameters plotted in figure 6.4) were made through equations similar to 6.4 and 6.5. 7 6 . 4 22 <D O 2 tc ^ ^ W . ERRORT ERRORT Wall* ed to Hot Flow mjmi » ERHORq Hot Flow tD WallflEd JJJ1 5 z) o 1 ***• 3, ERROR q Bed to W«ll -2 -3 _ -5 I 10 20 30 40 V S w • 50 ERHOR 0O2 Bed ID i HB Flow i 60 Global Convergence Number - IG 70 80 Figure 6.4: Global Convergence of Fully Coupled Run At convergence about two orders of magnitude drop occurred in these parameters. The results were not only compared to available experimental data but also to an overall kiln heat balance calculation. These comparisons will be covered in chapter 8 and give confidence to the solution obtained here. 118 7.0 INVESTIGATION OF AERODYNAMICS IN THE HOT END OF THE KILN Combustion in the rotary kiln turbulent diffusion flame is not controlled by a chemical or thermal limitation. It is dependent on the rate at which the core flow (primary air and fuel) mixes with the secondary stream. "Mixed is burnt" is the phrase often used to describe the flame conditions. This chapter investigates the accuracy of current modelling in predicting the flow patterns in the hot end. Various experimental conditions were modelled for the confined jet problem. Comparison to available experimental data was made for the core flow spread rate, center line velocity decay, mean velocity distributions and turbulent stresses. A range of possible error involved in using the standard K-e model [46] is estimated and used to indicate the effect on hot flow conditions on an industrial kiln. 7.1 Aerodynamics, Turbulence and Combustion The importance of aerodynamics in the rotary kiln was emphasized in chapter 1. Turbulence modelling is important since turbulent intensities are ultimately linked to where the fuel burns. If for example the turbulent shear stresses are underestimated, a lower radial momentum transfer results in a lower spread rate for the core flow. The flame would be longer than actual conditions and fuel would be predicted to burn further downstream into the kiln. 119 Figure 7.1 compares isothermal water modelflowdiagrams of Moles et al. [6] for a 1:24 scale industrial kiln with 2-D axisymmetric computations performed as part of the present thesis work. In the numerical computations uniform profiles were used as inlet boundary conditions. This holds for all cases described in this chapter. Boundary conditions for the inlet, outlet and wall are as those described in sections 4.1.1, 4.1.2 and 4.1.3 respectively. Curtet numbers and diameter ratios are given for each case computed. Further details may be found in the references given for each case. In the water model of Moles et al. photographs were taken in the central vertical plane of the kiln tube. Flow visualization was possible by using polystyrene tracer particles. The arrows used to portray the flow patterns in figure 7.1 for both the experimental data (upper half) and numerical results (lower half) indicate direction only, not velocity magnitude. It should be emphasized that this comparison is made for the sole purpose of showing the possible effect that incorrectly predicted core flow spread rates may have in the kiln hot end flowfield.Although uniform profiles were used as numerical inlet boundary conditions, the profiles on the left inlet boundary infigure7.1 depend on the flow path determined by the shape of the kilnfiringhood (which was taken into account experimentally). Under all experimental conditions an inlet vortex is formed below the primary inlet nozzle (the region covered by the numerical results in the figure). It is therefore emphasized that results shown here do not reflect the accuracy of the model in further calculations. 120 recirculation (a) Primary to Secondary Velocity Ratio = 15.3, Ct=0.73 (b) Primary to Secondary Velocity Ratio = 18.4, Ct=0.61 recirculation (c) Primary to Secondary Velocity Ratio = 23.2, Ct=0.51 Figure 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) Conditions in figure 7.1 were such that the total flow rate was the same in all cases. The velocity ratios between primary and secondary streams are indicated in the figure. Although not as critical for this comparison results given are essentially grid independent solutions. This in fact holds true for all results to follow. The approach in which grid independence is achieved is described in section 7.3. 121 The numerical results correctly predict the qualitative growth of the Craya-Curtet recirculation pattern as the Ct number is decreased from 0.73 in figure 7.1a to 0.51 in figure 7.1c. If however the recirculation zone position was under-estimated as shown in the figure, the spread rate would be too high and the results would ultimately underpredict the flame length. This can have consequences on the ability of the model to predict properly wall hot spots and heat rates to the bed and wall. A quantitative investigation on the prediction capability of the model follows in the next two sections. 7.2 Prediction Capability for the Isothermal Confined Jet Modelling of the experimental conditions of Razinsky et al. [77] were made for the confined jet. The mixing tube in Razinsky et al.'s experiments was 9.1 m (30 ft) in length with a 0.15 m (6 in) inner diameter. Two interchangeable core flow nozzles were used for a resulting secondary to primary diameter ratio of 6.0 and 3.0. Computations were made for 2-D axisymmetric conditions. Measurements on the spreading of the core flow and duct boundary layer were made from both the numerical and experimental velocity distributions. The duct boundary layer was defined as 95 % of the secondary stream velocity as shown in Figure 7.2. The position of 5 % excess velocity (center line minus secondary) was also chosen for the boundary of the core flow. 122 For J E T spread: Lj where U = U + 0.05(U -U ) s For DUCT B.L. spread : L where U = 0.95*U BL CL S S Figure 7.2: Measurement for Spreading of Jet Lj and Duct Boundary Layer L L B Results of the spread rate measurements are shown in Figure 7.3 for a secondary to primary diameter ratio of 6.0 and a Ct number of 2.6. Numerical solutions for Ct numbers of 0.7 and 4.5 are included for comparison. The predictions are generally in good agreement, however the spread rate is slightly under-predicted. This is also evident in the center line velocity decay comparison of Figure 7.4. In this figure the axial position is non-dimensionalized with the kiln tube radius and the center line velocity with its initial value (at primary jet exit plane). Due to the underprediction of the spread rate, the center line velocity distribution is slightly overpredicted. 80 . Radial Distance (mm) • 70 DUCT BOUNDARY LAYER \ ^ ^ ^ ^ ^ ^ ^ 60 50 - a=2.6 a=o.7 Ct=4.5 40 30 20 Experimental Results: • Duct Boundary Layer o Unattached Jet '- /jy^.,.-'-'''' UNATTACHED JET Razinsky+Brighton, 1971 10 n i i 0 500 1000 Distance downstream (mm) Figure 7.3: Spreading of Jet and Duct Boundary Layer (Exp. from Razinsky et ai. [77]) 0.2 - o.o t—•>—•—'—i—I—'—•—'—i—I—•—i—i—'—I—i—•—•—•—I—•—'—'—'—I 0 10 20 30 40 X/Rs Figure 7.4: Center Line Velocity Decay (Exp. from Razinsky et al. [77]) 50 124 Yule et al. [78] have carried out numerical computations of axisymmetric jet flows with the standard K - £ 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 - £ model. An additional term is included in the dissipation rate equation as follows, C.Gg p£V-^-V£ C£ 2 lP K C pe 2 (7.1) 4 K A* where in Cartesian tensor notation, (7.2) X = 0),jO) s , jk k 1 K du, duj (7.3) CO,, = 2£ 1 1 K du, " Tl s L = d y Xj + du. 1 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 s . The terms are nontj 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 C4 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 dr K£J (7.5) dx where u and u are the velocity components in the axial and radial directions x r respectively. Equation 7.5 reveals that vortex stretching in an axisymmetric jet occurs along the circumferential vorticity component co = co . This mechanism is not present xr g in the 2-D plane jet where % - 0 • 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 CENTER LINE VELOCITY Diameter Ratio = 22.3 1.0 Exp. Data [Yule+Damou, 1991] Std. k-e model Pope's correction [Ce 4=0.79] 0.8 Ucl/Uclo 0.6 0.4 0.2 h 0.0 1 0 * * • * * * * * • * 1 2 1 3 * • • • 1 4 • 1 • • 1 ' * ' ' 5 1 6 • • • • * • • • • 7 1 8 • • • • 1 9 • • • • 1 • ' • • 10 1 • • • • 11 1 * • • * 12 X/Rs Figure 7.5: Center Line Velocity Decay with Pope's Correction [81] (Exp. from Yule et al. [80]) 1 • • • • 13 1 • • ' 14 • 1 15 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 CENTER LINE VELOCITY Diameter Ratio=6.0 o Exp.Data [Razinsky+Brighton] Std. k-€ model Pope's correction [Ce 4=0.79] Ucl/Uclo 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 moreriskyit is to use Pope's correction. 128 7.3 Error Involved in Using the Standard K - £ Model Section 7.2 considered an investigation of the prediction capability for the isothermal confined jet. The standard K - £ 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 — e turbulence model. The range of error involved in using the standard K - £ 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 whetherfirstor 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 nondimensionalized with the square of the inlet center line velocity U . u v is evaluated by, 0 r uV = du dv^ 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 - £ 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 - £ 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 lOOOuv/Uo' (c)Z/Do= 1.875 10 15 lOOOuv/Uo 2 (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 10 lOOOuvflJo* 15 (b)Z/Do= 1.250 (c)Z/Do= 1.875 (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/D=0.625) as shown in figure o 7.10a. By axial position Z/D=2.5 the mean velocity distribution is more accurately 0 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 Ci value of 1.44 at an axial position of Z/D=0.625 as indicated in figure 7.1 la. By axial position Z/D=2.5 the turbulence stress o 0 would require a value of Ci 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 Ci . This range for Ci (1.30<C|<1.55) may now be used to evaluate the effect on a hot flow solution for a kiln. A generic industrial kiln has been modelled with a 2.6 m internal diameter (60 m length). A total of 25 MW is supplied with a Curtet number of 0.57. Excess 0 of 2.5 % 2 is present at the kiln exit. The radiation model used in this case is based on an optically thick medium. The radiation model was developed and implemented by Nowak [45]. An assumption of high dust levels in the hot flow region allows the use of the diffusion approximation for radiation [85]. The radiative source term in the energy equation (equation 3.1) takes the form, 0 =Vr v(oT ) 4 3a (7.7) R where a is the Rosseland mean absorption coefficient [86]. R The grid used for the generic kiln is shown in Figure 7.13. The number of control volumes used is 22x20x50 in the radial, circumferential and axial direction respectively. 136 (a) Vertical Central Plane (b) Cross Section Plane Figure 7.13: 3-D Grid Used for Generic Industrial Kiln Three solutions for values of d of 1.30, 1.44 and 1.55 are given in Figure 7.14. The temperature distribution on the central vertical plane is shown for the first 5 k i l n diameters. For the lowest value of Q (figure 7.14a) the higher turbulence levels (relative to the standard K-e model solution of figure 7.14b) results in the highest core flow spread rate. Since the flame length is shortened, the isotherms generally move upstream. The opposite effect is true for figure 7.14c. The difference in isotherm positioning between the three solutions is of the order of one kiln diameter. This suggests that the error involved in using the standard K-e prediction of wall hot spots etc. may be within one kiln diameter. model for the 137 i i i i i i i i i i i i i . i l 4 i 6 . 8 . . . i • • • • i 10 12 10 12 (a)C,=1.30 Level F E D C B A 9 8 TEMP 1391.25 1318.5 1245.75 1173 1100.25 1027.5 954.75 882 2 1 0 -1 -2 : : ' 4 6 (b)C,=1.44 2 1 0 h -1 -2 _L_i_ _L_i_ 4 6 10 12 (c)Ci=1.55 Figure 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) 7.4 Summary Remarks on Modelling of Lime Kilns The investigation outlined in previous sections only covers one of numerous topics of interest in lime kiln modelling. Investigation on the accuracy of other physical models (combustion, radiation) may also be of interest. 138 What may also be of equal or greater concern than the accuracy of the turbulence model is whether buoyancy calculations are included. Figure 7.15 for example compares solutions on the central vertical plane with and without buoyancy for the generic industrial kiln described above. The solution in figure 7.15a is axisymmetric. At about 6.6 m in the axial direction (2.54 kiln diameters) the temperature in the upper wall region is about 150 K hotter when buoyancy is included as compared to the solution with no buoyancy. 2h -2 —i—I—i—i—i—i—I—i—i—i—i—I—i—i—i—i—I—i—i—i—i 0 2 4 6 I 8 i i i i i 10 . . . . i . 12 (b) With Buoyancy Figure 7.15: Temperature (K) in Central Vertical Plane of Generic Industrial Kiln With and Without Buoyancy (Axes are in meters) (Legend as in Figure 7.14) Temperatures for both solutions in the cross sectional plane (at 6.6 m) are given in Figure 7.16. A line is drawn in these plots to estimate the volume occupied by a fictitious bed. 139 •1-0 -OS OX) OS IX) -1.0 -0.5 OX) 0.6 1.0 Figure 7.16: Temperature (K) in Cross Sectional Plane at Axial Position of 6.6 m for Generic Industrial Kiln (Axes are in meters) In the region where the bed surface would be present there exists as much as a 200 K difference between the two solutions. This error would be present if no buoyancy was included in the calculations (at least for the case considered here). Burner geometries with swirling flow has also not been considered here. This is an additional challenge for the prediction of lime kiln flows since the standard K-e model does not account for the streamline curvature effects on the turbulence structure. The streamline curvature effect may be accounted for by using the modified K-e model proposed by Launder et al. [87]. This has been used successfully for the high swirling hydrocyclone flows [88] and may be worthwhile exploring for lime kiln flows with burner inlet swirl. Limitations and improvements in the turbulence model as mentioned above is only part of the possible avenues of investigation. Influence of dust in the radiative properties for example has not been considered here. Nevertheless, an important 140 conclusion from the investigation carried out in this chapter is that the standard K-e model does a reasonable job in predicting rotary lime kiln flows. Any additional investigations must be weighed with the overall prediction capability of the model in predicting rotary kiln flows for design or optimization purposes. 141 8.0 RESULTS AND DISCUSSION Results for two coupled runs (hot flow/wall coupling and the fully coupled hot flow/bed wall case) are given in this chapter. Visualization of the 3-Dfieldsis included. The obtained information and its implication on rotary kiln physical phenomena is presented and discussed. 8.1 Hot Flow/Wall Coupled Results Coupled model results are described for the thermal measurements of Alyaser [61] on the UBC 0.41 m internal diameter pilot kiln. Alyaser's measurements include gas flow temperatures on the center-line and 60 mm above and below the center-line. Measurements were also made on the inner and outer walls along the length of the kiln. The kiln did not rotate and no bed was used in these experiments. Run #6 of Alyaser's experiments was modelled with a burner load of 77 kW and a primary (to total) air ratio of 0.4. Conditions for this run are given in Table 8.1. The Curtet number was evaluated using equation 1.5 based on a cross-sectional plane at the burner exit. Uniform flow conditions are assumed for this calculation and densities for the core flow and secondary stream are based on local temperatures. This Curtet number is only used for comparison with the fully coupled run results. 142 Natural Gas Flow 5.0 SCFM Burner Load • 77 kW Total Primary Air Flow 22.86 SCFM Total Secondary Air Flow 34.29 SCFM Excess Air 20% Ct number 0.65 Table 8.1: Pilot Kiln Trial Run #6 from Alyaser [61] Heat transfer rate imbalance data in the model are given in Table 8.2 and a comparison of results to experimentally determined heat transfer rates in Table 8.3. 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 less than 0.1 % of the total input flow rate. Heat transfer rates also compare favorably to the experimental values as shown in table 8.3. Model Hot Flow Wall Overall Imbalance Percent of Burner Load 0.9 kW 0.2 kW 0.7 kW 1.1 % 0.3 % 0.9 % Table 8.2: Model Heat Transfer Rate Imbalances for Hot Flow/Wall Coupling Heat loss through shell Heat carried with products of combustion Experiment [61] 3-D CFD calculation Difference 37 kW 36 kW 1 kW 43 kW 40 kW 3kW Table 8.3: Comparison of Heat Transfer Rates to Experimental Data of Alyaser [61] 143 Figures 8.1 and 8.2 show hot flow coupled results for the central vertical plane of the kiln. The effective flame length may be estimated from the fuel (CFLj) contours plot of figure 8.1. If for example a 1% or 2% mass fraction value is used as the flame boundary, the flame is shown to penetrate approximately 1 m into the kiln. This agrees with the burner manufacturers' expectations. Flame impingement on the upper wall refractories is also clearly visible from figure 8.1. _ 0 . 2 Fj c o •4=J o CO o CD >-0.2 i iiiI iiiiI -0.2 0 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 0.2 0.4 O.B 0.8 1 1.2 Axial Direction (rn) 1.4 1.6 1.8 >i 2 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 Axial Direction (rn) 1.2 1.4 1.6 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-affectedfieldas well as the expected difference between upper and lower inner wall temperatures. 8 % 1 2 3 Axial Position (m) 4 5 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 theflamezone 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 flowfieldis 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 offigure8.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 afixedwall 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 offigure8.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 Burner Load Total Primary Air Flow Total Secondary Air Flow Excess Air Ct number Limestone Feed Rate Wall Rotation 5.5 SCFM 88 kW 36.0 SCFM 54.0 SCFM 72% 0.73 55Kg/hr 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 Hot Flow Bed Wall Imbalance Percent of Burner Load 1.2 kW 1.4 kW 0.8 kW 1.4% 1.6% 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 ratesfromthis model to the CFD results is given in Table 8.6. The satisfactory comparison gives confidence in the computed 3-D fields. Heat carried with products of combustion Burner Thermal Load Heat Absorbed by Calcination Reaction Heat loss through shell Heat Balance Model 3-D CFD Result Difference 54.7 kW 55.1 kW 0.4 kW 88.1 kW 89.0 kW 0.9 kW 24.9 kW 25.7 kW 0.8 kW 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 Axial Position (m) 4 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 thefinalquarter 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, GWRe , effectively proportional to the inverse of a Froude 2 number) was made for both cases. The value of Gr /Re for the hot flow/wall coupled 2 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 (redirection) 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»3 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. 8 / o c l M - 1 i 3 o a* in 3 5 03 D 2 • • • 4 S s Aid a) CXra atari (m) Bed flow direction Figure 8.16: Solids Volume Fraction of C a C 0 in Bed 3 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 plugflowregion 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. 6r Ih fij a- Q_ n r 0 i i i i i 1 i i i i i i i i i i • • • 2 3 • i i i i • i • • 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 Bed Bulk Peclet Number - P E ' o w 1 Artificial Peclet N u m b e r - P E . 4 n. 3 c 3 z 5 Q_ i i i i I i i 1 2 1 1 1 1 1 1 • 1 3 • Axial Position (m) 1 1 • 1 • 1 4 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 10 was required to stabilize the scheme. 4 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. Y 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 172 temperature distribution is given in Figure 8.24. The inner wall temperature increases in the hotflowregion 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 900 - P 5 CO CD 800 f covered wall Axial direction 4.5 m 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: • 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-Newtonianfieldwas 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 • 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: • 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. • 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. • 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 • 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. • 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. • 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", 4 Int. Symp. on Comb., 789-796 (1953). th 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", 9 Int. Symp. on Comb., 7-20 (1963). th 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 101 Annual Meeting, Denver, CO. 151-161 (1972). st 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 118 Annual Meeting, Las Vegas, NV., 461-469 (1989). th 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 NonReacting 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):5576(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", 16 Int. Symp. on Comb., Cambridge, MA., 719-729 (1976). th 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 ForceARPA-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.,fifthedition (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. 82. POPE, S.B., "An Explanation of the Turbulent Round-Jet/Plane-Jet Anomaly", AIAA Journal, 16(3):279-281 (1978). 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 N for nodes / = 1 to / = 20 are as follows, t Corner nodes, / = 1-8 Ntf, ,C)= -(\ L V + tf )(\ )(\ i +m + £ C M i + W, + «*, -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 = \(l-Tl l\ 2 Ntt,nX) + &)(l + tf ) i (A.l.c) On sides, i = 13, 14, 15, 16 N,(€,T,,C) = -\(l-C Xl 2 (A.l.d) + & )Ll + 'm ) l l Description of unknown variables <p or Cartesian coordinates x, y and z are described by, 20 *(£,I7,0 = 5>,(£»7,£>, (A.2.a) I=I 20 x^,V,C) = ZNM,TJ,Ck (A.2.b) 20 AZ,T,£) = £N {4,n,c:)y 1 l 20 i=l where </>., x , y., and z. are values at nodes /'. ( ( l (A.2.c) 188 APPENDIX B Calculation of Surface Integrals Parametric equations x = x(^,n), y = y{^,n), and z = z(£, 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.... ................. (...•-••. i:: Figure B.l: Mapping of 3-D surface Tangent vectors P (on 77= constant curve) and P on (£ = constant curve) are x 2 described by, - dx ~ dy ~ dz - dx ~ dy ~ dz dn dn dn (B.l.a) r r (B.l.b) The cross-product P\~xP dictates the normal direction h to the surface and surface area 2 A is evaluated through, ^ = JJjy3 P fedn = jjdA x 2 (B.2) 189 = -(5*5), always Care must be taken so that the unit normal vector h = points P xP x 2 away from the element surface. For a known heat flux q = q the heat rate, Q, through the surface is, Q=\\qhdA P xP { = 2 j ^ x P ^ d n (B.3) If q is given by a convective flux (q = pC TV ) the heat rate is evaluated by, p = jjpC T[v-(P xP )\i&ri p l 2 (B.4) Finally, when q has the form of a diffusive flux (q = -KVT), Q = -jJKVT(P xP )dZdT l 2 1 dTd{y,z) dTd{z,x) , BTB{x,y) + Kdx d{£,r}) By B{£,n) Bz B^n) + where K v (B.5) 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 ( - °) 4 2 4 where S, and E are the parabolic functions shown in figure 4.9. Besides the angular 2 extent of the bed < 0 < 0 ) as shown in figure 4.10, three input variables are required 2 for the parabolic description of equation 4.40. Namely, the two peak values of granular temperature, !E and S m OTj and the position of 6\ (of figure 4.9). S, and E are 2 evaluated from. {e-djd-e ) 2 ( C 1 ) for d > e:, E _ (e-0\\d-e ) 2 n, 2 (6 E 2 -0\X n, (C.2.a) ~0 ) e mi 2 2 for 9 < e[, : =0 (C.2.b) 2 From the geometry, (c.3) e = -^ 6 m e =e;+^A (c.4) m 9 is evaluated from the "fractional range for s " orr t 2 &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, e A e° GT urface -1 (C.6) -1 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 <J that dictates the form of this function as shown in Gr Figure C. 1. Figure C l : Depth Distribution for Granular Temperature, E At 8 = 0 (active layer/plug flow interface) the value of S is zero. At 8 = A (free surface) E = S surface in table 4.1. . The four values required as input (S _, E , r and a ) are given m mi h0 GT 192 APPENDIX D Calcination Temperature Distribution for the Bed An exponential form for the calcination temperature in the bed is assumed. c T T T -T e i A (D.l) °< - 1 T e where A is the maximum bed thickness present. Calcination temperature T is a h c function of bed depth S . Values of the calcination temperature at the free surface, T , b and at the "bed bottom", , are given in table 4.1 along with the value of constant o . T The distribution is similar to that of figure C. 1 and its form depends on o . For known T T ^ and 7^ an increase in <r generally implies lower calcination temperatures and c r therefore higher values for the reaction rate. APPENDIX E Schematic Diagram of the UBC Pilot Kiln A simplified schematic of the UBC pilot kiln is shown in Figure E. 1. Figure E.l: Simplified Schematic of the UBC Pilot Kiln (not to scale) Details of the fire hood/burner arrangement in the hot end is shown in Figure E.2. Figure E.2: Simplified Schematic of the Fire Hood/Burner Arrangement in the UBC Pilot Kiln (not to scale) 194 APPENDIX F Overall Heat Balance Model Overall heat balance calculations were summarized in the spreadsheet given below. The first page is a summary of input and output data. HEAT BALANCE MODEL Temperatures in K Heat Rates in kW Enthalpies in kJ/Kgmol INPUT DATA m_CH4 m_CaC03 fract. CaC03 consumed % excess air Q_h,l 6.33 55.00 0.905 72.18 2.00 OUTPUT DATA (kW) kg/hr kg/hr (1-phi) % kW T hf.in TJif.out T top,C02 T b,in T_b,out KgMoles of Air = 288.7 835.0 1089.0 300.0 1089.0 K K K K K 3.4436 Kgmol / Kgmol CH4 Burner Thermal Load Q_hf,r 88.13 (complete combustion assumed) Bed Thermal Load Heat absorbed by React. 27.49 24.88 (equals Q_b,r if phMD.O) Q_b,r Kiln Shell Loss Q_w.l 15.82 Hot Flow Input Rate Hot Flow Output Rate QJif.in Q_hf.out 16.12 49.96 Bed Input Rate Bed Output Rate Q_b,in Q_b,out 5.27 11.51 Heat Rate of C02 (top) Q_top.co2 Heat Rate of C02 (hf_ou Q_b,co2 6.60 4.75 <== bulk estimate Heats of Formation Enthalpies at T (K) CH4 02 N2 C02 H20 (g) CaC03 CaO -74850 0 0 -393520 -241620 -1207000 -633559 Molecular Weights Kg/Kgmol CH4 02 N2 C02 H20(g) CaC03 CaO 16 32 28 44 18 100 56 a 11.981 30.372 28.590 45.076 30.879 115.000 64.400 b 31.663 2.093 1.884 4.396 5.149 0.000 0.000 h=aT+b(10A3)T*2+c(1O«5yr+d c -3.964 1.674 0.502 8.623 0.000 0.000 0.000 d 4979 000 -1120.000 -236.000 -7362.000 239.000 0.000 0.000 enthalpy at 10030.94 8678.468 8619.556 9354:61 9898.194 34270 19191.2 Q_b298 0.791 298 kW Heat Released in H.F. Reactants at T_hf,in species CH4 02 N2 288.7 moles 3.4436 12.947936 heat form. -74850 0 0 h at T h at 298 sensi. heat kJ/Kgmol CH4 9703.901 8402.684 8348.818 10030.94 -327.0367 8678.468 -275.784 8619.556 -270.7385 -75177.037 -949.68964 -3505.5048 Sum of Reactants => Products at T_hf,out species C02 02 N2 H20 -79632.231 kJ/Kgmol CH4 835.0 moles 1 1.4436 12.947936 2 heat form. -393520 0 0 -241820 h at T 34373.81 25900.39 25010.13 29612.98 h at 298 sensi. heat 9354.61 8678.468 8619.556 9898.194 25019.2 17221.92 16390.58 19714.78 Sum of Products => Net Heat Net Heat -495993.31 kJ/Kgmol CH4 -30999.582 kJ/Kg CH4 Heat Rate -54 484865 kW kJ/Kgmol CH4 -368500.8 24861.568 212224.13 -444210.43 -575625.54 kJ/Kgmol CH4 Heat Absorbed in Bed Reactants at T_b,in species 300.0 moles CaC03 heat form. 1 h at T -1207000 34500 h at 298 sensi. heat 34270 kJ/Kgmol CaC03 230 -1206770 Sum of Reactants => Products at T_b,out species 1089.0 moles C02 CaO CaC03 0.905 0905 0.095 (C02 at TJrf.out) heat form. h at T -393520 34373.81 -633559 70131.6 •1207000 125235 -1206770 kJ/Kgmol CaC03 835.0 h at 298 sensi. heat 9354.61 19191.2 34270 kJ/Kgmol CaC03 -333493 23 -527269.83 -106023.33 25019.2 50940.4 90965 Sum of Products => Net Heat Net Heat 239983.615 kJ/Kgmol CaC03 2399.83615 kJ/Kg CaC03 Heat Rate 36.6641634 kW -966786.38 kJ/Kgmol CaC03 Overall Balance Heat Rate Loss (other than wall) Heat Released in H.F. Heat Absorbed in Bed 2.00 -54.484865 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 Q_hf,r 88.13375 (complete combustion assumed) Bed Thermal Load and Q b,r Bed Thermal Load (max) Q b,r 179921 kJ/Kgmole CaC03 1799.21 kJ/Kg CaC03 27.487931 kW 24.876577 kW <== Enthalpy of reaction used in Bed Model (equals to Bed Thermal Load if f=1) Hot Flow Input Heat Rate Reactants at T_hf,in 288.7 species CH4 02 N2 moles h at T 1 3.4436 12.947936 kJ/Kgmol CH4 9703.901 8402.684 8348.818 9703.901 28935.48 108100 Sum of Reactants => Q_hf,in=> 146739.3 kJ/Kgmol CH4 9171.209 kJ/KgCH4 16.11932 kW Hot Flow Output Heat Rate Products at T_hf,out species C02 02 N2 H20 835.0 moles h at T 1 1.4436 12.947936 2 kJ/Kgmol CH4 34373.81 25900.39 25010.13 29612.98 34373.81 37389.8 323829.6 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 species 300.0 moles CaC03 hatT 1 kJ/Kgmol CaC03 34500 34500 Sum of Reactants => Q_b,in => 34500 kJ/Kgmol CaC03 345 kJ/Kg CaC03 5.270833 kW Bed Output Heat Rate Products at T_b,out species CaO CaC03 1089.0 moles 0.905 0.095 h at T kJ/Kgmol CaC03 70131.6 125235 63469.1 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 species C02 1089.0 moles 0.905 hatT 47730.31 kJ/Kgmol CaC03 43195.93 Sum of C02 Prod. => 43195.93 kJ/Kgmol CaC03 431.9593 kJ/Kg CaC03 Q_top,C02 =====> 6.599378 kW Hot Flow Output Heat Rate of Bed C02 Products at T_hf,out species C02 835.0 moles 0.905 hatT 34373.81 Sum of C02 Prod. => Q_b.C02 =====> kJ/Kgmol CaC03 31108.3 31108.3 kJ/Kgmol CaC03 311.083 kJ/KgCaC03 4.752656 kW
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mathematical modelling of lime kilns
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mathematical modelling of lime kilns Georgallis, Mike 2004
pdf
Page Metadata
Item Metadata
Title | Mathematical modelling of lime kilns |
Creator |
Georgallis, Mike |
Date Issued | 2004 |
Description | 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 κ-ε 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 CaCO3 and CaO. Carbon dioxide is released to the hot flow. 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 = ΩxR). 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. |
Extent | 14272182 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0103853 |
URI | http://hdl.handle.net/2429/15899 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2004-901890.pdf [ 13.61MB ]
- Metadata
- JSON: 831-1.0103853.json
- JSON-LD: 831-1.0103853-ld.json
- RDF/XML (Pretty): 831-1.0103853-rdf.xml
- RDF/JSON: 831-1.0103853-rdf.json
- Turtle: 831-1.0103853-turtle.txt
- N-Triples: 831-1.0103853-rdf-ntriples.txt
- Original Record: 831-1.0103853-source.json
- Full Text
- 831-1.0103853-fulltext.txt
- Citation
- 831-1.0103853.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0103853/manifest