- Email: [email protected]

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A dynamic delayed detached-eddy simulation model for turbulent ﬂows Chuangxin He a,b, Yingzheng Liu a,b,∗, Savas Yavuzkurt c a

Key Lab of Education Ministry for Power Machinery and Engineering. School of Mechanical Engineering, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China b Gas Turbine Research Institute, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China c Department of Mechanical and Nuclear Engineering, Pennsylvania State University, University Park, PA, United States

a r t i c l e

i n f o

Article history: Received 29 June 2016 Revised 3 December 2016 Accepted 20 January 2017 Available online 22 January 2017 Keywords: Dynamic DDES LES Separated ﬂow Turbulent wall heat transfer

a b s t r a c t The current study developed a new dynamic delayed detached-eddy simulation (dynamic DDES) model based on the k-ω SST model and the well-established dynamic k-equation subgrid-scale model. Instead of using a constant model coeﬃcient CDES in traditional DES formulations, the present model employs two coeﬃcients Ck and Ce , which are computed dynamically by taking into account the spatial and temporal variations of the ﬂow ﬁeld at the grid and test ﬁlter levels. A modiﬁcation on shielding function fd is proposed, with a spatial uniformization operator imposed on the velocity gradient to obtain a smooth and monotonous hybrid interface. A damping function ϕ d is introduced based on the local grid resolution and ﬂow condition to damp the Reynolds-averaged Navier–Stokes (RANS) region and achieve wall-modeled LES (WMLES) mode dynamically. The test of the model in developed channel ﬂow shows the log-layer mismatch (LLM) problem is signiﬁcantly improved with respect to the dynamic LES model and original DDES model. The use of the spatial uniformization operator and the damping function convincingly demonstrates the improvement in prediction of separated ﬂows, with the model coeﬃcients dynamically computed. The LES region is maximized at the limit of grid resolution and more turbulent vortical structures are resolved. The test in the ribbed channel ﬂow shows the present model has considerably better performance in prediction of the mean and turbulence velocity in the strong shear layer and the recirculation bubble. In addition, the simulation of impinging jet shows the model exhibits rapid switching from the RANS to LES under the ﬂow instabilities when the inﬂow does not include turbulence content. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Hybrid Reynolds-Averaged Navier–Stokes/Large-Eddy Simulation (RANS/LES) models have been developed to cope with the well-known shortcomings of RANS models in separated ﬂows and the high expense of LES models with high-Reynolds number wallbounded ﬂows. To this end, there exist two types of RANS-LES hybrid methodology, wall-modeled large-eddy simulation (WMLES) and detached-eddy simulation (DES). In WMLES models, RANS is applied only to a thin layer of the near-wall region where the wall distance is much smaller than the boundary layer thickness [1–3]. This technique captures most of the turbulence inside the boundary layer using LES, without resolving the smallest and highly Reynolds-number-dependent turbulent structures above the viscous sublayer. In DES models, a major part of the attached boundary layer is modeled by RANS, while LES is only applied in the ∗

Corresponding author. E-mail address: [email protected] (Y. Liu).

http://dx.doi.org/10.1016/j.compﬂuid.2017.01.018 0045-7930/© 2017 Elsevier Ltd. All rights reserved.

separated ﬂow regions [4–6]. In the vicinity of the wall, the grids are highly anisotropic, with large grid spacing in directions parallel to the wall. RANS is used to predict the averaged ﬂow quantities of the turbulent boundary layer, as the grid resolution is insufﬁcient for LES. For the original DES model (DES97) [5], however, the problematic behavior of Grid Induced Separation (GIS) was observed [7] due to the Modeled Stress Depletion (MSD) [8] when the grid was reﬁned to shift from RANS to LES without balancing the reduction of eddy viscosity by the resolved turbulence content. Accordingly, delayed DES (or DDES) was proposed to prevent this switch in the boundary layer region by using a generic formulation of the shielding function which depends on the eddy-viscosity, the local velocity gradient and the wall distance [8,9]. Two perspectives of the DES model are the WMLES capability in ﬂow conﬁgurations with enough turbulence in the inﬂow, and the rapid RANS-to-LES transition when a steady inlet boundary condition is imposed. However, simulation of the channel ﬂow [10,11] using the DES model identiﬁed a Log-Layer Mismatch (LLM) problem, which was also found in the DDES model. As Shur

C. He et al. / Computers and Fluids 146 (2017) 174–189

et al. [12] claimed, the simulations produced two logarithmic layers with different intercepts, resulted in under-prediction of the skin friction. In addition, DDES usually produced a thick RANS region in the ﬂows where WMLES is highly desirable. In this regard, the improved DDES (IDDES) was developed with a series of blending functions and a modiﬁed ﬁlter width [9,12]. This model improved LLM and ensures an automatic selection of DDES or WMLES mode depending on the turbulence content in the ﬂow. Motivated by the exploration of an alternate and simpler DES formulation which could overcome both the MSD and LLM, Reddy [13] recently modiﬁed the DDES formulation, in which the eddy viscosity was similar to the Smagorinsky formula far from the wall. With a hybrid ﬁlter width deﬁned, the problem of LLM was signiﬁcantly improved. When a steady inlet boundary condition is used, DES relies on the natural instability of the ﬂow to generate resolved turbulence. It is well known that the original DES model delays the RANS-to-LES transition in the free and separated shear layers [14], largely attributed to the conservative estimation of the ﬁlter width max which overpredicts the subgrid eddy viscosity and kills the turbulence content. This problem is more serious when the grids are strongly anisotropic, which usually occurs in the initial region of the free and separated shear layers. Chauvet et al. [15] modiﬁed the ﬁlter width deﬁnition by applying max only on the 2D plane which was perpendicular to the local vorticity vector. Further improvement were made by introducing solution-dependent kinematic measures to reduce the subgrid eddy viscosity in the initial region of the shear layers [16], and employing an alternative subgrid eddy viscosity formulation in the LES region which was expected to discern between quasi 2D and developed 3D ﬂows and decrease the eddy viscosity in the early shear layers [17]. In this contribution we introduce a dynamic DDES model formulation. This is designed to be suitable for application to a wide range of ﬂow conﬁgurations. In DDES models [6,8,9], the model coeﬃcient CDES was usually calibrated using the benchmark test of decaying isotropic turbulence (DIT). However, it was demonstrated that the calibrated value strongly depended on grid resolution [18]; a smaller CDES was required when increasing the mesh resolution. A constant CDES in the whole computational domain was not suitable, however, as the turbulence was inhomogeneous in wall-bounded ﬂows. In this regard, a dynamic version of DDES was proposed [19] by dynamically determining the model coeﬃcient CDES based on the local ﬂow, demonstrating relatively accurate prediction of the fully developed stationary channel ﬂow and rotating channel ﬂow in comparison with traditional DDES models (constant CDES ). However, in Yin’s model, the shielding function fd deﬁned by Spalart et al. [8] was used to hybridize RANS and LES based on spatial locations, which would introduce inaccurate prediction of complex ﬂows. Recall that even in the constantcoeﬃcient DDES models [9], signiﬁcant decrease of the shielding function was found in the free-stream region of the ﬂat-plate boundary layer, forming RANS islands surrounded by the LES zone. Gritskevich et al. [9] argued that this stemmed from the low values of strain rate magnitude but should not cause any problems as it happened outside the boundary layer. However, this effect may not be negligible in complex ﬂows when the model coeﬃcient is dynamically computed as smaller vortex structures are captured and give rise to wide variation of the velocity gradient, leading to dramatic fd decrease and unphysical eddy viscosity increase in ﬂow separation region. Although Yin’s model has an intrinsic capability to adjust the hybrid interface [19], the shielding of the RANS region is still somewhat conservative. Accordingly, it is highly desirable to establish a dynamic procedure to maximize the LES region in response to variations in local ﬂow structures and grid resolution. Inspired by the dynamic strategy [19], in this study, a new formulation of dynamic DDES model based on the k – ω SST turbulence RANS model [20] and the dynamic k – equation subgrid LES

175

model [21] is developed. The model coeﬃcients are dynamically computed along with a dynamic adjustment of the RANS area to achieve the WMLES. The same form of the k equation is used in both the RANS and LES regions, while the ω equation is only active in the RANS region. The shielding function fd in the constantcoeﬃcient DDES model [9] is modiﬁed by imposing a spatial uniformization operator on the velocity gradient to achieve a smooth hybrid interface. In addition, a damping function based on the local velocity gradient and grid resolution is applied in fd to activate the WMLES mode in various ﬂow conditions. Accordingly, the constant in the deﬁnition of the shielding function is recalibrated using the ﬂat-plate boundary layer. The open source code OpenFOAM (http://openfoam.org) has been used for all the present simulations. Gaussian ﬁnite volume integration with central differencing, or hybrid central-upwind differencing along with explicit correction based on the gradients, for interpolation, is used for spatial discretization of equations. A second-order backward difference scheme is selected for time integration. A PISO algorithm with two inner iterations is used for the velocity-pressure coupling in the N-S equations. Four test cases, i.e., developed channel ﬂow, ﬂow over periodic hills, ribbed channel ﬂow and impinging jet with heat transfer, are used to evaluate the model. 2. Model formulation 2.1. Hybrid method The two existing models that form the basis for the present dynamic DDES model are reproduced here for convenience, i.e., the k – ω SST RANS model [20], and the dynamic k – equation subgrid LES model [21]. The k – ω SST RANS model [20] is as follows:

∂ρ k μt ρ k3/2 ∇ k + Pk − √ + ∇ · (ρUk ) = ∇ · μ + ∂t σK3 k/Cμ ω

(1-a)

∂ρω μt ∇ k · ∇ω ∇ω +2(1 − F1 )ρ + ∇ · (ρ U ω ) = ∇ · μ + ∂t σω 3 σω 2 ω ω 2 +α3 Pk − ρβ3 ω (1-b) k

Here, the dissipation term in the k equation is rewritten for convenience of the DDES formulation. If we deﬁne the RANS length scale as

√

LRANS =

α1 k max(α1 ω, F2 S )

(2)

the turbulence eddy viscosity becomes

μt = LRANS k

(3)

The model constants are identical with those deﬁned in the original model [20]. The dynamic k – equation subgrid LES model [21] is as follows:

∂ρ k μt ρ k3/2 ∇ k + Pk − + ∇ · (ρUk ) = ∇ · μ + ∂t σK3 /Ce

(4)

When we deﬁne the LES length scale

LLES = Ck

(5)

the subgrid eddy viscosity becomes

μsgs = LLES k

(6)

In this model, the model coeﬃcients Ck and Ce are dynamically computed as

Ck =

1 Li j · Mi j 2 Mi j · Mi j

(7-a)

176

C. He et al. / Computers and Fluids 146 (2017) 174–189 -1

10

t=0 -2

10

E

t = 0.87 t=2 -3

10

Exp. by Comte-Bellot & Corrsin 3 64 3 32 3 16

-4

10

10

0

10

1

10

2

k Fig. 1. Energy spectra of the decaying isotropic turbulence calculated with the dynamic procedure.

Ce =

S2 − Sˆ2 (μ + μt ) K K 3/2 /

(2)

(7-b)

lles =

(9-b)

Ce

When the shielding function fd is introduced to hybridize the models, the DDES length scales are determined as

with

˙U Li j = U U −U i˙ j˙ i j˙

(8-a)

LDES = fd LLES + (1 − fd )LRANS

(10-a)

√

Mi j = −2 K Si j

(8-b)

ldes = fd lles + (1 − fd )lrans

(10-b)

˙U KK = U U −U i˙ i˙ i i˙

(8-c)

Based on the new deﬁned length scales, the dynamic DDES model is

In these equations, U is the Reynolds-averaged velocity tensor or the ﬁltered velocity tensor at the grid level in corresponding models. “ ˆ ” denotes the ﬁlter operator at the test level. Note that in the Smagorinsky model, the dynamically-determined model coeﬃcient is averaged either in the homogeneous direction [22] or along the ﬂow pathlines [23] to avoid numerical instability problem. In the present formulation, the average is not required due to the inherent smoothing of the transport equation. However, this can yield locally negative values of Ck and Ce , then clipping of the negative values is applied straightforwardly in the current dynamic procedure. A prior test of the dynamic procedure is illustrated in Fig. 1, where the energy spectra of the decaying isotropic turbulence are presented. In this assessment, the simulation is performed on grids with resolution 163 , 323 and 643 using the dynamic k–equation LES model (corresponding to DES with RANS deactivated). The experimental data of Comte-Bellot and Corrsin [24] is employed for comparison. The spectra at non-dimensional time t = 0 is contained by the initial ﬁeld, while good agreement of the spectra is found on different grids at t = 0.87 and 2. Note that the dynamic model needs a certain gird resolution for the test ﬁlter to be valid, the spectra on the 163 grid is slightly overpredicted. For more details about this dynamic procedure, please refer to [25]. Since the forms of the k equation in both models are identical except the length scale in the dissipation terms, we deﬁne two other length scales in the RANS and LES equations as

lrans

√ K = Cμ ω

(9-a)

∂ρ k μt ρ k3/2 + ∇ · (ρUk ) = ∇ · μ + ∇ k + Pk − ∂t σK3 ldes

(11-a)

∂ρω μt ∇ k · ∇ω ∇ω +2(1 − F1 )ρ + ∇ · (ρ U ω ) = ∇ · μ + ∂t σω 3 σω 2 ω ω 2 +α3 Pk − ρβ3 ω (11-b) k

and the eddy viscosity becomes

μt = LDES k

(12)

Hence, the dynamic DDES architecture is constructed by using the same equation in both RANS and LES branches. In the eddyresolving region, a standard LES model with dynamically calculated coeﬃcients Ck and Ce is achieved which accurately predicts the turbulence in the detached ﬂows. In the original DDES formulation, the ﬁlter width is deﬁned as the maximum length scale of the local grid hmax . To recover the dynamic LES model in the eddyresolving region, the cubic root of the cell volume V1/3 should be used. Inspired by Reddy et al. [13], a hybrid ﬁlter width is deﬁned as

= fdV 1/3 + (1 − fd )hmax

(13)

In the eddy-resolving region where fd = 1, this expression gives a ﬁlter width V1/3 for LES. Reddy mentions that the difference lies only in the region where 0 < fd < 1, which helps ameliorate the problem of LLM. The deﬁnition of the ﬁlter width to accelerate the RANS-to LES transition suggested by Mockett et al. [16] and Shur

C. He et al. / Computers and Fluids 146 (2017) 174–189

177

variation of the velocity gradients. This will give rise to a twisted hybrid interface, and undesired growth of the eddy viscosity in the far-wall region with the resolved eddies, as sketched in Fig. 2(a). Here, a spatial uniformization operator is imposed on the velocity gradients to obtain a smooth and monotonous connection between RANS and LES

Ui jUi j = Ui jUi j ·

Ui jUi j

Ui jUi j

14 (15)

The angle bracket “” denotes the average over the whole computational domain. A schematic example is shown in Fig. 2(b), where a sequence of random data is given. After the uniformization, a new set of data is obtained with a large shrinkage on the ﬂuctuation but with the variation tendency and average value remaining almost the same. The modiﬁed shielding function is deﬁned as

fd = 1 − tanh (Cd1 rd )Cd2

(16-a)

with

rd =

Fig. 2. (a): Sketch of the hybrid interface and (b): the effect of the spatial uniformization operator.

and Spalart [17] is not employed in the current version of the dynamic DDES model. Fortunately, we will also see the rapid RANSto-LES transition in the simulation of impinging jet (Figs. 20 and 21) using the cubic root of the cell volume in LES region. This is attributed to the much smaller eddy viscosity calculated by the dynamic procedure. In the S-A based DDES model [8], the shielding function is deﬁned as

fd = 1 − tanh (Cd1 rd )Cd2 (Cd1 = 8, Cd2 = 3 ) rd =

k/ω + υ κ 2 dw 2 Ui jUi j

(14-a) (14-b)

The two constants Cd1 and Cd2 are based on the test of the model in the ﬂat-plate boundary layer. These values ensure that the whole or at least a large portion of the boundary layer is simulated using the RANS model while the model smoothly turns to LES at the edge of the boundary layer. When this shielding function is applied in the SST-based DDES model, switching inside the boundary layer to LES is not completely eliminated. Gritskevich et al. [9] argues that Cd1 = 20 should be used to provide enough shielding for the RANS mode inside the boundary layer. A significant decrease of the shielding function was also noticed both in SST based and S-A based DDES, which was explained by low values of the velocity gradients near y/δ ≈ 1.2 where the area was outside of the boundary layer and expected to have no negative concerns (Fig. 2). Indeed, this behavior would not increase the eddy viscosity outside of the boundary layer and cause signiﬁcant problem for most of the external ﬂows as the ﬂow was almost uniform outside of the boundary layer. When the dynamic procedure is applied with the WMLES mode activated, smaller vortices near the wall will be resolved and captured by the local grid, resulting in large

k/ω + υ

κ

2d 2 w

(16-b)

Ui˙ j˙Ui˙ j˙

The constants in the fd need to be recalibrated as the ﬁeld of rd varies. We keep the same value of Cd2 = 3 and only calibrate Cd1 , which is associated with the thickness of the RANS region, using the ﬂow over a ﬂat plate. The grid used here is uniformly distributed in a streamwise direction such that the maximum grid size hmax near the wall gradually changes from hmax δ (Rex = 0) to hmax ≈ 0.2δ (Rex = 107 ) as the boundary develops. This corresponds to an “ambiguous” grid spacing [8] where MSD may occur. In this model, a much larger value of Cd1 should be used to provide enough shielding for the boundary layer. The curve in Fig. 3(a) with Cd1 = 80 shows that switching to LES is delayed to the border of the boundary layer near y/δ = 1, and is quite smooth and monotonous compared with the original formulation. In contrast, the shielding provided by the original formulation [9] with Cd1 = 8 is insuﬃcient with a slight decrease of fd near y/δ = 1.2. When the recommended value Cd1 = 20 [9] is used, the switching is delayed to y/δ = 1. However, the LES mode is not activated maturely even beyond y/δ = 2 with a substantial decrease of fd outside the boundary layer. This undesired behavior of fd is attributed to the hump of rd occurring in the range of 1 < y/δ < 1.5 (Fig. 3(b)). In addition, rd does not decrease to zero even at y/δ = 2. The modiﬁed formulation rd shows a smooth and monotonous decay inside the boundary layer, and remains zero beyond y/δ = 1 due to the spatial uniformization operator which increases the velocity gradients in the freestream region. When approaching the wall, rd grows rapidly and becomes much larger than that in the original formulation, because much smaller velocity gradients near the wall are obtained using the spatial uniformization operator. Fig. 4 shows the contours of the shielding function for both original SST based DDES [9] with Cd1 = 20 and the present dynamic DDES model. A large number of islands with substantial fd decay are seen outside the boundary layer for the original DDES model (Fig. 4(a)), while a mature switching to LES is obtained (Fig. 4(b)) when the spatial uniformization operator is applied. 2.2. RANS damping Excessively conservative shielding of the DES model would impair the turbulence resolving capability of the DDES model in separated ﬂow regions. Alternatively, damping of the shielding function is expected when the grid resolution is ﬁne enough for LES, especially in ﬂow with an abundance of turbulence. Similar to the

178

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 3. Comparison between DDES and dynamic DDES in ﬂat plate boundary layer: (a) shielding function fd , (b) rd quantity.

the maximum grid aspect ratio be 20 below which the damping can be activated. The value of ξ is set according to various tests to obtain a reasonable RANS region near the wall and the proper response to the grid reﬁnement. r gives the damping coeﬃcient between 0 and 1 based on the value of y+ . Generally, y+ relocal local mains signiﬁcantly large close to the wall due to high wall shear stress, and decays away from the wall. However, y+ does not local always decay smoothly and monotonously (which is compulsory when applied directly on the shielding damping) from the wall to the mainstream in complex ﬂows. Hence, the local average is used in the near-wall region. For convenience, we deﬁne

f = tanh Cdf rd and

ϕd =

Fig. 4. Contours of shielding function fd in ﬂat plate boundary layer computed by different models.

non-dimensional grid spacing on the wall (y+ = ∂dU/√∂νy ), the local w

non-dimensional grid spacing, which is deﬁned based on the velocity gradients, the cubic root of the cell volume and the molecular viscosity, serves as the criterion to judge whether the local grid is ﬁne enough for LES. In addition, another criterion to consider is the grid aspect ratio. Thus,

y+ local

hmax = max , ηhmin

ξV

Ui, j √ 1/3

υ

r = min 1, max y+ , 1 −1 local

(η = 20, ξ = 5)

(17-a)

(17-b)

In this formulation, y+ is deﬁned as the maximum value of local the grid aspect ratio (deﬁned by the maximum grid length scale hmax , the minimum length scale hmin and the constant η) and the local non-dimensional grid spacing. The value of η is set to allow

f · r f

3

Cdf = 60

(18-a)

2 (18-b)

In Eq. (18-a), rd is deﬁned in Eq. (16-b). When the wholedomain average operator “ ” is used, the damping function ϕ d calculates the approximate local average of r in the near-wall region speciﬁed by f. The constant Cdf = 60 is chosen to give a relatively thin layer of region to perform the average and thus provides a conservative damping function ϕ d which will be used in the new shielding function deﬁnition. The square in the expression of ϕ d is included to cope with the rapid increase of rd in the new shielding function when approaching the wall. In the present formulation, the variation of the domain average f · r with respect to the domain size is largely compensated by f Thus, the deﬁnition of ϕ d is applicable to various ﬂow conﬁgurations. Finally, the shielding function in the present dynamic DDES model is deﬁned as

fd = 1 − tanh (Cd1 rd ϕd )Cd2 (Cd1 = 80, Cd2 = 3 )

(19)

With this compact formulation, the damping strategy is more straightforward than IDDES model. In the boundary layer of the ﬂows without inﬂow turbulence content, y+ is signiﬁcantly large local and hence gives no damping on the shielding function until very ﬁne grid resolution is used to activate the WMLES mode. Fig. 5 shows the distribution of skin friction on the ﬂat plate predicted by different models. The grid is identical with the one used in the calibration of Cd1 . The original SST based DDES model [9] with Cd1 = 8 underpredicts the skin friction as the MSD occurs when the

C. He et al. / Computers and Fluids 146 (2017) 174–189

179

Fig. 5. Distribution of skin friction coeﬃcient Cf over the ﬂat plate computed by different models.

boundary layer is not fully shielded; Increasing Cd1 to 20 can solve this problem. The present dynamic DDES model is shown to adequately shield the boundary layer, giving excellent agreement with that obtained using the SST RANS model. 3. Test cases We turn to a comparative view of the performance of the present dynamic DDES model with respect to the original DDES model [9] and dynamic LES models. The previous results for ﬂow over the ﬂat plate with appropriate model constants (Fig. 3 and 5) are free of the MSD issue in the boundary layer simulation. Subsequently, four cases are tested for the WMLES capability of the present dynamic DDES model: developed channel ﬂow, ﬂow over periodic hills [26,27], ribbed channel ﬂow [28] and impinging jet with heat transfer [29]. These particular cases are chosen to eliminate the effect of the inlet boundary condition on the simulation results by using the periodic inlet/outlet boundary condition or the fully developed pipe ﬂow condition obtained by a precursor RANS simulation. In the last case, the turbulent Prandtl number model was used in the LES region derived based on the method proposed by Moin et al. [30] and hybridized with the constant hypothesis in the RANS region under the DDES framework. 3.1. Developed channel ﬂow This test is the most important for evaluation of the new model’s capability when applied to non-separated turbulent boundary layer ﬂows. As mentioned previously, both original DES [10] and DDES [8,9] exhibit unexpected LLM problems, containing two log-layers with different intercepts. The IDDES [12] corrects LLM by using a series of blending functions, but it is rather complicated. The current dynamic DDES model tends to avoid LLM by using a simpler formulation which is more akin to DDES rather than IDDES, but resolves as much turbulence as possible within the limit of grid resolution. The size of the computational domain used in the current test was Lx = 4H, Ly = 2H and Lz = 2H. Here H denoted the half-height of the channel. The grid resolution at different Reynolds numbers is shown in Table 1. Note that the dimensional grid size was kept the same for all the considered Reynolds numbers. The near-wall ystep was y/H = 5.2 × 10−4 and the stretching factor of the wallnormal step was 1.1. The resultant total grid points were around 1.1 million. The time step t was chosen to ensure that the maximum local Courant-Friedrichs-Lewy (CFL) number ≈ 1.2. The ﬂow was driven with a constant pressure gradient taken into account in the governing equation via a source term in the momentum equation. The periodic conditions were imposed both in the spanwise direction z and the streamwise direction x. Fig. 6 shows the non-dimensionalized velocity proﬁles along with the fd and rd quantities obtained using the present dynamic

Table 1 Grid resolution for developed channel ﬂow at different Reynolds numbers. Reτ

x+

y+w

z+

395 1200 2400 4500

15 25 48 88

0.4 0.7 1.2 2.3

13 21 40 73

DDES model. For ease of comparison, the non-dimensionalized velocity proﬁles predicted by the k-equation dynamic LES model [21] and the quantities predicted by the full version of the SSTbased IDDES model [9] are presented. This LES model is chosen due to its identical formulation with the dynamic DDES model in the eddy-resolving region. At the lowest Reynolds number, the grid + with y+ w < 1 and a small x and z supports the well-resolved LES. The results of both the dynamic LES and dynamic DDES model in Fig. 6(a) show good agreement with Reichardt’s correlation [31]. As for the dynamic DDES model, the damping function takes advantage of the grid and allows the LES mode to activate in the near wall region. We still have a thin RANS layer close to the wall due to the large rd value. However, this RANS layer, which only prevails in the viscous sublayer, is much thinner than that in Yin’s model [19], and is expected to be eliminated with further increase of the grid resolution. As the Reynolds number increases to Reτ = 1200, a small discrepancy (< 5%) between the LES result and Reichardt’s correlation is seen (Fig. 6(b)), while the problem of LLM is improved thanks to the RANS layer, which is ﬁve times thicker compared with that at the lower Reynolds number. At higher Reynolds numbers, i.e. Reτ = 2400 and 4500 (Fig. 6(c) & (d)) when y+ w > 1, the results of the LES model exhibit substantial disagreement (15%) with Reichardt’s correlation, as the grid resolution does not comply with the fully resolving of the turbulent boundary layer, resulting in a gross underprediction of skin friction coeﬃcient. Accordingly, the automatic adaptation mechanism is triggered in the dynamic DDES model with the enlargement of the RANS region close to the wall. The SST model with the corresponding wall function gives a good prediction of the mean velocity proﬁle in a large portion of the boundary layer. When the model coeﬃcients in the outer region of the boundary layer are dynamically updated and ﬁlter width is redeﬁned (Eq. (13)), the LLM problem is signiﬁcantly improved. The IDDES gives similar mean velocity proﬁles at all the Reynolds numbers compared with the dynamic DDES model, but the RANS-damping mechanism works in different ways. The quantity rdt (turbulent analogue of rd used in the shielding function definition) remains signiﬁcantly small while fd is only limited by the maximum size of the grid in the vicinity of the wall. Note that only fd is discussed for the hybrid length scale since fe = 0 in most of the region. This gives rise to the grid-only dependence of the

180

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 6. U + proﬁles in developed channel ﬂow at different Reτ computed by different models.

thin RANS region at Reτ = 395 and 1200, when the RANS region is much thicker than that given by the dynamic DDES model. At Reτ = 2400 and 4500, rdt begins to take effect and the thickness of the RANS region increases in response to the local ﬂow conditions, but the RANS region keeps considerably thinner than that given by the dynamic DDES. A further comparative view between the original DDES/IDDES model [9] and the present dynamic DDES model is made clear in simulation of the channel ﬂow at Reτ = 395 (Fig. 7). The original DDES model with Cd1 = 8 exhibits a large RANS region in the computational domain, as shown in Fig. 6(a), when the grid resolution meets the LES requirements. The disagreement of the velocity proﬁle in the log-layer is apparent, which arises in the eddy-resolving region. In addition, the rd quantity experiences non-monotonicity in the range of y+ = 10 – 20, where the switching to LES begins to take place. This situation is the primary mechanism for the formation of a twisted hybrid interface. The non-monotonicity of rd also presents in the range of y+ = 15 – 30 in Fig. 6(d); however, this is not a troublesome issue, as the RANS-to-LES switching occurs at y+ = 200. When the constant Cd1 = 20 is used as recommended by Gritskevich et al. [9], the original DDES model suffers from the LES-damping because almost all the domain is occupied by the RANS region as indicated by Fig. 7(b), even the ﬂow is initialized with abundant of turbulent ﬂuctuation. The RANS-to-LES switching commences at y+ = 100 and continues in the remainder the region, where the model is neither RANS nor LES. The large “grey area” results in the disagreement of the velocity proﬁle in the whole log-layer region. Fig. 7(c) and (d) show the resolved turbulence and the total turbulent kinetic energy computed by original DDES/IDDES and the present dynamic DDES model, along with the comparison with DNS data [32]. The total turbulent kinetic energy is deﬁned by

ktot = k +

1 2 u + v 2 + w 2 2

(20)

The underprediction of the resolved turbulence (Fig. 7(c)) computed by the original DDES model with Cd1 = 8 is due to the presence of a signiﬁcant RANS region near the wall. The discrepancy in the LES region is also large. The result of the original DDES model with Cd1 = 20 is not shown here since almost all the resolved turbulence is killed by the RANS. The result of the dynamic DDES model shows a good agreement of the solved turbulence with the DNS data, as the dynamic LES model achieves almost in the entire computational domain. It is worth noting that the thin RANS layer has little effect on the turbulence resolved near the wall. Thus, the total turbulent kinetic energy computed by the dynamic DDES model agrees well with the DNS data. When the original DDES model with Cd1 = 20 is used, the performance of the RANS model is poor, as shown in Fig. 7(d), where the total turbulent kinetic energy is substantially underestimated. Shifting Cd1 = 8 results in signiﬁcant overestimation of the total turbulent kinetic energy in the LES region. The resolved turbulence and the total turbulent kinetic energy computed by the IDDES model are also plotted in Fig. 7(c) and (d) for comparison, where a slight overprediction of the streamwise velocity ﬂuctuation and the total turbulent kinetic energy is seen at the peak. Fig. 8 shows the contours of the quantity y+ which is used as local the controller for the RANS damping. For each subﬁgure, the upper and lower borders represent the wall. At the lowest Reynolds number (Fig. 8(a)), nearly all the y+ quantities are less than 1 local except in the thin layer close to the wall, where the strain due to the shear is quite large. This gives a substantial damping to the shielding function and achieves the LES mode in almost all the ﬂuid domain. y+ increases signiﬁcantly with the increasing local Reynolds number, giving rise to weakened damping on the shielding function. It is notable that y+ shows strong spatial variation, local and several troughs are found even near the wall. Hence, it would be a wise choice to use a local average (Eq. (18)) to smooth the hybrid interface and eliminate small RANS islands.

C. He et al. / Computers and Fluids 146 (2017) 174–189

181

Fig. 7. Mean and turbulent quantities in developed channel ﬂow at Reτ = 395: (a) DDES with Cd1 = 8, (b) DDES with Cd1 = 20, (c) resolved turbulence compared with DNS, (d) total turbulent kinetic energy compared with DNS.

Fig. 8. Contours of y+ in developed channel ﬂow at different Reτ computed by dynamic DDES. (For each subﬁgure, the upper and lower border represents the wall). local

3.2. Flow over periodic hills This is a typical case showing the ﬂow separation from a smooth surface. The geometry and ﬂow conditions were described by Fröhlich et al. [33]. The size of the computational domain was 9H and 4.5H in the streamwise and spanwise directions respectively, where H was the hill height at the crest. The Reynolds number based on the hill height and the bulk velocity at the crest was 10,595. The grid used had 218 × 80 × 90 nodes in the streamwise, spanwise and wall normal directions. The periodic boundary conditions were enforced along the streamwise and spanwise directions. The ﬂow was driven by a pressure gradient source term which was

adjusted to sustain the required bulk velocity at the crest. A maximum local CFL number ≈ 1.2 was maintained throughout the entire domain. Fig. 9(a) and (b) respectively show the mean and ﬂuctuation of the streamwise velocity at x/H = 0.5, 2, 4, 6 and 8 with respect to a coordinate originating at the crest. It shows excellent agreement with the experimental data measured by Breuer et al. [27]. Fig. 9(c) also shows the skin friction Cf distribution compared with the LES data from [33]. Overall agreement is reached except for the region near the crest (x/H = 0). Reddy et al. [13] claimed that this discrepancy may be attributed to the streamwise grid spacing near the crest, which can be eliminated by further reﬁnement of the grid

182

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 9. Flow over periodic hills: (a) mean velocity proﬁles, (b) velocity root-mean-square proﬁles, (c) skin-friction coeﬃcient along the bottom wall.

Fig. 10. Contours of shielding function fd in ﬂow over periodic hills computed by different models.

resolution. The original DDES model with Cd1 = 20 [9] gives similarly accurate prediction except for the skin friction near x/H = 3 where is in the recirculation zone of the separated ﬂow. Fig. 10 shows the contours of the shielding function calculated by the original DDES model and the present dynamic DDES model. The original DDES model is quite conservative, as a large portion of the RANS region occurs in the ﬂow domain (Fig. 10(a)). This kills a large portion of the turbulence near the wall and in some situations will deteriorate the ﬂow ﬁeld simulation, as the RANS model is problematic in predicting ﬂows with massive separation. However, the dynamic DDES model has a helpful damping mechanism constraining the much thinner RANS region near the wall. In ad-

Fig. 11. Isosurface of Q = 0.05 colored by the vorticity in ﬂow over periodic hills computed by different models.

dition, the spatial uniformization operator imposed on the velocity gradients considerably smooths the hybrid interface. The beneﬁt of the present dynamic DDES model can also be examined by mean of the instantaneous isosurface of the vortex structures determined using the Q-criterion (Fig. 11). As the dynamic procedure is activated, the model can capture much smaller vortices in the

C. He et al. / Computers and Fluids 146 (2017) 174–189

183

Fig. 12. Ribbed channel ﬂow computed by different models: (a) mean velocity proﬁles, (b) normal-stress component u u , (c) normal-stress component v v , (d) shear-stress component u v .

ﬂow ﬁeld using the same mesh compared with the original DDES model. It can be seen that most of the vortex structures are killed by the original DDES model due to the large RANS region near the wall (Fig. 11(a)), while the small vortex structures near the upper wall are successfully captured by the dynamic DDES model (Fig. 11(b)). 3.3. Ribbed channel ﬂow Unlike the periodic hills conﬁguration, the ﬂow in the ribbed channel is separated at the leading edge of the rib, and reattaches onto the bottom wall behind it. This is a typical case in which the velocity gradients are strengthened on the top of the rib and in the separated shear layer, while becoming signiﬁcantly small a short distance away from the wall. In this situation, the spatial uniformization operator is highly desirable for the dynamic DDES model to deﬁne a smooth hybrid interface. The geometry and ﬂow conditions were described in [34]. The dimension of the computational domain was 7.2e, 5e and 2H in the streamwise, spanwise and wall normal directions respectively, where e was the height of the square rib and H = 2.5e was the half width of the channel. The Reynolds number based on the rib height and the bulk velocity above the rib was 37,200. The grid used had a total of 2,60 0,0 0 0 nodes in the whole computational domain. The periodic boundary conditions were enforced along the streamwise and spanwise directions. The ﬂow was driven by a pressure gradient source term which was adjusted to sustain the required bulk velocity above the rib. A maximum local CFL number ≈ 2.5 was maintained throughout the entire domain. The mean streamwise velocity and the second moment computed by the original DDES model with Cd1 = 20 [9] and the current dynamic DDES model are shown in Fig. 12; the experimental data measured by Drain and Martin [28] are extracted from [34]. The data at x/e = 0, 3.68, 4.82 and 6.8 with respect to the rib center are presented. Fig. 12(a) shows that the original DDES

model overestimates the mean streamwise velocity in the recirculation bubble and a slightly underestimation in the separated shear layer. However, the mean streamwise velocity determined by the dynamic DDES model agrees well with the measurement data, which is attributed to the resolved small eddies using the dynamically computed model coeﬃcients (Ck & Ce ). In Fig. 12(b) – (d), the present dynamic DDES model demonstrates relatively better predictions of the Reynolds stresses. The original DDES model’s larger over-prediction of the ﬂuctuation is seen in the separation and reattachment region, further demonstrating the advantage of the present dynamic DDES model in simulation of the strong shear and separated ﬂows. Fig. 13 shows the spatial distribution of the eddy (subgrid) viscosity, in which the hybrid interface is plotted for ease of understanding. The RANS region in Fig. 13(b) is thinner and much smoother than that in Fig. 13(a) due to the damping function and spatial uniformization operator in fd . Both models predict a smaller eddy viscosity in the freestream region, as the turbulence is not as strong as in the shear layer. The signiﬁcant difference is that the eddy viscosity computed by the present dynamic DDES model is much smaller than with the original DDES model. In addition, the spatial distribution of the eddy viscosity in Fig. 13(b) is more discrete resulting from the wide variation of the coeﬃcients Ck in space as demonstrated in Fig. 14(a). Recall that the coeﬃcient Ce is associated with the dissipation of kinetic energy. Accordingly, a large value of Ce (with maximum value larger than 2.0) is found in the freestream region, as shown in Fig. 14(b), where the turbulence is relatively weak compared to the shear layer. Near the rib and bottom wall, there exists another region of large Ce , which however has no effect on the ﬂow, as RANS prevails in this region. The inﬂuence of spatial uniformization on the velocity gradients is highlighted in Fig. 15, where the contour of the normalized velocity gradients computed by Yin’s model [19] without the spatial uniformization is presented. It is obvious that the velocity gradients are signiﬁcantly large near the wall and in the separated shear

184

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 13. Turbulent eddy viscosity in ribbed channel ﬂow computed by different models. (White curve represents the RANS/LES interface).

Fig. 15. Velocity gradient and RANS/LES interface in ribbed channel ﬂow computed by Yin’s dynamic DDES model.

is formed close to the surface and the WMLES mode is achieved. In the simulation of turbulent heat transfer in LES region, a closure of the turbulent Prandtl number Prt is employed, which is derived based on the assumption that the turbulent Prandtl number is unique at different ﬁlter levels [30]. Deﬁning the turbulent kinetic energy at the test level K as

K=

1 ˆ U ˙U ˙ − Ui˙Ui˙ + k 2 i i

(21)

the turbulent Prandtl number in the LES region is determined by

P rt =

Qj · Qj Pj · Q j

where Fig. 14. Model coeﬃcients in ribbed channel ﬂow computed by dynamic DDES: (a) Ck , (b) Ce .

layer. At a short distance away from the wall, however, several regions with very small velocity gradients are found, which is the main mechanism for the increase of rd and the resultant decrease of fd . As a consequence, a twisting coastline-like hybrid interface is formed, which results in inaccurate prediction (not shown here) and numerical instability in the present simulation. 3.4. Impinging jet with heat transfer In this simulation, the ﬂow keeps steady in the pipe, and switches from the RANS mode to the LES mode immediately out of the pipe. Where the jet impinges on the wall, a thin RANS region

(22)

√ √ ˆ ∂T ˆ ∂ T − Ck Q j = Ck K K ∂xj ∂xj Pj = Tˆ · Uˆ j − T · Uj

(23-a) (23-b)

In the DDES formulation, the shielding function is also used to hybridize the dynamically computed turbulent Prandtl number in the LES region and the constant one in the RANS region, giving

P rt = fd

Qj · Qj + (1 − fd )P rt0 Pj · Q j

(24)

Here, Prt0 = 0.85 is the constant turbulent Prandtl number used in the RANS region. The test case is based on the experimental conﬁguration investigated by Cooper et al. [29]. The Reynolds number based on the bulk velocity in the pipe and the pipe diameter was 23,0 0 0. The

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 16. Schematic representation of the computed domain for impinging jet.

oriﬁce-to-plate distance was H/D = 2. The computational domain was shown in Fig. 16. A circular jet issuing from a fully developed pipe ﬂow, entered a cylindrical ﬂuid domain with a diameter of 12D and a height of 3D, and impinged on the target plate located 2D from the pipe exit. This computational domain was deﬁned based on the study of Hadžiabdic´ and Hanjalic´ [35]. Other than imposing the inﬂow on the pipe exit directly, a section of pipe with length 0.5D was simulated along with the ﬂow domain to eliminate the effect of interpolation error on the ﬂow ﬁeld. Meanwhile, a small pipe thickness was considered for convenience of grid generation. The grid used here consisted of 117 × 135 × 300 points in height, radial and azimuthal direction respectively, and results in a total node number of around 4.8 million. To test the DDES model, a fully developed velocity, turbulent kinetic energy and turbulent eddy frequency distribution was obtained from a previous SST RANS simulation in a pipe, and was imposed on the inlet boundary without any turbulence content. A noslip condition for velocity and a constant heat ﬂux as the experimental condition was imposed on the target surface. As for the top and cylindrical surface, the deﬁnition of the boundary condition was not straightforward for such an open domain. Hadžiabdic´ and Hanjalic´ [35] noticed that there was no signiﬁcant difference in the results within the jet itself for different top boundary conditions. Hence, a free-slip condition for velocity and a constant pressure condition were employed in current computations. For the outﬂow boundary, the convective boundary conditions deﬁned by the hyperbolic convective equation were applied, which read

∂ψ + ∇ · (ρ U ψ ) = 0 ∂t

(25)

where U was the velocity calculated from the ﬂux at the previous time step. Here, ψ could be the velocity, temperature, turbulence kinetic energy and turbulent eddy frequency. The simulations presented here were performed by solving the N-S equations and the passive-scale transport equation. Both the original DDES model with Cd1 = 20 and the dynamic Smagorinsky subgridscale LES model with local average using the test ﬁlter were used for comparison. In order to avoid the unphysical ﬂuctuation results from the unboundness of the central differencing scheme in complex ﬂows, the hybrid central-upwind differencing along with explicit gradient correction were used for the surface interpolation of the convection term, which casts in the form

Uip = σ Ucentral + (1 − σ )Uupwind + (1 − σ )∇ U · r

(26)

where Uip was the velocity at the surface center, Ucentral and Uupwind were surface velocities respectively predicted by the central and upwind schemes. r was the vector from the upwind node to the surface center, and σ was the hybrid coeﬃcient (set to 0.75 in the present simulation). The time step t was chosen to ensure that the maximum local CFL number ≈ 2.8.

185

The time-averaged velocity modulus determined by the original DDES model, the dynamic LES model, the current dynamic DDES model and the experimental measurement [29] are compared in Fig. 17. The velocities are normalized using U0 , which denotes the bulk velocity in the pipe. Note that the node number of the grid used in the present simulation is much less than that used by Hadžiabdic´ and Hanjalic´ [35] (which was as much as 10 million), since the current intention is to demonstrate the dynamic DDES model on a relative coarse grid, as is most practical for industry applications. The dynamic Smagorinky LES model and the present dynamic DDES model provide similar results at the radial stations r/D = 0.5, 1.0 and 1.5 as shown in Fig. 17(a) – (c), where the discrepancy with the experimental data is not as pronounced as for the original DDES model. Although all the models overpredict the velocity in the stagnation region in range y/D = 0 – 0.25 as shown in Fig. 17(a), which results from the steady inlet condition used in this simulation, the agreement with the experimental results improves for the dynamic LES model and the present dynamic DDES model at larger radial stations (Fig. 17(b) and (c)), where the original DDES model shows relatively large deviations. Until the radial stations r/D = 2.0, 2.5 and 3.0 as shown in Fig. 17(d), (e) and (f) respectively, the situation improves greatly for the original DDES model. In addition, both the original DDES model and the present dynamic DDES model tend to more accurately predict the velocity near the wall at larger radial stations, while the dynamic LES model yields the overpredicted value near y/D = 0.05 at r/D = 2.5 and 3.0. This can be explained by the inadequacy of current grid resolution for LES, whereas the DDES models can take full advantage of the WMLES mode. Note that the overpredicted velocity at y/D > 0.2 is attributed to the grid resolution in the axial and azimuthal directions as declared by Hadžiabdic´ and Hanjalic´ [35]. A good DES model is expected to switch from the RANS mode to the eddy-resolving mode under the actuation of the ﬂow instability. Although the speed-up strategies for the RANS-to-LES transition introduced by Mockett et al. [16] and Shur and Spalart [17] are not employed in the present formulation, the ﬂow predicted by the dynamic DDES model has seen rapid transition to LES mode immediately out of the nozzle. As shown in Fig. 18, the turbulence intensity predicted by the dynamic DDES model is closer to the experimental data in the whole range. The dynamic LES model tends to have similar performance, but exhibits a slight overprediction near the wall, which is particularly obvious in Fig. 18(b), (c), (e) and (f). Neither of the two models capture the ﬂuctuation near the wall at r/D = 0.5, which is attributed to the steady inlet condition. As for the original DDES model, the ﬂuctuation is substantially underestimated at r/D = 0.5, 1.0 and 1.5 (Fig. 18(a) – (c)). This means that the switching from RANS to LES exhibited by the original DDES model is rather conservative. Even at a large radial station (Fig. 18(d) – (f)), the ﬂuctuation intensities predicted by the original DDES model show pronounced disagreement with the experimental data. For the time-averaged Nusselt number along the bottom surface (Fig. 19), the original DDES model has a substantial overshoot in the stagnation region (r/D < 2.0), while the results predicted by the dynamic LES model and the dynamic DDES model have relatively better agreement with the experimental measurement. All the models show underprediction at larger radial stations in the wall-jet region (r/D > 2.0), where the computed Nusselt number shows a similar monotonic decrease in the radial direction as in the measurement. Hadžiabdic´ and Hanjalic´ [35] identiﬁed the same problem, attributing it to insuﬃcient grid resolution in the walljet region (note that the total node number of their grid is up to 10 million). However, the dynamic DDES model provides a globally higher Nusselt number in the whole range due to the thin RANS region near the wall, which increases the eddy viscosity and hence the eddy diffusivity.

186

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 17. Mean velocity proﬁles in impinging jet at different locations: (a) r/D = 0.5, (b) r/D = 1.0, (c) r/D = 1.5, (d) r/D = 2.0, (e) r/D = 2.5, (f) r/D = 3.0.

Fig. 18. Velocity ﬂuctuation proﬁles in impinging jet at different locations: (a) r/D = 0.5, (b) r/D = 1.0, (c) r/D = 1.5, (d) r/D = 2.0, (e) r/D = 2.5, (f) r/D = 3.0.

C. He et al. / Computers and Fluids 146 (2017) 174–189

187

Fig. 19. Time-averaged Nusselt number in impinging jet computed by different models.

Fig. 21. Isosurface of Q-criterion = 5 ×105 colored by the axial velocity in impinging jet computed by different models. Fig. 20. Turbulent eddy viscosity in impinging jet computed by different models.

A further comparative view of different models can be inferred from contours of the eddy (subgrid) viscosity in Fig. 20. For ease of comparison, the same legend is used for all the models. It is obvious that the eddy viscosity in the original DDES model (Fig. 20(a)) is globally elevated to far beyond saturation, along with signiﬁcantly large eddy viscosity clustered within the jet, near the bottom surface and the outlet. Such substantial elevation in the eddy viscosity would kill a considerable portion of the resolved turbulence and suppress the switching from RANS mode to LES mode. As for the dynamic LES model (Fig. 20(b)) and the present dynamic DDES model (Fig. 20(c)), similar patterns are found except in the region in the pipe and near the bottom wall, where RANS mode prevails. This indicates that the dynamic DDES model can shield the steady internal pipe ﬂow, and switch rapidly to the dynamic LES mode out of the pipe by attenuating the modeled turbulent eddy viscosity. In the wall-jet region, the WMLES is activated with only a thin RANS region in the vicinity of the wall. Fig. 21 clearly demonstrates the ﬂow structures in terms of the isosurface of Q-criterion colored by the axial velocity component. The ring-like vortical structures near the wall shown in Fig. 21(a), which are almost continuous in the azimuthal direction, serve as an indication of high dissipation by the original DDES model. They are originally generated in the shear layer of the jet, but do not break down when impinging on the wall and propagating downstream. In contrast, abundant turbulence content is resolved by the dynamic LES model (Fig. 21(b)) and the present dynamic DDES model (Fig. 21(c)). A large number of ring-like vortical structures are generated in the shear layer of the jet, and then break down,

turning into smaller vortical structures immediately after the impingement. Fig. 22 shows contours of the Nusselt numbers on the target surface, which are strongly associated with the ﬂow structures. The original DDES model results in a ring pattern of Nusselt numbers, along with very large values in the stagnation region. As discussed previously, the Nusselt number in the stagnation region would be considerably overpredicted. Although a thin RANS region prevails near the wall, the present dynamic DDES model yields a Nusselt number pattern similar to the dynamic LES model. This means that the present dynamic DDES model is comparable with the dynamic LES model in the eddy-resolving region. 4. Concluding remarks A new dynamic DDES model is proposed based on the k-ω SST model and the well-established dynamic k-equation LES subgridscale model. The same form of k-equation is maintained in both the RANS and LES regions. The model coeﬃcients Ck and Ce (the two coeﬃcients correspond to CDES in the original DDES formulation) are computed dynamically by taking into account the local ﬂow ﬁeld at the grid and test ﬁlter levels. To obtain a smooth and monotonous hybrid interface, a spatial uniformization operator is imposed on the velocity gradient, which is used in the deﬁnition of the shielding function fd . Accordingly, the constant Cd1 in fd is recalibrated using the ﬂow over a ﬂat plate. A quantity y+ is local deﬁned analogously to the wall unit y+ , which serves as a criterion whether the local grid can support the LES. Thus, a damping function ϕ d deﬁned based on y+ is employed in the fd funclocal tion which damps the RANS region dynamically to achieve WMLES

188

C. He et al. / Computers and Fluids 146 (2017) 174–189

Fig. 22. Nusselt number distribution in impinging jet computed by different models.

mode. Four test cases, developed channel ﬂow, ﬂow over periodic hills, ribbed channel ﬂow and impinging jet with heat transfer, are used to evaluate the model. The test in simulation of ﬂat-plate boundary layer ﬂow shows the constant Cd1 = 80 provides enough shielding for the boundary layer. The damping procedure is not activated, as the y+ local remains signiﬁcantly large in the boundary layer. This leads to a similar prediction of wall friction coeﬃcient Cf as the k-ω SST model. When the spatial uniformization operator is applied, the decrease in fd outside the boundary layer, which occurs in the original DDES model, is eliminated. When the model is applied in the developed channel ﬂow, the thickness of the RANS region is dynamically adjusted depending on the grid resolution and local ﬂow conditions. The problem of LLM is improved and the model convincingly demonstrates considerable improvement over the dynamic LES model and original DDES model. When applied in separated ﬂows, the present dynamic DDES model exhibits better performance than the original DDES model and the dynamic LES model. The tests of the model on ﬂow over periodic hills and the ribbed channel show that a thinner RANS region is obtained and the hybrid interface is much smoother, and much smaller vortical structures are captured. In the simulation of ribbed channel ﬂow, with the dynamical computation of model coeﬃcients based on the local grid resolution and ﬂow conditions, the eddy viscosity is much smaller than the original DDES model, accurately predicting the ﬂow in both the shear layer and the recirculation bubble compared with that obtained by the original DDES model. In addition, the impinging jet ﬂow shows the present dynamic DDES model exhibits rapid switching from the RANS mode to LES mode when the inﬂow does not include any turbulent content, while the original DDES model kills most of the turbulence in the jet region and exhibits excessive dissipative behavior in the whole computational domain. Subsequently, implementation of the turbulent Prandtl number model within the DDES framework gives a better prediction of the turbulent wall heat transfer in comparison with the dynamic LES model. Future work will be directed to improving the present dynamic DDES model for more complex ﬂow conditions. The spatial uniformization operator imposed on the velocity gradients, which employs the domain average, plays a signiﬁcant role in obtaining the monotonous shielding function and smoothing the RANS-LES interface. Of particular interest is the generalization with different domain size in external ﬂows, which needs further effort. Extensive comparisons with the traditional DDES/IDDES and the dynamic LES models, and the alternative ﬁlter width deﬁnition which takes the grid isotropy into consideration, will also be considered in future studies. Conﬂict of interest The authors declare that they have no conﬂict of interest.

Acknowledgement The authors gratefully acknowledge ﬁnancial support from the National Natural Science Foundation of China (No.11372189) and the Chinese ‘111 Program of Introducing Talents of Discipline to Universities.

References [1] Piomelli U, Balaras E. Wall-layer models for large-eddy simulations. Annu Rev Fluid Mech 2002;34:349–74. [2] Foroutan H, Yavuzkurt S. A partially-averaged Navier–Stokes model for the simulation of turbulent swirling ﬂow with vortex breakdown. Int J Heat Fluid Flow 2014;50:402–16. [3] Kniesner B, Šari S, Mehdizadeh A, Jakirli S, Hanjali K, Tropea C, et al. Wall treatment in LES by RANS models: method development and application to aerodynamic ﬂows and swirl combustors. 2007. [4] Shur M, Spalart P, Strelets M, Travin A. Detached-eddy simulation of an airfoil at high angle of attack. Eng Turbul Model Exp 1999;4:669–78. [5] Spalart P, Jou W, Strelets M, Allmaras S. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Adv DNS/LES 1997;1:4–8. [6] Strelets M. Detached eddy simulation of massively separated ﬂows. AIAA, Aerospace Sciences Meeting and Exhibit, 39 th, Reno, NV; 2001. [7] Menter F, Kuntz M, Langtry R. Ten years of industrial experience with the SST turbulence model. Turbul Heat Mass Transf 2003;4. [8] Spalart PR, Deck S, Shur M, Squires K, Strelets MK, Travin A. A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theor Comput Fluid Dyn 2006;20:181–95. [9] Gritskevich MS, Garbaruk AV, Schütze J, Menter FR. Development of DDES and IDDES formulations for the k-ω shear stress transport model. Flow Turbul Combust 2012;88:431–49. [10] Nikitin N, Nicoud F, Wasistho B, Squires K, Spalart P. An approach to wall modeling in large-eddy simulations. Phys Fluids 20 0 0;12:1629–32 (1994-present). [11] Piomelli U, Balaras E, Pasinato H, Squires KD, Spalart PR. The inner–outer layer interface in large-eddy simulations with wall-layer models. Int J Heat Fluid Flow 2003;24:538–50. [12] Shur ML, Spalart PR, Strelets MK, Travin AK. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. Int J Heat Fluid Flow 2008;29:1638–49. [13] Reddy K, Ryon J, Durbin P. A DDES model with a Smagorinsky-type eddy viscosity formulation and log-layer mismatch correction. Int J Heat Fluid Flow 2014;50:103–13. [14] Spalart PR. Detached-Eddy Simulation. Fluid Mechanics 2009;41:203–29. [15] Chauvet N, Deck S, Jacquin L. Zonal detached eddy simulation of a controlled propulsive jet. AIAA J 2007;45:2458–73. [16] Mockett C, Fuchs M, Garbaruk A, Shur M, Spalart P, Strelets M, et al. Two non– zonal approaches to accelerate RANS to LES transition of free shear layers in DES. Notes on. Numer Fluid Mech Multidiscipl Des 2015;130:187–201. [17] Shur ML, Spalart PR. An enhanced version of DES with rapid transition from RANS to LES in separated ﬂows. Flow Turbul Combust 2015;95:1–29. [18] Mockett CA. Comprehensive study of detached eddy Simulation: univerlagtuberlin; 2009. [19] Yin Z, Reddy K, Durbin PA. On the dynamic computation of the model constant in delayed detached eddy simulation. Phys Fluids 2015;27(2):4–8. [20] Menter FR. Improved two-equation k-omega turbulence models for aerodynamic ﬂows; 1992. [21] Kim W-W, Menon S. A new dynamic one-equation subgrid-scale model for large eddy simulations. AIAA, aerospace sciences meeting and exhibit, 33 rd, Reno, NV; 1995. [22] Germano M, Piomelli U, Moin P, Cabot WH. A dynamic subgrid-scale eddy viscosity model. Phys Fluids A Fluid Dyn 1991;3:1760–5. [23] Meneveau C, Lund TS, Cabot WH. A Lagrangian dynamic subgrid-scale model of turbulence. J Fluid Mech 1996;319:353–85.

C. He et al. / Computers and Fluids 146 (2017) 174–189 [24] Comte-Bellot G, Corrsin S. Simple Eulerian time correlation of full-and narrow-band velocity signals in grid-generated,‘isotropic’turbulence. J Fluid Mech 1971;48:273–337. [25] Menon S, Kim W-W, Menon S. High Reynolds number ﬂow simulations using the localized dynamic subgrid-scale model. AIAA 960425, 34th AIAA aerospace sciences meeting: Citeseer; 1996. [26] Rapp C, Manhart M. Flow over periodic hills: an experimental study. Exp Fluids 2011;51:247–69. [27] Breuer M, Peller N, Rapp C, Manhart M. Flow over periodic hills–numerical and experimental study in a wide range of Reynolds numbers. Computers & Fluids 2009;38:433–57. [28] Drain L, Martin S. Two-component velocity measurements of turbulent ﬂow in a ribbed-wall ﬂow channel. In: International conference on laser anemometry—advanced and application; 1985. p. 99–112. [29] Cooper D, Jackson D, Launder BE, Liao G. Impinging jet studies for turbulence model assessment—I. Flow-ﬁeld experiments. Int J Heat Mass Transf 1993;36:2675–84.

189

[30] Moin P, Squires K, Cabot W, Lee S. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Phys Fluids A 1991;3:2746–57. [31] Reichardt H. Vollständige Darstellung der turbulenten Geschwindigkeitsverteilung in glatten Leitungen. ZAMM-J Appl Math Mech/Zeitschrift für Angewandte Mathematik und Mechanik 1951;31:208–19. [32] Moser RD, Kim J, Mansour NN. Direct numerical simulation of turbulent channel ﬂow up to Re= 590. Phys Fluids 1999;11:943–5. [33] Fröhlich J, Mellen CP, Rodi W, Temmerman L, Leschziner MA. Highly resolved large-eddy simulation of separated ﬂow in a channel with streamwise periodic constrictions. J Fluid Mech 2005;526:19–66. [34] Weihing P, Younis B, Weigand B. Heat transfer enhancement in a ribbed channel: development of turbulence closures. Int J Heat Mass Transfer 2014;76:509–22. [35] Hadžiabdic´ M, Hanjalic´ K. Vortical structures and heat transfer in a round impinging jet. J Fluid Mech 2008;596:221–60.