Flow #2: A Smooth Introduction to Smoothed Particle Hydrodynamics (SPH)

We are embarking on an epic journey. We’re on our way to build our own Machine Learning model that can predict the dynamics of a fluid of a simple (yet rich) physical scenario, as computed by a simulator based on Smoothed Particle Hydrodynamics (SPH).

In the previous post, we started off by explaining the basic physical principles behind the Mother of All Equations in Fluid Dynamics – the Navier-Stokes Equations (NSEs).

What’s next? In this post, we will get deeper into what SPH is all about — from its principles to technical methods. Then, it’s Fluid Dynamics in Particle Land. We will go through how particles of a fluid can be simulated à la SPH.

So, let’s flow!

Flow #2:
A Smooth Introduction to Smoothed Particle Hydrodynamics (SPH)

We can reflect on the point we originally raised in Flow #1 about our glass of water, containing an estimated \(10^{24}\) particles:

If we had infinite compute power and sophisticated CFD methods, we really could simulate the dynamics of the fluid down to the level of each individual particle.

A framework whose implementation can lead us towards this goal is Smoothed Particle Hydrodynamics (SPH).

Video 1: SPH simulations of various scenarios using SPlisHSPlasH. As we will see, the inclusion of custom bodies and complex geometries make this particle-based method very versatile. Credits: Jan Bender and Dan Koshier

Time to make all these notions more technically concise [1].
Here’s your smooth SPH 101:

Smoothed particle
what?

In essence, SPH refers to a range of mesh-free computational methods that solve dynamical relations between particles. The seminal work that introduced the foundations of SPH dates back to 1977 by Gingold & Monaghan and Lucy to Astrophysics problems of non-spherical stars and fission processes.

Due to its extremely versatility, the field quickly developed to extend to various other domains: think of Fluid Dynamics, Ballistics, Coastal Dynamics or even Fracture Mechanics!

Generally speaking, as long as phenomena can be modelled as a continuum, SPH can perform sophisticated computations to the dynamics of all particles used to represent this medium.

Fig. 1: (left) Rendering of output coming from a SPH simulation of the fluid flow that results from a rotating submerged turbine in the sea. (right) Density of particles making up a mechanical part. Credits: Koschier et al [1]

But what do we mean by a particle here?
Welcome to Particle Land!

The “Particles” in SPH
discretising the flow

SPH focuses on discrete entities - suggestively termed particles – and their interactions. These particles are not necessarily the actual physical atoms or molecules! In essence, these are nodal points where fluid properties are evaluated and updated in time and space.

The way these properties (such as their mass, temperature, pressure or any other quantity) change is expressed in accordance to equations governing the particle phenomena.

The aim of SPH is to solve how these features change at the particles’ positions only. This is the way we discretise our medium: a particle is a unit which stores, at a time instant \(t\), and position \(\mathbf{x}(t)\), features \(f(\mathbf{x}(t))\).

Video 2: Visualising the speed feature \(|\mathbf{v}(\mathbf{x}(t))|\) in a fluid of SPH particles close to a fluid-fluid interface. Each point represents the particle position and it is coloured according to how fast it is. Lighter (darker) colours show faster (slower) particles. Credits: David Carvalho / Inductiva

The “Smoothed” in SPH:
understanding particle averaging

In Particle Land, the notion of fluid does not go away. SPH relies precisely on particle averaging in order to compute local features of the fluid dynamics.

The way SPH estimates quantities locally is quite intuitive. The value of a function \(f\) at a position \(\mathbf{x}\) can be expressed with the integral representation:

\[f(\mathbf{x}) = \int f(\mathbf{x}') \delta(\mathbf{x} - \mathbf{x'}) d \mathbf{x'}\]

where \(\delta(\mathbf{x})\) is the Dirac delta distribution.

So far, this result does not help us much. This is an infinitely-spiked distribution at \(\mathbf{x}= \mathbf{0}\) (i.e. \(\delta(\mathbf{0}) = \infty\)) and 0 elsewhere. A way to approximate this unreachable result is to relax this steepness by introducing a more tractable distribution \(W(\mathbf{x}, h)\) [2]:

\[f(\mathbf{x}) \approx \int f(\mathbf{x}')W(\mathbf{x}-\mathbf{x}', h) \mathrm{d}\mathbf{x}'\]

Here, the \(W(\mathbf{x})\) distribution is the kernel function and \(h\) a smoothening length parameter such that \(\lim_{h \rightarrow 0} W(\mathbf{x}, h) = \delta(\mathbf{x})\). The role of \(h\) is to widen the distribution, making it differentiable.

Fig. 2: A SPH kernel is a distribution that defines the weighted average of the value of a function with respect to its neighbours. Credits: Wikimedia Commons / Jlcercos

Which kernels to use? A reasonable candidate could be the normal distribution \(\mathcal{N}(\mu, \sigma, \mathbf{x}) \propto e^{-(\mathbf{x}-\mu)/(2\sigma)}\) with mean \(\mu = 0\) and variance \(\sigma^2\).

However, its infinitely-extending tails never vanish, meaning all particle contributions would need to be computed. Computationally, it is simply undesirably expensive.

In practice, more sophisticated kernels are used instead [14]. For instance, bell-shaped splines for which \(W(\mathbf{x}-\mathbf{x}')=0\) for \(||\mathbf{x}-\mathbf{x}'|| \geq h\) are truncated kernels which have the desirable finite support.

Fig. 3: C² Wendland kernel [3] in 2 dimensions. It is compactly supported by the unit disc.

Okay – so far we have only slightly tweaked the integrand in the continuum integral representation of \(f(\mathbf{x})\).
We now need to discretise it à la SPH.

From continuous
… to discrete

Let’s get practical with SPH. For a given position \(\mathbf{x}(t)\), after applying a kernel \(W(\mathbf{x})\), a specific neighbourhood \(\mathcal{N}_{\mathbf{x}}\) arises, whose particles we know their values \(\mathbf{x}_j\) from \(j=1\) to some (hopefully large) \(j=N\).

To approximate the integral, we find that each finite volume \(d \mathbf{x}\) has a density \(\rho(\mathbf{x}) = dm / d \mathbf{x}\), whenever we know the particles masses \(m_j\).

We now sum over the neighbourhood \(\mathcal{N}(\mathbf{x})\) to obtain our local estimate of \(f(\mathbf{x})\). Shortening \(W_{ij} \equiv W(\mathbf{x}_i-\mathbf{x}_j)\):

\[f(\mathbf{x}_i) \approx \sum_j f_j W_{ij} \frac{m_j}{\rho_j} \ .\]

So there you have it: in SPH, functions are smoothened by applying a convolution with this specified kernel!

The recipe is clear: with a set kernel \(W_{ij}\) and knowing the particles positions \(\mathbf{x}_i\), the density can be computed as \(\rho(\mathbf{x}_i) = \sum_j W_{ij} m_j\) and then, any quantity \(f\) as shown above.

The “Hydrodynamics” in SPH:
The NSEs - Particle Edition

We need to find how to translate the NSEs into this particle-based world. Like we discussed before, incompressible Newtonian fluids particles must evolve like what the Navier-Stokes Equations (NSEs) dictate:

\[\underbrace{\rho \frac{D \mathbf{v}}{D t}}_{\rm 'total \ force'} = -\underbrace{\nabla P}_{\rm pressure} + \underbrace{\mu \nabla^2 \mathbf{v}}_{\rm viscosity} + \underbrace{\rho \mathbf{g}}_{\rm gravity} + \underbrace{\mathbf{F}^{\rm ext}}_{\rm external \ forces}\]

Solving the equation means we must know five quantities at the particle location \(\mathbf{x}(t)\): the pressure \(P\), density \(\rho\) and the velocity vector \(\mathbf{v} = (v_x, v_y, v_z)\) for all particles \(i\).

Once discretised à la SPH, and with the volume \(d \mathbf{x}_j = m_j/\rho_j\), the operators that we need could be approximated as [1]:

• Gradient: \(\nabla f_i \approx \sum_j f_j \nabla W_{ij}(m_j/\rho_j)\),
• Divergence: \(\nabla \cdot f_i \approx \sum_j f_j \cdot \nabla W_{ij}(m_j/\rho_j)\),
• Laplacian: \(\nabla^2 f_i \approx \sum_j f_j \nabla^2 W_{ij}(m_j/\rho_j)\)

Unfortunately, these simplistic SPH-based formulas lead to poor quality discretisation, leading to unstable simulations. For this reason, the SPH community contributes with better and computationally-efficient approximations [1,4] tailored to each particular contribution in the NSEs and physical regime.

Many of these are rather intricate [5]. For example, Divergence-Free SPH (DFSPH)[7] solvers for incompressible fluids or Weakly Compressible SPH [8] pressure solvers, as well as XSPH [9] viscosity solvers [10].

Operator
splitting

We need to be very careful in the way we discretize and evolve each differential operator. If our numerical updates fail to reasonably estimate the next steps, error buildup will inevitably misestimate the trial solution [5].

By applying operator splitting, the underlying governing equation (which, as we saw above, include terms pertaining to wildly different phenomena and scales) are decomposed in several sequential problems requiring their own individual techniques [6].

In our case, the viscosity and pressure contributions to the acceleration are computed in separate loops over the particles and the velocities partially updated mid-loop like this pseudo-code outlines:

Fig. 4: Pseudo-code outlining the typical loop of a Smoothed Particle Hydrodynamics simulation.

The idea is simple. Firstly, the density \(\rho_i\) is computed at each particle position. Then, each acceleration term in the NSEs is calculated using a suitable discretisation.

At this stage, all terms have been computed at time \(t_i\). In order to update the positions and velocities at \(t_{i+1}\), a time integration scheme must be applied (e.g. Euler integration) until the desired simulation duration is reached.

Handling
boundaries

So far we have mentioned on how to solve the NSEs in the interior of the domain. But how do we handle boundary conditions within the SPH framework? In other words, how do we include walls and other rigid bodies, as well as complex surfaces and model their interactions with the fluid?

The versatility in how boundaries can be included in a customisable way in the simulation is one of the strongest points in support of SPH as a CFD framework.

Video 3: SPH simulation of a wave energy harvester (a mechanical device to convert kinetic to electric energy) with the DualSPHysics simulator (in conjunction with the Project Chrono library for rigid-body interactions). Credits: DualSPHysics.

So, how is this implemented?

In particle-based methods [11], the boundaries are promoted to particles, although they are forced to be static. These boundary particles then contribute as any other neighbor when computing fluid properties within the SPH formulation.

How exactly? Think of a region \(\mathcal{F}\) where the fluid can flow close to a boundary region \(\mathcal{B}\) and of a particle \(i\) about to reach the boundary. The boundary effects need to be considered whenever any fluid particle neighbourhood \(\mathcal{N}_i = \{ x \in \mathbb{R}^3 \ | \ \left \| \mathbf{x}_i − \mathbf{x} \right \| ≤ h \}\) intersects with the boundary region \(\mathcal{B}\).

Fig. 5: Discretising a boundary interface with particles. When computing the SPH integral for particle \(i\), all particles within its neighbourhood \(\mathcal{N}_i\) contribute, both from the boundary region \(\mathcal{B}\) and fluid region \(\mathcal{F}\). Credits: Pedro Pereira / Inductiva.

When a particle approaches a boundary, the SPH integral can be split into the fluid and boundary regions:

\[f(\mathbf{x}_i) \approx \underbrace{\sum_{j^\mathcal{F}} f_{j^\mathcal{F}} W_{i j^\mathcal{F}} \left( \frac{m_{j^\mathcal{F}}}{\rho_{j^\mathcal{F}}}\right)}_{\rm Fluid} + \underbrace{\sum_{j^\mathcal{B}} f_{j^\mathcal{B}} W_{i j^\mathcal{B}} \left( \frac{m_{j^\mathcal{B}}}{\rho_{j^\mathcal{B}}}\right)}_{\rm Boundary}\]

Even though boundaries are incorporated in a consistent way, the added computational cost from both holding the boundary particles in memory and searching for neighbouring particles is high.

Also, determining a good boundary discretisation is not trivial. If the setup is too sparse, the boundary may not be sufficiently covered and fluid particles can penetrate it; if too dense, boundaries may become too bumpy and lead to artifacts in the fluid flow at the interface.

Not only that, inherent modelling deficiencies may introduce a lot of unphysical phenomena. For example, Modified Dynamic Boundary Conditions (mDBCs) have been used to correct artificial output from other boundary handling algorithms [15].

Fig. 6: Artefacts found in SPH simulations when handling the fluid-boundary interactions. Credits: A. English [15].

Towards supercomputation
(without a supercomputer)

Performing computations at the individual particle level is a heavy task and we should not look down upon this task, especially when scaling to larger batches of particles.

Things are moving fast though. As GPUs become more powerful and allow for multi-processing, the potential of SPH methods to push for more refined solutions becomes tangible.

For instance, using NVIDIA’s CUDA, GPU-enabled simulations using SPH have shown tremendous speed-ups (up to 2 orders of magnitude in comparison to their CPU-enabled counterparts) [13]. As stated there:

GPU implementation of SPH permits high resolution SPH modelling in hours and days rather than weeks and months on inexpensive and readily available hardware.

Video 4: SPH simulation of the interaction of a large wave with an oil rig using the MultiGPU implementation of DualSPHysics. In order to obtain the dynamics of a billion particles during 12 seconds, 64 Tesla 2090 GPUs were used for 92 hours. Credits: DualSPhysics

All of this can be performed at a fraction of the running costs associated to supercomputing clusters.

Better still, the availability of well-maintained open-source libraries — like SPlisHSPlasH, DualSPHysics, GPUSPH or AQUAgpuSPH — is democratizing the way high-quality flow data in highly-complex systems can be realistically obtained by virtually anyone!

In the next
post…

Well, the field of Smoothed Particle Hydrodynamics is barely 40 years old. From a computational point of view, it is still at its infancy. Now that we understand how the algorithms work and to what they apply to, we need to dive deeper into the computational side of how these SPH are deployed.

As you would probably expect, the mathematical structure of the data we need to know can be conceptualised in the language of… graphs.

Efficient graph algorithms provide a promising way to increase the finesse of SPH simulations. That’s what we will dive into.

Until then…
~~~ flow away! ~~~

---

References
& Remarks

[1] For an in-depth review of the SPH framework, check Koschier et al “Smoothed particle hydrodynamics techniques for the physics based simulation of fluids and solids.” (2020) and Monaghan’s “Smoothed particle hydrodynamics.” Annual review of astronomy and astrophysics 30.1 (1992): 543-574.

[2] The distributions, the kernels \(W(\mathbf{x},h)\) must also be normalised (i.e. its integral across the entire domain \(\Omega\) must be \(\int_{\Omega} W(\mathbf{x},h) d \mathbf{x} = 1\)), non-negative i.e. \(W(\mathbf{x},h) > 0 \ \forall \mathbf{x}\), symmetric i.e. \(W(\mathbf{x},h) = W(-\mathbf{x},h)\) and (twice) differentiable.

[3] An example is the Wendland Quintic C2 kernel: \(W(q) = \alpha_d \begin{cases} (1 - q/2)^4(2q + 1) & \text{if } 0 \leq q \leq 2\\ 0 & \text{if } q > 2 \end{cases} \ ,\) where \(q=||\mathbf{x}||/h\) and \(\alpha_d\) is a scaling constant dependent on the number of spatial dimensions \(d\).

[4] Cossins, Peter J., “Smoothed particle hydrodynamics.” arXiv preprint arXiv:1007.1245 (2010)

[5] We have already touched upon this issue in the context of solving the Heat Equation on a 2D plate, where a FTCS scheme was used to propagate the differential opertator solutions in space and time.

[6] An example of these is the modified Laplacian discretisation:

\[\nabla^2 f_i \approx 2(d + 2) \sum_j \frac{m_j}{\rho_j} \frac{(f_i - f_j) \cdot (\mathbf{x}_i - \mathbf{x}_j)}{||(\mathbf{x}_i - \mathbf{x}_j)||^2 + \epsilon}\nabla W_{ij},\]

where \(d\) is the number of spatial dimensions and \(\epsilon = 0.01h^2\) is added to avoid zero-division singularities.

[7] Bender, Jan, and Dan Koschier. “Divergence-free smoothed particle hydrodynamics.” Proceedings of the 14th ACM SIGGRAPH/Eurographics symposium on computer animation. 2015

[8] M. Becker and M. Teschner, “Weakly compressible sph for free surface flows” in Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation, 2007, pp. 209–217.

[9] Schechter, Hagit, and Robert Bridson. “Ghost SPH for animating water” ACM Transactions on Graphics (TOG) 31.4 (2012): 1-8.

[10] M. Weiler, D. Koschier, M. Brand, and J. Bender, “A physically consistent implicit viscosity solver for sph fluids,” in Computer Graphics Forum, vol. 37, no. 2. Wiley Online Library, 2018, pp. 145–155.

[11] Akinci, Nadir, et al. “Versatile rigid-fluid coupling for incompressible SPH.” ACM Transactions on Graphics (TOG) 31.4 (2012): 1-8.

[12] Dan Koschier, Jan Bender, Barbara Solenthaler, Matthias Teschner - “Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids”

[13] Hérault, Alexis, Giuseppe Bilotta, and Robert A. Dalrymple. “Sph on GPU with CUDA.” Journal of Hydraulic Research 48.sup1 (2010): 74-79.

[14] An extensive yet informal discussion on kernels can be found in:
Fasshauer, Gregory E. “Positive definite kernels: past, present and future.” Dolomites Research Notes on Approximation 4 (2011): 21-63.

[15] English, A., et al. “Modified dynamic boundary conditions (mDBC) for general-purpose smoothed particle hydrodynamics (SPH): Application to tank sloshing, dam break and fish pass problems.” Computational Particle Mechanics 9.5 (2022): 1-15.