Cavitation vortex dynamics of unsteady sheet/cloud cavitating flows with shock wave using different vortex identification methods *

Chang-chang Wang, Ying Liu, Jie Chen, Fu-yi Zhang, Biao Huang, Guo-yu Wang

School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China

Abstract: Cavitation is a complex multiphase flow phenomenon with an abrupt transient phase change between the liquid and the vapor, including multiscale vortical motions. The transient cavitation dynamics is closely associated with the evolution of the cavitation vortex structures. The present paper investigates the cavitation vortex dynamics using different vortex identification methods, including the vorticity method, the Q criterion method, the Omega method （Ω ）, the method and the Rortex method. The Q criterion is an eigenvalue-based criterion, and in the Ω method, the parameter is normalized, is independent of the threshold value and in most conditions Ω= 0.52. The Rortex method is based on an eigenvector-based criterion. Numerical simulations are conducted using the implemented compressible cavitation solver in the open source software OpenFOAM for the sheet/cloud cavitating flows around a NACA66 (mod) hydrofoil fixed at α=6°, σ=1.25 and Re =7.96 ×1 05. The flow is characterized by the alternate interactions of the re-entrant flow and the collapse induced shock wave. Results include the vapor structures and the vortex dynamics in the unsteady sheet/cloud cavitating flows, with emphasis on the vortex structures in the cavitation region, the cavity interface, the cavity closure, the cavity wakes, and the foil wakes with the shedding cavity. The comparisons of the various methods, including that the vorticity method, the Q criterion method, the Ω method, the method and the Rortex method, show the performances of different methods in identifying the cavitation vortex structures. Generally,during the attached cavity growth stage, the Q criteria can well predict the vortex structures in the cavitation region and at the foil trailing edge in the pure liquid region, while with the Ω method and the Rortex method, the vortex structures outside the attached cavity and on the foil pressure side can also be predicted. The method can well predict the vortex structures in the cavity closure region. During the re-entrant jet development stage, the vortex structures in the re-entrant jet region is weak. During the cavity cloud shedding stage, the vortex dynamics at the foil leading edge covered by newly grown cavity sheet is different from that during the attached cavity sheet growth stage. During the shock wave formation and propagation stage, strong vortex structures with both the size and the strength are observed owing to the cavity cloud shedding and collapse behavior. The influence of the small parameter ε in the Ω method on the cavitation vortex identification is discussed.

Key words: Cavitation vortex dynamics, vortex identification, sheet/cloud cavitating flow, re-entrant jet, shock wave, OpenFOAM

Introduction

The cloud cavitation is characterized by the periodic attached cavity growth and the large scale cavity cloud shedding, involving the multi-scale vortical motions from microbubble clusters in the micron dimension to the large scale cavity structure in the centimeter dimension[1-4]with different strength.The cavitating microvortices (e.g., the hairpin vortex)at the trailing edge of the attached cavity sheet and the large scale shedding vortical cavity cloud are known to be highly erosive[5-6]. Experiments[7-8]identified the re-entrant jet dominant mechanism and the shock wave dominant mechanism as the cause of the onset of the cloud cavitation, but with a limited understanding of the shock wave mechanism. In the shock wave dominant mechanism, the large scale cavity cloud collapse induced shock wave formation and propagation may serve as the important vorticity generation mechanism in the flow field and will greatly alter the cavitation dynamics[9]. Furthermore,with the violent mass transfer between the water and the vapor and thus the cavity structure changes in the process of the cavity cloud collapse, and its coherent interactions with the vortical structures, all make the cavitation vortex dynamics very complex. The cavitation vortex structure has attracted great attention in both scientific research and engineering fields[2,4,10-11].The understanding of the cavitation vortex dynamics will shed light on the cavitation instabilities,especially the shock wave mechanism.

In the unsteady cavitating flows, the cavitation instabilities are closely related with the formation of multiple scale cavitation vortical structures. It is indicated that the cavitation inception process is closely related with the local turbulence structures[12].On the other hand, the developed cavitation is also an important mechanism of the turbulence generation.Gopalan and Katz[3]applied the particle image velocimetry (PIV) and the high-speed photography to measure the flow structures in the closure region and the places downstream the sheet cavity. Results show that the vapor collapse in the closure region is the primary mechanism of the vorticity generation, and will significantly increase the turbulence, the momentum and the displacement thickness. The Reynolds stress could be increased by 25%-40% due to the cavitation. Iyer and Ceccio[13], Laberteaux and Ceccio[14]further pointed out that the baroclnic torques were responsible for the production of the vorticity during the vapor cloud collapse due to the alignment of the density gradient and the pressure gradient. However, owing to the high void fraction within the attached cavity sheet, the velocity measurement in the cavitation region is difficult, with only a limited understanding of the inner structures within the attached cavity. To further investigate the complex cavitation-vortex interactions, the accurate identification of the vortex structures is fundamental for the investigation of the cavitation vortex dynamics and various vortex identification methods, including the Eulerian local vortex identification methods[11-12,15-20](e.g., the vorticity threshold method, the Q- criterion method, the2λ method) and the Lagrange methods(e.g., the LCS)[21-22]are utilized for the cavitating flows. Ji et al.[10]conducted numerical simulations to investigate the cavitation shedding dynamics around a hydrofoil using the large eddy simulation method. The vortex structures are identified using the Q- criterion method withFurthermore, the cavitation vortex structures around a 3-D twisted hydrofoil are investigated in detail by using the Euler-Lagrangian method[19]. Huang et al.[23]numerically investigated the turbulent vortex-cavitation interactions in the unsteady sheet/cloud cavitating flows by using the large eddy simulation. The Q contour with the value of 0-105is used for the vortex structure visualization.The detailed vorticity generation and evolution mechanism is studied based on the budget analysis of the vorticity transport equation[2,24]. They concluded that the periodic formation, breakup, shedding, and collapse of the sheet/cloud cavities, along with the associated baroclinic and viscoclinic torques, are the important mechanisms for the vorticity generation and modification. Zhao et al.[25]conducted the Lagrangian analysis of the vortex dynamics by using the Lagrangian Coherent Structures (LCS) and the particle trajectory in the time-dependent cavitating flow around a clark-Y hydrofoil. Results show that during the re-entrant jet development, complex interactions exist between the LE-LCS and the TE-LCS, where the vortical structures inside the attached cavity are strongly unsteady. However, the research focus is mostly on the entire vortical region characteristics, with little attention paid to the mechanism of the well-defined cavitation vortical structures and their evolution (e.g., the U-shape shedding cavity cloud, the hairpin vortex at the trailing edge of the cavity sheet initiated in the attached cavity). To address the cavitation vortex dynamics and its relationship with the cavitation instabilities, especially in the process of the cavity cloud collapse induced shock wave under the shock wave mechanism, the development and the improvement of the cavitation vortex identification methods are essential for better understanding the cavitation instabilities.

The Q- criterion method is widely used in the cavitating flows in the Euler point of view for the vortex dynamics, however the choice of the critical Q value to identify certain cavitation vortex structures[26]is usually arbitrary based on experience.Traditionally, the vorticity method is used to describe vortices according to the Helmholtz's theorems.However, the vorticity-based definitions of the vortex cannot distinguish adequately the vortex from the shear flow, for example, in near wall regions without vortical motions but with vorticity[27-28]. Moreover, the point vortex has vortex structures but without vorticity.Consequently, the vortex identification methods were extensively explored to overcome these difficulties,and numerous vortex identification methods,including the Q criterion method[29], the Δ criterion method[30], thecriterion method[31], themethod[32], the Ω criterion method[33], etc.,were proposed. Zhang et al.[34]conducted a review of the vortex identification methods in applications,including the boundary layer transition, the wake flows, the atmosphere boundary layer, and the reversible pump turbine. They also proposed a method of determining the parameter in the Ω method and an interpretation of its physical meaning. Recently,Liu et al.[35-37]introduced a vector called the Rortex as a new vortex identification method, with the advantage of accurate prediction of both the local rotational direction and strength. However, these methods are based on the single phase flows. It is necessary to systematically validate the applications of these methods in the cavitating flows of liquid/vapor mixtures with the abrupt phase change,and to develop the identification criterion for the cavitation vortex structures.

The present work conducts numerical simulations of the time-dependent cloud cavitating flows around a NACA66 (mod) hydrofoil, by using the implemented compressible cavitation solver[38]based on the open source software OpenFOAM. The compressibility of both the liquid and the water is considered, to enable the cavity collapse induced shock wave to be captured.The cavitation instabilities are dominated by the alternate interactions of the re-entrant jet and the shock wave[38-39]. The objectives of the present paper are: (1) for the applications of different vortex identification methods including the vorticity method,the Q criterion method, the Ω method, the2λ method and the Rortex method in the cloud cavitating flows, (2) to investigate the cavitation vortex dynamics, with emphasis on the vortex structures in the process of the re-entrant jet development and the shock wave propagation, and (3) to discuss the influence of the small parameter ε on the performance of the Ω method.

1. Mathematical and physical model

1.1 Compressible Navier-Stokes equations

Based on the homogenous assumption that the two-phase flow is treated as a single mixture of the liquid and the vapor, the compressible Navier-Stokes equations including the mass equation, the momentum equation, and the energy equation, along with the void fraction transportation, are given below in the Cartesian coordinates.

where ρ, U , e and α are the density, the velocity vector, the internal energy and the void fraction,respectively, p , K and T are the pressure, the kinematic energy, and the temperature, κ=-??, αeffis the effective thermal diffusivity,is the relative velocity between the liquid phase and the vapor phase, as modeled according to Ref.[40],+m , m-are the mass transfer rates modeled by the cavitation model. More detailed description can be found in Ref. [38], the subscripts m, l and v refer to the mixture, liquid and vapor phases. The compressibility of the water and the vapor is considered by employing the equations of state, with the Tait equation of state for the water and the ideal gas equation of state for the vapor, respectively.

The pure liquid phase, the water, is described by the Tait EOS

The perfect gas law is used to describe the pure vapor phase

where the subscript v denotes the vapor-phase value and =461.6 J/(kg K)? is the gas constant.

The sonic speed in the water/vapor mixtures as a function of the void fraction is shown in Fig. 1. It can be found that the sonic speed predicted by the present model agrees well with the experiment data. Under the cavitation condition where the pressure is 2 338.16 Pa,the sonic speed is reduced to a value even below 3 m/s as shown by the arrow in Fig. 1, indicating a strong compressibility.

The transport proprieties of the mixture are given by the linear combination of the liquid and vapor phases as follows:

where ρ, μ are the density, and the dynamic viscosity, respectively.

Fig. 1 Sonic speed in water/vapor mixtures[38]

1.2 Mass transfer model and SST-SAS hybrid turbulence model

The transport-based Saito cavitation model[41]is used, where the temperature variance caused by the heat transfer along with the phase change is considered between the liquid and the vapor. The mass transfer rates measured by the source terms in Eq. (4) are expressed as follows:

wherecC is the model coefficient for the condensation from the vapor to the liquid andeC is the model coefficient for the evaporation from the liquid to the vapor. In the present study,is taken according to Ref. [41].

The RANS/LES hybrid turbulence model SAS is adopted, where the transition between the RANS simulation and the LES simulation is based on the local turbulence length other than the meshscale.Theshear stress transport (SST) turbulence is used for the RANS model and the WALE model is used for the LES. The second velocity gradient based von Karman length-scale, Lvk, is implemented into the ω equation as a sensor for detecting the flow unsteadiness.The SST-SAS[42]is dynamically scale-adaptive to the resolved structures in the flow field in order to capture more flow details, e.g., the vortex structures. The SAS term in the ω equation is expressed as follows

Numerical results are for the sheet/cloud cavitating flows with the re-entrant jet/shock wave alternate interaction around a NACA66 (mod) hydrofoil fixed at α=6° ,=5.33 m/s, which was experimentally investigated by Leroux et al.[39].According to Leroux et al.[39], the partial cavitating flows are characterized by the process of the alternate interactions of the re-entrant jet and the shock wave,caused by the cavity cloud collapse. The numerical simulation is conducted by using the implemented compressible cavitating flow solver[38]in the OpenFOAM software. The compressible cavitating flow solver is based on the finite volume method. In the simulation, the density variance is considered by employing the equations of state for both the liquid and the vapor, and the temperature variance is considered by solving the energy equation. The Saito cavitation model[41]and the SST-SAS turbulence model[42]are applied. More details about the numerical method, the computation domain and the boundary conditions can be found in Ref. [38]. In the previous study, the cavity cloud collapse induced shock wave is well captured by the solver[38].

The non-dimensional parameters used are defined as follows:

where c is the foil chord length, ν is the kinematic viscosity and poutis the outlet pressure.

2. Vortex identification methods

The characteristic equation of the velocity gradient tensor is given as follows:

where P, Q and R are the three invariants of the velocity gradient tensor. The velocity gradient tensor?u can be decomposed into two parts as the symmetric part (the deformation rate tensor) and the anti-symmetric part (the rotation rate tensor). In what follows, the definitions of different vortex identification methods based on the velocity gradient tensor are given.

2.1 Vorticity method

The vorticity has a clear mathematical definition as the anti-symmetry tensor of the velocity gradient tensor, e.g., the part b in Eq. (16).

The vortex is said to be identified in the region when the magnitude of the vorticity larger than a certain threshold,.

2.2 Q criterion method

The Q criterion is the second invariant of the velocity gradient tensor in Eq. (16), and is defined as[29]:

2.3 Ω method

The definition of Ω is as[33]:

where a, b are the squares of the Frobenius norm of A, B in Eq. (12). The physical meaning of Ω is how the vorticity overtakes the deformation.Assuming that both a, b are positive, the value of Ω is in the range of 0-1. When Ω=0, the flow is pure shear; when Ω=1, the flow is pure rigid body rotation.

2.4 2λ method

Based on the criterion that a pressure minimum should exist across a vortex core, the definition of2λ is as follows[31]

where, if one defines1λ ,2λ and3λ as the three real eigenvalues ofwith the order＞the second largest eigenvalue2λ ofis utilized as the vortex identification criterion. The region where λ2＜0 indicates the vortex region.

2.5 Rortex method

The Rortex vector[35-37]is defined as follows:

where U , V are in the XYZ- frame, and there exists a corresponding transformation between ?V ,?v[35-37]

where Q is a rotation matrix.

3. Vortex identification methods

3.1 Global multiphase structures of unsteady sheet/cloud cavitating flows

Figure 2 shows the time evolution of the normalized cavity volume,maxand the normalized cavity volume rates,cycle at=5.33 m/s. Also shown in Fig.2 is the bird view of the numerically predicted unsteady cavitation behaviors with the absolute pressure on both the foil surface and the side symmetry plane in eight typical instances. It can be found that the cavity volume shows unsteady characteristics in accordance with the transient cavity behaviors. The isosurface valueis used to draw the multiphase structures. The flow direction is from left to right as marked by the red arrow in picture 1. The comparisons between the computational results and the experimental data[39]are presented in Fig. 3,which shows a good agreement. The instances in Fig.)duringoneandaquartercavity, =1.25 2 are indicated by open circles. It should be noted that the pressure evolution on the foil suction side is strongly influenced by the cavity behaviors, with a reasonable deviation.

Fig. 2 (Color online) Time evolution of the normalized cavity volume and cavity volume with eight typical cavity structures marked by open symbols (isosurface, = 0.15 v α )during approximately one cavity cycle

Fig. 3 (Color online) Time evolution of pressure signals at x/c = 0.3c , 0.5c and 0.7c during one typical cavity shedding cycle

For the time-dependent partial cavitating flows around a NACA66 hydrofoil (mod) at =5.33 m/s

(1) Phase 1 attached cavity growth (picture 1,picture 2, picture 7 and picture 8): The cavity sheet is initiated at the foil leading edge and grows downstream, attaching on the foil surface all the time.The leading edge of the sheet is smoothly extended from the vapor fraction isosurface, while the trailing edge of the sheet is unsteady with a wavy cavity closure. The foil covered by the cavitation is acted by a low pressure. During this stage, the cavity volume gradually increases, indicating that the total mass transfer is from the liquid to the vapor. Because the specific volume of the vapor is much higher than that of the liquid, the total volume of the cavitating fluid increases abruptly in this process.

(2) Phase 2 re-entrant flow development (picture 3): When the cavity sheet length grows to a certain length about l / c=0.5, the cavity growth slows down and the unsteadiness of the cavity closure increases,when a re-entrant jet forms and develops. As shown in the cavity volume curve, it can be found that during this stage, the cavity volume reaches its maximum value and remains unchanged, indicating the relative balance between the evaporation rate and the condensation rate during the re-entrant movement beneath the cavity.

(3) Phase 3 cavity cloud shedding (picture 4):When the re-entrant jet reaches the foil leading edge,the interaction between the re-entrant jet front and the cavity interface will cause the cavity breakup. Under the action of the shear force of the main flow, the residual cavity sheet will roll up and shed downstream as a large scale cavity cloud as shown in picture 4. As illustrated in the cavity volume, the cavity volume is in a constant period, and a detailed study indicates that in this period, the shedding cavity cloud is just transported above the foil, showing the balance between the evaporation rate in the attached cavity sheet region and the condensation rate in the shedding cavity cloud.

(4) Phase 4 cavity cloud collapse and shock wave emission and propagation (picture 5, picture 6): When the cavity cloud is transported downstream into the high pressure region, a large pressure difference will cause the cavity cloud collapse. With the collapse of the cavity cloud, a shock wave is emitted as shown in picture 5 at the same time. Then, the shock wave propagates in the flow field both upstream and downstream. When the shock wave meets with the newly grown cavity sheet, the cavity sheet will be sharply destroyed due to the interaction between the shock wave and the cavity sheet. Then, the shock wave rebounds. As can be seen in the cavity volume,during the shock wave emission and propagation state,the cavity volume decreases quickly and in the rebound process, the cavity volume fluctuations can be observed, indicating a high mass transfer rate and a strong cavitation instability.

3.2 Vortex dynamics in time-dependent sheet/cloudcavitating flows

In the unsteady cavitating flows, in the stage of the attached cavity growth, plenty of liquid evaporates with the mass transfer from the liquid to the vapor.The cavitating fluid volume increases rapidly. As shown in Fig. 4(a), the cavity attaches to the foil leading edge (LE) and grows towards the trailing edge(TE) to a size of about 0.5c , where the cavity interface marked by the red dashed line indicates the contour line of. The 3-D view of the cavity structure is shown in picture 1 in Fig. 2. The vorticity Z (Zω ) contour on foil mid-plane is shown in Fig.4(b) together with the streamlines. The foil leading edge (FLE), the foil trailing edge (TEV), the cavity leading edge (CLE), the cavity trailing edge (CTE)and the cavity wake (CW) are marked. It can be observed that the vorticity concentrates in the cavitation region, the cavity wake region, at the foil trailing edge and in the foil wake. A vortex pair at the foil trailing edge labeled “TEV” and a foil wake sheet vortex in the foil wake “WVS” are seen in Fig. 4(b),with a similar magnitude as the upper vortex and the down vortex. The line pair shape of the WVS indicates the weak vortex in the foil wake region in this stage. A separated vortex pair shape TEV is observed at the foil trailing. The Q- criterion method in Fig. 4(c) can predict the vortex in the cavitation region and at the foil trailing edge. The Ω method gives contours distributed mainly in the cavitation region, outside the cavity interface, at the foil trailing edge, and in the foil wake region. It can be seen that

Fig. 4 (Color online) Void fraction distribution, mass transfer and comparisons of different vortex identification methods including the vorticity, and Rortex methods, at during the attached cavity growth. The red dashed line indicates the contour line of =0.2

compared with the Q - criterion method, the Ωmethod can predict the regular cavity interface-like vortex structures, the thickness of which is comparable with that of the attached cavity sheet,outside the cavity (OCV, =0.20 vα ), the strip vortex structures on the foil pressure side (FPV), and in the foil wake region (WVS). The thickness of the cavity shape-like vortex structures is larger than that of the strip vortex structures as shown in Fig. 4(d). It should be noted that the value of Ω in the cavitation region is smaller than that outside the cavity, indicating that in the cavitation region, the strain rate is comparable with the rotating rate. The regular large Ω value distribution is found outside the cavity, indicating the strong rigid body rotating. The2λ method is presented in Fig. 4(e), with contours mainly distributed at the foil leading edge, in the attached cavity region, at the cavity trailing edge, in the cavity wake region, and at the foil trailing edge. In the experiment, the vorticity distribution is identified at the cavity sheet rear part, and in the wake region[2-3],and it can be seen that the Q- criterion method, the Ω method and the2λ method can predict the vortex in the cavity wake. The cavity wake vortex structures and the positions predicted by different vortex identification methods are different. The positions of the vortex structures in the cavity wake region predicted by the2λ method are in the rear of those predicted by the Ω method, while the vortex structure predicted by the Q-criteria is thin and attaches on the foil suction side as shown in Fig. 4(c).Moreover, CTE are well predicted by the2λ method.The vortex structures in the foil wake region are identified by the Ω method, but not by the Qcriterion method and the2λ method. On the other hand, some noises can be seen in the predictions obtained by the Ω method as shown in Fig. 4(d). As shown in Fig. 4(f), the Rortex method ( R) can well predict the vortex structures both in the cavitation region and outside the cavity interface. All these vortex methods can predict the TEV, while the direction of the vortex pair is only predicted by the vorticity method. The vortex pair shape WVS in the foil wake is only identified by the vorticity method,the Ω method and the Rortex method. It should be noted that the Rortex is a vector, and to show the vortex strength to be comparable with that predicted by the Q , Ω methods, the strength of the(whereis presented in the present study.

Fig. 5 (Color online) Variance of water fraction (black), streamwise velocity ( u /U , red) and pressure (blue) along the mid-plane on foil suction side in three instances illustrating the re-entrant jet movement with water fraction contours inserted

3.2.2 Re-entrant jet development stage

To illustrate the re-entrant jet development clearly, Fig. 5 presents a series of snapshots to show the reverse flow beneath the attached sheet cavity along with the water fraction, the streamwise velocity ( u / U) and the pressure. The monitors are located along the foil suction side on the foil mid-plane, and extends to a size of 1.5c downstream the foil. The black lines in the water fraction contours indicate the velocity. It can be found that the re-entrant jet is mainly composed of liquid.With the re-entrant jet movement, the liquid jet moves upstream, with the decrease of the vapor fraction.

In the re-entrant jet development stage, the cavity is relatively long, and a liquid reverse flow can be found beneath the attached cavity. As shown in Fig. 6(a), the cavity attaches on the foil suction side and grows to cover the whole foil surface, with the re-entrant flow beneath the cavity arriving at 0.4c. The cavity interface is marked by red dashed lines, and the 3-D view is presented in picture 3 in Fig. 2. The purple lines indicate the region of the re-entrant jet strength ofand the reverse flow region is also labeled by purple lines/ =0 u U . From the streamlines in Fig. 6(b), it can be found that the streamlines roll up at the rear of the cavity, indicating a strong vortex. Vortex structures in the attached cavity is in the opposite direction with that inthereverseflowregion(e.g.,there-entrantjet) as labeled by “- ” and “+”. The trailing vortex pairs shape TEV and vortex sheet shape WVS can be seen in the foil wake region. It can be found that the vortex structures in the cavitation region can be identified by the Q- criterion method, while the vortex strength in the re-entrant jet region is weak and cannot be identified by the Q- criterion method, as shown in Fig. 6(c). It can be also found that the vortex structures concentrate in the front part and the rear part of the cavity while they are sparse in the cavity middle section as predicted by the Q- criterion method. The Ω method can identify the OCV, WVS,and the FPV as shown in Fig. 6(d). The WVS presents the line-pair shape as that during the attached cavity growth stage. Similar within the attached cavity stage,a high regular Ω value region can be found outside the cavity, indicating the strong rigid body rotating region outside the cavity, the thickness of which is in the order of the attached cavity thickness. The2λ method can identify the CTV and within the cavitation region as illustrated in Fig. 6(e). The Rortexmethodcombines the advantages of both the Q, Ω method.The OCV, FPV, and WVS are identified as shown in Fig. 6(f). It can be observed that during the re-entrant jet development stage, the vortex structures in the reverse flow region are weak, including the re-entrant jet region, mainly consisting of liquid. However, the vortex structures along the interface between the reverse flow region and the cavity are strong. The vortex structures within the cavitation and outside it,especially at the cavity leading edge and the trailing edge, are strong.

3.2.3 Cavity cloud shedding stage

Fig. 6 (Color online) Void fraction distribution, mass transfer and comparisons of different vortex identification methods, the Q , Ω and Rortex methods, at during the re-entrant jet development

In the cavity cloud shedding stage, the cavity structures are characterized by the newly attached cavity growth, and the large scale cavity cloud shedding as shown in Fig. 7(a), picture 4 in Fig. 2. As shown in Fig. 7(a), the cavity cloud sheds at the foil trailing edge. The cavity interface is outlined by thereddashedlinesinFigs.7(a)-7(c),7(f)andbythelight blue dashed lines in Figs. 7(d), 7(e). As illustrated in the vorticity distribution and the streamlines in Fig. 7(b), the streamlines roll up in the shedding cavity cloud region, and a shedding cavity cloud vortex pair (SCV) can be seen, indicating a strong vortex. In the newly attached cavity growth region, the negative vorticity can be found as labeled in Fig. 7(b). Unlike the line-pair shape of the WVS in the attached cavity growth and the re-entrant jet development, in this cavity cloud shedding process,the WVS presents a vortex pair shape, indicating a stronger vortex. As shown in Fig. 7(c), the Qcriterion method can identify the SCV, TEV and in the newly growth attached cavity. The Ω method can predict the vortex structures in the newly grown cavity (NCLV) and, the SCV, and the TEV, and the WVS. It should be noted that the vortex structures in the newly grown cavity sheet region are different from those during the first stage attached cavity sheet growth, showing a different mechanism in the cavity sheet growth, and thus a cavitation dynamics. TheWVS structure presents a vortex pair shedding shape,indicating that the transport of the large scale cavity cloud will enhance the vortex strength downstream.As can be seen in Fig. 7(d), identified by the Ω method, the FPV structure is connected with the TEV.The2λ method can identify the CW, the SCV and the TEV, as shown in Fig. 7(e). As can be seen in Fig.7, the vortex structures and positions in the cavity wake identified by the Q- criterion method, the Ω method, the2λ method and the Rortex method are different. The Rortex method combines the advantages of the Q, Ω methods, and can predict the NCLV, the TEV, the FPV, the WVS and the SCV as shown in Fig. 7(f). Unlike the prediction by the Ω method, when normalized, the Rortex method can not only predict the vortex structure, but also the vortex strength. It can be also found that the strength distributions predicted by the Q , R methods are similar, but different from those predicted by the Ω method. The Rortex method can identify, a strong vortex SCV and TEV.

Fig. 7 (Color online) Void fraction distribution, mass transfer and comparisons of different vortex identification methods, the Q , Ω and Rortex methods, at during the cavity cloud shedding

Fig. 8 (Color online) Time evolution of pressure signals on foil suction side at x/c = 0.3c , 0.4c , 0.6c , 0.8c , 1.0c ,1.2c and 1.4c during the shock wave formation and propagation process, with the flow structures in three instances

Fig. 9 (Color online) Void fraction distribution, mass transfer and comparisons of different vortex identification methods, the Q , Ω and Rortex methods, at t = t 0 +6204Tr ef during the cavity cloud collapse

3.2.4 Cavity cloud collapse induced shock wave propagation stage

Figure 8 shows the time evolution of the pressure signals along with the cavity structures during the shedding cavity cloud collapse induced shock wave formation and propagation. The dashed lines indicate the instances corresponding to the cavity structures on the right. It can be seen that the pressure peaks in the whole flow field are captured after the cavity cloud collapse. After the shock wave passes by, it is seen to rebound.

Figure 9(a) shows the flow structures during the cavity cloud collapse. A thin cavity sheet can be seen at the foil leading edge, and a collapsing cavity cloud can be observed in the high pressure region downstream the foil, which is also shown in picture 5 in Fig. 2. As shown in the vorticity contour and the streamlines in Fig. 9(b), especially, the streamlines in the foil wake region, the cavity cloud collapses,indicating the vortex structures there. Vortex pairs can be observed the TEV at the foil trailding edge and the SCV at the cavity cloud collapse region. Negative vorticity occurs at the newly grown attached cavity region as labeled in Fig. 9(b). It is indicated that the vorticity concentrates in the foil wake region and at the foil leading edge in the cavitation region, and at the trailing edge. It can be found that all vortex identification methods predict the TEV during the cavity cloud collapse stage, however, with different vortex structure and strength. Vortex directions can be predicted by the vorticity method. The Q- criterion method in Fig. 9(c), the Ω method, and the Rortex method can identify the NCLV, but not by themethod as indicated in Fig. 9(e). High value of Ω can be found at the foil leading edge where the new cavity sheet grows, with a cavity wake, showing a strong rigid body rotation around the newly grown cavity sheet. The SCV is identified by the vorticity method, the Q- criterion method, the Ω method,themethod, and the Rortex method. The vortex structure in the cavity wake region is predicted by thecriterion method, the Ω method, and the Rortex method with different vortex structure and strength.The strike vortex structure FPV is identified by the Ω method in Fig. 9(d) and the Rortex method in Fig.9(f). Unlike the cavity cloud shedding stage where the shedding cavity cloud is just above the foil suction side, the large scale vortex structures are found in the foil wake region, far larger than the cavity region. The large scale vortex structures at this stage are induced by the shedding voritical cloud. Moreover, the collapsing cavity cloud will generate the vorticity and enhance the vortex strength[2-3], leading to larger vortex structures than the cavitation region.

Fig. 10 (Color online) Comparisons of the Ω method with different ε values, 7 max ε =10- (b - a) , 6 max ε =10- (b - a) ,5 max ε =10- (b - a) and 4 max ε =10- (b - a) , during the re-entrant jet development in Fig. 6(d). The red dashed dot lines show the cavity interface = 0.2 v α and the pink dashed dot lines show the reverse flow region and re-entrant jet region/ =0 x u U , / = 0.56

3.3 Improvement of Ω method based on small parameterε

According to Liu et al.[33], when a, b in the Ω method is small enough as nearly zero, the value of Ω might be very large owing to the fact of being divided by a number nearly zero, thus resulting in noises or clouds, as shown in Fig 4(d). Therefore, to accurately identify the vortex structures in the field,the improvement of the Ω method is necessary, in particular when exploring new vortex structures without prior knowledge about the vortex. The small parameter method is one of the most popular method for the improvement of the Ω method. To link the small parameter to the flow field, the small parameter[27-28]is defined asmax. Then,the influence of the small parameter on the vortex structure identification performance of the Omega method is evaluated in the process of the re-entrant jet development and the cavity cloud collapse induced shock wave formation and propagation in the sheet/cloud cavitating flows.

3.3.1 Influence of ε on the Ω method in the process of re-entrant flow development

It is known that the re-entrant jet is one of the main cause of cavitation instabilities[1,19], and the shear force between the re-entrant jet and the cavity is complex and could generate bubbly vortical structures during the re-entrant jet movement. To investigate the coherent interactions of the vortex dynamics and the re-entrant jet dynamics of the cavitation characteristics, the influence of the small parameter on the vortex identification during the re-entrant jet development is important. Figure 10 shows the vortex identification performance of the Ω method with different small parameters (, ε=,and4(during the re-entrant jet development as shown in Fig. 6(d). With the increase of the small parameter,the vortex identified by the Omega method has different structures. As shown in Fig. 10(b), compared with that at, the FPV, OCV and WVS decrease in number and shrink at ε=While those along the interface between the re-entrant jet and the cavity change only a little. At, the line-shape vortex structures in the foil wake region and the strike vortex like FPV disappear, and the OCV with a strong rigid body rotating greatly decrease in number, while the vortex structures within the cavitation region change little, as shown in Fig. 10(c). The vortex structuresalongtheinterfacebetweenthere-entrantjetandthecavity become weak. With a further increase of the small parameter to, the vortex structures outside the cavity disappear, and the vortex structures within the cavitation region decrease in number as shown in Fig. 10(d). The vortex structures along the interface between the re-entrant jet and the cavity become further weak. It can be seen that the sensitivity of different vortex structures to the small parameter changes is different. The FPV, WVS are more sensible than the OCV. The vortex structure is the least sensible to the small parameter within the cavitation region, especially the CLV and the CTV. It can be seen that to identify certain vortex structures,the choice of a proper small parameter is critical. An improper selection of ε will generate noise, and it is suggested to increase the value of epsilon to remove the noise. On the other hand, to capture more vortex structures, it is suggested that the small parameter is not very large. In the present study,,which is proper.

3.3.2 Influence of ε on the Ω method in the process of cavity cloud collapse induced shock wave propagation

Fig. 11 (Color online) Comparisons of the Ω method with different ε values, 10 max ε =10- (b - a) , 9 max ε =10- (b - a) ,8 max a) , 7 max a)max a)max ε =10- (b - a) , during the cavity cloud collapse induced shock wave formation and propagation in Fig. 9(d)

In the shock wave mechanism, the cavity cloud collapse induced shock wave formation and propagation is the main cause of cavitation instabilities,b-which will induce large pressure fluctuations, as well as the noise emission. According to the experimental measurement in the closure region of the attached cavitation[3], the cavity collapse is an important source of the vorticity generation. It should be noted that the shock wave formation is caused by the large scale cavity cloud collapse and during the shock wave propagation, the interactions between the shock wave and the cavity will cause a local cavity collapse,resulting in a low void fraction[7]. Thus, in the process of the shock wave generation and propagation, the coherent interactions between the vortex dynamics and the shock wave dynamics are essential for understanding the physics and the mechanism of cavitation instabilities. In view of the fact that during the shock wave formation and propagation, the cavity structure changes abruptly due to the shock wave, the vortex dynamics is supposed to be different from that in other cavitation evolution stages, e.g., during the attached cavity growth, the re-entrant jet development,and the cavity cloud shedding. Consequently, the vortex identification is critical for the cavitation dynamics involved study. In what follows, the influence of the small parameter in the Ω method on the vortex identification performance is presented.Figure 11 shows the vortex structures identified by the Ω method using different small parameters,,,b-,and. Different vortex structures (the NCLV, the CW, the OCV, and the WVS) are labeled by yellow arrows as shown in Fig. 11(a). Several noises or clouds can be found as marked by red dashed circles, which are always connected with the walls. With the small parameter being decreased to, the noise structures shrink as shown in Fig. 11(b). At, the noises further decrease, and at the same time, the weak vortex structures in the foil wake region begin to shrink, while the vortex structure with a relative strong strength changes little as shown in Fig. 11(c).With the small parameter being further decreased to, the noise structures shrink a great deal and disappear. At the same time, the wake vortex structure in the foil wake region attenuates as indicated by the yellow dashed circles as shown in Fig.11(d), while the vortex structures at the foil leading edge and at the trailing edge change slightly. Atas shown in Fig. 11(e), the vortex structures at the foil leading edge and the trailing edge decrease a great deal in number and at ε=, besides the reduction of the vortex structures at the foil leading edge and the trailing edge,the vortex near the foil wake begins to shrink as shown in Fig. 11(f). It can be shown that the influences of the small parameter on the noise and the vortex structures are different, where a small change of the small parameter will much change the noise. On the other hand, a weak vortex structure is more sensible to the small parameter variation. To identify both the weak vortex and the strong vortex, the small parameter is critical and should be small enough. In the present study,is proper.However, the distinction between the noise and a weak vortex should be based on the understanding of the cavitating flow field.

Fig. 12 (Color online) Comparisons of different vortex identification methods including vorticity method, Q - criterion method, and Ω method

Fig. 13 (Color online) Comparisons of different vortex identification methods including vorticity method, Q -criterion method,Ω method and a, b, c ( c = b - a) distributions in Ω method in the process of re-entrant jet development

3.4 Quantitative analysis of vorticity method, Q -criterion method and Ω method

According to the definition, a close relationship can be found between the vorticity method, the Qcriterion method and the Ω method. To make quantitative comparisons between the vorticity method, the Q- criterion method and the Ω method,Fig. 12 shows the comparisons of different vortex identification methods, with x- axis representing, y- axis representing. It can be found that the results of the vorticity method are independent ofand linear withas shown in Fig.12(a). In the present study, the range of a and b is 0-1 000. The results of the Q-criterion are proportional to,as shown in Fig. 12(b). The results of the Ω method are non-linear with,and with the increase of,, the Ω value increases. From the results identified by the vorticity method,theQ-criterionmethod,andtheΩ method, it can be found that different methods can identify vortices in different regions in theplane, showing different identification performances for the vortex structures. The vortex strength can be evaluated by the vorticity method and the Q- criterion method based on, while with the Ω method, a normalized parameter is used, and difference between smallcan not be distinguished, therefore, the vortex strength can not be evaluated. Compared with the vorticity method,the Q- criterion method employs the stain rate sensor,with consideration of the effects of a on the vortex structures. The Ω method uses the normalized parameter, defined as the ratio of the rotation rate over the summation of the rotation rate and the strain rate.From the distribution of the high Ω value, it can be observed that the Ω method can identify both the strong vortex (high b value) and weak vortex (low b value). To further evaluate the vortex identification performance of different vortex identification methods, large in the cavitating flows, the distributions of a, b and c (c =b - a) within the cavitating flow field are discussed in the re-entrant jet development stage and the shock wave formation and propagation stage.

Fig. 14 (Color online) Comparisons of different vortex identification methods including vorticity method, Q -criterion method,Ω method and a, b, c distributions in Ω method in the process of shock wave formation and propagation

3.4.1 Re-entrant jet development

Figure 13 shows the distributions of Ω, a, b and cin the Ω method, the vorticity method and the Q- criterion method in the process of the re-entrant jet development. It can be found that high values of a, b are found within the cavitation region, at the cavity leading edge and the cavity trailing edge, where the vortex structures, the CLV,CTV, are localized, relatively low values of a, b are found outside the cavity, in the foil wake, on the foil pressure side, where the vortex structures, OCV,WVS and FPV, are positioned. With further evaluations of the vortex strength, combined with the vorticity method and the Q- criterion (and c), it can be found that the foil wake is the high vorticity concentration regionasshowninFig.13(e),butthevortex strength defined by the Q- criterion method(and c) is weak as shown in Fig. 13(f). The vortex structures within the cavitation region, at the cavity leading edge and the cavity trailing edge are strong with both high values of a, b, where the vorticity concentrates and a high Q value (and c) can be found. The vortex structures outside the cavity and on the foil pressure side are of low values of a, b,where high values of Ω a are observed. In the re-entrant jet region, low values of a, b, c, Ω and Q are found.

3.4.2 Shock wave formation and propagation

Figure 14 shows the distributions of Ω, a, b and c in the Ω method, the vorticity method and the Q- criterion method in the process of the shock wave formation and propagation. In this process, the large scale shedding cavity cloud will induce large scale vortex structuresWVS as can be predicted by the vorticity method and the Q- criterion method shown in Figs. 14(e), 14(f), respectively. It can be found that high values of a, b distribute mainly at the newly grown attached cavity sheet e.g., FLV, in the shedding cavity cloud region e.g., WVS. With the increase of the distance to the foil, the vortex strength attenuates as predicted by both the vorticity method and the Qcriterion method (and c) in Figs. 14(e), 14(f), where low values of a, b are found as shown in Figs.14(b), 14(c). Owing to the flow separation at the foil trailing edge, high vorticity can be found as shown in Fig. 14(e). High values of a, b and c, and Qcriterion are found in FTV. Outside the cavity region,low values of a, b are observed, as is similar to the case of the re-entrant jet development stage.

4. Discussions and conclusions

Numerical simulations of the unsteady sheet/cloud cavitating flows around a NACA66 (mod)hydrofoil are conducted based on the implemented compressible cavitation solver[37]in the open source software OpenFOAM. The foil is fixed with an angle of attack α=6° , σ=1.25, =5.33 m/sand Re. In the simulations, the set of governing equations for the compressible cavitating flows including mass, momentum and energy conservation equations are solved, with considerations of the cavitation fluid compressibility and the heat transfer between phases. The sheet/cloud cavitating flows include four phases: 1) the attached cavity growth, 2)the re-entrant jet development, 3) the cavity cloud shedding, and 4) the cavity cloud collapse induced shock wave formation and propagation, with a strong coherent coupling between the vortex dynamics and the cavitation dynamics. Different vortex identification methods including the vorticity method, the Qcriterion method, the Ω method, the2λ method,and the Rortex method are employed for the vortex structure identification in the transient cavitating flows. The main findings are as follows:

(1) The cavitation vortex structures are closely related with the transient cavity behaviors. During the attached cavity sheet growth stage, the vorticity concentrates in the CW and the WVS. The vortex structures within the cavitation region, the TEV, and the CW can be identified by the Q- criterion method,the Ω method, the2λ method, and the Rortex method. A regular strong rigid body rotating region the OCV can be identified by the Ω and Rortex methods, and a strip vortex structure is found on the FPV. The2λ method can well predict the vortex in the cavity closure region. The line-shape vortex structures in the foil wake can be identified by the vorticity method, the Ω method and the Rortex method. However, the vortex positions and structures in the cavity wake region identified by different methods are different, showing the difficulties of the vortex structure identification. It should be noted that compared with other vortex identification methods,the noise could be produced using the Ω method with an improper choice of the small parameter ε.

(2) During the re-entrant jet development stage,the flow structure is characterized by a reverse liquid flow beneath the cavity. Vortex structures distribute within the cavitation region, along the interface between the re-entrant jet and the cavity, and the CTV.Line-shape vortex structures the WVS can be identified. Within the re-entrant jet, mainly composed of liquid, the vortex structures are weak. The Qcriterion method can better predict the vortex structures within the cavitation region, and the2λ method in the CTV. The OCV and the FPV can be identified by both the Ω method and the Rortex method.

(3) During the cavity cloud shedding stage, the flow structure consists of the newly grown cavity sheet and the large scale shedding cavity cloud. Based on the vortex structure identification, it can be found that the NCLV are different from the CLV during the attached cavity sheet growth stage. The WVS have the vortex pair shape as identified by the Ω method and the Rortex method. The CW are complex and different vortex identification methods show different performances in both structures and positions. The shedding cavity cloud consists of strong vortex structures, along with the strong vortical structures outside the cavity cloud. Moreover, strong vortex can be found at the TEV, and in this process the strip vortex structures on the (FPV connects with the vortex structures at the TEV.

(4) During the cavity cloud collapse induced shock wave formation and propagation stage, strong vortex-pair shape vortex structures can be found in the WVS, with the size of the vortex far larger than that of the shedding cavity. The vortex structures, large in both size and strength in the foil wake region are supposed to be induced by the cloud cavity shedding behavior. Being exposed to a wetted flow, the flow separation bubble can be found at the foil trailing edge.The vortex structures at the foil leading edge, covered by the newly grown attached cavity sheet can be identified by the Q- criterion method, the Ω method and the Rortex method.

(5) Compared with the Q- criterion method and the2λ method, the Ω method and the Rortex method can identify the vortices outside the cavity, on the foil pressure side and in the foil wake region where the vortex strength in different regions is different. The Ω and Rortex methods can capture both strong vortices and weak vortices, or in other words, can capture more vortex structures, while the Q,2λ methods can only capture strong vortices but fail to identify weak vortices when a certain threshold is set. The small parameter ε can be used to improve the performance of the Ω method and reduce the noise or clouds. In the re-entrant jet movement process, with increasing the small parameter value, the vortex structures in the foil wake region, outside the cavity, and on the foil pressure side shrink gradually,following the cavitation region and the cavity leading edge and the trailing edge. In the shock wave formation and propagation process, with increasing the small parameter, the noise related with the wall and the weak vortex structures in the foil wake is reduced quickly. Compared with strong vortex structures, it can be found that the noise and the weak vortex are more sensitive to the small parameter.However, the distinctions between the noise and the weak vortex are based on the understanding of the cavitating flow field. To identify the weak vortex, the small parameter valueis suitable.Combined with the a, b and c analyses in the Ω method, it can be concluded that large values of a, b are found in the cavitation region, and small values of a, b are found in the OCV, the FPV, and the CW. For the unsteady sheet/cloud cavitating flows,the cavitation vortex identification in the cavity closure region and cavity wake region in both structures and positions is a challenge, with great discrepancies of the results obtained by different vortex identification methods.

Further investigations of the cavitation vortex dynamics should be conducted. Recently, the Ω -Liutex/Rortex method with the advantages of both the Ω and liutex/Rortex methods was proposed[43]. It is promising to employ the new vortex identification method for the unsteady cavitating flows. At the same time, the vortex structure study from the Lagranian point of view is also important. The Lagrangian methods, including the coherent structures (LCS) and the particle trajectory, will further promote the vortex dynamics investigation.

Acknowledgements

This work was supported by the Open Foundation of State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University, China (Grant No.1611), the Open Fund for the Key Laboratory of Fluid and Power Machinery, Ministry of Education, Xihua University.