Application of a pressure based CFD code with mass transfer model based on the Rayleigh equation for the numerical simulation of the cavitating flow around a hydrofoil with circular leading edge

The most common method for simulating cavitating flows is using the governing flow equations in a form with a variable density and treats both phases as incompressible in combination with a transport equation for the vapour volume fraction. This approach is commonly referred to as volume of fluid method (VoF). To determine the transition of the liquid phase to vapour and vice versa, a relation for the mass transfer is needed. Several models exist, based on slightly differing physical assumptions, for example derivation from the dynamics of single bubbles or large bubble clusters. In our simulation, we use the model of Sauer and Schnerr which is based on the Rayleigh equation. One common problem of all mass transfer models is the use of model constants which often need to be tuned with regard to the examined problem. Furthermore, these models often overpredict the turbulent dynamic viscosity in the two-phase region which counteracts the development of transient shedding behaviour and is compensated by the modification proposed by Reboud. In the presented study, we vary the parameters of the Sauer-Schnerr model with Reboud modification that we implemented into an OpenFOAM solver to match numerical to experimental data.


Introduction
Cavitation is a phenomenon which is still not completely understood and is, in turn, a topic of ongoing research.A simplified understanding of cavitation and a common modelling approach that is also used in the presented study is that cavitation occurs when the pressure in a flow drops below the vapour pressure of the liquid, leading to vaporisation and formation of a cavitation pocket.When this cavitation pocket collapses due to increasing pressure, microjets form which lead to sharp pressure peaks on adjacent surfaces.By this effect, cavitation can lead to erosive damage on e.g.pump impellers, turbines, ship propellers or in valves, but is also a source for unwanted effects like increased noise or vibration.However, cavitation can also be put to use to control the flow rate in fuel system injectors.
To model cavitation, several approaches exist.A wide spread approach -as it is computationally efficient and numerically robust -is the volume of fluid method (VoF), which adds a transport equation for the vapour void fraction to the set of governing equations.As the original volume of fluid method was developed to describe the transport of the vapour phase in conjunction with an interface reconstruction approach, but did not include inter-phase mass transfer (which is needed in case of cavitation), extensions to this model were proposed to include source and sink terms in the transport equation, for example by Sauer and Schnerr [1].
In the presented study, we summarise these extensions and show how to apply the model to simulate the cavitating flow around a hydrofoil.We make an effort to adjust the model parameters in a way that the simulation results match the shedding frequency of the experimental data.
a Corresponding author: christian.deimel@rub.de In former investigations, the same hydrofoil test case as investigated in the presented study has been simulated by Reboud [2] and Frobenius [3,4].However, Reboud uses a different approach to model the cavitating flow by employing a simple barotropic equation of state (density is directly dependant on pressure and vice versa) to obtain the mixture density.In the study of Reboud [2], two operation points are simulated, with results shown to be in good agreement to experimental values.
Frobenius, on the other hand, uses the Sauer-Schnerr mass transfer model for inter-phase mass transfer, but focusses on the evaluation of re-entrant jet velocities in [3].In [4], shedding frequencies for two different Reynolds numbers are evaluated, however it does not become clear how these frequencies were obtained from the numerical data.
Furthermore, studies applying OpenFOAM solvers and the Sauer-Schnerr mass transfer model to cavitating flows exist, but focus on different hydrofoils or internal flow test cases like a cavitating venturi nozzle.As a notable example we want to mention the study by Gosset et al. [5] who apply the Sauer-Schnerr mass transfer model with Reboud correction to the cavitating flow around a NACA0015 hydrofoil.
On the basis of the former studies we select the Sauer-Schnerr mass transfer model [1] for our investigation and adjust the model parameters to match experimentally observed frequencies for several operating points.Then, the model parameters can be applied to more complex flow problems like the simulation of hydraulic machinery where knowledge about the cavitation behaviour is of prime interest to determine the admissible operating range. 2 Numerical model

Two-phase flow modeling
A common assumption to model two-phase flow is not to model two distinct phases, but to describe the bulk liquid as a homogenous mixture of a liquid phase and a dispersed vapour phase.By this, the properties of the mixture can be related to the volume fraction of the vapour phase, leading to an equation for the mixture density and mixture dynamic viscosity: On a side note, equation 2 represents a strong simplification of the real physical behaviour but is a common assumption in two-phase flow modelling.The mixture has to satisfy the continuity equation, hence: By substituting equation 1 into equation 3, we obtain the non-conservative form of the continuity equation (i.e. the velocity field is no longer divergence free): Additionally, a transport equation for the vapour volume fraction is needed: In this equation, ṁ denotes the rate of inter-phase mass transfer.Again, by substituting equation 4 into equation 5, this can be rewritten as:

Sauer-Schnerr mass transfer model
To describe the vapour volume fraction α, Sauer and Schnerr [1] model the bulk liquid as a mixture of a liquid phase containing a large number of spherical bubbles: To derive an equation for the change in the vapour volume fraction, the volume change of the bubbles contained in a volume element is used [6], leading to To model the bubble growth, the Rayleigh-Plesset equation is used: By neglecting higher order derivatives, viscous terms and surface tension, this equation can be simplified to By assuming the bubble pressure to be equal to the saturation pressure of the liquid, applying the chain rule to equation 8 and substituting equation 10, this leads to Finally, substituting equation 11 into equation 6, the mass transfer rate for the Sauer-Schnerr model can be written as p − p sat ρ l for p > p sat (12) with the bubble radius R b derived from equation 7 To account for vaporisation and condensation to happen on different time scales, scaling factors F vap and F cond can be introduced into equation 12.For the original Sauer-Schnerr model, these model parameters equal For the initial number of bubbles per unit volume and the initial bubble radius Sauer and Schnerr chose values of In section 5, we will show how the model parameters can be modified to match simulation results to results obtained from experiments.

Turbulence model and parameters
The k − ω − S S T model as proposed by Menter [7] is used in our simulations.Freestream turbulence prescribed at the inlet is set to Tu = 5 %, so the turbulent kinetic energy can be calculated with the relation With k known, a value for the turbulent dissipation rate ω at the inlet can be calculated by using a characteristic length L, In our simulations, the characteristic length is set to an arbitrarily chosen value of L = 3.88 × 10 −3 m.

Turbulent viscosity modification
It is known from several studies (e.g.[8,9]) that RANS models are unable to correctly predict the transient shedding behaviour of cavitation without further modification, which is attributed to an overprediction of turbulent viscosity in the cavity closure region [10].Therefore, the formation of a re-entrant jet does not occur, leading to the cavity remaining attached.As the resulting quasi-steady behaviour is clearly unphysical in the regime of cloud cavitation that is considered in the present study, a simple correction of the turbulent viscosity has been suggested [8] to overcome this problem by defining the turbulent dynamic viscosity as a function of the mixture density: In this expression, f (ρ m ) is defined as Common choices for the exponent in this equation which are deemed to be sufficient for most applications are n = 3 or n = 10.After carefully comparing the results obtained with either choice, no influence of the parameter could be seen in the results, so for this study n = 3 was used in all numerical simulations.The hydrofoil investigated in our simulation is a circular leading edge-hydrofoil (CLE) with a chord length of l chord = 107.9mm and a thickness of d = 16 mm.The foil has been subject of several studies and PhD theses at the TU Darmstadt (see e.g.[11][12][13]), mainly focussing on the impact of cavitation erosion on different materials.

Experimental validation data
In the above mentioned studies, measurements were conducted with the foil mounted in the rectangular test section of a cavitation tunnel in a closed loop, and the experimental setup and results of [11][12][13] are summarised here.The upstream flow velocity in the experiments was set to u ∞ = 16 m s −1 which corresponds to a Reynolds numer of Re = 1.6 × 10 6 based on the hydrofoil chord length.Cavitation behaviour was analysed for different angles of attack ranging from α = 0 • to 10 • .To analyse dynamic wall pressure, pressure tappings were placed at four distinct positions on the suction side of the hydrofoil, connected to four pressure transducers inside the foil.The position of these pressure tappings is shown in figure 2. Transforming the time signals from the pressure transducers into the frequency domain, the period of the cavitation shedding can be determined.
Following an idea of Reboud et al. [2], we adopt this method for our numerical simulations.In the experiments, transient, fluctuating cloud cavitation could be observed for the considered operation conditions listed below.
During a shedding cycle, a re-entrant jet forms and leads to seperation of the cavity.The cavity is then transported downstream where it recondenses in areas of higher pressure.The whole process repeats itself, beginning with the formation of a new cavity at the leading edge of the foil (quasi-periodic behaviour) [11].The cavitation number was varied in the range of σ = 2 < σ < σ i in the experiments.Besides conducting measurements for different operation points based on the cavitation number, the angle of attack and freestream velocity of the flow were varied as well.Primarily, dynamic wall pressures were measured using the aforementioned pressure transducers, but visual analysis of the cavitating flow for selected operating points was conducted as well by means of high speed imaging.An exemplary cavitation cycle captured with this method is shown in figure 3. The extension of the computational domain corresponds to the channel height of h = 0.1 m as in the experiment.In accordance to the experimental data used for validation, the hydrofoil is placed in the computational domain at an angle of attack of 5 • .We created a block structured grid consisting solely of hexaeder cells.The grid features two O-grids wrapped around the hydrofoil to ensure good element quality in the critical region, see figure 4 for details.The fine grid G04 was generated by doubling the number of nodes per edge of each block of the coarser grid G02.The resolutions of these grids are shown in table 1.The y + values for the grid G02 are in the range of 7 < y + < 184 and 4 < y + < 101 for the fine grid G04 respectively.Although cavitation is a phenomenon that is generally three-dimensional, we know from preliminary studies of the cavitating flow around a NACA0015 hydrofoil that a two-dimensional simulation is sufficient for obtaining the shedding frequency.Hence, the span direction is only resolved with one cell.

Physical properties
For the physical properties used in the simulations, we select the values for water at 20 • C. At this temperature, the saturation pressure of water is 2318 Pa.
As heat transfer is not included in the used model (i.e. the simulation is isothermal) and these properties are known to be only weakly dependant on the pressure in the examined region, they are kept constant throughout the simulation.

Boundary conditions
The outlet pressure is prescribed by using a relation for the cavitation number σ, which is defined as The resulting outlet pressures for the cavitation numbers which were examined are listed in table 2. The inlet velocity is prescribed with a constant value of u ∞ = 13 m s −1 , which corresponds to a Reynolds number of Re = 1.3 × 10 6 based on the chord length.Note that this is different from the Reynolds number at which the measurements were conducted (Re = 1.6 × 10 6 ), however we know from [4,13] that there is only a weak dependency between Reynolds number and obtained shedding frequency for the considered test case.
The upper and lower walls of the simulation domain are modeled as slip walls, while at the hydrofoil surface a no-slip condition is employed.
For the turbulence quantities k, ω and ν t the default wall functions provided by OpenFOAM are used.For ω, this includes a continous blending approach between the viscous sublayer and the log-layer of the velocity profile close to the wall.
Regarding the parameter set for the Sauer-Schnerr mass transfer model, we chose a value of n 0 = 1 × 10 8 m −3 for the number of bubbles per unit volume, while the initial bubble radius was set to R b = 2.673 × 10 −5 m.

Numerical method
For our simulations, we use a solver derived from inter-PhaseChangeFoam included in OpenFOAM 2.1.x 1 .Inter-PhaseChangeFoam is a solver for two incompressible, isothermal immiscible fluids with phase-change.To solve the governing equations, this solver employs a coupled PISO-SIMPLE algorithm.The mass transfer models included in the standard distribution are the models by Merkle [14], Sauer and Schnerr [1], and Kunz [15], while only model [1] is utilised in the present study.However, additional models can be included easily by providing the appropriate source terms.For the implementation of the modification of the dynamic viscosity according to [8] and described in section 2.3.2, a re-implementation of the diffusive terms in the momentum equation has been done.
The convective terms in the momentum equation are discretised using the second order linearUpwind scheme.For time integration, an implicit first order Euler scheme is used due to stability reasons.To counterbalance the inaccuracies introduced by this scheme, a maximum CFLnumber of 0.2 is imposed, leading to time step sizes in the magnitude of 1 × 10 −6 s.
The volume fraction equation is discretised using the van Leer scheme [16] without interface compression [17] despite the fact that Schnerr and Sauer [1] recommend a first order upwind discretisation.For the convection of the turbulent kinetic energy and specific turbulent dissipation rate, the first order upwind method is used.Although second order discretisation is desirable for the turbulence quantities as well, this choice again was made due to stability reasons.

Frequency evaluation
The simulations have been run long enough that a quasiperiodic shedding develops which allows a reliable frequency analysis as well as frequency results that are independent of the simulated physical time.We employ a method analogous to the experimental strategy and sample the pressure at four distinct locations close to the hydrofoil surface.To obtain a shedding frequency from the time signals, a fast fourier transform of 1200 data points with a sampling frequency of f samp = 4000 Hz is performed.Furthermore, the FFT is used in combination with the Hann window function [18] to decrease spectral leakage.

Mass transfer model parameter variation
As we try to match our simulation results to the experimental data, we systematically vary the values for the constants F vap and F cond of the mass transfer model.Starting from the initial values of F vap = F cond = 1, we use combinations of powers of two in a range between 1 and 32 for both constants, varying one of them while the other one is held constant at its respective value.For each tested parameter set, the shedding frequency is evaluated with the method described.Most of these parameter sets either lead to quasi-stationary behaviour (i.e. the cavity remains attached) or do not show a distinct shedding frequency in the FFT analysis.The best results were obtained with a combination of F cond = 4 and F vap = 16, showing transient shedding behaviour with formation of a re-entrant jet and distinct frequency peaks in the range of the experimental values.

FFT frequency analysis
Examining the frequency spectra obtained from the simulations in detail, we find that in some cases more than one strong peak can occur.For σ = 2, a strong peak in the range of f ≈ 20 Hz exists.By visually examining the results , the actual shedding frequency is evaluated to be f ≈ 60 Hz, which corresponds to the second peak in the frequency spectrum for grid G04.We take into account the frequency analysis results of the higher values of σ and assume a continuous shape of the frequency curve over σ.Therefore, and due to the presence of the peak at f ≈ 60 Hz on the finer grid G04, we assume that this is the frequency that should be compared with the experimental data.The source of the strong low frequency component on both grids remains unclear.
For σ = 2.3, several peaks of almost equal amplitude exist.Again, by visually evaluating the shedding behaviour the first peak at f ≈ 110 Hz seems to be most feasible and is assumed to correspond to the experimental data.One observed shedding cycle is depicted in figure 9.
For σ = 2.5 results from both grids show good agreement, with a shedding frequency of f ≈ 140 Hz.The smaller peaks at a frequency of f ≈ 70 Hz are not feasible since we know from the experimental data that the shedding frequency increases with the cavitation number, hence we conclude that this is not the primary frequency.
For σ = 2.7, we can see two peaks in the spectrum for grid G02.Again, since the higher frequency peak occurs at the finer grid G04, and we consider the fine grid results to more reliable than the coarse grid results, we assume that the higher frequency is the shedding frequency in this case.
Shown in figure 10 are the results obtained from our simulations, together with the experimental data taken from the thesis of Hofmann [11].It can be seen that shedding frequencies are underpredicted by our model con- stants.However, the trend of continously increasing frequency with increasing σ in correspondence to the experiment was found.
All in all, our results not only emphasise the necessity of a study on different grids, but also the importance of not relying on only one tool (in this case the FFT) for evaluation.A visual examination of the results remains important but is impractical for parameter studies with a wider range of parameter variation.Integral values like the void fraction, integrated over the computational domain or the time evolution of the lift and drag coefficients could be included in the evaluation as well.

Conclusion and outlook
As we have shown in the previous section, we successfully applied the presented method to obtain shedding frequencies from results of numerical simulations.The Sauer-Schnerr mass transfer model in combination with the Reboud modification is deemed to be appropriate for the prediction of transient shedding behaviour in this kind of flow problem.Also, it has to be noted that we found the flow around this hydrofoil in particular to be more intricate in terms of evaluation than for example the cavitating flow around a NACA0015 hydrofoil, where good results were achieved in preliminary studies from solely evaluating the integrated vapour volume fraction in the computational domain with no need to apply more complex methods like a FFT.Although we did not reach the exact measurement data with our parameter set, this does not invalidate the frequency analysis method in general.Clearly, shedding frequencies are underpredicted by our model constants, but the general trend is in good agreement with the experimental data and leaves room for future improvement.Manual parameter variation helped to gain insight into the interaction of the model constants but remains a tedious task.To continue our work with mass transfer rate based models, an optimisation strategy that can be used for a sensitivity analysis on a larger scale (then, of course, automatised) has to be developed.As an objective function, a more refined frequency evaluation shall be considered, supported by using time averaged lift coefficient along the hydrofoil surface.
In conclusion, we are confident a better parameter set may be found with the help of the presented method which will not only suit this hydrofoil, but can be used for the simulation of more complex problems as well.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences,

Figure 2 .
Figure 2. Placement of pressure tappings along the suction side of the hydrofoil; distances are given in millimetres (redrawn according to [2])

Figure 4 .
Figure 4. Detailed view of the O-Grid around the hydrofoil

Table 1 .
Resolution of the computational grids used in the presented study

Table 2 .
Outlet pressure based on cavitation number