A COMPARISON OF THE EQUILIBRIUM AND THE DROPLETS BASED NON-EQUILIBRIUM COMPRESSIBLE PHASE CHANGE SOLVERS FOR CONDENSATION OF CARBON DIOXIDE INSIDE NOZZLES

In the current work, we simulate the condensation of supercritical CO2 during its high speed flow inside two different converging-diverging nozzles. We use the homogeneous equilibrium method and the classical nucleation theory based non-equilibrium phase change model for this purpose. The simulation results indicate significant influence of the nozzle inlet condition, nozzle shape and the fluid thermophysical behaviour on the nonequilibrium conditions prevailing inside the nozzles. We observe very low, ∼0.15 K, supercooling for the flow of CO2 inside the Claudio Lettieri nozzle compared to the supercooling of ∼3 K observed for the Berana nozzle. Very high nucleation rate (∼ 10 nucleation per m per second) is observed before the throat of the nozzles which remains confined to a very small axial distance. The nucleation rate takes much smaller values (∼ 10 nucleation per m per second) in rest of the nozzle. A maximum of 70 nano meter sized droplets with number densities of the order of 10 droplets per m are predicted inside the nozzles. Liquid mass fraction values between 0.2 to 0.4 are predicted by the solvers inside the nozzles. These results will be useful to the engineering community involved in the design and fabrication of CO2 based systems.


INTRODUCTION
There is recent surge in the use of CO2 for various industrial applications where it undergoes condensation during certain phases of the working cycle. Examples include the flow of CO2 inside compressors and turbines used in supercritical Brayton cycle (e.g., Lettieri et al. (2018); Lettieri et al. (2015); Ameli et al. (2018)) and the fracture induced depressurisation of pipelines used for carbon capture, transportation and storage (e.g., Brown et al. (2014); Nouri-Borujerdi and Ghazani (2017); Wen et al. (2019)). This phase change process is involved in other industrial applications as well which deal with high speed flow of vapors. The flow of steam through nozzles and turbine sections (e.g., Hill (1966); Moses and Stein (1978); Starzmann et al. (2018); Grübel et al. (2018); Xiao et al. (2017)) and the flow of refrigerants inside ejectors for pressure recovery (e.g., Palacz et al. (2015); Mazzelli et al. (2018); Giacomelli et al. (2019)) undergo similar phase transition. A common feature of such high speed flows with phase change is the existence of thermodynamic non-equilibrium conditions. The working fluids attain metastable thermodynamic states and suddenly undergo phase transition causing a condensation shock. The phase change process may be essential for the operation of the devices like condensing ejectors while having undesirable effects on the performance of steam turbines and supercritical CO2 compressors. The presence of shock waves, shear layers, boundary layer † Corresponding author. Email: ss.yadav@pilani.bits-pilani.ac.in separation, flow turbulence along with the phase change process make such flows difficult to analyse. There is a need to quantify the condensation behaviour with respect to the distribution of physical quantities like vapor supercooling, droplet nucleation rate, droplet diameter and number density, liquid mass fraction etc. Converging-diverging nozzles are simple devices with no moving parts and ideally suited for performing experimental or simulation based investigations of the phase change process. These are also used in most of the practical applications discussed before. Therefore, to quantify the phase change process of CO2, we simulate its flow inside converging-diverging nozzles. We use two different solvers available in Ansys CFX for this purpose, namely, the homogeneous equilibrium model based solver and the classical nucleation theory based non-equilibrium phase change solver. In the following, we briefly review the previous works related to the simulation of compressible phase change process focusing mainly on the works based on the homogeneous equilibrium and the non-equilibrium models.
The flow of CO2 through nozzles and radial compressors has gained attention recently due to its application in the supercritical Brayton cycle and in the carbon capture and storage based technique for atmospheric CO2 reduction. Such flows involve significant compressibility related effects and undergo phase change around leading edge of the the com-pressor blades. Here we summarize some of the most recent numerical works in this area. Ameli et al. (2017b) investigated the accuracy of various equations of state and the resolution of real gas property tables for flow inside supercritical radial compressors and turbines. For near critical inlet states, they observed that the Span and Wagner equation gives the most accurate results among the different equations used. Ameli et al. (2017a) compared the simulation based results on flow of supercritical CO2 through a converging-diverging nozzle, the experiments on which were provided by Berana et al. (2009). The thermophysical properties of CO2 in the metastable state were obtained using a bilinear extrapolation of the saturation state data. They observed that the flow becomes supersonic in the converging part of the nozzle just before the throat. They also highlighted the need to constrain the critical radius of the nucleating droplets to around a nano meter. Lettieri et al. (2018) experimentally investigated the metastable behaviour of CO2 during its high speed expansion inside a Laval nozzle. They obtained the Wilson line with the help of optical visualization and laser interferometry based techniques. It was observed that the computational data obtained with the extrapolated Span and Wagner equation of state agree within 2% of the experimentally observed values. Ameli et al. (2018) compared the numerical results on the flow of CO2 through a centrifugal compressor for different resolutions of the real gas property (RGP) table. They compared the CFD based distribution of entropy and density values over the compressor impeller with those obtained from NIST REFPROP (Lemmon et al. (2018)). They concluded that the errors in the simulated results decrease as the resolution of the RGP table increases but with increasing instability of the numerical solution procedure.
A significant amount of simulation based work has been performed on CO2 ejectors which are used for work recovery in supermarket refrigeration systems. Here we briefly discuss the important previous works related to ejector flow simulation. Smolka et al. (2013) expressed the heat conduction term in the energy equation in terms of enthalpy. This, combined with the pressure, was used to track the two phases under the homogeneous equilibrium approach. They observed a discrepancy of around 20% in the computed mass flow rates through the motive and the suction nozzles. Banasiak et al. (2014) analysed the entropy generation inside a two phase CO2 ejector based on the homogeneous equilibrium model. They concluded that the oblique shock train and flow turbulence contributed most towards the entropy generation inside the ejector's mixer and diffuser sections. Palacz et al. (2015) investigated the accuracy of homogeneous equilibrium method (HEM) for two phase flow through an ejector geometry. They concluded that the HEM model gives accurate results for operating conditions near and above the fluid critical point. The accuracy of the results decreased as the the inlet state was moved near to the saturation line. Palacz et al. (2017) compared the simulation data from the homogeneous equilibrium and homogeneous relaxation based phase change approaches with the experimental data on flow through an ejector. It was observed that both the phase change models predicted inaccurate mass flow rates for operating conditions far from the critical point. The HEM approach was found to be more accurate for operating conditions above the critical point. Using the homogeneous relaxation method, Haida et al. (2018) investigated the effect of relaxation time parameters on the distribution of pressure, velocity and vapor quality inside an ejector. An increase in the value of the coefficient θ0 decreased the vapor quality and the pressure inside the diverging portion of the motive nozzle. It also strongly influenced the shock formation and the flow process inside the mixing zone.
Recent simulation based works have considered the high speed condensation of vapor of one fluid inside another carrier gas. Using Ansys Fluent, Xiao et al. (2017) simulated the condensation of water vapor present in wet natural gas inside a nozzle. A maximum nucleation rate of 10 23 droplets per m 3 per second and droplet number density of 10 19 droplets per m 3 are observed inside the nozzle. Zheng et al. (2018) numerically investigated the condensation of CO2 during the high speed flow natural gas in a supersonic separator. A maximum nucleation rate of 10 21 droplets per m 3 per second and droplet diameter of 10 −7 m are observed inside the separator. Du and Hu (2020) investigated the effect of non-condensable gases on the phase change process of CO2 inside an evaporator. They theoretically showed that non-condensable gases convert the nucleation mechanism from homogeneous to heterogeneous and promote the phase change process.
From this brief literature survey, we observe that most of the simulation based works on condensation of CO2 have used the homogeneous equilibrium and the homogeneous relaxation based approaches. These approaches do not give an idea about the distribution of the nucleation rate, the number density and the diameter of the droplets. Although the classical nucleation theory has been widely applied for the condensation of steam, its use for the phase change of CO2 is more recent. In the current work, we simulate the phase change process of CO2 inside two different converging-diverging nozzles. We use two different compressible phase change solvers from the commercial flow solver Ansys CFX for this purpose. These are the homogeneous equilibrium based solver and the classical nucleation theory based non-equilibrium phase change solver. Based on the simulations, we predict the distribution of following physical quantities inside the nozzles: the supercooling levels attained by CO2; the nucleation rate; the droplet diameter and number density; and finally, the liquid mass fraction. This work is organized as follows: in Section 2, we describe the computational domains and the boundary conditions; in Section 3, we present the equations which govern the phase change of CO2 under the homogeneous equilibrium and the classical nucleation theory based non-equilibrium solvers; in Section 4, we present the details on the solver settings in Ansys CFX and the CO2 thermophysical properties; in Section 5, we discuss the results obtained by us using the simulation setup; and finally in Section 6, we conclude this work.

COMPUTATIONAL DOMAINS AND BOUNDARY CONDITIONS
In this section, we discuss the computational domains and the mesh related details for the two different converging-diverging nozzles used by us. The very first case considers the flow of CO2 inside a convergingdiverging nozzle used by Lettieri et al. (2018), the computational domain for which is shown in Fig. 1(a). A structured, body fitted, curvilinear mesh is generated inside the domain with smaller height of mesh elements near the walls. A total of 284739 hexahedral cells are generated in the domain using GMSH which is an open source mesh generator developed by Geuzaine and Remacle (2009). The wall mesh thickness is gradually varied normal to the walls, the non-dimensional wall distance y + varies from 8 near the throat region to 4 in the diverging potion under fully developed flow conditions. The domain shown in Fig. 1(b) investigates the flow of CO2 through a converging-diverging nozzle having straight walls and experimentally studied by Nakagawa et al. (2008), Berana et al. (2009). The mesh for this case has a total of 360477 hexahedral cells. The non-dimensional wall distance y + in this case varies from a maximum of 14 near the throat region to a minimum of 2 in the diverging portion of the nozzle under fully developed flow situation. In order to reduce the computational effort, we simulate the flow inside one fourth of the actual physical domains by using symmetry boundary conditions as shown in the figures. The various initial and boundary conditions at the inlet and outlet of the domains are shown in Table (1).

GOVERNING EQUATIONS
In this section, we describe the equations which govern the phase change process of a vapor under the homogeneous equilibrium and the classical nucleation theory based non-equilibrium models. The following discussion is based on Gerber (2008) and the Ansys CFX theory guide (ANSYS (2017)).

Equilibrium phase change solver
The homogeneous model assumes that the inter-phase transport processes occur at a much faster time scale compared to the fluid flow time scale so that there exists a mechanical, thermal and chemical potential equilibrium between the phases. The consequences of this assumption are the equality of the velocity, pressure, temperature and the Gibbs free enthalpy of the phases in the grid cells containing both the phases. The two phases are treated as a single fluid mixture with different thermophysical properties depending on the spatial distribution of the individual phases. Because of the mechanical and the thermal equilibrium between the phases, the homogeneous equilibrium model does not require an explicit formulation of the interfacial transport processes. In the following, the subscript l stands for the liquid phase, v stand for the vapor phase and m stands for their mixture. The mass conservation equation for the two-phase mixture is given by ∂ρm ∂t where ρm = α l ρ l + αgρg is the mixture density and α is the volume fraction of the phases with the constraint α l + αg = 1. The momentum conservation equation for the two phase mixture is given by where µm = α l µ l + αvµv is the viscosity of the mixture and S holds any other relevant momentum source term like that due to gravity. Ansys CFX solves for two different 'total enthalpy' based equations instead of solving a single mixture enthalpy equation. The equations for the total enthalpy of the liquid and the vapor phase are respectively given by In the above equations, Q l and Qg denote the heat transfer to the respective phases across the interface. The mass fraction of the vapor is calculated based on the enthalpy of the mixture using the following equation: This implies that a grid cell contains sub-cooled liquid if x < 0, superheated vapor if x > 1 and the two phase mixture if 0 ≤ x ≤ 1. The different thermophysical properties of the fluid are calculated based on the vapor mass fraction value.

Droplets based non-equilibrium phase change solver
The droplets based non-equilibrium phase change solver computes for the homogeneous nucleation and growth of droplets in a rapidly expanding fluid vapor. This model calculates the droplet diameter and the droplet number density as part of the solution. The mathematical formulation of the model involves equations for the continuous and the dispersed phase which are assumed to move at the same speed. In the following, the subscript c stands for continuous phase which is the vapor phase in our case while the subscript d stands for dispersed phase which is the condensed liquid phase in the form of droplets. The mass conservation equation for the continuous phase is given by while that for the dispersed phase is given by In the above equations, Sm represents the interfacial mass transfer term, J d is the critical nucleation rate andm d is the nucleated droplet mass based on the critical radiusR d . The constraint αc + α d = 1 applies for the volume fractions of the two phases. The analysis by Bakhtar et al. (2005) provides the following expression for the critical nucleation rate: Here qc is an empirical constant which is taken equal to one in this work, mm is the mass of a single molecule of CO2,Ĝb is the Gibbs number at the critical drop radius and is given bŷ In this equation, kB is the Boltzmann constant, ∆Ĝ is the change in the Gibbs free enthalpy of the supercooled vapor during the formation of a droplet of critical radius,R d , which are respectively given by . (11) In the above equations,R is the gas constant for CO2 vapor, Psat(Tc) is the saturation pressure at the vapor temperature, σ(Tc) is the surface tension coefficient which is a function of the continuous phase temperature, the expression for which is given later. Therefore, we get the expression for the Gibbs number at critical drop radius aŝ Further details on the theory and application in the context of steam condensation can be found in Bakhtar et al. (2005), Starzmann et al. (2018) and Grübel et al. (2018). The surface tension coefficient σ between the liquid and the vapor states of CO2 is a function of temperature and is implemented in the current work based on the following equation from Miqueu et al. (2000) σ(T ) = kBTcr NA Vcr 2/3 × (4.35 + 4.14ω)T 1.26 In this equation, NA is the Avogadro number, Vcr = 0.00214 m 3 /kg is the critical volume, ω = 0.239 is the acentric factor, Tr = 1 − T Tcr with Tcr = 304.13 K being the critical temperature of CO2. The variation of the surface tension coefficient with temperature based on Eq. (13) and the actual experimental data are shown in Fig. (2). Equation (13) is implemented as a CEL expression in Ansys CFX based on its reference guide, Ansys (2017).
The momentum conservation equation is similar to that for the homogeneous equilibrium model and is given by The continuous phase energy equation is given by while that for the dispersed phase is given by The equation governing the number density of the droplets is given by where the droplets move at the mixture velocity because of the no slip condition assumed between the phases. The droplet radius growth rate is given by where ξ is an empirical constant and Kn is the Knudsen number which is included to take into account the droplet size variation from the molecular to the continuum scales. For very small droplets (d ≤ 1µm), the dispersed phase energy equation, Eq. (16), is discarded in favour of the droplet temperature, T d , which is calculated using In this equation, Tsat(p) is the saturation temperature as a function of pressure p, Tsc is the supercooling level in vapor phase, R d is the droplet radius, and finally, R * d is the critical droplet radius. The term in Eqs. (6) and (7), related to the droplet growth due to mass transfer across the interface, is given by where β d = 3α d R d is the interfacial area available per unit volume. The source term in the equations for the total enthalpy of the phases, Eqs. (15), (16), is given by where ht,u is the upwinded total enthalpy whose value depends on the direction of mass transfer across the interface. The very small droplets (R d ≤ 1µ) are assumed to have a uniform temperature with all the heat transfer q d occurring in the continuous phase during evaporation and condensation with the heat flux calculated based on In the following section, we describe the solver settings used for the two compressible phase change solvers in Ansys CFX.

SOLVER SETTINGS AND THERMOPHYSICAL PROPERTIES
In this section, we briefly describe the various options used for the spatial and temporal discretization in Ansys CFX. The spatial discretization in Ansys CFX is based on control volume based finite element method under which the flow variables are stored at mesh nodes. The variation of the flow variables along the cell faces is given by geometric or parametric shape functions. The pressure and velocity are interpolated using a trilinear scheme while the pressure-velocity coupling is enforced using high resolution Rhie-Chow type scheme. The convective part of the mass, momentum and energy equations is discretized using high resolution scheme. The first order accurate implicit Euler method is used for temporal discretization. The volume fraction equation is solved in a coupled manner with the mass, momentum and energy equations with a residual and conservation target of 10 −4 . The vapor phase heat transfer is based on the total energy formulation. The surface tension coefficient is  calculated based on the vapor temperature. The particles based homogeneous nucleation model is used for the phase change process. The liquid phase heat transfer is calculated based on the small droplets assumption. The non-equilibrium solver requires fluid thermophysical properties in the subcooled state. This necessitates the use of a property database or an equation of state for the vapor phase amenable to extrapolation in the subcooled region. As per the Ansys solver theory guide (ANSYS (2017)), the IAPWS database for water and the Redlich-Kwong family of equation of state for CO2 are capable of such extrapolations but the latter's accuracy is unreliable near the critical point and for subcooled liquid states. For the equilibrium solver, we used the real gas properties of CO2 based on NIST REFPROP (Lemmon et al. (2018)). The temperature and pressure bounds for the RGP table are given in Table (2). The Soave-Redlick-Kwong equation of state, Soave (1972), is used for the thermophysical data required for non-equilibrium calculations, the table bounds for which are shown in Table (2). One has to use the expert parameter 'realeos liquid prop = 2' in Ansys CFX to force the solver to read the subcooled liquid properties based on the cubic equation of state. The parameters given in the Table ( For the non-equilibrium solver, a converged solution is obtained easily by adopting the following approach. First, the nucleation process is switched-off in the solver, allowing time steps of ∼ 10 −5 second to be used for attaining vapor supercooling in the domain. Afterwards, the nucleation is switched-on allowing phase change to be simulated at time steps of the order of 10 −9 second. Lastly, as mentioned in Ameli et al. (2017a), the critical radius of the droplets is calculated in Ansys CFX based on the Gibbs free enthalpy change of the continuous phase (∆Gc = g(Pc, Tc) − g(Psat, Tc)) using Near the critical point, the Gibbs free energy change becomes very small due to which the critical radius value become as high as 10 −3 m. In order to avoid such nonphysical values, an 'expert parameter' is used in the solver in order to constrain the critical radius to a maximum of 10 −8 m.
In the following section, we present the various results obtained by us based on the simulations performed with Ansys CFX.

RESULTS AND DISCUSSIONS
In this section, we discuss the various results obtained with Ansys CFX on the phase change process of CO2 inside the two converging-diverging nozzles. First we present the results on the flow of supercritical CO2 through a nozzle, the experiments on which were performed by Lettieri et al. (2018). We call this case as the Lettieri test case. Next we present the simulation results related to the flow of supercritical CO2 in a converging diverging nozzle with straight walls, the experiments on which were performed by Berana et al. (2009 Figure (3) shows the experimental pressure distribution along the center line of the Claudio Lettieri nozzle along with those based on the numerical simulations using the equilibrium and non-equilibrium phase change solvers in Ansys CFX. As can be seen, the equilibrium solver gives a better pressure distribution in the nozzle which may be due to the more accurate fluid thermophysical data coming from NIST REFPROP (Lemmon et al. (2018)). The non-equilibrium solver under-predicts the pressure in the diverging part of the nozzle which can be attributed mostly to the use of cubic equation of state for the thermophysical properties of CO2. Other parameters of the solver like the maximum critical radius, bulk nucleation tension factor, interphase transfer model may also be playing a role here. Figure (4) shows the distribution of the pressure and vapor temperature contours on the symmetry planes of the Lettieri nozzle based on the non-equilibrium solver. The fluid cools down as it expands inside the nozzle with a minimum temperature equal to 265 K at the exit of the nozzle. A similar distribution is obtained with the equilibrium solver as well. Figure (5) shows the distribution of the Mach number and the local speed of sound for the vapor phase of CO2 in the Claudio Lettieri nozzle. It can be seen that flow becomes supersonic only in the diverging portion of the nozzle while it is subsonic near the throat of the nozzle. This is in spite of the fact that the local speed of sound is lowest in Lettieri nozzle based on the equilibrium and the non-equilibrium solvers. It can be seen that the maximum liquid mass fraction at the nozzle outlet is around 0.3 based on the equilibrium solver while it is equal to 0.4 with the non-equilibrium solver. Fig. 7 Contours of the vapor supercooling and the nucleation rate inside the Lettieri nozzle based on the non-equilibrium solver. The nucleation rate is zero in the converging portion of the nozzle and becomes very high when the vapor becomes supercooled.
The non-equilibrium solver gives information on several other parameters which are discussed in the following. Figure (7) shows the distribution of the vapor supercooling and the nucleation rate on two mutually perpendicular, symmetry planes of the Lettieri nozzle. It can be noted that there exist very small (∼ 0.15 K) vapor supercooling inside the nozzle after the phase change process. This is due to the near critical condition (80 bar, 311 K) of CO2 at the inlet to the nozzle where the spinodal limit itself is near to the saturated vapor line. This probably hints at the low overall non-equilibrium conditions existing in the nozzle. The figure also shows that there is a very high droplet nucleation intensity (∼ 10 35 nucleation per m 3 per second) confined to a narrow portion of the nozzle length before the throat which falls steeply to a much smaller, almost constant value (∼ 10 7 nucleation per m 3 per second) in the diverging part of the nozzle. The line plot shown in Fig. (8) represent the same information along the Lettieri nozzle center line. It can be observed in the figure that the degree of liquid subcooling is also very small ∼ 0.15 K. Figures (9) and (10)    the droplet number density inside the Lettieri nozzle. It can be seen that the droplet size is the smallest near the nucleation zone which then increases steeply to around 65 nano meter in a very short axial distance. Similarly, starting from a negligible value at around 20 mm, the number density of the droplets increases steeply in a very small axial distance, and afterwards, decreases gradually along the nozzle length. Here we briefly discuss the results on the phase change process of CO2 inside the Berana nozzle, the experiments on which were performed by Berana et al. (2009). Figure (11) compares the simulation based pressure distribution in the Berana nozzle with the experimentally obtained values. It can be noticed that the equilibrium solver again predicts a better pressure distribution compared to that obtained with the non-equilibrium solver in spite of the better theoretical formulation of the latter. This can be attributed to two reasons. First, the equilibrium solver uses the REFPROP (Lemmon et al. (2018)) based more accurate thermophysical properties of CO2 in the form of RGP tables. Second, the non-equilibrium solver involves various empirical parameters in the theoretical formulation whose values are unknown before-hand. As can be seen in Fig. (11), the equilibrium solver predicts a shock wave located immediately downstream to the throat as revealed by the kink in the pressure distribution curve. Both the equilibrium and the non-equilibrium solvers under-predict the pressure values after the throat of the nozzle. It may be also be noted here that no experimental data is available for pressure distribution near the exit of the nozzle. Figure (12) shows the corresponding pressure contours along with the distribution of vapor temperature on the two symmetry planes inside the Berana nozzle. It should be noted that the vapor temperature increases as the throat is approached starting from a value of 318.15 K at the inlet boundary. After the throat, the temperature decreases along the diverging portion of the nozzle. Based on the non-equilibrium solver, Ameli et al. (2017a) simulated flow inside a two dimensional geometry of the Berana nozzle and observed that the flow becomes supersonic inside the converging portion of the nozzle, that is, before the throat of the nozzle. This is not the case Fig. 13 Contours of the local speed of sound and the Mach number inside the Nakagawa nozzle based on the non-equilibrium solver.
when a three dimensional computational domain is used as can be seen in Fig. (13) which shows a Mach number of around 0.65 at the throat. Information on the distribution of the liquid mass fraction of CO2 along the Berana nozzle center line can be gained from Fig. (14). The figure shows that the equilibrium solver predicts a higher value (∼ 0.2) of maximum liquid mass fraction distribution inside the nozzle compared to the nonequilibrium solver which gives a maximum value of around 0.16. Also, the liquid fraction starts increasing abruptly from a value of 0.1 under the equilibrium solver while it increases gradually from zero under the non-equilibrium solver. The liquid mass fraction decreases because of the viscous heating by the shock wave whose presence can be guessed from the pressure variation in Fig. (11).     Figure (15) shows the variation of supercooling / subcooling levels achieved by the phases of CO2 and the nucleation rate of the droplets along the center line of the Berana nozzle. In contrast to the Lettieri nozzle case, the supercooling level here is around 3 K which shows that non-equilibrium effects are stronger in this case. This may be due to the higher pressure and temperature values (90 bar, 318 K) at the inlet to the nozzle so that the state is farther from the critical point compared to that in the Lettieri case. However, it should be noted that the vapor supercooling and the liquid subcooling quickly becomes very small further downstream the nozzle showing that the flow equilibrates in a short distance (∼2 mm) after the throat. The nucleation rate starts with a very high value (10 35 nucleation per m 3 per second) and then falls steeply like in the Claudio Lettieri case. Afterwards, the nucleation process becomes almost constant further downstream in the nozzle with a value of around 10 6 nucleation per m 3 per second. Figure (16) shows the variation of drop diameter and their number density along the center line of the Berana nozzle . As can be seen, the droplet number density starts with a high value (∼ 10 20 droplets per m 3 ) at the throat followed by a sharp decrease before the kink. Afterwards, the droplet number density decreases gradual along the downstream of the nozzle. The size of the droplets increase gradually to 70 nano meter along the length of the nozzle starting from smaller values at the throat. It may be noted that the droplet diameter decreases while the droplet number density increases near the exit of the nozzle. This may be due to a shock wave existing in the diverging portion of the nozzle just before the outlet as can be seen from the Mach number contours in Fig. (13).

CONCLUSIONS
To summarize the current work, we compared two different compressible phase change solvers in Ansys CFX, namely, the equilibrium phase change solver and the classical nucleation theory based non-equilibrium solver for the phase change process of CO2 inside two different nozzles. The aim of the work was to gain an understanding about the capability and accuracy of the flow solvers for simulating the compressible, high speed flow of CO2 with phase change inside nozzles. Our experience is summarized as follows. The droplets based non-equilibrium solver needs thermodynamic data from the metastable states for simulating the nonequilibrium conditions inherent in such flows. The solver works well with cubic equations of state but requires fine tuning of several 'expert parameters' for obtaining correct results. The equilibrium solver is more robust and more accurate with respect to prediction of pressure inside the nozzles with REFPROP based real gas property data for CO2. Compared to experimental data, the solvers under-predicts the pressure distribution inside the nozzles. The simulation results indicate significant influence of the nozzle inlet condition, nozzle shape and fluid thermophysical behaviour on the non-equilibrium effects existing inside the nozzles. For example, we observe very low (∼ 0.15 K) supercooling for the flow of CO2 inside the Claudio Lettieri nozzle compared to the supercooling (∼ 3 K) observed for the Berana nozzle. The nucleation process inside the nozzles is observed to start abruptly with a very high rate (∼ 10 35 nucleation per m 3 per second) which remains confined to a very small axial distance and takes much smaller values (∼ 10 7 nucleation per m 3 per second) in rest of the nozzle. The non-equilibrium effects are also confined to a smaller portion of the nozzle and the vapor flow quickly equilibrates itself along the downstream of the nozzle. A maximum of 65 to 70 nano meter sized droplets and number densities of the order of 10 21 droplets per m 3 are predicted inside the nozzles. Liquid mass fraction values ranging from 0.2 to 0.4 are predicted inside the nozzles. In general, the solvers are able to simulate compressible, high speed phase change problems with acceptable accuracy which depends on several factors, among which, the use of real gas property database is the most important one. There is a need for a thorough investigation of the effect of different 'expert parameters' which control the behaviour of the non-equilibrium solver. The experience gained here is being used to study the flow of supercritical CO2 inside an ejector geometry.