Abstract
The objective of this paper is to evaluate the predictive capability of three popular transport equation-based cavitation models for the simulations of partial sheet cavitation and unsteady sheet/cloud cavitating flows around a stationary NACA66 hydrofoil. The 2D calculations are performed by solving the Reynolds-averaged Navier-Stokes equation using the CFD solver CFX with the k-ω SST turbulence model. The local compressibility effect is considered using a local density correction for the turbulent eddy viscosity. The calculations are validated with experiments conducted in a cavitation tunnel at the French Naval Academy. The hydrofoil has a fixed angle of attack of ° with a Reynolds number of Re = 750,000 at different cavitation numbers . Without the density modification, over-prediction of the turbulent viscosity near the cavity closure reduces the cavity length and modifies the cavity shedding characteristics. The results show that it is important to capture both the mean and fluctuating values of the hydrodynamic coefficients because (1) the high amplitude of the fluctuations is critical to capturing the extremes of the loads to ensure structural safety and (2) the need to capture the frequency of the fluctuations, to avoid unwanted noise, vibrations, and accelerated fatigue issues.
1. Introduction
In liquid flows, cavitation occurs if the local pressure drops below the saturated vapor pressure, which leads to the formation of vaporous bubble(s) to relieve the negative pressure [1]. Cavitation is commonly observed when a hydraulic machine operates in high speeds or under off-design conditions, and is still a limiting factor in its design and operation. The occurrence of unsteady cavitation in marine propulsion devices such as hydrofoils, propellers, and waterjets can lead to problems such as pressure pulsations, sudden changes in loads, vibrations, noise, and erosion. Moreover, unsteady cavitation is usually the primary physical phenomenon behind performance alterations in a hydraulic machinery [2–4]. For this reason, it is crucial to be able to accurately predict the development and evolution of cavitation, and the resultant impact on the performance.
Phenomenologically, cavitation often involves complex interactions between turbulent flow structures and phase-change dynamics with large variations in fluid density and pressure fluctuations [5]. These physical mechanisms are still not well understood because of the complex, multiscale, multiphase phenomenon. Moreover, there are significant computational challenges associated with the accuracy, stability, efficiency, and robustness of numerical algorithms for the simulation of unsteady, turbulent cavitating flows.
In the numerical modeling of cavitating flows, the selection of cavitation models plays a major role on the prediction of the onset, growth, breakup, and collapse of cavitation bubbles. In recent years, significant efforts have been made in the development of cavitation models; examples of recent review articles can be found in [6–8]. Most cavitation models assume the flow to be homogenous and isothermal and apply either a barotropic equation of state or a transport equation to solve for the variation of the mixture density in the multiphase flow.
In barotropic cavitation models, the local mixture density () is assumed to depend only on the local pressure: . Examples of barotropic cavitation models include those proposed by [9, 10]. A recent experimental finding by [11] has shown that vorticity production is an important aspect of cavitating flows, especially in the cavity closure region. Specifically, this vorticity production is a consequence of the baroclinic torque:
Clearly, if a barotropic equation of state, , is used, the gradients of density and pressure are always parallel, which leads to zero baroclinic torque. Hence, barotropic cavitation models will not be able to properly predict the dynamics of cavitating flows [12], particularly for cases with unsteady cavitation because of the importance of cavitation in vorticity production.
The other popular approach to simulate cavitating flows is transport equation models (TEMs), which solves an additional transport equation for either the mass or volume fraction of vapor, with appropriate source/sink term(s) to regulate the mass transfer between the liquid and vapor phases. Different source/sink terms have been proposed, including the popular models presented in [6, 13–15]. Kubota et al. [13] coupled the Rayleigh-Plesset equation to the URANS equations to solve for the local void fraction based on an assumed bubble radius, and hence it is appropriate for modeling of unsteady cavitation. Merkle et al. [16] and Kunz et al. [14] employed the artificial compressibility method with a special preconditioning formulation to solve the multiphase URANS equations, comprised of the mixture volume, mixture momentum, and constituent volume fraction equations. Singhal et al. [15] developed a “full-cavitation model” based on the phase-change rate derived from a reduced form of the Rayleigh-Plesset equation for bubble dynamics and local flow conditions. The apparent advantage of the Singhal model comes from the convective character of the equation. Senocak and Shyy [17] compared different TEM models to develop an interfacial dynamics-based cavitation model (IDM), which allows direct interpretation of the empirical parameters in existing transport equation-based models. In the TEM approach, the mixture density is a function of the transport process. Consequently, gradients of the density and the pressure are not necessarily parallel, suggesting that TEM can accommodate baroclinic vorticity generation.
The objective of this study is to evaluate the predictive capability of three popular transport-based cavitation models for the simulation of quasisteady and unsteady cavitating flow around a stationary hydrofoil. The aim is to improve the understanding of the influence of multiphase flow on the turbulent closure model, and the interplay between vorticity fields, cavity dynamics, and hydrodynamic performance. The numerical models are first presented. Three popular transport-based cavitation models (named herein as the Kubota model, the Merkle model, and the Singhal model for simplicity) and the - SST turbulence model are summarized, followed by presentation of the density correction model for the turbulent viscosity modification. The foil geometry, fluid mesh, and boundary conditions are presented, followed by a brief summary of the experimental setup of [18] used for validation of the computations. The results are then analyzed for an NACA 66 hydrofoil at a fixed angle of incidence () and Reynolds number (Re = 750,000) for a case with quasisteady partial sheet cavitation (), and a case with unsteady sheet/cloud cavitation (). To study the influence of cavitation model and its interaction with turbulence modeling, results are compared for three different cavitation models and for different turbulent viscosity modifications. Finally, the variation of the predicted cavity length, lift and drag coefficients over a range cavitation numbers are compared with experimental measurements.
2. Transport-Based Cavitation Models
2.1. Conservation of Mass and Momentum
The URANS equations, in their conservative form, for a Newtonian fluid without body forces and heat transfers are presented below along with the mass transport equation in the Cartesian coordinates [19]: where is the mixture density, is the liquid density, is the vapor density, is the vapor fraction, is the liquid fraction, is the velocity, is the pressure, is the mixture laminar viscosity, and are, respectively, the liquid and vapor dynamic viscosities, and is the turbulent viscosity. The subscripts (, , and ) denote the directions of the Cartesian coordinates. The source term and the sink term in (4) represent the condensation and evaporation rates, respectively.
2.2. Kubota Model
The Kubota model [13] assumes a constant nuclei density in the fluid domain. The growth and collapse of the bubble clusters are governed by the simplified Rayleigh-Plesset equation for single-bubble dynamics [20]. The cavitation process is governed by the mass transfer equation given in (4), and the source and sink terms are defined as follows: where is the nucleation volume fraction, is the bubble diameter, is the saturated liquid vapor pressure, and is the local fluid pressure. is the rate constant for vapor generated from the liquid in a region where the local pressure is less than the vapor pressure. Conversely, is the rate constant for reconversion of vapor back into liquid in regions where the local pressure exceeds the vapor pressure. As shown in (7) and (8), both the evaporation and condensation terms are assumed to be proportional to the square root of the difference between the local pressure and vapor pressure because the Kubota model was derived by assuming that the bubble dynamics are governed by the simplified Rayleigh-Plesset equation [20]. In this work, assumed values for the model constants are , m, , and , which are the default values in CFX [5] and are used because of their supposed general applicability.
It should be noted that (7) and (8) are different from the original model proposed by [13], and instead follow the form presented in [21]. Since the evaporation rate is much higher than the condensation rate, different coefficients are needed for condensation and evaporation terms, as shown in (7) and (8). In addition, is replaced by () in the vaporization term to account for the decrease of the nucleation site density with the increase of the vapor volume fraction.
2.3. Merkle Model
Several researchers have adopted the Merkle model proposed by [16] (e.g., see [12]), which has been presented in both the volume fraction form and the mass fraction form. It was derived primarily based on dimensional arguments for large-bubble clusters instead of individual bubbles. Consequently, the source and sink terms for the Merkle model shown in (9) are directly related to the pressure difference, , instead of the square root of the pressure difference as in the Kubota model:
In this work, the empirical factors are set to be and , which follows the constants used by Senocak and Shyy [17] for simulation of quasisteady cavitation around a NACA 66(MOD) hydrofoil. The mean flow time scale is defined as [12, 17].
2.4. Singhal Model
The Singhal mass transfer model proposed by [15] is originally known as the full-cavitation model. This model is also based on the reduced form of the Rayleigh-Plesset equation for bubble dynamics, but empirical constants have been added to take into account the interaction between water and vapor phases. In the derivation of the Singhal model, the bubble size is determined by the surface tension force, and the characteristic velocity, which reflects the effect of the local relative velocity between the liquid and vapor phases. Similar to the Kubota model, both evaporation and condensation terms are assumed to be proportional to the square root of the difference between the local pressure and vapor pressure. The Singhal model assumes that the mass transfer rate is related to the local turbulence kinetic energy, . Following [8], the equations below represent a modified form of the original Singhal model, which consider the effects of noncondensable gas. The evaporation and condensation terms are defined as follows:
In this work, the empirical factors, and (resp., the destruction and the production coefficients for the Singhal model), have the following values: m/s and m/s, and the surface tension coefficient is assumed to be N/m. is the vapor mass fraction, defined as . It should be noted that in the present form of the mass transfer equations, the empirical factors and have to be expressed with m/s units to be consistent with the dimensions of and in (4), which is in units of kg/(m3s1).
2.5. Turbulence Model and Local Compressibility Correction
The numerical results shown in this paper are performed using the commercial CFD code, CFX, to solve the URANS equations. The k- SST (shear stress transport) turbulence model is used, which combines the advantages of the original and models by using the model near the wall, and the model away from the wall [6]. The SST has been found to give good predictions of boundary layer detachment characteristics, and it has been validated for the prediction of subcavitating flow around the NACA66 hydrofoil shown in this paper [22].
It should be noted that the original URANS models were developed for fully incompressible single-phase flows and were not intended for flow problems involving compressible multiphase mixtures. As shown in [23], standard two-equation models (such as and models) overpredict the turbulent eddy viscosity in the rear part of the cavity, which results in overprediction of turbulent stresses, causing the reentrant jet in cavitating flows to lose momentum. Thus, the reentrant jet is not able to cut across the cavity sheet, which significantly modifies the cavity shedding behavior. For partial sheet cavitation in steady flow conditions, inability of the reentrant jet to reach the leading edge will lead to stabilization of the sheet cavity, which prevents the breaking and shedding of the sheet cavity. For unsteady flow conditions, inability of the reentrant jet to fully reach the leading edge will lead to frequent partial shedding of the cavity near the trailing edge, which tends to increase the cavity shedding frequency.
To improve numerical simulations by taking into account the influence of the local compressibility effect of two-phase mixtures on turbulent closure models, Coutier-Delgosha [23] proposed to reduce the mixture turbulent viscosity based on the local liquid volume fraction by substituting in (3) with :
The variation of the modified effective density, , with the vapor volume fraction, , for different values is shown in Figure 1.
In [23], they recommended to better simulate the reentrant jet behavior and vapor shedding process. In [24], the same model was applied for the simulation of supercavitating flow around a hydrofoil, and obtained favorable agreement between numerical and experimental results for . The same modification has also been applied by [25] for the NACA 66 hydrofoil considered in this paper; they used and obtained good agreements with experimental measurements.
3. Experimental Setup and Description
To evaluate the predictive capabilities of the three different cavitation models with different local compressibility corrections, numerical computations are compared with experimental measurements of the NACA 66 hydrofoil conducted at the cavitation tunnel at the Research Institute of the French Naval Academy. The test section has a cross-sectional area of 0.192 m2 and length of 1 m. The inflow velocity ranges between 3 m/s and 15 m/s, and the pressure in the test section ranges between 30 m bar to 3 bars. The tunnel inflow turbulence intensity, defined as / at the inlet of the test section, is about 2%. The foil has a uniform cross-section with an NACA 66 thickness distribution with a maximum thickness-to-chord ratio of 12%, and an NACA camber distribution with a maximum camber-to-chord ratio of 2%. The chord length is m and the span length is m. The hydrofoil is made of stainless steel, and behaved like a 2D rigid hydrofoil even though it was mounted using a cantilevered setup with a small gap (1 mm) between the free end of the hydrofoil and the test section wall. Pressure measurements were carried out using seventeen flush-mounted piezo-resistive transducers (Keller AG 2 MI PAA100-075-010) with a maximum pressure of 10 bars. The transducers locations were aligned along the chord on the suction side of the hydrofoil at mid-span, starting from the foil leading edge at a reduced coordinate of to the trailing edge at , with increments of 0.1 c. Lift and drag were measured using a resistive gauge hydrodynamic balance with a range up to 150 daN in lift and 15 daN in drag. Readers should refer to [18, 26] for additional details about the rigid hydrofoil experimental setup and results. The experimental results presented in this paper are taken from [18] for cases with partial sheet cavitation and from [26] for cases with unsteady sheet/cloud cavitation.
4. Numerical Setup and Description
To demonstrate and validate the numerical model, results are shown for the rigid rectangular NACA 66 hydrofoil described above. All the results shown in this paper corresponds to the hydrofoil fixed at and subject to a nominal free stream velocity of m/s, which yields a moderate-to-high Reynolds number of .
The density and dynamic viscosities of the liquid are taken to be kg/m3 and Pa·s, respectively, which correspond to fresh water at 25°C. The vapor density is kg/m3 and the vapor viscosity is Pa·s. The vapor pressure of water at 25°C is Pa. The 2D fluid domain is shown in Figure 2, which corresponds to the height of the experimental test section at the French Naval Academy. The computational domain has an extent of about 5 c upstream and 10 c downstream of the foil to simulate near-infinite boundary conditions at the inlet and outlet. A no-slip boundary condition is imposed on the hydrofoil surface, and symmetry conditions are imposed on the top and bottom boundaries of the tunnel. The inlet velocity is set to be m/s and the outlet pressure is set to vary according to the cavitation number, defined as , where is the tunnel pressure. A constant turbulent intensity of 2% is set at the inlet boundary and is equal to the experimentally measured turbulent intensity.
All cavitating runs have been initialized with a steady-state calculation using fully wetted models to avoid any vapor fraction at the initial time step. The tunnel pressure is then decreased progressively until the particular cavitation number is reached. Only 2D results are shown in this work. The 2D fluid mesh (shown in Figure 2) is composed of 120,000 elements with 50 structured elements across the foil boundary layer, which is selected to ensure , where is the thickness of the first cell from the foil surface, and is the wall frictional velocity. The regions outside the boundary layer have been discretized with unstructured triangular elements. Mesh refinements are performed at the foil leading edge, trailing edge, and in the wake region. The time integration scheme is a second-order backward Euler algorithm, and the spatial derivatives are computed using a second-order upwind scheme. Mesh convergence and time discretization have been studied in detail for the case of the rigid NACA 66 hydrofoil at a fixed angle of attack of in steady flow; see [22] for more details.
A temporal convergence study is carried out in the case of partial sheet cavitation, at an angle of attack of with . Figure 3 shows temporal convergence of the predicted surface pressure coefficients, , along with the measured values. The error bars represent the maximum and minimum pressure fluctuations. The numerical results were generated using the Merkle model with . The largest time step size used, s, yielded a much stronger pressure gradient in the cavity closure region compared to the measured values. The two finer time step sizes yielded very similar results and compared better with the measured values. Hence, s is chosen for the computations from here on, which gives an average CFL number of .
5. Results
5.1. Partial Sheet Cavitation
The influence of cavitation models and the effect of the turbulent viscosity are first analyzed for the case with quasisteady partial sheet cavitation (). Figure 4 shows the predicted vapor volume fraction and turbulent kinetic energy obtained using the Kubota, Merkle, and Singhal cavitation models with different values for the turbulent viscosity modification (see (11)). The vapor fraction contours for the three different cavitation models with and are, respectively, shown in the first and third rows in Figure 4; a vapor volume fraction of 0 (black) corresponds to pure water, whereas 1 (white) corresponds to pure vapor. The corresponding turbulent kinetic energy contours are shown in the second and fourth rows in Figure 4, where the lowest turbulent kinetic energy is shown in white, and the highest turbulent kinetic energy is shown in black.
(a) Kubota model
(b) Merkle model
(c) Singhal model
For (no turbulent viscosity modification), the three cavitation models yielded very different levels of turbulent kinetic energy in the cavity and at the cavity closure region, which in turn modified the vapor volume fraction and cavity shape, as shown in the first two rows of Figure 4. The Merkle model has the lowest level of turbulent kinetic energy in the vapor region, which leads to a longer sheet cavity compared to the other two models. Compared to the Merkle model, the Kubota model has a higher level of turbulent kinetic energy in the pure vapor region, and a shorter cavity length. It is observed that the production and destruction terms in the Merkle model (9) are related to the pressure difference, instead or the square root of the pressure difference for the Kubota and Singhal models (resp., (7), (8), and (10)). For an equivalent pressure level in subcavitating flow, the direct dependence of the production and destruction terms on the pressure difference leads to a longer cavity for the Merkle model. The Singhal model leads to the highest turbulent kinetic energy in the vapor region. This is because the vapor production and destruction terms for the Singhal model, as shown in (10), are directly related to the square root of the local turbulent kinetic energy, which creates a nonlinear feedback mechanism that increases the pressure gradient in the cavity closure region, and hence tends to destabilize the cavity.
For , the results for all three models in the bottom two rows of Figure 4 show a general decrease in the turbulent kinetic energy in the vapor region, and slightly more unsteady behavior at the cavity trailing edge, compared to . The results show that the density correction to the turbulent viscosity has different degrees of impact on the three cavitation models. The Merkle model predicts similar cavity lengths with and without the density modification, although the decrease of the turbulent kinetic energy allows the reentrant jet at the cavity trailing edge to move upstream, and hence increases the flow unsteadiness and leads to partial shedding of the cavity. The Kubota model shows a noticeable increase in cavity length and more unsteadiness at the cavity trailing edge with the modification. The Singhal model appears to be the most sensitive to the turbulent viscosity modification. Again, this is because of the dependence on the square root of the local turbulent kinetic energy in the production and destruction terms. When , the turbulent kinetic energy is completely suppressed in the cavity region, which allowed the cavity to grow longer but also lead to unsteady partial shedding of the vapor structure near the cavity trailing edge where the pressure gradient is high.
The increase of unsteadiness with increasing values for (11) is due to the development of a reentrant jet as a result of the reduced turbulent eddy viscosity level at the cavity closure. Figure 5 shows the predicted vapor fraction contours and the velocity vectors at the cavity closure region for obtained using the Merkle model with (no modification), , and . With (Figure 5(a)), the turbulent eddy viscosity in the cavity closure region is too high, which prevented the upstream motion of the reentrant jet, and hence the cavity is fully stable. With (Figure 5(b)), the reduction of the turbulent viscosity in the cavity closure region allows the reentrant jet to develop and move upstream, which lead to partial shedding of the cavity. With (Figures 5(c) and 5(d)), the turbulent eddy viscosity is reduced to near zero in the pure vapor region, which allowed the reentrant jet to reach the cavity leading edge. Subsequently, the cavity detached from the leading edge (Figure 5(d)) and then shed downstream as a cloud cavity. This complete breakup of the cavity from the leading edge and unsteady transformation to cloud cavitation obtained using are different from experimental observations by [18], whom reported partial sheet cavitation with behavior similar to that obtained using at .
| (a) , steady cavitation |
| (b) , quasisteady cavitation |
| (c) , unsteady cavitation, |
| (d) , unsteady cavitation, s |
As suggested in [7, 27], the baroclinic torque created by the mixture density and pressure gradients in the cavitating region are responsible for the alteration of the vorticity field, which is evident via the vorticity transport equation: where and are, respectively, the fluid vorticity and velocity vectors. is the laminar kinematic viscosity and is the turbulent kinematic viscosity.
The first term on the left-hand side (LHS) of (12) represents the rate of change of vorticity at a fixed point. The second term on the LHS represents the convection as well as the stretching and tilting of the vorticity due to velocity gradients. The first term on the right-hand side (RHS) of (12) is called the baroclinic torque, which is nonzero only in a variable-density flow field. The baroclinic torque corresponds to the generation of vorticity from unequal acceleration caused by unaligned pressure and density gradients, for example, the generation of a shear layer caused by the faster acceleration of the lighter fluid (vapor) relative to the heavier fluid (water). The second term on the right-hand side of (12) represents the laminar and viscous diffusion of the vorticity.
To examine the influence of cavitation on the production/alteration of vorticity around the case with quasisteady partial sheet cavitation, the strengths of and are compared in Figure 6 for the different cavitation models with different values. As expected, the baroclinic term modifies the vorticity field in regions with high density and pressure gradients, that is, along the liquid-vapor interface and near the cavity closure. The strength of the baroclinic torque increases with , particularly near the cavity closure region because the reduction of the turbulent eddy viscosity tends to promote the unsteady shedding of the cavities. Consequently, the vorticity field is modified, as evident via the contour of the vorticity transport/stretching term. The baroclinic term is found to be the lowest for the Merkle model, which shows a smoother cavity with consequently lower density and pressure gradients. The baroclinic term is higher for the Singhal and Kubota cavitation models, where the magnitude of the baroclinic torque is comparable to the vorticity transport/stretching term near the cavity closure, which will tend to increase the flow unsteadiness, as observed earlier in Figure 4.
| (a) Kubota model, |
| (b) Kubota model, |
| (c) Merkle model, |
| (d) Merkle model, |
| (e) Singhal model, |
| (f) Singhal model, |
Figures 7, 8, and 9 show the comparison of predicted wall pressure coefficients, obtained at , for the Kubota, Merkle, and Singhal cavitation models, respectively, along with the measured values from [18]. On each of the three figures, the mean and fluctuations of the predicted pressure coefficients obtained using , 3, and 10 are shown. The nearly constant pressure region with corresponds to the leading edge, partial sheet cavity on the suction side, which is followed by a sudden pressure jump at the cavity closure due to transition back to subcavitating flow. In general, the pressure fluctuations increased, particularly at the cavity closure region, with increasing values because of the decrease of turbulent eddy viscosity and the additional vorticity generation/modification by the baroclinic torque. The Kubota cavitation model (Figure 7) has very steady pressure behavior without modification (), which showed a cavitating region up to . The results are approximately the same for and 10, which showed a constant pressure region up to , and larger pressure fluctuations at the cavity closure. The Merkle model (Figure 8), on the other hand, showed negligible difference in the predicted pressure coefficients between the different values except for the increase in pressure fluctuations at the cavity closure for . Moreover, the predicted pressure distributions by the Merkle model are in good general agreement with experimental measurements. The Singhal model (Figure 9) shows the largest difference between the different values because of the explicit dependence on the turbulent kinetic energy, which changes with the turbulent eddy viscosity as well as vorticity generation/distortion by the baroclinic torque. Consequently, as increases, the cavity length (constant pressure region) predicted by the Singhal model decreases. A general unsteadiness is observed along the entire suction side pressure distribution. Moreover, the predicted pressure coefficients show greater discrepancy with experimental measurements compared to the Kubota and Merkle models. In general, the results obtained with the Merkle cavitation model with matched better with experiments.
5.2. Unsteady Sheet/Cloud Cavitation
To evaluate the predictability of the three different cavitation models for unsteady sheet/cloud cavitation, results are shown in this section for the case with , Re = 750,000, and . It should be noted that the Kubota model did not converge for this case with s. The authors believe that it is due to the high value of in (7), which makes the calculation unstable. Hence, results are shown in this case for the Merkle and Singhal models only.
At , Re = 750,000, and , unsteady sheet/cloud cavitation was reported in [26]. Figures 10 and 11 show the evolutions of the vapor fraction, turbulent kinetic energy, and vorticity for and with the Merkle and Singhal cavitation models, respectively. Globally, the results show significant reduction in the turbulent kinetic energy in the cavity region, particularly at the cavity closure, and longer cavity with increasing values because of the reduction of the turbulent eddy viscosity. For the Merkle model (Figure 10), it is observed that the level of vorticity is closely related to the vapor fraction contour and the distribution of turbulent kinetic energy. For , when the cavity develops, the decrease of turbulent kinetic energy lead to an increase of vorticity near the wall (Figure 10(b)), which allowed the reentrant jet to form and move upstream () while allowing the cavity to expand to the foil trailing edge; subsequently, the cavity breaks and lifts off at the at the leading edge () to form a large cavitation cloud, which then convects downstream. For the Merkle model, the cavity shedding frequency at is approximately 3.5 Hz, which agrees with the measured value, as shown in Table 1. On the contrary, for , the cavity partially breaks around the midchord (see the vapor fraction in Figure 10(a), at ), which limited the maximum cavity length to , and increased the shedding frequency to about 4.0 Hz, see Table 1. Moreover, because of the presence of cloud cavitation near the foil trailing edge (Figure 10(a) at and ), the overprediction of turbulent viscosity without the density correction () suppresses the formation of the counter-clockwise trailing edge vortex, which changes the vorticity patterns around the foil, as well as the convection of the cloud in the wake. When the modification is applied (), the level of turbulent kinetic energy in the cloud is lower, which allows the counter-clockwise trailing edge vortex to form and interact with the cloud cavity (Figure 10(b) at and ).
| (a) Merkle model, |
| (b) Merkle model |
| (a) Singhal model, |
| (b) Singhal model, |
As shown in Figure 11, the Singhal cavitation model is more sensitive to the turbulent viscosity modification. Without the modification (), the level of turbulent kinetic energy is very high near the cavity closure (see Figure 11(a)), which limited the length of the sheet cavity to approximately (), and the cavity partially breaks at (), followed by successive shedding of smaller vapor structures (). The cavity appears to be more unstable and transient than that predicted by the Merkle cavitation model, as shown by the scattered vorticity distribution near the wall (see Figure 11(a)). When the modification is applied (), the Singhal model also predicted an increase in the cavity length up to , although general cavity behavior is much more unstable than that predicted by the Merkle model. The cavity breaks around . However, it does not shed as one big cavity; instead, the aft end of the cavity is broken into multiple little cloud cavities, as shown in Figure 11(b). Therefore, the shedding frequencies predicted by the Singhal model (about Hz for and Hz for ) are much higher than those by the Merkle model and did not agree with experimental measurements, as summarized in Table 1. Again, the turbulent viscosity modification allowed the trailing edge vortex to develop earlier, which interacted with the cloud cavity as it sheds into the wake.
Figure 12 shows the baroclinic torque and the vorticity transport/stretching terms, at for both the Merkle and Singhal models. For both models, the baroclinic torque concentrates along the regions where the cavity breaks, although the magnitudes seem to be small compared to the vorticity transport/stretching term. Nevertheless, the influence of the cavity on the vorticity distribution is clearly evident by the significant modification of the vorticity transport/stretching term (Figure 12) and vorticity contours (Figures 10 and 11). Hence, the interdependency between cavitation generation/production and vorticity generation/modification due to the baroclinic torque is not negligible, particularly in terms of its effect on the cavity stability and shedding frequency.
| (a) Merkle model, |
| (b) Merkle model, |
| (c) Singhal model, |
| (d) Singhal model, |
Based on available experimental data presented in [26], an analysis of suction side loading can be calculated by summing the pressure coefficients along the suction side. The approximate suction side lift coefficient can be written as: where is the pressure coefficient at location , and is the nondimensional distance between two consecutive pressure transducers. A total of 10 pressure transducers were placed along the suction side at the mid-span of the foil [22]. For comparison purposes, the same procedure is applied to calculate the numerically derived suction side lift coefficients. It should be noted that since the flow is mostly attached along the pressure side, the changes to the total lift coefficient should be dominated by the suction side dynamics, which is represented by the suction side lift coefficient.
As shown in Figure 13, a fair comparison is observed between the experimental measurements and numerical computations obtained using the Merkle model with . The peaks generally correspond to the time instance when the maximum cavity length is observed ( in Figure 10(b)), and the unloading phase corresponds to the shedding of the cavity. The results show that the rate of the unloading phase predicted by the Merkle model is faster than the measured value. This difference could be due to inaccurate prediction of the convection/destruction of vapor in the cloud, which is by definition very transient and sensitive to the production and destruction rate of vapor, and the interaction between the cavity and the vorticity field through the baroclinic torque. Nevertheless, the shedding frequency predicted by the Merkle model with matched well with the measured frequency of Hz. As shown in Table 1, the cavity shedding frequency is higher for because the cavity is shorter and undergoes partial shedding at its trailing edge. Hence, it is important to take into account the influence of mixture density variation in the calculation of the turbulent viscosity in and near the cavity.
As shown in Figure 14, the Singhal model with shows a much higher shedding frequency, Hz, than the measured frequency of 3.5 Hz. The high frequency predicted by the Singhal model is because of the continuous presence of vapor along the hydrofoil surface (see Figure 11(b)), which significantly perturbs the wall pressure distributions and leads to the generation of vorticity through the baroclinic torque (see Figure 12(d)). As a consequence, the predicted suction side lift coefficient shows a much higher shedding frequency and much lower amplitude compared to the measured values. The results shown in Table 1 confirms that the shedding frequency is higher for the Singhal model if , without the turbulent viscosity modification.
It should be noted that the importance of capturing and comparing both the mean and fluctuating values is evident in Figures 13 and 14. As shown in the figures, although the mean suction side lift coefficient predicted by both the Merkle and Singhal models are approximately the same, and agreed well with experimental measurements, the amplitude and frequency of the fluctuations predicted by the Singhal model are very different from the Merkle model, which compared much better with experimental measurements.
Figure 15 shows the comparisons between the predicted vapor fraction contours obtained using the Merkle cavitation model with and photographs from the experiments presented in [26] for , Re = 750,000, and . The white vapor fraction represents the pure vapor and the black represents the pure water. The predicted flow streamlines are also shown to facilitate visualization of the flow dynamics. Similar to the comparison of the suction side loading, the numerical predictions obtained using the Merkle model with are in good agreement with the experiments. The cavity expanded up to the foil trailing edge at ( to ), followed by breakup near the foil leading edge due to the reentrant jet at . An interaction between the bubbly cavity mixture and the trailing edge vortex is observed at and , followed by complete convection of the cloud cavity into the wake at .
5.3. Comparisons of the Measured and Predicted Hydrodynamic Load Coefficients and Cavity Lengths
Figures 16, 17, and 18 show the measured and predicted normalized cavity lengths , lift coefficient , and drag coefficient over a range of cavitation numbers . The numerical results for all three models are obtained using , which compared the best with the experimental measurements, as shown in the previous section.
Figure 16 shows the predicted maximum and minimum cavity lengths, along with the measured maximum cavity lengths. Cavitation inception is observed at for both computations and experiments. For , the minimum and maximum values are equal, which suggest stable cavitation. It should be noted that all three cavitation models give very similar results for stable cavitation, and all agree well with the experimental measurements (star symbols). For , the cavity becomes unstable, as indicated by the diverging maximum and minimum values of the cavity length, and differences could be observed between the three cavitation models. The results suggest that the cavity breaks directly at the leading edge for the Kubota model, as indicated by the near-zero value of the minimum length represented by the open circles in Figure 16 for . However, experimental observations suggest the existence of partial sheet cavitation between , which is characterized by a partial break of the cavity with a stable vapor region at the cavity leading edge [18]. As shown in Figure 16, the Singhal and Merkle models agree well with the experimental observations, but not the Kubota model. As the cavitation number further decreases, the break of the cavity moves closer to the foil leading edge. At , unsteady sheet/cloud cavitation occurs, as shown in the numerical and experimental results reported in the previous section. It should be noted that the measured maximum cavity lengths agree well with the numerical values for all three cavitation models considered. However, as evident via the results shown in the previous section and in Table 1, even if the minimal and maximal values of the cavity length agree well between the Singhal and the Merkle cavitation models, the shedding frequencies are very different. It should be noted that the Kubota model did not converge due to excessive fluctuations for , and hence results are only shown for in Figures 16–18.
Figure 17 compares the measured and predicted lift coefficients (, where is the lift force) for different cavitation numbers. The lift coefficient increases slightly with the development of small, leading edge, and partial sheet cavitation (). The increase in lift is attributed to the suction effect created by the small leading edge cavitation bubble caused by the high pressure gradient. As the cavity grows with decreasing cavitation number (i.e., , ), the lift coefficient decreases. Both the slight increase in the lift coefficient for and subsequent drop in lift coefficient for are captured by the Merkle and Singhal models. The numerical predictions show large fluctuations in the lift coefficients for , as indicated by the error bars, due to the unsteady cavity dynamics. The numerical results show good agreement with the experimental measurements for , in the quasisteady cavitation and fully wetted flow regimes. The predicted lift coefficients by all three models are slightly less than the measured values for , in the partial sheet cavitation flow regime. In that range, the differences between the three cavitation models are relatively small, with the Merkle model showing the most favorable, and the Singhal model showing the least favorable comparisons. The Kubota model also predicts much higher fluctuations compared to the other two models in the unsteady cavitation flow regime, and the fluctuations grow with decreasing cavitation number.
Comparisons of the predicted and measured drag coefficients (, where is the drag force) are shown in Figure 18. All three cavitation models underpredicted the drag coefficients, but all three trends agree well with the experiments. It should be noted that the underprediction of the drag may be attributed to 3D effects, where downwash and wall boundary layer effects tend to increase the drag. The results show that the development of cavitation significantly increases the drag, particularly for unsteady cavitation. The combination of decreases in lift and increases in drag leads to a decrease in the efficiency, which is a classical issue for cavitation. It should be noted that the magnitudes of the fluctuations in the lift and drag coefficients are higher for the Merkle model than the Singhal model at , which is consistent with the results shown in Figures 13 and 14.
6. Conclusions
The predictive capability of three popular cavitation models is investigated for the prediction of cavitation flow around a NACA 66 hydrofoil at a fixed angle of incidence () at , for a case with partial sheet cavitation (), and a case with unsteady sheet/cloud cavitation (). Comparisons of the predicted and measured cavity lengths and hydrodynamic coefficients over a range of cavitation numbers are also shown. In addition to the application of different cavitation models, different levels of density correction to the turbulent viscosity are applied to investigate the interaction between cavitation and turbulent flows. The aim of this density correction is to take into account the local compressibility effects in the cavitating region, especially near the cavity closure, where standard turbulence models tend to overpredict the turbulent viscosity, which will modify the cavity shedding dynamics.
The general conclusions of the paper are as follows:
(i) the density modification helps to reduce the turbulent viscosity in the cavitating region, which decreases the turbulent kinetic energy in the cavity, as well as around the cavity closure. With the reduction of the turbulent viscosity, the reentrant jet is allowed to develop underneath the cavity and move upstream while allowing the cavity to expand toward the foil trailing edge. Consequently, the cavities tend to be slightly longer and more unsteady with the density modification. This lead to a better comparison with experimental measurements and observations for both the quasisteady cavitation case and the unsteady sheet/cloud cavitation case;
(ii) the density correction to the turbulent viscosity has varying degrees of impact on the different cavitation models. The Singhal model is the most sensitive to the turbulent viscosity modification, which already overpredicted the general unsteadiness due to the explicit dependence of the source and sink terms on the turbulent kinetic energy, the addition of the density correction will further increase the overall unsteadiness, and hence lead to greater discrepancies compared to experimental measurements. On the other hand, the Merkle model seems to be the least sensitive to the turbulent viscosity modification, and tends to be more stable than the Singhal and Kubota models, particularly for unsteady cavitating flows. Addition of the density correction for the Merkle model helps to avoid overprediction of the turbulent eddy viscosity, which allows the reentrant jet to properly develop, and hence leads to improved comparisons between predictions and measurements;
(iii) the baroclinic torque is responsible for the generation of additional vorticity due to shear flow created by the unequal acceleration between the lighter vapor and heavier liquid phases. Hence, in regions with high, unaligned density and pressure gradients, for example, along the liquid-vapor interface and the cavity closure region, the vorticity field is modified by the cavitation, which in turn changes the cavitation dynamics. The interdependency between the cavity and vorticity dynamics through the baroclinic torque is not negligible, particularly for unsteady cavitating flows. The magnitude of the baroclinic torque tends to increase with the power factor in the density correction for the turbulent viscosity (11), which in turn increases the flow unsteadiness, particularly near the cavity closure. New detailed experimental measurements of the mean and fluctuating velocity fields in the cavitating region are needed to investigate the validity of the density correction, and the importance of the baroclinic torque on the cavitating response of hydrofoils;
(iv) the results show that while predictions by all three cavitation models compare well with experimental measurements in the quasisteady partial sheet cavitation range, the predicted cavity shedding frequencies and load fluctuations predicted by the three cavitation models are very different in the unsteady sheet/cloud cavitation range. Moreover, the numerical stabilities of the cavitation models are very sensitive to the dependency of the condensation and evaporation rates on the power to the pressure difference term (1.0 in the Merkel model and 0.5 in the Kubota and Singhal models), which is probably attributed to the interdependency between the cavitation and vorticity dynamics through the baroclinic torque. The Singhal model is different from the Kubota and Merkel models in that it explicitly depends on the square root of the turbulent kinetic energy (see (10)), which increases the coupling effect between the vorticity distribution in the turbulent flow field and the cavity dynamics, and hence tends to destabilize the cavity. Consequently, the predicted cavity shedding frequencies was much higher than the measured values, and the amplitude of the load fluctuation was much lower than the measured values. Physically, as noted in [28], transition to turbulence will sweep away attached cavities. Hence, explicit dependence of the sink (vaporization) term on the turbulent kinetic energy does not agree with experimental observations presented in [28], nor with experimental measurements shown in this paper. Nevertheless, additional detailed experimental measurements of the mean and fluctuating velocity fields, as well as the time histories of the vapor fraction distribution in the cavity is needed to truly quantify the validity of the various cavitation models;
(v) in general, the Merkel model with a density correction of shows the best agreement with experimental measurements and observations for both steady and unsteady cavitating flows, particularly in terms of the cavity shedding frequency. The predicted shedding for unsteady sheet/cloud cavitation obtained using the Merkle model with matches with the measured value of 3.5 Hz. The Kubota model is not able to converge for the case of unsteady sheet/cloud cavitation, and the Singhal model predicts much higher cavity shedding frequency and much lower amplitude of fluctuations compared to the measured values;
(vi) the analysis showed that the destruction and production coefficients in the Singhal model are not dimensionless parameters, but instead should have units of velocity to be consistent with the units of the mass transport equation;
(vii) the results demonstrate the importance of capturing both the mean and fluctuating values of the hydrodynamic coefficients and cavity lengths for the case of unsteady sheet/cloud cavitation because (1) the fluctuations are in the same order as the mean values, and hence knowledge of the fluctuations is critical to capturing the extremes of the loads to ensure structural safety; (2) the need to capture the frequency of the fluctuations, which are critical to avoid unwanted noise, vibrations, and accelerated fatigue issues. Moreover, it is also critical to compare both the amplitude and frequencies of fluctuations in addition to the mean values when evaluating the accuracy of numerical models.
This paper helps to improve the understanding of the predictive capability of the three popular cavitation models. The impact of the density correction term on the different cavitation models and its influence on the cavitation dynamics, vorticity field, and the force coefficients have been investigated. Additional research is needed to better understand the influence of the local compressibility and baroclinic torque on the cavitation and wake dynamics, as well as turbulent vorticity fields, particularly for cases with spatially/temporally varying inflow, and for cases with rigid and/or elastic body motions. Such research is important because accurate prediction of the cavity and load fluctuations is critical when analyzing the hydroelastic stability, fatigue, noise, vibration, and erosion characteristics of hydraulic machineries.
Acknowledgments
The authors gratefully knowledge the Office of Naval Research (ONR) and Dr. Ki-Han Kim (program manager) for their financial support through Grant no. N0014-09-1-1204.