基于IIF模型的流體表面張力模擬

董蘭芳， 章 恒， 陳駿雄

(1. 中國科學(xué)技術(shù)大學(xué)計算機科學(xué)與技術(shù)學(xué)院，安徽 合肥 230027；2. 浙江大學(xué)流體動力與機電系統(tǒng)國家重點實驗室，浙江 杭州 310027)

基于IIF模型的流體表面張力模擬

董蘭芳1,2， 章 恒1， 陳駿雄1

(1. 中國科學(xué)技術(shù)大學(xué)計算機科學(xué)與技術(shù)學(xué)院，安徽 合肥 230027；2. 浙江大學(xué)流體動力與機電系統(tǒng)國家重點實驗室，浙江 杭州 310027)

基于光滑粒子動力學(xué)(SPH)方法模擬流體時流體表面張力的作用在固液、氣液交界處不可忽視，其影響模擬的準確性和視覺真實感。目前已有的表面張力模型如連續(xù)表面力(CSF)模型、粒子間相互作用力(IIF)模型都存在各自的缺陷。針對 IIF模型模擬表面張力時所產(chǎn)生的粒子非物理聚集、流體表面形狀不規(guī)則等現(xiàn)象，采用基于類Lennard-Jones勢函數(shù)的分子間聚斥力對表面張力建模，并定義了基于法向差的SPH張力修正項以解決IIF模型不能保證流體表面面積最小化問題。實驗結(jié)果表明，該方法能夠穩(wěn)定和正確地模擬兩相交界處的表面張力的效果。

計算機圖形學(xué)；流體模擬；光滑粒子動力學(xué)方法；表面張力；粒子間相互作用力模型

表面張力作為流體的一種常見和重要的物理性質(zhì)，隨著對流體細節(jié)和模擬精度要求的增強，在諸如計算機圖形學(xué)領(lǐng)域和工程數(shù)值計算領(lǐng)域?qū)Ρ砻鎻埩Φ哪M也越來越重視。從微觀層面說，表面張力的產(chǎn)生是由于表面區(qū)域的水分子受到各種方向分子間作用力的合力不為零，指向流體內(nèi)部。從宏觀表現(xiàn)上看，根據(jù)拉普拉斯定律，表面張力的作用就是最小化表面區(qū)域面積，如自然界中液滴的形成、聚合和分裂。

光滑粒子動力學(xué)方法(smoothed particle hydrodynamics, SPH)作為一種典型的拉格朗日流體模擬方法，由于具備相較于傳統(tǒng)基于網(wǎng)格方法的易編程、方便處理復(fù)雜邊界和大形變問題等優(yōu)點而被廣泛應(yīng)用于流體模擬領(lǐng)域[1-3]。目前基于SPH的表面張力模擬方法主要分為兩種：基于連續(xù)表面力模型(continuum surface force, CSF)和基于粒子間相互作用力模型(inter-particle interaction force, IIF)。CSF模型是將作用在流體表面上的力轉(zhuǎn)換為周圍流體體積力，最早由Morris[4]在2000年提出，并與基于網(wǎng)格的流體體積函數(shù) (volume of fluid, VOF)方法進行了對比，得到可信的結(jié)果。在Morris的基礎(chǔ)上，Müller等[5]采用顏色域函數(shù)對其進行簡化，Hu和Adams[6]則通過配置新的粘性項，解決了速度和剪切應(yīng)力在多相交界處的連續(xù)性問題。2011年強洪夫等[7]采用SPH方法中解決邊界核函數(shù)插值問題的校正平滑粒子方法(corrective smoothed particle method, CSPM)對CSF模型中表面張力計算進行修正，并對二維溶液中油滴的碰撞破碎過程進行了模擬。如文獻[7]所述，CSF模型模擬表面張力時在尖角、邊界粒子缺失嚴重的部位曲率計算誤差很大，需要花費額外的計算代價來進行修正。另外，CSF模型以一種非對稱的形式將力施加到SPH粒子上，這使得流體的動量不再守恒。

IIF模型是一種基于分子間相互作用力，從微觀尺度上模擬表面張力產(chǎn)生的模型。在流體內(nèi)部，分子受到的周圍領(lǐng)域內(nèi)分子在各種方向上作用力的合力為零，與此同時，處于流體表面的分子僅受到流體內(nèi)部分子的作用力而導(dǎo)致合力非零，宏觀表現(xiàn)為表面張力。SPH方法模擬表面張力時，將流體用粒子表示，應(yīng)用IIF模型思想則對每對粒子施加一對對稱的作用力[8-9]，在粒子分布近似均勻的前提下，對整個流體系統(tǒng)施加這一作用力可近似模擬表面張力效果。文獻[8]和[9]的不同之處在于采用不同的分子間作用力模型。IIF模型克服了CSF模型的缺點，不再需要計算曲率，因此編程簡單，計算效率高。同時由于對稱形式作用力的引入保證流體的動量守恒。然而實驗證明，僅僅引入粒子間相互作用力不能保證表面區(qū)域面積趨向最小，因為流體平衡狀態(tài)可以在非最小表面面積狀態(tài)下達到。

本文所基于IIF模型的表面張力模型，采用基于距離的類 Lennard-Jones勢函數(shù)為粒子間作用力建模。同時，為了解決IIF模型不保證表面面積最小問題，利用文獻[5]中的一些結(jié)論引入了一個張力修正項，通過對立方體液滴在去除重力下形變過程的模擬，得到了與文獻[7]中修正后的CSF模型幾乎一致的結(jié)果，證明了該方法的正確性。通過對碰撞破碎過程對比和多液滴合并的模擬證明該方法的有效性。

1 流體模擬框架

1.1 SPH基本形式

SPH方法中用粒子來承載諸如密度、質(zhì)量等各個物理量，連續(xù)的流體場被離散成大量獨立的粒子，流體連續(xù)的物理量根據(jù)周圍流體粒子相對應(yīng)的物理量插值得到。SPH方法中的核函數(shù)用于確定空間中周圍粒子對中心粒子插值的權(quán)重，一般來說，核函數(shù)的選擇使得離中心粒子越近的粒子插值權(quán)重越大。SPH基本離散化的形式[10]為：

其中，W為插值所用的光滑核函數(shù)；h為其核半徑；ri和rj分別為粒子i、j的位置；Aj為粒子j所對應(yīng)的物理量；Vj為粒子j的體積。

以當(dāng)前粒子為中心，h為半徑的球體范圍稱為當(dāng)前粒子的支持域。直觀上，SPH方法通過枚舉支持域內(nèi)所有粒子以用于插值，從而得到當(dāng)前粒子的相應(yīng)屬性值。

1.2 SPH流體運動方程

流體運動用經(jīng)典的Navier-Stokes方程來刻畫，如式(2)所示：

其中，ρ為流體的密度；u為速度矢量；p為壓強；μ為粘力系數(shù)；f為合外力項諸如重力和浮力等。

基于文獻[6]用SPH方法離散化Navier-Stokes方程，先將密度項根據(jù)粒子的分布計算出來，再用于計算壓力項和粘力加速度項，如式(3)～(5)所示。

為了從密度求出壓強，采用文獻[9]提出的弱可壓縮法(weakly compressible SPH，WCSPH)，模擬時間步長根據(jù)CFL條件[10]來估計。對于在流體與空氣交界處表面由于缺少鄰居粒子，導(dǎo)致密度估計失真而產(chǎn)生的拉伸不穩(wěn)定現(xiàn)象，采用Shepard filter來對流體密度進行平滑。對于在流體與固體交界處為了避免 sticking particles問題，采用文獻[11]中利用單層固體邊界粒子對周圍鄰居流體粒子進行密度校正的方法解決。

2 修正的IIF表面張力模型

流體的表面張力由分子間的內(nèi)聚力產(chǎn)生，然而，在采用單純的IIF模型模擬SPH流體是在宏觀層面上考慮粒子間的相互作用力，由于對稱的粒子作用力的引入，整個流體系統(tǒng)的動量是守恒的，流體系統(tǒng)的平衡態(tài)可以在流體的任意形態(tài)達到，這樣導(dǎo)致的結(jié)果就是流體表面區(qū)域面積的最小化難以保證，需要通過其他手段來進行修正。本文中表面張力模型中采用一個類Lennard-Jones勢函數(shù)為粒子間的聚斥力建模，同時考慮使得表面區(qū)域最小化的表面力項。

2.1 粒子間聚斥力模型

IIF模型中為粒子間施加一個基于粒子間距離的作用力，基本形式如式(6)所示，F(xiàn)cohesioni＜-j描述粒子j對粒子i的作用力。其中，mi和mj是粒子i和j的質(zhì)量，L(d)是一個核函數(shù)控制粒子間作用力的正負(分別對應(yīng)著排斥力和吸引力)，文獻[8]最早提出基于一個類余弦函數(shù)的粒子間聚斥力，但該函數(shù)在粒子距離 d趨于0時，函數(shù)值也趨于0，意味著粒子間的排斥力隨著距離的減小而減小，這會導(dǎo)致在模擬弱可壓縮流體時出現(xiàn)非物理聚集現(xiàn)象。文獻[9]采用一個SPH核函數(shù)來為粒子間吸引力建模，由于缺少排斥力項，該模型也會導(dǎo)致非物理聚集現(xiàn)象的發(fā)生。本文采用文獻[12]中模擬流體和可變形固體交互時流固間作用力所使用的類Lennard-Jones勢函數(shù)作用核函數(shù)：

其中，h為作用范圍；k為控制表面張力影響程度的系數(shù)；d0為靜止距離；d為當(dāng)前兩粒子間距離。

式(7)中不同作用范圍下L(d)的值隨粒子間距離d的變化而變化(圖1)，其中系數(shù)k值取0.001，d0值取4×10-3mm，h分別取1.5 d0, 2 d0, 2.5 d0。

圖1 不同作用范圍下L(d)隨粒子間距的變化

由圖1可見，當(dāng)d=d0時，粒子間作用力為0，當(dāng) h＞d＞d0時,函數(shù)值為負，表現(xiàn)出吸引力，當(dāng)0＜d＜d0時，函數(shù)值為正，變現(xiàn)為排斥力。另外,當(dāng)d趨于 0時，函數(shù)值趨于 k，這樣避免了文獻[8] 和[9]中的非物理聚集現(xiàn)象發(fā)生。

2.2 表面張力修正項

基于類Lennard-Jones勢函數(shù)的粒子間聚斥力可以有效地解決粒子非物理聚集的問題，并且保證動量守恒，但不能保證表面面積最小化。為了解決這個問題，本文引入一個表面張力的修正項：

其中，mi為當(dāng)前粒子i的質(zhì)量；ni和nj分別為粒子i處的法向和粒子j處的法向；λ為修正項系數(shù)。

該修正項的思想受文獻[13]中的分析啟發(fā)，從離散的觀點來說，表面張力趨向于將一個局部區(qū)域內(nèi)所有離散采樣點的法向變得一致。直觀地說，即撫平該區(qū)域的曲率。對于粒子 i，的意義在于作用于粒子i，使得i和j之間的法向差變小，推廣到整個計算域中所有流體粒子，項的施加將使得曲率大的地方(即流體表面)的法向差變小，即減小曲率，從而最小化表面面積。推廣到整個計算域上，粒子i的表面張力修正項用SPH形式表示有以公式：

其中，F(xiàn)fixi為粒子i的表面張力修正項。

對于一對粒子i和j而言：

粒子間的表面張力修正項也是對稱的，因此修正項的引入不破壞動量守恒條件。關(guān)于粒子法向的計算，采用文獻[5]的方法，即先為整個計算域定義一個顏色域函數(shù)，然后計算顏色函數(shù)的梯度，得到的梯度函數(shù)就是計算域上每處的未歸一化法向。用SPH的形式表示如下：

在CSF模型中，通過每個粒子處法向的散度來估計該處的曲率，然后定義一個大小與曲率相關(guān)，方向為法向的力，即該粒子處所受到的表面張力。整個過程需要對法向進行歸一化，而在流體內(nèi)部，對法向歸一化可能導(dǎo)致流體內(nèi)部的法向計算不準確，從而需要額外處理來修正。在本文模型中，由于避免了直接計算粒子處的曲率，從而也無需對法向進行歸一化，因此也就不需要進行額外的處理。

3 算法實現(xiàn)和參數(shù)選擇

本文模型中計算每個粒子的表面張力時需要考慮周圍鄰居粒子的貢獻，所以對表面張力的模擬是耦合在流體模擬框架中的。單獨來看，具體步驟如圖2所示。

圖2 表面張力模擬流程

關(guān)于 SPH流體的模擬，采用文獻[6]和[14]中所設(shè)計的幾個多項式核函數(shù)，針對計算不同物理量選取不同的核函數(shù)。SPH的光滑核半徑取值為1.5R，R是 SPH粒子的直徑，粒子半徑的選取與采樣流體所使用的粒子數(shù)有關(guān)，壓強的計算采用文獻[9]中的WCSPH方法，時間步長的選擇根據(jù)CFL條件[10]來估計，在實驗中時間步長設(shè)為0.005 s。數(shù)值積分采用簡單的leap-frog模式[15]。關(guān)于表面張力模型中涉及的2個主要參數(shù)：類Lennard-Jones勢函數(shù)中的聚斥力縮放系數(shù) k和表面張力修正項中修正系數(shù)λ，這兩個模型參數(shù)都是經(jīng)驗值。實驗表明，兩個參數(shù)的取值應(yīng)在同一數(shù)量級上，在實驗中將其保持與流體的粘度系數(shù)一致(0.001)，能夠在維持流體粒子的1%密度震蕩幅度的同時取得良好的可視化效果。

4 實驗算例和結(jié)果

實驗配備 2.8 GHz Intel Core i5-2300 四核CPU和GeForce GTS 450顯卡以及4 G內(nèi)存的硬件機器，流體粒子形式的渲染基于OpenGL平臺實時繪制。

4.1 立方體液體形變?yōu)橐?guī)則球體

算例1. 模擬了0.5×0.5×0.5(m3)的水柱在無重力的情況下由表面張力驅(qū)動，形變成規(guī)則球體的過程，其中流體粒子4 913個。在無重力等外力的情況下，流體穩(wěn)定狀態(tài)下的形狀主要由表面張力決定，表面張力使得流體各個地方的曲率趨于一致、表面面積趨于最小化，立方體液塊在表面張力的作用下將收縮震蕩，最終由于流體的粘性作用而穩(wěn)定成一個規(guī)則的球體。該算例的目的在于證明該表面張力算法的正確性。圖 3是液塊形變模擬過程的部分關(guān)鍵幀截圖信息(流體為三維表示，觀察方向為水平方向)。

圖3 立方體液塊形變過程

首先，液塊由立方體先是收縮，再迅速形變?yōu)橐粋€球體(圖 3(a)～(d))；然后，液塊由一個球體震蕩為一個具有6個尖角的極限狀態(tài)(圖3(e)～(h))；最后，流體由于粘力對動能的耗散作用使得液塊形最終穩(wěn)定在球體狀態(tài)下(圖3(i)～(l))。該結(jié)果與文獻[7]中采用修正的CSF模型在二維情況下得到的數(shù)值模擬結(jié)果類似。

4.2 液塊碰撞破碎和液滴合并

算例2. 模擬了在體積為1 m×1 m×1 m的固體容器中，體積為0.25 m×0.25 m×0.25 m的液塊在重力的作用下下落到與容器底碰撞破碎形成多個液滴，以及液滴在表面張力的作用下合并的過程，并與不施加表面張力情況下做出對比。其中粒子數(shù)為810個。圖4給出液塊在碰撞破碎后施加表面張力與不施加表面張力模擬過程對比(流體為三維表示，觀察方向為從容器頂部到底部的垂直方向)。

圖4 碰撞破碎過程比較

圖 4(a)和(b)自左向右分別是施加表面張力和不施加表面張力的液塊自下落開始第60、70、80、90幀的截圖信息。從2個模擬過程的對比中可以看出，在施加表面張力的情況下液塊在與容器底部碰撞后破碎形成液滴，如圖3(a)所示，而不施加表面張力時破碎后成網(wǎng)狀失真，并且不能形成分離的圓形液滴，如圖3(e)所示。另外，隨著流體的運動，施加表面的張力的模擬過程中還伴有液滴間的碰撞合并現(xiàn)象。

液體運動趨于穩(wěn)定后的比較如圖5所示(觀察方向為水平方向)，圖 5(a)為施加表面張力后形成的一個單獨液滴，包含3層SPH粒子，圖5(b)則為沒有施加表面張力的情況下，流體粒子無法聚攏，從而平鋪在容器底部，只有1層SPH粒子。

圖5 流體靜止?fàn)顟B(tài)比較

圖6為8幅流體粒子位置截圖，顯示碰撞破碎產(chǎn)生的多個獨立的運動液滴在表面張力的作用下逐漸合并為一個單獨的穩(wěn)定的圓形液滴的模擬過程，該過程發(fā)生液體在與容器壁發(fā)生碰撞后，多個分散的液滴向容器中央聚攏。

圖6 多液滴合并過程

4.3 計算性能對比

為了測試本文方法的計算性能，使用方形液滴自然變化和油滴互溶兩個算例，選取運行40幀后添加表面張力的效果與傳統(tǒng)的 IIF模型和 CSF模型進行對比，時間為每幀的計算時間(表1)。

表1 計算速度對比(s)

本文方法由于在計算表面張力修正項過程中直接使用顏色域梯度求解粒子法向，相較于傳統(tǒng)CSF模型省去了法向歸一化以及相應(yīng)修正處理過程，在單幀計算速度上略有提高，模擬的場景越復(fù)雜，效果越明顯。與傳統(tǒng)IIF模型相比，本文方法在保證模擬速度的基礎(chǔ)上，通過改進聚斥力模型和修正表面張力進一步提高了模擬精度。

5 結(jié) 論

本文針對傳統(tǒng)的表面張力模擬所采用的IIF模型中存在的缺陷，定義一種最小化表面面積的表面張力修正項，結(jié)合基于類L-J勢函數(shù)的粒子間聚斥力，為SPH方法模擬流體時表面張力的模擬建模。通過模擬立方體液塊僅由表面張力驅(qū)動時形變?yōu)榍蝮w的算例驗證了本方法的正確性。最后通過比較液塊在施加表面張力和不施加情況下碰撞破碎過程的對比以及多液滴合并算例驗證了本文方法的有效性。

[1] Solenthaler B, Pajarola R. Predictive-corrective incompressible SPH [J]. Acm Transactions on Graphics, 2009, 28(3): 341-352.

[2] 溫嬋娟, 歐嘉蔚, 賈金原. GPU通用計算平臺上的SPH 流體模擬[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2010, 22(3): 406-411.

[3] Zhu Y N, Bridson R. Animating sand as a fluid [J]. Acm Transactions on Graphics, 2006, 24(3): 965-972.

[4] Morris J P. Simulating surface tension with smoothed particle hydrodynam ics [J]. International Journal for Numerical Methods in Fluids, 2000, 33: 333-353.

[5] Müller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications [C]//Preceedings of Symposium on Computer Animation. Aire?La?Ville: Eurographics Association Press, 2003: 154-159.

[6] Hu X Y, Adams N A. A multi-phase SPH method for macroscopic and mesoscopic flow s [J]. Journal of Computational Physics, 2006, 213(2): 844-861.

[7] 強洪夫, 陳福振, 高巍然. 基于 SPH方法的表面張力修正算法及其應(yīng)用[J]. 計算力學(xué)學(xué)報, 2011, 28(3): 37-42.

[8] Tartakovsky A, Meakin P. Modeling of surface tension and contact angles with smoothed particle hydrodynam ics [J]. Physical Review E, 2005, 72(2): 026301.

[9] Becker M, Teschner M. Weakly compressible SPH for free surface flows [C]//Proceedings of the 2007 ACM SIGGRAPH/ Eurographics symposium on Computer Animation. Aire?La?Ville: Eurographics Association Press, 2007: 209-217.

[10] Monaghan J J. Smooth particle hydrodynamics [J]. Reports on Progress in Physics, 2005, 68(8): 1703-1759.

[11] Akinci N, Ihmsen M, Akinci G, et al. Versatile rigid-fluid coupling for incompressible SPH [J]. ACM Transactions on Graphics (TOG), 2012, 31(4): 62.

[12] Müller M, Schirm S, Teschner M, et al. Interaction of fluids with deformable solids [J]. Computer Animation and Virtual Worlds, 2004, 15(3/4): 159-171.

[13] He X W, Wang H M, Zhang F J, et al. Robust Simulation of Small-Scale Thin Features in SPH-based Free Surface Flows [EB/OL]. [2015-08-16]. http://life.kunzhou.net/ 2013/SPHsurfacetension.pdf, 2013.

[14] Desbrun M, Gascuel M P. Smoothed particles: a new paradigm for animating highly deformable bodies [M]. Vienna: Springer Vienna, 1996: 61-76.

[15] Kelager M. Lagrangian fluid dynam ics using smoothed particle hydrodynamics [D]. Denmark: University of Copenhagen, 2006.

Surface Tension Simulation for SPH Fluid Based on IIF Model

Dong Lanfang1,2, Zhang Heng1, Chen Junxiong1

(1. Department of Computer Science and Technology, University of Science and Technology of China, Hefei Anhui 230027, China; 2. The State Key Laboratory of Fluid Power Transm ission and Control, Zhejiang University, Hangzhou Zhejiang 310027, China)

Surface tension effect at fluid-air and fluid-solid interfaces should not be ignored when simulating fluid based on smoothed particle hydrodynam ics (SPH) method, it affects the accuracy of the simulation and visual realism. The existing surface tension models such as the continuum surface force (CSF) model and the inter-particle interaction force (IIF) model both have their own flaws. This paper focus on the IIF model which can causes non-physical particle clustering and the irregular shape of fluid surface, and surface tension is modeled by an attraction-repulsion force which is formalized by a Lennard-Jones like potential function, and a SPH tension correction term which based on the normal difference between neighboring particles is introduced to m inim ize the surface area and smooth the surface. Experimental results show that the proposed method can simulate the surface tension effect accurately and stably.

computer graphics; fluid simulation; smoothed particle hydrodynam ics method; surface tension; inter-particle interaction force method

TP 39

10.11996/JG.j.2095-302X.2016030302

A

2095-302X(2016)03-0302-06

2015-09-28；定稿日期：2015-11-13

浙江大學(xué)流體動力與機電系統(tǒng)國家重點實驗室開放基金項目(GZKF-201318)

董蘭芳(1970-)，女，安徽合肥人，副教授，碩士。主要研究方向為科學(xué)計算可視化、計算機動畫、視頻圖像智能分析。E-m ial：lfdong@ustc.edu.cn