A Gargiulo

T A Ozoroski

T Hallock

A Haghiri

R D Sandberg

M Visonneau

G Deng

E Guilmineau

D Geneau

T Jeans

M Bettle

M Kerkvliet ¶ ¶

S L Toxopeus

J E Duetsch-Patel

V Sundarraj

A Borgoltz

W J Devenport

C J Roy

K T Lowe

Computations of the BeVERLI Hill three-dimensional separating flow model validation cases

The BeVERLI (Benchmark Validation Experiment for RANS/LES Investigations) Hill project aims at producing a detailed experimental database of three-dimensional non-equilibrium turbulent boundary layers with various levels of separation while meeting the most exacting requirements of computational fluid dynamics validation as per Oberkampf and Smith [1]. A group of the Science and Technology Organization (STO) of the North Atlantic Treaty Organization (NATO) entitled NATO AVT-349 -"Non-Equilibrium Turbulent Boundary Layers in High Reynolds Number Flow at Incompressible Conditions" has recently considered the BeVERLI Hill case. Their goal is to advance the accuracy and range of prediction models for high Reynolds number non-equilibrium boundary layers. This highly collaborative and international group comprised of various academic, governmental, and industrial institutions has performed several Reynolds-averaged Navier-Stokes (RANS) simulations of the BeVERLI Hill using different grids, solvers, and turbulence models. The resulting solutions and available experimental data are presented in this paper to summarize and highlight key features, sensitivities, and the current predictive capability of RANS for the BeVERLI Hill case. The results suggest important sensitivities to Reynolds number, grid density, iterative converge,

III. Introduction

T he computational resources required to perform highly-resolved DNS and LES of full-scale aero-mechanical devices in turbulent flow at high Reynolds numbers remain out of reach and are not likely to be available for many decades [START_REF] Moin | Tackling Turbulence with Supercomputers[END_REF][START_REF] Spalart | Philosophies and Fallacies in Turbulence Modeling[END_REF][START_REF] Duraisamy | Status, Emerging Ideas and Future Directions of Turbulence Modeling Research in Aeronautics[END_REF]. As a result, the development of more efficient approaches, such as mean flow resolving RANS equations and hybrid RANS/LES methods, along with the turbulence models upon which these approaches rely, is at the center of current CFD and turbulence research. Industry's desire to increasingly utilize CFD for high-impact decisions [START_REF] Bolds-Moorehead | Aircraft Certification by Simulation[END_REF] further drives the demand for efficient and reliable predictions in more complex flow regimes than those for which RANS turbulence models are currently trusted.

The development of improved turbulence models demands validation data from benchmark experiments designed to accurately capture the physics and critical features of practically-relevant flows. Modelers must carefully inspect the source of validation data before use. Many experimental results indeed lack some level of critical completeness. Such experiments, for instance, provide an insufficient amount or completely omit uncertainty estimates, ignore inflow nonuniformities, or assume equivalence between as-designed and as-measured geometries. These simplified assumptions propagate as unknown uncertainties that directly impact the boundary conditions used to setup corresponding CFD computations [START_REF] Oberkampf | Verification and Validation in Scientific Computing[END_REF][START_REF] Rumsey | The NASA Juncture Flow Test as a Model for Effective CFD/Experimental Collaboration[END_REF].

An ongoing effort contributing to the mitigation of this problem and the acquisition of a benchmark validation dataset for CFD arises from a collaboration between the National Aeronautics and Space Administration (NASA) and Virginia Tech (VT) [START_REF] Lowe | Status of the NASA/Virginia Tech Benchmark Experiments for CFD Validation[END_REF][START_REF] Gargiulo | Examination of Flow Sensitivities in Turbulence Model Validation Experiments[END_REF][START_REF] Gargiulo | Flow Field Features of the BEVERLI Hill Model[END_REF]. The case analyzed is a three-dimensional wall-mounted hill configuration named BeVERLI Hill, which captures many features shared by applications involving three-dimensional flow separation. Similar studies involving flow over hills were conducted in the past by Simpson et al. [START_REF] Simpson | Study of Vortical Separation from an Axisymmetric Hill[END_REF], Byun and Simpson [START_REF] Byun | Structure of Three-Dimensional Separated Flow on an Axisymmetric Bump[END_REF], in the NASA Faith Hill experiments by Bell et al. [START_REF] Bell | Surface and Flow Field Measurements on the FAITH Hill Model[END_REF] and Husen et al. [START_REF] Husen | Global Luminescent Oil-Film Skin-Friction Meter Generalized to Three-Dimensional Geometry and Applied to FAITH Hill[END_REF], and more recently in another ongoing collaborative effort between the University of Washington and Boeing [START_REF] Williams | Experimental Study of a CFD Validation Test Case for Turbulent Separated Flows[END_REF][START_REF] Williams | Characterization of Separated Flowfield over Gaussian Speed-Bump CFD Validation Geometry[END_REF][START_REF] Robbins | Overview of Validation Completeness for Gaussian Speed-Bump Separated Flow Experiments[END_REF]. Other studies include two-dimensional hills with a range of shapes [START_REF] Webster | Turbulence Characteristics of a Boundary Layer over a Swept Bump[END_REF][START_REF] Greenblatt | Separation Control over a Wall-Mounted Hump[END_REF], axisymmetric aft body flows [START_REF] Bachalo | Transonic, Turbulent Boundary-Layer Separation Generated on an Axisymmetric Flow Model[END_REF][START_REF] Disotell | Design of an Axisymmetric Afterbody Test Case for CFD Validation[END_REF], and backward-facing ramps [START_REF] Simmons | Benchmark Smooth Body Flow Separation Experiments[END_REF]. The BeVERLI Hill project is one of the first to be undertaken with the intent of achieving the highest levels of completeness as defined in [START_REF] Oberkampf | Assessment Criteria for Computational Fluid Dynamics Model Validation Experiments[END_REF], including the complete documentation and uncertainty estimation of their database.

Recently, the Science and Technology Organization (STO) of the North Atlantic Treaty Organization (NATO) has integrated the BeVERLI Hill case into their technical activities. The NATO AVT-349 group focusing on "Non-Equilibrium Turbulent Boundary Layers in High Reynolds Number Flow at Incompressible Conditions" aims to advance the accuracy and range of prediction models for high Reynolds number non-equilibrium boundary layers. The group consists of members and collaborators around the globe from academia, government and non-government laboratories, and industry. The members of the NATO AVT-349 group are performing RANS computations of the BeVERLI Hill benchmark validation case as part of their computational activities. They are testing a large variety of standard linear, non-linear, and in-house turbulence models, different CFD solvers, and computational grids. In this paper, initial results from this ongoing collaborative effort available to date are presented to summarize and highlight key features, sensitivities, and the current predictive capability of RANS for the BeVERLI Hill case.

IV. Methodology

A. Geometry

The BeVERLI Hill geometry is defined in Gargiulo et al. [START_REF] Gargiulo | Examination of Flow Sensitivities in Turbulence Model Validation Experiments[END_REF] and represented in Fig. 1. It features a 5th-degree polynomial centerline profile, superelliptic corners, and a squared flat top of width 𝑠 = 0.093472 m. The hill width is 𝑤 = 0.93472 m and the height 𝐻 = 0.186944 m (aspect ratio = 𝑤/𝐻 = 5).

B. Simulation Cases

The current BeVERLI Hill RANS simulations were performed at flow conditions corresponding to hill-height-based Reynolds numbers of Re 𝐻 = 250, 000 and 650,000. The latter coincides with the target-scale Reynolds number selected in the original experiment performed at VT. One orientation of the hill with respect to the oncoming boundary layer was considered, namely the 45 • yaw angle orientation, as schematically illustrated in Fig. 2.

Fig. 2 BeVERLI Hill model yaw angle orientations with respect to the oncoming flow.

VT, being the primary institution leading the experimental efforts of the BeVERLI Hill project, has laid the groundwork for the corresponding RANS simulations. They proposed an initial modeling strategy, boundary conditions, reference conditions, and computational grids for the BeVERLI Hill case. This foundation was kept as consistent as possible amongst all participating members of the NATO AVT-349 group. However, solver-specific parameters and differences in boundary conditions, alternative modeling strategies, and in some cases different computational grids, were employed across institutions. The numerical methodology by VT is introduced next, followed by a detailed description of the configurations applied by the other participating institutions.

C. Virginia Tech

Numerical Solver

VT utilized two distinct CFD solvers for their simulations. First, the VT in-house solver SENSEI [START_REF] Jackson | Turbulence Model Implementation and Verification in the SENSEI CFD Code[END_REF][START_REF] Derlaga | SENSEI Computational Fluid Dynamics Code: A Case Study in Modern FORTRAN Software Development[END_REF] and second, the commercially available ANSYS Fluent solver [START_REF] Ansys | FLUENT Manual[END_REF]. In SENSEI, the inviscid flux discretization was performed using Roe's flux difference splitting [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes[END_REF] with MUSCL extrapolation to achieve a second-order accurate scheme. The viscous and turbulent fluxes were discretized to second-order accuracy using a central flux scheme after applying the Green-Gauss theorem. In ANSYS Fluent, the available density-based solver using an algebraic multigrid method (AMG) was used to solve the RANS equations numerically. The spatial discertization of the inviscid fluxes was performed using a second-order updwind scheme. Fluent achieves second-order accuracy using a multidimensional linear reconstruction approach [START_REF] Barth | The Design and Application of Upwind Schemes on Unstructured Meshes[END_REF]. The viscous and turbulent fluxes in Fluent were discretized to second-order accuracy using a central flux scheme and the Green-Gauss theorem.

Boundary and Reference Conditions

The computational domain employed by VT is represented in Fig. 3 and models the as-designed geometry of the BeVERLI Hill in the Virginia Tech Stability Wind Tunnel (VT SWT) test section, as depicted in Fig. 4, where the original experiment takes place. The VT SWT is a continuous, single return, subsonic wind tunnel. Its test section is nominally a rectangular cuboid 1.83 m by 1.83 m in cross section and 7.30 m long. The BeVERLI Hill is nominally mounted at the center of one of the walls of the test section. In the simulations, this wall coincides with the floor of the computational domain. A Cartesian coordinate system, the origin of which is located at the center of the hill in the plane of the wall hosting the geometry, is equally utilized in the experiment and the simulations. The 𝑥-axis coincides with the streamwise flow direction, the 𝑦-axis is in wall-normal direction, and the 𝑧-axis is chosen in spanwise direction to complete a right-handed coordinate system. The inflow plane of the computational domain is extended in the simulations to provide a similar incoming boundary layer as experimentally measured. The boundary layer thickness at 𝑥 = -1.83 m is 0.02 and 0.05 m for Re H = 250, 000 and 650,000, respectively. Although quite close, highly accurate matching of the experimental boundary layer parameters was not performed in the present work but is planned for future studies. The computational outflow plane location is selected to fall outside of the domain of influence of the hill. The boundary conditions as employed by VT for the SENSEI and ANSYS Fluent solvers for the hill-height-based Reynolds numbers tested, Re 𝐻 = 250, 000 and 650,000, are summarized in Fig. 3. The inlet boundary conditions are based off of values obtained from the original experiment. A Dirichlet condition specifying the stagnation temperature 𝑇 0,in and the stagnation pressure 𝑝 0,in is applied uniformly across the inflow domain. At the outlet, a Dirichlet condition fixes the static pressure 𝑝 out which is adjusted to match the reference pressure 𝑝 ref,exp (hereafter simply referred to as 𝑝 ref ) at a specified reference location of the experiment. The additional domain boundaries are modeled as smooth, adiabatic, viscous, non-permeable no-slip walls, and 𝑦 + values for all grids were below one. The flow of air within the system is modeled to be ideal, compressible, and viscous. The molecular dynamic viscosity 𝜇 is defined in accordance with Sutherland's law [START_REF] Sutherland | LII. The Viscosity of Gases and Molecular Force[END_REF]. Free-stream turbulence quantities are computed from the turbulence intensity 𝐼 and the turbulent viscosity ratio 𝜇 𝑡 /𝜇.

Reference values for the computation of normalized quantities, such as the pressure coefficient 𝐶 𝑝 or the friction coefficient 𝐶 𝑓 , are obtained at a specific common reference location between experiment and simulations. The common reference ensures a suitable comparison between the computational and experimental results. The reference pressure 𝑝 ref , to match between experiment and CFD, is calculated by averaging the pressure measurement from 7 spanwise-distributed surface static pressure probes (or ports) located 2.23 m upstream of the hill center on the wall opposite of the hill model. The exact 𝑥, 𝑦, and 𝑧 coordinates of the reference ports are summarized in Table 1. The reference ports are also depicted as green dots in Fig. 3. A projected copy of the reference ports on the wall hosting the BeVERLI Hill is also shown in

Turbulence Models

Two turbulence models were used by the VT investigators. The first, is the standard Spalart-Allmaras (hereafter referred to as SA) [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF][START_REF]Turbulence Modeling Resource -The Spalart-Allmaras Turbulence Model[END_REF] one-equation eddy-viscosity model. The second, is the the 2003 Menter 𝑘-𝜔 Shear Stress Transport (hereafter referred to as SST2003) two-equation eddy-viscosity model [START_REF] Menter | Ten Years of Industrial Experience with the SST Turbulence Model[END_REF][START_REF]Turbulence Modeling Resource -The Menter Shear Stress Transport Turbulence Model[END_REF]. Note, SENSEI uses the SA-neg [START_REF] Allmaras | Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model[END_REF] implementation of the SA model, while Fluent uses the standard formulation.

Computational Grids and Iterative Convergence

A family of five structured curvilinear grid levels was utilized by the VT team. Snapshots of the first four levels are illustrated in Fig. 5. The meshing process was performed using the commercial software Pointwise [START_REF] Pointwise | Pointwise Manual[END_REF]. The typical mesh strategy employed used a denser, high resolving allocation of grid nodes in the regions of and around the BeVERLI Hill as well as at the boundaries of the computational domain representing the walls of the VT SWT. The latter was done to fully resolve the boundary layer on each tunnel wall. Furthermore, the grids were constructed on one half of the computational domain and mirrored about the 𝑧 = 0 plane to ensure full spanwise symmetry. The five grid levels generated by VT are identified from finest to coarsest as Levels 1, 2, 3, 4, and 7. The grid levels are 465 × 353 × 409 (Level 1), 369 × 281 × 325 (Level 2), 293 × 223 × 259 (Level 3), 233 × 177 × 205 (Level 4), and 117 × 89 × 103 (Level 7) nodes in size. Each level is generated starting from the finest Level 1 grid using a non-integer refinement factor of 𝑟 =

D. University of Melbourne

Numerical Solver and Boundary Conditions

The RANS calculations of the BeVERLI Hill case carried out by the University of Melbourne used OpenFOAM [START_REF] Weller | A Tensorial Approach to Computational Continuum Mechanics using Object-oriented Techniques[END_REF], which is a set of object-oriented open source code written in C++ based on the finite volume method to perform CFD. The solver used in OpenFOAM was simpleFOAM, which simulates a steady and incompressible, viscous flow. All terms including gradient, divergence and Laplacian were discretized using a second-order accurate central scheme, except for the divergence of the mean velocity and the turbulence variables k and 𝜔, for which a second-order accurate scheme and a first-order-accurate upwind-based scheme were applied, respectively.

The simulations focused on the case of Re 𝐻 = 250, 000. For the boundary conditions, at the walls (top, side and bottom), no-slip boundary conditions were used. At the outlet, zero-gradient boundary conditions were used for all the variables except for the pressure, which was fixed. At the inlet, the boundary conditions were set as follows: a zero-gradient pressure condition, a uniform velocity vector normalized with the free-stream velocity (𝑈 = 1, 𝑉 = 0, 𝑊 = 0), and uniform turbulence variables 𝑘 and 𝜔 calculated from the free-stream turbulence intensity of 0.02 % (taken from the VT experiment) and 𝜇 𝑡 /𝜇 = 1.5.

Turbulence Models

For the baseline RANS calculations, the SST2003 model was used as introduced above. A part of the NATO AVT-349 efforts include data-driven development of turbulence closures that aim to improve steady RANS predictions of statistically three-dimensional fluid dynamics problems. Here, we use an in-house machine learning tool that is based on gene-expression programming (GEP) [START_REF] Weatheritt | A Novel Evolutionary Algorithm Applied to Algebraic Modifications of the RANS Stress-Strain Relationship[END_REF]. A key objective is to find novel turbulence models that generalize well for certain classes of flows, thus the model development performed in the current study is purposely not based on the BeVERLI Hill data but rather on a high-fidelity dataset of a finite square-mounted cylinder with a height-to-width ratio of 4 and Re 𝑑 = 11, 000 (based on the cylinder width d), described in more detail in a previous study [START_REF] Weatheritt | Improved Junction Body Flow Modeling through Data-driven Symbolic Regression[END_REF]. The original GEP-based approach is extended to not only develop a non-linear explicit algebraic Reynolds stress model, but also a turbulent-kinetic-energy-corrective term, as proposed by Schmelzer et al. [START_REF] Schmelzer | Discovery of Algebraic Reynolds-stress Models using Sparse Symbolic Regression[END_REF]. The models generated with the GEP optimization framework are shown in Appendix A. The newly developed data-driven models were then used to extend the SST2003 turbulence model implemented in OpenFOAM and tested a-posteriori with RANS calculations of the current BeVERLI Hill case at Re 𝐻 = 250, 000 and 45 • yaw angle orientation, thus the testing of the novel closures is conducted on a testing case "unseen" during training.

Computational Grids and Iterative Convergence

An earlier set of meshes presented in [START_REF] Gargiulo | Flow Field Features of the BEVERLI Hill Model[END_REF] and generated by the VT group were used. In particular, two different mesh resolutions were utilized: a coarse mesh (Level 2) and a fine mesh (Level 0), consisting of about 11,000,000 and 41,000,000 grid cells, respectively.

OpenFOAM utilizes a normalized form of solution residual, 𝑟 = 1 𝑛 𝑛 |𝒃 -𝐴𝒙|, over all computational nodes 𝑛. The convergence criterion for the pressure and velocity fields was set to reach a residual tolerance fixed at 1e-6 (with relative tolerance of 0.01). The tolerance for the turbulence 𝑘 and 𝜔 components was fixed at 1e-7 and 1e-11 (with relative tolerance of 0.01), respectively.

E. University of New Brunswick and DRDC

Numerical Solver and Boundary Conditions

The University of New Brunswick (UNB), in partnership with Defense Research and Development Canada (DRDC), utilized the commercially available ANSYS CFX solver for all simulations. The incompressible RANS equations are discretized using a finite volume/finite element method with co-located pressure-velocity coupling and a coupled multigrid solver [START_REF] Schneider | Control Volume Finite-Element Method for Heat Transfer and Fluid Flow using Colocated Variables-1. Computational Procedure[END_REF][START_REF] Raw | A Coupled Algebraic Multigrid Method for the 3D Navier-Stokes Equations[END_REF].

All simulations were performed for the case of a hill-height-based Reynolds number of Re 𝐻 = 650, 000. The boundary conditions were selected accordingly. A prescribed inlet velocity of 55 m/s and density of 1.1272 kg/m 3 were selected. The domain reference pressure was selected to match the experimental inlet total pressure assuming incompressible flow, and the inlet turbulence intensity was set to 3 %. At the outlet, the average static pressure was set based on a theoretically calculated pressure drop across the domain. All remaining boundaries in the CFD domain were modeled as no-slip walls.

Turbulence Models

To investigate the impact of the turbulence closure model on the predicted flow field two turbulence models were selected. First, the standard linear two-equation eddy-viscosity model, the SST2003 model, was selected. Second, the anisotropic seven-equation baseline Reynolds Stress model (BSLRSM) [START_REF] Launder | Progress in the Development of a Reynolds-stress Turbulence Closure[END_REF] was utilized.

Computational Grids and Iterative Convergence

To estimate the discretization error, simulations were performed on four structured multi-block grids for each orientation and turbulence model, resulting in 16 simulations in the final dataset. An older set of initial grids provided by VT, as described in [START_REF] Gargiulo | Flow Field Features of the BEVERLI Hill Model[END_REF], were modified to improve the final convergence, while maintaining the same overall grid topology. Specifically, upstream of the BeVERLI Hill the axial spacing was redistributed such that it grew uniformly from the base of the hill to the inlet without increasing the total control volumes in this region. Downstream the axial spacing was increased by a factor of 1.3 and then redistributed such that it grew uniformly from the downstream base of the hill and then kept at a constant value up to the outlet once the aspect ratio of the first control volume at the wall reached 10,000. This increased the overall grid size relative to the original VT grids by 39 %.

All simulations were completed with second-order spatial accuracy. Primary target residuals for all simulation cases were set to an 𝐿 ∞ norm of 1e-4 with a secondary auxiliary convergence setting to an 𝐿 2 norm of 1e-6. For some simulations added gradient and mass relaxation was necessary to achieve the desired convergence. This was most prominent for the BSLRSM, which is likely due to the added anisotropy.

F. Maritime Research Institute Netherlands -MARIN

Numerical Solver, Boundary Conditions, and Turbulence Models

ReFRESCO [START_REF] Vaz | Free-Surface Viscous Flow Computations: Validation of URANS code FRESCO[END_REF] is a CFD solver based on a finite volume discretization of the continuity and momentum equations written in strong conservation form. The solver uses a fully-collocated arrangement and a face-based approach that enables the use of cells with an arbitrary number of faces. Picard linearization is applied and segregated or coupled approaches are available with mass conservation ensured using a SIMPLE-like algorithm [START_REF] Klaij | SIMPLE-Type Preconditioners for Cell-Centered, Colocated Finite Volume Discretization of Incompressible Reynolds-Averaged Navier-Stokes Equations[END_REF] and a pressure-weighted interpolation technique to avoid spurious oscillations [START_REF] Miller | Use of a Pressure-Weighted Interpolation Method for the Solution of the Incompressible Navier-Stokes Equations on a Nonstaggered Grid System[END_REF]. A volume of fluid technique is used for multiphase flows and several alternative mathematical formulations can be used to solve turbulent flow. Thorough code verification is performed for all releases of ReFRESCO [START_REF] Eça | On Code Verification of RANS Solvers[END_REF]. Earlier applications of ReFRESCO in the NATO context can be found in [START_REF] Toxopeus | Viscous-Flow Calculations for KVLCC2 in Manoeuvring Motion in Deep and Shallow Water[END_REF][START_REF]Assessment of Stability and Control Prediction Methods for NATO Air and Sea Vehicles[END_REF][START_REF]Reliable Prediction of Separated Flow Onset and Progression for Air and Sea Vehicles[END_REF][START_REF]Evaluation of Prediction Methods for Ship Maneuvering and Control[END_REF][START_REF]Assessment of Prediction Methods for Large Amplitude Dynamic Manoeuvres for Naval Vehicles[END_REF][START_REF] Lopes | Assessing Numerical Aspects of Transitional Flow Simulations Using the RANS Equations[END_REF][START_REF]Flowfield Prediction for Manoeuvring Underwater Vehicles[END_REF]. In the present simulations, a single-phase, incompressible fluid was adopted. Simulations were performed with the segregated approach, which means that momentum equations, pressure-correction equation derived from mass conservation and turbulence quantities transport equations are solved sequentially. The momentum and turbulence equations were discretized using a second-order harmonic TVD scheme [START_REF] Waterson | Design Principles for Bounded Higher-Order Convection Schemes-A Unified Approach[END_REF].

The simulations were performed for both hill-height-based Reynolds numbers of Re 𝐻 = 250, 000 and Re 𝐻 = 650, 000. The computational domain and boundary conditions are chosen according to Fig. 3. At the inlet, Dirichlet boundary conditions are imposed for velocity and turbulence quantities, whereas the pressure follows a Neumann condition. The velocity and turbulence quantities are imposed according to the values given in the inset table of Fig. 3 for the two considered Reynolds numbers. At the outlet, zero normal derivatives for all flow quantities are imposed. The velocities at the outlet are scaled to ensure exact mass conservation between inflow and outflow. At the tunnel walls, a non-permeable no-slip wall condition is applied. A pressure reference point at the inlet is defined with a reference pressure of 𝑝 = 94,450 Pa for Re 𝐻 = 650,000 and 𝑝 = 94,220 Pa for Re 𝐻 = 250,000. Reference values for non-dimensionalization were obtained according to the specifications of Section IV.C.2.

The two-equation SST2003 model and the one-equation SA model were used as turbulence closures for the BeVERLI Hill RANS simulations from the MARIN group.

Computational Grids and Iterative Convergence

The structured grids provided by VT, which are discussed in Section IV.C.4, were used for the ReFRESCO calculations. The computations were conducted with a maximum number of 4000 iterations in steady mode, but were considered to be converged when the 𝐿 2 residuals dropped below 1e-8. For Re H = 250,000 on the coarse Level 7 grid, this convergence criterion was reached, see e.g. Fig. 6. This resulted in maximum changes of the flow quantities of less than 5e-5 for the SST2003 computation on grid Level 7. For the finer grids, stagnation occurred and the maximum number of iterations was reached. An example is shown in Fig. 7, indicating maximum changes of 4e-3 for the SA computation on grid Level 1. For SST2003 on the Level 1 mesh the changes in the velocity were less than 1e-2 and 1e-3 in the pressure.

G. CNRS-Centrale Nantes

Numerical Solver

The in-house solver ISIS-CFD developed by the CNRS-Centrale Nantes, also available as a part of the FINE™/Marine computing suite worldwide distributed by Cadence Design Systems, is an incompressible multiphase unsteady Reynoldsaveraged Navier-Stokes (URANS) solver mainly devoted to marine hydrodynamics. It is based on a fully-unstructured (face-based) finite volume discretization with specific functionalities needed for multiphase flows and industrial applications, see [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF][START_REF] Leroyer | Numerical Methods for RANSE Simulations of a Self-Propelled Fish-Like Body[END_REF]. Volume and surface integrals are evaluated according to second-order accurate approximations. The solver features several sophisticated turbulence models: apart from the classical two-equation 𝑘-𝜖 and 𝑘-𝜔 models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models, are available, see Deng et al. [START_REF] Deng | Comparison of Explicit Algebraic Stress Models and Second-Order Turbulence Closures for Steady Flows around Ships[END_REF], Duvigneau & Visonneau [START_REF] Duvigneau | On the Role Played by Turbulence Closures in Hull Shape Optimization at Model and Full Scale[END_REF] and [START_REF] Cécora | Differential Reynolds-Stress Modeling for Aeronautics[END_REF]. All models are available with wall-function or low-Reynolds near wall formulations. Hybrid RANS/LES turbulence models based on Detached Eddy Simulation (DES-SST, DDES-SST, IDDES) are also implemented and have been thoroughly validated on automotive flows characterized by large separations, see Guilmineau et al. [START_REF] Guilmineau | Assessment of Hybrid RANS-LES Formulations for Flow Simulation around the Ahmed Body[END_REF] and ships at steady drift [START_REF] Visonneau | Local Flow around a Surface Combatant at Various Static Drift Conditions: The Role Played by Turbulence Closures[END_REF]. Moreover, the solver accepts sliding and overset grids and features an anisotropic adaptive grid refinement functionality [START_REF] Wackers | Combined Refinement Criteria for Anisotropic Grid Refinement in Free-Surface Flow Simulation[END_REF][START_REF] Mozaffari | Average-Based Mesh Adaptation for Hybrid RANS/LES Simulation of Complex Flows[END_REF] applied to unstructured hexahedral meshes which will be used in this study.

Computational Grids, Boundary Conditions, and Turbulence Models

In the present study, the ISIS-CFD simulations used the structured grids provided by VT from Level 2 to Level 7. Both hill-height-based Reynolds number of Re 𝐻 = 250,000 and 650,000 were tested. For all meshes, at the inlet, the velocity and turbulence quantities are fixed while at the outlet the pressure is fixed and the other quantities have zero normal derivative. At the tunnel walls, a no-slip wall condition is used.

In the present study a single-phase fluid is used. The scheme used for the discretization of the convective fluxes in both the momentum equations and the equations for turbulence modelling is the AVLSMART [START_REF] Pržulj | Bounded Convection Schemes for Unstructured Grids[END_REF]. The turbulence is modeled with the help of the two-equation SST2003.

Iterative Convergence

All simulations using the VT grids are conducted in a steady mode. In steady mode, a reduction by at least 6 orders of magnitude of the non-linear residuals of the discrete momentum equation is required. Fig. 8 shows some examples of the convergence. With the Level 7 grid (see Fig. 8a) the convergence shows a stagnation, while with the Level 3 grid a reduction by 8 orders is obtained (see Fig. 8b).

V. Results and Discussion

Next, key results available from the BeVERLI Hill RANS simulations at 45 • yaw angle orientation and hill-heightbased Reynolds numbers of Re 𝐻 = 250,000 and 650,000 for a variety of numerical solvers, computational grids, and turbulence models are presented in order of institution. Important sensitivities of the numerical solutions to this vast set of computational configurations, the predictive capability of RANS for the BeVERLI Hill case, and key physical features of the flow over the BeVERLI Hill are highlighted. Lastly, a qualitative cross-comparison of the results is provided to conclude the present section.

A. Virginia Tech

The analysis of the VT computations focuses on the cases at Re 𝐻 = 250, 000 and the SA model. This particular case provides a direct cross-comparison between the SENSEI and the ANSYS Fluent solver computations, since many of the simulations for the higher Reynolds number or the SST2003 turbulence model cases are still ongoing and not completed on at least one of the respective solvers. Similarly, the results are discussed only for the Level 2 and Level 3 grids. Nevertheless, relevant sensitivities of the RANS solution to the choice of solver and key physical features of the flow over the BeVERLI Hill are identified.

Comparing the hill surface 𝐶 𝑝 distribution from the SENSEI and ANSYS Fluent RANS calculations along the centerline (𝑧 = 0 plane) and the centerspan (𝑥 = 0 plane) of the BeVERLI Hill to the VT experimental data (Fig. 9) shows that the numerical solutions closely follow the VT experimental data up to 𝑥/𝐻 < 1, after which the flow undergoes boundary layer separation. Flow separation is a well known weak point in the predictive capability of standard RANS linear eddy-viscosity models [START_REF] Levy | Summary of Data from the Fifth AIAA CFD Drag Prediction Workshop[END_REF][START_REF] Schuster | The Expanding Role of Applications in the Development and Validation of CFD at NASA[END_REF][START_REF] Fremaux | COMSAC: Computational Methods for Stability and Control[END_REF].

Furthermore, a closer look at the centerspan 𝐶 𝑝 distribution (Fig. 9b) reveals a similar picture, where the SENSEI computations very well predict the experimental data across the whole span. Notable differences, however, between the predictions of the SENSEI and ANSYS Fluent computations on the Level 3 grid are observed. The choice of solvers and numerical schemes appears to strongly affect the flow predictions, at least on the Level 3 grid. The SENSEI cases predict a more spanwise symmetric flow, similar to the VT experiment, while the ANSYS Fluent computations predict a spanwise asymmetric flow. Note, the simulations were performed on perfectly spanwise symmetric grids and boundary conditions. This asymmetry effect was observed in an earlier study [START_REF] Gargiulo | Flow Field Features of the BEVERLI Hill Model[END_REF] to occur both in the numerical simulations and the experiment of the BeVERLI Hill for the zero degree rotation angle. The flow over the BeVERLI Hill was hypothesized to be metastable and susceptible to small, destabilizing numerical or experimental asymmetries, such as experimental setup imperfections or asymmetries in the numerical solver scheme employed. This effect becomes more apparent when comparing contours of the 𝐶 𝑝 distribution over the hill surface with overlaid wall shear stress 𝜏 𝑤 lines, as depicted in Fig. 10. A clear spanwise asymmetry affecting the flow on the leeward (𝑥/𝐻 > 0) side of the hill is visible. In particular, the asymmetric flow topology in the separated region (wake) of the hill for the ANSYS Fluent solution (Fig. 10b) jumps to the eye. A pair of counter-rotating vortices at 𝑥/𝐻 = 1 and 𝑧/𝐻 = -2 and 𝑥/𝐻 = 1.5 and 𝑧/𝐻 = 0 is skewed towards the 𝑧/𝐻 < 0 region. In contrast, the flow topology predicted by the SENSEI solution, while not perfectly symmetric in the wake, is overall much more symmetric.

Comparing the 𝐶 𝑝 contours computed for the SENSEI solution on the Level 2 mesh, which was seen to be in good agreement with the experimental data by VT along the centerline and centerspan of the BeVERLI Hill, with a corresponding experimentally obtained contour plot, the quality of the SENSEI result is reinforced (see Fig. 11). The experimental contour was reconstructed through linear interpolation of the pressure data at the corresponding locations of the measurement, the latter which are depicted as black dots in Fig. 11b. It can also be noted that the pressure contour in the wake of the BeVERLI Hill, where the flow is separated, is poorly predicted, consistent with the previous observations.

B. University of Melbourne

As a first step, mesh-dependence is investigated for the RANS baseline case. Considering the mean streamwise velocity (Fig. 12) and the Reynolds stress profiles (Fig. 13 and14) for the baseline RANS runs on the Level 0 and Level 2 grid, it can be observed that for 𝑥/𝐻 ≤ 1 mesh-dependence is negligible, while in regions with 𝑥/𝐻 > 1 some differences between the solutions obtained with the two distinct grids are visible. Thus, further investigation of the impact of the computational grids on the RANS predictions is needed, and grids produced by other participants of this collaborative effort will be tested in future work.

Given the differences observed between the Level 0 and Level 2 grid for the baseline RANS model calculations, the evaluation of the GEP-developed models was conducted only on the finer Level 0 grid. When comparing with the VT experimental data, the GEP RANS improves the 𝐶 𝑝 prediction over the hill for some regions over the baseline RANS, in particular when comparing the 𝐶 𝑝 distribution along the centerspan (𝑥 = 0 plane) of the hill to the experimental data (Fig. 15b). However, the 𝐶 𝑝 distribution along the hill centerline (𝑧 = 0 plane) in the streamwise direction (Fig. 15a) does show an undeprediction of 𝐶 𝑝 on the hill itself and, as a consequence, a recovery downstream of the hill that overshoots the experimental values. This behavior can be explained by inspecting profiles of the GEP RANS predictions (Fig. 16), that show a significant over-prediction of the boundary layer thickness further upstream of the hill geometry, which is caused by a significant increase in the amount of the turbulent kinetic energy (TKE) production and its wall-normal range of impact. This significantly impacts the results, in that the flow over the hill does not separate due to the changed ratio of boundary layer thickness to the hill height. The observed overproduction of TKE has also been previously observed for RANS predictions of flow over a square-mounted cylinder [START_REF] Schmelzer | Discovery of Algebraic Reynolds-stress Models using Sparse Symbolic Regression[END_REF], but given the fixed separation point at the cylinder corners, the consequence of exaggerated TKE levels was not as severe. Thus, further research is needed on the data-driven model development in order to capture the physics pertinent to the flow physics of the current BeVERLI Hill case. Furthermore, notable differences in the prediction of the overall flow mechanisms between experiment and computations are observed, which is well visible in the 𝐶 𝑝 contour plots of Fig. 17. A spanwise asymmetric wake region is predicted by the SST2003 model for a computation with perfectly symmetric and nominal boundary conditions (Fig. 17a). In contrast, the experiment reveals a more symmetric 𝐶 𝑝 distribution (Fig. 17c), which better agrees with the GEP prediction (Fig. 17b). However, the presence of a local pressure maximum on the leeward side of the hill in the GEP contour indicates the prediction of fully attached flow. The experiment, on the other hand, shows that the flow is separated in the leeward region of the hill, indicated by the severe pressure drop in the corresponding 𝐶 𝑝 contour in that region. The GEP algorithm clearly produces an improvement in the prediction of the flow development up to separation as compared to the SST2003 model but, as stated above, the model must be further developed in future to capture important physical aspects of the BeVERLI Hill case.

C. University of New Brunswick and DRDC

The mesh-dependence is first investigated for the two turbulence models tested, the SST2003 and the BSLRSM model. Considering the 𝐶 𝑝 distribution along the centerline of the BeVERLI Hill, it can be noted that mesh-dependence is negligible for both turbulence models and all mesh levels on the windward side (𝑥/𝐻 < 0) of the hill (see Fig. 18a and18b). Over the hill top and on the leeward side (𝑥/𝐻 ≥ 0), differences are still visible between the two finest mesh levels, Level 0 and Level 1. Monotonic convergence towards the finest mesh solution is observed for both models in the fully attached flow regions of the hill (𝑥/𝐻 < 1). In the separated flow region (𝑥/𝐻 ≥ 1) this is not the case. Hereafter, given the above observations, the results will be evaluated and presented only on the finest Level 0 grid. When comparing the 𝐶 𝑝 distribution along the centerline and the centerspan of the BeVERLI Hill for the simulations on the Level 0 grid with the VT experiment (Fig. 19a,19b) the BSLRSM shows clear improvements as compared to the SST2003 model. This is particularly visible in the centerspan 𝐶 𝑝 distribution (Fig. 19b). Here, the local pressure minima in the centerspan 𝐶 𝑝 distribution are better predicted by the BSLRSM model. However, the general shape of the distributions appears to be better represented by the SST2003 model. Notably, both closure models predict a spanwise asymmetric flow despite the perfectly symmetric computational setup. The experimental data reflects the asymmetric behavior, suggesting a true physical nature of this effect. Following the interpretation provided in Section V.A, it can be noted that the susceptibility of the flow over the BeVERLI Hill to destabilizing asymmetries is well pronounced both in the simulations and the experiment. Comparing profiles of the mean streamwise velocity, 𝑈, along the domain centerline (Fig. 20) for the SST2003 and BSLRSM model, the predictions look consistent on the windward side (𝑥/𝐻 ≤ 1) but notable differences in the predictions on the leeward side (𝑥/𝐻 > 1), where the flow separates from the hill surface, are observed.

A qualitative comparison of the 𝐶 𝑝 contours and surface shear-stress-lines, 𝜏 𝑤 , from the simulations with the pressure data and fluorescent oil flow visualization measurement from VT, respectively, is represented in Fig. 21. Qualitatively, the 𝐶 𝑝 contours predicted by the SST2003 model more closely follow the corresponding experimental picture compared to the BSLRSM model. The flow topology represented by the surface shear-stress-lines is qualitatively consistent with the oil flow image for both closure models. The asymmetry of the flow in the wake of the BeVERLI Hill is correctly predicted by both closure models. However, from the experimental flow visualization image it is not possible to draw detailed conclusions about the wake topology. The simulations both predicted a pair of counter-rotating vortices skewed towards 𝑧/𝐻 < 0. One vortex is also clearly visible in the oil flow image. A more quantitative comparison of the 𝐶 𝑝 contours and surface shear stress-line against improved calibrated oil flow images will be an item of future studies on the BeVERLI Hill case.

D. Maritime Research Institute Netherlands -MARIN

Mesh-dependence for all tested closure models and Re 𝐻 is investigated in Fig. 22 in terms of the 𝐶 𝑝 distribution along the centerline of the BeVERLI Hill. In all cases, the solutions at 𝑥/𝐻 < 1 appear to monotonically converge to the finest Level 1 mesh solutions. In the aft section of the hill (𝑥/𝐻 > 1) substantial differences are observed between all mesh levels. Notably, the Level 1 mesh solution of the SA computations at Re 𝐻 = 250, 000 (Fig. 22a), shows the largest change, due to a transition from asymmetric solutions on the coarse grids to an almost symmetric one at the finest grid. Furthermore, the discretization uncertainty in the solutions was estimated. Double precision numbers were used in the computations, which means that round-off errors were negligible. The solutions were interpolated onto the hill centerline spanning -4.5 < 𝑥/𝐻 < 4.5. For each point on this line, the spatial discretization uncertainty 𝑈 𝐶 𝑝 for the pressure coefficent was determined, using the procedure from Eça and Hoekstra [START_REF] Eça | A Procedure for the Estimation of the Numerical Uncertainty of CFD Calculations based on Grid Refinement Studies[END_REF]. Fig. 23 shows the discretization uncertainty for both turbulence models and both Reynolds numbers. A net influence of the turbulence model is captured. On the windward side of the hill, the discretization uncertainty of SA is lower than SST2003 for the Re 𝐻 = 250, 000 case. At Re 𝐻 = 650, 000 this trend is strongly reversed in favor of the SST2003 model. Until 𝑥/𝐻 < 1, it appears that the uncertainty 𝑈 𝐶 𝑝 in the SA prediction are generally lower than those of the SST2003 prediction. After the hill, both turbulence models show a high uncertainty, due to the scatter in the results on the different grid levels: especially for SST2003, the two finest grid solutions show reasonable symmetry in the flow, while the three coarsest ones show large asymmetry to different sides. This is exemplified for the grid Levels 1 through 3 in terms of the surface 𝐶 𝑝 contours with overlaid shear stress lines in Fig. 24. The SA results show similar trends, but the asymmetry in the coarsest grids is less pronounced. Hereafter, given the above considerations, the results will be investigated only on the finest Level 1 grid. Qualitatively comparing the 𝐶 𝑝 contours on the finest grid level from the ReFRESCO computations with the VT experimental data (Fig. 25), we note that the largest deviations between CFD and experiment occur at Re 𝐻 = 650, 000. In this case, the simulations predict an almost perfectly symmetric flow over the BeVERLI hill for both turbulence models, which is in contrast with the asymmetric result obtained in the VT experiment. Although, some slight asymmetries are visible in the SA result when inspecting the flow on the leeward side of the hill. The simulation results at Re 𝐻 = 250, 000 appear to be more consistent with the experimental picture. The SA model predicts some slight asymmetries in the separated wake of the hill, whereas the SST2003 result indicates a more symmetric topology. The VT experiment clearly indicates a sensitivity of the flow over the BeVERLI hill to Reynolds number. In particular, the susceptibility of the flow to destabilizing asymmetries appears to increase with Reynolds number. The numerical solutions of the ReFRESCO simulations seem to be less susceptible to possible numerical asymmetries at both Reynolds numbers on the finest mesh. However, as previously observed, this result varies with decreasing mesh resolution. The same conclusions can also be drawn when considering a detailed view of the 𝐶 𝑝 distribution along the centerline and the centerspan of the hill (Fig. 26). In particular, the slight spanwise asymmetries in the SA results at both Re 𝐻 = 250, 000 and 650,000 become more apparent when looking at the spanwise distribution.

E. CNRS-Centrale Nantes

The ISIS-CFD computations with the SST2003 model are first checked for mesh-dependence. The 𝐶 𝑝 distribution along the hill centerline and centerspan are plotted in Fig. 27 for the Level 2 through 7 grids. Note, the level 2 result is not available for the Re 𝐻 = 650, 000 case. The centerline solution indicates that negligible mesh-dependence was achieved everywhere along the hill centerline for the Level 2 and 3 grids at both Reynolds numbers. The spanwise distribution shows a similar result. However, differences between the Level 2 and Level 3 mesh are still visible in the edge regions (𝑧/𝐻 ≈ ±1) of the hill.

Comparing the 𝐶 𝑝 distribution along the centerline and centerspan for the ISIS-CFD SST2003 computations to the VT experiment (Fig. 27), we note that the computations of the centerline distribution reasonably follow the experimental data for 𝑥/𝐻 < 0 at all grid levels and Reynolds numbers. Over the flat hill top, however, the simulations are seen to deviate from the experimental results. In this region, the flow is subjected to a series of rapid sign changes in the local pressure gradient from favorable (FPG) to adverse (APG), which typically gives rise to non-equilibrium boundary layer effects. The deviations over the hill top then largely affect the subsequent downstream behavior of the flow, where both the separation and recovery regions show underpredicted pressure levels. The spanwise 𝐶 𝑝 distribution reinforces the above observations but further reveals important mesh-dependent solution sensitivities. There is a fundamental mismatch between CFD and experiment in the qualitatively predicted physical behavior of the flow at Re 𝐻 = 250, 000. The finer, more converged meshes (Level 2 and 3) here seem to predict an asymmetric solution, whereas the experiment predicts a more symmetric picture. Interestingly, however, the coarsest Level 7 mesh better follows the symmetric experimental result. This is likely due to the stronger numerical dissipation of coarser meshes, which render the solution less susceptible to destabilizing asymmetries. In contrast, at Re 𝐻 = 650, 000 also the experiment indicates the presence of asymmetric flow, which is in agreement with the simulations.

A broader view at the overall flow topology in terms of the 𝐶 𝑝 contours and surface shear stress lines (see Fig. 28) reveals a complex wake behavior. The simulations predict two pairs of counter rotating vortices for both Reynolds numbers. At Re 𝐻 = 250, 000 the first pair is visible at around 𝑧/𝐻 = 0 and 𝑥/𝐻 = 1, the second pair consists of a vortex at about 𝑧/𝐻 = 2 and 𝑥/𝐻 = 0.75 and 𝑧/𝐻 = -2 and 𝑥/𝐻 = 1.5. Each of these vortex pairs is skewed asymmetrically.

F. Qualitative Cross-Comparison

A few points are worth noting when cross-comparing the above results across institutions (see Fig. 29 and30). First, the turbulent flow over the BeVERLI Hill, as seen in the experimental data, is found to be Reynolds number sensitive. At Re 𝐻 = 250, 000 the flow is found to be spanwise symmetric, in agreement with our physical intuition. At Re 𝐻 = 650, 000 the flow shows strong spanwise asymmetries from 𝑥/𝐻 >= 0 onwards. This is despite nominally symmetric boundary conditions in the experiment. A possible interpretation of this effect was provided in an earlier study by Gargiulo et al. [START_REF] Gargiulo | Flow Field Features of the BEVERLI Hill Model[END_REF]. The flow over the BeVERLI Hill was hypothesized to be inherently metastable and susceptible to small, destabilizing asymmetries. These asymmetries can originate from imperfection in the experimental setup and the corresponding boundary conditions. The asymmetric state would represent a more energetically advantageous state. The simulations also show that the particular topology of the wake asymmetry is sensitive to the computational setup. The number of nodal and saddle points contained in the separated wake are seen to vary across different solvers and meshes. This suggest that the level of numerical dissipation pertinent to a particular numerical scheme or mesh resolution notably impacts the results. This is particularly visible in the simulations at the higher Reynolds number case, which appear to be more sensitive to this effect.

Overall, the computational results and their scatter across institutions reveal relevant key features of the BeVERLI Hill that pose challenges to the predictive capability of current RANS models. This becomes apparent when inspecting Fig. 31. Note, in the absence of a rigorous validation of the results, including the estimation of numerical, modeling, and experimental errors and uncertainties, conclusions here can only be drawn qualitatively. The computational results seem to fairly consistently predict the attached flow over the BeVERLI Hill up to around 𝑥/𝐻 < 0. Although, the prediction of the highest and lowest pressure points seems to be sensitive to the particular choice of numerical discretization. A larger scatter of numerical results across institutions and deviations from the experimental data are observed over the flat top of the BeVERLI Hill. The repeated and rapid change in the sign of the pressure gradient in that region gives rise to complex non-equilibrium flow physics. The results in that region were seen to be particularly sensitive to the choice of turbulence model, hinting at a potential weakness and limitation of the tested turbulence closures. Discrepancies from the experimental data seem to then be amplified in the separated wake of the hill, knowingly a weakness of RANS models [START_REF] Levy | Summary of Data from the Fifth AIAA CFD Drag Prediction Workshop[END_REF][START_REF] Schuster | The Expanding Role of Applications in the Development and Validation of CFD at NASA[END_REF][START_REF] Fremaux | COMSAC: Computational Methods for Stability and Control[END_REF]. The three-dimensional curvature of the BeVERLI Hill geometry is likely responsible for the spanwise sensitivity of the flow, which in turn leads to the observed wake asymmetry.

All of these challenges render the case of the BeVERLI Hill an attractive case for the study of non-equilibrium flow physics. However, the susceptibility of the BeVERLI Hill case towards numerical asymmetries might arguably pose challenges to the appropriate validation of CFD. The present study is yet inconclusive on whether the asymmetry present or absent in the experiment is correctly caught by the simulations. The strong dependence of the asymmetry prediction on the grid density, iterative convergence level and discretization scheme might render it difficult to converge to the experimental solution. Additionally, it is important to ensure that the experimental solution is unambiguous and repeatable, and not prone to the same severe sensitivities recorded in the computations. To investigate the latter point, an ongoing effort is currently dedicated at reproducing and quantitatively cross-comparing the VT experiment across multiple experimental facilities around the globe. In [START_REF] Duetsch-Patel | The BeVERLI Hill Three-Dimensional Separating Flow Case: Cross-Facility Comparisons of Validation Experiment Results[END_REF] a qualitative comparison of the experimental solution has been completed and found to be consistent across facilities. Future work will be also dedicated to investigations with unsteady RANS (URANS), LES, and hybrid methods, which might open new avenues and interesting pathways for the BeVERLI Hill validation case.

VI. Conclusions

The BeVERLI Hill case has become an active subject of the NATO AVT-349 group's research activities. This international collaborative effort allows for extensive comparisons of RANS solutions of the BeVERLI Hill across all participating institutions and an array of CFD solvers, computational grids, and turbulence models, including RANS closures developed using data-driven approaches [START_REF] Weatheritt | A Novel Evolutionary Algorithm Applied to Algebraic Modifications of the RANS Stress-Strain Relationship[END_REF]. The present work presented key results available to date summarizing and highlighting the characteristic features and sensitivities as well as the current predictive capability of RANS for the complex flow of the BeVERLI Hill.

Six distinct solvers, three sets of computational grids, and four different turbulence closure models were employed to simulate the turbulent flow over the BeVERLI Hill at 45 • yaw angle orientation and Re 𝐻 = 250, 000 and 650,000. The results indicate important sensitivities to Reynolds number, grid density, solver and numerical scheme, and turbulence models. In particular, the experimentally observed wake asymmetry of the BeVERLI Hill was demonstrated to be also triggered purely by numerical effects, such as grid density, iterative convergence, or the particular choice of numerical schemes. The predictions of RANS were seen to be consistent in the fully attached flow regime over the BeVERLI Hill up to 𝑥/𝐻 < 0. Deviations from the experiment and scatter of results were correlated with the flat top and separated flow regions of the BeVERLI Hill. Three-dimensional curvature effects add to the complexity of this case.

The intricacies and complex physics of the BeVERLI Hill are an attractive case for CFD validation and for the advancement of turbulence models. However, the susceptibility of the BeVERLI Hill case towards numerical asymmetries might arguably render the appropriate use of this case in the validation of CFD models intricate. Pending a rigorous validation of the numerical results it is difficult to conclude whether the asymmetry present or absent in the experimental results is correctly caught by the simulations. The strong dependence of the asymmetry prediction on the grid density, iterative convergence level and discretization scheme might render it difficult to converge to the experimental solution. Additionally, it is to be determined, whether the experimental result is itself exactly reproducible. For this reason, ongoing cross-facility comparisons of the VT experiment across multiple experimental facilities around the globe are being conducted.

Future studies will be dedicated to quantitative and rigorous verification and validation of RANS results across institutions. This will include numerical error and uncertainty estimations and experimental uncertainty quantification. The release of additional experimental data from the VT team will support the improvement of boundary conditions matching between CFD and experiment, including the use of an as-tested hill geometry. Similarly, a final improved set of computational grids will be released and used consistently amongst all collaborators. Lastly, dedicated studies using URANS, LES, and hybrid methods will be conducted to explore alternative pathways for the BeVERLI Hill validation case.

Appendix A

For the present optimization, we choose to use gene expression programming (GEP), a symbolic regression algorithm that returns mathematical equations [START_REF] Ferreira | Gene Expression Programming: A New Adaptive Algorithm for Solving Problems[END_REF]. The specifics of applying GEP to turbulence models was introduced in detail in Weatheritt and Sandberg [START_REF] Weatheritt | A Novel Evolutionary Algorithm Applied to Algebraic Modifications of the RANS Stress-Strain Relationship[END_REF]. A non-linear constitutive stress-strain relationship is constructed using our GEP framework by adding non-linear terms to the linear Boussinesq approximation,

𝑢
Here, the Reynolds stress 𝑢 𝑖 𝑢 𝑗 is related to its trace 2𝑘, strain rate tensor 𝑆 𝑖 𝑗 and additional non-dimensional terms. 𝑇 𝑛 𝑖 𝑗 are basis functions [START_REF] Pope | A more General Effective-Viscosity Hypothesis[END_REF], formed from polynomials of non-dimensional strain 𝑠 𝑖 𝑗 = 𝜏 𝑆 𝑖 𝑗 and rotation rate 𝑤 𝑖 𝑗 = 𝜏 Ω 𝑖 𝑗 , with 𝜏 being the turbulence model timescale. The coefficients 𝜁 𝑛 in equation ( 1) are functions of the non-dimensional invariants I 𝑛 produced from 𝑠 𝑖 𝑗 and 𝑤 𝑖 𝑗 , of which there are five for a three-dimensional flow [START_REF] Pope | A more General Effective-Viscosity Hypothesis[END_REF].

In addition, following the approach of Schmelzer et al. [START_REF] Schmelzer | Discovery of Algebraic Reynolds-stress Models using Sparse Symbolic Regression[END_REF], a second model is developed to correct the production term in 𝑘 -𝜔 turbulence transport equations. The term 𝑅 (a local production term) is modelled as 𝑅 = 2𝑘 R𝑖 𝑗 𝜕 𝑗 𝑈 𝑖 .

The resulting GEP-derived non-linear models, trained using a high-fidelity dataset of a finite square-mounted cylinder with a height-to-width ratio of 4 and Re 𝑑 = 11, 000 (based on the cylinder width 𝑑), are: (3)

𝑎 𝑥,

Fig. 1

 1 Fig. 1 The BeVERLI Hill geometry. Shown are, on the left, a wall-mounted representation of the BeVERLI Hill along with a Cartesian coordinate system, the flow direction, and a schematic of the incoming boundary layer. On the right, a detail view of the BeVERLI Hill geometry, including the superelliptic corners (𝑪 1 -𝑪 4 ), the flat top, and the 5th-degree polynomial centerline profile 𝑷 5 (𝒙).

Fig. 3

 3 Fig. 3 Computational domain for the VT BeVERLI Hill RANS simulations with highlighted boundary conditions. The inlet (blue), outlet (purple), and floor (orange) plane are emphasized using different colors. The reference static pressure ports on the wall opposite to the BeVERLI Hill are represented as green dots. The table on the left summarizes the boundary conditions for the simulated hill-height-based Reynolds numbers.

Fig. 4 .

 4 Other reference quantities, such as the reference Mach number 𝑀 ref , static temperature 𝑇 ref , velocity 𝑉 ref , and density 𝜌 ref are computed from the reference pressure 𝑝 ref and the stagnation conditions 𝑝 0,in and 𝑇 0,in using standard compressible isentropic flow relations.

Fig. 4

 4 Fig. 4 Schematic of the wall-mounted BeVERLI Hill in the VT SWT test section. The Cartesian coordinate system utilized in both the experiment and the simulations is highlighted. Included is also a projected copy of the surface static pressure ports locations used in the original experiment. In particular, the location of the reference pressure ports is emphasized.

Fig. 5

 5 Fig. 5 VT systematically refined grid Levels 1 through 4 with refinement factor 𝒓 = 2 1/3 . The grids are ordered from top to bottom in decreasing level of refinement. The number of finite volume cells for each mesh level is indicated in parenthesis under each grid level label. A Cartesian coordinate system and the flow direction are included.

Fig. 6 Residual

 6 Fig. 6 Residual 𝑳 2 and 𝑳 ∞ norms for the MARIN ReFRESCO computations at Re 𝑯 = 250, 000 on the Level 7 grid and the SST2003 model.

  Fig. 7 Residual 𝑳 2 and 𝑳 ∞ norms for the MARIN ReFRESCO computations at Re 𝑯 = 250, 000 on the Level 1 grids and the SA model.

Fig. 8 Residual

 8 Fig. 8 Residual 𝑳 2 norms for the CNRS-Centrale Nantes ISIS-CFD computations at Re 𝑯 = 250, 000 on the Level 7 and Level 3 grids and the SST2003 model.

Fig. 9

 9 Fig. 9 Comparison of the BeVERLI Hill surface 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) and centerspan (𝒙 = 0 plane) from the ANSYS Fluent and SENSEI RANS SA model computations to the VT experiment at Re 𝑯 = 250, 000. Shown are the numerical results computed on the Level 2 and Level 3 grids. The Cartesian coordinates are normalized by 𝑯. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

Fig. 10

 10 Fig. 10 Comparison of the BeVERLI Hill surface 𝑪 𝒑 distribution contours with overlaid wall shear stress 𝝉 𝒘 lines from the ANSYS Fluent to the SENSEI RANS SA model computations at Re 𝑯 = 250, 000. Shown are the results computed on the Level 3 grid. A solid, white colored line traces the baseline profile of the BeVERLI Hill. The Cartesian coordinates are normalized by 𝑯.

Fig. 11

 11 Fig. 11 Comparison of the BeVERLI Hill surface 𝑪 𝒑 distribution contours with overlaid wall shear stress 𝝉 𝒘 lines from the SENSEI RANS SA model computations on the Level 2 grid with the VT experimental data at Re 𝑯 = 250, 000. A solid, white colored line traces the baseline profile of the BeVERLI Hill. The Cartesian coordinates are normalized by 𝑯.

1 Fig. 12

 112 Fig. 12 Mean streamwise velocity, 𝑼, profiles at different 𝒙-locations in the 𝒛 = 0 plane normalized by 𝑼 ref for the OpenFOAM baseline RANS SST2003 computations at Re 𝑯 = 250, 000. Shown is a comparison of the results on the Level 0 and Level 2 grids. The Cartesian coordinates are normalized by 𝑯.

1 Fig. 13

 113 Fig. 13 Streamwise normal Reynolds stress, 𝒖 2 , profiles at different 𝒙-locations in the 𝒛 = 0 plane normalized by 𝑼 ref for the OpenFOAM baseline RANS SST2003 computations at Re 𝑯 = 250, 000. Shown is a comparison of the results on the Level 0 and Level 2 grids. The Cartesian coordinates are normalized by 𝑯.

1 Fig. 14

 114 Fig. 14 Reynolds shear-stress, 𝒖𝒗, profiles at different 𝒙-locations in the 𝒛 = 0 plane normalized by 𝑼 ref for the OpenFOAM baseline RANS SST2003 computations at Re 𝑯 = 250, 000. Shown is a comparison of the results on the Level 0 and Level 2 grids. The Cartesian coordinates are normalized by 𝑯.

Fig. 15 BeVERLI

 15 Fig. 15 BeVERLI Hill surface 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) and centerspan (𝒙 = 0 plane) from the OpenFOAM RANS computations with the SST2003 and GEP models on the Level 0 grid as compared to the VT experiment at Re 𝑯 = 250, 000. The Cartesian coordinates are normalized by 𝑯. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

Fig. 16

 16 Fig. 16 Mean streamwise velocity, 𝑼, profiles normalized by the reference velocity 𝑼 ref for the incoming boundary layer of the BeVERLI Hill at 45 • yaw angle orientation and Re 𝑯 = 250, 000. Shown is a comparison between the results from the OpenFOAM baseline SST2003 and GEP RANS simulations for the Level 0 mesh. The spatial Cartesian coordinates are normalized by the hill height 𝑯.

Fig. 17

 17 Fig. 17 Comparison of the BeVERLI Hill surface 𝑪 𝒑 distribution contours from the OpenFOAM baseline RANS SST2003 and GEP computations at Re 𝑯 = 250, 000. Shown are the results computed on the Level 0 grid. The Cartesian coordinates are normalized by 𝑯.

Fig. 18 BeVERLI

 18 Fig. 18 BeVERLI Hill surface 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) from the ANSYS CFX RANS computations with the SST2003 and BSLRSM model on the grid levels 0 through 3 at Re 𝑯 = 650, 000. The Cartesian coordinates are normalized by 𝑯. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

Fig. 19 BeVERLI 1 Fig. 20

 19120 Fig. 19 BeVERLI Hill surface 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) and centerspan (𝒙 = 0 plane) from the OpenFOAM RANS computations with the SST2003 and BSLRSM models on the Level 0 grid as compared to the VT experiment at Re 𝑯 = 650, 000. The Cartesian coordinates are normalized by 𝑯. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

  Fig. 21 Comparison of the BeVERLI Hill surface 𝑪 𝒑 distribution contours from the ANSYS CFX RANS SST2003 and BSLRSM computations on the Level 0 mesh to the VT experimental data at Re 𝑯 = 650, 000. The numerical results feature and overlay of wall shear stress lines. The corresponding VT data consist of a measured pressure contour and fluorescent oil flow visualization. The spatial Cartesian coordinates are normalized by the hill height 𝑯. The perimeter of the hill is included in white color for the computational results.

  Fig. 22 Grid-independence study for the MARIN ReFRESCO RANS SA and SST2003 computations in terms of the BeVERLI Hill surface 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) at Re 𝑯 = 250, 000 and 650,000. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

Fig. 24 BeVERLIFig. 25

 2425 Fig. 23 𝑼 𝑪 𝒑 distribution along the centerline of the BeVERLI Hill for ReFRESCO computations at Re 𝑯 = 250, 000 and 650,000 with the SA and SST2003 models. The hill centerline profile is schematically included as black colored, dashed line.

Fig. 26 BeVERLI

 26 Fig. 26 BeVERLI Hill surface 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) and centerspan (𝒙 = 0 plane) from the ReFRESCO RANS computations with the SST2003 and SA models on the Level 1 grid as compared to the VT experiment at Re 𝑯 = 250, 000 and 650,000. The Cartesian coordinates are normalized by 𝑯. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

Fig. 27

 27 Fig. 27 Comparison of 𝑪 𝒑 distribution along the centerline (𝒛 = 0 plane) and centerspane (𝒙 = 0 plane) from ISIS-CFD RANS SST2003 computations on the grid levels 2 through 7 with the VT experiment at Re 𝑯 = 250, 000 and 650,000. A profile of the hill with the height scaled by a factor 𝒏 (i.e. 𝑯/𝒏) is included.

Fig. 28

 28 Fig. 28 Comparison of the BeVERLI Hill surface 𝑪 𝒑 distribution contours from the ISIS-CFD RANS SST2003 computations on the Level 2 and 3 mesh to the VT experimental data at Re 𝑯 = 250, 000 and Re 𝑯 = 650, 000, respectively. The numerical results feature and overlay of wall shear stress lines. The spatial Cartesian coordinates are normalized by the hill height 𝑯.

Fig. 30

 30 Fig. 30 Cross-comparison of BeVERLI Hill surface 𝑪 𝒑 distribution contours with overlaid wall shear stress, 𝝉 𝒘 , lines on the finest available mesh across institutions with the VT experiment at Re 𝑯 = 650, 000.

( a )Fig. 31

 a31 Fig. 31 Cross-comparison of BeVERLI Hill surface 𝑪 𝒑 distribution along the hill centerline and centerspan on the finest available mesh across institutions with the VT experiment at Re 𝑯 = 250, 000 and 650,000.

Table 1 Cartesian coordinates of reference static pressure ports. x [m] y [m] z [m]

 1

	-2.228 1.83 -0.6858
	-2.228 1.83 -0.4572
	-2.228 1.83 -0.2286
	-2.228 1.83	0
	-2.228 1.83	0.2286
	-2.228 1.83	0.4572
	-2.228 1.83	0.6858

  𝑖 𝑢 𝑗 = 2 3𝑘𝛿 𝑖 𝑗 -2𝜈 𝑡 𝑆 𝑖 𝑗 linear Bouss.

	+	2𝑘	10 ∑	𝜁	𝑔𝑒 𝑝 𝑛 𝑇 𝑛 𝑖 𝑗 .
			𝑛=1

𝑎 𝑥

𝑖 𝑗 :additional non-linear terms

  𝑔𝑒 𝑝 𝑖 𝑗= 𝑇1 𝑖 𝑗 (8.0𝐼 1 -0.399)+ 𝑇 2 𝑖 𝑗 (2.0𝐼 1 -3.43)-𝑇 4 𝑖 𝑗 (2.0𝐼 2 -𝐼 5 + 2.0)+ 𝑇 8 𝑖 𝑗 (𝐼 1 -𝐼 2 -𝐼 3 + 6.15)-𝑇 9 𝑖 𝑗 (𝐼 2 + 𝐼 3 -(𝐼 3 -0.15) (𝐼 5 -0.15) -1.178). 𝑖 𝑗 (4.0𝐼 2 -2𝐼 1 + 2.15)+ 𝑇 2 𝑖 𝑗 (𝐼 5 -𝐼 1 (1.911𝐼 3 -𝐼 1 + 𝐼 5 ) + 1.0)-𝑇 4 𝑖 𝑗 (2.0𝐼 2 -12.0)+ 𝑇 8 𝑖 𝑗 (𝐼 1 )+ 𝑇 9 𝑖 𝑗 (6𝐼 2 + 𝐼 1 𝐼 5 + 18.089).

	R𝑔𝑒𝑝 𝑖 𝑗	= 𝑇 1
		(2)

√2 to systematically coarsen and obtain each subsequent grid level down to the coarsest Level 7 grid. 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 order of magnitude from the initial levels.

Acknowledgments

The Virginia Tech team thanks NASA, in particular Dr. Michael Kegerise, Dr. Mujeeb Malik, and Dr. Christopher Rumsey, for their support under grant 80NSSC18M0146. We are extremely grateful for the contributions of Vidya Vishwanathan and Jarrod Banks to wind tunnel testing and the acquisition of experimental data. We would like to thank the Advanced Resource Computing division for the usage of HPC resources. We also thank the Virginia Tech Aerospace and Ocean Engineering Machine Shop, and particularly James Lambert and Cameron Hollandsworth, for their support in designing and manufacturing test hardware and instrumentation.

The University of Melbourne team acknowledges support from the U.S. Office of Naval Research (ONR) under NICOP Grant N62909-20-1-2046 with program monitors Dr. Ki-Han Kim (ONR), and Dr. Sung-Eun Kim (ONR Global, Tokyo).

The CNRS -Centrale Nantes team gratefully acknowledges GENCI (Grand Equipement National de Calcul Intensif) (Grant-A0082A00129, Grant-A0102A00129) for the use of HPC resources.

The research by MARIN is partly funded by the Dutch Ministry of Economic Affairs.

The computations of the different participating institution, however, have demonstrated that the asymmetry effect is also present in the numerical solutions and, therefore, triggered by numerical asymmetries. In particular, the VT, University of Melbourne, MARIN, and CNRS-Centrale Nantes computations show that an asymmetric solution can also be triggered at the lower Reynolds number case, contrasting the experimental result by VT. Note, only nominal and perfectly spanwise symmetric boundary conditions and grids were used in the computations. By exclusion, the numerical asymmetries which destabilize the flow must originate from the particular solver discretization and solution schemes utilized and iterative convergence achieved. However, the results clearly show that the numerical solution can be more or less susceptible to these numerical asymmetries depending on the numerical grid or the choice of turbulence model. For instance, the SST2003 model appears to be less sensitive and more often predicts a symmetric flow in tune with physical intuition.