劉漢濤 劉謀斌常建忠? 蘇鐵熊

1 引言

微觀或介觀尺度下多相流體流動包含復雜的流體動力學特性,對于現(xiàn)代微納米工程、生物工程、水資源及環(huán)境工程具有極為重要的作用.然而目前由于實驗研究和數(shù)值模擬的局限性,對其復雜的流動特性還未能很好地理解和掌握[1-3].

耗散粒子動力學最初是由Hoogerbrugge和Koelman[4,5]在研究懸浮時引入的,并結合了分子動力學、格子氣自動機和格子-波爾茲曼方法的優(yōu)點.Espanol和Warren[6],Marsh等[7]的工作奠定了耗散粒子動力學(DPD)方法在統(tǒng)計力學方面的基礎,DPD方法已成功應用于高分子溶液、膠體及多相流的模擬.與光滑粒子動力學(SPH)方法相比,DPD方法主要區(qū)別在于粒子間的相互作用力除耗散力外還包含了隨機力和保守力,隨機力引起粒子間的隨機振動,增加系統(tǒng)的動能,提高系統(tǒng)的溫度,這種振動只有在微觀尺度時才凸顯出來,這使DPD方法本質上成為一種介觀尺度上的模擬技術.國內采用DPD方法對簡單單相流體的流動及液滴的破碎和碰撞進行了模擬[8],國外DPD應用于多相系統(tǒng)的研究也只局限于相分離、不同相在界面處排斥的界面過程.Visser等[9]通過控制相互作用的粒子間耗散力近似不同粒子間黏度,對兩種不同黏度的流體系統(tǒng)進行了研究.模擬中采用將不同相的流體通過鏡面反射在分界面上嚴格分開,但實際上不同相的流體在界面上可以相互滲透,粒子可以通過界面進行擴散.

DPD方法一般采用傳統(tǒng)平方形式的權函數(shù)(1-r/rc)2,或為提高系統(tǒng)施密特數(shù),耗散力權函數(shù)和隨機力權函數(shù)的形式可以采用如下形式：

式中,rc為截距,粒子-粒子間的兩兩相互作用力局限在有限的截距內.耗散力權函數(shù)和隨機力權函數(shù)采用上述形式時,粒子間表現(xiàn)為純排斥性的作用力,導致DPD粒子在計算區(qū)域內因相互排斥而分散開來,從而占據(jù)整個計算區(qū)域.因此,基于這種純排斥性權函數(shù)的DPD方法實際上模擬的是氣體.雖然有的研究者利用傳統(tǒng)的DPD方法做了一些研究,這些研究都是基于模擬多組分系統(tǒng)性液-液和液-固交接面,而并非單一組分下帶有液-氣共存流體動力學特性.因此傳統(tǒng)的帶純排斥性權函數(shù)的DPD方法不能模擬帶有液-氣共存的多相流體或帶有自由表面的流動.

2 DPD方法及參數(shù)

文中DPD方法的基本思想、控制方程及數(shù)值方法可參見文獻[8].

文獻[10—12]采用光滑粒子動力學的分段光滑三次函數(shù),通過疊加不同作用強度和截距,構造了一種部分區(qū)域為正、部分區(qū)域為負的函數(shù),作為DPD模型中粒子間相互作用的保守力勢函數(shù),本文采用Liu[13]提出的四次方光滑函數(shù)：

式中,rc為截距,根據(jù)正則化或歸一化條件要求在一維、二維和三維空間中,αd分別取為 1/h,15/7πh,315/208πh3.

構造勢函數(shù)：

式中,W(r,rc1),W(r,rc2)為(2)式所表示的未正則化的四次函數(shù),rc1,rc2分別是權函數(shù)的截距,A,B為權函數(shù)常系數(shù),a為粒子間相互作用力系數(shù).圖1為取不同截距、不同系數(shù)時的勢函數(shù).

圖1 a=1,A=2,rc1=0.8,B=-1,rc2=1.0時四次方函數(shù)及新構造的勢函數(shù)U(r)

與三次函數(shù)相比,四次函數(shù)具有更高的精確性、穩(wěn)定性,并且是一條連續(xù)函數(shù),而不是分斷連續(xù)的函數(shù).

DPD模型中保守力權函數(shù)是勢函數(shù)的導數(shù)：

如圖2所示,可以看出,DPD粒子近時表現(xiàn)為排斥形式,距離趨遠時表現(xiàn)為吸引形式.這種遠程吸引近距排斥的勢函數(shù),使得模擬液氣共存的多相流動過程成為可能.本文利用新構造的保守力勢函數(shù),對流體流過十字型通道進行了計算.計算過程中采用的DPD參數(shù)如表1所示.

表1 模擬過程中采用的DPD參數(shù)

圖2 由U(r)函數(shù)構造的DPD粒子間相互作用的保守力

3 算例及分析

計算時首先將粒子注射到40×40×3的計算區(qū)域,粒子密度ρ=4,達到平衡后,只保留計算邊界一個DPD單位的粒子,凍結作為十字型固壁邊界,因而固壁粒子是隨機布置,與采用規(guī)則分布的固壁粒子相比,這種新的壁面能有效模擬復雜固壁,并準確實施壁面邊界條件[14].固壁粒子數(shù)為1655,模擬中DPD粒子從十字型通道左部隨機注射進去,注射速率是每100時間步100個粒子,流體與固體間保守力系為aw,aw/af=2,粒子在向右的外力g=0.01驅動下流動.

圖3為粒子在不同時間步的分布情況,若采用傳統(tǒng)的只包含排斥力的權函數(shù),DPD粒子會充滿整個計算區(qū)域,只能模擬類似氣體的流體,采用包含遠程吸引近距排斥的保守力權函數(shù)后,可以模擬液氣共存的兩相流動.可以看出,流體在流動過程中,流體與壁面形成一個較小的接觸角,在流過十字型通道分叉處后,以相同的流動速度和流動型式分別流入兩個通道,流動過程中接觸角度保持不變.在流體的表面張力、向右的流動重力以及流體與固壁相互作用下,得到了其流動過程及流動模式.圖4為采用Volume-of-f l uid(VOF)方法在不同時刻模擬的兩相流組分分布情況,模擬流體采用常溫下水的性質,同樣在向右的重力作用下流動,預先指定接觸角25°.無量綱的DPD方法在三個不同時刻中,采用的時間步間的比例為1：1.27：1.62,與VOF方法的真實時間步間的比例一致.與基于網(wǎng)格的方法可以預先指定接觸角不同,DPD方法中的接觸角由壁面與流體粒子間的相對位置計算得到.可以看出,兩幅圖中DPD方法與VOF方法得到的計算結果很符合.

每100時間步注入粒子數(shù)100個,而粒子間相互作用力比值為aw/af=5時,DPD模擬結果如圖5所示.流體具有很強的潤濕能力,此時流體與壁面間的接觸角變得很小,在流過十字型通道分叉口前,接觸角基本保持不變,在分叉口處,流體被黏附在兩側壁面上,分別經(jīng)兩個通道流出,流體與壁面的接觸角保持不變.

當增大流體與固體壁面間的作用力,降低注入粒子數(shù)時,采用每100時間步注入流體粒子20個,aw/af=10,x方向重力加速度g=0.01.如圖6,由于流體流速較低,流體與壁面間相互作用力很強,流體在整個十字型通道內都會形成薄膜流動.

圖3 DPD模擬中不同時刻粒子分布,每100時間步注入100個粒子,aw=2.0af,g=0.01 (a)3300;(b)4200;(c)5400

圖4 VOF方法模擬不同時刻兩相流動分布 (a)0.37;(b)0.47;(c)0.60

圖5 DPD模擬中不同時刻粒子分布,每100時間步注入100個粒子,aw=5.0af,g=0.01

圖6 DPD模擬中不同時刻粒子分布,每100時間步注入20個粒子,aw=10.0af,g=0.01

4 結論

1)由于四次函數(shù)具有較高的精確性、穩(wěn)定性,是一條連續(xù)函數(shù),采用四次方光滑函數(shù)構造了包含遠程吸引近距排斥的保守力勢函數(shù),并對十字型通道內的流動過程進行了計算,得到了其流動過程和流動模擬,并與基于網(wǎng)格的VOF計算結果進行了對比,計算結果一致;

2)采用不同的aw/af,粒子注入速率模擬了不同的多相系統(tǒng)界面和接觸線動力學特征,得到了其流動過程及模式.這種新的勢函數(shù)也可采用不同的排斥項系數(shù)A和不同的吸引項系數(shù)B以獲得不同的流體屬性,如流體表面張力、輸運性質等.

[1]Tindall J A,Kunkel J R 1999 Unsaturated Zone Hydrology for Scientists and Engineers(New Jersey：Prentice Hall)p273

[2]Chang C L 2006 J.Phys.D：Appl.Phys.39 1132

[3]Muhammad A,Kamal M A 2009 Appl.Math.Model.33 1933

[4]Hoogerbrugge P J,Koelman J M V A 1992 Europhys.Lett.19 155

[5]Koelman J M V A,Hoogerbrugge P J 1993 Europhys.Lett.21 363

[6]Espanol P,Warren P 1995 Europhys.Lett.30 191

[7]Marsh C A,Backx G,Ernst M H 1997 Phys.Rev.A 56 1676

[8]Chang J Z,Liu M B,Liu H T 2008 Acta Phys.Sin.57 3954(in Chinese)[常建忠,劉謀斌,劉漢濤2008物理學報57 3954]

[9]Visser D C,Hoefsloot H C J,Iedema P D 2006 J.Comput.Phys.214 491

[10]Liu M,Meakin P,Huang H 2006 Phys.Fluids 18 017101

[11]Liu M,Meakin P,Huang H 2007 J.Momput.Phys.2222 110

[12]Wang X L,Chen S 2010 Acta Phys.Sin.59 6778(in Chinese)[王曉亮,陳碩2010物理學報59 6778]

[13]Liu G R,Liu M B 2003 Smoothed Particle Hydrodynamics：A Meshfree Particle Method(Singapore：World Scientif i c)p58

[14]Liu M B,Chang J Z 2010 Acta Phys.Sin.59 7556(in Chinese)[劉謀斌,常建忠2010物理學報59 7556]