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

**Pedro Ferro Pereira**Co-author

**Fábio Cruz**Co-author

**David Carvalho**Co-author | Editor

**Ivan Pombo**Technical Reviewer

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

reallycould 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 “**P**articles” in S**P**H

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 “**S**moothed” in **S**PH:

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:

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]:

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^{2} 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 “**H**ydrodynamics” in SP**H**:

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:

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:

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.

## Recent posts from our blog

Maria Castro

The Fluid Cube dataset contains 100 fluid dynamics simulations of a fluid block with varying viscosity and initial shape, velocity and positions, flowing inside a unit cube domain.

Ivan Pombo

This dataset contains 1000 voltage measurements obtained through simulation of current propagation inside a circular two-dimensional domain for distinct electrical conductivity profiles.