Study of unsteady cavitation on NACA66 hydrofoil using dynamic cubic nonlinear subgrid-scale model

Abstract

In this article, we describe the use of a new dynamic cubic nonlinear model, a new nonlinear subgrid-scale model, for simulating the cavitating flow around an NACA66 series hydrofoil. For comparison, the dynamic Smagorinsky model is also used. It is found that the dynamic cubic nonlinear model can capture the turbulence spectrum, while the dynamic Smagorinsky model fails. Both models reproduce the cavity growth/destabilization cycle, but the results of the dynamic cubic nonlinear model are much smoother. The re-entrant jet is clearly captured by the models, and it is shown that the re-entrant jet cuts the cavity into two parts. In general, the dynamic cubic nonlinear model provides improvement over the dynamic Smagorinsky model for the calculation of cavitating flow.

Introduction

Cavitation is an unsteady flow phenomenon that usually has a negative or even destructive effect on hydraulic machinery embodied in pressure pulsation, erosion, and noise.^1,2 Cavitation occurs when the liquid pressure is lower than the saturated vapor pressure. In physical experiments, cavitation is produced by decreasing the pressure at the inlet or outlet of the flow field, leading to different types of cavitation that can be observed during this process. As the pressure decreases, an air cavity is stably attached to the surface of blades, followed by the development of cloud cavitation due to the destabilization of the cavity, resulting in the shedding of large vapor bubbles. Cloud cavitation is commonly associated with hydraulic machinery blades, and it has a highly dynamic flow pattern. ³ The intrinsic mechanism is often attributed to the re-entrant jet generated at the end of the cavity.^4–6

In order to simulate the unsteady characteristics of cavitation through numerical calculation, it is usually assumed that the vapor and liquid are mixed and that they share the same velocity and pressure with a no-slip velocity condition between the two phases, as proposed by Kubota et al. ⁷ Based on this assumption, Transport-Based Equation Modeling (TEM) is widely used in cavitating flows. With this method, an additional advection equation is solved. Many cavitation models have been used for this framework. For example, Singhal et al. ⁸ developed the full cavitation model; Kunz et al. ⁹ proposed a three-species formulation to account for the liquid, vapor, and the non-condensable gas. Apart from the advection equation, the other aspects related to TEM are the solution strategy for the advection equation and the turbulence model. ¹⁰

As a solution strategy, Volume of Fluid (VOF) method is a good technique with the advantage of the ability to capture the cavity interface accurately. A brief review of different VOF models can be found in Roohi et al. ¹⁰ The application of VOF in various cavitating flows has been shown to be successful in tracking the cloud and super-cavitation cavities.^10,11

For simulating turbulence, both Reynolds-averaged Navier–Stokes (RANS) and large eddy simulation (LES) models are widely used, although it seems that the results obtained with RANS modeling are not satisfactory without modification of the original formulations. For instance, Zhou and Wang ¹² studied a two-dimensional (2D) hydrofoil and found that the standard re-normalization group (RNG) k-ε model could not capture the cavity shedding unless the turbulent viscosity was modified to account for cavitation.

The theoretical basis of LES modeling is its ability to capture large vortical structures including cloud cavitation. Liu et al. ¹³ used the Smagorinsky–Lilly model to study the cavitation on hydrofoil and successfully predicted the existence of the re-entrant jet; Ji et al. ¹⁴ used the wall-adapting local eddy viscosity (WALE) model and also reproduced the cavity shedding. Both models needed no further modification, showing an obvious advantage over the regular RANS models. However, the simulated pressure pulsation, which is an important aspect of the cavitation process, tended to fluctuate more than the experiment.^14,15

A review of the literature shows that the LES models assume a linear relationship between the subgrid-scale (SGS) stress and the rate-of-strain, which may be the reason for the disagreement with experiment, as it has been figured out that the eigenvectors of the strain rate tensor and SGS stress tensor are not aligned but with a certain degree between each other in this linear relation.^16,17 A more reasonable relation to use would be the nonlinear tensorial polynomial constitutive form, which combines SGS stress with strain rate and vorticity.

In order to examine the performance of nonlinear models in simulating cavitating flow, a new nonlinear SGS model is proposed to simulate the cavitation around a hydrofoil. The results are compared with the experimental data obtained by Leroux et al. ⁶ and with the numerical results of the dynamic Smagorinsky model (DSM).

Governing equations

LES modeling

In the LES model, a low-pass filter (in frequency) is applied to the Navier–Stokes equations to get the following formulations (in Cartesian coordinates)

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$$

(1)

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

(2)

where the variables with a r bar over them represent the filtered quantities, u is the velocity, p is the pressure, $\mu$ is the laminar viscosity, and $\rho$ is the mixture density

$$\mu =(1-{\alpha }_l){\mu }_v+{\alpha }_l{\mu }_l$$

(3)

$$\rho =(1-{\alpha }_l){\rho }_v+{\alpha }_l{\rho }_l$$

(4)

where $\alpha$ is the volume fraction, and the subscripts v and l are the vapor and liquid, respectively. ${\tau }_{\mathrm{ij}}$ is the SGS stress, which is defined as

$${\tau }_{\mathrm{ij}}={\overline{u}}_i{\overline{u}}_j-\overline{u_i{u}_j}$$

(5)

DSM

DSM is a linear eddy viscosity model. The Smagorinsky model can be expressed as

$${\tau }_{\mathrm{ij}}^d=-2C_s{\overline{\Delta }}^2\left|\overline{S}\right|{\overline{S}}_{\mathrm{ij}}$$

(6)

where ${\overline{S}}_{\mathrm{ij}}$ is the rate-of-strain tensor, ${{\tau }_{\mathrm{ij}}}^d$ is the deviatoric part of the SGS stress tensor, $\overline{\Delta }$ is the filter width, and $\left|\overline{S}\right|={(2{\overline{S}}_{\mathrm{ij}}{\overline{S}}_{\mathrm{ij}})}^{1/2}$ is the characteristic filtered rate of strain.

In order to dynamically calculate the model coefficient $C_s$ , Germano ¹⁸ proposed to apply a second filter, called the test filter, to the SGS stress

$$\begin{array}{c}\multicolumn{1}{c}{T_{\mathrm{ij}}=({\tilde{\overline{u}}}_i{\tilde{\overline{u}}}_j-\widetilde{\overline{u_i{u}_j}})}\\ \multicolumn{1}{c}{L_{\mathrm{ij}}=T_{\mathrm{ij}}-{\tilde{\overline{\tau }}}_{\mathrm{ij}}}\end{array}$$

(7)

where “∼” is the test-filtering operation and $T_{\mathrm{ij}}$ is the stress at a test filter scale ( $\tilde{\overline{\Delta }}=2\overline{\Delta }$ ). Similar to equation (6)

$$T_{\mathrm{ij}}^d=-2C_s{\tilde{\overline{\Delta }}}^2\widetilde{\left|\overline{S}\right|}{\tilde{\overline{S}}}_{\mathrm{ij}}$$

(8)

where ${T_{\mathrm{ij}}}^d$ is the deviatoric part of $T_{\mathrm{ij}}$ . Combining equations (6) and (8)

$$T_{\mathrm{ij}}^d-{\tilde{\overline{\tau }}}_{\mathrm{ij}}^d=2C_s{\overline{\Delta }}^2\left|\overline{S}\right|{\overline{S}}_{\mathrm{ij}}-2C_s{\tilde{\overline{\Delta }}}^2\widetilde{\left|\overline{S}\right|}{\tilde{\overline{S}}}_{\mathrm{ij}}=C_s{M}_{\mathrm{ij}}$$

(9)

Using equation (9) to provide an approximation to $L_{\mathrm{ij}}^d$

$$L_{\mathrm{ij}}^d=C_s{M}_{\mathrm{ij}}$$

(10)

The use of least squares minimization of the error is suggested by Lilly ¹⁹ for the calculation of the value of $C_s$

$$C_s=\frac{M_{\mathrm{ij}}{L}_{\mathrm{ij}}}{M_{\mathrm{kl}}{M}_{\mathrm{kl}}}$$

(11)

It is well known that DSM suffers from numerical instability due to the large oscillation of the model coefficient. In this article, the eddy viscosity is bounded as

$${\nu }_r=max({\nu }_r,-\nu )$$

(12)

which means that the backscatter cannot be larger than the viscous dissipation ( $\nu$ is the kinetic viscosity).

Dynamic cubic nonlinear model

Kosovic ²⁰ proposed applying the concepts of nonlinear models of different levels of complexity for Reynolds stress closures to SGS model as they have been extensively studied by researchers and are much more developed. Therefore, Shih et al.’s method ²¹ is applied for the construction of the new model. The general form of the new model can be written as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\tau }_{\mathrm{ij}}^d=C_s{\overline{\Delta }}^2\left|\overline{S}\right|{\overline{S}}_{\mathrm{ij}}+C_1{\overline{\Delta }}^2({\overline{S}}_{\mathrm{ik}}{\overline{\Omega }}_{\mathrm{kj}}-{\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}})+C_2{\overline{\Delta }}^2}\\ \multicolumn{1}{c}{[\frac{1}{\left|\overline{S}\right|}({\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}^2-{\overline{S}}_{\mathrm{ik}}^2{\overline{\Omega }}_{\mathrm{kj}}+{\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{km}}{\overline{\Omega }}_{\mathrm{mj}}-\frac{1}{3}{\overline{\Omega }}_{\mathrm{kl}}{\overline{S}}_{\mathrm{lm}}{\overline{\Omega }}_{\mathrm{mk}}{\delta }_{\mathrm{ij}}+{\Pi }_s{\overline{S}}_{\mathrm{ij}}]}\end{array}$$

(13)

where ${\Pi }_s$ is the second principal invariant of $S_{\mathrm{ij}}$

$${\Pi }_s=\frac{1}{2}(S_{\mathrm{kk}}{S}_{\mathrm{mm}}-S_{\mathrm{kk}}^2)$$

(14)

Adapting the above model to the simple shear flow, it can be shown that equation (13) does not stand. Comparing with the Kosovic model, ²⁰ the term $({\overline{S}}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}-(1/3){\overline{S}}_{\mathrm{mn}}{\overline{S}}_{\mathrm{nm}}{\delta }_{\mathrm{ij}})$ is introduced into the new model and the final form appears as

$$\begin{array}{c}\multicolumn{1}{c}{{\tau }_{\mathrm{ij}}^d=C_s{\overline{\Delta }}^2\left|\overline{S}\right|{\overline{S}}_{\mathrm{ij}}+C_1{\overline{\Delta }}^2}\\ \multicolumn{1}{c}{[({\overline{S}}_{\mathrm{ik}}{\overline{\Omega }}_{\mathrm{kj}}-{\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}})-({\overline{S}}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}-\frac{1}{3}{\overline{S}}_{\mathrm{mn}}{\overline{S}}_{\mathrm{nm}}{\delta }_{\mathrm{ij}})]+C_2{\overline{\Delta }}^2}\\ \multicolumn{1}{c}{\left[\frac{1}{\left|\overline{S}\right|}({\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}^2-{\overline{S}}_{\mathrm{ik}}^2{\overline{\Omega }}_{\mathrm{kj}}+{\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{km}}{\overline{\Omega }}_{\mathrm{mj}}-\frac{1}{3}{\overline{\Omega }}_{\mathrm{kl}}{\overline{S}}_{\mathrm{lm}}{\overline{\Omega }}_{\mathrm{mk}}{\delta }_{\mathrm{ij}}+{\Pi }_s{\overline{S}}_{\mathrm{ij}})\right]}\end{array}$$

(15)

In order to dynamically calculate the three model coefficients (C_s, C ₁, and C ₂), the method proposed by Wang and Bergstrom ²² is applied. This is an extension of Germano’s method. The least square method is used

$$\left[\begin{array}{c}\hfill M_{\mathrm{ij}}{M}_{\mathrm{ij}}{M}_{\mathrm{ij}}{N}_{\mathrm{ij}}{M}_{\mathrm{ij}}{P}_{\mathrm{ij}}\hfill \\ \multicolumn{1}{c}{N_{\mathrm{ij}}{M}_{\mathrm{ij}}{N}_{\mathrm{ij}}{N}_{\mathrm{ij}}{N}_{\mathrm{ij}}{P}_{\mathrm{ij}}}\\ \multicolumn{1}{c}{P_{\mathrm{ij}}{M}_{\mathrm{ij}}{P}_{\mathrm{ij}}{N}_{\mathrm{ij}}{P}_{\mathrm{ij}}{P}_{\mathrm{ij}}}\end{array}\right]\left[\begin{array}{c}\hfill C_s\hfill \\ \multicolumn{1}{c}{C_1}\\ \multicolumn{1}{c}{C_2}\end{array}\right]=\left[\begin{array}{c}\hfill M_{\mathrm{ij}}{L}_{\mathrm{ij}}^d\hfill \\ \multicolumn{1}{c}{N_{\mathrm{ij}}{L}_{\mathrm{ij}}^d}\\ \multicolumn{1}{c}{P_{\mathrm{ij}}{L}_{\mathrm{ij}}^d}\end{array}\right]$$

(16)

where

$$M_{\mathrm{ij}}=4{\overline{\Delta }}^2\widetilde{\left|\overline{S}\right|}{\tilde{\overline{S}}}_{\mathrm{ij}}-{\overline{\Delta }}^2\widetilde{\left|\overline{S}\right|{\overline{S}}_{\mathrm{ij}}}$$

(17)

$$\begin{array}{c}\multicolumn{1}{c}{N_{\mathrm{ij}}=4{\overline{\Delta }}^2[({\tilde{\overline{S}}}_{\mathrm{ik}}{\tilde{\overline{\Omega }}}_{\mathrm{kj}}-{\tilde{\overline{\Omega }}}_{\mathrm{ik}}{\tilde{\overline{S}}}_{\mathrm{kj}})-({\tilde{\overline{S}}}_{\mathrm{ik}}{\tilde{\overline{S}}}_{\mathrm{kj}}-\frac{1}{3}{\tilde{\overline{S}}}_{\mathrm{mn}}{\tilde{\overline{S}}}_{\mathrm{nm}}{\delta }_{\mathrm{ij}})]}\\ \multicolumn{1}{c}{-{\overline{\Delta }}^2[(\widetilde{{\overline{S}}_{\mathrm{ik}}{\overline{\Omega }}_{\mathrm{kj}}}-\widetilde{{\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}})-(\widetilde{{\overline{S}}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}}-\frac{1}{3}\widetilde{{\overline{S}}_{\mathrm{mn}}{\overline{S}}_{\mathrm{nm}}}{\delta }_{\mathrm{ij}})]}\end{array}$$

(18)

$$\begin{array}{c}\multicolumn{1}{c}{P_{\mathrm{ij}}=4{\overline{\Delta }}^2}\\ \multicolumn{1}{c}{\left[\frac{1}{\widetilde{\left|\overline{S}\right|}}({\tilde{\overline{\Omega }}}_{\mathrm{ik}}{\tilde{\overline{S}}}_{\mathrm{kj}}^2-{\tilde{\overline{S}}}_{\mathrm{ik}}^2{\tilde{\overline{\Omega }}}_{\mathrm{kj}}+{\tilde{\overline{\Omega }}}_{\mathrm{ik}}{\tilde{\overline{S}}}_{\mathrm{km}}{\tilde{\overline{\Omega }}}_{\mathrm{mj}}-\frac{1}{3}{\tilde{\overline{\Omega }}}_{\mathrm{kl}}{\tilde{\overline{S}}}_{\mathrm{lm}}{\tilde{\overline{\Omega }}}_{\mathrm{mk}}{\delta }_{\mathrm{ij}}-\frac{1}{2}{\tilde{\overline{S}}}_{\mathrm{kk}}^2{\tilde{\overline{S}}}_{\mathrm{ij}})\right]}\\ \multicolumn{1}{c}{-{\overline{\Delta }}^2\left[\frac{1}{\left|\overline{S}\right|}(({\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{kj}}^2-{\overline{S}}_{\mathrm{ik}}^2{\overline{\Omega }}_{\mathrm{kj}}+{\overline{\Omega }}_{\mathrm{ik}}{\overline{S}}_{\mathrm{km}}{\overline{\Omega }}_{\mathrm{mj}}-\frac{1}{3}{\overline{\Omega }}_{\mathrm{kl}}{\overline{S}}_{\mathrm{lm}}{\overline{\Omega }}_{\mathrm{mk}}{\delta }_{\mathrm{ij}}-\frac{1}{2}{\overline{S}}_{\mathrm{kk}}^2{\overline{S}}_{\mathrm{ij}})\right]}\end{array}$$

(19)

We named the dynamic cubic nonlinear model (DCNM) as the new model contains a cubic term. In order to guarantee the numerical stability and save computation resources, the coefficients in equation (16) are expressed explicitly with the volume average applied to both the numerator and denominator of each coefficient’s expression.

Two-phase modeling

As discussed above, TEM is used and the mixture model is considered in this article. In the mixture model, the fluids are treated as a combination of vapor and liquid, and an additional transport equation is written as

$$\frac{\partial {\alpha }_l}{\partial t}+\frac{\partial \left({\alpha }_l{u}_i\right)}{\partial x_i}=\stackrel{·}{m}$$

(20)

where ${\alpha }_l$ is the volume fraction of liquid and $\stackrel{·}{m}$ is the mass transfer rate between the two phases. To express the mass transfer rate, the Kunz et al. model ⁹ is adopted

$${\stackrel{·}{m}}^{-}=\frac{C_{\mathrm{dest}}{\rho }_{v}min[0,p-p_v]{\alpha }_l}{{\rho }_l(0.5{\rho }_l{U}_{\infty }^2)t_{\infty }}$$

(21)

$${\stackrel{·}{m}}^{+}=\frac{C_{\mathrm{prod}}(1-{\alpha }_l){\alpha }_l^2}{{\rho }_l{t}_{\infty }}$$

(22)

where ${\stackrel{·}{m}}^{-}$ and ${\stackrel{·}{m}}^{+}$ represent condensation and evaporation, respectively; $C_{\mathrm{dest}}$ and $C_{\mathrm{prod}}$ are empirical constants; $U_{\infty }$ is the free-stream velocity; $t_{\infty }$ is the flow characteristic time defined as $t_{\infty }=L/U_{\infty }$ , where L is the body diameter; and $p_v$ is the saturation pressure.

In OpenFOAM, ²³ the VOF method is an improved version of CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes) as follows

$$\frac{\partial {\alpha }_l}{\partial t}+\frac{\partial \left({\alpha }_l{u}_i\right)}{\partial x_i}+\frac{\partial \left[U_i{\alpha }_l(1-{\alpha }_l)\right]}{\partial x_i}=0$$

(23)

The compression velocity is calculated as

$$\mathit{U}=min[C_l\left|\mathit{U}\right|,max\left(\left|\mathit{u}\right|\right)]\frac{\nabla {\alpha }_l}{|\nabla {\alpha }_l|}$$

(24)

where $C_l$ is a constant which controls the compression. The surface tension is evaluated by

$$f_{\sigma }=\sigma \kappa \nabla {\alpha }_l$$

(25)

where $\sigma$ is the surface tension coefficient, and the free surface curvature $\kappa$ is defined as

$$\kappa =-\nabla ·\left(\frac{\nabla {\alpha }_l}{|\nabla {\alpha }_l|}\right)$$

(26)

The surface tension is added to the momentum equation, and thus equation (2) is rewritten as

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

(27)

Numerical setup

The hydrofoil investigated in this study is an NACA66 series foil fixed in a square channel, with the chord length $c=0.15\mathrm{m}$ and span width $0.19\mathrm{m}$ . The specific parameters of the foil can be found in Leroux et al. ⁶ In order to reduce the calculation time, the span width is set to $0.3c$ , with the calculation domain as illustrated in Figure 1.

Figure 1.

Computational domain and boundary conditions.

The procedure proposed by Celik et al. ²⁴ is used to examine whether the grid is in the asymptotic range. To resolve the near-wall characteristics, $y^{+}$ is set to be located in the viscous sub-layer at walls. Three sets of grids with $8.8\times {10}^4$ , $2.2\times {10}^5$ , and $5.1\times {10}^5$ cells, respectively, are generated. DSM and the Kunz et al. model ⁹ are used to simulate the hydrofoil for the non-cavitating condition. The open-source code OpenFOAM ²³ is used for all the calculations. As stated above, a mixture method is employed for the two-phase modeling, and VOF is adopted to solve the additional transfer equation in the cavitation model. The PISO algorithm is used for the pressure–velocity coupling, and the time-step size is set to $1\times {10}^{-5}\mathrm{s}$ to ensure that the maximum Courant number $\mathrm{Co}<1$ .

The pressure coefficient is defined as

$$C_p=\frac{p_m-p_{\mathrm{out}}}{0.5\rho U_{\infty }^2}$$

(28)

where $p_m$ is the time-averaged ( $8\times {10}^4$ time-steps are used for averaging) pressure at the monitor points, which are equally spaced on the suction side of foil within $0.1c-0.9c$ from the leading edge; $p_{\mathrm{out}}$ is the pressure at outlet; and $U_{\infty }=5.33\mathrm{m}/\mathrm{s}$ is the free-stream velocity, which is set to be the inlet velocity. Both the fine-grid convergence index and the extrapolated value for $C_p$ are shown in Figure 2.

Figure 2.

(a) Fine-grid convergence index and (b) extrapolated value comparing to the result of fine grid.

From the error bar shown in Figure 2(a), it can be seen that the fine-grid convergence index of the grid is small with an acceptable maximum error below 0.07. Figure 2(b) shows that the extrapolated value is close to that of the fine grid (with the number of grid elements equal to $5.1\times {10}^5$ ). Hence, the fine grid is adopted for the modeling.

Results and discussions

According to the experimental results reported by Leroux et al., ⁶ when the cavitation number defined as $\sigma =(p_{\mathrm{out}}-p_{\mathrm{sat}})/(0.5\rho U_{\infty }^2)$ ( $p_{\mathrm{sat}}$ is the saturation pressure) reaches 1.25, the unsteady cavitating phenomenon will occur with the cavity periodically growing, shedding, and disappearing. This unsteady behavior will claim more requirements on the accuracy of the turbulence model. In the following calculations, DSM and DCNM are used under the cavitating condition of $\sigma =1.25$ .

Pressure coefficient and cavity shedding frequency

Figure 3 shows a plot of the time-averaged value of $C_p$ calculated along the suction side of the hydrofoil with the DSM and DCNM methods and the referenced experimental data. It can be seen that the pressure coefficients are overestimated at $x/c<0.5$ and underestimated at $x/c>0.5$ . This tendency is similar to the results obtained by Zhou and Wang, ¹² which can be clearly explained by the instantaneous pressure. When $x/c<0.5$ , the predicted pressure fluctuation shows a shorter duration of low pressure (as shown in Figure 6) which leads to the overestimation; when $x/c>0.5$ , the duration of high pressure is shorter (not shown), leading to the underestimation of the pressure coefficients. Despite the discrepancy with the experimental data, DCNM shows better agreement.

Figure 3.

Comparison of the pressure coefficients along the suction side.

Due to the intrinsic relationship between cavitation and pressure, it is convenient to determine the cavity shedding frequency from the pressure signal. As was done in Ji et al., ¹⁴ the power spectral density (PSD) is used in this article. Figure 4 shows the PSD profiles of two models, obtained in both cases by taking $7.3\times {10}^4$ samples at $x/c=0.6$ on the foil suction side. For both models, the cavity shedding frequency is 4.1 Hz, which is higher than the experimental value 3.625 Hz. In order to examine the models’ ability to capture the turbulence spectrum, the −5/3 slope is shown as an indication of the model accuracy. ¹⁴ It is clear that DCNM captures this slope accurately, while DSM fails, showing a much more flat tendency.

Figure 4.

PSD of DSM and DCNM at $x/c=0.6$ on the suction side: (a) DSM and (b) DCNM.

Cavity growth/destabilization cycle

The cavity growth/destabilization cycle can be clearly illustrated by the numerical results, as shown in Figure 5. The iso-surface of the vapor volume fraction of 0.1 is used to show the cavity shape. It is observed that the cavity basically experiences three stages: (1) the cavity occurs at the leading edge and develops downstream; (2) when the cavity length is maximum, the cavity is cut into a partial cavity and a cloud cavitation by a pressure perturbation (shown in Figure 6); and (3) when the cloud cavitation reaches the trailing edge, the cavitation suddenly disappears (not shown). Both the models correctly reproduce the whole cycle compared with the experimental results. It should be noted that three-dimensional characteristics are shown in the spanwise direction, which would probably lead to the overestimation of the shedding frequency due to the reduced span width used in the calculations.

Figure 5.

Cavity evolution in one cycle (time interval between two images: numerical data, 0.05 s; experimental data, 0.08 s): (a) DSM, (b) DCNM and (c) experiment.

Figure 6.

Comparison of instantaneous pressure at $x/c=0.4$ on the suction side: (a) DSM and (b) DCNM.

Three stages can be reflected in the temporal pressure signal on the suction side of the hydrofoil. Figure 6 shows a comparison of instantaneous pressure signals at $x/c=0.4$ on the suction side. In order to compare these values with the experimental data, it is necessary to choose a base point. Leroux et al. ⁶ found the major pressure perturbation occurs at approximately $t=0.411\mathrm{s}$ and cuts the cavity into two parts (stage 2 in Figure 5). Therefore, the points when the cavity is cut into two parts in DSM and DCNM are found accordingly and matched to $t=0.411\mathrm{s}$ . On the whole, the two models perform similarly: at $t=0.32-0.64\mathrm{s}$ , the results agree well with the experimental data and the major pressure perturbation is reproduced around $t=0.411\mathrm{s}$ , while at $t=0.2-0.32\mathrm{s}$ , the peak value of the pressure moves approximately 0.0645 s (for numerical results, the peak value is at approximately $t=0.31\mathrm{s}$ ; for experimental data, it is around $t=0.245\mathrm{s}$ ), which is equal to twice the phase difference between the numerical results and experimental data. This results from the step matching one cycle ( $t=0.32-0.64\mathrm{s}$ ) to the experimental data, and then the phase difference of two cycles is imposed on the previous cycle ( $t=0-0.32\mathrm{s}$ ). The major difference between the two results is that the data of DCNM tend to be much smoother than those of DSM. This may be the reason why DSM cannot reflect the right turbulence spectrum in Figure 4.

At stage 2, the cavity shedding is attributed to the re-entrant flow generated at the end of the cavity.^4–6 Figure 7 shows the velocity vector field around the foil at $t=0.411\mathrm{s}$ when the cavity is cut into two parts. In both images, the re-entrant flow is predicted at approximately $x/c=0.6$ , where the flow is cut by it. It should be noted that the re-entrant flow predicted by DSM is some distance away from the foil surface. Due to the pressure gradient, part of the main flow reverses direction and cuts the cavity from the bottom. ⁵ Therefore, the re-entrant flow should be exactly on the suction side as shown in Figure 7(b). This should be attributed to the insufficiency of DSM in the near-wall behavior, which is also documented in Yang et al. ²⁵

Figure 7.

Predicted re-entrant flow on the suction side: (a) DSM and (b) DCNM.

Model coefficients and SGS dissipation

From the above analysis, it is found that the instantaneous pressure of DCNM is smoother than the DSM’s. The underlying reason can be attributed to the model behavior during the calculation. The model coefficient is a representative quantity which will vary during the calculation. Figure 8 shows the instantaneous non-dimensional model coefficient, which is defined as $\Delta C=(C-\mathrm{mean}(C))/(\mathrm{mean}(C))$ at $x/c=0.4$ , where C represents the model coefficients of the models— $C_s$ in DSM and $C_s,C_1,C_2$ in DCNM.

Figure 8.

Comparison of the instantaneous non-dimensional model coefficient at $x/c=0.4$ .

As can be seen, the fluctuation of the model coefficient in DSM is much larger than in DCNM. This is a good indication that the pressure field calculated by DSM experiences more fluctuation.

Another index quantity is the non-dimensional SGS turbulent kinetic energy (TKE) production, which is defined as $P_r=-{\tau }_{\mathrm{ij}}^d{\overline{S}}_{\mathrm{ij}}/(U_{\infty }^4/\nu )$ , representing the TKE transferred to the residual motions. Figure 9 shows a comparison of the instantaneous SGS TKE production for the two models at $x/c=0.4$ . Most of the time, the SGS TKE production of DSM is negative, which means that the energy is transferred from the residual motions to the large-scale motions, namely, the backscatter. Therefore, the large-scale motions, including the pressure field fluctuation, tend to be influenced by the sub-grid scale motions. Additionally, the fluctuation of $P_r$ in DSM is also much larger than in DCNM, having an even larger effect on the large-scale motions.

Figure 9.

Comparison of the instantaneous non-dimensional SGS TKE production at $x/c=0.4$ .

Conclusion

In this article, we present the use of a linear model DSM and the new nonlinear model DCNM for the calculation of cavitating flow around an NACA66 series hydrofoil and the simulation results are compared with the experimental data. For both models, the cavity shedding frequency is overestimated, but the reason is still not clear. It is suspected that the reduced span width in the simulation may affect the shedding frequency as the cavitating flow shows three-dimensional characteristics. From the PSD analysis, it is found that DCNM can capture the turbulence spectrum while DSM fails to do so. The cavity growth/destabilization cycle, simulated by both models, duplicates the experimental data with little discrepancy.

Further investigation is done to explore the instantaneous pressure on the suction side. The results of DCNM are much smoother than in DSM, and this phenomenon may explain why DSM fails to capture the turbulence spectrum. Additionally, the re-entrant jet which cut the cavity into two parts is captured by the models. The instantaneous model coefficients and SGS TKE production are studied to reveal the underlying reason why DCNM predicts smoother pressure profile. It is found that both the model coefficient fluctuation and SGS TKE production are much larger for DSM. Additionally, significant backscatter of DSM is observed, which will affect the large-scale motions. Under these conditions, the pressure field is influenced in DSM, resulting in much fluctuation of the instantaneous pressure. In general, the nonlinear model does improve the ability to capture the turbulence, while it should be noted that the cavity shedding frequency is turbulence model independent. Further improvements of the calculation may focus on the cavitation model.