CFD analysis of flashing flow in two-phase geothermal turbine design


 A thermal power plant for the East African Rift countries is under study for combined energy and freshwater generation using geothermal water, available at above 500 kPa pressure and temperature exceeding 150°C. This article presents the computational fluid dynamics (CFD) model and analysis of the two-phase turbine used for power generation in this total flow thermal plant. Flash boiling was implemented using a two-fluid multiphase model with the thermal phase-change criteria for heat, mass, and momentum transfer in the CFD solver ANSYS CFX. Initially, flashing flow in a converging–diverging nozzle was validated. This stationary nozzle model was then extended to a curved rotating nozzle reaction turbine and the results of flow and power were evaluated against available test data at 400 kPa feed water pressure under subcooled condition of 117°C and a very low backpressure of 6 kPa. Flow through this turbine was predicted within 8% deviation. An overestimate in thermodynamic power by 30–50% was predicted at speeds below 4000 rpm, while at the design speed of 4623 rpm the deviation was less than 5%. Rotor torque and hence power estimate was found to be dependent on the bubble size, bubble number density, and heat transfer parameters prescribed in the CFD model. The vapour dryness fraction at turbine exit was close to an isentropic expansion vapour quality. The isentropic efficiency was 7.5–17% for the analysed speed range.


V:
velocity ω: angular velocity of the frame ρ: densitȳτ q : phaseq stress-strain tensor μ q : shear viscosity of phase q F q : external force V pq : interphase velocities − → V rel : relative velocity between phases pq : heat transfer coefficient at the phasic interface q q : heat flux h pq : interphase enthalpy A i : interfacial area Nu p : Nusselt number Re: Reynolds number J a T : modified Jakob number N b : bubble number density (BND) P t : turbine power x: vapour dryness fraction α: volume fraction − → V t : translational velocity of the frame k q : thermal conductivity of phase q P : pressure λ q : bulk viscosity of phase q R pq : interaction force between phases h q : specific enthalpy of phase q d q : bubble diameter of phase q m pq : mass transfer rate from phase p to phase q Q pq : heat exchange between phases p and q g: g r a v i t y T sat : saturation temperature P e: Péclet number P r: Prandtl number a p : thermal diffusivity of liquid phase N t : turbine speed T n : nozzle torque η i : isentropic efficiency Subscripts p, q: phase indices q: water vapour p: water liquid sat: saturation

Introduction
Geothermal hot brine deposits are a natural source of hightemperature and high-pressure fluids. Across the globe, techniques and systems to utilize the potential of geothermal wells are being explored. Output electrical energy or mechanical energy is found to be the preferred means of extracting the thermal energy contained in the geothermal well. At sites such as in Iceland, Africa, America, Italy, and Australia, several projects are underway. In the present investigation, a thermodynamic system for combined power generation and pure water cogeneration has been aimed for. This is due to the demands at sites in East Africa, Ethiopia, and Kenya where geothermal wells can be tapped, but the local population is in greater demand for drinkable water, as is the power requirement. In the context of mechanical power generation from geothermal fluids, an impulse-type steam turbine has been a classical approach. The systems employing binary fluids such as organic fluids have found use in low-temperature heat sources due to their flammability at higher temperature. In high-temperature sources, the geothermal water is initially flashed and steam so produced is used to drive the turbine. The total flow concept that directly uses the brine as the turbine's working fluid is the most efficient means of energy conversion but suffers from the backdrop of heavy mechanical erosion of the turbine and nozzle and mineral deposition, which makes the system inoperational in due course (Smith, 1993). Austin, Higgins, and Howard (1973) have presented the total flow concept for power recovery for application to the Salton Sea geothermal bed. Geothermal water was expanded from a pressure of 2200 psia at a temperature of 300 • C to 400 psia. Thermal energy was converted to kinetic head by expansion through a converging-diverging (CD) nozzle. The jet was then used to drive an impulse turbine. Similar concepts can be found in the literature that have used multiphase nozzles. Elliot (1982) has reported tests with water and nitrogen mixtures and with single-stage and two-stage impulse turbines. A special two-phase nozzle with mixer housing and liquid injection tube was designed to vary the mass ratio of liquid to vapour in the two-phase flow investigation. Fabris (1993Fabris ( , 2005 identified that the high inlet pressure to the two-phase nozzle entry, resulting from the centrifugal acceleration of the liquid flowing radially outwards, thereby subcooling it and delaying the flashing, was a deficiency of House's (1978) radial outflow reaction turbine configuration. As an improvement, he proposed a nozzle configuration as shown in Section 4, Fig. 5, to minimize the lateral forces on the accelerating and expanding fluid and to maximize the length of the nozzle in which expansion took place. Efficiencies of up to 50% have been reported from tests carried out on a turbine with this configuration (Fabris, 2005). It was claimed that efficiency could be further improved by additional stages (Comfort, 1978). More recently, Date, Khaghani, Andrews, and Akbarzadeh (2015) and Date, Vahaji, Andrews, and Akbarzadeh (2015) have presented experimental results of a curved nozzle total flow reaction turbine for low-temperature feed water application and reported a maximum isentropic efficiency of 25%. Abuaf, Wu, Zimmer, and Saha (1981) conducted an extensive experimental study on two-phase flow in a CD nozzle. The test data are available at various operating conditions in the form of pressure measurements across several points along the length of the nozzle. Vapour void fraction measurements have also been documented in this report and the data are highly suitable for conducting a computational model validation. A closedloop system to heat, pressurize, and control the water flow and condenser pressure and temperature was used in these experiments. The same nozzle has been analysed using the thermal phase-change model in Liao and Lucas (2015, 2018. They have extensively reported the influence of various modelling parameters, closure models, and grid dependence. Results have been validated in the form of pressure drop along the nozzle and vapour void fraction along the nozzle and across several transverse planes and compared with the BNL test data. A good match of the overall mass flow rate through the nozzle, pressure drop, and average void fraction at the exit was reported. In their simulations, the liquid-side heat transfer coefficient was formulated based on the Jacob number and Péclet number, while the vapourside heat transfer coefficient was formulated using the zeroresistance approach across a spherical bubble surface whose diameter was specified using BND. Grid refinement reported in this study showed an optimum node count of 744 800 for the nozzle geometry and the numerical solver used was ANSYS CFX in steady state. Water data for both liquid and vapour phases were specified using IAPWS-IF97 (Wagner et al., 2000). The SST k-omega turbulence model was used. Several test points and operating conditions were validated in this study. Particularly, discrepancies in the distribution of radial void fraction were reported. This literature provides a foundation for the modelling approach used in the current study. In order to evaluate the radial void fraction discrepancies in detail, Janet, Liao, and Lucas (2015) studied the influence of nucleation model using a bubble number transport equation instead of directly fixing a BND as was done by Liao and Lucas (2015). Small improvements in prediction were observed but were not satisfactorily validated with the BNL measurements. In addition to bubble number transport formulation, Janet et al. (2015) also considered the bubble breakup and coalescence (Edwards & O'Brien, 1970;Marsh & O'Mahony, 2008) in their model.
A thermal power plant for the East African Rift countries is under development for combined energy and freshwater generation using geothermal water available at >500 kPa pressure and temperature exceeding 150 • C. The present work deals with the implementation of the computational fluid dynamics (CFD) model and analysis of the two-phase turbine used for power generation of this total flow thermal plant. Here, a Eulerian-Eulerian multiphase model has been used in the numerical solver ANSYS CFX. Two-phase heat, mass, and momentum transfer has been modelled using the thermal phasechange model. A phase-change validation study has been conducted on published experimental data on a CD nozzle flow at various operating conditions (Abuaf et al., 1981;. Satisfactory results of flash boiling of liquid water to vapour were achieved by the model and quantitative data were validated with the literature. The model was then extended to the test turbine design reported in  and the results have been evaluated.

Governing Equations of the Two-Phase CFD Model
In the Eulerian-Eulerian approach of multiphase modelling, the primary and secondary phases are treated as interpenetrating continua. Phasic volume fraction α q represents the volume of a computational cell occupied by phase q (water vapour) and the conservation transport equations of mass, momentum, energy, and other scalars are satisfied by each phase individually. Conservation of mass: Continuity equation for phase q is Here, V q is the velocity of phase q and has three Cartesian components.ṁ pq is the mass transfer term from phase p (water liquid) to phase q (water vapour).ṁ qp is the mass transfer term from phase q (water vapour) to phase p (water liquid). S q is any additional source for phase q. Conservation of momentum: Momentum conservation for phase q is whereτ q is the phase q stress-strain tensor with components that are functions of shear and bulk viscosity. P is the pressure shared by all phases. μ q and λ q are the shear and bulk viscosity of phase q, respectively. F q represents the sum of external body forces such as Coriolis and centripetal accelerations (α q ρ q [ ω × ( V q − V t )]) and interaction forces such as lift, virtual mass, turbulent dispersion, etc. acting on phase q. ω is the angular velocity and V t is the translational velocity of the moving reference frame relative to stationary reference frame. R pq is the interaction force between phases. V pq and V qp are the interphase velocities.
Conservation of energy: Enthalpy conservation equation for phase q is Here, h q is the specific enthalpy of phase q. q q is the heat flux, S q is the enthalpy source, and Q pq is the intensity of heat exchange between phases p and q. h pq is the interphase enthalpy, i.e. vaporization enthalpy in case of evaporation given by equation (7).
pq is the heat transfer coefficient at the phasic interface, k q is the thermal conductivity of phase, Nu p is the Nusselt number, and d p is the vapour bubble diameter.

Thermal phase-change model -interphase mass transfer
The thermal phase-change model is based on thermal nonequilibrium between the phases. Heat, mass, and momentum transfer is governed by the local heat balance across the liquidvapour interface during phase transition. Such a formulation is suitable for modelling of phenomenon such as flash boiling occurring in CD nozzle flows. Within the Eulerian multiphase framework, the thermal phase-change model defines the mass transfer rates entirely based on the interphase heat transfer and assumes a mechanical equilibrium across the liquid-vapour interface . At the water liquid and water vapour interface, heat balance gives From the interface to water liquid phase, From the interface to water vapour phase, Enthalpy with subscript 's' is selected depending on evaporation or condensation: The interphase mass transfer in equations (1)-(3) is obtained from equation (4). Equations (5) and (6) are used to select the enthalpy for source in equation (3). Ifṁ pq > 0 (boiling), Ifṁ pq < 0 (condensation), As seen from equation (4), heat transfer coefficients p and q and interfacial area A i are required to be prescribed in the thermal phase-change model.

Interfacial area density
The particle model is adopted for interfacial transfer between two phases, which assumes that one of the phases ( p, here water) is continuous and the other (q, steam) is dispersed. The surface area per unit volume is then calculated by assuming that steam is present as spherical particles of mean diameter d q . Using this definition, the interfacial area density is Either mean diameter or the BND has to be specified. They are either prescribed based on experimental comparison or obtained by solving additional transport equation for conservation of BND (N b ). This transport equation has quantities that account for nucleation, bubble growth, break-up, or coalescence coupled with local flow conditions (Janet et al., 2015;Marsh & O'Mahony, 2008).

Interphase heat transfer
The heat flux across the spherical bubble surface that is transferred to water from the water-steam interface is Here, p is the overall heat transfer coefficient. These inputs are prescribed by specifying the interface Nusselt number correlation (Table 4). This is further discussed in Section 4.2.

Numerical settings of the flow solver
The coupled volume fraction algorithm of two-phase flow with implicit scheme of coupling velocity, pressure, and volume fraction equations has been used in this work in the flow solver ANSYS CFX. A high-resolution scheme (blend of first-order and second-order upwind schemes) was specified for discretization of advection terms and for the transient discretization a fully implicit second-order backward Euler scheme was used. The SST k-omega turbulence model was used as a homogeneous specification for both water and steam and the first-order upwind scheme was used for discretization. The steady-state solver was used with a timescale limit factor, which in case of the BNL nozzle flow case was 1e −4 s, while in case of the curved nozzle turbine a lower limit of 1e −5 s was required for stability of the solution. Special functions given by equations (10) and (12) were used to initialize the pressure and equations (11) and (13) for volume fraction in the domain as dependent on nozzle centreline location. For the BNL nozzle, the axial length and for the turbine the radius were used as the local parameters of which pressure and volume fraction were functions.
BNL nozzle pressure initialization function: BNL nozzle volume fraction initialization function: Turbine curved nozzle pressure initialization function: Turbine curved nozzle volume fraction initialization function: Liquid and vapour temperature initialization function (Wagner et al., 2000): 2 − 4 (n 9 + n 10 D) Liquid and vapour temperatures were initialized with equation (14) of saturation temperature defined by IAPWS-IF97 (Wagner et al., 2000) and local pressure (equations 10 and 12) was used as input to this equation in the factor β. Thus, initially a thermal equilibrium is forced on the system, and as numerical iterations progress, the two-phase flow pattern materializes. The convergence criterion was set to an rms level of 1e −4 for all equations and iterations were advanced until flow-averaged quantities such as feed water mass flow and nozzle exit vapour mass fraction stabilized. Abuaf et al. (1981) conducted an extensive experimental study on two-phase flow in a CD nozzle. A closed-loop system to heat, pressurize, and control the water flow and condenser pressure and temperature was used in these experiments. This section describes the CFD model of this symmetric CD nozzle and application of the thermal phase-change model to estimate the flash boiling and vapour void fraction distribution in the test section. In the literature, two-phase CFD modelling of the BNL nozzle has been analysed using the thermal phase-change model of ANSYS CFX solver in Liao and Lucas (2015). Same solver settings, water liquid and vapour-fluid properties using the IAPWS-IF97 database (Wagner et al., 2000), and numerical scheme have been applied in the current study in order to achieve consistent comparison. The Ranz-Marshall correlation (Ranz & Marshall 1952) for interphase heat transfer (Table 4) has been prescribed. Figure 1 shows the hexahedral computational mesh used for the 3D nozzle domain. Grid is refined closer to the throat region. Inlet boundary is defined using measured pressure, temperature, and saturated liquid condition. Outlet boundary is defined using measured pressure. Pressure and temperature initial conditions are specified using a geometric function and saturation temperature equation according to IAPWS-IF97 coefficients (Wagner et al., 2000). These initialization functions help in stability of the steady-state solver. About 15 000 iterations of the solver are computed to achieve residual imbalance below 1e −4 for all the transport equations. Three grid refinement levels were evaluated to achieve mass balance within 0.5% accuracy.

Results and validation
Results from the BNL nozzle experiments (Abuaf et al., 1981), the CFD results presented in the literature , and current computational model results have been compared here in order to validate the two-phase model. Results of pressure and vapour void fraction variation along the nozzle length have been quantitatively compared. Phase temperature variation in the domain has also been presented.
Experimental runs referenced as BNL309 and BNL296 have been computed in the current study. Table 1 presents the grid independence results evaluated for BNL309. Mesh sizes and boundary layer statistics have been reported. As the grid is refined from mesh 1 to mesh 3, the intake liquid mass flow rate through the nozzle tends to increase.
Mass imbalance is calculated between the intake and exit planes of the nozzle and it is seen that from mesh 1 to mesh 3 the imbalance reduces from ∼1% to <∼0.5%, which was considered to be an acceptable grid refinement level. Table 2 presents the operating conditions and results of the nozzle flow for the two test sets evaluated. In comparison to the measurements (ṁ; Abuaf et al., 1981), the CFD results reported by Liao and Lucas (2015) and the current computations are underestimated by ∼5%. Both the CFD models of Liao and Lucas (2015) (ṁ) and the current study (ṁ) closely predict the flow at both the operating conditions. Figure 2 is a plot of pressure variation and water vapour void fraction distribution along the centreline of the nozzle for BNL309. The nozzle inlet pressure is 555.9 kPa and saturation temperature is 422.25 K. Pressure steadily drops from 555.9 to ∼380 kPa in the converging section of the nozzle with minimum pressure occurring at the throat. Downstream of the throat, the pressure recovers to ∼402 kPa and then remains close to this magnitude in the diverging section of the nozzle. The nozzle exit and condenser pressure in this test was 402.5 kPa. Pressure variation follows close to that of measurements and the CFD results reported by Liao and Lucas (2015). In the converging section of the nozzle, water remains in liquid phase. At the throat, the flash boiling phenomenon initiates and in the diverging section of the nozzle phase change from liquid to vapour continues with vapour volume fraction increasing to ∼75% at the nozzle exit. Vapour void fraction follows close to that of measurements and CFD results reported by Liao and Lucas (2015). Similar validation of pressure and vapour void fraction was achieved in case of BNL296 operating conditions. Figure 3 is a contour plot of water liquid and vapour volume fraction distribution in the longitudinal plane of the nozzle. Radial variation of the vapour void fraction is observed in the contours. At the throat, flashing initiates with high vapour concentration at the nozzle centreline and gradually reducing radially outwards. In the diverging section, the vapour void fraction continues to increase to a value of 75% at the nozzle exit. Figure 4 is    a comparison of local saturation temperature (T sat ) according to IAPWS-IF97 (Wagner et al., 2000), local water liquid temperature, and local water vapour temperature in the longitudinal plane of the nozzle. Liquid at inlet is subcooled by ∼6 • C and hence in the converging section of the nozzle there is no flashing of the liquid. Just upstream of the throat when flashing initiates due to fall in saturation temperature, liquid temperature begins to drop. In the diverging section of the nozzle, liquid temperature continues to drop and the exit temperature is ∼418 K. The generated vapour is at local saturation temperature. From the nozzle throat up to the nozzle exit, it is seen that liquid temperature is dropping. This drop in temperature is due to the loss of latent heat (h pq ) in equation (3) transferred to the vapour phase. Results presented in Table 2 and the local pressure, temperature, and vapour void fraction distribution thus validate the twophase formulation of the thermal phase-change CFD model.  have presented experimental performance of a curved nozzle two-phase reaction turbine. Their curved nozzle design was based on the curvature calculation procedure described in Fabris (1993Fabris ( , 2005 with variable circular cross section such that pressure drop along the length of the diverging section is assumed to be linear and the change in relative velocity of the fluid along the length of the nozzle is also assumed to be linear. The two phases tend to separate within the curved nozzle due to the large lateral acceleration, and by designing the nozzle curvature using these assumptions, the lateral components of Coriolis acceleration, centripetal acceleration relative to motion, and centripetal acceleration with respect to turbine axis can be balanced. Figure 5 shows the layout of the curved nozzles inside the turbine disc and a mid-plane section of the turbine model. Throat of the curved nozzle is located at a radius of 59.15 mm. The full span of the nozzle from inlet to the exit is 240 • .  have reported test results from two operating feed water temperatures: for first test the average feed water temperature was maintained around 97 • C under local atmospheric pressure, and for the second test the average feed water temperature was maintained around 117 • C. For both the tests, the initial condenser pressure was maintained around 6 kPa. The maximum power output of the turbine was estimated to be around 1330 W with an isentropic efficiency of around 25%. The turbine operation was not continuous due to limits of the condenser to handle higher mass flow rates, which eventually raised the 6 kPa backpressure with time. The test lasted for ∼150 s of acceleration and further 100 s of deceleration of the turbine.

Two-Phase Reaction Turbine Analysis
The results from the CFD model of the turbine with feed water temperature of 117 • C have been discussed in this section.

Base turbine design and CFD model
The thermal phase-change, two-phase model established in the BNL nozzle study is also applied in case of the test turbine analysis. There are two main differences as compared to the BNL nozzle model. (i) A moving reference frame approach is used to apply the rotational speed to the turbine. (ii) The turbine has two symmetrically positioned curved nozzles. In order to reduce CFD domain, a 180 • sector of the rotor is modelled with periodic boundaries. Inlet to the nozzle is defined by the feed water pressure of 400 kPa and temperature of 117 • C. The exit of the turbine is open to a flash tank connected to condenser coil, which during the start of the experiments is vacuum conditioned to 6 kPa pressure. The available measurements  on the base turbine are listed in Table 3 and are used as boundary conditions for the CFD model at corresponding operating speed.
The inlet pressure (P in ) is maintained fixed, while the exit pressure (P out ) is specified as recorded during the measurements. Inlet feed water temperature (T in ) is specified with liquid volume fraction as 1 and vapour volume fraction as 0. Exit temperature is specified only for backflow definition. The isentropic efficiency listed in Table 3 is based on the gross shaft power and the isentropic enthalpy change due to the expansion from P in to P out . The CFD model uses a smaller flash tank domain up to a diameter of 225 mm, which is defined as pressure outlet boundary condition. CFD computations were performed over a range of turbine operating speeds from 1560 to 4623 rpm. Results from the measured data for feed water flow rate, shaft power, exit vapour dryness fraction, and isentropic efficiency have been used for the validation of the two-phase turbine model. Table 3: Measured performance data of the base two-phase turbine .   Figure 6 presents the computational domain of the twophase turbine and the main boundaries. The inlet is specified as pressure acting on an axial cross section of the inlet pipe. Only one of the curved nozzles is included in the CFD model as a 180 • sector of the turbine has been modelled. Consequently, there are periodic boundaries. In order to reduce the grid size further, a symmetry plane is defined at the mid-plane of the turbine. Due to this condition, gravity body force acting along the turbine rotational axis has been excluded from the two-phase model. A moving reference frame with steady-state solver applied angular velocity to the turbine.

Interphase heat transfer and turbine performance
Computational modelling of flashing flows is highly complex. Although the thermal phase-change model provides a frame-work within the Eulerian flow solver to include the mechanism of phase change and the accompanied heat, mass, and momentum transfer, it relies completely on the specification of interphase heat transfer. As seen in equations (1)-(3), the mass transfer and enthalpy change are closely coupled. In literature sources such as Janet et al. (2015) and Liao and Lucas (2015, 2018, it has been discussed and also in the present study it was experienced that a reliable correlation for estimate of interphase heat transfer coefficient in conjunction with the parameters of bubble size, growth, and nucleation is crucial. In Liao and Lucas (2018), test cases have been reported comparing predictions of bubble growth rate with test data and direct numerical simulation data for different Nusselt number correlations. The closeness of result was found to be dependent on the Jakob number and Reynolds number. Some of the analytical correlations for bubbles rising and translating in stagnant liquid with certain degree of superheat take into account the effect of conduction via an analogy with 1D unsteady thermal diffusion with moving boundary. Convection heat transfer is accounted via the penetration theory (Liao & Lucas, 2018). The effect of turbulence in experimental studies showed an increase in the Nusselt number with increase in turbulence intensity. Such effects have been used to modify the convection heat transfer factors using turbulence timescales. The Nusselt number is defined in terms of the Jakob number (J a T ) and Péclet number (P e) in these correlations. Other empirical correlations such as the Ranz-Marshall correlation (Ranz & Marshall, 1952) are based on experimental data such as spherical water droplets evaporating in an air stream. The Nusselt number is defined in these types of correlations in terms of the Reynolds number (Re) and Prandtl number (P r).
Though the Ranz-Marshall heat transfer correlation (equation 15) performed well in the case of the BNL nozzle study, a  Aleksandrov, Voronov, Gorbunkov, Delone, and Nechayev (1967) Nu p = ( 12 π 2 J a T 2 + 1 3π P e) 0.5 suitable Nu p correlation was not found for flashing flows in rotating channels; hence, multiple formulations have been examined for consistency and accuracy. Table 4 presents the correlations for Nusselt number that have been investigated in the current two-phase turbine study. From the literature survey, it was found that broadly the Nusselt number correlations for interfacial heat transfer are based on either Re and P r or P e and J a T as classified in Table 4. The Jakob number J a T directly accounts for the degree of liquid superheat. J a T represents the ratio of the sensible heat for a given volume of liquid to heat or cool through T sup = (T w − T sat ) in reaching its saturation temperature to the latent heat (L = h pq ) required in evaporating the same volume of vapour. Ranz-Marshall (Ranz & Marshall, 1952) type of correlations have been extensively used in flows with flashing and are based on the Reynolds number and Prandtl number. However, some other flashing studies , 2018 have identified Wolfert (1976) type of correlations (equation 18), using J a T , as superior for vapour bubbles in superheated liquid. BND is unknown for the turbine flow regime. BND controls the size variation (d v and A i in equations 1-4) of the flashed vapour according to equation (8) and hence is critical to achieve right flashing development. The BND value in the turbine nozzle analysis was selected based on the available literature on BNL CD nozzle flow as 5e +07 and 5e +10 m −3 and the results have been compared here between the two values. Figure 7 compares the estimate of turbine power between the two types of heat transfer correlations listed in Table 4. Thermodynamic power is calculated using nozzle torque in equation (20). Figure 7a compares the Ranz-Marshall type of correlations and Fig. 7b compares the Wolfert type of correlations. Within the family, these correlations are consistent and the predictions are close to each other. Similar observation was made with respect to the feed water flow rate estimate. However, there are differences when they are compared across the family. Figure 8a and b presents the comparison of turbine feed water flow rate between measurements  and CFD model for a presumption of BND = 5e +10 and 5e +07 m −3 , respectively. Results based on Ranz-Marshall and Wolfert correlations for Nu (Table 4) have also been compared. With a BND of 5e +10 m −3 (Fig. 8a), turbine flow at various operating speeds was underestimated within 10-18% deviation from measurements. However, with a BND of 5e +07 m −3 (Fig. 8b), turbine flow at various operating speeds was well estimated within 0.5-7.5% deviation from measurements.
In measurements , gross shaft power is expressed as the sum of shaft power extracted at the load dynamometer, frictional power loss, and the rate of change of kinetic energy at a given speed. This shaft power is equal to the power estimated by rotor torque from the CFD model. Figure 9a and b presents the comparison of turbine gross shaft power between measurements and CFD model for a presumption of BND = 5e +10 and 5e +07 m −3 , respectively. In addition, results based on Ranz-Marshall and Wolfert correlations for Nu (Table 4) have been compared. From Fig. 9b, with a BND of 5e +07 m −3 , the CFD model overpredicts the turbine power at operating speeds <4000 rpm in the range of 30-55%. At 4623 rpm, turbine power deviation in the range of 4.5% was observed. However, with a BND of 5e +10 m −3 , the overprediction was high, in the range of 20-90%. This suggested that a BND of 5e +07 m −3 or lower was more suitable for the two-phase regime formed inside this turbine nozzle domain. In the BNL nozzle validation by Liao and Lucas (2015, 2018 and other flashing flow experiments (Janet et al., 2015;Marsh & O'Mahony, 2008), it has been highlighted that an optimum BND specification is required at a given operating pressure and temperature. Similarly, the current turbine nozzle flow needs an optimum BND to be specified for each operating speed in order to achieve rotor torque and power closer to the measurements.
There was very small difference in flow prediction between the two types of Nu correlations ( Fig. 8a and b). With respect to turbine power prediction, however, the response from the two correlations was greatly different. Comparison of Fig. 9a and b shows that the power magnitude drops in case of the Ranz-Marshall correlation when the BND was changed from 5e +10 to 5e +07 m −3 . Also, the trend of power versus speed variation is largely deviated with a BND of 5e +07 m −3 . In case of the Wolfert correlation, there was a drop in magnitude of power, but the trend with speed was consistent.
The Ranz-Marshall correlation is an empirical correlation based on experimental data for spherical water droplets evaporating in an air stream. The Nusselt number here is defined in terms of the Reynolds number and Prandtl number. The degree of liquid superheat is missing in both these flow numbers and it is expected to be accounted by the empirical exponents and coefficients of the correlation. In the BNL nozzle case presented in Section 3, these factors worked reasonably well. However, in case of the turbine, upstream pressure is >400 kPa and downstream pressure close to 8 kPa. Thus, operating condition has a severe expansion ratio of 50. In the absence of a direct account for liquid superheat, as it exists in the form of Jakob number in the Wolfert correlation, the Ranz-Marshall correlation became highly sensitive to the size of the vapour bubble. Bubble size distribution controls the heat transfer interface area between the two phases. As such when BND was changed, the Ranz-Marshall correlation fails to reproduce the two-phase regime formed under this operating condition and resulted in an inconsistent turbine power variation with speed. If Fig. 9b is further studied, the deviation increases with speed. This is so because with increase in speed, the time available for bubble growth inside the nozzle reduces and this can lead to different bubble size distribution even though the operating pressure ratio is the same. The Wolfert correlation was thus found to be more acceptable choice for further investigation. Both BND and Nusselt number correlation factors require significant tuning at each operating condition or speed in case of the turbine, in order to achieve a lower deviation of thermodynamic power with test data over the full range of speed. However, since this will defeat the purpose of the model to be used as a design tool for improvement of the turbine geometry, in this study only the value of BND that provided least deviation (5%) at the design speed was used. At off-design speeds, the deviation is high and consistently increases as speed decreases.

Mesh refinement study
The grid independence study of the BNL nozzle (run 309 presented in Table 1) indicated that grid refinement tends to produce an improvement in flow prediction through the nozzle. Similar grid independence exercise was performed on the base turbine case at 1560 and 3050 rpm. This was also required to test the dependence on rotational speed, which was absent   in the BNL nozzle flows. Based on the analysis presented in Section 4.2, a BND of 5e +07 m −3 and the Wolfert correlation for interphase heat transfer were prescribed in this mesh refinement study. Mesh sizes and boundary layer statistics are reported for the turbine model in Table 5.
In this case, y + has been reported for 3050 rpm flow condition. Four mesh refinement levels were examined, named as mesh 1 to mesh 4. Mesh 1 is the reference mesh. Meshes 2, 3, and 4 have an increasing refinement in the curved nozzle domain by factors for 1.57, 4.61, and 14.3, respectively, over the base mesh 1 size. Mesh 1 has a node count of 0.1e+6, which was highly refined to 1.6e+6 in mesh 4. Meshes outside the curved nozzle, in the flash tank domain, and in the inlet pipe were not altered between these meshes, as the curved nozzle mesh sensitivity was required to be examined. Curved nozzle mesh refinement levels with cross-section grid factors at the throat area are shown in Fig. 10 for the four meshes.
The parameters for evaluation of mesh refinement were the feed water flow rate and the turbine power estimated at both the speeds. Figure 11a shows a comparison of feed water flow rate between measurement and CFD estimate for different mesh refinement levels. At 1560 rpm, the deviation in flow is in the range of 7.2-7.7%. At 3050 rpm, the deviation reduces and is in the range of 3.3-4.5%. At both the speeds, the deviation does not vary when the mesh is refined from mesh 1 to mesh 4, indicating that mesh 1 was sufficiently refined for consistency of flow rate estimate. Figure 11b shows a comparison of shaft power between measurement and CFD estimate for different mesh refinement levels. At 1560 rpm, the deviation in shaft power is in the range of 47-54%. At 3050 rpm, the deviation reduces and is in the range of 16-27%. Similarly to the feed water flow rate, at both the speeds the deviation in shaft power does not vary when the mesh is refined from mesh 1 to mesh 4. Thus, mesh 1 was sufficiently refined for consistency of flow rate as well as shaft power estimate at both the speeds.    Figure 12 presents the evolution of pressure along the axis of the curved nozzle in the turbine mid-plane. The inlet boundary to the nozzle was specified as 400 kPa. From inlet to the nozzle throat, a steady increase in the fluid pressure is observed. Even though this is a converging section of the nozzle, the pressure rises due to centrifugal force acting on the fluid. The peak pressure upstream of the nozzle goes up to 670 kPa and then severely drops down close to 45 kPa downstream of the nozzle throat. Using this information, the nozzle section could be redesigned in order to avoid the pressure rise, which could improve the specific power output of the turbine. Figure 12 also presents the variation of vapour mass fraction (equation 21) along the axis in the mid-plane. In the converging section, the pressure is high, liquid is subcooled, and hence flashing does not occur. Feed water remains in liquid state with no phase transition. Just downstream of the throat, a severe drop in local pressure is observed to ∼45 kPa and this initiates the flash boiling as local saturation temperature is lower than the liquid temperature. The vapour distribution in this region is observed only close to the nozzle wall and not in the core of the flow. Here (at 0.06 m), the vapour mass fraction is 0%. After the throat (close to 0.08 m), the inception of water vapour due to flash boiling is observed. Further downstream due to rotation of the nozzle and the curvature, the flash boiling continues on the trailing surface of the nozzle. At nozzle exit, the mass fraction is ∼7% at the centreline. Figure 13 presents the distribution of vapour dryness fraction in the turbine nozzle at 4623 rpm. Vapour dryness fraction is defined by equation (21). Flashing initiates just after the throat area, and although homogeneous nucleation is set in the CFD model, vapour generation is seen to be high, adjacent to the throat wall.

Liquid and vapour distribution
x = m q (m p + m q ) kg/kg, m p = liquid mass, m q = vapour mass.
This effect is due to liquid velocity being relatively high in the core of the flow compared to that at the throat wall, resulting in higher residence time at the throat wall for interphase heat and mass transfer to develop. From the throat region and in the diverging section of the nozzle, there is steady increase in vapour fraction. Downstream, the flash development is highly nonuniform and vapour mass fraction is high on the low-pressure surface, i.e. trailing surface of the nozzle. At the nozzle exit, average dryness fraction is 12%, which is close to the test data of 12.8%. Heavier liquid phase is present on the high-pressure surface, i.e. leading surface. This nature of distribution of vapour mass fraction and flashing regime is highly dependent on the turbine rotational speed, though the exit vapour dryness fraction was found to be in the range of 11-12% over the entire speed range. Figure 14 is a plot of the liquid, vapour, and local saturation temperatures along the nozzle axis in the mid-plane. For the vapour phase, energy equation (3) is not explicitly computed and the solver is set to model vapour temperature at local saturation temperature. This practice is valid in the two-phase regime where vapour is not expected to get into superheated state at any point. Accordingly, the variation of vapour temperature closely follows the local saturation temperature distribution as seen in Fig. 14. Liquid at inlet is in subcooled state at 117 • C, while the saturation temperature at this pressure of 400 kPa is close to 143.6 • C. Due to high subcooled state, the phase change does not initiate even though homogeneous nucleation has been defined throughout the nozzle channel. From inlet to the upstream of the throat, as centrifugal force tends to increase the liquid pressure, the local saturation temperature also increases to exceed 163 • C. In this region, the liquid temperature does not vary. Just after the throat, flashing initiates and drop in local pressure results in a high degree of superheat in the liquid. The difference between the liquid temperature and saturation temperature at any point in Fig. 14 is the degree of superheat. The generated vapour is at local saturation temperature. From the nozzle throat up to the nozzle exit, it is seen that liquid temperature is dropping. This drop in temperature is due to the loss of latent heat in equation (3) transferred to the vapour phase. However, still the liquid remains at a superheated state at the exit, indicating that further flashing will continue after the nozzle exit, in the flash tank. This domain is not modelled in the CFD set-up and could be one of the contributors to higher vapour dryness fraction being reported in the measurements as compared to CFD model estimate. Figure 5 shows the reference used for defining angular span of the curved nozzle. The inlet is at 0 • , throat at 120 • , and the exit is at a 240 • span. A specific torque variation is extracted by scanning the nozzle span with 1 • interval as plotted in Fig. 15 at 4623 rpm. The sum of the specific torque over a certain span of the nozzle is plotted as the cumulative torque variation in Fig. 15. Cumulative torque at any specific channel span indicates the net torque output if the nozzle was clipped to that length. From inlet to ∼30 • , the rotor torque is positive and increasing. After this, the nozzle curvature and surface area variation is such that the specific torque becomes negative. This results in a decrease in the net torque, up to the throat section. After the throat section, the specific torque again becomes positive and the cumulative torque starts to increase up to the nozzle exit. The base turbine design has a tangential exit section and pressure acting on the   trailing surface results in a negative specific torque. This further reduces the cumulative torque and net torque output from this turbine. The nature of torque variation presented in Fig. 15 is observed at all speeds of the turbine. Torque analysis provides information on redesign of the inlet and exit of the turbine to eliminate negative specific torque so that net power output and specific power of the turbine can be improved.

Flash generation and isentropic efficiency
An expansion process that follows isentropic path results in a maximum enthalpy drop and hence highest power could be potentially extracted. However, this process also results in least flashing and least vapour dryness fraction at the end of expansion, which reduces potential for freshwater generation. On the other hand, an isenthalpic expansion will result in highest vapour formation but no power could be extracted. In reality, the expansion process will follow a path somewhere in between the isentropic and isenthalpic processes. Figure 16 presents the variation of measured and CFD predicted vapour dryness fraction with rotational speed. It is seen that the CFD expansion process is approaching towards an isentropic expansion as the vapour dryness fraction at turbine exit is close to isentropic dryness fraction.
100%, h i = enthalpy at inlet, h o = enthalpy at outlet. Figure 16 also presents the variation of isentropic efficiency with rotational speed. Isentropic efficiency is calculated using equation (22). The CFD model overpredicted power at a lower speed. As a result, the CFD isentropic efficiency is higher than the measured data at these speeds. The influence of speed on the expansion process needs to be further investigated in terms of the BND specification and thermal boundary conditions.

Conclusions
A two-phase CFD model for calculation of flash boiling flows was applied using the thermal phase-change formulation in ANSYS CFX solver to nozzles. The BNL nozzle experimental data and other CFD literature results on the same test case were first used in the model set-up to validate the formulation. The two-phase CFD model so established was then adapted with the moving reference frame, periodic boundary conditions, etc., to be used for the design and development of geothermal reaction turbine with curved nozzle channels.
The CFD model of the two-phase turbine has been analysed with respect to the specification of heat transfer parameters of the Nusselt number correlation and prescribed vapour BND. Flash boiling that has been incorporated through the thermal phase-change approach relies on local thermal nonequilibrium with homogeneous nucleation. The influence of Nusselt number based on the Ranz-Marshall and Wolfert types of correlations was studied and the Wolfert correlation was found to be more consistent when BND was varied. A BND of 5e +07 m −3 was found to be more suitable for the vapour regime formed in the turbine nozzle at the design speed.
The test turbine results of flow and power were evaluated against available test data at 400 kPa feed water pressure under subcooled condition of 117 • C and a very low backpressure of 6 kPa. Flow through this turbine was evaluated within 8% deviation. An overestimate in thermodynamic power by 30-50% was predicted at speeds <4000 rpm, while at the maximum speed of 4623 rpm the deviation was <5%. Rotor torque and hence power estimate was found to be dependent on the prescription of BND and heat transfer correlation in the CFD model. The vapour dryness fraction at turbine exit was close to an isentropic expansion vapour quality. The estimated isentropic efficiency was in the range of 7.5-17%. The turbine channel will be redesigned using the implemented two-phase turbine model to achieve power output with reasonable efficiency at atmospheric backpressure condition.