The characteristics of unsteady cavitation around a sphere

Cavitation can be observed in many engineering systems, such as propeller, hydraulic machinery, and underwater bodies. The unsteady characteristics of cavitation have been the subjects of some research studies.^1–4 Kubota et al.¹ measured the unsteady flow structure around the cloud cavitation over a 2D (two-dimensional) hydrofoil. It can be seen that cloud cavitation exhibits obvious quasi-periodicity behaviors. Le et al.⁴ performed a large number of cavitation experiments on the flat-concave hydrofoil and analyzed the cavitation shedding characteristics from the perspective of the cavitation thickness. Furness and Hutton⁵ recorded the growth and collapse of fixed cavities in a two-dimensional convergent-divergent nozzle. They propose that the re-entrant jet is the main mechanism of sheet cavitation instability. Kawanami et al.⁶ further confirmed the effect of a re-entrant jet on cavity shedding by pressure measurement and studied the inhibitory effect of different obstacles on cavitation. Reisman et al.⁷ further discussed the effect of a shockwave on cavitation instability by their experimental or numerical studies.

Unsteady cavitation typically includes three processes: one is the appearance and growth of a cavity, the second is the instability and shedding of a cavity, and the third is the collapse of cloud cavitation and regrowth of the attached cavitation. Gosset et al.⁸ simulated the unsteady cavitating flow over a NACA0015 foil by using the implicit large eddy simulation (ILES) and the k − ω models. They captured the periodic phenomena related to the cavitation shedding. Their results showed that the cavitation shedding exhibits a single-frequency feature. The corresponding Strouhal number is about 0.2. Recently, Hsiao et al.⁹ employed a multiscale two-phase flow model to study the unsteady results and also showed a single-frequency feature for different angles of attack.

The cavitating flow over a sphere is more complicated because of the interaction between wake flow and cavitation. Brennen experimentally studied the stability of the cavity interface in the fully developed cavitating flow around a sphere.^10,11 Arakeri¹² discussed the effect of a laminar boundary-layer on the cavitation inception. A conditional semi-empirical approach was developed to predict the cavitation inception position on a smooth body. Recent work by Brandner et al.¹³ observed the cavitation over a sphere in a water tunnel under a turbulent condition (R_e = 1.5 × 10⁶). They found that the flow is dominant by a dual-frequency mode in the cases of small cavitation numbers.¹⁴ Pendar and Roohi¹⁵ performed numerical simulations under the same conditions as Brander’s. However, they did not analyze the frequency features associated with the periodic behaviors.

It is obvious that the cavitation over a 3D (three-dimensional) sphere is more complicated than over a 2D hydrofoil. The dual-frequency mode, found in the experiment,¹⁴ is interesting but not fully understood. In this paper, we numerically investigate the cavitating flow around a sphere, with special emphasis on the frequency analysis. The used cavitation number covers the cavitation states from the inception to supercavitation. The delayed detached eddy simulation based on Spalart-Allmaras (S-A DDES) turbulence model and the Sauer cavitation model are used to model the cavitating flow, then analyze the unsteady behavior of the cavitation, and finally discuss the frequency features and their reasons.

The whole numerical simulation is carried out on the numerical platform OpenFOAM, which is an open-source software library for solving partial differential equations. The two-phase flow solver of the OpenFOAM package, interPhaseChangeFoam, is employed. InterPhaseChangeFoam is a solver for two incompressible, isothermal immiscible fluids with phase-change and uses a VOF (volume of fluid) phase-fraction based interface capturing approach. The momentum and other fluid properties are of the “mixture,” and a single momentum equation is solved. The set of phase-change models provided are designed to simulate cavitation, but other mechanisms of phase-change are supported within this solver framework. Turbulence modeling is generic, i.e., laminar, Reynolds Averaged Simulation (RAS) method, or Large Eddy Simulation (LES) may be selected.

Assuming that the fluid is Newtonian and incompressible, the homogeneous mixture multiphase flow is governed by the incompressible Navier-Stokes (NS) equations as follows:

in which F is the quality force, p is pressure, μ is the mixture viscosity, and ρ is the mixture density. The physical properties of the mixture are weighted by the physical properties of each phase. The mixture density ρ and viscosity μ are computed as follows:

The volume fraction of liquid phase α_l and vapor phase α_v should satisfy the condition as follows:

The evolution of volume fraction is governed by a transport equation

where $ṁ$ represents the phase change rate.

In recent years, the transport equation-based model (TEM) and Equation of State (EOS) model have been used widely to calculate the density of mixture in cavitating flow. In an EOS model, a density-pressure dependency is needed to solve the rapid phase change process between two phases. Numerical methods for EOS models always refer to the density-based algorithm which is commonly employed in aerodynamics computations. TEM employs a phase transfer equation which reflects the mass conservation law. Typically, a pressure-based algorithm is used for the computation of cavitating flows with TEMs. There are several cavitation models based on the Rayleigh-Plesset equation of the spherical bubble theory, such as Singhal’s full cavitation model¹⁶ and Sauer’s volume fraction model.¹⁷ Sauer’s model can be given by

where

The initial number of bubbles (n) is set as 1.6 × 10¹³.

The Reynolds number of cavitating flow is generally high. When it is impossible to perform a direct numerical simulation (DNS), the governing equations need to be averaged or filtered and closed with the corresponding turbulence models. Here, we use the S-A DDES turbulence model in all cases. The S-A turbulence model, first proposed by Spalart and Allmaras,¹⁸ establishes a transport equation for the new eddy viscosity $ν̃$

and the relationship between the various parameters

where ν and ν_t are different viscosity coefficients, which in turn represent the molecular viscosity coefficient and fluid eddy viscosity coefficient. Here, S is chosen as the magnitude of the vorticity |ω|. The left side of Eq. (9) represents the time term and the convection term of the transport variable $ν̃$⁠, and the right side represents the production item, the diffusion term, and the destruction term, respectively. The destruction term represents the reduction of the stresses in the vicinity of solid walls and the dissipation of the turbulent kinetic energy in the near-wall region, which is directly related to the reciprocal of the minimal wall distance d in the Reynolds Average Navier-Stokes (RANS) mode. In the interPhaseChangeFoam solver, the corresponding constants are κ = 0.041, σ = 2/3, c_b1 = 0.1355, c_b2 = 0.622, c_$v$1 = 7.1, c_$w$1 = c_b1/κ² + (1 + c_b2)/σ, c_$w$2 = 0.3, c_$w$3 = 2.

The above is the standard RANS formulation of the S-A turbulence model. Spalart and Allmaras modified the model and proposed the detached eddy simulation (DES) method.¹⁹ They replaced the wall distance d with a new wall variable l_DES, where C_DES is constant (C_DES = 0.65) and Δ is the filter width in large eddy simulations. They satisfy the following relationships:

The basic idea of delayed detached eddy simulation (DDES) is to construct a function to determine the length variable l_DES in DES, which can delay the switching of the RANS mode to the LES mode in the boundary layer, and ensure that the RANS mode is used in the boundary layer.²⁰ The redefined length variable l_DES is

where ν represents the molecular viscosity, U_i,j is the velocity gradient, f_d is the delaying function, and r_d is a dimensionless length variable, which is the square of the ratio of the length scale of the model to the distance from the wall. When applying the S-A model, ν_t + ν can be replaced with $ν̃$⁠.

We simulate the cavitation flow in a domain like the test section of a water tunnel. Figure 1 shows a schematic of the computational domain and the boundary conditions imposed on the domain. The mean flow is along the y-axis, and the sphere is placed at the centerline of the water tunnel. The domain is 2.6 m long, and the diameter of the sphere D is 0.15 m. The upstream distance is 0.75 m, while the downstream distance is 1.85 m. Classical boundary conditions, imposed inlet velocity (u = 9 m/s), and fixed outlet static pressure are applied. We can change the reference static pressure p to adjust the cavitation number. The spherical surface uses the no-slip wall condition as well as the side of the domain. The Reynolds number is defined as R_e = uD/ν = 1.5 × 10⁶. The cavitation number is defined as

where p is the outlet static pressure, p_sat is the saturation pressure, and ρ is the water density. In all cases, the fluid properties are specified as ρ_l = 1000 kg/m³, ρ_$v$ = 0.02308 kg/m³, μ_l = 9 × 10⁻⁴ kg/(m s), μ_$v$ = 1 × 10⁻⁵ kg/(m s), and p_$v$ = 2300 Pa.

As shown in Fig. 2, the total number of grids is about 4 × 10⁶. In order to better study cavitation and wake flow, a special grid refinement is applied in the vicinity of the sphere and the centerline of the computing domain. For a typical simulation case at σ = 0.5, the mean value of y⁺ is 34.5 and the maximum value of y⁺ is 65.8 around the sphere. Considering the possible influence of the boundary layer of the wall on the flow past the sphere, the mesh at the wall is also refined and the mean value of y⁺ is 82.1.

First, the time-average properties of the flow are investigated. Due to the periodic behaviors of cavitation, the statistical results will be stable after several periods. We obtain the results by averaging more than 20 periods. In all cases, we chose the same isosurface with the liquid volume fraction of 0.9. Figure 3 shows the changes of cavitation morphology in a range of cavitation numbers. The first column of Fig. 3 presents the instantaneous experimental results.^13,14 The second presents the instantaneous simulation results, and the third shows the time-averaged simulation results. As the number of cavitation increases, the volume of the cavitation decreases, and the leading edge position of cavitation continues to move upstream. The simulation results fit well with the experimental data. We performed a statistical analysis of the time-averaged cavity length and the position of the cavity leading edge. To describe the inception position, we define a position angle β, which is shown in Fig. 2(d).

Figure 4 represents the inception position of the cavitation at various numbers. It can be seen that the position angle increases approximately linearly with the cavitation number. The maximum angle is about 95° at σ = 1.0. When the cavitation number comes to 0.3, the inception position angle drops to 68°. The simulated results are in good agreement with the experimental data.^21,13

In 1956, Garabedian²² proposed a formula to predict the length of the axisymmetric cavity

As shown in Fig. 5, we compare the numerical results with Garabedian’s formula. Although this formula aims at predicting the length of supercavity, we find that it is also applicable to the cases of partial cavitation. The cavity length increases continuously as the cavitation number decreases. The length change is relatively sensitive in the ranges of small cavitation numbers. In general, the cavity morphology is strongly affected by the cavitation number. When the cavitation number is large, the cavity is short. The cavitation inception position is located behind the peak-point on the sphere surface. With the cavity length increasing, the inception position moves upstream. This is an interesting flow pattern differing from the single-phase flow.

In this section, the transient properties of the flow are investigated. In Graff’s experiment, they measured the pressure signal at a downstream point on the sphere surface. This signal can be used to validate the transient result of simulation. Figure 6(a) shows the comparison of the surface pressure change over time at the same position with Graff’s experiment,¹⁴ where the time, t, has been nondimensionalized by t′ = Ut/D and the pressure, p, has been subtracted by the mean value and then divided by the standard deviation. As shown in Figs. 6(a) and 6(b), the numerical pressure evolution and Fast Fourier Transform (FFT) results fit well with the experimental data in terms of amplitude and period.

Figure 7 represents a typical cycle of cavity shedding at the cavitation number of 0.95. Cavitation keeps growing from time T₀ to T₄. During this process, a re-entrant jet occurs at the tail part of the cavity. This jet moves upstream and finally cuts off the cavity from the surface at time T₆. Then, the volume of cavitation decreases when moving downstream. This process shows the quasiperiodic behavior during the cavitation evolution. We try to find the frequency feature by the analysis of lift force curve. As shown in Fig. 8, after the Fast Fourier Transform (FFT), we obtain the frequency spectrum. However, the basic frequency is about 11 Hz. This frequency does not coincide with the period of cavitation shedding but corresponds to the frequency of wake flow. That means, the frequency at this condition is close to the frequency in single-phase flow, which is about 9 Hz [see Fig. 8(b)]. The effect of unsteady cavitation is weak since the cavity length is short.

As the cavitation number decreases, the length of the cavitation increases. Figure 9 shows a typical shedding cycle at the cavitation number of 0.75. Cavitation keeps growing from time t₀ to t₃. During this process, a re-entrant jet occurs at the tail part of the cavity. The re-entrant jet moves upstream and interacts with the surface of the cavity. Finally, the jet cuts off the cavity from the surface of the sphere at time t₆. The frequency spectrum analysis of the lift force is given in Fig. 10. Under this condition, two frequencies appear in the distribution. The low frequency is close to that in the case of σ = 0.95. The high frequency is about 35 Hz. It seems that the intensity level of cavitation has an influence on the flow instability. We will discuss the dual-frequency feature in detail in Sec. IV C.

In this section, we focus on the frequency analysis of the lift force in the range of simulated cavitation numbers. A typical St number is defined based on the diameter of the sphere. To describe the instability related to the cavitation, a St number based on the cavity length is also defined.

Figure 12(a) shows typical St numbers (St = fD/U) based on the diameter of the sphere at various cavitation numbers. Results show that the flow instability exhibits two modes. When the cavitation number is larger than 0.8, the cavitation has little effect on the flow. Just like single-phase flow, the St number exhibits a single-frequency mode. As the cavitation number is reduced, cavitation has a stronger effect on the flow. There is a step change around the cavitation number of 0.8. When the cavitation number drops below 0.8, the flow exhibits another dual-frequency mode. Then as the cavitation number decreases, both frequencies decreases continuously.

As shown in Fig. 10, the flow obviously exhibits a dual-frequency mode in frequencies of 11 Hz and 35 Hz at σ = 0.75. Figure 11 shows a typical process of a large-scale cavity shedding. Cavitation keeps growing from time t₀ to t₃. During this process, it can be clearly seen that a re-entrant jet occurs at the tail part of cavity and moves upstream. The re-entrant jet impinges on the cavity at time t₃ but only causes a small-scale cavity shedding. Simultaneously, the cavitation also keeps growing. The second re-entrant jet occurs at time t₄ and interacts with the cavity surface again at time t₅. The cavity is cut off at time t₆. Then, the shedding cavitation moves downstream and the attached cavitation regenerates at time t₇. With the large-scale cavity shedding, the lift force has a strong oscillation. We find that the frequency of large-scale cavity shedding corresponds to the high frequency in a dual-frequency mode. Combined with results shown in Fig. 12(a), the low frequency corresponds to the natural frequency of the sphere. For σ ≥ 0.8, the cavitation has a little effect on the flow. The frequency feature at this condition is close to the feature in single-phase flow and keeps a single-frequency mode. For σ < 0.8, the cavitation has a non-negligible influence on the flow. The flow exhibits a dual-frequency mode. The low frequency corresponds to the natural frequency of the sphere, while the high frequency reflects the dynamic evolution of cavitation. The natural frequency changes continuously because the wake flow is simultaneously affected by the shape of cavity.

For σ < 0.8, the cavitation has a non-negligible influence on the flow. To describe the instability related to the cavitation, a St number based on the cavity length in this range of cavitation number is used. As shown in Fig. 12(b), the frequency is almost stable at two constant numbers of 0.2 and 0.5. Combined with Garabedian’s theoretical formula for describing the cavity length,²² we proposed a formula to predict the frequencies for the dual-frequency mode

where St₁ is the stable St number based on the cavity length of dual-frequency mode.

The unsteady cavitating flow past a sphere was numerically investigated. The SA-DDES turbulence model and the Sauer cavitation model were used in the modeling of the flow. The used cavitation numbers cover the cavitation states from inception to supercavitation. The averaged and transient properties were analyzed and compared to the existing data. Results show that with the cavitation number decreasing, the cavity length increases, while the cavitation inception position moves upstream along the sphere surface. The frequency feature was analyzed by using the FFT approach. When the cavitation number is larger than 0.8, the cavitation shedding has little effect on the lift force. The flow exhibits a single-frequency mode. When the cavitation number drops below 0.8, the flow exhibits a dual-frequency mode. The low frequency corresponds to the natural frequency of the flow, while the high frequency reflects the dynamic evolution of cavitation. In this range of cavitation numbers, the natural frequency changes continuously because the wake flow is simultaneously affected by the shape of cavity. We found that the St numbers based on the cavity length are almost stable at two constant numbers of 0.2 and 0.5. Combined with the formula describing the cavity length, we proposed a formula to predict the frequencies for the dual-frequency mode.

This project was supported by the National Natural Science Foundation of China (No. 11772298), the State Key Program of National Natural Science of China (No. 91852204), and the Fundamental Research Funds for the Central Universities (No. 2018QNA4055).