{"responseHeader":{"status":0,"QTime":6,"params":{"q":"{!q.op=AND}id:\"195472\"","hl":"true","hl.simple.post":"","hl.fragsize":"5000","fq":"!embargo_tdt:[NOW TO *]","hl.fl":"ocr_t","hl.method":"unified","wt":"json","hl.simple.pre":""}},"response":{"numFound":1,"start":0,"docs":[{"ark_t":"ark:/87278/s6hd89gb","source_t":"Original in Marriott Library Special Collections, TP7.5 2012 .P33 2012","setname_s":"ir_etd","subject_t":"Cfd; Coal; Direct quadrature method of moments; Dqmom; Large eddy simulation; Particles","restricted_i":0,"department_t":"Chemical Engineering","format_medium_t":"application/pdf","identifier_t":"etd3/id/1783","date_t":"2012-08","mass_i":1515011812,"publisher_t":"University of Utah","description_t":"The Direct Quadrature Method of Moments (DQMOM) was implemented in the Large Eddy Simulation (LES) tool ARCHES to model coal particles. LES coupled with DQMOM was first applied to nonreacting particle-laden turbulent jets. Simulation results were compared to experimental data and accurately modeled a wide range of particle behaviors, such as particle jet waviness, spreading, break up, particle clustering and segregation, in different configurations. Simulations also accurately predicted the mean axial velocity along the centerline for both the gas phase and the solid phase, thus demonstrating the validity of the approach to model particles in turbulent flows. LES was then applied to the prediction of pulverized coal flame ignition. The stability of an oxy-coal flame as a function of changing primary gas composition (CO2 and O2) was first investigated. Flame stability was measured using optical measurements of the flame standoff distance in a 40 kW pilot facility. Large Eddy Simulations (LES) of the facility provided valuable insight into the experimentally observed data and the importance of factors such as heterogeneous reactions, radiation or wall temperature. The effects of three parameters on the flame stand-off distance were studied and simulation predictions were compared to experimental data using the data collaboration method. An additional validation study of the ARCHES LES tool was then performed on an air-fired pulverized coal jet flame ignited by a preheated gas flow. The simulation results were compared qualitatively and quantitatively to experimental observations for different inlet stoichiometric ratios. LES simulations were able to capture the various combustion regimes observed during flame ignition and to accurately model the flame stand-off distance sensitivity to the stoichiometric ratio. Gas temperature and coal burnout predictions were also examined and showed good agreement with experimental data.","rights_management_t":"Copyright © Julien Pedel 2012","title_t":"Large eddy simulations of coal jet flame ignition using the direct quadrature method of moments","id":195472,"publication_type_t":"dissertation","parent_i":0,"type_t":"Text","subject_lcsh_t":"Moments method (Statistics); Coal, Pulverized --Combustion -- Computer simulation; Eddies -- Computer simulation","dissertation_institution_t":"University of Utah","thumb_s":"/e6/37/e63757c220df84c6066e0a3d5efee54a1577a27c.jpg","oldid_t":"etd3 1783","author_t":"Pedel, Julien","metadata_cataloger_t":"CLR; car","format_t":"application/pdf","modified_tdt":"2017-12-11T20:33:36Z","dissertation_name_t":"Doctor of Philosophy","school_or_college_t":"College of Engineering","language_t":"eng","file_s":"/8b/69/8b699861a06ff4b7b817db47c62eef6fa8b62939.pdf","format_extent_t":"6,721,420 bytes","created_tdt":"2012-08-06T00:00:00Z","permissions_reference_url_t":"https://collections.lib.utah.edu/details?id=1288857","doi_t":"doi:10.26053/0H-JJAY-T800","_version_":1619873613070991360,"ocr_t":"LARGE EDDY SIMULATIONS OF COAL JET FLAME IGNITION USING THE DIRECT QUADRATURE METHOD OF MOMENTS by Julien Pedel A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah August 2012 Copyright © Julien Pedel 2012 All Rights Reserved The Univers i ty of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Julien Pedel has been approved by the following supervisory committee members: Philip J. Smith , Chair 06/01/2012 Date Approved Jeremy N. Thornock , Member 06/01/2012 Date Approved Jost O. L. Wendt , Member 06/01/2012 Date Approved Eric G. Eddings , Member 06/01/2012 Date Approved Thomas H. Fletcher , Member 06/01/2012 Date Approved and by JoAnn S. Lighty , Chair of the Department of Chemical Engineering and by Charles A. Wight, Dean of The Graduate School. ABSTRACT The Direct Quadrature Method of Moments (DQMOM) was implemented in the Large Eddy Simulation (LES) tool ARCHES to model coal particles. LES coupled with DQMOM was first applied to nonreacting particle-laden turbulent jets. Simulation results were com-pared to experimental data and accurately modeled a wide range of particle behaviors, such as particle jet waviness, spreading, break up, particle clustering and segregation, in differ-ent configurations. Simulations also accurately predicted the mean axial velocity along the centerline for both the gas phase and the solid phase, thus demonstrating the validity of the approach to model particles in turbulent flows. LES was then applied to the prediction of pulverized coal flame ignition. The stability of an oxy-coal flame as a function of changing primary gas composition (CO2 and O2) was first investigated. Flame stability was measured using optical measurements of the flame stand-off distance in a 40 kW pilot facility. Large Eddy Simulations (LES) of the facility provided valuable insight into the experimentally observed data and the importance of factors such as heterogeneous reactions, radiation or wall temperature. The effects of three parameters on the flame stand-off distance were studied and simulation predictions were compared to experimental data using the data collaboration method. An additional validation study of the ARCHES LES tool was then performed on an air-fired pulverized coal jet flame ignited by a preheated gas flow. The simulation results were compared qualitatively and quantitatively to experimental observations for different inlet stoichiometric ratios. LES simulations were able to capture the various combustion regimes observed during flame ignition and to accurately model the flame stand-off distance sensitivity to the stoichiometric ratio. Gas temperature and coal burnout predictions were also examined and showed good agreement with experimental data. Overall, this research shows that high-fidelity LES simulations combined with DQMOM can yield a deeper understanding of complex coal flames and their ignition mechanisms, indicate where experimental uncertainties lie and in the end, be a valuable tool for the design, retrofit and scale-up of oxy-coal boilers. iv Je souhaite dedier la realisation de ce memoire a mes parents qui ont toujours ete la pour moi et qui m'ont infailliblement soutenu. Je tiens par ces quelques mots a leur exprimer mon amour, ainsi qu'a ma tres chere soeur et a ma femme. CONTENTS ABSTRACT : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : iii LIST OF FIGURES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : viii LIST OF TABLES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xi NOMENCLATURE : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xii ACKNOWLEDGMENTS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xiv 1 INTRODUCTION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1 2 LARGE EDDY SIMULATION OF COAL PARTICLES IN TURBULENT COAXIAL JETS USING THE DIRECT QUADRATURE METHOD OF MOMENTS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 4 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3.1 Governing equations of the gas phase . . . . . . . . . . . . . . . . . . . 9 2.3.2 Governing equations of the dispersed phase . . . . . . . . . . . . . . . 9 2.3.3 Stokes number analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Experimental and simulation conditions . . . . . . . . . . . . . . . . . . . . . 21 2.4.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.2 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5.1 Single phase coaxial jet . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5.2 Stoke number analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5.3 Two-phase coaxial jet . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.7 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3 IGNITION OF COAXIAL TURBULENT DIFFUSION OXY-COAL JET FLAMES: EXPERIMENTS AND SIMULATIONS COLLABORATION: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 48 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.1 Governing equations of the gas phase . . . . . . . . . . . . . . . . . . . 56 3.4.2 Governing equations of the dispersed phase . . . . . . . . . . . . . . . 56 3.5 Validation and uncertainty quantification . . . . . . . . . . . . . . . . . . . . 68 3.5.1 Data Collaboration method . . . . . . . . . . . . . . . . . . . . . . . . 69 3.5.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.5.3 Surrogate model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.5.4 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.6 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.6.1 Ignition mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.6.2 Validation and uncertainty quantification analysis . . . . . . . . . . . . 87 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4 LARGE EDDY SIMULATION OF PULVERIZED COAL JET FLAME IGNITION USING THE DIRECT QUADRATURE METHOD OF MOMENTS : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 94 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.3 Numerical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.3.1 Governing equations of the gas phase . . . . . . . . . . . . . . . . . . . 97 4.3.2 Governing equations of the dispersed phase . . . . . . . . . . . . . . . 98 4.3.3 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.4.1 Qualitative analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.4.2 Quantitative analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 4.6 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5 CONCLUSION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 119 REFERENCES: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 121 vii LIST OF FIGURES Figure Page 2.1 DQMOM approximation of the NDF with 3 quadrature nodes . . . . . . . . . 12 2.2 Illustration of the different particle kinematic regimes as a function wave number for a fixed particle relaxation time (a) and the implications on LES modeling (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Schematic of the coaxial nozzle . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4 Gas x-velocity profile (m/s) for VR=0, no particles . . . . . . . . . . . . . . . 26 2.5 Gas x-velocity profile (m/s) for VR=1.0, no particles . . . . . . . . . . . . . . 27 2.6 Gas x-velocity profile (m/s) for VR=1.5, no particles . . . . . . . . . . . . . . 28 2.7 Comparison between LES and experimental data of gas averaged x-velocity along the centerline for different velocity ratios / No particles . . . . . . . . . 29 2.8 Location of the Nyquist cut-off and the unity Stokes wave number on the model energy spectrum for 25-micron particles (a) and Stokes number as a function of wave number for 25-micron particles (b) . . . . . . . . . . . . . . . 32 2.9 Location of the Nyquist cut-off and the unity Stokes wave number on the model energy spectrum for 70-micron particles (a) and Stokes number as a function of wave number for 70-micron particles (b) . . . . . . . . . . . . . . . 33 2.10 Instantaneous profiles of weights (# particle/m3) for 25-micron particles, VR=0 35 2.11 Instantaneous profiles of weights (# particle/m3) for 70-micron particles, VR=0 36 2.12 Time-averaged particle concentration at X = 15d [(a), (b)] and along the centerline [(c), (d)] for VR = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.13 Qualitative comparison of particle behavior for the 25-microns particles cases, LES screenshots of particle concentration on the left (Wsum), photographs from experiments on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.14 Qualitative comparison of particle behavior for the 70-micron particles cases, LES screenshots of particle concentration on the left (Wsum), photographs from experiments on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.15 Comparison between simulations and experimental data of time-averaged gas and particles velocities along the centerline for 25 m particles [(a) (c) (e)] and 70 m particles [(b) (d) (f)] . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.16 Radial profiles of particles number density taken at x = 15d . . . . . . . . . . 45 3.1 Average standoff distances with associated error . . . . . . . . . . . . . . . . . 55 3.2 DQMOM approximation of the NDF with 3 quadrature nodes . . . . . . . . . 60 3.3 Illustrative schematic of coal particle components and reactions. . . . . . . . . 62 3.4 Verification results of the Yamamoto devolatilization model with comparisons to results by the chemical percolation devolatilization model (CPD) . . . . . . 65 3.5 Integration of uncertainty in modeling parameters, boundary conditions and experimental data into quantification of an overall predictive output. . . . . . 69 3.6 A Box-Behnken design for three parameters . . . . . . . . . . . . . . . . . . . 74 3.7 Normalized weight (number of particles/m3) of each quadrature node at a given timestep of simulation 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.8 Gas temperature profiles (K) for simulations 1, 6 and 7 . . . . . . . . . . . . . 78 3.9 Gas temperature profiles (K) for simulations 8 and 9 . . . . . . . . . . . . . . 79 3.10 Profiles of temperatures and devolatilization rate at a given timestep of sim-ulation 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.11 Profiles of oxygen, raw coal mass fraction and char oxidation rate at a given timestep of simulation 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.12 Profiles of gas temperature (a) and particle raw coal mass ((b), (c) and (d)) at a given timestep of simulation 8 . . . . . . . . . . . . . . . . . . . . . . . . 83 3.13 Time-averaged values of temperatures and heat fluxes along the centerline for simulation 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.14 Time-averaged values of oxygen and raw coal mass fractions along the cen-terline for simulation 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.15 Surrogate model response surfaces for PO2;pri = 10:45% (a) and Twall = 1500K (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.16 Feasible set of parameter values (a) and model predictions compared to ex-perimental values (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 ix 4.1 DQMOM approximation of the NDF with three quadrature nodes . . . . . . . 101 4.2 Illustrative schematic of coal particle components and reactions. . . . . . . . . 103 4.3 Verification of the Yamamoto devolatilization model by comparison to the distributed activation energy model (DAEM) predictions . . . . . . . . . . . . 105 4.4 Gas temperature profiles (K) obtained by LES for various inlet stoichiometric ratios (SRs) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.5 Instantaneous profiles of each quadrature node weight (number of particles x 109/m3) for SR=0.36 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.6 Validation results for the flame stand-off distance . . . . . . . . . . . . . . . . 114 4.7 Comparison between experiments and LES results of gas temperature . . . . 115 4.8 Comparison between experiments and LES results of gas coal burnout (SR=0.36)116 4.9 Comparison between experiments and LES results of gas coal burnout (SR=0.22)117 x LIST OF TABLES Table Page 2.1 Experimental setup for the different velocity ratios . . . . . . . . . . . . . . . 22 2.2 Quadrature approximation for 25-micron particles (inlet values) . . . . . . . . 24 2.3 Quadrature approximation for 70-micron particles (inlet values) . . . . . . . . 25 3.1 Utah bituminous coal particle size distribution . . . . . . . . . . . . . . . . . 53 3.2 Average standoff distances (cm) for overall PO2=40%, secondary preheat 544 K 54 3.3 Chosen parameters and their ranges . . . . . . . . . . . . . . . . . . . . . . . 72 3.4 Parameter values for each simulation in the Box-Behnken design . . . . . . . 74 3.5 Operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.6 Quadrature approximation for 25-micron particles (inlet values) . . . . . . . . 75 3.7 Time-averaged stand-off distance for each simulation . . . . . . . . . . . . . . 88 3.8 Surrogate model coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.1 Operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.2 Coal analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.3 Quadrature approximation for 25-micron particles (inlet values) . . . . . . . . 109 NOMENCLATURE gas kinematic viscosity (m2/s) dp particle diameter (m) d primary nozzle diameter (m) DQMOM Direct Quadrature Method of Moments turbulence energy dissipation rate (m2/s3) g gas density (kg/m3) ug gas velocity (m/s) p particle relaxation time (s) hp particle enthalpy (J) internal coordinate vector wave number (m1) Kolmogorov length scale (m) LES Large Eddy Simulation n number density function (number of particles/m3/internal coordinate) N number of internal coordinates N number of quadrature nodes Nyquist cut-off length (m) up particle velocity in real space (m/s) vp particle velocity in internal coordinate space (m/s) c particle mass of raw coal (kg) h particle mass of char (kg) PO2;pri partial pressure of oxygen in the primary stream Qconv net convective and conductive flux to particle (J/s) Qrad Net radiative flux to particle (J/s) Tg gas temperature (K) Tp particle temperature (K) Tpri primary stream temperature (K) Twall wall temperature (K) f fluid characteristic time (s) U0 central jet maximum gas velocity (m/s) Ua annular jet maximum gas velocity (m/s) V R inlet velocity ratio w weight for quadrature node (number of particles/m3) h i abscissa for internal coodinate and quadrature node xiii ACKNOWLEDGMENTS I am very grateful to my family, friends and mentors who have supported me on this long, enlightening and life-changing experience. My mother and my father must share a large credit of this achievement and I cannot thank them enough for all they have done for me. I also would like to thank my sister for her support in difficult times and my wife for the patience and encouragement she gave me through all these years. Few students have the chance to have Dr. Philip Smith as an advisor. I was at first impressed with his technical insight and his ability to jump into multiple complex problems and provide helpful advice. As time passed, I came to realize that Phil not only had great technical knowledge but also extra-ordinary human qualities. Even though we had our share of disagreements, Phil always kept an open mind and never showed an ounce of arrogance. He gave me incredible freedom to explore the topics I wanted to and make my own mistakes. I greatly appreciated his continuous support and his loyalty, but above all, what truly makes Phil a role model for me is his ability to manage people by bringing out the best in them. He would systematically encourage people in their work and make them want to do better. This is something I hope I will remember and later apply in my life as I am convinced our societies would be much better off if there were more people like him. I am very grateful to Dr. Jeremy Thornock who spent a lot of time helping me and answering my questions. I have learned a lot at Jeremy's side during these past years but his critical thinking and self-discipline were the most inspiring to me. This work wouldn't have been possible without his patience and his guidance. I also would like to acknowledge Dr Jost Wendt, Kerry Kelly, Dr Jingwei Zhang, Dr Terry Ring and Dr Jennifer Spinti for their help and time. In addition, I want to thank Jery Schryver for her precious help regarding immigration issues. Finally, I would like to thank the supervisory committee for their time serving as com-mittee members. This research was permitted by the Department of Energy under Award Number DE-NT0005015. xv CHAPTER 1 INTRODUCTION Global warming, one of the largest environmental challenges of our time, is due to increased carbon dioxide levels in the atmosphere [66]. The world is currently depending on the use of fossil fuels, especially coal, for its energy supply, and will continue to do so for a long time. According to the Kyoto Protocol, by 2008-2012, Annex I countries have to reduce their greenhouse gas emissions by a collective average of 5% below their 1990 levels. Therefore, oxy-fuel combustion, as one of the promising technologies to reduce CO2 emissions, has attracted worldwide attention [1, 20, 48, 66]. Oxy-coal combustion, in which an O2=CO2 mixture replaces air, is one of the few possible capture technologies to enable CO2 sequestration for existing coal-fired boilers. Burning coal with relatively pure oxygen, together with recycled flue gases, can produce a highly concentrated (up to 95% CO2) flue gas stream, which makes carbon sequestration more economical. One issue of interest in rapidly implementing a strategy to retrofit existing air-fired burners or designing new ones is to understand how replacing air by a O2=CO2 mixture affects the kinetics, aerodynamics, pollutants formation, heat fluxes and flame ignition in the near burner region. Large Eddy Simulation has the potential to predict oxy-coal flame characteristics and to be an important tool in the retrofitting and designing processes. In a turbulent coal flame, accurate prediction of the particle behavior is critical for predicting the flame characteristics. Stokes number analyses suggest that LES has the capabilities to resolve all the turbulent length scales affecting coal particles in a coal flame, and could therefore lead to more accurate simulations. The main purpose of this research is to perform reliable Large Eddy Simulations (LES ) of oxy-coal flames. ARCHES, the LES tool developed at the University of Utah, initially only included gaseous phase simulations, so implementation of particle-modeling capabilities was 2 required. The Direct Quadrature Method of Moments (DQMOM) was chosen to represent particles in an Eulerian-Eulerian framework. This method allows the transport of multi-dimension probability density functions in a computationally cheap manner. The first step of this research project was thus the implementation of DQMOM in Arches in a stable and efficient way. The next one was to implement all the physical models related to coal particles combustion within the DQMOM approach and couple them with the gaseous phase. A non-exhaustive list of those models includes mass, momentum and heat transfer between the gas and solid phases, devolatization and char oxidation. As a general rule, simulations results should be compared to experimental data to assess the validity and accuracy of the models. This process, called Validation and Uncertainty Quantification, can be achieved be selecting relevant experimental data and comparing them with the predictions made by the simulations. This research has two main objectives. The first one is to show that DQMOM can be used in combination with LES to ac-curately model particles in complex turbulent flows. It was proposed to use Budilarto's data [14] to validate this approach. Budilarto investigated the effect of fluid aerodynamics on particle motion in the near field region of coaxial particle-laden jets. The fluid aerody-namics was modified by varying the inlet velocity ratio of the annular to central jet velocity at 0, 1.0 and 1.5. The dispersion of particles was enhanced with decreasing particle size and with increasing the velocity ratio greater than 1.0. The particle-phase radial velocity fluc-tuation was determined as an important parameter to characterize the particle-turbulence interaction and the spatial distribution of particles in coaxial laden jets. The investigation of particle size distribution effects showed that the radial spreading of the coarse particles was reduced by the presence of the fine particles. LES simulation of Budilarto's cases are presented in the first chapter and show that this LES combined with DQMOM is able to predict particle dispersion, particle clustering and gas-solid interactions. The second objective aims to show that LES coupled with DQMOM can be used to perform reliable simulations of oxy-coal flames. One issue of interest in developing oxy-coal technology is to understand and predict the effects of near burner zone environment, now 3 consisting of O2 and CO2 instead of N2, on flame stability and ignition. Flame stability is an important issue in designing burners because it affects the lifespan of the burner as well as pollutant formation and heat fluxes in the furnace. Experimental data on oxy-coal flame ignition that would meet the criteria for LES simulations are very limited. Zhang [112, 113] studied how composition of the coal transport medium (primary CO2 and O2 ) affects the observed ignition stand-off distance in coaxial turbulent diffusion, oxy-coal flames. Results show that flame stability is affected by primary PO2 , wall temperature, secondary preheat temperature, secondary PO2 and transport medium. The burner and the furnace geometry were chosen to be simple enough for LES simulations. It was therefore proposed to perform simulations of Zhang's cases and to compare predictions on the flame stand-off distance with the experimental data. To further assess the capabilities of LES coupled with DQMOM to predict coal flame ignition, another validation study was also performed for an air-fired pulverized coal flame. LES was applied to a coal jet flame ignited by a preheated gas flow. The simulation results were compared qualitatively and quantitatively to experimental observations for different inlet stoichiometric ratios. The details for each simulation performed in this thesis can be found in the input files available online at: https://software.crsim.utah.edu/svn/LES_Coal/Pedel. CHAPTER 2 LARGE EDDY SIMULATION OF COAL PARTICLES IN TURBULENT COAXIAL JETS USING THE DIRECT QUADRATURE METHOD OF MOMENTS Julien Pedel ; a, Jeremy N. Thornocka, Sean T. Smitha, Philip J. Smitha Submitted to the Journal of Computational Physics Corresponding author. Email: jpedel@gmail.com aInstitute for Clean and Secure Energy, University of Utah, 155 South 1452 East, Room 350, Salt Lake City, UT, 84112, USA 2.1 Abstract Modeling coal particles in turbulent coal flames is a challenging but nonetheless critical step for predicting flame characteristics. Traditional approaches such as Reynolds-Averaged Navier-Stokes (RANS) coupled with Lagrangian solvers are unable to correctly predict par-ticle concentrations and velocities in complex systems. This study shows how Large Eddy Simulations (LES) coupled with an Eulerian solver can address those issues. Previous Direct Numerical Simulations and Stokes number analysis suggest that LES has the capabilities to resolve all the turbulent length scales affecting coal particles in a coal flame and could therefore lead to more accurate simulation. The effects of fluid aerodynamics on the particle 5 motion in the near field region of nonreacting coaxial particle-laden jets were simulated using LES. The Direct Quadrature Method of Moments (DQMOM) was used to track the particle phase in an Eulerian framework. Simulation results were compared to experimental data and accurately modeled a wide range of particle behaviors, such as particle jet waviness, spreading, break up, particle clustering and segregation, in different configurations. Simu-lations also predicted the mean axial velocity along the centerline for both the gas phase and the solid phase with a maximum error of 12% relative to experimental data. This study therefore provides a solid validation of the LES with DQMOM approach to model particles in turbulent flows and justifies its use for coal flame simulations. 2.2 Introduction Coal is an important and indispensable energy resource for electricity production as its reserves are far more abundant than those of other fossil fuels. As environmental regulations get stricter, coal thermal power plants need to control pollutants emissions such as NOx, SOx and ash particles. Coal also plays an important role in global CO2 emissions, which will need to be reduced to limit the effects of climate change [66]. In order to meet those requirements, understanding the pulverized-coal combustion mechanisms and developing advanced combustion technology is necessary. However, the combustion of pulverized-coal is a complex phenomenon compared with that of gaseous or liquid fuels, since dispersion of coal particles, devolitilization and oxidation reactions take place simultaneously [45, 92]. Moreover, experimental measurements of pulverized-coal flame characteristics are extremely difficult [38,78,80,88,97,100,107]. Accordingly, the development of new combustion furnaces and burners is still empirical and requires a significant amount of time and money. In coal-combustion systems, finely pulverized coal particles are conveyed by a medium, air in general, that enters the combustor as a jet. The distribution of particles and its dis-persion affect the local environment surrounding each particle and therefore its combustion kinetics [84]. Although the reaction kinetics involved in coal combustion have been stud-ied extensively, the fundamental mechanics of coal trajectories and of particle-laden jets in 6 general are as yet poorly understood and have not been included in optimizing the design and operating flow condition for coal combustors. This issue can be addressed by using Computational Fluid Dynamics (CFD), which was shown to be a very useful and accurate modeling framework for the description of particulate systems, and for the description of turbulent flows in general [24, 29]. RANS (Reynolds-Averaged Navier-Stokes), LES (Large-Eddy Simulation) and DNS (Di-rect Numerical Simulation) are the three common CFD methods for turbulence modeling. DNS requires a very fine grid spacing and even though it has the highest numerical ac-curacy, its application to practical combustion systems is made impracticable by the high computer load required. RANS on the opposite end of the spectrum is widely used for prac-tical applications. There is however evidence that RANS can be inaccurate in predicting particle-laden flows [5, 86]. As a result, interest for LES is growing, as it directly solves the transport equations for the large eddies and models only the smallest eddies. Unsteady turbulent motions are evaluated and the number of model parameters is greatly reduced. LES has been convincingly shown to be superior to RANS in accurately predicting turbu-lent mixing and combustion dynamics [56,75]. Although LES requires a high computer load compared to RANS, its superiority in terms of prediction will likely make it the best choice for practical combustion applications as computational cost decreases. In order to describe the evolution of the dispersed particulate phase, CFD must be coupled with a population-balance equation (PBE). Several numerical approaches can be used to solve this equation, such as classes methods (CM) or Monte-Carlo methods (MCM) [77]. Classes methods, in which the internal coordinate (e.g., particle length or volume) is discretized into a finite series of bins, are the most popular [40, 53, 54] , whereas MCM [11] are well known for their ease of implementation. In order to get reasonable results, the CM method requires a large number of classes (e.g., 20-30), so it is a computationally expensive approach for CFD calculations. Although the Monte-Carlo method is theoretically applicable, especially for Lagrangian-Eulerian simulations, in order to reduce the statistical error, a very large number of particles must be used. Due to limitations on the computational 7 resources, the number of particles which can possibly be tracked is often too low to get accurate particle statistics for applications such as coal boilers. The method of moments (MOM) offers an attractive alternative where the PBE is tracked through its moments by integrating the internal coordinate [43]. The main advantage of MOM is that the number of scalars required is very small (i.e., usually 4-6), which makes the implementation in CFD codes feasible. However, due to the difficulties related with expressing transport equations in terms of the moments themselves, the method has been scarcely applied until recently. The so-called closure problem [62, 94, 98, 103] has been reviewed [19] and a new method known as quadrature method of moments (QMOM) has been proposed [64] , validated [8, 61, 63] and applied for studying a wide range of practical problems [62,94,98,101,102] . The main limitation of the classical QMOM is that it can treat only PBE tracking one property of the population of particles, such as particle mass, size, volume or area (that is, monovariate PBE). QMOM has then been extended to bivariate problems [103,108,109]. However, in a number of practical cases, it is interesting to describe the population of particles with multivariate PBE, where two or more properties of the population are simultaneously tracked. In order to address these issues, the direct quadrature method of moments (DQMOM) has been formulated and validated [26, 59]. DQMOM is based on the direct solution of the transport equations for weights and abscissas of the quadrature approximation and presents the advantage of being directly applicable to multivariate PBE. It has been shown to be a powerful approach for describing polydispersed solids undergoing segregation, growth, aggregation and breakage processes in the context of CFD simulations [26, 32, 60]. Pulverized coal flames are considered dilute systems as the volume of particles is less than 1% of the volume of the fluid. Aggregation and breakage phenomena are rare events and are therefore usually neglected. However, even in this context where drag is the main process affecting particles, the DQMOM equations are equivalent to the widely used Lagrangian particle method [23], but it has the advantage of precisely controlling the statistical noise in the lower order moments (e.g., particle number density, mass density, Sauter radius). For a given desired accuracy, this greatly reduces 8 the computational cost since it does not require a large number of particles [17, 18]. In previous work, the advantages of using DQMOM for treating particle populations with low Stokes numbers (i.e., the dispersed-phase velocity follows closely the velocity of the continuous phase) have been clearly demonstrated [61, 63, 98, 115]. However, coal particles, which typically range from 1 micrometer to 200 micrometers in pulverized coal flames, have a finite Stokes number and their velocity depends on the particle size. In quadrature methods, this dependency is accounted for by solving an Eulerian model where each quadrature node has its own velocity field. Fox et al. [32] have applied DQMOM to droplets spray with a finite Stokes number and compared the results with the multifluid method and a classical Lagrangian solver. Results show that DQMOM can provide a comparable accuracy with a significantly lower computational cost, making it a very good candidate for more complex two-phase combustion applications. The goal of this study is to show that DQMOM can be used in combination with LES to predict particle behavior in coal flames. Since few data on particles in coal flames are available, experimental data of nonreacting turbulent flows with particles were chosen. The study presents LES simulations of coaxial particle-laden jets and compares the results to experimental data [14]. Special interest is given to the following important characteristics of particle-laden jets: gas and particle velocity, particle size segregation and particle clus-tering. These phenomena affect operational considerations within a coal-fired boiler such as flame stability, pollutants formation or heat flux characteristics [84]. Simulation of coaxial particle-laden jets are performed with three different velocity ratios between the primary and secondary stream and two particle size distributions (25 microns and 70 microns). Sim-ulation results show predictions of both the gaseous and solid phase averaged velocity and are compared to experimental data. 2.3 Numerical modeling Simulations were performed using ARCHES, a LES tool resulting from a ten-year part-nership with the Department of Energy and the University of Utah. It is a massively parallel 9 code that solves conserved quantities (mass, momentum, energy, scalar) spatially and tem-porally in a turbulent flow field, allowing for detailed and accurate simulations of fires and flames [87]. The code is integrated into a C++ framework called Uintah which provides large-scale parallelization tools for physics components [9,16,73]. The code is maintained in a repository and distribution is freely available [41]. The governing equations and numerical methods for the gas phase and the dispersed phase are described in the following sections. Interactions between the two phases are taken into account. 2.3.1 Governing equations of the gas phase Since the particles considered are nonreacting, the gas-phase governing equations are solved using filtered continuity and momentum equations for incompressible flows in a finite volume formulation. The mass balance can be written as: @ ug;i @xi = 0 (2.1) where ug;i is the ith gas velocity component, and the momentum equations as: @ ug;i @t + @ @xj ( ug;i ug;j) = 1 g @ p @xi + 2 @ @xj Sij @ ij @xj + Sui (2.2) where g is the gas density, is the gas kinematic viscosity and Sij is the rate-of-strain tensor. ij is the residual stress tensor and is modeled using the dynamic subgrid scale model of Germano et al. [36]. The source term Su is a momentum source term that accounts for momentum transfer from the dispersed phase. 2.3.2 Governing equations of the dispersed phase The DQMOM method was implemented into ARCHES to solve the dispersed phase. A detailed description of the particles representation and of the DQMOM equations is given in the following paragraphs. 10 2.3.2.1 Number density function (NDF) The particle Number Density Function (NDF) describes the number of particles per volume as a function of the spatial location and of other independent variables. These inde-pendent variables are intrinsic characteristics of the particles (e.g., particle diameter, particle velocities, etc.) and are called internal coordinates. The vector of internal-coordinate ran-dom values is denoted by , and the internal coordinate sample space is denoted by . When the particle velocities are considered as internal coordinates, the random values are denoted by up, and the particle velocity sample space is denoted by vp. The full NDF is denoted n(v; ; x; t). At a given location in space and time (x0; t0), the total number of particles per volume is: np(x0; t0) = vp n(vp; ; x0; t0)d dvp (2.3) NDFs can be univariate or multivariate. Univariate NDFs are only functions of one internal coordinate, so the internal coordinate sample space has a single dimension. Mul-tivariate NDFs, however, are functions of multiple internal coordinates, so the internal co-ordinate sample space has N dimensions: n( ; x; t) = n( 1; : : : ; N ; x; t) (2.4) or, including the particle velocity as an internal coordinate, n(vp; ; x; t) = n(vp;1; vp;2; vp;3; 1; : : : ; N 3; x; t) (2.5) The particle Probability Density Function (PDF) pu denotes the probability of the velocity-scalar vector taking on a particular value. At a fixed point in space and time (x0; t0), the PDF is related to the NDF by: 11 n(vp; ; x0; t0) = np(x0; t0)pu (vp; ; x0; t0) (2.6) 2.3.2.2 NDF transport equation The PDF transport equation is given by [76]: @ @t (pu (vp; ; x; t)) + @ @xi (hvp;ijvp; i pu (vp; ; x; t)) = @ @vi (hAijvp; i pu (vp; ; x; t)) @ @ k (hGkjvp; i pu (vp; ; x; t)) (2.7) In this equation, the variable hvp;ijvp; i represents the entire velocity sample space and is thus a full distribution. The quantities hAijvp; i and hGijvp; i are conditional quantities that describe the \"velocity\" of the PDF in the phase space (vp; ). Ai is defined by: Ai = dvp;i dt (2.8) As the right side depends on vp and , this quantity is a distribution. The particular value of Ai depending on v and can be expressed as hAijvp; i. Likewise, Gk is defined by: Gk = d k dt (2.9) and is also a distribution, with the particular value of Gk expressed as a conditional quantity hGkjvp; i. Eq. (2.7) can be multiplied by the function np(x; t) and combined with the number balance equation to yield the NDF transport equation [77]: @ @t (n(vp; ; x; t)) + @ @xi (hvp;ijvp; i n(vp; ; x; t)) = @ @vi (hAijvp; i n(vp; ; x; t)) @ @ k (hGkjvp; i n(vp; ; x; t)) +h(vp; ; x; t) (2.10) 12 where h is a source term representing the birth and death of particles in the domain. This term is set to zero as particle breakage and aggregation effects are neglected. 2.3.2.3 Direct Quadrature Method of Moments The Direct Quadrature Method of Moments (DQMOM) approximates the NDF by a summation of multidimensional Dirac delta functions. Figure 2.1 gives an example for bivariate Gaussian distribution (i.e., two internal coordinates) of the NDF and its quadrature approximation with three quadrature nodes. The quadrature approximation can be regarded as an extension of the mono-kinetic assumption for the multifluid model. In the multifluid model, the closure of the system is obtained through the following assumptions [57]: Fig. 2.1: DQMOM approximation of the NDF with 3 quadrature nodes 13 - For a given particle size, at a given point (x,t), there is only one characteristic averaged velocity, - The velocity dispersion around the averaged velocity is zero in each direction at every point. This is equivalent to presuming that the NDF is described by: n(vp; ; x; t) = np( ; x; t) (vp hvpi) (2.11) The quadrature approximation extends the mono-kinetic assumption to all the internal coordinates: - For a given quadrature node , at a given point (x,t), there is only one characteristic averaged internal coordinate value h i , also called abscissa, - The internal coordinate dispersion around the averaged value h i is zero. With those assumptions, the NDF is written as : n(vp; ; x; t) = XN =1 w ( h i ) (vp hvpi (2.12) n(vp; ; x; t) = XN =1 0 @w NY 3 i=1 ( i h ii ) Y3 i=1 vp;i hvp;ii 1 A (2.13) where w is the number of particles per volume associated with the quadrature node . w , h i and hvpi depend on space and time, but the dependence is omitted for clarity of notation. The validity of the mono-kinetic assumption has been discussed for the multifluid method [17]; the hypothesis is valid for particles with low Stokes number and becomes inappropriate when Stokes number increases (St > 0.4) and particle trajectory crossing oc-curs. It is however important to note that this limitation can be partially overcome with DQMOM. In the multifluid model, only the particle size space is discretized, whereas DQ- 14 MOM allows us to discretize a multidimensional space. It is therefore possible to discretize the particle velocity space if the particle velocity is chosen as an internal coordinate and then to have several quadrature nodes with different velocities and otherwise identical internal coordinates values. As quadrature nodes can cross each other, the error associated with finite Stokes number particles crossing each others is reduced. This strategy was applied for large particle-size distribution and its validity is discussed in the second part of this work. 2.3.2.4 Quadrature-approximated NDF transport equation The quadrature-approximated NDF transport equation can be derived for a multivariate NDF by combining the NDF transport equation (Eq. (2.10)) with the quadrature approxi-mation (Eq. (2.13)). The quadrature node averaged phase space velocities hAii and hGii are introduced: hAii = Aijvp = hvpi ; = h i (2.14) hGii = Gijvp = hvpi ; = h i (2.15) The multivariate quadrature-approximated NDF transport equation is then given by a set of equations: @ @t (w ) + @ @xi hvp;ii w = a (2.16) @ @t (&n ) + @ @xi hvp;ii &n = bn (2.17) where &n = w h ni is the weighted abscissa for the nth internal coordinate. The terms on the right-hand side are determined by the following equations: 15 XN =1 2 4 YN j=1 j h ji + XN m=1 @ @ h mi 0 @ YN j=1 j h ji 1 Ah mi 3 5a XN =1 XN n=1 2 4 @ @ h ni 0 @ YN j=1 j h ji 1 A 3 5bn = S (2.18) where S is the sum of the phase space convective terms: S = XN =1 @ @vi 2 4hAii w 0 @ Y j vp;j hvp;ji 1 A 0 @ Y j j h ji 1 A 3 5 XN =1 @ @ i 2 4hGii w 0 @ Y j vp;j hvp;ji 1 A 0 @ Y j j h ji 1 A 3 5 (2.19) 2.3.2.5 Moment-transformed NDF transport equation The quadrature-approximated NDF transport equation (Eq. (2.18)) is a single equation but multiple unknowns need to be determined (weights and abscissas). In order to obtain a number of independent equations equal to the number of weights and abscissas, a set of independent moments are chosen and the moment transform of the quadrature-approximated NDF transport equation yields a set of independent equations, equal in number to the independent moments. A number of moments equal to (N + 1)N is required, where N is the number of internal coordinates and N the number of environments. This procedure results in the moment-transformed quadrature-approximated NDF transport equation [60]: XN =1 2 4 0 @ YN j=1 D kj j E 1 A XN m=1 km D km m E 0 @ YN j=1;j6=m D kj j E 1 A 3 5a + XN =1 XN n=1 2 4kn D kn1 n E 0 @ YN j=1;j6=n D kj j E 1 A 3 5bn = Sk (2.20) where: 16 Sk = XN =1 XN n=1 0 @ YN j=1;j6=n D kj j E 1 Akn D kn1 n E w hGni (2.21) This set of equations forms a linear system that can be solved to get the source terms of the weights and weighted abscissas transport equations (Eqs. (2.16) and (2.17)) a and bn . 2.3.2.6 Particle velocity model In this work, the particles considered are nonreacting glass beads. The aggregation/breakage phenomena are neglected and the only physical process modeled is the drag force, which can be described by the Stokes drag law. For a mesoscale size particle, other forces such as lift can be neglected and the momentum equation for the particle can be expressed by an ordinary differential equation. dup;i dt = fdrag p (ug;i up;i) + gi ( p g) p (2.22) where g is the gravity force acting on the particle and fdrag is the coefficient of the drag force, which has a close relationship with the particle Reynolds number [85]: fdrag = 8>>>>< >>>>: 1 Rep < 1 1 + 0:15Re0:687 p 1 < Rep < 1000 0:0183Rep Rep > 1000 (2.23) Rep = pdp jup ugj g (2.24) where g is the gas dynamic viscosity. In Eq. (2.22), p is the particle relaxation time: 17 p = pd2 p 18 g (2.25) The particle density p and particle diameter dp are assumed to be constant. This model can easily be implemented into the DQMOM formulation through the term: hAii = fdrag p ui;g hvp;ii + gi ( p g) p (2.26) The gas phase and the particle phase are then coupled through the source term Sui in the gas momentum equation (Eq. (2.2)): Sui = 1 g XN =1 w p fdrag p ui;g hvp;ii (2.27) 2.3.3 Stokes number analysis The major advantage of LES over RANS is that a large part of the turbulent kinetic energy spectrum (usually about 80%) is resolved as opposed to RANS where it is modeled. However, a subgrid scale model (SGS) is still needed to model the effects of the smallest eddies on the flowfield. The unclosed subgrid scale terms appear when the spatial filter is applied to the Navier-Stokes equations (Eq. (2.2)). A similar issue arises when filtering the NDF transport equation (Eq. (2.10)): @ @t (n(vp; ; x; t)) + @ @xi ( vp;in(vp; ; x; t)) = @ @vi hAijvp; in(vp; ; x; t) @ @ k hGkjvp; in(vp; ; x; t) +@ sgs;i @xi + @ sgs;ui @ui + @ sgs; k @ k (2.28) The subgrid scalar flux sgs;i represents the flux of the number density as a result of unresolved velocity fluctuations. Likewise, the subgrid scalar fluxes sgs;ui and sgs; k represent the subgrid flux of the number density in phase space (vp; ). 18 Few studies have focused on the effects of subgrid scale velocity fluctuations on particle motion in LES so far. Yei and Lei [106] used LES to investigate the motion of particles in isotropic and homogeneous shear flows. They generated particle statistics by neglecting the subgrid scale (SGS) effects on particles and showed that LES can successfully predict second-order statistics of the particle motion, such as root-mean-square velocity fluctua-tions, the mean square displacement and the dispersion coefficient. They also showed that particle dispersion was dominated by the large eddies. Later, Yand and Lei [105] studied the phenomenon of particle accumulation in low vorticity regions and found that it is controlled mainly by the small scales of turbulence with wave number k! corresponding to the max-imum of the dissipation (vorticity) spectrum. The small eddies with wave number greater than 2:5k! had no effect on the particle accumulation and the average settling velocity. Since the length scale l! = 1=k! is typically one order of magnitude greater than the Kolmogorov length scale , they concluded that LES was adequate for particle-laden flows and that the subgrid-scale effects can be neglected, provided that the cutoff wave number is greater than 2:5=l!. As a result, many LES studies of particle-laden flows neglected the SGS effects for various flow configurations [5, 10, 82, 99]. More recently, Shotorban and Mashayek [81] used a stochastic model for particle motion in LES and found that this assumption may break down when particles have small time constant and/or the filtered energy is significant. The error was however important for = 8 but very limited for = 4 , which tends to agree with Yang and Lei's results. The previous studies only considered a limited range of turbulence intensity or particle size. A general theoretical justification for direct use of the drag law in LES can be made by a Stokes number analysis. The particle Stokes number combines the effects of particle size, particle density and turbulence intensity, and simultaneously accounts for the variety of turbulent timescales. The Stokes number is defined as: St = p f (2.29) 19 where f is a characteristic time of the fluid. The difficulty of Stokes number analysis often lies in determining f . The turbulent kinetic energy (TKE) spectrum, denoted as EU, is convenient for this analysis since it breaks down the turbulent fluctuations into contributions of different sizes - with each size exhibiting a corresponding eddy turnover time. When the eddy turnover time is used as the fluid characteristic time in the Stokes number, and when each wave number in the turbulent-energy spectrum has an associated eddy turnover time, then a single particle of fixed radius and density will exhibit a unique Stokes number value for each different wave number. To illustrate this concept, the model energy spectrum for isotropic turbulence is shown in Figure 2.2a. This figure shows a dotted line corresponding to the wave number for which the eddy turnover time is equal to a particle's relaxation time - Stokes number equals one for that particle. For wave numbers much larger than the dashed line, the eddy turnover time is much smaller than the particle relaxation time and the particle behaves ballistically (with respect to the fluctuations at those wave numbers). On the other hand, for wave numbers much smaller than the dashed line, the eddy turnover time is much larger than the particle relaxation time and the particle behaves as a tracer element. It is only for a relatively narrow bandwidth that the turbulence interacts nontrivially with the particle. The implications of this concept on LES modeling are illustrated in Figure 2.2b. In these three plots, the model energy spectrum is shown with a hypothetical Nyquist cut-off. Three particle sizes are considered, one for each of the three plots, with the overlap bandwidth of nontrivial Stokes number illustrated in relation to the Nyquist cut off. The small particles, top plot, behave as tracers for all of the resolved scales. Therefore, the drag of particles in this unresolved regime should be modeled in the subgrid-scale flux. In the extreme of zero Stokes number, the subgrid-scale flux is often modeled using gradient diffusion. The large particles, bottom plot, are in a regime such that they are ballistic with regard to all the unresolved scales - their trajectories are essentially unaffected by the presence or absence of the unresolved scales. For this reason, particles in the resolved regime should be modeled using the drag directly - as done in direct numerical simulations. The 20 (a) (b) Fig. 2.2: Illustration of the different particle kinematic regimes as a function wave number for a fixed particle relaxation time (a) and the implications on LES modeling (b) 21 middle plot shows particles of intermediate size for which some of the nontrivial drag effects are resolved and some are not. Particles in this partial-resolution-of-drag regime present a challenging modeling problem. It will be shown in the results section, 2.5.2, that particles in the current study are in the resolved regime. Therefore, the drag law is used directly and the subgrid scale effects of turbulence on particles were neglected 2.4 Experimental and simulation conditions In this study, LES coupled with DQMOM is applied to simulate particle-laden coaxial jets and results are compared to experimental data obtained by Budilarto [14]. 2.4.1 Experimental setup In the experiments, particles were conveyed with air in the primary nozzle whereas the secondary nozzle was fed only with air (Figure 2.3). Air and particles were injected at room temperature in an 18x18 inch chamber. The experiments focused on the effects of the inlet velocity and the particle sizes on the flow-field and the particle motion in the near field region of the nozzle. The velocity ratio, defined as the ratio of the annular to central inlet velocity (V R = Ua=U0), was modified by varying the volumetric flow rate of the annular jet. Three different velocity ratios were studied: 0, 1.0 and 1.5. The Reynolds number based on the average central inlet velocity was 8,400. Table 2.1 lists the flow conditions for the experiments. Two particle-size distributions were used for this purpose, one with a mean diameter of 25 microns and the other with a mean diameter of 70 microns. Particles used were glass beads with a density of 2500 kg=m3. The particle mass loading was maintained at 0.5 for all velocity ratios. Laser Doppler Velocimetry and Phase Doppler Particle Analyzers were used to measure time-averaged velocities. Those experimental data were chosen for several reasons. First, measuring coal particles in coal flames is challenging and very few studies have been reported. Even though our final goal is to simulate coal flames, nonreacting systems offer a more reliable and accurate source of data for particle velocities, segregation and clustering effects. 22 Fig. 2.3: Schematic of the coaxial nozzle Table 2.1: Experimental setup for the different velocity ratios Case d (inch) D (inch) Central jet maximum inlet velocity U0 (m/s) Annular jet maximum inlet velocity Ua (m/s) Velocity ratio (VR) 1 0.56 1.254 11.7 0 0 2 0.56 1.254 11.7 11.6 1.0 3 0.56 1.254 11.7 18.0 1.5 Budilarto's data include extensive results in different operating conditions (3 velocity ratios) with two different particle sizes distributions, which allows us to test the model sensitivity and accuracy with respect to turbulence intensity and particle size. Then, the experiments offer a simple geometry with well-defined boundary conditions and are therefore suited for numerical simulations. Finally, the operating conditions are very close to those of pulverized coal flames in terms of velocities, turbulence intensity, particle mass loading, 23 particle sizes and volume fraction. In addition, coaxial burners are also used in coal flame experiments on a laboratory scale [72, 112, 113], so this study provides a first step towards the simulation of coaxial coal flames. 2.4.2 Simulation details The computational domain is a 0.28 m wide cube which is divided with a cartesian mesh into 8 millions cells, resulting in a 1.42 mm cellsize. A total of 9 simulations were performed; for each velocity ratio, simulations without particles, with 25-micron particles and with 70-micron particles were run on 512 CPUs for about 24 hours. Each simulation was performed for at least 15,000 timesteps, each timestep being about 25 s and at least 5,000 steps were used for time averaging. Measurements have shown that the common power law velocity profile of fully developed turbulent flow can be applied at the central inlet: U U0 = 1 2r d 1=7 (2.30) The dispersed phase was modeled using the DQMOM method with four internal co-ordinates: the particle diameter dp and the three particle velocity components vx, vy, vz. Particles were assumed to be perfect spheres of constant size and density. Special care must be taken when choosing the quadrature nodes to represent the full NDF. As discussed earlier, for small particles, the mono-kinetic assumption is applicable and all the particles of a same size can be assumed to have a unique velocity at a given location. On the other hand, larger particles have more inertia and their trajectory history may result in multiple velocities for same size particles at a given location. In this case, the mono-kinetic assumption would not be valid anymore. It is however possible to overcome this limitation with the DQMOM approach by discretizing the velocity phase space into several quadrature nodes. Taking those considerations into account, a different modeling approach is used to represent the 25-micron particle size distribution and the 70-micron distribution. 24 For the 25-micron distribution, the particle diameters range from a few microns up to 45 microns. The mass of the smallest particles is several orders of magnitude less than the mass of the largest particles. Since a 10-micron particle would behave significantly differently from a 40-micron particle, the size distribution was represented by three quadrature nodes in the DQMOM approach, each node corresponding to a different particle size. Table 2.2 shows the particle diameter associated with each quadrature node as well as their weight at the inlet. Even though the mean diameter is 25 microns, it can be noticed that there are more 14-micron particles than 25-micron particles. For the 70-micron distribution, the particle diameters range from 50 to 90 microns. The mass of all the particles is of the same order of magnitude, so the size distribution can be represented by a single particle size of 70 microns. However, these particles have more inertia than the small ones and, as they collide and bounce against the walls inside the nozzle, their initial velocity at the exit of the nozzle follows a distribution that should be taken into account. For the 25-micron particles, measured velocity fluctuations at the inlet were indeed small, so assuming that all the particles have the same inlet velocity is a good approximation. On the contrary, for the 70-micron particles, velocity fluctuations at the inlet were important and it is necessary to represent their velocity distribution. Based on the measured root mean square of velocity fluctuations at the inlet, a 2-D normal distribution in particle velocity components vp;y and vp;z was assumed. Fluctuations in vp;x were neglected since they are small compared to the mean velocity in this direction. The normal distribution was approximated by 5 quadrature nodes in the DQMOM approach, as shown in Table 2.3. Table 2.2: Quadrature approximation for 25-micron particles (inlet values) Quadrature node Particle diameter ( m) Weight (# particles/m3) 0 14 4.11E10 1 25 1.44E10 2 34 2.90E9 25 Table 2.3: Quadrature approximation for 70-micron particles (inlet values) Quadrature node Particle diameter ( m) vp;x (m/s) vp;y (m/s) vp;z (m/s) Weight (# particles/m3) 0 70 9.6 0 0 3.8E8 1 70 9.6 -0.55 0 2.3E8 2 70 9.6 0.55 0 2.3E8 3 70 9.6 0 -0.55 2.3E8 4 70 9.6 0 0.55 2.3E8 2.5 Results and discussion Simulations were performed for single-phase coaxial jets and for particle-laden coaxial jets to assess the uncertainties associated with turbulence modeling and particle modeling. 2.5.1 Single phase coaxial jet Simulations were performed for the three velocity-ratios cases without particles. The two main objectives are first to validate the modeling approach by ensuring that the gas mean velocities are predicted correctly and secondly to have a reference case to then understand the effect of particles on the gas dynamics. Figures 2.4, 2.5 and 2.6 show instantaneous profiles of the x-component of the gas velocity. For VR=0, the flow field is identical to a round jet, whereas for VR=1.0 and VR=1.5, the effect of the secondary stream on the flowfield can be observed. Simulations results were evaluated against Budilarto's experimental data by comparing the mean gas velocity along the centerline (Figure 2.7). Velocity profiles are predicted with a maximum error less than 10%. This error can be attributed to several factors, among which grid resolution, turbulence modeling or experimental measurements uncertainty. Results however predict the velocity profile and jet spreading accurately enough to validate the modeling approach. 26 Fig. 2.4: Gas x-velocity profile (m/s) for VR=0, no particles 27 Fig. 2.5: Gas x-velocity profile (m/s) for VR=1.0, no particles 28 Fig. 2.6: Gas x-velocity profile (m/s) for VR=1.5, no particles 29 0 2 4 6 8 10 12 14 16 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 x/d U g /U 0 VR = 0, LES VR = 0, Exp. VR = 1.0, LES VR = 1.0, Exp. VR = 1.5, LES VR = 1.5, Exp. Fig. 2.7: Comparison between LES and experimental data of gas averaged x-velocity along the centerline for different velocity ratios / No particles 2.5.2 Stoke number analysis Results from the single phase coaxial jet can be used to perform a Stokes number analysis a posteriori. For this analysis, we assume local isotropy of turbulence. The turbulent kinetic energy spectrum can be modeled by a simple model spectrum [76]: E( ) = C\"2=3 5=3fL( L)f ( ) (2.31) 30 where \" is the dissipation rate, is the wave number, L is the characteristic length scale of the large eddies and is the Kolmogorov length scale. fL and f are specified nondimensional functions: fL( L) = 0 @h L ( L)2 + cL i 1 A 5=3+p0 (2.32) f ( ) = exp h ( )4 + c4 i1=4 c (2.33) The constant values of C, p0 and are taken to be C = 1:5, p0 = 2, = 5:2. For high Reynolds number, the values of the constants cL and c can be approximated by cL 6:78 and c 0:40 [76]. For specified values of \", L and , the model spectrum is determined by Eqs. (2.31), (2.32) and (2.33). The dissipation rate can be estimated from the LES simulations by the simplified for-mula: \" = 15 @ug;1 @x1 2 (2.34) For the LES simulation of the single phase coaxial jet with the highest turbulence intensity, i.e., VR=1.5, the maximum of the dissipation rate along the centerline is calculated to be \" 40m2=s3. The Kolmogorov length scale can be estimated from this value: 3=\" 1=4 0:1mm (2.35) The characteristic length scale of the large eddies is obtained by: L =Re3=4 0:09m (2.36) 31 The resulting normalized model energy spectrum EU = E=( u2 ) is plotted on Figure 2.8a, where u (\" )1=4 is the Kolmogorov velocity. The relationship between the dissipa-tion energy and the energy spectrum (Eq. (2.37)) can be used to verify that the estimated values for , L and are consistent. \" = 1 0 2 2E( )d (2.37) The value of the integral in Eq. (2.37) returns \" = 39:4 J=s, which is close to the LES simulation value \" 40 J=s. The wave number = 1= corresponding to the Nyquist cut-off limit is indicated by a green line in Figure 2.8, where is the LES filter length ( = 1:4mm). The model spectrum indicates that more than 90% of the energy spectrum is resolved for the LES simulations performed in this study and can further be used to determine the eddy characteristic time f as a function of the wave number by assuming: f ( ) 1 [2k( )]0:5 (2.38) where k( ) E( ) is the turbulent kinetic energy associated with eddy of wave number . Eq. (2.38) combined with the Stokes number definition (Eq. (2.29)) allows us to plot the Stokes number as a function of the wave number for 25-micron and 70-micron particles, as shown on Figures 2.8b and 2.9b, respectively, ( p = 2500 kg=m3). The black dotted lines represent the wave number for which the particle Stokes number is equal to one. The top images of Figures 2.8 and 2.9 show the respective locations of the Nyquist limit and the unity Stokes number wavelength for the 25-micron particles (Figure 2.8a) and the 70-micron particles (Figure 2.9a). In both cases, the Nyquist cut-off wave number is higher than the unity Stokes wave number by a safe margin. It can be concluded that both particle sizes are in the resolved regime and that the specified mesh resolution for the LES simulations in this study allows us to resolve all the eddies affecting the particles motion. 32 (a) (b) Fig. 2.8: Location of the Nyquist cut-off and the unity Stokes wave number on the model energy spectrum for 25-micron particles (a) and Stokes number as a function of wave number for 25-micron particles (b) 33 (a) (b) Fig. 2.9: Location of the Nyquist cut-off and the unity Stokes wave number on the model energy spectrum for 70-micron particles (a) and Stokes number as a function of wave number for 70-micron particles (b) 34 2.5.3 Two-phase coaxial jet Simulation results of the coaxial jet for the three velocity ratios and the two particle size distributions are presented in this section. Results for the six simulations are analyzed both qualitatively and quantitatively. In general, the qualitative analysis of the particle behavior predicted by the simulations shows very good agreement with the experimental observations. Figures 2.10 and 2.11 show instantaneous snapshots of particle concentration for a given quadrature node in the XZ plane for the simple jet case (VR=0). Images for the 25-micron particles case are presented in Figure 2.10, whereas images for the 70-micron particles case are shown in Figure 2.11. For the 25-micron particle size distribution, particle segregation is observed: the smaller particles (quadrature node 0) disperse faster than the larger ones (quadrature node 2). The 14-micron particles follow closely the gas flowfield and form a wider jet than the 34-micron particles which remain mostly along the centerline. The 70- micron particles have more inertia and tend to keep their initial velocity. Figures 2.11a, 2.11b and 2.11c, respectively, show particle concentration for quadrature nodes 0, 3 and 4. Quadrature nodes 1 and 2 are not shown but would be similar to quadrature nodes 3 and 4 in the XY plane. For each quadrature node, little dispersion is seen around the main jet as particles are less affected by turbulence. The effect of their initial velocity is more prominent and is felt on their trajectory until they reach a semi-equilibrium with the gas phase. Figure 2.12 presents time-averaged particle concentrations at different locations for VR = 0. wsum is the sum of the weights of each quadrature node and represents the total number of particles and is defined as: wsum = XN =1 w (2.39) Figure 2.12a shows particle concentration in the radial direction at x = 15d for 25- micron particles. This plot confirms that particle dispersion decreases with particle size and that the dispersion of the three particle sizes is significantly different. 35 (a) Weight for quadrature node 0 (b) Weight for quadrature node 1 (c) Weight for quadrature node 2 Fig. 2.10: Instantaneous profiles of weights (# particle/m3) for 25-micron particles, VR=0 36 (a) Weight for quadrature node 0 (b) Weight for quadrature node 3 (c) Weight for quadrature node 4 Fig. 2.11: Instantaneous profiles of weights (# particle/m3) for 70-micron particles, VR=0 37 −2 −1 0 1 2 0 1 2 3 4 5 6 7 x 10 10 y/d W (# particles/m3) W sum W 0 W 1 W 2 (a) Number density at X = 15d, 25-micron particles −2 −1 0 1 2 0 2 4 6 8 10 x 10 8 y/d W (# particles/m3) W sum W 0 W 1 W 2 W 3 W 4 (b) Number density at X = 15d, 70-micron particles 0 5 10 15 0 2 4 6 8 10 12 x 10 10 x/d W (# particles/m3) W sum W 0 W 1 W 2 (c) Number density along centerline, 25-micron particles 0 5 10 15 0 2 4 6 8 10 12 14 x 10 8 x/d W (# particles/m3) W sum W 0 W 1 W 2 W 3 W 4 (d) Number density along centerline, 70-micron particles Fig. 2.12: Time-averaged particle concentration at X = 15d [(a), (b)] and along the centerline [(c), (d)] for VR = 0. 38 Figure 2.12b presents a similar plot for the 70-micron particles case. The particles with no initial radial velocity (w0) show very little dispersion, whereas the particles with an initial y-velocity (w1, w2) show a significant dispersion along the y-axis. As a result, the overall profile of wsum is wider than the profile of w0 and it can be concluded that particle dispersion for 70-micron particles occurs because of their initial velocity whereas for the 25- micron particles, it is due to gas phase interactions. Figures 2.12c and 2.12d show averaged particle concentration profiles along the centerline. For the 25-micron particles, a peak is observed before the particle concentration decreases. This peak is located further down the jet as the particle size increases. As a matter of fact, the peak for w0 in the 70-micron particles case is not shown on the figure but would be located at about x = 20d. To validate the LES with DQMOM approach to model particles-laden jets, we first compare qualitatively the particle behavior predicted by the simulations in the different configurations to that observed experimentally. Figure 2.13 presents instantaneous screen-shots of the total particle concentration wsum from LES simulations and photographs of particles taken during the Budilarto's experiment for the three velocity ratios cases with the 25-micron particle size distribution. For the velocity ratio of 0, Budilarto noted that the region with the highest particle number density was located along the centerline and its shape resemble a cone (Figure 2.13b). The 25-micron particles were dispersed symmetrically and gradually in the radial direction as the axial distance increased. He also observed that the conical region did not have a smooth interface, but its shape followed a \"pine-tree\" like structure where wedges of particles with backward slope were present. This effect resulted from the partial interaction of the particles with the large eddies: particles were thrown outwards at the leading edge of the large eddies, but the large eddies did not have enough energy to bring the particles into the high vorticity region, i.e., the center of the large ed-dies. As can be seen on Figure 2.13a, the LES simulation also predicts the conical and symmetrical shape and the \"pine-tree\" structures are present as well. For VR=1, Budilarto reported that for x<8d, there was no change in the cross-sectional area of the structures of the 25-micron particles and that the particle number density was approximately constant 39 (a) VR = 0 (b) VR = 0 (c) VR = 1.0 (d) VR = 1.0 (e) VR = 1.5 (f) VR = 1.5 Fig. 2.13: Qualitative comparison of particle behavior for the 25-microns particles cases, LES screenshots of particle concentration on the left (Wsum), photographs from experiments on the right. in this region (Figure 2.13d). The presence of the annular jet inhibited the spreading of the particles in the radial direction. At x>8d, the jet of particles started to become wavy as the turbulent structures from the outer shear layer interact with the particles. The LES simulation predicts the same behavior (Figure 2.13c); the jet of particles remains compact with a constant particle concentration until x=8d and then it becomes wavy. For VR=1.5, Budilarto pointed out that the particles jet started becoming wavy at the end of the po-tential core, which was located at x=4d (Figure 2.13f). This waviness was attributed to interactions between particles and the large eddies of the inner shear layer. The radial dis- 40 persion of particles was significantly enhanced at x>6d, as particles interacted with eddies of the length scale on the order of the jet diameter. Due to the presence of more energetic turbulent structures, the radial spreading of the particles in VR=1.5 was higher than that in VR=1.0. The more energetic eddies were able to drag the particles to follow their motion and to create large displacement in the radial direction, which caused the particle jet to break and become asymmetrical. For VR=0, the large eddies only had enough energy to displace particles near the jet edges whereas for VR=1.5, the large eddies are able to displace particles with no preference if they are at the center of the jet or near the edges. All those observations are in agreement with the LES simulation: the particles jet becomes wavy at x=4d and the waves grow at x>6d (Figure 2.13e). The jet eventually breaks up into clusters of particles as observed experimentally. Figure 2.14 presents the same results for the three cases with the 70-micron particle size distribution. For VR=0, Budilarto reported that the particles flowed downstream with a straight-line trajectory and were not affected by the presence of the turbulent structures in the flow (Figure 2.14b). The 70-micron particles spread in the radial direction and the spreading was attributed to the presence of radial velocity fluctuations at the nozzle inlet. The LES simulation results (Figure 2.14a) confirm those statements; the particle jet is little affected by turbulence and the spreading results mainly from the initial particle velocity, as shown previously in Figure 2.10. For VR=1.0, particles flowed with straight-line trajectories until x=9d. At x>9d, the turbulent structures formed in the outer shear layer began to affect the particles motion and small waves in the particles flow appeared (Figure 2.14d). The waviness was less pronounced than for the 25-micron particles case because of the particle inertia and did not increase the radial spreading. As can be seen on Figure 2.14c, LES results agree with the observations since small waves appear after x=9d and the jet spreading is limited. Finally, for VR=1.5, the particles structure became wavy at x>7d-8d because of the higher turbulence intensity (Figure 2.14f). A significant increase in the radial spreading was observed, exceeding other velocity ratios, as the particles were dragged and dispersed by the large eddies. The enhancement in the radial spreading for 41 (a) VR = 0 (b) VR = 0 (c) VR = 1.0 (d) VR = 1.0 (e) VR = 1.5 (f) VR = 1.5 Fig. 2.14: Qualitative comparison of particle behavior for the 70-micron particles cases, LES screenshots of particle concentration on the left (Wsum), photographs from experiments on the right. 42 the 70-micron particles was less important than that for the 25-micron particles. The LES simulation (Figure 2.14e) also predicts that the particles jet becomes wavy at x=7d and that the radial dispersion is greater compared to other velocity ratios but less important compared to the 25-micron particles case. Overall, simulations were able to capture very different particle behaviors, such as jet spreading, waviness, break up or particle clustering, in the six configurations studied. To further validate the LES coupled with DQMOM approach to model turbulent particle-laden flows, we now compare the time-averaged velocities along the centerline pre-dicted by the simulations to the experimental measurements. Figure 2.15 presents plots of the measured and predicted mean gas and particles velocities. For the simulation results, the mean particle velocity was calculated by averaging the first moment of the particle velocity v1 i over time, where: v1 i = PN =1 w hvii P N =1 w (2.40) Since the inertia of the 25-micron particles is smaller than that of the 70-micron parti-cles, they respond faster to the variations in the gas motion. This is clearly shown by the similarity in the centerline development of the axial mean velocity between the gas and the particles for the smaller particles (Figures 2.15a, 2.15c and 2.15e). In general, the develop-ment of the particles velocity along the centerline includes a first region where the gas-phase velocity is faster than the solid-phase velocity and, as the gas phase velocity decreases and the particles retain more momentum, the two curves cross each other and the particle ve-locity becomes higher than the gas velocity. The length of this region for the 25-micron particles is short because of their small slip velocity at the nozzle inlet. Since the 70-micron particles have a velocity at the nozzle exit that is about 20% lower than the gas velocity, the length of this region is greater, about 7d for VR=0 and 1.0 and is increased to 13d for VR=1.5 (Figures 2.15b, 2.15d and 2.15f). A comparison with the single phase results (Figure 2.7) also shows that the particles have a strong influence on the gas phase velocity. 43 0 5 10 15 5 6 7 8 9 10 11 12 x/d U c (m/s) Gas velocity, LES Gas Velocity, Exp. Particle velocity, LES Particle velocity, Exp. (a) VR = 0 0 5 10 15 5 6 7 8 9 10 11 12 x/d U c (m/s) Gas velocity, LES Gas velocity, Exp. Particle velocity, LES Particle velocity, Exp. (b) VR = 0 0 5 10 15 5 6 7 8 9 10 11 12 x/d U c (m/s) Gas velocity, LES Gas velocity, Exp. Particle velocity, Exp. Particle velocity, Exp. (c) VR = 1.0 0 5 10 15 5 6 7 8 9 10 11 x/d U c (m/s) Gas velocity, LES Gas velocity, Exp. Particle velocity, LES Particle velocity, Exp. (d) VR = 1.0 0 5 10 15 5 6 7 8 9 10 11 12 x/d U c (m/s) Gas velocity, LES Gas velocity, Exp. Particle velocity, LES Particle velocity, Exp. (e) VR = 1.5 (f) VR = 1.5 Fig. 2.15: Comparison between simulations and experimental data of time-averaged gas and particles velocities along the centerline for 25 m particles [(a) (c) (e)] and 70 m particles [(b) (d) (f)] 44 Particles indeed transfer some of their momentum to the gas phase and thus the gas velocity along the centerline decreases slower. For instance, for VR=0, the gas velocity magnitude at x=15d is less than 40% of its initial magnitude at the nozzle exit whereas it is more than 50% for both cases with particles. The LES simulations results are in good agreement with the experimental data; the maximum error over all the datapoints between the predicted velocities and the measured velocities for both phases is less than 12%. The uncertainty of the measured velocities can be assumed to be less than 5%, so the simulation predictions are slightly in disagreement with the experimental data. However, a comparison with Figure 2.7 shows that this disagreement can be attributed to the gas phase turbulence modeling since the maximum error for the single phase simulations was about 10%. The radial spreading of particles is an important effect in coaxial-laden jets but also in coal flames as smaller particles will disperse faster and ignite first. Figure 2.16 presents the time-averaged concentration of particles wsum predicted by LES simulations for the different velocity ratios along the radial direction at x = 15d. The predicted particle concentrations are compared to the particle data rate obtained from the LDV measurements by Budilarto. Values are normalized by the centerline value where the particle concentration is the highest. For every velocity ratio, the fine particles have a wider distribution whereas the coarse particle jet remains more compact. Overall, simulations give good predictions and show that the modeling approach is able to capture particle segregation, even though the spreading rate is over-predicted for VR = 0 and under-predicted for VR = 1.0 and VR = 1.5. The comparison with the experimental data is however difficult to interpret for several reasons. First, the LDV data rate is not an accurate measurement of the number of particles as the laser resolution is not high enough to distinguish two small particles from a large one. Then, the error due to the gas phase turbulence modeling may be responsible for the error in the prediction of the particle dispersion. Finally, the predicted distribution is highly sensitive to the initial choice of environments in the quadrature approximation, as shown in Figure 2.12, and a larger number of environments to represent the particle distribution in size and velocity would likely lead to more accurate results. For instance, for the 70-micron particles, 45 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 r/d W/W max 25 um, LES 25 um, Exp. 70 um, LES 70 um, Exp. (a) VR = 0 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 r/d W/W max 25 um, LES 25 um, Exp. 70 um, LES 70 um, Exp. (b) VR = 1.0 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 r/d W/W max 25 um, LES 25 um, Exp. 70 um, LES 70 um, Exp. (c) VR = 1.5 Fig. 2.16: Radial profiles of particles number density taken at x = 15d 46 the concentration of particles at r/d>1 tend to be under-predicted. This error could be reduced if the small fraction of particles with initial radial velocity greater than 0.55 m/s were taken into account. 2.6 Conclusion In this work, the Direct Quadrature Method of Moments was implemented in a mas-sively parallel LES code to model finite Stokes number particles in turbulent flows. The derivation of the DQMOM equations was shown for a multivariate NDF using the particle velocity components as internal coordinates. Equations were simplified by extending the mono-kinetic assumption used in the multifluid approach to all the internal coordinates and for each environment and by neglecting the subgrid-scale terms in the NDF filtered transport equations. This approach was tested by comparing simulation predictions to experimental results in nonreacting coaxial particle-laden jets. Simulations were performed for three different velocity ratios and two particle size distributions. Results showed very good agree-ment with the experimental observations, both qualitatively and quantitatively. The LES coupled with DQMOM approach was able to accurately model a variety of different particle behaviors, such as particle jet waviness, spreading, break up and particle clustering, for both fine and coarse particle size distributions. The development of the mean axial velocity along the centerline for the gas phase and the solid phase was also analyzed. The experimental values were predicted by the simulations with a maximum error below 12% and single phase simulations suggest that a significant part of this error can be attributed to the gas phase modeling. Segregation effects were also captured as small particles dispersed faster than the large ones. Moreover, simulations showed the importance of correctly discretizing the particle size phase-space for small particles and the particle velocity phase-space for large particles in the quadrature approximation. Overall, this study provides a solid validation of the LES with DQMOM approach to model particles in turbulent flows and justifies its further application to coal flame simulations. 47 2.7 Acknowledgements This material is based upon work supported by the Department of Energy under Award Number DE-NT0005015. The views and opinions of authors expressed herein do not neces-sarily state or reflect those of the United States Government or any agency thereof. CHAPTER 3 IGNITION OF COAXIAL TURBULENT DIFFUSION OXY-COAL JET FLAMES: EXPERIMENTS AND SIMULATIONS COLLABORATION Julien Pedel ; a, Jeremy N. Thornocka, Philip J. Smitha Submitted to Combustion and Flame Corresponding author. Email: jpedel@gmail.com aInstitute for Clean and Secure Energy, University of Utah, 155 South 1452 East, Room 350, Salt Lake City, UT, 84112, USA 3.1 Abstract The primary purpose of the study is to obtain physical insight into the stability of an oxy-coal flame as a function of changing primary gas composition (CO2 and O2). Flame stability was measured by Zhang et al. using optical measurements of the flame stand-off dis-tance in a 40 kW pilot facility. Large Eddy Simulations (LES) of the facility were performed using a multiscale simulation tool and provide additional insight into the experimentally observed data. The importance of factors such as heterogeneous reactions, radiation or wall temperature can be better understood thanks to simulations. The effects of three param-eters on the flame stand-off distance were studied. Simulation predictions were compared to experimental data using the data collaboration method. Overall, this study shows that high-fidelity LES simulations combined with experimental data can yield a deeper under- 49 standing of very complex coal flames, indicate where experimental uncertainties lie and be a valuable tool for design, retrofit and scale-up of oxy-coal burners. 3.2 Introduction Coal is an important energy resource for electricity production as its reserves are far more abundant than those of other fossil fuels. As environmental regulations get stricter, coal thermal power plants need to control pollutants emissions such as NOx, SOx and ash particles. Coal also plays an important role in global CO2 emissions, which will need to be reduced to limit the effects of climate change [66]. Oxy-coal combustion, in which an O2=CO2 mixture replaces air, is one of the few possible capture technologies to enable CO2 sequestration for existing coal-fired boilers. Despite the numerous research efforts that have been conducted on oxy-fuel combustion, several gaps of knowledge have been identified and need to be addressed in more detail to obtain a fundamental understanding of the changes between oxy-fuel combustion and conventional air-fired combustion [79, 95]. One area of interest is understanding how the replacement of air by a O2=CO2 mixture affects the flame ignition and stability for retrofit and new design of oxy-coal units [15]. Flame stabilization is a key technology to control NOx emissions because it increases the size of the fuel-rich flame region where NOx is reduced to nitrogen [6]. Furthermore, flame stability of oxy-fuel combustion strongly depends on oxygen concentration of the gas mixture [21] and can also affect the heat transfer performance of the plant. Thus, flame stabilization plays a very important role in controlling both NOx and CO2 emissions. Ignition of turbulent coal flames is a complex phenomenon involving multiphysics as-pects. The aerodynamics, kinetics and thermodynamics of both the gas and solid phase play an important role in understanding the flame characteristics. In the past, there have been numerous studies of single coal particle ignition mechanisms [7, 22, 25, 37, 49, 50, 67, 69, 89, 96, 110, 111, 114]. Both homogenous and heterogeneous mechanisms have been proposed to understand the ignition process. In the homogenous mechanism, the initial step is pyrolysis and subsequent ignition of the volatiles, followed by ignition of the char. The heterogeneous 50 reaction involves the direct attack of oxygen on carbon at the surface of the particle. The heterogeneous mechanism neglects the influence of volatile matter and depends on heat gen-eration by combustion on the coal surface, whereas the homogeneous mechanism involves the combustion of evolved volatile matter in the vicinity of the particle [37]. The domi-nant ignition mechanism is subject to controversy as Smoot and Smith [85] believed that reaction of oxygen on the coal surface was not a major factor in the ignition process while others identified surface oxidation as very important [22,25]. In addition, numerous factors, such as coal properties, particularly particle size and volatile content [7, 25, 49, 89], heating rate [25, 89], oxygen concentration [22, 67, 68, 89], pressure [69, 89, 110] and nature of the gas flow surrounding the particle [50, 67, 68] may all impact coal ignition. Oxy-coal combustion introduces an additional variable to those present in air-fired configurations, namely the concentration of O2 in the inlet oxidant mixture [15]. In practical turbulent diffusion flames, the overall inlet oxidant stream is split into at least two streams: a primary stream that transports the pulverized coal to and through the coal injector within the burner and a secondary stream that mixes with the primary coal jet downstream of the burner. In contrast to air combustion, the amount of oxygen within the primary jet in oxy-coal combustion is a variable, which can be set anywhere from zero up to any value that is considered safe. Previous work showed that flame stability is influenced by the PO2 in primary O2=N2 streams [15, 39, 51, 68, 96]. However, with the exception of the recent work of Molina and Shaddix [67,68], most studies have not considered the effect of elevated levels of CO2, as present in oxy-coal combustion systems, on the ignition of coal particles. The different properties of CO2 in comparison to N2 have been shown to cause differences in flame and furnace operation parameters such as ignition time and gas temperature profiles [58]. Molina and Shaddix [68] have suggested that differences between air- and oxy-fired coal particle ignition times and char burnout times are due to differences in molecular diffusion rates of O2 in N2 and CO2. Hu et al. [42] and Zeng et al. [110] have shown that the longer ignition delay of coal/char particles can be attributed to reduced O2 diffusion to the particle surface. 51 Although the reaction kinetics involved in coal particle ignition have been investigated extensively, most of the studies were done on single coal particles and on laminar flows. The single coal particle ignition mechanism cannot, however, sufficiently elaborate on the real turbulent flames due to certain limitations. First, the particles injected always contain a size distribution which should be taken into account to obtain correct kinetics [7]. Second, interactions between particles play a crucial role in the ignition and combustion behavior and pollutant formation and destruction. Annamalai et al. [3, 4] studied isolated coal/char particles and found that the particles with relatively low volatile matter, which ignite het-erogeneously when isolated, ignite homogeneously under cloud conditions. Finally, in real field boilers, with the mixing effect introduced by the large eddies in turbulence, single coal particle ignition mechanism must be coupled with turbulent mixing to understand or predict the ignition and the flame stability. Turbulence can easily wipe off the gradients of con-centrations, temperatures or velocities, which dramatically changes the ignition mechanism. The fundamental mechanics of coal particles in turbulent flows in general are as yet poorly understood [14, 24, 52] and have not been included in optimizing the design and operating flow condition for coal combustors. Smith et al. [84] showed that particle behaviors such as particle dispersion and particle clustering affected operational considerations within a coal-fired boiler such as flame stability, pollutants formation or heat fluxes. The gap between the fundamental research and the real industrial applications is thus still large and to understand how turbulent coal flames are stabilized, one must consider numerous complex physical processes and their interactions. This issue can be addressed by using Computational Fluid Dynamics (CFD), which was shown to be a very useful and accu-rate modeling framework for the description of particulate systems, and for the description of turbulent flows in general [24, 29]. Turbulence models based on the Reynolds-Averaged Navier-Stokes (RANS) equation are commonly used with gas and coal combustion models for the simulation. There is however evidence that RANS can be inaccurate in predicting particle-laden flows [5,86] and a turbulence model based on large eddy simulation (LES) has 52 been convincingly shown to be superior to RANS in accurately predicting turbulent mixing and combustion dynamics [56, 75]. This work is directed towards understanding the ignition mechanisms of coaxial tur-bulent diffusion, oxy-coal flames and how operating parameters such as the composition of the coal transport medium (primary CO2 and O2) or the furnace temperature affect the flame standoff distance. The focus is on interaction mechanisms between turbulent mixing in coaxial jets and coal particle ignition, rather than on each physical or chemical process taken individually. Flame stability was measured by Zhang et al. [113] using optical mea-surements of the flame stand-off distance in a 40 kW pilot facility. Experimental results suggest that flame stability is affected by primary PO2 , secondary preheat temperature, sec-ondary PO2 and the transport medium. Simulations of the facility are performed using a massively parallel LES simulation tool and provide additional insight into the experimental data. The importance of factors such as heterogeneous reactions, radiation, primary stream temperature or wall temperature is evaluated through the simulation results. A formal Validation and Uncertainty Quantification study is conducted with the Data Collaboration method [27, 33, 35]. The goal of this study is to show that LES can be a valuable tool for predicting and understanding oxy-coal flames characteristics and that combining experimen-tal data with high-fidelity modeling can guide experimentalists in their measurements by indicating where uncertainties lie. 3.3 Experiments This section briefly summarizes the experimental results published by Zhang et al. [112, 113] and the experimental conditions. A special attention is given to the uncertainties associated with the measurements. The experiments were carried out in the University of Utah 40 kW (100 kW maximum rating) down-fired oxy-coal combustor. The furnace consists of an oxy-fuel combustion chamber, followed by downstream-controlled temperature cooling to simulate practical fur-nace conditions [15]. It allows the systematic control of burner momentum and velocity 53 variables as well as wall temperatures. Quartz/sapphire windows allow nonintrusive optical diagnostics to measure the flame stability, quantified by the stand-off distance. The burner applied was a coaxial Type 0 (zero swirl) burner with a transport stream and pulverized coal in the center and a secondary stream in the annular sleeve. The coal employed was Utah bituminous pulverized coal. The mass average size of the particles is approximately 60 m (the coal particle size distribution is given in Table 3.1). The coal analysis is as follows: (1) Ultimate (wt%, daf): C 77.75%, H 5.03%, N 1.44%, S 0.45%, O 15.33%. (2) Proximate: moisture 3.03%, volatile matter 38.81%, fixed carbon 46.44%, ash 11.72%. (3) Higher (gross) heating value (HHV) = 27,286 kJ/kg. The flame stability was quantified by the observed stand-off distance, the distance between the burner tip and the visible ignition of the flame. An EPIX CMOS camera was set up at an exposure time of 8.3 ms and a frame rate of 30 fps to take sequences of flame images. Flame envelope/structure was defined at the above camera settings, which are consistent with the averaging process of the human eye. An image processing methodology was developed to quantify the flame stand-off distance from image data. The following operating parameters were fixed: coal feeding rate 4.84 kg/h, stoichio-metric ratio 1.15, wall temperature 1283 K, primary stream velocity 5.4 m/s and primary stream temperature 305 K. Both the primary and secondary streams were O2=CO2 mixtures and the overall oxygen concentration was maintained to 40 vol%. In this study, we are inter- Table 3.1: Utah bituminous coal particle size distribution Particle size ( m) Volume % 178 1.2 152 2.5 106 13.8 75 12.7 63 58.6 44 5.3 37 5.9 54 ested in a specific set of experiments where the secondary stream preheat temperature was set to 544 K and the primary PO2 varied from 0% to 20.9%. Variations in primary PO2 were obtained by transferring oxygen flows from the primary to the secondary stream, with only small effects on jet momentum. The averaged flame stand-off distance was measured with the CMOS camera for different values of primary PO2 . Each experiment was replicated five times. Results are presented in Table 3.2 and show that the flame stand-off distance decreases with increasing primary PO2 . Figure 3.1 shows the averaged values over the five replicates of the stand-off distance with associated error bars, where the upper and lower bounds correspond respectively to the maximum and minimum average value over the five replicates. This set of data was chosen as the main focus of this study because it is unclear how the primary PO2 affects the flame stand-off distance. High-fidelity computer simulations may provide insights on the ignition mechanisms and help design further experiments. This set of data will then later be used for the validation and uncertainty quantification analysis with the ARCHES LES tool provided in Section 3.6.2. Sources of uncertainty in the operating conditions were identified during those experi-ments. First, the wall temperature was measured with a thermocouple located behind the wall heater. The reported wall temperature value of 1283 K is therefore the value measured behind a 2 inches thick layer and not the wall temperature inside the furnace. The later one could indeed be much higher. An infrared camera was used to determine the wall tem-perature inside the University of Utah furnace in other experiments with similar operating conditions and measured wall temperatures at about 1800 K. In addition, thermocouple Table 3.2: Average standoff distances (cm) for overall PO2=40%, secondary preheat 544 K PO2;pri (vol%) 0 5.4 9.9 15.4 20.9 replicate 1 29.7 29.0 28.8 28.5 28.2 replicate 2 30.3 29.9 21.2 7.6 7.3 replicate 3 32.0 32.2 30.1 9.6 7.7 replicate 4 28.0 27.6 8.7 7.2 7.6 replicate 5 28.9 15.1 8.9 8.6 8.5 average 29.8 26.8 19.5 12.3 11.9 55 0 5 10 15 20 0.05 0.1 0.15 0.2 0.25 0.3 Primary P O 2 (%) Stand−off distance (m) Fig. 3.1: Average standoff distances with associated error measurements can be inaccurate and have up to 50K uncertainty. The primary stream was injected at room temperature and its temperature was assumed to be 305 K but was not measured directly. Operators noticed that the tip of the burner was glowing red after the experiments. The burner penetrates inside the furnace by about 10 cm and its tip is there-fore exposed to a high radiative flux. As a result, it is possible that the burner transfers heat to the primary stream and the coal particles, which would result in a higher primary stream temperature. 3.4 Numerical modeling Simulations were performed using ARCHES, a LES tool resulting from a ten-year part-nership with the Department of Energy and the University of Utah. It is a massively parallel code that solves conserved quantities (mass, momentum, energy, scalar) spatially and tem-porally in a turbulent flow field, allowing for detailed and accurate simulations of fires and flames [87]. ARCHES includes particle phase representation through the use of a moment 56 method, the direct quadrature method of moments (DQMOM). DQMOM combined with LES allows ARCHES to perform highly accurate simulations of coal combustion and gasi-fication applications. The code is integrated into a C++ framework called Uintah which provides large-scale parallelization tools for physics components [9,16,73]. ARCHES is main-tained in a repository and distribution is freely available [41]. The governing equations and numerical methods for the gas phase and the dispersed phase are described in the following sections. 3.4.1 Governing equations of the gas phase The gas phase is solved in a Eulerian manner by a finite volume method. Favre-filtered governing equations for continuity, momentum, enthalpy and scalars transport are solved. A second-order central difference scheme is used to calculate the convection and diffusion terms of the transport equations. For the continuity equation, the mass source term due to coal combustion is taken into account. Coal combustion models are described in Section 3.4.2. For the momentum equation, the Navier-Stokes equations are solved by taking the momentum source due to coal combustion into account. The subgrid scale (SGS) stress tensor is calculated by the dynamic subgrid scale model [36]. For the enthalpy equation, radiation, heat exchange between gas and coal particles, and heat transfer due to coal combustion are taken into account. The radiation source term is calculated by the discrete ordinate radiation method [2, 46, 47]. The gas-phase combustion is modeled through a mixture fraction approach and chemical equilibrium for the gas phase properties is assumed [93]. Spinti et al. [87] gave a more detailed description of the gas phase modeling in ARCHES. 3.4.2 Governing equations of the dispersed phase In order to describe the evolution of the disperse particulate phase, LES equations must be coupled with a population balance equation (PBE). The method of moments (MOM) 57 offers an attractive alternative to the traditional Lagrangian method where the PBE is tracked through its moments by integrating the internal coordinate [43]. The Direct Quadra-ture Method of Moments allows the description of population of particles with multivariate PBE, where two or more properties of the population are simultaneously tracked [26, 59]. DQMOM has been shown to be a powerful approach for describing polydisperse solids un-dergoing segregation, growth, aggregation and breakage processes in the context of CFD simulations [26, 32, 60] and to greatly reduce the computational cost compared to the La-grangian approach [17, 18]. In previous work, the advantages of using DQMOM for treating particle populations with low Stokes numbers (i.e., the dispersed-phase velocity follows closely the velocity of the continuous phase) have been clearly demonstrated [61, 63, 98, 115]. However, coal particles, which typically range from 1 micrometer to 200 micrometers in pulverized coal flames, have a finite Stokes number and their velocity depends on the particle size. In quadrature methods, this dependency is accounted for by solving an Eulerian model where each quadrature node has its own velocity field. Fox et al. [32] have applied DQMOM to droplets spray with finite Stokes number and compared the results with the multifluid method and a classical Lagrangian solver. Results show that DQMOM can provide a comparable accuracy with a significantly lower computational cost, making it a very good candidate for more complex two-phase combustion applications. The DQMOM method was implemented into ARCHES to solve the dispersed phase in an Eulerian manner. ARCHES coupled with DQMOM has been proven to be a reliable and accurate method to simulate particles in turbulent flows [74]. A brief description of the particles representation and of the DQMOM equations is given in the following paragraphs. A more detailed description of DQMOM implementation in ARCHES was previously given by Pedel et al. [74]. 58 3.4.2.1 Number density transport equation The particle Number Density Function (NDF) n describes the number of particles per volume as a function of the spatial location and of other independent variables. These independent variables are intrinsic characteristics of the particles (e.g., particle diameter, particle velocities, etc.) and are called internal coordinates. The internal coordinates are contained in the internal coordinate vector , and are dependent on space and time. At a given location in space and time (x0; t0), the total number of particles in a volume V is given by: np(x0; t0) = V n( ; x0; t0)d dV (3.1) The velocity of the NDF in both real and internal coordinate space can be described using a set of internal coordinates = [ 1; 2; :::], and external coordinates x = [x1; x2; x3]. For a single particle, the time derivatives of these variables may be written as dxi dt = vp;i (3.2) d j dt = Gj (3.3) where vp;i is the ith component of the particle velocity vector vp. The NDF transport equation is given by: @ @t (n( ; x; t)) + @ @xi (hvp;ij i n( ; x; t)) = @ @ j (hGj j i n( ; x; t)) + h( ; x; t) (3.4) In this equation, the variable hvp;ij i represents the velocity in real space of the NDF. This quantity is an average over particles having the same coordinates ( ; x; t) and is a distribution conditional to the value of the vector . Similarly, the quantity hGij i describes 59 the velocity of the NDF in the phase space. The term h is a source term representing the birth and death of particles in the domain. This term is set to zero as particle breakage and aggregation effects are neglected. 3.4.2.2 Direct Quadrature Method of Moments The Direct Quadrature Method of Moments (DQMOM) approximates the NDF by a summation of multidimensional Dirac delta functions. Figure 3.2 gives an example for bi-variate Gaussian distribution (i.e., two internal coordinates) of the NDF and its quadrature approximation with three quadrature nodes. In this approximation, the following assump-tions are made: - For a given quadrature node , at a given point (x,t), there is only one characteristic averaged internal coordinate value h i , also called abscissa, - The internal coordinate dispersion around the averaged value h i is zero. With those assumptions, the NDF is written as: n( ; x; t) = XN =1 0 @w YN i=1 ( i h ii ) 1 A (3.5) where both the weight w , which defines the number of particles per volume associated with the quadrature node , and the average internal coordinate value h i depend on space and time, but the dependence is omitted for clarity of notation. N is the number of quadrature nodes and N the number of internal coordinates. The validity of those assumptions has been discussed for particle velocity in the multi-fluid method [17]; the hypothesis is valid for particles with low Stokes number and becomes inappropriate when Stokes number increases (St > 0.4) and particle trajectory crossing oc-curs. It is however important to note that this limitation can be partially overcome with DQMOM. In the multifluid model, only the particle size space is discretized, whereas DQ-MOM allows us to discretize a multidimensional space. It is therefore possible to discretize the particle velocity space if the particle velocity is chosen as an internal coordinate and then 60 Fig. 3.2: DQMOM approximation of the NDF with 3 quadrature nodes to have several quadrature nodes with different velocities and otherwise identical internal coordinates values. As quadrature nodes can cross each other, the error associated with finite Stokes number particles crossing each others is reduced. DQMOM solves transport equations for the weights and weighted abscissas of the quadrature approximation. These transport equations are simply scalar transport equa-tions, and are easily implemented into a CFD code: @ @t (w ) + @ @xi hvp;ii w = a (3.6) @ @t (&n ) + @ @xi hvp;ii &n = bn (3.7) 61 where &n = w h ni is the weighted abscissa for the nth internal coordinate and hvp;ii is the average internal coordinate velocity. The right-hand side terms a and bn; are obtained by solving a linear system of equations. 3.4.2.3 Single coal particle description Following Smoot and Smith [85], a single coal particle can be characterized using several particle independent variables. These are denoted: - Raw coal cj - Char hj - Moisture wj - Ash aj - Particle diameter dpj - Particle enthalpy hpj - Particle velocity vector vpj The above quantities use the subscript j to denote the jth particle. The variable r can be used to denote reaction rates, so that rhj would be the net char reaction rate for the jth particle. Using this nomenclature, physical processes important to coal particles are depicted in Figure 3.3. The raw coal can react to form gaseous volatile matter in devolatilization reactions and solid char; the solid char can be oxidized to form more gaseous products. The ash mass is fixed, and ash is treated as inert. In this study, a total of seven internal coordinates were chosen to represent particle characteristics in the DQMOM formulation: - mass of raw coal c - mass of char h - particle diameter dp - particle enthalpy hp - particle velocity components vx, vy, vz 62 Fig. 3.3: Illustrative schematic of coal particle components and reactions. The particle diameter dp was assumed to be constant. To reduce the number of internal coordinates and the computational cost, the moisture in the particles was assumed to be evaporated before entering the domain and the water content was added to the inlet stream. 3.4.2.4 Particle velocity model In this work, the aggregation/breakage phenomena are neglected and the only physical process modeled is the drag force, which can be described by the Stokes drag law. For a mesoscale size particle, other forces such as lift can be neglected and the momentum equation for the particle can be expressed by an ordinary differential equation. dvp;i dt = fdrag p (ug;i vp;i) + gi ( p g) p (3.8) where g is the gravity force acting on the particle, g and p are respectively the gas and particle densities, ug is the gas velocity, and fdrag is the coefficient of the drag force, which has a close relationship with the particle Reynolds number [85]: 63 fdrag = 8>>>>< >>>>: 1 Rep < 1 1 + 0:15Re0:687 p 1 < Rep < 1000 0:0183Rep Rep > 1000 (3.9) Rep = pdp jup ugj g (3.10) where g is the gas dynamic viscosity. In equation (3.8), p is the particle relaxation time: p = pd2 p 18 g (3.11) 3.4.2.5 Devolatilization model In this study, an empirical devolatilization model proposed by Yamamoto et al. is used [104]. The approach was proven accurate and inexpensive by comparison to more sophisticated models such as the distributed activation energy model , FLASHCHAIN [71] or chemical percolation devolatilization (CPD) [31]. The devolatilization process is described as: (raw coal) rv ! Y (volatiles) + (1 Y ) (char) (3.12) where Y is a stoichiometric coefficient. The devolatilization rate of the model rc is based on a single reaction: rv = d c dt = Fk0exp (Ev=RTp) c (3.13) 64 F = exp X5 i=0 ciGi v ! (3.14) where k0 is the frequency factor, Ev is the activation energy and R is the gas constant. F is the modification factor and is expressed as a function of fraction devolatilized Gv = ( c;0 c)= c;0. Here c;0 is the initial mass of raw coal. The model parameters Y , k0, Ev and ci are determined by minimizing the difference between CPD and the present model at the heating rate of 100,000 K/s. Figure 3.4 shows results calculated by the Yamamoto model and CPD. The coal used here is the same as the coal used for the experiments. The results obtained by the present model show good agreement with CPD results. The maximum error of the fraction devolatilized is below 1.5% for both heating rates of 100,000 and 20,000 K/s. This model thus offers the advantage to predict the fraction devolatilized with the same accuracy as the CPD model but with the low computational cost of a single-rate reaction model. 3.4.2.6 Char oxidation model After the raw coal in a particle has devolatilized, it forms char. Historically, most kinetic data on burning coal char particles have been interpreted using an nth order Arrhenius model of char combustion. In this representation, the global surface reaction rate follows the expression q = ks (Tp) Pn O2;s (3.15) where n is the reaction order, Tp is the particle temperature, and ks is a temperature-dependent rate coefficient, assumed to follow an Arrhenius form: ks (Tp) = Asexp (Es=RTp) (3.16) 65 500 600 700 800 900 1000 1100 1200 1300 1400 1500 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Particle temperature (K) Fraction devolatilized CPD, 1.0E5 K/s CPD, 2.0E4 K/s Yamamoto, 1.0E5 K/s Yamamoto, 2.0E4 K/s Fig. 3.4: Verification results of the Yamamoto devolatilization model with comparisons to results by the chemical percolation devolatilization model (CPD) The parameter values were chosen from the experimental study of Murphy and Shaddix [70] for Eastern bituminous coal: As = 344 mol=s=m2=atmn , Es = 45:5 kJ=mol and n = 0:18. To be consistent with the experiments, the same equations used to determine these parameter values [70] are solved to calculate the surface reaction rate q. The partial pressure of oxygen at the char surface PO2;s is obtained by solving the gas-phase diffusion, PO2;s P = + PO2;1 P exp qdp 2CDO2;mix (3.17) In this equation, C = P=RTf is the average molar concentration of the gases in the boundary layer and = 1 2 where = CO2=CO= (1 + CO2=CO) is the fraction of carbon that becomes CO2. The CO2=CO production ratio at the particle surface is determined from the correlation CO2=CO = A0P 0 O2;sexp (B0=Tp) with the coefficients suggested by Tognotti et al. [91]. Given a trial value of q, Equation 3.17 must be solved iteratively because depends on the CO2=CO ratio, which in turn depends on PO2;s. 66 The char oxidation rate is then given by rh = d h dt = d2 pWCq (3.18) where Wc is the molecular weight of carbon. To account for the eventuality of heterogeneous ignition, the oxidation reaction was also applied to the raw coal before ignition. As shown in Sections 3.6.1 and 3.6.2, the oxidation rate was negligible before ignition, even though the activation energy value Es is significantly lower than oxidation activation energy values reported in previous studies at low temperature. There is indeed a relative agreement about the magnitude of the global intrinsic activation energy for carbon oxidation at a low temperature regime where direct measurement is possible; it falls between 105 and 180 kJ/mol with many values in the range 130 to 150 kJ/mol [30, 44, 83]. 3.4.2.7 Particle heat transfer The particle is heated by convection, radiation and reaction enthalpy changes: dhp dt = Qrad + Qconv + rjhj ; (3.19) where Qrad represents the net radiation flux to a particle, Qconv represents energy transfer due to convection and conduction between the gas and the particle, and rjhj represents both the amount of energy lost by the particles due to lost mass and the enthalpy released during devolatilization and oxidation reactions. Convection The convection term is expressed as: Qconv = Nu kg (Tg Tp) dp =2 e =2 1 (3.20) 67 where Tg is the gas temperature, kg is the gas thermal conductivity, and = WCqdpCpg=kg is a modified version of the Peclet number [70]. Cpg is the gas heat capacity. The Nusselt number Nu is calculated with the following correlation [55]: Nu = 2:0 + 0:65Re1=2 p Pr1=3: (3.21) The particle temperature Tp is obtained from the particle enthalpy by assuming Mer-rick's correlations for the particle heat capacity [65]. Radiation The radiative flux is given by: Qrad = Qincident Qemitted (3.22) where Qincident is the incident radiative flux to the particle and Qemitted is the radiative heat flux emitted by the particle. The absorption coefficient of the coal particles is assumed constant: Qabs = 0:8. The incident flux is obtained through the Discrete Ordinates radiation model whereas the emitted flux is calculated using the blackbody emissive power: Qemitted = Qabs d2 p T4 p : (3.23) Following [85], the coupling with the gas phase is done by modifying the gas absorption coefficient: kabs;g = kabs;g + kabs;p (3.24) with kabs;p = Xnp j=1 4Qabsd2 pj (3.25) 68 3.5 Validation and uncertainty quantification The common approach to model validation consists in comparing the predictions made by the model to available experimental data. Typically, such comparisons result in mixed outcomes: some show a good agreement and some do not. In the latter case, the apparent inconsistency between the model and the experiment is argued to imply either that the model is inadequate or that the experiment is incorrect. Before a decision can be made as to whether or not to use a model, it is essential for the decision maker to understand the uncertainty in the predictions made by the model. Simulation error or uncertainty arises from four types of uncertainties: • Model uncertainty, or error in the mathematical description of the physical process being modeled. This error includes uncertainty in the parameters used in these models. • Uncertainty in boundary conditions. This uncertainty derives from the uncertainty in the inputs and boundary conditions used for the simulation. • Numerical error (verification error). This error arises from the numerical approxima-tions used in producing the simulation (i.e., grid resolution, discretization error, convergence error, etc.) • Experimental error in validation data. The quantification of simulation error is achieved by comparing the simulation output to experimental observations; since the exper-imental measurements themselves are inexact, the simulation accuracy can not be known to any better degree of accuracy than the measurements themselves. Because of these inaccuracies, it is essential to simultaneously estimate the sensitivity of the assessed variability in the model output to changes in the assumptions concerning the model uncertainties, as well as to changes in the experimental measurements. This approach is illustrated in Figure 3.5. It tightly couples simulation and experi-mental data with formal verification and validation /uncertainty-quantification (V&V/UQ) methodologies. More specifically, given a bounded parameter space, this methodology seeks to find the minimum and maximum error by propagating the error/uncertainty in the input parameter space (priors) to the simulation output or posterior subject to the constraints 69 Fig. 3.5: Integration of uncertainty in modeling parameters, boundary conditions and ex-perimental data into quantification of an overall predictive output. of the experimental validation and verification uncertainty. This methodology provides the lower and upper bound on the predicted efficiency by requiring consistency between the simulation output and the measured experimental data. This works presents a formal validation and uncertainty quantification analysis of the ARCHES LES tool with the experimental data on oxy-coal flame ignition presented in Section 3.3. The Data Collaboration approach [27, 34, 35] is used to study the uncertainty associated with a couple of scenario parameters which were identified as potential sources of uncertainty. The following section first provides a brief description of the Data Collaboration method and then explains the methodology chosen to perform the analysis. Finally, a detailed description of the simulations is provided. 3.5.1 Data Collaboration method Recently, the Data Collaboration approach has been developed to provide quantitative and reliable uncertainty bounds on predictions from multiphysics, multiscale simulations [27, 34, 35]. Data Collaboration is a method that unites process models and associated admissible parameter values with experimental data and accompanying uncertainties. The 70 approach focuses on transferring the uncertainties of the experimental data into the model directly. Data obtained from various sources are integrated through models that depend on common parameters. This ensemble is called a dataset. Optimization-driven techniques can then be applied to dataset analysis, including performing dataset-based predictions [34] and determining dataset consistency [28]. Doing so allows one to harvest substantially more of the information content of the data [35] and determine more realistic bounds on model predictions [34]. 3.5.1.1 Dataset In the Data Collaboration approach [28], each experiment e is associated with a dataset unit Ue = (de; ue; le; Me) where de is the measured value, ue and le are respectively the upper bound and lower bound of the reported uncertainty and Me is a mathematical model of the experiment. The model depends only on a parameter vector to generate the prediction Me for ye. The true value of the experimental observable, ye, satisfies le ye de ue. A dataset is a collection of dataset units. The dataset is denoted as D and the set of indices e as E, i.e., D = fUege2E. Complex models usually have a large number of parameters. In the Data Collaboration method, a subset of model parameters which have a influence on the outcome is chosen. The parameters contained within this subset are called active variables. An individual active variable is designated as Xj and xj 2 R refers to a specific value of Xj . Individual dataset units may have different sets of active variables. Thus, Xj might be an active variable for one experiment but not another. We denote the list of active variables for experiment e as Xe. The total number of active variables is denoted as n and the vector x 2 Rn represents the active variable values. Associated with a vector x, xe are the values extracted from x that correspond to the active variables for experiment e. 71 3.5.1.2 Initial hypercube and feasible region The Data Collaboration approach assumes that prior information on the possible values of the dataset active variables is available. This prior information can be expressed as the confinement of possible values of the active variables to an n-dimensional \"hypercube,\" H = fx 2 Rn : j xj jg, where j and j are the lower and upper bounds on xj for j = 1; 2; :::; n. Each edge of the hypercube H represents the presumed interval of physically allowed values of the corresponding active variable. Some parameter values drawn from H may result in model predictions that lie outside the experimentally determined ranges. In other words, not every x 2 H predicts all exper-imental observations of the dataset within their specified uncertainties. The collection of parameter values that are both contained in the hypercube and satisfy le Me(xe)de ue for every e 2 E form the feasible region, FD. A point x that is not contained in FD can thus be eliminated from consideration as a possible value for the dataset active variables. An experiment may therefore eliminate portions of the hypercube H from consideration, thereby decreasing the uncertainty in the values of the parameters. 3.5.1.3 Dataset consistency The Data Collaboration approach also allows us to determine wether the data contained in the dataset are mutually consistent. A dataset D is said to be inconsistent if there is no single point x in the hypercube H that satisfies le Me(xe) de ue for all e 2 E . Otherwise, the dataset is consistent. In other words, the dataset is inconsistent if the feasible region is empty. The consistency of a dataset depends on the size of the hypercube. A consistent dataset may become inconsistent with a decrease in the intervals of active variables and/or in the level of experimental error. 3.5.2 Methodology In this study, we are interested in understanding the ignition phenomena in oxy-flames and the influence of operating conditions on the flame stability. The experimental data with 72 their associated uncertainties presented in Section 3.3 were selected for the validation study. Experiments suggest that the primary PO2;pri has an influence on the flame stand-off distance and was therefore chosen as a parameter of interest. Two other parameters were chosen: the wall temperature Twall and the temperature of the primary stream Tpri. Those parameters were selected because experiments show a high sensitivity to the wall temperature and the stream temperature [112, 113], but also because of their potential uncertainty. As described in Section 3.3, the wall temperature was measured at 1283 K but could be as high as 1800 K. The uncertainty range for Twall was thus assumed to be 1200 K - 1800 K. Similarly, the primary stream temperature could be higher than reported in the experiments since the burner tip is heated by radiation. We assumed that Tpri could be as high as the secondary stream temperature (544K). The range associated with each parameter is reported in Table 3.3. 3.5.3 Surrogate model The Data Collaboration method uses optimization procedures to determine the effects of input variables and associated uncertainties of system responses. This procedure requires a large number of function evaluations. Ideally, the function evaluations would be performed with ARCHES LES tool. However, LES simulations are computationally expensive and running hundreds of LES simulations is not feasible. In this situation, it is useful to have a simpler model, also called a \"surrogate model,\" that is much cheaper to evaluate and much simpler in form, such as a polynomial. These can be constructed based on limited information about the larger scale computer simulation model. Box and Wilson first proposed the response surface methodology (RSM) [13]. RSM uses polynomial surfaces to represent the response as a function of input variables. Polyno- Table 3.3: Chosen parameters and their ranges Parameter Range Wall temperature (Twall) 1200-1800 K Primary stream temperature (Tpri) 305-544 K Primary oxygen partial pressure (PO2;pri) 0-20.9 % 73 mials are very general and work well for many responses, and as a result, RSM has thrived as a widely used modeling technique for complex systems. In this work, a Box-Behn"}]},"highlighting":{"195472":{"ocr_t":[]}}}