Large Eddy Simulation of Turbulent Attached Cavitating Flows around Different Twisted Hydrofoils

: In this paper, the turbulent attached cavitating ﬂows around two different twisted hydrofoils, named as NACA0009 and Clark-y, are studied numerically, with emphasis on cavity shedding dynamic behavior and the turbulence ﬂow structures. The computational method of large eddy simulation (LES) coupled with a homogeneous cavitation model is applied and assessed by previous experimental data. It was found that the predicted results were in good agreement with that of the experiment. The unsteady cavity morphology of the two hydrofoils undergoes a similar quasi-periodic process, but has different shedding dynamic behavior. The scale of the U-type shedding structures forming on the suction surface of NACA0009 is larger than that of Clark-y. This phenomenon is also present in the iso-surface distributions of Q -criterion. Otherwise, the time-averaged cavity morphology is dramatically different for the two hydrofoils, and it is found that the attached location of the cavity is closely related to the hydrofoil geometry. The time ﬂuctuation of the lift force coefﬁcients is affected signiﬁcantly by the cavity shedding dynamics. Compared with NACA0009, the lift force of Clark-y shows more ﬂuctuation, due to its complicated shedding behavior. Further analysis of the turbulent structure indicates that the more violent shedding behaviors can induce higher levels of turbulence velocity ﬂuctuations.


Introduction
Cavitation is a very important flow phenomenon that influences the design and operation of hydraulic machines [1][2][3][4], such as pumps, ship propellers, and turbomachinery. A better understanding of the cavitating flow mechanism is substantially necessary for effective cavitation control in engineering designs [5].
Owing to the great efforts of many researchers in the past, knowledge of turbulent cavitating flows have been remarkably improved recently. Some attention has been focused on the time-evolution features of unsteady cavitating flows around various physical models. For example, Wang et al. [6] experimentally observed the unsteady cloud cavitating flow structures around a two-dimensional Clark-y hydrofoil; Barre et al. [7] studied the attached sheet cavitating flows in venture; Wu et al. [8] studied the transient characteristics of cloud cavitating flow around a flexible hydrofoil; Gopalan and Katz [9] investigated the flow structure in the closure region of attached cavitating flow in nozzles; Foeth [10] made numerous experiments to study the sheet/cloud cavitating structures around different

Numerical Methods and Validation
The present work applied the LES method, coupled with the cavitation model proposed by Zwart et al. [38], in order to simulate unsteady cavitating flows under the frame of a homogenous multiphase flow assumption. The details of the solver are given as follows.
In the homogenous multiphase flow, each fluid phase is assumed to share the same velocity and pressure field. The basic governing equations contain the mass and momentum conservation equations, shown as Equations (1) and (2): where u with different subscripts i and j presents the velocity component in different coordinate directions, ρ m is the mixture density, p is the pressure, and µ m is the mixture dynamic viscosity. For the LES approach, after a Favre-filtering operation to Equations (1) and (2), the LES basic equations are given as below: ∂ρ m ∂t where the variables with over-bars indicate the filtered quantities. Another new term τ ij , naming the sub-grid scale stresses, is added in Equation (4), and is defined as: A widely used modeling method is the eddy-viscosity model, which assumes that the sub-grid scale stresses τ ij are proportional to the modulus of the strain rate tensor S ij of the filter's large-scale flow. The following equation denotes this relationship: where µ t is the sub-grid scale turbulent viscosity, which is closed by the wall-adopting local eddy (WALE) viscosity model [39] in this work. The corresponding definitions in the WALE model are as follows: S d ij = g ij = ∂u i ∂x j , L s = min kd, C w V 1/3 (10) Energies 2018, 11, 2768 4 of 15 where L s is the mixing length for the sub-grid stresses, k is von Karman's constant, d represents the distance to the closest wall, and V is the volume of the computational cell. The default WALE constant C w is set to 0.5. For the cavitation process, the mass transfer equation is expressed as below: where α l is the water volume fraction. The source terms m + and m − denote the condensation and evaporation rates. In the present work, the source terms are expressed by the Zwart cavitation model. This model is derived from a simplified Rayleigh-Plesset equation, neglecting the second order derivative of the bubble radius, the surface tension term, and the viscous term. Then the source terms m + and m − are defined as: where F c and F e are two empirical coefficients for the condensation and vaporization processes, and α nuc and R B respectively represent the volume fraction and the radius of the nucleation sites. These empirical constants in the ANSYS CFX (18.0 Version, ANSYS, Canonsburg, PA, USA) software [40] are set to F c = 0.01, F e = 50, α nuc = 5 × 10 −4 , and R B = 1 × 10 −6 m. The Zwart cavitation model has been widely used to simulate cavitating flows, and more descriptions about this popular numerical method can be found in [41][42][43][44]. Brief introductions of the simulation cases are listed in Table 1. Case 1 is used to validate the numerical method again. The other two cases with the same flow condition are simulated to investigate the differences of unsteady cavitating flows between the two hydrofoils. To be consistent with the experiment setup of [10], in the simulation of Case 1, the inflow velocity is 6.97 m/s and the cavitation number is 1.07. The computational domain and the boundary conditions are similar to the description in the literature [37]. Half of the flow domain is computed because of the geometric symmetry. Thus, the domain is with the width of C 1 , the height of two C 1 , and the length of seven C 1 . Here, C 1 with the value of 0.15 m is the chord length of the Delft Twist-11 hydrofoil. Figure 1a,b respectively gives the observed and numerical unsteady cavitation structures around the Delft twisted hydrofoil in one typical cycle. Here, in the predicted results, the cavity morphology is expressed by the iso-surface of the water vapor volume fraction, with a value of 0.1. From Figure 1a,b, it is found that the experiments and simulations can be remarkably consistent for capturing the periodical evolution process of attached turbulent cavity around the hydrofoil, including the attached cavity growing up, breaking off, primary shedding, secondary shedding, and final collapse downstream. In addition, Table 2 gives the typical cavity shedding frequency f and the lift coefficient C L obtained by several researchers. Here, the predicted frequency f is obtained by a fast Fourier transform (FFT) analysis of the total vapor volume. The lift coefficient C L is defined as Equation (14). When compared with other documents, the data acquired in this work agree well with that of prior experiments. Consequently, the above analysis can reasonably indicate the present LES

Cases
CL f (HZ) EXP [10] 0.53 32.2 LES [45] 0.44 34 LES, fine [45] 0.45 -DES [45] 0.42 30 LES (in this work) 0. 45 34 Based on the numerical method mentioned above, this work will contribute to numerical studies of unsteady cavitating flows around two twisted hydrofoils, with special emphasis on the differences of cavitation characteristics between them. Figure 2 gives the geometry of two twisted hydrofoils with different sections, named NACA0009 and Clark-y, respectively. The two hydrofoils have the same chord length of C = 70 mm and span length of S = 70 mm. The only difference between the two twisted hydrofoils is the sectional profile. One has a symmetric NACA0009 section, and the other has an asymmetric Clark-y section. As shown in Figure 3, the sections of two hydrofoils vary with their angles of attack from 0° to 11° along the span direction, which is similar to the Delft twisted hydrofoil [46].    [45] 0.45 -DES [45] 0.42 30 LES (in this work) 0. 45 34 Based on the numerical method mentioned above, this work will contribute to numerical studies of unsteady cavitating flows around two twisted hydrofoils, with special emphasis on the differences of cavitation characteristics between them. Figure 2 gives the geometry of two twisted hydrofoils with different sections, named NACA0009 and Clark-y, respectively. The two hydrofoils have the same chord length of C = 70 mm and span length of S = 70 mm. The only difference between the two twisted hydrofoils is the sectional profile. One has a symmetric NACA0009 section, and the other has an asymmetric Clark-y section. As shown in Figure 3, the sections of two hydrofoils vary with their angles of attack from 0 • to 11 • along the span direction, which is similar to the Delft twisted hydrofoil [46]. hydrofoils with different sections, named NACA0009 and Clark-y, respectively. The two hydrofoils have the same chord length of C = 70 mm and span length of S = 70 mm. The only difference between the two twisted hydrofoils is the sectional profile. One has a symmetric NACA0009 section, and the other has an asymmetric Clark-y section. As shown in Figure 3, the sections of two hydrofoils vary with their angles of attack from 0° to 11° along the span direction, which is similar to the Delft twisted hydrofoil [46].    Figure 4 shows the setup of computational domain and mesh around the hydrofoil. From Figure 4a, the hydrofoil is placed horizontally in the domain, with a width of 0.5C, a height of 2.7C, and a length of 10C. It should be noted that only half of hydrofoil was considered, due to its geometric symmetry. The corresponding boundary conditions are listed in Table 3. According to the velocity of = 7 m/s and the pressure of 32,510.2 Pa, the cavitation number σ is equal to 1.2.
The grid distribution around the hydrofoil is shown in Figure 4b. An O-H type mesh was applied and refined near the foil surface. In order to find a balance between computation cost and accuracy, three sets of mesh with the same grid topology were tested. As listed in Table 4, the meshes with different spatial resolutions were evaluated by two variables, CL and CD. The time-averaged lift coefficient CL and drag coefficient CD of the NACA0009 hydrofoil are defined as below: where FL and FD are lift force and drag force, respectively; the other variables can be found above.   The upper, lower, and back side face Free-slip -Wall The foil surface No-slip -  Figure 4 shows the setup of computational domain and mesh around the hydrofoil. From Figure 4a, the hydrofoil is placed horizontally in the domain, with a width of 0.5C, a height of 2.7C, and a length of 10C. It should be noted that only half of hydrofoil was considered, due to its geometric symmetry. The corresponding boundary conditions are listed in Table 3. According to the velocity of U ∞ = 7 m/s and the pressure of 32,510.2 Pa, the cavitation number σ is equal to 1.2. The grid distribution around the hydrofoil is shown in Figure 4b. An O-H type mesh was applied and refined near the foil surface. In order to find a balance between computation cost and accuracy, three sets of mesh with the same grid topology were tested. As listed in Table 4, the meshes with different spatial resolutions were evaluated by two variables, C L and C D . The time-averaged lift coefficient C L and drag coefficient C D of the NACA0009 hydrofoil are defined as below: where F L and F D are lift force and drag force, respectively; the other variables can be found above.   Figure 4 shows the setup of computational domain and mesh around the hydrofoil. From Figure 4a, the hydrofoil is placed horizontally in the domain, with a width of 0.5C, a height of 2.7C, and a length of 10C. It should be noted that only half of hydrofoil was considered, due to its geometric symmetry. The corresponding boundary conditions are listed in Table 3. According to the velocity of = 7 m/s and the pressure of 32,510.2 Pa, the cavitation number σ is equal to 1.2.
The grid distribution around the hydrofoil is shown in Figure 4b. An O-H type mesh was applied and refined near the foil surface. In order to find a balance between computation cost and accuracy, three sets of mesh with the same grid topology were tested. As listed in Table 4, the meshes with different spatial resolutions were evaluated by two variables, CL and CD. The time-averaged lift coefficient CL and drag coefficient CD of the NACA0009 hydrofoil are defined as below: where FL and FD are lift force and drag force, respectively; the other variables can be found above.   The upper, lower, and back side face Free-slip -Wall The foil surface No-slip -Symmetry The front side face --  The upper, lower, and back side face Free-slip -Wall The foil surface No-slip -Symmetry The front side face -- In Table 4, the comparisons of C L and C D indicate that the small differences between the medium mesh and fine resolution mesh can be neglected. Besides, the y + of the medium mesh is nearly 1 for most of the hydrofoil surface. The definition of y + can be written as Equation (15). According to the previous comments for a valid LES method [37,44], a medium mesh with about 3.5 million nodes was selected as the final grid.
where ∆y is the distance to the nearest wall, u τ is the friction velocity at the nearest wall and υ is the kinematic viscosity. All the simulations were conducted using commercial CFD software named ANSYS CFX. The time-dependent cavitating flow simulation was performed with the initial condition of steady non-cavitation results. The high-resolution scheme was used for the convective term, and the second order backward Euler scheme was used for the transient term. The time step was set to 5 × 10 −4 s. Then, the computation of 800 time steps took about three days, using a computer with 64 GB memory and 24 parallel cores.

Results and Discussion
Figure 5a,b respectively gives some computational snapshots of cavity topology and the pressure distribution on the suction surface of NACA0009 and Clark-y in a typical cycle. Here, the iso-surface of water vapor volume fraction with a value of 0.1 was selected to depict the morphology of cavitation. The reference time is expressed by t 0 , and T is the period of unsteady cavity evolution. For the two hydrofoils, the time evolution of unsteady cavity morphology undergoes a similar process, including the attached cavity breaking off, and subsequently growing up again, as well as the primary shedding with a U-type vortex structure and the secondary shedding occurring at the two sides of the attached cavity tail. However, there are some distinct differences in the details of transient cavity morphology. Firstly, the fully developed attached cavity of NACA0009 with the regular convex closure completely covers the hydrofoil leading edge. In contrast, the attached cavity around the Clark-y hydrofoil with a larger scale is not covered the leading edge, especially at the two sides of the hydrofoil. Secondly, to focus on the primary shedding process, the U-type shedding vortex structure forms and immediately sheds downstream, with a larger scale for the NACA0009. However, for the Clark-y, the upstream part of the break-off cavity grows faster, and could catch up with the downstream part to join together again. This behavior is probably responsible for the final, smaller-scale structure of U-type shedding. closure of the attached cavity. However, by comparisons, it is found that the pressure of the NACA0009 is always larger than that of the Clark-y, except for the cavitation region. To focus on the leading edge, the pressure of Clark-y delays reaching the vapor pressure, which is probably attributed to its larger curvature at the leading edge. Consequently, it can be concluded that the body geometry substantially affects the pressure distributions, and then yields to the different attached cavity characteristics.  Figure 7 gives the comparisons of lift coefficients for the two hydrofoils. The time-dependent lift curves both show quasi-periodic characteristics, which is similar to the time evolution of cavity topology. However, compared with that of NACA0009, we observe that the transient lift coefficient signals of Clark-y present more small fluctuations over time, due to its more complicated cavity shedding process, as mentioned above. It is also found that the time-averaged value of the lift coefficient for Clark-y is much larger than that of NACA0009. Actually, the lift value is defined by the pressure difference of the pressure and suction surfaces of the hydrofoil. As for the cavitating  From Figure 6a, the time-averaged cavity topology for the two hydrofoils demonstrates the similar characteristics as that mentioned in the previous section. That is, for the NACA0009 the time-averaged cavity starts from the whole leading edge and curls into a convex closure. However, for the Clark-y, only a partial nearby area of the mid-plane of leading edge is covered by the cavity, which has a nearly straight closure. From the time-averaged pressure distributions, a larger area, with high pressure levels than that of the Clark-y, was found located behind the cavity colure of the NACA0009. It can be concluded that the high pressure restrains the growth of the cavity, in order to induce a smaller scale of cavity on the NACA0009. To further study the location difference of the cavity at the leading edge for the two hydrofoils, Figure 6b gives the geometry outline of a typical section and the corresponding pressure distributions for the two hydrofoils. The selected section is located at y/S = 0.1, marked by a red dashed line, as shown in Figure 6a. The profiles of pressure present a similar trend along the suction surface for the two hydrofoils, including decreasing dramatically at the leading edge of the vapor pressure and then increasing sharply behind the closure of the attached cavity. However, by comparisons, it is found that the pressure of the NACA0009 is always larger Energies 2018, 11, 2768 9 of 15 than that of the Clark-y, except for the cavitation region. To focus on the leading edge, the pressure of Clark-y delays reaching the vapor pressure, which is probably attributed to its larger curvature at the leading edge. Consequently, it can be concluded that the body geometry substantially affects the pressure distributions, and then yields to the different attached cavity characteristics.
Energies 2018, 11, x 9 of 15 flows, the pressure value of the suction surface is substantially affected by the cavitation area, which has low pressure. As a result, the larger cavitation region of Clark-y yields to lower pressure of suction surface, leading to the larger lift force. The cavity shedding behavior in unsteady cavitating flows is closely related to the reentrant flow and the attached cavity itself [47,48]. Then, in the present work, the velocity vector distributions on the suction surface were investigated during the cavity shedding process. In Figure 8, the cavity topology with iso-surface values of = 0.1 and the contour of u-velocity component are displayed to study the cavity shedding behaviors. It should be noted that the negative values of u-velocity can represent the reentrant flow behaviors. For the NACA0009 in Figure 8a, the reentrant flow invaded the attached cavity more widely with a radially-diverging pattern, which then caused the horseshoe vortex structure with larger-scale shedding from the attached cavity. In contrast, for the Clark-y in Figure 8b, it was found that the reentrant flow just invaded the center part of the suction surface, with a sharp corner head form. Moreover, after the first break-up, the upstream attached cavity developed quickly enough to overcome the collision of the weaker reentrant flow and join with the downstream part again. As a result, this conjunction behavior, as well as the reentrant flow, leads to the second breaking off at the downstream position and the smaller-scale of U-type shedding structure on the Clark-y.  Figure 7 gives the comparisons of lift coefficients for the two hydrofoils. The time-dependent lift curves both show quasi-periodic characteristics, which is similar to the time evolution of cavity topology. However, compared with that of NACA0009, we observe that the transient lift coefficient signals of Clark-y present more small fluctuations over time, due to its more complicated cavity shedding process, as mentioned above. It is also found that the time-averaged value of the lift coefficient for Clark-y is much larger than that of NACA0009. Actually, the lift value is defined by the pressure difference of the pressure and suction surfaces of the hydrofoil. As for the cavitating flows, the pressure value of the suction surface is substantially affected by the cavitation area, which has low pressure. As a result, the larger cavitation region of Clark-y yields to lower pressure of suction surface, leading to the larger lift force. The cavity shedding behavior in unsteady cavitating flows is closely related to the reentrant flow and the attached cavity itself [47,48]. Then, in the present work, the velocity vector distributions on the suction surface were investigated during the cavity shedding process. In Figure 8 Figure 8a, the reentrant flow invaded the attached cavity more widely with a radially-diverging pattern, which then caused the horseshoe vortex structure with larger-scale shedding from the attached cavity. In contrast, for the Clark-y in Figure 8b, it was found that the reentrant flow just invaded the center part of the suction surface, with a sharp corner head form. Moreover, after the first break-up, the upstream attached cavity The cavity shedding behavior in unsteady cavitating flows is closely related to the reentrant flow and the attached cavity itself [47,48]. Then, in the present work, the velocity vector distributions on the suction surface were investigated during the cavity shedding process. In Figure 8, the cavity topology with iso-surface values of α v = 0.1 and the contour of u-velocity component are displayed to study the cavity shedding behaviors. It should be noted that the negative values of u-velocity can  Figure 8a, the reentrant flow invaded the attached cavity more widely with a radially-diverging pattern, which then caused the horseshoe vortex structure with larger-scale shedding from the attached cavity. In contrast, for the Clark-y in Figure 8b, it was found that the reentrant flow just invaded the center part of the suction surface, with a sharp corner head form. Moreover, after the first break-up, the upstream attached cavity developed quickly enough to overcome the collision of the weaker reentrant flow and join with the downstream part again. As a result, this conjunction behavior, as well as the reentrant flow, leads to the second breaking off at the downstream position and the smaller-scale of U-type shedding structure on the Clark-y. As demonstrated by Gopalan and Katz [9], the bubbles at the closure region of the attached cavity presented hairpin-like structures. To more clearly demonstrate the turbulence flow structures during the shedding process, Figure 9 displays the iso-surface of the Q-criterion (Q = 4 × 10 5 ) for the two hydrofoils. It was found that the vortex shedding structures are closely related with the cavity morphology. The larger scale of the shedding cavity contributes to the larger scale of the U-type vortex and vice versa. During the cavity shedding process, the U-type vortex structures forming on the suction side are governed downstream by the main flow, and gradually collapse. The scale of the U-type vortex structure on the NACA0009 is larger than that of Clark-y. However, the downstream region of the cavity closure on the Clark-y contains more complicated turbulence structures, including numerous and correlative U-type structures.  As demonstrated by Gopalan and Katz [9], the bubbles at the closure region of the attached cavity presented hairpin-like structures. To more clearly demonstrate the turbulence flow structures during the shedding process, Figure 9 displays the iso-surface of the Q-criterion (Q = 4 × 10 5 ) for the two hydrofoils. It was found that the vortex shedding structures are closely related with the cavity morphology. The larger scale of the shedding cavity contributes to the larger scale of the U-type vortex and vice versa. During the cavity shedding process, the U-type vortex structures forming on the suction side are governed downstream by the main flow, and gradually collapse. The scale of the U-type vortex structure on the NACA0009 is larger than that of Clark-y. However, the downstream region of the cavity closure on the Clark-y contains more complicated turbulence structures, including numerous and correlative U-type structures. As demonstrated by Gopalan and Katz [9], the bubbles at the closure region of the attached cavity presented hairpin-like structures. To more clearly demonstrate the turbulence flow structures during the shedding process, Figure 9 displays the iso-surface of the Q-criterion (Q = 4 × 10 5 ) for the two hydrofoils. It was found that the vortex shedding structures are closely related with the cavity morphology. The larger scale of the shedding cavity contributes to the larger scale of the U-type vortex and vice versa. During the cavity shedding process, the U-type vortex structures forming on the suction side are governed downstream by the main flow, and gradually collapse. The scale of the U-type vortex structure on the NACA0009 is larger than that of Clark-y. However, the downstream region of the cavity closure on the Clark-y contains more complicated turbulence structures, including numerous and correlative U-type structures. The distributions of Reynolds stresses , , and on the mid-plane of hydrofoils are displayed in Figure 10, where and are the horizontal and vertical turbulent fluctuation velocity components on the mid-plane. As can be seen, the normal and shear Reynolds stresses of high magnitude are located just in the cavitation region. This indicates that the unsteady cavity The distributions of Reynolds stresses u u , u w , and w w on the mid-plane of hydrofoils are displayed in Figure 10, where u and w are the horizontal and vertical turbulent fluctuation velocity components on the mid-plane. As can be seen, the normal and shear Reynolds stresses of high magnitude are located just in the cavitation region. This indicates that the unsteady cavity evolutions greatly affect the turbulelnt features and can induce more turbulence velocity fluctuations. Similar conclusions were also reported by Laberteaux and Ceccio [21] and Ji et al. [24]. In more detail, the normal Reynolds stress u u in Figure 10a is dominant in the attached cavity region, due to the transient changes in the growing and shriking of the attached cavity. The Reynolds stresses u w and w w are predominant in the downstream and closure of the cavity, which is probably influenced by the shedding and the U-type vortex structure formation of the detached cavity. When compared with those of Clark-y, all the Reynolds stress components of NACA0009, as shown in Figure 10, are in higher levels and occupy more areas in the cavitation region. This again forcefully indicates that the time-evolution process of a cavity on NACA0009 shows more violent reentrant flows and a larger scale of U-type shedding structures. Reynolds stresses and are predominant in the downstream and closure of the cavity, which is probably influenced by the shedding and the U-type vortex structure formation of the detached cavity. When compared with those of Clark-y, all the Reynolds stress components of NACA0009, as shown in Figure 10, are in higher levels and occupy more areas in the cavitation region. This again forcefully indicates that the time-evolution process of a cavity on NACA0009 shows more violent reentrant flows and a larger scale of U-type shedding structures. NACA0009 Clark-y (a)

Conclusions
Numerical investigations of the turbulent attached cavitating flows around NACA0009 and Clark-y hydrofoils are respectively carried out by the method of LES, coupled with the Zwart cavitation model under the homogenous flow theory frame. The computational method was evaluated by the published experimental data and presents a fairly good ability to simulate the unsteady cavitating flows. The main conclusions are drawn and summarized below.
First, the U-type shedding characteristics of unsteady cavity are influenced by the interaction of the attached cavity and the re-entrant flow. The similar time evolution process of the turbulent attached cavity around the two hydrofoils presents different shedding traits. For NACA0009, it is easier to form the larger scale of the U-type vortex, and the side re-entrant flow can easily converge at the middle section to impinge the attached cavity together. However, for Clark-y, the break-off cavities tend to merge together again, to cause the final, smaller scale of the U-type shedding cavity. This behavior is probably because the weaker re-entrant flow can hardly cut off the attached cavity with a higher growth speed.

Conclusions
Numerical investigations of the turbulent attached cavitating flows around NACA0009 and Clark-y hydrofoils are respectively carried out by the method of LES, coupled with the Zwart cavitation model under the homogenous flow theory frame. The computational method was evaluated by the published experimental data and presents a fairly good ability to simulate the unsteady cavitating flows. The main conclusions are drawn and summarized below.
First, the U-type shedding characteristics of unsteady cavity are influenced by the interaction of the attached cavity and the re-entrant flow. The similar time evolution process of the turbulent attached cavity around the two hydrofoils presents different shedding traits. For NACA0009, it is easier to form the larger scale of the U-type vortex, and the side re-entrant flow can easily converge at the middle section to impinge the attached cavity together. However, for Clark-y, the break-off cavities tend to merge together again, to cause the final, smaller scale of the U-type shedding cavity. This behavior is probably because the weaker re-entrant flow can hardly cut off the attached cavity with a higher growth speed.
Second, the cavity attached characteristics are related to the pressure distributions. Compared with that of NACA0009, the time-averaged cavity attached on the Clark-y hydrofoil is much larger and is located more downstream from the leading edge, especially at the two sides of hydrofoil. Further analysis verifies this phenomenon is closely relevant to the different pressure distributions on the suction surface of the two hydrofoils.
Third, the turbulent fluctuations and vortex structures are closely related with the cavity shedding dynamics. Due to the more complicated shedding behaviors, the distributions of the Q-criterion iso-surface around the Clark-y hydrofoil shows more disorder with a much smaller scale; also, the time evolution of the lift coefficients shows more fluctuations for the Clark-y hydrofoil. In addition, the contours of the Reynolds stress components indicate that the shedding behaviors with more violent and stronger collision of reentrant flows can induce much higher levels of turbulence velocity fluctuations.
The turbulent cavitating flows around the twisted hydrofoils present some 3D effects, such as the U-type shedding structures and the reentrant flow behaviors. More investigations about these points are probably beneficial for the further study of turbulence-cavitation interactions. Actually, some experimental investigations about the unsteady cavitating flows around these two twisted hydrofoils have been conducted, and we expect the results can be shared in the near future. On the other hand, the verification of the numerical method is always a vital part in the CFD field. Due to the extensive use of the LES in cavitation research recently, it is necessary to discuss its reliability. This work tried to validate the LES method using the typical experimental data. However, some recent comments suggested that the previous verification and validation (V&V) method for LES has some limitations, and they provided a general framework for LES V&V [49,50]. In addition, it has been found that cavitation has a great influence on LES numerical error and modeling error [51]. Therefore, in our future study we are also trying to focus on the V&V investigation of the LES method.

Conflicts of Interest:
The authors declare no conflict of interest.