Numerical simulation of high-swirl flow in axial turbine stage

This article deals with numerical modelling of flow in an axial turbine stage with prismatic blades which are not equipped with shroud. In this case flow through a radial gap above and below ends of blades leads to generation of strong secondary vortices which significantly affect an efficiency of the stage. The main attention is paid to turbulence modelling. Besides a linear model based on Boussinesq approximation a nonlinear model based on higher order closure formula is used. Also an effect of a curvature correction is taken into account.


Introduction
The axial steam turbine of low power (around hundreds kilowatts to units megawatts) are typically designed in configuration with the drum-type rotor.It means that the rotor blades are carried directly by the rotor of large diameter.In case of low power turbines, in order to reduce a production costs a blades with free ends that are not equipped with the shroud are still in use.For more reduction of the production cost the prismatic blades are sometimes used as it is in the present work.In this configuration it is necessary to maintain the radial gap under the hub-end of the stator blades and above the tip-end of the rotor blades.However, flow through the radial gap generates an intensive secondary vortices which interacts with following blade wheel.This interaction leads to strong swirl of flow field in the peripheral parts of the blade span.This work focuses on numerical modelling of intensively swirl flow in case of the axial turbine stage using the URANS (Unstaedy Reynolds Averaged Navier-Stokes equations) framework.The main attention is paid on the turbulence model closure formula.
Presented work was carried out under the project: TA02021336 -Research and Development of Small Sized Turbine Stage (see [1]).Project deals with research of an aerodynamics and energy characteristics of the axial stage of small sized and low power steam turbine.
Also some previous works carried out under above mentioned project deal with numerical modelling of flow in the axial stage with blades that are not equipped with the shroud.In [2] there is discussed an effect of the radial gap dimension on steady flow in stator wheel without interaction with the rotor blades.Results of numerical and experimental investigation are compared in [3].There is also discussed an influence of the interface between stator and rotor part of the computational domain.Two different kinds of the interface are compared in [3].Paper [4] deals a e-mail: straka@vzlu.czwith clarification the mechanism of the rotor blade interaction with secondary vortices generated behind the radial gap under the hub-end of the stator blade.Also in that paper computational results are compared with experimental data.There is also discussed the possibility of reduction the computational domain by the method of scaling, which technique is applied also in present work.
Note that in [2][3][4] there the two-equation k − ω turbulence model with first order closure formula based on common Boussinesq approximation was used.However in general the turbulence models based on first order closure formula are not able to sufficiently well predict the secondary flow structures.Therefore present work is focussed on turbulence model based on nonlinear second order closure formula.

Geometric setup
Figure 1 shows the computational domain with marked inlet and outlet boundaries and sliding mesh interface between stator and rotor parts of the computational domain.It is seen that the computational domain contains only one stator and one rotor blade.Such a reduction of the actual number of blades is achieved by method of scaling.The scale factor of the stator blade is given by 2N s /(N s + N r ), where N s is real number of the stator blades and N r is real number of the rotor blades.The scale factor of the rotor blade is defined analogously.In our case the scale factor is 0.909 for the stator blade and 1.09 for the rotor blade.Note that scaling has no effect on static aerodynamics loading of the blades, because the pitch to chord ratio of blades is preserved.Scaling of the blade profile causes little change of the chord based Reynolds number, but it is minor problem in case of values of scale factors near 1.The method of scaling has, however, major impact on prediction of dynamic aerodynamic loading of the blades.Therefore this method is unusable in cases where such that prediction is a main scope.In our case we want to compare prediction of the rotor blade interaction with secondary vortices generated behind the radial gap under the hub-end of the stator blade using turbulence models with different closure formulas.After scaling, the chord of the stator blade is 29.0 mm and the chord of the rotor blade is 28.1 mm.The hub diameter is 330 mm and diameter of the outer casing is 430 mm.Between the hub-end of the stator blade and the hub-wal there is the radial gap of 1 mm.Also between the tip-end of the rotor blade and the outer casing there is the radial gap of 1.4 mm.The axial distance between the stator trailing edge and the rotor leading edge is 15 mm.The inlet boundary is placed 24 mm before the stator leading edge and the outlet boundary is placed 15 mm behind the rotor trailing edge.The sliding mesh interface between stator and rotor parts of the computational domain is placed in the middle of the axial gap between the stator and the rotor blades.

Flow conditions
Presented results were calculated for the total/static expansion ratio 1.4, the velocity ratio u/c is = 0.6 (u is the circumferential velocity at middle diameter and c is is the isentropic velocity at the hub diameter behind the rotor trailing edge).It corresponds with the inlet total pressure 120 kPa, the outlet pressure 85.7 kPa, speed of revolutions 7430 RPM, the total inlet temperature of 303 K.The isentropic outlet Mach number is M is = 0.7, while the outlet isentropic Reynolds number based on blade chord is Re is ≈ 3 × 10 5 .

Numerical aparatus
Flow through the axial turbine stage is modeled as unsteady, 3D, compressible, viscous, fully turbulent flow of the perfect gas.Simulation is performed using the inhouse numerical code based on solution of system of the RANS equation and two-equation k − ω turbulence model.Used numerical code [5,6] is focused on turbomachinery application.The interface between stator and rotor part of the computational domain is performed via the "sliding mesh" method, where communication between stationary stator domain and rotating rotor domain is based on interpolation.System of the governing equations reads: where The turbulent viscosity μ t and the Reynolds stress term τ t i j will be described later in Sect. 5. System of equations (1) to ( 5) we can formally rewrite in vector form as follows: where W is vector of the conservative variables, F is vector of the inviscid and viscous fluxes and Q is vector of the source terms.System (23) is discretized using the cellcentered finite-volume method on multi-block structured 02115-p.2mesh of hexahedral cells.Temporal discretization is performed using the second-orded backward Euler formula in implicit form, which is realized through a dual iterative process: where W n, 0 := W n , ΔW n, ν+1/2 = ΔW n, ν+1 − ΔW n, ν and W n+1 := W n, ν * , where W n, ν * is a steady solution in the dual iterative process.In scheme (24) there n means the time level, Δt is the time step, ν is index of the dual iteration, i is index of current cell in the computational mesh, j is index of the neighbouring cells and R ≈ ∂F j /∂x j is the numerical approximation of the inviscid and the viscous fluxes.The inviscid numerical fluxes are calculated using the exact solution of the 1D Rieman problem in normal direction to the cell edges.The viscous numerical fluxes are calculated using the central scheme using the Green-Gauss theorem on a dual cells.Higher order of accuracy in space is obtained using linear reconstruction with the Van Leer's slope limiter.The computational domain was discretized by structured multiblock mesh of 1.5 millions hexahedral cells which is displayed in figure 2.

Closure formula
As it was mentioned above, this work is focused on using of the turbulence model based on nonlinear second order closure formula.Among the many turbulence models suitable for prediction of secondary flows a quadratic explicit algebraic Reynolds stress model proposed by Rumsey and Gatski [7] was chosen.This model takes into account the production-to-dissipation rate ratio as well as the stress tensor anisotropyo.On the other hand, because the stress tensor is modeled via algebraic formula, this model has comparable computational cost as the linear eddy-viscosity turbulence models based on the Boussinesq approximation.Because the algebraic stress models are in general derived based on weak-equilibrium assumption, it partially lose sensitivity to streamline curvature.This may be problem in case of strong swirl flow field.This deficiency can be removed by modelling the convection of the anisotropy by some means.The curvature correction method based on the strain rate tensor, that is described in [8], is used in this work.This turbulence model will be referred here as EASM (Explicit Algebraic Stress Model).
Results calculated using models EASM are in this work compared with those from [4].Because the results published in [4] were computed using the turbulence model based on the first order closure formula, this model is also described here.Turbulence model based on the first order closure formula will be here referred as EVM (Eddy Viscosity Model).
Note that used turbulence models are based on two transport equations for turbulent kinetic energy k and specific dissipation rate of turbulent energy ω given by eq. ( 4) and ( 5).

EVM turbulence model
This model is based on the Boussinesq approximation.The Reynolds stress tensor τ t i j is modeled as a linear function of the strain rate tensor S i j : where the turbulent viscosity μ t = μ evm t is given by eq. ( 18).The model constants: β * = 0.09, β = 0.075, γ = 0.5532, σ k = 0.5, σ ω = 0.5 and σ d = 0.5 corresponds with the TNT k − ω turbulence model proposed by Kok [9].

EASM turbulence model
In this model the Reynolds stress tensor τ t i j is modeled as a quadratic function of the strain rate tensor S i j : It is clear that eq. ( 26) extends relation (25) about anisotropy part (term in square brackets).The turbulent viscosity μ t is now given as where the turbulent time scale τ is defined as where C τ = 6.0 is a model constant.The nominal level of the variable coefficient C μ in a zero-pressure-gradient loglayer is approximately 0.09.In general case value of C μ is obtained by solving following cubic equation: where The correct root to choose from this equation is the root with the lowest real part.See [7] for more details.Other parameters are gived by: where

Curvature correction
The rotation rate tensor Ω i j is defined now as: where and A 0 = −0.72,III S = S kl S lm S mk and ε i jk is the Levi-Civita symbol.
In eq. ( 38) is the material derivative of the strain rate tensor.The temporal derivative of the strain rate tensor is approximated using first order backward differencing scheme, while the convective part is discretized by standard manner of finite volume method.Following [10] a spatial filter for computed tensor Ω (r)  i j is applied to prevent some problems with numerical oscillations of DS /D t.

Computational results
The effect of flow through the radial gap under the stator blade is such that stream flowing from the radial gap generates (in interaction with main stream in blade channel) massive vortex above the stream and beside the vertical surface of the suction side as illustrates a scheme in Fig. 3.This upper vortex is in given regimes and geometrical configurations intensive enough to deform the stream flowing  from the radial gap which leads to separation of the stream from the wall and to generation of the other vortex under the stream.Both vortices then interact with the rotor blade so that they deforms in upward direction and the are carried on the suction side almost up to half of the span.This causes increase of the dissipation rate on the lower half of span on the suction side, as illustrates distribution of the dimensionless entropy index s = p/ρ κ (computed from the normalized pressure and density) on the rotor blade surface that is shown in figure 4.
The difference between the linear EVM model and the nonlinear EASM model, which takes the anisotropy into account, lies in prediction of forming of the upper and lower vortices.The EVM model predicts this phenomena much simpler than the EASM model.The EVM model predicts the lower vortex in parallel along the upper vortex up to the point of interaction with the leading edge of the rotor blade, while the EASM model predicts distortion of the lower vortex that is being wound around the upper vortex as is depicted schematically in figure 5 and in form of isosurface of (|Ω| 2 − |S | 2 ) in figure 6.This difference between EVM and EASM models has its origin in anisotropic part of the Reynolds stress tensor, which is ignored in the EVM model.
Another difference between both turbulent models lies in production of the turbulent kinetic energy.Figure 7 shows field of the turbulent kinetic energy in circumferential cross section at 10, 50 and 90 percent of span.The

Figure 3 .
Figure 3. Scheme of flow behind the radial gap under the stator blade.

Figure 4 .
Figure 4. Distribution of the dimensionless entropy index s = p/ρ κ on the rotor blade surface.

02115-p. 4 EPJFigure 5 .Figure 6 .
Figure 5. Scheme of formation of the secondary vortices behind the radial gap under the stator blade.

Figure 9 .
Figure 9. Angle of the velocity vector from the axial direction in tangential plane.Left column -EVM model, right column -EASM model, top -90% of span, middle -50_ of span, bottom -10% of span.