## 1. Introduction

Turbulent shear flows, including attached boundary layers and separated shear layers, are important application test cases for computational fluid dynamics prediction since they are observed for a wide array of engineering processes and systems. In particular, turbulence modeling for these systems is an important aspect of the simulation that is often responsible for predictive error. Reynolds-averaged Navier-Stokes (RANS) models are known to perform relatively well for turbulent boundary layer flow [1] due to the somewhat universal nature of wall-bounded turbulence. RANS models typically do not perform as well in regions of separated flow due to the presence of adverse pressure gradients, strong three-dimensionality, shear layer reattachment, and large-scale unsteadiness [2–8]. Large eddy simulation (LES) models theoretically provide greater accuracy than RANS approaches in these regions, but their application to attached wall-bounded flows suffers from significantly increased computational expense relative to RANS. As a consequence, the LES approach is still not widely used for industrial analysis and design for high Reynolds number flows, especially those with attached boundary layers [9]. Hybrid RANS/LES (HRL) [1] approaches offer a potentially attractive alternative to RANS or LES, since they attempt to combine the advantages of both RANS and LES modeling in an optimized manner that resolves both attached and separated flows effectively. Specifically, HRL models have the potential for greater accuracy than RANS and less expense than LES. Interest in HRL methods has, therefore, grown substantially over the past several years.

Hybrid RANS-LES models are categorized as either zonal or nonzonal. For zonal models, the RANS and LES models are separately employed in selected regions of the computational domain which are determined a priori. Development of effective methods for coupling the two model types at their interface remains a challenge and is an ongoing area of research [10, 11]. Nonzonal methods are generally simpler to implement and do not require the user to decide where LES or RANS is used in a particular simulation. In a nonzonal model, the eddy viscosity adopts a value representative of a RANS model in the near-wall region and a value representative of an LES subgrid stress (SGS) model in separated regions. Detached-eddy simulation (DES) [12] is probably the most commonly used nonzonal modeling methodology. Blending between RANS and LES model types in the original DES formulation is a function of the local grid size and has been shown to suffer from inaccuracy in attached boundary layers [13]. While several ad hoc modifications have been implemented to address these limitations, the modifications often only mitigate weaknesses while not completely resolving them. For example, to address the problem of reduced levels of eddy viscosity in attached boundary layers due to premature switching to LES mode, Spalart et al. [13] introduced a modified version of the baseline DES model denoted as Delayed DES (DDES). DDES adopts a definition for length scale different from that in the original baseline DES model, as a function of local mean flow and turbulence model variables. Shur et al. [14] introduced a DES version with an additional modification—improved delayed DES (IDDES)—in an attempt to eliminate the well known “log layer mismatch” issue that arises in classical DES and wall-modeled LES (WM-LES). The IDDES model behaves similar to DDES except that it replicates a WM-LES type model in boundary layer regions when resolved turbulent fluctuations are present.

The key challenge for nonzonal HRL modeling is that of effectively specifying the transition between RANS and LES modes in the simulation. Typically, such a transition is defined based on the value of the eddy viscosity that varies between the Reynolds stress model and the SGS value. Significantly, the Reynolds stress and subgrid stress are mathematically distinct and physically different; therefore, any method that switches between the two using a single parameter (i.e., eddy viscosity) is not straightforward. This is because the Reynolds stress is based on an ensemble averaging of all turbulent scales present in the flow field, while the subgrid stress is based on (typically spatial) filtering of all turbulence scales that are not directly resolved in the simulation for a given mesh size. The difficulty in defining zonal transition using a single variable for eddy viscosity has been identified as a major weakness for commonly used HRL models [15–17]. Additionally, many of the currently used HRL models incorporate the local grid size directly as a model variable. This necessitates that great care is taken when constructing grids for HRL modeling. For the best performance, the grid should be built with foreknowledge of the expected model behavior and design of the specific grid used as the method of forcing RANS-to-LES transition in the desired locations of the domain [15].

The transition from a purely modeled RANS stress to a resolved dominating LES stress has been recognized as a major concern [15]. The problem is exacerbated if separation occurs at a clearly defined location such as a sharp point (e.g., a backward-facing step) and the separating boundary layer lacks any initial level of fluctuating turbulence content. Paterson and Peltier [16] studied the issues for RANS-to-LES transition for cases where no geometrically imposed separation point exists. The authors note that a delay in the development of resolved fluctuations stress terms is present in the RANS-to-LES transition region upstream of the separation location; therefore, the (SGS) turbulent viscosity values obtain dominance over the RANS eddy viscosity values prematurely. This effect is the well known “modeled-stress depletion” described by Spalart et al. [13]. Nitikin et al. [17] also demonstrated the difficulties inherent in calculating the appropriate level of grid resolution for the “gray region” that lies between the RANS and LES modes, particularly for RANS-to-LES transition that occurs in the wall normal direction in boundary layer flows.

Previous researchers have attempted to resolve the RANS-to-LES transition issue [13, 14, 18–21] including efforts discussed above [13, 14]. Menter et al. [18, 19] introduced the concept of scale-adaptive simulation (SAS), which is significant as it provides the potential to develop turbulence models that are applicable in both RANS and LES modes without including any explicit grid-dependence. Hamba [20] proposed that rapid variation of the filter width near the interface of the RANS and LES zones is the primary reason for the log-layer mismatch in channel flow simulations and that the problem can be resolved by using an additional filtering operation. Piomelli et al. [21] have proposed the use of a pseudo-random forcing function as a backscatter model in the interface region as a means of resolving the underlying issues of the a transition layer between RANS and LES. It is important to note that most of these prior attempts should be viewed as ad hoc modifications rather than fundamental changes to the basic modeling approach. Celik [22], in a review paper, proposes that entirely new criteria are needed in order to effectively address the RANS-to-LES transition issue in HRL models in a more fundamental way.

The dynamic hybrid RANS-LES (DHRL) modeling methodology presented in this paper was specifically developed as an attempt to resolve the aforementioned weaknesses, including explicit grid dependence, that have been documented for most previous HRL models. Furthermore, it is assumed that these issues are fundamental in nature and unlikely to be effectively addressed using ad hoc modifications. The novel features of the DHRL modeling framework are that (1) it is a general framework which allows coupling of any particular RANS model with any particular LES model; (2) it does not include mesh size as a variable in the model formulation; (3) the blending between RANS and LES modes is enforced through the assumption of continuity in total turbulence energy production; and (4) it exactly reduces to the underlying RANS model in flow regions that exhibit numerically steady-state results.

In the present paper, a detailed investigation of the DHRL model is presented using the finite-volume based commercial CFD software Ansys FLUENT^{®}. The DHRL model was implemented into FLUENT using the user-defined function (UDF) capability available with that solver. For the present cases, the DHRL framework was used to integrate Menter’s SST k-ω model [23] as the RANS component, with monotonically integrated LES (MILES) [24] as the LES model component. To evaluate the performance of the DHRL model in both attached and separated flow regions, the test cases considered were turbulent channel flow matching the DNS case of Hoyas and Jimenez [25] and backward-facing step flow matching the experimental case of Driver and Seegmiller [26]. Simulation results using the DHRL model are compared with the corresponding DNS and experimental data, and with the results of companion simulations using other RANS and HRL turbulence models available in FLUENT.

## 2. Model formulation

The DHRL model formulation is summarized in this section. Readers are referred to Bhushan and Walters [27] for a more detailed description. Other models, including DES, DDES, and SST k-ω, are also summarized since they are used to perform companion simulations, and the results are compared with DHRL. In addition, the SST k-ω model is used as the RANS component of the DHRL model, while MILES is used as the LES component. The models used here were chosen for comparison purposes because they are widely used examples of HRL and RANS, respectively. Since they have been previously thoroughly documented, only a brief description and appropriate references are provided in Subsections 2.2–2.5.

### 2.1. Dynamic hybrid RANS-LES (DHRL) methodology

The description of the DHRL model presented here focuses on single-phase, incompressible, Newtonian flow with no body forces. Extension to compressible flow or flows with gravitational effects is straightforward. First, applying an (undefined) filtering operation to the momentum equation results in:

where *u*_{i} and

The turbulent stress term must be modeled in order to close the filtered momentum equation.

Most hybrid RANS-LES models, including the commonly used DES model, use a single scalar variable in the filtered momentum equation to model the turbulent stress using the Boussinesq hypothesis. This term is denoted as the eddy viscosity, and in theory obtains a value appropriate for a modeled Reynolds stress in the RANS region (typically near the wall) and a value appropriate for a modeled subgrid stress in the LES region (typically far from the wall).

As discussed above, blending the effects of ensemble-averaged velocity fields (Reynolds stress) and spatially filtered velocity fields (subgrid stress) using a single scalar variable introduces ambiguity into the model. The DHRL modeling methodology seeks to avoid this ambiguity. The model development begins with decomposition of the velocity field in such a way that the effects of ensemble-averaged velocity fields and spatially filtered velocity fields retain a distinct separation in the transitional or “gray” zones.

First is introduced a simulation-specific decomposition for the instantaneous velocity (*u*_{i}):

where *ū*_{i} is the mean (Reynolds-averaged) velocity, *u*″_{i} is the resolved fluctuating velocity, and *u*′_{i} is the unresolved or modeled fluctuating velocity. The Reynolds-averaged velocity and resolved fluctuating velocity can be obtained directly from the simulation. The modeled fluctuating velocity is defined using the turbulent stress term. Substitution of the decomposed instantaneous velocity (*u*_{i}) in Eq. (3) into Eq. (2), along with the assumption that resolved and unresolved velocity fluctuations are uncorrelated, yields the following expression for the residual stress:

The scale similarity concept can be used to model both of the terms on the right-hand side of Eq. (4). The result is an expression for the subfilter stress term as:

The first (both parts inside parenthesis) and the second terms on the right-hand side of Eq. (5) can be modeled as linear functions of the subgrid stress (SGS) and Reynolds stress, respectively. The SGS and Reynolds stresses can be obtained from any suitable SGS and RANS model. The proportionality constants α and β can be both spatially and temporally varying, but are assumed to be complementary, such that the residual stress term can be expressed as a weighted average of both the SGS and RANS stress as follows:

To determine the value of the blending coefficient α, a secondary filtering operation is applied, conceptually similar to the method of Lilly [28] for dynamic model coefficient evaluation for subgrid stress modeling. Based on the following assumptions:

the Reynolds-averaging operation can be applied to Eq. (2) as a secondary filter and combined with Eq. (7) to yield:

The Reynolds-averaged form of Eq. (6) is combined with Eq. (9) to eliminate

(10) |

The local value of the coefficient α is determined based on the relative contribution to turbulence production by the resolved scales, the mean (Reynolds-averaged) component of the subgrid model stress, and the RANS model stress. For stability, the value of α is limited in practice to lie between 0 and 1. Eq. (10) shows that the value of α becomes zero in regions with no resolved fluctuations, i.e., numerically steady-state flow. In those regions, therefore, a pure RANS model is recovered. In regions for which turbulence production by resolved velocity fluctuations is significant, the RANS stress contribution decreases, and an LES subgrid stress contribution appears in the momentum equation. This ensures a smooth variation of turbulent production through the transition between RANS and LES. If the resolved turbulent production in any region is sufficiently large, α increases to a value of 1, and a pure LES model is recovered.

One final key aspect of the DHRL methodology is the manner in which the RANS model component is computed. In contrast to most other hybrid methods, the DHRL approach computes the RANS terms based only on the Reynolds-averaged flowfield. For stationary flows, the velocity field used to compute all RANS terms can be obtained from a time-averaging operation that runs concurrent with the simulations. Other appropriate averaging methods can be adopted for nonstationary flows. For the current study, stationary flows are considered, and the RANS model is computed using the time-averaged flowfield in all of the turbulence model terms.

### 2.2. Detached-eddy simulation (DES)

The DES model is based on the one-equation eddy-viscosity RANS model developed by Spalart and Allmaras [29] (SA), which has been used extensively in applied CFD analysis. The SA model includes one additional transport equation for an eddy-viscosity variable,

where *C*_{d} is a variable model coefficient and *d* is the distance to the nearest wall. The wall distance also appears in other model coefficients.

The DES model replaces the wall distance with an effective length scale,

where ∆ is the local characteristic mesh spacing and *C*_{DES} is a model constant. The effect of the modification is to increase destruction of

### 2.3. Delayed detached-eddy simulation (DDES)

The DDES model was proposed as a modification to the baseline DES model in order to address the issue of grid-induced premature switching to the LES mode in attached boundary layers [13]. This anomalous activation of LES mode occurs in the DES model when a highly refined grid is used. In contrast to the baseline DES model, DDES redefines the length scale such that it depends on both the local grid size as well as the eddy viscosity:

The function *f*_{d} recovers a value of zero inside the boundary layer, ensuring that the RANS mode is activated; however, the value of *f*_{d} limits to one outside the boundary layer, thus recovering the underlying DES mode. The functional form of *f*_{d} is

where the variable *r*_{d} is obtained as:

### 2.4. Shear-stress transport (SST) k-ω

The SST k-ω model proposed by Menter [23] is conceptually based on the transport of the principal shear stress and was developed to improve prediction of flows with adverse pressure gradients. It has been used for a large number of practical RANS CFD simulations of complex turbulent flows [30]. In the SST model, the eddy viscosity from the baseline k-ω model is modified within the framework of the SST model as follows:

where *F*_{2} is a blending function, *a*_{1} is a constant, and *S* represents an invariant measure of the strain rate magnitude. *F*_{2} recovers a value of one in attached boundary layers, and a value of zero in free shear layers. The model improves prediction in adverse pressure gradient boundary layers by ensuring that production of turbulent kinetic energy is larger than dissipation.

Two transport equations are adopted in the SST model, one for the turbulent kinetic energy (*k*) and one for the specific turbulence dissipation rate (*ω*):

The function *F*_{1} plays a similar role as the blending function *F*_{2}, acting as an indicator variable for near-wall and farfield regions of the flow. Near the wall, *F*_{1} = 1, and the k-ω model form is obtained. In farfield regions, *F*_{1} = 0 and the model behaves similar to a k-ε model. Refer to Ref. [23] for further model details. For the DHRL model used in the present study, the SST model was used as the RANS component.

### 2.5. Monotonically integrated large eddy simulation (MILES)

MILES refers to an approach for implicit LES first documented by Boris et al. [31] for which the subgrid stress is assumed to be accurately approximated by the dissipation inherent in the discretization error for the convective terms. For upwind-biased schemes in particular, the numerical error can be shown mathematically to be similar in form to an explicit eddy-viscosity model for such an approach [24], and several studies have been documented in the literature in which practical LES solutions have been successfully obtained. For the present study, the MILES method is used for the reference LES simulations and also as the LES component in the DHRL implementation. Two different convective discretization schemes are used (these are discussed below). Both schemes are designed to preserve monotonicity through upwinding and flux limiting, which is appropriate for a MILES approach.

## 3. Numerical method

Simulations in the present study were all run using the commercial CFD code Ansys FLUENT. The DHRL methodology was implemented using the User-Defined Function subroutines available with that solver. Other turbulence models were all available as built-in modeling options in FLUENT. All of the test cases considered are incompressible, so the segregated pressure-based solver was used with the SIMPLE pressure-velocity coupling method [32]. Second-order implicit (three-point backward difference) discretization was used for the unsteady term in all time-dependent simulations including HRL and LES. A second-order upwind, linear reconstruction scheme with a Least-Squares gradient computation and conventional slope limiting [33] was used as the baseline discretization scheme for the convective terms in the momentum equations. Mass fluxes were computed using the momentum-weighted interpolation method of Rhie and Chow [34]. For some of the test cases, a less dissipative bounded central differencing spatial discretization scheme based on the Normalized Variable Diagram [35] was used to evaluate the influence of discretization scheme on the HRL results. A second-order centered discretization was used for the pressure gradient terms, and second-order central differencing was used for all diffusion terms.

Steady RANS simulations were run until a converged solution was obtained, with convergence determined based on stabilization of all flow variables to constant values as iterations increased, and a decrease in the L2 norm of all residuals by at least six orders of magnitude less from the initial value. For the hybrid RANS-LES (HRL) models, the converged steady-state RANS results were used as the initial conditions. The HRL model simulations were time-dependent and were run until a statistically steady-state flow field was obtained, as judged by time-averaged variables obtaining constant values with increasing time steps. The time-step size for unsteady simulations was selected to maintain maximum convective CFL number of approximately one, based on the smallest streamwise mesh dimension and the freestream velocity. For the time-dependent cases, the L2 norm of all residuals was verified to be reduced by at least three orders of magnitude during each time step. For both attached and separated flow test cases, an additional HRL simulation was run with the time-step size reduced by half, and results were compared to the default time step cases. No significant difference was observed for the resulting statistical quantities, so it was concluded that the use of a convective CFL number of approximately one is appropriate.

## 4. Attached turbulent flow test case

The selected attached flow test case is a fully developed turbulent channel flow, corresponding to the experimental study of Hoyas and Jimenez [25]. Results for this case using the DHRL model were previously obtained using a pseudo-spectral solver and a combination of the Spalart-Allmaras one-equation model [29] for the RANS component and the dynamic Smagorinsky model [28] for the LES component [27]. This paper investigates in more detail DHRL simulation using a general-purpose finite-volume flow solver.

The simulations were performed for the case of Reynolds number equal to 2003, where

*u*_{τ} is the wall friction velocity, *δ* is the channel half-height, and *ν* is the kinematic viscosity. The domain extent was 2*πδ* in the streamwise and *πδ* in the spanwise directions. Three different structured meshes were investigated, with streamwise × spanwise × wall-normal cell counts of 32 × 24 × 36 (Coarse), 64 × 48 × 72 (Medium), and 128 × 96 × 144 (Fine). The grid spacing was uniform in the streamwise and spanwise directions. Stretching was employed in the wall-normal direction so that the first near-wall cell centroid was located with *y*^{+} of 2, 1, and 0.5 for the Coarse, Medium, and Fine grids, respectively, and cells located at the channel centerline had aspect ratios of approximately unity.

Hybrid RANS-LES simulations were performed with three different models—the DHRL method described above, the original detached-eddy simulation (DES) model based on the Spalart-Allmaras RANS model [12], and the more recently proposed delayed detached-eddy simulation (DDES) model [13]. The DHRL model combined the k-ω SST model [23] as the RANS component and monotonically integrated LES (MILES) [24] as the LES component. For comparison purposes, steady RANS simulations were obtained using k-ω SST model alone, and LES simulations were obtained using the MILES method. Results for all models were evaluated based on comparison with the available DNS data [25]. Quantitative comparisons were made primarily for mean velocity and turbulent kinetic energy. Of particular interest was the effect of grid resolution level and discretization scheme for the different modeling approaches investigated. Unless otherwise noted, all results shown are plane-averaged, though all statistical quantities showed little variation in the streamwise and spanwise directions once the simulation had run sufficiently long for a stationary state to be reached.

**Figure 1** shows the mean velocity profile predicted by each of the models on the Coarse, Medium, and Fine grids. Note that the results with the SST model are only shown for the Coarse grid, since no noticeable differences were observed between that and the two finer grids. The MILES results are only shown for the Medium and Fine grids, since the results on the Coarse grids were comparable to laminar flow, and dramatically overpredicted the velocity in the middle portion of the channel. In fact, it is apparent from **Figure 1 (b)** that even results on the Medium grid are qualitatively incorrect, and only on the Fine grid does the characteristic log-layer behavior begin to be observed, though mean velocity of it is still significantly overpredicted.

For the Coarse grid, the SST, DDES, and DHRL models all show qualitative agreement with the DNS data, with the DDES model showing the best agreement, and both the SST and DHRL models showing a similar, slight overprediction in the log-layer region. The DES model shows the characteristic overshoot that has previously been noted in the literature, which is due to the performance of the model in the transition region between the (near-wall) RANS zone and the far-field (LES) region in the middle of the channel. As the grid is refined (**Figure 1 (b, c)**), the DES results agree more closely with the DNS data, and for the Fine grid, the log-layer mismatch is much smaller, though still apparent.

For all three grid sizes, the DDES model obtained a numerically steady-state result, effectively operating as a pure RANS model. This is not surprising, since the goal of the Delayed DES formulation is to inhibit the development of LES behavior in attached flow regions, while allowing its development in separated flow regions. As a consequence, the DDES mean velocity results were limited to a Spalart-Allmaras model equivalent, which is in good agreement with the reference data.

The predicted turbulent kinetic energy (*k*) profiles for the hybrid RANS-LES models are shown in **Figure 2**. Since the DDES model returned a steady RANS solution, no result is plotted. For the DHRL model, both resolved turbulence and total (resolved + modeled) turbulence is shown, where the total *k* is computed as:

For both the DES and DHRL models, the overall trends are correctly predicted, with a peak near the wall and decay to a minimum value at the channel centerline. Both models tend to underpredict the peak value of *k* near the wall as well as predicting the location of the peak to be too far from the wall, and for the Fine grid case show an underprediction throughout the channel. The DES result shows nontrivial sensitivity to the mesh refinement level, and surprisingly shows better agreement throughout most of the channel on the Coarse grid than on the Fine grid. The DHRL model shows relatively little dependence on the mesh resolution level, except for an improvement in the near-wall prediction as the grid is refined (as does the DES model).

The mesh sensitivity is further highlighted in **Figure 3**, which shows the change in the mean velocity profiles for each model as the mesh is refined from 27,648 total cells (Coarse) to 1,769,472 cells (Fine). For the DES model (**Figure 3 (a)**), the improvement with increasing mesh resolution is apparent, as the characteristic overshoot becomes smaller. The DDES model is not shown since it yielded a steady-state result that was insensitive to mesh refinement level for the grids considered here. The DHRL model (**Figure 3 (b)**) shows very little sensitivity to grid refinement level and shows better agreement with the DNS data than the DES model for all refinement levels. Also shown in **Figure 3 (b)** are the steady RANS results from the SST model, which is the RANS component for the DHRL hybridization. The results from DHRL agree quite closely with the SST for all mesh refinement levels, despite the fact that a significant portion of the turbulent kinetic energy is resolved as unsteady velocity fluctuations (**Figure 2**) in contrast to the steady SST result.

The effect of the choice of discretization scheme (second-order upwind versus NVD-based Bounded Central Difference) for the convection term is shown in **Figure 4**, for both the DES and DHRL models. While the BCD scheme has been shown to significantly reduce numerical dissipation on structured meshes [35, 36], the effect on the mean flowfield using hybrid RANS-LES models was found to be insignificant. For the DES model, as the numerical dissipation increases, the local shear stresses decrease in magnitude, resulting in lower overall magnitudes of the eddy viscosity. This decrease in model eddy viscosity helps to compensate for the higher numerical viscosity of the second-order upwind scheme. For the DHRL model, a similar effect occurs. As seen in Eq. (10), the blending coefficient *α* responds to changes in the resolved flowfield to adjust the RANS contribution accordingly.

The underlying method by which the DHRL model operates is highlighted in **Figure 5**, showing the distribution of the blending parameter *α* through the channel for all three grids investigated. Recalling that *α* = 0 corresponds to a pure RANS mode and *α* = 1 corresponds to a pure LES mode, it is apparent that both the Medium and Fine grid cases have a significant portion of the channel, up to about *y*+ = 1000, for which the model operates as pure MILES. The MILES region extends further toward the wall for the Fine grid case (*y*+ ~ 100) than for the Medium grid case (*y*+ ~ 300). In the near-wall region, *α* → 0 as *y*→0 and the RANS model is recovered. Likewise, near the channel centerline (*y*+ > 1000), *α* < 1 and some portion of the RANS model stress is introduced into the momentum equations. The behavior of the Coarse grid case is quite different. Since the mesh is too coarse to resolve a significant level of turbulent production, the value of *α* is quite low throughout the channel, and a significant portion of the RANS model stress is introduced. As seen in **Figure 1 (a)**, the effect is to minimize overshoot of the mean velocity profile in the middle portion of the channel and recover a result close to the steady RANS SST result.

**Figures 6** and **7** illustrate some of the features of the resolved velocity field in the hybrid RANS-LES simulations. **Figure 6** shows isosurfaces of Q-criterion on the Fine grid, for both the DES and DHRL models, to highlight the presence of vortical turbulent structures. Such structures are visible for both models; however, for an equivalent value of Q-criterion, the structures predicted by the DHRL model are more prevalent and of a smaller scale than those predicted by the DES model. This was found to be consistently true for all cases. The reason is that the DHRL model operates in a MILES mode with respect to the fluctuating velocity field and as such the dissipation of small-scale turbulent features is less than with the DES model. **Figure 7** shows qualitatively the effect of mesh refinement on the resolution of the fluctuating velocity field using the DHRL model by showing contours of the spanwise vorticity on the periodic boundary planes. It is apparent that for the Fine grid (**Figure 7 (b)**), a significant portion of the turbulent spectrum is resolved, similar to the Q-criterion contours in **Figure 6 (b)**. However, for the Coarse grid (**Figure 7 (a)**), only large-scale, low-frequency fluctuations are present in the flow. It is this lack of resolution of the turbulence field that results in the decrease of the blending parameter *α* (**Figure 5**) in the Coarse grid case, and the introduction of a relatively large RANS stress contribution into the momentum equations.

## 5. Separated turbulent flow test case

The separated flow test case is the backward facing step flow matching the experimental study of Driver and Seegmiller [26], which is a commonly used benchmark case for turbulence model verification and validation. Two different three-dimensional meshes were created using the commercial meshing tool Ansys GAMBIT^{®} in order to evaluate grid sensitivity issues common to most current HRL models [15], for example, the delay of instability and breakdown of separated shear layers. For both grids, the computational domain size extended from 4 step heights (H) upstream of the step to 32H downstream of the step, 16H vertically from the wall downstream of the step, and 6H total in the spanwise direction. A structured multiblock mesh was generated, with a baseline Coarse mesh containing 744,960 total cells, and a Fine mesh containing 7,946,400 total cells. An average *y*+ value less than one was enforced to satisfy the recommended *y*+ requirements for the RANS turbulence model and for accurate HRL modeling. Nearly isotropic quadrilateral cells were used in the LES region (from x/H = 0.0 to 10.0) for maximum accuracy. A 2D planar slice of the Fine mesh is shown for illustrative purposes in **Figure 8**.

Simulations were performed using the DHRL, DDES, and SST k-ω models, with all simulation parameters and boundary conditions matching the experimental study. Profiles of inlet flow variables including mean inflow velocity, turbulent kinetic energy, and specific dissipation rate for the DHRL and RANS model simulations, and modified eddy viscosity for the DDES model simulation were created to match the distributions in the experimental data. For the spanwise boundaries, a periodic boundary condition was applied in the transient HRL simulations, and a symmetry boundary condition was used for the steady RANS simulation.

**Figure 9** shows the spanwise-averaged mean wall static pressure distribution for both the Coarse mesh and the Fine mesh, respectively, using three different turbulence models. In the recirculation region (x/H > 0), the Coarse mesh results shown in **Figure 9 (a)** indicate that both the DHRL and RANS computations predict a smooth pressure decrease and resolve the negative peak pressure with reasonable accuracy compared with the experimental data. The DDES results in contrast show overprediction of the pressure decrease and an offset in the location of negative peak pressure. The pressure rise predicted by both the DHRL and RANS models in the separated flow region shows similar behavior to the experimental results, but the wall pressure in the DDES computation shows a delayed pressure recovery. The pressure distribution after flow reattachment is similar to the experimental results for all simulations, except for a small overprediction by the DDES model in the region 10 < x/H < 16. **Figure 9 (b)** shows that the Fine mesh computations compare reasonably well with the experimental data for all models used. Simulations using the Fine mesh show a substantial improvement relative to the Coarse mesh for the DDES model. The differences underscore the mesh sensitivity inherent in DDES. In contrast, the Coarse and Fine mesh results using the DHRL model are similar, indicating that the DHRL model is relatively insensitive to mesh resolution, which agrees with the results from the attached flow test case above.

**Figure 10** shows the spanwise-averaged mean wall skin friction distribution on both the Coarse mesh and Fine mesh. **Figure 10 (a)** indicates that the wall skin friction result on the Coarse grid in the separated flow region agrees closely with the experimental data for both the DHRL model and the RANS model. In particular, the negative peak *C*_{f} values are quite close to the experimental result. In contrast, the DDES model shows an overprediction of the size of the flow separation region along with an overprediction of skin friction. The reattachment point was determined in the experiments using a linear interpolation of the oil-flow laser skin-friction measurements and was found to be x/H = 6.38 [26]. The DHRL model results show an early prediction of flow reattachment and a reattachment point of approximately x/H = 5.60, an underprediction of 9%. The RANS model predicts the flow reattachment location at x/H = 6.30, which is closer to the experimental measurements. The DDES results show a delayed flow reattachment at x/H = 9.23. The grid is too coarse for the DDES model to accurately resolve the Reynolds stress contribution in either LES or RANS mode, leading to the delay in flow reattachment. The skin friction is predicted reasonably well by the DHRL model downstream of the reattachment point. The RANS model accurately predicts skin friction immediately downstream of the reattachment point but shows an underprediction in the region farther downstream. The differences in the DHRL model prediction, as compared to the SST model, can be attributed to the presence of resolved unsteady turbulent fluctuations in the flow that lead to more rapid mixing of the separated shear layer and transport of momentum in the wall-normal direction. This is also clearly seen in the mean velocity and turbulent kinetic energy profiles shown below.

The predicted skin-friction distribution for the Fine mesh, shown in **Figure 10 (b)**, indicates that results with all three models agree reasonably well with the experimental data. However, with regard to mesh sensitivity, some models show quantitative or even qualitative differences between the Coarse and Fine meshes. The DHRL model predictions are closest to being mesh insensitive, as the *C*_{f} calculations for both meshes are quite similar. The reattachment location (x/H = 5.70) calculated using the Fine mesh is less than 2% different than the result on the Coarse mesh. The RANS model predictions for both the Coarse and Fine meshes show very close agreement with one another in all regions, including the separation point. This is similar to the attached flow results above, for which the steady RANS results were effectively mesh independent. In contrast, a significant level of mesh sensitivity is shown by the DDES model predictions. Skin friction over most of the downstream wall shows significant differences between the Coarse and Fine grid results. Mesh refinement significantly improves the prediction of reattachment point (x/H = 6.31) for the DDES model.

Previous studies have shown that DDES model predictions have reported delayed shear layer breakup and limited resolution of the Reynolds stresses in separated flow regions, particularly on coarse grids [2, 37]. One key advantage of the DHRL model is that it inherently addresses these potential problems. **Figure 11** shows contours of instantaneous spanwise vorticity (*z*-vorticity) obtained from the DHRL and DDES model simulations, respectively, on the Coarse mesh. The DHRL model results indicate a better ability to resolve turbulent eddies in the separated shear layer and a more obvious breakup and transition to LES mode than the DDES model. The relatively poor predictions of wall static pressure and skin friction predictions by the DDES model shown above are partly explained by these results. **Figure 12** shows the contours of instantaneous spanwise vorticity (*z*-vorticity) on the Fine mesh for both models. DHRL and DDES are both able to resolve turbulent scales in the separation region and as a result are able to provide reasonably accurate prediction of *C*_{p} and *C*_{f}. With regard to mesh sensitivity, the DHRL model contours are similar on both the Coarse and Fine mesh computations except for the increased resolution of small scales on the Fine mesh. The DDES model contours, in contrast, show even qualitative differences between results on the Coarse and Fine meshes. **Figure 13** shows the spanwise vorticity (*z*-vorticity) contours predicted by the RANS simulations on the Coarse and Fine meshes, respectively. It is apparent that no turbulent scales are directly resolved, as expected. Rather, the effect of all turbulent scales on the mean flow is included via the eddy-viscosity term.

In order to minimize the dissipation errors due to the numerical method, low dissipation schemes are generally used for convective term discretization in LES and HRL simulations. The bounded central differencing (BCD) [35] scheme is an example of a low-dissipation scheme developed to provide stable results on both structured and unstructured grids. Fine mesh simulations using the BCD scheme were performed along with baseline second-order upwind simulations to investigate the effect of discretization on the HRL model results. **Figure 14** shows the distribution of pressure in the streamwise direction, obtained using the Fine mesh. Very little difference due to discretization scheme is observed for the DHRL results. The DDES simulations, in contrast, show more sensitivity to the scheme, with the BCD discretization method producing a slightly overpredicted negative *C*_{p} peak and a delayed pressure recovery relative to the second-order upwind results.

**Figure 15** shows skin friction predictions on the Fine mesh for DHRL and DDES using both discretization schemes. Again, there is little difference in DHRL simulations due to the numerical scheme, though slightly better agreement with experiments is observed in the BCD scheme predictions just downstream of the reattachment location. The reattachment point (x/H = 5.75) for the DHRL with BCD scheme is close to the predicted reattachment with the upwind scheme (x/H = 5.70). For the DDES model, the BCD scheme results show a streamwise offset for the *Cf* prediction in the separated flow region compared to the upwind scheme. The BCD scheme results also show a significantly delayed reattachment point (x/H = 6.80) relative to the upwind scheme (x/H = 6.31). To summarize, the effect of discretization scheme in the DHRL simulations appears to be minimal, with a slight improvement in the results when using the BCD scheme. The use of the BCD scheme for the DDES simulations appears to result in slightly less agreement with experimental data, which is surprising. The apparent reason for this behavior may be the violation of the convection boundedness criterion of the BCD method due to the very low sub-grid scale turbulent diffusivity produced by the DDES model simulations.

**Figure 16** shows mean streamwise velocity profiles at several streamwise measurement stations for simulations using the Fine mesh. In the separated flow region (stations x/H = 1.0, 1.5, 2.0, and 3.0), the profiles clearly show that characteristics including boundary layer growth and size of the separation bubble in the wall normal direction are well resolved by both the DHRL and RANS simulations in comparison to the experimental measurements. The DDES simulation shows a small underprediction of the negative velocity peak in the near wall region and an overprediction of the size of the separation bubble. Simulations using each of the models show good agreement with the experimental data far from the wall. Immediately upstream of the reattachment point (station x/H = 5.0), DHRL and RANS results show reasonably good agreement with the data in terms of the negative velocity peak and separation-bubble size, but the DDES model somewhat overpredicts the flow behavior. It is apparent that in the separated flow region near the wall, the DHRL simulation shows nearly RANS-like behavior, which is expected since the DHRL model tends to limit to the RANS result close to the wall unless the grid is highly refined. It is interesting to note that the experimental velocity profile at station x/H = 6.0 clearly shows that the flow has reattached, which is early relative to the experimentally determined reattachment location of x/H = 6.38 based on the oil-flow laser skin-friction measurements [26]. The reason for this apparent discrepancy in the experimental data is not clear. The DHRL results also predict reattached flow at station x/H =6.0, while the RANS and DDES simulations show separated flow. Farther downstream (x/H = 7.0), all models clearly show not only reattached flow but also an underprediction of the velocity very close to the wall. In the wake recovery region (x/H = 16.0 and 20.0), the two HRL models show good agreement with the experimental data while the RANS model results show an underprediction in the rate of wake recovery, presumably due to lower predicted levels of turbulent mixing.

**Figure 17** shows computed turbulent kinetic energy profiles (normalized by square of mean inlet velocity) obtained for the Fine mesh simulations, compared with experimental data, at the measurement stations identical to those used to compare the mean velocity profiles. In the separated shear layer just downstream of the step at x/H = 1.0, DHRL results show the sudden rise in turbulent kinetic energy seen in the experimental data, indicating relatively rapid breakup of the shear layer characteristic of DHRL simulations. In contrast, the DDES model results show a relatively small peak in turbulent kinetic energy, which is characteristic of the more delayed shear layer break up characteristic of that model. This is also reflected in the delayed mixing of mean momentum shown at the same location in **Figure 16**. RANS model results show a relatively rapid increase in modeled turbulent kinetic energy downstream of flow separation. Both the DHRL and RANS results show a steady increase in turbulent kinetic energy similar to the experimental data up to station x/H = 5.0. The DDES simulation shows a relatively slow increase at stations x/H = 1.5 and 2.0, but a more rapid increase at x/H = 5.0, characteristic of a delayed but sudden breakup of the shear layer. Near the reattachment location at stations x/H = 6.0 and 7.0, the computed turbulent kinetic energy profiles decay similar to the experimental results, but the peak value of turbulent kinetic energy is underpredicted by the DHRL and RANS models results and significantly overpredicted by the DDES results. Farther downstream (stations x/H = 16.0 and 20.0), results obtained with all three models show qualitative agreement with the experimental data.

## 6. Summary and conclusions

The recently proposed dynamic hybrid RANS/LES modeling framework (DHRL) has been investigated for two test cases, turbulent channel flow and flow over a backward facing step, in order to investigate the capability of the model to accurately predict both attached and separated turbulent flows. The test cases elucidated the behavior of the DHRL model in comparison with a more traditional hybrid RANS-LES model, including model response to mesh refinement level and changes in discretization scheme.

For the channel flow test case, the DES and DHRL models showed significant resolved velocity fluctuations for all mesh refinement levels and for both discretization schemes, on the order of the turbulence levels in the reference DNS data. In contrast, the DDES model yielded a steady RANS-type solution, which is consistent with its intended performance. The DES model exhibited the well-known log-layer mismatch in the mean velocity, which was significantly reduced for the DHRL model. Furthermore, the DES model showed significant mesh dependence, with results improving as the mesh was refined, as expected. The DHRL model showed almost no mesh sensitivity with regard to the mean flow results, although the resolved turbulent kinetic energy levels, as well as the range of resolved scales of fluctuating motion, increased with increasing mesh refinement. The DDES model and SST models showed no mesh dependence due to the fact that each yielded a steady-state result. Both the DES and DHRL models showed almost no sensitivity to discretization scheme for the attached flow case.

For the backward facing step test case, the DHRL results showed rapid breakdown of the shear layer from RANS to LES type flow even for the Coarse mesh simulations, and relatively small differences in results on both the Coarse and Fine meshes, as well as for two different discretization schemes. The DDES model results showed a delayed breakdown of the shear layer and substantial mesh dependence, as well as nontrivial dependence of the results on the discretization scheme. The mean flow results predicted using the DHRL model were similar to the SST model results in the separated flow region, but showed improved prediction downstream of the reattachment point. The DHRL results were also comparable to SST results and superior to DDES results in terms of predicted turbulent kinetic energy downstream of the step, although it should be noted that for the SST model turbulence was completely modeled while for the DHRL model turbulent kinetic energy was resolved as unsteady velocity fluctuations. A key benefit to the DHRL model is the potential to effectively address grid dependence issues present for most currently used HRL models in predictions of the mean flow, while providing for increasing resolution of unsteady turbulence fluctuations while operating in LES mode as the grid size is reduced.

Overall, the results indicate that the DHRL model performs as designed. The attached flow case demonstrates a smooth variation between the RANS and LES zones, perpendicular to the flow direction, controlled by the model blending parameter *α*. Likewise, the separated flow case showed a similar behavior, with transition between RANS and LES regions occurring in the streamwise and wall-normal directions. The lack of explicit grid dependence in the DHRL model terms resulted in minimal grid dependence of the mean flow results, in contrast with the DES and DDES models. Furthermore, the DHRL model was shown to agree closely with the RANS component (SST) results for the attached flow case, while resolving significant levels of turbulent fluctuations in the boundary layer, at least on well-refined grids. Results indicate that the DHRL model may provide a more flexible and universal method for hybrid RANS-LES than existing standalone models.