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

: 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 ﬂuid commercially adopted in powerplants: axial and, for smaller machines, radial inﬂow 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 ﬂow rate correspondingly increases along the machine channels), so that Radial Outﬂow 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 ﬁve di ﬀ erent working ﬂuids 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 ﬁxed set of operating conditions, di ﬀ erent kinematic parameters are optimized to obtain the most e ﬃ cient cascade conﬁguration. Moreover, as shown in the conclusions, the most e ﬃ cient blade proﬁle leads to higher friction losses, making further investigation regarding the best conﬁguration necessary.


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 efficiency 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 efficiency 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

Fluid Selection
In this work, a hydro-treated mineral oil, Paratherm NF, is the heat transfer medium taken as 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 water, which enters the condenser at room temperature (25 °C or 298 K) and is assumed to be available in the necessary quantity. The selection of the working fluid and the identification of the thermodynamic cycle was performed via a MATLAB code by using the CoolProp library [9] to calculate the fluid properties. A first selection of fluids from the CoolProp library was performed under two criteria: first, since the source is at low temperature, the fluid should have a low critical pressure [10,11]; second, the Ljungström turbine requires relatively high volume flow rates at its inlet and a well-identified variation in the density along the radius due to its peculiar cross-sectional area, and thus the fluid ought to have a low density at turbine inlet (i.e., high values of the coefficients in the Redlich-Kwong equation of state). According to these considerations the following fluids were considered: an alicyclic hydrocarbon, Cyclopentane (C5H10); an azeotropic mixture, SES36; a haloalkane refrigerant, R134a (CH2FCF3); an alkane refrigerant, R601 (C5H12); a hydrofluorocarbon refrigerant, R245fa (C3H3F5).
For each working fluid, as shown in Figure 2, a specific thermodynamic cycle was designed by fixing the top and bottom pinch points and specifying a 10 °C water temperature rise in the condenser. For a correct comparison, the thermal power input is the same for all cycles considered. The heat input and output in the boiler and the condenser were then matched by iteratively adjusting the working fluid mass flow rate or/and the evaporation temperature. Therefore, the turbine enthalpy drop and the pump work were determined with an assumed 0.8 efficiency for the pump and 0.83 (polytropic) for the turbine. The overall cycle efficiency is equal to:

Fluid Selection
In this work, a hydro-treated mineral oil, Paratherm NF, is the heat transfer medium taken as 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 water, which enters the condenser at room temperature (25 • C or 298 K) and is assumed to be available in the necessary quantity. The selection of the working fluid and the identification of the thermodynamic cycle was performed via a MATLAB code by using the CoolProp library [9] to calculate the fluid properties. A first selection of fluids from the CoolProp library was performed under two criteria: first, since the source is at low temperature, the fluid should have a low critical pressure [10,11]; second, the Ljungström turbine requires relatively high volume flow rates at its inlet and a well-identified variation in the density along the radius due to its peculiar cross-sectional area, and thus the fluid ought to have a low density at turbine inlet (i.e., high values of the coefficients in the Redlich-Kwong equation of state). According to these considerations the following fluids were considered: an alicyclic hydrocarbon, Cyclopentane (C5H10); an azeotropic mixture, SES36; a haloalkane refrigerant, R134a (CH2FCF3); an alkane refrigerant, R601 (C5H12); a hydrofluorocarbon refrigerant, R245fa (C3H3F5).
For each working fluid, as shown in Figure 2, a specific thermodynamic cycle was designed by fixing the top and bottom pinch points and specifying a 10 • C water temperature rise in the condenser. For a correct comparison, the thermal power input is the same for all cycles considered. The heat input and output in the boiler and the condenser were then matched by iteratively adjusting the working fluid mass flow rate or/and the evaporation temperature. Therefore, the turbine enthalpy drop and the pump work were determined with an assumed 0.8 efficiency for the pump and 0.83 (polytropic) for the turbine. The overall cycle efficiency is equal to:  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 sensible heat content of the working fluid before it enters the condenser. The percentage of thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6.   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 sensible heat content of the working fluid before it enters the condenser. The percentage of thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6.     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 sensible heat content of the working fluid before it enters the condenser. The percentage of thermal power recovered by the regenerator amounts to 7.98% of the whole thermal input for the Cyclopentane, 14.50% for the R601, 11.23% for the R134a, 15.99% for the R245fa and 21.42% for the SES36. Efficiencies and mass flow rates (indicative of expander size) are reported in Table 6.

The Optimization Algorithm
The standard formulation of the Euler work (2) shows that, for a Ljungström "stage", the only positive contribution to the specific work is given by the increase in the relative velocity, because u out is only slightly higher than u in : 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: • 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: • In all blade crowns, the ratio of the radial chord of the blade to the outlet radius is constant and equal to: • Due to the axial symmetry of the turbine, the absolute velocity at the inlet of the first row can be considered completely radial. The above assumptions suffice for a definition of kinematic efficiency. 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: That can be rewritten as: 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: The parameters to be considered to maximize the kinematic efficiency are, therefore, ρ r , χ, β out and β 0 . Due to the second and third assumptions it is not possible to maximize both η kin,i and η kin,1 with respect to ρ r and χ, thus only the η kin,i equation is considered here for the optimization.
It can be easily verified that the derivative of the kinematic efficiency (7) is linearly dependent on both χ and β out and thus it is not possible to maximize this equation with respect to these parameters. This conclusion is reinforced by a closer inspection of Equation (7), which is maximized either when the numerator of the second term tends to zero or the denominator to a maximum. Both cases leading to ideal turbine configurations: • χ = 1, meets the maximum of the efficiency for ρ r = cos(β out )/2, and causes the numerator to be null. In fact, in this case u in = u out and v in = v out so that the velocity triangle is similar to that of an impulse turbine; • β out = 0 causes the denominator to be at its maximum, in which case w in = 0 because the triangle collapses onto a line.
Therefore, the only parameter to be considered in the velocity triangles optimization is ρ r , which allows us to maximize the kinematic efficiency for different values of β out and χ, as shown in Figure 8. It can be easily verified that the derivative of the kinematic efficiency (7) is linearly dependent on both χ and βout and thus it is not possible to maximize this equation with respect to these parameters. This conclusion is reinforced by a closer inspection of Equation (7), which is maximized either when the numerator of the second term tends to zero or the denominator to a maximum. Both cases leading to ideal turbine configurations: • χ = 1, meets the maximum of the efficiency for ρr = cos(βout)/2, and causes the numerator to be null. In fact, in this case uin = uout and vin = vout so that the velocity triangle is similar to that of an impulse turbine; • βout = 0 causes the denominator to be at its maximum, in which case win = 0 because the triangle collapses onto a line.
Therefore, the only parameter to be considered in the velocity triangles optimization is ρr, which allows us to maximize the kinematic efficiency for different values of βout and χ, as shown in Figure 8. Considering the above, we can now proceed to the study of the isentropic efficiency, starting from its equation and the loss model adopted here: Considering the above, we can now proceed to the study of the isentropic efficiency, starting from its equation and the loss model adopted here: In order to obtain the value of w out,id , Soderberg [14] loss correlation was used: And after some rearrangement, the formula of the isentropic efficiency for the i-th row, with the exception of the first one, is obtained: While for the first row: 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 abovementioned prediction models. It was found that the adoption of different loss models does not noticeably influence the results, as shown in Figure 9. abovementioned prediction models. It was found that the adoption of different loss models does not noticeably influence the results, as shown in Figure 9. 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.
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 triangles. 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 .
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 triangles.
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.
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 triangles. 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. 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. 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 r in and β in that satisfy the continuity equation for all the l > l min .
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 (∆h i ).
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 in Figure 11.
Since the maximum efficiency is reached for b = 0 and β out = 0 and since these conditions are 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: where the leakage factor ξ was assumed here to be 0.98.
it is now possible to calculate the loss coefficient ξr and the enthalpy drop for the first stage (∆hi).
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 in Figure 11. Since the maximum efficiency is reached for b = 0 and βout = 0 and since these conditions are 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: where the leakage factor ξ was assumed here to be 0.98. Figure 11. Schematic representation of the design algorithm.
The above equation has two degrees of freedom. By fixing a proper range of values for β 0 and r in all possible combinations that meet the equation can be calculated. These ranges are constrained by two different conditions. β 0 must be limited in order to obtain w out,1 > w in,1 : while r in for the first row must be larger than the shaft radius that was calculated with the following formula: 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 coefficients 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.
As shown in Table 8 and, as expected, the χ parameter is strictly correlated to the efficiency of the turbine. A closer look to these results reveals that the less efficient thermodynamic cycles presents the higher isentropic efficiencies. The Ljungström turbine performs better with low enthalpy drop and high-volume flow rates (higher specific speed n s ). Moreover, since the cross-sectional area depends on the volume flow rate variation, not all the selected fluids lead to geometrical characteristics, such as 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.  In the present study, R601 was chosen as working fluid because, despite the slightly lower efficiency 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:

CFD Validation
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] and the number of blades was selected using a modified Zweifel's criterion [19,20], since ROTs have an axial blade shape. A MATLAB code based on the above parametric method was written to

CFD Validation
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] and the number of blades was selected using a modified Zweifel's criterion [19,20], since ROTs have an axial blade shape. A MATLAB code based on the above parametric method was written to 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. 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 threedimensional 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 threedimensional 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.    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.
After a satisfactory convergence was obtained with a residual norm of the order of 10 −5 , 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 twoand three-dimensional simulations, to compare the power output with the model prediction. 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.
After a satisfactory convergence was obtained with a residual norm of the order of 10 −5 , 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.    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 frame. 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 frame. 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: 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: Figure 17 clearly shows that, as the range of power increases, higher density working fluids are required to match the Ljungström turbine characteristics.

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

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 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 efficiency.

Conflicts of Interest:
The authors declare no conflict of interest.