多相流氣液邊界處理

韋 建, 朱曉臨, 張義群

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

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

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

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

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

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

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

1 SPH方法

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

(1)

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

(2)

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

(3)

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

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

(4)

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

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

(5)

(6)

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

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

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

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

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

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

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

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

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

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

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

3 壓力修正求解公式

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

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

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

(7)

(8)

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

(9)

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

(10)

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

(11)

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

(12)

(13)

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

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

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

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

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

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

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

(14)

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

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

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

4.2 氣液邊界壓力修正

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

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

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

圖4 氣液交界面模擬效果對(duì)比

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

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

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