Numerical investigation of high pressure condensing flows in supersonic nozzles

High-pressure non-equilibrium condensing flows are investigated in this paper through a quasi-1D Euler model coupled to the method of moments for the physical characterization of the dispersed phase. Two different numerical approaches, namely the so-called (a) the mixture and (b) continuum phase model, are compared in terms of computational efficiency and accuracy. The results are verified against experimental data of high-speed condensing steam measured at high pressure (100.7 bar). The analysis demonstrates that Model (b) markedly outperforms the mixture model in terms of computational cost, while retaining comparable accuracy. However, both models, in their original formulation, lead to considerable deviations in the nucleation onset prediction as well as an overestimation of the average droplet radius. A further investigation is then conducted to figure out the main physical parameters affecting the condensation process, i.e. the surface tension, the growth rate and the nucleation rate. It is eventually inferred that applying proper correction to these three quantities allows to obtain best fit with the experimental data. A final calculation is carried out to show the dependence of these three correcting factors from the thermodynamic conditions of the mixture.


Introduction
High-speed flows usually condense at non-equilibrium conditions, namely when the vapor reaches thermodynamic states below saturation without any actual formation of liquid droplets. The physical characterization of such mixtures is extremely relevant for a vast number of industrial applications. For example, non-equilibrium condensation is responsible of significant aerodynamic efficiency decay in steam turbines [1]. Recently, several studies have been devoted to investigate supercritical CO2 condensation [2] for applications in compressor stages (carbon capture and sequestration, enhanced oil recovery, power cycles). Moreover, rapid expansions from supercritical states through the liquid-vapor dome can be exploited to enhance the efficiency of Organic Rankine Cycles [3].
Especially for application at high reduced pressures, two-phase flow models for condensation have to meet several requirements. In particular, a complex equation of state must be employed to account for the vapor non ideality. Note that the classical nucleation theory makes use of the perfect gas assumption for the derivation of the nucleation and the growth rate. Additionally, the reliability of the transport models, and in particular of the surface tension, is very questionable when approaching non-ideal thermodynamic regions and this strongly affects the physics of 1 condensation. On another note, if the models are employed for design purposes, the simulations must be also computationally efficient. To the authors' knowledge, there is currently a lack of models complying with all these requirements.
Multiple models and numerical schemes potentially applicable to metastable condensation are described in detail in literature [4]. However, the great majority of them have been validated with steam-water experiments at low pressure [4,5,6]. It should be noted that the so called "high-pressure" validation is traditionally made with the experimental data by Bakhtar [7], and even in this case P 0 is lower than 35 bar, with a corresponding reduced pressure of only 0.15.
The aim of the present work is to investigate non-equilibrium condensing flows in convergingdiverging nozzles at pressures much higher than those traditionally considered for steam-water. Since the full droplet spectra is not of interest in this study, the prediction of the second phase is performed through the well-established method of moments [8] in terms of averaged properties. Furthermore, only homogeneous condensation is of interest in this study, as the presence of foreign contaminants does not commonly trigger heterogeneous condensation in high-speed flows within supersonic nozzles.
The numerical model [9] consists of a quasi-1D Euler solver supplemented by two additional transport equations for the description of the second phase. Two different formulations are used for the main flow, i.e. the so called mixture and continuum phase conservation laws [9]. The results are compared with the experimental observations by [1].
The paper is organised as follows: the first section describes the governing equations of the model. The second section presents the numerical scheme employed for their resolution. The third section shows the comparison with the experimental observations and debates the results. Finally, the last section illustrates three different dimensionless corrections to the original model aimed at increasing its accuracy at high-pressures.

Governing equations
The characterization of a quasi-1D condensing flow requires a minimum of 8 variables: two thermodynamic properties for each phase, the velocities, the liquid mass fraction and the average droplet radius. As presented in [9], in order to reduce the number of equations down to 5 the following assumptions are made: i) the phases are in mechanical equilibrium, ii) the vapor velocity is equal to the liquid one and iii) the droplets temperature is determined through a simplified capillarity model. These assumptions are largely employed in literature [5,6] and proved to be adequate for low pressure steam test cases.
The model adopted in this work consists of three Eulerian conservation laws and two additional transport equations taken from the moments theory [9]. All relations are reported in the next paragraphs.

Mass, momentum and energy balance
This work carries on the comparison between the two main approaches described in [9], i.e. the employment of the mixture conservation laws and the continuum phase balance equations. The paper refers to these models as (a) and (b) respectively hereafter. The Eulerian set of equations where ρ, v, P , e 0 , h 0 are the density, velocity, pressure, total energy and total enthalpy respectively, A c is the cross sectional area of the nozzle and S v is the source term accounting for mass exchange across the interface vapor-liquid. The work refers to the properties of the continuum phase and the mixture with the pedex c and m respectively. The interested reader can refer to [5] for the derivation of such equations. On the other hand, the relations for Model (a) are obtained from the previous system using mixture properties and setting the source term S v equal to 0. The great restriction of Model (b) is that the equations are derived assuming a negligible volume fraction of liquid [5]. The last assumption can be adequate for low-pressure steam expansions, due to the very high difference in specific volume between the liquid and the vapor phase. On the other hand, the densities of the two phases get closer when moving towards the critical point, and the approximation on the droplet volume can lead to inaccuracies in the final solution. The majority of literature [6,8] adopts Model (a) as the validity of the mixture equations is general and fluid-independent.
However, the study in [9] proved Model (b) to be up to six times more efficient in terms of computational time. The main reason is that Model (a) requires an additional iterative procedure to determine the mixture properties, significantly increasing the simulation cost [6,9]. This work carries out a second analysis on a high pressure test case to choose the optimum solver in terms of both accuracy and code efficiency.

Moments equations for the disperse phase
The formulation adopted by the majority of the literature was presented by Hill [8], and involves 4 moments equations instead of the two theoretically required. Even though this choice leads to an increase in the computational cost, [8] claims that such model can reach a higher accuracy in the final results. However, the recent work in [9] presented a new 2-equation approach that proved to give comparable results with respect to Hill's model for low pressure steam. Therefore, the present work adopts this last solution and tests the accuracy of such numerical approach at high-pressure. The additional equations employed are where µ 0 and µ 3 are the moments of order -0 and -3, defined as in which f is the radial derivative of the droplet number density [8], N is the number of droplets per unit of mixture mass, R and R * the average and the critical radius respectively, J the nucleation rate and G the growth rate.

Closure models
All the missing terms in Eq.(2) and (3) must be modelled in order to close the equation system. In particular, the nucleation rate J and the growth rate G are modelled as in [9]. The empirical parameter β that enters in the model for G is here set to 1. Finally, the critical radius is evaluated as where ∆G is free Gibbs energy variation due to the condensation process.

Thermodynamic model
The required thermodynamic properties for the vapor phase are retrieved by an in-house thermophysical library [10]. On the contrary, the liquid properties, i.e. density and enthalpy, are directly implemented following the models in [11]. For the sake of accuracy the Gibbs free energy variation also includes the liquid contribution. ∆G is defined as the sum of two contributions To maintain a fluid-independent approach, ∆G v is evaluated through the same commercial thermodynamic library adopted for the continuum phase, whereas ∆G l is found as in [11]. The speed of sound for the mixture is evaluated as in [9]. Furthermore, the final mixture properties are found by means of the quasi-Newton algorithm. Finally, the surface tension is evaluated as in [12]: this model was compared with experimental measurements in equilibrium conditions up to the critical temperature, showing a maximum deviation of 0.4 · 10 −3 N/m.

Numerical resolution
The conservation laws for both Model (a) and (b) are discretized by means of finite volumes schemes, whereas the moments equations are rewritten as in which where x inlet and x outlet are the inlet and outlet abscissae of the domain and ∆x is the cell dimension. The system (7) is discretized via the upwind numerical scheme in [9]. The code adopts a segregated approach, i.e. the droplets properties are kept as constant during the resolution of the Eulerian equations. Implicit time integration is adopted to reach convergence.

Validation
The model was already validated for low pressure by using the well-established Moore nozzle A test case. For the sake of brevity, the work refers to [9] for a more detailed analysis. Results show a well agreement with regard to static pressure trend along the nozzle as well as a deviation from the experimental diameter of less than 1%. In order to validate the model at high pressure, the present study considers two steam expansions from [1], namely the tests No.18C and 18B. The next paragraph reports the comparison with the measurements.

High pressure investigation
The first test-case considered is No.18C, with total inlet conditions P 0 = 100.7bar and T 0 = 615.2K (P r = 0.46, T r = 0.95 ). The simulations were made on a 1000-cell mesh with implicit time integration. The expansion of highly superheated steam at same total pressure [1] is taken as reference. Fig.1a shows the expansion in the P-T chart obtained with Model (a). Due to the very high speed of the steam, the onset of condensation manifests for a temperature 15 K lower than the corresponding T sat at same pressure. A rapid change in the steam temperature and pressure is observed close to the nozzle throat after inception of condensation. The reason behind this behaviour is here explained: i) the number of drops rapidly increases from 0 to 3.2 · 10 16 kg −1 , ii)  a large amount of energy in the form of latent heat is released and, as a consequence, iii) the total temperature of the fluid rises. The simplified equation for ideal gas [13] states that the static pressure increment depends on the total temperature and the cross sectional area variation. In this case, due to the latent heat released, the dominant term is the first one, thus, the expansion rate decreases with respect to the one at the nozzle inlet. Note that, depending on the test-case, the release of latent heat can even lead to a considerable increase of the static pressure, as clearly shown for the low pressure test case in [9]. After the metastable condensation started, the release rate of the latent heat decreases, and the total temperature variation becomes lower. The fluid is already in the divergent section of the nozzle, thus, the cross sectional area A c is increasing. Therefore, the vapor continues the expansion and the static pressure rapidly reduces. Furthermore, the effect of metastabilities vanishes with time, and the vapor phase tends to reach equilibrium conditions: the steam approaches the saturation line at the nozzle outlet.
Given the assumptions of the model, the liquid phase is in mechanical equilibrium with the steam, thus, both phases have the same static pressure. Moreover, the temperature of the droplets is determined through a capillarity model as a function of the average radius R. The resulting thermodynamic state is close to the saturation line from the condensation onset till the nozzle outlet. Therefore, the droplets can be considered as saturated liquid. Fig.1b depicts the experimental pressure field along the nozzle. Fig.2 reports the final solution obtained in terms of pressure and average radius. Both Model (a) and (b) provide similar results in terms of static pressure distribution (Fig.2a). The solution at the nozzle inlet and outlet is close to the measurements. However, a considerable deviation is observed in the condensation onset with respect to experimental observation. The characteristic pressure bump caused by the latent heat is observed in correspondence of a dimensionless static pressure of around 0.62 (62 bar) instead of 0.42 (42 bar) as for the experimental data (Fig.1b). Fig.2b depicts the average droplets radius along the nozzle. The trend for Model (a) is similar to the experimental curve, however the results obtained are almost twice the nominal values measured. The prediction for Model (b) is very close to the one of the first model (the difference is around 1%), even though the radius increases faster along the nozzle.

Comparison between model (a) and (b)
In analogy with [9], a comparison is carried out between Model (a) and (b) in terms of computational cost and numerical stability. A first analysis is made with explicit time integration and CFL = 1. The simulations are stopped after a decrease of 4 orders of magnitude of the residual vector components. A single-phase simulation with highly superheated vapor at same total pressure is taken as benchmark. Two different tests are made changing the initial motion field, i.e. (i) constant solution along the nozzle or (ii) superheated solution. Table 1 shows the physical time required by the benchmark, Model (a) and Model (b) respectively. For the first case, the comparison shows that the simulation time for Model (b) is two times higher than the single-phase one. On the other hand, it is worth pointing out that Model (a) requires a computational time that is more than 15 times greater with respect to the benchmark. As already observed by [9], the quasi-Newton algorithm adopted represents a large penalty for Model (a). Even for the second test, in which the initial motion field is close to the final solution, the presence of such iterative procedure leads to a time that is more than 2 times higher than Model (b). To test the numerical stability of the two models, a second analysis on the maximum CFL allowable is made. The tests are made on the same 1000-cell mesh with implicit time integration and constant CFL. The initial motion field is again the constant solution along the nozzle. Table 2 shows the final result.
Note that the maximum CFL allowable for all the models is relatively low, as it is kept as constant during the simulation. The introduction of a time-dependent CFL allows for a maximum value of more than 30 for both Model (a) and (b).

Discussion
One of the key properties affecting the characteristics of metastable condensation is the surface tension σ. The underestimation of this parameter leads to a substantial anticipation of condensation onset resulting in fairly large difference in pressure and temperature trends. The work in [14] introduced a new parameter σ * defined as in which r σ is an empirical correction incorporating the effect of droplets curvature. It is reported that in the case r σ > 1 the condensation onset is delayed when compared to r σ = 1. However, the change in the nucleation starting point comes along with a considerable increase in the droplet radius. To gain insight of this behaviour, a simulation was performed with r σ = 1.5. Fig.3 shows the pressure field and the radius obtained. Additionally, Fig.4 depicts the trend for the nucleation and the growth rate. It can be observed that J decreases of one order of magnitude (Fig.4a): Eq.(5) shows that an increase in σ * reflects on the critical radius value. Furthermore, the nucleation rate has an exponential dependence from R 2 * and σ. As a consequence, J falls down from 1.3 · 10 24 kg −1 s −1 to 1.8 · 10 23 kg −1 s −1 , and the number of drops N significantly decreases. At the same time, the higher critical radius required for the nucleation causes a delay in the condensation onset, thus, an increase in the steam degree of subcooling ∆T sub . Due to the linear dependence of the growth rate from ∆T sub , G reaches a value 1.5 times higher than that one obtained with r σ = 1. Ultimately, the combined effect of the droplets number reduction and the increase of G leads to exceedingly high values of the average radius (Fig.3b).  Multiple tests were carried out on the test case No.18C to empirically determine the values of all the three coefficients for the best data fit ( Table 3). The analysis showed that the variation of the nucleation rate slightly affects the location of the condensation onset (columns (4) and (5) in Table 3). In particular, the reduction of J causes a nucleation delay, in agreement with the behaviour observed in the previous section. However, the influence of this correction factor is rather limited: the change of r J in the range 1 to 0.3 shifts the pressure bump from 0.44 to 0.42. Moreover, the corresponding radius variation is also minimal, from 5.7 · 10 −8 m to 5.5 · 10 −8 m (the value is taken in correspondence of the last experimental point in [1]). As observed in columns (3) and (4) in Table 3, the growth rate has a minor influence on the pressure field with respect to the surface tension. However, it considerably affects the average radius (from 13 · 10 −8 m to 5.7 · 10 −8 m). The best fit with the experimental data is obtained for r σ = 1.4, r J = 0.3, r G = 0.23. Fig.5 shows the achieved pressure field and average radius.  It is worth pointing out that the corrections are approximatively equal to 1 for the lowpressure test case. In particular, the final solution presents no deviation in the condensation  starting point, and a very high accuracy on the radius ( [9]) with respect to the measurements. Additionally, the droplet dimension is comparable for both the low and high pressure test case (5.0 · 10 −8 m and 5.5 · 10 −8 m as average radius respectively at the nozzle outlet). Therefore, it can be inferred that r σ , r J , r G are not only a function of the droplet curvature as in [14], but they arguably show a dependence from the fluid thermodynamic conditions. To support this hypothesis, a second simulation was carried out on the test No.18B in [1], characterized by the same total pressure (100.7 bar) but a different total temperature (638.53 K) (Fig.1a). Table 4 reports the final results. As condensation starts in the same thermodynamic region at a slightly lower pressure and temperature, it is expected to find values of r σ , r J , r G closer to 1 with respect to the previous case 18C. This hypothesis was confirmed by the results, which gave final values of r σ = 1.35, r J = 0.35, r G = 0.26. A further investigation is required to gain insight of i) the influence of pressure and temperature on these parameters and ii) whether a correction of σ, J, G could be sufficient to increase the accuracy of the numerical model closer to the critical point of fluids.  Table 4: Determination of r σ , r J , r G : condensation pressure and average radius, Test No.18B

Conclusions
In this paper, high-pressure condensing steam flows were simulated using the method of moments. A comparison in terms of accuracy, computational cost and numerical stability was carried out between two alternative numerical approaches proposed in [9], employing the mixture or the continuum phase conservation laws, i.e. Model (a) and (b) respectively. The study showed the superiority of Model (b) compared to Model (a) in terms of robustness and computational efficiency even for high-pressure condensation. Nonetheless, both models were found quite inaccurate in terms of static pressure trend, due to a considerable deviation in the location of condensation onset. Furthermore, the average radius was found almost two times larger than the measurements. New correction factors were then proposed and implemented in the models to improve accuracy. In particular, the work introduced three correction factors for the surface tension, the nucleation rate and the growth rate model r σ , r J , r G . It was observed that the rise of the surface tension significantly affects the final solution, causing a delay in the condensation onset as well as an increase of the average radius. Additionally, the droplet dimension are significantly affected by the growth rate G. On the other hand, J is the parameter with less influence, as shown in Table 3. The values of r σ , r J , r G were determined empirically to obtain the best fit with the experimental data, reaching almost the same accuracy as in the lowpressure test case in [9]. A second analysis was carried out on the so-called No.18B (Table 4) to confirm the dependence of r σ , r J and r G on the fluid thermodynamic conditions. Future efforts will be devoted to establish a physical correlation for these three parameters exploiting the data in [1] aiming at improving the reliability of the model for metastable condensation of steam (and possibly other fluids) close to the critical point, while maintaining high computational efficiency.