高雷諾數(shù)下非混相Rayleigh-Taylor不穩(wěn)定性的格子Boltzmann方法模擬*

胡曉亮 梁宏? 王會利

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

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

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

1 引 言

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

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

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

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

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

2 數(shù)值方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 不同雷諾數(shù)下, 氣泡和尖釘振幅在初始階段的演化曲線, 其中數(shù)據(jù)點是統(tǒng)計結(jié)果, 實線則是擬合結(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量綱加速度演化曲線,紅色和藍色實線分別為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é) 論

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