Numerical Simulation Method and System for Cavitation Turbulence Based on Gas-Liquid Slip Velocity Correction

Abstract

This invention provides a numerical simulation method and system for cavitation turbulence based on gas-liquid slip velocity correction, belonging to the field of cavitation numerical simulation technology. It obtains the gas phase volume fraction distribution and liquid phase velocity field distribution in a fluid pipe, establishes the density gradient field of the gas-liquid interface region, constructs a dynamic characteristic function of the gas-liquid interface including a coupled expression of interface curvature and local vorticity terms, and uses an adaptive gas-liquid slip velocity correction model to correct the relative gas-liquid velocities, solving for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field. This significantly improves the accuracy and refinement of the gas-liquid interface dynamics description, enhances the adaptability of the slip velocity correction model to the nonlinear characteristics of the interface, and helps to more accurately simulate cavitation turbulence behavior, thereby improving the accuracy and reliability of cavitation turbulence numerical simulation and providing more scientific and effective technical support for cavitation prediction and control in complex engineering fields.

Description

Technical Field

[0001] This invention relates to the field of cavitation numerical simulation technology, and more specifically, to a cavitation turbulence numerical simulation method and system based on gas-liquid slip velocity correction. Background Technology

[0002] Cavitation is an important research subject in fluid mechanics, particularly in high-speed flow, liquid transport, and marine engineering. Cavitation typically occurs when the local pressure of a liquid drops below its saturated vapor pressure, causing the liquid to vaporize and form bubbles, which collapse as the pressure recovers. This process is accompanied by complex phase transitions, turbulence, and gas-liquid interface dynamics, resulting in highly nonlinear and multi-scale characteristics of cavitation. In recent years, with the development of numerical simulation techniques, numerical simulation of cavitation turbulence has become an important tool for studying cavitation mechanisms. Traditional cavitation models are mostly based on the Reynolds-averaged method (RANS) or the large eddy simulation (LES), combined with the volume fraction method (VOF) or the Euler-Lagrange method to describe the dynamic behavior of gas-liquid two-phase flow. However, these methods often rely on empirical models when dealing with the dynamic characteristics of the gas-liquid interface and slip velocity corrections, failing to fully reflect the interfacial nonlinearities and interphase coupling effects during the cavitation process.

[0003] While existing technologies have made significant progress in numerical simulation of cavitation turbulence, several shortcomings remain. First, current gas-liquid interface dynamics modeling methods provide a coarse description of the coupling between the interface's geometric properties and the fluid's local vorticity, making it difficult to accurately capture the influence of interface curvature changes on cavitation behavior. Second, traditional gas-liquid slip velocity models typically employ fixed empirical formulas, failing to fully consider the local dynamic characteristics of the gas-liquid interface and turbulence characteristic parameters, resulting in room for improvement in the accuracy and robustness of numerical simulation results. Furthermore, in solving the turbulent strain rate transport equations in the cavitation flow field, existing methods simplify the treatment of relative gas-liquid velocities, neglecting the feedback effect of nonlinear interface coupling on turbulence characteristics, thus limiting the accuracy of turbulence characteristic parameter predictions. These problems not only affect the reliability of numerical simulations of cavitation turbulence but also pose a serious challenge to the prediction and control of cavitation behavior in complex engineering applications. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention is proposed. This invention provides a numerical simulation method and system for cavitation turbulence based on gas-liquid slip velocity correction.

[0005] According to one aspect of the present invention, a numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction is provided, comprising:

[0006] The gas phase volume fraction distribution and liquid phase velocity field distribution in the fluid pipeline are obtained, and the density gradient field of the gas-liquid two-phase interface region is established based on the gas phase volume fraction distribution.

[0007] Using the density gradient field and liquid phase velocity field distribution, a dynamic characteristic function of the gas-liquid interface is constructed. The dynamic characteristic function includes a coupled expression of the interface curvature term and the local vorticity term.

[0008] Based on the numerical distribution of the dynamic characteristic function, an adaptive gas-liquid slip velocity correction model is used to correct the relative velocity of gas and liquid.

[0009] Substitute the modified gas-liquid relative velocity into the turbulent strain rate transport equation to solve for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field.

[0010] Furthermore, a density gradient field is established in the gas-liquid interface region, and a gas phase volume fraction distribution function is established by detecting and identifying bubble characteristics using acoustic emission, thereby calculating the gas-liquid two-phase mixing density.

[0011] The spatial rate of change of the mixed density is calculated using the central difference scheme. The interface region is identified by combining the magnitude and direction vector of the density gradient with the eigenvalues ​​of the second derivative tensor. Finally, a density gradient field characterizing the interface intensity and direction is constructed.

[0012] Furthermore, the phase distribution of the gas-liquid two-phase mixture density is determined based on the gas phase volume fraction. The two-phase region is weighted and the gas phase density is corrected by pressure. The mixture density is obtained through the weighted contribution value of the gas and liquid phases.

[0013] Furthermore, constructing the dynamic characteristic function of the gas-liquid interface includes:

[0014] The density Laplacian operator is obtained by performing divergence calculations on the density gradient field, and then the average curvature expression is constructed by combining the direction information of the density gradient.

[0015] The vorticity vector is obtained by calculating the curl of the velocity field;

[0016] The mean curvature and vorticity vector are projected onto the density gradient direction to obtain the interface normal curvature and normal vorticity components.

[0017] Dynamic characteristic functions of the gas-liquid interface are constructed based on normal curvature and normal vorticity components.

[0018] Furthermore, by calculating the Laplacian operator and normal vector of the density gradient field, the principal curvature is obtained on the interface tangent plane, and finally the mean curvature expression with mesh independence is obtained;

[0019] The vorticity vector is obtained by calculating the velocity gradient in each direction of the liquid phase velocity field, then obtaining the velocity field curl to get the vorticity vector, and finally using the time averaging method to obtain a stable characteristic vorticity.

[0020] Furthermore, the steps for correcting the relative velocity of gas and liquid include: coupling the dynamic characteristic function with the pressure gradient and vorticity field respectively, and combining the direction modulation factor to calculate the correction coefficients of the rising velocity and the lateral drift velocity of the bubble motion respectively, and using the correction coefficients to correct the relative velocity of gas and liquid.

[0021] Furthermore, the turbulent strain rate transport equation is shown below:

[0022] ;

[0023] in, For turbulent strain rate tensor, For time, For the average velocity component, For fluid density, For pressure, Kinematic viscosity, For the generation of relative motion turbulence, For interfacial stress, To add turbulent dissipation terms, This represents the partial derivative operator with respect to the spatial coordinate k in the direction of k. Let i represent the partial derivative operator with respect to the direction of spatial coordinate i. This represents the partial derivative operator with respect to the spatial coordinate j in the direction of j. It is the Laplace operator.

[0024] According to another aspect of the present invention, a numerical simulation system for cavitation turbulence based on gas-liquid slip velocity correction is provided, comprising:

[0025] The data acquisition and differentiation module is used to acquire the gas phase volume fraction distribution and liquid phase velocity field distribution in the fluid pipeline, and to establish the density gradient field of the gas-liquid two-phase interface region based on the gas phase volume fraction distribution.

[0026] The dynamic characteristic function construction module is used to construct a dynamic characteristic function of the gas-liquid interface using the density gradient field and liquid phase velocity field distribution. The dynamic characteristic function includes a coupled expression of the interface curvature term and the local vorticity term.

[0027] The correction module is used to correct the relative velocity of gas and liquid based on the numerical distribution of the dynamic characteristic function using an adaptive gas-liquid slip velocity correction model.

[0028] The parameter solving module is used to substitute the modified gas-liquid relative velocity into the turbulent strain rate transport equation to solve for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field.

[0029] Compared with existing technologies, the numerical simulation method and system for cavitation turbulence based on gas-liquid slip velocity correction provided by this invention obtains the gas phase volume fraction distribution and liquid phase velocity field distribution in a fluid pipe, establishes the density gradient field of the gas-liquid interface region based on the gas phase volume fraction distribution, further constructs a dynamic characteristic function of the gas-liquid interface including a coupled expression of interface curvature term and local vorticity term, and uses an adaptive gas-liquid slip velocity correction model to correct the relative velocity of gas and liquid, thereby solving for the turbulence characteristic parameters and pressure field distribution in the cavitation flow field. This significantly improves the accuracy and refinement of the gas-liquid interface dynamics description, enhances the adaptability of the slip velocity correction model to the nonlinear characteristics of the interface, and helps to more accurately simulate cavitation turbulence behavior, thus improving the accuracy and reliability of cavitation turbulence numerical simulation and providing more scientific and effective technical support for cavitation prediction and control in complex engineering fields. Attached Figure Description

[0030] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings:

[0031] Figure 1 This is a system block diagram of a numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to an embodiment of the present invention.

[0032] Figure 2 This is a comparison diagram of the turbulent viscosity distribution in the gas-liquid interface region between the traditional method and the present invention in the numerical simulation method of cavitation turbulence based on gas-liquid slip velocity correction according to an embodiment of the present invention. Detailed Implementation

[0033] Hereinafter, exemplary embodiments according to the present invention will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are merely some embodiments of the present invention, and not all embodiments of the present invention. It should be understood that the present invention is not limited to the exemplary embodiments described herein.

[0034] Figure 1 This is a system block diagram of a numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to an embodiment of the present invention. Figure 1 As shown, the numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction includes:

[0035] S1: Obtain the gas phase volume fraction distribution and liquid phase velocity field distribution in the fluid pipeline, and establish the density gradient field of the gas-liquid two-phase interface region based on the gas phase volume fraction distribution;

[0036] The process involves acquiring the gas phase volume fraction distribution and liquid phase velocity field distribution in a fluid pipeline, and establishing a density gradient field for the gas-liquid interface region based on the gas phase volume fraction distribution. Specifically, this includes: installing multiple pressure and velocity sensors at the inlet, middle, and outlet sections of the fluid pipeline; the pressure sensors collecting pressure data and the velocity sensors collecting velocity data; using acoustic emission detection technology to identify the acoustic signal characteristics of bubble bursting and merging based on the pressure data, and establishing a gas phase volume fraction distribution function; using a Doppler current meter to measure the instantaneous velocity of the liquid phase, acquiring the liquid phase velocity field distribution, which includes radial and axial velocity components; calculating the gas-liquid two-phase mixing density using the gas phase volume fraction distribution function, and calculating the rate of change of the mixing density in space using a central difference scheme, thus establishing a density gradient field for the gas-liquid interface region. This density gradient field characterizes the position and intensity distribution of the gas-liquid interface.

[0037] A gas phase volume fraction distribution function is established by identifying the acoustic signal characteristics of bubble bursting and merging using acoustic emission detection technology. Specifically, an acoustic emission sensor array is deployed at key locations in the fluid pipeline to collect acoustic signals generated during bubble movement. Wavelet transform processing is performed on the acoustic signals to obtain a time-frequency characteristic spectrum. When a bubble bursts, the acoustic signal exhibits high-frequency short-pulse characteristics; when bubbles merge, the acoustic signal exhibits low-frequency long-duration characteristics. Using the energy distribution characteristics in the time-frequency characteristic spectrum, a correspondence between acoustic signal intensity and bubble size is established. The acoustic signal energy during bubble bursting is proportional to the bubble size before bursting, and the acoustic signal energy during bubble merging is proportional to the total volume of the merging bubbles. Based on this correspondence, an acoustic signal triangulation algorithm is used to determine the spatial location of the bubbles. Combined with bubble size information, a volume fraction distribution function characterizing the spatial distribution of the gas phase is established, as shown in the following equation:

[0038] ;

[0039] in, Let be the gas phase volume fraction distribution function. Let be the acoustic signal energy value of the i-th bubble. Let be the volume of the i-th bubble. To calculate the total volume of the region, The acoustic signal attenuation coefficient, The distance from the bubble to the calculation point. Let be the distance from the i-th bubble to the j-th sensor. The number of bubbles identified, This represents the number of sensors.

[0040] Furthermore, the gas-liquid two-phase mixing density is calculated using the gas phase volume fraction distribution function. Specifically, the pure liquid phase density and pure gas phase density under standard conditions are obtained as reference parameters. Based on the gas phase volume fraction distribution function, the gas phase proportion at each calculation point in the flow field is determined. When the gas phase volume fraction is greater than zero and less than one, it indicates that the point is in the gas-liquid two-phase region. The mixing density is obtained using a gas-liquid two-phase weighted calculation method, where the weighting coefficient of the gas phase density is the value of the gas phase volume fraction distribution function, and the weighting coefficient of the liquid phase density is one minus the value of the gas phase volume fraction distribution function. If the gas phase volume fraction at the point is equal to zero, the mixing density is directly... The liquid phase density value is used. If the gas phase volume fraction at that point is equal to one, the gas phase density value is directly used. For the gas-liquid two-phase region, the influence of local pressure on the gas phase density is considered. A pressure correction coefficient is calculated based on the ideal gas law. The pressure correction coefficient is equal to the ratio of the local pressure to the standard atmospheric pressure. The pressure correction coefficient is multiplied by the gas phase density to obtain the corrected gas phase density. The corrected gas phase density is multiplied by the gas phase weighting coefficient to obtain the gas phase contribution value. The corrected liquid phase density is multiplied by the liquid phase weighting coefficient to obtain the liquid phase contribution value. The gas phase contribution value and the liquid phase contribution value are added together to obtain the gas-liquid two-phase mixing density.

[0041] Furthermore, a density gradient field for the gas-liquid two-phase interface region is established by calculating the rate of change of the mixing density in space using a central difference scheme. Specifically, the computational domain is divided into a structured grid, and the mixing density value of the gas-liquid two-phase mixture is obtained at each grid node. The density gradient is calculated based on the mixing density values ​​of adjacent grid nodes. The density gradient in the x-direction is obtained by dividing the difference in mixing density between the nodes before and after the node by the distance between the two nodes. The density gradients in the y-direction and z-direction are calculated using the same method. When the computational node is located at the boundary of the computational domain, the density gradient of the node is calculated using a forward difference or backward difference method. The magnitude and direction vector of the density gradient are calculated based on the density gradients in the x, y, and z directions, and the direction vector points in the direction of increasing mixing density. An adaptive curvature evaluation method is used to identify the interface region. The interface curvature is determined by calculating the eigenvalues ​​of the second derivative tensor of the local density field. When the curvature value indicates a significant change in the local density distribution, the region is marked as an interface region. The density gradient magnitude and direction vector are combined to form a density gradient field, where the density gradient magnitude characterizes the interface strength, and the direction vector characterizes the spatial orientation of the interface.

[0042] S2: Using the density gradient field and liquid phase velocity field distribution, a dynamic characteristic function of the gas-liquid interface is constructed. The dynamic characteristic function includes a coupled expression of the interface curvature term and the local vorticity term.

[0043] Based on the calculation of interface curvature characteristics using the density gradient field, the density Laplacian operator is first obtained through divergence calculation of the density gradient field, and then the average curvature expression is constructed by combining the direction information of the density gradient. Next, the local vorticity characteristics are calculated based on the liquid phase velocity field distribution, and the vorticity vector is obtained by curl calculation of the velocity field. Subsequently, the average curvature and vorticity vector are projected onto the density gradient direction to obtain the interface normal curvature and normal vorticity components. Then, a dynamic characteristic function of the gas-liquid interface is constructed based on the normal curvature and normal vorticity components. Finally, an adaptive weighting method is used to adjust the contribution ratio of the curvature term and the vorticity term.

[0044] Specifically, the density Laplacian operator is obtained through divergence calculation of the density gradient field, and the average curvature expression is constructed by combining the density gradient direction. The process involves: calculating the second-order partial derivatives in each direction using the central difference method based on the components of the density gradient field in three-dimensional space; summing the second-order partial derivatives in the x, y, and z directions to obtain the density Laplacian operator; normalizing the density Laplacian operator using the magnitude of the density gradient field to obtain the density change rate per unit density gradient direction; calculating the interface normal vector based on the density gradient field, and obtaining the unit normal vector by dividing the direction vector of the density gradient field by its magnitude; using the principal curvature calculation method in differential geometry, constructing an interface tangent plane based on the unit normal vector, calculating the directional derivative of the density field on the tangent plane to obtain the principal curvatures in two orthogonal directions on the interface; averaging the two principal curvatures to obtain the average curvature; constructing an interface shape factor based on the average curvature, and matching the average curvature with the local mesh scale to ensure mesh independence in curvature calculation. More specifically, the average curvature expression is shown below: ;

[0045] in, For the mean curvature, and For two principal curvatures, For density field, For density Laplacian operator, This represents the density gradient magnitude. It is a unit normal vector and Let be the second directional derivative of the density field in the direction of the normal vector. For local grid size, For the local radius of curvature, For the nabla operator.

[0046] Furthermore, based on the velocity components of the liquid phase velocity field in three-dimensional space, the velocity gradient in each direction is calculated using the central difference method, where the partial derivatives of the velocity in the x-direction with respect to y and z, the partial derivatives of the velocity in the y-direction with respect to x and z, and the partial derivatives of the velocity in the z-direction with respect to x and y are calculated. When the calculation node is located at the boundary, the velocity gradient is calculated using the one-sided difference method. The curl of the velocity field is calculated using the velocity gradient, where the x-component of the vorticity vector is equal to the partial derivative of the velocity in the z-direction with respect to y minus the partial derivative of the velocity in the y-direction with respect to z, the y-component is equal to the partial derivative of the velocity in the x-direction with respect to z minus the partial derivative of the velocity in the z-direction with respect to x, and the z-component is equal to the partial derivative of the velocity in the y-direction with respect to x minus the partial derivative of the velocity in the x-direction with respect to y. The vorticity intensity is calculated based on the vorticity vector, and the vorticity intensity is equal to the magnitude of the vorticity vector. The time averaging method is used to eliminate high-frequency fluctuations in the vorticity calculation. When the flow field is in a stable state, the average vorticity value within a specified time period is taken as the characteristic vorticity at that location.

[0047] Furthermore, a unit normal vector for the interface is constructed based on the density gradient field. This unit normal vector is obtained by dividing the density gradient vector by its magnitude, and it points in the direction of increasing density. The interface normal curvature is obtained by performing an inner product operation between the average curvature and the unit normal vector. The projection component of the vorticity in the interface normal direction is obtained by performing an inner product operation between the vorticity vector and the unit normal vector. This normal vorticity component characterizes the rotational intensity of the fluid around the normal axis near the interface. When the angle between the vorticity vector and the interface normal vector is close to 90 degrees, it indicates that the vortex motion mainly occurs in the tangential direction of the interface, at which point the normal vorticity component is close to zero. The vorticity vector is decomposed into normal and tangential components using a local coordinate transformation method. Dynamic characteristic parameters of the interface are constructed based on the interface normal curvature and the normal vorticity component to characterize the local deformation and rotational properties of the interface.

[0048] Finally, a dimensionless curvature parameter is constructed based on the absolute value of the interface normal curvature, obtained by multiplying the interface normal curvature by the characteristic length; a dimensionless vorticity parameter is constructed based on the normal vorticity component, obtained by multiplying the normal vorticity component by the characteristic time; the dimensionless curvature parameter and the dimensionless vorticity parameter are combined using an adaptive weighting coefficient; a dimensional correction function is constructed based on the Weber number and Reynolds number to adjust the influence of surface tension and viscous force on the dynamic characteristics of the interface. A larger Weber number indicates that inertial force is dominant, while a larger Reynolds number indicates that viscous dissipation is weak; the dimensionless parameter and the dimensional correction function are combined to construct the gas-liquid interface dynamic characteristic function, as shown in the following equation:

[0049] ;

[0050] in, This is the dynamic characteristic function of the gas-liquid interface. For characteristic length, The interface normal curvature (obtained by the inner product of the mean curvature and the unit normal vector). For characteristic time, It is the normal vorticity component (obtained by the inner product of the vorticity vector and the unit normal vector). Given the Weber number and We = ρU²L / σ (ρ is the density, U is the characteristic velocity, and σ is the surface tension coefficient), The critical Weber number. Given the Reynolds number and Re = UL / ν (ν is the kinematic viscosity). The critical Reynolds number is denoted as .

[0051] S3: Based on the numerical distribution of the dynamic characteristic function, an adaptive gas-liquid slip velocity correction model is used to correct the relative velocity of gas and liquid;

[0052] The bubble's rising velocity and lateral drift velocity are determined using a standard slip velocity model, considering only the balance between buoyancy and drag. Correction coefficients for both rising velocity and lateral drift velocity are constructed based on the values ​​of the dynamic characteristic function. An increase in the dynamic characteristic function value indicates increased interface instability, and both correction coefficients increase accordingly. The corrected rising velocity is obtained by multiplying the corrected rising velocity coefficients by the bubble's rising velocity reference value. These correction coefficients are determined through the coupling relationship between the dynamic characteristic function and the local pressure gradient. The corrected rising velocity correction coefficient reaches its maximum value when the angle between the pressure gradient and the gravity direction is small and the dynamic characteristic function value is large. Similarly, the corrected lateral drift velocity is obtained by multiplying the lateral drift velocity correction coefficient by the lateral drift velocity reference value. This correction coefficient is determined through the coupling relationship between the dynamic characteristic function and the local vorticity field. The lateral drift velocity correction coefficient reaches its maximum value when the local vorticity intensity is large and the dynamic characteristic function value is large. A gas-liquid relative velocity vector is constructed based on the corrected rising velocity and lateral drift velocity. This relative velocity vector is dynamically adjusted according to changes in the interface dynamics and local flow field characteristics.

[0053] Based on the numerical values ​​of the dynamic characteristic function, the ascent velocity correction coefficient and the lateral drift velocity correction coefficient are constructed respectively. Specifically, the ascent velocity correction coefficient is calculated based on the numerical values ​​of the dynamic characteristic function. A pressure coupling term is obtained by performing a dot product operation between the dynamic characteristic function and the local pressure gradient. The pressure coupling term represents the projection intensity of the pressure gradient in the gravity direction. A direction modulation function is constructed using the cosine value of the angle between the local pressure gradient and the gravity direction. The direction modulation function reaches its maximum value when the angle between the pressure gradient and the gravity direction is close to zero degrees. The ascent velocity correction coefficient is obtained by multiplying the pressure coupling term with the direction modulation function and mapping it through an exponential function. Based on the numerical values ​​of the dynamic characteristic function... The lateral drift velocity correction coefficient is calculated by performing a cross product operation between the dynamic characteristic function and the local vorticity field to obtain the vorticity coupling term, which characterizes the intensity of the tangential vorticity at the interface. The vorticity direction function is constructed using the cosine of the angle between the local vorticity vector and the interface normal vector. The vorticity direction function reaches its maximum value when the vorticity direction is perpendicular to the interface normal vector. The vorticity coupling term is multiplied by the vorticity direction function and mapped through a logarithmic function to obtain the lateral drift velocity correction coefficient. When the value of the dynamic characteristic function is large and the local flow field characteristics meet the corresponding conditions, the values ​​of the two correction coefficients increase accordingly, reflecting the modulation effect of the interface dynamic characteristics on the relative motion of gas and liquid, as shown in the following equations.

[0054] ;

[0055] ;

[0056] in, This is a correction factor for the rate of ascent. This is the lateral drift velocity correction factor. For dynamic characteristic functions, For local pressure gradient, For reference pressure, Let be the angle between the pressure gradient and the direction of gravity. This is the local vorticity vector. For reference vorticity, The unit normal vector of the interface. The magnitude of the local vorticity vector. Let be the magnitude of the unit normal vector of the interface.

[0057] S4: Substitute the modified gas-liquid relative velocity into the turbulent strain rate transport equation considering interfacial stress, and solve for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field.

[0058] First, the velocity gradient tensor is obtained by performing gradient calculations on the corrected gas-liquid relative velocities. This tensor contains the velocity shear distribution characteristics near the gas-liquid interface. A turbulent generation term for relative motion is constructed in the turbulent strain rate transport equation. This generation term is obtained by calculating the double dot product of the velocity gradient tensor and the Reynolds stress tensor. When calculating the interface stress term, the surface tension gradient stress is calculated using the interface curvature gradient, and the Marangoni stress is calculated using the interface temperature gradient. The resultant force of these two stresses is substituted into the turbulent strain rate transport equation as the interface stress term. An additional turbulent dissipation term is calculated based on the deformation rate of the gas-liquid interface. This dissipation term is proportional to the square of the interface deformation rate. By solving the corrected turbulent strain rate transport equation, the spatial distribution of turbulent kinetic energy and turbulent dissipation rate is obtained, and then the turbulent viscosity coefficient is calculated. The turbulent strain rate transport equation is shown below:

[0059] ;

[0060] ;

[0061] ;

[0062] ;

[0063] in, For time partial derivative operators, For turbulent strain rate tensor, For time, For the average velocity component, For fluid density, For flow field pressure, For fluid kinematic viscosity, For the generation of relative motion turbulence, For interfacial stress, To add turbulent dissipation terms, This is the corrected component of the relative gas-liquid velocity in the i-direction. This is the corrected component of the relative gas-liquid velocity in the j-direction. and Let Reynolds stress tensor be the stress tensor. The surface tension coefficient, For the interface curvature, For the interface Delta function, The Marangoni stress coefficient is... For temperature, This is the dissipation adjustment coefficient. This represents the amount of interface deformation. This represents the partial derivative operator with respect to the spatial coordinate k in the direction of k. Let i represent the partial derivative operator with respect to the direction of spatial coordinate i. This represents the partial derivative operator with respect to the spatial coordinate j in the direction of j. It is the Laplace operator. , , It represents the position component in the spatial coordinate system.

[0064] It should be noted that the turbulent viscosity coefficient is calculated based on the turbulent kinetic energy and turbulent dissipation rate obtained from the solution. Specifically, the following steps are taken: First, the local turbulent strain rate distribution is calculated using the turbulent strain rate transport equation. The generation and dissipation terms of the turbulent kinetic energy are then calculated based on the strain rate tensor. These generation and dissipation terms are substituted into the turbulent kinetic energy transport equation and the turbulent dissipation rate transport equation. The convection term is discretized using a second-order upwind scheme, and the diffusion term is discretized using a central difference scheme. When the gas-liquid interface deforms drastically, an additional source term caused by interface deformation is added to the turbulent kinetic energy equation. This source term is proportional to the square of the interface deformation rate. When a phase transition occurs, a phase transition term is added to the turbulent dissipation rate equation. An additional source term caused by the change is obtained, which is proportional to the product of the phase transition rate and the local turbulent time scale. The spatial distribution of turbulent kinetic energy and turbulent dissipation rate is obtained by solving the modified turbulent transport equations. The characteristic time scale is calculated based on the turbulent kinetic energy and turbulent dissipation rate, and the product with the turbulent kinetic energy is used as the benchmark value of the turbulent viscosity coefficient. When the gas-liquid interface exists, the interface shape function is constructed by using the ratio of the interface curvature to the local turbulent length scale, and the interface normal vector is introduced to calculate the turbulent anisotropy tensor. The benchmark value of the turbulent viscosity coefficient is corrected by using the interface shape function and the turbulent anisotropy tensor, and finally the turbulent viscosity coefficient considering the interface effect is obtained.

[0065] The turbulent viscosity coefficient is substituted into the momentum transport equation, and the pressure field is solved using the SIMPLE algorithm. The source term of the pressure correction equation is determined by the mass conservation residual. When the local pressure obtained is lower than the saturated vapor pressure, the region is identified as a cavitation region and its pressure is fixed at the saturated vapor pressure. Finally, the obtained pressure field is substituted into the dynamic characteristic function for calculation, thereby completing the coupled solution of the pressure field and the dynamic characteristics of the interface.

[0066] In summary, the numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to embodiments of the present invention is elucidated. It obtains the gas phase volume fraction distribution and liquid phase velocity field distribution in a fluid pipe, establishes the density gradient field of the gas-liquid interface region based on the gas phase volume fraction distribution, further constructs a dynamic characteristic function of the gas-liquid interface including a coupled expression of interface curvature and local vorticity terms, and uses an adaptive gas-liquid slip velocity correction model to correct the relative gas-liquid velocity, thereby solving for the turbulence characteristic parameters and pressure field distribution in the cavitation flow field. This significantly improves the accuracy and refinement of the gas-liquid interface dynamics description, enhances the adaptability of the slip velocity correction model to the nonlinear characteristics of the interface, and helps to more accurately simulate cavitation turbulence behavior, thus improving the accuracy and reliability of cavitation turbulence numerical simulation and providing more scientific and effective technical support for cavitation prediction and control in complex engineering fields.

[0067] According to another aspect of the present invention, a numerical simulation system for cavitation turbulence based on gas-liquid slip velocity correction is provided, comprising:

[0068] The data acquisition and differentiation module is used to acquire the gas phase volume fraction distribution and liquid phase velocity field distribution in the fluid pipeline, and to establish the density gradient field of the gas-liquid two-phase interface region based on the gas phase volume fraction distribution.

[0069] The dynamic characteristic function construction module is used to construct a dynamic characteristic function of the gas-liquid interface using the density gradient field and liquid phase velocity field distribution. The dynamic characteristic function includes a coupled expression of the interface curvature term and the local vorticity term.

[0070] The correction module is used to correct the relative velocity of gas and liquid based on the numerical distribution of the dynamic characteristic function using an adaptive gas-liquid slip velocity correction model.

[0071] The parameter solving module is used to substitute the modified gas-liquid relative velocity into the turbulent strain rate transport equation to solve for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field.

[0072] Here, those skilled in the art will understand that the specific operations of each step in the numerical simulation system of cavitation turbulence based on gas-liquid slip velocity correction have been referenced above. Figure 1 and Figure 2 The numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction has been described in detail, and therefore, its repeated description will be omitted.

[0073] It should be noted that, Figure 2The left-middle figure shows the distribution of turbulent viscosity coefficient calculated by the conventional method, the middle figure shows the distribution of turbulent viscosity coefficient calculated by the method of this invention, and the right figure shows a comparative profile of turbulent viscosity coefficient along the center line. The improved method proposed in this invention significantly optimizes the conventional method by introducing the interface effect into the calculation of turbulent viscosity coefficient. It can be clearly seen from the figures that in the gas-liquid interface region (the central region of the figure), the turbulent viscosity coefficient calculated by the method of this invention is significantly lower than that of the conventional method. This reduction effect is consistent with actual physical phenomena, as the presence of the interface suppresses the local turbulence intensity. A quantitative observation through the comparison of the center line profile shows that in the interface region, the method of this invention reduces the turbulent viscosity coefficient by approximately 35.93% compared to the conventional method. In regions far from the interface, the results of the two methods tend to be consistent. This distribution characteristic accurately reflects the local modulation effect of the interface on turbulent transport, thus enabling better prediction of turbulent characteristics in gas-liquid two-phase flow.

[0074] In summary, the numerical simulation system for cavitation turbulence based on gas-liquid slip velocity correction according to embodiments of the present invention is elucidated. It obtains the gas phase volume fraction distribution and liquid phase velocity field distribution in a fluid pipe, establishes the density gradient field of the gas-liquid interface region based on the gas phase volume fraction distribution, and further constructs a dynamic characteristic function of the gas-liquid interface including a coupled expression of interface curvature and local vorticity terms. An adaptive gas-liquid slip velocity correction model is used to correct the relative gas-liquid velocity, thereby solving for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field. This significantly improves the accuracy and refinement of the gas-liquid interface dynamics description, enhances the adaptability of the slip velocity correction model to the nonlinear characteristics of the interface, and helps to more accurately simulate cavitation turbulence behavior, thus improving the accuracy and reliability of cavitation turbulence numerical simulation and providing more scientific and effective technical support for cavitation prediction and control in complex engineering fields.

Claims

1. A numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction, characterized in that, include: The gas phase volume fraction distribution and liquid phase velocity field distribution in the fluid pipeline are obtained, and the density gradient field of the gas-liquid two-phase interface region is established based on the gas phase volume fraction distribution. Using the density gradient field and liquid phase velocity field distribution, a dynamic characteristic function of the gas-liquid interface is constructed. The dynamic characteristic function includes a coupled expression of the interface curvature term and the local vorticity term. Based on the numerical distribution of the dynamic characteristic function, an adaptive gas-liquid slip velocity correction model is used to correct the relative velocity of gas and liquid. Substitute the modified gas-liquid relative velocity into the turbulent strain rate transport equation to solve for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field. Constructing the dynamic characteristic function of the gas-liquid interface includes: The density Laplacian operator is obtained by performing divergence calculations on the density gradient field, and then the average curvature expression is constructed by combining the direction information of the density gradient. The vorticity vector is obtained by calculating the curl of the velocity field; The mean curvature and vorticity vector are projected onto the density gradient direction to obtain the interface normal curvature and normal vorticity components. Dynamic characteristic functions of the gas-liquid interface are constructed based on normal curvature and normal vorticity components; By calculating the Laplacian operator and normal vector of the density gradient field, the principal curvature is obtained on the interface tangent plane, and finally the mean curvature expression with mesh independence is obtained. The vorticity vector is obtained by calculating the velocity gradient in each direction of the liquid phase velocity field, then obtaining the velocity field curl to get the vorticity vector, and finally using the time averaging method to obtain a stable characteristic vorticity.

2. The numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to claim 1, characterized in that, A density gradient field is established in the gas-liquid interface region. Acoustic emission detection is used to identify bubble characteristics and establish a gas phase volume fraction distribution function, which is then used to calculate the gas-liquid two-phase mixing density. The spatial rate of change of the mixed density is calculated using the central difference scheme. The interface region is identified by combining the magnitude and direction vector of the density gradient with the eigenvalues ​​of the second derivative tensor. Finally, a density gradient field characterizing the interface intensity and direction is constructed.

3. The numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to claim 2, characterized in that, The gas-liquid two-phase mixing density is calculated based on the phase distribution of the calculation point determined by the gas phase volume fraction. The two-phase region is weighted and the gas phase density is corrected by pressure. The mixing density is obtained by weighting the contribution values ​​of the gas and liquid phases.

4. The numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to claim 3, characterized in that, The gas-liquid relative velocity is corrected by coupling the dynamic characteristic function with the pressure gradient and vorticity field respectively, and by combining the direction modulation factor to calculate the correction coefficients for the rising velocity and lateral drift velocity of the bubble motion, and then using the correction coefficients to correct the gas-liquid relative velocity.

5. The numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction according to claim 4, characterized in that, The turbulent strain rate transport equation is shown below: ； in, For turbulent strain rate tensor, For time, For the average velocity component, For fluid density, For pressure, Kinematic viscosity, For the generation of relative motion turbulence, For interfacial stress, To add turbulent dissipation terms, This represents the partial derivative operator with respect to the spatial coordinate k in the direction of k. Let i represent the partial derivative operator with respect to the direction of spatial coordinate i. This represents the partial derivative operator with respect to the spatial coordinate j in the direction of j. It is the Laplace operator.

6. A numerical simulation system for cavitation turbulence based on gas-liquid slip velocity correction, using the method described in any one of claims 1-5, characterized in that, include: The data acquisition and differentiation module is used to acquire the gas phase volume fraction distribution and liquid phase velocity field distribution in the fluid pipeline, and to establish the density gradient field of the gas-liquid two-phase interface region based on the gas phase volume fraction distribution. The dynamic characteristic function construction module is used to construct a dynamic characteristic function of the gas-liquid interface using the density gradient field and liquid phase velocity field distribution. The dynamic characteristic function includes a coupled expression of the interface curvature term and the local vorticity term. The correction module is used to correct the relative velocity of gas and liquid based on the numerical distribution of the dynamic characteristic function using an adaptive gas-liquid slip velocity correction model. The parameter solving module is used to substitute the modified gas-liquid relative velocity into the turbulent strain rate transport equation to solve for the turbulent characteristic parameters and pressure field distribution in the cavitation flow field.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, it implements the steps of the numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction as described in any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by the processor, it implements the steps of the numerical simulation method for cavitation turbulence based on gas-liquid slip velocity correction as described in any one of claims 1 to 5.

Images