光滑粒子流体动力学方法基本原理和求解固体力学问题的一般步骤

北京理工大学 | 李明健

1. 基本原理

光滑粒子流体动力学方法（Smoothed particle hydrodynamics，SPH）是计算流体力学、计算固体力学领域中一种常用的粒子类无网格数值模拟方法，它具有方便追踪自由界面、无网格畸变等优势。

SPH 方法采用核近似（Kernel approximation）和粒子近似（Particle approximation），对连续场变量和其梯度、散度在粒子上进行插值表示。下面介绍其基本思想和求解流程。

1.1 核近似

$\Omega$ $f(\boldsymbol r)$ $\boldsymbol r$ 为自变量的函数，该函数可用以下积分形式等效表示：

\begin{matrix} (1) & f (r) = \int_{Ω} f (r^{'}) δ (r - r^{'}) d r^{'} \end{matrix}

$\delta$ 是 delta 函数，其性质满足：

\begin{matrix} (2) & δ (r - r^{'}) = 0, r \neq r^{'} \end{matrix}

\begin{matrix} (3) & \int_{- \infty}^{\infty} δ (r - r^{'}) d (r - r^{'}) = 1 \end{matrix}

$W(\boldsymbol r - \boldsymbol r^\prime, h)$ $(\boldsymbol r - \boldsymbol r^\prime)$ $h$ ，其具有以下性质：

尖峰性（在零点处逼近 delta 函数的性质）：

\begin{matrix} (4) & lim_{h \to 0} W (r - r^{'}, h) = δ (r - r^{'}) \end{matrix}

归一性（用于质量等物理量时保证了守恒性）：
$\begin{matrix} (5) & \int_{- \infty}^{\infty} W (r - r^{'}) d (r - r^{'}) = 1 \end{matrix}$
紧支性（避免了全局求解，极大减少计算量）：

\begin{matrix} (6) & W (r - r^{'}) = 0, ∥ r - r^{'} ∥ > κ h \end{matrix}

对称性（通常是偶函数）、非负性（保证正确表示物理含义）、衰减性（距离增大时单调递减）、光滑性（保证数值稳定）。

这种光滑的加权函数，在 SPH 中称为核函数（kernel function），包括多种形式，下面给出两种三维问题常用的核函数：

（1）高斯核函数：

\begin{matrix} (7) & W = {\begin{cases} \frac{1}{π^{3 / 2} h^{3}} e^{- q^{2}}, & 0 ⩽ q ⩽ 3 \\ 0, & q > 3 \end{cases} q = \frac{∥ r - r^{'} ∥}{h} \end{matrix}

（2）三次样条核函数：

\begin{matrix} (8) & W = {\begin{cases} \frac{1}{π h^{3}} (1 - \frac{3}{2} q^{2} + \frac{3}{4} q^{3}), & 0 ⩽ q ⩽ 1 \\ \frac{1}{4 π h^{3}} (2 - q)^{3}, & 1 < q ⩽ 2 \\ 0, & q > 2 \end{cases} q = \frac{∥ r - r^{'} ∥}{h} \end{matrix}

$h$ $h$ 越大，则越平缓。取核函数以后，即可以将公式（1）写为以下形式：

\begin{matrix} (9) & f (r) \approx \int_{Ω} f (r^{'}) W (r - r^{'}, h) d r^{'} \end{matrix}

该式即为 SPH 中的积分表示，或称核近似。

1.2 粒子近似

下面，仍需要将核近似得到的公式进一步进行粒子近似，即用紧支域内的一组粒子来表示场量：

\begin{matrix} (10) & f (r) \approx \sum_{j = 1}^{N} V_{j} f (r_{j}^{'}) W (r - r_{j}^{'}, h) \end{matrix}

粒子近似 $\boldsymbol r_i$ $h$ $V_j$ $j$ $V_j = m_j / \rho_j$ ，类似地，对于矢量场，同样有：

\begin{matrix} (11) & f (r_{i}) \approx \sum_{j = 1}^{N} V_{j} f (r_{j}) W (r_{i} - r_{j}, h) \end{matrix}

$N$ $\Omega$ $\Omega$ 内对所有粒子求和，与在紧支域内的粒子求和是等价的。

对于标量场的梯度，可以表示为：

\begin{matrix} (12) & \nabla f (r_{i}) \approx \sum_{j = 1}^{N} V_{j} f (r_{j}) \nabla W (r_{i} - r_{j}, h) \end{matrix}

对于矢量场的散度，可以表示为：

\begin{matrix} (13) & \nabla \cdot f (r_{i}) \approx \sum_{j = 1}^{N} V_{j} f (r_{j}) \cdot \nabla W (r_{i} - r_{j}, h) \end{matrix}

这样，就把对场量的求导转化成对核函数的求导，使复杂偏微分方程的离散求解成为了可能，这与有限元法中将场量求导转化成对形函数求导的做法是类似的。

2. 固体力学问题求解

2.1 控制方程

在固体力学求解中，连续性方程（质量守恒）为：

\begin{matrix} (14) & \frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u) = 0 \end{matrix}

动量方程为：

\begin{matrix} (15) & \frac{\partial u}{\partial t} = \frac{1}{ρ} \nabla \cdot σ + g \end{matrix}

能量守恒方程（高速冲击中通常忽略热传导，主要考虑应力做功）：

\begin{matrix} (16) & \frac{\partial e}{\partial t} = \frac{1}{ρ} (- p \nabla \cdot u + s : \dot{ε}) \end{matrix}

$e$ $\boldsymbol s$ 为偏应力。

此外还需要满足几何方程、本构方程和状态方程。

2.2 求解流程

$N$ $i$ $m_i$ $\boldsymbol r_i$ $\boldsymbol u_i$ $\rho_i$ $\boldsymbol \sigma_i$ 等属性。随后，在粒子上对控制方程进行求解。每个求解时间步内的求解流程如下：

2.2.1 连续性方程：密度更新

首先更新粒子密度，通常有两种方式，第一种是通过核函数对粒子密度进行插值计算：

\begin{matrix} (17) & ρ_{i} = \sum_{j = 1}^{N} m_{j} W (∥ r_{i} - r_{j} ∥, h) \end{matrix}

根据核函数的归一性，该方程事实上隐含了质量守恒，这种处理简单高效，但当粒子分布不均匀时误差较大。

第二种方法是直接离散连续性方程（14），通过粒子速度场更新密度：

\begin{matrix} (18) & \frac{d ρ_{i}}{d t} = \sum_{j = 1}^{N} m_{j} (u_{i} - u_{j}) \cdot \nabla_{i} W_{i j} \end{matrix}

$\nabla_i W_{ij}$ $i$ 的梯度。随后再更新密度：

\begin{matrix} (19) & ρ_{i}^{n + 1} = ρ_{i}^{n} + Δ t \cdot \frac{d ρ_{i}}{d t} \end{matrix}

2.2.2 状态方程：压力计算

根据状态方程更新压力，对于高速冲击问题，可采用刚性气体状态方程（Stiffened Gas Equation of State），压力计算如下：

\begin{matrix} (20) & p_{i} = c_{0}^{2} (ρ_{i} - ρ_{0}) + (γ - 1) ρ_{i} e_{i} \end{matrix}

$p$ $c_0$ $\rho_0$ $\gamma$ $e$ 是比内能。随后更新声速：

\begin{matrix} (21) & c_{i} = \sqrt{c_{0}^{2} + (γ - 1) (e_{i} + \frac{p_{i}}{ρ_{i}})} \end{matrix}

$\gamma = 1$ ，即退化为等温状态方程（Isothermal Equation of State），即不考虑内能的变化对压力的影响。

2.2.3 几何方程：应变率/旋率计算

根据邻域内粒子上的速度，求解当前粒子上的速度梯度：

\begin{matrix} (22) & \nabla u_{i} \approx \sum_{j = 1}^{N} V_{j} (u_{j} - u_{i}) \otimes \nabla_{i} W_{i j} \end{matrix}

$\otimes$ 表示并矢（张量积），得到的速度梯度实际是一个二阶张量：

\begin{matrix} (23) & \nabla u_{i} = [\begin{array}{lll} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x} & \frac{\partial w}{\partial x} \\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} & \frac{\partial w}{\partial y} \\ \frac{\partial u}{\partial z} & \frac{\partial v}{\partial z} & \frac{\partial w}{\partial z} \end{array}] \end{matrix}

$u$ $v$ $w$ $\boldsymbol u$ $x$ $y$ $z$ 则为位矢的三个分量坐标。

根据速度梯度张量，即可计算应变率张量：

\begin{matrix} (24) & {\dot{ε}}_{i} = \frac{1}{2} [\nabla u_{i} + (\nabla u_{i})^{T}] \end{matrix}

以及计算旋率张量：

\begin{matrix} (25) & Ω_{i} = \frac{1}{2} [\nabla u_{i} - (\nabla u_{i})^{T}] \end{matrix}

2.2.4 本构方程：应力计算

随后根据应变率和旋率更新应力，对于线弹性材料，采用广义胡克定律得到偏应力变化率如下：

\begin{matrix} (26) & {\dot{s}}_{i} = 2 G [{\dot{ε}}_{i} - \frac{1}{3} tr ({\dot{ε}}_{i}) I] + s_{i} \cdot Ω_{i} + Ω_{i} \cdot s_{i} \end{matrix}

$G$ 是剪切模量，随后对偏应力进行更新：

\begin{matrix} (27) & s_{i}^{n + 1} = s_{i}^{n} + {\dot{s}}_{i} Δ t \end{matrix}

应力即为压力和偏应力的组合：

\begin{matrix} (28) & σ_{i} = s_{i} - p_{i} I \end{matrix}

2.2.5 能量方程：内能

考虑应力做功，将能量方程（16）离散为粒子形式：

\begin{matrix} (29) & \frac{d e_{i}}{d t} = \frac{1}{ρ_{i}} s_{i} : \dot{ε_{i}} + \sum_{j = 1}^{N} m_{j} (\frac{p_{i}}{ρ_{i}^{2}} + \frac{p_{j}}{ρ_{j}^{2}} + Π_{i j}) \cdot \nabla_{i} W_{i j} \end{matrix}

$\Pi$ 为人工粘性项，随后即可根据内能增量更新内能。

2.2.6 动量方程：速度/位移更新

将动量方程（15）离散为粒子形式：

\begin{matrix} (30) & \frac{d u_{i}}{d t} = \sum_{j = 1}^{N} m_{j} (\frac{σ_{i}}{ρ_{i}^{2}} + \frac{σ_{j}}{ρ_{j}^{2}} + R_{i j}) \cdot \nabla_{i} W_{i j} + g_{i} \end{matrix}

$\boldsymbol R_{ij}$ 是避免数值不稳定的人工应力。

该式即为粒子的加速度，采用显式时间积分，如向前欧拉、预估-校正、龙格库塔、蛙跳等格式，即可得到更新的速度和位置。

通过上述步骤，SPH 方法将连续介质力学方程转化为粒子系统的离散动力学问题，实现了固体力学行为的数值模拟。