SPH真实感流体交互模拟的改进算法

程志宇 ，徐国庆*，张岚斌 ，许 犇

1.武汉工程大学计算机科学与工程学院，湖北 武汉 430205；

2.智能机器人湖北省重点实验室（武汉工程大学），湖北 武汉 430205

真实感流体模拟是计算机图像学领域中的研究热点，不管是在电视特效和广告、电子游戏亦或是军用作战平台仿真、虚拟医学等领域中都有着广泛的应用。基于物理的流体模拟的计算模型应该包括流体每一时刻的运动方向、受力情况以及在运动过程中流体间的相互作用力等，同时为了更好的提高真实感要增强细节模拟［1-2］，展现实时流动的视觉效果。

现阶段，对于流体模拟主要采用欧拉网格法和拉格朗日粒子法两种方法来求解计算流体动力学中的控制方程。欧拉法是模拟流体在固定网格单元上的运动，将流体离散成空间中的固定点，空间点的物理属性包括密度、速度、压强等，通过固定点的属性变化，来模拟流体的运动过程，基于网格的方法可以非常有效地模拟蒸汽和烟雾这样的体积效应［3-4］。而基于粒子的拉格朗日法则是跟踪单个流体粒子的运动，并记录其物理属性变化，来完成对流体的模拟。拉格朗日粒子法在保持自由表面的质量或体积方面有明显优势。

基于上述分析，针对真实感流体具有的复杂动力学特性，提出一种在基于通用光滑粒子流体动力学（smoothed particle hydrodynamics，SPH）方法上建立的流体模型，利用粒子模拟流体运动过程，SPH初始化完成粒子创建，相邻粒子根据划分的网格单元进行交互。根据粒子的空间位置构建流体密度场计算粒子密度、压力、粘滞力、外力等，采用固定镜像粒子方法［5］处理交互边界问题［6-7］。通过移动立方体（marching cubes，MC）算法来重建流体网格［8-10］，实现对流体表面建模，光照计算实现反射和折射效果，完成对流体的高质量交互式渲染［11］，形成更逼真的动态演示效果。

1 流体建模算法

SPH方法是模拟流体流动的一种拉格朗日形式的无网格粒子法。它的理论基础是粒子方法，由于粒子的运动与流体的运动相似，可以用经典牛顿流体动力学的控制方程来模拟。与传统的流体模拟方法相比，采用SPH方法能够比较真实的模拟流体运动过程，对于边界不连续问题，大面积变形问题，模拟破裂合并以及飞溅等，SPH方法可以很好的解决。

它的核心思想是用一系列任意分布的粒子来表示问题域，粒子之间不需要任何连接，用积分形式来近似场函数，再通过局部区域（支持域）的相邻粒子所携带的计算数值叠加求和来取代场函数及其导数的积分表示式，实现进一步近似。在每一个时间步中都要进行粒子近似，因此当前粒子的性质取决于当前局部分布的粒子。

1.1 流体运动方程

基于拉格朗日的流体运动方程可由连续性方程（1）、运动方程（2）、压力状态方程（3）表示［12］。

其中：v为速度，p为压力，ρ为密度，μ为动力学黏性系数为流体高度，g为重力加速度，可引入人工黏性项∏来修正运动方程。

式（4）中：α是0到1之间的常系数，v是流体速度，r是粒子的位置矢量，η是保证分母不为零的参数，h是光滑核半径。显然，拉格朗日方法是针对于相对独立的各个流体粒子来求解运动方程，通过积分计算流体粒子下一时刻的物理属性值变化。

1.2 SPH基本方程

在SPH方法中，问题域是用具有质量和体积大小的有限粒子表示，通过对粒子物理属性插值计算区域任一点的对应值：

其中：A（r）是位置r的物理属性；N为粒子j在支持域中的粒子总数；mj是粒子j的质量；ρj是粒子j的密度；函数W为光滑核函数，h为光滑核半径，光滑核函数具有2个性质：首先必须是偶函数，满足W（-r）=W（r）；其次，必须为规整函数，满足

在SPH方法中，导数只影响光滑核。光滑属性函数A(r)的梯度以及拉普拉斯算子：

2 改进算法

2.1 粒子建模

流体方程计算时，首先在密度场应用式（5）计算粒子的求和密度ρ，使用Poly6核函数估计：

将式（6）和式（7）代入运动方程并进行对称化，得到作用在粒子上的压力和黏度力位于不同压强区的粒子间的相互作用力不等，所以采用相互作用粒子压强的算术平均值来优化计算式。

其中，基于计算精度和效率，估计计算选择核函数形式如下：

最后计算得出流体的加速度a，进行速度场的更新。

模拟流体动态效果的SPH流体粒子模型算法的基本流程如图1所示。

图1 SPH粒子系统建模流程Fig.1 Modeling flowchart of SPH particle system

此SPH模型适合用来求解高速碰撞过程中的流体大变形的问题，通过不断循环迭代修正密度差［13-14］，对流体的形状和运动状态都能模拟较真实效果，模拟具有有效性和高效性。

2.2 自由表面建模和渲染

流体表面重建是流体模拟的重要内容。在未实现自由界面的重构和渲染时，用SPH方法得到的是由一堆离散粒子组成的流体系统，流体的真实感以及运动连贯性欠缺。采用一种基于移动立方体（MC）算法来提取流体运动表面的真实效果，实现流体表面建模。MC算法的核心思想是在三维离散数据场中按一定规则划分体元，在每个体元内部根据设定值对每条边进行插值，通过线性插值得到的三角面片来逼近等值面，完成等值面提取。在动画模拟过程中图形处理器（graphics processing unit，GPU）加速进行实时渲染［15］，将流体表面渲染成透明的，并进行光照计算，实现了反射和折射效果。同时流体表面的光照和纹理渲染提高了流体场景的真实感。

3 结果与讨论

3.1 实验结果

实验平台是 Intel（R）Core（TM）i7-8550U CPU 1.8 GHZ，主存 16 GB，图形卡为 NVIDIA GeForce MX150，显卡内存为2 GB，开发环境Visual Studio 2017，图形的编程接口由OpenGL提供。在图形硬件上，实验完成了对复杂流体方程的求解，同时针对于交互结果进行了实时渲染展示。实验结果显示了本文方法模拟的有效性和高效性，实时场景渲染的真实感有所提高。

图2（a）显示了粒子系统初始形态，粒子数约为4 000，室外场景为背景，图2（b）是渲染后的流体系统，流体采用透明渲染，用MC算法构造了自由移动的流体表面的效果。其中每个粒子包含重力、压力、粘滞力、交互粒子间的引力与斥力作用。场景包括粒子的物理属性值计算、自由表面建模、光照计算和材质渲染等。

图2 SPH模型初始状态：（a）粒子模型，（b）流体模型Fig.2 Initial state of SPH models：（a）particle model，（b）fluid model

图3（a）展示了球体渲染模型在外部输入力的作用下运动的效果，粒子数约为4 000。图3（b）是透明渲染流体在外力作用下的运动效果。在模拟过程中，用户可以通过鼠标拖拽来调整流体的受力，其中每个粒子除了受到基础作用力外还受到外部输入的力的作用。在外力作用下流体表面发生剧烈变形，形成飞溅效果。

图3 SPH运动模型：（a）球体模型，（b）流体模型Fig.3 Models of SPH motion：（a）spherical model，（b）fluid model

实验将木块运动与实时流体模拟相结合，通过计算木块下落位置和速度来触发与流体的动态交互，使流体产生相应形变。图4是流体和木块交互的动态演示效果，包括木块下沉和上浮过程、场景包括相邻网格粒子搜索和复杂的边界交互计算。

图4 流体和木块交互模型：（a）木块下沉，（b）木块上浮Fig.4 Models of fluid and wood interaction：（a）wood sinking，（b）wood floating

3.2 对比分析

小尺度的流体模型渲染速度会受到粒子数和光滑半径的影响，实验过程中通过设置对照组进行对比实验，不同粒子数和光滑半径下的渲染速度关系如表1所示。实验结果表明光滑半径相同，随着粒子数增加，渲染速度变慢；粒子数不变，随着光滑半径的增加，渲染速度变慢。

表1 不同粒子数和光滑半径下的渲染帧率Tab.1 Rendering frame rates of different particle numbers and smoothing radii frames/s

4 结 语

本文主要讨论真实感流体的模拟问题，利用基于SPH的数值计算法进行粒子建模，再采用MC算法对粒子表面建模。实验结果表明，本文提出的对于真实感流体的模拟改进算法，综合利用体绘制和面绘制法实现流体三维模拟，真实感提高，优化效果明显。但对于数值方法的固有计算损耗、细节特征模拟以及硬件加速方面要做进一步改进，提高流体模拟的效率和精度。