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

Lucian Hanimann; Luca Mangani; Ernesto Casartelli; Damian M. Vogt; Marwan Darwish

doi:10.3390/ijtpp5030020

Hanimann, Lucian;Mangani, Luca;Casartelli, Ernesto;Vogt, Damian M.;Darwish, Marwan

2020-08-03 00:00:00

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 ﬂow 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 signiﬁcant 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 ﬂow 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 ﬁrst comparing against results obtained with a framework using perfect gas simpliﬁcations. 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 ﬂows 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 modiﬁcations 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 speciﬁc 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 ﬁrst 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 conﬁguration with the operating range varying from supersonic to transonic ﬂow 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 ﬂow 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 ﬂows. It is a measurement for the change in speed of sound through a change in density, its deﬁnition 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) ﬂuids, 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 ﬂow 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 ﬁnal 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 exempliﬁed 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 ﬂuids, each having its own justiﬁcation depending on the operting conditions and the working ﬂuid. 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 ﬂuids or refrigerants, well known virial equations are those of Benedict–Webb–Rubin (BWR) [15] and its modiﬁcations. 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. Modiﬁcations 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 ﬁnal 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 ﬁve 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 Deﬁnition 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 deﬁnitions 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 brieﬂy described here on the basis of internal energy. The deﬁnition for the change in internal energy is given in Appendix A and the ﬁnal 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 speciﬁc volume (B) using a given reference state e(T , p ). The ﬁnal 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 ﬂuid. This is always valid if the molecules are separated by an inﬁnite distance. In other words the density goes towards zero or speciﬁc volume to inﬁnity. 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 ﬁnal 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) ﬂows 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 simpliﬁcation 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 justiﬁed 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. Speciﬁc 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 ﬁnal 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 deﬁnition 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 ﬁnal 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 ﬁnal 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 ﬂow chart will also deﬁne 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 ﬁnalized showing the chosen strategies for the update of ﬂow 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 ﬂux modeled through Fourier ’s approach as: ¶T q = l (18) ¶x 5. Numerical Discretization of Governing Equations The governing equations for compressible ﬂuid 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 coefﬁcient 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 fulﬁll the conservation of mass. Equation (20) shows the continuity equation in its semi discrete form with U = u S being the velocity ﬂux 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 exempliﬁed 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 ﬁelds 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 ﬂux 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 ﬂux on the face, computed using adjacent cell values through linear interpolation, contributing to the coefﬁcients a , a and a in the coefﬁcient matrix A . This couples the continuity pu pv pw p equation to the velocity ﬁeld, 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 coefﬁcient 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 ﬁrst 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 ﬁnal 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 ﬁrst 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 coefﬁcient 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 ﬁnal 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 deﬁnition 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 ﬂows 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 ﬂux. Using Fourier ’s law to relate the heat ﬂux to the local temperature gradient with l being the conductivity of the working ﬂuid, 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 simpliﬁcation. To recover at least partly the implicit contribution a ﬁrst 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 deﬁnition of the speciﬁc 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 ﬁrst 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 simpliﬁcations 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 ﬁnd an iterative solution procedure to accurately compute these quantities. The algorithm used is a Newton–Raphson root ﬁnding 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 deﬁnition 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 ﬁnding method is explained as a reference. This procedure is later applied to a variety of solution strategies used for ﬁnding iterative solutions for thermophysical changes in real gas. The procedure starts with a ﬁrst 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 deﬁned 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 coefﬁcient 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 ﬁrst order linearization of the pressure with respect to the velocity ﬂux U is introduced as given in Equation (63). b indexes values evaluated at boundary faces. Again, the velocity ﬂux is deﬁned 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 ﬁeld 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 ﬂux of the new ¶U time step U . The value for the velocity ﬂux 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 ﬁnal deﬁnition 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 ﬂux . ¶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 ﬂuid particle when brought to rest isentropically are not explicitly available for the general case of a real gas. ¶ p To approximate the gradient the deﬁnition 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 deﬁnition 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 ﬂow 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 ﬁnding procedure can be done in analogy to Section 6.1.1. 7.2. Total Pressure and Total Temperature In contrast to the problem of speciﬁed total conditions for boundary conditions as described in Section 6.1, ﬁnding 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 ﬁnding method is used. Given the static values of temperature and speciﬁc 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 deﬁned, 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 conﬁguration 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 ﬂuid 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 ﬂuid 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 conﬁguration published by Moses and Stein [37]. This The second test case is a nozzle conﬁguration published by Moses and Stein [37]. 305 conﬁguration is frequently used to investigate droplet growth and real gas ﬂow physics, see e.g. This conﬁguration is frequently used to investigate droplet growth and real gas ﬂow physics, see 306 [38]. Since this investigation of real gas ﬂow physics lays the basis for further extension to multiphase e.g., [38]. Since this investigation of real gas ﬂow physics lays the basis for further extension 307 ﬂows, this test case provides an ideal reference for testing the robustness of the implementation. Fig. 6 to multiphase ﬂows, 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 ﬂow conﬁgurations 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 ﬂow conﬁgurations 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 ﬂuid. 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 reﬂect 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 oversimpliﬁcation 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 ﬂow 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 ﬂuid. 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. Inﬂuence of the Figure presence 7. Inﬂuence 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 ﬂow (black lines) to transonic ﬂow (gray lines) 322 standard if phase transition is detected. However, it needs to be mentioned that both solutions do until supersonic, chocking ﬂow regime (light gray lines). In Figure 8a a continuous line can be seen for 323 not entirely reﬂect 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 oversimpliﬁcation of the actual problem of the linearized equation system, while temperature is computed as a post-processing ﬁeld 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 ﬂow properties change rapidly however is a demonstration of the stability of the limited to 273.25 K which is the limit of deﬁned 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 ﬂow (black lines) to transonic ﬂow (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 ﬂow 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 ﬁeld 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 deﬁned 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 ﬂuid. 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 ﬂuid. 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 ﬂow is then corrected based on the non-dimensional mass ﬂow 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 modiﬁed 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 ﬂuid 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 ﬂuid. A signiﬁcant 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 signiﬁcant inﬂuence on the ﬁnal 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 speciﬁc 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 ﬁnd 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 efﬁciency (Eq. was the pressure ratio. This is due to the fact that the main objective of the measurement campaign 381 84) and work coefﬁcient (Eq. 85) against ﬂow coefﬁcient (Eq. 86). The circumferential velocity u is was to ﬁnd the surge lines to allow for a controlled operation. To validate the implemented 382 evaluated at the radius (D ) of the ﬁrst compressor outlet. procedure, the results are further compared with commercial code results by plotting polytropic efﬁciency Equation (84) and work coefﬁcient Equation (85) against ﬂow coefﬁcient Equation (86). The circumferential velocity u is evaluated at the radius (D ) of the ﬁrst 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 ﬂow to the stability Figure 10a shows the results for the whole operating range from choking ﬂow 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 Coefﬁcient 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 efﬁciency and reduced resources, numerical stability of the iterative process together with a small efﬁciency 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 Efﬁciency (dashed line) Work Coefﬁcient (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 ﬁciency e (dashed solutions line), obtained similar rwith esults the are obtained in-house for framework. ﬂow coefﬁcients The 403 investigation above the peak focuses efﬁciency on dif point. ferences Considering in the results peak when efﬁciency 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 ﬂow coefﬁcient. 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 ﬂow coefﬁcient 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 Coefﬁcient 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 efﬁciency (dashed line), similar results are obtained for ﬂow coefﬁcients Further investigation of the marked regions (IP1 in blue, IP2 in red) show the differences. 411 above the peak efﬁciency point. Considering peak efﬁciency, a difference of∼ 2% can be found. While IP1 corresponds to peak efﬁciency 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 ﬁxing 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 efﬁciency), to obtain the same ﬂow coefﬁcient, the perfect gas model needs a much smaller 414 state equation lower back pressure is necessary for the same volume ﬂow coefﬁcient. 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 ﬂow coefﬁcient ratio and ﬂow coefﬁcient 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 ﬂow 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 efﬁciency for both setups, IP2 are choking conditions. The operating points Comparing the exit ﬂow angles from Figure 12a the general orientation between real gas and 420 of the CFD simulation were set by ﬁxing the value of the pressure at the outlet. As can be seen from perfect gas in stage 1 do not show a signiﬁcant difference. This is independent whether choking 421 IP1 (peak efﬁciency), to obtain the same ﬂow coefﬁcient, the perfect gas model needs a much smaller or peak efﬁciency 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 ﬂow compr coefﬁ 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 ﬂow 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 efﬁciencyis for smaller the peak . efﬁciency point. 427 Comparing the exit ﬂow angles from Fig. 12a the general orientation between real gas and perfect gas 428 in stage 1 do not show a signiﬁcant difference. This is independent whether choking or peak efﬁciency 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 efﬁciency for the peak 434 efﬁciency point. Polytropic Efﬁciency (dashed line) Work Coefﬁcient (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 efﬁciency) at sections close to the impeller outlets. For stage 1 Figure 13 shows the results for IP1 (Peak efﬁciency) at sections close to the impeller outlets. 435 Fig. 13 shows the results for IP1 (Peak efﬁciency) 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 ﬂow 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 ﬂow physics are not captured accurately. of blade geometries because the ﬂow 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 ﬂow 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 modiﬁcations. 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 ﬂow as the quantities. direct algebraic A transonic relations bump for 448 conditions and iterative procedures for the update of various ﬂow quantities. A transonic bump 449 case was chosen to demonstrate the consistency of the modiﬁcations. 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 modiﬁcations. 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 conﬁguration. 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 ﬂow 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 ﬂow 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 ﬂow physics was mainly driven by L.H. as it is a part of his PhD. The ﬁnal target of the thesis is a framework, able at predicting two-phase ﬂow 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 ﬂow 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 ﬁnancial contribution provided by the Swiss National Science Foundation (SNF: Grant-Number 175900). Conﬂicts of Interest: The authors declare no conﬂict 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 Speciﬁc heat at constant pressure [J/kgK] c Speciﬁc 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 ﬂow [kg/s] MR Maxwell relation PG Perfect Gas p Pressure [Pa] q Heat [J/kg] R Universal gas constant [J/kgK] R Speciﬁc 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 ﬂux [m /s] u Velocity x-component [m/s] v Velocity y-component [m/s] w Velocity z-component [m/s] V Volume [m ] v Speciﬁc volume [1/kg] z Compressibility Greeks a Thermal expansivity [1/K] b Compressibility [1/Pa] G Fundamental derivative g Heat capacity ratio h Polytropic efﬁciency L Work coefﬁcient 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 coefﬁcient 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 ﬁrst and second law of thermodynamics, the change in speciﬁc 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 ﬁrst 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. Speciﬁc Heat Capacity The speciﬁc heat at constant volume is deﬁned as given in Equation (A27). ¶e c = (A27) ¶T Using the deﬁnition 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 speciﬁc 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 speciﬁc 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 ﬂuid 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 ﬂow simulations: Inﬂuence of thermodynamic models on ﬂow ﬁelds 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.; Pﬁtzner, 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-Inﬂow 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 ﬂow ﬁeld. 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 ﬂows 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 modiﬁed Redlich-Kwong equation of state. Chem. Eng. Sci. 1972, 27, 1197–1203. 18. McBride, B.J.; Zehe, M.J.; Gordon, S. NASA Glenn Coefﬁcients 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 ﬂows. 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-ﬂuid approach for two or more immiscible ﬂuids. 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 ﬂow 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 ﬂow solver. Appl. Math. Comput. 2011, 217, 5041–5054. 29. Darwish, M.; Moukalled, F. A fully coupled Navier-Stokes solver for ﬂuid ﬂow 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 ﬁnite-volume solver for incompressible two-phase ﬂow. 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 ﬂows 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 ﬂow past an airfoil with trailing edge separation. AIAA J. 1983, 21, 1525–1532. 34. Darwish, M.; Sraj, I.; Moukalled, F. A coupled ﬁnite volume solver for the solution of incompressible ﬂows on unstructured grids. J. Comput. Phys. 2009, 228, 180–201. 35. Demirdžic, ´ I.; Lilek, Ž.; Peric, ´ M. A collocated ﬁnite volume method for predicting ﬂows 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 ﬂuid 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

http://www.deepdyve.com/lp/multidisciplinary-digital-publishing-institute/real-gas-models-in-coupled-algorithms-numerical-recipes-and-bwS90WaiZt

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

Loading next page...

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

Journal

International Journal of Turbomachinery, Propulsion and Power – Multidisciplinary Digital Publishing Institute

Published: Aug 3, 2020

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 †

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 †

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies