Design and Parametric Analysis of a Supersonic Turbine for Rotating Detonation Engine Applications

: Pressure gain combustion is a promising alternative to conventional gas turbine technologies and within this class the Rotating Detonation Engine has the greatest potential. The Fickett–Jacobs cycle can theoretically increase the efﬁciency by 15% for medium pressure ratios, but the combustion chamber delivers a strongly non-uniform ﬂow; in these conditions, conventionally designed turbines are inadequate with an efﬁciency below 30%. In this paper, an original mean-line code was developed to perform an advanced preliminary design of a supersonic turbine; self-starting capability of the supersonic channel has been veriﬁed through Kantrowitz and Donaldson theory; the design of the supersonic proﬁle was carried out employing the Method of Characteristics; an accurate evaluation of the aerodynamic losses has been achieved by considering shock waves, proﬁle, and mixing losses. Afterwards, an automated Computational Fluid Dynamics (CFD) based optimization process was developed to ﬁnd the optimal loading condition that minimizes losses while delivering a sufﬁciently uniform ﬂow at outlet. Finally, a novel parametric analysis was performed considering the effect of inlet angle, Mach number, reaction degree, peripheral velocity, and blade height ratio on the turbine stage performance. This analysis has revealed for the ﬁrst time, in authors knowledge, that this type of machines can achieve efﬁciencies over 70%.


Introduction
Factors such as concerns about climate change, more sustainable mobility, and the target imposed upon the aviation sector in terms of pollutant emissions, have caused an incentive for ever more efficient engines.Given the great flexibility and responsiveness of gas turbines, they will continue to represent the core element of aeroengines for a long time, especially for medium/long haul.Their application will also benefit from the progressive shift towards non-fossil fuels, such as Sustainable Aviation Fuels or hydrogen [1].In addition to the engine core efficiency, weight is also an important task to be addressed by future technologies: compact and high specific work engines are expected to be preferred.These meet the need and potential of the energy sector, where gas turbines will play an important role in balancing the discontinuous supply by non-schedulable renewable energy sources, such as wind and solar power, providing an important prerequisite for the success of renewable technologies in the overall energy mix.In this regard, the adoption of zero-carbon combustion technologies, which are expected to be ever more common soon thanks to the use of carbon free fuels such as Hydrogen and Ammonia [2], can be accelerated.Therefore, it is necessary to invest in a renewal of gas turbine technologies that can fulfil the global need for higher efficiency and lower environmental impact both in the transport and in the energy conversion sectors.
Among the possible innovative gas turbine technologies, pressure gain combustion (PGC) has the greatest potential in terms of overall performance improvement and practical adoption.This technology utilizes a change in the thermodynamic cycle to achieve a significant efficiency boost for the same turbine inlet temperature.In fact, pressure gain combustion cycles exploit the more thermodynamically efficient processes of constant volume or detonation-based heat release (this latter well approximated by the Fickett-Jacobs cycle): the maximum pressure of these cycles rises during the combustion process with a further increase in compression work.As such, over 15% of thermodynamic efficiency can be potentially achieved respect to a Joule-Brayton cycle working at the same maximum temperature [3,4].
This potential gain justifies the rapidly growing interest in PGC, especially for its ready adaptability to hydrogen combustion and for the chance to realize efficient ultra-compact engines: among the different technologies the so-called rotating detonation engines are the most promising.In this kind of gas turbine, a detonation wave propagates around the perimeter of an annular combustion chamber to achieve the pressure increase.This supersonic combustion wave, with flame speeds around 2 km/s, allows for a lower number of compression stages, a more compact design with reduced combustor weight and a lower specific fuel consumption [5].
In addition to all issues related to the combustion chamber (out of this paper purpose), the transition duct connecting the combustion chamber to the turbine [6][7][8][9] and the turbine itself constitute important tasks to be addressed to have efficient cycles.
As for the turbine design, some papers are available in the open literature describing different issues.Paniagua et al. [10] report a preliminary design and a possible turbomachine stage arrangement, while Sousa et al. [11] describe guidelines for entropy minimization by working on the profile shape.Papers [12,13] report results and concepts related to the unsteadiness generated by the combustion chamber and the corresponding turbine frequency response.Liu et al. [14] present and discuss a CFD analysis of such turbines.
The size of the combustion chamber, the number of compression stages are expected to reduce with respect to actual ones, thus contributing to the overall engine weight reduction.This paper focuses on the turbine design too, because it constitutes a primary task for an efficient thermodynamic cycle: in fact, due to the high speed, high temperature and fluctuating flows released by the transition duct the turbine design is complex and requires different methodologies with respect to conventional gas turbine for Joule-Brayton cycles.
Among papers dealing with the turbine design, none of them proposes a parametric analysis on geometrical and fluid-dynamic parameters.The present paper proposes a novel parametric analysis on two geometrical parameters, namely solidity and blade height, on two flow features, such as inlet flow angle and inlet Mach number, and on the reaction degree being this latter one of the main issues describing the stage arrangement and affecting the stage efficiency.Authors are also aware about the criticism of the thermal load and cooling requirements on this kind of application even though they are not addressed in the present paper as they require a specific study.

Mean-Line Design
An original mean-line code was developed to perform a blade-to-blade design of a supersonic turbine.The code is based on the standard turbine theory and some specific features for supersonic blades are added: the Kantrowitz limit for a started blade row [15], an accurate evaluation and prediction of the shock losses and a Method of Characteristics to generate the supersonic profile shape.Before describing the logic applied to design the turbine, it is deemed necessary to briefly introduce the theoretical elements mentioned above.

Started or Unstarted Blade Row
A convergent passage, which is typical of turbines, behaves like a diffuser when subjected to a supersonic flow, thus reducing the velocity and increasing the pressure of the stream.In these conditions, it is critical to assess if the blade row is started or unstarted, where "started" refers to an oblique shock configuration at the cascade inlet which is less dissipative than the normal shock of the "unstarted" condition (Figure 1).Kantrowitz and Donaldson showed that the strength and the position of the shock depends on the geometry of the diffuser and on the upstream Mach number [15].Employing a conventional one-dimensional flow theory, the maximum permissible contraction ratio to have a self-started condition can be determined using Equation (1), where M n is the Mach number behind the normal shock wave.
Between the self-started and the unstarted region, there is a third intermediate one: in this region it is possible to achieve a started operating mode, but the diffuser must start from a self-started condition and then either the inlet Mach number should be reduced or a variable diffuser geometry should be employed (Figure 2).This latter concept is relevant in off design conditions.
As the blade turning influences the outlet area in a purely converging blade-to-blade channel, the self-starting limit affects also the maximum allowable turning as shown in Equation (2).At the same time, the blade turning impacts on the outlet tangential velocity magnitude and therefore on the work exchange.The blade height, contributing to the definition of the outlet discharge section, is also important for the amount of work extracted by the stage.

Evaluation and Prediction of Aerodynamic Losses
A profile working in a supersonic flow environment is less efficient than its subsonic counterpart due to the generation of a complex shock pattern in the channel.The capability to effectively estimate losses is a key requirement for a reliable mean-line code and this was achieved by taking into account the bow shock wave at the leading edge and the boundary layer (methodology also used in [16]) in addition to the losses related to the first oblique shock wave reflection and the mixing in the case of supersonic axial Mach number.
The detached bow shock wave generated from the leading edge is the main source of entropy production inside the turbine and it is responsible alone of almost two thirds of the total losses.The shape of the detached shock is predicted through the method developed and validated by Moeckel for two-dimensional and axially symmetric bodies [17].Once the detached shock wave has been determined, the entropy production is calculated locally by applying the oblique shock wave relations and then results are integrally averaged, as shown in Equation (3), where the subscript "up" and "l" refer, respectively, to the upper and lower blade shocks (Figure 1b).Notice that the intersection between shock waves is neglected as a simplifying assumption.
The detached shock wave losses are strongly dependent on the dimensionless ratio between the blade pitch and its thickness.This fact has physical meaning because a smaller thickness compared to the blade passage makes the leading-edge sharper and reduces the subsonic portion of the flow behind the bow shock, which is the main responsible of the losses.At the same time, a too high ratio cannot be used because the leading edge would be too small for structural and thermal integrity.
The first oblique shock wave, which is generated when the detached shock wave is reflected by the blade surface, was also determined and the associated losses were calculated to improve the accuracy in the calculation of the losses.Wave reflection is calculated by considering that the flow immediately adjacent to the wall must be parallel to it: by applying the classical theory of aerodynamics [18], the implicit Equation ( 4) is derived and solved numerically to find the slope of the reflected shock.The subsequent shock wave reflections were instead not considered because they are significantly weaker, and the assumptions required to evaluate them could seriously affect the reliability of the results.
The second most important source of losses in a supersonic blade row are boundary layer and mixing losses, having preliminarily neglected the shock-boundary layer interaction.The estimate is performed by a two-step procedure: First the boundary layer quantities for a compressible turbulent boundary layer are calculated using the procedure derived by Stratford and Beavers [19]; (ii) Then the overall loss coefficients are evaluated by applying the method described by Stewart for high velocity flows [20].
In the original methodology, Stewart considered flows with a subsonic axial Mach number; since in our application the axial Mach number is supersonic, the original approach was extended.Without repeating the whole demonstration, already reported in [20], here are reported only the analytical steps describing the novel extension to supersonic axial Mach numbers.
The solution of the second order equation gives two results and in the subsonic case the negative sign was selected to find the correct conditions after mixing.This is physically consistent because the increase in available area from the trailing edge to the mixed region reduces the velocity in subsonic conditions, but for a supersonic axial Mach number the channel behaves like a nozzle and not like a diffuser.Hence, differently from Stewart, the solution of the second order equation selected is the one with the positive sign Equation (6).
Once all the quantities have been determined, the total pressure ratio between the inlet and the completely mixed region at outlet is calculated through Equation (7).In literature there are mainly two methods based on MOC to design the profile of a supersonic blade row: corner-flow method by Shapiro [21] and vortex-flow method by Goldman [22].The Goldman one has been implemented in this mean-line code because it provides a smoother velocity distribution with the highest loading at mid-chord and less constraints on the flow turning.
A blade profile designed with the vortex-flow method is composed of the following fundamental parts, each of them having a specific function here briefly summarized for readers convenience (Figure 3): the inlet transition arcs convert a uniform parallel flow at the blade inlet into a vortex flow field, the concentric circular arcs preserve the vortex-flow structure while turning the flow with a constant surface Mach number and the outlet transition arcs reconvert the vortex flow into a uniform parallel flow.
The vortex flow method requires several inputs to define the blade profile: inlet Mach number M in , outlet Mach number M o , lower surface Mach number M l , upper surface Mach number M up , inlet flow angle, and the specific heat ratio of the working fluid.Although all the other inputs are given, M l and M up are selected to avoid issues related to flow separation: as suggested in [22], the maximum M up and the minimum M l are calculated to obtain a pressure recovery coefficient of 0.5.The actual M l and M up have been selected with a further 12% margin respect to those limit values.Once the dimensionless profile has been generated, two key geometrical parameters can be calculated: the dimensionless chord c * (c * = c/r) and pitch g * .In the mean-line code the chord is the only dimensional quantity assigned for the profile geometry, thus allowing for the calculation of r (the geometrical scale parameter).The geometries presented as examples in this work are for ground power production, hence the dimensionless ratios (c/g and g/thickness) were optimized to achieve the maximum efficiency.For aeronautical application, the turbine design requires a compromise between compactness and efficiency, which leads to different values of the dimensionless ratios and a shorter chord length.Once these quantities have been assigned, the dimensional profile, the pitch and the thickness can be finally calculated.
The profile generated by the Method of Characteristics is optimal from an aerodynamic point of view, but it is purely theoretical since zero thickness at leading and trailing edge is impossible to manufacture and has poor mechanical and thermodynamic performance, especially with a high velocity and high temperature flow.Hence, once an appropriate value of the thickness is selected according to mechanical calculations, the upper and lower surfaces are shifted, and the blade profile is closed by means of circular arcs.
A structural analysis was carried out proving that the blade can withstand without significant issues the static loadings; on the other hand, a detailed thermal analysis was not included in this work, because for this type of turbines this is an advanced topic and will be properly dealt in future studies.

Mean-Line Code: Description and Logic Explanation
The mean-line code implemented employs a multi-loop iterative approach to achieve more than one design objective, which is possible by changing quantities that would be usually fixed in a direct approach; while the stator blade height ratio H 2 /H 1 is an input selected by the user, the rotor blade height ratio H 4 /H 3 is iteratively adjusted to maintain the same axial velocity between rotor inlet and outlet, the latter being a common design feature of axial turbines.The absolute stage inlet flow angle is assumed equal to 0°: in fact, notwithstanding the flow entering the turbine is oscillating in pressure, temperature, and direction, its mass flow average is zero because fluctuations are symmetrical [8,23].By design the rotor absolute outlet flow angle is null to minimize the outlet Mach number and therefore to reduce the shock losses in the following stages.Finally, the code is implemented in a way that the user can select as an input the reaction degree X of the stage: this option does not only help the user to set a key parameter in the design of a turbine but has been extremely useful to perform the parametric analysis shown in Section 6.
Once stator/rotor losses have been calculated and the outlet total pressure is known, the outlet Mach number is found with the continuity equation and then all other quantities are determined by resorting to the ideal gas assumption.As soon as the stator and rotor design have been completed, the mean-line code checks if the Kantrowitz limit for a started blade row is satisfied.Then, it also verifies that the axial Mach number is higher than 1 to avoid the unique incidence problem [24], which is an unfavorable working condition considering the unsteadiness of the incoming flow.In conclusion the flowchart in Figure 4 shows the order and the logic of all the operations performed to design the rotor.

Computational Flow Model
To check the mean-line code validity and to deepen the understanding of the flow features inside the channel, CFD was run on a selected number of cases.Steady-state blade to blade and steady-state 3D simulations were performed on the commercial software Ansys CFX (19.1, Ansys, Canonsburg, U.S.), which implements a pressure-based implicit coupled solver with co-located grid arrangement.A high-resolution scheme was selected for the advection terms and the turbulence; this scheme belongs to the Total Variation Diminishing (TVD) algorithms, which have proven to be accurate in modelling shock wave reflections [25].The k−ω SST model was selected for turbulence closure, while the y + was maintained below 1 for proper boundary layer resolution.The working fluid is air which is assumed to be a thermally and calorically perfect gas.
Since the inlet is supersonic, all the relevant quantities are assigned at this boundary.The outlet is supersonic and the lateral boundaries are periodic.Adiabatic and no-slip conditions were selected for the blade walls, while adiabatic and free-slip conditions were set for the hub and shroud walls; this choice was made also for the 3D simulations, because at this stage the aim was to understand the effect of the radial variation of the quantities, working with a mesh that is suitable to investigate a large design space.Finally, a mixing plane interface was employed in the stator-rotor coupled simulations and all the computations were considered at convergence when the residuals are in the order of 10 −6 .

CFD Code Validation
Since the flow is supersonic with shocks and expansion waves in the domain, it is pivotal for the results reliability to assess if ANSYS CFX can accurately capture these distinct features with the settings reported above.Two test cases were selected to perform the code validation: the Sod shock tube [26] and the compression corner [27][28][29].Notice that a grid sensitivity analysis by way of consecutive refinements was performed for both test cases and all the reported results are grid independent.
In the shock tube problem, a diaphragm is initially separating a high-pressure gas on the left (x < 0 m) from a low-pressure gas on the right (x > 0 m).When the diaphragm is broken (t = 0 s), three new features appear in the domain: a shock wave and a contact surface propagating to the right and an expansion wave propagating to the left.The hypothesis for the analytical solution are that the problem is 1D, the tube is infinitely long and the flow is inviscid [18,26].Since ANSYS CFX works only with 3D domains, the geometry is modelled as a 10 m long tube, while the other two dimensions are significantly smaller and meshed only with one cell.Figure 5 shows a comparison between the CFD solution and the analytical one at t = 0.001 s: the two curves are almost identical for all the quantities proving that Ansys CFX is perfectly able to predict strength, position, and the wave speed of all three features.The 2D compression corner test case by Smits et al. (Figure 6) was selected because the shock/boundary-layer interaction is an advanced problem that contains altogether the difficulties of turbulence, compressibility, and viscous-inviscid interaction phenomena [28,29].Since the inlet Mach number is supersonic (M = 2.85), an oblique shock wave is generated when the wall turns into the flow domain; depending on the inlet and geometry condition, the interaction between the shock and the boundary layer can lead or not to a boundary layer separation.The computational domain is 3D, but one dimension is significantly smaller than the others and meshed only with one cell.The k−ω SST model was employed as turbulence model, while the y + was maintained below 1 for proper boundary layer resolution.The validation was performed for two ramp angles, respectively, 8 and 16 degrees.Figure 7 plots the surface pressure distribution for both cases and compares them to the experimental data of Smits et al. [28].For the 8-degree case, ANSYS CFX captures accurately the starting point of pressure increase, the maximum pressure after the shock and the overall trend.For the 16 degree case, CFD is correctly predicting the starting point of pressure increase and the maximum pressure, but there is a slight difference in the trend: in the CFD results there is a weak separation at the corner, which is not seen in the experimental data, but as stated by the authors the 16 degree case is at incipient separation and this aspect may be hidden by the uncertainty of the measurements [29].Finally, ANSYS CFX has been extensively used by our turbomachinery group and lots of effort has been dedicated over the years in properly validating the CFD results with the experimental data acquired on the test facilities present at the Laboratory of Fluid Machines (LFM) [30,31].

Mesh Independence Analysis
ANSYS TurboGrid was used to generate a structured mesh with hexahedral cells, high orthogonality, and low skewness.The mesh independence analysis was performed in two steps: first, the number of cells required in each layer for a 2D blade-to-blade grid independent solution was determined; then, the number of layers was increased in the radial direction till the independence condition was also satisfied for the 3D solution.
The grid convergence index method [32] was applied to perform the mesh sensitivity analysis and the standard deviation of the pressure distribution was selected as the reference quantity.For the blade-to-blade simulation (2 cells in radial direction), 100 k, 200 k, 400 k, 600 k, and 1600 k meshes were tested; the 600 k mesh can be considered as grid independent with a Grid Convergence Index (GCI) of 0.67%, calculated with the pressure distribution extracted at the outlet section mid-span.Then, 10, 30, and 50 layers were prescribed in the radial direction and the pressure distribution was tangentially averaged and extracted in the spanwise direction at half chord distance from the trailing edge; the 9 million cells mesh (30 layers) achieved the independence condition for both stator and rotor, with a GCI of 0.28% and 0.05%, respectively.

Mean-Line Code Results Validation through CFD in Design Condition
The physical parameters reported in Table 1 have been selected considering typical quantities for RDEs and high-pressure turbines [5,8,10,14,23,[33][34][35][36], while the geometrical parameters are chosen considering the blade structural feasibility.The validation procedure is accomplished by comparing the code outputs with the mass flow averages extracted from CFD simulations performed on the stator blade.The total entropy predicted by the mean line code is 44.17J/(kg K), while the value obtained from CFD is 43.87 J/(kg K), which results in a 0.7% deviation; pressure and temperature deviation are also within 1% at stator outlet.This result is very promising and assess the present mean line code as a primary tool in the preliminary turbine design.

Research of Optimal Loading by an Automated Optimization Process
Several studies have described optimal loading conditions for subsonic axial turbines, meanwhile there is a major knowledge gap regarding supersonic turbines.An original completely automated optimization process, which makes use of Python, MATLAB, and Ansys packages, was developed to find the optimal loading value that minimizes the losses while delivering a sufficiently uniform flow at outlet.
Zweifel's tangential force coefficient Equation ( 8) was selected as a driver parameter for the whole optimization process because of its capability to quantify the forces acting on a blade: its meaning is not related only to the blade solidity, but it also considers the flow deflection and the variation in the axial component of the velocity [37].
The search interval for the optimal loading condition was established after a series of preliminary CFD simulations.The upper limit was selected starting from the need to have sufficiently uniform flow at stator outlet to allow an efficient rotor operation.Blades cannot be too far apart because the outlet flow would be less uniform in magnitude and direction (Figure 9a).The lower limit is the solidity proposed by the MOC (Figure 9b): the associated C f t is quite low because the main objective of the vortex-flow method is to obtain a completely uniform flow at exit using a high number of blades.Then, it is necessary to consider that the original MOC designed a profile with infinitesimal thickness at leading and trailing edge, so the channel and behavior of the blade were modified when a finite thickness was assigned and, therefore, a new optimal value of the tangential force coefficient is expected.The optimal design should minimize losses while matching uniformity requirements.The selected optimization algorithm is the "golden section search method", which has a simple robust logic and converges towards a local minimum.This latter feature is a weak point of this method, but preliminary simulations confirmed that within the selected interval the trend of the entropy generated has only one minimum.
Effectively, for each C f t guessed by the algorithm, a high-fidelity evaluation process is performed: (i) the mean-line code (integrated with the MOC) generates the blade profile; (ii) blade-to-blade CFD simulations are set-up and run; (iii) data are extracted from CFD code; (iv) by post-processing, the entropy generation is calculated.
The entire process (described in Figure 10), including the optimization algorithm, has been fully automatized to significantly reduce the human effort and to easily repeat the whole operation with different inlet parameters.Notice that it is fundamental to run the mean-line code each time, because the search of the optimal value of C f t is performed while maintaining all the other dimensionless numbers fixed.Meshing the domain and the steps regarding the simulation are performed twice for each case: the first mesh is coarse and it is used to quickly initialize the whole domain and then the fine mesh (600k cells) is employed to have high quality data.The results of the optimization procedure applied to the stator blade row are shown in Figure 11; the smallest value of C f t is associated to the value of the solidity given by the MOC and the entropy generated in that condition is almost 18% higher than the optimal value; this improvement clearly justifies the novel optimization process developed for this purpose.Figure 11 also reports the mean-line code losses prediction, and they are very promising with a maximum error with respect to CFD lower than 8%, that decreases to less than 1% in the optimal working condition.For low C f t (higher number of blades and smaller pitch), the error increases because a higher number of oblique shock wave reflections are generated in the passage, all neglected but the first one by the mean-line code.At high C f t (lower number of blades), the increase in losses is explained by the small separation bubble formed behind the reflection point on the blade suction side (Figure 9a).It is also very interesting to observe that the recommended value of C f t in literature for conventional machines is 0.8, but the optimal values obtained here are quite different: this suggests that the use of correlations developed for conventional machines should be carefully assessed prior to their application on supersonic turbines.
Figure 12 shows the stator blade in optimal loading condition, which reduces both the wetted surface (and hence profile losses) and the shock wave losses.At the leading edge there is a detached bow shock wave, which is reflected only two times in the channel (labeled as A and B in Figure 12); the reflection in point B is the most critical because it generates a separation bubble, whose size and therefore loss content mainly depend on the blade overspeed and on the reflected shock strength.With a more loaded blade, a larger amount of flow is involved in separation and this leads to a notable increase in losses (Figure 9a); in optimal loading condition the separation bubble is considerably smaller and the boundary layer thickening in the interaction zone is limited.

Stator Off-design Analysis
Once the methodology has been assessed and the optimal solidity found, stator offdesign performances can be investigated.The two fluid dynamic parameters chosen to be varied are the inlet flow angle and the inlet Mach number.The reason for this choice is twofold: (a) the need to understand the behavior of the supersonic turbine in a wide range of different inlet conditions, (b) the fact that rotating detonation combustion chambers release high fluctuating flows, mainly on these two parameters.

Inlet Angle Analysis
In this set of off-design blade-to-blade simulations, all inlet boundary conditions except the inlet angle have been fixed.The angle has been varied between +/−16°with a step of 4°.Since the absolute Mach number is maintained constant, a secondary effect of the change in direction is a reduction in the mass flow rate and, consequently, a possible turbine unstarting for extreme value of the inlet angle.
Figure 13 reports the averaged value of the entropy production and total pressure ratio where it can be noted that curves are not "symmetric" with respect to the design condition, because of the shock waves lattice and its interaction with the pressure/suction side boundary layers.Entropy production more than doubles in the selected range thus evidencing how losses in the stator row are crucial for the stage performances.
As for the outlet flow direction, the stator shows good capability to deviate the flow close to the geometrical direction of 22°: Figure 13 also reports the outlet flow angle and its standard deviation calculated at the outlet section mid-span (reported on the figure as error bars).The flow-field is characterized by a medium level of variability at the outlet surface.The design condition shows the lowest deviation angle and a limited standard deviation.Meanwhile, in off-design condition the profile progressively loses the capability to impose the flow direction and the discrepancies with respect to the nominal condition are close to 4°in the worst scenario.These results allow to introduce some further considerations: since supersonic profiles do not perform well under negative incidence, in the frame of possible off-design conditions both an optimal robust design methodology and a slight positive incidence at the design point can be applied.By way of example, Mach number contours for the two cases of +/−12°are reported in Figure 14.In the case of −12°(Figure 14a), a separated region is found on the pressure side, due to the sharpness of leading edge.For the positive incidence case, (+12°Figure 14b), the shock pattern causes a flow detachment in the front part of the suction side that is subsequently partially recovered up to the last part of the suction side where a reflected shock causes the final flow detachment.The case of i = +16°, not shown for sake of brevity reports an unstarted situation on the stator with a strong normal shock upstream of the leading edge.

Inlet Mach Number Analysis
A second set of blade-to-blade simulations have been performed with variable inlet Mach number (range 1.6 to 2.8) but constant axial flow and static quantities (equal to the reference case).At first it has to be underlined that the change in the inlet Mach number, given all other quantities, modifies the mass flow rate and then, the following features can be noted (Figure 15): The outlet flow angle is reduced by few degrees when the speeds are lower than the design condition, while there is a low over-speed sensitivity; (b) The condition of uniformity of the outgoing flow is reduced as the speed increases; (c) Losses have a fairly linear trend around the design condition (M = 2): at higher Mach, the trend is more than linear since the reflected oblique shock affects the trailing edge region, favouring a premature flow separation.The analysis reported above for different stage inlet flow direction and magnitude demonstrate the complexity of the off-design operation of such turbines and the need for a proficient and robust design methodology.Preliminary evaluations of the stator outlet angle deviation, which in turn modifies the rotor inlet angle, are of paramount importance in the off-design stage analysis.

Effect of Reaction Degree and Peripheral Speed
In the context of improving the understanding and the applicability of such stages in ever more efficient gas turbines, a novel parametric analysis has been performed accounting for the effect of the reaction degree (X) and peripheral velocity (U) on the stage performance; the parametrical analysis is based on X and U as these two parameters hold a crucial position in the turbomachinery preliminary design.Such analysis made use of the mean line code and different blade height ratios have been also explored.Then, few selected cases have been also verified by means of CFD (Section 6.3).
The grid selected for these studies is defined in the following way: 12 equally spaced points for U from 260 m/s to 500 m/s and 20 equally spaced points for X from 0 to 0.55: overall, the grid consists in 240 points.Each case was studied in the condition of optimal solidity, this latter searched by the mean-line code thanks to its good agreement with CFD.
An in-depth description of the outcomes of the parametric analysis performed at constant stator blade height will be presented and then a comparison will be drawn at increasing values of H 2 /H 1 .

Constant Stator Blade Height
The link among the reaction degree, the peripheral velocity and the stator outlet flow angle is presented in Equation ( 9) to aid the results understanding, having imposed for the rotor constant axial velocity and axial outlet.
A different reaction imposes a different blade turning (low reaction requires high blade turning), that coupled to the peripheral speed allows to quantify the rotor inlet relative Mach number.The stator blade turning affects the area ratio and for high α 2 a significant deceleration on the stator may occur so that the Kantrowitz limit may be not satisfied leading to a normal shock upstream of the stator, too dissipative to be acceptable: the limit is reported on figures as red dashed lines.
Stator shock losses are primarily influenced by the primary shock and its first reflection which are not so affected by the reaction as, on the contrary, the boundary layer and mixing, even though these latter are less important.
An important source of loss that is strongly affected by the reaction degree is the one related to rotor shocks: this contribution is mainly driven by the relative rotor inlet Mach number.As shown in Figure 16a, the rotor inlet relative Mach number decreases with a reduction in the reaction or an increase in the peripheral speed (higher α 2 that leads to a smaller area at stator outlet, i.e., an increase in deceleration).Hence, rotor losses are strongly modulated by the variables U, X (Figure 16b).The boundary layer and mixing losses, as in conventional turbines, are mainly function of the outlet flow angle and the velocity.In the stator, losses increase closer to the starting limit because the flow angles are higher (Figure 17a); in the rotor, they decrease because, while the relative outlet flow does not vary significantly, the velocities are lower reducing especially the losses associated to mixing (Figure 17b).Considering that the increase in these losses in the stator is counterbalanced by the decrease in the same losses in the rotor, there is a significant cut in total losses near the starting limit because of a net reduction in shock losses for the rotor.
The Euler work obviously increases for lower X and higher U: Figure 18 describes this feature for different blade height ratios, here commented only for H 2 /H 1 = 1.The total-to-total efficiency (Figure 19a) is improving towards the Kantrowitz limit thanks to two effects: The amount of work extracted is increasing because the flow angles or U are larger; (ii) The total entropy production decreases.This is undoubtedly an unconventional result: for a subsonic turbine, the designer usually can either aim to maximize the efficiency or the amount of work extracted; here for a supersonic turbine no compromise between the two is required because the designer can theoretically create targets to extract the maximum work at the maximum efficiency.In conclusion, the X and U values should be carefully selected because the decrease in efficiency is strong, passing from 70% to 20% when lowering U and increasing X.

Variable Stator Blade Height
The blades with a variable height are designed with a constant hub radius, while the shroud radius is linearly varied to obtain the desired inlet to outlet blade height ratio.A variable stator blade height is a valuable option to extend the available design space limited by the starting condition: the reduction in area caused by larger flow angles is balanced by a taller blade at outlet section.In this way more pairs (U, X) are feasible, being characterized by larger flow deflections and therefore by a considerable increase in the amount of work extracted (Figure 18).Higher values of the Euler work carry two appealing advantages: (i) A lower total temperature at stage outlet reduces the cooling requirements on the next stages; (ii) A sizeable cut in NOx emissions, which are strongly dependent on high temperatures [38].To this end, the high mean flow velocity allows for a low permanence time at high temperature, thus preventing for high NOx production.
Notwithstanding the increase in the work exchange, M 3W is almost constant and low (even still supersonic) close to the starting limit even with an increasing stator blade height ratio.This result depends on the fact that the Kantrowitz limit is only function of the outlet-inlet area ratio and of inlet Mach number: the minimum area ratio is always the same and therefore the maximum possible deceleration will lead to the same M 3W close to the Kantrowitz line, which is shifting on the U-X diagram.From this, it is possible to conclude that in the most interesting design area, the shock wave losses generated by stator and rotor remain constant even when the stator blade height and the work are increasing.
The boundary layer and mixing losses are increasing with higher H 2 /H 1 ratio because the flow angles are larger for both the stator and the rotor.Hence, the stage total losses are increasing because the shock wave losses are almost constant, while the entropy generated by the boundary and mixing phenomena becomes higher.
Overall, the total-to-total efficiency is improving close to the Kantrowitz limit because the considerable increase in the turbine work overcomes the slight increase in losses (Figure 19b).By way of example, at constant blade height the maximum efficiency is 69% with 271 kJ/kg of work, while at H 2 /H 1 = 1.2 the maximum efficiency is 76% with 391 kJ/kg.Obviously, the design point cannot be at the edge of the Kantrowitz limit and a reasonable safety margin should be decided.To conclude, an increasing variable stator blade height is a beneficial option leading to an extension of the design space for U and X, thanks to a considerable increase in work extraction and higher efficiencies.

Detailed Study of Flow Three-Dimensionality for Multiple Design Points
Interesting results were revealed by the reaction degree and peripheral speed parametric analysis, but these data are generated by a mean-line code, which has been compared only to blade-to-blade simulations.Because one of the shortcomings of a mean-line code is the inability to consider the spanwise variation of the quantities, 3D coupled stator-rotor simulations with a mixing plane interface were performed in 4 meaningful design points (Figure 18 and Table 2).First, a detailed analysis will be provided for case D, which in authors opinion is the most interesting from a performance and flow features standpoint; then, the other points will be examined from a quantitative point of view.Design point D is directly at the edge of the Kantrowitz limit, but, as previously discussed (Section 6.2), in this region the turbine can extract a large amount of work with the highest efficiencies (Figure 19).From the results of the coupled 3D simulation, the turbine is extracting 36.02MW with an efficiency of 73.58%, which is in close agreement with the mean-line code prediction of 36.45MW with an efficiency of 74.47%.
Figure 20 displays the pressure gradient contour on the stator and rotor surfaces, and it is immediately clear that the shock structure has a complex 3D pattern.The first shock reflection line on the stator (1) is closer to the leading edge at the hub than at the shroud and this is the consequence of two factors: (a) the solidity at the hub is higher than the solidity at the shroud; (b) the hub wall is at constant radius, while the shroud wall is rising to achieve the desired blade height.The second and especially the third reflection lines on the stator suction side (2-3) exhibit a double shock structure, which is typical of a separation bubble: the upstream shock is the separation one, while the downstream shock is the reattachment one [39].It is significant to observe that at the stator last reflection (3), the shock at the hub is closer to the trailing edge reversing the initial trend; the reason behind this shock shifting is that the flow is gradually gaining tangential velocity and observing the radial equilibrium theory, the pressure at the hub is lower and the Mach number is higher.In the rotor blade, the intersection of the bow shock waves generates a Mach phenomenon, recognized by its typical structure composed of two slip lines and a normal shock between the two triple points (Figure 21).The Mach phenomenon manifests itself with high values of shock deflection for which the conventional solution with two oblique shock waves originating from the intersection point is not possible [39]; this is the reason why the Mach phenomenon appears only for design points at low reaction degree (A and D) and not for the higher ones (B and C).Even though this flow feature is undesired because it increases the pressure losses due to the normal portion of the shock, the low reaction turbines have still a better efficiency thanks to a significantly higher work that overcomes these losses.Table 3 reports all the relevant data for an in-depth comparison between the meanline code results and the solution of the coupled 3D CFD simulations.Considering the complexity of the problem, the mean-line code losses predictions are satisfactory, as shown by the entropy production values for both stator and rotor.Pressure, temperature, and Mach number at stator outlet are also accurately estimated and this is beneficial for a reliable rotor design.
Presently the main concern is the notable overturning in point B and C downstream of the stator (higher α 2 in CFD with respect to MLC), which leads to a higher power extraction and therefore higher efficiency compared to the mean-line code.This overturning results from the flow deviation induced by the last shock wave reflection on the suction side and by the complex flow structure built of shocks and expansion fans at trailing edge.This issue is emphasized for point B and C because most of the flow goes through the last reflected shock in the semi-vaned region, while in point A and D this happens only for a small portion; notice that the point D has also 3 reflections in the stator channel compared to the other three cases which have only 2. In conclusion, current work is focused on finding a solution to this issue, even though predicting with accuracy the position of the last shock reflection in a mean-line code is a difficult task.

Conclusions
This paper presents a supersonic turbine design along with its preliminary optimization.The research field of interest is related to the development of rotating detonation engines, for which conventional first stage turbines are unsatisfactory as they are unable to work at an acceptable efficiency.
A novel mean line code was developed to design supersonic turbines; this code includes additional features to assess the starting limit of a supersonic passage, to generate the supersonic profile through an implicit version of the MOC and to accurately predict the losses.For this final point, not only leading-edge bow shock waves and profile boundary layer were considered as recommended in literature, but our code also takes into account the first shock reflection and mixing losses by extending Stewart's theory for supersonic axial Mach numbers.The in-house mean-line code demonstrated excellent predictive capabilities during the validation process with CFD simulations.
Since no loading criterium, empirical or semi-empirical, exists for these pioneer turbines, an optimization procedure was developed and implemented to find the optimal loading condition for both stator and rotor.The optimization procedure, based on the golden-section search algorithm, iteratively and automatically evaluates the performance function by running the mean-line code, creating an adequate mesh and running the CFD simulation for the blade-to-blade section.In optimal loading condition the entropy generated by the stator is almost 18% lower compared to the original solidity proposed by the MOC (once a finite thickness is imposed), clearly justifying the optimization process developed for this purpose.
An extensive blade-to-blade simulation campaign was performed to assess the turbine performance in off-design condition.From this analysis it emerges that supersonic profiles exhibit a low over-speed sensitivity, while an optimal robust design should take into account that a slight positive incidence can be applied at design point to allow for a wider operating range.
Finally, a parametric analysis was carried out to assess the influence of the variable stator blade height, the reaction degree and the peripheral speed on the turbine stage performance.A variable stator blade height design mitigates the Kantrowitz limit, allowing to choose lower reaction degree and higher peripheral velocity designs; this extended design space is also characterized by a higher total-to-total efficiency up to 76% and it is important to underline that, previous to our work, efficiencies this high for this type of machines were not documented in literature.A second important and novel result of this analysis is that, opposite to common subsonic turbine design, the lower the reaction the higher the efficiency: this unconventional result comes from a reduction in rotor shock losses backed by a higher work generation when the reaction degree is lowered.The 3D coupled stator-rotor simulations were performed in 4 meaningful points to confirm these interesting outcomes; these simulations also display the full complexity of the flow field: the shocks exhibit a complex 3D structure from hub to shroud and Mach phenomena appear in the rotor channel at low reaction degree.

Figure 1 .
Figure 1.Comparison between unstarted (a) and started (b) condition for a supersonic blade row.

Figure 2 .
Figure 2. Contraction ratio for a started or unstarted blade row at varying inlet Mach number. p

Figure 3 .
Figure 3. Blade profile sector division in vortex-flow method.

Figure 4 .
Figure 4. Flowchart to illustrate the iterative process employed to design the rotor in the meanline code.

Figure 5 .
Figure 5.Comparison between the analytical and Ansys CFX solution at t = 0.001 s for the Sod shock tube problem.

Figure 6 .
Figure 6.Pressure field for the compression corner problem with a 16 degree ramp angle.

Figure 7 .
Figure 7.Comparison between the experimental solution [28] and Ansys CFX for the compression corner problem at 8°and 16°ramp angles.

Figure 8
Figure8compares the analytically determined shocks with a CFD simulation performed on the same supersonic profile and in the same flow condition.The density gradient (numerical Schlieren) is shown in the background, while in the foreground we can appreciate the detached shock wave and the first oblique shock wave reflection as calculated by the mean-line code.The predicted shock positions are really close to the ones obtained by CFD, which confirms the effectiveness of the Moeckel theory application.Even though the predicted reflected oblique shock wave is not exactly in the same location, slopes are very similar and thanks to this also the estimated entropy production.

Figure 8 .
Figure 8.Comparison between CFD numerical Schlieren in background and the analytically determined shocks in foreground (detached shock in yellow and reflected shock in purple in the online colored version).

Figure 9 .
Figure 9. (a) Minimum solidity conditions assumed in this paper.(b) Maximum solidity as proposed by MOC.

Figure 10 .
Figure 10.Flowchart to illustrate the procedure used to calculate the entropy production for each value of C f t given by the optimization algorithm.

Figure 11 .
Figure 11.Entropy produced by stator at variable C f t : in blue, values from CFD; in black, value from mean-line code results; in red, CFD-mean-line percentage difference.

Figure 12 .
Figure 12.Mach number field for the stator in optimal loading condition.

Figure 13 .
Figure 13.Entropy production and stator outlet angle with its standard deviation calculated at the outlet section mid-span for different inlet flow angles.

Figure 15 .
Figure 15.Stator entropy production and outlet flow angle with its standard deviation calculated at the outlet section mid-span for various inlet Mach numbers.

Figure 16 .
Figure 16.Relative Mach number at rotor inlet (a) and entropy production [J/(kg K)] due to shocks in the rotor (b) for constant stator blade height and different reaction degree and peripheral velocity.

Figure 17 .
Figure 17.Entropy production [J/(kg K)] in the stator (a) and in the rotor (b) due to boundary layer and mixing for constant stator blade height and different reaction degree and peripheral velocity.

Figure 18 .
Figure 18.Euler work [J/kg] performed by the stage for different stator blade height, different reaction degree and peripheral velocity.

Figure 19 .
Figure 19.Stage total to total efficiency for stator blade height ratios of 1 (a) and 1.2 (b) with different reaction degree and peripheral velocity.

Figure 20 .
Figure 20.The 3D shock structure visualization on stator suction side and rotor pressure side for design point D.

Figure 21 .
Figure 21.Mach phenomenon generated by the intersection of rotor bow-shock waves in design point D.

Table 1 .
Input parameters for the reference design stage.

Table 3 .
Detailed quantitative comparison between the results of the 3D coupled stator-rotor simulations and the mean-line code predictions.