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

劉旭，班曉娟，楊鳴遠，賀亮亮

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

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

劉旭，班曉娟，楊鳴遠，賀亮亮

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

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

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

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

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

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

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

1 SPH流體仿真

1.1 SPH方法

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

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

密度：流體運動滿足質量守恒，即滿足連續方程：

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

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

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

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

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

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

1.2 人工粘度

為了保證數值計算穩定性，在SPH流體仿真中通常加入粘性力，文獻[3]中對流體運動控制NS方程中的粘力項進行推導得到粘力計算公式：

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

其中，

式中，α是可調參數，ε=0.01，εh2的引入是為避免xij2=0。

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

1.3 流體表面張力

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

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

1.4 時間積分

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

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

2 基于質量校正的邊界處理

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

2.1 邊界泊松采樣

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

圖1剛體表面采樣結果

2.2 邊界粒子校正計算

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

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

文獻[11]對上式進行了改進。考慮到邊界粒子質量設置不合理或分布不均勻，將導致計算錯誤，使用下式估算邊界粒子質量：

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

因此，式(12)可改為：

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

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

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

式中的Πik同式(9)。

2.3 算法流程

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

3 實驗結果

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

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

圖2 噴泉實驗結果

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

手攪動水的實驗中流體粒子數約為1.17M，邊界粒子數63k，粘度系數α=0.05，表面張力系數κ=0.005，時間步長為0.000134s，如圖4所示。

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

圖4 手攪動水的模擬

4結論

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

[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].圖學學報,2015,36(1):111-116.

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

[24]Liu G R,Liu M B.光滑粒子流體動力學——一種無網格粒子法[M].韓 旭,楊 剛,強洪夫,譯.長沙:湖南大學出版社,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

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

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

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