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

2单流体SPH算法实现

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

2.1算法框架

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

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

Created with Raphaël 2.1.2 初始化粒子属性计算粒子密度ρ 计算粒子压强p 计算粒子受到的压力F_pressure 计算粒子的粘滞力F_viscosity 得到粒子的加速度更新粒子的位置

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

2.2 寻找邻域粒子

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

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

我们可以通过划分网格的方式搜索邻域粒子，如下图所示。以 h h <script type="math/tex" id="MathJax-Element-38">h</script>为网格长度划分网格，使得网格覆盖所有的粒子。这样搜索某一个粒子的邻域粒子时，只需要搜索器周围的9个网格，而无需遍历所有的粒子，大大增加的效率。

2.3算法实现

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

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

下面给出源代码：SPH

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