劉旭，班曉娟，楊鳴遠(yuǎn)，賀亮亮

（北京科技大學(xué)計算機與通信工程學(xué)院，北京100083）

一種SPH流體仿真邊界校正方法

劉旭，班曉娟，楊鳴遠(yuǎn)，賀亮亮

（北京科技大學(xué)計算機與通信工程學(xué)院，北京100083）

使用光滑粒子流體動力學(xué)方法進(jìn)行流體仿真，并提出一種邊界校正方法。使用快速泊松盤采樣算法對容器邊界進(jìn)行采樣，生成邊界粒子，對邊界粒子質(zhì)量進(jìn)行差值估算，計算邊界粒子對流體粒子的作用力，以此來仿真流體與邊界的相互作用。該方法可避免穿刺、滯留等現(xiàn)象的發(fā)生。通過實驗驗證了該算法的正確性。

流體仿真；光滑粒子流體動力學(xué)；邊界校正

計算機圖形學(xué)的一個重要目標(biāo)就是逼真地再現(xiàn)自然界各物體的運動及形態(tài)。流體作為自然界的重要組成部分，其運動復(fù)雜而多變，流體仿真一直以來都是圖形學(xué)工作者研究的難點與熱點。光滑粒子流體動力學(xué)((smoothed particle hydrodynamics, SPH)[1-2]方法最早用于解決三維開放空間的天體物理學(xué)問題，由于其去網(wǎng)格化對于仿真流體的潑濺、卷浪等現(xiàn)象具有天然優(yōu)越性，因此近年來越來越多的研究者開始研究流體仿真SPH方法。Müller等[3]首次開創(chuàng)性地提出了SPH流體仿真的算法框架。Becker和Teschner[4]對不可壓縮流體仿真進(jìn)行研究，提出了弱可壓縮SPH方法((weakly compressible SPH,WCSPH)。他使用Tait方程取代氣體狀態(tài)方程，利用對稱的壓力公式并提出了一種新的表面張力建模方法。Solenthaler和Pajarola[5]提出了一種預(yù)測矯正 SPH法(predictive-corrective incompressible SPH,PCISPH)，該方法可大大增加仿真的時間步長，從而提高了仿真系統(tǒng)的整體效率。邊界條件是SPH流體仿真的一個重要研究點，許多研究人員將剛體邊界使用粒子來表示，Monaghan等[6]通過給與剛體邊界粒子相鄰的流體粒子施加邊界力來解決。Ihmsen等[7]將流體粒子密度外推至邊界粒子，然后對流體重新計算。Schechter和Bridson[8]對其改進(jìn)，使用一種ghostSPH方法，得到了令人滿意的結(jié)果。

SPH方法流體仿真計算量及大，許多研究人員對于SPH加速問題進(jìn)行了研究。文獻(xiàn)[9]中使用自適應(yīng)大小粒子方法提高效率。Solenthaler和Gross[10]根據(jù)細(xì)節(jié)需要使用兩種大小的粒子對流體進(jìn)行模擬，減少了粒子數(shù)量從而提高了效率。文獻(xiàn)[11]提出了一種基于背景網(wǎng)格的并行化算法。文獻(xiàn)[12]針對大規(guī)模場景不宜使用背景網(wǎng)格的缺點，使用hash表來存儲粒子。文獻(xiàn)[13]考慮到Cache命中率問題，針對背景網(wǎng)格使用粒子Z排序方法，針對大規(guī)模場景的hash存儲提出了緊湊hash策略，該方法可提升2～3倍的運行效率。很多研究者也對使用GPU并行計算SPH進(jìn)行了研究[14-15]。更多的流體仿真知識可參考文獻(xiàn)[1]，更多SPH知識可參考文獻(xiàn)[2]。

流體的運動往往伴隨與其他物體的交互，因此使用SPH方法對流體和其他物體的交互進(jìn)行模擬也是一個研究重點。Mao和Yang[16]對不可混溶流體交互進(jìn)行研究，提出了一種不可混溶流體交互仿真方法。Du等[17]使用三角形網(wǎng)格表示可變形體，并通過持續(xù)碰撞檢測實現(xiàn)了一種流體與可變形體交互方法。Akinci等[18]提出了一種通用的流體剛體仿真方法，該方法使用一種矯正的密度評估方法來計算流體和剛體的相互作用力。SPH方法使用粒子對流體建模，在動畫制作中需要對流體自由表面進(jìn)行表面重構(gòu)以進(jìn)行后期渲染。Mathieu和Cani-Gascuel[19]使用Blobby球方法生成的表面，并施加“表面張力”來平滑表面，從而克服了不自然的凹凸效果；Müller等[3]用SPH方法代替三維高斯函數(shù)對空間標(biāo)量場進(jìn)行疊加；Zhu和Bridson[20]在Blobby球方法的基礎(chǔ)上對粒子附近的密度變化進(jìn)行補償，并進(jìn)一步對空間標(biāo)量場進(jìn)行了一次平滑，從而得到更光滑的表面。Yu和Turk[21]提出了一種用各向異性核函數(shù)(anisotropickernel)構(gòu)造表面的方法。陳沸鑌等[22-23]則進(jìn)一步用SPH方法模擬了剛體的破裂等現(xiàn)象。

但是由于剛體建模與流體建模不同，剛體表面由網(wǎng)格構(gòu)成且粒子分布較為稀疏。而基于SPH的流體模擬，流體由粒子構(gòu)成。這樣剛體與流體的接觸面上兩種粒子的數(shù)量差極大，不但會影響雙方受力的計算，還會產(chǎn)生液體粒子穿進(jìn)或穿過剛體的現(xiàn)象。因此本文將針對復(fù)雜場景流體運動進(jìn)行SPH仿真，使用泊松盤對場景邊界進(jìn)行采樣，此處將剛體視為場景的一部分，構(gòu)建場景邊界粒子，并估算其質(zhì)量，使用校正的公式來計算流體粒子與邊界的交互作用。

1 SPH流體仿真

1.1 SPH方法

SPH方法的核心思想是以離散化粒子的形式來表征連續(xù)的場，并對場物理量使用積分近似的方式進(jìn)行計算。場量A（x）可用下式近似計算：

其中，mj、ρj分別表示粒子質(zhì)量和密度；W （x-xj,h）為積分光滑核函數(shù)，一般滿足3個條件：歸一化條件、緊支持性和Dirac函數(shù)性；h為其積分半徑。核函數(shù)可取滿足條件的多個函數(shù)，文獻(xiàn)[2]中使用B-樣條函數(shù)，文獻(xiàn)[3]使用Spiky核函數(shù)并且在計算粘性力時設(shè)計了專門的核函數(shù)，此外，高斯函數(shù)也經(jīng)常被用做光滑核函數(shù)。本文使用三次B-樣條函數(shù)作為核函數(shù)，取h=4r，r為粒子半徑。

密度：流體運動滿足質(zhì)量守恒，即滿足連續(xù)方程：

其中，v表示粒子速度。在SPH流體仿真中，粒子密度使用插值進(jìn)行計算，由式(1)可知粒子密度計算公式：

該公式可保證系統(tǒng)質(zhì)量守恒。由該式可知粒子密度僅與當(dāng)前狀態(tài)周圍粒子分布有關(guān)，密度的計算是SPH流體仿真的基礎(chǔ)。只有當(dāng)周圍粒子呈現(xiàn)各向球狀均勻分布時，使用上式可得近似結(jié)果，而在流體邊界處粒子分布不滿足該條件，將導(dǎo)致粒子密度值偏小，導(dǎo)致邊界問題，這將在后一章進(jìn)一步討論。

壓力：流體運動滿足動量守恒，流體運動動量方程如下：

其中，P表示壓強，g表示外力場。通過推導(dǎo)[2]，可使用下式計算粒子所受壓力：

由該式可知，相鄰粒子相互作用力是對稱的，且能保證線動量和角動量守恒。

狀態(tài)方程：通過式(3)求得粒子密度之后，需要對粒子壓強進(jìn)行計算，以求得粒子壓力。文獻(xiàn)[3]使用氣體狀態(tài)方程P=k(ρ-ρ0)來計算粒子壓強，其中，k為常數(shù)， ρ0為常量密度。該方法將導(dǎo)致比較大的可壓縮性。本文使用文獻(xiàn)[4]中的Tait方程進(jìn)行密度計算：

1.2 人工粘度

為了保證數(shù)值計算穩(wěn)定性，在SPH流體仿真中通常加入粘性力，文獻(xiàn)[3]中對流體運動控制NS方程中的粘力項進(jìn)行推導(dǎo)得到粘力計算公式：

該式計算的粘力并不對稱。文獻(xiàn)[2]中提出了一種人工粘性力：

其中，

式中，α是可調(diào)參數(shù)，ε=0.01，εh2的引入是為避免xij2=0。

本文使用式(8)計算粘性力，該方法產(chǎn)生的人工粘力是對稱的，即可以保證數(shù)字計算的穩(wěn)定性，同時還可以模擬流體運動的粘性效果。

1.3 流體表面張力

表面張力是流體的一大重要特性，文獻(xiàn)[3-4]都對表面張力進(jìn)行了研究，本文采用文獻(xiàn)[4]的方法。該方法基于這樣一個事實，表面張力的微觀機理是分子間引力，通過在相鄰粒子間加入引力來對表面張力進(jìn)行模擬，如下式：

其中，κ為可調(diào)參數(shù)，用來控制表面張力的強弱。在流體內(nèi)部，由于周圍粒子分布均勻，各周圍粒子對粒子i產(chǎn)生的表面張力合力項為0，在流體表層則不為0。

1.4 時間積分

仿真的每一步，先根據(jù)粒子分布，計算出各粒子密度ρi，然后依據(jù)式(6)求得粒子壓強 Pi，最后分別計算粒子壓力FiP、粘性力Fiv和表面張力Fie，則每個粒子加速度為：

其中，g為外力場，通常為重力。之后需要對粒子進(jìn)行時間積分，常用的時間積分方法有：蛙跳法(leap-frog,LF)、預(yù)測校正法(predictor-corrector)、龍格-庫塔法(Runge-Kutta,RK)等，本文采用簡單易用的蛙跳法進(jìn)行實現(xiàn)，可參考文獻(xiàn)[25]。積分時間步長通過CFL(courant-friedrichs-lewy)條件[4]確定。

2 基于質(zhì)量校正的邊界處理

流體的邊界包括流體和空氣的交界面以及流體和容器壁的交界面，本文只討論流體和容器壁交界面的處理。為了簡化處理過程，將容器看成剛體。本文結(jié)合文獻(xiàn)[8]和文獻(xiàn)[11]的方法，通過引入邊界粒子來實現(xiàn)容器壁和流體的相互作用。首先將邊界采樣為邊界粒子，進(jìn)而通過計算流體粒子和邊界粒子之間的相互作用來計算容器壁和流體的相互作用。

2.1 邊界泊松采樣

為得到邊界粒子，需要對以網(wǎng)格形式存在的容器模型進(jìn)行采樣，本文對容器表面按泊松盤分布（Poisson disk distribution，一種藍(lán)色噪聲分布）進(jìn)行采樣，具體參考文獻(xiàn)[26]。注意，在整個仿真過程中對每個容器或剛體只進(jìn)行一次初始采樣，該采樣結(jié)果將在后續(xù)整個仿真過程中使用，因此，對容器或剛體的采樣并不會帶來額外的開銷。表面采樣結(jié)果如圖1所示。

圖1剛體表面采樣結(jié)果

2.2 邊界粒子校正計算

考慮到邊界粒子的影響，流體粒子的密度計算式(2)需要加入邊界粒子的加權(quán)求和結(jié)果：

式中，用fi來表示流體粒子i，bk表示邊界粒子k，用j迭代粒子if的所有流體鄰居粒子，用k迭代所有邊界鄰居粒子。

文獻(xiàn)[11]對上式進(jìn)行了改進(jìn)?？紤]到邊界粒子質(zhì)量設(shè)置不合理或分布不均勻，將導(dǎo)致計算錯誤，使用下式估算邊界粒子質(zhì)量：

其中，ρ0表示流體的常量密度，為邊界粒子所代表邊界區(qū)域體積的估算值。使用Ψbi（ρ0）取代邊界粒子質(zhì)量進(jìn)行計算可提高穩(wěn)定性。

因此，式(12)可改為：

容器壁和流體最重要的相互作用就是壓力，用下式計算邊界粒子對流體粒子的壓力產(chǎn)生的加速度：

式中，pfi＞0時取k=2。pfi＜0時，邊界粒子和流體粒子相互吸引，可對k進(jìn)行調(diào)整(0≤k≤2)以實現(xiàn)不同吸附效果，本文中取k=1。

為模擬容器壁與流體之間的摩擦力或?qū)崿F(xiàn)流體和剛體的交互，需要計算邊界粒子和流體粒子摩擦力，摩擦力的計算借鑒了人工粘度公式(8)：

式中的Πik同式(9)。

2.3 算法流程

本文流體仿真算法流程如下：

3 實驗結(jié)果

本文通過實驗驗證了算法的有效性，程序的運行平臺為Inteli5-3470（四核，3.20GHz，6MB Cache）、8GB內(nèi)存。仿真算法和表面重構(gòu)算法用C++語言采用多線程技術(shù)實現(xiàn)，仿真算法中的臨近粒子搜索使用空間背景網(wǎng)格進(jìn)行哈希查找。表面重構(gòu)采用文獻(xiàn)[15]的方法，使用各向異性核函數(shù)構(gòu)造顏色場，然后使用步進(jìn)立方體算法重構(gòu)表面，其中矩陣奇異值分解采用的是NIST的JAMA和TNT開源數(shù)學(xué)庫。使用OpenGL三維圖形庫實時顯示仿真和表面構(gòu)造的結(jié)果，并使用OpenCV庫錄制視頻。為了對復(fù)雜容器和場景進(jìn)行建模以及制作動畫，本文使用Blender軟件，后期高質(zhì)量流體效果的渲染使用光線跟蹤引擎Mitsuba。

為了驗證邊界處理算法的正確性，本文進(jìn)行了噴泉仿真實驗與城市洪澇實驗。噴泉實驗中流體粒子數(shù)約為100～200k，邊界粒子數(shù)150k，粘度系數(shù)α=0.08，表面張力系數(shù)κ=0.01，時間步長為0.0002s。圖2為噴泉實驗結(jié)果。

圖2 噴泉實驗結(jié)果

城市洪澇實驗中流體粒子數(shù)約為50～200k，邊界粒子數(shù)80k，粘度系數(shù)α=0.08，表面張力系數(shù)κ=0.005，時間步長為0.0005s。圖3為城市洪澇實驗結(jié)果。

手?jǐn)噭铀膶嶒炛辛黧w粒子數(shù)約為1.17M，邊界粒子數(shù)63k，粘度系數(shù)α=0.05，表面張力系數(shù)κ=0.005，時間步長為0.000134s，如圖4所示。

圖3 城市洪澇實驗不同時刻渲染結(jié)果

圖4 手?jǐn)噭铀哪M

4結(jié)論

本文通過實驗驗證了所使用的SPH流體仿真和邊界處理算法的有效性，并制作了噴泉動畫，實驗結(jié)果逼真地反應(yīng)了潑濺、吸附、表面張力、粘性等流體特征。但是算法計算量較大，今后工作將對如何提高算法效率進(jìn)行研究，此外還將研究流體剛體交互計算方法。

[1] Bridson R, Müller-fischer M. Fluid simulation: SIGGRAPH 2007 course notes video filesassociated w ith this course are available from the citation page[C]//ACM SIGGRAPH 2007Courses.New York,USA,2007:1-81.

[2]Monaghan J J.Smoothed particle hydrodynam ics[J]. Reportson Progress in Physics,2005,68(8):1703-1759.

[3]Müller M,Charypar D,Gross M.Particle-based fluid simulation for interactive applications [C]//ACM SIGGRAPH/Eurographics Symposium on Computer Animation.San Diego,CA,USA,2003:154-159.

[4]Becker M,Teschner M.Weakly compressible SPH for free surface flows[C]//ACM SIGGRAPH/Eurographics symposium on Computer animation.Sw itzerland,2007: 209-217.

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

[6]Monaghan JJ,Kos A,Issa N.Fluidmotion generated by impact[J].Journal of Waterway Port Coastal and Ocean Engineering,2003,129(6):250-259.

[7]Ihmsen M,Akinci N,Gissler M.Boundary handling and adaptive time-stepping for PCISPH[C]//Workshop in Virtual Reality Interactions and Physical Simulation, VRIPHYS(2010).Darmstadt,Germany,2010:79-88.

[8] Schechter H,Bridson R.Ghost SPH for animating water[J]. ACM Transcationson Graphics,2012,31(4):61:1-61:8.

[9]Adams B,Pauly M,Keiser R,et al.Adaptively sampled particle fluids[J].ACM Transactions on Graphics,2007, 26:48:1-48:8

[10]Solenthaler B,GrossM.Two-scale particle simulation[J]. ACM Transactionson Graphics,2011,30(4):81:1-81:8.

[11]Kalojanov J,Slusallek P.A parallel algorithm for construction of uniform grids [C]//In HPG '09: Proceedings of the 1st ACM Conference on High PerformanceGraphics.New York,USA,2009:23-28.

[12]Bell N,Yu Yizhou,Mucha P J.Particle-based simulation of granularmaterials[C]//ACM SIGGRAPH/Eurographics Symposium on Computer Animation,New York,USA, 2005:77-86.

[13]Ihmsen M,Akinci N,Becker M,et al.A parallel SPH implementation on multicore CPUs [J].Computer Graphics Forum,2011,30(1):99-112.

[14]Harada T,Koshizuka S,Kawaguchi Y.Smoothed particle hydrodynamics on GPUs [C]//In Proceedings of Computer Graphics International.Petropolis,RJ,Brazil, 2007:63-70.

[15]Goswam i P,Schlegel P,Solenthaler B,et al.Interactive SPH simulation and rendering on the GPU [C]// Proceedings of the 2010 ACM SIGGRAPH/Eurographics Symposium on Computer Animation.Madrid,Spain, 2010:55-64.

[16]Mao Hai,Yang Y.Fluid-fluid collision for particle-based immiscible fluid-fluid collision[C]//GI'06 Proceedingsof Graphics Interface 2006.Toronto Canada,2006:49-55.

[17]Du Peng,Tang M in,Meng Chang,et al.A fluid/cloth coupling method for high velocity collision simulation[C]// Proceedings of the 11th ACM SIGGRAPH International Conference on Virtual-Reality Continuum and its Applications in Industry.Singapore,2012:309-314.

[18]AkinciN,Ihmsen M,AkinciG,etal.Versatile rigid-fluid coupling for incompressible SPH[J].ACM Transactions on Graphics,2012,31:62:1-62:8.

[19]Mathieu D,Cani-Gascuel M.Active implicit surface for animation [C]//Graphics Interface.Vancouver,BC, Canada,1998:143-150.

[20]Zhu Yongning,Bridson R.Animating sand as a fluid[J]. ACM Transactionson Graphics,2005,24(3):965-972.

[21]Yu Jihun, Turk G. Reconstructing surfaces of particle-based fluids using anisotropic kernels[J].ACM Transactionson Graphics,2013,32(5):1-5,12.

[22]陳沸鑌,王長波,謝步瀛.可控性的雙層粒子剛體脆性破裂模擬動畫[J].圖學(xué)學(xué)報,2015,36(1):111-116.

[23]陳沸鑌,謝步瀛,冉修遠(yuǎn).基于細(xì)分粒子的剛體破裂動畫[J].圖學(xué)學(xué)報,2014,35(5):669-675.

[24]Liu G R,Liu M B.光滑粒子流體動力學(xué)——一種無網(wǎng)格粒子法[M].韓 旭,楊 剛,強洪夫,譯.長沙:湖南大學(xué)出版社,2005:1-433.

[25]Desbrun M,CaniM.Smoothed particles:a new paradigm for animating highly deformable bodies [C]//6th Eurographics Workshop on Computer Animation and Simulation'96.Poitiers,France,1996:61-76.

[26]Bridson R.Fast Poisson disk sampling in arbitrary dimensions [C]//In ACM SIGGRAPH Technical Sketches2007.San Diego,CA,USA,2007:1.

A Boundary Correction Method of SPH Fluid Simulation

Liu Xu, Ban Xiaojuan, Yang Mingyuan, He Liangliang

(Department of Computer Science and Technology, University of Science and Technology Beijing, Beijing 100083, China)

This article presents a boundary correction method in fluid simulation based on smoothed particle hydrodynamics. We sample boundaries as boundary particles with fast Poisson disk algorithm,interpolate mass of boundary particles, and calculate force between fluid particles and boundaryparticles as simulation of fluid-boundary interaction. The presented method can also avoid penetration and holdup. Experiments demonstrated conduct experiments to prove the validity of our method.

fluid simulation; smoothed particle hydrodynamics; boundary correction

TP391

A

2095-302X(2015)03-0462-06

2014-12-26；定稿日期：2015-02-08

國家自然科學(xué)基金資助項目(61272357，61300074)

劉旭(1987-)，男，北京人，博士研究生。主要研究方向為計算機圖形學(xué)、流體模擬及表面重構(gòu)。E-mail：liuxu.ustb@gmail.com

班曉娟(1970-)，女，遼寧朝陽人，教授，博士。主要研究方向為人工智能、計算機動畫。E-mail：Banxj@ustb.edu.cn