自适应分层定位流体仿真模型

梁欣鑫, 朱晓临, 刘洋洋, 闫志伟

(合肥工业大学 数学学院,安徽 合肥 230009)

自适应分层定位流体仿真模型

梁欣鑫, 朱晓临, 刘洋洋, 闫志伟

(合肥工业大学 数学学院,安徽 合肥 230009)

文章将自适应分层模型与3种现有的粒子方法相结合,得到3种改进的粒子方法,并对流体粒子进行大小粒子分层,可减少模拟使用的粒子,从而降低了原始方法的计算量;而通过使用较小粒子覆盖于流体表面,使得改进后的方法表面细节效果更佳。通过嵌入自适应分层模型,使得改进后的3种方法的计算效率得到了提高,表面细节也得到了保留甚至更佳。

流体模拟;自适应;分层

流体模拟是计算机图形学中重要的研究内容,特别是对液体的模拟。光滑粒子动力学(smoothed particle hydrodynamics,SPH)方法是一种用于液体模拟的主流粒子方法,文献[1-4]成功地将天体物理学中的SPH 方法应用到模拟连续固体力学和流体力学中,实现了流场中的激波强间断现象;但是由于SPH方法中目标粒子的密度对周围粒子数量和位置分布敏感,加上模型的非结构性导致了模拟流体的不可压缩性难以保证,表面细节效果不理想。文献[5]通过对三次B样条核函数的测试,指出了SPH方法不能保证张力的稳定性,可以通过放大有效应力来减小扰动。文献[6-7]中根据粒子的稠密程度通过SPH方法动态调节粒子的支集半径,虽然加快了SPH方法的模拟速度,但是降低了流体的表面细节效果。文献[8]在SPH方法中加入一个人工压力项，很好地解决了流体模拟中表面张力不稳定的问题,提升了模拟效果。文献[9]构建了一种可以对流体的自由表面交互模拟的SPH模型,此后SPH方法才正式广泛应用于模拟不可压缩流体以及多相流体模拟。文献[10]提出了一种在GPU通用平台上的SPH模拟方法,计算效率约为传统SPH方法的6倍,解决了SPH方法实时模拟的难题。文献[11]提出了一种基于SPH的非均匀粒子流体模拟方法,引入由粒子的位置、所受压力和黏性力3种因素决定的精细度参数,减少了同等场景下模拟需要的粒子数,提高了计算速度。

文献[12]提出了半隐式移动粒子(moving-particle semi-implicit,MPS)方法。该方法提出之初是为了解决不可压缩流体的模拟问题,用于模拟水体溃坝。MPS方法是一种半隐式方法,对流项是通过求解Poisson方程得到,而其他物理量则是通过显式方法得到。相比于SPH方法,MPS提出的时间还比较短,但因其能保证流体的不可压缩性,越来越多的学者关注MPS方法,将其应用到新场景或者对方法本身提出改进。

文献[13]使用MPS方法对流体的晃动问题进行了研究；文献[14]在模拟多相流体以及中等规模的流体时,结合网格对MPS进行改进；文献[15]对MPS方法中使用不同的核函数进行了讨论,从核函数的选取上增强了MPS方法的稳定性；文献[16]使用泰勒展开,对离散后粒子的局部坐标进行近似,确保了离散精度的一致性。

文献[17]提出了一种基于定位动态学(position based dynamics,PBD)的流体模拟方法,通过物体运动的限制方程来修正物体的位置,运算效率高,模拟流体运动具有很真实的动态感。文献[18]通过在薄层表面适当修改粒子间距来动态模拟的表层粒子重采样方法来模拟流体回落的过程,试验结果的实效性表明局部流体在适当修正位置后可以更好地体现出模拟流体的真实感。文献[19]结合SPH方法和PBD方法又提出了定位流体模拟(position based fluids,PBF)方法,通过粒子密度限制方程来校正位置,模型中人工添加的压力项很好地避免了粒子小范围飞溅或聚集的问题,模拟效果很好,适用于大规模的场景模拟。

文献[20]在模拟横向流体运动中引入了一种基于粒子几何特征尺寸的自适应采样算法。该算法允许在几何复杂的区域增加计算量来细化流体局部特征;通过粒子融合方式减少流体内部的粒子数来加速模拟;根据模拟视角的重要程度调节采样粒子密度来进一步提升细节模拟效果。但这些优点在粒子较少时效果并不明显。文献[21]在模拟水滴落入水杯的场景中,将水杯中模拟水的粒子分为大小不同的若干层,以提高模拟效率。该分层模拟的背景比较简单,相邻层之间的粒子仅需粒子融合即可保持模型的稳定。文献[22]提出了一种统一的粒子框架来模拟固体、液体之间相互作用的运动情形,很好地解决了传统粒子方法固液交互模拟中固体边界难以处理的问题。

针对模拟不可压缩流体往左右两侧溃坝的场景,本文将自适应分层模型与现有的3种粒子方法相结合,得到3种改进后的粒子方法。针对溃坝的液体采用分层技术:液体表面的粒子为小粒子,保证了液体表面细节的模拟效果;下层的粒子为大粒子,这些大粒子的使用保证了流体内部区域的整体结构,可以减少总体粒子数,提高模型运行速度。本文的改进方法对原始方法的计算流程并未进行改动,而是在保持原始方法计算流程的基础上,在方法的初始处加入分层,而后下次迭代计算进行之前,对流体粒子进行自适应分层。通过这种改进,本文有效地将3种经典的粒子方法进行改进,在保持表面细节的前提下,降低了原始方法的计算复杂度,改进后的方法具有更好的表面细节。本文的模拟背景比文献[21]复杂,而与文献[20]相比,本文的改进方法效率更高。

1 3种粒子方法简介

1.1 SPH方法

SPH方法的基本思想是将粒子的相关物理量用积分表述,然后通过核函数对其进行离散,进而对描述流体的微分方程进行数值离散。

SPH方法的基本公式如下:

(1)

其中,Ai为需要求解的物理量,例如密度、压强或黏性力等;h为支集半径;粒子i的支集表示以粒子i为中心,h为半径的圆形范围;X={x1,x2,…,xn}中xj(j=1,2,…,n)表示粒子i支集半径内粒子的位置;mj为粒子j的质量;W1(xi-xj,h)为一个核函数。在计算过程中用到的核函数W1、W2、W3参考文献[7]。通过(1)式可以得到影响粒子运动的压强力和黏性力,从而计算粒子的加速度,进而计算出粒子的位置变化。

SPH方法的流程如图1所示。

图1 SPH方法的流程

由图1可以看出,SPH方法整个过程都是显式求解的,求解过程并不能完全保证质量和动量的守恒,因此其存在精度不高的缺点,特别是自由边界与介质分界处的粒子,此区域的粒子能用来求解中心粒子的数目比流体内部粒子更少,相比而言精度也就更差。此外,SPH方法的表面张力稳定性较差,该方法本身对核函数的使用有较高的要求,比如要求紧支撑性,要求函数为高斯函数等。

1.2 MPS方法

MPS的基本思路与SPH方法类似,都是通过使用核函数来对粒子之间的相互作用进行建模。不同于SPH方法,MPS方法不离散对流项,其计算步骤分为显式和隐式2个部分。显式部分通过外力和边界条件对粒子的速度和位置进行直接修正;而隐式部分则是通过显式修正之后造成的粒子数密度变化,得到一个包含压力项的Poisson方程,求解方程得到压力项,再对粒子的速度和位置进行第2次修正,这样粒子的密度刚好不变,从而满足不可压缩条件。

在MPS方法中,为了描述粒子的疏密程度,定义粒子i处的粒子数密度为:

(2)

其中,h为支集半径;粒子j表示以粒子i为中心,h为半径的圆形范围内的粒子,xj(j=1,2,…,n)表示粒子i支集半径内粒子的位置;W(xi-xj,h)为一个核函数。粒子数密度实质上是密度的一种离散表达,不可压缩性要求流体的密度不变,相当于粒子数密度为常量,记为n0。

对描述流体运动的微分方程进行分解,先对粒子施加黏性力和外力,得到粒子的临时粒子数密度〈n*〉i以及临时速度与位置,然后得到压强的Poisson方程如下:

(3)

其中,Pn+1为下一时刻的压强。求解(3)式,得到新的压强后,对临时速度与位置进行修正。

MPS方法的流程如图2所示。

图2 MPS方法的流程

由图2可以看出,MPS方法只有在求解压强时才是隐式求解的，虽然隐式的方法保证了守恒,这样得出的解精度较好,但是求解Poisson方程变成整个算法中最为耗时的过程。

1.3PBF方法

PBF算法的的主要思想是:首先将流体模拟中的不可压缩性转化为关于位置的限制方程;然后通过对限制方程进行求解,得到新的位置;最后根据新的位置算出粒子的速度。该方法通过引入PBD的思想,先对粒子的位置应用外力进行第1次修正,然后找出粒子的相邻粒子,类似于SPH方法,显式地算出第2次修正的参数;然后进行第2次修正;最后使用新的位置来计算出粒子的速度。位置修正的方法如下:

(4)

其中,h为支集半径;粒子i的支集表示以粒子i为中心,h为半径的圆形范围;X={x1,x2,…,xn}中xj(j=1,2,…,n)为粒子i支集半径内粒子的位置;mj为粒子j的质量;W1(xi-xj,h)为一个核函数;λ为位置修正参数;scorr为加入的一个非负的人工压力项,用以避免粒子聚集成簇,其数值参数同文献[19]。在计算过程中用到的核函数W1、W2、W3参考文献[7]。

MPS方法的流程如图3所示。

图3 PBF方法的流程

PBF方法是一种高效的方法,与SPH方法类似,都是显式求解的；但是PBF方法的计算结果特别是表面张力也会对使用的核函数存在依赖性。此外,针对加入的人工压力项scorr,其设置一般是通过多次实验,选取适当值。

上述3种方法均有不同的优缺点,但在模拟大规模场景时,这3种粒子方法都需要大量的粒子,计算量很大。因此,本文引入自适应分层技术对上述3种方法进行改进。

2 本文模型

本文粒子模型的改进体现在粒子分层模型及其相适应的粒子融合分离机制2个方面,粒子的属性计算精度按照粒子速度大小自适应变化。

2.1 粒子融合分离机制

文献[20]在模拟流体和固体的碰撞场景时,按照模拟的区域判定粒子是否需要融合或分离,标记需要处理的粒子;然后根据粒子周围其他粒子的位置分布判定是否满足粒子融合或分离的条件。若满足,则依照标记处理;否则就更新粒子属性,直至下一时刻。这种简单的双重判定的方法使得“新粒子”的位置分布更合理，但是这种粒子融合分离机制在分离时要求粒子分布对称,在融合时要求适当的空间分配在周围粒子集中,这些比较难满足,循环多步后可能会造成标记粒子数量过多而影响粒子模型的稳定性。

本文的粒子模型首先在初始时刻T=0时先在一个40×20×20的水槽左半边平铺4层大粒子,每个大粒子半径为0.5,质量为1;然后在大粒子的上方和右边再分别放置6层和3列小粒子,小粒子的半径为rl=0.25(rl是实验模拟中最小的长度单位,也是文中的支集半径和搜索范围的距离衡量单位),质量为0.125。表层的小粒子可以提高模拟流体的表面细节,内部的大粒子支撑模拟流体的运动形态,同时减少了模拟中需要的粒子数,减小计算量。文中的粒子融合分离机制是通过计算目标粒子支集半径内异类粒子和同类粒子之间的数量比值γ来判定目标粒子是否仍处在目标粒子区域,其计算公式如下:

(5)

其中,N为目标粒子支集半径内不同粒子的数量;m为目标粒子支集半径内不同粒子的质量;X为目标粒子的种类；B为小粒子；A为大粒子。若大粒子支集半径内小粒子和大粒子的数量比超过了给定值γ,则判定该大粒子进入了小粒子区域,需要分离成2个小粒子。为了估算实验模拟中γ的数值,本文虚拟一个球形结构,球体内的上半部分均匀放置小粒子,下半部分放置大粒子。当球体的半径为大粒子的支集半径时,实验计算得到,大粒子和小粒子的比值为17/37,并以此作为大粒子分离时γ的临界参考值;当球体的半径为小粒子的支集半径时,这个比值为7/17,以此作为小粒子融合时γ的临界参考值。在估算出粒子分离在和融合的临界参考值后,本文还在实验中分别对上述的临界参考值进行了微调以检验粒子模型的稳定性,测得的结果显示大粒子分离时的稳定区间为[0.4,0.56],小粒子融合时的稳定区间为[0.38,0.48]。当所选用γ的参考值不在稳定区间内时,粒子模型在程序运行时粒子的结构层次混合,影响稳定性,降低模拟效果。这种新的粒子融合分离机制适用绝大多数的分层粒子模型的流体模拟场景。

2.1.1 粒子融合

粒子融合判定条件如下:

(1) 小粒子进入大粒子的区域,文中γ取值为0.42。

(2) 小粒子的支集半径内有其他同类粒子;若有多个,则取距离最近的一个。

融合的大粒子属性如下:

(1) 密度取粒子静止密度ρ=ρ0。

(2) 质量取2个较小质量粒子之和M=2m。

(3) 速度按照动量守恒定律

2mv=mv1+mv2⟹v=(v1+v2)/2。

(4) 新位置取2个质量较小粒子融合前位置的中点,即x=(x1+x2)/2。

简单的粒子融合二维图如图4所示。图4a判定小粒子满足融合条件,选取范围内最近粒子进行融合;图4b融合后的大粒子位置为融合前2个小粒子的中点。

图4 粒子融合二维图

2.1.2 粒子分离

粒子分离判定条件如下:

(1) 大粒子进入小粒子区域,文中γ取0.5。

(2) 该大粒子影响范围内没有其他同类粒子。

分离的2个小粒子属性如下:

(1) 密度取粒子静止密度ρ=ρ0。

(2) 质量取较大质量粒子的1/2,m=M/2。

(3) 速度按照运动学规律满足动量守恒和能量守恒

v1=v2=v。

(4) 2个新粒子的位置中点为分离前原粒子的位置,位置连线为原粒子速度方向的延长线,2个粒子距离为2.8rl。

简单的粒子分离二维图如图5所示。图5a判定大粒子满足分离条件;图5b分离成的2个小粒子位置在原大粒子的速度方向上,两者距离的中点即是分离大粒子的原位置。

图5 粒子分离二维图

2.1.3 孤立粒子

若一个小粒子的支集半径之内,全是比其大的粒子,则该粒子并不能进行融合,因此这种粒子就形成了孤立粒子;同理,对于大粒子亦如此。为了简便,以下本文仅对孤立的小粒子进行讨论。

对于孤立的小粒子,其作用域内实际上是较大的粒子存在的流体中部和下部区域,并且大粒子往往是动能较小的粒子,此时可通过对孤立粒子加一个推力,该推力的方向为粒子的运动速度方向,即

其中,λ的取值范围为[0.1,0.3],本文取λ=0.1。反之,对于小粒子区域中的孤立大粒子,同样给其加一个与速度方向相反的推力。

2.2 计算精度自适应变化

因为流体运动的局部连续性,模型中的粒子在速度较低时,其周围粒子的速度变化很小,其对周围粒子的作用力很微弱。因此,当粒子运动缓慢,通过降低消极粒子的支集半径或者减少计算的迭代次数并不会对下一时刻的数值计算产生较大影响。但是,当粒子运动速度加快时,活跃粒子的运动对周围粒子的位置开始敏感,因此扩大它的支集半径或者增加计算的迭代次数可以得到更精确的结果。本文根据粒子在时间间隔Δt内的位移大小将粒子分为活跃粒子、普通粒子和消极粒子3种粒子。3种粒子占总粒子数的百分比分别为10%、30%和60%,并通过输出实验数据确定时间间隔Δt内粒子位移的分界点。具体数据及每种粒子的支集半径迭代次数见表1所列。

表1 活跃粒子、普通粒子和消极粒子的分类及其支集半径

3 算法流程

本文分别将上述自适应分层方法与SPH、MPS、PBF方法相结合,建立自适应分层SPH方法、自适应分层MPS方法和自适应分层PBF方法。改进后的方法流程如图6～图8所示。

图6 自适应分层SPH方法流程

图7 自适应分层MPS方法流程

图8 自适应分层PBF方法流程

4 实验结果

本文的实验平台为Windows 7系统,程序开发的软件选用Microsoft Visual Studio 2008,三维图形的编程接口采用OpenGL。算法实现的硬件环境是Intel(R) Core(TM) i3-4130的CPU、4 G内存的PC机。原始3种方法和改进后方法的粒子模拟结果比较如图9所示。

图9 流体的粒子运动

在模拟中,原始方法和加入自适应分层改进后方法的帧率(帧率指1 s内的帧数)比较见表2所列。

表2 原始方法和本文方法的帧率比较

原始方法和改进后方法的模拟结果经表面重建及渲染后的比较如图10所示。

图10 渲染后流体运动

由表2可以看出,加入自适应分层改进之后的方法,其帧率的提高率至少提升了20%,其中MPS的提升最为明显,提高率达到了65.24%,而另外2种方法都不超过50%。此外,由图9及图10可得:分层后的粒子模型可以在模拟同等场景的流体时,减少约25%的粒子数;计算效率较原始方法快24%～66%;模拟流体表面细节效果可以近似达到全部用一种粒子进行计算得到的效果。

5 结 论

针对横向流动液体,本文提出使用嵌入自适应分层的改进的3种粒子方法,将粒子模型下层的大量粒子替换为体积更大的粒子,有效地减少了模拟中的粒子数;再通过一种较为简单的粒子融合和分离机制，很好地保证了分层模型的稳定性。文中改进后的方法相比于原始方法,可以自适应地调节粒子的计算精度,减少了消极粒子的计算量,提高了模拟运算速度；通过对表层高速粒子更精确的计算,能够更好地模拟水流的表面细节。改进后的方法特别在计算效率上有较明显的提高。

然而,由于本文模拟流体时的粒子数量偏少,分层模型在减少模拟粒子数量及对水流表面细节更精细地刻画方面的优势并没有完全体现出来。通过显卡的并行计算的GPU加速,文中模型的粒子数可以达到百万级时,本文模型的优势将会愈加突出,渲染后的模拟效果将更加逼真。

[1] GINGOLD R A,MONAGHAN J J.Smoothed particle hydrodynamics:theorya and application to nonspherical stars[J].Monthly Notices of the Royal Astronomical Society,1977,181:375-389.

[2] MONAGHAN J J.An introduction to SPH[J].Computer Physics Communications,1988,48(1):89-96.

[3] MONAGHAN J J.Smoothed particle hydrodynamics[J].Annual Review of Astronomy and Astrophysics,1992,30:543-574.

[4] MONAGHAN J J.Simulating free surface flows with SPH [J].Journal of Computational Physics,1994,110(2):399-406.

[5] SWEGLE J W,HICKS D L,ATTAWAY S W.Smoothed particle hydrodynamics stability analysis[J].Journal of Computational Physics,1995,116(1):123-134.

[6] JOHNSON G R,STRYK R A,BEISSEL S R.SPH for high velocity impact computations[J].Computer Methods in Applied Mechanics & Engineering,1996,139(1/2/3/4):347-373.

[7] OWEN J M,VILLUMSEN J V,SHAPIRO P R.Adaptive smoothed particle hydrodynamics[J].Astrophysical Journal Supplement,1995,116(2):155-209.

[8] MONAGHAN J J.SPH without a tensile instability [J].Journal of Computational Physics,2000,159(2):290-311.

[10] 温婵娟,欧嘉蔚,贾金原.GPU通用计算平台上的SPH流体模拟[J].计算机辅助设计与图形学学报,2010,22(3):406-411.

[11] 谭诗瀚,段茗,杨红雨.非均匀粒子流体模拟[J].计算机工程与设计,2011,32(8):2760-2763.

[12] KOSHIZUKA S,OKA Y.Moving-particle semi-implicit method for fragmentation of incompressible fluid[J].Nuclear Science and Engineering,1996,123(3):421-434.

[13] YOON H Y,KOSHIZUKA S,OKA Y.A particle-gridless hybrid method for incompressible flows[J].International Journal for Numerical Methods in Fluids,1999,30(4):407-424.

[14] PREMOE S,TASDIZEN T,BIGLER J.Particle-based simulation of fluids[J].Computer Graphics Forum,2003,22(3):401-410.

[15] ATAIE-ASHTIANI B,FARHADI L.A stable moving-particle semi-implicit method for free surface flows[J].Fluid Dynamics Research,2006,38(4):241-256.

[16] 张帅,楼一民,邢菲,等.一种改进的MPS粒子作用模型[J].计算力学学报,2013,30(1):124-129.

[20] ADAMS B,PAULY M,KEISER R,et al.Adaptively sampled particle fluids[J].ACM Transactions on Graphics (TOG),2007,26(3):1-8.

[21] HONG W,HOUSE D H,KEYSER J.Adaptive particles for incompressible fluid simulation[J].The Visual Computer，2008,24(7):535-543.

(责任编辑 闫杏丽)

Adaptive layered position based fluids simulation model

LIANG Xinxin, ZHU Xiaolin, LIU Yangyang, YAN Zhiwei

(School of Mathematics, Hefei University of Technology, Hefei 230009, China)

This paper modifies three methods which are based on particles with the adaptive layered technique. Through dividing the fluid into two kinds of particles, the modified methods use fewer particles and get lower computation complexity than the original methods. In addition, the modified methods even have better surface detail for using small particles to cover the surface. By adding the adaptive layered model, the modified methods can get lower computation complexity but the detail of surface is unchanged or even better than that of the original ones.

fluid simulation; adaptivity; hierarchy

2015-12-16；

2016-03-22

国家自然科学基金资助项目(61272024)

梁欣鑫(1991-),男,广西玉林人,合肥工业大学硕士生; 朱晓临(1964-),男,安徽池州人,博士,合肥工业大学教授,硕士生导师.

10.3969/j.issn.1003-5060.2017.02.026

TP301.6

A

1003-5060(2017)02-0277-07