{'fraction_non_alphanumeric': 0.035693210634599126, 'fraction_numerical': 0.01239214685571203, 'mean_word_length': 4.532692138004059, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bi-fidelity approximation for uncertainty quantification and sensitivity analysis of irradiated particle-laden turbulence

Efficiently performing predictive studies of irradiated particle-laden turbulent flows has the potential of providing significant contributions towards better understanding and optimizing, for example, concentrated solar power systems. As there are many uncertainties inherent in such flows, uncertainty quantification is fundamental to improve the predictive capabilities of the numerical simulations. For large-scale, multi-physics problems exhibiting high-dimensional uncertainty, characterizing the stochastic solution presents a significant computational challenge as many methods require a large number of high-fidelity solves. This requirement results in the need for a possibly infeasible number of simulations when a typical converged high-fidelity simulation requires intensive computational resources. To reduce the cost of quantifying high-dimensional uncertainties, we investigate the application of a non-intrusive, bi-fidelity approximation to estimate statistics of quantities of interest associated with an irradiated particle-laden turbulent flow. This method relies on exploiting the low-rank structure of the solution to accelerate the stochastic sampling and approximation processes by means of cheaper-to-run, lower fidelity representations. The application of this bi-fidelity approximation results in accurate estimates of the QoI statistics while requiring a small number of high-fidelity model evaluations.


Introduction
The ability to quantitatively characterize and reduce uncertainties, in conjunction with model verification and validation (V&V), plays a fundamental role in increasing the reliability of numerical simulations. These types of studies are commonly encompassed within the field of uncertainty quantification (UQ), which has attracted increasing attention in the modeling and simulation community. In this regard, the Predictive Science Academic Alliance Program (PSAAP) II at Stanford University [1], focuses on advancing the state-ofthe-art in large-scale, predictive simulations of irradiated particle-laden turbulence relevant to concentrated solar power (CSP) systems. To this end, physics-based models are developed and the numerical predictions are validated against data acquired from an in-house experimental apparatus designed to mimic a scaleddown particle-based solar energy receiver, and for which the quantification of uncertainties is of paramount importance. A significant challenge, and the scope of this work in particular, is investigating optimal UQ strategies for this complex, multi-physics flow when many sources of uncertainty are present.

Irradiated Particle-Laden Turbulent Flow
Turbulent flow laden with inertial particles, or droplets, in the presence of thermal radiation is encountered in a wide range of natural phenomena and industrial applications. For instance, it is well established that turbulence-driven particle inhomogeneity plays a fundamental role in determining the rate of droplet coalescence and evaporation in ocean sprays [2] and atmospheric clouds [3]. Another example is found when studying fires, in which turbulence, soot particles, and radiation are strongly interconnected resulting in very complex physical processes [4]. From an industrial point of view, important applications include the atomization of liquid fuels in combustion chambers [5], soot formation in rocket engines [6], and more recently, volumetric particle-based solar receivers for energy harvesting [7].
Even in the simplest configuration, e.g., homogeneous isotropic turbulence, particle-laden turbulent flow is known to exhibit complex interactions between the carrier and dispersed phases in the form of preferential concentration and turbulence modulation [8]. Preferential concentration is the phenomenon by which heavy particles tend to avoid intense vorticity regions and accumulate in regions of high strain rate, while turbulence modulation refers to the alteration of fluid flow characteristics in the near-field region of particle clusters as a result of two-way coupling effects, e.g., enhanced dissipation, kinetic energy transfer, or formation of wakes and vortexes. The physical complexity is further increased by the simple addition of solid walls as turbophoresis [9], i.e., tendency of particles to migrate towards regions of decreasing turbulence levels, becomes an important mechanism for augmenting the inhomogeneity in spatial distribution of the dispersed phase by driving particle accumulation at the walls.
As described above, characterization of particle-laden turbulent flow is a difficult problem; many experimental and numerical research studies have been devoted to this objective over the past decades, see, e.g., [10,11,12]. However, the problem of interest in this work involves, in addition to particle-flow coupling, heat transfer from the particles to the fluid through radiation absorption. The practical application motivating the study of this phenomenon is the improvement of energy harvesting in volumetric particle-based solar receivers. At present, most CSP technologies use surface-based collectors to convert the incident solar radiation into thermal energy. In this type of system, the energy is transferred to the working fluid downstream of the collection point via heat exchangers, typically resulting in large conversion losses at high temperatures. By contrast, volumetric solar receivers continuously transfer the energy absorbed by particles directly to the operating fluid as they are convected through an environment exposed to thermal radiation. This innovative technology is expected to increase the performance of CSP plants by avoiding the necessity of heat-exchanging stages, while requiring significantly high radiation-to-fluid energy transfer ratios. This requirement imposes a very complex design constraint as the physical mechanisms governing irradiated particle-laden turbulent flow are still not fully comprehended, and therefore is a topic of intense research. For example, Zamansky et al. [13] studied the interaction between radiation, particles and buoyancy in homogeneous isotropic turbulence (HIT) by means of point particle direct numerical simulation (PP-DNS) with the assumption of a constant radiation intensity in an optically thin environment. The study revealed how non-uniformities of the gas-particle mixture resulted in inhomogeneities in heat absorption, and therefore in spatial temperature variations that induced large-scale fluid motion by local gas expansion and buoyancy. The resulting baroclinic vorticity production triggered a feedback loop that generated larger non-uniformities by centrifuging the inertial particles leading to new non-uniformities in heat absorption. In the context of buoyancy effects in HIT, Frankel et al. [14] investigated the impact of heating on the settling of particles, finding that, contrary to non-heated particles which enhance small-and large-scale turbulent motions when their settling velocity is sufficiently high compared to the Kolmogorov velocity, the heating of particles resulted in a significant reduction of the mean settling velocity caused by rising buoyant plumes in the vicinity of particle clusters and which affected all scales of turbulence. However, in regimes relevant to particle-based solar receivers, the background turbulence dominates over the buoyancy effects as the Froude number of the systems are typically very large. This problem was studied by Pouransari and Mani [15] by considering a convective flow of gas and particles subject to radiation and in the absence of buoyancy effects. In their study, it is shown how different particle concentration patterns, resulting from different Stokes numbers, have a strong influence on the heat transfer rates between particles and gas. A similar flow configuration, but considering polydispersity of particles, has been recently analyzed by Rahmani et al. [16], finding also that the effective heat transfer rate between the two phases is significantly impacted by particle clustering in addition to polydispersity.
presented in [45], which can be used to determine if a given pair of low-and high-fidelity models will lead to an accurate BF approximation. In the present study, we adopt the BF approach and error estimate of [45] to illustrate their efficacy in UQ of a CSP systems as an instance of a complex, multi-physics engineering application.

Objectives and Organization of the Work
The system studied in this work is based on a small-scale apparatus [46] designed to enable the reproduction of the physical mechanisms encountered in volumetric particle-based solar receivers. Physics-based modeling of irradiated particle-laden turbulent flow -as detailed in Section 2 -and its numerical investigation and validation against data obtained from the experimental apparatus -description given in Section 3are difficult tasks that intrinsically require several model assumptions, selection of coefficient and parameter values, and characterization of initial and boundary conditions. These steps, even if performed carefully, result in potential sources of uncertainty that can impact the quantities of interest (QoIs). Examples of such uncertainties include the incomplete description of particle diameters [16] and thermal radiative properties [47], variability of the incident radiation and its complex interaction with boundaries, and structural uncertainty inherent in the approximations utilized, e.g., turbulence modeling [48]. In addition to the large number of uncertainty sources, accurate predictions of the complex interaction of particle-laden turbulent flow with radiative heat transfer demand the utilization of expensive HF numerical simulations. As an example, the cost of a medium-scale HF calculation of this problem requires approximately 500k core-hours per sample on the Mira supercomputer (ALCF) [49]. Therefore, if brute-force UQ techniques, e.g., MC simulation with O(10 3 ) samples, are to be performed, the total cost is of the order of 500M core-hours, thus motivating the need for accelerated UQ strategies. In this regard, the objective of this work is to investigate the BF approximation UQ strategy on large-scale, multi-physics applications based on the PSAAP II solar receiver. In particular, it can be shown that this BF approximation provides accurate estimates of the QoI statistics, while maintaining a reduced cost similar to that of LF models for simulations with O(10 3 ) samples.
The paper is organized as follows. In Section 2, the physical models and numerical method utilized to simulate irradiated particle-laden turbulent flow are described. The particle-based solar receiver studied is detailed next, in Section 3, in terms of computational setup, input uncertainties, and QoIs considered. The BF approximation strategy is presented in Section 4, as well as a brief discussion of the associated error bound. The performance of this BF approximation, with regards to both accuracy and cost is investigated in Section 5. From these results, comparisons are made between this approximation and alternative LF models. Finally, the work is concluded and future directions are proposed in Section 6.

Physics Modeling and Numerical Method
The problem under study comprises the physical mechanisms depicted in Figure 1. The instantaneous snapshot, extracted from a single HF simulation, corresponds to a slice of the normalized temperature increment, (T − T 0 )/T 0 , of fluid and particles in the x-y plane along the streamwise direction. The fluidparticle mixture flows from a development square duct, with statistically steady-state distribution of particles (clustering + turbophoresis) at temperature T 0 = 300 K, into a windowed section with the same dimensions but transversely illuminated by lasers from below. The result is that the micron-sized, nickel particles absorb thermal radiation as they are advected through the duct; predominantly near the walls where they accumulate due to turbophoretic effects. The thermal energy is then transferred to the fluid in the vicinity of particles by thermal exchange, and subsequently mixed into the bulk of the fluid by diffusion processes interacting with the cascade of turbulent scales. The equations describing this type of flow are continuity, Navier-Stokes in the low-Mach-number limit, conservation of energy assuming ideal-gas behavior, Lagrangian particle transport, and radiative heat transfer. The modeling of these three physical processes -turbulent flow, particle transport, and radiative heat transfer -and their couplings, are sequentially described in the subsections below. Figure 1: Irradiated particle-laden turbulent flow in a square duct (instantaneous, streamwise x-y plane snapshot). The quantity represented is the normalized temperature increment (T 0 = 300 K) of fluid and particles. It is interesting to note (i) the presence of elongated particle clusters resulting from preferential concentration, (ii) the accumulation of particles at the walls by turbophoretic effects, (iii) the absorption of radiation by particles (illuminated from below), and (iv) the transfer of energy from the particles to the fluid and the subsequent thermal mixing enhanced by turbulence.

Variable-Density Turbulent Flow
The volumetric particle-based solar receiver operates at atmospheric pressure conditions in which air, the carrier fluid, is assumed to follow the ideal-gas equation of state (EoS), P th = ρ g R air T g , where P th is the thermodynamic pressure, ρ g is the density, R air is the specific gas constant for air, and T g is the temperature. As indicated by the EoS, density varies with temperature. However, the Mach number of the flow Ma = u/c, with u the local flow velocity and c the speed of sound of the medium, is small (Ma < 0.03) for the range of velocities and temperatures considered. Therefore, the low-Mach-number approximation is utilized to separate the hydrodynamic part, p P th , from the total pressure, P tot = P th + p. This decomposition results in the following equations of fluid motion where u g is the gas velocity, ρ 0 is an ambient reference density, I is the identity matrix, g is the gravitational acceleration, and µ g and λ g are the dynamic viscosity [50] and thermal conductivity [51], respectively. Additionally, C v,g and C p,g are, respectively, the isochoric and isobaric specific heat capacities, and f T W C and S T W C are two-way coupling terms representing the effect of particles on the fluid and, respectively, approximated as point sources in the forms Here, m p = ρ p πd 3 p /6 and v p are the particle mass and velocity, respectively, u p is the gas velocity at the particle location, τ p = ρ p d 2 p /(18µ g ) is the particle relaxation time, d p is the particle diameter, and δ (x − x p ) is the Dirac delta function concentrated at the particle location x p . The expression for the fluid-particle convection coefficient is h = N uλ f /d p with N u the particle Nusselt number. In the problem studied, particles are assumed to be isothermal as the Biot number is small, Bi = hd p /λ p 1.

Lagrangian Particle Transport
The carrier fluid is transparent to the incident radiation. Hence, micron-sized nickel particles, i.e., the dispersed phase, are seeded into the gas to generate a non-transparent gas-particle mixture that absorbs and transfers the incident radiation from the particles to the gas phase. The diameters of the particles are several orders of magnitude smaller than the smallest significant (Kolmogorov) turbulent scale η, and the density ratio between particles and gas is ρ p /ρ g 1. As a result, particles are modeled following a Lagrangian point-particle approach with Stokes' drag as the most important force [52]. Their description in terms of position, velocity and temperature is given by where C v,p is the particle specific isochoric heat capacity, ω = Q s / (Q a + Q s ) is the scattering albedo with Q a and Q s the absorption and scattering efficiencies, respectively, I is the radiation intensity, σ is the Stefan-Boltzmann constant, and dΩ = sin θdθdφ is the differential solid angle. In the conservation equation for particle temperature, (8), the first term on the right-hand-side accounts for the amount of radiation absorbed by a particle, while the second term represents the heat transferred to its surrounding fluid. In the point-particle approximation, particle-wall and particle-particle interactions are typically described by one-dimensional collision models based on the balance of total momentum and energy. In the case of collisions involving two objects, A and B, the velocities after impact, v A and v B , are given by where u A and u B are the velocities of the objects before impact, m A and m B are the mass of the objects, and 0 ≤ C R ≤ 1 is the restitution coefficient. The limits of C R correspond to the cases in which the objects coalesce at impact (0, perfectly inelastic collision) and rebound with the same relative speed as before impact (1, perfectly inelastic collision); intermediate values represent inelastic collisions in which kinetic energy is dissipated. The above equations simplify to v A = −C R u A and v B = 0 when object B is a static wall.

Radiative Heat Transfer
In the problem under consideration, the flow and particle timescales are orders of magnitude larger than the radiation timescale, which is related to the speed of light. As a consequence, it can be assumed that the radiation field changes instantaneously with respect to temperature and particle distributions; i.e., radiation field is quasi-steady. Under this assumption, and considering that air is transparent at all wavelengths and that absorption and scattering are determined solely by the presence of particles and solid boundaries, the radiative heat transfer equation becomeŝ whereŝ is the direction vector, σ e = σ a + σ s is the total extinction coefficient with σ a and σ s the absorption and scattering coefficients, respectively, and Φ is the scattering phase function that describes the directional distribution of scattered radiation. The total extinction coefficient can also be defined in terms of absorption and scattering efficiencies as σ e = (Q a + Q s ) πd 2 p n p /4 with n p the local number density of particles. Moreover, assuming gray radiation Q a +Q s ≈ 1, which leads to ω ≈ Q s (see (8)), and as a result σ a ≈ Q a πd 2 p n p /4 and σ s ≈ Q s πd 2 p n p /4. Figure 2: Computational setup of the PSAAP II volumetric particle-based solar energy receiver. An isothermal periodic section (left domain) is utilized to generate fully developed particle-laden turbulent flow, which is used as inflow conditions for the second section (right domain) where the gas-particle mixture is irradiated perpendicularly to the flow direction from one the sides.

Numerical Method
The equations of fluid motion, (1)-(3), are solved following an Eulerian finite-volume discretization implemented in an in-house solver that is second-order accurate in space and suitable to non-uniform, Cartesian meshes. A fourth-order Runge-Kutta scheme is used for integrating the equations in time, together with a fractional-step method for imposing conservation of mass [53]. Integration in time of the Lagrangian position, velocity, and temperature of particles, (6)- (8), is fully coupled with the advancement of the flow equations to ensure fourth-order accuracy. The transfer of radiative heat, (11), is calculated by means of an in-house discrete ordinates method (DOM) interfaced to the flow solver via an Eulerian representation of the particles distribution.

Computational Setup and Physical Parameters
Numerical simulations of the volumetric particle-based solar receiver are performed on the computational setup depicted in Fig. 2. Two square duct domains, with dimensions 1.7L×W ×W (L = 0.16 m, W = 0.04 m) in the streamwise (x-axis) and wall-normal directions (y-and z-axis), are utilized to mimic the development and radiated sections of the experimental apparatus. The development section (left domain) is an isothermal, T 0 , periodic particle-laden turbulent flow generator that provides inlet conditions for the inflow-outflow radiated section (right domain). The solid boundaries of the development section (y-and z-sides) are considered smooth, no-slip, adiabatic walls. Regarding the radiated section, the same boundary conditions are imposed except for the radiated region in which the y-and z-boundaries are modeled as non-adiabatic walls accounting for heat fluxes due to the radiation energy absorbed by the glass windows.
The bulk Reynolds number of the gas phase at the development section is Re b = ρ g u b L/µ g = 20000, with u b the gas bulk velocity. As characterized from the experiments, the particle size distribution is approximated by 5 different classes (uniformly sampled) with Kolmogorov Stokes numbers St η = 6, 8, 10, 12, 17 (diameters listed in Table 1) and with a total mass loading ratio (MLR) of MLR = n p m p /ρ g ≈ 20%. Detailed values of the development section flow conditions and material properties are listed in Table 1. The gas-particle mixture, as depicted in Fig. 3, is volumetrically irradiated through an L × W glass window starting at ∆x = 0.1L from the beginning of the radiated section. The radiation source consists of an array of diodes mounted on a vertical support placed ∆y ≈ 3W from the radiated window and aligned with the streamwise direction of the flow. The diodes generate a total power of P ≈ 1 kW approximately uniform within a 18 • cone angle. Detailed configuration of the illuminated computational section (right). The test segment is 1.7L long and W width, with a radiated volume of L × W × W starting at 0.1L from the beginning of the section, and with a probe perpendicular to the flow and located 0.3L downstream.
The intervals of the input variables listed in Table 2 have been carefully characterized on the basis of information provided by the team responsible for conducting the experiments, and by taking into consideration results and conclusions extracted from published studies. In particular, the intervals of the particle restitution coefficients follow the trend observed in experimental investigations by Yang & Hunt [54] in which C R increases with Stokes number. The expression for Stokes' drag force correction and its coefficient interval is based on the theoretical work by Brenner [55]. The particle Nusselt number range is extracted from numerical studies of heated particles performed by Ganguli & Lele [56,57]. The intervals for particle absorption and scattering efficiencies are obtained from Mie scattering theory and take into account sensitivity to shape deformation as investigated by Farbar et al. [58]. The long development section of the experiment is modeled in the simulations by means of a periodic domain. This approach reduces the computational cost noticeably, but at the expense of not having direct control of the mass loading ratio. Consequently, the interval for this quantity is characterized based on preliminary comparisons with experimental data for a periodic particle-laden duct flow case without thermal radiation at steady-state conditions. Similarly, the heat fluxes from the walls to the fluid are not fully resolved in the simulations and are instead directly modeled as boundary conditions; their calculation would require the solution of a complex conjugate convective heat transfer problem. In order to properly represent the uncertainty associated with this approximation, comparisons against the experiment were performed for an irradiated turbulent flow case without particles. This was done in a manner that the thermal response of the system was fully driven by the fluxes from the walls to the fluid and the intervals for these uncertainties could be obtained. Finally, the variability of the incident radiation is based on the uncertainty in emission intensity of the laser diodes reported by the manufacturer and the uncertainty in absorptivity of the glass walls.
The performance of the BF approximation is focused on thermal QoIs at the probe location. As detailed in Fig. 3, the probe is located ∆x = 0.3L downstream from the radiated perimeter, and is perpendicular to the flow direction along the y-axis at z = W/2. Of particular interest in this study is the time-averaged, normalized increment of gas temperature along the y-axis profile, i.e., Q = ( T − T 0 )/T 0 = ∆T /T 0 , and the average heat flux over the plane at the probe location, i.e., Q = C p,g ρ g u g ∆T dS.

Simulation Strategy
In the experimental apparatus, fully developed conditions are achieved using a long duct (≈ 7m) with an aspect ratio of order hundred. To reduce the computational cost of simulating a long duct with inflowoutflow boundary conditions, the computational setup is divided in two domains as described in Section 3.1. First, randomly distributed particles are seeded into the development section with an initial fully developed turbulent velocity field. This system is then evolved in time for 20 flow through times (FTTs), defined as FTT = L/u b , to achieve fully developed turbulent particle-laden flow conditions as in the experiments. After 20 FTTs, the instantaneous Eulerian and Lagrangian solutions are copied into the second domain, radiative illumination is activated in the radiated region, and the total system (two domains) is integrated in time for 5 additional FTTs intended to flush the thermal transient (1 FTT) and collect statistics (4 FTTs). This procedure is repeated independently for each realization.

Description of the High-and Low-Fidelity Models
Three model fidelities have been designed to perform the UQ study: one HF model and two LF representations, denoted LF1 and LF2. All three models use the same description of uncertainty as described in Section 3.2. The HF model corresponds to a point particle direct numerical simulation (PP-DNS) with sufficient resolution (≈ 55M cells/section) to capture all the significant (integral to Kolmogorov) turbulent scales, while approximating the particles as Lagrangian points (≈ 15M particles/section) with nonzero mass. The flow grid is uniform in the streamwise direction with spacings in wall units equal to ∆x + ≈ 12, while stretched in the wall-normal directions with the first grid point at y + , z + ≈ 0.5 and with resolutions in the range 0.5 < ∆y + , ∆z + < 6. The radiative heat transfer equation is solved on a uniform DOM mesh of 270 × 160 × 160 grid points (≈ 7M cells) with 350 quadrature points (discrete angles) per element.
Based on this HF model, two LF models have been constructed by coarsening the Eulerian and Lagrangian resolutions, resulting in the LF1 and LF2 representations that are, respectively, ≈ 170× and ≈ 1300× cheaper to simulate (per sample) than the HF model. The flow and radiation meshes, and quadrature points are uniformly coarsened in each direction by a factor of 5 (LF1) and 10 (LF2), which additionally allows for larger time steps; a schematic of the different fidelity models for the turbulent flow phase is shown in Fig. 4. The discrete phase is coarsened by reducing the number of point particles in the calculation by 5× (LF1) and 10× (LF2). To preserve the dimensionless parameters of the problem, this is efficiently accomplished by grouping physical particles into parcels representing their total effect [59], viz. surrogate particles [60]. In the case of uniform parcels, the evolution of the N s surrogate particles can be described with the same set of Lagrangian equations used for the N p physical particles, (6)- (8), with the only modification of multiplying the physics coupling terms by the ratio W s = N p /N s as σ a ≈ W s Q a πd 2 p n p /4 and σ s ≈ W s Q s πd 2 p n p /4.

Construction of Bi-fidelity Approximation
The BF approximation of this work follows the approach of [43,45,44], which seeks to identify a low-rank representation of the HF QoI assisted by the realizations of its LF counterpart. While the QoIs in this work are scalar-valued, this BF approximation relies on vector-valued quantities u that describe the scalar-valued QoI. For example, when the QoI is the time-and spatial-averaged gas temperature along the probe line (see Fig. 3), the elements of u may be the time-averaged estimates of the gas temperature along this profile. The construction of the BF approximation is completed in four main steps, as illustrated in Fig. 5. In the first step, N independent samples of the inputs ξ, each denoted by ξ (i) , are generated according to their joint probability density function (PDF) and used to form N LF realizations of u, each referred to as u (i) L ∈ R m (see Fig. 5, orange boxes, top row). Here, m is the number of (spatial or temporal) degrees of freedom describing u L . These realizations are then organized in an m × N LF data matrix U L such that Typically N is large, as obtaining many LF samples is computationally feasible, and the value of m corresponds to the spatial degrees of freedom of the numerical solver. In the second step, we seek to identify a subset of realizations of u L of size r (and the associated input samples of ξ) that form a basis for the range space of U L . To this end, we perform a rank r N factorization of U L via matrix interpolative decomposition (MID) [61,62,63]. Specifically, the rank r approximationÛ L is obtained via column-pivoted (truncated) QR factorization where Q ∈ R m×r has r orthogonal columns, R 11 ∈ R r×r is an upper triangular matrix, R 12 ∈ R r×N −r , P ∈ R N ×N is a permutation matrix, and † indicates the pseudoinverse. It is straight-forward to show that the product QR 11 in (16) is equal to the left r columns of U L , denoted by and referred to as the column skeleton of U L . In (17), the column indices i k , k = 1, . . . , r, are identified from the permutation matrix P in (15) and indicate input samples ξ (i k ) corresponding to the LF realizations u (i k ) L (see Fig. 5, dark orange boxes, top row). From (16), setting C L = [I R † 11 R 12 ]P to be the LF coefficient matrix, the rank r MID approximation of U L is given aŝ indicating that U c L , or equivalently the LF realizations u L , act as a reduced basis for the range of the LF matrix U L .
In the third step, the HF model is simulated at the inputs ξ (i k ) corresponding to the LF reduced basis (see Fig. 5, dark blue boxes, bottom row). This results in the associated HF column skeleton Note that U c H is M × r, where M ≥ m as it corresponds to the number of HF spatial degrees of freedom. For HF models with a finer mesh resolution than the LF model it follows that M > m. However, as is the case in this work, data values may be extracted at equivalent coordinates by using extraction points independent of the grid tessellation, resulting in M = m. In the final step the BF approximation is formed by taking the product of the HF column skeleton in (19) and the LF coefficient matrix where the i-th column ofÛ H , denotedû H , is the BF approximation to u H at ξ (i) (see Fig. 5, blue boxes, middle row). Once formed, approximate QoI realizations are calculated directly from the BF realizationŝ u With regards to computational cost, the BF approximation requires N LF simulations, performing rankrevealing QR on U L with O(rmN ) floating point operations for basis identification and the coefficient matrix computation, and only r HF simulations to form U c H . Typically, obtaining HF solutions is a bottleneck, and thus limiting the number of HF simulations to small r is of great value and a fundamental component of this approximation.
Remark 1. By construction of MID, the columns ofÛ L with indices i k , k = 1, . . . , r, identified by the permutation matrix P in (15), are exactly the same as the corresponding columns of U L . Stated differently, an r × r sub-matrix of C L is a permutation of the identity matrix. This, therefore, implies that the columns of the BF matrixÛ H with indices i k (identified from MID of U L ) are the same as the r HF realizations forming the columns skeleton U c H .

Remark 2.
While the optimal rank for this BF approximation is not known a priori, the rank of the LF data matrix U L may indicate a good selection range. Theoretical results of [61,62,63] show the error of the rank r LF approximation in (18) is bounded by a scaling of the (r + 1)-th largest singular value of U L . This suggests that the decay of singular values of U L , specifically, where a significant drop occurs, to be a worthwhile consideration for the value of r.

Theoretical Error Estimation of Bi-Fidelity Approximation
This section briefly presents the theoretical results of [45] that provide a practical error bound on the spectral norm of U H −Û H . This error bound may be used with small cost to assess the suitability of a given pair of low-and high-fidelity models to produce an accurate BF approximation. Instances of successful BF approximation viaÛ H in (20) have been reported in a number of recent studies [41,42,43,44,45].
Following [45], the BF estimateÛ H is accurate as long as there exists an M × m matrix T , with bounded norm, such that U H ≈ T U L . For a given pair of low-and high-fidelity problems, the existence of such a mapping is not always guaranteed. However, [45] shows that such a mapping can be constructed when is small enough and that U L (similarly U H ) has fast decaying singular values. In (21), λ max (·) denotes the largest eigenvalue of a matrix, and U T H U H and U T L U L are the Gramians of the HF and LF matrices, respectively. The following theorem from [45] provides a bound on the BF approximation error. In particular we note the dependence of the error bound on τ and (τ ).
Theorem 4.1 (Theorem 1 of [45]). LetÛ H be the rank r BF approximation, as in (20), to the HF data matrix U H . LetÛ L be the rank r MID approximation, given in (18), to the LF data matrix U L , where C L is the corresponding coefficient matrix. For (τ ), as defined in (21), and · the spectral norm, it follows that the BF error may be bounded as where σ k and σ k+1 are the k-th and (k + 1)-th largest singular values of U L , respectively.
Evaluating the bound in (22) requires the computation of (τ ) for multiple values of τ . Following (21), this in turn requires access to the entire HF data matrix U H , which is not possible as U H is never generated in practice. As an alternative, estimates of (τ ) may be calculated using a subset of R N HF and LF samples viaˆ where the superscript R indicates the number of columns of U H and the corresponding columns of U L used to set U R H and U R L , respectively. Stated differently, estimates of (τ ) may be obtained using R HF samples, instead of N . To evaluate the remainder of the bound in (22), MID is applied to the LF data matrix to obtain values for U L −Û L and C L . Combining these estimates and minimizing over identified values of (τ,ˆ (τ )) and the singular values of U L results in an approximate bound. The numerical results of [45] show empirically that R slightly larger than the approximation rank r is sufficient to estimate the optimal pair (τ, (τ )). This therefore suggests the efficacy of (22) in estimating the BF approximation error or in suitability of a given pair of low-and high-fidelity models for BF modeling.

Using Bi-Fidelity Approximation to Estimate QoI Statistics
Our discussion so far has focused on the construction of N BF estimates of a vector-valued QoI using the corresponding N LF realizations, along with r selected HF realizations. We next turn our attention to how these N BF estimates may be utilized for the purpose of UQ or sensitivity analysis.
When the BF approximation achieves the desired accuracy, standard methods such as MC sampling, stochastic collocation, or sparse PCE may be employed onÛ H (more specifically, the columns ofÛ H ) instead of U H to estimate the moments and PDF of the QoI, and perform global sensitivity analysis. When the dimension d of the random inputs is not high, methods such as sparse PCE or stochastic collocation may be employed. Otherwise, MC sampling methods are preferred. Achieving small sampling errors, requires having access to large enough number of BF -and therefore LF -samples N . On the other hand, when the BF approximation does not meet the accuracy requirements but leads to estimates well correlated with the HF data, the BF approximation may serve as a control variate to MC in a single-level [26] or in a multilevel setting as in [36].
In the numerical results of Section 5, we use the 1 -minimization approach of [21,64] to build sparse PCEs of the QoIs as approximate maps between the random inputs ξ and the QoIs. The resulting PCEs are, in turn, used to estimate the statistics, here, histogram, of the QoIs and perform global sensitivity analysis. The QoI statistics can be either computed directly via the PC coefficients, e.g., for the mean and variance, or by sampling the PCE itself in a MC fashion, e.g., for histogram. We follow the latter approach in the experiments of Section 5 to generate histograms of the QoIs. For global sensitivity analysis, we perform variance decomposition to compute the so-called Sobol' indices [65], which provide a means to quantitatively describe the importance of input parameters, by calculating each parameter's contribution to the total variance of the output QoI. Following the work of [66], we compute the Sobol' indices directly from the PCE coefficients.

Numerical Results of Bi-Fidelity Approximation
To investigate the performance of the BF approximation, 256 LF2, 128 LF1, and 26 HF simulations were performed, such that the 26 HF simulations correspond to the first 26 simulations of LF1 and LF2, and the 128 LF1 simulations correspond to the first 128 of the LF2 simulations. From the two LF models, two BF approximations are formed: bi-fidelity 1 (BF1) approximation and bi-fidelity 2 (BF2) approximation. The BF1 approximation is formed from N = 128 LF1 samples and r HF samples, and the BF2 approximation is formed from N = 256 LF2 samples and r HF samples. The number of HF simulations is left as r, as the selection of this value will be discussed in the following results.
The motivation of this approximation is to form a BF model that accurately predicts the HF data, and the goal of these results is to investigate whether or not there is an improvement over the performance of the the LF models. For these tests, two primary time-averaged, i.e., · , thermal QoIs are considered: (i) heat flux through the plane at the probe location, Q = C p,g ρ g u g ∆T dS, (∆x = 0.3L downstream from the radiated perimeter) normalized by the mean QoI estimated from the HF data, and (ii) normalized temperature increment, ∆T /T 0 = ( T − T 0 )/T 0 , values along the profile at the probe location (∆x = 0.3L downstream from the radiated perimeter at z = 0.5W ), where focus is placed on spatially-averaged ∆T /T 0 , and point estimates at y/W = 0.5, y/W = 0.1, and y/W = 0.05.
For each QoI, five primary tasks are considered: (i) BF rank identification as to best optimize accuracy and computational cost of the approximation, (ii) BF error bound estimation, to verify the accuracy of the approximation for a fixed rank, (iii) QoI approximation via available data, (iv) estimation of statistics via sparse PCE of LF and BF models, and (v) cost analysis to compare approximate core-hours needed to obtain converged simulations from each model. For a given matrix U of LF or BF data, we report the relative spectral error, with respect to its HF counterpart U H , defined as For scalar QoIs, we report the relative 2 (root mean-square error), where Q (i) is a LF or BF simulated QoI and Q (i) H is the corresponding HF realization.

QoI #1: Heat Flux Through the ∆x = 0.3L Plane
To estimate the statistics of the heat flux QoI, the BF approximation is formed from realizations u of the heat flux values over the entire ∆x = 0.3L plane. For ease of comparison, these values are extracted as a 500 × 500 uniform grid of points for both HF and LF simulations (thus here m = M ). The scalar-valued QoI estimates Q are taken to be the heat flux determined from all of the elements of u, normalized by the mean QoI estimated from the HF data.
When forming the BF approximation, the first task is to identify the approximation rank. This is critical as it will dictate both the computational cost -corresponding to r HF realizations -and accuracy of the approximation. As a reduced cost relies on minimal number of HF samples, the selected rank must be as small as possible. However, excessively small approximation ranks may lead to inaccurate BF solutions. To aid in rank selection, consider the results of Fig. 6. Fig. 6 (a) provides the decay of singular values of LF and HF matrices. The magnitudes of these singular values decay most rapidly within the first six indices, indicating that six realizations of heat flux capture most of the information of the LF and HF data. Fig. 6 (b) displays the error bound estimate of the BF approximation matrix as a function of rank r. Note, the number of simulations to calculateˆ (τ ) is set to R = r + 2; a more thorough assessment of the error bound will be discussed shortly. The error bound estimate levels out for a conservative value of r = 6, which is the BF approximation rank we use for the rest of the results. For comparison, the LF1 and LF2 spectral errors are provided (see (24)), indicating that very few HF samples are needed to observe an improvement over the LF models. The presentation of the corresponding cost analysis is postponed until Section 5.3. While the previous results indicate that the rank r = 6 BF approximation more accurately describes the HF data, the calculated error must be verified by more thoroughly investigating the theoretical error bound estimates from Section 4.2. Specifically, the bound must be estimated for multiple values of R, the number of samples used to estimateˆ (τ ). Fig. 7 provides the error bound estimates as a function of R for the BF1 and BF2 models. Each point represents an error bound estimate calculated from R random columns of U H (out of 26 total columns) and the corresponding columns of U L . The solid line represents the average value of the points at each value of R. Numerical results of [45] suggest a value of R ≈ 2r will provide a true error bound. With rank r = 6, these results estimate the error bound to be 0.12 for both BF models. Recall from Fig. 6 (b) that this error bound estimate is smaller than that of either LF model, indicating that improvement in accuracy may be estimated without knowledge of the true BF error. To compare the performances of the LF and BF approximations with regards to the heat flux QoI, consider the results of Fig. 8, where 17 simulated values of the heat flux using the HF, LF and BF models are provided. Given the discussion of Remark 1, we exclude the heat flux values corresponding to the r HF realizations used to form the BF approximation. For each simulated value, the BF approximations are significantly more accurate than their LF counterparts. With regards to the relative 2 error, the BF QoIs are about 7× more accurate than the LF QoIs, with both BF1 and BF2 having errors of 0.03, LF1 an error of 0.20, and LF2 an error of 0.22. To estimate the moments and PDF, as well as perform sensitivity analysis, a sparse PCE surrogate is formed from available data, as discussed in Section 4.3. Table 3 provides the QoI mean and coefficient of variation (CoV) values determined from the approximate coefficients of each LF and BF surrogate model, as well as the relative error between each QoI mean and the HF QoI mean. Note the HF QoI mean is calculated directly from available data. These results show that the BF1 mean is about 100× more accurate than the LF mean values and the BF2 mean is about 10× more accurate than the LF mean values, where the BF1 model predicts within 0.15% of the HF QoI mean. The LF models, on the other hand, predict to within only 20 − 22% of the HF QoI mean. The CoV is provided for all models as well. The QoI variation for each fidelity model is comparable to that of the HF model, with that of BF1 most closely describing the HF data. Since the BF2 CoV is of the same order of magnitude as the respective error, it cannot be completely relied on as an estimate for the HF CoV. The BF1 error, on the other hand, is significantly smaller than the CoV, and thus the CoV estimate can be trusted, indicating that the BF1 approximation most closely represents the HF data.  To estimate the QoI PDFs, histograms of the LF and BF data are generated from 100, 000 sparse PCE samples. Fig. 9 provides these normalized histograms for the heat flux QoI to compare with available HF data. Fig. 9 (a) provides data derived from the LF1 and BF1 models, while Fig. 9 (b) provides the histograms from LF2 and BF2 data. We observe that the BF histograms more closely follow the histogram of the HF data. For the final set of results the total Sobol' indices are calculated from the PCE coefficients of the LF and BF heat flux samples. The relative contribution of each parameter to the variance of the estimated heat flux is displayed in Fig. 10 (a)-(d). The Sobol' indices determined from all four sets of PCE coefficients clearly indicate that the inputs representing the uncertainty in the heat fluxes from the radiated wall (ξ 12 ) and its opposite wall (ξ 13 ) contribute the most to the heat flux variation. Comparatively, the remaining uncertain parameters are best distinguished from the BF1 results of Fig. 10 (b) based on its accurate CoV estimate. In particular, the BF1 result indicates that inputs ξ i with i = 4 − 9, 11 are next in terms of contribution to the total QoI variance, with about 1% of the variance, and ξ i with i = 1 − 3, 10 have the least important contribution, with < 0.1% of the total variance. heat flux data. Starting with parameter ξ 1 at the top position, importance of each ξ i is provided in counterclockwise order with respect to increasing i, with corresponding description provided in Table 2. As determined from all four model surrogates, heat flux from the radiated wall (ξ 12 ) and the opposite wall (ξ 13 ) to the fluid contribute the most to the heat flux variability.

Model
As the results of this section have shown, both BF approximations are significantly more accurate than either LF approximation, with the sparse PCE surrogate of the BF1 approximation most accurately estimating the HF heat flux data. A discussion of cost comparisons between the five models may be found in Section 5.3.

QoI # 2: Spatially-Averaged and Point Estimates of ∆T /T 0 Along Profile at Probe Location
Next, time-averaged ∆T /T 0 = ( T − T 0 )/T 0 along the profile at the probe location is considered. The BF approximation employs realizations u of the change in temperature ∆T /T 0 along this profile. After investigating the estimation of the full temperature profile for all models, focus is placed on estimating the spatial mean of ∆T /T 0 along the profile, as well as ∆T /T 0 at three points along the profile, namely, y/W = 0.5, y/W = 0.1, and y/W = 0.05 at the probe location.
Similar to the heat flux QoI, an optimal rank of the BF approximation must be first identified. Fig. 11 (a) shows the decay of the singular values of the LF and HF matrices, indicating that an approximation rank r ≥ 4 accurately represents the LF1 and LF2 data. Fig. 11 (b) displays the BF approximation error bound and the calculated relative spectral errors (see (24)) for both LF data as a function of rank r. For these error bound estimates, R = r + 2 simulations are used to estimate the values ofˆ (τ ). As in QoI #1 above, the error bound estimate for BF1 levels out for the conservative value of r = 6 (notice the larger value of the error estimate for r = 5). We therefore set the BF1 approximation rank to r = 6. While the BF2 bound suggests a rank of r = 5, for the interest of a simpler presentation, we chose r = 6 for BF2 rank as well.  suggest larger values of R -compared to R = r + 2 used in the results of this section -lead to more accurate BF error estimates. In both cases for R = 15 the estimated error is larger than the corresponding LF model errors. Naively speaking, this implies that the BF approximation may not lead to significant error improvement over the LF models, unlike in Section 5.1 for the heat flux data. Fig. 13 (b) confirms this observation. To compare the abilities of the LF and BF models to reconstruct the ∆T /T 0 temperature profile at the probe location, the results from Fig. 13 are considered. Fig. 13 (a) displays the average ∆T /T 0 temperature profile derived from the LF, BF, and HF models. At most points along the profile, the mean BF1 and BF2 ∆T /T 0 are observed to be more accurate than the mean of either LF model. To quantify this error, Fig. 13 (b) provides the 2 error (see (25)) evaluated at each point along the profile. While the BF models are more accurate than the LF models at most points, the interior of the profile exhibits the greatest improvement. In the following, the focus is placed on four different QoIs of ∆T /T 0 : the spatial mean of ∆T along the profile, and ∆T /T 0 quantities at points y/W = 0.5, y/W = 0.1, and y/W = 0.05 along the profile. Figs 14 (a)-(d) provide the simulated values of the four ∆T /T 0 QoIs when using the HF, LF, and BF models. Note these do not include the data used for the BF approximation basis. Fig. 14 (a) and (b), which provide the mean value of ∆T /T 0 and the value of ∆T /T 0 at y/W = 0.5, respectively, show improved performance of the BF approximation compared to the LF models for all simulations. With regards to the accurate. On the other hand, in Fig. 14 (c) and (d), which provide the values of ∆T /T 0 at y/W = 0.1 and y/W = 0.05, respectively, error improvement does not appear to be significant. For both QoIs, the BF approximations have a smaller 2 error than their corresponding LF models; however, the gain is smaller than 2×. These results suggest that the approximation performs well for the mean and ∆T /T 0 QoIs near the interior of the profile, but accuracy decays for values closer to the walls. While prior temperature results were calculated directly from the simulated values, to estimate the moments and PDFs, as well as to perform sensitivity analysis, surrogates of the LF and BF models are formed via sparse PCE approximations. For each of the four temperature QoIs, the mean and CoV estimates determined from the sparse PCE coefficients are provided in Table 4 (a)-(d), as well as the relative error between each QoI mean and the HF QoI mean. For each subtable in Table 4, the BF mean estimates are more accurate than either of the LF estimates, by a factor of 3 − 11. Between the two BF approximations, BF2 consistently has a small error. Aside from Table 4 (b), in which BF1 is more accurate (2% error vs. 4% error), the BF2 errors are the smallest for each subtable. In terms of the QoI CoV estimates, the values of all models are comparable to that of the HF data. In Table 4 (a) and (b) the CoV results are similar to those of the heat flux results; specifically, since the BF CoV is the same order of magnitude of the respective errors, it is not necessarily distinguishable and thus can't completely be relied on as an estimate for the HF CoV. However, as the BF approximations exhibit small errors and CoV values, the mean estimates are reliable. On the other hand, the BF CoV estimates of Table 4 (c) and (d) are larger than the corresponding error, and therefore are more dependable than the estimates in  To estimate the PDFs of the four ∆T /T 0 QoIs, histograms of the LF and BF sparse PCE surrogates are formed using 25, 000 samples. These results are provided in Fig. 15, where the left column results are associated with LF1 and BF1 surrogates, and the right column results are associated with LF2 and BF2 surrogates. For comparison available simulated HF data are provided as well. With the exception of Fig. 15 (e) and (g), where the LF1 and BF1 results are closely overlaid, all of the histograms show the BF results more accurately approximating the HF data than the LF results. This significance is observed more so for Fig. 15 (a)-(d), which displays the QoI histograms of the spatial mean of ∆T /T 0 and ∆T /T 0 at the center of the profile. In Fig. 15 (c)-(h), which shows the ∆T /T 0 QoIs at specific points along the profile, the BF surrogates are consistently accurate. The LF surrogates are the least accurate in the center of the profile (Fig. 15 (c) and (d)), and improve for QoI point estimates near the wall (Fig. 15   The last result considered is global sensitivity analysis via Sobol' indices, as calculated from the sparse PCE coefficients. From these estimates, two sets of results are presented: comparisons between the four models and comparisons between the four QoIs. To compare the four models, sensitivity analysis is completed for the spatial mean ∆T /T 0 QoI. Fig. 16 provides the decomposition of important parameters from the (a) LF1, (b) BF1, (c) LF2, and (d) BF2 models. All pie charts suggest that heat flux from the radiated wall to the fluid (ξ 12 ) is the most important parameter affecting the QoI variance. However, it is important to note the results of Table 4 (a); specifically, the LF1 and LF2 errors are larger than the corresponding CoV estimates, indicating that the associated sensitivity analysis is not necessarily reliable. The BF1 and BF2 data, on the other hand, have an error that is smaller than the CoV, but on the same order of magnitude. As such, Fig. 16 (b) and (c) show that the heat flux from the radiated wall and opposite wall to the fluid (ξ 12 and ξ 13 , respectively) are the two most important input parameters. Further conclusions cannot be made with regard to the remaining parameters as the CoV estimates of Table 4 (a) are of the same order as the error.   Table 2.
Based on the moment estimations of Table 4 and histograms of Fig. 15, it is clear that the BF approximations provide improved information of the HF QoIs compared to the corresponding LF models. As a consequence, Fig. 17 provides the importance of input parameters of the four QoIs as determined by the BF2 sparse PCE coefficients. The BF2 model is selected as it consistently has a low error (see Table 4), and does so with lower cost than the BF1 surrogate model. Fig. 17 (a)-(d) shows how the importance of parameters changes with the four QoIs. From Fig. 17 (a)-(d), the heat flux from the radiated wall to the fluid (ξ 12 ) is the most important parameter for the mean ∆T /T 0 QoI, as well as the two point QoIs near the radiated wall (Figs. 17 (a), (b), and (c), respectively); in contrast, heat flux from the opposite wall to the fluid (ξ 13 ) is the most important parameter for the ∆T /T 0 QoI in the middle of the profile (see Fig. 17 (d)). This suggests that, over the whole profile, variations in the heat flux from the radiated wall will greatly affect the ∆T /T 0 values, but more so at points close to the radiated wall. For points further from this wall, the variations in this heat flux will play less of a role in the variations of ∆T /T 0 .
In terms of the remaining parameters, the results of Table 4 (a) and (b) indicate that there is lack of sufficient precision in the variance estimate to conclude the importance of the remaining parameters for Fig.  17 (a) and (b). However, the CoV estimates in 4 (c) and (d) are significantly larger than the corresponding error estimates for the BF2 models, allowing for further conclusions to be made from the sensitivity analysis results. Specifically, the data of Fig. 17 (c) and (d) show that the remaining parameters contribute equally to the variance of the ∆T /T 0 QoIs near the radiated wall.
.pdf (d) BF2, y/W = 0.05 Figure 17: Importance of input parameters from sparse PCE coefficients from the BF2 model for the ∆T /T 0 QoIs of (a) spatial mean along profile at probe location (b) y/W = 0.5 along profile at probe location, (c) y/W = 0.1 along profile at probe location, and (d) y/W = 0.05 along profile at probe location. Starting with parameter ξ 1 at the top position, importance of each ξ i is provided in counterclockwise order with respect to increasing i, with corresponding description provided in Table 2.
The ∆T /T 0 results of this section show that the BF approximations provide a more accurate reduced model than either of the LF models, with greatest improvements for mean ∆T /T 0 and ∆T /T 0 at the interior of the profile. In addition, theoretical error results guarantee that the BF approximations will be at least as accurate as the corresponding LF data. As will be shown in the following, the cost of this approximation is close to that of the LF models when many simulations are required.

Computational Cost Comparisons of the Five Models
The final component necessary to justify this BF approximation is to perform a cost comparison of all five models. Fig. 18 provides the approximate number of core-hours needed to generate N converged simulations from the HF, LF, and rank r = 6 BF models, extrapolated to large values of N . The number of simulations generated in this study are indicated by markers. For the BF1, the cost to generate N = 128 samples is 20× less expensive than the HF model, and for the BF2, the cost to generate N = 256 samples is 50× less expensive than the HF model. In comparison, the LF1 and LF2 models are 170× and 1300× less expensive, respectively; however, as shown in the results, they are poor approximations to the HF data. As the r HF simulations greatly affect the cost of the BF approximation, significant cost reduction is observed for larger values of N . For O(10 3 ) samples, which corresponds to values of N that are of interest in the context of the application studied in this work, the computational cost of the BF models more closely aligns with the cost of the LF models. While N = 10 3 HF simulations is approximately 500M core-hours, obtaining the equivalent number of simulations is approximately 6M and 3.5M core-hours for the BF1 and BF2 approximations, respectively. This drastic cost improvement, without a significant loss of accuracy that is observed with the LF models, makes the BF approximation a powerful tool to perform UQ for large-scale problems. Figure 18: Number of core-hours on Mira (ALCF) [49] to obtain N simulated values from each model. Markers provide number of core-hours for N = 128 and N = 256, to provide the cost of the simulated values generated for this study.

Conclusions
As the results of this work show, when solutions to large-scale parametric/stochastic problems exhibit a low-rank structure, i.e., lend themselves to a reduced basis representation, models with lower fidelities (cheaper to simulate) may be utilized to generate accurate approximations of the solution with significantly lower computational cost. An instance of such a BF approximation was utilized to quantify the uncertainty in thermal solutions of interest (time-averaged heat flux and temperature near the outflow boundary) of a particle-based solar receiver model involving turbulence, particle transport, and radiative heat transfer, an example of large-scale, multi-physics systems. The sources of uncertainty included particle properties (restitution coefficient, Nusselt number, absorption and scattering efficiencies), thermal boundary conditions, mass loading ratio, among others. A single HF simulation of the problem, consisting of DNS flow model, Lagrangian particle transport, and highly resolved DOM discretization of the radiative transfer equations, requires about half a million core-hours. This, in conjunction with the relatively large number of uncertain inputs, renders standard UQ techniques infeasible.
To tackle the computational cost issue, two lower-fidelity models of the problem were constructed based on coarsening the fluid and DOM grids, as well as by reducing the number of particles. These models were used to identify a reduced basis and an interpolation rule for the thermal QoIs of the HF model. As a result, only a small number of HF model realizations (e.g., six) at selected samples of uncertain parameters were needed to generate BF samples of the QoIs, which were in turn used for a PCE-based estimation of the QoI statistics and sensitivity analysis. All simulations were performed using the same computational code and in a black box fashion. The errors in predicting the QoIs and their statistics via these BF models were computed and compared to those of their LF counterpart. It was observed that the BF solutions, while requiring a small number of HF realizations, were considerably (as high as 100× for the mean of heat flux) more accurate than the LF estimates. The adopted BF strategy features an error bound which uses a small number of HF realizations (along with LF samples) to estimate the error with respect to the HF model. The efficacy of the bound in determining the number of required HF samples and estimating the BF error was also demonstrated on the thermal QoIs.
The construction of the BF models enabled us to efficiently carry out sensitivity analyses via PCE-based Sobol' indices. The results indicate that the thermal QoIs considered in this work are extremely sensitive to the heat fluxes from the walls to the gas-particle mixture resulting from the interaction between radiation and non-ideal transparent walls. The importance of this effect is amplified in wall-bounded, particle-laden turbulent flows as turbophoretic mechanisms tend to accumulate particles at the walls where the importance of non-ideal transmissivity is the highest. The analysis also extracted that uncertainty in particle properties, which is typically disregarded in computational studies of irradiated particle-laden turbulence, plays an important role in the solution outcome that may result in large variability in, otherwise, robust first-and second-order statistics.
Ongoing and future work focuses on exploring the performance of the BF approximation under more challenging physical conditions in terms of turbulent intensities, particle loading, and radiation intensity. Of notable interest is also the assessment of the approach in predicting higher order statistics, e.g., particles and temperature fluctuations, and more complex physical phenomena, such as particle clustering and turbulence modulation.