Continuous Eddy Simulation vs. Resolution-Imposing Simulation Methods for Turbulent Flows

Abstract

The usual concept of simulation methods for turbulent flows is to impose a certain (partial) flow resolution. This concept becomes problematic away from limit regimes of no or an almost complete flow resolution: discrepancies between the imposed and actual flow resolution may imply an unreliable model behavior and high computational cost to compensate for simulation deficiencies. An exact mathematical approach based on variational analysis provides a solution to these problems. Minimal error continuous eddy simulation (CES) designed in this way enables simulations in which the model actively responds to variations in flow resolution by increasing or decreasing its contribution to the simulation as required. This paper presents the first application of CES methods to a moderately complex, relatively high Reynolds number turbulent flow simulation: the NASA wall-mounted hump flow. It is shown that CES performs equally well or better than almost resolving simulation methods at a little fraction of computational cost. Significant computational cost and performance advantages are reported in comparison to popular partially resolving simulation methods including detached eddy simulation and wall-modeled large eddy simulation. Characteristic features of the asymptotic flow structure are identified on the basis of CES simulations.

1. Introduction

The characteristic features of computational methods for simulations of turbulent flows of practical relevance are well known. Reynolds-averaged Navier–Stokes (RANS) methods are computationally very efficient, but they suffer from their inability to correctly describe the physics of separated flows, which appear in most applications [1]. Large eddy simulation (LES) methods have a significantly higher predictive power, but they are often unaffordable because of the huge computational cost required for simulations of high Reynolds number ($Re$) turbulent flows involving complex geometries [1,2,3,4,5]. The most promising approach to overcome these issues is the design of hybrid RANS-LES methods that combine elements of RANS and LES equations [2,5,6,7,8,9]. In most applications, detached eddy simulation (DES) [10,11,12,13] and wall-modeled LES (WMLES) [14,15,16] methods are applied. However, such hybrid RANS-LES suffer from basic problems, too. There is significant uncertainty of predictions due to different models applied, regions where different models and grids are applied, different mesh distributions, and set-up options to manage the information exchange between flow regions [5]. Hence, all such predictions need validation data, which are often unavailable. In addition, although such hybrid RANS-LES are designed to reduce the cost of LES, it turns out that their computational cost can be very significant, even comparable to LES cost if performance deficiencies need to be addressed [17].

A closer look at concepts applied reveals that the predominant strategy to design computational methods for simulations of turbulent flows is to impose a desired flow resolution, usually by the setting of model viscosity. For example, a relatively large model viscosity aiming at zero flow resolution is applied in RANS methods, a relatively small model viscosity aiming at an almost complete flow resolution is applied in LES methods. In (partially) resolving simulation methods (including DES and WMLES), filter width $\Delta$ is used as characteristic model length scale in simulations on relatively coarse grids. The idea is to impose on the model side a resolution (model viscosity) ranging between corresponding RANS and LES resolutions (model viscosities). The problem is that the resolution imposed in this way does not need to determine the actual flow resolution (which can be measured as explained below). The actual flow resolution involves two ingredients, model contributions and resolved contributions, and the latter are not controlled by the model contribution away from RANS and LES regimes.

The latter has significant consequences.

• In this case, the hybrid RANS-LES model cannot properly handle the transition between RANS and LES regimes, which requires the model to increase (decrease) its contribution in response to a relatively low (high) actual flow resolution. Given the uncertainty of the actual flow resolution, there is implied uncertainty of simulation results. Hence, often unavailable validation data are required, and reliable predictions of very high-$Re$ turbulent flows are out of reach.
• The mismatch between the imposed and actual flow resolution introduces hybridization errors (see the discussion of minimal error methods below). To accomplish a desired simulation performance, finer grids (higher computational cost) are then needed to minimize hybridization errors [17].
• The imbalance between the damping of fluctuations (controlled by the modeled viscosity) and the resolved fluctuations can generate problems seen very often in hybrid RANS-LES, which require the stimulation of fluctuations to trigger the development of instantaneous turbulence or the damping of fluctuations to prevent the shifting from RANS to LES within the boundary layer.
• In particular, the modeled-stress depletion problem, the slow development of resolved motion after a RANS-to-LES transition, is caused by this imbalance [10,18]. This requires then empirical techniques to transfer energy from the modeled to the resolved scales and vice versa.

Hybrid RANS-LES methods that do not suffer from the mismatch of imposed and actual resolution were recently introduced as continuous eddy simulation (CES) methods [5,19,20,21,22,23,24,25]. Their characteristic feature is that CES methods enable the model to read and properly respond to the actual flow resolution. The key to accomplish the latter is a strict mathematical derivation via variational analysis. So far, the suitability of CES methods was only tested for periodic hill flow at a relatively high $Re$. This application includes relevant features such as separation, recirculation, and natural reattachment, but its complexity is not comparable with complex flow configurations seen in usual applications. The motivation of this paper is a comprehensive evaluation of benefits of CES methods for a moderately complex flow, the NASA wall-mounted hump flow that is frequently used for the evaluation of turbulence models. Specific questions addressed are the following ones.

• What are characteristic general differences between CES methods and simulation methods that impose flow resolution independent of the actual resolved motion?
• Do the conceptual advantages of CES methods imply specific practical advantages, e.g., in regard to reliability of simulation results and computational cost?
• Which insight can be obtained from the application of CES methods to extreme $Re$ regimes that cannot be reliably studied by any other approach?

The paper is organized in the following way. The modeling and computational approach is described in Section 2, and the hump flow simulation set up is described in Section 3. Section 4, Section 5 and Section 6 deal with the evaluation of CES model performance, the comparison with RANS and other hybrid RANS-LES, and the comparison with LES methods, respectively. Very high-$Re$ regimes are studied in Section 7, and conclusions are presented in Section 8.

2. Modeling and Computational Approach

The CES models considered and their computational implementation are introduced next.

2.1. CES Models Considered

CES hybridization can be applied to at least most popular turbulence models [19,24]. Typically, a RANS model involving a transport equation for the modeled kinetic energy k combined with a scale equation for the modeled dissipation rate $ϵ=k/\tau$ or turbulence frequency $\omega =1/\tau$ is considered, where $\tau$ refers to the modeled turbulence dissipation time scale. Variational analysis of these equations is applied, which covers the transition from RANS to LES by including partially resolving regimes. The latter enables it to establish relationships between model parameters and resolution measures (as, e.g., $L_{+}$, see below). Used in a fully resolving mode, it is relevant to note that these equations correspond to LES without any need to introduce an artificial length scale as given by filter width $\Delta$. An essential feature of this approach is the possibility to address the hybridization in several ways based on the k equation or the scale equation.

Here, the approach is applied in regard to a two-equation eddy viscosity model (EVM), in particular a $k-\omega$ model. The model is based on continuity and momentum equations,

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{U}}_i}{\partial x_i}=0,\frac{\stackrel{\sim }{D}{\stackrel{\sim }{U}}_i}{\stackrel{\sim }{D}t}=-\frac{1}{\rho }\frac{\partial \stackrel{\sim }{p}}{\partial x_i}+\frac{\partial \left(2\nu {\stackrel{\sim }{S}}_{ij}\right)}{\partial x_j}-\frac{\partial {\tau }_{ij}}{\partial x_j}.\end{array}$$

(1)

The sum convention is applied throughout this paper. The tilde notation refers to space-averaged variables. Specifically, expression $\stackrel{\sim }{D}/\stackrel{\sim }{D}t=\partial /\partial t+{\stackrel{\sim }{U}}_j\partial /\partial x_j$ represents the filtered Lagrangian time derivative. The velocity vector components are denoted as ${\stackrel{\sim }{U}}_i$, while $\stackrel{\sim }{p}$, $\rho$, and $\nu$ represent pressure, fluid density, and kinematic viscosity, respectively. Rate-of-strain tensor ${\stackrel{\sim }{S}}_{ij}$ is given by ${\stackrel{\sim }{S}}_{ij}=(\partial {\stackrel{\sim }{U}}_i/\partial x_j+\partial {\stackrel{\sim }{U}}_j/\partial x_i)/2$. In the momentum equation, the right-hand side (RHS) incorporates the unknown model stress tensor, ${\tau }_{ij}$. According to the EVM concept, ${\tau }_{ij}=2k{\delta }_{ij}/3-2{\nu }_t{\stackrel{\sim }{S}}_{ij}$ is applied here. Here, ${\delta }_{ij}$ denotes the Kronecker symbol, and ${\nu }_t$ represents the modeled viscosity of the flow. The latter is computed by ${\nu }_t=C_{\mu }k/\omega$, where $C_{\mu }$ is a model parameter with standard value $C_{\mu }=0.09$.

Three hybrid CES RANS-LES models are considered here; see

Table 1. Summary of hybrid CES models considered where $R=L_{+}^2$ [20].

Model	CES Hybrid Equations	Mode Control
CES-KOS:	$\begin{array}{l}\frac{Dk}{Dt}=P-ϵ+D_k,\\ \frac{D\omega }{Dt}=C_{{\omega }_1}{\omega }^2\left(\frac{P}{ϵ}-{\beta }^{*}\right)+D_{\omega }+D_{\omega c}\end{array}$	${\beta }^{*}=1+R(\beta -1)$
CES-KOK:	$\begin{array}{l}\frac{Dk}{Dt}=P-{\psi }_{\beta }ϵ+D_k,\\ \frac{D\omega }{Dt}=C_{{\omega }_1}{\omega }^2\left(\frac{P}{ϵ}-\beta \right)+D_{\omega }+D_{\omega c}\end{array}$	${\psi }_{\beta }=1+\beta -{\beta }^{*}$
CES-KOKU:	$\begin{array}{l}\frac{Dk}{Dt}={\psi }_{\beta u}^{-1}(P-{\psi }_{\beta u}^2ϵ+D_k),\\ \frac{D\omega }{Dt}=C_{{\omega }_1}{\omega }^2\left(\frac{P}{ϵ}-\beta \right)+D_{\omega }+D_{\omega c}\end{array}$	${\psi }_{\beta u}^2={\psi }_{\beta }$

. CES-KO refers to CES performed with the $k-\omega$ model. In particular, -KOS (-KOK, -KOKU) refers to CES-KO performed by modifying the scale equation (the k equation, the turbulence time scale in the k equation according to unified RANS-LES) [20]. Differences of these models are discussed in Section 4. The equations involve the production of modeled kinetic energy $P={\nu }_t{S}^2$, which involves the characteristic strain rate $S={\left(2{\stackrel{\sim }{S}}_{mn}{\stackrel{\sim }{S}}_{nm}\right)}^{1/2}$. The turbulent transport terms are given by

$$\begin{array}{c}{D}_k=\frac{\partial }{\partial x_j}\left[\left(\nu +{\nu }_t\right)\frac{\partial k}{\partial x_j}\right],D_{\omega }=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right],\\ D_{\omega c}=\frac{C_{\omega }}{k}(\nu +{\nu }_t)\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}.\end{array}$$

(2)

The model constants involved have values $C_{\omega 1}=0.49$, $\beta =1.63$, and ${\sigma }_{\omega }=1.8$ [26]. The integration of the $\omega$ equation through the viscous sublayer can cause numerical errors. To avoid this problem, $\omega =2\nu /d^2$ was used at the first cells above the wall, where d refers to the distance from the wall to the cell center of the first cell. These equations are applied in conjunction with the elliptic blending approach used in our previous work to account for wall damping [26] (in an RANS-mode for the entire flow field).

The essential new feature of the equations in Table 1 is the mode control relation, which relates the model parameters, ${\beta }^{*},{\psi }_{\beta },{\psi }_{\beta u}^2$, to flow resolution measures. These relations were obtained by variational analysis; see the beginning of this subsection. To simplify the analysis, cross-diffusion term $D_{\omega c}$ (which acts as a damping function) is neglected. In addition, the substantial derivatives of k and $\omega$ are neglected, too. It is possible to involve the latter terms [22], but this approach implies the need to consider more complex equations. Given the objective to demonstrate the basic advantages of the CES approach, no attempt is made to consider this option. In particular, ${\beta }^{*},{\psi }_{\beta },{\psi }_{\beta u}^2$ are determined as functions of $L_{+}=L/L_{tot}$, which is the turbulence length scale resolution ratio involving modeled (L) and total contributions ($L_{tot}$) [19]. The modeled contribution is calculated by $L={〈k〉}^{3/2}/〈ϵ〉$. The brackets refer to averaging in time. The total length scale is calculated correspondingly by $L_{tot}=k_{tot}^{3/2}/{ϵ}_{tot}$. In these relations, we have $k_{tot}=〈k〉+k_{res}$, where $k_{res}=\left(〈{\stackrel{\sim }{U}}_i{\stackrel{\sim }{U}}_i〉-〈{\stackrel{\sim }{U}}_i〉〈{\stackrel{\sim }{U}}_i〉\right)/2$. Correspondingly, ${ϵ}_{tot}$ is the sum of modeled and resolved contributions, ${ϵ}_{tot}=〈ϵ〉+{ϵ}_{res}$. Here, the resolved contribution is calculated by ${ϵ}_{res}=\nu \left(〈\partial {\stackrel{\sim }{U}}_i/\partial x_j\partial {\stackrel{\sim }{U}}_i/\partial x_j〉-〈\partial {\stackrel{\sim }{U}}_i/\partial x_j〉〈\partial {\stackrel{\sim }{U}}_i/\partial x_j〉\right)$.

The dependence of model parameters ${\beta }^{*},{\psi }_{\beta },{\psi }_{\beta u}^2$ on $L_{+}$ has essential implications. The current computational strategy for turbulent flow simulations is the use of models being unaware of the actual flow resolution, which has significant consequences as discussed in the introduction. In contrast to that, the CES approach informs the model via resolution measures as given by $L_{+}$ about the actual flow resolution. The latter enables the model to adjust its contribution to the simulation: the model contribution decreases (increases) if there is (there is not) a significant amount of resolved motion. This ability is the essential requirement of a functional RANS-LES swing. It is also worth noting that the CES approach minimizes the hybridization error [24]: there is no mismatch between the resolution imposed by the model and the actual resolution seen in simulations, which minimizes computational cost.

2.2. Computational Approach

The three hybrid CES models considered (the CES-KOS, CES-KOK, and CES-KOKU models) were implemented in the OpenFOAM open-source platform, a widely used computational fluid dynamics tool [27]. Numerical simulations were conducted using a finite-volume-based method, and a time-marching scheme was utilized to solve the unsteady turbulent flow equations on various numerical grids. The iterative solution process involved the application of both implicit and explicit time-step schemes. It should be noted that simulations employing an explicit Euler scheme suffer from instability issues. Therefore, to ensure stability and accuracy, the implicit backward-Euler time integration scheme, known for its second-order accuracy in time, was employed as one of the standard Newton sub-iteration methods in this study.

The spatial discretization of the diffusive and gradient terms in the governing equations was performed using a central difference scheme with Gaussian integration. For the convective part of the momentum equation, a second-order GammaV scheme was initially applied for a few time steps. Then, we switched to the bounded second-order van Leer scheme after achieving improved stability. The turbulent transport terms were also discretized using the second-order van Leer scheme, ensuring continuity and numerical stability. The consideration of the van Leer scheme in the CES framework is motivated by its ability to effectively capture discontinuities and shock waves in turbulent flow, particularly in wall-bounded separated flows. This is attributed to its ability to transport and handle steep gradients without generating spurious oscillations and unphysical negative mixing ratios [28,29].

• The remaining flow terms were discretized using the second-order limitedLinearV scheme. The pressure–velocity coupling was handled by the PISO algorithm, which dynamically adjusts the pressure gradient along the flow direction to maintain a constant mass flow rate [30]. The iterative solution process employed a preconditioned bi-conjugate gradient technique with diagonally incomplete LU preconditioning for all flow variables except pressure at each time step [31]. To solve the Poisson equation for pressure, an algebraic multi-grid solver was utilized.

3. Hump Flow Simulation Set-Up

We applied the CES and other models to moderately complex separated flow simulations: flow over a wall-mounted hump geometry, commonly referred to as the NASA hump model, which serves as a representative model for the upper surface of an airfoil. The experimental setup and design of the 2D streamwise-periodic hump flow are shown in Figure 1. The experimental studies of this flow, boundary conditions and the set-up of flow simulations are described next. Throughout the article, the standard values for hump chord length and free-flow velocity are $c=420$ mm and $U_{ref}=34.6$ m/s, respectively. Using these quantities, all variables shown below were made dimensionless.

3.1. Flow Configuration

Seifert and Pack [32] developed a wall-mounted hump model to investigate unsteady flow separation, reattachment and flow control at a high Reynolds number $Re=c{\rho }_{ref}{U}_{ref}/\mu \approx 936$ K based on chord length c and freestream velocity $U_{ref}$. Here, $\mu$ is dynamic viscosity and abbreviation $ref$ indicates the reference freestream conditions, which are determined at axial point $x/c=-2.14$. The model reflects the upper surface of a 20-thick Glauert–Goldschmied airfoil that was originally designed for flow control purposes in the early twentieth century. As a benchmark for comparison, we used the experiment conducted by Greenblatt et al. [33] without flow control. This benchmark case has been extensively documented on the NASA Langley Research Center’s Turbulence Modeling Resource webpage and has been widely used for evaluating different turbulence modeling techniques, as discussed in the 2004 CFD Validation Workshop. We see in Figure 1 a strongly convex region just before the trailing edge, which induces flow separation. The flow over the hump model exhibits significant separation at the trailing edge of the hump, accompanied by a free shear layer extending from the crest of the hump. At $x/c=0.65$, the flow over the hump model experiences a breakage, followed by reattachment at $x/c=1.1$ due to a pronounced adverse pressure gradient. The hump has a length of $c=420$ mm and a peak height of $53.7$ mm, with a predominantly two-dimensional flow and accounting for end-plate side walls. The estimated freestream Mach number is approximately $0.1$. At the high $Re=936$ K considered, accurately predicting the boundary layer’s reattachment and subsequent recovery to an equilibrium state poses significant challenges.

The discretization of the domain employed in the numerical simulation spans $L_x\times L_y\times L_z=6.14c\times 0.909c\times 0.2c$ in the streamwise (x), wall-normal (y), and spanwise (z) directions, respectively, with the origin of the coordinate system situated at the leading edge of the hump, thereby establishing an entrance point at $x/c=-2.14$ and an exit point at $x/c=4.0$. By locating the inlet patch upstream of the hump, the turbulent boundary layer is fully developed before reaching the hump region. For computational purposes, the top-wall boundary is set at $y/c=0.9009$. The span-wise boundary is normally constrained in the range of 0 to 0.4c [34,35]. In particular, You et al. [35] used $L_z=0.2c$. Thus, the present study uses a corresponding span-wise size of $0.2c$.

3.2. Boundary Conditions

Incoming boundary condition: To generate the necessary inflow, a RANS simulation of a channel flow with zero pressure gradient was performed using internal plane mapping techniques [36,37]. This precursor simulation produced inlet conditions matching experiments for velocity and turbulent kinetic energy [36]. Figure 2 (left) shows the precursor domain visually integrated with the main simulation. The precursor velocity contour was generated using a computational grid with 162 nodes in the flow direction, 62 along the walls, and additional 24 across the span. The recycling method by Simens et al. [38] was leveraged to effectively manage the mean velocity data. The mean flow velocity data was extracted from the recycling plane at $x/c=-6$ downstream of the hump leading edge, as recommended by Simens et al. [38]. However, the mapping plane was deliberately positioned at $x/c=-2.14$ to ensure that the velocity and turbulence characteristics attain a steady state by the time the flow reaches the mapping plane.

• Hence, the predicted velocity and turbulence characteristics (turbulence kinetic energy) were systematically recorded in samples at the mapping plane during each iteration over $250c/U_{ref}$ time units. This compiled a comprehensive flow library. Samples were then extracted, extrapolated, and fed in as inflow conditions at the main simulation’s inlet plane. Figure 2 (right) shows the time-averaged streamwise velocity profile at the inlet plane (mapping plane). This profile closely resembles the velocity profile reported by You et al. [35] at the same axial position.

Additional boundary conditions: The test case involves a 2D flow configuration confined within a channel bounded by two solid walls. The top wall of the channel is treated with no-stress and no-slip boundary conditions, ensuring that there is no shear stress and velocity parallel to the wall. On the other hand, the bottom wall is subjected to viscous adiabatic boundary conditions, accounting for the effect of viscosity without heat transfer. To align with experimental conditions, several simulations reported in previous workshops employed upper-wall symmetry or slip-wall boundaries, which better mimic the experimental setup. Regarding the exit boundary condition, the outflow patch is characterized by symmetric z-axis boundary conditions, ensuring a balanced flow through the boundary. Conversely, the inflow patch is treated with distinct convective boundary conditions, where the convective speed is determined based on the mean streamwise velocity at the outlet patch. This approach allows for maintaining the appropriate convective properties at the inflow boundary and ensures consistency with the downstream flow behavior.

3.3. Simulation Set-Up

A series of simulations were conducted on the hump model using four different grid resolutions: coarse, medium, fine and very fine. These grids, denoted as $G_1$, $G_2$, $G_3$, and $G_4$ respectively, are presented in

Table 2. The grids applied.

Run	Grid	Grid Points	Resolution	Inflow (Generation) Plane	Min. Wall Spacing	$\mathbf{\Delta }{\mathit{y}}^{+}$
$G_1$	Coarse	480 K	323 $\times 62\times 24$	plane mapping, $x/c=-2.14$	0.0001	1.25
$G_2$	Medium	960 K	$417\times 72\times 32$	plane mapping, $x/c=-2.14$	0.00001	0.75
$G_3$	Fine	1.7 million	$482\times 90\times 40$	plane mapping, $x/c=-2.14$	0.000005	0.25
$G_4$	Very fine	3.9 million	$723\times 135\times 40$	plane mapping, $x/c=-2.14$	0.000001	0.2

and Figure 3. Here, $\Delta y^{+}$ refers to the non-dimensional wall distance of the first grid point.

In the preliminary $G_1$ simulation, we employed a grid comprising $323\times 62\times 24$ nodes in streamwise, wall-normal, and spanwise directions, respectively. The mesh configuration was designed to concentrate nodes near the wall and around the hump model in order to accurately capture the intricacies of the boundary layer. The mesh was refined in the streamwise and spanwise directions to accommodate the overall flow dynamics. In order to prevent disturbances from propagating upstream, the grid was additionally and selectively refined near the outlet in the streamwise direction, creating a buffer zone. The streamwise $G_1$ mesh spacing corresponded to eight points per boundary layer thickness, denoted here as ${\delta }_{ref}$. In terms of wall units (derived from the wall shear and evaluated at the inlet), the $G_1$ mesh resolution perpendicular to the wall was determined as $\Delta y^{+}=1.25$. Using this initial mesh configuration, three additional levels of mesh refinement were obtained. The medium-mesh $G_2$ simulation involved $417\times 72\times 32$ nodes in the streamwise, wall-normal, and spanwise directions, while the fine-mesh $G_3$ simulation utilized $482\times 90\times 40$ nodes in the same directions. This yielded a resolution of 12 and 16 points per ${\delta }_{ref}$ for $G_2$ and $G_3$ respectively. The grid points in the wall-normal direction corresponded to $\Delta y^{+}=0.75$ for $G_2$ and $\Delta y^{+}=0.25$ for $G_3$. The hump region’s unique features required a more targeted grid concentration due to the presence of shear layers, separation, and reattachment. Hence, the $G_4$ grid has a higher mesh density concentrated on the hump surface in the near-wall area. This was achieved by refining the streamwise resolution of the $G_3$ grid by a factor of 1.5 in the inner region closer to the wall. Thus, the resulting $G_4$ grid had $723\times 135\times 40$ points, with a wall-normal spacing of $\Delta y^{+}=0.2$ and 20 points per ${\delta }_{ref}$. The spanwise direction features the same mesh distribution as the $G_3$ grid. The $G_4$ simulation results were only used to further support the main comparisons.

The iterative computations generated a cumulative maximum Courant–Friedrichs–Lewy (CFL) number of 0.4, which corresponds to $8.2\times {10}^{-3}c/a_{ref}\approx 0.1\times {10}^{-4}$ time-step for the Mach $=0.1$ case, where $a_{ref}$ is the typical speed of sound. The flow characteristics were studied at three Reynolds numbers: $Re=$ (936 K, 5 M, 10 M). Experimental data for comparisons were only available for $Re=936$ K. Before averaging span-wise turbulent quantities, all simulations were run for 500 flow-through times (FTT).

The flow resolution was determined via $0\le L_{+}=L/L_{tot}\le 1$, which characterizes the model contribution to the total characteristic length scale $L_{tot}$ of turbulent motions. In particular, $L_{+}=0$ corresponds to LES conditions where model contribution L is negligible compared to resolved contribution $L_{tot}$, whereas $L_{+}=1$ corresponds to RANS conditions where all the flow is modeled. CES simulations started with an initial value of $L_{+}=1$, equivalent to RANS modeling. Subsequently, the values of $L_{+}$ undergo dynamic updates by considering averaged values of modeled and resolved kinetic energies and dissipation rates, K, $k_{res}$, $ϵ$, and ${ϵ}_{res}$. The required variables were obtained by employing a temporal smoothing technique, beginning at the initial time and extending up to a maximum of 2000 FTT. In line with this, the resolution variables presented below represent averaged values obtained from a span of 575,000 time-steps, representing a total duration of 500 convective time units (1 CTU or FTT $=c/U_{\infty }$). After reaching the 1500 CTU mark, the mean quantities were further averaged over a period of 500 CTU to eliminate any transient effects, while also considering average values in the spanwise direction. Figure 4 illustrates the result of characteristic mean values of $L_{+}$ and $k_{+}=k/k_{tot}$. The variations correspond to expectation: a coarser grid and a higher $Re$ increase the relative amount of resolved motion. In correlation with corresponding periodic hill flow simulations [20], the grid and $Re$ variations cover almost resolving LES-type simulations via $Re=936$ K ($G_3$ case) and basically modeled RANS-type simulations via $Re=10$ M ($G_1$ case).

4. CES Model Performance

Characteristic CES features are described next. In particular, we consider the reflection of resolution variations and simulation performance. The influence of the CES version considered (see Table 1) was studied in preparation to these analyses. It was found that the use of CES-KOKU is less appropriate for the flow considered, probably because the CES-KOKU hybridization concept is not as simple and stable as CES-KOS and CES-KOK concepts. CES-KOS was found to provide the most accurate results; see the reference to corresponding CES-KOK characteristics in

Table 3. CES-KOS and CES-KOK results: Separation and reattachment points depending on the grids applied.

Methods	Separation Location ($\mathbf{x}/\mathbf{c}$)	Reattachment Location ($\mathbf{x}/\mathbf{c}$)	Bubble Length $\Delta \mathbf{x}/\mathbf{c}$	Error in Bubble LengthPrediction (%)
Exp. [33]	$0.665\pm 0.005$	$1.10\pm 0.005$	0.435	–
CES-KOS–$G_1$	0.6642	1.135	0.4708	8.2
CES-KOS–$G_2$	0.6669	1.1144	0.4475	2.9
CES-KOS–$G_3$	0.6737	1.1100	0.4363	0.3
CES-KOS–$G_4$	0.6680	1.1000	0.432	0.7
CES-KOK–$G_3$	0.6752	1.1333	0.4581	5.3

. Thus, the CES-KOS model (or simply the KOS model) is considered from now on.

4.1. Resolution Variations

Figure 5 and

Table 4. Range of velocity fluctuations for the cases described in Figure 5.

Fluc.	$\mathit{y}/\mathit{c}=0.01$	$\mathit{y}/\mathit{c}=0.05$	$\mathit{y}/\mathit{c}=0.1$
u	$-0.01\le u\le 0.01$	$-0.02\le u\le 0.02$	$0.04\le u\le 0.04$
v	$-0.01\le v\le 0.01$	$-0.02\le v\le 0.02$	$0.04\le v\le 0.04$
w	$-0.01\le w\le 0.01$	$-0.02\le w\le 0.02$	$0.04\le w\le 0.04$

show velocity fluctuations (u, v, and w) for different $Re$ and resolution conditions, in particular the $Re=5\mathrm{M}(G3,G1)$ and $Re=936\mathrm{K}(G3,G1)$ cases. The $Re=936$ K ($G_3$) case shows a variety of small-scale vortices which indicate a high level of resolution, i.e., this case is at least relatively close to being well resolved. On the other hand, the $Re=5$ M ($G_1$) case shows an instationary RANS-type regime: fluctuations are still present, but small-scale turbulence structures are not resolved anymore. The intermediate cases ($Re=936$ K ($G_1$) and $Re=5$ M ($G_3$)) represent transition regimes in between well-resolved simulations and instationary RANS simulations. Additional observations are the following ones. First, for the highly resolved case, the proximity to the wall has little impact on fluctuation patterns, the spatial extent and frequency of fluctuations is about the same. Second, grid coarsening and higher $Re$ cause continuous changes in fluctuations without sudden shifts. The most noticeable change is the spatial clustering of fluctuations, while the overall strength of fluctuations remains largely unaffected. Under instationary RANS conditions, fluctuations are still present (they do not vanish). In particular, they can be as strong as under almost resolving conditions.

Figure 6 shows local flow resolution variations $L_{+}=L/L_{tot}$, $k_{+}=k/k_{tot}$, ${ϵ}_{+}=ϵ/{ϵ}_{tot}$ caused by grid variations; corresponding variations due to $Re$ changes are presented below. $L_{+}$ profiles confirm expectations. We have $L_{+}^{G_3}\le L_{+}^{G_2}\le L_{+}^{G_1}$, indicating that coarser grids result in higher $L_{+}$ values due to a lower resolution. $L_{+}$ values become larger with smaller distance from the walls because of effective grid coarsening. Notably, the variation of $L_{+}$ in space is smooth, indicating the same resolution of physically equivalent regions without oscillations of $L_{+}$. It may be seen that ${ϵ}_{+}$ is close to unity except very close to the wall. The use of approximation ${ϵ}_{+}=1$ away from the walls results in $L_{+}=k_{+}^{3/2}$. Correspondingly, we see related variations of $L_{+}$ and $k_{+}$. Overall, the resolution features seen here are found to be very similar to corresponding features observed for periodic hill flows [20,26].

4.2. Flow Simulations

Figure 7 presents the influence of the mesh resolution on CES-KOS predictions of velocity profiles. It is observed that grid refinement can significantly affect the local variations of mean velocity profiles. At the initial station point ($x/c=0.65$), which is located just upstream of the separation region, the CES-KOS simulation using the $G_1$ grid fails to accurately estimate the velocity profile. This discrepancy arises from an underestimation of the adverse pressure gradient due to the cumulative blockage effect in the flow region. As a result, the boundary layer within the $G_1$ grid domain exhibits a considerable momentum, leading to early separation and, consequently, a larger recirculation zone. However, as the mesh density increases as given for $G_2$, $G_3$, and $G_4$ grids, the prediction of velocity profile improves significantly. Particularly, for the well-resolving $G_3$ and $G_4$ cases, the predicted velocity profile at the axial point $x/c=0.65$ closely matches the experimental result. The CES-KOS model used on $G_4$ effectively captures the momentum of the attached flow immediately after separation, indicating the importance of sufficiently fine spatial mesh resolutions along the thin shear layer for the accurate estimation of the velocity profile and identification of the reattachment point. Figure 8 shows the overall flow pattern depending on $G_1$ and $G_3$ grids in terms of streamline plots. This figure confirms the larger recirculation zone obtained by using the $G_1$ grid. Table 3 summarizes the separation and recirculation zone characteristics under different grid variations. Clearly, a significant reduction in the length of the separation bubble is observed as the grid becomes more refined. The separation bubble is quantified by

$${\delta }_{bubble}=\frac{1}{0.435}\left[{(x/c)}_{rea}-{(x/c)}_{sep}-0.435\right].$$

(3)

Here, ${(x/c)}_{rea}$ and ${(x/c)}_{sep}$ refer to reattachment and separation points, respectively, and 0.435 is the experimental bubble length.

The performance of the CES-KOS model is further evaluated by comparisons with measured Reynolds stresses, pressure and skin friction coefficients. Figure 9 shows such comparisons with the baseline experiment at axial points $x/c=(0.65,0.8)$ for the four grids considered. The modeled stress contributions are relatively small. For the fine $G_3$ and $G_4$ grids, the results show accurate CES-KOS predictions at $x/c=0.65$. At $x/c=0.8$, the CES-KOS implied stresses tend to overestimate the measured stresses, especially near the lower wall. Similar features were observed for corresponding periodic hill flow simulations [20]. In regard to the latter simulations, such discrepancies can be attributed to the fact that the experiments underestimated Reynolds stress components.

The last column in Figure 9 shows the influence of mesh refinements on $C_p$ and $C_f$ distributions compared to experimental data [33]. Here, the coefficients of pressure and skin friction [36] are calculated by

$$\begin{array}{c}\hfill C_p=\frac{p-p_{ref}}{(1/2){\rho }_{ref}{U}_{ref}^2},C_f=\frac{{\tau }_{wall}}{(1/2){\rho }_{ref}{U}_{ref}^2},\end{array}$$

(4)

where ${\tau }_{wall}$ is the viscous wall shear stress calculated from the wall surface (boundary). $C_p$ profiles demonstrate varying characteristics across different flow regions. Initially, in the range $0\le x/c<0.65$, the turbulent flow progresses towards the base of the hump model on the upstream side, resulting in a high-pressure zone. As the flow passes over the hump and moves downstream, the curvature of the streamlines leads to a linear decrease in $C_p$. The lowest surface pressure is observed at the origin of the hump (at $x/c=0$), followed by flow separation indicated by the inflection in $C_p$ profiles near $x/c=0.665\pm 0.005$. Reattachment occurs between $x/c=1.095$ and $1.105$, where the surface pressure gradually recovers to approximately the freestream value. Figure 9 shows that the estimated $C_p$ peak deviates from the experimental value by approximately $5\%$ for the $G_1$ grid and $3\%$ for the $G_2$ grid. This discrepancy is likely influenced by the applied grid type and the presence of end plates, which act as a barrier. However, refining the computational grid enhances the agreement between the experiment and simulation, with $G_3$ and $G_4$ grids closely resembling the suction peak profile within the separated zone. The CES-KOS simulation using $G_3$ and $G_4$ grids accurately captures the separation region and final pressure recovery, albeit slightly overestimating the pressure coefficient immediately after separation. Although the predicted reattachment point in the recirculation zone is slightly smaller than the experimental value, it does not significantly impact the size of the separation bubble. Similar findings have been reported in previous studies. The applied CES-KOS model provides better estimates of the pressure coefficient than the unsteady RANS model used by Capizzano et al. [39]. $C_f$ plots exhibit overall agreement across all simulations, with $G_3$ and $G_4$ grids achieving a higher level of accuracy. Thus, quality mesh density as given by $G_3$ and $G_4$ grids greatly improves both primary and secondary suction peak predictions, particularly within the recirculation zone where delayed flow reattachment leads to a slight downstream shift.

5. CES vs. Resolution Imposing Methods

The conceptual advantage of CES methods is the model’s ability to properly respond to variations of flow resolution. However, what is the relevance of this advantage to the performance of simulations? After preparing the discussion of this question by corresponding comparisons with RANS results in the next subsection, the latter question is addressed via comparisons of CES with PANS and PITM predictions. Instead of the CES setting $R=L_{+}^2$ (see Table 1), PANS applies a prescribed constant value of R. Alternatively, within PITM, an R setting is applied that is supposed to reflect $R=k_{+}$. More specifically, $R=C_{\Delta }{\Delta }_{+}^{2/3}$ is applied, where ${\Delta }_{+}=\Delta /L_{tot}$ and $C_{\Delta }=1.06$ [5,20]. Both PANS and PITM impose flow resolution, i.e., these methods are well appropriate to see the difference to the CES concept. In contrast to CES and PANS, the PITM concept involves filter width $\Delta$. For the latter, volume filter width $\Delta ={(\Delta x\Delta y\Delta z)}^{1/3}$ is applied. These comparisons are complemented in the third subsection by comparisons with one of the most popular hybrid RANS-LES methods: DES. The latter concept is based on the idea to impose LES resolution away from the wall to enable, e.g., much better simulations of separated turbulent flows.

5.1. CES vs. RANS

Two-dimensional steady RANS simulations (simply referred to as RANS in this subsection) are performed with the CES-KOS model by setting $L_{+}=1$. Comparisons of streamlines obtained by CES-KOS using the $G_3$ grid and RANS are shown in Figure 10 for $Re=936$ K. It may be seen that the RANS simulation significantly overpredicts the separation bubble. This is consistent with the known performance of standard RANS models, which have been reported to overestimate the separation bubble length by approximately $35\%$ [40], regardless of the mesh distribution and quality. The overestimation is caused by a low turbulence level in the separated shear layer [39,40,41]. In particular, Figure 10a,b correspond to a $0.3\%$ and a $16.1\%$ increase in the cross-sectional area of the separated bubble, respectively.

In the top row of Figure 11, we compare velocity profiles predicted by the CES-KOS model and the corresponding RANS model. It may be seen that the RANS simulation underestimates the flow velocity at the onset of separation and it predicts a delayed reattachment. The difference of accuracy is about $8\%$ at $x/c=0.65$ and is less than $5\%$ at $x/c=1.2$. Beyond $x/c=1.2$, the prediction remains consistent and comparable. We note that the CES-KOS model provides velocity profiles that deviate from the corresponding experimental data by no more than $3\%$ at every position. Figure 11 also compares streamwise stresses and Reynolds shear stresses. The RANS model significantly overestimates these stresses, particularly near the hump’s wall. On the contrary, the CES-KOS model offers a more precise alignment with measured normal and Reynolds shear stresses.

The accuracy of CES-KOS and RANS models in regard to predictions of pressure ($C_p$) and skin friction ($C_f$) coefficients is shown in the last column of Figure 11. In the reattachment zone, the RANS model fails to accurately capture the sudden decrease in the secondary peak of turbulent kinetic energy and the corresponding pressure coefficient profile. The skin friction plot illustrates that the RANS model predicts a significantly larger recirculation zone and late recovery due to low turbulent mixing in the freestream. In contrast, the CES-KOS model provides $C_f$ predictions that closely align with the measured data, exhibiting a notable increase in turbulent energy towards the end of the recirculation zone.

Table 5. Separation and reattachment locations obtained by the listed models. Notation–$G_3$ refers to simulations made by combining the simulation model with the corresponding settings of these methods.

Methods	Separation Location ($\mathbf{x}/\mathbf{c}$)	Reattachment Location ($\mathbf{x}/\mathbf{c}$)	Bubble Length $\mathbf{\Delta }\mathbf{x}/\mathbf{c}$	Error in Bubble Length Prediction (%)
Exp. [33]	$0.665\pm 0.005$	$1.10\pm 0.005$	0.435	–
CES-KOS–$G_3$	0.6737	1.1100	0.4363	0.3
RANS–$G_3$	0.645	1.150	0.505	16.1
RANS (Lardeau and Billard [42])	–	1.188–1.305	–	–
URANS (Capizzano et al. [39])	–	1.25	–	–
PANS-$G_3$	0.5833	1.2833	0.70	60.9
PITM-$G_3$	0.6301	1.2166	0.5865	34.8
DES (Probst et al. [43])	0.656	1.143	0.487	11.95
SST-DES (Ren et al. [44])	0.654	1.12	0.466	7.13

compares the separation and reattachment points obtained by CES-KOS and RANS models to the measured values, where other RANS results are included. Corresponding characteristics are also presented in Table 5. The CES-KOS model predicts a separation bubble length that closely matches the measured length, with separation identified at $x/c=0.6737$, slightly downstream of the experimental points ($x/c=0.665\pm 0.005$). Furthermore, the CES-KOS model accurately predicts reattachment in agreement with the measurements, while the RANS model estimates a somewhat higher surface pressure upstream of separation. The estimated variation in separation bubble length is approximately $0.3\%$ for the CES-KOS model and $16.1\%$ for the RANS model.

5.2. CES vs. PANS and PITM

Let us switch now to corresponding PANS and PITM simulation results, which were obtained by using corresponding PANS and PITM settings in the model applied. In regard to periodic hill flow simulations, it was found that PANS and PITM concepts to prescribe a desired flow resolution are hardly effective: such prescribed flow resolution does not really control the actual flow resolution measured by $L_{+}$ or $k_{+}$ [20]. On top of that, an unphysical behavior of flow resolution parameters was found close to boundaries. Corresponding features were found in regard to the NASA hump flow simulations presented here. An example is given for the PANS concept in Figure 12. The actual flow resolution measured by $k_{+}$ is found to be well above or below the prescribed flow resolution R, and the near-wall behavior of $k_{+}$ shows major differences from the prescribed R. The RHS plot of Figure 12 shows that $k_{+}$ is clearly affected by the grid in contrast to the idea of constant flow resolution.

Figure 13 presents the mean velocity streamlines obtained from CES-KOS, PANS, and PITM simulations on the $G_3$ grid. The CES-KOS model predicts a separation bubble approximately $0.3\%$ larger than the experiments. On the other hand, PITM and PANS predict separation bubbles $34.8\%$ and $60.9\%$ larger than the experiments, respectively; see Table 5.

In Figure 14, normalized velocity profiles are compared with the experimental results at two positions. PANS and PITM models properly predict streamwise velocity above the boundary shear layer region. However, they fail to accurately predict the velocities near the lower wall of the hump within the boundary layer region. In contrast, CES-KOS profiles provide reliable predictions of the normalized velocity profile both in the separation zone and after reattachment. Unlike PANS and PITM, the CES-KOS model accurately captures the rapid fluctuations of the shear layer and the mixing of the wake at the center of separation and downstream of reattachment, as well as the swift transition to a turbulent boundary layer. It shows significant improvements over PANS and PITM.

Figure 14 also shows the variations of the pressure and skin friction coefficients for the three different hybrid models. All employed models accurately estimate the mean pressure upstream of the hump. However, notable differences lie in the strength and location of the leeward low-pressure gradient, which directly influences the location of the separation region. PANS and PITM models predict significantly lower streamwise pressure gradients on the downstream side of the hump, leading to a delayed pressure recovery and a larger separation bubble. Conversely, the CES-KOS model predicts a moderate recovery of $C_p$ and a smaller separation bubble. Among all the models considered, the PANS model predicts the earliest separation on the hump surface and the largest recirculation zone. $C_f$ plots show remarkable performance deficiencies of both PANS and PITM models. In summary, the CES-KOS model provides the most accurate prediction of the reattachment zone, with reattachment occurring just beyond the base of the hump at $x/c=1.1100$.

A look at Figure 15, which shows separation and reattachment point predictions for the methods listed in Table 5, shows something interesting: PANS and PITM methods produce results that are not as good as RANS simulations, which is in contradiction to the idea of the design of hybrid RANS-LES methods.

5.3. CES vs. DES

We continue the comparison of CES results by comparing them with the results of DES (one of the most popular hybrid RANS-LES methods) obtained by Probst et al. [43]. DES hybrid RANS-LES are known to sensitively depend on simulation settings; performance shortcomings need to be often compensated by flow-dependent empirical adjustments. In the following, comparisons are presented with zonal improved delayed DES (IDDES) and a non-zonal approach based on a shear layer-adapted (SLA) definition of the subgrid length scale (the filter size $\Delta$). The zonal IDDES method is based on the shear stress transport (SST) model. It applies a synthetic turbulence generator at a (flow-dependent) RANS-LES interface positioned at $x/c=0.5$; therefore, this approach is referred to as zonal. Instead of using extra effort for a precursor simulation to provide proper inflow conditions, the idea of this approach is to introduce synthetic turbulence close to the separation point. In non-zonal DES approaches, the delay of transition to developed 3D turbulence occurs because instability is “blocked” by an excessive level of eddy viscosity in the initial region of the shear layers, which results from the convection of eddy viscosity from the attached upstream boundary layer treated by RANS [43]. In the SLA-DES version, this problem is addressed by empirical modifications of the subgrid length scale.

Figure 16 shows mean velocity, Reynolds shear stress, and skin friction coefficient $C_f$ and pressure coefficient $C_p$ in comparison with IDDES and DDES + ${\Delta }_{SLA}$ (the latter is referred to simply as DDES below) profiles from Probst et al. [43]. The mean velocity predictions are similar and in agreement with the experimental data. Figure 16 also contrasts the total shear stress profiles obtained by zonal and non-zonal DES approaches with corresponding CES-KOS profiles. It may be seen that DDES results at $x/c=0.65$ and IDDES results at $x/c=0.8$ reveal obvious shortcomings in comparisons with experimental data. The final two plots in the second row of Figure 16 show the profiles of skin friction and pressure coefficients compared to the experimental results. Evidently, all the hybrid models exhibit a similar pattern, an accurate capturing of the underlying data trends. Nevertheless, upon closer examination, slight variations become apparent between the profiles of the applied models. While the zonal IDDES approach induces a brief yet noticeable decrease in skin friction immediately downstream of the separation point, the DDES model produces a more sustained and continuous evolution of the skin friction coefficient beyond that region, ultimately matching the downstream experimental data with greater accuracy. The CES-KOS model demonstrates an enhanced accuracy in determining the separation and reattachment points in the separation region, as evidenced in Table 5. In the post-reattachment region, the DES models performed poorly, resulting in an overestimation of both the bubble length and the skin friction coefficient profile. Corresponding separation and reattachment characteristics are shown in Figure 15 in comparison with characteristics of other methods.

6. CES vs. LES Methods

We continue with the comparison of CES-KOS with LES methods, in particular LES on coarse grids (which are known to suffer from issues like significant overestimation of turbulent kinetic energy), WMLES, and WRLES. These comparisons are presented for the $Re=936$ K case by using, basically, $G_3$ and $G_4$ grids. The LES model applied is the constant coefficient Smagorinsky model. Implicit LES (iLES) of Avdis et al. [45] (involving $9.4$ M grid points) is partially involved. The deficiencies of both the constant coefficient Smagorinsky LES model and iLES are known, but they are widely used in applications. They are involved here to demonstrate their shortcomings for the flow considered.

6.1. CES vs. LES

Figure 17 and Figure 18 provide an immediate impression of differences between CES-KOS and LES. In particular, Figure 17 shows snapshots of velocity contours on the very coarse $G_1$ grid. It may be seen that the CES-KOS model generates much more fine-scale structures compared to the LES model, which acts merely in an almost RANS regime. Figure 18 shows corresponding velocity streamlines of the two models considered on the $G_3$ grid. The LES model produces a much smaller recirculation zone and an elongated separation bubble, which arises from the larger dissipation rate and subsequent fluid motion towards the wall. These effects lead to an early reattachment downstream, contributing to the observed difference in behavior between the two models.

The first row of Figure 19 shows velocity profiles at axial positions $x/c=(0.65,0.8,1.2)$ in the pre-separation, mid-separation, and post-separation regions. Both CES-KOS and LES yield similar velocity profiles in the attached region just before flow separation at $x/c=0.65$. The second row of Figure 19 shows profiles of Reynolds streamwise normal stress and turbulent shear stress at $x/c=0.8$. The LES model significantly underestimates the resolved turbulent shear stress due to the insufficient flow resolution. In addition, the LES model overestimates the intensity of Reynolds streamwise normal stress. This increase in normal stress corresponds to a significant decrease in the length of the separation bubble, and it contributes to enhanced momentum transport and early reattachment.

The last two plots of Figure 19 compare the profiles of the pressure and skin friction coefficients obtained on $G_3$ by also involving iLES of Avdis et al. [45]. All the models under consideration predict well the $C_p$ values on the upstream side ($x/c<0$) and the downstream side ($x/c>0.90$) of reattachment. The CES model accurately captures the minimum surface pressure at $x/c=0$, while the LES models consistently underestimate the peak negative pressure. LES and iLES indicate early separation occurring upstream at $x/c=0.6710$ and $x/c=0.658$, respectively. This results in a reduced separation bubble length, as indicated by the plateau in the $C_p$ plot. The underprediction of peak pressure and early separation is likely due to modeled stress depletion in this region, where the average cell size of the applied mesh is small. The CES-KOS model effectively deals with these challenges, providing reliable predictions of peak negative pressure and reattachment location, albeit with a slight delay in separation. Figure 19 also presents a corresponding comparison of skin friction coefficient profiles. The level of accuracy of these methods varies. The CES-KOS model accurately captures the observed peak in the skin friction profile, particularly in the region of $0.1<x/c<0.2$. This plateau likely results from relaminarization due to a favorable pressure gradient in that region. In contrast, LES overshoots the peak and exhibits early reattachment, which is attributed to an excessive eddy viscosity in the boundary shear layer and early transition induced by modeled stress depletion. Nevertheless, every model properly predicts the size of the separation bubble.

6.2. CES vs. WMLES, WRLES

Figure 20 presents a comparison of CES-KOS with results obtained by WMLES [16] and WRLES [46]. In particular, streamwise velocity profiles, streamwise Reynolds normal stress and shear stress profiles are shown. In regard to mean velocity, we see that the CES-KOS model and WRLES perform very well. Prediction deficiencies are observed with respect to WMLES predictions, particularly within the boundary layer region at all positions. The WMLES near-wall velocity gradient deviates from experimental values. This discrepancy may stem from the neglect of streamwise gradients ($du/dx$) in the equilibrium wall model applied, which (due to continuity constraints in nearly incompressible flows) could contribute to inaccuracies in the overall velocity profile. With respect to the stresses shown in Figure 20, we observe distinct differences between the predictions of methods considered. CES-KOS and WRLES provide realistic predictions. On the other hand (like RANS), WMLES significantly overpredicts Reynolds streamwise normal stress (by a factor of at least 3) and underpredicts Reynolds shear stress at $x/c=0.65$.

Figure 21 illustrates that all methods exhibit a reasonable agreement with the experimental pressure coefficient profiles. Predictions from WRLES match the experimental measurement profile downstream and the model is capable of mimicking the dominant features of the flow. However, within the reattachment region, the second wall pressure peak is underpredicted by WRLES compared to CES-KOS and WMLES. Figure 21 also displays the mean skin friction coefficient obtained by CES-KOS, WMLES, and WRLES simulations, demonstrating their agreement with experimental values. In the separation zone, from $0\le x/c\le 0.65$, WRLES underpredicts the skin friction coefficient, while WMLES overestimates the actual peak. In regard to post-reattachment, however, the $C_f$ profiles of WRLES and CES-KOS match relatively well, despite using different frameworks, mesh sizes, and grid resolutions.

Overall, WRLES and CES-KOS models provide the most accurate estimates of the separation bubble length, with reattachment occurring around $x/c=1.095$ and $x/c=1.1100$, respectively, slightly downstream of the hump (see

Table 6. Separation and reattachment locations obtained by the listed models.

Methods	Separation Location ($\mathbf{x}/\mathbf{c}$)	Reattachment Location ($\mathbf{x}/\mathbf{c}$)	Bubble Length $\Delta \mathbf{x}/\mathbf{c}$	Error in Bubble Length Prediction (%)
Exp. [33]	$0.665\pm 0.005$	$1.10\pm 0.005$	0.435	–
CES-KOS–$G_3$	0.6737	1.1100	0.4363	0.3
LES-$G_3$	0.6710	1.1018	0.4308	1.7
iLES (Avdis et al. [45])	0.658	1.079	0.421	3.2
WRLES (Uzun and Malik [34,46])	0.659	1.095	0.436	0.23
WMLES (Iyer and Malik [16])	0.637–0.655	1.035–1.105	0.398– 0.45	3.4–8.5

).

6.3. Cost Scalings

Computational cost differences of methods considered are addressed next by including DES in addition to CES-KOS, WMLES, and WRLES. By following Mokhtarpoor et al. [26], we consider two main cost factors, the number, N, of grid points applied and $N_t$, the number of time steps performed. Hence, simulation cost C can be calculated by

$$C=N{N}_t=TN/\Delta t.$$

(5)

Here, we apply $T=N_t\Delta t$, where T is the constant total physical simulation time and $\Delta t$ is the simulation time step. Both N and $\Delta t$ are known to vary with $Re$ according to $N={\alpha }_{1}R{e}^{{\beta }_1}$, $\Delta t={\alpha }_{2}R{e}^{-{\beta }_2}$, where ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$, and ${\beta }_2$ are constants [26].

For the flow considered, the dependence of N and $\Delta t$ can only be assessed for the CES-KOS model. For the given data, $N=(0.96,1.7,3.9)$ M and $\Delta t=(1,0.5,0.25){10}^{-5}$ for $Re=(0.936,5,10)$ M, respectively, the least-squares errors of N and $\Delta t$ of the CES-KOS model are minimized by the cost scalings given in

Table 7. Scaling factors and total cost estimates ($N/\Delta t$) for WRLES [47], WMLES [47], DES [15,43,48], and CES-KOS.

Methods	N	$\mathbf{\Delta }\mathit{t}$	$\mathit{C}/\mathit{T}=\mathit{N}/\mathbf{\Delta }\mathit{t}$
KOS	$N=492R{e}^{0.545}$	$\Delta t=0.0213R{e}^{-0.554}$	$N/\Delta t=$ 23,097 $R{e}^{1.099}$
WRLES	$N=0.00164R{e}^{1.86}$	$\Delta t=0.0393R{e}^{-0.857}$	$N/\Delta t=0.0417R{e}^{2.72}$
WMLES	$N=4.7R{e}^{1.0}$	$\Delta t=1.886\times {10}^{-6}R{e}^{-0.143}$	$N/\Delta t=2.492\times {10}^{6}R{e}^{1.143}$
DES	$N=2895R{e}^{0.545}$	$\Delta t=0.04R{e}^{-0.554}$	$N/\Delta t=$ 72,375 $R{e}^{1.099}$

. The latter scalings agree well with the scalings for the unified hybrid RANS-LES model presented by Mokhtarpoor et al. [26].

To enable a comparison with corresponding computational cost of WRLES and WMLES, we follow Yang et al. [47]. The Reynolds number scaling reported in Table 7 follows from Yang et al. [47]: their $R{e}_{L_x}$ corresponds to $Re$ considered here, and the $\Delta t$ scaling arises from the inverse scaling of $N_t$ considered by Yang et al. [47]. The combination with corresponding N and $\Delta t$ used for the $Re=0.936$ M case ($(N=210M,\Delta t=0.03\times {10}^{-5})$ for WRLES and $(N=4.4M,\Delta t=0.0264\times {10}^{-5})$ for WMLES) enables the calculation of factors of $Re$ scalings; see Table 7.

To the authors’ knowledge, there are no analyses of the scaling of computational cost of DES with $Re$ available so far. However, from a computational viewpoint, DES is very similar to unified RANS-LES methods characterized by $N\sim R{e}^{0.55}$ [26], which is equivalent to the CES-KOS N scaling observed here. In regard to DES $\Delta t$ scaling, there is no reason for the DES scaling to differ from the corresponding CES-KOS scaling. Thus, we apply the CES-KOS $Re$ scaling for DES, too. For the $Re=936$ K case considered, we know the DES cost: the latter is characterized by $N=5.2$ M and $\Delta t=2\times {10}^{-5}$. The DES cost data are used to determine the $Re$ scaling factors in Table 7.

An illustration of these cost scalings is given in Figure 22. Relevant observations are the following ones. Costs of $C/T={10}^{15}$ are currently feasible, as demonstrated by WRLES. A usual $Re$ range of many complex turbulent flows is ${10}^6\le Re\le {10}^9$. Hence, several decades are probably required before WRLES can be made applicable to many relevant turbulent flows. On the contrary, WMLES, DES, and CRES-KOS are applicable to ${10}^6\le Re\le {10}^9$. However, WMLES costs are more than two orders of magnitude higher than the corresponding CES-KOS costs, and CES-KOS costs are only 32% of the corresponding DES costs. The latter can be attributed to the fact that CES-KOS was derived by minimizing the hybridization error.

7. Very High-Re Simulations

Usually applied computational methods like RANS, LES, WMLES, and DES face significant problems with respect to predictions of high-$Re$ flows for which validation data are unavailable. The application of LES is restricted by resolution requirements that need to be satisfied; RANS methods are known to be unreliable, and hybrid RANS-LES like WMLES and DES suffer from the uncertainty of their predictions caused by unknown simulation settings needed to ensure agreement with validation data. In contrast, the CES-KOS model applied here incorporates a physics-based mechanism to respond to changes in resolved motion as given by an increasing $Re$. It is, therefore, of interest to study implications of increasing $Re$ on the flow structure. This is achieved next by extending the previous simulations at $Re=936$ K by simulations at $Re=5$ M and $Re=10$ M. The grids applied and time step settings for these simulations are described above.

7.1. $Re$ Effects on Flow Resolution

Figure 23 shows the profiles of resolution indicators $L_{+}$, $k_{+}$, and ${ϵ}_{+}$ at $x/c=1.2$ for the $Re=5$ M and $Re=10$ M cases. Relevant observations are the following ones. First, a coarser grid and a higher $Re$ increase the amount of modeled motion reflected by higher $L_{+}$ and $k_{+}$ values, i.e., the mode interaction mechanism does work correctly. In particular for the $Re=10$ M case, we have simulations that are truly in between RANS and LES regimes, i.e., far from well-resolved regimes. Second, very similar to the $Re=936$ K case, ${ϵ}_{+}$ is close to unity except the cases when it is very close to the walls. It is of interest to see slight deviations from unity for the $Re=10$ M case. Third, very similar to the experience obtained from periodic hill flow simulations [20], we see smooth variations of $L_{+}$, $k_{+}$, and ${ϵ}_{+}$ in the space without oscillations, i.e., physically equivalent regions are equivalently resolved. These features are fully consistent with corresponding features at $Re=936$ K seen in Figure 6.

7.2. $Re$ Effects on Flow Structures

Before considering the $Re$ effect on flow structures, it is of interest to clarify the influence of the grid for high-$Re$ cases. The latter question is addressed in terms of

Table 8. CES results: Summary of CES-KOS separation and reattachment points for high-$Re$ cases.

Cases	Re	Separation Location ($\mathbf{x}/\mathbf{c}$)	Reattachment ($\mathbf{x}/\mathbf{c}$)	Bubble Length $\mathbf{\Delta }\mathbf{x}/\mathbf{c}$	Error in Bubble Length Prediction (%)
Exp. [33]	$Re=936$ K	$0.665\pm 0.005$	$1.10\pm 0.005$	0.4350	–
	$Re=936$ K	0.6737	1.1100	0.4363	0.3
$G_3$	$Re=5$ M	0.6713	1.2167	0.5454	–
	$Re=10$ M	0.6641	1.2165	0.5524	–
	$Re=936$ K	0.6680	1.100	0.4320	0.7
$G_4$	$Re=5$ M	0.6725	1.2025	0.5300	–
	$Re=10$ M	0.6812	1.2245	0.5433	–

, which shows differences between CES-KOS $G_3$ and $G_4$ grid results with respect to separation and reattachment characteristics. Increasing $Re$ leads to a marginal increase in the separation region, consistent with the flattening effect of $Re$ on the size of the separation bubble, as previously reported by Uzun and Malik [34] and corroborated by Rumsey et al. [40] and Seifert and Pack [32]. We see that $G_4$ grid simulations imply slightly more consistent results in regard to predictions of separation an reattachment points, but these differences are indeed very small. The latter supports the suitability of using the $G_3$ grid for these $Re$ cases. The difference in separation bubble lengths between $Re=5$ M and $Re=10$ M cases is only $2.5\%$, indicating an almost asymptotic $Re$ regime.

The influence of increasing $Re=$ (936 K, 5 M, 10 M) on the flow structure was studied based on $G_3$ grid simulations. Figure 24 illustrates that higher $Re$ values have a negligible influence on the velocity profile of the flow (it is worth noting that the experimental data shown here as a reference do only apply to the $Re=936$ K case). Hence, the mean velocity profiles also indicate an asymptotic flow behavior.

Correspondingly, Figure 25 illustrates the effect of $Re$ on pressure and skin friction coefficients. While the overall wall pressure distribution remains largely unaffected, increasing $Re$ minimally reduces the primary maximum peak due to energized flow dynamics and increased turbulence. This is consistent with findings by Vatsa et al. [40]. The secondary pressure peak also flattens for higher $Re$, indicating low recovery and delayed reattachment. $C_f$ profiles show more noticeable differences across $Re$. The profiles show maxima at $x/c\approx 0.1-0.65$ that increase in magnitude over the hump surface as $Re$ increases. This is due to an increased intensity of turbulent fluctuations and larger-scale structures in the flow, leading to a stronger pressure gradient and higher skin friction coefficients. $C_f$ profiles vary little in the separation region across all $Re$ cases. Although the separation regions are nearly identical, the lowest-$Re$ case experiences early reattachment at $x/c\approx 1.1100$, while the reattachment points occur at approximately $x/c\approx 1.2167$ and $x/c\approx 1.2165$ for $Re=5$ M and 10 M, respectively. Overall, while increasing $Re$ shows negligible pressure coefficient effects, the skin friction coefficient and bubble length show Reynolds number sensitivity.

7.3. Asymptotic $Re$ Regime

The question about the asymptotic flow regime is certainly of interest. It matters to our understanding to see whether there are asymptotically stable regimes of wall-bounded turbulent flows involving flow separation. More specifically, this question relates to our understanding of which geometric conditions enable an asymptotically stable flow configuration [49], what is the corresponding concrete flow structure, and which $Re$ are needed to accomplish an asymptotic flow structure. Also, to provide comprehensive support for the validity and predictive power of simulation methods, evidence is needed for simulation predictions for an extended range of $Re$, including almost asymptotic $Re$ regimes. This leads to the question of what are almost asymptotic $Re$ regimes. As shown here, the computational cost of simulation methods depends significantly on $Re$. There is the question about the applicability of WRLES and experiments to provide validation data at required high $Re$.

The prediction of reattachment points and their convergence is well appropriate to address the question about the existence of an asymptotic $Re$ regime. Existing experience for the flow considered is summarized in

Table 9. Summary of reattachment points in Figure 26 for various methods. The upper and lower parts correspond to Mach numbers of 0.1 and 0.25, respectively.

Methods	Reference	Re	Reattachment Point (${\mathit{X}}_{\mathit{r}}$)
Exp.	Greenblatt et al. [33]	936 K	$1.10\pm 0.005$
RANS	Lardeau & Billard [42]	936 K	1.3050
LES	Avdis et al. [45]	936 K	1.0790
LES (WRLES)	Uzun and Malik [46]	936 K	1.0950
Hybrid (WMLES)	Iyer and Malik [16]	936 K	1.035–1.1050
Hybrid (DES)	Probst et al. [43], Radoslav & Roger [50]	936 K	1.1520, 1.1050
Hybrid (SST-DES)	Ren et al. [44]	936 K	1.1200
Hybrid (STM)	Woodruff [51]	936 K	1.0800–1.1900
Hybrid (KOS)		(936 K, 5 M, 10 M)	(1.1100, 1.2025, 1.2255)
Exp.	Vatsa et al. [40]	4.2 M	$1.21\pm 0.05$
URANS (FUN2D)	Vatsa et al. [40]	(4.2 M, 16 M)	(1.3330, 1.3100)
RANS (TLNS3D)	Vatsa et al. [40]	(4.2 M, 16 M)	(1.2700, 1.2550)
RANS (CFL3D)	Vatsa et al. [40]	(4.2 M, 16 M)	(1.3040, 1.2830)

and visualized in Figure 26. Without paying attention to CES-KOS results that are also shown, it may be seen that existing observations do not provide a conclusive answer to the question about the asymptotic regime. The large scatter of data does hardly allow conclusions about the asymptotic trend; in particular, no conclusion can be drawn about how the reattachment points vary with $Re$. It is of interest to note that the codes applied clearly have an influence on reattachment point predictions; see the RANS-type results in the last rows of Table 9 [40].

On the other hand, CES-KOS enables a clear conclusion about how reattachment points vary with $Re$. We obtain $x_r=1.229/(1+e^{-0.3857({10}^{-6}Re+4.842)})$, which is also shown in Figure 26. The CES-KOS predictions are well supported by corresponding RANS-type predictions. Something that matters in this regard is the fact that the CES-KOS result, $x_r=1.229/(1+e^{-0.3857({10}^{-6}Re+4.842)})$, is supported over a very large range of $Re$ values. This avoids problems of corresponding curve fits supported by narrow $Re$ ranges [52], which are inapplicable to larger $Re$ [20].

We note that CES-KOS results are also in meaningful agreement with corresponding periodic hill flow studies indicating the existence of an asymptotic $Re$ regime for $Re\ge 2.6$ M [20]. A discussion of similarities of conclusions about asymptotic $Re$ regimes predicted for canonical wall-bounded turbulent flows (channel flow, pipe flow, and zero-pressure gradient turbulent boundary layers) can be found elsewhere [53,54].

8. Conclusions

The paper reported the first application of CES methods to a moderately complex, relatively high-$Re$ turbulent flow: the NASA wall-mounted hump flow. Also, for the first time, we presented comprehensive comparisons of CES methods with RANS, LES, and several hybrid RANS-LES as PANS, PITM, DES, and WMLES. The following conclusions can be drawn in regard to questions Q1–Q3 posed in the introduction.

In regard to Q1, one of the most serious limitations of existing computational methods for turbulent flows of practical relevance is their dependence on data needed for validation. Practically, this excludes reliable predictions. Usual practice of hybrid RANS-LES is a choice of simulation settings that enables a good agreement with validation data, as reported here, for example, in regard to versions of DES and WMLES. WRLES does not suffer from the latter problem, but it also needs validation data to confirm the suitability of grids.

• From a general viewpoint, the advantage of CES methods demonstrated here was an excellent simulation performance independent from adjustments of simulation settings. The latter advantage is based on the model’s ability to adjust the model to changing flow conditions (the amount of resolved motion). A grid dependence of simulation results was found if extremely coarse grids are applied.
• The implementation of a sensibilization of the model to the degree of flow resolution is not the only requirement to enable a performance as demonstrated here in regard to CES methods. An equally important condition is the mathematically correct implementation of an appropriate model sensitivity. The latter was demonstrated here via comparisons with PANS and PITM methods, demonstrating a simulation performance that is outperformed by CES methods.
• Contrary to CES, other popular simulation methods impose a certain desired flow resolution: this applies to RANS, LES, and other hybrid RANS-LES. This inherent model inflexibility generates a sensitive dependence on simulation parameter settings, for which optimal solutions cannot be found for every flow. In addition, this concept does not ensure that the desired imposed resolution is actually realized. This can imply higher computational cost (finer grids) to compensate for performance deficiencies [17].

In regard to Q2, by considering both simulation performance and computational cost, what are the implications of the results reported here in regard to simulations of turbulent flows of practical relevance? The focus of the following is at least partially on resolving simulation methods.

• Table 7 speaks a clear language in regard to the use of WRLES: for decades to come, WRLES will be inapplicable to many relevant turbulent flows. The use of WRLES also faces questions from a simulation performance view point. Figure 21, for example, shows that the use of WRLES is no guarantee for excellent simulation results. The latter can be attributed to a significant conceptual issue of WRLES: the inclusion of the filter width as length scale, which can be unphysical.
• The most popular hybrid RANS-LESes are DES and WMLES. Both depend on a variety of simulation settings, as specified here, which need to be adjusted to (often unavailable) validation data. The computational costs of DES and WMLES are well above the CES cost; see the discussion related to Table 7. Performance-wise, CES-KOS was shown to perform better than DES versions and WMLES. These differences may be attributed to the CES feature to be a minimal error simulation method.
• PANS and PITM were introduced as alternatives to DES and WMLES. Figure 15 in conjunction with other related figures shows something interesting: the use of hybrid RANS-LES is no guarantee to obtain simulation results better than corresponding RANS prediction; the opposite may be the case. The underlying conceptual issue is the idea to impose a desired flow resolution. The latter does not only not determine the actual flow resolution, but the mismatch implied can deteriorate simulation results.

In regard to Q3, the question about the asymptotic flow structure is certainly relevant. This relates to an understanding of which geometric conditions enable an asymptotically stable flow configuration, what is the corresponding concrete flow structure, and which $Re$ are needed to accomplish an asymptotic flow structure. CES methods certainly offer advantages with respect to addressing such questions based on their ability to adjust the model to the changing flow resolution. The following observations were obtained here:

• The structural effects of increased $Re$ on the mean velocity field are relatively small; there are only very minor variations. This speaks for the consideration of an $Re$ range that is relatively close to an asymptotic regime. Correspondingly, the pressure coefficient distribution is little affected by increasing $Re$.
• On the other hand, the skin friction coefficient distribution and bubble length were identified as sensitive measures of Reynolds number dependence. Both are clearly affected by higher $Re$ in response to an increased intensity of turbulent fluctuations and larger-scale structures in the flow.
• Reattachment point predictions were used to specifically ask about the potential existence of an asymptotic $Re$ regime for the flow considered. CES-KOS predictions confirm that there is an asymptotic $Re$ regime. The $Re=10$ M case is at least very close to this asymptotic regime. This conclusion is well supported by corresponding RANS predictions. This observation is in meaningful agreement with corresponding periodic hill flow studies indicating the existence of an asymptotic $Re$ regime for $Re\ge 2.6$ M [20].