朱曉臨， 張義群， 郭清偉

(合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，安徽 合肥 230009)

基于物理的流體模擬方法已廣泛應(yīng)用于影視特效、生物醫(yī)療、災(zāi)難預(yù)防與營救等領(lǐng)域，而邊界處理及液體表面張力的計算一直是流體模擬研究的重要課題。

在光滑粒子流體動力學(xué)(smoothed particle hydrodynamics，SPH)方法中，傳統(tǒng)的流固邊界處理方法有3種：密度正則化法、邊界力法和虛粒子法。處理表面張力的方法為 2種：IIF方法(inter-particle interaction force)和CSF方法(continuum surface force)。

本文針對傳統(tǒng)的邊界力法和虛粒子法計算量大、對復(fù)雜邊界處理難以達(dá)到實時性的問題，提出改進(jìn)。

1997年，MORRIS等[1]提出了一種動態(tài)生成鏡像虛粒子的方法，該方法對于直邊界(邊界為平面)效果很好，但對于復(fù)雜邊界效果不佳；2002年，LIU等[2]將虛粒子法與邊界力法相結(jié)合，提高了SPH近似在邊界的精度，但需實時生成虛粒子，過程較為復(fù)雜；2006年，HU和ADAMS[3]將虛粒子看作另一種相，虛粒子關(guān)于邊界對稱鋪設(shè)插值點，其只對直邊界效果較好；2011年，MARRONE等[4]提出了一種只需在初始時刻生成虛粒子的方法，但其對復(fù)雜幾何形狀的場景，內(nèi)部插值點構(gòu)造復(fù)雜；2012年，SCHECHTER和 BRIDSON[5]提出 Ghost方法，用于流體自由表面和邊界條件的處理，虛粒子具有較穩(wěn)定的狀態(tài)分布，但需要更多的預(yù)處理和計算量；2012年，ADAMI等[6]將流體粒子的壓強(qiáng)賦值給虛粒子，利用狀態(tài)方程反推出虛粒子密度。該方法適用于具有復(fù)雜邊界的場景，但虛粒子的壓強(qiáng)可能為負(fù)值，導(dǎo)致流體穿透邊界；2015年，劉虎等[7]指出，該方法利用壓強(qiáng)反推密度是多余的。文獻(xiàn)[7]還將虛粒子看作流體粒子的擴(kuò)展，虛粒子有條件地參與控制方程的計算，對流體粒子的壓力、速度方向產(chǎn)生影響，但此方法對于復(fù)雜邊界效果不佳；2017年，ESLAMIAN和KHAYAT[8]用靠近邊界的2層粒子作為內(nèi)襯層代替虛粒子，避免了邊界力法施加人工力導(dǎo)致的動量不守恒問題和虛粒子法造成的不滿足連續(xù)性方程的問題，但該方法必須耦合其他壓力處理方法，否則會產(chǎn)生明顯的視覺壓縮，增加計算量。

針對上述方法計算量大、對于復(fù)雜邊界處理效果不佳的問題，本文提出了對稱區(qū)域邊界處理方法。對于靠近邊界的流體粒子，首先找到被邊界截斷的流體粒子的支持域，然后通過中心對稱，在流體粒子支持域中找到關(guān)于被截斷區(qū)域?qū)ΨQ的區(qū)域，在計算流體粒子密度時，增大該對稱區(qū)域內(nèi)的粒子權(quán)重，避免了邊界對粒子支持域截斷帶來的密度、壓力計算錯誤等問題。與傳統(tǒng)的邊界力法和虛粒子法相比，本文方法無需生成虛粒子，在保證逼真度的前提下滿足實時性要求，且隨著粒子數(shù)的增加，耗時增長也明顯比其他傳統(tǒng)方法慢，計算效率高的優(yōu)勢更明顯，更適合對復(fù)雜場景的模擬，同時避免了邊界處流體與邊界分離的現(xiàn)象。

本文的另一項工作是對傳統(tǒng)的 CSF方法進(jìn)行改進(jìn)。

1992年，BRACKBILL等[9]提出 CSF方法，通過將流體表面描述為一個具有有限厚度的過渡帶，并將表面張力轉(zhuǎn)化為體積力進(jìn)行計算，適用于任何形狀的表面，但其用一個光滑的色函數(shù)區(qū)分2種流體是沒必要的；2000年，MORRIS[10]將 CSF方法引入 SPH方法中計算表面張力，并采用光滑的色函數(shù)，但采用 SPH計算表面曲率時，散度的計算需要核函數(shù)的充分支持，而核函數(shù)在過渡帶內(nèi)是無法滿足的；2005年，TARTAKOVSKY和MEAKIN[11]針對CSF無法滿足表面張力的切向特性，提出了計算表面張力的公式，通過余弦函數(shù)定義吸引、排斥力來重現(xiàn)表面張力的影響。粒子運動具有較好的一致性，但計算得到的表面張力需要校準(zhǔn)，且在液滴的模擬過程中會出現(xiàn)收縮、膨脹交替的過程，不符合表面張力作用的實際情形。文獻(xiàn)[3]使用不連續(xù)的色函數(shù)解決了文獻(xiàn)[10]的曲率計算問題，修正了表面張力計算公式，避免了表面張力方向錯誤的問題，但計算量相對較大。2007年，HU和ADAMS[12]將CSF方法應(yīng)用于多相流的模擬；2010年，ADAMI等[13]提出一種新的散度近似方案，無須完整的核函數(shù)支持域，得到了精度較高的曲率，但其不完全遵守動量守恒；2017年，YANG等[14]將原始網(wǎng)格的每個頂點的表面張力強(qiáng)度轉(zhuǎn)移到其相鄰的粒子上，得到了較好的結(jié)果，但需多次迭代計算；2018年，KRIMI等[15]對文獻(xiàn)[13]提出的表面張力計算方法進(jìn)行改進(jìn)，提高了流體表面附近的穩(wěn)定性，考慮了表面張力的切向特性，并且可以應(yīng)用于兩相以上的流體模擬，但該方法依舊無法保證動量守恒。

針對 CSF方法將表面張力看作體積力進(jìn)行計算導(dǎo)致的曲率計算脫離表面形狀控制這一問題，本文對 CSF方法進(jìn)行了改進(jìn)，提出了基于表面粒子提取的表面張力計算方法：首先精確提取表面粒子，并通過表面粒子計算曲率，不但減小了傳統(tǒng)CSF方法計算曲率的誤差，而且減少了計算量，提高了計算速度。通過方形液滴實驗，驗證了該方法的有效性。

1 SPH方法

SPH方法是一種無網(wǎng)格插值法，通過核函數(shù)近似計算粒子的物理屬性；描述流體運動規(guī)律的N-S

方程[16]為

其中，ρ為密度；v為速度；t為時間；p為壓強(qiáng)；μ為粘性系數(shù)；g為重力加速度。式(1)等號右側(cè)依次為壓力項、粘性力項、外力項。粒子屬性通過離散求和公式[16]計算如下

其中，A(r)為位置r的目標(biāo)粒子的屬性值；N為目標(biāo)粒子鄰域內(nèi)粒子總數(shù)；mj，ρj，Aj，rj分別為目標(biāo)粒子鄰域內(nèi)粒子的質(zhì)量、密度、物理屬性值和位置；W(r-rj,h)為光滑核函數(shù)；h為光滑核半徑。

A(r)的梯度和拉普拉斯算子通過離散求和式[16]分別計算如下

其中，?為梯度算子；2?為拉普拉斯算子。

其中，r為目標(biāo)粒子與支持域內(nèi)粒子j的間距。

由式(2)得到流體粒子的密度公式[16]，即

對于弱可壓縮流體，采用泰特方程[17]計算壓強(qiáng)

其中，k為氣體常數(shù)；ρ0為流體靜止密度。

壓力和粘性力計算[16]如下

針對不同的物理量，采用 Muller的核函數(shù)，計算以下核函數(shù)：

(1) 密度核函數(shù)[16]

(2) 壓力核函數(shù)[16]

(3) 粘性力核函數(shù)[16]

其中，r為目標(biāo)粒子與支持域內(nèi)粒子j的間距。

2 對稱區(qū)域邊界處理方法

邊界力法是在邊界鋪設(shè)一層粒子并對流體粒子施加作用力，防止粒子穿透，該方法需要求解邊界力，足夠小的時間步長才能保證模擬的穩(wěn)定性。虛粒子法通過在邊界另一側(cè)鏡像生成流體粒子，并給其粒子賦密度和壓強(qiáng)來降低截斷誤差，該方法需要鋪設(shè)虛粒子，隨著粒子數(shù)的增加計算量增大；對于復(fù)雜邊界，難以確定虛粒子的位置，使得算法復(fù)雜度提高，因此傳統(tǒng)的邊界力法和虛粒子法難以滿足實時性要求。為此，本文提出一種對稱區(qū)域邊界處理方法。以下分直邊界、圓邊界、鋸齒邊界3種情形進(jìn)行討論。

(1) 直邊界。對靠近邊界的粒子i的支持域，首先找到被邊界截斷的區(qū)域2；再確定支持域中與區(qū)域2對稱的區(qū)域1；然后計算粒子的密度：計算區(qū)域3內(nèi)的粒子的密度時，按正常權(quán)重計算；計算區(qū)域1時，將其權(quán)重增大到原權(quán)重的2倍，如圖1所示。

圖1 直邊界的情形

圖1中，粗直線表示固壁邊界，虛線與粗線關(guān)于粒子i對稱，大圓表示粒子i的支持域，小圓圈表示支持域內(nèi)的其他流體粒子，區(qū)域 2表示支持域被邊界截斷的部分，區(qū)域1與區(qū)域2關(guān)于粒子i對稱。

對于直邊界，當(dāng)粒子碰到邊界時，若粒子速度為v，將v在邊界方向和與之垂直的方向進(jìn)行分解，于是反彈后的速度為

(2) 圓邊界。設(shè)圓形邊界所在的圓為以點O為圓心，R為半徑的圓O；對靠近邊界的粒子i，圓O關(guān)于i對稱的圓為圓O′。對粒子i的支持域，首先找到被邊界截斷的區(qū)域5；再確定支持域中與區(qū)域5對稱的區(qū)域4；然后計算粒子的密度：計算區(qū)域6內(nèi)的粒子的密度時，按正常權(quán)重計算；計算區(qū)域4內(nèi)的粒子的密度時，將其權(quán)重增大到原權(quán)重的2倍。如圖2所示。

圖2中，粗線圓為固壁邊界所在的圓，虛線圓與粗線圓關(guān)于粒子i對稱，中間的圓為粒子i的支持域，小圓圈表示支持域內(nèi)的其他流體粒子，區(qū)域5表示支持域內(nèi)被邊界截斷的部分，區(qū)域4與區(qū)域5關(guān)于粒子i對稱。

圖2 圓形邊界的情形

當(dāng)粒子到達(dá)邊界時會發(fā)生反彈，直邊界反彈處理比較簡單。對于圓形邊界，v在切線和垂直切線的方向進(jìn)行分解，得到切向和法向速度，即

其中，α為速度與法向的夾角；vt為切向速度；vn為法向速度，如圖3所示。

圖3 圓形邊界速度處理

轉(zhuǎn)換到直角坐標(biāo)系后，反彈后的粒子速度為

其中，θ為切線方向與水平方向的夾角。

(3) 鋸齒邊界。對靠近邊界的粒子i的支持域，首先找到被邊界截斷的區(qū)域8；再確定支持域中與區(qū)域8對稱的區(qū)域7；然后計算粒子的密度：計算區(qū)域9內(nèi)的粒子的密度時，按正常權(quán)重計算；計算區(qū)域7內(nèi)的粒子的密度時，將其權(quán)重增大到原權(quán)重的2倍，如圖4所示。

圖4 鋸齒邊界的情形

圖4中，粗直線表示固壁邊界，虛線與粗線關(guān)于粒子i對稱，大圓表示粒子i的支持域，小圓圈表示支持域內(nèi)的其他流體粒子，區(qū)域 8表示支持域被邊界截斷的部分，區(qū)域7與區(qū)域8關(guān)于粒子i對稱。

對于鋸齒邊界，其速度反彈的處理方法與圓形邊界速度反彈處理方法相同。模擬仿真效果見4.1節(jié)。

3 基于表面粒子提取的表面張力計算

CSF方法通過將流體表面描述為一個具有有限厚度的過渡帶，將表面張力轉(zhuǎn)化為體積力進(jìn)行計算。2種流體分別被賦予了不同的顏色值，如研究液滴的表面張力時，將液體粒子的色值設(shè)為1，將空氣粒子的色值設(shè)為0，根據(jù)SPH思想得到一個色函數(shù)，通過計算得到每個粒子相應(yīng)的色函數(shù)值，當(dāng)其大于閾值時計算表面張力。

在運用 CSF的過程中，雖然流體粒子整體分布是均勻的，但就局部而言，流體粒子分布并不均勻，致使表面曲率計算錯誤，從而導(dǎo)致表面張力的計算錯誤。

因為表面張力是一種表面力，曲率的計算只與表面的形狀有關(guān)。為此，本文對 CSF方法進(jìn)行了改進(jìn)，提出了一種基于表面粒子提取的表面張力計算方法。

提取表面粒子，首先要確定該粒子是表面粒子，下面給出判斷表面粒子的方法。以粒子i為向量起點做2個相互垂直的向量a和b，將粒子i的支持域內(nèi)的粒子j以粒子i為起點做向量xij，并與向量a和b的數(shù)量積的正負(fù)來判斷粒子i是否為表面粒子。若粒子i支持域內(nèi)所有粒子j滿足a·xij＜0或b·xij＜0，則說明從向量a到向量b的圓心角為90°的扇形區(qū)域內(nèi)沒有粒子，粒子i為表面粒子；否則，將向量a和b同時旋轉(zhuǎn)，當(dāng)旋轉(zhuǎn)到某個角度(設(shè)這個角度為α(0＜α＜2π))，粒子i支持域內(nèi)所有粒子j都滿足a·xij＜0 或b·xij＜0，則停止旋轉(zhuǎn)，此時粒子i為表面粒子；若向量a和b同時旋轉(zhuǎn)，直到轉(zhuǎn)回初始位置，粒子i支持域內(nèi)所有粒子j均不滿足a·xij＜0 且b·xij＜0，則粒子i不是表面粒子。如圖5所示。

圖5中，黑色實心小圓分別為粒子i的2種情況，虛線圓為粒子i的支持域。①粒子i位于圖中上方，此時i的支持域中從向量a到向量b的扇形區(qū)域內(nèi)無其他流體粒子，此時i為表面粒子；②粒子i位于圖中下方，此時i的支持域中，無論向量a和b如何旋轉(zhuǎn)，向量a與向量b之間的扇形區(qū)域內(nèi)均有其他流體粒子，此時i不是表面粒子。

圖5 表面粒子提取

下一步根據(jù)提取出的表面粒子計算曲率。將2種流體粒子的色值分別設(shè)置為1和0，其色函數(shù)[16]為

在CSF模型中，表面張力計算模型[16]為

其中，σ為表面張力系數(shù)；單位法向n= ?cs；曲率κ=-(? ·n) = -?2c。具體計算如式(3)～(6)。

計算曲率時，目標(biāo)粒子i鄰域內(nèi)的粒子是通過上一步表面提取得到的邊界粒子。模擬仿真效果見4.2節(jié)。

4 仿真實驗與分析

實驗在 Windows 7平臺上進(jìn)行，開發(fā)環(huán)境為visual studio 2013，實驗配置為Intel(R) Core(TM)i5-2400，3.10 GHz CPU，4 G 內(nèi)存，NVIDIA Geforce GT 620 顯卡。

4.1 二維潰壩流固邊界處理仿真實驗

為了驗證本文提出的對稱區(qū)域邊界處理方法的有效性，將邊界力法、動態(tài)虛粒子法、靜態(tài)虛粒子法及本文方法分別應(yīng)用于二維潰壩場景的直邊界處理(圖 6～8，表 1)，比較模擬效果，粒子數(shù)分別為1 156，2 809，5 329。再將上述3中傳統(tǒng)方法及本文方法分別應(yīng)用于二維潰壩場景的圓邊界處理(圖 9～11，表 2)及鋸齒邊界處理(圖 12～14，表 3)，比較模擬效果，粒子數(shù)分別為1 550，3 010，5 712。實驗時間步長取0.005 s，初始密度取1 000 kg/m3，粒子支持域半徑取0.04 m。

圖6 直邊界情形下，4種方法在第230幀時的模擬效果(粒子數(shù)為1 156)

圖7 直邊界情形下，4種方法在第230幀時的模擬效果(粒子數(shù)為2 809)

圖8 直邊界情形下，4種方法在第240幀時的模擬效果(粒子數(shù)為5 329)

表1 直邊界情形下，4種方法在不同粒子數(shù)及幀數(shù)時運行耗時表(s)

圖9 圓邊界情形下，4種方法在第170幀時的模擬效果(粒子數(shù)為1 550)

圖10 圓邊界情形下，4種方法在第160幀時的模擬效果(粒子數(shù)為3 010)

圖11 圓邊界情形下，4種方法在第120幀時的模擬效果(粒子數(shù)為5 712)

表2 圓邊界情形下，4種方法在不同粒子數(shù)及幀數(shù)時運行耗時表(s)

圖12 鋸齒邊界情形下，4種方法在第210幀時的模擬效果(粒子數(shù)為1 550)

圖13 鋸齒邊界情形下，4種方法在第196幀時的模擬效果(粒子數(shù)為3 010)

圖14 鋸齒邊界情形下，4種方法在第196幀時的模擬效果(粒子數(shù)為5 712)

表3 鋸齒邊界情形下，4種方法在不同粒子數(shù)及幀數(shù)運行耗時表(s)

由圖6～14可以看出，靜態(tài)和動態(tài)虛粒子法會出現(xiàn)由于虛粒子產(chǎn)生的壓強(qiáng)過大導(dǎo)致的流體與邊界分離等不穩(wěn)定現(xiàn)象；邊界力法由于時間步長的限制導(dǎo)致模擬速度緩慢；本文方法不僅不會出現(xiàn)上述問題，而且模擬效果穩(wěn)定、符合實際。由表1～3可以看出，本文方法在很大程度上減小了計算量，且隨著粒子數(shù)的增加，本文方法實時性的優(yōu)勢更加明顯，這為復(fù)雜邊界的處理提供了基礎(chǔ)。

4.2 形成圓形液滴的仿真實驗

在處理表面張力時，本文可精確提取表面粒子，通過表面粒子計算曲率，減小了傳統(tǒng) CSF方法計算曲率的誤差，提高了模擬效率。

在進(jìn)行圓形液滴形成的模擬實驗中，實驗粒子數(shù)為1 024，粒子支持域半徑為0.04 m，時間步長為0.02 s。圖15為表面粒子提取效果，其中紅色粒子為內(nèi)部粒子，藍(lán)色粒子為表面粒子。圖16為文獻(xiàn)[11](利用余弦函數(shù)計算表面張力)、傳統(tǒng)CSF方法及本文方法模擬圓形液滴形成過程的效果對比。

圖15 表面提取效果(藍(lán)色粒子為表面粒子，紅色粒子為內(nèi)部粒子)

圖16 3種方法模擬液滴形成不同時刻的模擬效果

圖16中3種方法運行到200幀時，耗時分別為4.5 s，4.2 s，3.5 s，可以看出，與文獻(xiàn)[11]的方法相比，本文方法在幀數(shù)為200時已經(jīng)得到了很好的效果，而文獻(xiàn)[11]方法的液滴還處于變化狀態(tài)。此外，本文方法的液滴在整個變化過程中一直處于擴(kuò)張狀態(tài)，而文獻(xiàn)[11]方法會出現(xiàn)收縮、膨脹交替的過程，不符合表面張力作用的實際情形；與傳統(tǒng)CSF方法相比，本文方法形成的圓形液滴更光滑?？傊疚姆椒M圓形液滴形成時，不但效果好，而且耗時短。

5 結(jié) 論

針對傳統(tǒng)的固壁邊界處理方法中的邊界力法、虛粒子法計算量大、對于復(fù)雜邊界處理難以達(dá)到實時性的問題，本文提出一種對稱區(qū)域邊界處理方法。對于靠近邊界的流體粒子，首先找到被邊界截斷的流體粒子的支持域，然后通過中心對稱，在流體粒子支持域中找到關(guān)于被截斷區(qū)域?qū)ΨQ的區(qū)域，在計算流體粒子密度時，對這個對稱區(qū)域內(nèi)的粒子增大權(quán)重，避免了邊界對粒子支持域截斷帶來的密度、壓力計算錯誤等問題。與經(jīng)典傳統(tǒng)方法相比，該方法無需生成虛粒子，減少了計算量，在保證逼真度的前提下達(dá)到實時性，且隨著粒子數(shù)的增加，該方法的耗時增長也明顯比其他傳統(tǒng)方法慢，計算效率高的優(yōu)勢更明顯，為復(fù)雜場景的模擬奠定了良好的基礎(chǔ)。此外，針對傳統(tǒng) CSF方法計算表面張力時曲率計算存在的誤差，本文進(jìn)行了修正，只通過邊界粒子計算邊界粒子的曲率，減小了計算量，提高了模擬效率，且模擬效果更好。此外，本文的對稱區(qū)域邊界處理方法可以推廣到任意復(fù)雜邊界，但邊界越復(fù)雜，本文方法的計算也越復(fù)雜；但比傳統(tǒng)的生成虛粒子的計算簡單快速。

本文在基于表面粒子提取的表面張力計算中，用到余弦函數(shù)，增大了計算量，計算效率有待進(jìn)一步提高。