Assessment of RANS turbulence models and Zwart cavitation model empirical coefficients for the simulation of unsteady cloud cavitation

The numerical simulation of unsteady cavitation flows is sensitive to the selected models and associated parameters. Consequently, three Reynolds Average Navier-Stokes (RANS) turbulence models and the Zwart cavitationmodelwere selected to assess their performance for the simulation of cloud cavitation on 2D hydrofoils. The experimental cavitation tests from a NACA65012 hydrofoil at different hydrodynamic conditions were used as a reference to tune the modeling parameters and the experimental tests from a NACA0015 were finally used to validate them. The effects of near wall grid refinement, time step, iterations andmesh elements were also investigated. The results indicate that the Shear Stress Transport (SST) model is sensitive to near wall grid resolution which should be fine enough. Moreover, the cavitation morphology and dynamic behavior are sensitive to the selection of the Zwart empirical vaporization, Fv , and condensation, Fc , coefficients. Therefore, a multiple linear regression approach with the single objective of predicting the shedding frequency was carried out that permitted to find the range of coefficient values giving the most accurate results. In addition, it was observed that they provided a better prediction of the vapor volume fraction and of the instantaneous pressure pulse generated by the main cloud cavity collapse. ARTICLE HISTORY Received 29 May 2019 Accepted 12 November 2019


Introduction
Cavitation is a topic of concern in the design and operation of a wide variety of hydraulic machinery and systems due to its negative effects like erosion, noise, vibrations and performance drop. Many research has been focused on the numerical simulation of the steady properties of the cavitation forms appearing for example in water jet nozzles (Chen, Hu, & Zhang, 2019), hydraulic turbines (Ayli, Celebioglu, & Aradag, 2016) and valves (Yuan, Song, & Liu, 2019). Broadly speaking, Computational Fluid Dynamic (CFD) models are becoming an indispensable tool to generate precise flow prediction and optimum design in many practical problems as demonstrated by recent works in very diverse disciplines such as architectural fields (Mou, He, Zhao & Chau, 2017), diesel engines (Akbarian, et al. 2018), heat exchangers (Ramezanizadeh, Nazari, Ahmadi & Chau, 2019), nanofluids (Ghalandari, Koohshahi, Mohamadian, Shamshirband & Chau, 2019), and in combination with machine learning methods (Mosavi, Shamshirband, Salwana, Chau & Tah, 2019).
More specifically, cavitation flows are turbulent and unsteady which make them complex fluid mechanics problems. For example, one of the most aggressive forms CONTACT Xavier Escaler xavier.escaler@upc.edu of cavitation is the cloud cavitation that forms on hydrofoils. This powerful and unstable type of cavitation generates strong vibrations on the hydrofoil that are prone to excite the structure and to erode the solid surface. Therefore, the unstable behavior of cavitation must be predicted during the design stage to guarantee a safe operation of many hydraulic machinery such as turbines and pumps. For that, accurate numerical simulations with CFD are needed which depend both on the turbulence and the cavitation models being used. Many works have addressed the significant influence of the turbulence model on the simulation of cavitation flows. For example, Bensow (2011) simulated unsteady cavitation on the Delft Twist11 foil with Large Eddy Simulation (LES), Detached Eddy Simulation (DES), and Reynolds-Averaged Navier-Stokes (RANS) turbulence models. He showed that the RANS models failed to capture the unsteady behavior unless the Reboud's correction (Reboud, Stutz, & Coutier, 1998) was used. This correction reduces the turbulent viscosity in the mixture of water and vapor by taking into account the compressibility effects. He also demonstrated that LES and DES can predict details of the cavitation dynamic behavior like the shedding frequency. Kim (2009) obtained a similar conclusion by simulating partial cavitation on a NACA0015 hydrofoil. Similarly, Pendar and  and Roohi, Pendar and Rahimi (2016) compared the RANS Shear Stress Transport (SST) and the LES models for cavitating and supercavitating flows, and they stated that LES provided the most accurate solutions. Kinzel, Lindau, Peltier, Kunz and Sankaran (2007) applied the DES and the RANS standard k-ε models to simulate the flow around a ventilated body and an ogive. They pointed out that DES can capture a much broader spectrum of the turbulent scales and the cavity dynamics, as well as predict better a range of cavitating flows. Sedlar, Ji, Kratky, Rebok and Huzlik (2016) employed LES, DES and Scale-Adaptive Simulation (SAS) to simulate the unsteady cavitating flow around a NACA2412 hydrofoil. They found that SAS and DES can predict better the dominant frequency of the cavity oscillation, which is overestimated by LES. Nevertheless, LES appeared to provide the best description of the vortex structures.
In spite of the superiority that the scale-resolving approaches (LES and DES) have demonstrated, it has been stated that these approaches become impractical in a fast-paced industrial context due to the significant requirements in computational capacity. Therefore, RANS models still remain the most widely used approach in industrial CFD for hydraulic machinery and systems. However, conclusive results regarding the influence of different RANS models on the cavitation simulation is quiet limited. For instance, Goncalvès (2011) used four different RANS models to simulate the cavitating flow in a Venturi test section. The numerical results showed that the SST model provided the best agreement with the experimental data.
Regarding the influence of the cavitation models, there are various computational approaches to simulate cavitation flows with different levels of complexity. The two-phase flow is often treated as a homogenous phase mixture of vapor and water consisting of a single fluid with varying density. Two strategies can be used to compute the density field which are based on either an equation of state (EOS) or on a volume fraction transport equation.
For example, Goncalves (2011) used a barotropic equation of state (EOS) that couples density directly with pressure to close the system. However, this model cannot take into account mass exchange and thermal transfer in the cavitation. Furthermore, it cannot capture all the vorticity production that is a fundamental mechanism at the cavity closure region (Gopalan and Katz, 2000;Senocak and Shyy 2002). Mani, Cervone and Hickey (2017) also demonstrated that this strategy is very sensitive to the turbulent closure model. Similarly, several authors have found that the morphology of the cavity is highly affected by the value of the minimum speed of sound (Pascarella, Salvatore & Ciucci, 2003;Gonçalves & Patella, 2009;Hejranfar, Ezzatneshan & Fattah-Hesari, 2015).
The approach based on the transport equation models (TEM) for the volume fraction including a cavitation source term permits to calculate the mass transfer between the vapor and water phases (Utturkar, Wu, Wang & Shyy, 2005). Several TEMs have been proposed which have been set as the default option in popular CFD software. For instance, the Zwart model (Zwart, Gerber, & Belamri, 2004) and the Kunz model (Kunz et al., 2000) are the native cavitation models in CFX R and OPENFOAM R , respectively. Meanwhile, Fluent R has adopted the Singhal cavitation model (Singhal, Athavale, Li & Jiang, 2002). However, the main drawback of these TEMs is that they are based on different source terms along with empirical constants. Besides, some ill-posed formulations have been detected in the works of several authors. For example, Pendar and Roohi (2016), Roohi, Pendar and Rahimi (2016), Senocak and Shyy (2002) and Utturkar, Wu, Wang and Shyy (2005) define the Kunz model mass transfer rate with inconsistent dimensions. Morevoer, Senocak and Shyy (2002) and Utturkar, Wu, Wang and Shyy (2005) also propose a different expression for the condensation term of the Singhal model using the product of water and vapor densities instead of the square value of the water density. Therefore, attention should be paid to carefully check the mathematical definitions of the cavitation models implemented by the researchers into their CFD codes in comparison with the original models being used as reference. The TEMs' empirical constants are necessary to tune the evaporation and condensation rates between the two phases which are not symmetric. Actually, the values for such empirical coefficients are in some way arbitrary and generally based on the studies carried out by the model's authors. Consequently, the use of the assumed default values could result in uncertain results for some cases. Therefore, it is crucial to assess the validity and the influence of these coefficients on the numerical results. In this sense, Vaidyanathan, Senocak, Wu and Shyy (2003) optimized the coefficients of the Kunz model based on the response surface method, and identified the best combination for attached cavities around a hemispherical projectile (Rouse & McNown, 1948) and the modified NACA66(MOD) hydrofoil (Shen & Dimotakis, 1989) at different cavitation numbers. Morgut, Nobile and Biluš (2011) and Bilus, Morgut and Nobile (2013) tuned the empirical coefficients of three different cavitation models (Zwart, Kunz and Singhal) with an optimization strategy based on sheet cavitation experiments carried out around  Dupont (1993) and around a NACA66(MOD) hydrofoil by Shen and Dimotakis (1989). The results demonstrated that the three cavitation models could provide similar levels of accuracy if optimized empirical coefficients were used. However, for the case of unsteady partial cavitation (sheet or cloud type), they also noted that further work was required to find out the adequate coefficients. Tseng and Wang (2014) modified the coefficients of the Zwart model into a dimensionless form, and determined a coefficient range also based on the experiments of a hemispherical projectile (Rouse & McNown, 1948) and of the NACA66(MOD) hydrofoil (Shen & Dimotakis, 1989), as well as the Clark Y hydrofoil (Wang, Senocak, Shyy, Ikohagi & Cao, 2001). They stated that their proposed range of values improved the generality and reduced the sensitivity of the numerical results to the cavitation model. As an example, Table 1 summarizes some research works with their recommended values for different cavitation models. Note here that the coefficients recommended by Tseng and Wang (2014) were obtained based on particular dimensionless coefficients. The recommended coefficients in Table 1 show that completely different values were proposed for the Kunz and Zwart models by different authors, although similar cavitation patterns were investigated by all of them. One possible reason for such discrepancies are that the turbulence models, which have a significant effect on the numerical results, are not aligned between them.
Based on the current state of art, the present work has been devoted firstly to investigate systematically and with detail the performance of different RANS models by comparing the simulated unsteady cavitation flows with the experimental results available under different operation conditions. Moreover, the sensitivity of each turbulence model to various numerical parameters has also been evaluated. Next, the influence of the empirical coefficients of the Zwart cavitation model on the simulation results of cloud cavitation dimensions, morphology, dynamic behavior and collapse process have be analyzed in detail to determine a range of optimal values to be used. It has to be noted that, unlike the previous works summarized in Table 1 in which most of the investigations have been based on steady state simulations, in our study we have focused on the shedding process of the cloud cavitation around hydrofoils.

Experimental results
The experimental investigations taken as a reference for validation purposes of our numerical results were carried out by Escaler, Farhat, Egusquiza and Avellan (2007) and Couty (2002) to determine the dynamic behavior and the intensity of erosive partial cavitation on hydrofoils. In particular, unstable cloud cavitation was tested in the High Speed Cavitation Tunnel at EPFL for various free stream velocities, U inf , on a NACA65012 hydrofoil. The tunnel test section was rectangular with dimensions 150 × 150 × 750 mm 3 . The hydrofoil had a chord length, c, of 100 mm and it was fixed with an incidence angle, α, of 6°.
Unsteady cloud cavitation on the hydrofoil suction side was generated by adjusting the inlet pressure, P in . The values of the corresponding cavitation numbers, σ , defined by Eq. 1 are indicated in Table 2 for 8 testing conditions that comprised free stream velocities, U inf , ranging from 15 to 30 m/s and two maximum cavity lengths, l, of 20% and 40% of the chord. The shedding frequency, f, was obtained by amplitude demodulation of the

Numerical model
Based on the homogeneous mixture flow assumption, the two phases of a turbulent cavitation flow, vapor and water, are assumed to have velocity and pressure equilibrium between them. Thus, the governing equations for the mixture quantities are the mass conservation equation (Eq. 3) and the momentum conservation equation (Eq. 4): where u and p are the mixture flow velocity and pressure, t is the time and μ t is the turbulent eddy viscosity. The mixture flow dynamic viscosity, μ m , and density, ρ m , are defined as by Eqs. 5 and 6: where ρ v and μ v are the vapor density and dynamic viscosity, respectively, ρ l and μ l are the water density and dynamic viscosity, respectively, and α v is the vapor volume fraction that is defined as the ratio of the vapor volume to the cell volume and obtained from an additional vapor mass conservation equation (Eq. 7): where R e and R c are the evaporation and condensation source terms, respectively, that account for the mass transfer rate between the water and vapor phases. The mass transfer rate calculation can be developed based on the dynamics of a spherical bubble in an infinite body of incompressible fluid governed by the Rayleigh-Plesset equation (Brennen, 1995) which is defined as Eq. 8: where R is the bubble radius, R 0 is the initial bubble radius, S is the liquid surface tension, P G0 is the pressure of non-condensable gas, P v is the saturated pressure, T ∞ and P ∞ are the temperature and the pressure at infinity, respectively, and T B is the temperature within the bubble.
On the left hand side of Eq. 8, the first term is the driving force determined by the conditions far from the bubble, the second term refers to the thermal effects which will play an important role when the temperature difference is large enough, and the third term refers to the noncondensable gas effect. On the right hand side of Eq. 8, the last two terms consider the influence of viscosity and surface tension, respectively. Since the bubble's growth and collapse are generally considered to be isothermal processes, thus Eq. 8 can be simplified to: Since in most cases the inertial forces are dominant, viscosity and surface tension do not play a significant role and the effects of non-condensable gas are neglected, Eq. 9 is further simplified to: This simplified form of the Rayleight-Plesset equation has been used to develop the cavitation models by Singhal (Singhal et al., 2002), Sauer (Schnerr and Sauer, 2001) and Zwart (Zwart et al., 2004) based on the relation between the bubble diameter and the vapor volume or mass fractions. Nevertheless, its use may affect the model accuracy if, for example, the second temporal derivative of bubble radius is ignored because then the initial bubble growth rate and the bubble collapse rate will be affected. In a similar way, the thermodynamic effect should also be taken into account because it will modify the pressure and temperature distribution between the bubble and the outer liquid and thus the mass transfer rate. Therefore, in order to predict the cavitation phenomenon more accurately, some empirical constants and other parameters such as the vapor bubble radius in the Zwart model, the noncondensable gas fraction in the Singhal model and the number of nuclei per liquid volume in the Sauer model need to be determined and tuned for each flow condition. More specifically, these cavitation models must also be significantly improved for some applications like cavitation in cryogenic liquids, in strongly viscous liquids and in flows with high gas content. In particular, for the Zwart model, the vaporization and condensation mass transfer rates are expressed as: where R B is the vapor bubble radius, α nuc is the nucleation site volume fraction, P v is the saturated water vapor pressure, p is the local fluid pressure, and F v and F c are the empirical coefficients for vaporization and condensation, respectively. The default model constants in ANSYS CFX R v16.2 software are R B = 10 −6 m, α nuc = 0.0005, The RANS models used to calculate μ t are able to reproduce the unsteady cavitation behavior if the compressibility effect of the mixture is taken into account using the method proposed by Reboud et al. (1998). This correction has been proved to reproduce the cloud cavitation shedding process because it is able to reduce the turbulent eddy viscosity in the mixture. Thus, it has been implemented in our simulations by introducing the function f(ρ) defined with Eq. 13: where the exponent n is recommended to have a value higher than 1. In our study, we have selected n = 10 after performing a sensitivity analysis. A brief mathematical description of the RANS twoequation models used in our simulations is given as follows, but more detailed formulation can be found in the CFX help manual (ANSYS, 2015). The k-ε model calculates μ t from the turbulent kinetic energy, k, and its dissipation rate, ε, with Eq. 14: where c μ = 0.09. The RNG model is based on renormalization group analysis of the Navier-Stokes equations. The transport equations for turbulence generation and dissipation are the same as those for the k-ε model (Eq. 10), but the model constant c μ = 0.085. Finally, the SST model, which improves the accuracy of prediction of the onset and the amount of flow separation under adverse pressure gradients, assumes that μ t is linked to k and to the turbulent frequency, ω, via Eq. 15: where F 2 is a function that equals 1 for boundary layer flows and 0 for free shear layers, S is an invariant measure of the strain rate and a 1 = 0.31.

Solution strategy
The computational domain and the selected coordinate system have been implemented in ANSYS CFX R version 16.2 based on the tunnel geometry and they are schematically plotted in Figure 1. The whole computational domain extends 2 chords upstream the hydrofoil leading edge and 4.5 chords downstream the trailing edge, and the thickness of the computational domain is 1 mm. The inlet boundary setup was defined with the corresponding normal velocity equal to U inf , a turbulent intensity of 1% and a water and vapor volume fractions equal to 1 and 0, respectively. The average static pressure was specified at the outlet boundary. The corresponding pressure value was obtained from preliminary simulations without the cavitation model activated. For that, the inlet boundary condition was defined as total pressure, with the static pressure being calculated according to the cavitation number, and the outlet boundary was defined as total mass flow rate. After convergence, the average pressure at outlet boundary was calculated. The obtained values for each case were then used as the outlet boundary condition to simulate cavitation. A no-slip wall condition was set for the top, bottom and hydrofoil surfaces. On the two lateral faces of the fluid domain, a symmetry condition was setup to simulate a 2D flow. Moreover, the vapor saturation pressure, the density and the dynamic viscosity of water and vapor were adjusted based on the experiment temperature of 17°C (P v = 2000 Pa, ρ v = 0.01389 kg/m 3 , μ v = 9.6 E −3 kg/m·s, ρ l = 998.7 kg/m 3 , μ l = 0.0011 kg/m·s).
The pressure-velocity direct coupling method was used to solve the governing equations. The high resolution scheme was used for the convection terms. The second order implicit time scheme was used for the transient term. Several successive iterations were set within each physical time step. A very small residual criterion of 10 −8 and a large iterative number were set to march the solution towards convergence in every time step. To

Numerical verification and validation
Some CFD uncertainties might be due to numerical and modeling errors caused by time and space discretization and by incomplete iterative and grid convergences. The y + has been defined with Eq. 16.
where τ ω is the wall shear stress, v is kinematic viscosity, y denotes is the distance between the first and second grid points off the wall. The effects of y + , time step, number of iterations and number of mesh elements have been evaluated only for the cavitation conditions corresponding to case 8 in Table 2. The hydrodynamic conditions of case 8 generate the most aggressive form of cloud cavitation in terms of erosion and vibrations. On the other hand, the modeling errors coming from the assumptions and approximations of the RANS turbulence model and the Zwart cavitation model have been assessed for all the cases in Table 2.

Mesh convergence study
A mesh convergence study was performed according to the Grid Convergence Method (GCI) provided by Fluids Engineering Division of ASME (Ismail et al., 2008) and taking into account that in CFX different turbulence models use different wall treatments. The SST model uses the so-called Automatic Wall Function that switches automatically from a typical wall function approach to a low-Re-number model by blending the wall value for ω between the logarithmic and the near wall formulation as the mesh is refined. Meanwhile, the k-ε and the RNG k-ε models use the so-called Scalable Wall Function in which y + is calculated as y + = max (y + , 11.06). Therefore, these two turbulence models do not resolve the viscous sublayer and they directly use the logarithmic relation to compute the near wall velocity.
Several mesh resolutions were tested with the three turbulence models. For that, the meshing tool ANSYS  ICEM R was used to create two structured meshes, named Grid1 and Grid2, with three different degrees of refinement as indicated in Table 3. All the meshes were built with similar topology and with the same recommended grid refinement factor r = √ 2. For Grid 1 meshes, the height of first layer was small enough so that the boundary layer could be solved with the low-Re-Number model. For Grid 2 meshes, the first layer of elements lied in the log-law region. As an example, Figure 2 shows the topology of G2.2 mesh with an O-grid around the hydrofoil embedded inside an H-grid.
The three turbulence models were applied to each grid to simulate case 8 from Table 2 corresponding to U inf = 30 m/s without considering cavitation. The lift coefficient, C L , defined by Eq. 17, where F L represents the lift force, c is the chord length and s is the span length, was computed to obtain the value of the fine-grid convergence indexes (GCI 21 and GCI 32 ) indicating whether calculations with additional grid refinement should be performed. Table 4 presents these indexes for each turbulence model and each mesh.
As shown in Table 4, the maximum GCI values for Grid 1 are 0.6% and 1.8% for the three turbulence models, indicating that the solution is well within the asymptotic range. For Grid 2, the GCI values are smaller for the k-ε and RNG models, while they are larger for the SST model reaching 9.7%. This is because the mesh resolution in Grid 2 is more consistent with the requirement of a scalable wall function. Based on these results, it was concluded that the medium grids G1.2 and G2.2 provided the necessary resolutions within the asymptotic range, and, considering the cost of simulation time, they were chosen for all the final simulations.

Reboud's correction for the turbulence models
To emphasize the importance of the Reboud's correction and the influence of its exponent value n, the SST turbulence model was selected to conduct the unsteady simulation for case 8 conditions with n = 1 (not recommended) and n = 10 (recommended). The value of n determines the rate of change of the modified effective density, f(ρ), with the water volume fraction, α l , as shown in Figure 3a.
For n = 1, the unstable behavior of the attached cavity is not correctly simulated and a quasi-steady behavior is predicted. However, a typical cloud cavitation behavior with a good shedding frequency and maximum cavity length is obtained with n = 10 as expected. The reason for this is the induced reduction of the eddy viscosity that is clearly overpredicted by the original turbulence model. In Figure 3b, the simulated mean vapor volume fraction and turbulent eddy viscosity obtained with n = 1 and n = 10 are plotted for comparison. It can be seen that, taking n = 1, very high eddy viscosity values are obtained in the cavity closure region which prevents the formation of the re-entrant jet. Whereas with n = 10, the eddy viscosity is almost zero on the hydrofoil extrados except the area close to the trailing edge thus allowing the formation of the re-entrant jet and reproducing the cavitation shedding process. The evaluation of the turbulence models for cavitation simulation is necessary in the particular field of hydraulic machinery because the Reynolds numbers are very high and the losses induced by viscosity are not negligible. Nevertheless, it must be noted that in some cases the unsteady cloud caviation behavior is mainly controlled by the inertial forces instead of the viscous forces. Consequently, inviscid solvers could also be used as demonstrated by the works of Budich, Neuner, Schmidt and Adams (2015) and Schenke and van Terwisga (2019) who modeled with good success the cavitation unsteady behavior over a sharp wedge and around a NACA0015 hydrofoil, respectively. They were able to predict the frequency of the shedding process, the re-entrant jet and the horse-shoe vortices. Thus, the application of inviscid solvers must be considered as a feasible alternative methodology in particular for inertia driven cavitating flows.

Sensitivity of the numerical parameters
In this section, a sensitivity analysis of the main parameters to simulate unsteady cavitation is presented based on case 8. In particular, attention has been given to y + , time step, number of iterations and number of mesh elements. For each computation, both the lift coefficient, C L , and the total vapor volume within the computational domain, V cav , during several cavitation cycles have been considered to assess the results. V cav has been defined with Eq. 18: where n denotes the number of mesh elements, α i is the vapor volume fraction within the element and V i is the element volume. The simulated transient values of these two physical quantities have been compared between different model setups in the time and the frequency domains. Mesh G1.2 with an average y + = 2 and mesh G2.2 with an average y + = 20 have been used to evaluate the influence of the near wall grid refinement on the numerical results for each turbulence model, as it can be seen on the shape of the time signals and their spectra plotted in Figure 4.
In Figure 4a, the results with the SST model show that V cav is significantly higher with the finer mesh. The evolution of C L over time also changes significantly depending on y + . The periodic frequency of the shedding process is of about 233 Hz for the finer mesh and of about 340 Hz for the coarser one, the former estimate being closer to the experimental one. Then, the influence of y + on the results predicted by the k-ε is not as significant as for those with the SST model as demonstrated by the similar results and frequency peaks with only a small difference of 30 Hz shown in Figure 4b. And finally for the RNG model, the shedding process is not well captured with the finer mesh as no clear frequency peak can be distinguished in the spectra in Figure 4c. Meanwhile, a dominant fluctuation is obtained around 300 Hz with the coarser mesh. In conclusion, the SST model appears to work well for smaller y + values. Meanwhile, the k-ε and the RNG models can only capture the cavitation periodic dynamic behavior for larger values of y + . Hence, mesh G1.2 was finally used with the SST model and mesh G2.2 was used with k-ε and RNG models.
Next, the numerical sensitivity to the rest of parameters -time step duration, number of iterations and number of mesh elements-was evaluated based on the calculated dominant frequency peak, f, for both C L and V cav results as indicated in Table 5. The last column represents the percent deviation of the frequency result obtained with the higher value of the sensitivity parameter relative to the result obtained with the lower value of the sensitivity parameter.
Two different time step durations were investigated with an approximate value of 1/100 and 1/200 times the measured shedding period in the experiment, respectively. As it can be seen, the influence on f is weak for all the turbulence models since the maximum error found is only around 3.4%. Based on this result, the time step for the rest of simulations was fixed to 0.00005 s. The influence of the number of iterations in every time step loop appears to be negligible for all the three turbulence models. The change of f is less than 2%. Therefore, a number of 20 iterations was considered sufficient in any case. To test the influence of the number of mesh elements, a new mesh with four times the original number of elements was built but with the same mesh topology and the same average y + . Based on the results presented in Table 5, the obtained f only suffers a slight change with less than 3% deviation. Thus, it was confirmed that the coarser meshes with less than 30,000 elements were suitable to carry out accurate predictions of the cavity dynamic behavior. In summary, and according to the slight effect of the checked numerical parameters, the finally selected values corresponded to those providing accurate numerical results and saving computational time and effort

Sensitivity of the dimension space
To investigate the effect of the number of computational domain dimensions, a 3D model was created according to the experimental configuration with a thickness of 150 mm corresponding to the tunnel test section and hydrofoil span wise size. A 3D mesh was created by extruding the 2D mesh G1.2 with 150 uniformly distributed elements in the cross direction Z, and hence nearly 4.5 million elements were needed to solve the problem. The same boundary conditions and numerical settings were fixed between 2D and 3D simulations, and no-slip wall conditions were used at the two sidewalls of the 3D domain instead of the symmetry conditions used for the 2D domain. The 3D cavitation simulation was only conducted with the SST turbulence model. Regarding the results for case 8, the 3D simulation with the SST model predicts a shedding frequency of 208 Hz meanwhile the 2D simulation predicts 233 Hz and the experiments indicates 226 Hz. Therefore, there is a difference of about 12% between the 3D and the 2D results, and the 3D simulation predicts a lower frequency than the 2D simulation. An explanation for such frequency reduction can be obtained from the comparison between the 2D and 3D numerical results shown in Figure 5. The top graphs show the maximum cavity length of a isosurface with vapor volume fraction α v = 0.2, the middle graphs show the instant of cavity break off and the bottom graphs show the velocity distribution at the cavity break off instant on a surface located 0.5 mm above the hydrofoil extrados. Note that for visualization convenience, the 2D results have been enlarged 150 times in span wise direction to be compared with the 3D results. It can be observed that the maximum cavity length predicted by the 3D simulation is slightly longer than that predicted by the 2D simulation. Moreover, the re-entrant jet structure in the 3D results is not uniform in span wise direction, as observed with the velocity distribution close to the hydrofoil extrados, and such irregularities reduce the average upstream velocity of the re-entrant jet. Therefore, both results explain why a longer time is needed by the re-entrant jet to break up the attached cavity and consequently the shedding frequency is reduced.
In spite of the differences observed between the 3D and the 2D simulation results, a significant similarity is found between them in general terms regarding the shedding process and the average pressure and velocity fields. In conclusion, given the objectives and the scope of the present study it has been decided to base our investigation on 2D simulation results and avoid the extremely high computational cost required for a full study with 3D models.

Sensitivity of the turbulence model
In Table 6, the f predicted by the three different RANS models has been compared with the shedding frequency measured experimentally for all the cavitation conditions described in Table 2, and the corresponding percent errors are also indicated. Note that the value of f was obtained from the time duration of at least ten periodic cycles or from the Fourier analysis of the simulated signals.
The obtained results show that f is always over predicted by all the RANS models with the only exception of the STT model for case 3. In general, the SST model provides more accurate results with maximum errors below 10% for all the cases. Instead, the k-ε and RNG models present larger errors for almost all cases. Further, the RNG model is not able to predict f for case 5. In conclusion, the SST model is the best one to simulate the unsteady periodic behavior of cloud cavitation with the current numerical parameters and a good boundary layer resolution. Therefore, all the final results presented in the following sections were obtained using the SST model and the previously stated model setup parameters.

Simulation of unsteady cavitation with empirical coefficients F v = 300 and F c = 0.03
The results presented in the previous section were calculated with the Zwart default empirical coefficients F v = 50 and F c = 0.01. For comparison, the numerical simulations with the SST model of all the cases listed in Table 2 were computed again using the Morgut et al. (2011) recommended empirical coefficients F v = 300 and F c = 0.03. Tables 7 and 8 show the calculated frequencies and maximum cavity lengths, respectively, simulated by both the default and the tuned coefficients, and the corresponding percent errors. Following the same criteria than for the experimental measurements, the predicted and measured cavity length, l, corresponds to its maximum size in chord wise direction during a shedding period before it detaches the hydrofoil surface.
These results indicate that the default coefficients provide more accurate values of f and l than the ones recommended by Morgut et al. (2011). For the default coefficients, the frequency error is less than 10% for all the cases, and the cavity length error is a little larger just for case 5 reaching around 17%. However, with Morgut's coefficients the model predicts values far from the experimental ones, especially for cases 1-4 with deviations from −24% to −38% in frequency estimates. This is due to an over prediction of l that, as indicated in Table 8, can reach a deviation of about +87.5%. On the contrary, for cases 5-8 the results with Morgut's coefficients are very accurate with errors below 6% in frequency estimation.
To show more details about the differences between the numerical results obtained with default and the Morgut's coefficients, the time history of V cav , C L and the local pressure at location x/c = 0.5 on the suction side of the hydrofoil calculated for each time step (0.00005 s) have been plotted in Figure 6 for case 8. In these plots, the horizontal axis corresponds to the dimensionless time, T*, which has been normalized relative to the characteristic cavitation shedding period T = 1/f and calculated as T* = t/T. The evolution of these values is only shown during two shedding cycles although a total of 10 cycles were simulated to achieve a stable and repetitive solution in all the cases. Then, Figure 7 shows the contour plots of vapor volume fraction at corresponding instants of time during one complete shedding cycle simulated with both sets of coefficients.
From Figure 7, it can be observed that the amount of vapor inside the attached cavity and the shed cloud increases significantly when the value of F v is increased  from 50 to 300. As a result, a large vapor volume fraction gradient is observed at the interface between the cavitation structures and the surrounding water. Moreover, the thickness of the cavitation interface is significantly reduced. In conclusion, the obtained results with F v = 300 appear to be closer to the observed cavitation morphology in the experiments by Escaler et al. (2007) and Couty (2002). These results also present similar vapor volume fractions to those measured with high-speed visualization and time-resolved X-ray densitometry measurements by Ganesh, Mäkiharju and Ceccio (2016) in a periodically shedding cavity forming from a wedge at the end of its growth.
On the other hand, the increase of F c from 0.01 to 0.03 permits to capture the instantaneous collapse of the shed cloud and the resulting pressure pulse on the hydrofoil surface, as it can be observed in Figure 6c. Simultaneously to the pressure peak, the C L suffers a sudden decrease. This pressure pulse could be the cause of the experimentally observed formation and propagation of a bubbly shock within the high void-fraction bubbly mixture in the separated cavity flow by Ganesh et al. (2016). They propose that the periodic cloud shedding could be induced by both a re-entrant jet and a shock-wave traveling upstream. Our simulations predict the re-entrant jet that pinches off the attached sheet of vapor and creates the cloud cavity that is convected downstream by the main flow. Unfortunately, our numerical model is based on an incompressible solver and thus it cannot capture the shock-wave formation.
Moreover, the final cloud collapse takes place behind the closure region of the main cavity using Morgut's coefficients, as it can be seen at instant 1.4T* in Figure 7b, which is in accordance with the region where the material erosion was measured in the experiments (Couty, 2002). It has to be noted that with the default coefficient F c = 0.01, the collapse takes place close to the trailing edge of the hydrofoil which is too far from the expected eroded region, as it can be observed at instant 2.0T* in Figure 7a.
In conclusion, the use of the Morgut's coefficients provided a more realistic simulation of the cavitation morphology and of the collapse process. Nevertheless, for some of the cases they failed to obtain the shedding frequency mainly because the length of the attached cavity was not the correct one. Therefore, it was decided to carry out a parametric study of the empirical coefficients so that the best combination of values predicting the dynamic behavior, the vapor content of the cavities and the pressure pulse on the hydrofoil surface due to the cloud collapse could be found for all the cases indicated in Table 2.

Parametric analysis of Zwart empirical coefficients
For the parametric analysis, the following range of values has been taken into account: • F v from 100 to 500. • F c from 0.02 to 0.10.
For each case indicated in Table 2, a series of transient simulations were carried out using different combinations of values for F v and F c , which comprised a total of 200 computational runs. As a result, a full factorial design space was obtained with 25 design points evenly spaced. Then, a response surface method was applied taking as the response variable the error between the simulated shedding frequency, f, and the measured one in the experiments, and as the independent variables the dimensionless values F * V = F V /50 and F * C = F C /0.01. To better fit the numerical results, a polynomial linear regression including the second order term and the interaction effects between the independent variables was obtained for each case. For example, the regression function obtained for case 8 is given as follows in Eq. 15: To show the accuracy of this function, Figure 8 permits to compare the frequency error of the CFD numerical results (bold) with those obtained with the function (cursive). It can be seen that the corresponding values are generally in good agreement. The maximum residual is found for the design point (300, 0.08) for which the frequency error predicted by CFD is −0.4% meanwhile it is 7.9% when predicted by the regression function. Hence, a total maximum residual of about 8.4% is found for this particular combination of values. Table 9 shows the statistical quality indicators of the regressions obtained for all cases. For all of them, the value of R 2 is close to 1, which indicates an accurate data fitting. The value of significance F is always less than 0.001, indicating that the model and dependent variable are statistically significant, and that the regression equation does have validity within the fitted data. Finally, it is noted that the maximum residual between the CFD and the regression model is less than 10% for all cases.
The obtained response surfaces are shown in Figure 9 for each cavitation case presented in Table 2 with their corresponding contour plots. It can be seen that for cases 1-4, the frequency error varies within a large range from −40% to 30%. On the other hand, the error range is significantly reduced for cases from 5 to 8 within the same design space. This observation suggests that the simulation of the dynamic behavior is more sensitive to the values of the empirical coefficients when the cavity length is shorter. Here it must be recalled that the first 4 cases correspond to a maximum length of 20% of the chord and the last 4 cases correspond to a maximum length of 40% of the chord.
In order to find a range of values that could predict with acceptable accuracy the shedding frequency for each cavitation condition including short and long attached cavity lengths, the frequency errors from cases 1 to 8 were averaged according to Eq. 20 and the average contour plot obtained has been plotted in Figure 10a.
In Figure 10b, the variance of the averaged frequency error for all cases has also been plotted. From theses graphs, a blue region with the smallest error range can be clearly identified at the upper right of the design space, that also presents the smallest variance. This region spans   from F V = 300 to 500 and from F C = 0.08 to 0.1. Therefore, any combination of empirical coefficients within these values should provide good estimates of the cavitation dynamic behavior independently of the cavitation size and hydrodynamic condition.
To validate the range of optimal empirical coefficients, Table 10 shows the frequency errors of the CFD results calculated with couples of coefficients within this region. The maximum frequency error found is of about 13% for case 6 and coefficients (300, 0.08). Moreover, the average value of the absolute frequency errors for each combination of coefficients spans from 5.1% to 7.0% as indicated in the bottom row, which correlates with the error range of the blue region in Figure 10a (4.9% to 7.4%). In conclusion, the feasibility and the accuracy of the results obtained by the response surface are demonstrated.
In addition, it can be observed that the optimal coefficients in Table 10 are significantly higher than the default empirical coefficients (50, 0.01) which underestimate the cavity vapor content and the intensity of the collapse process as discussed before. These results are in accordance with the recent works by Ghahramani, Arabnejad and Bensow (2019) and Schenke, Melissaris and van Terwisga (2019). Using different cavitation models, they found that the speed of the bubble collapse is significantly underestimated with low mass transfer coefficients. Moreover, they provide novel approaches which could be a future line of research for the current work.

Validation of the optimal empirical coefficients on a NACA0015 airfoil
To confirm the validity of the optimal range of coefficients found with the present study, the cloud cavitation behavior visualized and measured around a NACA0015 hydrofoil by Van Rijsbergen, Foeth, Fitzsimmons and Boorsma (2012) was simulated with the default coefficients (50, 0.01) and with two optimal combinations: (300, 0.08) and (500, 0.1). The SST model was used and a mesh with an average y + value less than 2 was created. The boundary conditions were set according to the Table 11. Shedding frequency measured and predicted with default and the optimal coefficients found in the current study for cloud cavitation on a NACA0015 hydrofoil. experiment description, and the same resolution strategy was followed as for the current NACA65012 model. Table 11 shows the measured and the predicted shedding frequency of the cloud cavitation around the NACA0015 airfoil. As it can be seen, the default coefficients overestimate the frequency by around 48.9%. Meanwhile, the shedding frequency is better predicted and the error is reduced to around 9.6% for the tuned couple of values F v = 300 and F c = 0.08.

Conclusion
In the present work, transient simulations of unsteady cloud cavitation on a 2D NACA65012 hydrofoil have been carried out systematically to assess the influence of the setup parameters, RANS turbulence models and Zwart cavitation model empirical coefficients on the numerical results. In summary, it has been concluded that: • The results are more sensitive to near wall mesh resolution than to time step, number of iterations and number of mesh elements. For the SST model, an average y + = 2 must be used, meanwhile for the k-ε and RNG models a coarser grid resolution is sufficient. • The SST model performs better than the k-ε and the RNG models when the Reboud's correction is used. • Increasing F v , larger amounts of vapor inside the cavities are obtained with thinner interfaces between vapor and water phases, which resembles more precisely the experimental observations. • Increasing F c , the instantaneous collapse of the shed cloud is captured and the induced pressure pulse on the hydrofoil surface at the main cavity closure region where erosion takes place is well simulated. • With a parametric study based on a response surface method and multiple linear regression, the optimal range of empirical coefficients to simulate the dynamic behavior of a wide range of cavitation conditions with different attached cavity lengths on 2D hydrofoils has been found. The recommended values are from 300 to 500 for F v , and from 0.08 to 0.1 for F c .
The present parametric study has been limited to a single objective corresponding to predict the cloud cavitation shedding frequency on 2D hydrofoils. Further research should incorporate multiple targets such as velocity and vapor volume fraction distributions, induced pressures and vibrations, that would be considered simultaneously in a multi-objective approach. However, this work will require advanced experimental measurements with more detailed results. In this sense, the compressible approach should be used to try to capture the bubbly shock propagation as a mechanism for sheet-to-cloud transition of partial cavities.
Since the current recommendations are only valid for cavitating 2D hydrofoils in high speed water tunnels, the next step would be to simulate the 3D unsteady cavitation phenomena in hydrofoils and other simple geometries to verify the validity of the present results. And finally, the models should be validated for cavitation in actual hydraulic machines such as turbines and pumps.