韋 建, 朱曉臨, 張義群

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

關(guān)于多相流動,文獻(xiàn)[1]認(rèn)為根據(jù)固體、液體、氣體中的任何兩相組合或三相組合,可以分為氣液兩相流、液液兩相流、氣固兩相流、固液兩相流、氣液固三相流等。近年來,氣液兩相流動數(shù)值模擬技術(shù)是計(jì)算流體動力學(xué)(computational fluid dynamics,CFD)領(lǐng)域研究的熱點(diǎn)和前沿課題之一。在圖形圖像模擬方面,氣液兩相流的交互模擬可以增強(qiáng)真實(shí)感,可廣泛用于游戲、動畫及影視制作等;在實(shí)際工程中,如針對水利工程中的問題,研究氣液交互也有重要作用,只有很好地把握波浪與結(jié)構(gòu)的相互作用關(guān)系,研究波的傳播、高度差引起的水跌、水底氣泡浮升時空氣體積因壓力所引起的變化、正在下降的水滴運(yùn)動,或者是搖晃密度不同流體所引起的界面流動等,才能更好地增加工程的穩(wěn)定性和可靠性。

目前,常用光滑粒子流體動力學(xué) (smoothed particle hydrodynamics,SPH)方法進(jìn)行液體模擬。 SPH積分時,鄰近邊界粒子的支持域會被問題域截?cái)?這種單邊影響會導(dǎo)致求解結(jié)果出現(xiàn)較大誤差,因此必須通過一些特殊方式對邊界進(jìn)行處理,且處理應(yīng)該盡可能小地對整個流場產(chǎn)生擾動。因此,在SPH方法中,邊界處理一直是研究的重點(diǎn)和難點(diǎn)。在應(yīng)用SPH方法模擬流體時,在氣液邊界處的流體粒子附近,流體粒子數(shù)目急劇減少,因此存在邊界截?cái)嗾`差問題,從而影響流體模擬的精度以及模擬的穩(wěn)定性。常見的處理氣液邊界的方法主要有2種,即虛粒子法和空氣粒子法。

虛粒子法是在靠近氣液邊界的流體粒子上方生成若干層虛粒子,用來填補(bǔ)流體粒子被邊界截?cái)嗟闹С钟?。文獻(xiàn)[2]應(yīng)用SPH方法模擬流體時,首次對氣液邊界粒子密度的計(jì)算方法進(jìn)行了改進(jìn),雖然模擬效果提升得不明顯,但是這為后續(xù)研究氣液邊界奠定了基礎(chǔ);文獻(xiàn)[3]在氣液邊界粒子上方鏡像生成多層虛粒子,參與表面張力的計(jì)算并作用于小型氣泡,但是該方法產(chǎn)生的虛粒子的密度值差異較大,且數(shù)值變換劇烈;文獻(xiàn)[4]提出泰特方程來降低不可壓縮性,并且利用高阻尼方程初始化粒子,使靠近自由表面的粒子能夠達(dá)到相對穩(wěn)定的密度值,但是該方法只適用于初始粒子形狀規(guī)整、分布均勻的情況;文獻(xiàn)[5]提出ghost方法,利用Poisson-disk算法生成的虛粒子具有穩(wěn)定的狀態(tài)分布,但是該算法的計(jì)算量較大;文獻(xiàn)[6]提出一種氣體壓力法來生成虛粒子,與ghost方法相比,該方法計(jì)算量減少、耗時短,但是氣液表面的精度和細(xì)節(jié)有所降低;文獻(xiàn)[7]先在流體模擬域內(nèi)初始均勻分布ghost粒子,且ghost粒子位置保持不變,再利用一個判定機(jī)制選出自由表面處的ghost氣體粒子,該方法一定程度上解決了動態(tài)生成ghost粒子困難的問題,但是在容器較大、流體運(yùn)動比較平緩的情況下,預(yù)先固定鏡像粒子將造成一定的浪費(fèi);文獻(xiàn)[8]提出一種空氣粒子均勻化分布的方法,將空氣粒子和流體粒子看成同一種粒子進(jìn)行處理,在所有的粒子之間另外施加一個余弦函數(shù)的力,在一定時間步空氣粒子可以達(dá)到足夠的均勻狀態(tài),但加入余弦函數(shù)意味著增加了很多的計(jì)算量,降低了模擬效率。

空氣粒子法是把空氣也看作一種流體,則氣液邊界也是一種液液邊界,只是兩相粒子的密度比為1∶1 000。因此,如何處理該密度比的多相流交互成為處理氣液邊界的關(guān)鍵。文獻(xiàn)[9]使用密度重整化方法,對空氣粒子使用大的人工表面張力和較高的波速,使速度場平滑,但是需要非常小的時間步長來維持模擬系統(tǒng)的穩(wěn)定,這導(dǎo)致模擬時間較長;文獻(xiàn)[10]提出一種改進(jìn)的多相SPH方法,并將其應(yīng)用于流體工程以估算多流體的密度,該方法能夠處理跨相交界面的密度不連續(xù)性,允許不連續(xù)的黏度,但是確保跨相交界面上的速度和應(yīng)力的連續(xù)性,該方法的缺點(diǎn)是增加了計(jì)算量,且添加表面張力后的效果不是很理想;文獻(xiàn)[11]改進(jìn)SPH方法中的密度計(jì)算公式,模擬了密度比為1∶1 000的空氣和水交互場景,但是該方法對多相流之間的波速比有較大的限制;文獻(xiàn)[12]使用(密度加權(quán))平均壓力來處理氣液邊界上的壓力梯度,但是邊角處容易出現(xiàn)粒子不穩(wěn)定的情況;文獻(xiàn)[13]提出一種排斥力來穩(wěn)定兩種流體之間的界面,將其應(yīng)用于密度比為1∶1 000的多相流模擬中,但是,該方法需要低密度粒子的數(shù)目是高密度粒子數(shù)的5～7倍,否則模擬效果較差,局限性較大;文獻(xiàn)[14]利用壓力泊松方程建模,使氣液邊界附近流體粒子的密度保持在一定范圍內(nèi),但是該方法使體積守恒的誤差比較大,且需要增加補(bǔ)償項(xiàng)來處理交界面的密度不連續(xù)問題;文獻(xiàn)[15-16]提出一種不可壓縮SPH和弱可壓縮SPH的耦合方法,對速度項(xiàng)應(yīng)用不可壓縮法,可壓縮法將壓力項(xiàng)作為界面處另一相的邊界條件,又在可壓縮相中加入人工黏度項(xiàng)用于界面的數(shù)值穩(wěn)定,該方法可用于模擬密度比為1∶1 000和有較大波速比的多相流,但是,該方法實(shí)施起來比較麻煩,模型中不確定參數(shù)較多,模擬難度較大;文獻(xiàn)[17]基于密度校正提出了一種新型WCSPH流體模型,該方法可以有效處理密度比為1∶1 000的氣液交互面密度連續(xù)性問題,能夠模擬與實(shí)際相符的波速比(4∶1)效果,但是該方法在氣體粒子與液體粒子的數(shù)目比比較大的情況下,空氣粒子會穿透水面,影響模擬效果。

與虛粒子法相比,空氣粒子法更符合現(xiàn)實(shí)情況(空氣也是粒子組成的),且在模擬多相流時,相數(shù)越多,空氣粒子法優(yōu)勢越大,即不用為每相流的氣液邊界動態(tài)生成虛粒子,減少了計(jì)算量,更容易使多相流的氣液邊界保持穩(wěn)定。

本文在處理氣液邊界時提取出界面粒子(氣液交界處“最薄”的2層粒子),與這些界面粒子產(chǎn)生交互的粒子標(biāo)記為“外交粒子”,這些“外交粒子”可用來加快界面粒子的鄰近粒子搜索,也可用來處理氣液邊界。通過實(shí)驗(yàn)?zāi)M時間對比發(fā)現(xiàn),相比于文獻(xiàn)[18],本文方法可以有效減少模擬時間,提高模擬效率。本文提出一種改進(jìn)的空氣粒子法壓力求解公式,將其應(yīng)用于上述的“外交粒子”來達(dá)到處理氣液邊界的目的。通過實(shí)驗(yàn)?zāi)M效果對比發(fā)現(xiàn),相比于文獻(xiàn)[17]方法,本文方法使得氣液交界面更光滑、清晰和無穿透現(xiàn)象。

1 SPH方法

SPH方法[19]中,在支持域Ω上,關(guān)于位置矢量r的屬性方程f定義為:

(1)

其中,r′為鄰近粒子位置矢量;δ(r-r′)為狄拉克函數(shù),即

(2)

在數(shù)值計(jì)算中用光滑函數(shù)W(r-r′,h)來取代狄拉克函數(shù),則f(r)為:

(3)

其中,W為平滑核函數(shù);h為核函數(shù)的有效支持長度。

采用粒子近似方法來對核函數(shù)插值離散,粒子i的場函數(shù)值可以通過核函數(shù)對該粒子支持域內(nèi)所有粒子函數(shù)值加權(quán)平均得到,即

(4)

其中,N為流體粒子i鄰近域內(nèi)的粒子總數(shù);j為粒子i的鄰近粒子;mj為粒子j的質(zhì)量;ρj為粒子j的密度。

對f(r)求梯度與求拉普拉斯運(yùn)算即是對核函數(shù)進(jìn)行梯度運(yùn)算與拉普拉斯運(yùn)算,即

(5)

(6)

2 界面粒子的鄰近粒子搜索

在用SPH方法模擬氣液多相流運(yùn)動時,一個必要的步驟是提取界面粒子,主要用來計(jì)算曲率(用于表面張力計(jì)算)和表面追蹤(用于表面重建);另一個必不可少的步驟是搜索每個粒子的鄰近粒子,這也是模擬中最耗時的步驟。在氣液多相流中,當(dāng)氣液交互較為劇烈時會產(chǎn)生大量的界面粒子,如圖1所示。圖1中,綠色和紅色粒子分別是A和B中的界面粒子。本文為界面粒子增加一個新功能,即幫助氣液多相流中所有相流的界面粒子搜索它們的鄰近粒子。

A—?dú)庀嗔?B—液相流圖1 氣液多相流中的界面粒子示意圖

在提取界面粒子時,將與界面粒子交互的粒子提取并存儲,這些粒子稱為“外交粒子”(與本身所在相流的界面粒子或者其他相流的界面粒子產(chǎn)生交互的粒子)。在計(jì)算界面粒子屬性需要搜索其鄰近粒子時,只需要在“外交粒子”集合中遍歷即可。

界面粒子一定是“外交粒子”,“外交粒子”不一定是界面粒子,如圖2所示。圖2中,綠色和藍(lán)色粒子都是“外交粒子”。

圖2 氣液多相流中的“外交粒子”示意圖

本文氣液多相流模擬的主要步驟如下:

(1) 初始化粒子位置,賦予粒子初始屬性。

(2) 提取界面粒子[18],在提取界面粒子過程中即可提取出“外交粒子”集合。

(3) 遍歷粒子,搜索每個粒子的鄰近粒子,計(jì)算目標(biāo)粒子的密度、壓力、速度及加速度等屬性。當(dāng)目標(biāo)粒子為界面粒子時,鄰近粒子只需要在“外交粒子”集合中搜索。

(4) 更新粒子位置,返回到步驟(2)。

3 壓力修正求解公式

利用文獻(xiàn)[17]方法處理密度比為1∶1 000的氣液交互面密度連續(xù)性問題時,在氣體粒子與液體粒子的數(shù)目比比較大的情況下,空氣粒子會穿透浪花,使得水波或浪花不正常斷裂,呈現(xiàn)出類似一串珍珠的樣子(應(yīng)當(dāng)只在浪的前端出現(xiàn)破碎、細(xì)小及飛濺的浪花),氣液交界面的局部存在粒子穿透與“毛糙”界面。出現(xiàn)上述情況的原因是文獻(xiàn)[17]模型中壓力求解公式在細(xì)節(jié)處的不適用,即在求解氣體粒子屬性時,特別是在氣液邊界處,容易產(chǎn)生數(shù)值不穩(wěn)定,產(chǎn)生較大速度和壓力,穿透了聚集在一起的液體粒子團(tuán)。

本文對“外交粒子”應(yīng)用一個新壓力求解公式(其他粒子仍使用文獻(xiàn)[17]中的壓力求解公式),以防止氣液交界面發(fā)生粒子穿透和氣液粒子混合影響交界面的清晰、光滑度。

“外交粒子”i的壓力求解公式為:

(7)

(8)

根據(jù)(7)式,粒子i受到粒子j的作用力為:

(9)

粒子j受到粒子i作用力為:

(10)

(7)式對支持域中有不同相粒子的“外交粒子”增加了壓力修正項(xiàng),使液體表面更光滑,具體分析如下。

(11)

則對這部分粒子不進(jìn)行壓力修正,(7)式化為:

(12)

(13)

(13)式可防止不同相粒子在界面出現(xiàn)穿透現(xiàn)象,并保持界面清晰和光滑,這與實(shí)際情況保持一致(與物理上的相斥流體界面現(xiàn)象相對應(yīng),即當(dāng)流體與另外一種與其不相容的流體相互接觸時,交界面上的力會保持界面的光滑性與清晰性)。

4 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)采用文獻(xiàn)[17-18]和本文的方法模擬多相流潰壩場景。實(shí)驗(yàn)在Window 7系統(tǒng)上,平臺為Visual Studio 2013。

4.1 界面粒子的鄰近粒子搜索

采用文獻(xiàn)[18]方法和本文通過提取“外交粒子”方法加快界面粒子中鄰近粒子搜索的方法模擬多相流(4.2中氣液交互場景),粒子總數(shù)為5 000、10 000、20 000、40 000時,2種方法的關(guān)鍵幀用時如圖3所示。

圖3 4種粒子總數(shù)下關(guān)鍵幀用時對比

4組實(shí)驗(yàn)?zāi)M總時長對比見表1所列。模擬幀數(shù)均為2 342幀,實(shí)驗(yàn)效率提高程度η計(jì)算公式為:

(14)

表1 4種粒子總數(shù)下模擬時長與實(shí)驗(yàn)效率對比

其中,t1為文獻(xiàn)[18]方法用時;t2為本文方法用時。

從4組實(shí)驗(yàn)結(jié)果對比可以看出,本文方法能夠提高模擬效率,增強(qiáng)模擬實(shí)時性。本文方法主要是加快對表面粒子中鄰近粒子的搜索,當(dāng)模擬場景中的粒子數(shù)較多且多相流之間交互劇烈(交互劇烈會產(chǎn)生大量的界面粒子)時,與文獻(xiàn)[18]方法相比,本文方法不僅提升了模擬效率,而且穩(wěn)定性更好,因此本文方法能更高效地處理復(fù)雜場景。

4.2 氣液邊界壓力修正

應(yīng)用文獻(xiàn)[17]方法和本文方法模擬氣液交互的場景。

氣液交界面模擬效果對比如圖4所示。圖4中,黃色是空氣粒子,紅色是流體粒子。

從圖4可以看出,在氣液交界面,文獻(xiàn)[17]方法存在粒子穿透與不光滑界面(用黑色圈標(biāo)記),本文方法界面更光滑、更清晰,且沒有穿透現(xiàn)象。

圖4 氣液交界面模擬效果對比

風(fēng)吹水面二維模擬效果對比如圖5所示。圖5a、圖5b中,從上至下是使用不同參數(shù)模擬不同“風(fēng)力”的模擬效果。

圖5 風(fēng)吹水面二維模擬效果對比

從圖5可以看出,采用文獻(xiàn)[17]方法出現(xiàn)水波不正常斷裂,呈現(xiàn)類似一串珍珠的樣子,與實(shí)際不符,模擬效果不佳;采用本文方法模擬的水波連續(xù)且光滑,與現(xiàn)實(shí)相符,模擬效果有明顯提升。

應(yīng)用本文方法進(jìn)行三維渲染實(shí)驗(yàn)的效果如圖6所示。

圖6 本文方法三維渲染實(shí)驗(yàn)效果

5 結(jié) 論

本文基于SPH方法提出一種多相流氣液邊界處理方法,解決了以下2個問題:① 針對氣液邊界交互粒子多,傳統(tǒng)SPH方法搜索鄰近粒子比較耗費(fèi)時間的問題,本文通過提取“外交粒子”來減少搜索界面粒子鄰近粒子的時間,該方法對于粒子總數(shù)多、交互劇烈的多相流模擬,效率有顯著提高;② 針對氣液交互面出現(xiàn)粒子穿透,交界面不光滑、“毛糙”等問題,本文對“外交粒子”應(yīng)用一種新的壓力計(jì)算公式達(dá)到處理氣液邊界的目的,使得交界面清晰、穩(wěn)定且沒有粒子穿透現(xiàn)象。多個仿真實(shí)驗(yàn)結(jié)果驗(yàn)證了本文方法的有效性。