表面張力對(duì)高雷諾數(shù)Rayleigh-Taylor不穩(wěn)定性后期增長(zhǎng)的影響*

黃皓偉 梁宏 徐江榮

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

采用多相流的相場(chǎng)格子Boltzmann方法數(shù)值研究了微通道內(nèi)高雷諾數(shù)單模Rayleigh-Taylor (RT)不穩(wěn)定性的后期演化規(guī)律， 重點(diǎn)分析表面張力對(duì)相界面動(dòng)力學(xué)行為以及氣泡與尖釘增長(zhǎng)的影響.數(shù)值實(shí)驗(yàn)表明，隨著界面張力的增大， 可以有效降低演化過(guò)程中相界面結(jié)構(gòu)的復(fù)雜程度， 并抑制不穩(wěn)定性后期相界面破裂形成離散液滴.另外， 增大表面張力可以先促進(jìn)后抑制氣泡振幅的增長(zhǎng)， 而當(dāng)表面張力較小時(shí)， 尖釘振幅增長(zhǎng)曲線之間并無(wú)明顯差別， 當(dāng)表面張力增大到一定值后， 它對(duì)尖釘振幅的抑制效果可明顯地被觀察到.進(jìn)一步， 根據(jù)不穩(wěn)定性速度增長(zhǎng)曲線， 將高雷諾數(shù)單模RT不穩(wěn)定性的演化劃分為線性增長(zhǎng)、飽和速度增長(zhǎng)、重加速、混沌混合四個(gè)發(fā)展階段.數(shù)值計(jì)算獲取氣泡與尖釘?shù)娘柡退俣确习缑鎻埩π?yīng)的勢(shì)流理論模型.另外還統(tǒng)計(jì)了不同表面張力和Atwood數(shù)下表征RT不穩(wěn)定性后期演化的氣泡與尖釘增長(zhǎng)率， 結(jié)果顯示氣泡與尖釘后期增長(zhǎng)率隨著表面張力的增大總體上呈現(xiàn)出先促進(jìn)后抑制的規(guī)律.最后， 從數(shù)值計(jì)算和理論分析兩方面研究了不同Atwood數(shù)下RT 不穩(wěn)定性發(fā)生的臨界表面張力， 發(fā)現(xiàn)兩者結(jié)果符合得很好， 并且臨界表面張力隨著流體Atwood數(shù)的增大而增大.

1 引 言

瑞利-泰勒(Rayleigh-Taylor， RT)不穩(wěn)定性現(xiàn)象是一種常見(jiàn)的流體不穩(wěn)定性問(wèn)題， 它發(fā)生于重力場(chǎng)作用下， 兩種不同密度的流體界面受到微小擾動(dòng)后.這類問(wèn)題不僅存在于自然現(xiàn)象中， 例如卷云的形成， 還在天體物理學(xué)、慣性約束聚變、地球氣候?qū)W等領(lǐng)域有著廣泛且重要的應(yīng)用， 也是多相流體界面、不穩(wěn)定性、湍流混合等基本問(wèn)題的理論基礎(chǔ)，因而吸引著許多學(xué)者對(duì)其開(kāi)展相關(guān)研究[1].RT不穩(wěn)定性的先驅(qū)工作歸功于英國(guó)流體力學(xué)大師 Rayleigh和Taylor[2，3]， 他們創(chuàng)造性地提出了線性穩(wěn)定性理論， 并證明了無(wú)黏流體的初始擾動(dòng)呈現(xiàn)指數(shù)增長(zhǎng)的規(guī)律.后來(lái)， Lewis[4]實(shí)驗(yàn)驗(yàn)證了線性增長(zhǎng)理論， 還報(bào)告了RT不穩(wěn)定性隨后進(jìn)入非線性發(fā)展階段.Sharp[5]對(duì)RT不穩(wěn)定性的早期研究工作做了很好的總結(jié)， 并定性地將RT不穩(wěn)定性的增長(zhǎng)劃分為四個(gè)發(fā)展階段.根據(jù)界面初始擾動(dòng)的不同， RT不穩(wěn)定性可以分為單模和多模問(wèn)題， 而本文主要關(guān)注單模RT不穩(wěn)定性的研究.2001年， Waddell 等[6]實(shí)驗(yàn)研究了二維單模 RT不穩(wěn)定性問(wèn)題， 指出擾動(dòng)振幅在初始階段以指數(shù)形式增長(zhǎng)， 即h=a1eγt+a2e-γt， 其中γ是線性增長(zhǎng)因子， 定義為其中g(shù)是重力加速度，At是Atwood數(shù)，k是波數(shù)，σ是流體間的表面張力，ρh和ρl分別為重輕流體的密度.進(jìn)一步， 他們還觀察到在不穩(wěn)定性演化后期氣泡與尖釘?shù)钠骄俣冉咏诔?shù).Glimm等[7]利用前追蹤方法模擬了二維單模RT不穩(wěn)定性問(wèn)題， 發(fā)現(xiàn)不穩(wěn)定性的增長(zhǎng)經(jīng)歷飽和速度階段后會(huì)出現(xiàn)重加速階段.后來(lái)， Wilkinson和Jacobs[8]對(duì)三維單模RT不穩(wěn)定性進(jìn)行了實(shí)驗(yàn)研究， 指出由于界面渦結(jié)構(gòu)間的相互作用， 氣泡與尖釘在后期階段被重加速， 導(dǎo)致演化速度高于Goncharov勢(shì)能模型[9]的理論值.2012 年， Ramaprabhu等[10]大規(guī)模模擬了不同Atwood數(shù)和雷諾數(shù)下三維單模RT不穩(wěn)定性的后期演化問(wèn)題， 首次將高雷諾數(shù)RT不穩(wěn)定性的發(fā)展歸結(jié)為線性增長(zhǎng)、飽和速度增長(zhǎng)、重加速、混沌混合四個(gè)發(fā)展階段， 并進(jìn)一步觀察到氣泡與尖釘在重加速階段的演化速度超過(guò)經(jīng)典勢(shì)流模型[9]的解析解， 而在后期的混沌混合階段， 發(fā)現(xiàn)氣泡與尖釘?shù)脑鲩L(zhǎng)會(huì)受到抑制.Wei和Livescu[11]利用直接數(shù)值模擬方法高精度研究了低Atwood數(shù)下二維單模RT不穩(wěn)定性問(wèn)題， 觀察到不穩(wěn)定性在高雷諾數(shù)條件下同樣經(jīng)歷了四個(gè)發(fā)展階段， 但不同于Ramaprabhu等[10]的結(jié)果， 他們發(fā)現(xiàn)氣泡的平均振幅在后期混沌階段具有二次增長(zhǎng)的規(guī)律.Lai等[12]采用離散 Boltzmann 方法研究了RT不穩(wěn)定性的熱動(dòng)非平衡效應(yīng).近期， Liang等[13-15]采用介觀格子Boltzmann (lattice Boltzmann， LB)方法模擬了二維及三維微通道內(nèi)單模RT不穩(wěn)定性后期演化問(wèn)題， 分析了Atwood數(shù)和雷諾數(shù)對(duì)不穩(wěn)定性的增長(zhǎng)階段和相界面動(dòng)力學(xué)行為的影響.

上述研究均忽略了表面張力對(duì)RT不穩(wěn)定性的影響， 已有研究表明RT不穩(wěn)定性在表面張力作用下可以顯示毛細(xì)波、收縮、破裂等獨(dú)特的界面動(dòng)力學(xué)行為[16].目前， 絕大多數(shù)的工作[16-20]集中于分析表面張力對(duì)初始擾動(dòng)為多模態(tài)的RT不穩(wěn)定性增長(zhǎng)， 而較少學(xué)者考慮表面張力對(duì)單模RT不穩(wěn)定性演化的影響.Daly[21]采用有限差分法數(shù)值研究了表面張力對(duì)單模RT不穩(wěn)定性早期增長(zhǎng)的影響， 獲取的不同界面張力下線性增長(zhǎng)因子與修正的線性穩(wěn)定性理論一致， 并且發(fā)現(xiàn)存在表面張力時(shí)可以抑制界面的卷吸行為.Zhang等[22]利用多相流LB方法模擬了單模RT不穩(wěn)定性問(wèn)題， 發(fā)現(xiàn)表面張力作用下可以減弱相界面的卷吸程度， 在演化后期誘導(dǎo)相界面發(fā)生夾斷生成離散液滴， 并且觀察到界面張力存在時(shí)可以減小氣泡與尖釘?shù)难莼俣?基于Goncharov勢(shì)流理論[9]， Young和Ham[23]分析界面張力對(duì)無(wú)黏流體的單模RT不穩(wěn)定性增長(zhǎng)的影響， 并提出了包含表面張力效應(yīng)的氣泡漸進(jìn)速度的解析模型，

進(jìn)一步， 他們還利用水平集數(shù)值方法模擬了表面張力作用下單模RT不穩(wěn)定性在飽和速度階段的增長(zhǎng)， 發(fā)現(xiàn)獲得的氣泡漸進(jìn)速度與提出的解析模型相一致.同樣從勢(shì)流理論出發(fā)， Sohn[24]分析了流體黏性和表面張力對(duì)單模RT不穩(wěn)定性發(fā)展的影響， 發(fā)現(xiàn)黏性和表面張力均抑制不穩(wěn)定性的發(fā)展， 減小氣泡的飽和增長(zhǎng)速度.結(jié)合黏性效應(yīng)， 氣泡在飽和速度階段的漸進(jìn)速度修正為

其中νh表示重流體的運(yùn)動(dòng)學(xué)黏性.夏同軍等[25]分別基于勢(shì)流模型和Zufiria模型[26]分析了表面張力對(duì)RT不穩(wěn)定性非線性增長(zhǎng)階段的影響， 推導(dǎo)出兩種模型下氣泡漸進(jìn)速度和漸進(jìn)曲率的解析表達(dá)式， 結(jié)果表明， 界面張力降低了氣泡的速度， 但對(duì)曲率沒(méi)有影響.另外發(fā)現(xiàn)， 基于勢(shì)流模型所得氣泡漸進(jìn)速度大于Zufiria模型的預(yù)測(cè)值.進(jìn)一步， 基于Zufiria理論模型， Li等[27]分析了流體黏性和表面張力對(duì)不穩(wěn)定性飽和速度階段的影響， 提出了包含兩種效應(yīng)的氣泡漸進(jìn)速度的解析式.Guo等[28]建立了含表面張力的單模RT不穩(wěn)定性的弱非線性模型， 分析了RT不穩(wěn)定性在表面張力作用下從線性增長(zhǎng)到非線性增長(zhǎng)的轉(zhuǎn)變.

綜上所述， 已有學(xué)者分析了表面張力對(duì)單模RT不穩(wěn)定性增長(zhǎng)的影響， 但這些研究聚焦于不穩(wěn)定性中氣泡增長(zhǎng)的中前期階段， 而對(duì)表面張力作用下不穩(wěn)定性的演化后期及尖釘增長(zhǎng)的描述還未有報(bào)道.另外， 上述研究所考慮的雷諾數(shù)均較小.本文將采用基于相場(chǎng)理論的多松弛LB方法研究冗長(zhǎng)微通道內(nèi)高雷諾數(shù)RT不穩(wěn)定性的演化規(guī)律， 著重分析表面張力對(duì)相界面動(dòng)力學(xué)行為及氣泡與尖釘后期增長(zhǎng)的影響.

2 數(shù)值模型

LB方法[29]是近幾十年發(fā)展起來(lái)的一種流體系統(tǒng)建模和模擬的介觀方法， 因其具有易于刻畫流體間與流固間的微觀相互作用、高性能并行計(jì)算及自動(dòng)追蹤相界面等優(yōu)點(diǎn)， 因而被廣泛用來(lái)研究多相流界面動(dòng)力學(xué)問(wèn)題.目前， 從流體間相互作用力的不同物理背景， 已提出了多種不同類別的多相流LB模型， 包括顏色模型、偽勢(shì)模型、自由能模型和相場(chǎng)LB模型.相比前三類模型， 相場(chǎng)LB模型因在相界面追蹤方面具有堅(jiān)實(shí)的物理學(xué)基礎(chǔ)而受到一些學(xué)者的廣泛關(guān)注， 感興趣的讀者可以參考近期關(guān)于多相流體的相場(chǎng)LB方法的綜述論文[30].相場(chǎng)模型中界面追蹤方程包括Cahn-Hilliard方程和Allen-Cahn方程， 基于Allen-Cahn模型的LB方法在求解序參數(shù)時(shí)具有更低的數(shù)值耗散， 適用于模擬大密度比兩相流問(wèn)題[31]， 而基于Cahn-Hilliard模型的LB方法在模擬界面大拓?fù)渥兓瘯r(shí)則更具有優(yōu)勢(shì).本文采用Liang等[13]提出的基于Cahn-Hilliard方程的相場(chǎng)LB模型， 該模型克服了已有模型中相界面求解精度低以及數(shù)值穩(wěn)定性差等問(wèn)題， 提高了LB方法在相界面捕獲方面的能力， 可顯式計(jì)算流場(chǎng)中的宏觀量， 并在求解復(fù)雜的流體界面RT不穩(wěn)定性問(wèn)題上顯示出較大的潛力[14，32].在相場(chǎng)理論中， 多相流體系統(tǒng)的總自由能泛函可表述為[33]

其中φ是序參數(shù)， 用于捕捉相界面;F(φ) 是體相區(qū)的自由能密度泛函， 具有雙勢(shì)阱形式F(φ)=βφ2(1-φ)2， 參數(shù)α和β依賴于表面張力和界面厚度對(duì)總的自由能泛函求變分運(yùn)算， 可以獲得化學(xué)勢(shì)的表達(dá)式為

進(jìn)一步， 假設(shè)序參數(shù)φ由流體速度u發(fā)生對(duì)流， 和化學(xué)勢(shì)μ產(chǎn)生擴(kuò)散， 從而得到描述相界面動(dòng)力學(xué)的經(jīng)典Cahn-Hilliard相場(chǎng)方程:

其中，M是遷移率.為了描述多相流體輸運(yùn)過(guò)程，相界面追蹤的Cahn-Hilliard方程需要耦合含外力的不可壓Navier-Stokes方程:

其中ρ是流體密度，p是流體動(dòng)力學(xué)壓力，ν是運(yùn)動(dòng)學(xué)黏性， 界面張力Fs取勢(shì)形式為Fs=μ?φ，G是系統(tǒng)中施加的外力.

本文的LB模型利用兩個(gè)分布函數(shù)fi和gi分別來(lái)求解Cahn-Hilliard和Navier-Stokes相互耦合的多相系統(tǒng)控制方程， 對(duì)應(yīng)的多松弛碰撞算子的演化方程可表示為[29]

其中，fi(x，t) ，gi(x，t) 是粒子分布函數(shù)，和是分布函數(shù)對(duì)應(yīng)的平衡態(tài)分布函數(shù)，ci是離散速度，M表示碰撞矩陣，Sf和Sg是對(duì)角松弛矩陣，F(xiàn)i(x，t) 和Gi(x，t) 分別為源項(xiàng)和外力項(xiàng)的分布函數(shù).為了恢復(fù)正確的宏觀方程， 平衡態(tài)分布函數(shù)可設(shè)定為[13]

且

其中，η是調(diào)節(jié)遷移率的自由參數(shù)，cs是聲速，ωi是權(quán)系數(shù).對(duì)于二維流動(dòng)， 本文采用流行D2Q9格子 模 型， 其 對(duì) 應(yīng) 的 權(quán) 系 數(shù)ωi為ω0=4/9 ，ω1-4=1/9 ，ω5-8=1/36 ， 離散速度ci定義為[34，35]

其中，c=δx/δt，δx和δt分別代表格子長(zhǎng)度和時(shí)間步長(zhǎng)，對(duì)于多相流LB模型， 通常設(shè)定c=δx=δt=1.另外， D2Q9 模型所對(duì)應(yīng)的碰撞矩陣M定義為[36]

以及松弛矩陣Sf和Sg可表示為

另外， 流體密度ρ可看成序參數(shù)φ的線性函數(shù):

通過(guò)Chapman-Enskog多尺度分析[13]， 可以證明本文采用的相場(chǎng)LB模型能夠正確地恢復(fù)到Cahn-Hilliard方程和Navier-Stokes方程， 并且遷移率M和運(yùn)動(dòng)學(xué)黏度ν的數(shù)學(xué)表達(dá)式為

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

本文關(guān)注冗長(zhǎng)的微管道內(nèi)多相流體界面單模擾動(dòng)的RT不穩(wěn)定性， 網(wǎng)格劃分為L(zhǎng)×W=8160×544.初始時(shí)刻， 在微管道的中心截面處施加波長(zhǎng)為W的微小擾動(dòng):

其中， 波數(shù)k定義為k=2π/W， 并設(shè)定序參量的初始分布為

其中， 界面厚度D固定為 4 個(gè)格子網(wǎng)格.需要指出的是， RT不穩(wěn)定性的界面演化可以用兩個(gè)無(wú)量綱參數(shù)來(lái)描述， 即雷諾數(shù)(Reynolds number，Re)和阿特伍德數(shù)(Atwood number，At)， 分別定義為:

本文著重分析高雷諾數(shù)RT不穩(wěn)定性的后期演化規(guī)律， 因而在數(shù)值模擬中雷諾數(shù)均設(shè)定為Re=104.為了表征重力效應(yīng)， 對(duì)微管道中輕重流體均施加一個(gè)豎直方向上的浮力:

其中， 重流體密度設(shè)定為ρh=1 ， 輕流體的密度可由給定的阿特伍德數(shù)確定， 重力加速度滿足關(guān)系式另 外， Peclet數(shù) (Pe)設(shè) 定 為50，其定義式為邊界條件設(shè)置如下: 左右邊界采用周期性邊界條件， 而上下邊界均采用無(wú)滑移半反彈邊界格式.半反彈格式將壁面設(shè)定在格點(diǎn)中線上， 假定粒子與壁面碰撞后速度發(fā)生逆轉(zhuǎn)， 即臨近壁面流體點(diǎn)xf的粒子分布函數(shù)設(shè)定為其中ci=-cˉi是指向壁面的粒子速度，表示碰撞后的粒子分布函數(shù)[39].具體而言， 當(dāng)xf位于上壁面時(shí)， 半反彈邊界格式實(shí)現(xiàn)的數(shù)學(xué)表達(dá)式為

而當(dāng)xf為下壁面時(shí)， 半反彈邊界條件可設(shè)置為

在本文的研究中， 選取特征長(zhǎng)度為擾動(dòng)的初始波長(zhǎng)W， 特征速度和特征時(shí)間分別給定為和下文統(tǒng)計(jì)的相關(guān)物理量均已被相應(yīng)特征值所無(wú)量綱化.

圖1給出了At=0.7 時(shí)， 四種不同表面張力下高雷諾數(shù)RT不穩(wěn)定性的相界面演化過(guò)程.可以發(fā)現(xiàn)， 對(duì)于不同的表面張力， 相界面在初始階段表現(xiàn)出相似的動(dòng)力學(xué)行為: 重流體和輕流體相互滲入形成了尖釘和氣泡結(jié)構(gòu).隨后在不同表面張力下， 尖釘和氣泡的增長(zhǎng)表現(xiàn)出明顯的差異.在表面張力很小， 僅為 5×10-6時(shí)， 尖釘持續(xù)下降并出現(xiàn)非線性效應(yīng)， 這是由于氣泡和尖釘界面上存在的速度差誘導(dǎo)產(chǎn)生了非線性 Kelvin-Helmholtz (KH)不穩(wěn)定性， 進(jìn)而使尖釘前端卷起形成一對(duì)渦.隨著兩個(gè)旋渦的持續(xù)增長(zhǎng)， 在卷起的尾端出現(xiàn)了二級(jí)渦， 并且混合區(qū)域內(nèi)的非線性效應(yīng)也越來(lái)越劇烈， 導(dǎo)致界面在多個(gè)位置發(fā)生卷起， 形成了復(fù)雜的界面拓?fù)浣Y(jié)構(gòu).在不穩(wěn)定性繼續(xù)發(fā)展過(guò)程中也伴隨著流體界面混沌的破裂， 最終導(dǎo)致混合區(qū)域存在著大量的離散液滴， 同時(shí)也觀察到相界面的增長(zhǎng)失去了左右對(duì)稱性.相界面的非對(duì)稱性發(fā)展同樣在高雷諾數(shù)下混相RT不穩(wěn)定性現(xiàn)象中被觀察到[40].隨著表面張力的增加， KH不穩(wěn)定性的出現(xiàn)開(kāi)始延后， 一級(jí)渦的產(chǎn)生相應(yīng)地延緩， 隨后產(chǎn)生的渦的數(shù)量也相應(yīng)地減少， 后期的結(jié)構(gòu)復(fù)雜度下降.不同于σ=5×10-6情形， RT不穩(wěn)定性在較大界面張力的增長(zhǎng)過(guò)程中，相界面圖案始終保持著左右對(duì)稱性， 整個(gè)演化過(guò)程中界面破裂的現(xiàn)象減少， 系統(tǒng)中離散的液滴數(shù)量也相應(yīng)減少， 尤其當(dāng)σ=1×10-2時(shí)， 不穩(wěn)定性的發(fā)展過(guò)程中未觀察到離散液滴.另外， 我們通過(guò)數(shù)值實(shí)驗(yàn)進(jìn)一步發(fā)現(xiàn)， 繼續(xù)增大表面張力， 擾動(dòng)會(huì)隨時(shí)間發(fā)生振蕩并逐漸恢復(fù)平穩(wěn)， RT不穩(wěn)定現(xiàn)象并未發(fā)生， 這表明存在一個(gè)臨界的表面張力的值， 超過(guò)臨界值后， RT不穩(wěn)定現(xiàn)象不會(huì)發(fā)生， 而小于該臨界值， RT 不穩(wěn)定性將會(huì)發(fā)生.下文將細(xì)致討論臨界表面張力的影響因素以及關(guān)系表達(dá)式.

氣泡和尖釘?shù)恼穹菃文T不穩(wěn)定研究中兩個(gè)重要的物理量， 因而， 本文統(tǒng)計(jì)了不同表面張力條件下隨時(shí)間變化的氣泡與尖釘振幅， 以定量比較界面張力的影響.圖2為氣泡與尖釘振幅在不同表面張力下隨時(shí)間的演化曲線.可以看出， 對(duì)于氣泡而言， 隨著表面張力的增加， 氣泡振幅的增長(zhǎng)并不是一直減緩的， 在表面張力很小的一段范圍內(nèi)，表面張力反而促進(jìn)了氣泡的增長(zhǎng).當(dāng)σ=5×10-4時(shí)， 氣泡振幅的增長(zhǎng)在發(fā)展前期受表面張力影響不大， 但是 在 后 期， 增 長(zhǎng)明 顯快 于σ=5×10-6和σ=5×10-5時(shí)的氣泡振幅.而表面張力增加到σ=5×10-3后， 氣泡振幅的增長(zhǎng)無(wú)論是前期還是后期均慢于σ=5×10-4時(shí)的增長(zhǎng).繼續(xù)增大表面張力到σ=1×10-2， 在前中期的抑制效果更加明顯.上述結(jié)果與文獻(xiàn)[23，24]中普遍認(rèn)為表面張力抑制氣泡的增長(zhǎng)有所不同， 原因可能是前人所考慮的界面張力范圍有限.而對(duì)于尖釘而言， 當(dāng)表面張力相對(duì)較小時(shí)， 振幅曲線沒(méi)有明顯的差異， 這表明當(dāng)表面張力足夠小時(shí)， 其對(duì)尖釘振幅的影響不再重要.但是當(dāng)表面張力增大到一定值后， 它對(duì)尖釘振幅的抑制效果則可以被明顯觀察到.我們還進(jìn)一步模擬了不同表面張力條件下中低 Atwood數(shù)(At=0.5，0.1)的RT不穩(wěn)定性問(wèn)題， 結(jié)果表明， 隨著界面張力的增大， 表面張力對(duì)于氣泡振幅的影響同樣是先促進(jìn)后抑制的， 說(shuō)明這是一個(gè)普遍的規(guī)律， 而表面張力可以減緩尖釘?shù)脑鲩L(zhǎng).

圖1 R e=104 ， A t=0.7 時(shí)， 表 面 張 力 對(duì) 高 雷 諾 數(shù) 非 混 相RT不 穩(wěn) 定 相 界 面 演 化 的 影 響 (a) σ =5×10-6 ; (b)σ=5×10-4 ; (c) σ =5×10-3 ; (d)σ=1×10-2Fig.1.Effect of the surface tension on the evolution of phase interface in the immiscible RT instability with R e=104 ， A t=0.7 :(a) σ =5×10-6 ; (b) σ =5×10-4 ; (c) σ =5×10-3 ; (d) σ =1×10-2.

圖2 表面張力對(duì)隨時(shí)間變化的氣泡和尖釘振幅的影響Fig.2.Influence of surface tension on the time variations of bubble and spike amplitudes.

接著對(duì)氣泡和尖釘?shù)脑鲩L(zhǎng)速度進(jìn)行定量分析.圖3給出了在不同表面張力下氣泡與尖釘速度隨時(shí)間變化的演化曲線.根據(jù)文獻(xiàn)[11，14]， 高雷諾數(shù)下的單模RT不穩(wěn)定性的發(fā)展經(jīng)歷了四個(gè)不同的階段， 包括線性增長(zhǎng)階段、飽和速度階段、重加速階段和混沌階段.當(dāng)表面張力足夠小時(shí)， 氣泡和尖釘速度在線性階段的發(fā)展十分相似， 增長(zhǎng)曲線基本重合， 即表面張力的效應(yīng)不顯著.而當(dāng)表面張力增大到一定值后， 線性階段的氣泡與尖釘速度減緩明顯.上述數(shù)值結(jié)果與線性穩(wěn)定性理論[6]分析一致:擾動(dòng)振幅在初始階段的增長(zhǎng)具有指數(shù)形式， 其線性增長(zhǎng)因子如(1)式所示， 這從理論上表明當(dāng)表面張力非常小時(shí)， 其值對(duì)初始階段擾動(dòng)增長(zhǎng)的影響可以忽略， 而當(dāng)表面張力較大時(shí)， 表面張力的效應(yīng)逐漸顯著， 增大其值可以有效地減弱擾動(dòng)的線性發(fā)展.緊接著線性增長(zhǎng)階段， 從圖3可以發(fā)現(xiàn)氣泡與尖釘?shù)脑鲩L(zhǎng)將以近似恒定的速度增長(zhǎng)， 這表明不穩(wěn)定性的發(fā)展進(jìn)入飽和速度階段.Goncharov[9]基于經(jīng)典的勢(shì)流理論模型預(yù)測(cè)了無(wú)黏流體在忽略表面張力條件的單模RT不穩(wěn)定性中氣泡與尖釘?shù)娘柡退俣龋?其無(wú)量綱形式分別表示為

Sohn[24]分析流體黏性和表面張力對(duì)不穩(wěn)定性增長(zhǎng)的影響， 并基于Goncharov的解析勢(shì)流模型， 提出了修正的氣泡漸進(jìn)速度的數(shù)學(xué)表達(dá)式((3)式).受Goncharov[9]工作的啟發(fā)， 耦合黏性和表面張力效應(yīng)的尖釘恒定速度的無(wú)量綱形式可以寫為

其中Bo為邦德數(shù)， 定義為Bo=ρhgW2/σ.根據(jù)Sohn[24]的解析模型， 我們發(fā)現(xiàn)氣泡與尖釘滿足σ≤5×10-3時(shí)漸進(jìn)速度的理論解幾乎差別不大，因而在圖3中僅給出σ=5×10-3和σ=1×10-2條件下氣泡與尖釘飽和速度的理論值.可以看出，氣泡漸進(jìn)速度的數(shù)值計(jì)算結(jié)果與理論解符合較好，并且發(fā)現(xiàn)在界面張力較大的情形下， 增大其值可以有效地降低氣泡的漸進(jìn)速度.而對(duì)于尖釘而言， 界面張力對(duì)尖釘飽和速度的解析解影響不大， 數(shù)值計(jì)算的尖釘漸進(jìn)速度在表面張力較小時(shí)與理論解一致， 而在表面張力較大時(shí)， 數(shù)值預(yù)測(cè)的結(jié)果略低于解析解， 這表明界面張力較大時(shí)可以輕微地減弱尖釘在飽和速度階段的增長(zhǎng).隨著不穩(wěn)定性的非線性強(qiáng)度持續(xù)增強(qiáng)， 氣泡與尖釘?shù)难莼俣葧?huì)超過(guò)各自對(duì)應(yīng)的飽和速度理論值， 這說(shuō)明不穩(wěn)定性的發(fā)展進(jìn)入了重加速階段.從圖3可以觀察到， 氣泡的再加速階段在小表面張力(σ=5×1 0-6)時(shí)并不明顯，在增大表面張力至σ=1×10-2后， 重加速階段則十分清晰.最后， 氣泡與尖釘速度隨時(shí)間上下波動(dòng)，顯示不穩(wěn)定的行為， 其值會(huì)被反復(fù)加速和減速，RT不穩(wěn)定性的演化進(jìn)入混沌混合階段.從圖3可以看出， 增大表面張力可以有效地抑制氣泡與尖釘速度在混沌混合階段的波動(dòng)程度.另外， 已有研究表明， 高雷諾數(shù)的單模RT不穩(wěn)定性的增長(zhǎng)在后期的混沌混合階段呈現(xiàn)出平均二次增長(zhǎng)的規(guī)律[11，14]，即氣泡與尖釘演化振幅隨時(shí)間變化的函數(shù)可表示為hs，b=αs，bAtgt2， 其中αs，b分別表示氣泡與尖釘?shù)暮笃谠鲩L(zhǎng)率系數(shù).值得注意的是， 氣泡與尖釘增長(zhǎng)率是RT不穩(wěn)定性研究中最為關(guān)心且重要的參數(shù)， 目前已提出了多種不同方法用于測(cè)量后期階段的氣泡與尖釘增長(zhǎng)率.為了顯示良好的收斂性， 本文采用Cabot等[41]提出的測(cè)量氣泡和尖釘?shù)暮笃谠鲩L(zhǎng)率的方法，表示氣泡與尖釘振幅對(duì)時(shí)間的一階導(dǎo)數(shù).圖4給出了不同Atwood數(shù)下氣泡與尖釘?shù)暮笃谠鲩L(zhǎng)率與表面張力的關(guān)系曲線.可以看出， 相同的Atwood數(shù)下， 尖釘?shù)暮笃谠鲩L(zhǎng)率遠(yuǎn)大于對(duì)應(yīng)的氣泡后期增長(zhǎng)率， 其效應(yīng)隨著Atwood數(shù)的增長(zhǎng)而越發(fā)顯著.而對(duì)于相同的界面張力， 尖釘后期增長(zhǎng)率會(huì)隨著Atwood數(shù)的增大而顯著遞增.另外還發(fā)現(xiàn)， 尖釘與氣泡的增長(zhǎng)系數(shù)隨著表面張力的增大總體上呈現(xiàn)出先遞增而后減小的趨勢(shì).

圖3 表面張力對(duì)隨時(shí)間演化的氣泡和尖釘增長(zhǎng)速度的影響.藍(lán)色和虛線分別表示勢(shì)能模型預(yù)測(cè)氣泡與尖釘速度在 σ =5×10-3 和 σ =1×10-2 時(shí)的解析解Fig.3.Influence of surface tension on the time evolutions of bubble and spike growth velocities.The blue and yellow dotted lines represent the analytical solutions of the bubble and spike velocities from potential flow model at σ=5×10-3 and σ =1×10-2.

圖4 不同Atwood數(shù)下氣泡與尖釘?shù)暮笃谠鲩L(zhǎng)率系數(shù)αb，s 隨表面張力的變化Fig.4.Variations of bubble and spike late-time growth coefficients α b，s with respect to the surface tension under different Atwood numbers.

最后討論RT不穩(wěn)定性現(xiàn)象發(fā)生的臨界表面張力.通過(guò)大規(guī)模的數(shù)值實(shí)驗(yàn)， 發(fā)現(xiàn)當(dāng)表面張力小于某個(gè)臨界閾值時(shí)， RT不穩(wěn)定性現(xiàn)象將會(huì)發(fā)生，而表面張力大于該臨界值時(shí)， 相界面是穩(wěn)定的，RT不穩(wěn)定性受到抑制.另外還發(fā)現(xiàn)， 臨界表面張力的值依賴于流體Atwood數(shù).臨界表面張力的范圍可以通過(guò)半分法來(lái)確定.對(duì)于固定流體Atwood數(shù)， 考慮一個(gè)較大界面張力的范圍 [σ0σ1] ，其中σ=σ0時(shí)， RT不穩(wěn)定性能發(fā)生， 而當(dāng)σ=σ1時(shí)， RT不穩(wěn)定性不能發(fā)生.接著， 考慮σ=(σ0+σ1)/2的情形， 如果通過(guò)數(shù)值模擬發(fā)現(xiàn) RT不穩(wěn)定性能發(fā)生， 則臨界表面張力坐落在 [ (σ0+σ1)/2σ1] 范圍內(nèi)， 否則臨界表面張力的范圍為 [σ0(σ0+σ1)/2].反復(fù)重復(fù)上述的數(shù)值實(shí)驗(yàn)， 最終可以確定臨界表面張力的數(shù)值預(yù)測(cè)的變化范圍.另外， 根據(jù)線性穩(wěn)定性理論， RT不穩(wěn)定性在初始階段的線性增長(zhǎng)因子如(1)式所示.要使RT不穩(wěn)定性能夠發(fā)生， 必須保證線性增長(zhǎng)因子中開(kāi)根號(hào)的值不小于 0， 因而可以推導(dǎo)出臨界表面張力的理論關(guān)系式為

這與利用文獻(xiàn)[24，27]提出的模型推導(dǎo)出的結(jié)果是相同的.圖5給出了LB方法數(shù)值計(jì)算獲得的不同Atwood數(shù)下的臨界表面張力以及理論模型所預(yù)測(cè)的臨界值.從圖5可以發(fā)現(xiàn)， 臨界表面張力會(huì)隨著Atwood數(shù)的增加而增大， 并且理論給出臨界值基本落在數(shù)值模擬所預(yù)測(cè)的范圍內(nèi)， 即數(shù)值計(jì)算結(jié)果與理論模型結(jié)果一致， 這也表明利用(30)式對(duì)臨界表面張力進(jìn)行計(jì)算是合理的.

圖5 不同 Atwood 數(shù)下的臨界表面張力Fig.5.Critical surface tensions at various Atwood numbers.

4 結(jié) 論

本文基于相場(chǎng)理論的LB方法模擬了冗長(zhǎng)微通道內(nèi)非混相流體的單模RT不穩(wěn)定性問(wèn)題， 考察了高雷諾數(shù)條件下表面張力對(duì)相界面動(dòng)態(tài)過(guò)程及擾動(dòng)增長(zhǎng)的影響.數(shù)值結(jié)果表明， Atwood數(shù)為0.7的單模RT不穩(wěn)定性在不同的表面張力下表現(xiàn)出明顯不同的界面動(dòng)力學(xué)特征: 低表面張力時(shí)， 不穩(wěn)定性演化后期流體界面出現(xiàn)混沌的破裂， 誘導(dǎo)系統(tǒng)中產(chǎn)生大量的離散液滴， 同時(shí)也觀察到相界面在后期演化過(guò)程中顯示非對(duì)稱性特征; 而隨著表面張力逐漸增大， 界面區(qū)域中產(chǎn)生離散液滴的數(shù)量減少， 在極大的表面張力作用下相界面甚至未出現(xiàn)破裂情形， 在整個(gè)演化過(guò)程中界面始終保持關(guān)于中軸線對(duì)稱.另外還分析了表面張力對(duì)氣泡與尖釘振幅增長(zhǎng)和演化速度的影響.在表面張力較小時(shí)， 表面張力對(duì)尖釘增長(zhǎng)的作用不顯著， 繼續(xù)增大表面張力可以明顯地抑制尖釘?shù)脑鲩L(zhǎng)， 而表面張力對(duì)氣泡增長(zhǎng)的影響則呈現(xiàn)出先促進(jìn)后抑制的規(guī)律.進(jìn)一步統(tǒng)計(jì)了不同表面張力和Atwood數(shù)下氣泡與尖釘演化后期的二次增長(zhǎng)系數(shù)， 可以發(fā)現(xiàn)對(duì)于固定的Atwood數(shù)， 尖釘?shù)暮笃谠鲩L(zhǎng)率遠(yuǎn)大于相應(yīng)的氣泡增長(zhǎng)率， 并且其效應(yīng)隨著Atwood數(shù)的增大而更加顯著.此外還發(fā)現(xiàn)， 氣泡與尖釘?shù)暮笃谠鲩L(zhǎng)率隨著表面張力的增大總體上呈現(xiàn)先增大而后遞減的規(guī)律.最后， 通過(guò)理論分析建立了RT不穩(wěn)定性現(xiàn)象發(fā)生的臨界表面張力的數(shù)學(xué)表達(dá)式， 揭示了臨界表面張力與流體Atwood數(shù)呈現(xiàn)遞增關(guān)系， 并且通過(guò)LB方法計(jì)算獲取的臨界表面張力與理論模型結(jié)果相符合.