Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Real Gas Models in Coupled Algorithms Numerical Recipes and Thermophysical Relations †

Real Gas Models in Coupled Algorithms Numerical Recipes and Thermophysical Relations † International Journal of Turbomachinery Propulsion and Power Article Real Gas Models in Coupled Algorithms Numerical Recipes and Thermophysical Relations 1, ‡,§ 1,§ 1 2 Lucian Hanimann * , Luca Mangani , Ernesto Casartelli , Damian M. Vogt and Marwan Darwish Competence Center for Fluid Mechanics and Hydro Machines, Lucerne University of Science and Arts, 6048 Horw, Switzerland; luca.mangani@hslu.ch (L.M.); ernesto.casartelli@hslu.ch (E.C.) Institut für Thermische Strömungsmaschinen und Maschinenlaboratorium, Universität Stuttgart, 70569 Stuttgart, Germany; damian.vogt@itsm.uni-stuttgart.de Departement of Mechanical Engineering, American University of Beirut, Beirut 110236, Lebanon; darwish@aub.edu.lb * Correspondence: lucian.hanimann@hslu.ch; Tel.: +41-41-349-34-58 † This paper is an extended version of our paper published in Proceedings of the European Turbomachinery Conference ETC13, Lausanne, Switzerland, 8–12 April 2019; Paper No. 218. ‡ Current address: Technik & Architektur, University of Applied Sciences Lucerne, Technikumstrasse 21, 6048 Horw, Switzerland. § These authors contributed equally to this work. Received: 13 May 2020; Accepted: 29 July 2020; Published: 3 August 2020 Abstract: In the majority of compressible flow CFD simulations, the standard ideal gas state equation is accurate enough. However, there is a range of applications where the deviations from the ideal gas behaviour is significant enough that performance predictions are no longer valid and more accurate models are needed. While a considerable amount of the literature has been written about the application of real gas state equations in CFD simulations, there is much less information on the numerical issues involved in the actual implementation of such models. The aim of this article is to present a robust implementation of real gas flow physics in an in-house, coupled, pressure-based solver, and highlight the main difference that arises as compared to standard ideal gas model. The consistency of the developed iterative procedures is demonstrated by first comparing against results obtained with a framework using perfect gas simplifications. The generality of the developed framework is tested by using the parameters from two different real gas state equations, namely the IAPWS-97 and the cubic state equations state equations. The highly polynomial IAPWS-97 formulation for water is applied to a transonic nozzle case where steam is expanded at transonic conditions until phase transition occurs. The cubic state equations are applied to a two stage radial compressor setup. Results are compared in terms of accuracy with a commercial code and measurement data. Results are also compared against simulations using the ideal gas model, highlighting the limitations of the later model. Finally, the effects of the real gas formulations on computational time are compared with results obtained using the ideal gas model. Keywords: real gas; numerics; compressor 1. Introduction CFD continues to be a critical tool in the design and development of turbomachinery equipment. It is used in the evaluation of the performance of complex multistage turbines at design and off-design points, and in the optimization of single blades. Furthermore, the continuous increase in available computational resource has allowed for a substantial improvement in the quality of simulations as compared to a decade ago. While the majority of compressible flows still uses the perfect gas Int. J. Turbomach. Propuls. Power 2020, 5, 20; doi:10.3390/ijtpp5030020 www.mdpi.com/journal/ijtpp Int. J. Turbomach. Propuls. Power 2020, 5, 20 2 of 32 model, it is becoming more common to work on applications where real gas effects are no longer negligible. Example where these situations occur include low pressure stages of steam turbines [1] or ORC turbines [2]. Many authors have investigated the impact of various real gas models, mainly focusing on possible achievements in the accuracy of the results and the impact on computational time due to the increased complexity of the thermophysical models, see, e.g., [2–7]. The later point, the increase in computational time, has, however, received less attention. While the complexity of real gas relations is probably the most obvious cause, other factors also play a role, such as the iterative nature of the solution procedures used to update thermophysical properties. This non-linearity is often addressed by using lookup tables, see, e.g., [5]. In this case the user needs to know in advance the range of the thermophysical states, which is not always available at the start of a simulation. This is further complicated by the fact that, prior to convergence, some intermediate solutions can run out of the tabulated range. For such situations, a stabilization criterion is needed. Additional issues arise concerning the resolution of the generated lookup table. A low resolution may lead to a faster numerical procedure, but can diminish the quality of the results. The high sensitivity of the results to lookup table resolution was described in [8]. Other techniques to address the weakness of lookup tables include the use of structured grids for the lookup tables and the use of various interpolation methods whether bi-linear or bi-cubic ([5]). In this paper the real gas relations are incorporated algebraically, i.e., without using lookup tables, into an in-house, pressure based, all-speed, fully coupled solver [9–11]. The main intention is therefore to contribute to the understanding of the necessary modifications considering a compressible all-Mach solver when using real gas state equations. Starting from the basic real gas relation, the conservation equations and boundary conditions and the derivation of the developed numerical techniques are presented in details with a main focus on retaining the performance and robustness of the coupled solver. Furthermore, the computational cost of the use of real gas relations is thoroughly investigated. Although the developed methods to deal with real gas relations are general, two specific relations are evaluated in more detail. It is for the reason that these models are later used for the validation and it is considered a contribution to help the reader to apply the procedure to his desired real gas model. The highlighted models are the cubic state equations, e.g., [12] and the IAPWS database for water at any state [13]. It will be shown that the increase in computational time often does not justify the problems arising from lookup table formulations. Replacement with tabulated data is, however, straight forward. Three test cases are used to validate and evaluate the presented numerical approach. In the first case, the consistency of the approach is demonstrated. This is accomplished by showing that the solution of additional iterative procedures needed for real gas formulations recovers the results of the perfect gas model when applied to a 2D transonic bump case. In the second test case, the IAPWS-97 formulation [13] is applied to a nozzle configuration with the operating range varying from supersonic to transonic flow regimes. The expansion for the transonic conditions yield to operating conditions in the two-phase region. Results are compared to those obtained from a commercial code. The last test case is an industrial type two stage compressor at one operational speed for the full range from choking flow to the stability limit, the cubic state equations are used as they are the manufacturer ’s choice for this application. Results are then compared against measurement data. Furthermore, for cases two and three, an evaluation of the computational cost of the real gas implementations is presented. 2. Thermodynamic Modeling of Real Gas Real gases are compositions of materials at thermophysical states under which its molecules are dense enough so that intermolecular forces can no longer be neglected. This is in contrast to the theory of ideal gases. For ideal gases, the molecules are sparsely distributed in space and these forces are not taken into account. It is this assumption that leads to the well known ideal gas state equation. Int. J. Turbomach. Propuls. Power 2020, 5, 20 3 of 32 The fundamental derivative of gas dynamics is a key parameter in understanding real-gas flows. It is a measurement for the change in speed of sound through a change in density, its definition is provided in Equation (1). r ¶c G = 1 + (1) c ¶r g + 1 While this expression reduces to a constant value of G = for perfect gas and is always greater than one, this does not apply for real gas. A special example are Bethe–Zel’dovich–Thompson (BZT) fluids, which are heavy polyatomic compositions. For these materials, G can even have negative values when approaching the critical point. The effect is an inverted gas dynamic behavior, where, contrary to perfect gas, the Mach number decreases as the flow expands. 2.1. General Comment on Real Gas Modeling The modeling of real gases in a numerical framework is accompanied with a variety of additional challenges. This is a direct result of the increased complexity of the state equations compared to the simple ideal gas state equations. The non-linearities introduced require iterative updates of a variety of quantities, which will be a major part of what follows. It also requires a profound analysis of all underlying equations and thermophysical relationships. Entropy, internal energy and enthalpy are a function of two variables for the general case of a real gas. Since cubic equations are used in the final validation (in which p = f (T, v)), the equations and partial derivatives are given using T, v as primary input variable. However, this does not limit the generality of the approach and just requires a careful derivation of what is exemplified below. p, T, v are interchangeable through the state equation, but depending on its formulation it may simplify the analysis if the input variables are changed to, e.g., T, p. 2.2. Real Gas State Equations Various models have been proposed to better account for intermolecular forces when encountering denser fluids, each having its own justification depending on the operting conditions and the working fluid. Analytical representation for superheated and two-phase steam is given by the model of Young [14] in virial form or for water at any state in the form given through IAPWS [13]. When modeling hydrocarbons, cryogenic fluids or refrigerants, well known virial equations are those of Benedict–Webb–Rubin (BWR) [15] and its modifications. Another group of real gas state equations is represented in the form of cubic state equations. The original formulation was introduced by Van der Waals [16] but many other models have emerged to overcome some of its drawbacks. Modifications have been applied to improve the phase transition, see Soave et al. [17], or the prediction of liquid phase density and accuracy in the critical region by Peng et al. [12]. 2.2.1. Cubic State Equations The cubic state equation, used for the final validation case, is based on the formulation provided by Peng et al. [12] and therefore given in Equation (2). T and p are temperature and pressure at the cr cr critical point and w being the accentric factor. R T aa p = (2) 2 2 v b v + 2bv b 2 2 T 0.45724R T s cr a = a(T) = 1 + k 1 a = (3) T p cr cr R T s cr b = 0.0778 k = 0.37646 + 1.54226w 0.26993w (4) cr Int. J. Turbomach. Propuls. Power 2020, 5, 20 4 of 32 2.2.2. IAPWS-97 The IAPWS-97 [13] formulations provide accurate state equations for the properties of water at any state. The equations are highly polynomial functions and provided for five different subregions described in Figure 1. Region 1 is the liquid phase and region 3 accounts for the transition from liquid to gaseous phase above the critical point. Region 2 and 5 are gaseous states and region 4 is the two-phase region below the critical point. Figure 1. IAPWS regions. While this set of equations allows for a highly accurate representation of the properties of water, they can lead to a massive increase in computational costs. Even by the use of provided backwards functions to eliminate at least partially the iterative procedures the increase in computational time is still considerable. The usual approach is therefore the use of lookup tables generated using IAPWS formulations. 2.3. Caloric State Equations and Definition of Entropy The numerical algorithm requires values of internal energy, entropy and enthalpy at a certain point. This can be the case during the update of variables that are not a direct outcome of a conservation equation such as temperature or density. However, it may also be required when updating boundary conditions, such as given stagnation conditions at the inlet. Using real gas state equations, the definitions of the caloric state equations for enthalpy and internal energy are functions of two intensive properties, e.g., temperature and pressure, and are no longer only a function of temperature alone. The derivation of generalized state equations is therefore briefly described here on the basis of internal energy. The definition for the change in internal energy is given in Appendix A and the final expression is shown in Equation (5). The integral form of this differential equation can be found using an integration along paths as shown in Figure 2. ¶ p de = c dT + T p dv (5) ¶T The value of e(T, p) is computed through integration along lines of constant temperature ( A, C) and constant specific volume (B) using a given reference state e(T , p ). The final expression for internal r r energy is thus given in Equation (6), formulations for enthalpy and entropy (see Appendix E) can be found using the same basic idea. Int. J. Turbomach. Propuls. Power 2020, 5, 20 5 of 32 Z Z Z v T v ¶ p ¶ p e(T, p) = T p dv + [c ] dT + T p dv +e(T , p ) (6) v r r ¶T ¶T v T v r r ¥ v v T T | {z } | {z } | {z } It is obvious that such a representation of the integral value can be understood as the sum of the ideal gas contribution plus a correction to adjust for real gas effects. As stated earlier, ideal gas behavior is encountered when intermolecular forces do not affect the properties of the fluid. This is always valid if the molecules are separated by an infinite distance. In other words the density goes towards zero or specific volume to infinity. Hence the integration along line B accounts for the ideal gas contribution, and integration along the lines A, C corrects the value to correspond to the real gas change in the state variable. The ideal gas heat capacity in term B can be described with a variety of models. The simplest is the assumption of calorically perfect gas (i.e., c = const.), other options are Sutherland models or NASA’s polynomial approach [18]. As stated in the introduction, the final validation of the coupled solution strategy for real gas CFD is based on an industrial sized two-stage radial compressor. It is the manufacturer ’s choice to use the Peng–Robinson [12] state equation. Appendix C therefore gives a detailed derivation of Equation (6) for this type of state equation. To provide an example of the generality of the approach the derivations for internal energy for the Benedict–Webb–Rubin (BWR) [15] state equation are given in addition in Appendix D. T, p A ¥ T , p r r Figure 2. Integration path. 2.4. Stagnation Conditions Evaluation of stagnation conditions given static values or vice versa is a necessity frequently encountered in engineering practice. The relation is usually based on the assumption of an isentropic state change, leading to the well known expressions for incompressible and compressible (assuming ideal gas) flows as given in Equations (7) and (8) exemplarilly for pressure. p = p + rjjujj (7) g1 g 1 p = p 1 + Ma (8) These relationships usually provide the basis in the derivation of stagnation conditions, however they do not apply to real gas. In the derivation of Equation (8) it is assumed that the isentropic exponent is a constant, a simplification which is not valid for real gas. Some authors ([19,20]) argue that the variations are, however, small and Equation (8) can be applied with minimal loss of accuracy. While this may be a justified assumption for a variety of applications, its generality is questionable. The generalized approach used in the present article is therefore explained in later sections when describing boundary conditions and the update of domain internal quantities. Int. J. Turbomach. Propuls. Power 2020, 5, 20 6 of 32 2.5. Specific Heat Capacity As already introduced when discussing the stagnation conditions, the heat capacity and heat capacity ratio do no longer follow the well known ideal gas relations. The heat capacity ratio can vary for different thermophysical states and requires a generalized approach. The derivation used in the implementation is therefore given in detail in Appendix F. The final expression is given in Equation (9) ¶ p ¶T c = c T (9) p v ¶ p ¶v 2.6. Mach Number The Mach number is a frequently used quantity for the post-processing and analysis of numerical results. In the following we will therefore provide a general formulation for the speed of sound. The definition of the Mach number is given by Equation (10). jjujj Ma = (10) The speed of sound for any compressible medium is given by Equation (11). ¶ p c = (11) ¶r Detailed derivation of the partial differential of pressure with respect to density under isentropic conditions can be found in Appendix G. The final expression is given in Equation (12). ¶ p ¶ p = v (12) ¶r c ¶v The partial differentials can again be found in Appendix B. 3. Numerical Environment Segregated pressure-based solution algorithms [21,22] based on the SIMPLE family of algorithms have been applied successfully over the past decades. The popularity of these algorithms is partly due to their low memory requirement based on the sequential method in solving the RANS equations. Another reason is the relatively simple implementation. On the downside, these algorithms suffer from convergence breakdown for large-scale problems. The reason for this issue is that the naturally strong pressure-velocity coupling is numerically decoupled when solving the momentum and continuity equation sequentially. Owing to the developments in computer architecture, coupled solution strategies for pressure-based algorithms have regained interest in recent years [23,24]. These algorithms solve the momentum and continuity equation together, i.e., coupled. The numerical framework used in this article is based on the coupled solution strategy presented in [9–11]. The coupled system of linearized momentum and continuity equations for each cell of the computational domain takes the form given in Equation (13). Int. J. Turbomach. Propuls. Power 2020, 5, 20 7 of 32 2 32 3 2 32 3 2 3 u p u p uu uv uw uu uv uw u a a a a u a a a a u b p p p p p nb c nb nb nb nb v p 6 76 7 6 v p76 7 6 7 vu vv vw vu vv vw v a a a a v a a a a v b 6 76 7 6 76 7 6 7 p p p p p nb nb nb nb nb 6 76 7 + 6 76 7 = 6 7 (13) w p w p wu wv ww wu wv ww w 4a a a a 54w 5 4a a a a 54w 5 4b 5 p p p p nb c nb nb nb nb pu pv pw p p pu pv pw p p p p a a a a p b a a a a p nb c p p p p nb nb nb nb | {z } | {z } p nb The final system of linearized equations is solved using algebraic multigrid algorithms on the basis of an additive correction approach [25,26] with a block-ILU smoother. Figure 3 highlights the main body of the coupled solution strategy. The blocks of the flow chart will also define the structure of the present article. In Section 5, the conservation equations are revisited with a focus on the difference when integrating real gas thermophysics. This is followed by a detailed analysis of boundary conditions mainly focusing on the implementation of total conditions, see Section 6. The theoretical part is finalized showing the chosen strategies for the update of flow properties in Section 7. n n1 Update Boundary Conditions t = t + Dt Coupled System Sequential Systems Momentum-Equation Energy Equation A p 3 A p Assemble Terms Turbulence Equations Compute D uu uv uw u p a a a a Assemble Total Conditions T, v = f (h, p) vu vv vw v p a a a a Property update wu wv ww w p Continuity Equation a a a a Correct U and p pu pw pw p p a a a a Assemble Terms Figure 3. Flow chart of the algorithm. 4. Governing Equations The above described numerical framework solves the RANS equations as: ¶r ¶ru Continuity Equation + = 0 (14) ¶t ¶x ¶ru u ¶t ¶ru i j ¶ p i j Momentum Equation + + f = 0 (15) ¶t ¶x ¶x ¶x j i j ¶rh ¶ p ¶ Energy Equation + ru h + q u t = 0 (16) i 0 i j i j ¶t ¶t ¶x t in the momentum and energy equation is the deviatoric stress tensor given as: i j ¶u ¶u j 2 ¶u i l t = m + m d (17) i j i,j ¶x ¶x 3 ¶x j i l q in the energy equation accounts for diffusive heat flux modeled through Fourier ’s approach as: ¶T q = l (18) ¶x 5. Numerical Discretization of Governing Equations The governing equations for compressible fluid dynamic analysis are the momentum, continuity, energy and turbulence equations. As described in Equation (13) the momentum and continuity equations are solved simultaneously to retain the pressure-velocity coupling, followed by a sequential solution of the energy and turbulence equations. A complete derivation of the governing equations, Int. J. Turbomach. Propuls. Power 2020, 5, 20 8 of 32 schemes and used solvers would go beyond the scope of this article. Implementation details of the numerical framework is best described in [9–11,27–32]. The main focus will be the integration of a generalized approach for real gas physics. Partial recapitulation of the derivations concerning conservation equations and boundary conditions will be given where necessary. 5.1. Momentum Equation The momentum equation remains unchanged compared to the derivation for perfect gas. It is, however, important in order to introduce the coupled solver framework and therefore given in Equation (19). ¶ru V F + (r u u ) + p t S = 0 (19) bi å f f j f i f f i j f i ¶t F accounts for body forces and S is the normal vector of the face with the length of its bi f i area. The implicit contributions of the single terms of the momentum equation to the coupled, block-structured coefficient matrix are assembled in what is marked as a red box in Figure 3. The coupling of the momentum equation to the pressure based continuity equation is dominated by the pressure gradient and assembled in a , a and a , shown in red in Equation (13). This is one of the u p v p w p major advantages of coupled solution strategies. The pressure gradient is treated implicitly and thus updated during inner iterations. Segregated solution strategies consider the pressure gradient in the momentum equation as a constant contribution to the source when solving the momentum equations. 5.2. Continuity Equation For pressure based numerical frameworks, the continuity equation is used to obtain an equation for the pressure which serves as a constraint to fulfill the conservation of mass. Equation (20) shows the continuity equation in its semi discrete form with U = u S being the velocity flux across a face. f f i f i dr V + r U = 0 (20) å f f dt The temporal derivative is expanded using the general assumption that the density can be expressed as a function of pressure and enthalpy. This is a procedure that has been described in [7]. The linearization of this term is therefore given as shown in Equation (21). Expressions have to be found for the thermophysical derivatives based on the chosen real gas state equation. This derivation will be exemplified in a later section. dr ¶r dh ¶r d p = + (21) dt ¶h dt ¶ p dt p h Equation (21) differs from perfect gas computations where it is usually assumed that the change of density with respect to time can be simply expressed as given in Equation (22). dr R T = (22) dt dt The second term of the continuity equation, the convective contribution, is discretized using Newton linearization, see Equation (23). Variables indexed with  are previous iteration values while variables without this index represent current iteration values. r U = r U + r U r U (23) f f f f f f f f Int. J. Turbomach. Propuls. Power 2020, 5, 20 9 of 32 To avoid checkerboarding pressure fields on collocated arrangements and to introduce the necessary implicit dependency on pressure together with a coupling of the continuity equation to the velocity, the velocity flux of the new time step (U ) is computed using the Rhie–Chow [33] interpolation technique. 0 1 B ! ! C B C ¶ p ¶ p B C U = U D S (24) B C f f f ,i,j f i B ¶x ¶x C j j @ f fA | {z } | {z } A B U is the velocity flux on the face, computed using adjacent cell values through linear interpolation, contributing to the coefficients a , a and a in the coefficient matrix A . This couples the continuity pu pv pw p equation to the velocity field, shown in green in Figure 3. The pressure gradient A is evaluated directly on the face , while term B is evaluated based on cell center values and then interpolated on the face. The term D is the inverted 3  3 coefficient sub-matrix A (given in green in Figure 3) after f p3 cell-to-face interpolation. The density r in Equation (23) is computed using a first order linearization, to obtain the necessary implicit dependency on pressure. This is in contrast to the derivation for ideal gas, where density is substituted with pressure by simply using the state equation. ¶r ¶r r = r + (h h ) + ( p p ) (25) f f f f f f ¶h ¶ p p, f h, f Inserting Equations (24) and (25) in Equation (23) and using Equation (21) in Equation (20) leads to the final continuity equation in pressure form shown in Equation (26). ¶r dh ¶r d p V + + r U f f ¶h dt ¶ p dt | {z } p h m ˙ 0 0 1 1 ! ! ¶ p ¶ p @ ¯ @ A A + r U D S å f f f ,i,j f i ¶x ¶x j j f f | {z } m ˙ ¶r ¶r + U (h h ) + ( p p ) r = 0 (26) å f f f f f f ¶h ¶ p p f h f | {z } m ˙ The term r U introduces an implicit coupling of the continuity equation to the velocity and pu pv pw p p will therefore contribute to a , a and a . Three terms are implicit contributions to a , namely c c c c ¶ p ¶r ¶r d p r D S , U p and V . å å f p f f i f f f f ¶x ¶ p ¶ p dt h f h To summarize, the difference in the derivation of the pressure equation for real gas compared to idealized gas is listed below: dr • Linearization of transient term dt • Linearization of the density of the actual time step r ¶r ¶r • Partial derivatives and ¶h ¶ p p h Int. J. Turbomach. Propuls. Power 2020, 5, 20 10 of 32 5.2.1. Thermophysical Partial Derivatives in the Continuity Equation The derivation of the partial derivatives given above is crucial for the stability of the solver. In this ¶r section the derivation of is given as an example. Expansion of the expression is first done ¶ p through triple product rule. ¶r ¶h ¶r = (27) ¶ p ¶ p ¶h h v p | {z } | {z } B1 B2 Part B1 can extended as follows ¶h ¶s = T + v dh = Tds + vd p (28) ¶ p ¶ p v v ¶v ¶v ¶s = T + v MR: = (29) ¶T ¶T ¶ p s s v ¶s ¶v = T + v Triple product rule (30) ¶T ¶s v T c b ¶s c ¶v b v v = T + v = and = (31) T a ¶T T ¶s a v T c b a = + (1 Ta)v c = c vT (32) v p a b In this derivation a is the thermal expansion coefficient and b is the compressibility as given in Equations (33) and (34). ¶ p 1 ¶T a = (33) v ¶ p ¶v 1 ¶v b = (34) v ¶ p Part B2 can be extended as follows ¶r ¶r ¶T = Product rule (35) ¶h ¶T ¶h p p p ¶r 1 = (36) ¶T c 1 ¶v 1 = (37) ¶T c a 1 = (38) v c The final expression can therefore be written as given in Equation (39) ¶r ¶h ¶r Ta a = = rb (39) ¶ p ¶ p ¶h c h v p Int. J. Turbomach. Propuls. Power 2020, 5, 20 11 of 32 5.3. Energy Equation The definition of the conservation equation for total internal energy is given in Equation (40). ¶re ¶ + ru e + u p + q u t dV = 0 (40) j 0 j j i i,j ¶t ¶x For compressible flows it is, however, more convenient to express the energy equation in terms of total enthalpy as h = e + p/r. Equation (40) is therefore rewritten as Equation (41). 0 0 ¶rh ¶ p ¶ + ru h + q u t dV = 0 (41) i 0 i j i j ¶t ¶t ¶x Compared to ideal gas state equations, the elementary changes are in the modeling of the heat flux. Using Fourier ’s law to relate the heat flux to the local temperature gradient with l being the conductivity of the working fluid, Equation (42). is derived. Z I ¶q ¶T dV = q n d A  q S = l S (42) i i å i f f i å f f i ¶x ¶x W ¶W i i f f This term would lead to a completely explicit contribution to the source of the energy equation, diminishing the stability of the iterative solution procedure. What is needed is an algebraic relation for the dependency of temperature in terms of enthalpy to obtain an implicit formulation of the heat conduction. For calorically perfect ideal gases, this poses no challenge, since enthalpy can be expressed as h h = c (T T ). A generalized approach, as considered in this publication, allows no such 2 1 p 2 1 simplification. To recover at least partly the implicit contribution a first order linearization of the enthalpy is introduced as given in Equation (43). ¶h ¶h h = h + (T T) +  ( p p) (43) ¶T ¶ p p T | {z } Expression A can be found in Appendix B and is simply the definition of the specific heat capacity c . B is zero because the pressure is kept constant during the solution of the energy equation. An expression for the temperature of the actual time step can therefore be found as given in Equation (45) h h T = T + (44) = T + (45) Substituting the temperature T in Equation (42) using Equation (45) allows to rewrite the Laplacian contribution as given in Equation (46). Term a introduces the implicit contribution providing the needed stability of the discretization. ¶T ¶T f ¶h l S = l S + S å f f i å f f i å f i ¶x ¶x c ¶x i i p f i f f f f f f | {z } ¶c h p + l S (46) å f f i ¶x p f | {z } b Int. J. Turbomach. Propuls. Power 2020, 5, 20 12 of 32 At convergence term a and b vanish due to their dependency on the correction of enthalpy, hence the heat conduction is modeled accurately. Part b can usually be neglected since the gradient of the heat capacity is often rather small. 6. Boundary Conditions Considering the implementation of boundary conditions for a coupled framework, it is crucial to identify any coupling between velocity and pressure as well as implicit dependencies on the variable that is solved for. A detailed derivation of the implemented boundary conditions can be found in [32,34]. While these articles give a detailed insight in the numerical discretization of a variety of boundary conditions, they assume perfect gas state equations. The focus of this section will therefore remain on the difference when integrating real gas state equations. 6.1. Total Conditions The use of total conditions at the inlet is a very common boundary condition. A detailed derivation is therefore given in what follows. As already shown, in a coupled framework, the momentum and continuity equation are assembled into one system of equations. In order to properly account for the boundary contribution, each term of the conservation equations has to be analyzed with respect to the constraint given by the user. This is described in detail in [32] for ideal gas. The derivation of total conditions in this section is given by first providing the algorithm for consistent total to static conversion. The algorithm uses the common assumption of an isentropic deceleration to stagnation conditions but without all the simplifications valid only for perfect gas. This is followed by the discretization of the Navier–Stokes equations itself. With the exception of the pressure gradient term, this is in close agreement to what is given in [32]. The main body of this section will therefore be a detailed revision of the discretization of the pressure gradient term when dealing with real gas state equations. 6.1.1. Isentropic Total-to-Static Conversion As already mentioned above, the computation of static values using given total conditions can no longer be obtained using well-known algebraic relations. It is, however, possible to find an iterative solution procedure to accurately compute these quantities. The algorithm used is a Newton–Raphson root finding method. The roots are the two underlying constraints of the total to static conversion. Root 1 as given in Equation (47) is the isentropic constraint, while root 2 as given in Equation (48) is simply using the definition of total enthalpy. f (x , x ) = s(T , v ) s(T, v) = 0 (47) 1 1 2 0 0 f (x , x ) = (h(T , v ) 0.5jjujj ) h(T, v) = 0 (48) 2 1 2 0 0 In the following the Newton–Raphson root finding method is explained as a reference. This procedure is later applied to a variety of solution strategies used for finding iterative solutions for thermophysical changes in real gas. The procedure starts with a first order Taylor expansion of the roots. n n n n ¶ f (x , x ) ¶ f (x , x ) 1 1 n+1 n+1 n n 1 2 n+1 n 1 2 n+1 n f (x , x ) = f (x , x ) + (x x ) + (x x ) (49) 1 1 1 2 1 2 1 1 2 2 ¶x ¶x 1 2 n n n n ¶ f (x , x ) ¶ f (x , x ) 2 2 n+1 n+1 n n 1 2 n+1 n 1 2 n+1 n f (x , x ) = f (x , x ) + (x x ) + (x x ) (50) 2 2 2 2 1 2 1 1 1 2 ¶x ¶x 1 2 From Equations (47) and (48) it is known that these equations are equal to zero. Hence we can reformulate them as shown in Equations (51) and (52). Int. J. Turbomach. Propuls. Power 2020, 5, 20 13 of 32 n n n n ¶ f (x , x ) ¶ f (x , x ) 1 1 n n 1 2 1 2 0 = f (x , x ) + Dx + Dx (51) 1 1 2 1 2 ¶x ¶x 1 2 n n n n ¶ f (x , x ) ¶ f (x , x ) 2 2 n n 1 2 1 2 0 = f (x , x ) + Dx + Dx (52) 2 1 2 1 2 ¶x ¶x 1 2 The corrections Dx and Dx are defined as: 1 2 n+1 n Dx = x x (53) n+1 n Dx = x x (54) 2 2 This is system of equation is then rewritten in matrix form as: 2 3 n n n n ¶ f (x , x ) ¶ f (x , x ) 1 1 " # " # 1 2 1 2 n n 6 7 Dx f (x , x ) 1 1 ¶x ¶x 1 2 1 2 6 7 = (55) n n n n n n 4 5 ¶ f (x , x ) ¶ f (x , x ) 2 2 Dx f (x , x ) 1 2 1 2 2 2 1 2 | {z } | {z } ¶x ¶x 1 2 | {z } x By inverting matrix A the system can be solved as shown in Equation (56) 2 3 n n n n ¶ f (x , x ) ¶ f (x , x ) " # " # 1 2 1 2 1 2 n n 6 7 Dx 1 f (x , x ) 1 ¶x ¶x 1 1 2 6 1 1 7 = (56) n n n n n n n n n n n n 4 5 n n ¶ f (x ,x ) ¶ f (x ,x ) ¶ f (x ,x ) ¶ f (x ,x ) ¶ f (x , x ) ¶ f (x , x ) Dx 1 2 2 1 f (x , x ) 1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 2 1 2 ¶x ¶x ¶x ¶x 1 2 1 2 | {z } | {z } ¶x ¶x 1 1 The variables x , x are then updated with the correction and the procedure repeated until certain 1 2 convergence criterion is achieved. Applied to the roots given in Equations (47) and (48) with x = T and x = p, the following derivatives have to be found using the chosen real gas state equation. n n n n ¶ f (x , x ) ¶h ¶ f (x , x ) ¶h 1 1 1 2 1 2 = and = (57) ¶x ¶T ¶x ¶v 1 2 v T n n n n ¶ f (x , x ) ¶s ¶ f (x , x ) ¶s 2 2 1 2 1 2 = and = (58) ¶x ¶T ¶x ¶v 1 2 v T ¶h ¶s ¶s can be found in Appendix B. is equal to and can be replaced with ¶T ¶T ¶v v v T ¶ p ¶h by using Maxwell relations. is found as given in Equation (59). ¶T ¶v v T ¶h ¶e ¶ p = + v + p h = e + pv (59) ¶v ¶v ¶v v T T ¶e In a second step, is derived from Equation (6) as given in Equation (60). ¶v ¶e ¶ p = T p (60) ¶v ¶T T v Hence Equation (59) can be written as given in Equation (61). The derivatives of the pressure are easily found through the chosen state equation. ¶h ¶ p ¶ p = T + v (61) ¶v ¶T ¶v v v T Int. J. Turbomach. Propuls. Power 2020, 5, 20 14 of 32 6.1.2. Pressure Gradient Term As already pointed out, the discretization of the Navier–Stokes equation for given total conditions is best described following what is written in [32]. When integrating real gas state equations, the only difference is in the discretization of the pressure gradient term of the momentum equation. Using Green–Gauss theorem to rewrite the volume integral value of the pressure gradient as a surface integral, this term is written as given in Equation (62). An expression for the value of the pressure at the boundary face is therefore needed. The requirement is a strong coupling to face adjacent cell values of pressure and velocity, leading to an implicit contribution to the coefficient matrix. This ensures a robust integration of boundary conditions into the linearized equation system. Z I ¶ p dV = pd A  p S (62) å f f i ¶x W ¶W Using the well known strong coupling between pressure and velocity, a first order linearization of the pressure with respect to the velocity flux U is introduced as given in Equation (63). b indexes values evaluated at boundary faces. Again, the velocity flux is defined as U = u S with S being a i f i f i face normal vector with the length of the area of the face. As for the continuity equation this leads to an implicit contribution to the cell velocity and an implicit coupling to the pressure field of the continuity equation. ¶ p p = p + (U U ) (63) b b b b ¶U ¶ p What is needed is an expression for the gradient and the value of the velocity flux of the new ¶U time step U . The value for the velocity flux is computed using Rhie–Chow interpolation, given in Equation (64). ! ! ! ¶ p ¶ p U = U D S (64) b b,i,j ib ¶x ¶x j j b c The gradient of the pressure at the boundary is approximated as given in Equation (65). ¶ p p p b b c = (65) ¶x jdj With p being the cell pressure value andjdj the normal distance from face to cell center. Inserting Equation (63) in Equation (64) leads to the final definition for the boundary pressure value given in Equation (66). ¶ p p = p + qU + qD p + qF U (66) c c b b ¶U DjS j f i D = (67) jdj ¶ p q = 1 + D (68) ¶U ¶ p ¶ p ˜ ˜ F = D p + D U + D S (69) f i b b ¶U ¶x b c ¶ p The only unknown is therefore the gradient of pressure with respect to the velocity flux . ¶U This is usually computed using the algebraic relations for the static to stagnation conversion, given in Equation (70) for incompressible cases and Equation (71) for compressible cases. Int. J. Turbomach. Propuls. Power 2020, 5, 20 15 of 32 p = p + 0.5rjjujj (70) g1 g 1 p = p 1 + M (71) As explained before, direct algebraic relations that describe the change in pressure of a fluid particle when brought to rest isentropically are not explicitly available for the general case of a real gas. ¶ p To approximate the gradient the definition given in Equation (72) is used. ¶U ¶ p ¶ p ¶u p + D p p ¶u = = lim (72) ¶U ¶u ¶U D p!0 u( p = p + D p) u( p) ¶U The static value of the pressure is computed using an iterative scheme given the total conditions based on T and p as described above. In a second step, the velocity is computed that corresponds to 0 0 a slightly increased pressure level p = p + D p, following an isentropic change. Hence the roots are: f (x , x ) = s(T , p ) s(T , v ) = 0 (73) 2 0 0 2 2 1 1 f (x , x ) = p p(T , v ) = 0 (74) 2 1 2 2 2 2 In addition to the derivatives given in the section dealing with total to static conversion the partial derivatives of the pressure are needed. ¶ f (x , x ) ¶ p ¶ f (x , x ) ¶ p 2 1 2 2 1 2 = = (75) ¶x ¶T ¶x ¶v 1 2 v T Using the values for T and v the value of the static enthalpy h = f (T , v ) is computed according 2 2 2 2 to Equation (76). The calculation of internal energy follows the previously explained integration along paths. h(T , v ) = e(T , v ) + p v (76) 2 2 2 2 2 2 The velocity magnitude corresponding to these updated values is further found through the definition of total enthalpy as given in Equation (77). jju( p )jj = 2.0(h(T , p ) h(T , v )) (77) 2 0 0 2 2 The value for the derivative of pressure with respect to velocity is then found by substituting these results into Equation (72). 7. Update of Flow Properties As shown in Figure 3, the outcome of the conservation equations for momentum, continuity and energy are updated values for velocity, pressure and total enthalpy. There are, however, many other flow properties that are not directly obtained but which have to be computed based on these known values. This may be because they are intrinsically necessary for the solution procedure, as, e.g., density or viscosity, or because they are valuable quantities for post-processing such as Mach-number or entropy. This section explains their computation when operating in regions, where real gas effects can no longer be neglected. 7.1. Update of Temperature and Density Given the velocity, the pressure and total enthalpy from the conservation equations, it is possible to update the temperature and density independent of the state equation. First, the static Int. J. Turbomach. Propuls. Power 2020, 5, 20 16 of 32 enthalpy is computed using the given total enthalpy h and the velocity based on the solution of the RANS equations. h = h 0.5u u (78) i i The value for temperature and density are computed using the Newton–Raphson method. The roots of the method are: f (x , x ) = h h(T, v) = 0 (79) 1 1 2 f (x , x ) = p p(T, v) = 0 (80) 2 1 2 The rest of the Newton–Raphson root finding procedure can be done in analogy to Section 6.1.1. 7.2. Total Pressure and Total Temperature In contrast to the problem of specified total conditions for boundary conditions as described in Section 6.1, finding stagnation values for static conditions is usually required for the remaining cells and boundary faces of the computational domain. As described, no explicit relation exists for the static to stagnation conversion, asking for an iterative procedure for the isentropic state change. Consequently, Newton–Raphson root finding method is used. Given the static values of temperature and specific volume, the two roots therefore are given in Equations (81) and (82). With x and x being 1 2 the requested stagnation values T and v . 0 0 f (x , x ) = s(T, v) s(T , v ) = 0 (81) 1 1 2 0 0 f (x , x ) = (h(T, v) + 0.5u u ) h(T , v ) = 0 (82) 2 1 2 i i 0 0 As written above, having the roots of the Newton–Raphson method defined, the solution is in close analogy to what is given in Section 6.1.1. With the stagnation values for T and v the total 0 0 pressure is easily reconstructed using the chosen state equation. 8. Validation and Results This section provides an assessment of the generality and accuracy of the implemented numerical solution procedure when dealing with real gas state equations. Initially, the well-known transonic bump of Peric [35] is used. This test case is used to test the consistency of the implemented procedures when perfect gas state equations are used. The expressions and iterative algorithms given above have to return the same solutions as when using the mentioned algebraic relations only valid for perfect gas. Validity of the implemented procedures is given through comparison of the RMS-residuals. The second test case is a transonic nozzle configuration operated using steam based on the IAPWS-97 formulation [13]. Results are compared to a commercial code. To demonstrate the generality of the approach to any type of state equations, cubic state equations are used for a two-stage radial compressor. Results are compared against commercial code and measurement data. 8.1. Consistency of the algorithms for Perfect Gas State Equation The transonic bump test case of Peric [35] is given in Figure 4. The thickness to chord ratio of the bump is 10% and the computational domain consists of 53,000 quadrilateral grid elements. The inlet boundary conditions is set to a total pressure of 1.4 bar and a total temperature of 783K, turbulence intensity is 0.1% and the eddy viscosity ratio is set 2. The outlet static pressure is 1bar. The working fluid is air with the turbulence modeled using the k-w-SST model of [36]. Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 18 of 32 289 8.1. Consistency of the algorithms for perfect gas state equation 290 The transonic bump test case of Peric [35] is given in Fig. 4. The thickness to chord ratio of the 291 bump is 10% and the computational domain consists of 53’000 quadrilateral grid elements. The inlet 292 boundary conditions is set to a total pressure of 1.4bar and a total temperature of 783K, turbulence 293 intensity is 0.1% and the eddy viscosity ratio is set 2. The outlet static pressure is 1bar. The working Int. J. Turbomach. Propuls. Power 2020, 5, 20 17 of 32 fluid is air with the turbulence modeled using the k− ω-SST model of [36]. Figure 4. Transonic Bump Figure 4. Transonic bump. 295 To prove the consistency of the implemented generalized approach, the results should reduce to To prove the consistency of the implemented generalized approach, the results should reduce 296 the same output as when directly using perfect gas assumptions, when deriving the fundamental to the same output as when directly using perfect gas assumptions, when deriving the fundamental 297 equations. Fig. 5a shows the comparison between the direct solution and the iterative procedure for equations. Figure 5a shows the comparison between the direct solution and the iterative procedure for 298 the ideal gas case. The small difference results from the fact that the generalized approach uses an the ideal gas case. The small difference results from the fact that the generalized approach uses an 299 iterative procedure to evaluate isentropic conditions and to update temperature and density, hence the iterative procedure to evaluate isentropic conditions and to update temperature and density, hence the 300 accuracy of the solution is limited by the convergence criterion. For reference the pressure distribution accuracy of the solution is limited by the convergence criterion. For reference the pressure distribution 301 along the bump surface is given in Fig. 5b. The solutions coincide for the two methods. along the bump surface is given in Figure 5. The solutions coincide for the two methods. ·10 1.2 −1 −2 −3 0.8 −4 0.6 −5 Ux Direct Pressure Direct Ux Iterative Pressure Iterative −6 0.4 0 20 60 80 2 40 1 1.2 1.4 1.6 1.8 CPU-Time [s] X-Coordinate [m] (a) RMS vs. cpu time (b) Pressure along bump surface Figure 5. Consistency results on transonic bump. 303 8.2. Validation using IAPWS-97 real gas state equations 8.2. Validation Using IAPWS-97 Real Gas State Equations 304 The second test case is a nozzle configuration published by Moses and Stein [37]. This The second test case is a nozzle configuration published by Moses and Stein [37]. 305 configuration is frequently used to investigate droplet growth and real gas flow physics, see e.g. This configuration is frequently used to investigate droplet growth and real gas flow physics, see 306 [38]. Since this investigation of real gas flow physics lays the basis for further extension to multiphase e.g., [38]. Since this investigation of real gas flow physics lays the basis for further extension 307 flows, this test case provides an ideal reference for testing the robustness of the implementation. Fig. 6 to multiphase flows, this test case provides an ideal reference for testing the robustness of the 308 shows the overall domain. implementation. Figure 6 shows the overall domain. 309 The inlet total pressure was set to 67760 Pa, while the outlet static pressure was varied. From an 310 outlet pressure of 58000 Pa the pressure was constantly lowered until choking flow configurations 311 were obtained. As for the transonic bump test case before, the k− ω SST turbulence model of [36] was 312 used. p = 67760 Pa p = 58000 Pa 0 Max T = 376.6 K Figure 6. Moses-Stein Nozzle. The inlet total pressure was set to 67,760 Pa, while the outlet static pressure was varied. From an outlet pressure of 58,000 Pa the pressure was constantly lowered until choking flow configurations RMS Residual I nlet Pressure [Pa] Outlet Int. J. Turbomach. Propuls. Power 2020, 5, 20 18 of 32 were obtained. As for the transonic bump test case before, the k-w SST turbulence model of [36] was used. Initially, the computational domain is thermophysically assigned to region 2 of the IAPWS standard, i.e., gaseous state. When lowering the outlet static pressure, the temperature drop due to Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 19 of 32 the acceleration in the nozzle theoretically leads to cells that are assigned to the two-phase region of the IAPWS standard, i.e., region 4. However, it was observed that the commercial code software restricts to region 2 of the IAPWS standard, when running single fluid. Figure 7a shows the centerline temperature along the nozzle geometry from inlet to outlet and demonstrates this behaviour. The dashed line is the solution of the in-house code, the continuous line shows the solution of the commercial code. Extrapolation of region 2 into the two-phase region 4, as done by the commercial p = 67760 Pa p = 58000 Pa 0 Max T = 376.6 K code, leads to considerably lower temperatures in the throat area. The in-house code applies the proper formulations of the IAPWS standard if phase transition is detected. However, it needs to be mentioned that both solutions do not entirely reflect the full reality. Although the properties of the water are clearly better described using the state equations for the two phase region 4, it is not equal to a full two-phase Figure 6. Moses-Stein Nozzle simulation. The assumption of instant thermal equilibrium for example is an oversimplification of the actual problem and prevents an accurate modeling of the phase change. Being able to switch the 313 Initially, the computational domain is thermophysically assigned to region 2 of the IAPWS standard, i.e. state equation in a region where the flow properties change rapidly, however, is a demonstration of 314 gaseous state. When lowering the outlet static pressure, the temperature drop due to the acceleration in the stability of the underlying methods. An option was implemented that allows a restriction of the 315 the nozzle theoretically leads to cells that are assigned to the two-phase region of the IAPWS standard, in-house code to stay within a prescribed region. In case a transition into the two-phase region is 316 i.e. region 4. However it was observed that the commercial code software restricts to region 2 of detected, the values of region 2 are extrapolated. The effect is presented in Figure 7b. The solutions 317 the IAPWS standard, when running single fluid. Fig. 7a shows the centerline temperature along the obtained with the in-house framework match the commercial code. nozzle geometry from inlet to outlet and demonstrates this behaviour. The dashed line is the solution 340 IC IC CC CC Corrdinate [m] Corrdinate [m] (a) Free IAPWS Regions (b) Fixed to Region 2 Figure 7. Influence of the Figure presence 7. Influence of Region of4. the In-house presence code of Region (IC) vs. 4. commercial code (CC). In-house code (IC) vs. commercial code (CC) The following validation was therefore carried out restricting the framework to region 2 of the 319 of the in-house code, the continuous line shows the solution of the commercial code. Extrapolation IAPWS standard. Figure 8 shows the temperature and pressure distribution along the symmetry plane 320 of region 2 into the two-phase region 4, as done by the commercial code, leads to considerably lower for three different outlet pressures. As can be seen the results are in close agreement with commercial 321 temperatures in the throat area. The in-house code applies the proper formulations of the IAPWS code results, over the complete range from subsonic flow (black lines) to transonic flow (gray lines) 322 standard if phase transition is detected. However, it needs to be mentioned that both solutions do until supersonic, chocking flow regime (light gray lines). In Figure 8a a continuous line can be seen for 323 not entirely reflect the full reality. Although the properties of the water are clearly better described the 15,000 Pa case for both numerical algorithms. The pressure curve in Figure 8b does not show this 324 using the state equations for the two phase region 4, it is not equal to a full two-phase simulation. The sharp limit. This is because the value of the pressure in the interior of the domain is a direct outcome 325 assumption of instant thermal equilibrium for example is an oversimplification of the actual problem of the linearized equation system, while temperature is computed as a post-processing field based on 326 and prevents an accurate modeling of the phase change. Being able to switch the state equation in a the values of pressure and enthalpy. In this post-processing step the values for the temperature are 327 region where the flow properties change rapidly however is a demonstration of the stability of the limited to 273.25 K which is the limit of defined temperatures of the IAPWS-97 formulation. 328 underlying methods. An option was implemented that allows a restriction of the in-house code to stay 329 within a prescribed region. In case a transition into the two-phase region is detected, the values of 330 region 2 are extrapolated. The effect is presented in Fig. 7b. The solutions obtained with the in-house 331 framework match the commercial code. 332 The following validation was therefore carried out restricting the framework to region 2 of the IAPWS 333 standard. Fig. 8 shows the temperature and pressure distribution along the symmetry plane for three 334 different outlet pressures. As can be seen the results are in close agreement with commercial code 335 results, over the complete range from subsonic flow (black lines) to transonic flow (gray lines) until Temperature [K] I nlet −0.05 0.00 0.05 0.10 Temperature [K] −0.05 0.00 0.05 Outlet 0.10 Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 20 of 32 336 supersonic, chocking flow regime (light gray lines). In Fig. 8a a continuous line can be seen for the 337 15000Pa case for both numerical algorithms. The pressure curve in Fig. 8b does not show this sharp 338 limit. This is because the value of the pressure in the interior of the domain is a direct outcome of the 339 linearized equation system, while temperature is computed as a post-processing field based on the 340 Int. values J. Turbomach. of pressur Propuls. e and Power enthalpy 2020, 5, 20 . In this post-processing step the values for the temperature are19 limited of 32 to 273.25K which is the limit of defined temperatures of the IAPWS-97 formulation. ·10 IC 58000 Pa IC 58000 Pa CC 58000 Pa CC 58000 Pa IC 40000 Pa IC 40000 Pa CC 40000 Pa CC 40000 Pa IC 15000 Pa IC 15000 Pa CC 15000 Pa CC 15000 Pa Corrdinate [m] Corrdinate [m] (a) Temperature distribution (b) Pressure distribution Figure 8. Allowing region 2 only. In-house code (IC) vs. commercial code (CC). Figure 8. Allowing region 2 only. In-house code (IC) vs. commercial code (CC) 8.3. Two-Stage Radial Compressor Final validation is carried out on the basis of a two-stage radial compressor with a refrigerant as a 342 8.3. Two-stage radial compressor working fluid. Due to the governmental requirements on the global warming potential (GWP), 343 Final validation is carried out on the basis of a two-stage radial compressor with a refrigerant as a replacement of the refrigerant became necessary. This led to a redesign of the components 344 a working fluid. Due to the governmental requirements on the GWP (Global Warming Potential), a accompanied with measurements herein used for the validation of the implemented numerical model. 345 replacement of the refrigerant became necessary. This led to a redesign of the components accompanied The main objective of the measurement campaign was to gain information on the surge line under 346 with measurements herein used for the validation of the implemented numerical model. The main various operating conditions to allow for a controlled operation of the compressor. A vaneless radial 347 objective of the measurement campaign was to gain information on the surge line under various diffuser follows each impeller while the return channel contains a single row of stationary blades to 348 operating conditions to allow for a controlled operation of the compressor. A vaneless radial diffuser reduce the swirl before the second impeller. Figure 9 shows the individual parts. 349 follows each impeller while the return channel contains a single row of stationary blades to reduce the swirl before the second impeller. Figure 9 shows the individual parts. Figure 9. Parts overview. 8.3.1. Measurement Setup Figure 9. Parts Overview The two stage compressor was installed in the laboratory of the Institute for Thermal Energy System and Process Engineering at HSLU in Lucerne. Restriction in the laboratory environment, however, led to operating conditions , which are not in line with the design. This concerns mainly the inlet temperature and pressure levels, which could not be tuned to design conditions. In order to account for this inconsistency, dynamic scaling was applied to match the operating conditions of the design. The rotational speed of the compressor was corrected using Equation (83), with f being the design frequency and f the corrected frequency of the laboratory setup. The mass Temperature [K] −0.05 0.00 0.05 0.10 Pressure [Pa] −0.05 0.00 0.05 0.10 Int. J. Turbomach. Propuls. Power 2020, 5, 20 20 of 32 flow is then corrected based on the non-dimensional mass flow as given Equation (83) for comparison with the design conditions. m ˙ f p M 0 0 1.5 f = q m ˙ = z z = (83) T rT R z p 8.3.2. Numerical Setup The manufacturers choice for the real gas state equations was the standard Peng–Robinson model as given in Equation (2). In order to account for the temperature dependency of the dynamic viscosity, the interacting sphere model given in [39,40] was applied. Furthermore, the thermal conductivity is modeled using the modified Eucken formulation as described in [41]. This decision is mainly based on comparability considerations concerning the CFD computations which have been carried out using a commercial tool during the measurement campaign. Of the available models in the commercial tool, the Peng–Robinson EOS showed the smallest error for the chosen working fluid and operating conditions. Table 1 shows the relative error comparing the chosen thermophysical models against tabulated values from REFPROP [42] over the operating range of the working fluid. A significant error can be seen in laminar conductivity and viscosity; however, this effect is negligible since dominated by the turbulent contributions. The error in the heat capacity was tested on its effect on the results and a change in the value in the range given in the table did not have a significant influence on the final results. Table 1. Relative error for Peng–Robinson and related models. Pressure [bar] Temperature [ C] r c g a l m 2 3 0.48% 3.15% 0.82% 0.13% 10.08% 11.46% 3 12 0.53% 4.04% 1.06% 0.11% 9.43% 11.60% 4 18 0.58% 4.92% 1.33% 0.09% 8.96% 11.63% 5 27 0.46% 5.19% 1.37% 0.02% 8.14% 11.30% 6 33 0.40% 5.64% 1.51% 0.10% 7.49% 10.99% 7 39 0.30% 5.95% 1.58% 0.20% 6.81% 10.53% 8 42 0.31% 6.66% 1.84% 0.24% 6.19% 10.13% 9 48 0.14% 6.69% 1.81% 0.38% 5.42% 9.48% 10 51 0.11% 7.25% 2.01% 0.45% 4.68% 8.89% 11 54 0.06% 7.75% 2.18% 0.51% 3.86% 8.20% 12 60 0.19% 7.39% 1.96% 0.68% 3.10% 7.39% The ideal gas contribution of any internal energy change (Path B in Figure 2) is evaluated using a fourth order polynomial for the specific heat capacity at constant pressure. Turbulence is accounted for by the use of the Menter-SST model [36] with automatic wall treatment. The mesh consists of a total of nine different regions connected with fourteen interfaces for general grid connections and six in-house developed mixing planes, connecting the regions with different periodicities. Details on the implementation of the mixing plane can be found in [43]. The computational domain consists of 4 mio cells with y values around 20. 8.3.3. Global Performance Comparison The complete operating range for a given rotational speed was computed using the above described real gas models by continuously increasing the backpressure until the stability limit was reached. In order to assess the quality of the results, the solutions are compared against commercial Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 22 of 32 373 Global Performance Comparison 374 The complete operating range for a given rotational speed was computed using the above 375 described real gas models by continuously increasing the backpressure until the stability limit was 376 reached. In order to assess the quality of the results, the solutions are compared against commercial 377 code results and measurement data. At the time, the only measurement data available for comparison Int. J. Turbomach. Propuls. Power 2020, 5, 20 21 of 32 378 was the pressure ratio. This is due to the fact that the main objective of the measurement campaign 379 was to find the surge lines to allow for a controlled operation. To validate the implemented procedure, code results and measurement data. At the time, the only measurement data available for comparison 380 the results are further compared with commercial code results by plotting polytropic efficiency (Eq. was the pressure ratio. This is due to the fact that the main objective of the measurement campaign 381 84) and work coefficient (Eq. 85) against flow coefficient (Eq. 86). The circumferential velocity u is was to find the surge lines to allow for a controlled operation. To validate the implemented 382 evaluated at the radius (D ) of the first compressor outlet. procedure, the results are further compared with commercial code results by plotting polytropic efficiency Equation (84) and work coefficient Equation (85) against flow coefficient Equation (86). The circumferential velocity u is evaluated at the radius (D ) of the first compressor outlet. Δh (s − s )(T − T ) p0 2 1 0,2 0,1 η = Δh = (h − h )− (84) p p0 0,2 0,1 Dh T (s s )(T T ) h p0− h 0,2 2 1 0,2 0,1 0,2 0,1 ln h = Dh = (h h ) (84) p p0 0,2 0,1 0,1 0,2 h h 0,2 0,1 ln 0,1 Δh Λ = (85) Dh ( 0u ) L = (85) (u ) 2 ˙ Φ = (86) V 2 u (D ) 2 2 F = (86) u (D ) 2 2 384 Figure 10a shows the results for the whole operating range from choking flow to the stability Figure 10a shows the results for the whole operating range from choking flow to the stability 385 limit. The legend entries IC and CC stand for the in-house code and commercial code respectively. It is limit. The legend entries IC and CC stand for the in-house code and commercial code, respectively. 386 evident that the solutions are in very close agreement, indicating the consistency of the underlying It is evident that the solutions are in very close agreement, indicating the consistency of the underlying 387 thermophysical models. Comparison against measurement data was only available in terms of static thermophysical models. Comparison against measurement data was only available in terms of static 388 pressure ratio, shown in gray. The differences in the results between numerical predictions and pressure ratio, shown in gray. The differences in the results between numerical predictions and 389 measurement setup are dominated by the limitations of the laboratory setup as explained in the section measurement setup are dominated by the limitations of the laboratory setup as explained in the section 390 covering the measurement setup. covering the measurement setup. 1.4 5 −3 0.8 1.2 −4 −5 0.7 −6 0.8 −7 IC Ux 0.6 −8 IC PR CC Ux 0.6 CC PR IC En 2 −9 Measurement CC En −10 0.5 0.4 10 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 Inlet Volume Flow Coefficient Iteration (a) Performance Map Overall System (b) Convergence History Figure 10. Two-stage compressor results. In-house Code (IC) vs commercial code (CC) Figure 10. Two-stage compressor results. In-house Code (IC) vs commercial code (CC). 391 History 392 In order to provide a reliable tool for the design of turbomachinery applications operating in real In order to provide a reliable tool for the design of turbomachinery applications operating in real 393 gas thermophysical regions, accurate numerical solutions are not enough. To allow for a high work gas thermophysical regions, accurate numerical solutions are not enough. To allow for a high work 394 efficiency and reduced resources, numerical stability of the iterative process together with a small efficiency and reduced resources, numerical stability of the iterative process together with a small 395 computational overhead are equally important. Investigations considering the convergence behavior computational overhead are equally important. Investigations considering the convergence behavior 396 have shown promising results as can be seen from Fig 10b. have shown promising results as can be seen from Figure 10b. 397 In order to examine the effects of real gas modeling compared to perfect gas thermophysical models, In order to examine the effects of real gas modeling compared to perfect gas thermophysical 398 cpu time was compared as well. An increase of about 50% compared when switching to real gas models models, cpu time was compared as well. An increase of about 50% compared when switching to real 399 could be observed for this test case setup. A value well within what is observed from commercial code gas models could be observed for this test case setup. A value well within what is observed from 400 results. commercial code results. 8.3.4. Real Gas Effects In this section, we only compare solutions obtained with the in-house framework. The investigation focuses on differences in the results when using real gas models instead of perfect gas assumptions. In order to do so, the characteristic was recomputed using standard models for Polytropic Efficiency (dashed line) Work Coefficient (dotted line) Static Pressure Ratio (continous line) RMS Residual Int. J. Turbomach. Propuls. Power 2020, 5, 20 22 of 32 perfect gas as Newtonian viscosity, constant heat capacity and Prandtl’s law for the conductivity Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 23 of 32 based on average values from the real gas solutions. The results are shown in Figure 11. Lines with square markers correspond to solutions of the perfect gas model (PG) and lines with circular markers correspond to the real gas solution (RG) based on Peng–Robinson state equation as already shown in 401 Real Gas Effects Figure 10a. 402 In With this respect section to polytr we only opic ef compar ficiency e (dashed solutions line), obtained similar rwith esults the are obtained in-house for framework. flow coefficients The 403 investigation above the peak focuses efficiency on dif point. ferences Considering in the results peak when efficiency using , areal differ gas ence models of 2% instead can be of found. perfect While this is surely not negligible, the differences are too small to draw any general conclusions. 404 gas assumptions. In order to do so, the characteristic was recomputed using standard models for Comparing the curves of static pressure ratio (continuous line), it is however clear that when using 405 perfect gas as Newtonian viscosity, constant heat capacity and Prandtl’s law for the conductivity perfect gas state equation lower back pressure is necessary for the same volume flow coefficient. 406 based on average values from the real gas solutions. The results are shown in Fig. 11. Lines with Using the perfect gas model, it was not possible to reach the same pressure levels as when using the 407 square markers correspond to solutions of the perfect gas model (PG) and lines with circular markers real gas model. The compressor stability limit was detected at much smaller levels of back pressure 408 correspond to the real gas solution (RG) based on Peng-Robinson state equation as already shown in and the flow coefficient dropped much faster for the perfect gas when increasing the back pressure. Fig. 10a. 1.4 0.8 1.2 IP 1 IP 2 0.7 0.6 0.8 IC RG 0.5 0.6 IC PG Measurement 0.4 0.4 0.6 0.7 0.8 0.9 1 Inlet Volume Flow Coefficient Figure 11. Performance Map: Real Gas (RG) vs. Perfect Gas (PG) Figure 11. Performance map: real gas (RG) vs. perfect gas (PG). 410 With respect to polytropic efficiency (dashed line), similar results are obtained for flow coefficients Further investigation of the marked regions (IP1 in blue, IP2 in red) show the differences. 411 above the peak efficiency point. Considering peak efficiency, a difference of∼ 2% can be found. While IP1 corresponds to peak efficiency for both setups, IP2 are choking conditions. The operating points 412 this is surely not negligible, the differences are too small to draw any general conclusions. Comparing of the CFD simulation were set by fixing the value of the pressure at the outlet. As can be seen from 413 the curves of static pressure ratio (continuous line), it is however clear that when using perfect gas IP1 (peak efficiency), to obtain the same flow coefficient, the perfect gas model needs a much smaller 414 state equation lower back pressure is necessary for the same volume flow coefficient. Using the perfect static pressure ratio between inlet and outlet. In fact, the outlet pressure of the perfect gas simulation 415 gas model, it was not possible to reach the same pressure levels as when using the real gas model. The is roughly 13% lower than for real gas. This is in contrast to IP2 (choking conditions) where pressure 416 compressor stability limit was detected at much smaller levels of back pressure and the flow coefficient ratio and flow coefficient are in close agreement. This can, however, be explained by Figure 12c. 417 dropped much faster for the perfect gas when increasing the back pressure. Under choking conditions, the flow is much closer to ideal gas behavior (Compressibility z = 1), 418 Further investigation of the marked regions (IP1 in blue, IP2 in red) show the differences. IP1 thus the effect of real gas modeling is smaller. 419 corresponds to peak efficiency for both setups, IP2 are choking conditions. The operating points Comparing the exit flow angles from Figure 12a the general orientation between real gas and 420 of the CFD simulation were set by fixing the value of the pressure at the outlet. As can be seen from perfect gas in stage 1 do not show a significant difference. This is independent whether choking 421 IP1 (peak efficiency), to obtain the same flow coefficient, the perfect gas model needs a much smaller or peak efficiency is considered. Increasing deviation can be found when comparing the results in 422 static pressure ratio between inlet and outlet. In fact, the outlet pressure of the perfect gas simulation stage 2, see Figure 12b. An explanation can be found when plotting the average compressibility factor 423 is roughly 13% lower than for real gas. This is in contrast to IP2 (choking conditions) where pressure in streamwise direction given in Figure 12c, showing the increasing deviation from ideal gas state 424 ratio through and the flow compr coeffi essor cient . This are in deviation close agrfr eement. om perfect Thisgas can modeling howeverassumptions be explainedbecome by Fig. mor 12c.eUnder and 425 choking conditions, the flow is much closer to ideal gas behavior (Compressibility z=1), thus the effect more dominant with increasing back pressure (denser gas), explaining the increased difference in 426 of polytr real opic gas modeling efficiencyis for smaller the peak . efficiency point. 427 Comparing the exit flow angles from Fig. 12a the general orientation between real gas and perfect gas 428 in stage 1 do not show a significant difference. This is independent whether choking or peak efficiency 429 is considered. Increasing deviation can be found when comparing the results in stage 2, see Fig. 12b. 430 An explanation can be found when plotting the average compressibility factor in streamwise direction 431 given in Fig. 12c, showing the increasing deviation from ideal gas state through the compressor. This 432 deviation from perfect gas modeling assumptions become more and more dominant with increasing 433 back pressure (denser gas), explaining the increased difference in polytropic efficiency for the peak 434 efficiency point. Polytropic Efficiency (dashed line) Work Coefficient (dotted line) Static Pressure Ratio (continous line) Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 24 of 32 Int. J. Turbomach. Propuls. Power 2020, 5, 20 23 of 32 Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 24 of 32 0.8 0.8 0.95 0.8 0.6 0.8 0.95 0.6 0.9 0.6 0.6 0.4 0.4 0.9 IP2 Ideal 0.4 0.4 IP2 Real 0.85 IP2 Ideal 0.2 0.2 IP1 Ideal IP2 Real 0.85 IP2 0.2 0.2 IP1 Real IP1 Ideal IP1 IP2 0 0 IP1 Real 0.8 IP1 0 20 40 0 20 40 60 80 0 0.5 1 1.5 2 0 0 0.8 ◦ ◦ 0 20 40 0 20 40 60 80 Flow Angle [ ] Flow Angle [ 0 0.5 1 1.5 2 Streamwise location Flow Angle [ ] Streamwise location (a) Stage 1 (b) Stage 2 (c) Compressibility (a) Stage 1 (c) Compressibility (b) Stage 2 Figure 12. Flow angle and compressibility Figure 12. Flow angle and compressibility. A normalized dynamic pressure was computed in order to compare the dynamic contribution A normalized dynamic pressure was computed in order to compare the dynamic contribution A normalized dynamic pressure was computed in order to compare the dynamic contribution between real and perfect gas computations as given in Eq. 87. between real and perfect gas computations as given in Eq. 87. between real and perfect gas computations as given in Equation (87). p p −−p ¯ p ¯ p p Ψ = (87) Ψ = (87) Y 0 = (87) p ¯− p ¯ p ¯− p ¯ p ¯0 p ¯ 435 Fig. 13 shows the results for IP1 (Peak efficiency) at sections close to the impeller outlets. For stage 1 Figure 13 shows the results for IP1 (Peak efficiency) at sections close to the impeller outlets. 435 Fig. 13 shows the results for IP1 (Peak efficiency) at sections close to the impeller outlets. For stage 1 436 the results are comparable, differences however become present in the second stage. This is especially For stage 1 the results are comparable, differences however become present in the second stage. This is 436 the results are comparable, differences however become present in the second stage. This is especially 437 true for the suction side (SS) of the splitter blade (SB), where the dynamic contribution is much lower. especially true for the suction side (SS) of the splitter blade (SB), where the dynamic contribution is 437 true for the suction side (SS) of the splitter blade (SB), where the dynamic contribution is much lower. much lower. (a) Real Gas, Stage 1 (b) Ideal Gas, Stage 1 (a) Real Gas, Stage 1 (b) Ideal Gas, Stage 1 (c) Real Gas, Stage 2 (d) Ideal Gas, Stage 2 Figure 13. Normalized Total Pressure Isolines, IP1 (c) Real Gas, Stage 2 (d) Ideal Gas, Stage 2 439 These differences in terms of compressor dynamics when switching between real and perfect Figure 13. Normalized total pressure isolines, IP1. 440 gas model explain the importance of an accurate representation of the thermophysical properties and 439 These differences in terms of compressor dynamics when switching between real and perfect 441 derivation of the conservation equations. Using perfect gas models would lead to errors in the design These differences in terms of compressor dynamics when switching between real and perfect 440 gas model explain the importance of an accurate representation of the thermophysical properties and 442 of blade geometries because the flow physics are not captured accurately. gas model explain the importance of an accurate representation of the thermophysical properties and 441 derivation of the conservation equations. Using perfect gas models would lead to errors in the design derivation of the conservation equations. Using perfect gas models would lead to errors in the design 443 CONCLUSIONS 442 of blade geometries because the flow physics are not captured accurately. of blade geometries because the flow physics are not captured accurately. 444 The present work summarizes the necessary steps during the derivation and implementation of 9. Conclusions 443 CONCLUSIONS 445 Navier-Stokes equations into a pressure based numerical framework when operating in non-ideal gs 446 conditions. Detailed analysis is given concerning the discretization of the conservation equations The present work summarizes the necessary steps during the derivation and implementation of 444 The present work summarizes the necessary steps during the derivation and implementation of 447 when Navier using –Stokes realequations gas state equations. into a pressur This e based is extended numerical furthermor framework e to the when development operating of in boundary non-ideal 445 Navier-Stokes equations into a pressure based numerical framework when operating in non-ideal gs 448 conditions and iterative procedures for the update of various flow quantities. A transonic bump gs conditions. Detailed analysis is given concerning the discretization of the conservation equations 446 conditions. Detailed analysis is given concerning the discretization of the conservation equations 449 case when was using chosen real to gas demonstrate state equations. the consistency This is extended of thefurthermor modifications. e to the The development added iterative of boundary solution 447 when using real gas state equations. This is extended furthermore to the development of boundary 450 pr conditions ocedures for and the iterative thermophysical procedur states es forreduce the update to the of same various results flow as the quantities. direct algebraic A transonic relations bump for 448 conditions and iterative procedures for the update of various flow quantities. A transonic bump 449 case was chosen to demonstrate the consistency of the modifications. The added iterative solution 450 procedures for the thermophysical states reduce to the same results as the direct algebraic relations for Span Span Compressibility Factor Compressibility Factor Inlet Inlet Impeller 1 Impeller 1 Diffuser 1 Crossover Diffuser 1 Bend Crossover Impeller 2 Bend Diffuser Impeller 2 2 Outlet Diffuser 2 Outlet Int. J. Turbomach. Propuls. Power 2020, 5, 20 24 of 32 case was chosen to demonstrate the consistency of the modifications. The added iterative solution procedures for the thermophysical states reduce to the same results as the direct algebraic relations for perfect gas. The real gas models were then tested using two different state equations, once using the IAPWS-97 formulation for water and once using the Peng–Robinson cubic state equation. Comparison against the commercial code results was carried out using the IAPWS-97 formulation on a transonic nozzle configuration. It was found that commercial code extrapolates the gaseous equations into the two-phase region. Applying the same treatment to the in-house code then allowed further comparison and a close agreement between the two frameworks was found. Cubic state equations were then applied to an industrial type compressor setup and compared against commercial code results and measurement data. It could be proven that the implemented procedure is not only able to accurately replicate commercial code results but also leads to a stability of the iterative scheme that goes beyond the capabilities of the used commercial tool. Comparison against perfect gas assumptions justify the increased complexity of the implementation. It could be demonstrated, based on overall performance prediction and local flow features, that the assumption of perfect gas shows its limitation under certain operating conditions, especially away from choking conditions. At moderate backpressure levels, however, i.e., high volume flow rate, the assumption of perfect gas might be a viable option. Especially during preliminary design phase, the possibility to use perfect gas assumption without great loss of accuracy can lower the computational costs and therefore increase the throughput. For a thorough analysis of the geometry, real gas state equations are, however, a necessity. Author Contributions: The extension of a pressure based coupled numerical algorithm to real gas flow physics was mainly driven by L.H. as it is a part of his PhD. The final target of the thesis is a framework, able at predicting two-phase flow physics for steam turbines. It is however based on the excessive work done by L.M. and M.D. in the development of such a coupled framework and thus also supported by these two authors. E.C. and D.M.V. provided valuable background information considering thermophysical relations and the various model assumptions when moving to real gas flow physics and their effect on turbomachinery design. All authors contributed to writing and editing this paper. All authors have read and agreed to the published version of the manuscript. Funding: The authors gratefully acknowledge the financial contribution provided by the Swiss National Science Foundation (SNF: Grant-Number 175900). Conflicts of Interest: The authors declare no conflict of interest. Nomenclature Abbreviations A 4  4 block matrix a Matrix element b Source term CC Commercial code c Speed of sound [m/s] c Specific heat at constant pressure [J/kgK] c Specific heat at constant volume [J/kgK] d Face to cell distance vector [m] D Inverted 4  4 block matrix e Internal energy [J/kg] F Body force [kg/ms ] bi f Frequency [1/s] h Enthalpy [J/kg] IC In-house code I P Investigation point Ma Mach number Int. J. Turbomach. Propuls. Power 2020, 5, 20 25 of 32 M Molar mass [kg/mol] m Mass flow [kg/s] MR Maxwell relation PG Perfect Gas p Pressure [Pa] q Heat [J/kg] R Universal gas constant [J/kgK] R Specific gas constant [J/kgK] RG Real Gas S Surface vector [m ] f i s Entropy [J/kgK] T Temperature [K] t Time [s] U Velocity flux [m /s] u Velocity x-component [m/s] v Velocity y-component [m/s] w Velocity z-component [m/s] V Volume [m ] v Specific volume [1/kg] z Compressibility Greeks a Thermal expansivity [1/K] b Compressibility [1/Pa] G Fundamental derivative g Heat capacity ratio h Polytropic efficiency L Work coefficient l Thermal conductivity [W/mK] m Dynamic viscosity [kg/ms] Y Normalized dynamic pressure [-] r Density [kg/m ] 2 2 t Stress tensor [kg/m s ] i j F Flow coefficient w Accentric factor Subscripts b Boundary value c Cell value cr Critical point values f Surface values ¥ Ideal gas state nb Neighbour matrix element p Diagonal matrix element r Reference state 0 Stagnation conditions Previous iteration value Superscripts n Time step index u Coupling to velocity x-component v Coupling to velocity y-component w Coupling to velocity z-component p Coupling to pressure Correction to given quantity Int. J. Turbomach. Propuls. Power 2020, 5, 20 26 of 32 Appendix A. Entropy Change Using the first and second law of thermodynamics, the change in specific internal energy can be described as given in Equation (A1). de = Tds pdv (A1) Using the fundamental derivative, ds can be rewritten as given in Equation (A2). ¶s ¶s ds = dT + dv (A2) ¶T ¶v v T Inserting Equation (A2) into Equation (A1) and using Maxwell relations, the general form given in Equation (A3) can be found. ¶ p de = c dT + T p dv (A3) ¶T Appendix B. Partial Thermal Derivatives Table A1. Constant Pressure. Derivative General ¶r ¶h vc ¶r ¶T v ¶T 1 ¶h c Table A2. Constant Volume. Derivative General ¶h c b + (1 Ta)v ¶ p ¶s ¶T ¶T b ¶ p a ¶h 2 c vT ¶T Table A3. Constant Temperature. Derivative General ¶v ¶s ¶ p bv ¶v ¶h aT1 ¶v T Int. J. Turbomach. Propuls. Power 2020, 5, 20 27 of 32 Table A4. Constant Enthalpy. Derivative General ¶r 2 Ta a rb ¶ p Table A5. Constant Entropy. Derivative General ¶ p c v p b c ¶r Appendix C. Real Gas Internal Energy Using Peng–Robinson State Equation As stated earlier in Equation (6) the internal energy of any real gas can be expressed as given in Equation (A4). Z Z Z v T v ¶ p ¶ p e(T, p) = T p dv + [c ] dT + T p dv +e(T , p ) (A4) v r r ¶T ¶T v T v r r ¥ v v T T | {z } | {z } | {z } ¶ p Using Peng–Robinson [12] state equation the derivative can be expressed as given in ¶T Equation (A7). ¶ p ¶ R T aa(T) = (A5) 2 2 ¶T ¶T v b v + 2bv b v v ¶a(T) a2a R ¶T = (A6) 2 2 v b v + 2bv b k T aa R T T s c = + (A7) 2 2 v b v + 2bv b ¶ p Hence the term T p can be written as given in Equation (A8). ¶T ¶ p aa(T)(1 + k) T p = (A8) 2 2 ¶T v + 2bv b Inserting Equation (A8) into term A of Equation (A4) leads to Equation (A9). " !# ¶ p aa(T )(1 + k) 2b T p dv = p tanh (A9) ¶T b + v r 2b This allows to rewrite the internal energy for any state T, v as given in Equation (A12) aa(T )(1 + k) 2b e(T, v) = p tanh (A10) b + v 2b r aa(T)(1 + k) 2b p tanh (A11) b + v 2b + c dT R (T T ) + e(T , v ) (A12) p s r r r r Int. J. Turbomach. Propuls. Power 2020, 5, 20 28 of 32 Appendix D. Real Gas Internal Energy Using Benedict–Webb–Rubin (BWR) State Equation As shown in Appendix C, the expression needed for the real gas contribution of the state change ¶ p is the term T p. BWR [15] state equation is given as: ¶T C c r 2 0,s s 2 3 6 2 g r p = rR T + B R T A r + (b R T a ) r + a a r + 1 + g r e (A13) s s s s s s s s 0,s 0,s 2 2 T T ¶ p Hence can be written is: ¶T ¶ p C c r 2 0,s 2 3 2 g r = rR + (B R + 2 )r + b R r 2 (1 + g r )e (A14) s 0,s s s s s 3 3 ¶T T T Hence the requested term reads as: ¶ p C c r 2 0,s 2 2 2 6 2 g r T p = r + A r + ar a a r 1 + g r e (A15) 0,s s s s 2 2 ¶T T T The remaining procedure just follows what is shown in Appendix C. Appendix E. Computation of Entropy Using the first law of thermodynamics and Maxwell relations, the entropy can be expressed as shown in Equation (A16). ¶ p c ds(v, T) = dv + dT (A16) ¶T T Integration along paths according to Figure 2 leads to Equation (A17). Z Z Z v T v ¶ p c ¶ p s(T, v) s(T , v ) = dv + dT R ln(T/T ) + dv (A17) r r s r T T ¶T T ¶T v T v r r ¥ v v Peng–Robinson This is an example of the implemented formulation when using Peng–Robinson [12] state equation. R a (T) s(T, v) s(T , v ) = dv (A18) r r T 2 2 v b v + 2bv b R a (T) + dv (A19) 2 2 v b v + 2bv b + dT R ln(T/T ) (A20) s r " !# a (T ) 2b s(T, v) s(T , v ) = R ln(v b) p tanh (A21) r r s b + v 2b " !# a (T) 2b + R ln(v b) p tanh (A22) b + v 2b + dT R ln(T/T ) (A23) s r r Int. J. Turbomach. Propuls. Power 2020, 5, 20 29 of 32 a (T) 2b s(T, v) s(T , v ) = R ln(v b) p tanh (A24) r r s b + v 2b a (T ) 2b R ln(v b) + p tanh (A25) s r b + v 2b r + dT R ln(T/T ) (A26) s r Appendix F. Specific Heat Capacity The specific heat at constant volume is defined as given in Equation (A27). ¶e c = (A27) ¶T Using the definition of internal energy of Equation (6) this equation can be rewritten as given in Equation (A30) ¶ ¶ p c = c T p dv (A28) v v ¶T ¶T ¶ ¶ p ¶ p = c T dv (A29) ¶T ¶T ¶T v T ¶ p = c T dv (A30) ¶T v T ¶ p The partial derivative can be obtained through differentiation of the state equation. ¶T Values for the ideal gas contribution (c ) can be found through various models, e.g., constant value or NASA’s polynomial expression [18]. To obtain an expression for the specific heat at constant pressure, the general statement of Equation (A32) is used. ¶q ¶e c c = (A31) p v ¶T ¶T p v ¶s ¶s = T (A32) ¶T ¶T p v Using the substantial derivative of the entropy and specific volume ¶s ¶s ds = dT + dv (A33) ¶T ¶v v T ¶v ¶v dv = dT + d p (A34) ¶T ¶ p p T Substituting Equation (A34) in Equation (A33) results in Equation (A35). When expressing the change in entropy through fundamental derivative a second expression is found as given in Equation (A36). " # ¶s ¶s ¶v ¶s ¶v ds = + dT + d p (A35) ¶T ¶v ¶T ¶v ¶ p v T p T ¶s ¶s ds = dT + d p (A36) ¶T ¶ p p T Int. J. Turbomach. Propuls. Power 2020, 5, 20 30 of 32 With Equations (A35) and (A36) it is possible to express the derivative of entropy with respect to temperature at constant pressure as given in Equation (A37). ¶s ¶s ¶s ¶v = + (A37) ¶T ¶T ¶v ¶T p v T p Substituting Equation (A37) in Equation (A32) allows to rewrite c as given in Equation (A39). ¶s ¶v c c = T (A38) p v ¶v ¶T T p ¶ p ¶v c = c + T Maxwell relation (A39) p v ¶T ¶T v p For the density based state equations it is, however, advantageous to rewrite those formulations using triple product rule. Details on the partial derivatives can be found in Appendix B. ¶v ¶v ¶ p = Triple product rule (A40) ¶T ¶ p ¶T p T v ¶ p ¶T c = c T  (A41) p v ¶ p ¶v Appendix G. Speed of Sound ¶ p ¶ p = v (A42) ¶r ¶v ¶ p ¶s = v Triple product rule (A43) ¶s ¶v v p ¶T ¶ p = v Maxwell relation (A44) ¶v ¶T s s ¶T ¶s ¶ p ¶s = v Triple product rule (A45) ¶s ¶v ¶s ¶T v T T p c ¶ p ¶T ¶s = v c = T (A46) c ¶T ¶v ¶T v p v p ¶ p = v Triple product rule (A47) c ¶v References 1. Starzmann, J.; Casey, M.M.; Mayer, J.F.; Sieverding, F. Wetness loss prediction for a low pressure steam turbine using computational fluid dynamics. Proc. Inst. Mech. Eng. Part A J. Power Energy 2014, 228, 216–231. 2. Colonna, P.; Rebay, S.; Harinck, J.; Guardone, A. Real-gas effects in ORC turbine flow simulations: Influence of thermodynamic models on flow fields and performance parameters. In Proceedings of the European Community on Computational Methods in Applied Sciences (ECCOMAS) CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, Egmond aan Zee, The Netherlands, 5–8 September 2006; Delft University of Technology: Delft, The Netherlands, 2006. 3. Pohl, S.; Jarczyk, M.; Pfitzner, M.; Rogg, B. Real gas CFD simulations of hydrogen/oxygen supercritical combustion. Prog. Propuls. Phys. 2013, 4, 583–614. Int. J. Turbomach. Propuls. Power 2020, 5, 20 31 of 32 4. Odabaee, M.; Sauret, E.; Hooman, K. CFD Simulation of a Supercritical Carbon Dioxide Radial-Inflow Turbine, Comparing the Results of Using Real Gas Equation of Estate and Real Gas Property File. Appl. Mech. Mater. 2016, 846, 85–90. 5. Boncinelli, P.; Rubechini, F.; Arnone, A.; Cecconi, M.; Cortese, C. Real Gas Effects in Turbomachinery Flows: A CFD Model for Fast Computations. In Proceedings of the ASME Turbo Expo 2003, Collocated with the 2003 International Joint Power Generation Conference, Atlanta, GA, USA, 16–19 June 2003; American Society of Mechanical Engineers: New York, NY, USA, 2003; pp. 1103–1112. 6. Travis, J.; Koch, D.P.; Xiao, J.; Xu, Z. Real-gas equations-of-state for the GASFLOW CFD code. Int. J. Hydrogen Energy 2013, 38, 8132–8140. 7. Lucas, C.; Rusche, H.; Schroeder, A.; Koehler, J. Numerical investigation of a two-phase CO ejector. Int. J. Refrig. 2014, 43, 154–166. 8. Ameli, A.; Afzalifar, A.; Turunen-Saaresti, T.; Backman, J. Effects of real gas model accuracy and operating conditions on supercritical CO compressor performance and flow field. J. Eng. Gas Turbines Power 2018, 140, 062603. 9. Mangani, L.; Buchmayr, M.; Darwish, M. Development of a Novel Fully Coupled Solver in OpenFOAM: Steady State Incompressible Turbulent Flows. Numer. Heat Transf. Part B Fundam. 2014, 66, 1–20. doi:10.1080/10407790.2014.894448. 10. Mangani, L.; Buchmayr, M.; Darwish, M. Development of a Novel Fully Coupled Solver in OpenFOAM: Steady State Incompressible Turbulent Flows in Rotational Reference Frames. Numer. Heat Transf. Part B Fundam. 2014, 66, 526–543. doi:10.1080/10407790.2014.894371. 11. Mangani, L.; Moukalled, F.; Darwish, M. An OpenFOAM pressure-based coupled CFD solver for turbulent and compressible flows in turbomachinery applications. Numer. Heat Transf. Part B Fundam. 2016, 69 413–431. 12. Peng, D.Y.; Robinson, D.B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. 13. Cooper, J.; Dooley, R. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam; IAPWS R7-97(2012); The International Association for the Properties of Water and Steam: Lucerne, Switzerland, 2007. 14. Young, J.B. An Equation of State for Steam for Turbomachinery and Other Flow Calculations. J. Eng. Gas Turbines Power 1988, 110, 1–7, doi:10.1115/1.3240080. 15. Benedict, M.; Webb, G.B.; Rubin, L.C. An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures I. Methane, ethane, propane and n-butane. J. Chem. Phys. 1940, 8, 334–345. 16. Van der Waals, J.D. Over de Continuiteit van den Gas-en Vloeistoftoestand. Ph.D. Thesis, Leiden University: Leiden, The Netherlands, 1873. 17. Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 1972, 27, 1197–1203. 18. McBride, B.J.; Zehe, M.J.; Gordon, S. NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species; E-13336; National Aeronautics and Space Administration: Washington, DC, USA, September 2002. 19. Baltadjiev, N.D. An Investigation of Real Gas Effects in Supercritical CO Compressors. Ph.D. Thesis, Massachusetts Institute of Technology: Cambridge, MA, USA, 2012. 20. Baltadjiev, N.D.; Lettieri, C.; Spakovszky, Z.S. An investigation of real gas effects in supercritical CO centrifugal compressors. J. Turbomach. 2015, 137, 091003. 21. Barton, I.E. Comparison of SIMPLE-and PISO-type algorithms for transient flows. Int. J. Numer. Methods Fluids 1998, 26, 459–483. 22. Anjorin, V.; Barton, I. Removal of temporal and under-relaxation terms from the pressure-correction equation of the SIMPLE algorithm. Int. J. Fluid Dyn. 2001, 5, 59–75. 23. Kissling, K.; Springer, J.; Jasak, H.; Schutz, S.; Urban, K.; Piesche, M. A coupled pressure based solution algorithm based on the volume-of-fluid approach for two or more immiscible fluids. In Proceedings of the V European Conference on Computational Fluid Dynamics ECCOMAS CFD, Lisbon, Portugal, 14–17 June 2010. 24. Cerne, G.; Tiselj, I.; Petelin, S. Upgrade of the VOF Method for the Simulation of the Dispersed Flow. In Proceedings of the ASME 2000 Fluids Engineering Division Summer Meeting, Boston, MA, USA, 11–15 June 2000. Int. J. Turbomach. Propuls. Power 2020, 5, 20 32 of 32 25. Hutchinson, B.; Raithby, G. A multigrid method based on the additive correction strategy. Numer. Heat Transf. Part A Appl. 1986, 9, 511–537. 26. Keller, S.; Cordazzo, J.; Hinckel, P.H.; Maliska, C.R. Additive correction multigrid method applied to diffusion problems with unstructured grids. In Proceedings of the 10th Brazilian Congress of Thermal Sciences and Engineering, Rio de Janeiro, RJ, Brazil, 29 November–3 December 2004. 27. Darwish, M.; Sraj, I.; Moukalled, F. A coupled incompressible flow solver on structured grids. Numer. Heat Transf. Part B Fundam. 2007, 52, 353–371. 28. Moukalled, F.; Aziz, A.A.; Darwish, M. Performance comparison of the NWF and DC methods for implementing high-resolution schemes in a fully coupled incompressible flow solver. Appl. Math. Comput. 2011, 217, 5041–5054. 29. Darwish, M.; Moukalled, F. A fully coupled Navier-Stokes solver for fluid flow at all speeds. Numer. Heat Transf. Part B Fundam. 2014, 65, 410–444. 30. Darwish, M.; Abdel Aziz, A.; Moukalled, F. A coupled pressure-based finite-volume solver for incompressible two-phase flow. Numer. Heat Transf. Part B Fundam. 2015, 67, 47–74. 31. Mangani, L.; Buchmayr, M.; Darwish, M.; Moukalled, F. A fully coupled OpenFOAM® solver for transient incompressible turbulent flows in ALE formulation. Numer. Heat Transf. Part B Fundam. 2017, 71, 313–326. 32. Darwish, M.; Mangani, L.; Moukalled, F. Implicit boundary conditions for coupled solvers. Comput. Fluids 2018, 168, 54–66. 33. Rhie, C.; Chow, W.L. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J. 1983, 21, 1525–1532. 34. Darwish, M.; Sraj, I.; Moukalled, F. A coupled finite volume solver for the solution of incompressible flows on unstructured grids. J. Comput. Phys. 2009, 228, 180–201. 35. Demirdžic, ´ I.; Lilek, Ž.; Peric, ´ M. A collocated finite volume method for predicting flows at all speeds. Int. J. Numer. Methods Fluids 1993, 16, 1029–1050. 36. Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994, 32, 1598–1605. 37. Moses, C.A.; Stein, G.D. On the growth of steam droplets formed in a Laval nozzle using both static pressure and light scattering measurements. J. Fluids Eng. 1978, 100, 311–322. 38. Starzmann, J.; Hughes, F.R.; White, A.J.; Grübel, M.; Vogt, D.M. Numerical investigation of boundary layers in wet steam nozzles. J. Eng. Gas Turbines Power 2017, 139, 012606. 39. Chung, T.H.; Lee, L.L.; Starling, K.E. Applications of kinetic gas theories and multiparameter correlation for prediction of dilute gas viscosity and thermal conductivity. Ind. Eng. Chem. Fundam. 1984, 23, 8–13. 40. Chung, T.H.; Ajlan, M.; Lee, L.L.; Starling, K.E. Generalized multiparameter correlation for nonpolar and polar fluid transport properties. Ind. Eng. Chem. Res. 1988, 27, 671–679. 41. Reid, R.C.; Prausnitz, J.M.; Poling, B.E. The Properties of Gases and Liquids; McGraw Hill Book Co.: New York, NY, USA, 1987. 42. Lemmon, E.W.; .; Bell, I.H.; Huber, M.L.; McLinden, M.O. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0; National Institute od Standards and Technology: Gaithersburg, MD, USA, 2018. doi:10.18434/T4JS3C. 43. Hanimann, L.; Mangani, L.; Casartelli, E.; Mokulys, T.; Mauri, S. Development of a novel mixing plane interface using a fully implicit averaging for stage analysis. J. Turbomach. 2014, 136, 081010. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Turbomachinery, Propulsion and Power Multidisciplinary Digital Publishing Institute

Real Gas Models in Coupled Algorithms Numerical Recipes and Thermophysical Relations †

Loading next page...
 
/lp/multidisciplinary-digital-publishing-institute/real-gas-models-in-coupled-algorithms-numerical-recipes-and-bwS90WaiZt

References (45)

Publisher
Multidisciplinary Digital Publishing Institute
Copyright
© 1996-2020 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). Terms and Conditions Privacy Policy
ISSN
2504-186X
DOI
10.3390/ijtpp5030020
Publisher site
See Article on Publisher Site

Abstract

International Journal of Turbomachinery Propulsion and Power Article Real Gas Models in Coupled Algorithms Numerical Recipes and Thermophysical Relations 1, ‡,§ 1,§ 1 2 Lucian Hanimann * , Luca Mangani , Ernesto Casartelli , Damian M. Vogt and Marwan Darwish Competence Center for Fluid Mechanics and Hydro Machines, Lucerne University of Science and Arts, 6048 Horw, Switzerland; luca.mangani@hslu.ch (L.M.); ernesto.casartelli@hslu.ch (E.C.) Institut für Thermische Strömungsmaschinen und Maschinenlaboratorium, Universität Stuttgart, 70569 Stuttgart, Germany; damian.vogt@itsm.uni-stuttgart.de Departement of Mechanical Engineering, American University of Beirut, Beirut 110236, Lebanon; darwish@aub.edu.lb * Correspondence: lucian.hanimann@hslu.ch; Tel.: +41-41-349-34-58 † This paper is an extended version of our paper published in Proceedings of the European Turbomachinery Conference ETC13, Lausanne, Switzerland, 8–12 April 2019; Paper No. 218. ‡ Current address: Technik & Architektur, University of Applied Sciences Lucerne, Technikumstrasse 21, 6048 Horw, Switzerland. § These authors contributed equally to this work. Received: 13 May 2020; Accepted: 29 July 2020; Published: 3 August 2020 Abstract: In the majority of compressible flow CFD simulations, the standard ideal gas state equation is accurate enough. However, there is a range of applications where the deviations from the ideal gas behaviour is significant enough that performance predictions are no longer valid and more accurate models are needed. While a considerable amount of the literature has been written about the application of real gas state equations in CFD simulations, there is much less information on the numerical issues involved in the actual implementation of such models. The aim of this article is to present a robust implementation of real gas flow physics in an in-house, coupled, pressure-based solver, and highlight the main difference that arises as compared to standard ideal gas model. The consistency of the developed iterative procedures is demonstrated by first comparing against results obtained with a framework using perfect gas simplifications. The generality of the developed framework is tested by using the parameters from two different real gas state equations, namely the IAPWS-97 and the cubic state equations state equations. The highly polynomial IAPWS-97 formulation for water is applied to a transonic nozzle case where steam is expanded at transonic conditions until phase transition occurs. The cubic state equations are applied to a two stage radial compressor setup. Results are compared in terms of accuracy with a commercial code and measurement data. Results are also compared against simulations using the ideal gas model, highlighting the limitations of the later model. Finally, the effects of the real gas formulations on computational time are compared with results obtained using the ideal gas model. Keywords: real gas; numerics; compressor 1. Introduction CFD continues to be a critical tool in the design and development of turbomachinery equipment. It is used in the evaluation of the performance of complex multistage turbines at design and off-design points, and in the optimization of single blades. Furthermore, the continuous increase in available computational resource has allowed for a substantial improvement in the quality of simulations as compared to a decade ago. While the majority of compressible flows still uses the perfect gas Int. J. Turbomach. Propuls. Power 2020, 5, 20; doi:10.3390/ijtpp5030020 www.mdpi.com/journal/ijtpp Int. J. Turbomach. Propuls. Power 2020, 5, 20 2 of 32 model, it is becoming more common to work on applications where real gas effects are no longer negligible. Example where these situations occur include low pressure stages of steam turbines [1] or ORC turbines [2]. Many authors have investigated the impact of various real gas models, mainly focusing on possible achievements in the accuracy of the results and the impact on computational time due to the increased complexity of the thermophysical models, see, e.g., [2–7]. The later point, the increase in computational time, has, however, received less attention. While the complexity of real gas relations is probably the most obvious cause, other factors also play a role, such as the iterative nature of the solution procedures used to update thermophysical properties. This non-linearity is often addressed by using lookup tables, see, e.g., [5]. In this case the user needs to know in advance the range of the thermophysical states, which is not always available at the start of a simulation. This is further complicated by the fact that, prior to convergence, some intermediate solutions can run out of the tabulated range. For such situations, a stabilization criterion is needed. Additional issues arise concerning the resolution of the generated lookup table. A low resolution may lead to a faster numerical procedure, but can diminish the quality of the results. The high sensitivity of the results to lookup table resolution was described in [8]. Other techniques to address the weakness of lookup tables include the use of structured grids for the lookup tables and the use of various interpolation methods whether bi-linear or bi-cubic ([5]). In this paper the real gas relations are incorporated algebraically, i.e., without using lookup tables, into an in-house, pressure based, all-speed, fully coupled solver [9–11]. The main intention is therefore to contribute to the understanding of the necessary modifications considering a compressible all-Mach solver when using real gas state equations. Starting from the basic real gas relation, the conservation equations and boundary conditions and the derivation of the developed numerical techniques are presented in details with a main focus on retaining the performance and robustness of the coupled solver. Furthermore, the computational cost of the use of real gas relations is thoroughly investigated. Although the developed methods to deal with real gas relations are general, two specific relations are evaluated in more detail. It is for the reason that these models are later used for the validation and it is considered a contribution to help the reader to apply the procedure to his desired real gas model. The highlighted models are the cubic state equations, e.g., [12] and the IAPWS database for water at any state [13]. It will be shown that the increase in computational time often does not justify the problems arising from lookup table formulations. Replacement with tabulated data is, however, straight forward. Three test cases are used to validate and evaluate the presented numerical approach. In the first case, the consistency of the approach is demonstrated. This is accomplished by showing that the solution of additional iterative procedures needed for real gas formulations recovers the results of the perfect gas model when applied to a 2D transonic bump case. In the second test case, the IAPWS-97 formulation [13] is applied to a nozzle configuration with the operating range varying from supersonic to transonic flow regimes. The expansion for the transonic conditions yield to operating conditions in the two-phase region. Results are compared to those obtained from a commercial code. The last test case is an industrial type two stage compressor at one operational speed for the full range from choking flow to the stability limit, the cubic state equations are used as they are the manufacturer ’s choice for this application. Results are then compared against measurement data. Furthermore, for cases two and three, an evaluation of the computational cost of the real gas implementations is presented. 2. Thermodynamic Modeling of Real Gas Real gases are compositions of materials at thermophysical states under which its molecules are dense enough so that intermolecular forces can no longer be neglected. This is in contrast to the theory of ideal gases. For ideal gases, the molecules are sparsely distributed in space and these forces are not taken into account. It is this assumption that leads to the well known ideal gas state equation. Int. J. Turbomach. Propuls. Power 2020, 5, 20 3 of 32 The fundamental derivative of gas dynamics is a key parameter in understanding real-gas flows. It is a measurement for the change in speed of sound through a change in density, its definition is provided in Equation (1). r ¶c G = 1 + (1) c ¶r g + 1 While this expression reduces to a constant value of G = for perfect gas and is always greater than one, this does not apply for real gas. A special example are Bethe–Zel’dovich–Thompson (BZT) fluids, which are heavy polyatomic compositions. For these materials, G can even have negative values when approaching the critical point. The effect is an inverted gas dynamic behavior, where, contrary to perfect gas, the Mach number decreases as the flow expands. 2.1. General Comment on Real Gas Modeling The modeling of real gases in a numerical framework is accompanied with a variety of additional challenges. This is a direct result of the increased complexity of the state equations compared to the simple ideal gas state equations. The non-linearities introduced require iterative updates of a variety of quantities, which will be a major part of what follows. It also requires a profound analysis of all underlying equations and thermophysical relationships. Entropy, internal energy and enthalpy are a function of two variables for the general case of a real gas. Since cubic equations are used in the final validation (in which p = f (T, v)), the equations and partial derivatives are given using T, v as primary input variable. However, this does not limit the generality of the approach and just requires a careful derivation of what is exemplified below. p, T, v are interchangeable through the state equation, but depending on its formulation it may simplify the analysis if the input variables are changed to, e.g., T, p. 2.2. Real Gas State Equations Various models have been proposed to better account for intermolecular forces when encountering denser fluids, each having its own justification depending on the operting conditions and the working fluid. Analytical representation for superheated and two-phase steam is given by the model of Young [14] in virial form or for water at any state in the form given through IAPWS [13]. When modeling hydrocarbons, cryogenic fluids or refrigerants, well known virial equations are those of Benedict–Webb–Rubin (BWR) [15] and its modifications. Another group of real gas state equations is represented in the form of cubic state equations. The original formulation was introduced by Van der Waals [16] but many other models have emerged to overcome some of its drawbacks. Modifications have been applied to improve the phase transition, see Soave et al. [17], or the prediction of liquid phase density and accuracy in the critical region by Peng et al. [12]. 2.2.1. Cubic State Equations The cubic state equation, used for the final validation case, is based on the formulation provided by Peng et al. [12] and therefore given in Equation (2). T and p are temperature and pressure at the cr cr critical point and w being the accentric factor. R T aa p = (2) 2 2 v b v + 2bv b 2 2 T 0.45724R T s cr a = a(T) = 1 + k 1 a = (3) T p cr cr R T s cr b = 0.0778 k = 0.37646 + 1.54226w 0.26993w (4) cr Int. J. Turbomach. Propuls. Power 2020, 5, 20 4 of 32 2.2.2. IAPWS-97 The IAPWS-97 [13] formulations provide accurate state equations for the properties of water at any state. The equations are highly polynomial functions and provided for five different subregions described in Figure 1. Region 1 is the liquid phase and region 3 accounts for the transition from liquid to gaseous phase above the critical point. Region 2 and 5 are gaseous states and region 4 is the two-phase region below the critical point. Figure 1. IAPWS regions. While this set of equations allows for a highly accurate representation of the properties of water, they can lead to a massive increase in computational costs. Even by the use of provided backwards functions to eliminate at least partially the iterative procedures the increase in computational time is still considerable. The usual approach is therefore the use of lookup tables generated using IAPWS formulations. 2.3. Caloric State Equations and Definition of Entropy The numerical algorithm requires values of internal energy, entropy and enthalpy at a certain point. This can be the case during the update of variables that are not a direct outcome of a conservation equation such as temperature or density. However, it may also be required when updating boundary conditions, such as given stagnation conditions at the inlet. Using real gas state equations, the definitions of the caloric state equations for enthalpy and internal energy are functions of two intensive properties, e.g., temperature and pressure, and are no longer only a function of temperature alone. The derivation of generalized state equations is therefore briefly described here on the basis of internal energy. The definition for the change in internal energy is given in Appendix A and the final expression is shown in Equation (5). The integral form of this differential equation can be found using an integration along paths as shown in Figure 2. ¶ p de = c dT + T p dv (5) ¶T The value of e(T, p) is computed through integration along lines of constant temperature ( A, C) and constant specific volume (B) using a given reference state e(T , p ). The final expression for internal r r energy is thus given in Equation (6), formulations for enthalpy and entropy (see Appendix E) can be found using the same basic idea. Int. J. Turbomach. Propuls. Power 2020, 5, 20 5 of 32 Z Z Z v T v ¶ p ¶ p e(T, p) = T p dv + [c ] dT + T p dv +e(T , p ) (6) v r r ¶T ¶T v T v r r ¥ v v T T | {z } | {z } | {z } It is obvious that such a representation of the integral value can be understood as the sum of the ideal gas contribution plus a correction to adjust for real gas effects. As stated earlier, ideal gas behavior is encountered when intermolecular forces do not affect the properties of the fluid. This is always valid if the molecules are separated by an infinite distance. In other words the density goes towards zero or specific volume to infinity. Hence the integration along line B accounts for the ideal gas contribution, and integration along the lines A, C corrects the value to correspond to the real gas change in the state variable. The ideal gas heat capacity in term B can be described with a variety of models. The simplest is the assumption of calorically perfect gas (i.e., c = const.), other options are Sutherland models or NASA’s polynomial approach [18]. As stated in the introduction, the final validation of the coupled solution strategy for real gas CFD is based on an industrial sized two-stage radial compressor. It is the manufacturer ’s choice to use the Peng–Robinson [12] state equation. Appendix C therefore gives a detailed derivation of Equation (6) for this type of state equation. To provide an example of the generality of the approach the derivations for internal energy for the Benedict–Webb–Rubin (BWR) [15] state equation are given in addition in Appendix D. T, p A ¥ T , p r r Figure 2. Integration path. 2.4. Stagnation Conditions Evaluation of stagnation conditions given static values or vice versa is a necessity frequently encountered in engineering practice. The relation is usually based on the assumption of an isentropic state change, leading to the well known expressions for incompressible and compressible (assuming ideal gas) flows as given in Equations (7) and (8) exemplarilly for pressure. p = p + rjjujj (7) g1 g 1 p = p 1 + Ma (8) These relationships usually provide the basis in the derivation of stagnation conditions, however they do not apply to real gas. In the derivation of Equation (8) it is assumed that the isentropic exponent is a constant, a simplification which is not valid for real gas. Some authors ([19,20]) argue that the variations are, however, small and Equation (8) can be applied with minimal loss of accuracy. While this may be a justified assumption for a variety of applications, its generality is questionable. The generalized approach used in the present article is therefore explained in later sections when describing boundary conditions and the update of domain internal quantities. Int. J. Turbomach. Propuls. Power 2020, 5, 20 6 of 32 2.5. Specific Heat Capacity As already introduced when discussing the stagnation conditions, the heat capacity and heat capacity ratio do no longer follow the well known ideal gas relations. The heat capacity ratio can vary for different thermophysical states and requires a generalized approach. The derivation used in the implementation is therefore given in detail in Appendix F. The final expression is given in Equation (9) ¶ p ¶T c = c T (9) p v ¶ p ¶v 2.6. Mach Number The Mach number is a frequently used quantity for the post-processing and analysis of numerical results. In the following we will therefore provide a general formulation for the speed of sound. The definition of the Mach number is given by Equation (10). jjujj Ma = (10) The speed of sound for any compressible medium is given by Equation (11). ¶ p c = (11) ¶r Detailed derivation of the partial differential of pressure with respect to density under isentropic conditions can be found in Appendix G. The final expression is given in Equation (12). ¶ p ¶ p = v (12) ¶r c ¶v The partial differentials can again be found in Appendix B. 3. Numerical Environment Segregated pressure-based solution algorithms [21,22] based on the SIMPLE family of algorithms have been applied successfully over the past decades. The popularity of these algorithms is partly due to their low memory requirement based on the sequential method in solving the RANS equations. Another reason is the relatively simple implementation. On the downside, these algorithms suffer from convergence breakdown for large-scale problems. The reason for this issue is that the naturally strong pressure-velocity coupling is numerically decoupled when solving the momentum and continuity equation sequentially. Owing to the developments in computer architecture, coupled solution strategies for pressure-based algorithms have regained interest in recent years [23,24]. These algorithms solve the momentum and continuity equation together, i.e., coupled. The numerical framework used in this article is based on the coupled solution strategy presented in [9–11]. The coupled system of linearized momentum and continuity equations for each cell of the computational domain takes the form given in Equation (13). Int. J. Turbomach. Propuls. Power 2020, 5, 20 7 of 32 2 32 3 2 32 3 2 3 u p u p uu uv uw uu uv uw u a a a a u a a a a u b p p p p p nb c nb nb nb nb v p 6 76 7 6 v p76 7 6 7 vu vv vw vu vv vw v a a a a v a a a a v b 6 76 7 6 76 7 6 7 p p p p p nb nb nb nb nb 6 76 7 + 6 76 7 = 6 7 (13) w p w p wu wv ww wu wv ww w 4a a a a 54w 5 4a a a a 54w 5 4b 5 p p p p nb c nb nb nb nb pu pv pw p p pu pv pw p p p p a a a a p b a a a a p nb c p p p p nb nb nb nb | {z } | {z } p nb The final system of linearized equations is solved using algebraic multigrid algorithms on the basis of an additive correction approach [25,26] with a block-ILU smoother. Figure 3 highlights the main body of the coupled solution strategy. The blocks of the flow chart will also define the structure of the present article. In Section 5, the conservation equations are revisited with a focus on the difference when integrating real gas thermophysics. This is followed by a detailed analysis of boundary conditions mainly focusing on the implementation of total conditions, see Section 6. The theoretical part is finalized showing the chosen strategies for the update of flow properties in Section 7. n n1 Update Boundary Conditions t = t + Dt Coupled System Sequential Systems Momentum-Equation Energy Equation A p 3 A p Assemble Terms Turbulence Equations Compute D uu uv uw u p a a a a Assemble Total Conditions T, v = f (h, p) vu vv vw v p a a a a Property update wu wv ww w p Continuity Equation a a a a Correct U and p pu pw pw p p a a a a Assemble Terms Figure 3. Flow chart of the algorithm. 4. Governing Equations The above described numerical framework solves the RANS equations as: ¶r ¶ru Continuity Equation + = 0 (14) ¶t ¶x ¶ru u ¶t ¶ru i j ¶ p i j Momentum Equation + + f = 0 (15) ¶t ¶x ¶x ¶x j i j ¶rh ¶ p ¶ Energy Equation + ru h + q u t = 0 (16) i 0 i j i j ¶t ¶t ¶x t in the momentum and energy equation is the deviatoric stress tensor given as: i j ¶u ¶u j 2 ¶u i l t = m + m d (17) i j i,j ¶x ¶x 3 ¶x j i l q in the energy equation accounts for diffusive heat flux modeled through Fourier ’s approach as: ¶T q = l (18) ¶x 5. Numerical Discretization of Governing Equations The governing equations for compressible fluid dynamic analysis are the momentum, continuity, energy and turbulence equations. As described in Equation (13) the momentum and continuity equations are solved simultaneously to retain the pressure-velocity coupling, followed by a sequential solution of the energy and turbulence equations. A complete derivation of the governing equations, Int. J. Turbomach. Propuls. Power 2020, 5, 20 8 of 32 schemes and used solvers would go beyond the scope of this article. Implementation details of the numerical framework is best described in [9–11,27–32]. The main focus will be the integration of a generalized approach for real gas physics. Partial recapitulation of the derivations concerning conservation equations and boundary conditions will be given where necessary. 5.1. Momentum Equation The momentum equation remains unchanged compared to the derivation for perfect gas. It is, however, important in order to introduce the coupled solver framework and therefore given in Equation (19). ¶ru V F + (r u u ) + p t S = 0 (19) bi å f f j f i f f i j f i ¶t F accounts for body forces and S is the normal vector of the face with the length of its bi f i area. The implicit contributions of the single terms of the momentum equation to the coupled, block-structured coefficient matrix are assembled in what is marked as a red box in Figure 3. The coupling of the momentum equation to the pressure based continuity equation is dominated by the pressure gradient and assembled in a , a and a , shown in red in Equation (13). This is one of the u p v p w p major advantages of coupled solution strategies. The pressure gradient is treated implicitly and thus updated during inner iterations. Segregated solution strategies consider the pressure gradient in the momentum equation as a constant contribution to the source when solving the momentum equations. 5.2. Continuity Equation For pressure based numerical frameworks, the continuity equation is used to obtain an equation for the pressure which serves as a constraint to fulfill the conservation of mass. Equation (20) shows the continuity equation in its semi discrete form with U = u S being the velocity flux across a face. f f i f i dr V + r U = 0 (20) å f f dt The temporal derivative is expanded using the general assumption that the density can be expressed as a function of pressure and enthalpy. This is a procedure that has been described in [7]. The linearization of this term is therefore given as shown in Equation (21). Expressions have to be found for the thermophysical derivatives based on the chosen real gas state equation. This derivation will be exemplified in a later section. dr ¶r dh ¶r d p = + (21) dt ¶h dt ¶ p dt p h Equation (21) differs from perfect gas computations where it is usually assumed that the change of density with respect to time can be simply expressed as given in Equation (22). dr R T = (22) dt dt The second term of the continuity equation, the convective contribution, is discretized using Newton linearization, see Equation (23). Variables indexed with  are previous iteration values while variables without this index represent current iteration values. r U = r U + r U r U (23) f f f f f f f f Int. J. Turbomach. Propuls. Power 2020, 5, 20 9 of 32 To avoid checkerboarding pressure fields on collocated arrangements and to introduce the necessary implicit dependency on pressure together with a coupling of the continuity equation to the velocity, the velocity flux of the new time step (U ) is computed using the Rhie–Chow [33] interpolation technique. 0 1 B ! ! C B C ¶ p ¶ p B C U = U D S (24) B C f f f ,i,j f i B ¶x ¶x C j j @ f fA | {z } | {z } A B U is the velocity flux on the face, computed using adjacent cell values through linear interpolation, contributing to the coefficients a , a and a in the coefficient matrix A . This couples the continuity pu pv pw p equation to the velocity field, shown in green in Figure 3. The pressure gradient A is evaluated directly on the face , while term B is evaluated based on cell center values and then interpolated on the face. The term D is the inverted 3  3 coefficient sub-matrix A (given in green in Figure 3) after f p3 cell-to-face interpolation. The density r in Equation (23) is computed using a first order linearization, to obtain the necessary implicit dependency on pressure. This is in contrast to the derivation for ideal gas, where density is substituted with pressure by simply using the state equation. ¶r ¶r r = r + (h h ) + ( p p ) (25) f f f f f f ¶h ¶ p p, f h, f Inserting Equations (24) and (25) in Equation (23) and using Equation (21) in Equation (20) leads to the final continuity equation in pressure form shown in Equation (26). ¶r dh ¶r d p V + + r U f f ¶h dt ¶ p dt | {z } p h m ˙ 0 0 1 1 ! ! ¶ p ¶ p @ ¯ @ A A + r U D S å f f f ,i,j f i ¶x ¶x j j f f | {z } m ˙ ¶r ¶r + U (h h ) + ( p p ) r = 0 (26) å f f f f f f ¶h ¶ p p f h f | {z } m ˙ The term r U introduces an implicit coupling of the continuity equation to the velocity and pu pv pw p p will therefore contribute to a , a and a . Three terms are implicit contributions to a , namely c c c c ¶ p ¶r ¶r d p r D S , U p and V . å å f p f f i f f f f ¶x ¶ p ¶ p dt h f h To summarize, the difference in the derivation of the pressure equation for real gas compared to idealized gas is listed below: dr • Linearization of transient term dt • Linearization of the density of the actual time step r ¶r ¶r • Partial derivatives and ¶h ¶ p p h Int. J. Turbomach. Propuls. Power 2020, 5, 20 10 of 32 5.2.1. Thermophysical Partial Derivatives in the Continuity Equation The derivation of the partial derivatives given above is crucial for the stability of the solver. In this ¶r section the derivation of is given as an example. Expansion of the expression is first done ¶ p through triple product rule. ¶r ¶h ¶r = (27) ¶ p ¶ p ¶h h v p | {z } | {z } B1 B2 Part B1 can extended as follows ¶h ¶s = T + v dh = Tds + vd p (28) ¶ p ¶ p v v ¶v ¶v ¶s = T + v MR: = (29) ¶T ¶T ¶ p s s v ¶s ¶v = T + v Triple product rule (30) ¶T ¶s v T c b ¶s c ¶v b v v = T + v = and = (31) T a ¶T T ¶s a v T c b a = + (1 Ta)v c = c vT (32) v p a b In this derivation a is the thermal expansion coefficient and b is the compressibility as given in Equations (33) and (34). ¶ p 1 ¶T a = (33) v ¶ p ¶v 1 ¶v b = (34) v ¶ p Part B2 can be extended as follows ¶r ¶r ¶T = Product rule (35) ¶h ¶T ¶h p p p ¶r 1 = (36) ¶T c 1 ¶v 1 = (37) ¶T c a 1 = (38) v c The final expression can therefore be written as given in Equation (39) ¶r ¶h ¶r Ta a = = rb (39) ¶ p ¶ p ¶h c h v p Int. J. Turbomach. Propuls. Power 2020, 5, 20 11 of 32 5.3. Energy Equation The definition of the conservation equation for total internal energy is given in Equation (40). ¶re ¶ + ru e + u p + q u t dV = 0 (40) j 0 j j i i,j ¶t ¶x For compressible flows it is, however, more convenient to express the energy equation in terms of total enthalpy as h = e + p/r. Equation (40) is therefore rewritten as Equation (41). 0 0 ¶rh ¶ p ¶ + ru h + q u t dV = 0 (41) i 0 i j i j ¶t ¶t ¶x Compared to ideal gas state equations, the elementary changes are in the modeling of the heat flux. Using Fourier ’s law to relate the heat flux to the local temperature gradient with l being the conductivity of the working fluid, Equation (42). is derived. Z I ¶q ¶T dV = q n d A  q S = l S (42) i i å i f f i å f f i ¶x ¶x W ¶W i i f f This term would lead to a completely explicit contribution to the source of the energy equation, diminishing the stability of the iterative solution procedure. What is needed is an algebraic relation for the dependency of temperature in terms of enthalpy to obtain an implicit formulation of the heat conduction. For calorically perfect ideal gases, this poses no challenge, since enthalpy can be expressed as h h = c (T T ). A generalized approach, as considered in this publication, allows no such 2 1 p 2 1 simplification. To recover at least partly the implicit contribution a first order linearization of the enthalpy is introduced as given in Equation (43). ¶h ¶h h = h + (T T) +  ( p p) (43) ¶T ¶ p p T | {z } Expression A can be found in Appendix B and is simply the definition of the specific heat capacity c . B is zero because the pressure is kept constant during the solution of the energy equation. An expression for the temperature of the actual time step can therefore be found as given in Equation (45) h h T = T + (44) = T + (45) Substituting the temperature T in Equation (42) using Equation (45) allows to rewrite the Laplacian contribution as given in Equation (46). Term a introduces the implicit contribution providing the needed stability of the discretization. ¶T ¶T f ¶h l S = l S + S å f f i å f f i å f i ¶x ¶x c ¶x i i p f i f f f f f f | {z } ¶c h p + l S (46) å f f i ¶x p f | {z } b Int. J. Turbomach. Propuls. Power 2020, 5, 20 12 of 32 At convergence term a and b vanish due to their dependency on the correction of enthalpy, hence the heat conduction is modeled accurately. Part b can usually be neglected since the gradient of the heat capacity is often rather small. 6. Boundary Conditions Considering the implementation of boundary conditions for a coupled framework, it is crucial to identify any coupling between velocity and pressure as well as implicit dependencies on the variable that is solved for. A detailed derivation of the implemented boundary conditions can be found in [32,34]. While these articles give a detailed insight in the numerical discretization of a variety of boundary conditions, they assume perfect gas state equations. The focus of this section will therefore remain on the difference when integrating real gas state equations. 6.1. Total Conditions The use of total conditions at the inlet is a very common boundary condition. A detailed derivation is therefore given in what follows. As already shown, in a coupled framework, the momentum and continuity equation are assembled into one system of equations. In order to properly account for the boundary contribution, each term of the conservation equations has to be analyzed with respect to the constraint given by the user. This is described in detail in [32] for ideal gas. The derivation of total conditions in this section is given by first providing the algorithm for consistent total to static conversion. The algorithm uses the common assumption of an isentropic deceleration to stagnation conditions but without all the simplifications valid only for perfect gas. This is followed by the discretization of the Navier–Stokes equations itself. With the exception of the pressure gradient term, this is in close agreement to what is given in [32]. The main body of this section will therefore be a detailed revision of the discretization of the pressure gradient term when dealing with real gas state equations. 6.1.1. Isentropic Total-to-Static Conversion As already mentioned above, the computation of static values using given total conditions can no longer be obtained using well-known algebraic relations. It is, however, possible to find an iterative solution procedure to accurately compute these quantities. The algorithm used is a Newton–Raphson root finding method. The roots are the two underlying constraints of the total to static conversion. Root 1 as given in Equation (47) is the isentropic constraint, while root 2 as given in Equation (48) is simply using the definition of total enthalpy. f (x , x ) = s(T , v ) s(T, v) = 0 (47) 1 1 2 0 0 f (x , x ) = (h(T , v ) 0.5jjujj ) h(T, v) = 0 (48) 2 1 2 0 0 In the following the Newton–Raphson root finding method is explained as a reference. This procedure is later applied to a variety of solution strategies used for finding iterative solutions for thermophysical changes in real gas. The procedure starts with a first order Taylor expansion of the roots. n n n n ¶ f (x , x ) ¶ f (x , x ) 1 1 n+1 n+1 n n 1 2 n+1 n 1 2 n+1 n f (x , x ) = f (x , x ) + (x x ) + (x x ) (49) 1 1 1 2 1 2 1 1 2 2 ¶x ¶x 1 2 n n n n ¶ f (x , x ) ¶ f (x , x ) 2 2 n+1 n+1 n n 1 2 n+1 n 1 2 n+1 n f (x , x ) = f (x , x ) + (x x ) + (x x ) (50) 2 2 2 2 1 2 1 1 1 2 ¶x ¶x 1 2 From Equations (47) and (48) it is known that these equations are equal to zero. Hence we can reformulate them as shown in Equations (51) and (52). Int. J. Turbomach. Propuls. Power 2020, 5, 20 13 of 32 n n n n ¶ f (x , x ) ¶ f (x , x ) 1 1 n n 1 2 1 2 0 = f (x , x ) + Dx + Dx (51) 1 1 2 1 2 ¶x ¶x 1 2 n n n n ¶ f (x , x ) ¶ f (x , x ) 2 2 n n 1 2 1 2 0 = f (x , x ) + Dx + Dx (52) 2 1 2 1 2 ¶x ¶x 1 2 The corrections Dx and Dx are defined as: 1 2 n+1 n Dx = x x (53) n+1 n Dx = x x (54) 2 2 This is system of equation is then rewritten in matrix form as: 2 3 n n n n ¶ f (x , x ) ¶ f (x , x ) 1 1 " # " # 1 2 1 2 n n 6 7 Dx f (x , x ) 1 1 ¶x ¶x 1 2 1 2 6 7 = (55) n n n n n n 4 5 ¶ f (x , x ) ¶ f (x , x ) 2 2 Dx f (x , x ) 1 2 1 2 2 2 1 2 | {z } | {z } ¶x ¶x 1 2 | {z } x By inverting matrix A the system can be solved as shown in Equation (56) 2 3 n n n n ¶ f (x , x ) ¶ f (x , x ) " # " # 1 2 1 2 1 2 n n 6 7 Dx 1 f (x , x ) 1 ¶x ¶x 1 1 2 6 1 1 7 = (56) n n n n n n n n n n n n 4 5 n n ¶ f (x ,x ) ¶ f (x ,x ) ¶ f (x ,x ) ¶ f (x ,x ) ¶ f (x , x ) ¶ f (x , x ) Dx 1 2 2 1 f (x , x ) 1 2 1 2 1 2 1 2 2 1 2 1 2 1 2 2 1 2 ¶x ¶x ¶x ¶x 1 2 1 2 | {z } | {z } ¶x ¶x 1 1 The variables x , x are then updated with the correction and the procedure repeated until certain 1 2 convergence criterion is achieved. Applied to the roots given in Equations (47) and (48) with x = T and x = p, the following derivatives have to be found using the chosen real gas state equation. n n n n ¶ f (x , x ) ¶h ¶ f (x , x ) ¶h 1 1 1 2 1 2 = and = (57) ¶x ¶T ¶x ¶v 1 2 v T n n n n ¶ f (x , x ) ¶s ¶ f (x , x ) ¶s 2 2 1 2 1 2 = and = (58) ¶x ¶T ¶x ¶v 1 2 v T ¶h ¶s ¶s can be found in Appendix B. is equal to and can be replaced with ¶T ¶T ¶v v v T ¶ p ¶h by using Maxwell relations. is found as given in Equation (59). ¶T ¶v v T ¶h ¶e ¶ p = + v + p h = e + pv (59) ¶v ¶v ¶v v T T ¶e In a second step, is derived from Equation (6) as given in Equation (60). ¶v ¶e ¶ p = T p (60) ¶v ¶T T v Hence Equation (59) can be written as given in Equation (61). The derivatives of the pressure are easily found through the chosen state equation. ¶h ¶ p ¶ p = T + v (61) ¶v ¶T ¶v v v T Int. J. Turbomach. Propuls. Power 2020, 5, 20 14 of 32 6.1.2. Pressure Gradient Term As already pointed out, the discretization of the Navier–Stokes equation for given total conditions is best described following what is written in [32]. When integrating real gas state equations, the only difference is in the discretization of the pressure gradient term of the momentum equation. Using Green–Gauss theorem to rewrite the volume integral value of the pressure gradient as a surface integral, this term is written as given in Equation (62). An expression for the value of the pressure at the boundary face is therefore needed. The requirement is a strong coupling to face adjacent cell values of pressure and velocity, leading to an implicit contribution to the coefficient matrix. This ensures a robust integration of boundary conditions into the linearized equation system. Z I ¶ p dV = pd A  p S (62) å f f i ¶x W ¶W Using the well known strong coupling between pressure and velocity, a first order linearization of the pressure with respect to the velocity flux U is introduced as given in Equation (63). b indexes values evaluated at boundary faces. Again, the velocity flux is defined as U = u S with S being a i f i f i face normal vector with the length of the area of the face. As for the continuity equation this leads to an implicit contribution to the cell velocity and an implicit coupling to the pressure field of the continuity equation. ¶ p p = p + (U U ) (63) b b b b ¶U ¶ p What is needed is an expression for the gradient and the value of the velocity flux of the new ¶U time step U . The value for the velocity flux is computed using Rhie–Chow interpolation, given in Equation (64). ! ! ! ¶ p ¶ p U = U D S (64) b b,i,j ib ¶x ¶x j j b c The gradient of the pressure at the boundary is approximated as given in Equation (65). ¶ p p p b b c = (65) ¶x jdj With p being the cell pressure value andjdj the normal distance from face to cell center. Inserting Equation (63) in Equation (64) leads to the final definition for the boundary pressure value given in Equation (66). ¶ p p = p + qU + qD p + qF U (66) c c b b ¶U DjS j f i D = (67) jdj ¶ p q = 1 + D (68) ¶U ¶ p ¶ p ˜ ˜ F = D p + D U + D S (69) f i b b ¶U ¶x b c ¶ p The only unknown is therefore the gradient of pressure with respect to the velocity flux . ¶U This is usually computed using the algebraic relations for the static to stagnation conversion, given in Equation (70) for incompressible cases and Equation (71) for compressible cases. Int. J. Turbomach. Propuls. Power 2020, 5, 20 15 of 32 p = p + 0.5rjjujj (70) g1 g 1 p = p 1 + M (71) As explained before, direct algebraic relations that describe the change in pressure of a fluid particle when brought to rest isentropically are not explicitly available for the general case of a real gas. ¶ p To approximate the gradient the definition given in Equation (72) is used. ¶U ¶ p ¶ p ¶u p + D p p ¶u = = lim (72) ¶U ¶u ¶U D p!0 u( p = p + D p) u( p) ¶U The static value of the pressure is computed using an iterative scheme given the total conditions based on T and p as described above. In a second step, the velocity is computed that corresponds to 0 0 a slightly increased pressure level p = p + D p, following an isentropic change. Hence the roots are: f (x , x ) = s(T , p ) s(T , v ) = 0 (73) 2 0 0 2 2 1 1 f (x , x ) = p p(T , v ) = 0 (74) 2 1 2 2 2 2 In addition to the derivatives given in the section dealing with total to static conversion the partial derivatives of the pressure are needed. ¶ f (x , x ) ¶ p ¶ f (x , x ) ¶ p 2 1 2 2 1 2 = = (75) ¶x ¶T ¶x ¶v 1 2 v T Using the values for T and v the value of the static enthalpy h = f (T , v ) is computed according 2 2 2 2 to Equation (76). The calculation of internal energy follows the previously explained integration along paths. h(T , v ) = e(T , v ) + p v (76) 2 2 2 2 2 2 The velocity magnitude corresponding to these updated values is further found through the definition of total enthalpy as given in Equation (77). jju( p )jj = 2.0(h(T , p ) h(T , v )) (77) 2 0 0 2 2 The value for the derivative of pressure with respect to velocity is then found by substituting these results into Equation (72). 7. Update of Flow Properties As shown in Figure 3, the outcome of the conservation equations for momentum, continuity and energy are updated values for velocity, pressure and total enthalpy. There are, however, many other flow properties that are not directly obtained but which have to be computed based on these known values. This may be because they are intrinsically necessary for the solution procedure, as, e.g., density or viscosity, or because they are valuable quantities for post-processing such as Mach-number or entropy. This section explains their computation when operating in regions, where real gas effects can no longer be neglected. 7.1. Update of Temperature and Density Given the velocity, the pressure and total enthalpy from the conservation equations, it is possible to update the temperature and density independent of the state equation. First, the static Int. J. Turbomach. Propuls. Power 2020, 5, 20 16 of 32 enthalpy is computed using the given total enthalpy h and the velocity based on the solution of the RANS equations. h = h 0.5u u (78) i i The value for temperature and density are computed using the Newton–Raphson method. The roots of the method are: f (x , x ) = h h(T, v) = 0 (79) 1 1 2 f (x , x ) = p p(T, v) = 0 (80) 2 1 2 The rest of the Newton–Raphson root finding procedure can be done in analogy to Section 6.1.1. 7.2. Total Pressure and Total Temperature In contrast to the problem of specified total conditions for boundary conditions as described in Section 6.1, finding stagnation values for static conditions is usually required for the remaining cells and boundary faces of the computational domain. As described, no explicit relation exists for the static to stagnation conversion, asking for an iterative procedure for the isentropic state change. Consequently, Newton–Raphson root finding method is used. Given the static values of temperature and specific volume, the two roots therefore are given in Equations (81) and (82). With x and x being 1 2 the requested stagnation values T and v . 0 0 f (x , x ) = s(T, v) s(T , v ) = 0 (81) 1 1 2 0 0 f (x , x ) = (h(T, v) + 0.5u u ) h(T , v ) = 0 (82) 2 1 2 i i 0 0 As written above, having the roots of the Newton–Raphson method defined, the solution is in close analogy to what is given in Section 6.1.1. With the stagnation values for T and v the total 0 0 pressure is easily reconstructed using the chosen state equation. 8. Validation and Results This section provides an assessment of the generality and accuracy of the implemented numerical solution procedure when dealing with real gas state equations. Initially, the well-known transonic bump of Peric [35] is used. This test case is used to test the consistency of the implemented procedures when perfect gas state equations are used. The expressions and iterative algorithms given above have to return the same solutions as when using the mentioned algebraic relations only valid for perfect gas. Validity of the implemented procedures is given through comparison of the RMS-residuals. The second test case is a transonic nozzle configuration operated using steam based on the IAPWS-97 formulation [13]. Results are compared to a commercial code. To demonstrate the generality of the approach to any type of state equations, cubic state equations are used for a two-stage radial compressor. Results are compared against commercial code and measurement data. 8.1. Consistency of the algorithms for Perfect Gas State Equation The transonic bump test case of Peric [35] is given in Figure 4. The thickness to chord ratio of the bump is 10% and the computational domain consists of 53,000 quadrilateral grid elements. The inlet boundary conditions is set to a total pressure of 1.4 bar and a total temperature of 783K, turbulence intensity is 0.1% and the eddy viscosity ratio is set 2. The outlet static pressure is 1bar. The working fluid is air with the turbulence modeled using the k-w-SST model of [36]. Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 18 of 32 289 8.1. Consistency of the algorithms for perfect gas state equation 290 The transonic bump test case of Peric [35] is given in Fig. 4. The thickness to chord ratio of the 291 bump is 10% and the computational domain consists of 53’000 quadrilateral grid elements. The inlet 292 boundary conditions is set to a total pressure of 1.4bar and a total temperature of 783K, turbulence 293 intensity is 0.1% and the eddy viscosity ratio is set 2. The outlet static pressure is 1bar. The working Int. J. Turbomach. Propuls. Power 2020, 5, 20 17 of 32 fluid is air with the turbulence modeled using the k− ω-SST model of [36]. Figure 4. Transonic Bump Figure 4. Transonic bump. 295 To prove the consistency of the implemented generalized approach, the results should reduce to To prove the consistency of the implemented generalized approach, the results should reduce 296 the same output as when directly using perfect gas assumptions, when deriving the fundamental to the same output as when directly using perfect gas assumptions, when deriving the fundamental 297 equations. Fig. 5a shows the comparison between the direct solution and the iterative procedure for equations. Figure 5a shows the comparison between the direct solution and the iterative procedure for 298 the ideal gas case. The small difference results from the fact that the generalized approach uses an the ideal gas case. The small difference results from the fact that the generalized approach uses an 299 iterative procedure to evaluate isentropic conditions and to update temperature and density, hence the iterative procedure to evaluate isentropic conditions and to update temperature and density, hence the 300 accuracy of the solution is limited by the convergence criterion. For reference the pressure distribution accuracy of the solution is limited by the convergence criterion. For reference the pressure distribution 301 along the bump surface is given in Fig. 5b. The solutions coincide for the two methods. along the bump surface is given in Figure 5. The solutions coincide for the two methods. ·10 1.2 −1 −2 −3 0.8 −4 0.6 −5 Ux Direct Pressure Direct Ux Iterative Pressure Iterative −6 0.4 0 20 60 80 2 40 1 1.2 1.4 1.6 1.8 CPU-Time [s] X-Coordinate [m] (a) RMS vs. cpu time (b) Pressure along bump surface Figure 5. Consistency results on transonic bump. 303 8.2. Validation using IAPWS-97 real gas state equations 8.2. Validation Using IAPWS-97 Real Gas State Equations 304 The second test case is a nozzle configuration published by Moses and Stein [37]. This The second test case is a nozzle configuration published by Moses and Stein [37]. 305 configuration is frequently used to investigate droplet growth and real gas flow physics, see e.g. This configuration is frequently used to investigate droplet growth and real gas flow physics, see 306 [38]. Since this investigation of real gas flow physics lays the basis for further extension to multiphase e.g., [38]. Since this investigation of real gas flow physics lays the basis for further extension 307 flows, this test case provides an ideal reference for testing the robustness of the implementation. Fig. 6 to multiphase flows, this test case provides an ideal reference for testing the robustness of the 308 shows the overall domain. implementation. Figure 6 shows the overall domain. 309 The inlet total pressure was set to 67760 Pa, while the outlet static pressure was varied. From an 310 outlet pressure of 58000 Pa the pressure was constantly lowered until choking flow configurations 311 were obtained. As for the transonic bump test case before, the k− ω SST turbulence model of [36] was 312 used. p = 67760 Pa p = 58000 Pa 0 Max T = 376.6 K Figure 6. Moses-Stein Nozzle. The inlet total pressure was set to 67,760 Pa, while the outlet static pressure was varied. From an outlet pressure of 58,000 Pa the pressure was constantly lowered until choking flow configurations RMS Residual I nlet Pressure [Pa] Outlet Int. J. Turbomach. Propuls. Power 2020, 5, 20 18 of 32 were obtained. As for the transonic bump test case before, the k-w SST turbulence model of [36] was used. Initially, the computational domain is thermophysically assigned to region 2 of the IAPWS standard, i.e., gaseous state. When lowering the outlet static pressure, the temperature drop due to Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 19 of 32 the acceleration in the nozzle theoretically leads to cells that are assigned to the two-phase region of the IAPWS standard, i.e., region 4. However, it was observed that the commercial code software restricts to region 2 of the IAPWS standard, when running single fluid. Figure 7a shows the centerline temperature along the nozzle geometry from inlet to outlet and demonstrates this behaviour. The dashed line is the solution of the in-house code, the continuous line shows the solution of the commercial code. Extrapolation of region 2 into the two-phase region 4, as done by the commercial p = 67760 Pa p = 58000 Pa 0 Max T = 376.6 K code, leads to considerably lower temperatures in the throat area. The in-house code applies the proper formulations of the IAPWS standard if phase transition is detected. However, it needs to be mentioned that both solutions do not entirely reflect the full reality. Although the properties of the water are clearly better described using the state equations for the two phase region 4, it is not equal to a full two-phase Figure 6. Moses-Stein Nozzle simulation. The assumption of instant thermal equilibrium for example is an oversimplification of the actual problem and prevents an accurate modeling of the phase change. Being able to switch the 313 Initially, the computational domain is thermophysically assigned to region 2 of the IAPWS standard, i.e. state equation in a region where the flow properties change rapidly, however, is a demonstration of 314 gaseous state. When lowering the outlet static pressure, the temperature drop due to the acceleration in the stability of the underlying methods. An option was implemented that allows a restriction of the 315 the nozzle theoretically leads to cells that are assigned to the two-phase region of the IAPWS standard, in-house code to stay within a prescribed region. In case a transition into the two-phase region is 316 i.e. region 4. However it was observed that the commercial code software restricts to region 2 of detected, the values of region 2 are extrapolated. The effect is presented in Figure 7b. The solutions 317 the IAPWS standard, when running single fluid. Fig. 7a shows the centerline temperature along the obtained with the in-house framework match the commercial code. nozzle geometry from inlet to outlet and demonstrates this behaviour. The dashed line is the solution 340 IC IC CC CC Corrdinate [m] Corrdinate [m] (a) Free IAPWS Regions (b) Fixed to Region 2 Figure 7. Influence of the Figure presence 7. Influence of Region of4. the In-house presence code of Region (IC) vs. 4. commercial code (CC). In-house code (IC) vs. commercial code (CC) The following validation was therefore carried out restricting the framework to region 2 of the 319 of the in-house code, the continuous line shows the solution of the commercial code. Extrapolation IAPWS standard. Figure 8 shows the temperature and pressure distribution along the symmetry plane 320 of region 2 into the two-phase region 4, as done by the commercial code, leads to considerably lower for three different outlet pressures. As can be seen the results are in close agreement with commercial 321 temperatures in the throat area. The in-house code applies the proper formulations of the IAPWS code results, over the complete range from subsonic flow (black lines) to transonic flow (gray lines) 322 standard if phase transition is detected. However, it needs to be mentioned that both solutions do until supersonic, chocking flow regime (light gray lines). In Figure 8a a continuous line can be seen for 323 not entirely reflect the full reality. Although the properties of the water are clearly better described the 15,000 Pa case for both numerical algorithms. The pressure curve in Figure 8b does not show this 324 using the state equations for the two phase region 4, it is not equal to a full two-phase simulation. The sharp limit. This is because the value of the pressure in the interior of the domain is a direct outcome 325 assumption of instant thermal equilibrium for example is an oversimplification of the actual problem of the linearized equation system, while temperature is computed as a post-processing field based on 326 and prevents an accurate modeling of the phase change. Being able to switch the state equation in a the values of pressure and enthalpy. In this post-processing step the values for the temperature are 327 region where the flow properties change rapidly however is a demonstration of the stability of the limited to 273.25 K which is the limit of defined temperatures of the IAPWS-97 formulation. 328 underlying methods. An option was implemented that allows a restriction of the in-house code to stay 329 within a prescribed region. In case a transition into the two-phase region is detected, the values of 330 region 2 are extrapolated. The effect is presented in Fig. 7b. The solutions obtained with the in-house 331 framework match the commercial code. 332 The following validation was therefore carried out restricting the framework to region 2 of the IAPWS 333 standard. Fig. 8 shows the temperature and pressure distribution along the symmetry plane for three 334 different outlet pressures. As can be seen the results are in close agreement with commercial code 335 results, over the complete range from subsonic flow (black lines) to transonic flow (gray lines) until Temperature [K] I nlet −0.05 0.00 0.05 0.10 Temperature [K] −0.05 0.00 0.05 Outlet 0.10 Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 20 of 32 336 supersonic, chocking flow regime (light gray lines). In Fig. 8a a continuous line can be seen for the 337 15000Pa case for both numerical algorithms. The pressure curve in Fig. 8b does not show this sharp 338 limit. This is because the value of the pressure in the interior of the domain is a direct outcome of the 339 linearized equation system, while temperature is computed as a post-processing field based on the 340 Int. values J. Turbomach. of pressur Propuls. e and Power enthalpy 2020, 5, 20 . In this post-processing step the values for the temperature are19 limited of 32 to 273.25K which is the limit of defined temperatures of the IAPWS-97 formulation. ·10 IC 58000 Pa IC 58000 Pa CC 58000 Pa CC 58000 Pa IC 40000 Pa IC 40000 Pa CC 40000 Pa CC 40000 Pa IC 15000 Pa IC 15000 Pa CC 15000 Pa CC 15000 Pa Corrdinate [m] Corrdinate [m] (a) Temperature distribution (b) Pressure distribution Figure 8. Allowing region 2 only. In-house code (IC) vs. commercial code (CC). Figure 8. Allowing region 2 only. In-house code (IC) vs. commercial code (CC) 8.3. Two-Stage Radial Compressor Final validation is carried out on the basis of a two-stage radial compressor with a refrigerant as a 342 8.3. Two-stage radial compressor working fluid. Due to the governmental requirements on the global warming potential (GWP), 343 Final validation is carried out on the basis of a two-stage radial compressor with a refrigerant as a replacement of the refrigerant became necessary. This led to a redesign of the components 344 a working fluid. Due to the governmental requirements on the GWP (Global Warming Potential), a accompanied with measurements herein used for the validation of the implemented numerical model. 345 replacement of the refrigerant became necessary. This led to a redesign of the components accompanied The main objective of the measurement campaign was to gain information on the surge line under 346 with measurements herein used for the validation of the implemented numerical model. The main various operating conditions to allow for a controlled operation of the compressor. A vaneless radial 347 objective of the measurement campaign was to gain information on the surge line under various diffuser follows each impeller while the return channel contains a single row of stationary blades to 348 operating conditions to allow for a controlled operation of the compressor. A vaneless radial diffuser reduce the swirl before the second impeller. Figure 9 shows the individual parts. 349 follows each impeller while the return channel contains a single row of stationary blades to reduce the swirl before the second impeller. Figure 9 shows the individual parts. Figure 9. Parts overview. 8.3.1. Measurement Setup Figure 9. Parts Overview The two stage compressor was installed in the laboratory of the Institute for Thermal Energy System and Process Engineering at HSLU in Lucerne. Restriction in the laboratory environment, however, led to operating conditions , which are not in line with the design. This concerns mainly the inlet temperature and pressure levels, which could not be tuned to design conditions. In order to account for this inconsistency, dynamic scaling was applied to match the operating conditions of the design. The rotational speed of the compressor was corrected using Equation (83), with f being the design frequency and f the corrected frequency of the laboratory setup. The mass Temperature [K] −0.05 0.00 0.05 0.10 Pressure [Pa] −0.05 0.00 0.05 0.10 Int. J. Turbomach. Propuls. Power 2020, 5, 20 20 of 32 flow is then corrected based on the non-dimensional mass flow as given Equation (83) for comparison with the design conditions. m ˙ f p M 0 0 1.5 f = q m ˙ = z z = (83) T rT R z p 8.3.2. Numerical Setup The manufacturers choice for the real gas state equations was the standard Peng–Robinson model as given in Equation (2). In order to account for the temperature dependency of the dynamic viscosity, the interacting sphere model given in [39,40] was applied. Furthermore, the thermal conductivity is modeled using the modified Eucken formulation as described in [41]. This decision is mainly based on comparability considerations concerning the CFD computations which have been carried out using a commercial tool during the measurement campaign. Of the available models in the commercial tool, the Peng–Robinson EOS showed the smallest error for the chosen working fluid and operating conditions. Table 1 shows the relative error comparing the chosen thermophysical models against tabulated values from REFPROP [42] over the operating range of the working fluid. A significant error can be seen in laminar conductivity and viscosity; however, this effect is negligible since dominated by the turbulent contributions. The error in the heat capacity was tested on its effect on the results and a change in the value in the range given in the table did not have a significant influence on the final results. Table 1. Relative error for Peng–Robinson and related models. Pressure [bar] Temperature [ C] r c g a l m 2 3 0.48% 3.15% 0.82% 0.13% 10.08% 11.46% 3 12 0.53% 4.04% 1.06% 0.11% 9.43% 11.60% 4 18 0.58% 4.92% 1.33% 0.09% 8.96% 11.63% 5 27 0.46% 5.19% 1.37% 0.02% 8.14% 11.30% 6 33 0.40% 5.64% 1.51% 0.10% 7.49% 10.99% 7 39 0.30% 5.95% 1.58% 0.20% 6.81% 10.53% 8 42 0.31% 6.66% 1.84% 0.24% 6.19% 10.13% 9 48 0.14% 6.69% 1.81% 0.38% 5.42% 9.48% 10 51 0.11% 7.25% 2.01% 0.45% 4.68% 8.89% 11 54 0.06% 7.75% 2.18% 0.51% 3.86% 8.20% 12 60 0.19% 7.39% 1.96% 0.68% 3.10% 7.39% The ideal gas contribution of any internal energy change (Path B in Figure 2) is evaluated using a fourth order polynomial for the specific heat capacity at constant pressure. Turbulence is accounted for by the use of the Menter-SST model [36] with automatic wall treatment. The mesh consists of a total of nine different regions connected with fourteen interfaces for general grid connections and six in-house developed mixing planes, connecting the regions with different periodicities. Details on the implementation of the mixing plane can be found in [43]. The computational domain consists of 4 mio cells with y values around 20. 8.3.3. Global Performance Comparison The complete operating range for a given rotational speed was computed using the above described real gas models by continuously increasing the backpressure until the stability limit was reached. In order to assess the quality of the results, the solutions are compared against commercial Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 22 of 32 373 Global Performance Comparison 374 The complete operating range for a given rotational speed was computed using the above 375 described real gas models by continuously increasing the backpressure until the stability limit was 376 reached. In order to assess the quality of the results, the solutions are compared against commercial 377 code results and measurement data. At the time, the only measurement data available for comparison Int. J. Turbomach. Propuls. Power 2020, 5, 20 21 of 32 378 was the pressure ratio. This is due to the fact that the main objective of the measurement campaign 379 was to find the surge lines to allow for a controlled operation. To validate the implemented procedure, code results and measurement data. At the time, the only measurement data available for comparison 380 the results are further compared with commercial code results by plotting polytropic efficiency (Eq. was the pressure ratio. This is due to the fact that the main objective of the measurement campaign 381 84) and work coefficient (Eq. 85) against flow coefficient (Eq. 86). The circumferential velocity u is was to find the surge lines to allow for a controlled operation. To validate the implemented 382 evaluated at the radius (D ) of the first compressor outlet. procedure, the results are further compared with commercial code results by plotting polytropic efficiency Equation (84) and work coefficient Equation (85) against flow coefficient Equation (86). The circumferential velocity u is evaluated at the radius (D ) of the first compressor outlet. Δh (s − s )(T − T ) p0 2 1 0,2 0,1 η = Δh = (h − h )− (84) p p0 0,2 0,1 Dh T (s s )(T T ) h p0− h 0,2 2 1 0,2 0,1 0,2 0,1 ln h = Dh = (h h ) (84) p p0 0,2 0,1 0,1 0,2 h h 0,2 0,1 ln 0,1 Δh Λ = (85) Dh ( 0u ) L = (85) (u ) 2 ˙ Φ = (86) V 2 u (D ) 2 2 F = (86) u (D ) 2 2 384 Figure 10a shows the results for the whole operating range from choking flow to the stability Figure 10a shows the results for the whole operating range from choking flow to the stability 385 limit. The legend entries IC and CC stand for the in-house code and commercial code respectively. It is limit. The legend entries IC and CC stand for the in-house code and commercial code, respectively. 386 evident that the solutions are in very close agreement, indicating the consistency of the underlying It is evident that the solutions are in very close agreement, indicating the consistency of the underlying 387 thermophysical models. Comparison against measurement data was only available in terms of static thermophysical models. Comparison against measurement data was only available in terms of static 388 pressure ratio, shown in gray. The differences in the results between numerical predictions and pressure ratio, shown in gray. The differences in the results between numerical predictions and 389 measurement setup are dominated by the limitations of the laboratory setup as explained in the section measurement setup are dominated by the limitations of the laboratory setup as explained in the section 390 covering the measurement setup. covering the measurement setup. 1.4 5 −3 0.8 1.2 −4 −5 0.7 −6 0.8 −7 IC Ux 0.6 −8 IC PR CC Ux 0.6 CC PR IC En 2 −9 Measurement CC En −10 0.5 0.4 10 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 Inlet Volume Flow Coefficient Iteration (a) Performance Map Overall System (b) Convergence History Figure 10. Two-stage compressor results. In-house Code (IC) vs commercial code (CC) Figure 10. Two-stage compressor results. In-house Code (IC) vs commercial code (CC). 391 History 392 In order to provide a reliable tool for the design of turbomachinery applications operating in real In order to provide a reliable tool for the design of turbomachinery applications operating in real 393 gas thermophysical regions, accurate numerical solutions are not enough. To allow for a high work gas thermophysical regions, accurate numerical solutions are not enough. To allow for a high work 394 efficiency and reduced resources, numerical stability of the iterative process together with a small efficiency and reduced resources, numerical stability of the iterative process together with a small 395 computational overhead are equally important. Investigations considering the convergence behavior computational overhead are equally important. Investigations considering the convergence behavior 396 have shown promising results as can be seen from Fig 10b. have shown promising results as can be seen from Figure 10b. 397 In order to examine the effects of real gas modeling compared to perfect gas thermophysical models, In order to examine the effects of real gas modeling compared to perfect gas thermophysical 398 cpu time was compared as well. An increase of about 50% compared when switching to real gas models models, cpu time was compared as well. An increase of about 50% compared when switching to real 399 could be observed for this test case setup. A value well within what is observed from commercial code gas models could be observed for this test case setup. A value well within what is observed from 400 results. commercial code results. 8.3.4. Real Gas Effects In this section, we only compare solutions obtained with the in-house framework. The investigation focuses on differences in the results when using real gas models instead of perfect gas assumptions. In order to do so, the characteristic was recomputed using standard models for Polytropic Efficiency (dashed line) Work Coefficient (dotted line) Static Pressure Ratio (continous line) RMS Residual Int. J. Turbomach. Propuls. Power 2020, 5, 20 22 of 32 perfect gas as Newtonian viscosity, constant heat capacity and Prandtl’s law for the conductivity Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 23 of 32 based on average values from the real gas solutions. The results are shown in Figure 11. Lines with square markers correspond to solutions of the perfect gas model (PG) and lines with circular markers correspond to the real gas solution (RG) based on Peng–Robinson state equation as already shown in 401 Real Gas Effects Figure 10a. 402 In With this respect section to polytr we only opic ef compar ficiency e (dashed solutions line), obtained similar rwith esults the are obtained in-house for framework. flow coefficients The 403 investigation above the peak focuses efficiency on dif point. ferences Considering in the results peak when efficiency using , areal differ gas ence models of 2% instead can be of found. perfect While this is surely not negligible, the differences are too small to draw any general conclusions. 404 gas assumptions. In order to do so, the characteristic was recomputed using standard models for Comparing the curves of static pressure ratio (continuous line), it is however clear that when using 405 perfect gas as Newtonian viscosity, constant heat capacity and Prandtl’s law for the conductivity perfect gas state equation lower back pressure is necessary for the same volume flow coefficient. 406 based on average values from the real gas solutions. The results are shown in Fig. 11. Lines with Using the perfect gas model, it was not possible to reach the same pressure levels as when using the 407 square markers correspond to solutions of the perfect gas model (PG) and lines with circular markers real gas model. The compressor stability limit was detected at much smaller levels of back pressure 408 correspond to the real gas solution (RG) based on Peng-Robinson state equation as already shown in and the flow coefficient dropped much faster for the perfect gas when increasing the back pressure. Fig. 10a. 1.4 0.8 1.2 IP 1 IP 2 0.7 0.6 0.8 IC RG 0.5 0.6 IC PG Measurement 0.4 0.4 0.6 0.7 0.8 0.9 1 Inlet Volume Flow Coefficient Figure 11. Performance Map: Real Gas (RG) vs. Perfect Gas (PG) Figure 11. Performance map: real gas (RG) vs. perfect gas (PG). 410 With respect to polytropic efficiency (dashed line), similar results are obtained for flow coefficients Further investigation of the marked regions (IP1 in blue, IP2 in red) show the differences. 411 above the peak efficiency point. Considering peak efficiency, a difference of∼ 2% can be found. While IP1 corresponds to peak efficiency for both setups, IP2 are choking conditions. The operating points 412 this is surely not negligible, the differences are too small to draw any general conclusions. Comparing of the CFD simulation were set by fixing the value of the pressure at the outlet. As can be seen from 413 the curves of static pressure ratio (continuous line), it is however clear that when using perfect gas IP1 (peak efficiency), to obtain the same flow coefficient, the perfect gas model needs a much smaller 414 state equation lower back pressure is necessary for the same volume flow coefficient. Using the perfect static pressure ratio between inlet and outlet. In fact, the outlet pressure of the perfect gas simulation 415 gas model, it was not possible to reach the same pressure levels as when using the real gas model. The is roughly 13% lower than for real gas. This is in contrast to IP2 (choking conditions) where pressure 416 compressor stability limit was detected at much smaller levels of back pressure and the flow coefficient ratio and flow coefficient are in close agreement. This can, however, be explained by Figure 12c. 417 dropped much faster for the perfect gas when increasing the back pressure. Under choking conditions, the flow is much closer to ideal gas behavior (Compressibility z = 1), 418 Further investigation of the marked regions (IP1 in blue, IP2 in red) show the differences. IP1 thus the effect of real gas modeling is smaller. 419 corresponds to peak efficiency for both setups, IP2 are choking conditions. The operating points Comparing the exit flow angles from Figure 12a the general orientation between real gas and 420 of the CFD simulation were set by fixing the value of the pressure at the outlet. As can be seen from perfect gas in stage 1 do not show a significant difference. This is independent whether choking 421 IP1 (peak efficiency), to obtain the same flow coefficient, the perfect gas model needs a much smaller or peak efficiency is considered. Increasing deviation can be found when comparing the results in 422 static pressure ratio between inlet and outlet. In fact, the outlet pressure of the perfect gas simulation stage 2, see Figure 12b. An explanation can be found when plotting the average compressibility factor 423 is roughly 13% lower than for real gas. This is in contrast to IP2 (choking conditions) where pressure in streamwise direction given in Figure 12c, showing the increasing deviation from ideal gas state 424 ratio through and the flow compr coeffi essor cient . This are in deviation close agrfr eement. om perfect Thisgas can modeling howeverassumptions be explainedbecome by Fig. mor 12c.eUnder and 425 choking conditions, the flow is much closer to ideal gas behavior (Compressibility z=1), thus the effect more dominant with increasing back pressure (denser gas), explaining the increased difference in 426 of polytr real opic gas modeling efficiencyis for smaller the peak . efficiency point. 427 Comparing the exit flow angles from Fig. 12a the general orientation between real gas and perfect gas 428 in stage 1 do not show a significant difference. This is independent whether choking or peak efficiency 429 is considered. Increasing deviation can be found when comparing the results in stage 2, see Fig. 12b. 430 An explanation can be found when plotting the average compressibility factor in streamwise direction 431 given in Fig. 12c, showing the increasing deviation from ideal gas state through the compressor. This 432 deviation from perfect gas modeling assumptions become more and more dominant with increasing 433 back pressure (denser gas), explaining the increased difference in polytropic efficiency for the peak 434 efficiency point. Polytropic Efficiency (dashed line) Work Coefficient (dotted line) Static Pressure Ratio (continous line) Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 24 of 32 Int. J. Turbomach. Propuls. Power 2020, 5, 20 23 of 32 Version July 9, 2020 submitted to Int. J. Turbomach. Propuls. Power 24 of 32 0.8 0.8 0.95 0.8 0.6 0.8 0.95 0.6 0.9 0.6 0.6 0.4 0.4 0.9 IP2 Ideal 0.4 0.4 IP2 Real 0.85 IP2 Ideal 0.2 0.2 IP1 Ideal IP2 Real 0.85 IP2 0.2 0.2 IP1 Real IP1 Ideal IP1 IP2 0 0 IP1 Real 0.8 IP1 0 20 40 0 20 40 60 80 0 0.5 1 1.5 2 0 0 0.8 ◦ ◦ 0 20 40 0 20 40 60 80 Flow Angle [ ] Flow Angle [ 0 0.5 1 1.5 2 Streamwise location Flow Angle [ ] Streamwise location (a) Stage 1 (b) Stage 2 (c) Compressibility (a) Stage 1 (c) Compressibility (b) Stage 2 Figure 12. Flow angle and compressibility Figure 12. Flow angle and compressibility. A normalized dynamic pressure was computed in order to compare the dynamic contribution A normalized dynamic pressure was computed in order to compare the dynamic contribution A normalized dynamic pressure was computed in order to compare the dynamic contribution between real and perfect gas computations as given in Eq. 87. between real and perfect gas computations as given in Eq. 87. between real and perfect gas computations as given in Equation (87). p p −−p ¯ p ¯ p p Ψ = (87) Ψ = (87) Y 0 = (87) p ¯− p ¯ p ¯− p ¯ p ¯0 p ¯ 435 Fig. 13 shows the results for IP1 (Peak efficiency) at sections close to the impeller outlets. For stage 1 Figure 13 shows the results for IP1 (Peak efficiency) at sections close to the impeller outlets. 435 Fig. 13 shows the results for IP1 (Peak efficiency) at sections close to the impeller outlets. For stage 1 436 the results are comparable, differences however become present in the second stage. This is especially For stage 1 the results are comparable, differences however become present in the second stage. This is 436 the results are comparable, differences however become present in the second stage. This is especially 437 true for the suction side (SS) of the splitter blade (SB), where the dynamic contribution is much lower. especially true for the suction side (SS) of the splitter blade (SB), where the dynamic contribution is 437 true for the suction side (SS) of the splitter blade (SB), where the dynamic contribution is much lower. much lower. (a) Real Gas, Stage 1 (b) Ideal Gas, Stage 1 (a) Real Gas, Stage 1 (b) Ideal Gas, Stage 1 (c) Real Gas, Stage 2 (d) Ideal Gas, Stage 2 Figure 13. Normalized Total Pressure Isolines, IP1 (c) Real Gas, Stage 2 (d) Ideal Gas, Stage 2 439 These differences in terms of compressor dynamics when switching between real and perfect Figure 13. Normalized total pressure isolines, IP1. 440 gas model explain the importance of an accurate representation of the thermophysical properties and 439 These differences in terms of compressor dynamics when switching between real and perfect 441 derivation of the conservation equations. Using perfect gas models would lead to errors in the design These differences in terms of compressor dynamics when switching between real and perfect 440 gas model explain the importance of an accurate representation of the thermophysical properties and 442 of blade geometries because the flow physics are not captured accurately. gas model explain the importance of an accurate representation of the thermophysical properties and 441 derivation of the conservation equations. Using perfect gas models would lead to errors in the design derivation of the conservation equations. Using perfect gas models would lead to errors in the design 443 CONCLUSIONS 442 of blade geometries because the flow physics are not captured accurately. of blade geometries because the flow physics are not captured accurately. 444 The present work summarizes the necessary steps during the derivation and implementation of 9. Conclusions 443 CONCLUSIONS 445 Navier-Stokes equations into a pressure based numerical framework when operating in non-ideal gs 446 conditions. Detailed analysis is given concerning the discretization of the conservation equations The present work summarizes the necessary steps during the derivation and implementation of 444 The present work summarizes the necessary steps during the derivation and implementation of 447 when Navier using –Stokes realequations gas state equations. into a pressur This e based is extended numerical furthermor framework e to the when development operating of in boundary non-ideal 445 Navier-Stokes equations into a pressure based numerical framework when operating in non-ideal gs 448 conditions and iterative procedures for the update of various flow quantities. A transonic bump gs conditions. Detailed analysis is given concerning the discretization of the conservation equations 446 conditions. Detailed analysis is given concerning the discretization of the conservation equations 449 case when was using chosen real to gas demonstrate state equations. the consistency This is extended of thefurthermor modifications. e to the The development added iterative of boundary solution 447 when using real gas state equations. This is extended furthermore to the development of boundary 450 pr conditions ocedures for and the iterative thermophysical procedur states es forreduce the update to the of same various results flow as the quantities. direct algebraic A transonic relations bump for 448 conditions and iterative procedures for the update of various flow quantities. A transonic bump 449 case was chosen to demonstrate the consistency of the modifications. The added iterative solution 450 procedures for the thermophysical states reduce to the same results as the direct algebraic relations for Span Span Compressibility Factor Compressibility Factor Inlet Inlet Impeller 1 Impeller 1 Diffuser 1 Crossover Diffuser 1 Bend Crossover Impeller 2 Bend Diffuser Impeller 2 2 Outlet Diffuser 2 Outlet Int. J. Turbomach. Propuls. Power 2020, 5, 20 24 of 32 case was chosen to demonstrate the consistency of the modifications. The added iterative solution procedures for the thermophysical states reduce to the same results as the direct algebraic relations for perfect gas. The real gas models were then tested using two different state equations, once using the IAPWS-97 formulation for water and once using the Peng–Robinson cubic state equation. Comparison against the commercial code results was carried out using the IAPWS-97 formulation on a transonic nozzle configuration. It was found that commercial code extrapolates the gaseous equations into the two-phase region. Applying the same treatment to the in-house code then allowed further comparison and a close agreement between the two frameworks was found. Cubic state equations were then applied to an industrial type compressor setup and compared against commercial code results and measurement data. It could be proven that the implemented procedure is not only able to accurately replicate commercial code results but also leads to a stability of the iterative scheme that goes beyond the capabilities of the used commercial tool. Comparison against perfect gas assumptions justify the increased complexity of the implementation. It could be demonstrated, based on overall performance prediction and local flow features, that the assumption of perfect gas shows its limitation under certain operating conditions, especially away from choking conditions. At moderate backpressure levels, however, i.e., high volume flow rate, the assumption of perfect gas might be a viable option. Especially during preliminary design phase, the possibility to use perfect gas assumption without great loss of accuracy can lower the computational costs and therefore increase the throughput. For a thorough analysis of the geometry, real gas state equations are, however, a necessity. Author Contributions: The extension of a pressure based coupled numerical algorithm to real gas flow physics was mainly driven by L.H. as it is a part of his PhD. The final target of the thesis is a framework, able at predicting two-phase flow physics for steam turbines. It is however based on the excessive work done by L.M. and M.D. in the development of such a coupled framework and thus also supported by these two authors. E.C. and D.M.V. provided valuable background information considering thermophysical relations and the various model assumptions when moving to real gas flow physics and their effect on turbomachinery design. All authors contributed to writing and editing this paper. All authors have read and agreed to the published version of the manuscript. Funding: The authors gratefully acknowledge the financial contribution provided by the Swiss National Science Foundation (SNF: Grant-Number 175900). Conflicts of Interest: The authors declare no conflict of interest. Nomenclature Abbreviations A 4  4 block matrix a Matrix element b Source term CC Commercial code c Speed of sound [m/s] c Specific heat at constant pressure [J/kgK] c Specific heat at constant volume [J/kgK] d Face to cell distance vector [m] D Inverted 4  4 block matrix e Internal energy [J/kg] F Body force [kg/ms ] bi f Frequency [1/s] h Enthalpy [J/kg] IC In-house code I P Investigation point Ma Mach number Int. J. Turbomach. Propuls. Power 2020, 5, 20 25 of 32 M Molar mass [kg/mol] m Mass flow [kg/s] MR Maxwell relation PG Perfect Gas p Pressure [Pa] q Heat [J/kg] R Universal gas constant [J/kgK] R Specific gas constant [J/kgK] RG Real Gas S Surface vector [m ] f i s Entropy [J/kgK] T Temperature [K] t Time [s] U Velocity flux [m /s] u Velocity x-component [m/s] v Velocity y-component [m/s] w Velocity z-component [m/s] V Volume [m ] v Specific volume [1/kg] z Compressibility Greeks a Thermal expansivity [1/K] b Compressibility [1/Pa] G Fundamental derivative g Heat capacity ratio h Polytropic efficiency L Work coefficient l Thermal conductivity [W/mK] m Dynamic viscosity [kg/ms] Y Normalized dynamic pressure [-] r Density [kg/m ] 2 2 t Stress tensor [kg/m s ] i j F Flow coefficient w Accentric factor Subscripts b Boundary value c Cell value cr Critical point values f Surface values ¥ Ideal gas state nb Neighbour matrix element p Diagonal matrix element r Reference state 0 Stagnation conditions Previous iteration value Superscripts n Time step index u Coupling to velocity x-component v Coupling to velocity y-component w Coupling to velocity z-component p Coupling to pressure Correction to given quantity Int. J. Turbomach. Propuls. Power 2020, 5, 20 26 of 32 Appendix A. Entropy Change Using the first and second law of thermodynamics, the change in specific internal energy can be described as given in Equation (A1). de = Tds pdv (A1) Using the fundamental derivative, ds can be rewritten as given in Equation (A2). ¶s ¶s ds = dT + dv (A2) ¶T ¶v v T Inserting Equation (A2) into Equation (A1) and using Maxwell relations, the general form given in Equation (A3) can be found. ¶ p de = c dT + T p dv (A3) ¶T Appendix B. Partial Thermal Derivatives Table A1. Constant Pressure. Derivative General ¶r ¶h vc ¶r ¶T v ¶T 1 ¶h c Table A2. Constant Volume. Derivative General ¶h c b + (1 Ta)v ¶ p ¶s ¶T ¶T b ¶ p a ¶h 2 c vT ¶T Table A3. Constant Temperature. Derivative General ¶v ¶s ¶ p bv ¶v ¶h aT1 ¶v T Int. J. Turbomach. Propuls. Power 2020, 5, 20 27 of 32 Table A4. Constant Enthalpy. Derivative General ¶r 2 Ta a rb ¶ p Table A5. Constant Entropy. Derivative General ¶ p c v p b c ¶r Appendix C. Real Gas Internal Energy Using Peng–Robinson State Equation As stated earlier in Equation (6) the internal energy of any real gas can be expressed as given in Equation (A4). Z Z Z v T v ¶ p ¶ p e(T, p) = T p dv + [c ] dT + T p dv +e(T , p ) (A4) v r r ¶T ¶T v T v r r ¥ v v T T | {z } | {z } | {z } ¶ p Using Peng–Robinson [12] state equation the derivative can be expressed as given in ¶T Equation (A7). ¶ p ¶ R T aa(T) = (A5) 2 2 ¶T ¶T v b v + 2bv b v v ¶a(T) a2a R ¶T = (A6) 2 2 v b v + 2bv b k T aa R T T s c = + (A7) 2 2 v b v + 2bv b ¶ p Hence the term T p can be written as given in Equation (A8). ¶T ¶ p aa(T)(1 + k) T p = (A8) 2 2 ¶T v + 2bv b Inserting Equation (A8) into term A of Equation (A4) leads to Equation (A9). " !# ¶ p aa(T )(1 + k) 2b T p dv = p tanh (A9) ¶T b + v r 2b This allows to rewrite the internal energy for any state T, v as given in Equation (A12) aa(T )(1 + k) 2b e(T, v) = p tanh (A10) b + v 2b r aa(T)(1 + k) 2b p tanh (A11) b + v 2b + c dT R (T T ) + e(T , v ) (A12) p s r r r r Int. J. Turbomach. Propuls. Power 2020, 5, 20 28 of 32 Appendix D. Real Gas Internal Energy Using Benedict–Webb–Rubin (BWR) State Equation As shown in Appendix C, the expression needed for the real gas contribution of the state change ¶ p is the term T p. BWR [15] state equation is given as: ¶T C c r 2 0,s s 2 3 6 2 g r p = rR T + B R T A r + (b R T a ) r + a a r + 1 + g r e (A13) s s s s s s s s 0,s 0,s 2 2 T T ¶ p Hence can be written is: ¶T ¶ p C c r 2 0,s 2 3 2 g r = rR + (B R + 2 )r + b R r 2 (1 + g r )e (A14) s 0,s s s s s 3 3 ¶T T T Hence the requested term reads as: ¶ p C c r 2 0,s 2 2 2 6 2 g r T p = r + A r + ar a a r 1 + g r e (A15) 0,s s s s 2 2 ¶T T T The remaining procedure just follows what is shown in Appendix C. Appendix E. Computation of Entropy Using the first law of thermodynamics and Maxwell relations, the entropy can be expressed as shown in Equation (A16). ¶ p c ds(v, T) = dv + dT (A16) ¶T T Integration along paths according to Figure 2 leads to Equation (A17). Z Z Z v T v ¶ p c ¶ p s(T, v) s(T , v ) = dv + dT R ln(T/T ) + dv (A17) r r s r T T ¶T T ¶T v T v r r ¥ v v Peng–Robinson This is an example of the implemented formulation when using Peng–Robinson [12] state equation. R a (T) s(T, v) s(T , v ) = dv (A18) r r T 2 2 v b v + 2bv b R a (T) + dv (A19) 2 2 v b v + 2bv b + dT R ln(T/T ) (A20) s r " !# a (T ) 2b s(T, v) s(T , v ) = R ln(v b) p tanh (A21) r r s b + v 2b " !# a (T) 2b + R ln(v b) p tanh (A22) b + v 2b + dT R ln(T/T ) (A23) s r r Int. J. Turbomach. Propuls. Power 2020, 5, 20 29 of 32 a (T) 2b s(T, v) s(T , v ) = R ln(v b) p tanh (A24) r r s b + v 2b a (T ) 2b R ln(v b) + p tanh (A25) s r b + v 2b r + dT R ln(T/T ) (A26) s r Appendix F. Specific Heat Capacity The specific heat at constant volume is defined as given in Equation (A27). ¶e c = (A27) ¶T Using the definition of internal energy of Equation (6) this equation can be rewritten as given in Equation (A30) ¶ ¶ p c = c T p dv (A28) v v ¶T ¶T ¶ ¶ p ¶ p = c T dv (A29) ¶T ¶T ¶T v T ¶ p = c T dv (A30) ¶T v T ¶ p The partial derivative can be obtained through differentiation of the state equation. ¶T Values for the ideal gas contribution (c ) can be found through various models, e.g., constant value or NASA’s polynomial expression [18]. To obtain an expression for the specific heat at constant pressure, the general statement of Equation (A32) is used. ¶q ¶e c c = (A31) p v ¶T ¶T p v ¶s ¶s = T (A32) ¶T ¶T p v Using the substantial derivative of the entropy and specific volume ¶s ¶s ds = dT + dv (A33) ¶T ¶v v T ¶v ¶v dv = dT + d p (A34) ¶T ¶ p p T Substituting Equation (A34) in Equation (A33) results in Equation (A35). When expressing the change in entropy through fundamental derivative a second expression is found as given in Equation (A36). " # ¶s ¶s ¶v ¶s ¶v ds = + dT + d p (A35) ¶T ¶v ¶T ¶v ¶ p v T p T ¶s ¶s ds = dT + d p (A36) ¶T ¶ p p T Int. J. Turbomach. Propuls. Power 2020, 5, 20 30 of 32 With Equations (A35) and (A36) it is possible to express the derivative of entropy with respect to temperature at constant pressure as given in Equation (A37). ¶s ¶s ¶s ¶v = + (A37) ¶T ¶T ¶v ¶T p v T p Substituting Equation (A37) in Equation (A32) allows to rewrite c as given in Equation (A39). ¶s ¶v c c = T (A38) p v ¶v ¶T T p ¶ p ¶v c = c + T Maxwell relation (A39) p v ¶T ¶T v p For the density based state equations it is, however, advantageous to rewrite those formulations using triple product rule. Details on the partial derivatives can be found in Appendix B. ¶v ¶v ¶ p = Triple product rule (A40) ¶T ¶ p ¶T p T v ¶ p ¶T c = c T  (A41) p v ¶ p ¶v Appendix G. Speed of Sound ¶ p ¶ p = v (A42) ¶r ¶v ¶ p ¶s = v Triple product rule (A43) ¶s ¶v v p ¶T ¶ p = v Maxwell relation (A44) ¶v ¶T s s ¶T ¶s ¶ p ¶s = v Triple product rule (A45) ¶s ¶v ¶s ¶T v T T p c ¶ p ¶T ¶s = v c = T (A46) c ¶T ¶v ¶T v p v p ¶ p = v Triple product rule (A47) c ¶v References 1. Starzmann, J.; Casey, M.M.; Mayer, J.F.; Sieverding, F. Wetness loss prediction for a low pressure steam turbine using computational fluid dynamics. Proc. Inst. Mech. Eng. Part A J. Power Energy 2014, 228, 216–231. 2. Colonna, P.; Rebay, S.; Harinck, J.; Guardone, A. Real-gas effects in ORC turbine flow simulations: Influence of thermodynamic models on flow fields and performance parameters. In Proceedings of the European Community on Computational Methods in Applied Sciences (ECCOMAS) CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, Egmond aan Zee, The Netherlands, 5–8 September 2006; Delft University of Technology: Delft, The Netherlands, 2006. 3. Pohl, S.; Jarczyk, M.; Pfitzner, M.; Rogg, B. Real gas CFD simulations of hydrogen/oxygen supercritical combustion. Prog. Propuls. Phys. 2013, 4, 583–614. Int. J. Turbomach. Propuls. Power 2020, 5, 20 31 of 32 4. Odabaee, M.; Sauret, E.; Hooman, K. CFD Simulation of a Supercritical Carbon Dioxide Radial-Inflow Turbine, Comparing the Results of Using Real Gas Equation of Estate and Real Gas Property File. Appl. Mech. Mater. 2016, 846, 85–90. 5. Boncinelli, P.; Rubechini, F.; Arnone, A.; Cecconi, M.; Cortese, C. Real Gas Effects in Turbomachinery Flows: A CFD Model for Fast Computations. In Proceedings of the ASME Turbo Expo 2003, Collocated with the 2003 International Joint Power Generation Conference, Atlanta, GA, USA, 16–19 June 2003; American Society of Mechanical Engineers: New York, NY, USA, 2003; pp. 1103–1112. 6. Travis, J.; Koch, D.P.; Xiao, J.; Xu, Z. Real-gas equations-of-state for the GASFLOW CFD code. Int. J. Hydrogen Energy 2013, 38, 8132–8140. 7. Lucas, C.; Rusche, H.; Schroeder, A.; Koehler, J. Numerical investigation of a two-phase CO ejector. Int. J. Refrig. 2014, 43, 154–166. 8. Ameli, A.; Afzalifar, A.; Turunen-Saaresti, T.; Backman, J. Effects of real gas model accuracy and operating conditions on supercritical CO compressor performance and flow field. J. Eng. Gas Turbines Power 2018, 140, 062603. 9. Mangani, L.; Buchmayr, M.; Darwish, M. Development of a Novel Fully Coupled Solver in OpenFOAM: Steady State Incompressible Turbulent Flows. Numer. Heat Transf. Part B Fundam. 2014, 66, 1–20. doi:10.1080/10407790.2014.894448. 10. Mangani, L.; Buchmayr, M.; Darwish, M. Development of a Novel Fully Coupled Solver in OpenFOAM: Steady State Incompressible Turbulent Flows in Rotational Reference Frames. Numer. Heat Transf. Part B Fundam. 2014, 66, 526–543. doi:10.1080/10407790.2014.894371. 11. Mangani, L.; Moukalled, F.; Darwish, M. An OpenFOAM pressure-based coupled CFD solver for turbulent and compressible flows in turbomachinery applications. Numer. Heat Transf. Part B Fundam. 2016, 69 413–431. 12. Peng, D.Y.; Robinson, D.B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. 13. Cooper, J.; Dooley, R. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam; IAPWS R7-97(2012); The International Association for the Properties of Water and Steam: Lucerne, Switzerland, 2007. 14. Young, J.B. An Equation of State for Steam for Turbomachinery and Other Flow Calculations. J. Eng. Gas Turbines Power 1988, 110, 1–7, doi:10.1115/1.3240080. 15. Benedict, M.; Webb, G.B.; Rubin, L.C. An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures I. Methane, ethane, propane and n-butane. J. Chem. Phys. 1940, 8, 334–345. 16. Van der Waals, J.D. Over de Continuiteit van den Gas-en Vloeistoftoestand. Ph.D. Thesis, Leiden University: Leiden, The Netherlands, 1873. 17. Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 1972, 27, 1197–1203. 18. McBride, B.J.; Zehe, M.J.; Gordon, S. NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species; E-13336; National Aeronautics and Space Administration: Washington, DC, USA, September 2002. 19. Baltadjiev, N.D. An Investigation of Real Gas Effects in Supercritical CO Compressors. Ph.D. Thesis, Massachusetts Institute of Technology: Cambridge, MA, USA, 2012. 20. Baltadjiev, N.D.; Lettieri, C.; Spakovszky, Z.S. An investigation of real gas effects in supercritical CO centrifugal compressors. J. Turbomach. 2015, 137, 091003. 21. Barton, I.E. Comparison of SIMPLE-and PISO-type algorithms for transient flows. Int. J. Numer. Methods Fluids 1998, 26, 459–483. 22. Anjorin, V.; Barton, I. Removal of temporal and under-relaxation terms from the pressure-correction equation of the SIMPLE algorithm. Int. J. Fluid Dyn. 2001, 5, 59–75. 23. Kissling, K.; Springer, J.; Jasak, H.; Schutz, S.; Urban, K.; Piesche, M. A coupled pressure based solution algorithm based on the volume-of-fluid approach for two or more immiscible fluids. In Proceedings of the V European Conference on Computational Fluid Dynamics ECCOMAS CFD, Lisbon, Portugal, 14–17 June 2010. 24. Cerne, G.; Tiselj, I.; Petelin, S. Upgrade of the VOF Method for the Simulation of the Dispersed Flow. In Proceedings of the ASME 2000 Fluids Engineering Division Summer Meeting, Boston, MA, USA, 11–15 June 2000. Int. J. Turbomach. Propuls. Power 2020, 5, 20 32 of 32 25. Hutchinson, B.; Raithby, G. A multigrid method based on the additive correction strategy. Numer. Heat Transf. Part A Appl. 1986, 9, 511–537. 26. Keller, S.; Cordazzo, J.; Hinckel, P.H.; Maliska, C.R. Additive correction multigrid method applied to diffusion problems with unstructured grids. In Proceedings of the 10th Brazilian Congress of Thermal Sciences and Engineering, Rio de Janeiro, RJ, Brazil, 29 November–3 December 2004. 27. Darwish, M.; Sraj, I.; Moukalled, F. A coupled incompressible flow solver on structured grids. Numer. Heat Transf. Part B Fundam. 2007, 52, 353–371. 28. Moukalled, F.; Aziz, A.A.; Darwish, M. Performance comparison of the NWF and DC methods for implementing high-resolution schemes in a fully coupled incompressible flow solver. Appl. Math. Comput. 2011, 217, 5041–5054. 29. Darwish, M.; Moukalled, F. A fully coupled Navier-Stokes solver for fluid flow at all speeds. Numer. Heat Transf. Part B Fundam. 2014, 65, 410–444. 30. Darwish, M.; Abdel Aziz, A.; Moukalled, F. A coupled pressure-based finite-volume solver for incompressible two-phase flow. Numer. Heat Transf. Part B Fundam. 2015, 67, 47–74. 31. Mangani, L.; Buchmayr, M.; Darwish, M.; Moukalled, F. A fully coupled OpenFOAM® solver for transient incompressible turbulent flows in ALE formulation. Numer. Heat Transf. Part B Fundam. 2017, 71, 313–326. 32. Darwish, M.; Mangani, L.; Moukalled, F. Implicit boundary conditions for coupled solvers. Comput. Fluids 2018, 168, 54–66. 33. Rhie, C.; Chow, W.L. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J. 1983, 21, 1525–1532. 34. Darwish, M.; Sraj, I.; Moukalled, F. A coupled finite volume solver for the solution of incompressible flows on unstructured grids. J. Comput. Phys. 2009, 228, 180–201. 35. Demirdžic, ´ I.; Lilek, Ž.; Peric, ´ M. A collocated finite volume method for predicting flows at all speeds. Int. J. Numer. Methods Fluids 1993, 16, 1029–1050. 36. Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994, 32, 1598–1605. 37. Moses, C.A.; Stein, G.D. On the growth of steam droplets formed in a Laval nozzle using both static pressure and light scattering measurements. J. Fluids Eng. 1978, 100, 311–322. 38. Starzmann, J.; Hughes, F.R.; White, A.J.; Grübel, M.; Vogt, D.M. Numerical investigation of boundary layers in wet steam nozzles. J. Eng. Gas Turbines Power 2017, 139, 012606. 39. Chung, T.H.; Lee, L.L.; Starling, K.E. Applications of kinetic gas theories and multiparameter correlation for prediction of dilute gas viscosity and thermal conductivity. Ind. Eng. Chem. Fundam. 1984, 23, 8–13. 40. Chung, T.H.; Ajlan, M.; Lee, L.L.; Starling, K.E. Generalized multiparameter correlation for nonpolar and polar fluid transport properties. Ind. Eng. Chem. Res. 1988, 27, 671–679. 41. Reid, R.C.; Prausnitz, J.M.; Poling, B.E. The Properties of Gases and Liquids; McGraw Hill Book Co.: New York, NY, USA, 1987. 42. Lemmon, E.W.; .; Bell, I.H.; Huber, M.L.; McLinden, M.O. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0; National Institute od Standards and Technology: Gaithersburg, MD, USA, 2018. doi:10.18434/T4JS3C. 43. Hanimann, L.; Mangani, L.; Casartelli, E.; Mokulys, T.; Mauri, S. Development of a novel mixing plane interface using a fully implicit averaging for stage analysis. J. Turbomach. 2014, 136, 081010. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

Journal

International Journal of Turbomachinery, Propulsion and PowerMultidisciplinary Digital Publishing Institute

Published: Aug 3, 2020

There are no references for this article.