A Validation/Uncertainty Quantification (V/UQ) Analysis for a 1.5 MW Oxy-coal Fired L1500 Furnace using a Swirling Boundary Condition

Conference Paper

A validation/uncertainty quantification analysis for a 1.5 MW oxy-coal fired L1500 furnace using a swirling boundary condition I Oscar H. Dı́az-Ibarra, Jennifer Spinti, Andrew Fry, Benjamin Isaac, Jeremy N. Thornock, Michal Hradisky, Sean Smith, Philip J. Smith1 Salt Lake City, UT1,1 , Institute for Clean and Secure Energy -University of Utah1,∗ Abstract The work described in this paper is part of the larger mission of the CarbonCapture Multidisciplinary Simulation Center (CCMSC) (http://ccmsc.utah.edu) at the University of Utah. This paper focuses on a validation/uncertainty quantification (VUQ) study performed on the 1.5 MW L1500 furnace, an oxy-coal fired facility located at the Industrial Combustion And Gasification Research Facility at the University of Utah. The L1500 is part of the overall project because it includes many of the physics present in full-scale boilers without the complications of multiple burners and very large scales. Experiments and simulations under oxy-coal combustion conditions with a swirling burner have been done in the L1500 furnace with Utah SUFCO coal in order to perform a VUQ analysis. A six-step VUQ framework is used for studying the impact of model parameter uncertainty on the quantity of interest (QOI) for the overall project, heat flux. Parameters from both the char oxidation and ash deposition models are examined. 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 2, an input/uncertainty I sensitivity analysis author Email address: [e-mail] (Institute for Clean and Secure Energy -University of Utah) URL: http://www.icse.utah.edu/ ( Salt Lake City, UT) ∗ Corresponding Preprint submitted to AFRC 2016 INDUSTRIAL COMBUSTION SYMPOSIUMJuly 29, 2016 #12;(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. The concept of an instrument model is also introduced. Keywords: LES, oxy combustion, swirl, heat flux, coal, instrument models, sensitivity analysis 1. Nomenclature ε Surface emissivity [−] qincident Incident radiation [W m−2 ] R Thermal resistance [W m−2 K −1 ] ki Thermal conductivity for layer i [W m−1 K −1 ] ∆xi Thickness for layer i [m] mc Coal off gas mass flow [kgs−1 ] mp Primary stream mass flow [kgs−1 ] ms Secondary stream mass flow [kgs−1 ] η Mixture fraction see equation [2] [−] Fp Mixture fraction see equation [3] [−] F Mixture fraction see equation [4] [−] ρ Gas density [kgm−3 ] φ Scalar u Velocity in the flow direction [ms−1 ] Ratio between Arches resolution and STAR-CCM+ resolution [−] q Radiative heat flux [W m−2 ] Ω Solid angel Nr Number of rays Ir radiative intensity in each ray ratio W m−2 θr θ View angle T Gas temperature [K] 2 #12;Gas absorption coefficient [m−1 ] ∆x Resolution [m] Io Radiative intensity from a wall [W m−1 ] εw Surface emissivity of wall [−] Tw Wall temperature [K] Stefan Boltzmann constant 5.670373x10−8 [W m−2 K −4 Heat removed by the cooling tubes [W ] Volumetric reaction rate of char consumed from oxidizer l reaction [kgm−3 s−1 ] A Particle surface area [m− 2] w Particle number density [#m−3 ] c Mixture molar concentration [kmolem−3 ] W Mixture molecular weight [kgkmole−1 ] kl Reaction rate coefficient for reaction l [ms−1 ] Char molecular weight [kgkmole−1 ] φl Stoichiometric coefficient ratio for species l [kmolechar kmolel ] kc Mass transfer coefficient [ms−1 ] Molar concentration of oxidizer l in the bulk [kmolem−3 ] rt Total volumetric reaction rate [kgm−2 s−1 ] Sh Sherwood number [−] Re Particle reynolds number [−] Sc Schmidt number [−] dp Particle diameter [m] Dom Mixture averaged diffusion coefficient of oxidizer [m−2 s] qnet Net heat flux [W M −2 ] External temperature [K] ddeposit Deposit thickness [m] kdeposit Deposit thermal conductivity [W m−1 K] Ash deposition velocity [kgs−1 ] Soot blowing time [s] k σ Qremoval rH,l WH cgO,l Tshell v tsb 3 #12;2. Introduction The Carbon Capture Multidisciplinary Simulation Center (CCMSC) (http://ccmsc.utah.edu) at the University of Utah is demonstrating the use of exascale computing with 5 verification, validation, and uncertainty quantification as a means of accelerating deployment of low cost, low emission, coal-fired power generation technologies. This effort employs a hierarchical validation approach to obtain simultaneous consistency between a set of selected experiments at different scales embodying the key physics components (large eddy simulations, multiphase flow, particle 10 combustion and radiation) of a full-scale, oxy-fired boiler. Figure 1 presents the CCMSC validation hierarchy. This paper presents validation and uncertainty quantification (VUQ) results for the 1.5 MW oxy-coal furnace brick in the laboratory scale validation level. This suite of oxycoal-fired experiments were conducted in the L1500, a 1.5 MW 15 furnace at the University of Utah. Details about the L1500 furnace and the experimental data can be found in Fry et al [1]. A first VUQ cycle for the 0% swirl condition was presented in 2015 at the AFRC conference [2]. In this paper, the preliminary phase of a second VUQ cycle at 100% swirl conditions is presented. It focuses on simulations of the 20 L1500 with a swirling burner, the collection/processing of the simulation data, and the results of a sensitivity analysis. 3. Description of VUQ Approach The National Institute of Statistical Sciences (NISS) group presented a framework for validation of computer models called Simulator Assessment and 25 Validation Engine (SAVE), which consists of a six-step procedure[3]. These VUQ tools have been explored by Schroeder [4]. In his dissertation, Schroeder presents a theoretical basis for a VUQ methodology that employs the six-step SAVE framework with a consistency analysis methodology referred to as boundto-bound consistency [5] applied in Step 5, model output analysis. A modified 30 version of the SAVE framework is used in this VUQ analysis. 4 #12;Figure 1: Validation 5hierarchy for CCMSC #12;In step 1, model output(s) are selected as evaluation criteria or quantities of interest (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 35 (model, scenario, and numerical) that may have an impact on 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. Depending on the QOIs, instrument models may be 40 required to process the simulation output for the sensitivity analysis; see 6.3. Parameters with high priority are selected as active while low priority parameters are fixed, and the active parameters are investigated further. Assuming that the uncertainty is a probability distribution, the uncertainty of the active parameters is then propagated through the model. 45 Step 3 is the collection of both experimental and simulation data. Step 3 is closely tied to step 1, selection of the QOIs, and step 4, surrogate model development. The data that is collected in Step 3 must allow the direct measurement or the calculation of the QOIs (through an instrument model) and the simulations that are performed are determined by the requirements of the surrogate model. 50 Step 4, model approximation or surrogate model development, is required when running the computer model is expensive. Step 5, analysis of model outputs, can be performed using various methodologies. The NISS group [3] framework is partially based on the Kennedy and O'Hagan Baseyian methodology [6]. The main task is the computation of the posterior distribution, which is the product 55 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, bound-to-bound consistency, was developed by Michael Frenklach and Andrew Packard at the University of California Berkeley [5]. The basic concept of this 60 consistency analysis boils down to comparing model outputs with experimental data. If their difference is less than the error in the experimental measurement, 6 #12;the simulation 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 experi65 mentalist must reevaluate the experimental methods and data to see if there are unaccounted for errors that might increase the uncertainty. Step 6 is feedback and feed-forward. In this step, a review of the I/U map is performed 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 70 reviewed to see if new data should be incorporated in a new VUQ cycle. This paper details steps 1 and 2. Future work will include the completion of the remaining steps (data collection, development of surrogate models, consistency analysis, and feedback/feed forward). This analysis was performed on the L1500 oxy-coal furnace brick seen in Figure 1. The experiments and simulations 75 are described next, followed by a description of how steps 1 and 2 are applied to this particular case. 4. L1500 Experiments The L1500 can be operated in air-fired or oxy-fired modes and can burn solid, liquid, or gaseous fuels. In oxy-fired mode, recycled flue gas (RFG) 80 is brought from the exit of the convective section back into the burners primary/secondary oxidant registers. Oxygen is supplied to the secondary and primary RFG streams just prior to entering the burner. More information about this reactor and burner can be found in [7] and [8]. A week-long experimental campaign was carried out in the L1500 where 85 both 0% swirl and 100% swirl operating conditions were tested. The L1500 reactor is 15.5m long with a 1x1m2 transversal area as shown in Figure 2. It is divided into 10 sections, and each section has several ports through which a variety of measurements can be taken. The reactor has eight sets of watercooled tubes that remove heat from the first four sections. Additionally, there 90 is a water-cooled steel grid at the furnace exit to reduce the temperature of the 7 #12;Figure 2: Photo of the L1500 reactor located at the Industrial Combustion And Gasification Research Facility [1]. combustion gases prior to entering the convection section. 4.1. Coal Characterization A Sufco Utah coal was used in the experimental campaign. The ultimate analysis for this coal is presented in Table 2. 95 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. As seen in Figure 3, data from both methods were fitted to a Rosin-Rammler distribution. From this 100 analysis, it was concluded that the PSD distribution could be approximated by three particle sizes: 15µm, 60µm, and 200µm, with mass weights of 57.4%, 26.2% and 16.4% respectively. Assuming 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 1300kg/m3 and the volume corresponding 105 to each particle size. 8 #12;Table 2: Ultimate analysis for the Sufco Utah coal. Data are averaged and normalized from the analysis of five samples. Sufco Coal % Mass C 66.89 H 4.51 N 1.17 S 0.36 O 13.6 Ash 7.88 H2O 5.58 HHV [J/kg] 27364.93 Figure 3: Experimental PSD and fitted Rosin-Rammler PSD. Sieving was performed for different times to determine the effect on the resulting PSD. Based on this analysis, only the 30-minute sieving data was used for the Rosin-Rammler fit. 9 #12;Table 3: L1500 operating conditions Stream Mass flow [kg/s] Temperature [K] Coal (ash and moisture free) 0.03534 338 Primary 0.07103 338 Inner Secondary 0.05899 533 Outer Secondary 0.23515 533 Mass Fraction Primary Inner and Outer Secondary O2 0.1684 0.2412 CO2 0.6464 0.6114 H2 O 0.1314 0.0965 SO2 0.0009 0.0009 N2 0.0529 0.0500 4.2. Operating Conditions The L1500 operating parameters are shown in Table 3. 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 110 flow rate, the gas temperatures for each of these streams, and the composition of the RFG. All of these parameters were recorded and controlled during the experimental campaign. These operating conditions were generally stable during the experiment. However, there was an air leak in the RFG stream as evidenced by the out- 115 let 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 coals moisture are included in the RFG composition shown in Table 3. In scoping simulations, it was determined that simulation results (specifi- 120 cally radiative heat flux and wall temperatures) were not sensitive to the shell temperature of the furnace. Therefore, constant values were used for each of the four sides of the furnace. Table 4 presents the values of shell temperature used 10 #12;Table 4: 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.2 in the simulations. These values were measured using a surface thermocouple during the experimental campaign. 125 The L1500 walls are composed of one layer of refractory material (Greencast w ), 94 Plus) with thickness ∆x = 0.2032m and thermal conductivity k = 2.36 mK and three layers of insulation material: Insboard 3000 (∆x = 0.0508m, k = w w 0.19 mK ); Insboard 2600 (∆x = 0.0254m, k = 0.1475 mK ); Insblock ∆x = w 0.0508m, k = 0.104 mK ). Using the physical properties of these materials ob130 tained from the manufacturer, the thermal resistance was computed using equation 1. The average thermal resistance for the wall is R = 1.02 mw2 K . R = P4 1 ∆xi i=1 ki (1) 5. L1500 Simulation Strategy The L1500 has a swirling burner that can operated in a range of modes from 0% swirl, where all-axial flow can be achieved, to 100% swirl, where flow with 135 a 60m/s tangential velocity and 40m/s axial velocity can be obtained. For the experiments used in this validation study, coal was fed with the primary stream at a rate of 0.0324 kg/s in order to produce 1 MW of heat, and the burner was set for 100% swirl. The simulation data were obtained from the coupling of Large Eddy Simu- 140 lation (LES) of the burner and of the main chamber (see Figure 4). Because the low-NOx burner is equipped with swirl blocks and has a complex geometry, 11 #12;the commercial software package STAR-CCM+ was used to build a mesh and perform numerical simulations of the burner. The average velocity profile at the burner tip was then extracted, filtered, and used as the boundary condition (e.g. 145 handoff file) for the Arches simulation of the main chamber. Arches is an LES research code developed by this paper's authors and other CCMSC researchers. 5.1. STAR-CCM+ Simulations The low-NOx burner was resolved with STAR-CCM+. This software could represent the complex geometry of the L1500 burner using an unstructured 150 mesh. A CAD file used in burner fabrication was imported into STAR-CCM+. The following process was used to fix the CAD geometry to the CFD model. First, the flange holes and swirl adjustment rod were removed from the CAD assembly and imported into STAR-CCM+. Next, the inlets were taped and the internal volume was extracted. Finally, a box measuring 1m x 1m x 4m with 155 a smaller box at the outlet was added to the geometry to represent the furnace itself (see Figure 5). The simulation was run with the operating conditions presented in Table 3. An LES approach was used for turbulence, and the subgrid scales were modeled with the WALE model [9]. The energy equation was solved with properties 160 of the gases computed assuming non-reacting, multi-component ideal gas mixtures. Only the gas phase was resolved. The burner was assumed adiabatic, a second-order implicit time solver with fixed time step of 5x10−5 was used, and a y + wall treatment approach was employed. The computational domain had approximately 47 million cells. The burner 165 simulations were run on a linux machine, Ash, owned by the Institute for Clean and Secure Energy and operated by the Center for High Performance Computing at the University of Utah. The simulations were run on 3,600 cores for two weeks, which is long enough to obtain 3 s of physical time. Data were averaged over an interval of 1 s to obtain a velocity profile. This interval was 170 much shorter than the averaging interval for furnace data because of the short residence time/high velocity through the burner. 12 #12;ARCHES Chamber STAR-CCM+ Burner 13 Figure 4: L1500 simulation coupling between Arches and STAR-CCM+ #12;14 Figure 5: STAR-CCM+ mesh. Burner and an extra volume are resolved. #12;5.2. Arches Simulations The simulations of the L1500's main chamber were performed with Arches, a component of the Uintah software suite [10]. Uintah is a solver for partial 175 differential equations on structured grids using hundreds to thousands of processors. In Arches, the turbulent flow is resolved by the filtered Navier-Stokes equations [11, 12, 13]. In this LES approach, the large scales are resolved, and the subgrid scales are modeled with the dynamic Smagorinsky model. For the 180 LES simulations described in this paper, first order discretization in time and a hybrid convection scheme (upwind and central differences) were used [14], . The solid (coal) phase was modeled using the direct quadrature method of moments (DQMOM) with three environments [15]. The internal coordinates were the raw coal mass, the char mass, the particle enthalpy, maximum par- 185 ticle temperature, and the particle velocity (in x, y, and z coordinates). In addition, the weights for the three environments were solved. An upwind discretization scheme was used to solve the equations for the internal coordinates and for the weights. The particle phase was coupled with the gas phase through source terms in the equations for momentum, heat and mass. More detail about 190 the DQMOM implementation in Arches is presented in [13]. The coal particle physics models included a new devolatilization model based on the Chemical Percolation Devolatilization (CPD) model [16] and the char oxidation model discussed below in the subsection 5.3. More detailed information about this new devolatilization model can be found in [4]. 195 The gas-phase reactions in the system were modeled using a three-stream mixture fraction approach. The three streams were the primary stream (mp ), the secondary stream (ms ), and the coal off gas (mc ). The mixture fractions based on these three streams are defined in equations [2], [3] and [4]. η= mc mp + ms + mc 15 (2) #12;Fp = mp mp + ms + mc (3) mp mp + ms (4) F = Transport equations were solved for η and Fp ; F was computed from the 200 other two mass fractions. A lookup table based on equilibrium chemistry assumptions was tabulated a priori as a function of these three mixture fractions and of the system enthalpy. The discrete ordinates (DO) method was used to compute radiation in the simulation. For the cases in this paper, S8 quadrature, representing 40 discrete 205 directions, was used. The radiation equations were solved every 20 iterations. The L1500 simulations were run on two linux machines: Ash, mentioned above, and Syrah, a machine operated by Lawrence Livermore National Laboratory. With a 15mm mesh resolution, the computational domain had 2,255,610 cells. The simulations were run on 217 cores long enough to obtain 30 s of phys- 210 ical time. On ash, this required 100 hours of run time. The simulations were run with the operating conditions presented in Table 3. 5.3. Char Oxidation The char oxidation model implemented in Arches includes heterogeneous reaction at the particle surface, mass transport of oxidizer from the bulk gas 215 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 equation [5]. rH,l = (Aw)2 cW kl WH φl kc cgO,l AwcW (kl + kc ) + rt (5) For the L1500 simulation, two global reactions were considered: one oxidation and one gasification reaction. The oxidation reaction with O2 as oxidizer −El 220 was Cs +O2 → CO2 . The reaction constant in Arrhenius form was kl = Al e RT . 16 #12;The two parameters in this equation were AO2 and EO2 . The gasification reaction with CO2 as oxidizer was Cs + CO2 → 2CO. As with the oxidation reaction, the reaction constant kl had two parameters, ACO2 and ECO2 . The mass transfer coefficient kc is obtained using a Sherwood number correlation 225 with a correction factor as shown in equation [6]. Sh = 2 + 0.6Re1/2 Sc1/3 = kc dp Dom (6) 5.4. Ash Deposition Arches uses equation [7] as a boundary condition for the walls. In this equation, Tw is the internal (hot face) temperature of wall and Tshell is the external (cold face) temperature. For the refractory wall, Tshell is the outside 230 wall temperature and for the cooling tubes, it is the water temperature inside the tubes. Also, 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 compute with equation [8]. In this equation, the 235 first term is the resistance produced by the refractory and insulation layers in the furnace wall or by the metal in the cooling tube. The second term is the resistance caused by the ash layer. The thermal conductivity of the deposit, kdeposit , is an input parameter, and deposit thickness, ddeposit , is computed with the ash deposition model represented by equation [9]. In this model, the 240 ash deposition velocity (v) is computed from a probability model [17] and tsb , the time since the last soot blowing, is an input parameter. This deposition model uses three regimes . The first regime is when Tw is lower than the ash fusion temperature Tslag ; the deposit thickness ddeposit is computed with the ash deposition model [9]. If Tw is equal to or greater than Tslag , ddeposit is 245 computed with equation [10]. If the computed value of ddeposit is greater than zero, then Tw is equal to Tslag and the model is in regime 2. If the computed value of ddeposit is less than zero, ddeposit = 0 and Tw is computed with equation [7] and the model is in the regime 3. 17 #12;qnet = (Tw − Tshell ) = ε(qincident − σTw4 ) R (7) Nlayer R= X ∆xi ddeposit + k k i deposit i=1 (8) ddeposit = vtsb (9) Tslag − Tshell (10) qnet max − Rwall Where qnet max is computed from equation [15] with Tw replaced with Tslag PNlayer ∆xi and Rwall replaced with i=1 ki . ddeposit = kdeposit 250 5.5. Handoff Plane Because of the difference in grid size between the STAR-CCM+ and Arches simulations, the velocity field computed in STAR-CCM+ needed to be filtered for use in Arches. This filtering process has seven steps. 255 • First is to extract the velocity vectors from STAR-CCM+ on a structured grid using nearest 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 simulation (15x10−3 m). These choices result in structured grid 260 of 419x 419 points with a resolution of 5x10−4 m. • Second is to map the STAR-CCM+ information to a 2-D array, adding zeros to the array where data was not extracted (outside of the burner). • Third is to create a 2-D array for the fraction of the primary stream, Fp , with the condition Fp = 1 for r ≤ 0.0387m (radius of the primary inlet) 265 and Fp = 0 for r > 0.0387m. • Fourth is to obtain density, ρ,from the equilibrium chemistry lookup table described above using Fp and assuming that there is no heat loss in the burner. 18 #12;• Fifth is to compute the components of the velocity in the flow (u) and the 270 tangential (v and w) directions. To do this, the mass flux at the inlet, ρu (units of kg m2 s ), is filtered using equation [11], where φ = 1 for u, φ = 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 STARCCM+. In order to obtain u, ρu is divided by ρ. For v and w, ρuv and 275 ρuw are divided by ρu. • Sixth is to 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 6. On the left are the three velocity components from the STAR-CCM+ simulation at the exit plane of the burner. On the 280 right are the filtered velocity components used in the Arches simulation. • Seventh is to create input (handoff) files for Fp and enthalpy at the exit plane of the burner to be used in the Arches simulation. These files are obtain using the lookup table. For the DQMOM inlet conditions, the velocity of the particles is assumed equal to the gas velocity and a constant 285 coal mass flow rate is assumed for each cell. (j+1)ratio (i+1)ratio (ρuφ)i,j = (ρuφ)m,n ratio2 n=j∗ratio m=i∗ratio X X (11) 5.6. Domain Size and Mesh Resolution In previous simulations of the entire furnace geometry at a resolution of 16mm for the 0% swirl case, it was determined that gas temperatures, velocities and chemical compositions were relatively constant after section 6 [2]. Also, 290 differences between simulation results from the full length (15.5m) and from simulations using a shortened domain ( 7m) were minimal. Hence, the computational geometry for the analysis that follows is 7m. The impact of mesh resolution was also tested by performing an Arches simulation with 12mm resolution. While the differences in the wall temperature 19 #12;Figure 6: Velocities at plane of the burner tip; left is STAR-CCM+ with a resolution of 0.5x10−3 m , right is Arches with a resolution of 15x10−3 m; ratio 30 . 20 #12;Figure 7: 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 12mm. 295 profiles for the 12mm and 16mm resolution cases were small, there were larger differences in gas temperatures near the burner. Therefore, a 15mm resolution has selected for the sensitivity analysis discussed in this paper and a resolution of 10mm for the yet-to-be completed consistency analysis. The shortened geometry of the L1500 with a 15mm resolution is presented in 300 Figure 7; this geometry includes the cooling tubes and the 10cm 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. 21 #12;6. VUQ Analysis 305 6.1. Quantities of Interest To perform a VUQ analysis, the quantity (or quantities) of interest (QOI) and the system parameters (scenario, model, numerical) that have a first order impact on the QOIs are identified. In this experimental dataset, the QOIs all relate to heat flux. The QOIs are heat flux measurements from three narrow 310 angle radiometers (sections 1, 2, and 3), five wall temperature measurements (sections 1, 2, 3, 4, 6), and heat removal from eight 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 serving as a heat flux gauge was installed in the port opposite to 315 the radiometer to measure total heat flux to the wall and 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 320 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 8). 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- 325 measured mass flow rates and inlet/outlet temperatures of the water flowing through the eight cooling tubes in the first four sections. 6.2. I/U Map for Char Oxidation and Ash Deposition Models For this study, we desired to know the impact of the two models described previously on the QOIs: the char oxidation model and the ash deposition model. 330 The parameters and their associated uncertainty ranges used in these models are presented in Table 5. In this table, the overall priority means the relative importance of that parameter on the QOIs. Because of the computational cost 22 #12;Figure 8: Thermocouple placement in furnace wall. Section 4 is shown but the placement is similar for all thermocouple measurements. of an LES simulation, we must select the most important parameters for further study in the sensitivity and consistency analyses. 335 In the ideal case a sensitivity analysis should be performed to assign a priority value to the parameters in the IU map. However, running a sensitivity analysis with 14 parameters is too expensive for a LES simulation, even with a small case such as the L1500 coarse mesh described above. If we use two simulations for each dimension and combine all the parameters, a total of 214 = 16384 340 simulations are needed. Consequently, we carried out a pre-selection of parameters to assign priorities. In Table 5, the most important parameters are assigned a priority of 6 and the least important a priority of 1. In the char oxidation model, the parameters are the kinetic constants for the 345 two reactions, an oxidation reaction and a gasification reaction (see subsection 5.3). Using 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 23 #12;6 3 3 1 E O2 ACO2 ECO2 Initial deposit thick- 24 6 3 ts Tshell 3 3 3 6 6 Ash density Tslag kdeposit wallemissivity i=1 ∆xi ki PNlayer 3 Erosion thickness ness 6 0.3 0.1 120666 17 58 min Priority Range AO 2 Parameter 1 1.8 137322 20 68 max 1516.5 1000 1.02 1 or 5 0 53700 1390 value s atm −1 O2 ] factor for oxidation reaction factor for gasification reaction Surface emissivity [-] Thermal conductivity for deposit [W m−1 K −1 ] Ash fusion temperature [K] for Sufco coal Constant [kg/m2 ] Thermal resistance in the L1500 walls [W m−2 s−1 ] temperature of L1500 Average temperature inside cooling tube or external surface ing tubes [s] Soot blowing time or time since ash was removed from cool- thing else Maximum thickness: 1cm for cooling tubes, 5cm for every- Deposit thickness [m] at the start of the simulation Activation energy for gasification reaction [kcal/mole] [gcm−2 s−1 atm−1 O2 ] Pre-exponential Activation energy for oxidation reaction [kcal/mole] [gcm −2 −1 Pre-exponential Nominal Basis/Comments Table 5: Input/uncertainty map for char oxidation and deposition model #12;parameters. The remainder of the parameters in Table 5 relate to the ash deposition 350 model. The initial deposit thickness, set to 0 for the L1500 simulations, is used in the simulation before the ash deposition model is turned on. While this parameter could have an effect on the outputs of the simulation if the simulation were not run for enough time to flush the initial condition, in this case the simulation time was 30 s, which is four times the residence time of 355 approximately 7 s. Therefore, a priority of 1 was given to this parameter. The erosion thickness, the maximum deposit thickness allowed in the simulation, was given a priority of 3. It is difficult to assign a value to the soot blowing time (ts ) 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 360 was assigned to this parameter. The shell temperature (Tshell ) is the exterior temperature for the L1500 wall and the inside tube temperature for the cooling tubes. This parameter was given a priority of 3 because the water temperature did not change much in the cooling tubes and previous analysis has shown that L1500 shell temperatures do not have a big impact on the internal furnace 365 temperature. PNlayer The thermal resistance ( i=1 ∆xi ki ) due to the refractory and the insu- lation 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 370 analysis. The effective thermal conductivity kdeposit and Tslag were analyzed in a VUQ study of a 15 MWth oxy-coal boiler. 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 the nominal value. For the kdeposit parameter, a priority of 6 was 375 assigned because of the lack of information about the thermal conductivity of the ash deposits on the cooling tubes. The wallemissivity was given a priority of 6, because there is almost no information about this parameter in this reactor. The active parameters for this analysis are the parameters with a priority 25 #12;of 6. This set includes two parameters related to char oxidation, AO2 and EO2 , 380 and three parameters related to ash deposition, kdeposit , wallemissivity , and ts . The other parameters in Table 5 are fixed, and the nominal value is used in the simulations. A sensitivity analysis is performed to reduce the active parameter set from five to three for the consistency analysis. The uncertainty ranges for the char oxidation parameters AO2 and EO2 , 385 which are the reaction parameters of the oxidation reaction C + O2 → CO2 , were taken from Smoot and Smith [18]. Smoot and Smith reported small variations in activation energy, EO2 , for three U.S. coals (17-20 kcal/mole). In the case of AO2 , the range of 58-68 gcm−2 s−1 atm−1 O2 was selected based on the reported value of 60 gcm−2 s−1 atm−1 O2 for Illinois No. 6 (similar coal type and 390 composition). The uncertainty ranges for the ash deposition model parameters were selected as follows. The range for kdeposit was bounded by the maximum and minimum values of the experimental data presented by Rezaei and coworkers [19]. The wallemissivity range was bounded by the maximum and minimum 395 values reported by Wall and coworkers [20, 21]. The parameter ts 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), 400 and the minimum value was the total time without soot blowing less the first day of experiments where the coal feed rate was lower. 6.3. Instrument models 6.3.1. Flux Measured by Radiometer Three narrow-angle radiometers were used to measure the radiative heat 405 flux in sections 1, 2 and 3 in the L1500. The simulation instrument model for the radiometer used a reverse-Monte Carlo ray tracing approach to compute the radiative heat flux by summing up the radiative intensities over all the rays comprising the field of view, θ of the radiometer as seen in equation [13] to 26 #12;compute the solid angle Ω. The radiative intensity in each ray, Ir , was computed 410 with equation [14]. For this analysis, Nr = 1; because of the coarseness of the computational mesh, one ray was sufficient to account for the narrow angle view of the radiometer [22]. In order to compute the radiative intensity with equation [14], the gas temperature T , the gas absorption coefficient k, and the mesh resolution ∆x were 415 obtained from the simulations. These values were extracted along one ray that extended from the wall opposite the radiometer to the radiometer lens. 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 420 emissivity, εw , were used; see equation [15]. q=Ω Nr 1 X Ir cos(θr ) Nr r=1 Ω = 2π(1 − cos(θ)) Ir = N −1 X j=0 j−1 N X X σ 4 ∆xki ) Tj exp(− ∆xki )(1 − exp(∆xkj )) + Io exp(− π i=0 i=0 σ Io = εw Tw4 π (12) (13) (14) (15) More detailed information about the instrument model can be found in [22, 23]. 6.3.2. Heat Removal by Cooling Tubes For each computational cell in the simulated cooling tubes, the heat removal 425 (Qremoval ) was computed with equation [16] and then added over the whole cooling tube to obtain the total heat removal. Thus, in equation [16], A corresponds to the surface area of the cell faces that are in contact with the gas, 27 #12;Tw is the surface temperature of the cooling tube, qincident is the incident heat flux, and ε is the surface emissivity of the tube. Qremoval = N X ε(qincident − σTw4 )A (16) n=0 430 6.3.3. Wall Temperature Measured by Thermocouple In the Arches simulation, it was assumed that the incident heat flux at the location of the wall thermocouple was not affected by the thermocouple; that is, the thermal conductivity of the computational cell at the thermocouple location was the wall thermal conductivity. For this reason, a thermcouple instrument 435 model is needed to correct for the heat lost due to differences in material properties between the ceramic sheath and the wall as shown schematically in Figure 8. However, for the sensitivty analysis described below, the wall temperature computed in the simulation at the thermocouple location was used as output in the sensitivity analysis described in section 7. 440 7. Sensitivity Analysis In the I/U map (Table 5), five variables were chosen as priority 6, which means these variables are active in this study. However, a five- dimensional study with Arches simulations is too expensive. Therefore, a sensitivity analysis was used to reduce the number of dimensions in this study. 445 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 [24]. UQTk uses the app pce_sens to compute total and main sensitivity indeces using a polynomial chaos (PC) surrogate model. 450 7.1. Generating PC Surrogate Models 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. UQTk app generate_quad generated 32 quadrature points of ξi = 28 #12;+/ − 0.57735027 with weights of w = 0.0625 for each dimension. The variable 455 ξi is mapped to physical input space for the Arches simulations using equation [17], where ai and bi are the bounds in the uncertainty interval presented in the I/U map (Table 5). λi = ai − bi bi − ai + ξi 2 2 for i = 1, ..., d. (17) A set of Arches simulations was run at the 32 quadrature points. All the simulations were run for more than 30 s; average fields were computed in an 460 interval of approximately 5 s. The radiometer instrument model (subsection 6.3.1) uses the gas temperature and the absorption coefficient from the Arches simulation to compute the radiative heat flux. The average field for these two variables was extracted at the radiometer location and read it into a python script that computed the ra- 465 diative heat flux. The radiative heat flux was the output used to generate the PC surrogate model. The heat removed by the cooling tubes (subsection 6.3.2) was computed from Arches simulation output using the visualziation tool Visit. In Visit, the total heat removal was calculated using the average fields of qincident and Tw . 470 This output is used to generate the PC surrogate model. The thermocouple instrument model (subsection 6.3.3) was not used in this sensitivity analysis. Rather, the average wall temperatures at the thermocouple locations were extracted from Arches and used to generate the PC surrogate model. 475 The coefficients of the PC surrogate model were then used to compute the main and the total sensitivity indeces using the UQTk app pce_sens. 7.2. Results In the Figure 9, the main sensitivity indeces for the radiometer instrument model are presented. The two most sensitive parameters with respect to the 480 radiometer measurements are wallemissivity and kdeposit ; both parameters are part of the ash deposition model. It is interesting that the main sensitivity 29 #12;1.2 Ao2 EO2 kdep wallemi Soot_blowing 1.0 Sensitivity 0.8 0.6 0.4 0.2 0.0 rad_1 rad_2 rad_3 Figure 9: Main sensitivity index for radiative heat flux measured by the radiometers index for ts is low by comparison. The main sensitivity indeces for the AO2 and EO2 are also low, meaning these parameters have little influence on the radiometer measurement. Thus, for the radiative heat flux measurement, the 485 five-dimensional (parameter) study can be reduced to a two-dimensional study. In Figure 10, the main sensitivity indeces for the heat removed by the cooling tubes are presented. The sensitivity indeces have a similar behavior to those obtained in the radiometer analysis. The wallemissivity and kdeposit parameters have bigger indeces while the sensitivity index for ts is small for the cooling 490 tubes. The sensitivities of the char oxidation parameters are small as well. Just 30 #12;1.2 Ao2 EO2 kdep wallemi Soot_blowing 1.0 Sensitivity 0.8 0.6 0.4 0.2 0.0 HR_1 HR_2 HR_3 HR_4 HR_5 HR_6 HR_7 Figure 10: Main sensitivity index for the heat removed by the cooling tubes 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 11, the main sensitivity indeces for the wall temperature at the thermocouple locations are presented. For all locations, the biggest sensitivity 495 index is for wallemissivity . However, the sensitivity index for kdeposit is low at location 1 (W HT1 ) compared with the much higher index values at the other locations. Thus, we conclude that the consistency analysis study can be reduced from five to two dimensions. The insensitivity to tsg (which determines deposit thickness) for all QOIs 31 HR_8 #12;1.2 Ao2 EO2 kdep wallemi Soot_blowing 1.0 Sensitivity 0.8 0.6 0.4 0.2 0.0 WT_1 WT_2 WT_3 WT_4 Figure 11: Main sensitivity index for wall temperature (thermocouple's location) computed with Arches 32 WT_5 #12;500 was unexpected, so the simulations were checked to look for an explanation. It was seen that the deposit thickness was zero for all the walls except for the cooling tubes. There are two reasons for this result. One, in all the simulations the wall temperature was higher than Tslag in sections 1 and 2, and for some simulations, wall temperatures exceeded Tslag for the entire length of the whole 505 reactor as presented in Figure 12. 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 temperature. The second reason is that for the 100 % swirl operating condition, most of the ash deposition is expected in sections 1 and 2 where the model is in 510 regime 3. In the Arches simulations, the only part of the furnace where the deposit thickness was greater than zero was on the cooling tubes. The influence of the cooling tubes deposits can be seen in the sensitivity index of the radiometers, where kdeposit has the second largest effect after wallemissivity and the where its 515 index increases with the section (see Figure 9). 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 10. In sections 1 and 4, HR1 , HR2 , HR7 , and HR8 are the only QOIs where the sensitivity 520 indeces of ts are greater than those ofkdeposit . In section 2 (HR3 and HR4 ), kdeposit has a bigger sensitivity index than wallemissivity . Finally, in section 3, kdeposit still has a strong influence. 8. Conclusions The work described in this paper represents the first part of a consistency 525 analysis study of experimental measurements taken in the L1500 oxy-fired pulverized coal combustor 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 33 #12;1900 Tslag Simulation Wall Temperature [K] 1800 1700 1600 1500 1400 1300 1200 0 1 2 3 4 Section No 5 Figure 12: Wall Temperature (thermocouple's location) computed with Arches, horizontal red line is the Tslag 34 6 7 #12;burner was computed using the commercial software Star-CCM+ and then pro530 vided as an input to the Arches simulations of the L1500 at the plane of the burner tip. The first two steps of the VUQ analysis are discussed in this paper: (1) selection of QOIs and (2) creation of an I/U map and selection of active parameters (those with a priority of 6). The active parameter selection requires 535 the development of instrument models to compute the QOIs from the simulation data and the performance of a sensitivity analysis to reduce the number of active parameters. From the sensitivity analysis, it can be concluded that of the set of five active parameters in the I/U map, the most sensitive parameters across three differ- 540 ent types of measurements (heat removal by cooling loops, heat flux measured by radiometers, and net heat flux at the wall) are kdeposit and wallemissivity . Therefore, for the consistency analysis, these two parameters will be used. 9. Acknowledgment This material is based upon work supported by the Department of En545 ergy, National Nuclear Security Administration, under Award Number(s) DENA0002375. 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, 550 including Ignacio Preciado, Teri Draper Snow, Eric Eddings and David Wagner. Finally, Oscar Dı́az-Ibarra would like to thank Dr. Debusschere and coworkers for their help with the Uncertainty Quantification Toolkit. References [1] A. Fry, J. Spinti, I. Preciado, O. Diaz-Ibarra, E. Eddings, Pilot-scale Inves555 tigation of Heat Flux and Radiation from an Oxy-coal Flame, in: American 35 #12;Flame Research Committee International Symposium, Salt Lake City, UT, U.S.A., 2015. [2] O. H. Diaz Ibarra, J. Spinti, A. Fry, B. Schroeder, J. N. Thornock, B. Issac, D. Harris, M. Hradisky, S. Smith, E. Eddings, P. J. Smith, A Vali560 dation/Uncertainty Quantification Analysis of a 1.5 MW Oxy-Coal Fired Furnace, AFRC 2015. [3] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. a. Cafeo, J. Cavendish, C.-H. Lin, J. Tu, A Framework for Validation of Computer Models, Technometrics 49 (2) (2007) 138-154. doi:10.1198/004017007000000092. 565 URL http://www.tandfonline.com/doi/abs/10.1198/ 004017007000000092 [4] B. B. Schroeder, Scale-bridging model development and increased model credibility, Ph.D. thesis, The University of Utah (2015). [5] T. Russi, A. Packard, M. Frenklach, Uncertainty quantification: Making 570 predictions of complex reaction systems reliable, Chemical Physics Letters 499 (1-3) (2010) 1-8. doi:10.1016/j.cplett.2010.09.009. URL http://dx.doi.org/10.1016/j.cplett.2010.09.009 [6] M. C. Kennedy, O. Anthony, Bayesian Calibration of Computer Models, Journal of the Royal Statistical Society . Series B (Statistical Methodology) 575 63 (2001) 425-464. doi:10.1111/1467-9868.00294. URL http://www.jstor.org/ [7] J. Ahn, R. Okerlund, A. Fry, E. G. Eddings, Sulfur trioxide formation during oxy-coal combustion, International Journal of Greenhouse Gas Control 5 (2011) S127-S135. doi:10.1016/j.ijggc.2011.05.009. 580 URL http://dx.doi.org/10.1016/j.ijggc.2011.05.009 [8] A. Fry, B. Adams, K. Davis, D. Swensen, S. Munson, 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 ControlFry, A., 36 #12;Adams, B., Davis, K., Swensen, D., Munson, S., & Cox, W. (2011). An 585 investigation into the likely impact of oxy-coal retrofit on fire-side corrosion behavior in utility boilers. International Journal 5 (2011) S179-S185. doi:10.1016/j.ijggc.2011.05.021. URL http://dx.doi.org/10.1016/j.ijggc.2011.05.021 [9] F. Nicoud, F. Ducros, Subgrid-Scale Stress Modelling Based on the Square 590 of the Velocity Gradient Tensor, Flow, Turbulence and Combustion 62 (3) (1999) 183-200. doi:10.1023/A:1009995426001. URL http://dx.doi.org/10.1023/A:1009995426001 [10] J. S. Justin Luitjens, John Schmidt, Alan Humphrey, J. Davison de St. Germain, Todd Harman, Jim Guilkey, Charles Reid, Dan Hinckley, 595 Jeff Burghardt, John M. Schreiner, Joseph Peterson, Jeremy Nicholas Thornock, Brian Leavy, Qingyu Meng, Jennifer Spinti, Chuck W, http://www.uintah.utah.edu/. URL http://www.uintah.utah.edu/ [11] P. Smith, R. Rawat, J. Spinti, S. Kumar, S. Borodai, A. Violi, Large Eddy 600 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, 2003. doi:doi:10.2514/6.2003-3697. URL http://dx.doi.org/10.2514/6.2003-3697 605 [12] R. Rawat, J. P. Spinti, W. Yee, 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. doi:10.1115/IMECE2002-33105. 610 URL https://www.scopus.com/inward/record. uri?eid=2-s2.0-78249257590&partnerID=40&md5= b8886f2bf89589bc62cccbe60641b6e9 37 #12;[13] J. Pedel, J. N. Thornock, P. J. Smith, Ignition of co-axial turbulent diffusion oxy-coal jet flames: Experiments and simulations collaboration, Combus615 tion and Flame 160 (6) (2013) 1112-1128. doi:10.1016/j.combustflame. 2013.01.022. URL http://dx.doi.org/10.1016/j.combustflame.2013.01.022 [14] H. K. Versteeg, W. Malaskekera, An Introduction to Computational Fluid Dynamics, Vol. M, 2007. doi:10.2514/1.22547. 620 [15] R. O. Fox, P. Vedula, Quadrature-Based Moment Model for Moderately Dense Polydisperse GasParticle Flows, Industrial & Engineering Chemistry Research 49 (11) (2010) 5174-5187. doi:10.1021/ie9013138. URL http://dx.doi.org/10.1021/ie9013138 [16] T. Fletcher, A. Kerstein, R. J. Pugmire, M. Solum, D. M. Grant, A 625 chemical percolation model for devolatilization: summary, . . . Report SAND92-8207 (1992) 1-66. URL http://www.researchgate.net/publication/253563635_A_ CHEMICAL_PERCOLATION_MODEL_FOR_DEVOLATILIZATION_SUMMARY/file/ e0b4952b4b428ab386.pdf 630 [17] A. Brink, D. Lindberg, M. Hupa, M. E. De Tejada, M. Paneru, J. Maier, G. Scheffknecht, A. Pranzitelli, M. Pourkashanian, A temperature-history based model for the sticking probability of impacting pulverized coal ash particles, Fuel Processing Technology 141 (2016) 210-215. doi:10.1016/ j.fuproc.2015.08.039. 635 URL http://dx.doi.org/10.1016/j.fuproc.2015.08.039 [18] L. D. Smoot, P. J. Smith, Coal Combustion and Gasification, 1st Edition, Vol. 22, PLENUM, New York, 1985, Ch. 4, pp. 77-110. [19] 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, Y. Endo, Thermal conductivity 640 of coal ash and slags and models used, Fuel 79 (13) (2000) 1697-1710. 38 #12;URL http://www.sciencedirect.com/science/article/ B6V3B-4106DWB-M/2/0af95ad9bf0075b7c993cfc8edb607be [20] 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 645 Science (1979) 1-29. [21] T. F. Wall, S. P. Bhattacharya, L. L. Baxter, G. Richards, J. N. Harb, The character of ash deposits and the thermal performance of furnaces, Fuel Processing Technology 44 (1995) 143-153. doi:10.1016/0378-3820(94) 00112-7. 650 [22] I. L. Hunsaker, D. J. Glaze, J. N. Thornock, P. J. Smith, A new model for virtual radiometers, Proceedings of the asme 2012 international heat transfer conference. [23] I. L. Hunsaker, Parallel-Distributed, reverse Monte-Carlo radiation in coupled, large eddy combustion Simulations, Ph.D. thesis, University of Utah 655 (2013). [24] B. Debusschere, K. Sargsyan, C. Safta, K. Chowd hary, The Uncertainty Quantification Toolkit ( UQTk ), in: R. Ghanem, D. Higdon, H. Owhadi (Eds.), Handbook of Uncertainty Quantification, 1st Edition, Springer International Publishing, 2017, Ch. The Uncert. 39