Application of a validation/uncertainty quantification (vuq) methodology at two scales: from modeling of char oxidation to simulation of a 1.5 mw coal-fired furnace

This dissertation presents two validation/uncertainty quantification (VUQ) studies, one on the 1.5 MWth coal fired furnace (L1500) and the other on the Reacting Particle and Boundary Layer (RPBL) char oxidation model. A six-step methodology is used in both cases. In Chapters 2 and 3, the VUQ for the L1500 furnace is presented; the quantities of interest (QOIs) are the heat removal by the cooling tubes and the wall temperature. In Chapter 2, the Arches simulation of the L1500 base case is described in detail. From the simulation, two models that impact the QOIs are selected for further analysis, the ash deposition model and the char oxidation model. An input and uncertainty (I/U) map is created from the parameters in these models and a sensitivity analysis is performed with the five parameters that have the greatest impact on the QOIs. From the sensitivity analysis, two parameters (thermal conductivity of the deposit and wall emissivity) are chosen for the next steps in the VUQ cycle. In Chapter 3, an updated version of the I/U map with two additional parameters, the coal feed rate and the swirl factor (this factor varies the tangential component of the velocity), is presented. The thermal conductivity of the deposit and wall emissivity are combined into one parameter, the effective thermal conductivity. These three active parameters are then used in the consistency analysis. The experimental uncertainty of the QOIs is estimated by adding the sampling and the systematic errors. Data collection for the simulation is done with 34 cases obtained by varying the three active parameters. For each experimental QOI, a Gaussian process (GP) surrogate model is built from the set of simulation data. The consistency analysis is performed with the GP surrogate models and the QOIs with their estimated uncertainties; consistency is achieved. Chapter 3 concludes with recommendations for reducing the uncertainty in the experimental measurements and a review of model assumptions. In Chapters 4 and 5, the VUQ for the RPBL model is presented; the QOIs are the particle temperature and velocity. In Chapter 4, the RPBL model formulation and the associated I/U map are given. One case is presented and explained in detail. A sensitivity analysis with nine parameters from the I/U map is performed, and five parameters (dp, ϕ_initial, Yc, ε_p, and r_inf/rp) are selected for the next steps. In Chapter 5, the priorities of the parameters in the I/U map are updated using the results from sensitivity analysis. To compute the experimental uncertainty associated with the QOIs, the main contribution is assumed to be the sampling error. The RPBL model is run with the five parameters to produce simulation data. A polynomial chaos (PC) surrogate model is built for the set of simulation data corresponding to each experimental QOI. The consistency analysis is performed with the PC surrogate models and the QOIs with their estimated uncertainties. Consistency is found for three different types of char and six different particle size ranges in two O2 environments, each with 6-14 QOIs. To conclude, Chapter 5 reviews the experimental measurements, analyzes what is learned about the parameters in the consistency analysis, and revisits the assumptions made in the RPBL formulation.

APPLICATION OF A VALIDATION/UNCERTAINTY QUANTIFICATION (VUQ) METHODOLOGY AT TWO SCALES: FROM MODELING OF CHAR OXIDATION TO SIMULATION OF A 1.5 MW COAL-FIRED FURNACE by Oscar Homero Dı́az Ibarra 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 2017 c Oscar Homero Dı́az Ibarra 2017 Copyright All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Oscar Homero Dı́az Ibarra has been approved by the following supervisory committee members: Jennifer P. Spinti , Chair(s) 8 Feb 2017 Date Approved Philip J. Smith , Member 8 Feb 2017 Date Approved Christopher R. Shaddix , Member 8 Feb 2017 Date Approved Eric G. Eddings , Member 8 Feb 2017 Date Approved Jost O. L. Wendt , Member 8 Feb 2017 Date Approved by Milind Deo , Chair/Dean of the Department/College/School of Chemical Engineering and by David B. Kieda , Dean of The Graduate School. ABSTRACT This dissertation presents two validation/uncertainty quantification (VUQ) studies, one on the 1.5 MWth coal fired furnace (L1500) and the other on the Reacting Particle and Boundary Layer (RPBL) char oxidation model. A six-step methodology is used in both cases. In Chapters 2 and 3, the VUQ for the L1500 furnace is presented; the quantities of interest (QOIs) are the heat removal by the cooling tubes and the wall temperature. In Chapter 2, the Arches simulation of the L1500 base case is described in detail. From the simulation, two models that impact the QOIs are selected for further analysis, the ash deposition model and the char oxidation model. An input and uncertainty (I/U) map is created from the parameters in these models and a sensitivity analysis is performed with the five parameters that have the greatest impact on the QOIs. From the sensitivity analysis, two parameters (thermal conductivity of the deposit and wall emissivity) are chosen for the next steps in the VUQ cycle. In Chapter 3, an updated version of the I/U map with two additional parameters, the coal feed rate and the swirl factor (this factor varies the tangential component of the velocity), is presented. The thermal conductivity of the deposit and wall emissivity are combined into one parameter, the effective thermal conductivity. These three active parameters are then used in the consistency analysis. The experimental uncertainty of the QOIs is estimated by adding the sampling and the systematic errors. Data collection for the simulation is done with 34 cases obtained by varying the three active parameters. For each experimental QOI, a Gaussian process (GP) surrogate model is built from the set of simulation data. The consistency analysis is performed with the GP surrogate models and the QOIs with their estimated uncertainties; consistency is achieved. Chapter 3 concludes with recommendations for reducing the uncertainty in the experimental measurements and a review of model assumptions. In Chapters 4 and 5, the VUQ for the RPBL model is presented; the QOIs are the particle temperature and velocity. In Chapter 4, the RPBL model formulation and the associated I/U map are given. One case is presented and explained in detail. A sensitivity analysis with nine parameters from the I/U map is performed, and five parameters (d p , φinitial , Yc , e p , and rin f rp ) are selected for the next steps. In Chapter 5, the priorities of the parameters in the I/U map are updated using the results from sensitivity analysis. To compute the experimental uncertainty associated with the QOIs, the main contribution is assumed to be the sampling error. The RPBL model is run with the five parameters to produce simulation data. A polynomial chaos (PC) surrogate model is built for the set of simulation data corresponding to each experimental QOI. The consistency analysis is performed with the PC surrogate models and the QOIs with their estimated uncertainties. Consistency is found for three different types of char and six different particle size ranges in two O2 environments, each with 6-14 QOIs. To conclude, Chapter 5 reviews the experimental measurements, analyzes what is learned about the parameters in the consistency analysis, and revisits the assumptions made in the RPBL formulation. iv To my parents, Oscar and Marı́a Elsomina CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv CHAPTERS 1. 2. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 1.2 1.3 1.4 1.5 Validation and uncertainty quantification methodology . . . . . . . . . . . . . . . . . . L1500 furnace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Char oxidation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Overview and organization of this dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 3 3 4 A VUQ ANALYSIS FOR THE L1500 FURNACE: SENSITIVITY ANALYSIS . . . 8 2.1 2.2 2.3 2.4 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Description of VUQ approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . L1500 experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 L1500 reactor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Swirl burner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Coal characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.4 Operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 L1500 simulation strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 STAR-CCM+ simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Arches simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2.1 Char oxidation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2.2 Ash deposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.3 Handoff plane boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.4 Domain size and mesh resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Analysis of quantities of interest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 I/U map for char oxidation and ash deposition models . . . . . . . . . . . . . . . . . . 2.8 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Computation of radiometer heat flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.2 Computation of heat removal by cooling tubes . . . . . . . . . . . . . . . . . . . . . 2.8.3 Computation of wall temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.4 Generating PC surrogate models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 10 11 13 13 14 14 15 15 16 16 18 19 20 21 21 22 24 24 25 26 26 26 28 2.10 Acknowledgment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3. A VUQ ANALYSIS FOR THE L1500 FURNACE: CONSISTENCY ANALYSIS . 40 3.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Description of VUQ approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Construction of input and uncertainty map . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Parameter selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Experimental data collection and uncertainty analysis . . . . . . . . . . . . . . . . . . 3.5.1 Experimental data collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Experimental measurement uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.1 Heat removal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.1.1 Systematic error for heat removal. . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.1.2 Sampling error for heat removal. . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.1.3 Total error for heat removal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.2 Wall temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.2.1 Systematic error for wall temperature. . . . . . . . . . . . . . . . . . . . 3.5.2.2.2 Sampling error for wall temperature. . . . . . . . . . . . . . . . . . . . . 3.5.2.2.3 Total error for wall temperature. . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.3 Oxygen molar fraction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2.3.1 Sampling error for O2 molar fraction. . . . . . . . . . . . . . . . . . . . . 3.5.2.3.2 Systematic error for O2 molar fraction. . . . . . . . . . . . . . . . . . . . 3.6 Simulation data collection and model development . . . . . . . . . . . . . . . . . . . . . 3.6.1 Wall model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Suite of simulations for consistency analysis . . . . . . . . . . . . . . . . . . . . . . . 3.6.3 Data extraction and postprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Construction of surrogate models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Analysis of model outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Feedback and feed-forward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9.2 Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9.3 Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Acknowledgment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 41 42 43 43 43 45 45 46 46 46 47 47 48 48 51 51 51 51 51 51 52 54 54 54 56 57 57 59 60 60 62 REACTING PARTICLE AND BOUNDARY LAYER MODEL: FORMULATION 73 4.1 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Development of the RPBL model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Species mass fraction equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Energy equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Solution methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Reactor description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Model verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.1 Input and uncertainty map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii 74 76 80 81 83 86 87 88 88 89 90 4.6.2 Base case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.3 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. 91 93 96 96 REACTING PARTICLE AND BOUNDARY LAYER (RPBL) MODEL: VUQ . . . . 108 5.1 5.2 5.3 5.4 5.5 5.6 Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Description of VUQ approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Selection of quantities of interest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Construction of input and uncertainty map . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Experimental data collection and uncertainty analysis . . . . . . . . . . . . . . . . . . 113 5.6.1 Reactor description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.6.2 Data collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.6.3 Measurement uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.7 Simulation data collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.7.1 Reactor model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.7.2 Data collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.8 Construction of surrogate models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.9 Analysis of model outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.9.1 Effect of particle size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.9.2 Effect of char type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.9.3 Effect of O2 concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.10 Feedback and feed-forward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.10.1 Particle size and shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.10.2 Other model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.10.3 Model assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.10.4 Uses for RPBL model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.11 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.12 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6. CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.1.1 L1500 furnace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.1.2 Char oxidation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 viii LIST OF FIGURES 1.1 Validation hierarchy for CCMSC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 Photo of the L1500 reactor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Combustion process for a coal particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Validation hierarchy for CCMSC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2 Six-step methodology with consistency analysis. The yellow box corresponds to the steps analyzed in this paper. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3 Drawing of the L1500 reactor located at the Industrial Combustion and Gasification Research Facility. a) Schematic of the L1500 multifuel combustion furnace. b) Schematic of the first four sections of the L1500 furnace. . . . . . . . . . 34 2.4 L1500 cooling tubes. a) Cooling tubes are located in the first four sections of the furnace. Photo was taken from the burner location looking toward the furnace exit. b) Dimensions of a set of cooling tubes (in inches). . . . . . . . . . . . . 34 2.5 Schematic of the swirl burner used in the L1500 furnace. Inlets are labeled. . . . 35 2.6 Swirl vanes in the burner at the 0% swirl and 100% swirl positions. . . . . . . . . . 35 2.7 Experimental PSD and fitted Rosin-Rammler PSD. Sieving was performed for different lengths of time to determine the effect on the resulting PSD. Based on this analysis, only the Beckman-Coulter diffraction data were used for the Rosin-Rammler fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.8 L1500 simulation coupling between Arches and STAR-CCM+ . . . . . . . . . . . . . . 36 2.9 Velocities at plane of the burner tip; left are STAR-CCM+ results with a resolution of 0.5x10−3 m, right are Arches results with a resolution of 15x10−3 m (ratio = 30). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.10 Shortened geometry for the L1500 simulations including the quarl, the eight sets of cooling tubes, and the step change in the reactor floor. Resolution is 15 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.11 Thermocouple placement in furnace wall. Section 4 is shown but the placement is similar for all thermocouple measurements. . . . . . . . . . . . . . . . . . . . . . . 37 2.12 Main sensitivity index for radiative heat flux measured by the radiometers. . . 37 2.13 Main sensitivity index for the heat removed by the cooling tubes. . . . . . . . . . . . 38 2.14 Main sensitivity index for wall temperature (thermocouple's location) computed with Arches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.15 Wall temperature at thermocouple locations computed with Arches. Horizontal red line is Tslag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 Validation hierarchy for CCMSC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.2 Schematic of the L1500 multifuel combustion furnace. . . . . . . . . . . . . . . . . . . . . 64 3.3 Schematic of the first four sections of the L1500 furnace. . . . . . . . . . . . . . . . . . . . 64 3.4 Six-step methodology with consistency analysis . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.5 Heat removal data in section 1 on February 27, 2015. The heat removal decreases continually in both 0% and 100% swirl modes due to ash deposition on the tube surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.6 Heat removal from cooling tubes with total error bars; see Eqn. (3.1). . . . . . . . . 65 3.7 Thermocouple placement in furnace wall. Section 4 is shown but the placement is similar for all thermocouple measurements. . . . . . . . . . . . . . . . . . . . . . . 66 3.8 Representation of the instrument model to correct thermocouple measurement in ceramic sheath to equivalent inside wall temperature. . . . . . . . . . . . . . 66 3.9 Procedure used to correct Ttc to Tw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.10 Thermocouple temperature and corrected wall temperature. . . . . . . . . . . . . . . . 67 3.11 Wall model verification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.12 Wall temperature at thermocouple locations; section 1 is closest to the burner. . 68 3.13 Suite of simulations that were run in Arches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.14 Verification of surrogate models for T1 , HR2 , and O2 molar fraction . . . . . . . . . 69 3.15 Surrogate models for T1 , HR2 and O2 molar fraction. . . . . . . . . . . . . . . . . . . . . . 70 3.16 Consistent sample space for all 14 QOIs in L1500 simulations. . . . . . . . . . . . . . . 70 3.17 Consistency range for heat removal (top), O2 molar fraction (middle), and wall temperature (bottom). Data are plotted by QOI number. The eight sets of cooling tubes are QOIs 1-8, the O2 mole fraction in the exhaust is QOI 9, and the wall temperatures in sections 1, 2, 3, 4, and 6 are QOIs 10-14. Consistency was obtained simultaneously across the 14 QOIs. . . . . . . . . . . . . . . . . . . . . . . . . 71 3.18 Effect of f s , ṁcoal = 275 [lb/h], k e f f = 3 [W/m/k ], Blue color is 300 K, red is 2100 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.19 Transient temperature effects in the refractory wall . . . . . . . . . . . . . . . . . . . . . . . 72 4.1 Representation of the particle and boundary layer. . . . . . . . . . . . . . . . . . . . . . . . 98 4.2 Schematic of optically accessible entrained flow reactor at Sandia National Laboratories. Figure adapted with permission from [34] . . . . . . . . . . . . . . . . . . . 100 4.3 Grid independence study for average particle temperature, d p = 80 µ m. . . . . . . 100 4.4 Grid independence study for particle velocity, d p = 80 µ m. . . . . . . . . . . . . . . . . . 101 4.5 Species mole fractions in the particle and boundary layer of 95 µm diameter char particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.6 Species mole fractions in 95 µm diameter char particle. . . . . . . . . . . . . . . . . . . . . 103 4.7 Carbon bulk density inside 95 µm diameter char particle. . . . . . . . . . . . . . . . . . . 104 x 4.8 Particle temperature (left) and particle and boundary layer gas temperature (right) for 95 µm diameter coal particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.9 Particle temperature sensitivity analysis for 50 µm diameter particle. . . . . . . . . 105 4.10 Particle velocity sensitivity analysis for 50 µm diameter particle. . . . . . . . . . . . . 105 4.11 Particle temperature sensitivity analysis for 80 µm diameter particle. . . . . . . . . 106 4.12 Particle velocity sensitivity analysis for 80 µm diameter particle. . . . . . . . . . . . 106 4.13 Particle temperature sensitivity analysis for 120 µm diameter particle. . . . . . . . 107 4.14 Particle velocity sensitivity analysis for 120 µm diameter particle. . . . . . . . . . . 107 5.1 Six-step methodology with consistency analysis . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.2 Schematic of optically accessible, entrained flow reactor at Sandia National Laboratory. Figure adapted with permission from [34] . . . . . . . . . . . . . . . . . . . . 131 5.3 Average temperature, velocity, and diameter with error bars for particles in the 75-90 µm size bin. The black dots are the experimental points, and the red error bars are computed from equation (5.1) with a confidence interval of 95 %. These data correspond to US char in an environment of O2 = 24 vol %, H2 O = 14 vol %, and balance CO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.4 Pseudocolor plots of gas temperature and gas velocity in the Sandia entrained flow reactor for O2 = 24 vol %, H2 O = 14 vol % , balance CO2 environment. . . . 133 5.5 Gas temperature profiles (single cell and averaged over 9 cells) along the center line of the entrained flow reactor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.6 Thermocouple temperature profiles along the center line of the entrained flow reactor. The profile is from position 0 [m] to position 0.254 [m]. . . . . . . . . 134 5.7 Comparison of spatially-averaged gas velocity profile from reactor simulation and experimental particle velocity profile for particles in the 53 - 63 µm size range. The average gas velocity was reduced by 3% to be in the uncertainty range of the small particle velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.8 Spatially-averaged mass fraction profiles for O2 , H2 O, and CO2 along the center line of the simulated entrained flow reactor. . . . . . . . . . . . . . . . . . . . . . . . 135 5.9 Verification plots for PC surrogate models. The red line is y = x. The positions correspond to an experiment with US char in an O2 = 24 vol%, H2 O = 14 vol% balance CO2 environment and bin size of 75-90 µm. . . . . . . . . . 136 5.10 Scatter plots for PC surrogate models; left column is particle temperature and right column is particle velocity. The positions (rows) correspond to an experiment with US char in an O2 = 24 vol%, H2 O = 14 vol% balance CO2 environment and bin size of 75-90 µm. Ranges for d p and φinitial are presented r in Table 5.1; Yc,intial =0.75 , rinpf = 85, and e p = 0.55. . . . . . . . . . . . . . . . . . . . . . . . . 137 5.11 Consistent sample for US char in an O2 = 24 vol%, H2 O =14 vol%, balance CO2 environment and 75-90 µm size bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 xi 5.12 Consistent range for US char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment and 75-90 µm size bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 5.13 Consistent range for I6 char in an O2 = 24 vol%, H2 O =14 vol%, balance CO2 environment for 53-63 µm, 63-75 µm, and 75-90 µm size bins. . . . . . . . . . . . . . . . 140 5.14 Consistent range for I6 char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment for 90-106 µm, 106-125 µm, and 125-150 µm size bins. . . . . . . . . . . 141 5.15 Consistent ranges for I6, US, and BT chars in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment for 90-106 µm size bin. . . . . . . . . . . . . . . . . . . . 142 5.16 Consistent range for US char (top) O2 = 24 vol% and (bottom) O2 = 36 vol%, both with H2 O = 14 vol% and balance CO2 environments for the 53-63 µm size bin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 xii LIST OF TABLES 2.1 Ultimate analysis for Utah Sufco coal. Data are averaged and normalized from the analysis of five samples taken from the bags of coal that were burned during the testing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2 L1500 operating conditions for 100 % swirl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.3 Shell temperatures averaged over all measurements made on a side . . . . . . . . . 30 2.4 Thermal properties obtained from the manufacturer for the L1500 walls . . . . . 30 2.5 Input/uncertainty map for char oxidation and deposition model . . . . . . . . . . . 31 3.1 Input/uncertainty map for consistency analysis . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.2 Input/uncertainty map for the thermocouple instrument model . . . . . . . . . . . . 62 3.3 Shell temperatures averaged over all measurements made on a side . . . . . . . . . 63 3.4 Parameter ranges after consistency analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.1 Surface reaction mechanism [31, 34] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.2 Input/uncertainty map for RPBL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.1 Input/uncertainty map for RPBL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.2 Parameter ranges prior to and after consistency analysis for US char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment and 75-90 µm size bin . . 128 5.3 Consistent parameter ranges for US, I6, and BT chars in an O2 = 24 vol%, r H2 O = 14 vol%, balance CO2 environment. Prior ranges are rinpf (50-120), φinitial (0.15 - 0.7), Yc,intial (0.5-1), and e p (0.1-1); see Table 5.1. . . . . . . . . . . . . . . . 129 5.4 Parameter ranges after consistency analysis for US, I6, and BT chars in an O2 r = 36 vol%, H2 O = 14 vol%, balance CO2 environment. Prior ranges are rinpf (50-120), φinitial (0.15 - 0.7), Yc,intial (0.5-1), and e p (0.1-1); see Table 5.1. . . . . . . . . 130 ACKNOWLEDGEMENTS I would like to acknowledge the support of many people during my PhD program: my advisor Jennifer Spinti for her assistance in the preparation of this dissertation and for her guidance throughout this research project; Benjamin Isaac, Jeremy Thornock, Michal Hradisky, and Sean Smith for their recommendations regarding simulations of the L1500 furnace; Phillip J. Smith for his recommendations and guidance; Christopher Shaddix and Ethan Hecht for their advice during my internship at Sandia National Laboratories; Eric Eddings for his advice during the first two years of my PhD program; Jost Wendt for his willingness to be part of my committee; Dana Overacker for his assistance in the experimental work during the first year of my PhD program; John Camilo Parra Alvarez and John Marcelo Fuertez Cordoba for all their discussions and recommendations to improve my research work; Minmin Zhou and Derek Harris for their contributions to the Arches code; Bert Debusschere and coworkers for their help with the Uncertainty Quantification Toolkit; and the team involved in the L1500 experimental campaign, including Ignacio Preciado, Teri Draper Snow, and David Wagner. Finally, I would like to thank Rachel Reynolds, Pete Bringhurst, Josh Cragg, Rose Reynolds, and Yarrow Reynolds-Bringhurst for all their support and kindness during the last four years that I have lived with them. This work was supported by the Department of Energy, National Nuclear Security Administration, under Award Number DE-NA0002375. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. CHAPTER 1 INTRODUCTION The work described in this dissertation is part of the University of Utah's Carbon Capture Multidisciplinary Simulation Center (CCMSC). The main objective of CCMSC is to use verification, validation, and uncertainty quantification (VVUQ) tools to accelerate deployment of low-cost, low-emission, coal-fired power generation technologies such as Ultra Super Critical (USC) and Advanced Ultra Super Critical (AUSC) boilers. To reach this goal, we must achieve simultaneous consistency with simulation and experimental data set across a range of scales that encompass the principal physics components of these types of industrial boilers including large eddy simulations, multiphase flow, particle combustion, and radiation. We are specifically computing consistency for the quantities of interest (QOIs) in the industrial scale system(s), which are incident heat flux and tube mid-wall temperature. This hierarchical validation approach for our application is presented in Figure 1.1. There are four scales represented in the hierarchy: bench, burner, pilot, and industrial. The work presented in this dissertation focuses on two bricks in the hierarchy: the 1.5 MWth oxy-coal furnace brick in the burner-scale validation level and the char oxidation brick in the bench-scale validation level. The 1.5 MWth oxy-coal furnace is included in the hierarchical validation for two reasons. First, most of the physical phenomena present in a real boiler are also present in this furnace. Second, this furnace is located at the University of Utah and is managed by the CCMSC experimental group. As a result, I was able to obtain detailed information related to the design and operating conditions. Char oxidation is included because CCMSC researchers have performed sensitivity analysis to find what phenomena have more influence on the QOIs, and char oxidation is one of those phenomena. It is a slow process that takes place in the entire domain of a coal-fired boiler. An accurate computation of char oxidation in the entire computational domain is required to predict the particle 2 temperatures (needed for radiative heat flux calculations) and the particle reaction rate. 1.1 Validation and uncertainty quantification methodology My research program uses a validation and uncertainty quantification (VUQ) methodology based on the six-step Simulator Assessment and Validation Engine (SAVE) framework [6] coupled with a bound-to-bound consistency analysis [75, 77] to compare experimental data with simulations output at two different scales. The steps of this methodology are selection of the QOIs, construction of an input and uncertainty map, collection of experimental and simulation data, construction of surrogate models, consistency analysis of model output and feedback and/or feed-forward. I applied this methodology to the L1500 furnace and to a char oxidation model. In this VUQ methodology, the experimental data are assumed to be the true values. However, these values have uncertainty ranges that are often overlooked. In this dissertation, I present methodologies for estimating the error in the experimental data, which can then be used in the consistency analysis. 1.2 L1500 furnace The 1.5 MWth oxy-coal furnace or L1500 furnace is located at the University of Utah's Industrial Combustion and Gasification Research Facility [28, 42]. A photo of the furnace is presented in Figure 1.2. The L1500 furnace has a low-NOx burner in which the tangential component of the gas velocity (swirl) can be varied. The burner can be operated in airfired or oxy-fired modes. In oxy-fired mode, the gas composition of the oxidant stream is determined by the relative flow rates of recycled flue gas and pure oxygen (O2 ). The furnace is controlled using an OPTO 22 system; all of the operating conditions and most of the measurements are saved every 1 s. Detailed measurements of gas phase species such as O2 , CO2 , CO, NOx , and SO2 are made along the length of the furnace and at the furnace exit. Measurements related to heat flux include heat removal by cooling tubes, heat flux measured by narrow-angle radiometers, and wall temperatures. I simulated the L1500 furnace with Arches, a component of the Uintah software suite [48]. Uintah is a framework for solving partial differential equations on structured grids using hundreds to thousands of processors. In Arches, the turbulent flow is resolved 3 by the Navier-Stokes equations [66-68, 71, 82] using a large eddy simulation (LES) approach. The direct quadrature method of moments (DQMOM) [25] is used to modeled the solid (coal) phase. Char oxidation and devolatilization models provide the particle contributions to the gas phase. An ash deposition model computes the effects of ash deposition on heat transfer. Gas composition and temperature are computed with an equilibrium approach. All Arches simulations described in this dissertation employ the discrete ordinates method [59] to compute radiation. 1.3 Char oxidation The reaction process of a coal particle is shown schematically in Figure 1.3. The reaction process begins when a coal particle is injected into a hot gas environment in the presence of an oxidizer (O2 , H2 , etc.). As the particle is heated by the gas, the moisture of the particle is evaporated (drying). The particle continues heating up, releasing the volatile material in the particle (devolatilization) [24, 41, 54, 73, 77]. After devolatilization is completed, the remaining solid is known as char. The amount of char depends on the coal and on environmental conditions such as heating rate and gas temperature. The last stage of the reaction process is char oxidation [18, 31, 34, 35, 80]. The reaction process ends when all the carbon has reacted and the only solid material left in the particle is ash. I initially chose the Surface Kinetics in Porous Particles (SKIPPY) model [3, 31] and the set of experimental data collected by Hecht [34] to perform a VUQ study for char oxidation. However, I found that the assumption in SKIPPY of a steady state particle was not adequate when comparing with the set of experimental data. Hence, using SKIPPY as the starting point, I developed the transient RPBL model (described in this dissertation) that accounts for changes in physical properties with burnout and has a time-dependent boundary condition. 1.4 Research objectives This dissertation addresses the following research objectives: • Perform a VUQ analysis on the 1.5 MW (L1500) oxy-coal furnace using the Arches simulation tool for the combustion simulations and experimental data collected during the February 2015 experimental campaign. 4 • Develop a char oxidation model that computes the transient-state conditions for a spherical, reacting, porous particle and its reacting boundary layer. • Perform a VUQ analysis of the char oxidation model using experimental data from Sandia National Laboratories [34]. 1.5 Overview and organization of this dissertation I applied the six-step VUQ methodology at two scales, burner scale (L1500 furnace) and bench-scale (RPBL char oxidation model). The results I obtained are presented as Chapters 2 through 5. In Chapters 2 and 3, I report the L1500 furnace results. In Chapter 2, I provide detailed information about a simulation of the L1500 using Arches. I also present a sensitivity analysis of five parameters selected from the ash deposition and char oxidation models on three different types of QOIs: heat removal by cooling tubes, wall temperatures, and incident heat flux as measured by narrow angle radiometers. These QOIs are related to the system QOIs of incident heat flux and tube mid-wall temperature. Using the sensitivity indices, I select the two most sensitive parameters, the deposit thermal conductivity (k deposit ) and deposit emissivity (ε w ). In Chapter 3, I complete the VUQ cycle with the four remaining steps. I also explain how the uncertainty was computed for two types of experimental measurements (wall temperature and heat removal). Using Gaussian response surfaces to construct surrogate models, I perform the consistency analysis. After reviewing the analysis, I give recommendations for improving the simulation tool and the experimental data collection. In Chapters 4 and 5, I report the results from the char oxidation model. In Chapter 4, I introduce the RPBL model, the transient char oxidation model that I developed. I describe in detail all the equations used in the model and present the results from applying the model to a specific case. I also perform a sensitivity analysis on nine parameters with respect to particle temperature and velocity and select five of these parameters for the consistency analysis. In Chapter 5, I describe the experimental data set for particle temperature and velocity collected by Hecht [34] in an entrained flow reactor at Sandia National Laboratories and calculate the uncertainty in the data. I also present an Arches simulation of the Sandia reactor from which I obtain bulk gas profiles for use as the boundary conditions in the char oxidation model. I then run the 5 char oxidation model and construct polynomial chaos surrogate models from the outputs. Comparing particle temperature and velocity data with model outputs, I find consistency with most of the data. Finally, I review the analysis and give recommendations for further improvements. 6 Figure 1.1. Validation hierarchy for CCMSC. 7 Figure 1.2. Photo of the L1500 reactor. Coal par(cle Drying Char oxida(on Devola(liza(on Figure 1.3. Combustion process for a coal particle. ash CHAPTER 2 A VUQ ANALYSIS FOR THE L1500 FURNACE: SENSITIVITY ANALYSIS Submitted to the Journal of Verification, Validation and Uncertainty Quantification, A validation/uncertainty quantification analysis for a 1.5 MW oxy-coal fired furnace: Part A Sensitivity Analysis. Oscar H. Dı́az-Ibarra, Jennifer Spinti, Andrew Fry, Jeremy N. Thornock, Michal Hradisky, Sean Smith, Philip J. Smith 9 This paper focuses on a validation/uncertainty quantification (VUQ) study performed on the 1.5 MWth L1500 furnace, an oxy-coal fired facility located at the Industrial Combustion and Gasification Research Facility at the University of Utah. As part of the study, both experiments and simulations under oxy-coal combustion conditions with a swirling burner were completed. A six-step VUQ framework is used for studying the impact of model parameter uncertainty on heat flux, the quantity of interest (QOI) for the project. This paper focuses on the first two steps of the framework. The first step is the selection of model outputs in the experimental and simulation data that are related to the QOI, heat flux. In step two, an input/uncertainty (I/U) map is developed and all the parameters are assigned a priority. A sensitivity analysis is performed on five parameters in order to reduce the number of parameters that must be considered in the remaining steps of the framework. 2.1 ε qincident R ki ∆xi mc mp ms η Fp F ρ φ u ratio q Ω Nr Ir θr θ T k ∆x Io εw Nomenclature Emissivity [−] Incident radiation [Wm−2 ] Thermal resistance [Wm−2 K −1 ] Thermal conductivity for layer i [Wm−1 K −1 ] Thickness of layer i [m] Coal off gas mass flow [kg s−1 ] Primary stream mass flow [kg s−1 ] Secondary stream mass flow [kg s−1 ] Mixture fraction see Eqn. (2.3) [−] Mixture fraction see Eqn. (2.4) [−] Mixture fraction see Eqn. (2.5) [−] Gas density [kg m−3 ] Scalar Velocity in the flow direction [m s−1 ] Ratio between Arches resolution and STAR-CCM+ resolution [−] Radiative heat flux [Wm−2 ] Solid angle Number of rays Radiative intensity in each ray [Wm−2 ] Ray polar angle [degrees] View angle [degrees] Gas temperature [K] Gas absorption coefficient [m−1 ] Resolution [m] Radiative intensity from a wall [Wm−1 ] Wall emissivityl [−] 10 Tw σ Qremoval r H,l A w c W kl WH φl kc g cO,l rt Sh Re Sc dp Dom qnet Tshell ddeposit k deposit v tsb Tslag mn xi wi Wall temperature [K] Stefan Boltzmann constant 5.67037x10−8 [Wm−2 K −4 ] Heat removed by the cooling tubes [W] Volumetric reaction rate of char consumed from oxidizer l reaction [kg m−3 s−1 ] Particle surface area [m−2 ] Particle number density [# m−3 ] Mixture molar concentration [kmole m−3 ] Mixture molecular weight [kg kmole−1 ] Reaction rate coefficient for reaction l [m s−1 ] Char molecular weight [kg kmole−1 ] Stoichiometric coefficient ratio for specie l [kmolechar kmolel−1 ] Mass transfer coefficient [m s−1 ] Molar concentration of oxidizer l in the bulk [kmole m−3 ] Total volumetric reaction rate [kg m−2 s−1 ] Sherwood number [−] Particle reynolds number [−] Schmidt number [−] Particle diameter [m] Mixture averaged diffusion coefficient of oxidizer [s m−2 ] Net heat flux [Wm−2 ] External temperature [K] Deposit thickness [m] Deposit thermal conductivity [Wm−1 K] Ash deposition velocity [kg s−1 ] Soot blowing time [s] Slagging temperature [K] Moments of the Rosin-Rammler distribution Particle sizes Weights 2.2 Introduction The Carbon Capture Multidisciplinary Simulation Center (CCMSC) at the University of Utah is demonstrating the use of exascale computing with verification, validation, and uncertainty quantification as a means of accelerating deployment of low-cost, low-emission, coal-fired power generation technologies [10]. This effort employs a hierarchical validation approach to obtain simultaneous consistency among a set of selected experiments at different scales. The key physics components of a full-scale, oxy-fired boiler, including large eddy simulations, multiphase flow, particle combustion, and radiation, are represented at these various scales. The CCMSC validation hierarchy is presented in Figure 2.1. This paper focuses on validation and uncertainty quantification (VUQ) results for the 1.5 MWth oxy-coal furnace (L1500) brick in the burner-scale validation level seen in Fig- 11 ure 2.1. There are multiple reasons for performing the analysis in the L1500. First, the L1500 has multiple access points for collecting a wide variety of data, including heat flux measurements, temperature measurements, heat removal rates, and gas phase species concentrations. This type of data is difficult to collect and/or obtain at the industrial scale. Second, the L1500 is located at the University of Utah, so the validation team has access to detailed information about the geometry and operating conditions. Also, because the simulation and experimental teams work closely together, reactor and/or data collection modifications can be performed if necessary. As a result, the L1500 experimental data are optimal for VUQ analysis. Third, the L1500 sits between a laboratory-scale and an industrial-scale reactor in size. At this scale, the wall boundary conditions have significant impact on the heat flux measurements. Therefore, the facility can provide data for studies on ash deposition, heat removal by a cooled surface, and the impact of burner input conditions. Fourth, many of the complex phenomena present in a full-scale, coal-fired furnace are replicated including gas and particle phase radiation, turbulent reacting flow, coal particle devolatilization and oxidation, and particle flow. In this paper, we present the first phase of a VUQ analysis of the L1500 furnace. For this study, we desired to know the impact of the char oxidation and the ash deposition models on heat flux measurements. 2.3 Description of VUQ approach Validation is the process of testing a model using experimental data as a reference. Uncertainty quantification (UQ) is the process of computing the accuracy with which the model calculates the experimental data or quantity of interest (QOI) [63]. It is also focused on mapping uncertainty in inputs to error in the QOIs. If the model is able to reproduce the experimental data, we say that the model represents the data and thus is validated. With uncertainty quantification, we say that the model computes the QOI with specific accuracy. In reality, experiments have non-negligible bias (in addition to zero-mean noise). We almost always find it necessary to correct the bias with an instrument model. There are several VUQ methodologies that have been published. The probability bounds analysis (PBA) approach was developed by Ferson and coworkers [20-22] and applied by 12 other researchers [5, 53, 63]. The approach developed by the American Society of Mechanical Engineers (ASME) is the standard for verification and validation in computational fluid dynamics and heat transfer [4]. The bound-to-bound consistency analysis approach was developed by Frenklach and Packard [75] and used by [46, 47, 63, 77]. A six-step framework for validation of computer models called Simulator Assessment and Validation Engine (SAVE) was developed by the National Institute of Statistical Sciences group [6]. More recently, Schroeder [77] explored these VUQ tools and presented a theoretical basis for a VUQ methodology that employs the six-step SAVE framework with the bound-tobound consistency analysis [75]. This analysis employs the modified version of the SAVE framework proposed by Schroeder and is presented in Figure 2.2. In step 1, model output(s) are selected as evaluation criteria or QOIs. This step ideally involves researchers from both the simulation and experimental teams so that the QOIs can be reasonably obtained given the available facilities and instrumentation. In step 2, a list of parameters (model, scenario, and numerical) that may impact the QOIs is created and refined. This list, which also includes the parameter uncertainties, is known as the input/uncertainty (I/U) map. A determination of the impact of each parameter (e.g. priority) on the QOIs is made based on prior knowledge and/or sensitivity analysis. Parameters with high priority are selected as active and varied during the analysis while low priority parameters are fixed. Assuming that the uncertainty is a probability distribution, the uncertainty of the active parameters is then propagated through the model. Step 3 is the collection of both experimental and simulation data and the computation of the QOIs based on the data. It also includes the estimation of experimental error. Step 4, model approximation or surrogate model development, is required for expensive model evaluation. There are several types of surrogate models, including Gaussian process (GP) [6], rational and quadratic surrogate [19], and polynomial chaos (PC) [13]. In this paper, we use a PC surrogate model of first order to compute the sensitivity indices. Step 5, analysis of model outputs is performed using various methodologies. The NISS group [6] framework is partially based on the Kennedy and O'Hagan Bayesian methodology [52]. The main task is the computation of the posterior distribution, which is the 13 product of the likelihood probability distribution (assumed to be a Gaussian distribution function of the discrepancy between the model and the experimental data) and the prior distribution functions for the parameters. Another methodology is bound-to-bound consistency [75]. The basic concept of this consistency analysis is the comparison of simulation outputs with experimental data. If their difference is less than the error in the experimental measurement, the simulation outputs and experimental data are consistent. If the data are not consistent, the simulation scientist must reassess whether the models and model parameter ranges are appropriate for the system being studied, and the experimentalist must reevaluate the experimental methods and data to see if errors not previously accounted for might result in increased uncertainty. Step 6 is feedback and feed-forward. In this step, a review of the I/U map is performed and based on the results, corrections are made to the model or new capabilities are added, uncertainty in the parameters is reevaluated, and the evaluation criteria are reviewed to see if new data should be incorporated in a new VUQ cycle. This paper details how steps 1 and 2 are applied to the L1500 oxy-coal furnace brick seen in Figure 2.1. A companion paper details the completion of the remaining steps. 2.4 L1500 experiments This work focuses on data collected during a week-long experimental campaign in the L1500 furnace in February of 2015. During this campaign, the furnace was oxy-fired with a Utah bituminous coal (Sufco). 2.4.1 L1500 reactor The refractory-lined L1500 reactor, shown in Figure 2.3, is 14.2 m long with a 1.05x1.05 m2 transversal area. It is divided into 10 sections. Each section has several ports through which a variety of measurements can be taken. The reactor has eight sets of water-cooled tubes that remove heat from the first four sections (see Figure 2.4). Additionally, there is a water-cooled steel grid at the furnace exit to reduce the temperature of the combustion gases prior to entering the convection section. The L1500 can be operated in air-fired or oxy-fired modes and can burn solid, liquid, or gaseous fuels [2, 27]. During oxy-firing, recycled flue gas (RFG) is brought from the 14 exit of the convective section back into the burner's primary/secondary oxidant registers. Oxygen is supplied to the secondary and primary RFG streams just prior to entering the burner. 2.4.2 Swirl burner The swirl burner used for these tests has a range of operating modes. A schematic of this burner is presented in Figure 2.5. There are two sets of swirl vanes in the burner, one set for the inner secondary stream, and the other for the outer secondary stream. These vanes are composed of two blocks: a fixed block and a moving block. To increase the tangential velocity components (swirl), the position of the moving block is changed. In Figure 2.6, the positions of the blocks for 0% and 100 % swirl are presented. This work presents results for the 100% swirl operating condition. We performed detailed burner simulations using STAR-CCM+ for this setting. The computed swirl numbers for the inner and outer secondary registers were 2.3 and 3.6, respectively. 2.4.3 Coal characterization A Utah Sufco coal was used in this experimental campaign. The ultimate analysis for this coal is presented in Table 2.1. To determine the particle size distribution (PSD) of the Sufco coal, the bags of coal to be burned during the experimental campaign were sampled at different depths. Both a sieving analysis and a Beckman-Coulter diffraction analysis were performed on the collected samples. The collected data from both analyses are presented in Figure 2.7. From the figure, we can see that the sieving data are close to the Beckman-Coulter diffraction data if long sieving times are used. Therefore, in Figure 2.7, we fitted only data from Beckman-Coulter diffraction analysis to the Rosin-Rammler distribution. The PSD distribution was approximated with three particle sizes: 15 µm, 60 µm, and 200 µm. The mass weights, 57.4 %, 26.2 %, and 16.4 %, were computed with the Eqn. (2.1). 2 ∑ xin wi = mn f or n = 0, 1, 2 (2.1) i =0 In this equation, mn are the moments of the Rosin-Rammler distribution presented in Figure 2.7. 15 2.4.4 Operating conditions The L1500 operating parameters for this VUQ study are shown in Table 2.2. The coal was fed with the primary stream. With the exception of the burner swirl setting (100 %), all of these parameters were controlled and recorded during the experimental campaign. These inputs were computed from the mass flow of the RFG to the primary and secondary registers, the oxygen mass flow to the primary and secondary inlets, the coal mass flow rate, the gas temperatures for each of these streams, and the composition of the RFG. These operating conditions were stable during the experiment. However, there was an air leak in the RFG stream as evidenced by the outlet CO2 concentration being lower than expected. To compensate for this air leak in the simulations, an overall mass balance was performed on the furnace. The leaked air and the coal moisture are included in the RFG composition shown in Table 2.2. Shell temperatures were measured during the experimental campaign using a surface thermocouple. These temperatures varied by approximately 100 K depending on location. In scoping simulations, we determined that simulation results (specifically radiative heat flux and wall temperatures) were not sensitive to the shell temperature of the furnace within the range of temperatures measured experimentally. Therefore, the averaged experimental values shown in Table 2.3 were used for each of the four sides of the furnace. The thermal properties of the L1500 walls are presented in Table 2.4. With these properties, the total thermal resistance was computed using Eqn. (2.2). The average total thermal resistance for the wall is R = 1.02 mW2 K . R= 2.5 1 ∆x ∑4i=1 ki i (2.2) L1500 simulation strategy The simulation data for the 100 % (burner) swirl experiment were obtained from the coupling of two Large Eddy Simulations (LES), one of the burner and the other of the main chamber (see Figure 2.8). Because the burner is equipped with swirl blocks and has a complex geometry, the commercial software package STAR-CCM+ was used to build a mesh and perform numerical simulations of the burner. The time averaged velocity profile at the burner tip was then extracted, filtered, and used as the boundary condition 16 (i.e. handoff file) for the simulation of the main chamber using Arches, an LES research code developed at the University of Utah. We did not perform a simulation of the entire reactor in STAR-CCM+ due to the lack of coal particle physics models (e.g. char oxidation, devolatilization, ash deposition, direct quadrature method of moments) required in our analysis. 2.5.1 STAR-CCM+ simulations We resolved the complex geometry of the swirl burner in STAR-CCM+ with an unstructured mesh. We used an LES approach for turbulence, and the subgrid scales were modeled with the WALE model [62]. To model the fluid flow in the boundary layer, we use cell sizes as small as 15 µm to fully resolve the viscous sublayer (with y+ values of approximately one). In the energy equation, the properties of the gases were computed assuming nonreacting, multicomponent ideal gas mixtures. We also assumed the walls of the burner were adiabatic. Only the gas phase was resolved; the particle phase was not included in these simulations because we were only interested in the gas velocity profiles at the tip of the burner. We assumed that the particles would not affect the gas velocity as there were no chemical reaction inside the burner and the particle loading was low in comparison to the primary gas feed rate. To solve the equations, we used a second-order implicit time solver with a fixed time step of 5x10−5 s. The burner simulations were run on a Linux cluster at the University of Utah. Each simulation required 3,600 cores for two weeks, which was long enough to obtain 3 s of physical time. Data were averaged over an interval of 1 s to obtain the velocity profile. This interval was much shorter than the averaging interval for L1500 furnace data because of the short residence time/high velocity through the burner. 2.5.2 Arches simulations The simulations of the L1500's main chamber were performed with Arches, a component of the Uintah software suite [48]. Uintah is a framework for solving partial differential equations on structured grids using hundreds to thousands of processors. Arches is a low-Mach, pressure projection-based, variable density code solving filtered (per the standard large eddy simulation approach) conservation equations for mass, momentum, and energy of the gas and solid phases [66-68]. For the gas phase, subgrid fluc- 17 tuations in the momentum equation are modeled using a dynamic Smagorinsky closure model. Convective terms in the momentum transport are discretized with a hybrid scheme (a blend between first-order upwind and central differencing) [92]. A first-order, explicit forward Euler approach is used to advance the equations forward in time. A pressure projection ensures that continuity is strictly satisfied for each time step. The solid (coal) phase was modeled using the direct quadrature method of moments (DQMOM) with three environments [25]. Seven internal coordinates, including raw coal mass, the char mass, the particle enthalpy, the maximum particle temperature, and the particle velocity (in x, y, and z coordinates), were used to describe the PSD. In addition, the weights and internal coordinates for the three environments were solved using a first order upwind discretization scheme. The particle phase was coupled with the gas phase through source terms in the equations for momentum, heat, and mass. More detail about the DQMOM implementation in Arches is presented in [66-68]. The coal particle physics models included both devolatilization and char oxidation models. The devolatilization model is detailed in [77]; it is based on the Chemical Percolation Devolatilization (CPD) model [24]. It assumes that the rate of devolatilization is a first order reaction rate. The ultimate volatile yield is a function of the particle temperature and the parameters in the ultimate volatile yield equation are fitted with the CPD model. The char oxidation model is discussed in section 2.5.2.1. The gas-phase reactions in the system were modeled using a three-stream mixture fraction approach. The three streams were the primary stream (m p ), the secondary (inner + outer) stream (ms ), and the coal off gas (mc ). During char oxidation or devolatilization coal off gas enters the gas phase with an elemental composition shown in Table 2.1. Gas composition and temperature are computed assuming equilibrium (with constant enthalpy and pressure). The mixture fractions based on these three streams are defined in Eqn. (2.3), Eqn. (2.4), and Eqn. (2.5). η= mc m p + ms + mc (2.3) Fp = mp m p + ms + mc (2.4) mp m p + ms (2.5) F= 18 Transport equations were solved for η and Fp ; F was computed from the other two mass fractions as 1 − Fp − η. A lookup table based on equilibrium chemistry assumptions was tabulated a priori as a function of three independent variables: η, Fp , and the linearized enthalpy [84]. The state space variables in the table include the gas temperature, mixture enthalpy, species mass fractions, mixture molecular weight, and mixture density. The discrete ordinates (DO) method was used to compute radiation in the simulation [59]. For the cases in this paper, S8 quadrature, representing 80 discrete directions, was used. The L1500 simulations were run on two Linux clusters, one at the University of Utah and the other at Lawrence Livermore National Laboratory. With a 15 mm mesh resolution, the computational domain had 2,255,610 cells. The time step was ∼ 6.5e − 05 [s]. The simulations were run on 217 cores for approximately 100 hours. This was long enough to obtain 30 s of physical time; on average, ∼ 15 s were required to reach steady state. We averaged the last 5 s of simulation data for our analysis. In this study, our focus is on the impact of the char oxidation and the ash deposition models on the QOIs. Therefore, we describe these model in the following sections. 2.5.2.1 Char oxidation The char oxidation model implemented in Arches includes heterogeneous reaction at the particle surface, mass transport of oxidizer from the bulk gas to the particle surface, and mass transport of devolatilization and oxidation products away from the particle surface. It computes the volumetric reaction rate of char consumed by the oxidizer in global reaction l using Eqn. (2.6). g r H,l ( Aw)2 cWk l WH φl k c cO,l = AwcW (k l + k c ) + rt (2.6) For the L1500 simulations, two global reactions were considered: one oxidation and one gasification reaction. The oxidation reaction with O2 as oxidizer is Eqn. (2.7). The reaction − El constant in Arrhenius form is k l = Al e RT . The two parameters in this equation are AO2 and EO2 . The gasification reaction with CO2 as oxidizer is Eqn. (2.8). As with the oxidation reaction, the reaction constant k l has two parameters, ACO2 and ECO2 . The mass transfer coefficient k c was obtained using a Sherwood number correlation with a correction factor as shown in Eqn. (2.9). 19 Cs + O2 → CO2 (2.7) Cs + CO2 → 2CO (2.8) Sh = 2 + 0.6Re1/2 Sc1/3 = 2.5.2.2 kc d p Dom (2.9) Ash deposition Arches uses a one-dimensional wall boundary condition as shown in Eqn. (2.10). In this equation, Tw is the internal (hot face) temperature of the wall and Tshell is the external (cold face) temperature. For the refractory wall, Tshell is the outside wall temperature and for the cooling tubes, it is the water temperature inside the tubes. R is the thermal resistance, ε is the emissivity of the ash deposit, and qincident is the incident heat flux. In this model, Tw is solved with a Newton solver every time step. The thermal resistance R is computed with Eqn. (2.11). In this equation, the first term is the resistance produced by the refractory and insulation layers in the furnace wall or by the metal in the cooling tubes. The second term is the resistance caused by the ash layer. The thermal conductivity of the deposit, k deposit , is an input parameter. qnet = ( Tw − Tshell ) = ε(qincident − σTw4 ) R Nlayer R= ∑ i =1 ddeposit ∆xi + ki k deposit (2.10) (2.11) This deposition model uses three regimes. In regime 1, Tw is lower than the ash fusion temperature, Tslag , and the deposit thickness, ddeposit , is computed with Eqn. (2.12). In this equation, the ash deposition velocity (v) is computed from a probability model [8] and tsb , the time since the last soot blowing event, is an input parameter. In regimes 2 and 3, Tw is equal to or greater than Tslag and ddeposit is computed with Eqn. (2.13). If the computed value of ddeposit is greater than zero, the model is in regime 2 and Tw = Tslag . If the computed value of ddeposit is less than zero, the model is in regime 3. In this regime, ddeposit = 0 and Tw is computed from Eqn. (2.10). ddeposit = vtsb ddeposit = k deposit Tslag − Tshell qnet max − Rw (2.12) (2.13) To compute ddeposit in Eqn. (2.13), qnet max is computed from Eqn. (2.10) with Tw replaced N with Tslag and Rw replaced with ∑i=layer 1 ∆xi ki . 20 Finally, Arches has a basic erosion model. If ddeposit computed with Eqn. (2.12) or Eqn. (2.13) is greater than the erosion thickness, ddeposit is set equal to this value. Thus, erosion thickness is the maximum deposit thickness allowed in the simulation. 2.5.3 Handoff plane boundary condition Because of the difference in grid size between the STAR-CCM+ and Arches simulations, the velocity field computed in STAR-CCM+ must be filtered for use in Arches. This filtering process has six steps. 1. Extract the velocity vectors from STAR-CCM+ on a structured grid using nearest neighbor interpolation. The number of points in the extracted grid depends on the desired ratio (an integer) between the STAR-CCM+ and Arches resolutions (here this value is 30) and the final resolution of the Arches simulations (15x10−3 m). These choices result in a structured grid of 419x419 points with a resolution of 5x10−4 m. 2. Map the STAR-CCM+ information to a 2-D array, adding zeros to the array where data were not extracted (outside of the burner). 3. Create a 2-D array for the fraction of the primary stream, Fp , with the condition Fp = 1 (Eqn. (2.4)) for r ≤ 0.0387m (radius of the primary inlet), and Fp = 0 for r > 0.0387m. The mixture fraction field for coal off gas, η (see Eqn. (2.3)), is zero at the inlet because we assumed there were no coal particles reactions in the burner. 4. Obtain density, ρ, from the equilibrium chemistry lookup table described above using Fp and assuming the burner is adiabatic. 5. Compute the components of the velocity in the flow (u) and the tangential (v and w) directions. To do this, the mass flux at the inlet, ρu (units of kg ), m2 s is filtered using Eqn. (2.14), where φ = 1 for the u component, φ = v for the v component, φ = w for the w component, and ρu is the product of the gas density and velocity in the flow direction extracted from STAR-CCM+. In order to obtain u, ρu is divided by ρ. For v and w, ρuv and ρuw are divided by ρu. ( j+1)ratio (i +1)ratio (ρuφ)i,j = (ρuφ)m,n ratio2 n= j∗ratio m=i ∗ratio ∑ ∑ (2.14) 21 6. Write the filtered velocities u, v, and w to an input (handoff) file that is read by Arches. This entire process is illustrated by the velocity fields shown in Figure 2.9. On the left are the three velocity components from the STAR-CCM+ simulation at the exit plane of the burner. On the right are the filtered velocity components used in the Arches simulation. The particle inlet conditions are computed separately because the STAR-CCM+ results are only for the gas phase. For the DQMOM inlet conditions, the velocity of the particles is assumed equal to the gas velocity and a constant coal mass flow rate is assumed for each cell. To compute the flux of particles at the burner tip, we assume that the particle velocities were the same as the gas velocity. These mass weights were converted to particles per cubic meter using a coal density of 1300 kg/m3 and the volume corresponding to each particle size. While the assumed density of the particle will affect the number of particles per cubic meter, the flow rate of coal is independent of this variable. 2.5.4 Domain size and mesh resolution In previous simulations of the entire furnace geometry at a resolution of 16 mm for the 0% swirl case, we determined that gas temperatures, velocities, and chemical compositions were relatively constant after section 6 [17]. Also, differences between simulation results from the full length (14.2 m) and from simulations using a shortened domain ( 7 m) were minimal. Hence, the computational geometry for the analysis that follows is 7 m. The shortened geometry of the L1500 with a 15 mm resolution is presented in Figure 2.10; this geometry includes the burner quarl, the cooling tubes, and the 10 cm step up in the bottom of the furnace between sections 4 and 5. The surface area of the cooling tubes in the computational mesh was adjusted to match the actual surface area of the tubes. 2.6 Analysis of quantities of interest To perform a VUQ analysis, the QOIs and the system parameters (scenario, model, numerical) that have a first order impact on the QOIs are identified. In this experimental data set, the QOIs all relate to heat flux. The QOIs are heat flux measurements from three narrow angle radiometers (sections 1, 2, and 3), five wall temperature measurements 22 (sections 1, 2, 3, 4, 6), and heat removal by eight sets of cooling tubes. During the experimental campaign, radiative heat flux was measured through the center port of sections 1, 2, and 3 using a narrow angle radiometer. A cold plate was installed in the port opposite to the radiometer to provide a known boundary condition for the radiometer measurements. In practice, the cold surface became coated with radiating particles, introducing uncertainty into the radiometer measurement. Wall temperature measurements were taken in sections 1, 2, 3, 4, 6, and 8 using Type B thermocouples encased in ceramic sheaths that were then inserted into small holes in the furnace ceiling located in the middle of each section (see Figure 2.11). Each sheath was inserted until it was flush with the inside wall of the furnace. The heat removal by the cooling tubes was determined by the experimentally-measured mass flow rates and inlet/outlet temperatures of the water flowing through the eight cooling tubes in the first four sections. 2.7 I/U map for char oxidation and ash deposition models Our focus in this study is on the impact of the char oxidation and the ash deposition models on the QOIs. Model parameters and their associated uncertainty ranges are presented in Table 2.5. In this table, the overall priority means the relative importance of that parameter on the QOIs. Because of the computational cost of an LES simulation, we must select the most important parameters for further study in the sensitivity and consistency analyses. In the ideal case, a sensitivity analysis would be performed to assign a priority value to the parameters in the I/U map. However, running a sensitivity analysis with 11 parameters is too expensive for an LES simulation, even with the L1500 coarse mesh (15 mm resolution) case described above. If we use two simulations for each dimension and combine all the parameters, a total of 211 = 2048 simulations would be needed. Consequently, we carried out a preselection of parameters to assign priorities. In Table 2.5, the most important parameters are assigned a priority of 6 and the less important a priority of 3. In the char oxidation model, the parameters are the kinetic constants for the two reactions, an oxidation reaction (Eqn. (2.7)) and a gasification reaction (Eqn. (2.8)). Using 23 prior knowledge about char oxidation, we assigned a priority of 6 to the oxidation reaction parameters and a priority of 3 to the gasification reaction parameters. The erosion thickness, the maximum deposit thickness allowed in the simulation, was given a priority of 3 because we set this value high enough that there was no limitation on the deposit thickness (1 cm on the cooling tubes and 5 cm on other surfaces) in the time frame of the simulation. It is difficult to assign a value to the soot blowing time (tsb ) because there was no soot blower in the L1500 and different conditions were tested throughout the week of experiments. For these reasons, a priority of 6 was assigned to this parameter. The shell temperature (Tshell ) is the exterior shell temperature for the L1500 wall and the inside tube temperature for the cooling tubes. This parameter was given a priority of 3 because (1) the water temperature did not change much in the cooling tubes between the inlets and outlets and (2) a one-dimensional heat transfer analysis through the furnace wall showed that Tshell had a negligible impact on the internal furnace temperature. N The thermal resistance (∑i=layer 1 ∆xi ki ) due to the refractory and the insulation material was computed using an average thermal conductivity obtained from the manufacturer's data. Because it could be well-characterized, a priority of 3 was assigned to this parameter. The ash density was kept constant for the analysis. The effective thermal conductivity k deposit and temperature parameter Tslag were analyzed in a VUQ study of a 15 MWth oxy-coal boiler [45] . This study concluded that the Tslag parameter was well-represented by the ash fusion temperature. Therefore, a priority of 3 was given to Tslag , and the ash fusion temperature for Sufco coal was used as the nominal value. For the k deposit parameter, a priority of 6 was assigned because of the lack of information about the thermal conductivity of the ash deposits on the cooling tubes. The ε w was given a priority of 6 because there is almost no information about this parameter in this reactor. The active parameters for this analysis were the parameters with a priority of 6. This set includes two parameters related to char oxidation, AO2 and EO2 , and three parameters related to ash deposition, k deposit , ε w , and ts . The other parameters in Table 2.5 were fixed at the nominal value in the simulations. The uncertainty ranges for the char oxidation parameters AO2 and EO2 , which are the reaction parameters of the oxidation reaction (Eqn. (2.7)), were taken from Smoot and Smith [83]. Smoot and Smith reported small variations (17-20 kcal/mole) in activation 24 energy, EO2 , for three U.S. coals: Wyoming subbituminous, Illinois No. 6 bituminous, and Pittsburgh No. 8 bituminous. For AO2 , the range of 58-68 g cm−2 s−1 atm−1 O2 was selected based on the reported value of 60 g cm−2 s−1 atm−1 O2 for Illinois No. 6 (similar coal type and composition). For the ash deposition model, the range for k deposit was bounded by the maximum and minimum values of the experimental data presented by Rezaei and coworkers [72] for coal ash and synthetic ash. The ε w range was bounded by the maximum and minimum values reported by Wall and coworkers [87, 94] for coal ash. The parameter tsb is the interval of time without ash removal. Since the L1500 did not have a soot blower at the time of the experimental campaign, this value was computed as the number of hours that coal was fed to the reactor. The maximum value was the total coal feed time for the campaign (at night, the furnace was fired with natural gas), and the minimum value was the total coal feed time for the campaign less the first day of experiments when the coal feed rate was lower. 2.8 Sensitivity analysis In the I/U map (Table 2.5), we gave five variables a priority of 6, which means these variables are active in this study. However, a full five-dimensional VUQ study with Arches simulations is too computationally expensive. Therefore, a sensitivity analysis was used to reduce the number of dimensions in this study from five to two for the consistency analysis. In order to perform the sensitivity analysis, each QOI was computed from simulation data as described in the following sections. 2.8.1 Computation of radiometer heat flux Three narrow angle radiometers were used to measure the radiative heat flux in sections 1, 2, and 3 in the L1500 experiments. To compute these heat fluxes from the simulation data, we used a reverse-Monte Carlo ray tracing approach that sums up the radiative intensities over all the rays comprising the field of view, θ, of the radiometer as seen in Eqn. (2.15). The solid angle Ω was computed using Eqn. (2.16). The radiative intensity in each ray, Ir , was computed with Eqn. (2.17). For this analysis, Nr = 1; because of the 25 coarseness of the computational mesh, one ray was sufficient to account for the narrow angle view of the radiometer [38]. In order to compute the radiative intensity with Eqn. (2.17), the gas temperature T, the gas absorption coefficient k, and the mesh resolution ∆x were obtained from the simulations. These values were extracted along one ray that extended from the wall opposite the radiometer to the wall where the radiometer tip was located. Ash buildup on the cooled targets flush with the walls opposite the radiometers resulted in a surface condition that was unknown. To compute Io , the intensity of the target, simulation values for wall temperature, Tw , and wall emissivity, ε w , were used; see Eqn. (2.18). q=Ω 1 Nr Nr ∑ Ir cos(θr ) Ω = 2π (1 − cos(θ )) Ir = N −1 ∑ j =0 (2.15) r =1 j −1 σ 4 T exp(− ∑ ∆xk i )(1 − exp(∆xk j )) π j i =0 (2.16) (2.17) N + Io exp(− ∑ ∆xk i ) i =0 Io = ε w σ 4 T π w (2.18) More detailed information about computation of the radiometer heat flux is found in [37, 38]. 2.8.2 Computation of heat removal by cooling tubes For each computational cell in the simulated cooling tubes, the heat removal (Qremoval ) was computed with Eqn. (2.19) and then added over the whole cooling tube to obtain the total heat removal. Thus, in Eqn. (2.19), A corresponds to the surface area of the cell faces that are in contact with the gas, Tw is the surface temperature of the cooling tube, qincident is the incident heat flux, and ε is the emissivity of the tube. N Qremoval = ∑ ε(qincident − σTw4 ) A (2.19) n =0 This computation was performed using the visualization tool VisIt [11]. The heat removal was calculated using the average fields of qincident and Tw . 26 2.8.3 Computation of wall temperature For the sensitivity analysis described below, we used the wall temperature computed in the simulation at the thermocouple location. 2.8.4 Generating PC surrogate models The sensitivity analysis is done using the Uncertainty Quantification Toolkit (UQTk), a set of C++ tools with a Python interface. It was developed by Debusschere and coworkers at Sandia National Laboratories [13]. UQTk uses the app pce_sens to compute total and main sensitivity indices using a polynomial chaos (PC) surrogate model. Using a PC surrogate model of order 1 with a full quadrature rule, a total of 32 Arches simulations (2n where n is the number of dimensions) were needed for this study. The UQTk app generate_quad generated 32 quadrature points of ξ i = +/ − 0.58 with weights of w = 0.0625 for each dimension. The variable ξ i is mapped to physical input space for the Arches simulations using Eqn. (2.20), where ai and bi are the bounds in the uncertainty interval presented in the I/U map (Table 2.5). λi = b − ai a i + bi + i ξi 2 2 for i = 1, ..., d (2.20) We ran a set of Arches simulations at the 32 quadrature points. All the simulations were run for approximately 30 s. Only the last 5 s of data were averaged for this analysis; simulations required at least 15 s to reach steady state, so we had a buffer of 10 s to ensure that steady state had been reached for all simulations. We then computed the radiative heat flux, heat removal, and wall temperatures and generated the PC surrogate model for each of these data points. 2.8.5 Results Our next step was to compute the main and the total sensitivity indices with the UQTk app pce_sens using the coefficients of the PC surrogate model. In Figure 2.12, the main sensitivity indices for the radiative heat flux are presented. The two most sensitive parameters with respect to the radiometer measurements are ε w and k deposit ; both parameters are part of the ash deposition model. It is interesting that the main sensitivity index for tsb is low by comparison. The main sensitivity indices for the AO2 and EO2 are also low, meaning these parameters have little influence on the radiometer measurement. Thus, 27 for the radiative heat flux measurement, the five-dimensional (parameter) study can be reduced to a two-dimensional study. Figure 2.13 shows the main sensitivity indices for the heat removed by the cooling tubes. The sensitivity indices have a similar behavior to those obtained in the radiometer analysis. The ε w and k deposit parameters have bigger indices while the sensitivity index for tsb is small for the cooling tubes. The sensitivities of the char oxidation parameters are small as well. Just as with the radiometer measurements, this analysis shows that for the heat removal measurements, the study can be reduced from five to two dimensions. In Figure 2.14, the main sensitivity indices for the wall temperature at the thermocouple locations are presented. For all locations, the biggest sensitivity index is for ε w . While the sensitivity index for k deposit is low at location 1 (WT1 ), its index values are much higher at the other locations. Thus, we conclude that the consistency analysis study can be reduced from five to two dimensions. The insensitivity to tsb (which determines deposit thickness) for all QOIs was unexpected, so the simulations were checked to look for an explanation. For all 32 simulations, the deposit thickness was zero for all the walls except for the cooling tubes after 30 s of simulation time. There are two reasons for this result. First, in all the simulations, the wall temperatures were higher than Tslag in sections 1 and 2, and for some simulations, wall temperatures exceeded Tslag for the entire length (7m) of the simulated reactor as presented in Figure 2.15. Consequently, at least the first two sections of the L1500 are in regime 3 of the ash deposition model where the ash deposit thickness is zero, meaning the deposit is slagging or evaporating due to high temperatures. Second, for the 100 % swirl operating condition, most of the ash deposition is expected in sections 1 and 2 where the model is in regime 3. The influence of the cooling tubes deposits can be seen in the sensitivity index of the radiometers, where k deposit has the second largest effect after ε w and where its index increases with the section (see Figure 2.12). This result means that the heat removal by the cooling tubes is being impacted by the ash deposition, which in turn impacts the gas temperature and thus the computation of the radiative heat flux. This impact of ash deposition is also seen in Figure 2.13. In sections 1 and 4, HR1 , HR2 , HR7 , and HR8 are the only QOIs where the sensitivity indeces of tsb are greater than those of k deposit . In section 28 2 (HR3 and HR4 ), k deposit has a bigger sensitivity index than ε w . Finally, in section 3 (HR5 and HR6 ), k deposit still has a strong influence. 2.9 Conclusions The work described in this paper represents the first part of a VUQ study of simulation and experimental measurements taken in the oxy-fired L1500 furnace at the University of Utah. For the experiments described in this paper, the burner was set at 100% swirl. The simulations performed in this study utilized a handoff strategy wherein the complex flow through the burner was computed using STAR-CCM+. Averaged scalar and velocity fields at the plane of the burner tip were then provided as inputs to the Arches simulations of the L1500. We discuss the first two steps of the VUQ analysis in this paper: (1) selection of QOIs and (2) creation of an I/U map. We perform a sensitivity analysis to reduce the number of active parameters. We conclude that of the five active parameters in the I/U map, the most sensitive parameters across three different types of measurements (heat removal by cooling loops, heat flux measured by radiometers, and wall temperature) are k deposit and ε w . Therefore, we will use these two parameters for the consistency analysis, as presented in part B. 2.10 Acknowledgment This material is based upon work supported by the Department of Energy, National Nuclear Security Administration, under Award Number DE-NA0002375. We would like to acknowledge the following groups and resources for their support of and contributions to this project: first, the Arches development team, especially Derek Harris and Minmin Zhou for their contributions to the Arches code; second, the team involved in the L1500 experimental campaign, including Ignacio Preciado, Teri Draper Snow, Eric Eddings, and David Wagner; third, the Center for High Performance Computing at the University of Utah and Lawrence Livermore National Laboratory for computing resources; fourth, Oscar Dı́az-Ibarra would like to thank Dr. Debusschere and coworkers for their help with the Uncertainty Quantification Toolkit. 29 Table 2.1. Ultimate analysis for Utah Sufco coal. Data are averaged and normalized from the analysis of five samples taken from the bags of coal that were burned during the testing. Sufco Coal C H N S O Ash H2 O HHV [ J/kg] % Mass 66.89 4.51 1.17 0.36 13.60 7.89 5.58 27364.93 Table 2.2. L1500 operating conditions for 100 % swirl Stream Mass flow [kg/s] Coal (ash and 0.03534 moisture free) Primary (m p ) 0.07103 Inner Secondary 0.05899 (ms ) Outer Secondary 0.23515 (ms ) Mass Fraction Primary (m p ) O2 CO2 H2 O SO2 N2 0.1684 0.6464 0.1314 0.0009 0.0529 Temperature [K ] 338 338 533 533 Inner and Outer Secondary (ms ) 0.2412 0.6114 0.0965 0.0009 0.0500 30 Table 2.3. Shell temperatures averaged over all measurements made on a side Location Shell Temperature [K] Quarl 434 Main chamber south side 362 Main camber north side 396 Main chamber bottom side 362 Main chamber top side 427 Table 2.4. Thermal properties obtained from the manufacturer for the L1500 walls W W Layer name ∆x [m] k [ mK ] σ [ mK ] Temperature range [K] Greencast0.2032 2.36 0.52 1143 94 Plus 1588 Insboard 0.0508 0.19 0.10 698 3000 1477 Insboard 0.0254 0.1475 0.0562 589-1255 2600 Insblock 0.0508 0.104 0.032 575-923 31 Table 2.5. Input/uncertainty map for char oxidation and deposition model Parameter Priority Range Nominal Basis/Comments value min max Char oxidation AO2 6 58 68 Pre-exponential factor for oxidation reaction Eqn. (2.7). [g cm−2 s−1 atm−1 O2 ] EO2 6 17 20 Activation energy for oxidation reaction [kcal/mole] ACO2 3 1390 Pre-exponential factor for gasification reaction Eqn. (2.8) [g cm−2 s−1 atm−1 O2 ] ECO2 3 53700 Activation energy for gasification reaction [kcal/mole] Ash deposition Erosion thick- 3 1 or 5 Maximum thickness: 1 cm for coolness ing tubes, 5 cm for everything else tsb 6 120666 137322 Soot blowing time or time since ash was removed from cooling tubes [s] Tshell 3 Average temperature inside cooling tube or external surface temperature of L1500 Nlayer ∆xi 3 1.02 Thermal resistance in the L1500 ∑ i =1 k i walls [Wm−2 s−1 ] Tslag 3 1516.5 Ash fusion temperature [K] for Sufco coal k deposit 6 0.1 1.8 Thermal conductivity of deposit [W m−1 K −1 ] εw 6 0.3 1 Wall emissivity [-] 32 Figure 2.1. Validation hierarchy for CCMSC 33 Figure 2.2. Six-step methodology with consistency analysis. The yellow box corresponds to the steps analyzed in this paper. 34 s ction 10 se a) To convection zone Burner Sample ports Cooling tubes b) Burner Sections 1 2 3 4 Figure 2.3. Drawing of the L1500 reactor located at the Industrial Combustion and Gasification Research Facility. a) Schematic of the L1500 multifuel combustion furnace. b) Schematic of the first four sections of the L1500 furnace. 8 6 7 4 5 3 2 Item No. Qty. Part No. 1 2 Section 1 Cooling Panel D a) 1 Description Material 1/2" Sch. 40 Pipe 316 SS D All Bends are 2.5" in Radius. C C b) B B 32" 16" 5" A 5" 5" 5" 5" 5" 6" 5" 7 MADE BY CK'D BY Andrew Fry 06/25/14 6" University of Utah, ICGRF 6" 870 South 500 West, Salt Lake City, UT 84101 801 587-1781 Section 1 Cooling Panel 26" 8 REVISION DATE 5" Fabricate 2 6 5 4 SIZE DWG. NO. REV. A 3 SCALE 1:12 CAD FILE: 2 Figure 2.4. L1500 cooling tubes. a) Cooling tubes are located in the first four sections of the furnace. Photo was taken from the burner location looking toward the furnace exit. b) Dimensions of a set of cooling tubes (in inches). SHEET 1 1 of 1 A 35 Primary and coal Inner Secondary Outer Secondary Figure 2.5. Schematic of the swirl burner used in the L1500 furnace. Inlets are labeled. 0 % swirl 100 % swirl Figure 2.6. Swirl vanes in the burner at the 0% swirl and 100% swirl positions. Figure 2.7. Experimental PSD and fitted Rosin-Rammler PSD. Sieving was performed for different lengths of time to determine the effect on the resulting PSD. Based on this analysis, only the Beckman-Coulter diffraction data were used for the Rosin-Rammler fit. 36 STAR-CCM+ ARCHES Burner Chamber Figure 2.8. L1500 simulation coupling between Arches and STAR-CCM+ Figure 2.9. Velocities at plane of the burner tip; left are STAR-CCM+ results with a resolution of 0.5x10−3 m, right are Arches results with a resolution of 15x10−3 m (ratio = 30). Cooling tubes Quarl Section 1 Section 2 Section 3 Step up Section 4 Section 5 Section 6 Figure 2.10. Shortened geometry for the L1500 simulations including the quarl, the eight sets of cooling tubes, and the step change in the reactor floor. Resolution is 15 mm. 37 Sheath Furnace Wall Insulation Material Section 4 Location of thermocouple measurement Gas Figure 2.11. Thermocouple placement in furnace wall. Section 4 is shown but the placement is similar for all thermocouple measurements. 1.2 Ao2 EO2 kdeposit tsb εw 1.0 Sensitivity 0.8 0.6 0.4 0.2 0.0 rad1 rad2 rad3 Figure 2.12. Main sensitivity index for radiative heat flux measured by the radiometers. 38 1.2 Ao2 EO2 tsb εw kdeposit 1.0 Sensitivity 0.8 0.6 0.4 0.2 0.0 HR 1 HR2 HR3 HR4 HR5 HR6 HR7 HR8 Figure 2.13. Main sensitivity index for the heat removed by the cooling tubes. 1.2 Ao2 EO2 tsb εw kdeposit 1.0 Sensitivity 0.8 0.6 0.4 0.2 0.0 WT1 WT2 WT3 WT4 WT5 Figure 2.14. Main sensitivity index for wall temperature (thermocouple's location) computed with Arches. 39 1900 Tslag Simulation Wall Temperature [K] 1800 1700 1600 1500 1400 1300 1200 0 1 2 3 4 Section No 5 6 7 Figure 2.15. Wall temperature at thermocouple locations computed with Arches. Horizontal red line is Tslag . CHAPTER 3 A VUQ ANALYSIS FOR THE L1500 FURNACE: CONSISTENCY ANALYSIS Submitted to the Journal of Verification, Validation and Uncertainty Quantification, A validation/uncertainty quantification analysis for a 1.5 MW oxy-coal fired furnace: Part B Consistency Analysis. Oscar H. Dı́az-Ibarra, Jennifer Spinti, Andrew Fry, Jeremy N. Thornock, Michal Hradisky, Sean Smith, Philip J. Smith 41 This paper focuses on a validation and uncertainty quantification (VUQ) study performed on the 1.5 MWth L1500 furnace, an oxy-coal fired facility located at the Industrial Combustion and Gasification Research Facility at the University of Utah. We performed this analysis with three types of data or quantities of interest (QOI) collected in the L1500 furnace: thermocouple temperature measurements at the wall, heat removal by cooling tubes, and oxygen molar fraction at the furnace exit. We estimated sampling error with raw experimental data and systematic error with instrument models for each type of measurement. We ran 34 large eddy simulation (LES) cases exploring the effect of three parameters on the QOIs: effective thermal conductivity of the ash, a burner swirl parameter, and the coal feed rate. From these cases, we created Gaussian response surface surrogate models that we used to perform a consistency analysis with the three QOIs. We found consistency between the experimental QOIs and the corresponding LES simulation data and determined that the burner inlet condition (swirl) had a significant impact on the QOIs. 3.1 b̄ bi cp error fs F HR k ṁ m̄ mi N qincident Qremoval T ue x ȳ yT β δ ∆ Nomenclature kg Mean intercept in calibration curves[ s ] Slope for a calibration curve i J ] Heat capacity [ kgK Error in experimental measurement Swirl parameter [−] Net heat flux [ mW2 ] Experimental heat removal [W] W Thermal conductivity [ mK ] kg lb Mass flow rate [ s or h ] kg Mean slope in calibration curves [ Hz s ] Slope for a calibration curve i Number of calibration curves Incident heat flux [ mW2 ] Simulation heat removal [W] Temperature [K] Experimental uncertainty Number of revolutions [Hz] Measured experimental value True value of measured variable Sampling or random error Systematic error Normalized discrepancy 42 ∆x ∆y ε ρ σ φhole Grid spacing in x direction [m] Grid spacing in y direction [m] Emissivity [−] kg Density [ m3 ] Stefan Boltzmann constant 5.670367x10−8 [Wm−2 K −4 ] Diameter of the thermocouple hole [m] Subscripts b: coal: deposit: ef f: i, j: in: ins: inter f ace : l: lower: out: shell: re f : r: std: sur f ace: tc: t: upper: w: water: 1: 2: Bottom Coal Ash deposit Effective Computational cell i,j Inlet Insulation Interface between thermocouple and wall Left Lower bound Outlet Outside wall Refractory Right Standard deviation Surface of the internal wall Thermocouple Top Upper bound Wall Liquid water Refractory layer Insulation layer 3.2 Introduction The Carbon Capture Multidisciplinary Simulation Center (CCMSC) at the University of Utah is demonstrating the use of exascale computing with verification, validation, and uncertainty quantification (VUQ) as a means of accelerating deployment of low-cost, lowemission, coal-fired power generation technologies [10]. This effort employs a hierarchical validation approach, presented in Figure 3.1, to obtain simultaneous consistency among a set of selected experiments and simulations at different scales with the ultimate objective of achieving predictive capability for a full-scale, coal-fired boiler, where experimental measurements are expensive and difficult to obtain. The key physics components in these systems include multiphase flow, turbulence, particle combustion, radiation, as well as the 43 numerical models which capture these phenomena. These components are represented at the various scales in Figure 3.1. This paper focuses on VUQ results for the 1.5 MWth oxy-coal furnace (L1500) brick in the burner-scale validation level seen in Figure 3.1. The L1500 has a transversal area of 1.05 x 1.05 m2 and is 14.2 m long. It is divided into 10 sections, as seen in Figure 3.2. A cutaway view through the middle of the burner and the first four sections of the furnace is shown in Figure 3.3. Four additional sets of cooling tubes (not shown) are located adjacent to the opposite interior wall. This experimental facility is discussed in greater detail in [2, 27, 28]. 3.3 Description of VUQ approach Our analysis employs a version of the Simulator Assessment and Validation Engine (SAVE) framework, originally developed by the National Institute of Statistical Sciences [6], and modified by Schroeder [77]; see Figure 3.4. Our focus is on data collected during a week-long experimental campaign in February of 2015. Based on overall programmatic goals related to incident heat flux, the quantities of interest (QOIs) in this analysis are five wall temperature measurements (sections 1, 2, 3, 4, 6), heat removal by eight sets of cooling tubes (one set each on north and south sides of sections 1-4), and oxygen (O2 ) mole fraction in the furnace exhaust. In a previous paper [14], we described steps 1 and 2 of the framework and identified two parameters that had a first order impact on the QOIs: thermal conductivity (k deposit ) and emissivity (ε w ) of the deposit. This paper reviews step 2 and then details steps 3 through 6. 3.4 Construction of input and uncertainty map The input and uncertainty (I/U) map is a list of parameters and their associated uncertainty ranges that are considered in the system VUQ analysis. The list may include numerical, scenario, and/or model parameters. An overall priority is assigned to each parameter that signifies its relative importance on the QOIs. 3.4.1 Parameter selection In our previous paper [14], we carried out a sensitivity analysis with five parameters selected from the I/U map. All parameters came from the char oxidation and the ash deposition models. Our conclusion was that the most sensitive parameters across the 44 different types of measurements were k deposit and ε w and that these parameters should be given the highest priority and made the active parameters in steps 3-6 of the VUQ analysis. However, because of the computational time required to perform a VUQ study with even 2-3 parameters, we performed some scoping tests to see if the sensitivity analysis had identified the appropriate active parameters. We found via qualitative comparisons that the wall temperature and heat removal values from the simulation were consistently higher than the experimental values. To determine the source of this inconsistency, we ran additional scoping simulations. First, we increased the ranges of the char oxidation parameters. We thought that by reducing the rate of char oxidation, a lower temperature in section 1 would result. However, this did not occur. Second, we analyzed the effect of wall thermal conductivity (k w ) on the wall temperature profile. While increasing the thermal conductivity did reduce the section 1 temperature, it also reduced the wall temperatures in the other sections, making those data points inconsistent. Third, we changed k w to be a function of temperature. Using an arctangent equation with three parameters, we produced a range of thermal conductivities in the furnace wall. The result was a lower wall temperature in section 1 without reducing wall temperatures in sections 3, 4, and 6. However, the thermal conductivity in the high temperature range was similar to that of stainless steel, which is unrealistic. Fourth, we compared a movie of our LES simulation with a video of the L1500 furnace operating in 100 % swirl mode. The L1500 video showed a cold zone in section 1 that was not present in the L1500 simulation movie. We hypothesized that our data inconsistency in section 1 resulted from not having a cold zone in the simulation. Based on this analysis, our objective was to reproduce this cold zone in the simulation. The L1500 inlet boundary condition was a handoff plane approach as described by Dı́az-Ibarra et al. [14]. This handoff plane included tangential and axial components of the velocity field exiting the burner. We reduced the tangential velocity components in the handoff plane and verified that we could produce the cold zone. As a result, we added a swirl parameter ( f s ) to our analysis. This parameter (range of 0-1) modifies the tangential velocity components of the inlet velocity field. A factor of 0 means the inlet has no tangential components (0% swirl) and a factor of 1 means it has the maximum tangential velocity possible (100% swirl) from the swirl burner. 45 The final parameter we added to our analysis was ṁcoal . While this parameter was controlled during the experiments, the coal feeding system was hampered by small leaks and did not have a high degree of accuracy. The range of this parameter is the maximum and minimum value in the 30-minute time period described in section 3.5.1. To reduce the size of our VUQ study, we combined the two most sensitive parameters from the sensitivity analysis, k deposit and ε w , into one parameter, which is the effective thermal conductivity of the ash deposit or k e f f . The range of this parameter is determined by a range of 0.1-1.8 for k deposit [72] and a range of 0.3-1.0 for ε w [87, 94]. Table 3.1 is an updated I/U map. 3.5 Experimental data collection and uncertainty analysis The experimental data or QOIs used in this analysis, temperature at the wall, heat removal by cooling tubes, and O2 mole fraction in the furnace exhaust, were collected in the L1500 furnace (see Figure 3.3) [2, 27]. The furnace was oxy-fired with a Utah Sufco coal (bituminous) and a mixture of recycled flue gas and O2 . Two burner swirl conditions were tested, 0 % (no swirl) and 100 % swirl (maximum swirl setting of the burner). Information relating to the coal characterization and to the operating conditions is found in Dı́az-Ibarra et al. [14]. Additional details about the L1500 experimental campaign are found in Fry et al. [28]. 3.5.1 Experimental data collection Figure 3.5 shows a plot of heat removal by the section 1 north and south side cooling tubes on February 27, 2015, during the period from 10 a.m. to 3 p.m. The burner swirl setting was switched from 0% to 100% swirl in this interval. The heat removal increased immediately after the switch but then decreased to pre-switch heat removal levels within two hours due to ash deposition on the tube surfaces. For this analysis, we selected the data interval from 3:00 - 3:30 p.m. (100 % swirl mode) because the furnace was stable and no intrusive measurements were being made. This interval is represented by the red box in Figure 3.5. While this figure only shows the heat removal data for two sets of cooling tubes, the other data sets (heat removal by other sets of cooling tubes, wall temperatures, O2 mole fraction) were extracted for this same time period. 46 As seen in Figure 3.5, the heat removal on the north side of section 1 is higher than that on the south side. We hypothesize that the burner geometry and the positions of the cooling tubes caused this behavior. We base this hypothesis on visual observations of the "flame" through optical access at several locations along the length of the furnace. The "flame" was not symmetric but drifted toward the north side. 3.5.2 Experimental measurement uncertainty There are generally two types of error in experimental data, sampling error (random) and systematic error (bias) [63]; see Eqn. (3.1). Sampling error (β) is due to the finite number of data points used to represent a population. Systematic error (δ) arises from the measurement instrument, which could result from a problem with the instrument itself or with the way the experimentalist uses the instrument. ȳ − y T = δ + β (3.1) We estimated the sampling error with Eqn. (3.2) using data saved every second over the 30-minute interval described in the previous section. In this equation, tα/2,v is a factor computed with the t-distribution and a given confidence interval, s is the standard deviation, and n is the number of data points (1800). s | β| ≤ tα/2,v √ n (3.2) We present the concept of instrument models as a method for analyzing potential sources of systematic error and quantifying the magnitude of those errors. All of the QOIs recorded by the data collection system were actually derived quantities. What was actually measured, for example a voltage signal, was then converted into the desired variable via a calibration curve or other model. The model, or set of models, we used to convert the data we collected (mV) into a value of interest (heat removal in W) is the instrument model. 3.5.2.1 Heat removal 3.5.2.1.1 Systematic error for heat removal. The heat removal instrument model consisted of five equations. In Eqn. (3.3), we computed the heat removal based on the heat capacity of water at 283 K (c p,water = 4192.1 J kgK ), the measured water mass flow rate (ṁwater ), and the temperature increase between the water entering (Tin,water ) and exiting 47 (Tout,water ) each set of cooling tubes as measured by type K thermocouples. Based on literature from the manufacturer [65], we used a value of ±1.1 K as an estimate of systematic error for each thermocouple temperature measurement. HR = c p,water ṁwater ( Tout,water − Tin,water ) (3.3) GEMS RotorFlow sensors [29] measure revolutions in hertz (Hz) to determine ṁwater indirectly. We employed a calibration procedure to convert Hz to kg/h (kilograms of water per hour) for each set of cooling tubes. From the linear calibration curves for the north and south cooling tubes in sections 3 and 4, we computed a single linear calibration curve from the mean slope (m̄) and mean intercept (b̄) of these four calibrations curves. We then used Eqn. (3.4) and Eqn. (3.5) to compute water mass flow with standard deviation for each recorded Hz measurement in the time interval of interest. Here, x is the measurement in Hz, N is the number of calibration curves, and ṁstd,water is the systematic error for the water mass flow rate. ṁmean,water = m̄x + b̄ ṁstd,water = 3.5.2.1.2 ∑iN=1 (mi − m̄) x + (bi − b̄) N (3.4) 1/2 (3.5) Sampling error for heat removal. We computed the sampling error associ- ated with the water mass flow rate and the two temperature measurements in Eqn. (3.3) using Eqn. (3.2) with a confidence interval of 95%. 3.5.2.1.3 Total error for heat removal. We computed the overall heat removal error from estimates of the lower and upper bounds given in Eqn. (3.6) and Eqn. (3.7). In these equations, errorṁ and error T are the total errors obtained by adding the sampling and systematic errors. These error bars are plotted in Figure 3.6 for the heat removal data from the eight sets of cooling tubes. ¯ water + errorṁ ) HRupper = c p,water (ṁ (3.6) ( Tout,water − Tin,water + errorT ) ¯ water − errorṁ ) HRlower = c p,water (ṁ ( Tout,water − Tin,water − errorT ) (3.7) 48 3.5.2.2 Wall temperature Temperature measurements (Ttc ) at the inside wall of the furnace were taken in sections 1, 2, 3, 4, and 6 using Type B thermocouples encased in ceramic sheaths. These sheaths were inserted into small holes in the furnace ceiling located in the middle of each section. Each sheath was inserted until it was flush with the inside wall of the furnace, as shown in Figure 3.7. Insulation material was then stuffed in the gap (interface) between the ceramic sheath and the refractory wall. In the L1500 simulations, we did not resolve the refractory wall as it would have been too computationally expensive. Rather, we computed a steady state solution for heat transfer through the wall, including the inside wall temperatures, assuming the properties of the refractory and insulation materials composing the wall. However, we could not directly compare these simulation wall temperatures to the temperatures measured experimentally because the heat transfer profile through the sheath and the insulation material surrounding the sheath was different than the heat transfer profile through the refractory wall. To compare the simulation and experimental data, we required an instrument model to correct the experimentally-measured wall temperatures. This correction accounted for differences in material properties between the ceramic sheath and the refractory wall, as shown schematically in Figure 3.7. 3.5.2.2.1 Systematic error for wall temperature. The instrument model solves the two-dimensional energy equation for a solid using the physical properties of the wall material, the thermocouple encased in the sheath, and the interface material around the sheath (see Figure 3.8). This model accounts for property variation with temperature in the refractory and insulation layers of the furnace wall; the properties of the thermocouple encased in the ceramic sheath and of the insulation material around the sheath are assumed to be constant over the temperature range of the analysis. We assume that the ash layer being deposited on the furnace wall has no impact on the heat transfer through the wall. This assumption is based on simulation results and visual observations showing that slagging conditions exist in sections 1 through 4. The model solves the set of ordinary differential equations (ODE) presented in Eqn. (3.8) using the LSODA [1, 9] ODE solver. Eqn. (3.8) is the discretized energy equation using the finite volume method. In this equation, Ti,j is the temperature of the solid at the center of cell i, j. The density ρi,j and 49 heat capacity c p(i,j) are computed based on the material in and the temperature of cell i, j. The thermal conductivities at cell faces (k l (i,j) , etc.) are computed with Eqn. (3.9). The thermal conductivity at the cell center, k i,j , is a function of the type of material (k tc , k inter f ace , kre f , k ins ) and the temperature (refractory and insulation layers only). Figure 3.8 is colored by kr(i,j) at 300 K. Finally, ∆y and ∆x are the cell size. ρi,j c p(i,j) Ti,j − Ti,j−1 − Ti,j + Ti,j+1 dTi,j = −k l (i,j) + kr(i,j) 2 dt ∆y ∆y2 Ti,j − Ti−1,j − Ti,j − Ti,j+1 − k b(i,j) + k t(i,j) ∆x2 ∆x2 k b(i,j) = 0.5(k i,j + k i−1,j ) (3.8) (3.9) k t(i,j) = 0.5(k i,j + k i+1,j ) k l (i,j) = 0.5(k i,j + k i,j−1 ) kr(i,j) = 0.5(k i,j + k i,j+1 ) The boundary conditions for Eqn. (3.8) were as follows. At the top (see Figure 3.8), we used a shell temperature of 290 K. We tested the model with shell temperatures ranging from 300 K to 500 K, and the results were insensitive to this parameter. At the bottom, we computed the net heat flux at the interface from Eqn. (3.10) where Tw is the temperature at the inside wall and ε is the emissivity. At the thermocouple location, we used ε = ε tc . At all other locations along the inside wall, we used ε = 1. We obtained the range of qincident values shown in Table 3.2 from a simulation of the L1500. We picked a wide range to cover all the possible qincident values along the length of the furnace. On the left and right sides of the computational domain, we applied a Neumann boundary condition. qnet = ε(qincident − σTw4 ) (3.10) The initial condition for the whole wall was 300 K. We ran the simulation for 60 h of physical time in order to reach steady state. The main goal of this instrument model was to correct the experimentally-measured temperature (Ttc ) to a surface refractory wall temperature (Tw ). We determined appropriate ranges for model parameters by performing a consistency analysis (described in section 3.8) with each Ttc (sections 1, 2, 3, 4, and 6) and the output of this instrument model, Ti,j , at the location labeled "Thermocouple" in Figure 3.8. 50 In Table 3.2, we present the I/U map for the thermocouple instrument model. We explored the four parameters with the highest priority of 6 ( k tc , k inter f ace , ε themo , and qincident ) and fixed the parameters with the lowest priority of 1 (diameter of the hole where the thermocouple sheath was placed, φhole , and the resolution of the computational domain, ∆x and ∆y). For k tc , the maximum value assumes the ceramic sheath is made from "OMEGATITE 450" [64], which has a thermal conductivity of 8.65 W mK W mK at 1073 K. The minimum value of 2 is the thermal conductivity of the refractory material (Greencast-pu Plus) at 1588 K. The material between the furnace wall and the thermocouple sheath is either air or insulation (see Figure 3.7). Since k inter f ace represents this material, we used a thermal conductivity range between that of air at 293 K and of the insulation. We assumed the insulation material had similar thermal conductivity to the insulation layers in the refractory wall. We did not have information regarding the range of ε tc , so we used a wide range of 0.1 - 1. The procedure to correct Ttc is presented in Figure 3.9. We identified a sample of 1296 cases within the ranges of the four active parameters from Table 3.2 using the Uncertainty Quantification Toolkit (UQTk) [13]. We then ran the instrument model for these 1296 cases, extracted a value Ti,j from the location labeled "Thermocouple" in Figure 3.8 (step (1) in Figure 3.9), and created a quadratic response surface (surrogate model) from these 1296 outputs (step (2) in Figure 3.9) [19]. Next, we computed the mean values for the thermocouple measurements using the experimental data from the selected time interval. For each measurement, we estimated the upper and lower bounds from the sampling error; see Eqn. (3.2). With the surrogate model and the upper and lower bounds of the thermocouple measurements, we performed a consistency analysis following the procedure described in section 3.8 for each thermocouple measurement. From the consistency analysis, we obtained a consistent sample space of the four parameters under study for each measurement (step (3) in Figure 3.9). We then extracted a value Ti,j from the location labeled "Wall temperature" in Figure 3.8 from the 1296 runs and created 1296 quadratic surrogate models (step (4) in Figure 3.9). Finally, we evaluated the Tw surrogate model within the consistent sample space obtained in step (3) to obtain the estimated uncertainty due to the systematic error (step (5) in Figure 3.9). The upper and lower bounds for Tw are the maximum and minimum values for Tw computed at each thermocouple location. 51 3.5.2.2.2 Sampling error for wall temperature. We computed the sampling error associated with Ttc using Eqn. (3.2) with a confidence interval of 95%. We obtained a value of ±2.3 K. Because this error is small relative to the systematic error, we ignored its contribution in the total error analysis. 3.5.2.2.3 Total error for wall temperature. The results of this analysis are presented in Figure 3.10. The systematic error due to the difference in physical properties between the wall and the thermocouple sheath is approximately ± 115 K. This error range is a prior estimation for the uncertainty of the wall temperature. In section 3.8, we perform a consistency analysis that reduces this uncertainty range. 3.5.2.3 Oxygen molar fraction. 3.5.2.3.1 Sampling error for O2 molar fraction. We computed the sampling error associated with O2 molar fraction using Eqn. (3.2) with a confidence interval of 90%. 3.5.2.3.2 Systematic error for O2 molar fraction. We assumed a systematic error of 0 for O2 molar fraction because we have not yet developed an instrument model to compute this error. 3.6 Simulation data collection and model development We collected simulation data for the L1500 using Arches, a component of the Uintah software suite [48]. Arches uses a large eddy simulation (LES) approach to resolve the turbulent flow field of the L1500. The solid (coal) phase is modeled using the direct quadrature method of moments with three environments [25]. We simulated the fluid flow through the complex geometry of the burner separately using the commercial software package STAR-CCM+. We then used handoff files to transfer data from the output of the burner simulation to the input for the Arches simulation; see Diaz-Ibarra [14] for details about this procedure. We previously determined that a computational domain encompassing the first seven sections of the L1500 delivered sufficient accuracy for this analysis [14]. Our mesh resolution for these simulations was 15 mm, resulting in a computational domain with 2,255,610 cells. Because the Uintah framework provides efficient parallelization, we performed the simulations required for this analysis on two high-performance computing platforms, one 52 at the University of Utah and the other at Lawrence Livermore National Laboratory. We ran a typical simulation to 30 s: 25 s to reach steady state and 5 s for data averaging (see Figure 3.12). Each simulation required approximately 750 processors for 96 hours to complete; see Dı́az-Ibarra [14] for additional details about the L1500 furnace simulations. 3.6.1 Wall model A new wall model has been tested and implemented in Arches for the refractory walls of the L1500. This model allows for both negative and positive net heat flux. Negative net heat flux occurs when a cold eddy passes by the wall surface after a hot eddy has passed by and heated up the wall. The main goal of this model is to compute the fire-side surface wall temperatures at steady state. The walls consist of layers of refractory and of insulation materials (see Figure 3.7). In this model, the relaxation coefficient is given by the properties of the wall (density and heat capacity). The wall model solves two energy equations, one for the refractory material and the other for the insulation material. The energy balance for the refractory material is shown in Eqn. (3.11). In this equation, Fsur f ace is the net heat flux from the combustion environment as computed by Eqn. (3.12), F1,2 is the net heat flux exchange between the refractory layer (1) and the insulation layer (2) as computed by Eqn. (3.13), and L1 = 21.6 cm is the thickness of the refractory layer. The thermal inertia term ρ1 c p,1 = 2.0e6 J ) m3 K is computed from density and heat capacity data obtained from the manufacturer [44]; In Eqn. (3.12), qincident is the incident heat flux at the wall computed in the simulation, T1 is the temperature of the refractory layer, and ε is emissivity. In Eqn. (3.13), T2 is the temperature of the insulation layer, k1 = 2.36 W mK is the thermal conductivity of the refractory layer, and L2 = 12.7 cm is the thickness of the insulation layer. Fsur f ace − F1,2 dT1 = dt L1 ρ1 c p,1 (3.11) Fsur f ace = ε(qincident − σT14 ) (3.12) F1,2 = k1 ( T1 − T2 ) 0.5( L1 + L2 ) (3.13) The energy balance for the insulation layer is shown in Eqn. (3.14). Here, F1,2 is computed with Eqn. (3.13), F2,shell is the net heat flux exchange between the insulation layer and the external wall as computed from Eqn. (3.15), and ρ2 c p,2 = 2.9e5 m3J K is the product of 53 the density and the heat capacity of the insulation layers based on data from the manufacturer [44]. In Eqn. (3.15), k2 = 0.15 W mk is the thermal conductivity of the insulation layers and Tshell is the shell temperature measured during the experimental campaign as shown in Table 3.3). F1,2 − F2,shell dT2 = dt L2 ρ2 c p,2 F2,shell = k2 ( T2 − Tshell ) 0.5L2 (3.14) (3.15) This model is an energy balance on the wall with only two control volumes. Thus, our main concern is whether the coarseness of the model affects the precision of the surface temperature. To test the sensitivity of the model to the number of cells, we developed a version of this model with an arbitrary number of cells. We present results from two-cell and 100-cell cases in Figure 3.11; both models were run for 12 h of physical time. We changed the time mean (µq ) and standard deviation (σq ) of qincident three times to emulate the different firing regimes during a typical day of testing (natural gas fired overnight, coal fired with 0 % and 100 % swirl). In the first 4 h, µq = 9.07e5 W/m2 and σq = 9.07e4 W/m2 . In the next 4 h, µq = 2.87e5 W/m2 and σq = 2.87e4 W/m2 . For the last 4 h, µq = 3.72e5 W/m2 and σq = 3.72e4 W/m2 . For each firing regime, we assumed a normal distribution of qincident and sampled from that distribution to obtain qincident at each time step. In Figure 3.11, the blue line (100 grid points) corresponds to the surface wall temperature, and the green line (2 grid points) corresponds to T1 (refractory layer temperature). There are some slight differences in model behavior before the cases reach steady state. However, once steady state is reached, the differences in surface wall temperature are negligible. Consequently, we implemented the simpler two-grid-point model, Eqn. (3.11) and Eqn. (3.14), in Arches. We solved the two-grid-point model using a first order scheme in time. However, Eqn. (3.11) and Eqn. (3.14) cannot be solved using the same time step as the LES simulation because it would take too long to reach steady state. Instead, we used a time step that was 5000 times larger than the time step of the LES simulation. Therefore, 1 s of the LES simulation corresponded to 5000 s in the wall model. To check whether this new model reached steady state, we extracted wall temperatures from simulation output at the five locations corresponding to the thermocouple measurement locations and plotted them in 54 Figure 3.12. After 15 s of LES simulation time, all the thermocouple locations have reached steady state. The cooling tubes used the Arches wall model described in Diaz-Ibarra [14]. 3.6.2 Suite of simulations for consistency analysis We ran a suite of 34 cases in Arches in which we varied the three parameters in Table 3.1. This suite of cases is graphically displayed in Figure 3.13 where each red point represents a simulation that was run for a least 30 s of simulation time. 3.6.3 Data extraction and postprocessing Using post-processing tools described in Diaz-Ibarra [14], we extracted data relating to the QOIs from the 34 simulations and averaged the data over the last 5 s of simulation time. Since heat removal (Qremoval ) by the cooling tubes is a derived quantity, we computed it from Eqn. (3.16) during postprocessing. In this equation, A corresponds to the surface area of the cell faces that are in contact with the hot gas, Tw is the surface temperature of the cooling tube, qincident is the incident heat flux on the tube, and ε is the surface emissivity of the tube. The averaged values of Qremoval in all the cells comprising a set of cooling tubes were summed up to obtain the total heat removal by that set of tubes. This procedure was repeated for all eight sets of cooling tubes. N Qremoval = ∑ ε(qincident − σTw4 ) A (3.16) n =0 To obtain the wall temperature simulation data, we extracted the averaged wall temperatures at the five thermocouple locations. For O2 molar fraction, we performed both time (last 5 s of simulation time) and spatial (exit plane of the simulation) averages of the mass fractions of O2 , CO2 , H2 O, SO2 , and N2 and then computed the O2 molar fraction on a dry basis (to match the experimental data) using this gas composition. 3.7 Construction of surrogate models The consistency analysis performed in section 3.8 requires thousands of function evaluations. Since the evaluation time of our function, an LES simulation, is on the order of days, we created surrogate models of the simulation outputs using Gaussian process 55 (GP) response surfaces [70, 76, 77]. Each GP surrogate model was created with the Arches simulation data obtained as described in section 3.6.3. The number of GP surrogate models was equal to the number of experimental measurements. We had 14 total measurements and thus 14 GP surrogate models: eight heat removal measurements, five corrected wall temperatures, and one O2 mole fraction measurement. In Figure 3.14, we present verification plots for three GP surrogate models: wall temperature in section 1 (T1 ), heat removal by the set of cooling tubes on the north side in section 1 (HR2 ), and O2 mole fraction in the exhaust. In these plots, the x axis is the surrogate output for the variable and the y axis is the Arches simulation output for the same variable. The green line is y = x. If the surrogate model is a good representation of the simulation data, all the x, y pairs should be close to y = x, which they are. Thus, the GP surrogate models do reproduce the output values from the Arches simulation. In Figure 3.15, we performed additional verification by evaluating the same three GP surrogate models with samples of 10,000 runs. In the first column, we used ṁcoal = 290 lb/hr and varied k e f f and f s over the ranges presented in Table 3.1. In the second column, we used f s = 0 and varied k e f f and ṁcoal over the ranges presented in Table 3.1. The GP surrogate model plots of T1 correspond to the first row in Figure 3.15. In the first column, the temperature increases with increasing f s until f s = 0.4. At this point, the temperature is approximately constant. In the second column, the temperature decreases with increasing k e f f , which means more energy is being removed from the reactor. The GP surrogate model plots of HR2 correspond to the second row in Figure 3.15. In the first column, we see that the heat removal is lower in the region f s ≤ 0.4 than in the region f s ≥ 0.4 , which corresponds to the lower temperatures seen in the T1 GP surrogate models. In the second column, increasing k e f f increases the heat removal. The plots of O2 molar fraction are shown in the third row of Figure 3.15. The scatter plots in both columns show that if f s increases, the O2 molar fraction also increases, but only by a small amount. On the other hand, if ṁcoal increases, the O2 molar fraction decreases because more coal means more O2 is consumed during combustion. We performed a similar verification analysis for the other GP surrogate models. Based on these qualitative analyses, we determined that all the surrogate models demonstrated adequate behavior. 56 3.8 Analysis of model outputs The VUQ methodology employed for this study is a consistency measure analysis referred as bound-to-bound consistency [75]. The basic concept of this consistency analysis is the comparison of model outputs with experimental data using Eqn. (3.17). In this equation, ym,e ( x ) is the simulation data point obtained from a simulation defined by the set of x parameters and ye is the experimental data point. If ∆ ≤ 1, the simulation data point using parameter set x is consistent with the experimental data point. If the simulation outputs for a parameter set x are consistent with all the experimental measurements, x is a consistent point in our analysis. ∆= |ym,e ( x ) − ye | ue (3.17) To perform the consistency analysis, we first computed ue and ye . Here, ue is the difference between the upper and lower bounds divided by two for each of the 14 measurements (ue = ye,upper −ye,lower ). 2 We estimated these bounds in section 3.5.2.2 for the thermocouple, section 3.5.2.1 for the heat removal, and section 3.5.2.3 for the O2 mole fraction measurements. We computed ye , the mean experimental value, by averaging the upper and lower bounds of each measurement (ye = ye,upper +ye,lower ). 2 Next, we created a random sample of 500,000 points within the parameter space defined by the ranges of the three parameters under study ( f s , ṁcoal , k e f f ) as presented in Table 3.1 and evaluated each one of the 14 GP surrogate models at all 500,000 points. We then applied the consistency analysis to the 14 measurements in our data set. For each point (set of parameter values) in the 500,000 point set, we extracted the maximum ∆, ∆max , of the 14 ∆ values that we computed. From this vector of ∆max values, we extracted the minimum ∆max , which is the most consistent point from the 500,000 points. If this minimum ∆max < 1, then we have found at least one consistent point. We reviewed the ∆max vector and selected all points for which ∆max < 1. From this group of points, we obtained a new range for the three parameters by finding the global maximum and minimum values of all the consistent points. This new range is presented in Table 3.4. If the minimum ∆max > 1, we need to improve our models, explore other parameters, or review ue . Our last step in the consistency analysis was to perform a refined consistency analysis 57 with the new (consistent) range. We created a random sample of 10,000 points within the consistent range and then reevaluated the 14 GP surrogate models at these 10,000 points. We next obtained a set of consistent points where ∆max < 1; these points are presented in Figure 3.16. These points are simultaneously consistent with all 14 QOIs. In this figure, the ranges for the three axes are the original ranges presented in Table 3.1 and the box containing the points is the consistent range. The color of each point corresponds to ∆max in the consistent samples. The consistent ranges for the three parameters are presented in Table 3.4. These ranges are the limits of the box that contains the consistent points (see Figure 3.16). All ranges were reduced with the help of the experimental data. Thus, the consistency analysis improves the simulation tool by refining (e.g. calibrating) the parameters that have a primary effect on the QOIs. It also improves the quality of the experimental data by reducing the uncertainty bars. We compare the experimental (red error bar), simulation (blue error bar), and consistency ranges (green error bar) for the wall temperature, heat removal, and O2 mole fraction in Figure 3.17. In all three sets of data, the consistency range is smaller than the experimental uncertainty range. The uncertainty in Tw was reduced from approximately ± 115 K to ± 16 - 26 K. The uncertainty in heat removal measurements prior to the consistency analysis was approximately ± 1.05e4 W to ± 1.24e4 W. After the consistency analysis, the uncertainty was reduced to approximately ± 0.08e3 W to ± 0.58e3 W. The consistency analysis reduced the uncertainty in all the experimental measurements by forcing simultaneous consistency of all 14 data points with simulation output. 3.9 Feedback and feed-forward In this step, we review our analysis, check assumptions, and make recommendations for the next VUQ cycle. We explored three parameters (k e f f , f s , and ṁcoal ) using 14 measurements (eight heat removal values, five wall temperatures, and one O2 mole fraction). 3.9.1 Parameters The parameter k e f f has a direct impact on the heat removal by the cooling tubes and on the wall temperatures as noted in the surrogate model scatter plots (see Figure 3.15). 58 This parameter represents the thermal conductivity and the surface emissivity of the ash layer on the cooling tubes as well as a correction for the changing external surface area of the cooling tubes. In the simulation, we use a nominal area to compute the heat removal. However, in the reactor, this surface area may be bigger because of the ash deposits on the surface. The consistency analysis produces a calibrated range for k e f f of 0.97-1.51 (see Table 3.4). This range indicates that using a constant value for this parameter in simulations is not adequate. Thus, we recommend developing a model that accounts for the variation in this effective thermal conductivity. This analysis shows the L1500 reactor can be used to develop models related to ash deposition that have large impacts on measurements related to heat flux such as heat removal and wall temperature. The swirl parameter, f s , also has a direct impact on the heat removal and wall temperature measurements. In Figure 3.15, the wall temperature in section 1 increases until it reaches a maximum value at f s = 0.4. To better understand this behavior, we present instantaneous plots of gas temperature for three Arches simulations in Figure 3.18; this figure represents a slice through the middle of the computational domain in the y plane. The vertical black lines show the position of the wall thermocouple in sections 1, 2, and 3. All three cases have the same ṁcoal and k e f f . For the case with f s = 1, the tangential velocity range is 0 - 75 m/s. When f s = 0.1, the tangential velocity range is 0 - 7.5 m/s. When f s = 0, the tangential velocity is zero. From the figures we see that when f s = 1, the gas temperature near the wall thermocouple in section 1, T1 , is higher than the same temperature for the other two cases. Therefore, if f s = 1, we have a short, wide zone of higher temperatures (red color in the figure) in section 1. If f s = 0, the long, narrow zone of higher temperatures extends through sections 1, 2, and 3. If f s has a value between 0 and 1, the behavior is between this two extremes as we can see when f s = 0.1. The experimental wall temperature data show a maximum temperature in section 2 (see Figure 3.10). Since f s = 1 produces a maximum temperature close to the burner in section 1, we need to vary f s if we want to obtain consistency. After performing the consistency analysis, we found that the range for f s was 0. - 0.1, which is close to the 0% swirl case. The low range for f s indicates a potential problem with the burner inlet condition since we were operating in 100 % swirl mode ( f s = 1 ). We performed a high-resolution LES 59 simulation of the burner given the fabrication drawings for the burner and the operating flow rates for the various inlet streams (oxidant and coal) to obtain velocity profiles at the burner exit plane [14]. Thus, we need to review the assumptions that were made in the simulation such as (1) no leaks in the burner, (2) the geometry in the fabrication drawings accurately represents the actual burner, and (3) the mechanical joints such as screws do not affect the fluid dynamics. From this analysis, we know what the tangential velocity has to be reduced in order to obtain consistency. The coal feed rate, ṁcoal , has a strong effect on the O2 molar fraction (see Figure 3.15, third column). It also has an effect on Tw and heat removal because increasing or decreasing ṁcoal increases or decreases the amount of energy in the reactor. The consistent range, 275 - 286 lb/h (see Table 3.4), is small and biased toward the lower feed rate. This behavior could be due to a coal leak during the experimental campaign. 3.9.2 Measurements The measured temperatures, Ttc , were corrected to Tw with the instrument model presented in section 3.5.2.2. While the variation in temperature measured by these thermocouples was small (standard deviation ± 2 K) over the 30-minute time period under analysis, we found that the uncertainty in Tw was ∼ ± 115 K. This case is a good example of the importance of estimating bias (systematic) errors in experimental measurements. The main cause of this systematic error was the difference in thermal conductivity among the thermocouple, the insulation material that was used in the gap between the thermocouple and the wall, and the refractory and insulation materials that made up the wall. To improve this measurement, we recommend that the thermocouple be imbedded in the same material of which the wall furnace is composed. The large uncertainty in heat removal (see in Figure 3.6) arises from several sources. The systematic error, corresponding to 78% - 85% of the total error, is due to calibration and to the small difference between the inlet and outlet water temperature (∼ 20 K). The error in two type K thermocouple measurements is ∼ 2.2 K, which is ∼ 10 % of the water temperature difference. The sampling errors are also substantial, representing 15%-22% of the total uncertainty. These errors were caused by ash deposits that were not controlled during the experimental campaign. 60 We assumed that the experimental data were taken when the L1500 was at steady state. This assumption is justified for the thermocouple-temperature data because the rate of temperature change was small; for the 30-minute time period analyzed in this study, the standard deviation was ∼2.3 K. However, the heat removal data shows transient behavior (see Figure 3.5), probably due to ash deposition. We tried to account for transient effects in the experimental error estimation. To reduce the effect of ash deposits on heat removal rates in future campaigns, we will replace four of the eight sets of cooling tubes with flat, water-cooled panels with an accompanying soot blowing system. With this system, we can correlate the effects of ash deposition with the time in the soot blowing cycle and thus obtain better estimates of the error due to the transient effects. 3.9.3 Models In the Arches simulations, we made the assumption that the temperature profile in the furnace wall could be computed with a steady state wall model. In order to review this assumption, we performed a simulation with the 2D wall model presented in section 3.5.2.2. We ran the model for 61 h with qincident = 5.20e5 W/m2 and then switched to qincident = 2.87e5 W/m2 for 6.3 h to emulate the switch from natural gas to coal in this experimental campaign. In Figure 3.19, we present three temperature profiles from this simulation: fire-side surface of the wall (0 cm), 1.27 cm (0.5 in) in the wall, and 10.16 cm (4 in) in the wall. The temperature at the surface reaches steady state in less than 2 h, while the temperature 1.27 cm (0.5 in) in the wall takes ∼ 5 h, and the temperature 10.16 cm (4 in) in the wall takes 6.3 h. These long transients are due to the thermal inertia of the refractory walls, which makes the heating process slow. Since the experimental wall temperature measurements were near the wall surface, we conclude that our assumption of a steady state wall condition is reasonable. However, for temperature measurements deeper in the refractory wall (e.g. farther from the hot face), the associated instrument model must include transient effects. 3.10 Conclusions We presented results from our VUQ analysis of an oxy-coal data set collected in the L1500 furnace. We employed a six-step methodology that was based on bound-to-bound 61 consistency analysis. Through this analysis, we reduced the uncertainty in the experimental data and in three simulation parameters: k e f f , f s , and ṁcoal . The consistency ranges for these parameters were: k e f f = 0.97-1.51 w mK , f s = 0.0-0.1, and ṁcoal = 275-283 lb h. The k e f f parameter, which includes both ε w and k deposit , had a large impact on wall temperatures in and heat removal from the L1500. Thus, we conclude that this reactor can be used to perform studies on ash deposition and to develop and validate models related to heat transfer to cool surfaces. The f s parameter, which affects the inlet condition, also had a large impact on the L1500 measurements, making the L1500 an adequate reactor for studying improvements to burner design and for testing swirl effects. The ṁcoal parameter had a small effect on wall temperature and heat removal measurements; its impact on the O2 mole fraction in the exhaust was much larger. This methodology requires uncertainty bounds on the experimental data that include both the sampling and systematic errors. We presented a procedure to estimate the systematic error in a thermocouple device through an instrument model and the application of a VUQ methodology. The instrument was the 2D representation of the energy equation for the wall of the L1500. With this procedure, we estimated an error of approximately ±115 K in the thermocouple measurements of wall temperature, which is much greater than the estimated random error of ∼ ±2.3 K. We then performed the consistency analysis and were able to reduce the experimental error from ∼ ± 115 K to ∼ ± 16-26 K. We also estimated the systematic error in heat removal data through cooling tubes with an instrument model. We found that the error due to calibration was ∼ 80 % of the total error. We then performed the consistency analysis and were able to reduce the experimental error from ∼ ± 1.06e4 W to ∼ ± 0.58e4 W. From the wall instrument model, we found that the surface wall temperature reaches steady state in a few hours for the conditions tested in the L1500. Therefore, performing an LES simulation with a steady state wall model is a reasonable approach. However, measurements of temperatures deeper in the wall need an instrument model in order to account for transient effects. 62 3.11 Acknowledgment This material is based upon work supported by the Department of Energy, National Nuclear Security Administration, under Award Number(s) DE-NA0002375. We would like to acknowledge the Arches development team, especially Derek Harris and Minmin Zhou for their contributions to the Arches code. We also wish to thank the team involved in the L1500 experimental campaign, including Ignacio Preciado, Teri Draper Snow, Eric Eddings, and David Wagner. Table 3.1. Input/uncertainty map for consistency analysis Parameter Range min max Scenario parameters f s [−] 0. 1. lb ṁcoal [ h ] 275. 310. Model parameters (ash deposition) k e f f [Wm−1 K −1 ] 0.1 5. Table 3.2. Input/uncertainty map for the thermocouple instrument model Parameter Priority Range min max W k tc [ mK ] 6 2.0 8.65 k inter f ace 6 0.0257 0.15 W [ mK ] ε tc [−] 6 0.1 1. qincident [ mW2 ] 6 1.88e5 5.95e5 Nominal Value φhole [m] 1 0.05 ∆x [m] 1 0.00635 ∆y [m] 1 0.01 63 Table 3.3. Shell temperatures averaged over all measurements made on a side Location Shell Temperature [K ] Quarl 434 362 Main chamber south side Main camber north side 396 Main chamber bottom side 362 Main chamber top side 427 Table 3.4. Parameter ranges after consistency analysis Parameter Prior range Consistent range w k e f f mK 0.1-5. 0.97-1.51 fs 0.0 1. 0.0- 0.1 ṁcoal lbh 275. - 310. 275. - 283. Figure 3.1. Validation hierarchy for CCMSC. Nominal value 1. 290. 64 ction 10 se s To convection zone South side Sample ports North side Burner Figure 3.2. Schematic of the L1500 multifuel combustion furnace. Cooling tubes Burner Sections 1 2 3 4 Figure 3.3. Schematic of the first four sections of the L1500 furnace. Selection of quantities of interest (QOIs) Step 1 Construction of input and uncertainty (I/U) map Step 2 Collection of experimental and simulation data Step 3 Construction of surrogate models Step 4 Analysis of model outputs - consistency analysis Step 5 Feedback and feedforward Step 6 Figure 3.4. Six-step methodology with consistency analysis 65 0% swirl setting 100% swirl setting Selected data Figure 3.5. Heat removal data in section 1 on February 27, 2015. The heat removal decreases continually in both 0% and 100% swirl modes due to ash deposition on the tube surfaces. Figure 3.6. Heat removal from cooling tubes with total error bars; see Eqn. (3.1). 66 Sheath Insulation layers Insulation Material Refractory layer Section 4 Location of thermocouple measurement Gas Figure 3.7. Thermocouple placement in furnace wall. Section 4 is shown but the placement is similar for all thermocouple measurements. Shell Temperature Sheath qincident Insulation layers Refractory layer Neumann Neumann Insulation material Wall temperature Thermocouple Figure 3.8. Representation of the instrument model to correct thermocouple measurement in ceramic sheath to equivalent inside wall temperature. 67 5 0.9 x 10 m2 hW i 0.8 qincident q incident (2) Surrogate model from Ttc 0.7 5 0.6 4 0.5 3 0.4 0.3 2 0.15 0.2 8 0.1 h W i kinterf ace k mK 4 0.05 2 y [m] T emperature [K] interface 0.1 6 h W i ktc kthermo mK (3) Consistency analysis 1580 Ttc 1560 Tw Ttc Tw 1540 Temperature [K] x [m] (1) 2D Model 1520 1500 1480 1460 (4) Surrogate model from Tw 1440 1420 1.8 1.85 1.9 No [−] 1.95 2 2.05 (5) Estimated uncertainty Figure 3.9. Procedure used to correct Ttc to Tw . 1700 Tw 1650 Ttc Temperature [K] 1600 1550 1500 1450 1400 1350 1300 0 1 2 3 4 Section [−] 5 6 7 Figure 3.10. Thermocouple temperature and corrected wall temperature. 68 Figure 3.11. Wall model verification. Selected data Steady state Figure 3.12. Wall temperature at thermocouple locations; section 1 is closest to the burner. Figure 3.13. Suite of simulations that were run in Arches. 69 Figure 3.14. Verification of surrogate models for T1 , HR2 , and O2 molar fraction 70 Figure 3.15. Surrogate models for T1 , HR2 and O2 molar fraction. Figure 3.16. Consistent sample space for all 14 QOIs in L1500 simulations. 71 Figure 3.17. Consistency range for heat removal (top), O2 molar fraction (middle), and wall temperature (bottom). Data are plotted by QOI number. The eight sets of cooling tubes are QOIs 1-8, the O2 mole fraction in the exhaust is QOI 9, and the wall temperatures in sections 1, 2, 3, 4, and 6 are QOIs 10-14. Consistency was obtained simultaneously across the 14 QOIs. 72 T1 T2 T3 fs = 0 fs = 0.1 fs = 1 Figure 3.18. Effect of f s , ṁcoal = 275 [lb/h], k e f f = 3 [W/m/k ], Blue color is 300 K, red is 2100 K qincident = 5.20e5 hW i m2 qincident = 2.87e5 hW i m2 Figure 3.19. Transient temperature effects in the refractory wall CHAPTER 4 REACTING PARTICLE AND BOUNDARY LAYER MODEL: FORMULATION In preparation for submission to Combustion and Flame, Reacting Particle and Boundary Layer (RPBL) model: Formulation. Oscar H. Dı́az-Ibarra, Jennifer Spinti, Philip J. Smith, Christopher Shaddix, Ethan Hecht. 74 The Reacting Particle and Boundary Layer (RPBL) model computes the transient-state conditions for a spherical, reacting, porous char particle and its reacting boundary layer. RPBL computes the transport of gaseous species with a Maxwell-Stefan multicomponent approach. Mass transfer diffusion coefficients are corrected to account for a nonstagnant bulk flow condition using a factor based on the Sherwood number. The homogeneous gas phase reactions are estimated with a syngas mechanism, and the heterogeneous reactions are calculated with a six-step reaction mechanism. Both reaction mechanisms are implemented in Cantera. RPBL computes carbon density (burnout) and uses the Bhatia and Perlmutter model to estimate the evolution of the specific surface area. RPBL solves two energy equations, one for the gas temperature and the other for the particle temperature. The physical properties of the particle are computed from the fractions of ash and carbon in the particle as well as the void fraction. The void fraction is computed assuming a constant diameter particle during the reaction process. RPBL solves a particle momentum equation in order to estimate the position of the particle in a specific reactor. In this paper, we present the model formulation and model outputs for a char particle in the optical entrained flow reactor at Sandia National Laboratories [34, 61]. We explore the sensitivity of two RPBL outputs, particle temperature and velocity, to nine model parameters. The two outputs are found to be sensitive to five of the nine parameters. 4.1 Nomenclature cp External surface area [m2 ] Area of the cell face (i − 1/2 or i + 1/2) [m2 ] Pre-exponential factor in units of [mol, cm2 , s] Molar concentration of the oxidizer in the gas phase at the surface of the particle [ mol ] m3 J Heat capacity [ kgK ] Do D dp E F F1 hk,i hsolid gas jk kr Multicomponent diffusion coefficient (mass basis) [ ms ] 2 Binary Maxwell-Stefan multicomponent diffusion coefficients [ ms ] Diameter of the particle [m] kJ Activation energy [ mol ] Heat flux [ mJ2 ] Drag coefficient factor [−] Gas enthalpy [ kgJ ] Internal particle heat transfer coefficient; hsolid gas = 1 [ W K] kg Mass diffusion flux relative to the mass average velocity [ m2 s ] Reaction constant for heterogeneous reaction [ ms ] Ap A A Cop 2 75 Kg kn m n nt nk Np N l P Pr Re R r rk Si s˙k Sgc Sh Sc t tr T V v vs Wc Wk ẇk,i xc xk Y ξp ρ σr φ τ µg ψ λ ε Index for last gas species Convection heat transfer coefficient [ mW2 ] Mass of ash [Kg] Apparent order of the reaction [−] kg Total mass flux relative to stationary coordinates [ m2 s ] kg Mass flux of species k relative to stationary coordinates [ m2 s ] Number of cells inside of particle Number of cells in whole domain Length scale in reactor [m] Pressure [ Pa] Prandtl number [−] Reynolds number [−] J Gas constant, 8.3144621 K mol Radius [m] kg Char reaction rate [ s ] Heat exchange source term [ mW3 ] mol Heterogeneous molar production rate of kth species [ m 2s ] m2 Specific surface area [ kg ] Sherwood number [−] Schmidt number [−] Time [s] Relaxation time [s] Temperature [K ] Volume [m3 ] velocity [ ms ] Stoichiometric coefficient [−] kg Molecular weight of carbon [ mol ] kg Molecular weight of specie k [ mol ] Molar production rate of species k [ kmol ] m3 s Conversion of carbon [−] Mole fraction of species k [−] Mass fraction [−] Particle area factor to account for internal surface burning [−] Density [kg/m3 ] 2 Specific surface area for heterogeneous reactions, σr = Sgc ρbulk,c [ m ] m3 Void fraction [−] Tortuosity [−] Dynamic viscosity [Pa s] Structure parameter [−] Thermal conductivity [ mWk ] Emissivity [−] Subscripts in f : g: p: Limit of the particle domain Gas Particle 76 ash: char: cond: h: inter f ace: c or cb: k: t: bulk or b: true: solid: initial: i: s: Ash in particle Char in particle Conduction Enthalpy Interface of particle Carbon Species k Total Bulk gas True density Solid particle Initial condition Computational cell Surface 4.2 Introduction Char oxidation is a slow process that takes place in the entire domain of a coal-fired boiler. An accurate computation of char oxidation in the entire computational domain is required to predict the O2 and CO concentrations, particle temperatures (needed for radiative heat flux calculations), and the carbon content of the fly ash. Because char oxidation is a complex process, with heterogeneous and homogeneous reactions, species transport, and the evolution of the physical properties of the porous matrix, a high-fidelity model of the process will have many parameters, some of which have high uncertainty. Experimental studies and modeling approaches for the oxidation of char particles have been research topics for the last seven decades [18, 23, 35, 36, 39, 40, 49, 55, 79-81, 83, 89-91, 93, 95]. The main purpose of all this experimental and modeling work is to compute the kinetic rates of char oxidation. One approach is the global model given in equation (4.1) where the char oxidation rate depends on the oxygen stoichiometric coefficient (vs ), the molecular weight of carbon (Wc ), a kinetic constant for the heterogeneous reaction rates (kr ), a particle area factor to account for internal surface burning (ξ p ), the external area of the particle (A p ), the gas concentration of the oxidizer (Cop ) at the surface of the particle, and the reaction order (n). Murphy and Shaddix [61] presented an experimental and modeling study where a kinetic constant was found for both a single-n th-order Arrhenius expression and for an nth-order Langmuir-Hinshelwood kinetic equation. n rk = vs Wc kr ξ p A p Cop (4.1) 77 While this kind of approach is easy to implement, the result may only apply to the specific conditions of the experiment. A more detailed approach is one where the intrinsic char reaction rate is computed. This approach is based on the specific surface area of the porous matrix; therefore, the internal diffusion and reaction of species are considered. Additionally, the porosity, specific surface area and bulk density need to be computed. There are two different methods for computing the pore diffusion. In the macroscopic method, an effective diffusion coefficient is defined in terms of the molecular diffusion and the Knudsen diffusion coefficient. This effective diffusion coefficient is multiplied by the ratio of porosity to tortuosity. The porosity is then obtained from the density assuming a shrinking core model [83]. In the microscopic model, the diffusion through a single pore is solved and then a statistical approach is used to predict the porosity in the particle [7, 83]. Both mass transport of species and heterogeneous chemical reactions are important phenomena that must be considered in char modeling. Researchers have proposed three zones to describe the behavior of the char reaction rate with particle temperature and size [83]. In Zone I, which occurs at low temperatures and for small particles, chemical reaction is the controlling step. The oxidizer is transported throughout the particle; therefore, the chemical reaction rate is uniform. In this zone, the diameter of the particle is constant and the bulk density is proportional to the mass loss rate [58, 83]. In Zone III, which occurs at high temperatures and for large particle, mass transport of species from the bulk to the surface of the particle is the controlling step. The reactions occur on the external surface of the particle, the bulk density is approximately constant, and the diameter of the particle reduces at a rate proportional to the one-third power of the mass loss. Finally, in Zone II, both mass transport and chemical reaction are important. In this zone, there is a partial penetration of oxidizer in the particle, leading to gradients of bulk density, specific area, and mass loss. Both diameter and bulk density are changing in this zone. In a real industrial boiler, most of the particles are burned under Zone II conditions. This three-zone approach is a simplified model to understand the char oxidation phenomenon. However, a more realistic approach requires that other phenomena besides mass transport and chemical reaction be considered, including char structure variations, particle size effects, radiation between particles, changes in specific surface area, fracturing of char, catalytic reactions of ash, evaporation of ash, and so on. 78 Erland et al. [18] developed a comprehensive model for char particle conversion for a volume containing multiple particles. Their model predicts the temporal variations of mass, bulk density, and diameter for a porous particle that is surrounded by a gas mixture of O2 and CO2 . A detailed nine-step heterogeneous reaction mechanism is used to account for the conversion of carbon. A homogeneous reaction mechanism (GRI-Mech 3.0) is used for the gas reactions. They assume that the heterogeneous reactions are occurring on the surface of the particle and then use a Thiele modulus strategy to correct for reactant penetration. Additionally, they take into account the particle-wall and particle-particle radiation effects. Equations for the mass of gas, the mass of each gas species, and the energy are solved in a hermetic volume with a fixed number of particles. In the energy equation, the heat of reaction as well as the convective and conductive heat fluxes are considered. For the particle phase, an equation for the mass of the particle is solved; this model assumes that there is no evaporation of ash. The mass loss of the particle is computed as the product of the carbon consumption per unit of surface area multiplied by the specific surface area. The specific surface area as a function of carbon conversion is calculated using an expression from Bhatia and Perlmutter [7]. Mitchell et al. [58] developed a detailed model that predicts the physical changes of a particle during the combustion process. In this approach, the particle is divided into a number of concentric, annular volume elements. For each element, the bulk density, mass loss rate, and specific area are computed. The model predicts the particle's burning rate, temperature, diameter, bulk density, and surface area. Differential equations for the carbon mass of the particle, the O2 concentration, and the fraction of the total sites having adsorbed oxygen atoms are solved. The model assumes an isothermal particle, and homogeneous gas reactions are not considered either inside of the particle or in the surrounding gas. However, a mass balance is used to compute the flux of O2 from the surrounding gas to the surface of the particle. A six-step heterogeneous reaction mechanism is used to compute the consumption of carbon. The main purpose of this model is to accurately predict the bulk density, diameter, and specific area of the particle during its lifetime. This model can be considered a comprehensive or high-fidelity model, where information about the diffusion process, the heterogeneous reactions, and physical properties of the particle can be obtained and further used to produce a low-fidelity model. 79 Singer and Ghoniem [80] developed a comprehensive, transient, one-dimensional char model. This model accounts for evolution of a multimodal pore structure, computes fluxes using a multicomponent approach, and incorporates reaction annealing, particle shrinkage, and ash fragmentation. The model uses three heterogeneous reactions and two homogeneous reactions, although the homogeneous reactions are confined to the boundary layer. The model has the option of a time-dependent boundary condition. Singer and Ghoniem compared their model output with experimental data for a synthetic char collected in an entrained flow reactor at Sandia National Laboratories [34]. This highly detailed model focuses on the porous matrix evolution and the transport of gas species in the matrix. Its predictions of particle temperature and conversion closely match the data with discrepancies in the position closest to the injection point. The Surface Kinetics in Porous Particles (SKIPPY) model solves differential equations for the mass fractions of species and for energy both inside the particle and in its surroundings. SKIPPY was developed by Brian Haynes and coworkers at the University of Sydney[3] and has subsequently been used by Molina and coworkers [60] and by Hecht and coworkers [31-34]. It is similar in approach to Mitchell and coworkers [58]. However, Mitchell uses an effective diffusion coefficient for oxygen while in SKIPPY, a MaxwellStefan approach is used and multicomponent mass transfer effects are considered. There are several differences between SKIPPY and the global reaction model shown in equation (4.1). In SKIPPY, the whole particle is resolved and carbon consumption throughout the particle is computed ( e.g. with a six-step surface reaction mechanism [34]). In this paper, we present the Reacting Particle and Boundary Layer (RPBL) model. This model solves conservation equations inside the particle and in its surroundings (boundary layer). It includes heterogeneous surface reactions, gas phase reactions, and the evolution of physical properties with carbon burnout. This model is a transient version of SKIPPY with some modifications. We use the same approach as SKIPPY for computing species transport. The source terms are equivalent, but RPBL uses Cantera [30] to compute them while SKIPPY uses CHEMKIN [74] and surface CHEMKIN [12]. RPBL allows a timedependent boundary condition for all the equations instead of the constant boundary used in SKIPPY. RBPL has two energy equations, one for the particle and one for the gas, while SKIPPY uses a combined energy equation for the gas and particle. Additionally, RPBL 80 has a momentum equation to compute the velocity of the particle separately from the gas velocity. RPBL also accounts for a flowing boundary condition by adding a correction to the stagnant boundary condition used in SKIPPY; see equation (4.9). SKIPPY uses Darcy's law to compute the pressure inside of the particle; in RPBL, we assume a negligible gradient of pressure. The reason we do not use Darcy's law is that it is a simplification of the momentum equation for steady state flow in a porous medium, and we are resolving a transient state. In the RPBL model, we focus our efforts on the detailed solution of both heterogeneous and homogeneous reactions. We use a simple approach to compute the void fraction and the evolution of the solid matrix. In contrast, Singer and Ghoniem [80] focus their efforts on the detailed resolution of the solid matrix and use a simpler reaction model than our approach. We think the evolution of the solid matrix is an important process in char oxidation. However, adding model complexity to RPBL increases the number of uncertain parameters. Consequently, in a companion paper, we perform a validation and uncertainty quantification study to see if we can represent the experimental data collected by Hecht [34] with the current formulation of RPBL. We will then decide whether to add complexity to the evolution of the solid matrix. 4.3 Development of the RPBL model The RPBL model is a one-dimensional model in the radial direction that computes the transient-state conditions for a spherical, reacting, constant-diameter, porous particle and its reacting boundary layer. The particle and its surrounding gas are considered to be continuous media and there is transport of gas species between them (see Figure 4.1). Therefore, the particle voids are full of gas. The porous particle surface has active carbon sites which react heterogeneously with gas species, including O2 , CO2 , and H2 O. These heterogeneous reactions are presented in Table 4.1. Additionally, homogeneous reactions can occur between the gas species. The model uses a syngas (H2 /CO) reaction mechanism with 11 species (O2 , H2 O2 , CO, CO2 , O, H, OH, HO2 , HCO, H2 , H2 O) developed by Ranzi and coworkers [69]. Because of the carbon molecules that are being gasified and leaving the porous matrix, the physical properties of the particle are changing during the reaction process. 81 4.3.1 Species mass fraction equations The mass fraction for a species k in cell i (Yk,i ) is computed from equation (4.2). This equation is obtained with a mass balance on a cell i of volume Vi . Here, ρt,i is total mass density computed with the ideal gas equation, the subindeces i − 1/2 and i + 1/2 correspond to values at the left and right faces of the cell i, A is the area of the cell face (i − 1/2 or i + 1/2), nk is the mass flux of species k relative to stationary coordinates at the cell face, and the last two terms on the right-hand side are source terms from the heterogeneous and homogeneous reactions, respectively. dYk,i 1 ( Ank )i+1/2 − ( Ank )i−1/2 + Vi ṡk,i σr,i Wk + Vi ẇk,i Wk φi =− dt Vi ρt,i (4.2) The value for nk in equation (4.2) is computed with equation (4.3) using a Maxwell-Stefan multicomponent approach [86]. Here, jk is the mass diffusion flux relative to the mass average velocity and nt is the total mass flux relative to stationary coordinates. nk = jk + Yk nt (4.3) The value for jk is computed with equation (4.4). The total density (ρt ) is computed at the cell face using a linear approximation between i − 1 and i for the left face (i − 1/2) and o is the matrix of diffusion between i and i + 1 for the right face (i + 1/2). The parameter Dk,l coefficients [D o ] computed using the Maxwell-Stefan multicomponent approach [86] as shown in equation (4.5). In this equation, the matrix [D ] is the inverse matrix of the matrix B (computed with equation (4.6)), [W ] is the diagonal matrix of mass fractions, [X ] is the diagonal matrix of molar fractions, and the [B ou ] matrix is computed from equation (4.7). The gradient of mass fraction ( ∂ (Y ) ∂r ) is computed with a linear approximation; the gradient at face (i − 1/2) uses node values at i − 1 and i while the gradient at face (i + 1/2) uses node values at i and i + 1. The void fraction (φ) is computed as the combustion process evolves (see equation (4.28)) and the tortuosity (τ) is assumed constant [26]. The ratio φ τ is only used inside the particle. Equation (4.4) computes K g − 1 mass diffusion fluxes where K g is the total number of gas species. The mass flux for the last species (jKg ) is computed with equation (4.8) because the sum of all jk has to be zero. jk = −ρt K g −1 ∑ l =1 o Dk,l ∂(Yl ) φ ∂r τ (4.4) 82 [D o ] = [B ou ][W ][X ]−1 [D ][X ][W ]−1 [B ou ]−1 (4.5) K g xk xl Bkk = + ∑ - k,Kg - kl D D l =1k 6=l 1 1 − Bkl = − xk - kl D - kKg D YKg xk ou Blk = δlk − Yl 1 − xKg Yk jKg = − K g −1 ∑ jk (4.6) (4.7) (4.8) k =1 - ) are obtained using The binary Maxwell-Stefan multicomponent diffusion coefficients (D - ∗ must be corrected because the gas surrounding Cantera. However, the Cantera output D the particle is flowing and the Cantera model assumes that it is stagnant. The correction is presented in equation (4.9). In this equation, Shk,l is the Sherwood number computed with equation (4.10). The Schmidt number (Sck,l ) is in turn computed using equation (4.11), where the dynamic viscosity (µ g ) and the density (ρ g ) of the gas are computed in each cell. The local Reynolds number Re is computed from the difference between the gas (v g ) and particle (v p ) velocities (see equation (4.12)). The correction factor (0.5Shk,l ) is 1 when the difference between gas and particle velocities is zero and increases as the velocity difference increases. ∗ - k,l = 0.5Shk,l D - k,l D Shk,l = 2 + Re1/5 Sc1/3 k,l Sck,l = Re = µg - k,l ρ gD 2r p |v g − v p |ρ g µg (4.9) (4.10) (4.11) (4.12) The value for the total mass flux, nt in equation (4.3) at the right face (i + 1/2) of each cell is computed with equation (4.13). This equation was derived using conservation of the total mass, assuming that dρ¯t dt ≈ 0. For cell i, the total mass flux at the left face is equal to the 83 right face mass flux for the cell i − 1. Thus, the right face mass flux is computed using the left face mass flux of the cell i − 1 beginning from i = 0 where ( Ant )−1/2 = 0. Kg ( Ant )i+1/2 = ( Ant )i−1/2 + Vi ∑ ṡk,i σr,i Wk (4.13) i =1 The production or destruction of species k in cell i is estimated with the two source terms in equation (4.2). The production or destruction at the carbon surface is represented by the term ṡk,i σr,i Wk . In this term, the molar production rate (ṡk,i ) for each species k is computed from a general surface kinetics formalism described in [50] and adapted in Cantera [30] that uses the surface reaction mechanism presented in Table 4.1; σr,i = Sgc,i ρbulk c,i is the product of the specific surface area (Sgc,i ) and the particle bulk density of carbon (ρbulk c,i ); and Wk is the molecular weight of species k. The second source term ẇk,i Wk φi is the production or destruction of species k due to gas phase reactions. The molar production rate (ẇk,i ) is computed from Cantera using the syngas reaction mechanism from Ranzi and coworkers [69] and the void fraction (φ) is computed from equation (4.28). The ρbulk c,i term is computed from equation (4.14), which is the carbon mass balance for cell i. Here, Wc is the molecular weight of carbon and ṡcb ,i is the molar destruction rate of bulk carbon. The Sgc,i term is computed with the Bhatia model [7, 58] presented in equation (4.15). In this model, Sgc,initial is the initial specific surface area, which is assumed to have a constant value throughout the particle. The conversion of carbon (xc ) is computed from the current and initial carbon bulk densities, xc,i = 1 − ρbulk c,i ρbulk c,initial . The structural parameter ψ has a value between 3 and 8 [18, 58]. d(ρbulk c,i ) = ṡcb ,i Wc σr,i dt Sgc,i = Sgc,initial (1 − xc,i )(1 − ρbulk c,i ρbulk c,initial )(1 − ψln(1 − xc,i ))1/2 (4.14) (4.15) The fixed boundary condition for equation (4.2) at rin f is the mass fraction of species k in the bulk gas. Symmetric boundary conditions are specified at r = 0, the center of the particle, jk,−1/2 = 0 and nt,−1/2 = 0. 4.3.2 Energy equations The gas temperature in cell i (Tg,i ) is computed from a gas phase energy balance on the cell as shown in equation (4.16). All the physical properties of the gas such as heat capacity 84 (c p g,i ), enthalpy (hk,i ), and thermal conductivity (λ g,i ) are computed with Cantera at the gas temperature and composition of cell i using the transport data from Ranzi et al [69]. The energy flux due to conduction (Fcond ) is computed with equation (4.17) where the gradient of gas temperature ( dTg dr ) is computed with a linear approach in the same way that the gradient of mass fraction is computed. The energy flux due to the transport of enthalpy Kg (Fh ) is computed with equation (4.18). The term ρt,i ∑k=1 hk,i dYk,i dt is a transient term due to the rate of change of the species mass fractions [51]. The heat exchanged between the gas inside the particle and the solid matrix of the particle (Si ) is computed with equation (4.19). We assume that the heat flux is proportional to the difference in temperature between the solid and the gas. The proportionality constant, hsolid gas , is a convection coefficient. Since the parameter hsolid gas is hard to measure or compute, we choose a value (hsolid gas = 1 [ W K ]) that forces the model to have equal particle and gas temperatures in each cell i. h dTg,i 1 − [( AFcond )i+1/2 − ( AFcond )i−1/2 ] − [( AFh )i+1/2 − = dt ρt,i c p g,i Vi Kg ( AFh )i−1/2 ] − Vi ρt,i Fcond = −λ g dTg dr ∑ k =1 hk,i (4.16) i dYk,i + Si dt (4.17) Kg Fh = ∑ nk hk (4.18) k =1 Si = hsolid gas ( Tp,i − Tg,i ) (4.19) The boundary condition at rin f is Tg,in f = Tg,bulk where Tg,bulk is the gas temperature in the bulk flow. In the center of particle (r = 0), Fcond,−1/2 = 0 and Fh,−1/2 = 0, which means the heat flux due to conduction or to the transport of enthalpy is zero at the left (−1/2) face in the first cell. The particle temperature in cell i (Tp,i ) is computed from a particle phase energy balance on the cell as shown in equation (4.20). In the equation, ρbulk p,i and c p p,i are the bulk density and heat capacity of the particle at cell i, respectively. The energy flux due to conduction (Fcond,p ) is computed with equation (4.21) where λ p is the thermal conductivity of the particle computed with equation (4.22). The thermal conductivity of the solid matrix (λs ) in equation (4.22) is assumed to be constant during the combustion process. The gradient of particle temperature dTp dr is computed with a linear approximation in the same 85 way that dTg dr is computed. The term hcb ,i ṡcb ,i σr,i is the production of energy due to the transformation of bulk carbon to the gas; hcb ,i is the enthalpy of solid carbon and ṡcb ,i is the carbon consumption rate. The Si term is computed with equation (4.19). i h dTp,i 1 = − [( AFcond,p )i+1/2 − ( AFcond,p )i−1/2 ] − hcb ,i ṡcb ,i σr,i Vi − Si dt ρbulk p,i c p p,i Vi (4.20) dTp dr (4.21) λ p,i = λs (1 − φi ) (4.22) Fcond,p = −λ p The ρbulk p,i value in equation (4.20) is computed with equation (4.23). It is assumed that no ash leaves the particle and that the ash mass in each cell (m ash,i ) is constant. The c p p,i value is computed from equation (4.24) where Yc,i is the mass fraction of carbon (see equation (4.25)) and the heat capacities of char (c p char,i ) and ash (c p ash,i ) are computed with Merrick's mathematical model [57]. ρbulk p,i = ρbulk c,i Vi + m ash,i Vi c p p,i = Yc,i c p char,i + (1 − Yc,i )c p ash,i Yc,i = ρbulk c,i Vi ρbulk c,i Vi + m ash,i (4.23) (4.24) (4.25) The boundary condition at r = 0 is Fcond,p,1 = 0. The boundary condition at r = r p is more complex because of heat exchanged between the gas and the particle and between the particle and its surroundings. We are modeling a particle in an entrained flow reactor at Sandia National Laboratories [34] as described in section 4.4. Thus, at r = r p , Fc,p = − Finter f ace . The value for Finter f ace is computed with equation (4.26), where k n is the convection coefficient computed with equation (4.27), Re is the Reynolds number, and Pr is the Prandtl number. The particle and gas temperatures at r = r p are Tp,r p and Tg,r p , respectively. The second term in equation (4.26) is the heat exchanged due to radiation from the wall of the reactor; Tw is the reactor wall temperature, and ε p is the emissivity of the particle surface. 4 Finter f ace = −k n ( Tp,r p − Tg,r p ) − ε p σ( Tp,r − Tw4 ) p kn = λg (2 + 0.6Re1/2 Pr1/3 ) 2r p (4.26) (4.27) 86 With ρbulk p,i defined, the void fraction (φ) in equation (4.2) is computed with equation (4.28); ρtrue p,i is the true density (see equation (4.29)) computed from the true density of carbon (ρtrue c ) and the true density of ash (ρtrue ash ). φi = 1 − ρbulk p,i ρtrue p,i 1 − Yc,i Y 1 = c,i + ρtrue p,i ρtrue c ρtrue ash 4.3.3 (4.28) (4.29) Momentum equation The composition and temperature of the bulk gas are boundary conditions for the species mass fraction equations (4.2) and the gas temperature equation (4.16). The bulk gas conditions need to be computed for a specific application. Our specific application is a set of char oxidation data from Sandia National Laboratories [34] collected in the entrained flow reactor described in section 4.4. We performed simulations of this reactor for different experimental configurations and then extracted from these simulations the bulk gas profiles for mass fraction, temperature, density, dynamic viscosity, and gas velocity along the particle trajectory. For additional details about the reactor simulations, see Dı́az-Ibarra [16]. To determine the appropriate bulk phase properties, we need to know where the particle is located in the reactor. Knowing its location, we can then interpolate the bulk gas boundary conditions. We do this by adding a momentum equation for the particle, presented in equation (4.30), to compute the velocity of the particle (v p ). The first term is the drag force contribution [56] and the second term is the buoyancy and gravity contribution. The bulk gas velocity (v g ) is imported from the reactor simulation output, F1 is computed with equation (4.31), which is a correlation to compute the drag coefficient [56], Re is computed with equation (4.12), the bulk gas density (ρ g,b ) and bulk gas dynamic viscosity (µ g ) are imported from the reactor simulation output, and tr is computed with equation (4.32). The total bulk density of the particle (ρbulk p,t ) is computed with equation (4.33), where Vt is total volume of the particle. In equation (4.30), the effect of the reduction of particle mass is accounted for through ρbulk p,t , which is decreasing with carbon consumption. The main intent of adding an equation for particle velocity is to compute the range of particle velocities seen in the experimental data of Hecht [34]. 87 ρ dv p F1 g,b = (v g − v p ) + g −1 dt tr ρbulk p,t (4.30) F1 = 1 + 0.15Re0.687 (4.31) tr = 4ρbulk p,t r2p 18µ g ρbulk p,t N ∑i 1 (ρbulk c,i Vi + m ash,i ) = Vt 4.3.4 Solution methodology (4.32) (4.33) We converted the above ordinary differential equations (ODEs) from time derivatives to length of the reactor derivatives using the chain rule and v p as shown in equation (4.34). dYk,i dY 1 = k,i dl dt v p d(ρbulk c,i ) d(ρbulk c,i ) 1 = dl dt vp (4.34) dTg,i dTg,i 1 = dl dt v p dTp,i 1 dTp,i = dl dt v p dv p dv p 1 = dl dt v p This set of ODEs is solved using the LSODA [1, 9] ODE solver. The number of ODEs depends on the number of species in the gas reaction mechanism and the number of cells in the computational domain. For the syngas (H2 /CO) reaction mechanism from Ranzi and coworkers [69], we need to solve equations (4.2) for the following 11 species: O2 , H2 O2 , CO, CO2 , O, H, OH, HO2 , HCO, H2 , and H2 O. This reaction mechanism has 33 reactions; we deleted species and reactions related to N2 . Thus, the number of ODEs we need to solve is 11 times N cells. Additionally, we must solve N ODEs for Tg (see equation (4.16)), Np (number of cells inside the particle) ODEs for Tp (see equation (4.20)), Np ODEs for ρbulk c,i (see equation (4.14)), and one ODE for v p (see equation (4.30)). In total, we must solve (12N + 2Np + 1) ODEs. We used a nonuniform grid size both inside and outside of the particle and set up the grid in such a way that r = r p is located at a cell face. We divided the volume inside the particle into Np cells of equal volume. Hence, the grid resolution is finer close to r = r p 88 and coarser near r = 0. We employed a log scale to locate the grid cells outside of the particle. The number of cells close to r = r p is large and decreases toward r = rin f . The total number of cells in the boundary layer is N − Np . 4.4 Reactor description Hecht [34] carried out a set of char oxidation experiments in an optically accessible, laminar, entrained flow reactor at Sandia National Laboratories (see Figure 4.2). Three types of coal were used to produce char, a high volatile bituminous coal (Illinois # 6 ), a western bituminous coal (Utah Skyline), and Black Thunder subbituminous coal (Power River Basin). The char was sieved into 6 narrow size bins of 53-63 µm, 63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm, and 125-150 µm. Twelve different environments were tested: O2 = 24, 36, 60 vol%, H2 O = 10, 14, 16 vol%, and the balance either N2 or CO2 . In these experiments, the particle size, velocity, and temperature of approximately 100 particles were obtained at 3-7 positions above the injection point (position 0 [m] in Figure 4.2. For the bulk gas condition in RPBL, we considered two different experimental environments: O2 = 24, 36 vol %, H2 O = 14 vol %, and balance CO2 . From computational fluid dynamics simulations of these two environments, we obtained bulk gas profiles for temperature, composition, and velocity [16]. The profiles were extracted from position 0 to position 0.254 m, the range of heights from which particle data were collected. 4.5 Model verification While we could not perform a verification study for all sets of equations used by RPBL, we did verify some components of the code. We used the method of manufactured solutions (MMS) [63] to verify the diffusion term in the species mass fraction equation (equation (4.2)). With MMS, a manufactured solution is applied to a differential equation. Because this solution is not the exact solution, the remainder needs to be coded as a source term. In this case, the differential equation was equation (4.2) with only the diffusion term. When this method was used, the solution converged to the manufactured solution. However, this verification was done with constant diffusion coefficients and only for the diffusion term. Due to the discontinuity between the particle and the gas, MMS could not be used when the solid and gas phases were coupled. 89 Next, we performed code-to-code comparisons between RPBL and SKIPPY [31, 34] to check for correct units and the correct use of the reaction mechanism. For these comparisons, we used the same surface reaction mechanism (see Table 4.1) and the same physical properties related to the carbon solid in both models. Then we compared numerical values for the source terms in the species equation (equation (4.2)) and the gas and particle phase energy equations (equation (4.16) and equation (4.20), respectively). We also ran a preliminary version of RPBL with constant properties (SKIPPY assumption) and compared the steady state solution with the SKIPPY solution; the results were the same. We checked grid independence for the two dependent quantities used in our sensitivity analysis, average particle temperature and particle velocity. We ran RPBL with 5 grid resolutions inside the particle: Np = 15, 20, 30, 50, 100. The same number of cells was used external to the particle for a total cell count of N = 30, 40, 60, 100, 200. We tested three particle sizes (diameter = 50 µ m, 80 µ m, 200 µ m). We obtained similar results for the three sizes, so we only present results for the 80 µ m case in Figure 4.3 and Figure 4.4. To generate these data, we ran RPBL and saved data every 0.0127 m from 0 (injection point) to 0.254 m above the injection point. We then computed the average temperature of the particle at each location using the nodes inside the particle. In this figure, the x axis is the position in the reactor, the y axis is the average particle temperature, and the dot colors correspond to the various grid resolutions. The average temperatures with Np = 15 are similar to those with Np = 100 with some small deviations at longer distances. We conclude that if we set Np = 20 (N = 40) or greater, we will obtain an accurate result for the average temperature. In Figure 4.4, we present the results of the particle velocity grid study. The effect of grid resolution is negligible. 4.6 Results In this paper, we present detailed results for a base case to show what type of data is produced from RPBL. We also present a sensitivity analysis that considers the effect of several uncertain model parameters on two quantities of interest measured by Hecht [34], particle temperature and particle velocity. The main objective of the sensitivity analysis is to select the most sensitive model parameters for a further validation analysis. For the base cases presented in this section, we ran RPBL with N = 60 cells (Np = 30), 90 which required the solution of 781 ODEs. For the all the cases in the sensitivity analysis, we ran RPBL with N = 40 cells (Np = 20). We ran the simulation until the particle reached a position of 0.254 m from the injection point. 4.6.1 Input and uncertainty map The RPBL parameters are presented in the form of the input and uncertainty (I/U) map in Table 4.2. Each parameter in the table is assigned a priority that represents the importance of the parameter on the quantities of interest (particle velocity and particle temperature). There are nine parameters with the highest priority of 6. The highest priority parameters all have ranges while the lower priority parameters are assigned nominal values. Model parameter values in Table 4.2 were selected from external sources and from nominal values used by SKIPPY [31]. The lower bound of 3 for τ is from Singer and Ghoniem [80]. The upper bound of 6 was modified from the SKIPPY value of 5 [31]. The rin f rp parameter has a value of 100 in SKIPPY [31], so we varied it from 50 to 120. The ψ parameter was determined from surface area data; we used the value of 3 from Mitchell and coworkers [58] as the lower bound and the value of 8 from Erland and coworkers [18] as the upper bound. While Singer and Ghoniem [80] and SKIPPY [31] used a value for e p of 0.85, we used a wider range (0.1 to 1.0) for our sensitivity analysis. For λ p , the range was bounded by the maximum and minimum values of the experimental data for ash deposits presented by Rezaei and coworkers [72]. We assumed that the nominal value for ρtrue c is 921 kg . m3 We based the nominal value for ρtrue ash on experimental correlations for the light components of ash [85]. The scenario parameters in Table 4.2 were obtained from a variety of sources. The d p range was based on experimental particle diameters reported by Hecht [34]. The parameters with a nominal value listed as "From reactor model" correspond to bulk gas parameters obtained from simulations of the entrained flow reactor at Sandia National Laboratories [34] described in section 4.4. The initial species mole fractions are based on the fact that the carrier gas in the particle feeding system system was CO2 . The nominal value used by SKIPPY [31] for φ was 0.4. We used a wider range of φinitial for the sensitivity analysis presented below. Due to the heterogeneous nature of the char particles, we chose 91 2 2 ] a wide range for Yc,intial . For Sgc,initial , we computed a value of 1e4 [ mkg ] from σr = 5.6e6 [ m m3 kg (σr = Sgc ρbulk c ) [31] and ρbulk c = 560 [ m3 ]. We added 20 % to this value to obtain the upper bound of Sgc,initial and removed 20 % to obtain the lower bound. The value for Tw reflects the approximate temperature of the quartz walls of the Sandia reactor. The pressure is the ambient pressure. 4.6.2 Base case For this particular case, we assigned the following values to these parameters: τ = 5, rin f rp = 53, ψ = 8, ε p = 0.96, λ p = 1.33, d p = 95 µm, φinitial = 0.18, Yc,initial = 0.98, and Sgc,initial = 8000 [kgc /m2 ]. For the other parameters listed in Table 4.2, we used the nominal value. The environment for this base case was O2 = 24 vol%, H2 O = 14 vol%, and balance CO2 . The initial conditions (t = 0) for the system of ODEs (see equation (4.34)) were as follows. For the cells internal to the particle, we assumed the particle was filled with CO2 because in the Sandia experiments, the carrier gas in the particle feeding system system was CO2 . This initial condition corresponds to the nominal values listed in Table 4.2. We computed the initial particle bulk density of carbon (ρbulk c,initial ) with equation (4.35) and assumed it was constant throughout the particle. In equation (4.35), we computed ρtrue p,initial from equation (4.29) with Yc,initial = 0.98. The initial temperature of the particle (800 K from the reactor model) was also assumed constant throughout the particle. The gas temperature in the cells inside the particle was assumed to be the same as the particle temperature. For the cells external to the particle, we assumed that xCO2 = 0.99 and that the initial gas temperature was equal to the bulk gas temperature at position 0 m (from the reactor model). ρbulk c,initial = ρtrue p,initial (1 − φinitial )Yc,initial (4.35) Even though the equations are cast in terms of mass fractions, we present the results here in molar fraction because it is easier to visualize low molecular weight species. The mole fraction results for the entire domain are presented in Figure 4.5. Because it is difficult to see the cells inside the particle cells in this figure, the mole fraction results for these cells are shown in Figure 4.6. Each plot in Figure 4.5 and Figure 4.6 corresponds to one of the 11 gas phase species. The x axis is r rp , which is the ratio of the distance from the 92 center of the particle to the radius of the particle. The y axis is the position of the particle in the reactor (distance above the injection point). The plots of O2 mole fraction are located in the first column and first row in Figure 4.5 and Figure 4.6. The mole fraction of O2 is ∼ 0 in the whole domain at position 0 m. From position 0 to ∼ 0.02 [m], O2 is transported from the bulk gas to the boundary layer and into the particle. In this initial stage, the particle is heating up and there are no chemical reactions (see CO plot in third column, first row and H2 plot in first column, fourth row in Figure 4.6; both are heterogeneous reaction products). Chemical reactions are initiated above position ∼ 0.02 [m]. From position ∼ 0.02 [m] to ∼ 0.03 [m], the O2 that has diffused into the particle is rapidly consumed (essentially zone III burning mode), so O2 mole fraction inside the particle drops back to 0 and reaction products such as CO and H2 can be seen. Above position ∼ 0.03 [m], carbon is consumed in the cells close to r rp = 1 (see Figure 4.7). As carbon bulk density is reduced, O2 diffuses further into the particle (zone II burning mode). The O2 in the boundary layer (Figure 4.5) shows a sharp gradient close to r rp = 1 due to the consumption of O2 inside the particle. The main products of the heterogeneous reaction mechanism, CO and H2 (see Table 4.1), show similar behavior. Both species are produced inside the particle and then react in the gas phase to produce species such as O (second column and second row ), OH (first column and third row), and CO2 (second column and second row). Consequently, the gas phase reactions have an impact on the heterogeneous production of CO and H2 . The particle carbon bulk density (ρbulk c ) is presented in Figure 4.7. In the early stages of particle heating (before ∼ 0.03 [m]), there is no carbon consumption. Then, from ∼ 0.03 [m] to ∼ 0.05 [m], there is a uniformly small reduction of ρbulk c throughout the particle due to the initiation of reactions with O2 . As char oxidation proceeds (above ∼ 0.05 [m]), the consumption of ρbulk c begins near rp rin f = 1 and then moves into the particle. The temperature inside the particle (Tp ) is presented on the left side of Figure 4.8. The temperature is relatively constant throughout the particle, so an assumption of an isothermal particle used by Mitchell and coworkers [58] is reasonable. This figure also illustrates the particle heating process. At position 0 [m], the entire particle is at a low temperature. From position 0 [m] to ∼ 0.02 [m], the particle is only heated by the gas. Above position ∼ 0.02 [m], the particle temperature increases rapidly, reaching a maximum of ∼ 1900 [K ] 93 at position ∼ 0.1[m]. This rapid increase indicates the onset of heterogeneous reactions. The particle temperature then decreases above position ∼ 0.12[m] due to the decreasing bulk gas temperature in the entrained flow reactor. The gas temperature (Tg ) both inside the particle and in its boundary layer is presented on the right side of Figure 4.8. The bulk gas temperature is the boundary condition for RPBL, and its effect on Tg far away from the particle is visible in this figure. The gas temperature is also affected by the particle. When the particle is cold (0 [m] - ∼ 0.02 [m]), the boundary layer has a temperature gradient from cold (particle) to hot (bulk gas). When the particle reaches a maximum temperature (0.08 [m] - 0.12 [m]), there is a large temperature gradient in the gas near the particle. Thus, under the conditions of this case, the boundary layer has an impact on the particle and the particle has an effect on the boundary layer. 4.6.3 Sensitivity analysis The primary objective of this sensitivity analysis is to determine the parameters that have the largest effect on the quantities of interest, particle temperature and particle velocity, for a subsequent validation analysis. The numerical, model, and scenario parameters required by RPBL are listed in Table 4.2 in the form of an I/U map, which includes the importance of the parameter to the quantities of interest (priority), the range of values for the parameter, and/or the nominal value of the parameter used in the simulations. In Table 4.2 there are nine parameters with the highest priority of 6; these nine parameters are the active parameters. We performed a sensitivity analysis with eight active parameters for three particle sizes: 50 µm, 80 µm, and 120 µm (particle size is the ninth parameter). We performed the sensitivity analysis using the Uncertainty Quantification Toolkit (UQTk), a set of C++ tools with a Python interface. It was developed by Debusschere and coworkers at Sandia National Laboratories [13]. UQTk uses the app pce_sens to compute total and main sensitivity indices using a polynomial chaos (PC) surrogate model. Using a first order PC surrogate model with a full quadrature rule, a total of 256 RPBL cases (2n where n is the number of dimensions or parameters) were needed for each particle size. The UQTk app generate_quad generated 256 quadrature points of ξ i = +/ − 0.577 with weights of w = 0.0625 for each dimension. The variable ξ i is 94 mapped to physical input space for the RPBL case using equation (4.36), where ai and bi are the bounds of the uncertainty intervals (parameter ranges) presented in the I/U map (Table 4.2). λi = b − ai a i + bi + i ξi 2 2 for i = 1, ..., n (4.36) RPBL produces data along the entire length of the reactor while the experimental data were taken at specific positions. Therefore, to create a PC surrogate model for each quantity of interest (particle temperature and particle velocity) at each experimental position, we extracted RPBL data at these positions, including particle velocity and Np temperatures in the particle. We then computed an average particle temperature from the Np temperatures. We selected the experiment with particles in the 53 − 60 µm range to represent the 50 µm diameter particle. In this experiment, data were taken at seven positions from 0.0375 m to 0.1125 m. For the 80 µm diameter particle, we selected the experiment with particles in the 75 − 90 µm range where data were taken at five positions from 0.05 m to 0.15m. For the 120 µm diameter particle, we selected data from the experiment with 106 − 125 µm diameter particles. These data were taken at six positions from 0.075 m to 0.2 m. Our next step was to compute the main sensitivity indices as defined by UQTk [13] at each experimental position for each quantity of interest. We used the UQTk app pce_sens with the coefficients of the PC surrogate model for each position. We present the main sensitivity indices for the 50 µm diameter particle in Figure 4.9 (particle temperature) and Figure 4.10 (particle velocity); this diameter represents small size particles. In these figures, the y axis is the main sensitivity index and the x axis is the position in the reactor. For particle temperature, the most sensitive parameter is φinitial ; this parameter is used to compute the density of the particle. A particle with a high φinitial is a light particle while a small φinitial indicates a heavy particle; in both cases, the carbon content is determined by Yc , the third most sensitive parameter. Since φinitial and Yc are controlling the amount of carbon in the particle, it makes sense that particle temperature has a high sensitivity to them, as the carbon density determines the duration of active char burning at high temperatures. The second most sensitive parameter is e p , which affects radiation heat loss of the particle. The high sensitivity is at least partially due to the wide uncertainty range of this parameter. The sensitivity index for this parameter is larger in the first and last 95 positions. In these positions, the particle temperature is lower and the effect of radiant heat loss is diminished. We expected high sensitivities for Sgc,initial and ψ, which are used in the Bhatia and Perlmutter specific surface area model [7]. However, their impact on the particle temperature is small. The velocity of the particle Figure 4.10 is only sensitive to φinitial and Yc . These parameters determine the bulk density of the particle, which in turn affects the gravity term in the momentum equation of the particle. The sensitivity analysis for the 80 µm diameter particle is presented in Figure 4.11 and Figure 4.12; this diameter represents medium size particles. For particle temperature, φinitial and Yc are the most sensitive parameters at some positions while e p is most sensitive for the 0.05 [m] and 0.075 [m] positions. The next most sensitive parameters are τ, rin f rp , and λ p . As with the 50 µm diameter particle, Sgc,initial and ψ have only a small impact on the particle temperature. The particle velocity sensitivities are similar to those for the 50 µm diameter particle; φinitial and Yc are the most sensitive parameters. The sensitivity analysis for the 120 µm diameter particle is presented in Figure 4.13 and Figure 4.14; this diameter represents large size particles. For particle temperature, the most sensitive parameter is e p followed by φinitial and Yc . The parameters rin f rp , τ, and λ p exhibit greater sensitivity for this particle size than for the smaller particle sizes, as expected based on the greater diffusional distance within the particle and greater flux of oxidation products from the particle surface. Similar to the 50 µm and 80 µm diameter particles, Sgc,initial and ψ have a small impact on the particle temperature. The particle velocity sensitivities also follow the trend of the other two particle sizes; φinitial and Yc are the most sensitive parameters. We learn the following from these sensitivity analyses at the three particle sizes. The particle velocities are only affected by φinitial and Yc . The three most sensitive parameters for particle temperature are φinitial , Yc , and e p . For the largest particles (120 µm), e p is the most sensitive parameter. The parameters rin f rp , τ and λ p have a discernible impact on particle temperature for medium and large particle sizes. However, Sgc,initial and ψ do not affect the particle temperature in the particle size range under study. 96 4.7 Conclusion We presented RPBL, a transient model for char oxidation in a spherical particle and in its boundary layer. This model accounts for mass transport of gas species between the particle and the boundary layer using a time dependent boundary condition. Heterogeneous and homogeneous reaction rates are computed with a six-step reaction mechanism and a syngas reaction mechanism, respectively. The burnout process is traced through the particle as it moves through the length of the entrained flow reactor at Sandia National Laboratories [61, 88]. We presented a base case for a 95 µm diameter char particle resolved by RPBL. We described in detail the mole fraction plots of the 11 species and the particle and gas temperature plots. Initially, the particle is cold and O2 diffuses inside of the particle. As reactions are initiated, the O2 mole fraction inside the particle drops to zero and a measurable O2 mole fraction is only seen near r rp = 1 (essentially zone III burning mode). As reactions progress, oxygen diffuses further into the particles (zone II burning mode). The model predicts a constant temperature throughout the particle. The time dependent boundary condition yields a boundary layer temperature that decreases with time (height in the reactor); this effect on the particle temperature is noted as well. We performed a sensitivity analysis for nine parameters on the particle temperature and velocity. We evaluated eight parameters directly; for the ninth parameter, d p , we performed the sensitivity analysis at three particle diameters. We performed the analysis in this way because the sensitivity to d p was so high that all other parameter were insignificant. For the particle temperature, φinitial , Yc , and e p are the most sensitive parameters. Three other parameters also have an impact: rin f rp , τ, and λ p . The velocity of the particle is sensitive to φinitial and Yc . Neither the particle temperature nor the velocity is sensitive to Sgc,initial and ψ. 4.8 Acknowledgements This material is based upon work supported by the Department of Energy, National Nuclear Security Administration, under Award Number DE-NA0002375. 97 Sandia National Laboratories is a multiprogram laboratory operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC0494AL85000. 98 Table 4.1. Surface reaction mechanism [31, 34] Reaction A[mol, cm2 , s] Cb + Cs + O2 → CO + C(O)s 3.3 × 1015 C(O)s + Cb → CO+ Cs 1.0 × 108 Cs + O2 → C(O2 )s 9.5 × 1013 C(O2 )s + Cb → CO2 +Cs 1.0 × 108 Cs + CO2 → CO +C(O)s variable Cs + H2 O → H2 +C(O)s variable E[kJ/mol] 167.4 0.0 142.3 0.0 251.0 222.0 Particle Bulk Flow rp rinf Boundary layer (Bulk conditions are inputs) Figure 4.1. Representation of the particle and boundary layer. 99 Table 4.2. Input/uncertainty map for RPBL Parameter Priority Range Nominal value min max Numerical Parameters Np 1 15 100 30 N 1 30 200 60 rtol 1 1e-4 Model Parameters τ [−] 6 3 6 rin f [−] 6 50 120 rp ψ [−] 6 3 8 e p [−] 6 0.1 1 W λ p [ mK ] 6 0.1 2 kg ρtrue c [ m3 ] 1 921 kg ρtrue ash [ m3 ] 1 2000 W hsolid gas [ m2 K ] 1 1 Scenario parameters d p [µm] 6 50 160 v g [ ms ] 3 From reactor model Tg [K ] 3 From reactor model O2,bulk [−] 3 From reactor model H2 Obulk [−] 3 From reactor model CO2,bulk [−] 3 From reactor model H2,bulk [−] 3 From reactor model CObulk [−] 3 From reactor model O2,initial [−] 3 1.00e-3 H2 Oinitial [−] 3 1.00e-3 H2,initial [−] 3 1.00e-3 COinitial [−] 3 1.00e-3 CO2,initial [−] 3 0.99 Tp,initial [K ] 3 From reactor model φinitial 6 0.15 0. 7 Yc,intial [−] 6 0.5 1 kg Sgc,initial [ m2 ] 6 8000 12000 Tw [K ] 3 500 kg Pressure [ m2 ] 3 1.00e5 100 Position 0 Char Feed Figure 4.2. Schematic of optically accessible entrained flow reactor at Sandia National Laboratories. Figure adapted with permission from [34] 2200 Average temperature [K] 2000 1800 1600 1400 Np = Np = Np = Np = Np = 1200 1000 800 0.00 0.05 0.10 0.15 Position [m] 0.20 15 20 30 50 100 0.25 Figure 4.3. Grid independence study for average particle temperature, d p = 80 µ m. 101 3.2 Particle velocity [m/s] 3.1 3.0 Np = Np = Np = Np = Np = 15 20 30 50 100 2.9 2.8 2.7 2.6 2.5 0.00 0.05 0.10 0.15 Position [m] 0.20 0.25 Figure 4.4. Grid independence study for particle velocity, d p = 80 µ m. 102 H2O2 CO 0.0000018 0.096 Position [m] 0.120 0.10 Molar fraction Position [m] 0.144 0.0000012 0.15 0.0000010 0.0000008 0.10 0 10 20 r rp 30 40 0.54 0.15 0.45 0.36 0.10 0.27 0.0000004 0.05 0.048 0.18 0.05 0.09 0.0000002 0.024 0.00 0.63 0.0000006 0.072 0.05 0.20 0.0000014 0.168 0.15 0.72 0.0000016 0.20 0.192 Molar fraction Position [m] 0.20 0.000 0.00 0 10 [−] 20 r rp CO2 30 40 0.0000000 0.00 0 10 [−] 20 r rp O 30 0.00 40 [−] H 0.00216 0.96 0.00280 0.87 0.20 0.20 0.00192 0.20 0.00245 0.00168 0.78 0.00175 0.00140 0.10 0.33 0 10 20 r rp 30 40 0.00120 0.00096 0.10 0.00072 0.00070 0.05 0.00048 0.05 0.00035 0.24 0.00 0.15 0.00105 0.42 0.05 0.00144 Position [m] 0.15 Molar fraction 0.51 0.10 Position [m] 0.60 Molar fraction Position [m] 0.00210 0.69 0.15 0.15 0.00 0 10 [−] 20 r rp OH 30 40 0.00024 0.00000 0.00 0 10 [−] 20 r rp HO2 0.20 30 [−] HCO ×10−7 5.6 0.0000144 0.20 0.0049 0.00000 40 0.0000162 0.0056 Molar fraction 0.216 Molar fraction O2 0.20 4.9 0.0000126 0.0000090 0.0000072 0.10 0.0021 0.0000036 0.05 0.0007 0.00 0 10 20 r rp 30 40 3.5 2.8 0.10 2.1 0.0000054 0.0014 0.05 0.15 1.4 0.05 0.7 0.0000018 0.0000 0.00 0 10 [−] 20 r rp H2 30 40 0.0000000 [−] Molar fraction 0.0028 0.10 4.2 0.0000108 0.15 Molar fraction Position [m] 0.0035 Position [m] 0.15 Molar fraction Position [m] 0.0042 0.00 0 10 20 r rp 30 40 0.0 [−] H2O 0.0135 0.130 0.0120 0.20 0.20 0.115 0.0105 0.15 0.085 0.070 0.10 0.055 0.0045 0.0030 0.05 0.040 0.05 0.025 0.0015 0.00 0 10 20 r rp 30 [−] 40 0.0000 Molar fraction 0.0060 0.10 Position [m] 0.0075 Molar fraction Position [m] 0.100 0.0090 0.15 0.00 0 10 20 r rp 30 40 0.010 [−] Figure 4.5. Species mole fractions in the particle and boundary layer of 95 µm diameter char particle. 103 O2 H2O2 CO 0.0000018 0.168 0.20 0.72 0.0000016 0.20 0.147 0.20 0.63 0.0000014 0.10 0.063 0.0000008 0.10 0.000 0.00 0.0 0.2 0.4 [−] r rp CO2 0.6 0.8 0.0018 0.15 0.0015 0.0012 0.10 0.42 0.8 0.15 0.00 0.0 0.2 0.4 [−] r rp 0.6 0.8 0.0024 0.10 0.0000108 0.15 0.0000090 0.0000072 0.10 0.0018 0.0000036 0.05 0.0006 0.6 0.8 0.0000 0.00 0.0 0.2 0.4 H2 r rp 0.6 0.8 2.8 2.1 0.0060 0.10 Position [m] 0.0075 0.0000000 0.7 0.00 0.0 0.115 0.15 0.085 0.070 0.10 0.055 0.0045 0.0030 0.05 0.040 0.05 0.025 0.0015 0.6 0.8 0.0000 1.4 0.05 0.100 Molar fraction Position [m] 3.5 0.10 0.130 0.20 0.0105 0.0090 4.9 0.15 H2O 0.0120 [−] 0.20 [−] 0.0135 0.15 ×10−7 0.0000018 [−] 0.20 [−] 0.0000054 0.0012 0.05 0.00000 0.8 4.2 Position [m] 0.0030 Molar fraction Position [m] 0.15 r rp 0.6 0.0000126 0.0036 r rp 0.4 5.6 0.0000144 0.20 0.4 0.2 HCO 0.0042 0.2 0.00 0.0 0.0000162 0.0048 0.00 0.0 0.00024 0.0000 HO2 0.20 r rp 0.00048 0.05 [−] 0.0054 0.4 0.00096 0.10 0.0003 OH 0.2 0.00120 0.00072 0.0006 0.05 0.24 0.6 0.00144 0.0009 0.33 0.05 0.00168 0.15 0.00 0.0 0.2 0.4 r rp 0.6 0.8 0.010 [−] Figure 4.6. Species mole fractions in 95 µm diameter char particle. 0.2 0.4 r rp 0.6 [−] 0.8 0.0 Molar fraction 0.10 0.00216 0.20 Position [m] 0.51 [−] 0.0021 Position [m] 0.60 0.00 0.8 0.00192 0.0024 0.20 Molar fraction Position [m] 0.69 0.15 0.00 0.0 r rp 0.6 H 0.78 r rp 0.4 0.0027 0.87 0.4 0.2 O 0.20 0.2 0.00 0.0 [−] 0.96 0.00 0.0 0.09 0.0000000 Molar fraction 0.8 0.0000002 Molar fraction r rp 0.6 0.18 0.05 Molar fraction Position [m] 0.4 0.36 0.10 Molar fraction 0.2 0.45 0.27 0.0000004 0.05 0.021 0.00 0.0 0.15 0.0000006 0.042 0.05 0.0000010 Molar fraction 0.084 0.54 0.0000012 0.15 Molar fraction Position [m] 0.105 Position [m] 0.15 Molar fraction Position [m] 0.126 104 kg Carbon bulk density [ m 3] 750 675 0.20 Position [m] 600 525 0.15 450 375 0.10 300 225 0.05 150 0.00 0.0 0.2 0.4 r rp 0.6 75 0.8 [−] Figure 4.7. Carbon bulk density inside 95 µm diameter char particle. Particle Boundary layer 1880 1880 1760 1760 0.20 1640 1640 1520 1280 0.10 Position [m] 1400 1520 0.15 Temperature [K] Position [m] 0.15 1400 1280 0.10 1160 1040 0.05 1160 1040 0.05 920 0.00 0.0 0.2 0.4 r rp 0.6 [−] 0.8 800 Temperature [K] 0.20 920 0.00 0 10 20 r rp 30 40 800 [−] Figure 4.8. Particle temperature (left) and particle and boundary layer gas temperature (right) for 95 µm diameter coal particle. 105 1.2 1.0 τ φinitial Yc,initial rinf rp Sgcinitial ψ ²p λp Sensitivity 0.8 0.6 0.4 0.2 0.0 0.0375 0.05 0.0625 0.075 0.0875 Position [m] 0.1 0.1125 Figure 4.9. Particle temperature sensitivity analysis for 50 µm diameter particle. 1.2 1.0 τ φinitial Yc,initial rinf rp Sgcinitial ψ ²p λp Sensitivity 0.8 0.6 0.4 0.2 0.0 0.0375 0.05 0.0625 0.075 0.0875 Position [m] 0.1 0.1125 Figure 4.10. Particle velocity sensitivity analysis for 50 µm diameter particle. 106 1.2 1.0 τ φinitial Yc,initial rinf rp Sgcinitial ψ ²p λp Sensitivity 0.8 0.6 0.4 0.2 0.0 0.05 0.075 0.1 Position [m] 0.125 0.15 Figure 4.11. Particle temperature sensitivity analysis for 80 µm diameter particle. 1.2 1.0 τ φinitial Yc,initial rinf rp Sgcinitial ψ ²p λp Sensitivity 0.8 0.6 0.4 0.2 0.0 0.05 0.075 0.1 Position [m] 0.125 0.15 Figure 4.12. Particle velocity sensitivity analysis for 80 µm diameter particle. 107 1.2 1.0 τ φinitial Yc,initial rinf rp Sgcinitial ψ ²p λp Sensitivity 0.8 0.6 0.4 0.2 0.0 0.075 0.1 0.125 0.15 Position [m] 0.175 0.2 Figure 4.13. Particle temperature sensitivity analysis for 120 µm diameter particle. 1.2 1.0 τ φinitial Yc,initial rinf rp Sgcinitial ψ ²p λp Sensitivity 0.8 0.6 0.4 0.2 0.0 0.075 0.1 0.125 0.15 Position [m] 0.175 0.2 Figure 4.14. Particle velocity sensitivity analysis for 120 µm diameter particle. CHAPTER 5 REACTING PARTICLE AND BOUNDARY LAYER (RPBL) MODEL: VUQ In preparation for submission to Combustion and Flame, Reacting Particle and Boundary Layer (RPBL) model: Validation and uncertainty quantification. Oscar H. Dı́az-Ibarra, Jennifer Spinti, Philip J. Smith, Christopher Shaddix, Ethan Hecht. 109 The Reacting Particle and Boundary Layer (RPBL) model solves species, energy, and carbon consumption equations for a char particle. We present a six-step validation and uncertainty quantification (VUQ) study of RPBL with experimental measurements of single particle temperatures and velocities. The quantities of interest (QOIs) in the analysis are particle temperature and velocity. We explore five RPBL model parameters (particle diameter, void fraction, carbon fraction of the particle, particle emissivity, and the nondimensional thickness of the boundary layer) to which the QOIs are most sensitive based on results from a previous sensitivity analysis. The experimental data were collected in an optically accessible, laminar, entrained flow reactor at Sandia National Laboratories. We selected data from three chars (obtained from Utah Skyline, Illinois # 6, and Black Thunder coals), two environments (O2 = 24, 36 vol %, H2 O = 14 vol %, and the balance CO2 ) and six size bins (53-63 µm, 63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm, and 125-150 µm). We performed 36 bound-to-bound consistency analyses to compare the simulation data with the experimental data. Because the consistency analysis required fast function evaluations of the simulation data, we constructed polynomial chaos surrogate models for the simulation data output from a suite of approximately 11,000 simulations. We found consistency for all 36 analyses. From these analyses, we learned that the particle diameter is the most important parameter for obtaining consistency. We were able to reproduce the particle temperature and velocity with the RPBL model. However, we did not reduce uncertainty in the parameters under study. 5.1 Nomenclature dp hsolid gas N Np P r rin f rp Particle diameter [µm] Convection coefficient, hsolid gas = 1 [ W K] Number of cells in entire domain Number of cells inside of particle Pressure [ Pa] Radius [m] Ratio of rin f to r p [−] Sgc T v Y e λp Specific surface area [ mkg ] Temperature [K ] Velocity [ ms ] Mass fraction [−] Particle emissivity [−] Particle thermal conductivity [ mWk ] 2 110 kg ρ σr τ φ ψ Density [ m3 ] Specific surface area for heterogeneous reactions, σr = Sgc ρbulk,c Tortuosity [−] Void fraction [−] Structure parameter [−] Subscripts ash: bulk: c: g: in f : initial: p: true: w: Ash in particle Bulk gas Carbon in particle Gas Limit of the particle domain Initial condition Particle True density Wall of the reactor 5.2 Introduction The Surface Kinetics in Porous Particles (SKIPPY) model was developed by Brian Haynes and coworkers at the University of Sydney [3], and subsequently used by Molina and coworkers [60] and by Hecht and coworkers [31-34]. Our original goal was to perform a validation and uncertainty quantification (VUQ) analysis with SKIPPY model outputs and experimental data collected by Hecht [34]. We performed the analysis for one data point (one location in the reactor), but we could not perform a simultaneous comparison using transient data that were collected along the length of the reactor. Consequently, we developed a new model form, the Reacting Particle and Boundary Layer (RPBL) model. This model, based on SKIPPY, computes the transient-state conditions for a spherical, reacting, porous char particle and its reacting boundary layer. While SKIPPY produces detailed data relative to the char oxidation process, it assumes the particle is at steady state. In contrast, RPBL is a transient model that includes a time-dependent boundary condition representing the char particle's movement through a specified reactor. The model formulation was previously described by Dı́az-Ibarra et al [15]. The RPBL model uses a Maxwell-Stefan multicomponent approach to compute the mass transport of species. A syngas (H2 /CO) mechanism is used to compute the homogeneous reaction rates; a six-step reaction mechanism is used to compute the heterogeneous reaction rates. RPBL uses a carbon density (burnout) equation to compute the evolution 111 of the carbon consumption. The physical properties of the solid matrix are computed from the carbon, ash, and void fractions. The void fraction is computed assuming a constant particle diameter throughout burnout. Two energy equations are solved, one for the particle and one for the gas. RPBL uses a time-dependent boundary condition that is obtained from the specific reactor where experimental data are collected, in this case an entrained flow reactor at Sandia National Laboratories [61, 88]. RPBL also has a momentum equation for computing the particle velocity along the length of the reactor. This paper focuses on a VUQ analysis of the RPBL model. The purpose of this analysis is to see if the RPBL model outputs can reproduce the experimental particle temperature and velocity data collected by Hecht [34]. The RPBL model is a high fidelity model because it solves differential equations to compute particle temperatures and gas phase composition in the particle and its boundary layer. However, this model includes many assumptions and is prescribed by a set of parameters. Therefore, it needs to be validated. Once validated, the RPBL model can be used to review assumptions of simpler models relating to the effect of the boundary layer, the effect of gas phase reactions, diffusion inside the particle, uniformity of particle temperature, and so on. 5.3 Description of VUQ approach In validation, we test a model using experimental data, otherwise known as the quantity of interest (QOI), as a reference [63]. We say that our model is validated if we can reproduce the experimental data. In uncertainty quantification (UQ), we explore uncertainty in the inputs of the model and determine the model's accuracy in calculating the experimental data. Our goal with the UQ analysis is to develop a model that computes the QOIs with the degree of accuracy required for the particular application. In this analysis, our VUQ methodology is a modified version of the Simulator Assessment and Validation Engine (SAVE) framework [6] as presented in Figure 5.1. This methodology, developed by Schroeder [77], consists of the six steps. In step 1, model output(s) are selected as evaluation criteria or QOIs. The input/uncertainty (I/U) map that is created in step 2 consists of a list of parameters (model, scenario, and numerical) and their uncertainties that may have an impact on the QOIs. Each parameter is assigned a priority; those with high priority are selected as active parameters and are investigated 112 further. Step 3, the collection of both experimental and simulation data, is closely tied to steps 1 and 4. The data that are collected must allow the direct measurement or the calculation of the QOIs. The requirements of the surrogate model determine the simulations that are performed. Step 4, surrogate model development, is required when running the computer model is expensive. In this work, we use a polynomial chaos (PC) surrogate model [13]. Various methodologies can be used to analyze model outputs in step 5 [52, 75]. We use a consistency analysis methodology referred to as bound-to-bound consistency [75]. The basic concept is that if the difference between model outputs and experimental data is less than the error in the experimental measurement, the simulation and experimental data are consistent. In step 6, feedback and feed-forward, the I/U map is reviewed based on the consistency analysis outcome, models are modified, parameter uncertainties are reevaluated, and the evaluation criteria are reviewed to see if new data should be incorporated. We present the application of this six-step methodology to the RPBL model output and the Hecht experiments [34] in the following sections. 5.4 Selection of quantities of interest The QOIs for this analysis are the particle temperature and velocity data for three char types (char were obtained from Illinois # 6 (I6), a high volatile bituminous coal, Utah Skyline (US), a western bituminous coal, and Black Thunder (BT), a subbituminous coal), two O2 mole fractions (O2 = 24, 36 vol %, H2 O = 14 vol %, balance CO2 ), and six size bins (53-63 µm, 63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm, and 125-150 µm). These QOIs were measured by Hecht [34] as detailed in section 5.6. 5.5 Construction of input and uncertainty map The I/U map is a list of parameters and their associated uncertainty ranges that are being considered in the VUQ analysis. An overall priority is assigned to each parameter that signifies its relative importance to the QOIs. The I/U map for the RPBL model is presented in Table 5.1. A detailed explanation of all the parameters in Table 5.1 and how their ranges were obtained is given in Dı́az-Ibarra et al [15]. In a previous sensitivity study [15], we identified five sensitive parameters for further 113 analysis: particle diameter (d p ), void fraction (φinitial ), carbon fraction of the particle (Yc ), particle emissivity (e p ), and the nondimensional thickness of the boundary layer ( 5.6 rin f rp ). Experimental data collection and uncertainty analysis 5.6.1 Reactor description The experimental data used in this analysis were collected in an optically accessible, laminar, entrained flow reactor at Sandia National Laboratories [61, 88]. The Sandia reactor is presented in Figure 5.2. This reactor was composed of a diffusion-flamelet-based Hencken burner, a quartz chimney, a particle feeding system, a collection probe, and an optical instrument for simultaneous in situ measurement of size, velocity, and temperature of a single particle. The optical instrument was fixed; the reactor moved up and down relative to the instrument to vary the position where the particle data were measured. For the particle feed system, a group of particles is fluidized with carrier gas in a test tube. A single fluidized particle then travels through a capillary tube that connects the particle feed system to the reactor. The particle is injected into the flame at position 0 (see Figure 5.2). As the particle travels along the length of the reactor, it is surrounded by hot gases from the Hencken burner. The optical instrument captures this moving, incandescent particle and LabView postprocesses the transmitted light. Using the signal recorded from the particle, the LabView software computes particle temperature with a Wien's law ratio of two signals [61]. The particle size is obtained by capturing a signal onto a coded aperture [88]. Velocity is measured using the time difference between two signals in the coded aperture [88]. A collection probe is used to sample particles to perform offline analyses such a fixed carbon or/and specific surface area. The optical instrument is not operated when the collection probe is in use. 5.6.2 Data collection Hecht [34] produced I6, US, and BT chars from their parent coals in a high temperature drop tube furnace with an environment of N2 (1473 K, 0.9 s residence time). The char particles were sieved into six narrow size bins of 53-63 µm, 63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm, and 125-150 µm. Hecht tested the sieved char samples in 12 different 114 environments in the entrained flow reactor. In each environment for each size bin, Hecht measured the size, velocity, and temperature of approximately 100 particles at 3-7 positions in the reactor. For each environment, Hecht also measured the gas temperature along the center line of the reactor without particles present with a type R thermocouple. For this VUQ analysis, we focused on the set of measurements from 36 different experiments comprising all six size ranges for the three chars in two environments, O2 vol % = 24 and 36, H2 O vol % = 14, and the balance CO2 . An experiment is defined as one char, one environment, and one particle size at all measurement positions in the reactor. The QOIs were the particle temperatures and velocities that were collected at various positions along the length of the reactor. 5.6.3 Measurement uncertainty Experimental error can be characterized as either sampling (random) error or systematic (bias) error [63]. As shown in equation (5.1), the difference between the true value and the experimentally measured value is the sum of the sampling (β) and systematic (δ) errors. Sampling error arises from uncontrolled conditions during the experiment and produces a distribution of QOIs (see Figure 5.3). We estimate this distribution from the finite number of data points that are collected. Systematic error results from problems with the measurement instrument or with the way the experimentalist uses the instrument. For this analysis, we assumed that the systematic error was much smaller than the sampling error. Hence, the total error in the system was represented by the sampling error. ȳ − y T = δ + β (5.1) At each measurement position in each of the 36 experiments, we estimated the sampling error with equation (5.2) using the temperature, velocity, and particle diameter data that were collected. In this equation, tα/2,v is a factor computed with the t-distribution and a given confidence interval (we used 95 %), s is the standard deviation, and n is the number of data points ( ∼ 100). s | β| ≤ tα/2,v √ n (5.2) Figure 5.3 shows plots of the error bars for one experiment with US char in a O2 = 24 vol %, H2 O = 14 vol %, balance CO2 environment for the 75-90 µm size bin. For all three plots, the 115 x axis is the position in the reactor above the injection point where the measurement was taken. From top to bottom, the y axis is particle temperature, velocity, and diameter. We used the particle temperature and velocity data with their associated error bars to compare with the RPBL model temperature and velocity data. We used the particle diameter as a prior range in our analysis (see section 5.9). 5.7 Simulation data collection The RPBL model solves a set of differential equations to compute species mass fractions, gas temperature, particle temperature, and particle velocity for a particle and its boundary layer moving through a specified reactor. As a particle travels the length of the reactor, it experiences different bulk gas conditions. The RPBL model uses the bulk gas conditions as boundary conditions for the species and temperature equations. For this reason, we required profiles of gas velocity, temperature, and composition in the Sandia entrained flow reactor [34]. 5.7.1 Reactor model The gas velocity and mass fraction profiles were not measured during the experiment, so the only way to estimate the data is to perform a simulation of the reactor given the experimental operating conditions. We used centerline temperature profiles measured with thermocouples to validate the results of the reactor simulation. We performed simulations of the reactor for the two environments under study using the Arches component of the Uintah software suite [48]. Uintah is a framework for solving partial differential equations on structured grids using hundreds to thousands of processors. Arches is a large eddy simulation (LES) research code developed by University of Utah researchers. It is a low-Mach, pressure projection-based, variable density code solving filtered (using the standard large eddy simulation approach) conservation equations for mass, momentum, and energy of the gas and solid phases [66-68, 71, 82]. For the simulations discussed in this paper, the gas-phase reactions were modeled using an equilibrium approach. The discrete ordinates method with S8 quadrature, representing 80 discrete directions, was used to compute radiation [59]. The simulations did not have particles. For the reactor wall, a quartz chimney, we used an emissivity of 0.93 and a 116 wall thermal conductivity that was a function of temperature [78]. We calculated the heat exchange between the reactor wall and the room where the reactor was located using a heat transfer coefficient of 4 W m2 K [43] and a room temperature of 298 K. The reactor had a cross-sectional area of 0.0508 m x 0.0508 m and a height of 0.4572 m. The computational mesh was ∼ 4 million cells. We ran the simulations on 360 cores using computational resources at Lawrence Livermore National Laboratories. Steady state was reached after 2 hours of run time. The gas temperature and gas velocity fields in the reactor for one environment (O2 = 24 vol %, H2 O = 14 vol %, balance CO2 ) are presented in Figure 5.4. This figure represents a slice through the middle of the computational domain in the y plane. The flow is clearly laminar. The carrier gas enters the reactor at position 0 on the x axis. The carrier gas velocity is higher than that of the surrounding gas while its temperature is lower. The wall has a large impact on the gas temperature as noted by the gradients of temperature near the wall. To validate this reactor model, we used center line gas temperature measurements. We extracted two profiles of gas temperature, one representing single cell values on the center line from the injection point to the exit of the reactor and the other representing an average of cells around and including the center line. Both profiles are presented in Figure 5.5. In the single cell gas profile, the initial temperature is 800 K and then rises rapidly; the averaged profile is more homogeneous. However, the effect of the injection gas is evident from the slight dip at 0.01 s. We corrected the center line gas temperatures extracted from the reactor model (Figure 5.5) to thermocouple (Ttc ) measurements using equation (5.3), which is an energy balance on the thermocouple. In this equation, the gas temperature (Tg ) was extracted from the reactor model, the thermal conductivity of the gas (k gas ) was computed from simulation outputs (gas composition and temperature), the Nusselt number was computed with a correlation [34] (see equation (5.4)), etc = 0.2 is the surface emissivity of the thermocouple, and Tw = 500 K is the reactor wall temperature [34]. In equation (5.4), the Reynolds number (Re) was computed from the diameter of the thermocouple (dtc = 2.54e − 5 m) and the velocity of the gas (extracted from the simulation output). The Prandtl number (Pr) was also computed from simulation outputs (gas composition and temperature). 117 Ttc = Tg − σetc dtc ( Tg4 − Tw4 ) Nuk gas Nu = 0.989Re0.33 Pr0.33 (5.3) (5.4) The results of this correction are shown in Figure 5.6. The center profile Ttc values are close to the experimental data between positions 0 m and 0.0254 m of the reactor. Both profiles are close to the experimental data for positions greater than 0.0635 m. Because the char particles can take different paths than just the center line, the spatially-averaged gas temperature profile is more representative of the environment a particle may experience as it accounts for the center line and eight other paths. Therefore, we used the average profile in the RPBL simulations. The initial particle temperature for all the RPBL simulations was 800 K, which is the gas temperature at the center line of position 0 m. With this initial temperature, we matched most of the thermocouple measurements near the injection point. We also validated the reactor model by comparing the spatially-averaged gas velocity profile from the reactor model with the velocities of the smallest US char particles (53 - 63 µm) in the experiment. The velocities of these small particles are a reasonable surrogate for the gas velocity because they rapidly approach the terminal (gas) velocity. We had to reduce the average gas velocity by 3% to be in the uncertainty range of the small particle velocity, as shown in Figure 5.7. One reason that we could not obtain a closer match between the two velocity profiles was that the burner geometry was not resolved in the simulation. The Hencken burner consists of a set of small tubes, each of which holds a small flame. We simulated the burner as a single large flame. Nevertheless, our correction is close to the particle velocity measurement precision of ± 2% reported by Murphy and Shaddix [61]. The spatially-averaged gas composition (averaged over the same nine points surrounding and including the center line) was also extracted from the reactor model, as shown in Figure 5.8. For the initial composition, we assumed that the gas inside the particle and in its boundary layer was 99% CO2 because that was the carrier gas. We used the average gas temperature, velocity, and mass fraction profiles from the reactor model as boundary conditions in the RPBL model at rin f , representing the bulk gas condition. 118 5.7.2 Data collection With the bulk gas profiles obtained from the reactor simulations of the two environments, our next step was to run the RPBL model based on the five parameters in our study and the requirements of the surrogate model. Given these five parameters and their ranges as presented in Table 5.1, we used the Uncertainty Quantification Toolkit (UQTk) [13] to produce a design of experiments. For the O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment, a polynomial chaos (PC) surrogate model of order 4 with a full quadrature rule required a total of 3125 cases to be run. For the O2 = 36 vol%, H2 O = 14 vol%, balance CO2 environment, a PC surrogate model of order 5 with a full quadrature rule required a total of 7776 cases to be run. We ran these RPBL model cases on a local Linux cluster at the University of Utah. We set the total number of cells to 60 and the number of cells inside the particle to 30, which meant we had to solve 781 ODEs. With 520 cores, it was possible to run all simulations required for the VUQ analysis in one day. To compare the QOIs from the RPBL model results with those from the experimental data, we defined a procedure for computing particle velocity and temperature. Because the velocity of the particle is resolved by the RPBL model, we only needed to extract the velocity at the positions where data were taken. For particle temperature, we computed the average of the 30 cells in the particle. Then, we extracted the average particle value at the positions where data were taken. 5.8 Construction of surrogate models We created PC surrogate models [13] for the consistency analysis with the RPBL data. We required a surrogate model because it took more than 1 hour of computational time to obtain output data (particle temperature and velocity) from the RPBL model for one set of parameters. This function evaluation time is too long for the consistency analysis tool, so we precomputed the particle temperature and velocity with the RPBL model, then built a surrogate model with these data. For each experiment, the number of measurement positions varied from 3 to 7. For example, in Figure 5.3, experimental data were taken at five positions. Hence, for this experiment, we built five PC surrogate models for particle temperature and five for ve- 119 locity. We used the experimentally-measured particle diameters as the prior range in the consistency analysis described in section 5.9. We used a PC surrogate model of order 4 for the O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment. In Figure 5.9, we present verification plots for ten PC surrogate models, five corresponding to the particle temperature (first column) and five to the particle velocity (second column) for the experiment presented in Figure 5.3. The rows are the position in the reactor. The x axis is the surrogate output for the variable and the y axis is the RPBL simulation output for the same variable. The surrogate model is a good representation of the simulation data because all the x, y pairs are close to y = x. Thus, the PC surrogate models reproduce the output values from the RPBL simulations. We performed one additional verification to determine if the PC surrogate models exhibited the behavior we expected: we evaluated the PC surrogate models for one experiment with a sample of 10,000 points. We held three parameters constant (Yc,intial =0.75, rin f rp = 85, e p = 0.55) and varied d p and φinitial using the ranges presented in Table 5.1. We present the results in Figure 5.10. The particle temperature PC surrogate models are in the first column and the particle velocity surrogate models in the second. Each row is a location in the reactor where data were taken for the 75-90 µm size bin. The particle temperature plots show how d p and φinitial affect the particle temperature. In the first row (position 0.049 m), particles smaller than 40 µ m with high φinitial values have temperatures of ∼ 1650 - 1700 K, which are close to the gas temperature. The maximum temperatures (∼ 1920 K) are observed for particles between 40 to 80 µ m, indicating that these particles are reacting. Particles bigger than 80 µ m are cooler than the gas, which means that these particles are being heated by the gas. Particles with low φinitial values (heavy particles) are cooler than those with high φinitial . In the second row (position 0.0735 m), particles smaller than 40 µ m with low φinitial (heavy particles) and particles larger than 80 µ m with a high φinitial (light particles) are close to the gas temperature (∼ 1700 K), which means these particles are composed primarily of ash. Maximum temperatures, indicating reacting particles, are observed for heavy particles between 60 to 100 µ m and for light particles between 70 to 120 µ m. Particles larger than 130 µ m with low φinitial are still heating up. In the third, fourth, and fifth rows, the particles show similar behavior to that in the first and second rows albeit with a different size of particle. The particle 120 temperature plots exhibit the behavior that we expect; small particles react earlier in the reactor than large particles. From the particle velocity plots, we can see that small particles have larger velocities than the large particles. For particles bigger than 60 µ m, the velocity decreases with decreasing φinitial . This result makes sense because the particle weight increases with decreasing φinitial . We plotted the PC surrogate models at different positions and with varying values for the other three parameters and observed similar behavior. We also checked the PC surrogate models for the O2 = 36 vol%, H2 O = 14 vol%, balance CO2 environment. In this case, we used a PC surrogate model of order 5. We concluded that since the PC surrogate models showed behavior that was similar to RBPL model results, these PC surrogate models were adequate for the consistency analysis. 5.9 Analysis of model outputs We employed a consistency measure analysis referred to as bound-to-bound consistency [75]. In this analysis, the model outputs and experimental data are compared using equation (5.5). ∆= |ym,e ( x ) − ye | ue (5.5) Here, ym,e ( x ) is the model data point defined by the set of x parameters, ye is the experimental data point, and ue is the experimental uncertainty. If ∆ ≤ 1, the model data point using parameter set x is consistent with the experimental data point. If the model outputs for a parameter set x are consistent with all the experimental measurements, x is a consistent point in our analysis. We carried out a consistency analysis for each of the 36 experiments performed by Hecht [34] as described in section 5.6. Here, we describe the consistency analysis for US char in the O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment and the 75-90 µm size bin, but all the analyses were similar. We first computed ue and ye for the particle temperature, velocity, and diameter. Here, ue is the sampling error (β) and ye is the mean experimental value computed for each position as discussed in section 5.6.3. The particle temperature and velocity are the QOIs that we compared with the PC surrogate models. The d p uncertainty range is used as a 121 prior estimate for the diameter of the particle in the next step. Second, we created a random sample of 500,000 points within the parameter space defined by the ranges of φinitial , Yc , e p , and rin f rp presented in Table 5.1. For the d p range, we chose the maximum upper bound and the minimum lower bound computed from the experimental data. Third, we built the PC surrogate models for the particle temperature and velocity measurements [13] and then evaluated each surrogate model at all 500,000 points. Fourth, we applied the consistency analysis to the particle temperature and velocity measurements to compute ∆ (see equation (5.5)). For each point (set of parameter values) in the 500,000 point set, we extracted the maximum ∆, ∆max , from the ∆ values at each measurement location. From this vector of ∆max values, we extracted the minimum ∆max , which is the most consistent point in the set of 500,000 points. If this minimum ∆max < 1, then we have found at least one consistent point. If the minimum ∆max > 1, we need to improve our models, explore other parameters, or review ue . Fifth, we reviewed the ∆max vector and selected all points for which ∆max < 1. From this group of points, we obtained new ranges for the five model parameters by finding the global maximum and minimum values of all the consistent points. The last step was a refined consistency analysis with the new consistent ranges. We created a random sample of 10,000 points within the consistent ranges and then reevaluated the PC surrogate models at these 10,000 points. We obtained a set of consistent points where ∆max < 1; these points are presented in Figure 5.11 for one experiment. In this figure, the ranges for the three axes are the original ranges presented in Table 5.1 and the box that contains the points is the consistent range. The color of each point corresponds to ∆max . Table 5.2 presents the prior and the consistent ranges for this particular consistency analysis. We compare the experimental range (red error bars), simulation (model) range (blue error bars), and consistency range (green error bars) for particle temperature and velocity in Figure 5.12. In both plots, the consistent range is equal to or smaller than the experimental uncertainty range. In addition, the model can reproduce the scatter of particle temperature in all positions. We hypothesize that the large uncertainty in the particle temperatures is due to the large range in particle diameter in a given size bin and the nonconstant physical 122 and chemical properties of a single particle. We are able to span the range of experimental particle velocities because we added a momentum equation for the particle rather than assuming the same velocity for gas and particle. We performed a similar consistency analysis for the other experiments and obtained consistency with all 36 experiments. The consistent ranges of all five parameters are shown in Table 5.3 for the O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment and in Table 5.4 for the O2 = 36 vol%, H2 O = 14 vol%, balance CO2 environment. There are two additional columns for d p . The nominal column corresponds to the size bins reported by Hecht [34] (obtained from sieving data). The prior column corresponds to the prior ranges used in the model. The prior range was taken from particle sizes measured in the entrained flow reactor for a given sieving size bin. It is clear from the table data that the particle size variability measured in the entrained flow reactor was much greater than the size range of the unburned char particles based on sieving. For all chars and all size bins, the consistent range for d p was reduced. With the exception of the 53-63 µ m and 63-75 µ m size bins for BT char in the O2 = 36 vol% environment (see Table 5.4), the nominal range was contained within the consistent range. For these two exceptions, measurements were taken close to position 0 where the effects of the injection point are significant. 5.9.1 Effect of particle size The results of six consistency analyses (representing the six bin sizes) for I6 char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment are presented in Figure 5.13 and Figure 5.14. Each row corresponds to a bin size. The first column is particle temperature and the second column is particle velocity. The x axis is the position and the y axis is particle temperature or particle velocity. In each position, the experimental error bars and consistency ranges are slightly offset with respect to the simulation range for visualization purposes. The actual position of the measurement corresponds to the position of the simulation range bar. From both figures, we can see that the RPBL model outputs and the experimental data are consistent for all the size bins. Thus, we can conclude that the model represents the data over the range of size bins for I6 char. We obtained a similar result for US char. For BT char, we obtained consistency for the six size bins, but the 123 consistency ranges were reduced, especially for the e p parameter (see Table 5.3). Comparing the plots in both figures, we can see that the measurements were taken at different positions depending on the size bin. For small size bins, particles were measured close to the injection position, while for larger size bins, particles were measured farther from the injection point. Small particles heat up faster than big particles and have a smaller amount of carbon than a bigger particle. Therefore, the char burnout time of a small particle is shorter than that of a big particle. All of these effects are captured by the model. 5.9.2 Effect of char type We compare consistency analysis results for US, I6, and BT chars in Figure 5.15. This figure shows error plots for the 90-106 µm size bin in the O2 = 24 vol%, H2 O= 24 vol%, balance CO2 environment. The first column is particle temperature and the second column is particle velocity; each row corresponds to a char type. Both I6 and US chars have similar temperature ranges while BT char has a higher particle temperature. The particle velocities are similar for the I6 and US chars but lower for BT chars at the first two positions. The consistent and experimental data ranges for particle temperature nearly overlap in all cases, indicating that the model is reproducing the range of the particle temperatures for all three chars. For BT char, we found consistency with all QOIs for all the size bins. However, the e p parameter was reduced toward the lower bound of the prior range. This trend can be explained by considering the particle temperatures in Figure 5.15; those for BT char are higher than those for US and I6 chars, so to find consistency, lower values of e p (i.e. reducing the amount of energy lost to the reactor wall) were required. 5.9.3 Effect of O2 concentration As part of our consistency analysis, we analyzed results from two O2 mole fraction environments, O2 = 24 and 36 vol% [34]. In Figure 5.16, we compare error bar plots in the two environments for the 53-63 µm size bin of a US char. The model accurately reproduces the measured particle temperatures at all positions for both environments. Also, higher temperatures are obtained with the O2 = 36 vol% environment than with the O2 = 24 vol% environment, which makes sense because the amount of O2 available for reaction is increased. We obtained similar results for all size bins with US, I6, and BT chars. 124 5.10 Feedback and feed-forward In this step, we review our analysis, check assumptions, and make recommendations for the next VUQ cycle. We explored five parameters (d p , φinitial , Yc,intial , rin f rp , e p ) using 36 experiments. The QOIs for our consistency analysis were particle temperature and particle velocity measured at various positions in the Sandia entrained flow reactor. With the RPBL model, we obtained simultaneous consistency across QOIs taken at multiple locations in 36 different experiments performing 36 consistency analyses, but we could not reduce the uncertainty of the five parameters. We believe that this scatter in the temperature and velocity data is due to the distribution of particle sizes and physical properties (especially φinitial and Yc,intial ) in the heterogeneous char samples. 5.10.1 Particle size and shape The most important parameter for obtaining consistency was d p . This parameter determines where the particle is located in the reactor and the amount of total carbon in the particle. The char used in this analysis was sieved into six different size ranges. However, no additional characterization of the particle distribution was done and only the nominal values were reported [34]. We used a prior range of particle sizes broader than these nominal sizes. After consistency analysis, the ranges were reduced for all size bins but were still much larger than the nominal ranges (see Table 5.3 and Table 5.4). We propose four possible explanations for the discrepancy. First, d p needs to be better characterized experimentally. We recommend that particle size distribution measurements be performed after the sieving. Second, the RPBL model assumes a constant particle diameter during the reaction process. Since we found consistency with all the experimental data, we think that this assumption is adequate. However, the particles may be increasing in size due to the reaction process. Third, the particles, assumed to be spherical in the RBPL model, may be nonspherical. The RPBL model could be improved by resolving a cylindrical particle in addition to a spherical particle. Fourth, d p is a distribution, not a range, and there may be bias in the experimental measurements such that the distribution is not sampled uniformly at each measurement location. With the model, we could use a distribution of d p rather than a range of d p as the prior. To obtain the particle phase size distribution, the Sandia reactor could be simulated using Arches [66-68] with the direct 125 quadrature method of the moments (DQMOM) [25]. 5.10.2 Other model parameters The parameters related to the particle morphology are φinitial and Yc,intial . These two parameters are also important for computing the particle temperature and velocity. The consistent range of φinitial was only reduced for the 53-63 µ m size bin experiment with BT char in the O2 = 36 vol% environment. However, looking at Figure 5.11, we can see that consistent range of d p changes depending on φinitial . Thus, if we could estimate this value with greater accuracy, we could find a smaller range for d p . The Yc,intial parameter has not been reduced in most of the experiments. While we use a range of values for φinitial and Yc,intial in this study, these parameters should be estimated from a devolatilization model prior to the onset of oxidation. 5.10.3 Model assumptions We assumed a gradient of pressure ∼ 0 in the RPBL model. There were no data available to test this assumption, but the model reproduced particle temperature and velocity behavior when the O2 mole fraction was increased or when the particle diameter was varied, so this assumption did not negatively impact RPBL model outputs. We used the same heterogeneous reaction mechanism for all three chars. We found that the US and I6 chars were well represented by this reaction mechanism. However, for the BT char, this reaction mechanism may not be adequate. We base our conclusion on the e p consistency range for BT char in Table 5.3 and Table 5.4. The low values of e p result in higher particle temperatures because the particle is losing less energy to the reactor wall. Particle temperature depends on the consumption of carbon, which is determined by the reaction mechanism. If we want to have one heterogeneous reaction mechanism for all types of char, we must explore parameters related to the reaction mechanism. 5.10.4 Uses for RPBL model Because RPBL simulation times are on the order of hours, it cannot be implemented in an LES code such as Arches [66-68]. However, this model can produce data for development and validation of simpler models. RPBL can also be used to produce surrogate or precomputed models for implementation in LES applications. 126 Given that we obtained consistency with the RPBL model and Sandia char oxidation data for three types of char, we might ask if this model can be used to perform predictions with another type of char or with another type of reactor. In order to use RPBL as a predictive tool, we need to eliminate as much of the bias (error) as possible from both the model and the experimental data. We did remove some bias in the SKIPPY model by adding new physics to create the RPBL model. However, some results indicate that other forms of bias must still be corrected. For example, the reaction mechanism reproduced the velocity and temperature data for I6 and US char particles, but for the BT char, the e p parameter was reduced in order to find consistency. Thus, we think that the model cannot be used to predict behavior from other char types until a heterogeneous reaction mechanism is developed that differentiates among types of chars. Also, we did not look for ways to identify and remove bias from the experimental measurements. We computed the experimental error assuming that systematic error was small in comparison with the sampling error. Only by removing experimental bias will we know if our model has the necessary physics to predict char behavior in other systems. 5.11 Conclusion We followed a six-step methodology to perform a VUQ analysis with the RPBL model output and with experimental char oxidation data collected by Hecht [34]. The QOIs in this analysis were the particle temperature and velocity. We selected five RPBL model parameters that showed higher sensitivity to the QOIs for our analysis: d p , φinitial , Yc , e p , and rin f rp . The uncertainty in the experimental data was computed using the approach of Oberkampf and Roy [63] from the ∼ 100 measurements made at each position in each experiment. We performed consistency analyses between experimental particle temperatures and velocities and those computed with the RPBL model for 36 experiments that corresponded to three types of char (US, I6, and BT), two gas compositions (O2 vol% = 24, 36; H2 O vol% = 14; balance with CO2 ), and six size bins (53-63 µm, 63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm, and 125-150 µm). Consistency was found with the 36 experimental data sets. From the consistency analysis, we learned that the US and I6 chars exhibited similar behavior with respect to particle temperature and velocity. However, BT char produced 127 higher temperatures than either US or I6 in both O2 environments that were analyzed. When the O2 mole fraction was increased, the experimentally-measured particle temperatures also increased; this behavior was captured by the RBPL model as was the effect of varying the particle size for all three chars. Our objectives with this analysis were to determine if the RPBL model outputs could reproduce the experimental data collected by Hecht [34] and if we could use the experimental data to estimate RPBL model parameters. From the analysis, we learned that the RPBL model does reproduce the experimental data with three char types in two oxygen environments. However, the most of the ranges of the five parameters we studied were not reduced in the consistency analysis. To answer the question of whether the RPBL model is validated and can thus be used as a predictive tool, we need to consider whether the bias has been removed from both the modeling the experimental data. In this analysis, we used the heterogeneous reaction mechanism from the original SKIPPY model. We think that model bias cannot be fully examined until we apply a reaction mechanism from a different source. On the experimental side, we did not account for any sources of error related to the experimental device. These sources of bias need to be removed and and new VUQ studies performed before the RPBL model can be used to perform predictions in other conditions with some known degree of accuracy. 5.12 Acknowledgements This material is based upon work supported by the Department of Energy, National Nuclear Security Administration, under Award Number DE-NA0002375. Sandia National Laboratories is a multi-program laboratory operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC0494AL85000. 128 Table 5.1. Input/uncertainty map for RPBL Parameter Range min max Model Parameters rin f 50 120 r p [−] e p [−] 0.1 1 Scenario parameters d p [µm] 30 160 φinitial 0.15 0. 7 Yc,intial [−] 0.5 1 Table 5.2. Parameter ranges prior to and after consistency analysis for US char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment and 75-90 µm size bin Parameter Prior range Consistent Nominal range value d p [µm] 36.9- 146.5 42.0-146.0 75-90 φinitial 0.15 - 0.7 0.15 - 0.7 Yc,intial [−] 0.5-1 0.5-1 rin f 50-120 50-120 r p [−] e p [−] 0.1-1 0.1-1 - 129 Table 5.3. Consistent parameter ranges for US, I6, and BT chars in an O2 = 24 vol%, H2 O = r 14 vol%, balance CO2 environment. Prior ranges are rinpf (50-120), φinitial (0.15 - 0.7), Yc,intial (0.5-1), and e p (0.1-1); see Table 5.1. rin f d p [µm] φinitial Yc,intial [−] e p [−] r p [−] Nominal 53-63 63-75 75-90 90-106 106-125 125-150 Prior 22-112 27-119 37-147 36-140 53-164 66-199 Nominal 53-63 63-75 75-90 90-106 106-125 125-150 Prior 28-118 30-119 29-150 30-130 67-144 73-165 Nominal 53-63 63-75 75-90 90-106 106-125 125-150 Prior 17-124 27-138 21-141 33-145 36-165 58-171 Utah Skyline (US) Consistent Consistent 36-112 50 -120 0.15-0.70 0.5 -1.00 34-120 50-120 0.15-0.70 0.5-1.00 42-142 50-120 0.15-0.70 0.5-1.00 38-148 50-120 0.15-0.70 0.5-1.00 70-160 50-120 0.15-0.70 0.5-1.00 74-160 50-120 0.15-0.70 0.5-1.00 Illinois # 6 (I6) Consistent Consistent 34-120 50 -120 0.15-0.70 0.5 -1.00 44-120 50-120 0.15-0.70 0.5-1.00 30-150 50-120 0.15-0.70 0.5-1.00 32-132 50-120 0.15-0.70 0.5-1.00 74-146 50-120 0.15-0.70 0.5-1.00 82-160 50-120 0.15-0.70 0.5-1.00 Black Thunder (BT) Consistent Consistent 60-126 50 -120 0.15-0.70 0.5 -1.00 70-120 50-115.92 0.16-0.70 0.64-1.00 70-144 50-120 0.15-0.70 0.5-1.00 88-148 50-120 0.15-0.70 0.5-1.00 102-160 50-120 0.15-0.70 0.5-1.00 96-160 50-120 0.15-0.70 0.5-1.00 0.1 -1 0.1-1 0.1-1 0.1-1 0.1-1 0.1-1 0.1 -1 0.1-1 0.1-1 0.1-1 0.1-1 0.1-1 0.1 -0.69 0.1-0.36 0.1-0.46 0.1-0.76 0.1-0.79 0.1-1 130 Table 5.4. Parameter ranges after consistency analysis for US, I6, and BT chars in an O2 r = 36 vol%, H2 O = 14 vol%, balance CO2 environment. Prior ranges are rinpf (50-120), φinitial (0.15 - 0.7), Yc,intial (0.5-1), and e p (0.1-1); see Table 5.1. rin f d p [µm] φinitial Yc,intial [−] e p [−] r p [−] Nominal 53-63 63-75 75-90 90-106 106-125 125-150 Simulation 25-127 32-126 41-149 37-159 40-192 43-203 Nominal 53-63 63-75 75-90 90-106 106-125 125-150 Simulation 31-117 22-129 39-163 32-167 47-196 42-189 Nominal 53-63 63-75 75-90 90-106 106-125 125-150 Simulation 23-138 24-134 27-155 28-172 32-179 43-176 Utah Skyline (US) Consistent Consistent 44-102 50 -120 0.15-0.70 0.5 -1.00 52-128 50-120 0.15-0.70 0.5-1.00 68-152 50-120 0.15-0.70 0.5-1.00 70-160 50-120 0.15-0.70 0.5-1.00 80-160 50-120 0.15-0.70 0.5-1.00 92-158 50-120 0.15-0.70 0.5-0.750 Illinois # 6 (I6) Consistent Consistent 46-118 50 -120 0.15-0.70 0.5 -1.00 54-130 50-120 0.15-0.70 0.5-1.00 64-160 50-120 0.15-0.70 0.5-1.00 72-160 50-120 0.15-0.70 0.5-1.00 86-160 50-120 0.15-0.70 0.5-1.00 84-160 50-120 0.15-0.70 0.5-1.00 Black Thunder (BT) Consistent Consistent 76-88 56-82 0.44 -0.67 0.91-1.00 76-136 50-120. 0.15-0.70 0.5-1.00 82-156 50-120 0.15-0.70 0.5-1.00 98-160 50-120 0.15-0.70 0.5-1.00 100-160 50-120 0.15-0.70 0.5-1.00 98-160 50-120 0.15-0.70 0.5-1.00 0.1 -1 0.1-1 0.1-1 0.1-1 0.1-1.0 0.1-0.55 0.1 -1 0.1-1 0.1-1 0.1-1 0.1-0.9 0.1-0.85 0.12-0.21 0.1-0.64 0.1-0.63 0.1-0.68 0.1-0.5 0.1-0.73 131 Selection of quantities of interest (QOIs) Step 1 Construction of input and uncertainty (I/U) map Step 2 Collection of experimental and simulation data Step 3 Construction of surrogate models Step 4 Analysis of model outputs - consistency analysis Step 5 Feedback and feedforward Step 6 Figure 5.1. Six-step methodology with consistency analysis Position 0 Char Feed Figure 5.2. Schematic of optically accessible, entrained flow reactor at Sandia National Laboratory. Figure adapted with permission from [34] 2100 2000 1900 1800 1700 1600 1500 Diameter[µ m] Particle velocity [m/s] Particle temperature [K] 132 3.2 3.0 2.8 2.6 2.4 2.2 2.0 1.8 350 300 250 200 150 100 50 0 CO2-O2-24-H2O-14-coal-US-size-75-90um 0.06 0.08 0.10 0.12 Position [m] 0.14 0.06 0.08 0.10 0.12 Position [m] 0.14 0.06 0.08 0.10 0.12 Position [m] 0.14 Figure 5.3. Average temperature, velocity, and diameter with error bars for particles in the 75-90 µm size bin. The black dots are the experimental points, and the red error bars are computed from equation (5.1) with a confidence interval of 95 %. These data correspond to US char in an environment of O2 = 24 vol %, H2 O = 14 vol %, and balance CO2 . 133 Gas temperature Gas velocity Figure 5.4. Pseudocolor plots of gas temperature and gas velocity in the Sandia entrained flow reactor for O2 = 24 vol %, H2 O = 14 vol % , balance CO2 environment. 1800 Averaged gas Center gas Temperature [K] 1600 1400 1200 1000 800 0.00 0.05 0.10 Position [m] 0.15 0.20 Figure 5.5. Gas temperature profiles (single cell and averaged over 9 cells) along the center line of the entrained flow reactor. 134 1700 Temperature [K] 1600 1500 1400 1300 Averaged themocouple Center themocouple Experimental thermcouple 1200 1100 0.00 0.05 0.10 Position [m] 0.15 0.20 Figure 5.6. Thermocouple temperature profiles along the center line of the entrained flow reactor. The profile is from position 0 [m] to position 0.254 [m]. 3.2 3.1 Averaged gas velocity Particles of 53 to 63µm Velocity [m/s] 3.0 2.9 2.8 2.7 2.6 2.5 0.00 0.05 0.10 0.15 Position [m] 0.20 0.25 0.30 Figure 5.7. Comparison of spatially-averaged gas velocity profile from reactor simulation and experimental particle velocity profile for particles in the 53 - 63 µm size range. The average gas velocity was reduced by 3% to be in the uncertainty range of the small particle velocity. 135 0.8 0.7 Mass fraction [-] 0.6 0.5 O2 CO2 H2 O 0.4 0.3 0.2 0.1 0.0 0.00 0.05 0.10 Position [m] 0.15 0.20 Figure 5.8. Spatially-averaged mass fraction profiles for O2 , H2 O, and CO2 along the center line of the simulated entrained flow reactor. 136 2200 1800 Simulation Simulation 2000 Particle temperature [K] at position 0.049 [m] y=x 1600 1400 2000 2200 2000 Particle temperature [K] at position 0.098 [m] y=x 1500 1600 1700 1800 Surrogate model 1900 2000 2100 2100 3.05 3.00 2.95 2.90 2.85 2.80 2.75 2.70 2.65 2.65 2100 3.05 3.00 2.95 2.90 2.85 2.80 2.75 2.70 2.70 Particle temperature [K] at position 0.1225 [m] y=x 1900 1800 1700 1600 1500 1500 1600 2100 Particle temperature [K] at position 0.147 [m] y=x 1800 1900 Surrogate model 2000 1900 Simulation Simulation 2000 1700 1800 1700 1600 1500 1500 1600 1700 1800 1900 Surrogate model 2000 2.95 2.90 2.85 2.80 2.75 2.70 2.65 2.60 2.55 2.55 3.00 2.95 2.90 2.85 2.80 2.75 2.70 2.65 2.60 2.60 Simulation 2100 2000 1900 1800 1700 1600 1500 1400 1400 Simulation 1600 1800 Surrogate model Particle temperature [K] at position 0.0735 [m] 2200 y=x 2100 2000 1900 1800 1700 1600 1500 1400 1300 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 Surrogate model 2100 Simulation 1400 Simulation Simulation Simulation 1200 1200 2.85 2.80 2.75 2.70 2.65 2.60 2.55 2.50 2.50 Particle velocity [m/s] at position 0.049 [m] y=x 2.55 2.60 2.65 2.70 Surrogate model 2.75 2.80 2.85 Particle velocity [m/s] at position 0.0735 [m] y=x 2.60 2.65 2.70 2.75 2.80 Surrogate model 2.85 2.75 2.80 2.85 Surrogate model 2.90 2.80 2.85 2.90 Surrogate model 2.95 2.90 2.95 Particle velocity [m/s] at position 0.098 [m] y=x 2.65 2.70 2.95 3.00 Particle velocity [m/s] at position 0.1225 [m] y=x 2.70 2.75 3.00 3.05 Particle velocity [m/s] at position 0.147 [m] y=x 2.75 2.80 2.85 2.90 Surrogate model 2.95 3.00 3.05 Figure 5.9. Verification plots for PC surrogate models. The red line is y = x. The positions correspond to an experiment with US char in an O2 = 24 vol%, H2 O = 14 vol% balance CO2 environment and bin size of 75-90 µm. 137 0.6 0.3 0.2 40 80 100 dp µm 120 140 0.3 0.2 40 60 80 100 dp µm 120 140 160 Particle temperature [K] at position 0.098 [m] 0.6 φinitial 0.5 0.4 0.3 0.2 40 60 80 100 dp µm 120 140 160 0.7 Particle temperature [K] at position 0.1225 [m] 0.6 φinitial 0.5 0.4 0.3 0.2 40 60 80 100 dp µm 120 140 160 Particle temperature [K] at position 0.147 [m] 0.6 0.5 φinitial 0.2 40 0.4 0.3 0.2 40 60 80 100 dp µm 120 140 160 1920 1860 1800 1740 1680 1620 1560 1500 1440 1950 1900 1850 1800 1750 1700 1650 1600 1950 1900 1850 1800 1750 1700 1650 1600 1950 1900 1850 1800 1750 1700 1650 1600 0.7 60 80 100 dp µm 120 140 160 Particle velocity [m/s] at position 0.0735 [m] 0.6 0.5 φinitial φinitial 0.5 0.7 0.3 Particle temperature [K] at position 0.0735 [m] 0.4 0.7 0.5 0.4 160 0.6 Particle velocity [m/s] at position 0.049 [m] 0.6 0.4 0.3 0.2 40 0.7 60 80 100 dp µm 120 140 160 Particle velocity [m/s] at position 0.098 [m] 0.5 0.4 0.3 0.2 40 60 80 100 dp µm 120 140 0.6 0.5 0.4 0.3 0.2 0.7 60 80 100 dp µm 120 140 160 Particle velocity [m/s] at position 0.147 [m] 0.6 0.5 0.4 0.3 0.2 40 60 80 100 dp µm 120 140 2.88 2.84 2.80 2.76 2.72 2.68 2.64 2.60 160 0.7 Particle velocity [m/s] at position 0.1225 [m] 40 2.79 2.76 2.73 2.70 2.67 2.64 2.61 2.58 2.55 2.96 2.92 2.88 2.84 2.80 2.76 2.72 2.68 0.6 φinitial 0.7 60 0.7 φinitial φinitial 0.5 0.4 1920 1840 1760 1680 1600 1520 1440 1360 φinitial Particle temperature [K] at position 0.049 [m] φinitial 0.7 160 3.00 2.96 2.92 2.88 2.84 2.80 2.76 2.72 3.03 3.00 2.97 2.94 2.91 2.88 2.85 2.82 2.79 2.76 Figure 5.10. Scatter plots for PC surrogate models; left column is particle temperature and right column is particle velocity. The positions (rows) correspond to an experiment with US char in an O2 = 24 vol%, H2 O = 14 vol% balance CO2 environment and bin size of r 75-90 µm. Ranges for d p and φinitial are presented in Table 5.1; Yc,intial =0.75 , rinpf = 85, and e p = 0.55. 138 Figure 5.11. Consistent sample for US char in an O2 = 24 vol%, H2 O =14 vol%, balance CO2 environment and 75-90 µm size bin. 139 QOI: Particle temperature 2200 2000 1800 1600 1400 1200 0 3.2 QOI: Particle velocity 3.1 Experimental range Simulation range Consistent range 1 2 3 4 5 6 8 9 10 11 Experimental range Simulation range Consistent range 3.0 2.9 2.8 2.7 2.6 2.5 5 6 7 Figure 5.12. Consistent range for US char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment and 75-90 µm size bin. 140 I6 53-63 µm 2000 1800 1600 0.05 0.10 2.7 2.6 2.5 1800 1600 0.05 0.10 0.15 Position [m] Experimental range Simulation range Consistent range 3.0 1400 0.20 I6 63-75 µm 3.1 Experimental range Simulation range Consistent range QOI: Particle velocity QOI: Particle temperature 2.8 2.3 0.20 I6 63-75 µm 2.9 2.8 2.7 2.6 0.05 0.10 0.15 Position [m] 2.5 0.20 I6 75-90 µm 2200 0.05 0.10 0.15 Position [m] 0.20 I6 75-90 µm 3.2 3.1 2000 QOI: Particle velocity QOI: Particle temperature 0.15 Position [m] 2000 1800 1600 1400 1200 2.9 2.4 2200 1200 Experimental range Simulation range Consistent range 3.0 1400 1200 I6 53-63 µm 3.1 Experimental range Simulation range Consistent range QOI: Particle velocity QOI: Particle temperature 2200 Experimental range Simulation range Consistent range 0.05 0.10 0.15 Position [m] 0.20 3.0 2.9 2.8 2.7 Experimental range Simulation range Consistent range 2.6 2.5 0.05 0.10 0.15 Position [m] 0.20 Figure 5.13. Consistent range for I6 char in an O2 = 24 vol%, H2 O =14 vol%, balance CO2 environment for 53-63 µm, 63-75 µm, and 75-90 µm size bins. 141 I6 90-106 µm 1800 1600 1400 Experimental range Simulation range Consistent range 0.05 0.10 0.15 Position [m] QOI: Particle velocity QOI: Particle temperature 1800 1700 1600 0.05 0.10 0.15 Position [m] 0.05 0.10 0.15 Position [m] 0.20 I6 106-125 µm Experimental range Simulation range Consistent range 3.5 3.0 2.5 1.5 0.20 I6 125-150 µm 3.3 2100 3.2 2000 QOI: Particle velocity QOI: Particle temperature 2.7 2.0 Experimental range Simulation range Consistent range 2200 1900 Experimental range Simulation range Consistent range 1600 1500 0.05 0.10 0.15 Position [m] 0.20 I6 125-150 µm Experimental range Simulation range Consistent range 3.1 3.0 2.9 2.8 2.7 1400 1300 2.8 4.0 1900 1700 2.9 4.5 2000 1800 3.0 2.5 0.20 2100 1400 3.1 2.6 I6 106-125 µm 2200 1500 Experimental range Simulation range Consistent range 3.2 2000 1200 I6 90-106 µm 3.3 QOI: Particle velocity QOI: Particle temperature 2200 0.05 0.10 0.15 Position [m] 0.20 2.6 0.05 0.10 0.15 Position [m] 0.20 Figure 5.14. Consistent range for I6 char in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment for 90-106 µm, 106-125 µm, and 125-150 µm size bins. 142 US 90-106 µm 3.1 2000 1800 1600 1400 1200 Experimental range Simulation range Consistent range 0.05 0.10 0.15 Position [m] QOI: Particle velocity 1400 Experimental range Simulation range Consistent range 0.05 2.7 2.6 0.05 0.10 0.15 Position [m] 0.10 0.15 Position [m] I6 90-106 µm Experimental range Simulation range Consistent range 3.1 3.0 2.9 2.8 2.7 2.5 0.20 3.2 2100 3.1 QOI: Particle velocity 2000 1900 1800 1700 1600 Experimental range Simulation range Consistent range 0.05 0.10 0.20 2.6 BT 90-106 µm 2200 QOI: Particle temperature 2.8 3.2 1600 1400 2.9 3.3 1800 1500 3.0 2.4 0.20 2000 1200 US 90-106 µm Experimental range Simulation range Consistent range 2.5 I6 90-106 µm 2200 QOI: Particle temperature 3.2 QOI: Particle velocity QOI: Particle temperature 2200 0.05 0.10 0.15 Position [m] 0.20 BT 90-106 µm Experimental range Simulation range Consistent range 3.0 2.9 2.8 2.7 2.6 0.15 Position [m] 0.20 2.5 0.05 0.10 0.15 Position [m] 0.20 Figure 5.15. Consistent ranges for I6, US, and BT chars in an O2 = 24 vol%, H2 O = 14 vol%, balance CO2 environment for 90-106 µm size bin. 143 Figure 5.16. Consistent range for US char (top) O2 = 24 vol% and (bottom) O2 = 36 vol%, both with H2 O = 14 vol% and balance CO2 environments for the 53-63 µm size bin. CHAPTER 6 CONCLUSIONS AND FUTURE WORK I studied two bricks of different scales in the CCMSC validation hierarchy. For each brick, I employed a six-step validation and uncertainty quantification (VUQ) methodology that is based on a bound-to-bound consistency analysis. Chapters 2 and 3 focused on the VUQ analysis of the 1.5 MWth oxy-coal (L1500) furnace and Chapters 4 and 5 focused on the new Reacting Particle and Boundary Layer (RPBL) char oxidation model. 6.1 6.1.1 Conclusions L1500 furnace I performed the VUQ analysis for three quantities of interest (QOIs), heat removal by cooling tubes, wall temperature, and O2 mole fraction, using experimental data from the L1500 furnace and simulation data from Arches simulations of the furnace. I selected two models for analysis in this study, the ash deposition model and the char oxidation model. I created an input and uncertainty map for these models and chose five model parameters with high priority for a sensitivity analysis. The experimental data I selected for this analysis were collected in the L1500 furnace in February of 2015 for the 100 % swirl burner operating condition. A STAR-CCM+ simulation of the burner was performed to compute the velocity profiles at the burner tip. These profiles were then used as an input boundary condition in the Arches simulation. The sensitivity analysis showed that the most sensitive parameters across the three different types of measurements were two ash deposition parameters, the thermal conductivity (k deposit ) and emissivity (ε w ) of the deposit. This analysis also showed that the char oxidation parameters had a small impact on the QOIs in comparison with the ash deposition parameters. I estimated the total error in the experimental data by adding the sampling and the systematic errors. The sampling error was computed from 30 minutes of data using a 145 t-distribution with a 95 % confidence interval [63]. I estimated the systematic error in the wall temperature with a 2D heat transfer model of the L1500 wall. I obtained a value of ∼ ± 115 [K ]. I estimated the systematic error in the heat removal data using calibration information for the water flow rates and the errors in the temperature measurements of the water inflow and out flow. I assumed that the total error for the O2 mole fraction was equal to the sampling error. I performed a consistency analysis with three parameters: k e f f (combination of k deposit and ε w ), burner swirl parameter ( f s ), and coal feed rate (ṁcoal ). Consistency was found between the simulation and all the QOIs. The consistent ranges for the parameters were: k e f f w , f s = 0.0-0.1, and ṁcoal = 275 - 283 lbh . The wall temperature uncertainty = 0.97-1.51 mK was reduced from ± 115 K to ∼ ± 16-26 [K ]. The heat removal uncertainty was reduced from ∼ ± 1.06e4 [W ] to ∼ ± 0.58e4 [W ]. The fact that the f s value was reduced by 90 % from the nominal value ( f s =1) to obtain consistency indicates a potential issue with the burner inlet condition. 6.1.2 Char oxidation I developed the new RPBL model for char oxidation. This model is based on SKIPPY [34] but has an added transient component. RPBL solves mass, energy, and momentum balance equations in a spherical char particle and its boundary layer. This model accounts for mass transport of gas species between the particle and the boundary layer using a time-dependent boundary condition. Heterogeneous and homogeneous reaction contributions are estimated with a six-step reaction mechanism and a syngas reaction mechanism respectively. The burnout process is traced through the particle and along the length of the entrained flow reactor at Sandia National Laboratories [61, 88]. For the VUQ analysis, I used experimental char oxidation data collected by Hecht [34]. I performed a sensitivity analysis on the particle temperature and velocity for nine parameters and selected five for the consistency analysis: particle diameter (d p ), void fraction (φinitial ), carbon fraction of the particle (Yc ), particle emissivity (e p ), and boundary layer thickness ( rin f rp ). The particle temperature and/or velocity were most sensitive to the first four parameters. I included the fifth parameter, compute this value. rin f rp , because I do not have a model to 146 I computed the uncertainty in the experimental data using the approach of Oberkampf and Roy [63] to calculate sampling error. The ∼ 100 measurements at each position were used to compute a particle temperature, velocity, and diameter with upper and lower bounds. The consistency analysis between particle temperatures and velocities computed with RPBL and the experimental values was performed for 36 groups of data that correspond to three types of char (char obtained from Utah Skyline (US), Illinois # 6 (I6), and Black Thunder (BT) coals), two gas compositions (O2 vol% = 24, 36; H2 O vol% = 14; balance CO2 ), and 6 size bins (53-63 µm, 63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm, and 125-150 µm). I found consistency for all the 36 groups of data. From the consistency analysis, I learned that the US and I6 chars have similar behavior with respect to particle temperature and velocity while BT char produces higher temperatures. When the O2 mole fraction is increased, the measured particle temperature also increases. This behavior is captured by the RPBL model, which also captures the effect of varying the particle size of the three chars. 6.2 Future work From the VUQ analysis performed on the L1500 furnace, I identified sources of uncertainty in the wall temperature and heat removal measurements. As a result, several changes have been made to the L1500 furnace. One, new thermocouples were installed in the top wall of every section to measure the wall temperature, and the material that was poured around the thermocouples had the same heat transfer properties as the materials in the L1500 walls. Two, the cooling tubes in sections 1 and 2 were replaced with cooling panels and a soot blowing system. With this system, better control of the ash deposits on the cooling panels is obtained and the simpler geometry is easier to simulate with Arches. I have the following recommendations from the L1500 VUQ analysis: • Preliminary studies with the 2D wall heat transfer model show that steady state conditions are reached after more than 40 h. This length of time is two orders of magnitude higher than the simulation time for an LES simulation (approximately 30 s). Transient wall effects on measurements of interest (mid-wall temperature) should be studied using a VUQ approach with the 2D wall heat transfer model. 147 • With the changes to the L1500, wall temperatures data are available along the entire length of the furnace. The entire furnace domain should be included in future simulations. • There are known leaks in the L1500 burner. The effect of leaks on the fluid dynamics in the furnace needs to be studied. In the Arches simulation presented in this dissertation, I included the effect of leaks only in the composition of oxidant. • The velocity profiles of the fluid flow at the burner exit need to be experimentally characterized. In addition, the effect of leaks between the vanes in the inner secondary and outer secondary streams should be evaluated. In the VUQ analysis of the RPBL model, I compared RPBL data and single particle experimental data from Sandia National Laboratories. In these experiments, only 100 particles were measured at each position, which is a relatively small sample given the heterogeneous nature of coal. For this reason, I chose a wide range of values for the initial mass fraction of carbon in the particle and for the initial void fraction in the VUQ analysis. I did not use the proximate analysis data for the initial mass fraction of carbon because char yield depends on the particle heating rate and the gas temperature, and these two conditions are different between the proximate analysis test and the char produced in the drop tube furnace. Therefore, my recommendations from the RPBL model VUQ analysis include the following: • The number of measurements at each position should be increased in order to reduce the sampling error. • Coal devolatilization that predicts void fraction and carbon mass fraction needs to be developed. • A better classification system for characterizing the size distribution of particles is needed. • In this dissertation, I did not compare RPBL outpout with offline measurements such as fixed carbon for several reasons. First, a collection probe was used to collect 148 particles for the offline analysis. This probe affects the profiles of gas temperature, velocity, and composition surrounding the particle and therefore needs to be accounted for in the reaction simulation. Second, these experiments used a higher flow rate of char than the single particle experiments to collect sufficient char for analysis. The RPBL model is a single particle model that does not account for the effect of other particles. For the offline measurements, a better approach would be use Arches with DQMOM for the particle phase to simulate the experiment. However, for this approach, a simpler model than RPBL must be used. REFERENCES [1] A. C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers." In Scientific Computing, 1983. [2] J. Ahn, R. Okerlund, A. Fry, and E. G. Eddings, Sulfur trioxide formation during oxycoal combustion, International Journal of Greenhouse Gas Control, 5 (2011), pp. S127- S135. [3] P. J. Ashman and B. Haynes, Improved techniques for the prediction of NOx formation from char nitrogen, (1999). [4] ASME, Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer ASME V & V 20-2009, tech. rep., 2009. [5] J. M. Aughenbaugh and C. J. J. Paredis, The Value of Using Imprecise Probabilities in Engineering Design, Journal of Mechanical Design, 128 (2006), p. 969. [6] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. a. Cafeo, J. Cavendish, C.H. Lin, and J. Tu, A Framework for Validation of Computer Models, Technometrics, 49 (2007), pp. 138-154. [7] S. K. Bhatia and D. D. Perlmutter, A random pore model for fluid-solid reactions: I. Isothermal, kinetic control, AIChE Journal, 26 (1980), pp. 379-386. [8] A. Brink, D. Lindberg, M. Hupa, M. E. De Tejada, M. Paneru, J. Maier, G. Scheffknecht, A. Pranzitelli, and M. Pourkashanian, A temperature-history based model for the sticking probability of impacting pulverized coal ash particles, Fuel Processing Technology, 141 (2016), pp. 210-215. [9] P. N. Brown and A. C. Hindmarsh, Reduced storage matrix methods in stiff ODE systems, Applied Mathematics and Computation, 31 (1989), pp. 40-91. [10] CCMSC, http://ccmsc.utah.edu/. [11] H. Childs, E. Brugger, B. Whitlock, J. Meredith, S. Ahern, D. Pugmire, K. Biagas, M. Miller, C. Harrison, G. H. Weber, H. Krishnan, T. Fogal, A. Sanderson, C. Garth, E. W. Bethel, D. Camp, O. Rübel, M. Durant, J. Favre, and P. Návratil, VisIt: An end-user tool for visualizing and analyzing very large data, High Performance Visualization-Enabling Extreme-Scale Scientific Insight, (2012), pp. 357- 372. [12] M. E. Coltrin, R. J. Kee, and F. M. Rupley, Surface CHEMKIN (Version 3.7): A FORTRAN package for analyzing heterogeneous chemical kinetics at a solid-surface-gas-phase interface, 1990. 150 [13] B. Debusschere, K. Sargsyan, C. Safta, and K. Chowd hary, The Uncertainty Quantification Toolkit ( UQTk ), in Handbook of Uncertainty Quantification, R. Ghanem, D. Higdon, and H. Owhadi, eds., Springer International Publishing, 1 ed., 2017, ch. The Uncert. [14] O. Diaz-Ibarra, J. Spinti, A. Fry, B. Isaac, J. N. Thornock, M. Hradisky, S. Smith, and P. J. Smith, A validation / uncertainty quantification analysis for a 1 . 5 MW oxy-coal fired L1500 furnace : Part A Sensitivity Analysis, In preparation for Journal of Verification, Validation and A validation (Chapter 2), (2017). [15] O. H. Diaz Ibarra, Reacting Particle and Boundary Layer (RPBL) model: Formulation, In preparation for submission to Combustion and Flame (chapter 4), 0 (2017). [16] O. H. Diaz Ibarra, Reacting Particle and Boundary Layer (RPBL) model: Validation and uncertainty quantification, In preparation for submission to Combustion and Flame (chapter 5), 0 (2017). [17] O. H. Diaz Ibarra, J. Spinti, A. Fry, B. Schroeder, J. N. Thornock, B. Isaac, D. Harris, M. Hradisky, S. Smith, E. Eddings, and P. J. Smith, A Validation/Uncertainty Quantification Analysis of a 1.5 MW Oxy-Coal Fired Furnace, AFRC 2015, (2015). [18] N. Erland, L. Haugen, R. E. Mitchell, and M. B. Tilghman, A comprehensive model for char particle conversion in environments containing O 2 and CO 2, Combustion and Flame, 162 (2015), pp. 1455-1463. [19] R. Feeley, Fighting the curse of dimensionality: A method for model validation and uncertainty propagation for complex simulation models, (2008). [20] S. Ferson, What Monte Carlo methods cannot do, Human and Ecological Risk Assessment: An International Journal, 2 (1996), pp. 990-1007. [21] S. Ferson and L. R. Ginzburg, Different methods are needed to propagate ignorance and variability, Reliability Engineering and System Safety, 54 (1996), pp. 133-144. [22] S. Ferson and J. G. Hajagos, Arithmetic with uncertain numbers: Rigorous and (often) best possible answers, Reliability Engineering and System Safety, 85 (2004), pp. 135-152. [23] M. A. Field, Measurements of the effect of rank on combustion rates of pulverized coal, Combustion and Flame, 14 (1970), pp. 237-248. [24] T. Fletcher, A. Kerstein, R. J. Pugmire, M. Solum, and D. M. Grant, A chemical percolation model for devolatilization: summary, . . . Report SAND92-8207, (1992), pp. 1- 66. [25] R. O. Fox and P. Vedula, Quadrature-Based Moment Model for Moderately Dense Polydisperse GasParticle Flows, Industrial & Engineering Chemistry Research, 49 (2010), pp. 5174-5187. [26] G. F. Froment, K. B. Bischoff, and J. De Wilde, Chemical Reactor Analysis and Design, 2011. 151 [27] A. Fry, B. Adams, K. Davis, D. Swensen, S. Munson, and W. Cox, An investigation into the likely impact of oxy-coal retrofit on fire-side corrosion behavior in utility boilers, International Journal of Greenhouse Gas Control, 5 (2011), pp. S179-S185. [28] A. Fry, J. Spinti, I. Preciado, O. Diaz-Ibarra, and E. Eddings, Pilot-scale Investigation of Heat Flux and Radiation from an Oxy-coal Flame, in American Flame Research Committee International Symposium, Salt Lake City, UT, U.S.A., 2015. [29] Gems Sensors and Control, http://www.gemssensors.com/Flow/Electronic-FlowSensors/Rotor-Flow/RFO-types-Flow-Sensor. [30] D. Goodwin, Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, 2015. [31] E. S. Hecht, C. R. Shaddix, M. Geier, A. Molina, and B. S. Haynes, Effect of CO2 and steam gasification reactions on the oxy-combustion of pulverized coal char, Combustion and Flame, 159 (2012), pp. 3437-3447. [32] E. S. Hecht, C. R. Shaddix, and J. S. Lighty, Analysis of the errors associated with typical pulverized coal char combustion modeling assumptions for oxy-fuel combustion, Combustion and Flame, 160 (2013), pp. 1499-1509. [33] E. S. Hecht, C. R. Shaddix, A. Molina, and B. S. Haynes, Effect of CO2 gasification reaction on oxy-combustion of pulverized coal char, Proceedings of the Combustion Institute, 33 (2011), pp. 1699-1706. [34] S. Hecht, Single Particle Studies of Pulverized Coal Oxy-Com Bustion, PhD thesis, University of Utah, 2013. [35] T. Holland and T. H. Fletcher, Global Sensitivity Analysis for a Comprehensive Char Conversion Model in Oxy-fuel Conditions, Energy & Fuels, 30 (2016), pp. 9339-9350. [36] J. Hong, W. C. Hecker, and T. H. Fletcher, Improving the Accuracy of Predicting Effectiveness Spherical Coordinates, Combustion, (2000), pp. 663-670. [37] I. L. Hunsaker, Parallel-Distributed, reverse Monte-Carlo radiation in coupled, large eddy combustion Simulations, PhD thesis, University of Utah, 2013. [38] I. L. Hunsaker, D. J. Glaze, J. N. Thornock, and P. J. Smith, A new model for virtual radiometers, Proceedings of the asme 2012 international heat transfer conference, (2012). [39] R. Hurt, J.-K. Sun, and M. Lunden, A Kinetic Model of Carbon Burnout in Pulverized Coal Combustion, Combustion and Flame, 113 (1998), pp. 181-197. [40] R. H. Hurt and J. M. Calo, Semi-Global Intrinsic Kinetics for Char Combustion Modeling, Combustion and Flame, 2180 (2001). [41] S. Iavarone, C. Galletti, F. Contino, L. Tognotti, P. J. Smith, and A. Parente, CFD-aided benchmark assessment of coal devolatilization one-step models in oxy-coal combustion conditions, Fuel Processing Technology, 154 (2016), pp. 27-36. [42] ICSE, http://www.icse.utah.edu/multifuel-furnace/, 2016. 152 [43] F. P. Incropera, D. P. DeWitt, T. L. Bergman, and A. S. Lavine, Fundamentals of Heat and Mass Transfer, vol. 6th, 2007. [44] H. International, www.thinkHWI.com. [45] B. J. Isaac, J. N. Thornock, S. T. Smith, Smith, and P. J., Predictive next-generation pulverized-coal boiler design, in 41st International Technical Conference on Clean Coal & Fuel Systems 2016: The Clearwater Clean Coal Conference, 2016. [46] A. Jatale, P. J. Smith, J. N. Thornock, S. T. Smith, and M. Hradisky, A Validation of Flare Combustion Efficiency Predictions From Large Eddy Simulations, Journal of Verification, Validation and Uncertainty Quantification, 1 (2015), p. 021001. [47] A. Jatale, P. J. Smith, J. N. Thornock, S. T. Smith, J. P. Spinti, and M. Hradisky, Application of a Verification, Validation and Uncertainty Quantification Framework to a Turbulent Buoyant Helium Plume, Flow, Turbulence and Combustion, 95 (2015), pp. 143- 168. [48] J. S. Justin Luitjens, John Schmidt, Alan Humphrey, J. Davison de St. Germain, Todd Harman, Jim Guilkey, Charles Reid, Dan Hinckley, Jeff Burghardt, John M. Schreiner, Joseph Peterson, Jeremy Nicholas Thornock, Brian Leavy, Qingyu Meng, Jennifer Spinti, Chuck W, http://www.uintah.utah.edu/. [49] O. Karlström, A. Brink, M. Hupa, and L. Tognotti, Multivariable optimization of reaction order and kinetic parameters for high temperature oxidation of 10 bituminous coal chars, Combustion and Flame, 158 (2011), pp. 2056-2063. [50] R. J. Kee, M. E. Coltrin, and P. Glarborg, Heterogeneous Chemistry, in Chemically Reacting Flow, John Wiley & Sons, Inc., 2005, pp. 445-486. [51] , The conservation Equations, in Chemically Reacting Flow, John Wiley & Sons, Inc., 2005, pp. 67-136. [52] M. C. Kennedy and O. Anthony, Bayesian Calibration of Computer Models, Journal of the Royal Statistical Society . Series B (Statistical Methodology), 63 (2001), pp. 425-464. [53] E. Kriegler and H. Held, Utilizing belief functions for the estimation of future climate change, International Journal of Approximate Reasoning, 39 (2005), pp. 185-209. [54] T. Maffei, a. Frassoldati, a. Cuoci, E. Ranzi, and T. Faravelli, Predictive one step kinetic model of coal pyrolysis for CFD applications, Proceedings of the Combustion Institute, 34 (2013), pp. 2401-2410. [55] D. J. Maloney, E. R. Monazam, K. H. Casleton, and C. R. Shaddix, Evaluation of char combustion models : measurement and analysis of variability in char particle size and density, Proceedings of the Combustion Institute, 30 (2005), pp. 2197-2204. [56] D. L. Marchisio and R. O. Fox, Computational Models for Polydisperse Particulate and Multiphase Systems, 2013. [57] D. Merrick, Mathematical models of the thermal decomposition of coal. 2. Specific heats and heats of reaction, Fuel, 62 (1983), pp. 540-546. 153 [58] R. E. Mitchell, L. Ma, and B. Kim, On the burning behavior of pulverized coal chars, Combustion and Flame, 151 (2007), pp. 426-436. [59] M. F. Modest, Chapter 17 - The Method of Discrete Ordinates (SN-Approximation), in Radiative Heat Transfer (Third Edition), M. F. Modest, ed., Academic Press, Boston, third edit ed., 2013, pp. 541-584. [60] A. Molina, A. F. Sarofim, W. Ren, J. F. Lu, G. X. Yue, J. M. Beér, and B. S. Haynes, Effect of boundary layer reactions on the conversion of char-N to NO, N2O, and HCN at fluidized-bed combustion conditions, Combust. Sci. Technol., 174 (2002), pp. 43- 71. [61] J. J. Murphy and C. R. Shaddix, Combustion kinetics of coal chars in oxygen-enriched environments, Combustion and Flame, 144 (2006), pp. 710-729. [62] F. Nicoud and F. Ducros, Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor, Flow, Turbulence and Combustion, 62 (1999), pp. 183-200. [63] W. L. Oberkampf and C. J. Roy, Verification and Validation for Scientific Computing, 2010. [64] Omega, http://www.omega.com/pptst/PTRA PTRM.html. [65] , http://www.omega.com/prodinfo/thermocouples.html, 2016. [66] J. Pedel, J. N. Thornock, and P. J. Smith, Large Eddy Simulation of Pulverized Coal Jet Flame Ignition Using, Energy & Fuels, 26 (2012), p. 6686. [67] J. Pedel, J. N. Thornock, and P. J. Smith, Ignition of co-axial turbulent diffusion oxy-coal jet flames: Experiments and simulations collaboration, Combustion and Flame, 160 (2013), pp. 1112-1128. [68] J. Pedel, J. N. Thornock, S. T. Smith, and P. J. Smith, Large eddy simulation of polydisperse particles in turbulent coaxial jets using the direct quadrature method of moments, International Journal of Multiphase Flow, 63 (2014), pp. 23-38. [69] E. Ranzi, A. Frassoldati, R. Grana, A. Cuoci, T. Faravelli, a. P. Kelley, and C. K. Law, Hierarchical and comparative kinetic modeling of laminar flame speeds of hydrocarbon and oxygenated fuels, Progress in Energy and Combustion Science, 38 (2012), pp. 468-501. [70] C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine learning., 2006. [71] R. Rawat, J. P. Spinti, W. Yee, and P. J. Smith, Parallelization of a large scale hydrocarbon pool fire in the uintah PSE, in ASME International Mechanical Engineering Congress and Exposition, Proceedings, vol. 6, University of Utah, Salt Lake City, UT, United States, 2002, pp. 49-55. [72] H. R. Rezaei, R. P. Gupta, G. W. Bryant, J. T. Hart, G. S. Liu, C. W. Bailey, T. F. Wall, S. Miyamae, K. Makino, and Y. Endo, Thermal conductivity of coal ash and slags and models used, Fuel, 79 (2000), pp. 1697-1710. [73] A. P. Richards and T. H. Fletcher, A comparison of simple global kinetic models for coal devolatilization with the CPD model, Fuel, 185 (2016), pp. 171-180. 154 [74] J. Robert, Chemkin-II: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics, Sandia National Laboratories Report, (1991), pp. SAND89--8009B. [75] T. Russi, A. Packard, and M. Frenklach, Uncertainty quantification: Making predictions of complex reaction systems reliable, Chemical Physics Letters, 499 (2010), pp. 1-8. [76] T. J. Santner, W. B., and N. W., The Design and Analysis of Computer Experiments, Springer-Verlag, 2003. [77] B. B. Schroeder, Scale-bridging model development and increased model credibility, PhD thesis, The University of Utah, 2015. [78] O. A. Sergeev, A. G. Shashkov, and A. S. Umanskii, Thermophysical properties of quartz glass, Journal of engineering physics, 43 (1982), pp. 1375-1383. [79] C. R. Shaddix, F. Holzleithner, M. Geier, and B. S. Haynes, Numerical assessment of Tognotti determination of CO2/CO production ratio during char oxidation, Combustion and Flame, 160 (2013), pp. 1827-1834. [80] S. L. Singer and A. F. Ghoniem, Comprehensive gasification modeling of char particles with multi-modal pore structures, Combustion and Flame, 160 (2013), pp. 120-137. [81] I. W. Smith, The combustion rates of coal chars: A review, Symposium (International) on Combustion, 19 (1982), pp. 1045-1065. [82] P. Smith, R. Rawat, J. Spinti, S. Kumar, S. Borodai, and A. Violi, Large Eddy Simulations of Accidental Fires Using Massively Parallel Computers, in 16th AIAA Computational Fluid Dynamics Conference, Fluid Dynamics and Co-located Conferences, American Institute of Aeronautics and Astronautics, June 2003. [83] L. D. Smoot and P. J. Smith, Coal Combustion and Gasification, vol. 22, PLENUM, New York, first ed., 1985, ch. 4, pp. 77-110. [84] J. P. Spinti, J. N. Thornock, E. G. Eddings, P. J. Smith, and a. F. Sarofim, Heat transfer to objects in pool fires, Dev. Heat Transfer, 20 (2008), pp. 69-136. [85] A. Strugal, Empirical relationships for the determination of true density of coal chars, 79 (2000), pp. 743-753. [86] R. Taylor and R. Krishna, Multicomponent Mass Transfer, vol. 60, 1995. [87] I. M. S. T.F wall, A. Lowe, L.J. Wibberley, Mineral matter in coal and the thermal performance of large boilers, Progress in Energy and Combustion Science, (1979), pp. 1- 29. [88] D. A. Tichenor, R. E. Mitchell, K. R. Hencken, and S. Niksa, Simultaneous in situ measurement of the size, temperature and velocity of particles in a combustion environment, Symposium (International) on Combustion, 20 (1985), pp. 1213-1221. [89] M. B. Tilghman and R. E. Mitchell, Coal and biomass char reactivities in gasification and combustion environments, Combustion and Flame, 2 (2015). [90] L. J. S. Tognotti, L., The products of the high temperature oxidation of a single char particle in an electrodynamic balance, 1990. 155 [91] H. Umetsu, H. Watanabe, S. Kajitani, and S. Umemoto, Analysis and modeling of char particle combustion with heat and multicomponent mass transfer, Combustion and Flame, 161 (2014), pp. 2177-2191. [92] H. K. Versteeg and W. Malaskekera, Chapter 5 - An Introduction to Computational Fluid Dynamics, in Fluid flow handbook. McGraw-Hill . . . , vol. M, 2007, pp. 134-178. [93] P. L. Walker, R. J. Foresti, and C. C. Wright, Surface Area Studies of CarbonCarbon Dioxide Reaction, Industrial & Engineering Chemistry, 45 (1953), pp. 1703-1710. [94] T. F. Wall, S. P. Bhattacharya, L. L. Baxter, G. Richards, and J. N. Harb, The character of ash deposits and the thermal performance of furnaces, Fuel Processing Technology, 44 (1995), pp. 143-153. [95] A. Zoulalian, R. Bounaceur, and A. Dufour, Kinetic modelling of char gasification by accounting for the evolution of the reactive surface area, Chemical Engineering Science, 138 (2015), pp. 281-290.