Uncertainty Quantification of Time-Dependent Quantities in a System With Adjustable Level of Smoothness

Fluctuating operating conditions are ubiquitous in applied environments and need to be considered when performing computational estimates of the relevant system characteristics [1]. This can be achieved by utilizing a technique termed uncertainty quantification (UQ). This is a methodology aiming to understand how predictions of a numerical model are affected by sources of uncertainty. One of the types of UQ problems is the forward propagation of uncertainty, which estimates the resultant uncertainty in model outputs propagated from uncertain inputs, in the form of low-order moments of the outputs, i.e., mean and variance. Within UQ, a given model output is referred to as the quantity of interest (QoI).

In this work, we consider computational fluid dynamics (CFD) models subject to uncertain inputs in which the input–output relation is governed by the Navier–Stokes equations, which are nonlinear partial differential equations (PDEs). For a review of the available UQ methods with a particular focus on CFD applications, the reader is referred to Refs. [2] and [3].

In UQ, the moments of the output distribution of QoI(s) may be computed via random sampling methods and nonrandom methods that are strategic in their sample selection. Random sampling methods include Monte Carlo (MC) [4], Quasi-Monte Carlo [5], and importance sampling [6]. Nonrandom methods include response surface [7,8] and polynomial chaos expansion (PCE)-based methods [9]. The choice of the UQ method depends on the system dynamics through which the uncertainty is to be propagated. The system dynamics under consideration refers to the smoothness of the underlying mathematical model responsible for calculating a QoI. For the purpose of our work, we understand smooth in the context of practical uncertainty quantification. A time-dependent QoI is considered to be smooth if its frequency is low enough to be surrogated with low-order nonsampling methods such as PCE [10].

If QoI exhibits smooth behavior, as is often demonstrated by a range of mathematical models, polynomial chaos expansion techniques are best suited to address such problems. This is because PCE is capable of accurately approximating smooth models using spectral representation of the uncertain quantities by expansion in an orthogonal polynomial basis. Evaluating a polynomial representation requires much less computational resources compared to the evaluation of an actual model, which makes the sampling of the parameter space far less computationally expensive [11]. When PCE is used to quantify uncertainty in a smooth system with finite variance, it provides exponential convergence at relatively low computational costs [12], and therefore PCE is the most commonly used UQ propagation technique when dealing with smooth processes [13].

The term convergence is used if an uncertainty propagation technique manages to determine the relevant statistical moments of a QoI that have converged and would no longer benefit from further increase in the sample size, providing an insight into the realistic system response when subjected to uncertainty. A synthetic example demonstrating converged and nonconverged UQ scenarios is given in Fig. 1.

Despite the clear advantages of the PCE method when used to quantify uncertainty in smooth processes, it fails to build an accurate surrogate model when applied to a system where time-dependent QoI(s) exhibit nonsmoothness and as a result, PCE calculates moments of a QoI that are not converged [14,15]. This occurs because as the nonsmooth QoI evolves, the degree of PCE's expansion polynomials must rapidly increase to capture the evolution of a QoI and maintain the accuracy of the expansion [16–20]. For this to occur, the system does not even have to be chaotic—it can be nonchaotic but of nonlinear nature [16]. As these types of systems evolve, their solutions start to develop their own random characteristics. Thus, the probability density function of a QoI evolves with time, and the PCE's orthogonal polynomials, which are optimal initially, lose their optimality at later instants. The increase in polynomial order only postpones the problem to a later instant in time. For a specified polynomial order, the error adds up and becomes unacceptable after some time. In short, PCE cannot build accurate surrogate models, including CFD-based models [21] that are nonsmooth. Clearly, if the surrogate fit is not representative of the underlying model, the PCE-based forward propagation of uncertainty will result in non-convered moments of a QoI.

This is restrictive because in the area of fluid mechanics, for example, flows and corresponding QoIs, such as velocity magnitude or pressure, develop nonsmoothness beyond a certain Reynolds number rendering PCE as ineffective. However, the ability to quantify uncertainty in nonsmooth simulations (such as turbulent CFD models) at manageable computational costs is recognized as an important area for development [22]. This has resulted in authors attempting to make PCE resistant to chaos, nonlinearity and the resulting nonsmoothness. To this end, Kantarakias et al. [13] coupled PCE with the multiple shooting shadowing [23] method and applied it to two standard chaotic systems, the Lorenz system and the Kuramoto–Sivashinsky equation. It was found that at the quadrature integration points of PCE, the resulting regularization lead to accurate results that matched well with the Monte Carlo simulations at significantly lower computational costs. Paulson et al. [24] developed a nonsmooth polynomial chaos expansion (nsPCE) method as an extension of traditional PCE that can handle the nonsmooth nature of computationally expensive dynamic models. The main idea behind nsPCE is to systematically partition the parameter space into two nonoverlapping elements on which the model response behaves smoothly. Separate PCE models are then fitted in both elements using a basis-adaptive sparse regression approach. The effectiveness of the method was demonstrated on a nonsmooth dynamic flux balance analysis model of an E. coli monoculture.

PCE's inability to build accurate surrogate models has been documented for unsteady flows, where each realization is slightly different in phase [16,19,21,25] and also in turbulent flows [26–29], which puts a severe limitation on applying the propagation technique in nearly any realistic fluid flow simulation. The only way to bypass this and apply PCE to a turbulent or unsteady fluid flow simulation is to build a surrogate based on specific QoIs that are smooth and steady despite the flow itself being turbulent or unsteady. Examples of smooth QoIs in non-smooth CFD simulations include statistical moments of decaying homogeneous isotropic turbulence [30] and root-mean-square values of velocity, temperature, and mixture fraction at specific domain locations [31]. Usually, an applied user would like to assess more typical QoIs such as temperature or pressure, which usually represent the general fluid-flow characteristics that are often turbulent. In summary, PCE cannot be applied to nonsmooth QoIs; thus, the following choice is presented when quantifying uncertainty in nonsmooth systems: use sampling methods such as Monte Carlo sampling or use PCE coupled with a mechanism that can somehow deal with nonsmooth system dynamics. However, the former method is computationally restrictive as it relies on the law of large numbers needing many model evaluations for statistical convergence and the latter is still in the area of novel research that is most likely to be unavailable to an applied user.

Thus, the aim of this work is to take a step back and try to quantify the relationship between the PCE's convergence behavior and the dynamics of the underlying system. Specifically, we want to understand whether there is a quantifiable relationship between the PCE's convergence metrics and the level of smoothness of quantatites of interest.

To this end, we perform an uncertainty quantification study of CFD-based combustion simulations with contrasting system dynamics: a real turbulent burner flame and a theoretical laminar flame front. This aims to showcase the effect of flow state and the underlying system dynamics on the convergence performance of the PCE technique. In both cases, the forward propagation of uncertainty in the injection speed of the fuel stream will be performed with respect to temperature, T, in the area of flame stabilization; thus, T is the QoI of this study. Setting the flame temperature as the QoI is a common direction of UQ analyses of combustion models [32].

The above-mentioned approach allowed to empirically show the existence of a power law scaling linking the convergence performance of PCE in relation to MC to the degree of smoothness of a QoI as measured by its standard deviation. We believe that the established metric will allow those requiring UQ to make a more informed decision about the uncertainty propagation technique to be used depending on the degree of smoothness exhibited by the problem.

The reason for using combustion models, specifically the burner model as the basis of this study, is due to the fact that such models present the opportunity to analyze the temperature response at different simulation time steps, each representing a unique model solution, thus representing a unique QoI response. In such a model, different simulation time steps correspond to different levels of temperature smoothness enabling one to analyze the effect of QoI dynamics on the PCE's convergence performance within the same system by performing UQ analysis at different simulation times. Furthermore, because the burner flame is a continually and naturally developing structure, there are no issues arising from the fact of capturing identical model solutions at different simulation time steps such as in a CFD model of a flow past a cylinder.

The focus of the study is to assess the convergence performance of PCE at varying levels of system smoothness. To analyze the change in the PCE's performance, its convergence metrics are compared against corresponding results obtained by means of MC sampling.

For the MC method, a number of independent samples, M is drawn randomly from probability-approximated uncertainty, ω and these are then propagated through the underlying model, f. Having obtained M model evaluations at the time of interest, t, the moments (i.e., mean or variance) of $f(t,ω)$ are computed, making it a very straightforward yet robust, assumption-free and easy to parallelize UQ technique. The higher M, the better the moments of a QoI will have converged. Mathematically, ω is a probability distribution of a given type that best describes the relevant uncertainty. An example of a probability-approximated ω would be a normal distribution with mean, μ and standard deviation, σ: $ω∼N(μ,σ2)$⁠. Other distributions can be employed if necessary.

Unlike the MC procedure, the PCE methodology is a nonsampling-based technique capable of approximating the model response at t and ω, $f(t,ω)$ as a spectral expansion in terms of orthogonal polynomials, $ϕn(ω)$ and expansion coefficients, $f̂n(t)$—a procedure in some ways analogous to Fourier decomposition. However, rather than employing sines and cosines as an orthogonal base for the decomposition, other functions are used. For example, if the underlying distribution of the uncertain parameters, ω is Gaussian, then the associated orthogonal polynomials are Hermite polynomials and if ω is uniform, they are Legendre polynomials [33]. In general, one can construct orthogonal polynomials for arbitrary distributions [34]. For practical computations, the spectral expansion, mathematically formulated as an infinite series is truncated to Nth order and rearranged for $f̂n(t)$ exploiting the orthogonality of the underlying basis. The rearranged expression takes the form of an integral that is numerically approximated using Gaussian quadrature [35], which generates a number of nodes, $xp$ needing model propagation and weights, $wp$ used in the integral approximation according to the quadrature degree, P. In theory, the higher the P, the higher the integral approximation accuracy will be and subsequently more converged moments are attained via the PCE method. Having obtained the expansion coefficients at the time of interest, $f̂n(t)$⁠, the moments of $f(t,ω)$ are computed solely based on the obtained coefficients. The steps involved in applying the algorithms of the two methods are summarized in Table 1, which demonstrates the difference in the implementation complexity of the two techniques.

The two nonintrusive UQ techniques considered, PCE and MC, are applied to two benchmark combustion cases. These are Sandia flame D and counter flow flame, which were modeled with CFD. Their nature and features are briefly summarized below.

Sandia flame D is the well-known air-piloted coaxial methane jet burner experiment conducted by the Sandia National Laboratories [36]. The turbulent flame operates at a high Reynolds number of 22,400 with little or no local extinction even at modest pilot. The mixing rates are sufficiently high such that the flame speed is limited by the diffusion rate with a well-defined reaction zone. Despite those good operating characteristics, the burner flame is not a demanding simulation with a simple geometry and is very well documented providing high quality experimental data. This made the case one of the most widely used benchmark problems in the validation of new numerical combustion models [37–40]. Briefly, the fuel stream is injected at 49.6 m/s from a round nozzle with a diameter, d of 7.2 mm and comprises 15% methane, and 85% air by mass. The burner is situated in a coflow with a bulk velocity of 0.9 m/s.

Counter flow flame is a combustion system in which pure methane gas and air move toward each other at equal low speeds of 0.1 m/s and upon sufficient mixing react around the center of the symmetrical domain, forming a laminar flame front. From the start, the computational domain is heated to $2000 K$ providing the necessary activation energy. openfoam's counter flow flame case is based on the work around validated axysymmetric counterflow flame simulations by Johnson et al. [41].

The numerical simulation of both cases was carried out using the open-source CFD tool—openfoam v7. In fact, both Sandia flame D and counter flow flame cases were simulated as part of the openfoam in-built combustion tutorial family and did not require any modifications except for a mesh refinement in the Sandia flame D case. For reference, the openfoam-specific names of these cases are “SandiaD_LTS” and “counterFlowFlame2DLTS_GRI_TDAC”. Both cases were simulated using the transient, compressible solver for nonpremixed combustion of gaseous fuels named reactingFoam. The reactingFoam solver works through the Favre-filtered Navier–Stokes and continuity equations to solve for momentum, species mass fraction and energy to simulate compressible, reacting flow. For a detailed description of the reactingFoam solver refer to Andersen et al. [42]. Briefly, the solver is based on the PIMPLE (pressure implicit split operator-semi-implicit method for pressure-linked equations) algorithm for pressure–velocity coupling and is stabilized by the local time stepping (LTS) scheme [43], which maximizes the individual time-step for each individual cell size. In practice, a single LTS iteration corresponds to a unique $Δt$ value; thus, the physical time of the simulation is equal to the sum of all individually calculated $Δt$ values for each LTS iteration.

When simulating Sandia flame D case, it was first ran for 1500 LTS iterations with chemistry deactivated to achieve a cold profile and then with chemistry activated until 8000 LTS iterations—the time of flame stabilization. The flame stabilization in the counter flow flame case is achieved after 1500 LTS iterations. In this study, the number of LTS iterations will be used as a measure of the simulation progress in time. For further reference, the simulation state at the LTS-time correspondent to flame stabilization is referred to as stabilized flow.

The submodels used in conjunction with the LTS-stabilized reactingFoam solver for both cases are shown in Table 2. The employed submodels used in the openfoam Sandia flame D simulation are those submodels that have been found to produce the best results in terms of their similarity to the experimental data [49].

The computational domain housing the resulting temperature field at the specified boundary conditions of the unmodified counter flow flame case at the time of flame stabilization (1500 LTS) is shown in Fig. 2. For this case, the forward propagation of uncertainty is performed with respect to the rectangular sampling zones 1 and 2, which cover most of the perpendicular area of the formed flame stabilization front.

The diameter (d)-normalized computational domain of Sandia flame D housing the temperature field, simulated using the described implementation approach, at the time of flame stabilization (8000 LTS) and at the specified boundary conditions is shown in Fig. 3(a). The computational domain and mesh were obtained from the original openfoam Sandia flame D tutorial. The domain is reduced to an axisymmetric 20 deg wedge, exploiting the flame's rotational symmetry to save on computational cost. The only modification made to the original mesh was that it was refined from 5170 cells to 25,850 cells. This is in line with work by Kadar [50] who found that for the reactingFoam-based simulation of Sandia flame D the results with 23,340 cells are optimal both in terms of estimation accuracy and computational cost. As shown in Fig. 3(b), the results obtained from the simulated Sandia flame D case are in good agreement with the experimental temperature, which was measured along the diameter(d)-normalized axis, x/d as well as in the normalized radial direction, r/d.

For Sandia flame D, the forward propagation of uncertainty is performed with respect to the 103 sampling zones. These are equally sized rectangles covering the entire flame stabilization area, as illustrated in Fig. 3(a). This number of 103 sample areas arises when trying to balance two factors: avoiding the excessive averaging of the QoI and computer storage limitations. In principle, it would be desirable to sample QoI results from every computational point of the domain (23,340 points) in order to avoid any QoI-averaging. However, this is prohibitive in terms of storage requirements. This is because for each sampling zone, the QoI results are recorded for every computational time-step (8000 time steps in total). This is needed as part of the upcoming analysis aimed at finding the time when the QoI has the highest standard deviation and with many model evaluations needing to be done as part of the MC and PCE analysis, it makes it prohibitive to sample QoI from each computational point.

The use of sampling zones, in the form of rectangles, lines, and circles, is a normal practice as it is more often that an applied user would like to estimate the resultant uncertainty in model outputs in a sizeable part of the computational domain rather than at a singular computational point.

Having obtained a validated CFD model of the Sandia flame D case and a CFD model of the unmodified counter flow flame case, the uncertainty quantification and subsequent convergence study may now be performed. First, an uncertainty needs to be introduced into the simulations. For both simulations, it was decided to propagate the uncertainty, approximated as normal probability distribution, in the injection velocity of the fuel stream: for Sandia flame D the fuel stream is termed as “inletCH4” and for counter flow flame “fuel inlet” as shown in Figs. 3(a) and 2. The uncertainty in the injection speed of the fuel stream in the Sandia flame D experiment [36] was reported to be $49.6±2 m/s$ and is propagated in this study. For the counter flow flame, a similar uncertainty of 5% of the mean value is propagated. In both cases, forward propagation of uncertainty in the injection speed of the fuel stream is performed at the respective times of the flame stabilization for the QoI of T. Mathematically, these uncertainty frameworks can be expressed as:

The convergence study will be based on the comparison of the MC and PCE obtained temperature means, T-means at the selected sampling zones. Both techniques are allocated the same computational allowance for effective technique comparison, which is measured in the number of model evaluations the technique requires to obtain a given T-mean value in a given sampling zone. The number of model evaluations ranges from 2 to 16. Since the UQ problems in question are 1D, the increasing quadrature order of the PCE method does not result in the exponential growth of the number of nodes needing model propagation; thus, for a given number of model evaluations in 1D setting, the number of MC runs, M is equal to the PCE's quadrature order, P (see Table 1).

The outlined convergence study, with respect to the defined Sandia flame D uncertainty framework, is performed for three sampling zones 27, 35, and 60 as shown in Fig. 4(a). Note, for these three sampling zones, the upcoming analyses will be done explicitly for demonstrative proposes and there is no specific reason for the choice of these zones apart from the fact they are from different parts of the computational domain. However, the overall analysis will be conducted for all 103 sampling zones. The convergence study, with respect to the defined counter flow flame uncertainty framework, is performed for sampling zones, 1 and 2, as shown in Fig. 4(b).

In this study, the convergence of a propagation technique is quantified as a variance of all the obtained T-mean values (15 in total). An uncertainty propagation technique with lower variance of the obtained T-means can be deemed as the better, more performant technique with higher convergence. To demonstrate the use of the introduced terms, consider the convergence plot for Sandia flame D's sampling zone 27 shown in Fig. 4(a): the variance of the PCE-obtained means was 5.8 and 1.6 for MC-obtained means, making PCE the less performant propagation technique in this particular case.

As evident from Fig. 4(a), the T-means calculated using both propagation techniques for all three selected sampling zones (27, 35, and 60) of Sandia flame D fluctuate significantly, even with PCE at higher orders of quadrature. Both MC and PCE demonstrate poor convergence properties. Moreover, the PCE-obtained T-means show no better convergence than the MC-obtained T-means at the same number of model evaluations.

With the counter flow flame simulation, the convergence behavior is opposite to that in Sandia flame D, with PCE being significantly more performant than MC, converging to a T-mean value in an instantaneous manner, requiring only two model evaluations (second order of quadrature) with PCE, for sampling zones 1 and 2 as can be seen in Fig. 4(b). As we can see for sampling zone 2, the variance of the PCE-obtained T-means is 2.1 and 427.8 with MC, which suggests that PCE calculates significantly more convergent statistics than MC given the same computational resources. The PCE's superior to MC convergence behavior is also observed for sampling zone 1.

This result prompts the question regarding the origin of the difference for the two cases in respect of their PCE convergence behavior. To establish this, a sensitivity study is carried out in which both models are subject to 16 fuel injection velocities correspondent to the nodes generated at the 16th order of Gauss quadrature according to the respective ω definitions and propagated as part of the PCE algorithm (see Fig. 5).

Analyzing Fig. 5, it is evident that in the Sandia flame D's case at 8000 LTS, a high degree of QoI nonsmoothness is observed. Although, the curves demonstrate predictable, related one-to-another behavior at first, past the point of reaching approximately 2000 K they diverge and eventually appear to be unrelated, fluctuating around the flame temperature, which can be described as chaotic behavior. As established, systems exhibiting chaos and, as a result, high level of QoI nonsmoothness cannot be approximated by the PCE's spectral expansion; thus, the PCE-UQ analysis shows poor results that are no better than that of small-sample MC.

In the counter flow flame case, the system does not exhibit nonsmooth QoI behavior at any point during the simulation. At the time of flame stabilization (1500 LTS), each quadrature node corresponds to a clearly defined, predictable temperature trajectory, which in turn results in each quadrature node contributing its clearly weighted role in the construction of an accurate PCE surrogate model and subsequent highly convergent PCE-UQ analysis as was seen in Fig. 4(b).

Within the CFD framework, the above differences in the system dynamics can be illustrated by plotting the standard deviation map of the QoI. To this end, the local standard deviation in temperature, σ_T is calculated based on 100 MC samples drawn from the respective ω definitions. The result of this procedure is illustrated in Fig. 6. The figure reveals that the counter flow flame case features a well-defined and relatively high σ_T front at and near the position of flame stabilization. On the contrary, the Sandia flame D case features a dispersed and relatively low σ_T front at and near the position of flame stabilization.

Based on the obtained maps of standard deviation, a conclusion can be drawn that the PCE UQ analysis of sampling zones housing higher standard deviation of the QoI significantly outperforms MC-based forward propagation of uncertainty. This is clearly evident from the analysis of the counter flow flame case, where PCE exhibits much better convergence in relation to MC than in the case of Sandia flame D (see Fig. 4) and features higher standard deviation of the QoI of temperature as can be seen in Fig. 6. Therefore, the analysis around the PCE's convergence and the standard deviation of QoI will be explored further. Note, this further analysis will only be performed for the case of Sandia flame D as its simulation suffers from the poor PCE convergence properties. There is no need to do any further analysis on the counter flow flame case as PCE demonstrated excellent convergence, which was expected as it is a laminar, smooth system. However, the counter flow flame case played an important role in setting the direction of the upcoming analysis by indicating that the PCE's convergence may improve with an increase in standard deviation of a QoI—in our case, temperature.

To conduct the analysis on the effect of the standard deviation of the QoI on the PCE's convergence, a method that increases σ_T found in sampling zones needs to be devised. If the increase in Sandia flame D's σ_T is achieved, it will allow to test the following hypotheses: the higher the standard deviation of the QoI in a sampling zone, the more performant PCE becomes relative to MC, in terms of the QoI's statistical convergence properties, when there is an uncertainty affecting that QoI.

A means to achieve the standard deviation increase in the Sandia flame D's sampling zones, without altering the model itself, nor the defined UQ framework, is to perform analyses at a different simulation time. A way to find the time, t at which a sampling zone experiences the highest standard deviation of the QoI, in our case $maxσT$⁠, is to find the time at which the temperature (QoI) in that sampling zone has the highest standard deviation. To this end, the local standard deviation in temperature, σ_T is calculated for all simulation time steps (and not just at the time of the flame stabilization) based on the same 100 MC ω-based samples used in the construction of the standard deviation maps shown earlier. The result of this procedure is shown in Fig. 7. Note that for each individual sampling zone, the time at which it experiences $maxσT$ is different as the flame development is dependent on the simulation time and will interact with each sampling zone on an individual basis. Therefore, the zone that has not been reached by the flame front at a given simulation time cannot physically experience $maxσT$⁠. For further reference, the flow state at the LTS-time correspondent to $maxσT$ is referred to as developing flow.

Applying the aforementioned strategy for finding the time (LTS time-step) at which a sampling zone experiences $maxσT$ to the previously analyzed sampling zones 27, 35, and 60 of Sandia flame D resulted in the data displayed in Fig. 8. For each analyzed zone, the data includes:

• Convergence plot at the time of flame stabilization, 8000 LTS, which was already shown in Fig. 4(a) (it is included for comparison)

• Convergence plot at a newly calculated LTS-time corresponding to $maxσT$⁠, which is termed as developing flow

• Standard deviation map of temperature at stabilized and developing flow states

• Temperature profile at stabilized and developing flow states

Referring to Fig. 8: for sampling zone 60 at the time of flame stabilization (8000 LTS), the variance of the PCE-obtained T-means was 53.5 and the variance of the MC-obtained T-means was 34.0 rendering PCE as ineffective in this case. However, at the time corresponding to the developing flow state (5432 LTS), the variance of the PCE-obtained T-means was 1236.4 and the variance of the MC-obtained T-means was 22037.4, making PCE the more performant propagation technique calculating significantly more convergent QoI statistics than MC. Therefore, at developing flow state, the situation for sampling zone 60 changes with PCE having significantly better convergence than MC, turning the PCE into a far more performant uncertainty propagation technique. Although, the variances of the T-means obtained via the two techniques are much higher at developing flow state than at fully stabilized flow; it is the relative difference in convergence between the techniques that is of interest. The same trend of significant improvement in the PCE's convergence in relation to MC when utilized during the developing flow states, applies to sampling zones 27 and 35.

Another trend that can be observed is that a sampling zone experiences $maxσT$ when the outer shear layer area of the developing flame interacts with that sampling zone. This consistent interaction indicates that the MC-based strategy of finding the flow state at which the sampling zone experiences $maxσT$ is devised correctly. Physically, what is happening is that a sampling zone (at time corresponding to the developing flow state) interacts with the outer shear layer area of the developing flame, which is the least turbulent part of a flame [51]. In other words, we are creating a time window during which the forward propagation of uncertainty is performed for a sampling zone when it interacts with the least turbulent (i.e., least nonsmooth/least chaotic) part of the flame. This results in an improved PCE convergence due to the smoother temperature-QoI response, which is quantified by the increase in its standard deviation.

To demonstrate the use of Eq. (1), let us calculate $ΔMC−PCE$ parameter for sampling zone 61 of Sandia flame D at stabilized and developing flow states (see Fig. 8):

The absolute difference is taken since the $ΔMC−PCE$ values for the stabilized flow data can be negative because at lower standard deviation of QoI, PCE often demonstrated slightly worse convergence than MC for the same number of model evaluations, as was demonstrated in Fig. 8. Since σ_T and $ΔMC−PCE$ values for the developing flow data are significantly larger than those for the stabilized flow, logarithmic representation is needed to show the $ΔMC−PCE$ versus σ_T data points for both flow scenarios together, which requires only positive data. Therefore, if $ΔMC−PCE$ is $≫0$⁠, MC-obtained T-means have higher variance than the PCE-obtained T-means, suggesting that the PCE is the more performant technique with better convergence and vice versa. Thus, the higher the $ΔMC−PCE$ parameter the more performant PCE is in relation to MC.

Having developed a metric quantifying the PCE's convergence performance in relation to MC, it can be used to observe the changes in the PCE's convergence performance with increasing standard deviation of the QoI found in a sampling zone. To this end, the $ΔMC−PCE$ parameter is plotted against every sampling zone's standard deviation of the QoI, σ_T, in double-logarithmic representation, at both stabilized and developing flow states (see Fig. 9). Figure 9 reveals that the relative convergence performance of PCE exhibits different behavior at the two analyzed flow states:

• Stabilized flow: σ_T is low as was seen in Fig. 6 and the PCE's convergence was similar (or even worse) to that of MC resulting in low $ΔMC−PCE$ values as was seen in Fig. 4(a). Thus, the data is clustered in the region of low σ_T and low $ΔMC−PCE$⁠. However, the data does demonstrate a weak trend where the performance of PCE in relation to MC improves with higher standard deviation of the QoI, which at the stabilized flow conditions is low and uniform but with some σ_T variability across the domain

• Developing flow: the performance of PCE also improves with increasing σ_T but with a much clearer dependency and reaching much higher values of σ_T (as was seen in Fig. 8—developing flow plots), which results in significantly better performance of PCE in relation to MC as defined by the $ΔMC−PCE$ parameter

Referring to Fig. 9, we see that the stabilized flow data shows no trend and is clustered in the region of relatively low σ_T and low $ΔMC−PCE$ values. This is observed because at the time of full flame stabilization, the temperature response is chaotic, which translates to the PCE's inability to construct an accurate surrogate of the QoI reducing the PCE's convergence properties, which is quantified by a high value of $Var(T−means PCE)$⁠, a component of the $ΔMC−PCE$ parameter (refer to Eq. (1)). On the other hand, low-sample MC will always be naturally associated with high values of $Var(T−mean MC)$⁠. Thus, in the scenarios, where QoI response is chaotic as quantified by its low standard deviation, $ΔMC−PCE$ values will be relatively low as there is no substantial difference between the convergence performances of the MC and PCE methods. However, this situation changes and PCE exhibits significantly better convergence behavior than MC if we are giving an opportunity for the defined sampling zones to interact with the QoI when it is nonchaotic and thereby smooth, which is the method used in obtaining the developing flow data. This occurs because PCE is capable of building high accuracy surrogates of systems with a nonchaotic response with a relatively few model evaluations resulting in much better convergence properties than those for low-sample MC. In these scenarios where QoI-response is smooth, the following is observed: $Var(T−means MC)≫Var(T−means PCE)$⁠, which means that the majority of the variance in the $ΔMC−PCE$ parameter is caused by the MC component. Thus, the plot of $ΔMC−PCE$⁠, where the PCE component is negligible, against σ_T should result in the square dependency as the standard deviation is the square root of the variance, which is observed in Fig. 9. What this tells us is that the best convergence performance PCE can possibly exhibit is proportional to the square of the QoI's standard deviation. Thus, technically, any deviation from the square dependency suggests that the QoI response is not completely smooth as there is some difficulty for it to be effectively surrogated with the PCE method. Clearly, in the case of Sandia flame D that deviation from the quare dependency is very small because each sampling zone was able to interact with the flame at a time when temperature (QoI) demonstrated smoother behavior, allowing for it to be effectively PCE-surrogated and compute QoI's moments of much higher convergence that those computed with MC. Thus, at developing flow, the plot of the MC-dominated $ΔMC−PCE$ values vs. standard deviation would naturally be a squared dependency.

The generality of the observed findings and specifically the value of the exponent in the derived power law scaling, for other fluid flows, still remains to be confirmed mathematically. Since this is not attempted, we cannot make any statements on the generality of the proposed power law scaling relationship. However, we present further empirical evidence to demonstrate that the power law scaling is not exclusive to the case of Sandia flame D. The steps involved in determining the power law scaling are repeated for two additional fluid flow cases referred to as DLR Flame A [52,53], a well-characterized turbulent jet flame and the turbulent pitzDaily flow scenario [54], which is a compressible, transient simulation with a preheated inlet stream based on a large eddy simulation (LES) turbulence model [55]. Their nature and features are briefly summarized below.

DLR Flame A is a well-documented burner test case with Re = 15,200 [52,53]. This test case is available in the openfoam tutorial family with the name of DLR_A_LTS. The case is simulated by means of the reactingFoam solver in conjunction with the Reynolds-averaged simulation (RAS) turbulence model. This burner test case is different to Sandia flame D in terms of burner geometry, fuel composition, fuel injection speed, the nozzle inner diameter (d = 8.0 mm) and the method of flame ignition. Analogous to Sandia flame D, UQ analysis of the DLR test case involved forward propagation of uncertainty in the injection speed of the fuel stream as a normal probability distribution in line with the experimentally reported data: $N(42.2,0.5)$⁠. The sampling mesh consisted of 102 zones. The domain showing the QoI of temperature, T at the time of flame stabilization (10,000 LTS), the sampling mesh and the final data structure along with the developing flow fit are shown in Fig. 10(a).

The turbulent pitzDaily case is a two-dimensional domain consisting of a short inlet, a backward facing step, and a converging nozzle as the outlet [54]. The case employs an LES turbulence model. The case is also part of the openfoam tutorial family and is a compressible fluid flow simulation solved by means of the rhoPimpleFoam solver. The QoI here was also temperature, which can be simulated as part of the rhoPimpleFoam solver. The temperature difference is induced by the hot inlet stream heated to 800 K, whereas the domain is at 292 K. The forward propagation of uncertainty was performed with respect to the injection speed of the hot stream, which is considered to be normally distributed: $N(10,0.4)$⁠. The mean of 10 m/s is the original inlet speed supplied with the tutorial and the standard deviation of 0.4 m/s represents 4% of the mean value, which is in line with the standard deviation of the uncertainty propagated through the Sandia flame D case. The sampling mesh consisted of 109 zones. The domain showing the QoI of temperature at the time when the hot temperature stream passes the backward facing step (this specific simulation time is shown for demonstrative purposes as at the time of quasi-steady-state, the domain is completely filled with hot temperature, which ineffectively communicates the nature of the fluid flow), the sampling mesh and the final data structure along with the “developing flow” fit are shown in Fig. 10(b).

It is evident that the $ΔMC−PCE∝σT∼2$ power law scaling holds for the two additional test cases presented. However, there is inevitable variability in the exponent values and the R² metrics between the cases. That variability arises from the processes associated with CFD and the nature of the fluid flows being simulated. In principle, the most important aspect which influences the exponent values and the R² metrics is the level of smoothness of a CFD-modeled QoI, which is measured by the standard deviation of that QoI. If a sampling zone features a smooth QoI that can be effectively PCE-surrogated, the fit should be close to that of a power law with the exponent of 2 and a high coefficient of determination. However, as per the premises of the study that is not a binary classification in terms of whether a QoI is capable of being surrogated or not. We state that there is a quantifiable scaling relationship between the PCE's ability to surrogate a QoI and the level of smoothness of that QoI, as quantified by its standard deviation. Thus, the exponent and R2 values of $ΔMC−PCE∝σTγ$ will always vary for different fluid flows as the level of QoI smoothness will always be unique to any particular simulation case.

Furthermore, to demonstrate that the scaling is independent of the sampling mesh, we also build the $ΔMC−PCE$ versus σ_T data structure for the case of Sandia flame D but with 53 sampling zones, instead of the original 103. The results obtained for the reduced number of sampling zones are included in Fig. 10(c). The comparison of Fig. 10(c) with Fig. 9 reveals that the power-law scaling has not been affected by a reduction in the number of sampling zones. In principle, the number of sampling zones should not affect the nature of the final data structure and the analysis could be performed for very few sampling zones as long as the sampling zones in question are able to demonstrate smooth QoI response of different levels with respect to the propagated uncertainty.

Another, important aspect that needs to be considered as it may heavily influence both the exponent value and lower the R² metric is the effects associated with the CFD modeling tools such as boundary conditions (BCs). These can sever the link between the propagated uncertainty samples and a QoI. For example, in the pitzDaily case (see Fig. 10(b)), the sampling zones were defined as such not to directly interact with a QoI at the walls of the domain because there, the BCs simply prescribe a value to a QoI as part of the CFD workings (referred to as Dirichlet boundary condition). In the DLR case, the sampling zones were defined as such not to interact with the part of the domain between $x/D=0$ to $x/D∼30$ as this is the area where a high temperature value is prescribed as a type of an initial condition to facilitate the ignition of the flame. The scenarios where QoI is prescribed, disregarding the input, weaken or even sever the link between the propagated values associated with the uncertain input and a QoI. Evidently, performing UQ analysis with respect to sampling zones that are facing interference in terms of input-QoI relation hinders the ability to calculate accurate moments by any UQ propagation technique including PCE, which, for the study's algorithms translates in significantly lower $ΔMC−PCE$ values increasing the variability of the power-law fit and affecting its exponent value. Note, these issues are not applicable to the case of Sandia flame D because there, the ignition occurs via a preheated pilot stream entering the domain via a separate inlet (see Fig. 3(a)) that does not impose its prescribed high temperatures onto any of the defined sampling zones (in both cases with 53 and 103 sampling zones) and there are no other conditions prescribing a value to a QoI anywhere in the domain.

This study has demonstrated that there exists a relationship between the convergence performance of polynomial chaos expansion (PCE) in relation to Monte Carlo (MC) sampling, as quantified by the comparative $ΔMC−PCE$ parameter, and the degree of smoothness of a quantity of interest (QoI), which is quantified by its standard deviation. The data analysis revealed a power-law scaling $ΔMC−PCE∝σTγ$⁠, with the exponent being very close to γ = 2 for the QoI of temperature and the flow scenarios investigated here.

The observed scaling relationship was derived by propagating the uncertainty in the injection speed of the fuel stream in the model of the benchmark turbulent combustion case, Sandia flame D by means of PCE and MC. At the time of flame stabilization, temperature, the study's QoI was exhibiting nonsmoothness resulting in poor PCE convergence. This was an expected result as PCE is known to have convergence problems when used to quantify uncertainty in nonsmooth systems. In fact, at the time of flame stabilization, per the defined $ΔMC−PCE$ parameter, it was deemed that the PCE's performance is no better than that of small sample MC.

Having performed an analogous uncertainty quantification (UQ) analysis in a laminar combustion model—counter flow flame with a smooth temperature response and as a result highly convergent PCE-UQ analysis—it was concluded that a quantifiable difference between the simulations that could be responsible for such contrasting convergence behavior of the PCE method was the standard deviation of QoI, temperature. It was observed that the higher the standard deviation of QoI, the more convergent model output statistics are produced by the PCE method in relation to MC. Following that observation, a method of increasing the temperature's standard deviation, without altering the Sandia flame D simulation itself, nor the defined UQ framework, was devised. This involved finding the simulation time at which a given sampling zone experienced the maximum standard deviation of temperature.This translated to an instancewhen a given sampling zone interacted with temperature in the outer shear layer area of the developing flame. This has worked and the obtained increase in the standard deviation of QoI-temperature significantly improved the PCE's convergence properties in relation to MC as quantified by the $ΔMC−PCE$ parameter. The rate of improvement in PCE's convergence as quantified by the increasing $ΔMC−PCE$ values scaled with the standard deviation of the study's QoI, temperature. These steps were repeated on further two test cases of DLR Flame A, another well-characterized turbulent jet flame and the LES-turbulent pitzDaily flow scenario demonstrating that the attained power law scaling of $ΔMC−PCE∝σT∼2$ is not exclusive to the case of Sandia flame D.

Based on the investigated flow scenarios, it is deemed important to perform a preliminary analysis aimed at identifying the underlying system dynamics, specifically the level of smoothness, which, as part of this study, was found to be quantifiable by the standard deviation of the system responses. Such an analysis allows to make a more weighted decision when choosing the uncertainty propagation technique because a quantified level of smoothness can tell us whether a system is capable or not of being approximated by surrogate techniques such as PCE. As part of this study, we have demonstrated on three different fluid flow cases that the PCE's ability to surrogate a system in relation to MC sampling, scales with the standard deviation of that system as a second-order power law, deviation from which suggests that a QoI is not entirely smooth. The devised relationship, whose applicability to other fluid flow systems still remains to be confirmed, may have the potential to help to make a more weighted decision when choosing the uncertainty propagation technique in fluid flow problems where the system dynamics are not so easily identifiable.

• The Engineering and Physical Sciences Research Council (EPSRC, UK) (Funding Data: 53889).

• Tata Steel via the Industrial CASE programme ((Funding Data: col 1178 GIPS 03132).