A Practical Real Gas Model in CFD

A PRACTICAL REAL-GAS MODEL IN CFD

Teemu Turunen-Saaresti, Jin Tang and Jaakko Larjola Lappeenranta University of Technology

Skinnarilankatu 34, 53850 Lappeenranta, Finland e-mail: teemu.turunen-saaresti@lut.fi web page: http://www2.et.lut.fi/lvt/en/

Key words: Real gas, CFD, least square fitting, nozzle

Abstract. Polynomial functions were developed to efficiently compute real-gas properties in a Navier-Stokes flow solver. Two-variable relative least square fitting is applied to generate these functions. The computational time requested is small and the precision of the polynomial function can be easily adjusted. In this article real-gas model is implemented for superheated toluene and the flow through a supersonic stator nozzle of an Organic Rankine Cycle turbine is simulated. The obtained results are compared to the measured ones.

1 INTRODUCTION

CFD (Computational Fluid Dynamics) is now widely used in the engineering and other scientific fields. Because of the use of the various mediums and processes occurring in the non-ideal gas region in the engineering applications, real-gas models describing the gas properties are needed. A sufficiently straightforward implementation of gas property models for the different mediums is necessary in order to maintain a sufficiently short lead-time of CFD calculations.

The gas properties can be calculated in CFD in several ways: using a look-up-table approach, an equation of the state or polynomial functions. Complex equations of the state can be used when the fluid does not behave as an ideal gas. In this article polynomial fittings are chosen. This approach is utilized because polynomial functions are easy to generate and implement into the code. Also the need for extra variable tables is avoided.

In the presented real-gas model, the temperature and the pressure are independent variables. Density, internal energy, dynamic viscosity and thermal conductivity are functions of the temperature and the pressure. The functions are described by polynomial equations. For a given property table, the coefficients of the polynomial equations are generated. Matlab® has been used to calculate the fitting coefficients. The polynomial equations are generated based on pvT data taken from reference1.

equations but this on the other hand increases the computational time.

The simulation of an ORC (Organic Rankine Cycle) supersonic nozzle is considered as an example in this article. The Navier-Stokes flow solver Finflo has been applied for the simulation. Finflo has been used to calculate flow fields in several turbomachinery applications2,3,4 The simulation results are compared with the experiment results.

2 NUMERICAL PROCEDURE

Reynolds-averaged Navier-Stokes equations were solved by a Finite-Volume method. Navier-Stokes equations can be written in Cartesian coordinates as follow

( )

( )

( )

where is the vector of the conservative variables and Q is the

source term. Here ρ is the density, ρu, ρv and ρw are the momentum components, E is the total internal energy, k is the kinetic energy of the turbulence, and ε is the dissipation of the kinetic energy of the turbulence. There are additional quantities such as the velocity

, pressure p, viscous stresses τ

(

)

Vr = r+ r+ r ij and heat fluxes qi in the flux vectors F(U),

G(U) and H(U). The total internal energy is defined as k V V e E =ρ +ρ ⋅ +ρ 2 r r ₍₂₎

where e is the specific internal energy. In a finite-volume technique the flow equations have the following integral form

( )

_∫

∫

∫

After performing the integrations, for a computational cell i the flow equations are

∑

where the sum is taken over the faces of the computational cell. The flux for the face is defined as H n G n F n Fˆ = _x + _y + _z (5)

The MUSCL- type approach is adopted for both primary flow variables and conservative turbulent variables and inviscid fluxes are evaluated using Roe’s flux difference splitting5. Details of the implementation can be found in reference6 and only essential features are presented here related to the real gas model implementation. The flux is calculated as

( )

where T is a rotation matrix used to transform the dependent variables to a local coordinate system normal to the cell surface. Now, only the Cartesian form F of the flux is needed and it is calculated from

(

)

[

( ) ( )

]

_∑

where Ul and Ur are the solution vectors on the different sides of the face, r(k) is a Jacobian matrix 1 − Λ = ∂ ∂ = R R U F A (8)

The corresponding eigenvalue is λ(k), and α(k) is the corresponding characteristic variable obtained from R-1ΔU, where ΔU = Ur – Ul. The flux components can be calculated easily but the calculation of the flux difference requires additional thermodynamic quantities. Different quantities can be utilized because the selection of the eigenvectors and characteristic variables is not unambiguous. This flow solver utilizes the derivates (∂ρ/∂h)p and (∂ρ/∂p)h.

Viscous fluxes are evaluated using thin-layer approximation. Turbulence is modeled using Chien’s k-ε turbulence model7. Convergence of the calculations is accelerated by a multi-grid method8. The y+ of the grid is maintained to be less than 1 in order to solve boundary layer properly. The flow solver Finflo is developed at the Helsinki University of Technology9. 2.1 Real gas model

One criterium for the real-gas model adopted in this study was the short generation time of the model for different real gases. Therefore the polynomial function approach was utilized. The temperature and the pressure are independent variables in this flow solver. Density, internal energy, dynamic viscosity and the thermal conductivity are required as a function of temperature and pressure:

(

)

(

)

(

)

(

)

The used solver also requires the derivates of the density described equations (10) and (11)

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ p p p T h T T c T p p ρ ρ ρ ρ ρ ρ 1 1 (11)

where the derivates (∂ρ/∂T)p and (∂ρ/∂p)T can be derived from equation (9) and heat capacity

can be calculated from equation

p p p T p T e c ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ = ρ ρ2 (12)

Polynomial fittings for the density, the internal energy, the dynamic viscosity and the thermal conductivity are constructed in the following functional form:

(

)

_∑

For example the fourth order polynomial function is shown in equation (14)

4 2 14 3 2 1 13 2 2 2 1 12 2 3 1 11 4 1 10 3 2 9 2 2 1 8 2 2 1 1 3 1 6 2 2 5 2 1 4 2 1 3 2 2 1 1 0 1 x a x x a x x a x x a x a x a x x a x x a x a x a x x a x a x a x a a y + + + + + + + + + + + + + + = (14)

The order of the equation (13) can be adjusted to change the accuracy of the fittings. Fourth-order functions are found to be accurate enough for the simulations presented here. Multiple regression was used to solve unknown coefficients (a0 … a19). Equation (15) shows the

criteria used in evaluating the fitting quality. It takes only few hours to generate the necessary functions with Matlab.

(

)

3.1 Validation against table data

3.2 Description of the test case

The test case used in this article is a supersonic ORC turbine nozzle where some experimental data is available from measurements using a pressure/temperature probe. The working fluid is toluene. The ORC turbine is a radial type. Figure 2(a) shows the computational grid used in the calculation without the probe. It represents one of the 20 stator vanes around the rotor. The flow enters the nozzle at the upper right and exits at the lower right towards the rotor. Mass flow and temperature are imposed at the inlet. The cyclic boundary condition is used in front of and behind the vane. Moreover, the pressure is imposed at the outlet. Two different cases are calculated. First the stator nozzle without the measuring probe and second the same nozzle geometry with the measuring probe. The measuring probe is placed between the stator trailing edge and the rotor leading edge. Calculation with the probe is conducted to obtain information about the pressure and temperature around the probe. Figure 2(b) shows the surface grid with the probe and figure 2(c) shows the detailed surface grid around the measuring probe. The measured values available are the temperature and the pressure at the inlet of the nozzle and the temperature measured with the probe and the static pressure measured at the wall around the probe. Table 1 shows the main geometrical and working fluid data. The grid sizes used in the simulations are shown in table 2.

(a) (b) (c)

Figure 2. (a) Surface grid of one nozzle. Every second grid line is visible. (b) Surface grid of one nozzle with measuring probe. (c) Detailed surface grid around the probe located at the

The real-gas model described is rather efficient. The real-gas code consumes a CPU time per cell 3.4% less than the corresponding ideal gas code. The reason for shorter CPU time might be the Sutherland’s formula which is used to calculate viscosity in the ideal gas code. The Sutherland’s formula is computationally expensive because it contains a division operation. Also the real-gas model uses the same amount of memory as the ideal gas code. However, the initial conditions were noticed to be a critical aspect to conduct successful runs. Trial and error method was used to choose the initial values.

Table 1. Main geometrical and working fluid data.

Rotor outer diameter R1 224.2 mm

Ratio of nozzle outlet and rotor inlet radius R1’/R1 1.05

Degree of reaction r 0.26

Height of nozzle b1’ 7.3 mm

Number of nozzles 20

Breadth of nozzle throat 0.90 mm

Mach number at nozzle exit 2.40

Mach number at rotor inlet 0.97

Mass flow qm 1.24 kg/s

Rotor speed (nominal) N 478 1/s

Table 2. Computational grid sizes. grid size

2D 13312 2D with probe 23040

3.3 Comparison of calculated and measured data

Table 3 shows the measured and calculated values. The calculated data is taken at the inlet, before the supersonic jets are mixing and at the outlet. Calculated values are shown in cases with and without the probe. The values taken from the computational results have been area averaged. The calculated values from the probe have been averaged around the probe surface.

Table 3. Measured and calculated temperature and pressure values. measured CFD with probe CFD without probe p [bar] T [K] p [bar] T [K] p [bar] T [K]

inlet 29.44 587.6 30.77 587.7 29.87 587.6

before mixing 0.59 473.0 0.57 473.4

probe 0.59 540 1.95 542.5

outlet 0.53 473.0 0.54 475.3

(a) (b)

Figure 3. (a) Computed temperature around the probe and (b) computed pressure around the probe.

4 CONCLUSIONS

The real gas model based on polynomial fits is implemented in a Reynolds averaged Navier-Stokes solver. The real gas model is generated using two-variable relative least square polynomial fittings which are made based on a given property table and a certain range of the working fluid. Matlab® is used to calculate fitting equations. Polynomial fitting are used because they can be generated fast and they are easy to implement to the code. The real gas model was noticed to consume 3.4% less CPU time per cell than corresponding ideal gas code.

REFERENCES

[1] Goodwin, R., D., 1989. Toluene Thermophysical Properties from 178 to 800 K at

Pressures to 1000 Bar, Journal of Physical and Chemical Reference Data, Vol. 18, No.4. [2] Turunen-Saaresti, T., Larjola, J.,2004. Unsteady Pressure Field in a Vaneless Diffuser of

a Centrifugal Compressor: An Experimental and Computational Analysis, Journal of Thermal Science, Vol. 13, No. 4.

[3] Hoffren, J., Talonpoika, T., Larjola, J., Siikonen, T., 2002. Numerical Simulation of Real-Gas Flow in a Supersonic Turbine Nozzle Ring, Journal of Engineering for Real-Gas Turbines and Power, vol. 124, p. 395-403.

[4] Turunen-Saaresti, T., Reunanen, A., Larjola, J., 2005. Effect of Pinch on the Performance of a Vaneless Diffuser in a Centrifugal Compressor, 6th European Conference on

Turbomachinery Fluid Dynamics and Thermodynamics, 7-11. March 2005, Lille, France. [5] Roe, P., L., 1981, Approximate Riemann Solvers, Parameter Vectors and Difference

Schemes, Journal of Computational Physics, n 43, p 357-372

[6] Siikonen, T., 1995. An Application of Roe’s Flux-difference Splitting for k-ε Turbulence Model, International Journal of Numerical Methods in Fluids, vol. 21, n. 11.

[7] Chien, K., 1982, Predictions of channel and boundary-layer flows with a low-Reynolds- number turbulence model, AIAA Journal, n. 20, p. 33-38.

[8] Jameson, A., Yoon, S., 1986. Multigrid solution of Euler Equations Using Implicit Schemes, AIAA Journal, Vol. 24, No. 11.