基于SPH统一框架的流体细节仿真模拟方法.pdf

2023 年第 7 期计算机与数字工程收稿日期：2023年1月13日，修回日期：2023年2月24日基金项目：国家自然科学基金项目（编号：61402205）；国家青年科学基金项目（编号：61902155）资助。作者简介：华磊，男，硕士，研究方向：计算机图形学、流体模拟。曾兰玲，女，博士，副教授，研究方向：计算机图形学。杨洋，男，博士，副教授，研究方向：计算机动画、人机交互。赵岩，男，博士，副教授，研究方向：计算机图形学。1引言流体模拟在计算机游戏、电影特效和军事训练等领域有广泛的需求和应用。经过相关学者和专家的不断努力，流体模拟已经取得了很多研究成果，并应用在不同的领域。但是现阶段研究方向都是针对单一流体的仿真模拟，因此如何通过统一框架模拟不同流体的运动状态非常具体挑战性。此外，在具体细节真实感方面，如火焰、烟雾的边缘细节优化还有待进一步提高。基于粒子的方法13和基于网格的方法4是目前流体模拟最流行的离散化不同层次的典型方案。基于网格的方法假设在无发散条件下具有不可压缩性，在拓扑转换上更有效率，但是无法精确有效的处理场的流动，容易丢失能量和几何合效果。SPH（Smoothed Particle Hydrodynamics）光滑颗粒的流体动力学5算法就是典型的拉格朗日粒子法，其基本原理是通过粒子模拟流体的运动规律，然后转换成网格进行流体渲染。如何在计算成本和真实细节之间做出的合理的选择是目前流体模拟的研究热点610，但该方面的工作基于的都是单一的流体性能方面的模拟。其中文献 1114 的研究工作是采用预先存储的数据，通过不同状态的组合生成最后的效果，文献 15 是通过加快邻域粒子的搜索速度来提高计算机模拟性能。针对多样化流体模拟研究方面的工作较少，因此，实现多样化流体模拟的统一框架具有实际应用意义。由于SPH可以有效地表示小的液体特性与低基于 SPH 统一框架的流体细节仿真模拟方法华磊曾兰玲杨洋赵岩（江苏大学计算机科学与通信工程学院镇江212000）摘要现有基于SPH方法对现实世界中存在的多种流体统一框架以及细节模拟研究较少。首先，针对现阶段缺少统一框架模拟不同流体细节仿真的问题，论文提出了在SPH框架下统一模拟流体运动的方法，通过优化后的核函数，模拟多样化的流体效果。其次，针对流体模拟中边界层细节问题，论文提出基于波动方程算法来模拟丰富的流体细节效果。实验结果表明，该统一框架和边界层细节优化模拟算法可以实现多样化的更加细节的流体效果。关键词流体；SPH；统一框架；压缩波中图分类号TM743DOI：10.3969/j.issn.1672-9722.2023.07.008Fluid Detail Simulation Method Based on SPH Unified FrameworkHUA LeiZENG LanlingYANG YangZHAO Yan（School of Computer Science and Communication Engineering，Jiangsu University，Zhenjiang212000）AbstractFirstly，to solve the problem of the lack of unified framework to simulate different fluid details，an optimized kernelfunction is presented to simulate diversified fluid effects under unified SPH framework.Secondly，for the detail problem in fluid simulation，a method is proposed to simulate rich fluid details based on the compression wave propagation model.Experimental resultsshow that the unified framework and detail simulation algorithm can achieve diversified and more detailed fluid effects as well as theless computational cost.Key Wordsfluids，SPH，unified framework，compressional waveClass NumberTM743总第 405期2023 年第 7 期计算机与数字工程Computer&Digital EngineeringVol.51No.71489第 51 卷内存成本等特点，本文研究如何利用SPH来实行多样化的流体效果，并通过SPH中的光滑核函数的优化，实现一个控制流体多样化的统一框架。细节模拟方面，由于SPH对粒子分布很敏感，在实现水波纹和烟雾波动等细节效果，可以选择在边界层粒子的位置进行细节化处理。目前的相关研究工作主要采用毛细管波和重力波模拟流体中的波纹细节效果，但毛细管波速度快、高频率难以模拟，而重力波在细节处理方粗糙，因此，文献 1617 研究在虚拟环境中模拟声波的传播。由于声波速度快，可以在流体的表面模拟声波，且不存在数值化问题等优点，本文采用声波，即压缩波对流体边界层细节进行优化处理。该工作背后的基本思想是基于SPH结合声波的方法，实现丰富多样化的3D流体细节效果。2相关工作流体模拟一直以来都是大规模图形模拟的重点研究方向，到目前为止，已经提出了许多流体模拟方法，并取得比较真实的结果。SPH是一种流体模拟算法，相对于其它的模拟算法，其特点就是简单快速，是典型的拉格朗日视角，它的基本原理就是通过粒子模拟流体的运动规律，可以有效地模拟小尺度的流体运动。现阶段在SPH上的主要研究工作是实现实时SPH流体模拟，主要采用并行加速计算方法。SPH在模拟各种流体效果的应用方面也取得了很多成果，如水体18、多相流体19和固液耦合20。Amada，et al21首次提出了基于图形处理器的SPH流体模拟，该方法在中央处理器上执行最近邻域搜索（Nearest Neighbor Search），并在图形处理器上求解N-S方程，以加速流体模拟。之后，Harada，et al22、Zhang，et al23将除渲染场景之外的所有步骤都放在GPU上，进一步提高了流体模拟的效率。现实生活中的很多流体边界层现象不能用单相流体来模拟，多流体模拟是解决这一现象的关键。Mller，et al24首先引入体积分数的概念来表示不相同的空间分布，实现混合流体的多流体模拟系统。Kim25、Misztal，et al26采用网格的方法，并假设混合流体不可压缩且不发散情况下，该方法在拓扑转换上获得较好的结果。面对混合现象由浓度差引起的扩散效应，Liu，et al27引入了一种新的多流体模拟方法，该方法能够真实地捕捉由于多种流体之间的相对运动、湍流相互作用和力变化分布而导致的复杂混合现象。湍流相互作用和力变化分布而导致的复杂混合现象。目前，主流方法是基于SPH的表面张力来优化流体的边界层模拟效果。早期工作大多集中在连续体表面力 CSF（Continuum Surface Force）方法上。CSF的主要缺点是在自由表面情况下严重依赖粒子样本。如果没有足够的粒子，会导致模拟结果数值不稳定性和不准确性的问题。此外，He等28提出在微水平上将表面张力模拟为粒子之间的分子力方法，实现流体边界模拟。最近，何小伟等29提出一种新的表面张力模型，该模型是基于扩散界面方法处理稀疏采样的自由表面。与其它表面张力相比，毛细管波尤其难以模拟，因为该方法需要一致和准确的表面张力估计。综上，本文提出基于SPH算法，并结合压缩波进行边界层细节优化的方法实现多流体细节模拟。3算法实现论文算法方案结构如图1所示。算法流程如下。计算粒子加速度开始决定流体样式通过密度光滑核函数计算粒子密度YN流体优化边界层粒子划分计算粒子物理属性更新边界层粒子粒子表面构建结束图1整体算法结构第一步：确定所要模拟的流体属性，选择本文优化的密度光滑核函数对应的属性因子m；第二步：光滑核函数代入到SPH算法，求得粒子的密度；第三步：是否选择优化流体效果，如果选择Y进入到第四步；否则直接进入到第六步；华磊等：基于SPH统一框架的流体细节仿真模拟方法14902023 年第 7 期计算机与数字工程第四步：通过流体的边界层进行边界层粒子划分；第五步：通过声波公式求出代更新的边界层密度；第六步：进行粒子的初始化 求粒子的物理属性（例如重力、粘滞力等），得出粒子的加速度；第七步：通过MC算法进行流体表面构建得到想要的流体模拟效果。3.1SPH等式SPH是最流行的拉格朗日方法。它能有效地模拟小尺度自由表面流体。不可压缩牛顿流体的运动由 N-S方程描述，而 SPH 方法中选取简易的N-S方程，分析得出作用在一个微粒上的作用力由三部分组成，如式（1）所示：F=Fexternal+Fpressure+Fviscosity（1）其中Fexternal称为外部力，一般就是重力，如式（2）所示：Fexternal=g（2）Fpressure是由流体内部的压力差产生的作用力，如式（3）所示：Fpressure=-p（3）Fviscosity是由粒子间的速度差引起的，跟流体的黏度系数以及速度差有关，如式（4）所示：Fviscosity=2u（4）根据式（2）式（4）带入到式（1），并带入公式F=a，可以得到：a=g-p+2u（5）加速度则为a=g-p+2u（6）以上就是粒子运动学的计算方法，其中u表示速度，表示密度，p表示压力，g为重力表示外力，表示粘度系数。SPH算法与其他的流体力学中的数学方法类似，同样涉及到“光滑核”的概念，因为流体中每个位置参与运算的值都是由周围一组粒子累加起来的，则该处某项属性Attribute的累加公式，如式（7）所示：Attribute(ri)=jAttributejmjjW(ri-rjh)（7）其中Attributej是要累加的某种属性，mj和j是周围粒子 j 的质量和密度，ri和 rj是粒子 i 和 j 的位置，h是光滑核半径，而W就是光滑函数。求解 N-S 方程需要满足动量守恒和质量守恒。由于SPH粒子的数量在产生后不会发生变化，所以SPH中的粒子i在其运动过程中可以满足质量守恒。SPH粒子i的加速度可以由式（5）计算得出：ai=duidt=Fii（8）其中Fi是等式（1）外力之和。ai为本文更新粒子i的位置的加速度，推算得到的加速度公式如式（6）所示。3.2基于SPH的优化由于基于SPH的原始方法以及扩展方法都无法实现一个统一的流体控制框架，针对这个问题，本文对SPH的光滑核函数进行优化，计算粒子密度的光滑核函数可以改写成式（9）所示：W=1|ri-rj-1 2mkh（9）其中h为核半径，ri-rj为半径内的j位置的粒子与圆心的距离，m为决定流体的属性因子，k为粒子间距离倍数。接下来根据式（7），用代替 Attribute可以得到式（10）：(ri)=jjmjjW(ri-rjh)=jmjW(ri-rjh)（10）计算密度中将不再使用之前光滑核函数称为Poly6的光滑核函数，而是使用式（9）。现将其带入到式（10）中，则 ri处的密度计算公式最终为式（11）：(ri)=jmj1|ri-rj-12mkh（11）根据式（7），压力计算中使用的光滑核函数称为Spiky函数，因此在ri处由压力产生的加速度部分为式（12）所示：aipressure=m45h6j(pi+pj2ij(h-r)2 ri-rjr)（12）其中r=|ri-rj。根据式（7），粘滞力计算中使用的光滑核函数称为Wviscosity，在ri处由粘滞力产生的加速度部分为式（13）所示：aiviscosity=m45h6jui-ujij(h-|ri-rj)（13）将式（12）、式（13）代入到公式中得到式（14）所示：1491第 51 卷a(ri)=g+m45h6j(pi+pj2ij(h-r)2)+m45h6jui-ujij(h-r)（14）通过修改本文采用的密度光滑核函数参数，对基本的流体仿真模型进行复杂形变，包括流体（火、烟）抖动形变和流体溅落形变。对基本流体仿真模型进行的复杂形变通过修改密度光滑核函数中的属性因子m、分母系数k实现，在可控区域内，m=1、0.5、h（核半径）是最接近现实流体的仿真构建，粒子显示结果如图2所示。（a）m=1粒子图（b）m=0.5粒子图（c）m=h粒子图图2统一框架呈现的粒子效果雏形3.3基于声波的细节处理本文在流体仿真细节优化处理中，使用声波的波动方程来更新密度变化，如式（15）所示：D2outerDt2=c22outer（15）其中c为波速，outer表示流体表面边界层中的粒子关于位置x和时间t的波强度的一个测量。c=k=k 0（16）式（16）所示，其中是表面张力，k是波数，0的流体密度。对于水，论文使用0=1.0*103kg/m3，对于气体相对应的可以直接使用空气的密度0=1.293kg/m3。外部粒子更新的粒子边界层主流体内部粒子图3边界层粒子划分尽管声波和水体的毛细管波虽然很多的行为相似，但是还有一定的区别：液体毛细管波只在表面传播，声波可以在液体的体积内传播。如果只是简单地求解式（16），模拟流体会得到不一样的模拟效果，尤其在模拟流体的粒子数目过多的情况下。为了实现更加逼真效果，针对上述问题，论文采用的解决方案是检测流体表面边界层的粒子，并对粒子进行范围的划分，让声波在不同边界层范围内进行传播，该方法可以实现较好的流体细节效果。因此，本文引入了流体的边界层划分算法，如式（17）所示：=5.0vLv（17）其中式（17）表示布拉修斯对平板层流边界层求得的精确解，其中v表示运动粘度，v表示的是主体流速，表示的是边界层厚度，边界层厚度指的就是壁面到流速与主体流速v相差不到1%的流体层厚度，如式（18）所示：v=（18）其中表示流体的动力粘度（粘性的度量，也称粘性系数），表示液体密度。根据式（14），可以求得每个粒子的加速度情况，因为速度是矢量，那么主流体的速度，如式（19）所示：v=ia(ri)*tn（19）其中公式t表示时间，n表示粒子数。如图3所示，接下来就是将式（19）中得到的主流体流速代入到式（17）中，得到了流体的边界层，将式（16）带入到式（15）可以得到outer，更新边界层中粒子的密度为outer。就可以将声波形成的细节效果带入到流体的表面。以水流为例，实现效果华磊等：基于SPH统一框架的流体细节仿真模拟方法14922023 年第 7 期计算机与数字工程与没有声波丰富细节明显的区别如图4所示，对比效果，图（a）细节化程度明显弱于图（b）采用的本文方法。且图（a）中文献 30 的实现方法需要构造一个密集的点集的拓扑变化，添加了算法的复杂度，因此本文中的方法更加容易实现。（a）文献方法30（b）本文方法图4（a）（b）表面波纹对比yyxx（a）文献方法30（b）本文方法图5粒子方法对比（a）原始方法（b）细节化粒子效果（c）细节化绘制效果图6丰富的细节化火焰为了进一步验证本文细节优化处理的有效性，通过粒子实验来进行验证，计算参数有：流体初始密度为=1000kg/m3，流体模拟粒子数为2200个，边界粒子数为 120个，粒子影响半径 r=0.02m。从图5中可以看出，在其他条件相同的情况下，图（a）相对图（b）的流体表面均匀，图（b）的凹凸效果更加明显，因此本文方法可以提高流体表面细节仿真效果。4实验结果本文实验是在一台配备英特尔酷睿i7-8700K3.70GHz、英伟达GTX1080显卡、16G内存，操作系统为Windows 10的电脑上测试文中的算法。整个实现基于 OpenGL，使用的是 MC 算法（MarchingCube Algorithm）进行流体液面绘制。这里给出了一些实验结果，这些结果显示了本文所提出的方法的有效性。经过实验不仅完成了整个基于物理的流体方程的推导，同时也将计算的结果实时绘制出来。实验结果验证了本文算法所对产生丰富的细节多样化流体效果，论文中的实验例子包括多样化的流体雏形效果（如图2），图2的实验得到了一个基于SPH不局限单一流体的框架。在图2多样化的流体雏形的框架基础上，添加粒子数，为粒子添加对应流体的颜色，结合MC的流体表面绘制，呈现小型流体波纹效果（如图3），火焰燃烧（图6）的流体效果，以及烟雾吹动的效果（图7）。这些例子说明了本文算法的正确性，以及类似现实世界中的流体的波动效果。在图 5 中，文中方法（图（b）、（c）比较Mllers（图（a）的方法在火焰上实现的对比，可以看出文中的方法优化流体细节方向这一块，使得流体显示的效果更加鳞次栉比。经过图7的实验，图（a）是将何晓伟等30实现毛细管波的方法来丰富烟雾的细节化，效果是近似现实效果，但计算时间长，图（b）、（c）实验展示本文的方法应用到烟雾模拟中的，图（b）中粒子细节效果更加丰富，再进行表面绘制得到的图（c）的结果，经过实现得到近似现实世界烟雾效果，在烟雾刚开始的起始阶段基于文中算法得到的细节效果好于文献 30 的方法。由表12中的数据结合图5和图7的模拟效果可以得出在相同的条件下，随着时间的增加，本文的最高速度与最低速度差明显大于Mllers的方法和文献 30 的方法，进一步得出本文方法速度波动大，流体模拟细节越显著。本文采用的算法实现多样化的流体效果，是局1493第 51 卷限小型流体，大型流体的效果会出现过度偏差细节化优化也会出现失真的效果，并且流体运动过长，细节方面也会出现精度问题。算法在细节优化步骤上是基于压缩波的，相对于液体的毛细管波等现实波纹细节方面，并不意味着物理上的精确。文中算法不指定流体的具体密度，是根据单一的密度来实现多样化的流体细节优化的效果。最后，基于文中算法的基础得到的结果在某种程度上是取决于表面的粒子颗粒上的。（a）原始方法（b）细节化粒子效果（c）细节化绘制效果图7丰富的细节化烟雾表1N=45000数据时间步长/s01*10-2510Mller方法最高速度与最低速度差00.019.4620.52He方法最高速度与最低速度差00.0314.7429.12本文方法最高速度与最低速度差00.0415.7336.24表2N=70000数据时间步长/s01*10-2510Mller方法最高速度与最低速度差00.1610.3723.24He方法最高速度与最低速度差00.3716.8632.52本文方法最高速度与最低速度差00.5519.7337.945结语论文提出了一种框架，基于SPH实现丰富细节的多样化的流体模拟。基于SPH的基础下优化光滑核函数实现了多样化的流体框架不再局限于单一的流体模拟，并且结合声波实现细节优化的方法。经过实验，表明了这是一个合理的假设。经过实验研究进一步表明了，在SPH框架声波更新的表面粒子可以为优化流体的细节提供很好的方向。本文中的方法可以实现多样化的流体模拟，但是将该方法应用于大型流体模拟是我进一步研究的工作。另外，在未来，希望基于以上算法的基础，进行理论的研究，实现多种多样流体控制。并且想进一步研究气液耦合方面的压强问题，以获得更精确的细节优化效果。更进一步的话，将通过机器学习来提高文中算法的性能，提高整体算法的模拟速度。参 考 文 献1Solenthaler B，Pajarola R.Predictive-corrective incompressible SPH J.ACM Trans.Gr.，2009，28（3）：1-6.2Charypar D，Gross M.Particle-based fluid simulation forinter-active applications J.In：AM Siggraph/Eurographics Symposium on Computer Ani-mation.Eurographics association，2003：154-159.3 Liu M.B.，Liu G.R.Smoothed particle hydrodynamics（SPH）：an overview and recent developments J.Arch.Comput.Methods Eng，2010，17（1）：25-76.4Kuznik F，Obrecht C，et al.LBM based flow simulation using GPUcomputing processor J.Compute，2010：59（7）：2380-2392.5MONAGHAN J J.Simulating free surface flows with SPHJ.Comput.Phys.，1994，110：399-406.6Yang M，Li X，Liu Y，et al.A novel surface tension formulation for SPH fluid simulation J.The Visual Computer，2016，33（5）：597-606.7Bender J，Koschier D.Divergence-Free SPH for Incompressible and Viscous Fluids J.IEEE Transactions on Visualization and Computer Graphics，2016，23（3）：1.8Yan Z，Zhu X，Liu Y，et al.An improved method of position based fluids J.Coll.Math.，2016，32（1）：38-43.9Boyd L，Bridson R.MultiFLIP for energetic two-phase fluid simulation J.ACM Trans.Gr.，2012，31（2）：1-12.10Huang T，Li D，Li L.Adaptive time-stepping particlefluid motion simulationC/In：International Conference on Computer Application and System Modeling.（下转第1616页）华磊等：基于SPH统一框架的流体细节仿真模拟方法1494第 51 卷International Symposium on Information，Computer，andCommunications Security Information，2009：218-227.13W Wu，Y Xiao，W Wang，et al.k-symmetry model foridentityanonymizationinsocialnetworks C/NewYork：EDBT 10：Proceedings of the 13th InternationalConference on Extending Database Technology，2010：111-122.14Ying X，Pan K，Wu X，et al.Comparisons of randomization and K-degree anonymization schemes for privacypreserving social network publishingC/New York：SNA-KDD 09：Proceedings of the 3rd Workshop on Social Network Mining and Analysis，2009：1-10.15Zhou B，Pei J.Preserving Privacy in Social NetworksAgainst Neighborhood Attacks C/Mexico：2008 IEEE24th International Conference on Data Engineering，Cancun，2008：506-515.16Boldi P，Bonchi F，Gionis A，et al.Injecting Uncertainty in Graphs for Identity Obfuscation J.Proceedings ofthe Vldb Endowment，2012，5（11）：1376-1387.IEEE，2010：117-121.11Jeong S H，Solenthaler B，Pollefeys M，et al.Data-driven fluid simulations using regression forestsJ.ACMTrans.Gr.，2015，34（6）：199.12Yang C，Yang X，Xiao X.Data-driven projection method in fluid simulation J.Comput.Anim.Vir.Worlds，2016：415-424.13Okabe M，Anjyor K，Onai R.Creating fluid animationfrom a single image using video databaseJ.Comput.Gr.，2011，30（7）：1973-1982.14Sato S，Morita T，Dobashi Y，et al.A data-driven approach for synthesizing high-resolution animation of fireJ.In：Proceeding of the Digital Production Symposium，2012：37-42.15刘恒星，范湘涛，刘健，等.WebGL技术下的SPH流体模拟方法 J.小型微型计算机系统，2017，38（10）：2406-2411.LIU Hengxing，FAN Xiangtao，LIU Jian，et al.SmoothedParticle Hydrodynamics Fluid Simulation Method Basedon WebGLJ.Journal of Chinese Computer Systems，2017，38（10）：2406-2411.16Jeschke S，Wojtan C.Water Wave Animation via Wavefront Parameter InterpolationJ.ACM Transactions onGraphics，2015，34（3）：1-14.17Jeschke S，Skrivan T，Muller-Fischer M，et al.Watersurface waveletsJ.ACM Transactions on Graphics（TOG），2018，37（4CD）：1-13.18Solenthaler B，Pajarola R.Density contrast SPH interfaces J.InProc.of SCA，2008：211-218.19Yu J，Wojtan C，Turk G，et al.Explicit mesh surfacesfor particle based fluids J.Comp.Graph.F orum（Eurographics），2012，31（4）：815-824.20Amada T，Imura M，Yasumuro Y，et al.Particle-basedfluid simulation on GPU J.In：Proc.of ACM Workshopon General-purpose Computing on Graphics Processors，2004：342-343.21Harada T，Koshizuka S，Kawaguchi Y.Smoothed particle hydrodynamics on GPUs C/In：Proc.of ComputerGraphics International，2007：63-70.22Zhang Y，Solenthaler B，Pajarola R.Adaptive samplingand rendering of fluids on the GPU J.In：Proc.of theFifthEurographics/IEEEVGTCConferenceonPoint-based Graphics，2008：137-146.23Muller M，Solenthaler B，et al.Particle-based fluid-fluid interaction J.ACM SIGGRAPH/Eurographics Symposium on Computer Animation，2005：237-244.24Kim B.Multi-phase fluid simulations using regional level sets J.ACM Trans Gr，2010，29（6）：175.25Misztal MK，Erleben K，Bargteil A，et al.Multiphaseflow of immiscible fluids on unstructured moving meshesJ.IEEE Trans Visual Comput Gr，2014，20（1）：4-16.26Liu S，Liu Q，Peng Q.Realistic simulation of mixing fluids J.Visual Comput，2011，27（3）：241-248.27Akinci N，Akinci G，Teschner M.Versatile surface tension and adhesion for SPH fluidsJ.ACM Trans.Graph.，2013，32（6）：1-8.28He X，Wang H，et al.Ro-bust simulation of sparselysampled thin features in SPH-based free sur-face flowsJ.ACM Trans.Graph.，2014，34（7）：1-7.29Mercier O，Beauchemin C，Thuerey N，et al.Surface turbulence for particle-based liquidsimulationsJ.ACMTrans.Graph.（SIGGRAPH Asia），2015，34：6.30Sheng YANG，Xiaowei HE，Huamin WANG，et al.Enriching SPH simulation by approximate capillary wavesP.Computer Animation，2016.（上接第1494页）彭擎宇等：基于三元闭包的不确定社交网络隐私保护算法1616