RANS of turbulent flow separation in a 3-D diffuser

International Journal of Heat and Fluid Flow 31 (2010) 820–832

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Numerical and physical aspects in LES and hybrid LES/RANS of turbulent flow separation in a 3-D diffuser S. Jakirlic´ a,b,d,*, G. Kadavelil a,b, M. Kornhaas a,c, M. Schäfer a,c, D.C. Sternel a,c, C. Tropea a,b,d a

Department of Mechanical Engineering, Technische Universität Darmstadt, Germany Institute of Fluid Mechanics and Aerodynamics, Petersenstr. 30, D-64287 Darmstadt, Germany c Institute of Numerical Methods in Mechanical Engineering, Dolivostr. 15, D-64293 Darmstadt, Germany d Center of Smart Interfaces, Petersenstr. 32, D-64287 Darmstadt, Germany b

a r t i c l e

i n f o

Article history: Received 20 December 2009 Received in revised form 14 April 2010 Accepted 12 May 2010 Available online 6 July 2010 Keywords: 3-D separation Wall-bounded flow LES Hybrid LES/RANS

a b s t r a c t An incompressible fully-developed duct flow expanding into a diffuser, whose upper wall and one side wall are appropriately deflected (with the expansion angles of 11.3° and 2.56° respectively), and for which reference experimental and DNS databases were provided by Cherry et al. (2008, 2009) and Ohlsson et al. (2009, 2010), was studied computationally by using a zonal hybrid LES/RANS (HLR) method, proposed recently by Kniesner (2008) and Jakirlic´ et al. (2009). In addition a complementary Large-Eddy Simulation (LES) method has been applied. The flow Reynolds number based on the height of the inlet channel is Reh = 10,000. The primary objective of the present investigation was the comparative assessment of the computational models in this flow configuration characterized by a complex 3-D flow separation being the consequence of an adverse-pressure gradient evoked by the duct expansion. The focus of the investigation was on the capability of different modelling approaches to accurately capture the size and shape of the 3-D flow separation pattern and associated mean flow and turbulence features. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction Configurations involving 3-D boundary-layer separation are among the most frequently encountered flow geometries in practice. Accordingly, the methods simulating them have to be appropriately validated using detailed and reliable reference databases. However, the largest majority of the experimental benchmarks being used for validating computational methods and turbulence models relates to internal, two-dimensional flow configurations, as e.g., flow in a 2-D diffuser (e.g., Obi et al., 1993), flow over a backward-facing step and a forward-facing step, or flow over fences, ribs, 2-D hills and 2-D humps mounted on the bottom wall of a plane channel. In these examples it is assumed that the influence of the side walls (according to Bradshaw and Wong (1972), the minimum aspect ratio – representing the ratio of the channel height to channel width – should be 1:10 in order to eliminate the influence of the side walls) is not felt at the channel midplane. Consequently, within a computational framework, the spanwise direction can be regarded as a homogeneous one, the fact enabling the application of the periodic boundary conditions (even 2-D computations when using the RANS approach). By doing so, the 3-D nature of the flow is completely missed: strong secondary mo* Corresponding author at: Institute of Fluid Mechanics and Aerodynamics, Petersenstr. 30, D-64287 Darmstadt, Germany. E-mail address: [email protected] (S. Jakirlic´). 0142-727X/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.ijheatfluidflow.2010.05.004

tion across the inlet section of the channel induced by the Reynolds stress anisotropy – which is, as generally known, beyond the reach of the eddy-viscosity RANS model group, complex 3-D separation patterns spreading over duct corners (corner separation and corner reattachment), etc. These circumstances were the prime motivation for the recent experimental study of the flow in a three-dimensional diffuser conducted by Cherry et al. (2008, 2009). Such a diffuser configuration has also a high practical relevance. It mimics a diffuser situated between a compressor and the combustor chamber in a jet engine. Its task is to decelerate the flow discharging from compressor over a very short distance to the velocity field of the combustor section. Typically a uniform inlet profile over the diffuser outlet is desirable. Such a flow situation is associated by a strong pressure increase. Cherry et al. provided a detailed reference database comprising the pressure distribution along the bottom non-deflected wall, three-component mean velocity field and the streamwise Reynolds stress component field within the entire diffuser section. Recently Ohlsson et al. (2009, 2010) have performed a complementary Direct Numerical Simulation using a massively parallel high-order spectral element code. The 3-D diffuser was meshed by approximately 172 million grid points. In addition to the mean velocity field, all six Reynolds stress components were evaluated, as well as the surface pressure distribution along the bottom wall. It should be noted that in the experimental investigations two three-dimensional

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

diffusers with the same fully-developed channel inlet but slightly different expansion geometries (the upper-wall expansion angle is reduced from 11.3° to 9°; the side-wall expansion angle is increased from 2.56° to 4°) were investigated, Cherry et al. (2008). Both diffuser flows are characterized by a three-dimensional boundary-layer separation, but the size and shape of the separation bubble exhibit a high degree of geometric sensitivity to the dimensions of the diffuser. The first diffuser configuration (Fig. 1), being also the flow geometry considered in the present work, has served as a test case of the 13th and 14th ERCOFTAC SIG15 Workshops on refined turbulence modelling, Steiner et al. (2009) and Jakirlic´ et al. (2010b). In addition to different RANS models, the LES and LES-related methods (different seamless and zonal hybrid LES/RANS models; DES – Detached Eddy Simulation) were comparatively assessed (http://www. ercoftac.org). The latter schemes, hybridizing the RANS and LES methods aimed at a reduction of spatial and temporal resolution, have recently experienced growing popularity in the Computational Fluid Dynamics (CFD) community. Their goal is to combine the advantages of both methods in order to provide a computational procedure that is capable of capturing the large-scale eddy structures with a broader spectrum and the bulk flow unsteadiness – as encountered in the flows involving separation, but at affordable costs. Interested readers are referred to a relevant review about hybrid LES/RANS methods from Fröhlich and von Terzi (2008). We mention here only their general classification into two main groups: the zonal, two-layer schemes – a RANS model resolving the near-wall region is bridged at a distinct interface with the conventional LES covering the outer layer (flow core) – and the seamless models, where a RANS-like model formulation, mimicking a sub-scale model, is applied in the entire flow domain. The computational model examined in the present work belongs to the former category. This zonal hybrid LES/RANS model couples an eddy-viscosity-based RANS model in the wall layer to an LES in the outer flow region. In addition to the exper-

821

imental and DNS database, the results obtained are comparatively assessed to a complementary LES simulation. 2. 3-D diffuser: case description The diffuser shape, dimensions and the coordinate system are shown in Fig. 1. The flow in the inlet duct (height h = 1 cm, width B = 3.33 cm) corresponds to fully-developed turbulence (enabled experimentally by a development channel being 62.9 channel heights long). The L = 15h long diffuser section is followed by a straight outlet part (12.5h long). Downstream of this the flow goes through a 10h long contraction into a 1 in. diameter tube. The curvature radius at the walls transitioning between diffuser and the straight duct parts are 6 cm. The bulk velocity in the inflow duct is Ubulk = Uinflow = 1 m/s in the x-direction resulting in the Reynolds number based on the inlet channel height of 10,000. The origin of the coordinates (y = 0, z = 0) coincides with the intersection of the two non-expanding walls at the beginning of the diffuser’s expansion (x = 0). The working fluid is water ( q = 1000 kg/m3 and l = 0.001 Pas). 3. Computational method The continuity and momentum equations governing the incompressible flow in the present 3-D diffuser configuration read:

@U i ¼0 @xi @U i @ðU i U j Þ 1 @P @ ¼ þ ð sm þ st Þ þ @xj @t q @xi @xj ij ij

ð1Þ

  Here, smij ¼ 2mSij (with Sij ¼ 0:5 @U i =@xj þ @U j =@xi being the rate of strain tensor) represents the viscous stress tensor, whereas the turbulent stress tensor stij is to be modelled. The rationale and the most

Fig. 1. Geometry of the 3-D diffuser considered (not to scale), Cherry et al. (2008).

822

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

important features of the turbulence models used in the present work are outlined in the following sub-sections. 3.1. Computational models 3.1.1. Hybrid LES/RANS (HLR) model The present hybrid LES/RANS formulation represents a zonal, two-layer hybrid approach with a RANS model covering the near-wall region and the LES model the remainder of the flow domain. Both methods share the same temporal resolution (time step). The equations governing the velocity (Eq. (1)) field operate as the Reynolds-averaged Navier–Stokes equations in the near-wall layer (with U i representing the ensemble-averaged velocity field) or as the filtered Navier–Stokes equations in the outer layer (with U i representing the spatially filtered velocity field). The turbulent stress tensor stij in Eq. (1) representing either the sub-grid-stress tensor (stij  sij ¼ 2mt Sij  skk dij =3) or the Reynolds-stress tensor (stij  ui uj ¼ 2mt Sij  2kdij =3) is expressed in terms of the mean strain tensor via Boussinesq’s relationship. The equations governing the velocity field in the hybrid LES/RANS framework read:

" !#  DU i 1 @p @ @U i @U j þ ðm þ m t Þ þ ¼ Dt q @xi @xj @xj @xi

ð2Þ

The coupling of the instantaneous LES field and the ensembleaveraged RANS field at the interface is realized via the turbulent viscosity, which makes it possible to obtain solutions using one system of equations. This means practically that the governing equations are solved in the entire solution domain irrespective of the flow sub-region (LES or RANS). Depending on the flow zone, the hybrid model implies the determination of the turbulent viscosity mt either from a near-wall k–e RANS model: mt = Clflk2/e or from the subgrid-scale (SGS) model in the LES formulation (standard Smagorinsky model was applied, 1963): mt ¼ mSGS ¼ ðC S DÞ2 jSj. The Smagorinsky constant CS takes the value of 0.1. D = (Dx  Dy  Dz)1/3 represents the filter width and jSj ¼ ðSij Sij Þ1=2 the strain rate modulus. The near-wall variation of the turbulent viscosity mt is obtained from a two-equation, k–e RANS model, implying solution of the transport equations for kinetic energy of turbulence k and its dissipation rate e. The near-wall and viscous damping functions and the relevant source terms are presently modelled in line with the proposal of Launder and Sharma (LS, 1974), with e representing a ‘‘homogeneous” part of the total viscous dissipation rate ðe ! ~e ¼ e  2m 1=2 1=2 ð@k =@xk Þð@k =@xk ÞÞ, taking zero value at the wall (~ewall ¼ 0). This near-wall model is a straightforward extension of the standard high-Reynolds number k–e model aimed at capturing the enhanced viscosity influence on turbulence in the immediate wall vicinity (viscous sub-layer and buffer layer). Solving the equations for k and e in the RANS sub-domain (being necessary for the turbulent viscosity determination) requires appropriate boundary conditions at the RANS-LES interface. The approach used presently is based on the equality of the RANS interface values for k and e with corresponding SGS values. This condition, originating from the equivalence of the total stresses at both sides of interface (the RANS total stress, similar as the LES one, comprises both resolved and RANS modelled parts) - sLES tot ¼ stot , reduces to the equality of their modelled contributions by assuming the continuity of the resolved fractions across the interface (ifce), Temmerman et al. (2005): 2mSGS Sij ¼ 2mt Sij . Finally, the condition applied at the interface implies the equality of the (modelled) turbulent viscosities:

mt;ifce jRANS-side ¼ mSGS;ifce jLES-side

esis assuming the equality between the production, dissipation and energy flux through the cutoff and the analogy to the eddy-viscosity-based RANS modeling in line with the proposal of Mason and Callen (1986):

kSGS ¼

ðC S DÞ2 jSj2 pffiffiffiffiffiffi ; Cl

eSGS ¼ ðC S DÞ2 jSj3

ð4Þ

with Cl = 0.09. The algorithm being used to numerically ensure that the condition expressed by Eq. (3) is fulfilled, resembles the wellknown procedure for setting the value of a variable at a computational node and is explained as follows. The RANS equations for k and e are solved in the entire flow field, but with the discretization coefficients taking zero values in the LES sub-region (Ak,SGS = 0). By appropriately manipulating the corresponding Pthe source  terms in P discretization equations k Ak þ SP SGS UP;SGS ¼ k Ak Uk;SGS þ SU;SGS by taking SU,SGS = 1030USGS and SP,SGS = 1030, with USGS  kSGS, eSGS, their numerical solution in the framework of the finite volume method provides the interface values of the kRANS and e RANS, being equal to the corresponding SGS values. By doing so, the boundary condition at the LES/RANS interface (ifce), expressed by Eq. (3), is implicitly imposed without introduction of any empirical ‘‘transitional” function. Both the RANS model and the SGS model are used in their original form without any further adjustment or modification. This is illustrated in Fig. 2 displaying the variation of the modelled turbulent viscosity across the interface at two streamwise positions x/h = 2 (inflow duct) and x/h = 10 (diffuser section). One can clearly see the damping of the RANS viscosity towards the one of the SGS model. A gradual adjustment of the turbulent viscosities at both interface sides to each other ensures their smooth transition. Fig. 2 also offers evidence that the adopted grid resolution for both HLR and LES is adequate (see corresponding discussion in Section 3.3). Unlike some seamless hybrid LES/RANS models, as e.g., Detached Eddy Simulation (DES, Spalart et al., 1997; Spalart, 2009), the present zonal model treats the LES/RANS interface explicitly. Its position in the flow field plays an important role with respect to the quality of the solution. The present hybrid LES/RANS model is ‘‘equipped” with an algorithm providing an adaptive, flow-dependent interface location, controlled by a parameter expressing the ratio of the modelled to the total kinetic energy of turbulence, averaged appropriately over all grid cells in the LES-region of interface (Jakirlic´ et al., 2009). However, this procedure is not applied in the present simulation. The interface position along all diffuser walls, that is the RANS layer thickness, is presently estimated on the basis of the introductory computations of the 3-D duct flow and the LES results (see Section 3.3 dealing with the computational details).

RANS

LES

ð3Þ

Because k and e are not provided within the LES sub-domain in the case of the sub-grid-scale model of Smagorinsky, their SGS values are presently estimated utilizing the local equilibrium hypoth-

Fig. 2. Variation of modelled turbulent viscosity across the RANS/LES interface at two selected streamwise locations in the diffuser configuration.

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

The present HLR model has been intensively validated in configurations of different geometrical complexity featuring flow separation (from sharp-edged and continuous curved surfaces) and strong heating under the conditions of constant and variable fluid properties. Also, different sub-grid-scale (SGS) model formulations in LES (besides the Smagorinsky model the transport model for the residual kinetic energy from Yoshizawa and Horiuti (1985) has been tested) were combined with different eddy-viscosity-based RANS models. More details about the present HLR method are given in Jakirlic´ et al. (2006, 2009) and Kniesner (2008).

3.1.2. Large-Eddy Simulation (LES) The sub-grid scale motion was modelled by the Smagorinsky (1963) formulation (mSGS ¼ C 2S D2 jSj) utilizing the dynamic determination of the model coefficient C g ¼ C 2S proposed by Germano et al. (1991). The determination of the coefficient Cg and the implementation of the dynamic procedure into the code FASTEST (Section 3.2) are briefly outlined in the following. A detailed description can be found in Ertem-Mueller (2003). b > D to the Application of a test-filter with the filter length D filtered momentum equations leads to the sub-test-filter stress tensor T test in analogy to the sub-grid stress tensor sSGS ij ij :

cc T test ¼ Ud i Uj  Ui Uj ij

s

SGS ij

1 test b2 b c T test  dij T test ij kk ¼ 2C g D j Sj Sij ¼ 2C g aij 3 1 2 SGS sSGS  dij sSGS ij kk ¼ 2C g D jSjSij ¼ 2C g aij 3

ð6Þ

^ are obtained by The (test-) filtered quantities denoted with ‘‘ðÞ” applying the test-filter to each control volume (CV) as the volumeweighted average of the CV value itself and its twenty-six (26) neighbours. Beside the cell volumes additional weighting factors are taken into account. The variable arrangement and the test-filter with corresponding weighting factors are depicted in Fig. 3. The resolved part of the turbulent stresses Lij is computed in accordance with Germano (1992):

cc Lij ¼ T test  sSGS ¼ Ud i Uj  Ui Uj ij ij

ð7Þ

This identity is used to determine the model coefficient Cg dynamically in conjunction with the instantaneous local flow conditions. Utilizing the model assumptions used for the Smagorinsky formulation results in the following expressions:

ð8Þ ð9Þ

Insertion of the above formulations in the identity (7) yields

dSGS Lij ¼ 2C g atest ij þ 2 C g aij

ð10Þ

In order to compute the model coefficient the following approximation SGS aSGS ¼ C g ad C gd ij ij

ð11Þ

is utilized, being equivalent to the assumption that Cg is constant within the volume of the test-filter (see Fig. 3). This leads to

  SGS Lij ¼ 2C g ad  atest ¼ 2C g M ij ij ij

ð12Þ

According to Germano et al. (1991) Eq. (12) can be contracted with the resolved strain rate Sij resulting in Lij Sij ¼ 2C g M ij Sij . Since Mij can become zero, Cg may be indeterminate. To overcome this the least-squares method is used, as proposed by Lilly (1992). Application of the square error minimization algorithm leads to the following relation defining the preliminary model coefficient C g :

ð5Þ

¼ Ui Uj  Ui Uj

823

C g ð~ x; tÞ ¼

Lij M ij 2M ij M ij

ð13Þ

The determination of the model coefficient is carried out for each SIMPLE (Section 3.2) iteration while computing the new time level tn+1. To obtain a smooth behaviour of Cg and to avoid its wellknown oscillations, a relaxation of the coefficient C g in time with the factor b 2 ð0; 1 is carried out utilizing its value from the previous time step at tn:

 nþ1  ~ C  x; t n Þ þ bC g ð~ x; t nþ1 Þ ¼ ð1  bÞC g ð~ g x; t

ð14Þ

Since C  g can become negative, leading to a negative effective viscosity and to unphysical behaviour and numerical instabilities, the negative values of Cg are clipped. Furthermore, an upper bound of Cg is introduced to ensure the model does not become too diffusive. This finally yields the model coefficient in the corresponding expression for the turbulent viscosity enabling a converged solution at the new time level tnþ1

(a) control volume with neighbouring nodes

(b) Test filter Δ

b with corresponding weighting factors (b). Fig. 3. Variable arrangement (a) and applied test-filter D

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

824

n

h

i o C g ð~ x; t nþ1 Þ ¼ min max C  g ;0 ;1

ð15Þ

3.2. Numerical method The computational results were obtained using the in-house code FASTEST (Flow Analysis by Solving Transport Equations Simulating Turbulence), which uses a finite volume method for blockstructured, body-fitted, non-orthogonal, hexahedral meshes (FASTEST-Manual, 2005). Block interfaces are treated in a conservative manner, consistent with the treatment of inner cell-faces. Cell centred (collocated) variable arrangement and Cartesian vector and tensor components are used. The equations are linearised and solved sequentially using an iterative method. The velocity–pressure coupling is ensured by the pressure-correction method based on the SIMPLE algorithm which is embedded in a geometric multi-grid scheme with standard restriction and prolongation and an ILU for smoothing, Durst and Schäfer (1996). To avoid decoupling of pressure and velocity on the collocated grid the selective interpolation method proposed by Rhie and Chow (1983) is applied. To enable simulations with a large number of control volumes the code is strictly parallelised. The parallelisation is based on a domain decomposition, which is directly related to the block structure of the spatial discretisation. By using an automatic partitioning tool, an optimal load balancing of the processors can be achieved, even

for complex geometries with hundreds of blocks, Sternel et al. (2004). Communication between the processors is organised via the Message Passing Interface (MPI) technique. The convective transport of all variables was discretised by a second-order, central differencing scheme for LES and DES. In the case of the HLR method some upwinding is used for k and e equations by applying the socalled ‘‘flux blending” technique. Time discretisation was accomplished applying the Crank-Nicolson scheme. 3.3. Computational details 3.3.1. LES The solution domain comprising a part of the development duct (5 h), the diffuser section (15 h), the outlet straight duct (12.5 h) as well as the convergent part (9 h), is meshed with 3.55 Mio. grid cells: Nx  Ny  Nz = 408  64  136, Fig. 4. Two simulations with and without the SGS model have been performed (only the results obtained by using the SGS model are shown here). The wall boundary layers are resolved with y+ values of O(1). The off-wall resolution is close to that used in HLR simulation, see Figs. 2, 5 and 7 and corresponding discussion. According to the experimental setup of Cherry et al. (2008) the fully-developed turbulent channel flow has been computed with respect to the inflow generation. These inlet data are generated by a simultaneously running periodic channel flow simulation of a channel (Nx  Ny  Nz = 48  64  136)

Fig. 4. Solution domain and computational mesh used in the LES simulation (every 5th grid line is plotted).

Fig. 5. Near-wall resolution expressed in wall units along the upper and lower walls in the central x–y plane (z/B = 1/2) in the HLR simulation.

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

825

Fig. 6. Interface position in the HLR simulation corresponding to y+  50.

with the same cross-section as the diffuser inlet; small figure in Fig. 4 (see also Fig. 8). To allow the flow through the diffuser to influence the flow field in the development channel a part (5 channel heights) of this channel has been modelled in front of the diffuser. The turbulent flow fields in a cross-section in the periodic channel are copied to this inlet location. We are aware of the fact that the 3-D streamwise-periodic channel of length 4.5h used for the inflow generation might be too short, keeping in mind the spatial extent of the characteristic vortex structures, which is in general larger (due to the secondary currents) than in a channel flow with the spanwise homogeneity. Furthermore, Nikitin (2008) argued that an auxiliary streamwiseperiodic simulation might not be suitable since it causes a spatial periodicity, which is not physical for turbulent flows. Let us recall that the solution domain in the DNS of Ohlsson et al. (2010) comprises the inflow development duct of 63h length, accounting even for the transition of the initially laminar inflow. The present simplification of the numerical setup is especially pertinent to the hybrid LES/RANS method, since the overall aim is to improve the efficiency (lower computational costs) and applicability to complex geometries. In order to achieve the same basis for mutual comparison of the presently employed LES and HLR, both methods used the same inflow conditions, i.e. the same inflow duct length. The emphasis of the present study was on the time-averaged flow field (see Section 4) which was in reasonable agreement

with the reference databases, retroactively justified the simplified inflow generation. We mention here also the comparable inflow generation of Schneider et al. (2010), successfully applied over a channel length of only 3 channel heights. The temporal resolution adopted corresponds to a non-dimensional time step of Dt = 0.0001Ubulk/h = 0.01 (normalized by the inlet channel parameters Ubulk = 1 m/s and h = 1 cm). This leads to the maximum CFL number being substantially smaller than one. As a comparison the value of the viscous time scale m=U 2s in the inflow region corresponds to a dimensionless time step of Dt = 0.028 (the value of friction velocity at the bottom wall of the inflow channel at z/B = 0.5 obtained by DNS was about 0.063 – this values can be extracted from Fig. 10-left; significantly lower values were documented in the diffuser section). Averaging was performed for approximately 80,000 time steps. The simulations were carried out on 16 IBM Power 5 CPUs with a load balancing efficiency of 100%. 3.3.2. HLR The inflow data were, in analogy to the present LES, generated by a precursor simulation of the fully-developed duct flow (with the auxiliary periodic duct length of 6.5h) using the respective models on the grid of Nx  Ny  Nz = 78  62  134 cells. The solution domain for the HLR simulation comprised a part of the development duct (5 h), the diffuser section (15 h) and the straight outlet duct (12.5 h), Fig. 1. At its outlet cross-section

Fig. 7. Profiles of the kinetic energy of turbulence at two selected streamwise locations (x/h = 2 – inlet duct and x/h = 10 diffuser section) obtained by the present HLR method (vertical lines denote the LES/RANS interface position).

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

826

Fig. 8. The instantaneous velocity field obtained by LES (small figure on the upper left illustrates the turbulent inlet data generation through a separate inflow duct simulation utilizing a simultaneously running periodic channel simulation).

(a) experiment:

(b) HLR:

0.1 m/s

0.1 m/s

Fig. 9. Velocity vectors in the y–z plane in the inflow duct.

the convective outflow conditions were applied. No-slip boundary conditions were applied at the walls. The grid applied contained Nx  Ny  Nz = 224  62  134 cells (approximately 1.86 Mio. grid cells in total; the number of grid cells in the cross-plane y–z as well as of those cells in the diffuser section correspond to the number applied in the LES simulation). A dimensionless time step of Dt ¼ 0:00028U bulk =h ¼ 0:028 was used for the computations, ensuring a CFL number less than unity throughout the entire solution domain (CFLmax  0.76). This is of the same order of magnitude as the viscous time scale m=U 2s (see discussion related to the LES-relevant computational details). The near-wall resolution corresponds to Dy+ < 0.8, Dx+  10–100 and Dz+  2–20 (along the x–y plane at z/B = 1/2) for the HLR simulation, Fig. 5.

Similar resolution was achieved on the side walls (not shown here). The final position of the interface, separating the near-wall (RANS; black area in Fig. 6) and off-wall (LES; gray area in Fig. 6) sub-regions, corresponds to the dimensionless distances from the bottom/upper and side walls y+, z+  50. This position was determined in accordance to the preliminary computations of the fully-developed flow in the inflow duct as well as the results obtained by the pure LES simulations. It should be noted that the symmetry plane in the central vertical plane of the inflow duct is at y+  300 (see position of the velocity maximum in Fig. 10-left). Accordingly, the wall layer thickness covered by the RANS model amounts to one sixth of the half-height of the inflow duct. Averaging was performed for approximately 80,000 time steps. The simulations were carried out on the same computer as LES but using only seven processors. One of the important issues when employing a hybrid LES/RANS method is the computational time used relative to the complementary LES simulation. Such an analysis could not be performed straightforwardly in the present work keeping in mind the differences in the grid size, model structure (HLR model employs the standard Smagorinsky model and LES its dynamic version) and number of processors. But for instance, in a fully-developed channel flow (the analysis was conducted by computing the flow at Res = 640), the present hybrid model needs about 10–30% more time (depending on the version of the low-Reynolds number, twoequation k–e model covering the RANS layer) if using the same grid, see Kniesner (2008). Further information about the suitability of the grid resolution can be gained from Figs. 2 and 7 by analysing the ratio of the turbulent SGS viscosity to molecular viscosity and the SGS contribution to total turbulence kinetic energy at the locations across the inflow duct (x/h = 2) and the interior of the diffuser (x/h = 10). Although extracted from the HLR simulation, these results could serve as the basis for discussion about the grid resolution assessment also for the LES simulation. Fig. 2 illustrates the SGS viscosity

Fig. 10. Semi-log plots of axial velocity component across the inflow duct (x/h = 2) and at a cross-section in the interior of the diffuser (x/h = 10).

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

Fig. 11. Surface pressure distribution at the bottom wall.

being substantially smaller than its molecular counterpart in the LES sub-region. Fig. 7 demonstrates clearly the contribution of the kSGS (obtained from Eq. (4)) to the total k-value being less than 10% everywhere in the LES sub-region (compatible with the one in the pure LES simulation); the RANS layer thickness corresponds approximately to 0.1h at both streamwise locations (the vertical lines denote the LES/RANS interface position). One can also see the fraction of the modelled turbulence in the RANS layer being between 20% and 25% of the total turbulence kinetic energy. The exception is the portion of the upper wall in the interior of the diffuser section being affected by separation, where almost the entire turbulence kinetic energy originates from the resolved fluctuations. 4. Results and discussion Fig. 8 displaying the instantaneous velocity field obtained by LES, provides a first impression about the flow topology. A slice through the computational domain and the periodic channel is depicted, showing the contours of the instantaneous streamwise velocity component at an arbitrary instant of time. The flow in the inflow duct complies with the fully-developed turbulence underlying approximately the equilibrium conditions (see Fig. 10 and corresponding discussion). Increased flow perturbations start already at the corner between the flat walls of the inflow duct and the sloped diffuser walls, causing the deceleration of the boundary layer and a deformation of the mean velocity profile. The boundary layers along the flared walls interact and separate finally in the corner region under the influence of the adverse-pressure gradient. The separation bubble grows and gradually occupies the entire upper diffuser wall. Fig. 8 illustrates clearly the unsteady nature of the diffuser flow, characterized by flapping motion of the separated shear layers and multiple separating and reattaching regions. A selection of the computational results along with the experimental and DNS data is displayed in Figs. 9–15. The presentation and analysis of the results obtained has been conducted with respect to flow in the inflow duct, the size and shape of the flow separation pattern and associated mean flow and turbulence features. 4.1. Inflow duct Unlike the flow through a circular pipe, the flow in a duct with rectangular cross-section is no longer unidirectional. It is characterized by a strong secondary motion with the velocity components being perpendicular to the axial direction, Fig. 9. This secondary flow transporting momentum into the duct corners is characterized by jets directed towards the duct walls bisecting each corner with associated vortices at both sides of each jet. This secondary current resembles Prandtl’s flow of the second kind

827

(possible only for turbulent flows) induced by the Reynolds stress anisotropy. Indeed, the Reynolds stress gradients cause the generation of forces which induce the normal-to-the-wall velocity components in the secondary flow plane. Accordingly, correct capturing of anisotropic turbulence in the inflow duct is an important prerequisite for a successful computation of the diffuser flow. Fig. 9 displays the time-averaged velocity vectors in the cross-plane y–z located at x/h = 2 obtained experimentally and computationally by the HLR model. Despite the lower intensity of the secondary motion, whose largest velocity has the magnitude corresponding approximately to <(1–2)% of the axial bulk velocity (Ubulk = 1 m/ s), its influence on the flow in the diffuser is significant. The qualitative agreement with respect to the secondary flow topology discussed above is obvious. Fig. 10-left depicts the semi-log plot of the axial velocity component across the central plane (z/B = 0.5) of the inflow duct at x/h = 2. The inflow conditions correspond clearly to those typical for an equilibrium flow. The velocity profile shape obtained by DNS follows closely the logarithmic law, despite a certain departure from it. This departure, expressed in terms of a slight underprediction of the coefficient B in the log-law (U+ = ln (y+)/ j + B; B = 5.2), can also be regarded as a consequence of the back-influence of the adverse-pressure gradient evoked by the flow expansion. The pressure coefficient evolution, displayed in Fig. 11, reveals an appropriate pressure increase already in the inflow duct (x/h 6 0). The present computational results exhibit a certain overprediction of the velocity in the logarithmic region. It resembles a typical outcome in the case of a somewhat coarser grid resolution. One would assume a sufficiently fine grid in this region according to the ratio of the turbulent viscosity to the molecular viscosity (Fig. 2) which is less than 1. On the other hand, the profile of the kinetic energy of turbulence (Fig. 7-left) reveals a somewhat higher level of the modelled turbulence in this attached flow region. All these issues contributed to an appropriate underprediction of the friction velocity Us, serving here for the normalization – U+ = U/Us. Whereas the LES returned a smooth profile, the velocity profile obtained by the present HLR method exhibits a certain bump in the region around the LES/RANS interface (y+  50). It represents an outcome typical for the application of a zonal hybrid LES/RANS method. The relatively high turbulent viscosity in the RANS layer strongly damps the fluctuations, suppressing their recovery in the LES sub-region until some distance behind the interface. Such an anomaly can be remedied by generating some additional fluctuations and superimposing them onto the LES ones (see Jakirlic´ et al. (2010a) for specific details). This so-called ‘‘forcing” procedure was not used in the present work. This is justified by the fact that the ‘‘buffer layer” associated with the interface region is limited only to the immediate wall vicinity in the inflow duct, and does not significantly affect other quantities in this region (see e.g., smooth behaviour of the kinetic energy profile in Fig. 7-left) or the remainder of the flow domain (see e.g., smooth evolution of the velocity and kinetic energy profiles at x/h = 10 across the interface in Figs. 7-right and 10-right; both profiles displaying a behaviour, typical of a flow affected by an adverse-pressure gradient – underprediction of the logarithmic law and strong enhancement of the turbulence intensity corresponding to a substantial mean flow deformation in the streamwise direction). 4.2. Diffuser section Fig. 11 illustrates the evolution of the surface pressure coefficient along the bottom non-expanding wall in the diffuser configuration, representing the most important integral characteristic. The surface pressure distribution along the bottom wall was evaluated from the expression C P ¼ ðp  pref Þ=ð0:5qU 2bulk Þ; the reference pressure was taken at the position x/L = 0.05. The pressure

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

x/h=15

x/h=12

x/h=8

x/h=5

x/h=2

828

Experiment

LES

HLR

Fig. 12. Comparison between experimentally and computationally obtained iso-contours of the axial velocity field in the cross-planes y–z at selected streamwise locations within the diffuser section (the thick line denotes the zero-velocity line).

curve exhibits a development typical of flow in diverging ducts. The pressure decrease in the inflow duct is followed by a steep pressure increase already at the very end of the inflow duct and especially at the beginning of the diffuser section. The transition from the initial strong pressure rise to its moderate increase occurs at x/L  0.3, corresponding to the position where about 9% of the entire cross-section is occupied by the flow reversal (see e.g., Fig. 13). Onset of the separation zone causes a certain contraction of the flow cross-section, leading to a weakening of the deceleration intensity and, accordingly, to a slower pressure increase. The region characterized by a monotonic pressure rise was reached in the remainder of the diffuser section. The results obtained by both computational methods display this correct tendency, but some quantitative differences are present. The result obtained by applying the HLR method exhibit very good agreement with the experimental data, apart from the slight underprediction in the flow

region affected by a steep pressure gradient. The LES results are typical of a higher pressure gradient in the region x/h = 3–5 (consistent with an enhanced flow deceleration in the lower half of the diffuser section, Fig. 14) resulting in an oveprediction of the CP value over the entire diffuser length. Fig. 12 shows the contour plots of the axial velocity component at five streamwise cross-sectional areas indicating the evolution of the flow separation pattern. The recirculation zone development displayed here can be analyzed in parallel with the quantitative information about the fraction of the diffuser cross-sectional area occupied by the reversed flow depicted in Fig. 13. The adversepressure gradient is imposed onto the intersecting boundary layers along flat walls upon entering the diffuser section. According to the experimental investigation the boundary layers along all walls are of comparable thickness. The separation starts immediately after the beginning of the diffuser section at x/L = 0 (x/h = 0), see.

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

Fig. 13. Fraction of the cross-sectional area occupied by the flow reversal.

Fig. 13. The onset of separation pattern is located in the upper-right diffuser corner, built up from the deflected side wall and the top wall, see e.g., the position x/h  2 (x/L = 0.13) in Fig. 12. Initial growth of this corner bubble reveals its spreading rate along two sloped walls being approximately of the same intensity, see position x/h = 5. As the adverse-pressure gradient along the upper wall outweighs significantly the one along the side wall due to substantially higher angle of expansion, 11.3° vs. 2.56°, the separation zone spreads gradually over the entire top wall surface, see position x/ h = 8. One notes a strong three-dimensional nature of the separation pattern. The maximum occupation of the diffuser cross-sectional area by the flow reversal, around 22% (Fig. 13), is documented at the position x/h = 12–15 (x/L = 0.8–1.0). The thickness of the flow reversal zone (its dimension in the normal-to-wall direction) is almost constant over the diffuser width in this region, resembling approximately a 2-D pattern. After this position the intensity of the back-flow weakens. The experimental results point out the reattachment region being located around x/L  1.4 within the straight outlet duct, Fig. 13. The results obtained by the LES and HLR methods exhibit qualitatively close agreement with the previously described spatial development of the three-dimensional back-flow region, both in size and shape. This is also the case when comparing with the DNS database (Fig. 3 in Ohlsson et al. (2010)). As in the experimental configuration, the separation starts in the upper-right corner. Whereas the LES simulation results in a corner bubble covering both deflected walls (streamwise position x/h = 5 in Fig. 12), in close agreement with the experimental findings with respect to shape but not size, the HLR model response to the adversepressure gradient onset is manifested in an asymmetrical corner bubble occupying mostly the upper wall. The size of the cross-sectional area occupied by the latter separation zone is smaller than the experimentally detected one. This outcome is reflected in the lower fraction of the back-flow region compared to the throughflow area within the first half of the diffuser section, Fig. 13. Whereas the corner bubble predicted by LES preserves its symmetrical shape growing equally in both coordinate directions in the vertical cross-plane y–z at x/h = 8 (x/L = 0.53; one notes the appearance of a smaller bubble in the upper left corner – disconnecting separation zone), the HLR simulation results in the back-flow zone similar in shape to the experimental one. The total amount of the cross-sectional area affected by the flow reversal is about 10% and agrees well with the experimental outcome at this position, Fig. 13. Both computational methods yield a separation zone spreading over the entire upper wall from the streamwise location x/h = 8–10 (x/L = 0.5–0.7). Contrary to the experimental findings the computationally obtained separation zone exhibits an asymmetric form with a non-constant thickness over the diffuser duct

829

width. The position of the maximum reverse flow cross-section (about 20%) is reached at x/h = 12–15 (x/L = 0.8–1.0) in good agreement with experiment. After this location the fraction of the crosssectional area covered by the separated flow follows closely the experimentally determined results. Analogous to the experiment the flow reattaches at position x=L  1:4—1:5 in the straight outlet part of the diffuser configuration. In summary it can be stated that despite certain deviations the topology of the separation region, with respect to its onset, shape and size, is reasonably well reproduced by both computational methods. Good overall agreement obtained by HLR and LES is further confirmed in Figs. 14 and 15 illustrating the axial velocity and streamwise stress component profiles at 14 streamwise locations situated in all characteristic flow regions – attached flow in the inflow duct, separation point, recirculation zone, reattachment and recovery region – within the diffuser configuration in three vertical cross-planes corresponding to the spanwise locations z/B = 1/4, z/B = 1/2 and z/B = 3/4 (B = 3.33 cm). Both the mean velocity and turbulence fields have been mapped three-dimensionally, enabling a detailed quantitative insight into the entire flow domain. In addition to the reference experiments, the present computational results are analysed along with the DNS database from Ohlsson et al. (2010). Fig. 14 displays the velocity field development typical for the flow in an expanding duct. The bulk flow exhibits deceleration, leading to an asymmetry of the velocity profile. The effect of the adverse-pressure gradient is especially visible in the flow region along the upper expanding wall. The velocity profile approaches gradually the form characterizing a separating flow, exhibiting regions of the zero velocity gradient (at separation and reattachment points) and profile inflection. The through-flow, that is the flow in the positive streamwise direction, is characterized by a spreading dictated by the pressure gradient arising from the geometry expansion. Accordingly, the position of the reduced velocity maximum is gradually shifted towards the upper wall, eventually reaching the center (y/h  2) of the straight outlet channel (with the height 4h). In this post-reattachment zone the velocity profile exhibits a fairly flattened form, being almost symmetric. The consequence of the velocity profile flattening is a continuous monotonic decrease of the wall shear stress. The velocity profile evolution is similar in all three longitudinal vertical planes. The specific differences are related to the vicinity of both bottom and upper walls, especially the latter upper wall, where the flow passes regions of three-dimensional flow reversal (see Fig. 12 and corresponding discussion). In Figs. 15, the development of the streamwise turbulence intensity is shown. The lowest turbulence intensity is situated in the region coinciding with the mean velocity maximum – flow zone with approximately zero velocity gradient – along the entire diffuser section. The Reynolds stress profiles exhibit their highest values in the regions with the most intensive flow deformation. These are the near-wall layer in the attached flow regions and the flow zone along the shear layer bordering the recirculation zone. The peak of the turbulence intensity originating from the boundary layer at the top inflow duct wall increases initially, after the strong rise in the pressure field is imposed (see pressure coefficient development in Fig. 11), and weakens slightly after flow transition to the second part of the diffuser section, being characterized by a decreasingly adverse-pressure gradient. The streamwise turbulence intensity in the outlet duct is uniformly distributed over the cross-section. The present computational results are in close agreement with both reference data sets, experimental and DNS. This holds for the individual velocity and Reynolds stress profiles but also for the flow development as a whole. However, some discrepancies should not be overlooked. The largest ones with respect to the velocity field are found in the region aligned with the mean dividing

830

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

Fig. 14. Evolution of the axial velocity profiles in the vertical plane x–y at three spanwise locations z/B = 1/4, 1/2 and 3/4.

streamline. This is especially the case with the LES simulations. Further deviations, particularly with respect to the Reynolds stress profiles, are noticeable in the core flow within the transitional region between the diffuser section and the straight outlet channel

(x/h P 15), where the flow is exposed to a relaxation. Both present simulations were performed on the grids with the cells clustered adequately upon approaching the solid walls. The wall-closest numerical nodes are situated deeply in the viscous sub-layer along

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

831

Fig. 15. Evolution of the turbulent streamwise stress component profiles in the vertical plane x–y at three spanwise locations z/B = 1/4, 1/2 and 3/4.

all duct walls. The adopted expansion factor led to appropriately fine near-wall grid resolution and a somewhat coarser grid in the flow core, the latter feature representing probably the reason for the afore mentioned deviation.

Similar deviations of the LES simulations from the reference data are observed in the region characterized by high turbulence intensity. Both computational methods complement each other in the core flow region. The basic differences originate from

832

S. Jakirlic´ et al. / International Journal of Heat and Fluid Flow 31 (2010) 820–832

the near-wall treatment. Whereas the hybrid LES/RANS method employs an elaborate near-wall RANS model which solves additional transport equations for k and e for the velocity and length scale determination in the turbulent viscosity formulation, the LES method extracts both representative scales from the velocity field and the grid spacing. As elaborated earlier, the wall boundary layer is covered by an irregular, highly non-uniform grid, the fact being particularly applicable with respect to the streamwise and spanwise resolutions. The LES models assume a cutoff wave number well situated in the inertial sub-range, a condition not fulfilled in the immediate wall vicinity characterized by lowReynolds numbers, where production and dissipative ranges lie almost immediately next to each other. This situation is consistent with a very narrow inertial region. Unlike the present LES, the near-wall RANS model works well under conditions of non-uniform, anisotropic grids with higher aspect ratio. Generally speaking, an appropriate grid refinement leading eventually to the improved agreement would be easily manageable; keeping in mind a fairly low flow Reynolds number (Reh = 10,000) a moderate increase in the number of grid cells would probably be sufficient to accomplish this task. However, as the primary objective was the validation of the hybrid LES/ RANS method designed to function under the condition of a coarser spatial resolution in such a complex flow configuration featured by three-dimensional separation, the presently used grid was intentionally adopted.

5. Conclusions Two eddy-resolving modelling approaches: Large-Eddy Simulation (LES) and a zonal hybrid LES/RANS scheme (HLR) were used to predict the flow in a 3-D diffuser, measured recently by Cherry et al. (2008, 2009) and computed by means of Direct Numerical Simulation by Ohlsson et al. (2009, 2010), aiming at a comparative analysis of their features and performance. The latter method, developed by the authors, merging a low-Reynolds-number k–e Reynolds-Averaged Navier–Stokes (RANS) model with Large-Eddy Simulation (LES) in the framework of a two-layer scheme, was the focus of the present investigation, keeping in mind that the basic motivation behind the HLR scheme is to obtain results on a coarser grid, whose quality is comparable with that of a well-resolved LES. In order to check its feasibility, an intensive validation has been performed. The flow configuration exhibited a 3-D recirculation pattern was thought to represent a suitable benchmark. Good results were obtained by both computational models with respect to the characteristics of the duct flow expanding into a diffuser section, the consequent separation flow region (onset, shape and size), the mean velocity field and associated integral parameters (pressure distribution), as well as the turbulence quantities. More importantly, a close agreement was achieved between the hybrid LES/RANS model and both experimental and DNS reference databases with respect to the near-wall behaviour of the timeaveraged quantities, demonstrating its high potential for computing such complex, wall-bounded turbulent flows.

Acknowledgements The work of G. Kadavelil and M. Kornhaas is financially supported by the Deutsche Forschungsgemeinschaft (DFG) within the German Collaborative Research Center (SFB 568) ‘‘ Flow and Combustion in Future Gas Turbine Combustion Chambers”. Our special thanks goes to the authors of the reference results J. Eaton and E. Cherry (experimental database) and J. Ohlsson (DNS database) for making their data available.

References Bradshaw, P., Wong, F.Y.F., 1972. The reattachment and relaxation of a turbulent shear layer. J. Fluid Mech. 52 (1), 113–135. Cherry, E.M., Elkins, C.J., Eaton, J.K., 2008. Geometric sensitivity of threedimensional separated flows. Int. J. Heat Fluid Flow 29, 803–811. Cherry, E.M., Elkins, C.J., Eaton, J.K., 2009. Pressure measurements in a threedimensional separated diffuser. Int. J. Heat Fluid Flow 30, 1–2. Durst, F., Schäfer, M., 1996. A parallel block-structured multigrid method for the prediction of incompressible flows. Int. J. Numer. Methods Fluids 22, 549–565. Ertem-Mueller, S., 2003. Numerical efficiency of implicit and explicit methods with multigrid for Large Eddy Simulation in complex geometries. Dissertation. Technische Universität Darmstadt. FASTEST-Manual, 2005. Institute of Numerical Methods in Mechanical Engineering, Department of Mechanical Engineering, Technische Universität Darmstadt, Germany. Fröhlich, J., von Terzi, D., 2008. Hybrid LES/RANS methods for the simulation of turbulent flows. Prog. Aerosp. Sci. 44, 349–377. Germano, M., 1992. Turbulence: the filtering approach. J. Fluid Mech. 238, 325–336. Germano, M., Piomelli, U., Moin, P., Cabot, W.H., 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A 3 (7), 1760–1765. Jakirlic´, S., Kniesner, B., Šaric´, S., Hanjalic´, K., 2006. Merging near-wall RANS models with LES for separating and reattaching flows. In: ASME Joint US–European Fluids Engineering Summer Meeting: Symposium on ‘‘DNS, LES and Hybrid RANS/LES Techniques”, Paper No. FEDSM2006-98039, Miami, FL, USA, July 17–20. Jakirlic´, S., Kniesner, B., Kadavelil, G., Gnirß, M., Tropea, C., 2009. Experimental and computational investigations of flow and mixing in a single-annular combustor configuration. Flow Turbul. Combust. 83 (3), 425–448. Jakirlic´, S., Kniesner, B., Kadavelil, G., 2010. A method for interface-turbulence forcing in hybrid LES/RANS simulations. In: 8th Int. ERCOFTAC Symp. on ‘‘Engineering Turbulence Modelling and Measurements”, Marseille, France, June 9–11. Jakirlic´, S., Kadavelil, G., Sirbubalo, E., von Terzi, D., Breuer, M., Borello, D., 2010. Report on 14th ERCOFTAC Workshop on Refined Turbulence Modelling. September 18, 2009. ‘Sapienza’ University of Rome, ERCOFTAC Bulletin, No. 84, 2010. Kniesner, B., 2008. Ein hybrides LES/RANS Verfahren für konjugierte Strömung, Wärme- und Stoffübertragung mit Relevanz zu Drallbrennerkonfigurationen (A hybrid LES/RANS method for conjugated flow, heat and mass transfer with relevance to swirl combustor configurations). PhD Thesis, Technische Universität Darmstadt. . Launder, B.E., Sharma, B.I., 1974. Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disc. Lett. Heat Mass Transfer 1, 131–138. Lilly, D.K., 1992. A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids A 4 (3), 633–635. Mason, P.J., Callen, N.S., 1986. On the magnitude of the subgrid-scale eddy coefficient in large-eddy simulation of turbulent channel flow. J. Fluid Mech. 162, 439–462. Nikitin, N., 2008. On the rate of spatial predictability in near-wall turbulence. J. Fluid Mech. 614, 495–507. Obi, S., Aoki, K., Masuda, S., 1993. Experimental and computational study of turbulent separating flow in an asymmetric plane diffuser. In: Proc. 9th Symp. on Turbulent Shear Flows, Kyoto, Japan, August 16–19, pp. 305.1–305.5. Ohlsson, J., Schlatter, P., Fischer, P.F., Henningson, D.S., 2009. DNS of threedimensional separation in turbulent diffuser flows. In: Advances in Turbulence XII. Proceedings of the 12th EUROMECH European Turbulence Conference, vol. 132. Springer Proceedings in Physics, Marburg. ISBN: 978-3-642-03084-0. Ohlsson, J., Schlatter, P., Fischer, P.F., Henningson, D.S., 2010. DNS of separated flow in a three-dimensional diffuser by the spectral-element method. J. Fluid Mech. 650, 307–318. Rhie, C.M., Chow, W.L., 1983. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation. AIAA J. 21, 1525–1532. Schneider, H., von Terzi, D., Bauer, H.-J., Rodi, W., 2010. Reliable and accurate prediction of three-dimensional separation in asymmetric diffusers using largeeddy simulation. ASME J. Fluids Eng. 132 (3), 031101–031111. Smagorinsky, J., 1963. General circulation experiments with the primitive equations: I. The basic experiment. Mon. Weather Rev. 91 (3), 99–164. Spalart, P.R., 2009. Detached-Eddy simulation. Annu. Rev. Fluid Mech. 41, 181–202. Spalart, P.R., Jou, W.-H., Strelets, M., Allmaras, S., 1997. Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach. In: 1st AFOSR Int. Conf. on DNS and LES, Columbus, OH, USA, August 4–8. Steiner, H., Jakirlic´, S., Kadavelil, G., Šaric´, S., Manceau, R., Brenn. G., 2009. Report on 13th ERCOFTAC Workshop on Refined Turbulence Modelling. September 25–26, 2008. Graz University of Technology, ERCOFTAC Bulletin, No. 79, pp. 24–29. Sternel, D.C., Junglas, D., Martin, A., Schäfer, M., 2004. Optimisation of partitioning for parallel flow simulation on block structural grids. In: Topping, B.H.V., Mota Soares, C.A. (Eds.), Proc. of the 4th Int. Conf. on Engineering Computational Technology, Stirling: Paper 93. Civil-Comp Press, Stirling. Temmerman, L., Hadziabdic, M., Leschziner, M.A., Hanjalic, K., 2005. A hybrid twolayer URANS-LES approach for large eddy simulation at high Reynolds numbers. Int. J. Heat Fluid Flow 26, 173–190. Yoshizawa, A., Horiuti, K., 1985. A statistically-derived subgridscale kinetic energy model for the large-eddy simulation of turbulent flows. J. Phys. Soc. Jpn. 54, 2834–2839.