A validation/uncertainty quantification (VUQ) study was 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. 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 heat flux: incident heat flux, heat removal by cooling tubes, and wall temperatures. We describe the experimental facility, the operating conditions, and the data collection process. To obtain the simulation data, we utilized two tools, star-ccm+ and Arches. The star-ccm+ simulations captured flow through the complex geometry of the swirl burner while the Arches simulations captured multiphase reacting flow in the L1500. We employed a filtered handoff plane to couple the two simulations. In step two, we developed an input/uncertainty (I/U) map and assigned a priority to 11 model parameters based on prior knowledge. We included parameters from both a char oxidation model and an ash deposition model in this study. We reduced the active parameter space from 11 to 5 based on priority. To further reduce the number of parameters that must be considered in the remaining steps of the framework, we performed a sensitivity analysis on the five parameters and used the results to reduce the parameter set to two.

## 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. 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 (LES), multiphase flow, particle combustion and radiation, are represented at these various scales. The CCMSC validation hierarchy is presented in Fig. 1. In this paper, we present the first phase of a validation and uncertainty quantification (VUQ) analysis of the 1.5 MWth oxy-coal furnace (L1500) brick in the burner-scale validation level. This study was part of a larger study encompassing the entire hierarchy that identified char oxidation and ash deposition as the two models with primary impact on the quantity of interest (QOI), heat flux. For this reason, the model parameters identified for this study are from these two models.

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. 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.

## Description of Validation and Uncertainty Quantification 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 QOI [1]. 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 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 approach was developed by Ferson et al. [2–4] and applied by other researchers [1,5,6]. 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 [7]. The bound-to-bound consistency analysis approach was developed by Frenklach and coworkers [8] and used by Oberkampf and Roy [1], Jatale et al. [9,10], and Schroeder [11]. 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 [12]. More recently, Schroeder [11] explored these VUQ tools and presented a theoretical basis for a VUQ methodology that employs the six-step SAVE framework with the bound-to-bound consistency analysis [8].

This analysis employs the modified version of the SAVE framework proposed by Schroeder and is presented in Fig. 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 [12], rational and quadratic surrogate [13], and polynomial chaos (PC) [14].

Step 5, analysis of model outputs, is performed using various methodologies. The NISS group [12] framework is partially based on the Kennedy and O'Hagan Bayesian methodology [15]. The main task is the computation of the posterior distribution, which is the 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 [8]. 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 includes both 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 reviewed to see if new data should be incorporated in a new VUQ cycle.

This paper describes how steps 1 and 2 are applied to the L1500 oxy-coal furnace. A companion paper details the completion of the remaining steps.

## L1500 Experiments

This work focuses on data collected during a week-long experimental campaign in February of 2015 during which the furnace was oxy-fired with a Utah bituminous coal (Sufco).

### L1500 Reactor and Swirl Burner.

The refractory-lined L1500 reactor, shown in Fig. 3, is 14.2 m long with a 1.05 × 1.05 m^{2} transversal area. It is divided into ten 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 Fig. 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 [16,17]. During oxy-firing, recycled flue gas (RFG) is brought from the 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.

The swirl burner used for these tests has a range of operating modes. A schematic of this burner is presented in Fig. 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 a fixed block and a moving block. To increase the tangential velocity components (swirl), the position of the moving block is changed. In Fig. 6, the positions of the blocks for 0% and 100% swirl are presented. This work presents results for the 100% swirl operating condition.

### Coal Characterization.

A Utah Sufco coal was used in this experimental campaign. The ultimate analysis for this coal is presented in Table 1.

Sufco coal | % Mass |
---|---|

C | 66.89 |

H | 4.51 |

N | 1.17 |

S | 0.36 |

O | 13.60 |

Ash | 7.89 |

H_{2}O | 5.58 |

Sufco coal | % Mass |
---|---|

C | 66.89 |

H | 4.51 |

N | 1.17 |

S | 0.36 |

O | 13.60 |

Ash | 7.89 |

H_{2}O | 5.58 |

To determine the particle size distribution (PSD) of the Sufco coal, we sampled bags of coal at different depths and then performed both a sieving and Beckman–Coulter diffraction analyses on the collected samples. The collected data from both analyses are presented in Fig. 7. Only data from the longest sieving time (30 min) are included because data at shorter sieving times were much different than the diffraction data. For this reason, we fitted only data from Beckman-Coulter diffraction analysis to the Rosin-Rammler distribution. During the experimental campaign, we minimized particle size uncertainty by selecting coal bags with similar PSDs to burn.

### Operating Conditions.

The L1500 operating parameters for this VUQ study are shown in Table 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.

Stream | Mass flow (kg/s) | Temperature (K) |
---|---|---|

Coal (ash and moisture free) | 0.03534 | 338 |

Primary (m)_{p} | 0.07103 | 338 |

Inner secondary (m)_{s} | 0.05899 | 533 |

Outer secondary (m)_{s} | 0.23515 | 533 |

Stream | Mass flow (kg/s) | Temperature (K) |
---|---|---|

Coal (ash and moisture free) | 0.03534 | 338 |

Primary (m)_{p} | 0.07103 | 338 |

Inner secondary (m)_{s} | 0.05899 | 533 |

Outer secondary (m)_{s} | 0.23515 | 533 |

Mass fraction | Primary (m)_{p} | Inner and outer secondary (m)_{s} |
---|---|---|

O_{2} | 0.1684 | 0.2412 |

CO_{2} | 0.6464 | 0.6114 |

H_{2}O | 0.1314 | 0.0965 |

SO_{2} | 0.0009 | 0.0009 |

N_{2} | 0.0529 | 0.0500 |

Mass fraction | Primary (m)_{p} | Inner and outer secondary (m)_{s} |
---|---|---|

O_{2} | 0.1684 | 0.2412 |

CO_{2} | 0.6464 | 0.6114 |

H_{2}O | 0.1314 | 0.0965 |

SO_{2} | 0.0009 | 0.0009 |

N_{2} | 0.0529 | 0.0500 |

These operating conditions were stable during the experiment. However, there was an air leak in the RFG stream as evidenced by the outlet CO_{2} concentration being lower than expected. We performed an overall mass balance and adjusted the RFG composition at the burner inlet until the outlet CO_{2} concentration matched the experimentally observed value. We then used this RFG composition to compute the composition of the primary and secondary inlet streams shown in Table 2. In this way, we added the effect of leaks in the composition.

We measured furnace shell temperatures with a surface thermocouple during the experimental campaign. 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 furnace shell temperature within the range of temperatures measured experimentally. Therefore, the averaged experimental values shown in Table 3 were used for each of the four sides of the furnace.

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 |

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 |

The cooling tube shell temperatures, representing the average of the inlet and outlet water temperatures in the cooling tubes, are presented in Table 4. The table also includes the temperature difference between the inlet and outlet flows (Δ*T*). In the L1500 simulations, we assume that the water temperature does not change along the length of the tube because Δ*T* is small (on the order of 20 K) and this small difference has a negligible effect on the heat exchange between the cooling tubes and the internal walls of the furnace.

Location | Shell Temperature (K) | ΔT (K) |
---|---|---|

Coil No 1 | 294 | 19.3 |

Coil No 2 | 294 | 20.2 |

Coil No 3 | 295 | 18.9 |

Coil No 4 | 296 | 29.3 |

Coil No 5 | 293 | 18.7 |

Coil No 6 | 297 | 26.6 |

Coil No 7 | 296 | 22.2 |

Coil No 8 | 293 | 18.2 |

Location | Shell Temperature (K) | ΔT (K) |
---|---|---|

Coil No 1 | 294 | 19.3 |

Coil No 2 | 294 | 20.2 |

Coil No 3 | 295 | 18.9 |

Coil No 4 | 296 | 29.3 |

Coil No 5 | 293 | 18.7 |

Coil No 6 | 297 | 26.6 |

Coil No 7 | 296 | 22.2 |

Coil No 8 | 293 | 18.2 |

*R*) using Eq. (2). The average value for $R=1.02(W/m2K)$.

Layer name | Δx (m) | Nominal value $(W/mK)$ | Range $(W/mK)$ | T range (K)_{w} |
---|---|---|---|---|

Greencast 94 Plus | 0.2032 | 2.03 | 2.03–2.09 | 1143–1588 |

Insboard 3000 | 0.0508 | 0.19 | 0.09–0.35 | 698–1477 |

Insboard 2600 | 0.0254 | 0.1475 | 0.09–0.22 | 589–1255 |

Insblock | 0.0508 | 0.104 | 0.07–0.15 | 478–923 |

Layer name | Δx (m) | Nominal value $(W/mK)$ | Range $(W/mK)$ | T range (K)_{w} |
---|---|---|---|---|

Greencast 94 Plus | 0.2032 | 2.03 | 2.03–2.09 | 1143–1588 |

Insboard 3000 | 0.0508 | 0.19 | 0.09–0.35 | 698–1477 |

Insboard 2600 | 0.0254 | 0.1475 | 0.09–0.22 | 589–1255 |

Insblock | 0.0508 | 0.104 | 0.07–0.15 | 478–923 |

### Selection of Quantities of Interest.

To perform a VUQ analysis, we first identify the QOIs. In this experimental dataset, the QOIs all relate to heat flux: incident flux measured by three narrow angle radiometers (sections 1–3), five wall temperature measurements (sections 1–4 and 6), and heat removal by the eight sets of cooling tubes. All data were collected at a frequency of 1 s and recorded by OPTO 22, a data acquisition hardware/software platform, for the duration of the experimental campaign.

During the experimental campaign, we measured radiative heat flux through the center port of sections 1–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, increasing the uncertainty of the radiometer measurement.

We measured wall temperature in sections 1–4, and 6, using type B thermocouples encased in ceramic sheaths that were inserted into small holes in the furnace ceiling located in the middle of each section (see Fig. 8). Each sheath was inserted until it was flush with the inside wall of the furnace.

We determined the heat removal by measuring the water mass flow rate with an RFO type electronic flow sensor and the inlet and outlet water temperatures with type K thermocouples for each of the eight sets of cooling tubes in the first four sections.

## L1500 Simulation Strategy

The simulation data for the 100% swirl experiment were obtained from the coupling of two LES simulations, one of the burner and the other of the main chamber (see Fig. 9). Because the burner is equipped with swirl blocks and has a complex geometry, we used the commercial software package star-ccm+ to build a mesh and perform numerical simulations of the burner. We then extracted and filtered the time-averaged velocity profile at the burner tip to use as the boundary condition (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 (DQMOM)) required in our analysis.

### 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 [18]. To model the fluid flow in the boundary layer, we used 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 burner tip. We assumed that the particles would not affect the gas velocity as there were no chemical reactions 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 5 × 10^{−5} s.

The burner simulations were run on a Linux cluster at the University of Utah. With approximately 30 × 10^{6} cells in the computational domain, we resolved the fluid flow down to the viscous sublayer (*y*^{+} < 3) with three grid points below *y*^{+} of 10. Each simulation required 3600 cores for two weeks in order 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.

Based on simulation results, the computed swirl numbers for the inner and outer secondary registers were 2.3 and 3.6, respectively.

### Arches Simulations.

We performed simulations of the L1500's main chamber with Arches, a component of the Uintah software suite [19]. Uintah is a framework for solving partial differential equations on structured grids using hundreds to thousands of processors.

Arches is low-Mach, pressure projection-based, variable density code solving filtered (per the standard LES approach) conservation equations for mass, momentum and energy of the gas and solid phases [20–22]. For the gas phase, subgrid fluctuations 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) [23]. 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 DQMOM with three environments [24]. Seven internal coordinates, including raw coal mass, char mass, particle enthalpy, maximum particle temperature, and 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 Refs. [20–22].

The coal particle physics models included both devolatilization and char oxidation models. The devolatilization model, detailed in Ref. [11], is based on the chemical percolation devolatilization model [25]. 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 chemical percolation devolatilization model. The char oxidation model is discussed in Sec. 5.1.

*m*), the secondary (inner + outer) stream (

_{p}*m*), and the coal off gas (

_{s}*m*). During char oxidation or devolatilization, coal off gas enters the gas phase with an elemental composition shown in Table 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 Eqs. (3)–(5).

_{c}Transport equations were solved for *η* and *F _{p}*;

*F*was computed from the other two mass fractions as $(F=Fp/1\u2212\eta )$. A lookup table based on equilibrium chemistry assumptions was tabulated a priori as a function of three independent variables:

*η*,

*F*and the linearized enthalpy [26]. The state space variables in the table include the gas temperature, mixture enthalpy, species mass fractions, mixture molecular weight, and mixture density.

_{p,}The discrete ordinates method was used to compute radiation in the simulation [27]. For the cases in this paper, we used S8 quadrature, representing 80 discrete directions.

We ran the L1500 simulations 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.5 × 10^{−5} s. The simulations were run on 217 cores for approximately 100 h. This was long enough to obtain 30 s of physical time. With a five-second averaging window, we compared the moving average to the average obtained for the entire interval from 10 to 30 s. After 15 s of simulation time, the relative error from averaging over a 5 s window dropped below 1%. We averaged over the last 5 s of simulation data for the analyses in this paper.

### 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 + was filtered for use in Arches by the following process.

- (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 between the star-ccm+ and Arches resolutions and the final resolution of the Arches simulations (15 × 10

^{−3}m). For a ratio of 30, the result is a structured grid of 419 × 419 points with a resolution of 5 × 10^{−4}m. - (2)
Create a two-dimensional array for the fraction of the primary stream, $Fp\xaf$, with the condition $Fp\xaf=1$ for

*r*≤*r*(radius of the primary inlet), and $Fp\xaf=0$ for_{p}*r*>*r*. The mixture fraction field for coal off gas,_{p}*η*, is zero at the inlet because we assumed there were no coal particle reactions in the burner. - (3)
Obtain density, $\rho \xaf$, from the equilibrium chemistry lookup table using $Fp\xaf$ and assuming the burner is adiabatic.

- (4)Compute the components of the velocity in the flow ($u\xaf$) and the tangential ($v\xaf$ and $w\xaf$) directions, by filtering the mass flux at the inlet,
*ρu*(units of $(kg/m2s)$), using Eq. (6), 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\xaf,\u2009\rho u\xaf$ is divided by $\rho \xaf$. For $v\xaf$ and $w\xaf,\u2009\rho uv\xaf$ and $\rho uw\xaf$ are divided by $\rho u\xaf$.$(\rho u\varphi )i,j\xaf=\u2211n=j*ratio(j+1)*ratio\u2211m=i*ratio(i+1)*ratio(\rho u\varphi )m,nratio2$(6) - (5)
Write the filtered velocities $u\xaf,\u2009v\xaf$, and $w\xaf$ to an input (handoff) file that is read by Arches.

This process is illustrated by the velocity fields shown in Fig. 10. 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 filtered mass flow rate of gas in Arches is equal to the experimental mass flow rate.

To compute the flux of particles at the burner tip, we assumed the particle velocities were equal to the gas velocity and the coal mass flow rate was constant for each cell. We converted these mass weights to particles per cubic meter using a coal density of 1300 kg/m^{3} 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.

### Domain Size and Mesh Resolution.

We performed prototyping runs with the 14.2 m domain using a 0% swirling condition. The gas temperature field from these simulations was relatively constant from 6 to 12 m as the coal was burned out by 6 m. Therefore, we reduced the computational domain to 7 m. In this way, we reduced the size of the computation by 50% and the time that we needed to run the simulation due to the reduced residence time; this reduction in computational load allowed us to perform more runs for the sensitivity analysis. To minimize the effects of the outlet boundary condition, we used a Neumann boundary condition for all the scalars and an outlet boundary condition for velocity and pressure. For the outlet boundary condition, we fixed the pressure at the outlet plane and then performed the pressure projection to enforce continuity across the outlet. Recent analysis with L1500 simulations shows errors of less than 2% in the QOIs are achieved with this methodology.

The shortened geometry of the L1500 with a 15 mm resolution is presented in Fig. 11; 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.

## Input/Uncertainty 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. Both are briefly described next.

### Char Oxidation Model.

*l*using the following equation:

_{2}as oxidizer is Eq. (8). The reaction constant in Arrhenius form is $kl=Ale(\u2212El/RT)$. The two parameters in this equation are $AO2$ and $EO2$. The gasification reaction with CO

_{2}as oxidizer is Eq. (9). As with the oxidation reaction, the reaction constant

*k*has two parameters, $ACO2$ and $ECO2$. The mass transfer coefficient

_{l}*k*was obtained using a Sherwood number correlation with a correction factor as shown in Eq. (10)

_{c}### Ash Deposition Model.

Arches uses the one-dimensional wall boundary condition shown in Eq. (11) where *T _{w}* is the internal (hot face) temperature of the wall,

*T*

_{shell}is the external (cold face) temperature,

*R*is the thermal resistance,

*ε*is the emissivity of the ash deposit, and

*q*

_{incident}is the incident heat flux. For the refractory wall,

*T*

_{shell}is the outside wall temperature and for the cooling tubes, it is the water temperature inside the tubes. In this model,

*T*is solved with a Newton solver every time-step.

_{w}*R*is computed with Eq. (12). 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. The deposit thickness,

*d*

_{deposit}, is computed from a deposition model

*T*is lower than the ash fusion temperature,

_{w}*T*

_{slag}, and

*d*

_{deposit}is computed with Eq. (13). In this equation, the ash deposition velocity ($v\xaf$) is computed from a probability model [28] and

*t*

_{sb}, the time since the last soot blowing event, is an input parameter. In regimes 2 and 3,

*T*is equal to or greater than

_{w}*T*

_{slag}and

*d*

_{deposit}is computed with Eq. (14). If the computed value of

*d*

_{deposit}is greater than zero, the model is in regime 2 and

*T*=

_{w}*T*

_{slag}. If the computed value of

*d*

_{deposit}is less than zero, the model is in regime 3. In this regime,

*d*

_{deposit}= 0 and

*T*is computed from Eq. (11)

_{w}### Input/Uncertainty Map Description.

Model parameters and their associated uncertainty ranges are presented in Table 6. In this table, the overall priority means the relative importance of that parameter on the QOIs. In the ideal case, a sensitivity analysis would be performed to assign priority values. However, running a sensitivity analysis with 11 parameters is too expensive for an LES simulation, even with the coarse mesh (15 mm resolution) case described above. If we were to use two simulations for each dimension and combine all the parameters, a total of 2^{11} = 2048 simulations would be needed. Consequently, we carried out a preselection of parameters by assigning priorities based on prior knowledge. In Table 6, the most important parameters are assigned a priority of 6 and the less important a priority of 3.

Range | |||||
---|---|---|---|---|---|

Parameter | Priority | min | max | Nominal value | Basis/comments |

Char oxidation | |||||

$AO2$ | 6 | 58 | 68 | Pre-exponential factor for oxidation reaction, Eq. (8)$(g\u2009cm\u22122\u2009s\u22121\u2009atm\u22121\u2009O2)$ | |

$EO2$ | 6 | 17 | 20 | Activation energy for oxidation reaction, (kcal/mole) | |

$ACO2$ | 3 | 1390 | Pre-exponential factor for gasification reaction Eq. (9)$(g\u2009cm\u22122\u2009s\u22121\u2009atm\u22121\u2009O2)$ | ||

$ECO2$ | 3 | 53,700 | Activation energy for gasification reaction (kcal/mole) | ||

Ash deposition | |||||

d_{max} | 3 | 1 or 5 | Maximum deposit thickness: 1 cm for cooling tubes, 5 cm for everything else | ||

t_{sb} | 6 | 120,666 | 137,322 | Soot blowing time or time since ash was removed from cooling tubes (s) | |

T_{shell} | 3 | External surface temperature - (1) for L1500 walls, see Table 3, (2) for cooling tubes, see Table 4 | |||

R | 3 | 1.02 | Thermal resistance $(\u2211i=1Nlayer(\Delta xi/ki))$ in the L1500 walls (Wm^{–2} s^{–1}) | ||

T_{slag} | 3 | 1516.5 | Ash fusion temperature for Sufco coal (K) | ||

k_{deposit} | 6 | 0.1 | 1.8 | Thermal conductivity of deposit (Wm^{–1}K^{–1}) | |

ε_{w} | 6 | 0.3 | 1 | Wall emissivity |

Range | |||||
---|---|---|---|---|---|

Parameter | Priority | min | max | Nominal value | Basis/comments |

Char oxidation | |||||

$AO2$ | 6 | 58 | 68 | Pre-exponential factor for oxidation reaction, Eq. (8)$(g\u2009cm\u22122\u2009s\u22121\u2009atm\u22121\u2009O2)$ | |

$EO2$ | 6 | 17 | 20 | Activation energy for oxidation reaction, (kcal/mole) | |

$ACO2$ | 3 | 1390 | Pre-exponential factor for gasification reaction Eq. (9)$(g\u2009cm\u22122\u2009s\u22121\u2009atm\u22121\u2009O2)$ | ||

$ECO2$ | 3 | 53,700 | Activation energy for gasification reaction (kcal/mole) | ||

Ash deposition | |||||

d_{max} | 3 | 1 or 5 | Maximum deposit thickness: 1 cm for cooling tubes, 5 cm for everything else | ||

t_{sb} | 6 | 120,666 | 137,322 | Soot blowing time or time since ash was removed from cooling tubes (s) | |

T_{shell} | 3 | External surface temperature - (1) for L1500 walls, see Table 3, (2) for cooling tubes, see Table 4 | |||

R | 3 | 1.02 | Thermal resistance $(\u2211i=1Nlayer(\Delta xi/ki))$ in the L1500 walls (Wm^{–2} s^{–1}) | ||

T_{slag} | 3 | 1516.5 | Ash fusion temperature for Sufco coal (K) | ||

k_{deposit} | 6 | 0.1 | 1.8 | Thermal conductivity of deposit (Wm^{–1}K^{–1}) | |

ε_{w} | 6 | 0.3 | 1 | Wall emissivity |

In the char oxidation model, the parameters are the kinetic constants for the oxidation reaction (Eq. (8)) and gasification reaction (Eq. (9)). The effect of the gasification reactions in oxy-combustion of pulverized coal has been studied by Hecht et al. [29,30]. While gasification reactions can make a measurable contribution to the carbon consumption, oxidation reactions are the main consumer of carbon. Hecht et al. [30] report literature values for relative rates of CO_{2} gasification of char with respect to oxidation in the range of 1–3 × 10^{−3} at 2000 K. Additionally, in the reaction mechanism used by Erland et al. [31], the pre-exponential factor (Arrhenius expression) for the CO_{2} gasification reaction is at least five-orders of magnitude lower than that of the O_{2} oxidation reaction. Based on this prior knowledge, we assigned a priority of 6 to the oxidation reaction parameters and a priority of 3 to the gasification reaction parameters.

The erosion thickness, *d*_{max}, 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 time frame of the simulation. It is difficult to assign a value to the soot blowing time (*t*_{sb}) 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 (*T*_{shell}) 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 *T*_{shell} had a negligible impact on the internal furnace temperature.

The thermal resistance $(\u2211i=1Nlayer(\Delta 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 is well-characterized, a priority of 3 was assigned to this parameter. The temperature parameter *T*_{slag} was analyzed in an unpublished VUQ study of a 15 MWth oxy-coal boiler performed by our research group. This study concluded that the *T*_{slag} parameter was well represented by the ash fusion temperature. Therefore, a priority of 3 was given to *T*_{slag}, 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 parameter *ε _{w}* was given a priority of 6 because there is no information about this parameter in the L1500.

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

*t*. The other parameters in Table 6 were fixed at the nominal value in the simulations.

_{sb}The uncertainty ranges for the char oxidation parameters $AO2$ and $EO2$ were taken from Smoot and Smith [32], who reported small variations (17–20 kcal/mole) in $EO2$ for three U.S. coals: Wyoming subbituminous, Illinois No. 6 bituminous, and Pittsburgh No. 8 bituminous. For $AO2$, the range of $58\u221268\u2009g\u2009cm\u22122\u2009s\u22121\u2009atm\u22121\u2009O2$ was selected based on the reported value of $60\u2009g\u2009cm\u22122\u2009s\u22121\u2009atm\u22121\u2009O2$ 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 et al. [33] for coal ash and synthetic ash. The *ε _{w}* range was bounded by the maximum and minimum values reported by Wall et al. [34,35] for coal ash. The parameter

*t*

_{sb}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 while the minimum value excluded the first day of experiments when the coal feed rate was lower.

## Sensitivity Analysis

Based on the I/U map (Table 6), there are five active variables. However, a full five-dimensional VUQ study with Arches simulations is too computationally expensive. Therefore, we performed a global sensitivity analysis to reduce the number of dimensions in this study to two.

### Computation of Quantity of Interests From Simulation Data.

Three narrow angle radiometers were used to measure the radiative heat flux in sections 1–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 Eq. (15). The solid angle Ω was computed using Eq. (16). The radiative intensity in each ray, *I _{r}*, was computed with Eq. (17) using the gas temperature,

*T*, the gas absorption coefficient,

*k*, and the mesh resolution, Δ

*x*, obtained from the simulations. Because of the coarseness of the computational mesh, one ray (

*N*= 1), extending from the wall opposite the radiometer to the wall where the radiometer tip was located, was sufficient to account for the narrow view angle of the radiometer [36].

_{r}*I*, we used simulation values for wall temperature,

_{o}*T*, and wall emissivity,

_{w}*ε*; see Eq. (18). More detailed information about the computation of the radiometer heat flux is found in Refs. [36] and [37].

_{w}*Q*

_{removal}) with Eq. (19) using time-averaged fields of

*q*

_{incident}and

*T*. Here,

_{w}*A*corresponds to the surface area of the cell faces that are in contact with the gas,

*T*is the surface temperature of the cooling tube,

_{w}*q*

_{incident}is the incident heat flux, and

*ε*is the emissivity of the tube. We then added values for the whole cooling tube to obtain the total heat removal. We performed this computation with the visualization tool VisIt [38] using data that were computed and saved in Arches

For the wall temperature QOI, we extracted the simulation wall temperature at the thermocouple location.

### Generating Polynomial Chaos Surrogate Models.

We performed the global sensitivity analysis using the Uncertainty Quantification Toolkit (UQTk), a set of C++ tools with a Python interface, developed by Debusschere et al. [14]. UQTk uses the app pce_sens to compute total and main sensitivity indices using a PC surrogate model.

*where*

^{n}*n*is the number of dimensions) were needed for this study. The UQTk app generate_quad generated 32 quadrature points of

*ξ*= ±0.58 with weights of

_{i}*w*= 0.0625 for each dimension. The variable

*ξ*is mapped to physical input space for the Arches simulations using Eq. (20), where

_{i}*a*and

_{i}*b*are the bounds in the uncertainty interval presented in the I/U map (Table 6).

_{i}We ran a set of Arches simulations at the 32 quadrature points. We then computed the radiative heat flux, heat removal, and wall temperatures from the simulations and generated the PC surrogate model for each of these data points.

### Results.

We computed the main and the total sensitivity indices with the UQTk app pce_sens using the coefficients of the PC surrogate model. In Fig. 12, the main sensitivity indices for the radiative heat flux are presented. The two most sensitive parameters with respect to the radiometer measurements are ash deposition parameters, *ε _{w}* and

*k*

_{deposit}. The main sensitivity index for

*t*

_{sb}is low by comparison. The main sensitivity indices for $AO2$ and $EO2$ are also low, meaning these parameters have little influence on the radiometer measurement. Thus, for the radiative heat flux measurement, the five-dimensional (parameter) study can be reduced to a two-dimensional study.

Figure 13 shows the main sensitivity indices for the heat removed by the cooling tubes. The sensitivity indices exhibit a similar behavior to those obtained in the radiometer analysis: the indices for *ε _{w}* and

*k*

_{deposit}are large while that for

*t*

_{sb}is small. The sensitivities of the char oxidation parameters are small as well. As with the radiometer measurements, the study can be reduced from five to two dimensions for the heat removal measurements.

In Fig. 14, the main sensitivity indices for the wall temperature are presented. For all locations, the biggest sensitivity index is for *ε _{w}*. While the sensitivity index for

*k*

_{deposit}is low at location 1 (

*WT*

_{1}), 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 *t*_{sb} (which determines deposit thickness in regime 1 of our deposit model) 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 *T*_{slag} in sections 1 and 2, and for some simulations, wall temperatures exceeded *T*_{slag} for the entire length (7 m) of the simulated reactor as presented in Fig. 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. 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 tube deposits is 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 Fig. 12). The heat removal by the cooling tubes is 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 Fig. 13. In sections 1 and 4,

*HR*

_{1},

*HR*

_{2},

*HR*

_{7}, and

*HR*

_{8}are the only QOIs where the sensitivity indices of

*t*

_{sb}are greater than those of

*k*

_{deposit}. In section 2 (

*HR*

_{3}and

*HR*

_{4}),

*k*

_{deposit}has a bigger sensitivity index than

*ε*. Finally, in section 3 (

_{w}*HR*

_{5}and

*HR*

_{6}),

*k*

_{deposit}still has a strong influence.

## 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. 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. The selection of QOIs requires an understanding of how both the experimental and simulation data were collected. The L1500 furnace was fired with a Utah bituminous coal (Sufco) during a week-long experimental campaign. For the experiments described in this paper, the burner was set at 100% swirl. The simulations performed in this study utilized a handoff strategy—the complex flow through the burner was computed using star-ccm+ and then averaged scalar and velocity fields at the plane of the burner tip were provided as inputs to the Arches simulations of the L1500. We analyze 11 parameters from both the char oxidation and ash deposition models for the I/U map. Using prior knowledge, we reduce the number of active parameters from 11 to 5. After performing a sensitivity analysis, 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 remainder of the VUQ study described in a companion paper.

## 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 et al. for their help with the Uncertainty Quantification Toolkit.

## Nomenclature

*A*=particle surface area (m

^{2})*c*=mixture molar concentration (kmole m

^{–3})- $cO,lg$ =
molar concentration of oxidizer

*l*in the bulk (kmole m^{–3}) *d*_{deposit}=deposit thickness (m)

*d*_{max}=erosion thickness (m)

*d*=_{p}particle diameter (m)

*D*=_{om}mixture averaged diffusion coefficient of oxidizer (m

^{–2}s^{−1})*F*=mixture fraction; see Eq. (5)

*F*=_{p}mixture fraction; see Eq. (4)

*I*=_{o}radiative intensity from a wall (W m

^{–1})*I*=_{r}radiative intensity in each ray (W m

^{–2})*k*=gas absorption coefficient (m

^{–1})*k*=_{c}mass transfer coefficient (m s

^{–1})*k*_{deposit}=deposit thermal conductivity (W m

^{–1}K^{–1})*k*=_{i}thermal conductivity for layer

*i*(W m^{–1}K^{–1})*k*=_{l}reaction rate coefficient for reaction

*l*(m s^{–1})*m*=_{c}coal off gas mass flow rate (kg s

^{–1})*m*=_{n}moments of the Rosin-Rammler distribution

*m*=_{p}primary stream mass flow rate (kg s

^{–1})*m*=_{s}secondary stream mass flow rate (kg s

^{–1})*N*=_{r}number of rays

*q*=radiative heat flux (W m

^{–2})*q*_{incident}=incident radiation (W m

^{–2})*q*_{net}=net heat flux (W m

^{–2})*Q*_{removal}=heat removed by the cooling tubes (W)

*R*=thermal resistance (W m

^{–2}K^{–1})*r*_{H}_{,}=_{l}volumetric reaction rate of char consumed by oxidizer reaction

*l*(kg m^{–3}s^{–1})*r*=_{p}radius of primary inlet (m)

*r*=_{t}total volumetric reaction rate (kg m

^{–2}s^{–1})- Re =
particle Reynolds number

- ratio =
ratio between Arches resolution and STAR-CCM+ resolution

- Sc =
Schmidt number

- Sh =
Sherwood number

*T*=gas temperature (K)

*t*_{sb}=soot blowing time (s)

*T*_{shell}=external surface temperature (K)

*T*_{slag}=slagging temperature of ash (K)

*T*=_{w}wall temperature (K)

*u*=velocity in the flow direction (m s

^{–1})- $v\xaf$ =
ash deposition velocity (kg s

^{–1}) *W*=mixture molecular weight (kg kmole

^{–1})*w*=particle number density (# m

^{–3})*W*=_{H}char molecular weight (kg kmole

^{–1})*w*=_{i}weights

*x*=_{i}particle size

*i*(*μ*m)- Δ
*x*=resolution (m)

- Δ
*x*=_{i}thickness of layer

*i*(m) *ε*=emissivity

*ε*=_{w}wall emissivity

*η*=mixture fraction; see Eq. (3)

*θ*=view angle (deg)

*θ*=_{r}ray polar angle (deg)

*ρ*=gas density (kg m

^{–3})*σ*=Stefan Boltzmann constant, 5.67037 × 10

^{–8}(W m^{–2}K^{–4})*ϕ*=scalar

*ϕ*=_{l}stoichiometric coefficient ratio for species $l\u2009[kmolechar\u2009kmolel\u22121]$

- Ω =
solid angle