一种拉格朗日粒子流体的高效表面重建方法

邵绪强， 刘 艳， 王新颖

(1. 华北电力大学控制与计算机工程学院，河北 保定 071003；2. 河北金融学院国际教育学院，河北 保定 071003)

一种拉格朗日粒子流体的高效表面重建方法

邵绪强1， 刘 艳2， 王新颖1

(1. 华北电力大学控制与计算机工程学院，河北 保定 071003；2. 河北金融学院国际教育学院，河北 保定 071003)

为了在几何表示层面捕获拉格朗日粒子流体的表面细节特征以进行真实感绘制，提出一种高效的基于八叉树的自适应流体表面重建方法。首先进行流体表面粒子的精确检测；然后根据结点中是否包含流体表面粒子，对流体粒子的包围盒进行基于八叉树的自适应剖分；最后只在包含流体表面粒子的叶子结点里建立隐式距离场。该方法在流体表面重建过程中的内存消耗和计算复杂度只取决于流体表面而不是流体体积，适合大规模粒子流体场景的表面提取和绘制。

拉格朗日粒子；表面重建；流体模拟；真实感绘制

波涛汹涌的大海、倾盆而下的暴雨、瞬间就要吞噬战舰的巨大漩涡等流体影视特效给好莱坞特效大片的观众带来了震撼的视觉享受。在欣赏这些精彩特效大片的同时，人们领略到了基于物理的流体模拟技术的高超魅力。由于真实物理模型的支撑，基于物理的流体模拟技术能逼真地模拟各种复杂流体现象，给人以震撼的视觉享受，因此吸引了国内外图形学研究学者们对其进行深入研究。

基于物理的流体模拟方法[1]主要分为欧拉网格方法(Eulerian grid method)和拉格朗日粒子方法(Lagrangian particle method)。欧拉网格方法把流体占据的计算空间按一定规则剖分为网格，然后分析被流体所充满的空间中每一个网格位置上流体的速度、压强、密度等参数随时间的变化。该方法可以较容易地保证流体的不可压缩特性，但由于模拟过程中要同时把固体边界网格离散化，不易于进行流固交互模拟，并且会丢失精度小于网

格分辨率的流体细节特征。而拉格朗日粒子方法从分析流体各个微团的运动着手，即研究每个流体微团的速度、压强、密度等描述流体运动的参数随时间的变化状况，这是一种以“实体”形式为流体建模的方法，建模的基本元素就是携带流场变量的粒子。与欧拉网格方法相比，拉格朗日粒子方法更直观，自动保证数值求解过程中的质量守恒，尤其易于进行固体边界处理和模拟流体细节。因此，在电子游戏和电影特效等实际流体模拟应用中，拉格朗日粒子方法变得越来越普遍。

为了得到逼真的流体动画，在通过数值求解流体运动方程而用粒子描述流体运动信息后，拉格朗日粒子方法需要重建出光滑完整的隐式流体表面以进行真实感绘制。拉格朗日粒子流体的隐式表面重建方法可分为密度场方法和距离场方法。对于密度场方法，Blinn[2]直接用一个连续性表面包裹所有流体粒子而形成密度云，并且这个表面随流体粒子的运动而改变。Müller等[3]把该方法引入到光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)流体中，并利用球形核函数重新定义密度场，使得表面重建结果更稳定，但存在局部区域出现鼓包的现象，从而不能重建出平滑的流体表面。为了使密度场方法能够重建出足够光滑的流体表面，Yu和Turk[4]基于粒子分布的主成分分析而用椭圆形核函数来插值计算密度场；为了解决密度场方法受限于粒子采样的均匀性，Zhu和 Bridson[5]提出了基于距离场的表面重建方法，通过直接计算流体粒子到表面的距离来表示隐式流体表面，但该方法在高曲率区域存在明显的瑕疵。Solenthaler等[6]通过动态调整粒子的半径使质心移动速度小于粒子半径衰减速度，从而改善了距离场方法在高曲率区域存在的瑕疵现象。然而，上述流体表面重建方法都是针对于流体体积进行的，当对大规模拉格朗日粒子流体场景进行表面重建时就会由于过高的内存和时间消耗而失败。针对此问题，本文给出了一种高效的基于八叉树的自适应流体表面重建方法。该方法首先进行流体表面粒子的精确检测；然后根据结点中是否包含流体表面粒子，对流体粒子的包围盒进行基于八叉树的自适应剖分；最后只在包含流体表面粒子的叶子结点中建立隐式距离场。因此，流体表面重建过程中的内存消耗和计算复杂度只取决于流体表面而不是流体体积，非常适合于大规模粒子流体场景的表面提取和绘制。

1 拉格朗日粒子流体模拟

1.1 不可压缩粒子流体模拟

确保不可压缩性的计算是SPH流体模拟最困难的部分，也是其关键研究的内容。微可压缩SPH流体[7](weakly compressible SPH，WCSPH)必须使用较小的时间步长以确保不可压缩性，严重影响计算效率。通过泊松方程计算压强的不可压缩方法[8-9](incompressible SPH，ISPH)可以使用较大时间步长，但在每个时间步内由于迭代求解线性方程组而导致计算量庞大。Ihmsen等[10]结合对称的SPH压力和连续方程离散方法来离散化压强泊松方程，并对压力进行实际计算，提高了求解的收敛速度。

为了能够在模拟过程中使用较大时间步长，本文采用Solenthaler和Pajarola[11]提出的预测修正不 可 压 缩 SPH 方 法 (predictive-corrective incompressible SPH，PCISPH)进行不可压缩流体模拟。PCISPH结合了WCSPH和ISPH的优点，允许使用较大时间步长，并简化了每个时间步内的计算复杂度。算法的核心思路是通过预测修正密度误差的方法迭代计算流体粒子的压强和压力，压力的改变会影响粒子的速度和位置，然后位置的相对变化会影响粒子密度，最后粒子的相对位置和密度影响压力大小和方向，在这样的迭代作用的循环中找到满足约束条件的压强为止。PCISPH的伪代码如图1所示。

图1 PICSPH伪代码

图1给出了PCISPH在每个时间步内的迭代计算过程。在迭代计算过程中，首先对每个粒子搜索其核函数紧支域内的邻居粒子信息。然后计算出除压力外的所有受力，包括粘性力、重力、表面张力、甚至漩涡力。进入预测修正循环之前，粒子的压强和压力初始化为 0。最后便进入PCISPH的压强计算的循环中。每次循环都要检查粒子的密度误差的大小是否超过阈值η，其中密度误差为。如果，则使用当前粒子的合力来预测速度和位置，根据预测出来的粒子速度和位置重新计算粒子的密度和密度误差，并将所有粒子的最大误差保存在密度误差中。在计算过程中，根据密度误差对流体的压强值进行累加

在预测下一步的密度时，应该根据预测的位置信息更新邻居粒子列表。考虑到计算效率，这里不再重新搜索邻居粒子，而是通过更新到现有邻居粒子距离的方式来近似。虽然这样的估计会带来误差，但是在最少迭代次数、时间步长、核函数影响范围衰减等约束条件的作用下，对整个流体运动模拟结果影响是可以忽略的。采用一个最小迭代次数的方法可显著减少误差的局部震荡，而且这给予粒子足够的时间来扩散预测信息。为了提高计算效率，一般使用最少迭代3次，而3次迭代已经足够将粒子密度误差约束到目标范围内。

从PCISPH的执行过程来看，邻居搜索是其最大的计算瓶颈，为了提高该部分的计算性能，可建立 KD-Tree以实现邻居例子的快速搜索。与WCSPH相比，PCISPH由于增加了迭代计算来计算压力而大幅度降低了性能，但是采用的时间步长却增加了30倍左右。

1.2 流固交互模拟

对于SPH粒子流体和运动固体边界的交互，本文结合文献[12]的惩罚力交互方法和文献[13]的自适应边界粒子采样方法实现双向流固交互模拟，计算过程分为3步：

(1) 在每个时间步内进行流固交互计算之前，根据初始流体粒子间距r0对任一表面三角形进行自适应边界粒子采样，如图 2所示。首先选取三角形顶点为边界粒子；然后以初始流体粒子间距r0对三角形的边进行边界粒子采样；确定三角形最短边es和最长边el，并选取向内的法线方向为扫描线方向对三角形内部进行边界粒子采样。自适应边界粒子采样可以更精确地表示固体表面的拓扑结构，从而实现更准确的固流交互。

图2 自适应收缩压粒子采样

(2) 考虑边界粒子的相对密度贡献[12-13]，计算任一流体粒子i的密度

其中，j为邻居流体粒子索引；k为属于固体边界的邻居粒子的索引；W为 SPH方法的核函数；为相对贡献度量函数； Vbk为第k个边界邻居粒子当前体积。式(2)通过φ() x度量边界粒子的相对贡献，解决了其非均匀分布特性。

(3) 对于一对相邻的流体粒子i和固体粒子j，计算 j对 i的包括法向力和切向力在内的交互力。利用 SPH方法求解流体方程得到的压力、粘滞力来分别近似法向力、切向力。并且，对i的作用力同样采用φ() x进行度量

其中，μ为物体表面磨擦系数；v为速度。

2 粒子流体表面重建

2.1 流体表面粒子检测

流体的表面可以由位于流体表面的粒子表示和跟踪。基于计算流体动力学(computational fluid dynamics，CFD)领域文献[14]方法，给出一个快速准确的流体表面粒子检测方法。对于任一流体粒子i，首先计算一个再归一化矩阵Bi

其中，Vj是第 j个邻居粒子的体积；W为再归一化的高斯核函数[15]； 40h= d，d0为初始流体粒子间距。

然后，计算矩阵Bi的最小特征值，并根据Bi和实验验证值的大小关系判定粒子 i是否属于SPH流体表面S

由式(6)可得到一个包含了SPH流体表面粒子的粒子集合。与文献[14]方法相比，本文方法通过设定miniλ 的阈值为0.75，把自由流体表面附近的流体粒子都确定为流体表面粒子集合，虽然引入了一些流体内部的粒子，但提高了检测准确度。

2.2 自适应距离场重建

为了降低粒子流体表面重建过程中的内存消耗和计算复杂度，基于文献[16]中流体表面粒子检测和隐式表面定义方法，给出一种基于八叉树[17]的自适应距离场构建算法，该方法只在包含流体表面粒子的叶子结点上计算距离场值。

在自适应距离场构建算法中，任一八叉树结点包含两个粒子集合：流体表面粒子集合和流体粒子集合。其中，流体表面粒子集合包含所有属于该结点的流体表面粒子，流体粒子集合包含所有属于该结点的流体粒子。本方法把八叉树的结点分为 3类：外部结点、内部结点和表面结点。对于外部结点来说，其粒子集合流体表面粒子集合和流体粒子集合均为空；内部结点的粒子集合流体表面粒子集合为空，而粒子集合流体粒子集合非空；而对于表面结点，其粒子集合流体表面粒子集合非空。八叉树距离场的递归创建过程如图3所示。

图3 基于八叉树的自适应距离场的建立

从图 3中左图可以看出，基于八叉树的自适应距离场的创建过程分为 3步：①先为整个流体模拟建立一个根结点。该根结点的物理空间为整个流体的最小立方体包围盒，流体表面粒子集合包含所有的流体表面粒子，流体粒子集合包含所有流体粒子。②从根结点开始，对每个八叉树结点进行递归剖分，直到其粒子集合流体表面粒子集合为空或者达到最大深度。③对于每个叶子结点，如果其流体表面粒子集合流体表面粒子集合非空，则根据下式计算其8个顶点的距离场值

半径； R= 4r； f∈ [0… 1 ]并通过下式进行计算

其中，EVmax为矩阵的最大特征值；γ通过文献[17]方法进行计算。图3右图给出了如何确定八叉树结点包含的粒子集合以避免粒子流体表面重建过程中产生的空洞。当为八叉树的任一结点建立粒子集合流体表面粒子集合和流体粒子集合时，本方法考虑落在比该结点物理空间大l的AABB包围盒内的所有流体表面粒子和流体粒子。在本方法中，l等于式(7)中的R。

为了保证基于八叉树的流体表面重建算法能够重建出仅由一层粒子构成的薄膜结构的表面，可给出了八叉树叶子结点边长和粒子半径的大小关系以避免表面重建失败，如图 4所示。在图 4中，用二维示意图的形式分析了八叉树叶子结点边长L和粒子半径r的大小关系，其中r值对于确定流体场景是固定的。图 4(a)给出的特例表明，

图4 粒子半径r与叶子结点边长L的大小关系

3 实例验证

将通过一系列 3D实验证明本文提出的针对拉格朗日粒子流体的表面重建算法。实验硬件配置为：Intel Xeon E5630 CPU，6 GB RAM 和NVIDIA Geforce GTX 680 GPU。编程环境为：Windows 7，Visual Studio 2008，OpenGL 2.0。在3D场景中，首先利用SPH粒子方法模拟了流体运动，然后利用本文给出的自适应流体表面重建方法创建流体的隐式表面，最后利用Marching Cubes算法[3]提取三角形网格面并利用开源渲染引擎Pov-Ray进行离线绘制。实验中用到的物理量以及取值如表1所示。

表1 漩涡模拟中物理参数及取值

3.1 表面重建效果分析

基于八叉树的自适应流体表面重建算法可以高效地重建出任一粒子集合的隐式表面，并捕获漩涡细节的几何结构以进行真实感绘制。该算法获取流体表面的过程如图5所示。在图5中利用基于八叉树的自适应流体表面重建算法对犰狳的采样粒子集合进行了表面重建。对于输入的犰狳网格模型，首先对其所占的空间进行规则粒子采样得到12.1 k粒子(图5(a))，采样距离为0.018 m。接着利用本文给出的表面粒子检测方法检测出粒子集合的表面粒子；图5(b)给出了利用本方法进行表面粒子检测得到的结果，表面粒子数为3.5 k。然后，根据检测到的流体表面粒子，建立八叉树自适应距离场(图 5(c))。最后图 5(d)给出了Marching Cubes算法提取的三角形网格面的绘制结果。

图5 自适应流体表面重建算法

图 6给出了初始形状为犰狳模型的流体在重力作用下的运动模拟结果。通过使用本文提出的粒子流体表面重建方法建立距离场，并利用开源绘制引擎Pov-Ray对Marching Cubes算法提取的三角形网格面进行真实感绘制。从绘制结果可以看出，本文的方法可以光滑地重建出具有复杂形状的粒子流体表面，清晰地表现飞溅液滴等细小流体特征。

图6 犰狳模型流体的表面重建

图7给出了刚性圆柱体以10 m/s的速度在水面移动产生的漩涡细节模拟结果。从模拟结果可以看出，本文的表面重建算法可以重建出刚性圆柱体后面产生的丰富漩涡细节特征，增强了流固交互模拟的真实感。

图 8给出了多个具有不同密度的刚性物体和水坝的双向交互模拟效果。从效果图可以看出，本文的表面重建算法可以重建出固流交出中的飞溅液滴这一细节特征，有效提高了流固交互模拟的真实感。

图7 刚性圆柱和流体单相交互的流体表面重建

图8 水坝冲击不同密度物体的流体表面重建

3.2 表面重建性能分析

使用 3种不同方法：CAVW_2007[6]、CGF_2012[18]和本文方法对不同粒子间距采样后的犰狳模型进行了表面重建，并比较了内存消耗。对比结果如表 2所示。表中给出了在不同粒子半径 r和 Cube大小下不同重建方法的内存消耗。内存消耗通过在重建过程中需要计算距离场值的立方体个数(#Cube)来进行衡量。CAVW_2007方法对流体空间进行规则均匀划分，在每一个Cube里都要计算距离场，因此内存消耗是最大的；CGF_2012方法则只在流体表面附近的 Cube里计算距离场，大幅度降低了内存消耗，但由于不能准确地检测出流体表面粒子和设定r和L的大小关系，所以不能重建出水膜的表面；针对本文方法同样大幅度降低了内存消耗，虽然在八叉树的自适应创建和访问引入了部分计算开销，但整体计算性能与 CGF_2012方法持平，且能重建出由单层粒子组成的水膜结构的表面。

表2 流体表面重建算法的内存消耗对比

本文利用基于八叉树的自适应表面重建方法对所有三维场景进行流体表面重建，并统计了性能数据，如表 3所示。在表 3中，#FP、#SP和#FPsurface分别代表流体粒子数、固体粒子数以及流体表面粒子数。#Cube3和#Cube1分别代表本方法和 CAVW_2007方法在流体表面重建过程中需要计算距离场值的Cube数量。表示本文方法与CAVW_2007方法的运行时间比值。#Triangles表示利用Marching Cubes算法提取三角形网格面后得到的三角形数量。从表 3给出的统计数据可以看出，在三维流体模拟过程中，流体表面粒子数量通常只占流体粒子总数的很小一部分(＜20%)。

与 CAVW_2007方法相比，基于八叉数的自适应表面重建算法在大幅度降低内存消耗的同时，显著提高了粒子流体表面重建的计算效率，适合大粒子规模下的流体表面的重建。

表3 基于八叉树的自适应表面重建算法的内存消耗

本文的流体表面重建算法只在整个流体模拟的计算过程中引入了较小的额外计算量。还给出了上述三维场景模拟过程中每个时间步内的平均计算时间，并对本文的表面重建方法进行时间性能分析。时间统计包括：物理模拟计算时间Tps和基于八叉树的自适应表面重建时间 Tode，统计结果如表4所示。从表4数据可以看出，模拟过程中主要的计算瓶颈是复杂的物理模拟部分，占据了总模拟时间的 70%以上。而自适应表面重建部分的计算时间包括流体表面粒子检测时间和八叉树自适应距离场创建时间，平均占据了总模拟时间的10%左右。

表4 时间性能分析(s)

通过以上性能数据的统计分析可以得到，本文的自适应流体表面重建方法可以大幅度降低拉格朗日粒子流体表面重建过程中的内存消耗和计算时间。但本文方法在创建基于八叉树的距离场时，对于每个叶子节点都要计算其 8个顶点的距离场值，存在计算冗余。后续工作中将研究八叉树叶子节点的顶点共享问题，降低计算冗余。

4 结 论

本文给出的基于八叉树的自适应流体表面重建方法可以高效地捕获拉格朗日粒子流体的表面细节特征以进行真实感绘制，该方法通过只在流体表面建立自适应距离场，使得流体表面重建过程中的内存消耗和计算复杂度，非常适合于虚拟现实应用中大规模粒子流体场景的表面提取和绘制。

[1] 柳有权. 基于物理的计算机动画及其加速技术的研究[D]. 北京: 中国科学院软件研究所, 2005.

[2] Blinn J F. A generalization of algebraic surface drawing [J]. ACM Transactions on Graphics (TOG), 1982, 1(3): 235-256.

[3] Müller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications [C]//In Proceedings of the 2003 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: CAM Press, 2003: 154-159.

[4] Yu J H, Turk G. Reconstructing surfaces of particle-based fluids using anisotropic kernels [C]//In Proceedings of the 2010 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: CAM Press, 2010: 217-225.

[5] Zhu Y N, Bridson R. Animating sand as a fluid [J]. ACM Transactions on Graphics (TOG), 2005, 24(3): 965-972.

[6] Solenthaler B, Schlafli J, Pajarola R. A unified particle model for fluid-solid interactions [J]. Computer Animation and Virtual Worlds, 2007, 18(1): 69-82.

[7] Becker M, Teschner M. Weakly compressible SPH for free surface flows [C]//In Proceedings of the 2007 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: CAM Press, 2007: 209-217.

[8] Ando R, Thurey N, Tsuruno R. Preserving fluid sheets with adaptively sampled anisotropic particles [J]. IEEE Transactions on Visualization and Computer Graphics (TVCG), 2012, 18(8): 1202-1214.

[9] He X W, Liu N, Li S, et al. Local poisson SPH for viscous incompressible fluids [J]. Computer Graphics Forum, 2012, 31(6): 1948-1958.

[10] Ihmsen M, Cornelis J, Solenthaler B, et al. Implicit incompressible SPH [J]. IEEE Transactions on Visualization and Computer Graphics (TVCG), 2014, 20(3): 426-435.

[11] Solenthaler B, Pajarola R. Predictive-corrective incompressible SPH [J]. ACM Transactions on Graphics (TOG), 2009, 28(3): 40-46.

[12] Akinci N, Ihmsen M, Akinci G, et al. Versatile rigid-fluid coupling for incompressible SPH [J]. ACM Transactions on Graphics (TOG), 2012, 31(4): 621-628.

[13] Akinci N, Cornelis J, Akinci G, et al. Coupling elastic solids with SPH fluids [J]. Computer Animation and Virtual Worlds, 2013, 24(3-4): 195-203.

[14] Marrone S, Colagrossi A, Le Touze D, et al. Fast free-surface detection and level-set function definition in SPH solvers [J]. Journal of Computational Physics, 2010, 229(10): 3652-3663.

[15] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics [J]. Journal of Computational Physics, 2003, 191(2): 448-475.

[16] Shao X Q, Zhou Z, Zhang J S, et al. Realistic and stable simulation of turbulent details behind objects in smoothed-particle hydrodynamics fluids [J]. Computer Animation and Virtual Worlds, 2015, 26(1): 79-94.

[17] 张文胜, 解 骞, 钟 瑾, 等. 基于八叉树邻域分析的光线跟踪加速算法[J]. 图学学报, 2015, 36(3): 339-344.

[18] Akinci G, Ihmsen M, Akinci N, et al. Parallel surface reconstruction particle-based fluids [J]. Computer Graphics Forum, 2012, 31(6): 1797-1809.

An Efficient Surface Reconstruction Method for Lagrangian Particle Fluid

Shao Xuqiang1, Liu Yan2, Wang Xinying1

(1. School of Control and Computer Engineering, North China Electric Power University, Baoding Hebei 071003, China; 2. International Education College, Hebei Finance University, Baoding Hebei 071003, China)

In order to simulate small-scale geometrical features of Lagrangian particle fluid surface, this paper presents an efficient octree-based adaptive fluid surface reconstruction approach. First, this approach implements the more accurate detection of fluid surface particles. Then, the approach adaptively decomposes the fluid particles’ AABB bounding boxes including fluid surface particles based on octree. Finally, the implicit distance field is only constructed in the leaf nodes including fluid surface particles. The memory consumption and computation complexity of fluid surface reconstruction depend on fluid surface not fluid volume, so our method is very suitable for large scale fluid scenes.

Lagrangian particle; surface reconstruction; fluid simulation; realistic rendering

TP 391

10.11996/JG.j.2095-302X.2016050607

A

2095-302X(2016)05-0607-07

2016-03-15；定稿日期：2016-05-03

国家自然科学基金项目(61502168)；中央高校基本科研业务费专项资金(2015MS126)；河北省自然科学基金项目(F2016502069)

邵绪强(1982–)，男，山东泰安人，讲师，博士研究生。主要研究方向为虚拟现实、计算机图形学。E-mail：shaoxuqiang@163.com