- Email: [email protected]
Numerical simulation of circumferentially averaged flow in a turbine
Applied Mathematics and Computation xxx (2015) xxx–xxx
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical simulation of circumferentially averaged flow in a turbine Jirˇí Fürst, Jaroslav Forˇt, Jan Halama ⇑, Jirˇí Holman, Jan Karel, Vladimír Prokop, David Trdlicˇka FME CTU Prague, Karlovo nám. 13, 121-35 Prague, Czech Republic
a r t i c l e
i n f o
Keywords: Turbine Finite volume method Loss model
a b s t r a c t The paper refers about the development of a fast computational code, which should be able to provide an approximate information about the three-dimensional flow field in a multistage turbine. The code is based upon the solution of circumferentially averaged Euler equations coupled with the thermodynamic, geometry and loss prediction models. The computational domain is the meridional cut of a turbine. The Euler equations are solved by a finite volume solver with the AUSM type flux. Initial tests showed, that developed solver is able to predict well radial distributions of flow parameters upstream and downstream considered blade cascades at a fraction of CPU time compared to fully three-dimensional simulations. Ó 2015 Elsevier Inc. All rights reserved.
1. Introduction The paper refers about the development of a fast computational code intended for early stages of turbine design. There exist different methods ranging from quasi 1D solvers [1] to solvers based on the circumferentially averaged Euler [2] or Navier–Stokes [3] equations. The aim of this work is to develop a code, which is able to provide an approximate information about the three-dimensional flow field in a multistage turbine. It is based upon the solution of steady circumferentially averaged flow. The computational domain is two-dimensional. It is defined by the meridional cut of a turbine, i.e. the real shapes of hub and tip casing are parts of the domain boundary. The shapes of the pressure and the suction surfaces of the blade, which cannot be described by the definition of the computational domain, are included in the form of additional parameters. The variable thickness of blades is taken into account using the channel width parameter. The shape of the blades is approximated by the shape of the midplane, which is located in the middle between the suction and the pressure surfaces of the blade. This plane is implemented in the form of normal vectors. Similarly to [2] the flow is modeled by the circumferentially averaged Euler equations, which are coupled with geometry and loss models. The solution domain is discretized by a structured grid. The flow solver uses a finite volume method with AUSM (Advection Upstream Splitting Method) type flux [7]. Initial tests show, that developed solver is able to predict well the radial distributions of flow parameters upstream and downstream the considered blade cascades at short CPU time compared to fully three-dimensional simulations. The used finite volume method has advantage over previous methods based on a streamlines and stream functions, which did not work well for transonic speeds and they did not generally guarantee the conservation of mass, momentum and energy.
⇑ Corresponding author. E-mail address: [email protected] (J. Halama). http://dx.doi.org/10.1016/j.amc.2015.01.086 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
2
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
2. Flow model Consider the system of cylindrical coordinates r; u and z, which is fixed to a single blade row of axial turbine. The z axis coincides with the axis of the cascade. The cascade can rotate around its axis with angular velocity x. The Euler equations in the considered system of coordinates read
e @H f 1 @ðr e e FÞ 1 @G @W e þ þ þ ¼ B; @t r @r r @ u @z 3
2
q 7 6 6 qv r 7 7 6 f 7 W ¼6 6 qv u 7; 7 6 4 qv z 5 e 2
ð1Þ
3
2
7 6 6 qv 2r þ p 7 7 6 e 7 F ¼6 6 qv r v u 7; 7 6 4 qv r v z 5 v r ðe þ pÞ
6 6 6 e G¼6 6 6 4
2
qv r
qv u qv u v r qv 2u þ p qv u v z
v u ðe þ pÞ
3
2
7 7 7 7; 7 7 5
6 6 6 e H¼6 6 6 4
qv z qv z v r qv z v u qv 2z þ p
v z ðe þ pÞ
3 7 7 7 7; 7 7 5
3
0
7 6 vu p 6 q r þ r þ qr x2 þ 2qxv u 7 7 6 v v r u e 7; B¼6 q r 2qxv r 7 6 7 6 5 4 0 2
0 where the symbol t denotes time, q density, v r ; v u and v z velocity components, e the total energy per unit volume and p the pressure. We consider the steady solution, therefore we might take a single blade passage (i.e. blade-to-blade channel together with sufficiently long upstream and downstream parts) as the solution domain D. The ranges for the coordinates z and r are defined by the domain Dzr (projection of D onto zr plane). The range for the coordinate u is limited by u1 ðr; zÞ and u2 ðr; zÞ. The idea of circumferential averaging is equivalent to the finite volume discretization of D with a single layer of cells in the tangential direction. Consider an arbitrary cell of this discretization defined by V Dzr and u1 ðr; zÞ < u < u2 ðr; zÞ. The integral form of the Eq. (1) then reads
ZZ
Z V
u2
u1
! ZZ e @H f 1 @ðr e e FÞ 1 @G @W rdudrdz ¼ þ þ þ @t r @r r @ u @z
Z
u2
e udrdz: Brd
ð2Þ
u1
V
Swapping the derivative and the integration with respect to u and using the mean value theorem for the left hand side of the Eq. (2) yields
ZZ V
@ðbrWÞ @ðbrFÞ @ðbrHÞ þ þ þ ðF 1 ; G1 ; H1 Þ~ n1 ðF 2 ; G2 ; H2 Þ~ n2 drdz; @t @r @z
where bðr; zÞ ¼ u2 ðr; zÞ u1 ðr; zÞ is the width of blade channel and V is the finite volume of the discretization of Dzr , see e and H f; e e with respect to u. The vector F i F; G examples in the Figs. 2 or 10(b). The vectors W; F; G and H are averages of W e for u ¼ u and is the vector F i
@ ui 1 @ ui ~ ; ; ; ni ¼ r @z @r
ð3Þ
is the normal vector to the surface u ¼ ui ðr; zÞ. The surfaces u ¼ ui ðr; zÞ coincide partly with the blade surface. Therefore we T need to apply the non-permeability condition in the form ðF i ; Gi ; Hi Þ~ ni ¼ ½0; pi~ ni ; 0 on the part of u ¼ ui ðr; zÞ surfaces and for n1 ¼ ~ n2 and p1 ¼ p2 . The non-permeabilthe remaining part to apply the periodicity conditions, i.e. ðF 1 ; G1 ; H1 Þ ¼ ðF 2 ; G2 ; H2 Þ, ~ ity condition yields
3 3 2 0 0 7 7 6 6 6 r@ u2 =@r 7 6 r@ u1 =@r 7 7 7 6 6 7 7 6 ðF 2 ; G2 ; H2 Þ~ n2 ðF 1 ; G1 ; H1 Þ~ n1 ¼ p2 6 6 1 7 p1 6 1 7: 7 7 6 6 4 r@ u2 =@z 5 4 r@ u1 =@z 5 2
0
ð4Þ
0
The above substitution (4) can be used also in the case of periodicity conditions, where both sides turn to zero. The circumferentially averaged integral form of governing equations is therefore
ZZ
@ðbrWÞ dzdr þ @t V
ZZ
div ðbrH; brFÞdzdr ¼ V
ZZ
brQdzdr;
ð5Þ
V
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
3
with
3 0 7 6 6 @b=@r 7 7 p6 0 7 Q ¼Bþ 6 7; b6 7 6 4 @b=@z 5 0 2
f; e e The set of Eqs. (5) have to be closed by the therF and H. and the u-averages W; F and H have the same components as W modynamic model, i.e. pressure as a function of W. The equations are further coupled with a geometry model described in the next section as well as with a loss prediction model in order to approximate the three-dimensionality of flow and the dissipation of kinetic energy respectively. 3. Computational algorithm Computational algorithm is based on explicit finite volume method for the system (5) with additional terms due to geometry and loss prediction models. The finite volume method uses AUSM type flux. The example of algorithm with the first order discretization in time reads nþ1
ðbrW ÞV
n
¼ ðbrW ÞV
Dt X
lðVÞ
n n n ðbrHÞl Dr l ðbrFÞl Dzl þ ðbrQ ÞV þ SV ;
ð6Þ
l
where subscript l is used for parts of boundary of V and superscript n for the time level. Both geometry and loss prediction T
f tot ; 0 with external force ~ f tot ¼ ~ f b þ~ f d , which is non-zero models are included in the form of additional term SV ¼ ½0; q~ b ~ between the blades. The force f models the change of flow direction due to blades. It is perpendicular to the midplane uðr; zÞ ¼ ðu1 ðr; zÞ þ u2 ðr; zÞÞ=2 þ ucor , where ucor includes incidence and deviation corrections from the loss prediction f d models the deceleration of flow due to dissipation of kinetic energy. The computational domain, i.e. a model. The force ~ domain in the zr plane, is discretized by structured quadrilateral grid. The coupling between finite volume method and loss model has following steps: Grid lines j ¼ const are considered as streamlines. The solution W ni;j for each j ¼ const line gives an input data for the loss prediction model, which returns the total pressure loss recomputed further as the entropy rise Dsloss j . The constant loss model without ucor correction and AMDC-KO loss model according to the book [4] are implemented. Computation of the entropy rise Dsnj related to the numerical solution W nj , i.e. the amount of numerical dissipation. f dj related to Dsloss Update of force D~ Dsnj . When the entropy rise of the numerical solution and the loss model equilij f d will become constant. brates, the force ~ f dj along j ¼ const line between leading and trailing edge. Redistribution of ~ ~ Update of force Df bi;j proportional to the component of the velocity perpendicular to the midplane. When the velocity perf b will stay constant. The force ~ f b is applied between leading and pendicular to the midplane goes to the zero, the force D~ i;j
trailing edges to mimic the guidance of the flow by the blades. Evaluation of unknowns in ghost cells – implementation of common boundary conditions for inlet, outlet, hub wall and tip wall boundaries. f bi;j and ~ f di;j are relaxed to ensure the robustness of the by the finite volume method. External forces D~ Computation of W nþ1 i;j whole iterative algorithm. 4. Case 1: flow in a high pressure core stage of a gas turbine The geometry of the first test case is rather simple with short blades and moderate 3D shaping, see the Fig. 1. The computational grid with highlighted positions of leading and trailing edges is plotted in the Fig. 2. The positions of the leading and the trailing edges are important since they define the region where the geometry and the loss prediction models are switched on. The Fig. 3 shows the contours of the relative channel width b. The values of b are ranging from zero to unity, where zero would theoretically mean fully closed channel and unity the part without blades. The geometry is available in [5]. The discussed results correspond to the cold air test [5] with the stage inlet total pressure 101:3 kPa, the stage inlet total temperature 288:2 K and the stage outlet pressure 22:8 kPa. The flow at the stage inlet has axial direction and the rotational speed of the rotor cascade is 8081 RPM. The dissipation of kinetic energy of flow has been modeled by the constant loss model with 5:5% loss for each cascade (design loss [5]). The models of flow incidence and deviations have not been used. The contours of the relative Mach number (note the discontinuity between the absolute and the relative frame) are plotted
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
4
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
Fig. 1. Example of the shape of the midplane for the high pressure core stage of a gas turbine.
Fig. 2. Computational domain and grid for the case of high pressure core stage of a gas turbine.
Fig. 3. Relative width of blade channel for the case 1.
in the Fig. 4. There is only a weak gradient of solution in the radial direction due to relatively short blades and moderate 3D shaping. The streamlines in the Fig. 5 show non-zero radial velocity in the rotor. The qualitative comparison of flow angles at the mid-height is available in the Table 1. The absolute flow angle a1 downstream the stator and the relative flow angle b2 downstream the rotor obtained by the current method are equal to the design geometry angles. It means, that the geometry model works correctly. One could expect, that a deviation model will move the obtained values towards the values from experiment. The relative flow angle b1 upstream the rotor is influenced by several factors, therefore the difference with respect to the design and experiment might be bigger. Moreover the original paper [5] notes the slight difference between
Fig. 4. Relative Mach number contours in the meridional plane for the case 1.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
5
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
Fig. 5. Streamlines in the meridional plane for the case 1.
Table 1 The absolute and the relative flow angles downstream the stator (a1 and b1 ) and the relative flow angle downstream the rotor (b2 ) at the mid-height and the total mass flow.
Design [5] Experiment [5] 3D Euler [6] Current method
a1
b1
b2
_ m
73:0 72:2 72:7 73:0
50:8 48:3 48:0 52:2
59:6 56:4 55:0 59:7
3:71 3:86 3:95 3:77
the design geometry and the real geometry of the stator blade. The total mass flow estimated by the current method fits within the range given by the values from design and the experiment. The results of the three-dimensional simulation [6] are also included in the Table 1. The effect of the geometry model, i.e. the shape of midplane between the blades and the width of flow channel is demonstrated by the pressure distributions between the leading and the trailing edges at three different radial positions. These pressure distributions, which are denoted by symbols in the Figs. 6–8, should approximate the
Fig. 6. Distribution of normalized pressure along stator (left) and rotor (right) blade profiles at the hub section. Full line denotes the solution from 6, symbols denote the results of current method.
Fig. 7. Distribution of normalized pressure along stator (left) and rotor (right) blade profiles at the mid-height. Full line denotes the solution from [6], symbols denote the results of current method.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
6
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
Fig. 8. Distribution of normalized pressure along stator (left) and rotor (right) blade profiles at the tip section. Full line denotes the solution from [6], symbols denote the results of current method.
Fig. 9. Radial distributions of relative flow angle and relative Mach number upstream the stator cascade, in the gap between stator and rotor and downstream the rotor.
Fig. 10. Geometry and grid for the case 2.
pressure distributions [6] along the blade profiles, which are plotted by the full lines. Since there is no pressure neither suction side of the blade in the current simulation, the distribution of pressure achieved by the current method should be understood as an average. Pressure distributions in the Figs. 6–8 show, that current method is able to approximate well the
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
7
Fig. 11. Mach number contours for the case 2.
expansion through the blade channel. The differences in the vicinity of the rotor leading edge are most probably caused by differences in the relative flow angle b1 , discussed already for the Table 1. The examples of radial distributions of flow angles and relative Mach number are plotted in the Fig. 9. 5. Case 2: flow in the last stator cascade of a steam turbine The second test case is the flow in the last stator cascade from of a low pressure part of a steam turbine. The geometry of the cascade has significant 3D shaping, see the Fig. 10. The grid of the middle section from 3D Euler simulation in the Fig. 10(a) corresponds to the shape of midplane for the current method. The comparison of Mach number contours in the middle section of 3D Euler simulation [8], see the Fig. 11(a), with the results of current method used without loss prediction model, see the Fig. 11(b) shows how the current method approximates transonic flow field in more complex geometry. Stronger gradients in the case of 3D Euler simulation in the outlet part close to the hub (the right bottom part of the domain in Fig. 11(a)) are due to the presence of trailing edge shock waves. These shock waves cannot be modeled by the current method (the geometry model takes into account only direction and contraction of flow channel), therefore these gradients are missing in the Fig. 11(b). 6. Conclusions Current method based on the iterative coupling of circumferentially averaged Euler equations with the geometry and the loss prediction models can approximate the three-dimensional flow field in a turbine. The dissipative effects are included strictly by the loss prediction model, therefore the computational grid does not need to be refined along the hub and tip casings. The implementation of loss prediction model compensates the numerical dissipation of a used finite volume method. The complete simulation for a turbine stage is completed in a relatively short time (several minutes using common computer). The method is suitable for the simulation of flow in a multistage configuration during the early steps of turbine design, where the information about mass flow and radial distributions of flow angles, pressure and velocity are essential. Acknowledgments This work has been supported by the Grant TACR TE01020020.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
8
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
References [1] [2] [3] [4] [5]
O. Leonard, O. Adam, A quasi-one-dimensional cfd model for multistage turbomachines, J. Therm. Sci. 17 (1) (2008) 7–20. G. Fahua, M.R. Anderson, Cfd-based throughflow solver in turbomachinery design systems, in: ASME Turbo Expo 2007, ASME, Montreal, 2007. J.F. Simon, O. Leonard, Modeling of 3-d losses and deviations in a throughflow analysis tool, J.Therm. Sci. 16 (3) (2007) 208–214. R.H. Aungier, Turbine Aerodynamics: Axial-Flow and Radial-Inflow Turbine Design and Analysis, ASME Press, 2006. T.P. Moffit, E.M. Szanca, W.J. Whitney, F.P. Behning, Design and cold-air test of single uncooled core turbine with high work output, Tech. Rep. NASA Technical paper 1680, Lewis Research Center, Cleveland, Ohio (May 1980). [6] T. Arts, Calculation of the three-dimensional, steady, inviscid flow in a transonic axial turbine stage, in: 29th International Gas Turbine Conference and Exhibit Amsterdam, ASME, 1984. paper 84-GT-76. [7] M.S. Liou, C.J. Steffen, A new flux splitting scheme, J. Comput. Phys. 107 (1) (1993) 23–39. [8] J. Forˇt, J. Fürst, J. Halama, K. Kozel, Numerical simulation of 3D transonic flow through cascades, Math. Bohem. 126 (2) (2001) 353–361.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical simulation of circumferentially averaged flow in a turbine Jirˇí Fürst, Jaroslav Forˇt, Jan Halama ⇑, Jirˇí Holman, Jan Karel, Vladimír Prokop, David Trdlicˇka FME CTU Prague, Karlovo nám. 13, 121-35 Prague, Czech Republic
a r t i c l e
i n f o
Keywords: Turbine Finite volume method Loss model
a b s t r a c t The paper refers about the development of a fast computational code, which should be able to provide an approximate information about the three-dimensional flow field in a multistage turbine. The code is based upon the solution of circumferentially averaged Euler equations coupled with the thermodynamic, geometry and loss prediction models. The computational domain is the meridional cut of a turbine. The Euler equations are solved by a finite volume solver with the AUSM type flux. Initial tests showed, that developed solver is able to predict well radial distributions of flow parameters upstream and downstream considered blade cascades at a fraction of CPU time compared to fully three-dimensional simulations. Ó 2015 Elsevier Inc. All rights reserved.
1. Introduction The paper refers about the development of a fast computational code intended for early stages of turbine design. There exist different methods ranging from quasi 1D solvers [1] to solvers based on the circumferentially averaged Euler [2] or Navier–Stokes [3] equations. The aim of this work is to develop a code, which is able to provide an approximate information about the three-dimensional flow field in a multistage turbine. It is based upon the solution of steady circumferentially averaged flow. The computational domain is two-dimensional. It is defined by the meridional cut of a turbine, i.e. the real shapes of hub and tip casing are parts of the domain boundary. The shapes of the pressure and the suction surfaces of the blade, which cannot be described by the definition of the computational domain, are included in the form of additional parameters. The variable thickness of blades is taken into account using the channel width parameter. The shape of the blades is approximated by the shape of the midplane, which is located in the middle between the suction and the pressure surfaces of the blade. This plane is implemented in the form of normal vectors. Similarly to [2] the flow is modeled by the circumferentially averaged Euler equations, which are coupled with geometry and loss models. The solution domain is discretized by a structured grid. The flow solver uses a finite volume method with AUSM (Advection Upstream Splitting Method) type flux [7]. Initial tests show, that developed solver is able to predict well the radial distributions of flow parameters upstream and downstream the considered blade cascades at short CPU time compared to fully three-dimensional simulations. The used finite volume method has advantage over previous methods based on a streamlines and stream functions, which did not work well for transonic speeds and they did not generally guarantee the conservation of mass, momentum and energy.
⇑ Corresponding author. E-mail address: [email protected] (J. Halama). http://dx.doi.org/10.1016/j.amc.2015.01.086 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
2
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
2. Flow model Consider the system of cylindrical coordinates r; u and z, which is fixed to a single blade row of axial turbine. The z axis coincides with the axis of the cascade. The cascade can rotate around its axis with angular velocity x. The Euler equations in the considered system of coordinates read
e @H f 1 @ðr e e FÞ 1 @G @W e þ þ þ ¼ B; @t r @r r @ u @z 3
2
q 7 6 6 qv r 7 7 6 f 7 W ¼6 6 qv u 7; 7 6 4 qv z 5 e 2
ð1Þ
3
2
7 6 6 qv 2r þ p 7 7 6 e 7 F ¼6 6 qv r v u 7; 7 6 4 qv r v z 5 v r ðe þ pÞ
6 6 6 e G¼6 6 6 4
2
qv r
qv u qv u v r qv 2u þ p qv u v z
v u ðe þ pÞ
3
2
7 7 7 7; 7 7 5
6 6 6 e H¼6 6 6 4
qv z qv z v r qv z v u qv 2z þ p
v z ðe þ pÞ
3 7 7 7 7; 7 7 5
3
0
7 6 vu p 6 q r þ r þ qr x2 þ 2qxv u 7 7 6 v v r u e 7; B¼6 q r 2qxv r 7 6 7 6 5 4 0 2
0 where the symbol t denotes time, q density, v r ; v u and v z velocity components, e the total energy per unit volume and p the pressure. We consider the steady solution, therefore we might take a single blade passage (i.e. blade-to-blade channel together with sufficiently long upstream and downstream parts) as the solution domain D. The ranges for the coordinates z and r are defined by the domain Dzr (projection of D onto zr plane). The range for the coordinate u is limited by u1 ðr; zÞ and u2 ðr; zÞ. The idea of circumferential averaging is equivalent to the finite volume discretization of D with a single layer of cells in the tangential direction. Consider an arbitrary cell of this discretization defined by V Dzr and u1 ðr; zÞ < u < u2 ðr; zÞ. The integral form of the Eq. (1) then reads
ZZ
Z V
u2
u1
! ZZ e @H f 1 @ðr e e FÞ 1 @G @W rdudrdz ¼ þ þ þ @t r @r r @ u @z
Z
u2
e udrdz: Brd
ð2Þ
u1
V
Swapping the derivative and the integration with respect to u and using the mean value theorem for the left hand side of the Eq. (2) yields
ZZ V
@ðbrWÞ @ðbrFÞ @ðbrHÞ þ þ þ ðF 1 ; G1 ; H1 Þ~ n1 ðF 2 ; G2 ; H2 Þ~ n2 drdz; @t @r @z
where bðr; zÞ ¼ u2 ðr; zÞ u1 ðr; zÞ is the width of blade channel and V is the finite volume of the discretization of Dzr , see e and H f; e e with respect to u. The vector F i F; G examples in the Figs. 2 or 10(b). The vectors W; F; G and H are averages of W e for u ¼ u and is the vector F i
@ ui 1 @ ui ~ ; ; ; ni ¼ r @z @r
ð3Þ
is the normal vector to the surface u ¼ ui ðr; zÞ. The surfaces u ¼ ui ðr; zÞ coincide partly with the blade surface. Therefore we T need to apply the non-permeability condition in the form ðF i ; Gi ; Hi Þ~ ni ¼ ½0; pi~ ni ; 0 on the part of u ¼ ui ðr; zÞ surfaces and for n1 ¼ ~ n2 and p1 ¼ p2 . The non-permeabilthe remaining part to apply the periodicity conditions, i.e. ðF 1 ; G1 ; H1 Þ ¼ ðF 2 ; G2 ; H2 Þ, ~ ity condition yields
3 3 2 0 0 7 7 6 6 6 r@ u2 =@r 7 6 r@ u1 =@r 7 7 7 6 6 7 7 6 ðF 2 ; G2 ; H2 Þ~ n2 ðF 1 ; G1 ; H1 Þ~ n1 ¼ p2 6 6 1 7 p1 6 1 7: 7 7 6 6 4 r@ u2 =@z 5 4 r@ u1 =@z 5 2
0
ð4Þ
0
The above substitution (4) can be used also in the case of periodicity conditions, where both sides turn to zero. The circumferentially averaged integral form of governing equations is therefore
ZZ
@ðbrWÞ dzdr þ @t V
ZZ
div ðbrH; brFÞdzdr ¼ V
ZZ
brQdzdr;
ð5Þ
V
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
3
with
3 0 7 6 6 @b=@r 7 7 p6 0 7 Q ¼Bþ 6 7; b6 7 6 4 @b=@z 5 0 2
f; e e The set of Eqs. (5) have to be closed by the therF and H. and the u-averages W; F and H have the same components as W modynamic model, i.e. pressure as a function of W. The equations are further coupled with a geometry model described in the next section as well as with a loss prediction model in order to approximate the three-dimensionality of flow and the dissipation of kinetic energy respectively. 3. Computational algorithm Computational algorithm is based on explicit finite volume method for the system (5) with additional terms due to geometry and loss prediction models. The finite volume method uses AUSM type flux. The example of algorithm with the first order discretization in time reads nþ1
ðbrW ÞV
n
¼ ðbrW ÞV
Dt X
lðVÞ
n n n ðbrHÞl Dr l ðbrFÞl Dzl þ ðbrQ ÞV þ SV ;
ð6Þ
l
where subscript l is used for parts of boundary of V and superscript n for the time level. Both geometry and loss prediction T
f tot ; 0 with external force ~ f tot ¼ ~ f b þ~ f d , which is non-zero models are included in the form of additional term SV ¼ ½0; q~ b ~ between the blades. The force f models the change of flow direction due to blades. It is perpendicular to the midplane uðr; zÞ ¼ ðu1 ðr; zÞ þ u2 ðr; zÞÞ=2 þ ucor , where ucor includes incidence and deviation corrections from the loss prediction f d models the deceleration of flow due to dissipation of kinetic energy. The computational domain, i.e. a model. The force ~ domain in the zr plane, is discretized by structured quadrilateral grid. The coupling between finite volume method and loss model has following steps: Grid lines j ¼ const are considered as streamlines. The solution W ni;j for each j ¼ const line gives an input data for the loss prediction model, which returns the total pressure loss recomputed further as the entropy rise Dsloss j . The constant loss model without ucor correction and AMDC-KO loss model according to the book [4] are implemented. Computation of the entropy rise Dsnj related to the numerical solution W nj , i.e. the amount of numerical dissipation. f dj related to Dsloss Update of force D~ Dsnj . When the entropy rise of the numerical solution and the loss model equilij f d will become constant. brates, the force ~ f dj along j ¼ const line between leading and trailing edge. Redistribution of ~ ~ Update of force Df bi;j proportional to the component of the velocity perpendicular to the midplane. When the velocity perf b will stay constant. The force ~ f b is applied between leading and pendicular to the midplane goes to the zero, the force D~ i;j
trailing edges to mimic the guidance of the flow by the blades. Evaluation of unknowns in ghost cells – implementation of common boundary conditions for inlet, outlet, hub wall and tip wall boundaries. f bi;j and ~ f di;j are relaxed to ensure the robustness of the by the finite volume method. External forces D~ Computation of W nþ1 i;j whole iterative algorithm. 4. Case 1: flow in a high pressure core stage of a gas turbine The geometry of the first test case is rather simple with short blades and moderate 3D shaping, see the Fig. 1. The computational grid with highlighted positions of leading and trailing edges is plotted in the Fig. 2. The positions of the leading and the trailing edges are important since they define the region where the geometry and the loss prediction models are switched on. The Fig. 3 shows the contours of the relative channel width b. The values of b are ranging from zero to unity, where zero would theoretically mean fully closed channel and unity the part without blades. The geometry is available in [5]. The discussed results correspond to the cold air test [5] with the stage inlet total pressure 101:3 kPa, the stage inlet total temperature 288:2 K and the stage outlet pressure 22:8 kPa. The flow at the stage inlet has axial direction and the rotational speed of the rotor cascade is 8081 RPM. The dissipation of kinetic energy of flow has been modeled by the constant loss model with 5:5% loss for each cascade (design loss [5]). The models of flow incidence and deviations have not been used. The contours of the relative Mach number (note the discontinuity between the absolute and the relative frame) are plotted
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
4
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
Fig. 1. Example of the shape of the midplane for the high pressure core stage of a gas turbine.
Fig. 2. Computational domain and grid for the case of high pressure core stage of a gas turbine.
Fig. 3. Relative width of blade channel for the case 1.
in the Fig. 4. There is only a weak gradient of solution in the radial direction due to relatively short blades and moderate 3D shaping. The streamlines in the Fig. 5 show non-zero radial velocity in the rotor. The qualitative comparison of flow angles at the mid-height is available in the Table 1. The absolute flow angle a1 downstream the stator and the relative flow angle b2 downstream the rotor obtained by the current method are equal to the design geometry angles. It means, that the geometry model works correctly. One could expect, that a deviation model will move the obtained values towards the values from experiment. The relative flow angle b1 upstream the rotor is influenced by several factors, therefore the difference with respect to the design and experiment might be bigger. Moreover the original paper [5] notes the slight difference between
Fig. 4. Relative Mach number contours in the meridional plane for the case 1.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
5
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
Fig. 5. Streamlines in the meridional plane for the case 1.
Table 1 The absolute and the relative flow angles downstream the stator (a1 and b1 ) and the relative flow angle downstream the rotor (b2 ) at the mid-height and the total mass flow.
Design [5] Experiment [5] 3D Euler [6] Current method
a1
b1
b2
_ m
73:0 72:2 72:7 73:0
50:8 48:3 48:0 52:2
59:6 56:4 55:0 59:7
3:71 3:86 3:95 3:77
the design geometry and the real geometry of the stator blade. The total mass flow estimated by the current method fits within the range given by the values from design and the experiment. The results of the three-dimensional simulation [6] are also included in the Table 1. The effect of the geometry model, i.e. the shape of midplane between the blades and the width of flow channel is demonstrated by the pressure distributions between the leading and the trailing edges at three different radial positions. These pressure distributions, which are denoted by symbols in the Figs. 6–8, should approximate the
Fig. 6. Distribution of normalized pressure along stator (left) and rotor (right) blade profiles at the hub section. Full line denotes the solution from 6, symbols denote the results of current method.
Fig. 7. Distribution of normalized pressure along stator (left) and rotor (right) blade profiles at the mid-height. Full line denotes the solution from [6], symbols denote the results of current method.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
6
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
Fig. 8. Distribution of normalized pressure along stator (left) and rotor (right) blade profiles at the tip section. Full line denotes the solution from [6], symbols denote the results of current method.
Fig. 9. Radial distributions of relative flow angle and relative Mach number upstream the stator cascade, in the gap between stator and rotor and downstream the rotor.
Fig. 10. Geometry and grid for the case 2.
pressure distributions [6] along the blade profiles, which are plotted by the full lines. Since there is no pressure neither suction side of the blade in the current simulation, the distribution of pressure achieved by the current method should be understood as an average. Pressure distributions in the Figs. 6–8 show, that current method is able to approximate well the
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
7
Fig. 11. Mach number contours for the case 2.
expansion through the blade channel. The differences in the vicinity of the rotor leading edge are most probably caused by differences in the relative flow angle b1 , discussed already for the Table 1. The examples of radial distributions of flow angles and relative Mach number are plotted in the Fig. 9. 5. Case 2: flow in the last stator cascade of a steam turbine The second test case is the flow in the last stator cascade from of a low pressure part of a steam turbine. The geometry of the cascade has significant 3D shaping, see the Fig. 10. The grid of the middle section from 3D Euler simulation in the Fig. 10(a) corresponds to the shape of midplane for the current method. The comparison of Mach number contours in the middle section of 3D Euler simulation [8], see the Fig. 11(a), with the results of current method used without loss prediction model, see the Fig. 11(b) shows how the current method approximates transonic flow field in more complex geometry. Stronger gradients in the case of 3D Euler simulation in the outlet part close to the hub (the right bottom part of the domain in Fig. 11(a)) are due to the presence of trailing edge shock waves. These shock waves cannot be modeled by the current method (the geometry model takes into account only direction and contraction of flow channel), therefore these gradients are missing in the Fig. 11(b). 6. Conclusions Current method based on the iterative coupling of circumferentially averaged Euler equations with the geometry and the loss prediction models can approximate the three-dimensional flow field in a turbine. The dissipative effects are included strictly by the loss prediction model, therefore the computational grid does not need to be refined along the hub and tip casings. The implementation of loss prediction model compensates the numerical dissipation of a used finite volume method. The complete simulation for a turbine stage is completed in a relatively short time (several minutes using common computer). The method is suitable for the simulation of flow in a multistage configuration during the early steps of turbine design, where the information about mass flow and radial distributions of flow angles, pressure and velocity are essential. Acknowledgments This work has been supported by the Grant TACR TE01020020.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086
8
J. Fürst et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx
References [1] [2] [3] [4] [5]
O. Leonard, O. Adam, A quasi-one-dimensional cfd model for multistage turbomachines, J. Therm. Sci. 17 (1) (2008) 7–20. G. Fahua, M.R. Anderson, Cfd-based throughflow solver in turbomachinery design systems, in: ASME Turbo Expo 2007, ASME, Montreal, 2007. J.F. Simon, O. Leonard, Modeling of 3-d losses and deviations in a throughflow analysis tool, J.Therm. Sci. 16 (3) (2007) 208–214. R.H. Aungier, Turbine Aerodynamics: Axial-Flow and Radial-Inflow Turbine Design and Analysis, ASME Press, 2006. T.P. Moffit, E.M. Szanca, W.J. Whitney, F.P. Behning, Design and cold-air test of single uncooled core turbine with high work output, Tech. Rep. NASA Technical paper 1680, Lewis Research Center, Cleveland, Ohio (May 1980). [6] T. Arts, Calculation of the three-dimensional, steady, inviscid flow in a transonic axial turbine stage, in: 29th International Gas Turbine Conference and Exhibit Amsterdam, ASME, 1984. paper 84-GT-76. [7] M.S. Liou, C.J. Steffen, A new flux splitting scheme, J. Comput. Phys. 107 (1) (1993) 23–39. [8] J. Forˇt, J. Fürst, J. Halama, K. Kozel, Numerical simulation of 3D transonic flow through cascades, Math. Bohem. 126 (2) (2001) 353–361.
Please cite this article in press as: J. Fürst et al., Numerical simulation of circumferentially averaged flow in a turbine, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.01.086