Structural aspects of the attached turbulent boundary layer flow over a hill

Three-dimensional turbulent boundary layers under strong pressure gradients and curvature are the rule in real-world flow applications but are typically not well predicted by turbulence models due to isotropic eddy viscosity or equilibrium assumptions. Validation-quality data in complex 3D flows is necessary for continued efforts to improve simulation accuracy. The Benchmark Validation Experiments for RANS/LES Investigations (BeVERLI) Hill bump model, designed specifically for validation experiments, was tested in the Virginia Tech Stability Wind Tunnel to collect validation experiment data on the three-dimensional (3D) boundary layer flow over a 3D hill. Laser Doppler velocimetry measurements on the bump model were used to study the mean flow and turbulence structure and evaluate the impact of pressure gradient and curvature upon the total shear stress in the boundary layer and evaluate the impact of pressure gradient and curvature upon the total shear stress behavior in the near-wall region. From analysis of the BeVERLI Hill flow, including the boundary layer just upstream of the hill, and comparison with the 3D flow around a wing-body junction of Ölçmen et al. (Exp Fluids 31:219–228, 2001), it is shown that none of the stations studied exhibit a constant shear stress region over any significant region of the boundary layer.


Introduction
Three-dimensional turbulent boundary layers (3D TBLs) are of extensive interest for practical engineering applications, but have been less frequently studied than two-dimensional turbulent boundary layers (2D TBLs). Due to the challenges that three-dimensional boundary layers present to accurate computational modeling, their continued study is relevant today. Some well-known features of 3DTBLs include the lag between the shear-stress direction vector and the velocitygradient vector direction, a reduction in near-wall turbulent kinetic energy and Reynolds stresses, and a reduction in the turbulence parameter A 1 = √ 3 ) below its typical 2D value of 0.15, to as low as 0.03 to 0.05 in the case of Anderson and Eaton (1989). Because of this, 3D TBLs present a particular challenge to modeling, as the stress/strain misalignment invalidates the assumptions behind any turbulence models that assume an isotropic eddy viscosity. This limits the utility of some popular Reynolds-Averaged Navier-Stokes (RANS) computational fluid dynamics (CFD) turbulence models for flows of engineering interest, even as CFD is increasingly used for high-impact decisions in vehicle design. Improved modeling of 3D flows depends upon a better understanding of 3D boundary layers and turbulence to improve the assumptions that go into turbulence models.
Unlike RANS, large-eddy simulation (LES) resolves the larger turbulence scales and models the small scales that cannot be resolved on the grid. As a result, LES requires a sub-grid scale (SGS) model for the effects of the small scales on the resolved flow. Per Bose and Park (2018), the total number of grid points required for wall-resolved and wall-modeled LES are N WR ≈ Re respectively. Because of this, wall-resolved LES is infeasible for higher Reynolds number applications, requiring wallmodeled LES. Like RANS turbulence models, existing wall models are largely built on equilibrium assumptions, which are not valid in the vast majority of applicable 3D flows. These assumptions can be built into the model directly, for example through assuming an equilibrium stress or velocity distribution, or indirectly through assumptions regarding eddy viscosity and inner-layer Reynolds stress (Bose and Park 2018).
Many different specific definitions for "equilibrium" have been used in literature. In this work, equilibrium describes an idealized state in which all flow quantities are self-similar (in which the profiles of a particular quantity collapse to a single curve with the proper normalization) for a given set of consistent scaling parameters. In this state, the normalized flow quantities are no longer a function of streamwise position in the flow. Equilibrium two-dimensional boundary layers typically have a region near the wall where the total shear stress is nearly constant, approximately 0.95 ≤ total ∕ w ≤ 1 for sufficiently high Reynolds numbers. Devenport and Lowe (2022) demonstrate the derivation of 2D zero-pressure gradient boundary layer equations that produce this behavior, showing that as the viscous shear stress decreases rapidly outside the viscous sublayer, the turbulent shear stress increases to balance out the total shear stress. In nonzero pressure gradient flows, a linear gradient of the total shear stress will be present instead.
Equilibrium approximations to the thin-boundary layer equations are often made to relax strict assumptions on the velocity profile that would otherwise be required for LES (Bose and Park 2018). These approximations also lead to the assumption of a constant total stress region in the near-wall region of the boundary layer. However, most flows of engineering interest involve strong pressure gradients in both the streamwise and lateral direction and other non-equilibrium phenomena, which will invalidate these assumptions.
The flow over a wall-mounted bump is a common test case to the goal of improved computational fidelity (e.g. Byun and Simpson 2006;Bell et al. 2012), as the simple geometry generates a flow field with pressure gradient, curvature effects, and three-dimensional boundary layer separation and reattachment. All these features are generally challenging for CFD and are very common in practical applications. This work presents results from experiments studying the three-dimensional turbulent attached and separated flow over the BeVERLI (Benchmark Validation Experiment for RANS/LES Investigations) Hill. Previous publications of this flow case include Gargiulo et al. (2020), Gargiulo et al. (2021), Gargiulo et al. (2022), andDuetsch-Patel et al. (2022), among others. The flow over this model experiences three-dimensional pressure gradients of varying sign and magnitude, skewing, three-dimensional curvature, and three-dimensional flow effects.
We study the evolution of the boundary layer from upstream of the hill to over the hill itself, specifically the impact of curvature, 3D pressure gradients, and skewing upon the total shear stress distribution in the boundary layer. Results are also compared RANS simulations and with the 3D TBL flow over a flat wall around a wing-body junction of Ölçmen et al. (2001), which also experiences three-dimensional pressure gradients and flow skewing. It will be shown that streamwise and lateral pressure gradients and skewing significantly impact the flow, reducing the turbulent shear stress throughout the boundary layer and resulting in no constant shear stress region for any of the stations studied. Any model that relies on equilibrium shear stress assumptions will not be able to accurately predict turbulence in threedimensional turbulent boundary layers.

Stability wind tunnel
Experiments were conducted in the stability wind tunnel (SWT) at Virginia Tech, a continuous, single-return subsonic wind tunnel, with a 1.85 × 1.85 × 7.32 m test section. A schematic of the facility is shown in Fig. 1. The SWT is capable of reaching a maximum flow speed of approximately 80 m/s, corresponding to a Reynolds number per meter of 5 × 10 6 . The facility features two interchangeable test sections, one of which has highly modular side walls to facilitate model and instrumentation installation , as well as an instrumented contraction liner to track the surface static pressure through the contraction . Measurements over the BeV-ERLI Hill were performed in a configuration as depicted in Fig. 2.
Extensive characterization has been conducted for the facility, including measurements of boundary layer and cross-sectional flow variations. An example of this characterization is measurements of the inflow boundary layer made using a boundary layer rake consisting of 30 Pitot probes to measure total pressure and an additional two static ports on the wall for static pressure. Measurements obtained for two different experimental builds in February 2020 and May 2021 were shown to be in good agreement (Duetsch-Patel et al. 2022). The inflow boundary layer profile at three points across the span of the inflow is shown in Fig. 3, and the inflow boundary layer parameters are shown in Table 1. As shown in Table 1, measurable spanwise non-uniformity is exhibited in the boundary layer parameters. This non-uniformity in the inflow boundary layer is a critical boundary condition to be included in CFD simulations of this validation case.

BeVERLI hill model
The bump geometry is defined in detail by Gargiulo et al. (2020) and has a design width w of 0.93472 m and height H of 0.1869 m. A schematic of the hill geometry is shown in Fig. 4. The 2D centerline and centerspan cross-section of the 3D geometry is defined by a fifth-degree polynomial mirrored about x 1 = 0 , with a flat top region in the center connecting the two polynomials to form a cruciform. These profiles form cylindrical regions s = 0.0935 m wide. The edges of these perpendicular cylindrical sections are connected at a constant height using a fourth order superellipse. Because of these features, the model exhibits 90 degree rotational symmetry, allowing for different orientations of the model with respect to the incoming flow. For these experiments, the bump was oriented as shown in Fig. 2 with one of the pairs of rounded corners directly aligned with the approaching freestream. This orientation is referred to as the 45 • yaw case. Coordinate systems labeled X i indicate global coordinates, with the origin at the center of the hill as shown in Fig. 2, while coordinate systems labeled lowercase x i indicate local wall-shear stress coordinates, with the origin at the wall-location of the profile. The hill was mounted on  the wind tunnel side wall to facilitate access to embedded instrumentation. Two CNC-milled models were used: One with 135 pressure taps installed in the surface and one with slots for clear windows for laser Doppler velocimetry (LDV) measurements along the fifth-degree polynomial curves. 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 and a measurement accuracy of approximate 8 μ m for the size of the scanned volume. These scans found that the manufactured geometry matched the design geometry to within ±0.40 mm and ±0.30 mm, respectively.
The hill was installed 6.88 m downstream of the 3.18 mm boundary layer trip in the wind tunnel contraction, and 1.83 m downstream of where 2D inflow measurements were captured. Measurements were made with the freestream flow adjusted to ensure hill height-based Reynolds numbers Re H = U ∞ H of 250,000 ± 5000 , 325,000 ± 5000 , and 650,000 ± 5000 . Measurements collected at Re H = 250, 000 are the focus of this study. The ratio of the bump height relative to the incoming boundary layer at the inflow along the centerline, H∕ inflow , is 3.3.

Instrumentation
Pressure measurements on the model surface were collected using 1.6 mm (0.063 in) diameter embedded steel tubes, which were connected via tubing to pressure scanners. Ninety-five taps were connected to three DTC ESP 32 HD scanners with a range of 2.5 psi and a manufacturerrated accuracy of 0.06% full-scale. The remaining taps were connected to an Esterline 9816/98RK pressure scanner with range of 2.5 psi and a rated accuracy of 0.05% full scale. The locations of these taps on the hill are indicated in Fig. 5. Measurements at four symmetric yaw angles were used to compute an averaged surface pressure contour. Comprehensive uncertainty estimations in the style of Aeschliman and Oberkampf (1998) and Rhode and Oberkampf (2017) found an uncertainty in static pressure coefficient measurements on the bump, including random uncertainty, flowfield nonuniformity uncertainty, and instrumentation and model geometry uncertainty, to be ±0.025 at Re H = 250,000 at 95% confidence.
A specialized embedded three-component LDV system, originally developed by AUR, Inc. (Lowe et al. 2014), was used to take measurements of the turbulent boundary layer over the bump surface and upstream of the model. The probe was mounted inside the bump, with five laser beams passing through the 1.5 mm anti-reflective curvature-fitting acrylic windows, as shown in Fig. 6. Three beams are generated with a continuous wave, frequency-doubled, diode pumped As-designed geometry of the BeVERLI Hill. Highlighted blue regions indicate the cylindrical regions of the hill, defined by the fifth-degree polynomial geometry. A flat top (highlighted in red) connects the cruciform formed by these four cylindrical fifth-degree polynomial regions, and the corners between these cruciforms are connected by superelliptic curves. The hill exhibits 90 • rotational symmetry ND:YVO4 solid state Coherent Verdi V6 laser with a wavelength of 532 nm, while two beams are generated with a continuous wave, single-line mode, optically pumped semiconductor Coherent Genesis MX laser with a wavelength of 488 nm. The collimating lenses and receiving fiber were replaced from the original configuration to reduce the measurement volume size and allow for measurements nearer to the surface.
Three measurement volumes were created from the intersection of the five beams to obtain three components of velocity. IntraAction Corporation acousto-optic modulators were used to frequency shift one beam out of each pair forming a measurement volume to generate a moving interference fringe pattern. This moving fringe pattern recovers directional information on the movement of the particle through the measurement volume that would be lost if a static fringe pattern were used. Two of the three green Verdi V6 beams were frequency-shifted by + 40 MHz and − 80 MHz, respectively, while one of the two Genesis MX beams was frequency-shifted by + 60 MHz. The Verdi V6 was operated at approximately 0.8-1.2 W and the Genesis MX was operated at approximately 0.25 to 0.50 W during testing. The laser beam power emitted from the probe was approximately 70-100 mW after passing through table optics and fibers.
Velocity components are measured perpendicular to the local interference fringes generated by the intersection of a given pair of laser beams, as shown in Fig. 7. Each component is a function of the detected Doppler frequency ( f D ) generated by the light scattered by the particle as it passes through the fringes and the fringe spacing ( Δx ), as shown in Eq. 1. The spacing of interference fringes created by the intersection of two laser beams is a function of the beam wavelength ( b ) and the angle between the beams ( ).
These velocity components are in a non-orthogonal coordinate system due to the non-orthogonal orientation of the fringe and measurement volume intersections, and thus must be converted to orthogonal coordinates. The vector in the fringe-perpendicular direction can be computed by measuring the beam direction vectors, � ⃗ v b1 and � ⃗ v b2 for a given pair. From the view of Fig. 7, the cross-product of these two vectors will yield a vector pointing out of the page. The cross-product of this vector with the bisector of � ⃗ v b1 and � ⃗ v b2 , Schematic of the LDV probe embedded in the BeVERLI Hill. The probe was traversed along the plane of the cylindrical section of the hills from station to station, as shown in a. The probe was rotated at each station as shown in b to ensure that the probe was oriented normal to the surface. The local coordinate x 2 at each station is normal to the local surface curvature Fig. 7 Intersection of two beams to form one measurement volume, measuring one component of velocity. Δx indicates the spacing between interference fringes generated by the intersection, and � ⃗ b 1,2 is the bisector of the beam direction vectors, � ⃗ v b1 and � ⃗ v b2 � ⃗ b 1,2 , results in a vector perpendicular to the fringe direction. The specific direction of the vector (pointing towards the top of the page or bottom of the page in Fig. 7) depends on the sign of the frequency shift that generates the moving fringe pattern. This transformation can be completed for each of the three measurement volumes, resulting in the transformation between the non-orthogonal, fringe-perpendicular velocity components, C i , and orthogonal surface-normal components, U i , as shown in Eq. 2. Here, � ⃗ d g,N is laser beam direction for the green Verdi V6 laser for N frequency shift (0 for unshifted, 40 for + 40 MHz, and 80 for − 80 MHz), � ⃗ d b,N is laser beam direction for the blue Genesis MX laser for N frequency shift (0 for unshifted and 60 for + 60 MHz), and � ⃗ b 0,N is the vector bisecting the two laser beam directions. The components in the second row are negated due to the negative frequency shift for this component. Extracting the orthogonal velocity components requires the inverse transformation matrix to be calculated.
Once the data is in local surface-normal coordinates, an additional transformation into wall-shear stress coordinates is completed via a rotation matrix applied to the data, as shown in Eq. 3. This is completed by rotating about the surface-normal x 2 -axis by W = arctan(U 3 ∕U 1 ) w at the nearestwall profile point.
The Doppler frequency signals are obtained using AURStudio from AUR, Inc, which uses Fourier domain processing for Doppler frequency estimates from photodetector signals. One block of samples (ranging from 5000 to 15,000 samples) was taken at each measurement point above the wall. Near-wall points were collected over approximately one minute of sampling time to prevent seed particles from collecting on the acrylic windows, and thus typically had fewer samples than points taken above x 2 ≈ 0.5 mm. Data further above the surface was collected over approximately one to three minutes of sampling time. During post-processing, outlying points from histograms were removed. Mean velocity results required a minimum of 500 samples to be accepted, while Reynolds stresses required a minimum of 1000 samples. Most points had a minimum of 4000 valid samples after processing. The impact of the varying number of samples for accepted datasets is included in the uncertainty via statistical uncertainty analysis. Nominal x 2 -heights at each point in the profile were corrected by fitting the nearest-wall data to U + 1 = x + 2 for data in the linear sublayer. The measurement volume has approximate dimensions of 63 μm × 63 μm × 50 μm (approximately 2-4 viscous units).
Measurements were achievable at an approximate minimum and maximum height above the wall of 0.1 mm and 30 mm (four viscous units and 1400 to 2000 viscous units), respectively, though some stations were able to capture data below this height. The flow was seeded using a mineral oil smoke machine (MDG MAX 3000 APS) with particles ranging in size from 0.5 to 1 μ m, resulting in Stokes numbers small enough to capture the full range of relevant scales.
Data was collected at 13 locations on the bump surface, as well as at the inflow to the test section and on the flat wall just upstream of the hill. For this study, four locations on and upstream of the hill are of interest at Re H = 250,000 , as shown in Table 2. These specific stations were chosen for this analysis because they are all in attached flow regions and have at least one point in the profile that is in or very near the edge of the linear sublayer, allowing for a direct calculation of w from the velocity gradient. The ratio of the centerline inflow boundary layer thickness to radius of curvature magnitude at the stations of interest on the hill is shown in Fig. 8. The sign of the curvature changes around halfway up the height of the hill, indicating a change from concave curvature ( inflow ∕R > 0 ) to convex curvature ( inflow ∕R < 0 ). The locations of all LDV data stations, including those not included in this particular study, are indicated in Fig. 9 by black filled circles.
Rigorous uncertainty quantification encompassing many uncertainty sources (e.g. article geometry uncertainty effects and repeatability of measurements) is critical for the BeV-ERLI Hill experiments to ensure the results are suitable for CFD validation. Comprehensive uncertainty quantification has been conducted on the LDV measurements, encompassing uncertainty from instrumentation sources including statistical convergence, signal processor resolution, precision in beam alignment, laser coherence, and the precision of the rotation matrix to transform the data from the nonorthogonal probe coordinates to an orthogonal coordinate system. In addition, the Error Sampling Method of Smith and Oberkampf (2014) has been applied to the data, comparing results across repeated runs to isolate sources due to uncertainty in the linear traverse motion and variations in flow conditions between runs.
The uncertainty components of these two sources (instrumentation error, dominated by bias error sources, and error sampling, dominated by random error sources) have been combined via a root-sum square to arrive at the total uncertainty interval for the mean velocities and turbulence components in global tunnel coordinates. For data given in wall-shear stress coordinates, the additional uncertainty due to the additional coordinate system transformations have also been incorporated into the results, resulting in comprehensive uncertainty estimates for the mean velocities and Reynolds stresses. Uncertainties for derived components, such as shear stress and turbulence production components, was completed via propagation of uncertainties.
Each profile has its own set of uncertainties at each point given the specific setup and conditions for that profile, and, as such, error bars are shown for each dataset and are of different magnitudes at different stations. Representative 95% confidence uncertainties are given in Table 3.

Results
Laser Doppler velocimetry data at the four stations of interest are supplemented by surface static pressure contours, oil flow visualization, and RANS CFD simulation results to analyze the global flowfield and the boundary layer history leading up to stations W2, W3, and W5. In addition, experimental data from the three-dimensional turbulent boundary layer flow on the flat wall around a wing-body junction of Ölçmen et al. (2001) is compared with the experimental data on the BeVERLI hill under similar skewing conditions. Pressure gradient and flow history effects are of primary focus, while effects of the hill curvature will be studied in future work. Figure 9a shows the mean freestream pressure coefficient contours over the entire hill at Re H = 250,000 . As shown, the pressure gradient changes sign multiple times over the surface in both the streamwise and spanwise directions as the flow away from the centerline is accelerated around the sides of the model. These strong streamwise and spanwise pressure gradients and 3D curvature induce strong near-wall crossflow. LDV profile locations are indicated by the black circles. The four stations of interest for this study are labeled W1, W2, W3, and W5. Figure 9b shows a focus view of the stations of interest on the hill, including wall-shear stress direction (from the mean flow angle of the nearest LDV measurement to the wall) indicated by black vectors, and

(b) A cross-section view of the hill cur-
vature along the grey dashed line in (a). δ inflow /R is colored according to the contour in (a). FA is the mean flow angle closest to the wall relative to the global x 1 axis, inflow ∕R is the ratio of the inflow tunnel boundary layer thickness to the local radius of curvature at each station, and the mean flow direction x + 2 ≈ 450 indicated by blue vectors. The maximum x + 2 captured at each profile is given in Table 2. The flow begins largely two-dimensional in an adverse pressure gradient (APG) at station W1 of X 1 = 1.3 . As the flow approaches the hill, the APG increases in strength up to the stagnation region on the front of the model, after which it is accelerated over the windward face of the hill. Along the sides, stations W2 and W3 experience similar pressure gradients and pressure gradient histories. At these stations, the flow direction next to the wall is skewed nearly 40 • relative to the upstream flow direction, indicating strong deformation under the lateral pressure gradients present upstream of these stations. At W2, the direction of � ⃗ = ( X 1 , 0, X 3 ) is only just slightly angled off the horizontal axis, indicating that the streamwise pressure gradient effects should dominate at this location. The spanwise pressure gradient impact becomes much more significant at W3, as � ⃗ tilts significantly in the positive X 3 -direction as the spanwise pressure gradient begins to dominate.
In contrast, station W5 experiences an entirely different pressure gradient and pressure gradient history. The local streamwise pressure gradient remains favorable, as the nearwall flow direction remains nearly perpendicular to the local pressure coefficient contours. However, the lateral pressure gradients have induced milder skewing. Here, the flow direction next to the wall is yawed −4.4 • relative to the upstream flow direction and � ⃗ indicates that the streamwise and spanwise pressure gradients are of a similar magnitude, although the spanwise component continues to dominate. The magnitudes of X 1 and X 3 are also the largest at this station at − 4.4 and 5.4, respectively.
The surface static pressure coefficient contour at Re H = 250,000 is shown overlaid with oil f low visualization images at Re H = 250,000 and Re H = 650,000 in Fig. 10 and with streamlines predicted by a k − SST simulation of the flow in ANSYS Fluent at Re H = 250,000 . In the oil flow visualization, gravity effects due to the hill's position on the side wall of the facility can be seen, especially in the lower Reynolds number case, when compared with the streamlines predicted by RANS CFD on the windward face of the hill. At Re H = 250,000 , the separated flow in the wake was not strong enough to prevent the oil mixture from flowing towards the floor under gravity (with the gravity vector pointing in the −X 3 direction). This also likely affected the windward side of the bump as well, as the windward streamlines are noticeably tilted more strongly in the direction of the floor than the streamlines at Re H = 650,000 and in the CFD. The effect of gravity on the oil flow visualization is supported by the similarity of the mean flow direction in the sublayer to the streamlines predicted by CFD. ±0.00020

Mean velocities and turbulence
LDV measurements were collected in surface normal coordinates, as shown in Fig. 2, and transformed into local wall-shear stress coordinates for analysis. In this coordinate system, the x 1 -axis points in the estimated direction of the wall-shear stress vector at the wall via the mean flow angle closest to the surface, indicated in Fig. 9b by black arrows. The x 2 -axis is perpendicular to the local surface tangent plane, and the x 3 -axis completes a right-handed coordinate system. As shown in Fig. 11, the flow begins largely 2D with near-zero components in U 2 and U 3 at W1. However, the flow experiences vastly different physics in different regions over the BeVERLI Hill, with varying levels of skewing and streamwise acceleration depending on the location on the model, as discussed. Though the U 2 and U 3 components are near-zero near the wall due to the wall-shear stress coordinate system, the U 3 components for all stations on the bump increase significantly further away from the surface due to the flow skewing. The U 2 components are near-zero for all profiles, although all profiles on the hill experience a small, increasingly negative component of U 2 as a function of height above the surface, especially at Station W5. This is an artifact of the strong tilt of the local surface normal ( 53 • relative to the X 1 -axis) wall-normal coordinate rotation. Near the wall, the flow is nearly tangential to the surface. Further above the wall, the body of the hill acts as a constriction, accelerating the flow and compressing the streamlines together to squeeze past the body. In the wallnormal coordinate system, this appears as a −U 2 component, as the streamlines passing through this profile are accelerated strongly in the outer region and squeezed together. This results in streamlines that are no longer parallel to the local surface curvature and appear to have a negative wall-normal velocity component.
The behavior in the Reynolds stresses in Figs. 12 and 13 is even more varied between stations than the mean velocities. Select Reynolds stress components are shown in Fig. 13 normalized on the friction velocity for comparison with equilibrium 2D TBL behavior.
stresses are approximately zero within uncertainties at W1. The enhanced outer peak in u �2 1 ∕U 2 ∞ is also consistent with an adverse pressure gradient turbulent boundary layer, like that of Harun et al. (2013). The shear stress u � 1 u � 2 ∕u 2 appears to be reduced compared to the 2D ZPG case, unlike the typical increase in Reynolds shear stress that would be expected in an APG TBL, as seen when compared to the zero-pressure gradient (ZPG) direct numerical simulation (DNS) data of Schlatter and Örlü (2010) in Fig. 13.
In contrast, all three stations on the bump have significant components of both normal and shear stresses. Stations W2 and W3, located near each other and in a region of similar pressure gradient and pressure gradient histories, experience very similar turbulence behavior, as well as significant skewing of the near-wall flow relative to the tunnel X 1 -axis and local concave curvature (see Fig. 8). While concave curvature generally has a destabilizing effect on a 2D TBL, Baskaran et al. (1990) found that the combined effects of mild convex and concave curvature and mean flow threedimensionality were small compared to the impact of curvature and three-dimensionality when applied individually, so it is likely that the difference in local pressure gradient and the different pressure gradient history of the flow in this region is a larger influence. The behavior between the stations primarily differs in u �2 1 ∕U 2 ∞ . The inner-region peak in u �2 1 ∕U 2 ∞ is not captured at station W2, and the behavior The gravity vector in is in the −X 3 direction, indicated by the � ⃗ g vector symbol in the top left corner above the x + 2 peak of u ′2 1 is lower in magnitude than at nearby station W3.
Station W5 experiences very different turbulence behavior than the other stations due to the different flow history, pressure gradients, and local curvature in this region. Here, the flow most notably experiences a sharp increase in the peak of u �2 1 ∕U 2 ∞ , an increase in the magnitude of u �2 3 ∕U 2 ∞ , and an inflection point in u � 1 u � 2 ∕U 2 ∞ in the log layer. As noted in Fig. 9a, there is an extended reach upstream of W5 which experiences a strong favorable pressure gradient. These turbulence reductions seen at W5 could be due to the strong flow acceleration sustained over a substantial distance as the flow is accelerated over the windward face of the hill. Insight into the significant difference between the turbulence at Station W5 and the other stations can be provided by analyzing the Reynolds stress production. The production tensor, P ij , is shown in Eq. 4 (Pope 2000).
Due to the nature of the LDV measurements, only surface-normal gradients were captured. The production is thus approximated for the Reynolds stresses by assuming that only surface-normal mean velocity gradients are significant and neglecting mean velocity gradients in the x 1 -and x 3 -directions, as shown in Eq. 5.
For all stations, the production of u ′ 1 u ′ 2 is dominated by the gradient of U 1 , as this component is more than an order of magnitude larger than the velocity gradients in the other coordinate directions. The difference in shear stress behavior between stations W2 and W3 and station W5 can be attributed to higher magnitude u ′2 2 behavior and the inflection midway through the profile, which is not present in the other stations. The production throughout the profiles are shown in Fig. 14. (5a) Fig. 11 Components of the mean velocity in local wall-shear stress coordinates, normalized by the freestream velocity. Symbols as in Table 2 Page 11 of 16 38 The behavior between profiles begins to differ in the other stress components due to the significant U 2 ∕ x 2 parameter at W5, where U 2 ∕ x 2 is of the same order of magnitude as U 3 ∕ x 2 . At W2 and W3, due to the smaller magnitude of U 2 ∕ x 2 , the −u �2 2 U 3 x 2 term dominates in the production of u ′ 2 u ′ 3 . In contrast, at W5, the term containing U 2 ∕ x 2 is non-negligible, and thus contributes to a larger magnitude u ′ 2 u ′ 3 profile at W5 as the production is more significant at this station, as shown in Fig. 13b. The significant U 3 ∕ x 2 and u ′ 2 u ′ 3 components together then contribute to enhance the production of u �2 3 ∕U 2 ∞ at W5, resulting in a larger component of u �2 3 ∕U 2 ∞ in comparison to the other stations.

Total shear stress
In the wall-shear stress coordinate system, the total shear stress is composed of the viscous component, U 1 ∕ x 2 , Fig. 12 Components of the Reynolds stress in local shearstress coordinates, normalized by the freestream velocity. Symbols as in Table 2 which is maximum in the viscous sublayer, and the turbulent component, The wall-parallel momentum equations can be simplified to focus upon the behavior of the flow in the limit of being in the viscous sublayer by ignoring the inertial terms and Reynolds shear stresses and assuming that the x 1 -axis is in the wall-shear direction, resulting in Eqs. 6 and 7. Fig. 13 Selected Reynolds stresses in wall-shear stress coordinates, normalized by u and compared to the 2D equilibrium DNS data to the DNS data of Schlatter and Örlü (2010). Symbols as in Table 2 Fig. 14 Reynolds shear stress production, estimated as in Eq. 5 In this wall-shear stress coordinate system, an FPG in the x 1 -direction indicates that the flow is being accelerated in the wall-shear direction. Skewing is caused by a lateral pressure gradient perpendicular to the local wall-shear direction, inducing a nonzero U + 3 component. If the pressure gradient is assumed to be approximately constant as a function of height in the viscous sublayer, then an FPG in the wall-shear direction will cause the viscous stress to reduce more rapidly from the wall in wall units than a ZPG or APG due to the influence of the pressure term.
The turbulent shear stress is also very sensitive to flow conditions. The presence of mean streamwise vorticity is known to impact the quasi-coherent turbulent structures that are responsible for the production of u ′ 1 u ′ 2 and u ′ 2 u ′ 3 Johnston and Flack (1996) and lead to the reduction of these shear stress components. However, Flack and Johnston (1998) have also found that turbulent increases outside of the buffer layer in regions of increasing three-dimensionality. This can also be seen at some stations in the data of Ölçmen et al. (2001). For both experiments, although u � 1 u � 2 ∕u 2 decreased in magnitude as the three-dimensionality of the flow increased, an increase in the magnitude of u � 2 u � 3 ∕u 2 compensated for this behavior and led to an increase in turbulent outside of the buffer layer.
The viscous and turbulent shear stress components of the total shear stress in the BeVERLI Hill stations of interest are normalized by w , computed via w = u 2 , in Fig. 15 and compared with the 2D ZPG equilibrium boundary layer DNS of Schlatter and Örlü (2010). The uncertainty of u is estimated as ±0.025 m/s, leading to an uncertainty in w of approximately 6%. As shown in Fig. 15, stations W2 and W3 both experience a region below x + 2 ≈ 200 where the turbulent shear stress is significantly reduced and not able to balance out the viscous stress it drops outside of the linear sublayer. This results in a total stress distribution that is far from the equilibrium 0.95 ≤ ∕ w ≤ 1 that is typically present in 2D ZPG boundary layers and is assumed in some LES equilibrium-stress models, such as those of Piomelli et al. (1989) and Marusic et al. (2001).
Above x + 2 ≈ 200 , there is a sharp increase in the turbulent stress, as was also observed by Flack and Johnston (1998) (though in their case, this was observed above x + 2 ≈ 50 ). There is, therefore, a region below x + 2 ≈ 200 that is dominated by viscous and wall-effects, and a region above this point of the boundary layer where the effects of the crossflow stress u ′ 2 u ′ 3 begin to contribute more significantly as the flow turns away from the local wall-shear stress direction. At Stations W2 and W3, this significant drop in turbulent shear stress is expected, as a reduction in turbulent shear stress is often seen in 3D flows, such as the experimental flow of Bradshaw and Pontikos (1985) and the DNS simulations of Moin et al. (1990), where the significant near-wall skewing and transverse strain suppresses the pressure-strain term in the turbulent transport.
In contrast, at station W5, the flow behavior is quite different due to the different pressure gradient and flow history in this region. The shear stress behavior is depressed throughout the full profile below the equilibrium values and experiences an initial rise to a local maximum, decrease to a local minimum, before rising again towards the end of the LDV profile. The decreased u �2 2 ∕u 2 profile relative to stations W2 and W3 is consistent with what would be expected from analysis of the production of each turbulence component, as the production of P 12 and P 22 are linked through the u ′2 2 component. The turbulent shear stress component in this region is dominated by the contribution from u ′ 1 u ′ 2 , which suggests that the same mechanisms that led to a lowered u � 1 u � 2 ∕u 2 are contributing to the strongly nonequilibrium turbulent shear stress behavior.
The profiles are compared in Fig. 17 with the 3D TBL on a flat wall around a wing-body junction from Ölçmen et al. (2001) in regions that experience similar pressure gradient and skewing behavior but without wall curvature. A schematic of the local freestream and wall-shear stress directions in their experiment is shown in Fig. 16. As shown, the flow experiences increasing skewing as the near-wall flow feels the impact of the airfoil presence more strongly than the freestream flow. After Station O3, the magnitude of the skewing relative to the freestream flow begins to decrease, as the wall-shear stress and freestream flow directions become nearly identical. Of the nine flow stations, four are  Table 2. Open symbols indicate viscous stress components, U + 1 ∕ x + 2 , and filled symbols indicate turbulent shear stress components, √ of particular interest for comparison to the skewed flow on the BeVERLI Hill. At Stations O3 and O4 of Ölçmen et al, the flow is extremely skewed ( > 30 • from the upstream flow direction), and streamwise pressure gradient is favorable as it passes around the side of the wing-body junction. This global flowfield behavior is similar to that of stations W2 and W3 of the BeVERLI flow, and, as shown, results in similar turbulent shear stress behavior, as shown in Fig. 17a. However, the data of Ölçmen et al is missing in the region in which Stations W2 and W3 begin experiencing a sharp increase in the turbulent shear stress, so it is not possible to identify if this behavior is also present in the juncture flow.
A similar trend follows for Ölçmen et al's stations O6 and O7 and the hill's station W5, as seen in Fig. 17b. At both stations, the flow is much less skewed than at previous stations as it recovers past the junction and is also experiencing a stronger streamwise FPG than further upstream in the flow, similar to W5. However, while the shapes of the profiles are similar, the magnitudes are very different. Notably, stations O6 and O7 experience a region through the buffer layer where the turbulent shear stress nearly matches the 2D DNS, before a sudden drop in turbulent shear stress above the buffer layer, and a subsequent increase above x + 2 ≈ 200 . The behavior at W5 follows a similar trend, but severely depressed in magnitude. The pressure gradient magnitudes at W5 are significantly stronger than those at O6 and O7 and are sustained over a substantial distance. This suggests that the strong streamwise and spanwise pressure gradients could be a source of the substantial turbulence reductions seen at W5, which then in turn reduce the turbulent shear stress distribution. The station at W5 is also located in a region of convex curvature on the hill, and curvature effects of the BeVERLI Hill may also be a cause behind the difference in magnitude between these cases.

Conclusion
LDV measurements of the flow upstream of and on the three-dimensional BeVERLI Hill model were used to analyze the fluid physics at those stations. As shown, the selected stations experience significantly different mean flow and turbulence behavior due to differing pressure gradients, radii of curvature, and history effects at each station. Of these effects, local pressure gradient and skewing appear to be the strongest influences on the turbulence at each station. Curvature effects due to the 3D hill geometry were not specifically examined in this study and will instead be a focus of future work.
Analysis of the BeVERLI Hill data and comparison with the 3D TBL on a flat wall around a wing-body junction from Ölçmen et al. (2001) show there is no region in the flow over the BeVERLI Hill where equilibrium assumptions regarding the shear stress distribution are valid. This is especially severe in strongly skewed flow, such as at stations W2 and W3. In accelerated regions of the flow with milder flow skewing, such as at station W5, the turbulent stresses were significantly depressed throughout the full profile relative to the 2D equilibrium  . 17 Comparison of turbulent and viscous shear stress behavior of selected stations on the BeVERLI Hill with the 2D equilibrium DNS data of Schlatter and Örlü (2010) and the 3D flat-wall TBL of Ölçmen et al. (2001). a Stations W1, W2, and W3 of the BeVERLI Hill compared with stations O3 and O4 of Ölçmen et al. (2001). b Station W5 of the BeVERLI Hill compared with stations O6 and O7 of Ölçmen et al. (2001) (a) BeVERLI symbols as in Table 2, : Station O4. Table 2, : Station O7.

(b) BeVERLI symbols as in
case. Similar behavior seen in the BeVERLI Hill flow was also observed in the flat wall flow around a wing-body juncture of Ölçmen et al. (2001) in regions with similar skewing and pressure gradient behavior. The comparisons match very well with stations W2 and W3, located in regions with milder curvature and strong flow skewing. The behavior at W5 follows the same trend as the compared stations but with a significantly lowered magnitude, which may be due a combination of differences in the magnitude of the pressure gradients between the two compared flows and curvature effects. Turbulence models that rely on this equilibrium shear stress assumption will not be able to accurately predict the turbulence behavior of 3D flow.