Reynolds Stress Perturbation for Epistemic Uncertainty Quantification of RANS Models Implemented in OpenFOAM

Reynolds-averaged Navier-Stokes (RANS) models are widely used for the simulation of engineering problems. The turbulent-viscosity hypothesis is a central assumption to achieve closures in this class of models. This assumption introduces structural or so-called epistemic uncertainty. Estimating that epistemic uncertainty is a promising approach towards improving the reliability of RANS simulations. In this study, we adopt a methodology to estimate the epistemic uncertainty by perturbing the Reynolds stress tensor. We focus on the perturbation of the turbulent kinetic energy and the eigenvalues separately. We first implement this methodology in the open source package OpenFOAM. Then, we apply this framework to the backward-facing step benchmark case and compare the results with the unperturbed RANS model, available direct numerical simulation data and available experimental data. It is shown that the perturbation of both parameters successfully estimate the region bounding the most accurate results.


Introduction
The motion of a fluid in the turbulent regime is fully described by the Navier-Stokes equations. A numerical solution encompassing all spatial and temporal scales is referred to as direct numerical simulation (DNS). Due to the significant computational cost of DNS, approximations like large-eddy simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) have been developed and are more widely used, in particular for practical engineering problems. RANS models use ensemble averages of the physical quantities, thereby decreasing the computational cost of the simulations. Due to the assumptions that are required to achieve closure, all RANS closures naturally introduce structural, or epistemic, uncertainty.
A systematic approach of the epistemic uncertainty quantification (EUQ) in RANS models, focusing on the Reynolds stress tensor, was first proposed by Emory et al. [1]. They introduced perturbation on the Reynolds stress tensor using the barycentric map proposed by Banerjee et al. [2], as a mean to visualize the degree of anisotropy. The same method was used later in a simulation by Gorle et al. [3] comparing RANS with LES results for an under-expanded jet in a supersonic cross flow. Another contribution, proposed by Gorle et al. [3] and further developed by Gorle and Iaccarino [4], was to introduce the perturbation not only in the momentum equations but also in the turbulent scalar fluxes in the scalar transport equation.
Gorle et al. [5] further developed the methodology with the goal of defining the uncertainty in the separation region. They introduced the idea of a marker, which identifies the regions of the

Epistemic Uncertainty
RANS models decompose the quantities into an averaged partū and a fluctuating part u such that Substituting (1) into the Navier-Stokes equations and applying Reynolds averaging leads to the RANS equations (see e.g., [10][11][12]) ∂ tūi +ū j ∂ jūi = ∂ ip + ν∂ j ∂ jūi − ∂ j R ij (2) and the Reynolds-averaged incompressibility constraint where t is time,ū i is the mean velocity in the i−direction, where i = 1, 2, 3 respectively represents the streamwise (x), wall-normal (y) and lateral (z) directions,p is the mean kinematic pressure, ν is the kinematic viscosity, and is the Reynolds stress tensor. The Reynolds stress tensor is an extra term that arises as a result of the Reynolds decomposition and contains the fluctuating parts of the velocity. The closure problem in RANS amounts to finding approximations for R ij that do not include the fluctuating part, since that information is not available in an actual RANS simulation. Yet, any model for R ij will be an approximation that introduces an additional epistemic uncertainty. Most popular RANS closure models are based on the Boussinesq turbulent-viscosity hypothesis (or eddy-viscosity hypothesis) that approximates the Reynolds stress tensor as a function of the stretching tensor (mean strain-rate tensor S ij = (∂ iūj + ∂ jūi )/2), the turbulent kinetic energy (k = u i u i /2) and the turbulent viscosity (effect of turbulent eddies on the flow ν t ) such that

Decomposition of the Reynolds Stress Tensor
The Reynolds stress tensor can be decomposed into components that establish its amplitude, shape and orientation. This is done by defining the turbulence anisotropy tensor a ij as Equation (6) can be normalized by dividing by 2k, yielding The eigenvalues of R ij are non-negative due to the condition of realizability established by Speziale et al. [13] and the Cauchy-Schwartz inequality holds for the off-diagonal components. Hence, the values are constrained within the intervals and In view of (7), the Reynolds stress tensor can be written as The eigendecomposition of the normalized anisotropy tensor yields (see e.g., [14]) where • Λ nl is a diagonal tensor containing the eigenvalues of the anisotropy tensor in an order such that v ij is a tensor containing the eigenvectors of the anisotropy tensor in the same order as the eigenvalues.

Perturbation of the Reynolds Stress Tensor
In Section 2.2, the Reynolds stress tensor is rearranged into components that define the amplitude (k), shape (Λ) and orientation (v). Perturbation of these components yields an estimate on the epistemic uncertainty. The perturbation yields where • v * = v + ∆v is the perturbation on the orientation of the Reynolds stresses, • Λ * = Λ + ∆Λ is the perturbation on the anisotropy of the Reynolds stresses, • and k * = [k/n k , n k k], is the amplitude of the perturbation of the turbulent kinetic energy. It is established as a range with a minimum and a maximum k * .
In this article, we focus on the perturbation of the turbulent kinetic energy k and the eigenvalues Λ. The perturbation of k is realized through the parameter n k ≥ 1 that determines the limits of perturbation. The maximum perturbation corresponds to k * = n k k and the minimum to k * = k/n k . The perturbation of the eigenvalues is realized using the barycentric map proposed by Banerjee et al. [2] and illustrated in Figure 1.
x 1c x 2c x 3c are functions of the eigenvalues associated with the Reynolds stress tensor. The limiting states are normalized such that The limiting state of the tensor can be represented in a two-dimensional coordinate system with coordinates or where B −1 x 1c = (2/3, −1/3, −1/3) T , B −1 x 2c = (1/6, 1/6, −1/3) T and B −1 x 3c = (0, 0, 0) T . Once the coordinates of the anisotropy tensor are located in the barycentric map, the perturbation is based on two parameters, the direction of the perturbation x (t) and the magnitude of perturbation δ B . The direction of the perturbation defines the chosen corner in the barycentric map and magnitude describes the distance of displacement to the chosen corner in the barycentric map in a range of [0, 1] as shown in Figure 1, The perturbed eigenvalues are calculated as [6] Here, it is noteworthy to mention that, in the RANS simulation of two-dimensional flows using eddy-viscosity models, at least one of the eigenvalues (λ l ) of the normalized anisotropy tensor (b ij ) is zero. Hence, all the points in the flow domain will be located along the line representing the plane shear flow (see Figure 1) [6].

Results and Analysis
We applied the previously explained methodology to the backward-facing step at Re h = 5100 (Re is dependent on the step height h, the inlet velocity u 0 and the kinematic viscosity ν, Re h = u 0 h/ν), for which DNS [15] and experimental data [9] are available. The backward facing step offers a simple two-dimensional case, similar to the canonical channel flow but adds a certain component of complexity in the physics right after the expansion appears.
The calculation was performed using the open source software OpenFOAM. The implementation of the Reynolds stress perturbation in OpenFOAM is explained in detail in the Appendix A. The boundary and initial conditions of the case were the same as the ones used by Le and Moin [15] and Jovic and Driver [9], including the mean velocity profile imposed at the inlet by Spalart [16]. They are summarized in Table 1. The k-ω SST model is chosen as RANS closure. The transport equations and their initial conditions were calculated following Ferziger and Peric [17] and Versteeg and Malalasekera [18]. Figure 2 shows the dimensions of the flow domain and its corresponding areas.  The topology and mesh resolution is similar to the one described in [19]. Three different mesh resolutions were considered to assess the sensitivity of the results to the grid. A relevant parameter to be taken into account, in this case, is the dimensionless wall distance y + which is defined as Keeping this value below one (y + < 1) allowed the simulation to capture the physics of the viscous sub-layer. Figure 3 shows the mesh topology. It can be seen that the grid gets finer when it approaches the bottom of the domain and also when it gets closer to the step, which is the most relevant part of the domain. Table 2 summarizes the topology and refinement used in each of the blocks of the domain. The mesh is designed such that the size of the first and last cells can be computed as a function of the grading and the number of cells. The grading is the ratio of the size of the last cell divided by the size of the first cell.  Table 2. Mesh refinement and topology for each of the three types of refinements applied. The Cells columns represent the amount of cells in each block in the x-and y-direction, respectively. The Grading columns represent the refinement in each block in the x-and y-direction, respectively.  Figure 4 illustrates the zoom into the corner where the step starts and the expansion begins. It can be seen that the size of the cells does not change sharply. It changes smoothly in both directions and it gets finer towards the lower wall. In Figure 5, all three refinements were compared for the pressure coefficient, friction coefficient, mean velocity and Reynolds stresses at x/h = 4, x/h = 6 and x/h = 10. The comparisons between different mesh refinements show the grid convergence for the intermediate and the fine grids. The rest of the work is based on the fine mesh described in Table 2. While the focus on the present article are the structural uncertainties in the RANS closures, it is worth mentioning that numerical uncertainties can make up a significant contribution to the overall model error. For a summary of different uncertainty contributions in computational fluid dynamics refer to the ASME standard [20] or more specifically to Roache [21] and Roache [22] in which the estimation of the numerical uncertainty associated with finite grid sizes is addressed. The following subsections show the behaviour of: pressure coefficient, friction coefficient, mean velocity profile and Reynolds stress computed at several positions in the domain (perturbed and unperturbed). The corresponding DNS and experimental data are also included for comparison.

Pressure Coefficient
The pressure coefficient was computed at the bottom of the domain throughout the stream-wise direction as where p o is the reference wall static pressure, p is the mean pressure computed at any location in the domain and u 0 is the upstream freestream reference velocity. We can see in Figure 6 that the results are bounded in a well defined area, thus decreasing the level of uncertainty. We can also see that a small amount of perturbation does not capture properly the physics of the problem. When the amount of perturbation increases, the bounds widen and they are able to cover better the discrepancy between the RANS model and the reference DNS and experimental data.

Friction Coefficient
The friction coefficient with τ w the wall shear stress was computed at the bottom of the domain, throughout the stream-wise direction. The results seen in Figure 7 show similar behavior as the one seen in Figure 6 with the exception that the turbulent kinetic energy perturbation for n k = 2 almost entirely defines the bounds of discrepancy between RANS and DNS and experimental results.

Mean Velocity in the x-Direction
The mean streamwise velocityū/u 0 throughout the y-axis was computed at x/h = 4, x/h = 6 and x/h = 10, respectively. The comparison between the models with perturbed kinetic energy, perturbed eigenvalue and the non perturbed model in Figure 8 show the pattern seen in previous parameters. The higher the amount of perturbation, the wider the bounds, thus the proportion of data covered in the bounds is higher. These results show the capability of the framework to capture the uncertainty in the simulation results.

Reynolds Shear Stress
Reynolds shear stress was computed vertically at x/h = 4, x/h = 6 and x/h = 10, respectively, and normalized with respect to the upstream reference velocity. As can be seen from Figure 9, there is some discrepancies between the unperturbed RANS model and the reference data from experiments and DNS. By introducing the structural uncertainties to the RANS, the model can bound the experimental data. It is also observed that the model uncertainty is not uniformly distributed in the domain. In some regions of the flow, the model uncertainty is larger due to the shortcoming of the model assumptions in capturing the physics of the flow. Also, it can be seen that increasing the amount of perturbation in both the magnitude and shape widen the gray regions that envelop the baseline results. For instance, increasing the amount of perturbation for turbulent kinetic energy from n k = 1.5 to n k = 2 widens the uncertainty bound (gray region) by a factor of about 2. A relatively similar trend is observed when the injected uncertainty into the shape of the stress tensor increases from δ B = 0.1 to δ B = 0.25.

Summary and Conclusions
In this article, we applied the framework, originally proposed by Emory et al. [6] for quantifying the structural uncertainties in RANS models. The quantification consisted in bounding the regions, where the most accurate results are certain to be located. This methodology focuses on the perturbation of the Reynolds stress tensor in the momentum equation. The Reynolds stress tensor is decomposed into components that represent the amplitude (turbulent kinetic energy), the shape (eigenvalues) and the orientation (eigenvectors). This investigation focused only on the influence of the amplitude and shape perturbations.
This perturbation is carried out by the turbulent kinetic energy and the eigenvalues of the Reynolds stress tensor separately and is applied to the backward-facing step case using the open-source software OpenFOAM. In this investigation, the following quantities were monitored: pressure and friction coefficients at the wall, mean streamwise velocity field and the Reynolds stress tensor. The investigation shows promising results. The results showed that for a specific amount of perturbation, the model provides a range of values where the physics of the case are well represented. An interesting feature of this research is that the perturbation of the turbulent kinetic energy shows results as good as the perturbation of the eigenvalues. The eigenvalues have well-established limits of the perturbation through the barycentric map, while the turbulent kinetic energy has not. This implies that the amount of the perturbation applied to the turbulent kinetic energy is not chosen systematically, thus it is left for future investigation.
In view of the results of this article, it is recommended for future work to focus on applying the perturbation to other parameters of the Reynolds stress such as the eigenvectors or a combination that optimizes the range that captures the most precise solution [7]. Recent studies [23,24] suggested that machine learning can be used as a powerful tool to predict discrepancies in the magnitude, anisotropy and orientation of the Reynolds stress tensor. It would be useful to continue applying this methodology to a new set of cases to monitor its behavior. simulations are run again until the results are converged. The corresponding flowchart is also plotted in Figure A4. x 1c x 2c x 3c x 1c x 2c x 3c x 1c x 2c x 3c (1) (2) (3) Figure A3. Locations of the chosen points shown in Figure A2 in the barycentric map. The points were perturbed towards the three limiting states and the amount of perturbation is δ B = 0.25. Figure A4. A flowchart describing the implementation of Reynolds stress perturbation in OpenFOAM.