Abstract
Two cavitation models with thermodynamic effects were established based on the Rayleigh-Plesset equation to predict accurately the cavitation characteristics in the high-temperature fluid. The evaporation and the condensation coefficient of the cavitation model were corrected. The cavitation flow of NACA0015 airfoil was calculated using the modified cavitation model, where the influence of the thermodynamic effects of airfoil cavitation was analyzed. The result showed that the pressure coefficient distribution and the bubble volume fraction simulated have the same tendency of Zwart-Gerber-Belamri model’s result. According to the experimental data, the two models provide more accurate results. At the room temperature, the values of obtained by the two improved models are approximately equal. The difference between the two models’ results increases gradually with the temperature increasing, but it is still small. The simulation results are consistent with the experimental data when the evaporation coefficient is 10 and 1. When the evaporation coefficient is 1, the bubble growth is inhibited, the volume fraction becomes lower, and the cavitation area becomes flat. As the temperature increases, the cavitation area and the bubble volume fraction at airfoil front edge become larger, showing that the temperature plays a “catalytic” role in the cavitation process.
1. Introduction
Cavitation belongs to multiphase flow category. It is a phase transition phenomenon when the internal fluid pressure of the pump is lower than the saturation pressure [1–3]. At present, two main types of numerical models (homogeneous and the nonhomogeneous model) can solve the multiphase flow. The thermal dynamic mechanism plays an important role in the cavitation process; it has already attracted the attention of researchers. In the theoretical analysis aspect, the rule of bubble growth in water was deduced [4]. A method to model the thermal dynamical process of cavitation was proposed using thermal conduction in bubble surface. Low-temperature medium was applied to compare the temperature drops in cavitation flows where the simulation results and experimental data are basically the same [5]. Considering the thermal effect, the cavitation phenomenon in an inducer was studied. As the temperature was increasing, the length of the cavity is reduced while cavitation performance was improved [6, 7]. The influence of the dimensionless thermodynamic coefficient on the critical cavitation number in inducer and its cavitation states under different temperature was studied. When the thermodynamic coefficient is lower than 0.54, the critical cavitation number decreases along with the increase of thermodynamic coefficient; and vice versa: when it is higher than 0.54, the critical cavitation number is not related to the thermodynamic coefficient [8].
The cavitation experiments at different temperatures in an inducer and NACA0015, respectively, were carried out. When the temperature was increased under the same cavitation number, the pressure pulsation was reduced and the cavity thickness and length on the hydrofoil were increased [9, 10]. Rotating cavitation at three different temperatures and the distribution of cavity lengths and temperatures were investigated [11]. The cavitation flow characteristics in fluorinated ketone of NACA0015 hydrofoil at different speeds and different attack angles were illustrated [12]. For the cavitation in low temperature, the dynamic characteristics of liquid carbon dioxide cavitation near the critical point were studied. The vapor kinetic equation considering the thermodynamic effect and quantitatively predicting the cavitation of liquid carbon dioxide near the critical point was also established. It was found that the thermodynamic boundary has an important effect on the bubble dynamics [13]. A cavitation model was set up to capture the heat transfer process during the collapse of the cavity bubbles, which provides a theoretical basis for calculating the local temperature, pressure, and heat transfer parameters during the collapse of the cavitation [14]. Also, a cavitation model considering the thermodynamic effects based on the Zwart cavitation model was set up and applied for simulating the cavitation state of the hydrothermal manifold at different temperatures where the experiments were carried out for calibrating parameters [15]. A heat transfer model of a bubble based on the Fourier heat transfer law was used to calculate the effect of thermodynamics on the NACA0015 airfoil; then the experiments were carried out to verify the accuracy [16]. Based on the Plesset-Zwick formula and the effect of temperature on bubble growth, a cavitation model was created to introduce the concept of the thermal boundary of the bubbles, numerical simulation about cavitation in NACA0015 hydrofoil, and the subminiature pump [17]. Also, a cavitation model based on Singhal full cavitation model including thermals effect was investigated. The model was verified by tested pressure coefficient of NACA66 hydrofoil in hot water when [18]. The vapor pressure based on Singhal cavitation model and the thermal effect on cavitation under high-temperature conditions was studied [19]. A secondary development of technology was applied to introduce modified Kunz and Merkle cavitation model and physical parameters of nitrogen vary with the temperature of the flow field. The cavitation flow in symmetrical rotors and hydrofoils nitrogen was simulated [20, 21]. The mixed fluid model based on Rayleigh-Plesset fluid equations and thermals theory was applied to analyze the rotating cavitation phenomena in the two-dimensional cascades of inducer. Briefly, the research on cavitation under thermals effect has made some achievements, but the accuracy and applicability of the model have not been improved yet, and a large deviation still exists between the numerical and experimental results [22].
In this study, based on the Rayleigh-Plesset equation, the influence of thermodynamic effects is introduced in the evaporation and condensation coefficient of the cavitation model. Two cavitation models considering thermodynamic effects are proposed and were added to the calculation by commercial software CFX secondary development. The cavitation flow of the hydrofoil airfoil and the NACA0015 airfoil was calculated by adapting the Zwart-Gerber-Belamri cavitation model and the modified cavitation model. By contrast and analysis between the simulation result and experimental data, the applicability of the modified cavitation model in high-temperature cavitation flow was validated.
2. Geometric Model and Meshing
The geometry simulated was consistent with the model in the experiment mentioned before [9]. The specific parameters of the NACA0015 airfoil are chord length 115 mm, wingspan 80 mm, and adjustable attack angle of 4°, 5°, 6°, and 8°. During the test, three pressure holes were evenly arranged at 100 mm in front and 140 mm at the rear of the airfoil. At the same time, 10 pressure holes were arranged on the airfoil suction surface while two pressure holes were arranged on the airfoil pressure surface, as shown in Figure 1.
In the numerical calculation, the average + of the airfoil suction plane and pressure surface are 21.2 and 18.5, respectively. The whole, the head, and the tail grids are shown in Figure 2. To exclude the effect of the grid number on the simulation accuracy, grid independence study is carried out at the beginning of the numerical investigation. The flow field around the NACA0015 hydrofoil was simulated using meshes with five different grid numbers, and the lift force calculated from different mesh is shown in Figure 3. When the grid number is larger than 1 million the lift force is almost unchanged, and the error is still below 0.3%. Consequently, the grid number for the later investigation is selected as 1 million.
The boundary conditions are consistent with those in the experiment. The inlet boundary is set as velocity inlet, and the outlet boundary adapted static pressure under forced outflow. The airfoil suction and pressure surface are set as nonslipping solid wall and the roughness of the airfoil surface was ignored. The convergence accuracy was set to , where the cavitation number and the pressure coefficient are defined as follows:where is the reference pressure, is the reference velocity, and is the vaporization pressure of water.
3. Governing Equations and Simulation Method
3.1. Governing Equations
The homogeneous equilibrium flow model assumes that the medium is a mixture with the same velocity and pressure field. The governing equations are mainly set according to the continuity equation, the momentum equation, and the energy equation.where , is the time, is the mixing medium velocity, is the density of the mixed medium, is the pressure of the mixed medium, is the laminar viscosity coefficient, is the turbulent viscosity coefficient, is the specific heat capacity, is the temperature, is the mass fraction of water vapor, is the latent heat of vaporization, is the liquid volume fraction, and is the enthalpy of the substance.
3.2. The Selection of Turbulence Model
For the turbulence closure to solve in (3) and (4), a turbulent model is necessarily solved in combination with the governing equations. According to the experimental investigations, turbulence is found to have a significant effect on the cavitation intensity [23, 24]. Thus, to make a better prediction for cavitating flows, it is important to simulate the flow field accurately and deal with the turbulent terms appropriately. The common turbulent models usually approximate all the length scales, which may create higher eddy viscosity field. Due to the numerical dissipation, some large periodic turbulence flows have the possibility of being ignored [25]. Consequently, filter functions are proposed to limit the eddy viscosity to a reasonable degree [26, 27].
In the current research, the Rayleigh-Plesset cavitation model is solved in combination with four common turbulence models (standard k-ε turbulence model, standard k-ω turbulence model, RNG k-ε turbulence model, and SST turbulence model) and also two different filter turbulent models (k-ε FBM turbulence model and RNG k-ε FBM turbulence model). Simulation results are compared to select the suitable turbulent model for predicting the cavitation flow field. NACA0015 airfoil at 5° attack angle and 1.5 cavitation number at 25°C is simulated.
As shown in Figure 4, the pressure coefficient distributions of the airfoil obtained by numerical calculation using typical turbulence model and the filter turbulence models have the same tendency. There are two mutations between and 0.4. When = 0.4~1.0, the pressure coefficients simulated by different turbulence models are not so different and well agree with the experimental data. In the airfoil head cavitation area, where = 0~0.25, there is a big difference between the predicted value and tested values. The cavity length simulated using the filtered turbulence model is shorter than the typical turbulence model and the cavity length of the k-ω turbulence model is the longest one among the typical turbulence models, and the cavity length of the SST turbulence model is the shortest one. By comparing the simulation results using different turbulence models with the experiment data, it can be judged that the calculated result using RNG k-ε turbulence model is the closest to the experimental value. In the position where chord length ratio = 0.25~0.4, the value of calculated using the filter models is lower than the result of the k-ε model and RNG k-ε model. It is indicated that the simulation result is sensitive to the filter function of turbulent viscosity and more robust filter function or more suitable filter scale is necessary for the simulation of turbulent cavitating flows. The result of the common turbulence model without filter functions is similar, while the k-ε based models are closer to the experimental data than the k-ω based models, and RNG k-ε model is more accurate than standard k-ε models as it replaced the constant with more detailed function . Therefore, the RNG k-ε turbulence model is chosen in this study.
3.3. The Modified Cavitation Models
The common cavitation models are mostly based on the bubble growth equation proposed by Rayleigh-Plesset (R-P) [28]. The R-P model provides the basic equations for controlling the vapor generation and condensation and describes the growth of bubbles in the liquid; the details are as follows:where is the bubble radius, is the bubble pressure, is the liquid pressure, is the surface tension coefficient, and is liquid kinematic viscosity.
Assuming that the bubble growth process is a spherical model, the bubble volume change rate is as follows:
The mass change rate of the bubble is
If there are bubbles per unit volume, the volume fraction can be expressed as
Then the total interphase mass transfer rate (evaporation process) per unit is
Similarly, the condensation coefficient is deduced as follows:
The thermodynamic parameter is a function of temperature, the Taylor series of vaporization pressure is expanded, and the higher order terms (including the quadratic term), the surface tension term, and the viscosity term in the R-P (5) are ignored. The reduced model is as follows:where is the pressure difference term and II is the thermodynamic term.
If the thermodynamic effect is not taken into account, the pressure difference between liquid phase and the vapor phase drives the interphase mass transfer; then ; that is, . The bubble radius growth rate without thermodynamic effect is reduced as
However, in the variable temperature medium, the latent heat of vaporization absorbed during the vaporization process forces the surrounding liquid to form a temperature gradient outside the bubble, which results in a temperature difference between the bubble and the outer liquid. The temperature difference will affect the growth of the bubble. Therefore, in order to predict the cavitation phenomenon more accurately, the existing cavitation model must be properly modified by adding the source term considering the thermodynamic effect.
Method 1. According to the function of pressure and temperature, the relationship of is given directly [29]:Then dR/dt was obtained:In the formula, .
Method 2. was obtained by using the Clapeyron formula:Then dR/dt was obtained:Equations (14) and (16) are put into (9) and (10), respectively, to obtain two kinds of the cavitation models considering thermodynamic effects (Thermal Cavitation Model, TCM). These two cavitation models are named TCM1 and TCM2, respectively. The specific expression is as follows:where is the factors of thermodynamic effects and the two forms of are and , respectively.
The results show that turbulent kinetic energy may have an important effect on cavitation. The above models use the modified method to reduce the effect of turbulent kinetic energy on cavitation [19]. The formula is as follows:where is the local saturated vapor pressure and is the local turbulent kinetic energy.
4. Validation and Modification of Cavitation Model considering Thermodynamic Effects
4.1. Comparison and Analysis of TCM1 and TCM2
The surface pressure coefficient distributions of the NACA0015 airfoil predicted by the three cavitation models (the Default cavitation model, TCM1, and TCM2) at 25°C are shown in Figure 5. It can be seen that, with different cavitation models, the position of important turning point on the suction surface of the airfoil is almost the same as that in the experimental result, which indicates that the thermodynamic effect of cavitation at this temperature is not prominent. Therefore, the numerical simulation of cavitation under room temperature can ignore the effect of thermodynamic effects. However, on the airfoil leading edge, the numerical simulation results predicted by the three cavitation models are quite different from the test experiment data, and the error of the pressure distribution coefficient of the suction surface at the chord length ratio 0.1 is about 15%. The pressure coefficient distribution and bubble volume fraction distribution map are very similar in two kinds of modified cavitation models. Also, the region and the location where bubble occurred in these two models are almost the same.
(a) Comparison of NACA0015
(b) Comparison of vapor volume fractions
Comparing (13) and (15), it is found that the main difference between TCM1 and TCM2 is in the definition of . TCM1 is obtained based on the derivative of the continuous function , and the function describing relationship between vaporization pressure and the temperature can be deduced by polynomial fitting. In TCM2, the term is obtained based on the Clapeyron equation approximation.
The specific values of the two improved cavitation models TCM1 and TCM2 at different temperatures are shown in Table 1. The values of obtained in room temperature of two improved models are basically the same. With the increase of temperature, the difference between the two shows a gradually increasing trend, but the overall gap is not large.
It can be seen from the above analysis that the improved thermodynamic cavitation models TCM1 and TCM2 are of the same essence and only different in the expression form. The effect of cavitation thermodynamics at room temperature is negligible. When the temperature is under 100°C, the difference in values between TCM1 and TCM2 is small. Therefore, at 70°C, the cavitation model TCM2 is selected for verification.
4.2. Modification and Validation of TCM2
Figure 6 shows the comparison chart of pressure coefficient distribution on NACA0015 airfoil surface at 70°C of modified cavitation model TCM2 and the test value. It can be seen from the figure that the simulated values in the airfoil rear area are better identified with the test value, but in the airfoil head area where chord ratio ranges from 0.1 to 0.4 the difference is large. Therefore, the cavitation model TCM2 still needs to be corrected.
Evaporation and condensation coefficients both have an effect on cavitation, so the evaporation and condensation coefficients need to be corrected after the cavitation model is modified. The default evaporation and condensation coefficients are 50 and 0.01. According to experience, the other two sets of coefficients are defined as 10, 0.002 and 1, 0.0002, respectively. Figure 7 shows the comparison chart of the airfoil pressure distribution coefficient of cavitation model TCM2 at 70°C under three evaporation and condensation coefficients. It can be seen from the figure that when the evaporation coefficient is 10 and 1, the simulation value agrees well with the experimental value. With the decrease of evaporation and condensation coefficient, the pressure coefficient distribution curve becomes smooth. However, when the evaporation coefficient is 1, several turning points of some important chord ratios are not shown, the bubble growth is suppressed, the bubble volume fraction becomes lower, and the shape of the bubble area becomes flat. Therefore, the cavitation model TCM2 is chosen in the following cavitation simulation, and the evaporation and condensation coefficient are set as 10 and 0.002, respectively, which is the same as in [17].
5. Cavitation Thermodynamic Effects in Hydrofoil
Based on the cavitation model TCM2 with evaporation and condensation coefficients of 10 and 0.002, respectively, the cavitation conditions in hydrofoil airfoil and thermodynamic effects in NACA0015 airfoil were simulated, respectively.
5.1. The Cavitation in Hydrofoil at 25°C and 50°C
The hydrofoil airfoil is first simulated, and its boundary conditions are set as follows: the inlet velocity is 16.91 m/s and the outlet static pressure is 51957 Pa. Airfoil working conditions and grid are shown in Figure 8; the airfoil grid average y+ = 89.2.
The interphase mass exchange rate distribution diagram of airfoils at 25°C and 50°C is shown in Figure 9. The interphase mass exchange rate (vapor to liquid) reflects the speed of interphase mass exchange, which indirectly reflects the generation and development of bubbles. It can be seen from Figure 9 that, at different temperatures, the regions where bubble is generated and disappears are similar to each other. The bubble is mainly generated downstream of the middle of airfoil chord on the suction surface and disappears after it reaches the high-pressure zone. However, compared to the interphase mass exchange rate distribution at 25°C, the bubble generation region at 50°C is wider, and the area of the vaporization region and the bubble volume fraction are larger. At 50°C, the process of the bubble to become smaller or disappear is more obvious and the area of collapse region becomes larger, which both will cause important consequence. It can be seen that the effect of thermodynamics on the cavitation process of hydrofoil airfoil is obvious, and the higher the temperature, the more obvious the development trend of cavitation.
5.2. The Cavitation in NACA0015 Hydrofoil at 25°C, 50°C, and 70°C
Numerical simulations were carried out at 25°C, 50°C, and 70°C, respectively, based on the NACA0015 airfoil, at 8-degree attack angle and 2.5 cavitation number. The distribution of bubble volume fraction and the isosurface of 0.5 is shown in Figure 10. It can be seen from Figure 10 that the degree of cavitation in NACA0015 airfoil at different temperatures is different and cavitation region occurs mainly in the wing head where chord ratio ranges from 0.01 to 0.25. With the increase of temperature, the bubble volume fraction in the cavitation region becomes higher and the cavitation thermodynamic effect becomes more obvious, showing that the temperature plays a “catalytic” role in the cavitation process.
(a) Bubble volume fraction
(b) Isosurface with 0.5 bubble volume fraction
Figure 11 shows the pressure coefficient distribution of airfoil surface at 25°C, 50°C, and 70°C with 8-degree attack angle. It can be seen from the figure that the positions of the important turning points on three curves are similar. The pressure distribution curve can be divided into three regions. The first region is the cavitation generation area, where the pressure coefficient value of suction surface () decreases slowly along the chord. The second region corresponds to the cavitation weakening process, and the value decreases rapidly in this region. In the third region, the cavitation effect almost ends and the pressure decreases slowly along the chord length ratio. Compared with the calculation results at three different temperatures, it can be found that the three regions of the pressure coefficient curve have obvious demarcation points at 70°C, and these demarcation points correspond to the position where the chord length ratio is and 0.18, respectively, while at 25°C and 50°C, that curve drops comparatively gently and there is no obvious demarcation point to distinguish three different regions. At 70°C, the first region corresponds to the area where the chord length ratio is about = 0~0.1, and the length of the region where cavitation occurred is greater than that in 25°C and 50°C situation. In the area where chord length ratio is in the range from 0 to 0.07, the pressure coefficients at three temperatures are relatively close, while in the area where chord length ratio is in the range from 0.07 to 0.18, the pressure coefficients at three temperatures show a large difference. When the chord length ratio ranges from 0.07 to 0.14, the lowest value of the pressure in the airfoil surface occurs at 70°C while the highest value occurs at 25°C. Then the pressure at 70°C rises rapidly, while the pressure at 25°C rises slowly. When the chord length ratio reaches 0.14–0.18, the highest value of the pressure in the airfoil surface occurs at 70°C while the lowest value occurs at 25°C. When the chord length ratio is larger than 0.18, the pressure values in the airfoil surface at three different temperatures all increase gradually and the difference between the three temperatures decreases.
Figure 12 shows the distribution of vacuolar volume fraction along the chord length ratio in the suction surface at three temperatures. It can be seen that when the chord length ratio ranges from 0 to 0.14, the highest value of vapor volume fraction occurs at 70°C while the smallest value occurs at 25°C. But when the vapor volume fraction decreases downstream along the chord, the higher the temperature, the faster the decrease speed, which indicates that, in the area where cavitation disappears, higher temperature promotes condensation process.
In summary, the bubbles in the growth process will transfer downstream along with the flow field; in the area where chord length ratio is in the range from 0.07 to 0.18, the difference of pressure coefficient under different temperature is most obvious. The higher the temperature, the more serious the cavitation. In cavitation inception stage (where chord length ratio = 0~0.07), the pressure coefficient at three temperatures is basically the same, which indicates that cavitation has occurred at all temperatures in this region. In the cavitation development stage (where chord length ratio = 0.07~0.1), the effect of cavitation on the pressure was significantly different at each temperature and the higher temperature means larger cavity length. But in the cavitation disappearance stage (where chord length ratio = 0.1~0.18), the pressure is restored much faster; when the cavitation effect disappears, the pressure in the airfoil surface rises slowly at all temperatures, and the difference between them is narrowed gradually.
6. Conclusions
To select the suitable turbulent model for cavitation prediction, four common turbulent models and two filter turbulent models are solved in combination with Rayleigh-Plesset cavitation model, and their accuracies are compared. Based on the Rayleigh-Plesset equation, the source considering thermodynamic effects is added to the Zwart cavitation model using the Taylor formula and the Clapeyron formula. Two cavitation models considering thermodynamics effects, TCM1 and TCM2, are proposed and tested by comparing the result with the experimental data of NACA0015. To improve the accuracy of the currently proposed cavitation model, the evaporation and condensation coefficient are further corrected. The following are found from the current theoretical and numerical research:(i)The accuracy of a simulation result for cavitation flows is affected by the turbulent model in large degree. The cavity length simulated using the filtered turbulence model is shorter than that of the typical turbulence model, and the result of RNG k-ε turbulence model is most close to the experimental value.(ii)According to the simulation results using TCM1 and TCM2 models, the fluid temperature has an effect on the cavitation degree around the hydrofoil, and the cavitation becomes more intensive with the increase of the temperature. The results calculated using these two cavitation models are similar, the difference between them shows a gradually increasing trend with the increase of temperature, but the overall gap is small.(iii)The evaporation and condensation coefficient has a great effect on the simulation result of the cavitation flow under different temperature, and the value of 10 and 0.002 for evaporation and condensation coefficients are most suitable for the currently proposed cavitation model. The simulation result using the modified cavitation model agrees well with the experimental result, which provides a new method to simulate the cavitation flow considering the thermodynamic effects.
Nomenclature
| Reference pressure [Pa = N/m2] | |
| : | Vaporization pressure [Pa = N/m2] |
| Density [kg/m3] | |
| Cavitation number = | |
| Pressure coefficient = | |
| : | Reference velocity [m/s] |
| Time [s] | |
| Position [m] | |
| Velocity [m/s] | |
| Temperature [K] | |
| Mass fraction of water vapor | |
| Latent heat of vaporization | |
| Volume fraction | |
| Enthalpy [kJ/kg] | |
| Laminar viscosity | |
| Turbulent viscosity coefficient | |
| Unit tensor | |
| Pr: | Prandtl number |
| Bubble radius | |
| Surface tension coefficient | |
| Liquid kinematic viscosity. |
| Mixture | |
| Dimension of tensors | |
| Liquid | |
| Vapor. |
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was sponsored by the National Natural Science Foundation of China (no. 51679111, no. 51579118, and no. 51409127), PAPD, Six Talents Peak Project of Jiangsu Province (no. HYZB-002), Key R&D Projects in Jiangsu Province (BE2015119, BE2015001-4, BE2016319, and BE2017126), the Natural Science Foundation of Jiangsu Province (no. BK20161472), Science and Technology Support Program of Changzhou (no. CE20162004), and Scientific Research Start Foundation Project of Jiangsu University (no. 13JDG105).