Systematic Validation Study of an Unsteady Cavitating Flow over a Hydrofoil Using Conditional Averaging: LES and PIV

Abstract

We present results of Large-eddy simulations (LES) modeling of steady sheet and unsteady cloud cavitation on a two-dimensional hydrofoil which are validated against Particle image velocimetry (PIV) data. The study is performed for the angle of attack of 9° and high Reynolds numbers $R{e}_C$ of the order of ${10}^6$ providing a strong adverse pressure gradient along the surface. We employ the Schnerr–Sauer and Kunz cavitation models together with the adaptive mesh refinement in critical flow regions where intensive phase transitions occur. Comparison of the LES and visualization results confirms that the flow dynamics is adequately reproduced in the calculations. To correctly match averaged velocity distributions, we propose a new methodology based on conditional averaging of instantaneous velocity fields measured by PIV which only provides information on the liquid phase. This approach leads to an excellent overall agreement between the conditionally averaged fields of the mean velocity and turbulence intensity obtained experimentally and numerically. The benefits of second-order discretization schemes are highlighted as opposed to the lower-order TVD scheme.

1. Introduction

Cavitation often occurs in marine engineering applications and hydraulic machinery systems where the local pressure drops below the saturation vapor pressure of an operating liquid giving rise to a phase transition [1,2,3,4,5]. Under certain flow conditions, cavitation is strongly governed by hydrodynamic instabilities featuring unsteady flow dynamics. Enhanced noise, increased mechanical vibrations, and intensive material erosion are among harmful consequences of the cavitation development. In addition, these instabilities can trigger the resonance with separate elements or the whole duct of a hydraulic unit, leading to its irreversible damage or failure.

For careful modeling, prediction, and control of cavitating flows representing a complex multiphase, multiscale phenomenon, one benefits the most from simultaneous application of advanced experimental and numerical tools complementing each other. Modern experimental techniques used in cavitation research are typically nonintrusive and allow full-field measurements with a high spatial (and sometimes temporal) resolution. For example, high-speed imaging is usually utilized to track spatio-temporal evolution of cavitation structures and to determine their integral parameters such as cavity shape, dimensions, advection velocity, and shedding frequency [6,7,8]. Particle image velocimetry (PIV) or Particle tracking velocimetry (PTV) are applied to measure velocity distributions around a test body and evaluate statistical characteristics of the flow [9,10]. Visualization and PIV methods are often employed together and are capable of extracting accurate and reliable quantitative data.

The mathematical framework of numerical calculations typically addresses three issues concerning numerical discretization and simulations of cavitation and turbulence. The turbulence modeling is connected to selection of a certain type of averaging applied to the Navier–Stokes equations to reduce computational time. Without simplifications, we end up with Direct numerical simulations (DNS) when all spatio-temporal scales are resolved on a chosen computational mesh. In case of time (ensemble) averaging, the Reynolds-averaged Navier–Stokes equations (RANS) can be derived where all temporal scales are modeled leading in practice to low computational expenses [11]. A now mature framework of Large-eddy simulations (LES) appears to be a reasonable compromise between DNS and RANS relying on a low-frequency spatial filtering [12] and is well suited to describe unsteady flows with strong vortex shedding.

A suitable way to treat cavitation is to apply the homogeneous mixture model [13] where both liquid and vapor are considered as the same phase with variable density and viscosity fields. This approach is typically referred to as Volume of fluid (VOF) that requires an additional transport equation for the vapor volume fraction field. The source terms corresponding to evaporation and condensation processes are modeled by a number of ways [14,15,16,17], where the phase-change rate expressions are typically derived from a reduced form of the Rayleigh–Plesset equation for bubble dynamics [1]. Most simulations in the literature employ these models together with RANS using both steady and unsteady formulations [18,19,20,21,22,23,24,25,26,27,28,29,30]. Although a good agreement with experiments is typically obtained for some regimes, one has to tune the model reducing, for instance, the eddy viscosity in certain flow regions. A promising overall strategy is to couple the VOF and LES to describe unsteady cavitating flow [31,32]. Recent advances in high-performance computing led to a rapidly increasing number of studies focusing on turbulent flows over hydrofoils [33,34,35,36,37], bluff bodies [38,39], and other geometries [40,41,42,43,44,45,46] at high Reynolds numbers.

The choice of a modeling strategy and its further success are tightly connected to the coupling of a particular model and numerical accuracy of the time-dependent problem discretization. A cavitating flow features large density variations leading to additional numerical issues. The order of temporal and spatial discretization is not a key factor within the RANS framework, while LES typically requires at least second-order accurate schemes, such as central differences typically combined with upwind differences including TVD schemes [47,48,49,50]. Moreover, for both RANS and LES, one has to employ the Favre averaging [51] for variable density flows to keep the governing equations in a standard form, although this results in extra unclosed terms. Generally, this aspect is left unnoticed in literature when numerical and experimental results are compared.

In this paper, we focus on a cavitating flow over a 2D hydrofoil at a moderate angle of attack and high Reynolds numbers. In Section 2 we describe visualization and measurement techniques used to study spatio-temporal evolution of cavitation structures and to measure flow velocity. Further, in Section 3 details on numerical simulations for cavitating regimes of the flow are given and a number of aforementioned issues are addressed, such as the effects of the cavitation model and discretization scheme with a detailed comparison of the computations and experiments. For a more accurate matching of the numerical and experimental data, we propose an improved methodology, which is based on the conditional averaging of velocity fields taking into consideration only the liquid phase. Section 4 is the primary one where we present our numerical and experimental results and, based on the mean flow velocity and turbulence characteristics, comparatively analyze the accuracy and quantitatively assess the effectiveness of the suggested methodology of modeling. Finally, in Section 5 we summarize our findings and draw concluding remarks.

2. Experimental Setup and Flow Regimes

The experimental part of this research was performed in the cavitation tunnel in the Institute of Thermophysics SB RAS described previously [52]. The test section of the tunnel is a $1.3$ m long channel with a rectangular cross-section of $0.08\times 0.25$ m${}^2$. We selected a hydrofoil section with a 2D symmetric profile resembling a scaled-down model of guide vanes of a Francis turbine (to a scale of 1:13.26) as a test object, see [53]. The chord length C is 100 mm and the aspect ratio is $L/C=0.8$, where $L=80$ mm is the foil span that is equivalent to the width of the test section, see Figure 1. The centerline of the rotating axis of the profile coincides with the geometric center of the hydrofoil. The test model is supported as a cantilever on one of the test channel sidewalls adjoining tightly to both of them, with the distance from the inlet cross-section to the hydrofoil center equal to $7.5C=750$ mm.

We consider two flow regimes for the attack angle $\alpha =9°$ corresponding to transitional sheet and unsteady cloud cavitation regimes with the flow parameters provided in

Table 1. Flow regimes considered in present study.

	${\mathit{U}}_0$ [m/s]	${\mathit{p}}_0$ [kPa]	${\mathit{p}}_{\mathit{i}\mathit{n}}$ [kPa]	$\mathit{R}{\mathit{e}}_{\mathit{C}}$ [-]	$\mathit{\sigma }$ [-]
I	$10.47$	$105.8$	$85.3$	$1.32\times {10}^6$	$1.86$
II	$9.45$	$115.3$	$102.4$	$1.19\times {10}^6$	$2.49$

. When a cavitation sheet is quasisteady, its closure is almost stable and only minor oscillations of its length are detectable (regime II). In this regime, a vapor trail is often visible alongside the model surface and in the hydrofoil wake due to irregular detachments of tiny vapor clouds from the cavity closure [25,26,54,55]. Unsteady cloud cavitation reveals unstable behavior characterized by quasiperiodic large-scale vapor structures shedding (regime I). In fact, cloud cavitation is one of the most common forms of cavitating flow unsteadiness [5,56,57,58].

The Reynolds and cavitation numbers are defined as

$$R{e}_C=\frac{{\rho }_l{U}_{0}C}{{\mu }_l},\sigma =\frac{p_0-p_v}{{\rho }_l{U}_0^2/2},$$

(1)

where $p_0$ and $p_v$ are the reference static pressure at the test section sidewall gauged by a diaphragm strain-gauge pressure transducer near the inlet cross-section at a distance of $7C=700$ mm upstream from the leading edge of the hydrofoil and the saturation vapor pressure of the operating liquid (distilled water), ${\rho }_l$ and ${\mu }_l$ are the density and dynamic viscosity of the operating liquid. Subscripts l and v denote the liquid and vapor phase, respectively. The reference velocity $U_0$ corresponds to the incoming flow measured upstream of the hydrofoil by PIV in the central vertical longitudinal section over the entire height of the test channel close to the inlet. Both dimensionless quantities ($R{e}_C$ and $\sigma$) are controlled by changing $U_0$. At water temperature of 30 °C, ${\rho }_l$ = 995.61 kg/m${}^3$, ${\mu }_l=7.978\times {10}^4$ Pa × s and $p_v=4.24\times {10}^3$ Pa. Concentration of the dissolved oxygen in the operating liquid is monitored by a Mettler Toledo InPro 6850i sensor (measurement range from 6 ppb to saturation, relative error 1%) connected to a M400 multiparameter transmitter and maintained almost constant (about $7.4$ mg/L) for both flow regimes.

We used a high-speed Photron FASTCAM SA5 camera to visualize the overall cavitating flow and determine integral characteristics of cavitation structures and a PIV method to measure planar velocity distributions around the hydrofoil. Sufficient details on the measurement instrumentation employed for the PIV implementation, processing, and validation procedures to evaluate velocity vector fields and assessment of measurement uncertainties for various statistical characteristics are available in our recent study [59]. The processing routine applied to raw PIV images accounts for local tracer concentration. Velocity vectors are calculated in image areas where seeding particles are suspended corresponding merely to the liquid phase since the vapor phase is free of tracers. This is clearly visible in Figure 2, where velocity vectors are present only in the flow regions without cavitation structures while the regions occupied by the vapor phase do not contain velocity vectors. These measurement data features based only on the liquid phase require an appropriate conditional averaging applied to the numerical results for a consistent and systematic validation as described below.

3. Modeling and Computational Details

In this section, we describe the modeling strategy assuming isothermal flow conditions and neglecting the effects of the surface tension.

3.1. Governing Equations

The Navier–Stokes and VOF equations for a fluid with variable density and viscosity are as follows [51]:

$$\begin{array}{c}\hfill \frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\sigma }_{ij}}{\partial x_j},\end{array}$$

(2)

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_j\right)}{\partial x_j}=0,\end{array}$$

(3)

$$\begin{array}{c}\hfill \frac{\partial \left({\rho }_v\alpha \right)}{\partial t}+\frac{\partial \left({\rho }_v\alpha u_j\right)}{\partial x_j}=R,\end{array}$$

(4)

where

$$\begin{array}{c}\hfill {\sigma }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_k}{\partial x_k}\right)\end{array}$$

(5)

is the viscous stress tensor. The VOF framework determs $\rho$ and $\mu$ as follows:

$$\begin{array}{c}\hfill \rho =\alpha {\rho }_v+(1-\alpha ){\rho }_l,\mu =\alpha {\mu }_v+(1-\alpha ){\mu }_l.\end{array}$$

(6)

Thus, the scalar field $\alpha$ representing the vapor volume fraction varies from 0 for the liquid phase to 1 for the vapor.

3.2. Cavitation Modeling

The term R in Equation (4) is responsible for phase transitions. We employ the Schnerr–Sauer model [16] as a basic approach and compare results to the Kunz model [15]. The Schnerr–Sauer model is formulated in the following way:

$$\begin{array}{ccc}\hfill & & R=\frac{{\rho }_v{\rho }_l}{\rho }\alpha (1-\alpha )\frac{3}{R_b}\mathrm{sign}(p_v-p)\sqrt{\frac{3}{2}\frac{|p_v-p|}{{\rho }_l}},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill & & R_b={\left(\frac{3}{4\pi }\frac{1}{n_0}\frac{\alpha }{1-\alpha }\right)}^{1/3}.\hfill \end{array}$$

(8)

The only empirical parameter is $n_0$ corresponding to bubble concentration per a unit volume of liquid and determining the mean radius of a cavitation microbubble which is set to be $n_0=1.6\times {10}^{13}$ according to the original publication [16].

The Kunz model involves more empiricism:

$$\begin{array}{c}\hfill R=\frac{C_1\alpha {\rho }_v^2\min [p-p_v,0]}{{\rho }_l\left(0.5{\rho }_l{U}_0^2\right)t_0}+\frac{C_2(1-\alpha ){\alpha }^2{\rho }_v}{t_0{\rho }_l},\end{array}$$

(9)

where $t_0=C/U_0$ is the flow characteristic time, and $C_1=C_2=1.0\times {10}^3$ are empirical constants.

3.3. Large-Eddy Simulation Framework

In line with the common routine, we introduce Favre-averaged quantities convenient for variable-density flows [51,60], which are defined as follows:

$$\begin{array}{c}\hfill {\widehat{u}}_i=\widetilde{\rho u_i}/\tilde{\rho },\end{array}$$

(10)

where the tilde denotes local spatial filtering within the LES framework. Applying spatial filtering to Equations (2)–(4) after some manipulations, we arrive at the system of equations:

$$\begin{array}{c}\hfill \frac{\partial \left(\tilde{\rho }{\widehat{u}}_i\right)}{\partial t}+\frac{\partial \left(\tilde{\rho }{\widehat{u}}_i{\widehat{u}}_j\right)}{\partial x_j}=-\frac{\partial \tilde{p}}{\partial x_i}+\frac{\partial {\widehat{\sigma }}_{ij}}{\partial x_j}-\frac{\partial {\widehat{\tau }}_{ij}}{\partial x_j},\end{array}$$

(11)

$$\begin{array}{c}\hfill \frac{\partial \tilde{\rho }}{\partial t}+\frac{\partial \left(\tilde{\rho }{\widehat{u}}_j\right)}{\partial x_j}=0,\end{array}$$

(12)

$$\begin{array}{c}\hfill \frac{\partial \left({\rho }_v\tilde{\alpha }\right)}{\partial t}+\frac{\partial \left({\rho }_v\tilde{\alpha }{\widehat{u}}_j\right)}{\partial x_j}=\tilde{R}-\frac{\partial {\widehat{q}}_j}{\partial x_j},\end{array}$$

(13)

where

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{ij}=\tilde{\mu }\left(\frac{\partial {\widehat{u}}_i}{\partial x_j}+\frac{\partial {\widehat{u}}_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial {\widehat{u}}_k}{\partial x_k}\right).\end{array}$$

(14)

We express the continuity Equation (12) in the form:

$$\begin{array}{c}\hfill \frac{\partial {\widehat{u}}_j}{\partial x_j}=-\frac{1}{\tilde{\rho }}\frac{d\tilde{\rho }}{dt}=-\frac{{\rho }_v-{\rho }_l}{\tilde{\rho }}\frac{d\tilde{\alpha }}{dt},\end{array}$$

(15)

where the substantial derivative and Equation (6) for $\tilde{\rho }$ were used to obtain this last relation. The subgrid-scale stresses in Equation (11) contain the input from the convective and viscous terms due to the nonlinearity:

$$\begin{array}{c}\hfill {\widehat{\tau }}_{ij}=\left(\tilde{\rho }\widehat{u_i{u}_j}-\tilde{\rho }{\widehat{u}}_i{\widehat{u}}_j\right)+\left({\tilde{\sigma }}_{ij}-{\widehat{\sigma }}_{ij}\right).\end{array}$$

(16)

A similar term appears in Equation (13):

$$\begin{array}{c}\hfill {\widehat{q}}_j={\rho }_v\widetilde{\alpha u_j}-{\rho }_v\tilde{\alpha }{\widehat{u}}_j.\end{array}$$

(17)

Below we summarize a list of modeling assumptions:

• In a standard manner we ignore the commutation errors within the LES framework caused by the inhomogeneous mesh [12].
• The term ${\widehat{\tau }}_{ij}$ in Equation (11), see Equation (16), is modeled employing a subgrid-scale model without distinguishing between the effects of convective and viscous contributions.
• We add an additional term $\partial \left({\rho }_v\tilde{\alpha }(1-\tilde{\alpha }){\widehat{v}}_j\right)/\partial x_j$ to the left-hand-side of Equation (13), where $\widehat{\mathit{v}}={\widehat{\mathit{u}}}_v-{\widehat{\mathit{u}}}_l$ is the relative velocity vector between the two phases, which is also called the compression velocity [61]. It is sometimes described as an artificial term added to compress the interface, but in fact is directly derived from the mass conservation equation of the vapor phase. The compression velocity is modeled as:

  $$\min [C_{\alpha }|\widehat{\mathbf{u}}|,\max (|\widehat{\mathbf{u}}\left|\right)]\frac{\nabla \tilde{\alpha }}{|\nabla \tilde{\alpha }|},$$

  (18)

  which is based on the maximum velocity in some region of the interface, $C_{\alpha }=1$ is the constant that controls intensity of the interface compression.
• The term $\tilde{R}\left(\alpha \right)$ in Equation (13) is nonlinear; thus, it has to be modeled. We ignore this fact employing the relation $\tilde{R}(p,\rho ,\alpha )=R(\tilde{p},\tilde{\rho },\tilde{\alpha })$.
• The term ${\widehat{q}}_j$ in Equation (13), see Equation (17), is typically neglected to avoid additional diffusion in the equation for $\tilde{\alpha }$.

The subgrid stresses ${\widehat{\tau }}_{ij}$ are approximated using Dynamic k–equation subgrid-scale model (DKSGS) [62]. An additional transport equation is written as follows:

$$\frac{\partial \tilde{\rho }{k}_{sgs}}{\partial t}+\frac{\partial \tilde{\rho }{\widehat{u}}_j{k}_{sgs}}{\partial x_j}=-\rho {\widehat{\tau }}_{ij}{\widehat{S}}_{ij}-ϵ-\frac{\partial }{\partial x_j}\left(\tilde{\rho }(\nu +{\nu }_t)\frac{\partial k_{sgs}}{\partial x_j}\right),$$

where $k_{sgs}=\widehat{{\mathit{u}}^2}-{\widehat{\mathit{u}}}^2$ is the subgrid-scale (SGS) kinetic energy, ${\widehat{S}}_{ij}=(\partial {\widehat{u}}_i/\partial x_j+\partial {\widehat{u}}_j/\partial x_i)/2$ is the strain-rate tensor. The term ${\widehat{\tau }}_{ij}$ is expressed in the following way of the SGS eddy viscosity ${\nu }_t$:

$${\widehat{\tau }}_{ij}=-2{\nu }_t\widehat{S_{ij}}+\frac{2}{3}{\delta }_{ij}{k}_{sgs},$$

(19)

where ${\nu }_t=C_{\nu }\sqrt{k_{sgs}}\tilde{\Delta }$ and $C_{\nu }$ is an unknown coefficient. The dissipation term $ϵ$ is expressed according to the scaling arguments:

$$ϵ=C_{ϵ}\tilde{\rho }{k}_{sgs}^{3/2}/\tilde{\Delta },$$

(20)

where $\tilde{\Delta }$ is the characteristic width of the spatial filter and $C_{ϵ}$ is also an unknown coefficient. The Germano–Lilly [63,64] procedure is applied for the calculations to find the coefficients $C_{\nu }$ and $C_{ϵ}$.

3.4. Numerical Details

LES of the cavitating flow is performed using the computational code OpenFOAM [65] based on the finite volume method (FVM) and interPhaseChangeDyMFoam solver. Within the FVM approach, the computational domain is divided into polyhedral cells with the equations of motion written in the integral form. For the convective term discretization in the momentum equation, see Equation (11), we use two different schemes. The basic one is the second order accurate Linear upwind differencing scheme (LUDS) [47], see the results described in Section 4.4. The results obtained with LUDS are compared to simulations with MUSCL as well as the TVD approach (Total variation diminishing) with MINMOD limiter [49], see Section 4.4. The diffusion terms are discretized using second-order central differences. The van Leer scheme [66] is applied for the second term in Equation (13) to guarantee boundedness of the scalar field. For time discretization we use the Crank–Nicolson scheme [67]. The PISO scheme consisting of a predictor and corrector steps is used for velocity-pressure coupling [68].

The computational domain is shown in Figure 3, where x, y and z correspond to the streamwise, transversal, and spanwise directions, accordingly. An unstructured computational mesh contains around $10.1\times {10}^6$ nodes. In the simulation, the dynamic mesh refinement procedure [69] is employed leading to an increase in the number of cells to $20\times {10}^6$ nodes. Dynamic refinement is applied based on the local volume fraction of vapor $\alpha$ by dividing the hexagonal volumes into 8 subvolumes for the cells with $0.4<\tilde{\alpha }<0.6$, see the inset in Figure 4. For detailed information about mesh refinement, the reader is referred to Section 4.1. Nondimensional computational time step $\mathrm{\Delta }t=1.2\times {10}^{-5}$, which corresponds to the CFL number of $0.51$. On the walls we set the no-slip boundary condition. The inflow velocity profile is set in agreement with the experiments [59], while the convective boundary condition is put on the outflow. The inflow boundary conditions for the subgrid-scale kinetic energy $k_{sgs}$ and the volume vapor fraction $\tilde{\alpha }$ are set to fixed values close to zero. As for the pressure field, the Neumann conditions are specified at all boundaries, while the pressure at the inlet cross-section is extrapolated according to the pressure value $p_0$ measured upstream. To do that, the velocity field known from the experiment is integrated over the area of the inlet cross-section of the channel. Then, we find the change in pressure along the channel from the cross-section where $p_0$ is gauged to the inlet cross-section of the computational domain at a distance of $1.65C$ from the leading edge of the hydrofoil, Figure 3. The pressure $p_{in}$ obtained in this way is ultimately set at the inlet cross-section of the computational domain, see Table 1.

3.5. Conditional Averaging

A reasonable way to adequately compare LES and PIV data for unsteady flow conditions is to perform conditional averaging of velocity fields for certain values of $\tilde{\alpha }$ to account only for the liquid phase characteristics as in PIV, see Section 2. The conditional averaging of velocity fields is defined as:

$$\begin{array}{c}\hfill {\langle {\widehat{u}}_i\rangle }_{{\tilde{\alpha }}_{th}}=\langle {\widehat{u}}_i|\tilde{\alpha }<{\tilde{\alpha }}_{th}\rangle ,\end{array}$$

(21)

where averaging $\langle .\rangle$ is performed only for samples with $\tilde{\alpha }<{\tilde{\alpha }}_{th}$ in the same spatial position, where ${\tilde{\alpha }}_{th}$ is a fixed threshold. Obviously, comparison should be performed for real velocity fields rather than for Favre-averaged values. However, the influence of spatial filtering for both LES and PIV appears to be of minor significance as we can see from the direct comparison below.

The value of ${\tilde{\alpha }}_{th}$ was varied in a set of preliminary simulations showing a negligible influence on the results for small vapor fractions, thus, ${\tilde{\alpha }}_{th}=0.25$ is adapted. Figure 5 shows the effect of conditional averaging on the mean velocity for regime I. The described procedure makes it possible to obtain the velocity field very close to the experimental one.

4. Results and Discussion

We compare PIV and LES results for the two flow regimes, i.e., regime I ($R{e}_C=1.32\times {10}^6$) and regime II ($R{e}_C=1.19\times {10}^6$), see Table 1. The regime with the higher Reynolds number is considered first as it corresponds to unsteady cloud cavitation, which may be more suitable for the chosen LES model. Regime II for transitional sheet cavitation could be expected to be more challenging albeit delivering similarly good results as described below.

4.1. Mesh Refinement

The dynamic mesh refinement procedure was used to alter the mesh density throughout the domain during the calculation. We considered two cases of the mesh named ‘coarse’ and ‘fine’ to investigate the effect of the mesh refinement on capturing the cavitation structures. As mentioned above, a coarse mesh contained around $10.1\times {10}^6$ computational nodes. The number of nodes along the upper side of the hydrofoil was 380. In y-direction the number of computational nodes was 430 from the hydrofoil upper side to the upper wall of the domain and the average distance to the first node from the surface was equal to $2.3\times {10}^{-5}C$. The dynamic mesh refinement procedure increased a number of nodes along the upper side of the hydrofoil to 620 as well as in y direction to 550 leading to a decrease in the average distance to the first node almost twice. As for a typical large cavitation cloud with the transverse size of around $0.3C$, the number of computational nodes was ∼140 in the y-direction for the coarse mesh, while for the fine mesh it was ∼260. The number of cells increased nearly twice in the region where the cavitation was present. We also considered a mesh with an increased number of computational nodes in the region near the upper side of the hydrofoil. The refined region represented a volume of the size of ($\mathrm{\Delta }x/C\times \mathrm{\Delta }y/C\times \mathrm{\Delta }z/C$) = (1 × 0.395 × 0.8) built around the center of the foil. In this region, each computational volume was divided into 8 subvolumes leading to $26.5\times {10}^6$ nodes of computational mesh. The dynamic mesh refinement procedure led to a significant improvement in the accuracy of calculations, see Figure 6, demonstrating the refinement effect on the time-averaged profiles of the streamwise velocity component. Compared to the mesh with the dynamic refinement, the static refinement delivered similar results. Nevertheless, the dynamic refinement was chosen due to the large number of computational nodes for the static refinement.

Figure 7 shows typical moments of the cavity detachment obtained with and without the dynamic mesh refinement procedure and with the static mesh refinement over the upper side of the hydrofoil. Simulations without the refinement perform poorly in capturing the attached vapor clouds because the mesh size appears to be insufficient for a large number of fine structures to be resolved. The patterns of cavitation structures obtained using static refinement and dynamic refinement have no visual differences. However, due to the fact that the mesh with the static refinement has a larger number of computational nodes and consumes more computational resources, the mesh with the dynamic refinement is used in this paper.

4.2. Visualization and Flow Dynamics

Snapshots during high-speed recording were taken at a sampling rate of 20 kHz, i.e., 20 thousand images correspond to 1 second of real time lapse. This made it possible to directly compare experimental photographs of the cavity with isosurfaces of the instantaneous vapor fraction obtained in the LES calculations. In the simulation, the exposure time of each frame is 1/25,000 of a second. Figure 8 shows a set of photographs of the unsteady cloud cavitating flow together with an isosurface of the instantaneous vapor fraction provided by LES with the Schnerr–Sauer model for regime I. This regime is characterized by a progressive growth of a cavitation sheet followed by its quasiperiodic detachments from the hydrofoil surface, cloud cavity formation, and shedding. This flow unsteadiness is caused by the development of a re-entrant jet moving upstream underneath the attached cavity due to an adverse pressure gradient. According to the cloud size, we distinguish partial, see Figure 8c, and full detachment, see Figure 8f. Immediately after the cloud detachment, the cavity closure front obtains a U-shaped contour. In this case, the attached cavity grows again in the central part of the channel, while close to the sidewalls its length does not change substantially. This behavior is apparently pronounced both in the experiment and simulation, see a–f in Figure 8. The reduced frequency (Strouhal number) of cloud cavitation $St=fC/U_0$, where f is the cavity shedding frequency, is well reproduced by LES with $S{t}_{LES}\approx 0.20$ that is very close to $S{t}_{PIV}\approx 0.19$. The higher dominant frequency $S{t}_{LES}\approx 0.58$ in LES is slightly lower than that in the experiment with $S{t}_{Exp}\approx 0.61$. Such good agreement of the $St$ values for regime I confirms that the cavitation Schnerr–Sauer model is appropriate and the level of spatial resolution is sufficient for simulating unsteady cloud cavitation.

Contrary to regime I, regime II features quasisteady flow dynamics with uniform detachments of small-scale clouds from the sheet cavity closure in the spanwise direction. Similar to the previous case, Figure 9 shows two characteristic instants for quasi-steady cavitating flow conditions for which some low-amplitude variations in the length of the cavitation sheet are nevertheless observed in the calculations. According to the LES data, there exists a weak quasiperiodicity with the non-dimensional frequency $S{t}_{LES}\approx 0.23$, although visualization does not reveal any fluctuations of the cavity. This issue still remains unclear and will be addressed in detail in the future. Despite this small discrepancy, the Schnerr–Sauer model visually predicts regime II quite well.

4.3. Comparison of Cavitation Models

To analyze fields of the time-averaged velocity and its fluctuations, we use an improved methodology for a proper comparison of the PIV and LES data based on conditional averaging that accounts only for the liquid phase distinguished by the threshold ${\tilde{\alpha }}_{th}=0.25$, see Equation (21), and then match it against the time-averaged results corresponding to ${\tilde{\alpha }}_{th}=1$. The results are averaged over the total duration of the LES calculations $t_{LE{S}_{total}}\approx 66C/U_0$. The PIV measurements were performed with a high spatial resolution, which evidently affects the sampling rate and the total time of the experiment. Each experimental vector corresponds to a window of $8\times 8$ pixels or a real flow area of $0.47\times 0.47$ mm. PIV data were recorded with a frequency of 4 double frames per second. The high spatial resolution allowed us to make a detailed comparison of the velocity profiles in the selected cross-sections. Values of the mean velocity and its fluctuations were evaluated based on the statistics of 5000 realizations (instantaneous velocity fields), which corresponds to a 20-min series of consecutive measurements in real time or 1255.2$C/U_0$. Figure 10 compares the averaged streamwise velocity profiles for both regimes in the five cross-sections shown in Figure 2, namely $x=0.1C$, $0.3C$, $0.5C$, $0.7C$ and $0.9C$. The conditional averaging indeed significantly improves the LES results. Note that the Schnerr–Sauer model typically demonstrates results which are slightly closer to the PIV data in most of the positions compared to the Kunz model. The results for cross-sections 1 and 2 are in excellent agreement, but a small discrepancy is pronounced in positions 3, 4 for regime I and 3, 5 for regime II which is discussed below.

The visual data shown in Figure 8 and Figure 9 demonstrate that, while for regime I the cavitation is basically present in the middle part of the foil, for regime II vapor bubbles are highly concentrated at the trailing edge. This observation explains the discrepancy mentioned above between the PIV and LES results in positions 3, 4 for regime I, and 3, 5 for regime II in Figure 10 representing the most challenging flow regions for modeling as well as the large difference between the time- and conditionally averaged results in these areas. In cross-section 5 for regime I, all the profiles collapse due to intermittent appearances of cloud cavities which periodically pass through this region with a low frequency, leading to a negligible contribution of the vapor phase to the statistical characteristics. Position 5 for regime II demonstrates that the simulation is capable of reproducing a small recirculation zone near the hydrofoil trailing edge, while the time-averaged results point at a relatively high positive velocity due to the vapor phase.

As for the velocity fluctuations averaged conditionally and over the whole time interval, their comparison against the PIV data is shown in Figure 11. Again, it is visible that the conditional averaging significantly improves the LES calculations. The Schnerr–Sauer and Kunz models demonstrate similar results in all cross-sections showing good overall agreement for regime I, except for position 4 where the squared amplitude of the streamwise fluctuations is overpredicted by a factor of 2 compared to the experiments. For cross-sections 2 and 3 of regime I, where the cavitation is largely present, the use of the conditional averaging is essential to reduce the time-averaged value of the velocity fluctuations and match with the experiments. For regime II, the results in the foil midchord, see position 3 in Figure 11, show that LES overestimates the velocity fluctuations compared to PIV, while close to the trailing edge, see position 5, the curves converge for all cases.

As for the conditional averaged shear stress obtained by the Schnerr–Sauer and Kunz models, their profiles are presented in Figure 12 in comparison with the PIV data. For the cross-sections 2 and 3 the conditional averaging reduces the time-averaged value of the shear stress. In positions 4 and 5, the numerical and experimental results are close to each other. The values of the shear stresses for regime II are a bit overestimated but, closer to the trailing edge, they are in good agreement with the PIV measurements. The overall results of the PIV and LES comparison are excellent for most flow regions due to the use of the conditional averaging.

4.4. Comparison of Convective Schemes

Above we compared the two cavitation models using the second-order accurate LUDS scheme to discretize the convective term in the momentum equation, see Equation (11). In this subsection, LUDS is compared with a more dissipative MINMODE and less dissipative MUSCL schemes to highlight the importance of second-order schemes for multiphase LES. Figure 13 shows together profiles of the conditionally averaged velocity and its fluctuations simulated with the Schnerr–Sauer model and the ones measured in the experiment. In the first two cross-sections, slight differences between the curves for LUDS, MINMOD, and MUSCL are clearly seen, which leads to larger deviations downstream. Significant overestimation of the mean velocity in positions 4 and 5 in Figure 13 points out at an enhanced diffusion of the momentum in the spanwise direction for the MINMOD scheme. On the contrary, distributions of the velocity fluctuations exhibit larger deviations in the first three cross-sections, while in positions 4 and 5 they are qualitatively similar. Note that the most delicate near-wall characteristics are poorly predicted by MINMOD, indicating the importance of second-order schemes to describe cavitating flow satisfactorily. The comparison of LUDS and MUSCL showed that there is only a slight discrepancy between these models.

5. Conclusions

We presented results of Large-eddy simulations of transitional sheet and unsteady cloud cavitation flow regimes around a 2D hydrofoil at the angle of attack of 9° and high Reynolds numbers $R{e}_C=1.19\times {10}^6$ and $1.32\times {10}^6$ in comparison with the PIV data [53]. We employed finite-volume simulations with the adaptive mesh refinement in critical flow regions where intensive phase transitions occur to increase the overall accuracy together with the Schnerr–Sauer and Kunz cavitation models. The comparison of LES and visualization results confirmed that the flow dynamics in the unsteady cloud cavitation regime was correctly reproduced in the numerical simulations including the shedding frequencies of partial and full detachment of cloud cavities and typical U-shaped front of the cavitation sheet closure after the cavity breakup.

To compare averaged velocity fields, we take into account that the measured PIV data provided the information only for the liquid phase. Thus, we proposed a methodology based on the conditional averaging of instantaneous velocity fields including only relatively small vapor fractions. This approach led to excellent overall agreement of distributions of the conditionally averaged velocity over the leading edge of the hydrofoil for both flow conditions. A slight difference between the LES and PIV data for the unsteady cloud cavity regime was observed in the midchord region above the hydrofoil, where the cavitation was very intensive. For the transitional sheet cavitation regime, we observed some discrepancies near the trailing edge of the hydrofoil where the cavitation modeling is quite challenging. The profiles of velocity fluctuations and the main component of shear stresses averaged conditionally over the whole time interval showed that the proposed methodology significantly improves the LES results compared to the PIV data. Next, the importance of the second-order LUDS discretization schemes was highlighted in comparison with corresponding results for the lower-order TVD MINMOD scheme. The comparison of LUDS and MUSCL schemes showed a slight discrepancy. Nevertheless, the results obtained with the LUDS model are in a good agreement with PIV.

Further application of the verified models will allow to study dynamic characteristics of the cavitating flow in more detail with a particular attention to lower-Reynolds-number transitional regimes typically representing a challenge for the considered cavitation models.