Smooth particle hydrodynamic model pressure intensity smoothing method based on distance and direction information

Abstract

The invention discloses a smoothed particle hydrodynamic model pressure intensity smoothing method based on distance and direction information, and the method comprises the steps: constructing a space correlation kernel function based on the distance between particles and the smooth length, and enabling the space correlation kernel function to be used for quantifying the space correlation degree between the particles; on the other hand, a directional correction kernel function considering the included angle between the particle position vector and the gravity field direction is constructed to reduce excessive smoothness in the gravity direction, so that the linear distribution rule of hydrostatic pressure in the gravity direction is maintained. And performing product combination on the two types of kernel functions to form a composite kernel function capable of reflecting the spatial distribution characteristics and the direction characteristics of the fluid at the same time. By means of the method, efficient correction of the water phase density field in a calculation area can be achieved, the density gradient characteristic under the theoretical hydrostatic pressure effect is kept while numerical oscillation is restrained, then the calculation stability and result accuracy of the model are obviously improved, and the method can be popularized and applied to various types of weakly compressible smoothed particle hydrodynamic models (WCSPH) and has good application prospects. And a technical support is provided for water disaster prevention and control research.

Description

Technical Field

[0001] The present invention belongs to the field of computational fluid dynamics and relates to an innovative smoothing method for smooth particle fluid dynamics models. Specifically, the invention relates to a Shepard mean pressure smoothing method based on distance and direction information, which can be used to improve the computational stability and accuracy of weakly compressible smooth particle fluid dynamics models. Background Art

[0002] Computational fluid dynamics (CFD) plays an irreplaceable role in water disaster prediction and prevention. By numerically solving the Navier-Stokes equations, fluid dynamics models can accurately simulate the hydrodynamic behavior of complex environments (such as rainstorms, floods, dam-break surges, debris flows, landslides, and tsunamis), revealing the physical mechanisms of disaster formation and evolution. With the development of CFD technology, fluid dynamics models have provided a crucial scientific basis for water disaster prediction and prevention. Current CFD models can be classified according to the governing equations (such as the Navier-Stokes equations, shallow water equations, and Saint-Venant equations), the physical processes of interest (pure fluid models, fluid-solid coupling models), and the discretization methods (grid models, meshless models). Among them, meshless models have attracted much attention due to their natural advantages in simulating violently fluctuating water flows (such as the strongly nonlinear free surface flows common in water disasters). Common meshless methods include the Moving Particle Semi-implicit (MPS) method, the Material Point Method (MPM), the Lattice Boltzmann Method (LBM), and the Smoothed Particle Hydrodynamics (SPH) model. The SPH model has attracted widespread attention due to its ease of application. Based on the pressure calculation method, SPH models are further divided into incompressible SPH (ISPH) and weakly compressible SPH (WCSPH) models. In the ISPH model, the fluid is assumed to be incompressible, so its density is always constant. Fluid pressure is obtained by solving the Pressure Poisson Equation (PPE), which is computationally expensive. In the WCSPH model, the fluid is assumed to be weakly compressible, so the fluid density is a variable. Fluid pressure is obtained by solving the explicit Equation of State (EOS). The WCSPH model has been widely developed and applied due to its ease of parallel computation, which significantly improves computational efficiency. However, the WCSPH model suffers from numerical oscillations, which affect computational stability and accuracy. Therefore, employing appropriate pressure smoothing methods is essential for the WCSPH model and is a crucial prerequisite for accurately predicting flood disasters. Summary of the Invention

[0003] Aiming at the numerical oscillation problem existing in the existing weakly compressible smoothed particle hydrodynamics model (WCSPH), the present invention proposes a pressure smoothing method for the smoothed particle hydrodynamics model based on distance and direction information to eliminate numerical oscillation and improve calculation stability and accuracy.

[0004] This method innovatively introduces a dual-kernel function coupling mechanism: first, a distance kernel function with particle distance and smoothing length as variables is constructed to characterize the spatial correlation between particles; simultaneously, a direction kernel function with the angle between the particle relative position vector and the gravity vector as a variable is established to weaken the smoothing in the direction of gravity, thereby preserving the linear distribution characteristics of the hydrostatic pressure. By multiplying these two kernel functions, a composite kernel function is constructed that comprehensively reflects the spatial distribution characteristics and directional information of the fluid. This composite kernel function can effectively smooth the water phase density distribution field within the computational domain, correct the water phase density, and maintain the density gradient distribution characteristics under the action of theoretical hydrostatic pressure while eliminating numerical oscillations, thereby significantly improving the computational stability and numerical accuracy of the model.

[0005] The proposed pressure smoothing method based on distance and direction information can be coupled to other WCSPH models to eliminate numerical oscillations and improve computational stability and accuracy. This method was coupled to a water-sand two-phase flow WCSPH model (Shi, H., Si, P., Dong, P. Yu, X. (2019). A two-phase SPH model for massive sediment motion in free surface flows. Advanced Water Resources, 129, 80-98.) to verify the use and effectiveness of the pressure smoothing method.

[0006] The present invention provides a pressure smoothing method for a smooth particle fluid dynamics model based on distance and direction information, which includes the following steps:

[0007] Step 1: Correct the water phase density based on distance and direction information

[0008] (1) For the target particle "a", search for the adjacent particle "b" in its kernel function domain (see Figure 1 Neighboring particles are defined as all particles within a radius of 2h from the target particle "a." Here, h is the smoothing length, which is set to 1.3Δ. Δ is the initial SPH particle spacing, which depends on a combination of considerations, including sensitivity analysis of the target case and convergence criteria.

[0009] The quotation marks in particle "a" and particle "b" indicate that ab is a code, not a variable, and the quotation marks are used to distinguish it from the variable.

[0010] (2) Calculate the relative distance r between the target particle "a" and the adjacent particle "b" ab , and calculate the particle distance kernel function W accordingly ab .

[0011]

[0012] Where r ab is the relative distance between the target particle "a" and the adjacent particle "b"; h is the smoothing length.

[0013] (3) Calculate the cosine value of the angle between the particle's relative distance vector and the gravity vector, and calculate the direction kernel function W based on it θ .

[0014]

[0015]

[0016] Where r ab is the relative distance vector between the target particle "a" and the adjacent particle "b"; g is the gravity vector; e is a natural constant; σ is the coefficient that controls the degree of smoothness, and it is recommended to take σ = 0.2.

[0017] (4) Multiply the two kernel functions to obtain the composite kernel function W ab W θ .

[0018] (5) Correction for water phase density.

[0019] According to the Shepard weighted average principle, the composite kernel function is used as the weighting function to obtain the corrected water phase density of particle "a".

[0020]

[0021] Where, the superscript “filtered” represents the corrected value, the subscript “a” represents the target particle, and the subscript “b” represents the neighboring particle of the target particle; ρ w is the density of the water phase, (ρ w ) b represents the water phase density of the adjacent particle “b” before updating; W ab is the distance kernel function; W θ is the direction kernel function; V b is the volume of the adjacent particle “b”, V b =(m w ) b / (α w ρ w ) b ,(m w )b is the mass of the aqueous phase of particle “b”, (α w ρ w ) b =(α w ) b ×(ρ w ) b , is the mass density of the aqueous phase of particle "b", where (α w ) b is the water phase volume fraction of particle "b"; the summation symbols in the numerator and denominator indicate that the summation of all adjacent particles "b" of the target particle "a" is performed to obtain the corrected water phase density of the target particle "a"

[0022] Step 2: Apply the pressure smoothing scheme of step 1 to the single-phase flow WCSPH, two-fluid WCSPH, or water-sand two-phase flow WCSPH model to perform pressure smoothing:

[0023] First, according to the initial variable field, the control equations and state equations of each model are solved according to the SPH model strategy, and the time iterative calculation is performed according to the prediction and correction method to obtain the basic variables of each phase carried by all SPH particles until the n×Nth s Numerical iteration steps. According to step 1 (1) to (5), the water phase density is corrected, and the corrected water phase density is substituted for the value before correction (i.e., the variable is updated). The updated water phase density is then brought into the control equation and the iterative calculation is continued. That is, every N s The water phase density is corrected and updated by numerical iteration steps, and the cycle continues until the calculation is completed. Figure 2 .

[0024] Furthermore, in the above method, step 2 couples the pressure smoothing scheme based on distance and direction information to a water-sand two-phase flow WCSPH model (Shi, H., Si, P., Dong, P. Yu, X. (2019). A two-phase SPH model for massive sediment motion in free surface flows. Advanced Water Resources, 129, 80-98.) to perform pressure smoothing. The specific process is as follows:

[0025] (2-1) First, according to the conventional SPH strategy, the time iteration calculation is performed according to the prediction and correction method, see Figure 2 , the control equations are as follows (5) to (8), and the state equation is as follows (9):

[0026]

[0027] In the formula, the variables in the iterative calculation refer to the variables of the target particle "a", u w is the water phase velocity; u s is the sediment phase velocity; p w is the water phase pressure; τ w is the water phase shear force; F is the volume force; g is the acceleration due to gravity; τ s is the sediment phase shear force; p s is the sediment phase pressure; α w is the volume fraction of the water phase; α s is the volume fraction of sediment phase; ρ w0 is the density of water phase at standard atmospheric pressure, ρ w0 =1000kg / m 3 ; c0 is the numerical sound velocity, which is generally taken as 10 times the maximum velocity of the water phase; ζ is the model parameter, and ζ = 7 is recommended.

[0028] (2-2) in the n×Nth s The water phase density term is corrected according to the above step one, that is, steps (1) to (5) in step one are implemented.

[0029] (2-3) Based on the corrected water phase density of the target particle “a”, the mass density terms of the water phase and the sediment phase of the target particle “a” are updated.

[0030]

[0031] Wherein, the subscript "a" represents the variables corresponding to particle a; (α w ρ w ) a =(α w ) a ×(ρ w ) a , is the water phase mass density of target particle “a”; (α w ) a is the water phase volume fraction of the target particle “a”, defined as (α w ) a =(m w ) a / (ρ w ) a V a ;(m w ) a is the mass of the aqueous phase of the target particle "a"; (ρ w ) a is the water phase density of target particle “a”; V a is the volume of the target particle "a"; (α s ) a is the volume fraction of the sediment phase of the target particle “a”, defined as (αs ) a =(m s ) a / ρ s V a ;(m s ) a is the sediment phase mass of target particle “a”; ρ s is the sediment phase density, which depends on the physical properties of the sediment and is a constant.

[0032] (2-4) The updated water phase mass density term (α w ρ w ) filtered and sediment phase density term Replace α in the original control equation and state equation respectively w ρ w and α s , continue the normal iterative calculation.

[0033] (2-5) Solve the control equations (5) to (8) and the state equation (9) according to the conventional SPH strategy, and perform time-based iterative calculations according to the prediction and correction method to obtain the basic variables of the two phases. s In the numerical iteration step, the water phase density is modified again according to steps (1) to (5) in step 1, and the variables and iterative calculations are updated according to steps (2-1) to (2-5). s The water phase density is corrected and updated once in each iteration, and the cycle continues until the calculation is completed. s The smaller it is, the higher the correction frequency is. Generally, N is recommended. s Take 10 to 50.

[0034] Furthermore, in the above method, in step (2), the particle distance kernel function W ab Relative distance r ab The relationship between the smooth length h is as follows Figure 3 (a) W ab Follow r ab decreases with the increase of r ab >2h, then W ab =0.

[0035] Furthermore, in the above method, in step (3), the directional kernel function W θ The relationship between the angle and the vector is as follows Figure 3 (b) (σ=0.2). When the vector angle θ is close to π / 2 and 3π / 2, the directional kernel function W θ Reaching the maximum value of 1.0, the distance vector r abPerpendicular to the gravity vector. This is equivalent to strengthening the influence of the adjacent particle "b" on the target particle "a" in the direction perpendicular to gravity. When the vector angle θ is close to 0, π and 2π, the direction kernel function W θ Reaching a minimum value of 3.7×10 -6 , at this time the distance vector r ab Parallel to the gravity vector. This effectively reduces the influence of neighboring particle "b" on target particle "a" in the direction parallel to gravity. This maintains the natural gradient of hydrostatic pressure in the direction of gravity and avoids numerical errors caused by oversmoothing in the direction of gravity.

[0036] Furthermore, in the above method, in step (4), the composite kernel function W ab W θ The relationship with the independent variable is as follows Figure 3 (c). When r ab / h=0,W ab ×W θ It reaches its maximum value at θ = π / 2 and θ = 3π / 2, when the adjacent particle "b" and the target particle "a" coincide in space. ab When / h is a non-zero constant, W ab W θ As θ approaches 0, π, and 2π, it decreases; W ab W θ As θ approaches π / 2 and 3π / 2, it increases. When the vector angle θ is constant, W ab W θ Follow r ab / h increases and decreases.

[0037] Furthermore, in the above method, when the particle exceeds the calculation domain of the composite kernel function, that is, r ab >2h, at this time W ab =W θ =W ab W θ =0.

[0038] Compared with the prior art, the present invention has the following beneficial effects:

[0039] 1. The method of the present invention proposes a pressure smoothing method for a smoothed particle hydrodynamic model based on distance and direction information, which can achieve efficient smoothing of the water phase density field in the calculation area. While suppressing numerical oscillations, it maintains the density gradient characteristics under the action of theoretical hydrostatic pressure, thereby significantly improving the pressure calculation stability and result accuracy of the model. It can be extended to various types of weakly compressible smoothed particle hydrodynamic models (WCSPH), providing technical support for water disaster prevention and control research.

[0040] 2. The method of the present invention can preserve the linear distribution of hydrostatic pressure in the direction of gravity, avoiding the destruction of the linear distribution of hydrostatic pressure due to excessive smoothing, thereby ensuring the stability and accuracy of the calculation.

[0041] 3. The method of the present invention can avoid the water surface elevation calculation error caused by the destruction of the linear distribution of hydrostatic pressure due to excessive smoothing. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 Schematic diagram of the kernel function calculation domain;

[0043] Figure 2 Flowchart of implementation steps;

[0044] Figure 3 is the relationship diagram between kernel function and independent variable. (a) is the distance kernel function W ab ; (b) Directional kernel function W θ ; (c) Composite kernel function W ab W θ .

[0045] Figure 4 Comparison of the calculation results of three pressure smoothing schemes. (ac) is the calculation result without considering the pressure smoothing scheme; (df) is the calculation result considering only the distance kernel function W ab The calculation results of the pressure smoothing scheme (only W is used in Equation 4) ab is the weighting function); (gi) is the composite kernel function W ab W θ Results of the pressure smoothing method. DETAILED DESCRIPTION

[0046] The pressure smoothing method of the smoothed particle fluid dynamics model based on distance and direction information of the present invention will be further described below in conjunction with specific embodiments.

[0047] Example 1

[0048] A water tank with a length of 0.4m, a width of 0.2m and a height of 0.4m is filled with water with a depth of h. w = 0.3m of tap water. When the water is completely still, the pressure distribution at this time conforms to the theoretical hydrostatic pressure, that is, p w =ρ w gh. It presents a triangular distribution in the direction of gravity, the relative pressure at the water surface is 0Pa, and the relative pressure at the bottom of the tank is p w =ρ w gh w =2943Pa.

[0049] The pressure smoothing method of the smoothed particle hydrodynamic model based on distance and direction information proposed in this paper is coupled to a water-sand two-phase flow WCSPH model (Shi, H., Si, P., Dong, P. Yu, X. (2019). A two-phase SPH model for massive sediment motion in free surface flows. Advanced Water Resources, 129, 80-98.) to simulate the pressure distribution in the water tank. Since there is only pure water and no sediment in this case, α w =1.0 and α s = 0.0. The model only calculates the two-dimensional pressure distribution in the length-height plane, so periodic boundary conditions are used in the width direction. Solid boundary conditions are used on the side walls and bottom of the tank. Based on sensitivity analysis, the initial spacing of the SPH particles in this case is Δ = 5 mm. Considering the convergence condition, the calculation step size is Δt = 10 -5 s. The scheme without considering pressure smoothing, the scheme only considering distance kernel function and the scheme considering composite kernel function are coupled to the water-sand two-phase flow WCSPH model respectively. According to the initial variable field and the conventional SPH strategy (Shi, H., Yu, X., & Dalrymple, RA (2017). Development of a two-phase SPH model for sediment laden flows. Computer Physics Communications, 221, 259-272.), the control equations (Equations (5) to (8) and state equation (9)) are solved. The prediction and correction method (Monaghan, JJ (1989). On the problem of penetration in particle methods. Journal of Computational physics, 82 (1), 1-15.) is used for time iteration. Every 20 iteration steps (N s =20) performs a water phase density correction process (steps (1) to (5)), and updates the variables and iterates the calculation accordingly (steps (2-1) to (2-5)), and the calculation results of the pressure distribution in the water tank can be obtained.

[0050] The results are as follows Figure 4 As shown in Figure 2, (a) to (c) represent the calculation results without considering the pressure smoothing scheme; (d) to (f) represent the calculation results when only the particle distance kernel function W is considered in formula (4). abis the calculation result of the pressure smoothing scheme of the weighted function; (g) to (i) represent the calculation results of the composite kernel function pressure smoothing method based on distance and direction information. It can be clearly seen that when no pressure smoothing scheme is used, the calculated water phase pressure distribution is obviously chaotic and seriously deviates from the theoretical linear hydrostatic pressure distribution. When the kernel function W that only considers the particle distance is used, the calculated water phase pressure distribution is obviously chaotic and seriously deviates from the theoretical linear hydrostatic pressure distribution. ab When the pressure smoothing scheme is used, the water pressure distribution is smooth and basically conforms to the linear distribution, but the water surface pressure exceeds 0Pa, causing the SPH particles on the water surface to rise erroneously, resulting in a large error in the calculation of the free water surface. In addition, the pressure at the bottom of the water tank is significantly lower than the theoretical value. At t=20s, the relative error in the calculated pressure at the bottom of the water tank is approximately -29.9%. On the contrary, the use of the composite kernel function pressure smoothing method can effectively eliminate the numerical oscillation of the WCSPH model, and the calculated pressure distribution is smooth and in good agreement with the theoretical hydrostatic pressure. Because the directional information between the distance vector and the gravity vector is taken into account, the linear distribution of the hydrostatic pressure is well preserved, avoiding excessive smoothing in the direction of gravity. In addition, the composite kernel function pressure smoothing method proposed in this invention also ensures that the relative pressure at the free water surface is near 0Pa, avoiding errors in the calculation of the free water surface. In summary, the invention proposes a pressure smoothing method for a smoothed particle fluid dynamics model based on distance and direction information, which improves the computational stability and accuracy of the WCSPH model, with significant results.

Claims

1. A pressure smoothing method for a smooth particle fluid dynamics model based on distance and direction information, characterized in that: Includes the following: Step 1: Correct the water phase density based on distance and direction information (1) For the target particle "a", the neighboring particles "b" in the kernel function domain are retrieved. The neighboring particles are defined as all particles within a radius of 2h of the target particle "a", where h is the smoothing length, h = 1.3Δ, and Δ is the initial SPH particle spacing, which depends on comprehensive considerations such as the sensitivity analysis of the target case and the convergence criterion; (2) Calculate the relative distance r between the target particle "a" and the adjacent particle "b" ab , and calculate the particle distance kernel function W accordingly ab , Where r ab is the relative distance between the target particle "a" and the adjacent particle "b"; h is the smoothing length; (3) Calculate the cosine value of the angle between the particle's relative distance vector and the gravity vector, and calculate the direction kernel function W based on it θ , Where r ab is the relative distance vector between the target particle "a" and the adjacent particle "b"; g is the gravity vector; σ is the coefficient controlling the smoothness, σ = 0.2; (4) Multiply the two kernel functions to obtain the composite kernel function W ab W θ ; (5) Correction of water phase density: According to Shepard's weighted average principle, the composite kernel function is used as the weighting function to obtain the corrected water phase density of particle "a". Where, the superscript "filtered" represents the corrected value, the subscript "a" is the target particle, and the subscript "b" is the neighboring particle of the target particle; ρ w is the density of the water phase, (ρ w ) b represents the water phase density of the neighboring particle "b" before updating; W ab is the distance kernel function; W θ is the direction kernel function; V b is the volume of the adjacent particle "b", V b =(m w ) b / (α w ρ w ) b ,(m w ) b is the mass of the aqueous phase of particle "b", (α w ρ w ) b =(α w ) b ×(ρ w ) b , is the mass density of the aqueous phase of particle "b", where (α w ) b is the water phase volume fraction of particle "b"; the summation symbols in the numerator and denominator indicate that the summation of all adjacent particles "b" of the target particle "a" is performed to obtain the corrected water phase density of the target particle "a" Step 2: Apply the pressure smoothing scheme of step 1 to the single-phase flow WCSPH, two-fluid WCSPH, or water-sand two-phase flow WCSPH model to perform pressure smoothing: First, according to the initial variable field, the control equations and state equations of each model are solved according to the SPH model strategy, and the time iterative calculation is performed according to the prediction and correction method to obtain the basic variables of each phase carried by all SPH particles until the n×Nth s Numerical iteration steps; according to step 1 (1) to (5), the water phase density is corrected, and the corrected water phase density is substituted for the value before correction, that is, the variable is updated, and the updated water phase density is brought into the control equation, and the iterative calculation is continued, that is, every N s The water phase density is corrected and updated through numerical iteration steps, and the operation is repeated until the calculation is completed.

2. The method according to claim 1, characterized in that Step 2: The specific process of coupling the pressure smoothing scheme based on distance and direction information into the water-sand two-phase flow WCSPH model to perform pressure smoothing is as follows: (2-1) First, according to the conventional SPH strategy, the time iteration calculation is performed according to the prediction correction method. The control equations are as follows: (5) to (8), and the state equation is as follows: In the formula, the variables in the iterative calculation refer to the variables of the target particle "a", u w is the water phase velocity; u s is the sediment phase velocity; p w is the water phase pressure; τ w is the water phase shear force; F is the volume force; g is the acceleration due to gravity; τ s is the sediment phase shear force; p s is the sediment phase pressure; α w is the volume fraction of the water phase; α s is the volume fraction of sediment phase; ρ w0 is the density of water phase at standard atmospheric pressure, ρ w0 =1000kg / m 3 ; c0 is the numerical sound velocity, which is 10 times the maximum velocity of the water phase; ζ is the model parameter, which is ζ = 7; (2-2) in the n×Nth s Numerical iteration steps are performed to correct the water phase density term according to the above step 1, i.e., steps (1) to (5) in step 1 are implemented; (2-3) Based on the corrected water phase density of the target particle "a", the mass density terms of the water phase and sediment phase of the target particle "a" are updated; Wherein, the subscript "a" represents the variables corresponding to particle a; (α w ρ w ) a =(α w ) a ×(ρ w ) a , is the mass density of the target particle "a" in the water phase; (α w ) a is the volume fraction of the water phase of the target particle "a", defined as (α w ) a =(m w 0 a (ρ w 0 a V a ;(m w ) a is the mass of the aqueous phase of the target particle "a"; (ρ w ) a is the water phase density of target particle "a"; V a is the volume of the target particle "a"; (α s ) a is the volume fraction of the sediment phase of the target particle "a", defined as (α s ) a =(m s ) a / ρ s V a ;(m s ) a is the sediment phase mass of target particle "a"; ρ s is the sediment phase density, which depends on the physical properties of the sediment and is a constant; (2-4) The updated water phase mass density term (α w ρ w ) filtered and sediment phase density term Replace α in the original control equation and state equation respectively w ρ w and α s , continue the conventional iterative calculation; (2-5) Solve the control equations (5) to (8) and the state equation (9) according to the conventional SPH strategy, and perform time-based iterative calculations according to the prediction and correction method to obtain the basic variables of the two phases until the next n×N s In the numerical iteration step, the water phase density is modified again according to steps (1) to (5) in step 1, and the variables and iterative calculations are updated according to steps (2-1) to (2-5). s The water phase density is corrected and updated once in an iterative step, and the cycle continues until the calculation is completed. s The smaller it is, the higher the correction frequency is. s Take 10 to 50.

Images