Smoothed Particle Hydrodynamics

1. Notes
2. Kernel
- 2.1. Cubic Spline 2D
- 2.2. Cubic Spline 2D - Modified

(This page is a work in progress)

SPH is a mesh-free method for discretization of functions and partial differential operators.

Functions are discreized into samples equipped with kernel function \(W\)

Some common myths about SPH are:

SPH particle represents
- a physical atom/molecule
- a droplet
- a grain (e.g. in sand simulation)
SPH is better than Eulerian approaches
Eulerian approaches are better than SPH
"Proper" engineering CFD methods are grid based
SPH is only 0th-order consistent

1. Notes

https://www.youtube.com/watch?v=Xby1yDiDgVE ranks opensource SPH software. SPHinXsys which is second on the list has good documentation that goes to the formulation it uses in detail. It also support multibody physics, so rigid and elastic solids interact with fluids seemlessly.

Paper References for SplishSplash: https://splishsplash.readthedocs.io/en/latest/bibliography.html
First rule of SPH: Interpret the Kernel as Gaussian
Second rule of SPH: Rewrite formula with density placed inside the operator (Smoothed Particle Hydrodynamics - Monaghan 1992.pdf: Page 3)
- use lagrange equation to derive discreized versions. you get energy conservation for free. https://www.youtube.com/watch?v=QZN_Kj8cTP4&t=458s

2. Kernel

Any kernel function works. Some necessary properties are:

Integral over the support area should be 1
Should be symmetric
Should be monotonically decreasing as distance increases

Good properties to have:

Smooth derivative

Smooth second derivative is also good to have but it doesn't matter much because we don't use that directly. Instead we use gradient of \(W\) to approximate higher order derivatives.

Sometimes different Kernels are used for different cases (e.g. for Cohesion, Adhesion forces).

2.1. Cubic Spline 2D

Smoothed Particle Hydrodynamics - Monaghan 1992.pdf: Page 12

\begin{equation*} W(r) = \frac {10} {7 \pi h^2} \begin{cases} 1 - \frac 3 2 q^2 (1 - q/2) & \textrm{ for } 0 \leq q \leq 1 \\ \frac 1 4 (2 - q)^3 & \textrm { for } 1 \leq q \leq 2 \\ 0 & \textrm { otherwise} \end{cases} \end{equation*}

where, \(q = r/h\)

If we integrate this function we get total area 1:

\begin{align*} 1 &= \int_0^1 \int_0^1 W(q) \diff x \diff y \\ &= \int_0^{2\pi} \int_0^1 W(q) r \diff r \diff \theta \\ &= 2\pi \int_0^{2h} W(r) r \diff r \end{align*}

< Collapse code block> Expand code block

(defun W (r h)
  "Cubic Spline Kernel 2D
r = 0 to 2h"
  (let ((q (/ r h)))
    (*
     (/ 10 7 pi (expt h 2))
     (cond ((<= 0 q 1)
            (- 1 (* 3/2 (expt q 2) (- 1 (/ q 2)))))
           ((<= 1 q 2)
            (* 1/4 (expt (- 2 q) 3)))
           (t 0.0d0)))))

(defun integrate (f s e step)
  (loop for x from s to e by step
        summing (* step (funcall f x))))

;; Test W
(let ((h 0.25))
  (* 2 pi (integrate (lambda (x) (* x (W x h)))
                     0 (* 2 h) 0.001)))

0.9999976896687713d0

2.2. Cubic Spline 2D - Modified

When using cubic spline 2d kernel given by monaghan, we have to consider neighbhours that are within \(\hbar = 2h\) radius of the particle. i.e. the support radius is \(2h\). to avoid confusing support radius (\(\hbar\)) with the kernel parameter \(h\), we can modify the above kernel as follows: (koschier 2019 pg. 3)

\begin{equation*} W(r) = \frac {40} {7 \pi \hbar^2} \begin{cases} 1 - 6 q^2 (1 - q) & \textrm{ for } 0 \leq q \leq 1/2 \\ 2 (1 - q)^3 & \textrm { for } 1/2 \leq q \leq 1 \\ 0 & \textrm { otherwise} \end{cases} \end{equation*}

where, \(q = r/\hbar\) and the above formula is obatined by changing the kernel parameter \(h\) to \(\hbar/2\) and substituting \(q\) with \(2q\) in the original formula.

Again the total area is equal to 1.

< Collapse code block> Expand code block

(defun W (r h)
  "Cubic Spline Kernel 2D
r = 0 to h"
  (let ((q (/ r h)))
    (*
     (/ 40 7 pi (expt h 2))
     (cond ((<= 0 q 1/2)
            (- 1 (* 6 (expt q 2) (- 1 q))))
           ((<= 1/2 q 1)
            (* 2 (expt (- 1 q) 3)))
           (t 0.0d0)))))

(defun integrate (f s e step)
  (loop for x from s to e by step
        summing (* step (funcall f x))))

;; Test W
(let ((h 0.5))
  (* 2 pi (integrate (lambda (x) (* x (W x h)))
                     0 h 0.001)))

0.9999976896687713d0

The derivative of kernel function is

\begin{equation*} W'(r) = \frac {40} {7 \pi \hbar^3} \begin{cases} - 6 (2q - 3q^2) & \textrm{ for } 0 \leq q \leq 1/2 \\ -6 (1 - q)^2 & \textrm { for } 1/2 \leq q \leq 1 \\ 0 & \textrm { otherwise} \end{cases} \end{equation*}

< Collapse code block> Expand code block

(defun dW/dr (r h)
  "First derivative of kernel function (Cubic Spline 2D)
  r = 0 to +h+"
  (let ((q (/ r h)))
    (* (/ 40 7 pi (expt h 3))
       (cond ((<= 0 q 1/2)
              (* -6 (- (* 2 q) (* 3 (expt q 2)))))
             ((<= 1/2 q 2)
              (* -6 (expt (- 1 q) 2)))
             (t 0d0)))))

(defun diff-at (x h step)
  (- (/ (- (W (+ x step) h)
           (W x h))
        step)
     (dw/dr x h)))

;; Test dW/dr
(let ((h 0.5))
  (values
   (diff-at (* 0.25 h) h 0.001)
   (diff-at (* 0.81 h) h 0.001)))

-0.04312213131925091d0, 0.033100545837216444d0

Smoothed Particle Hydrodynamics

Table of Contents

1. Notes

2. Kernel

2.1. Cubic Spline 2D

2.2. Cubic Spline 2D - Modified