## Abstract

This work investigates the importance of verification and validation (V&V) to achieve predictive scale-resolving simulations (SRSs) of turbulence, i.e., computations capable of resolving a fraction of the turbulent flow scales. Toward this end, we propose a novel but simple V&V strategy based on grid and physical resolution refinement studies that can be used even when the exact initial flow conditions are unknown, or reference data are unavailable. This is particularly relevant for transient and transitional flow problems, as well as for the improvement of turbulence models. We start by presenting a literature survey of results obtained with distinct SRS models for flows past circular cylinders (CCs). It confirms the importance of V&V by illustrating a large variability of results, which is independent of the selected mathematical model and Reynolds number. The proposed V&V strategy is then used on three representative problems of practical interest. The results illustrate that it is possible to conduct reliable V&V exercises with SRS models and evidence the importance of V&V to predictive SRS of turbulence. Most notably, the data also confirm the advantages and potential of the proposed V&V strategy: separate assessment of numerical and modeling errors, enhanced flow physics analysis, identification of key flow phenomena, and ability to operate when the exact flow conditions are unknown or reference data are unavailable.

## 1 Introduction

Scale-resolving simulation (SRS) models are nowadays widely applied in computation of turbulent flow problems of practical interest, e.g., vehicle aerodynamics and hydrodynamics [1–7], offshore engineering [8–13], propeller design [14–19], materials mixing [20–25], and combustion [26–31]. These formulations are characterized by their ability to resolve a fraction of or all turbulence scales present in a flow. This property unleashes the potential of SRS methods to obtain high-fidelity predictions of complex flows not amenable to modeling with the Reynolds-averaged Navier–Stokes (RANS) equations due to their limitations representing transient phenomena, onset and development of turbulence, instabilities, and coherent structures. Among all SRS formulations, hybrid [32–37], bridging [38–44], and large-eddy simulation (LES) [45–47] are engineering practitioners' most usual modeling strategies.

The further establishment of SRS methods in many areas of science and engineering is inevitably dependent on the ability to obtain predictive computations, i.e., simulations with a quantified and acceptable level of numerical, input, and modeling uncertainty [48,49]. Yet, the estimation of these sources of computational uncertainty is difficult in SRS and, as such, often neglected. There are four primary reasons for this:

*Complexity*—SRS models calculate instantaneous variables which are closely dependent on the physical resolution or range of scales resolved by the model. This hampers the development of accurate and robust SRS turbulence closures by making the comparison and interpretation of some quantities difficult, as it is not possible to match the physical resolution of the reference experiments or simulations.*Cost*—SRS calculations are usually computationally more intensive than RANS owing to the numerical requirements (spatiotemporal grid resolution, simulation time, iterative convergence criterion, etc.) necessary to accurately resolve a fraction of the turbulence field. This makes detailed verification and validation (V&V) exercises time consuming and resource intensive.*Dynamic modeling*—The physical resolution of most SRS formulations is defined based on the grid resolution and the simulated flow properties to optimize the use of the available computational resources. However, such modeling option makes the formulation's physical resolution vary upon grid refinement. This turns SRS models grid dependent, precluding the evaluation of modeling and numerical errors separately, a requirement for reliable V&V exercises. Also, dynamic SRS formulations lead to commutation errors, since the filtering operator of the model does not commute with spatial and temporal differentiation [50].*Numerical setup*—A precise validation exercise requires numerical simulations and reference experiments performed on the same flow conditions, i.e., domain, boundary conditions, material properties, and initial or inflow conditions. For SRS, it is not enough to know the mean values of these quantities; initializing the model requires both the time-resolved large-scale motions and the detailed statistics of the small scales. However, there are many cases where some of this information is not available from the reference experiments. This is particularly important for transient and transitional flows and may significantly affect proper comparisons between simulations and reference experiments.

To illustrate the role of V&V exercises [48,51–53] to numerical prediction and turbulence modeling, Fig. 1 depicts the outcome of a literature survey [49,54] that compiles results for the time-averaged drag coefficient obtained for the flow around a circular cylinder (CC) at different Reynolds numbers, Re. One set of experimental measurements is included for comparison [55], and the numerical results are colored by category of model: direct numerical simulation (DNS), implicit large-eddy simulation (ILES), LES, the bridging partially averaged Navier–Stokes (PANS) equations, hybrid, and RANS. Note that except for RANS, all these formulations are SRS models. The numerical results exhibit an enormous range of values, which does not appear to improve for any specific mathematical model. Nor does it converge with physical resolution, from RANS where turbulence is completely modeled, to LES/ILES where most turbulent scales are resolved. Considering the case $Re=3900$, the collected $C\xafD$ results vary from 0.6 to 1.8, representing a maximum variation of 191% (taking $C\xafD=0.6$ as the reference). Overall, Fig. 1 emphasizes the need for V&V exercises to further establish SRS methods in engineering and generate confidence in their numerical results when reference solutions are unavailable.

This study investigates the importance of V&V exercises to achieve predictive SRS of turbulence. This is a necessary step to diminish the variability of SRS results in Fig. 1 and gain confidence on such formulations. Toward this end, we propose a novel but simple V&V approach based on grid and physical resolution refinement studies that can be used even when the initial flow conditions are not sufficiently characterized, or reference results are unavailable. This is crucial for transient and transitional flow simulations and for the improvement of the accuracy of turbulence models. The proposed method is utilized to ascertain the simulations' fidelity of three representative test-cases involving transition to turbulence: the flow around a CC at $Re=3.9\xd7103$, the Taylor–Green vortex (TGV) at $Re=3000$ [56], and the Rayleigh–Taylor (RT) flow [57,58] at Atwood number $At=0.5$. The first problem is a statistically unsteady flow with well-characterized flow conditions and available experimental results and uncertainties, enabling the estimation of validation uncertainties. The second test-case is a transient flow with exact flow conditions and available reference DNS results. Yet, the numerical uncertainty of these DNS results has not been estimated. The third problem is a transient flow without reference experimental results due to the complexity of characterizing the exact experimental RT flow conditions (see literature survey of Pereira et al. [59]). These three test-cases constitute the ideal validation space for the new V&V method. All computations are based on the PANS equations [23,42,60] at constant physical resolution to enable the assessment of numerical and modeling errors separately and prevent commutation errors.

## 2 Verification and Validation Strategy

Therefore, unlike an error, uncertainty is always a positive quantity.

Here, $Ei(\varphi )$ is the input or parameter error caused by inexact values of boundary conditions, material properties, or initial conditions, $Em(\varphi )$ is the modeling error defined as the difference between the exact solution of the mathematical model and physical truth, and $En(\varphi )$ is the numerical error which can be divided into a discretization, iterative, round-off, and statistical (unsteady flows like the CC case) component [48,51–53].

The estimation of $Uv(\varphi )$ can be quite challenging for problems in which the initial flow conditions are not well characterized (large $Ui(\varphi )$). Also, the calculation of $Uv(\varphi )$ is not possible in problems where reference experiments are unavailable. This is often the case for such problems as transitional flow, high-speed combustion, climate prediction, and full-scale problems. In addition to these difficulties, the governing equations of most SRS models directly depend on the grid resolution. This feature prevents the separate and reliable quantification of numerical and modeling errors ($En(\varphi )$ and $Em(\varphi )$) until all flow scales are resolved in a DNS.

Therefore, the numerical prediction of complex flows with SRS models raises the following question: how can one conduct reliable V&V exercises to evaluate the accuracy of the computations? We propose the following strategy. Let us start by considering an SRS formulation whose governing equations do not depend on the grid resolution. This requires a model that uses an explicit type of filtering operator. Hence, one can evaluate modeling and numerical errors separately, and the conditions necessary to derive SRS models and assure scale-invariance [61] are satisfied. It also prevents commutation errors [50]. In this work, we use the PANS equations [42] supplemented with turbulence closures for incompressible single fluid [60,62] and variable-density flow [23]. We emphasize that we could have selected other SRS formulations as long as they use an explicit closure to represent the unresolved flow scales that can operate at constant physical resolution. This might require model modifications. For example, traditional SRS models utilizing a local characteristic grid size $\Delta SRS$ to set the physical resolution (e.g., detached eddy simulation (DES), very large simulation (VLES), extra large eddy simulation (XLES), LES, etc.) can be modified to either use a uniform $\Delta SRS$ or the $\Delta SRS$ distribution of the coarsest grid used in the study. This strategy removes the governing equations dependence from the grid cell size, is simple to implement, and enables decoupling numerical from modeling errors. Alternatively, one can use explicit filtering. Yet, the first approach is better suited for practical simulations of turbulence.

where *k _{u}* is the unresolved or modeled turbulence kinetic energy, and

*k*is the total. Note that

_{t}*f*= 0 is equivalent to DNS because all turbulence scales are resolved, whereas

_{k}*f*= 1 is equivalent to RANS since the turbulence closure models all turbulence scales. For a better physical understanding of the impact of

_{k}*f*on the predicted flow fields, Fig. 2 illustrates how the computed velocity field evolves with the relative filter length size,

_{k}*n*.

*n*equal to one is a direct numerical simulation (

*f*= 0) where all flow scales are resolved, leading to a highly unsteady velocity field with steep gradients that increase the cost of numerical simulations substantially. In contrast,

_{k}*n*= 349 models most of the turbulence field ($fk\u22480.96$), which decreases the cost of the computations.

This class of SRS formulations provides the background to propose the new V&V strategy based on grid and physical resolution refinement exercises. It enables the evaluation of the numerical accuracy of the simulations through grid refinement exercises performed at constant *f _{k}*, and the assessment of the modeling accuracy of the calculations through physical resolution refinement studies since the solutions converge toward the “truth” as $fk\u21920$. Yet, the solutions are expected to converge at coarser physical resolutions than the limiting case of

*f*= 0 (DNS for a continuous medium). This is shown below and in numerous studies [12,59,60,64–68], pointing out the method's potential to be applied to complex new problems without experimental solutions. The latter feature is a unique feature of the new V&V method, which can be crucial to many complex flow problems.

_{k}and/or the physical resolution, *f _{k}*. In relation (7),

*N*

_{1}and

*N*are the number of cells of the finest grid available (

_{i}*i*= 1) and a selected grid

*i*. The data can be assessed and visualized in three ways: (i) fixing

*f*and plotting $\varphi (ri)$; (ii) fixing

_{k}*r*and depicting $\varphi (fk)$; and (iii) plotting $\varphi (fk,ri)$. This last approach provides a more comprehensive analysis, since it enables the assessment of the numerical and modeling accuracy, their interaction and interdependency, and the interpretation of the predicted flow physics and its dependence on

_{i}*f*(see Sec. 4). Figure 3 shows that the solutions of $\varphi $ converge both with the grid and physical resolutions as $ri\u21920$ and $fk\u21920$. Note that the smallest

_{k}*r*shown throughout this work is equal to one, corresponding to the solution obtained on the finest grid tested.

_{i}*r*= 0 would require the extrapolation of the solution to an infinitely fine grid (e.g., Refs. [51] and [69]). From the interpretation of Fig. 3, it is possible to infer that the present technique can be utilized when the experimental flow conditions are poorly characterized or unknown or in cases in which the reference solution is unavailable. In such cases, one can use the extrapolated solution ($ri\u21920$ and $fk\u21920$) as reference. It is important to mention that the former procedure could be “generalized” to cases of grid dependent SRS models by considering the cross-derivatives with respect to

_{i}*r*and

_{i}*f*, $(\u22022\varphi )/(\u2202ri\u2202fk)$.

_{k}## 3 Test-Cases and Numerical Details

### 3.1 Flow Problems.

We use three representative flow problems to investigate the importance of V&V to SRS of turbulence and evaluate the potential of the proposed V&V strategy. The first problem is a statistically unsteady flow with well-characterized flow conditions and available experimental results and uncertainties, enabling the estimation of validation uncertainties. The second test-case is a transient flow with exact flow conditions and available reference DNS results. Yet, the numerical uncertainty of these DNS results has not been estimated. The third problem is a transient flow without reference experimental results due to the complexity of characterizing the exact experimental initial flow conditions. These three test-cases constitute the ideal validation space for the new V&V method and are described below.

#### 3.1.1 Circular Cylinder.

The flow past a CC at $Re=3900$ is an archetypal problem commonly used to assess the prediction of transition, turbulence, and separation around bluff-flows in the subcritical regime [79,80]. This problem is typical of model-scale applications and unmanned vehicles. The flow is statistically unsteady and characterized by a laminar boundary-layer, flow separation leading to a free shear-layer, onset and development of turbulence in the free shear-layer, and finally a turbulent wake. Also, the flow physics is closely dependent on two coherent structures: the Kelvin–Helmholtz in the free shear-layer and the large-scale vortex-shedding in the wake. The first is responsible for transition to turbulence in the free shear-layer. The spatial development of the flow in this regime comprises four major steps, proper representation of which determines the accuracy of any numerical simulation: $(i)$ onset of the Kelvin–Helmholtz instability in the free shear-layer; $(ii)$ spatial development of the Kelvin–Helmholtz rollers; $(iii)$ breakdown to high-intensity turbulence; and $(iv)$ turbulent free shear-layer roll-up, leading to vortex-shedding. A comprehensive numerical analysis of this problem is given in Refs. [60], [64], and [65].

In this work, we use the experimental measurements of Parnaudeau et al. [81] as reference to evaluate the modeling error. These have been conducted in a wind-tunnel with a square cross section of width equal to $23.3D$ (*D* is the cylinder's diameter), and a uniform inflow characterized by a turbulence intensity smaller than $I=0.2%$. The Reynolds number, $Re\u2261VoD/\nu $, is equal to 3900 (*V _{o}* is the inflow velocity, and

*ν*is the kinematic viscosity), and the sampling period used to compute the flow statistics is $\Delta TVo/D=2.08\xd7105$. The quantity of interest extracted from this study is the time-averaged velocity field, $V\xafi$. This has been measured using particle image velocimetry, and the reported experimental uncertainties do not exceed 1.0%. We also used the experiments of Norberg [82,83] to assess the predictions of the time-averaged drag coefficient, $C\xafD$, root-mean-square lift coefficient, $CL\u2032$, and pressure coefficient on the cylinder's surface, $C\xafp(\theta )$. Compared to Parnaudeau et al. [81], these experimental studies are conducted with two slightly different setups: $Re=4000,\u2009215.5D\xd780.0D$ cross section, $I<0.06%$, and $\Delta TVo/D=2.10\xd7104$ for Ref. [82], whereas $Re=4400,\u2009314.1D\xd7105.0D$ cross section, and $\Delta TVo/D=3.20\xd7104$ for Ref. [83]. The reported experimental uncertainties are 0.01 for $C\xafp(\theta )$, 1.0% for $C\xafD$, and 4.0% for $CL\u2032$.

The numerical simulations try to replicate the experimental apparatus of Parnaudeau et al. [81] closely. The computational domain is shown in Fig. 4(a). The domain is 50*D* long, with the inlet at 10*D* from the inflow, 24*D* wide, and 3*D* thick, with cylinder axis centered. The impact of the cylinder's length ($L=3D$) has been ascertained and guarantees acceptable input uncertainties [54]. At the inlet ($x1/D=\u221210$), the prescribed boundary conditions are constant velocity (*V _{o}*) aligned with the streamwise direction

*x*

_{1}, and the pressure is extrapolated from the interior of the domain. The turbulence quantities are set to result in a turbulence intensity $I=0.2%$ and a ratio between turbulent and molecular kinematic viscosity equal to $10\u22123$. At the outlet ($x1/D=40$), all streamwise derivatives are set equal to zero, whereas no-slip and impermeability conditions are applied on the cylinder's surface. Also, the normal pressure gradient and the turbulence quantities are equal to zero. The exception is the specific dissipation which is given at the center of the nearest-wall cell by $\omega =80\nu d\u22122$ (

*d*is the nearest-wall distance). The pressure and transverse derivatives of the remaining dependent quantities are equal to zero at $x2/D=\xb112$, and symmetry conditions are prescribed at $x3/D=0$ and 3. The simulations are initialized from a 200 time units ($\Delta TVo/D$) computation and run for a minimum of 500 time units. The flow statistics are calculated using a sampling period $\Delta TVo/D\u2265350$ [49].

#### 3.1.2 Taylor–Green Vortex.

*V*is the initial velocity magnitude. The corresponding pressure field is obtained from solving the Poisson equation

_{o}Here, *P _{o}* and

*ρ*are the pressure and density magnitude at

_{o}*t*= 0. After

*t*= 0, these structures interact and deform, leading to the formation of vortex sheets. Afterward, the structures get closer and roll-up, causing the onset of turbulence due to a complex vortex reconnection process between pairs of counter-rotating vortices [68]. In the next instants, the remaining vortical structures breakdown, and turbulence further develops. This leads to rapid dissipation of the flow kinetic energy. Since the TGV is an isolated system, the flow kinetic energy decays in time.

The flow configuration of the simulations replicates that used in the DNS studies of Brachet et al. [84] and Drikakis et al. [85] at $Re\u2261VoLo/\nu =3000$. The simulations are conducted in a cubical computational domain with width $L=2\pi Lo$, where periodic conditions are applied on all boundaries. The initial Mach number is set $Mao=0.28$, and the simulations run until $t=20Vo/Lo$. The initial thermodynamic and flow properties are the following: vortex velocity $Vo=104$ cm/s, box size $Lo=1.00$ cm, density $\rho o=1.178\xd710\u22123$ g/cm^{3}, mean pressure $Po=105$ Pa, dynamic viscosity $\mu =3.927\xd710\u22123$ g/(cm s), heat capacity ratio $\gamma =1.40$, turbulence kinetic energy $ko=10\u22127$ cm^{2}/s^{2}, and dissipation length scale $So=6.136\xd710\u22123$ cm. The DNS studies of Brachet et al. [84] at $Ma=0$ (incompressible) and Drikakis et al. [85] at $Mao=0.28$ are used to evaluate the modeling accuracy of the simulations. These do not report numerical uncertainties. We emphasize that the differences in $Mao$ between our simulations and Brachet et al. [84] computations may lead to discrepancies.

#### 3.1.3 Rayleigh–Taylor.

The RT [57,58] is a benchmark test-case of variable-density flow. This problem is a fundamental test-case for flows with material mixing and turbulence generation by baroclinicity. The flow initially consists of two fluids of different densities separated by a perturbed interface. These materials are at rest, with the dense medium (*ρ _{h}*) located on top of the light material (

*ρ*) as shown in Fig. 4(c). Immediately after

_{l}*t*= 0, the heavy fluid accelerates downward by the action of a body force, whereas the light fluid moves upwards. The interface perturbations create a misalignment between the density gradient and the pressure, which induces the RT instability. This generates structures: penetration of heavy fluid into the light fluid (called bubbles) and, conversely, spikes by penetration of light fluid into the heavy fluid. The boundaries of these structures experience shearing, resulting in a Kelvin–Helmholtz instability, which triggers the onset of turbulence. Afterward, turbulence further develops, enhancing the mixing of the two materials. The density difference is characterized by the Atwood number, $At=(\rho h\u2212\rho l)/(\rho h+\rho l$).

The perturbations' wavelengths range between modes 30 and 34 ($30\u2264n2+m2\u226434$), and their amplitude is featured by a maximum standard deviation equal to 0.04*L _{o}*. In Eq. (12), the modes

*m*and

*n*are selected to include the most unstable mode of the linearized problem, and

*r*

_{1}and

*r*

_{3}are random numbers between 0 and 1. The present numerical experiments use the ideal gas equation of state, and the initial temperature is set to maintain the flow $Ma<0.10$. The initial thermodynamic and flow properties are defined as follows: $\mu l=0.002$ g/(cm s), $\mu h=0.006$ g/(cm s), $\rho l=1.0$ g/cm

^{3}, $\rho h=3.0$ g/cm

^{3}, $\gamma l=\gamma h=1.40$,

*g*= −980 cm/s

^{2}, $ko=10\u22126cm2/s2,\u2009So=10\u22126cm$, and Schmidt and Prandtl numbers equal to one.

### 3.2 Numerical Settings.

The CC flow is calculated with the incompressible flow solver refresco [87]. This is a community-based open-usage flow solver optimized for maritime applications. refresco solves multiphase incompressible viscous flows using the filtered continuity and Navier–Stokes equations, complemented with turbulence closures. The equations are discretized using a finite volume approach with cell-centered collocated variables, in strong-conservation form, and a pressure-correction equation based on the semi-implicit method for pressure linked equations (SIMPLE) algorithm is used to ensure mass conservation. Time integration is performed implicitly with a second-order backward scheme. The spatial discretization schemes are also second-order accurate, and we use a second-order upwind-based scheme to discretize the convective terms of the governing equations. Four spatiotemporal grid resolutions are used, ranging from $1.04\xd7106$ to $4.55\xd7106$ cells and time-steps between $8.53\xd710\u22123$ and $5.21\xd710\u22123$ time units. These lead to a refinement ratio *r* of approximately 1.64. At each implicit time-step, the nonlinear system for velocity and pressure is linearized with Picard's method, and a segregated approach is adopted for the solution of all transport equations. The simulations run on double precision, and the iterative convergence criterion [49] guarantees that the maximum normalized residual of all dependent quantities is always smaller than $10\u22125$.

The implementation is face-based, which permits grids with elements consisting of an arbitrary number of faces. refresco was initially developed as a RANS solver [88], where the Reynolds-stress tensor was modeled through linear turbulent viscosity closures. Pereira and co-workers [54,89] extended the available mathematical models and turbulence closures to PANS, XLES, delayed DES (DDES), improved DDES (IDDES), RANS-explicit algebraic Reynolds stresses model (EARSM), and RANS-Reynolds stresses model (RSM) methods. The code is parallelized using message passing interface (MPI) and subdomain decomposition and runs on Linux workstations and high performance computing (HPC) clusters. The code is currently being developed, verified, and tested at MARIN (the Netherlands) in collaboration with several universities [87].

*c*is the speed of sound, and $\Delta x$ is the grid size. The CFL number is set equal to 0.45 for the TGV and 0.50 for the RT. The code can utilize an adaptive mesh refinement algorithm for following waves, especially shock-waves and contact discontinuities, but this option is not used in this work to prevent hanging-nodes [94] and, as such, the simulations use orthogonal uniform hexahedral grids. These have 128

^{3}, 256

^{3}, and 512

^{3}(and 1024

^{3}for $fk=0.00$) cells for the TGV and $642\xd7192,\u20091282\xd7384,\u2009and\u20092562\xd7768$ (and $5122\xd71536$ with CFL $=0.8$ for $fk=0.00$) cells for the RT.xrage models miscible material interfaces and high convection-driven flows with a van Leer limiter [95], without artificial viscosity, and no material interface treatments [96,97]. The solver uses the assumption that cells containing more than one material are in pressure and temperature equilibrium as a mixed cell closure. The effective kinematic viscosity in multimaterial problems [24] is defined as

where *n* is the material index, *n _{t}* is the number of materials, and

*f*is the volume fraction of material

_{n}*n*. For the RT flow, the diffusivity $D$ and thermal conductivity

*κ*are defined by imposing the Schmidt ($Sc\u2261\nu /D$) and Prandtl ($Pr\u2261cp\mu /\kappa $) numbers equal to one.

## 4 Results

### 4.1 Circular Cylinder.

We initiate the discussion of the results with the flow past a circular cylinder in the subcritical regime. Figure 5 presents change with grid, *r _{i}*, and physical resolution,

*f*, refinement in three time-averaged quantities: the drag coefficient, $C\xafD$, base pressure coefficient, $C\xafpb$, and recirculation length, $L\xafr$, and one root-mean-square quantity, RMS lift coefficient, $CL\u2032$. As expected, the simulations converge both with numerical resolution ($ri\u21920$) and physical resolution ($fk\u21920$). The results approach the experimental values with varying fidelity depending on the quantity. Yet, it is still possible to observe small discrepancies between the converged numerical results and the experiments. These discrepancies are caused by differences between the numerical and experimental apparatus, i.e., input errors. The experimental flow at the inlet is difficult to replicate because the spectral properties of the turbulence at the inlet were not experimentally measured (only the averaged quantities are available). The length and times scales of disturbances are crucial for transitional problems and particularly important for numerical predictions of the recirculation region length [98–101]. The differences in aspect ratio (length) of the cylinder [102–104], surface roughness [98,105,106], or Re [55] can also have a minor impact to the differences between numerical and experimental measurements. These results reaffirm the importance of the input error and also show the potential of

_{k}*r*and

_{i}*f*refinement studies to create reference solutions that can be used to analyze complex flows. Without all the information in Fig. 5, one could argue that an SRS model at $fk=0.40$ is “better” calculating $L\xafr$ than formulations at finer physical resolutions and draw modeling conclusions based on such idea. However, such a conclusion would be the result of insufficient data and not properly identifying error cancellation between input and modeling/numerical errors. In this case, the input error decreases $L\xafr$, while the modeling/numerical error increases it.

_{k}Figure 5 also shows that the selected mathematical model can accurately predict the present flow problem provided an adequate value of *f _{k}* and sufficient grid resolution. The results show only small discrepancies between experiments and simulations at $fk<0.50$. This is also seen in the values of $Ec(\varphi )$ and $Uv(\varphi )$ presented in Table 1 for three

*f*(see Refs. [60] and [65] for the comprehensive analysis). These quantities are obtained using the V&V20 metric [48]. For example, the results show that the values of $CL\u2032$ and $\u2212C\xafpb$ decrease from 0.664 and 1.41 at $fk=1.00$ (RANS) to 0.095 and 0.86 at $fk=0.25$, where the experimental values are equal to 0.096 and 0.86, respectively. This constitutes a reduction of $Ec(C\u2032L)$ from 591.8% to −1.1% and $Ec(C\xafpb)$ from 59.7% to −1.8%.

_{k}Model | $\Phi $ | $C\xafD$ | $CL\u2032$ | $\u2212C\xafpb$ | $L\xafr$ |
---|---|---|---|---|---|

$fk=1.00$ | $\varphi $ | 1.25 | 0.664 | 1.41 | 0.51 |

$Ec(\varphi )$ | +27.1 | +591.8 | +59.7 | −63.8 | |

$Uv(\varphi )$ | ±1.5 | ±57.6 | ±16.2 | ±49.7 | |

$fk=0.50$ | $\varphi $ | 1.04 | 0.284 | 1.05 | 1.12 |

$Ec(\varphi )$ | +5.7 | +195.7 | +19.3 | −25.6 | |

$Uv(\varphi )$ | ±32.2 | ±549.7 | ±61.7 | ±62.0 | |

$fk=0.25$ | $\varphi $ | 0.93 | 0.095 | 0.86 | 1.73 |

$Ec(\varphi )$ | −5.4 | −1.1 | −1.8 | +14.4 | |

$Uv(\varphi )$ | ±3.9 | ±55.9 | ±6.4 | ±20.2 | |

Expt. | $\varphi e$ | 0.98 | 0.096 | 0.88 | 1.51 |

U_{e} | ±0.01 | ±0.004 | ±0.01 | ±0.02 |

Model | $\Phi $ | $C\xafD$ | $CL\u2032$ | $\u2212C\xafpb$ | $L\xafr$ |
---|---|---|---|---|---|

$fk=1.00$ | $\varphi $ | 1.25 | 0.664 | 1.41 | 0.51 |

$Ec(\varphi )$ | +27.1 | +591.8 | +59.7 | −63.8 | |

$Uv(\varphi )$ | ±1.5 | ±57.6 | ±16.2 | ±49.7 | |

$fk=0.50$ | $\varphi $ | 1.04 | 0.284 | 1.05 | 1.12 |

$Ec(\varphi )$ | +5.7 | +195.7 | +19.3 | −25.6 | |

$Uv(\varphi )$ | ±32.2 | ±549.7 | ±61.7 | ±62.0 | |

$fk=0.25$ | $\varphi $ | 0.93 | 0.095 | 0.86 | 1.73 |

$Ec(\varphi )$ | −5.4 | −1.1 | −1.8 | +14.4 | |

$Uv(\varphi )$ | ±3.9 | ±55.9 | ±6.4 | ±20.2 | |

Expt. | $\varphi e$ | 0.98 | 0.096 | 0.88 | 1.51 |

U_{e} | ±0.01 | ±0.004 | ±0.01 | ±0.02 |

We need to emphasize that $CL\u2032$ is a quantity that does not permit direct comparisons of simulations at different physical resolutions (consistent with the first challenge involved in V&V of SRS identified in Sec. 1). Its magnitude will increase upon *f _{k}* refinement, as more of the fluctuations which contribute to $CL\u2032$ are resolved, and this behavior is independent of the modeling accuracy of the simulation. The present results show that the modeling error of simulations at $fk=1.00$ is larger than this effect since $CL\u2032$ grows as $fk\u21921.00$.

A very significant result of Fig. 5 from a V&V perspective is the minima observed on the solutions convergence with *f _{k}* on the two coarsest grids ($ri\u22651.35$). This behavior occurs at $fk\u22640.50$ and leads to significantly larger comparison errors at the finest physical resolution ($fk=0.15$) than at 0.25. It is caused by insufficient grid resolution to resolve all scales at $fk=0.15$ and demonstrates that the advantages and potential of SRS models can only be achieved with sufficient numerical resolution. The proposed V&V strategy identifies this behavior because it clearly distinguishes between the numerical errors and the modeling errors. Importantly, the method does not require reference data to make this distinction.

*f*). As observed in Table 1, the validation uncertainty is highest at $fk=0.50$. This can be explained in terms of the sensitivity of the flow to the effective Reynolds number, $Ree$,

_{k}*ν*is an effective viscosity, which can be decomposed as

_{e}with *ν*, *ν _{t}*, and

*ν*the molecular, turbulent, and numerical kinematic viscosities, respectively. Note that

_{n}*ν*comes from the numerical error. As demonstrated in Ref. [64], $Ree$ can determine the ability of a mathematical model to capture the key instabilities and coherent structures driving this flow. Since

_{n}*f*dictates

_{k}*ν*(and

_{t}*ν*indirectly), simulations at values of

_{n}*f*close to the threshold at which the model begins to resolve the key instabilities and coherent structures of the flow are highly sensitive to

_{k}*ν*. Alternatively, near a critical physical Reynolds number, the correct physical behavior is recovered only if both the turbulent and numerical viscosities are sufficiently small. Once again, this can be identified by the proposed V&V strategy. These results also represent crucial information for the utilization and development of SRS turbulence models, enabling a better understanding of the predicted flow physics and phenomena which are difficult to model. Finally, it is essential to mention that simulations at $fk=0.75$ can lead to slightly larger comparison errors than those at $fk=1.00$. This is a result of the calibration issues discussed in Sec. 1 that are not overcome by

_{n}*f*at coarse physical resolutions [11].

_{k}The trends of Fig. 5 and Table 1 are also observed in predictions of local flow quantities. This is shown in the following figures, which show the time-averaged pressure coefficient on the cylinder's surface, $C\xafp(\theta )$, and the streamwise velocity profiles, $V\xaf1$, in the near-wake at different grid resolutions for constant $fk=0.25$ (Fig. 6) and different values of *f _{k}* using the finest grid resolution (Fig. 7). All solutions converge upon the grid and physical resolution refinement and tend to the experimental measurement. Once again, the differences between simulations at different physical resolutions on the finest grid depend on the ability of the model to accurately resolve the instabilities driving the flow physics at a given

*f*. Such behavior is also observed for other SRS formulations [11,65].

_{k}### 4.2 Taylor–Green Vortex.

Next, we analyze the TGV results. This is a transient problem and, as such, V&V exercises are significantly more challenging and highly dependent on time. A flow misrepresentation at instant *t* will affect the subsequent instants and can compromise the overall simulation. This renders the accurate simulation and evaluation of the TGV (and RT) flow more complex than for the CC.

Figure 8 depicts the variation of the total (resolved and modeled) kinetic energy, *k*, with the grid, *r _{i}*, and physical,

*f*, resolution at instants before, during, and after the onset of turbulence (the solution at

_{k}*f*= 0 on $N=10243$ is not included). As discussed in Ref. [68], the period $t\u226412$ is crucial to the flow dynamics and simulation accuracy. At

_{k}*t*= 5, all simulations converge upon grid refinement ($ri\u21920$) and are independent of

*f*. This stems from the fact that the flow is laminar so that all simulations predict negligible magnitudes of turbulent stresses [68]. Later at

_{k}*t*= 7, when the flow is supposed to undergo transition to turbulence, the computations become dependent on

*f*and

_{k}*r*. Yet, it is observed that the solutions converge upon the refinement of these parameters, and the discrepancies between simulations at $fk\u22640.35$ are small. Also, it is possible to infer that varying physical and numerical resolution cause similar variations on the solutions of

_{i}*k*. At later times, the impact of

*f*relative to

_{k}*r*grows significantly, but the solutions still converge upon

_{i}*f*and

_{k}*r*refinement.

_{i}Figure 9 presents the temporal evolution of *k* and respective numerical uncertainty, $Un(k)$, at different physical resolutions, *f _{k}*. $Un(k)$ is estimated with the method proposed by Eça and Hoekstra [69] using the three finest grids available: $N=1283$, 256

^{3}, and 512

^{3}for all

*f*, and $N=2563$, 512

_{k}^{3}, and 1024

^{3}for $fk=0.00$. In Fig. 9, the black line indicates

*k*(

*t*), and the height of the red area the value of $Un(k)$ at a given instant. We assume the discretization error as the main contributor to the numerical error. It is important to reiterate that $Un(k)$ is the uncertainty associated with the numerical solution procedure, and so the red envelope in Fig. 9 is the range in which we expect the exact solution to the physical model to fall. The comparison error, $Ec(k)$, which is the difference between the numerical solution and the reference solution, can be large or small, independent of the numerical uncertainty. For example, for

*f*= 1, the error is large (Fig. 8), because the model performs poorly, but the numerical uncertainty at late time can be small, if the numerical solution is well converged. This means the numerical solution accurately represents the mathematical model, but that the model does not accurately capture the physics. Conversely, for

_{k}*f*= 0, at late time, the modeling error is small, but the numerical uncertainty is expected to be larger (for the same grid) because of the high requirements for both numerical and statistical convergence of a fully resolved turbulent flow. For a better physical and numerical interpretation of the results, recall that the onset of turbulence is expected to occur at $t\u22487$.

_{k}Referring to the simulations at $fk\u22650.50$, the results exhibit reduced $Un(k)$ for most of the simulation time. Taking the case at $fk=1.00$ (RANS), the estimates of $Un(k)$ do not exceed 3.3% of the predicted *k*. The exception lies in the period between *t* = 6 and 10, in which $Un(k)$ can reach 10.8%. Such values of $Un(k)$ are mostly caused by the difficulties of predicting the onset of turbulence with a RANS closure [68]. RANS models overpredict turbulence during transition, leading to a more rapid dissipation of *k*, which increases the flow gradients steepness and, consequently, the simulations' local numerical requisites. This illustrates the importance of a careful numerical and modeling interpretation of V&V exercises.

On the other hand, the refinement of *f _{k}* to 0.35 and 0.25 reduces the predictions' numerical uncertainty significantly. $Un(k)$ is negligible at early times and does not exceed 8.7% at late times when the flow features high-intensity turbulence and

*k*is significantly smaller than at

*t*= 0. However, the further refinement of physical resolution ($fk=0.0$) increases $Un(k)$ at

*t*> 11.0 considerably. This is the period when the flow is characterized by high-intensity turbulence and the largest numerical requirements since $fk=0.00$ aims to resolve all flow scales. For this reason, for simulations at $fk=0.00$, even the grid with $N=10243$ cells still has higher numerical uncertainties than those at $fk=0.25$ on $N=5123$. The simulation at $fk=0.00$ on the finest grid leads to $(Un(k))max=8.4%$ at

*t*= 13, whereas $Un(k)$ does not exceed 2.9% at this instant for $fk=0.25$.

*k*and its dissipation,

*ε*,

obtained at different *f _{k}* on the finest grid resolution available. The results of

*ε*are compared to the DNS of Brachet et al. [84] ($Ma=0,\u2009DNS1$) and Drikakis et al. [85] ($Mao=0.28,\u2009DNS2$). These reference studies do not provide numerical uncertainties. Figure 10(a) shows that all simulations converge monotonically upon

*f*refinement, and those at $fk\u22640.35$ exhibit small discrepancies, which do not exceed 0.0047. The exception in monotonic convergence is observed at

_{k}*t*> 11 at $fk\u22640.35$. This result is attributed to numerical uncertainty, especially for the case at $fk=0.00$ (Fig. 9). The results also show that all computations are independent of

*f*until $t\u2248tc=6$. This is the time when simulations at different

_{k}*f*lead to non-negligible levels of turbulence stresses.

_{k}The results for *ε* depicted in Fig. 10(b) exhibit similar tendencies. Until $t\u2248tc$, all computations are independent of *f _{k}* and lead to minor comparison errors. In contrast, the simulations become highly dependent on

*f*at $t>tc$, and their solutions converge upon

_{k}*f*refinement toward the reference data. Also, solutions of simulations at $fk\u22640.35$ are significantly less dependent on

_{k}*f*, and the resulting comparison errors are small. Comparing the two DNS data sets, our results converge toward those of Drikakis et al. [85] ($DNS2$). This stems from the fact that our simulations match the $Ma$ of such DNS, highlighting the role of input errors and $Ma$ [107] to the comparisons. Regarding simulations at $fk\u22650.50$, these lead to large comparison errors. For instance, the magnitude of the peak of dissipation predicted with $fk=1.00$ is 52.0% larger than that reported by Drikakis et al. [85].

_{k}In summary, Figs. 8–10 allow us to evaluate the simulations' modeling accuracy through *r _{i}* and

*f*refinement studies. These converge toward the reference solutions, leading to small comparison errors. The results also reaffirm that the proposed V&V method can be confidently applied when reference data are unavailable or the flow conditions are unknown since it can generate a reference solution ($ri\u21920$ and $fk\u21920$). The interpretation of the results also indicates the existence of two groups of simulations: those at $fk\u22650.50$ and $fk<0.50$. This can be explained by the critical value of

_{k}*f*that allows the SRS model to accurately resolve the key instabilities and coherent structures of the flow, not amenable to modeling. For TGV, these are the vortex reconnection phenomena responsible for the onset of turbulence [68].

_{k}### 4.3 Rayleigh–Taylor.

We conclude this section with a brief presentation of the RT flow, as these results exhibit behaviors similar to those of the CC and TGV. A complete discussion of the current results is given in Refs. [59] and [68]. The RT is a variable-density flow highly dependent on the initial conditions. Hence, we conduct simulations at different values of *f _{k}* and compare the results against the case at

*f*= 0 to guarantee the same initial conditions for all values of

_{k}*f*[108]. Figure 11 presents the temporal evolution of the mixing-layer height,

_{k}*h*, and the molecular mixing parameter,

*θ*, predicted with different values of

*f*on the finest grid resolution available for all

_{k}*f*($N=2562\xd7768$). The second quantity quantifies the homogeneity of the mixing-layer, with

_{k}*θ*= 0 and

*θ*= 1 referring to a heterogeneous and homogeneous mixture, respectively.

The results show that the simulations converge monotonically toward the reference solution upon physical resolution refinement, $fk\u21920$. Also, whereas simulations at $fk\u22640.35$ are in good agreement with the reference solution, those at $fk\u22650.50$ lead to large discrepancies. As discussed in Refs. [23] and [59], this is caused by the ability of simulations at $fk<0.50$ to accurately predict the instabilities and coherent structures governing the flow dynamics. This requires resolving the RT instability, laminar spikes and bubbles, and the Kelvin–Helmholtz instability responsible for the onset of turbulence in the RT flow. Once again, the grid and *f _{k}* refinement studies are critical to draw such conclusions.

Finally, Fig. 12 presents the numerical uncertainty of *h* obtained with *f _{k}* = 0 (most demanding case) using the three finest grid resolutions—$N=1282\xd7384,\u20092562\xd7768$, and $5122\xd71536$ cells. The solutions of

*h*for the three grids are also included. The results of Fig. 12(a) indicate that the numerical uncertainty of the simulations is relatively small and does not exceed 16.8% of the predicted value during the simulated time. Note that this represents an upper limit for $Un(h)$, since the coarsest grid used to estimate the numerical uncertainty does not possess sufficient resolution to resolve the interface perturbations at

*t*= 0. This idea is confirmed in Fig. 12(b), which compares the solutions of

*h*for the three grids. Whereas the two finest grids lead to nearly identical solutions, the coarsest one leads to large discrepancies. Such outcome demonstrates the importance of (adequately resolving) the initial flow conditions, and the adequacy of the grid with $2562\xd7768$ cells for this study.

## 5 Conclusions

This work investigated the importance of V&V to the credibility and further establishment of SRS of turbulence. A novel but simple V&V strategy based on grid and physical resolution studies is proposed that permits the separate assessment of numerical and modeling errors. It can be used when reference data are unavailable since it can estimate a reference solution from the physical resolution refinement studies. Also, it permits performing physical resolution refinements using the same initial flow conditions, avoiding the common difficulties generated by undercharacterized or unknown experimental conditions.

To ascertain the former aspects, we simulate three benchmark transitional flows using different grids and physical resolutions: the flow past a circular cylinder at $Re=3900$, the TGV at $Re=3000$, and the variable-density RT flow at $At=0.50$. The results confirm the importance of V&V to SRS of turbulence, indicate that it is possible to perform reliable V&V exercises with SRS formulations, and show the potential of the proposed V&V strategy. It enables an easier identification of cases experiencing error canceling, illustrates that insufficient grid resolution can suppress the modeling advantages of SRS of turbulence, permits a better flow analysis, and the identification of key flow phenomena. These properties can be used to develop better mathematical models and understand the reasons why models fail or succeed.

We expect the proposed V&V strategy to be a valuable contribution to further develop turbulence modeling, allow users attaining higher-quality computations, and promote the use of V&V in practical simulations of turbulence. Also, it can be a valuable tool for flows without available reference data or undercharacterized initial and flow experimental conditions. Naturally, performing grid and physical resolution refinement studies is computationally intensive, but the results of this study and Fig. 1 illustrate the risks of not assessing the accuracy of the simulations.

## Acknowledgment

We would like to thank the reviewers for their suggestions that improved our paper. Los Alamos National Laboratory (LANL) is operated by TRIAD National Security, LLC for the U.S. DOE NNSA.

## Funding Data

Laboratory Directed Research and Development program of Los Alamos National Laboratory (Project No. 20210289ER).

## References

*15th International Conference on Numerical Methods in Fluid Dynamics*(Lecture Notes in Physics, Vol. 490),

^{4}and 10

^{5}