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.

Keywords

Large eddy simulation nonlinear subgrid-scale model hydrofoil cavitation re-entrant jet

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 ρ}{\partial t} + \frac{\partial (ρ u_{i})}{\partial x_{i}} = 0

(1)

\frac{\partial ρ {\bar{u}}_{i}}{\partial t} + \frac{\partial}{\partial x_{j}} (ρ {\bar{u}}_{i} {\bar{u}}_{j}) = - \frac{\partial \bar{p}}{\partial x_{i}} + \frac{\partial}{\partial x_{j}} (μ \frac{\partial {\bar{u}}_{i}}{\partial x_{j}}) + \frac{\partial τ_{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, $μ$ is the laminar viscosity, and $ρ$ is the mixture density

μ = (1 - α_{l}) μ_{v} + α_{l} μ_{l}

(3)

ρ = (1 - α_{l}) ρ_{v} + α_{l} ρ_{l}

(4)

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

τ_{ij} = {\bar{u}}_{i} {\bar{u}}_{j} - \bar{u_{i} u_{j}}

(5)

DSM

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

τ_{ij}^{d} = - 2 C_{s} {\bar{Δ}}^{2} | \bar{S} | {\bar{S}}_{ij}

(6)

where ${\bar{S}}_{ij}$ is the rate-of-strain tensor, ${τ_{ij}}^{d}$ is the deviatoric part of the SGS stress tensor, $\bar{Δ}$ is the filter width, and $| \bar{S} | = {(2 {\bar{S}}_{ij} {\bar{S}}_{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{matrix} T_{ij} = ({\tilde{\bar{u}}}_{i} {\tilde{\bar{u}}}_{j} - \tilde{\bar{u_{i} u_{j}}}) \\ L_{ij} = T_{ij} - {\tilde{\bar{τ}}}_{ij} \end{matrix}

(7)

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

T_{ij}^{d} = - 2 C_{s} {\tilde{\bar{Δ}}}^{2} \tilde{| \bar{S} |} {\tilde{\bar{S}}}_{ij}

(8)

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

T_{ij}^{d} - {\tilde{\bar{τ}}}_{ij}^{d} = 2 C_{s} {\bar{Δ}}^{2} | \bar{S} | {\bar{S}}_{ij} - 2 C_{s} {\tilde{\bar{Δ}}}^{2} \tilde{| \bar{S} |} {\tilde{\bar{S}}}_{ij} = C_{s} M_{ij}

(9)

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

L_{ij}^{d} = C_{s} M_{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_{ij} L_{ij}}{M_{kl} M_{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

ν_{r} = max (ν_{r}, - ν)

(12)

which means that the backscatter cannot be larger than the viscous dissipation ( $ν$ 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{matrix} τ_{ij}^{d} = C_{s} {\bar{Δ}}^{2} | \bar{S} | {\bar{S}}_{ij} + C_{1} {\bar{Δ}}^{2} ({\bar{S}}_{ik} {\bar{Ω}}_{kj} - {\bar{Ω}}_{ik} {\bar{S}}_{kj}) + C_{2} {\bar{Δ}}^{2} \\ [\frac{1}{| \bar{S} |} ({\bar{Ω}}_{ik} {\bar{S}}_{kj}^{2} - {\bar{S}}_{ik}^{2} {\bar{Ω}}_{kj} + {\bar{Ω}}_{ik} {\bar{S}}_{km} {\bar{Ω}}_{mj} - \frac{1}{3} {\bar{Ω}}_{kl} {\bar{S}}_{lm} {\bar{Ω}}_{mk} δ_{ij} + Π_{s} {\bar{S}}_{ij}] \end{matrix}

(13)

where $Π_{s}$ is the second principal invariant of $S_{ij}$

Π_{s} = \frac{1}{2} (S_{kk} S_{mm} - S_{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 $({\bar{S}}_{ik} {\bar{S}}_{kj} - (1 / 3) {\bar{S}}_{mn} {\bar{S}}_{nm} δ_{ij})$ is introduced into the new model and the final form appears as

\begin{matrix} τ_{ij}^{d} = C_{s} {\bar{Δ}}^{2} | \bar{S} | {\bar{S}}_{ij} + C_{1} {\bar{Δ}}^{2} \\ [({\bar{S}}_{ik} {\bar{Ω}}_{kj} - {\bar{Ω}}_{ik} {\bar{S}}_{kj}) - ({\bar{S}}_{ik} {\bar{S}}_{kj} - \frac{1}{3} {\bar{S}}_{mn} {\bar{S}}_{nm} δ_{ij})] + C_{2} {\bar{Δ}}^{2} \\ [\frac{1}{| \bar{S} |} ({\bar{Ω}}_{ik} {\bar{S}}_{kj}^{2} - {\bar{S}}_{ik}^{2} {\bar{Ω}}_{kj} + {\bar{Ω}}_{ik} {\bar{S}}_{km} {\bar{Ω}}_{mj} - \frac{1}{3} {\bar{Ω}}_{kl} {\bar{S}}_{lm} {\bar{Ω}}_{mk} δ_{ij} + Π_{s} {\bar{S}}_{ij})] \end{matrix}

(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

[\begin{matrix} M_{ij} M_{ij} M_{ij} N_{ij} M_{ij} P_{ij} \\ N_{ij} M_{ij} N_{ij} N_{ij} N_{ij} P_{ij} \\ P_{ij} M_{ij} P_{ij} N_{ij} P_{ij} P_{ij} \end{matrix}] [\begin{matrix} C_{s} \\ C_{1} \\ C_{2} \end{matrix}] = [\begin{matrix} M_{ij} L_{ij}^{d} \\ N_{ij} L_{ij}^{d} \\ P_{ij} L_{ij}^{d} \end{matrix}]

(16)

where

M_{ij} = 4 {\bar{Δ}}^{2} \tilde{| \bar{S} |} {\tilde{\bar{S}}}_{ij} - {\bar{Δ}}^{2} \tilde{| \bar{S} | {\bar{S}}_{ij}}

(17)

\begin{matrix} N_{ij} = 4 {\bar{Δ}}^{2} [({\tilde{\bar{S}}}_{ik} {\tilde{\bar{Ω}}}_{kj} - {\tilde{\bar{Ω}}}_{ik} {\tilde{\bar{S}}}_{kj}) - ({\tilde{\bar{S}}}_{ik} {\tilde{\bar{S}}}_{kj} - \frac{1}{3} {\tilde{\bar{S}}}_{mn} {\tilde{\bar{S}}}_{nm} δ_{ij})] \\ - {\bar{Δ}}^{2} [(\tilde{{\bar{S}}_{ik} {\bar{Ω}}_{kj}} - \tilde{{\bar{Ω}}_{ik} {\bar{S}}_{kj}}) - (\tilde{{\bar{S}}_{ik} {\bar{S}}_{kj}} - \frac{1}{3} \tilde{{\bar{S}}_{mn} {\bar{S}}_{nm}} δ_{ij})] \end{matrix}

(18)

\begin{matrix} P_{ij} = 4 {\bar{Δ}}^{2} \\ [\frac{1}{\tilde{| \bar{S} |}} ({\tilde{\bar{Ω}}}_{ik} {\tilde{\bar{S}}}_{kj}^{2} - {\tilde{\bar{S}}}_{ik}^{2} {\tilde{\bar{Ω}}}_{kj} + {\tilde{\bar{Ω}}}_{ik} {\tilde{\bar{S}}}_{km} {\tilde{\bar{Ω}}}_{mj} - \frac{1}{3} {\tilde{\bar{Ω}}}_{kl} {\tilde{\bar{S}}}_{lm} {\tilde{\bar{Ω}}}_{mk} δ_{ij} - \frac{1}{2} {\tilde{\bar{S}}}_{kk}^{2} {\tilde{\bar{S}}}_{ij})] \\ - {\bar{Δ}}^{2} [\frac{1}{| \bar{S} |} (({\bar{Ω}}_{ik} {\bar{S}}_{kj}^{2} - {\bar{S}}_{ik}^{2} {\bar{Ω}}_{kj} + {\bar{Ω}}_{ik} {\bar{S}}_{km} {\bar{Ω}}_{mj} - \frac{1}{3} {\bar{Ω}}_{kl} {\bar{S}}_{lm} {\bar{Ω}}_{mk} δ_{ij} - \frac{1}{2} {\bar{S}}_{kk}^{2} {\bar{S}}_{ij})] \end{matrix}

(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 α_{l}}{\partial t} + \frac{\partial (α_{l} u_{i})}{\partial x_{i}} = \overset{\cdot}{m}

(20)

where $α_{l}$ is the volume fraction of liquid and $\overset{\cdot}{m}$ is the mass transfer rate between the two phases. To express the mass transfer rate, the Kunz et al. model⁹ is adopted

{\overset{\cdot}{m}}^{-} = \frac{C_{dest} ρ_{v} min [0, p - p_{v}] α_{l}}{ρ_{l} (0.5 ρ_{l} U_{\infty}^{2}) t_{\infty}}

(21)

{\overset{\cdot}{m}}^{+} = \frac{C_{prod} (1 - α_{l}) α_{l}^{2}}{ρ_{l} t_{\infty}}

(22)

where ${\overset{\cdot}{m}}^{-}$ and ${\overset{\cdot}{m}}^{+}$ represent condensation and evaporation, respectively; $C_{dest}$ and $C_{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 α_{l}}{\partial t} + \frac{\partial (α_{l} u_{i})}{\partial x_{i}} + \frac{\partial [U_{i} α_{l} (1 - α_{l})]}{\partial x_{i}} = 0

(23)

The compression velocity is calculated as

U = min [C_{l} | U |, max (| u |)] \frac{\nabla α_{l}}{| \nabla α_{l} |}

(24)

where $C_{l}$ is a constant which controls the compression. The surface tension is evaluated by

f_{σ} = σ κ \nabla α_{l}

(25)

where $σ$ is the surface tension coefficient, and the free surface curvature $κ$ is defined as

κ = - \nabla \cdot (\frac{\nabla α_{l}}{| \nabla α_{l} |})

(26)

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

\frac{\partial ρ {\bar{u}}_{i}}{\partial t} + \frac{\partial}{\partial x_{j}} (ρ {\bar{u}}_{i} {\bar{u}}_{j}) = - \frac{\partial \bar{p}}{\partial x_{i}} + \frac{\partial}{\partial x_{j}} (μ \frac{\partial {\bar{u}}_{i}}{\partial x_{j}}) + \frac{\partial τ_{ij}}{\partial x_{j}} + f_{σ}

(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 m$ and span width $0.19 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.3 c$ , 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} s$ to ensure that the maximum Courant number $Co < 1$ .

The pressure coefficient is defined as

C_{p} = \frac{p_{m} - p_{out}}{0.5 ρ 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.1 c - 0.9 c$ from the leading edge; $p_{out}$ is the pressure at outlet; and $U_{\infty} = 5.33 m / 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 $σ = (p_{out} - p_{sat}) / (0.5 ρ U_{\infty}^{2})$ ( $p_{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 $σ = 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 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 s$ . On the whole, the two models perform similarly: at $t = 0.32 - 0.64 s$ , the results agree well with the experimental data and the major pressure perturbation is reproduced around $t = 0.411 s$ , while at $t = 0.2 - 0.32 s$ , the peak value of the pressure moves approximately 0.0645 s (for numerical results, the peak value is at approximately $t = 0.31 s$ ; for experimental data, it is around $t = 0.245 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 s$ ) to the experimental data, and then the phase difference of two cycles is imposed on the previous cycle ( $t = 0 - 0.32 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 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 $Δ C = (C - mean (C)) / (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} = - τ_{ij}^{d} {\bar{S}}_{ij} / (U_{\infty}^{4} / ν)$ , 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.

Footnotes

Academic Editor: Moran Wang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors would like to acknowledge the financial support given by the Ministry of Education, Science and Technology Project (Grant No. 113010A), the Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20130008110047), and the National Natural Science Foundation of China (Grant No. 51209206).

References

Zhang

. A general approach for rectified mass diffusion of gas bubbles in liquids under acoustic excitation. J Heat Trans: T ASME 2014; 136: 042001.

Roohi

Pendar

Rahimi

. Simulation of three-dimensional cavitation behind a disk using various turbulence and mass transfer models. Appl Math Model. Epub ahead of print 2 July 2015. DOI: 10.1016/j.ijheatmasstransfer.2013.10.019.

Arndt

REA

. Some remarks on hydrofoil cavitation. J Hydrodyn 2012; 24: 305–314.

Pham

Larrarte

Fruman

. Investigation of unstable sheet cavitation and cloud cavitation mechanisms. J Fluid Eng: T ASME 1999; 121: 289–296.

Callenaere

Franc

Michel

. The cavitation instability induced by the development of a re-entrant jet. J Fluid Mech 2001; 444: 223–256.

Leroux

Astolfi

Billard

. An experimental study of unsteady partial cavitation. J Fluid Eng: T ASME 2004; 126: 94–101.

Kubota

Kato

Yamaguti

. A new modeling of cavitating flows: a numerical study of unsteady cavitation on a hydrofoil section. J Fluid Mech 1992; 240: 59–96.

Singhal

Athavale

. Mathematical basis and validation of the full cavitation model. J Fluid Eng: T SAME 2002; 124: 1–8.

Kunz

Boger

Stinebring

. A preconditioned Navier–Stokes method for two-phase flows with application to cavitation. Comput Fluids 2000; 29: 849–875.

10.

Roohi

Zahiri

Passandideh-Fard

. Numerical simulation of cavitation around a two-dimensional hydrofoil using VOF method and LES turbulence model. Appl Math Model 2013; 37: 6469–6488.

11.

Passandideh-Fard

Roohi

. Transient simulations of cavitating flows using a modified volume-of-fluid (VOF) technique. Int J Comput Fluid D 2008; 22: 97–114.

12.

Zhou

Wang

. Numerical simulation of cavitation around a hydrofoil and evaluation of a RNG k-ε model. J Fluid Eng: T ASME 2008; 130: 011302.

13.

Liu

. LES numerical simulation of cavitation bubble shedding on ALE 25 and ALE 15 hydrofoils. J Hydrodyn 2009; 21: 807–813.

14.

Luo

Arndt

REA

. Large Eddy Simulation and theoretical investigations of the transient cavitating vortical flow structure around a NACA66 hydrofoil. Int J Multiphas Flow 2015; 68: 121–134.

15.

Hidalgo

Luo

Escaler

. Numerical investigation of unsteady cavitation around a NACA 66 hydrofoil using OpenFOAM. In: 27th IAHR symposium on hydraulic machinery and systems, Montreal, QC, Canada, 22–26 September 2014. IOP Publishing.

16.

Tao

Katz

Meneveau

. Statistical geometry of subgrid-scale stresses determined from holographic particle image velocimetry measurements. J Fluid Mech 2002; 457: 35–78.

17.

Horiuti

. Roles of non-aligned eigenvectors of strain-rate and subgrid-scale stress tensors in turbulence generation. J Fluid Mech 2003; 491: 65–100.

18.

Germano

. Turbulence: the filtering approach. J Fluid Mech 1992; 238: 325–336.

19.

Lilly

. A proposed modification of the Germano subgrid-scale closure method. Phys Fluids A: Fluid 1992; 4: 633–635.

20.

Kosovic

. Subgrid-scale modelling for the large-eddy simulation of high-Reynolds-number boundary layers. J Fluid Mech 1997; 336: 151–182.

21.

Shih

Zhu

Liou

. Modeling of turbulent swirling flows. NASA technical memorandum 113112, 1997, http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19970034949.pdf

22.

Wang

Bergstrom

. A dynamic nonlinear subgrid-scale stress model. Phys Fluids A: Fluid 2005; 17: 035109-1–035109-15.

23.

OpenCFD. OpenFOAM^®, The opens source CFD toolbox, user guide. Berkshire: OpenCFD Ltd, 2012.

24.

Celik

Ghia

Roache

. Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. J Fluid Eng: T ASME 2008; 130: 078001.

25.

Yang

Cui

. Large eddy simulation of rotating turbulent channel flow with a new dynamic global-coefficient nonlinear subgrid stress model. J Turbul 2012; 13: 1–20.