A Numerical Study on the Shedding Frequency of Sheet Cavitation

VII International Conference on Computational Methods in Marine Engineering

MARINE 2017

M. Visonneau, P. Queutey and D. Le Touzé (Eds)

A NUMERICAL STUDY ON THE SHEDDING FREQUENCY OF SHEET CAVITATION

Abstract. The last decades there is a strong interest in predicting cavitation dynamics as it is a prerequisite in order to predict cavitation erosion. Industrial applications require accurate results in an acceptable time span and as a result there is a focus on large scale dynamics. In this paper the RANS equations are used to investigate the shedding frequency of sheet cavities in two-dimensional simulations. First a verification study is made for the NACA 0015 in 6 degrees angle of incidence. A grid sensitivity study is conducted in wetted flow and in steady (non-shedding) cavitating condition (σ=1.6). Then an investigation is conducted in order to capture the shedding frequency. The results show that only when a correction for turbulent viscosity at the cavity-water interface is used it was possible to capture the shedding frequency as found in other numerical studies. Furthermore, a validation study is conducted on a NACA66-312 α=0.8 for two different angles of attack. The obtained results are compared and validated with the experimental data from Leroux et al. They indicate that the 2D shedding frequency predicted by the numerical simulations is in good agreement with the frequency obtained in the experiment.

2 MATHEMATICAL MODELS

2.1 Governing equations

The equations solved are the Reynolds Averaged Navier-Stokes (RANS) equations, where each instantaneous quantity can be split into time-averaged and fluctuating components. An incompressible segregated flow model is selected solving the integral conservation equations of mass and momentum in a sequential manner combined with the SIMPLE pressure-velocity coupling algorithm. From the basic principles of conservation of mass and momentum, the governing equations are written as follows:

	${\frac {\partial \rho }{\partial t}}+\nabla \cdot \left(\rho {\mathit {\boldsymbol {u}}}\right)=$$0$	(1)
	${\frac {\partial }{\partial t}}\left(\rho {\mathit {\boldsymbol {u}}}\right)+\nabla \cdot \left({\mathit {\boldsymbol {\rho uu}}}\right)=$$-\nabla p+\rho f+\nabla \cdot \tau$	(2)

where u is the velocity tensor, ρ is the fluid density, ${\textstyle \,\nabla p}$ the pressure gradient, ${\textstyle \rho f}$ the exterior force density per unit mass and τ the viscous part of the stress tensor.

The multiphase model used is the Volume of Fluid (VOF) method. Within the VOF approach, the fluid is treated as a single continuum, assuming a no-slip condition between liquid and vapor phases, with varying properties in space according to its composition ${\textstyle {\,a}_{v}}$. The volume fraction of the components is determined from the condition:

${a}_{v}+{a}_{l}=1$

(3)

while density and viscosity are defined as:

$\rho ={a}_{l}{\rho }_{l}+{\alpha }_{v}{\rho }_{v}\quad \,\,and\quad \,\mu ={a}_{l}{\mu }_{l}+$${a}_{v}{\mu }_{v}$

(4)

2.2 Turbulence modeling

In this study the SST k-ω turbulence model developed by Menter [5] is used in order to fully resolve the boundary layer. This approach effectively blends a k-ε model in the far field with a k-ω model near the wall. Reboud and Delannoy [6] showed the important role of the re-entrant jet on the cavity break-off cycle. However, the use of this turbulence model leads to very strong turbulent viscosity in the cavity wake hindering the re-entrant jet formation. It is stated by Reboud et al. [7] that this effect, which is not representative of the real behaviour, has been analysed to be related to the hypothesis of homogeneous flow and its no-slip condition between the two phases. That no-slip condition behaves as an artificial increase of dissipation.

That problem has been treated by an empirical reduction of turbulence dissipative terms in the two-phase regions, by modifying the turbulent viscosity [7]:

${\mu }_{t}=f\left(\rho \right){C}_{\omega }{\frac {k}{\omega }};\quad f\left(\rho \right)=$${\rho }_{\nu }+{\frac {{\left({\rho }_{m}-{\rho }_{\nu }\right)}^{n}}{{\left({\rho }_{l}-{\rho }_{\nu }\right)}^{n-1}}};\quad \,n\gg 1$

(5)

where ${\textstyle {\rho }_{\nu }}$ is the vapor density, ${\textstyle {\rho }_{l}}$ the liquid density and ${\textstyle {\rho }_{m}}$ the mixture density. For the constant ${\textstyle n}$ a recommended value ${\textstyle n=}$$10$ has been used. This modification improves the numerical simulations by taking into account the influence of the local compressibility effects of the vapor/liquid mixture on the turbulent structure.

2.3 Cavitation modeling

The conservation equation that describes the transport of vapour is similar to the mass conservation for liquid and is described by:

${\frac {\partial {a}_{v}}{\partial t}}+\nabla \cdot \left({a}_{v}{\mathit {\boldsymbol {u}}}\right)=$${S}_{{a}_{v}}$

(6)

In equation (6), ${\textstyle {S}_{{a}_{v}}}$ represents the source of volume fraction of vapor. In order to close the system of equations formed by the RANS equations and the transport equation for the vapor, a cavitation model should be introduced for the source term of the volume fraction of vapor. The cavitation model used is the model proposed by Schnerr-Sauer (2000) based on a simplified Rayleigh-Plesset equation, which neglects the influence of bubble growth acceleration, as well as viscous and surface tension effects:

${\frac {dR}{dt}}=sign\left({p}_{\nu }-p\right){\sqrt {\frac {2\left|{p}_{\nu }-p\right|}{3{\rho }_{l}}}}$

(7)

Where ${\textstyle {p}_{\nu }}$ is the saturation pressure, ${\textstyle p}$ is the local pressure around the bubble and ${\textstyle {\rho }_{l}}$ is the fluid density. According to this rate the source term in equation (6) is defined as:

${S}_{{a}_{v}}={\frac {4\pi {R}^{2}{n}_{0}}{1+\left({\frac {4}{3}}\pi {R}^{3}\right){n}_{0}}}{\frac {dR}{dt}}\,$

(8)

3 CASE DESCRIPTION

Two hydrofoils are going to be calculated using RANS equations in StarCCM+. The one case is the NACA 0015, on which a verification study is conducted and the second one is the NACA66-312 α=0.8 as an interest to validate the results with experimental data. A velocity inlet boundary is used for the upstream flow and a pressure outlet is defined on the outlet boundary. The pressure on the outlet for every condition is derived by the cavitation number:

$\sigma ={\frac {{P}_{out}-{P}_{\nu }}{\rho {U}_{in}^{2}/2\,\,}}\,\,$

(9)

3.1 NACA 0015 computational domain

For the verification study the NACA 0015 hydrofoil is used at 6 deg angle of attack with a chord length ${\textstyle c=}$$200$ mm and a computational domain as used in the VIRTUE WP4 Workshop [8]. The length of the domain is extended 2 and 4 chord lengths upstream and downstream of the foil respectively and in the y direction the wall is 285 mm away from the center of the foil (1400x570).

Table 1: Grid features for the NACA 0015

Grid	# Cells	y+	Level
G1	25,059	1.2261	Coarse
G2	49,931	0.8407	Medium
G3	77,081	0.6808	Fine
G4	100,569	0.5829	Very fine

Figure 1: Grid visualization of the whole domain (top) the refinements around the foil (bottom left) and the prism layers (bottom right).

The 2D computational grid is derived from a 3D mesh using trimmed hexahedral cells with local refinements and prism layers along the walls (Figure 1). The grid generation process is driven by specifying a base mesh size, relative to which all the sizes are defined (cell size in various regions, prism layer near wall thickness etc.). Finer meshes of the same topology are then automatically created by just reducing the base size. The only parameter that was kept constant was the prism layer total thickness so as to include the boundary layer in any case. The geometrical similarity was controlled then by changing the number of the prism layers. A grid sensitivity study is conducted following the approach as described in [9] using four different grids as indicated in Table 1.

The nominal speed during the simulations is 6 m/s corresponding to a Reynolds number based on the foil chord length equal to 1.09x10⁶. A pressure that corresponds to two different cavitation numbers σ=1.6 and σ=1.0 is applied on the outlet resulting on a steady and an unsteady condition respectively. A no-slip condition is applied on the foil and a slip condition on the top and bottom wall (see Table 2). The number of inner iterations per time step, the time step and the order of temporal discretization are defined after a sensitivity study in the unsteady cavitating condition (σ=1.0). In the end, 40 inner iterations and a time step corresponding to Courant number 0.75 are selected using a second order temporal discretization scheme. The Courant number is defined as the product of the inlet velocity and the time step over the cell size in the x-direction:

$Courant={\frac {{V}_{in}\cdot \Delta t}{\Delta x}}\,\,$

(10)

Table 2: Boundary conditions and flow characteristics for the two hydrofoils

Boundary Conditions	NACA 0015		NACA66mod α=0.8
Velocity inlet	6 m/s		5.33 m/s
Angle of incidence	6 deg		6 deg	8 deg
Pressure Outlet	Pout =20.9 kPa (σ=1.0)	Pout =31.7 kPa (σ=1.6)	Pout =16.5 kPa (σ=1.0)	Pout =20.5 kPa (σ=1.28)
Turbulent Viscosity	1%		1%
Turbulent Viscosity Ratio	10		10
Reynolds Number	1.09x106		0.8x106
Water Temperature	24 °C		20 °C

3.2 NACA 66 computational domain

A second case has been used in order to validate the results with experimental data. The selected geometry is a modified NACA 66 series with a chord length ${\textstyle c=}$$0.150\,$ m, which can be referenced as “N1.1-mod. NACA 66(mod.)-312, ${\textstyle \,\alpha =}$$0.8$” tested by Leroux et al. [10]. The relative maximum thickness was ${\textstyle \tau =}$$12\%$ at ${\textstyle 45\%}$ from the leading edge and the relative maximum camber was ${\textstyle 2\%}$ at ${\textstyle 50\%}$ from the leading edge. The exact coordinates of the tested hydrofoil can be found in [11]. In the experiment the hydrofoil was tested in four different angles of incidence, however simulations only for two angles are performed (see Table 7), one for a low shedding frequency (6 deg) and one for a higher shedding frequency (8 deg).

The same approach as in the NACA 0015 case has been used in order to validate the numerical set-up. The same grid topology is applied with a small refinement around the hydrofoil due to the thinner shape (especially at the leading edge) leading to a coarse grid with 28,593 cells (Figure 2). A nominal free stream velocity of 5.33 m/s is applied to the inlet corresponding to a Reynolds number based on the foil chord length equal to 0.8x10⁶. A pressure that corresponds to cavitation number equal to σ=1.0 and σ=1.28 is applied on the outlet for angle of incidence 6 and 8 deg respectively. No-slip condition on the foil and slip condition on the walls are applied. The number of inner iterations per time step, the courant number and the order of temporal discretization are identical as in the NACA 0015 case. All the details for both foils can be found on Table 2.

Figure 2: Refinement visualization around the NACA66 foil (left) and the applied computational domain (right)

4 RESULTS

4.1 NACA 0015

The flow around the frequently used NACA 0015 hydrofoil is investigated in three conditions: wetted flow, steady cavitating flow (σ=1.6) and unsteady cavitating flow (σ=1.0). The results in wetted flow are compared with experimental data [8] and the ones in cavitating flow only with other numerical works [13-18].

4.1.1 Wetted flow

A grid sensitivity study is conducted in wetted flow for the drag (Figure 3) and lift (Figure 4) coefficients and the results are compared with experimental data available from the VIRTUE Workshop [8]. The results are shown in detail in Table 3. An error estimation has been made by using an approach with power series expansion proposed by Eca and Hoekstra [9]. The expansions are fitted to the data in the least-squares sense.

Figure 3: Drag coefficient for different grid density and comparison with the experimental value.

Figure 4: Lift coefficient for different grid density and comparison with the experimental value.

Table 3: Drag and Lift coefficient values for all the grids and comparison with experimental values (Δexp). The uncertainty (U_φ) of the values is also shown.

Grid	CD	Uφ	% Δexp	CL	Uφ	% Δexp
Experiment	0.014			0.658
G1 (Coarse)	0.01459	8.47%	4.20%	0.67018	2.40%	1.85%
G2 (Medium)	0.01435	3.48%	2.47%	0.66838	1.62%	1.58%
G3 (Fine)	0.01430	2.36%	2.16%	0.66838	1.69%	1.58%
G4 (Very Fine)	0.01432	3.05%	2.28%	0.66703	1.05%	1.37%

The results show a good agreement with the experimental data and the assessment of the uncertainty shows that there is a small sensitivity of the drag coefficient to the grid density, however already with the medium mesh (G2) the solution is quite consistent.

4.1.2 Steady cavitating flow

In the steady cavitating condition a steady sheet cavity is expected that covers approximately the 20% of the foil. The time step and the number of inner iterations per time step (I_N) are investigated in the coarse mesh by comparing the pressure distribution, the vapor volume fraction and the time history of the total vapor volume. A grid sensitivity study for the lift and drag coefficient is conducted for this condition as well (Table 4).

Table 4: Drag and Lift coefficient values in cavitating flow for all the grids and the computed uncertainty (U_φ).

Grid	CD	Uφ	CL	Uφ
G1 (Coarse)	0.01916	15.67%	0.6352	0.24%
G2 (Medium)	0.01849	6.02%	0.6361	0.68%
G3 (Fine)	0.01852	6.00%	0.6358	0.51%
G4 (Very Fine)	0.01838	3.61%	0.6348	0.35%

Table 5: Numerically obtained frequencies from various sources (6 deg angle of incidence and σ=1.0) [12].

Author	Vinlet(m/s)	Shedding Frequency
Koop [13]	12	≈24
Sauer [14]	12	≈11
Schnerr et al. [15]	12 12	≈11.18 (incompressible) ≈9 (compressible)
Oprea [16]	6	≈14
Hoekstra & Vaz [17]	6	≈15.4
Ziru Li [18]	6	≈11.4

Figure 5: Pressure distribution along the foil for different I_Nper time step (σ=1.6).

Figure 6: Vapor volume fraction along the foil for different time steps (σ=1.6).

Figure 7: Total vapor volume for different time steps with 5 and 100 inner iterations per time step.

Figure 8: Sheet cavity visualization obtained by the coarse (top) and the very fine mesh (bottom).

It is concluded that a time step according to a courant number even higher than one can be sufficient (see Figure 6) and the number of inner iterations does not affect the development of the cavity although the vapor volume might not be fully converged in every time step (for instance despite the fact that using a time step corresponding to courant number 0.75 the total vapor volume is converged after 20 iterations per time step, the results between 5, 20 and 100 inner iterations are identical, Figure 5). The uncertainty assessment is in line with the results in wetted flow showing dependence to the mesh (stronger this time) for the drag coefficient. In addition to that, a slight grid sensitivity regarding the shape at the trailing edge of the sheet cavity is observed (see Figure 8).

4.1.3 Unsteady cavitating flow

An investigation on the impact of the number of inner iterations per time step is conducted first. Using 100 inner iterations and a time step corresponding to Courant number 0.75 the solution shows that a number of 40 inner iterations per time step is sufficient (Figure 9). An unsteady periodic cycle is predicted giving a shedding frequency of about 3.6 Hz (Figure 11). However, according to other numerical studies such a frequency seems to be very low (Table 5). To this end the modification for the turbulent viscosity is applied. The results show that a higher frequency can be achieved with a second order temporal discretization scheme or with lower time step. Eventually a shedding frequency of about 13.6 Hz is computed, using 40 inner iterations, courant number 0.75 and a second order temporal discretization scheme (Figure 11) together with Reboud’s correction for eddy viscosity.

The instantaneous images of the volume fraction are shown in Figure 10. First, as the cloud cavities from the previous cycle are moving downstream, a sheet cavity starts to grow at the leading edge in combination with some cavities growing at the trailing edge (steps 1-3). The re-entrant jet is formed and moves towards the leading edge as the bubbly cloud from the previous cycle collapses (steps 4-5). Then the sheet cavity starts to shed (step 6) and as it becomes smaller and smaller it continues shedding (steps 6-8) till it completely disappears (step 9). Finally all the shedding parts are combined into a bubbly cloud and move downstream to the trailing edge as the sheet cavity starts to grow again at the leading edge and the new cycle starts (steps 10-12).

A grid sensitivity study for the shedding frequency has also been conducted. The results are shown in Table 6. With every mesh a shedding frequency between 13 and 14 Hz has been computed. The high uncertainty of the solution can be explained by the unsteadiness and randomness of the shedding.

Figure 9: Convergence of the residuals and the total vapor volume in every time step.

Figure 10: Total Vapor volume, drag and lift coefficient during a typical shedding cycle.

Figure 11: Total vapor volume in time and frequency domain for NACA 0015 with (right) and without (left) Reboud’s correction for eddy viscosity.

4.2 NACA66 (mod.)-312

The NACA66 hydrofoil as described before is used to validate the computational set-up. The results are compared with experimental data obtained by Leroux et al. for two different conditions. The experimental and the computational obtained data are shown in Table 7.

Table 6: Grid sensitivity study on the shedding frequency and assessment of the uncertainty for each mesh for the NACA 0015.

Grid density	Shedding (Hz)	Uφ
G1 (Coarse)	13.59	15.40%
G2 (Medium)	13.92	22.44%
G3 (Fine)	13.50	12.32%
G4 (Very Fine)	13.23	7.39%

Table 7: Experimental frequency and Strouhal number based on chord length as measured by Leroux et al. V_ref=5.33m/s, Re = 0.8 x 10⁶[10].

	Experiment			CFD
α	f (Hz)	σ	stc	f (Hz)	σ	stc
5.5	2.88	0.88	0.081	-	-	-
6.0	3.50	0.99	0.098	3.60	1.00	0.101
7.0	4.50	1.13	0.127	-	-	-
8.0	18.00	1.28	0.507	18.29	1.28	0.515

Figure 12: Total vapor volume in time and frequency domain for 6 deg (top) and 8 deg (bottom) angle of incidence.

Figure 13: Experimental-numerical comparison on the NACA66 for 6 deg angle of incidence. Experimental images are computed from an average of three instantaneous periods (Δt=1/50 s) and compared with computed instantaneous void fraction images with the same period.

The results show that in both cases a frequency similar to this in the experiment is obtained. There is a difference less than 3% for the low frequency case and 2% for the high frequency case. A comparison between the computations and the experimental data of the foil in 6 deg angle of incidence is shown in Fig. 13. As illustrated by Leroux et al. two steps can be identified during a typical shedding cycle: The first step consists of the growth of the sheet cavity (Figure 13 a-e) till it is slowed down and counterbalanced by the shedding of vapor structures (secondary clouds) in the wake (Figure 13 f-i).After the shedding of secondary clouds, the detachment of a large vapor cloud (main cloud) occurs (Figure 13 j). It is followed by the roll-up and convection of the main cloud (Fig. 13 k) together with the growth of the residual cavity. The second step occurs just after the cavity break-off. Indeed, the growth of the residual cavity is abruptly stopped at nearly the sametime the main cloud of vapor collapses (Figure 13 l), and the residual cavity almost entirely disappears (Figure 13 m). Then the cavity starts to grow again.

Similar cavitation dynamics are calculated by the simulations. The growth of the cavity and the secondary clouds are captured as well as the detachment of the large vapor cloud and the sudden vanishing of the cavity after the collapse of the cloud. Discrepancies can only be observed on the growth of the residual cavity, where a larger expansion of the residual cavity is predicted in the computations (the same behaviour was also predicted in the computations by Leroux et al).

5 CONCLUSIONS AND RECOMMENDATIONS

In this study an attempt was made to verify the incompressible RANS solver in StarCCM+ in cavitating flow. Despite the three-dimensional character of cavitation dynamics a first investigation was conducted on the grid and numerical (time step, inner iterations etc.) sensitivity with the intension to predict the shedding frequency using two-dimensional domain. For the current computational set-up and the tested conditions the following conclusions are drawn:

• When a steady sheet cavity is predicted the effect of the time step and the number of the inner iteration on the results are negligible. However, the grid density had a slight impact on the shape of the sheet cavity.

• In the unsteady condition, on the other hand, the time step and the number of inner iterations per time step seem to play an important role on the prediction of the shedding frequency. A higher frequency was captured only when the correction for the turbulence viscosity in areas with higher vapor volume was applied. Without the correction the effect of the re-entrant jet could not be captured thoroughly, leading to a “delayed” shedding and consequently to a lower shedding frequency.

• Furthermore, it should be noted that the number of iterations needed per time step changes for different time steps and order of temporal discretization, so it is suggested that they are selected in such a way that convergence of the total vapor volume per time step is achieved.

• A grid independent solution has been reached; even the coarsest mesh was capable of capturing the dynamic shedding in a high frequency after application of Reboud’s eddy viscosity correction.

• As a second step, an effort to validate the model was made comparing the numerical results with experimental data. Good agreement was obtained and the shedding frequency was accurately predicted. Discrepancies can only be observed in the maximum total volume per cycle.

• Further computations on a three-dimensional domain are recommended to investigate possible alterations on the cavitation dynamics and the shedding frequency.

• It is finally recommended to investigate the possible erosion mechanisms and the capability of predicting the potentially erosive cavitation implosions using incompressible URANS solver.

REFERENCES

[1]	M. Bijlard and N. Bulten, "RANS simulations of cavitating azimuthing thrusters," in Fourth International Symposium on Marine Propulsors, Austin, Texas, USA, (2015).
[2]	G. Bark, N. Berchiche and M. Grekula , "Application of principles for observation and analysis of eroding cavitation - The EROCAV observation handbook," Edition 3.1, (2004).
[3]	R. Fortes-Patella, J. L. Reboud and L. Briancon-Marjollet, "A phenomenological and numerical model for scaling the flow aggressiveness in cavitation erosion," EROCAV Workshop, Val de Reuil, (2004).
[4]	T. J. C. V. Terwisga, P. A. Fitzsimmons, Li Ziru and E. J. Foeth, "Cavitation Erosion - A review of physical mechanisms and erosion risk models," in Proceedings of 7th International Symposium on Cavitation, CAV2009, Ann Arbor, Michigan, USA, August 2009.
[5]	F. Menter, "Two-equation eddy-viscosity turbulence modeling for engineering applications," AIAA Journal, vol. 32, pp. 1598-1605, (1994).
[6]	J. L. Reboud and Y. Delannoy, "Two-phase flow modelling of unsteady cavitation," in 2nd Int. Symposium on Cavitation , Tokyo , (1994).
[7]	J.-L. Reboud, B. Stutz and O. Coutier, "Two-phase flow structure of cavitation: experiment and modelling of unsteady effects," in Third International Symposium on Cavitation, Grenoble, France, (1998).
[8]	European Comission Project FP6, http://www.virtual-basin.org/: Virtual Towing tank Utility in Europe, (2005-2008).
[9]	L. Eca and M. Hoekstra, "A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies," Journal of Computational Physics, vol. 262, pp. 104-130, (2014).
[10]	J.-B. Leroux, O. Coutier-Delgosha and J. A. Astolfi, "A joint experimental and numerical study of mechanisms associated to instability of partial cavitation on two-dimensional hydrofoil," Phys. Fluids 17, 052101, (2005).
[11]	J.-B. Leroux, J. A. Astolfi and J. Y. Billard, "An Experimental Study of Unsteady Partial Cavitation," ASME J. Fluids Eng., vol. 126, pp. 94-101, (2004).
[12]	Z. Li, M. Pourquie and T. Terwisga, "Assessment of Cavitation Erosion With a URANS Method," Journal of Fluids Engineering, Vols. 136, 041101, (2014).
[13]	A. H. Koop, "Numerical Simulation of Unsteady Three-Dimensional Sheet Cavitation," Ph.D. thesis, University of Twente, Enschede, The Netherlands, (2008).
[14]	J. Sauer, "Instationär kavitierende strömungen - Ein neues modell, basierend auf fron capturing (VoF) und blasendynamik," Ph.D. Thesis, Karlsruhe University, Karlsruhe, Germany, (2000).
[15]	G. H. Schnerr, S. J. Schmidt, I. H. Sezal and M. Thalhamer, "Shock and Wave Dynamics of Compressible Liquid Flows With Special Emphasis on Unsteady Load on Hydrofoils and Cavitation in Injection Nozzles," in Proceedings of 6th International Symposium on Cavitation, Wageningen, The Netherlands, (2006).
[16]	A. Oprea, Prediction of Tip Vortex Cavitation for Ship Propellers, PhD Thesis, University of Twente, Enschede, The Netherlands: Uitgeverij POXPress, (December 2013).
[17]	M. Hoekstra and G. Vaz, "FreSCo Exercises for NACA0015 Foil," VIRTUE WP4 Workshop, Contribution from MARIN, (2008).
[18]	Z. Li, Assessment of Cavitation Erosion with Multiphase Reynolds-Averaged Navier-Stokes Method. PhD Thesis, Technische Universiteit Delft, (2012).
[19]	J.-P. Franc and J.-M. Michel, Fundamentals of Cavitation, volume 76, Dordrecht: Kluwer Academic Publishers, (2004).