2单流体SPH算法实现
经过前一章的介绍,知道了SPH算法的原理,这一章我们介绍SPH算法的代码具体实现
2.1算法框架
SPH算法的思想是用粒子来模拟流体,其中粒子承载了各种属性(如 位置、速度、加速度、密度、压强等),通过不断更新粒子的属性(其实最终的目的就是为了更新粒子的位置),达到流体动态模拟的效果。而粒子的属性受其周围粒子的影响,一般通过光滑核函数 ,来计算粒子的各种属性。
所以,SPH算法迭代的流程如下:
算法框架如下:
1. 初始化粒子的位置、速度、加速度、静态密度。
2. 根据光滑核函数计算粒子的插值密度
3. 根据公式
计算粒子的压强
4. 根据光滑核函数以及以及前二步得到的密度
和压强
,计算粒子的压力
5. 根据光滑核函数以及粒子的粘度系数
和粒子的速度
计算粒子的粘滞力
6. 累加压力
、粘滞力
和重力
得到粒子的受力
,除于密度
,进而得到粒子的加速度
7. 更新粒子的速度
、位置
8. 返回第二步,一直迭代。
2.2 寻找邻域粒子
我们知道,粒子的各种属性是由光滑核半径
之内的邻域粒子的属性计算得到的,那如何找到这些邻域粒子呢?
最常规的思路是,计算两个粒子之间的距离
,如果
,则代表两个粒子相互影响。
然而我们发现,用这种方式去搜索所有粒子的邻域粒子的话,其时间复杂度为
,可见这种方式效率不高。那如何改进呢?
我们可以通过划分网格的方式搜索邻域粒子,如下图所示。以 为网格长度划分网格,使得网格覆盖所有的粒子。这样搜索某一个粒子的邻域粒子时,只需要搜索器周围的9个网格,而无需遍历所有的粒子,大大增加的效率。
2.3算法实现
写了一个简单的示例程序,运行效果如下:
该demo中粒子的绘制是用OSG写的。如果只是为了研究SPH算法,其实不需要太关注粒子绘制的问题,OSG绘制粒子只需要给SPH提供一个函数接口——输入是所有粒子的坐标,输出就是OSG用自己的渲染线程把粒子绘制出来。所以,每次更新完粒子的位置都调用一次这个函数就可以了。
下面给出源代码:SPH
程序是用vs2010开发的,OSG库的版本是3.4.1,要跑起来得要配置一下OSG。另外,SPH算法都写在sph.h和sph.cpp这两个文件中。