Transonic Flow of Wet Steam — Numerical Simulation

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.


Introduction
A flow in a nozzle or a turbine cascade is typically accelerated and thus subjected to a pressure and temperature 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 40 K 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 compressible 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.

Flow model
The flow of a mixture is described by the Euler or Navier-Stokes equations, which are complemented by the transport equations for moments of liquid phase.The complete set of equations reads where  and where ρ denotes mixture density, u x and u y mixture velocity components, e mixture total energy per unit volume, w mass fraction of liquid phase (wetness), τ shear stress, q heat flux, η viscosity, T temperature.The droplet size spectrum is described by three moments [7] where N denotes the total number of droplets per unit mass of mixture, and r i is the radius of the i-th droplet.The average droplet radius r is defined as 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]: where the specific heat ratio γ is considered as a local function of temperature.The condensation process consists of two different 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]: where σ denotes the water surface tension, m v is the vapor molecule mass, ρ v = (1 − w)ρ the vapor density, ρ l the liquid density, k B the Boltzmann constant, T v the vapor temperature, β the surface tension correction coefficient according to [10] (β = 1.328p 0.3 cor ± 0.05, where p cor [bar] denotes the pressure at the intersection of the isentropic expansion from reservoir conditions with the steam saturation line) and r c the critical radius, which is  The vapor pressure p v is considered equal to the pressure 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]: Term ( 8) also takes into account the evaporation.If T v > T s we set r c = 0.

Problem formulation
Consider the integral form of system (1) for any subset V ⊂ D, where D ⊂ R 2 is the solution domain.The solution W : D → R 8 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 T 0 and p 0 , the flow direction, w = 0, Q 0 = 0, Q 1 = 0 and Q 2 = 0.The non-permeability condition (u x , u y ) 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 supplemented by proper Neumann's conditions.The nonpermeability condition along the walls is, of course, replaced by the no-slip condition u x = 0, u y = 0 and

Numerical method
The properties of the source term make time integration difficult.The nucleation rate can change very rapidly with respect to space as well as time variables.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) splitting method, see [15], i.e. we solve equations (i)-(iii) successively, instead of directly solving the system (1) A finite volume method is used to solve the convection-diffusion part (ii) and the explicit twostage 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' W n i,j and time step ∆t by the symbol FV(W n i,j , ∆t) and one step of the Runge-Kutta method as RK(W n i,j , ∆t).Then one step of the full algorithm reads , 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 integration.
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 W 0 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').The current numerical algorithm with several finite volume methods (Lax-Wendroff, CFFV method, SRNH method) has been successfully validated for the case of inviscid transonic flow in a convergentdivergent nozzle, for details see [6].
Figures 2-7 show numerical results for two-dimensional two-phase flow in turbine cascade SE1050 obtained by the inviscid and laminar (Re = 1.5•10 6 ) 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 comparison 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 p 01 = 36730 P a, the inlet total temperature T 01 = 340 K and the mean outlet pressure p outlet /p 01 = 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 gradient 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 number of larger-size droplets, unlike the larger number

Conclusions
The physical model considered here takes into account condensation as well as evaporation of the droplets.Numerical tests have shown that the proposed numerical method has sufficient robustness.Strictly local character of the integration of the source term enables 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 onephase flow computation for the same grid and equivalent boundary conditions.The method has been validated for the Barschdorff nozzle [6].The numerical 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.

Figure 1 :
Figure 1: Example of the computational domain.
a) inviscid flow model (b) laminar flow model
(a) inviscid flow model (b) laminar flow model
(a) inviscid flow model (b) laminar flow model

Figure 6 :
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.
(a) inviscid flow model (b) laminar flow model