Simulation of the Flow in a Ribbed Rotating Duct with a Hybrid k-ω RANS/LES Model

The objective of this work is to verify the capabilities of a hybrid k-ω RANS/LES model for simulation of the unsteady three-dimensional flow in a ribbed duct subjected to system rotation. The Reynolds number is 15,000 and the rotation number is 0.3, both based on hydraulic diameter and bulk velocity. A correction term for system rotation is introduced into the originating k-ω RANS model. Simulation results in the mid-span section are compared with experimental data by Coletti et al. (Exp. Fluids 52:1043–1061, 2012). The comparison is complemented by analysis of the flow features in cross-sections. It is demonstrated that the hybrid k-ω RANS/LES model produces an accurate simulation of the rotating ribbed duct flow. Results are compared with those by the originating time-accurate k-ω RANS model. The k-ω RANS model is not accurate concerning secondary features in the longitudinal mean flow recirculation patterns and the secondary flow in cross-sections, but it reproduces quite well the time-averaged longitudinal flow.


Introduction
Study of the flow in a ribbed rotating duct is relevant for technical applications, such as the flow in blade passages of radial compressors and in internal cooling channels of rotating S. Kubacki slawomir.kubacki@meil.pw.edu.pl 1 blades in axial turbine stages of gas turbine engines. The rotating duct with ribs is also a challenging test case for validation of turbulence models and LES methods. Recent papers with a discussion of the flow phenomena and a summary of the broad literature on experimental and numerical techniques for analysis of flow in ribbed rotating ducts are, e.g., by Narasimhamurthy and Andersson [1] and by Elfert et al. [2].
The basic phenomena of rotating-duct flow were already described by Johnston [3]. A first feature is amplification or damping of turbulence in the boundary layers. If the Coriolis force points towards a wall, pressure builds up and turbulence is enhanced in wall vicinity. Because the magnitude of the Coriolis force is proportional to the velocity magnitude, motion of a fluid element with instantaneous high longitudinal velocity, towards the low velocity zone at a wall, is amplified by the Coriolis force. Similarly, motion of a fluid element with instantaneous low longitudinal velocity, away from wall vicinity towards the centre of the channel, is amplified as well. Consequently, turbulence is enhanced. The boundary layer at the pressure side is then said to be destabilised, which means that it becomes more turbulent. If the Coriolis force points away from a wall, this wall becomes the suction side and turbulence is damped. The boundary layer is then said to be stabilised, which means that the flow tends to be more laminar. A second Coriolis force induced phenomenon is vortex flow patterns in cross-sections, called secondary flow. The pressure difference between the trailing (pressure side) and leading (suction side) walls of the rotating duct, caused by the Coriolis force, is approximately the same in the core of the flow and in the side-wall boundary layers. But the reduced velocity in the side-wall boundary layers means a reduced Coriolis force. The local unbalance between the Coriolis force and the transverse pressure gradient causes a transversal flow in the side-wall boundary layers, from the pressure side to the suction side.
The first experimental techniques for analysis of rotating ribbed ducts were by laser-Doppler velocimetry (LDV). A well-known example is the work by Iacovides et al. [4], who measured the velocity field in a longitudinal section of a rotating U-bend with orthogonal square ribs on the leading and trailing walls, staggered to each other. Mean and fluctuating velocity components were recorded along lines perpendicular to the ribbed walls. A recent example of measurements by particle image velocimetry (PIV) is by Elfert et al. [2], who analysed secondary flow patterns in a two-pass ribbed channel representative for the front part of an internal blade cooling system. Such an analysis is possible with PIV, but practically unrealisable with LDV. With the cited measurements, the experimental facility was stationary with respect to the rotating duct, as it is in the vast majority of experiments. This means that illumination of the measurement point or the measurement plane has to be synchronised to a chosen circumferential position of the rotating test model and that the number of recordings per passage of the model is very limited. Further, the relative velocity of the internal flow is obtained by subtracting the peripheral speed from the measured absolute velocity and there is uncertainty on the peripheral speed and on the position of the measurement plane. The consequence is limited temporal and spatial resolution and rather large uncertainties, resulting in inaccurate detection of small-scale unsteadiness and high velocity gradients. Much more accurate measurements are possible with the PIV system rotating together with the duct, but such measurements have only been realised up to now by two teams. Visscher et al. [5] performed PIV measurements with the equipment rotating together with a long smooth duct. But the duct has a large aspect ratio (width to height is 5:1), and the Reynolds number is much lower than with flow through cooling channels of gas turbine blades. The objective was fundamental study of the effect of the Coriolis force on turbulence enhancement or damping in a nominally 2D channel. The first, and probably still the only, PIV measurements of flow in a ribbed rotating duct, relevant for internal cooling channels of gas turbine blades, with the equipment rotating together with the duct, were realised by Coletti and Arts [6] and Coletti et al. [7]. The cross-section of the duct is approximately square. There are square ribs on one wall, orthogonal to the main flow direction and parallel to the rotation axis, with rib size to channel height ratio of 0.1. The Reynolds number is 15,000 and the rotation number is 0.3, both based on hydraulic diameter and bulk velocity. Rotations in clockwise and counter-clockwise sense were investigated. With clockwise rotation, Kelvin-Helmholtz rollers are produced in the separated shear layers behind the ribs on the suction side of the channel due to damping of the turbulence by the Coriolis force. This leads to a much longer separation bubble behind a rib than with a stationary channel. With counter-clockwise rotation, a small separation bubble forms at the leading edge of a rib which reattaches very quickly at the top of the rib. Further downstream, a three-dimensional turbulent motion is created due to amplification of the turbulence by the Coriolis force on the pressure side of the channel. This leads to a shorter length of the separation bubble behind a rib compared to the non-rotating case. Coletti et al. [7] recorded mean-flow streamlines in the mid-span plane of the channel and profiles of mean and fluctuating streamwise velocity components on 6 positions in the mid-span plane. They also recorded the streamwise evolution of the mean wall-normal velocity component on a line just above the ribs. The availability of these high-quality data is a motivation for us to use them for verification of the qualities of a hybrid k-ω RANS/LES turbulence model. Narasimhamurthy and Andersson [1] performed direct numerical simulation (DNS) of flow in a spanwise periodic part of a nominally 2D channel with square ribs on corresponding positions on both walls, orthogonal to the main flow direction. They analysed the stationary channel and two rotation speeds. But up to now, there are no DNS of flows in rotating ribbed ducts with side walls, which could be used for verification of turbulence models. Most computational studies employ the steady RANS technique. An example is by Iacovides [8] with a differential stress model (DSM) and a linear eddy viscosity model (EVM) of k-ε type, combined with a one-equation model near the wall. No terms for system rotation are employed in the transport equations of the k-ε model. The analysis is on a streamwise periodic module of a square duct with square orthogonal ribs on leading and trailing walls, once with the ribs on corresponding positions and once with staggered ribs. For the staggered ribs, simulations are compared with experimental results of the ribbed Ubend measured by Iacovides et al. [4]. The conclusion is that the mean longitudinal flow is well predicted by both RANS approaches, but that the Coriolis force driven secondary motion and the turbulent quantities in the flow over the ribs are much better reproduced by the seven-equation DSM model. The strength of the cross flow and the turbulence intensity level determine critically the heat transfer in a rotating ribbed duct. So, for accurate prediction of heat transfer, it is crucial that Coriolis force induced flow phenomena are well captured by the turbulence model. In practical blade cooling applications, the turbulence model is further challenged by significant effects due to rotational buoyancy. The conclusion concerning EVM-type RANS models is similar for more modern variants. For instance, the cases of Iacovides [8] were analysed by Raisee et al. [9] by steady RANS with a cubic nonlinear eddy-viscosity k-ε model (NLEVM) and a linear k-ε model (EVM). The non-linear model is sensitised to streamline curvature and to system rotation. Both models predict the mean longitudinal flow field quite well, but are not very accurate concerning the strength of the secondary flow and the turbulence quantities, although the NLEVM performs better than the EVM. A second example employing the experimental data of Iacovides et al. [4] on the duct with staggered ribs is by Saha and Acharya [10]. They tested the Kato-Launder k-ε model in unsteady form (URANS), which is a model with no specific terms for system rotation. They also performed large eddy simulations (LES) with a dynamic Smagorinsky subgrid model, but for half the rotation number of the experiments and a much lower Reynolds number. Their conclusion is that the URANS technique, in spite of the timedependent formulation, produces only a mild unsteadiness, unless the rotation number is very high. The foregoing studies on RANS turbulence models are all somewhat hindered by lack of precise experimental data, but, nevertheless, the conclusion is that RANS eddy viscosity models, in principle, cannot reproduce accurately the secondary flow and the turbulence intensity in rotating ducts with ribs and, consequently, cannot accurately predict heat transfer.
Abdel-Wahab and Tafti [11] employed LES with the dynamic Smagorinsky subgrid model for flow in a streamwise periodic module of a square rotating duct with square orthogonal ribs at corresponding positions on the leading and trailing walls, for various rotational speeds and bulk Reynolds number Re=20,000. The flow conditions correspond to those of flow in rotating cooling channels in real aero-engines. They found that the peak values of the turbulent kinetic energy differ by a factor of two in the boundary layers on the suction and pressure sides of the channel and that the amplitude of the small-scale fluctuations close to the pressure side changes significantly with the rotational speed. The LES data of Abdel-Wahab and Tafti [11] were used by Viswanathan and Tafti [12,13] as reference data for testing a hybrid RANS/LES turbulence model of DES-type (Detached Eddy Simulation), based on the k-ω SST model. They showed that the DES results are quite close to the LES results, although the grid for DES was a factor of 2 coarser in each direction. They also studied the qualities of the originating k-ω RANS model on the same grid as for DES and concluded that the RANS model seriously underpredicts the turbulence intensity level, especially at the trailing edge side of the channel, i.e. where turbulence is enhanced by the Coriolis force, and, consequently, the heat transfer rate.
In the current work, numerical results by a hybrid k-ω RANS/LES model and the originating RANS k-ω model are compared with the experimental results by Coletti and Arts [6] and Coletti et al. [7]. Coletti et al. [7] observed that the Coriolis force affects the dynamics of both large and small scales, leading to complicated flow patterns with recirculation zones and secondary vortex motions. As discussed above, the applicability of a hybrid RANS/LES technique of DES type to flow in a ribbed rotating duct was already investigated by Viswanathan and Tafti [12,13], but comparison was made with LES results. Since LES involves modelling, we argue that it has added value to compare with high-quality experimental data. We complement the comparison with the experimental results [7] by analysis of the flow features in cross-sections of the channel. Basic secondary flow features can, in principle, already be understood from experimental analysis by PIV, as e.g. by Elfert et al. [2] and from DES simulations, as e.g. by Viswanathan and Tafti [12], but we argue that a systematic discussion helps in understanding the flow features that are important for enhancing heat transfer in ribbed rotating ducts.
The hybrid k-ω RANS/LES model belongs to the class of unified DES-type models [14]. The originating RANS model is the newest version of the k-ω model by Wilcox [15], with the correction for frame rotation by Hellsten [16] applied to the destruction term in the ωequation. From the discussion above, we know that such an eddy-viscosity model cannot accurately predict the effects of the Coriolis force in a rotating duct. Under the assumption of local equilibrium, the eddy viscosity reduces to a Smagorinsky subgrid viscosity with constant C s = 0.1. The functioning of the hybrid model is based on this feature, which ensures that a proper subgrid model is active in LES mode. The hybrid k-ω model is the last version of a model that was gradually developed in earlier work and tested on plane and round impinging jet flows [17][18][19]. First and second phase results were already presented at conferences [20,21]. The present paper concludes the research.

k-ω RANS and Hybrid k-ω RANS/LES Model Equations
For flow of a constant-density fluid in a duct with periodicity in x-direction (direction 1), the continuity and momentum equations, in a frame subjected to rotation in z-direction (direction 3), are where σ ij = 2νS ij and τ ij = 2ν t S ij − 2 3kδ ij are the components of the molecular viscous stress tensor and the modelled turbulent stress tensor, with ν the kinematic molecular viscosity, ν t the turbulent viscosity and S ij = 1 2 ∂U i /∂x j + ∂U j /∂x i − 1 3 (∂U k /∂x k ) δ ij the components of the shear rate tensor. P denotes the effective pressure, which is a part of the pressure combined with contributions from the potential energies associated to gravity and centrifugal force. ∂P ∂x 1 is the mean pressure gradient in the streamwise direction of the periodic domain, used to impose a prescribed mass flow rate. ε ijk is the Levi-Civita symbol and is the angular velocity around the z-axis of the relative frame with respect to an inertial frame. The transport equations of the hybrid k-ω RANS/LES model are where k is the turbulent kinetic energy and ω the specific dissipation rate. The motivation for the formulation of the dissipation term in Eq. 3 is that the dissipation is ε = β * kω = k 3/2 /L t , where the turbulent length scale is L t = k 1/2 /(β * ω). So, it means that in the dissipation term, the turbulent length scale L t is replaced by a length scale proportional to the grid size (C DES ) when the model turns into LES mode. The length scale in the kequation (3), is the maximum size = max( x , y , z ), where x , y , z denote the distances between the centres of the cell faces in x, y and z directions, following Strelets [22]. The closure coefficients and functions of Wilcox [15] are β * = 0.09, α = 0.52, β = β 0 f β , β 0 = 0.0708, σ = 0.5, σ * = 0.6, σ do = 0.125, with * ij = 1 2 ∂U i /∂x j − ∂U j /∂x i − ε ij 3 the components of the rotation rate tensor as seen in the inertial frame.
The eddy viscosity is defined according to Kok et al. [23] by where LES = ( x y z ) 1/3 . In Eq. 6, the RANS eddy viscosity can be expressed by ν t = k/ω = β * k 1/2 L t . It means that also in the eddy viscosity expression, the turbulent length scale L t is replaced by a length scale proportional to the grid size (C DES LES ). The chosen grid size is the cube root of the cell volume, so the typical LES grid size. Both grid size measures, and LES , are multiplied with a tuning constant C DES . The justification for using different grid scales in the k-equation (3) and in the eddy viscosity formula (6) comes from the observation that under local equilibrium (production of k equal to dissipation of k), the eddy viscosity reduces in LES mode to a Smagorinsky subgrid viscosity (S is the magnitude of the shear rate tensor): The role of the term ( / LES ) 1/4 is to increase the eddy viscosity on high aspect ratio cells, with respect to the value obtained by the cube root grid size in all turbulence length scale substitutions. We follow here the model approach by Scotti et al. [24], which improves the predictive qualities of LES on anisotropic grids. In Eq. 7, the Smagorinsky factor C s = (β * ) 3/4 C DES is set to the customary value 0.1, which results into C DES = 0.6086. With a hybrid simulation, the first role of the RANS model is to describe near-wall flows. An eddy viscosity model is sufficient for this purpose. The second role is to generate a subgrid model in the core of the flow, which is equivalent to a typical LES subgrid model. The present eddy viscosity RANS model turns into a Smagorinsky model ( 7) on fine enough grids (see later Fig. 5). So, the eddy viscosity model suffices and it does not pay off to use a complicated RANS model in a hybrid simulation.
The hybrid k-ω RANS/LES model belongs to the class of unified DES models [14]. The term DES refers to the switch from RANS behaviour to LES behaviour by modification of the destruction term in the k-equation in order allow breakdown into small-scale eddies of large-scale structures produced by flow instability in shear layers away from walls. The approach is named Detached Eddy Simulation (DES), according to the principle that large "detached" eddies, in shear layers away from walls, are resolved while small "attached" eddies, in wall bounded shear layers, are modelled. It is a double-substitution type, which means that the grid size measure is introduced in LES mode in both the destruction term of the k-equation and in the eddy viscosity formula. It is called a unified model, which means that there is no a priori definition of interfaces between RANS zones and LES zones, but that these follow from the model itself.
The factor for frame rotation and streamline curvature of Hellsten [16] is applied to the destruction term in the ω-equation (4): where Ri is the Richardson number, C rot is a constant equal to 3.6, * is the magnitude of the rotation rate tensor as seen in the inertial frame andS is the magnitude of the shear rate tensor. The motivation for using a correction factor for system rotation in the destruction term of the ω-equation (4), and not in the production term of the k-equation (3), which is also possible, is that the correction term is only allowed in the RANS zones of the hybrid model. In LES zones the role of such a correction term should be much reduced since most of the flow physics is resolved there. The present hybrid model uses the k-equation for determination of the subgrid turbulent kinetic energy in LES mode. The ω-equation decouples from the k-equation in LES mode, apart from the small influence of the diffusion terms. This procedure ensures that a true Smagorinsky model is recovered in LES mode under local equilibrium conditions. At walls, k and ω are set to where u τ = (τ w /ρ) 1/2 , τ w = μS, S R = (200/k + s ) 2 and k + s is a dimensionless roughness height. The walls are assumed to be hydraulically smooth, which means that the dimensionless roughness height may be at maximum 5 [15]. The precise value is not sensitive. The value is set to 4. Results do not change by lowering the value to 1. Figure 1 shows the computational domain for counter-clockwise rotation, the coordinate system, the rotation axis and the components of the Coriolis force for positive x− and y−velocity components. We remark that in the recirculation zone behind a rib, the Coriolis force points away from the ribbed surface in wall vicinity due to reversed flow direction. The flow is simulated as being periodic on a domain of two streamwise modules of the periodic duct, in order to allow detection of flow structures larger than one module. The other boundaries are walls. The ribs have a square cross-section with H = 8 mm. The size of the computational domain is 160 ×83 ×75 mm (L x ×L y ×L z ) in x, y and z-directions [7].

Computational domain and grids
In Fig. 1, the ribbed duct rotates in counter-clockwise sense around the z-axis with the rotation number Ro = D h /U 0 set to 0.3 (D h is the hydraulic diameter and U 0 is the bulk streamwise velocity, both for a cross-section without ribs). The Reynolds number based on  U 0 and D h is Re = 15,000. The clockwise and counter-clockwise rotations of the experiments are simulated by keeping the coordinate frame and the sense of the angular velocity as in Fig. 1, but by turning the computational domain by 180°around the x-axis in order to obtain the clockwise rotation of the experiments. This way, the pressure side is always at the bottom and the suction side at the top. A list of test cases is given in Table 1.
Block-structured grids have been generated. Figure 2 shows longitudinal and cross sections through a basic grid (a) and a fine grid (b). The basic and fine grids have about 2.1 and 6.1 million (M) cells. The grids have been refined close to the ribbed walls. Table 2   lists the maximal values of the non-dimensional grid sizes, in wall units, on the ribbed surface, side walls and smooth surface opposing the ribbed surface. Also, the maximal values of the ratio of the cell sizes in x, y and z directions to the Kolmogorov length scale η are given. The Kolmogorov length scale, η, is calculated from the dissipation rate, obtained by assuming production equal to dissipation in Eq. 3, following You et al. [25]: η = (ν 3 /ε) 1/4 with ε = ε SGS . Particular attention was paid to obtain a good quality grid in vicinity of the ribbed surface. The values of + max are 13 and 12 on basic and fine grids, respectively. Larger values of + max are used on the side walls and the smooth surface opposing the ribbed surface ( + max = 20-35), but these are typical values used in LES computations. The maximal values of d + are 1.4 and 1.0 near the ribbed surface, and 3.5 and 2.0 at other walls in the basic and fine grid simulations, respectively. Grid densities were chosen based on the LES computations of Abdel-Wahab and Tafti [11], who used a grid with about 2.1M cells for the flow in a rotating ribbed duct at Re = 20,000, but considering only one periodic module. Here, the grid has about 2.1 M cells with a lower Reynolds number (Re = 15,000), but for calculations on two modules. We demonstrate in Section 3.3 that the grid with 2.1M cells is fine enough to resolve the majority of the flow dynamics in the rotating duct.

Algorithmic choices
The computations were performed with the ANSYS FLUENT CFD code, version 14. The unsteady turbulence equations (Eqs. [3][4] were implemented with the user-defined scalar functionality and the source terms were added to the momentum equations (2) to account for the Coriolis force. This way, the grid generation and numerical discretisation methods from the package are used, but we determine ourselves the exact form of the turbulence model equations. The mean pressure gradient ∂P /∂x 1 in Eq. 2 was adjusted at each time step, such that a prescribed mass flow rate was enforced. Due to the splitting of the pressure gradient in a mean component (in the streamwise direction) and a deviating part, the deviating pressure term P becomes periodic. With the hybrid RANS/LES model, the bounded central differencing scheme was applied to the convective terms in the momentum equations and the second order upwind scheme to the convective terms in the kand ω-equations. For RANS, the second order upwind scheme was used for discretisation of the convective terms in all equations. For temporal discretisation, the second-order backward Euler implicit scheme was applied for both the RANS and hybrid RANS/LES model simulations. An implicit time stepping technique was chosen to guarantee stability for a CFL-number larger than unity ( Table 1). The time step was, however, small enough so that the CFL-number was smaller than 2 in the whole flow field, in order to ensure that the dissipation due to the time stepping remained small. The basic simulations were performed with time step t = 6·10 −5 sec (which corresponds to a dimensionless time step size t · U 0 /D h = 2·10 −3 ). This time step is the same for the unsteady RANS simulations and the hybrid model simulations. With the fine grid (H3 and H4-cases in Table 1), the time step was reduced to t = 4·10 −5 sec ( t · U 0 /D h = 1.5·10 −3 ) so that the maximum value of the CFL-number was again less than 2. At each time step, inner iteration steps were applied to lower the scaled residuals of the momentum and the transport equations below 10 −5 . The governing equations were sequentially solved with the pressure-correction PISO method and the pressure staggering option was used for the pressure-velocity coupling.
Time-accurate 3-D RANS simulations and hybrid model simulations were initialised with steady RANS results obtained on a 2-D grid (x − y plane in Fig. 1). The 2D RANS results were copied to the x − y planes of the 3D computational domain. After an initial transient stage of 20 through-flow time periods, statistics of selected quantities were collected for the next 70 periods. A through-flow period is the time needed for fluid particles moving at bulk velocity U 0 to travel the streamwise distance L x . The choice for 20 and 70 periods was made in order to ensure that statistics were sufficiently converged. This was verified by comparing averaged values taken during 40, 50 and 70 periods after an initial phase of 20 periods. The computing times for the unsteady RANS and the hybrid simulations are approximately the same. The numbers of inner iterations necessary to obtain the imposed convergence level in a time step are similar (about 20).  Fig. 3). In wall units the RANS zone extends up to d + =20. With the fine grid (6.1M), the thickness of the RANS zones is smaller (not shown). With clockwise rotation, RANS zones also appear in the central part of the duct, due to almost complete absence of turbulence. This results formally in a very small length scale of the turbulence. With the DES-formulation, such a zone is then interpreted as RANS. But, there is no consequence to it because the eddy viscosity generated by the RANS model is almost zero. The RANS zones in the central part are smaller with the fine grid simulations (not shown). Figure 4 shows instantaneous fields of the ratio of modelled to molecular viscosity, ν t /ν, by the hybrid model. The viscosity ratio is much lower than unity in a large part of the flow field, already on the basic grid. Values just above unity are obtained in small parts of the flow field. So, the contribution of the subgrid model to the turbulent fluctuations is already very small on the basic grid.

Characteristics of the DES model
The role of the subgrid model is further illustrated in Fig. 5 by instantaneous fields of the Smagorinsky factor C s = (β * ) 3/4 C DES , calculated by Eq. 7, from the eddy viscosity produced by the turbulence model on the basic grid. In the LES zones, the averaged value is about 0.1. The maximum value is around 0.2 and the minimum value is around 0.05. These values prove that the turbulence model functions as a proper subgrid model for LES, already on the basic grid. So, the basic grid is fine enough and the time step is small enough for true LES. This can also be concluded from the values of the ratio of the grid size measures to the Kolmogorov length scale in Table 2. According to Pope [26], such ratios should be between 8 and 60, which they satisfy. Figure 6 shows the energy spectra of the resolved x-velocity component at position x/H=1, y/H=1.2 (according to the distances in Fig. 1) in the shear layer behind the rib for clockwise and counter-clockwise rotations using the hybrid RANS/LES model on the basic (H1, H2) and fine (H3, H4) grids. There is no dominant frequency. A similar conclusion was drawn by Coletti and Arts [6] by analysing spectra in the shear layer behind the rib at the same position. With clockwise rotation, a somewhat larger energy level is reproduced in the fine grid simulation in the range 0.3 <fH/U 0 <3. But both the basic and fine grid results are very much comparable at low (fH/U 0 <0.3) and high (fH/U 0 >3) frequencies. The spectra are very much comparable for all frequencies fH/U 0 up to 10 on the basic and fine grids for  counter-clockwise rotation. The correspondence of the spectra means that numerical dissipation at the small-scale level (high frequency) is very small on the basic grid. Further, we remark that a very large range of frequencies is resolved, already on the basic grid. Figure 7a shows the vortex structures near the suction (top) and pressure (bottom) surfaces for clockwise rotation, obtained by hybrid model simulation on the fine grid, visualised by iso-surfaces of the q-criterion (values ±10 5 ), coloured by streamwise velocity. Figure 7b is a visualisation of the hairpin-like vortices (S1, S2 and S3) impacting the front part of the rib. The figure demonstrates that the turbulent flow around the ribs is almost fully resolved. The very low contribution of the subgrid model is also clear from Fig. 4. The resolution on the basic grid is somewhat less (results not shown). Figure 8 shows the vortex structures for the counter-clockwise rotation. The conclusion is the same.   Figure 9 shows profiles of mean and fluctuating streamwise velocity at distance x/H =4 from the rib in the mid-span plane of the duct (see Fig. 1) obtained with the hybrid model on the basic and fine grids. Quantities in Fig. 9 were calculated by time-averaging and by space-averaging over the spanwise distance of 1mm, which is the thickness of the laser sheet in the measurements by Coletti and Arts [6]. The width of the averaging strip is 2 and 3 cells on the basic and fine grids, respectively. Averaging was also done over the first and second module of the periodic duct. The small differences between the results on the two grids show again that the resolution of the basic grid is sufficient and that the time step is small enough for accurate resolution of the flow features. Similar levels of agreement are observed at the other distances from the rib (x/H=0, 2 and 6; results not shown). Figure 10 shows contour plots of the instantaneous streamwise velocity component in the x − y plane at mid-span, obtained by the hybrid model on the basic grid for clockwise (a) and counter-clockwise (b) rotation. A first observation is that the ribs cause obstruction to the flow. The consequence is low streamwise velocity in the vicinity of the ribbed surface and positioning of the core of the flow towards the smooth surface of the duct. A second observation is that there is damping of the turbulence in the shear layers behind the ribs for clockwise rotation and amplification for counter-clockwise rotation. With clockwise rotation, the result is a low turbulence level all over the duct. With counter-clockwise rotation, the turbulence from the shear layers spreads over the whole duct width. The spreading is caused by intense cross-flow, discussed in a later section. These observations are confirmed by the profiles at x/H =4 behind the ribs, shown in the earlier Fig. 9. The difference in turbulence level and consequently in turbulent mixing, results in a much broader low velocity zone with clockwise rotation.  Figure 11 shows streamlines of the mean flow in the mid-span plane in the vicinity of a rib, obtained by the hybrid model. The observation is that there is a much larger recirculation zone behind a rib with clockwise rotation. The recirculation zones resemble very well the flow patterns measured by Coletti et al. [7]. In particular, the small recirculation zones in the corners of the ribs and the walls are accurately represented, as well as the separation zone on the upper rib surface for counter-clockwise rotation [27]. The main recirculation zone behind the ribs is accurately obtained for counter-clockwise rotation on both basic and fine grids ( Table 3). The difference between measured (x/H =1.53) and computed (x/H =1.45 and 1.40) streamwise location of the centre of primarily recirculation bubble is less than 9 %. The differences are less than 5 % in transverse y-direction. With clockwise rotation, the numerically obtained rear part of the primary recirculation zone accords very well with the experiments. The front zone differs somewhat. In the corner behind a rib, the numerically obtained intensity of the flow recirculation is slightly less than in the experiments. The consequence is a difference in steamwise location of the centre of the main recirculation   shear stress. Improved results are obtained on the fine grid. There, the position of the vortex centre becomes x/H =2.3 (not shown). The relative error is reduced to 18 %. Figure 12 shows that the large recirculation patterns in the longitudinal flow are also quite well reproduced by the time-accurate RANS, although much less accurately than with the hybrid model. The secondary features are not fully captured. With RANS, the centres of the primary recirculation zones behind the ribs are obtained at the streamwise distances x/H =3.50 (relative error 80 %) and x/H =2.30 (relative error 50 %) for clockwise and counter-clockwise rotations (Table 3). With clockwise rotation, the difference may be explained by the too low turbulence level in the shear layer behind the rib (see Fig. 13a in the next subsection and Fig. 26a in the later section on turbulent shear stress), resulting in too low mixing in the front part of the recirculation zone. With counter-clockwise rotation the difference may be explained by a too large stress level produced in the accelerating zone in front of the ribs (see the later Fig. 26b), leading to entrainment of the rear part of the separation bubble. Despite a significant difference in the centre location, the overall size of the recirculation zones is only slightly larger than by the hybrid model.

Resolved and modelled fluctuations
Some understanding of the different behaviour of the k-ω RANS and the hybrid k-ω models may be gained from the difference in resolved and modelled fluctuations. Figure 13 shows  profiles of resolved and total fluctuating streamwise velocity near the ribbed wall at distance x/H =4 from a rib, obtained on the basic grid. With the RANS model (R1), the fraction of the resolved turbulence is 60 % to 90 % for the clockwise rotation (Fig. 13a). The level of total fluctuating velocity is too low compared to the measurements. However, as shown in Fig. 14 in the next subsection, this does not have a strong negative effect on the mean velocity profiles which are close to the measured profiles and the profiles obtained by the hybrid model. With the hybrid model (H1), the fraction of resolved turbulence is very high, about 80 % (Fig. 13a). The 80 % level is normally used as a criterion for fully resolved LES [26]. With counter-clockwise rotation (Fig. 13b), the total fluctuating velocity values by RANS and hybrid RANS/LES are comparable and near the experimental level, but the fractions of resolved and modelled turbulence are much different. With the RANS model, a very high fraction of the turbulence is modelled (only 5-10 % is resolved) on the pressure side of the duct. The correct turbulence level by RANS comes from the correction term for system rotation by Hellsten [16], who calibrated this term for flows in straight rotating channels. With the hybrid model (H2), similarly as for clockwise rotation, the fraction of resolved turbulence is about 80 %, which means that reliable LES is reached in the LES zones. The results also show that the role of the correction term for system rotation is secondary for the hybrid model, since most of the turbulence is resolved.

Profiles of mean and fluctuating velocity
Profiles of mean and total (resolved plus modelled) fluctuating streamwise velocity by the RANS and hybrid RANS/LES models are shown in Figs. 14 and 15, along lines perpendicular to the ribbed surface at mid-span at streamwise distances x/H = 0, 2, 4 and 6 behind the ribs. The results are obtained on the basic grid with 2.1 million cells. Similarly as before, the flow statistics are determined by time-averaging and by space-averaging over a strip of 1 mm. Figure 14 shows that there are no big differences between mean velocity profiles obtained by the RANS and hybrid models and the experimental ones for clockwise rotation. We note a too weak flow reversal at x/H =2 both by the RANS and hybrid models (Fig. 14b). Too low values of the fluctuating streamwise velocity component are obtained by RANS  Fig. 15), while the fluctuating velocities are accurately predicted by the hybrid model. The differences between the results of the RANS and the hybrid models are larger for counterclockwise rotation. The time-accurate RANS produces a too large velocity gradient at the walls just after the reattachment (Fig. 14c and d). This is mainly due to the too large turbulent shear stress predicted by RANS in between the ribs (see the later Fig. 26b). Results by the hybrid RANS/LES model are much better for mean and fluctuating velocity components (Figs. 14 and 15). The peak value of the fluctuating velocity component (Fig. 15a and b) is much better reproduced by the hybrid model at y/H =1 than by RANS. This is due to the ability of the hybrid model to resolve the streamwise flow unsteadiness there. Figure 16 shows the streamwise evolution of the mean wall-normal velocity component at distance y/H =0.1 from the top of the ribs, which is at the height of the separated shear layer behind a rib, for clockwise and counter-clockwise rotation. The positive V /U 0 in Fig. 16a the negative V /U 0 in Fig. 16b denote the y-velocity component directed towards the ribbed wall. Overall, the hybrid model results agree better with the experiments than the RANS model results. Both the hybrid and RANS models produce a somewhat too small variation of the y-velocity component over the primary recirculation bubble at x/H <5, for clockwise rotation (Fig. 16a). This shows that both the RANS and hybrid models have difficulties in reproducing the weak cross-flow unsteadiness in the separated shear layer  (Figs. 11a and 12a) this results in a somewhat too long size of the separation bubble returned by both the RANS and hybrid models. For counterclockwise rotation (Fig. 16b), the RANS model misses the peak value around x/H =0.4. This is due to difficulties by RANS in reproducing the separation bubble on the top of a rib owing to a too large turbulent shear stress in the accelerating flow in front of a rib (see the later Fig. 26b). Figures 17 and 18 show the mean and fluctuating velocities obtained by the hybrid model on the basic and fine grids. The fine grid results do not differ much from the basic grid results. It shows, again, that sufficient resolution is obtained on the basic grid. The agreement between measurements and simulations is particularly good close to the wall at distances x/H =4 and 6. The mean flow characteristics are somewhat less accurately reproduced at x/H =0 and 2 for the clockwise rotation ( Fig. 17a and b), but this has no big effect on the results downstream ( Fig. 17c and d).

Conclusion concerning the longitudinal flow field
The hybrid k-ω RANS/LES model produces very realistic mean and fluctuating longitudinal flow velocities. The time-accurate k-ω RANS model produces also quite good results   for the time-averaged longitudinal flow, if we accept that RANS produces a somewhat too strong flow reattachment for the counter-clockwise rotation (Fig. 14b and c). For clockwise rotation, the k-ω RANS model results for the mean longitudinal flow are very similar to the results of the hybrid model and they are close to experiments. The prediction of the fluctuating velocity in the longitudinal flow is much less accurate by the k-ω RANS model than by the hybrid k-ω RANS/LES model, but the predicted fluctuating velocity levels are not fully erroneous. The better results of the hybrid k-ω RANS/LES model are clearly due to the much higher level of resolved turbulence.

Hybrid RANS/LES results
Transversal flow patterns are shown in Fig. 19 by mean velocity vectors projected in z-y planes at x/H =0 and 4, obtained with the hybrid model for clockwise and counterclockwise rotations. The velocity vectors are displayed in each 4-th node. For clarity, the reference velocity magnitude is the bulk velocity U 0 divided by 80 for clockwise rotation (Fig. 19a) and divided by 10 for counter-clockwise rotation (Fig. 19b). Much stronger cross-flow is obtained for counter-clockwise rotation. The right panels of Fig. 19 show a magnified view of the flow near the right side wall (at x/H =4). Strong flow towards the suction side is observed in the boundary layers at the side walls. The boundary layer thickness is about z/H =0.4 and z/H =1.2 at y/H =5 for clockwise and counter-clockwise rotations, respectively. So the side-wall boundary layers are much thicker for counterclockwise rotation. The main reason for the thicker boundary layers with counter-clockwise rotation is the much higher turbulence level in the whole duct. A second reason is the superposition of two phenomena that feed transversal flow into the side-wall boundary layers. First, there is the obstruction effect of the ribs. Because the streamwise velocity magnitude is much lower in the side-wall boundary layers than in the centre of the duct, but pressure does not change much in spanwise direction in the approximately parallel flow approaching a rib, the stagnation pressure in the centre of the duct is larger than at the side-walls. The consequence is flow migration above a rib towards the side-walls and formation of streamwise vortices similar to a leg of a horseshoe vortex at the base of a bluff body placed on a flat plate. The result is two vortices above a rib and downstream of a rib near the side-walls with rotation sense corresponding with transversal flow in the side-wall boundary layers away from the rib. These vortices are magnified by the movement towards the ribbed wall in the centre of the duct downstream of a rib. The generation of rib-induced vortices was clearly demonstrated by PIV measurements in cross-sections of a stationary duct with square orthogonal ribs on one side by Casarsa and Arts [28] and by LES of the flow in the same ribbed duct by Lohasz et al. [29]. The second phenomenon is the Coriolis force induced transversal flow in side-wall boundary layers from the pressure side wall to the suction side wall. With counter-clockwise rotation, the rib-induced and Coriolis-induced transversal flows are in the same sense. This explains the strong transversal flow in the side-wall boundary layers. The consequence is two clearly formed streamwise vortices in the corners of the side-walls and the suction-side wall (Fig. 19b). With clockwise rotation, the ribinduced and Coriolis-induced transversal flows have opposite senses. This leads to a much weaker transversal flow from the pressure-side wall to the suction-side wall in the side-wall boundary layers. The consequence is a much lower velocity level of the cross-flow and more complex streamwise vortex patterns. Two rib-induced corner vortices can be identified in Fig. 19a near the rib, two Coriolis-induced vortices near the side walls and four streamwise vortices at the pressure-side wall. Two of these are generated by the flow movement in the core towards the pressure side and two are Coriolis force induced in the corners with the side-walls.
With clockwise rotation (Fig. 19a), fluid elements are attracted towards the suction side for y/H >8. This effect is caused by the Coriolis force pointing towards the wall in near-wall vicinity inside the large separation zone behind the ribs (negative streamwise velocity resulting in positive Coriolis force component 2 U). With counter-clockwise rotation (Fig. 19b), the centres of the Coriolis-induced vortices are near the suction side (y/H =8.5) and much larger vortex cells are produced with respect to clockwise rotation. At x/H =4, a strong flow towards the pressure side downstream of the recirculation zone behind the rib is noticeable with counter-clockwise rotation. The smaller secondary flow vortices produced at the pressure side with clockwise rotation do not appear with counter-clockwise rotation. This is due to strong turbulence in the shear layer behind the rib which leads to damping of the secondary flow motion.
Further insight may be gained by analysis of instantaneous cross-flow fields, as shown in Fig. 20. The much more intense cross-flow with counter-clockwise rotation is clear in the plots of instantaneous velocity vectors in z − y planes shown together with contour plots of streamwise velocity. With clockwise rotation, the unsteady flow motion on the suction side of the duct is separated from the perturbed flow motion on the pressure side. The transverse flow is almost unperturbed between y/H =4 and 8. There is some reduction of instantaneous streamwise velocity at the side walls and between the vortex cells on the pressure side. This reduction results in local movement of fluid elements towards the suction side in the side-wall boundary layers and between the vortex cells due to local reduction of the Coriolis force there (-2 U). In the time-averaged secondary flow, the latter effect is very weak, however. This explains that the Coriolis force induced vortex cells produced near the pressure side remain in the vicinity of the pressure side (Fig. 19a). With clockwise rotation, vortex cells are visible near the suction side in the instantaneous velocity field (Fig. 20a), but they do not lead to clear vortices in the time-averaged velocity field (Fig. 19a). Apparently, they are mixed out by the weak turbulence in the shear layers behind the ribs. The cross-flow is much more intense with counter-clockwise rotation (Fig. 20b). Strong and large vortex structures are produced in the shear layer behind the rib on the pressure side. The flow is highly three-dimensional. There are large zones with low instantaneous streamwise velocity. In these zones, fluid particles are pushed instantaneously towards the suction side (reduced Coriolis force). Conversely, fluid elements with higher streamwise velocity are pushed towards the pressure side. The vortex structures perturb the side-wall boundary layers already near the suction side of the duct and their activity increases strongly towards the pressure side due to their increased size. So, strong instantaneous transversal velocities inside the whole duct and thick side-wall boundary layers are caused with counter-clockwise rotation by the large highly time-dependent vortex structures in the shear layer behind a rib on the pressure side of the duct. Small-scale turbulent motions result by breakdown of the large vortices and spread over the full duct section. This confirms the observation by Coletti et al. [7] that the Coriolis force affects both the large-and small-scale dynamics. The strong mixing at the pressure side with counter-clockwise rotation explains why no vortices occur in the time-averaged cross-flow velocity field near the pressure side and the main vortices are pushed towards the suction side (Fig. 19b). Figure 21 shows the time-averaged velocity vectors in the cross-sections x/H =0 and 4, obtained by the RANS model. The reference velocity magnitude is the bulk velocity U 0 divided by 80. Figure 21a shows that the main cross-flow patterns of the RANS and hybrid models are similar for clockwise rotation. It means that the RANS model, similarly to the hybrid model, is successful in representing the turbulence suppression on the suction side for clockwise rotation. There is some difference in the vortex patterns on the pressure side, but, as we showed earlier by Fig. 14, this does not affect the mean velocity profiles near the ribbed surface which remain close to the profiles obtained by the hybrid model.

RANS results
The counter-clockwise rotation is much more challenging for RANS (Fig. 21b). In this case, the RANS model does not reproduce the secondary flow patterns near the side walls. As we demonstrate by the later Fig. 26b, this is due to overprediction of turbulent shear stress by the RANS model on the pressure side of the duct. The consequence is that the RANS model does not produce noticeable flow unsteadiness in the shear layer behind the rib and that the flow is largely two-dimensional. The almost complete absence of resolved fluctuations is demonstrated in Fig. 22b. The flow unsteadiness is also very small for clockwise rotation (Fig. 22a). As a result, the side-wall boundary layers stay very thin for counter-clockwise rotation due to lack of perturbations by three-dimensional unsteadiness. The thickness of the boundary layers is even smaller than that with clockwise rotation. The flow inside the boundary layers also experiences much smaller acceleration. Overall, a much weaker cross-flow is obtained by RANS for counter-clockwise rotation (take into account the different magnitude of the reference velocity in Figs. 19b and 21b). The above results indicate that the hybrid model is much more successful in reproducing the secondary flow.

Conclusion concerning secondary flow
The hybrid k-ω RANS/LES model produces a strongly time-dependent secondary flow. Although we do not have experimental results on the secondary flow, based on the good representation of fluctuating velocity components in the longitudinal flow, discussed in the previous section, we can deduce that the representation of the secondary flow by the hybrid model is basically correct. The k-ω RANS model produces a nearly steady secondary flow, which cannot be correct. 6 Turbulent Shear Stress

Hybrid RANS/LES results
The turbulent shear stress is analysed in order to further illustrate the effect of the Coriolis force on the turbulent flow dynamics. Figure 23 shows contours of total Reynolds shear stress (τ Bouss -<u'v'>)/U 2 0 (modelled plus resolved) in the x − y plane at z/H =0 for clockwise and counter-clockwise rotations, obtained by the hybrid model on the basic grid (take into account the much different scales in Fig. 23a and b). The modelled shear stress τ Bouss is approximated with the Boussinesq assumption, τ Bouss =2<ν t S xy >, where S xy =1/2(∂u/∂y+∂v/∂x). Turbulence is suppressed in the separated boundary layers on the suction side for clockwise rotation (Fig. 23a) and enhanced in the separated boundary layers on the pressure side for counter-clockwise rotation (Fig. 23b). With clockwise rotation, the turbulent shear stress just behind a rib is somewhat too small with respect to the measurements [7]. In the experiments, the peak Reynolds shear stress level of -<u'v'>/U 2 0 = -0.0040 is also obtained immediately behind the ribs. The difference in shear stress means that computationally not enough turbulence is generated just behind the ribs. It explains the difference in the measured and computed centres of the main recirculation zone behind a rib (Fig. 11a). With counter-clockwise rotation, the predicted and measured turbulent shear stresses agree well with each other. The foregoing global comparison is visual, based on the contour plots published by Coletti et al. [7].
A quantitative comparison of the turbulent shear stress is presented in Fig. 24 by profiles at the positions x/H =0, 2, 4 and 6 in the mid-span plane for counter-clockwise rotation (the experimental data were privately provided by Dr. F. Coletti). The correspondence between  Figure 25 shows contour plots of the total turbulent shear stress (τ Bouss -<u'v'>)/U 2 0 in cross-sections at three streamwise distances x/H =0, 4 and 8. A very low shear stress level is obtained with clockwise rotation on both suction and pressure sides of the duct in all planes. Two peaks of shear stress are visible on the pressure side (Fig. 25a). They are induced by the streamwise-oriented vortices near the pressure side (see Fig. 19a). The turbulence is transported away from the pressure side at z/H = -2 and 2 into the side-wall boundary layers. There is only weak turbulent diffusion in these layers. A much higher level of turbulent shear stress is produced on the pressure side for counter-clockwise rotation (take into account the much different scales in Fig. 25a and b). The turbulence produced in the shear layer behind a rib on the pressure side is transported over the full duct section. This transport is caused by the strong Coriolis-induced streamwise-oriented vortices (Fig. 20b). Strong turbulent diffusion is observed in the side-wall boundary layers.
The above results show that the Coriolis-induced secondary flow features have a strong effect on the turbulent flow dynamics inside the duct, for counter-clockwise rotation. With clockwise rotation, the turbulent shear stress produced beneath the roll cells at the pressure side is weak and it is only weakly transported towards the centre of the duct. A much higher level of turbulent shear stress is produced on the pressure side with counter-clockwise rotation due to the presence of the ribs. This stress is also much further transported towards the centre. Figure 26 shows contour plots of the turbulent shear stress in the mid-plane obtained by the RANS model. With clockwise rotation, the turbulent shear stress level on the pressure side (Fig. 26a) is similar to that of the hybrid model (Fig. 23a). Too low turbulent shear stress is produced on the suction side. It means that the correction term for system rotation overestimates somewhat the damping of turbulence on the suction side for clockwise rotation. This explains, therefore, the somewhat too large size of the primary recirculation zone obtained with the RANS model (Fig. 12a). A much larger level of turbulent shear stress is produced on the pressure side with counter-clockwise rotation (Fig. 26b). The zone with extremely large shear stress is not in agreement with the measurements [7]. The peak value of the turbulent shear stress should be immediately behind the ribs, similarly to what is obtained by the hybrid model, but not in the accelerating flow region between ribs. This means that the k-ω RANS model largely overpredicts the turbulent shear stress on the pressure side of the duct. The large level is caused by reduced turbulent dissipation (ω) on the pressure side due to the correction term for system rotation. So, it means that the correction term overreacts in the zone of accelerating flow. The correction term is meant to bring the level of fluctuating velocity to the correct value at the pressure side (Figs. 13b and 15), but with an eddy-viscosity model the unavoidable consequence is a large turbulent shear stress. Obviously, this large turbulent shear stress is not physically correct. But the zones of high turbulent shear stress by RANS are rather far away from the primary recirculation zones behind the ribs. This limits the effect on the mean and fluctuating velocity profiles in the zone behind the ribs, which are quite well captured by RANS (Figs. 14 and 15).

Conclusion concerning turbulent shear stress
The hybrid k-ω RANS/LES model reproduces quite accurately the turbulent shear stresses of the experiments, while the turbulent shear stresses by the k-ω RANS model are very erroneous. In particular, for counter-clockwise rotation, the k-ω RANS model produces a much too high turbulent shear stress level at the pressure side of the duct in the accelerating flow approaching the ribs. The correction term for system rotation, which is meant to bring the modelled turbulence intensity to the correct level, has this overprediction as a side effect. The result is that the flow upstream of the ribs is represented as almost steady with RANS, both for clockwise and counter-clockwise rotation. This explains the almost complete absence of unsteadiness in the secondary flow obtained by k-ω RANS (Fig. 22).

Conclusions
The hybrid k-ω RANS/LES model produces a very good simulation of the flow through a ribbed duct subjected to system rotation. The hybrid model allows for successful prediction of primary and secondary flow details, both for mean velocity and for fluctuating velocity. With the hybrid model most of the turbulence is resolved, so the role of the underlying turbulence model and the correction term for system rotation is secondary. Addition of the correction term for system rotation, which was originally tuned for smooth rotating ducts, allows obtaining realistic mean flow characteristics in the mid-span plane by the timeaccurate k-ω RANS technique. Also the fluctuating velocity components in the longitudinal flow are quite realistic by the k-ω RANS model. But the k-ω RANS model is strongly erroneous on the Coriolis-induced secondary flow features for counter-clockwise rotation. This means that the k-ω RANS model is unreliable for heat transfer prediction. The correction factor for system rotation does not help for this aspect. It brings the level of the fluctuating velocity to the correct value for the longitudinal flow, but a side effect is that it damps almost all fluctuations in the secondary flow near ribs where turbulence is enhanced by the Coriolis force. The strong quality difference between the results for secondary flow demonstrates the large benefit by turning the time-accurate k-ω RANS model into a hybrid k-ω RANS/LES model. With the same space and time resolutions, this change is done without any supplementary need for computer memory and computing time. components of shear rate tensor S ij = 1 2 ∂U i /∂x j + ∂U j /∂x i − 1 3 (∂U k /∂x k ) δ ij (s −1 ) grid size based on distances between centres of cell faces: max( x , y , z ) LES grid size based on cube-root of the cell volume: ( x y z ) 1/3

Nomenclature
x , y , z distances between centres of cell faces in x, y and z directions ε ij k Levi-Civita symbol ν t turbulent viscosity (m 2 /s) ω specific dissipation (s −1 ) angular velocity of the duct with respect to an inertial frame (rad/s) * magnitude of rotation rate tensor in the inertial frame: 2 * ij * ij (s −1 ) * ij components of vorticity tensor as seen in the inertial frame * ij = 1 2 ∂U i /∂x j − ∂U j /∂x i − ε ij 3