Numerical Analysis of a Two-Phase Turbine: A Comparative Study Between Barotropic and Mixture Models

Unlocking the potential of low-temperature energy resources, such as geothermal energy, solar power, and waste heat from thermal engines and industrial processes, can significantly advance our efforts to reduce the emissions of the power generation sector. Over the years, significant attention has been dedicated to advancing low-temperature power generation technology, including developments in heat exchanger and turbomachinery designs as well as working fluid selection and cycle optimization [1]. Instead of limiting the focus to incremental improvements of existing technology, adopting innovative concepts enabled by two-phase turbines offers the possibility to achieve considerably higher performance gains [2]. In two-phase turbines, a pressurized liquid or liquid-dominated mixture undergoes a large expansion ratio, resulting in a substantial increase in gas volume fraction. A prime example leveraging this technology is the trilateral Rankine Cycle, which enables a better utilization of sensible heat sources due to a better match between the heat source and the working fluid temperature profiles in the primary heat exchanger [3]. Thermodynamic comparative analyses indicate that the Trilateral Rankine Cycle could achieve a higher power output than the conventional Rankine cycle, potentially increasing the net power by up to 30–40% depending on the temperature glide of the heat source [4,5]. However, the broad adoption of this concept has been limited by the absence of turbines that can efficiently and reliably expand a pressurized liquid into the two-phase region [6].

Besides their prospects for power generation, two-phase turbines also bring substantial advantages to refrigeration processes. In most refrigeration processes, cooling is achieved by expanding a liquid–vapor mixture in a Joule–Thomson valve. While this throttling process accomplishes the desired cooling effect, it also introduces a significant source of irreversibility into the system. Replacing the valve with a two-phase turbine offers two main benefits. First, the mixture's expansion in the turbine produces shaft work, which can either be converted to electrical power or used directly to drive the compressor [6]. Second, this expansion reduces the enthalpy of the mixture, resulting in an exit state with higher liquid content. This enhancement in liquid content increases the cooling capacity in refrigeration cycles [6], and directly boosts product yield in liquefaction processes, including liquefied natural gas production [7] and liquid air energy storage systems [8].

The aforementioned benefits were quantified by Hays and Brasz [6], who analyzed the effect of using two-phase turbines within vapor-compression refrigeration cycles. Their thermodynamic analysis predicted a reduction in compressor power ranging from 5% to 25%, depending on working conditions and the isentropic efficiency assumed for the two-phase turbine. In addition, the authors evaluated the benefits of using a two-phase turbine in the context of an industrial chiller with a cooling capacity of 500 tons (approximately 1760 kW). The conventional Joule–Thomson valve was replaced by a single-stage axial turbine, designed to deliver a nominal power output of about 15 kW. Performance tests indicated that, at the design point, the nozzle efficiency reached 84%, while the rotor efficiency was between 60% and 65%, resulting in an overall turbine efficiency of approximately 50–55%. Notably, the rotor efficiency was about 10% points lower than the value predicted during the design phase, highlighting a considerable gap between the model predictions and actual performance. Even if the turbine efficiency was lower than the design target, the integration of the two-phase expander was successful from the system's perspective. Specifically, replacing the Joule-Thomson value reduced the compression power and increased cooling capacity, leading to an overall improvement in the chiller system efficiency of 8%. Furthermore, the authors noted a substantial potential to improve the turbine efficiency through a better understanding of loss mechanisms and development of more accurate predictive tools.

Despite their prospects to improve the efficiency of power generation and refrigeration systems, the development of two-phase turbines has historically been considered technically challenging and economically unviable. This perception is largely due to the limitations of operating in the two-phase region, including the risk of erosion and cavitation, as well as the comparatively lower efficiency resulting from additional loss mechanisms that are absent in single-phase turbines.

Cavitation and erosion, which are commonly identified as risks in two-phase turbomachinery, may potentially be mitigated with an appropriate design of the turbine's operating conditions and blade geometries. Erosion, a process characterized by material wear resulting from liquid droplet impact on turbine blades, can be mitigated by keeping the fluid velocity below the material's threshold impact velocity [9]. Furthermore, erosion in two-phase turbines tends to be less severe than in condensing flows (e.g., last stages of steam turbines) because the liquid phase tends to form a flowing layer that protects the blades [10]. With regard to cavitation, a phenomenon occurring as vapor bubbles violently collapse where local pressure exceeds the fluid's saturation pressure, the risk is notably lower in two-phase turbines than in hydraulic pumps. This reduced risk is primarily due to the natural decrease in pressure along the expansion path, which inherently lessens the likelihood of cavitation. While empirical data on erosion and cavitation in two-phase turbines remain scarce, the study by Welch and Boyle offers a positive indication of the technical feasibility of two-phase turbines [11]. They tested a 7-kW impulse turbine for a trilateral Rankine cycle at a geothermal power plant and found no signs of erosion or cavitation in the rotor after 150 h of operation. While these initial findings are encouraging, additional experimental evidence is required to assess the risks under different operating conditions and over extended periods of time.

The comparatively low efficiency of two-phase turbines has been another factor limiting the broader adoption of this technology. In general, the efficiency of turbines is influenced by several inter-related flow mechanisms, commonly known as losses. These include viscous dissipation in blade boundary layers and wakes, tip-leakage and end-wall losses, as well as losses due to entropy generation in shock waves [12]. In the case of two-phase turbines, we encounter additional loss sources caused by the interaction between the liquid and gas phases. While the underlying flow physics in two-phase turbines are not yet fully understood, certain mechanisms have been identified as key contributors to increased losses. One such mechanism is an increased dissipation caused by drag forces arising from the relative velocity (also known as slip velocity) between liquid and gas phases [13]. Additionally, the different density of the phases causes flow segregation and formation of liquid films on the pressure surface of the blades, which is associated with higher trailing edge losses and lower turbine efficiency [14].

In addition to the specific loss mechanisms, the distinct thermophysical behavior of two-phase mixtures introduces further complexities into the design and operation of these turbines. The presence of gas bubbles in a liquid creates a mixture with unusual characteristics: it is easily compressed because of the presence of a gas phase, whereas it maintains a relatively high density because of the dominant mass of the liquid. One important consequence of this behavior is that the effective speed of sound of the mixture is significantly lower than that of each separate phase [15]. As a result, compressible flow phenomena, such as choking and shock waves, may occur at relatively modest velocities, making the design of efficient two-phase turbines more challenging.

In this context, computational fluid dynamics (CFD) stands as an essential tool for the design and analysis of two-phase turbines, enabling the simulation and visualization of complex flow configurations. This is crucial to predict the performance across different operation conditions and geometry configurations, enhance our understanding of the flow physics, and guide the design process toward turbines with higher efficiency. While CFD is an established tool for single-phase turbomachinery [16,17], its application for multiphase flow simulations is comparatively less mature. To date, the majority of multiphase CFD studies in turbomachinery have focused on the analysis of condensation in the final stages of steam turbines [18–20], and on cavitation in hydraulic pumps and turbines [21]. However, these application areas are not representative of the flow conditions occurring in two-phase turbines because condensation and cavitation are usually confined to specific regions of the flow, while two-phase turbines experience a substantial increase of volume fraction of gas along the expansion as well as significant two-phase interactions throughout the entire bulk of the flow.

To the best of our knowledge, the only study modeling a two-phase turbine using CFD is the work of Rane and He [22]. They utilized a pressure-based Euler–Euler solver for simulating the expansion of a geothermal brine within a radial outflow turbine characterized by a 100% degree of reaction. This unconventional turbine architecture has the advantage that it does not rely on impinging jets to drive the rotor, which precludes erosion problems typically associated with jets of geothermal brine. However, this type of turbine inherently suffers from low efficiency because of the adverse pressure gradient caused by centripetal acceleration and the high friction losses in the long, curved flow paths necessary to counterbalance lateral Coriolis forces. These factors lead to low isentropic efficiencies, typically ranging from 10% to 40%, making this turbine configuration unsuitable for achieving the high efficiency levels required for successful application in an industrial context. Simulation results were compared with experimental data. A deviation in mass flowrate within 8% and an overestimation of shaft power by 30–50% was reported. Notably, the model showed high sensitivity to tuning parameters, which significantly influenced the simulation outcomes.

In response to the identified research gaps, specifically the absence of an established CFD modeling methodology suitable for two-phase turbines, our study focuses on addressing this challenge. Our objective is to identify an accurate and robust CFD modeling approach for predicting the performance of two-phase turbines. To achieve this, we conduct a comparative analysis of two multiphase models: the mixture model and the barotropic model, each offering different levels of fidelity and complexity. Our evaluation is centered around the single-stage impulse turbine examined by Elliot [2], where we will compare our model predictions on nozzle and rotor performance with experimental data. This work represents, to the best of our knowledge, the first time that a two-phase turbine configuration suitable for industrial application is evaluated via three-dimensional (3D) CFD simulations and validated against experimental data.

The rest of the paper is organized as follows: First, we introduce the test case used for the simulations (Sec. 2), setting the context for the subsequent sections. After that (Sec. 3), the mixture and barotropic two-phase flow models are described. The simulation results are organized in two parts: an analysis of the nozzle using two-dimensional (2D) axisymmetric simulations is presented first (Sec. 4), followed by the 3D simulations of the full turbine stage (Sec. 5). The paper concludes (Sec. 6) with a summary of the key findings, highlighting several areas for further modeling improvements.

Most experimental studies on two-phase expansion primarily focus on nozzles, leaving a notable gap in experimental data concerning rotating machinery. Among the few experimental studies considering two-phase turbines, a majority either provide insufficiently detailed data for accurate replication via CFD simulations [6,11,23], or they examine unconventional turbine configurations yielding efficiencies unsuitable for industrial use [24,25]. Given this scarce landscape, the key subject of our analysis is the single-stage impulse turbine analyzed by Elliot [2]. This test case was selected because (1) it includes detailed performance data across different angular speeds and working fluid liquid-to-gas ratios, and (2) it is representative of typical flow conditions in two-phase turbines (i.e., large expansion ratio and significant increase of the volume fraction of the gas phase along the expansion).

As illustrated in Fig. 1, the Elliot turbine setup comprises a single converging–diverging nozzle followed by a rotor row with 151 blades. In this configuration, the working fluid, a mixture of liquid water and gaseous nitrogen at room temperature, undergoes an expansion from 2000 kPa down to atmospheric pressure within the nozzle. This expansion process generates a high-velocity jet that drives the rotor and produces torque at the shaft. The nozzle, specifically designed for two-phase applications, features a 10 deg converging angle and a 2.5 deg diverging angle. These are smoothly joined through a throat section with a width of 6.55 mm designed with a gradual area change to reduce velocity differences between the phases.

Our study encompasses two series of experimental tests based on Elliot's work. First, we analyze the performance of the uncut axisymmetric nozzle under various water-to-nitrogen mass ratios in Sec. 4. These initial simulations provide the foundation for the validation of the barotropic and mixture models in two dimensions. Subsequently, for the complete two-phase turbine testing, Elliot modified the nozzle exit to an elliptical shape by cutting it at a 20 deg angle relative to the axis. The simulation results for the complete turbine stage, conducted at a single liquid-to-gas ratio and different angular speeds, are presented in Sec. 5. This replication analysis offers insights into the level of accuracy achievable with the barotropic and mixture models in the context of 3D rotating machinery.

Selecting an appropriate modeling approach is essential to successfully predict the performance of two-phase turbines and achieve a satisfactory balance among accuracy, computational speed, and convergence robustness. Within this context, the spectrum of modeling approaches includes the two-fluid, mixture, and barotropic models, each offering distinct advantages and limitations. The two-fluid model is often considered more advanced, as it solves separate transport equations for each phase and provides a detailed representation of phase interactions [26]. However, this level of detail incurs a higher computational cost and potentially less robustness. Furthermore, the inherent uncertainties in the closure relations and parameters used in the two-fluid model can compromise its accuracy, rendering it less practical compared to simpler approaches [13]. By contrast, both the mixture and barotropic models solve a single set of transport equations for the two-phase mixture [27]. Although these methods may not capture phase interactions as effectively as the two-fluid model, they offer greater computational efficiency and robustness. This can be particularly beneficial in turbomachinery applications, where the primary focus is the assessment of overall performance, rather than the fluid dynamics of the phases and their interactions.

The mixture model solves the Reynolds-Averaged Navier–Stokes equations for the two-phase mixture, treating it as a single-phase fluid with averaged properties and incorporates an additional transport equation for the volume fraction of the dispersed phase. Additional details about the mathematical formulation of the mixture model can be found in Ref. [13].

Regarding the thermophysical modeling of the fluid properties, we make several assumptions that are valid for the temperature and pressure ranges of the cases considered in this work. First, the amount of nitrogen dissolved in the liquid and of water evaporated in the gas is negligible. Therefore, the two phases were modeled as pure fluids with no interphase mass transfer. In addition, due to the modest pressure levels, incompressible liquid and ideal gas models were used for the water and nitrogen, respectively. Moreover, the temperature along the expansion is nearly constant because the heat capacity of water is much higher than that of nitrogen. Consequently, we assume that the viscosity of both phases is constant.

The barotropic fluid model is a simplified mathematical representation of the fluid’s thermophysical behavior based on the assumption that all properties depend only on pressure [28]. This approach is well suited to describe the expansion process in turbines, where the evolution of fluid properties can be approximated as an isentropic process or, more generally, as a process characterized by a polytropic efficiency. Given that the fluid properties follow a polytropic process, solving an energy transport equation to determine the fluid’s thermodynamic state is unnecessary. Consequently, the flow problem only requires the solution of the mass and momentum equations for the mixture, along with a turbulence model. This makes the barotropic model inherently more robust and computationally efficient than the mixture model [27].

As a sample calculation, Fig. 2 depicts the density, gas volume fraction, temperature, and speed of sound of the mixture as a function of pressure for the different water-to-nitrogen ratios (⁠$R=y1/y2$⁠) examined in Elliot's experiments. As expected, the density of the mixture falls within that of the water and nitrogen phases. Unlike the expansion of an incompressible liquid, the density of the two-phase mixture experiences a significant variation with pressure because of the expansion of the nitrogen bubbles. This effect is also manifested as an increase of the void fraction (i.e., volume faction of gas) along the expansion.

Additionally, Fig. 2 illustrates the dependence of the speed of sound on both pressure and liquid-to-gas mass ratio. Notably, the speed of sound of the mixture is significantly lower than that of the individual liquid and gas phases. This distinctive behavior of two-phase mixtures is a direct consequence of the highly nonlinear character of Eq. (5), and it has important implications in the design of two-phase nozzles and turbines, as compressible flow phenomena like shock waves and choking can occurs at relatively modest velocities.

The calculations shown in Fig. 2 were performed under the assumption of a polytropic efficiency $ηp=100%$⁠. In order to explore the impact of varying polytropic efficiencies on fluid properties, we conducted a sensitivity analysis for the case with 68.00 water-to-nitrogen mass ratio. Figure 3 shows the density deviation as a function of pressure, this analysis indicates that the density predicted by the barotropic model is remarkably insensitive to changes in polytropic efficiency. Specifically, the relative difference in density between the extreme scenarios of isentropic (⁠$ηp=100%$⁠) and isenthalpic (⁠$ηp=0%$⁠) expansions is less than 0.5% across the entire range. This insensitivity can be attributed to the expansion process being nearly isothermal, as the liquid water has a much higher heat capacity than the gaseous nitrogen. Given the minimal impact of the polytropic efficiency, all barotropic model fluid properties were based on isentropic calculations.

Steady-state simulations of the convergent-divergent nozzle were executed using Ansys Fluent's 2D axisymmetric pressure-based solver [30]. The mixture model simulations were conducted using the built-in mixture model. With regards to the barotropic model, Fluent does not have native functionality for the model described in Sec. 3.2. Therefore, the fluid property calculations for the barotropic model were implemented externally. The density, viscosity, and speed of sound of the mixture along an isentropic expansion were precalculated. These values were then fitted with eighth order polynomials, achieving a relative deviation smaller than 0.01%. The fitted polynomials were implemented using the Fluent Expression Language.

The turbulence model selected for this study was the shear stress transport (SST) k − ω model [31]. This model is widely recognized for its accuracy in turbomachinery applications, and it proved to be more effective than the k − ε model in predicting flow separation in the inactive rotor blade passages. An additional advantage of the SST k − ω model implemented in Fluent is that it features automatic wall functions, allowing for accurate turbulence predictions near wall surfaces without necessitating precise control over the y⁺ distribution.

For the boundary conditions, the mixture model simulations require specifying total pressure, total temperature, and volume fraction at the inlet. In contrast, the barotropic model simulations only require the total pressure at the inlet. For both models, a turbulence intensity of 5% and a turbulence viscosity ratio of 10 were applied at the inlet. The static pressure was specified at the exit, while zero velocity and roughness were applied at the walls.

The convergence criterion for the simulations was that the root-mean-square residual for each equation should be lower than 10⁻⁶. Additionally, the exit velocity, inlet mass flow, and outlet mass flow were independently monitored as supplementary criteria for convergence.

A structured grid with quadrilateral elements was utilized for the simulations. The grid independence study for the 2D nozzle axisymmetric nozzle is detailed in Table 1. The middle columns compare the relative change of mass flowrate for each grid compared to the finest grid, for both the barotropic and mixture models. Additionally, the last column reports the average y⁺ value of each grid. Variations in mass flow are deemed negligible for meshes finer than grid 5, which was subsequently selected for the 2D simulations. To further assess the impact of near-wall mesh refinement on our results, supplementary simulations were conducted using three variants of grid 5, each maintaining the same number of elements but featuring distinct average y⁺ values (17, 23, 34). These tests confirm that modifications to the near-wall mesh refinement do not significantly alter the solution, indicating that the resolution of grid 5 close to the walls is adequate.

Table 2 and Fig. 4 show the exit velocity for the barotropic model, mixture model, and experimental data at various liquid-to-gas mass-based mixture ratios. The deviation in exit velocity ranges between 5.78% and 10.80% for the barotropic model and 7.54% and 14.11% for the mixture model, and it increases as the mixture ratio decreases in all cases. The experimental mass flowrate is only available for the case with a mixture ratio of 68.00, and the measured value is 3.657 kg/s. Both models accurately predict this mass flowrate, with the barotropic model and mixture model deviating by −0.19% and 1.56%, respectively. As the mass flowrate in a choked nozzle is determined by the conditions upstream of the throat, the very good agreement between simulation and experimental results suggests that both models accurately predict the performance of the converging section of the nozzle. However, the significant deviation in exit velocity implies that neither model precisely captures the behavior of the diverging section, where the flow becomes supersonic. This discrepancy may be caused by the unaccounted velocity difference between phases, particularly in the diverging section immediately following the throat, where the flow experiences significant acceleration. The difference in velocity arises because the gas phase, having lower density and inertia than the liquid phase, accelerates more easily. This leads to a higher velocity of the gas relative to the liquid, thereby increasing frictional losses due to drag between the phases and ultimately resulting in lower exit velocities. Furthermore, Elliot [2] estimated a gas-to-liquid velocity ratio of 1.46 using a one-dimensional, two-fluid model, supporting our reasoning.

Figure 5 illustrates the distribution of velocity and pressure within the nozzle for the case with a mixture ratio of 68.00. The flow undergoes acceleration owing to the pressure differential between the inlet and exit, and the maximum velocity is attained near the nozzle exit. Although both the barotropic and mixture models yield similar results, the mixture model predicts a slightly higher velocity in the divergent part of the nozzle. Additionally, the barotropic model predicts that the nozzle is slightly over-expanded, resulting in a small shock wave near the exit that aligns the outlet pressure with the boundary condition.

Table 3 provides further insights into the type of expansion occurring inside the nozzle. This table summarizes the barotropic model results at all mixture ratios, including the average velocity and pressure just upstream of the exit. It can be observed that the nozzle is slightly under-expanded for the lowest mixture ratio, while higher liquid ratios are associated with overexpansion.

The pressure profiles for various mixture ratios are illustrated in Fig. 6. In the converging section, all pressure profiles are nearly identical, exhibiting only negligible changes in volume fraction. By contrast, near the throat, the changes in volume fraction are significant, which is reflected in the steep gradient of the pressure profiles. For cases with lower mixture ratios, the changes in pressure and volume fraction near the throat are less steep. Beyond the throat, the flow undergoes a gradual expansion in the diverging section.

Figure 7 illustrates the fluid domain for the 3D Reynolds-Averaged Navier–Stokes simulations of the full turbine stage. Both the barotropic model and mixture model simulations were performed in Ansys CFX [32] due to its better convergence behavior in turbomachinery simulations compared to Ansys Fluent. The boundary conditions for the 3D simulations mirrored those of the 2D cases. At the nozzle inlet, the mixture model requires the specification of total pressure, total temperature, and volume fraction, whereas the barotropic model necessitates only the total pressure. Consistent with the 2D simulations, both models utilize the SST k − ω turbulence model with automatic wall functions and a specified inlet turbulence intensity of 5% and a turbulence viscosity ratio of 10. A frozen rotor boundary condition was used at the interface between the rotor and nozzle. At the rotor exit, an outflow condition enforcing average pressure and radial equilibrium was applied. No-slip conditions were applied to the nozzle walls, as well as to the rotor blades and hub and shroud surfaces in the rotating frame of reference. The no-slip condition at hub and shroud is consistent with Elliot's turbine configuration [2], featuring a shrouded rotor with no tip gap.

For numerical discretization, a first-order upwind scheme was used for both flow and turbulence models. The mesh independence study, detailed in Table 4 and discussed below, confirms that the discretization error is small despite using first-order methods. Regarding termination criteria, simulations were deemed converged when the root-mean-square residual for each equation fell below 10⁻⁷. In addition, the nozzle exit velocity and rotor torque were monitored as supplementary convergence criteria.

The 3D simulations use a hexahedral-dominated swept mesh for the nozzle and a multiblock structured mesh for the rotor domain. Special attention was given to achieving high mesh quality near the sharp corner caused by the oblique cut at the exit of the nozzle. The blade geometry was reconstructed based on the data from the Elliot report [2], and the rotor's mesh was generated using Ansys Turbogrid [33]. The grid independence study, detailed in Table 4, considers three grid levels: coarse, medium, and fine, with a refinement factor of approximately 1.5 between successive grids. The mesh sensitivity analysis, focusing on three key performance indicators (mass flowrate, nozzle exit velocity, and rotor torque), suggests that the relative changes in flow quantities between the medium and fine grids are minimal. Consequently, the medium grid was selected for all subsequent simulations, offering a significant reduction of computational cost while maintaining the accuracy of the results. As seen in Table 4, the average $y+$ value across the solid surfaces of the flow domain is approximately 8, indicating that the mesh resolution near the walls is adequate for capturing the influence of boundary layers.

In addition to the boundary conditions listed in Sec. 5.1, periodic boundary conditions were applied on the lateral ends of the rotor domain. Although this approach is not fully representative of the nonperiodic nature of the actual turbine geometry, it is supported by practical and computational considerations. In the examined turbine configuration, the majority of the torque is generated by the rotor blades directly exposed to the nozzle flow, termed the active blade passages. By contrast, the rotor segment outside the path of the nozzle jet, called inactive blade passages, has a minimal effect on torque generation but contributes to windage losses due to the friction with the surrounding environment. Considering the distinct behavior of active and inactive passages, the torque generation can be accurately modeled using a computational domain bounded by periodic boundary conditions that includes the active passages. This strategic approach circumvents the need for a full annulus simulation, drastically reducing the required execution time and random-access memory usage. The effects of windage losses associated with the inactive blade passages are further evaluated through separate CFD simulations, as elaborated in Sec. 5.5.

Ideally, the number of blade passages of the computational domain should be an integer submultiple of the total number of passages. The turbine considered has 151 blades, a prime number. To implement periodic boundary condition for scenarios involving 10, 15, and 25 blade passages, we considered a total of 150 blades. Out of the potential options, 15 was ultimately selected because this choice includes all the blade passages directly impinged by the jet and at least two adjacent passages on either side as a safety margin. To ensure that the predicted torque was not influenced by the number of passages considered, we carried out a separate simulation with 25 blade passages. This comparison indicates that increasing the number of inactive passages did not alter the flow field within the active segment (the results are not included in the paper for brevity), thereby confirming the validity of this modeling approach.

Figure 8 depicts the contours of absolute velocity inside the nozzle. The velocity profile is similar to that observed in the 2D axisymmetric simulations. The flow velocity increases slowly in the converging section. A significant change in velocity is observed near the constant-area throat section due to a change in the volume fraction. Subsequently, the flow continues to expand in the diverging section, attaining its peak velocity at the nozzle exit. Due to the nozzle being cut at a 20 deg angle from the axis, the flow experienced additional acceleration in the diverging channel situated between the circular and elliptical ends, reaching its maximum velocity at the elliptical end of the nozzle. The mixture model predicts a slightly higher exit velocity than the barotropic model. The mass averaged absolute velocity at the nozzle's exit for the mixture model and barotropic model is 104.40 m/s and 103.19 m/s, respectively.

Figures 9 and 10 compare the relative velocity and pressure contours of the barotropic model and the mixture model at a 50% span location and design angular speed. The results suggest that the nozzle directs the majority of the flow through nine blade passages, spanning from the 3rd to the 11th blades, where the most significant changes in flow properties take place. By contrast, the remaining blades are largely unaffected by the nozzle, exhibiting minimal variation in flow properties. In the case of the mixture model, the results indicate that the absolute velocity at the nozzle exit is marginally higher than that in the barotropic model, leading to an increased relative velocity at the rotor inlet. Despite this difference in inlet velocity, the flow pattern inside the rotor remains remarkably similar in both cases. The strength of the leading-edge bow shock, passage shock and expansion fans is slightly higher in the mixture model due to the increased inlet velocity.

Figure 11 compares the void fraction between the barotropic model and the mixture model at a 50% span location, showing remarkably similar distributions between the two. This shows that the barotropic model, for which the void fraction is a function of the pressure only, can correctly represent the flow evolution for the fluid of interest. Within the active blade passages, the lower vapor fraction found after the shockwave and on the pressure side is a consequence of the higher pressure reducing the volume of the nitrogen. The void fraction is overall slightly higher for the mixture model, which explains the higher relative velocity observed.

In this section, we analyze the main features of the flow field. Initially, our focus is centered on the flow characteristics at the design angular speed, where we explore several key aspects, including the spanwise distribution of the flow, shock system at the rotor inlet, blade loading distributions, formation of secondary flows at the endwalls, and main sources of loss (i.e., entropy generation). Following the examination of the design speed, we expand our analysis to assess the impact of the angular velocity on the flow field. To this aim, we analyze the flow at two additional angular speeds: one below and one above the design speed. All CFD flow visualizations presented herein are derived from simulations based on the barotropic model.

The formation of the leading-edge shock wave system is illustrated in Fig. 12, where the relative Mach number contours are shown at five different span locations. At a 5% span, the blade passages between 4 and 8 exhibit slightly higher flow compared to the other blade passages, even if the blades do not directly receive the jet from the nozzle at lower span locations. Nevertheless, there is a gradual increase in flow velocity toward the admission arc, attributed to the redistribution of the flow. This enhancement can be observed at the 25% span location. The rotor blade sections directly exposed to the nozzle flow experience either supersonic or transonic flow, observable at 50%, 75%, and 90% span. At the 50% span location, a leading-edge shock and passage shock form at blades 4–9. Owing to the leading-edge thickness, the leading-edge shock manifests as a detached bow shock that eventually coalesces with the passage shock. This interaction is clearly depicted in the right panel of Fig. 12. The strength of the leading-edge shock increases from blades 4 to 9 in correlation with the rise in flow velocity. The flow undergoes compression as it passes through the leading-edge shock, leading to a reduction in the Mach number. Notably, a small flow separation is discernible on the blade suction surface around the midchord, which is caused by the interaction of the shock with the blade surface boundary layer. The presence of an adverse pressure gradient across the passage shock induces a thickening of the boundary layer, prompting flow separation. The flow reattaches quickly downstream the shock-boundary layer interaction. As the inlet Mach number increases, the leading-edge shock angle progressively becomes more perpendicular to the flow. At the 75% span, a stronger flow separation is observed, primarily attributed to the higher intensity of the shock boundary layer interaction. This phenomenon is linked to the presence of a stronger passage shock, resulting in elevated losses and flow blockage. At the 90% span location, a mild leading-edge and passage shock forms at blades 4 and 5.

Figure 13 illustrates the relative Mach number contours at outlet of the rotor. In active blade passages, the flow pattern resembles that of a conventional turbine, characterized by secondary flows and corner flow separation near the hub and tip. As the rotor is a shrouded rotor with no tip gap and thus there is no tip leakage vortex. Flow redistribution occurs within the blade passage as it passes through. Due to the high incidence and insufficient flow energy inside the inactive blade passages, the flow is unable to follow the large flow turning; hence, it separates and forms a large separation vortex covering the entire blade passage from hub to tip (blades 11-15 and 1-2). In active blade passages, significant reverse flow is observed near the hub and shroud region. The flow separation near the hub is more intense, which can be primarily attributed to the relative location between nozzle and rotor blades. Since the nozzle flow impinges on the rotor above the midspan, it channels a relatively larger mass flow closer to the shroud, helping to mitigate the formation of corner separation vortices.

Figure 14 illustrates the blade loading distributions at 50% and 75% spans. To maintain clarity and avoid visual clutter, only selected blades are presented in the figures. The blade loading for inactive blades is found to be negligible. The results suggest that at the 50% span, a noteworthy variation in blade loading becomes apparent starting from blade 4. Subsequently, there is a consistent increase in blade loading from the blade 4 onwards, reaching its peak for blades 6 and 7, followed by a significant reduction beyond blade 10 as flow recedes. Examining blade 4 specifically, the loading remains relatively modest until 20% chord, after which it experiences a noticeable increase toward the trailing edge. On the pressure surface of blades 5–9, the pressure gradually increases from the leading edge, reaching a maximum near the midblade chord, and subsequently decreasing toward the trailing edge. Both peaks and valleys are observed on the suction surface of the blade. The presence of a passage shock over the suction surface causes a steep increase of pressure between 40% and 50% of the axial chord. The pressure drops to its lowest point prior to the passage shock, followed by a sharp increase in pressure across the shock. The flow separation due to shock boundary layer interaction reduces the diffusion capability of the turbine blades, impacting the overall performance and power output of the turbine. As the passage shock moves upstream from blade 4 to blade 9, the locus of the discontinuity also shifts upstream due to the change in position of the passage shock. A second valley in the pressure on the suction surface is identified around 70% of the chord location. This is attributed to the small flow separation and reattachment over the suction surface at this chord location.

The discussion regarding blade loading at 50% span also holds true for the 75% span. At 75% span, a noticeable change in loading is observed for blade 3 (not shown in Fig. 14). The locus of maximum pressure on the pressure surface shifts further downstream and is observed at around 60% of the chord length. Additionally, a significant reduction in loading is observed on the leading-edge portion of the blade (between the leading edge and the 50% span), in contrast to the observations at the 50% span location. This reduction in loading is attributed to lower diffusion and acceleration on the pressure and suction sides, respectively, resulting in an overall decrease in blade loading in the leading-edge portion. Notably, the second valley observed at the 50% span location vanishes at the 75% span. At this span, the trailing half of the rotor blade experiences a higher load.

Figure 15 provides a closer examination of the losses within the rotor passages, focusing on the entropy and entropy gradients at blade midspan. The entropy contours are indicative of the accumulated losses within the flow, while the entropy gradient contour represents the local generation of loss. As already depicted in Fig. 13, the inactive blade passages are entirely filled by the separation vortex, leading to a modest amount of entropy generation. As the flow in these passages is almost stagnant, entropy produced by viscous dissipation in the vortices accumulates and ultimately results in high entropy levels. As for the active blade passages, the main sources of loss can be attributed to three key phenomena. First, the losses at the leading-edge shock system, comprising both the bow shock and passage shock. Second, the flow separation regions caused by the shock-boundary layer interaction at about 40% of the blade chord. Lastly, substantial entropy generation occurs within the wakes and due to the mixing of the flow downstream of the blades. The cumulative effect of the loss mechanisms in the rotor blade passages is referred to as “rotor loss” in the quantitative loss analysis presented in Sec. 5.5.

Transitioning from the detailed analysis at the design angular speed, we now explore the impact of off-design operation on the flow field; Fig. 16 shows the relative Mach number contours at 50% span for three different rotational velocities. As demonstrated in Fig. 16(a), increasing the peripheral speed while keeping the absolute velocity of the jet constant reduces the incoming relative velocity. Consequently, a weaker leading-edge shock system is observed when the angular speed is 36% higher than the design speed. This weaker leading-edge and passage shock results in milder shock boundary layer interactions, preventing flow separation and eliminating the separation regions near the midchord. Despite this, a weaker shock forms over the suction surface of the active blades at around 60% of the blade chord, triggering flow separation in the rear part of the suction surface. The situation is reversed when the angular speed is lower than the design value. As illustrated in Fig. 16(c), the strength of the leading-edge shock and passage shock increases and moves upstream when the rotational speed is 49% of the design value. Due to the stronger passage shock, the adverse pressure gradient experienced by the boundary layer is more intense, resulting in a large flow separation bubble over the suction surface near the midchord. This phenomenon is clearly illustrated in the right part of Fig. 16.

In the preceding sections, we examined the flow field within the nozzle and rotor domains, identifying sources of loss from a qualitative perspective. Building on this foundation, this section aims to provide a quantitative analysis of all losses within the turbine, breaking down the impact of each loss mechanism on the turbine's overall performance.

The CFD simulations conducted with the barotropic and mixture models captured several key loss mechanisms, which can be grouped as (1) frictional losses within the nozzle, (2) rotor losses in the blade passages, and (3) kinetic energy at the rotor exit that is not converted to work. In addition to these, there are two significant loss mechanisms not modeled by the CFD simulations: (1) windage loss resulting from friction between the rotor and its surroundings, and (2) pumping loss (also known as end sector loss) caused by the unsteady emptying and filling of the rotor blade passages at the extremities of the active sector. Capturing the windage and pumping losses necessitates full annulus unsteady CFD simulations. However, such simulations, while technically feasible, are prohibitively expensive in terms of computational resources due to the large number of blade passages. Given that the simulations conducted in this study are steady-state simulations of a rotor sector, it becomes necessary to introduce correction terms to estimate the impact of these two mechanisms on turbine performance.

The direct computation of the windage loss is challenging due to two primary factors. First, it necessitates a full annulus simulation, which entails a prohibitive computational cost. Moreover, the experiments conducted by Elliot revealed that the inactive passage are almost entirely filled by the gas phase (nitrogen), a condition not accurately captured by the CFD simulations in this study due to their inability to model flow segregation. Consequently, the simulated fluid density within the inactive passages is significantly higher than the actual gas phase density. Since the windage loss is proportional to the fluid density in the inactive passages, performing full annulus simulations with either the mixture or barotropic models would lead to a significant overestimation of the windage losses. Considering these limitations, we evaluated the windage losses indirectly. First, we estimated the windage torque by means of a steady-state, incompressible simulation, employing nitrogen as the working fluid and reducing the pressure at the inlet of the nozzle to atmospheric pressure. This method mirrors Elliot's experimental approach, in which windage torque was determined by stopping the flow from the nozzle, while driving the rotor at the desired rotational speed with an electric motor. Utilizing an incompressible simulation without the energy equation is necessary to prevent divergence caused by the accumulation of frictional losses in the absence of net flow. The windage torque estimated with this procedure was finally converted into an equivalent enthalpy loss (⁠$Δhwindage)$ using Eq. (7).

where the corrected enthalpy change can be identified with the actual enthalpy change in Eq. (8).

The outlined calculation procedure was employed to estimate the overall turbine performance and the contribution of each loss mechanism at the different angular speeds measured by Elliot [2]. Figure 17 presents the breakdown of the losses obtained for the mixture model expressing the losses as total-to-static efficiency penalties. The losses within the nozzle are insensitive to rotor angular speed, consistent with expectations for a supersonic nozzle where downstream disturbances cannot propagate upstream. By contrast, rotor losses decrease as angular speed increases, a finding that is consistent with the insights obtained from the qualitative analysis of Fig. 16. Although rotor losses decrease at higher angular speeds, the condition of minimum overall loss occurs at the nominal speed of 1652 rpm. This is due to the increase of kinetic energy, windage, and pumping losses at higher speeds. Notably, the windage torque predicted through the indirect CFD procedure increases quadratically with rotational speed, which is consistent with expectations from dimensional considerations.

Figure 18 compares the experimental torque values with simulation data obtained from the barotropic and mixture models, both before and after applying corrections for pumping and windage losses. These corrections are derived by converting the enthalpy changes calculated from Eq. (16) into torque values through Eq. (7). Initially, both models significantly overestimate torque, with deviation of about 28% at design speed and narrowing to approximately 15% at half the design speed. Including the correction terms significantly enhances the model accuracy, reducing the deviation with respect to the experimental values below the ±10% uncertainty margin of the experiments. In particular, the torque deviation at the design speed is about 8%, and it narrows further to approximately 3% at half of the design speed. These findings underscore the critical influence of pumping and windage losses on the performance of the two-phase turbine configuration analyzed in this work.

Despite these improvements, a systematic overprediction of torque persists, which cannot be explained by the effects of pumping and windage losses alone. This discrepancy is likely associated with the overestimation of nozzle exit velocity. As reported in Table 2, the barotropic and mixture models predict velocities that are 7.23% and 9.15% higher than the experimental values, respectively. This overestimation of nozzle exit velocity leads to higher kinetic energy at the inlet of the rotor, which in turn contributes to the torque overprediction. As discussed in Sec. 4.2, the discrepancy in exit velocity may be caused by the effect of slip, particularly in the diverging section of the nozzle where the flow experiences a significant acceleration. Therefore, refining the models to include the influence of slip stands as a promising avenue for enhancing the accuracy of nozzle exit velocity predictions, which, in turn, could further align the torque with the experimental data.

This work evaluated the mixture and barotropic models for the simulation of a two-phase turbine, comparing numerical predictions at different liquid-to-gas mass ratios and angular speeds against experimental data from a single-stage impulse turbine operating with a working fluid consisting of a mixture of liquid water and gaseous nitrogen. The following conclusion are drawn from the results of this study.

Both the barotropic and mixture models demonstrate remarkable accuracy in predicting the mass flowrate, deviating by only −0.19% and 1.56%, respectively for the case with 68.00 liquid-to-gas mass ratio. Given that the mass flowrate in a choked nozzle is only affected by the conditions upstream the throat, the good agreement between simulations and experimental data indicates that both models accurately predict the flow in the converging section of the nozzle.

However, both models exhibit a deviation in the nozzle exit velocity, ranging from 5.78% to 10.80% for the barotropic model and from 7.54% to 14.11% for the mixture model, with the discrepancy increasing at lower mixture ratios. This discrepancy is likely caused by the effect of slip (velocity difference between the phases), especially in the portion of the diverging section close to the throat, where a rapid acceleration of the flow takes place. A proper understanding and modeling of the effect of slip effect on nozzle performance requires further analysis.

The internal flow within the turbine is complex and inherently three-dimensional. The active rotor passages are characterized by a supersonic inlet and leading-edge bow shocks that merge with the passage shocks. The interaction between the entry shock system with the boundary layer on the suction surface of the blades leads to boundary layer thickening and, in some cases, to severe flow separation. Conversely, within the inactive blade passages, the flow undergoes complete separation, forming multiple recirculation zones and resulting in negligible blade loading. The complexity of the flow phenomena, compounded by the two-phase effects, calls for a robust modeling approach. The CFD solver used in this research demonstrates its robustness, successfully converging for both the mixture and barotropic models. The capability of more sophisticated multiphase models, such as the two-fluid Eulerian model, to achieve reliable convergence within the demanding context of two-phase turbines remains to be evaluated.

The predictions both from the barotropic and mixture models are remarkably consistent. The differences in key metrics such as torque and mass flowrate are below 1% across all cases, with the barotropic model being closer to the experimental values. Moreover, the barotropic model is less computationally expensive and more robust as it only requires the solution of the mass and momentum equations along with a turbulence model. By contrast, the mixture model requires the inclusion of the energy equation and a transport equation for the void fraction. The combination of simplicity and accuracy positions the barotropic model as a promising method for simulating two-phase flows in turbomachinery.

While both models successfully capture the approximately linear relationship between torque and angular speed, they significantly overpredict the torque values. This discrepancy primarily arises from the inability of the blade sector, steady-state CFD simulations to capture the windage friction and the unsteady effects associated with partial admission. By introducing correction terms to account for these additional losses, the deviation against the experimental data is substantially reduced and falls within the uncertainty range of ±10% across all cases. Despite this adjustment, a persistent overprediction of torque remains, which is likely caused by the overestimation of nozzle exit velocity.

The simulations conducted in this study highlight various limitations of the current modeling approach, leading to opportunities for further refinement and exploration. Incorporating the effect of slip in the simulation can potentially enhance the prediction of nozzle exit velocity and rotor torque. Employing a correlation for pumping loss has been shown to be advantageous in refining the torque predictions from steady-state simulations. However, given that the existing correlations were originally developed for single-phase turbines, a calibration process is needed to adapt these corrections for two phase turbines.

The modeling approaches examined in this research hold promising potential for advancing the simulation of two-phase turbomachinery. Both the barotropic and mixture models prove computationally robust in handling the large expansion ratio in the nozzle and the complex supersonic flow within the rotor. The predictions of both models are remarkably similar and show good agreement with the experimental data used for validation. Furthermore, this work identifies key areas where model improvements are necessary to increase the accuracy of the predictions. We believe that the present work constitutes a significant initial step toward the realization of a reliable and computationally efficient CFD model for the simulation of two-phase turbines.

• “Sustainable large-scale energy storage in Egypt” by the Ministry of Foreign Affairs of Denmark and administrated by Danida Fellowship Centre (Project No. 21-M13-DTU; Funder ID: 10.13039/501100001731).

• European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie (Grant Agreement No. 899987; Funder ID: 10.13039/100010663).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.