Design of supercritical carbon dioxide centrifugal compressors

Design of supercritical carbon dioxide centrifugal compressors PhD Thesis Seville, March 2014 PhD Candidate: Benjamín Monge Brenes Supervisor: Tomás Sánchez Lencero Co-supervisor: David Sánchez Martínez This page left intentionally blank AGRADECIMIENTOS (ACKNOWLEDGEMENTS) Al final de mi periodo de doctorado, me enfrento probablemente a la tarea más complicada de las que he realizado durante estos 4 años: los agradecimientos. Y es que es imposible expresar en una sola página mi gratitud hacia quienes me han brindado esta oportunidad de desarrollo profesional, pero sobre todo; y es lo mejor que me llevo; personal. En primer lugar deseo manifestar también mi agradecimiento a la Consejería de Economía, Innovación, Ciencia y Empleo de la Junta de Andalucía por haber financiado el proyecto de excelencia “Desarrollo de turbocompresores centrífugos de dióxido de carbono supercrítico para implementación en sistemas de producción de potencia limpios y con alta eficiencia” (referencia TEP09-5300) en la convocatoria 2009. Esta tesis se ha desarrollo en el seno de dicho proyecto”. Y por supuesto, agradecer a la Universidad de Sevilla el haberme concedido la beca de Plan Propio que he gozado y al Grupo de Máquinas y Motores Térmicos por confiar en mí como alumno aspirante. Especial agradecimiento merecen el director de dicho grupo, el Prof. Dr. D. Tomás Sánchez Lencero y el Prof. Dr. D. David Sánchez Martínez, tutor y co-tutor respectivos de esta tesis. La guía basada en su experiencia y conocimiento ha sido la clave en muchas ocasiones para el desarrollo de esta tesis. Por supuesto no me puedo olvidar del Prof. Dr. D. Ricardo Chacartegui Ramírez, cuyo asesoramiento y ayuda han sido fundamentales en el proceso de aprendizaje de mi doctorado, además de enriquecedores en otros trabajos conjuntos. Sin salir del terreno universitario, quiero expresar mi reconocimiento a la Cranfield University (Reino Unido) y en especial a los Prof. Pilidis y Savill, con quienes compartí un periodo breve, pero productivo y muy determinante para la elaboración de esta tesis. No puedo terminar mis agradecimientos en la Universidad de Sevilla sin nombrar a quienes me han aguantado en el día a día durante estos años: Los compañeros, pero sobre todo, amigos de los laboratorios. En especial a José Muñoz de Escalona, gracias a cuyo trabajo el mío ha resultado más fácil y fructífero. Además de sentirme un privilegiado por formar parte de este inmenso grupo humano, considero haber hallado con ellos algo tan valioso como esta propia tesis doctoral. Y ahora quiero mostrar mi agradecimiento a las personas que me han hecho posible encontrar esta posibilidad. Esas personas que con su apoyo, esfuerzo material e inmaterial y sobre todo su amor, le han dado y le siguen dando sentido y motivación a mi vida. Hablo por supuesto de mi familia: De la única mujer de mi vida, Lola, con quien llevo más de 14 años de vida compartidos y aún sigue impresionándome por las muestras de amor, cariño y sobre todo comprensión hacia mí. De mis padres, Manolo y Consuelo, y de mi hermana, Natalia, quienes me han inculcado los valores de humildad, compromiso y amor por la familia de los que espero nunca desprenderme. Gracias a que ellos me han procurado siempre lo mejor, he podido llegar hasta este punto. Hablo también de mis fieles amigos de toda la vida, porque ellos también han aportado a mi currículum para lograr este objetivo. Especialmente mi primo político, el único, el inigualable Tirado con su primo y su hermano. Y por último, esa otra gran parte de mi vida: Mi Banda, donde puedo desarrollar mi pasión por la música y la docencia a la vez. A toda mi familia… …a los que ya no están, a los que están y a los que estarán. This page left intentionally blank - i - Contents CHAPTER 1. Introductory concepts ..................................................................................................... 1 1.1 Scope of the work ...................................................................................................................... 1 1.2 Main constraints of CSP technologies ....................................................................................... 2 1.3 Open simple Brayton cycle as an alternative for solar power production ................................ 4 1.4 Open recuperative Brayton cycle as an alternative for solar power production ...................... 5 1.5 Why supercritical carbon dioxide Brayton cycle for solar power production? ......................... 7 1.6 The fluid: Carbon Dioxide ........................................................................................................ 13 1.7 The turbomachinery issue ....................................................................................................... 15 1.8 International research and commercial interest in S-CO2 technology.................................... 16 1.8.1 Czech Republic institutes .................................................................................................... 16 1.8.2 Tokyo Institute of Technology ............................................................................................ 16 1.8.3 SANDIA National Laboratories ............................................................................................ 17 1.8.4 National Renewable Energy laboratory .............................................................................. 18 1.8.5 Argonne National Laboratory ............................................................................................. 18 1.8.6 The University of Queensland, Australia ............................................................................ 18 1.8.7 International diffusion and investments............................................................................. 19 1.8.8 Grupo de Máquinas y Motores Térmicos de la Universidad de Sevilla (GMTS) ................. 19 CHAPTER 2. Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery...................................................................................................................... 20 2.1 Basics of CFD ............................................................................................................................ 20 2.1.1 Pre-processing .................................................................................................................... 20 2.1.2 The solver............................................................................................................................ 21 2.1.3 Post-processing ................................................................................................................... 21 2.2 Modelling material properties: equations of state.................................................................. 21 2.3 Modelling physics: turbulent flow ........................................................................................... 23 2.3.1 Turbulent flow treatment in CFD ........................................................................................ 26 2.3.2 Turbulence model applied to S-CO2 turbomachinery ......................................................... 27 2.4 Additional numerical considerations for S-CO2 simulation ..................................................... 30 2.5 Uncertainties associated to CFD simulation with S-CO2 .......................................................... 31 2.6 Summary and conclusions ....................................................................................................... 33 CHAPTER 3. Analysis of S-CO2 flows in elementary diffusion processes ............................................ 34 3.1 Definition and scope of the analysis ........................................................................................ 35 3.1.1 Diffuser operation ............................................................................................................... 35 3.1.2 Parameters of influence in conical diffuser performance .................................................. 37 - ii - 3.1.3 Presentation of results and range of study ......................................................................... 39 3.1.4 The tool: definition, validation and verification ................................................................. 40 3.2 Effect of geometry ................................................................................................................... 49 3.3 Effect of first order parameters ............................................................................................... 61 3.3.1 Effect of aerodynamic blockage.......................................................................................... 61 3.3.2 Effect of Mach and Reynolds numbers ............................................................................... 64 3.4 Effect of second order parameters .......................................................................................... 65 3.4.1 Effect of inlet distortion (velocity profile) ........................................................................... 65 3.4.2 Effect of swirl component ................................................................................................... 66 3.4.3 Effect of turbulence ............................................................................................................ 67 3.5 Summary and conclusions ....................................................................................................... 68 3.6 Further calculations towards more reliable turbulence modelling ......................................... 69 CHAPTER 4. One-dimensional model of supercritical carbon dioxide centrifugal compressors ......... 72 4.1 Review of one-dimensional centrifugal compressor models .................................................. 73 4.2 Model of performance of supercritical CO2 centrifugal compressors: direct problem ........... 77 4.2.1 Fundamentals of the model ................................................................................................ 78 4.2.2 Inducer model ..................................................................................................................... 79 4.2.3 Impeller model .................................................................................................................... 82 4.2.4 Vaneless space model ......................................................................................................... 92 4.2.5 Vaned diffuser model ......................................................................................................... 96 4.2.6 Volute model .................................................................................................................... 101 4.3 Validation of the model ......................................................................................................... 103 4.3.1 Validation for supercritical carbon dioxide ....................................................................... 103 4.3.2 Validation for air ............................................................................................................... 105 4.4 Design strategy for supercritical CO2 centrifugal compressors: Inverse problem ................. 107 4.4.1 Cycle definition ................................................................................................................. 107 4.4.2 Compression system design ............................................................................................. 109 4.4.3 Stage design ...................................................................................................................... 112 4.5 Summary and conclusions ..................................................................................................... 119 4.6 Justification of the research .................................................................................................. 120 CHAPTER 5. Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery ... 123 5.1 Geometry generation ............................................................................................................ 123 5.2 The computational model ..................................................................................................... 125 5.2.1 The mesh .......................................................................................................................... 125 5.2.2 Other numerical aspects of the computational model ..................................................... 130 5.3 CFD results ............................................................................................................................. 132 - iii - 5.4 Comparison between 3D-CFD and 1D results ....................................................................... 140 5.5 Summary and conclusions ..................................................................................................... 142 CHAPTER 6. Conclusions and future developments ........................................................................ 144 6.1 Conclusions from the comparative analysis of air and S-CO2 in conical diffusers ................. 144 6.1.1 Effect of geometry ............................................................................................................ 145 6.1.2 Effect of first order parameters ........................................................................................ 145 6.1.3 Effect of second order parameters ................................................................................... 145 6.2 Conclusions from S-CO2 turbomachinery design ................................................................... 146 6.3 Original contributions ............................................................................................................ 147 6.4 Future developments ............................................................................................................ 147 CHAPTER 7. References .................................................................................................................. 149 - iv - Figures Figure 1. Scheme of a Concentrating Solar Power (CSP) system. ................................................................. 2 Figure 2. Influence of receiver temperature (𝑇𝑟𝑒𝑐) on system efficiency (𝜂𝑠𝑦𝑠𝑡) for different geometric concentration ratios (𝐶𝑅), considering an irradiation of 𝜙=1000 W/m2, an optical efficiency of 𝜂𝑜𝑝𝑡=0.8, an ambient temperature of 𝑇𝑎𝑏𝑚=20 ºC and emissivity and absorptivity of 𝛼= 𝜀=0.85. ........................................................................................................................................................ 3 Figure 3. Temperature- Entropy diagram of a conventional open Brayton cycle. ....................................... 4 Figure 4. Effects of pseudo-presure ratio (𝛿) and temperature ratio (𝜃) on cycle efficiency (𝜂, in black) and specific work (𝑊, in red) for an ideal gas with 𝛾=1.4, compressor inlet temperature of 𝑇01=300 K and compressor and turbine isentropic efficiencies of 80 and 90 % respectively. ..................... 5 Figure 5. Temperature- Entropy diagram of an open simple recuperative Brayton cycle. .......................... 6 Figure 6. Effects of pseudo-presure ratio (𝛿) and temperature ratio (𝜃) on cycle efficiency (𝜂, in black) and specific work (𝑤, in red) for an ideal gas with 𝛾=1.4 in a recuperative cycle with 𝜀𝑟𝑒𝑐=0.85, compressor inlet temperature of 𝑇01=300 K and compressor and turbine isentropic efficiencies of 80 and 90 % respectively. ........................................................................................................... 7 Figure 7. Temperature- Entropy diagram of a Supercritical Brayton cycle .................................................. 7 Figure 8. Supercritical CO2 compressor impeller. Designed and manufactured by Barber Nichols Inc. for SANDIA National Laboratories facilities (Wright, et al., September 2010). ................................... 8 Figure 9. Recuperator effectiveness in different conditions in a recuperative supercritical carbon dioxide Brayton cycle for thermal equilibrium (Dostal, et al., March, 2004). ............................................. 8 Figure 10. Cumulative efficiency drops due to different equipment for both the simple cycle (in blue) and the recompression one (in red). ................................................................................................... 10 Figure 11. S-CO2 cycle configurations studied by Angelino in 1968 (Angelino, July 1968). ....................... 11 Figure 12. S-CO2 cycle configurations studied by Angelino in 1969 (Angelino, 1969). .............................. 12 Figure 13. Efficiency comparison amongst the supercritical cycle configurations studied by Angelino (Angelino, July 1968). ................................................................................................................... 13 Figure 14. Generalised compressibility factor diagram. ............................................................................. 13 Figure 15. Compressibility factor variation vs. density ratio with respect to the inlet conditions in compression processes for air water and supercritical CO2. ........................................................ 14 Figure 16. Carbon dioxide properties vs. temperature as a function of pressure. ..................................... 15 Figure 17. Experimental loop at Běchovice research institute. .................................................................. 16 Figure 18. Radial impeller of the S-CO2 cycle experimental facility at the Tokyo Institute of Technology. ...................................................................................................................................................... 17 Figure 19. Flow diagram of the closed compression test loop of SANDIA National Laboratories. ............. 17 Figure 20. Scheme of a S-CO2 experimental loop at Argonne National Laboratory for heat transfer studies. Not in operation nowadays. ......................................................................................................... 18 Figure 21. Scheme of a geo-thermosiphon based on S-CO2 cycle, object of study at the QGECE. ............ 19 Figure 22. Stages in pre-processing. ........................................................................................................... 21 Figure 23. Results of the simulations of the SANDIA S-CO2 radial compressor by (Pecnik, et al., 2012). The points (cells) where vapour-liquid equilibrium conditions are found are shown in red. ............. 23 Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 2 than 1000ºC with reduced heat losses. However, it must also be acknowledged that a technology enabling the current turbine inlet temperature levels of modern gas turbines (1500ºC) is not envisaged in the near future and, therefore, auxiliary firing will have to be implemented if high efficiencies are sought. In this context, the reference supercritical carbon dioxide power cycle considered in this work may contribute to fill the gap between steam and gas turbines, offering competitive thermodynamic performance, modularity and applicability in renewable energy systems but, currently, at a prohibitively high cost. 1.2 Main constraints of CSP technologies Transforming solar power into electricity requires an initial process of concentration, whose aim is generating high quality energy available for the prime mover or, widely, system that carries out the different steps of the energy conversion process: solar to heat to mechanical to, eventually, electrical energy. A higher concentring capacity brings about a higher receiver temperature and consequently a higher efficiency potentially attainable by the power conversion system. It is therefore an important parameter for the characterisation of solar facilities the so-called geometric concentration ratio, 𝐶𝑅= 𝐴𝑟𝑒𝑐/𝐴𝑐𝑜𝑙, which is simply the ratio between receiver aperture area and collector aperture area. Figure 1. Scheme of a Concentrating Solar Power (CSP) system. The heat transferred to the fluid per unit of receiver surface (𝑄𝑎𝑏𝑠/𝐴𝑎𝑏𝑠) is given by applying the energy conservation law (Romero-Álvarez & Zarza, May 2007): 𝑄𝑎𝑏𝑠/𝐴𝑎𝑏𝑠=𝛼·𝐶𝑅·𝜙·𝜂𝑜𝑝𝑡−𝜎·𝜀(𝑇𝑟𝑒𝑐 4−𝑇𝑎𝑚𝑏 4) [Equation 1.1] This heat is thus a function of the absorptivity (𝛼) and emissivity (𝜀) of the receiver, the temperatures of the hot source (𝑇𝑟𝑒𝑐) and cold sink (𝑇𝑎𝑚𝑏) and the specific thermal power (per unit area) impinging the receiver from the collector field, 𝜙·𝜂𝑜𝑝𝑡, where 𝜙 is the normal direct irradiation [W/m2] and 𝜂𝑜𝑝𝑡 the optical efficiency of the facility, which accounts for the losses taking place from the solar energy striking the solar field to that hitting the receiver surface. The receiver efficiency is thus expressed as: 𝜂𝑎𝑏𝑠= Power absorbed Power received =𝑄𝑎𝑏𝑠/𝐴𝑎𝑏𝑠 𝐶𝑅·𝜙·𝜂𝑜𝑝𝑡=𝛼−𝜎·𝜀(𝑇𝑟𝑒𝑐 4−𝑇𝑎𝑚𝑏 4) 𝐶𝑅·𝜙·𝜂𝑜𝑝𝑡 [Equation 1.2] The power absorbed by the receiver is converted into electrical energy (in a series of energy conversion processes) by a power block operating between a hot source and a cold sink. As known, given a hot and cold temperatures for these source and sink, the practical heat to mechanical energy conversion efficiency is limited by that of the Carnot Cycle: 𝜂𝐶𝐴𝑅𝑁𝑂𝑇=(𝑇𝑟𝑒𝑐−𝑇𝑎𝑚𝑏)/𝑇𝑎𝑚𝑏. Thus, the total ideal efficiency of the receiver and (ideal) power block is expressed as the product of both, 𝜂𝑠𝑦𝑠𝑡=𝜂𝑟𝑒𝑐· Hot source, Trec [K] Cold source, Tamb [K] Compression Expansion SOLAR POWER, ϕ [W/m2] Acol [m2] Arec [m2] Qabs [W] We [W] Qamb [W] Optical losses, ηopt Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 3 𝜂𝐶𝐴𝑅𝑁𝑂𝑇, making it evident that for given irradiation and optical efficiency, increasing the receiver temperature (peak temperature of the working fluid) has counteracting effects on the global efficiency of the plant. On one hand, it increases the thermal efficiency of the power cycle while, on the other, it also increases the radiation losses from the receiver (i.e. heat radiated from the receiver surface to the surroundings). This combined effect is shown in figure 2 where, for a concentration ratio, an optimum receiver temperature is found in terms of global efficiency. For temperatures higher than the optimum one, the receiver radiation losses are not compensated for by the higher cycle efficiency. Figure 2. Influence of receiver temperature (𝑇𝑟𝑒𝑐) on system efficiency (𝜂𝑠𝑦𝑠𝑡) for different geometric concentration ratios (𝐶𝑅), considering an irradiation of 𝜙=1000 W/m2, an optical efficiency of 𝜂𝑜𝑝𝑡=0.8, an ambient temperature of 𝑇𝑎𝑏𝑚=20 ºC and emissivity and absorptivity of 𝛼= 𝜀=0.85. For the most developed solar technologies, geometric concentration ratios vary between 30 and 80 for parabolic trough and between 200 and 1000 for central receiver plants, which limits the system efficiency to values below 60 % even if the cycle works ideally, as shown in figure 2. Moreover, the real power blocks operate very far from Carnot efficiency so their global performance is further reduced. Consequently, in terms of project economics, the Levelised Cost Of Electricity (LCOE in €/kWh) measuring the cost of each unit electrical energy obtained, all lifelong capital and operating expenditures of the plant into account, departs from that of the hypothetical ideal plant. LCOE values around 0.12-0.15 €/kWh for parabolic trough (Turchi, et al., 2010) and 0.09-0.12 €/kWh for central receiver plants (Prior, 2011) have been reported. These account for all the expenditures and investment costs that need to be faced along the plant lifetime operation, 𝐿𝐶𝑂𝐸=(∑𝐼𝑛𝑣𝑒𝑠𝑡𝑚𝑒𝑛𝑡+ ∑𝐸𝑥𝑝𝑒𝑛𝑠𝑒𝑠)/∑𝑃𝑟𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛, and it is therefore an important parameter to appraise the economic feasibility of a project. In a simple model, one can express:  The capital investment as the product of the specific capital cost and the power installed: ∑𝐼𝑛𝑣𝑒𝑠𝑡𝑚𝑒𝑛𝑡=𝑐𝑖𝑛𝑣·𝑊󰇗,  The operational expenditures as the sum of financial (debt services), operating and fuel cost, the latter being null for a purely solar application (no fuel hybridisation): ∑𝐸𝑥𝑝𝑒𝑛𝑠𝑒𝑠= 𝑋𝑑𝑒𝑏𝑡·∑𝐼𝑛𝑣𝑒𝑠𝑡𝑚𝑒𝑛𝑡+𝑘1·𝑐𝑂𝑃&𝑀·∑𝑃𝑟𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛,  And the production as a function of the lifetime (𝑁 in years), plant load factor (𝑓𝑙𝑜𝑎𝑑) and rated output (𝑊󰇗): ∑𝑃𝑟𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛=8760·𝑘2·𝑁·𝑓𝑙𝑜𝑎𝑑·𝑊󰇗. Manipulating the previous expressions, LCOE can be expressed as: LCOE= 𝑐𝑖𝑛𝑣·(1+𝑋𝑑𝑒𝑏𝑡) 8760·𝑘2·𝑁·𝑓𝑙𝑜𝑎𝑑·𝑊󰇗+𝑘1·𝑐𝑂𝑃&𝑀 [Equation 1.3] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 300 500 700 900 1100 1300 System effiency, ηsyst Receiver tempeature, Trec (K) CR = 1 CR = 50 CR = 10 CR = 100 CR = 200 CR = 500 CR = 1000 CR = 2000 Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 4 Where 𝑐𝑖𝑛𝑣 and 𝑐𝑂&𝑀 are the specific investment and operational/maintenance costs in €/kW and €/kWh respectively, 𝑘1 and 𝑘2 are constant coefficients that bring both maintenance costs and production to their present values at the initial moment of the project and 𝑋𝑑𝑒𝑏𝑡 the ratio between the total financial expenses and the initial investment. The new proposal can only influence the specific costs in a partial way, since a large portion of costs is related not to the power block but to the solar field (where there would be a reduction only if the efficiency of the power block were higher). It is acknowledged here that costs depend largely on the maturity of the technology and hence commercial development is necessary for these costs to be accurate; nevertheless, attention must also be paid to other aspects that permit reducing the LCOE: simplicity and compactness of equipment and layouts, possibility of integrating energy storage systems and good partial load performance are some of them. The last limitation to highlight regarding solar power plants based on Rankine cycles is the dependence on water, necessary for cooling (note that air cooled condenser reduce the consumption of water but at the expense of a lower efficiency and therefore larger solar field) and also for cycle make up. Typical yearly water consumptions for Rankine power plants are around 2 L/kWh (Walker, et al., 2013), establishing an important handicap for the location of these plants in desert areas, where the solar resource is highest but hardly any water is available. 1.3 Open simple Brayton cycle as an alternative for solar power production The use of Brayton instead of Rankine cycles could offer some advantages such as no water consumption, a plant footprint considerably reduced thanks to its higher specific power, a priori similar thermal efficiency to state-of-the-art Rankine-based plants. However, both efficiency and specific work are preconditioned by the working fluid and operating parameters: pressure ratio and TIT mainly. Figure 3. Temperature- Entropy diagram of a conventional open Brayton cycle. Figure 3 shows the four characteristic stagnation points of a conventional Brayton cycle in a temperature-entropy diagram. The heat injected per unit mass is 𝑄=ℎ03−ℎ02 and the gross specific shaft work 𝑊=(ℎ03−ℎ04)−(ℎ02−ℎ01), both in J/kg. The efficiency is the quotient of both (𝜂=𝑊/𝑄) and it can be expressed as shown in [equation 1.4] by using the definition of the pseudo-pressure ratio 𝛿=𝜌 (𝛾−1)/𝛾 and the temperature ratio 𝜃=𝑇03/𝑇01, where 𝜌 is the pressure ratio and it has been assumed ideal gas and no pressure drops. 𝜂=𝜃𝛿𝜂𝐶𝜂𝑇−1 𝜃−1 𝛿−1𝜂𝐶−1 [Equation 1.4] 𝑊=𝐶𝑝𝑇01(𝛿−1) 𝜂𝐶(𝜃𝛿𝜂𝐶𝜂𝑇−1) [Equation 1.5] The equations above are obtained for an ideal cycle under the following assumptions:  The working fluid in the entire cycle is air which is also assumed to behave ideally. 01 02 03 04 01 02 03 04 T s Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 5  There are no pressure drops.  The fuel flow rate is negligible with respect to air one.  Non-isentropic behaviour of compressor and turbine. This simple model yields two well-known general conclusions for Brayton cycles (Figure 4):  There exist two different optimum pressure ratios for a given temperature ratio: one for highest efficiency and another for highest specific work, the former being always higher than the latter. The choice of one of these options depends on the application.  The achievable efficiency is upper limited by the temperature ratio, increasing in a directly proportional manner. Figure 4. Effects of pseudo-presure ratio (𝛿) and temperature ratio (𝜃) on cycle efficiency (𝜂, in black) and specific work (𝑊, in red) for an ideal gas with 𝛾=1.4, compressor inlet temperature of 𝑇01=300 K and compressor and turbine isentropic efficiencies of 80 and 90 % respectively. It is worth noting that the curves in Figure 4 do only apply to the rated operating conditions; i.e. they are not valid for varying operating conditions at part-load or off-design operation. Thus, for the characteristic range of TIT in solar applications, cycle efficiency is restricted to values under 30 % for temperature ratios not higher than 3.5 (TIT<1050 K), which is under the expectation of the new proposal. This can be improved by implementing several modifications in the conventional Brayton cycle lay-out, or by using a combustion chamber in series with the solar receiver to increase the temperature from receiver outlet to turbine inlet. However, this second option contradicts the renewable nature of any solar facility and it is therefore discarded. There exits two other options worth considering though: the first one is to make the cycle recuperative (Section 1.4) and the second one is the possibility to change the working fluid in order to operate near its critical points, giving place to supercritical Brayton cycles (Section 1.5). 1.4 Open recuperative Brayton cycle as an alternative for solar power production The recuperator located between compressor and heater recovers part of the heat content of the turbine exhaust which is used to heat up the compressor delivery air (Figure 5). Consequently, the heat required to reach a certain TIT decreases and therefore, the cycle efficiency is improved. 0 50 100 150 200 250 300 350 400 0 0.1 0.2 0.3 0.4 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Specific work, w(kJ/kg) Cycle efficiency, η Pseudo-pressure ratio, δ θ=5.0 θ=4.5 θ=4.0 θ=3.5 θ=3.0 θ=2.5 θ=2.0 θ=5.5 Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 6 Figure 5. Temperature- Entropy diagram of an open simple recuperative Brayton cycle. Points 02p and 04p represent the high and low pressure outlet streams of the recuperator and, based on them, the recuperator effectiveness can be defined As 𝜀𝑟𝑒𝑐=(ℎ02𝑝−ℎ02)/𝑄𝑚𝑎𝑥 or 𝜀𝑟𝑒𝑐= (ℎ04−ℎ04𝑝)/𝑄𝑚𝑎𝑥, since due to energy conservation (ℎ02𝑝−ℎ02)=(ℎ04−ℎ04𝑝). The maximum heat is simply calculated as 𝑄𝑚𝑎𝑥=ℎ02−ℎ04 for substances with constant specific heat such as ideal gases or liquids. However, if the specific heat varies strongly from across each side of the heat exchanger, as it happens with supercritical fluids; the maximum heat is generally calculated as exposed by (Dostal, et al., March, 2004): 𝑄𝑚𝑎𝑥=𝑀𝐼𝑁{(ℎ04−ℎ04∗),(ℎ02∗−ℎ02)} [Equation 1.6] Where point 02* is at same pressure as 02p and at the temperature of point 04, and equivalently point 04* is at same pressure as than 04p and at the temperature of 02. Therefore, these points imply thermal equilibrium at one of the two recuperator ends (that with lower specific heat). If the effectiveness of the recuperator is now combined with the aforementioned nondimensional parameters 𝛿 and 𝜃, the new cycle efficiency is defined as in [equation 1.7] from which [equation 1.4] is obtained for the non-recuperated cycle (𝜀𝑟𝑒𝑐=0). On the contrary, the equation yielding the specific work of the cycle does not vary with respect to the previous expression [Equation 1.5] as long as the simplifying assumptions are still considered. 𝜂= 𝜃𝛿𝜂𝐶𝜂𝑇−1 𝜃−1 𝛿−1𝜂𝐶−1−𝜀𝑟𝑒𝑐(𝜃−𝛿2) 𝛿(𝛿−1)𝜂𝐶 [Equation 1.7] A parametric representation of this expression is presented in figure 6, where a recuperator effectiveness equal to 85 % is considered (which is a technological value for air recuperator effectiveness according to literature review (Muñoz de Escalona, et al., 2010)). In addition to a cycle efficiency enhancement with respect to the simple cycle (Figure 4), another effect of recuperation is the displacement of the pressure ratio for maximum efficiency towards lower values, because increasing the pressure ratio in such configuration implies (for a given TIT) reducing the recuperative potential as the temperature gap between 04 and 02 decreases. As shown, it is possible to exceed 35 % of efficiency with this level of recuperation (Figure 6), a value that is close to the expectations for a gas turbine based application even though another alternative is still studied in the next section: the supercritical Brayton cycle. 01 02 03 04 02p 04p T s 01 02p 03 04 04p 02 Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 7 Figure 6. Effects of pseudo-presure ratio (𝛿) and temperature ratio (𝜃) on cycle efficiency (𝜂, in black) and specific work (𝑤, in red) for an ideal gas with 𝛾=1.4 in a recuperative cycle with 𝜀𝑟𝑒𝑐=0.85, compressor inlet temperature of 𝑇01=300 K and compressor and turbine isentropic efficiencies of 80 and 90 % respectively. 1.5 Why supercritical carbon dioxide Brayton cycle for solar power production? Non-conventional fluids in supercritical or transcritical cycles are usually linked to closed or externally heated configurations as shown in Figure 7. A cooling system is therefore needed which is a non-desirable feature for the new proposal as long as it involves water consumption. Nevertheless, the motivation of using supercritical cycles is mainly twofold: Figure 7. Temperature- Entropy diagram of a Supercritical Brayton cycle  On one hand, the real gas behaviour provides the fluid with a lower compressibility factor that makes the compression process less energy-consuming. In fact, even with a lower expansion work at the turbine, the net power output of still high as it is compensated for by the reduction in compression work. The cycle thus becomes the only practical alternative to improve cycle efficiency under the peak temperature limitation found in solar applications (without significantly affecting the layout of the plant).  On other hand, the highly pressurised conditions in the system reduce the footprint of the power block. In fact, the graphical comparison shown in (Dostal, et al., March, 2004) exposes very clearly the size reduction attainable with supercritical CO2 turbomachinery in comparison to air (Brayton) and water (Rankine) turbines. The first supercritical CO2 turbomachinery design (Figure 8) came to demonstrate the surprising size reduction that supercritical technologies can make possible in relation to conventional systems. This is a 0 50 100 150 200 250 300 350 400 0 0.1 0.2 0.3 0.4 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Specific work, w(kJ/kg) Cycle efficiency, η Pseudo-pressure ratio, δ θ=5.0 θ=4.5 θ=4.0 θ=3.5 θ=3.0 θ=2.5 θ=2.0 θ=5.5 εrec=0.85 01 02 03 04 02p 04p T s 01 02p 03 04 04p 02 Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 8 very important aspect, since it should bring about cost reductions in a future scenario of technical maturity and, at the same time, the expected compactness can make it very suitable for modular applications (so common in renewable and waste energy applications). Figure 8. Supercritical CO2 compressor impeller. Designed and manufactured by Barber Nichols Inc. for SANDIA National Laboratories facilities (Wright, et al., September 2010). One of the first individuals to research the supercritical cycle was Feher (Feher, 1968). This author demonstrated the need to adopt recuperative configurations to reach high efficiencies, since the heat recovery between high and low pressure streams can be performed with high effectiveness. However, recuperation finds a particular problem in supercritical cycles called the pinch problem, which is associated to the enormous differences in isobaric specific heat that can exist between the cold and hot ends and the high and low pressure sides of the recuperator. Feher also exposed it very clearly in an enthalpytemperature diagram, which was later manipulated by Dostal (Dostal, et al., March, 2004) to present the effect of this circumstance on the effectiveness of the recuperator (Figure 9). Even if the minimum temperature difference is set to 0ºC at one recuperator end, there is still a region, called the pinch-point region by Dostal, where the effectiveness is lower and does not attain 100%. Figure 9. Recuperator effectiveness in different conditions in a recuperative supercritical carbon dioxide Brayton cycle for thermal equilibrium (Dostal, et al., March, 2004). This particular performance of the working fluid in the heat exchangers at low temperature (close to the critical point) brings about a subsequent penalisation on cycle efficiency. Therefore, compact and highly efficient recuperators are mandatory in order not to incur in an excessive deterioration of the global efficiency of the system though this issue does not represent a technological barrier today. In effect, Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 9 printed circuit heat exchangers (PCHE) can achieve efficiencies of up to 95% easily (Muñoz de Escalona, et al., 2010), in particular for CO2 which is a very advantageous fluid for supercritical cycles thanks to its critical conditions (31ºC, 73.8 bar), Table 1. Fluid Critical temperature (K) Critical pressure (bar) Air 132.5 37.9 Water (H2O) 647.1 220.6 Ammonia (NH3) 405.4 113.3 Carbon dioxide (CO2) 304.1 73.8 Perfluoropropane (C3F8) 345.0 2.6 Sulfur dioxide (SO2) 430.6 7.9 Sulfur hexafluoride (SF6) 318.7 3.8 Xenon (Xe) 289.7 5.8 Propane (C3H8) 369.9 42.5 n-Butane (C4H10) 425.1 38.0 Isopentane (C5H12) 460.4 33.8 Table 1. Critical properties of different fluids (Lemmon, et al., 2010). There is a second aspect that needs to be studied in addition to the selection of the working fluid. Depending on the pressure range, (Brayton) power cycles can be differentiated in transcritical, when the critical pressure of the working fluid is met somewhere in the compression process, and supercritical, when the compression process starts at a pressure above the critical point 1 . In the transcritical layout, a new classification is needed depending on whether the working fluid is in liquid state at the inlet to the compression device (pump in this case) or if it is in gaseous state (compressor). The former case is usually termed as the supercritical condensation cycle which was extensively studied by Prof. Angelino in his seminal work (Angelino, July 1968) which where through a complete evaluation of several layouts this author presented a body of knowledge so as to which supercritical cycle is more efficient and where the inefficiencies are located. Figures 11 and 12 are taken from the aforecited reference and show the different layouts studied by Angelino. Angelino pointed out that most of the cycle’s irreversibility is found in the recuperator, as a result from the fact that heat transfer occurs between streams with different heat capacities (Angelino, July 1968). Thus, this author proposed different layouts aiming at reducing these sources of efficiency drops. The most important contribution of that work is the recompression cycle (Cycle A in Figure 11), which manages to reduce the difference between heat capacities in this process. To that aim, this layout comprises two compressors (main and recompression) and two recuperators (high and low temperature recuperator, HTR and LTR). Differences in specific heats between both sides of the HTR are negligible, since the real behaviour of the gas tends to disappear at elevated temperature. Therefore, similar heat capacities will be obtained as long as mass flow rates are kept equal in both sides of the HTR. However, in the LTR, due to more noticeable differences in the specific heat between streams, the hot outflow of the LTR is split: a fraction is directed to recompression and the remaining flows to the condenser, thus attaining similar heat capacities in both sides of the LTR. Irreversibility reductions are verified in Figure 10, where cumulative efficiency drops are represented for the different equipment of the cycle and their consequent losses, taking the Carnot 1 This concept can be extended to temperature (supercritical cycles being those operating at a higher than critical temperature) though this is not common. Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 10 efficiency between maximum and minimum cycle temperatures (𝜂𝐶𝐴𝑅𝑁𝑂𝑇=1−𝑇𝑚𝑖𝑛/𝑇𝑚𝑎𝑥) as a reference. Same conditions and cycle parameters as in (Angelino, July 1968) are considered for this calculation and it is confirmed that Angelinos’s recompression cycle manages to reduce the overall cycle irreversibility noticeably. Figure 10. Cumulative efficiency drops due to different equipment for both the simple cycle (in blue) and the recompression one (in red). Starting off from this layout, other cycle configurations were suggested by Angelino with the objective of further reducing irreversibilities as well as avoiding other technical problems. Cycle layout B is the pre-cooling cycle and makes the turbine exhaust pressure independent from condensation pressure. This layout comprises an additional cooling process from points 9 to 10, which provides additional waste heat available for integration opportunities. Layout C belongs to the reheating cycle, and presents another alternative to layout A which reduces the stresses in the high temperature devices. The pre-compression cycle (layout D) is a last option similar to the pre-cooling cycle but without additional cooling. In this case, the pre-compression enhances the potential to recuperate heat in the LTR due to the temperature increase of point 10, however this cycle does not contain any recompression. 100 150 200 250 300 5 10 15 20 25 30 35 40 45 Turbine Inlet Pressure [bar] hcycle - hCARNOT [%] Qheat & Qcool Qheat & Qcool Recompressor Pump Compressor Turbine Turbine LTR HTRRecuperator Pressure Drops Pressure Drops Simple cycle Recompression cycle Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 11 Figure 11. S-CO2 cycle configurations studied by Angelino in 1968 (Angelino, July 1968). Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 18 1.8.4 National Renewable Energy laboratory In 2010, the National Renewable Energy laboratory (NREL) started working in this field (Ma & Turchi, 2011) and nowadays, the Head of the CSP department at NREL, Dr. Craig Turchi, is world-widely acknowledged as an active leader of the S-CO2 community. NREL is participating in the development of a 16 M$ Project financed by the Department Of Energy (DOE) of the US in the scope of the Sunshot initiative (NREL, 2012) (see section 1.8.7). 1.8.5 Argonne National Laboratory The Argonne National Laboratory (ANL) from the U.S. Department of Energy is actively working on the S-CO2 cycle oriented to nuclear applications. Heat transfer studies were initially done in the facility shown in Figure 20, which has two separated S-CO2 loops thermally linked by the heat exchanger: the lower pressure loop operates at around 80 bar and the higher pressure loop at up to 200 bar. CO2/CO2 and water/CO2 heat transfer phenomena were studied aimed at developing design techniques for the recuperator and the cooler of the S-CO2 power cycle. Figure 20. Scheme of a S-CO2 experimental loop at Argonne National Laboratory for heat transfer studies. Not in operation nowadays. After these heat transfer studies, ANL is currently developing an experimental project on Sodium- Cooled Fast Reactors (SFR) in a consortium with Commissariat a l’Energie Atomique (CEA) of France and Japan Nuclear Cycle Development Institute (JNC) (Chang, et al., 2005). This type of advanced nuclear reactor is nowadays conceived as one of the most promising options in the mid-long term beyond 2020 and this project would cover all the experimental stages regarding equipment for the S-CO2 cycle oriented to nuclear applications. The Idaho National Laboratory is also contributing to this project by investigating the material compatibility for the so-called Very High Temperature Reactors (VHTR) (Chang, et al., 2006). 1.8.6 The University of Queensland, Australia In 2009, the Queensland government founded the Queensland Geothermal Energy Centre of Excellence (QGECE) with an initial investment of 15 M$. Amongst their different initiatives, the development of turbines for supercritical fluids, CO2 amongst them (Haghshenas Fard, et al., 2010), is worth mentioning. In this case, the application studied consists of an advanced geothermal system based on a CO2 geo-thermosiphon that extracts heat from rock layers at a depth of 5 km, thus heating up a fluid that is later expanded in the turbine, cooled down and pumped again towards the hot zone in the subsoil. Introductory concepts Design of supercritical carbon dioxide centrifugal compressors 19 Figure 21. Scheme of a geo-thermosiphon based on S-CO2 cycle, object of study at the QGECE. 1.8.7 International diffusion and investments In terms of dissemination of research, important achievements must also be highlighted such as the Supercritical CO2 Power Cycle Symposium, organised every two years by the S-CO2 community, as well as the creation of the technical committee number 34, Supercritical CO2 power within the International Gas Turbine Institute (IGTI) of ASME. In addition to that, the ASME-IGTI decided to take on the organisation of the 2nd International Symposium on ORC Power Systems, which resulted from a recent European initiative driven by Professor Piero Colonna (Technical University of Delft, Netherlands). This event took place last October in Rotterdam. Nowadays, it seems that this intense investigation of the S-CO2 cycle has provided sufficient knowledge and gathered enough interest to start thinking about addressing its commercial deployment, in order to which important governmental support plans are being put forth in this field. In this sense, the SunShot Initiative promoted by the Department of Energy of the United States has fixed the objective of the reducing the cost of electricity to $0.06/kWh for Concentrating Solar Power (CSP) technologies, for which 21 projects totalling an investment of $54.7 million have been put in place for. Amongst this, $8 million are being devoted to the construction of a S-CO2 power cycle by the National Renewable Energy Laboratory, NREL (NREL, 2012). 1.8.8 Grupo de Máquinas y Motores Térmicos de la Universidad de Sevilla (GMTS) In this scope, GMTS is trying to contribute to the supercritical carbon dioxide community, not only with theoretical and numerical studies, but also by participating in experimental projects such as the one developed in consortium with ALTRAN TECHNOLOGIES, S.A. and aimed at designing a wind tunnel with pressurised CO2 and a S-CO2 modular design where both simple and recuperative cycle will be possible to analyse as well as radial compressor stages in stand-alone loops. Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 20 CHAPTER 2. Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Computational Fluid Dynamics (CFD) can be defined as a mathematical technique to obtain more or less realistic solutions to real problems involving physical-chemical fluid phenomena. The realism of the solution depends not only upon the validity of the models employed, but also on the skills from the user, solely acquired by experience on simulation. Although experiments must be seen as the most creditable method to get data, CFD is increasingly becoming more widely implemented in industry whether as part of the design processes or for research, as it exhibits advantages with respect to high-quality experimental facilities. The following are worth noting:  Although investment costs of a high-performance computing facility (both for acquiring CFD software and computing hardware) are not negligible, they can be drastically reduced with respect to test rigs. Moreover, a reduction of more than an order of magnitude must be expected for large systems (for instance, gas turbine engines)  State-of-the-art numerical codes are versatile and allow for the implementation of multiple physical phenomena (aerodynamic, hydrodynamic, chemical reactions, heat transfer, etc.), the so-called multi-physics, thus recreating the real operating conditions of the experimental facility but in a virtual environment.  Tests of the system in hazardous operating conditions can be developed without personal or material danger.  Thanks to the great level of detail with which results can be obtained, CFD allows obtaining data where, perhaps, instrumentation could never have been installed. Based on these advantages, which are even more important in the scenario of economic crisis that we are currently going through, CFD simulation in ANSYS FLUENT® has been employed as the principal tool to develop this project, since supercritical carbon dioxide experimental facilities are nowadays extremely costly. Nonetheless, the few experimental results available in the field are considered as references for validation. 2.1 Basics of CFD It is well-acknowledged that commercial CFD software, including ANSYS FLUENT®, employs the finite volume method (Versteeg & Malalasekera, 2007), which forms part of the finite differences techniques and consists of solving the flow conservation equations in the differential control volumes that compose the entire continuum. This is carried out by the solver and it is just a part of the problem since CFD also involves other two additional stages: pre-processing and post-processing which take place and after the solver completes the calculations. 2.1.1 Pre-processing Pre-processing comprises all the stages previous to the calculation, accordingly defining the necessary inputs, models and methods to use by the solver. Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 21  Geometry: creation of the computational domain by using CAD software.  Mesh generation: discretization of the geometry into elementary cells.  Selection of physical-chemical models in accordance with the phenomena of interest: turbulence, radiation, chemical reaction, multiphase flow, etc.  Material: specification of properties.  Definition of the boundary conditions.  Setting up the numerical aspect to be used by the solver in a way that they do not have an effect on the final solution but, on the contrary, only help to improve the computation process. However, this can become crucial for convergence in highly complex problems such as the precise case of the majority of supercritical carbon dioxide simulations. Figure 22. Stages in pre-processing. 2.1.2 The solver Resolution is obviously the heart of CFD and, as said before, the software employed uses the finite volume method, whose algorithm is based on the following steps: (i) Integration of the conservation equations on each cell, (ii) discretisation of equations and (iii) resolution of the algebraic equation system in an iterative manner (for more details on solver operation (Versteeg & Malalasekera, 2007), (Ferziger & Peric, 2002)). 2.1.3 Post-processing Once a solution has converged to an acceptable level, all the relevant magnitudes of the flow field are available to be screened in a number of formats: numerical reports, contour and vector plots, 2D graphs, streamlines and even animation when the case is time-dependent. While experimental facilities require high-quality instrumentation that are specific for each magnitude and only acquire data for a particular location, post-processing a CFD solution is so powerful that permits getting any data at any location within the flow domain. Moreover, this is done without perturbing the flow as the instrument would do. On the negative side tough, it must be noted that the fidelity of these outputs is strongly dependent on how the pre-processing is done so that an incorrect pre-processing brings about unreal post-processing outputs. 2.2 Modelling material properties: equations of state The real behaviour of the fluid through the compressor (see section 1.6) manifests as low compressibility factors (see Figure 15), meaning that ideal gas models do not describe the state of the fluid accurately. In this context, real gas effects with compressibility factors slightly lower than unity can be modelled accurately enough with surrogate ideal gas models (Takagi, et al., 2010), i.e. modifying the molar mass of the gases in order to adjust density to its actual value. Should a stronger real gas behaviour occur, this could be modelled by a cubic Equation Of State (EOS) for certain ranges of pressure and temperature, but still they fail to evaluate the properties of the fluid in the vicinity of the critical point. Modelling the flow in this region usually requires a dedicated EOS (for each fluid). Several attempts to model carbon dioxide have been made since 1970. A compendium of them is presented in (Span & Wagner, 1996) and summarised in table 2. Geometry creation Mesh generation MODELS Turbulence Radiation Reaction Model selection Material definition Boundary conditions URF = 0.7 P ϵ [1, 2] bar T ϵ [0, 100] ºC ... Numerics Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 22 Authors Year Temperature range (K) Pressure range (MPa) Bender 1970 216-1076 0-50 Altunin and Gadetskii 1971 215-1300 0-300 Stein 1972 Non published Non published Starling et al. 1972 243-413 0-48 Meyer-Pittroff 1973 200-1273 0-60 Angus et al. 1976 220-1100 0-100 Huang et al. 1985 216-423 0-310 Ely 1986 216-1023 0-300 Ely et al. 1987 216-1023 0-300 Pitzer and Schreiber 1988 230-1030 0-100 Ely et al. 1989 216-1023 0-316 Pitzer and Sterner 1994 220-2000 0-10000 Table 2. Equations of state developed for carbon dioxide. Extracted from (Span & Wagner, 1996). Since 1996, it is worldwide accepted that the most precise EOS for modelling carbon dioxide is the one developed by Span and Wagner (Span & Wagner, 1996).This equation is based on Helmholtz energy (𝐴) as a function of density (𝜌) and temperature (𝑇) in their dimensionless forms: 𝜙=𝐴/𝑅𝑇, 𝛿= 𝜌/𝜌𝑐 and 𝜏=𝑇𝑐/𝑇 respectively. As a function of these three variables and its partial derivatives, Span & Wagner presented all the functions and coefficients necessary to calculate the following list of thermodynamic properties: pressure, entropy, enthalpy, internal energy, isochoric and isobaric heat capacity, saturated liquid heat capacity, speed of sound, Joule-Thompson coefficient, fugacity and second and third virial coefficients. Nonetheless, the EOS solves part of the problem, as transport properties need to be modelled as well, especially viscosity and conductivity, both of which participate in the CFD simulation through the momentum and energy equations. There exists consensus also for this and thus complementary information to Span & Wagner’s EOS are needed; for instance the works by Fenghour et al. (Fenghour, et al., 1998) and Vesovic et al. (Vesovic, et al., 1990). This set of models is recommended by the National Institute of Standard and Technology (NIST) in its Standard Reference Data (SRD) (NIST, Nov. 2012) for modelling carbon dioxide properties. This approach is employed by REFPROP®, the reference software for the simulation of fluid properties developed by this institution. The implementation of this set of equation into ANSYS Fluent® can be complex and tedious and thus some researchers have opted for parametric (or look-up) tables which have been employed mostly in proprietary codes. Amongst this, the works by Colonna’s team stands out (Pecnik, et al., 2012). Parametric tables constitute powerful tool given the considerable reduction in elapsed time that can be achieved given that they only involve polynomial interpolations in 2D. Fortunately enough for this work, ANSYS Fluent®, in addition to offering the possibility to use parametric tables by means of User Defined Real Gas Model (UDRGM) functions, also contains a database of 83 pure fluids modelled according to NIST (as included in REFPROP®). One of these fluids available is carbon dioxide which is modelled satisfactorily Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 23 in the ranges 200-100 K, 0-100 MPa 3 . This fluid model is activated in ANSYS Fluent by the following text command: define/user-defined/real-gas-model/nist-real-gas-model/yes/”co2.fld” NIST Real Gas Models (NRGM) has nonetheless limitations in ANSYS Fluent®, amongst which the impossibility to model two-phase flow must be cited. It must be noted this feature may have a significant impact on the project since condensation is a phenomenon likely to appear in supercritical fluid turbomachinery, as discussed by (Pecnik, et al., 2012). Even though it might appear partial load operation, condensation is a priori a non-desirable phenomenon in turbomachinery design (as it is the case in an LP steam turbine). In consequence, from a design standpoint, no two-phase flow will need to be simulated as long as the design space is constrained to condensation-free regions. Figure 23. Results of the simulations of the SANDIA S-CO2 radial compressor by (Pecnik, et al., 2012). The points (cells) where vapour-liquid equilibrium conditions are found are shown in red. In order to explore how the EOS affects the results, section 3.1.4.4 shows the evaluation of different options to model S-CO2 in a conical diffuser (two-dimensional simulations). Even though small differences were found in the estimates of the representative performance parameters of the diffuser, such as blockage factor and pressure rise coefficient, much larger differences were observed in fluid properties like density and speed of sound, and consequently in Mach number. Therefore, it was concluded that cubic EOS are not appropriate for the purpose of designing turbomachinery where fluid properties play a crucial role. Consequently, the NRGM for carbon dioxide is employed for simulation. 2.3 Modelling physics: turbulent flow As the majority of industrial applications of interest, conventional turbomachinery problems involve turbulence flows because of the typically high values of Reynolds number, which are still higher for supercritical carbon dioxide due to its high density in comparison with air (Figure 16). Therefore, turbulence is the main physical phenomenon to be concerned about in the CFD simulation of interest for this study. In nature, turbulent flow is non-steady and heterogeneous in the three spatial coordinates; in fact, Wilcox pointed out that “it is characterized by the presence of a large range of excited length and 3 These are the ranges of validity for the thermal conductivity model employed by (Vesovic, et al., 1990), which is the most restrictive one in this case. Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 24 time scales”. A simpler manner to describe turbulent flow is by the Reynolds decomposition, whereby each fluid dynamic variable at a certain position and time 𝜙(𝑥,𝑡) is decomposed in a mean value 𝜙(𝑥,𝑡) and a fluctuating term 𝜙′(𝑥,𝑡). Figure 24.Illustration of a typical flow variable distribution in turbulent flow. In the heterogeneous and unsteady nature of turbulent flows, there coexist eddies of different scales that interact between them and also with the mean flow. The largest eddies have characteristic length and velocity in the same order of magnitude of the problem itself, let them be 𝑙0 and 𝑣0 respectively, and they are dominated by inertia. Therefore the characteristic time at this scale is 𝜏0= 𝑙0/𝑣0. Kolmogorov’s Theory (Pope, 2000) postulates that:  Large eddy velocities are in the order of magnitude of the turbulent intensity: 𝑣0∝𝑣′=(2𝑘/3)1/2 [Equation 2.1]  Large eddies dissipate their energy along their characteristic time, for which the dissipation rate is: 𝜀∝𝑣02/2 𝜏0 [Equation 2.2] Therefore the large eddy length scale, termed as the integral scale, is derived as: 𝑙0∝𝑘3/2 𝜀 [Equation 2.3] The Largest eddies interact with the mean flow, where the energy is produced and transferred to successively smaller scales until the energy is delivered to the smallest turbulence scales ~ 0.1-0.01 mm (Kolmogorov scale) (Versteeg & Malalasekera, 2007). Kolmogorov scales are dominated by viscosity, and consequently energy is finally transformed into thermal energy due to the action of viscous friction. Accordingly, kinetic energy flows from the energy containing to the dissipation range, passing through an intermediate range (Taylor scale) where inertial forces are predominant. In this energy flow; called energy scale (Figure 25), there must exist a balance so that the energy extracted by large eddies from the mean flow is dissipated in thermal form by the smallest ones. t Φ(x,t) Φ(x,t) Φ’(x,t) Φ(x,t) = Φ(x,t)+Φ’(x,t) Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 25 Figure 25. Energy cascade energy established in turbulent flows. Evaluating the Kolmogorov scale for a certain problem can be done by assuming that Reynolds number is approximately unity, since viscosity is dominant in this region. Kolmogorov’s velocity scale can be written as 𝑣𝜂=𝜈/𝜂, which together with the definition of dissipation rate, 𝜀𝜂=𝑣𝜂2/𝜏𝜂, and time scale, 𝜏𝜂=𝜂/𝑣𝜂, allow expressing the characteristic magnitudes for the Kolmogorov scale as functions of the fluid properties (𝜈=𝜇/𝜌) and the nature of the problem (𝜀) (Tennekes & Lumley, 1972). 𝜂=(𝜈3𝜀 ⁄)1/4 𝑣𝜂=(𝜀𝜈)1/4 𝜏𝜂=(𝜈𝜀 ⁄)1/2 [Equation 2.4] From which ratios between integral and Kolmogorov scales can be deduced in a certain problem. 𝜂𝑙0=𝑅𝑒𝐿−3/4 𝜏𝜂 𝜏0=𝑅𝑒𝐿−1/2 𝑣𝜂 𝑣0=𝑅𝑒𝐿−1/4 [Equation 2.5] The turbulent spectrum scales were also established by Kolmogorov’s theory and experimentally verified by several researchers (see figure 6.14 in (Pope, 2000)). The usual distribution of energy with respect to the eddy length scale (𝜆), expressed by means of the commonly employed wave number 𝜅= 2𝜋/𝜆, is represented in figure 26 (Pope, 2000), where it is possible to observe: MEAN FLOW L, T, V Integral scale l0, τ0, v0 Kolmogorov scale η, τη, vη k k Taylor scale Inertial range Dissipation range Energy containing range Production of energy Dissipation of energy Energy transfer to smaller eddies Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 26 Figure 26. Spectrum distribution of the turbulent energy.  The majority of the turbulent energy is contained in the largest eddies.  The integral length scale, 𝑙0, has the highest energy content.  In the Taylor scales, energy decreases linearly in a double logarithm system with a slope of -5/3 (Kolmogorov’s 5/3 Law).  There exists a length scale that separates inertial from dissipation ranges, below which the energy content decreases sharply. 2.3.1 Turbulent flow treatment in CFD Considering a temporal interval sufficiently longer than the longest time scale of the problem, it is possible to average the fluid motion equations giving place to the Reynolds-Averaged Navier-Stokes (RANS) equations, presented next:  Mass conservation: ∇·𝑈 󰇍 󰇍  =0 [Equation 2.6]  Momentum conservation: 𝜕𝑈 󰇍 󰇍  𝜕𝑡+∇·(𝑈 󰇍 󰇍  𝑈 󰇍 󰇍  )=−1𝜌∇𝑝 󰇍 󰇍 󰇍 󰇍 󰇍  +𝜈Δ𝑈 󰇍 󰇍  −∇·𝑣′ 󰇍 󰇍 󰇍  𝑣′ 󰇍 󰇍 󰇍  +f𝑚 [Equation 2.7]  Energy conservation: 𝜌𝑐(𝜕𝑇 𝜕𝑡+𝑈 󰇍 󰇍  ·∇T 󰇍 󰇍 󰇍 󰇍 󰇍  )=𝑘Δ𝑇+𝜙𝑈 󰇍 󰇍  +𝜙𝑣′ 󰇍 󰇍 󰇍 󰇍  −𝜌𝑐∇·𝑣′ 󰇍 󰇍 󰇍  𝑇′+𝑄 [Equation 2.8] This results in a system with 5 equations and 8 unknowns: velocities (𝑈 󰇍 󰇍  ), pressure (𝑝) and temperature (𝑇), the Reynolds stress tensor (𝜏𝑖𝑗=𝑣′ 󰇍 󰇍 󰇍  𝑣′ 󰇍 󰇍 󰇍  ), the turbulent heat transport (𝑣′ 󰇍 󰇍 󰇍  𝑇′) and the viscous dissipation due to turbulence (𝜙𝑣′ 󰇍 󰇍 󰇍 󰇍  ) (Zamora Parra, 2008). The concept of turbulent viscosity (𝜇𝑡) introduced by Boussinesq (Boussinesq, 1877) allows to close the problem. Turbulent viscosity relates the stresses in the flow with mean flow velocity gradients and is dependent upon the nature of the problem: 1.E-15 1.E-13 1.E-11 1.E-09 1.E-07 1.E-05 1.E-03 1.E-01 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 Spectral energy, E(m3/s2) Length scale, λ (m) l0 η Dissipation range Inertial range Length scale, λ (m) Spectral energy, E (m3/s2) Considerations about the CFD simulation of supercritical carbon dioxide in turbomachinery Design of supercritical carbon dioxide centrifugal compressors 27 𝜏𝑖𝑗=𝜇𝑡(𝜕𝑈𝑖 𝜕𝑥𝑗+𝜕𝑈𝑗 𝜕𝑥𝑖) [Equation 2.9] Although 𝜇𝑡 is a function of the mean flow velocities, the problem is not still closed. For that, a turbulence model needs to be implemented. In decreasing order of complexity, turbulence models can be classified into:  Direct Numerical Simulation (DNS) is the absence of turbulence model, meaning that Navier-Stokes equations are directly solved unsteadily and capturing the whole spectrum of scales. Very fine meshes are thus required in order to capture the flow phenomena in the Kolmogorov scale, which therefore poses a maximum cell size limit. The computational effort required by DNS is within the capacity of only some computers all over the world.  Large-Eddy Simulation (LES) does include turbulence models but only for the smaller scales while the interaction of larger eddies is simulated. This is in consequence inherently unsteady and three-dimensional. LES was developed by Smagorinsky in 1963 (Smagorinsky, 1963) for studying the dynamics of general circulation in the atmosphere. In 1970, Deardorff (Deardorff, 1970) applied this proposal to Pouseuille flow at high Reynolds number, where wall-affected flows were analysed by advanced turbulence models. After that, LES was profusely extended to other industrial applications because of the good agreement to reality and the lower computational cost in comparison with DNS. LES philosophy is based on a low-pass filter that is applied to fluid motion equation so that Sub- Grid Scale (SGS) models are applied for scales smaller than the cell sizes, whereas above that the equations are solved directly. Although lower than DNS, the computational cost of performing LES is important, even more for S-CO2 as it is characterised by a Reynolds number higher than air. Based on this, it is clear that LES is out of the scope of this project.  RANS equations, as explained before, result from time averaging the Navier-Stokes equations. Thus, steady simulations can be performed easily and simplifications such as symmetric or periodic conditions can be applied. The system composed by equations 2.6- 2.8 has 2 degrees of freedom and therefore RANS models require the implementation of additional turbulence models with one (Spallart-Allmaras (Spalart & Allmaras, 1992)), two (𝑘−𝜀 (Launder & Spalding, 1972) and 𝑘−𝜔 (Wilcox, 1994)) or even more equations. Theoretically, if DNS and LES results were time-averaged along long enough intervals of time they would coincide with results of the RANS model. 2.3.2 Turbulence model applied to S-CO2 turbomachinery Taking into account that no doubts exist regarding the EOS, the choice of a suitable turbulence model has the largest impact on the results, apart from numerical aspects such as discretisation, interpolation schemes and other solver features. Given the lack of experimental results in the field (only some results of the centrifugal compressor at SANDIA NL publically available (Wright, et al., September 2010)), an examination of the various turbulence models available in ANSYS Fluent® has been performed in order to help make the choice. It is also true that other numerical simulations have been carried out in the field of supercritical carbon dioxide turbomachinery ( (Pecnik, et al., 2012), (Rinaldi, et al., 2013) and (Munroe, et al., 2009)) and heat exchangers (Van Abel, et al., 2011), which offer a sort of state-of-the-art in the field. These works confirm that two-equation models are the most common choice since, as stated by Wilcox “They are, in fact, the simplest complete model of turbulence” (Wilcox, 1994). Within the usual constraints in computational capacity, whether use 𝑘−𝜀 or 𝑘−𝜔 models is the most frequent question to answer. The 𝑘−𝜔 model was initially proposed by Wilcox (Wilcox, 1988), but it was not till its Shear-Stress-Transport (SST) version was created by Menter (Menter, 1994), that its use was profusely extended across the industrial world (Menter, et al., 2003). Older than 𝑘−𝜔, 𝑘−𝜀 models, originally introduced by Launder & Spalding (Launder & Spalding, 1972) are the alternative. These are especially suitable for fully turbulent (high Reynolds number) non-separated flows. However, its main Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 34 CHAPTER 3. Analysis of S-CO2 flows in elementary diffusion processes Turbocompressors are sophisticated machines that increase the total pressure of a certain inlet flow of gas by converting mechanical work at the rotor shaft into total enthalpy of the fluid at the exit. This energy conversion is a two-step process taking place in the rotor and stator rows of the machine. Thus, the rotor converts shaft work into kinetic energy (expressed in terms of absolute velocity) and enthalpy (i.e. pressure), delivering a flow at high speed to the static blade row where the gas decelerates and kinetic energy is converted into enthalpy giving place to a flow at high pressure and low/moderate velocity. The transformation of kinetic energy into pressure is called diffusion, after which the static blade row of a radial compressor is called the diffuser, and in subsonic adiabatic flow it requires increasing cross sectional area ducts. Based on this rationale, the intersection between a stream surface in an axial compressor and the channels enclosed within two adjacent blades presents a two-dimensional divergent cross sectional area pitch-wise. However, this geometry must not be mistaken for the reduction in channel height spanwise which accommodates to the increasing density of the gas at constant mass flow rate. In an axial compressor, the higher the outlet to inlet area ratio, the higher the static pressure rise. This is all presented in the figure below (Bölcs & Tsamourtzis, 1991). Figure 31. Details of turbomachinery geometry (Bölcs & Tsamourtzis, 1991). The corresponding hub (bottom) and shroud (top) surfaces must increase their cross section area, both for rotor and stator, especially in axial machines, where the increment of fluid total pressure through the rotor is much lower than in radial machines. In the latters, centrifugal force allows for achieving higher rotor pressure ratio, bringing about marked density changes and so determining smaller outlet section areas. Nevertheless, independently on machine type (axial or radial) and blade row (rotor or stator), the relative flow to the blade is always diffused, accordingly decreasing relative velocity form inlet to outlet. The afore discussed fundamentals apply both to axial and radial compressors. Nevertheless, in the latter, shaft work is converted into enthalpy (i.e. pressure) by virtue of the so called centrifugal forces. In effect, one feature that sets radial machines apart from axial ones is the noteworthy variation of blade speed. This velocity, which is assumed constant in virtually all the stages in an axial turbomachinery (in particular those stages with a hub-to-tip ratio higher than 0.7-0.75), almost triples from the inlet to the outlet of a centrifugal compressor impeller, thus allowing for a much larger energy transfer from the shaft to the fluid. This differential feature will eventually lead to a much higher stage pressure ratio. Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 35 Nevertheless, prior to addressing the whole centrifugal compressor problem, it is deemed convenient to begin with an analysis of elementary diffusers, since diffusion plays a crucial role in the turbomachinery. This is further confirmed by Japikse and Baines (Japikse & Baines, 1998) who pointed out that “a change of 0.01 in pressure recovery can be equivalent to a few tenths of a point of stage efficiency”. Conical diffusers involve simpler physical phenomena in comparison with the complex fluid dynamics of a full compression stage, what can help to comprehend how supercritical carbon dioxide flows behave in diffusion processes, which differences exhibit with respect to a conventional fluid (for instance air) and, in the light of these results, how the standard design rules can be customised in order to get highly efficient machines with the new working fluid. In this regard, an extensive study on conical diffusers is carried out in this chapter with the endeavour of estimating the impact of the main dimensionless parameters of influence in accordance with the available literature. 3.1 Definition and scope of the analysis The analysis of diffusers was initiated by Gibson in 1910 (Gibson, 1910), who basically analysed the influence of some geometrical parameters on the performance of conical diffusers working with water for hydraulic applications. This work gave place to an intense activity in this topic between the 50’s and 80’s most of which considered air as the working fluid. Unexpectedly, some of the parameters that had been disregarded previously proved to have a noteworthy effect on system performance whereas others which were presumed to be dominant showed up to hardly have any influence, for instance Mach number as demonstrated in (Dolan & Runstadler, 1973). As a consequence of this research, a vast amount of information was generated, mainly in the form of performance maps that resulted very useful for turbomachinery design (Japikse, 1996) and contributed significantly to the development of gas turbines and aero engines. In spite of a huge database of air diffuser performance being currently available, there is very little information of such devices operating on other less common fluids. This raises the following questions: is this database applicable to supercritical carbon dioxide diffusers? How large is the impact of the real gas behaviour on the performance of such diffuser? These are the two main questions that this chapter aims to answer. For this purpose, the parameters of influence are firstly presented in order to define the boundaries of the study. As a second step, the numerical tool on which the work bases is presented and validated against the available information for air diffusers. Finally, the results are presented for the different parameters considered and a number of conclusions about how to design efficient turbomachinery for S-CO2 are drawn. 3.1.1 Diffuser operation In spite of its usually simple geometry, the operation of a diffuser can become very complex due to the irreversibilities caused by diverging walls. As opposed to nozzles, where converging walls tend to practically remove the outflow boundary layer, in diffusers even isentropic flows at the inlet can become very distorted and unsteady at the outlet depending on geometry and inlet flow conditions. Several works done during the 50’s and 60’s focused on studying different operation regimes of diffusers and on the impact of flow features on performance: for instance inlet turbulence (Fox & Kline, 1962) (Moore & Kline, 1958) and aerodynamic blockage (Runstadler & Dolan, 1973) (Reneau, et al., 1967). Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 36 Figure 32. Operating regions of a straight line diffuser 6 . (Smith, 1975) Even though figure 4 refers to channel diffusers, it is useful to illustrate all the possible regimes in which any diffuser can operate; this includes conical ones. The desired performance is no appreciable stall and it is achieved for moderate divergence angles (2𝜃) and non-dimensional lengths (𝐿/𝑤1). On the contrary, jet flow appears when divergence is very high. In between these, two different stall phenomena are found. First, when divergence angle is increased starting from no stall operation, there is a certain and reproducible point where transient separation takes place (line AA). Bubbles are created locally on the diffuser walls, generating large and appreciable eddies which make the flow transient. For a given nondimensional length, there is a certain and also reproducible divergence angle for which maximum transitory effects are noticed. If the divergence angle is further increased, a second transition to a new steady operation is observed (line BB) which is characterised by a two-dimensional steady flow completely detached from the wall. A visible stream tube travels through the diffuser with a nearly constant cross sectional area while a stagnant fluid fills the rest of the diffuser. Figure 33. Diffuser operating in two-dimensional stall region. When increasing divergence angle, the separation between two-dimensional stall and jet flow is represented by line CC in figure 32figure 4. However, if divergence angle is decreased from jet flow operation, transition to two-dimensional stall doesn’t happen in the same line, but for a smaller angle in line DD, giving place to a hysteresis phenomenon. 6 The geometrical parameters participating in this figure are: the divergence angle (2θ) and the length to inlet width ratio (L/w1). Linked to these two parameters, it could be used the area ratio, which is defined for straight line diffusers by AR=1+ (L/w1)tanθ. 1.5 2 3 4 6 8 10 15 20 30 40 60 80 100 2θ 1.5 12 3 4 6 8 10 15 20 30 40 60 L/w1 NO APPRECIABLE STALL LARGE TRANSITORY STALL LINE OF MAXIMUM UNSTEADINESS INCREASING UNSTEADINESS FULLY DEVELOPED TWO-DIMENSIONAL STALL HYSTERSIS ZONEJET FLOW A A B C D C D B Stall, stagnant fluid Nearly constant area flow Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 37 3.1.2 Parameters of influence in conical diffuser performance Dimensional analyses provide useful information about which and how many parameters have an effect on the behaviour of a certain physical problem/system, the extent of which is then revealed by experiments. For diffusers, this has regularly been done over the years by a number of researchers. A complete review of this is presented in (Japikse & Baines, 1998), from where the following list of parameters is derived: Figure 34. Diffuser geometry.  Non-dimensional length, which is defined as the ratio from the length to the throat diameter of a conical diffuser.  Divergence angle, i.e. 2𝜃, or equivalently area ratio between outlet and inlet. These geometric parameters are related as shown by [Equation 3.1] for conical diffusers. 𝐴𝑅=[1+2(𝐿/𝐷𝑡ℎ)tan𝜃]2 [Equation 3.1]  Aerodynamic blockage at the throat, defined as the fraction of cross sectional area that is occupied by the boundary layer. It can also be related with the ratio between actual and ideal mass flow rates and is directly linked to the boundary layer displacement thickness. 𝐵𝑡ℎ=1−𝐴𝑒𝑓𝑓 𝐴=1−𝑚󰇗 𝑚󰇗𝑖𝑑=4𝛿∗ 𝐷𝑡ℎ(1−𝛿∗ 𝐷𝑡ℎ) [Equation 3.2] Amongst all fluid dynamic parameters subsequently described, aerodynamic blockage is the most relevant one. Boundary layer displacement thickness was already investigated in the 60’s by Stratford and Tubbs (Stratford & Tubbs, 1965), who found that it constrains the maximum pressure recovery in a diffuser (first with theoretical models and afterwards by comparison against experiments). Nearly a decade later, Dolan and Runstadler (Dolan & Runstadler, 1973) completed a detailed study on the effect of blockage and others aerodynamic parameters on diffuser performance. This is considered a seminal work in the topic and will be used here for validation purposes.  Both Mach and Reynolds numbers are commonly taken into account in any fluid dynamic problem. However, although in the early stages of diffuser research they were considered critical, it was later confirmed experimentally that they have a rather soft impacts. As a matter of fact, virtually no influence of Reynolds number was observed beyond 105 in (Dolan & Runstadler, 1973), confirming the existence of a critical value of Re around that value, beyond which Reynolds number does not affect.  Apart from that blockage, other features of the inlet velocity distribution have been studied in literature: swirl and distortion. Swirl refers to the angular velocity and is usually proportional to the radius. It modifies the development of both core and boundary layer flows through the diffuser, largely improving its performance as shown in figure 35. On the other hand, distortion affects to the axial component and it manifests as a momentum L Dth Dout 2θ (1) (2) Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 38 defect/excess at certain locations of the diffuser throat. Three main studies were done in this field: (Srinivasan & Lakshmi Narasimhan, 1978 ) , (Sajben, et al., 1974) and (Shimizu, et al., 1982), where performance improvements were sometimes observed, even though not as noticeable as for swirl. Figure 35. Comparison of conical diffuser performance without (a) and with (b) swirl effect found by (Mc Donald, et al., 1971)  The influence of turbulence in conical diffuser performance can also be analysed by comparing the works available in literature even though it needs to be extrapolated from the results for other geometries. In this regard, the systematic study done by Hoffman for channel diffusers ( Hoffman, 1981) is worth mentioning, figure 36. Figure 36. Turbulence intensity and turbulence length scale effects in channel diffusers ( Hoffman, 1981) As noticed, there is a critical value of the non-dimensional turbulence length scale beyond which this parameter has no effect (~8 in figure 36). And similarly, this figure suggests the existence of a critical turbulence intensity after which its effects disappear. Additionally, Hoffman pointed out that in addition to intensity and length scale having an impact on diffuser performance, the turbulence structure, i.e. whether it is isotropic or not, influences it.  The real behaviour of CO2 near the critical point suggests that the influence of the compressibility factor (𝑍=𝑝𝑣/𝑅𝑇) be studied. It must be highlighted that this is an 1.0 0.8 0.6 0.4 0.2 00.5 2.0 4.0 8.0 (a) PERFORMANCE COEFFICIENT , Cp Symbol Divergence Angle, 2θ 4.0 º 6.0 º 8.0 º 12.0 º 15.8 º 31.2 º AR - 1 IDEAL PERFORMANCE Cp,i 1.0 0.8 0.6 0.4 0.2 00.5 2.0 4.0 8.0 (b) PERFORMANCE COEFFICIENT , Cp Symbol Divergence Angle, 2θ 4.0 º 6.0 º 8.0 º 12.0 º 15.8 º 31.2 º AR - 1 IDEAL PERFORMANCE Cp,i 0.583 (0.713) 0.584 (0.709) 0.658 (0.751) 0.631 (0.736) Cp = 0.620 (Cp = 0.730) Cp = 0.680 (Cp = 0.770) 0.703 (0.798) 0.674 (0.785) 0.717 (0.795) Cp AT 2θ = 20º (Cp AT 2θ = 12º) ROD AXIS IN Z DIRECTION : :ROD AXIS IN Y DIRECTION 0.678 (0.773) 0.687 (0.782) 0.616 (0.724) 0.682 (0.778) 0.716 (0.788) Cp = 0.718 (Cp = 0.786) 0.702 (0.792) 0.713 (0.793) 0.835 (0.737) 0.888 (0.754) 0.898 (0.773) 46 8 10 12 2 3 4 5 6 7 λg / δ* um(%) v’2/ Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 39 additional parameter never considered before because of the ideal behaviour of most working fluids in turbomachinery applications (𝑍=1). In contrast, the fluid of interest in this works exhibits a very variable compressibility factor depending on the static pressure and temperature across the compressor.  The parameters shown up to now are the most relevant ones because of the large impact that they have on diffuser performance as well as their variability from one diffuser application to another. Apart from them, some authors studied different aspects related to some specific diffuser applications such as the rotation of the diffuser, oscillatory inflows, etc… This effects that will not be considered in this work. To sum up, after dimensional analysis, diffuser performance can be expressed as a function of 11 parameters: area ratio (𝐴𝑅), non-dimensional length (𝐿/𝐷𝑡ℎ), aerodynamic blockage (𝐵𝑡ℎ), Mach (𝑀𝑡ℎ) and Reynolds numbers (𝑅𝑒𝑡ℎ), dimensionless parameters for swirl (𝐶𝑠𝑤𝑖𝑟𝑙) and distortion (𝐶𝑑𝑖𝑠𝑡), turbulence intensity (𝑇𝑢𝑡ℎ), non-dimensional turbulence length scale (𝑙𝑡ℎ/𝛿∗) and a parameter identifying the turbulence structure (𝐶𝑡,𝑒𝑠𝑡), i.e. its anisotropy. The compressibility factor (𝑍𝑡ℎ) is also added. 𝐶𝑝=𝑓(𝐴𝑅,𝐿/𝐷𝑡ℎ,𝐵𝑡ℎ,𝑀𝑡ℎ,𝑅𝑒𝑡ℎ,𝐶𝑠𝑤𝑖𝑟𝑙,𝐶𝑑𝑖𝑠𝑡,𝑇𝑢𝑡ℎ,𝑙𝑡ℎ/𝛿∗,𝐶𝑡,𝑒𝑠𝑡,𝑍𝑡ℎ) [Equation 3.3] Given that the experiments are carried out by numerical simulation in ANSYS FLUENT®, and owing to the limited computational capacity at the disposal of this work, no analyses were done in terms of turbulence structure. In other words, the RANS approach is compulsory and this does not allow this type of analysis which would require Large Eddy Simulation (LES). Hence only the influence of the remaining parameters of influence was analysed. 3.1.3 Presentation of results and range of study Once the key parameters have been exposed, the outputs and the format in which they are presented are described in this section. The most common diffuser performance indicator is the static pressure rise coefficient (𝐶𝑝), defined as the ratio from static pressure rise to inlet dynamic head. This coefficient does not provide complete information about how efficiently the diffuser operates though. In effect, it is possible to achieve a high value of 𝐶𝑝 in a rather irreversible and thus inefficient process (i.e. prohibitive total pressure losses) as long as diffuser geometry allows that level of pressure recovery. Therefore, the diffuser effectiveness (𝜂𝐷), defined as the ratio between the actual and ideal values of 𝐶𝑝 is needed to fully assess the optimality of the device, where the ideal is 𝐶𝑝 obtained in isentropic conditions and only depends on area ratio: 𝐶𝑝,𝑖=1−𝐴𝑅−2. An additional performance parameter considered in this work, which is equivalent to the previous one, is the total pressure loss coefficient (𝐾𝐷), defined as the ratio between total pressure losses and inlet dynamic head. 𝐶𝑝=𝑝2−𝑝1 𝑝01−𝑝1 [Equation 3.4] 𝜂𝐷=𝐶𝑝/𝐶𝑝,𝑖 [Equation 3.5] 𝐾𝐷=𝑝01−𝑝02 𝑝01−𝑝1 [Equation 3.6] These parameters were measured for different fluid dynamic conditions in 30 diffuser geometries defined by red points in figure 37: 𝐴𝑅 = 3 - 8, 𝐿/𝐷𝑡ℎ = 6 - 20 and 2𝜃 = 4 - 10º. The effect of geometry was plotted by curves of constant 𝐶𝑝 in the 𝐿/𝐷𝑡ℎ- 𝐴𝑅 plane. Coordinates were selected to transform constant angle curves into lines, hence a double-logarithmic scale was chosen to represent the effect of geometry: log(𝐿/𝐷𝑡ℎ) in abscissa and log(√𝐴𝑅−1) in ordinate [equation 3.1]. Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 40 Figure 37. Typical performance map of a diffuser, showing the effect of geometry. This is the standard map employed to represent the effects of geometry along with the so called first order parameters: 𝐵𝑡ℎ, 𝑀𝑡ℎ and 𝑅𝑒𝑡ℎ. The influence of second order parameters, which are those related to velocity profile and turbulence, are available in closely related MSc thesis (López Florenciano, 2013) developed by fellow student at the GMTS - Department of Energy Engineering under the supervision of Prof. D. Sánchez. The main results of this complementary work are summarised in section 3.4. 3.1.4 The tool: definition, validation and verification The heart of the tool employed for the analysis is a two-dimensional model created in ANSYS FLUENT® which comprises the k-ω SST model for turbulence, the Span and Wagner equation of state for supercritical CO2 (Span & Wagner, 1996) and a mesh that accomplishes the required refinement level with 𝑀𝐴𝑋(𝑦+)≤1. This model acts as a slave and runs through journal files 7 (*.jou), that are previously generated by a Matlab script, which is the master block commanding the previous one (see figure 38). The main loop controls every variable by modifying a certain input parameter of the model, therefore the following relations between manipulating and control variables can be established: 7 A journal file contains a sequence of ANSYS FLUENT® commands, arranged as they would be typed interactively in the program (ANSYS, 2011). 6 8 10 12 14 16 18 20 2 3 4 5 6 7 8 9 L/Dth AR 2θ = 4º 2θ = 5º 2θ = 6º 2θ = 7º 2θ = 8º 2θ = 9º 2θ = 10º Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 41 Control variable Manipulating variable Aerodynamic blockage at throat Parameter 𝛿 in the inlet flux profile Mach number at throat Static pressure at outlet boundary Reynolds number at throat Throat diameter Swirl coefficient Angular velocity profile Distortion parameter Type of inlet flux profile Turbulence intensity Turbulence intensity at inlet boundary Dimensionless turbulence length scale Turbulence length scale at inlet boundary Compressibility factor Stagnation temperature at inlet boundary Wall 𝑦+ Size of first element in prismatic layer Table 3. Control and manipulating variables of the tool. 3.1.4.1 Simplified model (in EES) For given influence parameters, a simplified model developed in EES calculates the geometry of the diffuser and the isentropic flow conditions at the inlet and outlet. Geometric results are used for geometry and mesh generation with ICEM CFD® while isentropic flow conditions are useful for initialising both the flow and boundary conditions (see section 3.1.4.4), some of which are manipulating variables that need to be changed along the calculation loop (Figure 38) in order to keep the control variables within the desired ranges. As shown below, the model has been built by imposing mass and energy conservation along with no entropy generation. ℎ0=ℎ+𝑣2/2 𝑚󰇗=𝜌𝑣𝜋4𝐷2 𝑠0=𝑠 [Equation 3.7] The Span and Wagner equation of state for CO2 closes the model when applied for static and stagnation conditions both at inlet and outlet together with the definitions of Mach and Reynolds numbers. 𝑀=𝑣/𝑆𝑆 𝑅𝑒=𝜌𝑣𝐷 𝜇 [Equation 3.8] Where the properties of interest are: enthalpy (ℎ), entropy (𝑠), density (𝜌), speed of sound (𝑆𝑆) and viscosity (𝜇). Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 42 Figure 38. Block diagram of the numerical tool. Given the geometrical parameters and fluid dynamic conditions at the throat, it is possible to calculate the dimensions of the diffuser as well as the thermodynamic conditions (pressures and temperatures) at inlet and outlet. Therefore, geometry and mesh can be constructed. Simulations require defining initial values, which are taken from the outputs of the simplified model, and boundary conditions, some of them specified from the beginning and others obtained from the simplified model. Once the simulation is run, wall y+ is checked to be lower than one. If this is not accomplised, a mesh adaption is performed and the simulation launched again. If wall y+ <1, the throat Mach number is adjusted by modifying the static pressure outlet (the lower this pressure, the higher the Mach number). Then, if the Reynolds number obtained is far from the target value, the problem is scaled up/down in order to modify the throat diameter. Finally, once all parameters have converged, Fluent is open for last time to read the values of interest. Fluent END Geometry Isentropic flow Mesh Initial values Boundary conditions Mesh Generation Simplified model EES ICEM udf generation MAX(y+) ≤ 1? Desired Mth? Desired Reth? Fluent Mesh Adaption Pout modification Mesh Scaling Fluent Fluent Matlab CFD Simulation Inlet flux profile Tuth, lth/δ* T0th, Tout, Pout NO NO NO YES YES YES Bth, Cswirl, Cdist vth, ρth INPUTS: AR, L/Dth, Bth, Mth, Reth, Cswirl, Cdist, Tuth, lth/δ*, Zth AR, L/Dth, Mth, Reth, Zth Read case Fluent Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 43 3.1.4.2 Flux profile adjustment As presented in section 3.1.4.4, a mass flow inlet boundary is imposed at the diffuser throat, for which a flux profile (mass flow rate per unit area, 𝜌𝑣) can be specified in lieu of a total mass flow rate. This allows changing the velocity profile thus playing with blockage, swirl and distortion independently. Though some variations were done for distortion analysis (section 3.4.1), the basic flux profile for blockage modifications is composed by a constant value in the core flow and a parabola in a virtual boundary layer 8 with zero flux in the wall and peak flux in the vertex (Figure 39). Figure 39. Flux profile at diffuser throat, imposed as boundary condition. Density (𝜌), core velocity (𝑈∞) and throat radius (𝑟𝑡ℎ) are known when the influence parameters are given, so the following conditions are required to determine the parabola constants (𝑎, 𝑏 and 𝑐) and the virtual thickness 𝛿𝑣.  Null flux on the wall because of the non-slip conditions fixed in the CFD model. 𝑎·𝑟𝑡ℎ 2+𝑏·𝑟𝑡ℎ+𝑐=0 [Equation 3.9]  Continuity of the flux function at the virtual boundary layer edge. 𝑎·(𝑟𝑡ℎ−𝛿𝑣)2+𝑏·(𝑟𝑡ℎ−𝛿𝑣)+𝑐=𝜌·𝑈∞ [Equation 3.10]  Maximum flux at the virtual boundary layer edge, [𝑑𝑓/𝑑𝑟]𝑟=𝑟𝑡ℎ−𝛿𝑣=0, which results in: 2𝑎·(𝑟𝑡ℎ−𝛿𝑣)2+𝑏=0 [Equation 3.11]  Specified value of the virtual aerodynamic blockage, which is directly related to the virtual thickness. 𝛿𝑣=𝑟𝑡ℎ(1−√1−𝐵𝑡ℎ,𝑣) [Equation 3.12] 8 The term “virtual” is used herein because δ will not coincide with the displacement boundary layer thickness obtained after simulation, as explained in section 3.1.4 later. flux (kg/s·m2) r (m) f = ρU∞ rth rth - δv 0 f = ar2 + br + c Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 50 obtained for a 1 % turbulence intensity and non-dimensional turbulence length scale equal to 10. Blockage and Mach number are set to the same values for air and supercritical CO2 whereas, on the contrary, Reynolds number differs from one fluid to another due to the very different fluid properties, density in particular (more than 400 times higher for CO2). Nevertheless, as stated in (Dolan & Runstadler, 1973), the effect of 𝑅𝑒 on diffuser performance is negligible in the range of study. The maps in figure 46 comprise curves for S-CO2 (simulation) in red and air (experiments) in green. The particular (fluid dynamic) boundary conditions appear at the top of each single map for CO2 and on a green square at the bottom for air. Both 𝐶𝑝 surfaces intersect in the dashed purple line, which divides every map into two regions: (i) one where air has better performance (shadowed in green) and (ii) one where supercritical carbon dioxide performs better (shadowed in red). Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 51 Figure 46. Comparative performance maps for CO2 and air. Results for CO2 are numerical whilst air data are taken from (Dolan & Runstadler, 1973). 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.68 0.72 0.76 0.8 0.82 0.84 0.86 0.88 0.78 0.76 0.74 0.72 0.70 Bth = 3 % ; Mth = 0.2 ; Reth = 1.2·105 Air + S-CO2 + 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.620.66 0.7 0.74 0.78 0.8 0.8 0.82 0.84 Bth = 6 % ; Mth = 0.2 ; Reth = 1.2·105 0.72 0.70 0.68 0.660.64 0.62 S-CO2 + Air + 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.54 0.6 0.64 0.68 0.71 0.74 0.76 0.76 0.78 0.78 Bth = 9 % ; Mth = 0.2 ; Reth = 1.2·105 0.68 0.66 0.640.60 0.58 0.56 0.54 Air + S-CO2 + 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.5 0.56 0.6 0.64 0.67 0.7 0.72 0.72 0.74 0.74 Bth = 12 % ; Mth = 0.2 ; Reth = 1.2·105 0.64 0.62 0.60 0.580.56 0.54 0.52 S-CO2 + Air + 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.680.720.75 0.78 0.81 0.81 0.84 0.86 0.88 Bth = 3 % ; Mth = 0.4 ; Reth = 2.27·105 0.78 0.76 0.74 0.72 0.70 Air + S-CO2 + 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.60.64 0.68 0.72 0.75 0.75 0.78 0.78 0.8 0.8 0.82 0.70 0.68 0.66 0.64 0.62 Bth = 6 % ; Mth = 0.4 ; Reth = 2.27·105 Air + S-CO2 + 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.50.540.58 0.62 0.66 0.7 0.7 0.73 0.73 0.75 0.75 0.66 0.64 0.62 0.60 0.58 0.56 Bth = 9 % ; Mth = 0.4 ; Reth = 2.27·105 S-CO2 + Air + 0.68 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.45 0.5 0.54 0.58 0.62 0.62 0.66 0.66 0.69 0.69 0.72 0.72 0.64 0.62 0.60 0.58 0.560.54 0.52 Bth = 12 % ; Mth = 0.4 ; Reth = 2.27·105 S-CO2 + Air + Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 52 By observing the previous maps, the first conclusion is that CO2 performs better than air in a wider region of the range considered, this region being larger when blockage decreases. The shapes of the curves are also different for each fluid. The green lines present a wide region where the area ratio hardly affects the performance (constant 𝐶𝑝 are almost vertical) while this zone does not even exist for CO2in many cases. This means that for certain combinations of operating conditions, the cone diffuser operating with air seems to be more resistant to high divergence angles and thus to flow detachment. In other words, high 𝐶𝑝’s can be sought by increasing both the area ratio and the non-dimensional length. On the contrary, S-CO2 is less resistant to flow detachment and hence needs an increment of dimensionless length at limited area ratio in order to improve the performance. The previous comparison shows experimental data for air and numerical results for carbon dioxide. Apart from this comparison, simulations have also been carried out with air as working fluid (see maps in figure 49), whose validation is shown in section 3.1.4.5. For these simulations, similitude regarding turbulence parameters and equal blockage factor and Mach number are considered and, indeed, it comes to show the improvement in performance when CO2 is used as well as the fact that varying the area ratio can be favourable for air but rarely for carbon dioxide. Additional information about this effect will later be shown in section 3.3.1. Figure 47. Comparison between air (blue) and S-CO2 (red): effect of area ratio onto diffuser performance at constant non-dimensional length, 𝑀𝑡ℎ= 0.2 and 𝐵𝑡ℎ= 3 %. This result evidences that the supercritical fluid adapts better to longer diffusers with lower divergence angles, which was actually to be expected given the much higher Reynolds number characteristic of S-CO2. As latter explained in section 3.3.1, the dominant effect of inertia over friction incurs a lower capability of the flow to attach to the wall and hence the flow is very sensitive to an increase in cross sectional area. In order to further explain the observed behaviour, a comparison between air and S-CO2 flows at 𝑀𝑡ℎ=0.2 and 𝐵𝑡ℎ=6 % are shown for three different geometries: cases 3, 20 and 25. Cases 3 and 20 have the same non-dimensional lengths, cases 20 and 25 have the same area ratios and cases 3 and 25 have similar divergence angles. From the comparison of diffusers with the same non-dimensional length (Figure 48Figure 61 top), it is concluded that the wider the divergence angle, the faster the boundary layer growth in the entry region (near the throat) for S-CO2. This evidences that air adapts better to area changes. However, when diffusers with same area ratios but different lengths are compared (Figure 48 centre), it is observed that the aerodynamic blockage tends to converge to a common value as the flow develops along the diffuser. Furthermore, the blockage factor is even lower for S-CO2 if the diffuser is long enough, as observed for case 25. Therefore, near the throat, carbon dioxide separates from the walls because of inertia but, as velocity decreases and the significance of viscosity forces increases, the velocity profiles become more uniform in such a way that the boundary layer at the outlet is thinner for the supercritical fluid if the diffuser is long enough. Finally, a comparison of diffusers with similar divergence angles (Figure 48 bottom) shows that the aerodynamic blockage develops similarly streamwise. When the length is increased, the region where viscous forces are more important is proportionally longer and, in consequence, the velocity profiles are more homogeneous and the blockage factor is reduced. On the negative side, this increases the effect of wall friction though this is attenuated for the supercritical fluid thanks to the higher Reynolds number of CO2 with respect to air. This explains why high non-dimensional lengths benefit CO2 in terms of diffuser performance. 0.65 0.70 0.75 0.80 0.85 0.90 0 2 4 6 8 10 Cp AR S-CO2 vs Air , M=0.2 , Bth=3% L/D=6 L/D=8 L/D=10 L/D=12 L/D=14 L/D=16 L/D=18 L/D=20 L/D=6 L/D=8 L/D=10 L/D=12 L/D=14 L/D=16 L/D=18 L/D=20 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 53 Figure 48. Development of blockage (dashed lines) and velocity profile (solid symbols) along various diffusers for air (blue) and S-CO2 (red). Constant non-dimensional length (top), constant area ratio (middle) and similar divergence angle (bottom). Figure 49. Air diffuser maps obtained numerically. A last set of figures is now presented to confirm the beneficial effect of non-dimensional length on the performance of supercritical carbon dioxide diffusers. In effect, even though blockage factor somehow transforms the shape of constant-𝐿/𝐷𝑡ℎ curves, constant-𝐴𝑅 curves in the figures below always exhibit the same trends, definitely confirming that increasing the non-dimensional length brings about better diffuser performance. 0 10 20 30 40 50 60 70 0 5 10 15 20 010 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 Blockage (%) Radius (mm) Length (mm) ( ) Case 20: AR = 7, L/Dth = 10 ( ) Case 3: AR = 3, L/Dth = 10 (blue) AIR (red) S-CO2 0 10 20 30 40 50 60 70 0 5 10 15 20 020 40 60 80 100 120 140 160 180 200 220 240 260 Blockage (%) Radius (mm) Length (mm) ( ) Case 20: AR = 7, L/Dth = 10 ( ) Case 25: AR = 7, L/Dth = 20 (blue) AIR (red) S-CO2 0 5 10 15 20 25 30 35 0 5 10 15 20 020 40 60 80 100 120 140 160 180 200 220 240 260 Blockage (%) Radius (mm) Length (mm) ( ) Case 25: AR = 7, L/Dth = 20 (blue) AIR (red) S-CO2 ( ) Case 3: AR = 3, L/Dth = 10 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for AIR Bth = 3% ; Mth = 0.2 ; Reth = 6.0·104 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.68 0.7 0.72 0.74 0.74 0.76 0.76 0.78 0.8 0.82 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for AIR Bth = 6% ; Mth = 0.2 ; Reth = 6.0·104 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.60.64 0.680.7 0.72 0.740.76 0.78 0.8 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for AIR Bth = 9% ; Mth = 0.2 ; Reth = 6.0·104 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.54 0.6 0.64 0.680.7 0.72 0.720.74 0.76 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for AIR Bth = 12% ; Mth = 0.2 ; Reth = 6.0·104 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.5 0.54 0.58 0.62 0.66 0.68 0.68 0.7 0.72 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 54 Figure 50. Pressure rise coefficient vs. area ratio (left) and non-dimensional length (right) for S- CO2 at 𝐵𝑡ℎ= 3 %, 𝑀𝑡ℎ= 0.2, 𝑅𝑒𝑡ℎ= 7.4·106. Figure 51. Pressure rise coefficient vs. area ratio (left) and non-dimensional length (right) for S- CO2 at 𝐵𝑡ℎ= 3 %, 𝑀𝑡ℎ= 0.4, 𝑅𝑒𝑡ℎ= 1.6·107. Figure 52. Pressure rise coefficient vs. area ratio (left) and non-dimensional length (right) for S- CO2 at 𝐵𝑡ℎ= 15 %, 𝑀𝑡ℎ= 0.2, 𝑅𝑒𝑡ℎ= 1.6·107. In order to complete this section on the effect of geometry, all the maps obtained numerically for carbon dioxide are presented below. 0.65 0.70 0.75 0.80 0.85 0.90 2 4 6 8 Cp AR M=0.2, Re=7.4e6, Bth=3% L/D=6 L/D=8 L/D=10 L/D=12 L/D=14 L/D=16 L/D=18 L/D=20 0.65 0.70 0.75 0.80 0.85 0.90 4 6 8 10 12 14 16 18 20 22 Cp L/D M=0.2, Re=7.4e6, Bth=3% AR=3 AR=4 AR=5 AR=6 AR=7 AR=8 0.65 0.70 0.75 0.80 0.85 0.90 2 4 6 8 Cp AR M=0.4, Re=1.6e7, Bth=3% L/D=6 L/D=8 L/D=10 L/D=12 L/D=14 L/D=16 L/D=18 L/D=20 0.65 0.70 0.75 0.80 0.85 0.90 4 6 8 10 12 14 16 18 20 22 Cp L/D M=0.4, Re=1.6e7, Bth=3% AR=3 AR=4 AR=5 AR=6 AR=7 AR=8 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 2 4 6 8 Cp AR M=0.2, Re=1.6e7, Bth=15% L/D=6 L/D=8 L/D=10 L/D=12 L/D=14 L/D=16 L/D=18 L/D=20 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 4 6 8 10 12 14 16 18 20 22 Cp L/D M=0.2, Re=1.6e7, Bth=15% AR=3 AR=4 AR=5 AR=6 AR=7 AR=8 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 55 Figure 53. Diffuser maps for S-CO2 at 𝑀𝑡ℎ= 0.2, 𝑅𝑒𝑡ℎ= 7.4·106, 𝑇𝑢𝑡ℎ= 1 %, 𝑙𝑡ℎ/𝛿𝑡ℎ= 10, 𝑍𝑡ℎ= 0.3. 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.68 0.72 0.76 0.8 0.82 0.84 0.86 0.88 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.620.66 0.7 0.74 0.78 0.8 0.8 0.82 0.84 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.54 0.6 0.64 0.68 0.71 0.74 0.76 0.76 0.78 0.78 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.5 0.56 0.6 0.64 0.67 0.7 0.72 0.72 0.74 0.74 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 15% ; Mth = 0.2 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.440.48 0.52 0.56 0.6 0.64 0.68 0.7 0.7 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 56 Figure 54. Diffuser maps for S-CO2 at 𝑀𝑡ℎ= 0.2, 𝑅𝑒𝑡ℎ= 1.0·107, 𝑇𝑢𝑡ℎ= 1 %, 𝑙𝑡ℎ/𝛿𝑡ℎ= 10, 𝑍𝑡ℎ= 0.3. 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.2 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.68 0.72 0.76 0.8 0.82 0.84 0.86 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.2 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.620.66 0.7 0.74 0.780.8 0.82 0.84 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.2 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.54 0.6 0.64 0.68 0.71 0.74 0.760.78 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.2 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.5 0.56 0.6 0.64 0.67 0.7 0.72 0.74 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 15% ; Mth = 0.2 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.440.48 0.52 0.56 0.6 0.64 0.68 0.7 0.7 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 57 Figure 55. Diffuser maps for S-CO2 at 𝑀𝑡ℎ= 0.2, 𝑅𝑒𝑡ℎ= 1.6·107, 𝑇𝑢𝑡ℎ= 1 %, 𝑙𝑡ℎ/𝛿𝑡ℎ= 10, 𝑍𝑡ℎ= 0.3. 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.2 ; Reth = 1.6·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.64 0.68 0.72 0.76 0.80.82 0.85 0.88 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.2 ; Reth = 1.6·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.8 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.2 ; Reth = 1.6·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.2 ; Reth = 1.6·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 15% ; Mth = 0.2 ; Reth = 1.6·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 58 Figure 56. Diffuser maps for S-CO2 at 𝑀𝑡ℎ= 0.4, 𝑅𝑒𝑡ℎ= 7.4·106, 𝑇𝑢𝑡ℎ= 1 %, 𝑙𝑡ℎ/𝛿𝑡ℎ= 10, 𝑍𝑡ℎ= 0.3. 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.680.720.75 0.78 0.81 0.81 0.84 0.86 0.88 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.60.64 0.68 0.72 0.75 0.75 0.78 0.78 0.8 0.8 0.82 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.50.540.58 0.62 0.66 0.7 0.7 0.73 0.73 0.75 0.75 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.45 0.5 0.54 0.58 0.62 0.62 0.66 0.66 0.69 0.69 0.72 0.72 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 15% ; Mth = 0.4 ; Reth = 7.4·106 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.38 0.44 0.5 0.54 0.58 0.58 0.62 0.62 0.64 0.64 0.66 0.66 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 59 Figure 57. Diffuser maps for S-CO2 at 𝑀𝑡ℎ= 0.4, 𝑅𝑒𝑡ℎ= 1.0·107, 𝑇𝑢𝑡ℎ= 1 %, 𝑙𝑡ℎ/𝛿𝑡ℎ= 10, 𝑍𝑡ℎ= 0.3. 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 3% ; Mth = 0.4 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.68 0.72 0.75 0.78 0.81 0.81 0.84 0.84 0.86 0.88 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 6% ; Mth = 0.4 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.55 0.620.67 0.7 0.73 0.76 0.76 0.79 0.79 0.82 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 9% ; Mth = 0.4 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.50.55 0.6 0.64 0.68 0.72 0.72 0.75 0.75 0.78 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 12% ; Mth = 0.4 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.420.5 0.55 0.59 0.63 0.63 0.67 0.67 0.7 0.7 0.73 6 8 10 12 14 16 18 20 3 4 5 6 7 8 L/D AR Cp for Bth = 15% ; Mth = 0.4 ; Reth = 1.0·107 ; Tuth = 1% ; lth/th = 10 ; Zth = 0.3 2 = 10º 2 = 9º 2 = 8º 2 = 7º 2 = 6º 2 = 5º 2 = 4º 0.40.45 0.5 0.55 0.6 0.6 0.63 0.630.66 0.66 0.69 Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 66 𝜆=∫𝑣 𝑑𝐴 𝐴𝑢 ∫𝑣 𝑑𝐴 𝐴𝑏 [Equation 3.17] The next bullet points along with the results in Table 5 summarise the results with respect to the effect of distortion as obtained in the aforementioned thesis (López Florenciano, 2013):  Regardless of the operating fluid, the diffuser has an amplifying effect on distortion, evidenced by an increase in 𝜆 as the fluid flows through the diffuser.  This amplifying effect is stronger for air and it is found to be accentuated with the asymmetry of the velocity profile.  Supercritical carbon dioxide exhibits a slightly higher pressure rise coefficient, the differences between fluids being emphasised for those inlet profiles for which the radial velocity gradients are higher. Air Carbon dioxide 𝑫𝟏/𝑫𝟐 𝑪𝒑 𝝀𝒊𝒏 𝝀𝒐𝒖𝒕 𝑪𝒑 𝝀𝒊𝒏 𝝀𝒐𝒖𝒕 1.5/0.8 0.502 1.065 1.608 0.517 1.094 1.401 1.5/0.5 0.509 1.054 1.479 0.528 1.072 1.307 1.5/0.3 0.520 1.036 1.397 0.541 1.056 1.207 1.5/0.0. 0.540 1.006 1.021 0.569 1.042 1.044 Table 5. Summary of the effect of distortion on diffuser performance for 𝐴𝑅=4, 𝐿/𝐷𝑡ℎ=12 (López Florenciano, 2013). 3.4.2 Effect of swirl component Swirl modifies the flow structure inside the diffuser both in the core flow and in the boundary layer, bringing about important beneficial effects on pressure recovery, as well as in other engineering fields such as fuel-air mixing in reciprocating engines or homogenisation in coal combustion in boilers. These effects have also been evaluated numerically in (López Florenciano, 2013), making use of the flexibility of ANSYS Fluent for 2D axisymmetric swirl simulation. To this aim, the angular component of velocity is modelled to vary proportional to the radius by means of a udf and this is then added to the axial velocity distribution described in section 3.4.1 before. Although complex parameters have historically been defined to quantify the intensity of swirl [Equation 3.18] (Mc Donald, et al., 1971), the results are here presented as a function of the maximum swirl angle at the throat, Figure 66. 𝑚=𝐴𝑛𝑔𝑢𝑙𝑎𝑟 𝑚𝑜𝑚𝑒𝑛𝑡𝑢𝑚 𝐴𝑥𝑖𝑎𝑙 𝑚𝑜𝑚𝑒𝑛𝑡𝑢𝑚 =∫𝑣𝑚𝑣𝜃𝑟2 𝑑𝑟 𝑅 0𝑅∫𝑣𝑚 2 𝑟𝑑𝑟 𝑅 0 [Equation 3.18] 𝜃𝑚𝑎𝑥 =max[atan( 𝑣𝜃 𝑣𝑎𝑥𝑖𝑎𝑙)] [Equation 3.19] As opposed to the effect of inlet distortion, only soft differences were observed in the influence of swirl for each fluid. A very close look at figure 66 would only evidence that the swirl angle for peak Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 67 pressure rise is slightly higher for supercritical carbon dioxide than for air (24.5 vs. 23.5) as so is the improvement in 𝐶𝑝 (0.10 vs, 0.09), though the patterns are just alike, representing an improvement of 10 %. Figure 66. Summary of the effect of swirl on diffuser performance for the geometry 𝐴𝑅=5, 𝐿/𝐷𝑡ℎ=14 at 𝐵𝑡ℎ=6%, 𝑀𝑡ℎ=0.2, 𝑅𝑒𝑡ℎ=7.4·106 (López Florenciano, 2013). Additionally, the aforementioned effect of the swirl into the boundary layer flow structure causes in some cases even a null velocity in the axis at diffuser outlet. 3.4.3 Effect of turbulence There are various features of turbulence that have an effect on diffuser performance: intensity, length scale and structure (i.e. anisotropy), the first two of which were found to present critical values beyond which no further influence is to be expected in a channel diffuser, Figure 5 in ( Hoffman, 1981). Unfortunately, no systematic studies on turbulence have been done in the field of conical diffusers to date, even if Klein (Klein, 1981) carried out a complete compendium on the effect of turbulence amongst other parameters. In this work, the author concluded that “The effect of increasing inlet turbulence intensity is to increase pressure recovery”, which is true until a certain limit because, as reported in ( Hoffman, 1981), there seems to exists a value of intensity beyond which diffuser performance stabilises as shown in figure 67. Apart from this, this figure shows a slightly higher improvement of 𝐶𝑝 when 𝐿/𝐷𝑡ℎ is increased for S-CO2. This comes about because of the beneficial effect that higher lengths have on S- CO2 diffuser performance. Figure 67. Effect of turbulence intensity at the throat for S-CO2 (left) and air (right) for different geometries and Reynolds numbers (López Florenciano, 2013). 0.00 0.03 0.06 0.09 0.12 0 3 6 9 12 15 18 21 24 27 30 Cp-Cp swirl max = 0º Max Swirl Velocity Angle (º) Cp-Cpswirl=0 (Air) Cp-Cpswirl=0 (SCO2) Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 68 A critical value of non-dimensional turbulence length scale is easily identified in figure 68, which is in agreement with the results presented by Hoffman ( Hoffman, 1981) for channel diffusers. This critical value is around 100, meaning that the characteristic length scale of the turbulence phenomenon is in the order of magnitude of the throat diameter. Regarding the effect of fluid nature, no important differences are worth noting. Figure 68. Effect of non-dimensional turbulence length scale at throat for S-CO2 (left) and air (right) (López Florenciano, 2013). 3.5 Summary and conclusions Conical diffuser operation has been comparatively analysed for S-CO2 and atmospheric air by CFD simulation in ANSYS Fluent®. For that, an integral tool has been created by combining the CFD software with a mesh generator (ICEM CFD®) and also Matlab® and EES®. This tool permits to modify any parameter of influence in conical diffuser operation, not unlimited but till a certain extent. Nonetheless, for the range of analysis, the tool has been validated for air, since no experimental results are available for S-CO2 in this field. Therefore, an analysis of uncertainties has been performed for S-CO2 by considering different EOS and turbulence models in addition to apply the GCI analysis in order to determine numerical uncertainties due to the grid, the only source of errors. Although at the beginning it seems to exist huge differences between atmospheric air and S-CO2 in regard to conical diffuser operation, the development of this study has shown also the existence of important similitudes between both fluids. From a geometrical standpoint, the main difference to highlight is diffuser performance is enhanced only by lowering the divergence angle if S-CO2 is the working fluid, while for air it is possible to increase the pressure rise coefficient by increasing the area ratio. This is herein attributed to the much higher Reynolds number that characterizes S-CO2, as the higher inertia of this flow penalizes it adaption to divergent walls. However, this tendency is affected by the blockage factor, in such a way that this difference disappear at a 12 % blockage. It means that, at this level of analysis, the influence of fluid nature disappear in considerably distorted flows. It is another similitude the fact that blockage is the parameter with stronger influence. Moreover, as for air, Mach and Reynolds number have virtually no effect on the pressure rise coefficient. From the study on inlet velocity profiles, it was observed that: (i) air experiences a more pronounced distortion that the supercritical fluid, and by other side (ii) for those profiles with higher gradient of radial velocity, the operation with S-CO2 will be improved at higher extent in relation to air because of the lower weight of viscous forces. Finally, in regard to turbulence and swirl parameters, there weren’t observed important differences between both fluids. Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 69 3.6 Further calculations towards more reliable turbulence modelling Under the lack of experimental results in the field for the fluid of interest, the numerical approach employed in this project was compared against another one based on a more complex turbulence model, in contrast to the RANS turbulence models herein considered. DNS was discarded whereas focus was put on LES. However, since S-CO2 diffuser involves high Reynolds number and wall bounded flows, computational efforts are very strong in comparison with free flow simulations. In consequence, the Wall- Modeled Large Eddy Simulation (WMLES) (Shur, et al., 2008) approach was used for the comparison. WMLES belongs to the ELES (Embedded Large Eddy Simulation) approach, where two different zones are distinguished in the computational mesh: one where turbulent scales are modelled by using RANS models, and another one where turbulent scales are solved by Large Eddy Simulation. It is therefore a hybrid method that permits using 2-equation models in the near wall region, while the core flow can be directly solved by LES. Thus, three diffuser cases were run on WMLES computational models, requiring a total of around between 35,000 and 55,300 CPU hours per case, which would require 1 year and 7 months of just computational time if these calculations were performed with the GMTS computers and the academic licenses available. Therefore, this work was subcontracted to ANSYS, Inc, in parallel to the development of other tasks of this project. In general terms, the main conclusion was that results with RANS and WMLES were very similar, showing noticeable differences near the end of the cone, where stall flow becomes more likely. Next, details of ANSYS WMLES calculation are provided for the case 1 only, being the remaining cases very similar to the on presented here with results qualitatively identical. Case 𝑨𝑹 𝑳/𝑫𝒕𝒉 𝑴𝒕𝒉 𝑹𝒆𝒕𝒉 𝑩𝒕𝒉 Comments 1 3 8 0.4 7.4·106 0.03 Max. 𝐶𝑝 @ low 𝐵𝑡ℎ 2 5 8 0.4 7.4·106 0.12 Max. 𝐶𝑝 @ high 𝐵𝑡ℎ 3 3 10 0.2 7.4·106 0.03 High 2𝜃 (near stall) Table 6.Cases analysed in the ANSYS WMLES study. Six profiles were measured along each diffuser at 0, 20, 40, 60, 80 and 100 % of their lengths and velocity, pressure and temperature profiles were extracted and plotted in the table below, both for RANS (solid lines) and WMLES (dashed lines) calculations. As observed, results start differing apart from 80 % of the length, i.e. where stall flow may happen. However, mean values of pressures are very similar in both cases, offering very similar diffuser performance coefficients, 𝐶𝑝. Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 70 RANS WMLES Velocities Pressures Temperatures Table 7.RANS vs. WMLES for case 1. Figure 69. Velocity magnitude (in m/s) profiles comparison for case 1.RANS vs. WMLES comparison. Analysis of S-CO2 flows in elementary diffusion processes Design of supercritical carbon dioxide centrifugal compressors 71 Figure 70. Static pressure (in Pa) profiles comparison for case 1. RANS vs. WMLES comparison. NOTE: Operating pressure of reference: 7.5·106 Pa. Figure 71. Static temperature (in ºC) profiles comparison for case 1. RANS vs. WMLES comparison. Two main conclusions were extracted from this analysis:  Mean values of interest, i.e. static pressure to calculate afterwards performance coefficient; are very similar for both approaches. Therefore, diffuser maps obtained with the RANS-2D tool are considered numerically satisfactory, being still pending the corresponding empirical validation.  On the other hand, local values are noticeably different between both models, especially apart from the 60 % of the diffuser length. In consequence, measurements related to that such as peaks of velocity, wall shear stresses, etc. are not reliable in the 2D tool. From the light of these results, it could be deduced that RANS simulation provides satisfactory results for cases operating in the non-appreciable stall region (Figure 32), not being an accurate tool for predicting stall flow. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 72 CHAPTER 4. One-dimensional model of supercritical carbon dioxide centrifugal compressors In spite of the powerful simulation tools existing nowadays, one-dimensional techniques still constitute an essential part in the design process of a turbomachinery (Aungier, 1995). As an example, Dr. D. Japikse recommends the flow diagram of figure 72 for turbomachinery design (Japikse, 1996), where everything starts with one-dimensional (also called quasi-two-dimensional) evaluations. Moreover, because of their simplicity, these tools allow for quick evaluations of how a specific design would perform at given conditions, whether it be on or off-design operation. In addition to conservation laws, 1D models require experimental correlations for the evaluation of key performance parameters such as pressure loss coefficients and deviation angles amongst others. Nevertheless, it seems that, as far as supercritical carbon dioxide is concerned, there is not a body of knowledge which ensures the successful design of turbomachinery. This observation is based on the fact that the expected design efficiencies are not being reached by the very few existing machines at SANDIA National Laboratories, Bechtel Marine Propulsion Corporation and Tokyo Institute of Technology. In fact, a number of researchers have pointed out the need to develop a new design theory specific for supercritical carbon dioxide turbomachinery (Ulizar & Pilidis, 2000) (Wang, et al., 2004). In consequence, new proposals are included in this chapter aiming to contribute to fulfilling this gap. A literature review is first presented in section 4.1. Then the one-dimensional performance analysis is presented in detail (direct problem) and validated in the following section. Once the validity of the tool is confirmed to the extent possible (validation against experimental data), the equations of the performance model are re-organised to solve the indirect problem, which is presented in section 4.4. Finally, the main conclusions are exposed together with a summary of the chapter. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 73 Figure 72. Design strategy proposed by D. Japikse (Japikse, 1996). 4.1 Review of one-dimensional centrifugal compressor models A one-dimensional model is a numerical tool which solves for the thermo-fluid dynamic state of the gas/liquid at the interface between the different elements comprised by the turbomachinery under analysis. As such, the flow is characterised by the mean values of the relevant properties at each station of interest. Obviously, conservation laws of mass, momentum and energy, along with the equation of state must be taken into account by the model. Nevertheless, there exist several approaches to model the physics inside each compressor element, all of them potentially valid as long as its comparison against experiment shows good agreement. An author to highlight in the field of radial turbomachinery is Dr. David Japikse, whose work is a remarkable contribution to illustrate the physical phenomena that occur inside each turbomachinery One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 74 component. In addition, he has provided methods and tools both for solving the direct (analysis) and the indirect (design) problems. All this work has been managed through Concepts ETI Inc., created in 1980, which is nowadays one of the leading software developers for turbomachinery design worldwide. And not only this, the books by Dr. Japikse and co-workers, for instance (Japikse, 1996), are very illustrative and provide valuable graphical information about different aspect of turbomachinery. Focusing on one-dimensional models, Dr. Japikse pays particular attention to the two-zone model originally developed by (Dean & Senoo, 1960), which is based on splitting the flow at the impeller exit into an isentropic core flow (jet) and a non-isentropic distorted flow (wake) where all the impeller losses occur. There are two main parameters that characterise this model: the ratio of secondary to total mass flow rates 𝜒 and the secondary regime area fraction 𝜀, the dependency between which is correlated in figure 73. The model is closed by assuming a value of 𝜒 depending on the application, typically between 0.02 and 0.25 for efficient and poor design cases respectively. Implementing the two-zone model requires the assumption of 𝜒 which is well controlled for standard turbomachinery with conventional fluids but incurs large uncertainty when it comes to using supercritical carbon dioxide (since no experimental data exist). Concluding with this model, the flow is mixed after the impeller by applying mass, momentum and energy conservation. Therefore, as the author himself states, “the concept of using a blockage parameter at the rotor exit is fundamentally inconsistent”. Figure 73. Two-zone model results (Japikse, 1996). The Two-Elements-In-Series (TEIS) model is an alternative approach to compressor design (Japikse, 1996). This model assumes that the impeller channel can be modelled as two simple elements (see Figure 74) in series, the first one of which might be a nozzle or a diffuser, depending on the flow regime (high or low flow respectively), whilst the second one is always a diffuser. The performance of each separated device is affected by the parameters exposed in chapter 3, making it evident that the geometry and flow features at the inlet to both elements need to be identified in order to calculate their performances, whether by static pressure rise coefficients (𝐶𝑝,𝑎 and 𝐶𝑝,𝑏) or equivalent efficiencies (𝜂𝑎 and 𝜂𝑏). Again, experience plays a crucial role for closing this model as the elements’ performance parameters are fixed according to sizes and applications. From the point of view of supercritical carbon dioxide, serious doubts arise once more due to the limited database available. In this sense, chapter 3 tries to contribute to overcoming this lack of knowledge so that a complete set of information for the new fluid of interest is available in the future. In contrast to the two-zone model, the impeller outflow in the TEIS model is characterised by a blockage factor, as it is the case in diffuser theory. 1.00 0.90 0.80 0.70 0.60 0.50 0.40 0.30 0.20 0.10 0.00 0.00 0.20 0.40 0.60 0.80 1.00 NONROTATING, SINGLE PASSAGE DIFFUSERS ECKARDT DATA <3 ROTORS>. D2 = 0.4M CETI DATA <9 ROTORS>. D2 = 0.06M TO 0.13M χ ε χ = ε2 χ = 0.5ε2 - 0.05ε One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 75 Figure 74. Two-Elements-In-Series (TEIS) conceptual model. Slip is another fundamental feature of a radial impeller and it is conceptually similar to deviation in axial machinery. The formation of a bubble on the suction side at the outlet makes the streamlines deviate from the pressure to the suction side, bringing about a difference between blade and flow angles at the trailing edge of the blade. This fact has been investigated by several authors amongst which Dr. Japikse suggests using the one by Eckardt (Eckardt, 1980), expressed in [Equation 4.1]. 1−𝜎𝑠𝑙𝑖𝑝=𝑣2𝑢 𝑢2(1 𝜎′𝑠𝑙𝑖𝑝−1) [Equation 4.1] Similar to Japikse (Japikse, 1996), Cumpsty’s work (Cumpsty, 1989) makes a complete review of the models available for compressor design, where a compendium of the different approaches from onedimensional concepts such as the Jet and Wake and other loss-based models, to two-dimensional blade- to-blade and hub-to-shroud calculations is gathered. Additionally, this book contains design recommendations that applicable to conventional turbomachinery and also an important set of impeller flow measurements. Therefore, this reference emerges as an important source of information to find physical validation of the results of this project rather than as a specific guide for radial machinery analysis or design. Another work of relevance in the field of radial turbomachinery is the one developed by Ronald H. Aungier (Aungier, 2000). Throughout his carrier, R. H. Aungier has created methodologies both for turbomachinery analysis and design with methods based on loss models rather than two-zone models as it was the case for Japikse, approaches shared by several authors in the field (Yahya, 2005), (Balje, 1981), (Schobeiri, 2005). The main flaw of this approach is nevertheless that is virtually impossible to measure all loss terms in each turbomachinery element individually what, in turn, makes it impossible to validate them if not from a global standpoint. However, as it can be later seen, most loss models rely on the physics of the flow which makes them less dependent on fluid nature. This becomes a considerable advantage for the particular case of this project as the fluid is not conventional. Aside from conservation laws, Aungier’s models are based on the calculation of several loss contributions to finally calculate a total pressure drop coefficient for each element. These are all one-zone models, so no distinction between secondary and primary zones is made. Instead, aerodynamic blockage factor are considered at each station. As opposed to Japikse, who adopted Eckardt’s approach, Aungier developed a somewhat new slip factor expression by modifying the original work by Wiesner (Wiesner, 1967). His book (Aungier, 2000) is a compendium of all his work on radial compressors and constitutes a complete guide the aerodynamic design of centrifugal compressors, from the initial concept to the final three-dimensional geometry optimisation. Although more design-oriented, the work presented by A. Whitfield and N.C. Baines (Whitfield & Baines, 1990) also provides guidelines for the performance analysis of both for turbines and compressors. The approach is similar to Aungier’s but, in lieu of pressure loss coefficients, Whitfield and Baines’s proposal is based on entropy gain contributions, defined as 𝜎=ln(−∆𝑠/𝑅𝑔) and calculated by empirical Throat Element “a” at low flow regime Element “a” at high flow regime Element “b” Element “b”Element “a” Outlet OutletThroat Actual impeller channel Impeller in TEIS model One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 82 Figure 80. Flow diagram for the analysis of the inducer. 4.2.3 Impeller model Following the philosophy of the project, the fluid dynamic model of the impeller is based on pressure loss coefficients, as it is for the inducer. The flow is analysed in the relative frame of reference, for which loss coefficients apply to the relative dynamic pressure across the element. Impeller losses are classified as internal and parasitic, the latter not existing in stationary elements as diffusers and volute. Internal losses refer to the irreversibilities originated due to the non-ideal behaviour of the flow, which produces an increase in entropy or, equivalently, a total pressure losses. In contrast, parasitic losses have a mechanical nature, reducing the total enthalpy increase of the fluid with respect to the mechanical work produced by the shaft. There are three different sources of parasitic losses: disk friction (𝐼𝐷𝐹), recirculation flow (𝐼𝑅) and leakage flow (𝐼𝐿) (see figure 81). Inducer loss coefficients and the ones explained in sections 4.2.3.1 to 4.2.3.5 refer to internal losses, while parasitic losses are modelled in sections 4.2.3.6 to 4.2.3.8 by means of work input coefficients (non-dimensionalised with respect to 𝑢32). Total inlet, mass flow rate and Inducer geometry Section 1: [T1, P1, v1]=f(T01, P01, m, A1) Contraction losses, ϖcont Section 2: [T2, P2, v2]=f(P02, h02, m, A2) Energy conservation & Pressure losses Initialisation: No losses Loss calculation: ϖCH, ϖinc, ϖDIF Admissible εinc, εDIF ? Choking limit calculation NO Throat: [Tth, Pth, vth]=f(P0th, h0th, m, Ath) Energy conservation & Pressure losses END YES One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 83 Figure 81. Parasitic losses (Japikse & Baines, 1997) . 4.2.3.1 Wall friction losses Friction losses are evaluated with the model developed by (Schlichting, 1979), which provides a friction coefficient 𝑐𝑓 as a function of the characteristic Reynolds number of the duct (𝑅𝑒=𝜌𝑤𝑑𝐻/𝜇, averaged between inlet and outlet) and relative roughness (𝜀/𝑑𝐻) to be later applied in: 𝜛𝑓𝑟=4𝑐𝑓𝐿𝐵 𝑑𝐻(𝑤 𝑤2)2 [Equation 4.15] Where 𝐿𝐵 is the impeller channel length and 𝑑𝐻 the mean hydraulic diameter between impeller throat and exit. The mean relative velocity 𝑤 is computed as 𝑤=√(𝑤22+𝑤32)2 ⁄, being 𝑤≥ √(𝑤𝑡ℎ 2+𝑤32)2 ⁄ always. Two additional aspects must be noted with respect to the flow areas and the shape of the impeller channel (figure 82):  The inlet and outlet sections are assumed normal to the relative flow: 𝐴𝑖𝑛,𝑐ℎ𝑎𝑛𝑛𝑒𝑙= 𝐴2cos𝛽2 and 𝐴𝑜𝑢𝑡,𝑐ℎ𝑎𝑛𝑛𝑒𝑙=𝐴3cos𝛽3.  In order to account for the influence of the splitter blades on these losses, 𝐴3 is is calculated with respect to the total number of blades, 𝑧𝐹𝐵+𝑧𝑆𝐵, when calculating hydraulic diameter. Recirculation Disk friction Disk friction Leakage One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 84 Figure 82. Impeller channel geometry. This friction loss term considers the effect of the mean velocity of the fluid flow through the impeller channel. Therefore no local effects are considered, for which dedicated models are presented subsequently. 4.2.3.2 Aerodynamic loading losses Blade loading losses come about because of the deflexion of streamlines inside the impeller channels, which causes the diffusion of the relative flow. Some authors describe them as “the momentum loss due to boundary layer build-up” (Jansen, 1967). Given the much higher density of the fluid in comparison with air, and based on the fact that these losses stem from the change in angular momentum of the flow, a significant contribution of this loss term is expected in S-CO2 turbomachinery. Moreover, even if Its effect could be reduced by increasing the number of blades or reducing the blade curvature (i.e. flow deflexion), this would on the other hand penalise the performance of the impeller due to increased friction increase. These losses are evaluated with the model proposed in (Aungier, 2000), which requires a simple model for the meridional (stream-wise) velocity distribution along the pressure and suction sides of the blades. A simple (close to ideal) velocity profile can be obtained from a simple irrotational flow analysis and is depicted in figure 83. Figure 83. Blade loading model. Figure adapted from (Aungier, 2000). The velocity distribution shown in figure 83 yields the following maximum velocity difference: A2 θ r A3Aout,channel Ain,channel β 2 β 3 Impeller channel SUCTION SURFACE VELOCITY PRESSURE SURFACE VELOCITY Δw wmax w2 w3 ξ One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 85 ∆𝑤=2𝜋𝑑2 𝑢2 𝐼𝐵 𝑧𝑒𝑓𝑓 𝐿𝐵 [Equation 4.16] Which can be implemented into the equivalent diffusion ratio developed by Lieblein (Lieblein, 1965), which is a means to evaluate the aerodynamic loading of the blades: 𝐷𝑒𝑞=𝑤𝑚𝑎𝑥 𝑤3 ; 𝑤𝑚𝑎𝑥=(𝑤2+𝑤3+Δ𝑤)/2 [Equation 4.17] 𝐼𝐵 is the blade work input coefficient (described later in section 4.2.3.9) while the effective number of blades is calculated by having into account the splitter blades as well: 𝑧𝑒𝑓𝑓=𝑧𝐹𝐵+𝑧𝑆𝐵· (𝐿𝑆𝐵/𝐿𝐹𝐵). Therefore, the blade loading loss coefficient can be calculated as: 𝜛𝐵𝐿=1 24(∆𝑤 𝑤2)2 [Equation 4.18] This latter correlation accounts for the blade loading losses which are eventually caused by the pressure gradient in the blade-to-blade direction. Nevertheless, this is not the only pressure gradient in an impeller channel even though it is the most important. There exists another pressure gradient from hub to shroud whose effect must be taken into account by an additional loss term. This second contribution to the aerodynamic loading losses is taken into account by 𝜛𝐻𝑆 which depends on the mean streamline curvature, 𝜅𝑚=(𝛼𝐶3−𝛼𝐶2)/𝐿𝐵, mean channel width, 𝑙=(𝑙2+𝑙3)/2, and velocity profile. 𝜛𝐻𝑆=16(𝜅𝑚𝑙𝑤 𝑤2)2 [Equation 4.19] 4.2.3.3 Mixing losses Two separate contributions to mixing are consideredat the impeller outlet even though, in contrast to two-zone models, no secondary and primary zones are distinguished in the approach presented so far. Nevertheless, blockage is still modelled and hence the outflow is characterised by a tip blockage 𝐵3, directly related with the well-known distortion factor: 𝜆= 1 1−𝐵3 [Equation 4.20] At non-singular operating conditions of the compressor, the boundary layer is almost removed in the inducer thanks to the flow acceleration in the convergent geometry. Hence, an almost negligible boundary layer is expected at the inlet to the blading which then builds up in the impeller channel as the relative flow is diffused. Different phenomena are the cause of this build-up: (i) wall friction, (ii) diffusion of the relative flow, (iii) hub-to-shroud aspect ratio of the impeller channel and also (iv) relative importance of the clearance gap; all these effects contribute to increasing the outlet blockage. The following model proposed by Aungier (Aungier, 1995) summons all of them. 𝐵3=𝜛𝑓𝑟𝑃02−𝑃2 𝑃03−𝑃3√𝑤2 𝑑𝐻 𝑤3 𝑏3+[0.3+𝑏32 𝐿𝐵 2]𝐴𝑅 2 𝜌3 𝑏3 𝜌2 𝐿𝐵+𝛿𝐶𝐿 2 𝑏3 [Equation 4.21] Where the channel area ratio is calculated as 𝐴𝑅=𝐴3cos𝛽3/(𝐴2cos𝛽𝑡ℎ) and 𝛿𝐶𝐿 is the clearance gap size. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 86 The first mixing process is that of the distorted flow (included in the displacement boundary layer 12 ) and the core flow stream, causing pressure losses that are modelled with an abrupt expansion coefficient as proposed in (Benedict, et al., 1966). 𝜛𝜆=[(𝜆−1) 𝑣𝑚3 𝑤2]2 [Equation 4.22] The second mixing process takes place downstream of the trailing edge where the wake mixes with the channel flow. As in the previous case, this pressure loss coefficient is modelled as an abrupt expansion for which the size of the wake must be estimated. A diffusion limit of 𝐷𝑒𝑞=2 is assumed for stall-free operation in such a way that the separation velocity is calculated as follows: 𝑤𝑠𝑒𝑝={𝑤3 𝑖𝑓 𝐷𝑒𝑞≤2 𝑤3𝐷𝑒𝑞 2 𝐼𝑓𝐷𝑒𝑞>2 [Equation 4.23] Mixing involves meridional velocities only (Aungier, 2000), which are calculated before (𝑣𝑚3,𝑤𝑎𝑘𝑒) and after mixing (𝑣𝑚3,𝑚𝑖𝑥). 𝑣𝑚3,𝑤𝑎𝑘𝑒=√𝑤𝑠𝑒𝑝 2−𝑤𝑢3 2 𝑣𝑚3,𝑚𝑖𝑥=𝑣𝑚3 𝐴3 𝜋 𝑑3 𝑏3 [Equation 4.24] The destruction of absolute meridional velocity yields the wake mixing loss coefficient: 𝜛𝑤𝑎𝑘𝑒=[𝑣𝑚3,𝑤𝑎𝑘𝑒−𝑣𝑚3,𝑚𝑖𝑥 𝑤2]2 [Equation 4.25] 4.2.3.4 Clearance losses In open impellers, there exists a flow of gas through the gap between impeller tip and shroud from the pressure to the suction side of the blades. This portion of fluid is consequently driven by the impeller whereas it does not flow out from it, thus implying an additional pressure loss. Regarding S-CO2 turbomachinery, and given its reduced dimensions, the relative size of the gap can be noteworthy and this loos is likely to gain importance. The evaluation of clearance losses requires estimating both the pressure difference and the mass flow rate through the gap, to which aim angular momentum equilibrium is imposed: 12 Remember that the displacement boundary layer thickness (δ∗) is directly related with the aerodynamic blockage factor (See chapter 3). One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 87 Figure 84. Gap flow model (Aungier, 2000). 𝜏=𝑧𝑒𝑓𝑓 𝐹𝐵×𝑟 󰇍  [Equation 4.26] What, after mathematical rearrangement, yields the pressure difference across the gap: Δ𝑝𝐶𝐿=𝑚󰇗(𝑟3 𝑣𝑢3−𝑟2 𝑣𝑢2) 𝑧𝑒𝑓𝑓 𝑟 𝑏 𝐿𝐵 [Equation 4.27] On The other hand, the following flow velocity through the gap is proposed by Aungier (Aungier, 2000): 𝑣𝐶𝐿=0.816√2Δ𝑝𝐶𝐿/𝜌3 [Equation 4.28] Where the coefficient 0.816 comes from the simple hydraulic model shown in Figure 84 which comprises an abrupt contraction (from the pressure side to the gap) followed by an abrupt expansion (from the gap to the suction side). Let 𝐶𝑅𝐶𝐿 (defined <1) be the contraction ratio that the fluid experiences when flowing through the abrupt contraction and expansion. It is acknowledged, and it has already been said (see [Equation 4.7]), that the loss coefficient for both elements is (1−𝐶𝑅𝐶𝐿)2. Therefore, assuming that the density of the fluid remains constant, the application of Bernoulli’s equation yields 0.816=1[√2(1−𝐶𝑅𝐶𝐿)] ⁄, resulting in a contraction ratio equal to 𝐶𝑅𝐶𝐿=0.1334. In other words, this means that the flow area across the gap represents around 13% of the meridional section of the impeller channel. Nevertheless, in order to provide a more general approach, a variable coefficient depending on blade heights (𝑏2 y 𝑏3) and clearance gap (𝛿𝐶𝐿) is used in lieu of the constant 0.816 value: 𝑣𝐶𝐿=√Δ𝑝𝐶𝐿/𝜌3 1−𝐶𝑅𝐶𝐿 ; 𝐶𝑅𝐶𝐿=12(𝛿𝐶𝐿 𝑏2+𝛿𝐶𝐿+𝛿𝐶𝐿 𝑏3+𝛿𝐶𝐿) [Equation 4.29] The leakage flow through the gap is then: 𝑚󰇗𝐶𝐿=𝜌3 𝑧𝑒𝑓𝑓 𝑠 𝐿𝐵 𝑣𝐶𝐿 [Equation 4.30] And finally the clearance loss coefficient is: 𝜛𝐶𝐿=2 𝑚󰇗𝐶𝐿 Δ𝑝𝐶𝐿 𝑚󰇗 𝜌2 𝑤22 [Equation 4.31] τ F B= Δ pCL·AB r2 r3 rvCL Mechanical model Hydraulic model One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 88 4.2.3.5 Supersonic flow losses Although the inlet and outlet relative velocities are in the subsonic range, it is possible to find local relative velocities in the supersonic region, which form shock waves incurring additional losses. These losses are evaluated with the following loss coefficient (Aungier, 2000): 𝜛𝑐𝑟={0.4 [𝑤𝑚𝑎𝑥−𝑤∗ 𝑎2]2 𝑖𝑓 𝑤𝑚𝑎𝑥≥𝑤∗ 0 𝑖𝑓 𝑤𝑚𝑎𝑥<𝑤∗ [Equation 4.32] Where 𝑤∗ is the relative velocity at sonic conditions calculated in the inducer analysis (see page 81). 4.2.3.6 Disk friction losses The disk friction work input coefficient is calculated with the classical model of Daily and Nece (Daily & Nece, Mar, 1960) (Daily & Nece, Sept, 1960), based on the plot shown in figure 85 which gives the torque coefficient (𝐶𝑀0) as a function of the characteristic Reynolds number (𝑅𝑒𝑑𝑖𝑠𝑘=𝜌 𝜔 𝑟𝑑𝑖𝑠𝑘 2/𝜇) and the relative roughness (𝜀/𝑟𝑑𝑖𝑠𝑘) of the disk. 𝐶𝑀=2𝜏 𝜌𝜔2𝑟5 [Equation 4.33] Figure 85. Daily & Nece’s model to estimate the torque coefficient for disk friction calculations (Aungier, 2000). The final torque coefficient (defined in [Equation 4.33]) results from implementing some corrections given by Aungier into the values provided by Daily and Nece (Aungier, 2000). A first correction is done to account for the leakage flow according to the following expressions: 𝐶𝑀 ′=𝐶𝑀0(1−𝐾)2 (1−𝐾0)2 𝐾0=0.46 (1+𝛿𝑑𝑖𝑠𝑘 𝑟𝑑𝑖𝑠𝑘) ; 𝐾=𝐾0+𝐶𝑞(1.75 𝑣𝑢3 𝑢3−0.316)𝑟𝑑𝑖𝑠𝑘 𝛿𝑑𝑖𝑠𝑘 ; 𝐶𝑞=𝑚󰇗𝐿 𝑅𝑒𝑑𝑖𝑠𝑘 1/5 2𝜋 𝜌3 𝑟𝑑𝑖𝑠𝑘 2 𝑢3 [Equation 4.34] One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 89 This corrected coefficient is split into two due to the fact that Daily and Nece’s model yields the torque coefficient at both sides of the wheel (disk and cover). These partial coefficients are given by Aunger based on experimental data: {𝐶𝑀,𝑑𝑖𝑠𝑘=0.75 𝐶𝑀 ′ 𝐶𝑀,𝑐𝑜𝑣𝑒𝑟=0.75 𝐶𝑀 ′ 𝐿𝐵1−(𝑟2𝑠/𝑟3)5 𝑟3−𝑟2 𝑓𝑜𝑟 𝑐𝑙𝑜𝑠𝑒 𝑖𝑚𝑝𝑒𝑙𝑙𝑒𝑟𝑠 [Equation 4.35] Finally, disk friction losses are numerically characterised by the following work input coefficient: 𝐼𝐷𝐹=(𝐶𝑀,𝑑𝑖𝑠𝑘+𝐶𝑀,𝑐𝑜𝑣𝑒𝑟)𝜌3 𝑢3 𝑟32 2𝑚󰇗 [Equation 4.36] 4.2.3.7 Recirculation flow losses When the relative flow diffusion gets high, the risk of stall taking place must be taken into account as it would lead to partial flow recirculation back to the impeller channel. Therefore a consistent parameter for evaluating this situation is needed, which is usually taken as Lieblein’s equivalent diffusion ratio [Equation 4.17] being higher than 2 (Lieblein, 1959). This into account, the following work input coefficient due to recirculation losses is obtained: 𝐼𝑅=(𝐷𝑒𝑞 2−1)(𝑤𝑢3 𝑣𝑚3−2tan𝛽3)≥0 [Equation 4.37] The two conditions needed for recirculation to take place are (i) diffusion ratio higher than the threshold value, and (ii) deviation at impeller exit higher than: tan𝛼3′>2tan𝛽3 (as 𝑤𝑢3/𝑣𝑚3=tan𝛼3′) [Equation 4.38] 4.2.3.8 Leakage losses The design of efficient labyrinth seals for high pressure turbomachinery is anything but a simple task. A thorough analysis is usually needed to balance two opposing effects. Reducing the gaps in the labyrinth leads to lower leakage flows but, at the same time, increases the potential rub and tear of the fins, hence increasing mechanical losses and reducing the reliability and useful life of the component. On top of this, too narrow gaps make it more demanding to operate the unit during (cold) start-up and shutdown when thermal strains play an important role. This into account, any labyrinth seal will incur a certain leakage flow which reduces the performance of the turbomachine. And depends on the geometry of the seal and pressure ratio across it (in this particular case, the pressure ratio across the seal is the ratio from impeller exit to atmosphere). Based on this twofold dependence, there are two possible ways to calculate the leakage flow, depending on whether the impeller is open of covered. For open impellers, Aungier (Aungier, 1995) assumed that “half of the clearance gap leakage flow is re-entrained into the blade passage flow and re-energized by the impeller”. 𝐼𝐿,𝑜𝑝𝑒𝑛=𝑚󰇗𝐶𝐿 𝑣𝐶𝐿 2 𝑢3 𝑚󰇗 [Equation 4.39] However, for covered impeller, it is necessary to estimate the leakage mass flow rate, for which the classical paper by Egli (Egli, 1935) is used as proposed by Aungier (Aungier, 2000). One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 90 𝑚󰇗𝐿=𝜋 𝑑𝑠𝑒𝑎𝑙 𝛿𝑠𝑒𝑎𝑙 𝐶𝑟 𝐶𝑡 𝐶𝑐 𝜌3 √𝑅𝑇 [Equation 4.40] Where the seal contraction ratio (𝐶𝑟), seal throttling coefficient (𝐶𝑡) and seal carryover coefficient (𝐶𝑐) depend on the geometry of the seal and the number of fins. Figure 86.Labyrinth seals geometrical definition. Thus, for covered impeller, the leakage work input coefficient is calculated by: 𝐼𝐿,𝑐𝑜𝑣𝑒𝑟𝑒𝑑=𝐼𝐵𝑚󰇗𝐿 𝑚󰇗 [Equation 4.41] 4.2.3.9 The blade work input coefficient The blade work input coefficient is the non-dimensional enthalpy change of the fluid across the impeller. This is well established by Euler’s equation as the variation of 𝑣𝑢𝑢 from the inlet to the outlet. For ideally guided flow (𝛽3=𝛼3′), the tangential component of the absolute velocity at impeller exit is: 𝑣𝑢3 𝑖𝑑=𝑢3−𝑤𝑢3=𝑢3−𝑚󰇗 𝜌3𝐴3tan𝛽3 [Equation 4.42] As a consequence of the pressure difference between both sides of the blade, the flow deviates from the direction of the blade trailing edge towards the suction side. This deviation effect is accounted for by the slip factor (𝜎), which is defined as “the ratio between the actual blade work input and the value that would exist if the flow were perfectly guided by the blades” (Aungier, 2000). Another non-ideal effect is caused by the aerodynamic blockage at impeller exit, generally considered in radial machinery by means of the distortion factor (see section 4.2.3.3). Consequently, the blade work input coefficient is computed by the equation: 𝐼𝐵=𝜎(1−𝑚󰇗 𝜆 tan𝛽3 𝜌3 𝐴3 𝑢3)−𝑢2 𝑣𝑢2 𝑢32 [Equation 4.43] The slip factor is calculated based on Wiesner’s equation (Wiesner, 1967) with the correction implemented by Aungier (Aungier, 1995). As a matter of fact, Aungier’s correction makes adopts Wiesnser’s factor if the meanline radius ratio 𝜖 is below its limiting value and a modified expression when this limiting value is exceeded: 𝜎𝑊𝐼𝐸𝑆𝑁𝐸𝑅=1−√cos𝛽3sin𝛼𝐶2/𝑧𝑒𝑓𝑓 0.7 [Equation 4.44] 𝜖𝑙𝑖𝑚=𝜎𝑊𝐼𝐸𝑆𝑁𝐸𝑅−𝜎∗ 1−𝜎∗ ; 𝜎∗=sin(37𝑜−0.2𝛽3) [Equation 4.45] rseal δseal tseal Pseal FIN AXIS One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 91 𝜎={𝜎𝑊𝐼𝐸𝑆𝑁𝐸𝑅 𝐼𝑓 𝜖≤𝜖𝑙𝑖𝑚 𝜎𝑊𝐼𝐸𝑆𝑁𝐸𝑅[1−(𝜖−𝜖𝑙𝑖𝑚 1−𝜖𝑙𝑖𝑚)√9𝑜−𝛽3/10] 𝐼𝑓 𝜖>𝜖𝑙𝑖𝑚 [Equation 4.46] 4.2.3.10 Algorithm for the calculation of impeller performance Taking the solution of the inducer flow as an input, it is possible to calculate the ideal discharge of the impeller if null deviation and isentropic flow are considered. Then, a first estimate of the flow field is obtained, which is later used to calculate the various losses and other coefficients that are necessary to calculate the real impeller flow. Moreover, as opposed to the inducer and other stationary elements, the impeller analysis is done in the relative reference frame which gives place to the following equation for the conservation of rothalpy and the calculation of total pressure losees: ℎ03 ′=ℎ02 ′+(𝑢32−𝑢22)/2 𝑃03 ′=𝑃03𝑖𝑑 ′−𝑓𝑐 (𝑃02 ′−𝑃2)∑𝜛𝑖 2→3 [Equation 4.47] The correction factor 𝑓𝑐 is employed to account for the fact that the losses are applied at the impeller outlet (i.e., they are applied to the total pressure at the outlet from this component) but they are defined as a fraction of the inlet dynamic head. In accordance, the factor mitigates errors caused by high pressure differences between impeller inlet and outlet and, even though the one employed by Aungier is 𝑓𝑐=(𝜌03 ′𝑇03 ′)(𝜌02 ′𝑇02 ′)⁄ (Aungier, 1995), it is herein redefined as: 𝑓𝑐=𝑃03 ′/𝑃02 ′ [Equation 4.48] The factor translates the relative dynamic pressure at the inlet (𝑃02 ′−𝑃2) to comparable terms at the outlet, a fraction of which is later subtracted from the ideal relative total pressure at this section. It is noted here that for adiabatic and isentropic processes involving ideal gases, Aungier’s original factor and the factor proposed herein are ust the same since 𝑃03 ′/𝑃02 ′=(𝜌03 ′𝑇03 ′)(𝜌02 ′𝑇02 ′)⁄ . Nevertheless, since supercritical carbon dioxide is not an ideal gas, the ratio of relative total pressures is preferred. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 98 4.2.5.1 Incidence losses Following the philosophy applied for the impeller, incidence losses are the consequence of two effects: the difference between the actual flow angle and the optimum one (in terms of minimum loss), and the abrupt contraction caused by the blade leading edges (i.e, the effect of blade thickness). In addition, if the diffuser is stalled (𝑣4>𝑣4𝑠), a supplementary term must be introduced in this loss contribution. 𝜛𝑖𝑛𝑐,𝐷𝐼𝐹= { 𝜛𝑖𝑛𝑐0,𝐷𝐼𝐹+0.8(𝑣4−𝑣4∗ 𝑣4)2 𝐼𝑓 𝑣4≤𝑣4𝑠 𝜛𝑖𝑛𝑐0,𝐷𝐼𝐹+0.8{((𝑣4 𝑣4𝑠)2−1)(𝑣𝑡ℎ,𝐷𝐼𝐹 𝑣4)2+(𝑣4𝑠−𝑣4∗ 𝑣4𝑠 )2} 𝐼𝑓 𝑣4>𝑣4𝑠 [Equation 4.68] Where the minimum incidence losses in the diffuser is given by: 𝜛𝑖𝑛𝑐0,𝐷𝐼𝐹=0.8(𝑣4∗−𝑣𝑡ℎ 𝑣4)2+(𝑧𝐷𝐼𝐹𝑡4 2𝜋𝑟4)2 [Equation 4.69] 4.2.5.2 Choking Choking losses in the diffuser are based on the same model employed for the impeller (see section 4.2.2, [Equation 4.12]). The calculation of the flow conditions at the throat is thus necessary, for which incidence losses are taken into account. 4.2.5.3 Wall friction losses 𝜛𝑓𝑟,𝐷𝐼𝐹=4𝑐𝑓𝐿𝐵,𝐷𝐼𝐹 𝑑𝐻,𝐷𝐼𝐹(𝑣 𝑣4)2(𝑑𝐻,𝐷𝐼𝐹 2𝛿 )0.25 [Equation 4.70] As observed in the equation above, friction losses are calculated as in the impeller [Equation 4.15] though affected by the factor (𝑑𝐻,𝐷𝐼𝐹 2𝛿 ⁄)0.25, which includes a correction due to the fact that boundary layers in diffusers are not fully-developed usually. 4.2.5.4 Mixing losses As a consequence of the mixing process downstream of the diffuser trailing edges, a pressure loss arises. As for the impeller, this is a function of the meridional wake and mixed-flow velocities, 𝑣𝑚5,𝑚𝑖𝑥 and 𝑣𝑚5,𝑚𝑖𝑥 respectively. 𝜛𝑚𝑖𝑥,𝐷𝐼𝐹=(𝑣𝑚5,𝑤𝑎𝑘𝑒−𝑣𝑚5,𝑚𝑖𝑥 𝑣4)2 [Equation 4.71] where the wake velocity is estimated from a separation velocity (𝑣𝑆𝐸𝑃) and the mixed-flow velocity is given by continuity: 𝑣𝑚5,𝑤𝑎𝑘𝑒=√𝑣𝑆𝐸𝑃 2−𝑣𝑢5 2 [Equation 4.72] 𝑣𝑚5,𝑚𝑖𝑥=𝑣𝑚5 𝐴5 𝜋 𝑑5 𝑏5 [Equation 4.73] One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 99 The separation velocity depends on the divergence angle through the parameter 𝐶𝜃, defined in section 4.2.5.5. 𝑣𝑆𝐸𝑃=𝑣4 1+2𝐶𝜃 [Equation 4.74] 4.2.5.5 Aerodynamic blockage at vaned diffuser outlet An empirical correlation presented in (Aungier, 1988) is used for the calculation of blockage [Equation 4.75]. It was seen in chapter 3 that flow at the outlet from the diffuser is strongly dependent upon throat conditions, especially the aerodynamic blockage. However, according to the model employed, the blockage at diffuser outlet does not depend upon the blockage at diffuser inlet, the rationale being that the acceleration of the flow from the inlet to the throat brings about a negligible boundary layer thickness at the latter section. 𝐵5=[𝐾1+𝐾2(𝐶𝑅2−1)]𝐿𝐵,𝐷𝐼𝐹 𝑙5,𝐷𝐼𝐹 [Equation 4.75] Where 𝐶𝑅=12(𝑣𝑚4cos𝛽5 𝑣𝑚5cos𝛽4+1) 𝐾1=0.2(1− 1 𝐶𝜃 𝐶𝐿 ) 𝐾2=2𝜃𝐶 125𝐶𝜃 (1−2𝜃𝐶 22𝐶𝜃) [Equation 4.76] Where the width of the flow passage is calculated as 𝑊=(2𝜋𝑟cos𝛽)/𝑍𝐷𝐼𝐹−𝑡 and 𝐶𝜃 and 𝐶𝐿are correction coefficients associated to the equivalent diffuser divergence angle (2𝜃𝐶) and blade loading parameter (𝐿), which are calculated as follows: 2𝜃𝐶=2atan[𝑙5𝑏5𝑏4 ⁄−𝑙4 2𝐿𝐵,𝐷𝐼𝐹 ] 𝐿= Δ𝑣 𝑣4−𝑣5 =2π(𝑟4𝑣𝑢4−𝑟5𝑣𝑢5)(𝑧𝐷𝐼𝐹𝐿𝐵,𝐷𝐼𝐹) ⁄ 𝑣4−𝑣5 [Equation 4.77] It is worth noting that a pronounced deterioration of diffuser performance has been observed experimentally when the loading parameter exceeds 1/3 or the divergence angle is higher than 11º. Therefore, coming back to correction coefficients, they are calculated as: 𝐶𝜃=max(1,2𝜃𝐶/11) 𝐶𝐿=max(1,3𝐿) [Equation 4.78] 4.2.5.6 Deviation at the vaned diffuser outlet The calculation of the outflow angle is based on Howell’s work on axial blade cascades (Howell, 1947) and thus requires a transformation of the coordinate system. At the minimum-loss incidence, the deviation is: One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 100 𝛿∗=𝜃[0.92 (𝑎/𝑐)2+0.02 𝛽5] √𝜎−0.02𝜃 [Equation 4.79] Where the cascade parameters (relative location of maximum camber, 𝑎/𝑐, solidity, 𝜎, and camber angle, 𝜃) are: 𝑎𝑐 ⁄=2−(𝛽−𝛽4)(𝛽5−𝛽4)⁄ 3 𝜎=𝑧𝐷𝐼𝐹(𝑟5−𝑟4)(2𝜋𝑟4cos𝛽) ⁄ 3 𝜃=𝛽4−𝛽5 [Equation 4.80] But diffusers do not usually operate at the minimum-losses incidence, and therefore the effect of incidence on deviation must be computed by means of the following derivative: 𝜕𝛿 𝜕𝑖=exp{[(1.5−𝛽4/60)2−3.3]𝜎} [Equation 4.81] Finally, the glow angle at the outlet from the diffuser is calculated, at any incidence angle, as: 𝛼5=𝛽5−𝛿=𝛽5−𝛿∗−𝜕𝛿 𝜕𝑖(𝛽4−𝛼4) [Equation 4.82] 4.2.5.7 Algorithm for the analysis of vaned diffusers As it was the case for the impeller, the vaned diffuser is initialised with ideal discharge conditions, i.e. no deviation and no-losses. Once a flow field is obtained, it is possible to compute losses and hence update the flow field at the end of each iteration. Convergence criteria are set for static pressure and velocity at diffuser outlet. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 101 Figure 91. Flow diagram for the analysis of vaned diffusers. 4.2.6 Volute model The high pressure fluid delivered at the diffuser outlet is collected by the volute, which still may produce a certain pressure rise thanks to flow diffusion. This process brings about three different contributions to the total loss of the compressor: Figure 92. Volute geometry definition  The “destruction” of the meridional component of velocity at diffuser exit, since the geometry of the volute imposes a tangential velocity in the full-collection plane (Figure 92). Diffuser inlet from vaneless space resolution No deviation, isentropic & enthalpy conservation Ideal discharge: [T5id, P5id, v5id]=f(h05, s05, m, A5) Initialization Real discharge: [T5, P5, v5]=f(h05, s05, m, A5) Loss calculation: ϖfr,DIF, ϖinc,DIF, ϖCH,DIF, ϖmix,DIF Other calculations: B5, checking stall Admissible εP5, εv5 ? END YES NO Energy conservation & Pressure losses r5 r6Afcp One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 102 𝜛𝑚=(𝑣𝑚5/𝑣5)2 [Equation 4.83]  With regard to the tangential component, additional losses are generated when the angular momentum is not conserved, which is evaluated by the so-called sizing parameter: 𝑆𝑃=𝑟5 𝑣𝑢5 𝑟6 𝑣𝑢6 [Equation 4.84] 𝜛𝑈= { 12𝑟5 𝑣𝑢5 2 𝑟6 𝑣52(1− 1 𝑆𝑃2) 𝐼𝑓 𝑆𝑃≥1 𝑟5 𝑣𝑢5 2 𝑟6 𝑣52(1− 1 𝑆𝑃 )2 𝐼𝑓 𝑆𝑃<1 [Equation 4.85]  And, finally, friction losses exist and they are evaluated with: 𝜛𝑓𝑟,𝑉=4 𝑐𝑓(𝑣6 𝑣5)2𝐿𝑉 𝑑𝐻,𝑉 𝐿𝑉=𝜋(𝑟5+𝑟6)/2 𝑑𝐻,𝑉=√4 𝐴6/𝜋 [Equation 4.86] 4.2.6.1 Algorithm for the analysis of the volute The velocity at the outlet from the volute is initialised with a default value, for instance 𝑣5/2, which allows to calculate the three loss terms that participate in the volute. Then, the volute outlet flow is updated and the convergence (velocity at the volute exit in this case) criterion checked. Figure 93. Flow diagram for volute analysis. Volute inlet from vaned diffuser resolution Initialization of v6 Real discharge: [T6, P6, v6]=f(h06, s06, m, A6) Loss calculation: ϖfr,VOL, ϖm,VOL, ϖU,VOL Admisible εv6 ? END YES NO Energy conservation & Pressure losses One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 103 4.3 Validation of the model The model presented in earlier sections is now validated both for supercritical carbon dioxide and air. The radial compressor existing at SANDIA National Laboratories (Wright, et al., September 2010) is taken for validation with carbon dioxide, while the NASA CC3 compressor (Tan, 2003) is used for air. 4.3.1 Validation for supercritical carbon dioxide Table 9 presents the geometrical parameters that define the supercritical CO2 compressor at SANDIA NL, which has been used for validation. Impeller data are explicitly indicated in the document (Wright, et al., September 2010) while for the diffuser some parameter need to be calculated from graphical information available in the same document. Additional experimental data is found in (Vilim, 2010) as well. The geometry of the diffuser geometry is estimated by measuring the dimensions in figure 5-35 in (Wright, et al., September 2010), where the component can be seen completely (see figure 94). It is noted though that the diffuser modelled is of the vaned type (i.e. aerodynamic aerofoils) whereas the one used in the validation compressor is of the channel type; nevertheless, in order to minimise the effect of this difference, the diffuser geometry was set up with the objective of keeping the divergence angle and the area ratio identical: ~8.5 º and ~2 respectively. Figure 94. Full S-CO2 turbo-compressor at SANDIA NL (Wright, et al., September 2010). One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 104 IMPELLER GEOMETRY Number of full blades 6 Number of splitter blades 6 Impeller inlet radius at hub 2.537585 mm Impeller inlet radius at shroud 9.372047 mm Impeller exit radius 18.68170 mm Blade angle of the impeller leading edge at hub 17.88 º Blade angle of the impeller leading edge at mean radius 37.13 º Blade angle of the impeller leading edge at shroud 50.00 º Blade angle of the impeller trailing edge -50.00 º (backward) Angle between streamlines and shaft at impeller inlet 0 º Angle between streamlines and shaft at impeller exit 90 º Full blade length 25.0 mm Splitter blade length 12.5 mm Axial length of the impeller 15.9 mm Blade thickness at impeller leading edge 0.762 mm Blade thickness at impeller trailing edge 0.762 mm Blade height at impeller leading edge 1.7 mm Clearance gap at impeller tip 0.254 mm VANED DIFFUSER GEOMETRY Number of vanes 17 Diffuser inlet radius 18.5 mm Diffuser exit radius 26.0 mm Blade angle at diffuser inlet 71.50 º Blade angle at diffuser exit 42.44 º Blade height at diffuser inlet 1.8 mm Blade height at diffuser exit 1.8 mm Diffuser channel length 10.6 mm Blade thickness at diffuser inlet 0.0 mm Blade thickness at diffuser exit 3.35 mm Table 9. Geometry employed in the simulation of the SANDIA NL compressor used for validation with S-CO2. Once the geometry of the component is defined, the total conditions at the inlet, mass flow rate and shaft speed are set for each validation case, whereas static pressure at impeller exit and total pressure and temperature at diffuser outlet are taken as the reference parameters to be compared against experiments. The results of the validation are presented in table 10 below, where last three points are taken from (Vilim, 2010) and the remaining ones from (Wright, et al., September 2010). One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 105 N (rpm) 𝑻𝟎𝟏 (K) 𝑷𝟎𝟏 (bar) 𝒎󰇗 (kg/s) 𝑷𝟑 (bar) 𝑷𝟎𝟓 (bar) 𝑻𝟎𝟓 (K) Exp. Mod. Dev. (%) Exp. Mod. Dev. (%) Exp. Mod. Dev. (%) 10000 305.5 76.76 0.454 76.76 77.40 0.83 79.79 77.74 1.28 - - - 20000 305.5 76.76 0.771 78.54 79.49 1.21 80.69 81.11 0.53 - - - 28000 305.5 76.76 1.134 82.11 82.06 -0.06 85.33 85.19 -0.16 - - - 39000 305.6 77.11 1.451 85.68 87.88 2.54 92.82 94.55 1.86 - - - 49000 306.3 78.54 1.816 94.25 95.46 1.28 106.39 106.18 -0.20 - - - 55000 306.4 78.90 2.043 99.96 100.53 0.57 113.53 114.41 0.78 - - - 56000 306.6 78.26 2.088 101.04 101.54 0.49 114.96 115.85 0.77 - - - 60000 306.9 79.97 2.225 102.11 105.69 3.51 121.39 122.37 0.81 - - - 64900 307.9 82.11 2.406 108.53 111.98 3.18 129.24 131.59 1.82 - - - 64384 308.71 82.86 2.860 106.7 108.94 2.10 119.4 108.94 4.30 323.82 324.17 0.11 29888 306.78 79.20 1.315 82.64 84.86 2.69 85.68 88.01 2.72 310.09 4 310.05 6 -0.01 59584 308.33 82.24 2.609 102.60 104.89 2.23 112.28 118.5 5.54 321.64 321.88 4 0.08 Table 10. Comparison of the experimental and analytical performances for supercritical carbon dioxide operation. Given the complexity of the model, a mean deviation of 1.67 % with a peak at 5.54 % for the relevant parameters being screened, along with a deviation lower than 1 % in delivery temperature, are deemed satisfactory. Therefore, the model is considered to be suitable for its implementation into a design and optimisation tool to be used in supercritical carbon dioxide operation. 4.3.2 Validation for air The suitability of the model to simulate the performance of radial compressors operating on air is also checked. To this aim, the NASA CC3 compressor is considered, for which data is available in (Tan, 2003) and (Landon Tarr, 2008). In this case, the comparison between the experimental and simulated performance is done for the pressure ratio and efficiency vs mass flow rate curves at constant-speed (60, 70, 80, 90 and 100 % the design value, 21798 rpm). The results of this comparison are shown graphically in figure 95, and the geometrical data is available in Table 11. As it was the case for S-CO2, the diffuser is of the channel type and is modelled by modifying the trailing edge angle in order to keep the divergence angle in the actual case, while the leading edge angles of impeller and diffuser are set to have null incidence at the design point (21789 rpm and 4.54 kg/s). One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 106 IMPELLER GEOMETRY Number of full blades 15 Number of splitter blades 15 Impeller inlet radius at hub 41 mm Impeller inlet radius at shroud 105 mm Impeller exit radius 73 mm Blade angle of the impeller leading edge at hub 24.8695 º Blade angle of the impeller leading edge at mean radius 39.5336 º Blade angle of the impeller leading edge at shroud 49.8897 º Blade angle of the impeller trailing edge -50.00 º (backward) Angle between streamlines and shaft at impeller inlet 0 º Angle between streamlines and shaft at impeller exit 90 º Full blade length 197 mm Splitter blade length 137.9 mm Axial length of the impeller 133.89 mm Blade thickness at impeller leading edge 7.58 mm Blade thickness at impeller trailing edge 7.58 mm Blade height at impeller leading edge 17 mm Clearance gap at impeller tip 0.81 mm VANED DIFFUSER GEOMETRY Number of vanes 24 Diffuser inlet radius 232.76 mm Diffuser exit radius 450.813 mm Blade angle at diffuser inlet 72.1222 º Blade angle at diffuser exit 30 º Blade height at diffuser inlet 17 mm Blade height at diffuser exit 17 mm Diffuser channel length 288.794 mm Blade thickness at diffuser inlet 0.0 mm Blade thickness at diffuser exit 49.99 mm Table 11. . Geometry employed in the simulation of the CC3 compressor at NASA Glenn Research Centre with air. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 107 Figure 95. Comparison between numerical and experimental performance data of the CC3 compressor at NASA Glenn Research Centre (Tan, 2003) . Having into account the difficulties founds when defining the geometry and accordingly the inherent uncertainty of this process, in the author considers that the results presented in figure 95 are acceptable, in particular for operating conditions close to the design point, i.e. far from choking (performance drop at high mass flow) and stall. The fact that this model has been developed for design purposes makes these results even more satisfactory and thus it is concluded that this research can move on to the next stage. 4.4 Design strategy for supercritical CO2 centrifugal compressors: Inverse problem The inverse or design problem consists of calculating the geometry that yields the desired performance a compressor given a set of operating conditions, i.e. inlet stagnation point, mass flow rate and pressure ratio for a specific gas, whether it be a pure substance or a mixture. Although the equations involved are the same as in the direct problem (see section 4.2), they have to be solved through a different algorithm, which gives place to the design strategy. This section presents a design strategy which accounts for general rules in radial turbomachinery design but also includes specific considerations for supercritical fluids. From the concept of the power cycle to the final design of the compressor, the process can be divided into three steps:  Firstly, the operating cycle is design for given net power, which means determining the total pressure and temperature at the inlet and outlet sections of each component and also the mass flow rate. To this aim, component efficiency and other performance parameters need to be guessed.  The compression system can be designed once the inlet conditions, pressure ratio and mass flow rate are known. A target efficiency is also fixed by selecting appropriate values of specific speed (𝑁𝑠) and diameter (𝐷𝑠), which eventually yields shaft speed and impeller diameter for each compression stage.  Finally, the geometry for each stage is obtained by solving the inverse problem. 4.4.1 Cycle definition According to the solar-oriented application of the project, the design concept is based on an S- CO2 power cycle for a 10 MW central receiver solar plant. The turbine inlet temperature is set to 1100 K, which is an affordable value for this type of technology (central receiver) even though the existing plants still work at a lower temperature (in the order of 850 K for a multi-megawatt application). The maximum pressure of the cycle is on the other hand limited to 250 bar which brings about a much more demanding 1 2 3 4 5 1 1.5 2 2.5 3 3.5 4 4.5 5 Total-to-toal pressure ratio Mass flow rate (kg/s) PR vs PR_t_t 100 % 90 % 80 % 70 % 60 % 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1 1.5 2 2.5 3 3.5 4 4.5 5 Total-to-total efficiency Mass flow rate (kg/s) eta vs. eta_t_t 100 % 90 % 80 % 70 % 60 % One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 114 Figure 100. Absolute and relative Mach numbers at the inducer throat of the first stage of a compression system with 2 (left) and 3 (right) stages. The two-stage system does not allow for designing the stage with a Mach number lower than 0.62 in the throat, which is risky taking into account that AMC = 0.689 (red line). Additionally, the relative Mach number for this alternative is always over 0.95. In contrast, when a three-stages system is adopted, it is possible to employ relative Mach numbers at the throat lower than the AMC. Therefore the threestages option is finally selected. 4.4.3.2 Impeller design Once the inlet section of the impeller has been sized in the previous section, the design of the entire impeller implies calculating the blade angle at impeller exit (𝛽3), for which the desired stage pressure ratio is achieved. it must be noted that even if the impeller exit radius (𝑟3) plays an important role in this issue as well, it is not considered as a degree of freedom for it is determined by the desired specific diameter (section 4.4.2). It must be also noted that pressure losses occur downstream the impeller which have an impact on stage pressure ratio. Therefore, coefficients of performance of both diffusers and volute must be initially guessed prior to closing the design of the impeller. In few words, given the mutual influence between the different elements of the compressor, this draft design will later have to be refined to account for the performance of the remaining of the machine. Summing up, impeller design provides the impeller exit angle and requires the following inputs:  The impeller exit radius (𝑟3), imposed by the specific diameter. The total-to-total efficiency (𝜂𝑡𝑡) of the stage and the total pressure loss coefficient of the system formed by diffuser and volute, [Equation 4.93], must be guessed in order to get the desired conditions at impeller exit. The development of the design process requires to update this parameter for convergence. 𝐾𝑑𝑖𝑓𝑓−𝑣𝑜𝑙=𝑃03−𝑃06 𝑃03 [Equation 4.93]  The ratio of meridional velocities between impeller inlet and outlet is a design choice commonly set to 1. The design loop for the impeller is as follows: 1. A value of 𝐾𝑑𝑖𝑓𝑓−𝑣𝑜𝑙 is assumed and kept constant throughout the entire impeller design process. 2. A value of 𝜂𝑡𝑡 is initially guessed. 0.045 0.05 0.055 0.06 0.065 0.07 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 radius (m) Mach ( ) Absolute Mach number at throat Relative Mach number at throat Maximum Mach established 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.4 0.5 0.6 0.7 0.8 0.9 radius (m) Mach ( ) Absolute Mach number at throat Relative Mach number at throat Maximum Mach established One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 115 3. The absolute total conditions at impeller exit are defined by pressure and enthalpy, calculated by means of the following system of equations: { ℎ03=ℎ02+Δℎ0𝑠 𝜂𝑡𝑡 𝑃03=𝑃01 𝑃𝑅𝑠𝑡𝑎𝑔𝑒 1−𝐾𝑑𝑖𝑓𝑓−𝑣𝑜𝑙 [Equation 4.94] Where Δℎ0𝑠 is the isentropic total enthalpy change of the stage, calculated from the desired total-to-total pressure ratio. The equation for 𝑃03 is derived from the proper definition of 𝐾𝑑𝑖𝑓𝑓−𝑣𝑜𝑙 taking into account that 𝑃𝑅𝑠𝑡𝑎𝑔𝑒=𝑃06/𝑃01. 4. 𝑣𝑢3 is calculated with Euler’s equation while the meridional component (𝑣𝑚3) is a design choice, usually 𝑣𝑚3=𝑣𝑚2. 5. The performance parameters of the impeller parameter are assumed ideal, i.e. 𝜎=𝜆= 1 and ∑𝜛𝑖2→3 =0. 6. The impeller blade exit angle is calculated by 𝛽3=atan[(𝑣𝑢3/𝜎−𝑢3)(𝜆𝑣𝑚3)⁄ ] 7. The velocity diagram at impeller exit is solved and, consequently, the static conditions. 8. By continuity, the blade height at impeller exit (𝑏3) is calculated. 9. The impeller is thus completely defined 16 . Therefore, it is simulated with the model presented in section 4.2.3 in order to update slip and distortion factors as well as losses in it. In consequence, it is necessary to iterate from points 6 to 10 to achieve convergence. 10. Once the impeller is consistently solved, the new outlet conditions can be calculated since 𝐾𝑑𝑖𝑓𝑓−𝑣𝑜𝑙 is still the same. The total-to-total efficiency is then updated, giving place to an external iteration loop from points 3 to 10. 4.4.3.3 Vaneless space design There usually exists a vaneless space between impeller and diffuser with the objective of homogenising the flow and produce a certain flow diffusion at the same time. The design of this part comprises the definition of the blade height (𝑏4), blade angle (𝛽4) and radius (𝑟4) at the leading edge of the diffuser blades. [Equation 4.95] is an empirical correlation proposed in (Aungier, 1988) to size the radius: 𝑟4 𝑟3=1+90−𝛼4 360 +𝑀32 15 [Equation 4.95] This space provides the necessary path downstream of the impeller to reduce the distortion of the flow coming out from the latter element, thus reducing the locally high Mach number to the extent possible. An appropriate design of the angle 𝛽4 should guarantee a wide enough stable operating range, i. e. stall-free operation. Figure 101 shows how incidence and inlet flow angle in stall conditions vary with diffuser leading edge angle for different Mach numbers and considering 20 vanes in all cases. These curves are result from solving [Equation 4.64] to [Equation 4.66] for the various cases. As a general rule of thumb, it is concluded that the lower the Mach number, the higher the stall margin, which justifies the interest of diffusing the flow before going into the diffuser. However, this parameter is not a design choice as it comes imposed by impeller design. In contrast, the diffuser leading edge angle is here a design result that is defined by incidence (typically around -1º (Aungier, 2000)). The flow angle at this point is controlled by blade height, which is a design choice constraint to: 16 The impeller geometry is calculated by a simple model where the blade camber line is assumed circular. More details of the method are provided in (Aungier, 2000). One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 116 𝑏4≤𝑏3 [Equation 4.96] Figure 101 Effect of the diffuser leading angle and the Mach number at diffuser inlet onto the stall incidence (blue) and flow incidence angle (red) for a diffuser with 20 vanes. Empirical observation (Aungier, 2000) shows that the most favourable diffuser leading angle designs are found for 𝛽4 between 68 and 74 º and for the radius ratio, 𝑟4/𝑟3, between 1.06 and 1.12. This region is marked with dashed lines in Figure 101. The design loop is presented below: 1. The meridional velocity at the diffuser leading edge (𝑣𝑚4) is guessed. 2. [Equation 4.95] is solved along with the momentum conservation equation (𝑟3𝑣𝑢3= 𝑟4𝑣𝑢4) and velocity diagrams, yielding the flow angle 𝛼4. An internal loop is necessary in order to solve this system of equations. 3. Mass and energy conservation is solved simultaneously with the equation of state and the total pressure loss equation 17 , yielding the static conditions. 4. Loss terms are calculated. In this initial estimation, only two loss terms are considered: one due to wall friction and another due to the geometric contraction at the diffuser leading edge (Benedict, et al., 1966). The latter is modelled as: 𝜔𝑐𝑜𝑛𝑡,𝑉𝐿=cos2𝛼3(1−𝑏3/𝑏4)2 [Equation 4.97] 5. In parallel, 𝑣𝑚4 is updated. 6. The geometry is fully defined when iterating from points 2 to 5. Then, this component is simulated in detail with the vaneless space model presented in section 4.2.4. 4.4.3.4 Vaned diffuser design Given that the leading edge is defined in previous design loops, designing the vaned diffuser comprises the calculation of the outlet and throat geometries only. Thus, there are three degrees of freedom that need to be specified for the full geometrical definition of the diffuser (Aungier, 2000). Thus, as opposed to other elements, the diffuser is designed by a sort of a trial-and-error process where, after defining those three degrees of freedom, the geometry is computed and the diffuser is modelled as described in section 4.2.5. Once this is completed, the parameters are checked to be within acceptable ranges and, if not, they are updated and the process is initiated again. Following the indications in (Aungier, 2000), the most commonly employed geometrical specifications are used here: radius ratio (𝑅=𝑟5/𝑟4), area ratio (𝐴𝑅) and number of blades (𝑧𝐷𝐼𝐹). 17 In the first iteration no losses are considered. 0 10 20 30 40 50 60 70 80 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 Stall flow angles (  s) and incidences (is) Vaned difuser leading edge angle, 4 is(M4=0.5) is(M4=0.6) is(M4=0.7) is(M4=0.8)  s(M4=0.5)  s(M4=0.6)  s(M4=0.7)  s(M4=0.8) One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 117 Amongst all the combinations of this parameters, the one that provides acceptable values for the divergence angle (2𝜃𝑐), diffuser loading (𝐿) and design parameter (𝐸), [Equation 4.61] is adopted as the final design. From a designer standpoint, the following ranges are recommended based on experience:  A divergence angle lower than 11º: 2𝜃𝑐≤11º  A loading parameter not higher than 1/3: 𝐿≤1/3  And a design parameter in the range 1.5-1.7, which indicates that using the vaned diffuser really improves the performance with respect to a vaneless diffuser: 1.5≤𝐸≤1.7. As opposed to the impeller, the blade camber line is not assumed circular. Instead, the method developed in (Aungier, 1988), and also described in (Aungier, 2000) is used. 4.4.3.5 Volute design There are two different methods to design the volute:  The simple area schedule (SAS) method is based on varying the cross sectional area linearly with the angular coordinate 𝜃. The design is accordingly defined by the full-collection plane area.  More common and more precise, the conservation of mass and angular momentum (CAM) approach considers area variations in such a way that these two conservation laws are satisfied. When angular momentum conservation (𝑟𝑣𝑢=𝑟5𝑣𝑢5) and continuity (𝜌𝐴𝑣= 𝜃𝜌5𝑟5𝑏5𝑣𝑚5) are satisfied, the cross-sectional area at a certain angle 𝜃 of the volute is: 𝐴= 𝜃𝑟5𝑏5/tan𝛼5. Therefore, the general equation for volute sizing is: 𝐴=𝑆𝑃 𝜃 𝑟5 𝑏5/tan𝛼5 [Equation 4.98] Where 𝑆𝑃 is the sizing parameter that controls the design, already used in section 4.2.6. It is desirable that the sizing parameter be in the range [1, 1.2], 1 being the optimum value which means that momentum is conserved. Taking 𝑆𝑃 as the input, the second approach is implemented in this work as follows: 1. The radius of the full-collection plane is initially guessed, 𝑅𝑓𝑐𝑝. 2. The geometry of the volute is completed with the simple expressions: {𝑟6=𝑟5+𝑅𝑓𝑐𝑝 𝐴𝑓𝑐𝑝=𝜋𝑅𝑓𝑐𝑝 [Equation 4.99] 3. The volute is modelled as described in section 4.2.6. 4. The sizing parameter is calculated by [Equation 4.84]. 5. 𝑅𝑓𝑐𝑝 is modified if 𝑆𝑃 has a value different to the desired one. 4.4.3.6 Design overview and results The content dedicated to stage design and presented in previous sections is finally summarised here. figure 102 and Table 16 present the result of a valid design for the first compression stage with 𝑁𝑠 and 𝐷𝑠 equal 0.6 and 4 respectively, which has been obtained by the methodology developed herein. A sketch of this impeller is represented in figure 103. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 118 Figure 102. First stage geometry obtained with the one-dimensional design tool. IMPELLER GEOMETRY Number of full blades 20 Number of splitter blades 20 Impeller inlet radius at hub 25 mm Impeller inlet radius at shroud 50 mm Impeller exit radius 103.6 mm Blade angle of the impeller leading edge at hub 33.0 º Blade angle of the impeller leading edge at mean radius 44.3 º Blade angle of the impeller leading edge at shroud 52.4 º Blade angle of the impeller trailing edge -1.49 º (backward) Angle between streamlines and shaft at impeller inlet 0 º Angle between streamlines and shaft at impeller exit 90 º Full blade length 72.5 mm Splitter blade length 36.3 mm Axial length of the impeller 66.1 mm Blade thickness at impeller leading edge 1 mm Blade thickness at impeller trailing edge 1 mm Blade height at impeller leading edge 6.9 mm Clearance gap at impeller tip 2 % of the span Table 16. Geometry of the impeller obtained for the S-CO2 cycle design. One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 119 Figure 103. Loop for stage design. 4.5 Summary and conclusions A literature review in the field of centrifugal compressor design has confirmed the importance of using one-dimensional tools in the preliminary stages of this process. After this evidence, a 1D tool has been developed which is able to provide a tentative detailed design of all the elements in such a turbomachinery: inducer, impeller, vaneless space, vaned diffuser and volute. Data input Initial Guess: ηtt, Kdiff-vol 1. Inducer design 1.1. Guess r2s. 1.2. Geometry of (2). 1.3. Velocity triangle at (2). 1.4. Throat geometry. 1.5. Throat flow àMth’. 1.6. Look for the desing point. 2. Inducer simulation 3. Impller desired conditions 3.1. h03 = f(ηtt) 3.2. P03 = f(Kdiff-vol) 4. Velocity triangle at impeller exit 4.1. Euler equation à vu3 4.2. Design choice à vm3 5. Impellerer design 5.1. Guess σ , λ , Σϖ . 5.2. Velocity triangle at (3) à β3 . 5.3. Static conditions at (3). 5.4. Blade height b3 by continuity. 5.5. Whole impeller geometry. 5.6. Impeller simulation. 5.7. Updating σ , λ , Σϖ . 5.8. From 5.2 to 5.7 till convergence. 6. Updating ηtt for the same Kdiff-vol Converged ηtt? INDUCER DESIGNIMPELLER DESIGN VANELESS SPACE DESIGN NO 7. Vaneless space design 7.1. Guess vm4 . 7.2. Eqn. 4.95 and momentum conservationà α4 . 7.3. Static conditions at (4). 7.4. Losses estimation. 7.5. Updating vm4 . 7.6. From 7.2 to 7.5 till convergence. 8. Vaneless space simulation DIFFUSER DESIGN 9. Diffuser geometry calculation 10. Diffuser simulation YES Acceptable design? NO YES 11. Volute design 11.1. Guess Rfcp . 11.2. Volute geometry. 11.3. Volute simulation. 11.4. Compute SP . 11.5. From 11.1 to 11.4 till gettting SP desired. VOLUTE DESIGN 12. Updating ηtt, Kdiff-vol Converged ηtt, Kdiff-vol? End YES NO One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 120 This 1D model of performance is presented in section 4.2 and its application to the direct problem is validated in section 4.3, both for supercritical carbon dioxide and air. The centrifugal compressor existing in the experimental S-CO2 loop at SANDIA NL was taken as reference for the first validation, whilst the NACA CC3 unit operating at the Glenn Research Center was considered for air. Good agreement was obtained for both cases and thus a tool based on the similar principles was implemented for the inverse problem (producing the geometry that yields a certain performance with the highest possible efficiency). This design process for the inverse problem is affected by two aspects that are original for the application under analysis. These features must be bear in mind as they do not apply in standard turbomachinery using air:  The Acceleration Margin to Condensation (AMC) becomes an important design parameters, in particular for the first stages. It is a way of measuring how far the flow is from meeting the saturation line (i.e. where condensation would take place theoretically), which would incur additional losses in addition to mechanical damage due to the impact of liquid droplets on the walls of the flow passage.  Since high efficiency is desired, the compression system usually comprises multiple stages because the rotating speed of the shaft is limited by the AMC. A high enthalpy rise requires high rotating speeds, which would violate the AMC constraint. Therefore, the total pressure rise must be split into several stages. Finally, the geometry of the first compression stage is presented at the end of the section 4.4 and will be analysed in the next chapter. 4.6 Justification of the research Once the design process together with all related concepts have been presented, it is time to vindicate the research project presented in this dissertation. It must be acknowledged first of all that even though supercritical carbon dioxide is an original working fluid, the corresponding compressor raising its pressure is still a piece of turbomachinery and, as such, its design can be approached following conventional methodologies. Nevertheless, this being said, it must also be bear in mind that the unconventional behaviour of carbon dioxide in the vicinity of the supercritical point brings about some special features into this process. This introduces substantial differences when designing these machineries with respect to standard ones. On one hand, regarding the operation of the exiting test rigs, there exist some operational constraints derived from the fluid nature. Thus, the condensation of the supercritical fluid in the inducer one represents an important inconvenient, since the experimentation with some existent prototypes has proven physical damage in the inducer blades when operating inside the biphasic region (below the saturation dome). Additionally, given the lower speed of sound of S-CO2, choke is encountered earlier (at a lower fluid velocity) than in air turbomachinery. These aspects impose restrictions onto shaft speed, which generally imply larger impeller diameter for the same pressure ratio (in order to achieve high enough blade speeds at impeller outlet). However, in the light of the results presented in chapter 3, the risk of stall in the diffuser channels is higher for S-CO2, unless the (dimensionless) length of the channel is also increased. A compromise must be thus achieved between these and other opposing effects. In summary, the design space of a supercritical carbon dioxide compressor is modified because of the fluid nature. Moreover, the performance of this turbomachinery is expected to vary with respect to similar equipment operating on air or other conventional fluid. A sample of this is presented below in order to make it evidence to what extent the design space is influenced by the nature of the fluid. Let the pressure ratio of a single stage radial compressor be considered as a function of the reduced pressure and temperature at the inlet (i.e. the thermodynamic distance from the critical point). The following assumptions are made: One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 121  The ranges considered for the reduced inlet conditions are equivalent to assuming compressibility factors in the range from 0.24 to 0.87 approximately.  The ratio from throat Mach number to AMC is set to 0.8 in order to avoid that avoid condensation during the initial expansion.  The mass flow rate is kept constant at 73.04 kg/s whereas the specific speed and diameter are set to 0.7 and 4 respectively in order to increase efficiency as much as possible. Figure 104. Effect of reduced inlet conditions on stage pressure ratio. The results shown in figure 104 confirm the following:  An anomalous behaviour is observed for reduced pressures over 1.05 when the reduced temperature gets close to unity. This is due to the sharp density rise with respect to temperature (see figure 16). Based on the definition of specific speed itself, the enthalpy rise decreases when density increases (Δℎ0𝑠∝𝜌−2/3), and pressure is related with it. These points are characterised by very low values of the compressibility factor (around 0.25).  The previous behaviour is not observed for the remaining curves as reducing the inlet temperature would give place to equilibrium condensation (coexistence of liquid and vapour in the subcritical region) rather than to a fluid with transition properties.  As a main conclusion, it is observed that effect of the supercritical nature of carbon dioxide is translated into a limitation of the stage pressure ratio down to values around 1.2-1.7. These pressure ratios are much lower than expected from ideal gas compressors and they tend to increase with inlet temperature even if this influence vanishes as temperature increases. As a matter of fact, there is no further increase in pressure ratio for reduced temperatures higher than 2 regardless of the reduced pressure at the inlet. For higher temperatures, pressure ratio is limited by choke in the throat.  Having into account that reducing the pressure ratio affects the specific work of the cycle, it would be necessary to increase the mass flow rate or, a rather elegant alternative, to employ multistage compressor, incurring in higher costs. From a thermoeconomic standpoint, it is concluded that the working fluid sets new boundaries for the design space where the designer must make a decision between opposing effects. On one hand, using single stage compressors imply a reduction on cycle pressure ratio which impacts the performance of the system negatively. This translates into lower efficiency (regardless of power output) and higher costs for a given output (larger equipment for larger flow rates). On the other hand, should the designer 1 1.5 2 2.5 3 3.5 4 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 Reduced Stagnation Temperature, T0,r PR1stage P0,r = 1.02P0,r = 1.02 P0,r = 1.05P0,r = 1.05 P0,r = 0.98P0,r = 0.98 P0,r = 0.95P0,r = 0.95 P0,r = 0.90P0,r = 0.90 P0,r = 1.10P0,r = 1.10 P0,r = 1.20P0,r = 1.20 One-dimensional model of supercritical carbon dioxide centrifugal compressors Design of supercritical carbon dioxide centrifugal compressors 122 target higher efficiencies, it would be necessary to employ multiple stage turbomachinery thus incurring much higher costs. This very simple example is just to confirm the interest of this research, which is deemed relevant for the S-CO2 community. Throughout the document, the main differences between air and supercritical carbon dioxide are discussed and design recommendations are provided whenever possible. Furthermore, the tool presented is expected to enable the development of an optimisation tool based on which design engineers will hopefully save time and money when developing turbomachinery for a particular application. Even if this tool exceeds the scope of this work, this research paves the way for its successful development in the near future. Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 123 CHAPTER 5. Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Before proceeding to prototype manufacturing, it is absolutely necessary to perform accurate aerodynamic simulations in order to ensure that the compressor will integrate correctly into the power cycle. One-dimensional tools like those presented in the previous chapter are useful for the rapid evaluation and global geometry definition of a number of design alternatives, but they provide limited information of the flow features within the machine stations and thus fully three-dimensional calculations are required. As illustrated in figure 72, it is also usual to employ simulations of intermediate complexity for optimisation, streamline curvature methods being the most widespread (Boyer, 2001) followed by inviscid-3D and turbulent-3D simulations. This chapter is aimed at fully-3D simulations of the impeller that resulted from the 1D design. It means that the exact geometry of the impeller is run by CFD (using ANSYS Fluent®) considering the particular features of supercritical carbon dioxide exposed in chapter 2, including boundary layer effects. Prior to that, a coarser impeller model is run, which excludes those elements that add considerable computational burden, i.e. clearance gap and boundary layer refinement. 2D and inviscid-3D simulations are hence omitted in this project even if they are mentioned at the end of the chapter as the subject of future research at GMTS. These 2D simulations are faster than 3D CFD models and with acceptable precision, thus being excellent choices when developing optimisation studies. In fact, they are catalogued as the “backbone of turbomachinery design” by Denton and Dawes (Denton & Dawes, 1999). 5.1 Geometry generation One-dimensional design provides blade angles and radii at the hub and shroud for the inlet and outlet stations of each compressor element, whereas other dimensions such as blade lengths and impeller axial length are estimated based on simple considerations proposed in (Aungier, 2000). It is worth noting that generating the full three-dimensional shape of the impeller and other elements requires extensive work, involving complex mathematical treatment of the geometry. It is nevertheless possible today to circumvent this hurdle thanks to specialised software such as the one employed in this project, ANSYS BladeGen®; which is able to create the 3D geometry apart from the parameters that follow:  Hub and shroud radii at inlet: 𝑟2ℎ and 𝑟2𝑠.  Radius and blade height at impeller exit: 𝑟3 and 𝑏3.  Impeller axial length: 𝐿𝑎𝑥𝑖𝑎𝑙  Number of full and splitter blades: 𝑧𝐹𝐵 and 𝑧𝑆𝐵  Blade angles at hub, mean-radius and shroud in the inlet station, as well as the blade angle at impeller exit: 𝛽2ℎ, 𝛽2, 𝛽2𝑠 and 𝛽3.  Length ratio between splitter and full blade: 𝐿𝑆𝐵/𝐿𝐹𝐵.  Thickness at leading and trailing edge: 𝑡2 and 𝑡3. Bezier curves (Casey, 1983) are the most common option employed to define the contours of hub and shroud and is here used (this is possible in ANSYS BladeGen®), Figure 105. The inlet flow is assumed in the axial direction whereas the outflow is contained in a plane normal to the shaft (purely radial machine), after which the default geometry provided by the software is assumed valid. This Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 130 Figure 111. Representation of rough wall model in ANSYS Fluent® (ANSYS, 2011). 𝐾𝑆+=𝜌 𝐾𝑆 𝑢∗ 𝜇 [Equation 5.4] Figure 112. Wall y+ obtained in a full impeller model with smooth (a) and rough (b) wall boundary conditions. Roughness grain height 25 μ m. 5.2.2 Other numerical aspects of the computational model In addition to the main remarks given in chapter 2 for the simulation of S-CO2 flows, this section presents the features of the numerical model that is specific for the compressor model. First of all, a Single Reference Frame (SRF) approach is used, meaning that the fluid domain is linked to a relative reference frame to which a rotational speed is imparted. Consequently, wall boundaries are set to moving rather than stationary, as required by ANSYS Fluent® when a frame motion is used. Moreover, zero velocity with respect to the relative reference frame is set for blades and hub, whilst for the shroud zero velocity is imposed with respect to the absolute frame. Finally, wall boundaries are completely defined with a mean roughness grain height of 25 μm and adiabatic condition. Inflow and outflow boundaries are defined as mass flow inlet and pressure outlet respectively, default choices for compressible flow. Moreover, the mass flow rate at the inlet boundary is set to 𝑚󰇗/𝑧𝐷𝐼𝐹 in a direction normal to the inlet boundary as only one channel of the impeller is simulated. The definition of this inlet flow is completed by the static inlet temperature obtained from the 1D calculation and the following turbulence parameters: Turbulence intensity = 1.92%, Turbulence length-scale = 7 mm 19 . In contrast, only the static pressure is required at the outlet, even though other static properties like temperature and turbulence must be specified in the eventual case of backflow (reverse flow). 19 These turbulence parameters have been set based on ANSYS Fluent® recommendations [ANSYS UserGuide]: Tu= 0.16Re−1/8 and lTu=0.07dH. y+ MAX (a) Blade leading edge 1.116 Blade pressure face 0.900 Blade suction face 0.757 Blade trailing edge 0.532 Blade tip 0.811 Hub 0.780 Shroud 0.829 Splitter leading edge 1.134 Splitter pressure face 0.858 Splitter suction face 0.639 Splitter trailing edge 0.554 Splitter tip 0.699 y+ MAX (b) Blade leading edge 1461 Blade pressure face 1189 Blade suction face 850 Blade trailing edge 620 Blade tip 1198 Hub 981 Shroud 1058 Splitter leading edge 1486 Splitter pressure face 995 Splitter suction face 764 Splitter trailing ege 614 Splitter tip 978 Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 131 Periodic contours close the fluid domain externally but there are other boundary conditions that need to be defined inside: the interfaces. In effect, when clearance gap is modelled, there exist surfaces that separate the pressure and suction sides of the blades and require the corresponding mesh interfaces. Figure 113. Impeller mesh including boundary layer refinement and clearance gap modelling. The Coupled scheme is employed for the pressure-velocity coupling algorithm whereas first order upwind interpolation schemes have been used for each variable. Other choices regarding pressurevelocity coupling algorithms have been unsuccessfully tested, arriving at non-converged solutions. On the other hand, fou schemes have been used because of its computational robustness, since higher orders increase the elapsed times considerably. As mentioned before, the final solution is obtained following a step-by-step methodology where the mean mesh size is fixed at the reference value (Mesh 3 in the GCI analysis) and the boundary layer is progressively refined until wall y+ lower or near 1 is achieved. Each mesh is initialised with the results from the immediately coarser mesh except for the first one, for which the Hybrid Initialization method implemented in ANSYS Fluent® is applied. This method requires both inlet and outlet pressure boundaries. Therefore, a total pressure inlet equal 75 bar was set up only for the hybrid initialization process, later changed to perform the simulation with the mass flow inlet boundary. Once initialized, simulations were carried out as follows: 1. Minimum pressure and temperature are limited to values above the critical point, 75 bar and 305 K respectively. 2. Simulations are started with flow equations only, i.e. continuity and momentum. 3. When the residuals show a clearly descending pattern, the Turbulence equations are activated. 4. The energy equation is activated when the difference of mass flow rates between inlet and outlet is considerably low: ~10−3 kg/s. 5. At this point, there exist a large number of cells with their pressure and temperature set to the minimum values imposed in 1, but this number of cells starts to decrease at a certain point in the iterative process. When the minimum number of limited cells, both for pressure and temperature, is nearly achieved, the aforesaid limits are modified to 30 bar and 200 K, values of pressure and temperature that must not be achieved in the computational models if they are accurate enough. Mass Flow Inlet Pressure Outlet Periodic Moving walls Mesh Detail: Boundary layer Mesh Detail: Clearance gap Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 132 5.3 CFD results Once the uncertainties associated with the mesh are, if not removed, at least, reduced and known, the numerical model of the impeller is calibrated with a main purpose: defining general guidelines to be applied in future simulations of S-CO2 turbomachinery. Therefore, the two main inaccuracies committed in the model employed for the GCI analyses should be now quantitatively evaluated. These are: (i) the boundary layer refinement and (ii) the clearance gap modelling, both of them previously omitted. Accordingly, the following analysis of four cases in parallel has been done:  Case 1: Without gap modelling and without boundary layer refinement.  Case 2: With gap modelling but without boundary layer refinement.  Case 3: Without gap modelling but with boundary layer refinement  Case 4: With gap modelling and with boundary layer refinement Four different meshes are generated in TurboGrid® with a common mean size (that of Mesh 3 in the CGI analysis, section 5.2.1.1)and the boundary layer is refined by reducing the y+ offset down to 1 where necessary. Regarding turbulence models and according to the recommendation in chapter 2, two different alternatives are used: 𝑘−𝜔 SST and 𝑘−𝜀 models depending on whether the boundary layer refinement is implemented or not. Regarding the inflow/outflow boundary conditions, these are just the same for all four cases at expected. Table 20 summarises the mass-weighted averaged values obtained in this 4-case analysis performed with ANSYS Fluent®. The values reported are measured at the inlet and outlet sections, the static pressure at the outlet being omitted for it is a boundary condition. Case 1 2 3 4 𝑇𝑖𝑛 (K) 309.60 309.69 309.51 309.65 𝑃𝑖𝑛 (bar) 71.17 71.61 70.66 71.44 𝑇𝑜𝑢𝑡 (K) 336.53 336.00 337.08 336.06 𝑣𝑖𝑛 (m/s) 55.99 55.22 56.88 55.53 𝑣𝑜𝑢𝑡 (m/s) 161.49 163.06 161.44 163.85 𝑤𝑖𝑛 (m/s) 83.77 83.27 84.37 83.47 𝑤𝑜𝑢𝑡 (m/s) 74.29 74.90 72.07 76.30 𝑣𝑚,𝑖𝑛 (m/s) 55.88 55.12 56.78 55.43 𝑣𝑚,𝑜𝑢𝑡 (m/s) 68.28 69.45 66.42 70.67 𝜌𝑖𝑛 (kg/m3) 221.43 224.49 217.96 223.27 𝜌𝑜𝑢𝑡 (kg/m3) 276.49 278.87 273.98 278.77 Table 20. Mass-weighted averaged values for the four-case analysis. Apart from the reported values, it is interesting to define the impeller operation by means of its main performance parameters (see chapter 4): total pressure loss coefficient (∑𝜛𝑖2→3 ), total-to-total isentropic efficiency of the impeller (𝜂𝑡𝑡,𝑖𝑚𝑝) and slip (𝜎) and distortion (𝜆) factors. The simple model presented below is been created for this analysis:  The static thermodynamic properties are calculated with the averaged temperatures and pressures directly: enthalpy ℎ=ℎ(𝑇,𝑃) and entropy 𝑠=𝑠(𝑇,𝑃).  The absolute and relative stagnation enthalpies are then calculated with the absolute and relative velocities: ℎ0=ℎ+𝑣2/2 and ℎ0′=ℎ+𝑤2/2; the corresponding entropies remain the same: 𝑠=𝑠0=𝑠0′. Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 133  Stagnation temperatures and pressures are calculated with the corresponding enthalpyentropy pairs. These values are presented in table 21, where the inlet total temperature (𝑇0,𝑖𝑛) is the same in all cases (it is a boundary condition). Based on these variables/properties, the aforelisted performance parameters can be calculated from their definitions directly:  The total-to-total isentropic efficiency of the impeller is calculated as: 𝜂𝑡𝑡,𝑖𝑚𝑝=ℎ0,𝑜𝑢𝑡,𝑠−ℎ0,𝑖𝑛 ℎ0,𝑜𝑢𝑡−ℎ0,𝑖𝑛 [Equation 5.5] Where the isentropic outlet enthalpy (ℎ0,𝑜𝑢𝑡,𝑠) is obtained for the actual total pressure at the outlet and the inlet entropy: ℎ0,𝑜𝑢𝑡,𝑠=ℎ(𝑃0,𝑜𝑢𝑡,𝑠0,𝑖𝑛).  Calculating the distortion factor requires the aerodynamic blockage factor (𝜆= 1(1−𝐵)⁄ ), which is computed with the following equation directly: 𝑚󰇗=𝜌 𝑣𝑚 𝐴 (1−𝐵) [Equation 5.6]  The slip factor results from the resolution of both the actual (measured) and ideal velocity diagrams at the impeller exit, both sharing the same meridional velocity and differing in the absolute and relative flow angles (note that absolute flow angle and blade angle are the same in the ideal diagram) (Figure 114). Figure 114. Ideal (transparent) vs. actual/measured velocity diagrams at impeller exit.  Finally the total pressure loss coefficient is calculated: 𝑃0,𝑜𝑢𝑡 ′=𝑃0,𝑜𝑢𝑡,𝑖𝑑 ′−𝑃0,𝑜𝑢𝑡 ′ 𝑃0,𝑖𝑛 ′(𝑃0,𝑖𝑛 ′−𝑃𝑖𝑛)∑𝜛𝑖 2→3 [Equation 5.7] where the ideal discharge has the same entropy as the inlet and the same total enthalpy as the outlet. The following conclusions are drawn from the results obtained from the model and presented in Table 21:  No large differences are observed between cases, neither in terms of averaged variables nor performance parameters, and this situation remains when statistical information is taken into account (mean 𝜇, variance 𝜎2 and standard deviation 𝜎 have been calculated for each variable in the group of four cases object of the study and included in Table 21). It w3 w3,id v3 v3,id u3 Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 134 can be observed that typical variations are 1% of mean values for all the variables except for the pressure loss coefficient and the distortion factor. The former has a descendent trend when the complexity of the model increases whilst the latter does not follow any apparent trend.  When it comes to considering whether to include the gap (case 2) or the boundary layer refinement (case 3), the first option is slightly preferred as its results are a bit more similar to the complete model in terms of performance parameter with a hardly noticeable exception for the pressure loss coefficient.  Including the gap seems to have an important effect on flow distortion. In particular, eliminating the gap from the complete model (case 3 respect to case 4) brings about a deviation higher than 7.5% in 𝜆, while this deviation downs to 1.6 % if the gap is modelled and the boundary layer refinement eliminated (case 2 respect to case 4). This fact suggests that case 3 exhibits a more pronounced difference with respect to case 4 than the one exhibited by case 2 in regard to the velocity profile at impeller exit. Case 1 2 3 4 𝜇 𝜎2 𝜎/𝜇 (%) 𝑇0,𝑖𝑛 (K) 313.15 313.15 313.15 313.15 313.15 0.0000 0 𝑃0,𝑖𝑛 (bar) 74.70 75.09 74.25 74.94 74.745 0.1010 0.425 𝑇0,𝑜𝑢𝑡 (K) 362.49 362.51 362.94 362.81 362.6875 0.0373 0.053 𝑃0,𝑜𝑢𝑡 (bar) 141.36 142.64 140.95 143.04 141.9975 0.7508 0.610 ∑𝜛𝑖 2→3 0.7282 0.7078 0.6932 0.5976 0.6817 0.0025 7.352 𝜂𝑡−𝑡 89.12 89.84 89.57 90.99 89.88 0.4768 0.768 𝜎 0.8414 0.8491 0.8486 0.8459 0.84625 9.32·10-6 0.361 𝜆 1.220 1.252 1.176 1.273 1.23025 0.0013 2.972 Table 21. Stagnation points and performance parameters for the four-case analysis. In order to thoroughly analyse this differences, a comparison in terms of velocity profiles at the impeller outlet is made. The outcome of this comparison is summarised in figure 115, yielding the following surprising observations: Figure 115. Absolute velocity profiles at impeller outlet for the four cases under analysis.  The most pronounced difference, both quantitatively and qualitatively, is found between cases 3 and 4. Looking at case 3 in detail, a velocity drop is observed in the wake downstream of the pressure side of the splitter which is not found in the other cases. This 158 160 162 164 166 168 170 172 174 0 0.2 0.4 0.6 0.8 1 Velocity magnitude (m/s) Pitchwise coordinate (-) CASE 1 CASE 2 CASE 3 CASE 4 Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 135 difference comes about because of the following numerical issue experienced by the cells that are adjacent to the blade tip and shroud in case 3 (where there is no interstitial gap in the mesh):: Figure 116. Detail of a cell sharing faces both with shroud and splitter surfaces. Figure 116 summarises the boundary conditions of the aforecited cell, namely cell A. On one hand, the velocity of cells A is influenced by the velocity imparted by the splitter wall, i.e. 𝑣𝑐𝑒𝑙𝑙 𝐴~𝑟·𝜔, whereas in the upper face of the cell the velocity is fixed to zero given the stationary wall boundary. Therefore, the velocity gradient in the tangential direction (which affects the wall shear stress) yields: 𝜏𝑤𝑎𝑙𝑙≅𝜇𝑣𝑐𝑒𝑙𝑙 𝐴 Δ𝑦𝐴/2 [Equation 5.8] On the other hand, the pressure and shear stress produced by the fluid (𝑃𝑓𝑙𝑢𝑖𝑑 and 𝜏𝑓𝑙𝑢𝑖𝑑) surrounding the wall-free sides of cell A are related by the following simplified equation as long as the tangential velocity of the fluid is assumed constant (stationary motion in the tangential coordinate): 𝜏𝑤𝑎𝑙𝑙+𝜏𝑓𝑙𝑢𝑖𝑑+𝑃𝑓𝑙𝑢𝑖𝑑Δ𝑦𝐴=𝑃𝑓𝑎𝑐𝑒Δ𝑦𝐴 [Equation 5.9] The main effect of a boundary layer refinement is that Δ𝑦𝐴 is largely reduced what in turns brings about a n inverse increase in 𝑃𝑓𝑎𝑐𝑒 to compensate for this effect. This is the reason why a pronounced pressure peak is observed in Figure 117. It must be noted that this numerical issue does not take place in case 1 where the gap is not present either. The reason is that Δ𝑦𝐴 is much larger due to the absence of mesh refinement. And neither does it occur in case 4 since type A cells are not found in the mesh due to the interstitial gap being modelled. In summary, it is the combination of boundary layer refinement and absence of interstitial gap which brings about the occurrence of the discussed numerical instability.  In consequence, where the boundary layer is refined, the impact of implementing the gap is more visible since the phenomena around the blades are solved for rather than being implemented in the turbulence model as an overall effect. The conclusion is thus that whenever the boundary layer is refined, the gap must be also implemented. On the contrary, if the boundary layer is not refined, it is less important whether or not the Shroud (stationary wall) Splitter (moving wall) vCell A τ wall τ fluid Pfluid Pface Cell A ΔyA Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 136 interstitial gap is also included. This is of course as far as local velocities are concerned. In terms of the global parameters, implementing the gap but not refining the boundary layer seems to be the most leveraged solution.  This explanation is further confirmed by the static pressure distribution on the blade walls (Figure 117), which show that the high pressures typically found on the blade tips are more pronounced in case 3 (see coloured scales on the left). This pressure peak in case 3 is behind the velocity drop of the wake in figure 115. Figure 117. Pressure contours for cases 1, 3 and 4, where pressure peaks on the blade tips can be observed in case 3.  The most difficult question to answer is whether or not the added computational burden of case 4 is worth the gain in accuracy. This is disputable and the best way to provide an answer is to ask ourselves the type of results on which we are interested. Thus, in terms of global parameters (averaged values), one would say that case 2 is most balanced. On the contrary, should the local flow phenomena at impeller outlet be of interest, then it would be necessary to select the most complex model including interstitial gap and boundary layer refinement.  In terms of turbulence model, it seems that the Realizable 𝑘−𝜀 model suffices to provide satisfactory results when no boundary layer refinement is implemented whereas the more complex 𝑘−𝜔 SST model must be used when the boundary layer is refined around the blades and in the wakes. Graphical comparisons between the four cases are additionally provided in the next figures, comprising contours of entropy, pressure and velocity, streamlines and velocity vectors. CASE 1 CASE 3 CASE 4 Peaks of pressure at blade tips Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 137 Figure 118. Contours of entropy at different meridional positions along the impeller channel. The previous figure presents those zones where flow irreversibilities are created, generally near the walls and in the gap, and more in particular in the shroud near the impeller exit. The next figure shows pressure contours that are very similar qualitatively, with the exception of case 3 because of the aforediscussed effect, which is again visible in Figure 117. Figure 119. Static pressure contours in hub, blades and splitters. Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 138 The next figures do not reveal any appreciable difference between the four cases, as reported in the previous quantitative analysis. Figure 120. Absolute velocity contours on a surface located at midspan. Figure 121. Relative velocity contours on a surface located at midspan. Three-dimensional CFD simulation of supercritical carbon dioxide turbomachinery Design of supercritical carbon dioxide centrifugal compressors 139 Figure 122. Streamlines starting from the inlet boundary coloured by absolute velocity. Figure 123. Vectors of absolute velocity on a surface located at midspan.