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

Learn More →

Optimal Design of a Ljungström Turbine for ORC Power Plants: From a 2D model to a 3D CFD Validation

Optimal Design of a Ljungström Turbine for ORC Power Plants: From a 2D model to a 3D CFD Validation International Journal of Turbomachinery Propulsion and Power Article Optimal Design of a Ljungström Turbine for ORC Power Plants: From a 2D model to a 3D CFD Validation Umberto Coronetta * and Enrico Sciubba Department of Mechanical and Aerospace Engineering, University of Roma Sapienza, Via Eudossiana 18, 00184 Roma, Italy; enrico.sciubba@uniroma1.it * Correspondence: umberto.coronetta@gmail.com Received: 23 April 2020; Accepted: 14 July 2020; Published: 20 July 2020 Abstract: In the last few years, waste-energy recovery systems based on the Organic Rankine Cycle (ORC) have gained increased attention in the global energy market as a versatile and sustainable technology for thermo-electric energy conversion from low-to-medium temperature sources, up to 350 C. For a long time, water has been the only working fluid commercially adopted in powerplants: axial and, for smaller machines, radial inflow turbines have been the preferred expanders since their gulp capacity matches the -T curve of water steam. The density of most organic compounds displays extremely large variations during the expansion (and the volume flow rate correspondingly increases along the machine channels), so that Radial Outflow Turbines (ROTs) have been recently considered instead of traditional solutions. This work proposes a two-dimensional inviscid model for the stage optimization of a counter-rotating ROT, known as the Ljungström turbine. The study starts by considering five di erent working fluids that satisfy both the gulp requirements of the turbine and the hot source characteristics. On the basis of a limited number of geometric assumptions and for a fixed set of operating conditions, di erent kinematic parameters are optimized to obtain the most ecient cascade configuration. Moreover, as shown in the conclusions, the most ecient blade profile leads to higher friction losses, making further investigation regarding the best configuration necessary. Keywords: turbine CFD; Ljungström turbine; Organic Rankine Cycle 1. Introduction The research and development of the Organic Rankine Cycle (ORC), based energy conversion systems, is mainly aimed at the accurate prediction of the performance of the expander, which strongly influences the eciency of the entire plant. Moreover, based on the characteristics of energy sources exploited by the plant, some expanders are more suitable than others, leading to a higher eciency of the overall process. With the exception of the smallest ORC power plants (few kW), where volumetric expanders of scroll and screw type can be used, whenever higher power outputs are needed, and consequently high expansion ratios and/or high volume flow rates are necessary, dynamic turbines are the only feasible expanders [1]. Tocci et al. [2], in their review of ORC for small scale applications, argue that the selection of the expander type in the range of power 20–70 kW represents an open question: on the one side, the size of volumetric expanders increases exponentially with an increment of output power [2], on the other, traditional radial inflow solutions may lead to unrealistic rotational speeds [3,4]. In light of the above, the purpose of this paper is to investigate the suitability of a special class of Ljungström turbines, specifically designed for the very large expansion rates of organic fluids in the range of power for which a preferred choice is not yet defined and consolidated. To do so, we will start Int. J. Turbomach. Propuls. Power 2020, 5, 19; doi:10.3390/ijtpp5030019 www.mdpi.com/journal/ijtpp Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 2 of 16 In light of the above, the purpose of this paper is to investigate the suitability of a special class of Ljungström turbines, specifically designed for the very large expansion rates of organic fluids in Int. J. Turbomach. Propuls. Power 2020, 5, 19 2 of 17 the range of power for which a preferred choice is not yet defined and consolidated. To do so, we will start from the study of the velocity triangles of a Ljungström type radial outflow turbine that exploit from the s low study -T tof hermal sou the velocity rces and pro triangles of duces s a Ljungström ignificantype t power out radial outflow put [5] and t turbine hen cont that exploits inue tolow-T their opti thermal mizasour tion. ces and produces significant power output [5] and then continue to their optimization. The Ljun The Ljungström gström tturbine urbine o owes wes it its s n name ame tto o it its s in inve vent ntors, ors, the Ljun the Ljungström gström brothers, who brothers, who in in 1908 1908 conceived, d conceived, designed esigned and pr and pr oduced oduced a stat a statorless orless steam t steam urbine, turbine, as sho as shown wn in Fig in u Figur re 1, cons e 1, consisting isting of twof o count two counter er-rota-r ting co otating axi coaxial al disks disks fitted w fitteditwith h seve several ral con concentric centric bla blade de crowns crowns [6[]. 6 ].The in The inlet let sect section ion is is located located rad radially ially ne near ar the sh the shaft aft while the while the outlet is at outlet is at th the e outer r outer radius adius at the exit at the exit of the last blade of the last bladecrown. crown. This arrangement compensates for the increase in the volume of the fluid during the expansion by a This arrangement compensates for the increase in the volume of the fluid during the expansion by a correspondin corresponding g incr increase ease in t in the he flow flow are area a along along the flo the flow w path. path. Figure 1. Ljungström Radially Outward Turbine (ROT): Fixed and Counter-rotating arrangements [7,8]. Figure 1. Ljungström Radially Outward Turbine (ROT): Fixed and Counter-rotating arrangements [7,8]. 2. Materials and Methods 2. Materials and Methods 2.1. Fluid Selection In this work, a hydro-treated mineral oil, Paratherm NF, is the heat transfer medium taken as the 2.1. Fluid Selection hot source. It enters the boiler at 140 C (413 K) with a mass flow rate of 4 kg/s. The cold source is water, In this work, a hydro-treated mineral oil, Paratherm NF, is the heat transfer medium taken as which enters the condenser at room temperature (25 C or 298 K) and is assumed to be available in the the hot source. It enters the boiler at 140 °C (413 K) with a mass flow rate of 4 kg/s. The cold source is necessary quantity. The selection of the working fluid and the identification of the thermodynamic water, which enters the condenser at room temperature (25 °C or 298 K) and is assumed to be cycle was performed via a MATLAB code by using the CoolProp library [9] to calculate the fluid available in the necessary quantity. The selection of the working fluid and the identification of the properties. A first selection of fluids from the CoolProp library was performed under two criteria: first, thermodynamic cycle was performed via a MATLAB code by using the CoolProp library [9] to since the source is at low temperature, the fluid should have a low critical pressure [10,11]; second, calculate the fluid properties. A first selection of fluids from the CoolProp library was performed the Ljungström turbine requires relatively high volume flow rates at its inlet and a well-identified under two criteria: first, since the source is at low temperature, the fluid should have a low critical variation in the density along the radius due to its peculiar cross-sectional area, and thus the fluid pressure [10,11]; second, the Ljungström turbine requires relatively high volume flow rates at its inlet ought to have a low density at turbine inlet (i.e., high values of the coecients in the Redlich–Kwong and a well-identified variation in the density along the radius due to its peculiar cross-sectional area, equation of state). According to these considerations the following fluids were considered: an alicyclic and thus the fluid ought to have a low density at turbine inlet (i.e., high values of the coefficients in hydrocarbon, Cyclopentane (C5H10); an azeotropic mixture, SES36; a haloalkane refrigerant, R134a the Redlich–Kwong equation of state). According to these considerations the following fluids were (CH2FCF3); an alkane refrigerant, R601 (C5H12); a hydrofluorocarbon refrigerant, R245fa (C3H3F5). considered: an alicyclic hydrocarbon, Cyclopentane (C5H10); an azeotropic mixture, SES36; a For each working fluid, as shown in Figure 2, a specific thermodynamic cycle was designed by haloalkane refrigerant, R134a (CH2FCF3); an alkane refrigerant, R601 (C5H12); a hydrofluorocarbon fixing the top and bottom pinch points and specifying a 10 C water temperature rise in the condenser. refrigerant, R245fa (C3H3F5). For a correct comparison, the thermal power input is the same for all cycles considered. The heat input For each working fluid, as shown in Figure 2, a specific thermodynamic cycle was designed by and output in the boiler and the condenser were then matched by iteratively adjusting the working fixing the top and bottom pinch points and specifying a 10 °C water temperature rise in the fluid mass flow rate or/and the evaporation temperature. Therefore, the turbine enthalpy drop and the condenser. For a correct comparison, the thermal power input is the same for all cycles considered. pump work were determined with an assumed 0.8 eciency for the pump and 0.83 (polytropic) for the The heat input and output in the boiler and the condenser were then matched by iteratively adjusting turbine. The overall cycle eciency is equal to: the working fluid mass flow rate or/and the evaporation temperature. Therefore, the turbine enthalpy drop and the pump work were determined with m  an L as sumed L 0.8 efficiency for the pump and 0.83 w f turb pump = , (1) cyc . (polytropic) for the turbine. The overall cycle efficiency is equal to: m  DH oil oil Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 3 of 16 Int. J. Turbomach. Propuls. Power 2020, 5, 19 3 of 17 = , (1) Figure 2. Temperature–entropy diagram for the fluids considered in this study fluids, in the same scale Figure 2. Temperature–entropy diagram for the fluids considered in this study fluids, in the same as the usual water Mollier diagram. scale as the usual water Mollier diagram. Tables 1–5 present the thermodynamic properties of the considered fluids at the relevant stations Tables 1–5 present the thermodynamic properties of the considered fluids at the relevant stations in the process. Figures 3–7 clearly show that, in all of the cycles, a regenerator recovers a portion of the in the process. Figures 3–7 clearly show that, in all of the cycles, a regenerator recovers a portion of sensible heat content of the working fluid before it enters the condenser. The percentage of thermal the sensible heat content of the working fluid before it enters the condenser. The percentage of power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Cyclopentane, thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the SES36. Eciencies Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the and mass flow rates (indicative of expander size) are reported in Table 6. SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6. Table 1. Cyclopentane thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 73,593 378,072 2.04 Pump (2) 313.00 73,593 18,013 725.38 Economizer (3) 313.06 250,885 17,768 725.52 Evaporator (4) 353.00 250,885 61,970 682.63 Superheater (5) 353.00 250,885 426,794 6.43 Turbine (6) 373.00 250,885 458,217 6.01 Regenerator (7) 341.04 73,593 416,033 1.86 at the inlet of each equipment. Table 2. R601 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Figure 3. Cyclopentane T–S diagram. Equipment [K] [Pa] [J/kg] [kg/m ] Table 1. Cyclopentane thermodynamic parameters. Condenser (1) 313.00 115,093 363,518 3.35 Pump (2) 313.00 115,093 9010 605.85 1 1 1 1 Temperature Pressure Enthalpy Density Economizer (3) 313.09 366,619 9425 606.13 Equipment [K] [Pa] [J/kg] [kg/m ] Evaporator (4) 353.00 366,619 108,917 562.19 Superheater (5) 353.00 366,619 427,177 10.11 Condenser (1) 313.00 73,593 378,072 2.04 Turbine (6) 373.00 366,619 468,490 9.35 Pump (2) 313.00 73,593 −18,013 725.38 Regenerator (7) 348.94 115,093 430,071 2.96 Economizer (3) 313.06 250,885 −17,768 725.52 at the inlet of each equipment. Evaporator (4) 353.00 250,885 61,970 682.63 Superheater (5) 353.00 250,885 426,794 6.43 Turbine (6) 373.00 250,885 458,217 6.01 Regenerator (7) 341.04 73,593 416,033 1.86 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2020, 5, 19 4 of 17 Table 3. R134a thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 1012,509 419,363 49.87 Pump (2) 313.00 1012,509 256,185 1147.37 Economizer (3) 314.11 2624,797 257,590 1155.48 Evaporator (4) 353.00 2624,797 322,105 929.39 Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 3 of 16 Superheater (5) 353.00 2624,797 428,830 154.37 Turbine (6) 373.00 2624,797 460,255 121.11 Regenerator (7) 333.93 1012,509 442,132 43.82 = , (1) at the inlet of each equipment. Table 4. R245fa thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 249,412 435,245 13.95 Pump (2) 313.00 249,412 252,838 1297.13 Economizer (3) 313.55 1648,462 253,916 1300.84 Evaporator (4) 385.44 1648,462 360,253 1038.29 Superheater (5) 385.44 1648,462 482,315 98.34 Turbine (6) 405.44 1648,462 509,011 84.42 Regenerator (7) 354.94 249,412 476,033 11.89 at the inlet of each equipment. Table 5. SES36 thermodynamic parameters. Figure 2. Temperature–entropy diagram for the fluid 1 s 1considered in this stu 1 dy fluids, in the 1 same Temperature Pressure Enthalpy Density Equipment scale as the usual water Mollier diagram. [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 116,523 389,772 8.78 Tables 1–5 present the thermodynamic properties of the considered fluids at the relevant stations Pump (2) 313.00 116,523 233,461 1333.94 in the process. Figures 3–7 clearly show that, in all of the cycles, a regenerator recovers a portion of Economizer (3) 313.31 857,341 234,016 1335.99 Evaporator (4) 386.78 857,341 327,876 1104.87 the sensible heat content of the working fluid before it enters the condenser. The percentage of Superheater (5) 386.78 857,341 442,402 63.85 thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Turbine (6) 406.78 857,341 466,484 56.91 Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the Regenerator (7) 370.35 116,523 439,563 7.21 SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6. at the inlet of each equipment. Figure 3. Cyclopentane T–S diagram. Figure 3. Cyclopentane T–S diagram. Table 1. Cyclopentane thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 73,593 378,072 2.04 Pump (2) 313.00 73,593 −18,013 725.38 Economizer (3) 313.06 250,885 −17,768 725.52 Evaporator (4) 353.00 250,885 61,970 682.63 Superheater (5) 353.00 250,885 426,794 6.43 Turbine (6) 373.00 250,885 458,217 6.01 Regenerator (7) 341.04 73,593 416,033 1.86 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 4 of 16 Int. J. Turbomach. Propuls. Power 2020, 5, 19 5 of 17 Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 4 of 16 Figure 4. R601 T–S diagram. Table 2. R601 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 115,093 363,518 3.35 Pump (2) 313.00 115,093 9,010 605.85 Economizer (3) 313.09 366,619 9,425 606.13 Evaporator (4) 353.00 366,619 108,917 562.19 Superheater (5) 353.00 366,619 427,177 10.11 Turbine (6) 373.00 366,619 468,490 9.35 Regenerator (7) 348.94 115,093 430,071 2.96 at the inlet of each equipment. Figure 4. R601 T–S diagram. Figure 4. R601 T–S diagram. Table 2. R601 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 115,093 363,518 3.35 Pump (2) 313.00 115,093 9,010 605.85 Economizer (3) 313.09 366,619 9,425 606.13 Evaporator (4) 353.00 366,619 108,917 562.19 Superheater (5) 353.00 366,619 427,177 10.11 Turbine (6) 373.00 366,619 468,490 9.35 Regenerator (7) 348.94 115,093 430,071 2.96 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 5 of 16 Figure 5. R134a T–S diagram. Figure 5. R134a T–S diagram. Table 3. R134a thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 1012,509 419,363 49.87 Pump (2) 313.00 1012,509 256,185 1147.37 Economizer (3) 314.11 2624,797 257,590 1155.48 Evaporator (4) 353.00 2624,797 322,105 929.39 Superheater (5) 353.00 2624,797 428,830 154.37 Turbine (6) 373.00 2624,797 460,255 121.11 Regenerator (7) 333.93 1012,509 442,132 43.82 at the inlet of each equipment. Figure 6. R245fa T–S diagram. Figure 5. R134a T–S diagram. Figure 6. R245fa T–S diagram. Table 3. R134a thermodynamic parameters. Table 4. R245fa thermodynamic parameters. 1 1 1 1 1 1 1 1 Temperature Pressure Enthalpy Density Temperature Pressure Enthalpy Density Equipment Equipment [K] [Pa] [J/kg] [kg/m ] [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 1012,509 419,363 49.87 Condenser (1) 313.00 249,412 435,245 13.95 Pump (2) 313.00 1012,509 256,185 1147.37 Pump (2) 313.00 249,412 252,838 1297.13 Economizer (3) 314.11 2624,797 257,590 1155.48 Economizer (3) 313.55 1648,462 253,916 1300.84 Evaporator (4) 353.00 2624,797 322,105 929.39 Evaporator (4) 385.44 1648,462 360,253 1038.29 Superheater (5) 353.00 2624,797 428,830 154.37 Superheater (5) 385.44 1648,462 482,315 98.34 Turbine (6) 373.00 2624,797 460,255 121.11 Turbine (6) 405.44 1648,462 509,011 84.42 Regenerator (7) 333.93 1012,509 442,132 43.82 Regenerator (7) 354.94 249,412 476,033 11.89 at the inlet of each equipment. at the inlet of each equipment. Figure 7. SES36 T–S diagram. Table 5. SES36 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 116,523 389,772 8.78 Pump (2) 313.00 116,523 233,461 1333.94 Economizer (3) 313.31 857,341 234,016 1335.99 Evaporator (4) 386.78 857,341 327,876 1104.87 Superheater (5) 386.78 857,341 442,402 63.85 Turbine (6) 406.78 857,341 466,484 56.91 Regenerator (7) 370.35 116,523 439,563 7.21 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 5 of 16 Figure 6. R245fa T–S diagram. Table 4. R245fa thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 249,412 435,245 13.95 Pump (2) 313.00 249,412 252,838 1297.13 Economizer (3) 313.55 1648,462 253,916 1300.84 Evaporator (4) 385.44 1648,462 360,253 1038.29 Superheater (5) 385.44 1648,462 482,315 98.34 Turbine (6) 405.44 1648,462 509,011 84.42 Regenerator (7) 354.94 249,412 476,033 11.89 Int. J. Turbomach. Propuls. Power 2020, 5, 19 6 of 17 at the inlet of each equipment. Figure 7. SES36 T–S diagram. Figure 7. SES36 T–S diagram. Table 6. Eciency, mass flow rate and expansion ratio for the considered fluids. Table 5. SES36 thermodynamic parameters. m ˙ 1 1 1 1 Fluid Temperature Pressure Enthalpy Density [–] [kg/s] [–] Equipment [K] [Pa] [J/kg] [kg/m ] Cyclopentane 9.63% 1.35 3.41 Condenser (1) 313.00 116,523 389,772 8.78 R601 9.79% 1.50 3.19 Pump (2) 313.00 116,523 233,461 1333.94 R134a 11.23% 4.03 2.59 R245fa 16.34% 1.50 6.61 Economizer (3) 313.31 857,341 234,016 1335.99 SES36 15.46% 1.50 7.36 Evaporator (4) 386.78 857,341 327,876 1104.87 Superheater (5) 386.78 857,341 442,402 63.85 2.2. The Optimization Algorithm Turbine (6) 406.78 857,341 466,484 56.91 Regenerator (7) 370.35 116,523 439,563 7.21 The standard formulation of the Euler work (2) shows that, for a Ljungström “stage”, the only at the inlet of each equipment. positive contribution to the specific work is given by the increase in the relative velocity, because u is out only slightly higher than u : in 2 2 2 2 2 2 v v + u u + w w in out in out out in L = , (2) Eul The novel method for the mean-line design proposed here begins with the following two positions, first suggested by Kearton [12] and Shepherd [13]: All blades, except those in the innermost ring, have the same cross-section, the same profile and the same stagger angle and, as a consequence, also have the same outlet angle: = const., (3) out In all blade crowns, the ratio of the relative outlet steam velocity to the peripheral velocity at the outlet edge of the blade ring is constant and equal to: out = < 1, (4) out In all blade crowns, the ratio of the radial chord of the blade to the outlet radius is constant and equal to: b r in = 1 = < 1, (5) r r out out Due to the axial symmetry of the turbine, the absolute velocity at the inlet of the first row can be considered completely radial. Int. J. Turbomach. Propuls. Power 2020, 5, 19 7 of 17 The above assumptions suce for a definition of kinematic eciency. Since the relative velocity provides the only positive contribution to the Euler ’s work, the absolute velocity and the peripheral velocity are to be considered as losses, and this leads to the following formula: 2 2 2 2 v v + u u out in out in = 1 , (6) kin 2 2 w w out in That can be rewritten as: 2 2 1   1 + 2 2 cos r r out = 1 , (7) kin,i 2 2 1  (1 + 4 4 cos ) r r out This equation provides additional useful information because it is similar for every i-th blade ring except for the first one, where it takes the following form: 2 + 1 2 cos r r out cos = 1 , (8) kin,1 cos The parameters to be considered to maximize the kinematic eciency are, therefore,  , , r out Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 7 of 16 and . Due to the second and third assumptions it is not possible to maximize both  and 0 kin,i kin ,1 with respect to  and , thus only the  equation is considered here for the optimization. r kin,i It can be easily verified that the derivative of the kinematic efficiency (7) is linearly dependent It can be easily verified that the derivative of the kinematic eciency (7) is linearly dependent on on both χ and βout and thus it is not possible to maximize this equation with respect to these both  and and thus it is not possible to maximize this equation with respect to these parameters. out parameters. This conclusion is reinforced by a closer inspection of Equation (7), which is maximized This conclusion is reinforced by a closer inspection of Equation (7), which is maximized either when either when the numerator of the second term tends to zero or the denominator to a maximum. Both the numerator of the second term tends to zero or the denominator to a maximum. Both cases leading cases leading to ideal turbine configurations: to ideal turbine configurations: • χ = 1, meets the maximum of the efficiency for ρr = cos(βout)/2, and causes the numerator to be = 1, meets the maximum of the eciency for  = cos( )/2, and causes the numerator to be r out null. In fact, in this case uin = uout and vin = vout so that the velocity triangle is similar to that of an null. In fact, in this case u = u and v = v so that the velocity triangle is similar to that of an out out in in impulse turbine; impulse turbine; • βout = 0 causes the denominator to be at its maximum, in which case win = 0 because the triangle = 0 causes the denominator to be at its maximum, in which case w = 0 because the triangle out collapses onto a line. in collapses onto a line. Therefore, the only parameter to be considered in the velocity triangles optimization is ρr, which allows us to Thereforma e, th ximiz e only e th parameter e kinematic to effi be consider ciency for diff ed in erent the velocity values of triangles βout and optimization χ, as shown in Fig is  , which ure 8. allows us to maximize the kinematic eciency for di erent values of and , as shown in Figure 8. out (a) (b) Figure 8. Kinematic eciency: (a) Variation of the kinematic eciency with  for di erent and Figure 8. Kinematic efficiency: (a) Variation of the kinematic efficiency with rρr for different outβout and fixed ; (b) Variation of the kinematic eciency with  for di erent  and fixed . fixed χ; (b) Variation of the kinematic efficiency with rρr for different χ and fixed outβout. Considering the above, we can now proceed to the study of the isentropic efficiency, starting from its equation and the loss model adopted here: (9) ℎ = , In order to obtain the value of wout,id, Soderberg [14] loss correlation was used: = 1 , (10) And after some rearrangement, the formula of the isentropic efficiency for the i-th row, with the exception of the first one, is obtained: = , (11) , , While for the first row: = , (12) The Soderberg loss coefficient is adopted for all the ensuing calculations, but different models were tested as well, such as Ainley and Mathieson [15] and Craig and Cox [16] prediction models. With little modifications in the code, the loss coefficient ξr can be evaluated for any of the Int. J. Turbomach. Propuls. Power 2020, 5, 19 8 of 17 Considering the above, we can now proceed to the study of the isentropic eciency, starting from its equation and the loss model adopted here: 2 2 2 2 w w u u out,id in out in Dh = , (9) ss In order to obtain the value of w , Soderberg [14] loss correlation was used: out,id w = w 1 +  , (10) out,id out r And after some rearrangement, the formula of the isentropic eciency for the i-th row, with the exception of the first one, is obtained: cos r r 2 out 2 1 + = ! , (11) is,i p 1+ +4 4 1+ cos 1+ r,i1 r r r,i1 out r,i 1 +  1 2 2 r r While for the first row: cos r out r =   , (12) is,1 1+ r,1 1 +  1 2 2 cos r 0 The Soderberg loss coecient is adopted for all the ensuing calculations, but di erent models were Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 8 of 16 tested as well, such as Ainley and Mathieson [15] and Craig and Cox [16] prediction models. With little modifications in the code, the loss coecient  can be evaluated for any of the abovementioned abovementioned prediction models. It was found that the adoption of different loss models does not prediction models. It was found that the adoption of di erent loss models does not noticeably influence noticeably influence the results, as shown in Figure 9. the results, as shown in Figure 9. (a) (b) Figure Figure 9. 9. Isentr Isentropic eff opic eiciency: ciency: ( (a a)) Variation Variation of of the i the isentr sentropi opicce eff ciency iciency wit withh ρr for for di differe erent nt βout and and r out fixed ; (b) Variation of the kinematic eciency with  for di erent  and fixed . fixed χ; (b) Variation of the kinematic efficiency with ρr for different χ and fixed βout. r out Since  depends on the fluid dynamic and geometric conditions that are di erent in each row, Since ξr depends on the fluid dynamic and geometric conditions that are different in each row, it is not possible to choose a single value of  for all rows that optimizes  . This justifies the choice r is,i it is not possible to choose a single value of ρr for all rows that optimizes ηis,i. This justifies the choice of a velocity triangles optimization based on the selection of the  that maximizes the  . of a velocity triangles optimization based on the selection of the ρr that maximizes the ηkin,i kin,i. As shown in Figure 10, for a fixed geometry (fixed  and ) and by varying the velocity, out As shown in Figure 10, for a fixed geometry (fixed χ and βout) and by varying the velocity, and and consequently  , the maximum of the isentropic eciency is close to the range of the maxima consequently ξr, the maximum of the isentropic efficiency is close to the range of the maxima of the kinematic efficiency. This supports the choice of an optimization procedure based on the velocity triangles. Figure 10. Comparing the variation of ηkin and ηis with ρr for fixed χ and βout and different flow velocity. A MATLAB code was written to obtain a complete mean line analysis, and then the flow in the resulting geometry was simulated via a CFD analysis to validate the model. The inputs of the algorithm are shown in Table 7. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 8 of 16 abovementioned prediction models. It was found that the adoption of different loss models does not noticeably influence the results, as shown in Figure 9. (a) (b) Figure 9. Isentropic efficiency: (a) Variation of the isentropic efficiency with ρr for different βout and fixed χ; (b) Variation of the kinematic efficiency with ρr for different χ and fixed βout. Since ξr depends on the fluid dynamic and geometric conditions that are different in each row, it is not possible to choose a single value of ρr for all rows that optimizes ηis,i. This justifies the choice of a velocity triangles optimization based on the selection of the ρr that maximizes the ηkin,i. Int. J. Turbomach. Propuls. Power 2020, 5, 19 9 of 17 As shown in Figure 10, for a fixed geometry (fixed χ and βout) and by varying the velocity, and consequently ξr, the maximum of the isentropic efficiency is close to the range of the maxima of the kinematic efficiency. This supports the choice of an optimization procedure based on the velocity of the kinematic eciency. This supports the choice of an optimization procedure based on the triangles. velocity triangles. Figure 10. Comparing the variation of  and  with  for fixed  and and di erent flow velocity. r out kin is Figure 10. Comparing the variation of ηkin and ηis with ρr for fixed χ and βout and different flow velocity. A MATLAB code was written to obtain a complete mean line analysis, and then the flow in the A MATLAB code was written to obtain a complete mean line analysis, and then the flow in the resulting geometry was simulated via a CFD analysis to validate the model. The inputs of the algorithm resulting geometry was simulated via a CFD analysis to validate the model. The inputs of the are shown in Table 7. algorithm are shown in Table 7. Table 7. Inputs and constraints of the algorithm. m [kg/s] Working Fluid Mass Flow Rate  [MPa] Allowable Stress at the Shaft 80 all p [Pa] Inlet pressure ! [rpm] Maximum speed 4400 max in p [Pa] Outlet pressure  [–] Blade encumbrance 0.85 out T [K] Inlet temperature l [m] Minimum blade height 0.01 in min b [m] Radial chord 0.01 min [deg] Minimum relative outlet angle 18 out By setting a trial rotational speed ! = ! the shaft diameter can be found that sets a lower max limit to the inlet radius at the entrance of the first row. Thence we obtain a matrix of all possible combinations of r and that satisfy the continuity equation for all the l > l . in in min The abovementioned matrix leads to another matrix of all possible values of  . By taking the kin maximum value we can also find the optimal values of  and . With all these parameters in place, it is now possible to calculate the loss coecient  and the enthalpy drop for the first stage (Dh ). n rows Dh = Dh , (13) i tot i = 1 If (13) is not satisfied it is possible to design additional rows and if an exact match is not achieved, a modified rotational speed ! = ! D! is selected and the process is iteratively repeated, as shown max in Figure 11. Since the maximum eciency is reached for b = 0 and = 0 and since these conditions are out clearly unphysical, a minimum value must be arbitrarily specified. In addition, the change of the cross-sectional area along the radius requires the designer to specify a minimum value of the axial blade length at the inlet: w f h = , (14) min 2!r  tan in b where the leakage factor  was assumed here to be 0.98. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 9 of 16 Table 7. Inputs and constraints of the algorithm. [kg/s] Working Fluid Mass Flow Rate τall [MPa] Allowable Stress at the Shaft 80 pin [Pa] Inlet pressure ωmax [rpm] Maximum speed 4400 pout [Pa] Outlet pressure δb [–] Blade encumbrance 0.85 Tin [K] Inlet temperature lmin [m] Minimum blade height 0.01 bmin [m] Radial chord 0.01 βout [deg] Minimum relative outlet angle 18° By setting a trial rotational speed ω = ωmax the shaft diameter can be found that sets a lower limit to the inlet radius at the entrance of the first row. Thence we obtain a matrix of all possible combinations of rin and βin that satisfy the continuity equation for all the l > lmin. The abovementioned matrix leads to another matrix of all possible values of ηkin. By taking the maximum value we can also find the optimal values of ρr and χ. With all these parameters in place, it is now possible to calculate the loss coefficient ξr and the enthalpy drop for the first stage (∆hi). ∑ ∆ℎ =∆ℎ , (13) If (13) is not satisfied it is possible to design additional rows and if an exact match is not achieved, a modified rotational speed ω = ωmax − ∆ω is selected and the process is iteratively repeated, as shown Int. J. Turbomach. Propuls. Power 2020, 5, 19 10 of 17 in Figure 11. Figure 11. Schematic representation of the design algorithm. Figure 11. Schematic representation of the design algorithm. The above equation has two degrees of freedom. By fixing a proper range of values for and r 0 in Since the maximum efficiency is reached for b = 0 and βout = 0 and since these conditions are all possible combinations that meet the equation can be calculated. These ranges are constrained by clearly unphysical, a minimum value must be arbitrarily specified. In addition, the change of the two di erent conditions. must be limited in order to obtain w > w : out, 0 1 in ,1 cross-sectional area along the radius requires the designer to specify a minimum value of the axial blade length at the inlet: 0 1 B C B 1 C B C B C 0 < < cos  ∙ , (15) 0 B r C @ 2 2 2 A ℎ = , (14) 1 +  2 cos + r r out r where the leakage factor ξ was assumed here to be 0.98. while r for the first row must be larger than the shaft radius that was calculated with the in following formula: 3 16P d  K , (16) sha f t where K = 1.3 is a safety factor, P [MPa] the power at the shaft and ! [rpm] its rotational speed [17]. Subsequently the code designs the first row using the inputs provided by the optimization function. It starts from the velocity triangles, calculates the isentropic heat drop using the loss coecients and finally derives the thermodynamic quantities at the outlet of the row. The quantities thus derived are used as inputs for the following stage, leading to an iterative process. The proposed algorithm leads to a choice in the number of stages as a function of the maximum angular speed. w f h = , (17) min 2!r  tan . in b As shown in Table 8 and, as expected, the  parameter is strictly correlated to the eciency of the turbine. A closer look to these results reveals that the less ecient thermodynamic cycles presents the higher isentropic eciencies. The Ljungström turbine performs better with low enthalpy drop and high-volume flow rates (higher specific speed n ). Moreover, since the cross-sectional area depends on the volume flow rate variation, not all the selected fluids lead to geometrical characteristics, such as Int. J. Turbomach. Propuls. Power 2020, 5, 19 11 of 17 axial blade length, which correspond to the allowed manufacturing limits. As shown in Table 9, the only suitable fluids are Cyclopentane and R601, the only ones leading to an axial blade length acceptable for the imposed constraints. Table 8. Crucial design parameters for di erent working fluids. ! Power Fluid n of Rows n kin is [rpm] [kW] Cyclopentane 6 0.587 3648 56 0.914 0.836 0.925 R601 6 0.586 4016 56 0.897 0.826 0.909 R134a 6 0.434 3714 61 0.815 0.729 0.821 R245fa 4 0.499 9171 34 0.698 0.590 0.688 SES36 6 0.350 4541 32 0.783 0.675 0.792 calculated at first stage. Table 9. Axial blade length (expressed in (m)) at the inlet of each row for di erent working fluids. Fluid 1st Row 2nd Row 3rd Row 4th Row 5th Row 6th Row Cyclopentane 0.010 0.011 0.011 0.011 0.011 0.012 R601 0.010 0.010 0.010 0.009 0.010 0.010 R134a 0.011 0.008 0.006 0.005 0.004 0.003 R245fa 0.010 0.006 0.003 0.003 SES36 0.010 0.007 0.005 0.004 0.004 0.004 In the present study, R601 was chosen as working fluid because, despite the slightly lower eciency of its cycle, it is much less flammable than Cyclopentane. Moreover, its condenser pressure is about 1 bar against 0.7 bars for Cyclopentane, making the manufacturing easier. The design parameters thus obtained, as shown in Table 10, make it possible to draw the velocity triangles of the best performing configuration, as shown below in Figure 12: Table 10. Output parameters for the selected fluid. Parameter Unit 1st Row 2nd Row 3rd Row 4th Row 5th Row 6th Row – 0.615 0.850 0.845 0.839 0.834 0.830 is n – 0.54 0.31 0.29 0.27 0.26 0.26 d – 3.72 4.77 5.13 5.43 5.66 5.77 u [m/s] 42 46 51 56 61 68 in u [m/s] 46 51 56 61 68 74 out v [m/s] 29 54 60 66 72 79 in v [m/s] 54 60 66 72 79 87 out w [m/s] 51 30 33 36 40 44 in w [m/s] 96 106 116 128 141 155 out – 0.7 0.64 0.64 0.64 0.64 0.64 in – 0.64 0.64 0.64 0.64 0.64 0.64 out – 0.00 0.98 0.98 0.98 0.98 0.98 in – 0.98 0.98 0.98 0.98 0.98 0.98 out h [m] 468,500 466,000 461,000 456,000 449,000 441,000 in h [m] 466,000 461,000 456,000 449,000 441,000 431,000 out T [K] 373 371 368 364 360 355 in T [K] 371 368 364 360 355 349 out s [J/kgK] 1,343 1,347 1,349 1,352 1,355 1,360 in p [Pa] 366,619 366,040 292,470 247,050 201,090 156,400 in p [Pa] 366,040 292,470 247,050 201,090 156,400 115,050 out [kg/m ] 9.35 8.54 7.42 6.27 5.11 3.99 in [kg/m ] 8.54 7.42 6.27 5.11 3.99 2.95 out Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 11 of 16 The design parameters thus obtained, as shown in Table 10, make it possible to draw the velocity triangles of the best performing configuration, as shown below in Figure 12: Table 10. Output parameters for the selected fluid. Parameter Unit 1st Row 2nd Row 3rd Row 4th Row 5th Row 6th Row ηis – 0.615 0.850 0.845 0.839 0.834 0.830 ns – 0.54 0.31 0.29 0.27 0.26 0.26 ds – 3.72 4.77 5.13 5.43 5.66 5.77 uin [m/s] 42 46 51 56 61 68 uout [m/s] 46 51 56 61 68 74 vin [m/s] 29 54 60 66 72 79 vout [m/s] 54 60 66 72 79 87 win [m/s] 51 30 33 36 40 44 wout [m/s] 96 106 116 128 141 155 ϕin – 0.7 0.64 0.64 0.64 0.64 0.64 ϕout – 0.64 0.64 0.64 0.64 0.64 0.64 ψin – 0.00 0.98 0.98 0.98 0.98 0.98 ψout – −0.98 −0.98 −0.98 −0.98 −0.98 −0.98 hin [m] 468,500 466,000 461,000 456,000 449,000 441,000 hout [m] 466,000 461,000 456,000 449,000 441,000 431,000 Tin [K] 373 371 368 364 360 355 Tout [K] 371 368 364 360 355 349 sin [J/kgK] 1,343 1,347 1,349 1,352 1,355 1,360 pin [Pa] 366,619 366,040 292,470 247,050 201,090 156,400 pout [Pa] 366,040 292,470 247,050 201,090 156,400 115,050 ρin [kg/m] 9.35 8.54 7.42 6.27 5.11 3.99 Int. J. Turbomach. Propuls. Power 2020, 5, 19 12 of 17 ρout [kg/m] 8.54 7.42 6.27 5.11 3.99 2.95 Figure 12. Velocity triangles for the first two rows. Figure 12. Velocity triangles for the first two rows. 3. Results 3. Results CFD Validation CFD Validation The configuration proposed in Figure 12 was simulated via a viscous turbulent CFD analysis, to The configuration proposed in Figure 12 was simulated via a viscous turbulent CFD analysis, to validate the proposed model. The blade shape was designed using Dunham’s parametric method [18] Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 12 of 16 validate the proposed model. The blade shape was designed using Dunham’s parametric method and the number of blades was selected using a modified Zweifel’s criterion [19,20], since ROTs have an [18] and the number of blades was selected using a modified Zweifel’s criterion [19,20], since ROTs axial blade shape. A MATLAB code based on the above parametric method was written to define the define the shape of the blades by interpolating over 1000 points chordwise; then the blade shape was have an axial blade shape. A MATLAB code based on the above parametric method was written to shape of the blades by interpolating over 1000 points chordwise; then the blade shape was imported in imported in SolidWorks and simulated by Ansys Fluent. SolidWorks and simulated by Ansys Fluent. For the CFD simulation a RANS model with a k–ε turbulence model was used. For this For the CFD simulation a RANS model with a k–" turbulence model was used. For this preliminary preliminary design validation, the viscous sublayer was not resolved. Both two- and three- design validation, the viscous sublayer was not resolved. Both two- and three-dimensional simulations dimensional simulations were performed. The relative motion was simulated by the sliding mesh were performed. The relative motion was simulated by the sliding mesh method in the two-dimensional method in the two-dimensional case and by the Multiple Reference Frame (MRF) model in the three- case and by the Multiple Reference Frame (MRF) model in the three-dimensional case. dimensional case. A Conformal Mesh was generated using a growth rate of 1.045 in order to have a mesh continuity A Conformal Mesh was generated using a growth rate of 1.045 in order to have a mesh continuity as gradual as possible with triangular elements, except in the near-wall region where rectangular as gradual as possible with triangular elements, except in the near-wall region where rectangular elements provide better accuracy. A mesh sensitivity analysis was performed (see Figures 13 and 14) elements provide better accuracy. A mesh sensitivity analysis was performed (see Figures 13 and 14) to check the independency of the results from the mesh generated. The objective function for the mesh to check the independency of the results from the mesh generated. The objective function for the sensitivity was the total entropy. mesh sensitivity was the total entropy. Figure 13. Mesh quality before refinement. Figure 13. Mesh quality before refinement. Figure 14. Mesh quality after refinement. The above settings led to an acceptable number of nodes for an accurate two-dimensional study (200,862 nodes and 387,202 elements) and a y+ from 40 to 150 along the blade surfaces. By the same line of reasoning, we obtained 2.6 million hexahedral elements for the 3D study. −5 After a satisfactory convergence was obtained with a residual norm of the order of 10 , different thermodynamic quantities were extracted from the numerical solution and compared with the analytic model, as shown in Figure 15. The output torque was calculated at each time step of the two- and three-dimensional simulations, to compare the power output with the model prediction. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 12 of 16 define the shape of the blades by interpolating over 1000 points chordwise; then the blade shape was imported in SolidWorks and simulated by Ansys Fluent. For the CFD simulation a RANS model with a k–ε turbulence model was used. For this preliminary design validation, the viscous sublayer was not resolved. Both two- and three- dimensional simulations were performed. The relative motion was simulated by the sliding mesh method in the two-dimensional case and by the Multiple Reference Frame (MRF) model in the three- dimensional case. A Conformal Mesh was generated using a growth rate of 1.045 in order to have a mesh continuity as gradual as possible with triangular elements, except in the near-wall region where rectangular elements provide better accuracy. A mesh sensitivity analysis was performed (see Figures 13 and 14) to check the independency of the results from the mesh generated. The objective function for the mesh sensitivity was the total entropy. Int. J. Turbomach. Propuls. Power 2020, 5, 19 13 of 17 Figure 13. Mesh quality before refinement. Figure 14. Mesh quality after refinement. Figure 14. Mesh quality after refinement. The above settings led to an acceptable number of nodes for an accurate two-dimensional study The above settings led to an acceptable number of nodes for an accurate two-dimensional study (200,862 nodes and 387,202 elements) and a y+ from 40 to 150 along the blade surfaces. By the same (200,862 nodes and 387,202 elements) and a y+ from 40 to 150 along the blade surfaces. By the same line of reasoning, we obtained 2.6 million hexahedral elements for the 3D study. line of reasoning, we obtained 2.6 million hexahedral elements for the 3D study. 5 After a satisfactory convergence was obtained with a residual norm of the order of 10 , di erent −5 After a satisfactory convergence was obtained with a residual norm of the order of 10 , different thermodynamic quantities were extracted from the numerical solution and compared with the analytic thermodynamic quantities were extracted from the numerical solution and compared with the model, as shown in Figure 15. The output torque was calculated at each time step of the two- and Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 13 of 16 three-dimensional simulations, to compare the power output with the model prediction. analytic model, as shown in Figure 15. The output torque was calculated at each time step of the two- and three-dimensional simulations, to compare the power output with the model prediction. (a) (b) Figure Figure 15. 15. Contours Contours of of static staticenthalpy enthalpyfor for a si a single ngle tim timeestep: step: (a (a )) 2D mode 2D model;l; ( (bb ))3D 3D model. model. Figure 16 shows the output power for both rows in all the analyzed cases. The continuous Figure 16 shows the output power for both rows in all the analyzed cases. The continuous lines lines represent the model prediction, while the dashed lines represent the 2D configuration and the represent the model prediction, while the dashed lines represent the 2D configuration and the dotted dotted lines are their averages. The asterisks represent the 3D configuration which was analyzed in a lines are their averages. The asterisks represent the 3D configuration which was analyzed in a single single frame. frame. Figure 16. Output power as a function of time. The isentropic efficiency of the whole turbine is about 0.82; however, to make a comparison with the CFD results, only the first two rows are to be considered. Therefore, the isentropic efficiency decreases to 0.759 because of the low efficiency of the first row. For the 2D configuration, the efficiency for the first two rows, calculated on the basis of the average power and enthalpy drop, is about 0.762—for the 3D configuration, about 0.751. The reduction in efficiency from the 2D to the 3D configuration is caused by the increase in the frictional losses, in particular in the hub and tip boundary layers. In fact, since the Soderberg loss model takes into account 3D losses, the slightly higher value of the efficiency in the 2D analysis compared with the “analytical” efficiency can be attributed to the neglection of wall friction (on the disks), which are instead considered in the 3D analysis. Now that the effectiveness of the proposed design algorithm has been demonstrated via the above described CFD simulations, it is possible to formulate an educated assessment of the suitability of the Ljungström turbine for different low-T thermal sources by using different working fluids in different ranges of power: Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 13 of 16 (a) (b) Figure 15. Contours of static enthalpy for a single time step: (a) 2D model; (b) 3D model. Figure 16 shows the output power for both rows in all the analyzed cases. The continuous lines represent the model prediction, while the dashed lines represent the 2D configuration and the dotted lines are their averages. The asterisks represent the 3D configuration which was analyzed in a single Int. J. Turbomach. Propuls. Power 2020, 5, 19 14 of 17 frame. Figure 16. Output power as a function of time. Figure 16. Output power as a function of time. The isentropic eciency of the whole turbine is about 0.82; however, to make a comparison with the CFD results, only the first two rows are to be considered. Therefore, the isentropic eciency The isentropic efficiency of the whole turbine is about 0.82; however, to make a comparison with decreases to 0.759 because of the low eciency of the first row. For the 2D configuration, the eciency the CFD results, only the first two rows are to be considered. Therefore, the isentropic efficiency for the first two rows, calculated on the basis of the average power and enthalpy drop, is about decreases to 0.759 because of the low efficiency of the first row. For the 2D configuration, the 0.762—for the 3D configuration, about 0.751. The reduction in eciency from the 2D to the 3D efficiency for the first two rows, calculated on the basis of the average power and enthalpy drop, is configuration is caused by the increase in the frictional losses, in particular in the hub and tip boundary about 0.762—for the 3D configuration, about 0.751. The reduction in efficiency from the 2D to the 3D layers. In fact, since the Soderberg loss model takes into account 3D losses, the slightly higher value configuration is caused by the increase in the frictional losses, in particular in the hub and tip of the eciency in the 2D analysis compared with the “analytical” eciency can be attributed to the boundary layers. In fact, since the Soderberg loss model takes into account 3D losses, the slightly neglection of wall friction (on the disks), which are instead considered in the 3D analysis. higher value of the efficiency in the 2D analysis compared with the “analytical” efficiency can be Now that the e ectiveness of the proposed design algorithm has been demonstrated via the above attributed to the neglection of wall friction (on the disks), which are instead considered in the 3D described CFD simulations, it is possible to formulate an educated assessment of the suitability of the analysis. Ljungström turbine for di erent low-T thermal sources by using di erent working fluids in di erent Now that the effectiveness of the proposed design algorithm has been demonstrated via the Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 14 of 16 ranges of power: above described CFD simulations, it is possible to formulate an educated assessment of the suitability Figure 17 clearly shows that, as the range of power increases, higher density working fluids are of the Fig Ljung ure 17 clearly ström turbine showfor d s that, as the r ifferent low-T ther ange of power mal sou increases, rces by h uis ghe ingr di den fferent sity wor work king ing fluids fluids ar in e required to match the Ljungström turbine characteristics. req different r uired to m anga es of power: tch the Ljungström turbine characteristics. Figure 17. Di erent power output for di erent working fluids. Figure 17. Different power output for different working fluids. 4. Conclusions To the best of our knowledge, there is not a commercially available Ljungström blade design software, and therefore the present study proposes the use of the kinematic efficiency as a new relevant design parameter in the Ljungström design. However, as discussed in section 2, such a choice opens up a lot of points worthy of further discussion and investigation. The main reason is that this procedure does not take into account the centrifugal forces (absent in axial flows), and therefore the “optimality” of the blade profile proposed by the analytical model must be verified by an experimental validation. In the case studied here, this validation was obtained by a series of 2D time- dependent and by a steady state 3D CFD simulations. Though a full picture of the time evolution of the 3D flow is not available, nevertheless, some useful considerations can be drawn. As clearly shown in Table 9, the axial blade length is a crucial parameter in Ljungström design. By looking at Equation (17), it is clear that once all parameters, except for density, are fixed, in fact the inlet radius of each row is determined by the assumption of a fixed radius ratio, while the flow coefficient is fixed by the assumption of a constant outlet angle. In light of the above, it is also clear that the axial blade length depends on the mass flow rate and on the enthalpy drop. Following Shepherd and Kearton, it was assumed that all blades, except those in the innermost ring, have the same cross-section and they are all set similarly in the rotor. Some authors question the adequacy of these assumptions, used in the past for the design of large scale steam-powered ROT, since they lead to geometrical optimization instead of maximizing the isentropic efficiency [7]. According to these authors, if one removes the assumptions of constant βout and χ, it is possible to exploit two additional degrees of freedom for the axial blade length Equation (17). As a consequence, the formula of the kinematic efficiency changes its structure so that a multi-dimensional optimization is required. Some existing designs split the turbine into a high pressure and a low pressure stage to avoid the reduction in axial blade length in the downstream blade crowns of the turbine. This alternative solution consists of a small reduction in the mass flow rate at the inlet (first stage) and an additional injection in the low pressure wheel; however, this is not in contrast with our model, although further investigations could be useful. Clearly all the proposed changes should be analyzed in more depth because none of the above is supported by calculations. The method we propose is based on kinematic efficiency maximization, while isentropic efficiency determines the turbine performance. Even though ηkin does not take into account fluid dynamic losses, it is maximized with respect to ρr, the ratio between inlet and outlet Int. J. Turbomach. Propuls. Power 2020, 5, 19 15 of 17 4. Conclusions To the best of our knowledge, there is not a commercially available Ljungström blade design software, and therefore the present study proposes the use of the kinematic eciency as a new relevant design parameter in the Ljungström design. However, as discussed in Section 2, such a choice opens up a lot of points worthy of further discussion and investigation. The main reason is that this procedure does not take into account the centrifugal forces (absent in axial flows), and therefore the “optimality” of the blade profile proposed by the analytical model must be verified by an experimental validation. In the case studied here, this validation was obtained by a series of 2D time-dependent and by a steady state 3D CFD simulations. Though a full picture of the time evolution of the 3D flow is not available, nevertheless, some useful considerations can be drawn. As clearly shown in Table 9, the axial blade length is a crucial parameter in Ljungström design. By looking at Equation (17), it is clear that once all parameters, except for density, are fixed, in fact the inlet radius of each row is determined by the assumption of a fixed radius ratio, while the flow coecient is fixed by the assumption of a constant outlet angle. In light of the above, it is also clear that the axial blade length depends on the mass flow rate and on the enthalpy drop. Following Shepherd and Kearton, it was assumed that all blades, except those in the innermost ring, have the same cross-section and they are all set similarly in the rotor. Some authors question the adequacy of these assumptions, used in the past for the design of large scale steam-powered ROT, since they lead to geometrical optimization instead of maximizing the isentropic eciency [7]. According to these authors, if one removes the assumptions of constant and , it is possible to out exploit two additional degrees of freedom for the axial blade length Equation (17). As a consequence, the formula of the kinematic eciency changes its structure so that a multi-dimensional optimization is required. Some existing designs split the turbine into a high pressure and a low pressure stage to avoid the reduction in axial blade length in the downstream blade crowns of the turbine. This alternative solution consists of a small reduction in the mass flow rate at the inlet (first stage) and an additional injection in the low pressure wheel; however, this is not in contrast with our model, although further investigations could be useful. Clearly all the proposed changes should be analyzed in more depth because none of the above is supported by calculations. The method we propose is based on kinematic eciency maximization, while isentropic eciency determines the turbine performance. Even though  does not take into kin account fluid dynamic losses, it is maximized with respect to  , the ratio between inlet and outlet radii. The latter also plays a major role in the axial blade length Equation (17). In the end, the shorter the blade length, the larger the wall friction losses and the lower the isentropic eciency. Author Contributions: Conceptualization, U.C. and E.S.; Formal analysis, U.C.; Investigation, U.C.; Methodology, U.C.; Supervision, E.S.; Validation, E.S.; Writing—original draft, U.C.; Writing—review and editing, E.S. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Conflicts of Interest: The authors declare no conflict of interest. Int. J. Turbomach. Propuls. Power 2020, 5, 19 16 of 17 Nomenclature Symbol Description Unit of Measure Symbol Description Unit of Measure b radial chord [m] V absolute flow velocity [m/s] c specific heat [J/(kgK)] W relative flow velocity [m/s] d specific diameter wf working fluid d shaft diameter [m] wo working oil shaft axial blade length/ [m] h relative inlet angle enthalpy [J/kg] id ideal expansion ratio in inlet  blade encumbrance l blade chord [m] " deviation angle L work [J/kg]  cycle’s eciency cyc L Euler ’s work, [J/kg] [J/kg]  isentropic eciency Eul is l length scale, [m] [m]  kinematic eciency t kin m mass flow rate [kg/s]  turbine eciency turb n specific speed  leakage factor out outlet  Soderberg loss coecient p pressure [Pa]  density [kg/m ] blade to relative outlet velocity P power [kW] ratio r radius [m]  blade solidity Re Reynolds number  allowable stress [MPa] all s entropy [J/(kgK)] ' flow coecient s-s static-to-static  inlet to outlet radii ratio T temperature [T] load coecient U blade peripheral velocity [m/s] ! angular velocity [rpm] References 1. Quoilin, S.; van den Broek, M.; Declaye, S.; Dewallef, P.; Lemort, V. Techno-economic survey of Organic Rankine Cycle (ORC) systems. Renew. Sustain. Energy 2013, 22, 168–186. [CrossRef] 2. Tocci, L.; Pal, T.; Pesmazoglou, I.; Franchetti, B. Small Scale Organic Rankine Cycle (ORC): A Techno-Economic Review. Energies 2017, 10, 413. [CrossRef] 3. Kang, S.H. Design and experimental study of ORC (Organic Rankine Cycle) and radial turbine using R245fa working fluid. Energy 2012, 41, 514–524. [CrossRef] 4. Pei, G.; Li, Y.; Ji, J. Performance evaluation of a micro turbo-expander for application in low-temperature solar electricity generation. J. Zhejiang Univ. 2011, 12, 207–213. [CrossRef] 5. Palumbo, C.F.; Barnabei, V.F.; Preziuso, E.; Coronetta, U. Design and CFD analysis of a Ljungstrom turbine for an ORC cycle in a waste heat recovery application. In Proceedings of the ECOS 2016, Portorož, Slovenia, 19–23 June 2016. 6. Ljungström, F. The Development of the Ljungström Steam Turbine and Air Preheater. Proc. Inst. Mech. Eng. 1949, 160, 211–223. [CrossRef] 7. Casati, E.; Vitale, S.; Pini, M.; Persico, G.; Colonna, P. Centrifugal Turbines for Mini-Organic Rankine Cycle Power Systems. J. Eng. Gas Turbine Power 2014, 136, 122607. [CrossRef] 8. Pini, M.; Persico, G.; Casati, E.; Dossena, V. Preliminary design of a centrifugal turbine for ORC applications. In Proceedings of the First International Seminar on ORC Power Systems, Delft, The Netherlands, 22–23 September 2011. 9. Bell, I.H.; Wronsky, J.; Quoilin, S.; Lemort, V. Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library CoolProp. Ind. Eng. Chem. Res. 2014, 53, 2498–2508. [CrossRef] [PubMed] 10. Franchetti, B.; Pesiridis, A.; Pesmazoglou, I.; Sciubba, E.; Tocci, L. Thermodynamic and technical criteria for the optimal selection of the working fluid in a mini-ORC. In Proceedings of the ECOS 2016, Portorož, Slovenia, 19–23 June 2016. 11. Drescher, U.; Brüggemann, D. Fluid selection for the Organic Rankine Cycle (ORC) in biomass power and heat plants. Appl. Therm. Eng. 2007, 27, 223–228. [CrossRef] Int. J. Turbomach. Propuls. Power 2020, 5, 19 17 of 17 12. Kearton, W.J. Steam Turbine Theory and Practice: A Textbook for Engineering Students; I. Pitman & Sons: London, UK, 1958. 13. Shepherd, D.G. Principles of Turbomachineryi; Macmillan Pub Co.: New York, NY, USA, 1961. 14. Soderberg, C.R. Unpublished Notes; Gas Turbine Laboratory, Massachusetts Institute of Technology: Cambridge, MA, USA, 1949. 15. Ainley, D.G.; Mathieson, G.C.R. An Examination of the Flow and Pressure Losses in Blade Rows of Axial-Flow Turbines; HM Stationery Oce: London, UK, 1951. 16. Craig, H.R.M.; Cox, H.J.A. Performance Estimation of Axial Flow Turbines. Proc. Inst. Mech. Eng. 1971, 185, 407–424. [CrossRef] 17. Sciubba, E. Lezioni di Turbomacchine; Euroma La Goliardica: Rome, Italy, 2001. 18. Dunham, J. A Parametric Method of Turbine Blade Profile Design. In Proceedings of the ASME International Gas Turbine Conference and Products Show, Zurich, Switzerland, 30 March–4 April 1974. 19. Zweifel, O. Optimum Blade Pitch for Turbo-Machines with Special Reference to Blades of Great Curvature. Eng. Dig. 1946, 7, 358–360. 20. Pini, M.; Persico, G.; Casati, E.; Dossena, V. Preliminary design of a centrifugal turbine for ORC applications. J. Eng. Gas Turbines Power 2013, 134, 042312. [CrossRef] © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) license (http://creativecommons.org/licenses/by-nc-nd/4.0/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Turbomachinery, Propulsion and Power Multidisciplinary Digital Publishing Institute

Optimal Design of a Ljungström Turbine for ORC Power Plants: From a 2D model to a 3D CFD Validation

Download PDF
17 pages

Loading next page...
 
/lp/multidisciplinary-digital-publishing-institute/optimal-design-of-a-ljungstr-m-turbine-for-orc-power-plants-from-a-2d-udNBJAHBW4

References (22)

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

Abstract

International Journal of Turbomachinery Propulsion and Power Article Optimal Design of a Ljungström Turbine for ORC Power Plants: From a 2D model to a 3D CFD Validation Umberto Coronetta * and Enrico Sciubba Department of Mechanical and Aerospace Engineering, University of Roma Sapienza, Via Eudossiana 18, 00184 Roma, Italy; enrico.sciubba@uniroma1.it * Correspondence: umberto.coronetta@gmail.com Received: 23 April 2020; Accepted: 14 July 2020; Published: 20 July 2020 Abstract: In the last few years, waste-energy recovery systems based on the Organic Rankine Cycle (ORC) have gained increased attention in the global energy market as a versatile and sustainable technology for thermo-electric energy conversion from low-to-medium temperature sources, up to 350 C. For a long time, water has been the only working fluid commercially adopted in powerplants: axial and, for smaller machines, radial inflow turbines have been the preferred expanders since their gulp capacity matches the -T curve of water steam. The density of most organic compounds displays extremely large variations during the expansion (and the volume flow rate correspondingly increases along the machine channels), so that Radial Outflow Turbines (ROTs) have been recently considered instead of traditional solutions. This work proposes a two-dimensional inviscid model for the stage optimization of a counter-rotating ROT, known as the Ljungström turbine. The study starts by considering five di erent working fluids that satisfy both the gulp requirements of the turbine and the hot source characteristics. On the basis of a limited number of geometric assumptions and for a fixed set of operating conditions, di erent kinematic parameters are optimized to obtain the most ecient cascade configuration. Moreover, as shown in the conclusions, the most ecient blade profile leads to higher friction losses, making further investigation regarding the best configuration necessary. Keywords: turbine CFD; Ljungström turbine; Organic Rankine Cycle 1. Introduction The research and development of the Organic Rankine Cycle (ORC), based energy conversion systems, is mainly aimed at the accurate prediction of the performance of the expander, which strongly influences the eciency of the entire plant. Moreover, based on the characteristics of energy sources exploited by the plant, some expanders are more suitable than others, leading to a higher eciency of the overall process. With the exception of the smallest ORC power plants (few kW), where volumetric expanders of scroll and screw type can be used, whenever higher power outputs are needed, and consequently high expansion ratios and/or high volume flow rates are necessary, dynamic turbines are the only feasible expanders [1]. Tocci et al. [2], in their review of ORC for small scale applications, argue that the selection of the expander type in the range of power 20–70 kW represents an open question: on the one side, the size of volumetric expanders increases exponentially with an increment of output power [2], on the other, traditional radial inflow solutions may lead to unrealistic rotational speeds [3,4]. In light of the above, the purpose of this paper is to investigate the suitability of a special class of Ljungström turbines, specifically designed for the very large expansion rates of organic fluids in the range of power for which a preferred choice is not yet defined and consolidated. To do so, we will start Int. J. Turbomach. Propuls. Power 2020, 5, 19; doi:10.3390/ijtpp5030019 www.mdpi.com/journal/ijtpp Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 2 of 16 In light of the above, the purpose of this paper is to investigate the suitability of a special class of Ljungström turbines, specifically designed for the very large expansion rates of organic fluids in Int. J. Turbomach. Propuls. Power 2020, 5, 19 2 of 17 the range of power for which a preferred choice is not yet defined and consolidated. To do so, we will start from the study of the velocity triangles of a Ljungström type radial outflow turbine that exploit from the s low study -T tof hermal sou the velocity rces and pro triangles of duces s a Ljungström ignificantype t power out radial outflow put [5] and t turbine hen cont that exploits inue tolow-T their opti thermal mizasour tion. ces and produces significant power output [5] and then continue to their optimization. The Ljun The Ljungström gström tturbine urbine o owes wes it its s n name ame tto o it its s in inve vent ntors, ors, the Ljun the Ljungström gström brothers, who brothers, who in in 1908 1908 conceived, d conceived, designed esigned and pr and pr oduced oduced a stat a statorless orless steam t steam urbine, turbine, as sho as shown wn in Fig in u Figur re 1, cons e 1, consisting isting of twof o count two counter er-rota-r ting co otating axi coaxial al disks disks fitted w fitteditwith h seve several ral con concentric centric bla blade de crowns crowns [6[]. 6 ].The in The inlet let sect section ion is is located located rad radially ially ne near ar the sh the shaft aft while the while the outlet is at outlet is at th the e outer r outer radius adius at the exit at the exit of the last blade of the last bladecrown. crown. This arrangement compensates for the increase in the volume of the fluid during the expansion by a This arrangement compensates for the increase in the volume of the fluid during the expansion by a correspondin corresponding g incr increase ease in t in the he flow flow are area a along along the flo the flow w path. path. Figure 1. Ljungström Radially Outward Turbine (ROT): Fixed and Counter-rotating arrangements [7,8]. Figure 1. Ljungström Radially Outward Turbine (ROT): Fixed and Counter-rotating arrangements [7,8]. 2. Materials and Methods 2. Materials and Methods 2.1. Fluid Selection In this work, a hydro-treated mineral oil, Paratherm NF, is the heat transfer medium taken as the 2.1. Fluid Selection hot source. It enters the boiler at 140 C (413 K) with a mass flow rate of 4 kg/s. The cold source is water, In this work, a hydro-treated mineral oil, Paratherm NF, is the heat transfer medium taken as which enters the condenser at room temperature (25 C or 298 K) and is assumed to be available in the the hot source. It enters the boiler at 140 °C (413 K) with a mass flow rate of 4 kg/s. The cold source is necessary quantity. The selection of the working fluid and the identification of the thermodynamic water, which enters the condenser at room temperature (25 °C or 298 K) and is assumed to be cycle was performed via a MATLAB code by using the CoolProp library [9] to calculate the fluid available in the necessary quantity. The selection of the working fluid and the identification of the properties. A first selection of fluids from the CoolProp library was performed under two criteria: first, thermodynamic cycle was performed via a MATLAB code by using the CoolProp library [9] to since the source is at low temperature, the fluid should have a low critical pressure [10,11]; second, calculate the fluid properties. A first selection of fluids from the CoolProp library was performed the Ljungström turbine requires relatively high volume flow rates at its inlet and a well-identified under two criteria: first, since the source is at low temperature, the fluid should have a low critical variation in the density along the radius due to its peculiar cross-sectional area, and thus the fluid pressure [10,11]; second, the Ljungström turbine requires relatively high volume flow rates at its inlet ought to have a low density at turbine inlet (i.e., high values of the coecients in the Redlich–Kwong and a well-identified variation in the density along the radius due to its peculiar cross-sectional area, equation of state). According to these considerations the following fluids were considered: an alicyclic and thus the fluid ought to have a low density at turbine inlet (i.e., high values of the coefficients in hydrocarbon, Cyclopentane (C5H10); an azeotropic mixture, SES36; a haloalkane refrigerant, R134a the Redlich–Kwong equation of state). According to these considerations the following fluids were (CH2FCF3); an alkane refrigerant, R601 (C5H12); a hydrofluorocarbon refrigerant, R245fa (C3H3F5). considered: an alicyclic hydrocarbon, Cyclopentane (C5H10); an azeotropic mixture, SES36; a For each working fluid, as shown in Figure 2, a specific thermodynamic cycle was designed by haloalkane refrigerant, R134a (CH2FCF3); an alkane refrigerant, R601 (C5H12); a hydrofluorocarbon fixing the top and bottom pinch points and specifying a 10 C water temperature rise in the condenser. refrigerant, R245fa (C3H3F5). For a correct comparison, the thermal power input is the same for all cycles considered. The heat input For each working fluid, as shown in Figure 2, a specific thermodynamic cycle was designed by and output in the boiler and the condenser were then matched by iteratively adjusting the working fixing the top and bottom pinch points and specifying a 10 °C water temperature rise in the fluid mass flow rate or/and the evaporation temperature. Therefore, the turbine enthalpy drop and the condenser. For a correct comparison, the thermal power input is the same for all cycles considered. pump work were determined with an assumed 0.8 eciency for the pump and 0.83 (polytropic) for the The heat input and output in the boiler and the condenser were then matched by iteratively adjusting turbine. The overall cycle eciency is equal to: the working fluid mass flow rate or/and the evaporation temperature. Therefore, the turbine enthalpy drop and the pump work were determined with m  an L as sumed L 0.8 efficiency for the pump and 0.83 w f turb pump = , (1) cyc . (polytropic) for the turbine. The overall cycle efficiency is equal to: m  DH oil oil Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 3 of 16 Int. J. Turbomach. Propuls. Power 2020, 5, 19 3 of 17 = , (1) Figure 2. Temperature–entropy diagram for the fluids considered in this study fluids, in the same scale Figure 2. Temperature–entropy diagram for the fluids considered in this study fluids, in the same as the usual water Mollier diagram. scale as the usual water Mollier diagram. Tables 1–5 present the thermodynamic properties of the considered fluids at the relevant stations Tables 1–5 present the thermodynamic properties of the considered fluids at the relevant stations in the process. Figures 3–7 clearly show that, in all of the cycles, a regenerator recovers a portion of the in the process. Figures 3–7 clearly show that, in all of the cycles, a regenerator recovers a portion of sensible heat content of the working fluid before it enters the condenser. The percentage of thermal the sensible heat content of the working fluid before it enters the condenser. The percentage of power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Cyclopentane, thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the SES36. Eciencies Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the and mass flow rates (indicative of expander size) are reported in Table 6. SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6. Table 1. Cyclopentane thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 73,593 378,072 2.04 Pump (2) 313.00 73,593 18,013 725.38 Economizer (3) 313.06 250,885 17,768 725.52 Evaporator (4) 353.00 250,885 61,970 682.63 Superheater (5) 353.00 250,885 426,794 6.43 Turbine (6) 373.00 250,885 458,217 6.01 Regenerator (7) 341.04 73,593 416,033 1.86 at the inlet of each equipment. Table 2. R601 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Figure 3. Cyclopentane T–S diagram. Equipment [K] [Pa] [J/kg] [kg/m ] Table 1. Cyclopentane thermodynamic parameters. Condenser (1) 313.00 115,093 363,518 3.35 Pump (2) 313.00 115,093 9010 605.85 1 1 1 1 Temperature Pressure Enthalpy Density Economizer (3) 313.09 366,619 9425 606.13 Equipment [K] [Pa] [J/kg] [kg/m ] Evaporator (4) 353.00 366,619 108,917 562.19 Superheater (5) 353.00 366,619 427,177 10.11 Condenser (1) 313.00 73,593 378,072 2.04 Turbine (6) 373.00 366,619 468,490 9.35 Pump (2) 313.00 73,593 −18,013 725.38 Regenerator (7) 348.94 115,093 430,071 2.96 Economizer (3) 313.06 250,885 −17,768 725.52 at the inlet of each equipment. Evaporator (4) 353.00 250,885 61,970 682.63 Superheater (5) 353.00 250,885 426,794 6.43 Turbine (6) 373.00 250,885 458,217 6.01 Regenerator (7) 341.04 73,593 416,033 1.86 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2020, 5, 19 4 of 17 Table 3. R134a thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 1012,509 419,363 49.87 Pump (2) 313.00 1012,509 256,185 1147.37 Economizer (3) 314.11 2624,797 257,590 1155.48 Evaporator (4) 353.00 2624,797 322,105 929.39 Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 3 of 16 Superheater (5) 353.00 2624,797 428,830 154.37 Turbine (6) 373.00 2624,797 460,255 121.11 Regenerator (7) 333.93 1012,509 442,132 43.82 = , (1) at the inlet of each equipment. Table 4. R245fa thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 249,412 435,245 13.95 Pump (2) 313.00 249,412 252,838 1297.13 Economizer (3) 313.55 1648,462 253,916 1300.84 Evaporator (4) 385.44 1648,462 360,253 1038.29 Superheater (5) 385.44 1648,462 482,315 98.34 Turbine (6) 405.44 1648,462 509,011 84.42 Regenerator (7) 354.94 249,412 476,033 11.89 at the inlet of each equipment. Table 5. SES36 thermodynamic parameters. Figure 2. Temperature–entropy diagram for the fluid 1 s 1considered in this stu 1 dy fluids, in the 1 same Temperature Pressure Enthalpy Density Equipment scale as the usual water Mollier diagram. [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 116,523 389,772 8.78 Tables 1–5 present the thermodynamic properties of the considered fluids at the relevant stations Pump (2) 313.00 116,523 233,461 1333.94 in the process. Figures 3–7 clearly show that, in all of the cycles, a regenerator recovers a portion of Economizer (3) 313.31 857,341 234,016 1335.99 Evaporator (4) 386.78 857,341 327,876 1104.87 the sensible heat content of the working fluid before it enters the condenser. The percentage of Superheater (5) 386.78 857,341 442,402 63.85 thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Turbine (6) 406.78 857,341 466,484 56.91 Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the Regenerator (7) 370.35 116,523 439,563 7.21 SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6. at the inlet of each equipment. Figure 3. Cyclopentane T–S diagram. Figure 3. Cyclopentane T–S diagram. Table 1. Cyclopentane thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 73,593 378,072 2.04 Pump (2) 313.00 73,593 −18,013 725.38 Economizer (3) 313.06 250,885 −17,768 725.52 Evaporator (4) 353.00 250,885 61,970 682.63 Superheater (5) 353.00 250,885 426,794 6.43 Turbine (6) 373.00 250,885 458,217 6.01 Regenerator (7) 341.04 73,593 416,033 1.86 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 4 of 16 Int. J. Turbomach. Propuls. Power 2020, 5, 19 5 of 17 Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 4 of 16 Figure 4. R601 T–S diagram. Table 2. R601 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 115,093 363,518 3.35 Pump (2) 313.00 115,093 9,010 605.85 Economizer (3) 313.09 366,619 9,425 606.13 Evaporator (4) 353.00 366,619 108,917 562.19 Superheater (5) 353.00 366,619 427,177 10.11 Turbine (6) 373.00 366,619 468,490 9.35 Regenerator (7) 348.94 115,093 430,071 2.96 at the inlet of each equipment. Figure 4. R601 T–S diagram. Figure 4. R601 T–S diagram. Table 2. R601 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 115,093 363,518 3.35 Pump (2) 313.00 115,093 9,010 605.85 Economizer (3) 313.09 366,619 9,425 606.13 Evaporator (4) 353.00 366,619 108,917 562.19 Superheater (5) 353.00 366,619 427,177 10.11 Turbine (6) 373.00 366,619 468,490 9.35 Regenerator (7) 348.94 115,093 430,071 2.96 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 5 of 16 Figure 5. R134a T–S diagram. Figure 5. R134a T–S diagram. Table 3. R134a thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 1012,509 419,363 49.87 Pump (2) 313.00 1012,509 256,185 1147.37 Economizer (3) 314.11 2624,797 257,590 1155.48 Evaporator (4) 353.00 2624,797 322,105 929.39 Superheater (5) 353.00 2624,797 428,830 154.37 Turbine (6) 373.00 2624,797 460,255 121.11 Regenerator (7) 333.93 1012,509 442,132 43.82 at the inlet of each equipment. Figure 6. R245fa T–S diagram. Figure 5. R134a T–S diagram. Figure 6. R245fa T–S diagram. Table 3. R134a thermodynamic parameters. Table 4. R245fa thermodynamic parameters. 1 1 1 1 1 1 1 1 Temperature Pressure Enthalpy Density Temperature Pressure Enthalpy Density Equipment Equipment [K] [Pa] [J/kg] [kg/m ] [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 1012,509 419,363 49.87 Condenser (1) 313.00 249,412 435,245 13.95 Pump (2) 313.00 1012,509 256,185 1147.37 Pump (2) 313.00 249,412 252,838 1297.13 Economizer (3) 314.11 2624,797 257,590 1155.48 Economizer (3) 313.55 1648,462 253,916 1300.84 Evaporator (4) 353.00 2624,797 322,105 929.39 Evaporator (4) 385.44 1648,462 360,253 1038.29 Superheater (5) 353.00 2624,797 428,830 154.37 Superheater (5) 385.44 1648,462 482,315 98.34 Turbine (6) 373.00 2624,797 460,255 121.11 Turbine (6) 405.44 1648,462 509,011 84.42 Regenerator (7) 333.93 1012,509 442,132 43.82 Regenerator (7) 354.94 249,412 476,033 11.89 at the inlet of each equipment. at the inlet of each equipment. Figure 7. SES36 T–S diagram. Table 5. SES36 thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 116,523 389,772 8.78 Pump (2) 313.00 116,523 233,461 1333.94 Economizer (3) 313.31 857,341 234,016 1335.99 Evaporator (4) 386.78 857,341 327,876 1104.87 Superheater (5) 386.78 857,341 442,402 63.85 Turbine (6) 406.78 857,341 466,484 56.91 Regenerator (7) 370.35 116,523 439,563 7.21 at the inlet of each equipment. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 5 of 16 Figure 6. R245fa T–S diagram. Table 4. R245fa thermodynamic parameters. 1 1 1 1 Temperature Pressure Enthalpy Density Equipment [K] [Pa] [J/kg] [kg/m ] Condenser (1) 313.00 249,412 435,245 13.95 Pump (2) 313.00 249,412 252,838 1297.13 Economizer (3) 313.55 1648,462 253,916 1300.84 Evaporator (4) 385.44 1648,462 360,253 1038.29 Superheater (5) 385.44 1648,462 482,315 98.34 Turbine (6) 405.44 1648,462 509,011 84.42 Regenerator (7) 354.94 249,412 476,033 11.89 Int. J. Turbomach. Propuls. Power 2020, 5, 19 6 of 17 at the inlet of each equipment. Figure 7. SES36 T–S diagram. Figure 7. SES36 T–S diagram. Table 6. Eciency, mass flow rate and expansion ratio for the considered fluids. Table 5. SES36 thermodynamic parameters. m ˙ 1 1 1 1 Fluid Temperature Pressure Enthalpy Density [–] [kg/s] [–] Equipment [K] [Pa] [J/kg] [kg/m ] Cyclopentane 9.63% 1.35 3.41 Condenser (1) 313.00 116,523 389,772 8.78 R601 9.79% 1.50 3.19 Pump (2) 313.00 116,523 233,461 1333.94 R134a 11.23% 4.03 2.59 R245fa 16.34% 1.50 6.61 Economizer (3) 313.31 857,341 234,016 1335.99 SES36 15.46% 1.50 7.36 Evaporator (4) 386.78 857,341 327,876 1104.87 Superheater (5) 386.78 857,341 442,402 63.85 2.2. The Optimization Algorithm Turbine (6) 406.78 857,341 466,484 56.91 Regenerator (7) 370.35 116,523 439,563 7.21 The standard formulation of the Euler work (2) shows that, for a Ljungström “stage”, the only at the inlet of each equipment. positive contribution to the specific work is given by the increase in the relative velocity, because u is out only slightly higher than u : in 2 2 2 2 2 2 v v + u u + w w in out in out out in L = , (2) Eul The novel method for the mean-line design proposed here begins with the following two positions, first suggested by Kearton [12] and Shepherd [13]: All blades, except those in the innermost ring, have the same cross-section, the same profile and the same stagger angle and, as a consequence, also have the same outlet angle: = const., (3) out In all blade crowns, the ratio of the relative outlet steam velocity to the peripheral velocity at the outlet edge of the blade ring is constant and equal to: out = < 1, (4) out In all blade crowns, the ratio of the radial chord of the blade to the outlet radius is constant and equal to: b r in = 1 = < 1, (5) r r out out Due to the axial symmetry of the turbine, the absolute velocity at the inlet of the first row can be considered completely radial. Int. J. Turbomach. Propuls. Power 2020, 5, 19 7 of 17 The above assumptions suce for a definition of kinematic eciency. Since the relative velocity provides the only positive contribution to the Euler ’s work, the absolute velocity and the peripheral velocity are to be considered as losses, and this leads to the following formula: 2 2 2 2 v v + u u out in out in = 1 , (6) kin 2 2 w w out in That can be rewritten as: 2 2 1   1 + 2 2 cos r r out = 1 , (7) kin,i 2 2 1  (1 + 4 4 cos ) r r out This equation provides additional useful information because it is similar for every i-th blade ring except for the first one, where it takes the following form: 2 + 1 2 cos r r out cos = 1 , (8) kin,1 cos The parameters to be considered to maximize the kinematic eciency are, therefore,  , , r out Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 7 of 16 and . Due to the second and third assumptions it is not possible to maximize both  and 0 kin,i kin ,1 with respect to  and , thus only the  equation is considered here for the optimization. r kin,i It can be easily verified that the derivative of the kinematic efficiency (7) is linearly dependent It can be easily verified that the derivative of the kinematic eciency (7) is linearly dependent on on both χ and βout and thus it is not possible to maximize this equation with respect to these both  and and thus it is not possible to maximize this equation with respect to these parameters. out parameters. This conclusion is reinforced by a closer inspection of Equation (7), which is maximized This conclusion is reinforced by a closer inspection of Equation (7), which is maximized either when either when the numerator of the second term tends to zero or the denominator to a maximum. Both the numerator of the second term tends to zero or the denominator to a maximum. Both cases leading cases leading to ideal turbine configurations: to ideal turbine configurations: • χ = 1, meets the maximum of the efficiency for ρr = cos(βout)/2, and causes the numerator to be = 1, meets the maximum of the eciency for  = cos( )/2, and causes the numerator to be r out null. In fact, in this case uin = uout and vin = vout so that the velocity triangle is similar to that of an null. In fact, in this case u = u and v = v so that the velocity triangle is similar to that of an out out in in impulse turbine; impulse turbine; • βout = 0 causes the denominator to be at its maximum, in which case win = 0 because the triangle = 0 causes the denominator to be at its maximum, in which case w = 0 because the triangle out collapses onto a line. in collapses onto a line. Therefore, the only parameter to be considered in the velocity triangles optimization is ρr, which allows us to Thereforma e, th ximiz e only e th parameter e kinematic to effi be consider ciency for diff ed in erent the velocity values of triangles βout and optimization χ, as shown in Fig is  , which ure 8. allows us to maximize the kinematic eciency for di erent values of and , as shown in Figure 8. out (a) (b) Figure 8. Kinematic eciency: (a) Variation of the kinematic eciency with  for di erent and Figure 8. Kinematic efficiency: (a) Variation of the kinematic efficiency with rρr for different outβout and fixed ; (b) Variation of the kinematic eciency with  for di erent  and fixed . fixed χ; (b) Variation of the kinematic efficiency with rρr for different χ and fixed outβout. Considering the above, we can now proceed to the study of the isentropic efficiency, starting from its equation and the loss model adopted here: (9) ℎ = , In order to obtain the value of wout,id, Soderberg [14] loss correlation was used: = 1 , (10) And after some rearrangement, the formula of the isentropic efficiency for the i-th row, with the exception of the first one, is obtained: = , (11) , , While for the first row: = , (12) The Soderberg loss coefficient is adopted for all the ensuing calculations, but different models were tested as well, such as Ainley and Mathieson [15] and Craig and Cox [16] prediction models. With little modifications in the code, the loss coefficient ξr can be evaluated for any of the Int. J. Turbomach. Propuls. Power 2020, 5, 19 8 of 17 Considering the above, we can now proceed to the study of the isentropic eciency, starting from its equation and the loss model adopted here: 2 2 2 2 w w u u out,id in out in Dh = , (9) ss In order to obtain the value of w , Soderberg [14] loss correlation was used: out,id w = w 1 +  , (10) out,id out r And after some rearrangement, the formula of the isentropic eciency for the i-th row, with the exception of the first one, is obtained: cos r r 2 out 2 1 + = ! , (11) is,i p 1+ +4 4 1+ cos 1+ r,i1 r r r,i1 out r,i 1 +  1 2 2 r r While for the first row: cos r out r =   , (12) is,1 1+ r,1 1 +  1 2 2 cos r 0 The Soderberg loss coecient is adopted for all the ensuing calculations, but di erent models were Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 8 of 16 tested as well, such as Ainley and Mathieson [15] and Craig and Cox [16] prediction models. With little modifications in the code, the loss coecient  can be evaluated for any of the abovementioned abovementioned prediction models. It was found that the adoption of different loss models does not prediction models. It was found that the adoption of di erent loss models does not noticeably influence noticeably influence the results, as shown in Figure 9. the results, as shown in Figure 9. (a) (b) Figure Figure 9. 9. Isentr Isentropic eff opic eiciency: ciency: ( (a a)) Variation Variation of of the i the isentr sentropi opicce eff ciency iciency wit withh ρr for for di differe erent nt βout and and r out fixed ; (b) Variation of the kinematic eciency with  for di erent  and fixed . fixed χ; (b) Variation of the kinematic efficiency with ρr for different χ and fixed βout. r out Since  depends on the fluid dynamic and geometric conditions that are di erent in each row, Since ξr depends on the fluid dynamic and geometric conditions that are different in each row, it is not possible to choose a single value of  for all rows that optimizes  . This justifies the choice r is,i it is not possible to choose a single value of ρr for all rows that optimizes ηis,i. This justifies the choice of a velocity triangles optimization based on the selection of the  that maximizes the  . of a velocity triangles optimization based on the selection of the ρr that maximizes the ηkin,i kin,i. As shown in Figure 10, for a fixed geometry (fixed  and ) and by varying the velocity, out As shown in Figure 10, for a fixed geometry (fixed χ and βout) and by varying the velocity, and and consequently  , the maximum of the isentropic eciency is close to the range of the maxima consequently ξr, the maximum of the isentropic efficiency is close to the range of the maxima of the kinematic efficiency. This supports the choice of an optimization procedure based on the velocity triangles. Figure 10. Comparing the variation of ηkin and ηis with ρr for fixed χ and βout and different flow velocity. A MATLAB code was written to obtain a complete mean line analysis, and then the flow in the resulting geometry was simulated via a CFD analysis to validate the model. The inputs of the algorithm are shown in Table 7. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 8 of 16 abovementioned prediction models. It was found that the adoption of different loss models does not noticeably influence the results, as shown in Figure 9. (a) (b) Figure 9. Isentropic efficiency: (a) Variation of the isentropic efficiency with ρr for different βout and fixed χ; (b) Variation of the kinematic efficiency with ρr for different χ and fixed βout. Since ξr depends on the fluid dynamic and geometric conditions that are different in each row, it is not possible to choose a single value of ρr for all rows that optimizes ηis,i. This justifies the choice of a velocity triangles optimization based on the selection of the ρr that maximizes the ηkin,i. Int. J. Turbomach. Propuls. Power 2020, 5, 19 9 of 17 As shown in Figure 10, for a fixed geometry (fixed χ and βout) and by varying the velocity, and consequently ξr, the maximum of the isentropic efficiency is close to the range of the maxima of the kinematic efficiency. This supports the choice of an optimization procedure based on the velocity of the kinematic eciency. This supports the choice of an optimization procedure based on the triangles. velocity triangles. Figure 10. Comparing the variation of  and  with  for fixed  and and di erent flow velocity. r out kin is Figure 10. Comparing the variation of ηkin and ηis with ρr for fixed χ and βout and different flow velocity. A MATLAB code was written to obtain a complete mean line analysis, and then the flow in the A MATLAB code was written to obtain a complete mean line analysis, and then the flow in the resulting geometry was simulated via a CFD analysis to validate the model. The inputs of the algorithm resulting geometry was simulated via a CFD analysis to validate the model. The inputs of the are shown in Table 7. algorithm are shown in Table 7. Table 7. Inputs and constraints of the algorithm. m [kg/s] Working Fluid Mass Flow Rate  [MPa] Allowable Stress at the Shaft 80 all p [Pa] Inlet pressure ! [rpm] Maximum speed 4400 max in p [Pa] Outlet pressure  [–] Blade encumbrance 0.85 out T [K] Inlet temperature l [m] Minimum blade height 0.01 in min b [m] Radial chord 0.01 min [deg] Minimum relative outlet angle 18 out By setting a trial rotational speed ! = ! the shaft diameter can be found that sets a lower max limit to the inlet radius at the entrance of the first row. Thence we obtain a matrix of all possible combinations of r and that satisfy the continuity equation for all the l > l . in in min The abovementioned matrix leads to another matrix of all possible values of  . By taking the kin maximum value we can also find the optimal values of  and . With all these parameters in place, it is now possible to calculate the loss coecient  and the enthalpy drop for the first stage (Dh ). n rows Dh = Dh , (13) i tot i = 1 If (13) is not satisfied it is possible to design additional rows and if an exact match is not achieved, a modified rotational speed ! = ! D! is selected and the process is iteratively repeated, as shown max in Figure 11. Since the maximum eciency is reached for b = 0 and = 0 and since these conditions are out clearly unphysical, a minimum value must be arbitrarily specified. In addition, the change of the cross-sectional area along the radius requires the designer to specify a minimum value of the axial blade length at the inlet: w f h = , (14) min 2!r  tan in b where the leakage factor  was assumed here to be 0.98. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 9 of 16 Table 7. Inputs and constraints of the algorithm. [kg/s] Working Fluid Mass Flow Rate τall [MPa] Allowable Stress at the Shaft 80 pin [Pa] Inlet pressure ωmax [rpm] Maximum speed 4400 pout [Pa] Outlet pressure δb [–] Blade encumbrance 0.85 Tin [K] Inlet temperature lmin [m] Minimum blade height 0.01 bmin [m] Radial chord 0.01 βout [deg] Minimum relative outlet angle 18° By setting a trial rotational speed ω = ωmax the shaft diameter can be found that sets a lower limit to the inlet radius at the entrance of the first row. Thence we obtain a matrix of all possible combinations of rin and βin that satisfy the continuity equation for all the l > lmin. The abovementioned matrix leads to another matrix of all possible values of ηkin. By taking the maximum value we can also find the optimal values of ρr and χ. With all these parameters in place, it is now possible to calculate the loss coefficient ξr and the enthalpy drop for the first stage (∆hi). ∑ ∆ℎ =∆ℎ , (13) If (13) is not satisfied it is possible to design additional rows and if an exact match is not achieved, a modified rotational speed ω = ωmax − ∆ω is selected and the process is iteratively repeated, as shown Int. J. Turbomach. Propuls. Power 2020, 5, 19 10 of 17 in Figure 11. Figure 11. Schematic representation of the design algorithm. Figure 11. Schematic representation of the design algorithm. The above equation has two degrees of freedom. By fixing a proper range of values for and r 0 in Since the maximum efficiency is reached for b = 0 and βout = 0 and since these conditions are all possible combinations that meet the equation can be calculated. These ranges are constrained by clearly unphysical, a minimum value must be arbitrarily specified. In addition, the change of the two di erent conditions. must be limited in order to obtain w > w : out, 0 1 in ,1 cross-sectional area along the radius requires the designer to specify a minimum value of the axial blade length at the inlet: 0 1 B C B 1 C B C B C 0 < < cos  ∙ , (15) 0 B r C @ 2 2 2 A ℎ = , (14) 1 +  2 cos + r r out r where the leakage factor ξ was assumed here to be 0.98. while r for the first row must be larger than the shaft radius that was calculated with the in following formula: 3 16P d  K , (16) sha f t where K = 1.3 is a safety factor, P [MPa] the power at the shaft and ! [rpm] its rotational speed [17]. Subsequently the code designs the first row using the inputs provided by the optimization function. It starts from the velocity triangles, calculates the isentropic heat drop using the loss coecients and finally derives the thermodynamic quantities at the outlet of the row. The quantities thus derived are used as inputs for the following stage, leading to an iterative process. The proposed algorithm leads to a choice in the number of stages as a function of the maximum angular speed. w f h = , (17) min 2!r  tan . in b As shown in Table 8 and, as expected, the  parameter is strictly correlated to the eciency of the turbine. A closer look to these results reveals that the less ecient thermodynamic cycles presents the higher isentropic eciencies. The Ljungström turbine performs better with low enthalpy drop and high-volume flow rates (higher specific speed n ). Moreover, since the cross-sectional area depends on the volume flow rate variation, not all the selected fluids lead to geometrical characteristics, such as Int. J. Turbomach. Propuls. Power 2020, 5, 19 11 of 17 axial blade length, which correspond to the allowed manufacturing limits. As shown in Table 9, the only suitable fluids are Cyclopentane and R601, the only ones leading to an axial blade length acceptable for the imposed constraints. Table 8. Crucial design parameters for di erent working fluids. ! Power Fluid n of Rows n kin is [rpm] [kW] Cyclopentane 6 0.587 3648 56 0.914 0.836 0.925 R601 6 0.586 4016 56 0.897 0.826 0.909 R134a 6 0.434 3714 61 0.815 0.729 0.821 R245fa 4 0.499 9171 34 0.698 0.590 0.688 SES36 6 0.350 4541 32 0.783 0.675 0.792 calculated at first stage. Table 9. Axial blade length (expressed in (m)) at the inlet of each row for di erent working fluids. Fluid 1st Row 2nd Row 3rd Row 4th Row 5th Row 6th Row Cyclopentane 0.010 0.011 0.011 0.011 0.011 0.012 R601 0.010 0.010 0.010 0.009 0.010 0.010 R134a 0.011 0.008 0.006 0.005 0.004 0.003 R245fa 0.010 0.006 0.003 0.003 SES36 0.010 0.007 0.005 0.004 0.004 0.004 In the present study, R601 was chosen as working fluid because, despite the slightly lower eciency of its cycle, it is much less flammable than Cyclopentane. Moreover, its condenser pressure is about 1 bar against 0.7 bars for Cyclopentane, making the manufacturing easier. The design parameters thus obtained, as shown in Table 10, make it possible to draw the velocity triangles of the best performing configuration, as shown below in Figure 12: Table 10. Output parameters for the selected fluid. Parameter Unit 1st Row 2nd Row 3rd Row 4th Row 5th Row 6th Row – 0.615 0.850 0.845 0.839 0.834 0.830 is n – 0.54 0.31 0.29 0.27 0.26 0.26 d – 3.72 4.77 5.13 5.43 5.66 5.77 u [m/s] 42 46 51 56 61 68 in u [m/s] 46 51 56 61 68 74 out v [m/s] 29 54 60 66 72 79 in v [m/s] 54 60 66 72 79 87 out w [m/s] 51 30 33 36 40 44 in w [m/s] 96 106 116 128 141 155 out – 0.7 0.64 0.64 0.64 0.64 0.64 in – 0.64 0.64 0.64 0.64 0.64 0.64 out – 0.00 0.98 0.98 0.98 0.98 0.98 in – 0.98 0.98 0.98 0.98 0.98 0.98 out h [m] 468,500 466,000 461,000 456,000 449,000 441,000 in h [m] 466,000 461,000 456,000 449,000 441,000 431,000 out T [K] 373 371 368 364 360 355 in T [K] 371 368 364 360 355 349 out s [J/kgK] 1,343 1,347 1,349 1,352 1,355 1,360 in p [Pa] 366,619 366,040 292,470 247,050 201,090 156,400 in p [Pa] 366,040 292,470 247,050 201,090 156,400 115,050 out [kg/m ] 9.35 8.54 7.42 6.27 5.11 3.99 in [kg/m ] 8.54 7.42 6.27 5.11 3.99 2.95 out Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 11 of 16 The design parameters thus obtained, as shown in Table 10, make it possible to draw the velocity triangles of the best performing configuration, as shown below in Figure 12: Table 10. Output parameters for the selected fluid. Parameter Unit 1st Row 2nd Row 3rd Row 4th Row 5th Row 6th Row ηis – 0.615 0.850 0.845 0.839 0.834 0.830 ns – 0.54 0.31 0.29 0.27 0.26 0.26 ds – 3.72 4.77 5.13 5.43 5.66 5.77 uin [m/s] 42 46 51 56 61 68 uout [m/s] 46 51 56 61 68 74 vin [m/s] 29 54 60 66 72 79 vout [m/s] 54 60 66 72 79 87 win [m/s] 51 30 33 36 40 44 wout [m/s] 96 106 116 128 141 155 ϕin – 0.7 0.64 0.64 0.64 0.64 0.64 ϕout – 0.64 0.64 0.64 0.64 0.64 0.64 ψin – 0.00 0.98 0.98 0.98 0.98 0.98 ψout – −0.98 −0.98 −0.98 −0.98 −0.98 −0.98 hin [m] 468,500 466,000 461,000 456,000 449,000 441,000 hout [m] 466,000 461,000 456,000 449,000 441,000 431,000 Tin [K] 373 371 368 364 360 355 Tout [K] 371 368 364 360 355 349 sin [J/kgK] 1,343 1,347 1,349 1,352 1,355 1,360 pin [Pa] 366,619 366,040 292,470 247,050 201,090 156,400 pout [Pa] 366,040 292,470 247,050 201,090 156,400 115,050 ρin [kg/m] 9.35 8.54 7.42 6.27 5.11 3.99 Int. J. Turbomach. Propuls. Power 2020, 5, 19 12 of 17 ρout [kg/m] 8.54 7.42 6.27 5.11 3.99 2.95 Figure 12. Velocity triangles for the first two rows. Figure 12. Velocity triangles for the first two rows. 3. Results 3. Results CFD Validation CFD Validation The configuration proposed in Figure 12 was simulated via a viscous turbulent CFD analysis, to The configuration proposed in Figure 12 was simulated via a viscous turbulent CFD analysis, to validate the proposed model. The blade shape was designed using Dunham’s parametric method [18] Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 12 of 16 validate the proposed model. The blade shape was designed using Dunham’s parametric method and the number of blades was selected using a modified Zweifel’s criterion [19,20], since ROTs have an [18] and the number of blades was selected using a modified Zweifel’s criterion [19,20], since ROTs axial blade shape. A MATLAB code based on the above parametric method was written to define the define the shape of the blades by interpolating over 1000 points chordwise; then the blade shape was have an axial blade shape. A MATLAB code based on the above parametric method was written to shape of the blades by interpolating over 1000 points chordwise; then the blade shape was imported in imported in SolidWorks and simulated by Ansys Fluent. SolidWorks and simulated by Ansys Fluent. For the CFD simulation a RANS model with a k–ε turbulence model was used. For this For the CFD simulation a RANS model with a k–" turbulence model was used. For this preliminary preliminary design validation, the viscous sublayer was not resolved. Both two- and three- design validation, the viscous sublayer was not resolved. Both two- and three-dimensional simulations dimensional simulations were performed. The relative motion was simulated by the sliding mesh were performed. The relative motion was simulated by the sliding mesh method in the two-dimensional method in the two-dimensional case and by the Multiple Reference Frame (MRF) model in the three- case and by the Multiple Reference Frame (MRF) model in the three-dimensional case. dimensional case. A Conformal Mesh was generated using a growth rate of 1.045 in order to have a mesh continuity A Conformal Mesh was generated using a growth rate of 1.045 in order to have a mesh continuity as gradual as possible with triangular elements, except in the near-wall region where rectangular as gradual as possible with triangular elements, except in the near-wall region where rectangular elements provide better accuracy. A mesh sensitivity analysis was performed (see Figures 13 and 14) elements provide better accuracy. A mesh sensitivity analysis was performed (see Figures 13 and 14) to check the independency of the results from the mesh generated. The objective function for the mesh to check the independency of the results from the mesh generated. The objective function for the sensitivity was the total entropy. mesh sensitivity was the total entropy. Figure 13. Mesh quality before refinement. Figure 13. Mesh quality before refinement. Figure 14. Mesh quality after refinement. The above settings led to an acceptable number of nodes for an accurate two-dimensional study (200,862 nodes and 387,202 elements) and a y+ from 40 to 150 along the blade surfaces. By the same line of reasoning, we obtained 2.6 million hexahedral elements for the 3D study. −5 After a satisfactory convergence was obtained with a residual norm of the order of 10 , different thermodynamic quantities were extracted from the numerical solution and compared with the analytic model, as shown in Figure 15. The output torque was calculated at each time step of the two- and three-dimensional simulations, to compare the power output with the model prediction. Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 12 of 16 define the shape of the blades by interpolating over 1000 points chordwise; then the blade shape was imported in SolidWorks and simulated by Ansys Fluent. For the CFD simulation a RANS model with a k–ε turbulence model was used. For this preliminary design validation, the viscous sublayer was not resolved. Both two- and three- dimensional simulations were performed. The relative motion was simulated by the sliding mesh method in the two-dimensional case and by the Multiple Reference Frame (MRF) model in the three- dimensional case. A Conformal Mesh was generated using a growth rate of 1.045 in order to have a mesh continuity as gradual as possible with triangular elements, except in the near-wall region where rectangular elements provide better accuracy. A mesh sensitivity analysis was performed (see Figures 13 and 14) to check the independency of the results from the mesh generated. The objective function for the mesh sensitivity was the total entropy. Int. J. Turbomach. Propuls. Power 2020, 5, 19 13 of 17 Figure 13. Mesh quality before refinement. Figure 14. Mesh quality after refinement. Figure 14. Mesh quality after refinement. The above settings led to an acceptable number of nodes for an accurate two-dimensional study The above settings led to an acceptable number of nodes for an accurate two-dimensional study (200,862 nodes and 387,202 elements) and a y+ from 40 to 150 along the blade surfaces. By the same (200,862 nodes and 387,202 elements) and a y+ from 40 to 150 along the blade surfaces. By the same line of reasoning, we obtained 2.6 million hexahedral elements for the 3D study. line of reasoning, we obtained 2.6 million hexahedral elements for the 3D study. 5 After a satisfactory convergence was obtained with a residual norm of the order of 10 , di erent −5 After a satisfactory convergence was obtained with a residual norm of the order of 10 , different thermodynamic quantities were extracted from the numerical solution and compared with the analytic thermodynamic quantities were extracted from the numerical solution and compared with the model, as shown in Figure 15. The output torque was calculated at each time step of the two- and Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 13 of 16 three-dimensional simulations, to compare the power output with the model prediction. analytic model, as shown in Figure 15. The output torque was calculated at each time step of the two- and three-dimensional simulations, to compare the power output with the model prediction. (a) (b) Figure Figure 15. 15. Contours Contours of of static staticenthalpy enthalpyfor for a si a single ngle tim timeestep: step: (a (a )) 2D mode 2D model;l; ( (bb ))3D 3D model. model. Figure 16 shows the output power for both rows in all the analyzed cases. The continuous Figure 16 shows the output power for both rows in all the analyzed cases. The continuous lines lines represent the model prediction, while the dashed lines represent the 2D configuration and the represent the model prediction, while the dashed lines represent the 2D configuration and the dotted dotted lines are their averages. The asterisks represent the 3D configuration which was analyzed in a lines are their averages. The asterisks represent the 3D configuration which was analyzed in a single single frame. frame. Figure 16. Output power as a function of time. The isentropic efficiency of the whole turbine is about 0.82; however, to make a comparison with the CFD results, only the first two rows are to be considered. Therefore, the isentropic efficiency decreases to 0.759 because of the low efficiency of the first row. For the 2D configuration, the efficiency for the first two rows, calculated on the basis of the average power and enthalpy drop, is about 0.762—for the 3D configuration, about 0.751. The reduction in efficiency from the 2D to the 3D configuration is caused by the increase in the frictional losses, in particular in the hub and tip boundary layers. In fact, since the Soderberg loss model takes into account 3D losses, the slightly higher value of the efficiency in the 2D analysis compared with the “analytical” efficiency can be attributed to the neglection of wall friction (on the disks), which are instead considered in the 3D analysis. Now that the effectiveness of the proposed design algorithm has been demonstrated via the above described CFD simulations, it is possible to formulate an educated assessment of the suitability of the Ljungström turbine for different low-T thermal sources by using different working fluids in different ranges of power: Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 13 of 16 (a) (b) Figure 15. Contours of static enthalpy for a single time step: (a) 2D model; (b) 3D model. Figure 16 shows the output power for both rows in all the analyzed cases. The continuous lines represent the model prediction, while the dashed lines represent the 2D configuration and the dotted lines are their averages. The asterisks represent the 3D configuration which was analyzed in a single Int. J. Turbomach. Propuls. Power 2020, 5, 19 14 of 17 frame. Figure 16. Output power as a function of time. Figure 16. Output power as a function of time. The isentropic eciency of the whole turbine is about 0.82; however, to make a comparison with the CFD results, only the first two rows are to be considered. Therefore, the isentropic eciency The isentropic efficiency of the whole turbine is about 0.82; however, to make a comparison with decreases to 0.759 because of the low eciency of the first row. For the 2D configuration, the eciency the CFD results, only the first two rows are to be considered. Therefore, the isentropic efficiency for the first two rows, calculated on the basis of the average power and enthalpy drop, is about decreases to 0.759 because of the low efficiency of the first row. For the 2D configuration, the 0.762—for the 3D configuration, about 0.751. The reduction in eciency from the 2D to the 3D efficiency for the first two rows, calculated on the basis of the average power and enthalpy drop, is configuration is caused by the increase in the frictional losses, in particular in the hub and tip boundary about 0.762—for the 3D configuration, about 0.751. The reduction in efficiency from the 2D to the 3D layers. In fact, since the Soderberg loss model takes into account 3D losses, the slightly higher value configuration is caused by the increase in the frictional losses, in particular in the hub and tip of the eciency in the 2D analysis compared with the “analytical” eciency can be attributed to the boundary layers. In fact, since the Soderberg loss model takes into account 3D losses, the slightly neglection of wall friction (on the disks), which are instead considered in the 3D analysis. higher value of the efficiency in the 2D analysis compared with the “analytical” efficiency can be Now that the e ectiveness of the proposed design algorithm has been demonstrated via the above attributed to the neglection of wall friction (on the disks), which are instead considered in the 3D described CFD simulations, it is possible to formulate an educated assessment of the suitability of the analysis. Ljungström turbine for di erent low-T thermal sources by using di erent working fluids in di erent Now that the effectiveness of the proposed design algorithm has been demonstrated via the Int. J. Turbomach. Propuls. Power 2019, 4, x FOR PEER REVIEW 14 of 16 ranges of power: above described CFD simulations, it is possible to formulate an educated assessment of the suitability Figure 17 clearly shows that, as the range of power increases, higher density working fluids are of the Fig Ljung ure 17 clearly ström turbine showfor d s that, as the r ifferent low-T ther ange of power mal sou increases, rces by h uis ghe ingr di den fferent sity wor work king ing fluids fluids ar in e required to match the Ljungström turbine characteristics. req different r uired to m anga es of power: tch the Ljungström turbine characteristics. Figure 17. Di erent power output for di erent working fluids. Figure 17. Different power output for different working fluids. 4. Conclusions To the best of our knowledge, there is not a commercially available Ljungström blade design software, and therefore the present study proposes the use of the kinematic efficiency as a new relevant design parameter in the Ljungström design. However, as discussed in section 2, such a choice opens up a lot of points worthy of further discussion and investigation. The main reason is that this procedure does not take into account the centrifugal forces (absent in axial flows), and therefore the “optimality” of the blade profile proposed by the analytical model must be verified by an experimental validation. In the case studied here, this validation was obtained by a series of 2D time- dependent and by a steady state 3D CFD simulations. Though a full picture of the time evolution of the 3D flow is not available, nevertheless, some useful considerations can be drawn. As clearly shown in Table 9, the axial blade length is a crucial parameter in Ljungström design. By looking at Equation (17), it is clear that once all parameters, except for density, are fixed, in fact the inlet radius of each row is determined by the assumption of a fixed radius ratio, while the flow coefficient is fixed by the assumption of a constant outlet angle. In light of the above, it is also clear that the axial blade length depends on the mass flow rate and on the enthalpy drop. Following Shepherd and Kearton, it was assumed that all blades, except those in the innermost ring, have the same cross-section and they are all set similarly in the rotor. Some authors question the adequacy of these assumptions, used in the past for the design of large scale steam-powered ROT, since they lead to geometrical optimization instead of maximizing the isentropic efficiency [7]. According to these authors, if one removes the assumptions of constant βout and χ, it is possible to exploit two additional degrees of freedom for the axial blade length Equation (17). As a consequence, the formula of the kinematic efficiency changes its structure so that a multi-dimensional optimization is required. Some existing designs split the turbine into a high pressure and a low pressure stage to avoid the reduction in axial blade length in the downstream blade crowns of the turbine. This alternative solution consists of a small reduction in the mass flow rate at the inlet (first stage) and an additional injection in the low pressure wheel; however, this is not in contrast with our model, although further investigations could be useful. Clearly all the proposed changes should be analyzed in more depth because none of the above is supported by calculations. The method we propose is based on kinematic efficiency maximization, while isentropic efficiency determines the turbine performance. Even though ηkin does not take into account fluid dynamic losses, it is maximized with respect to ρr, the ratio between inlet and outlet Int. J. Turbomach. Propuls. Power 2020, 5, 19 15 of 17 4. Conclusions To the best of our knowledge, there is not a commercially available Ljungström blade design software, and therefore the present study proposes the use of the kinematic eciency as a new relevant design parameter in the Ljungström design. However, as discussed in Section 2, such a choice opens up a lot of points worthy of further discussion and investigation. The main reason is that this procedure does not take into account the centrifugal forces (absent in axial flows), and therefore the “optimality” of the blade profile proposed by the analytical model must be verified by an experimental validation. In the case studied here, this validation was obtained by a series of 2D time-dependent and by a steady state 3D CFD simulations. Though a full picture of the time evolution of the 3D flow is not available, nevertheless, some useful considerations can be drawn. As clearly shown in Table 9, the axial blade length is a crucial parameter in Ljungström design. By looking at Equation (17), it is clear that once all parameters, except for density, are fixed, in fact the inlet radius of each row is determined by the assumption of a fixed radius ratio, while the flow coecient is fixed by the assumption of a constant outlet angle. In light of the above, it is also clear that the axial blade length depends on the mass flow rate and on the enthalpy drop. Following Shepherd and Kearton, it was assumed that all blades, except those in the innermost ring, have the same cross-section and they are all set similarly in the rotor. Some authors question the adequacy of these assumptions, used in the past for the design of large scale steam-powered ROT, since they lead to geometrical optimization instead of maximizing the isentropic eciency [7]. According to these authors, if one removes the assumptions of constant and , it is possible to out exploit two additional degrees of freedom for the axial blade length Equation (17). As a consequence, the formula of the kinematic eciency changes its structure so that a multi-dimensional optimization is required. Some existing designs split the turbine into a high pressure and a low pressure stage to avoid the reduction in axial blade length in the downstream blade crowns of the turbine. This alternative solution consists of a small reduction in the mass flow rate at the inlet (first stage) and an additional injection in the low pressure wheel; however, this is not in contrast with our model, although further investigations could be useful. Clearly all the proposed changes should be analyzed in more depth because none of the above is supported by calculations. The method we propose is based on kinematic eciency maximization, while isentropic eciency determines the turbine performance. Even though  does not take into kin account fluid dynamic losses, it is maximized with respect to  , the ratio between inlet and outlet radii. The latter also plays a major role in the axial blade length Equation (17). In the end, the shorter the blade length, the larger the wall friction losses and the lower the isentropic eciency. Author Contributions: Conceptualization, U.C. and E.S.; Formal analysis, U.C.; Investigation, U.C.; Methodology, U.C.; Supervision, E.S.; Validation, E.S.; Writing—original draft, U.C.; Writing—review and editing, E.S. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Conflicts of Interest: The authors declare no conflict of interest. Int. J. Turbomach. Propuls. Power 2020, 5, 19 16 of 17 Nomenclature Symbol Description Unit of Measure Symbol Description Unit of Measure b radial chord [m] V absolute flow velocity [m/s] c specific heat [J/(kgK)] W relative flow velocity [m/s] d specific diameter wf working fluid d shaft diameter [m] wo working oil shaft axial blade length/ [m] h relative inlet angle enthalpy [J/kg] id ideal expansion ratio in inlet  blade encumbrance l blade chord [m] " deviation angle L work [J/kg]  cycle’s eciency cyc L Euler ’s work, [J/kg] [J/kg]  isentropic eciency Eul is l length scale, [m] [m]  kinematic eciency t kin m mass flow rate [kg/s]  turbine eciency turb n specific speed  leakage factor out outlet  Soderberg loss coecient p pressure [Pa]  density [kg/m ] blade to relative outlet velocity P power [kW] ratio r radius [m]  blade solidity Re Reynolds number  allowable stress [MPa] all s entropy [J/(kgK)] ' flow coecient s-s static-to-static  inlet to outlet radii ratio T temperature [T] load coecient U blade peripheral velocity [m/s] ! angular velocity [rpm] References 1. Quoilin, S.; van den Broek, M.; Declaye, S.; Dewallef, P.; Lemort, V. Techno-economic survey of Organic Rankine Cycle (ORC) systems. Renew. Sustain. Energy 2013, 22, 168–186. [CrossRef] 2. Tocci, L.; Pal, T.; Pesmazoglou, I.; Franchetti, B. Small Scale Organic Rankine Cycle (ORC): A Techno-Economic Review. Energies 2017, 10, 413. [CrossRef] 3. Kang, S.H. Design and experimental study of ORC (Organic Rankine Cycle) and radial turbine using R245fa working fluid. Energy 2012, 41, 514–524. [CrossRef] 4. Pei, G.; Li, Y.; Ji, J. Performance evaluation of a micro turbo-expander for application in low-temperature solar electricity generation. J. Zhejiang Univ. 2011, 12, 207–213. [CrossRef] 5. Palumbo, C.F.; Barnabei, V.F.; Preziuso, E.; Coronetta, U. Design and CFD analysis of a Ljungstrom turbine for an ORC cycle in a waste heat recovery application. In Proceedings of the ECOS 2016, Portorož, Slovenia, 19–23 June 2016. 6. Ljungström, F. The Development of the Ljungström Steam Turbine and Air Preheater. Proc. Inst. Mech. Eng. 1949, 160, 211–223. [CrossRef] 7. Casati, E.; Vitale, S.; Pini, M.; Persico, G.; Colonna, P. Centrifugal Turbines for Mini-Organic Rankine Cycle Power Systems. J. Eng. Gas Turbine Power 2014, 136, 122607. [CrossRef] 8. Pini, M.; Persico, G.; Casati, E.; Dossena, V. Preliminary design of a centrifugal turbine for ORC applications. In Proceedings of the First International Seminar on ORC Power Systems, Delft, The Netherlands, 22–23 September 2011. 9. Bell, I.H.; Wronsky, J.; Quoilin, S.; Lemort, V. Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library CoolProp. Ind. Eng. Chem. Res. 2014, 53, 2498–2508. [CrossRef] [PubMed] 10. Franchetti, B.; Pesiridis, A.; Pesmazoglou, I.; Sciubba, E.; Tocci, L. Thermodynamic and technical criteria for the optimal selection of the working fluid in a mini-ORC. In Proceedings of the ECOS 2016, Portorož, Slovenia, 19–23 June 2016. 11. Drescher, U.; Brüggemann, D. Fluid selection for the Organic Rankine Cycle (ORC) in biomass power and heat plants. Appl. Therm. Eng. 2007, 27, 223–228. [CrossRef] Int. J. Turbomach. Propuls. Power 2020, 5, 19 17 of 17 12. Kearton, W.J. Steam Turbine Theory and Practice: A Textbook for Engineering Students; I. Pitman & Sons: London, UK, 1958. 13. Shepherd, D.G. Principles of Turbomachineryi; Macmillan Pub Co.: New York, NY, USA, 1961. 14. Soderberg, C.R. Unpublished Notes; Gas Turbine Laboratory, Massachusetts Institute of Technology: Cambridge, MA, USA, 1949. 15. Ainley, D.G.; Mathieson, G.C.R. An Examination of the Flow and Pressure Losses in Blade Rows of Axial-Flow Turbines; HM Stationery Oce: London, UK, 1951. 16. Craig, H.R.M.; Cox, H.J.A. Performance Estimation of Axial Flow Turbines. Proc. Inst. Mech. Eng. 1971, 185, 407–424. [CrossRef] 17. Sciubba, E. Lezioni di Turbomacchine; Euroma La Goliardica: Rome, Italy, 2001. 18. Dunham, J. A Parametric Method of Turbine Blade Profile Design. In Proceedings of the ASME International Gas Turbine Conference and Products Show, Zurich, Switzerland, 30 March–4 April 1974. 19. Zweifel, O. Optimum Blade Pitch for Turbo-Machines with Special Reference to Blades of Great Curvature. Eng. Dig. 1946, 7, 358–360. 20. Pini, M.; Persico, G.; Casati, E.; Dossena, V. Preliminary design of a centrifugal turbine for ORC applications. J. Eng. Gas Turbines Power 2013, 134, 042312. [CrossRef] © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Journal

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

Published: Jul 20, 2020

There are no references for this article.