Smoothed-particle hydrodynamics (SPH) is a computational method used for simulating the mechanics of continuum media, such as solid mechanics and fluid flows. It was developed by Gingold and Monaghan and Lucy in 1977, initially for astrophysical problems. It has been used in many fields of research, including astrophysics, ballistics, volcanology, and oceanography. It is a meshfree Lagrangian method (where the co-ordinates move with the fluid), and the resolution of the method can easily be adjusted with respect to variables such as density.
Method
Advantages
By construction, SPH is a meshfree method, which makes it ideally suited to simulate problems dominated by complex boundary dynamics, like free surface flows, or large boundary displacement. The lack of a mesh significantly simplifies the model implementation and its parallelization, even for many-core architectures.
SPH can be easily extended to a wide variety of fields, and hybridized with some other models, as discussed in Modelling Physics. As discussed in section on weakly compressible SPH, the method has great conservation features. The computational cost of SPH simulations per number of particles is significantly less than the cost of grid-based simulations per number of cells when the metric of interest is related to fluid density (e.g., the probability density function of density fluctuations). This is the case because in SPH the resolution is put where the matter is.
Limitations
Setting boundary conditions in SPH such as inlets and outlets and walls is more difficult than with grid-based methods. In fact, it has been stated that "the treatment of boundary conditions is certainly one of the most difficult technical points of the SPH method". This challenge is partly because in SPH the particles near the boundary change with time. Nonetheless, wall boundary conditions for SPH are available. The computational cost of SPH simulations per number of particles is significantly larger than the cost of grid-based simulations per number of cells when the metric of interest is not (directly) related to density (e.g., the kinetic-energy spectrum). Therefore, overlooking issues of parallel speedup, the simulation of constant-density flows (e.g., external aerodynamics) is more efficient with grid-based methods than with SPH.
Examples
Fluid dynamics
Smoothed-particle hydrodynamics is being increasingly used to model fluid motion as well. This is due to several benefits over traditional grid-based techniques. First, SPH guarantees conservation of mass without extra computation since the particles themselves represent mass. Second, SPH computes pressure from weighted contributions of neighboring particles rather than by solving linear systems of equations. Finally, unlike grid-based techniques, which must track fluid boundaries, SPH creates a free surface for two-phase interacting fluids directly since the particles represent the denser fluid (usually water) and empty space represents the lighter fluid (usually air). For these reasons, it is possible to simulate fluid motion using SPH in real time. However, both grid-based and SPH techniques still require the generation of renderable free surface geometry using a polygonization technique such as metaballs and marching cubes, point splatting, or 'carpet' visualization. For gas dynamics it is more appropriate to use the kernel function itself to produce a rendering of gas column density (e.g., as done in the SPLASH visualisation package).
One drawback over grid-based techniques is the need for large numbers of particles to produce simulations of equivalent resolution. In the typical implementation of both uniform grids and SPH particle techniques, many voxels or particles will be used to fill water volumes that are never rendered. However, accuracy can be significantly higher with sophisticated grid-based techniques, especially those coupled with particle methods (such as particle level sets), since it is easier to enforce the incompressibility condition in these systems. SPH for fluid simulation is being used increasingly in real-time animation and games where accuracy is not as critical as interactivity.
Recent work in SPH for fluid simulation has increased performance, accuracy, and areas of application:
B. Solenthaler, 2009, develops Predictive-Corrective SPH (PCISPH) to allow for better incompressibility constraints M. Ihmsen et al., 2010, introduce boundary handling and adaptive time-stepping for PCISPH for accurate rigid body interactions K. Bodin et al., 2011, replace the standard equation of state pressure with a density constraint and apply a variational time integrator R. Hoetzlein, 2012, develops efficient GPU-based SPH for large scenes with Fluids v.3 N. Akinci et al., 2012, introduce a versatile boundary handling and two-way SPH-rigid coupling technique that is completely based on hydrodynamic forces; the approach is applicable to different types of SPH solvers M. Macklin et al., 2013 simulates incompressible flows inside the Position Based Dynamics framework, for bigger timesteps N. Akinci et al., 2013, introduce a versatile surface tension and two-way fluid-solid adhesion technique that allows simulating a variety of interesting physical effects that are observed in reality J. Kyle and E. Terrell, 2013, apply SPH to Full-Film Lubrication A. Mahdavi and N. Talebbeydokhti, 2015, propose a hybrid algorithm for implementation of solid boundary condition and simulate flow over a sharp crested weir S. Tavakkol et al., 2016, develop curvSPH, which makes the horizontal and vertical size of particles independent and generates uniform mass distribution along curved boundaries W. Kostorz and A. Esmail-Yakas, 2020, propose a general, efficient and simple method for evaluating normalization factors near piecewise-planar boundaries Colagrossi et al., 2019, study flow around a cylinder close to a free-surface and compare with other techniques
Astrophysics
Smoothed-particle hydrodynamics's adaptive resolution, numerical conservation of physically conserved quantities, and ability to simulate phenomena covering many orders of magnitude make it ideal for computations in theoretical astrophysics.
Simulations of galaxy formation, star formation, stellar collisions, supernovae and meteor impacts are some of the wide variety of astrophysical and cosmological uses of this method.
SPH is used to model hydrodynamic flows, including possible effects of gravity. Incorporating other astrophysical processes which may be important, such as radiative transfer and magnetic fields is an active area of research in the astronomical community, and has had some limited success.
Solid mechanics
Libersky and Petschek extended SPH to Solid Mechanics. The main advantage of SPH in this application is the possibility of dealing with larger local distortion than grid-based methods. This feature has been exploited in many applications in Solid Mechanics: metal forming, impact, crack growth, fracture, fragmentation, etc.
Another important advantage of meshfree methods in general, and of SPH in particular, is that mesh dependence problems are naturally avoided given the meshfree nature of the method. In particular, mesh alignment is related to problems involving cracks and it is avoided in SPH due to the isotropic support of the kernel functions. However, classical SPH formulations suffer from tensile instabilities and lack of consistency. Over the past years, different corrections have been introduced to improve the accuracy of the SPH solution, leading to the RKPM by Liu et al. Randles and Libersky and Johnson and Beissel tried to solve the consistency problem in their study of impact phenomena.
Dyka et al. and Randles and Libersky introduced the stress-point integration into SPH and Ted Belytschko et al. showed that the stress-point technique removes the instability due to spurious singular modes, while tensile instabilities can be avoided by using a Lagrangian kernel. Many other recent studies can be found in the literature devoted to improve the convergence of the SPH method.
Recent improvements in understanding the convergence and stability of SPH have allowed for more widespread applications in Solid Mechanics. Other examples of applications and developments of the method include:
Metal forming simulations. SPH-based method SPAM (Smoothed Particle Applied Mechanics) for impact fracture in solids by William G. Hoover. Modified SPH (SPH/MLSPH) for fracture and fragmentation. Taylor-SPH (TSPH) for shock wave propagation in solids. Generalized coordinate SPH (GSPH) allocates particles inhomogeneously in the Cartesian coordinate system and arranges them via mapping in a generalized coordinate system in which the particles are aligned at a uniform spacing.
Numerical tools
Interpolations
The Smoothed-Particle Hydrodynamics (SPH) method works by dividing the fluid into a set of discrete moving elements , referred to as particles. Their Lagrangian nature allows setting their position by integration of their velocity as:
These particles interact through a kernel function with characteristic radius known as the "smoothing length", typically represented in equations by . This means that the physical quantity of any particle can be obtained by summing the relevant properties of all the particles that lie within the range of the kernel, the latter being used as a weighting function . This can be understood in two steps. First an arbitrary field is written as a convolution with :
The error in making the above approximation is order . Secondly, the integral is approximated using a Riemann summation over the particles:
where the summation over includes all particles in the simulation. is the volume of particle , is the value of the quantity for particle and denotes position. For example, the density of particle can be expressed as:
where denotes the particle mass and the particle density, while is a short notation for . The error done in approximating the integral by a discrete sum depends on , on the particle size (i.e. , being the space dimension), and on the particle arrangement in space. The latter effect is still poorly known.
Kernel functions commonly used include the Gaussian function, the quintic spline and the Wendland kernel. The latter two kernels are compactly supported (unlike the Gaussian, where there is a small contribution at any finite distance away), with support proportional to . This has the advantage of saving computational effort by not including the relatively minor contributions from distant particles.
Although the size of the smoothing length can be fixed in both space and time, this does not take advantage of the full power of SPH. By assigning each particle its own smoothing length and allowing it to vary with time, the resolution of a simulation can be made to automatically adapt itself depending on local conditions. For example, in a very dense region where many particles are close together, the smoothing length can be made relatively short, yielding high spatial resolution. Conversely, in low-density regions where individual particles are far apart and the resolution is low, the smoothing length can be increased, optimising the computation for the regions of interest.
Discretization of governing equations
For particles of constant mass, differentiating the interpolated density with respect to time yields
where is the gradient of with respect to . Comparing this equation with the continuity equation in the Lagrangian description (using material derivatives),
it is apparent that its right-hand side is an approximation of ; hence one defines a discrete divergence operator as follows:
This operator gives an SPH approximation of at the particle for a given set of particles with given masses , positions and velocities .