A numerical method for variable temperature cavitation considering turbulence correction and energy-phase coupling

Abstract

The present application provides a kind of temperature variation cavitation numerical calculation method considering turbulence correction and energy-phase coupling, it is related to cavitation numerical calculation technical field, the present application obtains fluid unit temperature, current pressure, mixed phase density, flow velocity vector, turbulent kinetic energy and turbulent dissipation rate;Static saturated vapor pressure is established by fluid unit temperature, combined with the compensation of the generation of dynamic critical vaporization pressure of turbulence pressure fluctuation characteristics;Based on turbulent kinetic energy and turbulent dissipation rate, establish the reference item of moving turbulent viscosity, perform product substitution operation on it using density attenuation characteristics to generate local effective turbulent viscosity;According to the deviation of current pressure and dynamic critical vaporization pressure, establish phase change driven pressure difference, output phase change mass transfer rate and convert into phase change energy source term;Subsection analysis temperature evolution sensitivity, generate energy damping factor by mixed phase bulk heat capacity reduction process;Combined with flow velocity vector, joint iteration is solved and flow field distribution data is output, so as to improve the precision of temperature variation cavitation numerical calculation.

Description

Technical Field

[0001] This invention relates to the field of cavitation numerical calculation technology, specifically to a variable-temperature cavitation numerical calculation method that considers turbulence correction and energy-phase transition coupling. Background Technology

[0002] As various hydraulic machines evolve towards higher speeds and wider temperature ranges, fluids in low-pressure regions are prone to phase changes, leading to temperature-variable cavitation, which in turn causes vibration and performance degradation. High-precision numerical calculations of this cavitation flow field are a prerequisite for optimizing the structural design of hydraulic machinery. Traditional isothermal cavitation numerical methods often assume a constant flow field temperature, neglecting the localized temperature drop caused by the heat absorption during fluid vaporization. This approach suffers from a serious lack of understanding of the underlying physical mechanisms, making it difficult to accurately predict thermodynamic effects.

[0003] Existing numerical methods for variable-temperature cavitation based on computational fluid dynamics typically compare fixed theoretical parameters with averaged variables of the flow field to delineate phase boundaries. However, when dealing with complex flow fields containing strong turbulence, the existing numerical solution architecture for variable-temperature cavitation is severely inadequate in capturing transient fluctuations in strong turbulence, leading to passivation of the predicted vaporization critical point. Furthermore, due to the limitations of fixed basic phase transition coefficients, it is often necessary to supplement these with a large number of empirical parameters for manual optimization, using these as a benchmark to extrapolate cavitation morphology. Overall, these variable-temperature cavitation calculation methods mainly rely on a single turbulent viscosity assumption and a synergistic path specified by static mass transfer parameters to achieve numerical calculations of cavitation.

[0004] In real-world variable-temperature cavitation environments, the violent collapse and formation of bubbles are often accompanied by extremely high interphase mass transfer rates, leading to localized temperature evolution and thermodynamic distortions dominated by latent heat release and absorption. However, existing numerical solutions struggle to balance the highly sensitive nonlinear relationship between phase change rates and temperature evolution, resulting in numerical singularities during the energy source term iteration. Furthermore, the vapor-liquid mixing region often exhibits variations in interfacial drag caused by density steps, frequently leading to abnormal expansion of drag terms in the solution matrix, making the numerical calculation process highly unstable. This results in severely distorted predictions, making it difficult to perform objective and stable numerical calculations of variable-temperature cavitation solely based on the fundamental governing equations.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a numerical calculation method for variable-temperature cavitation that considers turbulence correction and energy-phase transition coupling, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling, the specific steps of which include: Step 1: Obtain the fluid element temperature, current pressure, mixed phase density, velocity vector, turbulent kinetic energy, turbulent dissipation rate, volume fraction of each phase, and element characteristic scale of each fluid element within the preset grid at the current moment; Step 2: Establish the static saturated vapor pressure based on the thermal state mapping between the fluid unit temperature and the preset physical property constants, construct the turbulent pressure pulsation term by cross-modulation of the mixed phase density and turbulent kinetic energy, and perform dynamic compensation on the static saturated vapor pressure to generate the dynamic critical vaporization pressure; Step 3: Construct a density decay function based on the evolution characteristics of mixed phase density and preset pure phase density, establish a motion turbulent viscosity reference term based on turbulent kinetic energy and turbulent dissipation rate, and perform a product substitution operation on the motion turbulent viscosity reference term using the density decay function to generate local effective turbulent viscosity. Step 4: Based on the deviation between the current pressure and the dynamic critical vaporization pressure, establish the phase change driving pressure difference, use turbulent kinetic energy to perform mass transfer amplification on the preset basic phase change coefficient, combine the volume fraction of each phase and the characteristic scale of the unit to perform phase state logic branch determination, and output the phase change mass transfer rate. Step 5: Based on the preset latent heat of vaporization, the phase change mass transfer rate is converted into a phase change energy source term, and the volumetric heat capacity of the mixed phase is established in combination with the volume fraction of each phase. The temperature evolution sensitivity of the phase change energy source term is analyzed in segments, and the temperature evolution sensitivity is reduced by using the volumetric heat capacity of the mixed phase to generate an energy damping factor. Step 6: Inject the local effective turbulent viscosity and phase change mass transport rate into the momentum conservation equation and cavitation transport equation, respectively. Combine the velocity vector to establish the spatial convection relationship. Inject the phase change energy source term into the energy conservation equation and use the energy damping factor to implement diagonal constraints. Solve the equations iteratively and output the flow field distribution data.

[0008] Furthermore, the volume fraction of each phase consists of the volume fraction of the liquid phase and the volume fraction of the vapor phase; The steps for constructing the preset mesh and obtaining the unit feature scale include: Based on the differences in physical flow regime within the target watershed, the watershed space is divided into a high-shear zone near the wall, a phase transition evolution zone, and a core mainstream zone; In the high shear zone near the wall, multiple layers of wall-attached mesh are arranged and the normal thickness of the first layer of wall-attached mesh is defined. In the phase transition evolution zone, a dense mesh with gradually decreasing size is arranged. In the core mainstream zone, a large-scale space-filling mesh is arranged, which together constitute the preset mesh. Extract the mesh volume of each fluid element in the preset mesh, and perform three-dimensional spatial dimensionality reduction feature extraction on the mesh volume to generate the element feature scale corresponding to each fluid element.

[0009] Furthermore, the preset physical property constant is composed of a preset first constant, a preset second constant, a preset third constant, and a preset reference pressure; The specific steps for generating the dynamic critical vaporization pressure are as follows: The temperature of the fluid element is summed with a preset third constant to construct a temperature offset denominator. The preset second constant is divided by the temperature offset denominator to obtain a nonlinear curvature offset term. The difference between the preset first constant and the nonlinear curvature offset term is calculated to construct a thermal state index term. Using the numerical value of 10 as the base and the thermal state index as the exponent, an exponential mapping operation is performed to obtain a reference mapping value. The reference mapping value, the preset reference pressure, and the numerical value 133.322 are then multiplied together to obtain the static saturated vapor pressure. The preset fluctuation coefficient, mixed phase density, and turbulent kinetic energy are multiplied together to construct the turbulent pressure fluctuation term; The turbulent pressure fluctuation term is used as a dynamic compensation increment and superimposed with the static saturated vapor pressure to generate the dynamic critical vaporization pressure. The specific calculation formula is as follows: In the formula, This is the corrected dynamic critical vaporization pressure. This is the static saturated vapor pressure. To preset the reference pressure, These are the first preset constant, the second preset constant, and the third preset constant, respectively. For the fluid element temperature, To preset the pulsation coefficient, For the mixed phase density, It is turbulent kinetic energy.

[0010] Furthermore, the preset pure phase density is composed of a preset pure vapor phase density and a preset pure liquid phase density; The specific steps for generating the local effective turbulent viscosity are as follows: The difference between the mixed phase density and the preset pure vapor phase density is calculated and calibrated as the first difference, and the difference between the preset pure liquid phase density and the preset pure vapor phase density is calculated and calibrated as the second difference. The first difference is divided by the second difference to obtain the phase evolution benchmark term. Using the phase evolution benchmark term as the base and the preset decay exponent as the exponent, a power function mapping operation is performed to construct a density decay function; Divide the square of the turbulent kinetic energy by the turbulent dissipation rate to obtain the turbulent micro-characteristic ratio. Perform a multiplication operation on the turbulent micro-characteristic ratio and the preset turbulence constant to establish the kinematic turbulent viscosity reference term. The density decay function is multiplied and substituted with the kinematic turbulent viscosity reference term to generate the local effective turbulent viscosity. The specific calculation formula is as follows: In the formula, For the local effective turbulent viscosity, For the mixed phase density, To preset the pure vapor phase density, To preset the pure liquid phase density, To preset the attenuation index, As a preset turbulence constant, For turbulent kinetic energy, This represents the turbulent dissipation rate.

[0011] Furthermore, the preset basic phase change coefficient is composed of a preset basic evaporation coefficient and a preset basic condensation coefficient; The specific steps for establishing the phase change driving pressure difference based on the deviation between the current pressure and the dynamic critical vaporization pressure, and using turbulent kinetic energy to perform mass transfer amplification on the preset basic phase change coefficient, include: Extract the preset pure liquid phase density from the preset pure phase density; The absolute difference between the dynamic critical vaporization pressure and the current pressure is calculated to establish the phase change driving pressure difference. The phase change driving pressure difference is divided by the preset pure liquid phase density and multiplied by two-thirds of the value, and then the square root operation is performed to generate the bubble dynamics reference term. The sum of the turbulent kinetic energy and the preset anti-divergence reference kinetic energy is used as the denominator, and the turbulent kinetic energy is used as the numerator to perform a ratio calculation. The result of the ratio calculation is multiplied by the preset mass transfer amplification factor and then superimposed with a value of one to construct the mass transfer amplification factor. The preset basic evaporation coefficient and the mass transfer amplification factor are multiplied and modulated to generate the corrected evaporation coefficient. The specific calculation formula is as follows: In the formula, To correct the evaporation coefficient, To preset the basic evaporation coefficient, To preset the mass transfer amplification factor, For turbulent kinetic energy, The preset reference kinetic energy is used to prevent divergence.

[0012] Furthermore, the volume fraction of each phase includes the volume fraction of the vapor phase; The specific steps for calculating the phase change mass transfer rate are as follows: The vapor phase volume fraction is extracted from the volume fraction of each phase, and the preset pure vapor phase density is extracted from the preset pure phase density; Based on the magnitudes of the current pressure and the dynamic critical vaporization pressure, a piecewise solution is performed: When the current pressure is lower than the dynamic critical vaporization pressure and the vapor phase volume fraction is less than the preset vapor phase cutoff threshold, the difference between value one and the vapor phase volume fraction is calculated. This difference is then multiplied by the preset condensation nucleus volume fraction, the preset pure vapor phase density, and value three. The resulting product is then divided by the unit characteristic scale to obtain the evaporation interface geometric mapping term. The corrected evaporation coefficient, the evaporation interface geometric mapping term, and the bubble dynamics reference term are then multiplied to calculate the positive phase change mass transfer rate. When the current pressure is higher than the dynamic critical vaporization pressure and the vapor phase volume fraction is greater than zero, the vapor phase volume fraction, the preset pure vapor phase density and the numerical value are multiplied together, and the resulting product is divided by the unit feature scale to obtain the condensation interface geometric mapping term. The preset basic condensation coefficient, the condensation interface geometric mapping term and the bubble dynamics reference term are multiplied together and the negative number is taken to calculate the negative phase change mass transfer rate. When any of the above phase constraint conditions are not met, the fluid unit is determined to be in the phase stability region, and the phase change mass transfer rate is set to zero. The specific calculation formula is as follows: In the formula, For phase change mass transfer rate, To correct the evaporation coefficient, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. To preset the vapor phase cutoff threshold, This refers to the preset pure vapor phase density within the preset pure phase density. For unit feature scale, The phase change driving pressure difference This refers to the preset pure liquid phase density within the preset pure phase density. It is the dynamic critical vaporization pressure. Due to current pressure, This is the preset basic condensation coefficient.

[0013] Furthermore, the volume fraction of each phase also includes the volume fraction of the liquid phase; The specific steps for converting the phase change mass transfer rate into a phase change energy source term based on a preset latent heat of vaporization, and for determining the mixed phase volumetric heat capacity in conjunction with the volume fraction of each phase, include: The preset latent heat of vaporization and the phase change mass transfer rate are multiplied, and the negative is taken to convert the phase change mass transfer rate into a phase change energy source term. The specific calculation formula is as follows: In the formula, For phase transition energy source term, To preset the latent heat of vaporization, This refers to the phase change mass transfer rate; The volume fraction of the liquid phase is extracted from the volume fraction of each phase. The liquid phase heat capacity term is obtained by performing a multiplication operation on the preset pure liquid phase constant pressure specific heat capacity, the preset pure liquid phase density, and the liquid phase volume fraction. The vapor phase heat capacity term is obtained by performing a multiplication operation on the preset pure vapor phase constant pressure specific heat capacity, the preset pure vapor phase density, and the vapor phase volume fraction. The liquid phase heat capacity term and the vapor phase heat capacity term are summed to establish the mixed phase volume heat capacity.

[0014] Furthermore, the preset physical property constants include a preset second constant and a preset third constant; The specific steps for generating the energy damping factor are as follows: Extract the preset second and third constants from the preset physical property constants; The fluid element temperature is summed with the preset third constant. The square of the summation result is used as the denominator of the partial derivative. The product of the preset second constant and the natural logarithm of the value of 10 is calculated and used as the numerator of the partial derivative. The numerator of the partial derivative is divided by the denominator of the partial derivative to obtain the derivative nonlinear coefficient. Then, it is multiplied with the static saturated vapor pressure to construct the absolute derivative term. Calculate the product of the phase change driving pressure difference and the preset pure liquid phase density and two-thirds of the value, and take the square root of the product and find its reciprocal as the velocity derivative bias term. Segmented analysis is performed based on the magnitude of the current pressure and the dynamic critical vaporization pressure: When the difference between the dynamic critical vaporization pressure and the current pressure is greater than the preset anti-divergence pressure difference threshold, and the vapor phase volume fraction is less than the preset vapor phase cutoff threshold, the preset latent heat of vaporization, the corrected evaporation coefficient, one-third of the values ​​of the evaporation interface geometric mapping term, the velocity derivative bias term, and the absolute derivative term are multiplied together and the negative is taken to analyze the temperature evolution sensitivity of the phase change energy source term under the evaporation condition. The temperature evolution sensitivity is divided by the mixed phase volume heat capacity to perform order reduction processing and generate the energy damping factor under the evaporation condition. When the difference between the current pressure and the dynamic critical vaporization pressure is greater than the preset anti-divergence pressure difference threshold, and the vapor phase volume fraction is greater than zero, the preset latent heat of vaporization, the preset basic condensation coefficient, one-third of the value of the condensation interface geometric mapping term, the velocity derivative bias term, and the absolute derivative term are multiplied together and the negative is taken to analyze the temperature evolution sensitivity of the phase change energy source term under condensation conditions. The temperature evolution sensitivity is divided by the mixed phase volume heat capacity to perform order reduction processing and generate the energy damping factor under condensation conditions. When any of the above operating conditions are not met, the energy damping factor is converged and assigned a value of zero. The specific calculation formula is as follows: In the formula, It is the energy damping factor. For temperature evolution sensitivity, To preset the latent heat of vaporization, To correct the evaporation coefficient, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. This refers to the preset pure vapor phase density within the preset pure phase density. For unit feature scale, This represents the liquid phase volume fraction within each phase volume fraction. This refers to the preset pure liquid phase density within the preset pure phase density. To preset the specific heat capacity of pure liquid phase at constant pressure, To preset the specific heat capacity of the pure vapor phase at constant pressure, The phase change driving pressure difference This is the static saturated vapor pressure. For the fluid element temperature, These are the second and third preset constants in the preset physical property constants. It is the dynamic critical vaporization pressure. Due to current pressure, To preset the vapor phase cutoff threshold, To preset the basic condensation coefficient, The preset anti-divergence pressure difference threshold is set.

[0015] Compared with the prior art, the beneficial effects of the present invention are: This invention generates a dynamic critical vaporization pressure based on the static saturated vapor pressure determined by the fluid element temperature, combined with compensation for turbulent pressure fluctuations. It then uses a kinematic turbulent viscosity benchmark term and a density decay function to perform a product substitution operation to generate a locally effective turbulent viscosity. This invention overcomes the limitations of traditional methods that rely excessively on simple comparisons of static saturated vapor pressure. By introducing density decay characteristics, it corrects the non-physical turbulent viscosity dissipation bias caused by density steps in the mixed phase. Through deep physical coupling of turbulent kinetic energy and turbulent dissipation rate with the critical threshold of phase transition, it realistically recreates the interphase shear and cavitation evolution state under strong turbulent conditions, significantly improving the accuracy of numerical calculations. This invention also establishes the phase change driving pressure difference based on the deviation between the current pressure and the dynamic critical vaporization pressure, and uses the mixed phase volumetric heat capacity to reduce the order of temperature evolution sensitivity to generate an energy damping factor. This mechanism solves the drawbacks of traditional fixed phase change coefficients and heavy reliance on empirical parameters for manual tuning, effectively resolving the highly sensitive coupling conflict between phase change mass transfer rate and temperature evolution. Finally, in the joint iterative solution process, the nonlinear numerical singularities generated by the energy source term are effectively suppressed, achieving high convergence and numerical stability in the numerical calculation of variable-temperature cavitation flow fields. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a schematic diagram illustrating the evolution of the response relationship between the phase change energy source term and the temperature of the fluid element under variable-temperature cavitation flow field conditions. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0018] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0019] Example: Please see Figures 1-2 The present invention provides a technical solution: A numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling, the specific steps of which include: Step 1: Obtain the fluid element temperature, current pressure, mixed phase density, velocity vector, turbulent kinetic energy, turbulent dissipation rate, volume fraction of each phase, and element characteristic scale of each fluid element within the preset grid at the current moment.

[0020] Cavitation phase transitions are often accompanied by the initial formation and collapse of bubbles, as well as the intense release and release of latent heat, resulting in extreme non-uniformity of the flow field. Discretizing using conventional uniform grids cannot accurately capture the turbulent viscous substructure near the wall and the thermodynamic distortions in the core region. Therefore, based on the differences in physical flow regimes within the target computational domain, the domain space is divided into a near-wall high-shear zone, a phase transition evolution zone, and a core-mainstream zone. In the near-wall high-shear zone, multiple layers of wall-attached grids are deployed with a defined normal thickness for the first layer. In the phase transition evolution zone, a gradually decreasing density grid is deployed, and in the core-mainstream zone, a large-scale space-filling grid is deployed. This division of the flow domain space further improves the accuracy of capturing local transient features. In this context, since the near-wall region where cavitation occurs is accompanied by extreme velocity shear gradients and temperature jumps, if the first-layer mesh is too large, the flow details of the innermost viscous sublayer will be lost. Therefore, the mesh size is gradually reduced according to the spatial span of the velocity shear gradient in the near-wall boundary layer until the phase transition characteristic parameters reach mesh-independent convergence. In this embodiment, the normal thickness is set to 0.01 mm to block the occurrence of non-physical numerical diffusion in the boundary layer.

[0021] To provide spatial parameters to eliminate interference from mesh heterogeneity, the mesh volume of each fluid element is extracted. Spatial dimensionality reduction feature extraction is performed by calculating the mesh volume to the power of one-third, generating the element feature scale corresponding to each fluid element. Simultaneously, the fluid element temperature, current pressure, mixed phase density, velocity vector, turbulent kinetic energy, turbulent dissipation rate, volume fraction of each phase, and element feature scale are acquired at the current moment. The volume fraction of each phase consists of the liquid phase volume fraction and the vapor phase volume fraction.

[0022] Step 2: Establish the static saturated vapor pressure based on the thermal state mapping between the fluid unit temperature and the preset physical property constants. Construct the turbulent pressure pulsation term by cross-modulating the mixed phase density and turbulent kinetic energy, and perform dynamic compensation on the static saturated vapor pressure to generate the dynamic critical vaporization pressure.

[0023] Cavitation processes in variable-temperature environments are accompanied by intense heat transfer, and local temperature differences in the flow field can lead to nonlinear drift of the vaporization phase transition boundary. Based on the thermal state mapping between the fluid element temperature and preset physical property constants (preset first constant, preset second constant, preset third constant, and preset reference pressure), the pure thermodynamic phase transition threshold at a specific temperature node can be accurately reproduced. To ensure that this pure thermodynamic phase transition threshold accurately matches common hydraulic machinery working media, this embodiment uses specific physical property constants derived from fitting the evaporation characteristics of real water media over a wide temperature range. Constants A, B, and C from the standard Antoine equation are extracted and used to replace the preset first constant, preset second constant, and preset third constant. The values ​​of the preset first constant are set to 8.07131, the preset second constant to 1730.63, and the preset third constant to -39.724 to avoid distortion of the phase transition boundary caused by pure mathematical extrapolation.

[0024] This embodiment uses the sum of the fluid unit temperature and a preset third constant to construct a temperature offset denominator, establishing a reference for the offset of the thermodynamic absolute temperature scale. The sum of the fluid unit temperature and the preset third constant only becomes zero when the absolute temperature of the fluid unit drops to 39.724 Kelvin (approximately -233 degrees Celsius). At this point, the water medium has already completely solidified into solid ice, so the sum can never be zero. This further establishes the offset reference for the absolute temperature scale and, based on the real-world physical laws of liquid phases, prevents the risk of a division-to-zero singularity. Subsequently, a preset second constant is divided by this temperature offset denominator to obtain a nonlinear curvature bias term, which is then subtracted from a preset first constant to construct a thermal state exponent term. This reduces the complexity of the fluid phase change latent heat evolution curve to an algebraic polynomial that can be executed by exponential operations. Mapping is performed using a numerical value of 10 as the base and the thermal state exponent as the exponent. Multiplication is then performed with the numerical value 133.322 and a preset reference pressure to finally obtain the static saturated vapor pressure. Wherein, the value 133.322 is the standard physical dimension conversion coefficient from millimeters of mercury to Pascals, further eliminating the distortion of the solution matrix caused by dimensional mismatch; to ensure the conservation of dimensions in algebraic operations, this embodiment sets the preset reference pressure to the value 1. The specific calculation formula for the static saturated vapor pressure is as follows: In the formula, This is the static saturated vapor pressure. To preset the reference pressure, These are the first preset constant, the second preset constant, and the third preset constant, respectively. The temperature of the fluid element.

[0025] The static saturated vapor pressure is used to characterize the absolute pressure threshold required for a fluid medium to undergo a phase change under purely thermodynamic conditions. A larger absolute pressure threshold indicates more intense molecular thermal motion within the fluid, making it easier for molecules to break free from intermolecular bonds and vaporize significantly. In this case, only a small pressure drop in the surrounding flow field is needed to induce cavitation. Conversely, a smaller absolute pressure threshold indicates a highly stable liquid phase, making vaporization difficult. The static saturated vapor pressure is dominated by the temperature of the independent fluid unit. Increased temperature imparts higher kinetic energy to the liquid molecules, causing the static saturated vapor pressure to rise nonlinearly. This embodiment utilizes a preset first constant, a preset second constant, and a preset third constant to construct a thermal state exponential term, accurately reproducing the nonlinear physical curvature of the latent heat evolution of the medium. Subsequently, a dimensionless basis is provided through exponential mapping of the numerical value of 10 and a preset reference pressure, transforming the empirical state into a real pressure dimension. Finally, a forced dimensional transformation from millimeters of mercury to Pascals is performed by multiplying by the numerical value of 133.322 to ensure absolute alignment of pressure magnitudes and completely eliminate the discrete solution distortion caused by dimensional mismatch.

[0026] Real fluid machinery often involves intense high-speed shearing and eddy evolution. Furthermore, the macroscopic time-averaged pressure calculation framework filters out transient low-pressure extremes within microscopic vortices, leading to significant hysteresis errors in predicting cavitation initiation in high-shear regions. Therefore, this embodiment performs a multiplication operation on a preset fluctuation coefficient, mixed phase density, and turbulent kinetic energy to construct a turbulent pressure fluctuation term. To ensure this dynamic pressure is reasonably converted into macroscopic pressure drop extremes, this embodiment sets the preset fluctuation coefficient to 0.195 based on the statistical principle of normal distribution extremes under the isotropic turbulence assumption, ensuring the objective physical basis for parameter selection. When the target calculation scenario is a strongly anisotropic turbulent field such as a strong swirling flow or a narrow-gap flow, a conventional linear proportional fine-tuning can be performed within a reasonable empirical range of 0.15 to 0.39 based on the estimated degree of anisotropy of the flow field. That is, the greater the field shear rate and the stronger the turbulent anisotropy, the more positively the coefficient should be increased to compensate for the transient pressure drop deeper in the vortex core. This multiplication step further quantifies the promoting effect of turbulent vortices on cavitation-induced cavitation. Subsequently, the turbulent pressure fluctuation term is used as a dynamic compensation increment and superimposed with the static saturated vapor pressure to generate the dynamic critical vaporization pressure. This allows for the sensitive detection of early cavitation initiation phenomena caused by microscopic vortices in high turbulence regions of the flow field, even before the macroscopic uniform pressure has dropped to the static vaporization point. The specific calculation formula for the dynamic critical vaporization pressure is as follows: In the formula, This is the corrected dynamic critical vaporization pressure. This is the static saturated vapor pressure. For the fluid element temperature, To preset the pulsation coefficient, For the mixed phase density, It is turbulent kinetic energy.

[0027] The dynamic critical vaporization pressure is used to characterize the actual phase change triggering boundary, which integrates macroscopic thermodynamics and microscopic fluid dynamics. A higher dynamic critical vaporization pressure indicates that the fluid is in an extremely unstable high-shear state, and cavitation can be induced earlier even under higher ambient pressure. Conversely, a lower dynamic critical vaporization pressure indicates a stable local flow pattern and difficulty in vaporization. The static saturated vapor pressure establishes the purely thermodynamic phase change bed, while the mixed phase density, turbulent kinetic energy, and a preset pulsation coefficient are multiplied to construct a turbulent pressure pulsation term. This turbulent pressure pulsation term essentially converts the mechanical kinetic energy pulsation within the microscopic vortex into an equivalent macroscopic pressure drop amplitude. An increase in the thermodynamic baseline or an intensification of turbulent pulsation will positively raise this triggering boundary. In real hydraulic machinery, the microscopic vortex center generated by high-speed shearing easily forms a transient low pressure, and macroscopic flow field calculations tend to smooth out this low-pressure extreme value, leading to a delay in phase change prediction. In this embodiment, the turbulent pressure fluctuation term is superimposed on the static saturated vapor pressure as a dynamic compensation increment, which reasonably raises the critical judgment threshold. This allows the flow field to keenly capture the early cavitation induced by microscopic turbulence even when the macroscopic uniform pressure has not dropped to the static vaporization point.

[0028] Step 3: Construct a density decay function based on the evolution characteristics of the mixed phase density and the preset pure phase density, establish a reference term for the motion turbulent viscosity based on turbulent kinetic energy and turbulent dissipation rate, and perform a product substitution operation on the reference term for the motion turbulent viscosity using the density decay function to generate the local effective turbulent viscosity.

[0029] In the numerical evolution of variable-temperature cavitation flow fields, classical turbulence equations often overestimate the turbulent viscosity of the vapor-liquid mixed phase region. The preset pure phase density consists of a preset pure vapor phase density and a preset pure liquid phase density, and according to real physical laws, the preset pure liquid phase density is always greater than the preset pure vapor phase density. In this embodiment, based on standard thermodynamic properties and actual ambient temperature, the preset pure liquid phase density is set to 998.2 kg / m³, and the preset pure vapor phase density is set to 0.0173 kg / m³ to eliminate non-physical viscosity decay distortion caused by arbitrary density boundary values. During the actual numerical iterative solution of the flow field, the preset pure liquid phase density and preset pure vapor phase density can be replaced with the dynamic pure phase density corresponding to the current fluid unit temperature based on the real-time fluid unit temperature and the actual fluid medium property table.

[0030] This embodiment calculates and labels the difference between the mixed phase density and the preset pure vapor phase density as a first difference, and calculates and labels the difference between the preset pure liquid phase density and the preset pure vapor phase density as a second difference. Dividing the first difference by the second difference yields a phase evolution benchmark term, which is used to forcibly map the complex mixture density range to a range of 0 to 1, thus eliminating the interference of differences in absolute density values ​​between different fluid media. Subsequently, a power function mapping operation is performed using the phase evolution benchmark term as the base and a preset decay index as the exponent to construct a density decay function. Specifically, to accurately suppress overdeveloped viscous substrates, this embodiment sets the preset decay index to 10 based on the empirically recommended value from classical computational fluid dynamics. This is used to reflect the extremely strong nonlinear curvature characteristics of the viscosity decay process in the region of intense phase change. The preset decay index can be adjusted according to the unsteady intensity of cloud cavitation shedding within the target watershed. Within the empirical range, routine trial and error adjustments should be made. That is, if the periodic shedding of large-scale cavitation clusters is not fully captured (i.e., insufficient viscosity suppression), the index should be appropriately increased positively to enhance the targeting resistance suppression effect.

[0031] Finally, the square of the turbulent kinetic energy is divided by the turbulent dissipation rate to obtain the turbulent micro-characteristic ratio, which further reflects the characteristic time and spatial scale of the local eddies. The turbulent micro-characteristic ratio and a preset turbulence constant are multiplied together to establish the kinematic turbulent viscosity reference term. In this embodiment, the preset turbulence constant is set to 0.09 by retrieving the standard universal constant from the classical fluid dynamics governing equations, accurately restoring the basic macroscopic momentum diffusivity before cavitation occurs. The density decay function and the kinematic turbulent viscosity reference term are multiplied and substituted to generate the local effective turbulent viscosity, as shown in the following formula: In the formula, For the local effective turbulent viscosity, For the mixed phase density, To preset the pure vapor phase density, To preset the pure liquid phase density, To preset the attenuation index, As a preset turbulence constant, For turbulent kinetic energy, This represents the turbulent dissipation rate.

[0032] The locally effective turbulent viscosity is used to characterize the true macroscopic momentum exchange capacity of a fluid undergoing a drastic phase transition, after deducting non-physical interfacial resistance. A higher locally effective turbulent viscosity indicates a pure phase state with strong momentum dissipation; a lower locally effective turbulent viscosity closely matches the actual physical phenomenon of extremely low viscosity inside cavitation bubbles and at interfaces. Specifically, the kinetic viscosity benchmark increases when turbulent kinetic energy increases or turbulent dissipation rate decreases. Simultaneously, as the mixed phase density decreases, its first difference from the preset pure vapor phase density decreases, and the obtained phase evolution benchmark decreases proportionally. This embodiment constructs a density decay function using the phase evolution benchmark and a preset decay index, and uses it to completely replace the mixed phase density in traditional standard viscosity calculations. This adaptively reduces artificially generated excessive viscous damping, implements targeted resistance suppression at interfaces, and ultimately ensures accurate capture of the transient stripping characteristics of cloud-like cavitation.

[0033] Step 4: Based on the deviation between the current pressure and the dynamic critical vaporization pressure, establish the phase change driving pressure difference, use turbulent kinetic energy to perform mass transfer amplification on the preset basic phase change coefficient, combine the volume fraction of each phase and the characteristic scale of the unit to perform phase state logic branch determination, and output the phase change mass transfer rate.

[0034] In this embodiment, the preset basic phase change coefficient is composed of a preset basic evaporation coefficient and a preset basic condensation coefficient. To accurately establish the basic rate boundary of phase change in a purely thermodynamic dimension, this embodiment sets the preset basic evaporation coefficient to 50 and the preset basic condensation coefficient to 0.01 based on the convergence experience of historical cavitation simulations of similar hydraulic machinery. In actual fluid engineering extensions, conventional forward or reverse proportional fine-tuning can be performed by comparing the actual macroscopic cavitation experimental morphology of the target hydraulic machinery under standard operating conditions. Since the volume expansion or contraction of cavitation bubbles is strictly controlled by the mechanical pressure imbalance on both sides of the bubble, this embodiment extracts a preset pure liquid phase density from the preset pure phase density, calculates the absolute difference between the dynamic critical vaporization pressure and the current pressure to establish the phase change driving pressure difference, divides the phase change driving pressure difference by the preset pure liquid phase density, multiplies it by two-thirds, and then performs a square root operation to generate a bubble dynamics benchmark term to eliminate secondary interferences such as surface tension and viscosity. Because the internal flow field of real fluid machinery is not static, strong turbulent shearing can cause bubbles to break up and split, leading to a geometric increase in the actual contact area between the vapor and liquid phases. Therefore, this embodiment uses the sum of the turbulent kinetic energy and the preset anti-divergence reference kinetic energy as the denominator, and the turbulent kinetic energy as the numerator to perform a ratio operation. Based on the principle of anti-singularity in numerical calculations, the preset anti-divergence reference kinetic energy is set to 0.001, meaning its absolute value is required to be much greater than the current hardware's machine truncation error and much less than one-thousandth of the macroscopic average turbulent kinetic energy of the flow field, to avoid non-physical interference with the macroscopic flow state. The obtained ratio is multiplied by a preset mass transfer amplification factor and then summed to construct the mass transfer amplification factor. This embodiment, based on hydraulic empirical calibration, sets the preset mass transfer amplification factor to 1.0. This allows for conventional adaptive scaling of the amplification factor, taking into account the actual shear tearing degree of the flow field, without causing oscillations in the overall numerical residuals. The macroscopic flow field turbulence intensity is converted into a gain weight for the microscopic phase change rate. Finally, the preset basic evaporation coefficient and the mass transfer amplification factor are multiplied and modulated to generate a corrected evaporation coefficient. The specific calculation formula is as follows: In the formula, To correct the evaporation coefficient, To preset the basic evaporation coefficient, To preset the mass transfer amplification factor, For turbulent kinetic energy, The preset reference kinetic energy is used to prevent divergence.

[0035] The modified evaporation coefficient is used to characterize the actual evaporation intensity benchmark after superimposing microscopic turbulent shear effects. A larger modified evaporation coefficient indicates that the flow field is subjected to extremely strong shear, causing severe breakup of the liquid phase, rapid expansion of the microscopic bubble group, and highly active gas-liquid mass transfer; a smaller modified evaporation coefficient indicates a stable local flow regime. Turbulent kinetic energy characterizes the microscopic mechanical strength of local eddies; a higher value indicates a stronger breaking effect on the fluid medium. In this embodiment, the sum of turbulent kinetic energy and a very small preset anti-divergence reference kinetic energy is used as the denominator, and turbulent kinetic energy is used as the numerator to perform a ratio calculation. The resulting ratio is multiplied by a preset mass transfer amplification factor, and then a numerical value is superimposed to construct a mass transfer amplification factor. Finally, a product modulation operation is performed between the preset basic evaporation coefficient providing a purely thermodynamic bed and this factor. The modified evaporation coefficient is nonlinearly positively correlated with turbulent kinetic energy. Because real strong turbulent shear causes bubble breakup and splitting, resulting in a geometric increase in the contact surface area between the two phases, conventional equations often underestimate this phase change rate. This embodiment objectively quantifies the gain weight of the flow field on microscopic mass transfer through turbulent kinetic energy, and cleverly utilizes the preset anti-divergence reference kinetic energy at the denominator to ensure the convergence stability of the underlying calculation under high shear boundary conditions.

[0036] In this embodiment, to provide a high-fidelity vapor-liquid mass exchange source term, this embodiment extracts the initial saturated dissolved gas content of the target fluid under the current ambient temperature and pressure according to Henry's Law and performs equivalent physical conversion, thereby setting the preset condensation nucleus volume fraction. In this embodiment, the preset condensation nucleus volume fraction is set to 0.0005. This embodiment uses spatial mass conservation extremum derivation of the multiphase flow continuity equation and combines it with rounding errors caused by grid discretization to set the range of the preset vapor phase cutoff threshold to... In this embodiment, the preset vapor phase cutoff threshold is set to 0.99. This embodiment extracts the vapor phase volume fraction from the volume fractions of each phase and extracts the preset pure vapor phase density from the preset pure phase density, performing piecewise calculations based on the values ​​of the current pressure and the dynamic critical vaporization pressure. When the current pressure is lower than the dynamic critical vaporization pressure and the vapor phase volume fraction is less than the preset vapor phase cutoff threshold, the difference between value one and the vapor phase volume fraction is calculated. This difference is then multiplied by the preset condensation nucleus volume fraction, the preset pure vapor phase density, and value three. The resulting product is then divided by the unit feature scale to obtain the evaporation interface geometric mapping term. The specific calculation formula is as follows: In the formula, For the geometric mapping term of the evaporation interface, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. This refers to the preset pure vapor phase density within the preset pure phase density. The unit feature scale.

[0037] The evaporation interface geometric mapping term characterizes the effective physical interface size that allows for gas-liquid contact and mass transfer during the evaporation stage. A larger value indicates a large specific surface area of ​​microbubbles and abundant mother liquor, resulting in highly active gas-liquid mass transfer; a smaller value indicates a scarcity of liquid phase available for vaporization and limited mass transfer. This formula quantifies the specific surface area characteristics of microscopic spherical cavities by using the quotient of the numerical value and the unit characteristic scale. A smaller unit characteristic scale indicates a geometrically larger phase change contact area. The difference between the numerical value and the vapor phase volume fraction characterizes the proportion of currently unvaporized usable liquid phase space, which, along with the preset condensation nucleus volume fraction, positively constrains the initial nucleation scale of the phase change. Finally, a multiplication operation is performed with the preset pure vapor phase density to further achieve precise separation of the actual effective evaporation area boundary.

[0038] The modified evaporation coefficient, the evaporation interface geometric mapping term, and the bubble dynamics reference term are multiplied together to calculate the positive phase change mass transfer rate, which reflects the physical reality that the less liquid space there is, the scarcer the mother liquor available for evaporation.

[0039] When the current pressure is higher than the dynamic critical vaporization pressure and the vapor phase volume fraction is greater than zero, the vapor phase volume fraction, the preset pure vapor phase density, and the numerical value are multiplied together, and the resulting product is divided by the unit feature scale to obtain the condensation interface geometric mapping term. The specific calculation formula is as follows: In the formula, For the geometry mapping term of the condensation interface, This represents the vapor phase volume fraction within the total volume fraction of each phase. This refers to the preset pure vapor phase density within the preset pure phase density. The unit feature scale.

[0040] The condensation interface geometric mapping term characterizes the effective physical contact interface size available for vapor condensation during the collapse phase of the fluid. A larger value indicates a large local cavitation surface area and highly active condensation mass transfer; a smaller value indicates that the bubbles gradually disappear and mass transfer tends to stagnate. This formula accurately quantifies the specific surface area of ​​microscopic spherical bubbles by using the quotient of the numerical value and the unit characteristic scale. A smaller unit characteristic scale indicates a geometrically enlarged condensation contact surface. The vapor phase volume fraction characterizes the actual residual vapor percentage in the current space, positively dominating the absolute base available for condensation. Finally, it is multiplied by the preset pure vapor phase density to complete the physical dimensions required for mass transfer. This formula restores the true contraction characteristics of the effective condensation interface nonlinearly decaying as the bubble volume shrinks, ensuring that negative mass transfer converges smoothly with the bubble's collapse.

[0041] The preset basic condensation coefficient, the condensation interface geometric mapping term, and the bubble dynamics benchmark term are multiplied together and their opposites are taken to calculate the negative phase change mass transfer rate, ensuring that collapse condensation only occurs when cavitation bubbles actually exist in physical space.

[0042] When any of the above phase constraint conditions are not met, the fluid unit is determined to be in the phase stability range, and the phase change mass transfer rate is assigned to zero, further eliminating false mass transfer caused by numerical oscillations. The specific calculation formula is as follows: In the formula, For phase change mass transfer rate, To correct the evaporation coefficient, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. To preset the vapor phase cutoff threshold, This refers to the preset pure vapor phase density within the preset pure phase density. For unit feature scale, The phase change driving pressure difference This refers to the preset pure liquid phase density within the preset pure phase density. It is the dynamic critical vaporization pressure. Due to current pressure, This is the preset basic condensation coefficient.

[0043] The phase change mass transfer rate is used to characterize the absolute physical exchange mass of liquid vaporization or vapor condensation per unit time and volume. A larger absolute value of the phase change mass transfer rate indicates more intense cavitation expansion or collapse within the flow field, and highly active gas-liquid mass transfer; a smaller absolute value indicates a more stable local medium phase state. The relative magnitude of the current pressure and the dynamic critical vaporization pressure defines the mass transfer direction. The bubble dynamics benchmark term, constructed from the phase change driving pressure difference and the preset pure liquid phase density, determines the radial motion velocity of the bubble wall, which is positively correlated with the transfer rate. The unit characteristic scale equivalently defines the specific surface area of ​​the two-phase contact, which is inversely correlated with it. Under evaporation conditions, the evaporation interface geometric mapping term is constructed from the difference between the numerical value and the vapor phase volume fraction, the preset condensation nucleus volume fraction, and the preset pure vapor phase density, and then multiplied together with the corrected evaporation coefficient. Under condensation conditions, the condensation interface geometric mapping term is constructed from the vapor phase volume fraction and the preset pure vapor phase density, and multiplied together with the preset basic condensation coefficient. This embodiment further restores the phase change dynamics law driven by pressure difference, and introduces a preset vapor phase cutoff threshold to perform space capacity limit interception, completely blocking the non-physical divergent solution that still occurs after the liquid is exhausted and then evaporates beyond the limit.

[0044] Step 5: Based on the preset latent heat of vaporization, the phase change mass transfer rate is converted into a phase change energy source term, and the volumetric heat capacity of the mixed phase is established in combination with the volume fraction of each phase. The temperature evolution sensitivity of the phase change energy source term is analyzed in segments, and the temperature evolution sensitivity is reduced in order using the volumetric heat capacity of the mixed phase to generate an energy damping factor.

[0045] In this embodiment, variable-temperature cavitation not only manifests as a change in the spatial volume of the vapor and liquid phases, but is also accompanied by intense latent heat release. This heat exchange reverses the local temperature field of the fluid. To accurately construct a cross-dimensional bridge between mass conservation and energy conservation, this embodiment multiplies the preset latent heat of vaporization with the aforementioned acquired phase change mass transfer rate, and takes the negative value to convert it into a phase change energy source term. This is used to equivalently convert simple hydrodynamic mass transfer into a heat source or heat sink in the thermodynamic dimension; taking the negative value is to conform to the actual physical thermodynamic direction, that is, when liquid evaporation generates positive mass transfer, the fluid particles will absorb heat from the outside, and the phase change energy source term will be negative, thus leading to local cooling; conversely, when cavitation bubbles condense and collapse, they release latent heat, leading to local heating. The specific calculation formula for the phase change energy source term is as follows: In the formula, For phase transition energy source term, To preset the latent heat of vaporization, This refers to the phase change mass transfer rate.

[0046] The phase change energy source term characterizes the absolute total amount of heat generated or absorbed by the fluid due to vapor-liquid phase change per unit time and unit volume. A larger absolute value of the phase change energy source term indicates more intense energy transfer and thermodynamic state reversal in the local flow field; a smaller absolute value means a relatively stable local phase state with no significant heat exchange. The preset latent heat of vaporization characterizes the inherent energy required for phase change per unit mass of medium, determining the energy baseline scale under equal mass transfer. The phase change mass transfer rate characterizes the actual phase change flow rate, positively controlling the output scale of the energy source term. The phase change energy source term uses the product of the preset latent heat of vaporization and the phase change mass transfer rate because classical thermodynamic laws define the linear conduction nature of latent heat of phase change. Furthermore, the negative sign before the product is to strictly match the actual heat transfer direction in the physical environment: when liquid evaporation generates positive mass transfer, fluid particles absorb heat, causing the source term to be negative and thus driving local cooling; conversely, when cavitation bubbles collapse and condense, they release latent heat, driving local heating. The introduction of a negative sign ensures that the formula can correctly map the objective thermodynamic laws of evaporation endothermy and condensation exothermy.

[0047] Since the fluid elements within the flow field are often in a discrete state of vapor-liquid coexistence, the thermophysical properties of the pure phase medium cannot truly reflect the transient heat storage capacity of the mixing region. Therefore, this embodiment extracts the liquid phase volume fraction from each phase volume fraction, and multiplies the preset pure liquid phase constant-pressure specific heat capacity, preset pure liquid phase density, and liquid phase volume fraction to obtain the liquid phase heat capacity term; similarly, multiplying the preset pure vapor phase constant-pressure specific heat capacity, preset pure vapor phase density, and vapor phase volume fraction together to obtain the vapor phase heat capacity term; finally, the sum of these two terms establishes the mixed phase volumetric heat capacity. This is used to quantify the true thermal inertia of the current fluid element at a specific vapor-liquid ratio. The larger the mixed phase volumetric heat capacity, the smaller the temperature fluctuation caused by the absorption of an equal amount of latent heat in that region; conversely, in the vapor-rich region, even a small heat release can trigger a drastic temperature jump. Specifically, for the target working medium, standard property scatter data for the corresponding temperature range are extracted, and a property call curve dynamically changing with temperature is generated through high-order polynomial fitting; during the numerical iterative solution process, the fitted curve is interpolated in real time based on the current fluid element temperature. This embodiment uses pure water as a reference under specific operating conditions, and sets the preset latent heat of vaporization value as follows: The preset value of the constant pressure specific heat capacity of the pure liquid phase is set to 4182, and the preset value of the constant pressure specific heat capacity of the pure vapor phase is set to 1918.

[0048] In high-dimensional numerical solutions of non-isothermal flow fields, the phase change energy source term is dominated by temperature nonlinearity. Explicit integration can easily lead to non-physical oscillations and divergences in the temperature variable. This embodiment extracts a preset second constant and a preset third constant. The sum and square of the fluid element temperature and the preset third constant are used as the denominator of the partial derivative. The product of the preset second constant and the natural logarithm of the value is used as the numerator of the partial derivative. Dividing the two yields the derivative nonlinearity coefficient, which is then multiplied by the static saturated vapor pressure to construct the absolute derivative term. Subsequently, the reciprocal of the square root of the product of the phase change driving pressure difference and the preset pure liquid phase density and two-thirds of the value is calculated as the velocity derivative bias term, used to reflect the sensitivity of the microbubble wall interface motion velocity to the local pressure difference.

[0049] To prevent division singularities and non-physical derivative oscillations from being induced in the square root operation of the denominator of the velocity derivative bias term when the minimum pressure difference approaches zero, this embodiment evaluates the truncation error level of the double-precision floating-point operation used by the current fluid computing platform and sets the threshold to a small positive number slightly larger than the local pressure difference extreme value corresponding to the truncation error. In this embodiment, the preset anti-divergence pressure difference threshold is set to 1.0. Subsequently, this embodiment strictly performs piecewise analysis based on the difference between the current pressure and the dynamic critical vaporization pressure. When the difference between the dynamic critical vaporization pressure and the current pressure is greater than the preset anti-divergence pressure difference threshold, and the vapor phase volume fraction is less than the preset vapor phase cutoff threshold, the preset latent heat of vaporization, the corrected evaporation coefficient, one-third of the value of the evaporation interface geometric mapping term, the velocity derivative bias term, and the absolute derivative term are multiplied together and their negative values ​​are taken to analyze the temperature evolution sensitivity of the phase change energy source term under the evaporation condition. The temperature evolution sensitivity is divided by the mixed phase volumetric heat capacity to perform order reduction processing and generate the energy damping factor under the evaporation condition.

[0050] When the difference between the current pressure and the dynamic critical vaporization pressure is greater than the preset anti-divergence pressure difference threshold, and the vapor phase volume fraction is greater than zero, the preset latent heat of vaporization, the preset basic condensation coefficient, one-third of the value of the condensation interface geometric mapping term, the velocity derivative bias term, and the absolute derivative term are multiplied together and their negatives are taken to analyze the temperature evolution sensitivity of the phase change energy source term under condensation conditions. The temperature evolution sensitivity is divided by the mixed phase volumetric heat capacity to perform order reduction processing and generate the energy damping factor under condensation conditions.

[0051] When any of the above operating conditions are not met, the energy damping factor is converged and assigned a value of zero. The specific calculation formula is as follows: In the formula, It is the energy damping factor. For temperature evolution sensitivity, To preset the latent heat of vaporization, To correct the evaporation coefficient, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. This refers to the preset pure vapor phase density within the preset pure phase density. For unit feature scale, This represents the liquid phase volume fraction within each phase volume fraction. This refers to the preset pure liquid phase density within the preset pure phase density. To preset the specific heat capacity of pure liquid phase at constant pressure, To preset the specific heat capacity of the pure vapor phase at constant pressure, The phase change driving pressure difference This is the static saturated vapor pressure. For the fluid element temperature, These are the second and third preset constants in the preset physical property constants. It is the dynamic critical vaporization pressure. Due to current pressure, To preset the vapor phase cutoff threshold, To preset the basic condensation coefficient, The preset anti-divergence pressure difference threshold is set.

[0052] The temperature evolution sensitivity characterizes the slope of the response of the heat released or absorbed during the instantaneous phase change of the flow field to future temperature changes. A larger absolute value indicates a more severe physical response of the flow field to temperature distortion; a smaller absolute value indicates lower sensitivity of heat transfer to temperature fluctuations. This sensitivity is positively dominated by the corrected evaporation coefficient or the preset basic condensation coefficient, as well as the interface geometry mapping term, determining the basic scale involved in the phase change. The velocity derivative bias term, constructed by the phase change-driven pressure difference and the preset pure liquid phase density inversely, reflects that the smaller the pressure difference, the more sensitive the bubble interface motion is to temperature. The absolute derivative term obtained by differentiating the static saturated vapor pressure with respect to the preset second and third constants accurately captures the nonlinear characteristics of vapor pressure rising with temperature. Combined with the preset latent heat of vaporization, this piecewise calculation fits the real fluid mass transfer boundary, objectively restoring the true physical evolution of the phase change heat gradient.

[0053] Furthermore, the energy damping factor characterizes the standardized energy bias correction step size of the underlying computational matrix after eliminating the interference of differences in the thermal inertia of the medium. A larger absolute value means a stronger numerical correction is required for the calculation; a smaller value indicates a smoother current state. This factor is positively correlated with the temperature evolution sensitivity (the numerator) and is significantly suppressed by the mixed-phase volumetric heat capacity (the denominator). Introducing division for order reduction essentially imbues the pure mathematical partial derivative with a real physical thermal inertia constraint, forcing the computational step size to compromise with the local real heat storage capacity, acting as a physical damper, and completely suppressing the false temperature spikes induced by the low heat capacity region of the vapor-rich phase, ensuring absolute convergence of the full-domain calculation to the real operating conditions.

[0054] Table 1 is an example table of comprehensive data on the energy-phase transition strong coupling and damping order reduction evolution of the variable-temperature cavitation flow field at 25 spatial nodes.

[0055] Table 1: Comprehensive Data Table of Energy-Phase Transition Strong Coupling and Damped Order Reduction Evolution Table 1 shows the accurate quantification capability of this method for latent heat throughput and local temperature response in a variable-temperature cavitation flow field. The phase change mass transport rates at nodes 1 to 6 are in a slightly negative range, corresponding to the phase change energy source term. Within a certain range, the energy is released gradually, and the energy damping factor is clamped at a low level, indicating that the flow field is in an initial stable state with localized minor condensation and lack of violent thermodynamic driving force. In nodes 7 to 13, the phase change mass transfer rate reaches the evaporation peak of 45.839 at node 10, causing the phase change energy source term to drop sharply to the endothermic limit of -103.459, resulting in a deep nonlinear step drop in the fluid unit temperature. At this time, the energy damping factor automatically triggers a negative extremum, forcibly applying the thermal inertia constraint of the real fluid to suppress the false temperature oscillations that are prone to occur in the low heat capacity vapor phase region. In nodes 14 to 25, the phase change mass transfer rate reaches the negative collapse peak of -28.913 at node 16. The phase change energy source term released by the collapse of the cavitation bubble drives the fluid unit temperature to rise again. As the cavitation bubble completely disappears, the phase change energy source term and the energy damping factor converge rapidly to zero in sync. This scheme effectively eliminates the non-physical numerical interference generated by conventional isothermal calculations and accurately restores the strongly coupled spatial mapping law of mass transfer and thermodynamic evolution inside the real complex flow field.

[0056] Step 6: Inject the local effective turbulent viscosity and phase change mass transport rate into the momentum conservation equation and cavitation transport equation, respectively. Combine the velocity vector to establish the spatial convection relationship. Inject the phase change energy source term into the energy conservation equation and use the energy damping factor to implement diagonal constraints. Solve the equations iteratively and output the flow field distribution data.

[0057] In this embodiment, the precise physical closure of the governing equations is the underlying foundation for achieving high-fidelity numerical simulation. This embodiment injects the locally effective turbulent viscosity into the viscous diffusion term of the momentum conservation equation, replacing the traditionally constantized macroscopic kinematic viscosity. Due to the extremely strong gas-liquid slip and viscous bed reconstruction in the cavitation region, injecting locally effective turbulent viscosity can accurately strip away non-physical excessive interfacial resistance. Simultaneously, the phase change mass transport rate is injected as a core source term into the cavitation transport equation, further defining the true flow rate of the gas-liquid two-phase mass conversion per unit time, ensuring the dynamic mass conservation of the multiphase medium's volume proportion in the spatial grid.

[0058] In the spatial evolution dimension, due to the accompanying high-speed shearing and transient migration of fluids caused by temperature-variable cavitation, it is necessary to use velocity vectors to establish the spatial convection relationships of various conserved variables between adjacent grids. First, the velocity vector at the interface between adjacent grids is extracted, and a dot product operation is performed between it and the outward normal area vector of the interface to obtain the absolute mass flux, the sign of which defines the true windward direction of the fluid. Second, the upstream and downstream topological relationships are established based on the sign of the absolute mass flux: when the flux is positive, the inner grid node of the interface is taken as the upstream reference; when it is negative, the outer adjacent grid node is taken as the upstream reference. Subsequently, a second-order upwind discretization scheme is adopted, with the value of the conserved variables at the upstream reference as the origin, and combined with the directional distribution of the velocity vector and the grid spatial gradient tensor, Taylor series expansion is used to interpolate and extrapolate the conserved variables to the interface. Finally, the obtained interface conservation variable values ​​are multiplied by the absolute mass flux, and transformed into the surface integral of the net convection flux across the grid boundary according to the Gaussian divergence theorem. This is then directly substituted into the divergence integral terms of each conservation equation to ensure the physical continuity of mass migration. Specifically, each conservation equation covers the cavitation transport equation, momentum conservation equation, energy conservation equation, and turbulence equation, with the vapor volume fraction, velocity vector component, fluid element temperature, and turbulent kinetic energy and turbulent dissipation rate as the corresponding conservation variables.

[0059] To address the cross-dimensional thermodynamic coupling characteristics, this embodiment directly injects the phase change energy source term into the temperature gradient calculation core of the energy conservation equation. Since variable-temperature cavitation phase change is accompanied by intense latent heat release, simple fluid dynamics equations cannot describe the thermodynamic distortion of the medium. Directly injecting this energy source term establishes a mathematical bridge between mechanical work and thermodynamic heat transfer, forcing the fluid temperature field to respond in real-time to the local temperature drops or sudden rises caused by the heat absorption and release during phase change.

[0060] In the discrete solution of partial differential equations, drastically fluctuating phase transition energy source terms can easily disrupt the diagonal dominance of the underlying algebraic equation matrix, inducing divergence in numerical computation. Therefore, this embodiment utilizes the energy damping factor to implement strict matrix diagonal constraints. This energy damping factor is directly superimposed as a bias increment onto the main diagonal elements of the Jacobian matrix, forcibly increasing the absolute value weight of the main diagonal of the underlying iterative equations. This effectively suppresses non-physical oscillations within the implicit integration step size from the underlying algebraic operation structure, greatly improving the convergence stability of the algebraic equation solution under harsh conditions.

[0061] Finally, given the highly nonlinear spatiotemporal dependence between momentum, mass, and energy, this embodiment performs a joint iterative solution to the above conservation equations. By repeatedly iterating within a unified time step, the flow field distribution data is calculated and output after the residuals satisfy the convergence criterion.

[0062] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0063] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0064] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0065] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling, characterized in that, The specific steps include: Obtain the fluid element temperature, current pressure, mixed phase density, velocity vector, turbulent kinetic energy, turbulent dissipation rate, volume fraction of each phase, and element characteristic scale of each fluid element within the preset grid at the current moment; The static saturated vapor pressure is established based on the thermal state mapping between the fluid unit temperature and the preset physical property constants. The turbulent pressure pulsation term is constructed by cross-modulation of the mixed phase density and turbulent kinetic energy. Dynamic compensation is performed on the static saturated vapor pressure to generate the dynamic critical vaporization pressure. A density decay function is constructed based on the evolution characteristics of mixed phase density and preset pure phase density. A reference term for kinematic turbulent viscosity is established based on turbulent kinetic energy and turbulent dissipation rate. The density decay function is used to perform a product substitution operation on the reference term for kinematic turbulent viscosity to generate local effective turbulent viscosity. The phase change driving pressure difference is established based on the deviation between the current pressure and the dynamic critical vaporization pressure. The mass transfer amplification of the preset basic phase change coefficient is performed using turbulent kinetic energy. The phase state logic branch determination is performed in combination with the volume fraction of each phase and the characteristic scale of the unit, and the phase change mass transfer rate is output. Based on the preset latent heat of vaporization, the phase change mass transfer rate is converted into a phase change energy source term. The volumetric heat capacity of the mixed phase is established by combining the volume fraction of each phase. The temperature evolution sensitivity of the phase change energy source term is analyzed in segments. The mixed phase volumetric heat capacity is used to perform a reduction process on the temperature evolution sensitivity to generate an energy damping factor. The local effective turbulent viscosity and phase change mass transport rate are respectively injected into the momentum conservation equation and the cavitation transport equation. The spatial convection relationship is established by combining the velocity vector. The phase change energy source term is injected into the energy conservation equation, and the diagonal constraint is implemented by using the energy damping factor. The equations are solved iteratively and the flow field distribution data is output.

2. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 1, characterized in that: The volume fraction of each phase consists of the volume fraction of the liquid phase and the volume fraction of the vapor phase; The steps for constructing the preset mesh and obtaining the unit feature scale include: Based on the differences in physical flow regime within the target watershed, the watershed space is divided into a high-shear zone near the wall, a phase transition evolution zone, and a core mainstream zone; In the high shear zone near the wall, multiple layers of wall-attached mesh are arranged and the normal thickness of the first layer of wall-attached mesh is defined. In the phase transition evolution zone, a dense mesh with gradually decreasing size is arranged. In the core mainstream zone, a large-scale space-filling mesh is arranged, which together constitute the preset mesh. Extract the mesh volume of each fluid element in the preset mesh, and perform three-dimensional spatial dimensionality reduction feature extraction on the mesh volume to generate the element feature scale corresponding to each fluid element.

3. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 1, characterized in that: The preset physical property constants consist of a preset first constant, a preset second constant, a preset third constant, and a preset reference pressure; The specific steps for generating the dynamic critical vaporization pressure are as follows: The temperature of the fluid element is summed with a preset third constant to construct a temperature offset denominator. The preset second constant is divided by the temperature offset denominator to obtain a nonlinear curvature offset term. The difference between the preset first constant and the nonlinear curvature offset term is calculated to construct a thermal state index term. Using the numerical value of 10 as the base and the thermal state index as the exponent, an exponential mapping operation is performed to obtain a reference mapping value. The reference mapping value, the preset reference pressure, and the numerical value 133.322 are then multiplied together to obtain the static saturated vapor pressure. The preset fluctuation coefficient, mixed phase density, and turbulent kinetic energy are multiplied together to construct the turbulent pressure fluctuation term; The turbulent pressure fluctuation term is used as a dynamic compensation increment and superimposed with the static saturated vapor pressure to generate the dynamic critical vaporization pressure. The specific calculation formula is as follows: In the formula, This is the corrected dynamic critical vaporization pressure. This is the static saturated vapor pressure. To preset the reference pressure, These are the first preset constant, the second preset constant, and the third preset constant, respectively. For the fluid element temperature, To preset the pulsation coefficient, For the mixed phase density, It is turbulent kinetic energy.

4. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 1, characterized in that: The preset pure phase density is composed of the preset pure vapor phase density and the preset pure liquid phase density; The specific steps for generating the local effective turbulent viscosity are as follows: The difference between the mixed phase density and the preset pure vapor phase density is calculated and calibrated as the first difference, and the difference between the preset pure liquid phase density and the preset pure vapor phase density is calculated and calibrated as the second difference. The first difference is divided by the second difference to obtain the phase evolution benchmark term. Using the phase evolution benchmark term as the base and the preset decay exponent as the exponent, a power function mapping operation is performed to construct a density decay function; Divide the square of the turbulent kinetic energy by the turbulent dissipation rate to obtain the turbulent micro-characteristic ratio. Perform a multiplication operation on the turbulent micro-characteristic ratio and the preset turbulence constant to establish the kinematic turbulent viscosity reference term. The density decay function is multiplied and substituted with the kinematic turbulent viscosity reference term to generate the local effective turbulent viscosity. The specific calculation formula is as follows: In the formula, For the local effective turbulent viscosity, For the mixed phase density, To preset the pure vapor phase density, To preset the pure liquid phase density, To preset the attenuation index, As a preset turbulence constant, For turbulent kinetic energy, This represents the turbulent dissipation rate.

5. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 4, characterized in that: The preset basic phase change coefficient is composed of the preset basic evaporation coefficient and the preset basic condensation coefficient; The specific steps for establishing the phase change driving pressure difference based on the deviation between the current pressure and the dynamic critical vaporization pressure, and using turbulent kinetic energy to perform mass transfer amplification on the preset basic phase change coefficient, include: Extract the preset pure liquid phase density from the preset pure phase density; The absolute difference between the dynamic critical vaporization pressure and the current pressure is calculated to establish the phase change driving pressure difference. The phase change driving pressure difference is divided by the preset pure liquid phase density and multiplied by two-thirds of the value, and then the square root operation is performed to generate the bubble dynamics reference term. The sum of the turbulent kinetic energy and the preset anti-divergence reference kinetic energy is used as the denominator, and the turbulent kinetic energy is used as the numerator to perform a ratio calculation. The result of the ratio calculation is multiplied by the preset mass transfer amplification factor and then superimposed with a value of one to construct the mass transfer amplification factor. The preset basic evaporation coefficient and the mass transfer amplification factor are multiplied and modulated to generate the corrected evaporation coefficient. The specific calculation formula is as follows: In the formula, To correct the evaporation coefficient, To preset the basic evaporation coefficient, To preset the mass transfer amplification factor, For turbulent kinetic energy, The preset reference kinetic energy is used to prevent divergence.

6. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 5, characterized in that: The volume fraction of each phase includes the volume fraction of the vapor phase; The specific steps for calculating the phase change mass transfer rate are as follows: The vapor phase volume fraction is extracted from the volume fraction of each phase, and the preset pure vapor phase density is extracted from the preset pure phase density; Based on the magnitudes of the current pressure and the dynamic critical vaporization pressure, a piecewise solution is performed: When the current pressure is lower than the dynamic critical vaporization pressure and the vapor phase volume fraction is less than the preset vapor phase cutoff threshold, the difference between value one and the vapor phase volume fraction is calculated. This difference is then multiplied by the preset condensation nucleus volume fraction, the preset pure vapor phase density, and value three. The resulting product is then divided by the unit characteristic scale to obtain the evaporation interface geometric mapping term. The corrected evaporation coefficient, the evaporation interface geometric mapping term, and the bubble dynamics reference term are then multiplied to calculate the positive phase change mass transfer rate. When the current pressure is higher than the dynamic critical vaporization pressure and the vapor phase volume fraction is greater than zero, the vapor phase volume fraction, the preset pure vapor phase density and the numerical value are multiplied together, and the resulting product is divided by the unit feature scale to obtain the condensation interface geometric mapping term. The preset basic condensation coefficient, the condensation interface geometric mapping term and the bubble dynamics reference term are multiplied together and the negative number is taken to calculate the negative phase change mass transfer rate. When any of the above phase constraint conditions are not met, the fluid unit is determined to be in the phase stability region, and the phase change mass transfer rate is set to zero. The specific calculation formula is as follows: In the formula, For phase change mass transfer rate, To correct the evaporation coefficient, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. To preset the vapor phase cutoff threshold, This refers to the preset pure vapor phase density within the preset pure phase density. For unit feature scale, The phase change driving pressure difference This refers to the preset pure liquid phase density within the preset pure phase density. It is the dynamic critical vaporization pressure. Due to current pressure, This is the preset basic condensation coefficient.

7. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 6, characterized in that: The volume fraction of each phase also includes the volume fraction of the liquid phase; The specific steps for converting the phase change mass transfer rate into a phase change energy source term based on a preset latent heat of vaporization, and for determining the mixed phase volumetric heat capacity in conjunction with the volume fraction of each phase, include: The preset latent heat of vaporization and the phase change mass transfer rate are multiplied, and the negative is taken to convert the phase change mass transfer rate into a phase change energy source term. The specific calculation formula is as follows: In the formula, For phase transition energy source term, To preset the latent heat of vaporization, This refers to the phase change mass transfer rate; The volume fraction of the liquid phase is extracted from the volume fraction of each phase. The liquid phase heat capacity term is obtained by performing a multiplication operation on the preset pure liquid phase constant pressure specific heat capacity, the preset pure liquid phase density, and the liquid phase volume fraction. The vapor phase heat capacity term is obtained by performing a multiplication operation on the preset pure vapor phase constant pressure specific heat capacity, the preset pure vapor phase density, and the vapor phase volume fraction. The liquid phase heat capacity term and the vapor phase heat capacity term are summed to establish the mixed phase volume heat capacity.

8. The numerical calculation method for variable-temperature cavitation considering turbulence correction and energy-phase transition coupling according to claim 7, characterized in that: The preset physical property constants include a preset second constant and a preset third constant; The specific steps for generating the energy damping factor are as follows: Extract the preset second and third constants from the preset physical property constants; The fluid element temperature is summed with the preset third constant. The square of the summation result is used as the denominator of the partial derivative. The product of the preset second constant and the natural logarithm of the value of 10 is calculated and used as the numerator of the partial derivative. The numerator of the partial derivative is divided by the denominator of the partial derivative to obtain the derivative nonlinear coefficient. Then, it is multiplied with the static saturated vapor pressure to construct the absolute derivative term. Calculate the product of the phase change driving pressure difference and the preset pure liquid phase density and two-thirds of the value, and take the square root of the product and find its reciprocal as the velocity derivative bias term. Segmented analysis is performed based on the magnitude of the current pressure and the dynamic critical vaporization pressure: When the difference between the dynamic critical vaporization pressure and the current pressure is greater than the preset anti-divergence pressure difference threshold, and the vapor phase volume fraction is less than the preset vapor phase cutoff threshold, the preset latent heat of vaporization, the corrected evaporation coefficient, one-third of the values ​​of the evaporation interface geometric mapping term, the velocity derivative bias term, and the absolute derivative term are multiplied together and the negative is taken to analyze the temperature evolution sensitivity of the phase change energy source term under the evaporation condition. The temperature evolution sensitivity is divided by the mixed phase volume heat capacity to perform order reduction processing and generate the energy damping factor under the evaporation condition. When the difference between the current pressure and the dynamic critical vaporization pressure is greater than the preset anti-divergence pressure difference threshold, and the vapor phase volume fraction is greater than zero, the preset latent heat of vaporization, the preset basic condensation coefficient, one-third of the value of the condensation interface geometric mapping term, the velocity derivative bias term, and the absolute derivative term are multiplied together and the negative is taken to analyze the temperature evolution sensitivity of the phase change energy source term under condensation conditions. The temperature evolution sensitivity is divided by the mixed phase volume heat capacity to perform order reduction processing and generate the energy damping factor under condensation conditions. When any of the above operating conditions are not met, the energy damping factor will converge and be assigned a value of zero. The specific calculation formula is as follows: In the formula, It is the energy damping factor. For temperature evolution sensitivity, To preset the latent heat of vaporization, To correct the evaporation coefficient, To preset the volume fraction of condensation nuclei, This represents the vapor phase volume fraction within the total volume fraction of each phase. This refers to the preset pure vapor phase density within the preset pure phase density. For unit feature scale, This represents the liquid phase volume fraction within each phase volume fraction. This refers to the preset pure liquid phase density within the preset pure phase density. To preset the specific heat capacity of pure liquid phase at constant pressure, To preset the specific heat capacity of the pure vapor phase at constant pressure, The phase change driving pressure difference This is the static saturated vapor pressure. For the fluid element temperature, These are the second and third preset constants in the preset physical property constants. It is the dynamic critical vaporization pressure. Due to current pressure, To preset the vapor phase cutoff threshold, To preset the basic condensation coefficient, The preset anti-divergence pressure difference threshold is set.

Images