CN109726496B - IISPH-based implementation method for improving simulation efficiency of incompressible water - Google Patents

IISPH-based implementation method for improving simulation efficiency of incompressible water Download PDF

Info

Publication number
CN109726496B
CN109726496B CN201910011938.6A CN201910011938A CN109726496B CN 109726496 B CN109726496 B CN 109726496B CN 201910011938 A CN201910011938 A CN 201910011938A CN 109726496 B CN109726496 B CN 109726496B
Authority
CN
China
Prior art keywords
particle
formula
density
pressure
divergence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910011938.6A
Other languages
Chinese (zh)
Other versions
CN109726496A (en
Inventor
艾明晶
李锋
郑爱玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201910011938.6A priority Critical patent/CN109726496B/en
Publication of CN109726496A publication Critical patent/CN109726496A/en
Application granted granted Critical
Publication of CN109726496B publication Critical patent/CN109726496B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention relates to a realization method for improving simulation efficiency of an incompressible water body based on an IISPH (Implicit incompressible Smooth particle hydrodynamics) method. Firstly, introducing a speed non-divergence model into IISPH, and correcting divergence errors of a speed field in each time step of the IISPH method; secondly, the density constant model and the speed non-divergence model share a calculation factor through the improved density constant model, and efficient simulation of the incompressible water body is achieved. By introducing the non-divergence speed model, the divergence error is corrected, and the problem that the density error is continuously increased along with time is solved; by improving the density constant model, the calculation factors of the speed non-divergence model and the density constant model are shared, and the problem of redundant calculation is solved. The method can achieve the real-time simulation speed of 30.10 frames/second under the 17k particle scale, and compared with a PCISPH (Predictive-Corrective Incompression SPH) method, the simulation efficiency is improved by 24.17%; compared with the IISPH method, the simulation efficiency is improved by 16.21%.

Description

IISPH-based implementation method for improving simulation efficiency of incompressible water
Technical Field
The invention relates to the field of computer graphics, in particular to a realization method for improving simulation efficiency of incompressible water based on IISPH (Implicit incompressible Smooth particle hydrodynamics).
Background
Limited by the computing performance and efficiency of computers, water simulation technology is in the theoretical research stage for a long time. In recent years, with the improvement of computer hardware performance and computational efficiency, water simulation is realized, and the development of computer graphics further promotes the development of a water simulation technology based on physics, so that the water simulation technology gradually becomes one of the most concerned research directions. The Smooth Particle Hydrodynamics (SPH) method is a meshless method for calculating the overall motion state of a fluid (including a water body) by solving a particle group kinetic equation, and is widely applied to various fields such as movies, games, advertisements, military affairs and the like because of the advantages of conservation of mass, easiness in capturing water bloom and bubbles and the like.
The traditional SPH method solves the pressure on water body particles by using an ideal gas state equation, so that the water body has obvious compressibility, and obvious unnatural oscillation occurs on the surface of the water body when large-scale water body and large-scale movement are simulated, and obvious visual distortion is caused, so that the realization of the incompressibility of the water body has important significance for improving the simulation reality of the water body.
The existing thinking for realizing the incompressibility of the water body is mainly divided into two categories: explicit incompressible methods and implicit incompressible methods. The explicit incompressible method solves the Pressure by directly solving a liquid state Equation or a Pressure Poisson Equation (PPE), and the implicit incompressible Equation mostly solves the Pressure by means of predictive correction.
1. Explicit incompressible method
The explicit incompressible method mainly obtains the pressure applied to the water body particles by directly solving a liquid state equation or a pressure Poisson equation, thereby realizing the incompressibility of the water body.
In 2007, Becker et al proposed a slightly Compressible SPH (WCSPH) method (reference 1: Monaghan, J.J.Simuling free surface flows with SPH [ J ]. Journal of Computational Physics,1994,110(2): 399-. The WCSPH ensures the incompressibility of the water body simulation, and obviously improves the reality of the water body simulation. However, the stiffness of the water body dominates the CFL (Courant-Friedrichs-Levy) condition, the large stiffness coefficient only allows the use of a small time step, the calculation cost is increased along with the reduction of compressibility, and therefore the WCSPH method cannot effectively simulate complex water body movement.
In the same year, Adams et al proposed an Incompressible SPH (ISPH) method (reference 2: X.Y.Hu, N.A.Adams.an in compressible multi-phase SPH method [ J ]. Journal of Computational Physics,2007,227(1): 264-. According to the method, the pressure Poisson equation is directly solved through the conjugate gradient, so that the pressure applied to the particles is obtained. The ISPH method firstly calculates other forces except the pressure as resultant force, secondly calculates the intermediate speed of the particles through the resultant force, and then disperses the Poisson equation to calculate the pressure applied to the particles. Unlike WCSPH, ISPH allows the use of a large time step, but this method requires the solution of a complex formula and equation set in each time step, and the calculation overhead is large, and efficient simulation of an incompressible water body cannot be realized.
In 2015, Kang et al introduced the non-Divergence Condition into the ISPH method (ref. 3: Kang N, Sagong D. Incompressible SPH using the Divergence-Free Condition [ J ]. Computer Graphics Forum,2015,33(7):219-228) which implemented the non-Divergence Condition when solving the momentum equations. Ensuring that the volume is incompressible while maintaining a divergence-free velocity field, the incompressible being performed simultaneously on the position level and on the velocity level. Compared with ISPH, the method can greatly accelerate the convergence rate of pressure solving and reduce the iteration times of a solver, thereby improving the simulation efficiency of the water body. However, in the simulation process, the phenomenon that the particle density at the free surface of the water body is underestimated, so that particle aggregation and the like are generated easily occurs.
2. Implicit incompressible method
The implicit incompressible method mainly solves the pressure Poisson equation indirectly to obtain the pressure on the water body particles, thereby realizing the incompressibility of the water body.
In 2009, solenthler et al proposed a Predictive-Corrective based SPH (PCISPH) method (ref.4: solenthler B, Pajarola r. Predictive-Corrective SPH [ J ]. Acm Transactions on Graphics,2009,28(3):1-6) which is divided into two phases: a prediction stage and a correction stage. In the prediction phase, the method uses other forces than the pressure as a resultant force, and calculates the properties of the particle, such as the speed, the position and the like, of the next time step, namely the predicted physical properties of the particle according to the resultant force. An iterative correction loop is then performed, adjusting the velocity of the particles by applying pressure to the particles when the error in the particle density is greater than a given error value, and exiting the correction loop if the average density error is within the given error range. And finally, recalculating the particle attribute of the next time step by using the solved pressure and other forces. The method avoids the calculation cost of directly solving the pressure Poisson equation, but the indirect pressure solving needs to be carried out through multiple iterations, and the calculation cost is still increased to a certain extent.
In 2012, He et al proposed a method for solving Local Poisson SPH (LPSPH) of an incompressible water body (reference 5: X.He, N.Liu, S.Li, et al, Local Poisson SPH for a visual in-compliance fluid [ J ]. Computer Graphics Forum,2012,31(6): 1948-1958). Firstly, converting a Poisson equation into an integral form; and then, a discretization method is applied to convert the continuous integral equation into discrete summation, so that the calculation cost of directly solving the global pressure Poisson equation is avoided.
In 2014, Ihmsen et al proposed a pressure projection formula, namely an Implicit Incompressible SPH (Impulse incorporated SPH, IISPH) method (refer to document 6: Ihmsen M, Cornelis J, Solenthler B, et al. Impulse incorporated compressed SPH [ J ]. IEEE Transactions on Visualization and Computer Graphics,2014,20(3):426-435.), which combines the symmetric SPH pressure of the continuity equation with the SPH dispersion to obtain a discrete form of the pressure Poisson equation, and solves the pressure using the relaxed Jacobi algorithm. This method improves the convergence speed of the solver compared to previous projection schemes, and furthermore, corrects the density deviation based on speed rather than position, which improves the stability of the time integration to some extent. However, the velocity field calculated by the method has divergence errors, so that the average iteration number of density correction is increased, and the simulation efficiency is reduced.
In 2014, Ihmsen et al proposed IISPH-LFIP (Implicit Incompression SPH-Fluid additive Particle, IISPH-FLIP) method (reference 7: Cornelis J, Ihmsen M, Peer A, et al. IISPH-FLIP for compressing fluids [ J ]. Computer Graphics Forum,2014,33(2): 255-262.). The method adopts an implicit incompressible method to perform pressure projection, and adopts a Fluid implicit Particle (FLIP) solver to perform boundary processing, so as to realize the simulation of the incompressible water body. This new combination solves two problems of the existing SPH and FLIP solver, namely the conservation of mass in the FLIP and the efficiency in the SPH. The IISPH-FLIP solver provided by the method can simulate the incompressible water body, and the quantized and imperceptible density deviation is less than 0.1%.
The explicit incompressible water simulation method directly solves a liquid state equation or a pressure Poisson equation, and the calculation complexity is high; implicit incompressible water simulation methods such as PCISPH and IISPH adopt a prediction correction method, but the number of iterations in the solution process is too large, and each iteration needs to perform a time-consuming adjacent particle search step, so that the calculation amount is increased sharply along with the increase of the number of particles. In summary, although the existing incompressible water simulation method realizes the incompressibility of the water body, the cost is too large when the incompressibility is applied, and the simulation efficiency is reduced.
Disclosure of Invention
The invention provides a realization method for improving the simulation efficiency of an incompressible water body based on IISPH, which comprises the steps of firstly, introducing a speed non-divergence model into an IISPH method, correcting the divergence error of the speed in the IISPH method, and avoiding the density error from increasing continuously along with the increase of time, thereby reducing the average iteration times of correcting the density error; and then improving the density constant model of the IISPH, and solving the density variation through the density variation rate, so that the improved density constant model can use the calculation factor solved by the speed non-divergence model, thereby avoiding redundant calculation and improving the simulation efficiency of the incompressible water body. The method specifically comprises the following steps:
the method comprises the following steps: introducing a speed non-divergence model to correct divergence errors of the IISPH method
The Navier-Stokes (N-S) equation is a kinematic equation used to describe the conservation of momentum of viscous incompressible fluids. This equation is primarily used to describe the motion of liquids and gases, and essentially follows newton's second law, which can also be expressed as: the momentum change of an infinitely small volume of fluid is the sum of its effects of gravity, viscous forces, pressure, and other forces. The N-S equation describes the law of fluid motion in nature, which is based on the nature of physics making it well suited for describing natural phenomena such as bodies of water, smoke, air flow, etc.
The three-dimensional expression form of the N-S equation is shown as the formula (1) and the formula (2).
Figure BDA0001937652730000031
Figure BDA0001937652730000032
Where u represents the velocity of the particles, ρ represents the density of the particles, p represents the pressure of the particles, i.e., the pressure exerted per unit area, and f represents the external force (e.g., gravity, etc.) to which the fluid is subjected. The operator · is a dot product of the spatial vector,
Figure BDA00019376527300000415
representing partial derivatives of spatial vectors over spatial dimensions, i.e.
Figure BDA0001937652730000041
Also known as gradient, operator
Figure BDA00019376527300000416
Is the laplacian operator.
Equation (1) describes the conservation of momentum when the fluid is in motion, and the right equation represents all the force field quantities of the fluid. The first term on the right of the equation represents the pressure difference between the particles, the second term represents the viscous force between the particles, and the third term is the external force (e.g., gravity). The physical significance is as follows: the movement of water body particles is determined by pressure, viscous force and external force.
The formula (2) is called as an incompressible condition, and the formula indicates that the velocity divergence of the water body in the moving process is 0, namely the inflow and outflow of the water body in a unit volume does not change in a certain time step, so that the incompressible property of the water body in the moving process is ensured.
The IISPH corrected particle velocity has divergence error, i.e. is not satisfied
Figure BDA0001937652730000042
The invention introduces a velocity non-divergence model and uses the model to calculate a pressure term, thereby realizing the correction of divergence errors of particle velocities, and the magnitude of the pressure borne by a particle i is determined by an equation (3).
Figure BDA0001937652730000043
Wherein
Figure BDA0001937652730000044
Which is indicative of the pressure to which the particle i is subjected,
Figure BDA0001937652730000045
which represents the volume of the water body particles,
Figure BDA0001937652730000046
indicating a pressure gradient. Introduction of
Figure BDA0001937652730000047
The calculation formula of the pressure gradient, which represents the stiffness coefficient for correcting the velocity v of the particle i, is shown in the formula (4).
Figure BDA0001937652730000048
WhereinmjRepresenting the mass of the neighboring particle j,
Figure BDA0001937652730000049
the gradient of the kernel function is represented by,
Figure BDA00019376527300000410
i.e. the stiffness coefficient meeting the condition of no divergence in speed.
The particles satisfy conservation of momentum during movement. The stress analysis of the particle i shows that the pressure applied to the particle i is the sum of the pressures of the adjacent particles j, and the pressure between the particles is a pair of interaction forces. Because the interaction forces are equal in magnitude and opposite in direction, the pressure exerted on the particle i and the pressure exerted on the adjacent particle j by the particle i satisfy
Figure BDA00019376527300000411
Wherein
Figure BDA00019376527300000412
The pressure of the particle i against the particle j is expressed, which means that the pressure exists as an inter-particle force in the form of an internal force, and the sum of the pressures of all the particles is 0. From the whole view of the water body, the pressure does not generate the change of the kinetic energy of the water body, and the law of conservation of the kinetic energy is satisfied. The pressure of the particle i to the adjacent particles is related to the distance between the particles, and the solving formula is shown in formula (5).
Figure BDA00019376527300000413
Wherein xjRepresenting the three-dimensional coordinates of the particle j, the solver must correct the velocity of the particle by changing the pressure value so that the velocity of the particle satisfies the non-divergence. The rate of change of the density with respect to time is shown in the formula (6).
Figure BDA00019376527300000414
Wherein WijRepresenting a smooth kernel function, and ensuring that the water body is not compressible, namely the density does not change, the change rate of the density obtained by the formula (6) relative to the time is 0. V in formula (6)i-vjWhich represents the relative velocity difference between particles caused by the pressure difference, the pressure to which different particles are subjected differs, thus causing the velocity difference between particles. The amount of speed change caused by the pressure is as shown in equation (7).
Figure BDA0001937652730000051
Wherein Δ viExpression (8) is obtained by substituting expression (7) for expression (6) indicating the amount of speed change due to pressure.
Figure BDA0001937652730000052
The pressure on the particles can be calculated by the formula (3), and the rigidity coefficient can be solved by combining the formula (6) with the density change rate of 0; by correcting the particle velocity using this stiffness coefficient, the divergence of the velocity can be made 0, and the rate of change of the density with respect to time can be made 0. Substituting formula (4) for formula (3) to obtain
Figure BDA0001937652730000053
Then, the compound is substituted with the formula (8) to obtain the formula (9).
Figure BDA0001937652730000054
Solving the rigidity coefficient according to the formula (9)
Figure BDA0001937652730000055
As shown in equation (10).
Figure BDA0001937652730000056
Wherein
Figure BDA0001937652730000057
And alpha isiIs only related to the particle position.
Equation (10) is a velocity non-divergence model that solves for a stiffness coefficient that allows velocity to satisfy non-divergence
Figure BDA0001937652730000058
The rigidity coefficient is used for calculating pressure, and the density change rate of the water body is guaranteed not to fluctuate in the movement process, so that the incompressibility of the water body is achieved.
Finally, the resultant force to which the particle is subjected is calculated according to equation (11), and the velocity of the particle is recalculated using the resultant force.
Figure BDA0001937652730000059
Wherein
Figure BDA00019376527300000510
Which represents the resultant force to which the particle i is subjected,
Figure BDA00019376527300000511
representing the resultant of forces other than pressure,
Figure BDA00019376527300000512
representing the pressure to which the particle i is subjected to the neighboring particle j.
And step two, improving the density constant model to enable the density constant model and the speed non-divergence model to share a calculation factor so as to correct the density error of the IISPH.
In the velocity non-divergence model, the particle velocity is corrected using a rigid coefficient satisfying the velocity non-divergence to ensure that the density change rate is less than a given divergence error value etadiv. The alpha factor needs to be calculated when the non-divergence velocity correction is carried out, and the density constant model provided by the invention reuses the alpha calculation factor to further correct the density of the water body particles.
Pass density of the inventionCalculating density variation according to the variation rate and the time step, and obtaining the density value rho 'at the moment of t + delta t'iAs shown in formula (12).
Figure BDA0001937652730000061
Substituting formula (8) for formula (12) by the density change rate of the velocity non-divergence model, and combining with the constant density condition rhoi=ρ0Obtaining the density error rho 'in delta t time step'i0As shown in formula (13).
Figure BDA0001937652730000062
Where ρ is0The initial density of the particles is represented by the formula (3) or (4) in place of the formula (13).
Figure BDA0001937652730000063
Figure BDA0001937652730000064
The rigidity factor is obtained from the formula (14)
Figure BDA0001937652730000065
The solving formula of (2) is shown in formula (15).
Figure BDA0001937652730000066
Wherein
Figure BDA0001937652730000067
In order to meet the requirement of a rigid coefficient with constant water density, alpha is a calculation factor in the velocity non-divergence model, and the calculation factor is only calculated once in the same time step, so that the time overhead caused by repeated calculation is reduced.
Equation (15) is a constant density model that solves for a constant density stiffness coefficient
Figure BDA0001937652730000068
The rigidity coefficient is used for calculating the pressure on the particles, and the density variation of the water body is ensured not to fluctuate in the moving process, so that the incompressibility of the water body is realized.
From the equation (11), the resultant force exerted on the particle after the density-constant model is corrected is shown in the equation (16).
Figure BDA0001937652730000069
The velocity of the particle is corrected by using the resultant force applied to the corrected particle, and the particle velocity v ' at the time t + delta t after correction can be obtained by combining Newton's second law 'iAs shown in (17).
Figure BDA00019376527300000610
Calculated according to the above
Figure BDA00019376527300000611
And the velocity v of the particle of formula (17)iCorrected to obtain a corrected particle velocity v'i
The method combines the speed non-divergence model formula (11) with the improved density constant model formula (15), and realizes the correction of the speed divergence error and the density error, thereby reducing the density correction iteration times, reducing the redundant calculation and effectively improving the simulation efficiency of the incompressible water body. Test results show that the simulation efficiency of the method provided by the invention is superior to that of IISPH and PCISPH methods when the incompressible water body is simulated.
Drawings
FIG. 1 is a flow chart of an implementation method for improving simulation efficiency of incompressible water based on IISPH.
FIG. 2 is a diagram of an improved density invariant model algorithm.
FIG. 3 is a velocity non-divergence model algorithm diagram.
Fig. 4 is a design of an incompressible water system module.
FIG. 5 is a graph comparing simulation efficiencies of PCISPH, IISPH, and the method of the present invention.
FIG. 6 is a comparison of multi-thread efficiency for the method of the present invention.
Fig. 7 is a three-dimensional scene water drop falling simulation and effect diagram.
Detailed Description
The technical solution and the result of the present invention will be further described in detail with reference to the accompanying drawings and examples.
The invention provides a speed non-divergence incompressible water body simulation method, which comprises the steps of firstly, introducing a speed non-divergence model into an IISPH method, correcting a divergence error of a speed in the IISPH method, and avoiding that a density error is continuously increased along with the increase of time, so that the average iteration times of correcting the density error is reduced; and then improving the density constant model of the IISPH, and solving the density variation through the density variation rate, so that the improved density constant model can use the calculation factor solved by the speed non-divergence model, thereby avoiding redundant calculation and remarkably improving the simulation efficiency of the incompressible water model.
The method comprises the following steps: and (3) introducing a speed non-divergence model into the IISPH, improving an IISPH density constant model, and constructing a non-compressible water body simulation framework, as shown in figure 1.
Step 1.1, performing adjacent particle search to obtain an adjacent particle set N of each water body particle in the current time stepi(t)。
Step 1.2, set N of neighboring particles according to each particlei(t) calculating a resultant force other than the pressure to which the particles are subjected, and calculating the predicted velocity and the predicted density of the particles based on the resultant force.
Step 1.3, combining the density constant model to correct the particle speed, firstly calculating a shared calculation factor, and calculating a rigidity coefficient according to the calculation factor, wherein the rigidity coefficient can enable the density error to be smaller than a given valueA density error value of (a); then, the velocity of the particles is corrected by using the rigidity coefficient so that the variation | ρ'i(t+Δt)-ρ0| less than a given density error value ηρ
Step 1.4, calculating a rigidity coefficient by using the shared calculation factor calculated in the step (1.3) in combination with the particle velocity obtained in the velocity non-divergence model correction step (1.3), wherein the rigidity coefficient can enable the density change rate error to be smaller than a given divergence error value etadivThen, the rigidity coefficient is used to correct the velocity of the particles so that the density change rate is less than a given divergence error value etadiv
And step 1.5, traversing all the particles, obtaining the particle speed according to the calculation, updating the speed to the particle speed of the next time step, and updating other attributes of the particles.
Step two: the density error was corrected in conjunction with a modified constant density model, as shown in FIG. 2.
Step 2.1, traverse all particles, combine its adjacent particle set Ni(t), calculating the average density variation.
Step 2.2, if the average density variation is less than the given density error value (p)avg0ρ) Or the iteration number iter is not less than 2, the density constant model is exited, otherwise, step 2.3 is executed.
Step 2.3, for each particle, traverse its set of neighboring particles Ni(t) calculating the corresponding stiffness coefficient
Figure BDA0001937652730000071
And calculating the pressure to which the particle and the adjacent particles are subjected according to the rigidity coefficient.
Step 2.4, correcting the particle velocity according to the pressure calculated in step 2.3, and using the corrected particle velocity, executing steps 2.1-2.4 in a circulating manner until the density variation is smaller than the given density error (rho)avg0ρ)。
Step three: the divergence error is corrected in conjunction with the velocity non-divergence model, as shown in FIG. 3.
Step 3.1, traverse all particles, combine its neighboring particle set Ni(t), the density change rate of the particles is calculated.
Step 3.2, if the density change rate is less than or equal to the given divergence error value
Figure BDA0001937652730000081
Or the iteration number is not less than 1, the speed non-divergence model is exited, otherwise, the step 3.3 is executed.
Step 3.3, for each particle, traverse its set of neighboring particles Ni(t) calculating the corresponding stiffness coefficient
Figure BDA0001937652730000082
And calculating the pressure to which the particle and the adjacent particles are subjected according to the rigidity coefficient.
Step 3.4, correcting the particle velocity according to the pressure calculated in step 3.3, and using the corrected particle velocity, executing steps 3.1-3.3 in a circulating manner until the velocity divergence error is smaller than the given divergence error value
Figure BDA0001937652730000083
Step four: and (4) constructing an incompressible water body motion scene by using the physical models solved in the second step and the third step, and outputting and displaying a result.
The incompressible water body model solved in the second step and the third step is only a result obtained based on mathematical physics solution, is visually described and is composed of a series of particles, and the actual water body scene cannot be vividly presented. Therefore, the model also needs to be correspondingly rendered, so that a water body with reality is obtained.
All tests of the invention were carried out on a Dall desktop, with the CPU model Intel Core i7-4790@3.6GHz and the graphics card model AMD RadeonTMR5 240。
Fig. 4 is a design diagram of a module of an incompressible water body system, and the whole system is divided into 4 modules: the device comprises an initialization module, a prediction correction module, a boundary processing module and a rendering module. The initialization module mainly receives parameters set by a user and generates particle initial information including attributes such as particle position, pressure, density, speed and the like.
The prediction correction module integrates a speed non-divergence model and a density constant model, is executed once in each time step, generates the particle speed of the next time step, and is mainly divided into two parts: a preprocessing section and a velocity calculation section. The preprocessing part mainly calculates and stores data only related to the positions of the particles, the data of the preprocessing part does not change in the time step, and repeated calculation can be avoided through a pre-calculation mode. The velocity calculating section corrects the predicted velocity of the particle using a velocity non-divergence model and a density constant model so that the corrected velocity satisfies velocity non-divergence and density constant.
The boundary processing module is used for correcting the speed corrected by the prediction correction module, so that the particles are prevented from penetrating the solid boundary in the moving process. Firstly, detecting whether particles collide with a boundary or not, and if so, correcting the particle speed through a collision response model based on repulsive force so as to ensure the impermeability of the container boundary; if no collision occurs, the collision response model does not need to be executed.
The rendering module is mainly used for realizing water body rendering according to the particle position information, and comprises the steps of anisotropic nuclear surface construction, depth smoothing, illumination application, sky box drawing and the like, so that the smoothness of the water body surface and the reality of the water body rendering are improved.
The invention introduces the speed non-divergence model into the IISPH, improves the density constant model of the IISPH, and can more efficiently simulate the incompressible water body motion scene by the implementation method for improving the simulation efficiency of the incompressible water body based on the IISPH.
The PCISPH algorithm, the IISPH algorithm and the method are respectively used for realizing the simulation of a water drop falling model, the model consists of 17074 particles, the test comprises 12 independent experiments, and the average frame rate acceleration ratio of 12 results is obtained through calculation. Table 1 is a comparison table of simulated frame rates of drop models, and fig. 5 is a corresponding line graph. The data in table 1 are calculated, and the average frame rate of the PCISPH algorithm is 24.24fps, the average frame rate of the IISPH algorithm is 25.91fps, and the average frame rate of the incompressible water acceleration simulation method is 30.10 fps. Compared with the PCISPH algorithm, the computational efficiency of the method is improved by 24.17%, and compared with the IISPH algorithm, the computational efficiency of the method is improved by 16.21%.
TABLE 1 Water drop model simulation frame rate comparison table
Figure BDA0001937652730000091
The invention uses multithreading to improve the efficiency of searching adjacent particles, kernel functions and gradient solving. In order to verify the effectiveness and efficiency of the multi-thread implementation of the program of the present invention, 1, 4, and 8 threads are respectively started to simulate 7104 particles, and a simulated frame rate diagram as shown in fig. 6 is obtained. As can be seen from fig. 6, when 1 thread is started, the real-time simulation requirement cannot be met; when 4 threads are started, the simulation efficiency is obviously improved, and the real-time requirement is met; when 8 threads are started, the simulation efficiency is further improved, so that the effectiveness and the high efficiency of the multithreading method are verified.
Fig. 7 depicts a three-dimensional scene water drop falling simulation effect diagram, fig. 7(a) shows an initial state of a water drop, fig. 7(b) and fig. 7(c) show the effects of multiple rebounds after the water drop falls and collides with the ground, and fig. 7(d) shows the effects of gradual restoration and stabilization of a water body, and it can be seen from the diagram that turbulence phenomena such as particle penetration boundary, particle adhesion boundary, particle aggregation at the boundary and the like do not occur in the moving process of the water body.

Claims (1)

1. An implementation method for improving simulation efficiency of incompressible water based on IISPH is characterized by comprising the following implementation steps:
the method comprises the following steps: introducing a speed non-divergence model to correct divergence errors of the IISPH method:
the IISPH corrected particle velocity has divergence error, i.e. is not satisfied
Figure FDA0002650837700000015
By introducing a speed non-divergence model and using the model to calculate a pressure term, the divergence error of the particle speed is corrected, and the pressure on the particle i is determined by the formula (1);
Figure FDA0002650837700000011
wherein Fi pWhich is indicative of the pressure to which the particle i is subjected,
Figure FDA0002650837700000016
which represents the volume of the water body particles,
Figure FDA0002650837700000017
representing a pressure gradient; introduction of
Figure FDA0002650837700000018
A calculation formula of a pressure gradient, which is a rigidity coefficient for correcting the velocity v of the particle i, is represented by the following formula (2);
Figure FDA0002650837700000012
wherein m isjRepresenting the mass of the neighboring particle j,
Figure FDA0002650837700000019
the gradient of the kernel function is represented by,
Figure FDA00026508377000000110
namely the rigidity coefficient meeting the condition of no divergence of the speed;
the particles satisfy the conservation of momentum in the process of moving, and the stress analysis on the particles i can know that the pressure borne by the particles i is the sum of the pressures of the adjacent particles j, and the pressure among the particles is a pair of interaction forces due to the phaseThe interaction forces are equal in magnitude and opposite in direction, so that the pressure on the particle i and the pressure on the adjacent particle j by the particle i satisfy
Figure FDA00026508377000000111
Wherein
Figure FDA00026508377000000112
The pressure of the particle i to the particle j is expressed, which means that the pressure exists as an inter-particle acting force in the form of an internal force, and the sum of the pressures of all the particles is 0; from the whole view of the water body, the pressure does not generate the change of the kinetic energy of the water body, and the law of conservation of the kinetic energy is satisfied; the pressure of the particle i to the adjacent particles is related to the distance between the particles, and the solving formula is shown in formula (3),
Figure FDA0002650837700000013
wherein xjRepresenting the three-dimensional coordinates of the particle j, the solver must correct the velocity of the particle by changing the pressure value so that the velocity of the particle satisfies the non-divergence property, the rate of change of the density with respect to time is as shown in equation (4),
Figure FDA0002650837700000014
wherein WijRepresenting a smooth kernel function, and ensuring that the water body is not compressible, namely the density is not changed, wherein the change rate of the density obtained by the formula (4) relative to time is 0; v in the formula (4)i-vjRepresenting the relative velocity difference among the particles, wherein the velocity difference is caused by the pressure difference, and the pressure of different particles is different, so that the velocity of different particles is different; the amount of speed change caused by the pressure is as shown in equation (5);
Figure FDA0002650837700000021
wherein Δ viRepresenting the amount of speed change due to pressure, and substituting formula (5) for formula (4) to obtain formula (6);
Figure FDA0002650837700000022
the pressure on the particles can be calculated by the formula (1), and the rigidity coefficient can be solved by combining the formula (4) with the density change rate of 0; by correcting the particle velocity using the rigidity coefficient, the divergence of the velocity can be made 0, and the rate of change of the density with respect to time can be made 0; substituting formula (2) into formula (1) to obtain
Figure FDA0002650837700000023
Then substituting the formula (6) to obtain a formula (7);
Figure FDA0002650837700000024
solving the rigidity coefficient according to the formula (7)
Figure FDA0002650837700000027
As shown in formula (8);
Figure FDA0002650837700000025
wherein
Figure FDA0002650837700000028
And alpha isiThe size of (a) is related only to the particle position;
equation (8) is a velocity non-divergence model that solves for stiffness coefficients
Figure FDA0002650837700000029
And using the stiffness coefficient to calculate the pressure, therebyEnsuring that the density change rate of the water body does not fluctuate in the movement process, namely ensuring the incompressibility of the water body;
finally, calculating the resultant force applied to the particle according to the formula (9), and recalculating the velocity of the particle by using the resultant force;
Figure FDA0002650837700000026
wherein
Figure FDA00026508377000000210
Which represents the resultant force to which the particle i is subjected,
Figure FDA00026508377000000211
representing the resultant of forces other than pressure,
Figure FDA00026508377000000212
representing the pressure to which particle i is subjected to adjacent particle j;
step two, improving the density constant model to enable the density constant model and the speed non-divergence model to share a calculation factor so as to correct the density error of IISPH:
in the velocity non-divergence model, the particle velocity is corrected using a rigid coefficient satisfying the velocity non-divergence to ensure that the density change rate is less than a given divergence error value etadiv(ii) a Alpha factors need to be calculated when non-divergence speed correction is carried out, the proposed density constant model multiplexes the alpha calculation factors, and the density of water body particles is further corrected;
calculating the density variation through the density variation rate and the time step, and then calculating the density value rho 'at the t + delta t moment'iAs shown in formula (10);
Figure FDA0002650837700000031
the formula (6) is replaced by the formula (10) by solving the density change rate of the speed non-divergence model, and the density constant condition (rho) is combinedi=ρ0) Obtaining the density error rho 'in delta t time step'i0As shown in formula (11);
Figure FDA0002650837700000032
where ρ is0The initial density of the particles is expressed by substituting the formulas (1) and (2) into the formula (11);
Figure FDA0002650837700000033
Figure FDA0002650837700000034
the rigidity factor is obtained from the formula (12)
Figure FDA0002650837700000036
As shown in formula (13);
Figure FDA0002650837700000035
wherein
Figure FDA0002650837700000037
In order to meet the requirement of a rigidity coefficient with constant water density, alpha is a calculation factor in a speed non-divergence model, and the calculation factor is only calculated once in the same time step, so that the time overhead caused by repeated calculation can be reduced;
equation (15) is a density-invariant model obtained by solving the stiffness coefficient
Figure FDA0002650837700000038
And the pressure calculated by the rigidity coefficient is used for ensuring that the density variation of the water body does not fluctuate in the movement process, namely ensuring the water bodyIncompressibility;
from the formula (11), the resultant force exerted on the particle after the density-constant model is corrected is shown as the formula (14);
Figure FDA0002650837700000041
the velocity of the particle is corrected by using the resultant force applied to the corrected particle, and the particle velocity v ' at the time t + delta t after correction can be obtained by combining Newton's second law 'iAs shown in (15);
Figure FDA0002650837700000042
calculated according to the above
Figure FDA0002650837700000043
And the velocity v of the particle of formula (15)iCorrected to obtain a corrected particle velocity v'i
Step three: and (3) constructing an incompressible water body motion scene by using the physical models solved in the first step and the second step, outputting and displaying the result:
and realizing water body rendering by adopting a screen space rendering technology, thereby displaying the fluid effect of the incompressible water body motion scene.
CN201910011938.6A 2019-01-07 2019-01-07 IISPH-based implementation method for improving simulation efficiency of incompressible water Active CN109726496B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910011938.6A CN109726496B (en) 2019-01-07 2019-01-07 IISPH-based implementation method for improving simulation efficiency of incompressible water

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910011938.6A CN109726496B (en) 2019-01-07 2019-01-07 IISPH-based implementation method for improving simulation efficiency of incompressible water

Publications (2)

Publication Number Publication Date
CN109726496A CN109726496A (en) 2019-05-07
CN109726496B true CN109726496B (en) 2020-12-11

Family

ID=66298146

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910011938.6A Active CN109726496B (en) 2019-01-07 2019-01-07 IISPH-based implementation method for improving simulation efficiency of incompressible water

Country Status (1)

Country Link
CN (1) CN109726496B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112560326B (en) * 2019-09-26 2023-07-18 腾讯科技(深圳)有限公司 Method and device for determining pressure field
CN115906596B (en) * 2022-11-18 2024-03-22 上海索辰信息科技股份有限公司 Wall oil film calculation method

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102867094A (en) * 2012-09-19 2013-01-09 西安交通大学 Establishment method for free surface flow model in moving particle semi-implicit algorithm
CN105956262A (en) * 2016-04-28 2016-09-21 清华大学 Multi-component solid and fluid simulation method and system based on SPH (Smoothed Particle Hydrodynamics) method
JP2017059444A (en) * 2015-09-17 2017-03-23 トヨタ自動車株式会社 Electrode simulation method and device for all-solid battery, and method of manufacturing electrode for all-solid battery
CN107908918A (en) * 2017-10-19 2018-04-13 新疆大学 The SPH method for numerical simulation of grains of sand surge start in a kind of flat sand bed
CN108269299A (en) * 2017-01-04 2018-07-10 北京航空航天大学 A kind of viscous fluid modeling method based on SPH method approximate solutions

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150242545A1 (en) * 2014-02-21 2015-08-27 Junghyun Cho Method of Simulation of Moving Interfaces using Geometry-Aware Volume of Fluid Method
CN106777662B (en) * 2016-12-12 2020-04-10 西安交通大学 Aircraft fuel tank oil mixing characteristic optimization method based on smooth particle fluid dynamics
US20180239848A1 (en) * 2017-02-21 2018-08-23 Livermore Software Technology Corporation Numerical Blast Simulation Methods and Systems Thereof

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102867094A (en) * 2012-09-19 2013-01-09 西安交通大学 Establishment method for free surface flow model in moving particle semi-implicit algorithm
JP2017059444A (en) * 2015-09-17 2017-03-23 トヨタ自動車株式会社 Electrode simulation method and device for all-solid battery, and method of manufacturing electrode for all-solid battery
CN105956262A (en) * 2016-04-28 2016-09-21 清华大学 Multi-component solid and fluid simulation method and system based on SPH (Smoothed Particle Hydrodynamics) method
CN108269299A (en) * 2017-01-04 2018-07-10 北京航空航天大学 A kind of viscous fluid modeling method based on SPH method approximate solutions
CN107908918A (en) * 2017-10-19 2018-04-13 新疆大学 The SPH method for numerical simulation of grains of sand surge start in a kind of flat sand bed

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"An incompressible multi-phase SPH method";X.Y. Hu etc.;《Journal of Computational Physics》;20071110;第227卷(第1期);全文 *
"Numerical investigation of Newtonian and non-Newtonian multiphase flows using ISPH method";A. Zainali etc.;《Computer Methods in Applied Mechanics and Engineering》;20130228;第254卷;全文 *
"单轴压缩条件下岩石破坏的光滑粒子流体动力学数值模拟";周小平等;《岩石力学与工程学报》;20150531;第 34 卷(第增 1期);全文 *

Also Published As

Publication number Publication date
CN109726496A (en) 2019-05-07

Similar Documents

Publication Publication Date Title
Schneiders et al. An efficient conservative cut-cell method for rigid bodies interacting with viscous compressible flows
Yokoi A density-scaled continuum surface force model within a balanced force formulation
Chentanez et al. Real-time Eulerian water simulation using a restricted tall cell grid
Kim et al. Practical animation of turbulent splashing water
Mullen et al. Energy-preserving integrators for fluid animation
Zhang et al. A PPPM fast summation method for fluids and beyond
US8803887B2 (en) Computer graphic system and method for simulating hair
Ando et al. A particle-based method for preserving fluid sheets
US9984489B2 (en) Fluid dynamics framework for animated special effects
Wang et al. Hierarchical optimization time integration for cfl-rate mpm stepping
Jones et al. Deformation embedding for point-based elastoplastic simulation
Weiler et al. Projective fluids
CN109726496B (en) IISPH-based implementation method for improving simulation efficiency of incompressible water
Thürey et al. Interactive free surface fluids with the lattice Boltzmann method
CN112861445A (en) Simulation method of grid-free numerical model
EP3674962A1 (en) Computer simulation of physical fluids on a mesh in an arbitrary coordinate system
Liu et al. High-order particle method for solving incompressible Navier–Stokes equations within a mixed Lagrangian–Eulerian framework
Kou et al. Immersed boundary method for high-order flux reconstruction based on volume penalization
Winchenbach et al. Constrained neighbor lists for SPH-based fluid simulations.
Ibayashi et al. Simulating liquids on dynamically warping grids
US20090040219A1 (en) System and method for surfacing of particle systems
Liu et al. Turbulent details simulation for SPH fluids via vorticity refinement
Bolduc et al. A high-order entropically-damped artificial compressibility approach on moving and deforming domains
Weaver et al. Fluid Simulation by the Smoothed Particle Hydrodynamics Method: A Survey.
Nikitin et al. A splitting method for numerical simulation of free surface flows of incompressible fluids with surface tension

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant