Assessment of Cavitation Models for Compressible Flows Inside a Nozzle

This study assessed two cavitation models for compressible cavitating flows within a single hole nozzle. The models evaluated were SS (Schnerr and Sauer) and ZGB (Zwart-Gerber-Belamri) using realizable k-epsilon turbulent model, which was found to be the most appropriate model to use for this flow. The liquid compressibility was modeled using the Tait equation, and the vapor compressibility was modeled using the ideal gas law. Compressible flow simulation results showed that the SS model failed to capture the flow physics with a weak agreement with experimental data, while the ZGB model predicted the flow much better. Modeling vapor compressibility improved the distribution of the cavitating vapor across the nozzle with an increase in vapor volume compared to that of the incompressible assumption, particularly in the core region which resulted in a much better quantitative and qualitative agreement with the experimental data. The results also showed the prediction of a normal shockwave downstream of the cavitation region where the local flow transforms from supersonic to subsonic because of an increase in the local pressure.


Introduction
Cavitation in confined flows is a highly complex phenomenon observed in many applications and can affect their performance considerably; the occurrence of cavitation in fuel supply systems of aircraft and automobiles, cryogenic pipelines and pipelines of nuclear power plants can lead to unexpected degradation in the system performance and damage to the components. The bubbly cavitating fluid behaves like a compressible fluid and has an unusual characteristic; the presence of even a small fraction of gas/vapor can considerably reduce the local sonic speed of liquid-gas mixture compared to the sonic speed within its constituents [1][2][3][4]. This unusual phenomenon occurs because the mixture phase is easily compressed owing to the compressibility of gas bubbles, however, it remains dense due to the dominant mass of the liquid [5]. Therefore, the bubbly flows in nozzles frequently become locally supersonic which leads to choking and the development of local shockwaves that impairs the performance [6]. Hence, because of the importance of an understanding of internal flows in nozzle design, several studies were conducted on flow inside the orifice which revealed complex bubbly flow behavior as described above. Typical nozzle geometries were sharp-edged entrance orifices [6][7][8], converging-diverging nozzles [9][10][11], and I-Channel [12]. Besides, many theoretical studies have been conducted to realize this phenomenon. Dorofeeva et al. [13] proposed a 1D barotropic model for cavitation prediction in the converging-diverging nozzle. The model assumed the fluid to be a homogenous mixture of liquid and cavitation vapor; the liquid was considered to be incompressible and the cavitation vapor to be an ideal gas. The model predicted a sharp rise in pressure or "shock" just downstream of the cavitation clouds where vapor transforms to liquid because of higher local pressure.

Flow Configuration
Duke et al. [15] performed X-ray radiography of a model single-hole nozzle to measure the cavitation vapor distribution. The test geometry was a polycarbonate nozzle with a throat diameter of D = 0.5 mm, a length/diameter ratio (L/D) of 5, and inlet to nozzle diameter ratio of 5 as shown in Figure 1. The working fluid was commercial surrogate gasoline with a density equal to 777.97 kg/m 3 , the viscosity of 8.48 ×10 −4 kg ms , and the vapor pressure of 580 Pa at 25 • C. The working fluid was delivered via a piston accumulator system pressurized with an inert gas. The fluid temperature and pressure were monitored immediately upstream, and downstream of the nozzle and a turbine flow meter was used to monitor the flow rate. The injection pressure P inj and back pressure P back were varied independently to obtain the desired flow rate. Operating condition of the experiment is listed in Table 1, and the same operating condition was used for the CFD simulations in this study. delivered via a piston accumulator system pressurized with an inert gas. The fluid temperature and pressure were monitored immediately upstream, and downstream of the nozzle and a turbine flow meter was used to monitor the flow rate. The injection pressure and back pressure were varied independently to obtain the desired flow rate. Operating condition of the experiment is listed in Table 1, and the same operating condition was used for the CFD simulations in this study.  During the X-ray radiography experiment, the 7-BM beamline of the advanced photon source was with the focused beam acting as a microprobe, allowing the void fraction to be probed on the small area along the line of sight. They passed the nozzle through the fixed beam in a raster pattern at 100 transverse positions (y) and 19 streamwise positions (x), the raster scan grid can be seen in Figure 2a. The time-averaged values were interpolated onto a series of contour plots, as shown in Figure 2b which shows cavitation originating near the sharp entrance of the nozzle. At half of the nozzle length, the vapor layer attached to the wall begins to collapse while a considerable void cloud was recorded at the core of the nozzle which then extends to the nozzle outlet. It was argued [15] that the presence of void at the centerline is more likely due to convection of vapor from the wall because of negatively inward pressure gradient than due to isolated nucleation events as described by Bauer et al. [36]. Nonetheless, the time-averaged measurements of vapor quantity achieved an uncertainty of 2% implying a steady flow regime. Duke [15] has also performed imaging using conventional Xray, the result of which is shown in Figure 2c which is a time-averaged X-ray image of the nozzle inlet which also described the above-stated phenomenon. For more details about the experimental work, please refer to the original work [15].  During the X-ray radiography experiment, the 7-BM beamline of the advanced photon source was with the focused beam acting as a microprobe, allowing the void fraction to be probed on the small area along the line of sight. They passed the nozzle through the fixed beam in a raster pattern at 100 transverse positions (y) and 19 streamwise positions (x), the raster scan grid can be seen in Figure 2a. The time-averaged values were interpolated onto a series of contour plots, as shown in Figure 2b which shows cavitation originating near the sharp entrance of the nozzle. At half of the nozzle length, the vapor layer attached to the wall begins to collapse while a considerable void cloud was recorded at the core of the nozzle which then extends to the nozzle outlet. It was argued [15] that the presence of void at the centerline is more likely due to convection of vapor from the wall because of negatively inward pressure gradient than due to isolated nucleation events as described by Bauer et al. [36]. Nonetheless, the time-averaged measurements of vapor quantity achieved an uncertainty of 2% implying a steady flow regime. Duke [15] has also performed imaging using conventional X-ray, the result of which is shown in Figure 2c which is a time-averaged X-ray image of the nozzle inlet which also described the above-stated phenomenon. For more details about the experimental work, please refer to the original work [15].  [15]. (a) The raster scan grid where X-ray radiography measurements of void were taken, (b) time-averaged X-ray radiography measurement of the void at CN = 11.2, (c) time-integrated X-ray images of nozzle inlet at CN = 11.2. (Reproduced with permission).

Mixture Model
The mixture model solves one set of conservation equation for mass, momentum, and energy and the volume fraction equation for the second phase with the addition of the turbulence closure model. The model assumes that velocity, temperature, and pressure between the phases are equal. This assumption is based on the idea that the difference between three potential variables would facilitate mass, momentum, and energy transfer between phases so fast that equilibrium would be achieved within an extremely small time scale (almost instantly). Therefore, the resulting equations  [15]. (a) The raster scan grid where X-ray radiography measurements of void were taken, (b) time-averaged X-ray radiography measurement of the void at CN = 11.2, (c) time-integrated X-ray images of nozzle inlet at CN = 11.2. (Reproduced with permission).

Mixture Model
The mixture model solves one set of conservation equation for mass, momentum, and energy and the volume fraction equation for the second phase with the addition of the turbulence closure model. The model assumes that velocity, temperature, and pressure between the phases are equal. This assumption is based on the idea that the difference between three potential variables would facilitate mass, momentum, and energy transfer between phases so fast that equilibrium would be achieved within an extremely small time scale (almost instantly). Therefore, the resulting equations resemble a pseudo-fluid with mixture properties and the equation of state which links phases, so mixture thermodynamics properties are obtained [18,19].
The steady-state version of the model equations are presented here. The following formulation is used to calculate mixture density: where α k is the volume fraction of the phase k. Since in the present study, the multiphase system is comprised of two phases, hence, the mixture density is computed using the following equation: Since, Hence, Equation (2) becomes: Similarly, the viscosity of the mixture can also be calculated as: where µ l is the molecular viscosity of liquid and µ v is the molecular viscosity of the vapor. The continuity and momentum equations are as follows: where → v is the mixture velocity, p is the pressure, and = τ is the mixture shear stress tensor. Similar formalisms apply to energy conservation and turbulence models [17,18,21,26]; equations are excluded here for brevity.

Schnerr and Sauer (SS) Model
The SS model [23] treats the bubbly flow as a homogeneous vapor-liquid mixture. The model assumes the vapor phase to consist of numerous mini spherical bubbles. The model uses the Rayleigh relation [21] for the description of bubble growth and collapse. .
where R B is the radius of the bubble, p v is the saturation vapor pressure, and p ∞ is the far-field pressure, however for the practical purpose p ∞ is taken to be the same cell center pressure [26]. The model defines the vapor volume fraction using the following equation: where n B is defined as the bubble number density per unit volume of pure liquid. Using Equation (8), the following equation is derived for the bubble radius: The transport equation for the vapor fraction is: where R is the net phase change rate. The phase transfer rate terms are as follows: When p ∞ ≥ p v , condensation The vapor volume fraction is dependent on the n B which is a constant and assuming no change in the nucleation site per unit volume of pure liquid. The originally recommended value of n B is 10 13 /m 3 which is used in this study in initial simulations, voften n B is required to be adjusted to match the experimental data as shown by Yuan et al. [37].

Zwart, Gerber, and Belamri (ZGB)
The ZGB model [27] also assumes the fluid to be a homogeneous liquid-vapor mixture, and the model similarly assumes vapor to consist of tiny spherical bubbles. The ZGB model also uses the Rayleigh relation (Equation (8) for the description of bubbles growth and collapse in vapor production (cavitation) and destruction (condensation) terms [27]. The following equation defines the vapor volume fraction: where R B is again the bubble radius, N B is the bubble concentration per unit volume of the fluid mixture. The transport equation for the vapor fraction is formerly identical to the SS model as: where R is again the net phase change rate and is defined using the following equation: The phase transfer rate terms are as follows: When p ∞ ≥ p v , condensation Zwart et al. [27] reported that the model works well for condensation but produces physically incorrect results (and numerically unstable) when employed for cavitation. This is because the model assumes no interaction between the cavitation bubbles, which is only possible in the early stages of cavitation when cavitation bubbles grow from a nucleation site. However, they argued that when the vapor volume increases, the nucleation site density must decrease accordingly. Hence, Zwart et al. proposed to replace α v with α nuc (1 − α v ) in the evaporation term R e , where α nuc is the nucleation site volume fraction. Hence, final forms of rate equations are: When When p ∞ ≥ p v , condensation The values of model constants are, α nuc = 5 × 10 −4 , F vap = 50 and F cond = 0.01.

Tait Equation
Tait establishes a nonlinear relationship between the density of the liquid and its pressure under isothermal conditions. The Tait equation can be expressed in terms of pressure and density using the following relationship: where a and b are coefficients which can be defined by assuming that bulk modulus is a linear function of pressure, the values of a and b are based on the reference state values of pressure, density, and bulk modulus. n is known as density exponent, having a similar role as a ratio of specific heats [38]. The simplified Tait equation as implemented in the ANSYS Fluent platform has the following form [39]: where and where p 0 is the reference liquid pressure, ρ o is the reference liquid density which is the density at the reference pressure p 0 , n is density exponent, for which the value 7.15 is used which corresponds to weakly compressible materials such as liquids [38], K 0 is the reference bulk modulus at the reference pressure p 0 , p is the liquid pressure (absolute), ρ is the liquid density at the pressure p, and K is the bulk modulus at the pressure p. The reference values that are considered in this study are of gasoline at the average temperature of 25 • C [40] and listed in Table 2. The speed of sound in the liquid phase is calculated using the following equation:

Ideal Gas Law
The vapor density is modeled using the ideal gas law which has the following form: where R u is the universal gas constant, MW is the molecular weight, and T is the temperature. The speed of sound in the gas phase is calculated using the following equation: where γ is the specific heat ratio, and it is defined by the following equation: where R is the specific gas constant R u MW , c p is the specific heat capacity at the constant pressure, and c v is the specific heat capacity at the constant volume.

Wallis Sonic Speed in the Two-Phase System
The speed of sound in the mixture phase is computed using the Wallis Model [41], the model assumes local thermodynamic equilibrium between phases and has been implemented in the following form in ANSYS Fluent [42]: where a m is the speed of sound in the liquid-vapor mixture. The readers can find the complete derivation of Equation (29) in [43].

Numerical Methods
The governing equations are solved using commercial CFD solver ANSYS Fluent v14.5 which uses the finite volume method [30]. Acknowledging that the flow regime is steady [15], steady-state simulations have been performed, and for pressure-velocity coupling, the SIMPLE algorithm is employed. Since the nozzle geometry is axis-symmetrical, a 2D axis-symmetrical flow domain (as shown in Figure 3) has been chosen for the CFD simulations and cylindrical form of governing equations (Equations (6), (7), and (11)) are solved. The 2D grid contained quad type elements. The operating boundary conditions used to simulate the flow are the same as those used in the experiment as specified in Table 1. The grids of the main flow region (the red box shown in Figure 3) were locally refined as shown in Figure 4. The influence of the number of cells in grids was checked by comparing predicted vapor volume fraction and velocity as shown in Figure 5a,b, respectively. The first grid comprised 89,189 cells or an approximate average cell size of 16.87 microns, the second grid contained 121,526 cells and hence had an average cell size of about 14.45 microns, and the third grid had 273,563 cells with average cell size was approx.10.06 microns. The cell ratio between the first grid and the second grid was 1.16 and between the second and third grid was 1.43. The authors achieved a good compromise between second and third grids in terms of results of mean axial velocity and vapor volume fraction, as shown in Figure 5. However, after the third refinement, the authors did not achieve y+ < 1 which is required to run low-Re turbulence models [30]. Hence grids were refined only in the near-wall region to ensure y+ < 1 was achieved. Hence, there is not much difference between the third grid and the final grid in terms of the cell quantity.          Tables 1 and 2. The effect of convergence was checked by changing the convergence criteria of continuity, momentum, k, and epsilon from 1 × 10 −4 to 1 × 10 −6 which showed no observable changes on the solution.   In the grid independence test high Reynolds number standard k-epsilon model [32] was used, and near-wall turbulence was simulated using the EWT (enhanced wall treatment method), which is a near-wall modeling approach that blends linear (laminar) and logarithmic (turbulent) laws-of-the-wall using a joining function [44]. For cavitation modeling, the ZGB cavitation model has been used with a mixture multiphase model. The operating condition and physical properties of the working fluid are listed in Tables 1 and 2. The effect of convergence was checked by changing the convergence criteria of continuity, momentum, k, and epsilon from 1 × 10 −4 to 1 × 10 −6 which showed no observable changes on the solution.

Influence of Turbulence Modelling
In this section, impacts of four different turbulence models on cavitation vapor predictions are assessed. Turbulence models evaluated here are the Standard k-epsilon, Realizable k-epsilon, Launder and Sharma k-epsilon, and SST k-omega. The density of both vapor and liquid phases are assumed to be constant. This assumption is only justified if the volume fraction is low [4], but has advantages as it allows the volume fraction transport equation (see Equations (10) and (15)) to become just a volume conservation equation. Thus it helps solutions to become more numerically stable, and therefore enables researchers to assess models in the early stages of assessments, and hence, it has been used by many researchers like [27,39,[45][46][47]. Rigorous evaluations were carried out to select an appropriate turbulence model for this study with full details given by [48], and a summary of some of these simulation tests with comparison to experimental data is given below in Figures 6 and 7. In this section, impacts of four different turbulence models on cavitation vapor predictions are assessed. Turbulence models evaluated here are the Standard k-epsilon, Realizable k-epsilon, Launder and Sharma k-epsilon, and SST k-omega. The density of both vapor and liquid phases are assumed to be constant. This assumption is only justified if the volume fraction is low [4], but has advantages as it allows the volume fraction transport equation (see Equations (10) and (15)) to become just a volume conservation equation. Thus it helps solutions to become more numerically stable, and therefore enables researchers to assess models in the early stages of assessments, and hence, it has been used by many researchers like [27,39,[45][46][47]. Rigorous evaluations were carried out to select an appropriate turbulence model for this study with full details given by [48], and a summary of some of these simulation tests with comparison to experimental data is given below in Figures 6 and 7.
The total quantity of vapor integrated over the nozzle cross-section at different axial locations (x/L) for different models has been compared with experiment and is represented in Figure 6. Before discussing the quantitative comparison, it would be useful to know how the vapor is distributed within the nozzle. Observing the experimental X-ray radiography image, Figure 7a, shows that vapor is generated at the sharp entrance of the nozzle orifice, as expected. The vapor then disperses from near the wall and mixes with the flow as it travels further downstream. Also, the accumulation of vapor along the nozzle axis downstream of the entrance can be seen clearly. It was argued [15] that the observed asymmetric vapor distribution has been attributed to small machining defects which may have provided additional nucleation sites to the nozzle. On comparison of the predicted results using for different turbulence models shown in Figure  6, it is clear that the realizable k-epsilon and standard k-epsilon models outperformed the other two models, and that the realizable k-epsilon model predictions give closer agreement with the experimental data than the standard k-epsilon model, especially toward the exit of the nozzle. The results show the realizable k-epsilon model overpredicts the vapor fraction near the entrance and becomes similar as x/L increases, while the standard k-epsilon model underpredicts in the beginning of the nozzle and the difference becomes much wider with downstream distance. This can be correlated with the larger recirculation zone predicted using the realizable k-epsilon model at the beginning of the nozzle (x/L = 0.0167) as shown fully in [48] from the mean axial velocity profiles; the larger recirculation zone corresponds to higher flow acceleration into the nozzle, which in turn Figure 6. Comparison of predicted total vapor volume fraction, obtained from four turbulence models, with experimentally measured total vapor volume fraction across the nozzle at different axial positions along the nozzle. Error bars in the experimental data show uncertainty of around 2% [15].
reaches the nozzle axis, and that the cavitation structure extends beyond the exit of the nozzle, Figure  7c. These features can be seen better in Figure 8 where the cavitation contours for both models are plotted and clearly show the extension of cavitation into the circular expansion region of the orifice for the realizable k-epsilon model (Figure 8b) while it terminates at the exit of orifice for standard kepsilon, Figure 8b. It is also clear from Figures 7 and 8 that the intensity of the cavitation is higher with the realizable k-epsilon model in line with the results obtained in Figure 6.  The total quantity of vapor integrated over the nozzle cross-section at different axial locations (x/L) for different models has been compared with experiment and is represented in Figure 6. Before discussing the quantitative comparison, it would be useful to know how the vapor is distributed within the nozzle. Observing the experimental X-ray radiography image, Figure 7a, shows that vapor is generated at the sharp entrance of the nozzle orifice, as expected. The vapor then disperses from near the wall and mixes with the flow as it travels further downstream. Also, the accumulation of vapor along the nozzle axis downstream of the entrance can be seen clearly. It was argued [15] that the observed asymmetric vapor distribution has been attributed to small machining defects which may have provided additional nucleation sites to the nozzle.
On comparison of the predicted results using for different turbulence models shown in Figure 6, it is clear that the realizable k-epsilon and standard k-epsilon models outperformed the other two models, and that the realizable k-epsilon model predictions give closer agreement with the experimental data than the standard k-epsilon model, especially toward the exit of the nozzle. The results show the realizable k-epsilon model overpredicts the vapor fraction near the entrance and becomes similar as x/L increases, while the standard k-epsilon model underpredicts in the beginning of the nozzle and the difference becomes much wider with downstream distance. This can be correlated with the larger recirculation zone predicted using the realizable k-epsilon model at the beginning of the nozzle (x/L = 0.0167) as shown fully in [48] from the mean axial velocity profiles; the larger recirculation zone corresponds to higher flow acceleration into the nozzle, which in turn reduces the pressure more (or more saturated pressure as shown in [48]), and hence more cavitation leading to the larger prediction of vapor volume.
In addition, the qualitative comparison of the vapor fraction predictions with experiment, Figure 7, confirms the conclusion made in Figure 6 that the realizable k-epsilon and standard k-epsilon models outperform the other two models; the graphs for other two models are not plotted here but can be found in [48]. However, the conclusion drawn between the realizable k-epsilon and standard k-epsilon models is not as clear as that of Figure 6. The comparison of the results in Figure 7a-c, indicate a slightly better agreement between the experiment and the standard k epsilon model than the realizable k epsilon model. The CFD results suggest that vapor is generated at the sharp entrance of the nozzle orifice for both models. However, the generated vapor with the standard k-epsilon models is dispersing from the vicinity of the wall toward axis and mixing with the liquid as the flow is converted downstream, but the cavitation structure disappears at the orifice exit, Figure 7b. Although similar dispersion of the vapor can be seen with realizable k-epsilon model, the vapor never reaches the nozzle axis, and that the cavitation structure extends beyond the exit of the nozzle, Figure 7c. These features are further confirmed when the cavitation contours for both models are plotted outside the orifice by [48], where the results clearly show the extension of cavitation into the circular expansion region of the orifice for the realizable k-epsilon model, while it terminates at the exit of the orifice for the standard k-epsilon model. It is also clear from Figure 7 that the intensity of the cavitation is higher with the realizable k-epsilon model in line with the results obtained in Figure 6.
The observed difference in the behavior of standard k-epsilon and realizable k-epsilon models can be attributed to the turbulent variable C µ used in the realizable k-epsilon model while it is treated as a constant in the standard k-epsilon model. Further analysis made by Kumar [48] showed clearly that the C µ had a direct influence on the turbulent viscosity µ t and the whole flow field so that a lower turbulent viscosity is predicted using the realizable k-epsilon model than the standard k-epsilon. Therefore, the predicted larger recirculation region, mentioned above, with the realizable k-epsilon at the orifice entrance can be attributed to the lower turbulent viscosity that resulted in more saturated pressure region, and thus more intense cavitation. On the other hand, the higher turbulent viscosity predicted with the standard k-epsilon model has resulted in more diffusion between phases, see Figure 7a,b. Therefore, from these analyses, it can be stated that the different turbulent viscosity formulation in the realizable k-epsilon model has resulted in a more accurate prediction of the recirculation region and flow reattachment and led to a more accurate prediction of vapor volume fraction in the present case.
In general, the full analysis given in [48] showed that the volume fraction has been underpredicted using the Launder and Sharma low-Re k-epsilon turbulence model because it has the same transport equations as the standard k-epsilon model, as well as the same turbulent viscosity formulation. However, this model has additional source terms in the transport equations for k-epsilon as well as damping functions that are active only close to the wall allowing the model to predict k-epsilon in the viscous sublayer, hence, the model does not use wall functions. Nevertheless, the results indicate that the model's low-Re approach has failed to completely capture the flow regime for the present case. The model has predicted a smaller recirculation region at the axial position x/L = 0.0167 near the entrance of the nozzle, which resulted in the smaller low-pressure region and therefore lesser cavitation; the analysis also showed that as the fluid moves downstream the cavitation is largely seen to stay in the vicinity of the wall and collapses well before the nozzle exit; around two-third of the orifice length.
The SST k-omega model uses features of the standard k-epsilon in the far-field region and the standard k-omega model [30] in the near-wall region; the k-omega terms are integrated throughout the viscous sublayer hence the model does not require wall functions. In addition, the model also has a different formulation for turbulent viscosity than the k-epsilon family of models. The model, in general, well overpredicts the generation of vapor, Figure 6, but the vapor tends to remain in the vicinity of the wall as fluid moves downstream and enters the circular expansion region of the nozzle. On comparing the mean axial velocity profiles, at x/L = 0.0167, it formed the largest recirculation region (compare to other three models) leading to the largest vapor prediction, Figure 6. However, the cavitation with this model did not diffuse and mixes with liquid like other three cases and remained close to the wall, which can be attributed to different turbulent viscosity formulation in the SST k-omega model that give very low turbulent viscosity leading to lower diffusion. Overall, the performed analysis concludes that the realizable k-epsilon model gave the best overall results for the present case, therefore it is chosen in the present study in later simulations. Readers might refer to [48] to gain more information about the above simulation and analysis on the selection of turbulence model for this study.

Cavitation Model Assessment
As mentioned before, two cavitation models have been assessed, the SS model and the ZGM model. The models are analyzed quantitatively and qualitatively by comparing the predicted results with the experimental data. The quantitative comparison of predicted results with experimental measurements of [15] is presented in Figure 8 and shows the total volume fraction of vapor along the nozzle axis. The comparison indicates that the SS model gives a good agreement with experimental data up to x/L = 0.6 with a difference of up to 7%, while the ZGB model predicts a higher vapor volume within the same range with differences up to 20%. Both models underpredict the volume faction at the exit, in particular, the SS models. This will be explained later when considering the assessment of evaporation and condensation of the two models. To obtain a better understanding of vapor fraction distribution within the nozzle, the qualitative spatial comparison of the predicted vapor fraction contours for both models with X-ray radiography measurements is made and are shown in Figure 9a-c, which suggests that both models have failed to capture the flow features adequately in the core of the nozzle along the axis, in particular at the downstream near the exit where the experimental results show considerable void while they are absent in both model predictions. This can be better realized by the results of both models presented in Figure 10a,b, where they show that the cavitation initiates at the sharp entrance of the orifice after which it converts downstream with the liquid. Both models predict that the significant proportion of Figure 8. Comparison of predicted total vapor volume fraction using two cavitation models with experimentally measured total vapor volume fraction across the nozzle at different axial positions along the nozzle. Error bars in the experimental data show uncertainty of around 2% [15].
To obtain a better understanding of vapor fraction distribution within the nozzle, the qualitative spatial comparison of the predicted vapor fraction contours for both models with X-ray radiography measurements is made and are shown in Figure 9a-c, which suggests that both models have failed to capture the flow features adequately in the core of the nozzle along the axis, in particular at the downstream near the exit where the experimental results show considerable void while they are absent in both model predictions. This can be better realized by the results of both models presented in Figure 10a,b, where they show that the cavitation initiates at the sharp entrance of the orifice after which it converts downstream with the liquid. Both models predict that the significant proportion of vapor remains near the wall with the uniform annular expansion with downstream distance (x/L) so that at the exit, the vapor bubbles cover almost 2/3 of the diameter. After the exit, the vapor bubbles start to reduce considerably so that they are diminished at distances 1.0D and 0.5D downstream of the exit for the SS and ZGB model, respectively; the reason for this is the vapor transforms back to liquid because of higher pressure within the circular expansion region. As aforementioned, both models have failed to predict vapor along the centerline as found in experimental results. To further analyze the predictions, it was decided to plot and compare the evaporation and condensation rate terms of both models. It should be noted that the evaporation and condensation rates in the SS model is non-linear and is proportional to α v (1 − α v ) and has an unusual property; it approaches zero when α → 0 and α → 1 and is maximum in between. Hence, the SS model requires a particular value of vapor fraction at the inlet boundary to initiate cavitation. The value used in the present case for the SS model is α v = 1 × 10 −5 as recommended by [26], which is not the case for the ZGB model. The vaporization and condensation rates in the ZGB model are linear and start from the initial value.   Hence, the ratio of evaporation and condensation rates of both models for the range of 0.005 ≤ α v ≤ 0.995 is plotted and is shown in Figure 11a,b. The evaporation rate graph in Figure 11a shows that the evaporation rate of the ZGB model is higher than that of the SS model in the range of 0.005 < α v < 0.35 after which it becomes smaller than the SS model. The condensation rate ratio graph in Figure 11b also suggests a higher condensation rate for ZGB model in a much small range, i.e., 0.995 > α v > 0.975, following which the condensation rate is greater in the SS model in the remaining of the α v range. From this comparison, it can be stated that the higher vapor volume fraction predicted in the ZGB model was because of the higher initial evaporation rate in the ZGB model over the SS model and lesser overall condensation rate of ZGB model. Hence, the ratio of evaporation and condensation rates of both models for the range of 0.005 ≤ ≤ 0.995 is plotted and is shown in Figure 11a,b. The evaporation rate graph in Figure 11a shows that the evaporation rate of the ZGB model is higher than that of the SS model in the range of 0.005 < < 0.35 after which it becomes smaller than the SS model. The condensation rate ratio graph in Figure 8b also suggests a higher condensation rate for ZGB model in a much small range, i.e., 0.995 > > 0.975, following which the condensation rate is greater in the SS model in the remaining of the range. From this comparison, it can be stated that the higher vapor volume fraction predicted in the ZGB model was because of the higher initial evaporation rate in the ZGB model over the SS model and lesser overall condensation rate of ZGB model. This study also assesses the influence of model constants on cavitation results associated with liquid quality. One of the essential factors that affect the liquid quality is the bubble nuclei population. In real liquids, nuclei may comprise small gas bubbles and impurities [47]. Brennen [4] explained that the bubble nuclei of larger than critical value could assist cavitation more substantially while smaller ones would stay in the liquid. Both models contain constants that need to be adjusted to account for the liquid quality. One such constant in the SS model is , which expresses the bubble concentration per unit volume of pure liquid. The ZGB model also contains many constants such as bubble radius This study also assesses the influence of model constants on cavitation results associated with liquid quality. One of the essential factors that affect the liquid quality is the bubble nuclei population. In real liquids, nuclei may comprise small gas bubbles and impurities [47]. Brennen [4] explained that the bubble nuclei of larger than critical value could assist cavitation more substantially while smaller ones would stay in the liquid. Both models contain constants that need to be adjusted to account for the liquid quality. One such constant in the SS model is n B , which expresses the bubble concentration per unit volume of pure liquid. The ZGB model also contains many constants such as bubble radius R B , nucleation site volume fraction α nuc , F vap , and F cond that need tuning for optimum physical results.
In the initial simulations, the SS model constant n B was varied and ZGB model constant α nuc is adjusted. The results presented in Figure 12 indicate that the cavitation vapor volume in these models can be adjusted by changing model constants. The SS model gave the best result when the value of the model constant n B is 10 13 /m 3 . The ZGB model slightly overpredicts the vapor volume when the value of the model constant α nuc is set to 0.0005 however, when α nuc is set to 0.0001, the model to some extent underpredicts the same. The ZGB model contains other empirical tunable constants such as R B , F vap , and F cond which are scalars and can be adjusted to allow users to further optimize predictions. However, the impacts of varying these constants are not considered here to avoid reader divergence since changing these constants leads to many permutations and combinations leading the present investigation to become merely a curve-fitting exercise. Also, in the above simulations, both phases are considered incompressible which is not physically accurate [4]. Therefore, it is necessary to evaluate both models by considering both phases to be compressible as described in the next section.

Influence of Compressibility
The liquid compressibility is modeled using the Tait equation (Equation (21)), and the ideal gas law (Equation (25)) is applied to model the vapor compressibility. The speed of sound in the bubbly phase is modeled using the Wallis equation (Equation (28)). Figure 13 shows the influences of vapor compressibility by comparing the predicted vapor faction from different models with the experimental. It is evident that the inclusion of the vapor compressibility makes significant improvement and predicts a closer vapor distribution to that of the experimental data, Figure 13a,c,e. However, when the vapor was assumed to be an ideal gas, the results obtained from a fully converged solution using the SS model were non-physical and deviated significantly from experimental results particularly in the core region of the flow as shown in Figure 13c. This is mainly believed to be due to the overprediction of condensation rate (as previously shown in Figure 11). Therefore, the SS model was abandoned at this stage, and that the results obtained using only the ZGB model are

Influence of Compressibility
The liquid compressibility is modeled using the Tait equation (Equation (21)), and the ideal gas law (Equation (25)) is applied to model the vapor compressibility. The speed of sound in the bubbly phase is modeled using the Wallis equation (Equation (28)). Figure 13 shows the influences of vapor compressibility by comparing the predicted vapor faction from different models with the experimental. It is evident that the inclusion of the vapor compressibility makes significant improvement and predicts a closer vapor distribution to that of the experimental data, Figure 13a,c,e. However, when the vapor was assumed to be an ideal gas, the results obtained from a fully converged solution using the SS model were non-physical and deviated significantly from experimental results particularly in the core region of the flow as shown in Figure 13c. This is mainly believed to be due to the overprediction of condensation rate (as previously shown in Figure 11). Therefore, the SS model was abandoned at this stage, and that the results obtained using only the ZGB model are considered for further analysis in this study. The quantitative comparison of predicted vapor volume fraction using the ZGB model presented in Figure 14 clearly shows an increase in the vapor volume fraction when vapor is assumed to be the ideal gas. This is because the ideal gas assumption allowed the vapor volume to change with the change in the local pressure. Therefore, the vapor volume expanded and its density decreased as the pressure decreased. The comparison also shows better agreements achieved with the experimental data when using the ideal gas assumption for cavitation vapor. In addition, another simulation with The quantitative comparison of predicted vapor volume fraction using the ZGB model presented in Figure 14 clearly shows an increase in the vapor volume fraction when vapor is assumed to be the ideal gas. This is because the ideal gas assumption allowed the vapor volume to change with the change in the local pressure. Therefore, the vapor volume expanded and its density decreased as the pressure decreased. The comparison also shows better agreements achieved with the experimental data when using the ideal gas assumption for cavitation vapor. In addition, another simulation with the liquid compressibility was modeled using the Tait equation and the results are superimposed in Figure 14, which clearly show that the liquid compressibility assumption made no difference in the cavitation results, implying that liquid compressibility has an insignificant impact in this case as expected for an upstream-downstream pressure ratio of 10 bars [4]. In this case study, during the assessment of the SS model, the value of n B was increased by a factor of ten from 10 10 /m 3 to 10 13 /m 3 . The results of n B = 10 10 /m 3 are only presented as the results for a subsequent increase in n B did not show any significant physical changes and improvements. While assessing the ZGB model, α nuc was set equal to 0.0001 and was not further changed as it provided fairly good agreement with experimental measurements and with the vapor compressibility model (ideal gas law). The other model constants were not changed.
Fluids 2020, 5, x FOR PEER REVIEW 18 of 24 Figure 14, which clearly show that the liquid compressibility assumption made no difference in the cavitation results, implying that liquid compressibility has an insignificant impact in this case as expected for an upstream-downstream pressure ratio of 10 bars [4]. In this case study, during the assessment of the SS model, the value of was increased by a factor of ten from 10 10 / 3 to 10 13 / 3 . The results of = 10 10 / 3 are only presented as the results for a subsequent increase in did not show any significant physical changes and improvements. While assessing the ZGB model, was set equal to 0.0001 and was not further changed as it provided fairly good agreement with experimental measurements and with the vapor compressibility model (ideal gas law). The other model constants were not changed. As was shown in Figure 13, the qualitative comparison of predicted vapor fraction with experimental measurements, presented in Figure 13a,e, clearly indicates a much better agreement when vapor is assumed to be an ideal gas than incompressible, Figure 13d, and therefore has been adopted to simulate the flow. Prediction results using the ideal gas model for vapor compressibility over the whole flow domain shows that the vapor generates at the sharp entrance of the nozzle orifice. As the fluid travels further downstream, the vapor separates from the wall and moves toward the nozzle center as it mixes with the liquid flow and reaches the nozzle axis (centerline) at around 2.5 diameters downstream of the nozzle entrance as it can be seen in Figure 15a. When the vapor is assumed to be incompressible, the predicted vapor contours are similar to the previous results presented in Figure 9c with no presence of cavitation around the center of the nozzle. As was shown in Figure 13, the qualitative comparison of predicted vapor fraction with experimental measurements, presented in Figure 13a,e, clearly indicates a much better agreement when vapor is assumed to be an ideal gas than incompressible, Figure 13d, and therefore has been adopted to simulate the flow. Prediction results using the ideal gas model for vapor compressibility over the whole flow domain shows that the vapor generates at the sharp entrance of the nozzle orifice. As the fluid travels further downstream, the vapor separates from the wall and moves toward the nozzle center as it mixes with the liquid flow and reaches the nozzle axis (centerline) at around 2.5 diameters downstream of the nozzle entrance as it can be seen in Figure 15a. When the vapor is assumed to be incompressible, the predicted vapor contours are similar to the previous results presented in Figure 9c with no presence of cavitation around the center of the nozzle. Moreover, on further examining the flow field of numerical simulation, in which vapor assumed as the ideal gas and liquid density modeled using the Tait equation, the vapor volume fraction contours in Figure 15a shows that after the nozzle exit the vapor fraction diminishes radially with downstream distance so that all vapor transforms into a liquid at a distance 0.75 away from the exit on the centreline. At the same position ( = 0.75 ) the pressure contours in Figure 15b show a sudden surge in pressure (up to 3.7 bars), which can be seen more clearly on the inserted zoomed picture of the radial pressure distribution. This sudden sharp pressure rise downstream of the nozzle exit can also be seen clearly in Figure 16a at the same location where the pressure profile on-axis shows a sudden rise with a correspondence sharp drop in velocity at the same location as shown in Figure 16b. These sharp rise and fall in pressure and velocity, respectively, in bubbly flows are usually associated with the presence of a normal shock wave [18]. Moreover, on further examining the flow field of numerical simulation, in which vapor assumed as the ideal gas and liquid density modeled using the Tait equation, the vapor volume fraction contours in Figure 15a shows that after the nozzle exit the vapor fraction diminishes radially with downstream distance so that all vapor transforms into a liquid at a distance 0.75D away from the exit on the centreline. At the same position (x = 0.75D) the pressure contours in Figure 15b show a sudden surge in pressure (up to 3.7 bars), which can be seen more clearly on the inserted zoomed picture of the radial pressure distribution. This sudden sharp pressure rise downstream of the nozzle exit can also be seen clearly in Figure 16a at the same location where the pressure profile on-axis shows a sudden rise with a correspondence sharp drop in velocity at the same location as shown in Figure 16b. These sharp rise and fall in pressure and velocity, respectively, in bubbly flows are usually associated with the presence of a normal shock wave [18]. To further investigate the formation of the shockwave, local sonic speeds in the liquid phase, , the vapor phase, , and the mixture phase (cavitation region), and the flow velocity magnitude on the nozzle axis are plotted in Figure 17, which show that local sonic speed in the mixture phase is substantially reduced and becomes significantly smaller than the flow velocity within the nozzle orifice region (−1.3 < < 2.5 ) and up to a distance 0.75 (or = 3.25 ) from the nozzle exit. Figure 17 also indicates that local sonic speed in the mixture phase becomes considerably lower than local sonic speeds of vapor and liquid phases, implying that the local flow in the mixture phase becomes supersonic in the cavitation region. Downstream of the cavitation region at = 3.25 (or 1.63 mm), the fluid transforms back to the liquid state because of the higher local pressure ( Figure  15b and Figure 16a). Thus, the local flow becomes subsonic from supersonic and hence, as a consequence, a normal shockwave is formed.  To further investigate the formation of the shockwave, local sonic speeds in the liquid phase, a l , the vapor phase, a v , and the mixture phase (cavitation region), a m and the flow velocity magnitude on the nozzle axis are plotted in Figure 17, which show that local sonic speed in the mixture phase is substantially reduced and becomes significantly smaller than the flow velocity within the nozzle orifice region (−1.3D < x < 2.5D) and up to a distance 0.75D (or x = 3.25D) from the nozzle exit. Figure 17 also indicates that local sonic speed in the mixture phase becomes considerably lower than local sonic speeds of vapor and liquid phases, implying that the local flow in the mixture phase becomes supersonic in the cavitation region. Downstream of the cavitation region at x = 3.25D (or 1.63 mm), the fluid transforms back to the liquid state because of the higher local pressure (Figures 15b and  16a). Thus, the local flow becomes subsonic from supersonic and hence, as a consequence, a normal shockwave is formed. To further investigate the formation of the shockwave, local sonic speeds in the liquid phase, , the vapor phase, , and the mixture phase (cavitation region), and the flow velocity magnitude on the nozzle axis are plotted in Figure 17, which show that local sonic speed in the mixture phase is substantially reduced and becomes significantly smaller than the flow velocity within the nozzle orifice region (−1.3 < < 2.5 ) and up to a distance 0.75 (or = 3.25 ) from the nozzle exit. Figure 17 also indicates that local sonic speed in the mixture phase becomes considerably lower than local sonic speeds of vapor and liquid phases, implying that the local flow in the mixture phase becomes supersonic in the cavitation region. Downstream of the cavitation region at = 3.25 (or 1.63 mm), the fluid transforms back to the liquid state because of the higher local pressure ( Figure  15b and Figure 16a). Thus, the local flow becomes subsonic from supersonic and hence, as a consequence, a normal shockwave is formed.

Conclusions
A cavitating flow in a submerged orifice was used to evaluate the predictive capability of different cavitation models for compressible flows. The current flow geometry in a single nozzle has been chosen because there was useful detailed experimental data, which were used to assess the models. The performance of different turbulence models was also evaluated for the same flow geometry. The simulation results were compared with experimental data for validation of models and selection of an appropriate model to evaluate the cavitating flow accurately. The following are a summary of the main findings and contributions of the current research study:

1.
Four turbulence models were evaluated, the standard k-epsilon, realizable k-epsilon, Launder and Sharma k-epsilon, and SST k-omega. The results showed that turbulent viscosity profoundly affected the results, and that the realizable k-epsilon model was found to provide the most accurate results which is believed to be due to the formulation of turbulent viscosity.

2.
The performance of the SS and ZGB cavitating models were analyzed and compared. The comparison of results showed higher vapor volume fraction when using the ZGB model because of the higher initial evaporation rate in the ZGB model, and lesser overall condensation rate with the ZGB model. The analyses also suggested that using incompressible liquid and vapor assumption, the models failed to capture the flow physics accurately.

3.
Simulations of the models with both liquid and vapor phase to be compressible showed that the SS model failed to predict the flow physics with a weak agreement with experimental data particularly in the nozzle core region, which is again believed to be due to the overestimation of condensation rate. However, the prediction results of the ZGB model were in better agreement with the experimental data and showed the ideal gas assumption allowed the cavitation vapor volume to change by the change in the local pressure. As a result, the cavitation vapor volume has been improved and increased compared to the incompressible assumption, particularly on the core region. In addition, a normal shockwave was predicted downstream of the cavitation region where the local sonic speed becomes subsonic from supersonic.
Overall, the chosen ZGB cavitation model with inclusion of the compressibility provided very good agreements with the experimental data across the nozzle particularly in the core of the nozzle along the axis and outside the nozzle exit. The outcome of this study will provide appropriate tools for engineers to further improve the calculation of cavitating flows and the designs of the related machines.  vapor volume fraction α nuc nucleation site per volume fraction γ specific heat ratio µ molecular viscosity of the mixture µ l molecular viscosity of the liquid µ v molecular viscosity of the vapor ρ 0 reference density of liquid ρ density of mixture ρ l density of the liquid ρ v density of the vapor = τ mixture shear stress tensor