Numerical modelling of turbulent transition in complex geometries

The work deals with numerical simulation of laminar-turbulent transition in transonic flows in turbine cascades. The 3D cascade geometry as well as 2D model cascade in a wind tunnel is simulated. The γ-ζ transition model is based on empirical criteria for start of the transition. The implementation of the model is discussed including reformulation of the criterion for transition on separation bubble.


Introduction
Mathematical models based on the Favre-averaged Navier-Stokes equations are well accepted tools in predicting turbulent flows, however, the accuracy is influenced by turbulence model and often by its capability to predict transition to turbulence. Common two-equation eddy-viscosity turbulence models usually predict too early start of transition and then the transition is too fast. This problem is further emphasized by over-prediction of the turbulent energy production on the leading edge of the blade which has its origin in the eddy-viscosity assumption. Some ad hoc remedies of the later problem has been proposed e.g. limiter by Medic, Durbin [1]. However the transition still requires explicit triggering. Considering transition models based on transport equations which seem more general than algebraic ones, the models contain equation for an intermittency variable and also for other auxiliary variable or variables. More recent examples are 3-equation model by Walters and Cokljat [2] or 2-equation model by Menter, Langtry [3]. Also 1-equation model is proposed by Durbin [4]. These models have "local" form enabling easy implementation especially on parallel computers. However they also share disadvantage of containing transition criteria implicitly, or as Menter, Langtry model require evaluating criteria in every point of the flow-field which is time consuming and ambiguous especially when moving walls are present. Although there are new attempts to improve upon these models, their reliability remains to be confirmed. In this work we apply a more conservative γ-ζ model of Lodefier, Dick [5] and Kubacki et al [6] instead. The model contains transition criteria explicitly and any new criterion can be added easily. The empirical transition criteria are reliable. The downside is that the model distinguishes free-stream and boundary layer and thus is not local. Nevertheless we show that also for non-trivial geometry when using multi-block grids the model can be blocka Corresponding autor: petr.louda@fs.cvut.cz local and does not require case-specific input under assumption that the whole thickness of boundary layer is contained in one block, at least in region where transition occurs. The model is applied here in domains with non-trivial boundary as 2D cascade in wind tunnel or 3D periodic cascade The implementation does not require any case-specific user input.

Mathematical model
The mathematical model of turbulent flow is based on Favreaveraged Navier-Stokes (NS) equations, see e.g. Wilcox [7]. The system consisting of continuity, 3 momentum and energy equations can be written in 3D in Cartesian coordinates as where V is control volume, n i outer unit normal vector components of its surface, t time, ρ density, u i velocity vector components, E total energy per unit volume, H = E + p/ρ is total enthalpy and p static pressure. The normal velocity u c = u i n i . Summation convention is used for repeated indices. Equation of state for perfect gas is prescribed in the form with the ratio of specific heats γ = 1.4 and k being turbulent energy. The molecular stress tensor and heat flux vector respectively are assumed in the form where δ i j is the Kronecker delta. The dynamic viscosity µ and Prandtl number satisfy The effect of turbulent fluctuations is present by the Reynolds stress tensor τ i j and turbulent heat flux q t i , which need to be modeled. An eddy viscosity model assumes where µ t is eddy viscosity and the turbulent Prandtl number is set Pr t = 0.91. In a k-ω model, the eddy viscosity µ t ∼ ρk/ω. The k-ω system can be written where the turbulent production P k = τ i j ∂u i ∂x j , the α, β, β * , σ k , σ ω are model coefficients and CD a cross-diffusion term, CD ∼ (∂k/∂x i )(∂ω/∂x i ). In this work the SST variant of k-ω model is use, see Menter [8].

Transition modelling
The two-equation turbulence models for RANS alone do not predict the transition to turbulence with satisfactory precision. Therefore some intermittency function multiplying turbulent stresses or merely turbulence production is introduced. In this work the γ-zeta intermittency model of Dick et al. [9] is chosen. It has been developed with emphasis on transition phenomena occurring in turbomachinery as is the bypass transition, transition on separation bubble or wake induced transition. For better prediction of the later case the intermittency is split into near-wall part γ and free-stream part ζ each governed by its transport equation with source terms where µ t is eddy viscosity, µ ζ viscosity for ζ (see [9]), β γ growth parameter, U γ mean velocity scale, n wall-normal direction. The Prandtl numbers σ γ = σ ζ = 1. The F s is switching function equal either 0 in laminar boundary layer or 1 if a criterion for turbulent flow is satisfied. The transition criteria can be chosen according to expected transition mechanism. The source term S ζ guarantees that ζ vanishes in the boundary layer, so that both intermittencies are complementary in the turbulent flow. Since the source terms are not directly linked, the total intermittency (γ + ζ) has to be limited to < 0, 1 > explicitly. The boundary conditions are γ = 0, ζ = 1 in the inlet, ∂γ/∂n = ∂ζ/∂n = 0 in the outlet and ∂γ/∂n = 0, ζ = 0 on the wall. The γ-ζ transition model relies on meeting various transition criteria in switching source term for near-wall intermittency γ. For bypass transition the Mayle and Abu-Ghannam, Shaw criteria are used, see [10,6]. The criteria are functions of boundary layer parameters as momentum thickness,freestream velocity and turbulence intensity. These have to be computed preferably in an automatic manner. The boundary layer edge detection algorhitm is taken from the original article [9]. The boundary layer edge is defined as distance from wall where the vorticity reaches 1% of its maximum, found on the wall in the attached flow. In separated flow the local maximum farthest from the wall is taken. This thickness is further increased by 30 % to reach the free-stream.
The criteria for transition on separation bubble applicable in internal aerodynamics moreover depend explicitly on the upstream history of the flow. Here the Mayle criterion in the form is used, where Re st is the Reynolds number from the distance from the point of separation x s , Re Θs is the Reynolds number from momentum thickness at the point of separation and T u e free-stream turbulence intensity. Obviously, the tracking of separation is very cumbersome in 2-dimensional case (although it has also been done by the authors) and plain ambiguous and thus impossible in 3D case. Therefore the following re-formulation is proposed, based on the observation of the simulated separation bubbles. To an acceptable degree of precision, one can assume that the separation thickness δ 0 starts growing linearly with th distance from separation x s . Then where δ 0 is the distance of zero velocity streamline from the wall and the slope in streamwise direction ∂δ 0 ∂s is evaluated on standard stencil parallel to wall. The non-local formulation persists only in the wall-normal direction which is acceptable. The empirical coefficient A = 0.9 is not very critical. If the slope ∂δ 0 ∂s is negative, the maximum thickness of the bubble has been reached. It is assumed that in this part the transition started at latest. In the attached state the criterion for attached boundary layer takes over. Further step to application of the transition model in separated flow is necessary. Since the source term in γ-equation is always non-negative, once the intermittency is captured in the recirculation it can not diminish except by diffusion which, however, has mostly opposite effect. The separation is then predicted as too short. Therefore a numerical sink term is switched-on inside the recirculation bubble.
In 3D cases a flow parallel to a corner can be present. The vorticity magnitude here has no more maxima on the wall. Therefore the boundary layers on both walls are observed and if they are no more well defined, the average of free-stream parameters from both walls are used for the whole corner area. The intersection of 3 walls is excluded, not least because of the lack of physically valid transition criteria. Such corner in principle can be dealt with by placing arbitrarily small grid block into the corner and setting the intermittency production there ad hoc, e.g. to zero.

Numerical method
The above mathematical model is solved numerically by an in-house Fortran 90 code implementing implicit finite volume upwind solver. For spatial discretization we use a cell centered finite volume method with quadrilateral (in 2D) or hexahedral (in 3D) finite volumes composing multiple blocks of structured grid. The numerical inviscid flux is computed by the AUSMPW+ splitting [11]. The higher order of accuracy is achieved by linear interpolation in the direction of grid lines with e.g. van Leer limiter. The discretization of diffusive flux is central. The approximation of cell face derivatives needed in diffusive terms uses octahedral dual finite volumes constructed over each face of primary volume -the vertices are located in vertices of primary face and in centers of adjacent primary volumes. For time discretization, the implicit backward Euler scheme is employed where the steady residual at new time level is approximated by linear extrapolation. The Jacobi matrices of the flux are obtained as derivatives of discrete expressions for flux with respect to nodal values from the stencil of implicit operator. We chose 7-point stencil, which leads to block 7-diagonal system of linear equations (not considering boundary conditions). The size of a block equals to the number of coupled equations. For more details see [12]. The turbulence model equations are solved in the same manner, however, decoupled from the RANS equations. The decoupling saves memory and computational time since there is no need for matrix-vector multiplication if the source terms Jacobi matrix is chosen diagonal.
Some remarks concerning transition model have already been presented. In parallel implementation it is natural to distribute the work block-wise, although it cannot be optimal in all cases. Therefore we require the model be at least blocklocal. This is satisfied if the whole thickness of a boundary layer is contained in same block. It can be satisfied by sufficiently thick O-grid around the blade.

Computational results
In this part we show application of the model to 2D and 3D flow through turbine cascades. In 2D we consider flow through model cascade consisting of six mid-section blades in a wind tunnel. The boundary approximates wind tunnel facility of the Institute of Thermomechanics AS CR and is shown in Fig. 1. To improve periodicity of the flow around cascade, internal wall is located behind top blade. In measurements, carried out by Luxa et al. [13,14] the internal wall is perforated plate with staggered rows of circular holes. Obviously, in 2D a model is needed. In present simulations it is simplified to nonpermeable line of zero thickness. The isentropic outlet Mach number of the tunnel is 1.120, the isentropic Reynolds number for the blade 1.2 · 10 6 .
The computational grid consist of approx. 310 000 finite volumes divided into 192 blocks. Each blade is discretised by 260 grid points. We note that the model of transition does not require any case specific user input other than merely activating the model. The Fig. 2 shows simulations of the experimental cascade in the wind tunnel, the Fig. 3 same cases with internal wall. The internal wall greatly improves the periodicity of the flow through cascade, although only impermeable internal wall has been applied. It can be expected that with perforated wall the flow pattern will be similar but with weaker shock wave reflexions. The transition causes maximum Mach numbers little bit higher than in fully turbulent case. Both cases with internal wall are very similar except that with transition model the boundary layer on the suction sides is laminar up to the interaction with the incident shock wave. In more detail, the middle channel is shown in Figs. 4 and 5. Fully turbulent simulation shows simple reflection of the shock wave from turbulent boundary layer. With the model of transition the boundary layer in front of the interaction is still laminar and small separation exists under the shock wave. The deformation of boundary layer causes formation of separation shock at beginning of the recirculation bubble and reattachment shock at its end. In between expansion waves are present. A better resolution of the interaction is possible on finer grid but for the purpose of the whole tunnel study we stop here. The corresponding near wall intermittency γ is shown in Fig. 6, where red color denotes turbulent boundary layer (γ = 1).
Considering comparison with experiment finding corresponding regimes appears difficult. The reference is isentropic Mach number in a certain point in the settling chamber. This is different for each tunnel simulation although the outlet isentropic Mach number is same (M 2is = 1.120) and also the mass flow is practically same (44 kg/s/m). The reference Mach number varied between 1.244 and 1.320 in mentioned 4 simulations. Therefore we show two measured regimes in Fig. 7 for two Mach numbers. The configuration is without internal wall and should correspond to simulation from Fig. 4. The computed shock wave-boundary layer interaction is correct but the angle of the shock wave seems to correspond to higher reference Mach number.
As an example of 3D simulation, the SE1050 mid-section turbine cascade is chosen. The outlet Mach number M 2is = 1.2 and Re 2is = 1.5 · 10 6 . One period of the cascade is solved and the blade is prismatic with symmetry condition in the midplane. The blade has 230 × 95 grid points on its surface. The inlet angle is design angle. A small recirculation bubble exists where incident shock reaches suction side. This is shown on streamwise velocity component in blue in Fig. 8 whereas the interaction with incident shock is turbulent as shown on the isolines of Mach number in the same figure. Indeed the boundary layer on the largest part of the suction side is turbulent as shown on the near wall intermittency in Fig. 9. The boundary layer on the pressure side is laminar.

Conclusions
The work presented simulations of 2D and 3D turbulent flow through turbine cascades with eddy-viscosity turbulence model complemented with the γ-ζ model of transition to turbulence. The mathematical model is solved by implicit AUSM finite volume method on multi-block structured grids. The implementation of transition model does not rely on explicit prescription of boundary layer edge and is adaptive as long as the whole thickness of boundary layer is contained in one block, which is typically O-grid around the blade consisting of several blocks in tangential direction. Also the treatment of corners is automatic. The transition mechanisms include bypass transition and transition on the separation bubble. The transition criterion for the later mechanism has been re-formulated in order to be local in streamwise direction. The physical correctness of transition prediction in the flow in convex corner however still needs to be confirmed by measurement. The results are shown for 2D turbine cascade in wind tunnel where the transition model enables to capture measured interaction of shock wave with laminar boundary layer. Further the model has been applied to 3D periodic turbine cascade.