SPH算法的理论和实践（2）

2单流体SPH算法实现

经过前一章的介绍，知道了SPH算法的原理，这一章我们介绍SPH算法的代码具体实现

2.1算法框架

SPH算法的思想是用粒子来模拟流体，其中粒子承载了各种属性（如位置、速度、加速度、密度、压强等），通过不断更新粒子的属性（其实最终的目的就是为了更新粒子的位置），达到流体动态模拟的效果。而粒子的属性受其周围粒子的影响，一般通过光滑核函数 $A (r_{i}) = \sum A_{j} \frac{m_{j}}{ρ_{j}} W (\vec{r} - {\vec{r}}_{j}, h)$ ，来计算粒子的各种属性。

所以，SPH算法迭代的流程如下：

算法框架如下：
1. 初始化粒子的位置、速度、加速度、静态密度。
2. 根据光滑核函数计算粒子的插值密度 $ρ (r_{i})$
3. 根据公式 $p_{i} = k (ρ (r_{i}) - ρ_{0})$ 计算粒子的压强
4. 根据光滑核函数以及以及前二步得到的密度 $ρ (r_{i})$ 和压强 $p_{i}$ ，计算粒子的压力 $F_{i}^{p r e s s u r e}$
5. 根据光滑核函数以及粒子的粘度系数 $μ_{i}$ 和粒子的速度 ${\vec{u}}_{i}$ 计算粒子的粘滞力 $F_{i}^{v i s c o s i t y}$
6. 累加压力 $F_{i}^{p r e s s u r e}$ 、粘滞力 $F_{i}^{v i s c o s i t y}$ 和重力 $F_{i}^{e x t e r n a l}$ 得到粒子的受力 $F_{i}$ ，除于密度 $ρ (r_{i})$ ，进而得到粒子的加速度 $\vec{a}$
7. 更新粒子的速度 ${\vec{u}}_{i}$ 、位置 ${\vec{r}}_{i}$
8. 返回第二步，一直迭代。

2.2 寻找邻域粒子

我们知道，粒子的各种属性是由光滑核半径 $h$ 之内的邻域粒子的属性计算得到的，那如何找到这些邻域粒子呢？

最常规的思路是，计算两个粒子之间的距离 $d = {\vec{r}}_{i} - {\vec{r}}_{j}$ ，如果 $d < h$ ，则代表两个粒子相互影响。
然而我们发现，用这种方式去搜索所有粒子的邻域粒子的话，其时间复杂度为 $O (n^{2})$ ，可见这种方式效率不高。那如何改进呢？

我们可以通过划分网格的方式搜索邻域粒子，如下图所示。以 $h$ 为网格长度划分网格，使得网格覆盖所有的粒子。这样搜索某一个粒子的邻域粒子时，只需要搜索器周围的9个网格，而无需遍历所有的粒子，大大增加的效率。

2.3算法实现

写了一个简单的示例程序，运行效果如下：

该demo中粒子的绘制是用OSG写的。如果只是为了研究SPH算法，其实不需要太关注粒子绘制的问题，OSG绘制粒子只需要给SPH提供一个函数接口——输入是所有粒子的坐标，输出就是OSG用自己的渲染线程把粒子绘制出来。所以，每次更新完粒子的位置都调用一次这个函数就可以了。

下面给出源代码：SPH

程序是用vs2010开发的，OSG库的版本是3.4.1，要跑起来得要配置一下OSG。另外，SPH算法都写在sph.h和sph.cpp这两个文件中。