A one-dimensional analytical calculation method for obtaining normal shock losses in supersonic real gas flows

The calculation of isentropic flow and normal shock waves of real gases are important, especially in the preliminary design of turbo-machinery and test rigs. In an ideal gas, the relations for one-dimensional isentropic flow and normal shock waves are well known and can be found in standard textbooks. However, for fluids exhibiting strong deviations from the ideal gas assumption universal relations do not exist due to complex equations of state. This paper presents a analytical method for the prediction of isentropic real gas flows and normal shock waves, based on the Redlich-Kwong (RK) equation of state. Explicit expressions based on a series expansion for describing isentropic flow of Novec™ 649 are compared to Refprop data and ideal gas equations. For moderate pressures the RK method is in very good agreement with the Refprop data, while the ideal gas equations fail to predict the real gas behaviour. The same observations are made for normal shock calculations, where both real gas methods yield very close results. Especially the predicted stagnation pressure losses across a shock wave are in excellent agreement.

649 are compared to Refprop data and ideal gas equations. For moderate pressures the RK method is in very good agreement with the Refprop data, while the ideal gas equations fail to predict the real gas behaviour. The same observations are made for normal shock calculations, where both real gas methods yield very close results. Especially the predicted stagnation pressure losses across a shock wave are in excellent agreement.

Nomenclature
Latin Symbols

Introduction
Despite the rapid advancements in computational fluid dynamics (CFD), a design and optimization approach solely based on CFD simulations is still time consuming. Therefore onedimensional methods, especially in the preliminary design phase, are important. In an ideal gas, the relations for one-dimensional isentropic flow and normal shock waves are well known and can be found in standard textbooks. However, for fluids exhibiting strong deviations from the ideal gas assumption, universal relations do not exist. It is the purpose of this paper to develop a method for the description of one-dimensional isentropic flow and normal shock losses in supersonic real gas flows, using the analytical Redlich-Kwong (RK) equation of state (EoS). Fluids exhibiting real gas behaviour are of special interest for organic Rankine Cycle (ORC) applications. An ORC process is a thermodynamic cycle, based on the classical steam Rankine cycle that utilizes organic fluids instead of water as working fluid. The key parameter in terms of organic fluid properties is a low boiling point, making the usage of low-temperature heat sources possible. In recent years, research has focused on so called dry fluids, with an overhanging saturation curve, which is favourable to avoid blade erosion [1].
The reliable design of ORC components requires both numerical and experimental work. For this reason, the Laboratory for Thermal and Power Engineering at Muenster University of Applied Sciences, Germany is building up a closed loop wind tunnel for organic vapors, aiming at the investigation of real-gas effects in nozzles and axial turbine blade cascades [2]. However, to determine the working point of the test rig, the overall losses at operating conditions have to be estimated. The dominating contribution to the total losses occurs in the test section. Details on its design can be found in [3,4]. To get a better quantitative estimate about the pressure losses and hence the required pressure rise to drive the wind tunnel, a normal shock wave in the test section was considered. The computational method to calculate these losses for an organic vapour are developed and tested in this contribution. Another important technical application for normal shock waves is the correction of total pressure measurements from supersonic Pitot-tubes.
The paper is structured as follows. After a brief literature review on existing real gas flow models, the present approach is described in detail. First, the thermodynamic real gas model is presented, followed by methods to describe isentropic flow and normal shock calculations. The results are then presented and discussed.

Literature Review
A number of models for describing the steady one-dimensional isentropic flow of a real gas can be found in the literature. An excellent historical review of isentropic flow models has been published by Sullivan [5]. Reimer [6] and Johnson [7] presented isentropic relations based on modified ideal gas relations that accurately describe fluid behaviour. However, in both cases only slight deviations from ideal gas behaviour occurred, which can be attributed to the fluids and the range of operating conditions selected for the study. Kouremenos and Kakatsios [8] chose a similar approach and defined explicit relations approximating real gas behaviour, with a similar mathematical form as the ideal gas equations. Later, they extended this work to the description of thermodynamic properties using the RK EoS [9]. One of the few experimental works has been carried out by Bier et al. [10]. They investigated the flow of the refrigerant R22 in an axially symmetric de Laval nozzle for Mach numbers up to M = 2.2 . The experimental results were compared to one-dimensional computations using an advanced EoS as well as an approximation procedure, that could reproduce all measured values with deviations within a few percent.
A different approach on isentropic flow makes use of the fundamental derivative Γ as defined by Thompson [11]. The fundamental derivative of gas dynamics is a quantitative measure of the variation of the speed of sound with respect to density in an isentropic change of state [12]. Thompson and Sullivan [13] reported a method for predictions of the sonic conditions in a real gas utilizing the fundamental derivative. The results are relatively simple relations based on a Taylor series expansion. However, the authors remark, that while the application of the fundamental derivative leads to a simple and effective description of isentropic real gas flows, the calculation of Γ itself is not necessarily simple. Colonna et al. [12] investigated the possibility of computing the fundamental derivative using different EoS. They found non-analytical or reference EoS to be superior to analytical formulations. Especially in the vicinity of the critical point, analytical EoS fail to predict the correct physical behaviour.
One of the few analytical approaches to calculate normal shock waves of real gases was published by Kouremenos [14]. He derived a non-linear equation using the RK EoS, but this first study lacks information on the results and accuracy of the method. Later on, Kouremenos et al. [15] presented a modified calculation procedure and investigated the influence of different EoS and stagnation conditions on normal shock waves in steam. They conclude, that the RK EoS provides accurate solutions for moderate pressures, but becomes increasingly inaccurate for higher pressures and temperatures.

Thermodynamic Real Gas Model
The description of the thermodynamic behaviour of a real gas requires a suitable EoS. The most accurate type are non-analytic or fundamental EoS [16], explicit in Helmholtz energy. Equations in this form are capable of describing the thermodynamic behaviour of pure substances with an accuracy that probably exceeds measurement uncertainties [16]. Unfortunately, this type of EoS has to be solved numerically and cannot be used to derive simple working equations.
For the present work the widely used Redlich-Kwong EoS [17] has been selected. The equation in the form p(T, v) is represented by Herein, p, v, T are pressure, specific volume, and temperature, and R is the specific gas constant R = R univ /M m . The constants α and β are functions of pressure and temperature at the critical point (p c , T c ) and are given by If the compressibility factor Z is close to unity, the ideal gas EoS is reproduced by Eq 1 [18]. A full description of the thermodynamic behaviour using a cubic EoS requires additional analytical expressions for enthalpy and entropy. However, it can be shown, that a fluid model consisting of the three equations p(T, v), h(T, v), and s(T, v) is thermodynamically equivalent to a fundamental equation of state [18]. For a thermal EoS with the independent variables temperature and volume, the enthalpy of a real gas is defined as [18] where u(T, v) describes the specific internal energy, given by Herein, u ig is the specific internal energy of the corresponding ideal gas state (p ig → 0). The first integral in Eq (4) describes the temperature dependency of u for an ideal gas state, whereas the second integral accounts for deviations from ideal behaviour. The specific entropy of a real gas can be dealt with in a similar manner, resulting in the expression where s ig is the specific entropy of the corresponding ideal gas state (p ig → 0). The selection of the reference point is arbitrary. In the present work, internal energy u ig and entropy s ig were set to zero at T ig = 273.15 K and p ig = 1 Pa. For the speed of sound a holds the exact relation [18] Equation 6 shows that the speed of sound can be derived from the thermal EoS p(T, v) and the real gas isochoric specific heat capacity c v (T, v): Herein, the ideal gas specific heat capacity c ig v (T ) is approximated by a second-order polynomial The constants A, B, and C can be taken from the literature or from a fit procedure. For this paper, Refprop data was used. Equation 8 was fitted by means of a least-squares method and reproduces the exact values within 0.15 %. After integration Eq (3) and Eq (5) become: .
(10) Thus, Eq (1,9,10) represent a closed thermodynamic real gas model, which has been used for the following investigations.

Isentropic Flow Relations
From a practical point of view, simple working formula for the description of isentropic changes in a real gas are desirable. A method for obtaining relations in the form p/p 0 , T /T 0 , v/v 0 , and A * /A in terms of Mach number is outlined in the following.
The governing equations of steady one-dimensional isentropic flow are the continuity equation and the energy equation where c is the velocity, v the specific volume, and h and s the specific enthalpy and entropy. In the following, an isentropic expansion process proceeding from a known set of stagnation conditions p 0 , T 0 , v 0 along an isentrope s 0 = const. to an arbitrary point p 1 , T 1 , v 1 is considered. Furthermore, the stagnation or reservoir conditions are assumed to be known. Consequently, the isentrope s 0 passes through this point. To describe this isentrope we calculate the fluid properties along this path for an arbitrary number of points n i . The local Mach number M i at each point is defined by Herein, the local speed of sound a i can be obtained from Eq 6 and the local velocity c i is related to the enthalpy by Start Compute: Generate linearly spaced vector:  Pressure, temperature, specific volume and area ratios are then fitted to a polynomial in terms of Mach number by means of a least-squares fit. The present work has shown, that simple ideal-gas-like formulations, as used for instance by Kouremenos and Kakatsios [8], fail to predict real gas behaviour of fluids like Novec TM 649. Therefore the equations were extended by means of a series expansion. A sixth-order polynomial was found to be sufficient to describe real gas behaviour in the present case. The polynomials F (M) are of the generic form: Herein, ζ i describe the fitted constants. The computational procedure is summarised in Fig 1. The correct value of the lower pressure limit for the fit p lb has to be determined iteratively, based on a prescribed maximum Mach number M max . It is worth noting, that this algorithm is independent of the EoS and can be used with cubic as well as fundamental EoS.

Normal Shock Calculations
A shock wave is a discontinuity in a very small region that occurs between two boundaries in supersonic flow. Normal shock waves are perpendicular to the flow direction and transform a supersonic flow to a subsonic flow regime. A schematic sketch of a normal shock wave is shown in Fig 2 (left). Figure 2 (right) shows the idealised process in the h-s-diagram. Shock waves are dominated by viscous and heat transfer effects. Therefore the assumption of isentropic flow is no longer valid. As a result from the second law of thermodynamics the entropy across a shock wave has to increase, i.e. s 02 > s 01 .
The description of a normal shock wave involves five variables, i.e. static pressure, temperature, specific volume, enthalpy, and velocity. Therefore in addition to the thermal EoS p = p(T, v), enthalpy equation h = h(T, v), continuity of mass Eq (11), and conservation of energy Eq (12), a fifth equation is necessary to compute the changes across a normal shock. This is achieved by introducing the conservation of momentum: Herein, p describes the static pressure, c is the velocity, and v the specific volume.
For an ideal gas, this set of equations can be solved in closed form, resulting in the normal shock or Rankine-Hugoniot equations [19]. For non-ideal fluids this is usually impossible due to the complexity of the EoS. A way of finding an analytical solution is by means of cubic EoS.  Figure 2. Schematic sketch of an idealised normal shock wave (left) and changes across a normal shock wave in the specific enthalpy versus entropy plane (right).
Using Eq (1) and Eq (10) for thermal and enthalpy EoS, in principle allows for a solution of one analytical expression containing only one unknown, i.e. v 2 = f (T 2 ). However, the present study has shown that the resulting expression becomes quite cumbersome and not very reliable to solve. Far better results were obtained from a system of two non-linear equations with two unknowns T 2 and v 2 that can be solved numerically using a root finding algorithm, such as the Newton-Raphson method [20]. Supplying the corresponding ideal gas values as initial data, results in a quick convergent behaviour. Using the conservation of momentum in combination with the RK EoS (Eq (1)) to eliminate the unknown pressure p 2 results in: Similarly the energy equation, Eq (12), in combination with the enthalpy equation, Eq (10), can be used to obtain the second expression: The numerical solution of Eq (17) and Eq (18) provides T 2 and v 2 . All other properties downstream from the shock wave can be obtained from the RK EoS (Eq (1)) and enthalpy and entropy equations (Eq (10,9)).

Results and Discussion
The calculation procedures presented in this paper are tested with the fluorinated ketone Novec TM 649, which is the working fluid for the closed loop wind tunnel [2]. Figure 3 shows an h-s-diagram for Novec TM 649 [21]. Of special interest for the present study are stagnation conditions in the superheated region with pressure and temperature in the range of p 0 = 4...6 · 10 5 Pa and T 0 = 370...420 K. This range was chosen because it represents the possible operating conditions of the planned closed loop wind tunnel. The fluid possesses the desired overhanging saturation curve, as mentioned in Sec 1. Deviations from ideal gas behaviour can already be observed for moderate pressures, i.e. at 5 · 10 5 Pa and 390 K, the compressibility factor is already as low as Z = 0.829.
In the first step, accuracy and consistency of the RK real gas model were compared with the reference EoS (McLinden et al. [22]), which is implemented in Refprop [21]. For stagnation conditions of 5 · 10 5 Pa and 390 K an exemplary comparison between both EoS is presented in Fig 4. The RK predictions show deviations mostly lower than 10 %. The biggest errors occur in the prediction of the specific volume, with a maximum deviation of 20.9 %. Errors increase for higher pressures, near the saturation line, and in the vicinity of the critical point. However, for the range of interest in this work, the RK EoS is a reasonable approximation of the real fluid behaviour. It should also be noted, that some initial tests with other cubic EoS, such as the formulations by Aungier [23], Soave [24], and Peng and Robinson [25] showed no significant improvement over the RK EoS.
In the next step, the isentropic relations were fitted to the fluid properties according to the procedure described in Sec 4 and Fig 1. Isentropic expansions using the RK and Refprop EoS were compared to each other over a wide range of stagnation conditions, up to the critical point region. Table 1 shows the maximum error for the RK EoS with respect to the corresponding Refprop data. The results show, that the prediction of the specific volume is afflicted with the highest deviations up to 40 % near the critical point, followed by the static pressure with deviations up to 18 %. The lowest errors were observed for the temperature with maximum deviations around 2 %. Furthermore, the results indicate, that the RK EoS becomes less reliable for higher pressures and temperatures, especially for expansions from stagnation conditions above p 0 = 12 · 10 5 Pa, errors higher than 10 % are being observed.
Deviations were also found to increase in the vicinity of the saturation curve. As shown in Tab 1, for a stagnation pressure of 5 · 10 5 Pa three different stagnation temperatures were considered. For T 0 = 378 K, the point closest to the saturation curve, the deviations are larger than for the two points with higher temperatures, that are farther away from the saturation curve. A similar behaviour, albeit with higher errors, has been observed for higher pressures. For the range of interest in this work (p 0 = 4...6·10 5 Pa, T 0 = 370...420 K), the RK EoS has been found to deliver re-      Finally, the method for calculating the normal shock was tested, using the following procedure: For a given set of stagnation conditions p 0 , T 0 the fluid properties for a specified Mach number M 1 are computed, using the isentropic relations described in Sec 4. The next step is the normal shock calculation. In case of the RK EoS the procedure described in Sec 5 is followed, whereas for the Refprop case three governing equations (Eq (11,12,16)) in combination with the EoS by McLinden et al. [22] were solved numerically, using a Newton-Raphson method. Subsequently, the stagnation conditions behind the shock are determined from the EoS and the results are checked for consistency, i.e. the stagnation enthalpy has to remain constant h 02 = h 01 , whereas the entropy has to increase s 02 > s 01 . Figure 6 shows the variation of M 2 , p 1 /p 2 , p 02 /p 01 , T 1 /T 2 and v 2 /v 1 in terms of upstream Mach number M 1 . The results for the RK EoS were obtained from the procedure described above, while the solutions with Refprop were computed by directly solving the three governing equations (Eq (11,12,16)) in combination with the EoS by McLinden et al. [22]. A stable convergent behaviour was achieved by setting the corresponding ideal gas solution as initial values for the iterative solution. Again, excellent agreement was found between RK and Refprop predictions, while the ideal gas formulations fail to predict the non-ideal fluid behaviour. An additional parameter in the description of normal shock waves is the strength of a shock wave P , which is defined as the ratio of the increase in static pressure to the initial static pressure upstream from the shock wave [19]: The results shown in Fig 7 (right) indicate, that normal shock waves calculated from real gas formulations are weaker than the solutions obtained from the ideal gas equations. However, the absolute loss in stagnation pressure for the real gas model is up to 40 % higher than the corresponding losses for an ideal gas. Of major interest is also the consistency of the calculations.  Figure 6. Variation of M 2 , p 1 /p 2 (left) and T 1 /T 2 , v 2 /v 1 (right) in terms of M 1 for normal shock waves in Novec TM 649. Stagnation conditions were set at p 0 = 5 bar and T 0 = 390 K.  Tests with different EoS, i.e. using Refprop for the description of isentropic flow and the RK EoS for normal shock calculations, resulted in inconsistencies and non-physical solutions.

Conclusion and Outlook
A method for calculating one-dimensional isentropic flow and the normal shock of a real gas is presented. A thermodynamic real gas model based on the cubic RK EoS accurately predicts the fluid properties of Novec TM 649 for moderate pressures, but becomes less accurate for higher pressures and near the critical point. Simple working formulae based on a sixth-order polynomials are sufficient to describe the isentropic flow of real gases. The agreement with Refprop data is excellent, while the ideal gas formulation fails to predict real gas behaviour. A system of two non-linear equations, using the RK EoS, is solved to calculate normal shock waves. The results are compared to Refprop and ideal gas data. While both real gas models are in very good agreement, the ideal gas shows strong deviations. Compared to the ideal gas model, both real gas models predict higher absolute pressure losses, but lower values for the shock strength.
This paper shows, that a analytical real gas model based on the RK EoS accurately predicts isentropic flow and normal shocks for moderate pressures. Therefore it can be used as an alternative to a purely empirical approach using Refprop. Although, the prediction of the specific volume by the RK model is afflicted with relatively high errors, the normal shock losses are in very good agreement with Refprop data. The model has been found to be reliable for moderate pressures and temperatures, higher errors occur near the critical point region. Also the behaviour at supercritical conditions has to be investigated.
The results of the present paper suggest the following areas for future research: So far only Novec TM 649 has been considered. An extension to other fluids is necessary to further investigate the reliability of the method. While initial tests with other cubic EoS have not shown significant improvement over the RK EoS, the investigation of other cubic EoS (e.g. Soave [24]) could further improve the method. Subsequently, an extension to oblique shock waves will be of interest.