Numerical Solution of Transonic Wet Steam Flow in Blade-to-Blade Cascade with Non-equilibrium Condensation and Real Thermodynamics

We present an engineering approach to mathematical modeling and numerical solution of 2D invis- cid transonic flow of wet steam in a steam turbine cascade channel of penultimate stage at rotor tip section in full Eulerian framework. Our flow model consists of the Euler system for the mixture (dry steam + homogeneously dispersed water droplets) and transport equations for moments of droplet number distribution function known as method of moments. Thermodynamic properties of vapor steam are provided by set of IAPWS equations. For equation of state for vapor phase valid both in superheated and wet (meta-stable) region we adopted recently developed equation in CFD formulation for low pressures provi1ded by Hrub´ y et al. (9), (8), (10). For extrac- tion of vapor parameters from the mixture ones we implemented simple relations in polynomial form describing thermodynamic properties of saturated liquid state. Nucleation model is resorting to modified classical nucle- ation theory. Linear droplet growth model is implemented for calculation of liquid sources. Numerical method is simple: cell-centered finite volume approach, 1st-order AUSM+ scheme for spatial derivatives, symmetrical fractional step method for separation of convection and condensation part, explicit 2-stage 2nd-order Runge- Kutta method for time integration. Geometry of blade profile and experimental results are provided by Bakhtar's work (22), (23). Results were obtained for one subsonic inlet/subsonic outlet regime and gave quite reasonable accordance with experiment.


Introduction
By the wet steam flow we mean just one of the various regimes of multiphase flow -the flow in which a large number of tiny spherical liquid droplets of various size are homogeneously dispersed in a continuous gas phase consisting only of pure vapor. Steam turbines are mechanical devices that transfer thermal energy from expanding steam to electrical power. It is expected that they will continue to play their prominent role in the ensuing decades in this respect. Nucleation process is accompanied with the irreversible (latent) heat release which results in additional total pressure loss and influences thermodynamic and flow parameters. Rapid expansion of initially superheated dry steam leads to the crossing of equilibrium saturation curve, but the first droplets occur only after certain supercooling is attained (non-equilibrium process). Physical theory, experimental measurements and numerical simulations are more complicated in the presence of the second phase. In practice many nucleation models rely on almost 80 years old Becker-Döring classical nucleation theory which with some added corrections gives in case of steam quite reasonable results. Accurate simulations require real thermodynamic concept even if low pressure conditions prevail, as condensation models are very sensitive to the thermodynamic parameters. Our approach is based on the homogeneous mixture non-equilibrium model which approximates pressure of the mixture by the vapor pressure and neglects slip between phases. a e-mail: vladimir.hric@fs.cvut.cz b e-mail: jan.halama@fs.cvut.cz 2 Mathematical models

Mixture model
Equations describing dynamics of fluid flow are based on conservation laws of mass, momentum and total energy in Eulerian frame of reference. We consider that there is no velocity slip between phases and that the phases share the same pressure field. Therefore, instead of balancing conservative properties for both phases separately -so-called two-fluid model [7], we treat both phases as a one mixture -so-called mixture (or single-fluid) model. Detailed assessment of single-and two-fluid model can be found e.g. in [6]. Provided that viscosity, heat conduction and gravitational effects are negligible, governing equations for mixture dynamics can be expressed in the form of inviscid Euler system of equations where ρ denotes mixture mass density, u is velocity vector, E is total energy of mixture, p is static pressure, H = E + p/ρ is total mixture enthalpy and I is identity matrix. The system needs closure in the form of vapor phase equation of state for pressure. Benefit of the mixture approach is reduction of computational demands, but a main drawback is that vapor properties are not directly accessible from the mixture ones. For this purpose we extract vapor density and vapor internal energy in the following way Further, we can consider thermodynamic states of liquid phase to be close to saturation liquid boundary. Therefore, we implemented simple relations in polynomial form with temperature as independent variable for liquid properties like density ρ , enthalpy h and speed of sound w , such that they fit liquid saturation data. Specific internal energy of liquid phase e is considered to comprise bulk part only (energy of surface is not accepted). More elegant and thermodynamically consistent method is outlined in [9]. It is assumed that at lower pressures (p < 1 bar) internal energy stored in liquid phase differs not too much from enthalpy of liquid phase, hence we adopted simplification e ≈ h .

Equation of state for vapor phase
Most of CFD codes still uses only perfect gas relation for pressure p = ρRT and internal energy e = RT/(γ − 1), sometimes with modified specific gas constant R and ratio of specific isobaric to specific isochoric heat capacity κ = c p /c v . However, modeling of thermodynamic states in meta-stable region requires using more general equation of state. State-of-the-art equation of state for steam valid both in vapor and meta-stable region provides IAPWS in so-called fundamental formulation either by relation for dimensionless Helmholtz free energy φ = φ(δ, τ) with 16 terms in total [17], [18] or relation for dimensionless Gibbs free enthalpy γ = γ(π, τ) with 23 terms in total [4], [19]. The main disadvantage of IAPWS equations of state is the fact that they require solution of one or two nonlinear algebraic equations in order to calculate pressure from density and internal energy, which prolongs calculation time. Therefore, at the Institute of thermomechanics ASCR developed so-called CFD formulation of equation of state [8] which provides function for dimensionless entropy with independent variables scaled density δ and scaled internal energy ξ η = η(δ, ξ) = η • (δ, ξ) + η r (δ, ξ) (5) in which the ideal part of entropy η • and the residual part of entropy η r are given by relations with constants c i and b 1,2,3 . Scaling of independent variables is defined as follows with constants β 1,2 and constant critical parameters of steam denoted by subscript (·) c . More details and percentage deviations from experimental data can be found in [8], [10].
Here we give only derived relation for pressure, temperature and speed of sound (for vapor phase) where coefficients are given as follows where (·) δ,ξ mean partial derivatives. The foregoing equation of state in CFD formulation is applicable only for lower pressures. The low pressure boundary is determined by the value of entropy s ≥ 7.5 kJ kg −1 K −1 . Below the entropy limit, deviations from experimental data grow significantly.

Evolution of liquid phase
Evolution of liquid phase can be described by mass conservation law for droplets in phase space (x, r). Instead of expressing balance equation in the general discrete kinetic form which takes interaction between all monomers (clusters containing one molecule) into account, we use its reduced continuum version for balance of droplet number continuous distribution function f = f (x, r, t). Then, the total number of droplets in unit volume with radii smaller thanr at the location x and at the time t is given by integral r 0 f (x, r, t) dr. For much more details about droplet size distributions we refer to [1]. This reduction assumes EPJ Web of Conferences 02025-p.2 that droplets with very small number of molecules are neglected. Then, balance equation takes the form The first term on the RHS is related to the droplet growth as a result of capturing new molecules from the surrounding vapor and the second source concerns creating of new droplets at the critical radius, δ is Dirac delta function. Denucleation source and droplet growth owing to liquid compressibility is not included. Coagulation and segregation of molecules are not respected as well. Because the distribution function f is not known a priori and we do not intend to model its shape, we reduce the balance equation further by its multiplication by r k and subsequent integration ∞ 0 (·) dr. With the definition of the k-th moment of the droplet number distribution function and applying integration by parts in droplet growth source term with boundary conditions we get the efficient so-called moment version of balance equation for moments In our CFD code we utilized only the first four moments, i.e. k = 0,1,2,3. This is sufficient in practical applications for basic modeling of polydispersed droplet spectra. In the classical method of moments droplet growth law (equation forṙ) cannot be used in the general form with respect to radial coordinate r. We decided to provide its linear approximation in the vicinity of average droplet radius, hencė r = a + br, where constants a, b have to be determined properly and are functions of vapor properties only. For average droplet radius we used Hill's definitionr = μ 2 /μ 0 . With the mentioned closure method we can rewrite moment equations as follows where k = 0 : 3. The last equation (k = 3) can be transformed to the equation for wetness factor by its multiplying by factor 4πρ /3 and using definition for wetness fraction ρy = 4πρ μ 3 /3 and utilizing ρ = const.
Moments of droplet number distribution function μ k are occasionally expressed by means of so-called liquid moments Q k and holds

Creating of new water droplets
Rapid expansion of initially superheated steam into the metastable region causes its supercooling, but steam remains dry until from the energetic point of view it is more advantageous to start nucleation process. Driving force for nucleation is so-called supersaturation ratio S = p/p sat or equivalently supercooling ΔT = T sat − T v . Condensation (phase change) process-initiating work can be expressed by means of thermodynamic potential Gibbs free enthalpy G. The change in free enthalpy during condensation process that poses the minimal work required to form a spherical water droplet is derived for general equation of state in [3] as follows where σ is surface tension of water and is modeled by IAPWS release [20], which allows its applicability for temperatures below 273.15K (possibility of its variation with droplet radius is not considered), r is droplet radius, R is specific gas constant of water, is difference in residual part of specific free enthalpy between actual and saturated state at constant vapor temperature for which we utilized IAPWS-97 formulation [19]. Then, we can write where τ = T c /T , π = p (MPa) and n i , I i , J i are given constants. The radius at which the change in free enthalpy attains maximum is called critical radius r c (system is in unstable equilibrium) and it results from applying condition for local extremum of function ΔG with respect to radius Clusters with radii larger than critical tend to grow contrary to clusters with radii smaller than critical which tend to decay (system tends to minimize its free enthalpy). Nucleation rate J is the rate at which number of new cores of droplets (nuclei) in unit volume is spontaneously generated from supersaturated vapor steam at the onset of homogeneous condensation. Because we consider steam to be perfectly pure, no foreign particles exist in the flow field an therefore we neglect heterogeneous condensation altogether. Model for nucleation rate is for engineering applications frequently provided by classical kinetic nucleation theory which assumes capillarity approximation and works with macroscopic quantities. Classical steady-state nucleation rate of Becker and Döring is expressed as where so-called nucleation barrier ΔG c = 4πr 2 c σ/3 is given by substituting critical radius (27) into Eq. (25). Group k B T v represents thermal energy with k B Boltzmann constant, k B = 1.3806488e −23 J K −1 and K is kinetic prefactor in which m 1 = M/N A is mass of single steam molecule, M = 18.015268e −3 kg mol −1 is molecular mass of water and N A = 6.02214129e 23 mol −1 is Avogadro constant. Classical nucleation rate J cl is frequently corrected by Courtney's and Kantrowicz's correction factors which take into account effects of partial pressure of the clusters and temperature difference between clusters and vapor phase. According [21] joint Courtney's and Kantrowicz's correction is where α c is convective heat transfer coefficient for critical cluster and is defined later, Eq. (38). Then, corrected nucleation rate is given by

Growth of existing water droplets
Temperature of water droplets differs from surrounding vapor phase temperature the more the faster expansion takes place. Faster expansion imply faster droplet growth, hence there is less time to diminish temperature difference. The rate of change in radial direction of water droplet can be determined by balance of energy over a surface of spherical droplet (absence of mass transfer). Only a part of released latent heat by condensing vapor steam is used for droplet heating. Simplified energy balance equation can be written as follows where the term on LHS is the rate of local specific latent heat released to the droplet and the RHS term is a portion of heat flowing toward vapor by convection mechanism between droplet and surrounding vapor steam. Heat transfer by conduction throughout droplet is neglected. Substituting m = 4πρ r 3 /3 for mass and A = 4πr 2 for surface of single droplet we arrive at the default equation for droplet growth Important quantity in determining droplet growth is Knudsen number Kn defined as a ratio of mean free path of vapor molecules¯ v to the droplet diameter where˜ v is given by kinetic theory. Regime with Kn 1 (Kn > 4.5) is considered as free molecular and with Kn 1 (Kn < 0.01) as continuum, while in between there is the transition regime. Since size of droplets in steam turbines covers quite wide range, providing droplet growth model covering entire range of Knudsen number is useful from practical point of view.
Gyarmathy postulated that temperature of sub-micron water droplets can be estimated utilizing capillarity approximation in the following way Hence, critical droplets have temperature of surrounding vapor and the more the droplets are larger then their critical size the closer is their temperature to saturation temperature. For the convective heat transfer coefficient he presented empirical relation [5]. Here we adopted its slightly modified form by utilizing Langmuir droplet model which defines interface radius r i = r + β˜ v . We distinguish following regimes r ≤ r i . . . continuum regime, r > r i . . . f ree molecular regime.
Choice β = 0.75 is frequently suggested. Gyarmathy-Langmuir heat transfer coefficient takes the form We rewrote it as follows by utilizing relation for Prandtl number Pr = μc p /λ and isobaric specific heat capacity for perfect gas c p = κR/(κ − 1). Substituting κ = 1 + 1/n • 3 ≈ 1.332633 from [18] we arrive at α = 1 where c = 1.032663666979 e −3 . In the continuum regimẽ r, relation (38) is in accordance with conduction theory for which α = λ/r. For our implementation we used only temperature dependent ideal part of vapor thermal conductivity λ v = λ v (T v ) and the same holds for vapor phase dynamic viscosity μ v = μ v (T v ), both stated in [4]. With Gyarmathy and Langmuir assumptions final form of droplet growth takes the forṁ where A, B depend only on vapor parameters Another (probably better) growth model presents e.g. White and Young in [11] and Dykas and Wróblewski in [7]. Coefficients of linear approximation of droplet growth g(r) ≈ a + br are given as follows a =ṙ(r) −ṙ (r)r, b =ṙ (r), whereṙ (r) is radial derivative of droplet growth evaluated at average radiusr and in our CFD code is implemented exactly.

Final form of governing equations
Following set of partial differential equations in conservative form is solved where indexes i = 1 : 4 form Euler equations for mixture and indexes i = 5 : 8 form equations for liquid phase, W i are conservative variables, F i are components of flux vector in x-direction, G i are components of flux vector in y-direction and S i are liquid phase sources. In more detail for any subsetD ⊂ D and boundary conditions prescribed at boundaries of solution domain.

Calculation of vapor parameters from conservative parameters
In order to calculate conservative vector W n+1 in new time level (n + 1) we have to determine static pressure, vapor temperature and sound speed of mixture from the conservative vector in previous time level W n . Pressure and vapor temperature are calculated from vapor density and vapor specific internal energy utilizing equation of state for but in our case it is calculated with frozen composition, which results in vanishing the last term in denominator. Note that we calculated liquid density, liquid specific internal energy, cf. (4), and liquid sound speed from vapor temperature in previous time level T n v , otherwise we would had to solve non-linear algebraic equations simultaneously.

Boundary conditions for subsonic inlet and subsonic outlet
Due to equation of state in CFD formulation, no iteration is needed in order to calculate static pressure and/or vapor temperature in inner cells, but there is different situation in case of fluxes calculations through inlet and outlet boundary. Inlet ghost cells work with first order extrapolation of static pressure p from the first neighboring cell entropy s and enthalpy h equal to their stagnant values and calculated from prescribed total pressure p stag and total temperature T stag during preprocessing stage prescribed flow angle α assumption of dry steam: W 4 = W 5 = W 6 = W 7 = 0 Therefore, numerical solution of system of two nonlinear algebraic equations (from given pressure and entropy calculate density and temperature) is necessary provided we intend to keep the same state equation. Newton-Raphson method is implemented with appropriate stop criterion. Properties of vapor phase are equal to their mixture values, as we consider dry steam at inlet. Similar situation is in outlet ghost cells, where we solve one nonlinear equation -vapor temperature is calculated from prescribed average static back pressure and vapor density which is extracted from first order extrapolation of mixture density and wetness factor from the first neighboring cell.
In the inviscid case we implemented slip condition on the wall boundaries. Normal component of velocity vector is zero, hence fluxes are composed from pressure part only. Static pressure on the wall is approximated by the value of static pressure in the first neighboring cell.
On periodic boundaries we implemented point-to-point periodicity.

Solution procedure
In order to tackle stiff character of the system due to different time scales of convection and condensation, we split the system into two subsystems by fractional step method. Because we implemented this method in symmetrical form [12], conservative variables in new time level (n+1) are calculated from the ones in old time level (n) by following successive steps

Time step
Δt in inviscid case is estimated as in [2] and we used CFL = 0.75. Time integration in case of condensation part is repeated k-times with reduced time step, where k is given as and condensation time scale τ c is estimated by Wegener formula [14] τ Convective fluxes of conservative variables through cell interfaces are calculated applying mid-point rule and AUSM+ scheme [16]. Only constant distribution of variables in cell volume is considered, that is why only 1st order spatial accuracy is attained. The scheme uses arithmetic averaged speed of sound in mixture w from left and right cells. For time integration we implemented explicit 2nd-order 2-stage Runge-Kutta method in case of both parts.

Calculation of liquid sources
We distinguished four cases in calculation of liquid sources S. Our algorithm is as follows: Calculate supersaturation ratio S = p/p sat (T v ) and average droplet radiusr where we decided for W 7 min = 1 e −12 . Then proceed with the following branching We set r min = 2 e −10 m.

Test case and remarks
We followed the validation test case no. 27 in [22], [23]. So, we prescribed at inlet boundary p stag = 0.0999 MPa, T stag = 360.8 K, α = 38 • , and at outlet boundary average pressure ratio p stag /p out = 2.34.
Only one period of blade cascade was meshed. We used unstructured mesh configuration with quadrangle elements. Free tool Gmsh [13] was used for this purpose.  Implementation of linear model of droplet growth gave good results in case of flow in nozzles (with slightly worse convergence comparing to constant droplet growth). However its using in case of blade-to-blade cascades resulted in very different results comparing to experiment data and results obtained by several authors. We blame unrealistic expansion in the vicinity of sharp trailing edge for that. Therefore, we simplified our model by consideration of constant droplet growth (b = 0) in evaluating liquid sources. Figure 4 presents distribution of supercooling ΔT = T sat − T v and figure 5 supersaturation ratio S = p/p sat of expanding steam in the blade cascade with some representative contours of constant value. Saturation contour S = 1 is showed as well. Region with S < 1 represents superheated region where evaporation of droplets occurs. This region is initiated by shock wave behind which vapor temperature crosses saturation temperature. Figure 6 depicts some contours of nucleation rate J. Unphysical strong expansion at the trailing edge is the reason of maximum nucleation rate (J > 1e24 m −3 s −1 , region MAX in the figure). Static pressure is presented in figure. 7. Figure 8 shows Mach number in the domain. Maximum local mach number is 1.28. Contour with M = 1 is presented as well (thick blue line in the figure). Calculated surface-averaged droplet radius in the domain is in figure 9. The biggest droplets with  radius about 0.03 μm are located in the middle of the channel. In the region behind the shock wave droplets are partly shrunk. Wetness factor is in figure 10 and has maximum value about 0.058.

Results with commentary
Along the blade surface we present following: Figures  11, 12 shows supercooling and supersaturation, respectively. Nucleation rate is presented in figure 13. Ratio of static pressure to the total inlet pressure is plotted in figure 14. Figure 15 presents Mach number. Figure 16 and figure 17 plots droplet radius and wetness factor, respectively. The distribution of static pressure on the pressure side is in very good accordance with experiment [22], [23]. The most interesting is the suction part of the blade. First nucleation occurs in the initial part of the suction surface as a result of       Figure 13. Nucleation rate on blade surface curved blade profile. Second nucleation part occurs in the middle of channel, where pressure knee is visible. However, the magnitude of pressure rise is smaller and location is shifted upstream comparing to experiment data. Small location of trailing edge on the pressure side exhibits extrema of nearly all parameters resulting from unphysical over-expansion stemming from using inviscid model.

Conclusion
The main contribution of the paper is implementation of iteration-free vapor equation of state in CFD formulation. We based our work on the full Eulerian approach which combines the Euler inviscid system for the whole mixture and the method of moments for liquid phase. Presented flow model was validated on one selected regime EFM 2014 02025-p.9 [ P :HWQHVV IDFWRU \ Figure 17. Wetness factor on blade surface of transonic wet steam flow in the steam turbine cascade channel of penultimate stage at the rotor tip section provided by Bakhtar. Taking limitations which arise from inviscid approach, our numerical results are in satisfied accordance with the experiment data and numerical simulations provided by the others authors. Our implementation of real thermodynamics of steam was successful. Attention requires implementation of one additional term in liquid sources resulting from linear approximation of droplet growth, which is our goal for future work.