Particle-Based Fluid Simulation with SPH

Particle-Based Fluid Simulation with SPH
Published 05/23/2016

===========================================

Table of Contents

Introduction
The Equations of Fluid Motion
- A Brief Derivation…
Smoothed Particle Hydrodynamics
Extensions
Citation

## Introduction

Fluid simulation is a popular topic in computer graphics that uses computational fluid dynamics to generate realistic animations of fluids such as water and smoke. Most of these simulation techniques are simply numerical approximations of the Navier-Stokes equations, which describe the motion of fluids.

There are two primary approaches when numerically approximating the Navier-Stokes equations (or any problem involving flow fields): Lagrangian, or so-called "particle-based" methods, and Eulerian, or "grid-based" methods. Taking a Lagrangian perspective, we observe an individual discretized fluid parcel as it moves through space and time. Contrastingly, Eulerian solvers focus on the motion of a fluid through specific locations in space as time passes. One can think of modeling these Lagrangian parcels as particles and these Eulerian fixed locations in space as grid cells.

Though both Lagrangian and Eulerian solvers are widely used in computer graphics and computational fluid dynamics to great effect, a notable property of particle-based solvers is that they are generally not as spatially restricted since they do not rely entirely on an underlying Eulerian simulation grid. This, as well as their comparative ease of implementation, makes Lagrangian solvers uniquely well suited for applications such as video games, which often involve large, dynamic environments.

One of the most popular methods for the simulation of particle-based fluids is Smoothed Particle Hydrodynamics (SPH), first developed by Gingold and Monaghan in 1977 for astrophysics simulations. Müller first applied this same technique to real-time fluid simulation for games, and SPH has since grown into the method of choice for simulators with applications in engineering, disaster prevention, high-fidelity computer animation, and, of course, video games and other real-time effects.

## The Equations of Fluid Motion

Before discussing the Smoothed Particle Dynamics formulation of the Navier-Stokes equations in depth, it is necessary to build some amount of intuition by means of a derivation from more basic physics. I suggest a strong knowledge of multi-variable calculus, and a working familiarity with differential equations.

The Navier-Stokes equations are given by:

\rho(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}) = -\nabla p + \eta\nabla^2\mathbf{u}+\rho\mathbf{g}

subject to the incompressibility constraint

\nabla\cdot\mathbf{u}=0

where $\rho$ is the density of the fluid, $\mathbf{u}$ is the velocity vector field, $p$ is the pressure, and $\mathbf{g}$ is a vector of external forces (e.g. gravity).

### A Brief Derivation…

This can easily be understood from a more reasonable foundation of Newtonian physics using the continuity of mass, momentum, and energy. First, recall the definition of the material or Lagrangian derivative, which, as opposed to the Eulerian derivative with respect to a fixed position in space, quantifies the change in a field following a moving parcel. The material derivative is defined as the nonlinear operator

\frac{D}{Dt}\equiv\frac{\partial}{\partial t}+\mathbf{u}\cdot\nabla

where $\mathbf{u}$ is the flow velocity.

Recall the time-dependent classical relation between force and acceleration $\mathbf{F}=m\mathbf{a}$ . Applying the material derivative, we have

m\frac{D\mathbf{u}}{Dt}=\mathbf{F}^{\text{total}}

Which we can then expand to the forces acting on a particle

\rho\frac{D\mathbf{u}}{Dt}=\sum{\mathbf{F}}=\mathbf{F}^{\text{pressure}}+\mathbf{F}^{\text{viscosity}}+\rho\mathbf{g}

noting that the discretization of density limits to mass as

\rho=\lim_{\nabla V\rightarrow L^3}\frac{\Delta m}{\Delta V}

where the volume is approaching the size of the domain.

The contribution of the fluid to the total force on each particle is modeled in two parts: pressure and viscosity, with $\mathbf{g}$ again representing any external force such as gravity. Modeling the force contribution of pressure as the gradient of the pressures, and the force contribution of viscosity proportional to the divergence of the gradient of the velocity field, we have

\rho\frac{D\mathbf{u}}{Dt} = \rho(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}) = -\nabla p + \eta\nabla\cdot\nabla\mathbf{u}+\rho\mathbf{g}

We wish to model incompressible fluids, or fluids in which the material density is constant within each infinitesimal volume that moves with the flow velocity. Examples of incompressible fluids include liquids such as water, which are of great interest in computer graphics and physics simulation.

Following from this statement, we constrain the volume change of any arbitrary "chunk" of fluid $\Omega$ , motivated by the conservation of mass

\frac{d}{dt}\text{volume}(\Omega)=0

This volume change can be measured by integrating the normal component of the fluid velocity field along the boundary of our "chunk":

\iint_{\partial\Omega}\mathbf{u}\cdot\mathbf{n}=0

Gauss' Theorem gives

\iiint_{\Omega}\nabla\cdot\mathbf{u}=0\\ \therefore\nabla\cdot\mathbf{u}=0

While a formal derivation of the incompressible Navier-Stokes equations is much more involved, this should help to motivate some physical understanding of their origins and validity. Further reading could help show the fundamentals of pressure and viscosity forces, which can be better understood from the stress in the system (see the Stokes hypothesis or Cauchy's equation of motion).

## Smoothed Particle Hydrodynamics

The basic idea of SPH is derived from the integral interpolation, similar to Kernel Density Estimation. In essence, SPH is a discretization of the fluid into discrete elements, particles, over which properties are "smoothed" by a kernel function. This means that neighboring particles within the smoothing radius affect the properties of a given particle, such as pressure and density contributions—a surprisingly intuitive way of thinking about fluid dynamics simulation.

This cornerstone kernel function gives the approximation of any quantity $A$ at a point $\mathbf{r}$

A_s(\mathbf{r}) = \int A(\mathbf{r}')W(\mathbf{r}-\mathbf{r}',h)d\mathbf{r}' \approx \sum_{j}A_j\frac{m_j}{\rho_j}W(\mathbf{r}-\mathbf{r}_j,h)

where $m_j$ and $\rho_j$ are the mass and density of the j^th particle, respectively, and $W(\mathbf{r},h)$ is a radially symmetric smoothing kernel with length $h$ having the following properties:

W(-\mathbf{r},h)=W(\mathbf{r},h))\\ \int W(\mathbf{r})d\mathbf{r}=1

Note that we are summing over all other particles $j$ . Applying this to the above discussed Navier-Stokes equations, we first must determine a way to define the fluid density based upon neighboring particles, though this follows trivially from the previous approximation, substituting $\rho$ for the arbitrary quantity $A$

\rho_i = \rho(\mathbf{r}_i)=\sum_j m_j\frac{\rho_j}{\rho_j}W(\mathbf{r}_i-\mathbf{r}_j,h) = \sum_j m_jW(\mathbf{r}_i-\mathbf{r}_j,h)

Due to the linearity of the derivative $\nabla$ , the spatial derivtive (gradient) of any quantity can be easily obtained as follows

\nabla A(\mathbf{r})=\nabla\sum_jm_j\frac{A_j}{\rho_j}W(\mathbf{r}-\mathbf{r}_j,h)=\sum_jm_j\frac{A_j}{\rho_j}\nabla W(\mathbf{r}-\mathbf{r}_j,h)

The same is true for the Laplacian

\nabla^2 A(\mathbf{r})=\sum_jm_j\frac{A_j}{\rho_j}\nabla^2 W(\mathbf{r}-\mathbf{r}_j,h)

### SPH and Navier-Stokes

Following from the previous discussion of the Navier-Stokes equations, it is evident that the sum of the force density fields on the right hand side of the equation (pressure, viscosity, and external) give the change in momentum $\rho\frac{D\mathbf{u}}{Dt}$ of the particles on the left hand side. From this change in momentum, we can compute the acceleration of particle $i$

\mathbf{a}_i=\frac{d\mathbf{u}_i}{dt}=\frac{\mathbf{F}_i}{\rho_i}

We are therefore interested in approximating the force terms of the Navier-Stokes equations for pressure and viscosity with SPH as follows

\mathbf{F}^{\text{pressure}}_i=-\nabla p(\mathbf{r}_i)=-\sum_jm_im_j(\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2})\nabla W(\mathbf{r}_i-\mathbf{r}_j,h)\\ \mathbf{F}^{\text{viscosity}}_i=\eta\nabla^2\mathbf{u}(\mathbf{r}_i)=\eta\sum_jm_j\frac{\mathbf{u}_j-\mathbf{u}_i}{\mathbf{\rho}_j}\nabla^2W(\mathbf{r}_i-\mathbf{r}_j,h)

The pressure $p$ at a particle is calculated using some equation of state relating density to a defined rest density, usually the Tait equation or the ideal gas equation, such as

p=k(\rho-\rho_0)

where $\rho_0$ is the rest density and $k$ is some defined gas constant dependent on the temperature of the system.

Note that these formulations can change between different methodologies of implementation and more advanced techniques. For example, the direct application of SPH to the pressure term $-\nabla p$ gives

\mathbf{F}^{\text{pressure}}_i=-\nabla p(\mathbf{r}_i)=-\sum_jm_j\frac{p_j}{\rho_j}\nabla W(\mathbf{r}_i-\mathbf{r}_j,h)

which is not symmetric. Consider the case of two particles, wherein particle $i$ would use the pressure only of particle $j$ to compute its pressure force and vice versa. The pressures at two particle locations are not equal in general, therefore the force is not symmetric. The above SPH pressure force equation uses what has become a canonical symmetrization (weighted sum). In Müller's original paper, for reference, the symmetrization is performed using an arithmetic mean

\mathbf{F}^{\text{pressure}}_i=-\nabla p(\mathbf{r}_i)=-\sum_jm_j\frac{p_i+p_j}{2\rho_j}\nabla W(\mathbf{r}_i-\mathbf{r}_j,h)

This again has to be addressed in the viscosity term, which naively yields the asymmetric relation

\mathbf{F}^{\text{viscosity}}_i=\eta\sum_jm_j\frac{\mathbf{u}_j}{\rho_j}\nabla^2W(\mathbf{r}_i-\mathbf{r}_j,h)

This is addressed using velocity differences as a natural approach due to the viscosity force's dependence only on velocity differences and not absolute velocity which can be thought of as looking at the neighbors of particle $i$ from $i$ 's own moving frame of reference.

### Surface Tension

So far, we have given a brief justification for the Navier-Stokes from basic physical principles and have discussed the use of SPH as a Lagrangian discretization scheme, allowing us to model pressure and viscosity forces. A crucial missing component for a believable simulation, however, is surface tension. Surface tension is the result of attractive forces between neighboring molecules in a fluid. On the interior of a fluid, these attractive forces are balanced, and cancel out, whereas on the surface of a fluid, they are unbalanced, creating a net force acting in the direction of the surface normal towards the fluid, usually minimizing the curvature of the fluid surface. The magnitude of this force depends on the magnitude of the current curvature of the surface as well as a constant $\sigma$ depending on the two fluids that form the boundary.

Surface tension forces, while not present in the incompressible Naiver-Stokes equations, can be explicitly modeled in a number of different ways. Again, the most common and perhaps easiest to grasp approach is presented by Müller et al., which revises the Navier-Stokes equations by adding another term based upon the above physical principles for surface tension

\rho(\mathbf{r}_i)\frac{\mathbf{u}_i}{dt}=-\nabla p(\mathbf{r}_i)+\eta\nabla^2\mathbf{u}(\mathbf{r}_i)+\rho(\mathbf{r}_i)\mathbf{g}(\mathbf{r}_i)\boxed{-\sigma\nabla^2 c_s(\mathbf{r}_i)\frac{\nabla c_s(\mathbf{r}_i)}{|\nabla c_s(\mathbf{r}_i)|}}

where $c_s(\mathbf{r})$ is a color field $1$ at particle locations and $0$ everywhere else, given in a smoothed form by

c_s(\mathbf{r})=\sum_jm_j\frac{1}{\rho_j}W(\mathbf{r}-\mathbf{r}_j,h)

note that the gradient of the smoothed color field $\mathbf{n}=\nabla c_s$ yields the surface normal field pointing into the fluid and the divergence of $n$ measures the curvature of the surface

\mathbf{\kappa}=\frac{-\nabla^2 c_s}{\mathbf{n}}

Thus, we have the following force density

\mathbf{F}^{\text{surface}}=\sigma\kappa\mathbf{n}=-\sigma\nabla^2 c_s\frac{\mathbf{n}}{|\mathbf{n}|}

formed by distributing the surface traction among particles near the surface by multiplying a normalized scalar field which is non-zero only near the surface.

### Kernel Functions

SPH Kernels

Müller states "(s)tability, accuracy and speed of the SPH method highly depend on the choice of the smoothing kernels." The design of kernels for use in SPH is an area of active research and varies widely by implementation, but there are a number of core principles which are generally adhered to, well exemplified by Müller's original proposed kernels: poly6, spiky, and viscosity. The figure above shows these kernels (in order), with the thick line corresponding to their value, the thin line to their gradient towards the center, and the dashed line to their Laplacian.

Kernels with second-order interpolation errors can be constructed as even and normalized. Kernels that are $0$ with vanishing derivatives at the boundary are conductive to stability. Justification of this is left as an exercise to the reader.

Müller's poly6 kernel

W_{\text{poly6}}(\mathbf{r},h)=\frac{315}{64\pi h^9}\begin{cases} (h^2-r^2)^3 &\quad 0\leq r\leq h\\ 0 &\quad \text{otherwise} \end{cases}

satisfies the design constraints and can be evaluated without computing square roots in distance computations since $r$ only appears squared.

The spiky kernel

W_{\text{spiky}}(\mathbf{r},h)=\frac{15}{\pi h^6}\begin{cases} (h-r)^3 &\quad 0\leq r\leq h\\ 0 &\quad \text{otherwise} \end{cases}

originally proposed by Desburn for use in SPH deformable body simulation solves the problem of particle clumping when poly6 is used for pressure force computation by necessitating a non-zero gradient near the center.

Finally, the viscosity kernel

W_{\text{viscosity}}(\mathbf{r},h)=\frac{15}{2\pi h^3}\begin{cases} -\frac{r^3}{2h^3}+\frac{r^2}{h^2}+\frac{h}{2r}-1 &\quad 0\leq r\leq h\\ 0 &\quad \text{otherwise} \end{cases}

has the following properties

\begin{aligned} W(|\mathbf{r}|=h,h)&=0\\ \nabla W(|\mathbf{r}|=h,h)&=\mathbf{0}\\ \nabla^2W(\mathbf{r},h)&=\frac{45}{\pi h^6}(h-r) \end{aligned}

We desire a viscosity kernel with a smoothing effect only on the velocity field. For other kernels, there can be a non-zero negative Laplacian of the smoothed velocity field for two close particles, creating forces that increase particle relative velocities, creating artifacts and instability especially when the number of particles is low. The given viscosity kernel has a positive Laplacian everywhere and satisfies the desired constraints.

## Extensions

As discussed, vanilla (i.e. Müller-like) SPH formulations have a number of limitations, notably in the calculation of pressure and viscosity, artificial surface tension, and the relatively arbitrary choice of kernel functions. Most importantly, but beyond the scope of this brief introduction, are the problems with enforcing incompressibility and maintaining simulation stability with large timesteps with techniques such as PCISPH. Further extensions include using different equations of state for pressure calculation (see WCSPH), a discussion of applicable numerical integration methods, adaptive timestepping, and a myriad of performance optimizations (the astute reader might already see a need for sorting particles into a grid based upon the kernel support radius).

We plan to release a number of tutorials in the near future detailing both basic and advanced SPH implementations, as well as discussions of the mathematics behind some of the techniques used in bleeding-edge particle-based fluid solvers.

## Citation

Cited as:

Schuermann, Lucas V. (May 2016). Particle-Based Fluid Simulation with SPH. Writing. https://lucasschuermann.com/writing/particle-based-fluid-simulation

@misc{schuermann2016sph,
  title   = {Particle-Based Fluid Simulation with SPH},
  author  = {Schuermann, Lucas V.},
  year    = {2016},
  month   = {May},
  url     = "https://lucasschuermann.com/writing/particle-based-fluid-simulation"
}