Ocean Engineering 156 (2018) 424–434

Contents lists available at ScienceDirect

Ocean Engineering
journal homepage: www.elsevier.com/locate/oceaneng

Numerical investigation of cavitation vortex dynamics in unsteady

cavitating flow with shock wave propagation
Changchang Wang a, Qin Wu a, b, *, Biao Huang a, Guoyu Wang a
a
School of Mechanical Engineering, Beijing Institute of Technology, Beijing, 100081, China
b
Department of Thermal Engineering, Tsinghua University, Beijing, 100084, China

A R T I C L E I N F O A B S T R A C T

Keywords: The objective of this paper is to study the cavitation vortex dynamics in compressible turbulent cavitating flow,
Cavitating flows around a NACA66 hydrofoil. The simulations are conducted based on the open source software OpenFOAM,
Cavitation vortex dynamics solving the compressible governing equations together with the Tait state equation for water and ideal gas state
Compressible cavitating solver equation for vapor. The Saito cavitation model is used to model the cavitation phase change process and SAS SST
Shock wave turbulence model is adopted to account for the turbulence effects. The numerical results showed a good agree-
OpenFOAM
ment with the experiments. The cavity evolution presents quasi-periodic behaviors, and the alternative re-entrant
flow movement and shock wave dynamics play an important role in cavitation vortex dynamics. Strong coherent
relationship between cavitation behaviors and vortex dynamics is illustrated. Based on the budget analysis of
vorticity transport equation, it shows that during the attached cavity growth, the re-entrant flow development and
the cloud cavity being shed stages, the baroclinic torque term dominates the vorticity transport process. During
the shock wave propagation stage, the vortex dilatation term dominates in the attached cavity sheet region.
Besides, the baroclinic torque term has strong influence on vorticity transport characteristics.

1. Introduction coherent interaction between the shock phenomena and the multiphase
compressible flow dynamics.
Hydrodynamic cavitation is generally undesirable phenomena due to Numerous researchers have experimentally measured the pressure
its severe damage to the hydraulic systems, such as pumps, turbines, pulses in the cavitating flows (Avellan and Farhat, 1989; Le et al., 1994;
marine propellers and rocket propulsion inducers in fluid machinery and Mckenney and Brennen, 1994; Reisman et al., 1994), and investigated
spillway in hydropower station (Knapp et al., 1970; Brennen, 1995; the shock wave propagation characteristics in the bubbly mixture
Wang et al., 2001). The very destructive effects, such as pressure fluc- (Noordij and Wijngaarden, 1974; Mørch, 1980, 1981; Hanson et al.,
tuations, variation, noise and erosion, are suffered in the condition of 1981). Reisman and Brennen (1996) conducted experiments to investi-
sheet/cloud cavitation, which is characterized by periodic breakup of gate the pressure impulses emitted in the cloud cavitation systematically,
attached cavity sheet, shedding and collapse of cloud cavity (Soyama using piezo-electric transducers on both stationary and oscillating hy-
et al., 1992; Joseph, 1995; Reisman, 1996; Chen et al., 2015; Wang et al., drofoils. They observed that pressure pulses are highly related with the
2017b; Wu et al., 2015, 2018). Many researchers have conducted in- large low void fraction region, which is supposed to be caused by the
vestigations for the shedding mechanism of the cloud cavity, among cavity collapse event. The impulse intensity was enhanced with the
which, the shock wave propagation during the cloud cavity collapse is increasing of inlet speed and the decreasing of cavitation number. After
supposed to be the main region for the cloud cavitation instabilities that, they further concluded that the shock wave dynamics led to the
(Reisman et al., 1998; Leroux et al., 2005; Budich et al., 2015; Wang damage and noise in cavitating flows. Arndt et al. (2001) applied a
et al., 2017a). It has also been pointed out that a marked decrease in the combined experimental and numerical method to investigate the
sound speed phenomena appears in liquid/gas mixture by Wallis (1967) sheet/cloud cavitation dynamics around a NACA0015 hydrofoil. They
and Brennen (1995), and thus it is very likely that the local liquid/gas concluded that at low value of σ /2α, the bubbly shock mechanism is the
flows reaches to transonic or supersonic speed even when the bulk flow cause of the cloud cavitation generation. Leroux et al. (2004) employed
velocity is relatively low. Hence, it is necessary to understand the the wall-pressure transducers to investigate partial cavity instabilities

* Corresponding author. School of Mechanical Engineering, Beijing Institute of Technology, Beijing, 100081. China.
E-mail address: wuqin919@163.com (Q. Wu).

https://doi.org/10.1016/j.oceaneng.2018.03.011
Received 11 October 2017; Received in revised form 7 February 2018; Accepted 3 March 2018

0029-8018/© 2018 Elsevier Ltd. All rights reserved.

C. Wang et al. Ocean Engineering 156 (2018) 424–434

around a NACA66 hydrofoil. They showed that at the condition of cavitation region is a high vorticity region and the cavity shedding pro-
σ ¼ 1.25 at attack of angle AoA ¼ 6 , the cloud cavity collapse induced cess is strongly related to the interaction between the cavity and vortex
shock wave would significantly cause the newly growth attached sheet structures. Although much efforts has been made to investigate the
cavity being shorten and collapsed, and thus resulting in the cavitation vortex dynamics in cavitating flows (Huang et al., 2013a, b; Ji et al.,
cycle increase and the cavitation Strouhal number reduction. Recently, 2014, 2015; Wang et al., 2016; Chen et al., 2017; Long et al., 2017a),
Ganesh et al. (2016) used high speed video and time-resolved X-ray there exist very limited studies of the cavitation vortex dynamics in such
densitometry to investigate the sheet to cloud cavitation transition complex cavitating flow accompanied by the re-entrant flow and the
mechanism in a venturi section. The shock wave propagation in the shock wave propagation.
attached cavity was well observed and it was concluded that the shock The present study applied the compressible cavitating flow solver to
wave propagation mechanism dominates the periodic large scale cloud study the unsteady sheet/cloud cavity dynamics around a NACA 66 hy-
cavity shedding process. Wu et al. (2016a) used hydrophone measure- drofoil. The Tait state equation for water and ideal gas state equation for
ments accompanied with high-speed videos and time-resolved X-ray vapor are employed to consider the compressibility effects. The objective
densitometry to investigate the partial cavity shedding dynamics around of this paper are to (1) illustrate the unsteady flow structures during the
a NACA0015 hydrofoil. They observed that with the reduction of the unsteady sheet/cloud cavity evolution, (2) study the cavitation vortex
cavitation number, the cavity dynamics changed from the pinch-off from dynamics, especially under the influence of the cloud cavity collapse
the rear of the cavity caused by the re-entrant jet, to a propagation bubbly induced shock wave propagation.
shockwave mechanism caused by the collapse of the shed cloud cavity.
With the development of the numerical technique, most past simu- 2. Physical and mathematical model
lations adopt the incompressible method in the cavitation calculation
(Bensow and Bark, 2010; Li et al., 2008; Huang et al., 2014; Wu et al., 2.1. Governing equations
2016b; Yu et al., 2017). Considering that the shock wave phenomena in
complex cavitating flows is highly related to the compressible charac- The set of governing equations based on homogeneous multiphase
teristics, it is necessary to employ the compressible solver into the cavi- flows strategy consist of the three-dimensional compressible Navier-
tating flow numerical simulations (Kunz et al., 2001; Sezal et al., 2009; Stokes equations (continuity, momentum and energy equations) along
Haren et al., 2016; Gnanaskandan and Mahesh, 2016; Egerer et al., with a transport equation for the void fraction of vapor:
2016). The solving of compressible multiphase equations is a challenging
problem, carrying a significant computational burden and requires ∂ρm
þ r⋅ðρm UÞ ¼ 0 (1)
careful treatment, in which the selection of time step is critical to guar- ∂t
antee the convergence and should satisfy the maximum Courant number


∂ðρm UÞ 2
criteria. Venkateswaran et al. (2002) applied a preconditioned þ r⋅ðρm UUÞ ¼ rp þ r⋅ μm rU þ ðrUÞT  ðr⋅UÞI (2)
time-marching algorithm to investigate the compressibility effects of ∂t 3
cavitating flow. The results showed that the method considering the  
compressibility can improve the cavitation dynamics compared with the ∂ ρm Cp T  
þ r⋅ ρm UCp T ¼ r⋅ðκrTÞ (3)
previous incompressible computations. Saito et al. (2007) developed the ∂t
compressible two-phase Navier-Stokes equations and applied them to
investigate the unsteady behaviors of cloud cavitation around a NACA ∂ρl αl
þ r⋅ðαl ρl UÞ ¼ mþ þ m (4)
0015 hydrofoil. The method accurately captured the U-shape cloud ∂t
cavity which is supposed to be related to the cavity collapse event.
Schnerr et al. (2008) developed the compressible cavitating flow solver ρm ¼ αρv þ ð1  αÞρl (5)
to investigate the collapse induced shock wave dynamics with nano-
second time step. The pressure peak up to 230 bar was calculated as the μm ¼ αμv þ ð1  αÞμl (6)
shock wave was induced by the cloud cavity collapse. However, inves-
In the above equations, ρm is the mixture density, ρl is the liquid
tigation on the cavitation vortex dynamics concerning with the cavitation
density, ρv is the vapor density, α is the volume fraction, U is the velocity
compressible characteristics, especially in the complex cavitating flows
tensor, I is the unit tensor, αl is the liquid fraction, αv is the vapor fraction,
along with the shock wave dynamics, has been rarely conducted. The
μm is the mixture dynamic viscosity, μl is the liquid dynamic viscosity, μv
turbulence model plays an important role in high Reynolds number
is the vapor dynamic viscosity, Cp is the specific heat, T is the tempera-
cavitating flows. Although the RANS approach has been widely used to
ture, k is the thermal conductivity, mþ is the evaporation rate, m is the
simulate turbulent cavitating flows, it has the weakness of the
condensation rate, respectively. The dynamic viscosity of water
over-prediction of turbulent eddy viscosity in multiphase mixture region
(μl ¼ 1.0018  103 Pa s) and vapor (μv ¼ 9.7271  106 Pa s) are
(Coutier-Degosha et al., 2003; Reboud et al., 2008; Huang et al. 2013a,
assumed constant at 293.15 K due to the weak thermal effects in room
2013b; Long et al., 2017b). The LES approach requires much more
temperature. The subscripts m, l and v denote the mixture, liquid and
computational resources (Wang and Ostoja-Stazewski, 2007; Ji et al.,
vapor respectively.
2015; Huang et al., 2014). Thus, some hybrid RANS-LES turbulence
In order to consider the compressibility within the cavitating flows,
models have been developed (Ducoin et al., 2012; Huang et al., 2013b;
the state of equation for water/vapor is used to update the density field in
Chen et al., 2016). Different from the grid-based hybrid turbulence mode,
the flow. The following Tait equation of state for water (Core, 1948) is
the RANS and LES transition of the SST SAS (Scale-adaptive simulation)
employed
turbulence model is based on local turbulence structures von Karman
length-scale. Considering the multi-scale characteristics in both time and
N
pl þ B ρl
space domain in unsteady cavitating flows, the SST SAS can solve more ¼ (7)
plsat þ B ρlsat
smaller turbulence eddy structures and better capture the unsteady in-
stabilities involved in the multi-scale cavitation dynamics. Most of the where psat ¼ 2338.6 Pa and ρsat ¼ 998.16 kg/m3 are the saturation pres-
experimental and numerical researches have mainly focused on the sure and saturation density of liquid water at 293.15 K according to
cavity patterns and pressure characteristics during the shock wave for- National Institute of Standards and Technology (NIST) data.
mation and propagation in cloud cavitating flow. However, the B ¼ 3.06  108 Pa and N ¼ 7.1 are the fitted constants.

425
C. Wang et al. Ocean Engineering 156 (2018) 424–434

2.2. The ideal gas equation of state is applied for vapor Where Pk ¼ τij ∂Ui =∂xj is the production of turbulent kinetic energy, and
CDkω ¼ ∂k=∂xj ∂ω=∂xj is the turbulent cross-diffusion term, and ε ¼ β⋅kω
is the dissipation rate. F1 is the blending function between the ω equation
pv ¼ ρv Rv Tv (8)
ε equation and F2 is another blending function. The turbulent viscosity is
pffiffiffiffiffiffiffiffiffiffi
where subscript v denotes vapor-phase value and Rv ¼ 461.6 J/(kg.K) is given by νT ¼ a1 k=maxða1 ω;SF2 Þ, S ¼ Sij Sij . The SAS term is defined as
gas constant. In the present study, the non-condensable gas is ignored in following
the gas phase. "
2
 #
2 L 2ρk 1 ∂k ∂k 1 ∂ω ∂ω
QSAS ¼ max ρζ2 S  CSAS max 2 ; ;0
2.3. Saito cavitation model LvK σΦ k ∂xj ∂xj ω2 ∂xj ∂xj
(14)
The cavitation model used in the present work was developed by

Saito et al. (2007). This cavitation model is a transport-equation based

∂U=∂y

Lvk ¼ κ

2
(15)
cavitation model that ∂ U=∂y

maxððp  pv Þ; 0Þ pffiffiffi 1=4

m_  ¼ Cc α2 ð1  αÞ2 pffiffiffiffiffiffiffiffiffiffiffiffiffi ; if p > pv (9) Where ζ2 ¼ 3.51, σ Ф ¼ 2/3, CSAS ¼ 2.0, L ¼ k=ðcμ ⋅ωÞ is the length
2 π Rg T
2
scale of the modeled turbulence. When the term ρζ 2 S2 LLvK is larger
ρl maxððpv  pÞ; 0Þ

m_ þ ¼ Ce α2 ð1  αÞ2 pffiffiffiffiffiffiffiffiffiffiffiffiffi ; if p < pv (10)
ρg 2π Rg T than CSAS 2σρΦk max k12 ∂∂xkj ∂∂xkj ; ω12 ∂∂ω ∂ω
xj ∂xj , the QSAS term will be active.

Here, α is the vapor fraction, ρl is the liquid density, ρv is the vapor

density, p is the local pressure, pv is the saturated vapor pressure, T is the 3. Numerical setup and validation
temperature and Rg ¼ 461.6 J/kgK is the gas constant. Cc is the rate co-
efficient of condensation when local pressure exceeds the saturated vapor The numerical simulation is conducted through the open source
pressure, and Ce is the rate coefficient of evaporation when local pressure software platform OpenFOAM (Jasak, 1996), which is based on finite
below the saturated vapor pressure. In the present work, Cc ¼ Ce ¼ 0.1 is volume method. The present volume of fraction (VoF) based compress-
taken according to the study by Saito et al. (2007). ible cavitating flow solver is developed by implementing the phase
The temperature-dependent saturation vapor pressure is given ac- change algorithm into the native pressure-based compressible two-phase
cording to the empirical formula (Schmidt and Grigull, 1980) as solver compressibleInterFoam of the open source package

 OpenFOAM-4.0. The phase change is implicitly treated as a source to
Tc 
pv ¼ pc exp  7:85823θ þ 1:83991θ1:5  11:7811θ3 þ 22:6705θ3:5 implement into the compressible phase volume fraction equation and
T compressible Pressure Poisson equation. The detailed description of the

 discretization of the governing equations above and the pressure-based
 15:9393θ4 þ 1:77516θ7:5
numerical method can be found in Miller et al. (2013).
(11)
3.1. Test case
with θ ¼ 1-T/Tc. The subscript “c” denotes temperature, pressure and
density at the critical point. For water pc ¼ 22.064 MPa, ρc ¼ 322 kg/m3,
In this paper, a two-dimensional foil of the NACA66 series by Leroux
Tc ¼ 647.14 K.
et al. (2004) is selected. The geometry coordinates of the foil can be
found in the reference (Leroux et al., 2004). It is characterized by a chord
2.4. Scale-adaptive-simulation (SAS) SST turbulence model length of c ¼ 0.150 m, the relative maximum thickness is τ ¼ 12% at 45%
from the leading edge, and the relative maximum chamber is 2% at 50%
In the present paper, the open source software OpenFOAM is used, from the leading edge. The experiment data shows that under the con-
and the compressible URANS equations are solved with the SAS SST dition of attack angle is 6 and the inlet flow velocity is U ¼ 5.33 m/s, at
turbulence model (Menter et al., 2003). The SAS SST turbulence model is the cavitation number is σ ¼ ðpout  pv Þ=ð0:5ρl U 2 Þ ¼ 1:25, lower
a hybrid RANS-LES model, where in the attached flow region like Strouhal number around 0.1 exists and the shock wave dynamics are
boundary layer, the SST model acts, showing RANS-like behaviors and in dominated in the unsteady cavitation instabilities.
the unsteady regions, the SAS is active, showing LES-like behaviors. It is In the present study, the large scale cavity cloud collapse emitted
not sensitive to the local length scale derived from the computation grid, shock wave induced newly attached sheet cavity being shorten phe-
so that can better capture unsteady cavitation instabilities characterized nomena and the pressure peaks are focused. In the simulation, the static
by the multi-scale dynamics. The von Karman length-scale, Lvk, based on pressure is adjusted to vary according to the cavitation number. The
the second velocity gradients is implemented into the ω equation as a water temperature is at 293.15 K.
sensor for detecting unsteadiness, which allows the SAS model to be
dynamically scale-adaptive to the resolved structures in the unsteady
3.2. Computational domain and boundary conditions
RANS calculation.
The turbulent kinetic energy (k) equation and the turbulent frequency
The 3D fluid domain is shown in Fig. 1 (a), which corresponds to the
(ω) equation are:
height of the experimental test section (Leroux et al., 2004). It should be
 
Dk ∂ ∂k noted that due to the large computational cost to use the real spanwise
¼ Pk  ε þ ðν þ σ k νT Þ (12)
Dt ∂xj ∂xj size of the experimental hydrofoil, the spanwise length is chosen as 0.3 c
in the present paper, the same domain treatment as Sagaut (2002) and Ji
  et al. (2015). The computational domain has an extent of about 2 c up-
Dω γ ∂ ∂ω σ ω2
¼ Pk  βω2 þ ðν þ σ ω νT Þ þ 2ð1  F1 Þ CDkω þ QSAS stream and 6 c downstream of the foil to simulate far-field boundary
Dt νT ∂xj ∂xj ω |ffl{zffl}
Addition term
conditions at the inlet and outlet, respectively. The 3D calculation mesh
(13) near the foil surface is shown in Fig. 1 (b), the total grid nodes are
approximately 2.5 million with 80 nodes uniformly distributed along the

426
C. Wang et al. Ocean Engineering 156 (2018) 424–434

Fig. 1. Computational domain and mesh details.

spanwise. Sufficient refinement near the foil surface was generated, with shown in Table 1. The relative error is within 3%, showing the good
the value of y þ at the first grid point within 1–2. The non-reflect accuracy of the present numerical results in predicting the typical cloud
boundary conditions are used for inlet and outlet to avoid the pressure cavitation dynamics characterized by low Strouhal number.
reflection. A no-slip boundary condition is imposed on the hydrofoil
surface, top and bottom tunnel boundaries, and symmetry conditions are 4. Results and discussion
imposed on the front and back boundaries. The Euler scheme for time
discretization and Gauss upwind scheme for spatial discretization are 4.1. Global multiphase structure associated with sheet/cloud cavitating
used in the present simulation. In order to reach a high accuracy in time, flows
the maximum Courant number (maxCo) is kept below 0.4 and thus the
time step 1  105 s is chosen. The total cavitation simulation time is Fig. 3 shows the comparisons between the experimentally observed
0.85 s, which includes three complete cavitation cycles to guarantee cavity behaviors (left, Leroux et al., 2004), the numerically predicted
simulation convergence. vapor volume fraction iso-surface with the absolute pressure contours on
symmetry plane and velocity distribution on foil surface (middle), and
the Q iso-surface (Q¼(ΩijΩij-SijSij)/2 ¼ 70000, by Hunt et al. (1988)
3.3. Numerical validation

Fig. 2 shows the time evolution of the cavity volume during three Table 1
cavitation cycles and the power spectral density (PSD) distribution based Comparisons of the measured (Exp.) and predicted (Num.) cavitation evolution
on Fast Fourier Transform (FFT) analysis of cavity volume. It can be frequency and the Strouhal number based on the foil chord length. σ ¼ 1.25,
found that the sheet/cloud cavitation presents quasi-periodic behaviors α ¼ 6 .
and the cavitation cycle is estimated to be about 3.52 Hz according to the Exp. (Leroux et al. (2004)) Num.
PSD analysis. The comparisons between the experimentally measured f 3.625 3.52 (2.9%)
and numerically predicted cavitation frequency (f) and the correspond- St ¼ fc/U∞ 0.102 0.099 (2.9%)
ing cavitation Strouhal number (St) based on the foil chord length are

Fig. 2. (a) The time evolution of cavity volume, and (b) PSD analysis cavity volume evolution.

427
C. Wang et al. Ocean Engineering 156 (2018) 424–434

Fig. 3. Comparisons of the experimentally observed cavitation pattern (left, Leroux et al. (2004)), numerically predicted vapor fraction isosurface (αv ¼ 0.15) and
vorticity distribution on symmetry plane (middle) and Q critier isosurface (Q ¼ 70000, right). Time interval between two consecutive images is 1/7 Tref.
(Q¼(ΩijΩij-SijSij)/2).

428
C. Wang et al. Ocean Engineering 156 (2018) 424–434

during one typical cavitation cycle. In the typical flow condition, a When the upstream propagating shock wave meets the newly attached
complex flow pattern, characterized by a lower Strouhal number, is ob- cavity sheet, it will cause the thin cavity sheet collapse. It is supposed the
tained. The leading edge (LE) and trailing edge (TE) of the foil are shock wave propagation and rebound is responsible for the reduction of
labeled, respectively. The flow direction is from right to left as indicated cavitation evolution Strouhal number. The detailed description of the
by the red arrow line in Fig. 3 (a). Meanwhile, the temporal-space evo- shock wave propagation process will be given following.
lution of the vapor volume fraction, absolute pressure and velocity dis- Fig. 6 shows the typical instances during the large-scale cloud cavity
tribution at the specific positions (as shown in Fig. 4) during one typical collapse process, namely U-shape cloud cavity formation, U-shape cloud
cavity evolution cycle are extracted and shown in Fig. 5. The positions of cavity head collapse, U-shape cloud cavity legs collapse and shock wave
the monitors are 1 mm above the foil surface and extended to 1.5 foil emission stages, respectively. It can be observed that the newly cavity
chord length. sheet grows to beyond x/c ¼ 0.3 before the cloud cavity collapse event
It is observed that the sheet/cloud cavitation has a distinct quasi- takes place. After the collapse event takes place, the emitted shock wave
periodic patterns. The unsteady cavity behaviors can be depicted as the propagates rapidly, when it arrives at the cavity sheet region, it will make
following four stages: (1) the growth of the attached cavity sheet in Fig. 3 the cavity sheet collapse. After that, the flows induced low pressure re-
(a), (b), (c), (d), and (i), (2) the development of re-entrant flow and gion grows towards the downstream from the foil leading edge to trailing
attached cavity sheet breakup in Fig. 3 (e), (3) attached cavity sheet edge. The absolute pressure fluctuations during the shock wave forma-
being rolled up and cloud cavity shed downstream in Fig. 3 (f) and (g), tion and propagation process at positions x/c ¼ 0.1–0.9 along the foil
(4) cloud cavity collapse and shock wave generation and propagation in surface at the mid-plane are plotted in Fig. 7 (a). The average void
Fig. 3 (h). During the first stage, the attached cavity sheet begins to form fraction is averaged using the void fraction at start point and end in-
at the foil leading edge and grows downstream towards the trailing edge, stances during shock wave propagation (αv ¼ αv,t1þαv,t2). Two shock
as shown in Fig. 3 (a), (b), (c) and (d). The Q-isosurface has the similar wave propagation events are observed, event I (cloud cavity collapse
shape as the attached cavity shape, indicating that the cavitation region is induced shock wave propagation) and event II (it is supposed to be
a high vorticity region. A vortex pair can also be found concentrate at the happened by the leading edge collapse induced shock wave propaga-
foil trailing edge. The velocity vector on the foil surface shows the ex- tion). According to the pressure signals, the average absolute pressure
istence of the reverse vortex there as indicated by the red curved arrows propagation speed along with the average void fraction and corre-
as shown in Fig. 3 (b). With the attached cavity development, the Q- sponding Wallis' sound speed (ρ1c2 ¼ ραcv2 þ 1 αv
ρ c2
) is shown in Fig. 7 (b). The
v v l l
isosurface becomes broken, illustrating small vortex structures within the
threshold value 25 kPa is used to calculate the cloud cavity collapse
attached cavitation region. As presented in Fig. 3(b), it can be found that
induced shock wave propagation speed. It can be found that lowest speed
the velocity magnitude distribution on the foil surface shows that the
of shock wave propagation in the cavitation region reaches 35.7 m/s
velocity within the attached cavity region is relatively lower than that
between 0.1 c to 0.2 c with the average void fraction around αv ¼ 0.13. It
outside of the attached cavity region, indicating the relatively stable
can be found that the vapor fraction has a significant effects on the
characteristics within the cavity. As shown by the velocity vectors on
pressure propagation speed. The cloud cavity collapse induced shock
region covered by the attached cavity, it can be concluded that cavitating
wave propagation speed in event I has the lower propagation speed than
flows are highly turbulent with many multi-scale vortical structures. As
that during the leading edge sheet cavity collapse induced shock wave
can be seen in Fig. 5 (c), the re-entrant flow appears at t/Tref ¼ 0.053, and
propagation in event II, where after the event I, the void fraction reduces
then develops with the attached cavity growth. During the second stage,
due to shock wave propagation, which can be observed by the average
when the attached cavity sheet reaches to a certain length, a relatively
void distribution in Fig. 7(b). During the leading edge sheet cavity
high speed re-entrant flow, as shown in Fig. 5 (c) with the blue arrow,
collapse induced shock wave propagation in event 2 just after the cloud
generates at the rear of the attached cavity sheet, and moves upstream.
cavity collapse induced newly cavity sheet collapse, the propagation has
When the re-entrant flow arrives at the leading edge, it will cut off the
the fast propagation speed due to collapse created lower vapor fraction. It
attached cavity sheet and cause the attached cavity sheet breakup. The Q-
should be noted that the Wallis's sound speed variation trend is in
isosurface is more broken, with small vortex structures generated. During
agreement with the void fraction variation, however, the Walls' sound
the re-entrant flow movement, higher and lower velocity streaks as
speed is smaller than the shock wave propagation speed based on the
shown in Fig. 3 (e) can be found. During the third stage, the broken cavity
pressure signals. The reasons for the discrepancy may be that firstly the
sheet is being rolled up and shed downstream, and eventually a U-shape
average shock wave speed calculation is based on the propagation pro-
cloud cavity forms. The U-shape cloud cavity legs consist of high vorticity
cess, and along the propagation path, the void fraction is variable, while
as shown in the Q-surface in Fig. 3 (g). During the fourth stage, when the
the Wallis' sound speed calculation is based on the average void fraction
cloud cavity being transported into high pressure region, it will collapse
on start and end points. Secondly, the Wallis's sound speed formula ig-
and the shock wave propagation and rebound phenomena along with
nores the non-equilibrium effects, heat transfer and mass transfer, which
high absolute pressure propagation takes place as indicated by the white
are important in the cavity cloud collapse process.
arrows in Fig. 5 (b). It can be found that with the shock wave rebound
process may be more than one time with the shock wave magnitude
attenuation. The shock wave will propagate both upstream and down- 4.2. The cavitation vortex dynamics analysis
stream, originating from the cloud cavity collapse position. The averaged
shock propagation speed from the cavity cloud collapse position to the Fig. 8 shows the three-dimensional isosurfaces of vorticity X, vorticity
foil leading edge is about 49.3 m/s based on the temporal-spatial diagram Y and vorticity Z at typical transient cavity evolution stages, the devel-
in Fig. 5(b), which reflects the sharp reduction in speed of sound in opment of attached cavity sheet, the movement of re-entrant flow, the
vapor/water mixture compared with that in pure water for 1450 m/s. shedding of cloud cavity and the cloud cavity collapse induced shock
wave emittion and propagation, respectively. The black isosurface in-
dicates the vorticity value of 500 and the white isosurface indicates the
vorticity value of 500. It can be found that the isosurface of vorticity Z has
the largest distribution region than that of vorticity X and vorticity Y,
indicating the vorticity Z dominates the cavitation vortex dynamics, and
the high-vorticity region in all three directions are all highly related to
the cavity regions. The trailing edge vortex is mainly composed of
vorticity Z as shown by the vorticity Z isosurface. During the attached
Fig. 4. NACA66 section and the monitor positions. cavity sheet growth stage, as shown in Fig. 8 (a), the vorticity X

429
C. Wang et al. Ocean Engineering 156 (2018) 424–434

Fig. 5. Numerically predicted temporal-space evolution of the predicted (a) water vapor fraction (αv), (b) absolute pressure and (c) reverse u-velocity at various
sections (the positions is reported in ordinate, from the foil leading edge x/c ¼ 0 to foil wake x/c ¼ 1.5) along the foil surface. The red vertical arrows indicate the first
shock wave propagation process. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6. Bird view of predicted isosurface of vapor fraction 0.15 and pressure distribution on the foil surface and side plane in the process of large scale cloud cav-
ity collapse.

distributes in the region of attached cavity and vorticity Y distributes D!

ω  ! rρ  rp
¼ ð!
ω ⋅rÞ  !
ω r⋅ V þ m 2 þ ðνm þ νt Þr2 !
ω (16)
streaks along the streamwise direction. During the re-entrant flow Dt ρm
development stage, as shown in Fig. 8 (b), the high vorticity region en-
In the transport equation, the first term on the right hand side (RHS)
larges, characterized by the larger vorticity X distribution region and the
is the vortex stretching term, representing the stretching and tilting of a
longer vorticity Y streaks. During the cloud cavity shedding stage, the
vortex by velocity gradients. The second term on the RHS is the vortex
cloud cavity mainly consists of vorticity X and vorticity Y as shown in
dilatation term, representing the volumetric expansion/contraction. This
Fig. 8 (c). During the cloud cavity collapse induced shock wave propa-
describes the role of fluid compressibility effects on the vorticity trans-
gation stage, vorticity Y can be found within the collapse region in the
port. The third term on the RHS is the baroclinic torque, representing the
foil wake as indicated by the red dashed circle. And the newly attached
role of misaligned pressure and density gradients. The last term on the
cavity sheet are composed of vorticity X and vorticity Z as indicated by
RHS is the viscosity diffusion term, indicating the rate at which the
the yellow dashed circles.
vorticity changes due to vorticity viscosity diffusion. It should be noted
From above, it can be obtained that the unsteady flow structures in
that the viscous diffusion term has a much smaller effects on the vorticity
the cavitating flow process is different, especially during the four cavi-
transportation than other three terms in high Reynolds number flows
tation development stages. In order to investigate the differences of
(Huang et al., 2013b; Ji et al., 2014).
cavitation-vortex interaction in the cavitation evolution, especially in the
The predicted contours of vortex stretching term, vortex dilatation
process of the re-entrant flow development and the shock wave propa-
term, and baroclinic torque term in vorticity Z transport equation, which
gation, the vorticity transportation equation in the compressible turbu-
dominates the vortex dynamics, at the mid-plane of the hydrofoil, as well
lent flow region is employed as follows:
as the ratio of baroclinic term and vortex stretching term and the vortex

430
C. Wang et al. Ocean Engineering 156 (2018) 424–434

Fig. 7. (a) Numerical pressure signals during the cloud cavity collapsing for stage 4 on foil surface at x/c ¼ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9, and (b)
propagation speed of shock wave based on pressure signals along with the average void fraction and corresponding Wallis' sound speed shown by gray lines (The
monitor locations are shown in Fig. 4.).

Fig. 8. Isosurfaces of (a) vorticity X, (b) vorticity Y and (c) vorticity Z at different cavitation evolution stages. Also shown is the (x, y) axis system used. (Black
isosurface for ω ¼ 500 and white isosurface for ω ¼ 500, foil surface is colored by the velocity magnitude).

431
C. Wang et al. Ocean Engineering 156 (2018) 424–434

Fig. 9. Budget of vortex transport equation at the mid-plane of hydrofoil for vortex stretching term (left), vortex dilatation term (middle), and baroclinic torque term
(right), and the corresponding ratio of baroclinic torque to the vortex stretching and vortex dilatation.

432
C. Wang et al. Ocean Engineering 156 (2018) 424–434

dilatation term are presented in Fig. 9 at four typical cavitating shedding baroclinic term dominate. The ratio of baroclinic term and the
instances, namely the attached cavity growth, the re-entrant flow vortex stretching term/vortex dilatation term further indicates,
movement, cloud cavity being shedding and cloud cavity collapse the dominant role of baroclinic term in the cavity growth, devel-
induced shock wave propagation process, respectively, to study the opment and shedding process. During the cavity collapsed
cavitation-vortex interaction more detailed. induced shock wave propagation process, the budget of vorticity
From the contours of the ratio of baroclinic term and the vortex transport as well as the ratio of baroclinic torque term to the
stretching term and the vortex dilatation term respectively, it can be vortex stretching/vortex dilatation shows that the vortex dilata-
clearly seen that the baroclinic term mainly dominates the vorticity tion term plays the major role in the vortex transport in the
transportation process. In the process of attached cavity growth in Fig. 9 attached cavity sheet region.
(a), the vortex stretching term and vortex baroclinic term distribution is
stronger than that the vortex dilatation term. The ratio of the baroclinic Regarding future work, further research on unsteady cavitation vor-
term and vortex stretching term shows the dominant influence of baro- tex structures will be conducted from the Lagrangian viewpoint. The
clinic term on the vorticity transport process. In the re-entrant flow Lagrangian Coherent Structures (LCS) method has the advantage that the
process, as shown in Fig. 9 (b), it can be observed that the vortex the LCS method could define transient vortex structure boundaries without
baroclinic torque term plays a main role in the vortex transport process relying on a threshold and could reveal more vortical structures than
than the vortex stretching term and the vortex dilatation term, as well as classical Eulerian fields.
that indicated by the ration of baroclinic term and the vortex stretching
term and the vortex dilatation term. In the cavity cloud shedding process Acknowledgement
in Fig. 9 (c), the vortex dilatation term and the baroclinic term are
stronger than the vortex stretching term, and from the ratio distribution, The authors gratefully acknowledge support by the National Post-
it can be observed that the baroclinic term dominates. In the cavity doctoral Program for Innovative Talents (Grant No: BX201700126), the
collapsed induced shock wave propagation process in Fig. 9 (d), the National Foundational of China (NSFC, Grant Nos: 51239005 and
vortex dilatation term plays the major role in the vortex transport, 51679005), National Natural Science Foundation of Beijing (Grant No:
indicated by both the budget of vorticity transport equation distribution 3172029), and the Open Foundation of State Key Laboratory of Ocean
and the budget ratio. It may be the reason that the re-entrant flow is Engineering (Shanghai Jiao Tong University, China) (Grant No: 1611).
induced by the adverse pressure gradient, which will significantly in-
fluence the baroclinic torque term, while the shock wave is highly related References
with the compressibility effects, and thus the vortex dilatation term is
mainly influenced by the shock wave propagation. Arndt, R.E.A., Song, C.C.S., Kjeldsen, M., et al., 2001. Instability of partial cavitation: a
numerical/experimental approach. In: Symposium on Naval Hydrodynamics.
Avellan, F., Farhat, M., 1989. Shock wave pressure generated by cavitation vortex
5. Discussion and conclusions collapse. In: ASME Int. Symp. On Cavitation Noise and Erosion in Fluid Systems, USA.
Bensow, R.E., Bark, G., 2010. Implicit LES predictions of the cavitating flow on a
propeller. J. Fluid Eng. 132, 041302.
In this paper, the unsteady cavitating flow around a NACA66 hy- Brennen, C.E., 1995. Cavitation and Bubble Dynamics. Oxford University Press.
drofoil fixed at the angle of attack AoA ¼ 6 at σ ¼ 1.25 has been Budich, B., Schmidt, S.J., Adams, N.A., 2015. Numerical investigation of a cavitating
numerically simulated, and the unsteady cavitating flow structures and model propeller including compressible shock wave dynamics. In: Fourth
International Symposium on Marine Propulsors.
the cavitation vortex dynamics, especially in the process of the re-entrant Chen, G., Wang, G., Hu, C., Huang, B., et al., 2015. Observations and measurements on
flow movement and the cloud cavity collapse induced shock wave for- unsteady cavitating flows using a simultaneous sampling approach. Exp. Fluid 56,
mation and propagation, have been analyzed. The numerical simulation 1–11.
Chen, T., Huang, B., Wang, G., Zhao, X., 2016. Numerical study of cavitating flows in a
is performed using compressible cavitation solver, closured by the Tait wide range of water temperatures with special emphasis on two typical cavitation
equation of state for water and ideal gas equation of state for vapor to dynamics. Int. J. Heat Mass Tran. 101, 886–900.
consider both the water, vapor and water/vapor mixture compressibility, Chen, Y., Chen, X., Li, J., Gong, Z., Lu, C., 2017. Large eddy simulation and investigation
on the flow structure of the cascading cavitating shedding regime around 3D twisted
based on the OpenFOAM platform. Several flow features involved tran-
hydrofoil. Ocean. Eng. 129, 1–19.
sient cavitation behaviors, accompanied with the re-entrant flow and Core, R.H., 1948. Underwater Explosion. Princeton Univ. Press, Princeton, NJ.
shock wave formation and propagation are addressed. Coutier-Degosha, O., Fortes-Patella, R., Rebound, J.L., 2003. Evaluation of the turbulence
model influence on the numerical simulations of unsteady cavitation. J. Fluid Eng.
125, 38–45.
(1) The cavity evolution presents quasi-periodic behaviors, and the Ducoin, A., Huang, B., Young, Y.L., 2012. Numerical modeling of unsteady cavitating
alternative re-entrant flow movement and shock wave dynamics flows around a stationary hydrofoil. Int. J. Rotating Mach. 2012.
play an important role in cavitation vortex dynamics. The un- Egerer, C.P., Schmidt, S.J., Hickel, S., Adams, N.A., 2016. Efficient implicit LES method
for the simulation of turbulent cavitating flows. J. Comput. Phys. 316, 453–469.
steady cavitation behaviors can be divided into four stages: (a) the Ganesh, H., M€akiharju, S.A., Ceccio, S.L., 2016. Bubbly shock propagation as a
growth of the attached cavity sheet, (b) the development of re- mechanism for sheet-to-cloud transition of partial cavities. J. Fluid Mech. 802, 37–78.
entrant flow and attached cavity sheet breakup, (c) attached Gnanaskandan, A., Mahesh, K., 2016. Large eddy simulation of the transition from sheet
to cloud cavitation over a wedge. Int. J. Multiphas. Flow 83, 86–102.
cavity sheet being rolled up and cloud cavity shed downstream Hanson, I., Kedrinskii, V.K., Mørch, K.A., 1981. On the dynamics of cavity clusters.
and (d) cloud cavity collapse and shock wave generation and J. Appl. Phys. 15, 1725–1734.
propagation. Among which, the cloud cavity collapse process can Haren, W., Heged} us, F., Kai, Y., Koch, S., Meyerer, B., Neu, W., Teubner, U., 2016. Shock
wave emission during the collapse of cavitation bubbles. Shock Waves 26, 385–394.
be further depicted into three stages: (a) the formation of the U- Huang, B., Ducoin, A., Young, Y.L., 2013a. Physical and numerical investigation of
shape cloud cavity; (b) the collapse of U-shape cavity head, and (c) cavitating flows around a pitching hydrofoil. Phys. Fluids 25, 102109.
the collapse of U-shape cavity legs. The collapse of the U-shape Huang, B., Young, Y.L., Wang, G., Shyy, W., 2013b. Combined experimental and
computational investigation of unsteady structure of sheet/cloud cavitation. J. Fluid
cavity legs, is the main source of shock wave emittion.
Eng. 135, 071301.
(2) The cavitation vortex dynamics show strong coherent relationship Huang, B., Zhao, Y., Wang, G., 2014. Large eddy simulation of turbulent vortex-cavitation
with the transient cavity behaviors. The vorticity Z dominates the interactions in transient sheet/cloud cavitating flows. Comput. Fluids 93, 113–124.
cavitation vortex dynamics. During the attached cavity growth Hunt, J.C.R., Wray, A.A., Moin, P., 1988. Eddies, stream, and convergence zones in
turbulent flows. In: Center for Turbulence Research Report, CTR-S 88. Center for
stage, the vortex stretching term and vortex baroclinic term dis- Turbulence Research, Stanford, CA.
tribution is stronger than that the vortex dilatation term. During Jasak, H., 1996. Error Analysis and Estimation for the Finite Volume Method with
the re-entrant flow process, the vortex the baroclinic torque term Applications to Fluid Flows. Ph.D. Thesis. University of London.

plays a main role in the vortex transport process. During the cavity
cloud shedding process, the vortex dilatation term and the

433
C. Wang et al. Ocean Engineering 156 (2018) 424–434

Ji, B., Luo, X.W., Arndt, R.E.A., Wu, Y.L., 2014. Numerical simulation of three Reisman, G.E., 1996. Dynamics, Acoustics and Control of Cloud Cavitation on Hydrofoils.
dimensional cavitation shedding dynamics with special emphasis on cavitation- Ph.D. thesis. California Institute of Technology.
vortex interaction. Ocean. Eng. 87, 64–77. Reisman, G.E., Wang, Y., et al., 1998. Observations of shock waves in cloud cavitation.
Ji, B., Luo, X.W., Arndt, R.E.A., Peng, X., Wu, Y., 2015. Large eddy simulation and J. Fluid Mech. 335, 255–283.
theoretical investigations of the transient cavitating vortical flow structure around a Reisman, G.E., McKenney, E.A., Brennen, C.E., 1994. Cloud cavitation on an oscillating
NACA66 hydrofoil. Int. J. Multiphas. Flow 68, 121–134. hydrofoil. In: Proc. 20th INR Symp. On Naval Hydrodynamics, pp. 78–89.
Joseph, D.D., 1995. Cavitation in a flowing liquid. Phys. Rev. E 51, 1649–1650. Reisman, G.E., Brennen, C.E., 1996. Pressure pulses generated by cloud cavitation. In:
Knapp, R.T., Daily, J.W., Hammitt, F.G., 1970. Cavitation. McGrwa Hill. Fluids Eng. Division Conference, vol. 236, pp. 319–328.
Kunz, R.F., Lindau, J.W., Billet, M.L., et al., 2001. Multiphase CFD Modeling of Sagaut, P., 2002. Large Eddy Simulation for Incompressible Flows. Springer.
Development and Supercavitating Flows, Multiphase Cfd Modeling of Development & Saito, Y., Takami, R., Nakamori, I., et al., 2007. Numerical analysis of unsteady behavior
Supercavitating Flows. of cloud cavitation around a NACA0015 foil. Comput. Mech. 40, 85–96.
Le, Q., Franc, J.P., Michel, J.M., 1994. Parial cavities: pressure pulse distribution around Schmidt, E., Grigull, U., 1980. Properties of Water and Stream in SI-Units, fourth ed.
cavity closure. ASME J. Fluid Eng. 115, 249–254. Springer, Oldenbourg, Berlin.
Leroux, J.B., Astolfi, J.A., Billard, J.Y., 2004. An experimental study of unsteady partial Schnerr, G.H., Sezal, I.H., Schmidt, S.J., 2008. Numerical investigation of three-
cavitation. J. Fluid Eng. 126, 94–101. dimensional cloud cavitation with special emphasis on collapse induced shock wave
Leroux, J.B., Coutier-Delgosha, O., Astolfi, J.A., 2005. A joint experimental and numerical dynamics. Phys. Fluids 20, 040703.
study of mechanisms associated to instability of partial cavitation on two- Sezal, I.H., Schmidt, S.J., Schnerr, G.H., Thalhamer, M., Fӧrster, M., 2009. Shock and
dimensional hydrofoil. Phys. Fluids 17, 052101. wave dynamics in cavitating compressible liquid flows in injection nozzles. Shock
Li, X., Wang, G., Zhang, M., et al., 2008. Structure of supercavitating multiphase flows. Waves 19, 49–58.
Int. J. Therm. Sci. 47, 1263–1275. Soyama, H., Kato, H., Oba, R., 1992. Cavitation observations of severely erosive vortex
Long, X., Cheng, H., Ji, B., Arndt, R.E.A., 2017a. Numerical investigation of attached cavitation arising in a centrifugal pump. In: Proceedings of the Third IMechE
cavitation shedding dynamics around the clark-Y hydrofoil with the FBDCM and an International Conference on Cavitation.
integral method. Ocean. Eng. 137, 247–261. Venkateswaran, S., Lindau, J.W., Kunz, R.F., et al., 2002. Computational of multiphase
Long, Y., Long, X., Ji, B., Huai, W., Qian, Z., 2017b. Verification and validation of URANS mixture flows with compressibility effects. J. Comput. Phys. 180, 54–77.
simulations of the turbulent cavitating flow around the hydrofoil. J. Hydrodyn. 29 Wallis, G., 1967. One-dimensional Two-phase Flow. McGraw-Hill, New York.
(4), 610–620. Wang, C., Huang, B., Wang, G., et al., 2017a. Unsteady pressure fluctuation characteristics
McKenney, E.A., Brennen, C.E., 1994. On the dynamics and acoustics of cloud cavitation in the process of breakup and shedding of sheet/cloud cavitation. Int. J. Heat Mass
on an oscillating hydrofoil. In: Proc. ASME Symp. On Cavitation and Gas-liquid Flows Tran. 114, 769–785.
in Fluid Machniery and Devices, vol. 190, pp. 195–202. Wang, G., Ostoja-Stazewski, M., 2007. Large eddy simulation of a sheet/cloud cavitation
Menter, F.R., Kuntz, M., Bender, R., 2003. A scale-adaptive simulation model for on a NACA0015 hydrofoil. Appl. Math. Model. 31, 417–447.
turbulent flow predictions. In: AIAA 2003-0767, 41st Aerosppace Science Meeting Wang, G., Senocak, I., Shyy, W., Ikohagi, T., Cao, S., 2001. Dynamics of attached
and Exhibit, Teno, Nevada. turbulent cavitating flows. Prog. Aero. Sci. 37, 551–581.
Miller, S., Jasak, H., Boger, D., Paterson, F., Nedungadi, A., 2013. A pressure-based, Wang, G., Wu, Q., Huang, B., 2017b. Dynamics of cavitation-structure interaction. Acta
compressible, two-phase flow finite volume method for underwater explosions. Mech. Sin. 33, 685–708.
Comput. Fluids 87, 132–143. Wang, Y., Abdel-Maksoud, M., Wang, K., Song, B., 2016. Prediction of tip vortex
Mørch, K.A., 1980. On the collapse of cavity cluster in flow cavitation. In: Proc. First Int. cavitation inception with low-order panel method. Ocean. Eng. 125, 124–133.
Conf. On Cavitation and Inhomogenieties in Underwater Acoustics, Springer Sries in Wu, J., Ganesh, H., Ceccio, S.L., 2016a. Acoustic emission due to partial cavity shedding
Electrophysics, vol. 4, pp. 95–100. on a hydrofoil. J. Acoust. Soc. Am. 140, 3330–3331.
Mørch, K.A., 1981. Cavity cluster dynamics and cavitation erosion. In: Proc. ASME Wu, Q., Huang, B., Wang, G., Gao, Y., 2015. Experimental and numerical investigation of
Cavitation and Polyphase Flow Forum, pp. 1–10. hydroelastic response of a flexible hydrofoil in cavitating flow. Int. J. Multiphas. Flow
National Institute of Standards and Technology (NIST), U.S. Department of Commence. 74, 19–23.
http://webbook.nist.gov/chemistry/. Wu, Q., Huang, B., Wang, G., 2016b. Lagrangian-based investigation of the transient flow
Noordij, L., Wijngaarden, L., 1974. Relaxation effects, caused by relative motion, on shock structures around a pitching hydrofoil. Acta Mech. Sin. 32, 64–74.
waves in gas-bubble/liquid mixtures. J. Fluid Mech. 66, 115–143. Wu, Q., Huang, B., Wang, G., Cao, S., 2018. The transient characteristics of cloud
Reboud, J.L., Stutz, B., Coutier-Degosha, O., 2008. Two-phase flow structure of cavitation cavitating flow over a flexible hydrofoil. Int. J. Multiphas. Flow 99, 162–173.
experiment and modeling of unsteady effects. In: 3rd International Symposium on Yu, C., Wang, Y., Huang, C., Wu, X., Du, T., 2017. Large eddy simulation of unsteady
Cavitation (CAV1998), Grenoble, France. cavitating flow around a highly skewed propeller in nonuniform wake. J. Fluid Eng.
139, 041302 .

434