Cavitation is an undesirable phenomenon in hydraulic systems, as it causes erosion and noise. The main difficulty in cavitation prediction when using Fluent software is lack of an openly accessible tool for implementation of a freely chosen homogeneous cavitation model. In this paper the main challenge is to make such a tool, user defined function (UDF). The second challenge is to use a qualitative method in the assessment of the results of verification process. Three cavitation models are verified in Fluent 14.5: Singhal *et al.*, Schnerr & Sauer and Zwart *et al.* The verification is based on the benchmark example from the *Cavitation Modeling* tutorial. Three methods of the algorithms verification are used: analysis of the convergence history of volume fraction, comparison of vapour volume fractions and statistical analysis of these data. The original achievements are not only the verified codes but also statistical analysis based on the computer methods of image analysis performed using two correlation coefficients: the first based on the cavitation intensity, and the second based on the changes of the cloud shape. The results of the analyses do not give any reasons to reject the UDFs. The appendix contains the analysed codes (available with the online version of this paper).

## INTRODUCTION

The design of hydraulic systems and their elements require experimental measurements and numerical simulations. Scientists interested in spillways and dams first emphasized the influence of the aeration phenomenon on the considered flow (Castillo *et al.* 2014, 2016; Bayon-Barrachina & Lopez-Jimenez 2015). A similar tendency is visible also in closed hydraulic systems, such as pipe systems. Aeration is seen as a way to avoid the undesirable cavitation phenomenon (Najafi *et al.* 2012). For that reason, in numerical simulations, a multiphase flow with two phases, i.e. water and air, is considered. Cavitation in computational fluid dynamics (CFD) calculations and during the design stage is usually neglected. It should be emphasized that aeration helps to avoid cavitation, but it does not give any guarantee that cavitation will not occur. Additionally, there are some devices requiring cavitation in order to work correctly, e.g. converging-diverging nozzles (Ashrafizadeh & Ghasemmi 2015). Complex design works in hydraulic systems should also consider the cavitation phenomenon.

Cavitation is a phenomenon that consists in appearance, growth and disappearance of bubbles containing vapour of a given liquid. The bubbles grow in areas where the static local fluid pressure drops below the saturated vapour pressure. Next, the bubbles rapidly decrease and implode in areas of the increased pressure. The threshold value for pressure to start the cavitation process could be significantly lower than the saturated vapour pressure for homogeneous liquids or higher for liquids with a large amount of gases (Bagieński 1998).

Cavitation has been a topic of research since 1894 (Reynolds 1894). Thorneycroft & Barnaby (1895) subordinated this phenomenon to the wear of the surface of a screw propeller. The negative influence of cavitation, e.g. erosion and vibration, also started to be known in elements of hydro-turbines and parts of hydraulic systems. Experimental investigations are the first way to learn more about this danger for machines. To have a comprehensive picture of this phenomenon, a complex diagnostic system should be used. The measurement techniques of cavitation can be divided into two categories: contact and contactless. The contactless methods predominate in experimental investigations. To these methods we can classify high-speed imaging, particle image velocimetry (PIV), X-ray attenuation and optoelectronic systems. The oldest and most popular is the high-speed imaging method. A large number of frames per second, i.e. even millions, allows registering of even the smallest development stages of the cavitation area (Kravtsova *et al.* 2014). The PIV technique has been applied in the investigation of cavitating flow since 1995 (Tassin *et al.* 1995). The main restriction for potential users is the huge cost of this system. The X-ray is one more option for capturing images of the cavitation area. Data acquisition rates during X-ray measurements are high (more than 1 kHz) but they do not achieve the frequency values available using a high-speed camera. Stutz & Legoupil (2003) presented the first experimental results of cavitating flow investigations using X-ray. The double optical probe is an example of the contact method. The aim of these measurements is to extract the local void fraction from a liquid-vapour mixture (Stutz & Reboud 1997). The optoelectronic system registering the shape and intensity of cavitation cloud is one of the newest achievements in the field of investigations of cavitating flows. Vapour bubbles disperse the laser light and through it, and the signal measured via digital oscilloscope changes. The method is simple and does not require any large costs (Lipiński & Niedźwiedzka 2016; Niedźwiedzka & Lipiński 2016).

*et al.*2000; Schnerr & Sauer 2001; Singhal

*et al.*2002; Zwart

*et al.*2004), while in the barotropic approach a correlation between mixture density and pressure is used (Delannoy & Kueny 1990; Ventikos & Tzabiras 2000; Liu

*et al.*2004).

The paper focuses on the homogeneous approach. Among the numerical software for calculations of cavitating flows, there are many tools for simulations using homogeneous models. One of these tools is Fluent, part of Ansys. In this software, there are two possibilities to implement the homogeneous models: directly from interface and using user defined functions (UDFs). Only three models are available directly from interface: Schnerr & Sauer (2001), Singhal *et al.* (2002) and Zwart *et al.* (2004). In the literature, the number of homogeneous models is greater. To have more possibilities in numerical simulations in Fluent, application of the UDFs is necessary. There are few works concerning numerical simulations of cavitating flows with implementation of homogeneous models using the UDFs (Susan-Resiga *et al.* 2003; Bernad *et al.* 2004; Zhang *et al.* 2008; Cao *et al.* 2009). Additionally, these works do not refer to any of the existing homogeneous models, but to the Rayleigh equation. There is no information on the verification of the implementation tool, so it is impossible to assess the confidence level of the whole modelling process. The lack of a verified openly accessible tool for implementation of homogeneous cavitation models in Fluent software motivated our work.

The main goal of the presented article is to verify the UDFs defined using the *Define Cavitation Rate* macro for the three above-mentioned homogeneous models*.* The AIAA (American Institute of Aeronautics and Astronautics) Guide (Oberkampf *et al*. 1998) gives two possibilities to verify computational physics codes: exact analytical solution and benchmark solution. In the case of Fluent software, a list of benchmark examples is accessible in Ansys Fluent Tutorial Guide (2012). Assuming that benchmark solutions are characterized by a high accuracy, this way seems to be the best solution. In the article, the benchmark example from the *Modeling Cavitation* tutorial is used for the verification process in Fluent. Oberkampf *et al.* (2004) define two main principles, which should be followed when determining the assessment of the agreement level of the homogeneous cavitation models defined in *Define Cavitation Rate* macro with the same models available directly from the interface. Firstly, verification and validation are ‘processes of determining’. It means that both methods have an ongoing character. The expected finish point of solution is not specified. For that reason, it is impossible to speak about the veracity of the analysed conceptual or computerized models, or about its rejection. Secondly, the level of accuracy should be correctly understood in relation to an engineering perspective. The most popular definition of accuracy describes it as a correctness measure. It is worth emphasizing that from an engineering point of view, the exact solution is not required. It is sufficient to have a meaningful comparison. Whereas the main challenge for the authors is to make a verified tool for implementation of homogenous cavitation models in Fluent software, the second challenge is to use a qualitative method of the assessment of the results of the verification process. Consequently, the original achievement is not only the code itself, but also use of the computer methods of image analysis in the statistical analysis of the results of the verification process. The quantitative analysis is made using two correlation coefficients. A positive evaluation of the prepared tool allows using it for numerical simulations in any specified hydraulic systems and other applications. This information is necessary for future implementation of other homogeneous cavitation models.

## THEORETICAL BACKGROUND

### Single-fluid method based on the homogeneous approach

The crucial moment for the homogeneous approach was 1992. In this year, Kubota *et al.* (1992) presented the first homogeneous model called Local Homogeneous Model (LHM). The instability of this model was inspiration for other scientists to look for new ways in the formulating of the source terms. Merkle *et al.* (1998) showed the next widespread known model. Their mass transfer rates formulation was not the only response for the pioneering equation. Since then, new versions of source terms appear regularly. Several works (Senocack & Shyy 2001; Frikha *et al.* 2008; Goel *et al.* 2008; Žnidarčič *et al.* 2011; Goncalves & Charrière 2014) tried to set the ideas in order, but because of the passage of time, an update of information should be regularly made. The last work with an overview of homogeneous models comes from Niedźwiedzka *et al.* (2016). This work contains not only the final forms of source terms, but also detailed descriptions including advantages and disadvantages as well as examples of application areas.

*ρ*is density (kg/m

^{3}),

*t*is time (s), is velocity (m/s), is spherical stress tensor (Pa), is viscous molecular stress tensor (Pa), is turbulent Reynolds stress tensor (Pa), is intensity of the mass forces source (N/m

^{3}),

*e*is energy (J), is molecular heat flux (kg/s

^{3}), is turbulent heat flux (kg/s

^{3}), is intensity of the energy source (J/(m

^{3}/s)) (Sobieski 2011). The additional transport equation, which expresses the changes in the fraction of the chosen constituent (liquid or vapour), appears most commonly as in the expression shown below: where

*α*is the volume fraction source term for evaporation (1/s),

^{+}*α*is the volume fraction source term for condensation (1/s) and

^{−}*α*is the vapour volume fraction (−). It considers vapour as the constituent and analyses the changes in its volume. In the form of mass fraction, the equation has the following form: where is the mass source term for condensation (kg/(m

_{v}^{3}·s)), is the mass source term for evaporation (kg/(m

^{3}·s)) and

*ρ*is the vapour density (kg/m

_{v}^{3}). For the second possible constituent liquid, the above-mentioned equations are formulated as follows: where

*α*is the liquid volume fraction (−) and

_{l}*ρ*is the liquid density (kg/m

_{l}^{3}).

*p*is saturated vapour pressure (Pa). In the calculations, the densities of liquid and vapour are assumed to be incompressible, and only the mixture density varies.

_{sat}*et al.*2004) use the Rayleigh (1917) equation: where

*R*is the bubble radius (m), to the description of the bubble dynamic. Plesset & Prosperetti (1977) presented a new more comprehensive form of this equation: where

*σ*is surface tension (N/m),

*μ*is liquid dynamic viscosity (Pa·s). The authors of the third model (Singhal

_{l}*et al.*2002), inspired by this formula, analysed the validity of consideration of its selected parts for numerical simulations of cavitation phenomenon, and proposed their own transformations, the most useful from their point of view.

### Mass transfer models

*et al.*(2002) and Zwart

*et al.*(2004). The Schnerr & Sauer (2001) model considers only quantitative values of the physical parameters and so it is widely used in numerical simulations of cavitating flows. Their source terms are formulated as follows: where

*ρ*is mixture density (kg/m

_{m}^{3}).

*et al.*(2002) model, called Full Cavitation Model. The name comes from the large amount of physical parameters that are taken into consideration through the expression of source terms in the following way: where

*C*is condensation constant (−),

_{p}*C*is evaporation constant (−),

_{d}*k*is local turbulent kinetic energy (m

^{2}/s

^{2}),

*f*is vapour mass fraction (−),

_{v}*f*is non-condensible gases mass fraction (−).

_{g}*et al.*(2004) model, the empirical constants appear, as in the Singhal

*et al.*(2002) model. In the evaporation rate, the vapour volume fraction (

*α*) is replaced with the product of the nucleation site of volume fraction (

_{v}*α*) and the remaining fluid volume fraction (1 −

_{nuc}*α*)

_{v}*.*The source terms are formulated as follows:

## TEST CASE

*Modeling Cavitation*(Ansys Fluent Tutorial Guide 2012) is used (see Figure 3). It is a simple sharp-edged orifice with the total length of 48 mm. The inlet diameter is equal to 23.04 mm and the outlet diameter is equal to 8 mm. The geometry is axisymmetric. The inlet is 16 mm in length. In the narrower part of the orifice, the average velocity increases and consequently, the static pressure decreases. It is in accordance with the Bernoulli's equation:

where *v _{i}* is inlet velocity (m/s),

*p*is inlet static pressure (Pa),

_{i}*g*is acceleration (m

^{2}/s) and

*h*is elevation of the point above a reference plane (m),

*v*is outlet velocity (m/s),

_{o}*p*is outlet static pressure (Pa). For that reason, cavitation is expected in these areas.

_{o}Boundary conditions are defined as follows: the inlet pressure is equal to 5 × 10^{5} Pa and the outlet pressure to 9.5 × 10^{4} Pa. The geometry of the analysed model is axisymmetric and two-dimensional (2-D) while the numerical calculations are performed for steady state. The chosen turbulence model is the k-*ɛ* turbulence model with standard wall functions. It is described using two parameters: turbulent kinetic energy with the value of 0.02 m^{2}/s^{2} and turbulent dissipation rate with the value of 1 m^{2}/s^{3}. The Reynolds number for the orifice is equal to 153,232. The detailed description of the setup and solution stage is presented in the *Cavitation Modeling* tutorial (Ansys Fluent Tutorial Guide 2012). During the implementation of the UDFs, in the solution methods, for pressure–velocity coupling a simple scheme is selected.

## METHODS IN NUMERICAL ALGORITHM VERIFICATION

The best way to determine the degree of compliance of the codes is comparison of the results obtained by numerical simulations. The most important results, which should be analysed, are distributions of cavitating areas with simultaneous consideration of its intensity. To have a comprehensive picture of the work of the examined computer algorithm, it is necessary to consider the process of numerical calculations. This information gives the analysis of convergence of the iterative procedure.

*X*and

*Y*are images dimensions (equal for both compared images),

*A*and

_{xy}*B*are values of cavitation intensity for pixels with

_{xy}*x*and

*y*coordinates in images

*A*and

*B*(images being compared), respectively. This coefficient gives information on how correlated are intensities of cavitation in images being compared. Typical values of these coefficients are in the range of −1 to 1, where 1 is the total positive correlation, 0 is no correlation, and −1 is total negative correlation.

*x*and

*y*coordinates) of images are obtained as: The value of this coefficient should be interpreted just as values of the coefficient.

## RESULTS AND DISCUSSION

### Analysis of the convergence history of volume fraction on default interior

*et al.*2002; Zwart

*et al.*2004), including distinction in implementation method: interface and the UDFs (see Figure 5). The first remark concerns the number of iterations which should be made until the calculations achieve the required convergence status. For the calculations made using models from interface (Schnerr & Sauer 2001; Zwart

*et al.*2004), the number of iterations is larger than for calculations made using the same models implemented through the UDFs. For the third model (Singhal

*et al.*2002), both implementation methods require a large number of iterations to achieve the expected convergence status.

The aim of the analysis of the diagram is the assessment of the correctness of the prepared codes. The similar level of volume fraction at the final stage of numerical calculations is the basis for the decision that there is no reason to reject the UDFs for the considered pairs of cavitation models (i.e. implemented directly from interface and using the UDFs). In two cases, the Schnerr & Sauer (2001) and Zwart *et al.* (2004) models, this assumption is satisfied. For the third model, doubts appear. The maximum value of volume fraction on default interior for the calculations made using the cavitation model available from the interface is about three times smaller than for the calculations using the UDF. The shape of the convergence history of calculations using the Singhal *et al.* (2002) model from the interface also stands out from other convergence curves. Noticing this large difference allows expecting also a significant difference in the final vapour volume fraction in the orifice. Taking into consideration that in other calculations the convergence status at a similar level is achieved, it can indicate that in this algorithm an error can be involved. Taking an explicit stand on this issue requires further investigations.

### Analysis of the vapour volume fraction in the orifice

*et al.*2004) do not give any reasons to reject the UDFs. The total length of vapour cloud is similar for the results achieved in both analyses: using the UDFs and the models available direct from interface. The only difference is the intensity of vapour cloud, which can be analysed more precisely using differential images. In both models using the UDFs, close to the rapidly changing diameter, the intensity of the vapour cloud is larger than in models available from the interface (Figures 6 and 8). The intensity of the cloud in the results obtained with models from the interface is larger in both cases in the internal middle and outlet parts. In the case of the results of numerical simulations of cavitating flow using the Singhal

*et al.*(2002) model, the large difference in the convergence history of volume fraction reflects in the contours of the vapour volume fraction. The vapour cloud achieved in the numerical calculations using the model from the interface is shortened by half in relation to the results achieved using the UDF. Because of the large similarity between the results of other models, the results achieved using the UDF for the Singhal

*et al.*(2002) model can be considered as correct.

### Statistical analysis of the vapour volume fraction in the orifice

The above presented results are statistically interpreted using the correlation coefficients presented above under ‘Methods in numerical algorithm verification’. Table 1 shows obtained values of the correlation coefficient. Analysis of Table 1 shows that there are no significant differences between the obtained values for five out of six implementations. The achieved values are 0.956 for the Schnerr & Sauer (2001) model and 0.942 for the Zwart *et al.* (2004) model. It means the results are consistent in more than 90%, even more than 95% for the Schnerr & Sauer (2001) model. One exception is the interface implementation of the Singhal *et al.* (2002) model, which is also visible in Figure 10. The results are correlated only in about 50%. The results obtained with the Singhal *et al.* (2002) model using the UDF are very good related to the results obtained using the Schnerr & Sauer (2001) and Zwart *et al.* (2004) models. The compatibility level is more than 95%.

. | Schnerr & Sauer (2001) interface . | Schnerr & Sauer (2001) UDF . | Singhal et al. (2002) interface
. | Singhal et al. (2002) UDF
. | Zwart et al. (2004) interface
. | Zwart et al. (2004) UDF
. |
---|---|---|---|---|---|---|

Schnerr & Sauer (2001) interface | – | 0.956 | 0.474 | 0.960 | 0.989 | 0.949 |

Schnerr & Sauer (2001) UDF | 0.956 | – | 0.518 | 0.972 | 0.945 | 0.993 |

Singhal et al. (2002) interface | 0.474 | 0.518 | – | 0.495 | 0.465 | 0.519 |

Singhal et al. (2002) UDF | 0.960 | 0.972 | 0.495 | – | 0.955 | 0.971 |

Zwart et al. (2004) interface | 0.989 | 0.945 | 0.465 | 0.955 | – | 0.942 |

Zwart et al. (2004) UDF | 0.949 | 0.993 | 0.519 | 0.971 | 0.942 | – |

. | Schnerr & Sauer (2001) interface . | Schnerr & Sauer (2001) UDF . | Singhal et al. (2002) interface
. | Singhal et al. (2002) UDF
. | Zwart et al. (2004) interface
. | Zwart et al. (2004) UDF
. |
---|---|---|---|---|---|---|

Schnerr & Sauer (2001) interface | – | 0.956 | 0.474 | 0.960 | 0.989 | 0.949 |

Schnerr & Sauer (2001) UDF | 0.956 | – | 0.518 | 0.972 | 0.945 | 0.993 |

Singhal et al. (2002) interface | 0.474 | 0.518 | – | 0.495 | 0.465 | 0.519 |

Singhal et al. (2002) UDF | 0.960 | 0.972 | 0.495 | – | 0.955 | 0.971 |

Zwart et al. (2004) interface | 0.989 | 0.945 | 0.465 | 0.955 | – | 0.942 |

Zwart et al. (2004) UDF | 0.949 | 0.993 | 0.519 | 0.971 | 0.942 | – |

The coefficients showing correlations between the UDF and interface implementations are shown in bold.

*r*focuses on the shape of the cavitation cloud. The basis for this coefficient are images obtained by using Equation (18). Through this modification, the contour of the cloud is emphasized and the inner part of the cloud fades into the background. In other words, intensity of contours visible in Figure 12 is linked with dynamics of changes in borders of cavitation cloud. The final images for each implementation are shown in Figure 12. Because of the symmetric character of the obtained images, only the upper half of each image is shown.

_{s}Table 2 shows values of the correlation coefficient, obtained based on the analysis of images shown in Figure 12. The analysis of the cloud shape confirms the results of the correlation analysis with the use of the first coefficient. The best-correlated pair is the model of Schnerr & Sauer (2001) – 54%, the worst is the model of Singhal *et al.* (2002) – 26%. Similarly, like in the first analysis, the results obtained with the Singhal *et al.* (2002) model using the UDF are very good related to the results obtained using the Schnerr & Sauer (2001) model and the Zwart *et al.* (2004) model. Interesting is the fact that cloud shapes are better correlated when grouped by method of implementation, not when grouped by model. This leads to the interesting conclusion that the method of model implementation in Fluent plays a significant role. It is the best visible for the Schnerr & Sauer (2001) and Zwart *et al.* (2004) models implemented using the UDFs. In the analysis using the first correlation coefficient, the compatibility is 99% and using the second coefficient 86%.

. | Schnerr & Sauer (2001) interface . | Schnerr & Sauer (2001) UDF . | Singhal et al. (2002) interface
. | Singhal et al. (2002) UDF
. | Zwart et al. (2004) interface
. | Zwart et al. (2004) UDF
. |
---|---|---|---|---|---|---|

Schnerr & Sauer (2001) interface | – | 0.541 | 0.309 | 0.452 | 0.706 | 0.547 |

Schnerr & Sauer (2001) UDF | 0.541 | – | 0.266 | 0.535 | 0.501 | 0.857 |

Singhal et al. (2002) interface | 0.309 | 0.266 | – | 0.261 | 0.320 | 0.262 |

Singhal et al. (2002) UDF | 0.452 | 0.535 | 0.261 | – | 0.456 | 0.497 |

Zwart et al. (2004) interface | 0.706 | 0.501 | 0.320 | 0.456 | – | 0.480 |

Zwart et al. (2004) UDF | 0.547 | 0.857 | 0.262 | 0.497 | 0.480 | – |

. | Schnerr & Sauer (2001) interface . | Schnerr & Sauer (2001) UDF . | Singhal et al. (2002) interface
. | Singhal et al. (2002) UDF
. | Zwart et al. (2004) interface
. | Zwart et al. (2004) UDF
. |
---|---|---|---|---|---|---|

Schnerr & Sauer (2001) interface | – | 0.541 | 0.309 | 0.452 | 0.706 | 0.547 |

Schnerr & Sauer (2001) UDF | 0.541 | – | 0.266 | 0.535 | 0.501 | 0.857 |

Singhal et al. (2002) interface | 0.309 | 0.266 | – | 0.261 | 0.320 | 0.262 |

Singhal et al. (2002) UDF | 0.452 | 0.535 | 0.261 | – | 0.456 | 0.497 |

Zwart et al. (2004) interface | 0.706 | 0.501 | 0.320 | 0.456 | – | 0.480 |

Zwart et al. (2004) UDF | 0.547 | 0.857 | 0.262 | 0.497 | 0.480 | – |

The coefficients showing correlations between the UDF and interface implementations are shown in bold.

## CONCLUSIONS

This paper presents the results of verification of the codes (UDFs) intended for implementation of homogeneous cavitation models obtained using Fluent software. Considered are three cavitation models available direct from Fluent interface: Schnerr & Sauer (2001), Singhal *et al.* (2002) and Zwart *et al.* (2004) models. The degree of compliance of the codes is assessed using three methods. The first method compares the final vapour void fraction in the considered area without any quantitative methods, the second with using of the statistical analyses and the third compares the convergence of the iterative procedure. The following concluding remarks can be made:

Three chosen verification methods give similar results.

Generally, the number of iterations is larger for calculations made using interface models than for calculations made by using the same models implemented through the UDFs, i.e. the UDF models are less computationally demanding.

The best-correlated pair is the Schnerr & Sauer (2001) model. In the analysis, the use of the first correlation coefficient gives the compatibility on the level of 96%, while use of the second coefficient gives the level of 54%.

The results of the correlation analysis for the Zwart

*et al.*(2004) model are also correct. The compatibility using the first coefficient is 94% and using the second coefficient 48%.All models correlate very well with the others with the exception of the Singhal

*et al.*(2002) model implemented directly from the interface.A better correlation is achieved through grouping by the method of implementation than by models.

There are no reasons to reject the UDFs for the Schnerr & Sauer (2001) and Zwart

*et al.*(2004) models. In case of the UDF for the Singhal*et al.*(2002) model, where the analyses give low compatibility values, because of the large similarity between the results of other models, the results can be considered as correct.