Lecture 13 Fluid Simulation 01| The Particle-based (Lagrangian) Methods

Lagrangian Method = Particle-Based

# Incompressible Fluid Dynamics

usually compressible for explosion / shock wave / … ; incompressible for slower dynamics

# Forces for Incompressible Fluids

m\vb{a} = \vb{f} = \vb{f}_{\text{ext}} +\vb f_{\text{pres}} +\vb f_{\text{visc}}

# Incompressible Navier–Stokes Equations

The Navier-Stokes Equation ($p = k (\rho_0 - \rho)$)

\left\{\begin{array}{ll} \underbrace{\rho \frac{\mathrm{D}v}{\mathrm{D}t}} _{ma/V} = \underbrace{\rho g}_{mg /V} -\underbrace{\grad p}_{\text{adv}} + \underbrace{\mu \laplacian v}_{\text{diff}} & \text{(Momentum Equation)} \\ \div\ v = 0 \Leftrightarrow \frac{\mathrm{D}\rho}{\mathrm{D}t} = \rho (\div\ v ) = 0\ &\text{(Divergence free -- Mass conserving)} \end{array}\right.

# The Spatial Derivative Operators

• Gradient: \grad: \R^1 \rightarrow \R^3

  \grad s = \left[\frac{\partial s}{\partial x}, \frac{\partial s}{\partial y}, \frac{\partial s}{\partial z} \right]^{\mathrm{T}}

• Divergence: $\div$ : \R^3 \rightarrow \R^1

  $\div\ v = \frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}$

• Curl: $\curl $ : \R^3 \rightarrow \R^3

  \curl\ v = \left[\frac{\partial v_{z}}{\partial y}-\frac{\partial v_{y}}{\partial z}, \frac{\partial v_{z}}{\partial x}-\frac{\partial v_{x}}{\partial z}, \frac{\partial v_{y}}{\partial x}-\frac{\partial v_{x}}{\partial y}\right]^{\mathrm{T}}

• Laplace: \Delta = \laplacian = \div \grad: \R^n \rightarrow \R^n (Diffusion)

  \laplacian s = \div(\grad s) = \frac{\partial^{2} s}{\partial x^{2}}+\frac{\partial^{2} s}{\partial y^{2}}+\frac{\partial^{2} s}{\partial z^{2}}

# Time Discretization

# Operator Splitting

# General Steps

Separate the N-S Eqn into 2 parts and use half the timestep to solve each part (Advection-Projection)

• Step 1: Input $v_n$, output $v_{n+0.5}$: (Explicit, everything is known)
  • \rho \frac{\mathrm{D}v}{\mathrm{D}t} = \rho g + \mu \laplacian v
• Step 2: Input $v_{n+0.5}$, output $v_{n+1}$: (Implicit: $\rho$ and \grad p are unknown)
  • \rho \frac{\mathrm{D}v}{\mathrm{D}t} = -\grad p
  • $\div\ v = 0$

# Time Integration

Given $x_n$ & $v_n$

• Step 1: Advection / external and viscosity force integration
  • Solve: \mathrm{d}v = g + v\laplacian v_n
  • Update: $v_{n+0.5} = v_n + \Delta t \ \mathrm{d}v$
• Step 2: Projection / pressure solver (Poisson Solver)
  • Solve: \mathrm{d}v = -\frac{1}{\rho} \grad(k(\rho - \rho_0)) and $\frac{\mathrm{D} \rho}{\mathrm{D}t} = \div(v_{n+0.5}+\mathrm{d}v) = 0$
  • Update: $v_{n+1} = v_{n+0.5} + \Delta t \ \mathrm{d}v$
• Step 3: Update position
  • Update: $x _{n+1} = x_n + \Delta t\ v_{n+1}$
• Return $x_{n+1}$, $v_{n+1}$

# Integration with the Weakly Compressible (WC) Assumption

Storing the density $\rho$ as an individual var that advect with the vel field

\frac{\mathrm D v}{\mathrm D t}=g-\frac{1}{\rho} \nabla p+v \nabla^{2} v\ ;\quad \cancel{\div\ v = 0}\ \text{(weak compressible)}

# Time Integration

Symplectic Euler Integration: In this case step 1 and 2 can be combined since they are both explicit (no order diff)

• Step 1: Advection / external and viscosity force integration
  • Solve: \mathrm{d}v = g + v\laplacian v_n
  • Update: $v_{n+0.5} = v_n + \Delta t \ \mathrm{d}v$
• Step 2: Projection / pressure solver (Use the current-point-evaluated $\rho$)
  • Solve: \mathrm{d}v = -\frac{1}{\rho} \grad(k(\rho - \rho_0))
  • Update: $v_{n+1} = v_{n+0.5} + \Delta t \ \mathrm{d}v$
• Step 3: Update position
  • Update: $x _{n+1} = x_n + \Delta t\ v_{n+1}$
• Return $x_{n+1}$, $v_{n+1}$

# Spatial Discretization (Lagrangian View)

• Previous knowledge using Lag. View: Mesh-based simulation (FEM)
• Today: Mesh-free simulation => a simular example: marine balls

# Basic Idea

Cont. view -> Discrete view (using particles):

\frac{d v_{i}}{d t}=\underbrace{g-\frac{1}{\rho} \nabla p\left(x_{i}\right)+\nu \nabla^{2} v\left(x_{i}\right)}_{a_i}\ , \text{ where }\nu = \frac{\mu}{\rho_0}

Time integration (Symplectic Euler):

$v_i = v_i + \Delta t\ a_i$ & $x_i = x_i + \Delta t\ v_i$

But still need to evaluate $\rho(x_i)$, \grad p(x_i) & \laplacian v(x_i)

# Smoothed Particle Hydrodynamics (SPH)

# Dirac Delta

# Trivial Indentity

The Dirac function only tends to infinity at 0 but equals to 0 otherwise. And its overall integral is 1.

f(r) = \int^{\infty}_{-\infty} f(r') \delta (r-r')\ \mathrm{d}r\\ \delta(r)=\left\{\begin{array}{l} +\infty, \text { if } r=0 \\ 0, \text { otherwise } \end{array} \text { and } \int_{-\infty}^{\infty} \delta(r)\ \mathrm{d} r=1\right.

# Widen the Dirac Delta

$f(r) \approx \int f\left(r^{\prime}\right) W\left(r-r^{\prime}, h\right) d r^{\prime}, \text { where } \lim _{h \rightarrow 0} W(r, h)=\delta(r)$

$W(r, h)$ - kernel funciton:

• Symmetric: $W(r, h) = W(-r, h)$
• Sum to unity: $f(r) \approx \int f\left(r^{\prime}\right) W\left(r-r^{\prime}, h\right) d r^{\prime}, \text { where } \lim _{h \rightarrow 0} W(r, h)=\delta(r)$
• Compact support: $W(r, h ) = 0, \text{ if } |r|>h$ (Sampling radius - $h$)

e.g. $W(r, h) = \left\{\begin{aligned}&\frac{1}{2h},\text{ if } |r|<h\\ &0, \text{otherwise} \end{aligned} \right.$ (Error: $\mathcal{O}(h^2)$, can be reduced by decreasing the sampling radius $h$)

# Finite Probes: Summation

f(r) \approx \int f\left(r^{\prime}\right) W\left(r-r^{\prime}, h\right)\ \mathrm{d} r^{\prime} \approx \sum_{j} V_{j} f\left(r_{j}\right) W\left(r-r_{j}, h\right)

# Smoother Kernel Function

Still use the summation: $f(r) \approx \sum_{j} V_{j} f\left(r_{j}\right) W\left(r-r_{j}, h\right)$. But “trust” the closer probes (the orange & yellow ones) better => Smoother $W$ function

Use a \C^3 cont. spline: e.g. (The figure below shows the 1st / 2nd / 3rd order der)

$W(r, h)=\sigma_{d} \begin{cases}6\left(q^{3}-q^{2}\right)+1 & \text { for } 0 \leq q \leq \frac{1}{2} \\ 2(1-q)^{3} & \text { for } \frac{1}{2} \leq q \leq 1 \\ 0 & \text { otherwise }\end{cases}\ ; \quad \text{with } q=\frac{1}{h}\|r\|, \sigma_{1}=\frac{4}{3 h}, \sigma_{2}=\frac{40}{7 \pi h^{2}}, \sigma_{3}=\frac{8}{\pi h^{3}}$

# Smoothed Particle Hydrodynamics (SPH)

# Discretization

1D:

$f(r) \approx \sum_{j} V_{j} f\left(r_{j}\right) W\left(r-r_{j}, h\right)=\sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) W\left(r-r_{j}, h\right)$

2D:

Par stored in every particles

• Intrinsic quantities:
  • $h$: support radius
  • $\tilde{h}$: particle radius -> $V$: particle volume
• Time varying quantities:
  • $\rho$: density
  • $v$: velocity
  • $x$: position

# Evaluate 2D fields using the SP

$f(r_1) \approx \frac{m_{2}}{\rho_{2}} f\left(r_{2}\right) W\left(r_{1}-r_{2}, h\right) + \frac{m_{3}}{\rho_{3}} f\left(r_{3}\right) W\left(r_{1}-r_{3}, h\right) + \frac{m_{4}}{\rho_{4}} f\left(r_{4}\right) W\left(r_{1}-r_{4}, h\right)$

# SPH Spatial Derivatives

The operators will affect only on the kernel func ($W$)

• $\nabla f(r) \approx \sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) \nabla W\left(r-r_{j}, h\right)$
• \nabla\cdot \vb{F}(r) \approx \sum_{j} \frac{m_{j}}{\rho_{j}} \vb{F}\left(r_{j}\right) \cdot \nabla W\left(r-r_{j}, h\right)
• \nabla \times \boldsymbol{F}(r) \approx-\sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) \times \nabla W\left(r-r_{j}, h\right)
• $\nabla^2 f(r) \approx \sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) \nabla^2 W\left(r-r_{j}, h\right)$

# Improving Approximations for Spatial Derivatives

• $f(r) \approx \sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) W\left(r-r_{j}, h\right)$
• $\nabla f(r) \approx \sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) \nabla W\left(r-r_{j}, h\right)$
• Let $f(r)\equiv 1$: (Red line)
  • $1 \approx \sum_{j} \frac{m_{j}}{\rho_{j}} W\left(r-r_{j}, h\right)$ in the graph: not exactly 1.0
  • 0 \approx \sum_{j} \frac{m_{j}}{\rho_{j}} \grad W\left(r-r_{j}, h\right) not equal = 0 exactly (even if increase sampling points)

# Anti-symmetric Form

Better for projection / divergence / …

• => Since \grad f(r) \equiv \grad f(r) \cdot 1:

  • \grad f(r) = \grad f(r) \cdot 1 + f(r) \cdot \grad 1
  • Or equivalently: \grad f(r) = \grad f(r) - f(r) \cdot \grad 1 => the 2 grads can be evaluated using SPH

• \grad f(r) \approx \sum_{j} \frac{m_{j}}{\rho_{j}} f\left(r_{j}\right) \grad W\left(r-r_{j}, h\right)-f(r) \sum_{j} \frac{m_{j}}{\rho_{j}} \grad W\left(r-r_{j}, h\right) \approx \sum_{j} m_{j} \frac{f\left(r_{j}\right)-f(r)}{\rho_{j}} \nabla W\left(r-r_{j}, h\right) (Anti-symmetric form, green line)

  

# Symmetric Form

Better for forces

• A more gen case:
  • \grad\left(f(r) \rho^{n}\right)=\grad f(r) \cdot \rho^{n}+f(r) \cdot n \rho^{n-1} \grad\rho
  • \Rightarrow \grad f(r)=\frac{1}{\rho^{n}}\left(\grad\left(f(r) \cdot \rho^{n}\right)-f(r)\cdot n \cdot \rho^{n-1} \grad\rho\right)
• \grad f(r) \approx \sum_{j} m_{j}\left(\frac{f\left(r_{j}\right) \rho_{j}^{n-1}}{\rho^{n}}-\frac{n f(r)}{\rho}\right) \grad W\left(r-r_{j}, h\right)
• Special Case: when $n=-1$:
  • \grad f(r) \approx \rho \sum_{j} m_{j}\left(\frac{f\left(r_{j}\right)}{\rho_{j}^{2}}+\frac{f(r)}{\rho^{2}}\right) \grad W\left(r-r_{j}, h\right) (Symmetric form)

# Implementation Details (WCSPH)

# Simulation Pipeline

• For i in particles:

  • Search for neighbors j -> Apply the support radius $h$

    

• For i in particles:

  • Sample the $v$ / $\rho$ fields using SPH
    • Desity: $\rho_{i}=\sum_{j} \frac{m_{j}}{\rho_{i}} \rho_{j} W\left(r_{i}-r_{j}, h\right)=\sum_{j} m_{j} W_{i j}$ ($W_{ij}$ is for $r_i - r_j$)
    • Viscosity: \nu\grad^{2} v_{i}=v \sum_{j} m_{j} \frac{v_{j}-v_{i}}{\rho_{j}} \grad^{2} W_{i j} (Asymmetric form)
    • Pressure Gradient: $-\frac{1}{\rho_{i}} \nabla p_{i}=-\frac{\rho_{i}}{\rho_{i}} \sum_{j} m_{j}\left(\frac{p_{j}}{\rho_{j}^{2}}+\frac{p_{i}}{\rho_{i}^{2}}\right) \nabla W_{i j}, \text { where } p=k\left(\rho_{j}-\rho_{0}\right)$
  • Compute $f$ / $a$ using N-S Eqn
    • \frac{\mathrm{d} v_{i}}{\mathrm{d} t}=g-\frac{1}{\rho_{i}} \grad p_{i}+v \grad^{2} v_{i}

• For i in particles (Symplectic Euler)

  • Update $v$ using $a$
    • $v_i = v_i + \Delta t \cdot \frac{\mathrm{d}v_i}{\mathrm{d}t}$
  • Update $x$ using $v$
    • $x_i = x_i + \Delta t \cdot v_i$

# Boundary Conditions

# Problems of Boundaries

• Insufficient samples

  • Free Surface: lower density and pressure => Generate outward pressure

    

    -> Sol: Clamp the negative pressure (everywhere), $p = \max (0, k(\rho- \rho_0))$ (only be compressed, cannot be streched)\

  • Solid Boundary: lower density and pressure => Fluid leakage (outbound velocity)

    -> Sol: $p = \max (0, k(\rho-\rho_0))$

    for leakage:

    

    1. Refect the outbound veclocity when close to boundary
    2. Pad a layer of SP underneath the boundaries: $\rho_{\mathrm{solid}} = \rho_0$ & $v_{\mathrm{solid}} = 0$ => increases numerical viscosity

# Neighbor Search

Navie search: $\mathcal{O}(n^2)$ time

=> Use background grid: Common grid size = $h$ (same support radius in SPH)

(Each neighbor search takes 9 grids in 2D and 27 grids in 3D)