Smoothed-Particle Hydrodynamics

Basic Ideas

Smoothed Particle Hydrodynamics (SPH) is a Lagrangian mesh-free particle method whose core idea is to discretize a continuous fluid medium into a collection of particles carrying physical attributes (such as mass, density, velocity, pressure, etc.). The interactions between these particles approximate the governing equations of fluid motion, such as the Navier–Stokes equations. Unlike traditional grid-based Eulerian methods (e.g., finite difference or finite volume methods), SPH does not rely on a fixed computational mesh; instead, it tracks the motion of each fluid particle to simulate the evolution of the entire fluid. This unique Lagrangian nature grants SPH inherent advantages when dealing with complex geometric boundaries, free surfaces, large deformations, and multiphase flows.

Why Particle System?

SPH adopts a particle-based computational framework primarily due to the unique advantages of particle discretization in describing fluid motion.

In SPH, each particle carries physical information such as mass, position, velocity, density, and pressure, and updates itself through interactions with neighboring particles via kernel functions. The particle’s density is obtained by a weighted sum of neighbor masses, while accelerations are determined by pressure gradients, viscous forces, and external forces.

Since particles move with the fluid, SPH naturally incorporates convection, avoiding the numerical diffusion and instability found in Eulerian discretization of convective terms — simplifying equation solving while improving stability and accuracy.

SPH can directly impose boundary conditions using “boundary particles,” significantly reducing the cost of mesh generation and adaptation, providing high flexibility when handling complex geometries and dynamic boundaries.

SPH excels at simulating free-surface flows, such as waves or splashing droplets. The free surface naturally emerges from the distribution of the outermost particles, enabling the method to capture complex interface evolutions such as breaking and merging without extra treatment.

Kernel Function: From Continuous to Discrete

The kernel function smooths or interpolates discrete particle information into continuous space. The choice and properties of the kernel function directly affect the accuracy, stability, and computational efficiency of SPH.

Dirac Delta Function

The core idea of SPH is to approximate continuous fields (such as density or velocity fields) using discrete particles. A kernel function measures how much a particle at $x_i$ influences the field value at some point $x$, with the influence decreasing with distance $r$ and controlled by a smoothing length $h$.

Mathematically, SPH is based on approximations of the Dirac delta function. The $\delta$-function $\delta(x - x’)$ is a generalized function that is infinite at $x = x’$, zero elsewhere, and integrates to 1 over the whole space. Any continuous function $f(x)$ can be expressed through convolution with the delta function:

\[f(x) = \int f(x') \delta(x - x') \,\mathrm{d} x'.\]

In SPH, the kernel function $W_h(\mathbf{x} - \mathbf{x}’)$ is used as an approximation of the delta function, so we expect the following properties:

\[W(\mathbf{0})=1,\quad W(\infty)=0,\quad W(r)\downarrow \text{ as } r\uparrow,\quad \int W=1.\]

Thus, the basic SPH formulation follows:

\[\begin{align} f(x) & = \int f(x')\delta(x-x') \,\mathrm{d}x' \notag\\ & = \int f(x')W_{h}(x-x') \,\mathrm{d}x' +o(h) \notag\\ & \approx \int f(x') W_h(x - x') \,\mathrm{d}x' \notag \end{align}\]

The approximation error is $O(h^2)$.

Calculus on Particle Interpolation

In the discretized particle system, the continuous fluid is replaced by $N$ particles. A particle $i$ has position $x_i$, mass $m_i$, and local density $\rho_i$.

Mathematically, the core SPH operation is approximating a continuous field $f(\mathbf{x})$ by a weighted sum of particle attributes:

\[\begin{align} f(\mathbf{x}) & \approx \int f(x') W_h(x - x') \,\mathrm{d}x' \notag \\ & = \int \frac{f(\mathbf{x}')}{\rho(\mathbf{x}')}W_{h}(\mathbf{x}-\mathbf{x'})(\rho(\mathbf{x}') d\mathbf{x}') \notag \end{align}\]

In the discrete view, the mass element $\rho(x’) \,\mathrm{d}x’$ is replaced by particle mass $m_i$, yielding:

\[f(\mathbf{x}) \approx \sum_{i} \frac{f(\mathbf{x}_i)}{\rho_{i}} W_{h}(\mathbf{x}-\mathbf{x}_{i}) m_{i}. \tag{1}\]

Gradient

Using interchangeability of differentiation and integration:

\[\begin{align} \nabla f(\mathbf{x}) & \approx \nabla\left[\int f(\mathbf{x}') W_h(\mathbf{x} - \mathbf{x}') \,\mathrm{d}\mathbf{x}'\right] \notag\\ & = \int f(\mathbf{x}') \nabla_\mathbf{x}\left[W_h(\mathbf{x} - \mathbf{x}') \right] \,\mathrm{d}\mathbf{x}' \notag\\ & \approx \sum_{i} \frac{f(\mathbf{x}_i)}{\rho_{i}} \nabla W_{h}(\mathbf{x}-\mathbf{x}_{i}) m_{i} \tag{2} \end{align}\]

Divergence

For a vector field $\mathbf{F}(\mathbf{x})$, its divergence is:

\[\begin{align} \nabla_\mathbf{x} \cdot [\mathbf{F}(\mathbf{x}') W_h(\mathbf{x} - \mathbf{x}')] & = \nabla_\mathbf{x} W_h(\mathbf{x} - \mathbf{x}') \cdot \mathbf{F}(\mathbf{x}') \notag \end{align}\]

Thus:

\[\begin{align} \nabla\cdot\mathbf{F}(\mathbf{x}) & \approx \int \nabla W_h(\mathbf{x} - \mathbf{x}') \cdot \mathbf{F}(\mathbf{x}')\,\mathrm{d}\mathbf{x}' \notag\\ & \approx \sum_{i} \frac{\mathbf{F}(\mathbf{x}_i)}{\rho_{i}} \cdot\nabla W_{h}(\mathbf{x}-\mathbf{x}_{i}) m_{i} \tag{3} \end{align}\]

Curl

Similarly, for curl:

\[\begin{align} \nabla_\mathbf{x} \times [\mathbf{F}(\mathbf{x}') W_h(\mathbf{x} - \mathbf{x}')] & = \nabla_\mathbf{x} W_h(\mathbf{x} - \mathbf{x}') \times \mathbf{F}(\mathbf{x}') \notag \end{align}\]

Thus:

\[\begin{align} \nabla\times\mathbf{F}(\mathbf{x}) & \approx \int \nabla W_h(\mathbf{x} - \mathbf{x}') \times \mathbf{F}(\mathbf{x}')\,\mathrm{d}\mathbf{x}' \notag\\ & \approx \sum_{i} \frac{\mathbf{F}(\mathbf{x}_i)}{\rho_{i}} \times\nabla W_{h}(\mathbf{x}-\mathbf{x}_{i}) m_{i} \tag{4} \end{align}\]

Implicit Estimation

In practice, gradient estimates suffer from particle irregularity, making
$\sum_i (m_i/\rho_i)\nabla_x W_h(x-x_i)\neq 0$, causing a ghost force and yielding an $O(h^2)$ error even for constant functions.

A better method couples $f$ and $\rho$ using:

\[\nabla(f\rho^n) =nf\rho^{n-1}\nabla \rho + \rho^n\nabla f \Rightarrow \nabla f=\rho^{-n}[\nabla(f\rho^n)-nf\rho^{n-1}\nabla \rho]\]

Substituting into Eq.(1):

\[\nabla f(\mathbf{x}) \approx \rho^{-n}(\mathbf{x}) \sum_{i} m_{i}(f(\mathbf{x}_{i})\rho(\mathbf{x}_{i})^{n-1}-nf(\mathbf{x})\rho(\mathbf{x})^{n-1})\nabla W_{h}(\mathbf{x}-\mathbf{x}_{i}) \tag{5}\]

For $n=\pm 1$, symmetry appears:

\[\begin{align} n&=1 & \nabla f&=\rho^{-1} \sum_{i} m_{i}(f-f_{i})\nabla W_{h} \tag{6}\\ n&=-1 & \nabla f&=\rho \sum_{i}m_{i} \left( \frac{f}{\rho^2}+ \frac{f_{i}}{\rho_{i}^2}\right) \nabla W_{h} \tag{7} \end{align}\]

Similarly, divergence estimates:

\[\begin{align} n&=1 & \nabla \cdot \mathbf{F}&=\rho^{-1} \sum_{i} m_{i}(\mathbf{ F}-\mathbf{F}_{i})\cdot\nabla W_{h} \tag{8}\\ n&=-1 & \nabla \cdot \mathbf{F} &=\rho \sum_{i}m_{i} \left( \frac{\mathbf{F}}{\rho^2}+ \frac{\mathbf{F}_{i}}{\rho_{i}^2}\right) \cdot\nabla W_{h}\tag{9} \end{align}\]