Smoothed Particle Hydrodynamics

Kernel

SPH defines "smoothness" using a function called a "kernel". When a particle position is given, the kernel function spreads out any values stored in a the particle nearby, starting from the center point of the particle, the function fades out to zero as the distance from the center reaches the kernel radius.

$\displaystyle W_{std(r)} = \frac{315}{64{\pi}h^3} \left\{ \begin{array}{lr} (1 - \frac{r^2}{h^2})^2 & 0 \leq r \leq h \\ 0 & otherwise \end{array} \right.$

A valid kernel function should integrate to 1

$\displaystyle \int W(r) = 1$

Interpolation

The basic idea of SPH Interpolation is to measture any physical quantities at any given location by looking up the nearby particles.

$\displaystyle \phi(x) = m\sum_{j} \frac{\phi_{j}}{\rho_{j}}W(x - x_j) $

where $x, m, \phi, \rho$ and $W(r)$ are location for the Interpolation, mass, quantity that we want to interpolate, density and kernel function, respectively. Subscript $j$ represents the $j$th neighbouring particle

Density

From the interpolation equation above density is required to calculate any required qauntity, to calculate density we use the interpolation equation to get:

$\displaystyle \rho(x) = m\sum_{j} W(x - x_j) $

Differential Operators

We now have the basic tools to perform mathematical computation in the SPH world. To compute fluid dynamics, however we need the differential operators

Gradient

$\displaystyle \nabla\phi(x) = m \sum_{j} \frac{\phi_{j}}{\rho_{j}} \nabla{W(|x - x_{j}|)}$

This implementation of the gradient, however, is not symmetric. It means the gradient calculated from two nearby particles with respect to each other can be different. The implementation below resolve the symmetry problem

$\displaystyle \nabla\phi(x) = \rho_{i} m \sum_{j} \left( \frac{\phi_{i}}{\rho_{i}^2} + \frac{\phi_{j}}{\rho_{j}^2} \right) \nabla{W(|x - x_{j}|)}$

Laplacian

$\displaystyle \nabla^2 \phi{(x)} = m \sum_{j} \left( \frac{\phi_{j} - \phi_{i}}{\rho_{j}} \right) \nabla^2 W(x - x_{j})$

Dynamics

The pressure gradient, viscosity and gravity forces are the key components for implementing a fluid solver. The presure gradient force results in fluid flow from higher to lower pressure area and viscosity defines how think the fluid is.

$\displaystyle \mathbf{a} = \mathbf{g} - \frac{\nabla{p}}{\rho} + \mu \nabla^2 v$

SPH simulator takes the following steps:

  1. Measure the density with particles' current position
  2. Compute the pressure based on the density
  3. Compute the presure gradient force
  4. Compute the viscosity force
  5. Compute the gravity and other external forces
  6. Peform time integration

Pressure Gradient Force

Pressure is highly related to density - higer density generates higer pressure. to calculate pressure we are going to need to use the Equation of state

Equation of State

The equation of state or EOS, describes the relationship between state variables. In this case, we need to map density to pressure

$\displaystyle p = \frac{\kappa}{\gamma} \left( \frac{\rho}{\rho_{0}} - 1 \right)^\gamma $

Here $p$ is pressure, $\kappa$ is eosScale, $\gamma$ is eosExponent, $\rho$ is density, and $\rho_{0}$ is targetDensity

eosScale can be calculated using the formula below:

$\displaystyle \kappa = \rho_{0} \frac{c_{s}}{\gamma}$

where $c_{s}$ is speed of sound in the fluid.

Computing Pressure Gradient

To compute the pressure gradient, we can use the symmetric version of the gradient equation

$\displaystyle \mathbf{f}_{p} = -m \frac{\nabla{p}}{\rho}$

we can apply the symmetric gradient operator to get

$\displaystyle \mathbf{f}_{p} = -m^2 \sum_{j} \left( \frac{p_{i}}{\rho_{i}^2} + \frac{p_{j}}{\rho_{j}^2} \right) \nabla{W(x - x_{j})}$

Viscosity

The equation for the viscosity force can be written as:

$\displaystyle \mathbf{f}_{\upsilon} = m \mu \nabla^2 \mathbf{u}$

based on the Laplacian operator the equation above can be rewritten as:

$\displaystyle \mathbf{f}_{\upsilon}(x) = m^2 \mu \sum_{j} \left( \frac{u_{j} - u_{i}}{\rho_{j}} \right) \nabla^2 W(x - x_{j})$