自適應分層定位流體仿真模型

梁欣鑫, 朱曉臨, 劉洋洋, 閆志偉

(合肥工業大學 數學學院,安徽 合肥 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