Comparative Performance Assessment between Incompressible and Compressible Solvers to Simulate a Cavitating Wake

1. Introduction

Vortical flow structures and cavitation are two common phenomena encountered in hydraulic machinery. As the fluid passes the fixed vanes, the guide vanes or the runner blades, it creates complex wake flow patterns that might develop in what is known as a vortex street. The resulting alternating shed vortices induce force fluctuations on the solids and provoke vortex-induced vibrations which may induce severe fatigue damage. In addition, vaporous cavities will be initiated and persist inside the center of the alternating shed vortices if the ambient pressure is low enough [1]. As these cavities move to the high-pressure region and collapse, they might also induce severe material erosion on adjacent surfaces. It is well known that the occurrence and development of cavitation can strongly alter the dynamic behavior of vortex street flows [2,3,4,5]. For example, it has been reported that cavitation tends to enhance the amplitude of the vortex-induced vibrations [4]. Therefore, further exploration of the effect of cavitation on the vortex street flow is essential to mitigate the damage caused by the cavitation and/or vortex-induced vibrations.

Currently, the computational fluid dynamic (CFD) method has been proved to provide detailed insights into complex flow fields, particularly in cases involving cavitation. Various numerical approaches, each with different levels of complexity, are available to simulate the cavitating flow. For a comprehensive review of the cavitation modeling approaches, see Niedźwiedzka et al. [6] and Folden and Aschmoneit [7]. Among these numerical models, the transport equation model (TEM) is one of the most widely used which includes an additional transport equation for the vapor phase volume fraction and a source term that accounts for mass transfer between phases. Without considering fluid compressibility, TEM-based simulations can effectively capture the transition from sheet to cloud cavitation, driven by a liquid re-entrant jet across different geometries [8,9,10], aligning well with experimental observations [11,12,13,14,15]. Conversely, the transition from the sheet cavitation to the cloud cavitation can be dominated by the condensation shock [16,17,18], which is highly related to the fluid compressibility of the mixture. Coupled with the equation of state for the liquid and vapor, TEM can successfully capture the condensation shock and its impact on the transition from sheet cavitation to cloud cavitation[10,18,19,20,21].

For the developed cavitation within the vortex street flows, bubble-vortex interactions alter the growth and collapse dynamics of cavitation [22,23,24,25,26]. This makes the roles for fluid compressibility in cavitating vortex street flows distinct from those observed in unsteady cloud cavitation. Previous numerical investigations have highlighted the impact of fluid compressibility on cavitating vortex street flows, particularly at extremely low cavitation numbers. For example, Brandao and Mahesh [27] have numerically captured the condensation shock on the surface of the circular cylinder where cavitation is highly developed. Additionally, Kim et al. [28] confirmed that fluid compressibility can alter the morphology of fully developed cavities when the cavitation number is less than half of the cavitation inception number. However, the effects of fluid compressibility on the dynamics of vortex street flows have received limited attention, especially at intermediate cavitation numbers.

The aim of this study is to demonstrate the influence of fluid compressibility on the dynamics of a cavitating vortex street by numerically resolving the cavitating flow field using both incompressible and compressible cavitation solvers. This paper is composed of five sections. Following this introduction, the second section explains the governing equations. The third section details the corresponding numerical treatment. In the fourth section, the implemented numerical solvers are validated, and the results from both the incompressible and compressible solvers are compared. Finally, the fifth section presents the conclusions and summarizes the results.

2. Governing Equations

2.1. Mass Conservation Equation

The equation for the conserved mass without external mass source is given by:

$$\partial \rho /\partial t+\nabla ·\left(\rho \mathit{u}\right)=0$$

(1)

here, $\rho$ is the density, $\mathit{u}$ is the velocity, and $t$ is the time. This equation is valid for both incompressible and compressible flows.

2.2. Momentum Conservation Equation

In an inertial reference frame, the conserved momentum equation can be written as follows:

$$\partial \left(\rho \mathit{u}\right)/\partial t+\nabla ·\left(\rho \mathit{u}\mathit{u}\right)=-\nabla p+\nabla ·\mathit{\tau }$$

(2)

where $p$ is the pressure and $\mathit{\tau }$ is the stress tensor which is defined as:

$$\mathit{\tau }=\mu \left[\nabla \mathit{u}+{\nabla \mathit{u}}^T-2/3\left(\nabla ·\mathit{u}\right)\mathit{I}\right]$$

(3)

here, $\mu$ is the fluid molecular viscosity and $\mathit{I}$ is the unit tensor.

2.3. Realizable $k$-Epsilon Delayed Detached Eddy Simulation (DDES) Model

In the realizable $k$-epsilon DDES model, the governing equation for the turbulence kinetic energy, $k$, can be written as follows:

$$\partial \left(\rho k\right)/\partial t+\nabla ·\left(\rho \mathit{u}k\right)=P_k-Y_k^{\mathrm{*}}+\nabla ·\left[\left(\mu +{\mu }_t/{\sigma }_k\right)\nabla k\right]$$

(4)

where $\epsilon$ is the dissipation, $P_k$ is the generation term of $k$, ${\mu }_t$ is the turbulent viscosity and ${\sigma }_k$ is the turbulent Prandtl number.Furthermore, the dissipation of turbulent kinetic energy, $Y_k^{\mathrm{*}}$, is given by:

$$Y_k^{\mathrm{*}}={\displaystyle \frac{\rho k^{\frac{3}{2}}}{l_{des}}}$$

(5)

where

$$l_{des}=min\left(l_{rke},l_{les}\right);l_{rke}={\displaystyle \frac{k^{\frac{3}{2}}}{\epsilon }};l_{les}=C_{des}{\mathsf{\Delta }}_{max}$$

(6)

here, ${\mathsf{\Delta }}_{max}$ is the maximum length of the cell edge and $C_{des}=0.61$ is the model constant. If $l_{rke}$= $l_{les}$, $l_{des}$ is replaced by the following functions:

$$l_{des}=l_{rke}-f_{y}max\left(0,l_{rke}-C_{des}{\mathsf{\Delta }}_{max}\right)$$

(7)

$$f_y=1-tanh\left[{\left(20.0r_y\right)}^3\right]$$

(8)

$$r_y={\displaystyle \frac{\left(\mu +{\mu }_t\right)/\rho }{\sqrt{\mathit{u}:\mathit{u}}\kappa y^2}}$$

(9)

where y represents the wall distance and κ = 0.41.

2.4. Cavitation Modelling

The volume fraction transport equation for the vapor volume fraction, ${\alpha }_v$, is given by:

$$\partial \left({\rho }_v{\alpha }_v\right)/\partial t+\nabla ·\left({\rho }_v{\alpha }_v\mathit{u}\right)=\dot{m}$$

(10)

The $\rho$ and $\mu$, of the fluid mixture are defined, respectively, as:

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

(11)

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

(12)

where ${\rho }_v$ and ${\rho }_l$ are the water and vapor densities, respectively, ${\mu }_v$ and ${\mu }_l$ are the water and vapor dynamic viscosities, respectively.

The mass transfer rate between the liquid and vapor due to the cavitation is modelled with a source term,$\dot{m}$, presented in equation (10). If the effects of viscosity, non-condensable gas, surface tension, and second-order derivation are neglected, the Rayleigh-Plesset equation can be simplified and written as:

$$\dot{R}=\sqrt{{\displaystyle \frac{2\left(P_{ref}-p_v\right)}{3{\rho }_l}}}$$

(13)

where $p$ is the saturated vapor pressure and $R$ is the bubble radius.

Assuming that the distribution of the nucleation inside the flow field is uniform and independent of the flow structures, then the relationship between the bubble radius, $R$, and ${\alpha }_v$ is given by:

$${\alpha }_v=n{\displaystyle \frac{4\pi }{3}}{R}^3$$

(14)

where $n$ is the number of nucleation per unit volume. Then, $\dot{m}$ is calculated by:

$$\dot{m}={\rho }_v{\dot{\alpha }}_v=n{\rho }_{v}3{\displaystyle \frac{4\pi }{3}}{R}^2\dot{R}$$

(15)

By submitting $n$ given by equations (13) and (14) into equation (15), $\dot{m}$ can be expressed as:

$$\dot{m}={\displaystyle \frac{3{\rho }_v{\alpha }_v}{R}}\sqrt{{\displaystyle \frac{2\left(P_{ref}-p_v\right)}{3{\rho }_l}}}$$

(16)

In the vaporization process, the nucleation site density must decrease as the vapor volume fraction increases. To model this process, $\alpha$ needs to be replaced with $\left(1-{\alpha }_v\right){\alpha }_{\mathrm{nuc}}$ in equation (16). Then, the final ZGB cavitation model [29] is presented as follows:

$$\dot{m}=\left\{\begin{array}{c}{F}_c{\displaystyle \frac{3{\alpha }_v{\rho }_v}{R_0}}\sqrt{{\displaystyle \frac{2}{3}}{\displaystyle \frac{\left(p-p_{\mathrm{v}}\right)}{{\rho }_l}}} p>p_{\mathrm{v}}\\ F_v{\displaystyle \frac{3{\rho }_v\left(1-{\alpha }_v\right){\alpha }_{\mathrm{nuc}}}{R_0}}\sqrt{{\displaystyle \frac{2}{3}}{\displaystyle \frac{\left(p_{\mathrm{v}}-p\right)}{{\rho }_l}}} p<p_{\mathrm{v}}\end{array}\right.$$

(17)

where $p_{\mathrm{v}}$ is the saturated vapor pressure. The initial value of the bubble radius is $R=1\mu \mathrm{m}$, the nucleation site of the volume fraction ${\alpha }_{\mathrm{nuc}}=5\times {10}^{-4}$. Here, the optimal empirical condensation and vaporization coefficients have been selected as $F_c=0.001$ and $F_v=50.0$, respectively.

2.5. Equation of State

The density of the compressible liquid phase can be computed using the Tait equation [20,30]:

$${\rho \left(p\right)}_l={\rho }_{l,sat}{\left({\displaystyle \frac{p+B}{p_v+B}}\right)}^{1/N}$$

(18)

and the corresponding sound speed is defined as

$${c\left(p\right)}_l=\sqrt{{\displaystyle \frac{\partial p}{\partial \rho }}}=\sqrt{{\displaystyle \frac{N}{{\rho }_{l,sat}}}\left(p_v+B\right){\left({\displaystyle \frac{p+B}{p_v+B}}\right)}^{\left(N-1\right)/N}}$$

(19)

where $B$ is the water bulk modulus, $3.1\times {10}^8$Pa, ${\rho }_{l,sat}$ is the liquid density at $p_v$, 998.18 kg/m3, and $N$ is a model constant equal to 7.15.

The density of the compressible vapor phase using the polytropic equation of state [31] can be written as:

$${\rho \left(p\right)}_v={\left(p/C_v\right)}^{1/n}$$

(20)

and the corresponding sound speed is:

$${c\left(p\right)}_v=\sqrt{{\displaystyle \frac{\partial p}{\partial \rho }}}=\sqrt{C_{v}n{\left(p/C_v\right)}^{\left(n-1\right)/n}}$$

(21)

where the model constant $C_v$ can be estimated from a density of 0.017 kg/m³ at $p_v$. The polytropic exponent $n$ for an adiabatic assumption is 1.4.

2.6. Sponge Layer Conditions

Up to date, although the use of the pressure-based method has achieved notable success in a wide range of compressible cavitating flows [32,33], several unresolved numerical difficulties remain, particularly concerning the determination of the inlet and outlet boundary conditions in CFD simulations. Typically, such calculations are conducted on a truncated domain of the entire system and the standard inlet/outlet boundary conditions may fail to allow some of the flow features to leave the computational domain as physically expected. For that reason, artificial treatments of the truncated domains are required to avoid pressure waves reflected at the boundaries. Therefore, the sponge layer has been adopted for its flexibility and simplicity [34]. The set of formulas of the sponge layer based on the pressure and velocity for the pressure-based method is given as follows:

$$\partial \rho /\partial t+\nabla ·\left(\rho \mathit{u}\right)=\mathit{\sigma }\left(\mathit{x}\right){\left.{\displaystyle \frac{\partial \rho }{\partial p}}\right|}_{p=p_{ref}}\left(p_{ref}-p\right)$$

(22)

$$\partial \left(\rho \mathit{u}\right)/\partial t+\nabla ·\left(\rho \mathit{u}\mathit{u}+p\mathit{I}-\mathit{\tau }\right)=\mathit{\sigma }\left(\mathit{x}\right)\rho \left({\mathit{u}}_{ref}-\mathit{u}\right)$$

(23)

where the subtitle “ref” refers to a target value of the flow variable and $\mathit{\sigma }\left(\mathit{x}\right)$ is the function of the damping coefficient, which is defined as:

$$\mathit{\sigma }\left(\mathit{x}\right)={\displaystyle \frac{3{{\alpha }_{pol}\left(\mathit{x}/L\right)}^{{\beta }_{pol}}}{2{∆t}^{{\gamma }_{pol}}}}$$

(24)

here, ${\alpha }_{pol}=1.0$, ${\beta }_{pol}=3.0$, ${\gamma }_{pol}=1.0$ and x is the distance from the leading edge of the sponge layer and $L$ is the length of the sponge layer, as shown in Figure 1.

3. Numerical Method

In this section, the governing equations are discretised using the Finite Volume Method (FVM) over the collocated unstructured mesh. More attention will be paid to the numerical algorithm for the pressure-based multiphase solver with/without consideration of fluid compressibility using the predefined macro “DEFINE_LINEARIZED_MASS_TRANSFER” inside the ANSYS Fluent^® platform.

3.1. Incompressible Mixture/VOF Model

Without consideration of fluid compressibility, the fluid densities are constant. The numerical algorithm for the mass conservation and volume fraction transport equation with linearized source terms is detailed in the following subsection.

3.1.1. Volume Continuity Equation

To guarantee mass conservation and numerical stability, the pressure-correction equation is based on the total volume continuity instead of the mass conservation equation.

$${\displaystyle \frac{1}{{\rho }_v}}\left[\partial \left({\rho }_v{\alpha }_v\right)/\partial t+\nabla ·\left({\rho }_v{\alpha }_v\mathit{u}\right)-\dot{m}\right]+{\displaystyle \frac{1}{{\rho }_l}}\left[\partial \left({\rho }_l{\alpha }_l\right)/\partial t+\nabla ·\left({\rho }_l{\alpha }_l\mathit{u}\right)+\dot{m}\right]=0$$

(25)

Integrating equation (25) over a control volume, the discretized form of the volume continuity is given by:

$$\left\{\begin{array}{c}{\displaystyle \frac{1}{{\rho }_v}}\left({\displaystyle \frac{{\rho }_v{\alpha }_v-{\rho }_v{{\alpha }_v}^0}{∆t}}dV+\sum _f{\rho }_v{\alpha }_v{\mathit{u}}_f{\mathit{A}}_f\right)+\\ {\displaystyle \frac{1}{{\rho }_v}}\left({\displaystyle \frac{{\rho }_l{\alpha }_l-{\rho }_l{{\alpha }_l}^0}{∆t}}dV+\sum _f{\rho }_l{\alpha }_l{\mathit{u}}_f{\mathit{A}}_f\right)\end{array}\right\}=\left\{\begin{array}{c}\left({\displaystyle \frac{1}{{\rho }_v}}-{\displaystyle \frac{1}{{\rho }_l}}\right)\left[{\dot{m}}^{\mathrm{*}}-{\left({\displaystyle \frac{\partial \dot{m}}{\partial p}}\right)}^{\mathrm{*}}\left(p^{\mathrm{*}}-p_{\mathrm{v}}\right)\right]+\\ \underset{ds/dp}{\underbrace{{\left({\displaystyle \frac{\partial \dot{m}}{\partial p}}\right)}^{\mathrm{*}}\left({\displaystyle \frac{1}{{\rho }_v}}-{\displaystyle \frac{1}{{\rho }_l}}\right)}}\left(p-p_{\mathrm{v}}\right)\end{array}\right\}$$

(26)

where $dV$ is the volume of the control cell, ${\mathit{A}}_f$ is the face area vector within the control volume, ${\rho }_v{\alpha }_v{\mathit{u}}_f$ and ${\rho }_l{\alpha }_l{\mathit{u}}_f$ are the vapor and liquid mass fluxes, respectively. On the collocated grid, the phase mass flux is computed using a Rhie–Chow interpolation.

More importantly, the term $ds/dp$ is the linearized mass transfer related to the pressure for the pressure correction equation, which can enhance the numerical stability of the cavitation simulation.

3.1.2. Second phase Fraction Equation

In ANSYS Fluent^®, only the transport equation for the second phase volume fraction is solved using the VOF or mixture model. Normally, the vapor phase is set to be the second phase. Thus, the vapor volume fraction is obtained with the vapor phase continuity equation as presented in equation (10). Correspondingly, the discretized form of the equation can be written as follows:

$${\displaystyle \frac{1}{{\rho }_v}}\left({\displaystyle \frac{{\rho }_v{\alpha }_v-{\rho }_v{{\alpha }_v}^0}{∆t}}dV+\sum _f{\rho }_v{\alpha }_v{\mathit{u}}_f{\mathit{A}}_f\right)={\displaystyle \frac{\dot{m}}{{\rho }_v}}dV$$

(27)

Assuming that the mass source term $\dot{m}$ can be rewritten as the function of ${\alpha }_v$ and ${\alpha }_l$:

$${\displaystyle \frac{\dot{m}}{{\rho }_v}}={\dot{m}}_0+S_l{\alpha }_l-S_v{\alpha }_v$$

(28)

With linearized the mass source term $\dot{m}$, the equation (27) can be rewritten as:

$${\displaystyle \frac{1}{{\rho }_v}}\left({\displaystyle \frac{{\rho }_v{\alpha }_v-{\rho }_v{{\alpha }_v}^0}{∆t}}dV+\sum _f{\rho }_v{\alpha }_v{\mathit{u}}_f{\mathit{A}}_f\right)=\left(S_c+S_p{\alpha }_v\right)dV$$

(29)

where

$$\left\{\begin{array}{c}{S}_c={\dot{m}}_0+S_l\\ S_p=-\left(S_v+S_l\right)\end{array}\right.$$

(30)

here, $S_p$ is the linear part of the source term $\dot{m}$, and $S_c$ is the part of $\dot{m}$ which cannot be linearized.

3.2. Compressible Mixture/VOF Model

In this section, the robust numerical algorithm for multiphase flows to handle fluid compressibility is introduced [32]. With the SIMPLE method, velocity and density fields can be decomposed into two parts as follows:

$$\left\{\begin{array}{c}{\rho }_v=\left({{\rho }_v}^{\mathit{*}}+{{\rho }_v}^{\mathit{\text{'}}}\right)={{\rho }_v}^{\mathit{*}}+\left({\displaystyle \frac{\partial {{\rho }_v}^{\mathit{*}}}{\partial p}}\right)p^{\text{'}}\\ {\rho }_l=\left({{\rho }_l}^{\mathit{*}}+{{\rho }_l}^{\mathit{\text{'}}}\right)={{\rho }_l}^{\mathit{*}}+\left({\displaystyle \frac{\partial {{\rho }_l}^{\mathit{*}}}{\partial p}}\right)p^{\text{'}}\\ u=\left({\mathit{u}}^{\mathit{*}}+{\mathit{u}}^{\mathit{\text{'}}}\right)\end{array}\right.$$

(31)

where “$\mathrm{*}$” indicates the tentative values and “$\mathrm{\text{'}}$” indicates the variable correction.

3.2.1. Volume Continuity Equation

Unlike the incompressible pressure-correction equation, the phase masses for vapor and liquid, ${\rho }_v{\alpha }_v$ and ${\rho }_l{\alpha }_l$, respectively, and the phase flux terms for vapor and liquid, ${\rho }_v{\alpha }_v\mathit{u}$ and ${\rho }_l{\alpha }_l\mathit{u}$, respectively,are computed with the following expressions:

$$\left\{\begin{array}{c} {\rho }_v{\alpha }_v={{\rho }_v}^{\mathbf{*}}{\alpha }_v+\left({\displaystyle \frac{\partial {{\rho }_v}^{\mathbf{*}}}{\partial p}}\right)p^{\mathrm{\text{'}}}{\alpha }_v\\ {\rho }_l{\alpha }_l={{\rho }_l}^{\mathbf{*}}{\alpha }_l+\left({\displaystyle \frac{\partial {{\rho }_l}^{\mathbf{*}}}{\partial p}}\right)p^{\mathrm{\text{'}}}{\alpha }_l\end{array}\right.$$

(32)

$$\left\{\begin{array}{c} {\rho }_v{\alpha }_{v}u=\left({{\rho }_v}^{\mathbf{*}}+{{\rho }_v}^{\mathbf{\text{'}}}\right){\alpha }_v\left({\mathit{u}}^{\mathbf{*}}+{\mathit{u}}^{\mathbf{\text{'}}}\right)=\left\{\begin{array}{c}{{\rho }_v}^{\mathbf{*}}{\mathit{u}}^{\mathbf{*}}{\alpha }_v+\left({\displaystyle \frac{\partial {{\rho }_v}^{\mathbf{*}}}{\partial p}}\right)p^{\mathrm{\text{'}}}{\mathit{u}}^{\mathbf{*}}{\alpha }_v+\\ {{\rho }_v}^{\mathrm{*}}{\mathit{u}}^{\mathbf{\text{'}}}{\alpha }_v+\underset{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}-\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}}{\underbrace{{{\rho }_v}^{\mathbf{\text{'}}}{\alpha }_v{\mathit{u}}^{\mathbf{\text{'}}}}}\end{array}\right\}\\ {\rho }_l{\alpha }_{l}u=\left({{\rho }_l}^{\mathbf{*}}+{{\rho }_l}^{\mathbf{\text{'}}}\right){\alpha }_l\left({\mathit{u}}^{\mathbf{*}}+{\mathit{u}}^{\mathbf{\text{'}}}\right)=\left\{\begin{array}{c}{{\rho }_l}^{\mathbf{*}}{\mathit{u}}^{\mathbf{*}}{\alpha }_l+\left({\displaystyle \frac{\partial {{\rho }_l}^{\mathbf{*}}}{\partial p}}\right)p^{\mathrm{\text{'}}}{\mathit{u}}^{\mathbf{*}}{\alpha }_l+\\ {{\rho }_l}^{\mathrm{*}}{\mathit{u}}^{\mathbf{\text{'}}}{\alpha }_l+\underset{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}-\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}}{\underbrace{{{\rho }_l}^{\mathbf{\text{'}}}{\alpha }_l{\mathit{u}}^{\mathbf{\text{'}}}}}\end{array}\right\}\end{array}\right.$$

(33)

Normally, the high-order correction terms in equation (33) are omitted due to their small magnitude compared with the rest of terms and their negligible influence on the resolved field. Then the discretised form of the volume continuity can be written as:

$$\left\{\begin{array}{c}{\displaystyle \frac{1}{{\rho }_v}}\left\{\begin{array}{c}{\displaystyle \frac{{{\rho }_v}^{\mathrm{*}}{{\alpha }_v}^{\mathrm{*}}+\left(\frac{\partial {{\rho }_v}^{\mathbf{*}}}{\partial p}\right)p^{\mathrm{\text{'}}}-{{\rho }_v}^0{{\alpha }_v}^0}{∆t}}dV+\\ \sum _f\left({{\rho }_{v,f}}^{\mathbf{*}}{\mathit{u}}^{\mathbf{*}}{\alpha }_v+{\left({\displaystyle \frac{\partial {{\rho }_v}^{\mathbf{*}}}{\partial p}}\right)}_f{p}^{\mathrm{\text{'}}}{\mathit{u}}^{\mathbf{*}}{\alpha }_v+{{\rho }_{v,f}}^{\mathrm{*}}{\mathit{u}}^{\mathbf{\text{'}}}{\alpha }_v\right){\mathit{A}}_f\end{array}\right\}+\\ {\displaystyle \frac{1}{{\rho }_l}}\left\{\begin{array}{c}{\displaystyle \frac{{{\rho }_l}^{\mathrm{*}}{{\alpha }_l}^{\mathrm{*}}+\left(\frac{\partial {{\rho }_l}^{\mathbf{*}}}{\partial p}\right)p^{\mathrm{\text{'}}}-{{\rho }_l}^0{{\alpha }_l}^0}{∆t}}dV+\\ \sum _f\left({{\rho }_{l,f}}^{\mathbf{*}}{\mathit{u}}^{\mathbf{*}}{\alpha }_l+{\left({\displaystyle \frac{\partial {{\rho }_l}^{\mathbf{*}}}{\partial p}}\right)}_f{p}^{\mathrm{\text{'}}}{\mathit{u}}^{\mathbf{*}}{\alpha }_l+{{\rho }_{l,f}}^{\mathrm{*}}{\mathit{u}}^{\mathbf{\text{'}}}{\alpha }_l\right){\mathit{A}}_f\end{array}\right\}\end{array}\right\}=\left\{\begin{array}{c}\left({\displaystyle \frac{1}{{\rho }_v}}-{\displaystyle \frac{1}{{\rho }_l}}\right)\left[{\dot{m}}^{\mathrm{*}}-{\left({\displaystyle \frac{\partial \dot{m}}{\partial p}}\right)}^{\mathrm{*}}\left(p^{\mathrm{*}}-p_{\mathrm{v}}\right)\right]+\\ +\underset{ds/dp}{\underbrace{{\left({\displaystyle \frac{\partial \dot{m}}{\partial p}}\right)}^{\mathrm{*}}\left({\displaystyle \frac{1}{{\rho }_v}}-{\displaystyle \frac{1}{{\rho }_l}}\right)}}\left(p-p_{\mathrm{v}}\right)\end{array}\right\}$$

(34)

where ${{\rho }_{v,f}}^{\mathit{*}}{\mathit{u}}^{\mathit{*}}{\alpha }_v$ is the phase mass flux at the cell face.

The above equation not only couples the pressure and velocity but also involves the interaction between the flow field and the mass transfer source term. Note that the terms ${\displaystyle \frac{\partial {{\rho }_v}^{\mathrm{*}}}{\partial p}}$ and ${\displaystyle \frac{\partial {{\rho }_l}^{\mathrm{*}}}{\partial p}}$ are related to the effect of the vapor and liquid compressibility on the pressure correction, respectively. By default, these terms are approximated by the upwind scheme. The compressibility-related term in the transient term and mass transfer source term can be easily linearized, which can significantly enhance the numerical stability without affecting the final solution.

3.2.2. Second Phase Fraction Equation

With consideration of the compressible effect of the second phase, the discretised form of the volume fraction equation can be written as follows when the second phase is the vapor:

$${\displaystyle \frac{1}{{\rho }_v}}\left({\displaystyle \frac{{\rho }_v{\alpha }_v-{{\rho }_v}^0{{\alpha }_v}^0}{∆t}}dV+\sum _f{\rho }_{v,f}{\alpha }_v{\mathit{u}}_f{\mathit{A}}_f\right)=\left(S_c+S_p{\alpha }_v\right)dV$$

(36)

And the left-hand side of the equation can be rewritten as:

$${\displaystyle \frac{{\alpha }_v-{{\alpha }_v}^0}{∆t}}dV+\sum _f{\alpha }_v{\mathit{u}}_f{\mathit{A}}_f=\left\{\begin{array}{c}\left(S_c+S_p{\alpha }_v\right)dV+\\ \left\{\underset{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}}{\underbrace{-\left({\displaystyle \frac{{\rho }_v-{{\rho }_v}^0}{{\rho }_v}}{\displaystyle \frac{{{\alpha }_v}^{0}dV}{∆t}}\right)+\left(\sum _f{\displaystyle \frac{{\rho }_v-{\rho }_{v,f}}{{\rho }_v}}{\alpha }_v{\mathit{u}}_f{\mathit{A}}_f\right)}}\right\}\end{array}\right\}$$

(37)

where the last term is the expansion source term accounting for the effect of fluid compressibility on the volume conservation equation. It is composed of two parts: one corresponding to the transient term and the other one to the advection term.

3.3. Pressure Limits

As we know, the pressure-based solver can provoke an unbounded pressure field, e.g. negative pressure field. For the compressible flow, it is necessary to limit the pressure since the obtained negative pressure field is out of the predefined range for the fluid material. However, limiting the pressure is a challenging numerical issue, which will lead to a decrease in the convergence rate and to numerical instability. In the case of the multiphase compressible flows, the situation will be more severe. As mentioned by Li and Vasquez [32], it is difficult to directly limit the pressure above zero in the region where only a negligible amount of gas is present. In such a case, the limitation of the pressure must consider the effect of the local flow characteristic. According to the pressure-limited method proposed by Li and Vasquez [32], the fluid density is computed using the barotropic law with a pre-described pressure limit when the local pressure turns negative. For instance, the Tait equation for the compressible liquid density with the pressure-limited can be rewritten as:

$${\rho \left(p\right)}_l={\rho }_{l,sat}{\left[{\displaystyle \frac{max\left(p,p_{lim}\right)+B}{p_v+B}}\right]}^{1/N}$$

(38)

and the corresponding sound speed with the pressure-limited is defined as

$${c\left(p\right)}_l=\sqrt{{\displaystyle \frac{\partial p}{\partial \rho }}}=\sqrt{{\displaystyle \frac{N}{{\rho }_{l,sat}}}\left(p_v+B\right){\left[{\displaystyle \frac{max\left(p,p_{lim}\right)+B}{p_v+B}}\right]}^{\left(N-1\right)/N}}$$

(39)

where $p_{lim}$ is the pre-described pressure limit.

4. Validation

4.1. Case1: 1-D Two-Phase Time-Dependent Test Case

In this section, the compressible cavitation model is used to describe the effects of a density reduction below the density value of the saturated liquid induced by two symmetrical expansion waves as in the numerical setup defined by Schmidt [35]. The domain consists of a 1-D tube with a length of 1 m. At time $t$= 0 s, the tube is full of pure liquid water, and the pressure is constant with a value of $p$= 0.9 bar. The velocity field is assumed to jump at $x$= 0.5 m from the left velocity, $u_L$= - 10 m/s, to the right velocity, $u_R$= 10 m/s, to force the phase change from liquid to vapor to occur. Furthermore, the domain has been divided into 1000 cells and the time integration has been performed using the first-order implicit scheme with a time step $∆t$ = 1∙10^-7 s.

To verify the implementation of the compressible cavitation model, the current numerical results have been compared with the ones obtained by Schmidt [35]. Figure 2 presents the pressure, the velocity, the vapor volume fraction, and the sound speed distributions along the tube at time $t$= 1.5∙10^-4 s using the present compressible cavitation model and the results obtained by Schmidt [35]. Without experiencing numerical oscillations in regions with high gradients, the current implemented solver is proved to be numerically stable. Furthermore, the distribution of the flow quantities along the tube is almost identical between the current ones and the those extracted from the reference simulation. The precise agreement between both sets of results demonstrates the validity of the implemented compressible cavitation model.

4.2. Case2: Cavitating Flow over a Circular Cylinder

The present simulations of the non-cavitating and cavitating flows over a stationary cylinder have been performed with a large 2D computational domain of dimensionless dimensions relative to the cylinder diameter, D, in the horizontal and vertical ranges -50 ≤ $x/D$ ≤ 50 and -50 ≤$y/D$≤ 50, respectively, to avoid the effect of the domain boundaries, as depicted in Figure 3. The cylinder has been located in position (0, 0) and the number of grid elements in the circumferential and radial directions are 480 and 220, respectively. The symbol, $\theta$, represents the angle around the cylinder surface measured from its front stagnation point. The inflow condition of $u=u_{ref}$ and the pressure condition of $p=p_{ref}$ has been set at the boundaries located at $x/D$=-50 and $x/D$=50, respectively, while the symmetry condition has been used at the top and bottom boundaries.

The drag and lift coefficients, denoted as $C_D$ and $C_L$, have been defined as:

$$C_D={\displaystyle \frac{F_x}{\frac{1}{2}{\rho }_l{U_{ref}}^{2}D}}$$

(40)

$$C_L={\displaystyle \frac{F_y}{\frac{1}{2}{\rho }_l{ref}^{2}D}}$$

(41)

where $F_x$ and $F_y$ are the streamwise and transverse components of the force acting on the cylinder surface, respectively, and $U_{ref}$ is the free-stream velocity.

The unsteady vortex shedding case has been calculated for liquid water at $Re$=200 as in previous studies [36,37,38,39,40,19,41,42]. Vortices are generated on the cylinder surface and detach alternatively from its upper and lower surfaces which results in fluctuations of $C_D$ and $C_L$ as shown in Figure 4. The time-average value of $C_D$, denoted as $C_{D,av}$, the maximum value of $C_L$, denoted as $C_{L,max}$, and the $St$ value are listed in Table 1 and compared with published data obtained by different authors using different meshes and methods. It is confirmed that the present results match well with the reference values.

Simulations are extended to cavitation conditions at $\sigma$= 1.0. As shown in Figure 5, the oscillating flow around the cylinder results in $C_D$ and $C_L$ fluctuations. The corresponding $C_{D,av}$, $C_{L,max}$, and $St$ values are listed in Table 2 and compared with previous results obtained by other authors. As the liquid pressure decreases to $p_v$, cavitation occurs around the lateral surface of the cylinder and then the vaporous bubbles detach from the cylinder surface and travel within the shed vortices. The positive velocity divergence caused by the appearance of cavitation leads to the vorticity dilatation, thus elongating and weakening the shed vortices [19]. Therefore, the computed shedding frequencies are reduced to $St$= 0.17, which are comparable to the corresponding predictions of 0.16 reported by Seo et al. [38] and Gnanaskandan and Mahesh [19], and of 0.177 reported by Hong and Son [42]. Moreover, Figure 6 compares the instantaneous vorticity contour colored with vapor volume fraction at $\sigma$=1.0 obtained by Gnanaskandan and Mahesh [19] and the present method, showing similar results.

5. Results

5.1. Computational Domain and Boundary Conditions

Figure 7 shows the dimensions of the computational domain to simulate the flow around a wedge which are the same as the ones used in the experimental setup presented by Wu et al. [5] . The cross-section of the tunnel is rectangular and its width is 4D, where D is the height of the wedge base which equal to 0.019 m. The span of the wedge is 4D. For the incompressible cavitation case, the inlet is located at 13.5D upstream from the wedge trailing edge plane and the outlet is at 19.5D downstream from it. For the compressible cavitating simulation the fluid domain is extended to implement the sponge layers at the inlet and outlet boundaries which are located 25.9D upstream and 36.6 D downstream of the trailing edge, respectively.

The flow around the wedge has been simulated using the DDES model to capture the globally unstable flow characteristics. Such flow has the feature that the turbulence in the separated zone is independent of the turbulence within the incoming attached boundary layer [43]. This flow feature can require to reduce significantly the mesh resolution. Following the suggestions provided by ANSYS Fluent^® [43], more than 20 cells per characteristic length, D, are sufficient to resolve the main flow instability. Furthermore, an isotropic cell (Δx = Δy = Δz) has been used in the area near the trailing edge of the wedge to avoid the numerical error resulting from cells with too large aspect ratios. To assess the effects of mesh refinement on the calculations, three meshes with different refinements have been tested as listed in Table 3. The topology of mesh M1 around the wedge is depicted in Figure 8(a) and 8(b).

For the non-cavitating case, the inlet boundary condition has been set as a uniform inflow velocity, $U_{ref}$ =6 m/s, and the outlet type boundary condition has been set as shown in Figure 9(a). For the cavitation case at $\sigma$= 1.9, the corresponding total pressure has been specified at the inlet boundary condition and the fixed mass flowrate boundary condition has been fixed at the outlet. To absorb the pressure wave generated by the compressible fluid and avoid reflections, sponge layer boundaries have been added at the inlet and outlet areas [34], as depicted in Figure 9(b). The specific values of the different boundary conditions for the two operating conditions without and with cavitation are detailed in Table 4.

All the simulations have been run with ANSYS Fluent^® using the Bounded Central Difference (BCD) advected scheme and a time step of $\Delta t$ =5×10^-5 s where the corresponding advected Courant number is below 1 behind the wedge.

5.2. Verification and Validation of the Non-Cavitating Case

Figure 10 shows the spectra of the simulated lift coefficient calculated using FFT for the three different grid refinements at non-cavitation conditions with the compressible model. It can be seen that for all the meshes the same maximum frequency peak is obtained with similar amplitude even when the coarse mesh is used. Therefore, to save computational resources, the mesh M1 has been selected and used for the following analysis.

As presented in Figure 10, the main frequency peak corresponding to the vortex shedding frequency is found around 91 Hz which y is in good agreement with the value of 90 Hz obtained by Wu et al. [5].

The numerically obtained $C_p$ on the surface of the wedge is plotted in Figure 11. Here, the position at $y=\pm 0.5D$ corresponds to the two vertices at the truncated trailing edge. It can be seen that, at the centerline of the base of the wedge, the value of $C_p$ = -1.5 is exactly the same as the experimental value reported by Wu et al. [5]. Therefore, it can be concluded that, in the non-cavitating case, the current numerical solution of the unsteady flow field provides a reasonable resemblance to the flow conditions measured during the experiment carried out by Wu et al. [5].

5.3. Assessment of the Compressible Cavitation Model

Using the validated incompressible and compressible modelling approaches presented in the previous sections, the cavitating wake behind the wedge at $U_{ref}$=6 m/s and $\sigma$= 1.9, has been calculated numerically to demonstrate the effects of fluid compressibility on the dynamics of the cavitating vortex shedding behind the wedge.

5.3.1. Pressure on the Wedge Surface

Figure 12 compares the mean values of $C_p$ on the wedge surface using both the incompressible and the compressible cavitation solvers and it can be seen that the results are exactly the same. Moreover, , the numerical $C_p$ of -1.4 at the centreline of the base of the wedge is in good agreement with the experimental value of -1.45 obtained by Wu et al. [5].

5.3.2. Unsteady Loads on the Wedge Surface

For comparison, Figure 13 shows the time histories of $C_L$ and their spectra with the incompressible and compressible cavitation solvers for $\sigma$= 1.9. It can be seen that both solvers predict the same value of 210 Hz for the vortex shedding frequency behind the wedge. Such result is in good agreement with the dominant frequency reported by Wu et al. As shown in Figure 13(b), the shape and amplitudes of the spectra obtained with the incompressible and compressible cavitation solver are quite similar for the lower frequencies below 800 Hz. However, at high frequencies above 1000 Hz, the spectra obtained with the incompressible cavitation solver present lower amplitude values than with the compressible cavitation solver. It is interesting to note that similar results have also been observed in the investigation of cloud cavitation around the hydrofoil conducted by Wang et al. [20].

5.3.3. Cavitation Structures

Figure 14 shows visualization of the 3D cavitation structures inside the shed vortices using iso-surface plots with the numerical results obtained with both the incompressible and the compressible cavitation solvers. And they are compared in the same figure with the photographs taken during the experimental tests carried out by Wu et al. [5]. It can be seen that the first few pairs of spanwise vortex cores (identified with the $Q$ criterion) are filled with vapor, and that these cavity structures (represented with iso-surfaces of ${\alpha }_v$ = 0.05) are advected downstream. Furthermore, the relative positions of the vortices and the spacings obtained with the incompressible and compressible cavitation models are quite similar and they align well with the experimental results.

Figure 15 shows the simulated instantaneous field of ${\alpha }_v$in the cavitating wake at different instants of time obtained with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers at $\sigma$ = 1.9 as well as the experimental (Exp.) images reported by Wu et al. [5] using time-resolved X-ray densitometry. All of them show clearly the periodically shed vortices at the vortex shedding frequency. The evolution of the cavitation structures with time confirms that the periodic shedding is governed by the alternating shedding vortex street behind the wedge. In each shedding period, a cavity starts to form and fill the center of the shed vortex. Then, the cavitating vortex is advected from the attached boundary layer and shed from the trailing edge. The alternating cavitating shedding vortices form the vortex street behind the wedge. The comparison between numerical results using incompressible and compressible cavitation solvers and the experimental results indicates that both cavitation solvers can capture the main flow characteristics in terms of the cavities morphology and their relative position as it can be seen in Figure 15(a) and Figure 15(c).

Figure 16 shows the time average void fraction field calculated with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers and the experimental results. It can be seen that the simulations with two different cavitation solvers produce almost indistinguishable mean and RMS void fraction fields behind the wedge. The shape of mean cavity structures obtained numerically is comparable to the pattern observed in the experimental image. Note that, there is a downstream offset in the mean and RMS void fraction fields between the numerical and the experimental results that may be due to the low numerical resolution of the locally unstable flow structures within the detachment area.

6. Conclusions

In this study, the numerical results predicted by with an incompressible and a compressible cavitation model have been compared to assess the effects of the fluid compressibility on the characteristics and dynamics of the cavitating flow behind a wedge. For validation, the experimental results obtained by Wu et al. [5] for a cavitating wake behind a wedge have been taken as a reference. It has been found that both solvers capture the dynamics of the cavitating vortex street and give similar results in terms of mean pressure, lift and drag forces as the experimental ones. The main difference has been found on the frequency content of the fluid force fluctuations:

• For low frequencies up to 800 Hz, both solvers predict similar amplitudes.
• For higher frequencies above 1000 Hz, the compressible solver predicts larger amplitudes.

Moreover, regarding the prediction of the dominant vortex shedding frequency and the instantaneous and mean void fraction fields, the simulations performed with both solvers have provided almost identical results.

To summarize, it has been confirmed that the proposed compressibility solver is able to predict the cavitating vortex shedding behind a wedge.