黃劍霖 宋廣毅 王靜竹*,*?,1) 王一偉,**

* (中國科學院力學研究所流固耦合系統(tǒng)力學重點實驗室,北京 100190)

? (中國科學院大學未來技術(shù)學院,北京 100190)

** (中國科學院大學工程科學學院,北京 100190)

?? (廣東空天科技研究院,廣州 511458)

引言

Rayleigh-Taylor (RT)不穩(wěn)定性是一種液體或氣體界面上出現(xiàn)的不穩(wěn)定現(xiàn)象.當兩個不同密度的流體在接觸面上形成密度梯度時,通過重力或其他外部力的作用,界面發(fā)生擾動和不規(guī)則變形,并產(chǎn)生波浪、褶皺或渦旋等[1-2].RT 不穩(wěn)定性現(xiàn)象在自然界和工程應用中廣泛存在[3-5].例如,在海洋中不同密度的海水層之間的相互作用可能導致海浪的形成和傳播[6-7].在火山噴發(fā)過程中,熔巖和氣體之間的界面會出現(xiàn)RT 不穩(wěn)定性,導致噴發(fā)的物質(zhì)和形狀的變化[8-9].在慣性約束聚變中,兩種不同密度的激光束界面可能引RT 不穩(wěn)定性,這可能導致激光能量的損失和散射,從而影響聚變反應的效率[10-11].在航空航天領(lǐng)域,飛行器在大氣中飛行時,機翼表面的氣流分離可能導致RT 不穩(wěn)定性,從而影響飛行器的性能和穩(wěn)定性[12].因此,理解與定向控制RT 不穩(wěn)定性對于工程應用來說是至關(guān)重要的.

Taylor 等[13-14]提出平面幾何界面線性不穩(wěn)定性分析理論,結(jié)合實驗總結(jié)了RT 不穩(wěn)定性的發(fā)展的幾個階段: (1)線性階段,擾動振幅遠小于擾動波長,擾動隨時間指數(shù)增長,可以用線性理論描述;(2)變形階段,擾動振幅發(fā)展,界面開始表現(xiàn)出不對稱特征[15];(3)規(guī)則非線性階段,低密度流體形成空泡進入高密度流體,高密度流體形成尖釘進入低密度流體中.Birkhoff[16]在其基礎上總結(jié)出了規(guī)則非線性后的不穩(wěn)定性發(fā)展階段;(4)不規(guī)則非線性發(fā)展階段,界面發(fā)展受到速度差異等多種作用耦合影響;(5)湍流混合階段,界面兩側(cè)流體會發(fā)生摻混.

同時,Bell[17]和Plesset[18]分別給出了柱面和球面幾何下的RT 不穩(wěn)定性演化方程,建立了運動幾何效應影響下的不穩(wěn)定性理論.這種界面運動軌跡影響界面擾動發(fā)展的效應被稱為Bell-Plesset (BP)效應.BP 效應下界面不穩(wěn)定性在工程中同樣廣泛存在,比如預混燃燒[19]、爆炸[20-21]和液體噴霧[22],可以簡化成內(nèi)部體積振蕩誘導液滴界面變形這一模型.Obreschkow 等[23]通過高速攝像機觀察了微重力條件下球形液滴中空泡振蕩誘導界面擾動、變形等現(xiàn)象.呂明等[24]基于數(shù)值模擬研究了液滴內(nèi)空泡生長的規(guī)律,提出了空泡潰滅的3 個歷程,包括快速潰滅期、緩慢潰滅期以及穩(wěn)定期.Zeng 等[25]基于實驗和數(shù)值模擬研究了液滴內(nèi)空泡潰滅產(chǎn)生的射流,他們提出內(nèi)部空泡振蕩產(chǎn)生的徑向加速度誘導了界面RT 不穩(wěn)定性的形成.為了在更大參數(shù)空間評估界面穩(wěn)定性,推導了包含液體黏性的球面RT 不穩(wěn)定性擾動增長解析模型,結(jié)果發(fā)現(xiàn)發(fā)現(xiàn)空泡潰滅后產(chǎn)生的表面射流會隨著液體黏性增加而減少.Wang等[26]鑒于實驗觀測手段難以捕捉三維不穩(wěn)定性的發(fā)展,將模型簡化至柱坐標系,結(jié)合實驗、數(shù)值模擬和理論全面研究了板間液滴下空泡潰滅誘導的RT 不穩(wěn)定性現(xiàn)象,他們總結(jié)了對于不同液滴初始半徑和空泡初始內(nèi)壓工況下液滴出現(xiàn)的破碎、通氣和穩(wěn)定3 種形態(tài)特征(如圖1 所示),同時建立了考慮RT 不穩(wěn)定性和空泡振蕩的分析模型,用于分析擾動的增長和液滴形態(tài)的評估.Zhang 等[27]對不同空泡和液滴半徑比的飛濺情況進行分析,結(jié)果表明隨著半徑比的增加,射流飛濺效果會更加顯著,在這種情況下空泡的潰滅時間會變小.Rosselló等[28]對自由落體液滴內(nèi)空泡潰滅的界面現(xiàn)象展開研究,發(fā)現(xiàn)空泡潰滅誘導液滴表面出現(xiàn)RT 不穩(wěn)定性和Rayleigh-Plateau 不穩(wěn)定性,并提出了不穩(wěn)定性產(chǎn)生的機制: 空泡快速膨脹或收縮加速液層導致液滴表面波紋產(chǎn)生,波紋在界面進一步發(fā)展從而觀察到不穩(wěn)定性現(xiàn)象.

圖1 不同半徑二維液滴內(nèi)3 種形態(tài)Fig.1 Observations of three distinct deformation characteristics for 2-D water droplets with different diameters

在工程應用中,液體黏性是影響不穩(wěn)定性產(chǎn)生和發(fā)展的關(guān)鍵因素,對應用的效率和可靠性有著重要的影響.Zeng 等[25]提出了考慮黏性的球形液滴內(nèi)空泡誘導RT 不穩(wěn)定性形成的機理,但是在黏性對液滴特征形態(tài)和不穩(wěn)定性評估的影響規(guī)律尚未研究清楚.本文采用直接數(shù)值模擬和高精度界面捕捉方法,研究了液體黏性影響下柱狀液滴內(nèi)空泡誘導不穩(wěn)定性的產(chǎn)生與發(fā)展,以期為慣性約束核聚變等工程應用提供參考.

1 數(shù)值模擬方法

1.1 計算模型

本文基于開源框架OpenFOAM 建立直接數(shù)值模擬方法,開展液體黏性對內(nèi)空泡振蕩誘導柱狀液滴界面不穩(wěn)定性的影響規(guī)律研究.

計算模型和邊界條件如圖2 所示,計算域內(nèi)液滴和空泡分別用藍色和白色表示.為了模擬空泡的振蕩,初始空泡初始設置成高溫高壓狀態(tài).模擬過程忽略軸向的變化而重點關(guān)注徑向的不穩(wěn)定性現(xiàn)象.計算域為三維矩形結(jié)構(gòu),選取1/4 對稱結(jié)構(gòu)開展計算,其大小為12 mm×12 mm×1 mm,右面和下面為對稱面,上面和左面為出口條件,前面和后面為空邊界.全區(qū)域采用結(jié)構(gòu)性網(wǎng)格均勻劃分,最小網(wǎng)格尺寸Δx=20 μm,其中初始空泡內(nèi)包含20 個網(wǎng)格.

圖2 計算域及網(wǎng)格劃分示意圖Fig.2 Schematic of computational domain and grid structure

1.2 控制方程

數(shù)值模擬中,通過多相可壓縮求解器compressible-InterIsoFoam 直接求解可壓縮的納維-斯托克斯方程,假設氣相和液相均為可壓縮的互不相溶牛頓流體[25],考慮兩相間的熱量傳輸[26,29],其控制方程為ρ

連續(xù)性方程

動量方程

能量方程

式中,t為時間,ρ為混合物密度,U為速度矢量,p為壓力,τ為應力張量,σ為表面張力系數(shù)為0.07 N/m,κ為界面曲率,δs為特征函數(shù),在界面處網(wǎng)格為1,非界面為0,n為界面單位法向量,e為比內(nèi)能,K為比機械能,β為熱擴散系數(shù),液相熱擴散系數(shù)為1.4 ×10-7m2/s,氣相熱擴散系數(shù)為1.8×10-5m2/s.

考慮到氣相和液相之間會發(fā)生熱交換,其狀態(tài)方程采用理想氣體和理想流體進行建模求解[26,30],表達式分別為

式中,ρg為氣相密度為1.29 kg/m3,ρl為液相密度為1000 kg/m3,R為常數(shù)為8.31 J/(mol·K),T為溫度,pl0為空泡初始內(nèi)壓.

兩相界面通過volume of fluid (VoF) 方法捕捉,其核心是引入流體體積分數(shù)α,定義為對應相體積與網(wǎng)格體積的比值(氣相為0,液相為1),通過求解相分數(shù)方程求解實現(xiàn)兩相界面捕捉[31-32].相分數(shù)方程定義為

式中,α為體積分數(shù),Ur為界面處兩相之間的相對速度,作用于人工壓縮確保界面尖銳.基于流體體積分數(shù),兩相系統(tǒng)中的密度和黏度定義為

式中,μ為液相的黏度,具體取值在后文討論,μg為氣相的黏度為1.0×10-5Pa·s.

考慮到不穩(wěn)定性發(fā)展對于精細界面捕捉的要求,本文采用了isoAdvector 的幾何VoF 方法.該方法采用了有效等值面的思想,對等值面處網(wǎng)格點切割并重新進行跨面連接,形成單元內(nèi)部新的等值面.在迭代過程中,該方法假設界面在子時間步長的間隔期內(nèi)穩(wěn)定移動,來分析計算通過網(wǎng)格表面的體積分數(shù)通量,過程中沒有對單元形狀作出假設,因此可以適用于任意形狀的網(wǎng)格[33].

求解過程采用PIMPLE 算法[34]求解瞬態(tài)的可壓縮納維-斯托克斯方程,采用自適應時間步長基于庫朗數(shù)進行調(diào)整,最大庫朗數(shù)設定為0.5.每個時間步分為外循環(huán)和內(nèi)循環(huán),外循環(huán)首先求解相分數(shù)方程,然后順序求解連續(xù)性方程、動量方程和能量方程,內(nèi)循環(huán)是壓力修正過程,利用溫度通過狀態(tài)方程更新各相密度和混合物密度實現(xiàn)壓力修正,達到求階精度后結(jié)束循環(huán),并步進到下一時間步,直至達到預定求解時間步停止求解.求解過程時間采用Euler 格式進行離散,相分數(shù)方程對流項采用GaussvanLeer 格式進行離散,動量方程對流項采用Gauss線性格式進行離散,其余項均為Gauss 迎風格式.

1.3 數(shù)值模擬驗證

為了驗證數(shù)值模擬的準確性,本小節(jié)討論了數(shù)值模擬結(jié)果的網(wǎng)格無關(guān)性,并與參考文獻[26]中的實驗結(jié)果進行對比.根據(jù)實驗的空泡最大半徑和潰滅時間,確定數(shù)值模擬空泡的初始條件為Rb0=0.2 mm,Pb0=45 MPa,Tb0=1500 K;液滴的初始條件與實驗一致,初始液滴半徑Rd0=6 mm,初始溫度Td0=300 K.網(wǎng)格無關(guān)性選取了粗、中、細3 套網(wǎng)格,最小尺寸分別為40,20 和10 μm,具體參數(shù)見表1.表1 給出了3 套網(wǎng)格中數(shù)值模擬空泡最大半徑Rmax的結(jié)果,并與實驗進行了對比,同時給出不同網(wǎng)格對擾動波數(shù)n和最大擾動幅值ηmax的計算結(jié)果,最大擾動幅值為擾動最大位置(Rηmax(t))與最小值位置(Rηmin(t))的差:ηmax=0.5(Rηmax(t)+Rηmin(t)),結(jié)果顯示空泡最大半徑、擾動波數(shù)和最大擾動幅值誤差逐級減小.圖3(a)展示了液滴半徑(Rd)、空泡半徑(Rb)及擾動幅值(η)變化曲線對比結(jié)果,液滴半徑通過當前時刻擾動最大位置(Rηmax(t))和最小值位置(Rηmin(t))的平均值計算得到Rd(t)=0.5(Rηmax(t)+Rηmin(t)).其中圓形和方形標記分別代表實驗的空泡半徑和液滴半徑,線條為數(shù)值模擬結(jié)果,藍色虛線、橙色實線和綠色點線分別對應了粗網(wǎng)格、中網(wǎng)格和細網(wǎng)格.結(jié)果表明3 套網(wǎng)格計算結(jié)果逐步收斂,粗、中和細網(wǎng)格空泡最大半徑與實驗結(jié)果的誤差分別為0.56%,0.43%和0.37%.因此采用中網(wǎng)格進行計算是合適并且準確的.

表1 粗中細網(wǎng)格參數(shù)細節(jié)Table 1 Details for coarse,medium and fine grids

圖3 實驗與數(shù)值模擬結(jié)果對比Fig.3 Comparisons of experiment and numerical simulations

圖3(b)展示了4 個典型時刻液滴和空泡形態(tài)的實驗與數(shù)值模擬對比結(jié)果,其中左半為數(shù)值模擬結(jié)果,右半為實驗結(jié)果,氣液界面為α=0.5 等值線.過程中,t=0.1～ 0.45 ms 為空泡的膨脹階段,t=0.45～0.85 ms 為收縮階段.從t=0.1 ms 開始,高溫高壓的空泡向外膨脹,并帶動液滴膨脹,此時液滴界面仍然保持光滑.隨后在t=0.3 ms 液滴界面觀察到擾動的產(chǎn)生(如黃色方框局部放大圖所示),此時不穩(wěn)定性開始出現(xiàn),隨后空泡繼續(xù)膨脹在t=0.45 ms 達到最大半徑,液滴界面逐漸明顯.收縮階段,隨著空泡體積迅速收縮,液滴界面擾動快速增長,在t=0.7 ms時能觀察到擾動明顯變大,該階段內(nèi)不穩(wěn)定性迅速發(fā)展.整個發(fā)展過程實驗與數(shù)值模擬對比良好,本文建立的數(shù)值模擬方法可以準確描述空泡振蕩和液滴界面擾動演化過程.下文基于該空泡初始條件,改變液體黏性(μ)和液滴的初始半徑(Rd0)討論液體黏性對柱狀液滴界面不穩(wěn)定性的影響.

2 結(jié)果與討論

柱狀液滴內(nèi)空泡潰滅誘導不穩(wěn)定性產(chǎn)生,隨著不穩(wěn)定性的發(fā)展液滴會展現(xiàn)出穩(wěn)定、通氣和破碎3 種不同的特征形態(tài),在液體黏性影響下空泡和液滴的行為都會發(fā)生改變.

2.1 液滴破碎形態(tài)

圖4 展示了液滴初始半徑Rd0=1.5 mm 不同液體黏性條件下液滴與空泡的行為,其中圖4(a)～圖4(c) 分別對應液體黏性為μ=0.001,0.05,0.1 Pa·s,該工況下液滴表現(xiàn)為破碎特征形態(tài).圖4(a)觀察可得,隨著空泡膨脹,液滴界面開始向外運動,在膨脹初期t=0.04 ms 觀察到液滴表面初始擾動在產(chǎn)生.隨后擾動繼續(xù)增長,在空泡和液滴膨脹到最大半徑前與外界大氣形成通路.此時處于膨脹階段,擾動穿透后空泡和液滴由于慣性繼續(xù)向外運動,液滴在擾動增長和空泡膨脹的共同作用下發(fā)生破碎(t=0.15 ms),不再保持完整形態(tài).這種在膨脹階段擾動穿透空泡界面的情況被稱為液滴的破碎形態(tài).

圖4 液體黏性對Rd0=1.5 mm 破碎形態(tài)液滴與空泡形態(tài)的影響Fig.4 Influence of liquid viscosity on the droplet and bubble at a splashing state at Rd0=1.5 mm

圖4(b)展示了液體黏性變大后的液滴和空泡行為.當液體黏性μ=0.05 Pa·s 時,擾動在t=0.05 ms時產(chǎn)生,且擾動幅值變小.擾動穿透空泡界面仍處于膨脹階段,隨后液滴發(fā)生破碎.

2.2 液滴通氣形態(tài)與通氣發(fā)生時間

圖5 展示了液滴初始半徑Rd0=2.5 mm 不同液體黏性條件下液滴與空泡的行為,其中圖5(a)～圖5(c) 分別對應液體黏性為μ=0.001,0.05,0.1 Pa·s,該工況下液滴初始半徑相較于圖4 變大,液滴表現(xiàn)為通氣特征形態(tài).圖5(a)觀察可得,膨脹階段空泡帶動液滴向外膨脹,在t=0.2 ms 擾動開始產(chǎn)生,其相比于破碎形態(tài)(圖4)更晚,并且此時空泡界面與液滴界面距離更大.空泡在t=0.25 ms膨脹到最大半徑,隨后進入收縮階段.在收縮階段擾動幅值繼續(xù)增大,擾動穿透空泡界面并在空泡和外部大氣形成通路(t=0.35 ms),定義該時刻為通氣發(fā)生時間(tv),此時空泡內(nèi)部壓力小于外部大氣壓,氣體從外部進入內(nèi)部發(fā)生通氣.這種在收縮階段擾動穿透空泡界面的情況被稱為液滴通氣形態(tài).

圖5 液體黏性對Rd0=2.5 mm 通氣形態(tài)液滴與空泡形態(tài)的影響Fig.5 Influence of liquid viscosity on the droplet and bubble at a ventilating state at Rd0=2.5 mm

圖5(b)和圖5(c)展示了液體黏性變大后的液滴和空泡行為.液體黏性增大后擾動產(chǎn)生時間延后,μ=0.05,0.1 Pa·s 分別在t=0.22,0.23 ms 觀察到擾動的產(chǎn)生.當空泡膨脹到最大半徑,擾動幅值明顯小于μ=0.001 Pa·s 黏性工況.隨著擾動發(fā)展其通氣發(fā)生時間也延后,分別是t=0.38,0.4 ms.隨著液體黏性增大,抑制了擾動的發(fā)展,使得需要更長時間才能穿透空泡界面,形成通氣形態(tài).

在液滴通氣形態(tài)下,提取了空泡半徑和擾動幅值作為關(guān)鍵特征量開展分析.圖6(a)展示了液滴通氣特征形態(tài)下(Rd0=2.5 mm,μ=0.001 Pa·s),液體黏性對空泡半徑的影響.觀察發(fā)現(xiàn),隨著液體黏性增大空泡在膨脹階段行為幾乎相同,這是由于膨脹階段空泡膨脹會擠壓液滴,在液滴內(nèi)部形成高壓一起向外膨脹,此時液體黏性并不起主要作用,液體黏性的影響主要體現(xiàn)在收縮階段.在收縮階段,Rd0=2.5 mm,μ=0.001 Pa·s 條件下空泡的最小半徑為0.69 mm,當液體黏性增大至μ=0.05,0.1 Pa·s 時最小半徑分別增大了179%和238%.收縮階段,空泡內(nèi)部低壓會在外部大氣壓的作用下收縮,相比與高溫高壓的初始空泡此時大氣壓的壓差并不大,液體黏性作用逐漸顯著,形成一個與壓差方向相反的力,因此空泡在液體黏性的作用下會更難收縮.液體黏性增加使空泡收縮減慢,最小半徑變大.

圖6 不同液體黏性下,Rd0=2.5 mm 液滴通氣形態(tài)空泡半徑與擾動幅值發(fā)展Fig.6 Temporal variations of bubble radius and perturbation amplitude for different liquid viscosity at a ventilating state at Rd0=2.5 mm

圖6(b) 中展示了液滴通氣特征形態(tài)下(Rd0=2.5 mm,μ=0.001 Pa·s),液體黏性對擾動幅值發(fā)展歷程的影響,定義擾動幅值為η(t),其通過當前時刻液滴界面擾動的最大位置和最小位置計算得到

觀察發(fā)現(xiàn),擾動幅值隨著發(fā)展過程逐漸增大,但在膨脹階段和收縮階段行為有所差異.在膨脹階段(t＜ 0.25 ms)擾動幅值較小并且隨著空泡膨脹緩慢增長,明顯受到液體黏性抑制,液體黏性較小時(μ=0.001 Pa·s)更早觀察到擾動.收縮階段(t＞ 0.25 ms)開始后,擾動隨著空泡收縮快速增長,這是由于空泡收縮帶動液滴內(nèi)界面向內(nèi)運動造成的.當擾動增長至穿透空泡界面,液滴發(fā)展到通氣形態(tài).隨著液體黏性增加擾動的整體幅值均減小,但整體趨勢仍保持一致.因此可得,液體黏性增大擾動幅值會逐漸變小,但并不改變擾動的發(fā)展趨勢.

圖7 展示了液體黏性對不同液滴初始半徑通氣發(fā)生時間(tv)的影響.由圖7 可知,液滴初始半徑相同時,通氣發(fā)生時間隨著液體黏性增大逐漸增加,對于Rd0=2.5 mm 的工況從μ=0.001～ 0.1 Pa·s 增加了15.9%.這是由于液體黏性增大擾動幅值變小,其需要更長時間的增大以形成通氣形態(tài).相比之下,液滴初始半徑的增加也會增加通氣開始的時間,并且其效果比液體黏性的更加明顯.對于μ=0.1 Pa·s 和Rd0=1.8 mm 的工況,發(fā)生通氣僅需要0.22 ms 而Rd0=2.5 mm 則需要0.39 ms,增加了77.27%.液滴初始半徑增大達到通氣形態(tài)需要的擾動幅值更大,因此擾動需要更長的發(fā)展時間,對于繼續(xù)加大液滴初始半徑則液滴形態(tài)會變成穩(wěn)定,此時擾動在最大值時無法發(fā)生通氣.

圖7 液體黏性對不同液滴初始半徑通氣發(fā)生時間的影響Fig.7 Influence of liquid viscosity on the start time of ventilation for different droplet initial radius

2.3 液滴穩(wěn)定形態(tài)與擾動幅值

圖8 展示了液滴初始半徑Rd0=6 mm 不同液體黏性條件下液滴與空泡的行為,其中圖8(a)～ 圖8(c)分別對應液體黏性為μ=0.001,0.05,0.1 Pa·s,此時液滴初始半徑更大,該工況下液滴表現(xiàn)為穩(wěn)定特征形態(tài).圖8(a)觀察可得,隨著空泡膨脹,液滴界面開始向外運動,在t=0.3 ms 觀察到液滴表面初始擾動在產(chǎn)生,其擾動產(chǎn)生時間相比于通氣形態(tài)(圖5)更晚,并且此時空泡界面與液滴界面距離更大.隨后在t=0.45 ms 時膨脹到最大半徑,同時擾動逐漸變大.在空泡收縮過程中液滴界面擾動繼續(xù)增長,在t=0.85 ms 空泡收縮到最小半徑,此時擾動并未穿透空泡界面,空泡內(nèi)部與外界氣體并未形成通路,該液滴形態(tài)成為穩(wěn)定特征形態(tài).

圖8 液體黏性對Rd0=6 mm 穩(wěn)定形態(tài)液滴與空泡形態(tài)的影響Fig.8 Influence of liquid viscosity on the droplet and bubble at a stable state at Rd0=6 mm

圖8(b)和圖8(c)展示了液體黏性變大后的液滴和空泡行為.隨著液體黏性增加,出現(xiàn)擾動時間同樣延后,μ=0.05,0.1 Pa·s 分別在t=0.35,0.4 ms 觀察到擾動的產(chǎn)生.空泡在t=0.45 ms 膨脹到最大半徑隨后開始收縮,液體黏性增大的條件下,收縮階段擾動幅值的增加明顯變緩,在t=0.85 ms觀察到擾動的幅值隨著液體黏性增加逐漸變小,空泡最小半徑則變大.液滴黏性增加,擾動幅值變小,液滴保持穩(wěn)定形態(tài).

在液滴穩(wěn)定形態(tài)下,同樣對關(guān)鍵特征量開展分析.圖9(a) 中展示了液滴穩(wěn)定特征形態(tài)下(Rd0=6 mm,μ=0.001 Pa·s),液體黏性對空泡半徑變化的影響.與通氣形態(tài)下(圖6(a))空泡行為相似,液體黏性的影響主要發(fā)生在收縮階段.液體黏性較小時(μ=0.001 Pa·s)空泡收縮會更快并且半徑會更小,在收縮到最小半徑時為0.56 mm.當液體黏性增大到μ=0.05,0.1 Pa·s 時,其最小半徑變?yōu)榱?.1 和1.9 mm,分別增大了96%和239%.相比于液滴通氣形態(tài),該工況下空泡周期更長,這是由于液滴變大,空泡膨脹需要驅(qū)動水的質(zhì)量增加導致膨脹和收縮都變得更慢.

圖9 不同液體黏性下,Rd0=6 mm 液滴穩(wěn)定形態(tài)空泡半徑與擾動幅值發(fā)展Fig.9 Temporal variations of bubble radius and perturbation amplitude for different liquid viscosity at a stable state at Rd0=6 mm

圖9(b) 中展示了液滴穩(wěn)定特征形態(tài)下(Rd0=6 mm,μ=0.001 Pa·s),液體黏性對擾動幅值發(fā)展歷程的影響.在膨脹階段(t＜ 0.45 ms)液體黏性較小時(μ=0.001 Pa·s)更早觀察到擾動,此后擾動幅值隨著空泡膨脹緩慢增加.收縮階段(t＞ 0.45 ms)開始后,擾動隨著空泡收縮快速增長,并在t=0.85 ms 到達最大值.與通氣形態(tài)下(圖6(b))液滴和空泡行為不同,該條件下液滴初始半徑更大不再發(fā)生通氣,擾動幅值會隨著空泡第2 次膨脹而減小,從而出現(xiàn)先增加后減小的趨勢.此時擾動存在最大幅值,在液體黏性小時(μ=0.001 Pa·s)可達到1.85 mm.

圖10 展示了液滴穩(wěn)定特征形態(tài)下,液體黏性對不同液滴初始半徑最大擾動幅值(ηmax)的影響.由圖可知,當液滴初始半徑相同時,擾動最大值隨著液體黏性增大呈現(xiàn)出逐漸減小的趨勢.對于Rd0=4 mm 的工況,μ=0.001 Pa·s 時ηmax=5.36 mm 而在μ=0.1 Pa·s 時僅為ηmax=1.85 mm,減小了66.2%.這也是由于液體黏性增大后空泡收縮變慢最小半徑變大,更緩慢的變化導致了不穩(wěn)定性擾動增長變慢,最大值變小.對于不同液滴初始半徑其基本趨勢保持一致,但隨著液滴初始半徑增大,其擾動幅值也減小,對于μ=0.001 Pa·s 的工況,Rd0=7 mm 時僅為ηmax=2.23 mm,減小了58.4%.這是由于液滴初始半徑增大后空泡膨脹和收縮變得更加緩慢,因此擾動的增長會減慢.相比之下,從Rd0=4 mm 增大到Rd0=5 mm 最大擾動幅值平均減少了29.3%,而從Rd0=6 mm 增大到Rd0=7 mm 最大擾動幅值平均僅減少了15.0%,證明該效果在液滴初始半徑較小時更明顯.

圖10 液滴穩(wěn)定形態(tài)下,液體黏性對不同液滴初始半徑最大擾動幅值的影響Fig.10 Influence of liquid viscosity on maximum perturbation amplitude for different droplet initial radius in a stable state

2.4 液體黏性影響下的液滴界面不穩(wěn)定性相圖

基于液體黏性和液滴初始半徑對液滴內(nèi)空泡誘導不穩(wěn)定性形態(tài)的影響,本小節(jié)在μ∈[0.001,100]Pa·s,Rd0∈[1.5,7] mm 兩個參數(shù)區(qū)間上開展了一系列的數(shù)值模擬,并歸納得到液滴不穩(wěn)定性相圖,如圖11 所示,其中藍色為破碎形態(tài),橙色為通氣形態(tài),綠色為穩(wěn)定形態(tài),圖11 中藍色圓形、橙色方形和綠色三角形為數(shù)值模擬得到的結(jié)果.由圖11 可得,在本文的空泡參數(shù)下,小液滴(Rd0＜ 2 mm)條件下,液體黏性較小時(μ=0.001 Pa·s)通氣發(fā)生在膨脹階段,此時液滴為破碎形態(tài).隨著液體黏性增加,擾動幅值減小,通氣發(fā)生時間延后.在μ＞ 0.1 Pa·s 時,空泡與外界大氣連通延后至收縮階段,液滴形態(tài)從破碎變?yōu)橥?如μ=0.1,1 Pa·s).隨著液體黏性持續(xù)增大(μ=10 Pa·s),擾動減小至不能穿透空泡界面,液滴變?yōu)榉€(wěn)定形態(tài).Rd0=1.5 mm 的典型液體黏性的液滴形態(tài)在下方黑色方塊中展示,隨著液體黏性增加,液滴形態(tài)從破碎轉(zhuǎn)變?yōu)橥庠俎D(zhuǎn)變?yōu)榉€(wěn)定.中液滴(2 mm ＜Rd0＜ 3 mm)的特征形態(tài)受液體黏性影響,在液體黏性較小時(μ=0.001,0.1 Pa·s)擾動隨著發(fā)展穿透空泡界面,液滴為通氣形態(tài),隨著液體黏性增大最大擾動幅值減小至小于空泡界面和液滴界面的距離,在μ＞ 2 Pa·s 時不再通氣,液滴從通氣轉(zhuǎn)變?yōu)榉€(wěn)定形態(tài)(μ=1,10 Pa·s).Rd0=2.5 mm 的典型液體黏性的液滴形態(tài)在下方紫色方塊中展示,液體黏性增加液滴形態(tài)從通氣轉(zhuǎn)變?yōu)榉€(wěn)定.隨著液滴初始半徑變大液滴變得更穩(wěn)定,大液滴(Rd0＞ 3 mm)的特征形態(tài)不再隨液體黏性發(fā)生變化,液體黏性增加最大擾動幅值減小,液滴維持穩(wěn)定形態(tài).綜上,液滴初始半徑增大和液體黏性增大均能使液滴更加穩(wěn)定,液滴會從破碎變?yōu)橥?進而變成穩(wěn)定形態(tài).

圖11 不同初始半徑與黏度下的液滴不穩(wěn)定性相圖Fig.11 Phase diagram of the interfacial instability of droplets with different initial radius and liquid viscosity

3 結(jié)論

在內(nèi)空泡振蕩下液滴會發(fā)生破碎、通氣及保持穩(wěn)定3 種特征形態(tài),本文基于開源OpenFOAM 框架,采用直接數(shù)值模擬和isoAdvector 幾何VoF 方法對液滴界面不穩(wěn)定性開展研究,討論了不同液滴初始半徑下液體黏性對液滴界面擾動、通氣時間等關(guān)鍵特征量的影響規(guī)律,總結(jié)了液體黏性與液滴初始半徑參數(shù)下的液滴不穩(wěn)定性相圖.

(1)液體黏性同時影響空泡振蕩和液滴界面擾動演化.對于空泡發(fā)展,液體黏性主要影響空泡的收縮運動,液體黏性的增加會使空泡收縮變得更慢,且最小半徑會更大;對于擾動發(fā)展,液體黏性會抑制擾動前期增長,隨著液體黏性增加擾動幅值減少,但擾動的發(fā)展趨勢不發(fā)生改變.

(2)液體黏性和液滴初始半徑共同影響了該過程的關(guān)鍵特征量趨勢.對于通氣形態(tài),隨空泡發(fā)展擾動幅值增大至穿透空泡界面發(fā)生通氣,液體黏性或液滴初始半徑增加,通氣發(fā)生時間增加;對于穩(wěn)定形態(tài),擾動幅值呈現(xiàn)先增加后減小趨勢,隨著液體黏性或液滴初始半徑增加,最大擾動幅值減小.

(3)液體黏性和液滴初始半徑顯著改變柱狀液滴的特征形態(tài).在文章討論的參數(shù)范圍內(nèi),隨著液滴初始半徑增加,液滴形態(tài)從破碎先轉(zhuǎn)變?yōu)橥庠俎D(zhuǎn)變?yōu)榉€(wěn)定.隨著黏性增加,小液滴(Rd0＜ 2 mm)的形態(tài)從破碎轉(zhuǎn)變?yōu)橥膺M而變成穩(wěn)定;中液滴(2 mm ＜Rd0＜ 3 mm)的形態(tài)從通氣轉(zhuǎn)變?yōu)榉€(wěn)定,不出現(xiàn)破碎形態(tài);大液滴(Rd0＞ 3 mm)的形態(tài)不隨液體黏性改變,保持穩(wěn)定形態(tài).液滴形態(tài)會隨著液體黏性或初始半徑增加趨于穩(wěn)定.

本文重點分析了液滴初始半徑和液體黏性對液滴界面不穩(wěn)定性發(fā)展的影響,在今后研究中,計劃基于穩(wěn)定性分析方法從理論上開展系統(tǒng)參數(shù)研究,推導得出液體黏性作用下的擾動模數(shù)及其增長特性.