于蔚然，周 訓(xùn)，李維仲

(大連理工大學(xué) 海洋能源利用與節(jié)能教育部重點(diǎn)實(shí)驗(yàn)室，大連 116024)

1 引 言

在自然界及工程界普遍存在的氣液兩相流動(dòng)現(xiàn)象中，兩相流因其復(fù)雜的流體系統(tǒng)和氣液相界面的不斷變換，是一項(xiàng)具有挑戰(zhàn)性的研究工作，因此該問題備受科學(xué)研究者的關(guān)注。氣流通過界面剪切力驅(qū)動(dòng)的液膜流動(dòng)是一種典型的兩相流動(dòng)形式，也是許多傳熱傳質(zhì)設(shè)備中常見的流動(dòng)形態(tài)。在航空發(fā)動(dòng)機(jī)的噴霧裝置中，帶有預(yù)膜板的氣動(dòng)霧化噴嘴在進(jìn)行霧化的過程中，高速流動(dòng)的氣體會(huì)通過界面剪切力驅(qū)動(dòng)液膜在預(yù)膜板上流動(dòng)并產(chǎn)生波浪形的氣液相界面變化，且氣液層有較大的密度比，同時(shí)氣體雷諾數(shù)也明顯高于液體。研究表明，氣流的剪切作用對(duì)液膜鋪展及流動(dòng)形態(tài)有著重大影響。由于氣液在預(yù)膜板上的流動(dòng)形態(tài)直接影響預(yù)膜板唇邊后部的液體霧化質(zhì)量，因此，探尋不同氣體速度對(duì)液膜流動(dòng)形態(tài)的影響規(guī)律，對(duì)于后期液膜破碎機(jī)理的研究有著重要意義。針對(duì)氣體速度對(duì)液膜流動(dòng)的影響，已有研究結(jié)果顯示[1,2],氣流的加速會(huì)導(dǎo)致液膜的傳播速度增強(qiáng)，進(jìn)而影響液膜界面的變化。本文將針對(duì)該問題進(jìn)行深入的探尋并總結(jié)相應(yīng)的規(guī)律。

關(guān)于液膜流動(dòng)問題，已有諸多理論分析[3,4]、實(shí)驗(yàn)[5,6]及數(shù)值模擬對(duì)其進(jìn)行相關(guān)研究。但對(duì)于相界面的追蹤，數(shù)值模擬方法的優(yōu)越性尤為明顯。目前對(duì)于多相流的數(shù)值模擬方法較為常用的是VOF和Level Set方法，此類方法已成熟地運(yùn)用在液膜流動(dòng)、液滴運(yùn)動(dòng)和溝流等模擬[7,8]上。這類流動(dòng)屬于自由表面的流動(dòng)，周圍氣體對(duì)流動(dòng)的影響微乎其微，往往在數(shù)值模擬中較好處理。但是由于較高氣流剪切作用下氣液相界面會(huì)出現(xiàn)復(fù)雜的拓?fù)湫巫?，VOF和Level Set難于準(zhǔn)確追蹤界面位置，且現(xiàn)有方法很難處理氣液界面處兩相大密度比問題，因而找到一種可以合理處理相關(guān)問題的方法顯得尤為重要。格子玻爾茲曼LBM(Lattice Boltzmann method)[9]，一種基于分子動(dòng)理論的模擬流體流動(dòng)的數(shù)值方法,在經(jīng)過近些年的不斷發(fā)展和完善后，目前已廣泛應(yīng)用于多相流、湍流、相變傳熱和微通道驅(qū)替[10-12]等多種復(fù)雜物理現(xiàn)象的模擬。發(fā)現(xiàn)LBM在追蹤具有較大拓?fù)湫巫兊南嘟缑嫜莼^程中有著獨(dú)特的優(yōu)勢(shì)。針對(duì)用于多相流的LBM模型，現(xiàn)有模型大多只能處理中小密度比的兩相流體流動(dòng)[13-16]。在涉及大密度比的多相LBM模型中，Liang等[17]提出的基于相場(chǎng)理論的LBM模型在降低假速度和保持質(zhì)量守恒方面有著出色表現(xiàn)。從理論上講，該模型的界面捕捉方程采用的是AC(Allen-Cahn)方程，由于該方程具有比CH(Cahn-Hilliard)方程低階的擴(kuò)散項(xiàng)，因此其在計(jì)算序參數(shù)及密度場(chǎng)等方面具有比基于CH方程模型[18]更高的數(shù)值精度和穩(wěn)定性，同時(shí)在界面捕獲中也有相應(yīng)的優(yōu)勢(shì)，因而選擇該模型處理本文提出的問題不失為良好的選擇。但在使用Liang等[17]模型進(jìn)行計(jì)算時(shí)，由于基于相場(chǎng)理論的LBM模型在兩相界面附近存在密度變化，該處的連續(xù)性方程不滿足不可壓縮性條件。因而本文在上述模型中引入額外的界面力來消除此影響。

概括來看，本文首先對(duì)Liang等[17]提出的LBM模型進(jìn)行修正，從而提升模型計(jì)算準(zhǔn)確性。進(jìn)而將該模型用以模擬計(jì)算本文所研究的高密度比氣液兩相流問題，探尋氣體速度對(duì)液膜流動(dòng)形態(tài)的影響，最后總結(jié)規(guī)律。

2 數(shù)學(xué)方法

2.1 基于相域理論的LBM兩相流模型

2.1.1 控制方程

本文所采用模型包含兩組LB方程，一組用于求解界面捕捉的AC方程，一組用于求解不可壓縮的Naiver-Stokes(NS)方程。其中，AC方程[19]為

(1)

式中n為界面處的單位法向量，

(2)

而對(duì)于式(2)的λ，其定義式為

(3)

對(duì)于不可壓縮流體，帶外力項(xiàng)[20]的NS方程可表示為

·u=0

(4a)

(4b)

式中F為外力項(xiàng)總和，詳見下文。

2.1.2 AC方程的LBM演化方程

基于BGK假設(shè)，AC方程的LBM演化方程可寫為

(5)

(6)

式中ωi為權(quán)系數(shù)，其數(shù)值取決于所選擇的格子玻爾茲曼模型。對(duì)于本文所選D2Q9模型，ωi的取值為ω0=4/9,ω1 - 4=1/9,ω5 - 8=1/36，離散速度ci為

(7)

在式(5)中，外力源項(xiàng)Fi的定義式為

(8)

(9)

(10)

(11)

2.1.3 NS方程的LBM演化方程

含外力項(xiàng)的NS方程的格子Boltzmann-BGK演化方程可以表示為[22]

δtGi(x,t)

(12)

(13)

(14)

式(12)的外力源項(xiàng)Gi與現(xiàn)有LB模型中[23-25]的表達(dá)式不同，本文采用更為簡(jiǎn)單的形式：

(15)

式中F為總的外力項(xiàng)。

F=Fs+Fa+G

(16)

式(16)為計(jì)算中可能出現(xiàn)的體積力，F(xiàn)s為表面張力。

(17)

式中β和k的取值與界面厚度以及表面張力系數(shù)σ有關(guān)。

(18)

基于相場(chǎng)理論的LBM模型在兩相界面附近存在密度變化，該處的連續(xù)性方程不滿足不可壓縮性條件，為了消除此影響，本文引入一個(gè)額外的界面力Fa=qau來提高模型的準(zhǔn)確性，其中參數(shù)qa的計(jì)算式為

(19)

在此基礎(chǔ)上，宏觀速度和壓力的計(jì)算公式修正如下：

(20)

(21)

另外，通過Chapman-Enskog展開，流體運(yùn)動(dòng)粘度由式(22)確定。

(22)

要注意的是，在兩相流體系中，運(yùn)動(dòng)粘度系數(shù)是在變化的，為了使其在界面處平滑過渡，通常的處理辦法是將運(yùn)動(dòng)粘度系數(shù)當(dāng)作序參數(shù)的線性或逆線性函數(shù)[15,26]?？紤]到更簡(jiǎn)單的情況，本文采取了線性形式計(jì)算該參數(shù)，

(23)

式中υl和υg分別為液氣兩相的運(yùn)動(dòng)粘度系數(shù)。

3 模型驗(yàn)證

(24)

模擬過程中，計(jì)算域上下壁面為無滑移邊界，左右壁面為周期邊界，網(wǎng)格量設(shè)定為256×1024。

圖2為Re =2048，At=0.5時(shí)，由R-T不穩(wěn)定性引發(fā)的兩相界面隨時(shí)間的演化過程。可以看出，由于初始時(shí)刻兩相界面存在擾動(dòng)，在重力作用下，密度大的流體往下沉(spike)，密度小的流體向上浮(bubble)。當(dāng)兩種流體間的剪切速度差足夠大時(shí)，則會(huì)觸發(fā)Kelvin-Helmholtz不穩(wěn)定性現(xiàn)象，從而使相界面處出現(xiàn)卷曲變形。圖2相界面的演化過程與文獻(xiàn)[27,28]的計(jì)算結(jié)果基本一致。

通過以上結(jié)果證實(shí)了該模型在低密度比下捕獲R-T不穩(wěn)定性特征的能力，但對(duì)于修正后的模型在計(jì)算高雷諾數(shù)下具有大密度比的問題，其可靠性沒有得到充分的驗(yàn)證，因此接下來本文將進(jìn)一步對(duì)該模型進(jìn)行驗(yàn)證計(jì)算。文獻(xiàn)[17,29]指出，在R-T不穩(wěn)定性發(fā)生初期，界面擾動(dòng)呈指數(shù)型增長(zhǎng)，增長(zhǎng)規(guī)律為

a=a0eα t

(25)

圖1 R-T不穩(wěn)定性示意圖

圖2 Re=2048,At=0.5時(shí)R-T不穩(wěn)定性現(xiàn)象中界面形態(tài)隨時(shí)間的演化

4 模擬結(jié)果與討論

為了探究在高速氣流驅(qū)動(dòng)下液膜在水平預(yù)膜板表面的流動(dòng)狀態(tài)，構(gòu)建如圖4所示的兩相流體系，氣體層和液體層厚度分別為Hg和Hl，氣液初始速度分別為Ug和Ul，且Ug>Ul。液體進(jìn)口速度給定為

ux(y)=Ultanh(y/δ)

(26)

(27)

圖3 不同k*下增長(zhǎng)率α*的計(jì)算結(jié)果與文獻(xiàn)[17]和

解析解[29]的結(jié)果對(duì)比

Fig.3 Comparison of the dimensionless growth rate obtained from the present model with ref.[17] and analytical data[29]

圖4 液膜流動(dòng)計(jì)算域

式中y為計(jì)算域豎直方向各點(diǎn)縱坐標(biāo)值。氣相初速度給定為

(28)

圖5給出了在Reg=5592.377(表1中 Case 3)時(shí)液膜在預(yù)膜板表面鋪展形態(tài)的變化。隨著演化的進(jìn)行，可以看出在氣流的剪切作用下，液膜表面逐漸形成波浪形突起并向右推進(jìn)，直至離開計(jì)算區(qū)域。液膜表面形成波浪形突起是由于氣液兩相間的速度差造成的Kelvin-Helmholtz不穩(wěn)定性引起的，這與Baptiste等[30]對(duì)該速度比下的液膜流動(dòng)形態(tài)描述相符。

圖6分別描繪了表1中Case 1～3不同氣流速度下氣液相界面的變化?？梢婋S著氣體速度的增大，氣體驅(qū)動(dòng)液膜形態(tài)改變的進(jìn)程也在加快，流體界面波動(dòng)幅度增大。事實(shí)上，隨著入口氣體速度的增加，氣液間的剪切速度差也隨之增加，較高的速度差會(huì)導(dǎo)致更強(qiáng)的Kelvin-Helmholtz不穩(wěn)定性現(xiàn)象，從而在氣液界面掀起更高振幅的波浪。為了更好地比較氣流速度對(duì)氣液相界面波動(dòng)的影響，本文統(tǒng)計(jì)了計(jì)算域中液膜的厚度峰值A(chǔ)隨時(shí)間的變化。如圖7所示，縱坐標(biāo)A/Hl代表各時(shí)刻液膜厚度峰值與初始液膜厚度的比。隨著演化的開始，氣液間的剪切速度差引起的Kelvin-Helmholtz不穩(wěn)定性現(xiàn)象使相界面出現(xiàn)擾動(dòng)，液膜表面出現(xiàn)表面波，液膜厚度峰值有所提升，并隨著演化的進(jìn)行不斷波動(dòng)。當(dāng)氣體速度較小時(shí)，液膜在氣體剪切力作用下產(chǎn)生較小振幅的表面波并向前推進(jìn)直至流出計(jì)算域，待氣液兩相體系穩(wěn)定后氣液界面會(huì)出現(xiàn)周期性微小波動(dòng)，液膜厚度基本保持穩(wěn)定。隨著氣體速度的提升，液膜表面出現(xiàn)更為頻繁且大振幅的波動(dòng)。當(dāng)Reg=5592.377，演化進(jìn)行到一定階段后，可以明顯看到規(guī)律性的波浪出現(xiàn)直至離開，從而出現(xiàn)圖7液膜厚度峰值的周期性波動(dòng)。這也印證了本文所述，氣體速度的提高增大了氣液相界面的不穩(wěn)定性和波動(dòng)性。

圖5 Case 3中氣液界面隨時(shí)間的演化

圖6 Case 1～3在t=2.5 ms時(shí)氣液界面狀態(tài)

圖8描述了不同氣體速度下在t=25 ms時(shí)刻液膜沿預(yù)膜板流動(dòng)的速度，并與液膜入口速度Ul進(jìn)行比較?？梢钥闯觯耗ち鲃?dòng)速度在經(jīng)歷輕微增長(zhǎng)后沿著板面逐漸下降趨于穩(wěn)定。這是因?yàn)樵跉庖洪_始相互作用的入口段，液膜受到氣流剪切力的作用速度增高，而后期隨著液膜不斷向后鋪展,在流體粘滯力的影響下又逐漸降低。氣體速度很大程度影響液膜速度，液膜速度隨氣體驅(qū)動(dòng)速度的增大而增大。最后本文討論了氣體速度對(duì)平均液膜厚度的影響。平均液膜厚度是流體流動(dòng)的一個(gè)重要參數(shù)，二維流動(dòng)中的平均液膜厚度相當(dāng)于預(yù)膜板上的液體持液量，反應(yīng)液體流量的變化。本文平均液膜厚度定義為H=Q/L,L為液膜板長(zhǎng)度，初始時(shí)刻Qo=HlL。如圖9所示，與本文液膜表面狀態(tài)相對(duì)應(yīng)，氣體雷諾數(shù)的增高引起波液膜表面的波浪形變化，計(jì)算域內(nèi)的液膜平均厚度也隨之發(fā)生相應(yīng)的波動(dòng)。可以看出，隨著氣體速度的升高，液膜的平均厚度隨之降低[30]，這是因?yàn)樵诟咚贇饬髯饔孟?，相同的時(shí)間里液膜會(huì)以更快速度向后鋪展，減少了液膜在預(yù)膜板后段的堆積，因而固定長(zhǎng)度段的平均液膜厚度相應(yīng)降低，這與圖8氣液間剪切力的增大會(huì)促進(jìn)液膜的流動(dòng)結(jié)論相一致。

圖7 Case 1～3中液膜峰值隨時(shí)間的變化

圖8 Case 1～3中Reg對(duì)液膜流動(dòng)速度的影響

圖9 Case 1～3中Reg對(duì)平均液膜厚度的影響

5 結(jié) 論

本文對(duì)Liang等[17]提出的基于相場(chǎng)理論的LBM模型進(jìn)行修正，通過引入一個(gè)額外的界面力來消除兩相界面間由于密度變化而導(dǎo)致無法滿足體系不可壓縮條件的影響，并驗(yàn)證了其實(shí)用性。而后使用該模型對(duì)于不同氣體速度下液膜流動(dòng)的形態(tài)進(jìn)行模擬，并總結(jié)了相應(yīng)的規(guī)律。

(1) 經(jīng)過修正后的模型能更好地模擬較大密度比條件下的兩相流問題，且對(duì)于氣液相界面變化的捕捉更為準(zhǔn)確。

(2) 氣體雷諾數(shù)增高時(shí)，氣體驅(qū)動(dòng)液膜在預(yù)膜板上的流動(dòng)速度隨之增高，同時(shí)會(huì)引起氣液相界面間出現(xiàn)更大振幅的表面波。

(3) 高雷諾數(shù)氣體驅(qū)動(dòng)下液膜流動(dòng)速度較快，因而預(yù)膜板上液膜平均厚度降低，液膜堆積現(xiàn)象減弱。