A stretch-stable sph fluid simulation method

Abstract

The present application relates to a kind of stretch stable SPH fluid simulation method, can effectively weaken stretch instability problem in SPH fluid simulation process.The present application introduces two new particles in SPH fluid simulation process, namely: as the sampling point of pressure field virtual particle, and with speed attribute and carry fluid mass real particle, and based on these two particles, it is constructed based on the projection method of SPH fluid simulation technology method.In addition, the present application provides two virtual particle generation methods, to efficiently generate virtual particles without increasing too much computing expense.The present application can effectively weaken stretch instability problem in the original SPH fluid simulation method, significantly improve the stability, visual quality and accuracy of fluid simulation.The present application not only can simulate the form of fluid in stretch state, such as thin layer, fine stream, splash and droplet, but also can realize the stable simulation of incoherent fluid.

Description

Technical Field

[0001] This invention belongs to the field of computer graphics, specifically relating to a method for simulating tensile-stabilized SPH fluids. Background Technology

[0002] Smoothed Particle Hydrodynamics (SPH) is a widely used fluid simulation method. This method has important applications in film and television special effects, game development, computer-aided engineering, and marine engineering. After many years of development, this method can now achieve rapid and accurate simulations of various fluid materials and multiphase flow scenarios on computers. However, since its invention in 1977, the SPH fluid simulation method has been plagued by the problem of tensile instability. This makes it difficult for SPH to accurately simulate the tensile state of fluids and to realistically simulate the detailed features of free surface flows such as thin layers, filaments, and splashes. This problem severely limits the application of SPH-based simulation methods.

[0003] The main reasons for the tensile instability problem in the SPH fluid simulation method are as follows: (1) The particle position update calculation in the SPH method is realized from the Lagrange perspective, but the related calculation of the interaction between particles must use the SPH kernel function from the Euler perspective. This will cause a group of particles with negative pressure to gradually increase the adsorption force as the distance between them decreases, making it impossible for the particles to maintain a stable tensile state; (2) The particles in the SPH method are both sampling points of physical fields such as velocity and sampling points of their differential fields. The overlap of these two sampling points will introduce the "zero-energy modes" and affect the accuracy of the calculation, which will further aggravate the tensile instability problem; (3) The calculation accuracy of the differential fields commonly used in the SPH method is not high enough, which will also introduce non-physical disturbances and exacerbate the tensile instability problem. Summary of the Invention

[0004] The technical problem solved by this invention is to overcome the shortcomings of the prior art and provide a technical method that can weaken the tensile instability problem in SPH fluid simulation, significantly improve the stability of the SPH fluid simulation method, and efficiently improve the visual quality of SPH fluid simulation without introducing too much computation.

[0005] The specific technical solution of the present invention is as follows: A method for simulating tensile-stabilized SPH fluid, comprising the following steps:

[0006] Step 1: In the initialization phase before simulation, real particles representing the fluid are generated. Each real particle carries part of the fluid's mass and has density and velocity properties.

[0007] Step 2: Apply viscous and volume forces to the real particles to update their velocity;

[0008] Step 3: Generate the position of virtual particles based on the position of real particles. Virtual particles only carry pressure properties and cannot move.

[0009] Step 4: Find the neighboring real particles and neighboring virtual particles for each real particle and virtual particle;

[0010] Step 5: Discretize the Poisson pressure equation using virtual particles, real particles, and their neighboring real particles and neighboring virtual particles to obtain a set of linear equations for solving the pressure of virtual particles, and solve for the pressure of virtual particles.

[0011] Step 6: Discretize the fluid Euler equations using virtual particles, real particles, and their neighboring real particles and virtual particles. Then calculate the pressure gradient of each real particle based on the pressure of the virtual particles, and update the velocity of the real particles using the pressure gradient of the real particles.

[0012] Step 7: Delete all virtual particles, and then update the positions of the real particles based on their velocities;

[0013] Step 8: Repeat steps 2-7 until the fluid simulation task is completed.

[0014] Furthermore, in step 3, virtual particles are generated based on the positions of real particles. There are two ways to generate the positions of virtual particles. The following are the implementation steps of the first method for generating virtual particle positions:

[0015] (1) Perform a deep copy of the position data of the real particles, and use the data after deep copying as the initial position data of the virtual particles;

[0016] (2) Correct the initial position of the virtual particles so that the position distribution of the virtual particles tends to be uniform;

[0017] The following are the implementation steps for the second method of generating virtual particle positions:

[0018] (1) Define orthogonal, equally spaced spatial points in the simulation space as anchor points for virtual particle generation;

[0019] (2) Generate candidate points for virtual particles at anchor points near each real particle;

[0020] (3) Convert all candidate points of virtual particles into Morton codes and store them in an array A;

[0021] (4) Sort array A;

[0022] (5) Perform deduplication on array A, that is: check whether each element in array A is the same as its neighboring element. If they are the same, remove the duplicate elements.

[0023] (6) Convert each element in the deduplicated array A into spatial coordinates, which are the positions of the virtual particles.

[0024] Furthermore, in step 5, the pressure Poisson equation is discretized using virtual particles, real particles, and their neighboring real particles and neighboring virtual particles, wherein the continuous form of the pressure Poisson equation is:

[0025]

[0026] In the above formula This refers to the velocity vector field of the fluid after the application of viscous and volume forces. The pressure field of the fluid; For the fluid density field, The time step is the time interval between two simulated frames;

[0027] The discretized left-hand side of the pressure Poisson equation takes the following form:

[0028]

[0029] The discretized source terms of the pressure Poisson equation, i.e., the right-hand side terms of the pressure Poisson equation, are in the following form:

[0030]

[0031] In the above two formulas, and The label for the virtual particle. and The designation for real particles; The first derivative of the kernel function in the SPH method; and virtual particles and real particles The distance; and Virtual particles and real particles The volume of space it occupies; virtual particles The pressure; For real particles that may not satisfy the divergence-free velocity condition; virtual particles With real particles The relative positions between them; The density of the fluid in its static state. The SPH support domain radius is the neighborhood radius divided by the smoothness length. The velocity at the virtual particle's position.

[0032] Furthermore, in step 6, the Euler equations for the fluid are discretized using virtual and real particles. The discretized Euler equations are in the following specific form:

[0033]

[0034] in, For real particles that satisfy the incompressibility condition speed; real particles With virtual particles The relative positions between them; real particles With virtual particles The relative distance between them.

[0035] Furthermore, in step 4, finding the neighboring particles of each particle specifically includes:

[0036] ① The virtual particle neighbors of each virtual particle;

[0037] ②The real particle neighbors of each virtual particle;

[0038] ③ The virtual particle neighbors of each real particle;

[0039] ④ The real particle neighbors of each real particle.

[0040] Furthermore, the velocity at the virtual particle's position The calculation formula is:

[0041]

[0042] In the above formula, This refers to the kernel function in the SPH method.

[0043] Furthermore, to avoid the accumulation of density errors during the simulation, the source term, i.e., the right-hand side term, of the pressure Poisson equation is compensated using the following formula:

[0044]

[0045] in, The source terms of the Poisson equation before compensation, These are the compensated source terms of the Poisson equation. The density value calculated using the real particle neighborhood of the virtual particle is: , real particles quality This is an empirical coefficient, with a value of 1.

[0046] Compared with existing SPH fluid simulation methods, this invention has the following advantages:

[0047] (1) In this invention, the virtual particles that store pressure and velocity divergence, and the real particles that store velocity and pressure gradients can be separated from each other, that is, the sampling points of the physical field and its differential field no longer coincide, thus effectively weakening the influence of the "zero energy" problem.

[0048] (2) The real particles in this invention do not carry pressure properties, so the adsorption force between real particles will not increase incorrectly as the position decreases. Therefore, visual artifacts such as non-physical aggregation of particles are not easily generated when simulating fluid stretching.

[0049] (3) When the fluid presents a thin layer, filament, etc., since there are sampling points of uniform and sufficient pressure field inside and outside the fluid - virtual particles, the real particles can stably and relatively accurately calculate the pressure and update their own velocity, thereby ensuring the uniform distribution of real particles and stably maintaining the thin layer, filament, etc. form of real particles.

[0050] (4) In this invention, the virtual particles storing pressure are evenly and sufficiently distributed, which can effectively improve the convergence speed of the iteration and improve the stability of the simulation calculation.

[0051] (5) The present invention can stably simulate fluids while retaining negative pressure. The presence of negative pressure will keep the velocity and direction of particle motion close. Therefore, the present invention can stably simulate inviscid fluids.

[0052] (6) In this invention, all steps, including virtual particle generation, are easy to implement in parallel. With the help of hardware devices such as GPUs, the time cost of the fluid simulation process can be further reduced.

[0053] (7) Compared with the SPH method which uses "stress point", there is no topological relationship or strict positional dependence between the virtual particle position and the real particle position in the method presented in this paper, so it is very suitable for simulating large deformation materials such as fluids. Attached Figure Description

[0054] Figure 1 This is a schematic diagram of the steps in the method for realizing fluid incompressibility in this invention.

[0055] Figure 2 This is a schematic diagram of virtual particles obtained using virtual particle generation methods (I) and (II).

[0056] Figure 3 This is a schematic diagram illustrating the implementation steps of the virtual particle generation method (I) in this invention.

[0057] Figure 4 This is a schematic diagram illustrating the implementation steps of the virtual particle generation method (II) in this invention. Detailed Implementation

[0058] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be described in detail below with reference to specific implementations and accompanying drawings, but this does not constitute a limitation on the present invention.

[0059] This invention provides a dual-particle method to address the stretch instability problem in SPH fluid simulation. It introduces two new particle types into the SPH fluid simulation method: virtual particles serving as pressure field sampling points, and real particles possessing velocity attributes and carrying fluid mass. Based on these two particle types, a fluid simulation technique based on the SPH projection method is constructed. Simultaneously, this invention provides two virtual particle generation methods to efficiently generate virtual particles without significantly increasing computational overhead. This invention effectively overcomes the stretch instability problem of the original SPH fluid simulation method, significantly improving the stability and visual quality of fluid simulation. It can realistically simulate thin layers, streams, splashes, and droplets of viscous fluids under stretching conditions on a computer, while ensuring that the SPH method can stably simulate inviscid fluids. Furthermore, this invention can also effectively improve the simulation quality of highly viscous fluids.

[0060] 1. The hardware platform of the method of this invention uses an Intel i7-8700K six-core CPU with a main frequency of 3.7GHz and an NVIDIA GeForce GTX 2080 graphics card with 8GB of video memory. The system program is written in C++, with CUDA language used for acceleration of the parallel computing part, and is compiled and executed with the help of Microsoft Visual Studio 2017. Open source libraries such as OpenGL and Freeglut are used in the development process.

[0061] 2. This invention uses the Euler equations for fluids and the divergence-free condition of velocity as the dynamic basis for the fluid simulation process, and its specific form is as follows.

[0062]

[0063]

[0064] By discretizing the time derivative of velocity in Euler's formula for fluids and substituting it with the divergence-free condition for velocity, we can obtain the Poisson equation for fluid pressure.

[0065]

[0066] 3. In the initialization phase before simulation, real particles representing the fluid are first generated. Each real particle carries part of the fluid's mass and has density and velocity properties.

[0067] 4. After the simulation starts, viscous and volume forces need to be applied to the real particles in the initial stage of each frame of calculation, and the velocity of the real particles needs to be updated.

[0068] 5. This invention uses the SPH method to discretize each physical field. Virtual particles and real particles are used in the discretization process. Virtual particles only carry pressure properties and cannot move. Virtual particles need to be generated before performing incompressibility calculations and deleted after the incompressibility calculations are completed. Figure 2 This is a schematic diagram of two different methods for generating virtual particles. Figure 3 The following are the implementation steps of virtual particle generation method (I), which is used to generate virtual particles with a relatively uniform distribution. The implementation method of virtual particle generation method (I) is as follows:

[0069] ① Obtain the position of the real particle ( Figure 3 Step 1 in the middle).

[0070] ② Copy the position data of real particles as the initial position of virtual particles, and find and store the virtual particle neighbors of each virtual particle. Figure 3 Step 2 in the middle).

[0071] ③ Apply the PBF (Position Based Fluid) algorithm or ParticleShifting algorithm to virtual particles to make their positions closer to those of real particles and their distribution more uniform. Figure 3 Step 3 in the middle.

[0072] 6. Figure 4 The implementation steps of virtual particle generation method (II) are as follows: This method can generate orthogonal, equally spaced, and sufficiently numerous virtual particles (such as...) inside and near the fluid. Figure 2 (As shown in the right figure). (Refer to...) Figure 4 The specific implementation method of virtual particle generation method (II) is as follows:

[0073] ① Obtain the position of the real particles, and at the same time define "anchor points" that are evenly spaced in the simulation space. The "anchor points" do not need to be explicitly generated, and their spacing can be 1 to ~1.8 times the sampling spacing of the real particles. Figure 4 (Step 1 and Step 2).

[0074] ② For each real particle, the nearest "anchor points" within the support domain radius are used as virtual particle candidate points, and the coordinates of these virtual particle candidate points are stored contiguously in memory. If two real particles are close together, there may be points with the same coordinates among their neighboring virtual particle candidate points, i.e., duplicate virtual particle candidate points. These duplicate points need to be deleted in subsequent operations.

[0075] ③ Convert the stored candidate point coordinates into Morton codes and store them in array A, and delete the original virtual particle candidate point coordinate data. Figure 4 Step 3 in the middle.

[0076] ④ Sort the array A containing the Morton codes of the candidate virtual particles. Since the Morton codes derived from the coordinates of the same position are also the same, the Morton codes of the candidate points at the same position will be stored in adjacent positions in array A after sorting. Figure 4 Step 4 in the middle.

[0077] ⑤ Remove duplicate elements from array A. The steps are as follows: Check if each element in array A is the same as its adjacent elements, and store the results in an additional array B. Array A and array B have the same size. If an element in array A has no duplicates or is the first of several consecutive identical elements, store a 1 in the corresponding position in B; otherwise, store a 0. After completing the search, count the number of elements with a value of 1 in array B, denoted as n. Create a new array C with size n, and using array B as a reference, store the non-duplicate elements from array A in array C. After completing all the above operations, delete arrays A and B. Figure 4 Step 5 in the middle.

[0078] ⑥ Inversely convert the Morton codes in array C to spatial coordinates, and use these spatial coordinates to generate virtual particles. Figure 4 Step 6 in the middle.

[0079] The above-mentioned method for generating virtual particles is easy to parallelize on multi-threaded computing such as GPUs.

[0080] 7. Find the neighboring real particles and neighboring virtual particles of each real particle and virtual particle;

[0081] 8. This invention uses the SPH method to discretize each physical field. Virtual particles and real particles are used in the discretization process. The specific implementation method is as follows.

[0082] ① Since this invention contains two different types of particles—virtual particles and real particles—the discretization of the physical field in this invention is calculated using the following four equations:

[0083]

[0084]

[0085]

[0086]

[0087] In the above formula, This is the SPH kernel function; For a certain type of physical field; , , and Calculate real particles using real particle neighbors respectively Physical quantities, using virtual particle neighbors to calculate real particles Physical quantities, using virtual particle neighbors to calculate virtual particles Physical quantities, and the calculation of virtual particles using virtual particle neighbors. The physical quantity.

[0088] ② Virtual particles not only carry pressure, but also properties such as velocity divergence and the Laplace's expression for pressure. The Laplace's discretization formula for pressure at the virtual particle's location is:

[0089]

[0090] The formula for discretizing the velocity divergence at the virtual particle's position is:

[0091]

[0092] ③Real particles not only carry velocity but also pressure properties. The pressure gradient discretization method in pressure is as follows:

[0093]

[0094] The final formula for calculating the velocity update of real particles is:

[0095]

[0096] 9. Figure 1 This method for realizing fluid incompressibility can be divided into three steps:

[0097] ① Generate virtual particles based on the positions of real particles ( Figure 1 (Step 1 and Step 2)

[0098] ② Use virtual and real particles to discretize the Poisson equation for pressure and solve for the fluid pressure. Figure 1 Step 3 in the middle).

[0099] ③ Use virtual particle pressure to calculate pressure force and update the velocity and position of real particles. Figure 1 Step 4 in the middle.

[0100] Based on the aforementioned velocity divergence discretization method and pressure Laplace's discretization method, the Poisson equation for pressure can be transformed into a system of linear equations. To ensure that this system of linear equations is nonsingular, free surface boundary conditions (i.e., Diraclee boundary conditions) need to be introduced into the Poisson equations. One possible approach to applying free surface boundary conditions is to use semi-analytical free surface boundary conditions, the specific implementation of which is as follows:

[0101] First, it is necessary to detect whether the virtual particle is located outside the fluid. A virtual particle is considered to be outside the fluid if it meets the following conditions:

[0102]

[0103]

[0104] To solve the coefficient matrix of the linear equation system of pressure, Diagonal elements can also be considered as Virtual particles The attributes; among which, It is calculated from virtual particles within the fluid that have sufficient neighboring particles, i.e. In the above formula It is a relatively small constant, whose value is generally less than the sampling distance of real particles. times.

[0105] When the virtual particle is outside the fluid, the pressure Laplace discretization method in the particle's equation must be adjusted to the following equation:

[0106]

[0107] When virtual particles are located inside the fluid, the pressure Laplace discretization method is as follows:

[0108]

[0109] After correction, the linear equations for solving pressure are transformed into a symmetric positive definite linear system, and the pressure can be quickly solved using the conjugate gradient method or the Jacobi method.

[0110] 10. Although the linear system of pressure can be solved quickly after the above corrections to the coefficients, to ensure the stability of the simulation, the following two corrections are needed for the source terms of the pressure Poisson equation (the velocity divergence term on the right-hand side of the Poisson equation):

[0111] ① Eliminate fluid volume loss (density error). Although solving the above pressure Poisson equation can obtain a velocity field that satisfies the incompressibility condition, during simulation, calculation errors may still introduce density errors, leading to fluid volume loss. Therefore, in the implementation of this invention, the source terms in the original pressure Poisson equation need to be compensated according to the fluid density error. The compensation method is as follows:

[0112]

[0113] in, The source terms of the Poisson equation before compensation, and It is the source term of the Poisson equation after compensation. The density value calculated using the real particle neighborhood of the virtual particle is: . This is an empirical coefficient, and its value is generally 1.

[0114] ② Velocity divergence calculation correction. To avoid significant errors in velocity divergence calculations near the fluid due to the lack of real particle neighbors, the velocity divergence can be discretized as follows:

[0115]

[0116]

[0117] 11. Delete all virtual particles, then update the positions of the real particles based on their velocities, and perform simulation calculations for the next frame.

Claims

1. A simulation method for tensile-stabilized SPH fluid, characterized in that, Includes the following steps: Step 1: In the initialization phase before simulation, real particles representing the fluid are generated. Each real particle carries part of the fluid's mass and has density and velocity properties. Step 2: Apply viscous and volume forces to the real particles to update their velocity; Step 3: Generate the position of virtual particles based on the position of real particles. Virtual particles only carry pressure properties and cannot move. Step 4: Find the neighboring real particles and neighboring virtual particles for each real particle and virtual particle; Step 5: Discretize the Poisson pressure equation using virtual particles, real particles, and their neighboring real particles and neighboring virtual particles to obtain a set of linear equations for solving the pressure of virtual particles, and solve for the pressure of virtual particles. Step 6: Discretize the fluid Euler equations using virtual particles, real particles, and their neighboring real particles and virtual particles. Then calculate the pressure gradient of each real particle based on the pressure of the virtual particles, and update the velocity of the real particles using the pressure gradient of the real particles. Step 7: Delete all virtual particles, and then update the positions of the real particles based on their velocities; Step 8: Repeat steps 2-7 until the fluid simulation task is completed; In step 3, virtual particles are generated based on the positions of real particles. There are two ways to generate the positions of virtual particles. The following are the implementation steps of the first method for generating virtual particle positions: (1) Perform a deep copy of the position data of the real particles, and use the data after deep copying as the initial position data of the virtual particles; (2) Correct the initial position of the virtual particles so that the position distribution of the virtual particles tends to be uniform; The following are the implementation steps for the second method of generating virtual particle positions: (1) Define orthogonal, equally spaced spatial points in the simulation space as anchor points for virtual particle generation; (2) Generate candidate points for virtual particles at anchor points near each real particle; (3) Convert all candidate points of virtual particles into Morton codes and store them in an array A; (4) Sort array A; (5) Perform deduplication on array A, that is: check whether each element in array A is the same as its neighboring element. If they are the same, remove the duplicate elements. (6) Convert each element in the deduplicated array A into spatial coordinates, which are the positions of the virtual particles.

2. The tensile-stabilized SPH fluid simulation method according to claim 1, characterized in that: In step 5, the pressure Poisson equation is discretized using virtual particles, real particles, and their neighboring real particles and neighboring virtual particles. The continuous form of the pressure Poisson equation is: In the above formula, This refers to the velocity vector field of the fluid after the application of viscous and volume forces. The pressure field of the fluid; For the fluid density field, The time step is the time interval between two simulated frames; The discretized left-hand side of the pressure Poisson equation takes the following form: The discretized source terms of the pressure Poisson equation, i.e., the right-hand side terms of the pressure Poisson equation, are in the following form: In the above two formulas and The label for the virtual particle. and The designation for real particles; The first derivative of the kernel function in the SPH method; and virtual particles and real particles The distance; and Virtual particles and real particles The volume of space it occupies; virtual particles The pressure; For real particles that may not satisfy the divergence-free velocity condition; virtual particles With real particles The relative positions between them; The density of the fluid in its static state. The SPH support domain radius is the neighborhood radius divided by the smoothness length. The velocity at the virtual particle's position.

3. The tensile-stabilized SPH fluid simulation method according to claim 2, characterized in that: In step 6, virtual particles and real particles are used to discretize the Euler equations of the fluid. The specific form of the discretized Euler equations is as follows: in, For real particles that satisfy the incompressibility condition speed; real particles With virtual particles The relative distance between them.

4. The tensile-stabilized SPH fluid simulation method according to claim 1, characterized in that: Step 4, which involves finding the neighboring particles of each particle, specifically includes: ① The virtual particle neighbors of each virtual particle; ②The real particle neighbors of each virtual particle; ③ The virtual particle neighbors of each real particle; ④ The real particle neighbors of each real particle.

5. The tensile-stabilized SPH fluid simulation method according to claim 2, characterized in that: velocity at the virtual particle position The calculation formula is: In the above formula, This refers to the kernel function in the SPH method.

6. The tensile-stabilized SPH fluid simulation method according to claim 2, characterized in that: To avoid the accumulation of density errors during the simulation, the source term, i.e., the right-hand side term, of the pressure Poisson equation is compensated using the following formula: in, The source terms of the Poisson equation before compensation, These are the compensated source terms of the Poisson equation. The density value calculated using the real particle neighborhood of the virtual particle is: , real particles quality This is an empirical coefficient, with a value of 1.

Images