A dynamic delayed detached-eddy simulation model for turbulent flows

A dynamic delayed detached-eddy simulation model for turbulent flows

Computers and Fluids 146 (2017) 174–189 Contents lists available at ScienceDirect Computers and Fluids journal homepage: www.elsevier.com/locate/com...

5MB Sizes 1 Downloads 9 Views

Computers and Fluids 146 (2017) 174–189

Contents lists available at ScienceDirect

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

A dynamic delayed detached-eddy simulation model for turbulent flows 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 flow 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 coefficient CDES in traditional DES formulations, the present model employs two coefficients Ck and Ce , which are computed dynamically by taking into account the spatial and temporal variations of the flow field at the grid and test filter levels. A modification 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 flow 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 flow shows the log-layer mismatch (LLM) problem is significantly 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 flows, with the model coefficients 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 flow 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 flow instabilities when the inflow 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 flows and the high expense of LES models with high-Reynolds number wallbounded flows. 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.compfluid.2017.01.018 0045-7930/© 2017 Elsevier Ltd. All rights reserved.

separated flow 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 flow quantities of the turbulent boundary layer, as the grid resolution is insufficient 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 refined 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 flow configurations with enough turbulence in the inflow, and the rapid RANS-to-LES transition when a steady inlet boundary condition is imposed. However, simulation of the channel flow [10,11] using the DES model identified 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 flows where WMLES is highly desirable. In this regard, the improved DDES (IDDES) was developed with a series of blending functions and a modified filter width [9,12]. This model improved LLM and ensures an automatic selection of DDES or WMLES mode depending on the turbulence content in the flow. Motivated by the exploration of an alternate and simpler DES formulation which could overcome both the MSD and LLM, Reddy [13] recently modified the DDES formulation, in which the eddy viscosity was similar to the Smagorinsky formula far from the wall. With a hybrid filter width  defined, the problem of LLM was significantly improved. When a steady inlet boundary condition is used, DES relies on the natural instability of the flow 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 filter 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] modified the filter width definition 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 flows 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 flow configurations. In DDES models [6,8,9], the model coefficient 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 flows. In this regard, a dynamic version of DDES was proposed [19] by dynamically determining the model coefficient CDES based on the local flow, demonstrating relatively accurate prediction of the fully developed stationary channel flow and rotating channel flow in comparison with traditional DDES models (constant CDES ). However, in Yin’s model, the shielding function fd defined by Spalart et al. [8] was used to hybridize RANS and LES based on spatial locations, which would introduce inaccurate prediction of complex flows. Recall that even in the constantcoefficient DDES models [9], significant decrease of the shielding function was found in the free-stream region of the flat-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 flows when the model coefficient 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 flow 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 flow 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 coefficients 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 constantcoefficient DDES model [9] is modified 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 flow conditions. Accordingly, the constant in the definition of the shielding function is recalibrated using the flat-plate boundary layer. The open source code OpenFOAM (http://openfoam.org) has been used for all the present simulations. Gaussian finite 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 flow, flow over periodic hills, ribbed channel flow 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 define 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 defined 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 define the LES length scale

LLES = Ck 

(5)

the subgrid eddy viscosity becomes

 μsgs = LLES k

(6)

In this model, the model coefficients 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 defined length scales, the dynamic DDES model is

In these equations, U is the Reynolds-averaged velocity tensor or the filtered velocity tensor at the grid level in corresponding models. “ ˆ ” denotes the filter operator at the test level. Note that in the Smagorinsky model, the dynamically-determined model coefficient is averaged either in the homogeneous direction [22] or along the flow 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 field, 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 filter 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 define 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 coefficients Ck and Ce is achieved which accurately predicts the turbulence in the detached flows. In the original DDES formulation, the filter width is defined 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 filter width is defined as

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

(13)

In the eddy-resolving region where fd = 1, this expression gives a filter 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 definition of the filter 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 fluctuation but with the variation tendency and average value remaining almost the same. The modified shielding function is defined 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 defined 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 flat-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 significant problem for most of the external flows as the flow 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 field 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 flow over a flat 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 insufficient 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 modified 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 flow regions. Alternatively, damping of the shielding function is expected when the grid resolution is fine enough for LES, especially in flow 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 flat 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 refinement. r gives the damping coefficient between 0 and 1 based on the value of y+ . Generally, y+ relocal local mains significantly 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 flows. Hence, the local average is used in the near-wall region. For convenience, we define



f = tanh Cdf rd and

ϕd =

Fig. 4. Contours of shielding function fd in flat 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 defined 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 fine 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 defined as the maximum value of local the grid aspect ratio (defined 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 defined 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 specified 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 definition. 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 definition of ϕ d is applicable to various flow configurations. Finally, the shielding function in the present dynamic DDES model is defined 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 flows without inflow turbulence content, y+ is significantly large local and hence gives no damping on the shielding function until very fine grid resolution is used to activate the WMLES mode. Fig. 5 shows the distribution of skin friction on the flat 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 coefficient Cf over the flat 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 flow over the flat 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 flow, flow over periodic hills [26,27], ribbed channel flow [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 flow 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 flow This test is the most important for evaluation of the new model’s capability when applied to non-separated turbulent boundary layer flows. 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 flow 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 profiles along with the fd and rd quantities obtained using the present dynamic

Table 1 Grid resolution for developed channel flow 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 profiles 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 five 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 coefficient. 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 profile in a large portion of the boundary layer. When the model coefficients in the outer region of the boundary layer are dynamically updated and filter width is redefined (Eq. (13)), the LLM problem is significantly improved. The IDDES gives similar mean velocity profiles 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 significantly 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 + profiles in developed channel flow 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 flow 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 flow 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 profile 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 flow is initialized with abundant of turbulent fluctuation. 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 profile 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 defined 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 significant 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 significant 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 fluctuation 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 subfigure, 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 fluid domain. y+ increases significantly 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 flow 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 flow at different Reτ computed by dynamic DDES. (For each subfigure, the upper and lower border represents the wall). local

3.2. Flow over periodic hills This is a typical case showing the flow separation from a smooth surface. The geometry and flow 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 flow 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 fluctuation 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 refinement of the grid

182

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

Fig. 9. Flow over periodic hills: (a) mean velocity profiles, (b) velocity root-mean-square profiles, (c) skin-friction coefficient along the bottom wall.

Fig. 10. Contours of shielding function fd in flow 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 flow. 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 flow domain (Fig. 10(a)). This kills a large portion of the turbulence near the wall and in some situations will deteriorate the flow field simulation, as the RANS model is problematic in predicting flows 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 flow over periodic hills computed by different models.

dition, the spatial uniformization operator imposed on the velocity gradients considerably smooths the hybrid interface. The benefit 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 flow computed by different models: (a) mean velocity profiles, (b) normal-stress component u u , (c) normal-stress component v v , (d) shear-stress component u v .

flow field 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 flow Unlike the periodic hills configuration, the flow 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 significantly small a short distance away from the wall. In this situation, the spatial uniformization operator is highly desirable for the dynamic DDES model to define a smooth hybrid interface. The geometry and flow 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 flow 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 coefficients (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 fluctuation 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 flows. 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 significant 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 coefficients Ck in space as demonstrated in Fig. 14(a). Recall that the coefficient 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 flow, as RANS prevails in this region. The influence 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 significantly 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 flow computed by different models. (White curve represents the RANS/LES interface).

Fig. 15. Velocity gradient and RANS/LES interface in ribbed channel flow 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 filter levels [30]. Defining 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 coefficients in ribbed channel flow 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 flow 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 configuration 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.

orifice-to-plate distance was H/D = 2. The computational domain was shown in Fig. 16. A circular jet issuing from a fully developed pipe flow, entered a cylindrical fluid 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 defined based on the study of Hadžiabdic´ and Hanjalic´ [35]. Other than imposing the inflow on the pipe exit directly, a section of pipe with length 0.5D was simulated along with the flow domain to eliminate the effect of interpolation error on the flow field. 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 flux as the experimental condition was imposed on the target surface. As for the top and cylindrical surface, the definition of the boundary condition was not straightforward for such an open domain. Hadžiabdic´ and Hanjalic´ [35] noticed that there was no significant 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 outflow boundary, the convective boundary conditions defined by the hyperbolic convective equation were applied, which read

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

(25)

where U was the velocity calculated from the flux 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 filter were used for comparison. In order to avoid the unphysical fluctuation results from the unboundness of the central differencing scheme in complex flows, 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 coefficient (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 flow 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 flow 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 fluctuation near the wall at r/D = 0.5, which is attributed to the steady inlet condition. As for the original DDES model, the fluctuation 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 fluctuation 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] identified the same problem, attributing it to insufficient 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 profiles 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 fluctuation profiles 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 significantly 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 flow, 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 flow 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 flow 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 coefficients Ck and Ce (the two coefficients correspond to CDES in the original DDES formulation) are computed dynamically by taking into account the local flow field at the grid and test filter levels. To obtain a smooth and monotonous hybrid interface, a spatial uniformization operator is imposed on the velocity gradient, which is used in the definition of the shielding function fd . Accordingly, the constant Cd1 in fd is recalibrated using the flow over a flat plate. A quantity y+ is local defined analogously to the wall unit y+ , which serves as a criterion whether the local grid can support the LES. Thus, a damping function ϕ d defined 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 flow, flow over periodic hills, ribbed channel flow and impinging jet with heat transfer, are used to evaluate the model. The test in simulation of flat-plate boundary layer flow shows the constant Cd1 = 80 provides enough shielding for the boundary layer. The damping procedure is not activated, as the y+ local remains significantly large in the boundary layer. This leads to a similar prediction of wall friction coefficient 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 flow, the thickness of the RANS region is dynamically adjusted depending on the grid resolution and local flow 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 flows, the present dynamic DDES model exhibits better performance than the original DDES model and the dynamic LES model. The tests of the model on flow 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 flow, with the dynamical computation of model coefficients based on the local grid resolution and flow conditions, the eddy viscosity is much smaller than the original DDES model, accurately predicting the flow in both the shear layer and the recirculation bubble compared with that obtained by the original DDES model. In addition, the impinging jet flow shows the present dynamic DDES model exhibits rapid switching from the RANS mode to LES mode when the inflow 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 flow conditions. The spatial uniformization operator imposed on the velocity gradients, which employs the domain average, plays a significant 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 flows, which needs further effort. Extensive comparisons with the traditional DDES/IDDES and the dynamic LES models, and the alternative filter width definition which takes the grid isotropy into consideration, will also be considered in future studies. Conflict of interest The authors declare that they have no conflict of interest.

Acknowledgement The authors gratefully acknowledge financial 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 flow 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 flows 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 flows. 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 flows. 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 flows; 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 flow 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 flow in a ribbed-wall flow 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-field 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 flow 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 flow 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.