Numerical Investigation of the Cavitation Effects on the Vortex Shedding from a Hydrofoil with Blunt Trailing Edge

Abstract

Vortex cavitation can appear in the wake flow of hydrofoils, inducing unwanted consequences such as vibrations or unstable behaviors in hydraulic machinery and systems. To investigate the cavitation effects on hydrofoil vortex shedding, a numerical investigation of the flow around a 2D NACA0009 with a blunt trailing edge at free caviation conditions and at two degrees of cavitation developments has been carried out by means of the Zwart cavitation model and the LES WALE turbulence model which permits predicting the laminar to turbulent transition of the boundary layers. To analyze the dynamic behavior of the vortex shedding process and the coherent structures, two identification methods based on the Eulerian and Lagrangian reference frames have been applied to the simulated unsteady flow field. It is found that the cavitation occurrence in the wake significantly changes the main vortex shedding characteristics including the morphology of the vortices, the vortex formation length, the effective height of the near wake flow and the shedding frequency. The numerical results predict that the circular shape of the vortices changes to an elliptical one and that the vortex shedding frequency is significantly increased under cavitation conditions. The main reason for the frequency increase seems to be the reduction in the transverse separation between the upper and lower rows of vortices induced by the increase in the vortex formation length.

1. Introduction

In recent decades, the interaction effects between cavitation and vortical flow structures have attracted great attention in academic research and engineering practice. The high vorticity intensity in the center of a vortical flow generates a low-pressure region which is prone to cavitation inception [1,2]. Furthermore, this type of cavitation is affected by vorticity generation mechanisms [3]. Once a cavitation bubble is formed and develops inside the vortex core, the unsteady vortex dynamic behavior, e.g., vortex shedding frequency and vortex formation length, tend to change significantly [4]. As cavitation is further developed, it may give rise to vibrations, noise and material erosion which is an actual risk for the safe operation of hydraulic machinery and systems [1].

To distinguish vortical flow cavitation and the cavitation appearing in a shear flow induced by a high-speed submerged jet, so-called jet cavitation, the cavitation in the vortices generated in the wake of the flow around a bluff body wake is called wake cavitation [2]. Many experimental and numerical research works can be found treating wake vortices behind bluff bodies in single-phase flows, like behind a circular cylinder [5,6,7] and a 2D wedge [8,9], which permit explaining the alternate vortex shedding mechanism. Although some researchers have also paid attention to the cavitating flow behind a circular cylinder [10,11,12,13] and a 2D wedge [3,4,14,15,16,17,18], there is a lack of investigations devoted to understanding the cavitation effects on the dynamics of wake vortices, especially behind 2D symmetric hydrofoils.

Young and Holl [14] were the first to experimentally measure the vortex shedding frequency behind a 2D triangular wedge. They found that the vortex shedding frequency remains constant when cavitation appears at cavitation number σ_i, but as the cavitation number σ further decreases below σ_i and the cavitation size increases, the vortex shedding frequency suffers an abrupt rise near ½ of σ_i. After the vortex shedding frequency reaches this maximum value and σ is further decreased, the individual cavitation vortices grow until the point that they suffer a coalescence and it becomes impossible to measure any shedding frequency. Rao and Chandrasekhara [11] and Ramamurthy [12,15,19] also observed and confirmed similar results both in a circular cylinder and triangular wedge geometry under different area ratios of the model to the test section. Based on high-speed imaging results, Belahadji [4] found that cavitation cannot be considered as irrelevant for the flow around a 2D wedge since, as soon as cavitation happens, it can dramatically change the flow dynamics and lead to an increase in the vortex shedding frequency. However, this behavior is not often mentioned in research devoted to the observation and characterization of wake cavitation behind 2D symmetric hydrofoils. Ausoni et al. [20] conducted an experimental research to determine the cavitation effects on the flow behind a truncated NACA 0009 hydrofoil at high Reynolds numbers. With a high-speed camera and a laser vibrometer, they confirmed that the vortex shedding frequency increased gradually with the development of cavitation, although no mechanism for the phenomena was thoroughly described.

An important parameter that determines the vortex shedding frequency at the high Reynolds numbers is the vortex formation length as proposed by Gerrard [5]. In a single-phase flow, Gerrard insisted that the larger the vortex formation region, the lower the vortex shedding frequency. However, when cavitation occurs, the situation becomes more complex. For example, Kumar et al. [21] observed and confirmed similar results at subcritical Reynolds numbers on the cavitating flow over a circular cylinder employing high-speed imaging. In contrast, Ramamurthy et al. [12] observed an increase in the vortex formation length when σ decreased in both a 2D wedge and a circular cylinder. Despite those works, few studies can be found related to the vortex formation length behind hydrofoils. Ausoni et al. [20] provided measurements of the vortex formation length for different free stream velocities at cavitation free regimes. They concluded that the vortex shedding frequency is reduced as the vortex formation length increases both in the natural boundary layer transition and tripped transition setups, which agrees well with Gerrard [5]. Nevertheless, no measurements of the vortex formation length for cavitation conditions were reported in their study.

Some numerical investigations regarding the vortex shedding flow around hydrofoils without cavitation should be mentioned. Ausoni et al. [20] carried out a numerical investigation of the vortex shedding behind a NACA 0009 hydrofoil with a truncated trailing edge for several free stream velocities and high Reynolds numbers under no cavitation conditions. Their numerical results showed that the boundary layer transition from laminar to turbulent regimes predicted with the SST transition turbulence model provided more consistent results with the experimental measurements than the standard SST turbulence model. Lee et al. [22] and Chen et al. [23] used 2D LES and 3D LES turbulence models, respectively, to investigate the same hydrofoil geometry. They stated that the choice of a turbulent model with the capability of predicting the boundary layer laminar to turbulence transition is a key point to obtain an accurate prediction of the vortex shedding frequency.

Under cavitation conditions, the most commonly used cavitation model is the homogeneous mixture model, which is derived from the simplified Rayleigh–Plesset bubble dynamics equation. The coupling of the cavitation model with RANS and LES turbulence models has been found to capture with good accuracy the cavitating flow phenomena in the wake of bluff bodies. For example, Kim et al. [24] used a combination of URANS with a homogeneous mixture cavitation model to investigate the unsteady cavitating flow over a 2D triangular wedge. Further, Gnanaskandan et al. [13,17] applied the LES turbulence model coupled with a compressible homogeneous mixture cavitation model to investigate the cavitating wake flow over a circular cylinder and a 2D wedge. In summary, these numerical investigations confirmed the capability of the homogeneous mixture cavitation models to predict the vortex shedding frequency and the vortex dynamics behind a bluff body. Futhermore, their results demonstrated that the patterns and processes of vortex shedding were well captured for both geometries. More specifically, it was found that LES can offer a better performance in the prediction of the intermittent vortices compared with RANS. Thanks to the latest works of the research groups at the Swiss Federal Institute of Technology of Lausanne [25,26,27,28], the Technical University of Munich [29,30,31], the Chalmers University of Technology [32,33,34], the Delft University of Technology [35] and the University of London [36,37,38], great progress has been made in the simulation of cavitation in vortical flows with LES. These works have increased the understanding of the cavitation effects on vortices/eddies dynamic behavior, ranging from the large scales to the intermittent ones. Based on their findings, the LES WALE model has been applied in the present simulations for its demonstrated performance when predicting the laminar–turbulence boundary layer transition and the cavitating vortical flows.

Although there is no agreed definition of what is a “vortex” among the fluid mechanics community [39], a vortex may be derived from a snapshot of the velocity field and its gradient under the Eulerian frame of reference. The most used criteria for vortex tracking are the Q criterion [40], the vorticity magnitude [41] and the ∆ criterion [42], among others. All of these Eulerian vortex criteria have been widely applied in fundamental vortex dynamic research [39] as well as in industrial engineering applications [2,4,43]. On the other hand, a vortex is just one type of “coherent structure” in which the fluid particles share a common orbit axis from the Lagrangian point of view. The term “coherent structure” can be first traced back to Brown and Roshko [44]. In 1991, Robinson [45] gave a more precise meaning for the coherent structure as a “three-dimensional region of the flow over which at least one fundamental flow variable (velocity component, density, temperature, etc.) exhibits a significant correlation with itself or with another variable over a range of space and/or time that is significantly larger than the smallest local scales of the flow”. Under the Lagrangian frame of reference, the coherent structures can be identified by tracking the fluid parcel trajectories over a reasonable finite time. Based on this method, Haller [46,47] firstly used the finite-time Lyapunov exponent (FTLE) as a criterion for the identification of Lagrangian coherent structures (LCS). After the FTLE scalar field was established, the diagnostic and analysis of the vortex dynamics using LCS identification approaches became active, and they served in the studies of vortex ring pinch-off flow [48,49,50], the wake of pitching panels [51], circular cylinders [52], wings [53], vortex break downs [54], turbines draft tube vortexes [55], ocean current eddies [56] and hurricanes [57]. A comprehensive review and a rigorous assessment of the LCS diagnostic approaches and their application was conducted by Haller in 2014 [58].

Inspired by the success of LCS analysis in single-phase flows, some researchers have applied it to study two-phase cavitating flows. Tseng et al. [59] used LCS analysis to track the trajectories of the fluid particles in the cavitating flow around a hydrofoil. They discovered three different LCS groups on the hydrofoil surface and they enhanced the understanding of the interactions between the cavitation bubbles and the vortex region. Wang et al. [60] used LCS analysis to study the ventilated cavitation over a bluff body and established a correlation between the vortex shedding frequency and the vortex formation region. Hence, these works demonstrate that LCS analysis is a powerful tool to uncover the physics of the vortex structures in cavitating flows.

In this study, the Zwart homogeneous cavitation model [61] and the LES wall-adapting local eddy-viscosity turbulence model [62] are used to simulate the cavitating vortex shedding flow behind the truncated trailing edge of a 2D NACA 0009 hydrofoil. The numerical setup has been validated by comparison with existing experimental measurements. Based on the simulation results, it has been possible to find out the correlation between the degree of cavitation and the vortex shedding frequency and formation length. Eulerian and Lagrangian vortex identification criteria have been applied to the numerically predicted flow velocity field behind the trailing edge in order to analyze the evolution of the cavitating vortex shedding in comparison with the cavitation free conditions. As a result, the effects of the cavitation on the vortex shedding dynamic behavior have been identified and analyzed. The novelty of the present work lies in the fact that few attention has been given to investigate the cavitation effects on the Von Kármán vortex shedding dynamic behavior behind hydrofoils, which are widely used in hydraulic machinery. Moreover, it permits proving that the use of both Eulerian and Lagrangian vortex identification methods is a valuable tool to assess the dynamics and the morphology of the primary Von Kármán vortex shedding in cavitation conditions.

2. Methods

2.1. Experimental Setup Description

The experimental tests taken as a reference to validate our models were carried out at the EPFL high-speed cavitation tunnel inside a rectangular test section with the dimensions 150 × 150 × 750 mm with an NACA 0009 hydrofoil. The hydrofoil chord length, C, and the span, SP, were 100 and 150 mm, respectively. For additional details, see Ausoni’s thesis [20].

The hydrofoil was truncated at 90% of its original chord length and the resulting trailing edge thickness b was 3.22% of the new chord length, b = 3.22 mm. The hydrofoil was fixed at both sides and its span was equal to the width of the test section, 150 mm. Hydrofoil vibrations and flow velocity profiles were measured at different sections during the tests. The vortex-induced vibrations on the hydrofoil surface were measured with a laser vibrometer. The vibration signals permitted to determine hydrofoil natural frequencies and the vortex shedding frequency. The velocity profiles were measured using a laser Doppler velocimetry (LDV) technique. High-speed camera visualization techniques were used to show the cavitating flow characteristics. Further details of the experiments can be found in Ausoni [20].

2.2. Governing Equations

The numerical simulations were conducted with the commercial CFD solver code ANSYS CFX^® version 18.1. The turbulence model selected to close the unsteady Navier–Stokes equations was the LES wall-adapting local eddy-viscosity (WALE) model [62] because it can eliminate the additional false eddy viscosity of the Smagorinsky model [63] in the laminar flow regions and give a more precise prediction of the laminar–turbulence transition.

The homogeneous multiphase model was used which neglects the velocity differences between the liquid and vapor media. Thus, the corresponding governing equations are given by

$$\frac{\partial {\rho }_m}{\partial t}+\frac{\partial \left({\rho }_m{u}_j\right)}{\partial x_j}=0$$

(1)

$$\frac{\partial \left({\rho }_m{u}_i\right)}{\partial t}+\frac{\partial \left({\rho }_m{u}_i{u}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left({\mu }_m\frac{\partial u_i}{\partial x_j}\right)$$

(2)

where u_i is the velocity components and p is the pressure. The mixture density, ρ_m, and the mixture viscosity, µ_m, are expressed by

$${\rho }_m={\rho }_l{\alpha }_l+{\rho }_g{\alpha }_g$$

(3)

$${\mu }_m={\mu }_l{\alpha }_l+{\mu }_g{\alpha }_g$$

(4)

where ρ, µ and α are the density, viscosity and volume fraction, respectively, and subscripts l and g indicate liquid and gas, respectively.

Applying the implicit filtering operation on the unsteady Navier–Stokes equations, one can obtain

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho {\overline{u}}_j\right)}{\partial x_j}=0$$

(5)

$$\frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial t}+\frac{\partial \left(\rho {\overline{u}}_i{\overline{u}}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial {\overline{u}}_i}{\partial x_j}\right)-\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(6)

where ${\overline{u}}_j$ is the filtered velocity components and τ_ij is the sub-grid scale (SGS).

The Zwart–Gerber–Belamri cavitation mass transfer model [61] was used because of its better robustness and generalization. The corresponding transport equation of the vapor volume fraction, α_g, is given by

$$\frac{\partial {\rho }_g{\alpha }_g}{\partial t}+\frac{\partial \left({\rho }_g{\alpha }_g{u}_j\right)}{\partial x_j}=\dot{m}$$

(7)

$$\dot{m}=\{\begin{array}{c}{C}_p\frac{3{\alpha }_g{\rho }_v}{R}\sqrt{\frac{2}{3}\frac{\left(p-p_{\mathrm{v}}\right)}{{\rho }_l}},p>p_{\mathrm{v}}\\ -C_d\frac{3{\rho }_g\left(1-{\alpha }_g\right){\alpha }_{\mathrm{nuc}}}{R}\sqrt{\frac{2}{3}\frac{\left(p_{\mathrm{v}}-p\right)}{{\rho }_l}},p<p_{\mathrm{v}}\end{array}$$

(8)

where $\dot{m}$ is the mass transfer rate between the water liquid and vapor; and p_v is the saturated vapor pressure with a constant value of 2000 Pa. The constants in Equations (7) and (8) are given by the initial value of the bubble radius R = 1 µm, the nucleation site of the volume fraction α_nuc = 5 × 10⁻⁴ and the empirical condensation and vaporization coefficients C_d = 0.01 and C_p = 50.0, respectively.

2.3. Mesh and Numerical Setup

A 2D computational domain, as shown in Figure 1a, was considered because it has been assumed that the span-wise effects are not significant for the present study. The inlet surface was defined as semicircular and located at 1.2C upstream the hydrofoil leading edge. The outlet was located at 3.0C downstream the trailing edge.

The grid of the fluid domain close to the hydrofoil is shown in Figure 1b. In order to find a mesh-independent solution, four (4) different meshes with different refinement levels and dimensions were tested. As indicated in

Table 1. Details of the checked meshes and results obtained for the mesh independency study.

Name	Number of Elements	Max y+ [–]	Resolution in Near Wake [mm]	Outlet Boundary Location	Shedding Frequency [Hz]	St [–]	Deviation [%]
MC	29525	3.5	0.6	3C	1550	0.250	10.7
MM	49420	2.8	0.4	3C	1500	0.242	7.1
MF_1	71200	2.8	0.2	3C	1450	0.233	3.6
MF_2	313790	0.38	0.2	6C	1450	0.233	3.6

, the number of elements and the mesh resolution were increased progressively from mesh MC to mesh MF_1 while keeping the same domain dimensions. Another mesh, named MF_2, was also checked with the outlet boundary wall located at a distance of 6C from the trailing edge that was double the other three (3) meshes. The meshes MC, MM and MF_1 have a first layer height equal to 1 × 10⁻⁶ m, corresponding to a maximum y+ of 2.8 and an average y+ lower than 1. Meanwhile, mesh MF_2 has a first layer height equal to 1 × 10⁻⁷ m, a maximum y+ of 0.38 and an average y+ lower than 0.1.

A uniform inflow at 20 m/s was set at the inlet according to Ausoni [20]. Moreover, constant pressure was applied at the outlet boundary based on the corresponding cavitation number, σ$=\frac{p_{\infty }-p_{\mathrm{v}}}{1/2\rho V_{\mathrm{inlet}}^2}$, where $p_{\infty }$ is the reference pressure, which is set to equal the pressure located at the outlet boundary. The cavitation number range, σ/σ_i, has been taken from 1.3 to 0.4, where σ_i denotes the inception cavitation number.

2.4. Vortex Identification Methods

In order to identify the vortices being shed in the wake flow of the hydrofoil both at free cavitation and cavitation conditions, two criteria were applied based on the Eulerian and Lagrangian frames of reference, respectively. The first one is the Eulerian Q criterion which was proposed by Hunt et al. [41] and applied by Rockwood et al. [52], who found that there are no significant differences with other Eulerian criteria such as the λ₂ and Δ. So, the advantages of the Q criterion are its simplicity and its consistency with the rest of the Eulerian criteria. The second one is the finite-time Lyapunov exponent (FTLE) used by Haller [46,47]. The FTLE is a scalar derived from the fluid particle trajectories which measures the tendency of neighbor particles to separation over a finite time.

The Eulerian Q criterion is defined as the difference between the magnitude of the rotation rate, Ω, and of the strain tensor, S, given by

$$Q=\frac{1}{2}\left[{\parallel \mathsf{\Omega }\parallel }^2-{\parallel S\parallel }^2\right]$$

(9)

where Ω and S are derived from the velocity gradient tensor ∇u with $\mathsf{\Omega }=\frac{1}{2}\left(\nabla \mathrm{u}-\nabla {\mathrm{u}}^{*}\right)$ and $S=\frac{1}{2}\left(\nabla \mathrm{u}+\nabla {\mathrm{u}}^{*}\right)$. Q can be proved to be Galilean invariant. Where Q > 0, a vortex can be identified indicating that Ω is larger than S. However, in practice, only the condition Q > 0 cannot identify with enough precision the vortex region. Unfortunately, a threshold value is required for the vortex identification which is somehow subjective.

The Lagrangian FTLE criterion can be obtained from the velocity field data obtained from numerical or experimental results over the domain of interest. The fluid particle trajectories, which are also called the flow map x(t;x₀,t₀), are generated from the velocity field u(x,t) via differential equations by

$$\dot{\mathit{u}}=\left(\mathit{x},t\right) \mathit{x}\left(t_0;{\mathit{x}}_{\mathbf{0}},t_0\right)=x_{\mathbf{0}} \mathit{x}\in \mathcal{D} \mathrm{t}\in \left[t_0,t_1\right]$$

(10)

Then, considering an arbitrary point $\mathit{x}\in \mathcal{D}$ at the initial time t₀, when being advected by the flow, the new position at the time t₁ = t₀ + T is denoted by

$$\mathit{x}=F_{t_0}^{t_1}\left(x_{\mathbf{0}}\right)$$

(11)

Further, considering the trajectory of the virtual flow particle at y and t₀, denoted by y = x + δx(t₀), and assuming that δx(t₀)→ 0, this perturbation displacement after a time interval T is given by

$$\mathit{\delta }\mathit{x}\left(t_0+T\right)=F_{t_0}^{t_0+T}\left(\mathit{y}\right)-F_{t_0}^{t_0+T}\left(\mathit{x}\right)=\frac{d{F}_{t_0}^{t_0+T}\left(\mathit{x}\right) }{d\mathit{x}}\mathit{\delta }\mathit{x}\left(t_0\right)+\mathcal{O}\left({\parallel \mathit{\delta }\mathit{x}\left(t_0\right)\parallel }^2\right)$$

(12)

The symmetry matrix ∆, called the Cauchy–Green deformation tensor, is defined as

$$\Delta ={\left(\frac{d{F}_{t_0}^{t_0+T}\left(\mathit{x}\right) }{d\mathit{x}}\right)}^{*}\frac{d{F}_{t_0}^{t_0+T}\left(\mathit{x}\right) }{d\mathit{x}}$$

(13)

Then, the maximum stretching ratio, $\underset{\mathit{\delta }\mathit{x}\left(t_0\right)}{\max }\parallel \mathit{\delta }\mathit{x}\left(t_0+T\right)\parallel$, that occurs between points x and y is associated with the maximum eigenvalues of ∆:

$$\underset{\mathit{\delta }\mathit{x}\left(t_0\right)}{\max }\parallel \mathit{\delta }\mathit{x}\left(t_0+T\right)\parallel =e^{{\sigma }_{t_0}^T\left(\mathit{x}\right)\left|T\right|}\parallel \overline{\mathit{\delta }\mathit{x}\left(t_0\right)}\parallel$$

(14)

where

$${\sigma }_{t_0}^T\left(\mathit{x}\right)=\frac{1}{\left|T\right|}\ln \sqrt{{\mathsf{\lambda }}_{max}(\Delta )}$$

(15)

which is aligned with the eigenvector associated with λ_max(Δ), and ${\sigma }_{t_0}^T\left(\mathit{x}\right)$ represents the FTLE at $\mathit{x}\in \mathcal{D}$ and t₀ with a finite integration time in which T can take both positive and negative values. When T is positive, the FTLE scalar is called forward-time FTLE and it represents the degree of repulsion of the flow particles at t₀. When T is negative, the FTLE scalar is called backward-time FTLE and it represents the degree of attraction of the flow particles at t₀. From the ridges of the backward- and forward-time FTLEs, the vortex boundary of the complex flow at t₀ can be extracted. In fact, the computational cost to obtain the FTLE is much higher than to calculate the traditional Eulerian Q vortex criterion. However, FTLE does not require any threshold selection and the boundary of the vortex derived from the FTLE scalar field can objectively provide additional information of the vortex dynamics.

3. Results

3.1. Validation of the Numerical Simulation

In order to investigate the mesh independence, the lift and drag coefficients, C_L and C_D, of the truncated hydrofoil were calculated from the instantaneous lift and drag forces, F_y and F_x, acting on the hydrofoil with the following expressions:

$$C_L=\frac{F_y}{1/2\rho V_{\mathrm{inlet}}^2{A}_{ref}}$$

(16)

$$C_D=\frac{F_x}{1/2\rho V_{\mathrm{inlet}}^2{A}_{ref}}$$

(17)

where V_inlet is the free stream velocity with a value of 20 m/s, ρ is the water density with a value of 1000 kg/m³ and A_ref is the reference area calculated as the hydrofoil chord length times the span length.

The four (4) meshes previously mentioned were tested and the calculated shedding frequencies and the corresponding Strouhal numbers are indicated in Table 1. The Strouhal number, St, was calculated as

$$St=\frac{fb}{V_{inlet}}$$

(18)

where f is the vortex shedding frequency obtained with a fast Fourier transform (FFT) of the C_L time history. The shedding frequency experimentally obtained by Ausoni [20] is about 1400 Hz, corresponding to an St equal to 0.225. The percent deviations between the numerical estimates and the experimental value demonstrate that mesh MF_1 is a good compromise to model vortex shedding.

The ratio between the average mesh resolution Δ and Taylor scales for mesh MF_1 is plotted in Figure 2. The calculation of Δ and the Taylor scale is based on

$$\Delta =\sqrt{\Delta \mathrm{x}\Delta \mathrm{y}}$$

(19)

$$\mathrm{Taylor} \mathrm{scales}=C_{\mu }^{0.75}\frac{k^{1.5}}{\epsilon }$$

(20)

where Δx and Δy are the mesh sizes in x and y directions, respectively, C_µ is a constant value of 0.09, k is the turbulence kinetic energy and ε is the turbulence dissipation rate. The distributions of k and ε have been obtained from a steady-state simulation with mesh MF_1. Close to the hydrofoil surface and trailing edge area, the maximum ratio is less than 3, which suggests that the grid resolution has a similar order to the Taylor scales. Therefore, it is confirmed that mesh MF_1 with a resolution of 0.2 mm meets the adequate conditions required by the LES turbulence model.

The hydrofoil mean pressure coefficient, $C_p=\frac{p-p_{\infty }}{1/2\rho V_{\mathrm{inlet}}^2}$, was plotted along the chord as shown in Figure 3a where the dimensionless position x/L = 0 corresponds to the center of the hydrofoil chord line. The numerical results were compared with the results obtained with the vortex panel method by Bouziad [64] and they show a very good agreement around 90% of the chord length. Some discrepancy mainly lays on the trailing edge area which may result from the Kutta thickness trailing edge conditions involved in the vortex panel method. Figure 3b shows the good agreement between experimental measurements and the current numerical results of the mean normalized horizontal and vertical velocity profiles, u_x,mean/V_inlet and u_y,mean/V_inlet, along a vertical line located 2.0 mm downstream the trailing edge. The vertical locations were normalized by t_max which is the maximum thickness of the hydrofoil. For example, the maximum value of u_y,mean/V_inlet predicted by the the simulation is 0.26 and the experimental one is 0.24 with a percent deviation of around 8%.

The time histories of C_L and C_D are shown in Figure 4 for one free cavitation condition at σ/σ_i = 1.3, and for two cavitation conditions at σ/σ_i = 0.6 and σ/σ_i = 0.4. As it has been marked in Figure 3a, all the signals present an irregular evolution at the beginning of the run until the Von Kármán vortex shedding street behind the trailing edge is fully developed. After 0.01 s, a regular oscillation is achieved for both quantities in all the cases. Hence, this time was used as a reference starting time to obtain the mean values of the physical quantities of interest, e.g., velocity and pressure, over 10 cycles.

To validate the unsteady simulation of the vortex shedding dynamic behavior, FFT was applied to identify the main vortex shedding frequency. Figure 5 shows the spectra of C_L for all the simulated cases which permit identifying from the maximum frequency peak the corresponding shedding frequency. Under free cavitation conditions, the frequency is equal to 1450 Hz, but when cavitation occurs at σ/σ_i = 0.6, the frequency rises up to 1600 Hz and at σ/σ_i = 0.4, the frequency reaches a value of 1650 Hz.

Figure 6 permits comparing the current numerical simulation results with the experimental measurements in terms of St numbers. It can be confirmed that the numerical results are able to capture the linear increase in St as the cavitation grows in the wake flow with the reduction in the cavitation number.

3.2. Vortex Formation Length

The vortex formation length is a magnitude related to the vortex formation process at the near wake of the bluff body which is important for the overall evolution of the vortex shedding. There are several definitions of the vortex formation length and thorough notes about the definition methods are given by Griffin [65]. In the present work, the vortex formation length was found based on the maximum of the wake velocity fluctuation on the wake centerline, as proposed by Gerrard [5] for a bluff body wake at high Reynolds numbers. The same method has also been used in the works of Lee et al. [22] and Gnanaskandan et al. [13].

Figure 7 shows the root mean square value of the fluctuating velocity along the horizontal line behind the hydrofoil trailing edge normalized by the inlet flow velocity, $u_{\mathrm{rms}}^{\prime }/V_{\mathrm{inlet}}$, at free cavitation and cavitation conditions. When there is no cavitation for σ/σ_i = 1.3, the position of the peak value of $u_{\mathrm{rms}}^{\prime }/V_{\mathrm{inlet}}$ is located at 1.8 mm downstream of the trailing edge and this value is almost the same as the one obtained by Ausoni [20]. As the cavitation number further decreases, the position of the maxima of $u_{\mathrm{rms}}^{\prime }/V_{\mathrm{inlet}}$ increases gradually from 2.7 at σ/σ_i = 0.6 to 3.3 mm at σ/σ_i = 0.4. However, it is observed that the amplitude of $u_{\mathrm{rms}}^{\prime }/V_{\mathrm{inlet}}$ is generally reduced as the cavitation number decreases, which suggests that the velocity fluctuations decay due to the cavitation development.

3.3. Vortex Shedding Morphology and Dynamic Behavior

Figure 8 shows the time history of C_L over 6 × 10⁻⁴ s for the different cavitation conditions. Five typical instants of time were marked (from t₀ to t₄) to compare the vortex shedding process among the three flow conditions. For each instance, the instantaneous Eulerian Q criterion was derived from a snapshot of the velocity field and its gradient. Moreover, the contour of the backward FTLE scalar was also calculated considering as the initial time the one corresponding to each of the selected time instants. It should be mentioned here that the finite integration time, T, was chosen equal to −0.0015 s, which corresponds to almost two periods of the vortex shedding process. In fact, no significant differences are found among the results if larger integration times are considered.

In Figure 9, both the evolution of the Eulerian Q criterion on the instantaneous velocity field and the backward FTLE scalar levels display the alternating vortex shedding cycle behind the trailing edge at free cavitation conditions. At instant t₀, both the closer upper vortex to the trailing edge, marked with letter “A”, and the closer lower vortex, marked with letter “B”, are observed immediately behind the hydrofoil at the vortex formation region. At this time, vortex A has a larger size than vortex B. From t₀ to t₁, vortex B continues its growth and induces enough velocity to the wake flow so that the link between the upper boundary shear layer and vortex A is cut off. Thus, vortex A reaches its maximum size at instant t₁. Then, the upper vortex A is shed from the vortex formation area, which corresponds to the minimum value of the lift coefficient. From t₁ to t₄, a new upper vortex appears, marked with the letter “C”, and it grows as it is fed by the upper boundary shear layer. When it reaches sufficient size to cut off the previous lower vortex B, the process already described for vortex A is repeated by vortex C. In general, all the shed vortices take a circular shape and the two rows of vortices are clearly separated in the transverse direction. All these observations of the vortex shedding evolution are in agreement with the results for a bluff body obtained by Wang et al. [60].

As cavitation occurs, the instantaneous contour of the vapor volume fraction, the instantaneous Eulerian Q criterion and the contours of the backward FTLE scalars are plotted in Figure 10 and Figure 11 for two degrees of cavitation.

At σ/σ_i = 0.60, the pressure inside the vortex center decreases to the vapor pressure, p_v, and the vortices cavitate at their cores, as it can be seen in Figure 10a where the graphs show the evolution of the vapor volume fraction. Firstly, the upper cavitating vortex A is being formed at instant t₀. Then, vortex A is stretched and reaches its maximum size at instant t₁ when C_L is minimum. From t₁ to t₂, the vapor also fills a long filament that connects the cavitating upper vortex A with the trailing edge region. After that moment, vortex A is completely detached from the trailing edge and the lower cavitating vortex B is still attached to the trailing edge and begins to be stretched. Further, a new upper cavitating vortex C begins to be formed. From t₂ to t₃, vortex B continues to be stretched and its size keeps increasing. Finally and just before instant t₄, vortex B is cut off by vortex C and completely shed from the trailing edge area.

At σ/σ_i = 0.40, the instantaneous evolution of the vapor volume fraction presented in Figure 11a permits observing that the cavitating vortices are more stretched than at σ/σ_i = 0.60. From t₁ to t₂, it is observed that the upper cavitating vortex A is connected to the upper and lower cavitating vortices C and B, respectively, by a thin cavitation filament. Then, at t₃, the lower cavitating vortex B separates from the upper vortex C, and vortex A is completely shed downstream. At t₄, it can be seen how two vortices coexist in contact with the hydrofoil trailing edge which corresponds to vortex C and to the rear part of vortex B that occupies the very near wake and provokes the boundary layers coming from the extrados and the intrados to require a longer distance to create shed new vortices.

In Figure 10b or Figure 11b, the instantaneous contour plots of the Eulerian Q criterion show the evolution of the vortex shedding and the simultaneous velocity field under cavitation conditions. Compared with the free cavitation condition presented in Figure 9a, it is observed that the shedding vortex pattern and shapes of the vortices near the trailing edge are significantly affected by the occurrence of cavitation. For example, at instant t₁, the dimensions of vortex B tend to increase as the cavitation number decreases. Moreover, the vortices are elongated in the main flow direction and they cannot be entirely cut off by the opposite vortex. Consequently, their rear part remains attached to the blunt trailing edge instead of being shed with the front part.

In Figure 10c or Figure 11c, the instantaneous contour plots of the backward FTLE scalar show the evolution of the vortex shedding from the Lagrangian viewpoint. The obtained results are in agreement with the previous observation and thus confirm that the morphology and shape of the vortices continue to be affected by the increasing level of cavitation as well as their formation length.

In summary, at the free cavitation condition, the two rows of vortices are well separated in the transverse direction and each vortex depicts a well-defined circular shape. However, under cavitation conditions, the tranverse distance between the two rows is significantly reduced and the vortices change their shape towards an elliptical one mainly elongated in the transverse direction. As cavitation increases, the sizes of the vortices decrease and they are almost in line with the horizontal axis of the hydrofoil. Therefore, all the results seem to indicate that the effective height of the near wake formation region is significantly reduced which explains the fact that the shedding frequency is increased if the Strouhal number is assumed to be invariant.

4. Conclusions

The characteristics of the cavitating wake flow of a 2D NACA0009 with a blunt trailing edge have been numerically investigated and compared with the free cavitation conditions by means of the LES WALE turbulence model and the Zwart cavitation model. The numerical setup has been validated in comparison with experimental observations. To understand the effects of cavitation on the vortex shedding process, two vortex identification criteria based on the Eulerian and Lagrangian frames of reference have been applied to the unsteady flow velocity field. The following conclusions have been obtained:

• Vortices shape: it has been observed that the circular shape of the vortices changes to an elliptical one as they are stretched in the streamwise direction with the development of cavitation.
• Vortices transverse separation: it has been observed that the transverse separation between the upper and lower rows of vortices is reduced with the development of cavitation.
• Vortex formation length: it has been numerically found that the vortex formation length increases by around 83% with the development and growth of cavitation in the wake flow in comparison with the free cavitation condition when the cavitation number is reduced to 40% of the incipient cavitation number.
• Vortex shedding frequency: it has been numerically found that the vortex shedding frequency increases by 14% with the development and growth of cavitation in the wake flow in comparison with the free cavitation condition when the cavitation number is reduced to 40% of the incipient cavitation number.