光滑粒子流体动力学_基于SPH（光滑粒子流体动力学）算法的流体仿真

实现的效果很粗糙，没有添加渲染和表面绘制，现在只做到了点和线的程度。这个月考试比较多，所以做的时间也没有很多。

先放一下最终的效果图。

碰撞检测的粒子碰撞后速度的计算还有有些问题的，碰撞检测做的比较简单。有时间了再单独修改碰撞检测方面的内容，先把流体仿真实现了再说。

没有表面构建，也看不出来具体是不是有什么问题（但是粗略的看除了上边几个粒子不动了，也没有什么太大的问题，可能就是碰撞时候出的问题？有空再多看看相关的资料），构建表面时候还是会整体进行修改的。

在介绍SPH算法之前，首先我们先了解一些数学原理。

偏导数

若z=f(x,y)，则z对x的偏导数为：

同理对y的偏导数只要将x改为y即可。

哈密顿算子

哈密顿算子∇在流体力学中非常重要，这里介绍一下哈密顿算子，所谓“算子”，就是那种不能单独存在，必须和其他符号放在一起的一种数学符号，例如微分中的那个“d”。哈密顿算子的定义如下：

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

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

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

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

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

拉普拉辛算子

拉普拉辛算子∇2是二阶微分算子，有时也可写作Δ或者∇∙∇

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

接下来是粒子的受力分析。

SPH算法的基本思想就是将连续的流体看成许多个相互作用的粒子，通过粒子之间的相互作用形成了复杂的流体运动。每个粒子都有自己的基本属性，包括位置、速度、质量、所受的重力、压力、粘力等等。并且除质量外都要每一帧进行重复计算。

对于每个粒子，都遵循基本的牛顿第二定律。

F = ma；

在SPH算法里，流体的质量是由流体单元的密度决定的，所以一般用密度代替质量：

在这里，粒子所受的力我们只研究三个，分别是重力、压力以及粘力。

试想一下在水管中流动的液体，进水口区域的压力一定会比出水口区域大，所以液体才会源源不断的流动，数值上，它等于压力场的梯度，方向由压力高的区域指向压力低的区域。

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

将上述三个力结合起来，就得到

加速度则为

光滑核函数：

sph算法中每个粒子的光滑核函数我们可以理解为：在光滑核半径内，粒子受到的力受到粒子之间距离的影响，距离越近受到的影响越大。具体形式如下图：

我们将流体看成多个粒子的集合，每个粒子都受到周围粒子的影响，

我们看系统中的某个粒子，他到临近的粒子的距离分别为r1，r2......rn，则所观察粒子受到的影响即为每个粒子带来的影响之和。这个属性可以是任何需要叠加的属性，比如密度、压力、粘度等。

其中Aj是要累加的某种属性（如密度，压力，粘度）。其中m是粒子的质量，ρ是粒子的密度。r是粒子的位置，h是光滑核半径，W就是光滑核函数。

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

密度

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

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

其中

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

由于所有粒子的质量相同都是m，所以在3D情况下，所研究粒子的密度为：

压力

同样，我们用p代替A，就会求得我们所需要的压力公式。

由于位于不同压强区的两个粒子之间的作用力不等，所以计算中一般使用双方粒子压强的算术平均值代替单个粒子由压力产生的作用力的计算公式为：

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

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

在3D情况下，

Kspiky = 15/（πh6），即

由此我们可以整理出公式中压力产生的加速度部分：

粘度

公式最后一部分，由粘度产生的作用力：

这个公式同样有“不平衡”的问题，考虑到公式中的速度其实并不是绝对速度，而是粒子间的相对速度，所以这个公式的正确写法应该是：

其中的光滑核函数形式如下：

在3D情况下，Kviscosity = 15/(2πh3）,由此我们可以得到：

由此可得到的粘度部分加速度：

到这里，根据我们已经获得的数据，我们就可以计算最后的加速度了：

我们所需要的理论，都已经介绍完了，在搞懂这些之后就可以编写代码，实现我们自己的流体仿真了。

具体的代码我已经上传到我的git上，项目地址：

wangfeng70117/fluid_simulation​github.com