于蔚然，周 訓，李維仲

(大連理工大學 海洋能源利用與節(jié)能教育部重點實驗室，大連 116024)

1 引 言

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

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

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

2 數(shù)學方法

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

2.1.1 控制方程

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

(1)

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

(2)

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

(3)

對于不可壓縮流體，帶外力項[20]的NS方程可表示為

·u=0

(4a)

(4b)

式中F為外力項總和，詳見下文。

2.1.2 AC方程的LBM演化方程

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

(5)

(6)

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

(7)

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

(8)

(9)

(10)

(11)

2.1.3 NS方程的LBM演化方程

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

δtGi(x,t)

(12)

(13)

(14)

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

(15)

式中F為總的外力項。

F=Fs+Fa+G

(16)

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

(17)

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

(18)

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

(19)

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

(20)

(21)

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

(22)

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

(23)

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

3 模型驗證

(24)

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

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

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

a=a0eα t

(25)

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

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

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

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

ux(y)=Ultanh(y/δ)

(26)

(27)

圖3 不同k*下增長率α*的計算結(jié)果與文獻[17]和

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

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

圖4 液膜流動計算域

式中y為計算域豎直方向各點縱坐標值。氣相初速度給定為

(28)

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

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

圖5 Case 3中氣液界面隨時間的演化

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

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

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

圖8 Case 1～3中Reg對液膜流動速度的影響

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

5 結(jié) 論

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

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

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

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