Numerical implementation and verification of the Zwart-Gerber-Belamri cavitation model with the consideration of thermodynamic effect

In the world of pumps, turbopumps, and industrial applications, the absorption and release of latent heat in high-temperature zones play a pivotal role in the cavitation process. This phenomenon exerts a profound influence on system performance, especially when handling hot and cryogenic fluids, which are particularly susceptible to thermodynamic influences. The pressing need arises for a robust computational approach to address the intricacies of thermodynamic effects on cryogenic cavitation. The study at hand introduces a modification to the Zwart-Gerber-Belamri cavitation model, incorporating thermodynamic considerations to simulate quasi-steady cavitating flow around a NACA0015 hydrofoil. The SST k-ε turbulence model and homogeneous mass transfer cavitation model are employed to account for thermal effects, while the Clausius-Clapeyron equation is utilized to adjust the saturated vapor pressure within the cavitation model. Comparing the results with experimental data from Cervone et al. [1], especially in the thermal domain, reveals congruence in the estimated pressure and temperature drop (ΔT) within the cavity under varying free stream temperature conditions. Notably, thermodynamic effects exert a significant influence on cavitation dynamics during the phase-change process, potentially hindering or delaying cavitation in hot and cryogenic fluids. This enhanced cavitation model offers a marginally improved prediction of cavitation in water.


Introduction
When the pressure inside a liquid falls below its vapor pressure, a complicated phenomenon known as cavitation takes place.It can cause problems in fluid machinery, leading to noise, vibration, erosion, and reduced performance [2].While cavitation in water at room temperature is commonly considered an isothermal process, this is not true for thermosensitive fluids like hot water.The study of cavitation in thermosensitive fluids involves heat transfer, thermodynamics, fluid mechanics, and combustion, with a focus on the phase-change process of fluids.The thermodynamic effect of cavitation in thermosensitive fluids is influenced by the low liquid-to-vapor density ratio.When cavitation occurs, the fluid undergoes a phase transition from liquid to vapor, absorbing heat from the surrounding fluid.This causes a decrease in temperature within the cavitation zone [3].The temperature also affects properties such as vapor pressure, resulting in significant changes to cavitation dynamics.
Thermosensitive fluids include a variety of substances like cryogenic fluids, refrigerants, fluoroketones, and high-temperature water.Cavitation in thermosensitive fluids presents both challenges and opportunities for studying the complexities of thermo-fluids.It has practical importance

Governing equations and numerical methods
In Cartesian coordinates, under homogeneous fluids modelling, the governing equations for cavitating flows consist of the conservative version of the continuity, momentum, energy, and cavitation model equations.These equations are provided below: ) 3 (1 ) (1 ) The The variables in this equation are denoted as ρm, u, and p, which, respectively, stand for the density, velocity, and pressure of the mixed species.The coordinate axes are denoted by the subscripts i, j, and k along with the letter x.The terms μm and μt stand for the mixture's dynamic and turbulent viscosities, respectively, whereas λm and λt stand for the thermal conductivity and turbulent thermal conductivity.The following symbols stand for temperature, latent heat, specific heat, volume fraction, transfer source term, condensation source term, evaporation source term, and Kronecker delta, in that order: T, L, Cpm, α, m  , m +  , m −  and δij.

Extensional Zwart-Gerber-Belamri (EZGB) model
The Rayleigh-Plesset (R-P) equation is utilized by the Zwart-Gerber-Belamri cavitation model [18] to describe the evaporation and condensation phases.Moreover, the Zwart cavitation model develops a volume fraction transport equation for the vapor phase, which can be represented as: The surface tension and second derivative of time with respect to bubble radius effects are ignored in the following estimation, which leads to the R-P equation.In the equation, pv(T∞) stands for the saturation vapor pressure at a local thermodynamic state, S stands for surface tension, p indicates the local fluid pressure and RB stands for the spherical bubble radius.
The term m  for the transfer source can be written as follows [22]: ( ) Temperature changes have a dramatic impact on material properties and phase transition in thermosensitive fluids, hence the thermal effect must be considered [2].
In this study, the local saturation vapor pressure is determined through the coupling of pressure and temperature.Experimental findings have also highlighted the significant impact of turbulence on cavitating flow.The simulated saturation vapor pressure is typically lower than the actual saturation vapor pressure and is often adjusted using the Clausius-Clapeyron equation, followed by correction through turbulent kinetic energy.Hence, the saturated vapor pressure is expressed as follows [17]: The R-P equation incorporates thermal effects by introducing a term that considers vapor pressures associated with the bubble and liquid temperatures and is rewritten as follows; There are three elements on the right side of the equation: the incondensable gas term, the pressure term (representing mechanical forces), and the thermal effect term.Based on references Brennen [23] and Franc et al. [24], the equation incorporates Pv(Tb) to represent the vapor pressure at the bubble temperature and Pv(T∞) to represent the vapor pressure at the temperature of the liquid away from the bubble.
In the case of isothermal conditions (Tb = T∞), these vapor pressures are equal.However, when is assumed not to be isothermal, the equation accounts for the temperature difference between the upstream liquid and the bubble due to the phase transition.Consequently, the vapor pressure at the bubble temperature is lower than the vapor pressure at the upstream liquid temperature outside the cavitation zone.For cases where there are slight temperature differences, such as in simulated scenarios, the Taylor expansion method can be utilized to calculate the temperature difference accurately.As a result, the right-hand side of the R-P equation exhibits the following structure.
The temperature difference (Tb -T∞) results from a balance between the latent heat involved in condensation or evaporation and the heat exchanged between the liquid and the bubble.( ) In this equation, ρv represents the vapor density, L denotes the liquid latent heat of evaporation, and hb signifies the bubble convective heat transfer coefficient.By substituting Equation ( 15) into Equation ( 13) and applying the Taylor expansion, the R-P equation is transformed into the following form; The relationship between the temperature difference ΔT between the bubbles and the surrounding liquid and the heat flow Q per unit area A, transported through convection or phase transition, is characterized by the heat transfer coefficient, abbreviated as hb.
Selecting an appropriate value of hb is crucial for accurately predicting simulation outcomes [25].When hb is set to a high value, it can cause an overestimation of the cavity volume in the prediction.Conversely, when hb is too low, it can lead to a significant underestimation of the predicted cavity volume.The equation provided by Oresta et al. [26], which is shown as Eq. ( 18) below, is used to locally determine the convective heat transfer coefficient, hb, for the purposes of the present study.
( ) The results from Ref. [27] are fitted to determine the parameter "n," which in this study is considered to be 2.65.Reb can be calculated using the formula 2RB|v -u|/λl, where "v" and "u" stand for the bubble and liquid velocities, respectively.In addition, the terms cpl, µl and δl stand for the specific heat, dynamic viscosity, and thermal conductivity of the liquid, respectively.The formula for calculating the thermal diffusivity of a liquid, λl is λl = δl/(ρcp), where cp is the specific heat at constant pressure and ρ is the density of the liquid.The bubble radius is indicated by RB.Ja is referred to as the sensible heat to latent heat ratio, and Nub,0 is a low bubble Péclet number's Nusselt number.While "Peb" stands for the bubble Péclet number, "Pec" stands for the crossover Péclet number.The volume change of bubbles brought on by evaporation or condensation is constrained by a large latent heat value, which is indicated by a small value of the Jakob number (Ja).
Two terms are used to describe thermal effects: the term for the latent heat energy source mL  and the term for convective heat transfer, Chb.Convective and condensation heat transfer at the liquid-vapor interface are closely correlated, necessitating this adjustment.The pressure-based mass transfer word is specifically subtracted from the convective heat transfer term, and the latent heat energy source term mimics the mass transfer rate but excludes the thermal term.The thermal term is given as follows; ( ) With the inclusion of thermal effects in cryogenic cavitating flow, the EZGB model can now be described as follows: ( ) ( ) where the saturation temperature, Tsat(p), is a function of the local pressure.
The empirical coefficients Fvap and Fcond play a crucial role in accurately simulating thermal cavitating flows.Sun et al. [28] re-evaluated these coefficients for cryogenic cavitating flows, setting Fvap to 5.0 and Fcond to 0.01.This calibration led to more precise predictions of temperature and pressure.Therefore, in all subsequent simulations, we have set the values of Fvap and Fcond to 5.0 and 0.01, respectively.Additionally, a nucleation site volume fraction of αnuc = 5×10 -4 has been assigned, and the bubble radius RB has been fixed at 1×10 -6 m.

Turbulence model
The standard k-ε model and its improved RNG k-ε version are highly accurate in simulating fully developed turbulent flow on far wall surfaces.However, their precision is limited when it comes to simulating flow separation and large backpressure gradients.To address this limitation, Menter [29] introduced the two-equation SST k-ω model, a modification of the BSL k-ω model.Strong curvature flow and separated flow simulation accuracy have both been greatly improved by this new model.
The k-ω SST (Shear Stress Transport) turbulence model, which combines the benefits of the k-ε and k-ω models, is the one we used in the present study.In contrast to how it is applied, the k-ω model is implemented close to the wall [30].The typical issue with k-ω models, which is their susceptibility to inlet free-stream turbulence parameters, is resolved by the k-ε SST model by switching to k-ω behavior in the free-stream.Since the model successfully predicts boundary layer detachment characteristics, it is the model of choice for real-world engineering flow computations.The equation for the model is expressed as follows: For the turbulent kinetic energy equation, For the turbulence frequency equation, For the turbulent viscosity equation, where F1 and F2 are blending functions; S represents the shear tensor; Pk stands for the turbulent kinetic energy term generated by viscous force, and Pkb represents the term for turbulent kinetic energy created by buoyancy.The constant coefficients are: a1 = 0.31, σk = 2, β′ = 0.09, σω = 2, α = 5/9, β = 0.075, and σω2 = 1/0.856.

Numerical setup and description
An established cavitation experiment performed by Cervone et al. [1] is used to validate the numerical results in this paper.Experimental research was conducted in the Cavitating Pump Rotor Dynamic Test Facility (CPRTF) lab at various temperatures of water surrounding the NACA0015 hydrofoil and the performance of cavitation was examined in relation to the thermodynamic effect.

NACA0015 hydrofoil geometry model
The computational domain for the NACA0015 hydrofoil experiment has a hydrofoil chord length of 115 mm.The angle of attack is 5 o , and the experimental section features pressure measurement holes at the inlet, exit, along the suction surface, and along the pressure surface.The 2D model of the hydrofoil is shown in Figure 2, positioned in a 0.12m high channel.The time step for accurate simulation of transient cavitating flow is 7.19×10 -5 s, based on recommendations by Coutier-Delgosha et al. [31].A mesh independence study reveals that a medium resolution mesh with around 151,107 nodes provides sufficient accuracy (Table 1), and a slightly increased time step of 1.0×10 -4 s does not significantly impact the results.

Non-cavitation calculation
For this study, we also aim to understand how the hydrofoil performs in non-cavitating normal operating conditions.By studying the non-cavitating behavior of the hydrofoil, we gained insights into the lift and drag forces, flow patterns, and overall efficiency of the hydrofoil.Secondly, the transient calculation process is first started from steady non-cavitating flow field to avoid any vapor fraction at the initial time step.We also validated the numerical approach by examining the steady non-cavitating flow around the hydrofoil.Figure1(a) shows the pressure coefficient distribution around the hydrofoil without cavitation.The predicted pressure distribution along the suction side of the hydrofoil matches the experimental data remarkably well, with no significant discrepancies observed between the predicted and measured results.Figure 3. Mesh generation for hydrofoil.

Numerical methodology
User-defined functions (UDF) were incorporated into the simulations using the ANSYS FLUENT CFD program.The continuity, momentum, and energy equations were solved using the pressure-based solver.
For convective fluxes, a second-order upwind differencing scheme was employed, while a central differencing scheme was applied for diffusive fluxes.The calculation of pressure at cell faces utilized the body-force-weighted scheme [32].
The coupling between velocity and pressure was achieved using the SIMPLEC algorithm at double precision [33,34].Second-order backward Euler algorithms were used for time integration while second-order upwind algorithms were used for spatial derivatives [35].After the convergence of the initial steady non-cavitating flow, unsteady cavitating flow simulations were initialized.Properties that are dependent on temperature exhibit different variations as the temperature changes.Figure 4 displays the selected properties as a function of temperature, with Cl and ρv exhibiting a more rapid change with temperature than Cv and L. Furthermore, as temperature increases, the slope of the saturation curve for ρv becomes steeper, indicating that the properties undergo varying changes with changes in temperature.

Results and discussion
The study utilizes the extensional Zwart-Gerber-Belamri (EZGB) model to investigate the thermo-fluid cavitation behaviors in water.The flow conditions and results for cavitating flows at different free stream temperatures are presented in Table 2. Significant variations in the results are observed.The modified numerical model successfully captures the primary characteristics of the pressure coefficient profile, as depicted in Figure 5, using the thermal energy solution.The local saturated vapor pressure falls with decreasing temperature, resulting in smaller temperature reductions via evaporative cooling in the EZGB model, which results in smaller cavity sizes.Evaporation absorbs latent heat, causing the temperature to drop slightly below the reference fluid temperature.This leads to a negative coefficient of pressure (Cp) in the cavity region, as shown in Figure 5. Temperature variations in thermo-fluid cavitating flow can significantly affect material properties, resulting in notable changes in pressure fields and cavity structures.It can be observed from Figure 5(c) that there is a significant difference between the three plot compared, at higher temperatures, the EZGB model clearly predicts the experimental points better than the EZGB model in Xu et al. [21] and the ZGB model especially at the cavitation closure region.The introduction of the thermodynamic parameter Σ [36], which is defined by ( ) allows for a precise discussion of the thermodynamic effects.The parameter Σ, which is called the Thermodynamic parameter and is measured in units of m/s 3/2 , is a key factor in determining the dynamic response of an individual spherical bubble.It is computed using the Rayleigh-Plesset equation, together with the energy balance and heat diffusion equations [36].The parameter is employed to forecast the variation in temperature as the diameter of the bubble changes.
In above equations ρl [kg/m³] is the liquid density, ρv [kg/m³] is vapor density at mainstream pressure and temperature.K is the thermal conductivity and L [J/kg] is latent heat.Cpl [J/(kg K)] and αl [m²/s] are isobaric specific heat and thermal diffusivity of water (functions in an unsteady heat transfer process).C [m] is chord length of hydrofoil.The cavity size is affected by the density ratio between liquid and vapor phases.Lower density ratios require increased vaporization to maintain a constant cavity size.A greater temperature drop during phase change (ΔT) leads to an increase in latent heat (L) and a decrease in specific heat (Cv), resulting in an elevated thermodynamic parameter Σ. Figure 6 illustrates the variation in Σ for water concerning normalized temperature (T/Tc, where Tc is the critical temperature).As temperature rises, Σ experiences an increment due to temperature-dependent material properties, indicating more pronounced thermodynamic effects at higher temperatures under equivalent free-stream conditions.The thermodynamic parameter Σ indicates the strength of thermodynamic effects in different fluids.Liquid hydrogen exhibits higher Σ values than nitrogen, suggesting more pronounced thermodynamic effects.At low temperatures, water has negligible thermodynamic effects, but at around 350K, its effects become similar to nitrogen, indicating that hot water cavitating flow can experience significant thermodynamic effects.The intensity of these effects is nearly the same for water, nitrogen, and hydrogen at specific temperatures.Raising the water temperature results in elevated Σ values, emphasizing the significance of thermodynamic effects in both hot water cavitating flow and cryogenic fluids.These effects are reliant on thermophysical properties rather than the type of fluid, given equivalent freestream conditions.Hence, accurately predicting the intensity of thermodynamic effects in thermo-fluid cavitating flows necessitates careful consideration of the critical thermodynamic parameter Σ and a comprehensive evaluation of physical properties.
The EZGB cavitation model enhances the accuracy of numerical pressure coefficient predictions compared to experimental data at different free stream temperatures.Additionally, compared to the ZGB model at the same temperature, it displays a larger region of low pressure on the hydrofoil suction surface.As water gets closer to its critical temperature, it becomes more sensitive to changes in temperature, especially if there is a pressure gradient inside the cavitation closure region.The revised thermal term of the EZGB model primarily influences the cavitation rate around the closure region, leading to improved pressure predictions.This modification enables a greater release of the heat of vaporization, thereby restraining cavitation intensity more effectively at higher temperatures.
Figure 7 illustrates how the thermodynamic processes and the thermosensitivity of water at high temperatures play a role in cavitation phenomena and how the EZGB cavitation model agrees with experimental findings.This study emphasizes the influence of thermodynamic influences on cavitation behaviour by modelling water cavitation at varied free stream temperatures.

Conclusion
The present study investigates the influence of thermodynamic effects on cavitation and the thermosensitive behavior of water at high temperatures.Specifically, the study extends the Zwart cavitation model to incorporate cavitating flows around a NACA0015 hydrofoil in water at different temperatures.The results demonstrate that the extended Zwart-Gerber-Belamri (ZGB) model provides a modest enhancement in predicting cavitation characteristics in hot water cavitating flow by considering thermodynamic effects.
Furthermore, the study highlights that thermodynamic effects are predominantly governed by thermophysical properties rather than the type of fluid under equivalent freestream conditions.In this context, the EZGB model demonstrates improved accuracy in predicting the numerical pressure coefficient compared to the experimental pressure coefficient across various freestream temperatures.Notably, the model reveals a wider region of low pressure on the suction surface of the hydrofoil in comparison to the ZGB cavitation model at the same temperature.
The limitations of this study might include the over simplified assumptions of the EZGB model.For instance, the model's assumptions about uniform properties and idealized bubble behavior might not fully capture the complex dynamics of cavitating flows.Also, the study primarily focuses on a limited range of temperatures in water.Extending the analysis to a wider temperature range and exploring different fluids could provide a more comprehensive understanding of the impact of thermodynamic effects on cavitation.Thirdly, the model's use of the Rayleigh-Plesset equation and the assumption of spherical bubbles might not accurately represent the actual shapes and behavior of cavities in various flow conditions.In the future performing detailed experiments under a wider range of conditions and comparing them rigorously with model predictions would help validate and refine the EZGB model.Also Applying the model to practical engineering systems, such as marine propulsion or hydraulic machinery, and comparing results with real-world performance data can validate the applicability of the model and provide insights for design and optimization.Thirdly, conducting sensitivity analyses to identify the most influential parameters in the model can guide efforts to improve its accuracy and reliability.

Figure 1 .Figure 3
Figure 1.(a) The pressure coefficient on the suction side of the NACA0015 hydrofoil under non-cavitating conditions for water flow.(b) The distribution plot of y + .

Figure 4 .
Figure 4. Variation of physical properties along saturation curve.

Figure 6 .
Figure 6.The parameter Σ for water, liquid nitrogen, and liquid hydrogen expressed as a function of normalized temperature T/Tc.

Figure 7 .
Figure 7. Predicted results of (a) temperature, and (b) mass transfer rate in water at 298K, 323K, and 343K.

Table 2 .
Setting of parameters at different temperatures.