Generalized wall-modeled large eddy simulation model for industrial applications Open Access

Turbulent flows are present in multiple industrial applications such as aircraft and vehicle aerodynamics, turbomachinery, medical devices, and weather prediction. The understanding of turbulent flow is important, since turbulence produces effects that include the increase in friction losses and heat transfer, the delay of separation in flows exposed to adverse pressure gradients, and the noise generation. The complete flow description is obtained with the direct numerical simulations (DNS), in which the entire turbulence spectrum is resolved. Estimates derived from Kolmogorov's hypotheses show that, for a turbulent flow above a flat plate of length L, the number of grid points is $N DNS ∼ R e L 37 / 14$⁠.¹ 

The computational cost of DNS is reduced with the large eddy simulation (LES) models, which solve part of the turbulence spectrum and model the remainder. The dynamics of the large-scale motions, which are affected by the flow geometry and are not universal, are computed explicitly. The small scales have a universal character and follow the Kolmogorov hypotheses, and they are represented by simple models. As a general rule, LES simulations require a minimum of 80% of the turbulent kinetic energy to be resolved,² which is restrictive in wall-bounded flows at high Reynolds numbers; for example, in the turbulent flow above a flat plate, the number of grid points is $N LES ∼ R e L 13 / 7$⁠.¹ The scaling problem near the wall is reduced in the wall-modeled LES (WMLES) models, where a mixing-length model or a wall function is used up to the inner part of the logarithmic layer. For a WMLES simulation of a turbulent flow above a flat plate, the number of grid points is $N WMLES ∼ R e L$⁠.¹ The accuracy of current WMLES models is limited by the log-layer mismatch (LLM), which is the deviation of the velocity from the log-law, especially where the model switches from the wall function (or the mixing-length model) to LES.³ Multiple solutions have been proposed, including the addition of forcing terms in the momentum equations,⁴ although the constants of the forcing terms are not universal. This is partially solved in dynamic approaches that are based on control systems,⁵ although at the expense of the increase in the complexity and the numerical stability. The error of the LLM is proportional to the dimensionless wall distance $y ̃ +$ of the first near-wall cell, and mesh refinement is not always possible due to the requirement of using cells of low aspect ratio for resolved turbulence. If a mixing-length model is used near the wall,⁶ the use of high aspect ratio cells is limited due to the reduced width in which the Reynolds-averaged Navier–Stokes (RANS) part is applicable.

In the Reynolds-averaged Navier–Stokes (RANS) models, the computational cost does not depend on the Reynolds number. The entire energy spectrum is modeled, and since the variables of the equations are mean quantities, temporal resolution is not required in the simulations. RANS models have shown great accuracy in wall-bounded attached flows, where the law of the wall establishes a reference for the derivation of the model constants. Two-equation models are currently the reference for industrial simulations, since they provide the optimum balance between computational cost, stability, and accuracy. However, mean quantities are not well predicted in free shear flows⁷ and separated flows,⁸ which constitute an important part of industrial applications. Furthermore, due to their time-averaged formulation, RANS models are usually not suitable for predicting time-dependent quantities, for example, in acoustics simulations,⁹ fluid–structure interactions,¹⁰ and multiphase flows.¹¹ There are particular cases—namely, massively separated flows—where unsteady RANS simulations (URANS) provide unsteady solutions, and specifically part of the turbulence spectrum, although the accuracy is limited in these situations.¹² 

In wall-bounded flows at high Reynolds numbers, the computational cost of LES is unaffordable for practical engineering applications. On the contrary, the cost in RANS is several orders of magnitude lower, and the results given by these models for boundary layers without detachment are often acceptable. In free shear flows, much less computational power is required to perform LES simulations. A compromise solution in terms of accuracy and computational cost is given by the hybrid RANS/LES models, where RANS is used in the boundary layers, and LES in the rest of the domain. Both WMLES with the mixing-length approach and hybrid RANS/LES models share the idea of blending RANS and LES models, which leads to the use of both names interchangeably by some authors. In the present article, WMLES is assumed to be a model where the blending takes place at a $y ̃ +$ value in the inner part of the log-layer (i.e., $10 < y ̃ + < 50$⁠). On the contrary, hybrid RANS/LES models are considered those where transition takes place at higher values of $y ̃ +$ up to the edge of the boundary layer.

A reference hybrid RANS/LES model is the detached eddy simulation (DES) model. The two-equation DES model introduces the grid spacing $Δ = max ( Δ x , Δ y , Δ z )$ in the dissipation term of the turbulent kinetic energy equation,¹³ where Δ_i refers to the grid spacing in the i direction. The model switches to LES if the turbulence length scale $l t > C DES Δ$ for a constant $C DES ≈ 1$⁠; therefore, it reduces to a RANS model in coarse meshes, and to LES in fine meshes. However, for intermediate conditions, the transition from RANS and LES can occur without control. If the original mesh satisfies $l t < C DES Δ$⁠, then the model stays in RANS [see Fig. 1(a)]. If the mesh is refined so that $l t > C DES Δ$ in some regions of the boundary layer, then the model switches to LES. Since l_t reaches its maximum value in the middle of the boundary layer,¹⁴ a real situation can be the one shown in Fig. 1(b), where δ is the thickness of the boundary layer. This is a stable flow because there is no phenomenon that triggers an instability, so no vortex structures are generated, leading to early separation, especially in flows exposed to an adverse pressure gradient. This numerical behavior is called grid-induced separation (GIS).¹⁵ If the mesh is further refined, the switching of the model from RANS to LES formulations occurs closer to the wall; that is, the width of the LES zone is extended, while there is no certainty that the mesh is fine enough for a LES simulation. Thus, in DES, mesh refinement does not guarantee a more accurate prediction.

An improvement of DES called delayed-DES (DDES) was published later¹⁵ that introduces a shielding function in the formulation. While in DES the effect of GIS is noticeable from $Δ ≈ δ$⁠, in DDES, it is reduced to $Δ ≈ 0.25 δ$⁠.¹⁴ However, the impact is much more drastic than in DES, since for $Δ < 0.25 δ$⁠, small changes in the mesh size produce significant variations in the solution. As a consequence, DES/DDES can provide less accurate solutions than RANS.¹⁶ Another important aspect of DES/DDES is that the transition from RANS to LES takes place in a distance whose size is not negligible; in particular, in DDES, it has been shown that it is around $0.3 δ$ in a flat plate.¹⁷ This is called gray area in the literature, since the turbulence model is not well defined there. Gray areas delay the generation of turbulent structures and generate spurious buffer layers¹⁸ and the LLM.¹⁹ 

Variants of the DDES model have been proposed, based on the use of the blending functions of the $k − ω$ SST model,²⁰ which reduce the GIS down to $Δ ≈ 0.1 δ$⁠, but the gray areas are not eliminated. More elaborated formulations involve the use of additional complex functions that depend on flow variables that act as sensors.²¹ The model cited solves the problem of the gray areas, although the GIS is not addressed. Moreover, it includes multiple additional adjusting constants, and it has not been verified that the values proposed are universal. The improved-DDES (IDDES) model is a modification of DES developed with the capability to run in WMLES mode,⁶ but with limited shielding against GIS. For meshes that are not suitable for WMLES, the accuracy of DDES and IDDES is similar,²² and recent studies have shown that IDDES does not solve completely the LLM.²³ 

Another weakness shared by all DES family is the reduced accuracy of the LES part. By assuming an equilibrium of the source terms of the RANS turbulence variables, it can be demonstrated that the LES formulation corresponds to a Smagorinsky model; specifically, for the original constant $C DES = 0.61$ (Ref. 13) used in DDES with the $k − ω$ SST model as the RANS model,²⁴ the corresponding Smagorinsky constant outside the boundary layers is $C s ≈ 0.175$⁠. This value is too dissipative for free shear flows, for which $C s ≈ 0.1$ has been shown to be optimal.²⁵ DDES also inherits limitations from the Smagorinsky model, namely, the need to adjust the model constant based on the specific flow being simulated. An additional deficiency is the persistent presence of non-zero turbulent viscosity in standard laminar regions, including laminar shear layers, where it should ideally vanish.

Similar to the DES family, the scale-adaptive simulation (SAS) model is a formulation that takes RANS equations as reference, and a second equation for the turbulence scale is reformulated.²⁶ This results in the inclusion of the second velocity space derivative in the RANS equations, so that the scale equation is able to adjust to resolved scales in the flow. SAS has been shown to provide performances superior to that of RANS (URANS) models in flows with strong instabilities.^27,28 However, for a given mesh suitable for LES far from walls, when the flow instability is not enough to enable the LES model, SAS solves turbulent structures of a size similar to those of URANS,²⁹ and the model can be less accurate than URANS.³⁰ Thus, in these cases, gray areas are produced, since the model does not solve the equations neither in URANS nor in LES mode. While in the DES family the location, width, and conditions of the gray areas are reasonably known, in SAS that knowledge does not exist. Other hybrid RANS/LES models available in the literature solve the partially averaged Navier–Stokes (PANS) equations.³¹ The difference between RANS and PANS is that the latter contains modified parameters in the turbulence dissipation and diffusion terms, so that, under some restrictions, the required level of modeled-to-resolved turbulent kinetic energy and dissipation rate can be set. PANS does not include additional shielding functions, so the presence of gray areas, the LLM, and the GIS are potential problems. The exact transition point from RANS to LES can be defined with the Zonal RANS/LES methods, where RANS and LES are solved in different parts of the domain, and a synthetic turbulence generator transfers modeled RANS turbulence to resolved LES turbulence.³² However, there are not many industrial applications in which the location of the RANS/LES interface can be accurately defined.

The stress-blended eddy simulation (SBES) approach has been developed by Ansys to overcome the deficiencies of the original DES model and its variants.¹⁴ The basis is the blending of the RANS and LES models at the stress level, which is equivalent to the blending of the turbulent viscosities if both models use the Boussinesq turbulent viscosity assumption. All RANS and LES models can be combined, which is especially relevant in the LES case, since more consistent models than the Smagorinsky model can be used. In addition, motivated by the higher time scale of the RANS turbulence quantities in comparison with the rest flow variables, in SBES the equations of RANS variables can be solved every 5–10 flow time steps. This benefit is obtained due to the decoupling of the RANS and LES turbulent viscosities and leads to a computational cost similar to the WMLES/LES models, which is lower than the one of DES/DDES/IDDES that solve the equations of the RANS variables at all time steps. Ansys claims that the blending function is defined, such that the RANS model covers the entire boundary layer, with no GIS problem, and with gray areas of negligible size. It has been verified that the model switches to WMLES mode depending on the mesh size and the turbulence content of the flow. However, the details of the formulation are not published, and it is available only in Ansys Fluent/CFX. SBES has shown similar or greater accuracy than IDDES/DDES in a mixing layer,¹⁴ automotive applications,³³ and multiphase flows.³⁴ 

Hybrid RANS/LES models can also be applied to transitional flows with transition RANS models, namely, the $γ − R e θ t$⁠,³⁵ which has been validated in turbomachinery applications.³⁶ The model has shown to be capable of predicting some of the physics of the laminar transition in further studies on wing profiles³⁷ and wind turbines.^38,39 However, in some cases, discrepancies with experimental data have been encountered near the separation regions. Transitional flows have been studied with LES models, although significant errors have been observed in the early stages of separation bubbles.⁴⁰ This is because the residual stress models do not detect the entire laminar regions in medium and coarse meshes, so that the resolution needs to be increased to levels that are close to DNS simulations.⁴¹ Alternatives have been proposed using sensors to detect the transition point,⁴² although this methodology implies modeling and the addition of uncertainties in general cases. Since the LES approach combines the high CPU cost and the modeling requirement, an interesting approach is to use a transition RANS model as a basis for a hybrid RANS/LES model in these cases. However, as noted in Ref. 43, the number of publications on this topic is reduced.

In summary, in the industrial context, DNS simulations are not an option in the near future, and LES is limited to free shear flows due to the high computing power required to perform these simulations. Given its limited accuracy, RANS is an option just for preliminary simulations and/or non-complex turbulent flows. Thus, WMLES simulations for moderate Reynolds numbers and hybrid RANS/LES models for high Reynolds numbers are stated as reference for industrial flows in the near future.^44,45 In the hybrid RANS/LES case, provided there is no GIS, the requirement of solving the outer layer of the boundary layer in RANS can be dropped, since it is not universal and the mesh requirements are of the same order as in the freestream region. Thus, both WMLES and hybrid RANS/LES can be unified to a model that solves in RANS the viscous sublayer and/or the log-layer, depending on the simulation requirements.

The main goal of this work is to develop a generalized WMLES formulation (GWMLES) that extends the classical WMLES approach. The model is divided into two separated formulations: GWMLES_sub, where turbulence is modeled up to the inner part of the log-layer, and $GWMLES log$⁠, where the viscous sublayer and the entire log-layer are modeled. The GWMLES_sub aims to improve current WMLES models, while the purpose of the $GWMLES log$ is to reduce the computational cost of the DES family, as well as to eliminate their shortcomings that lead to solution deterioration (see Fig. 2). The GWMLES model is based on the SBES approach of blending the turbulent viscosities, but with two key differences based on the information provided by Ansys: first, that in $GWMLES log$ the RANS model does not cover the outer layer of the boundary layer, which is not universal, and second, that the selection between $GWMLES log$ and GWMLES_sub is done by the user and not depending on the mesh size. The remaining requirements for the blending function are shared with SBES to overcome the deficiencies of the DES family, namely, the dissipative LES part, GIS, and the gray areas. This paper is organized as follows. First, the methodology is presented in Sec. II, including the details of the model. Then, in Sec. III, the results given by the model are compared with DDES, IDDES, SBES, and URANS, as well as with experimental data and DNS/LES/WMLES results available in the literature; in particular, in four test cases that cover three key industrial applications: a free shear flow, two turbulent wall-bounded flows with both GWMLES submodels, and a laminar–turbulent transitional flow. Finally, the conclusions and future work are drawn in Sec. IV.

The model equations are listed in Sec. II A, which are closed with the expression of the blending function of the model that is developed in Sec. II A 1. The numerical implementation is discussed in Sec. II B, including, in Sec. II B 3, the definition of quality estimators used in numerical simulations for model validation.

The parameters of GWMLES derived from the previous analyses are given in Table I.

Ansys Fluent 2023R1 is used,⁵⁰ which incorporates Eq. (5) by default, and the $f GWMLES$ blending function defined in Eq. (7) is implemented with a user-defined function (UDF).

The poly-hexcore method of Ansys Fluent is used for meshing, which creates a mesh consisting of octrees in the bulk region, and keeps a layered hexahedral mesh in the boundary layers. These elements are connected with general polyhedral elements (see Fig. 5), which have been used successfully in LES simulations of jets⁵¹ and external aerodynamics.⁵² The computational mesh is generated with different cell sizes Δ ranging from $Δ min$ to $Δ max$⁠, and the finest cell size is placed near the most relevant walls of the computational domain (see Fig. 5).

The use of cells with an aspect ratio greater than 1 has been shown to deteriorate the accuracy of scale-resolving simulations.^53–55 Thus, the mesh is generated so that $f GWMLES = 1$ (RANS) in the hexahedral cells with an aspect ratio greater than 1, and LES is considered in the rest of the domain. For flows at high Reynolds numbers, the instantaneous profile of $f GWMLES$ is expected to be time-dependent due to velocity fluctuations. To take this into account and avoid that LES is solved in cells of high aspect ratios, the transition ratio of the boundary layer cells is adapted such that the aspect ratio of the outermost hexahedral cells is close to 1 (see Fig. 5). It should be noted that the mean (time-averaged) contour of the function $f GWMLES$ is obtained after a precursor RANS simulation; that is, no unsteady simulations are needed to adapt the mesh to achieve the right numerical setup. The selection between GWMLES_sub and $GWMLES log$ is done so that $f GWMLES = 0$ (LES) in a computational cell if and only if the cell aspect ratio is close to 1; in particular, based on the cited references, values lower than 2 are considered acceptable. The set of constants of Table I is not optimal for all cases; for a given mesh, $GWMLES log$ may lead to solving in RANS large regions suitable for LES, while GWMLES_sub may lead to solving in LES a region that should be solved in RANS (with the corresponding risk of GIS). In those cases, the modification of the constant $c sub$ is proposed. The effect of this modification will be assessed in one of the test cases presented in Sec. III.

It is verified that, at the end of the iterative process that is done every time step, the globally scaled residuals are $< 5 × 10 − 4$ (continuity), $< 10 − 5$ (momentum and turbulent variables), and $< 10 − 8$ (energy). These values are lower than those by default in Ansys Fluent of $< 10 − 3$ (continuity, momentum, and turbulent variables) and $< 10 − 6$ (energy). Simulations are started at t = t₀ from a precursor RANS simulation (steady-state solver with second-order upwind discretization in the convective terms), and the time step dt is selected such that the maximum cell convective Courant number Co < 0.5. The same time step is used in all simulations of the mesh independence studies. In SBES and GWMLES models, the equations of the RANS turbulence variables (k and ω) are solved every 5 time steps. Simulations are finished at t = t₁ when the flow statistics reach a statistically steady state (see Algorithm 1).

Simulations are performed in an eight-node cluster, with each node consisting of a 16-core AMD EPYC 7302 processors, which gives a total number of 128 cores. For each test case, the computational cost is obtained taking the following flow parameters as references:

• Dimensionless time step, $dt ̂$⁠: the ratio between the simulation time step dt and the characteristic time $t c = L c / v ∞ , d t ̂ = d t / t c ,$ where L_c is the flow characteristic length.

• Dimensionless domain length, $L ̂ Ω$⁠: the ratio between the length of the simulation domain $L Ω$ and the characteristic length scale, $L ̂ Ω = L Ω / L c$⁠.

• Number of iterations per time step, $N iter$⁠: for a given dimensionless time step, it is defined as the number of iterations performed by the numerical scheme before reaching the convergence conditions shown in previous paragraphs. As previously noted, the number of equations solved is not the same in all time steps for the case of SBES and GWMLES models; therefore, the value given is assumed to be the average value.

• Iteration time $t iter$⁠: computational time required to perform one iteration.

• Number of flow times during transition $N flow tr$⁠: number of times that the mean flow passes through the computational domain before the extraction of flow statistics starts.

• Number of flow times during data extraction, $N flow dat$⁠: number of times that the mean flow passes through the computational domain, while the flow statistics are extracted.

• Simulation time $t sim$⁠: time required to perform the entire simulation, including the extraction of flow statistics. It is calculated as $t sim = ( L ̂ Ω / d t ̂ ) ( N flow t r + N flow dat ) N iter t iter$⁠.

Three parameters related to the accuracy of the simulations will be considered.

• Resolution quality, $Q r$⁠: this is an estimate of the resolved turbulence spectrum² 

  $Q r = k ¯ k ¯ + k r ,$

  (10)
  where $k ¯ = ( u ′ 2 + v ′ 2 + w ′ 2 ) / 2$ is the resolved turbulent kinetic energy, and k_r is the residual kinetic energy obtained from the kinetic energy transport model⁵⁹ 

  $k r = μ t 2 C k 2 ρ 2 Δ 2$

  (11)

  for a value of the constant $C k = 0.05$ derived from simulations and experimental data.⁵⁹ It is worth noting that, with additional equilibrium assumptions, the value of C_k can be also obtained from the Smagorinsky model.⁶⁰ For a typical range of the Smagorinsky constant $C s ∈ [ 0.1 , 0.2 ]$⁠, the corresponding range of C_k is $∈ [ 0.05 , 0.12 ]$⁠, which means that the value chosen from Ref. 59 gives a conservative prediction of the resolution quality.

• Number of computational cells per integral length scale, N: since there is no universal method for estimating the resolution quality, another method will be considered that consists of computing the integral length scale, which is obtained as follows:² 

  $L i i = 1 R i i ( 0 , x , t ) ∫ 0 ∞ R i i ( e i r , x , t ) d r ,$

  (12)
  where $R i i ( r , x , t ) = ⟨ v i ( x , t ) v i ( x + r , t ) ⟩$ are the components of the two-point correlation tensor, with $e i$ referring to the unit vector in the x_i coordinate direction, and $⟨ v i ⟩$ to the time-averaged value of v_i. In this work, the spanwise integral length scale L will be calculated; therefore, i = 3, and $L = L 33$⁠. The number of computational cells per integral length scale N is used as indicator for the quality of a LES simulation;⁶¹ in particular, in a subdomain $S 1 ⊂ Ω$ of the computational domain Ω,

  $N S 1 = L S 1 Δ .$

  (13)

  An approximation of $N / Δ$ can be obtained for the particular case of a slight variation of the von Kármán spectrum.⁶² The assumption that four cells are needed for a particular wavelength—which is considered as the characteristic length of the eddies—leads to a value of $N ≈ 4$ for $Q r = 0.5$⁠, and $N ≈ 10$ for $Q r = 0.8$⁠.

• Q criterion, Q: it is defined as the second invariant of $∇ v$⁠,⁶³ 

  $Q = 1 2 ( | | W | | 2 − | | S | | 2 ) ,$

  (14)

  where $W$ is the vorticity tensor, and $| | W | | = [ tr ( W W t ) ] 1 / 2$ and $| | S | | = [ tr ( S S t ) ] 1 / 2$⁠. The iso-surfaces of Q allow the visualization of the turbulence structures.

• Taylor microscale, $λ$⁠: the length scale that represents the onset of the turbulence dissipation range, which is defined from the two-point correlation function² 

  $λ i i = 2 − ∂ 2 R i i ∂ x i 2 | ( 0 , x , t ) .$

  (15)
  The calculation of the Taylor microscale for meshes that differ more than one order of magnitude from the Kolmogorov scale can lead to large uncertainties, since the number of points close to the $r = 0$ is reduced. An approximation for the Taylor microscale will be used to estimate the cutoff limit of the turbulence spectrum of numerical solutions² 

  $λ = 10 ν k ¯ ε ¯ ,$

  (16)

  where $ϵ ¯ = ν | | S | | 2$ is the resolved turbulent dissipation rate.

• Normalized error, $ϵ$⁠: in each test case, the numerical result $g num$ will be compared with the reference experimental or DNS data $g ref$ in a subdomain $S 2 ⊂ Ω$ of the computational domain Ω. The normalized error ϵ is defined as follows:

  $| | ϵ ( g ) | | S 2 = | | ϵ ( g ) | | L 2 ( S 2 ) = | | g num − g ref | | L 2 ( S 2 ) | | g ref | | L 2 ( S 2 ) .$

  (17)

The accuracy of the GWMLES model will be compared hereafter with other models such as SBES, DES family, and URANS, as well as with results from LES/WMLES simulations available in the literature. DNS results or experimental data will be taken as reference in each case. First, the subsonic round jet is studied in Sec. III A to verify the requirement that the blending function does not activate in turbulent regions far away from walls, which is equivalent to show that it gives a LES solution in a free shear flow; second, the periodic channel in Sec. III B serves as reference for the validation of the classical WMLES mode, which is given by GWMLES_sub; third, the flow around the Ahmed body will be studied in Sec. III C, which is a complex turbulent wall-bounded flow used to test the $GWMLES log$ model; finally, the capabilities of using a transition RANS model in the $GWMLES log$ approach will be evaluated by studying the flow around a circular cylinder around the critical regime in Sec. III D.

The subsonic round jet has been extensively studied in the literature to obtain the acoustic properties of jets,⁶⁴ which are key to the design of future aircraft engines.⁶⁵ Taking experimental data as Ref. 66, simulation results will be compared with LES results obtained by the present authors, as well as with WMLES results available on the literature.⁶⁷ 

The axisymmetric domain with respect to the x axis is shown in Fig. 6. The Mach number M = 0.5; therefore, the pressure at the outlet of the convergent nozzle is equal to the ambient pressure. The model assumes that the turbulence at the nozzle outlet is negligible compared to that produced by the Kelvin–Helmholtz instability. Thus, the nozzle is not part of the computational domain, and no turbulence fluctuations are defined at the inlet. The dimensions of the computational domain have been deduced after performing a domain independence study.

The ambient pressure and temperature are $p ∞ = 101 325$ Pa and $T ∞ = 298.15$ K, while the Reynolds number $R e D = 2.8 × 10 5$ referred to D = 25.4 mm and the ambient density, which is obtained with the ideal gas law. At the inlet, the turbulence boundary conditions are referred to the RANS turbulence variables, so that the turbulent kinetic energy and the turbulent specific dissipation rate can be calculated. The turbulence intensity $I = 2 k / 3 / u ∞$ is obtained from Ref. 66. Like in the rest of the validation cases, it has been verified that the value of the inlet turbulent viscosity does not have an impact on the results. The boundary condition for the turbulent specific dissipation rate is the asymptotic solution in the viscous sublayer.

Three different meshes are considered, with $Q r$ [Eq. (10)] ranging from 0.5 to 0.8 in the region inside the edge of the shear layer. Table II shows that the estimation N = 10 [Eq. (13)] for $Q r = 0.8$ presented in Sec. II B 3 is met, so that both methods show that the resolution of the fine mesh is acceptable for a LES simulation. The normalized error decreases as the mesh size decreases. The differences between the fine and the medium mesh are lower than those between the medium and the coarse mesh; therefore, a convergent behavior is achieved.

In the fine mesh, the normalized error in the fluctuations is lower than 15%, which is related to the $Q r$ value that is reasonably uniform along the jet [see Fig. 7(b)]. Performing a Richardson extrapolation does not lead to $ϵ → 0$ for $Δ → 0$⁠, which is explained by two reasons: first, neglecting the turbulence content coming from the nozzle flow is expected to lead to errors that are important when $Q r → 1$⁠; second, as shown in Sec. III A 2 b, the experimental curve is the mean of several measurements; thus, for $Q r → 1$⁠, the uncertainty of the measurements dominates against the normalized error related to the simulations.

Comparisons with experimental data will be shown for the fine mesh in all cases. Figure 8 shows the results of the mean and fluctuating axial velocities given by different models. The mean differences between $GWMLES log$ and SBES are lower than 2%, and closer to experimental data than DDES in Fig. 8(d). Even if $f GWMLES log$ is activated far from the jet region [see Fig. 7(a)], LES results (WALE) are shown to verify that the regions with $f GWMLES log > 0$ do not affect the results. The data available in the literature are for a WMLES case that includes the nozzle region, for a fine mesh consisting of $∼ 10 2$ M elements. As expected, its accuracy is higher than the one of the simulations performed in this work. URANS gives a reasonable prediction of the mean flow, although this model significantly underpredicts the fluctuations. The normalized errors (Table II) of the URANS solution obtained with the fine mesh are larger than those given by $GWMLES log$ with the coarse mesh in the velocity fluctuations; therefore, for free shear flows, using the $GWMLES log$ model is recommended over URANS even for coarse meshes.

Using a different residual stress model in $GWMLES log$ (WALE) and DDES (Smagorinsky) does not lead to important variations in the results. However, as previously mentioned, the model constant used in the DDES model is too dissipative in free shear flows. This is translated into the resolution of coarser turbulence structures, in particular near the nozzle (see Fig. 9).

The increased dissipation produced by DDES can be verified by looking at the turbulence spectrum at $x / D = ( 10 , 0 , 0 )$ [see Fig. 10(a)], which is outside of the potential core of the jet. Both turbulence models give almost the same spectrum for high wavenumbers. However, the cutoff limit is higher in $GWMLES log$⁠, which is an important feature for acoustics calculations. This can be further verified quantitatively by calculating the Taylor microscale [see Eq. (15)]. The higher values of λ close to y = 0 in the case of $x / D = 2$ are because this zone is inside the potential core of the jet, where the development of the turbulence structures is limited. Analyzing the Q criterion, the turbulence spectrum, and the Taylor microscales leads to the same conclusion: $GWMLES log$ solves finer turbulence structures than DDES.

It has been verified that, for the same mesh, $GWMLES log$ requires 30% lower computational cost than DDES, and 10% more computing time than LES (see Table III). The reason is that, in $GWMLES log$⁠, RANS equations are solved every five time steps, while in DDES they are solved at all time steps, and no additional equation is solved for LES. Considering that the number of cores (128) is moderate, the simulation time is reasonable for an industrial framework.

The periodic channel flow is a benchmark case that provides extensive information on near-wall turbulence without the requirement of simulating a large and/or complex domain. DNS results will be taken as Refs. 68 and 69, and the accuracy of the GWMLES_sub model will also be compared with WMLES results obtained by other authors.⁷⁰ The influence of the constant $c sub$ on the performance of the model will be studied by considering two simulation approaches: the baseline with $c sub = 500$ (⁠ $GWMLES sub 500$⁠) and a variation with $c sub = 1000$ (⁠ $GWMLES sub 1000$⁠).

In what follows, the air is assumed to be incompressible. The domain is a rectangular periodic channel (see Fig. 11), with the dimensions validated in Ref. 71 against larger domains. In the longitudinal (x) direction, the mass flow rate $m ̇$ is defined from the freestream velocity $v ∞$⁠, while in the spanwise direction (z) the mass flow rate vanishes. In the no-slip walls, the wall-parallel velocity is given by the solution in the viscous sublayer.⁴⁸ The Reynolds numbers $R e = 2 ρ u ∞ h / μ = 8.7 × 10 4$ and $R e τ = ρ u τ h / μ = 2000$⁠, where $u τ$ is the friction velocity.

The results of the mesh independence study are shown in Table IV. The values given by $GWMLES sub 500$ converge toward a value whose normalized error is lower than 10% in all variables. This error is lower than 20% in the coarse mesh, even if the resolution quality is only 60%. This is attributed to the fact that this test is one of the classical benchmarks that are used to validate numerical methods and turbulence models, with no particular complex flow features such as recirculation or separation zones. The value of the (time-averaged) $y ̃ w + = 1$ is uniform due to the periodicity of the flow, and 10 cells normal to the wall are defined where $f GWMLES = 1$⁠. Following the trends shown before, solving 80% of the turbulence spectrum corresponds to defining approximately ten computational cells per integral length scale. It can be concluded that both mesh quality indicators are related, with a relationship $N ( Q r = 0.8 ) ≈ 10$⁠. Thus, in the following mesh independence studies, only the resolution quality will be specified.

As expected from the reasoning followed in Sec. III A 1, most of the domain is solved in LES mode [see Fig. 12(a)]. Details of the transition from RANS to LES will be given below. The resolution quality is uniform in all meshes considered [see Fig. 12(b)].

The comparison between velocity statistics in the boundary layer given by different simulations is illustrated in Fig. 13. The agreement with DNS data is excellent in the GWMLES_sub and SBES models, while IDDES predicts a 50% lower thickness of the log-layer [see Fig. 13(a)]. The WMLES results⁷⁰ are obtained with a wall function approach and provide great agreement in the log-layer, although errors up to 20% in the mean dimensionless x-velocity are noticed in the viscous sublayer and the buffer layer.

The deviations between different methods are more pronounced in the fluctuating variables. GWMLES_sub, SBES, and IDDES results differ up to 20% in Fig. 13(b), although they are significantly more accurate than WMLES [see Figs. 13(b)–13(d)]. The mesh is coarser in the simulation of Ref. 71 (⁠ $Δ min + = 79$⁠), but this is partially balanced by the fifth-order spectral difference method that is used. Independently of the mesh and the numerical method, the issues of the wall function approach have been remarked in the introduction. This is relevant, since the fluctuations in the buffer layer are of the same order of magnitude as those in the log-layer. Even if RANS models are calibrated for steady cases, it is demonstrated that the model developed in this work is able to capture the velocity fluctuations that are produced in the viscous sublayer and the log-layer, while the wall function approach leads to errors that are as high as 50% in Fig. 13(b). Solving the viscous sublayer and the buffer layer in RANS leads to a higher computational cost than with the wall function-based approach. Nevertheless, in these zones, the cell aspect ratio is much higher than 1; therefore, both mesh sizes are of the same order of magnitude.

URANS simulations have also been performed, leading to a deviation lower than 1% in the mean velocity, and with negligible fluctuations. Starting from a precursor RANS simulation, the high dissipation of URANS leads to no fluctuations. The same conclusion is extracted after sufficient time if the URANS simulation is initialized from a GWMLES_sub solution. This could be anticipated before performing the simulation, since all parts of the domain are inside the boundary layer, with no relevant flow instabilities. Given that in GWMLES_sub, the normalized error of the mean values obtained with the coarse mesh is lower than 4% (see Table IV), and that URANS predicts negligible fluctuating variables, GWMLES_sub provides more accurate results than URANS in the periodic channel flow for all meshes considered.

The increase in $c sub$ defined in Eq. (7) leads to a RANS-LES transition at higher $y ̃ +$ (see Fig. 14). The transition region of the mean value of the blending function contains more than two cells. However, as illustrated in Fig. 12(a), the fluctuations in the function are important because of the high resolution quality, which means that the time-averaged value is smoothed with respect to the instantaneous value; that is, in Fig. 14, the transition region is of the order of $y ̃ +$⁠, while Fig. 12(a) shows that the instantaneous transition region is much thinner. The mean SBES blending function is essentially the $GWMLES sub 500$ function translated approximately by $y ̃ + = 5$⁠. In $GWMLES sub 1000$⁠, the $y ̃ +$ in which RANS-LES transition takes place is twice the one of $GWMLES sub 500$⁠, with a lower slope of the function. No LLM is noticed in Fig. 13(a), and the velocity fluctuations are underpredicted, which is attributed to the higher dissipation as a consequence of the extended RANS region. The great accuracy of the GWMLES_sub, SBES, and IDDES models is also proved in the mean skin friction coefficient, which is predicted with an error lower than 5% in all models (see Table V). This study confirms the capabilities of the proposed model to improve the current wall function and mixing-length approaches for WMLES.

The IDDES model has been particularly calibrated for periodic channel flows.⁶ Moreover, there are no transition regions, which means that the transition from modeled-to-resolved turbulence is not as relevant as in other more complex wall-bounded flows. This can be verified in the size of the resolved turbulence structures, in which no qualitative differences can be observed between GWMLES_sub and IDDES (see Fig. 15). Quantitative results derived from the distribution of the Taylor length scale lead to the same conclusion (see Fig. 16).

The computational cost for GWMLES_sub (same as for SBES) is around 35% lower than the one for DDES (Table VI). The reduction is greater than in the subsonic round jet case, since in compressible flow RANS equations represent $2 / 7 = 28 %$ of the total number of equations, while in incompressible flow without heat transfer the ratio is $2 / 6 = 33 %$⁠. The computational cost is of the same order as the one of the simulations of the subsonic round jet (see Table III). The higher cell count and lower time step, as well as the large number of time steps required for the turbulence structures to be fully developed, are compensated by the reduced numbers of iterations per time step and dimensionless length of the computational domain.

The Ahmed body is a model of a generic road car⁷² (see Fig. 17). The flow shares some features with the real cars such as a large three-dimensional recirculation zone at the back of the geometry, and the interaction of two counter-rotating vortices that are generated in the slant edges. For slant angles lower than the critical angle of 30°, the flow separates in the beginning of the slant, and it reattaches downstream, leading to complex flow features. Taking experimental data as Refs. 72–75, simulation results will be compared with the results obtained from LES and WMLES simulations available in the literature.^76,77

The computational domain is a 3/4 open wind tunnel with the dimensions of the cross section obtained from Ref. 73 (see Fig. 18). The lengths upstream and downstream of the Ahmed body have been obtained after performing a domain independence study. The Reynolds number $R e L = 2.7 × 10 6$⁠, and the turbulence intensity $u ′ / u$ at the inlet, which is used to calculate the turbulent kinetic energy and the turbulent specific dissipation rate, is obtained from Ref. 73.

Three different meshes are considered, with $Q r$ in the slant region ranging from 0.5 to 0.8. Coarser meshes lead to numerical instabilities when the low-dissipation discretization methods mentioned in Sec. II B 2 are used. Table VII shows that the error in the velocities of the boundary layer decreases as the mesh size decreases, and a convergent behavior is achieved.

The drag parameters—with $⟨ C d , f ⟩$ referring to the friction component of the drag coefficient—are not compared with experimental data; the reason is that there is no consensus based on different experiments. In the normal direction of the boundary layers, 20 cells are defined, with $y ̃ w + < 2$ in the turbulent boundary layers [see Fig. 19(a)]. These conditions imply small differences in the quantities that strongly depend on the RANS model, such as the friction component of the drag coefficient. In contrast, there are separation and reattachment processes, and the drag coefficient and the fluctuating values are mainly determined by the accuracy of the LES region. This is the reason why these variables are more sensitive to mesh refinement. In the fine mesh, $Δ min < 0.1 δ$⁠, even though, as noted in Ref. 20, the flow is not exposed to an adverse pressure gradient on the roof; therefore, GIS is not a problem for the present geometry.

It should be noted that the normalized error of the fluctuating velocity of the fine mesh is still 26%. Two universal reasons are related to the use of a RANS model in the boundary layers, as well as the limited resolution of the turbulence spectrum. Even if $Q r ≈ 0.8$ in the slant, there are some regions in the wake where the value is $0.6 − 0.7$ [see Fig. 19(c)]. In addition, in the experimental data available, there are only ten data points in the boundary layer, which are only spaced by 4 mm. This means that the accuracy of the measuring technique could have an important role in the comparison, as well as the manufacturing tolerances of the body. In Sec. III C 2 c, the results will be compared with data available on the literature, so that a better picture can be obtained in regard to the accuracy of the model proposed.

The mean profiles of the x-velocity and the $GWMLES log$ blending function in different positions of the roof of the Ahmed body are presented in Fig. 20. The equilibrium solutions in the viscous sublayer $⟨ u sub + ⟩ = y ̃ +$ and the log-layer $⟨ u log + ⟩ = C + + log ( y ̃ + ) / κ$ for values of the von Kármán constant $κ = 0.41$ and $C + = 5$ (Ref. 48) are plotted as reference. Figure 20(a) shows that the boundary layer is not fully developed in the first three reference points. The transition point from RANS to LES depends on the x-coordinate [see Fig. 20(b)], since the log-layer thickness is proportional to the distance to the leading edge when the flow is not fully developed. At x = −250 mm, the profile obtained from the numerical simulations matches the law of the wall. The mesh size in the roof region is $Δ ≈ 0.05 δ$ (half of the size used in the slant region), and the transition from RANS to LES takes place in around two cells. This is the expected result, given that the reference considered in Sec. II A 1 was 1 cell for $Δ ≈ 0.1 δ$⁠. In the fully developed case (x = −250 mm), the transition point is located near the end of the log-layer, as predicted.

Figure 21 compares the streamwise velocity statistics in the slant region of the symmetry plane (see Fig. 17) with experimental data, LES, and WMLES simulations. The results of LES and WMLES simulations are those obtained with the fine mesh used in the corresponding work. The differences between $GWMLES log$⁠, SBES, and DDES are lower than 5% at almost all points, especially in the mean velocity. These models are more accurate than WMLES in the mean values, although the deviation from the experimental data is more pronounced than in the LES simulations. Even if the accuracy of all models compared is similar in the fluctuating velocity, at x = −203 mm, there are some fluctuations predicted by the $GWMLES log$⁠, SBES, and LES models that are not present in the DDES model and in the experimental data. URANS fails to determine the velocity profiles, with negligible temporal fluctuations; indeed, the normalized errors (Table VII) are greater than those obtained with $GWMLES log$ and the coarse mesh.

Figure 21(a) allows the detection of the separation and reattachment point from different models. In the attached flow, the streamwise velocity component is positive, while the opposite occurs when the flow is separated. Except from URANS, in which the flow is separated along the entire slant, the separation at the beginning of the slant and the subsequent reattachment before the end are predicted by all models. However, $GWMLES log$⁠, SBES, and DDES models overpredict the width of the separated zone, while LES and WMLES provide the opposite behavior. It should be noted that a positive streamwise velocity component is not a sufficient condition for the flow to be attached, and it has been verified that in these cases, the vertical (z) velocity component is negative, as required to deduce that the flow is attached.

The effect of the prediction of the transition point will be analyzed in the drag parameters. Among the models that give a reasonable prediction of the velocity statistics, since the separated region is the largest in the DDES results, it is expected that DDES provides the highest drag coefficient value. In contrast, the lowest should be given by the LES model given the overprediction of the attached region by this model. This is indeed shown in Table VIII. As previously anticipated, the value of the drag coefficient is not clear taking previous articles as references. The quality of the test sample has a strong influence on the drag coefficient;⁷⁵ it has been shown that small roundings in the upper edge of the slant lead to a decrease in the drag coefficient, and, in particular, the drag coefficient does not depend on the radius of curvature in the cases studied. That is, there is a reference value in the body with a sharp edge, 0.356, which is similar to the one obtained in Ref. 74. The drag coefficient reduces to $≈ 0.297$ with a rounded edge, which is closer to the one of Ref. 72. Another key parameter is the friction component of the drag. The values given by $GWMLES log$⁠, SBES, and DDES models differ by less than 10% and agree with the value given in Ref. 72. LES results provide a low value of the friction component, which may be a reason for the discrepancy in $⟨ C d ⟩$ with respect to Refs. 74 and 75. Since friction effects are dominant in the attached flow, the friction component given by URANS is in line with the one given by $GWMLES log$⁠, SBES, and DDES. Given the large deviations of the URANS result with respect to experimental data in the velocity statistics, it is expected that the friction coefficient predicted by URANS deviates significantly, although this is not the case. Furthermore, flow variables need to be considered to explain this behavior.

Even if the drag coefficient is not obtained in the experiments taken as reference in Fig. 21 (Ref. 73), the pressure coefficient distributions in the slant and the rear surface of the Ahmed body are available. Figure 22 compares numerical and experimental results in the slant and shows that simulation values provide an accurate overall picture of the distribution. However, errors up to $≈ 20 %$ are obtained in the numerical simulations (see Table IX). As expected from the comparison of velocity statistics, the prediction of URANS is of reduced accuracy, since this model does not predict clearly two of the main flow physics: the separation and subsequent reattachment, as well as the two counter-rotating vortices, which are extremely weak [see Fig. 22(a)]. It has been verified that the normalized errors obtained with the coarse mesh in GWMLES $log$ are lower than those of URANS in the fine mesh, which leads to the same conclusion as in the previous cases: GWMLES provides more accurate results than URANS even in coarse meshes.

The contours of the pressure coefficient are less accurate at first sight in the rear of the body, which corresponds to the plane x = 0 (see Fig. 23). However, the area-averaged value is highly accurate in all cases, even in URANS (see Table IX); the latter is considered to be a coincidence, since the pointwise contour of URANS [Fig. 23(b)] deviates significantly from the experimental data [Fig. 23(a)]. Assuming that there is no error in the remaining surfaces, from the results of Table IX, it is deduced that $GWMLES log$ underpredicts the drag obtained in the experiments of Ref. 73. This enforces the assumption that the correct drag coefficient is 0.34, so that the hybrid RANS/LES models would provide a better accuracy than the LES model for a much lower computational cost. In that situation, $GWMLES log$ model is more accurate than the SBES, which can be due to solving the outer layer in LES mode instead of in RANS.