{'fraction_non_alphanumeric': 0.02437851435335898, 'fraction_numerical': 0.006936223734832791, 'mean_word_length': 4.840445068596738, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantifying model form uncertainty in Reynolds-averaged turbulence models with Bayesian deep neural networks

Data-driven methods for improving turbulence modeling in Reynolds-Averaged Navier-Stokes (RANS) simulations have gained significant interest in the computational fluid dynamics community. Modern machine learning algorithms have opened up a new area of black-box turbulence models allowing for the tuning of RANS simulations to increase their predictive accuracy. While several data-driven turbulence models have been reported, the quantification of the uncertainties introduced has mostly been neglected. Uncertainty quantification for such data-driven models is essential since their predictive capability rapidly declines as they are tested for flow physics that deviate from that in the training data. In this work, we propose a novel data-driven framework that not only improves RANS predictions but also provides probabilistic bounds for fluid quantities such as velocity and pressure. The uncertainties capture both model form uncertainty as well as epistemic uncertainty induced by the limited training data. An invariant Bayesian deep neural network is used to predict the anisotropic tensor component of the Reynolds stress. This model is trained using Stein variational gradient decent algorithm. The computed uncertainty on the Reynolds stress is propagated to the quantities of interest by vanilla Monte Carlo simulation. Results are presented for two test cases that differ geometrically from the training flows at several different Reynolds numbers. The prediction enhancement of the data-driven model is discussed as well as the associated probabilistic bounds for flow properties of interest. Ultimately this framework allows for a quantitative measurement of model confidence and uncertainty quantification for flows in which no high-fidelity observations or prior knowledge is available.


Introduction
Over the past decade, with the exponential power increase of computational hardware, computational fluid dynamics has become an ever more predominate tool for fluid flow analysis. The Reynolds-averaged Navier-Stokes (RANS) equation provides an efficient method to compute time-averaged turbulent flow quantities making RANS solvers a frequently selected CFD method. However, it is common knowledge that RANS simulations can be highly inaccurate for a variety of flows due to the modeling of the Reynolds stress term [1]. Although over recent years Large Eddy Simulations (LES) or Direct Numerical Simulations (DNS) have become more accessible, these methods still remain out of the scope of practical engineering applications. For example, design and optimization tasks require repeated simulations with rapid turnaround time requirements for which RANS simulations are the choice modeling tool. Thus improving the accuracy of RANS simulations and providing measures of their predictive capability remains essential for the CFD community.
Turbulence models seek to resolve the closure problem that is brought about from the time averaging of the Navier-Stokes equations. While CFD and computational technology has made significant strides over the past decade, turbulence models have largely become stagnate with many of today's most popular models being developed over two decades ago. Many of the most widely used turbulence models use the Boussinesq assumption as the theoretical foundation combined with a set of parameters that are described through one or more transport equations. In general, these turbulence models can be broken down into families based off the number of additional partial differential equations they introduce into the system. For example, the Sparlart-Allmaras model [2] belongs to the family of single equation models. While the Sparlart-Allmaras has been proven to be useful for several aerodynamic related flows [3], its very general structure severely limits the range of flows that it is applicable. In the two-equation family, models such as the k-model [4,5] and the k-ω model [6] provide better modeling for a much larger set of flows even though their limitations are well known. In all the aforementioned models, a set of empirically found constants are used for model-calibration thus resulting in potentially poor performance for flows that were not considered in the calibration process. This combined with empirical modeling of specific transport equations, such as the equation, result in a significant source of model form uncertainty. While many have proposed more complex approaches such as using different turbulence models for different regions of the flow [7] or using a turbulence model with additional transport equations [8], these methods still rely heavily on empirical tuning and calibration. Thus model form uncertainty introduced by turbulence models continues to be one of the largest sources of uncertainty in RANS simulations.
This work aims to improve turbulence modeling for RANS simulations using machine learning techniques that also allow us to quantify the underlying model error. While the use of machine learning methods in CFD simulations can be traced back to over a decade ago [9], recently there has been a new wave of integrating innovative machine learning algorithms to quantify and improve the accuracy of CFD simulations. Earlier work in quantifying the uncertainty and calibration of turbulence models focused on treating model parameters as random variables and sampling via Monte Carlo to obtain a predictive distribution of outcomes [10,11]. Rather than constraining oneself to a specific model, an alternative approach was to directly perturb components of the anisotropy term of the Reynolds stress [12]. Lately, the use of machine learning models has been shown to provide an efficient alternative to direct sampling. In general, the integration of machine learning with turbulence models can be broken down into three different approaches: modeling the anistropic term of the Reynolds stress directly, modeling the coefficients of turbulence models and modeling new terms in the turbulence model. Tracey et al. [13] explored the use of kernel regression to model the eigenvalues of the anisotropic term of the Reynolds stress. Later, Tracey et al. [14] used a single layer neural network to predict a source term in the Sparlart-Allmaras turbulence model. Similarly, Signh et al. [15] have used neural networks to introduce a functional corrective term to the source term of the Sparlart-Allmaras turbulent model for predicting various quantities over airfoils. Zhang et al. [16] investigated the use of neural networks and Gaussian processes to model a correction term introduced to the turbulence model. Ling et al. [17] considered deep neural networks to predict the anisotropic tensor using a neural network structure with embedded invariance [18]. Ling et al. [19] additionally proposed using random forests to improve RANS predictions for a flow with a jet in a cross flow.
While the above works have managed to improve the accuracy of RANS simulations, uncertainty quantification has largely been ignored. Arguably, the integration of black box machine learning models increases the importance of uncertainty quantification in the context of quantifying the error of the improved turbulence model but also quantifying the uncertainty of the machine learning model itself. This is largely due to the significant prediction degradation of these proposed machine learning models for flows that vary from the training data in either fluid properties or geometry [13,17]. Past literature has clearly shown that data-driven methods are not exempt from the conflicting objectives of predictive accuracy versus flow versatility seen in traditional turbulence modeling.
Several works have taken steps towards using machine learning to provide uncertainty quantification analysis of RANS simulations. For example, Xiao et al. [20] proposed a Bayesian data-driven methodology that uses a set of high-fidelity observations to iteratively tune an ensemble of Reynolds-stress fields and other quantities of interest. While proven to work well with for even sparse observational data, this work is limited to a single flow with which the machine learning model was trained explicitly on. Wu et al. [21] used the Mahalanobis distance and kernel density estimation to formulate a method to predict the confidence of a data-driven model for a given flow. While this allows the potential identification of regions of less confidence after training, this is is limited to the prediction of the anisotropic stress and fails to provide any true probabilistic bounds.
For machine learning methods to be a practical tool for reliably tuning RANS turbulence models, transferability to flows with different geometries and fluid properties is important. Additionally, quantifying the model uncertainty is critical for assessing both the accuracy and confidence of the machine learning model and of the resulting predicted quantities of interest.
The novelty of our work is the use of a data-driven model with a fully-Bayesian Deep Learning framework to provide the means of improving the accuracy of RANS simulations and allow for the quantification of the model-form uncertainty arising in the model. This uncertainty is then propagated to the quantities of interest, such as pressure and velocity. The focus of our work will not be on application on flows that are the same or similar to those in the training set, but rather to flows defined by different geometries and fluid properties. We aim to take a much more practical and expansive view of using these innovative machine learning models for improved turbulence modeling. The specific novel contributions of this work are fourfold: (a) the use of a Bayesian deep neural network as a model to predict a tuned Reynolds stress field, (b) introducing a stochastic data-driven RANS algorithm that allows us to calculate probabilistic bounds for any flow field quantity, (c) assessment of the data-driven model on flows that are geometrically different from the training simulations and (d) comparison of both performance and confidence of the datadriven model across several Reynolds numbers. This paper is structured as the following: In Section 2, we review the governing equations and motivation for this work. In Section 3, the proposed data-driven frame work is discussed in detail. We discuss the invariant machine learning model in Section 3.1, its extension to the Bayesian paradigm in Section 3.2 and the stochastic data-driven RANS methodology to propagate uncertainty from the Bayesian datadriven model to quantities of interest in Section 3.3. In Section 4, various implementation details are reviewed including information regarding flow data used, training techniques and integration in the selected CFD solver. Section 5 details the results of applying this model to the two test flows at three different Reynolds numbers. Results for a flow over a backwards step and over as wall mounted cube are presented in Section 5.1 and Section 5.2, respectively. Finally discussion and conclusion are provided in Section 6.

Governing Equations
As previously mentioned, the difficulty of RANS is the fundamental closure problem that is introduced when the Navier-Stokes equations are averaged with respect to time. The RANS momentum equation is as follows: As always, the challenge is to close this equation by approximating the Reynolds stress (R-S) term u i u j . u i indicates a fluctuation velocity defined as u i (x, t) − u i (x, t) in which · indicates time-averaged or mean value. The turbulent viscosity theory, originally developed by Boussinesq [22], proposes a form of the R-S that is mathematically analogous to the stress-strain rate of a Newtonian fluid: where k, ν t are the turbulent kinetic energy (TKE) and turbulent viscosity, respectively. Assuming that the flow is incompressible results in the following: This representation is used not because of its accuracy but instead due to the simplifications that result when it is substituted into the RANS equation. This form is known as the Boussinesq eddy viscosity assumption.

RANS Turbulence Models
The context of this work is focused on the k − turbulence model [4,23,24] which is the most commonly used closure model for RANS simulations to date [1]. Starting with the Boussinesq eddy viscosity assumption, the k − model approximates the effective viscosity ν t in terms of the turbulent kinetic energy k and the turbulent dissipation rate with the R-S given as follows: where C µ is one of five model constants of this model. Through manipulation of the Navier-Stokes equations, the kinetic energy can be derived precisely for the case of high Reynolds number. On the other hand, the standard transport equation for the turbulent dissipation, , should be thought of as an empirical fit [1]. For this work, we will use the standard k − model for fully turbulent, incompressible flow [6]: The five constants C µ , C 1 , C 2 , σ k , σ are tunable parameters whose optimal values depend on the flow under consideration. For this work, we use the values originally proposed by Launder et al. [23] obtained by data fitting over various turbulent flows: The advantages of the k − model are its numerical robustness, computational efficiency, easy implementation and general validity for fully-turbulent flows. However, with this versatility comes some significant drawbacks including poor accuracy for complex fluid flows, and for problems with flow separation and sharp pressure gradients [7,25]. Core assumptions such as the formulation of the turbulent dissipation equations, the turbulent model constants and even the Boussinesq approximation provide large sources of uncertainty for the k − model. Converged simulations using the k − model with the parameters discussed above will be referred to as baseline RANS simulations. Ultimately, we seek to improve the prediction of a baseline simulation through the proposed data-driven framework.

Data-Driven Framework
In this work, our goal is to introduce a data-driven model to increase the accuracy of a given RANS simulation and to provide uncertainty bounds for quantities of interest thus capturing the error of the turbulence model. The proposed framework is illustrated in Fig. 1 which, in a broad sense, shares similar characteristics to earlier works on data-driven turbulence models [17,20,15]. However, we introduce several novel modifications to the process. We break this framework down into two key phases: the training of a model using a set of pre-existing flow data and the prediction stage for which the model is sampled from to produce fluid flow responses.
The training data that is driving our model is a small library of different fluid flows that attempt to capture different fluid physics. Ideally each training flow should bring new information for the model to learn thus increasing its potential predictive capability. For each unique flow, there is a low-fidelity control RANS solution and a time-averaged high-fidelity LES solution. The objective of this model is to learn the mapping from some baseline RANS flow input information to a turbulent property yielding an improved R-S field matching that of the corresponding high-fidelity simulation. This turbulent property could be tuned model coefficients, model correction terms or components of the R-S directly. For the scope of this work, we will focus on modeling the R-S tensor directly but this framework can extend to other approaches. An error or loss function that quantifies the discrepancy between the predicted R-S and the true high-fidelity field is used to update the model in an iterative process. We select a Bayesian neural network to serve as this model. Its formulation is discussed in Section 3.1 with a Bayesian extension presented in Section 3.2. The methods and techniques used to train the model are outlined in Section 4.
Once the model has been trained, it can be used as a regression model to sample predicted R-S fields for a given reference RANS solution. This process starts with a baseline RANS simulation whose flow field will serve as the input into the calibrated model. From this model, a set of turbulent properties are sampled that correspond to a predicted high-fidelity representation of the R-S field. For each predicted field, an independent forward simulation is completed in which the R-S is held constant and the remaining state variables are relaxed around the predicted field from their baseline values to updated perturbed values. We refer to this process of executing an ensemble of forward simulations as stochastic data-driven RANS (SDD-RANS). The forward simulations for different samples of the R-S can then be used to compute statistical bounds for quantities of interest as discussed in Section 3.3.
Remark 1. While we choose to use LES as the high-fidelity method for the training data, this can be substituted with higher accuracy methods such as DNS or even a combination assuming the turbulent statistics are consistent. On the other hand, the training and testing RANS simulations should maintain consistency with regards to the selected baseline turbulence models since different turbulence models can lead to significantly different results. This would potentially cause conflicting training data (e.g. same input different outputs).

Invariant Neural Network
As previously mentioned, in the scope of this work, we will predict the R-S field directly by the anisotropic component shown in Eq. (2). This approach has been used by multiple earlier works [17,20]. These works consider an explicit representation of the R-S component in particular in the form of a constant field. The remaining fluid flow quantities (mean velocity and pressure) are then propagated forward by solving a numerical system around this constrained R-S field. Such an approach does not constrain our work to model specific assumptions. We note that explicit R-S approaches can often result in significant prediction error for for flows when the Reynolds number approaches 5000 and above [26,27]. Additionally, small errors in the predicted R-S field in an explicit representation can be amplified leading to instabilities. We accept this as an open problem and while alternative implicit approaches have been proposed [27], we leave the discussion of such methods to future works. Rather we show how our framework can provide a confidence estimation for when the selected method starts to fail.
We select a Bayesian neural network to map the baseline RANS flow to a highfidelity R-S field due to the impressive performance of neural networks for highdimensional supervised learning tasks [28]. For the underlying neural network model, we choose the neural network proposed by Ling et al. [17], illustrated in Fig. 2. This neural network predicts the anisotropic tensor of the R-S using the symmetric and antisymmetric tensor components of the velocity gradient tensor. Through use of tensor invariants, the neural network is able to achieve both Galilean invariance as well as invariance to coordinate transformations. This makes such a model attractive for predictions of flows that deviate in geometry from the training data. Here, we briefly review the fundamentals of this neural network for completeness of the presentation.
The theoretical foundation of this invariant neural network is the non-linear eddy viscosity model developed in [29]. In this model, the anisotropic tensor of the R-S is expressed as a function a(s, ω) of the normalized mean rate-of-strain tensor s and rotation tensor ω: where both s and ω are scaled by the TKE and turbulent dissipation. Through application of the Cayley-Hamilton theorem, it can be shown that every second-order anisotropic tensor can be expressed in the following form: ITr s 2 , ITr(ω 2 ), ITr(s 2 ω 2 ), where T λ is one of 10 independent, symmetric tensor functions and G λ are the respective coefficients in the linear model which can be each expressed as functions of the five invariants I 1 , · · · , I 5 . For complete details on the invariants, tensor functions and the derivation, we refer the reader to [29]. The neural network model proposed by Ling et al. [17] models the anisotropic term by using the linear combination in Eq. (13). As illustrated in Fig. 2, rather than using the components of the symmetric and antisymmetric tensors (s and ω) directly, the invariants and tensor basis functions in Eq. (14) and Eq. (15) are used instead. To enforce invariance to coordinate transformations, the neural network is used to learn the tensor basis coefficients G λ which are functions of the five invariants in Eq. (14). These predicted coefficients, G λ , can then be used with the tensor basis functions, T λ , to produce the anisotropic tensor a. Thus while the model predicts the anisotropic tensor a given the symmetric and antisymmetric tensors of the velocity gradient, the function of the basis coefficients with respect to the five invariants is what is being learned. If a model uses inputs with specific invariant properties, the model has the same invariance properties as well [30]. This allows the neural network to be (a) Galilean invariant due to the use of the rate-of-strain and rotation tensors which are functions of the velocity gradient; and (b) invariant to coordinate transformations through the use of the invariant inputs I i . Additionally, since this eddy-viscosity model is the most general formulation, this neural network model does not share any of the limitations of other simpler eddy viscosity models that place restrictive assumptions on the form of the anisotropic term.
Remark 2. This neural network formulation is trained on entirely local (point-wise) information. The key advantage of a spatially local model is that it extends very easily to training flow data provided on non-uniform meshes which are essential in practical CFD simulations. Approaches such as convolution neural networks require training data on a uniform mesh following an image-to-image like regression approach [31]. However, similarly to turbulent eddy-viscosity models, this approach implies that the R-S mean convection D u i u j /Dt is governed entirely by local quantities (e.g. k, , ∂ u i /∂x j ). This is a questionable assumption for flows that exhibit strong inhomogeneity [1]. A model that incorporates spatial correlations would likely be more descriptive, physically robust and easier to train. Figure 2: Invariant, fully connected (some connections omitted for clarity), neural network architecture proposed by Ling et al. [17]. The circles indicate scalar values and the rectangles represent 3 × 3 second-order tensors.

Bayesian Neural Network
Traditionally neural networks are not designed to yield predictive statistics, however multiple recent works explore Bayesian reformulations of neural networks. Older techniques for obtaining Bayesian statistics include the placement of distributions over network weights and sampling with Monte-Carlo methods to approximate statistical bounds [32,33] as well as ensemble methods [34,35]. More recently, methods involving stochastic variational inference have brought a new wave a Bayesian neural network techniques [36,37,38,39]. In this work, we choose to use Stein's variational gradient decent (SVGD) recently proposed by Liu et al. [39,40] that approximates a variational distribution through a set of particles. SVGD is a non-parametric algorithm of similar form as standard gradient decent.
We follow closely the work of Zhu and Zabaras [31] in which SVGD is successfully applied to deep convolutional neural networks used for surrogate modeling. For the invariant neural network architecture discussed previously, we will use the following representation: where the input x = {s, ω} consists of the strain and rotation tensors s and ω along with the neural networks parameters w which include weights and biases. In the equation above, we have defined the neural network model as a function that has absorbed the calculation of the invariants, tensor basis functions and the linear combination detailed in Eqs. (13)-(15). Thus we will refer to the function f as the invariant neural network model for convenience. We wish to treat the neural network's learnable parameters, w, as random variables. Due to the potentially large number of weights in a fully connected neural network, we assume that the weights have a probability density function of a fully factorizable zero mean Gaussian and Gamma-distributed precision α: Additionally, output-wise noise is added onto the predicted output to represent inherent uncertainty within the model's formulation or uncertainty that cannot be reduced with more training data. This results in an additional noise term to the likelihood function of the neural network. We assume that the noise takes the form of a zero mean Gaussian with precision β that is a learnable parameter. We assume that β is Gamma distributed and that the training data have small noise variance, where a 1 = 2, b 1 = 2 · 10 −6 . To optimize the parameters in the neural network, SVGD minimizes the KL divergence between the true parameter posterior, p(w|D), for which p(w|D) is the unnormalized posterior and K is the log normalization constant that is not required to be computed during optimization. For the given neural network, we prescribe a Gaussian likelihood function and the priors discussed previously: Rather than attempting to recover a parametric form of the variational distribution, SVGD describes q(w) by a particle approximation. Namely, a set of N neural networks each representing a particle . Thus the objective is now for the empirical probability measure, µ N , to converge in distribution towards the true measure of the posterior ν, To minimize the KL divergence, we assume that q(w) is from a class of distributions that can be obtained through a set of smooth transforms. A small perturbation function, resembling that of standard gradient decent, is used to update the particles: where is the step size and φ(w j ) is the direction of the update that lies in a function space F for the t-th iteration. It is now a matter of finding the optimal direction to permute the particles in which should be chosen such that the gradient of the KL divergence is maximized, where T µ denotes the updated empirical measure of the particles. Liu et al.. [39] identify connections between the function φ and Stein's method and show that: in which T p is known as the Stein's operator. Assuming that this function space F lies within a reproducing kernel Hilbert space H with positive kernel k(w, w ), the optimal direction has the closed form: In this work, we choose to use the standard radial basis function kernel for k(w, w ). This formulation results in a simple two line update procedure in which the optimal decent direction for all particles is calculated with Eq. (29) and then updated by Eq. (26). Monte Carlo approximations can then be used to find the predictive mean and variance: For complete details on SVGD, we direct the reader to the original work by Liu et al. [39,40] along with the work of Zhu and Zabaras [31].

Uncertainty Quantification with SDD-RANS
We now wish to propagate this uncertainty obtained for the anisotropic term to the fluid properties such as pressure or velocity. We use a stochastic system approach for which the model parameters in a system of PDEs are considered as random variables. This methodology has been used extensively in the past for model calibration, prediction and selection [41,42,43,10]. Consider a dynamical system defined by the model output q (θ, u(θ)), and the system response R I which is dependent on some information I. u(θ) are state variables that evolve with the dynamical system and θ represents a set of model parameters with probability density p I (θ). Using the theorem of total probability, the expected response of the system can be written as: Traditionally the true form distribution p I (θ) from which the model parameters are sampled from is largely not known. However, under the assumption that samples can be drawn from the parameter distribution, the expected response as well as the respective variance can be approximated by vanilla Monte Carlo simulation (MCS): To extend this to the problem of interest and motivate SDD-RANS, let us consider the system of functions q as the RANS equations and the state variables u(θ) to be the fluid's velocity, pressure and all other derived properties. As previously discussed, we will be taking an explicit representation of the tuned R-S in which a tuned R-S field is predicted and held constant while the other state variables are propagated forward. Thus we are able to view a predicted R-S field, more specifically the anisotropic tensor a * , as a model parameter θ. Rather than sampling the anisotropic term directly as the stochastic parameter in the MCS, we instead recall the following representation of the invariant neural network model in Eq. (16): In the SVGD formulation of the neural network, the weights w are in fact the stochastic parameters in the model. We can modify the MCS such that we sample the weights of the Bayesian neural network: The SVGD algorithm intrinsically provides samples of the marginal prior p(w i |D).
Given that SVGD uses a particle representation, each sample is a particle (or invariant neural network) used during training. This combination of using a Bayesian data-driven model with a stochastic model parameter is why we have name this process stochastic data-driven RANS (SDD-RANS). With SDD-RANS, we have opened up the ability to obtain sample statistics for all flow quantities through traditional MCS. This allows for the quantification of uncertainty regarding our data-driven model beyond the R-S itself.
Remark 3. Note that here we use a single neural network to calculate a R-S prediction of the entire flow domain. In principle the neural network's inputs are spatially independent between mesh nodes. This potentially allows for each point to be treated independently of its neighbors when predicting a * (i.e. neural network weights can be re-sampled at each point). However, using a single neural network for an entire field allows for a much more spatially smooth result which is an important for simulation stability and convergence.

Framework Implementation
We use this Bayesian framework in the system of RANS equations by setting the R-S term as the stochastic parameter that is sampled from the predictive distribution obtained through the Bayesian neural network. We summarize the offline training process: • The training data consist of both baseline RANS and high-fidelity data for a set of different flows that attempt to capture different flow types.
• The underlying machine learning model is an invariant neural network that uses local fluid quantities to predict the anisotropic term of the R-S.
• We extend this invariant model to the Bayesian paradigm by using SVGD in which a set of neural networks approximate the posterior p(w|D) by particle representation.
• The parameters in each particle (or neural network) are optimized by minimizing the KL divergence between a particle variational approximation and the posterior of the parameters.
• An iterative algorithm, resembling the form of standard gradient decent, updates the parameters of each particle until convergence.
With the Bayesian neural network trained, one can make predictions to new flows: • For the flow of interest, a base-line RANS solution is obtained and the corresponding invariants and tensor functions at each mesh point are calculated.
• Each neural network used during training with SVGD is used to predict a corresponding high-fidelity R-S field.
• For each predicted field, the R-S is then constrained to the predicted values and a forward execution of the constrained system updates the remaining state variables.
• An equivalent number of state variable samples are then obtained for which probabilistic bounds can be calculated.

CFD Methods
For obtaining the training and test flows, the open source CFD platform OpenFOAM (Open source Field Operation And Manipulation) [44,45] is used. OpenFOAM is a widely accepted CFD package that contains a vast number of solvers for incompressible, compressible and multi-phase flows along with pre-and post-processing utilities.
For the baseline RANS simulations the steady-state, incompressible solver simple-Foam was used which employs the semi-implicit method for pressure linked equations (SIMPLE) algorithm [46] to solve both the momentum and pressure equations. The high-fidelity LES simulations used the pimpleFoam transient solver that combines both the PISO (Pressure Implicit with Split Operator) [47] and SIMPLE algorithms to solve the pressure and momentum equations. The Smagorinsky subgrid-scale model [48] with Van-Driest style damping was used for all LES flows. Both the baseline RANS and LES domains are discretized by second-order methods. Each training and testing flow is outlined in Table 1, and all meshes are non uniform such that the mesh density increases around the feature of interest. All simulations were run with a CFL number below 0.3 for numerical accuracy.

Machine Learning Implementation
To train the neural network, the Python machine learning library PyTorch [58] was used. The software and data used in this work would become available to the public upon publication of the paper at https://github.com/cics-nd/rans-uncertainty. The details of the network architecture used are given in Table 3. The specific number of hidden layers and hidden nodes can likely be varied to yield improved results, however the following was found to work adequately for the problem of interest. The Leaky Rectifier function was used as the activation function as opposed to the standard Rectifier function to prevent too many nodes from becoming zero during training. Additionally the number of nodes in the hidden layers is tapered at the end of the network to prevent weights from being too small, which improved training performance. Compared to other potential network models, we found that the model selected originally by Ling et al. [17] proved to be exceptionally difficult to train. This is reflected in the original work by the extremely low learning rate used of 2.5e−6. During training we also found only very low learning rate could be used for training to be stable. To increase training performance and efficiency, we also used the following techniques: • To increase training efficiency rather than using every single mesh point, a subset of training points is selected. In this work, we use only 1e4 total training points that are evenly distributed among all training flows (i.e. 2e3 points for each of the five test flows). Although the entire flow domain does contain relevant information, in practice most flows contain a large bulk region and different training flows may contain similar flow structures making many training points redundant. The training points selected are those that have the largest predictive variance for the anisotropic component of the R-S tensor computed in Eq. (33).
• Points with the largest variance will change as training progresses, thus every 15 epochs the variance for all flow fields is calculated and the training points updated which is illustrated in Fig. 3. While this causes training to be more noisy, re-sampling helps prevent the model from focusing too much on the bulk flow and work more in the areas at which the prediction quality is poorer. Without re-sampling and simply using ever data point available, we often found that the neural network would over smooth the entire prediction field in an attempt to cater towards the large number of data points in the bulk region.
• Training points are shuffled randomly and mini-batched every epoch such that data from multiple flows can reside in a single mini-batch. This avoids the network from learning the physics of each flow independently.
• Since the invariant inputs to the neural network tended to vary strongly in magnitude, they are re-scaled by a log operation to improve training. In addition, the tensor basis functions were normalized by the L 2 norm of the matrix: in which sgn is the signum function and is a tunable parameter which is set to be 1.0.
To quantify the training accuracy, we use a testing data set of 2000 randomly selected data points from the two test flows in Table 2. For both the test and training data we plot the mean squared prediction error (MSPE) along with the mean negative log probability (MNLP) for each training flow. The sharp discontinuities are due to the re-sampling of training points every 15 epochs. The training process on a single NVIDIA P100 GPU took approximately 4.5 wall-clock hours.

Constrained R-S Simulation
To integrate the sampled R-S field into OpenFOAM, a small modification is made to the simpleFoam such that the R-S is now a constant field in the momentum RANS equation. Since the R-S field is held constant, the calculation of the TKE and turbulent dissipation is no longer needed. An important issue is the handling of boundary conditions. This includes the treatment of domain boundaries as well as areas in which wall functions may be used. We address these issues using two different methods. First, the use of the baseline RANS TKE as a scaling factor of the anisotropic term shown in Eq. (11) allows for many turbulent boundary conditions to be satisfied. Second, to address areas in which wall functions may be used, we take inspiration from hybrid LES/RANS methods and introduce a blending function proposed by Xiao et al. [60]: where a dd is the data-driven prediction of the anisotropic tensor, a rans is the baseline RANS anisotropic tensor, λ = k/ , d is the distance from the wall and α 1 is a tunable parameter. This function allows a smooth transition between the use of the baseline R-S near the wall and the data-driven prediction in the bulk flow. In the original work, it is suggested that the selection of α 1 be a value that achieves Γ = 0.5 somewhere is the log region. We found the value of 0.05 worked well for our test cases. As previously mentioned, often for larger Reynolds numbers explicit R-S approaches can result in convergence difficulty for the propagated fluid properties. Ultimately this is not a large concern of ours since we view these new fluid fields simply as perturbations towards a more accurate solution. Thus for each test case, the momentum and pressure equations are propagated forward for the same number of time-steps as the baseline simulation. As we will show in the subsequent sections, numerical instabilities in the predicted field allow for a good indicator of the confidence we have in the data-driven turbulence model.

Numerical Results
The use of data-driven models for test simulations whose domain is similar or identical to the training data is a frequent occurrence in the literature but does not correctly assess a data-driven models performance. Since our selected neural network has already been shown to work adequately for similar flows in [17], our test cases are selected to deviate significantly from the training flows in both flow geometry and Reynolds number. We have selected the two test flows detailed in Table 2: flow over a backwards step and flow around a wall mounted cube both at three different Reynolds numbers. The geometry for each flow can be seen in Fig. 4. For the first test flow, we focus on how this framework can provide a confidence measurement and indicate when a data-driven simulation can be trusted. In the second test flow, we focus more rigorously on how our framework can be used for UQ for the fluid quantities of interest. Ultimately, we wish to assess the predictive performance of the data-driven model for these geometrically different flows and use the proposed stochastic data-driven RANS algorithm to calculate probabilistic bounds on flow state variables by conducting UQ on the data-driven turbulence model.

Backwards Step
In the first test case, we select a backwards facing step at three different Reynolds numbers. As illustrated in Fig. 4a, this flow features a constant velocity inlet channel followed by a backwards facing step of height h. For this flow, we select the inlet channel to be the same height as the step. The no-slip walls are on both the top and bottom faces and the z direction is periodic. On the x − y plane, we place the origin at the corner of the step. The flow features of interest are the recirculating regions that appear not only directly after the step but also on the upper channel wall down stream which is seen in the LES simulations in Figs of the anisotropic terms of both the baseline RANS and the expected values of the Bayesian neural network are shown in Fig. 5 for Reynolds number 500. The region just downstream of the step, where the recirculation regions occur, is plotted. While it is clear that both RANS and the data-driven model vary heavily from the true LES solution, we note that the data-driven approach is able to improve the prediction over the baseline RANS. The data-driven model is able to correct the sign for the tensor terms a 11 and a 22 as well as predict non-zero values for the a 33 term. We note that the baseline RANS predicts a zero field for the a 33 term due to the linear eddy-viscosity assumption. The stream-wise velocity contours of the flow for the LES, baseline RANS and the expected velocity prediction of stochastic data-driven RANS (SDD-RANS) are depicted in Figs. 6-8. To keep plot labels uncluttered, we refer to the expected values as just SDD-RANS. For each Reynolds number, we also show the standard deviation field corresponding to the samples from SDD-RANS. We look first at the lowest Reynolds number of 500 for which the model produced the best prediction for. While being the Reynolds number farthest from the training data in Table 1, the stochastic model was able to successfully predict the appearance of the second recirculation region The baseline RANS simulation only predicted a single the eddy behind the step. The second eddy however is predicted to be much closer to the inlet than the target LES solution. We presume that this is likely due to the network being trained on higher Reynolds number flows, and thus over predicting the magnitude of the turbulence in the flow. This over estimation in magnitude is clearly seen in Fig. 5. Figure 6: Normalized stream-wise mean velocity contours for Reynolds number 500. The top is the LES solution, below is the baseline RANS prediction followed by the data-driven mean field. Lastly is the standard deviation field of the data-driven prediction.
As the Reynolds number is increased, a very clear degradation in the predictive performance of SDD-RANS performance is seen. For both Reynolds number 2500 and 5000, SDD-RANS does not yield any prediction improvement over the baseline RANS simulation. As briefly mentioned before, this is a common observation for explicit R-S treatments. While the training of the neural network model could be improved to yield better estimates of the R-S, this does not necessarily guarantee an improved prediction [26]. In the high Reynolds number test cases (see , the use of SDD-RANS versus a standard deterministic data-driven approach is extremely beneficial since SDD-RANS can provide a measure of variance for QoI. As the R-S is increased, minor deviations in the anisotropic term are amplified resulting in potentially starkly different flow predictions [27]. SDD-RANS accurately captures these phenemona. As the Reynolds number increases, so does the standard deviation indicating a lack of model confidence. This is most clearly seen in Reynolds number 5000 (Fig. 8) where the magnitude of the stream-wise velocity standard deviation in some areas is almost the magnitude of the inlet velocity. The extremely high variance provides an excellent indicator in this case for when the model is likely to under perform. With a deterministic data-driven model such indicators would not be present allowing for no interpretable information on prediction confidence without observed high-fidelity data. Figure 8: Normalized stream-wise mean velocity contours for Reynolds number 5000. The top is the LES solution, below is the baseline RANS prediction followed by the data-driven mean field. Lastly is the standard deviation field of the data-driven prediction.
Another beneficial outcome of SDD-RANS variance estimation is the ability to identify regions of the flow that are likely to be less accurate. To illustrate this feature, the predicted stream-wise velocity is plotted against the true LES velocity in Fig. 9. In general, as the points deviate from the ideal line the variance associated with that local point also increases. This is most clearly seen with Reynolds number 500, where SDD-RANS is able to indicate the precise regions that deviate the greatest from the high-fidelity result.

Wall Mounted Cube
The second test case is a flow around a wall mounted cube with height h as shown in Fig. 4b. Unlike the majority of the flows that have been tested by data-driven models in the literature [13,15,16,17,20] as well as our training flows, this test flow contains an obstacle that is not semi-infinite. This means that flow with this geometry cannot be modeled by a two-dimensional RANS simulation as was the case for all previously considered flows. Additionally, similar to the backwards step, none of our training flows contain a geometry that is similar to this. As a result, we consider this flow an excellent test to investigate the limits of SDD-RANS in generalizing to a true 3D test case.
The set-up of this flow consists of an uniform inlet velocity and two channel walls normal to the y-axis. The cube is placed slightly down stream of the inlet. The features of interest is primarily the recirculation region behind the cube itself. Additionally, as the Reynolds number increases, flow separation occurs on the sides of the cube. As will be shown in the subsequent figures, this flow separation is often non-existent for the baseline RANS predictions. Even though the mean flow is symmetrical about the x − y plane in the middle of the channel (z = 3.5H), we simulate the entire cube in order to observe non-symmetrical behavior in predictions. The stream-wise velocity contours on the plane of symmetry for the LES, baseline RANS and the expected value of SDD-RANS are depicted below in Fig. 10. Again, for the sake of brevity and keeping the plot labels uncluttered, we will refer to the expected value of the data-driven model simply as SDD-RANS unless stated otherwise. Similar to the backwards step, as the Reynolds number increases, the standard deviation of the SDD-RANS prediction also increases. However, the magnitude of the variance is significantly smaller than that of the backwards step for higher Reynolds numbers which reflects the more reasonable predictions. Again the predicted streamwise velocities of both the baseline RANS and SDD-RANS are plotted against the high-fidelity LES solution in Fig. 11 where the SDD-RANS variance is again able to provide a a good measure for simulation reliability. To take a closer look at the performance of SDD-RANS, stream-wise velocity profiles are plotted for both Reynolds number 500 and 5000 in Figs. 12 and 13, respectively. For each, a plot containing the resulting SDD-RANS velocity field samples are shown as well as the respective standard deviation. As expected the variance in the velocity samples increases with the increased Reynolds number, however this only occurs in the recirculation region where the instantaneous flow is turbulent. In the bulk region above the recirculation zone, the variance remains small.  A significant improvement by SDD-RANS for both Reynolds numbers is the prediction of the detached flow on top of the cube (x/h = 4) which the baseline RANS fails to predict entirely. For both, improvements in the recirculation region prediction are present however the lower Reynolds number achieves the best predictions by SDD-RANS. The higher Reynolds numbers appear to have an over estimation of the length of the recirculation which are again regions of larger variance. Since this obstacle is not semi-infinite, the velocity contours for the horizontal plane at y = 0.5H are shown in Fig. 14. Similarly, velocity profiles are plotted for Reynolds number 2500 in Fig. 15. In general, we can see the same trends as previous results for which the variance increases with Reynolds number. We note that SDD-RANS is able to predict the presence of the detached flow on the side of the cube as well. While some asymmetry exists in the SDD-RANS predictions, the predictions overall retain a general symmetric profile. It is likely that increasing the number of samples of R-S fields would further improve the symmetry and smoothness of predictions. Figure 14: Normalized stream-wise mean velocity contours for Reynolds numbers 500, 2500 and 5000 on the y = 0.5H plane. The top is the time averaged LES solution, below is the baseline RANS prediction followed by the SDD-RANS expected velocity. Lastly is the standard deviation field of the data-driven prediction. This framework allows probabilistic bounds to be calculated for other fluid properties such as pressure, drag, shear stress, etc. For example, two pressure profiles along the face of the wall mounted cube are plotted in Figs. 16 and 17. In general we see that for most faces, SDD-RANS is able to provide an improved prediction compared to the baseline RANS. SDD-RANS partially corrects the unphysical pressure drop that occurs on the edge of the leading face in the baseline RANS simulation.  Figure 17: Normalized mean surface pressure profile on the plane y = 0.5H for Re= 5000.

Conclusions
As the CFD community continues to investigate the use of machine learning tools for data-driven modeling, the need to accurately quantify the induced uncertainty from the use of such models will be essential. In this work, we have presented a novel framework that allows for the quantification of such model form uncertainty. To satisfy invariant properties, we use the neural network architecture originally proposed by Ling et al. [17]. Using Stein's variational gradient decent and following the work of Zhu and Zabaras [31], we extended this invariant neural network model to a fully Bayesian deep neural network to allow us to compute the distribution of the anisotropic R-S stress tensor for a given baseline solution. To propagate the uncertainty of this model to fluid flow quantities of interest, a stochastic datadriven RANS algorithm is proposed that utilizes standard Monte Carlo simulation. The integrated framework was rigorously investigated on two flows to observe its generalization property.
In the presented implementation of this framework, we found that the invariant neural network used to model the anisotropic tensor proved difficult to train and yield satisfactory predictions for unseen flows and geometries. From our studies, we hypothesize that using just the five invariant inputs does not provide enough descriptive information to accurately map from the coarse to high-fidelity flow physics. Although the network contains desired invariant proprieties, other flow quantities would most likely need to be used as model inputs to improve the prediction quality. Thus a critical area to be investigated is the development of more accurate Reynolds stress representations identifying the important local-and non-local variables that influence its values. The potential use of spatial correlations and information at neighboring nodes (non-local models) may prove to be extremely beneficial. Such an approach, while difficult to implement for non-uniform grids, can be easily applied in the context of convolutional neural networks that are capable of mapping many high-dimensional inputs to multi-outputs of high-dimensionality.
SDD-RANS can provide a good indication of when the model is performing poorly, as seen in the backwards step test flow. The produced results will be of little use without a physically accurate calculation of the Reynolds stress. While a number of other models in the literature can yield much better training predictions, most of these models remain only useful to a small family of flows resembling those in the training dataset. The generalization property of these models remains an open problem for the data-driven community.
With improvements in the representation of the tuned Reynolds stress in the RANS equations, we believe that this framework can provide extremely beneficial information for data-driven models. While the most obvious application is its use to assess a given model's predictive confidence, the developed framework can also be used to identify locations of potentially lower accuracy. This could be useful for identifying areas that may require finer mesh resolutions or high-fidelity simulations. Additionally, the use of a Bayesian neural network not only provides a method for sampling model parameters but also predictive bounds for the quantities of interest. This can be extremely useful in cases where the training data is limited. Even though the use of the Bayesian neural network and SDD-RANS requires more computational time than deterministic data-driven approaches, we found that when compared to high-fidelity simulations the computational cost remains low.
Finally, the use of a Bayesian neural network opens up the potential of implementing experimental design techniques by investigating the impact of training data on the quality of the model's predictions. This can range from the assessment of a limited data case or how specific training flows at various Reynolds numbers impact a specific test case prediction.