胡曉亮 梁宏? 王會(huì)利

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

2) (武漢紡織大學(xué)數(shù)學(xué)與計(jì)算機(jī)學(xué)院, 武漢 430200)

本文采用相場(chǎng)格子Boltzmann方法研究了豎直微通道內(nèi)中等Atwoods數(shù)流體的單模Rayleigh-Taylor不穩(wěn)定性問(wèn)題, 系統(tǒng)分析了雷諾數(shù)對(duì)相界面動(dòng)力學(xué)行為以及擾動(dòng)在各發(fā)展階段演化規(guī)律的影響.數(shù)值結(jié)果表明高雷諾數(shù)條件下, 不穩(wěn)定性界面擾動(dòng)的增長(zhǎng)經(jīng)歷了四個(gè)不同的發(fā)展階段, 包括線性增長(zhǎng)階段、飽和速度階段、重加速階段及混沌混合階段.在線性增長(zhǎng)階段, 我們計(jì)算獲得的氣泡與尖釘振幅符合線性穩(wěn)定性理論, 并且線性增長(zhǎng)率隨著雷諾數(shù)的增加而增大.在第二個(gè)階段, 我們觀察到氣泡與尖釘將以恒定的速度增長(zhǎng), 獲得的尖釘飽和速度略高于Goncharov經(jīng)典勢(shì)能模型的解析解[Phys.Rev.Lett.2002 88 134502], 這歸因于系統(tǒng)中產(chǎn)生了多個(gè)尺度的旋渦, 而渦之間的相互作用促進(jìn)了尖釘?shù)脑鲩L(zhǎng).隨著橫向速度和縱向速度的差異擴(kuò)大,氣泡和尖釘界面演化誘導(dǎo)產(chǎn)生的Kelvin–Helmholtz不穩(wěn)定性逐漸增強(qiáng), 從而流體混合區(qū)域出現(xiàn)許多不同層次的渦結(jié)構(gòu), 加速了氣泡與尖釘振幅的演化速度, 并在演化后期階段, 導(dǎo)致界面發(fā)生多層次卷起、劇烈變形、混沌破裂等行為, 最終形成了非常復(fù)雜的拓?fù)浣Y(jié)構(gòu).此外, 我們還統(tǒng)計(jì)了演化后期氣泡與尖釘?shù)臒o(wú)量綱加速度, 發(fā)現(xiàn)氣泡和尖釘?shù)恼穹诤笃诔尸F(xiàn)二次增長(zhǎng)規(guī)律, 其增長(zhǎng)率系數(shù)分別為0.045與0.233.而在低雷諾條件下, 重流體在不穩(wěn)定性后期以尖釘?shù)男问较蛳逻\(yùn)動(dòng)而輕流體以氣泡的形式向上升起.在整個(gè)演化過(guò)程中, 界面變得足夠光滑, 氣泡與尖釘在后期的演化速度接近于常數(shù), 未觀察到后期的重加速與混沌混合階段.

1 引 言

瑞利-泰勒 (Rayleigh-Taylor, RT)不穩(wěn)定性是一種經(jīng)典而又古老的流體界面不穩(wěn)定性現(xiàn)象.當(dāng)重流體置于輕流體之上并在相界面處施加一個(gè)微小的擾動(dòng), 兩相流體的界面是不穩(wěn)定的, RT不穩(wěn)定性現(xiàn)象將會(huì)發(fā)生.RT不穩(wěn)定性問(wèn)題廣泛存在于天體物理(卷云、超星系爆炸、蟹狀星云等)[1]、地球物理(鹽丘、油氣礦的形成)[2]、工程界(混合過(guò)程、慣性約束聚變)[3], 也在核能約束聚變熱能利用中扮演著重要角色.另外, 我國(guó)著名核物理學(xué)家賀賢土院士指出有效抑制RT不穩(wěn)定性的后期湍流混合是慣性約束聚變過(guò)程中實(shí)現(xiàn)點(diǎn)火成功的關(guān)鍵[4].因此, RT不穩(wěn)定性問(wèn)題在流體力學(xué)領(lǐng)域具有重要的理論價(jià)值和廣泛的工程應(yīng)用, 從而長(zhǎng)期以來(lái)吸引許多學(xué)者對(duì)其開(kāi)展了相關(guān)研究.早在1883年, 著名學(xué)者Rayleigh[5]在研究云的分層問(wèn)題中第一次描述了RT不穩(wěn)定性現(xiàn)象.后來(lái), Taylor[6]在原子能武器的研究中進(jìn)一步發(fā)現(xiàn)了RT不穩(wěn)定性現(xiàn)象.他們對(duì)RT不穩(wěn)定性進(jìn)行理論分析并提出了描述擾動(dòng)發(fā)展的經(jīng)典線性增長(zhǎng)理論: 當(dāng)界面初始擾動(dòng)的振幅小于它的波長(zhǎng)時(shí), 擾動(dòng)將變得不穩(wěn)定, 其振幅將以指數(shù)形式增長(zhǎng).后來(lái), 一些學(xué)者相繼分析了流體壓縮性、流體黏性、界面張力對(duì)擾動(dòng)增長(zhǎng)率的影響,并提出了修正的線性增長(zhǎng)理論[7–9].Lewis[10]第一次通過(guò)實(shí)驗(yàn)研究了RT不穩(wěn)定性問(wèn)題, 驗(yàn)證了Taylor的線性增長(zhǎng)理論, 并初步而定性地描述了不穩(wěn)定性中擾動(dòng)的發(fā)展階段.后來(lái), Sharp[11]進(jìn)一步研究了RT不穩(wěn)定性問(wèn)題, 并較為完善地劃分了不穩(wěn)定性發(fā)展的四個(gè)階段: 初始階段, 擾動(dòng)振幅符合經(jīng)典的線性增長(zhǎng)理論, 將以指數(shù)形式增長(zhǎng), 直至增長(zhǎng)到初始擾動(dòng)波長(zhǎng)的 1 0?40% ; 緊接著, 界面擾動(dòng)增長(zhǎng)轉(zhuǎn)變?yōu)榉蔷€性發(fā)展階段, 表現(xiàn)為重流體滲透到輕流體中呈現(xiàn)出蘑菇的形狀, 稱(chēng)為尖釘, 而輕流體上浮至重流體中, 稱(chēng)為氣泡; 在輕重流體相互滲透的過(guò)程中, 兩相系統(tǒng)的非線性強(qiáng)度逐漸增強(qiáng), 蘑菇頭部發(fā)生擠壓而尾部出現(xiàn)卷吸現(xiàn)象, 這表明產(chǎn)生了第二類(lèi)不穩(wěn)定性現(xiàn)象, 即開(kāi)爾文-赫姆霍茲(Kelvin-Helmholtz, KH)不穩(wěn)定性; 最后, 受兩類(lèi)不穩(wěn)定性的雙重影響, 氣泡和尖釘發(fā)生破裂行為,流體界面表現(xiàn)為后期的混沌狀態(tài).

自Sharp[11]的研究工作以來(lái), 許多學(xué)者對(duì)RT不穩(wěn)定性的各個(gè)發(fā)展階段開(kāi)展了研究, 包括單模態(tài)RT不穩(wěn)定性和隨機(jī)擾動(dòng)模態(tài)的RT不穩(wěn)定性[12–15].本文主要關(guān)注相界面初始擾動(dòng)為單模態(tài)的RT不穩(wěn)定性問(wèn)題.Waddell等[16]實(shí)驗(yàn)研究了低Atwood數(shù)下二維單模RT不穩(wěn)定性, 其結(jié)果表明擾動(dòng)振幅在初始階段以指數(shù)的形式增長(zhǎng), 獲得的線性增長(zhǎng)率與經(jīng)典的線性理論分析相一致, 還進(jìn)一步發(fā)現(xiàn)尖釘和氣泡的平均速度在演化后期接近于常數(shù).Goncharov[17]通過(guò)理論分析證實(shí)了這一點(diǎn),并定量地給出理想流體的尖釘和氣泡恒定的飽和速度

其中 ub是氣泡的速度, us是尖釘?shù)乃俣? g 是重力加速度, At是 Atwood 數(shù), k 是波數(shù), Cg對(duì)于三維情形取1, 而對(duì)于二維情形則取為3.Wilkinson和Jacobs[18]對(duì)三維單模RT不穩(wěn)定性問(wèn)題進(jìn)行了實(shí)驗(yàn)研究, 觀察到氣泡與尖釘在演化后期的增長(zhǎng)速度均高于Goncharov[17]的勢(shì)能理論模型的解析解,并將這種與理論預(yù)測(cè)不一致現(xiàn)象歸因于渦結(jié)構(gòu)間的相互作用.后來(lái), 一些學(xué)者還分析了渦結(jié)構(gòu)相互作用、流體黏性、表面張力對(duì)第二階段中氣泡演化速度的影響, 并提出了包含這些物理因素的理論表達(dá)式.

單模態(tài)RT不穩(wěn)定性的上述研究還僅限于擾動(dòng)發(fā)展的前三個(gè)階段, 而對(duì)于不穩(wěn)定性的后期階段, 非線性效應(yīng)往往越來(lái)越強(qiáng)烈, 流體界面則會(huì)顯示出非常復(fù)雜而劇烈的拓?fù)渥兓? 因而難以通過(guò)實(shí)驗(yàn)方法測(cè)量演化后期的相關(guān)物理統(tǒng)計(jì)量, 也難以通過(guò)理論分析方法來(lái)求解.近幾十年來(lái), 隨著計(jì)算機(jī)計(jì)算技術(shù)與數(shù)值方法的快速發(fā)展, 數(shù)值模擬已經(jīng)成為一種基本的研究手段, 在RT不穩(wěn)定性研究方面發(fā)揮了越來(lái)越重要的作用.He等[19]采用一個(gè)雙分布等溫格子 Boltzmann (Lattice Boltzmann, LB)方法模擬了三維單模RT不穩(wěn)定性問(wèn)題, 分析了雷諾數(shù)和Atoowd數(shù)對(duì)不穩(wěn)定性界面結(jié)構(gòu)的影響.Glimm等[20]基于歐拉方程的前追蹤方法研究了二維單模RT不穩(wěn)定性問(wèn)題, 發(fā)現(xiàn)不穩(wěn)定性的發(fā)展經(jīng)歷了飽和恒定速度階段后會(huì)出現(xiàn)一個(gè)加速階段, 這個(gè)階段后來(lái)被Wilkinson和Jacobs[18]實(shí)驗(yàn)所發(fā)現(xiàn).Celani[21]利用相場(chǎng)方法模擬了二維單模RT不穩(wěn)定性問(wèn)題, 并在初始階段獲得了與線性理論相符合的擾動(dòng)增長(zhǎng)率.需要指出的是, 大多數(shù)的數(shù)值研究是用于驗(yàn)證所發(fā)展的計(jì)算流體力學(xué)方法的正確性與有效性的, 研究階段也局限于不穩(wěn)定性發(fā)展的中前期.Ramaprabhu等[22]模擬了三維混相流體的單模RT不穩(wěn)定性現(xiàn)象, 考察了Atwood數(shù)和雷諾數(shù)對(duì)氣泡和尖釘振幅在后期階段的演化規(guī)律.他們發(fā)現(xiàn)低Atwood數(shù)的不穩(wěn)定性在高雷諾數(shù)情形下經(jīng)歷了指數(shù)增長(zhǎng)、飽和速度增長(zhǎng)、重加速、混沌混合等四個(gè)發(fā)展階段, 并進(jìn)一步指出在重加速階段,氣泡和尖釘速度已經(jīng)超過(guò)經(jīng)典勢(shì)能模型所預(yù)測(cè)的理論解, 而在后期階段, 氣泡和尖釘速度會(huì)出現(xiàn)急劇的下降.Wei等[23]采用直接數(shù)值模擬方法研究了低Atwoods數(shù)下二維混相流體的單模RT不穩(wěn)定性問(wèn)題, 分析了雷諾數(shù)對(duì)氣泡和尖釘振幅增長(zhǎng)的影響.他們同樣地觀察到混相不穩(wěn)定性在高雷諾數(shù)時(shí)經(jīng)歷了一系列的發(fā)展階段, 但不同于Ramaprabhu等[22]的結(jié)果, 氣泡演化速度在后期的混沌階段出現(xiàn)了隨時(shí)間波動(dòng)的現(xiàn)象, 并且其平均加速度接近于常數(shù), 這表明氣泡在演化后期具有二次增長(zhǎng)的規(guī)律.Liang等[24,25]利用基于相場(chǎng)理論的多相流LB方法模擬二維及三維長(zhǎng)微管道內(nèi)低Atwoods數(shù)下非混相流體的單模RT不穩(wěn)定性問(wèn)題, 詳細(xì)地考察了雷諾數(shù)對(duì)不穩(wěn)定性各發(fā)展階段擾動(dòng)增長(zhǎng)和流體界面動(dòng)態(tài)行為的影響.數(shù)值結(jié)果表明高雷諾數(shù)時(shí), 不穩(wěn)定性也經(jīng)歷了四個(gè)發(fā)展階段, 包括線性增長(zhǎng)階段、飽和速度增長(zhǎng)階段、重加速階段及混沌混合階段.在后期階段, 相界面發(fā)生多層次卷起、劇烈變形、混沌破裂等行為, 形成了非常復(fù)雜的拓?fù)浣Y(jié)構(gòu), 并且氣泡振幅顯示二次增長(zhǎng)的規(guī)律.而低雷諾數(shù)時(shí), 相界面的演化變得相對(duì)光滑,后期的發(fā)展階段也相繼難以到達(dá).最近, Hu等[26]數(shù)值研究了低Atwoods數(shù)混相流體的單模RT不穩(wěn)定性問(wèn)題, 發(fā)現(xiàn)氣泡速度在重加速階段后會(huì)出現(xiàn)不停地加速與減速, 并將這種現(xiàn)象歸因于旋渦強(qiáng)度的變化.

綜上所述, 已有許多學(xué)者[16–25]對(duì)單模RT不穩(wěn)定性問(wèn)題開(kāi)展了研究, 豐富了人們對(duì)不穩(wěn)定性發(fā)展規(guī)律的認(rèn)識(shí).然而, 絕大多數(shù)的工作[16–21]局限于不穩(wěn)定性發(fā)展的前期階段, 并且所考慮的雷諾數(shù)均較小.盡管一些學(xué)者[22–26]對(duì)高雷諾數(shù)下不穩(wěn)定性演化的后期階段開(kāi)展了部分研究, 但考慮的兩種流體為相互混溶的[22,23,26], 所關(guān)注的流體間Atwood數(shù)也太小[23–26].鑒于此, 本文將采用基于相場(chǎng)理論的多松弛LB方法[24]研究中等Atwood數(shù)下非混相流體RT不穩(wěn)定性的演化規(guī)律, 著重分析雷諾數(shù)對(duì)相界面動(dòng)力學(xué)行為和擾動(dòng)后期增長(zhǎng)的影響.

2 數(shù)值方法

基于分子動(dòng)理學(xué)理論的LB方法[27]是近十幾年發(fā)展起來(lái)的介觀數(shù)值方法, 它不再基于宏觀連續(xù)介質(zhì)模型, 而通過(guò)描述流體粒子分布函數(shù)的演化再現(xiàn)復(fù)雜流動(dòng)的宏觀行為, 從而相比傳統(tǒng)數(shù)值方法具有一些獨(dú)特的優(yōu)勢(shì), 比如易于處理復(fù)雜物理邊界、程序?qū)崿F(xiàn)簡(jiǎn)單且天然并行、直觀刻畫(huà)多相流體間微觀相互作用.目前, 從流體組分間相互作用力的不同物理背景出發(fā), 已經(jīng)提出了多種描述多相流體輸運(yùn)的LB模型, 包括顏色模型、偽勢(shì)模型、自由能模型、基于相場(chǎng)理論的LB方法.而基于相場(chǎng)理論的LB方法在描述多相流體界面動(dòng)力學(xué)方面具有清晰的物理機(jī)制, 可用于模擬多相系統(tǒng)中界面具有拓?fù)渥兓牧鲃?dòng), 因而受到了許多學(xué)者的廣泛關(guān)注, 并已成功地應(yīng)用于軸對(duì)稱(chēng)多相流[28,29]、三相流[30,31]、復(fù)雜微通道內(nèi)液滴動(dòng)力學(xué)[32–34]、液滴潤(rùn)濕固體表面[35]等復(fù)雜多相流問(wèn)題.本文將采用Liang等[24]提出的相場(chǎng)LB模型研究長(zhǎng)微通道中非混相RT不穩(wěn)定性的后期演化規(guī)律, 該模型相比前人相場(chǎng)LB模型可正確地恢復(fù)Cahn-Hilliard和Navier-Stokes的耦合方程, 并且速度和壓力場(chǎng)可通過(guò)分布函數(shù)進(jìn)行顯示計(jì)算.另外, 為了提高模型在高雷諾數(shù)時(shí)的數(shù)值穩(wěn)定性, 我們還在相場(chǎng)LB模型的演化過(guò)程中采用多松弛的碰撞算子[36].本文相場(chǎng)LB模型利用兩個(gè)獨(dú)立的分布函數(shù) fi和 gi,其對(duì)應(yīng)的多松弛演化方程可表述為:

其中 fi(x,t) 是描述粒子在空間 x 和時(shí)間t時(shí)的序參數(shù)分布函數(shù), gi(x,t) 是密度分布函數(shù),為兩分布函數(shù)所對(duì)應(yīng)的平衡態(tài)分布函數(shù),M 是分布函數(shù)空間到矩空間的變換矩陣, Sf, Sg為松弛矩陣, Fi(x,t) 、 Gi(x,t) 分別為源項(xiàng)和外力項(xiàng)的分布函數(shù), δx和 δt為空間和時(shí)間步長(zhǎng).為了匹配宏觀控制方程, 平衡態(tài)分布函數(shù)和可設(shè)定為[24]

其中 ? 是序參數(shù), μ 是化學(xué)勢(shì), η 為調(diào)節(jié)遷移率的自由參數(shù), ρ 是密度, p 是流體動(dòng)力學(xué)壓力, u 是流體速度, ωi為權(quán)系數(shù), cs為聲速.根據(jù)相場(chǎng)理論, 范德瓦爾斯 (van der Waals)流體的化學(xué)勢(shì) μ 可表述為序參數(shù)的相關(guān)函數(shù), 具體定義為

其中參數(shù)k, β 與界面厚度 (D) 、表面張力 (σ) 相關(guān),

另外, ωi和的選取依賴(lài)于格子速度的離散模型.針對(duì)本文所研究的二維微通道流動(dòng), 我們采用最有代表性的D2 Q9格子離散速度模型[37,38], 其權(quán)系數(shù) 和 聲 速 為 ω0=4/9 , ω1?4=1/9 , ω5?8=1/36 ,離散速度 ei為

且采用規(guī)則的正方形格子, 其對(duì)應(yīng)的變換矩陣為[36]

同樣地, 在多松弛的 D2Q9 格子模型, fi, gi所對(duì)應(yīng)的松弛矩陣分別定義為:

其中 I 是單位矩陣, Fs是界面張力, G 是外力項(xiàng).在相場(chǎng)模型中, 表面張力取為勢(shì)能形式 Fs= μ?? ,發(fā)現(xiàn)能夠有效地減少相界面處的虛假速度.

在本模型中, 流體的宏觀量可以通過(guò)求解粒子分布函數(shù)的各階矩獲得, 具體的計(jì)算式如下[39]:

而宏觀量密度(ρ)可以看成是關(guān)于序參數(shù)(?)的一個(gè)簡(jiǎn)單的線性插值函數(shù),

其中 ρh和 ρl分別為重質(zhì)流體和輕質(zhì)流體的密度.

最后, 通過(guò)Chapman-Enskog多尺度理論分析[24,39], 可以證明本文采用的基于相場(chǎng)理論的多松弛LB模型可以正確地恢復(fù)界面追蹤的Cahn-Hilliard控制方程,

和描述流體動(dòng)力學(xué)的不可壓Navier-Stokes方程組

其中 遷移率M、流體運(yùn)動(dòng)性黏性 ν 與松弛 因子關(guān)系可分別表示為:

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

本文將采用基于相場(chǎng)理論的多松弛LB模型研究二維長(zhǎng)管道內(nèi)非混相流體的單模RT界面不穩(wěn)定性問(wèn)題, 并詳細(xì)地考察雷諾數(shù)對(duì)界面不穩(wěn)定性發(fā)展的影響, 以及定量分析氣泡與尖釘?shù)恼穹S著雷諾數(shù)的變化規(guī)律.為了研究不穩(wěn)定性的后期演化規(guī)律, 本文考慮的物理問(wèn)題為一個(gè)長(zhǎng)方形的微通道, 其高度和寬度分別為 H 和 W, 且 H /W=8 .初始時(shí)刻, 密度較大的重流體置于另一種輕質(zhì)流體的上方, 且給流體界面處一個(gè)微小的擾動(dòng), 在重力場(chǎng)的作用下, 擾動(dòng)會(huì)逐漸發(fā)展, 并最終達(dá)到混沌混合狀態(tài).在長(zhǎng)方形微通道的中心截面處, 施加一個(gè)波長(zhǎng)為W具有余弦函數(shù)的微小初始擾動(dòng),

其中 k =2π/W 是擾動(dòng)波數(shù).為了使序參數(shù)變量在相界面處光滑連續(xù)的變化, 設(shè)定序參數(shù) ? (x,y) 的初始分布為

為了表征二維RT不穩(wěn)定性的演化特征, 引了兩個(gè)常用的無(wú)量綱參數(shù), 即雷諾數(shù)(Re)和阿特伍德數(shù)(A t), 分別定義為[12,13]

其中g(shù)是重力加速度, ν 是流體運(yùn)動(dòng)性黏性.在本文的研究中, 我們考慮中等的阿特伍德數(shù)下界面不穩(wěn)定性的演化規(guī)律, 將輕重流體的密度分別設(shè)定為 1和3, 其對(duì)應(yīng)的 Atwoods數(shù)為 0.5, 其他的物理參數(shù)選取如下: W =256 , D =4 , σ =5×10?5,為了演化計(jì)算, 需要選擇合適的邊界格式處理物理邊界處的粒子分布函數(shù), 本文對(duì)上下固壁采用無(wú)滑移的半反彈邊界條件, 左右邊界應(yīng)用周期邊界條件.另外, 標(biāo)記氣泡和尖釘?shù)恼穹鶠镠b和 Hs, 分別定義為氣泡與尖釘?shù)那岸说綄?duì)應(yīng)初始位置的距離, 因此氣泡和尖釘?shù)恼穹诔跏紩r(shí)刻的值為0.進(jìn)一步, 我們將管道寬度W和選為特征長(zhǎng)度和特征時(shí)間, 還定義了氣泡和尖釘無(wú)量綱的演化速度, 也被稱(chēng)為氣泡和尖釘?shù)腇roude數(shù)[23],

其中 ub和 us為氣泡和尖釘前端的速度, 可以由氣泡和尖釘振幅計(jì)算獲得.除了特別聲明, 本文接下來(lái)給出的長(zhǎng)度、速度、時(shí)間(τ)等相關(guān)物理量均已被相應(yīng)的特征值所無(wú)量綱化.

圖1描述了四種典型的不同雷諾數(shù)下非混相RT不穩(wěn)定性中相界面的演化過(guò)程.從圖1中可以發(fā)現(xiàn), 對(duì)于不同的雷諾數(shù), 不穩(wěn)定性的擾動(dòng)在初始階段顯示相似的界面動(dòng)力學(xué)特行為: 重流體在重力作用往下運(yùn)動(dòng)而輕流體向上浮起, 即輕重流體之間相互滲透, 從而形成了氣泡和尖釘圖案.緊接著,流體界面在不同的雷諾數(shù)下展示出顯著不同的動(dòng)力學(xué)特征.對(duì)高 Re 情形 (Re = 10000), 尖釘繼續(xù)向下運(yùn)動(dòng)并逐漸地向上卷起, 形成了旋轉(zhuǎn)方向相反的兩個(gè)旋渦, 這是KH不穩(wěn)定性出現(xiàn)并作用于相界面的結(jié)果.隨著時(shí)間的演化, 兩個(gè)旋渦不斷地增長(zhǎng),在卷起的尾端處形成了一對(duì)二級(jí)旋渦.隨后, 在多個(gè)旋渦相互作用下, 不穩(wěn)定性系統(tǒng)的非線性效應(yīng)越來(lái)越劇烈, 尖釘卷起的長(zhǎng)度也越來(lái)越長(zhǎng), 并逐漸地靠近進(jìn)而接觸中軸線附近的流體界面.在高流體界面剪切力作用下, 中軸線附近的界面在多處位置出現(xiàn)卷起與變形行為.最終, 流體界面發(fā)生了混沌的破裂, 在系統(tǒng)中形成了許多離散的小液滴.另外,我們還觀察到流體界面在整個(gè)演化過(guò)程中始終保持關(guān)于中軸線對(duì)稱(chēng).當(dāng)Re減少至2048時(shí), 同樣觀察到尖釘發(fā)生卷起行為, 在尾端也形成了一對(duì)二級(jí)渦, 并最終導(dǎo)致相界面在多個(gè)位置發(fā)生卷吸、變形和破裂, 形成較為復(fù)雜的結(jié)構(gòu).然而, 相比 Re =10000的情形, 在演化后期, 系統(tǒng)中相界面的混沌程度減弱了.當(dāng)Re數(shù)進(jìn)一步降低至50時(shí), 重流體的尖釘往下運(yùn)動(dòng), 經(jīng)過(guò)一段演化時(shí)間, 在尾端發(fā)生卷吸現(xiàn)象, 也形成了一對(duì)旋轉(zhuǎn)方向相互相反的旋渦, 但與高雷諾數(shù)情形相比, 界面卷起發(fā)生的時(shí)刻推遲了, 旋渦的卷吸幅度也相應(yīng)地減弱了.最后,形成的旋渦隨著時(shí)間演化而不斷地發(fā)展, 伴隨著尾端卷起的部分也越來(lái)越長(zhǎng).在整個(gè)演化過(guò)程中, 未觀察到二次旋渦卷吸和界面后期發(fā)生破裂的現(xiàn)象.當(dāng) Re 充分小 (Re = 5)時(shí), 卷吸現(xiàn)象不再發(fā)生, 重流體將以尖釘?shù)姆绞讲粩嗟叵蛳逻\(yùn)動(dòng), 界面也變得足夠光滑, 未出現(xiàn)混沌破裂等復(fù)雜拓?fù)洮F(xiàn)象, 這是由于強(qiáng)黏性作用使流體之間的剪切層保持穩(wěn)定, 流動(dòng)在整個(gè)過(guò)程表現(xiàn)為層流狀態(tài).

上面討論了雷諾數(shù)對(duì)單模RT不穩(wěn)定性中相界面動(dòng)力學(xué)行為的影響, 而氣泡與尖釘振幅及演化速度是描述RT不穩(wěn)定性問(wèn)題中另外兩個(gè)非常重要的物理量.為了進(jìn)一步顯示雷諾數(shù)的效應(yīng), 我們定量地分析了不同雷諾數(shù)下氣泡和尖釘振幅、運(yùn)動(dòng)速度隨時(shí)間的演化規(guī)律.圖2分別給出了氣泡與尖釘在不同雷諾數(shù)下隨時(shí)間變化的演化曲線.從圖中可以發(fā)現(xiàn), 對(duì)所有的Re數(shù)情形, 氣泡和尖釘?shù)恼穹S著時(shí)間演化而不斷增大.當(dāng)Re逐漸增大時(shí),可以觀察到同一時(shí)刻所獲得的尖釘振幅也越大, 而當(dāng)Re增大至足夠大時(shí), 雷諾數(shù)對(duì)尖釘振幅的影響將不再顯著.在不可壓流體的RT不穩(wěn)定性中, 尖釘?shù)倪\(yùn)動(dòng)特征在理論上由單位質(zhì)量的浮力和黏性耗散力之間的競(jìng)爭(zhēng)關(guān)系決定[40],

圖1 雷諾數(shù)對(duì)非混相 RT 不穩(wěn)定性中相界面演化圖案的影響 (a) Re=10000 ; (b) Re=2048 ; (c) Re=50 ; (d) Re=5Fig.1.The effect of the Reynolds number on the evolution of interfacial patterns in the immiscible RT instability: (a) Re=10000 ;(b) Re=2048 ; (c) R e=50 ; (d) Re=5 .

其中 δ ρ =(ρh? ρl)/2 , ρˉ =(ρh+ρl)/2 , Fd是黏性耗散力, 定義為重力方向上的動(dòng)量耗散率Fd= ??/ν, ? 是能量耗散率.對(duì)于無(wú)黏時(shí)間尺度,Sreenivasan[41]理論分析給出了能量耗散率的數(shù)學(xué)表達(dá)式, ? =Cν3/W , 其中 C 是常數(shù).根據(jù)上述分析, 在雷諾數(shù)從小增大的過(guò)程中, 黏性耗散力在不斷地減少, 從而理論上可以獲得更大的尖釘運(yùn)動(dòng)速度, 以及更大的尖釘?shù)恼穹? 而當(dāng)雷諾數(shù)足夠大時(shí),黏性耗散力已充分小, 浮力相比黏性耗散力更占統(tǒng)治地位, 從而繼續(xù)增大雷諾數(shù)對(duì)尖釘振幅增長(zhǎng)的影響將不再顯著.在本文模擬中, 我們確實(shí)觀察到與上述的理論分析相一致的數(shù)值結(jié)果.另外, 從圖2中還可以發(fā)現(xiàn), 當(dāng)雷諾數(shù)在一定范圍內(nèi), 增大雷諾數(shù)可以有效地促進(jìn)氣泡振幅的增長(zhǎng), 而當(dāng)雷諾數(shù)充分大時(shí), 氣泡振幅的增長(zhǎng)在后期會(huì)隨著雷諾數(shù)的增大而呈現(xiàn)出一種遞減的趨勢(shì).這是由于當(dāng)雷諾數(shù)足夠大時(shí), 在浮力驅(qū)動(dòng)的不穩(wěn)定性演化中后期, 誘導(dǎo)產(chǎn)生了KH不穩(wěn)定性.受兩類(lèi)不穩(wěn)定性的共同作用, 系統(tǒng)中出現(xiàn)了許多不同尺度的旋渦, 這些渦效應(yīng)在一定程序上減緩了氣泡向上運(yùn)動(dòng).

圖2 雷諾數(shù)對(duì)無(wú)量綱化的氣泡與尖釘隨時(shí)間演化振幅的影響Fig.2.The effect of the Reynolds number on the dimensionless bubble and spike amplitudes.

進(jìn)一步, 我們還定量地統(tǒng)計(jì)了不同雷諾數(shù)下氣泡和尖釘隨時(shí)間演化的運(yùn)動(dòng)速度 ub, us, 其大小是通過(guò)氣泡和尖釘?shù)恼穹鶎?duì)演化時(shí)間差分計(jì)算獲得的, 并根據(jù)(25)式無(wú)量綱化得到氣泡與尖釘?shù)腇roude數(shù).圖3描述了不同雷諾數(shù)下氣泡與尖釘隨時(shí)間演化的Froude數(shù).從圖3中可以發(fā)現(xiàn), 氣泡與尖釘?shù)倪\(yùn)動(dòng)速度在不同雷諾數(shù)下表現(xiàn)出顯著不同的演化規(guī)律.雷諾數(shù)越高, 尖釘往下運(yùn)動(dòng)的速度越快,而當(dāng)雷諾數(shù)充分大時(shí), 尖釘?shù)难莼俣葘?duì)雷諾數(shù)的變化不再敏感, 這與上述(26)式理論分析的結(jié)果是相統(tǒng)一的.雷諾數(shù)對(duì)氣泡運(yùn)動(dòng)速度的影響則顯現(xiàn)出先促進(jìn)而后抑制的規(guī)律, 這是由于當(dāng)雷諾數(shù)充分高時(shí), 在不穩(wěn)定性的演化后期產(chǎn)生了許多不同尺度的旋渦, 這些旋渦的運(yùn)動(dòng)減弱了氣泡往上運(yùn)動(dòng)的速率而使之發(fā)生旋轉(zhuǎn).Ramaprabhu等[22]數(shù)值研究了混相流體的RT不穩(wěn)定性的后期演化規(guī)律, 他們也發(fā)現(xiàn)當(dāng)雷諾數(shù)足夠大時(shí), 繼續(xù)增大雷諾數(shù)會(huì)減小氣泡往上的演化速度.另外, 我們進(jìn)一步對(duì)氣泡與尖釘速度演化圖進(jìn)行分析, 將高雷諾數(shù)下非混相流體的單模RT不穩(wěn)定性的演化歸結(jié)于四個(gè)發(fā)展階段, 包括線性增長(zhǎng)階段、飽和速度增長(zhǎng)階段、重加速階段、混沌混合階段.在初始階段, 不穩(wěn)定性的擾動(dòng)發(fā)展符合線性穩(wěn)定性理論, 其振幅的增長(zhǎng)具有指數(shù)形式的規(guī)律[16],

圖3 雷諾數(shù)對(duì)無(wú)量綱化的氣泡和尖釘演化速度的影響Fig.3.The effect of the Reynolds number on the dimensionless bubble and spike velocities.

其中H代表氣泡或者尖釘?shù)恼穹? a1與 a2是擬合參數(shù), γ 是線性增長(zhǎng)因子.圖4給出了不同雷諾數(shù)下隨時(shí)間演化的氣泡與尖釘振幅的數(shù)值模擬結(jié)果以及曲線擬合結(jié)果, 可以發(fā)現(xiàn)氣泡和尖釘?shù)某跏荚鲩L(zhǎng)確實(shí)符合線性穩(wěn)定性理論, 并且獲得的線性增長(zhǎng)因子隨著雷諾數(shù)的增大而增大.緊接著線性增長(zhǎng)階段, 氣泡與尖釘將以近似恒定的速度增長(zhǎng), 這表明不穩(wěn)定性的發(fā)展進(jìn)入飽和速度增長(zhǎng)階段.Goncharov[17]解析分析了單模RT不穩(wěn)定性的非線性增長(zhǎng)區(qū)域,提出了經(jīng)典的勢(shì)能理論模型以預(yù)測(cè)氣泡與尖釘?shù)娘柡驮鲩L(zhǎng)速度, 其表達(dá)式如(1)式所示.進(jìn)一步根據(jù)(1)式, 可以推導(dǎo)出氣泡與尖釘在飽和速度階段所對(duì)應(yīng)的無(wú)量綱Froude數(shù)分別為0.325與0.564.從圖3可以發(fā)現(xiàn), 高雷諾數(shù)下尖釘在飽和速度階段的Froude數(shù)略高于勢(shì)能模型的解析解, 這是由于在實(shí)際模擬中, 界面在該階段發(fā)生卷吸行為, 產(chǎn)生了許多不同尺度的旋渦, 這些渦效應(yīng)會(huì)促進(jìn)尖釘?shù)陌l(fā)展, 而勢(shì)能模型的理論解未包含渦效應(yīng).Goncharov[17]通過(guò)數(shù)值模擬同樣驗(yàn)證了尖釘實(shí)際演化速度高于勢(shì)能模型的解析解.另外, 我們發(fā)現(xiàn)在高雷諾數(shù)條件下, 繼續(xù)增大雷諾數(shù)會(huì)減少氣泡演化速度, 從而導(dǎo)致氣泡飽和速度小于勢(shì)能模型的解析解.接下來(lái), 各尺度的渦結(jié)構(gòu)之間相互作用逐漸增加, 使得氣泡和尖釘Froude數(shù)高于經(jīng)典勢(shì)能模型的理論解, 這預(yù)示著不穩(wěn)定發(fā)展進(jìn)入了重加速階段.重加速階段不能持續(xù)地發(fā)展下去, 在演化后期,氣泡與尖釘?shù)腇roude數(shù)變得不穩(wěn)定, 開(kāi)始隨著時(shí)間波動(dòng), 這表明界面的演化進(jìn)入了混沌混合階段.為了揭示混沌混合階段不穩(wěn)定性的發(fā)展規(guī)律, 我們計(jì)算了氣泡與尖釘增長(zhǎng)的無(wú)量綱加速度,其中表示氣泡與尖釘振幅對(duì)時(shí)間的兩階導(dǎo)數(shù), 實(shí)際通過(guò)對(duì)氣泡與尖釘振幅關(guān)于時(shí)間的二階差分計(jì)算獲得.圖5給出了高雷諾下氣泡與尖釘無(wú)量綱加速度隨時(shí)間的演化曲線.從圖可以發(fā)現(xiàn), 氣泡與尖釘加速度在演化后期不穩(wěn)定性,分別繞著常數(shù)0.045與0.233上下波動(dòng), 預(yù)示著后期氣泡與尖釘平均加速度是一個(gè)常數(shù), 并表明RT不穩(wěn)定性的后期發(fā)展呈現(xiàn)出二次增長(zhǎng)的規(guī)律.當(dāng)雷諾數(shù)足夠低時(shí), 在不穩(wěn)定性的整個(gè)演化過(guò)程中, 不能觀察到重加速階段與混沌混合階段, 氣泡與尖釘在后期階段將以恒定的飽和速度增長(zhǎng).另外, 我們發(fā)現(xiàn)低雷諾數(shù)下氣泡與尖釘?shù)娘柡退俣鹊陀诮?jīng)典的勢(shì)能模型的理論解, 其原因在于勢(shì)能模型考慮的是理想無(wú)黏性流體的不穩(wěn)定性現(xiàn)象, 未考慮流體黏性對(duì)演化速度的影響.

圖4 不同雷諾數(shù)下, 氣泡和尖釘振幅在初始階段的演化曲線, 其中數(shù)據(jù)點(diǎn)是統(tǒng)計(jì)結(jié)果, 實(shí)線則是擬合結(jié)果Fig.4.The curves of the early-time bubble and spike amplitudes with different Reynolds numbers, where the data points and solid lines are the statistical and fitting results.

圖5 高雷諾數(shù)下, 氣泡和尖釘?shù)臒o(wú)量綱加速度演化曲線,紅色和藍(lán)色實(shí)線分別為0.045和0.233Fig.5.The curves of dimensionless bubble and spike accelerations at a high Reynolds number, and the red and blue solid lines are 0.045 and 0.233.

4 結(jié) 論

本文基于相場(chǎng)理論的LB方法模擬了一個(gè)長(zhǎng)微管道內(nèi)非混相流體的RT不穩(wěn)定性問(wèn)題, 該方法可以準(zhǔn)確地追蹤相界面動(dòng)力學(xué)行為, 以及采用多松弛碰撞模型可以很好地處理高雷諾數(shù)流動(dòng).本文著重分析雷諾數(shù)對(duì)中等Atwoods數(shù)的不穩(wěn)定性演化中界面動(dòng)力學(xué)行為和擾動(dòng)增長(zhǎng)規(guī)律的影響.研究發(fā)現(xiàn)在高雷諾數(shù)情形下, RT不穩(wěn)定性的發(fā)展先后經(jīng)歷線性增長(zhǎng)階段、飽和速度增長(zhǎng)階段、重加速階段以及混沌混合階段.在線性階段, 重流體與輕流體之間相互滲透, 擾動(dòng)的增長(zhǎng)符合經(jīng)典的線性穩(wěn)定性理論.緊接著, 不穩(wěn)定性的演化進(jìn)入了飽和速度階段.在該階段中, 輕重流體形成了氣泡與尖釘?shù)慕Y(jié)構(gòu), 氣泡與尖釘將以恒定的速度增長(zhǎng), 并隨著時(shí)間演化在尖釘尾端處出現(xiàn)卷吸行為.隨著橫向速度和縱向速度的差異逐漸擴(kuò)大, 氣泡與尖釘?shù)难莼T導(dǎo)產(chǎn)生了KH不穩(wěn)定性, 進(jìn)而使系統(tǒng)中出現(xiàn)了許多不同尺度的渦結(jié)構(gòu), 從而加速了氣泡與尖釘?shù)难莼俣?重加速階段不能一直持續(xù)下去, 氣泡與尖釘?shù)难莼俣仍诤笃陔A段會(huì)發(fā)生波動(dòng)現(xiàn)象, 這表明不穩(wěn)定性的發(fā)展進(jìn)入了混沌混合階段.在混沌混合階段, 界面會(huì)發(fā)生多層次卷起、劇烈變形、混沌破裂,形成了非常復(fù)雜的界面拓?fù)浣Y(jié)構(gòu), 系統(tǒng)中也產(chǎn)生了許多離散的小液滴.另外, 我們統(tǒng)計(jì)了演化后期的氣泡與尖釘加速度, 其無(wú)量綱的值分別圍繞著常數(shù)0.045與0.233上下波動(dòng), 這表明不穩(wěn)定性在演化后期具有二次增長(zhǎng)的規(guī)律.而當(dāng)雷諾數(shù)足夠小,擾動(dòng)的增長(zhǎng)變得非常緩慢, 流體相界面也變得足夠光滑, 未出現(xiàn)卷吸和破裂行為, 也未觀察到后期的重加速與混沌混合階段, 氣泡與尖釘將以恒定速度一直演化下去.