Development and Experimental Validation of Real Fluid Models for CFD Calculation of ORC and Steam Turbine Flows

The article describes an interpolation–analytical method of reconstruction of the IAPWS-95 equations of state and the modified Benedict–Webb–Rubin equations of state with 32 terms (mBWR32). The method enables us to provide the thermodynamic closure in 3D computational fluid dynamics (CFD) calculations of turbomachinery flows with real working media, such as steam and Organic Rankine Cycle (ORC) fluids. The described approach allows for the sufficient accuracy of 3D flow calculations and does not require a significant increase in computational cost over perfect gas calculations. The method is validated against experimental data from measurements and compared with computational results from the model using the Tammann equation of state. Three turbine blading systems are considered—a multi-stage configuration from a low-pressure cylinder of a large-power steam turbine and two ORC microturbines working with organic media HFE7100 and R227ea. The calculation results obtained using the described method of approximation of the IAPWS-95 and mBWR32 equations exhibit satisfactory agreement with the experimental data, considering pressures, temperatures and enthalpies in key sections, as well as turbine power and efficiency in a wide range of changing thermodynamic parameters. In contrast, the Tammann equation of state provides acceptable results only for relatively small changes of thermodynamic parameters.


Introduction
Computational fluid dynamics (CFD) methods are widely used in various fields of science and technology: power engineering, aviation and astronautics, chemical industry, oil and gas industry, etc. A leading approach in CFD for gas dynamics and turbomachinery applications is based on numerical integration of the Reynolds-averaged Navier-Stokes (RANS) equations. In order to close RANS equations and establish a relation between the thermodynamic quantities, equations of state are used, the choice of which depends essentially on the reliability of the obtained computational results.
Currently, the most common equations of state used in 3D calculations are the perfect gas, Tammann, and Van der Waals equations [1][2][3][4]. In many cases, their use is justified and provides acceptable results. A significant increase in the accuracy of the obtained results can be achieved with an individual selection of the equations' constants, as appropriate for the actual range of change of thermodynamic parameters.
However, when the processes occur in a wide range of thermodynamic conditions (as in multi-stage steam turbines or in ORC turbines with a large pressure ratio), or when a phase transition takes place, the simulation accuracy can be significantly reduced [4].
Then, it is necessary to apply more complex equations of state that are polynomials with a large number of terms. Examples of such equations are the IAPWS-95 equations of state [5] and various forms or modifications of the Benedict-Webb-Rubin equations [6][7][8][9]. The IAPWS-95 equations are used to describe the thermodynamic properties of water and steam, while the Benedict-Webb-Rubin equations of state are applied to a wider range of working fluids and are among basic equations of the USA National Institute of Standards and Technology [10]. Unfortunately, direct use of these equations in 3D flow calculations is currently impossible, because in this case, the computing processor time increases by 1.5-2 orders of magnitude.
In response to this question, in the present paper, an interpolation-analytical method to represent thermodynamic functions in complex state equations such as IAPWS-95 equations of state and Benedict-Webb-Rubin equations of state is proposed. The use of this method in 3D CFD calculations allows us to ensure sufficient computational accuracy and, on the other hand, does not require a significant increase in computational costs over perfect gas computations. Here, one of the most accurate modifications of the Benedict-Webb-Rubin thermal equation of state that has 32 terms is used. A method for the determination of constants of the modified Benedict-Webb-Rubin equation of state with 32 terms, on the basis of the available fields of thermodynamic values, is described. These constants of the mBWR32 equation are given for fluids HFE7100 [11] and R227ea [12] working in the investigated ORC microturbines.
Three-dimensional (3D) CFD calculations of the flow parts of a low-pressure cylinder (LPC) of a large power steam turbine and two ORC microturbines are carried out with the help of the software package IPMFlow [13]. The results of CFD calculations obtained with the proposed reconstruction of the IAPWS-95 and mBWR32 equations are compared with the available experimental data from the considered turbines and also with the computational results obtained using the Tammann equation of state. The main purpose of these studies is the validation of the proposed approach of approximating complex thermodynamic equations of state rather than a detailed analysis of the thermodynamic processes and flow phenomena occurring in the investigated turbomachines. In addition, admissible intervals of change of thermodynamic parameters, for which the use of Tammann equation of state can yield reasonable results, are established.

Flow Solver
Three-dimensional (3D) CFD calculations are performed using the IPMFlow software package [13], which is the evolution of earlier software packages FlowER and FlowER-U [14]. It implements the following elements in the mathematical model: Reynoldsaveraged unsteady Navier-Stokes equations [15,16], Menter's k-ω SST two-equation turbulence model [17], the finite volume method, and the implicit quasi-monotone high-order ENO-scheme [18]. Such an approach to provide 3D CFD is typical for calculations of turbine flow parts [19,20]. In most cases for steady flows, the computational domain includes one blade-to-blade channel of each blade row, with periodic boundary conditions. At the boundaries between adjacent blade rows, values of thermodynamic parameters are averaged in the circumferential direction and exchanged. To speed up the convergence process, a quasi-multigrid algorithm [21] and an individual time step are used. Usually, a block-structured H-type grid is used for calculations. The domain of tip leakage is meshed, whereas the domain of shroud leakage is not meshed, and the so-called sink-source approach is used. This approach enables extractions and injections of the fluid streams from and into the main flow domain, so it can be applied to over-the-shroud leakages, inter-disk leakages, and technological extractions [22]. In order to obtain computational results with acceptable accuracy, it was found sufficient to use about 500,000 cells (elementary volumes) in one blade-to-blade channel in each row and a mesh resolution near walls that provided y + values below 2. Every time, final computations were proceeded by grid independence checks.

Equations of State
The equation of perfect gas and the caloric equations of state obtained on its basis [23,24] are most widely known and most often used in 3D CFD calculations. Using the perfect gas equation, it is possible to carry out CFD calculations in flow parts with a relatively short expansion/compression line, such as in a single stage of a high or medium pressure steam turbine cylinder in the superheated steam zone, in one stage of an air compressor of a gas turbine engine, or in a single-stage industrial fan, etc. The constants of the equation of state in such calculations are usually determined by the "average" thermodynamic values between the inlet and outlet.
The Tammann thermal equation of state is quite simple but at the same time more accurate: It differs from the perfect gas equation by the presence of an additional constant. The caloric equations of state obtained on the basis of the thermal equation of state are described in detail in the paper [25]. The extensive experience of the authors has shown that for the calculation of turbomachinery flows, the constants for the Tammann equation should be determined from values of total temperature, pressure and density at the inlet to the flow domain, as well as isentropic static values of temperature, pressure and density at the outlet of the computational domain. Then, the constants R, p 0 , and γ can be determined from Equation (1) and the following relation: The IAPWS-95 equation is nowadays the most accurate equation of state that describes the properties of water and steam. Its detailed description and the corresponding caloric equations of state are presented in [5].
The modified Benedict-Webb-Rubin thermal equation of state with 32 terms has the following form [26][27][28]: The above equation and its corresponding forms have a very complex shape. The caloric equation and equations describing thermodynamic functions corresponding to the BWR32 equation are enclosed and described in detail in Appendix A.

Determination of Constants of the Benedict-Webb-Rubin Equation of State
Usually, constants for equations of state, including mBWR32, are determined based on the experimental data. In the literature, one can find information about constants of simple equations of state for various types of working fluids. Open information about values of constants for the mBWR32 equation of state is available only for a few working fluids. However, there are various software packages, for example REFPROP [10] that allows us to calculate the array of fields of thermodynamic functions y i for a fairly wide range of working fluids. Based on them, it is possible to determine the constants for the working fluids of our interest. The gas constant R is determined as a ratio of the universal gas constant to the molecular weight of the considered working fluid. The remaining constants γ and G(i) are determined using the least squares method to assure the smallest square deviation of the dimensionless unknown function from the array base point values: where f i -the required thermodynamic function of the mBWR32 equation of state at point i of the array; y i -value of the thermodynamic function at point i, obtained using REFPROP software package; n-dimension of the base points array. Problem (2) can be solved in the following way: if γ is assumed as known and constant, then the condition (2) can be replaced by the condition: Expression (3) forms a system of 37 linear equations with respect to 37 unknowns G (j) for the thermodynamic functions: pressure, Helmholtz free energy, entropy and the partial derivative of pressure to density at constant temperature. The linear system of Equation (3) is solved by the Gauss method with dominant diagonal terms. The accuracy of calculations is set at a quadruple precision with 32 characters. Such a large mantissa is needed to maintain the required high accuracy. The global search for the solution of the problem defined by Equation (2) is carried out by varying γ in the range: where ρ * is the value of density at the critical point. Constants are found for the simultaneous fulfilment of condition (2) for the following thermodynamic functions: pressure, the Helmholtz free energy, entropy, and the partial derivative of pressure to density at constant temperature. The values of the constants of the mBWR32 equation, obtained using the method described above, for water vapour, HFE7100, and R227ea as working fluids in the investigated turbines are summarized in Appendix B. As initial data, thermodynamic arrays obtained using the REFPROP [10] program were taken at (1) 1100 points in the entire range of variation and 73 points on the saturation line for temperature from 280 to 800 K and pressure up to 96 MPa for water vapour, (2) 610 points in the entire range of variation and 35 points on the saturation line for temperature from 300 to 600 K and pressure up to 3 MPa for HFE-7100, and (3) 800 points in the entire range of variation and 45 points on the saturation line for temperature from 150 to 470 K and pressure up to 8.6 MPa for R227ea working fluid.

Method of Interpolation-Analytical Representation of Thermodynamic Functions
As mentioned earlier, the direct use of complex equations of state (IAPWS-95, mBWR32, etc.) in numerical algorithms for 3D viscous flow calculations leads to an increase in computational time by 1.5-2 orders of magnitude, which is unacceptable. To reduce computational cost, an interpolation-analytical method of representation of thermodynamic functions is used. This method was first applied to take into account the thermodynamic properties of water and steam in 3D calculations based on the IAPWS-95 equation of state in [29]. According to this approach, the required thermodynamic functions are determined from the following dependencies: where Z t = Z t (ρ,p), Z ρ = Z ρ (h,p), Z u = Z u (ρ,p), Z p = Z p (ρ,u), Z Cv = Z Cv (ρ,p), Z Cp = Z Cp (ρ,p), and Z S = Z S (ρ,p) are the dimensionless compressibility coefficients for the corresponding thermodynamic functions. These are determined by interpolation from pre-calculated arrays of the base points. This approach has significant advantages over the "direct" interpolation of thermodynamic functions, which is used in some software packages [30,31]. This is due to the fact that the range of variation of the dimensionless compressibility coefficients is much narrower, and unlike the required thermodynamic functions, the dimensionless compressibility coefficients are monotonic functions of the assumed dependent variables. Due to this, arrays of base points of much smaller dimensions can be stored to ensure the acceptable interpolation accuracy. Furthermore, to reduce the dimensions of the arrays without sacrificing the accuracy, the independent variables of pressure and density are considered in a logarithmic scale. To interpolate the compressibility coefficients, polynomials of the third degree are used. Values of the dimensionless compressibility coefficients are defined in the base points as: where the corresponding values of p, ρ, T, u, h, a, c p , c v and S are calculated using the thermal equations of state (IAPWS-95, mBWR32, etc.). A more detailed description of the interpolation method has been presented in the work of Rusanov et al. [29].

Low-Pressure Cylinder of 360 MW Steam Turbine
The steam flow in the low-pressure cylinder (LPC) of a 360 MW turbine was investigated experimentally and numerically. A view of the investigated LPC is shown in Figure 1. Thermodynamic parameters were measured in span-wise distribution in axial diffusers behind stages 3, 4, and 5 by probes 1, 2, and 3, respectively and in the exit diffuser by probe 4, as shown in Figure 1. Details of the experimental study and measurement method were described in earlier publications [4,32]. A numerical study of this LPC using various equations of state was first presented in [4]. In that paper, CFD results obtained using the equation of state for a perfect gas with constant and variable specific heats were compared and validated against experimental  A numerical study of this LPC using various equations of state was first presented in [4]. In that paper, CFD results obtained using the equation of state for a perfect gas with constant and variable specific heats were compared and validated against experimental data. It was shown that specifying the heat capacities as linear functions of temperature brought the obtained results closer to the experiment; however, this approach significantly worsened the stability of the computational code. The results of calculations of this LPC using the IAPWS-95 equation of state and their comparison with computational results based on the perfect gas equation and Van der Waals equation of state were shown in [29]. A fundamental improvement in the coincidence of experimental data and computational results while using the IAPWS-95 equation of state was illustrated.
In this study, the LPC is investigated numerically with the help of CFD code IPMFlow, making use of three equations of state-the Tammann equation, the IAPWS-95 equation, and the mBWR32 equation. To carry out the calculations, the following boundary conditions are assumed-the total pressure and temperature at the inlet, and the static pressure at the outlet (see Table 1). The calculations take into account two recuperative extractions after stages 3 and 4, with the mass extraction of 6 and 5.5 kg/s respectively, as well as inter-disk and over-the-shroud leakages [4]. The fifth stage rotor blade is unshrouded, so there is a meshed radial gap available for tip leakage flow.  The constants of the Tammann equation of state are determined from values of total thermodynamic parameters at the inlet and isentropic ones at the outlet ( Table 1). The values of constants of the Tammann equation and the isentropic enthalpy difference between the inlet and outlet of the flow part obtained on its basis are presented in Table 2. The isentropic enthalpy drop obtained from the Tammann equation of state is 9.2% greater than that obtained from the IAPWS-95 and mBWR32 equations of state.  Table 3 shows a comparison of the calculation results with experimental data, including pressures, temperatures, enthalpies, and flow rates at the outlet of the subsequent stages 3, 4, and 5. It can be seen that the differences between the experimental data and calculation results obtained using the IAPWS95 and mBWR32 equations of state are insignificant. The greatest difference is observed for the static temperature behind the third stage, −1.2%; in the other cases, the differences are much smaller than 1%, which is in favour of the presented computational method. The differences between the experimental data and calculation results obtained with the help of the Tammann equation of state are much more significant and, in some cases, reach 15%. However, it should be noted that the results obtained using the Tammann equation are more accurate compared to the perfect gas equation with constant heat capacities and similar in quality to the results obtained using the perfect gas equation with variable specific heats [4,29].

Radial ORC Turbine with HFE7100 Working Fluid
A 2.5 kW ORC micro power plant working on organic fluid HFE7100 supplied from the wood chips boiler was developed at the Institute of Fluid Flow Machinery, Polish Academy of Sciences (IMP PAN) as a model cogeneration unit dedicated for household applications. A general view of the micro-power plant is shown in Figure 4. A detailed description of the plant with its operational experience and the results of experimental studies can be found in the papers [33][34][35].

Radial ORC Turbine with HFE7100 Working Fluid
A 2.5 kW ORC micro power plant working on organic fluid HFE7100 supplied from the wood chips boiler was developed at the Institute of Fluid Flow Machinery, Polish Academy of Sciences (IMP PAN) as a model cogeneration unit dedicated for household applications. A general view of the micro-power plant is shown in Figure 4. A detailed description of the plant with its operational experience and the results of experimental studies can be found in the papers [33][34][35].  The main element of this cogeneration unit is a multi-stage radial microturbine. The turbine consists of four stages, where the first two stages are centripetal (radial inward), and the other two are centrifugal (radial outward). The turbine has small dimensions; the outlet edges of the fourth-stage rotor blades are located at a diameter of 74.5 mm (relative to the turbine rotation axis). The stator and rotor blades from all stages have tip gaps near the meridional contours of size 0.15 mm. The first stage is designed with partial admission supply (degree of partial admission − 0.5). The design rotational velocity is equal to 20,000 rpm. A view of the individual elements and meridional section of the ORC microturbine flow part is illustrated in Figure 5. The main element of this cogeneration unit is a multi-stage radial microturbine. The turbine consists of four stages, where the first two stages are centripetal (radial inward), and the other two are centrifugal (radial outward). The turbine has small dimensions; the outlet edges of the fourth-stage rotor blades are located at a diameter of 74.5 mm (relative to the turbine rotation axis). The stator and rotor blades from all stages have tip gaps near the meridional contours of size 0.15 mm. The first stage is designed with partial admission supply (degree of partial admission − 0.5). The design rotational velocity is equal to 20,000 rpm. A view of the individual elements and meridional section of the ORC microturbine flow part is illustrated in Figure 5.
Numerical studies of microturbine flow and comparison of computational results with the experimental data are carried out for three operational regimes. Boundary conditions for these regimes, including inlet total temperatures and pressures, as well as outlet static pressures, rotational speeds, and isentropic parameters are provided in Table 4. The calculations were made using two equations of state-the Tammann and mBWR32 equation. The Tammann equation constants are determined from total parameters at the inlet and isentropic parameters at the outlet as in Table 4. Values of the Tammann equation constants and the isentropic enthalpy difference between the inlet and outlet of the flow part obtained on its basis are presented in Table 5. It can be seen that the maximum difference between the isentropic enthalpy drop obtained from the Tammann equation of state is 8% compared to that found from the REFPROP [10] software. Numerical studies of microturbine flow and comparison of computational results with the experimental data are carried out for three operational regimes. Boundary conditions for these regimes, including inlet total temperatures and pressures, as well as outlet static pressures, rotational speeds, and isentropic parameters are provided in Table 4. The calculations were made using two equations of state-the Tammann and mBWR32 equation. The Tammann equation constants are determined from total parameters at the inlet and isentropic parameters at the outlet as in Table 4. Values of the Tammann equation constants and the isentropic enthalpy difference between the inlet and outlet of the flow part obtained on its basis are presented in Table 5. It can be seen that the maximum difference between the isentropic enthalpy drop obtained from the Tammann equation of state is 8% compared to that found from the REFPROP [10] software.    The comparisons of calculated results and experimental data are shown in Table 6. In the experiment, the power and, accordingly, the efficiency were measured at the generator terminals. The power and efficiency of the microturbine calculated in CFD studies were corrected taking into account the generator efficiency whose average value was determined at 85%. As it can be seen from the presented comparison, the results obtained using the mBWR32 equation of state are much closer to the experimental data than those obtained using the Tammann equation. For the mBWR32 equation, the maximum discrepancy between the experiment and computation in mass flow is 6.5%, in power, it is -2.7% and in efficiency, it is -4.1%. At the same time, for the Tammann equation of state, those discrepancies are 10.9%, 7.4%, and 6.2%, respectively.

Axial ORC Turbine with R227ea Working Fluid
An ORC plant working with R227ea was developed at the West Pomeranian University of Technology in Szczecin as a model installation for geothermal applications. This ORC unit is supplied by hot water from the district heating network. A general view of the ORC power plant is presented in Figure 6, whereas its more detailed description is given in the papers [35,36]. The ORC cycle is equipped with an axial single-stage microturbine with a degree of partial admission equal to 1/9. Turbine blades have low heights: 9.7 mm-stator blade and 11.4 mm-rotor blade. A view on the ORC turbogenerator and turbine flow path is shown in Figure 7. The small size of the turbine and a high degree of partial admission make its CFD calculation a non-trivial task. The ORC cycle is equipped with an axial single-stage microturbine with a degree of partial admission equal to 1/9. Turbine blades have low heights: 9.7 mm-stator blade and 11.4 mm-rotor blade. A view on the ORC turbogenerator and turbine flow path is shown in Figure 7. The small size of the turbine and a high degree of partial admission make its CFD calculation a non-trivial task. The ORC cycle is equipped with an axial single-stage microturbine with a degree of partial admission equal to 1/9. Turbine blades have low heights: 9.7 mm-stator blade and 11.4 mm-rotor blade. A view on the ORC turbogenerator and turbine flow path is shown in Figure 7. The small size of the turbine and a high degree of partial admission make its CFD calculation a non-trivial task. CFD investigations of the microturbine flow and comparison of computational results with the experimental data were carried out for four operational modes. Boundary conditions, including inlet total temperatures and pressures as well as outlet static pressures, rotational speeds, and isentropic parameters are provided in Table 7. Calculations CFD investigations of the microturbine flow and comparison of computational results with the experimental data were carried out for four operational modes. Boundary conditions, including inlet total temperatures and pressures as well as outlet static pressures, rotational speeds, and isentropic parameters are provided in Table 7. Calculations were conducted for two equations of state-the Tammann and mBWR32 equations. The Tammann equation constants determined from the total parameters at the inlet and isentropic parameters at the outlet, together with the isentropic enthalpy difference between the inlet and outlet of the flow part, are presented in Table 8. The maximum difference for the isentropic enthalpy drop obtained from the Tammann equation of state is 3% compared to that found from the REFPROP [10] software.  Table 9 shows the comparison of calculation results and experimental data. In the experiment, the power was measured at the generator terminals. Thus, the power and efficiency of the microturbine calculated in CFD studies were corrected, taking into account the generator efficiency whose average value was determined this time at 91%. Table 9. Comparison of experimental and calculated microturbine flow rate, power and efficiency for given operational regimes, as defined in Table 7. From the presented results, it can be seen that, in contrast to the previous examples, there is no explicit advantage of using the mBWR32 equation of state. Generally, there is a satisfactory agreement between the experimental data and calculation results obtained from both equations of state. According to the authors, this is due to a relatively small drop in thermodynamic parameters in the considered flow part, which leads to an insignificant difference in isentropic enthalpy drops obtained during the calculations using the mBWR32 and Tammann equations of state (less than 3%). It can be concluded that for small isentropic enthalpy drops, the Tammann equation of state satisfactorily describes the thermodynamic properties of the working fluid.

Conclusions
The method of interpolation-analytical reconstruction of the complex IAPWS-95 and mBWR32 equations of state was presented to take into account the real properties of working fluids in 3D CFD calculations. The method for determination of the constants of the modified Benedict-Webb-Rubin equation of state with 32 terms was proposed. The constants of the mBWR32 equation were calculated for the HFE-7100 and R227ea working fluids.
The results of the 3D flow calculation in the LPC of a steam turbine and ORC microturbines with the HFE7100 and R227ea working fluids were presented. Calculations were performed using the Tammann equation and the proposed interpolation-analytical method of representation of the IAPWS-95 and mBWR32 equations of state. The comparison of the calculation results with experimental data showed satisfactory agreement in the case of the IAPWS-95 and mBWR32 equations for a wide range of change of thermodynamic parameters. It was found that the Tammann equation of state provides acceptable results for relatively small drops in thermodynamic parameters in the flow part, when the isentropic enthalpy drop obtained during the calculations using the Tammann equation differs from the exact value by not more than 3-4%.
The method is able to improve the quality of flow modelling in thermal turbomachinery, especially in terms of the flow of fluids, the behaviour of which differs significantly from that of an ideal gas. This functionality was not available in earlier versions of the in-house IPMFlow software.
The function f i is a component of the Helmholtz free energy containing non-exponential terms, whereas f v contains exponential terms. I 1 -I 6 are density integrals of the exponential terms. An additional polynomial consisting of five terms is introduced into Equation (A2), which leads to an increase in the number of constants G from 32 to 37. The introduction of additional components increases the accuracy in determination of the Helmholtz free energy and other thermodynamic functions.
In view of the form of Equations (2) and (A2), the required thermodynamic functions take the form given below: Internal energy: Isobaric heat capacity:  Entropy: Sonic speed: Table A1 shows values of the constants of the mBWR32 equation, which were obtained using the method described in Sections 2.2 and 2.3 and the equations in Appendix A. Table A1. The values of the constants of the mBWR32 equation.

Header
Water