劉程 梁宏

(杭州電子科技大學(xué)物理系,杭州 310018)

基于多組分相場(chǎng)理論提出了一類模擬三相流體流動(dòng)的軸對(duì)稱格子Boltzmann 模型.該模型利用兩個(gè)粒子分布函數(shù)來(lái)捕捉三種不同流體之間的相界面,另一個(gè)粒子分布函數(shù)來(lái)求解流體動(dòng)力學(xué)方程以獲得流場(chǎng)信息.為了刻畫坐標(biāo)變換引起的軸對(duì)稱效應(yīng),巧妙地設(shè)計(jì)了演化方程中平衡態(tài)分布函數(shù)和外力項(xiàng)分布函數(shù),從理論上保證本文模型可以正確恢復(fù)三相流體系統(tǒng)的宏觀控制方程,并且軸對(duì)稱效應(yīng)產(chǎn)生的源項(xiàng)中不包含任何復(fù)雜的梯度項(xiàng),從而比現(xiàn)有的軸對(duì)稱格子Boltzmann 模型更加簡(jiǎn)單高效.首先通過(guò)模擬一系列軸對(duì)稱多相流的基準(zhǔn)算例,包括靜態(tài)的雙液滴、液體透鏡的擴(kuò)展和二元流體Rayleigh-Plateau 不穩(wěn)定性,來(lái)驗(yàn)證本文模型的有效性與正確性.接下來(lái),利用該模型研究了三相流體的Rayleigh-Plateau 不穩(wěn)定性的增長(zhǎng)過(guò)程,定量分析了波數(shù)和液柱半徑比對(duì)復(fù)合液體線程破裂過(guò)程中界面動(dòng)力學(xué)行為、界面破裂時(shí)間以及生成子液滴尺寸的影響.可以發(fā)現(xiàn)復(fù)合的液體線程在波數(shù)較大時(shí)破裂生成一個(gè)復(fù)合主液滴和衛(wèi)星液滴,而在波數(shù)較小時(shí)可以生成更多數(shù)量的衛(wèi)星液滴,這導(dǎo)致復(fù)合主液滴和衛(wèi)星液滴的尺寸隨著波數(shù)的增加呈現(xiàn)先增大而后減少的趨勢(shì).另外,我們發(fā)現(xiàn)內(nèi)部流體優(yōu)先于中間流體發(fā)生界面破裂,并且流體界面的破裂時(shí)間均隨著波數(shù)的減小而增加.最后發(fā)現(xiàn),增大液柱半徑比可以促進(jìn)內(nèi)部液體線程的破裂,而延緩中間液體線程發(fā)生破裂,并且復(fù)合主液滴的尺寸隨著液柱半徑比的增加而增大,而復(fù)合衛(wèi)星液滴的尺寸對(duì)液柱半徑比的變化不顯著.

1 引言

三相流體系統(tǒng)廣泛存在于自然界和工程應(yīng)用中,如水污染物的輸運(yùn)、強(qiáng)化采油中油水氣的傳輸過(guò)程、微流控裝置中復(fù)合液滴的生成技術(shù)等.由于涉及多種流體界面的遷徙、變形、破裂及合并等現(xiàn)象,三相流體的界面動(dòng)力學(xué)行為十分復(fù)雜.隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,數(shù)值模擬已成為研究三相流體流動(dòng)問(wèn)題的有效手段.它可以方便地提供實(shí)驗(yàn)研究中不易測(cè)量的物理量,如在相界面變化過(guò)程中某一時(shí)刻的界面形狀、界面處流體的速度、壓力及密度分布等.根據(jù)不同的數(shù)值方法,一些學(xué)者對(duì)三相流體系統(tǒng)進(jìn)行了研究,并在笛卡爾坐標(biāo)系下建立了多種數(shù)值模型,如Smith 等[1]開發(fā)了一種水平投影方法,用于在水平集公式中處理多個(gè)連接的運(yùn)動(dòng),并對(duì)三相流系統(tǒng)進(jìn)行了研究;Bonhomme 等[2]發(fā)展了三相流的流體體積函數(shù)法,用于研究氣泡穿過(guò)兩相流體界面的動(dòng)力學(xué)行為;Kim[3]基于相場(chǎng)理論,提出了包含表面張力的三組分不混溶流體流動(dòng)模型.

格子Boltzmann (lattice Boltzmann,LB) 方法[4]是近三十年來(lái)發(fā)展起來(lái)的流體系統(tǒng)模擬方法,通過(guò)描述粒子分布函數(shù)的演化再現(xiàn)流體的宏觀行為.由于其微觀本質(zhì)和介觀特點(diǎn),LB 方法相比傳統(tǒng)數(shù)值方法具有算法簡(jiǎn)單、天然并行、易于處理流體間的微觀相互作用等優(yōu)點(diǎn),被廣泛用于模擬多相流體[5,6]以及三相流體輸運(yùn)問(wèn)題.Spencer 等[7]基于兩相流LB 方法中顏色模型,重新設(shè)計(jì)了三相流體間的顏色梯度,從而建立了三組分流體的顏色模型.Leclaire 等[8]在顏色梯度的LB 模型中引入了三個(gè)碰撞算子,并使用高階離散算子計(jì)算顏色梯度,從而提出了數(shù)值穩(wěn)定性更好的三相流LB 模型.Yu 等[9]基于Leclaire 等[8]開發(fā)的微擾算子,推導(dǎo)出了描述不同流體之間相互作用的界面力公式,又應(yīng)用了Spencer 等[7]提出的重著色算法來(lái)維持界面并確保流體的不混溶性,提出了用于模擬具有表面張力的不混溶三相流體的顏色梯度LB 模型.Semprebon 等[10]基于自由能泛函理論,推導(dǎo)出流體間和流固間的相互作用,從而建立了表面張力可以獨(dú)立調(diào)節(jié)的三相自由能LB 模型.W?hrwag 等[11]將非理性流體狀態(tài)方程耦合至三相自由能泛函中,并采用了熵的碰撞算子以提高模型的穩(wěn)定性,從而發(fā)展了滿足熱力學(xué)一致性、穩(wěn)定性好的三相流LB 模型.Liang 等[12,13]基于多組分的相場(chǎng)理論,建立了模擬非混相三相流體輸運(yùn)的LB 模型.該模型巧妙地設(shè)計(jì)了平衡態(tài)分布函數(shù)和外力項(xiàng)分布函數(shù),從而理論上保證可以正確地恢復(fù)多組分的Cahn-Hilliard 方程和不可壓Navier-Stokes 方程.

上述模型的提出促進(jìn)了LB 方法在三相流領(lǐng)域的發(fā)展.然而,需要指出的是,這些模型均基于笛卡爾直角坐標(biāo)系建立的.而在三相流體系統(tǒng)中,存在著相界面顯示為軸對(duì)稱特性的問(wèn)題,如兩種不同液滴的對(duì)頭碰撞、氣泡穿過(guò)液-液兩相界面的上升動(dòng)力學(xué)過(guò)程、三相Rayleigh-Plateau 不穩(wěn)定性問(wèn)題等.當(dāng)采用LB 方法來(lái)模擬此類軸對(duì)稱三相流動(dòng)問(wèn)題,最直接的方式是利用三維的三相流LB 模型和合適的曲邊界處理格式.然而,這種方式?jīng)]有利用軸對(duì)稱流動(dòng)的特點(diǎn): 具有軸對(duì)稱特性的流體系統(tǒng)可以將三維流動(dòng)轉(zhuǎn)化為子午面上的二維問(wèn)題,大大減少了計(jì)算量和數(shù)據(jù)儲(chǔ)存.為了充分利用軸對(duì)稱流動(dòng)的特點(diǎn),一些學(xué)者考慮圓柱坐標(biāo)系下建立了有效的多相流LB 模型.然而,絕大多數(shù)的軸對(duì)稱多相模型是基于兩相流情形[14],針對(duì)圓柱坐標(biāo)系下的三相流LB 模型的研究鮮有報(bào)道.Haghani 等[15]從Cahn-Hilliard 相場(chǎng)理論出發(fā),建立了軸對(duì)稱坐標(biāo)下的三相流LB 模型,但該模型中由軸對(duì)稱效應(yīng)引入的源項(xiàng)包含許多復(fù)雜的梯度項(xiàng),這導(dǎo)致模型的算法復(fù)雜且梯度項(xiàng)的離散也會(huì)帶來(lái)額外的數(shù)值誤差.另一方面,針對(duì)三相流體Rayleigh-Plateau 不穩(wěn)定性的研究工作也較少.Yang 等[16]利用軸對(duì)稱的相場(chǎng)算法,對(duì)三相流體的Rayleigh-Plateau 不穩(wěn)定性問(wèn)題開展了初步研究,定性分析了黏性比、界面張力和初始擾動(dòng)對(duì)復(fù)合線程的界面動(dòng)力學(xué)的影響.

鑒于上述分析,本文將基于多組分的相場(chǎng)理論,提出一個(gè)簡(jiǎn)單高效的軸對(duì)稱三相流LB 模型.該模型利用了合適的平衡態(tài)分布函數(shù)和外力項(xiàng)分布函數(shù),從理論上保證可以正確地描述軸對(duì)稱三相流體流動(dòng),并且由軸對(duì)稱效應(yīng)引入的源項(xiàng)不包含任何的梯度項(xiàng),從而比已有軸對(duì)稱三相流LB 模型更加簡(jiǎn)單.通過(guò)模擬靜態(tài)的雙液滴、液體透鏡的擴(kuò)展和二元流體Rayleigh-Plateau 不穩(wěn)定性算例來(lái)驗(yàn)證本文LB 模型的準(zhǔn)確性,并以此模擬了三相流體Rayleigh-Plateau 不穩(wěn)定性的增長(zhǎng)過(guò)程,詳細(xì)分析了波數(shù)和內(nèi)外液柱半徑比對(duì)復(fù)合液體線程的界面動(dòng)力學(xué)行為、界面破裂時(shí)間及生成子液滴尺寸的影響機(jī)制.

2 數(shù)值模型

2.1 宏觀控制方程

本文模型建立在一個(gè)含有三種不互溶流體組成的不可壓縮系統(tǒng),并引入?1,?2,?3三個(gè)序參數(shù)來(lái)描述三種流體在混合物所占的體積分?jǐn)?shù),且滿足如下的約束條件:

根據(jù)多組分流體的相場(chǎng)理論[17,18],三相流體系統(tǒng)的總的自由能泛函可以表述為

其中,Ω為三相系統(tǒng)所處的流域,F(?1,?2,?3) 為體相區(qū)的自由能,D為界面厚度的特征尺度,λi為與表面張力有關(guān)的物理參數(shù).在相場(chǎng)理論[19]中,化學(xué)勢(shì)定義為自由能泛函對(duì)序參數(shù)的變分.然而,針對(duì)三相流體系統(tǒng),為了滿足序參數(shù)的守恒條件,需要在自由能的變分項(xiàng)中引入一個(gè)拉格朗日乘數(shù)β,從而化學(xué)勢(shì)μi可以表述為

序參數(shù)?i的時(shí)間演化由流體速度的對(duì)流和化學(xué)勢(shì)的擴(kuò)散來(lái)描述,其宏觀控制方程可以寫成多組分Cahn-Hilliard 方程的形式:

其中,u為流體速度,Mi為遷移率.假設(shè)遷移率滿足約束條件M1λ1=M2λ2=M3λ3=M0,并對(duì)方程(4)中所有i進(jìn)行求和,可以得到關(guān)于變量S的方程:

其中S=?1+?2+?3.根據(jù)方程(1)給出的守恒條件,S=1 顯然是方程(5)的解,可以推導(dǎo)出

其中,λT定義為

將β的表達(dá)式(6)代入方程(3),可以得到化學(xué)勢(shì)μi的表達(dá)式為

為了使相場(chǎng)模型滿足減少一致性條件以及具有很好的適定性屬性,Boyer 等[17,18]理論構(gòu)造了如下的體相區(qū)自由能:

其中,λi是與表面張力有關(guān)的參數(shù),

其中,σ12,σ13,σ23分別表示三相體系中兩種流體之間的表面張力.將(9)式代入(8)式,最終可以獲得化學(xué)勢(shì)的表達(dá)式:

其中,p是壓力,ρ是流體密度,ν是運(yùn)動(dòng)學(xué)黏性系數(shù),G是外力,Fs是表面張力.另外,為了減少相界面處的虛假速度,本文選用勢(shì)能形式的表面張力

為了獲得軸對(duì)稱三相流體動(dòng)力學(xué)的宏觀控制方程,進(jìn)行了從笛卡爾坐標(biāo)系到柱坐標(biāo)系的變換

其中,x=rcosθ,y=rsinθ,r,z,θ分別表示徑向坐標(biāo),軸向坐標(biāo)和方位角.經(jīng)過(guò)一系列的代數(shù)計(jì)算,可以得到柱坐標(biāo)系下的宏觀控制方程:

2.2 軸對(duì)稱多相流的格子Boltzmann 模型

本文需要利用兩個(gè)粒子分布函數(shù)的LB 演化方程來(lái)捕捉三相流體間的相界面,另一個(gè)粒子分布函數(shù)的LB 方程來(lái)求解流場(chǎng)動(dòng)力學(xué)特性.采用最簡(jiǎn)單的單松弛碰撞模型,追蹤三相流體相界面的LB 演化方程可以表述為[12]

其中,(x,t) 是粒子在t時(shí)刻在x位置的序參數(shù)分布函數(shù),為對(duì)應(yīng)的平衡態(tài)分布函數(shù),τi是無(wú)量綱松弛時(shí)間,δt是時(shí)間增量,(x,t) 是外力項(xiàng)分布函數(shù).將(16)式的最右端項(xiàng)吸收至對(duì)流項(xiàng),從而巧妙設(shè)計(jì)了序參數(shù)的平衡態(tài)分布函數(shù):

其中,ωk和ek分別是權(quán)重系數(shù)和離散速度,cs是聲速,η是引入的一個(gè)自由參數(shù),用來(lái)調(diào)整遷移率系數(shù).采用最具有代表性的D2Q9 離散速度模型[20],其對(duì)應(yīng)權(quán)系數(shù)的取值為ω0=4/9,ω1-4=1/9,,離散速度ek定義為

式中c=δx/δt表示格子速度,δx代表格子步長(zhǎng).為了能夠準(zhǔn)確地恢復(fù)軸對(duì)稱的多組分Cahn-Hilliard方程,還構(gòu)造了全新的外力項(xiàng)分布函數(shù):

在本模型中,通過(guò)計(jì)算粒子分布函數(shù)的零階矩,來(lái)獲得序參數(shù)?i的宏觀量

另外,序參數(shù)?3的值可以通過(guò)守恒條件(1)得到,即?3=1-?1-?2,并且流體的密度ρ可以看成序參數(shù)?i的線性插值函數(shù):

其中,ρ1,ρ2,ρ3分別表示三相流體系統(tǒng)所對(duì)應(yīng)的相密度.利用Chapman-Enskog 多尺度分析技術(shù)[21],可以證明本模型可以準(zhǔn)確地恢復(fù)到三相流體界面追蹤的軸對(duì)稱Cahn-Hilliard 相場(chǎng)方程以及遷移率Mi和松弛因子之間的關(guān)系:

由(24)式可知,當(dāng)遷移率非常小的情況下,無(wú)量綱松弛時(shí)間τi將會(huì)接近 0.5,這可能會(huì)導(dǎo)致LB 方法的數(shù)值不穩(wěn)定,因此引入自由參數(shù)η是非常有必要的,且不會(huì)增加任何額外的計(jì)算成本.

對(duì)于三相不可壓縮流體的建模,需要將多組分Cahn-Hilliard 相場(chǎng)方程與Navier-Stoke 流動(dòng)方程相耦合.基于單松弛碰撞模型,描述流體動(dòng)力學(xué)的的LB 演化方程可以寫成[12]

其中,gk(x,t) 為描述流場(chǎng)的粒子分布函數(shù),(x,t)為平衡態(tài)分布函數(shù),τg為與流體黏度有關(guān)的松弛因子,Gk(x,t) 為外力項(xiàng)分布函數(shù).考慮到軸對(duì)稱效應(yīng)帶來(lái)的影響,為了正確的恢復(fù)連續(xù)性方程,本文采用了如下修正的平衡態(tài)分布函數(shù)[14]:

另外,需要考慮將外力項(xiàng)引入介觀LB 方法中所產(chǎn)生的格子離散效應(yīng).Liang 等[22]近期考慮了格子離散效應(yīng),設(shè)計(jì)了一類準(zhǔn)確并且簡(jiǎn)化的外力項(xiàng)分布函數(shù).在軸對(duì)稱坐標(biāo)系下,流場(chǎng)的外力項(xiàng)分布函數(shù)可以表示為

通過(guò)計(jì)算粒子分布函數(shù)的一階矩和零階矩,可以獲得流體流動(dòng)的速度

和流體動(dòng)力學(xué)壓力p為

利用多尺度理論分析手段[22],可以證明本模型可以恢復(fù)正確的軸對(duì)稱不可壓Navier-Stokes 方程,并且運(yùn)動(dòng)性黏性與松弛因子之間的關(guān)系為

在數(shù)值模擬中,需要采用合適的差分格式對(duì)模型中時(shí)間導(dǎo)數(shù)項(xiàng)和空間導(dǎo)數(shù)項(xiàng)進(jìn)行離散運(yùn)算,本模型采用顯式的歐拉差分格式計(jì)算時(shí)間導(dǎo)數(shù)項(xiàng),

以及各向同性的二階中心差分格式計(jì)算梯度算子和拉普拉斯算子,

其中χ表示任意變量.另外,為了避免軸對(duì)稱LB 模型由于涉及r-1的項(xiàng)而在r=0 處產(chǎn)生奇異現(xiàn)象,將第一層網(wǎng)格放在r=0.5δx處,并對(duì)軸向邊界處應(yīng)用對(duì)稱邊界條件[14].

3 數(shù)值結(jié)果和討論

本節(jié)首先通過(guò)模擬幾個(gè)典型的軸對(duì)稱多相流問(wèn)題來(lái)驗(yàn)證本文提出的介觀LB 模型,包括靜態(tài)的雙液滴、液體透鏡和二元流體Rayleigh-Plateau 不穩(wěn)定性問(wèn)題.接著,利用該模型來(lái)研究三相流體的Rayleigh-Plateau 不穩(wěn)定性問(wèn)題,詳細(xì)分析波數(shù)和內(nèi)外液柱的半徑比對(duì)液線破裂過(guò)程中界面運(yùn)動(dòng)過(guò)程以及產(chǎn)生的復(fù)合液滴尺寸的影響.

3.1 靜態(tài)的雙液滴

靜態(tài)的雙液滴作為三相流的基準(zhǔn)算例,首先被用來(lái)驗(yàn)證本文所建立的軸對(duì)稱三相流 LB 模型.在靜態(tài)的雙液滴測(cè)試中,將計(jì)算網(wǎng)格設(shè)置為Nz×Nr=400×100,其中尺寸相同的兩個(gè)液滴靜止放置于計(jì)算區(qū)域中,液滴的半徑為R=50,兩液滴的中心坐標(biāo)分別是 (z1,r0)=(100,0) 和 (z2,r0)=(300,0).初始的序參數(shù)分布設(shè)定為

而序參數(shù)?3可以由方程(1)得到,其他物理參數(shù)設(shè)置為: 三相系統(tǒng)密度ρ1=10,ρ2=5,ρ3=1,界面厚度D=4,運(yùn)動(dòng)學(xué)黏性ν=0.1,界面張力σ12=σ13=σ23=0.1,松弛時(shí)間τ1=τ2=0.8,遷移率M0=0.01.粒子分布函數(shù)fi和gi的邊界條件設(shè)置如下: 左右邊界運(yùn)用周期邊界條件,上壁面采用無(wú)滑移反彈邊界條件,下邊界則采用對(duì)稱性邊界條件.圖1 給出了穩(wěn)態(tài)時(shí)三個(gè)序參數(shù)在計(jì)算區(qū)域中的分布.可以看出,序參數(shù)的穩(wěn)態(tài)分布符合初始條件給定的雙曲分布,避免了三相流的數(shù)值方法可能會(huì)產(chǎn)生非物理的虛假相.進(jìn)一步統(tǒng)計(jì)了序參數(shù)沿著界面方向的穩(wěn)態(tài)分布以及相應(yīng)的解析解,并將結(jié)果列于圖2.從圖中可以觀察到數(shù)值預(yù)測(cè)的序參數(shù)分布與解析解能很好地符合,這表明本文提出軸對(duì)稱三相流的LB 模型能夠準(zhǔn)確捕獲界面的分布.

圖1 靜態(tài)雙液滴測(cè)試中三個(gè)序參數(shù)的穩(wěn)態(tài)分布 (a) ?1 的分布;(b) ?2 的分布;(c) ?3 的分布Fig.1.Steady-state distributions of three order parameters in the static double-droplet test: (a) Distribution of ?1 ;(b) distribution of ?2 ;(c) distribution of ?3.

圖2 靜態(tài)雙液滴測(cè)試中序參數(shù)沿水平線的穩(wěn)態(tài)分布Fig.2.Steady distributions of the order parameters across the horizontal line in static double-droplet test.

3.2 液體透鏡的擴(kuò)展

液體透鏡的擴(kuò)展問(wèn)題被廣泛用于測(cè)試所提出的三相流LB 模型,然而絕大多數(shù)工作限制于二維笛卡爾坐標(biāo)系[9,10,12],針對(duì)三維或柱坐標(biāo)空間的模擬鮮有報(bào)道.本節(jié)將模擬三維液體透鏡的擴(kuò)展來(lái)驗(yàn)證所提出的三相流軸對(duì)稱LB 模型.如圖3 所示,圓形的液體透鏡初始時(shí)刻放置于另外兩種不混溶的流體之間,在表面張力作用下會(huì)發(fā)生變形,并最終達(dá)到穩(wěn)定的形狀.根據(jù)Neumann定律[23],液體透鏡在平衡時(shí)的形狀由三個(gè)表面張力系數(shù)決定:

圖3 液體透鏡相關(guān)幾何參數(shù)示意圖Fig.3.Schematic diagram of relevant geometric parameters of liquid lens.

其中,θ1,θ2表示接觸角.另外,透鏡的理論長(zhǎng)度d可以表示為

其中,A是初始時(shí)刻圓形透鏡的面積.液體透鏡的高度可以解析地表示為

該問(wèn)題的初始條件設(shè)置如下: 半徑R=100 的圓形透鏡坐落在Nz ×Nr=500×250 計(jì)算區(qū)域內(nèi),其中心點(diǎn)的坐標(biāo)為 (zc,rc)=(250, 0),從而序參數(shù)的初始分布可以設(shè)定為

其他相關(guān)的物理參數(shù)給定為ρ1=10,ρ2=1,ρ3=5,界面厚度D=4,運(yùn)動(dòng)學(xué)黏度ν=0.1,松弛時(shí)間τ1=τ2=0.8,遷移率M0=0.01.邊界條件采用與靜態(tài)的雙液滴測(cè)試相同的邊界處理格式.圖4 給出了三種不同表面張力比σ12:σ13:σ23=0.6:0.6:1,1:1:1, 1:4/3:1 時(shí),液體透鏡到達(dá)系統(tǒng)平衡狀態(tài)具有的界面形狀.從圖中可以發(fā)現(xiàn),液體透鏡在不同的表面張力比時(shí)具有不同穩(wěn)態(tài)界面形狀,其界面形狀與文獻(xiàn)報(bào)道的數(shù)值結(jié)果[12]定性上一致.進(jìn)一步測(cè)量了液體透鏡在平衡時(shí)的接觸角θ1和θ2,并將數(shù)值預(yù)測(cè)的計(jì)算結(jié)果與相應(yīng)的解析解總結(jié)于表1中.從表1 可以發(fā)現(xiàn),LB 方法所預(yù)測(cè)接觸角的數(shù)值結(jié)果與解析解的最大相對(duì)誤差不超過(guò)1.7%,即數(shù)值結(jié)果與理論解相符合.本文還計(jì)算了透鏡的長(zhǎng)度d和高度 (h1,h2) 以及相應(yīng)的解析解,并將結(jié)果列于表2 中.從表2 可以發(fā)現(xiàn),當(dāng)表面張力比為(0.6:0.6:1) 時(shí),長(zhǎng)度d的數(shù)值解與解析解存在著較大誤差,而當(dāng)表面張力比取其他情形,數(shù)值預(yù)測(cè)的結(jié)果與解析解相比大體上可以接受.這種誤差較高的結(jié)果也存在于前人的軸對(duì)稱或三維的LB 方法模擬[15,24]中.可以通過(guò)采用更精細(xì)的網(wǎng)格或更高階的差分格式來(lái)計(jì)算梯度和拉普拉斯項(xiàng)以降低計(jì)算誤差.

表1 不同表面張力比的液體透鏡的平衡接觸角θ1 和 θ2Table 1.Equilibrium contact angles θ1 and θ2 of liquid lens at different surface tension ratios.

表2 不同表面張力比液體透鏡的長(zhǎng)度d 和高度 h1,h2Table 2.Length d and heights h1,h2 of liquid lens at different surface tensions.

圖4 不同界面張力比液體透鏡的平衡形態(tài) (a) σ12 :σ13 :σ23=0.6:0.6:1 ;(b) σ12 :σ13 :σ23=1:1:1 ;(c) σ12 :σ13 :σ23=1:4/3:1.Fig.4.Equilibrium morphology of liquid lens at different interfacial tension ratios: (a) σ12 :σ13 :σ23=0.6:0.6:1 ;(b) σ12 :σ13 :σ23=1:1:1 ;(c) σ12 :σ13 :σ23=1:4/3:1.

3.3 Rayleigh-Plateau 不穩(wěn)定性

Rayleigh-Plateau 不穩(wěn)定性問(wèn)題廣泛存在于日常生活和實(shí)際工程中,如水龍頭中流出的射流、噴墨打印和基因芯片排列等.在不穩(wěn)定性的演化過(guò)程中,其界面動(dòng)力學(xué)行為呈現(xiàn)出軸對(duì)稱的性質(zhì),因此,本文提出的軸對(duì)稱多相流LB 模型適合用于模擬經(jīng)典的Rayleigh-Plateau 不穩(wěn)定性問(wèn)題.根據(jù)文獻(xiàn)調(diào)研可知[14],已有的研究多數(shù)著眼于兩相情形,而本文的軸對(duì)稱三相流模型LB 模型理論上也滿足減少一致性條件,因而,我們首先模擬二元流體的Rayleigh-Plateau 不穩(wěn)定性問(wèn)題來(lái)驗(yàn)證數(shù)值模型的正確性.計(jì)算的網(wǎng)格分辨率設(shè)定為Nz×Nr=λ×200,其中λ表示波長(zhǎng).序參數(shù)的初始分布設(shè)置如下:

其中,R=60 為液柱的半徑,d=0.1Rcos(2πz/λ)為兩相界面處施加的微小擾動(dòng).為了考察波數(shù)k(k=2πR/λ) 對(duì)不穩(wěn)定性界面動(dòng)力學(xué)行為的影響,本文考慮三種不同的波長(zhǎng),分別設(shè)定為λ=600,1200 和1800,其他物理參數(shù)為ρ1=10,ρ3=1,σ13=0.1,D=4,M0=0.1,τ1=τ2=1.0,ν=0.1.另外選擇毛細(xì)時(shí)間為特征時(shí)間,且如下的演化時(shí)間均被特征值無(wú)量綱化.圖5 給出了不同波數(shù)下Rayleigh-Plateau 不穩(wěn)定性的界面動(dòng)力學(xué)行為.從圖中可以發(fā)現(xiàn),界面處不穩(wěn)定性的擾動(dòng)隨著時(shí)間逐漸發(fā)展,導(dǎo)致中間區(qū)域的液體線程變得越來(lái)越纖細(xì),而兩端的液體尺寸逐漸增大.經(jīng)過(guò)一段時(shí)間的演化,液柱界面在兩個(gè)夾點(diǎn)發(fā)生破裂,形成了一段細(xì)長(zhǎng)的液體燈絲和一個(gè)離散的主液滴.接下來(lái),不穩(wěn)定性在不同的波數(shù)下顯示出不同的界面動(dòng)力學(xué)特征.當(dāng)波數(shù)k=0.63 或0.31 時(shí),液體燈絲在表面張力的作用下逐漸收縮,并最終形成一個(gè)衛(wèi)星液滴.而當(dāng)波數(shù)較小k=0.21 時(shí),液體燈絲最初也隨時(shí)間而逐漸收縮,并伴隨著燈絲兩側(cè)的體積逐漸增大,從而在界面表面又形成了二次Rayleigh-Plateau 不穩(wěn)定性,導(dǎo)致了液體線程發(fā)生二次破裂,最終釋放出三個(gè)衛(wèi)星液滴.本文計(jì)算獲得不同波數(shù)下Rayleigh-Plateau 不穩(wěn)定性的界面動(dòng)力學(xué)行為與前人的研究結(jié)果[14,25]在定性上是相符合的.

圖5 不同波數(shù)下兩相Rayleigh-Plateau 不穩(wěn)定性的界面演化過(guò) 程 (a) k=0.63 ;(b) k=0.31 ;(c) k=0.21.Fig.5.Interfacial evolution processes of two-phase Rayleigh-Plateau instability with different wavenumbers:(a) k=0.63 ;(b) k=0.31 ;(c) k=0.21.

接下來(lái)采用本文提出的軸對(duì)稱三相流LB 模型研究三相流體的Rayleigh-plateau 不穩(wěn)定性問(wèn)題,重點(diǎn)分析波數(shù)和內(nèi)部液體與中間液體半徑比對(duì)三相流體界面動(dòng)力學(xué)行為、界面破裂時(shí)間及生成子液滴尺寸的影響.根據(jù)三相系統(tǒng)中流體所處的位置,將三種流體分別稱為內(nèi)部液體、中間液體和外部液體.計(jì)算網(wǎng)格設(shè)定為Nz ×Nr=λ×2 00,邊界條件與上一個(gè)算例相同.序參數(shù)的初始分布給定如下:

3.3.1 波數(shù)的影響

首先考察波數(shù)對(duì)三相流體Rayleigh-Plateau不穩(wěn)定性演化的影響,并通過(guò)調(diào)節(jié)初始擾動(dòng)的波長(zhǎng)λ來(lái)獲得不同的波數(shù)k=2πRI/λ,其中波長(zhǎng)λ的值分別為600,900,1200,1500,1800.圖6給出了不同波數(shù)下三相Rayleigh-Plateau 不穩(wěn)定性的界面演化過(guò)程.從圖中可以觀察到,不穩(wěn)定性在初始階段表現(xiàn)出相似的界面動(dòng)力學(xué)特征: 界面隨時(shí)間發(fā)生變形,其擾動(dòng)的振幅逐漸增大,導(dǎo)致中間區(qū)域的液體線程變得越來(lái)越細(xì)薄,而兩端的液體體積逐漸增大.另外可以發(fā)現(xiàn),內(nèi)部液體的表面擾動(dòng)也隨時(shí)間逐漸增長(zhǎng),并經(jīng)歷了更大的界面變形,從而在波數(shù)較大(k=0.42 或0.21)時(shí),優(yōu)先發(fā)生界面破裂形成了液體燈絲和主液滴結(jié)構(gòu),而在波數(shù)較小(k=0.14)時(shí),在界面表面上形成了多模態(tài)的擾動(dòng).隨后,當(dāng)波數(shù)k=0.42 時(shí),內(nèi)部的液體燈絲在表面張力作用下首先逐漸收縮,在界面表面上形成了二次擾動(dòng),并在二次Rayleigh-Plateau 不穩(wěn)定性的作用下,形成的液體燈絲會(huì)再次發(fā)生斷裂,分離出兩個(gè)衛(wèi)星液滴.衛(wèi)星液滴隨著中間液體的斷裂,分別向位于兩端的主液滴移動(dòng)并與主液滴的內(nèi)部流體發(fā)生合并,從而最終形成了位于兩端的復(fù)合主液滴和由中間液體斷裂形成的衛(wèi)星液滴,并且復(fù)合主液滴表現(xiàn)為中間液體包裹內(nèi)部液體的形態(tài).當(dāng)波數(shù)k=0.21 時(shí),內(nèi)部的液體燈絲會(huì)在多個(gè)位置發(fā)生二次破裂,從而形成了更多數(shù)量的衛(wèi)星液滴.緊接著,兩測(cè)的衛(wèi)星液滴向主液滴方向移動(dòng),與主液滴中內(nèi)部流體接觸而合并,最終形成了復(fù)合的主液滴.而中間液體界面斷裂后形成的液體燈絲受表面張力作用逐漸收縮,致使余下的衛(wèi)星液滴向中間聚攏并融合成單一液滴,并最終形成了復(fù)合的衛(wèi)星液滴.當(dāng)波數(shù)充分小(k=0.14)時(shí),多模態(tài)的擾動(dòng)隨后逐漸發(fā)展,并使中間流體和內(nèi)部液體的燈絲在多個(gè)位置發(fā)生斷裂,最終生成了一個(gè)復(fù)合主液滴和三個(gè)復(fù)合的衛(wèi)星液滴,這表明隨著波數(shù)的減小,三相流體Rayleigh-Plateau 不穩(wěn)定性可以誘導(dǎo)生成更多的復(fù)合子液滴,這與二元Rayleigh-Plateau 不穩(wěn)定性的結(jié)果[26]在定性上是一致的.另外還統(tǒng)計(jì)了波數(shù)對(duì)界面破裂時(shí)間的影響,發(fā)現(xiàn)不同波數(shù)(k=0.42,0.28, 0.21, 0.17, 0.14)下內(nèi)部流體線程破裂的無(wú)量綱時(shí)間分別為10.6,15.4,20.0,23.8,26.3,而相應(yīng)的中間流體線程的無(wú)量綱破裂時(shí)間為21.9,27.5,34.4,38.8,40.0,這表明在三相Rayleigh-Plateau不穩(wěn)定性的演化中,內(nèi)部流體優(yōu)先于中間流體發(fā)生界面破裂,并且內(nèi)外流體界面的破裂時(shí)間均隨著波數(shù)的減小而增加.

圖6 不同波數(shù)下三相Rayleigh-Plateau 不穩(wěn)定性的界面演化過(guò)程 (a) k=0.42 ;(b) k=0.21 ;(c) k=0.14Fig.6.Interfacial evolution processes of three-phase Rayleigh-Plateau instability with different wavenumbers:(a) k=0.42 ;(b) k=0.21 ;(c) k=0.14.

進(jìn)一步,研究了波數(shù)對(duì)復(fù)合液體線程破裂生成子液滴尺寸的影響,統(tǒng)計(jì)結(jié)果如圖7 所示,其中R*表示生成子液滴的無(wú)量綱半徑,其值由生成的復(fù)合液滴的半徑與內(nèi)部液體的初始半徑RI的比值得到.由于在低波數(shù)情形,不穩(wěn)定性可以誘導(dǎo)生成多個(gè)衛(wèi)星液滴,因此本文僅考慮系統(tǒng)中尺寸最大的復(fù)合主液滴和尺寸最小的衛(wèi)星液滴.從圖7 可以發(fā)現(xiàn),復(fù)合主液滴和衛(wèi)星液滴的尺寸均隨著波數(shù)的增大呈現(xiàn)出先增長(zhǎng)而后減少的趨勢(shì).這是由于在波數(shù)非常小時(shí),不穩(wěn)定性的波長(zhǎng)較大,液體線程發(fā)生破裂時(shí)產(chǎn)生了更多數(shù)目的子液滴,從而根據(jù)液體體積守恒,導(dǎo)致生成復(fù)合主液滴和衛(wèi)星液滴的尺寸較小.當(dāng)波數(shù)較大時(shí),液體線程發(fā)生破裂生成一個(gè)主液滴和衛(wèi)星液滴,其尺寸隨著波數(shù)的增加而減小,該結(jié)果與二元流體Rayleigh-Plateau 不穩(wěn)定性的結(jié)果[14,25,26]一致.

圖7 不同波數(shù)下復(fù)合液體線程破裂生成主液滴和衛(wèi)星液滴的尺寸Fig.7.The sizes of main and satellite droplets originating from the breakup of compound liquid threads with different wavenumbers.

3.3.2 內(nèi)外液柱半徑比的影響

本節(jié)研究?jī)?nèi)部液體與中間液體半徑比對(duì)三相流體Rayleigh-Plateau 不穩(wěn)定性演化的影響.計(jì)算區(qū)域設(shè)置為Nz ×Nr=1200×200,其中內(nèi)部液體的半徑RI的值固定為40,而中間液體的半徑RM的值分別取 60,70,80,90,100.圖8 給出了不同內(nèi)外液柱半徑比下三相流體Rayleigh-Plateau 不穩(wěn)定性的界面演化過(guò)程.從圖中可以發(fā)現(xiàn),三元Rayleigh-Plateau 不穩(wěn)定性在不同的液柱半徑比下表現(xiàn)出不同的相界面動(dòng)力學(xué)過(guò)程.當(dāng)液柱半徑比為RM/RI=1.5 時(shí),不穩(wěn)定性的界面擾動(dòng)在初始階段隨時(shí)間逐漸發(fā)展,導(dǎo)致中間區(qū)域的線程變薄而兩端部分不斷增大.內(nèi)部流體經(jīng)歷更大的界面變形,從而優(yōu)先于中間流體發(fā)生破裂,形成了主液滴和液體燈絲.緊接著,液體燈絲的末尾也發(fā)生斷裂釋放出兩個(gè)衛(wèi)星液滴.衛(wèi)星液滴向主液滴方向移動(dòng),與主液滴發(fā)生合并.與此同時(shí),中間流體的線程也隨時(shí)間的演化而越來(lái)越薄,并發(fā)生破裂形成了復(fù)合的主液滴和液體燈絲.在復(fù)合液體燈絲收縮的過(guò)程中,相界面處形成了微小的擾動(dòng),進(jìn)而發(fā)生二次不穩(wěn)定性現(xiàn)象,導(dǎo)致相界面發(fā)生破裂生成了三個(gè)衛(wèi)星液滴.由于慣性作用,兩側(cè)的復(fù)合衛(wèi)星液滴逐漸靠近中間位置的復(fù)合衛(wèi)星液滴,在表面張力的作用下最終合并為一個(gè)體積更大的復(fù)合衛(wèi)星液滴.當(dāng)增大液柱半徑比RM/RI至2 或2.5 時(shí),可以發(fā)現(xiàn)內(nèi)部流體仍能優(yōu)先于中間流體發(fā)生破裂,形成了主液滴和液體燈絲.然而,不同于低液柱半徑比情形,液體燈絲隨后逐漸收縮,受到中間流體的界面變形的影響,發(fā)生破裂形成了被中間流體包裹的多個(gè)子液滴.最外側(cè)的兩個(gè)子液滴向兩側(cè)運(yùn)動(dòng),并與主液滴接觸后發(fā)生融合.緊接著,中間流體的界面也發(fā)生破裂,形成了復(fù)合的主液滴和包裹著多個(gè)子液滴的復(fù)合液體燈絲.在表面張力作用下,復(fù)合液滴燈絲逐漸收縮,多個(gè)子液滴同時(shí)也向中間位置移動(dòng)并發(fā)生合并,最終形成了一個(gè)復(fù)合的衛(wèi)星液滴.此外,還研究了內(nèi)外液柱半徑比對(duì)三相Rayleigh-Plateau 不穩(wěn)定性界面破裂時(shí)間的影響,結(jié)果表明內(nèi)部流體線程在不同半徑比下(RM/RI=1.5, 1.75,2.0, 2.25, 2.5)發(fā)生破裂的無(wú)量綱時(shí)間分別為 21.0,20.5,20.0,19.4,18.1,而中間流體線程破裂所對(duì)應(yīng)的無(wú)量綱時(shí)間為26.3,30.6,34.4,39.4,44.0.可以發(fā)現(xiàn)內(nèi)部液體線程仍然比中間液體線程更早發(fā)生破裂,并且增大液柱半徑比可以促進(jìn)內(nèi)部液體線程的破裂,而延緩中間液體線程發(fā)生破裂,這是由于中間流體在較大的液柱半徑比下占比更多的體積,需要更大的界面變形直至破裂.

圖8 不同內(nèi)外液柱半徑比下三相Rayleigh-Plateau 不穩(wěn)定性的界面演化過(guò)程 (a) RM/RI=1.5 ;(b) RM/RI=2 ;(c) RM/RI=2.5Fig.8.Interfacial evolution of three-phase Rayleigh-Plateau instability with different radius ratios of the inner and outer liquid columns: (a) RM/RI=1.5 ;(b) RM/RI=2 ;(c) RM/RI=2.5.

最后,統(tǒng)計(jì)了不同液柱半徑比下,液體線程破裂生成復(fù)合的主液滴和衛(wèi)星液滴的尺寸,結(jié)果如圖9 所示.從圖中可以發(fā)現(xiàn),復(fù)合主液滴的尺寸隨著液柱半徑比的增加而增大,而復(fù)合衛(wèi)星液滴的尺寸隨著液柱半徑比的增大變化不顯著,這表明中間流體體積增大的部分隨著不穩(wěn)定性界面的變形均往復(fù)合的主液滴方向運(yùn)動(dòng).

圖9 不同液柱半徑比下復(fù)合液體線程破裂生成復(fù)合主液滴和衛(wèi)星液滴的尺寸Fig.9.The sizes of the compound main and satellite droplets originating from the breakup of compound liquid threads with different radius ratios.

4 結(jié)論

本文基于多組分Cahn-Hilliard 相場(chǎng)理論,建立了一類簡(jiǎn)單高效準(zhǔn)確的軸對(duì)稱三相流LB 模型.該模型利用了三個(gè)粒子分布函數(shù): 兩個(gè)序參數(shù)分布函數(shù)追蹤三相流體的相界面,另一個(gè)密度分布函數(shù)求解流體的速度場(chǎng)和壓力場(chǎng)分布.通過(guò)引入修正的平衡態(tài)分布函數(shù)和合適的源項(xiàng),本文模型可以正確地描述軸對(duì)稱三相流體系統(tǒng),并且不同于前人的軸對(duì)稱三相LB 模型,本文模型中由軸對(duì)稱效應(yīng)引起的源項(xiàng)不包含任何復(fù)雜的梯度項(xiàng),從而在算法實(shí)現(xiàn)上更加簡(jiǎn)單高效.為了驗(yàn)證所提出的軸對(duì)稱三相流LB 模型,模擬了一系列經(jīng)典的軸對(duì)稱多相流算例,包括靜態(tài)的雙液滴、液體透鏡的擴(kuò)展和二元流體Rayleigh-Plateau 不穩(wěn)定性.數(shù)值結(jié)果顯示本文模型可以準(zhǔn)確地捕捉三相流體的相界面,預(yù)測(cè)的不同界面張力比下液體透鏡的穩(wěn)態(tài)形狀和Rayleigh-Plateau 不穩(wěn)定性中界面動(dòng)力學(xué)過(guò)程與解析解或文獻(xiàn)報(bào)道數(shù)據(jù)一致.進(jìn)一步采用提出的LB 算法研究了三相流體Rayleigh-Plateau 不穩(wěn)定性問(wèn)題,詳細(xì)分析了波數(shù)和液柱半徑比對(duì)不穩(wěn)定性的界面演化過(guò)程、線程破裂時(shí)間以及子液滴尺寸的影響.可以發(fā)現(xiàn)隨著波數(shù)的減少,不穩(wěn)定性可以誘導(dǎo)生成更多的復(fù)合子液滴.不同于兩相Rayleigh-Plateau 不穩(wěn)定性情形,液體線程破裂形成子液滴的尺寸隨著波數(shù)的增加先增大而后減少.另外,增大液柱半徑比可以促進(jìn)復(fù)合主液滴尺寸的增長(zhǎng),而對(duì)復(fù)合衛(wèi)星液滴尺寸的影響較小.