Home About the Journal Latest Work Current Issue Archive Special Issues Editorial Board


2021,  3 (2):   118 - 128

Published Date:2021-4-20 DOI: 10.1016/j.vrih.2021.01.003


SPH method has been widely used in the simulation of water scenes. As a numerical method of partial differential equations, SPH can easily deal with the distorted and complex boundary. In addition, the implementation of SPH is relatively simple, and the results are stable and not easy to diverge. However, SPH method also has its own limitations. In order to further improve the performance of SPH method and expand its application scope, a series of key and difficult problems restricting the development of SPH need to be improved.
In this paper, we introduce the idea of Helmholtz decomposition into the framework of smoothed particle hydrodynamics (SPH) and propose a novel velocity projection scheme for three-dimensional water simulation. First, we apply Helmholtz decomposition to a three-dimensional velocity field and decompose it into three orthogonal subspaces. Then, our method combines the idea of spatial derivatives in SPH to obtain a discrete Poisson velocity equation. Finally, the conjugate gradient (CG) is utilized to efficiently solve the Poisson equation.
The experimental results show that the proposed scheme is suitable for various situations and has higher efficiency than the current SPH projection scheme.
Compared with the previous projection scheme, our solution does not need to modify the particle velocity indirectly by pressure projection, but directly by velocity field projection. The new scheme can be well integrated into the existing SPH framework, and can be applied to the interaction of water with static and dynamic obstacles, even for viscous fluid.


1 Introduction
Since it was first proposed in 1977[1], the smoothed particle hydrodynamics (SPH) method has been widely used in the simulation of various types of fluids. The SPH method, which yields the numerical solutions of partial differential equations, can easily deal with distorted and complex boundaries. In addition, SPH is relatively simple and provides stable results that are difficult to diverge[2]. After decades of development, the SPH method has made remarkable achievements in the simulation of fluids[3,4], elastic solids[5], foam[6], snow[7], multiphase flow[8-10], and fluid control[11-14].
The standard SPH (SSPH)[15,16] uses an equation of state (EOS) to calculate pressure, which is suitable for compressible fluids. The weakly compressible SPH[17] imposes strict time-step restrictions on weakly compressible fluids, which limits its performance. Predictive-corrective incompressible SPH (PCISPH)[18] calculates the corresponding pressure value based on SSPH and EOS iterative prediction and correction of particle position. The PCISPH can deal with a larger time step and improve computational efficiency. It is based on the incompressible density of particles. The IISPH (Implicit Incompressible SPH) approach[19] allows a greater time step and a higher convergence rate, which further improves the computational efficiency, but its accuracy and stability remain to be improved.
In terms of details, SPH is not accurate in the calculation of fluid dynamic pressure, which leads to the low accuracy of SPH in simulating local phenomena of the fluid boundary layer and vortex. Regarding operation speed, for each calculation unit, SPH requires the interaction information of 30-50 surrounding particles, which is much larger than the information required by grid-based approaches. For example, the finite volume method[20] generally requires information on 5-13 adjacent elements. In addition, the additional search time of neighboring particles makes the fluid simulation of SPH on a small scale more time-consuming. In terms of mathematics, SPH lacks a complete error analysis system. Once the SPH simulation is started, each particle will be in a disordered but non-random state, which makes the error analysis of SPH very difficult. Therefore, improving the accuracy and efficiency of the SPH method remains a difficult problem.
In this paper, a new SPH projection scheme is proposed to simulate water. The idea of Helmholtz decomposition is applied to a three-dimensional velocity field to eliminate the velocity divergence introduced by the pressure term. Compared with the previous projection scheme, the proposed scheme does not need to modify the particle velocity indirectly through pressure projection, but directly projects the velocity field. In addition, because of the decoupling of the pressure and velocity fields, the conjugate gradient (CG) can be used to solve the Poisson equation efficiently. Various experimental results show that the new scheme proposed in this paper can be combined with the existing SPH framework, and is suitable for the interaction between water and static and dynamic obstacles, and can even be extended to viscous fluids. Concurrently, compared with the most advanced SPH projection schemes, this scheme has higher computational efficiency.
2 Background
Before introducing the proposed scheme, we first briefly explain the background knowledge and concepts related to our method in this section, including the SPH and the Helmholtz decomposition methods.
2.1 Smoothed Particle Hydrodynamics (SPH)
In 1999, the SPH method was first introduced into the field of fluid animation by Stam[21], while the initial application of SPH in the field of water simulation dates back to 2003[16]. The basic idea of SPH is to treat a continuous fluid as a collection of discrete fluid particles, and that the complex fluid motion is caused by the interaction between the fluid and rigid particles. The main process of SPH in water simulation involves the following steps:
(1) Establish the spatial index for all particles.
(2) Compute the pressure, density, viscous force, and external force for all particles, and then evaluate the acceleration of the particles based on the action of these forces.
(3) Compute the speed and position of the next moment according to the particle's acceleration and time step (time integral), then handle the necessary collisions.
(4) Render all particles.
2.1.1 SPH interpolation and smooth kernel function
Because SPH is based on discrete particles, it is necessary to interpolate the physical quantities of particles in a certain neighborhood for the calculation of a physical quantity at any position in the space:
is the physical quantity at an arbitrary position
is the known physical quantity at the position
of the neighboring particle, and
is a kernel function of the form
is called the smooth kernel radius; therefore,
is also called the smooth kernel function.
2.1.2 Spatial derivatives
Based on the differentiable smooth kernel function, the SPH method can avoid the problems of mesh winding and distortion in the grid-based method; thus, the spatial derivatives can be computed efficiently. The SPH method generally uses the approximation of the spatial derivatives:
2.1.3 Force analysis of particles
The force acting on the particle directly affects the movement of the particle, thereby affecting the dynamic behavior of the fluid as a whole; hence, an accurate analysis of the force on the particle is essential. The force acting on a particle generally consists of three parts, written as
, and
are the external force, pressure gradient force, and viscous force, respectively, and
is the viscosity coefficient.
The acceleration calculation of each particle follows Newton's second law, i.e.,
, and the acceleration of a particle can be obtained by substituting Equation (7) into Newton's second law:
2.2 Helmholtz decomposition
In 1858, Helmholtz proposed the decomposition theorem of flow field velocity (hereinafter referred to as the Helmholtz decomposition theorem)[22]. The Helmholtz decomposition theorem states that at an appropriate asymptotic infinity, a vector field can be decomposed into two parts:
The divi-free part represents the rotation, which can be expressed as the curl of the vector potential function.
The curl-free part represents translation and compression/expansion, which can be expressed as the gradient of a scalar potential function.
This means that when the divergence and curl of the vector field are known, the vector field can be uniquely constructed. Its inverse process, i.e., the process of decomposing the vector field into non-divergence and non-curvature components, is called Helmholtz decomposition, which can be described as
is a smooth vector field on the domain
with a smooth boundary
is a scalar potential function, and
is the vector potential function. By definition,
is curl-free (
) and
is div-free (
). Equation (9) can be abbreviated as
, where
is the divergent component and
is the rotational component that is tangent to the boundary
Chorin and Denaro et al., respectively, gave the boundary conditions of the Helmholtz decomposition and proved the existence, orthogonality, and uniqueness of the decomposition[23,24]. Specifically, the Helmholtz decomposition must satisfy one of the following boundary conditions:
where n is the outer normal vector of the boundary, and
are the normal and tangential components of the boundary, respectively.
The above derivation only considers the situation in infinite space, i.e., the situation where the harmonic field is zero. For bounded regions, there may be non-zero harmonic components. In this case, the Helmholtz decomposition will split the vector field into three orthogonal parts:
where h is the harmonic component, and its divergence and curl are both zero. In this case, the Helmholtz decomposition needs to satisfy the boundary conditions expressed by Equation (10).
3 Helmholtz decomposition-based projection scheme
This section proposes a projection scheme based on the Helmholtz decomposition of the three-dimensional velocity field. The advantage of this scheme lies in the decoupling of the pressure field and the velocity field when solving. Specifically, because the calculation guarantees the orthogonality of the decomposition, the calculation error of one item will not be reflected in the calculation of the other item. This means that it is more effective to use the Helmholtz decomposition of the velocity field for projection than to solve the Navier-Stokes equations in which velocity and pressure are coupled. The kernel process of our scheme is shown in Figure 1.
3.1 Particle velocity prediction
At the beginning of each time step, our scheme calculates the particle acceleration based on the external force, and updates the particle velocity with this acceleration. This process is called particle velocity prediction. First, the smooth kernel function is used to approximate the pressure and viscous force on the particle:
can be approximated as
is calculated according to the EOS:
, where parameters
are used to scale the pressure proportionally, called the stiffness coefficient, and
is the static density of the fluid.
Then, the current particle acceleration
can be evaluated by substituting the calculated pressure
and viscous force
into Equation (8). Finally, the particle velocity at each time step can be computed according to the particle acceleration:
is the size of the time step.
3.2 Helmholtz decomposition of the 3D velocity field
Owing to the existence of the pressure term, the velocity field calculated in the previous step contains divergence, which needs to be eliminated when handling incompressible flows. In this study, the projection operation is based on the principle of the Helmholtz decomposition. First, the velocity field
is decomposed into three parts in the form of Equation (11):
Since we assume that the fluid velocity field is irrotational, the following limitations exist:
Therefore, the rotational component
in Equation (15) is only a vector field with zero everywhere, and thus Equation (15) can be simplified as
To ensure incompressibility, the harmonic components
need to be retained to eliminate divergence:
By applying the divergence operator to both sides of Equation (16), after deformation, we can obtain
The term
is zero for the harmonic vector field. Therefore, Equation (17) can be simplified as a Poisson equation:
. After solving this Poisson equation, we can obtain the scalar field
, and then we substitute it into Equation (16), and the non-divergence velocity field
can be obtained.
3.3 Poisson equation solution based on the conjugate gradient
The Poisson equation can be solved in different ways, such as the Jacobi iteration[20], successive over-relaxation[25], and CG[26]. Among them, CG can solve Poisson equations with complex boundary conditions in an efficient interactive way, so it is an ideal solution.
For the simplified Poisson equation
, we use the smooth kernel function to approximate and obtain the discrete Poisson equation:
Equation (18) can be written as a linear equation system
. To solve this system, we first set the initial conditions
ϕ0=0, q0=r0=b
Then, we perform the following iterative calculations:
αl=rlTrlqlTAql, ϕl+1=ϕl+αlql, rl+1=rl-αlAql, ql+1=ql+rl+1Trl+1rlTrlql
The termination condition of the iteration is
, where
is the iteration termination threshold.
We can obtain
v=v* -φ
after iteration, and the divergence-free velocity field
can be calculated as
v=v* -ϕ
Thus far, the projection scheme based on the Helmholtz decomposition has been completely introduced.
4 Results
In this section, we discuss the experimental results of our scheme and compare the results of the different approaches. All experiments were performed on a personal computer equipped with a 2.6 GHz Intel i7-6700HQ CPU, 32 GB memory, and NVIDIA GeForce GTX 1060 GPU. We integrated our projection scheme into the SPH framework proposed by Koschier et al.[27]. For boundary handling, we used the scheme proposed by Akinci et al.[28], surface tension was treated using the method proposed by He et al.[29], and the liquid surface was reconstructed using the method proposed by Akinci et al.[30]. The images were rendered using Arnold.
For all the following scenes, we compared our scheme with IISPH[19] and DFSPH (Divergence-free SPH)[31]. The radius of the liquid particles in all scenes was 0.025 m, and the maximum volume error rate was set to 0.01%. In addition, the liquids in all the scenes were simulated in a cubic water tank. Considering the lighting factor, we did not render the water tank.
4.1 Two-dimensional scene
Figure 2 shows the rendering result in a two-dimensional scene, where the gray liquid is water and the blue liquid is viscous fluid, separated by obstacles. This scene contains 1726 liquid particles. The efficiency comparison between the proposed method and other approaches in this scene is shown in Table 1. For the given time steps, the proposed method is superior in average iteration times and simulation time (only slightly inferior to DFSPH when the time step is 5ms), and the advantage of our method is most evident when the time step is 2ms.
Comparison of the efficiency of each method under different time steps in a two-dimensional scene. For DFSPH, the avg. iter. column represents the number of iterations of the constant density solver and the divergence-free solver
Δt/msIISPHDFSPHOur scheme
avg. iter.comp. time/msavg. iter.comp. time/msavg.iter.comp. time/ms
4.2 Double emitter
Figure 3 shows the rendering result in a double-emitter scene, where the left emitter starts to drain from the 1st frame and stops at the 52nd frame, and then the right one starts to drain. This scene contains 1938 liquid particles at the beginning, and the number of particles increases with time. The efficiency comparison between the proposed method and other approaches in this scene is shown in Table 2. The efficiency of our method is higher than that of IISPH, but compared with DFSPH, our scheme only has advantages when the time step is 0.5ms.
Comparison of the efficiency of each method under different time steps in double emitter scene
Δt/msIISPHDFSPHOur scheme
avg. iter.comp. time/msavg. iter.comp. time/msavg. iter.comp. time/ms
4.3 Double dam break with sphere
Figure 4 shows the rendering result in a double-dam-break scene with a floating green sphere interacting with water. This scene contains 11492 liquid particles. The efficiency comparison between the proposed method and other approaches in this scene is shown in Table 3. The computational overhead of our method is lighter than that of IISPH and is similar to that of DFSPH.
Comparison of the efficiency of each method under different time steps in double dam break with sphere scene
Δt/msIISPHDFSPHOur scheme
avg. iter.comp. time/msavg. iter.comp. time/msavg. iter.comp. time/ms
4.4 Spinner
Figure 5 shows the rendering result in a spinner scene, where a red cuboid that continuously rotates clockwise exists at the bottom of the scene. This scene contains 5054 liquid particles. As shown in Table 4, compared to IISPH and DFSPH, our scheme requires fewer iterations and achieves higher efficiency at certain time steps.
Comparison of the efficiency of each method under different time steps in the spinner scene
Δt/msIISPHDFSPHOur scheme
avg. iter.comp. time/msavg. iter.comp. time/msavg. iter.comp. time/ms
4.5 Limitations
As shown in the above experiment, the projection scheme proposed in this chapter can not only be easily integrated into the existing SPH framework, but also performs well in two-dimensional scene simulations, fixed obstacle interaction scenes, and moving obstacle interaction scenes. In terms of computation overhead, our scheme achieves a noticeable improvement in average iterations and computation time compared to IISPH. In addition, the efficiency and performance of our scheme are also comparable to those of DFSPH.
However, our method has certain limitations. First, our method assumes that the three-dimensional velocity field is irrotational, so that the divergence-free component can be ignored when performing the Helmholtz decomposition. This means that our scheme cannot deal with fluid curl, so it is not suitable for simulating the liquid effects caused by vorticity. In addition, according to Ihmsen et al.[19], the application CG cannot clamp the negative pressure between iterations to eliminate the cohesion effect. In fact, clamping the negative pressure in the CG results in an extremely unstable pressure field. Finally, when using CG to solve the linear equation system, i.e.,
, the coefficient matrix
is required to be symmetric. This means that the CG is only suitable for simulating fluids with uniform density. Therefore, a solution to the Poisson equation that is suitable for fluids with non-uniform density is still needed in the future.
5 Conclusion and future work
This paper proposed a projection scheme based on the Helmholtz velocity field decomposition to simulate water. This paper first gave the Helmholtz decomposition form of the three-dimensional velocity field and then combined it with the SPH method to give the approximate discrete form of the velocity field Poisson equation. Finally, this paper used the CG to efficiently solve the discrete Poisson equation. Multiple sets of experiments have indicated that compared with the most advanced SPH projection schemes, our method can compute the harmonic velocity field with higher efficiency. The proposed scheme is also suitable for scenes interacting with various types of obstacles and can be extended to viscous fluid simulation.
For future research, we would like to improve the stability of our method to allow a larger time step. In addition, this study assumes that the velocity field is irrotational, so we plan to take curl into consideration. Finally, it would be interesting to promote our scheme to simulate more general phenomena, such as multiphase flow, fluid-rigid coupling, and even transformation between fluid and rigid.



Gingold R A, Monaghan J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181(3): 375–389 DOI:10.1093/mnras/181.3.375


Gissler C, Peer A, Band S, Bender J, Teschner M. Interlinked SPH pressure solvers for strong fluid-rigid coupling. ACM Transactions on Graphics, 2019, 38(1): 5 DOI:10.1145/3284980


Bao K, Zhang H, Zheng L L, Wu E H. Pressure corrected SPH for fluid animation. Computer Animation and Virtual Worlds, 2009, 20(2/3): 311–320 DOI:10.1002/cav.299


Huang C, Zhu J, Sun H Q, Wu E H. Parallel-optimizing SPH fluid simulation for realistic VR environments. Computer Animation and Virtual Worlds, 2015, 26(1): 43–54 DOI:10.1002/cav.1564


Peer A, Gissler C, Band S, Teschner M. An implicit SPH formulation for incompressible linearly elastic solids. Computer Graphics Forum, 2018, 37(6): 135–148 DOI:10.1111/cgf.13317


Bender J, Koschier D, Kugelstadt T, Weiler M. Turbulent micropolar SPH fluids with foam. IEEE Transactions on Visualization and Computer Graphics, 2018, 25(6): 2284–2295 DOI:10.1109/tvcg.2018.2805110


Gissler C, Henne A, Band S, Peer A, Teschner M. An implicit compressible SPH solver for snow simulation. ACM Transactions on Graphics, 2020, 39(4): 36 DOI:10.1145/3386569.3392431


Koschier D, Bender J, Solenthaler B, Teschner M. Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids. 2020


Ren B, Li C F, Yan X, Lin M C, Bonet J, Hu S M. Multiple-fluid SPH simulation using a mixture model. ACM Transactions on Graphics, 2014, 33(5): 171 DOI:10.1145/2645703


Yan X, Jiang Y T, Li C F, Martin R R, Hu S M. Multiphase SPH simulation for interactive fluids and solids. ACM Transactions on Graphics, 2016, 35(4): 79 DOI:10.1145/2897824.2925897


Feng G, Liu S G. Detail-preserving SPH fluid control with deformation constraints. Computer Animation and Virtual Worlds, 2018, 29(1): e1781 DOI:10.1002/cav.1781


Zhang X Y, Liu S G. Parallel SPH fluid control with dynamic details. Computer Animation and Virtual Worlds, 2018, 29(2): e1801 DOI:10.1002/cav.1801


Zhang X Y, Liu S G. SPH haptic interaction with multiple-fluid simulation. Virtual Reality, 2017, 21(4): 165–175 DOI:10.1007/s10055-017-0308-1


Zhang X Y, Liu S G. SPH fluid control with self-adaptive turbulent details. Computer Animation and Virtual Worlds, 2015, 26(3/4): 357–366 DOI:10.1002/cav.1637


Monaghan J J. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 1992, 30(1): 543–574 DOI:10.1146/annurev.aa.30.090192.002551


Müller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications. In: ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2003, 154–159


Becker M, Teschner M. Weakly compressible SPH for free surface flows. In: ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2007, 209–217


Solenthaler B, Pajarola R. Predictive-corrective incompressible SPH. ACM Transactions on Graphics, 2009, 28(3): 40 DOI:10.1145/1531326.1531346


Ihmsen M, Cornelis J, Solenthaler B, Horvath C, Teschner M. Implicit incompressible SPH. IEEE transactions on Visualization and Computer Graphics, 2013, 20(3): 426–435


Acharya S, Baliga B R, Karki K, Murthy J Y, Prakash C, Vanka S P. Pressure-based finite-volume methods in computational fluid dynamics. Journal of Heat Transfer, 2007, 129(4): 407–424 DOI:10.1115/1.2716419


Stam J. Stable fluids. In: SIGGRAPH, 1999,121–128


Bhatia H, Norgard G, Pascucci V, Bremer P T. The Helmholtz-hodge decomposition—A survey. IEEE Transactions on Visualization and Computer Graphics, 2013, 19(8): 1386–1404 DOI:10.1109/tvcg.2012.316


Chorin A J, Marsden J E. A mathematical introduction to fluid mechanics. New York, NY, Springer US, 1990 DOI:10.1007/978-1-4684-0364-0


Maria Denaro F. On the application of the Helmholtz-Hodge decomposition in projection methods for incompressible flows with general boundary conditions. International Journal for Numerical Methods in Fluids, 2003, 43(1): 43–69 DOI:10.1002/fld.598


Zapata M A U, van Bang D P, Nguyen K D. Parallel SOR methods with a parabolic-diffusion acceleration technique for solving an unstructured-grid Poisson equation on 3D arbitrary geometries. International Journal of Computational Fluid Dynamics, 2016, 30(5): 370–385 DOI:10.1080/10618562.2016.1234045


Helfenstein R, Koko J. Parallel preconditioned conjugate gradient algorithm on GPU. Journal of Computational and Applied Mathematics, 2012, 236(15): 3584–3590 DOI:10.1016/j.cam.2011.04.025


Koschier D, Bender J, Solenthaler B, Teschner M. Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids. In: Eurographics, 2020, 1–41


Akinci N, Ihmsen M, Akinci G, Solenthaler B, Teschner M. Versatile rigid-fluid coupling for incompressible SPH. ACM Transactions on Graphics, 2012, 31(4): 1–8 DOI:10.1145/2185520.2185558


He X W, Wang H M, Zhang F J, Wang H, Wang G P, Zhou K. Robust simulation of sparsely sampled thin features in SPH-based free surface flows. ACM Transactions on Graphics, 2014, 34(1): 1–9 DOI:10.1145/2682630


Akinci G, Ihmsen M, Akinci N, Teschner M. Parallel surface reconstruction for particle-based fluids. Computer Graphics Forum, 2012, 31(6): 1797–1809 DOI:10.1111/j.1467-8659.2012.02096.x


Bender J, Koschier D. Divergence-free smoothed particle hydrodynamics. In: Proceedings of the 14th ACM SIGGRAPH / Eurographics Symposium on Computer Animation. Los Angeles California, New York, NY, USA, ACM, 2015,147–155 DOI:10.1145/2786784.2786796