Acta Polytechnica Vol. 52 No. 6/2012 Transonic Flow of Wet Steam — Numerical Simulation Jan Halama Department of Technical Mathematics, Faculty of Mechanical Engineering, Czech Technical University in Prague, Karlovo nám. 13, 121 35, Prague Corresponding author: jan.halama@fs.cvut.cz Abstract The paper presents a numerical simulation of the transonic flow of steam with a non-equilibrium phase change. The flow of steam is approximated by a mixture model complemented by transport equations for moments. Proper formulation of the problem consists of domain definition, a complete set of equations, and appropriate choice of initial and boundary conditions. This problem is then solved numerically by a numerical code, that has been developed in-house. The code is based on a fractional step method and a finite volume formulation. Important issues related to numerical solution are discussed. Results for flow in a turbine are presented. Keywords: two-phase flow, fractional step method, homogeneous nucleation. 1 Introduction A flow in a nozzle or a turbine cascade is typically accelerated and thus subjected to a pressure and tem- perature drop, which can in the case of steam flow initiate condensation. This pressure and temperature drop is very fast in the case of transonic flow and condensation starts later at a certain sub-cooling (a typical sub-cooling is around 30 to 40K below the saturation temperature). The condensation in the turbine decreases the thermal efficiency and causes erosion of the blade. Experimental tests of steam flow are demanding. This fact motivates the development of numerical methods. We consider here a compress- ible flow of wet steam, which is a mixture of vapor and condensed droplets. The flow of the mixture is approximated by inviscid or laminar flow models. We further consider homogeneous condensation, simple convection of droplets by vapor, low-level wetness (negligible volume of droplets) and common pressure for both vapor and droplets. 2 Flow model The flow of a mixture is described by the Euler or Navier-Stokes equations, which are complemented by A B Γ Γ Γ Γ Γ Γi p p w w o Figure 1: Example of the computational domain. the transport equations for moments of liquid phase. The complete set of equations reads ∂tW + ∂x(Fc − Fv) + ∂y(Gc −Gv) = Q, (1) where W =  ρ ρux ρuy e ρw ρQ2 ρQ1 ρQ0  , Q =  0 0 0 0 4 3πr 3 cρlJ + 4ρπQ2ṙρl r2 cJ + 2ρQ1ṙ rcJ + ρQ0ṙ J  , Fc =  ρux ρu2 x + p ρuxuy (e+ p)ux ρwux ρQ2ux ρQ1ux ρQ0ux  , Fv =  0 τxx τxy uxτxx + uyτxy − qx 0 0 0 0  , Gc =  ρux ρuyux ρu2 y + p (e+ p)uy ρwuy ρQ2uy ρQ1uy ρQ0uy  , Gv =  0 τxy τyy uxτxy + uyτyy − qy 0 0 0 0  124 Acta Polytechnica Vol. 52 No. 6/2012 (a) inviscid flow model (b) laminar flow model Figure 2: Mach number isolines, single-phase flow case. and τxx = η 4 3∂xux − 2 3η∂yuy, τxy = τyx = η (∂yux + ∂xuy) , τyy = η 4 3∂yuy − 2 3η∂xux, qx = −λ∂xT, qy = −λ∂yT, (2) where ρ denotes mixture density, ux and uy mixture velocity components, e mixture total energy per unit volume, w mass fraction of liquid phase (wetness), τ shear stress, q heat flux, η viscosity, T tempera- ture. The droplet size spectrum is described by three moments [7] Q0 = N, Q1 = N∑ i=1 ri, Q2 = N∑ i=1 r2 i , (3) where N denotes the total number of droplets per unit mass of mixture, and ri is the radius of the i-th droplet. The average droplet radius r is defined as r = { 0 if w ≤ 10−6,√ Q2/Q0 if w > 10−6, (4) where the limit value 10−6 for the wetness is set to avoid division by a number close to zero. The system of equations is closed by the equation for the pressure according to [14] or [11]: p = (γ − 1)(1− w) 1 + w(γ − 1) [ e− 1 2ρ(u 2 x + u2 y) + ρwL ] , (5) where the specific heat ratio γ is considered as a local function of temperature. The condensation process consists of two differ- ent phenomena. The first is the appearance of new droplets (nucleation), when the vapor temperature drops sufficiently below the saturation temperature. The number of new droplets per unit volume and per second is approximated by term (6) according to [2]: J = √ 2σ πm3 v · ρ 2 v ρl · exp ( −β · 4πr2 cσ 3kBTv ) , (6) where σ denotes the water surface tension, mv is the vapor molecule mass, ρv = (1−w)ρ the vapor density, ρl the liquid density, kB the Boltzmann constant, Tv the vapor temperature, β the surface tension correc- tion coefficient according to [10] (β = 1.328p0.3 cor±0.05, where pcor [bar] denotes the pressure at the intersec- tion of the isentropic expansion from reservoir con- ditions with the steam saturation line) and rc the critical radius, which is rc = 2σ ρlRvTv ln(pv/ps) . (7) 125 Acta Polytechnica Vol. 52 No. 6/2012 (a) inviscid flow model (b) laminar flow model Figure 3: Mach number isolines, two-phase flow case. The vapor pressure pv is considered equal to the pres- sure of the mixture from Eq. (5). The second phenomenon is the growth of an existing droplet, which is approximated by the time derivative of the radius of the droplet, see [17]: ṙ = λv(Ts − Tv) Lρl(1 + 3.18 ·Kn) · r − rc r2 [m s−1], Kn = νv · √ 2πRvTv 4rpv [1]. (8) Term (8) also takes into account the evaporation. If Tv > Ts we set rc = 0. 3 Problem formulation Consider the integral form ∂t ∫∫ V W dV = − ∫∫ V (∂xF + ∂yG−Q) dV (9) of system (1) for any subset V ⊂ D, where D ⊂ R2 is the solution domain. The solution W : D → R8 has to fulfill the integral form (9) for any subset V ⊂ D and proper boundary conditions along boundary ∂D. The boundary consists of the following parts: inlet Γi, outlet Γo, non-permeable wall Ωw and periodical boundary Γp (we expect periodicity in the vertical direction), see the example in Fig. 1. We use the following boundary conditions for the inviscid flow model. The velocity component, which is perpendicular to the inlet boundary, is subsonic for all considered flow cases. Therefore according to the 1D theory of characteristics ‘the number of unknowns minus one’ parameters have to be set. We prescribe constant values of reservoir conditions T0 and p0, the flow direction, w = 0, Q0 = 0, Q1 = 0 and Q2 = 0. The non-permeability condition (ux, uy)~n = 0, where ~n denotes the unit vector normal to the boundary, is considered along the walls. The velocity component, which is perpendicular to the outlet boundary, is also subsonic for all considered cases, i.e. according to the 1D theory of characteristics one parameter has to be specified — we prescribe a mean outlet static pressure value. The computational domain for turbine flow calculations consists of one blade passage, i.e. we consider the spatial periodicity of the solution. The boundary conditions for viscous flow are similar. The system of the above mentioned conditions is sup- plemented by proper Neumann’s conditions. The non- permeability condition along the walls is, of course, replaced by the no-slip condition ux = 0, uy = 0 and 126 Acta Polytechnica Vol. 52 No. 6/2012 (a) inviscid flow model (b) laminar flow model Figure 4: Wetness isolines, two-phase flow case. the adiabatic wall condition ∂T/∂~n = 0. 4 Numerical method The properties of the source term make time integra- tion difficult. The nucleation rate can change very rapidly with respect to space as well as time vari- ables. Numerical simulation thus requires fine spatial discretization in the nucleation zone and sufficiently small time step to achieve accurate prediction of the number of new droplets. The required time step can be one or two orders smaller than the time step appropriate for convection. We therefore separate convection and nucleation phenomena by a splitting method based on the symmetrical (or Strang) split- ting method, see [15], i.e. we solve equations (i)–(iii) successively, instead of directly solving the system (1) ∂tW = Q, (i) ∂tW = −∂xFc + ∂xFv − ∂yGc + ∂yGv, (ii) ∂tW = Q (iii) A finite volume method is used to solve the convection-diffusion part (ii) and the explicit two- stage Runge-Kutta method is applied in each grid point to compute condensation parts (i) and (iii). Let us denote one step of a finite volume method with ‘initial data’ Wn i,j and time step ∆t by the sym- bol FV(Wn i,j ,∆t) and one step of the Runge-Kutta method as RK(Wn i,j ,∆t). Then one step of the full algorithm reads W(0) i,j = Wn i,j , W(k+1) i,j = RK ( W(k) i,j , ∆t 2N ) , k = 0, . . . ,N − 1, W(N+1) i,j = FV ( W(N ) i,j ,∆t ) , (10) W(k+1) i,j = RK ( W(k) i,j , ∆t 2N ) , k = N+1, . . . , 2N , Wn+1 i,j = W(2N+1) i,j , with N = ∆t/τ , where the step ∆t should satisfy the stability condition of the finite volume method for equation (ii), and step τ should be small enough to ensure sufficient accuracy of the source term integra- tion. Step ∆t is used globally (the same value for all finite volumes), while step τ is used locally (smaller steps are applied especially within the nucleation zone, where the nucleation rate changes its magnitude by many orders). As initial conditions W0 i,j we usually prescribe the converged steady solution of the problem with the same geometry and boundary conditions, but without the source term (condensation is ‘switched off’). 127 Acta Polytechnica Vol. 52 No. 6/2012 (a) inviscid flow model (b) laminar flow model Figure 5: Average droplet radius, two-phase flow case. The current numerical algorithm with several fi- nite volume methods (Lax-Wendroff, CFFV method, SRNH method) has been successfully validated for the case of inviscid transonic flow in a convergent- divergent nozzle, for details see [6]. 5 Results and Discussion Figures 2–7 show numerical results for two-dimen- sional two-phase flow in turbine cascade SE1050 ob- tained by the inviscid and laminar (Re = 1.5·106) flow models and the Lax-Wendroff finite volume method for the convection-diffusion part. The geometry of turbine cascade SE1050 is quite often used by many authors, and various numerical results, mainly for air flow, are available in the literature, see e.g. [4] or [3]. The geometry of the SE1050 cascade is provided by the QNET network, or it can be found e.g. in [14]. Although the laminar flow model used for high Reynolds number flow is an artificial case, in to our experience inviscid, laminar and turbulent models for air flow yield a similar pressure field. We consider two cases: the single-phase flow, which omits the source term Q in Eq. (1), i.e. it takes into account only convection and diffusion without condensation, and two-phase flow, which takes into account the full set of equations (1) with the source term Q. A com- parison between those two cases shows the effect of the addition of latent heat to the flow. The results of the single-phase flow case are used as the initial data for the computation of two-phase flow. We consider the same boundary conditions for all computed cases. The inlet flow angle is equal to 19.3°from the axial direction, the inlet total pressure p01 = 36730 Pa, the inlet total temperature T01 = 340 K and the mean outlet pressure poutlet/p01 = 0.423. Figure 2 shows the results in the form of Mach number isolines. The results for single-phase inviscid and single-phase laminar flow models have a similar structure due to the very thin boundary layer for the laminar flow model. The difference is mainly in reflection of the right running trailing edge shock wave. The latent heat released by condensation slows down the supersonic flow and it changes the expansion (different position of the throat, the structure of the shock wave). The wetness isolines in Fig. 4 show a higher gradi- ent in the nucleation zone for the inviscid flow model. The isolines of the average droplet radius are plotted in Fig. 5. The inviscid model predicts bigger droplets. The total number of droplets per unit volume is given in Fig. 7. The inviscid model yields a smaller num- ber of larger-size droplets, unlike the larger number 128 Acta Polytechnica Vol. 52 No. 6/2012 (a) inviscid flow model (b) laminar flow model Figure 6: Vapor sub-cooling, two-phase flow case. of smaller droplets for the laminar flow model. The total amount of liquid phase (wetness) at the outlet is similar for both inviscid and laminar two-phase flow models, and it approaches the equilibrium wetness at the domain outlet. The isolines of sub-cooling in Fig. 6 also show faster nucleation for the inviscid flow model. 6 Conclusions The physical model considered here takes into account condensation as well as evaporation of the droplets. Numerical tests have shown that the proposed numer- ical method has sufficient robustness. Strictly local character of the integration of the source term en- ables the strictly local refinement of the time steps for condensation (i.e. more cycles of the explicit Runge- Kutta 2-stage method). This means that N differs from point to point. The computation of two-phase flow usually takes twice as much CPU time as one- phase flow computation for the same grid and equiv- alent boundary conditions. The method has been validated for the Barschdorff nozzle [6]. The numeri- cal results presented here for steam flow in the SE1050 turbine cascade are physically expectable, and show the ability of our method to model condensation and evaporation in more complex flow fields. Acknowledgements This work has been supported by MSM CR Research Plan No. 212200009. References [1] D. Barschdorff. Verlauf der Zustandgroessen und gasdynamische Zuammenhaenge der spontanen Kondensation reinen Wasserdampfes in Lavaldue- sen. Forsch. Ing.-Wes., 37(5), 1971. [2] R. Becker, W. & Döring. Kinetische Behand- lung der Keimbildung in übersättingten Dämpfen. Ann. d. Physik, 24(8), 1935. [3] V. Dolejší. Anisotropic mesh adaptation tech- nique for viscous flow simulation. East-West Journal of Numerical Mathematics, 9(1):1–24, 2001. [4] V. Dolejší, M. Feistauer, J. Felcman, A. Kliková. Error estimates for barycentric finite volumes combined with nonconforming finite elements ap- plied to nonlinear convection-diffusion problems. Aplications of Mathematics, 47(4):301–340, 2002. [5] S. Dykas, K. Goodheart, G. H. Schnerr. Nu- merical study of accurate and efficient modelling 129 Acta Polytechnica Vol. 52 No. 6/2012 (a) inviscid flow model (b) laminar flow model Figure 7: Isolines of log10 Q0, two-phase flow case. for simulation of condensing flow in transonic steam turbines. In 5th European conference on Turbomachinery, Prague, 2003, pp.751–760. [6] J. Halama, F. Benkhaldoun, J. Fořt. Flux schemes based finite volume method for inter- nal transonic flow with condensation. Interna- tional Journal for Numerical Methods in Fluids, 65(8):953–968, 2011. [7] P. G. Hill. Condensation of water vapor during supersonic expansion in nozzles, part 3. Journal of Fluid Mechanics, 3:593–620, 1966. [8] R. J. LeVeque. Numerical methods for conserva- tion laws. Birkhäuser, 1999. [9] R. H. Ni. A multiple grid scheme for solving Euler equations, AIAA Journal 20(1), 1981. [10] V. Petr, M. Kolovratník. Heterogeneous Effects in the Droplet Nucleation Process in LP Steam Turbines. Proceedings of 4th European Confer- ence on Turbomachinery, Firenze, Italy, 2001. [11] M. Šejna, M. J. Lain. Numerical modelling of wet steam flow with homogeneous condensation on unstructured triangular meshes. Journal ZAMM, 74(5):T375–T378, 1994. [12] P. Sopuch. Kinetics of phase change vapor-liquid and its numerical simulation. Doctoral thesis, IT CAS CR, Prague, 1996 (in Czech). [13] M. Šťastný, P. Šafařík. Boundary layer effects on the transonic flow in a straight turbine cascade. ASME Paper 92-GT-155, 1992. [14] M. Šťastný, M. Šejna. Condensation effects in transonic flow through turbine cascade. Proceed- ings of The 12th international conference on the properties of water and steam, New York, Begel House, 2005, pp.711–719. [15] G.Strang. On the construction and comparison of difference schemes SIAM Journal of Numerical Analysis, 5:506–517, 1968. [16] S. M. Stringer, K. W. Morton. Artificial viscosity for the cell vertex method. Report no. 96/08, Oxford University Computing Laboratory, 1996. [17] J. Valha. The flow of wet steam in the through- flow part of a steam turbine. Doctoral thesis, SVUSS Běchovice, 1988, (in Czech). [18] J. B. Young, M. Moheban. A time marching method for the calculation of blade-to-blade non- equilibrium wet steam flows in turbine cascades. I. Mech. E. C76/84:89–99, 1984. 130