Assessing IDDES-Based Wall-Modeled Large-Eddy Simulation (WMLES) for Separated Flows with Heat Transfer

: The Wall-modeled Large-eddy Simulation (WMLES) methods are commonly accompanied with an underprediction of the skin friction and a deviation of the velocity proﬁle. The widely-used Improved Delayed Detached Eddy Simulation (IDDES) method is suggested to improve the prediction of the mean skin friction when it acts as WMLES, as claimed by the original authors. However, the model tested only on ﬂow conﬁgurations with no heat transfer. This study takes a systematic approach to assess the performance of the IDDES model for separated ﬂows with heat transfer. Separated ﬂows on an isothermal wall and walls with mild and intense heat ﬂuxes are considered. For the case of the wall with heat ﬂux, the skin friction and Stanton number are underpredicted by the IDDES model however, the underprediction is less signiﬁcant for the isothermal wall case. The simulations of the cases with intense wall heat transfer reveal an interesting dependence on the heat ﬂux level supplied; as the heat ﬂux increases, the IDDES model declines to predict the accurate skin friction.


Introduction
Wall-bounded turbulent flows play a key role in industrial applications in addition to their importance in relatively simple configurations for academic research. The main feature of wall-bounded turbulent flows is the existence of the turbulent boundary layer with the multi-scale nature throughout, which becomes gradually distinct with increasing Reynolds number [1,2]. Consequently, with increasing Reynolds number, several constraints on the grid resolution are imposed on the simulation domain if using Large-Eddy Simulation (LES), in which all of the large energetic scales throughout the boundary layer need to be resolved. This could lead to serious challenges of modeling when flow separation, drag, or heat transfer are of interest [3]. In recent years, wall-modeled LES has gained momentum as a high fidelity tool for industrial designs [4]. In WMLES, the large-scale unsteady, energetic turbulent configurations in the outer region of the boundary layer are resolved, meaning structures with low energy and universal behavior to be modeled. This enables a competitive computational cost compared with other high fidelity approaches. These WMLES approaches generally fall into one of two categories: (i) Methods that model the wall shear stress directly, and (ii) methods that switch to a Reynolds-Averaged Navier-Stokes (RANS) model in the inner layer [5]. The second category includes hybrid LES/RANS and Detached Eddy Simulation (DES). The hybrid LES/RANS and DES methods are popular in both academic and industrial applications. However, hybrid LES/RANS and DES accompanied with a deviation of the velocity profile named "Log-Layer Mismatch" (LLM), resulting in the inaccurate prediction of skin friction that is on the order of 15%, which is too low [6]. It was reported that in the switching location, the amount of turbulent viscosity contribution to the shear stress is too small, while energetic eddies have not yet been developed [7]. This could lead to challenges of modeling when flow separation, drag, or heat transfer are of interest [3]. This problem was first reported by Nikitin et al. [6] who used the standard DES as a wall model in which the switching point from RANS to LES is located in the logarithmic region.
The Improved Delayed DES (IDDES), which is a widely-used DES-like method pioneered by Shur et al. [8], to improve the deficiency of the hybrid RANS/LES model. In subsonic wall-bounded turbulent flows, the good performance of IDDES acting as WMLES has been validated (e.g., see Shur et al. [8,9]). However, some studies [10][11][12][13] showed that the IDDES model does not thoroughly resolve the deficiency. When using IDDES as a WMLES in supersonic boundary layer flows, Peterson and Candler [10] reported that the IDDES model did not resolve LLM. Gieseking et al. [11] reported a severe underprediction of skin friction at about 25%. Peterson and Candler [10] concluded that the reason attributed to the fact that the IDDES model was calibrated using low-speed canonical flows. However, Mockett et al. [10] showed that LLM could occur even in IDDES of a low-speed channel flow. The deficiency of the IDDES model in low-speed channel flows was also reported by Saini et al. [13]. Mockett et al. [10] argued that the problem might be due to the numerical sensitivity of the IDDES model. Guseva et al. [14] assessed the IDDES model combined with the shear-layer adapted subgrid model in separated flows and observed that it considerably mitigates the delay of transition from fully-modeled to partially-resolved turbulence in the separated shear layers compared to the standard DES.
As mentioned earlier, LLM could be a critical issue when flow separation, drag, or heat transfer are of interest. Xu et al. [15] reported WMLES results of low-speed turbulent flows in a plane channel with heat transfer and compared the results with Direct Numerical Simulations (DNSs) of low-speed plane channel flow. Jakirlic et al. [16] assessed a zonal hybrid LES/RANS method in backward-facing steps with heat transfer. To the authors' knowledge, the evaluation of WMLES performance on flows with heat transfer is limited to a few investigations only. Thereby, it is of interest to develop an assessment of the IDDES model acting as WMLES in the simulation of separated flows with heat transfer. This paper aims to systematically assess the IDDES model in terms of its ability to predict heat transfer. We limit our discussion to low-speed flows, i.e., flows with Mach numbers well below one. For these flows, the momentum equation and the energy equation are only weakly coupled. Two configurations of backward facing steps with an isothermal wall and walls with heat fluxes are considered, and the results are compared to previously-published LES and experimental data.

IDDES-Based WMLES
A hybrid RANS/LES model solves unsteady 3D Navier-Stokes equations on an embedded near-wall mesh with a RANS-type closure. Recently, significant progress has been reached in building-hybrid RANS/LES models for high Reynolds number massively separated flows and more general cases that include attached and mildly separated regions [17]. For the former, one of the widely-used approaches is DES and its modification, Delayed DES (DDES) [18]. A general feature of these two methods is that the major part of the attached boundary layer is treated by RANS, while LES is applied only in the separated flow regions. In contrast, in the second group of RANS/LES hybrids, which offers wallmodeling in LES of high Reynolds number flows, RANS is used only in a much thinner near-wall region, in which the wall distance is much smaller than the boundary layer thickness but is still potentially very large in the wall unit. Shur et al. [8] pioneered the IDDES model, which builds a single set of formulas both for natural DDES and WMLES applications so that different regions inside a single simulation over a complex geometry can each be treated by a very effective model. They concluded that IDDES significantly extends the application area of DDES by allowing the activation of RANS and LES in different flow regions, giving an admittedly intricate but well-balanced and strong numerical approach to complicated turbulent flows at high Reynolds numbers. The concept of the IDDES model couples RANS and LES approaches via the introduction of the following blended RANS-LES length-scale: where l RANS and l LES are the RANS and LES length scales, respectively. The non-zero function f e corrects the RANS behavior, and empirical blending-function f B varies from 0 to 1 to provide rapid switching of the model from RANS mode within the wall distance.
More details about this model can be found in Shur et al. [8].

Governing Equations
The following filtered compressible Navier-Stokes equations are solved that is designed to resolve the outer layer scales in WMLES [19]: where τ ij =ρ u i u j −ρũ iũj , the subgrid-scale stress tensor, and q j =ρC p u j T −ρC pũjT , the subgrid-scale heat flux, representing the effect of the small scales of turbulence on the large ones. The set of equations is closed by setting ∂ũ k ∂x k δ ij and by the modified averaged Equation of State, P =ρRT. Note that µ t and Pr t are the turbulent eddy viscosity and turbulent Prandtl number, respectively.
Here the SST-based IDDES model is used to model the eddy viscosity. The governing equations of the SST IDDES model are based on solving transport equations for a specific turbulent dissipation rate (ω) and turbulent kinetic energy, which read as follows [20]: where l IDDES is presented in Equation (1). Therefore, the eddy viscosity can be calculated as: The details of coefficients are used in Equations (5)-(7) and the IDDES length scales can be found in Reference [20].
A WMLES that solves the fully-coupled RANS equations requires an auxiliary walllayer grid to be provided ( Figure 1). The inner layer RANS is solved in a grid that is refined in the wall-normal direction (in RANS context) and is embedded in the outer layer (wall-modeled) coarse LES mesh. At each time step, the wall-shear stress τ w and wall-heat flux q w calculated from the wall-model solution are used as Neumann boundary conditions for the LES. The solutions from the LES are transferred to the inner-layer RANS model at the corresponding point. The matching point in the LES mesh where the RANS top boundary meets the LES mesh is not necessarily at the first off-wall LES nodes where subgrid modeling and numerical errors are large because of the poor near-wall grid resolution in the LES mesh [3].

Errors in Wall Modeling
Next, we discuss the numerical and subgrid-modeling errors in the first few grids off the wall. Since the wall model takes as input the solution from the LES at the first grid point, it means that the wall model is fed by inaccurate information. Therefore even a perfect wall model would not be able to accurately predict skin friction. The second error is the well-known LLM. The inner log-layer, provided by the RANS model and the outer log layer computed by LES are mismatched. It is reported that the total shear stress (µ + µ t ) ∂Ũ ∂y and heat flux −C p ( µ t Pr t + µ Pr ) ∂T ∂x j between the LES and RANS has discrepancies at the matching location.

Flow Solver
An in-house finite-volume code is used to carry out WMLES. The code can capture discontinuities such as shockwaves while preserving the energy in turbulent regions [21]. Therefore, it is appropriate for the LES type of turbulent flow simulations. A skew-symmetric form of a collocated second order scheme with low numerical dispersion and dissipation is implemented to model the turbulent region in the combination of a semi-discrete, central difference scheme to capture discontinuities, with a switch based on the Ducros sensor [22]. This sensor allows the hybrid scheme to capture the discontinuity and predict the accurate decay of turbulent kinetic energy in turbulent regions, making the model suitable for the LES of compressible turbulent flows. The governing equations are integrated in time using the original implicit Euler scheme. The code was validated for numerous cases of compressible and incompressible flows [23][24][25][26]. The IDDES model is implemented in the solver to study the boundary layer separation and reattachment in transonic airfoils [27]. In addition, the IDDES-based WMLES approach is used to assess the WMLES performance in hypersonic flows with heat transfer [28].

Flow Configuration
Two different flow configurations over backward-facing steps are considered. Figure 2a shows the first configuration in which the turbulent boundary layer at the edge of the step separates (fixed separation point) and reattaches at the bottom wall forming the recirculation region behind the step. The flow Reynolds number based on the step height (H = 0.038 m) is Re H = 28,000. The upstream flow conditions correspond to a developing channel flow of height 4H. The boundary layer developing at the lower wall separates at the step edge transforming into the free shear layer within the expanded channel section of height 5H. The boundary layer thickness at the step edge is approximately equal to δ = 1.0658 H. This configuration is similar to the experimental database is provided by Vogel and Eaton et al. [29]. The lower wall in the channel after the expansion is heated by the heat flux q w = 270 W/m 2 . In the spanwise direction z + = 12 and the wall normal direction y + = 1.5. The streamwise spacing is set to reach 10 ≤ x + ≤ 14.7. This domain was meshed with 312 × 120 × 84 grid cells, resulting in about 3 million cells in total. For sufficient confidence in the computational results, the number of the girds is determined by a mesh convergence study. Figure 2b shows the second configuration. The flow Reynolds number based on the step height (H = 0.041 m) is Re H = 5540. The upstream conditions correspond to fully-developed flow in a channel of height 2H (inflow was generated by performing precursor channel flow calculations). An isothermal wall and two strong wall heat fluxes (q w = 1 kW/m 2 , 3 kW/m 2 ) are considered for the second configuration. All other walls are treated as adiabatic. Digital filter-based inlet boundary condition is used to generate synthetic turbulence-alike time series [30]. A non-reflecting boundary condition is imposed to attempts to reconstruct this kind of non-reflective scheme allows the sound wave flow smoothly out of the domain through the outlet. Two LES simulations serving as a computational reference in the present work have been performed by Akselvoll and Moin [31] and Keating et al. [32]. In the spanwise direction z + = 12 and in the wall normal direction y + = 1.5. The streamwise spacing is set to reach 10 ≤ x + ≤ 14.7.

Results and Discussion
Some selected results obtained by the IDDES model for both configurations are shown and discussed in this section. The reference experimental database (Kasagi and Matsunaga [33]; Vogel and Eaton [29]) are used for comparisons.

First Configuration
The results of the reattachment length obtained by several computational methods and experiment are summarized in Table 1. The IDDES model slightly overpredicts the mean reattachment point compared to the experimental and LES results. The iso-surface of Q colored with the instantaneous x-velocity is shown in Figure 3. The turbulence coherent Q-criterion, was originally introduced by Jeong and Hussain [34] to visualize a vortex, and it is defined as Q = 1 2 (||Ω|| 2 −||S|| 2 ) where Ω and S are rotation and strain rate tensors, respectively. The coherent structures in the undisturbed boundary layer present classic streamwise elongated hairpin vortices. The Kelvin-Helmholtz instability at the edge of the step breaks the upcoming longitudinal structures. Near the reattachment line, some large-scale coherent structures in the boundary layer can be identified. In the recovery region, detached structures from the mixing layer are observed.

Experiment [29]
6.70H  Figure 4 shows the evolution of the mean axial velocity profiles at selected streamwise locations. After the step, in the recirculation zone, the inner velocity is reduced with the thickening of the low-speed inner layer. At x/H = 3.0, 4.7, and 5.93, the negative value of velocity can be observed, indicating the appearance of the reverse flow. After the reattachment, the boundary layer redevelops and finally recovers close to an equilibrium state near the exit of the computational domain. It can be seen that IDDES profiles go under faster development compared to LES results however, they exhibit good agreement with the experimental results. Figure 5 depicts the development of the profiles of the mean temperature at several streamwise locations. Some deviations of the computed temperature field from the reference experiment are noticeable only immediately behind the step in the separation zone (locations x/H = 0.33 and 1.67).
The evolution of the spatial distribution of the friction coefficient, Figure 6. The skin friction exhibits a clear plateau region with negative values, which indicates the existence of the separation bubble with reverse flows. The skin friction decreases with increasing the adverse pressure gradient and reaches the minimum value. After reaching the minimum value, C f begins to increase, and the reattachment point can also be identified by C f = 0. Further downstream, the increase of C f slows down, which indicates the recovery of the turbulent boundary layer towards another equilibrium state. The IDDES predicted C f follows the LES and experimental results closely. However, the underprediction of C f is apparent for the wall with heat flux, especially in the separation region. The underprediction is more noticeable compared to LES-predicated C f . This underprediction which is reported in many cases of using the IDDES model is due to the the log-layer mismatch at the switching location. The St values show good agreements with LES and experimental results with an underprediction both in the recirculating and recovery regions. The reason can be attributed to the mean temperature profiles are shown in Figure 5. Since IDDES goes under a faster temperature recovery, the value of the IDDES-predicted St number is lower in this region.

Second Configuration
The second configuration represents the isothermal wall and also the strong heat flux of the wall, which can affect the fluid properties. The comparison of the mean streamline patterns in the recirculation region for the isothermal case, Figure 7, shows that the IDDES model yields the reattachment length of x R /H = 6.85, which is slightly larger compared to the experimental one x R /H = 6.51. A discrepancy also is apparent in the corner bubble. The core of the recilculation region moved forward for the IDDES model compared to the analogous experiment, and it resides near the reattachment point. The discrepancy is seen between the IDDES and the experiment in the scale and location of the vortex structures for the recirculation zone can be attributed to the predicted properties of the boundary layer in the vicinity of the reattachment location. The C f evolution obtained by the IDDES model agrees reasonably with experimental results with respect to the reattachment length as shown in Figure 8. In contrast, the LES result reports the reattachment length to be slightly shorter. Furthermore, the negative peak value does not suffer from the significant underprediction for the isothermal wall case. The instantaneous and mean velocity for the isothermal wall are shown in Figure 9. Figure 10 shows the mean velocity and Reynolds shear stress in various streamwise locations. It can be seen that the results from the IDDES model and experimental data have very good agreement in all distinct regions of the backward-facing step with an isothermal wall.  It was discussed that the IDDES model suffers from the underprediction of the skin friction for the case of the wall with a mild heat transfer rate (first configuration). The simulations of the cases with an intense wall heat transfer are performed for two different wall heat fluxes, q w = 1 and 3 kW/m 2 . The C f evolution in the streamwise direction displayed in Figure 11a reveals a very interesting dependence on the heat flux level supplied. Directly comparing the mean C f and St (Figure 11b) with the experimental data reveals that as the heat flux increases, the IDDES model declines to predict the accurate C f . The St prediction also follows the same trend. However, it is apparent that the strong temperature variation does not significantly influence the reattachment length. Figure 12 shows that the trend of the wall temperature variations agree very well with the experimental reference data. The position of the peaks of the temperature is well predicted by the IDDES model. It can be seen that the temperature is overpredicted by the computational results in the separation zone for both cases; however, it suffers from an underprediction in the recovery region. The overall prediction is quite acceptable for the case of the wall with the lower heat transfer rate (q w = 1 kW/m 2 ).

Conclusions
This study is conducted to assess the performance of the IDDES model acting as a WMLES in separated flows with heat transfer. Two configurations of backward-facing steps with an isothermal wall and walls with heat transfer are considered, and the results are compared to previously-published LES and experimental data. Qualitatively, the results obtained by the IDDES model with respect to the reattachment lengths, C f and Stanton number distributions, fluid flow, and thermal fields follow the reference experiment and reference LES results closely. For the case of the wall with heat flux, the skin friction and Stanton number are underpredicted by the IDDES model however, the underprediction is less significant for the isothermal wall case. The simulations of the cases with intense wall heat transfers reveal an interesting dependence on the heat flux level supplied, as the heat flux increases, the IDDES model declines to predict the accurate C f . A possibility for future work is to examine the IDDES model for a broader range of Reynolds numbers, including high speed flows, to assess the performance of the IDDES model in high Reynolds number flows with various heat transfer configurations. A possibility for future work is improving the error of the wall modeling by investigating the blending function between RANS and LES to find the optimal switching location for the IDDES model.