Three-dimensional unsteady stator-rotor interactions in high-expansion organic Rankine cycle turbines

Organic Rankine cycle (ORC) systems are a readily available technology to convert thermal energy from renewable- and waste heat sources into electricity. However, their thermal performance is relatively low due to the low temperature of the available heat sources, but more importantly, due to the low ef ﬁ ciency of the employed expander. Designing the turboexpander is exceptionally challenging, because the ﬂ ow ﬁ eld is highly supersonic and unsteady, and since the expansion takes place in the highly non-ideal dense-vapor region. In this work, we perform unprecedented three-dimensional unsteady simulations of several high-expansion cantilever ORC turbines to highlight distinctive loss mechanisms. The simu- lations indicate strong unsteady effects in the rotor blade passage, as a result of unsteady propagating shock waves interacting with viscous wakes and boundary layers. Moreover, the ﬂ ow ﬁ eld in the rotor blade passage is strongly affected by three-dimensional secondary ﬂ ow features and a sharp expansion in the shroud region at the inlet of the rotor blade. These span-wise mechanisms and unsteady ﬂ ow in- teractions introduce irreversible losses which must be taken into account for designing highly ef ﬁ cient ORC expanders. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). considering three-dimensional numerical domains. Most computational ﬂ uid dynamic (CFD) studies of ORC


Introduction
ORC power systems are a viable alternative to convert low-tomedium grade heat sources at temperatures between 120 and 500 + C [1] to electrical power. The thermal energy available from biomass combustion, solar radiation, geothermal reservoirs, or waste heat from industrial processes, is arguably immense. Still, the low temperature of these heat sources is the main drawback for their conversion into electricity. While the operating principle of a steam Rankine cycle is comparable, the working fluid of an ORC is an organic compound with a high molecular-weight. As a result, the fluid can be selected to best match the temperature of the heat source, adding a degree of freedom to the system design, yielding a higher thermal efficiency compared to a steam Rankine cycle for the same conditions [1]. In terms of overall efficiency, the expander is the most critical component because of its direct impact on the overall system performance [2].
The expansion process of an ORC is distinct from its Rankine or Brayton cycle counterpart. First, the expansion takes place close to the working fluid's critical point, in the dense-vapor region, where the ideal gas assumption is invalid. Consequently, complex equations of state (EoS) are necessary to accurately describe organic fluids [3e7]. Second, an ORC expander features a small specific enthalpy drop due to the molecular complexity of the organic fluid. Therefore, the designer can choose a few numbers of stages (one or two) without the expander experiencing efficiency drop or high rotational speed [1]. As a consequence, large pressure ratios per stage can be achieved, as large as 100 in a number of applications. However, the large expansion ratio and the low speed of sound of organic fluids d relative to air or steam d lead to supersonic flows with inevitable shock waves and complex flow features. For these two reasons, mean-line flow models and diagrams derived from experimental data (e.g., Balje diagram), which are commonly used for steam and gas turbines, are not highly reliable for designing ORC turbines [8].
Several efforts in improving the predictive capability and understanding the flow features have been conducted, which are summarized in Table 1. This literature survey is expanded in the following paragraphs, structured in terms of simulation type, for instance, time-independent/dependent simulations or estimations considering three-dimensional numerical domains.
Most computational fluid dynamic (CFD) studies of ORC turbomachinery are under the assumption of steady-state with a mixing plane boundary condition at the interface between the stator and the rotor [10e14,17]. Because of the averaging at the interface and the assumption of periodicity in the tangential direction, only one blade passage is simulated per rotor d independent of the number of blades of the turbine stage d making the simulation computationally less expensive. However, because of the mixing plane assumption, the unsteady information between the stator and the rotor is lost. For instance, the impact of the viscous wake emanating from the trailing edge of the stator over the rotor blade can not be quantified [22]. An alternative is to use the frozen rotor approximation. Still, these simulations are physically inconsistent due to the absence of unsteadiness. These types of simulations are unable to accurately estimate unsteady effects that bring insights on otherwise unmeasured loss mechanisms [15,16,18e21].
To evaluate large expansion ratio ORC turbines (PR > 10) operating with an organic compound with a relatively low speed of sound, it is of paramount importance to consider time-resolved unsteady simulations. Such machines have a transonic/supersonic flow at the exit of the stator. This flow regime enhances the already inherently unsteady flow between the stator and rotor due to the formation of shock waves, which interact with the boundary layer and enhance flow mixing. Several studies that compared both steady and unsteady simulations of an ORC turbine stage with a pressure ratio in the order of 15 are [15,16,20,21]. All these studies reported a significant performance drop in the unsteady simulation if compared to the steady computation. In particular, Rubechini et al. [21] recognized a nearly linear relationship of the relative error between the two computations (steady and unsteady) with the nozzle discharge Mach number. To the authors' knowledge, Rinaldi et al. [18] is the first available open literature, which modelled the stator/rotor interaction of an ORC turbine with a PR greater than 100. From this study, the unsteady losses resulted in a 1.4% drop in efficiency. It is obvious from these studies that unsteady effects influence the predicted performance of an ORC turbine with a large pressure ratio.
Three-dimensional (3D) effects are another flow feature in ORC turbines that are frequently overlooked, as only a few studies have accounted for them in their CFD simulations [11,13,14,16]. Most papers consider a two-dimensional (2D) or a quasi-threedimensional (Q3D) flow, see Ref. [17,18]. An article that does consider the full 3D geometry of a high-expansion ORC turbine is Harinck et al. [14]. This scientific paper utilized a commercial software, ANSYS CFX, to perform a steady-state simulation of a cantilever ORC turbine. Marconcini et al. [20] investigated the unsteady interaction of a radial inflow ORC turbine (PR ¼ 20) by performing a full annulus simulation. More recently, Vitale et al. [11] performed a 3D simulation of an optimized rotor blade of a radial outflow ORC turbine. Other authors have also performed 3D simulations of ORC turbines but mainly consider standard radial inflow turbines with an expansion ratio of less than 11, see for example [12,13,16].
To the authors' knowledge, the present study is the first CFD investigation of a ORC turbine stage with a pressure ratio larger than 100 which accounts for the combination of real gas effects, unsteady stator/rotor interaction, and three-dimensional effects. We consider two different rotor blade geometries, while the geometry of the stator nozzle is left unchanged. We refer to the rotor geometry as (1) the "old" blade design, which closely resembles the geometry from an ORC manufacturer, and (2) a "new" design, which has been obtained by an indirect method that minimizes the change of flow turning and the cross-sectional area in the rotor passage. The novelty of this work is that it shows for the first time a detailed analysis of the three-dimensional unsteady phenomena (shock waves, viscous wakes, and shockwave-boundary layer interaction) in high-expansion ORC turbines.
The structure of the present article is as follows: first, an overview of the system is presented with a brief description of the two ORC expander geometries. Later, we introduce the numerical infrastructure, which includes the computational domain, and the numerical methods and models. We discuss the results of the CFD simulations in section 4; which consists of the unsteady and threedimensional effects. Subsequently, we make a quantitative comparison between the two geometries and the numerical domains in terms of the entropy generation and turbine performance. The final section of this article presents conclusions and future works.

System: ORC expander
In the present study, a cantilever turbine (radial inflow-radial outflow turbine) is studied for thermodynamic conditions representative of a commercial ORC unit. The electrical power output of the ORC turbogenerator, manufactured by Tri-O-Gen BV in the Netherlands, is in the range of 150e200 kW. Fig. 1 (a) shows a photograph of Tri-O-Gen's turboexpander which consists of 18 stator nozzles and 43 rotor blades. The expander has a PR > 100 and operates with toluene as working fluid with inlet conditions close to the critical point. The reduced temperature and pressure are in the range of T red ¼ ½0:9; 1:0 and p red ¼ ½0:75; 0:85, respectively.
Moreover, the rotational speed is in the range of u ¼ ½400; 500 Hz; for the current study, we apply a rotational speed of 430 Hz (z26 krpm). Table 1 Literature on numerical simulations of ORC rotors/turbines with information related to the working fluid, the expander type, the pressure ratio (PR), the simulation approach, and the flow regimen between the stator and rotor.

Turbine stage geometry
As mentioned earlier, we consider two rotor blade geometries and maintain the stator nozzle constant for two investigated turbine stages. The investigated rotor blade geometries consist of an old rotor blade, similar to the current blade of Tri-O-Gen's ORC system, and an optimized blade geometry [23]. Both of these rotor blade geometries are analyzed in the present article using unsteady 3D CFD simulations with the numerical domain depicted in Fig. 1(b). The design methodology for the stator and the rotor blades is based on an in-house developed design methodology and parametrization, and is briefly discussed hereafter.
Stator: First, the diverging section of the stator vane is designed using the Method of Characteristics (MOC) for an axial configuration for fluids near the critical point [5]. This supersonic axial nozzle is then converted into a radial configuration using conformal mapping to conserve the area ratio and the outflow angle. The reader is referred to Refs. [7] for more details on the stator vane design methodology.
Old rotor blade: The geometry of the old blade is determined using a direct approach (without any aerodynamic optimization) with ten geometrical parameters. We investigate numerically a blade that closely resembles the turboexpander geometry from Tri-O-Gen. In span-wise direction, the blade thickness varies due to strength and vibrational constraints in the existent machine. New rotor blade: We use an inverse design approach coupled with an optimization procedure to generate the new rotor blade shape [23]. The first step is purely geometrical and follows the idea to design a rotor blade passage that has a smooth turning flow and a steady increase in the cross-sectional area. This methodology is then applied in an aerodynamic optimization, where the design parameters are the blade angles, stagger angle, number of blades, among others. The objective of the optimization is to minimize the outlet total enthalpy, which essentially results in an increased power output of the turbine stage. We did not aim to maximize the turbine's total-to-static efficiency because the outflow of the domain is supersonic: the outlet boundary static pressure is not fixed. Therefore, the algorithm could inadvertently decrease the pressure ratio of the turboexpander by optimizing the efficiency. The fitness of each blade is estimated by means of Q3D simulations using SU2 (CFD solver) [24] with a mixing plane approximation between the stator and rotor; the mixing plane surface includes a non-reflecting boundary condition as presented by Giles [25]. Finally, the optimal blade shape is mechanically analyzed, which resulted in the conclusion that the blade thickness in span-wise direction can remain constant. For further information on the new blade design methodology, the reader is referred to Ref. [23]. Tri-O-Gen B.V. is planning to manufacture this blade and test it in the near future.
Cantilever turbines, which are usually of low reaction [26], require an increase of the cross-sectional area due to the expansion of the medium. However, an increase of the cross-sectional area in the blade passage is difficult to accommodate due to the high volume flow ratio of organic fluids and due to the radial inflow configuration. The blade height distribution can hence be used as a free parameter to control the cross-sectional area distribution in the rotor passage.
In Fig. 2, we depict the resulting blade profile and height distribution from the two blade design methodologies. This figure compares the new blade d including the shape optimization d with the old one. The blade thickness for the new blade has increased and the number of rotor blades has increased as well from 43 to 47. However, to achieve a steady increase in the crosssectional area, the blade height distribution has also increased, such that the trailing edge height of the new blade is almost 15% larger than the one of the old blade.

Numerical domain
Since this study aims to investigate three-dimensional unsteady effects in a high-expansion ORC turbine, we perform two types of simulations: Q3D and fully 3D simulations. Note, to have a reasonable comparison, the grid topology and number of mesh cells from the Q3D mesh is also used for the 3D mesh at mid-span and then projected towards the respective end walls.
In terms of the Q3D geometry, the cross-sectional area increases towards the machine axis to capture the same expansion as for the 3D geometry. We depict the Q3D numerical domain as a function of the radius for both rotor blades in the right plot of Fig. 2.
The 3D simulations account for the whole blade passage in span-wise direction, which is discretized by 50 cells, clustered near the walls at the hub and the shroud, achieving a y þ < 1 for most of the hub and shroud; however, close to the rotor leading edge at the shroud the simulations achieved a y þ max < 10. Even though a shroud leakage is present in the real machine, the tip-clearance is out of the scope in this research; therefore, we simulate shrouded blades.
In all simulations, the nozzle and blade boundary layers are discretized with O-type structured meshes. Their elements are also clustered at the wall, controlled by a hyperbolic tangent function, to achieve a y þ z1. To ensure mesh-independent results, we performed a Q3D simulation on a fine grid, which was determined using a grid refinement strategy as done in our previous works for the same turbo-expander configuration [7,18,23]; the grid used has z 45,000 elements while the fine grid z 68,000 elements. Grid convergence was evaluated in terms of turbine performance, mass flow rate, and outlet conditions, resulting in a relative difference of < 0.2% d for all parameters d between the fine grid and the mesh used in our simulations.
Several approximations are available in the literature to alleviate the computational demand of a 3D unsteady simulation. The most common are the phase-lagged boundary conditions [27], harmonic balance [28], and changing the blade count. While phase lagged can take a long time to converge, the harmonic balance has drawbacks in handling discontinuities at boundaries. We chose to scale the rotor blade count as it represents a trade-off between accuracy and computational cost. The number of rotor blades was decreased to 36, such that the numerical domain consists of two rotor blades per one stator nozzle, as illustrated in Fig. 1(b). The blade count reduction changes the solidity, which leads to inaccuracies in the blade loading and loss mechanism.

Solver and models
An in-house CFD code is used, which discretizes the compressible Navier-Stokes equation using the finite volume formulation. This in-house code has been extensively validated for turbomachinery flows in previous studies [18,29e31] and more details about the implemented numerical methods and models can be found in Refs. [32]. In this work, all simulations are performed with a second-order accurate spatial and temporal discretization scheme.
The fluid properties for toluene are obtained with the multiparameter EoS of Lemmon and Span [33] and tabulated using a look-up table (LUT) approach as discussed in Rinaldi et al. [18]. The ranges for temperature and density are (T ¼ ½310; 390 K) and (r ¼ ½0:005; 145:0 kg m À3 ), respectively, using 400 points in each direction. The LUT tabulates pressure, entropy, enthalpy, viscosity, thermal conductivity, and the speed of sound. The computational time with the LUT is drastically reduced, more than 20 times compared to the direct evaluation of the EoS, and the maximum interpolation error is below 0:01% [18].
Turbulence has been modelled with the one-equation turbulence model of Spalart-Allmaras (SA) [34], in accordance to previous authors [16e18, 20,21]. Moreover, in our previous study on turbulent flows [35], the SA model gave accurate results in fullydeveloped channels with non-ideal gases if compared to direct numerical simulations.

Boundary conditions
The boundary conditions are specified as follows. Uniform distributions of total pressure, total temperature, eddy viscosity and velocity angle are prescribed at the inlet of the computational domain. At the outflow of the rotor domain the static pressure is prescribed, as long as the radial outflow velocity is subsonic, otherwise, Neumann boundary conditions are specified for all quantities. Standard periodic boundary conditions are applied along the circumferential direction. At all walls no-slip boundary conditions are specified for the velocity and the eddy viscosity is set to zero. Adiabatic boundary conditions are applied for the wall temperature, which includes the stator vane, the rotor blades, and the turbine's hub and shroud of the 3D simulations. For the Q3D simulation, a symmetry boundary is set at the top and bottom of the domain.
To model the unsteady stator-rotor interaction, we use a fullyconservative flux-assembling technique. An auxiliary mesh is generated at the interface between the stator and rotor by intersecting the elements of the two non-conformal grids. The conservation of the flow quantities is ensured by calculating the fluxes for the elements of one mesh and adding their contribution to the respective control volumes of the other mesh. During the unsteady simulation, the auxiliary mesh is updated each time step to obtain a new auxiliary mesh topology. The reader is referred to Ref. [36] for more details on the flux-conserving technique. The non-matching mesh technique was extended in this research to handle 3D geometries, like the cantilever ORC turbine stage.

Flow field analysis of the CFD simulations
This section reports the outcome of all unsteady simulations for the old and the new blade geometry. For each geometry a Q3D and 3D simulation has been performed, which results in four simulations in total. We first examine the unsteadiness of the flow field, followed by an analysis of the three-dimensional effects.
The number of time steps per period (two-blade pass) is 60,000 and 80,000 for the Q3D and 3D simulation, respectively. An

Unsteady flow interaction
The time evolution of the relative Mach number for the Q3D simulations of the old and new rotor blades are depicted in Figs. 3 and 4, respectively. We have limited our analysis to half the angular period. In these plots, and most contour plots in this paper, the isobars are depicted as continuous black lines.

Instantaneous flow field
The flow features in the stator passage are equivalent for both blade geometries and are typical of a radial inflow supersonic stator vane. Because of the large pressure ratio, the flow is accelerated to sonic conditions at the throat of the nozzle and later becomes supersonic in the divergent section. A highly supersonic flow is reached at the end of the nozzle, with a Mach number of around 2.7. As the flow enters the free expansion region between the stator and the rotor, the supersonic flow can not maintain its direction due to the flow coming from the adjacent nozzle. Therefore, two oblique shock waves, ST 1 and ST 2 in Figs. 3(d) and 4(a), emanate from the stator trailing edge (STE). The shock wave ST 1 directly enters the rotor passage, while ST 2 hits the neighboring stator wall. This shock wave is then reflected (R) towards the rotor as well, as can be seen in the pressure gradient contour in Figs. 6(a) and 7(a). Both shock waves disturb the flow field downstream. However, these are weak oblique shocks and do not generate strong recompression in the semi-bladed region between the stator and rotor. Thus the design of the stator blade achieves the desired expansion ratio.
The old rotor blade suffers from several flow phenomena that increase losses, see Fig. 3. A bow shock is formed at the rotor leading edge (RLE) because the flow is supersonic in the relative frame of reference (relative Mach number is approximately 1.25). Moreover, the flow becomes detached at the suction side (SS) of the blade inducing a recirculation bubble and another oblique shock wave. The flow detachment is likely produced by the strong curvature of the suction side and an excessive incidence angle at the RLE. Because of the increase of the cross-sectional area in the blade passage, the flow continues to accelerate; therefore, two oblique shocks emanate from the rotor trailing edge (RTE). Two obvious solutions to reduce losses in the rotor passage are (1) a thicker blade that would reduce the flow recirculation on the SS, and (2) a higher rotational speed to have a subsonic flow in the relative frame of reference at the inlet of the rotor. The former solution is considered by the new rotor blade as seen in Fig. 4.
The new blade reduces the separation bubble at the SS of the blade. The flow detachment is likely to be produced by the strong curvature of the SS of the blade. Moreover, the new blade has a more uniform expansion, if compared to the old blade, along the rotor passage, as depicted by the isobars in Fig. 4. Qualitatively, the new blade considerably improves the flow field inside the highexpansion ORC turbine; even though the flow field is highly supersonic downstream of the nozzle.

Time evolution
Figs. 3 and 4 also reveal unsteady nature of the flow between the stator and the rotor, and in the rotor blade passage. To avoid redundancy, we will only discuss the unsteady shock wake interaction for the Q3D simulation.
Following three unsteady structures within the old rotor blade passage one can observe the viscous wake (VW) of the STE, the bow shock (BS) at the RLE, and the oblique shock (OS) emanating from the flow separation of the SS of the blade; all named in Fig. 3(a). The blade rotation influences the shock waves (ST 1 and ST 2 ) and VW from the STE, generating a highly unsteady flow between stator and  of the blade is also unsteady, e.g., at t ¼ 2=12 and t ¼ 3= 12, we observe shock waves induced by the separation bubble of the SS of blade 1. These shock waves are created by the interaction of the VW with the flow separation; this is observed by following the VW that impinges at blade 2 RLE at t ¼ 0=12 and continues to disturb the flow at the PS of blade 2 until t ¼ 3=12. Between t ¼ 3= 12 and t ¼ 4=12, the OS at the SS of blade 1 is formed as a consequence of the flow detachment. This OS interacts with the shock waves induced by the separation and the BS at the RLE of blade 2 (t ¼ 5= 12); this is how the complex shock-shock interaction, already discussed above, is formed. It is important to note that this highly nonuniform flow field and these unsteady interactions can only be captured with unsteady simulations. Accordingly, unsteady simulations are necessary to account for their associated loss mechanism in the stage performance. The time-averaged, maximum, and minimum values of the pressure distribution on the old rotor blade are displayed in Fig. 5 to quantitatively analyse the unsteadiness at the rotor blade surface. Ideally, the blade loading along the surfaces should be as high and uniform as possible. However, we interpret large oscillations of the blade loading as the minimum and maximum pressure largely deviate from the time-averaged values; the maximum pressure differs from the time-averaged value by around 1 bar at the RLE pressure side. This pressure oscillations are directly related to the bow shock at the RLE, while the flow detachment of the rotor blade is visible between a blade chord distance of 0.5 and 0.9. A significant variation is seen at this location between the maximum and minimum pressure loading. As such, the forces acting on the blade are highly unsteady and can be included in the mechanical and vibrational analysis of the rotor blade loading.  The highly unsteady flow field of the new rotor design is illustrated in Fig. 4. The oblique shock waves stemming from the STE (ST 1 and ST 2 ) and its reflection (R) intensifies the non-uniformity of the flow in the region between stator and rotor and in the rotor passage. Because the rotor inlet flow is still supersonic in the relative frame of reference, shock waves are produced in the rotor blade passage. Three shock waves are depicted in the rotor blade passage named in Fig. 4(d): (1) a bow shock (BS) produced at the RLE because the flow is supersonic in the relative frame of reference (Ma rel z1:3), (2) a separation shock wave (SW) generated by the flow detachment at the suction side of the blade, and (3) the reflected shock (R) from the tail of SW impinging at the PS of the rotor blade. These three shock waves rotate with the blade developing a complex shock structure at the PS, which is highly unsteady, see from t ¼ 0=12 to t ¼ 2=12. Subsequently, at t ¼ 3= 12, R and the tail of SW interact with the BS until this shock-shock interaction merges into a single shock at t ¼ 4=12. Near the reattachment region at the SS, an expansion fan (E) is being generated, see Fig. 4(d); the flow recompresses in this region. The flow continues to accelerate after the separation bubble achieving supersonic conditions again. Therefore, a fish-tail shock waves emanate from the RTE.

Qualitative comparison between the Q3D and 3D simulations
There are distinct similarities and differences between the flow field of the Q3D and 3D (at mid-span) simulations for both rotor blade designs, which we will highlight by plotting the instantaneous pressure gradient (Figs. 6 and 7) and the time-averaged entropy contours of the Q3D and 3D simulations at mid-span ( Fig. 8(a) and (d) for the old blade, and Fig. 8

(b) and (g) for the new blade).
Concerning the similarities, both simulations (Q3D and 3D) depict equivalent phenomena, e.g., the viscous wake from the STE, the flow separation at the SS of the blade (compare Fig. 8(a) and (d) of the old blade and Fig. 8(b) and (g) of the new blade), and complex shock-shock interaction at the rotor passage (see Fig. 6(aeb) and Fig. 7(aeb)). No qualitative differences are seen in the stator flow field between the two simulation type.
Concerning the differences, the 3D simulation at mid-span shows a thicker wake and larger thermodynamic entropy increase within the recirculation bubble, a consequence of the spanwise effects, and a change in the strength of the bow shock at the RLE. For the old blade, the bow shock wave is stronger than for the new design in the 3D simulations. Both rotor blades suffer from span-wise effects, which we will investigate in more detail in the following sub-section. Fig. 8(c)-(e) and 8(f)e(g) describe the time-averaged entropy contours of the 3D simulation for both rotor designs at three spanwise locations: 10% (near the hub), 50% (at mid-span), and 90% (near the shroud). Large differences in the flow field are observed in the span-wise direction at the rotor passage. Three important flow features that change along the span-wise direction are: (1) the size

Old rotor blade.
The strongest three-dimensional effects are seen within the rotor passage. Contrary to the mid-span, the entropy increase at the hub is not restricted to the recirculation bubble, compare for example the rotor passage in Fig. 8(c) and (d).
The flow field near the blade shroud has the highest difference if compared to the mid-span. The flow detachment changed location relative to the chord distance of the blade; the separation bubble at the shroud is closer to the RLE, see Fig. 8(e). Moreover, the shape of the flow separation at the shroud is distinct. Fig. 10(a) depicts the reason for these differences: a secondary flow caused by the strong flow detachment at the suction side and the high flaring angle of the blade height distribution. The high flaring angle is a consequence of the radial inflow configuration, and the high volume flow ratio. Therefore, the height distribution needs to increase rapidly in the rotor passage to compensate for these two factors. In this particular case, the high flaring angle in the radial-to-axial bend drives low momentum fluid towards the shroud because of flow recirculation and the meridional stream-wise curvature in the span-wise direction. A sharp increase in the averaged entropy d due to the secondary flow d is seen just at the inlet of the rotor passage for the shroud section. The 3D effects are actively modifying the flow field inside the high-expansion ORC turbine; therefore, the designer needs to consider these span-wise effects during the design phase by adapting both the blade profile and height distribution.

New rotor blade.
Just as in the old blade design, the threedimensional effects become apparent in the rotor passage as shown by the time-averaged entropy contour of Fig. 8(f and g) and the instantaneous relative Mach number contour in Fig. 9 at different span-wise locations. But unlike in the old blade design, there is a sudden expansion at the inlet of the rotor passage near the shroud, check the rotor inlet of Figs. 8(f) and 9(e). This pressure change is related to the height distribution, see right plot of Fig. 2, which has an abrupt blade height increase at this radius. The bow shock at the RLE is a three-dimensional structure as depicted in Fig. 7(c). The shock wave structures within the rotor passage are stronger at mid-span, given the absence of near-wall effects. Compared to the 25% span-wise location, there is barely any difference in the relative Mach number contour, compare Fig. 9 (b) and (c). However, there is a significant disparity with the rest of span-wise locations, see Fig. 9(a), (d) and (e). For the flow near the hub (10%) and the shroud (90%), the bow shock at the RLE is much weaker, which can be attributed to near-wall effects. Furthermore, this 3D bow shock wave affects the blade loading; there is a sharp decrease in pressure at the PS in Fig. 11 near the inlet of the rotor.
There is flow separation in all span-wise contours in the rotor passage, but differ due to the 3D effects close to the casing and shroud: At 90% span, the flow detachment is barely present. The entropy increase due to near-wall effects from the shroud have a larger influence in the flow field at this region. As mentioned already, the pressure at the inlet of the rotor passage is less (z0:3Dbar) at the shroud than at other span-wise locations. This phenomenon is also clearly seen in the blade loading near the shroud at a scaled radius of 0.97 (Fig. 11). Near the hub (10%), the separation creates a sharp increase in thermodynamic entropy and has a distinct shape, different from the mid-span and near the shroud, see Fig. 8(h); a secondary flow d illustrated by Fig. 10(b) d causes this sharp entropy increase at the hub.
The flow in the rotor passage suffers from a secondary flow because of a significant increase in shroud height. The sharp Fig. 9. Snapshot of the relative Mach number contour with isobars at five span-wise locations (10%, 25%, 50%, 75%, and 90%) at t ¼ 4=12 for the 3D simulation of the new rotor blade. expansion at the shroud of the rotor inlet drives fluid to bend from radial to the axial direction; still, not all the fluid manages to turn, especially near the hub, as depicted by Fig. 10(b). The new blade leads to reduce viscous losses by decreasing the separation bubble at the SS of the blade; still, strong span-wise effects are generated in the rotor passage.

Quantitative analysis of the turbine stages
We now quantitatively compare the two rotor blades, utilizing both numerical domains: Q3D and 3D, by analyzing the accumulative rate of entropy generation (Fig. 12) and the time-averaged performance ( Table 2). In Fig. 12, we depict the accumulated rate of entropy generation based on span percentage to have a fair comparison between the Q3D and the mid-span section of the 3D. The rate of entropy generation per unit volume, see Ref. [37], is calculated as: where l and t ij are the thermal conductivity and the shear stress, respectively.

Accumulative rate of entropy generation
The accumulative entropy generation analysis of the ORC turboexpander can help us to clearly identify the sources of losses. Furthermore, by comparing the mid-span simulation with its 3D counterpart, we can also discover which potential losses are concealed by a mid-span design approach. This section consists of four parts. First, the Q3D simulations are compared, followed by a comparison between the Q3D and 3D simulations at mid-span for the stator and rotor, respectively. And finally, the 3D simulations for both geometries are analyzed.

Q3D simulations
By examining the accumulative entropy generation for the two Q3D simulations in the right plot of Fig. 12, we arrive at two   12. Accumulative rate of entropy generation per unit volume as function of the radius of both geometries. The left plot depicts the hub (mass averaging between 0% and 35% span), mid-span (35%e65% span), and shroud (65%e100% span) sections of the 3D simulations. The right plot illustrates the Q3D and 3D simulations; the latter considers the entropy generation across the whole span and only the mid-span section (35% to 70% span). Table 2 Time-averaged performance comparison between the old and new blade for the high-expansion ORC turbine for a Q3D and 3D simulation. We measure efficiency in terms of isentropic efficiency h is ¼ ðh in Àhout Þ=ðh in Àh out;is Þ and the total-to-static efficiency h ts ¼ ðh in;0 Àh out;0 Þ=ðh in;0 Àh out;is Þ with h as enthalpy. Moreover, p S=R is the static pressure at the stator/rotor interface while p nozzleout is the nozzle's design static outlet pressure. We averaged in time (for one period) and in space using area averaged for the pressure and mass averaged for the rest. conclusions. (1) The strong bow shock of the old blade disturbs the flow field upstream, as depicted in Fig. 6(a), resulting in a significant increase of the entropy generation between the throat and STE as compared to the new blade; this demonstrates the importance of taking into account unsteady effects, as steady-state simulations can not capture this phenomenon. (2) Reducing the flow detachment lowers the entropy generation in the rotor passage as the slope of the accumulative entropy generation line between the RLE and RTE has decreased significantly for the new blade if compared to the old one.

Q3D/3D stator
In terms of the stator nozzle, there is no significant qualitative or quantitative difference between the Q3D and 3D simulations. The flow field and the accumulative entropy generation from the Q3D simulations agree well with the ones from the 3D simulations for the stator nozzle. Therefore, we conclude that the nozzle can be designed using a mid-span approach because there are no additional 3D effects (apart from the boundary layers at hub and shroud) that affect the stator performance. Furthermore, entropy generation before the throat is negligible due to the low velocity.

Q3D/3D rotor
For the rotor, the comparison between the Q3D and 3D simulations prove that designing a rotor blade using a mid-span approach can potentially conceal additional losses. Both of the rotor designs have substantial entropy generation due to secondary span-wise flow effects that cannot be captured by Q3D simulations. The 3D simulations exhibit two different span-wise mechanisms that can introduce additional losses: 1. Span-wise effects close to the rotor inlet can significantly alter the PR over the free expansion region between the stator and rotor. In the case of the new blade, the sharp expansion at the rotor inlet is modifying the pressure between the stator and rotor (see Table 2), leading to an increase in entropy generation in that region; this is not captured by the Q3D simulation (see Q3D concerning the 3D at mid-span for the new blade in Fig. 12).
Concerning the old blade, however, a 3D effect (the secondary flow) occurs further downstream of the rotor passage because of the more conservative increase of the shroud height. It, therefore, does not affect the free expansion region between the stator and the rotor. 2. 3D flow effects can increase the entropy generation locally. For the old blade design, the accumulation of entropy generation is equivalent for both simulation types (Q3D and 3D) across the turbine stage until the rotor passage, as depicted by the right plot of Fig. 12. However, in the rotor passage of the old blade, the 3D simulation has a higher accumulation of entropy generation where the secondary flow is located. Fig. 12 shows that the turbine with the new blade has a higher entropy generation, than the one with the old blade, at the exit of the stage. This can be explained using the mechanisms described in section 5.1.3. First, due to the first span-wise mechanism, the PR over the free expansion region increases for the new rotor blade. Extra entropy is generated near the hub and shroud compared to the old blade (see left Fig. 12), leading to a larger entropy accumulation at the inlet of the rotor. Second, due to the second spanwise mechanism, there is a steep increase in entropy production in the rotor passage near the shroud of the old blade, a loss mechanism that is considerably less for the new rotor. Furthermore, because the entropy production at mid-span is almost similar for both blade geometries, we can conclude that the reduction of entropy generation near the shroud of the new geometry is the leading cause for a higher entropy generation in the rotor passage compared to the old blade.

3D simulations
The entropy generation analysis demonstrates the importance of considering the span-wise and unsteady flow effects, and in particular, its influence over the pressure between the rotor and stator (p S=R ). The unsteady stator-rotor interaction strongly modifies the free expansion region, e.g., there is a quantitative difference in entropy generation between both rotor blades for the Q3D simulations. In the 3D simulation of the new blade, p S=R changes significantly from its Q3D counterpart (z0:16Dbar) due to the sharp expansion at the shroud of the RLE. This pressure is of paramount importance because it affects the flow field in the free expansion region (oblique shock waves and viscous wake from the STE) and, therefore, the uniformity of the flow entering the rotor blade passage. Moreover, the region between the stator and rotor can be the source of additional (mixing) losses that are not accounted for in a steady-state Q3D simulation, as seen in the right plot of Fig. 12. Even though the rotor performance does not directly influence the stator nozzle design, special care must be taken with the pressure between the stator and rotor as it affects the oblique shock wave from the STE and the flow uniformity at the rotor inlet.

Turbine performance
The 3D calculations give a worse performance of the highexpansion ORC turbines compared to the Q3D counterpart. We have measured the performance of the ORC turbine employing the power output, and the isentropic and total-to-static efficiency, as displayed by Table 2. On the one hand, for the old design, there is a 3:15D% and 1:19D% drop in the 3D simulation relative to isentropic and total-to-static efficiency of the Q3D simulation, respectively, resulting in less power output: 15.3 DkW. On the other hand, the new turbine stage also has a drop in performance for the 3D simulation: 3.57D%, À3.30D%, and À20.9 DkW for the total-tostatic, isentropic efficiency and power output, respectively. The additional loss mechanisms, quantifiable by a 3D simulation, decrease the performance of the high-expansion cantilever ORC turbine stage, proving the importance of the span-wise effects in such a machine.
The steady-state optimization that produced the new rotor blade improves the performance of the cantilever ORC turbine. There is an increment on the power output d if compared to the old blade d for both types of simulations: 18.38 DkW and 11.97 DkW for Q3D and 3D, respectively. In terms of total-to-static efficiency, the improvement of the new rotor blade is of þ6.81D% and þ4.22D% for Q3D and 3D, respectively. The ORC turbine does not have such an improvement in terms of isentropic efficiency ( < þ0.55D%), a result that matches the entropy generation analysis of the 3D simulation for both geometries (right plot of Fig. 12). The reason for such an improvement of power output and total-tostatic efficiency of the high-expansion ORC turbine is because the new blade is capable of extracting more kinetic energy from the flow, as illustrated by Table 2.

Conclusion and future works
This article presents detailed unsteady numerical simulations via three-dimensional calculations of high-expansion radial inflow ORC turbines, which operates in the dense-vapor region. To account for the unsteady stator-rotor interaction, a conservative flux assembling technique for non-matching three-dimensional meshes is applied. We consider two rotor blade shapes, a geometry that closely resembles the turboexpander of an ORC manufacturer which we call old blade, and a new blade derived from an indirect design methodology and an aerodynamic optimization. The simulations indicate strong unsteady effects, especially in the rotor blade passage, for both geometries. Because of the highly supersonic flow at the stator exit, unsteady shock waves emanate from the trailing edge of the stator and interact downstream with a bow shock at the rotor leading edge and the viscous structures at the suction side of the blade. In terms of three-dimensional effects, the stator blade can be designed at mid-span without any consequence in the third-direction; this area of the ORC turbine has no threedimensional effects apart from the boundary layer at the stator's hub and shroud. However, both rotor blades suffer from strong three-dimensional effects, e.g., a secondary flow inside the rotor passage and a sharp expansion in the shroud region at the inlet of the rotor. These span-wise mechanisms introduce additional irreversible losses: in the free expansion region between the stator and rotor, and (locally) in the rotor passage. The high-expansion cantilever turbine designer needs to consider these loss mechanisms as they decrease the overall stage performance (À3D% in isentropic efficiency for our particular case). In a broader sense, our detailed simulation results are valuable for designers in (1) ORC applications with different working fluids under similar thermodynamic conditions (reduced pressure and temperature) and (2) supersonic turbomachinery where there is a strong unsteady shock-wave interaction between the stator and the rotor.
Considering each geometry in more detail, the flow field of the old blade geometry indicates a large recirculation bubble at the suction side of the blade and strong three-dimensional effects; a secondary flow is generated at the rotor leading edge close to the shroud of the blade. Moreover, the blade loading of this geometry features sharp fluctuations d as high as 1 bar d caused by the unsteady shock wave interactions. On the other hand, the new blade has a more smooth pressure distribution, decreases the kinetic energy wasted at the outflow, and reduces the secondary flow in the rotor passage. As a consequence, an increase of the power output (þ10 DkW) is predicted by the simulation of the new highexpansion ORC turbine if compared to the old one. Still, the new blade suffers from another span-wise effect (substantial expansion at the shroud near the rotor leading edge), which modifies the pressure between the stator and rotor. Even though the stator nozzle design is independent of the unsteady and threedimensional effects, the free expansion area between the stator and rotor is being modified by the rotor blade's profile and height distribution. These findings can be used for a stator redesign by changing the nozzle's back-pressure.
In a future publication, we will present a comparison between our simulations and experimental measurements. The ORC manufacturer is currently instrumenting both the old and new highexpansion turbine designs, which will give us the unique opportunity to validate our detailed three-dimensional unsteady simulations and to quantify the impact of the reduced blade count approximation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.