html <span data-en="Smoothed Particle Hydrodynamics (SPH)" data-i18n="1" data-zh="Smoothed Particle Hydrodynamics (SPH)">Smoothed Particle Hydrodynamics (SPH)</span>
Language
01 / 12
SPH Basic Concept -- designed by Yi Zhan

Smoothed Particle Hydrodynamics(SPH)

Smoothed Particle Hydrodynamics represents continuum fields by particles and evaluates field variables, gradients, and divergence operators through local kernel-weighted interactions.

Lagrangian mesh-free method Kernel approximation Particle approximation Navier–Stokes equations
“在不可逆的熵增时空里,每一对粒子相遇都是物理学奇迹”
02 / 12
Starting Point

The original Dirac delta idea

Consider any continuum field variable $\phi(\mathbf{x})$, such as density, pressure, or one velocity component. The value at a target position $\mathbf{x}$ can be written as a weighted integral over the whole domain:

(1)$$\phi(\mathbf{x})=\int_{\Omega}\phi(\mathbf{x}')\,\delta(\mathbf{x}-\mathbf{x}')\,\mathrm{d}\mathbf{x}'$$

Here, $\delta(\mathbf{x}-\mathbf{x}')$ acts as an ideal weighting function: it assigns zero weight to all positions except $\mathbf{x}'=\mathbf{x}$, where the contribution is concentrated.

Def.$$\delta(\mathbf{x}-\mathbf{x}')=0\quad(\mathbf{x}'\ne\mathbf{x}),\qquad \int_{\Omega}\delta(\mathbf{x}-\mathbf{x}')\,\mathrm{d}\mathbf{x}'=1$$

This is the key idea: a field value can be represented through a weighted contribution of surrounding field information. The delta function is the limiting case of a perfect selector.

Ideal point selection

The delta function is not a regular finite-valued function. It is a generalized function with unit integral and infinitely concentrated support at the target point.

03 / 12
Why Replace δ?

Why the delta function is not directly usable

The exact identity is mathematically elegant, but it is not directly suitable for numerical continuum mechanics.

  • Singular amplitude:the delta function has an infinite peak at the target point.
  • No finite neighborhood:it does not provide a practical way to use neighboring particles.
  • Discrete sensitivity:if no particle is located exactly at the target point, the point selection becomes numerically fragile.
  • Differential operators:fluid mechanics requires gradients and divergence terms, which need smooth and stable approximations.
SPH therefore replaces the ideal selector $\delta$ with a smooth, normalized, finite-support kernel $W$.

From ideal selection to local averaging

Dirac delta $\delta$
ideal but singular
Kernel $W(\mathbf{r},h)$
smooth, non-negative, finite support
Particle summation
practical SPH approximation
$$\delta(\mathbf{x}-\mathbf{x}')\;\longrightarrow\;W(\mathbf{x}-\mathbf{x}',h)$$
04 / 12
Kernel Approximation

Kernel approximation: a smooth weighting function

Replacing the delta function by a kernel gives the SPH kernel approximation:

(2)$$\langle \phi(\mathbf{x})\rangle=\int_{\Omega}\phi(\mathbf{x}')\,W(\mathbf{x}-\mathbf{x}',h)\,\mathrm{d}\mathbf{x}'$$

The smoothing length $h$ controls the size of the local support domain. A larger $h$ uses more neighboring information and gives a smoother estimate; a smaller $h$ gives a sharper, more local estimate.

Kernel conditions$$W\ge 0,\qquad \int_{\Omega}W(\mathbf{r},h)\,\mathrm{d}\mathbf{r}=1,\qquad W(\mathbf{r},h)=0\ \text{for}\ |\mathbf{r}|>\kappa h$$

The 3D sketch shows the kernel as a positive weight surface. The peak is located at the target particle, and the value decays outward.

Interactive 3D kernel surface

0.52

Drag the plot to rotate the view. The kernel value is always non-negative and rises upward from the support plane.

05 / 12
Particle Approximation

Kernel-weighted particle summation

After the kernel integral is discretized, the field value at the target particle $a$ is obtained by adding the weighted contribution of each neighboring particle $b$.

(3)$$\int_{\Omega}\phi(\mathbf{x}')W(\mathbf{x}-\mathbf{x}',h)\,\mathrm{d}\mathbf{x}'\approx \sum_b \phi_b W_{ab}\Delta V_b$$
(4)$$\boxed{\langle \phi_a\rangle=\sum_b \underbrace{\frac{m_b}{\rho_b}\phi_b W_{ab}}_{\text{contribution of particle }b}}$$

The interaction panel explicitly shows the accumulation process: select a particle $b_k$, compute its contribution, and add it to the partial sum.

Step-by-step kernel summation

current particleb = --
current term--
partial sum--
Contribution: $(m_b/\rho_b)\phi_b W_{ab}$

Blue particles are outside the support; green particles have already been added; orange is the current particle being accumulated.

06 / 12
Gradient Approximation

Gradient as accumulated pairwise contributions

For a scalar field such as pressure or density, the SPH gradient is obtained by summing the contribution from each neighboring particle.

(5)$$\nabla \phi(\mathbf{x})\approx \int_{\Omega}\phi(\mathbf{x}')\nabla W(\mathbf{x}-\mathbf{x}',h)\,\mathrm{d}\mathbf{x}'$$
(6)$$\boxed{\langle \nabla \phi_a\rangle=\sum_b \underbrace{\frac{m_b}{\rho_b}(\phi_b-\phi_a)\nabla_a W_{ab}}_{\text{vector contribution of particle }b}}$$

Each particle contributes a small vector. The final gradient is the vector sum of these local contributions.

Step-by-step gradient accumulation

current particleb = --
current vector--
partial gradient--
Contribution: $(m_b/\rho_b)(\phi_b-\phi_a)\nabla_a W_{ab}$

The orange arrow is the current vector term; the green arrow is the accumulated gradient at particle $a$.

07 / 12
Divergence Approximation

Divergence as accumulated velocity projection

For velocity divergence, each neighboring particle contributes a scalar obtained by projecting the velocity difference onto the kernel-gradient direction.

(7)$$\left\langle \nabla\!\cdot\!\mathbf{v}_a \right\rangle = \sum_b C_b, \qquad C_b = \frac{m_b}{\rho_b}\left(\mathbf{v}_b-\mathbf{v}_a\right)\!\cdot\!\nabla_a W_{ab}$$

If the partial sum is positive, the local particle configuration behaves like expansion; if negative, it behaves like compression.

Step-by-step divergence accumulation

current particleb = --
current scalar--
partial divergence--
Contribution: $(m_b/\rho_b)(\mathbf{v}_b-\mathbf{v}_a)\cdot\nabla_a W_{ab}$

Velocity arrows are prescribed for illustration; the displayed number shows how each neighbor changes $\nabla\cdot\mathbf{v}_a$.

08 / 12
Governing Equations

From accumulated operators to Navier–Stokes terms

The same particle-by-particle accumulation is inserted into the Lagrangian Navier–Stokes equations. The continuity equation uses velocity divergence, while the momentum equation uses pressure-gradient and viscous terms.

(8)$$\frac{\mathrm{D}\rho}{\mathrm{D}t}=-\rho\nabla\cdot\mathbf{v},\qquad \frac{\mathrm{D}\mathbf{v}}{\mathrm{D}t}=-\frac{1}{\rho}\nabla p+\nu\nabla^2\mathbf{v}+\mathbf{g}$$
For weakly compressible SPH, density is evolved by continuity, pressure is obtained from an equation of state, and acceleration is then accumulated from pressure, viscosity, and body-force terms.

Step-by-step NS term assembly

current particleb = --
current term--
partial result--
Select continuity or pressure acceleration to see how each neighbor contributes to the governing equations.

This panel connects the abstract operators to the actual terms used in the SPH time-marching loop.

09 / 12
Closure

Equation of state for closure

In weakly compressible SPH, pressure is not obtained from a Poisson equation. Instead, it is computed from density using an equation of state.

(11)$$\boxed{p=B\left[\left(\frac{\rho}{\rho_0}\right)^\gamma-1\right]},\qquad B=\frac{c_0^2\rho_0}{\gamma}$$

Here, $\rho_0$ is the reference density, $c_0$ is the artificial speed of sound, and $\gamma$ is commonly taken as $7$ for water-like fluids.

This closes the equation system: once $\rho$ is evolved from continuity, $p$ can be calculated, and then the momentum equation can update velocity.

Closed WCSPH loop

Continuity equation
update $\rho_a$
Equation of state
compute $p_a=p(\rho_a)$
Momentum equation
update $\mathbf{v}_a$ and $\mathbf{x}_a$
10 / 12
SPH Form of NS: Mass and Pressure

SPH discretization of mass and pressure terms

Density summation

(12)$$\rho_a=\sum_b m_b W_{ab}$$

This reconstructs density from the local mass distribution.

Continuity equation

(13)$$\frac{\mathrm{D}\rho_a}{\mathrm{D}t}=\sum_b m_b(\mathbf{v}_a-\mathbf{v}_b)\cdot\nabla_a W_{ab}$$

This form evolves density according to pairwise velocity differences.

Pressure from EOS

(14)$$p_a=\frac{c_0^2\rho_0}{\gamma}\left[\left(\frac{\rho_a}{\rho_0}\right)^\gamma-1\right]$$

Pressure-gradient force

(15)$$\left(\frac{\mathrm{D}\mathbf{v}_a}{\mathrm{D}t}\right)_p=-\sum_b m_b\left(\frac{p_a}{\rho_a^2}+\frac{p_b}{\rho_b^2}\right)\nabla_a W_{ab}$$

The symmetric form helps preserve pairwise force balance.

11 / 12
SPH Form of NS: Viscosity and Time Evolution

SPH discretization of viscous and body-force terms

Viscous diffusion

(16)$$\left(\frac{\mathrm{D}\mathbf{v}_a}{\mathrm{D}t}\right)_\nu=\sum_b m_b\frac{4\nu\,\mathbf{r}_{ab}\cdot\nabla_a W_{ab}}{(\rho_a+\rho_b)(r_{ab}^2+\eta^2)}(\mathbf{v}_a-\mathbf{v}_b)$$

This is a common laminar-viscosity form. $\eta^2$ is a small regularization term to avoid singularity.

Complete acceleration

(17)$$\frac{\mathrm{D}\mathbf{v}_a}{\mathrm{D}t}=\left(\frac{\mathrm{D}\mathbf{v}_a}{\mathrm{D}t}\right)_p+\left(\frac{\mathrm{D}\mathbf{v}_a}{\mathrm{D}t}\right)_\nu+\mathbf{g}$$

Particle advection

(18)$$\frac{\mathrm{D}\mathbf{x}_a}{\mathrm{D}t}=\mathbf{v}_a$$

Definitions

$$\mathbf{r}_{ab}=\mathbf{x}_a-\mathbf{x}_b,\qquad r_{ab}=|\mathbf{r}_{ab}|,\qquad \nabla_a W_{ab}=\frac{\partial W(\mathbf{x}_a-\mathbf{x}_b,h)}{\partial \mathbf{x}_a}$$
12 / 12
Take-home Messages

The complete conceptual chain

Dirac delta
ideal selection
Kernel $W$
smooth weighting
Particle summation
SPH operators
Navier–Stokes
mass + momentum
EOS closure
$p=p(\rho)$

1. Field reconstruction

SPH reconstructs a continuum field from local kernel-weighted particle information.

2. Operators

Gradient and divergence operators are obtained by differentiating the kernel, not the noisy particle field directly.

3. Fluid solver

The discretized continuity equation, EOS, and momentum equation form a closed WCSPH time-marching system.

SPH Advantages

  • Explicit and highly parallelizable:local particle interactions make SPH well suited for large-scale GPU acceleration.
  • No explicit free-surface tracking:naturally handles breaking waves, splashing, droplets, and violent free-surface deformation.
  • Mesh-free Lagrangian description:suitable for large deformation, large displacement, and strongly nonlinear fluid–structure interaction.

SPH Limitations

  • Relatively low consistency:basic SPH approximations are often close to first-order accuracy without correction.
  • Complex boundary treatment:solid walls, open boundaries, and pressure boundaries usually require dedicated treatments.
  • Particle disorder and noise:irregular particle distributions may induce pressure oscillations, numerical noise, and artificial dissipation.
Core idea$$\boxed{\text{SPH} = \text{kernel-weighted field reconstruction} + \text{particle-based Navier--Stokes evolution}}$$