Towards Uncertainty Quantification of LES and URANS for the Buoyancy-Driven Mixing Process between Two Miscible Fluids—Differentially Heated Cavity of Aspect Ratio 4

Abstract

Numerical simulations are subject to uncertainties due to the imprecise knowledge of physical properties, model parameters, as well as initial and boundary conditions. The assessment of these uncertainties is required for some applications. In the field of Computational Fluid Dynamics (CFD), the reliable prediction of hydrogen distribution and pressure build-up in nuclear reactor containment after a severe reactor accident is a representative application where the assessment of these uncertainties is of essential importance. The inital and boundary conditions that significantly influence the present buoyancy-driven flow are subject to uncertainties. Therefore, the aim is to investigate the propagation of uncertainties in input parameters to the results variables. As a basis for the examination of a representative reactor test containment, the investigations are initially carried out using the Differentially Heated Cavity (DHC) of aspect ratio 4 with $Ra=2\times {10}^9$ as a test case from the literature. This allows for gradual method development for guidelines to quantify the uncertainty of natural convection flows in large-scale industrial applications. A dual approach is applied, in which Large Eddy Simulation (LES) is used as reference for the Unsteady Reynolds-Averaged Navier–Stokes (URANS) computations. A methodology for the uncertainty quantification in engineering applications with a preceding mesh convergence study and sensitivity analysis is presented. By taking the LES as a reference, the results indicate that URANS is able to predict the underlying mixing process at $Ra=2\times {10}^9$ and the variability of the result variables due to parameter uncertainties.

1. Introduction

Numerical approximations and mathematical models are used to solve technical issues in computational science and engineering. Each of these numerical simulations is inherently subject to uncertainties, and in some applications, their consideration is essential. An important representative of this category in the field of Computational Fluid Dynamics (CFD) is the reliable prediction of hydrogen distribution and pressure build-up in nuclear reactor containment during an accident scenario, which is of utmost importance with regard to maintaining the integrity of the containment and thus preventing the release of radioactive substances [1,2]. An illustration of the primary coolant system of a nuclear power plant and of the prevailing processes during a severe nuclear accident scenario is shown in Figure 1. Figure 1c shows consecutive mechanisms within the reactor pressure vessel. Decay heat production of the core leads first to heating and then to evaporation of the water in the reactor pressure vessel. After the core is uncovered, the fuel begins to overheat because the heat transfer from the fuel to the steam is less than to the liquid water. The cladding tubes of the fuel rods consist of a zirconium alloy (Figure 1d), which reacts with the surrounding water vapour at temperatures above $1100 \mathrm{K}$. Zirconium oxide and hydrogen are formed within an exothermic reaction. Hence, additional energy is released, which heats the fuel rods and increases the reaction rate and hydrogen production. Furthermore, large quantities of hydrogen are additionally produced after the failure of the reactor pressure vessel during nuclear meltdown of the fuel when the concrete decomposes (Figure 1c). After pressure relief of the reactor cooling circuit or the melting through the reactor pressure vessel with the associated Molten Core–Concrete Interaction (MCCI), the hydrogen is released into the containment. An ignitable hydrogen air–steam gas mixture can be formed from the hydrogen, the atmospheric oxygen present in the containment, and the evaporated water. If the ignition of this hydrogen air–steam gas mixture occurs, the integrity of the containment is jeopardized since the design pressure of the containment may be exceeded. As a result, radioactive material would be released into the environment. Therefore, for safety, it is essential to understand the transport and mixing processes of hydrogen in the containment [3,4,5].

By means of Computational Fluid Dynamics (CFD), the flow and transport processes prevailing in the containment, which are primarly driven by buoyancy effects, can be calculated numerically. Numerous investigations have been conducted to predict the formation, stability, and erosion of a local flammable gas layer. For validation of the CFD models under application relevant conditions, the THAI test series is part of a comprehensive concept [6]. The acronym THAI stands for Thermal-hydraulics Hydrogen, Aerosols and Iodine. The THAI facility is a unique, technical-scale plant built for experimental research in nuclear reactor containment safety, which is often used for the validation of CFD codes.

In this work, the focus is on the buoyancy-driven mixing process of two miscible fluids. Since hydrogen is a highly flammable gas, helium is used in the experiments. The stratification and mixing phenomena are primarily governed by the density differences in the atmosphere. Both helium and hydrogen are characterized by a much smaller molecular mass compared to air. Since the difference in the molecular weights of these light gases causes only a very small contribution to the atmospheric density differences, similar volumetric concentrations are achieved. Different physical properties, such diffusivity, thermal conductivity, or heat capacity, are of minor impact. Therefore, the buoyancy-driven mixing process with helium instead of hydrogen can provide appropriate knowledge about the flow present in the containment [6].

In the corresponding THAI TH22 experiment, which examines the erosion by natural convection of a stratified helium layer, an unstable stratification is established by cooling the upper third and by heating two lower thirds of the cylindrical vessel walls. The experiment is comprise of three consecutive phases. The first is the initial phase in which natural circulation with air is established, as shown in Figure 2a. This is followed by the injection phase (Figure 2b), in which helium is injected into the reactor test containment through a nozzle and the circulation is suppressed. After the injection process has been completed and a stable stratification has formed, the erosion phase of the helium layer begins, as can be seen in Figure 2c. Once the circulation forms again, the stratification is eroded rapidly, as shown in Figure 2d.

On the basis of the natural convection experiments MISTRA NATHCO (CEA, France) and THAI TH22 (Becker Technologies, Germany), CFD analysis revealed a significant influence of the specified initial and boundary conditions (e.g., initial gas and wall temperatures) on the mixing process of air and helium [7]. However, the definition of initial and boundary conditions are subject to uncertainties in CFD simulations, since measurement errors in experiments occur or a lack of required input variables, which cannot be measured with justifiable effort, exists [8,9]. Thus, the configuration is not sufficiently specified precisely. This leads to results that also contain uncertainties. For their evaluation though, insufficient experience exists [10,11]. For the exploration and further development of Uncertainty Quantification (UQ) methods, the investigation of a single-phase mixing process by means of CFD is recommended [12].

Therefore, the aim of this work is to investigate the propagation of uncertainties in initial and boundary conditions as well as in material values for this application. As a basis for later examination of the THAI containment (Figure 3b), the investigations are initially carried out using a test case from the literature [13] with a superimposed mixing process, shown in Figure 3a. This allows for gradual development of guidelines to quantify the uncertainty of natural convection flows in large-scale industrial applications, which is the primary purpose of this work. Noting that uncertainty quantification requires a large number of Large Eddy Simulation (LES) runs to be performed, the deliberate choice of performing engineering LES requires compromises in terms of mesh resolution and numerical techniques. Appropriate parameters have been selected based on the results in Section 3.1 as a balance between accuracy, robustness, and computing efficiency. The geometry corresponds to a tall cavity with a height-to-width ratio of 4, for which the left wall has a higher temperature compared to the right side, and thus, natural convection is formed within the cavity. The configuration is generally referred to as a Differentially Heated Cavity (DHC). This simplified test case, which reflects similar physics as those expected inside the reactor containment, enables a simple and clear comparison of the results in the course of method development, since interpretation of the changing physical effects that result from varying the initial and boundary conditions is easier. Moreover, computations of the simplified test case are more resource-efficient and allow for different analysis methods to be applied.

A few related studies are reported in the literature and will be briefly summarized in the following. An uncertain single-phase turbulent mixing process in the presence of density gradients was investigated within a flow channel, where two co-flowing water streams were initially separated by the splitter plate. Sucrose was used to increase the density of one water stream to achieve density differences between the two streams of up to 10%. Various methodologies for uncertainty quantification were applied, and final assessment of the UQ methodologies was based on experimental results. Turbulent mixing in the presence of density gradients is a typical situation encountered in many reactor issues and is of practical interest to nuclear reactor safety [14,15,16,17,18].

Uncertain thermo-fluid flow within a cubic DHC was studied by Le Maître et al. in the Boussinesq [19] and non-Boussinesq limits [20] implementing the variable-density low-Mach-number equations using an intrusive Polynomial Chaos Expansion (PCE) approach. A steady laminar recirulating flow regime was investigated. The cold wall temperature was presumed as uncertain, and arising uncertainties in the mean velocity field were analysed. Rayleigh–Bénard convection (RBC) in the Boussinesq limit using PCE was analysed in [21]. For this configuration, the heated bottom wall was considered uncertain. In the present work, the application of uncertainty quantification is extended to the buoyancy-driven turbulent mixing process between two miscible fluids within the DHC. Integral result quantities are used to evaluate the transient mixing process. Developed methods can be transferred to other applications.

Uncertainty can be classified into aleatoric and epistemic uncertainties. Aleatoric uncertainty describes the natural randomness in a process. For example, turbulence is characterized by a three-dimensional flow field with an apparently randomly varying component in time and space. Hence, turbulent flow has aleatoric character. Epistemic uncertainty is defined as any lack of knowledge or information in any phase or activity of the modeling process, e.g., by inaccurate measurement or neglect of certain effects by models [22]. Unsteady Reynolds-Averaged Navier–Stokes (URANS) computations, due to the low computational effort, have a high degree of modelling and thus possibly high epistemic uncertainties, which makes it less reliable compared to Large Eddy Simulation (LES). In URANS, the complete turbulence spectrum is modeled and the Navier–Stokes equation is perfectly deterministic [23]. Uncertainties that originate from modeling of turbulence are epistemic since the uncertainty is based on incomplete knowledge [22]. Therefore, a dual approach is applied, in which the LES method is used as a reference for the URANS computations to be able to make statements on the uncertainty of the URANS modelling.

By direct comparison of the results between LES and URANS, it can be assessed to what extent a URANS-based sensitivity and uncertainty analysis can reflect the variability of a result variable as a function of an input variable. Furthermore, the parameters are revealed, which have major influences on the mixing process of two miscible fluids driven by buoyancy effects.

2. Numerical Methods and Flow Configuration

In Section 2.1, the governing equations are introduced while the flow configurations and numerical techniques are detailed in Section 2.2 and Section 2.3.

2.1. Governing Equations

The extensions of the computational domain are height H, width W, and depth D. The height-to-width aspect ratio and the width-to-depth aspect ratio are ${\phi }_{HW}=H/W$ and ${\phi }_{WD}=W/D$, respectively. Besides the initial and boundary conditions, the physics are completely defined by the Prandtl number $Pr=\nu /\alpha$; the height-to-width aspect ratio ${\phi }_{HW}$, which is equal to 4; and the Rayleigh number $Ra=\frac{g\beta \mathsf{\Delta }T{H}^3}{\nu \alpha }$, which is equal to $2\times {10}^9$. Periodic flow is assumed in the third spatial direction. For the mesh convergence study, the DHC was considered filled with a single incompressible Newtonian viscous fluid of kinematic viscosity $\nu$ and thermal diffusivity $\alpha$. The Prandtl number $Pr=0.71$ corresponds to air. The modelled URANS/LES equations of an incompressible fluid for continuity, momentum, and temperature considering the Boussinesq approximation are

$$\begin{array}{cc}\hfill & \nabla ·\mathbf{u}=0,\hfill \\ \hfill & \frac{\partial \mathbf{u}}{\partial t}+\nabla ·\left(\mathbf{u}\mathbf{u}\right)=-\frac{1}{{\rho }_0}·\left(\nabla p-\rho \mathbf{g}\right)+\nabla ·\left(2{\nu }_{eff}D\left(\mathbf{u}\right)\right),\hfill \\ \hfill & \frac{\partial T}{\partial t}+\nabla ·\left(\mathbf{u}T\right)=\nabla ·\left({\alpha }_{eff}\left(\nabla T\right)\right),\hfill \end{array}$$

(1)

where $\mathbf{u}$ is the velocity field, $\rho$ is the density field, p is the static pressure field, T is the temperature field, and $\mathbf{g}=\left(0,g,0\right)$ is the gravitational acceleration vector. The density $\rho$ in the gravitational term is expressed by the linear function of the temperature T: $\rho \approx {\rho }_{\beta }\left(1-\beta \left(T-T_{\beta }\right)\right)$. $\beta$ is the volumetric thermal expansion coefficient. We denote the reference density by ${\rho }_{\beta }$ at the reference temperature $T_{\beta }$. The rate of strain tensor is defined as $D\left(\mathbf{u}\right)=\frac{1}{2}\left(\nabla ·\mathbf{u}+{\left(\nabla ·\mathbf{u}\right)}^T\right)$. The effective kinematic viscosity ${\nu }_{eff}$ is the sum of the molecular and turbulent or subgrid-scale viscosity regarding URANS and LES, respectively. By using the gradient flux approach with the turbulent Prandtl number $P{r}_t=0.85$, the effective thermal diffusivity ${\alpha }_{eff}$ results as the following sum of laminar and turbulent or subgrid-scale thermal diffusivities: ${\alpha }_{eff}=\frac{\nu }{Pr}+\frac{{\nu }_t}{P{r}_t}$.

For the sensitivity analysis and uncertainty quantification, we considered a low Mach number flow of two Newtonian viscous fluids in the cavity. The Prandtl numbers $Pr=0.71$ and $Pr=0.66$ used correspond to air and helium, respectively. The material values of air are used to define the Rayleigh number $Ra=\frac{g\beta \mathsf{\Delta }T{H}^3\rho }{\mu \alpha }$. The continuity, momentum, energy, and mass transfer equations take the form

$$\begin{array}{cc}\hfill & \frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{u}\right)=0,\hfill \\ \hfill & \frac{\partial \rho \mathbf{u}}{\partial t}+\nabla ·\left(\rho \mathbf{u}\mathbf{u}\right)=-\nabla p+\rho \mathbf{g}+\nabla ·\left(2{\mu }_{eff}D\left(\mathbf{u}\right)\right)-\nabla ·\left(\frac{2}{3}{\mu }_{eff}\left(\nabla ·\mathbf{u}\right)\right),\hfill \\ \hfill & \frac{\partial \rho h}{\partial t}+\frac{\partial \rho K}{\partial t}+\nabla ·\left(\rho \mathbf{u}h\right)+\nabla ·\left(\rho \mathbf{u}K\right)-\frac{\partial p}{\partial t}=\nabla ·\left(\rho {\alpha }_{eff}\left(\nabla h\right)\right)+\rho \mathbf{u}·g+S_D,\hfill \\ \hfill & \frac{\partial \rho Y_i}{\partial t}+\nabla ·\left(\rho \mathbf{u}{Y}_i\right)=\nabla ·\left(\rho D_{eff}\left(\nabla Y_i\right)\right),\hfill \end{array}$$

(2)

where h is the enthalpy, $K=\frac{1}{2}{\left|\mathbf{u}\right|}^2$ is the kinetic energy of the system, $Y_i$ is the mass fraction of the ith species from the set of gas species indices given by $N=\{1,2\}$, and $S_D$ accounts for the enthalpy transport due to diffusive mass transport and the associated correction of the heat conduction:

$$S_D=\sum _{i\in N}\nabla ·\left(\rho \left(D_{eff}-{\alpha }_{eff}\right)c_{p,i}T\left(\nabla Y_i\right)\right).$$

The effective dynamic viscosity ${\mu }_{eff}$ is the sum of the molecular and turbulent or subgrid-scale viscosity regarding URANS and LES, respectively. h is the sum of the internal energy per unit mass e and the kinematic pressure $h=e+\frac{p}{\rho }$. The effective molecular diffusivity $D_{eff}$ is the sum of the molecular and turbulent or subgrid-scale molecular diffusivity: $D_{eff}=D+D_t$. The molecular diffusivity is assumed to be constant. According to the gradient flux approach, the thermal diffusivity results from ${\alpha }_{eff}=\frac{\mu }{\rho ·Pr}+\frac{{\nu }_t}{P{r}_t}$ and the molecular diffusivity results from $D_{eff}=D+\frac{{\nu }_t}{S{c}_t}$ with the turbulent Schmidt number $S{c}_t=0.85$. Mixture properties ${\phi }_m$ are computed from the individual specie properties ${\phi }_i$ and species mass fractions $Y_i$: ${\phi }_m={\sum }_{i\in N}{Y}_i{\phi }_i$. A comprehensive model description can be found in [24].

2.2. Case Definition

For the present work, different case definitions form the basis of the investigations and the evaluation of the results. For validation of the results of the mesh convergence study, the DHC with aspect ratio 4, which is filled with air, was analysed. A schematic sketch is shown in Figure 4a and is referred to as the default case. Figure 4b,c represent a test case that can be considered as an intermediate step between the THAI containment and the DHC. Furthermore, important variables for normalization and parameter variation are introduced.

The non-slip boundary condition is imposed on the velocity at the four enclosing walls at $x=0$, $x=W$, $y=0$, and $y=H$. The cavity is subject to a temperature difference $\mathsf{\Delta }T={\overline{T}}_{left}-{\overline{T}}_{right}$. Thermal radiation is neglected. The height H of the cavity is 1 m. A 2D simulation is considered with URANS. With LES, the third spatial dimension is taken into account due to the three-dimensional character of turbulent flow. Hence, the flow field is assumed to be periodic in the z-direction. For the mass fraction, a zero-gradient at the enclosing walls is defined: ${\left.\frac{\partial Y_i}{\partial \mathbf{n}}\right|}_{\partial \mathsf{\Omega }}=0$, where $\mathbf{n}$ denotes the wall-normal unit vector. For LES, a zero gradient is defined for the turbulent viscosity and thermal diffusivity: ${\left.\frac{\partial {\nu }_t}{\partial \mathbf{n}}\right|}_{\partial \mathsf{\Omega }}=0$, ${\left.\frac{\partial {\alpha }_t}{\partial \mathbf{n}}\right|}_{\partial \mathsf{\Omega }}=0$. The Wall-Adapting Local Eddy-viscosity (WALE) model with $C_w=0.5$ is applied for modelling viscous subgrid-scale effects [25]. For URANS, wall functions are applied for the turbulent viscosity and thermal diffusivity.

The k-omega Shear Stress Transport (SST) turbulence model [26] is used with included buoyancy terms, implemented in analogy with ANSYS CFX [27].

As initial conditions, temperature $T_0=\overline{T}=298.15 \mathrm{K}$ and pressure $p_0=1\mathrm{bar}$ are applied. In the mesh convergence study, the cavity is filled with air and the top and bottom walls are adiabatic: $\frac{\partial T}{\partial y}\left(y=0\right)=0$ and $\frac{\partial T}{\partial y}\left(y=H\right)=0$. The temperature over the left and right wall is constant. Since this configuration has been investigated with Direct Numerical Simulation (DNS) [13], the results from the LES can be verified. This default case is shown in Figure 4a.

As shown in Figure 4b, the same configuration is used as the reference case for the sensitivity analysis and uncertainty quantification, with the difference that the upper third contains 40 $\mathrm{vol}\%$ of helium, which is uniformly distributed. This corresponds to the uniform reference mole fraction $X_0=0.4$. The material values for air and helium at $\overline{T}=298.15 \mathrm{K}$ and $p=1\mathrm{bar}$ listed in

Table A1. Applied material properties of air and helium.

Gas	Mole Weight M	Isobaric Heat Capacity ${\mathit{C}}_{\mathit{p}}$	Dynamic Viscosity $\mathit{\mu }$	Prandtl Number $\mathit{Pr}$
air [28]	$28.96 \mathrm{g}{\mathrm{mol}}^{-1}$	$1006.5 \mathrm{J} {\mathrm{K}}^{-1}$	$1.845\times {10}^{-5}$$\mathrm{Pa}$$\mathrm{s}$	$0.707$
helium [46]	$4.0 \mathrm{g}{\mathrm{mol}}^{-1}$	$5193 \mathrm{J} {\mathrm{K}}^{-1}$	$1.985\times {10}^{-5}$$\mathrm{Pa}$$\mathrm{s}$	$0.664$

are applied. From this, the molecular diffusion coefficient, which was assumed to be constant, was also derived using Fuller’s method [28] with $D_0=6.904\times {10}^{-5}{\mathrm{m}}^2/\mathrm{s}$. The remaining properties are derived by the ideal gas law at $\overline{T}=298.15 \mathrm{K}$ with the gas constant R: $\beta =\frac{1}{\overline{T}}$, $\rho =\frac{pM}{R\overline{T}}$. To investigate the influence of the initial and boundary conditions, different parameters are varied, which were chosen to be representative of possible uncertainties in the THAI-TH22 experiment. An illustration of the parameters is provided in Figure 4c.

Constant temperature specification or adiabatic conditions at the enclosing walls provide approximations for the real prevailing conditions in THAI containment. Insulations are still thermally conductive and heat transfer between the heating fluid and the gas is associated with a spatially variable temperature. Therefore, the variation in the thermal boundary conditions includes the wall temperature difference between the left and right walls, the wall-tangential temperature gradient at the left and right walls, and the temperatures at the top and bottom walls of the cavity. The wall temperature difference is defined by

$$\begin{array}{ccccc}\hfill {\overline{T}}_{left}& =\overline{T}+\frac{\mathsf{\Delta }T}{2},\hfill & \hfill & & \hfill \mathsf{\Delta }T=\vartheta \mathsf{\Delta }{T}_0=\vartheta 21.431 \mathrm{K},\\ \hfill {\overline{T}}_{right}& =\overline{T}-\frac{\mathsf{\Delta }T}{2},\hfill & \hfill & & \hfill ({\overline{T}}_{left}+{\overline{T}}_{right})/2=\overline{T}=298.15 \mathrm{K}=const.,\end{array}$$

(3)

where $\mathsf{\Delta }{T}_0$ is the reference temperature difference, $\vartheta$ describes the relative proportion of the reference temperature difference, and $\mathsf{\Delta }T$ describes the actual temperature difference under consideration. The average temperature of the boundary is kept constant. The wall-tangential temperature gradient is defined by

$$\begin{array}{ccccc}\hfill T\left(x=0,\psi \right)=\frac{\varphi \mathsf{\Delta }T}{H}\psi +{\overline{T}}_{left}& ,\hfill & \hfill & & \hfill \psi :=y-\frac{H}{2},\\ \hfill T\left(x=W,\psi \right)=\frac{\varphi \mathsf{\Delta }T}{H}\psi +{\overline{T}}_{right}& ,\hfill & \hfill & & \hfill \frac{\partial T}{\partial y}\left(x=0\vee x=W\right)=\frac{\varphi \mathsf{\Delta }T}{H},\end{array}$$

(4)

where $\varphi$ indicates the relative proportion of the temperature change due to the temperature gradient over the entire height to the characteristic temperature difference $\mathsf{\Delta }T$. This temperature gradient also occurs in a real application within the heating and cooling zones. The defined gradient increases linearly in the vertical direction and is identical on the left and right walls to keep the local Rayleigh number constant. Since an adiabatic wall is an idealized condition, a temperature profile is also defined for the upper and lower walls to investigate the influence of the associated convective heat transfer on the mixing process. To maintain consistency in the temperature field in the corners and edges, a parabolic profile is applied that is continuous and fully defined by the inflection point and the near wall temperature. A polynomial order of two was chosen for simplicity. The temperature boundary layer width ${\delta }_T$ at the top and bottom walls is assumed to be 20% of the total width W. The definition of the temperature field at the top and bottom walls is

$$\begin{array}{cccccc}\hfill {\delta }_T=x_l=0.2W,& \hfill & \hfill T_{top}={\overline{T}}_{left}-\theta \mathsf{\Delta }T,& \hfill & \hfill (T_{top}+T_{bottom})/2& =\overline{T}=298.15 \mathrm{K}=const.,\hfill \\ \hfill x_r=0.8W,& \hfill & \hfill T_{bottom}={\overline{T}}_{right}+\theta \mathsf{\Delta }T,& \hfill & \end{array}$$

$$\begin{array}{cc}\hfill T\left(x,y=H\right)& =\left\{\begin{array}{cc}\frac{T\left(x=0,y=H\right)-T_{top}}{{\delta }_T^2}{\left(x-x_l\right)}^2+T_{top},\hfill & \mathrm{if}x\le x_l\hfill \\ \frac{T\left(x=W,y=H\right)-T_{top}}{{\delta }_T^2}{\left(x-x_r\right)}^2+T_{top},\hfill & \mathrm{if}x\ge x_r\hfill \\ T_{top},\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \\ \hfill T\left(x,y=0\right)& =\left\{\begin{array}{cc}\frac{T\left(x=0,y=0\right)-T_{bottom}}{{\delta }_T^2}{\left(x-x_l\right)}^2+T_{bottom},\hfill & \mathrm{if}x\le x_l\hfill \\ \frac{T\left(x=W,y=0\right)-T_{bottom}}{{\delta }_T^2}{\left(x-x_r\right)}^2+T_{bottom},\hfill & \mathrm{if}x\ge x_r\hfill \\ T_{bottom},\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(5)

where $\theta$ is the relative proportion of the characteristic temperature difference $\mathsf{\Delta }T$. In addition, there is uncertainty in the actual build-up of the helium stratification after the injection process. Therefore, the initial helium stratification is changed by variation in the initial mole fraction difference $\mathsf{\Delta }X$. This is the deviation in the mole fraction at the beginning of the helium stratification. The mole fraction of helium is linearly distributed over the remaining height defined by the equations

$$\begin{array}{cccc}\hfill h:=\frac{2H}{3},& \hfill & \hfill \upsilon :=y-h=y-\frac{2H}{3},& \\ \hfill \mathsf{\Delta }X=\chi X_0,& \hfill & \hfill X\left(\upsilon \right)=\frac{2\mathsf{\Delta }X}{H-h}\upsilon +\left(X_0-\mathsf{\Delta }X\right)=\frac{2\chi X_0}{H-h}\upsilon +X_0\left(1-\chi \right),\end{array}$$

(6)

where $\chi$ is the relative proportion of the constant mole fraction $X_0$ of the reference case. Due to its significant influence with respect to the mixing process, the molecular diffusion coefficient D in Equation (2) is also taken into account with

$$\begin{array}{cccc}\hfill D_{eff}& =D+D_t,\hfill & \hfill D& =\xi D_0,\hfill \end{array}$$

(7)

where $\xi$ is the relative proportion of the reference molecular diffusivity $D_0$.

2.3. Numerical Framework and Spatial Discretization

The open-source C++ toolbox OpenFOAM v.1906 [29,30] was utilized for solving the nonlinear set of governing equations in a finite-volume framework. The calculations for the incompressible flow (Equation (1)) within the mesh convergence study were performed by using the buoyantBoussinesqPimpleFoam solver. The equations for the low Mach number flow (Equation (2)) when considering the mixing process were solved by the self-defined buoyantMixingFoam solver. The pressure–velocity coupling was addressed by using the PIMPLE algorithm. It is ensured that the normalized residuals fall below the value ${10}^{-4}$. For the mesh convergence study, the convective and diffusive fluxes were evaluated by second-order linear upwind and linear schemes, respectively. For the sensitivity and uncertainty analysis, the convective momentum flux was evaluated by second-order linear upwind. The remaining convective fluxes and diffusive fluxes were evaluated by the limited linear scheme. The convective flux of the mass fraction was discretized by the limited linear scheme that was bounded between 0 and 1. Temporal advancement was achieved by blending between the Euler and Crank–Nicolson scheme with a value of 0.9, which is a good compromise between accuracy and robustness. Here, 0 and 1 correspond to Euler and Crank–Nicolson, respectively. For the mesh convergence study, it has been ensured that the CFL number is always below the value of 0.5. For the sensitivity and uncertainty analysis, the maximum CFL number of 0.9 was chosen as a compromise between accuracy and stability. The spatial grid resolution has to be fine enough to resolve most of the turbulent fluctuations for the LES. In Figure 5, the applied mesh refinement strategy is shown. The mesh is refined linearly starting from the central planes of the cavity in the direction of the walls with a constant factor. Due to the center-point-symmetrical flow conditions, the mesh was also defined point symmetrically around the cavity center. For the correct resolution of the wall-boundary layer, the dimensionless horizontal and vertical normal wall distances $x_{\perp }^{+}=x/l_{{\tau }_x}=1$ and $y_{\perp }^{+}=y/l_{{\tau }_y}=1$ are applied. $l_{{\tau }_x}=\nu /u_{{\tau }_x}$ and $l_{{\tau }_y}=\nu /u_{{\tau }_y}$ denote the viscous length at the vertical and horizontal walls, respectively. $u_{{\tau }_x}$ and $u_{{\tau }_y}$ are the shear velocities at the vertical and horizontal walls, respectively. Within a mesh convergence study the dimensionless wall tangential cell sizes $\mathsf{\Delta }{x}_{\Vert }^{+}=\mathsf{\Delta }{x}_c/l_{{\tau }_y}$ and $\mathsf{\Delta }{y}_{\Vert }^{+}=\mathsf{\Delta }{y}_c/l_{{\tau }_x}$ at the central planes of the cavity are stepwise refined. A sufficient length in the periodic direction ${\phi }_{WD}=W/D=1$ is applied to ensure that turbulence fluctuations are uncorrelated at a separation of one half-period [13]. The mesh in the periodic direction is uniformly distributed with $\mathsf{\Delta }{z}^{+}=\mathsf{\Delta }z/l_{{\tau }_x}=20$. $l_{{\tau }_x}$ is applied for the definition of $\mathsf{\Delta }{z}^{+}$ because $l_{{\tau }_x}$ takes on smaller values than $l_{{\tau }_y}$ and therefore provides a more strict condition. The mesh was built with the parameters listed in

Table A2. Physical and spatial discretization parameters.

$\mathit{Ra}$	${\mathit{N}}_{\mathit{x}}$	${\mathit{N}}_{\mathit{y}}$	${\mathit{N}}_{\mathit{z}}$	$\frac{\mathbf{\Delta }{\mathit{x}}_{\mathit{c}}}{\mathbf{\Delta }{\mathit{x}}^{+}}$	$\frac{\mathbf{\Delta }{\mathit{y}}_{\mathit{c}}}{\mathbf{\Delta }{\mathit{y}}^{+}}$	${\mathit{\phi }}_{\mathbf{WD}}$
$2\mathrm{e}9$	44	220	34	$24.439$	$9.209$	$1.0$

.

3. Results and Discussion

Before presenting the results of the sensitivity and uncertainty analysis in Section 3.3 and Section 3.4, respectively, the mesh convergence study is discussed in Section 3.1 and the quantities of interest are introduced in Section 3.2.

3.1. Mesh Convergence Study

First, a mesh convergence study for LES at $Ra=2\times {10}^9$ was performed to determine a mesh that satisfactorily balances accuracy and computing resources. The underlying mesh refinement strategy is shown in Figure 5c. The mesh is linearly refined towards the wall with a constant factor. The normal wall distance of the first cells $\mathsf{\Delta }x$ and $\mathsf{\Delta }y$, which correspond to the dimensionless wall distances $x_{\perp }^{+}=1$ and $y_{\perp }^{+}=1$, are kept constant. This ensures correct resolution of the wall-boundary layer. The wall tangential cell dimensions $\mathsf{\Delta }{x}_c$ and $\mathsf{\Delta }{y}_c$ at the central planes of the cavity are gradually refined based on values for the dimensionless tangential cell size $\mathsf{\Delta }{x}_{\Vert }^{+}$ and $\mathsf{\Delta }{y}_{\Vert }^{+}$ from best practice guidelines that can be found in the literature [31]. The recommendations are ${\mathsf{\Delta }}_{\perp }^{+}\approx 2$, ${\mathsf{\Delta }}_{\Vert }^{+}\approx 50-100$, and ${\mathsf{\Delta }}_{{\Vert }^{*}}^{+}\approx 15-20$ for the wall-normal, streamwise, and spanwise grid resolutions. Through the stepwise refinements of $\mathsf{\Delta }{x}_{\Vert }^{+}$ and $\mathsf{\Delta }{y}_{\Vert }^{+}$ through $\mathsf{\Delta }{x}_c$ and $\mathsf{\Delta }{y}_c$, the required number of cells in the boundary layer and the required wall tangential resolution can be determined, which are crucial for correct calculation. It is possible to vary the two parameters independently of each other in order to determine the appropriate wall tangential and wall normal resolution. To keep the computational effort reasonable, the parameters were changed simultaneously, assuming that the required dimensionless mesh resolution is identical on the enclosing walls. An appropriate mesh for the large parametric analysis balancing accuracy and computational effort was finally determined by using the local Nusselt number as a convergence criterion: $N{u}_{local}=-\left(H{\left.(\partial T/\partial x)\right|}_{wall}\right)/\mathsf{\Delta }T$. For this purpose, the time-averaged Nusselt number profile and the local standard deviation were investigated. The results were obtained by averaging the profiles on the left and right walls and by averaging in the z-direction for the 3D LES simulation. For time averaging, the time $t/\left(\left(H^2/\alpha \right)R{a}^{-1/2}\right)=800$ [13] was used. Afterwards, the LES results were compared with the DNS results using the mean absolute error $\left(\mathrm{MAE}\right)$ [32,33]. The DNS was defined as a reference. avg and sdev denote the average value and standard deviation, respectively. ${\langle ·\rangle }_{z,t}$ denotes the temporal and spatial average in the periodic direction. ${\left(·\right)}^{\prime }$ is the fluctuating component according to $\phi ={\langle \phi \rangle }_t+{\phi }^{\prime }$ with a flow property $\phi$ and its temporal mean component ${\langle \phi \rangle }_t$. From this, the evaluation is described by the following expressions:

$$\begin{array}{cc}\hfill \underset{\mathrm{LES},\mathrm{DNS}}{\mathrm{MAE}}\left(\mathrm{avg}\left(Nu\right)\right)& =\frac{1}{H}{\int }_{y=0}^H\left|{\langle N{u}_{LES}\rangle }_{z,t}-{\langle N{u}_{DNS}\rangle }_{z,t}\right|\mathrm{d}y,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \underset{\mathrm{LES},\mathrm{DNS}}{\mathrm{MAE}}\left(\mathrm{sdev}\left(Nu\right)\right)& =\frac{1}{H}{\int }_{y=0}^H\left|{\left({\langle {\left(N{u}_{LES}^{\prime }\right)}^2\rangle }_{z,t}\right)}^{\frac{1}{2}}-{\left({\langle {\left(N{u}_{DNS}^{\prime }\right)}^2\rangle }_{z,t}\right)}^{\frac{1}{2}}\right|\mathrm{d}y.\hfill \end{array}$$

(9)

Figure 6a–d clearly shows that, with increasing refinement, the LES results are in very good agreement with the DNS results. Figure 6a,b shows the normalized MAE of the LES from the DNS results regarding the average and standard deviation of $Nu$, which gradually decreases when the mesh is refined. It can be seen in Figure 6a,b that, for $\mathsf{\Delta }{x}_{\Vert }^{+},\mathsf{\Delta }{y}_{\Vert }^{+}=10$, the difference is very close to zero. The time-averaged local $Nu$ number profile in Figure 6c already shows good results with a coarse resolution of the domain. There are minor differences near the top and bottom of the cavity. A detail plot within Figure 6c shows the gradual convergence of the LES to the DNS profile. However, for a coarse grid, the transition point, which corresponds to the point with maximum fluctuations and thus the largest standard deviation, is not predicted properly, as can be seen in Figure 6d. In addition, no recirculation zone is formed. The formation of the recirculation zone, as previously described by Trias et al. [13], only takes place in the case of fine grids.

Due to the fact that the recirculation zone does not develop and the transition occurs further downstream for coarse grids, the fluid is diverted by the horizontal walls instead and causes fluctuations when it hits the vertical walls.

To define a mesh that is accurate enough to reflect the global characteristics compared to the DNS, an upper limit for the relative deviation between LES and DNS was defined: $\underset{\mathrm{LES},\mathrm{DNS}}{\mathrm{MAE}}\left(\mathrm{avg}\left(Nu\right)\right)/{\langle \mathrm{avg}\left(N{u}_{DNS}\right)\rangle }_y\le {ϵ}_{MAE}\approx {10}^{-2}$ with ${\langle \mathrm{avg}\left(N{u}_{DNS}\right)\rangle }_y=66.63$ [13]. ${\langle ·\rangle }_y$ denotes the spatial average in the y-direction. This condition is fulfilled from a grid resolution of $\mathsf{\Delta }{x}_{\Vert }^{+},\mathsf{\Delta }{y}_{\Vert }^{+}=30$. Therefore, this mesh resolution was selected for the following investigations. The same spatial resolution was chosen for the LES and URANS to ensure that the spatial discretization error is comparable and provides a direct comparison of the methods independent of the discretization techniques. When considering the LES results, the structure and spatial resolution of the boundary layer were investigated, as can be seen in Figure 7. The thickness of the velocity boundary layer ${\delta }_u$ was determined by the occurring inflection point $\left(\mathrm{ip}\right)$ in the velocity profile, which is characteristic for natural convection. The temperature boundary layer ${\delta }_T$ was defined by a gradient criterion $\left(\mathrm{grad}\right)$: $\partial T/\partial x\ge {10}^{-2}{\left.\partial T/\partial x\right|}_{wall}$. These criteria made it possible to investigate the velocity and temperature boundary layer thickness over the entire height in a consistent manner and to determine the required number of cells in the boundary layer. The maximum root-mean-square (rms) of the velocity, which takes place after reaching the maximum velocity, is approximately in the area of the maximum velocity gradient, characterized by the inflection point. Therefore, the definition by the inflection point and the rms criterion are in good agreement at the middle height plane, as shown in Figure 7a. The situation is analogous for the temperature boundary layer. To achieve comparable results to the rms criterion, the gradient criterion can be adjusted accordingly. From the results the mean cell number was determined with the following expression: ${\overline{n}}_{\delta }=\frac{1}{H}\sum \mathsf{\Delta }{y}_i{n}_{\delta }$, where $\mathsf{\Delta }{y}_i$ is the vertical cell length of the corresponding cell, $n_{\delta }$ is the respective number of cells in the boundary layer, and ${\overline{n}}_{\delta }$ is the mean number of cells in the boundary layer over the entire height H. With the previously defined criteria, this results in a velocity and temperature boundary layer resolution of ${\overline{n}}_{{\delta }_{u/ip}}=7.87$ and ${\overline{n}}_{{\delta }_{T/grad}}=10.10$ respectively. The determination of the required wall normal resolution with ${\overline{n}}_{{\delta }_{u/ip}}$ and ${\overline{n}}_{{\delta }_{T/grad}}$ and the required wall tangential resolution $\mathsf{\Delta }{x}_{\Vert }^{+}$ and $\mathsf{\Delta }{y}_{\Vert }^{+}$ with the present mesh refinement strategy enables the transfer of results to problems with similar physics such as the THAI containment.

3.2. Quantities of Interest

Before the assessment of uncertainties and analyzing sensitivities, the result quantities, which adequately describe the underlying physics of the mixing process in the cavity, have to be defined. In the context of uncertainty quantification, these are referred to as Quantities of Interest (QoI). For this purpose, integral scalar quantities provide a plain description of the transient profiles of the mixing process. In this way, summarizing statements can be made about the complete transient. Since the flow is driven by buoyancy effects, the spatial averaged Nusselt number $Nu$ over the respective walls is of paramount importance. In this way, the convective heat transfer in the cavity is evaluated and the transients during the mixing process can be examined more closely. It is determined by the expression

$$Nu=\frac{1}{A_{wall}}{\int }_{A_{wall}}\frac{{\left.\frac{\partial T}{\partial \mathbf{n}}\right|}_{wall}H}{\mathsf{\Delta }T}\mathrm{d}A,$$

(10)

where $\mathbf{n}$ denotes the wall-normal unit vector. The effects on the convection mechanisms are described by global kinetic energy. $E_k$ is the quotient of the global kinetic energy by a reference kinetic energy $\frac{{\alpha }^2}{H^2}Ra$, which contains the material properties of air:

$$E_k={\left(\frac{{\alpha }^2}{H^2}Ra\right)}^{-1}·\frac{1}{M}{\int }_M\frac{1}{2}{u}^2\mathrm{d}m.$$

(11)

m denotes the mass and M denotes the total mass of the fluid domain. Due to the thermal boundary conditions, natural convection develops within the cavity and the helium stratification is eroded. This mixing process continues until a homogeneous mixture finally forms. The achievement of this state can be quantified by the mixture uniformity ${\sigma }_X$, which is the volume-weighted standard deviation of the mole fraction X from the homogeneous equilibrium state mole fraction $\overline{X}$ over the whole fluid domain. The initial mixture uniformity ${\sigma }_{X_0}$ is the highest occurring standard deviation during the mixing process. When varying the initial mole fraction difference $\mathsf{\Delta }X$, the initial mixture uniformity ${\sigma }_{X_0}$ changes. Therefore, for normalization, the initial mixture uniformity ${\sigma }_{X_{0/ref}}$ of the reference case, which corresponds to the initial conditions of a stratification with constant mole fraction, is used. Following the definition of Danckwerts [34], this gives the expression for the segregation intensity, described by the Equation (12):

$$I=\frac{{\sigma }_X^2}{{\sigma }_{X_{0/ref}}^2}=\frac{\frac{1}{V}{\int }_V{\left(X-\overline{X}\right)}^2\mathrm{d}V}{\frac{1}{V}{\int }_V{\left(X_{0/ref}-\overline{X}\right)}^2\mathrm{d}V}$$

(12)

From the segregation intensity I, the definition of the mixing intensity $\mathcal{M}=1-\sqrt{I}$ [35] can also be derived for the interpretation of the results. In the reference case, $I=1$ describes a completely inhomogeneous mixture and $I=0$ characterizes a completely homogeneous mixture and vice versa for the mixing intensity $\mathcal{M}$. Since the initial state changes, when the mole fraction difference $\mathsf{\Delta }X$ is varied, the initial segregation intensity I is also different and the mixture transients are consistently captured with all occuring changes. Finally, a criterion for the mixing time can be derived from I. A completely homogeneous mixture is characterized by $I=0$, and therefore by definition of an upper bound $I\le {ϵ}_I={10}^{-3}$, which I has to fall below, the achievement of this state can be quantified. When considering mass transfer processes, the Fourier number $Fo=D_{ref}t/H^2$ enables a dimensionless description of time, where $D_{ref}=D_0$ is the diffusion coefficient of the reference case. Hence, the time, when the homogeneous state is reached, can be described by $F{o}_{ϵ}=D_{ref}{t}_{ϵ}/H^2$. The time-dependent quantities described above are combined to the result variables vector $\mathbf{R}$, where the bold notation indicates a vector or matrix.

$$\mathbf{R}=\left(\begin{array}{cccccc}N{u}_{left}& N{u}_{right}& N{u}_{bottom}& N{u}_{top}& E_k& I\end{array}\right)$$

(13)

Together with $F{o}_{ϵ}$, additional scalar quantities that describe the mixing process were derived. For this purpose, the integral mean value $\tilde{\mathbf{R}}$ of the described quantities $\mathbf{R}$ (Equation (13)) is formed over the respective mixing time in Equation (14).

$$\tilde{\mathbf{R}}=\frac{1}{F{o}_{ϵ}}{\int }_0^{F{o}_{ϵ}}\mathbf{R}\left(Fo\right)\mathrm{d}Fo$$

(14)

Next to the actual integral mean value over the mixing process, this also provides a measure of the shape of the mixing transient because transients, which can be mapped on each other by linear stretching or compression of the time coordinate, have an identical integral mean value. Therefore, occurring deviations also indicate a change in the shape of the mixing transient. Because different transients can have the same integral mean value $\tilde{\mathbf{R}}$, the mean absolute deviation $\widehat{\mathbf{R}}$ from the reference case over the reference mixing time $F{o}_{ϵ/ref}$ was therefore defined in Equation (15).

$$\widehat{\mathbf{R}}=\frac{1}{F{o}_{ϵ/ref}}{\int }_0^{F{o}_{ϵ/ref}}\left|\mathbf{R}\left(Fo\right)-{\mathbf{R}}_{ref}\left(Fo\right)\right|\mathrm{d}Fo$$

(15)

All differences in the mixing transient are captured over the mixing time of the reference case and provide a good measure for the variability. At the same time, it can be determined how big the difference between the adiabatic configuration and a configuration with temperature specification is. Finally, all variables under investigation are summarized in a matrix in Equation (16).

$$\underline{\mathbf{R}}=\left(\begin{array}{cc}\tilde{\mathbf{R}}& F{o}_{ϵ}\\ \widehat{\mathbf{R}}& 0\end{array}\right)=\left(\begin{array}{ccccccc}{\widetilde{Nu}}_{left}& {\widetilde{Nu}}_{right}& {\widetilde{Nu}}_{bottom}& {\widetilde{Nu}}_{top}& {\tilde{E}}_k& \tilde{I}& F{o}_{ϵ}\\ {\widehat{Nu}}_{left}& {\widehat{Nu}}_{right}& {\widehat{Nu}}_{bottom}& {\widehat{Nu}}_{top}& {\widehat{E}}_k& \widehat{I}& 0\end{array}\right)$$

(16)

The results for the reference case are different for LES and URANS. The reference values ${\underline{\mathbf{R}}}_{ref/LES}$ are applied to the LES results and the reference values ${\underline{\mathbf{R}}}_{ref/URANS}$ are applied to the URANS results. As listed in Equation (17), the mean integral values of the reference case ${\tilde{\mathbf{R}}}_{ref}$ are applied for the normalization of the later results, since the mean absolute deviation ${\widehat{\mathbf{R}}}_{ref}$ from the reference case itself is zero.

$${\underline{\mathbf{R}}}_{ref}=\left(\begin{array}{cc}{\tilde{\mathbf{R}}}_{ref}& F{o}_{{ϵ}_{ref}}\\ {\tilde{\mathbf{R}}}_{ref}& 0\end{array}\right)$$

(17)

3.3. Sensitivity Analysis

For the investigation of uncertainties, identification of the most influential parameters is of great importance, since uncertainties occurring in these parameters mainly dominate the uncertainty in the results. Variance-based methods also enable quantification of the proportion of the result uncertainty caused by individual parameters. Therefore, uncertainty analysis is closely related to sensitivity analysis. For assessment of the sensitivities, the Morris method is applied [36]. $\mathbf{q}$ denotes the parameter vector with the elements $q_i$. The defined parameter space of the computational model is explored by r-trajectories, for which the initial points are gained by random sampling, and then one parameter after the other is changed step by step. The respective trajectory is described by the index j. Then, for each trajectory, the elementary effects ${\mathbf{d}}_{q_i}$ are calculated, which represent the gradient when varying the ith parameter using the forward difference. From this, the mean and modified mean of the elementary effects ${\mu }_{q_i}$ and ${\mu }_{q_i}^{*}$ as well as the standard deviation of the elementary effects ${\sigma }_{q_i}$ are determined. The modified mean indicates the overall effect and the standard deviation indicates nonlinear or interaction effects between the input parameters. For evaluation of the sensitivity, both parameters are important because a low value of ${\mu }_{q_i}^{*}$ combined with a considerable value of ${\sigma }_{q_i}$ means that there is a great impact in certain regions of the parameter space and that the corresponding parameter should be taken into account for the subsequent analysis. The expressions for the calculation [36] are

$$\begin{array}{cccccc}\hfill {\mathbf{d}}_{q_i}^{\left(j\right)}=\frac{\underline{\mathbf{R}}\left({\mathbf{q}}^{\left(j\right)}+{\mathsf{\Delta }}_{\mathsf{\Theta }}·{\mathbf{e}}_{q_i}\right)-\underline{\mathbf{R}}\left({\mathbf{q}}^{\left(j\right)}\right)}{\mathsf{\Delta }},& \hfill & \hfill {\mu }_{q_i}=\frac{1}{r}\sum _{j=1}^r{\mathbf{d}}_{q_i}^{\left(j\right)},& \hfill & \hfill \mathsf{\Delta }=\frac{p}{2\left(p-1\right)},& \\ \hfill {\sigma }_{q_i}=\sqrt{\frac{1}{r-1}\sum _{j=1}^r{\left({\mathbf{d}}_{q_i}^{\left(j\right)}-{\mu }_{q_i}\right)}^2},& \hfill & \hfill {\mu }_{q_i}^{*}=\frac{1}{r}\sum _{j=1}^r\left|{\mathbf{d}}_{q_i}^{\left(j\right)}\right|,& \hfill & & \hfill \end{array}$$

(18)

where ${\mathbf{e}}_{q_i}$ is the ith unit vector and the parameter step is ${\mathsf{\Delta }}_{\mathsf{\Theta }}=\mathsf{\Delta }\mathsf{\Theta }=\frac{p}{2(p-1)}\mathsf{\Theta }$ with the input space uniformly partitioned into the even number of levels p. $\mathsf{\Theta }$ denotes the parameter interval width vector and $\mathsf{\Delta }$ describes the scalar parameter step. Global sensitivity analysis requires a large number of trajectories r per parameter. However, a large number of parameters in the mixing process are involved, which might be subject to uncertainties. Therefore, the total number of parameters was divided into three parameter vectors to maintain clarity for the further steps, such that the vector of input parameters is given by $\mathbf{q}=\left(\begin{array}{ccc}{\mathbf{q}}_{\mathrm{I}}& {\mathbf{q}}_{\mathrm{II}}& {\mathbf{q}}_{\mathrm{III}}\end{array}\right)$ with:

$$\begin{array}{cc}\hfill {\mathbf{q}}_{\mathrm{I}}& =\left(\begin{array}{ccc}\vartheta & \chi & \varphi \end{array}\right),{\mathbf{q}}_{\mathrm{I}\mathrm{I}}=\left(\begin{array}{cc}\xi & \theta \end{array}\right),\hfill \\ \hfill {\mathbf{q}}_{\mathrm{I}\mathrm{I}\mathrm{I}}& =\left(\begin{array}{ccccccccc}{T}_0& c_{p/Air}& {\mu }_{Air}& P{r}_{Air}& c_{p/He}& {\mu }_{He}& P{r}_{He}& S{c}_t& P{r}_t\end{array}\right).\hfill \end{array}$$

(19)

First, the most important parameters are extracted by preselection with a one-at-a-time (OAT) approach. Starting from an initial point ${\mathbf{q}}_{oat}^0$ in parameter space, the parameters were changed along each orthogonal dimension with two trajectories $r=2$, respectively. The initial point was defined as follows:

$$\begin{array}{cc}\hfill {\mathbf{q}}_{\mathrm{I}}^0& =\left(\begin{array}{ccc}1.0& 0.0& 0.1\end{array}\right),{\mathbf{q}}_{\mathrm{I}\mathrm{I}}^0=\left(\begin{array}{cc}1.0& 0.2\end{array}\right),\hfill \\ \hfill {\mathbf{q}}_{\mathrm{I}\mathrm{I}\mathrm{I}}^0& =\left(\begin{array}{ccccccccc}{T}_0^0& c_{p/Air}^0& {\mu }_{Air}^0& P{r}_{Air}^0& c_{p/He}^0& {\mu }_{He}^0& P{r}_{He}^0& S{c}_t^0& P{r}_t^0\end{array}\right).\hfill \end{array}$$

(20)

${\mathbf{q}}_{\mathrm{I}\mathrm{I}\mathrm{I}}^0$ takes the initial values from Section 2.2 and according to Table A1. The most influential parameters were determined by assuming an identical parameter interval, which has the width of 20 percent with respect to the corresponding nominal value in each parameter. The associated values for $\mathsf{\Theta }$ and $\mathsf{\Delta }$ result according to

Table 1. Parameter variation values for the one-at-a-time (OAT) approach.

$q_{oa{t}_i}$	$\vartheta$	$\chi$	$\varphi$	$\xi$
$[q_{oa{t}_i}^a,q_{oa{t}_i}^b]$	$[0.9,1.1]$	$[0.0,0.2]$	$[0.0,0.2]$	$[0.9,1.1]$
${\mathsf{\Theta }}_i$	$0.2$	$0.2$	$0.2$	$0.2$
$\mathsf{\Delta }$	$\pm 0.5$	$+0.5/+1.0$	$\pm 0.5$	$\pm 0.5$
$q_{oa{t}_i}$	$\theta$	$T_0$	$c_{p/Air}$	${\mu }_{Air}$
$[q_{oa{t}_i}^a,q_{oa{t}_i}^b]$	$[0.1,0.3]$	$T_0^0+[-0.1,0.1]\mathsf{\Delta }{T}_0$	$[0.9,1.1]c_{p/Air}$	$[0.9,1.1]{\mu }_{Air}$
${\mathsf{\Theta }}_i$	$0.2$	$0.2\mathsf{\Delta }{T}_0$	$0.2c_{p/Air}$	$0.2{\mu }_{Air}$
$\mathsf{\Delta }$	$\pm 0.5$	$\pm 0.5$	$\pm 0.5$	$\pm 0.5$
$q_{oa{t}_i}$	$P{r}_{Air}$	$c_{p/He}$	${\mu }_{He}$	$P{r}_{He}$
$[q_{oa{t}_i}^a,q_{oa{t}_i}^b]$	$[0.9,1.1]P{r}_{Air}$	$[0.9,1.1]c_{p/He}$	$[0.9,1.1]{\mu }_{He}$	$[0.9,1.1]P{r}_{He}$
${\mathsf{\Theta }}_i$	$0.2P{r}_{Air}$	$0.2c_{p/He}$	$0.2{\mu }_{He}$	$0.2P{r}_{He}$
$\mathsf{\Delta }$	$\pm 0.5$	$\pm 0.5$	$\pm 0.5$	$\pm 0.5$
$q_{oa{t}_i}$	$S{c}_t$	$P{r}_t$
$[q_{oa{t}_i}^a,q_{oa{t}_i}^b]$	$[0.9,1.1]S{c}_t$	$[0.9,1.1]P{r}_t$
${\mathsf{\Theta }}_i$	$0.2S{c}_t$	$0.2P{r}_t$
$\mathsf{\Delta }$	$\pm 0.5$	$\pm 0.5$

. Through this definition, the parameters are varied symmetrically around the defined initial point with the assumed standard deviation ${\sigma }_{\mathcal{Q}}$ in the parameters according to Section 3.4. The standard deviation is representative for the total dispersion around the mean, and in this way, an estimation for the mean impact of single parameters is achieved.

From the results, the modified mean ${\mu }_{q_i}^{*}$ can be obtained. The standard deviation ${\sigma }_{q_i}$ was not taken into account here, as the number of trajectories is small. Here, the main criterion for the evaluation of the sensitivity was the mixing time measured by the dimensionless Fourier number $F{o}_{ϵ}=D_{ref}{t}_{ϵ}/H^2$. It can be seen in Figure 8 that the material values such as the heat capacity $c_{p/Air}$ and $c_{p/He}$, the dynamic viscosities ${\mu }_{Air}$ and ${\mu }_{He}$, as well as $P{r}_{He}$ have no great influence on the mixing time for $Ra=2\times {10}^9$. Additionally, the model coefficients of the gradient flux approach $S{c}_t$ and $P{r}_t$ do not show major impacts on the mixing time. On the other hand, $P{r}_{Air}$ causes a medium impact, since the main proportion of the mixture is air. However, it is assumed that there is a small uncertainty in the prediction of material values. The thermal boundary conditions via $\vartheta$, $\varphi$, and $\theta$ show great influence on the mixing time, as the mixing process is primarily based on buoyancy effects. The diffusion coefficient via $\xi$ also leads to a large deviation. The change in the helium stratification through $\chi$ also causes relevant differences in the mixing time.

On the basis of these findings, the five most important parameters were selected. These are investigated next by Global Sensitivity Analysis (GSA) using the Morris method according to Equation (18). Consequently, the parameter vector, which contains the prioritized parameters, gives the following expression:

$${\mathbf{q}}_{morris}=\left(\begin{array}{cc}{\mathbf{q}}_{\mathrm{I}}& {\mathbf{q}}_{\mathrm{I}\mathrm{I}}\end{array}\right)=\left(\begin{array}{ccccc}\vartheta & \chi & \varphi & \xi & \theta \end{array}\right)$$

(21)

The input space is uniformly partitioned into four levels $p=4$. For the parameter intervals in

Table 2. Parameter variation values for the Morris method.

$q_{morri{s}_i}$	$\vartheta$	$\chi$	$\varphi$	$\xi$	$\theta$
$[q_{morri{s}_i}^a,q_{morri{s}_i}^b]$	$[0.8,1.2]$	$[0.0,0.4]$	$[0.0,0.3]$	$[0.8,1.2]$	$[0.0,0.4]$
${\mathsf{\Theta }}_i$	$0.4$	$0.4$	$0.3$	$0.4$	$0.4$

scaled to $[0,1]$, this results in the scalar parameter step $\mathsf{\Delta }=2/3$. A random initial point ${\mathbf{q}}_{morris}^0$ is chosen from the set $\{0,1/(p-1),2/(p-1),\cdots ,1\}$ for each parameter. Starting from this point, the parameters are changed in an arbitrary order. The assessment of the sensitivity was conducted by $r=15$ trajectories. In Table 2, the intervals under investigation are listed. The parameter intervals were defined according to the probability density functions according to Equation (31) in Section 3.4. With ${\mu }_{\mathcal{Q}}$ and ${\sigma }_{\mathcal{Q}}$ denoting the expectation and standard deviation of the random input variables $\mathcal{Q}$, the interval was defined by $[{\mu }_{\mathcal{Q}}-2{\sigma }_{\mathcal{Q}},{\mu }_{\mathcal{Q}}+2{\sigma }_{\mathcal{Q}}]$ as it contains the majority of all possible parameter values. The aim of this study was to capture local interaction effects and possible nonlinearities using a reasonable number of trajectories. With a constant number of levels p, a wide interval width would lead to larger parameter steps, which cause less accurate prediction of local effects. For this reason, the parameter intervals were chosen in this way.

All previously described result variables according to Section 3.2 are examined. The integral mean Nusselt number $\widetilde{Nu}$ measures the convective heat transfer between the fluid and the enclosing walls during the mixing process, which significantly influences the mixing process. Heat is supplied on the left wall and dissipated on the right wall. Cold gas coming from the right wall is heated by the lower wall and hot gas coming from the left wall is cooled by the upper wall. ${\widetilde{Nu}}_{left}$ and ${\widetilde{Nu}}_{right}$, which are shown in Figure 9a,b, exhibit the biggest change through the temperature specification of the associated walls via $\vartheta$ because this simultaneously causes a change in the wall adjacent temperature gradient. The remaining parameters show medium impact, which is comparable for these parameters, when the values of the respective opposite walls are considered averaged with each other. LES and URANS show almost linear behavior regarding the integral mean Nusselt number in Equation (14), as can be seen in Figure 9e,f. Both approaches show good agreement. Additionally, for ${\widetilde{Nu}}_{bottom}$ and ${\widetilde{Nu}}_{top}$, which are shown in Figure 9c,d, the corresponding temperature specification via $\theta$ causes the greatest influence. The influence of $\chi$ is greater at the top wall than at the bottom wall because the helium is initially placed at the top and thus causes an influence during the mixing process. $\xi$ leads to the smallest impact for this case. In comparison to the left and right walls, ${\sigma }_{q_i}$ has larger values for the bottom and top walls. There is good correspondence between LES and URANS for both ${\mu }_{q_i}^{*}$ and ${\sigma }_{q_i}$ at the bottom and top walls.

$\widehat{Nu}$ is a quantity for measuring any occurring temporal local deviation in the convective heat transfer from the reference case and thus provides information about the influence of the parameters over the whole mixing transient. The results depicted in Figure 10 show a similar trend to that for $\widetilde{Nu}$. Therefore, the conclusions from $\widetilde{Nu}$ can be transferred to the results for $\widehat{Nu}$. An important difference exists for the measure ${\sigma }_{q_i}$, which captures nonlinear effects. The target function $\widehat{\mathbf{R}}$ according to Equation (15) evaluates the enclosed area between the mixing transients. Hence, nonlinear effects arise in the case of a strongly nonlinear profile of the transient if it is shifted through a shorter mixing time. In addition, longer mixing transients are cut off when the mixing time of the reference case is reached, which also results in nonlinearities. This fact explains the higher values in ${\sigma }_{q_i}$. This approach was chosen because the fixed integration interval enables consistent normalization. An integration interval with the mixing time of the respective case would cause identical nonlinear effects with no consistent normalization opportunity. Defining a long integration interval that contains all mixing times has the disadvantage that steady-state values are taken into account, which are irrelevant for the mixing process and therefore falsify the results. The evaluation according to Equation (15) has minor drawbacks and therefore offers a good basis for the evaluation.

The thermal boundary conditions via $\vartheta$, $\theta$, and $\varphi$ have great impact on the integral mean global kinetic energy during the mixing process $\widetilde{E_k}$, as can be seen in Figure 11a. The remaining parameters show minor impact. $\widetilde{E_k}$ also behaves almost linearly when the parameters are changed, as can be seen in Figure 11c. In contrast to ${\tilde{E}}_k$, the other parameters also show noticeable influence for the deviation in kinetic energy from the reference case $\widehat{E_k}$, shown in Figure 11b. This is due to the fact that, although the integral mean kinetic energy $\widetilde{E_k}$ remains the same by changing these parameters, the mixing time changes relative to the reference case and this leads to the measured deviation. ${\sigma }_{q_i}$ in Figure 11d correlates to ${\mu }_{q_i}^{*}$, since the target function itself is nonlinear, and consequently, larger change yields greater associated nonlinearity.

The duration of the entire mixing process is measured by the mixing Fourier number $F{o}_{ϵ}$. As can be seen in Figure 12a, the associated modified mean ${\mu }_{q_i}^{*}$ is greatest for $\vartheta$ and $\xi$. ${\mu }_{q_i}^{*}$ has a comparable order of magnitude for the other parameters, but $\varphi$ and $\chi$ have slightly higher ${\mu }_{q_i}^{*}$ than $\theta$. In Figure 12d, the nonlinearities or interaction effects are largest for $\vartheta$, $\varphi$, and $\xi$. When $\varphi$ is changed, larger interaction behavior or nonlinearities occur with URANS compared to LES, as can be seen from ${\sigma }_{q_i}$. In summary, LES and URANS show comparable sensitivity with regard to the mixing time.

The change in shape or rather the integral mean segregation intensity $\tilde{I}$ according to ${\mu }_{q_i}^{*}$ in Figure 12b is largely determined by $\chi$ because the initial distribution of the helium mole fraction varies. This is reflected in the mixing behavior and the shape of the mixing transient. The remaining parameters show medium impact, and the shape of the mixing transient is mostly retained. The conclusions from $F{o}_{ϵ}$ can largely be applied to the mean absolute deviation in the segregation intensity $\widehat{I}$. $\vartheta$ and $\xi$ have the greatest influence, as they also have the greatest influence on $F{o}_{ϵ}$. Hence, the change in these parameters shifts the mixing transients regarding time and consequently increases $\widehat{I}$. The larger values for ${\sigma }_{q_i}$ in Figure 12f are due to the nonlinearity effects of the target function $\widehat{\mathbf{R}}$, which has already been discussed.

Since a large number of parameters enormously increases the required computing resources for the investigation of uncertainties, three of the remaining parameters will be examined further. The mixing Fourier number $F{o}_{ϵ}$ was selected as the main criterion to pick the most relevant parameters. The characteristic temperature difference via $\vartheta$ has the greatest importance for the buoyancy-driven mixing process. $\chi$ and $\varphi$ have slightly higher impact than $\theta$. For this reason, $\vartheta$, $\chi$, and $\varphi$ were chosen for the uncertainty quantification. The diffusion coefficient via $\xi$ has also great importance for the mass transfer. For the application to a binary mixture, however, it is assumed that the diffusion coefficient can be described sufficiently accurate by a dynamic model and is subject to minor uncertainties. Therefore, this parameter is also neglected in the following considerations.

3.4. Uncertainty Quantification

Based on the preceding sensitivity analysis in Section 3.3, the propagation of uncertainties in one initial condition, which is represented by the change in the initial mole fraction difference of the helium stratification through $\chi$, and two thermal boundary conditions, which are represented by the change in the characteristic temperature difference through $\vartheta$ and the change in the wall tangential temperature gradient through $\varphi$, were investigated.

The open-source software Dakota 6.10 [37] was used as the uncertainty quantification framework. For the evaluation with LES and URANS, non-intrusive Polynomial Chaos Expansions (PCE) were applied because of the high convergence rate of the stochastic results with an increasing number of simulation runs. Therefore, very accurate results can potentially be obtained even with a small number of calculations. The random input variables $\mathcal{Q}:\mathsf{\Omega }\to \mathsf{{\rm Y}}\subset {\mathbb{R}}^n$ are functions that map events $\omega \in \mathsf{\Omega }$ from the sample space $\mathsf{\Omega }$ to realizations $\mathbf{q}\in \mathsf{{\rm Y}}$. PCE is a spectral method in which random response functions $\underline{\mathcal{R}}\left(\omega \right)$ are described by suitable multidimensional orthogonal polynomials ${\mathsf{\Psi }}_j\left(\mathcal{Q}\right)$ as a function of the random input variables $\mathcal{Q}$. With a sequence ${\left\{{\mathcal{Q}}_i\left(\omega \right)\right\}}_{i=1}^{\infty }$ of random variables, the infinite expansion results in expression Equations (22) and (23) [38,39]:

$$\begin{array}{cc}\hfill \underline{\mathcal{R}}\left(\omega \right)& =\underline{\mathbf{R}}\left(\mathcal{Q}\right)=a_0{B}_0+\sum _{i_1=1}^{\infty }{a}_{i_1}{B}_1\left({\mathcal{Q}}_{i_1}\right)+\sum _{i_1}^{\infty }\sum _{i_2}^{i_1}{a}_{i_1{i}_2}{B}_2\left({\mathcal{Q}}_{i_1},{\mathcal{Q}}_{i_2}\right)+\dots \hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \underline{\mathcal{R}}\left(\omega \right)& =\sum _{j=0}^{\infty }{\alpha }_j{\mathrm{\Psi }}_j\left({\mathcal{Q}}_1,{\mathcal{Q}}_2,\dots \right)\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\underline{\mathcal{R}}}^P\left(\omega \right)& =\sum _{j=0}^P{\alpha }_j{\mathrm{\Psi }}_j\left({\mathcal{Q}}_1,{\mathcal{Q}}_2,\dots ,{\mathcal{Q}}_n\right)=\sum _{j=0}^P{\alpha }_j{\mathrm{\Psi }}_j\left(\mathcal{Q}\right)\hfill \end{array}$$

(24)

There is a one-to-one correspondence between the PCE coefficients $a_{i_1{i}_2\cdots i_n}$ and ${\alpha }_j$ and between the multidimensional orthogonal polynomials $B_n\left({\mathcal{Q}}_{i_1},{\mathcal{Q}}_{i_2},\cdots ,{\mathcal{Q}}_{i_n}\right)$ and ${\mathrm{\Psi }}_j\left(\mathcal{Q}\right)$. In Equation (22), each additional set of nested summations designates a collection of polynomials with increasing order. Term-based indexing in Equation (23) instead of order-based indexing simplifies the expression. Finally, a limited number of random variables n and order truncation leads to Equation (24) with P summation terms. According to [39], the orthogonal polynomials are generated numerically by using Gauss–Wigert [40], discretized Stieltjes [41], Chebyshev [41], or Gramm–Schmidt [42] approaches. The Gauss points and weights are computed by the Golub–Welsch [43] tridiagonal eigensolution. This allows us to define arbitrary probability density functions for the input variables and eliminates the need to induce additional nonlinearity through variable transformations. The PCE coefficients ${\alpha }_j$ are estimated here by using spectral projection. The orthogonality property of the polynomials helps to extract each coefficient. The following expression [38], which contains the inner product ${\langle ·,·\rangle }_{\varrho }$ on $\mathsf{{\rm Y}}$ with the weight ${\varrho }_{\mathcal{Q}}\left(\mathbf{q}\right)$, gives the coefficients by

$${\alpha }_j=\frac{1}{{\gamma }_j}\mathbb{E}\left[\underline{\mathbf{R}}\left(\mathcal{Q}\right){\mathrm{\Psi }}_j\left(\mathcal{Q}\right)\right]=\frac{1}{{\gamma }_j}{\langle \underline{\mathbf{R}},{\mathsf{\Psi }}_j\rangle }_{\varrho }=\frac{1}{{\gamma }_j}{\int }_{\mathsf{{\rm Y}}}\underline{\mathbf{R}}\left(\mathbf{q}\right){\mathrm{\Psi }}_j\left(\mathbf{q}\right){\varrho }_{\mathcal{Q}}\left(\mathbf{q}\right)\mathrm{d}\mathbf{q},$$

(25)

where ${\varrho }_{\mathcal{Q}}\left(\mathbf{q}\right)={\prod }_{i=1}^n{\varrho }_{{\mathcal{Q}}_i}\left(q_i\right)$ is the joint probability density (weight) function. The inner product ${\gamma }_j={\langle {\mathrm{\Psi }}_j,{\mathrm{\Psi }}_j\rangle }_{\varrho }$ can be computed analytically. For solving the multidimensional integral in Equation (25), the discrete projection, which is also termed pseudospectral, is applied. The multidimensional integral can be approximated by a tensor product of one-dimensional quadrature formulas. A one-dimensional quadrature operator with the level $l\in {\mathbb{N}}_{+}$, the quadrature points $q_l^k=\left[q_l^1,\dots ,q_l^{K_l}\right]$, and a function $\mathit{f}\in C^{\alpha }$ gives the following expression [38]:

$${\mathcal{U}}_l^{\left(1\right)}\circ \mathit{f}\left(q\right)=\sum _{k=1}^{K_l}\mathit{f}\left(q_l^k\right)w_l^k$$

(26)

For the multivariate case $n>1$ and a multi-index $\mathbf{l}=\left(l_1,\dots ,l_n\right)\in {\mathbb{N}}_{+}^n$ with $\left|\mathbf{l}\right|={\sum }_{i=1}^n{l}_i$, the full tensor product quadrature formula results in Equation (27) [38].

$$\begin{array}{cc}\hfill {\mathcal{Q}}_{\mathbf{l}}^{\left(n\right)}\circ \mathit{f}\left(\mathbf{q}\right)& =\left({\mathcal{U}}_{l_1}^{\left(1\right)}\otimes \cdots \otimes {\mathcal{U}}_{l_n}^{\left(1\right)}\right)\circ \mathit{f}\left(\mathbf{q}\right)\hfill \\ \hfill & =\sum _{k_1=1}^{K_{l_1}}\cdots \sum _{k_n=1}^{K_{l_n}}\mathit{f}\left(q_1^{k_1},\dots ,q_n^{k_n}\right)\left(w_{l_1}^{k_1}\otimes \cdots \otimes w_{l_n}^{k_n}\right)\hfill \end{array}$$

(27)

However, for evaluation of this full tensor product, a very large number of function evaluations is required. Therefore, multi-dimensional integration by the Smolyak sparse grid method [44] according to Equation (28) is performed, which tremendously reduces the number of quadrature points while a high accuracy is preserved. The sparse grid quadrature rule is defined by the following expression [38,39]:

$$\begin{array}{c}\hfill \mathcal{A}\left(m,n\right)\circ \mathit{f}\left(\mathbf{q}\right)=\sum _{m+1\le \left|\mathbf{l}\right|\le m+n}{\left(-1\right)}^{m+n-\left|\mathbf{l}\right|}\left(\begin{array}{c}n-1\\ m+n-\left|\mathbf{l}\right|\end{array}\right)·\left({\mathcal{U}}_{l_1}^{\left(1\right)}\otimes \cdots \otimes {\mathcal{U}}_{l_n}^{\left(1\right)}\right)\circ \mathit{f}\left(\mathbf{q}\right)\end{array}$$

(28)

The expression before the tensor product is a binomial coefficient, which is defined as follows:

$$\left(\begin{array}{c}n-1\\ m+n-\left|\mathbf{l}\right|\end{array}\right)=\frac{\left(n-1\right)!}{\left(m+n-\left|\mathbf{l}\right|\right)!\left(\left|\mathbf{l}\right|-m-1\right)!}$$

(29)

The dimension independent maximum sparse grid level m controls the number of function evaluations and the associated accuracy of the PCE. The PCE is built through a linear combination of separate tensor polynomial chaos expansions for each underlying tensor quadrature grid [45]. Summation of the expansion terms is conducted with the Smolyak combinatorial coefficient in Equation (28). This improves accuracy in the coefficient estimation and preserves the consistency of the PCE [39].

In the present investigation, a sparse grid level $m=2$ is considered and closed fully nested Clenshaw–Curtis points are applied for quadrature. Univariate and bivariate effects in $\underline{\mathbf{R}}$ are modeled with orthogonal polynomials of highest-order $a=4$, resulting from the PCE construction according to [45]. Therefore, together with the GSA in Section 3.3, which shows the approximately linear or low nonlinear behavior of $\underline{\mathbf{R}}$, the results are assumed to be sufficiently accurate.

For the subsequent uncertainty analysis, a three-dimensional parameter space $n=3$ with the parameter vector ${\mathbf{q}}_{sparsegrids}$ contains the prioritized parameters in the following expression:

$${\mathbf{q}}_{sparsegrids}={\mathbf{q}}_{\mathrm{I}}=\left(\begin{array}{ccc}\vartheta & \chi & \varphi \end{array}\right)$$

(30)

The initial and boundary conditions are often not exactly known and are therefore subject to uncertainties. Next to the initial conditions, the boundary conditions enable the unambiguous solution of a partial differential equation regarding space and time. However, in case of the Navier–Stokes equations, inherent uncertainties in the boundary conditions arise due to the interaction of turbulent fluctuations with the boundary. For this reason, a normal distribution for uncertainty in the characteristic temperature difference through $\vartheta$ was defined according to the central limit theorem, which states that many independent random effects lead to normal distribution. The uncertainty of the wall tangential temperature gradient through $\varphi$ was defined to follow a log-normal distribution, which also in the course of the central limit theorem represents the distribution that results from the product of many positive independent random variables. The temperature gradient represents the temperature profile at the boundary that is established by the counterflow heat exchanger in the application case and since the supplied heat from the heat exchanger is always positive, the assumption of a log-normal distribution was made. In addition, the injection process creates large uncertainties in the actual build-up of helium stratification due to the turbulent flow of the gas stream through the injection nozzle. Hence, for the initial condition through $\chi$, a half normal distribution was also assumed because, if the stratification is stable, the density decreases in the vertical direction. The random variables for the realizations $\vartheta$, $\chi$, and $\varphi$ are $\mathsf{\Gamma }$, $\mathsf{\Lambda }$, and $\mathsf{\Phi }$. The random vector of parameters is $\mathcal{Q}=(\mathsf{\Gamma },\mathsf{\Lambda },\mathsf{\Phi })$. The standard deviation of $\mathsf{\Gamma }$ and $\mathsf{\Phi }$ was assumed to be ${\sigma }_{\mathcal{Q}}=0.1$. For $\mathsf{\Lambda }$, the standard deviation is defined with ${\sigma }_{\mathcal{Q}}=0.2$ since a large uncertainty in the structure of the stratification is assumed. $\mathcal{N}\left(\mu ,{\sigma }^2\right)$ and $\mathcal{LN}\left(\mu ,{\sigma }^2\right)$ denotes a normal distribution and a log-normal distribution with expectation $\mu$ and variance ${\sigma }^2$. $\mathcal{TN}\left(\mu ,{\sigma }^2,a,b\right)$ denotes a truncated normal distribution with a and b as the lower and upper bounds. The corresponding probability density functions $f_{\mathcal{Q}}\left(\mathit{q}\right)$ of the uncertain input parameters are then as follows:

$$\begin{array}{cccccc}\hfill \mathrm{\Gamma }& \sim \mathcal{N}\left(1,0.1^2\right)\hfill & \hfill \mathrm{\Lambda }& \sim \mathcal{TN}\left(0,0.2^2,0,1\right)\hfill & \hfill \mathsf{\Phi }& \sim \mathcal{LN}\left(0.1,0.1^2\right)\hfill \end{array}$$

(31)

After the propagation of these uncertainties, the probability density functions (PDF) $f_{\underline{\mathcal{R}}}\left(\underline{\mathbf{R}}\right)$, the expansion mean ${\mu }_{\underline{\mathcal{R}}}=\mathbb{E}\left(\underline{\mathcal{R}}\right)$, and the standard deviation ${\sigma }_{\underline{\mathcal{R}}}=\sqrt{\mathrm{var}\left(\underline{\mathcal{R}}\right)}$ for LES and URANS were determined. In Appendix B, the first-, second- and total-order Sobol indices are listed in

Table A3. First-order Sobol indices $S_{q_i}$ for the random response $\underline{\mathcal{R}}$.

	LES			URANS
$\underline{\mathbf{R}}$	${\mathit{S}}_{\mathbf{\vartheta }}$	${\mathit{S}}_{\mathbf{\chi }}$	${\mathit{S}}_{\mathbf{\varphi }}$	${\mathit{S}}_{\mathbf{\vartheta }}$	${\mathit{S}}_{\mathbf{\chi }}$	${\mathit{S}}_{\mathbf{\varphi }}$
${\widetilde{Nu}}_{left}$	0.8550	0.1083	0.0359	0.8392	0.1161	0.0441
${\widetilde{Nu}}_{right}$	0.8786	0.0473	0.0734	0.8471	0.0543	0.0978
${\widetilde{Nu}}_{bottom}$	0.0310	$7.4810\times {10}^{-4}$	0.9676	0.0319	$1.8460\times {10}^{-4}$	0.9673
${\widetilde{Nu}}_{top}$	0.0649	0.2382	0.6948	0.0558	0.2254	0.7175
${\widehat{Nu}}_{left}$	0.7741	0.1619	0.0048	0.7301	0.2086	0.0032
${\widehat{Nu}}_{right}$	0.7979	0.1106	0.0049	0.7179	0.1657	0.0050
${\widehat{Nu}}_{bottom}$	0.0213	0.0012	0.9768	0.0308	$2.4910\times {10}^{-4}$	0.9684
${\widehat{Nu}}_{top}$	0.1896	0.2295	0.5723	0.1962	0.2242	0.5736
${\tilde{E}}_k$	0.9454	0.0023	0.0493	0.9788	$3.3030\times {10}^{-4}$	0.0201
${\widehat{E}}_k$	0.8612	0.0147	0.0115	0.8508	0.0258	0.0067
$F{o}_{ϵ}$	0.8940	0.0343	0.0702	0.9027	0.0348	0.0624
$\tilde{I}$	0.0061	0.9094	0.0771	0.1270	0.7063	0.1624
$\widehat{I}$	0.6429	0.2137	0.0238	0.6275	0.2420	0.0190

,

Table A4. Second-order Sobol indices $S_{q_i{q}_j}$ for the random response $\underline{\mathcal{R}}$.

	LES			URANS
$\underline{\mathbf{R}}$	${\mathit{S}}_{\mathbf{\vartheta }\mathbf{\chi }}$	${\mathit{S}}_{\mathbf{\vartheta }\mathbf{\varphi }}$	${\mathit{S}}_{\mathbf{\chi }\mathbf{\varphi }}$	${\mathit{S}}_{\mathbf{\vartheta }\mathbf{\chi }}$	${\mathit{S}}_{\mathbf{\vartheta }\mathbf{\varphi }}$	${\mathit{S}}_{\mathbf{\chi }\mathbf{\varphi }}$
${\widetilde{Nu}}_{left}$	$5.3040\times {10}^{-4}$	$2.5950\times {10}^{-4}$	$4.6230\times {10}^{-5}$	0.0252	$4.4470\times {10}^{-4}$	$9.5290\times {10}^{-4}$
${\widetilde{Nu}}_{right}$	$4.2520\times {10}^{-4}$	$6.0350\times {10}^{-4}$	$8.2280\times {10}^{-5}$	0.0275	$4.9870\times {10}^{-4}$	0.0010
${\widetilde{Nu}}_{bottom}$	$4.2090\times {10}^{-5}$	0.0011	$8.2600\times {10}^{-5}$	$4.1120\times {10}^{-5}$	0.0010	$6.7580\times {10}^{-6}$
${\widetilde{Nu}}_{top}$	0.0010	$2.6460\times {10}^{-4}$	$2.8610\times {10}^{-4}$	0.0168	$4.9460\times {10}^{-4}$	$9.8280\times {10}^{-4}$
${\widehat{Nu}}_{left}$	0.0581	0.0016	0.0013	0.0861	0.0013	0.0025
${\widehat{Nu}}_{right}$	0.1028	0.0044	0.0046	0.1345	0.0043	0.0112
${\widehat{Nu}}_{bottom}$	$3.4510\times {10}^{-5}$	0.0014	$6.8360\times {10}^{-5}$	$2.3570\times {10}^{-4}$	$9.5770\times {10}^{-4}$	$3.6120\times {10}^{-5}$
${\widehat{Nu}}_{top}$	0.0011	0.0017	0.0052	0.0288	0.0052	0.0086
${\tilde{E}}_k$	$1.6030\times {10}^{-4}$	0.0023	$3.8430\times {10}^{-5}$	0.0604	0.0153	0.0036
${\widehat{E}}_k$	0.0635	0.0663	$4.6140\times {10}^{-4}$	0.0796	0.0207	0.0258
$F{o}_{ϵ}$	$3.3040\times {10}^{-4}$	$5.0210\times {10}^{-4}$	$1.0890\times {10}^{-4}$	0.1769	0.0053	0.0101
$\tilde{I}$	0.0028	0.0042	0.0017	0.0912	0.0070	0.0052
$\widehat{I}$	0.0991	0.0233	$8.1460\times {10}^{-4}$	0.1743	0.0052	0.0193

and

Table A5. Total-order Sobol indices $S_{T_{q_i}}$ for the random response $\underline{\mathcal{R}}$.

	LES			URANS
$\underline{\mathbf{R}}$	${\mathit{S}}_{{\mathit{T}}_{\mathbf{\vartheta }}}$	${\mathit{S}}_{\mathit{T}\mathbf{\chi }}$	${\mathit{S}}_{\mathit{T}\mathbf{\varphi }}$	${\mathit{S}}_{{\mathit{T}}_{\mathbf{\vartheta }}}$	${\mathit{S}}_{\mathit{T}\mathbf{\chi }}$	${\mathit{S}}_{\mathit{T}\mathbf{\varphi }}$
${\widetilde{Nu}}_{left}$	0.8557	0.1088	0.0362	0.8398	0.1164	0.0445
${\widetilde{Nu}}_{right}$	0.8792	0.0476	0.0738	0.8479	0.0545	0.0984
${\widetilde{Nu}}_{bottom}$	0.0316	$8.0080\times {10}^{-4}$	0.9682	0.0326	$1.9480\times {10}^{-4}$	0.9679
${\widetilde{Nu}}_{top}$	0.0664	0.2399	0.6958	0.0569	0.2263	0.7182
${\widehat{Nu}}_{left}$	0.8326	0.2201	0.0067	0.7864	0.2657	0.0060
${\widehat{Nu}}_{right}$	0.8791	0.1935	0.0140	0.8159	0.2725	0.0231
${\widehat{Nu}}_{bottom}$	0.0220	0.0013	0.9775	0.0313	$2.6020\times {10}^{-4}$	0.9689
${\widehat{Nu}}_{top}$	0.1946	0.2365	0.5776	0.1999	0.2293	0.5768
${\tilde{E}}_k$	0.9481	0.0028	0.0520	0.9795	$4.0750\times {10}^{-4}$	0.0209
${\widehat{E}}_k$	0.9730	0.0737	0.0659	0.9664	0.0819	0.0684
$F{o}_{ϵ}$	0.8953	0.0353	0.0708	0.9028	0.0349	0.0624
$\tilde{I}$	0.0113	0.9156	0.0804	0.1312	0.7095	0.1637
$\widehat{I}$	0.7622	0.3171	0.0403	0.7389	0.3337	0.0388

respectively. Thus, the uncertainties in the responses can be apportioned to the uncertainties in the input parameters. $f_{\underline{\mathcal{R}}}\left(\underline{\mathbf{R}}\right)$ was determined by sampling on the PCE approximation considering the input probability density functions $f_{\mathcal{Q}}\left(\mathit{q}\right)$. The mean ${\mu }_{\underline{\mathcal{R}}}$, standard deviation ${\sigma }_{\underline{\mathcal{R}}}$, and Sobol indices were computed analytically from the PCE coefficients ${\alpha }_j$.

When comparing the mean values and the standard deviations in

Table 3. Expansion mean ${\mu }_{\underline{\mathcal{R}}}$ and standard deviation ${\sigma }_{\underline{\mathcal{R}}}$ for the random response $\underline{\mathcal{R}}$.

	LES		URANS
$\underline{\mathbf{R}}$	${\mathbf{\mu }}_{\underline{\mathcal{R}}}$	${\mathbf{\sigma }}_{\underline{\mathcal{R}}}$	${\mathbf{\mu }}_{\underline{\mathcal{R}}}$	${\mathbf{\sigma }}_{\underline{\mathcal{R}}}$
${\widetilde{Nu}}_{left}$	47.7100	2.4600	48.6900	2.3070
${\widetilde{Nu}}_{right}$	−47.7800	2.2790	−48.6900	2.1500
${\widetilde{Nu}}_{bottom}$	14.9900	2.2560	15.2400	2.3140
${\widetilde{Nu}}_{top}$	−9.5170	2.2030	−9.9850	2.1100
${\widehat{Nu}}_{left}$	5.2100	2.4790	4.6910	2.1640
${\widehat{Nu}}_{right}$	4.2710	2.2060	4.1000	1.8390
${\widehat{Nu}}_{bottom}$	15.3600	2.2520	15.5400	2.3200
${\widehat{Nu}}_{top}$	9.5870	2.1930	10.0400	2.1270
${\tilde{E}}_k$	$4.7310\times {10}^{-4}$	$6.9620\times {10}^{-5}$	$4.2430\times {10}^{-4}$	$5.4710\times {10}^{-5}$
${\widehat{E}}_k$	$1.2910\times {10}^{-4}$	$8.2110\times {10}^{-5}$	$1.0430\times {10}^{-4}$	$5.9700\times {10}^{-5}$
$F{o}_{ϵ}$	0.0579	0.0054	0.0606	0.0062
$\tilde{I}$	0.4263	0.0118	0.4482	0.0118
$\widehat{I}$	0.0444	0.0304	0.0451	0.0302

, it is noticeable that the results for ${\mu }_{\underline{\mathcal{R}}}$ and ${\sigma }_{\underline{\mathcal{R}}}$ largely coincide for LES and URANS. With $\langle ·\rangle$ and $\left|·\right|$ denoting the mean and absolute value, respectively, the difference can be quantified in summary by the mean relative deviation of URANS from LES regarding the statistical moments in Table 3 with $〈\left|{\mu }_{\underline{\mathcal{R}}}^{URANS}-{\mu }_{\underline{\mathcal{R}}}^{LES}\right|/{\mu }_{\underline{\mathcal{R}}}^{LES}〉$ and $〈\left|{\sigma }_{\underline{\mathcal{R}}}^{URANS}-{\sigma }_{\underline{\mathcal{R}}}^{LES}\right|/{\sigma }_{\underline{\mathcal{R}}}^{LES}〉$, which are $0.055$ and $0.092$ respectively. From this, good agreement of both approaches becomes visible. However, when considering the probability density functions, correspondence and deviations between LES and URANS become clear in more detail. A part of the interpretation of the results is based on the fact that, if the first-order Sobol index $S_{q_i}$ is large with respect to an input variable and the response function additionally shows linear behavior regarding this parameter, then the result variable approximately shows a similar distribution type ${\mathcal{I}}_{\mathcal{Q}}$, as defined for the input variable. Considering the affine model for the approximation of a single response (Equation (32)), for which the variance is predominantly determined by one input variable, the shape of the probability density function after propagation remains approximately the same whereas the expectation and variance change. From this, conclusions can be made and the plausibility of the results can be checked.

$$\underline{\mathcal{R}}\approx A\mathcal{Q}+B\sim {\mathcal{I}}_{\mathcal{Q}}\left(A{\mu }_{\mathcal{Q}}+B,A^2{\sigma }_{\mathcal{Q}}^2\right)$$

(32)

Starting with ${\widetilde{Nu}}_{left}$ and ${\widetilde{Nu}}_{right}$ in Figure 13a,b, it is noticeable that the shape of the PDF is similar to a normal distribution. This is the case because the variance is predominantly caused by the temperature specification at the left and right walls via $\vartheta$, as can be seen from the first-order Sobol indices $S_{q_i}$ in Table A3, and shows a nearly linear behavior that becomes visible due to the low interaction behavior through the second-order Sobol indices $S_{q_i{q}_j}$ in Table A4 and a low ${\sigma }_{q_i}$ in Figure 9e,f. The distribution for ${\widetilde{Nu}}_{bottom}$ tends towards a log-normal distribution. As can be seen in Table A3, the variance is predominantly determined by the wall-tangential temperature gradient via $\varphi$, which was specified with a log-normal distribution. For this reason, a similar distribution for the result originates. For ${\widetilde{Nu}}_{top}$, the effects of $\chi$ also become visible. Increasing $\chi$ enlarges the amount of helium at the top wall. Due to the lower density of helium, this attenuates the erosion at the top wall. This results in a lower convective heat transfer in this region. This creates the superimposed distribution based on the effects of $\varphi$ and $\chi$, shown in Figure 13d. When comparing LES and URANS, the probability density distributions in Figure 13 and the statistical moments in Table 3 are in very good agreement. If taking LES as a reliable reference, the absolute value of the Nusselt number for URANS is slightly overestimated on the left, right, and top walls.

The distributions for ${\widehat{Nu}}_{left}$ and ${\widehat{Nu}}_{right}$ in Figure 14a,b, approximately resemble the shape of a chi-distribution because the response function evaluates the absolute difference between the mixing transients. Analogous to $\widetilde{Nu}$, the variance primarily arises through the variance in the temperature at the left and right walls, which is normally distributed and dominates the distribution of the response. This can be seen in Table A3 by means of the first-order Sobol indices $S_{q_i}$. The distributions of URANS and LES are in good agreement. The distributions for the top and bottom walls in Figure 13c,d and Figure 14c,d are nearly the same as that for $\widetilde{Nu}$ because the top and bottom walls are adiabatic for the reference case, and this results in the integral deviation from the abscissa.

In Figure 15a, a normal distribution also arises for the PDF of $\widetilde{E_k}$, as again the temperature on the left and right walls via $\vartheta$ causes the major variance, as can be seen with $S_{q_i}$ in Table A3. The global kinetic energy ${\tilde{E}}_k$ is underestimated with URANS and leads to longer mixing times $F{o}_{ϵ}$ in Figure 16a since the mixing process is dominated by the arising natural convection. The distribution for LES is wider than for URANS because LES captures most of the turbulent structures. The larger amount of the standard deviation also explains the larger deviation from the reference case ${\widehat{E}}_k$ for LES in Figure 15b. By taking the LES as reference, the results for ${\widehat{E}}_k$ indicate that URANS is able to predict the variability of $E_k$ at $Ra=2\times {10}^9$ due to the parameter uncertainties.

In Figure 16a, the shape of the PDF for $F{o}_{ϵ}$ resembles a normal distribution due to the large influence of $\vartheta$ for both LES and URANS, as can be seen with $S_{q_i}$ in Table A3. Low values for ${\sigma }_{q_i}$ in Figure 12d and $S_{q_i{q}_j}$ in Table A4 indicate the approximately linear behavior of $F{o}_{ϵ}$. In comparison to URANS, LES predicts shorter mixing times and causes a smaller width of the distribution, which also can be measured by the smaller mean value ${\mu }_{\underline{\mathcal{R}}}$ and standard deviation ${\sigma }_{\underline{\mathcal{R}}}$ in Table 3. The distributions for $\tilde{I}$ in Figure 16b tend towards larger values, as can be seen especially for the LES results. This tendency arises from the truncated normal distribution for the linear stratification parameter $\chi$, which increases the segregation intensity and causes the largest proportion of variance, as indicated by $S_{q_i}$ in Table A3. The variance component including the smaller values mainly arises from the remaining parameters $\vartheta$ and $\varphi$. In contrast to the larger integral mean ${\mu }_{\underline{\mathcal{R}}}$ of the segregation intensity $\tilde{I}$ in URANS, the results indicate that LES has a higher mixing intensity $\tilde{\mathcal{M}}$ whereas the width of the distributions is almost identical. The shorter mixing time measured by $F{o}_{ϵ}$ for LES and the clear difference in $\tilde{I}$ between LES and URANS arises since the more resolved convection in LES, which might not be represented sufficiently accurate by the turbulence model in URANS, results in a more accurate prediction of mass transfer, since anisotropic large-scale eddies have a large influence on the mixing process. The deviation from the reference case $\widehat{I}$, which provides a good measure for the variability of the results when the parameters change, shows very good agreement between LES and URANS.

4. Conclusions

In this paper, a methodology for the assessment of uncertainties of Large Eddy Simulation (LES) and Unsteady Reynolds-Averaged Navier–Stokes (URANS) for the buoyancy-driven mixing process between two miscible fluids at $Ra=2\times {10}^9$ was presented, starting from the preceding mesh convergence study and sensitivity analysis to the actual uncertainty quantification. For this purpose, an initial stratification with 40 $\mathrm{vol}\%$ helium was defined in the upper third of the Differentially Heated Cavity (DHC) in addition to air. From the investigation, the parameters that have significant influence on the mixing process were identified and it was found that there is good agreement between LES and URANS for $Ra=2\times {10}^9$ regarding the final stochastic results.

The mesh convergence study provides an approach for meshing the fluid domain by means of linear expansion factors, which allows for systematic determination of the required wall normal resolution with the number of cells in the boundary layer ${\overline{n}}_{{\delta }_{u/ip}}$ and ${\overline{n}}_{{\delta }_{T/grad}}$ and the required wall tangential resolution with the dimensionless cell widths $\mathsf{\Delta }{x}_{\Vert }^{+}$ and $\mathsf{\Delta }{y}_{\Vert }^{+}$, which can be transferred to other engineering applications underlying similar physics. The applied convergence criteria $\underset{\mathrm{LES},\mathrm{DNS}}{\mathrm{MAE}}\left(\mathrm{avg}\left(Nu\right)\right)/{\langle \mathrm{avg}\left(N{u}_{DNS}\right)\rangle }_y\le {ϵ}_{MAE}\approx {10}^{-2}$ containing the relative mean average error (MAE) between the results for the local Nusselt number received from LES and DNS ensures sufficiently accurate reflection of the global characteristics. This condition is fulfilled from a wall tangential resolution of $\mathsf{\Delta }{x}_{\Vert }^{+},\mathsf{\Delta }{y}_{\Vert }^{+}=30$ and boundary layer resolution of ${\overline{n}}_{{\delta }_{u/ip}}=7.87$ and ${\overline{n}}_{{\delta }_{T/grad}}=10.10$.

As a prelimary step for the uncertainty analysis, the number of relevant parameters was reduced by using a one-at-a-time (OAT-) method, in which all parameters were varied step by step starting from an initial point. It was found that the thermal boundary conditions, such as the temperature of the enclosing walls and the wall tangential temperature gradient at the vertical walls, are of essential importance, as the mixing process is primarily driven by buoyancy effects. In addition, the initial build-up of the helium stratification shows significant influence on the mixing process. Additionally, the diffusion coefficient for mass transfer causes large changes. Therefore, these parameters were examined in more detail by global sensitivity analysis using the Morris method to show more precise statements about the sensitivity of the response function to the different parameters and possible nonlinear effects or interaction behavior. In summary, it can be seen that the sensitivity behavior of LES and URANS is qualitatively comparable with moderate quantitative deviations at $Ra=2\times {10}^9$.

The propagtion of the uncertainty in selected parameters was eventually investigated with Polynomial Chaos Expansions (PCE). For the determination of the expansion-coefficients, the spectral projection approach was applied. The multi-dimensional integration was conducted by the Smolyak sparse grid method. The probability density functions (PDF) and Sobol indices elucidate slight differences between LES and URANS. The heat transfer, which is measured by the Nusselt number, is well predicted with URANS in comparison to LES. However, the global kinetic energy and associated convective mass transport shows differences between LES and URANS. This results in different mixing times and mixing intensity. If LES can be regarded as a suitable reference, URANS with the currently used turbulence model predicts a slightly incorrect mixing behavior and small deviations in the LES for the variability of the result variables exist. The deviations are probably caused by the greater limitation of degrees of freedom of the URANS equations, which complicates the accurate prediction of internal buoyancy-driven mixing processes. Turbulence has a decisive influence on the mass transport and deviations directly imply effects on the natural convection flow. According to Section 3.4, the mean relative deviation between LES and URANS for the statistical moments is small. While the PDFs indicate that discernable differences arise between the uncertainty analysis based on URANS and LES, there is still reasonable agreement in most of the statistical moments. In summary, the resource-efficient URANS is suitable as a lower-fidelity model for the prediction of the input uncertainties propagation at $Ra=2\times {10}^9$ for the investigated mixing process based on the reference analysis conducted with LES. Future work will extend the investigations to higher Rayleigh numbers comparable to the THAI application as well as the THAI large-scale facility itself.