【SPH模型入门】很适合于新手的一篇文章

转自：

https://thecodeway.com/blog/?p=83

https://thecodeway.com/blog/?p=139

https://thecodeway.com/blog/?p=161

https://thecodeway.com/blog/?p=204

SPH算法简介（一）: 数学基础

SPH(Smoothed Particle Hydrodynamics)算法是一种流体模拟算法，特点是简单快速，可以用在例如游戏这样的实时的交互软件中。SPH算法虽然简单，但要完全搞明白其中的原理和实现方法，也不是易事，写这个系列希望能全面介绍一下相关的内容，如果你搜索到这里，可以仔细看一下这个系列，希望能帮到你。
烟雾、海浪、水滴…，这些司空见怪的自然现象其实有着非常复杂的数学规律，对于流体的研究，有两种完全不同的视角，分别是欧拉视角和拉格朗日视角。欧拉视角的坐标系是固定的，如同站在河边观察河水的流动一样，用这种视角分析流体需要建立网格单元，还会涉及到有限元等复杂的工程方法，一般用在离线的应用中。而拉格朗日视角则将流体视为流动的单元，例如将一片羽毛放入风中，那么羽毛的轨迹可以帮我们指示空气的流动规律。

SPH算法是典型的拉格朗日视角，它的基本原理就是通过粒子模拟流体的运动规律，然后再转换成网格进行流体渲染。

介绍一下SPH算法涉及到的相关数学概念：

标量场和矢量场

如果空间区域内一点M，都有一个确定的数量f(M)，则称这个空间区域内确定了一个标量场，如果空间区域内任意一点M，都有一个确定的向量F(M)，则称这空间区域内确定了一个矢量场。例如，液体中的密度就是标量场，而速度就是矢量场。

偏导数

对于多元函数z=f(x,y)，定义z在(x0,y0)处相对于x的偏导数为

$\frac{\partial z}{\partial x}=\lim_{\Delta x\rightarrow 0}\frac{f(x_0+\Delta x,y_0)-f(x_0,y_0)))}{ \Delta x}$ (1.1)

例如，定义 $z=x^2+2xy+y^3$ ,那么 $\frac{\partial z}{\partial x}=2x+2y$ , $\frac{\partial z}{\partial y}=2x+3y^2$

哈密顿算子

哈密顿算子∇在流体力学中是如此重要，以至于很多地方将这个符号作为流体力学的标志，所以这里要着重介绍一下，所谓“算子”，就是那种不能单独存在，必须和其他符号放在一起的一种数学符号，例如微分中的那个“d”。哈密顿算子的定义如下：

$\nabla\equiv \overrightarrow{x} \frac{\partial }{\partial x}+\overrightarrow{y} \frac{\partial }{\partial y}+\overrightarrow{z} \frac{\partial }{\partial z}$ (1.2)

哈密顿算子有很多有趣的特性，它本身虽然并不是一个矢量，但很多运算确实可以把它视为一个矢量，例如把它作用在一个标量场A=f(x,y,z)上，那么

$\nabla A =\overrightarrow{x} \frac{\partial f}{\partial x}+\overrightarrow{y} \frac{\partial f}{\partial y}+\overrightarrow{z} \frac{\partial f}{\partial z}$ (1.3)

这个运算可以视为一个矢量和标量的乘法，得到的∇ A是一个矢量场，称为A的“梯度”，顾名思义，梯度的含义就是标量场A在某处变化快慢和方向，比如一个标量场H(x,y)是一座高山在(x,y)处的高度，则H的梯度是该高山在某处陡峭的程度，并且方向指向高处。

上面两个图中，标量场是黑白的，黑色表示大的数值，而其相应的梯度用蓝色箭头表示

而如果把哈密顿算子作用在一个矢量场 $\overrightarrow{A}$ 上，得到的∇ A称为矢量场A的“散度”，散度的计算和矢量的点积运算相似，得到的是一个标量场。

$\nabla \cdot \overrightarrow{A} =(\overrightarrow{x} \frac{\partial }{\partial x}+\overrightarrow{y} \frac{\partial }{\partial y}+\overrightarrow{z} \frac{\partial }{\partial z}) \cdot (\overrightarrow{x} A_x+\overrightarrow{y} A_y+\overrightarrow{z} A_z)=\frac{\partial A_x}{\partial x}+\frac{\partial A_y}{\partial y}+\frac{\partial A_z}{\partial z}$ (1.4)

散度的意义也很明显，就是描述一个矢量场“发散”的程度，例如下面的两个矢量场，左边的有很大的散度，而右边的散度为0

两个散度差别很大的矢量场

拉普拉辛算子

拉普拉辛算子∇2 $\nabla^2$ 是二阶微分算子，有时也可写作△或者∇⋅∇

$\nabla^2\equiv\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}+\frac{\partial^2 }{\partial z^2}$ (1.5)

例如对于A=f(x,y,z)

$\nabla^2 A\equiv\frac{\partial^2 A}{\partial x^2}+\frac{\partial^2 A }{\partial y^2}+\frac{\partial^2 A }{\partial z^2}$ (1.6)

SPH算法简介（二）: 粒子受力分析

SPH算法的基本设想，就是将连续的流体想象成一个个相互作用的微粒，这些例子相互影响，共同形成了复杂的流体运动，对于每个单独的流体微粒，依旧遵循最基本的牛顿第二定律。

$m\overrightarrow{a}=\overrightarrow{F}$ (2.1)

这是我们分析的基础，在SPH算法里，流体的质量是由流体单元的密度决定的，所以一般用密度代替质量

$\rho \overrightarrow{a}=\overrightarrow{F}$ (2.2)

这里的的作用力F的量纲发生变化，正常情况下，“力”的量纲 $dim F=MT^{-2}L$ ，而在这里 $dim F=MT^{-2}L^{-2}$ 后面的分析都是用这个量纲的“作用力”，这一点一定要注意。作用在一个微粒上的作用力由三部分组成

$\overrightarrow{F}=\overrightarrow{F}^{external}+\overrightarrow{F}^{pressure}+\overrightarrow{F}^{viscosity}$ (2.3)

其中 $\overrightarrow{F}^{external}$ 称为外部力，一般就是重力

$\overrightarrow{F}^{external}=\rho \overrightarrow{g}$ (2.4)

$\overrightarrow{F}^{pressure}$ 是由流体内部的压力差产生的作用力，试想一下在水管中流动的液体，进水口区域的压力一定会比出水口区域大，所以液体才会源源不断的流动，数值上，它等于压力场的梯度，方向由压力高的区域指向压力低的区域。

$\overrightarrow{F}^{pressure}=-\nabla p$ (2.5)

$\overrightarrow{F}^{viscosity}$ 是由粒子之间的速度差引起的，设想在流动的液体内部，快速流动的部分会施加类似于剪切力的作用力到速度慢的部分，这个力的大小跟流体的粘度系数μ以及速度差有关

$\overrightarrow{F}^{viscosity}=\mu \nabla^2\overrightarrow{u}$ (2.6)

带入公式2.2，可以得到

$\rho\overrightarrow{a}=\rho\overrightarrow{g}- \nabla p+\mu \nabla^2\overrightarrow{u}$ (2.7)

加速度形式：

$\overrightarrow{a}=\overrightarrow{g}- \frac{\nabla p}{\rho}+\frac{\mu \nabla^2\overrightarrow{u}}{\rho}$ (2.8)

实际运算过程中，有时还要考虑表面张力的影响，所谓表面张力大家应该并不陌生，肥皂泡、毛细管等有趣的物理现象都跟表面张力有关，这个力可以简单理解为流体试图减小表面而产生的力。

表面张力是由于界面层流体分子受力不均衡产生的

由于表面张力只涉及到表层的粒子，所以计算方法和上面的有所不同，这部分会在以后的章节介绍。
经过上面的分析，我们基本上搞清楚了SPH粒子的运动计算方法，下节我们将正式开始介绍SPH算法的关键部分，如何通过光滑核函数计算粒子运动规律。

SPH算法简介（三）: 光滑核函数

和其他流体力学中的数学方法类似，SPH算法同样涉及到“光滑核”的概念，可以这样理解这个概念，粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，这种随着距离而衰减的函数被称为“光滑核”函数，最大影响半径为“光滑核半径”。

光滑核函数一般具有的形态

反过来不难理解，尽管我们将流体视为一个个分散的粒子，但流体毕竟是连续充满整个空间的，流体中每个位置参与运算的值都是由周围一组粒子累加起来的。

设想流体中某点 $\overrightarrow{r}$ （此处不一定有粒子）,在光滑核半径h范围内有数个粒子，位置分别是 $\overrightarrow{r_0},\overrightarrow{r_1},\overrightarrow{r_2},...,\overrightarrow{r_j}$ ，则该处某项属性A的累加公式为：

$A(\overrightarrow{r})=\sum_{j}A_j \frac{m_j}{\rho_j}W(\overrightarrow{r}-\overrightarrow{r_j},h)$ (3.1)

其中 $A_j$ 是要累加的某种属性， $m_j$ 和 $\rho_j$ 是周围粒子的质量和密度， $\overrightarrow{r}$ 是该粒子的位置，h是光滑核半径。函数W就是光滑核函数。
光滑核函数两个重要属性，首先一定是偶函数，也就是W(−r)=W(r)，第二，是“规整函数”，也就是∫W(r)dr=1

SPH推导过程

我们假设流体中一个位置为 $\overrightarrow{r_i}$ 的点，此处的密度为 $\rho_i$ 、压力为 $p_i$ 、速度为 $\overrightarrow{u}(r_i)$ ，那么我们可以根据上一篇的公式2.8，可以推导出此处的加速度 $\overrightarrow{a}(r_i)$ 为

$\overrightarrow{a}(r_i)=\overrightarrow{g}- \frac{\nabla p(r_i)}{\rho (r_i)}+\frac{\mu \nabla^2\overrightarrow{u}(r_i)}{\rho(r_i)}$ (3.2)

对于SPH算法来说，基本流程就是这样，根据光滑核函数逐个推出流体中某点的密度，压力，速度相关的累加函数，进而推导出此处的加速度，从而模拟流体的运动趋势，下面我们逐个来分析

密度

根据公式3.1，用密度ρ代替A，可以得到

$\rho({r_i})=\sum_{j}\rho_j \frac{m_j}{\rho_j}W(\overrightarrow{r_i}-\overrightarrow{r_j},h)=\sum_{j}{m_j}W(\overrightarrow{r_i}-\overrightarrow{r_j},h)$ (3.3)

计算使用的光滑核函数称为Poly6函数，具体形式为:

其中是一个固定的系数，根据光滑核的规整属性，通过积分计算出这个系数的具体值，在2D情况下，在极坐标中计算积分：

3D情况下，在球坐标中计算:

由于所有粒子的质量相同都是m，所以在3D情况下， $\overrightarrow{r_i}$ 处的密度计算公式最终为：

压力

根据上一节的结论，在位置 $r_i$ 之处的由压力产生的作用力的计算公式为

不过不幸的是，这个公式是“不平衡”的，也就是说，位于不同压强区的两个粒子之间的作用力不等，所以计算中一般使用双方粒子压强的算术平均值代替单个粒子的压力 $r_i$ 之处的由压力产生的作用力的计算公式为

对于单个粒子产生的压力p,可以用理想气体状态方程计算

其中 $\rho_0$ 是流体的静态密度， $K$ 是和流体相关的常数，只跟温度相关。
压力计算中使用的光滑核函数称为Spiky函数

在3D情况下，

将公式3.12带入3.9，可以整理出公式3.2中压力产生的加速度部分

粘度

SPH算法简介（四）:Hello,SPH

上几节，我们推导出一大推复杂无比的公式，似乎有点纸上谈兵，这节来点真的，写一个可以运行的SPH系统，下面就是SPH基本的运算流程

初始化粒子，为每个粒子赋上初始位置
根据公式3.7计算每个粒子的密度
根据公式3.10计算每个粒子的压强
根据公式3.18计算每个粒子的加速度
根据临界条件调整加速度
根据加速度计算每个粒子的速度变化
根据速度计算粒子位置的变化
绘制粒子
回到步骤2

下面有个简单的示例程序，运行效果如下

        这个程序基本上没有怎么考虑效率，只是让系统跑起来，所以比较适合拿来对照公式学习，按照惯例，放出源代码和可执行程序
        源码下载:fluid_source.zip(395KB) [https://thecodeway.com/blog/?p=204]
        可执行程序下载: fluid.zip(120KB)
        SPH还有很多细节值得讨论，比如表面张力、并行计算、构建网格、真实材质的水渲染等，这些部分我会抽时间再写一些东西出来介绍。