Strategies for Computational Fluid Dynamics Validation Experiments

To date, turbulence modeling and Reynolds-averaged Navier–Stokes (RANS) equations represent the most sensible computational fluid dynamics (CFD) approach to predicting turbulent flows in practical engineering applications. High-fidelity CFD methods, such as direct numerical simulation (DNS) and large eddy simulation (LES), are forecast to remain impractical for many decades due to their prohibitive cost [2]. Even wall-modeled LES, although continuously gaining traction [3–5], is still too expensive for routine practical use. However, the known weaknesses of RANS models for predicting many classes of flows, notably involving surface roughness, pressure gradients, surface curvature, flow separation, and large-scale unsteadiness (see, e.g., [6] for a review), can be very problematic. These limitations and the industry's desire to increasingly utilize CFD for high-impact decisions [7] require improved turbulence models.

Further advancing turbulence models demands validation data from high-fidelity experiments and DNS designed to capture the physics and critical features of practically-relevant flows [8]. Beyond validating a model, these data enable the identification of model weaknesses, the calibration of model parameters, the training of data-driven models, and novel insight into the physics of turbulent flows. While its cost restricts DNS to simple flow cases, experiments, although still expensive, can access turbulent flows over complex geometries that challenge the predictive capability of current models.

However, designing and conducting CFD validation experiments demands utmost precision, attention to detail, financial resources, and adequate instrumentation to generate appropriate data. Unfortunately, many experimental results lack critical completeness, such as providing incomplete or entirely omitting uncertainty estimates, ignoring inflow non-uniformities, or assuming equivalence between as-designed and as-measured geometries [9]. These simplified assumptions propagate as unknown uncertainties directly impacting the boundary conditions used to setup corresponding CFD computations. As a result, conclusions regarding a model's predictive capability, accuracy, and weaknesses will be inaccurate, flawed, or in some cases, wrong.

Despite the turbulence modeling community being long aware of these shortcomings [10], only more recently have researchers actively attempted to mitigate them. The community has dedicated effort to establish standards and means for conducting, assessing, and documenting experiments specifically with model validation in mind [1]. These guidelines have created a considerable paradigm shift in the design and execution of fluid dynamics experiments. Indeed, model validation experiments adhering to these new standards or any flow experiment of such high fidelity have, to our knowledge, never been conducted before. However, several exceptionally well-conducted experiments focusing on model validation have paved the way into this new era of fluid dynamics experimentation, notably the experiments on NASA's wing juncture flow [11,12], NASA's two-dimensional (2D) hump [13,14], Sandia National Laboratories' revisitation of the Bachalo-Johnson hump [15], Byun and Simpson's bump [16,17], Boeing's speed bump [18,19], and Notre Dame's 2D ramp [20].

Compared to previous experiments, the design process of validation experiments is marked by conspicuous synergies between modelers and experimentalists and by increasingly stringent precision requirements. Performing high-fidelity model validation experiments is a new and developing practice, and the increasing demand for complex flow data is making this field proliferate, incentivizing researchers to take on the execution of such experiments. Consequently, to ensure the success of future model validation experiments and their conformity to the new quality standards, it is critical to establish adequate design procedures, risk mitigation strategies, uncertainty quantification (UQ) techniques, and procedures for the usage, documentation, and storage of validation data.

In this work, we present several strategies for the design, execution, and risk mitigation of model validation experiments and discuss their effectiveness in maximizing the outcomes and success of these experiments. Our arguments are presented based on the experiences collected and lessons learned during the Benchmark Validation Experiment for RANS/LES Investigations, or BeVERLI in short. Project BeVERLI is a research program, collaborative with NASA Langley Research Center, and part of the NATO Applied Vehicle Technology program, AVT-349—Non-Equilibrium Turbulent Boundary Layers in High Reynolds Number Flow at Incompressible Conditions². This project aims to produce detailed experimental data sets of smooth-wall flows with various levels of separation that, for the first time, meet the highest levels of model validation experiment completeness (MVEC) per Oberkampf and Smith [1]. Thereby, incompressible, high Reynolds number flow of a plane turbulent boundary layer over a three-dimensional (3D) hill, termed the BeVERLI Hill, is investigated. The project is an international collaboration that brings government labs, industry, and academia together to work on the validation experiment and corresponding simulations. The main validation experiment is carried out at the Virginia Tech (VT) Stability Wind Tunnel, and as of September 2022, three wind tunnel entries conducted in February 2020, May 2021, and October 2021 have been completed. The VT team for these entries consisted of six members, three for experiments and three for computations, with one member contributing to both. They had shared access to archived and version-controlled data, using log files for tracking changes. While currently available upon request to the present authors, the data will be publicly accessible in the future. Additional experiments with scaled versions of the same hill geometry are also being conducted at the University of Toronto Institute for Aerospace Studies (UTIAS) and SINTEF Ocean in Norway [21]. Similarly, several groups have computed the main BeVERLI Hill geometry for several flow conditions, numerical solvers, grids, and turbulence models using RANS [22].

In what follows, we outline the experimental and computational setup for the primary validation experiment and simulations conducted at VT, including approaches to selecting a suitable flow case for CFD model validation and an overview of the achieved MVEC level as of September 2022. We then present and discuss the proposed strategies for model validation experiments, and conclude with a summary and thoughts on future work.

A flow and its model geometry should exhibit particular qualities to be suited for CFD validation. Selecting a model geometry entails nuanced reasoning, branching into different aspects of a validation case. We suggest reducing the discussion to two attributes with primary and secondary evaluation criteria, as listed in Table 1.

First, the type of CFD models and corresponding weaknesses a validation study should address are considered. Clearly, if a CFD model's predictive performance is lacking, e.g., for separating flows over curved surfaces, the validation geometry should involve curvature and generate appreciable adverse pressure gradients. In contrast, considering the target CFD model type is less obvious. Different flows are more or less suited for simulation using the different CFD methods (RANS, LES, DNS, and hybrid methods) and, consequently, for validating a specific CFD model. However, the features making a flow suitable are, at the same time, simple and deceiving and should be evaluated with care. For example, in (steady) RANS model validation, the uniqueness and stability of a flow field's long-time or ensemble average are typically desired qualities. But, even in steady validation, it is crucial not to disregard flow cases that do not exhibit these qualities. Flows exist and need to be modeled in real-life applications that do not obey these restrictions, and while such flows might not be entirely suited for RANS, we must consider that their modeling accuracy is currently unknown. Moreover, clarifying a model's predictive capability, even against inapt flow characteristics, leads to improved models and solution strategies.

The second attribute considers a flow's complexity and realizability as an experiment. Related to the flow's realizability is the accessibility to measurements of flow quantities of interest (QoIs) over the model geometry. The selected geometry should be simple enough to be manufactured, placed in a test facility, and accessed with measurement instrumentation. It is beneficial if the model allows the control of experimental parameters to access various flow regimes, e.g., different degrees of flow separation, representing distinct validation flow cases. In addition, having the ability to integrate instrumentation on the surface of the model, e.g., optical access for measurement techniques like particle image velocimetry (PIV) and laser Doppler velocimetry (LDV), static pressure probes, or surface microphones, is a plus.

BeVERLI was conceived to capture aerodynamic flow with the same range of physics present in high-lift wings, typically subjecting the flow to separation and non-equilibrium effects under the simultaneous action of pressure gradients and curvature. BeVERLI further focuses upon the study of incompressible, high Reynolds number, smooth-wall turbulent boundary layers (TBLs) and steady RANS/LES model validation. RANS models are considered adequate for many classes of flows, such as fully attached equilibrium TBLs [23,24]. They are less adequate or even unacceptable for flows with massive separation [23,25,26] and non-equilibrium flows with streamline curvature and strong pressure gradients [6]. A 3D hill geometry was selected for project BeVERLI to study the development of an initially 2D high Reynolds number TBL as it departs from equilibrium over the hill, experiencing relatively strong pressure gradients and curvature, which ultimately leads to separation on the hill's leeward side. The simplicity of the hill geometry allows further alteration of these parameters, e.g., by yawing the hill about its vertical axis to change its position relative to the undisturbed oncoming flow. Two yaw angle orientations, later referred to as the 0 deg and the 45 deg yaw cases, are intended to generate massive and moderate flow separation, respectively. The hill's aspect ratio and the ratio of the approaching TBL's thickness to the hill's height were similarly adjusted to generate the desired level of flow separation. Other geometrical features of the BeVERLI Hill, such as its cylindrical or flat top region, as detailed later, were designed to integrate optical access windows for PIV and LDV.

As explained later, for certain hill-height-to-approach-boundary-layer-thickness ratios and height-based Reynolds numbers, the 0 deg and 45 deg yaw cases exhibit reflectional symmetry breaking and bi-modality in both the simulations and experiments. That is, the long-time-averaged flow field over the BeVERLI Hill takes on distinct spanwise asymmetric modes. The spanwise symmetric mean flow, expected in an experiment with symmetric boundary conditions, cannot be sustained. Instead, an asymmetric, presumably more stable, mean flow field is observed, featuring a spanwise shifted wake on the leeward side of the hill. We hypothesize that the stability of the flow is susceptible to small system asymmetries as found, e.g., in as-tested experimental boundary conditions or CFD solver algorithms. Asymmetries may also originate within the flow as shear layer instabilities from the interaction of the distinctly developing boundary layer over the top and flanks of the hill. The asymmetric wake is steady for certain scale ratios of the 45 deg yaw case, but an unsteady, bi-modal wake is found that chaotically switches sides spanwise for the 0 deg yaw case. Naturally, given its peculiarities, it would be tempting to disregard the symmetric orientation BeVERLI cases entirely for RANS validation. However, while the unsteady wake bi-modality of the 0 deg yaw case might be unsuitable for RANS, it is important to recognize that the steady asymmetry found in the 45 deg yaw case is a desired characteristic. Symmetry breaking occurs in many relevant engineering applications, e.g., for the flow around an Ahmed body [27–30], and represents a challenge for current turbulence models. Disregarding the symmetric orientation cases solely due to asymmetry, without considering all aspects of this symmetry breaking feature, would be a missed opportunity. Moreover, the 0 deg yaw case may represent a valid candidate for unsteady RANS (URANS) validation.

The BeVERLI Hill geometry is shown in Fig. 1. The 3D hill has a design length, w, of 0.9347 m and a height, H, of 0.1869 m, resulting in an aspect ratio, w/H, of 5. To construct the geometry, a Cartesian coordinate system, x_i, is used. The hill base is in the x₁–x₃ plane, and its height is erected in the x₂-direction. The 2D centerline (⁠$x3=0$ plane) and centerspan (⁠$x1=0$ plane) cross section of the 3D hill geometry are constructed by mirroring a fifth-degree polynomial, $P5(x1)$ and $P5(x3)$⁠, about the $x1=0$ and $x3=0$ plane, respectively, and joining the original and mirrored polynomial by a flat curve through the center. The fifth-degree polynomial is precisely defined by Gargiulo et al. [31] and provides enough parameters to ensure zero slope and curvature at the base of the hill and the flat top. The polynomial and flat curves are represented as dashed lines in Fig. 1. The resulting cruciform is symmetrically extruded about the centerline and centerspan to form a squared flat top and cylindrical surfaces of width $s=$ 0.0935 m, surrounding the dashed lines of the cross-sectional profiles as shaded surfaces in Fig. 1. The hill corners, labeled C₁ through C₄ (with C₂ shaded and outlined by dotted lines in Fig. 1), are created using lofted surfaces. These surfaces are constrained by the polynomial edges of the perpendicular cylindrical surfaces adjacent to each corner and a segment of a parametric fourth-order superellipse at the hill's base, the latter which is also precisely defined by Gargiulo et al. [31]. These three constraining curves are delineated for corner C₂ in Fig. 1, using dotted lines. To accurately define each corner's loft surface, a 3D computer model of the BeVERLI Hill was created using Solidworks 3D CAD software. The CAD model of the BeVERLI Hill, reflecting the precise definition of each corner, can be downloaded from the link reported in the footnotes³.

Two BeVERLI Hill models were CNC-milled out of epoxy tooling foam as shown in Fig. 2. The first features 135 pressure probes installed in the model's surface to measure mean and unsteady surface static pressure. The second one features slots with anti-reflective, clear windows for LDV measurements along the fifth-degree polynomial curves. A circular base 1.1684 m in diameter was added to both models to facilitate their installation and rotation in the VT SWT. Both models were painted using a black base coat and a glossy clear coat to support optical measurements, including PIV, oil film interferometry (OFI), and oil flow visualization. Similarly, two of the four dedicated LDV windows were painted to allow overlapping OFI at LDV measurement locations. The models were scanned after manufacturing, but before paint was applied to the surface with a GOM ATOS 5 blue light laser scanner with a resolution of 8 megapixels. These scans found that the pressure-tapped and the LDV model matched the as-designed geometry to within ±0.40 mm and ±0.30 mm, respectively.

Important to note is that a different initial model of the BeVERLI Hill (not shown here), made entirely from optically clear acrylic via vacuum forming, was used for the first wind tunnel entry in February 2020. This model aimed to provide unrestricted optical access for PIV and LDV measurements. However, as explained later, it deviated significantly from the intended BeVERLI geometry and presented other notable challenges. Consequently, it was replaced by the two precise CNC-milled models for subsequent tests.

The primary BeVERLI Hill experiment is performed at the Virginia Tech Stability Wind Tunnel (VT SWT), a continuous, single return, subsonic wind tunnel, with a test section that is nominally a rectangular cuboid 1.85 m by 1.85 m in cross section and 7.32 m long. A schematic of the facility is shown in Fig. 3. Flow speeds up to 80 m/s at turbulence intensities as low as 0.01%–0.03% can be achieved. The facility features an adaptable hybrid test section with two interchangeable walls: a smooth-wall aerodynamic configuration and an anechoic configuration with Kevlar windows as side walls. The smooth, hard wall configuration was used for the present study [32]. Lastly, the contraction upstream of the test section features an instrumented liner to track the static pressure along its surface [33] and a 3.18 mm boundary layer trip.

Surface and field measurements were performed over the BeVERLI Hill in a configuration as illustrated in Fig. 4. The hill was mounted at the center of the test section's port-side wall. A Cartesian coordinate system, x_i, whose origin is located at the center of the BeVERLI Hill in the plane of the tunnel wall, is prescribed to locate points within the configuration. The x₁-axis of the coordinate system points in the freestream direction, the x₂-axis is in the tunnel wall-normal direction, and the x₃-axis is spanwise and completes the coordinate system in the right-handed sense. The BeVERLI Hill is located 6.88 m downstream of the 3.18 mm boundary layer trip in the wind tunnel contraction, and 3.40 m downstream of the end of the wind tunnel contraction, coinciding with the start of a transition structure that connects the contraction to the test section.

The BeVERLI Hill can freely rotate about its x₂-axis, allowing different orientations of the model with respect to the oncoming flow in the VT SWT. The main hill orientation considered in the present study is as shown in Fig. 4, with one of the pairs of rounded corners directly aligned with the approaching freestream. This orientation is referred to as the 45 deg yaw case. However, a 0 deg yaw case, with the cylindrical surfaces of the hill directly aligned with the approaching freestream, was additionally measured. Measurements were made with the freestream flow adjusted to achieve the three hill-height-based Reynolds numbers, $ReH=H·U∞/ν$⁠, of 250,000, 325,000, and 650,000. However, the present study focuses on the 250,000 and 650,000 cases.

Additionally, extensive measurements of the boundary conditions of the facility have been conducted, including detailed measurements of the oncoming boundary layer and its cross-sectional uniformity. The corresponding thickness, $δ0.99$⁠, displacement thickness, $δ*$⁠, and momentum thickness, θ, of the oncoming boundary layer at $(x1, x2, x3)$ = $(−1.83 m, 0, 0)$⁠, or $(x1/H, x2/H,x3/H)$ = $(−9.78, 0, 0)$⁠, as measured using a 32-probe boundary layer rake, are summarized in Table 2 for $ReH≈$ 250,000, 325,000, and 650,000.

Mean static pressure, $P¯$⁠, and instantaneous pressure, P, measurements on the surface of the BeVERLI Hill were collected using 1.6 mm (0.063 in) diameter steel tubes embedded in the model, which were connected via Tygon tubing to scanners for the acquisition. Exactly 95 of the 135 total pressure taps were connected to three DTC ESP 32HD scanners with a range of 2.5 psi and a manufacturer-rated accuracy of 0.06% full-scale. The remaining taps were connected to an Esterline 9816/98RK pressure scanner with a rated accuracy of 0.05% full-scale. A layout of the pressure taps distribution over a RANS contour plot of the mean pressure coefficient, $Cp=(P¯−Pref)/(12ρrefUref2)$⁠, on the 45 deg yaw case hill surface at $ReH≈250,000$ is shown in Fig. 5. The tap layout features a denser agglomeration of taps along the centerline, centerspan, and diagonals of the hill. One half of the hill has a larger number of taps to enable the reconstruction of mean pressure over the entire surface of the hill by interpolating the data from symmetrically equivalent hill orientations, e.g., 45 deg, 135 deg, 225 deg, and 315 deg.

The mean velocity, $U¯i$⁠, and turbulence statistics of the boundary layer over the BeVERLI Hill were measured using a boundary layer rake, PIV, and LDV. The corresponding measurement locations are presented in Fig. 6 over a RANS C_p contour plot and surface streamlines of the 45 deg yaw case hill at $ReH≈250,000$⁠. Notably, the oncoming boundary layer is measured using all three techniques at (⁠$x1/H,x3/H$⁠) = (⁠$−9.78,0$⁠) (diamond symbol in Fig. 6). Similarly, at (⁠$x1/H,x3/H$⁠) = (⁠$−0.8,−0.8$⁠) (star symbol in Fig. 6), PIV and LDV are employed simultaneously.

Particle image velocimetry data were acquired using a time-resolved, two-dimensions-three-components (2D3C) system. The data were measured along specified x₁–x₂ flow planes except at (⁠$x1/H,x3/H$⁠) = (⁠$−0.8,−0.8$⁠), where PIV and LDV are employed simultaneously and the PIV plane is aligned diagonally with the LDV windows. The PIV planes were distributed to roughly capture the passage of three streamlines with distinct straining histories, which are highlighted as dashed curves in Fig. 6. The PIV system's components and configuration are shown in Fig. 7 with corresponding labels (1)–(9). The system consists of a dual-head, high-repetition-rate frequency-doubled Nd:YAG laser (2) of 532 nm wavelength from Photonics Industries International, Inc., allowing repetition rates up to 35 kHz. The laser was securely mounted under the test section of the VT SWT on a specifically designed optical table. A precisely calibrated articulating optical arm (5) from the company LaVision directed the laser from under the test section into it at the desired measurement location through an optically clear acrylic window mounted on the tunnel side-wall directly opposite of the BeVERLI Hill model. An assembly of collimating optics and a 50 mm cylindrical lens attached at the end of the laser delivery arm (7) generated the required PIV laser sheet with an approximate thickness of 1.5 mm. Two Phantom v2512 high-speed cameras (9) acquired two sets of 24,000 double-frame images at a frame rate of 1 kHz and 12.5 kHz. The low-rate, long-time PIV data (slow set) is utilized to achieve better convergence of statistical quantities, whereas the high-rate, short-time PIV data (fast set) provides time-resolved information. The cameras were mounted under the test section floor of the VT SWT and imaged the measurement location through additional floor-mounted, optically clear acrylic windows. Depending on the working distance, a Nikon camera lens of 300 mm fixed focal length and a Sigma lens of 105 mm fixed focal length with a Nikon 2.0x teleconverter were used. All the equipment were connected through a LaVision PTUx (programmable timing unit) (4) to a computer (6) running the LaVision processing software DaVis. The flow in the VT SWT was seeded using an MDG MAX 3000 APS fog machine and MDG neutral fluid, providing seeding particles roughly 0.5–0.7 μm in size, which resulted in Stokes numbers small enough to capture the full range of relevant scales. In order to calculate the particle displacements needed for generating the PIV vector field using the DaVis software, the frames of each double-frame image underwent a multipass cross-correlation process. This process involved two initial passes with a 64 × 64 pixels interrogation window, followed by two final passes with either a 16 × 16 pixels window or, depending on the image quality of a specific measurement plane, a 24 × 24 pixels window. The final passes were executed with 75% window overlap. The resulting PIV spatial resolution ranged from approximately 1.8–2.5 mm, with the closest data points obtained at a distance of 0.3–1.5 mm from the surface of the hill.

Laser Doppler velocimetry was employed to acquire near-wall point-wise measurements of the boundary layer along the dedicated windows on the BeVERLI Hill model. An embedded three-component system was utilized with the same flow seeding method of PIV. The system's probe was mounted inside the BeVERLI Hill. Five corresponding measuring laser beams passed through dedicated 1.5 mm anti-reflective curvature-fitting acrylic windows mounted on the surface of the hill and crossed within a measurement volume of approximately 50 μm in diameter at the measurement location of interest over the hill surface. This allowed measurements up to a minimum and maximum height above the hill surface of 0.1 mm and 30 mm, respectively. The LDV system is described in its entirety by Duetsch-Patel et al. [34].

The modular total pressure boundary layer rake consists of 30 Pitot probes aligned in the direction of the flow. The probes are arranged with a linear spacing of ten probes that are flush against each other closest to the wall and twenty logarithmically spaced probes further from the wall. The rake is encased in a 3D printed (PA 614—40% glass-filled Nylon 12) NACA 0012 fairing, with a chord of $0.1143 m$⁠. The probes extend $0.1016 m$ from the leading edge of the fairing and connect to DTC ESP 32HD scanners with a range of 20 in H₂O for the acquisition.

Direct surface shear stress measurements were acquired using OFI. The data were collected primarily at the LDV locations and upstream of the BeVERLI Hill. Two low pressure sodium vapor lamps emitting monochromatic light at 589 nm wavelength were used to generate the required light interference fringe pattern between silicon oil of different viscosities and the hill surface. The hill surface paint finish was selected for optimal light reflectivity while suppressing secondary rays. The fringes were imaged using a Sony Alpha a6500 mirrorless digital camera and processed using the matlab code developed by Naughton [35] and Liu [36] to obtain the desired shear stress information on the surface of the BeVERLI Hill. Details concerning the OFI setup, system calibration, instrumentation parameters, and the processing of the fringe images are detailed by Sundarraj et al. [37].

In addition to the above quantitative measurement techniques, surface oil flow visualization was used to acquire a broad, qualitative picture of the flow over the BeVERLI Hill. The hill areas of interest were lightly painted with oil in streamwise streaks, making multiple passes over each segment of the model. Two different oil mixtures, one of which was fluorescent, were used for the measurements. The non-fluorescent mixture consisted of kerosene, titanium oxide, and oleic acid (10:4:1 volume ratio), whereas the fluorescent mixture consisted of AeroShell 100 mineral oil aviation oil and kerosene (2:1 volume ratio). The VT SWT was then turned on such that the flow moved the oil over the hill. The tunnel was run until motion of the oil was no longer visible and the flow had reached a steady-state. Both pictures and videos were taken using a Sony Alpha a6500 mirrorless digital camera and a 18–135 mm zoom lens mounted behind the tunnel wall opposite to the model.

For the RANS simulations discussed in the present study, the commercially available ANSYSFluent solver [38] was used. The density-based solver with algebraic multigrid method was selected to solve the RANS equations. The spatial discretization of the inviscid fluxes is performed using a second-order upwind scheme. Fluent achieves second-order accuracy using a multidimensional linear reconstruction approach [39]. The viscous and turbulent fluxes in Fluent are discretized to second-order accuracy using a central flux scheme and the Green-Gauss theorem.

The computational domain for modeling the BeVERLI Hill in the VT SWT, as depicted in Fig. 8, is limited to the test section and excludes other tunnel circuit components. However, the inflow plane of the computational domain is extended to provide a similar oncoming boundary layer as experimentally measured. Similarity is assessed based on integral boundary layer parameters, compared between simulation and experimental results in Table 3. Notably, the TBL thickness is determined based on 95% of $U∞$ to account for numerical errors and improve comparison reliability. The indicated uncertainties correspond to 95% confidence intervals for the experiment and encompass the scatter caused by the choice of turbulence model and numerical grid level in the simulation. Additionally, the outflow plane location is placed outside of the domain of influence of the hill.

The boundary conditions are as follows: a Dirichlet condition specifying the stagnation temperature, $T0,in$⁠, and the stagnation pressure, $P0,in$⁠, is applied uniformly across the inlet plane. The stagnation conditions are obtained from the experiment. At the outlet, a Dirichlet condition fixes the static pressure, $Pout$⁠, which is manually adjusted to match a measured reference pressure, $Pref$⁠, at the specified reference location of the experiment. The additional domain boundaries are modeled as smooth, adiabatic, viscous, and non-permeable no-slip walls. The flow of air within the system is modeled to be ideal, compressible, and viscous. The absolute dynamic viscosity, μ, is defined in accordance with Sutherland's law [40]. Free-stream turbulence quantities are determined using the turbulence intensity, I, obtained from calibration data of the VT SWT via a hot wire probe, along with an estimated turbulent viscosity ratio, $μt/μ$⁠. The boundary and important reference conditions for $ReH≈$ 250,000 and 650,000 are summarized in Table 4.

Reference values for the computation of normalized quantities, such as C_p or the friction coefficient, $Cf=τw/(12ρrefUref2)$⁠, are obtained at a specific common reference location for both the experiment and simulations. The reference pressure to match between experiment and CFD is calculated by averaging the pressure measurement from seven spanwise-distributed surface static pressure ports located $2.23 m$ upstream of the hill center on the wall opposite of the hill model. The exact Cartesian coordinates of the reference ports are summarized in Table 5. The reference probes are labeled and depicted as small spheres in Fig. 8. Other reference quantities, e.g., $ρref$ and $Uref$⁠, are computed from the reference pressure and the stagnation conditions using standard compressible isentropic flow relations.

Two turbulence models were considered for this study: the standard Spalart–Allmaras (SA) one-equation eddy-viscosity model [41] and the Menter k–ω shear stress transport (SST) two-equation eddy-viscosity model [42]. The standard SA model was chosen for the BeVERLI Hill's curved geometry over the version with rotation/curvature correction to facilitate comparison with solutions from different solvers used by various groups working on the case, including non-commercial academic development codes. This decision ensured consistent comparisons due to the wider availability of the standard SA model.

A systematic family of four structured curvilinear grid levels was generated. The commercial software Pointwise [43] was used for meshing. The grids feature a denser, high-resolving allocation of nodes in the regions of the BeVERLI Hill and the tunnel walls to fully resolve the boundary layer and ensure a normalized wall distance of $x2+=uτx2/ν≤1$ on all grid levels. Furthermore, the grids were constructed on one-half of the computational domain and mirrored about $x3=0$ to ensure spanwise symmetry. The grid levels are identified from finest to coarsest as Levels 1 through 4 and are, respectively, $557×353×385$ (Level 1), $443×279×305$ (Level 2), $351×221×241$ (Level 3), and $279×177×193$ (Level 4) nodes in size. Each level is generated starting from the finest Level 1 grid and then systematically coarsened using a non-integer refinement factor of $r=21/3$ to obtain each subsequent grid level down to the coarsest Level 4 grid. An analysis of the geometric similarity of these grids is given in Appendix  B, showing they were correctly systematically refined. The iterative convergence criterion for both the mean flow and turbulent equations was set to achieve a relative iterative residual reduction of at least 8 orders of magnitude from the initial levels. All equations were run using first-order-accurate schemes for the first 3000 solver iterations before switching to second-order-accurate schemes. Using this protocol, the simulations typically converged between 22,000 and 28,000 solver iterations for all grid levels. The simulations were run on a high-performance computing cluster featuring 308 nodes, each running two AMD EPYC 7702 chips (128 cores/node), to take advantage of code parallelization. For the solution to reach convergence within 12–24 hours, 500 cores were used for the Level 1 mesh, 300 for Level 2, 180 for Level 3, and 120 for Level 4.

The grid dependence of the simulations was analyzed to assess the numerical uncertainty caused by discretization error (DE). The detailed DE analysis is provided in Appendix  C of the paper. The analysis revealed that the computed flow quantities of interest, such as the hill surface C_p and C_f distributions, exhibit sufficiently low DE on both the Level 2 and Level 1 meshes for the purpose of this study. Additionally, the DE analysis showed that the SA model yields lower uncertainty due to DE compared to the SST model at $ReH≈250,000$⁠, but higher uncertainty at $ReH≈650,000$⁠. The results of the DE analysis are used in this paper to incorporate error bars in the presented computations whenever applicable.

An internal assessment of BeVERLI's MVEC as of September 2022 was conducted, following Oberkampf and Smith's guidelines [1]. While the full assessment presentation is beyond this paper's scope, Appendix  A offers a summary and a breakdown of the MVEC score for each corresponding attribute suggested by Oberkampf and Smith [1]. The results highlight the effectiveness of the strategies discussed in the subsequent sections. The majority of attributes achieved the highest completeness level 3, with the remaining categories scoring between levels 2 and 3, predominantly leaning toward 3. A comprehensive NASA Technical Memorandum (TM) is currently underway, which, among other relevant items, will include the detailed MVEC assessment and plans to achieve completeness level 3 in all categories. An important focus is placed on improving documentation for the experimental facility and boundary conditions, aligning with the development of “Common Research Wind Tunnels for CFD Verification and Validation” under the NATO AVT-387 panel [44]. Additionally, the prospect to conduct an MVEC assessment through external reviewers is considered. The TM's publication is anticipated for the end of 2023.

Synergistically supporting the experiments of a validation study with corresponding CFD simulations, whether at a reduced scale for risk-reduction or the actual target scale, is vital. The purpose of CFD simulations in a validation study is twofold. First, parametric studies can be performed before the validation experiment to investigate the flow under a vast and diverse range of conditions that would otherwise not be accessible or which are too expensive to measure experimentally. Key flow features, many of which go unseen or are unresolved in the experiments, and their associated sensitivities can be studied over a selected geometry while altering its parameters or other boundary conditions. As a result, the flow's suitability for the validation study can be investigated early, and the model's geometry and experimental boundary conditions can be optimized, adapted, or entirely changed, avoiding opportunity costs, e.g., arising from the manufacturing of a suboptimal model geometry. Additionally, measurement locations and QoIs can be chosen methodically and coordinated efficiently across various measurements, allowing the design of a comprehensive test matrix with repeated and redundant measurements that enables extensive data UQ. Similarly, information on the required scaling and range of instrumentation and sensors can be extracted to adjust, expand, or change the available measurement equipment. The computed flow fields can even be utilized to estimate aerodynamic loads to mitigate challenges associated with installing model interfacing structures and instrumentation within a testing facility.

Second, while the CFD results cannot be considered the ground truth before rigorous validation, observing qualitatively similar trends between CFD and experiments increases the confidence and likelihood that the measurements are of good quality. Conversely, the CFD framework and its boundary conditions can be qualitatively validated, resulting in adjustments to the framework, such as improved computational grids and boundary conditions that more accurately represent the validation experiment.

Reynolds-averaged Navier–Stokes simulations were a determining tool in project BeVERLI. As alluded to previously, the flow over the BeVERLI Hill displays symmetry breaking at symmetric model orientations. The long-time averaged flow field has multiple solutions for the 0 deg and 45 deg yaw cases, a hypothesized metastable symmetric and a stable asymmetric mode. This characteristic would have gone unnoticed for a long time without the combined use of experiments and simulations, leading to inadequate preparation and failure of the experiment. The symmetry breaking characteristics of the BeVERLI Hill manifest as follows. In the RANS simulations, despite perfectly symmetric boundary conditions and meshes, the 0 deg and 45 deg yaw cases converge to a series of distinct spanwise asymmetric solutions instead of the intuitively expected symmetric flow. The asymmetry is visible in Fig. 9, showing C_p contour plots with overlaid surface shear stress (C_f) lines from the 45 deg yaw case RANS solutions on the Level 1 grid at $ReH≈$ 250,000 and 650,000. The asymmetry most prominently affects the wake region of the hill. As visible in Fig. 9, the flow is sensitive to the choice of turbulence model and the Reynolds number. Additionally, as reported by Gargiulo et al. [22], the flow exhibits sensitivity to the numerical solver, the computational grids, and the iterative solution convergence. In the experiments, the asymmetry is present for both the 0 deg and 45 deg yaw cases. However, for the 45 deg yaw case it is only observed at $ReH≈$ 650,000. More importantly, as will be detailed later, the long-time behavior of the asymmetry for the 0 deg case is considerably different than for the 45 deg yaw case. Figure 10 shows images from oil flow visualization over the surface of the 0 deg and 45 deg yaw case hill at $ReH≈650,000$⁠, clearly indicating the presence of the wake asymmetry in the experiments.

Regarding the nature of the symmetry breaking, we hypothesize that the intuitively expected symmetric flow over the BeVERLI Hill is metastable and susceptible to small system asymmetries, e.g., due to asymmetries in the experimental boundary conditions or as present in particular CFD solver algorithms. Algorithmic asymmetries can occur, for example, if a numerical scheme has a preferred sweep direction through a computational grid. Alternatively, the asymmetries may be inherent to the flow and manifest themselves as shear layer instabilities from the interaction of the distinctly developing boundary layers over the top and sides of the hill meeting on the leeward side. Once triggered by asymmetries, the flow departs from its symmetric state and settles into a stable and more energetically efficient asymmetric mode. While the asymmetry itself does not break the integrity of the validation experiment and represents an exciting feature that occurs in engineering problems of practical relevance and that challenges current turbulence models, its stability has important implications. An intermittent switching of the wake asymmetry between different modes favoring distinct sides of the hill would be undesirable for (steady) RANS model validation. Indeed, if the time scale at which the switching events occur is smaller than the time scale of the experiment, a stable mean flow field cannot be measured, jeopardizing the data's repeatability and consistency. Unsteady surface static pressure measurements, detailed in a later section of the present paper, indicate that the wake asymmetry is stable for the 45 deg yaw case. In contrast, an unsteady, bi-modal wake is found that chaotically switches sides in the spanwise direction for the 0 deg yaw case. Thus, the 45 deg yaw case represents a perfectly valid RANS model validation case despite the asymmetry, whereas the 0 deg yaw case is inadequate for RANS, both because of the challenges involved with measuring the mean flow field of this case and because of the fundamental limitations of steady RANS models. However, the 0 deg case remains valuable for development and validation of modeling in other simulation approaches, such as URANS or wall-modeled LES.

Recognizing the symmetry breaking nature of the BeVERLI Hill occurred in steps and required the combined use of simulations and experiments. The success of BeVERLI is, among other things, attributable to this critical discovery, leading to adequate experimental preparation, appropriate experimental design, and substantial improvements in both the experimental setup and the computations. Early risk-reduction experiments by Gargiulo et al. [31] revealed a relatively symmetric flow over the BeVERLI Hill. These experiments examined the 0 deg and 45 deg yaw cases, with $ReH$ up to 325,000. Note, at $ReH≈$ 325,000, tested at both reduced and full scale, the hill-height-to-TBL-thickness ratio, $H/δ$⁠, varies due to the difference in model and testing facility scale, yielding a value of 1.5 for the reduced- and 3.2 for the full-scale case. The 0 deg yaw case exhibited symmetry breaking at the full-scale but not at the reduced scale for $ReH≈$ 325,000, indicating sensitivity to $H/δ$⁠. Slight flow asymmetries in the risk-reduction experiments, while present, were attributed to the testing facility. Although some PIV data were available during the risk-reduction phase, the asymmetry analysis primarily relied on (steady) static surface pressure and oil flow visualization data, the latter which are shown in Fig. 11. Nevertheless, isolating the critical asymmetry feature ultimately required adding RANS simulations.

Because the asymmetric solutions were unexpected in the simulations, the CFD framework was checked for any anomaly in the nominal boundary conditions, particularly concerning the symmetry of the numerical grids. No appreciable asymmetry was identified, but this discovery incentivized a remeshing of the numerical grids, building them on one-half of the computational domain in the spanwise direction and mirroring them about the $x3=0$ plane to warrant and preserve symmetry. The subsequent continued presence of asymmetry in the simulations strengthened the hypothesis of symmetry breaking but could not exclude a numerical artifact. The simulations became crucial in helping correctly reprioritize the goals of the subsequent full-scale validation experiment. Experimentally capturing and investigating the presence of asymmetry became a primary target. The successful observation of asymmetry in the experiment, was a significant achievement and effectively confirmed the hypothesis of symmetry breaking. Subsequent in-depth boundary condition perturbation studies and, as detailed later, unsteady static surface pressure measurements and cross-institutional experiments and simulations, revealed crucial sensitivities and the bi-modality of the asymmetry for the 0 deg yaw case, leading to the shift of focus from the 0 deg to the 45 deg yaw case.

Furthermore, while the flow asymmetry has guided many of the design choices of project BeVERLI, other critical risks were mitigated through the combined use of experiments and simulations. While letting the asymmetry go unnoticed, the early risk-reduction experiments drew attention to a number of flow features, including a Reynolds-number-dependent separation region on the windward side of the hill close to the hill top for the 0 deg yaw case (labeled “S” in Fig. 11(a)). This phenomenon is due to Reynolds-number-dependent relaminarization due to the rapid acceleration of flow on the windward face. At low Reynolds number, the acceleration parameter becomes large enough to dramatically reduce turbulence in the boundary layer, which makes the flow susceptible to separation in the adverse pressure gradient at the top of the hill. Turbulence model predictions in accelerating flows have been identified by AVT-349 as a key area in need of further improvement for vehicle aerodynamics predictions with RANS. The windward separation region had the potential risk of drastically changing the downstream flow characteristics and obstructing flow features that the hill geometry was meant to produce naturally. Corresponding simulations of the BeVERLI Hill at the reduced scale of the risk-reduction experiment similarly identified a windward separation region for the 0 deg yaw case, albeit at the windward foot of the hill, and showed that it would disappear at high Reynolds numbers and extend only into the laminar sublayer when present (refer to [31] and [45]). Subsequent reduced-scale risk-reduction experiments confirmed the CFD observations and demonstrated the separation region to wash out entirely at the higher Reynolds numbers of interest to the validation experiment (compare Figs. 11(a) and 11(b)), effectively eliminating any risk. This finding also incentivized the addition of direct surface shear stress measurements using OFI to diagnose flow separation over the BeVERLI Hill better.

During the BeVERLI wind tunnel entries, qualitative and quantitative comparisons of the measured data with CFD were decisive in obtaining immediate feedback to adjust the measurements and test matrix. Experimental contour plots of hill surface C_p, as shown in Fig. 12, were compared qualitatively to their CFD equivalents (see, e.g., Fig. 9), receiving immediate feedback on the plausibility as well as quality of the experimental data. The creation of comprehensive C_p contours, like Fig. 12, depends on the interpolation scheme choice, e.g., cubic or linear, due to the discrete and unevenly distributed nature of experimental pressure tap layouts. Simultaneous comparative CFD data played a crucial role in assessing the impact and quality of different interpolation schemes, aiding in the selection of the optimal scheme (here linear) for the most representative C_p contour. The measurements of Fig. 12 also gave an immediate view of the flow's symmetry, showing that the wake asymmetry for the 45 deg yaw case is primarily triggered at higher Reynolds numbers in the experiment. As a result, the test matrix was adapted during the experiment to include repeated measurements at the higher tested Reynolds numbers, including focused unsteady surface pressure measurements using precisely located pressure taps at spanwise symmetric locations. Similarly, the pressure contours supported the placement of the PIV and LDV measurement locations as illustrated in Fig. 6.

More quantitative line plot comparisons of the surface pressure coefficient along the centerline and centerspan of the BeVERLI Hill to CFD conversely were useful in providing direct feedback regarding the turbulence models utilized in the computations. As illustrated in Figs. 13(a) and 13(c), both the k–ω SST and the SA model fail to predict flow separation occurring around $x1/H=1.5$⁠, which is a well-known and acknowledged weakness of standard eddy viscosity models. Perhaps less expected, however, is to see variation in the model predictions over the hill top, where the flow is fully attached. The variation is more pronounced at $ReH$ of 250,000, indicative that this area is connected to Reynolds number dependent physical characteristics. This valuable finding identified the hill top as a region containing complex physics capable of challenging turbulence models. Consequently, many of the measurements, such as PIV, were focused in that particular area. The spanwise C_p distribution indicates that both tested RANS models are capable of predicting flow asymmetry, but there are appreciable variations in the shape and form of the asymmetry with turbulence model and numerical grid. Sometimes the RANS models predict the asymmetry when it does not occur in the experiments and vice versa. This leaves the important question open whether the occurrence of asymmetry for a particular turbulence model at any given flow condition is the manifestation of the true physics of the problem or triggered purely by numerical conditions. Are the models succeeding in predicting the true underlying physical problem, but the numerical grids are yet too coarse to capture the physics responsible for features like the wake asymmetry, or are they predicting a different physical problem that is driven by the numerics embedded in a particular CFD solver? Regardless of a definitive answer to the preceding question, these comparisons show the value of the BeVERLI Hill data in identifying weaknesses and improving turbulence models and in identifying the role of numerical error and grid sensitivity.

It is noteworthy that the analysis of asymmetry relied on steady and unsteady surface static pressure and oil flow visualization data, as well as mean flow quantities extracted from RANS CFD data. This powerful combination of surface data and RANS CFD data provided a holistic view of the flow asymmetry, enabling the identification and understanding of its presence, behavior, sensitivities, and initial explanations for its physical origins. However, to gain a more comprehensive understanding, especially pertaining the dynamics of the asymmetry, it is necessary to include additional field mean and instantaneous velocity and turbulence data from techniques like PIV and LDV. As of September 2022, the available PIV and LDV data, primarily localized and unilateral to the hill, did not provide further insights beyond what the surface data already offered. To address remaining questions about the flow asymmetry, future investigations in BeVERLI will consider additional PIV and LDV data, especially spanning both sides of the hill, as well as volumetric data from measurement techniques such as particle tracking velocimetry (PTV), the latter which could also provide reliable estiamtes of field pressure. Furthermore, involving researchers who can provide high-fidelity simulations (LES or even DNS) is expected to bring further clarity to the flow asymmetry.

Incorporating time-dependent measurements is highly advisable, certainly necessary to reach high levels of experimental completeness per Oberkampf and Smith's guidelines in Ref. [1], whether or not capturing unsteady flow characteristics is a primary target of the validation study. The repeatability and consistency of the data in a validation experiment are critical, and using time-dependent measurements allows assessing these qualities in more depth. At a minimum, they should be used to confirm that the flow is not unsteady in the mean.

Time-dependent measurements were another important asset of project BeVERLI. Initially, only mean static pressure data were acquired on the hill surface, and instantaneous pressure measurements were a later addition. The wake asymmetry raised many concerns. As previously discussed, while the asymmetry itself is not a problem, its stability is crucial. An appreciable likelihood existed that the wake would be bi-modal and rapidly switch between different asymmetric modes, each favoring a distinct side of the hill. If the time scale at which the switching events happened was smaller than the measurement time scale, a stable mean flow field could never be measured, jeopardizing data repeatability and consistency. The first concerns were raised after observing the asymmetric wake experimentally for the first time in February 2020 for the 0 deg yaw case, which up to that point represented the main validation case. As detailed later, the mean pressure measurements indicated a long-time switching of the wake for equivalent symmetric orientations of the 0 deg yaw case, which was, however, attributed to the appreciable geometric deviations from nominal of the initial acrylic model used at that time. Notwithstanding these conclusions, it could not be excluded that the wake switching instead represented an inherent unsteadiness of the 0 deg yaw case and that the 45 deg yaw case was not similarly affected.

Introducing the current, more precise hill models (Fig. 2) and unsteady surface pressure measurements for the subsequent May 2021 test revealed vital information regarding the wake asymmetry. Figures 14(a) and 14(b) show unsteady surface pressure measurements for the 0 deg and 45 deg yaw case, respectively, as a moving-average-subtracted pressure coefficient, $Cp−Cp,avg$⁠. At each orientation the data are measured in the wake and presented at three specific pressure port locations: one on the centerline and two symmetrically about the $x3=0$ plane. Each of the three ports is associated with a different curve in Fig. 14. Additionally, in Fig. 14, the data are accompanied by a top-view outline of the hill with the pressure port locations marked as dots. Furthermore, circles outline the locations of the pressure ports associated with the presented data. The data sets were acquired for a total duration of 60 s at a sampling frequency of 300 Hz. Because the pressure ports were not calibrated, and their frequency response was therefore unknown, these results are qualitative in nature. Nevertheless, the additional benefit that these unsteady measurements provide is irrefutable. The results for the 0 deg yaw case indicate that the asymmetric wake is unstable and changes sides chaotically as the unsteady pressure on either side of the hill can be seen to be almost exactly out of phase (see Fig. 14(a)). The switch occurs at roughly a time scale of 5 s, which is clearly below any facility-driven phenomenon. Additional cross-facility experiments [21], performed in SINTEF Ocean's water tunnel in Norway, have further ruled out any facility biases. A corresponding Strouhal number, St, based on the hill height and the freestream velocity was computed for both the Virginia Tech and the SINTEF experiment, resulting in a closely matching, very low value of 0.003. On the other hand, the 45 deg yaw case does not display any switching of the wake asymmetry. Because of these results, the 0 deg yaw case was immediately excluded as a primary validation case in favor of the stable 45 deg yaw case for the remainder of the May 2021 and the subsequent October 2021 test.

Before conducting the validation experiment, it is essential to estimate as-tested rather than nominal boundary conditions. Quantifying deviations of as-tested conditions from nominal allows interpreting and distinguishing anomalies or off-nominal behavior from actual physical phenomena in the collected data. It further offers indications of the provenance of these anomalies, enabling their mitigation or, at the very least, raising awareness. As a result, the risk of accidentally eliminating false, apparent outliers in the data is minimized. Furthermore, running supporting CFD simulations using as-tested boundary conditions yields vital information about the flow's sensitivity to boundary condition variations, which can in turn be used to inform the design of the validation experiment. The use of as-tested boundary conditions in the simulations can also lead to potentially more insightful and accurate predictions closer to the actual experiment.

Scans of the manufactured model geometry are an excellent tool for accomplishing the above. More importantly, they provide means of gauging the quality of a manufactured model, the utmost purpose of which is the accurate reproduction of the as-designed geometry. As described by Rhode and Oberkampf [46], physical departures of the model surface coordinates from the as-designed geometry, whether due to errors in the model fabrication, damage to the surface, distortion and bending from pressure and temperature loads, or variations in surface roughness that could lead to flow field asymmetries and discontinuities, directly impact measurements as model and instrumentation uncertainties. Model scans support the quantification of these uncertainties and the assessment of a flow field's sensitivity to model geometry deviations. In addition, manufacturing models for a large-scale validation experiment is generally tied to considerable investments, and catching flaws in a model early on is essential for the experiment's success, especially when prototyping is not an option.

Model scans have led to the decision to discard the first manufactured model of the BeVERLI Hill in favor of the two current, more precise CNC-machined models (see Fig. 2), despite the significant investment. The first model was manufactured entirely from clear acrylic using a vacuum-forming technique. The reason behind constructing the hill entirely out of clear acrylic was to obtain virtually unlimited optical access for PIV and LDV. As reported by Gargiulo et al. [47], model scans performed using two different scanners (a laser and a contact-probe scanner) independently showed regions of significant deviation of the acrylic model from the as-designed geometry, within a range of up to ±6 mm. More importantly, some of the most significant perturbations were asymmetrically distributed across the hill surface. These large excursions in the acrylic model geometry could be directly related to variations observed in the BeVERLI Hill surface pressure data, including changes in the spanwise location of the wake asymmetry. As a result, it was no longer possible to conclude definitively whether the symmetry breaking predicted by the nominal CFD simulations was the same phenomenon observed in the experiments. Upon receiving the two current high-precision CNC-machined models previously presented in Fig. 2, the two models were scanned, indicating much lower deviations within ±0.4 mm from the as-designed geometry. These highly precise CNC models showed no relationship to the observed pressure variations within experimental uncertainties in subsequent experiments. The continued presence and consistency of symmetry breaking and bi-modality in the experiments could then be confidently categorized as a genuine feature of the BeVERLI Hill. The findings also provided further insight into the sensitivity of the flow, showing that the model geometry is not a dominant source of asymmetry.

A validation experiment should be conducted at multiple testing facilities whenever possible. In doing so, facility biases that alter the nominal boundary conditions and skew the experimental results, such as damages in the wind tunnel circuit or any anomalies in the wind tunnel components, can be identified. Key sensitivities to boundary conditions can be efficiently studied using cross-facility experiments. Similarly, allowing CFD simulations to be performed collaboratively between multiple groups can speed up the study of sensitivities of the numerical solutions. Different numerical solvers, turbulence models, computational grids, and numerical schemes can be explored in a remarkably shorter time than when performing simulations alone.

The discussed symmetry breaking of the BeVERLI Hill had an appreciable potential of being affected by the primary testing facility at Virginia Tech. In particular, the wake bi-modality observed for the 0 deg yaw case and its associated characteristic frequency could be driven by the specific testing facility and either occur at a different time scale or be completely absent in a different facility. To investigate this hypothesis, the BeVERLI Hill was tested at multiple facilities around the globe as part of the NATO AVT-349 program, namely, at UTIAS in Canada and SINTEF Ocean in Norway. These studies are detailed in Ref. [21] and, as discussed previously, demonstrated that not only was the symmetry breaking present in all those facilities for the hill orientations and Reynolds numbers tested at Virginia Tech but the unsteady, bi-modal wake oscillations observed for the 0 deg yaw case were nearly identical, demonstrating that the different experiments described a real feature of the BeVERLI Hill.

The BeVERLI Hill was also studied numerically as part of NATO AVT-349 collaborations, as detailed in Ref. [22]. Five groups simulated the case across six solvers, three sets of numerical grids, and seven RANS closure models, including physics-based and data-driven models. This large-scale numerical collaboration identified significant sensitivities of the RANS solutions to the computational grids, the turbulence models, and the solver. While not rigorously isolated, the iterative convergence appeared to additionally have a significant impact on the solutions.

Like synergism, redundancy and repetition of measurements enhance a validation experiment. They can be even seen as a requirement directly impacting the completeness of a validation experiment, as defined by Oberkampf and Smith [1]. Since no measurement is free of error, and no single measurement technique is best for all flow applications and ranges of parameters, redundant measurements should be performed whenever possible. Indeed, they should be considered if there is a suspicion that a measurement technique is of questionable applicability under the tested conditions. Having multiple diagnostic techniques allows to confirm the uncertainty analysis performed for each measurement, ultimately leading one to learn new things about the uncertainty of the measurements, especially bias errors, which are very difficult to identify. On the other hand, the repeatability of measurements can strengthen, confirm, and improve hypotheses formulated from the data about a particular flow mechanism or physical behavior. It further provides rigorous means of quantifying the measurement uncertainty and identifying flow sensitivities.

BeVERLI's experimental run matrix was designed to take advantage of these principles in multiple ways. For instance, the surface pressure data were acquired repeatedly at each condition from a run matrix combining three Reynolds numbers and two hill orientation groups (0 deg, 90 deg, 180 deg, 270 deg and 45 deg, 135 deg, 225 deg, 315 deg). At each run, a specific hill orientation would be tested and kept constant while the Reynolds number would be cycled from the lowest (⁠$ReH≈$ 250,000) to the highest (⁠$ReH≈$ 650,000) test value. Each run would be completed after one upwards and one downwards Reynolds number sweep. This Reynolds number cycling protocol enabled the acquisition of data at repeated conditions, capturing data variance related to the facility and instrumentation. Each orientation was tested several times during an entire measurement campaign to also capture data variance due to model/instrumentation uncertainty, hourly or daily atmospheric variation, and the reproducibility of the experimental setup. Because of the abundance of repeated pressure data in project BeVERLI, rigorous UQ using advanced methods, e.g., as proposed by Aeschliman and Oberkampf [10], became possible.

Sensitivities of the flow over the BeVERLI Hill to model misalignment, the as-manufactured model geometry deviations, and other boundary conditions could be successfully isolated using the repeated mean surface pressure data. Studying complex and sensitive flow phenomena, such as the wake asymmetry, could only be conducted thanks to redundancy and repetition of measurements. The consistency between the repeated measurements and the equivalent observations of asymmetry in redundant measurements increased the confidence that the observed wake asymmetry was an inherent characteristic of the BeVERLI Hill. Moreover, sensible hypotheses could be established about the origin and mechanism of this phenomenon thanks to the abundance of consistent and repeated data.

The benefit of redundant measurements is shown using velocity measurements from the PIV–LDV-overlap station located at $(x1, x3)=(−0.15 m,−0.15 m)$⁠. In Fig. 15, the PIV and LDV data are compared alongside corresponding RANS results. The data are presented in a local Cartesian coordinate frame (prime symbol), whose origin is on the hill's surface at $(x1, x3)=(−0.15 m,−0.15 m)$⁠. The $x1′$-axis of the coordinate system is locally tangential to the LDV windows and points toward the hill top, the $x2′$-axis is aligned with the local wall-normal direction, and the $x3′$-axis completes the system in the right-handed sense. Presented are wall-normal profiles of the mean velocity component in the $x1′$-direction, $U¯1′$⁠, and the Reynolds normal stress component in the $x2′$-direction, $u22¯′$⁠. Note that the ordinate of the plots in Fig. 15 is scaled logarithmically to better represent the data close to the hill surface. An estimate of the uncertainty due to random error and bias error from various sources, $ULDVRandom & Bias$⁠, detailed in [34], is provided for the LDV data and represented as 95% confidence interval error bars about the LDV data in the corresponding figure. For the PIV data, only estimates of the uncertainty due to random error, $UPIVRandom$⁠, according to Ref. [48], are provided. These estimates are represented as a 95% confidence interval shaded surface about the PIV data. Random error for PIV is typically small due to the large amount of independent samples (⁠$>10,000$⁠) available; thus, it is only moderately visible in Fig. 15. The estimation of uncertainty due to bias error for PIV is knowingly a very difficult and involved effort (see, e.g., [49]) and remained inaccessible for the present study. However, as discussed next, in several instances, e.g., when redundant data were available, valuable information about the PIV bias error could be extracted.

The $U¯1′$ data show how PIV, primarily focused in the outer layer of the flow due to the limited spatial resolution of the system, is nicely complemented by the near-wall LDV data, which in turn have limited accuracy in the outer layer due to a reduced number of available samples. The reasonably smooth transition between the two data sets indicates that the two systems yield consistent mean velocity measurements. Discrepancies are observed for the outer points and inner points of the LDV and PIV data,respectively. However, the PIV data still fall within the estimated LDV uncertainty bounds, leading to the following implications. Assuming that the PIV data provide more accurate mean flow measurements in the discrepancy region, where LDV is subjected to a larger random error from the reduced number of samples, the results imply that the uncertainty analysis for LDV provides reasonable uncertainty bounds and that the PIV data closest to the wall still offers reliable mean velocity measurements.

This estimate is represented as a hatched surface in Fig. 15. These findings can be used as a valuable basis to compare future quantitative bias error estimations to assess their plausibility and to gauge bias error where a rigorous estimate is unavailable.

Another benefit of redundant measurements is apparent when incorporating additional measurements across multiple wind tunnel entries. In the case of BeVERLI, shear stress measurements using OFI were included in later entries. To ensure validity and quality, the OFI shear stress measurements were compared to redundantly extracted values from LDV data and RANS simulations at overlapping locations. LDV-based shear stress estimation utilized the surface normal velocity gradient (⁠$∂U¯1/∂x2$⁠) computed with the nearest data points to the hill surface. Figure 16 shows the comparison of these measurements at three overlapping locations along the dedicated acrylic LDV windows: location 1 (⁠$x1/H, x3/H$⁠) = (−1.42, −1.42), location 2 (−1.29, −1.29), and location 3 (−0.80, −0.80). These measurements significantly overlap within the indicated uncertainty range, providing strong evidence of the adequate performance of the OFI measurement and enhancing confidence in the obtained OFI measurements' accuracy and reliability.

We presented strategies and tools for designing and conducting CFD model validation experiments based on the successes and lessons learned during project BeVERLI. These strategies aim to promote the highest levels of MVEC per Oberkampf and Smith [1]. The key findings from this study are:

1. The proposed strategies effectively mitigate significant risks in validation experiments, such as inapt flow features, opportunity costs, and anomalies in the experimental setup and boundary conditions.

2. Applying the exemplified principles supports the adequate design of validation experiments, including an experimental run matrix for the highest levels of data completeness, documentation, and UQ.

3. The combined application of the suggested strategies enables the systematic identification, diagnostic, UQ, and physical characterization of critical flow features, as-tested boundary conditions, and sensitivities.

4. The proposed ideas apply to estimating critical and difficult-to-obtain bias error uncertainties of different measurement systems, e.g., the underprediction of high-order statistical moments of the velocity field data from PIV due to spatial filtering effects, and the quality assessment of uncertainty estimates.

5. Synergistic CFD simulations are essential in guiding the experiment design, identifying and characterizing critical flow features and sensitivities, monitoring measurement consistency, and studying turbulence models' predictive capability and weaknesses.

The presented strategies led to key achievements for BeVERLI, including correctly identifying and characterizing the symmetry-breaking and bi-modal features of the hill flow, representing a significant risk for the experiment. The 45 deg yaw case produces steady asymmetry, while the 0 deg case exhibits unsteady bi-modality. The symmetry breaking and bi-modality are thought to be a consequence of the flow's inherent metastability and a genuine feature of the hill geometry, not an experimental setup artifact. It is unknown whether the symmetry breaking is triggered by small boundary condition asymmetries or an inherent flow mechanism, like a shear layer instability on the hill's leeward side from the interaction of the distinctly developing boundary layer over the hill's top and flanks. This open question will be investigated in the future. The strategies made BeVERLI achieve high initial MVEC levels, reflected in the completeness of data, documentation, and UQ, including bias error estimates often omitted in similar studies. CFD simulations of the BeVERLI Hill enabled a rigorous characterization of the flow and established the baseline capabilities and weaknesses of linear eddy viscosity models, specifically the SA and SST model.

Looking ahead, project BeVERLI will conduct additional wind tunnel tests that focus on non-symmetric hill orientations and involve comprehensive measurements such as pressure, PIV, LDV, and OFI. Additional volumetric data through PTV will be included. The BeVERLI database, currently available upon request to the present authors, will be made publicly accessible in collaboration with NASA. It will contain reduced data, raw data, experimental facility details, boundary condition information, and other relevant components. The data archive will likely be hosted on a cloud service and accessible through an official website for convenient access. The NASA TM publication, expected by the end of 2023, will outline the detailed archival procedures and plans. Additionally, a public blind CFD model validation challenge⁴ has been launched for the 30 deg yaw case. Participants are asked to submit simulations of the case using their choice of solver and turbulence models. The boundary conditions and numerical grids have been provided. Upon submission, additional validation data for surface pressure, skin friction, and velocity field will be provided.

The authors thank NASA, in particular Dr. Michael Kegerise, Dr. Mujeeb Malik, and Dr. Christopher Rumsey, for their support under grant 80NSSC18M0146 and 80NSSC22M0061. We thank BeVERLI collaborators Vignesh Sundarraj, Tom Hallock, and Thomas Ozoroski for their assistance with wind tunnel testing and thoughtful discussion on the data presented. We also thank Stability Wind Tunnel engineer, Bill Oetjens, for his assistance and expertise during wind tunnel testing, and the Aerospace and Ocean Engineering Machine Shop staff for their assistance in the fabrication and instrumentation of test equipment.

• National Aeronautics and Space Administration (NASA) (Grant Nos. 80NSSC18M0146 and 80NSSC22M0061).

• Langley Research Center (Award Nos. 80NSSC18M0146 and 80NSSC22M0061; Funder ID: 10.13039/100006199).

C_f =

surface skin friction coefficient

C_p =

pressure coefficient

H =

hill height (m)

I =

turbulence intensity (%)

k =

turbulent kinetic energy (⁠$m2s−2$⁠)

P =

pressure (Pa)

T =

temperature (K)

U_i =

velocity components (m s^–1)

$uτ$ =

friction velocity (m s^–1)

x_i =

Cartesian coordinates (m)

Greek Symbols

μ =

absolute viscosity (kg m^–1 s^–1)

μ_t =

turbulent eddy viscosity (kg m^–1 s^–1)

ν =

kinematic viscosity (m² s^–1)

ρ =

density (kg m^–3)

ω =

specific turbulence dissipation rate (s^–1)

τ_w =

skin friction (N m^{– 2})

Dimensionless Groups

M =

Mach number

Re =

Reynolds number

Superscripts and Subscripts

0 =

stagnation

$∞$ =

free stream value

in, out =

inlet, outlet

ref =

reference value

Table 6 provides an overview of BeVERLI's MVEC assessment. The assessment evaluates the completeness of the BeVERLI model validation experiment's documentation up until September 2022 across the six categories, or attributes, suggested by Oberkampf and Smith [1]: Experimental facility (EF), analog instrumentation and signal processing (AI&SP), boundary and initial conditions (B&IC), fluid and material properties (F&MP), test conditions (TC), and system response quantities (SRQ). The primary goal of this assessment is to determine the extent to which the necessary information is provided for the modeling of the test case. To determine each attribute's completeness level, which ranges from 0 (lowest) to 3 (highest), Oberkampf and Smith [1] define a set of descriptors that must be fulfilled at each level. These descriptors serve as queries to evaluate the amount and quality of descriptive documentation, measured data, and corresponding UQ. In order to achieve a particular completeness level for an attribute, all descriptors for that level must be fulfilled, in addition to all descriptors from lower completeness levels. If only some descriptors at a certain level have been fulfilled, a fractional score is assigned to the lower completeness level instead of a whole number score. In Table 6, the highest completeness level achieved for each attribute, where all descriptors were satisfied, is labeled with the letter A (assessed). If the descriptors were only partially fulfilled at a certain level, the number of satisfied descriptors is indicated as x out of y or x/y, where x and y are real numbers. Descriptors are considered either fully satisfied (scored as 1) or partially satisfied (scored as 0.5). The final fractional score, indicating the completeness level for each attribute, is provided with two decimal places of precision.

must hold, where $rnominal,i$ is the expected nominal grid refinement factor of grid level i. On the other hand, if the three effective grid refinement factors are unequal, the grids are not (exactly) geometrically similar. The effective refinement factors for the grids of the present study and their error relative to the nominal value, $εk,i=(rk,i−rnominal,i)/rnominal,i$⁠, for k = 1, 2, 3, are shown in Table 7. While the different effective grid refinement factors are not exactly equal at every grid level, they are found well within 0.5% of their corresponding nominal value, demonstrating the successful systematic refinement of the grids.

where $r=h2/h1$ is the grid refinement factor.

Note, this error estimate approaches infinity as $p→0$⁠, resulting in an unrealistically large error estimate. Therefore, the OOA should be limited to some small positive number, p_l. Phillips and Roy [54] recommend limiting the minimum order of accuracy to 0.5. In most practical cases $p=p̂$ and $0.5<p<pf$⁠. Further note, p_f represents the formal order of accuracy of a scheme only when the solution is smooth. The presence of singularities or discontinuities reduces the formal order of accuracy.

where N is the number of Richardson nodes. Note, because the formal order of the numerical schemes used in this study is only second-order (p_f = 2) for entirely smooth solutions, but the flow over the BeVERLI Hill is not due to its separation characteristics and recirculation regions, the formal order of accuracy is expected to be between $pf/(pf+1)=2/3$ and p_f = 2. Thus, a value of $pf=1.5$ was selected in the present study.

which uses the solution on grid Level 1 and Level 2, denoted f₁ and f₂, respectively. Due to the oscillatory convergence of the BeVERLI solutions a conservative factor of safety, F_s = 3, was chosen in the present study.

The grid study for the presented RANS simulations is summarized in Figs. 17 and 18 in terms of C_p and C_f hill surface distributions of the 45 deg yaw case for $ReH≈$ 250,000 and 650,000, respectively. The resulting OOA for these distributions are summarized in Table 8, with the corresponding $UDECp$ and $UDECf$ distributions plotted in Figs. 20 and 21 for $ReH≈$ 250,000, and in Figs. 22 and 23 for $ReH≈$ 650,000. Furthermore, representative uncertainty values in terms of the root-mean-square (rms) and maximum (max) of the distributions across all grid levels are reported in Table 9.

Notably, at $ReH≈650,000$⁠, the behavior of the Level 1 SA solution (see Fig. 18) differs from the solution on the other grids. As clearly visible from the hill surface contour plots in Fig. 19, the asymmetric flow solution predicted on the coarse grid levels is contrasted by the nearly symmetric solution on the Level 1 grid. At $ReH≈$ 250,000 the SST model has a higher DE in C_p and C_f than the SA model. Overall, the Level 2 and the Level 1 mesh provide reasonable error levels at $ReH≈$ 250,000. For the higher $ReH≈$ 650,000 case, the SST model shows lower uncertainty values than the SA model across the board, owing to the consistency of the solutions. However, again, the Level 2 and the Level 1 mesh are deemed adequate for the analysis at $ReH≈650,000$⁠.