Study on the Influence of Thermodynamic Effects on the Characteristics of Liquid Nitrogen Cavitating Flow around Hydrofoils

1. Introduction

Liquid cryogenic media exhibit unique characteristics, such as ultra-low temperatures, low thermal conductivity, extremely low liquid-to-vapor density ratio, and highly unstable critical states. Fluid machinery employed in the transportation of cryogenic media must withstand severe environmental conditions while maintaining excellent cavitation performance and operational stability. For instance, in order to achieve higher thrust and specific impulses, cryogenic liquid fuel is commonly utilized as propellant in liquid-fueled rocket engines. The cavitation performance of turbine pumps employed for cryogenic fuel transportation directly impacts the overall propulsion performance and manufacturing costs of the engine [1,2]. In the offshore liquid cargo ship transportation system, cavitation occurrence in liquid cargo pumps can induce structural vibrations, noise, performance degradation, increased transportation costs, and compromised transportation safety. Hence, understanding the unique properties of cavitation in cryogenic fluids and developing accurate prediction methods for such phenomena pose urgent challenges that need to be addressed.

Due to the unique properties of cryogenic media in fluid machinery and cryogenic engineering, the cavitation process induces a significant temperature decrease in the cavitation region. This temperature decrease causes a substantial reduction in the local saturated vapor pressure of the liquid, thereby impeding further cavitation development. This phenomenon, resulting from temperature variations, is commonly known as the thermodynamic effect [3]. Initially proposed to quantitatively characterize the magnitude of cavitation heat effects and predict cavitation performance in pumps, the B factor theory was introduced. The B factor is defined by Stahl et al. [4], represents the ratio of vapor volume to liquid volume during the evaporation process. Subsequently, various methods have been proposed by researchers to estimate the B factor. However, cryogenic cavitation flow involves complex coupling mechanisms among vortices and bubbles of different scales and frequencies, as well as the influence of temperature. As a result, it becomes an extremely challenging flow phenomenon, and the results obtained from theoretical analysis often deviate significantly from actual situations.

Numerical simulation technology is extensively employed in various fields, including aerospace, aviation, shipping, and energy, significantly reducing the design cycle for relevant products. Numerical simulation has also been applied to the study of cryogenic cavitation. However, due to the highly unstable nature of cryogenic cavitation, researchers often modify existing turbulence and cavitation models to accurately capture unstable flow characteristics. For instance, Rodio et al. [5] modified the heat transfer coefficient in the two-fluid model for cryogenic cavitation. Ahuja et al. [6] adjusted the empirical coefficient in the Merkle cavitation model to simulate cavitation flow of liquid hydrogen and liquid nitrogen. The full cavitation model was adapted by Zhang X.B. et al. [7,8] to accommodate numerical calculations of cryogenic cavitation flow. Sun T.Z. et al. [9] corrected the evaporation and condensation coefficients in the Zwart model, resulting in numerical results that were more consistent with Hord experiment results. The classical model was modified by Zhang et al. [10,11], who adjusted the empirical coefficients of the Zwart model and incorporated the heat source term. Li W.G. et al. [12] proposed a new cavitation model that considers thermodynamic effects and is applicable for numerical calculations of various cavitation mechanisms. Moreover, due to the complexity of cavitation mechanisms, the evolution of bubbles is closely related to multiscale vortex motion. Therefore, the development of more accurate turbulence models for cavitation numerical calculations is crucial. The Partially-Averaged Navier-Stokes (PANS) model, with the ability to transition from Reynolds-averaged turbulence models to direct numerical simulation, holds certain advantages [13,14]. Despite the advantages of numerical simulation, the availability of experimental data for verifying the accuracy of cryogenic cavitation simulations is limited. This is mainly due to the significant influence of equipment and technology on cryogenic cavitation experiments. Notably, Hord conducted a series of cryogenic cavitation experiments from 1972 to 1974 with the support of NASA, and the results have been widely used to validate numerical calculations [15,16,17,18]. Kazuki et al. [19] performed cryogenic cavitation experiments in a water tunnel and verified the inhibitory effect of thermal effects on cavitation development. It was confirmed by Cervone et al. [20,21] that the size of cavitation bubbles can be influenced by water temperature through induced wheel and hydrofoil experiments. Chen et al. [22,23] conducted experimental research on the characteristics of cryogenic nitrogen cavitation flow in contraction-expansion tubes.

This paper focuses on investigating the cavitation flow around a hydrofoil at low temperatures. The DDES turbulence model and Sauer-Schnerr cavitation model are employed to conduct numerical simulations on the Hord experiment hydrofoil, allowing for the verification of numerical calculation accuracy. On this basis, numerical simulations are carried out, and the thermodynamic effect, coupled motion characteristics of cavitation and vortices in the unsteady cavity flow of liquid nitrogen around the NACA 0015 hydrofoil are analyzed.

2. Numerical calculation methods and verification

2.1. Bubble number density correction calculation model

When conducting numerical simulations of cavitation flow using ambient temperature water as the medium, the typical setting for the bubble number density in the Sauer-Schnerr cavitation model is 10¹³. However, liquid nitrogen, being a cryogenic fluid, exhibits significant differences in physical properties compared to ambient temperature water. Therefore, it is necessary to modify the bubble number density in the cavitation model for liquid nitrogen.

To address the modeled stress loss issue in Detached Eddy Simulation (DES), the Delayed Detached Eddy Simulation (DDES) method is employed. DDES introduces a delay function to reconstruct the DES length scale, considering the grid scale and vortex viscosity field [24,25]. This approach enables more accurate capture of flow structures in the primary cavitation region, making DDES the preferred turbulence model for this study.

The numerical simulations are validated against the Hord hydrofoil liquid nitrogen flow experiments. In order to save computational resources, a two-dimensional hydrofoil grid is used for the calculations. The inlet width of the hydrofoil is 25.4 mm, and the total length of the flow passage is 177.2 mm. The computational domain is shown in Figure 1. In order to ensure y⁺≤10, the computational domain is divided into quadrilateral grids, and the grid refinement technique is applied near the wall surfaces. The inlet is set as a velocity inlet, the outlet as a pressure outlet, and the walls are assumed to be adiabatic and have no slip boundary conditions. During the calculations, the outlet pressure is adjusted continuously to match the measured inlet pressure values from the Hord experiments. The averaged results from different time steps of the numerical calculations are compared with the experimental values.

Figure 1. Schematic diagram of computational domain.

Figure 1. Schematic diagram of computational domain.

At the operating condition of 290C, a grid independence analysis is being conducted on the computational domain to ascertain whether the pressure difference at the inlet and outlet stabilizes as the number of grids increases. For this purpose, five grid partitioning schemes have been designed, and the corresponding number of grids for each scheme is presented in Table 1.

Table 1. Different mesh numbers division schemes.

Grid division schemes	grid number	pressure difference ΔP(Pa)
1	95,844	89258
2	176,548	115879
3	295,858	148658
4	534,048	167155
5	715,489	169518

Table 1. Different mesh numbers division schemes.

Grid division schemes	grid number	pressure difference ΔP(Pa)
1	95,844	89258
2	176,548	115879
3	295,858	148658
4	534,048	167155
5	715,489	169518

Figure 3 illustrates the results of the grid independence verification. It can be observed that the inlet-outlet pressure difference remains relatively stable when the number of grids exceeds 534,048. Considering the trade-off between computational efficiency and accuracy, a total of 534,048 grids were selected for the calculations. The grid distribution across the computational domain is depicted in Figure 2.

Figure 2. (a) Schematic diagram of mesh; (b) Wall Yplus for Hord hydrofoil.

Figure 2. (a) Schematic diagram of mesh; (b) Wall Yplus for Hord hydrofoil.

Figure 3. Grid independence verification.

Figure 3. Grid independence verification.

The numerical calculations for this study involve solving the control equations, which include the continuity equation, momentum equation, and energy equation, as described in the following text:

$$\frac{\partial {\rho }_m}{\partial t}+\frac{\partial }{\partial x_j}\left({\rho }_m{u}_j\right)=0$$

(1)

$$\frac{\partial }{\partial t}\left({\rho }_m{u}_i\right)+\frac{\partial }{\partial x_j}\left({\rho }_m{u}_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\mu +{\mu }_t\right)\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right]$$

(2)

$$\frac{\partial \left({\rho }_m{c}_{p}T\right)}{\partial t}+\frac{\partial }{\partial x_j}\left({\rho }_m{u}_j{c}_{p}T\right)=\nabla \cdot \left(k_{eff}\nabla T\right)-\frac{\partial \left({\rho }_m{f}_{v}L\right)}{\partial t}-\frac{\partial \left({\rho }_m{u}_j{f}_{v}L\right)}{\partial x_j}$$

(3)

where, ρ_m=α_lρ_l+α_vρ_v represents the mixture density of the two phases; u is velocity; μ is dynamic viscosity of the fluid; μ_t is turbulent viscosity; k_eff is thermal conductivity coefficient; L is latent heat of vaporization; C_p is specific heat capacity at constant pressure; i and j represent the coordinate directions, while v and l represent the vapor and liquid phases, respectively.

Sauer-Schnerr cavitation model [26]:

$${\dot{S}}_e=\frac{3{\alpha }_v\left(1-{\alpha }_v\right)}{R_B}\frac{{\rho }_v{\rho }_l}{{\rho }_m}{\left[\frac{2}{3}\frac{P_v\left(T\right)-P}{{\rho }_l}\right]}^{1/2},P<P_v\left(T\right)$$

(4)

$${\dot{S}}_c=-\frac{3{\alpha }_v\left(1-{\alpha }_v\right)}{R_B}\frac{{\rho }_v{\rho }_l}{{\rho }_m}{\left[\frac{2}{3}\frac{P-P_v\left(T\right)}{{\rho }_l}\right]}^{1/2},P\ge P_v\left(T\right)$$

(5)

R_B (Bubble radius) can be expressed by the following formula：

$$R_B={\left(\frac{{\alpha }_v}{1-{\alpha }_v}\frac{3}{4\pi }\frac{1}{n}\right)}^{1/3}$$

(6)

where, α_v is vapor volume fraction; ρ_l, ρ_v, ρ_m is liquid density, vapor density and mixture density; R_B is bubble radius; n is bubble number density in the liquid.

The DDES turbulent transport equation based on the SST k-ω model is:

$$\frac{\partial (\rho k)}{\partial t}+\frac{\partial (\rho k{{\displaystyle u}}_i)}{\partial {{\displaystyle x}}_i}=\frac{\partial }{\partial x}\left[\left(\mu +\frac{{{\displaystyle \mu }}_t}{{{\displaystyle \sigma }}_{k3}}\right)\frac{\partial k}{\partial {{\displaystyle x}}_j}\right]+{{\displaystyle P}}_k-\rho {{\displaystyle k}}^{3/2}/{{\displaystyle l}}_{DDES}$$

(7)

$$\begin{array}{l}\frac{\partial (\rho \omega )}{\partial t}+\frac{\partial (\rho \omega {{\displaystyle u}}_i)}{\partial {{\displaystyle x}}_i}=\frac{\partial }{\partial {{\displaystyle x}}_j}\left[\left(\mu +\frac{{{\displaystyle \mu }}_t}{{{\displaystyle \sigma }}_{\omega 3}}\right)\frac{\partial \omega }{\partial {{\displaystyle x}}_j}\right]+{{\displaystyle \alpha }}_3\frac{\omega }{k}{{\displaystyle P}}_k-{{\displaystyle \beta }}_3\rho {{\displaystyle \omega }}^2\\ +2(1-{{\displaystyle F}}_1)\rho \frac{1}{\omega {{\displaystyle \sigma }}_{\omega 2}}\frac{\partial k}{\partial {{\displaystyle x}}_j}\frac{\partial \omega }{\partial {{\displaystyle x}}_j}\end{array}$$

(8)

$${{\displaystyle \mu }}_t=\rho \frac{{{\displaystyle a}}_{1}k}{\max ({{\displaystyle a}}_1\omega ,S{F}_2)}$$

(9)

where the mixed functions F₁ and F₂ of the SST k-ω model can be represented as:

$${{\displaystyle F}}_1=\tanh ({{\displaystyle \xi }}^4)$$

(10)

$$\xi =\min \left[\max \left(\frac{\sqrt{k}}{{{\displaystyle \beta }}^{*}\omega {{\displaystyle d}}_{\omega }},\frac{500\mu }{\rho {{{\displaystyle d}}_{\omega }}^2\omega }\right),\frac{4\rho k}{{{{\displaystyle D}}_{\omega }}^{+}{{\displaystyle \sigma }}_{\omega 2}{{{\displaystyle d}}_{\omega }}^2}\right]$$

(11)

$${{{\displaystyle D}}_{\omega }}^{+}=\max \left[2\rho \frac{1}{{{\displaystyle \sigma }}_{\omega 2}}\frac{1}{\omega }\frac{\partial k}{\partial {{\displaystyle x}}_j}\frac{\partial \omega }{\partial {{\displaystyle x}}_j},{{\displaystyle 10}}^{-10}\right]$$

(12)

$${{\displaystyle F}}_2=\tanh \left({{\displaystyle \eta }}^2\right)$$

(13)

$$\eta =\max \left\{\frac{2{{\displaystyle k}}^{1/2}}{{{\displaystyle \beta }}^{*}\omega {{\displaystyle d}}_{\omega }}\right.,\left.\frac{500\mu }{\rho {{{\displaystyle d}}_{\omega }}^2\omega }\right\}$$

(14)

where, d_ω is the distance from the calculation point to the wall surface; P_k is the turbulence generation term caused by viscous forces, defined the same as the DES model; α₁ = 5/9, β₁ = 0.075, k₁ = 1.176, σ_ω1 = 2, α₂ = 0.44, β₂ = 0.0828, σ_k2 = 1, σ_ω2= 1/0.856, a₁ = 0.31, β ^* = 0.09.

$${{\displaystyle l}}_{DDES}={{\displaystyle l}}_{RANS}-{{\displaystyle f}}_d\max \left(0,{{\displaystyle l}}_{RANS}-{{\displaystyle l}}_{LES}\right)$$

(15)

$${{\displaystyle l}}_{RANS}=\frac{{{\displaystyle k}}^{1/2}}{{{\displaystyle \beta }}^{*}\omega }$$

(16)

$${{\displaystyle l}}_{LES}={{\displaystyle C}}_{DES}\Delta$$

(17)

$$C_{DES}=F_1{C}_{DES1}+\left(1-F_1\right)C_{DES2}$$

(18)

$$\Delta =\max \left\{\Delta x,\Delta y,\left.\Delta z\right\}\right.$$

(19)

$$f_d=1-\tanh \left[{\left(C_{d1}{r}_d\right)}^{C_{d2}}\right]$$

(20)

$$r_d=\frac{v_t+v}{\sqrt{\frac{1}{2}\left(S^2+{\Omega }^2\right)}{k}^2{d}_{\omega }{}^2}$$

(21)

where, f_d is the delay function; r_d is the delay factor; S is the value of the strain rate tensor; Ω is the value of the curl tensor; Δ is the maximum side length of the unit; Constant k=0.41, C_DES1 =0.78, C_DES2=0.61, C_d1=8, C_d2=3.

2.2. Discussion on the results of bubble number density correction

Initially, numerical simulations were performed for the swirling flow experiment of liquid nitrogen at 290C, considering different bubble number densities. The results, depicted in Figure 4, reveal notable differences in pressure and temperature distributions compared to the experimental values when the bubble number densities were set at 10⁷ and 10⁹. However, with a bubble number density of 10⁸, the numerical calculations exhibit the closest agreement with the experimental results.

Figure 4. Comparing simulated and experimental pressure and temperature drops. (290C operating condition).

Figure 4. Comparing simulated and experimental pressure and temperature drops. (290C operating condition).

Figure 5 and Figure 6 display the results obtained by numerically simulating two operating conditions, 293A and 296B, using the modified bubble number density. By comparing the measured values from the experiments with the simulated values at the same measurement points in Table 2, Table 3, Table 4 and Table 5, it can be observed that for the 293A operating condition, the maximum relative error in pressure drop is 11.4%, and the maximum relative error in temperature drop is 0.43%. For the 296B operating condition, the maximum relative error in pressure drop is 12.4%, and the maximum relative error in temperature drop is 0.22%. It is feasible to adjust the bubble number density in the Sauer-Schnerr cavitation model to 10⁸ and use the modified model for numerical calculations of cryogenic cavitation flow.

Table 2. Comparing simulated and experimental pressure drops. (293A operating condition).

Distance x(mm)	Test pressure value Pexp(Pa)	Simulation pressure value Psim(Pa)	Relative error |Pexp-Psim|/Pexp
1.651	89400	97225.8	8.8%
5.08	90000	98865.9	9.9%
10.16	98700	105640.3	7.0%
19.05	240700	214538.8	10.9%
29.21	328900	291525.5	11.4%

Table 2. Comparing simulated and experimental pressure drops. (293A operating condition).

Distance x(mm)	Test pressure value Pexp(Pa)	Simulation pressure value Psim(Pa)	Relative error |Pexp-Psim|/Pexp
1.651	89400	97225.8	8.8%
5.08	90000	98865.9	9.9%
10.16	98700	105640.3	7.0%
19.05	240700	214538.8	10.9%
29.21	328900	291525.5	11.4%

Table 3. Comparing simulated and experimental temperature drops. (293A operating condition).

Distance x(mm)	Test temperature value Texp(K)	Simulation temperature value Tsim(K)	Relative error |Texp-Tsim|/Texp
2.032	76.35	76.53	0.23%
6.985	76.62	76.42	0.25%
13.335	77.05	76.72	0.43%
22.86	77.49	77.50	0.01%
32.512	77.48	77.57	0.12%

Table 3. Comparing simulated and experimental temperature drops. (293A operating condition).

Distance x(mm)	Test temperature value Texp(K)	Simulation temperature value Tsim(K)	Relative error |Texp-Tsim|/Texp
2.032	76.35	76.53	0.23%
6.985	76.62	76.42	0.25%
13.335	77.05	76.72	0.43%
22.86	77.49	77.50	0.01%
32.512	77.48	77.57	0.12%

Table 4. Comparing simulated and experimental pressure drops. (296B operating condition).

Distance x(mm)	Test pressure value Pexp(Pa)	Simulation pressure value Psim(Pa)	Relative error |Pexp-Psim|/Pexp
1.651	269800	282665	4.8%
5.08	275700	300560.2	9.0%
10.16	327900	368529.1	12.4%
19.05	472000	454998.1	3.6%
29.21	481300	470346.9	2.3%

Table 4. Comparing simulated and experimental pressure drops. (296B operating condition).

Distance x(mm)	Test pressure value Pexp(Pa)	Simulation pressure value Psim(Pa)	Relative error |Pexp-Psim|/Pexp
1.651	269800	282665	4.8%
5.08	275700	300560.2	9.0%
10.16	327900	368529.1	12.4%
19.05	472000	454998.1	3.6%
29.21	481300	470346.9	2.3%

Table 5. Comparing simulated and experimental temperature drops. (296B operating condition).

Distance x(mm)	Test temperature value Texp(K)	Simulation temperature value Tsim(K)	Relative error |Texp-Tsim|/Texp
2.032	86.57	86.73	0.18%
6.985	87.22	87.02	0.22%
13.335	88.08	88.02	0.07%
22.86	88.47	88.49	0.02%
32.512	88.43	88.49	0.07%

Table 5. Comparing simulated and experimental temperature drops. (296B operating condition).

Distance x(mm)	Test temperature value Texp(K)	Simulation temperature value Tsim(K)	Relative error |Texp-Tsim|/Texp
2.032	86.57	86.73	0.18%
6.985	87.22	87.02	0.22%
13.335	88.08	88.02	0.07%
22.86	88.47	88.49	0.02%
32.512	88.43	88.49	0.07%

Figure 5. Comparing simulated and experimental pressure and temperature drops. (293A operating condition).

Figure 5. Comparing simulated and experimental pressure and temperature drops. (293A operating condition).

Figure 6. Comparing simulated and experimental pressure and temperature drops. (296B operating condition).

Figure 6. Comparing simulated and experimental pressure and temperature drops. (296B operating condition).

3. Unsteady cavitation flow around a hydrofoil with liquid nitrogen

Considering the fundamental role of the NACA0015 hydrofoil as a design element in numerous turbomachinery blades, the same computational approach is employed to simulate the unsteady cavitation flow around the NACA0015 hydrofoil. This numerical simulation accounts for the variations in key parameters (density, specific heat capacity, dynamic viscosity, heat transfer coefficient) of the gas-liquid two-phase flow with temperature. The objective is to investigate the mechanisms underlying unsteady cryogenic cavitation flow in liquid nitrogen.

The computational model is centered around the NACA0015 hydrofoil, featuring a chord length (c) of 50mm and an angle of attack (α) of 7.5°, as depicted in Figure 7. The computational domain is discretized using quadrilateral grids. To accurately simulate the cavitation flow near the hydrofoil boundary, the wall region encompassing the entire hydrofoil is suitably refined to maintain that y⁺≤5. The key dimensionless parameters considered in the calculations encompass the incoming cavitation number (σ_∞), lift coefficient (C_l), and drag coefficient (C_d):

$${\sigma }_{\infty }=(P_{in}-P_v)/0.5{\rho }_l{u}_{in}^2$$

(22)

$$C_l=F_y/0.5{\rho }_l{u}_{in}^{2}c$$

(23)

$$C_d=F_x/0.5{\rho }_l{u}_{in}^{2}c$$

(24)

where, P_in is inlet pressure; P_v is saturation vapor pressure corresponding to inlet temperature; u_in is inlet velocity; ρ_l is liquid density corresponding to inlet temperature; F_x and F_y are drag and lift on the hydrofoil, respectively; c is chord length.

Figure 7. Computational domain for NACA0015 hydrofoil.

Figure 7. Computational domain for NACA0015 hydrofoil.

In order to assess the grid independence of the computational domain, the stability of the lift coefficient of the hydrofoil with respect to the increasing number of grids was examined. To achieve this, five different grid division schemes were devised, and the corresponding number of grids for each scheme is presented in Table 6.

Table 6. Different grid schemes for NACA0015 hydrofoil.

Grid division schemes	grid number	lift coefficient Cl
1	102,987	0.02298
2	215,875	0.02088
3	387,655	0.01880
4	477,732	0.01839
5	624,875	0.01828

Table 6. Different grid schemes for NACA0015 hydrofoil.

Grid division schemes	grid number	lift coefficient Cl
1	102,987	0.02298
2	215,875	0.02088
3	387,655	0.01880
4	477,732	0.01839
5	624,875	0.01828

Figure 9 presents the outcomes of the grid independence verification, revealing that the lift coefficient exhibits no significant variation when the number of grids reaches 477,732 or higher. To strike a balance between computational cost and accuracy, 477,732 grids were selected for the calculation. The grid distribution across the computational domain is illustrated in Figure 8.

Figure 8. (a)Computational domain mesh diagram of NACA0015 hydrofoil; (b)Wall surface Yplus diagram of NACA0015 hydrofoil.

Figure 8. (a)Computational domain mesh diagram of NACA0015 hydrofoil; (b)Wall surface Yplus diagram of NACA0015 hydrofoil.

Figure 9. Grid independence verification of NACA0015 hydrofoil.

Figure 9. Grid independence verification of NACA0015 hydrofoil.

During the calculation, the boundary conditions were specified as velocity inlet, pressure outlet, and all wall surfaces were modeled as adiabatic walls. To improve the convergence of the numerical calculation, a combination of the SST k-ω turbulence model and the pseudo-transient algorithm was initially employed for steady-state calculations. The results obtained from the steady-state calculation were subsequently utilized as the initial values for the numerical simulation of unsteady cavitating flow. The corresponding time step was set to 1×10^-5s.

Figure 10 displays the temporal evolution of the lift and drag coefficients of the hydrofoil for two scenarios: isothermal and thermodynamic effects. The calculation parameters encompass an inlet velocity (u_in) of 23m/s, inlet temperature (T_in) of 77K, fluid density (ρ_l) corresponding to the inlet temperature of 807.392 kg/m³, and the saturated vapor pressure of liquid nitrogen (P_v) corresponding to the inlet temperature of 106606.4 Pa. The cavitation number of the inflow was calculated by Equation (22): σ_∞=0.616.

Observing the figure, it becomes evident that the lift and drag coefficients of the hydrofoil exhibit distinct periodic characteristics, irrespective of whether thermodynamic effects are taken into account. Upon considering thermodynamic effects, both the lift coefficient (C_l) and the drag coefficient (C_d) exhibit increased amplitudes. The fluctuation amplitudes of both coefficients are notably higher in the isothermal cavitation case. However, after considering the thermodynamic effects, the periodic variation of the two coefficients become more complicated. In the subsequent sections, a detailed analysis will be provided to elucidate the underlying factors contributing to the complex periodic variations of thermodynamic effects of cavitation.

Figure 10. The variation of lift and drag coefficients of the NACA0015 hydrofoil with time under isothermal and non-isothermal conditions.

Figure 10. The variation of lift and drag coefficients of the NACA0015 hydrofoil with time under isothermal and non-isothermal conditions.

3.1. Analysis of thermodynamic effects on cavitation bubble growth and vortex shedding characteristics

3.1.1. Isothermal cavitation

From the time-domain plot of the drag coefficient, it is evident that the cavitation flow around the NACA0015 hydrofoil exhibits significant periodic characteristics. The initial time is designated as t₀, and subsequent times increase in Δt/12 increments. As depicted in Figure 11, at time t₀, there is an inverse pressure gradient at the trailing edge of the hydrofoil that results in fluid backflow. As the inverse pressure gradient propagates upstream, it induces small-scale vortices. By t₀+2Δt/12, the inverse pressure gradient reaches the leading edge of the hydrofoil, causing the cavitation region to detach as a whole. During this detachment process, small-scale vortices combine into larger-scale vortices, which subsequently shed from the trailing edge due to the influence of cavitation bubble clusters. Starting from t₀+6Δt/12, new attached cavities began to appear on the hydrofoil surface, and these cavities continued to grow until the upstream inverse pressure gradient urged them to separate again.

Figure 11. Instantaneous gas phase volume cloud map obtained from DDES simulation, time interval of 1ms.

Figure 11. Instantaneous gas phase volume cloud map obtained from DDES simulation, time interval of 1ms.

Figure 12. Instantaneous vorticity cloud map obtained from DDES simulation, time interval of 1ms.

Figure 12. Instantaneous vorticity cloud map obtained from DDES simulation, time interval of 1ms.

3.1.2. Thermodynamic effect cavitation

When accounting for the thermodynamic effect, the cavitation effect noticeably diminishes under the same inflow cavitation number. The thermodynamic effect assumes a significant inhibitory role in cryogenic cavitation. Owing to the weakened cavitation effect, the cavitation region displays a distinctive "mist-like" structure that differs from isothermal cavitation. With the incorporation of the thermodynamic effect, the extent of the inverse pressure gradient is limited to the middle section of the hydrofoil, resulting in the detachment of the cavitation region commencing from the middle. In a similar vein, vorticity originates from the middle part of the hydrofoil. Initially, due to the disturbance of inverse pressure gradient, small-scale eddy current appeared. Subsequently, as the vortices progress downstream, they amalgamate and consolidate into large-scale vortices due to the interaction with cavitation bubbles, eventually shedding at the trailing edge.

Figure 13. Instantaneous gas phase volume cloud map obtained from DDES simulation, with a time interval of 1ms.

Figure 13. Instantaneous gas phase volume cloud map obtained from DDES simulation, with a time interval of 1ms.

Figure 14. Instantaneous vorticity cloud map obtained from DDES simulation, with a time interval of 1ms.

Figure 14. Instantaneous vorticity cloud map obtained from DDES simulation, with a time interval of 1ms.

3.2. Influence of thermodynamic effects on the pressure field and velocity field

3.2.1. Isothermal cavitation

From Figure 15, it is evident that at time t₀, the entire suction surface of the hydrofoil is enveloped by the cavitation region, with the upper part of the surface completely submerged within the low-pressure region. At t₀+2Δt/12, the low-pressure region begins to rupture, causing the cavities on the upper part of the suction surface to fragment. Following the fragmentation, small cavitation clusters move downstream with the main flow, continuously coalescing. By t₀+6Δt/12, the small cavitation clusters amalgamate to form a large cavitation cluster that eventually sheds at the hydrofoil's trailing edge. Subsequently, the low-pressure region reappears on the suction surface, initiating a new cycle of evolution until the primary cavitation region once again engulfs the entire surface.

Figure 15. Instantaneous pressure distribution cloud map obtained from DDES simulation, with a time interval of 1ms.

Figure 15. Instantaneous pressure distribution cloud map obtained from DDES simulation, with a time interval of 1ms.

Based on existing research [27,28], the evolution of hydrofoil cavitation is closely linked to the phenomenon of re-entrant jet, and the progression of re-entrant jet is highly influenced by the range of inverse pressure gradients. In the subsequent analysis, the evolution of re-entrant jet within cavitation flow will be assessed in conjunction with velocity vectors. As depicted in Figure 16, the velocity vector field at the hydrofoil's trailing edge at t₀ is presented. The incoming flow from the distant field bifurcates into the main flows along the suction surface and pressure surface upon traversing the hydrofoil's leading edge. These two main flows converge at the hydrofoil's trailing edge, forming a localized zone of high pressure and low velocity. This convergence region serves as the catalyst for re-entrant jet. Following the occurrence of the re-entrant jet, it is propelled by the inverse pressure gradient, entering the lower-pressure cavity from the high-pressure convergence region, and progresses upstream.

Figure 16. Velocity vector distribution at the trailing edge of the hydrofoil at time t_0.

Figure 16. Velocity vector distribution at the trailing edge of the hydrofoil at time t_0.

Figure 17 illustrates the distribution of velocity vectors along the suction surface of the hydrofoil at t₀+3Δt/12. At this specific moment, the re-entrant jet originating from the high-pressure region at the hydrofoil's trailing edge propagates towards the leading edge, triggering the overall detachment of the primary cavitation region. From the figure, it can be observed that as the re-entrant jet infiltrates the primary cavitation region and progresses upstream, it interacts with the mainstream flow on the suction surface, resulting in the formation of a shear layer characterized by substantial velocity gradients. This shear layer, in turn, contributes to the formation of small-scale and large-scale vortex structures. These vortices serve as the primary mechanism driving the fragmentation of the primary cavitation region.

Figure 17. Velocity vector distribution at the middle section of the hydrofoil at time t₀+3Δt/12.

Figure 17. Velocity vector distribution at the middle section of the hydrofoil at time t₀+3Δt/12.

3.2.2. Thermodynamic effect cavitation

This statement highlights the impact of thermodynamic effects on the hydrofoil's evolution cycle. As illustrated in the Figure 18, the front 1/3 of the hydrofoil's suction surface consistently experiences coverage by a low-pressure region. Beyond this low-pressure region at the leading edge, large-scale fragmented areas of low pressure emerge, which are intricately tied to the evolution of vortices.

Figure 18. Transient pressure distribution contour map obtained from DDES simulation, with a time interval of 1ms.

Figure 18. Transient pressure distribution contour map obtained from DDES simulation, with a time interval of 1ms.

Figure 19 presents the distribution of velocity vectors on the suction surface of the hydrofoil at time t₀, considering thermodynamic effects. Upon examination of the figure, it becomes apparent that when thermodynamic effects are taken into account, the adverse pressure gradient can only develop up to a specific point in the middle section of the hydrofoil, whereas in isothermal cavitation, it extends all the way to the leading edge. As the adverse pressure gradient serves as the driving force for the re-entrant jet, its inability to develop up to the leading edge under thermodynamic effects prevents the overall detachment of the main cavitation region. Instead, it can only progress up to the middle section of the hydrofoil, leading to partial detachment of the main cavitation region.

When thermodynamic effects are considered, the upstream development of the re-entrant jet exhibits distinct characteristics compared to isothermal cavitation. Under the influence of thermodynamic effects, at the same free stream cavitation number, the intensity of cavitation is significantly suppressed, resulting in increased liquid-phase content within the main cavitation region. Moreover, the re-entrant jet encounters higher resistance during its upstream development. In the figure, it can be observed that the interaction between the re-entrant jet and the mainstream flow on the suction surface, leading to the formation of large-scale vortices, takes place when considering thermodynamic effects. These large-scale vortices then interact with the upstream mainstream flow on the suction surface, consequently giving rise to the formation of localized re-entrant jets.

Figure 19. Velocity vector distribution at the middle section of the hydrofoil at time t_0.

Figure 19. Velocity vector distribution at the middle section of the hydrofoil at time t_0.

4. Conclusions

In this study, the Sauer-Schnerr cavitation model was modified by adjusting the bubble number density to improve the agreement with experimental results conducted by Hord. Based on this modification, numerical simulations were carried out to investigate the influence of thermodynamic effects on the unsteady cavitation flow around a two-dimensional NACA0015 hydrofoil in liquid nitrogen.

• Through comparison and analysis of the numerical simulations with the experimental results, it was determined that the best agreement between the simulation and experimental data was achieved when the bubble number density in the Sauer-Schnerr cavitation model was set to 10⁸.
• The validated numerical model was then employed to simulate the cavitation flow around the NACA0015 hydrofoil in liquid nitrogen without considering thermodynamic effects. It was observed that the small-scale vortices induced by the upstream development of the re-entrant jet were the primary cause of fragmentation within the main cavitation region. The shedding motion of the bubbles contributed to the integration of these vortices.
• Subsequently, the validated numerical model was utilized to simulate the cavitation flow around the NACA0015 hydrofoil in liquid nitrogen, this time taking into account thermodynamic effects. It was observed that the cavitation effect was significantly diminished, resulting in a "mist-like" structure of the cavitation region. The incomplete development of the re-entrant jet upstream was identified as the fundamental reason for the inability of the cavitation cloud to detach as a whole.
• Irrespective of whether thermodynamic effects were considered or not, the re-entrant jet originated from the high-pressure, low-velocity region formed by the interaction of mainstream flows at the hydrofoil's trailing edge. When thermodynamic effects were incorporated, the upstream development of the re-entrant jet faced greater obstacles, leading to the formation of larger-scale vortices compared to isothermal cavitation.