油滴撞擊油膜層內氣泡的變形與破裂過程的數值模擬?

周劍宏 童寶宏? 王偉 蘇家磊

1 引 言

隨著工業化水平的提高,人們對于機械產品的安全性、可靠性、節能性以及服役壽命等提出了更高的要求.據統計,在工業生產中約50%的設備損壞是由磨損引起的.潤滑作為一種減小摩擦和降低磨損的有效手段,普遍應用于工業生產中.在對旋轉工作的機械零部件和設備進行潤滑時將不可避免地涉及油滴撞擊油膜的現象.例如,在對航空發動機的軸承進行噴油潤滑時,為防止潤滑油泄漏引起焦化和著火等現象,常采用空氣系統對軸承室進行密封,使軸承室內部處于復雜的油-氣混合狀態,具體表現為油滴和油膜間的相互作用行為[1];在對齒輪、渦輪等減速器裝置進行飛濺潤滑時,機箱中旋轉的零部件將潤滑油霧化成油滴帶到摩擦副上形成自動潤滑,部分油滴將會與機箱中的油膜發生碰撞行為;在對高速電主軸進行油-氣潤滑時,零部件表面受連續油滴撞擊并形成潤滑油膜層,隨時間變化,后續油滴常與沉積的油膜層發生作用.

目前,在學術期刊上較常見的是非油類液滴與液膜的碰撞行為研究,其研究內容主要分為兩個方面:一方面是關注液滴撞擊液膜的流動行為;另一方面則集中考察氣泡夾帶現象.

針對液滴撞擊液膜的流動行為已開展了大量的研究工作,且具有較高的成熟度[2?11].Rioboo等[2]和Okawa等[3]實驗研究了液滴撞擊液膜的過程,并考察了液滴撞擊速度、液滴直徑和液膜厚度對撞擊結果的影響.研究結果表明,在不同的條件下,液滴撞擊液膜后會出現3種不同的運動形態:鋪展、皇冠水花和飛濺.郭加宏等[4]對液滴撞擊液膜后的水花流動行為進行了可視化實驗研究,發現液體黏性、表面張力、液膜厚度和液滴直徑等因素均對水花流動有一定的影響.宋云超等[5]通過數值模擬對液滴撞擊濕潤壁面的現象進行了研究,發現隨著液滴撞擊速度的增大,會依次出現4種不同的運動形態:黏附鋪展、波動運動、皇冠幾何體運動以及飛濺運動.梁剛濤等[6]采用耦合水平集-體積分數(CLSVOF)方法數值模擬了液滴撞擊液膜的運動過程,發現了液滴頸部射流的產生機理,并研究了黏度和表面張力影響水花形態轉變的規律.

隨著液膜厚度的增加,可以忽略壁面形貌特征對液滴撞擊過程的影響.Bisighini[7]提出液滴撞擊深液膜(液池)后會出現融合、彈坑狀水花和飛濺3種不同的現象.Rein[8]實驗研究了韋伯數(We)對聚結和飛濺現象之間的過渡過程產生影響的規律.Blanchette等[9]通過數值模擬研究了低速液滴撞擊深液膜的過程,發現在某些情況下,液滴不會與液膜完全融合,而會產生一個子液滴,并將該現象稱為部分融合.Chen等[10]對液滴撞擊液膜中的部分融合現象進行了可視化實驗,認為這些現象主要受慣性和界面張力的影響.Ray等[11]采用CLSVOF方法數值研究了液滴撞擊液膜過程中由部分融合到完全融合的變化規律.

當液滴撞擊液膜形成彈坑狀水花時,容易形成氣泡夾帶現象,這一現象已引起學者們的廣泛關注.在海洋和氣候科學中,氣泡夾帶現象有利于氣液界面氣體進行分子交換;而在液體食品和藥劑加工過程中,氣泡夾帶現象則會導致泡沫產生,嚴重影響產品質量和生產速率.由此可見,確定夾帶發生的條件以及確定這種現象在數量上的發生概率是極為重要的.

相關研究顯示,液滴撞擊液池后會夾帶大小不同的氣泡[12?21].通常最小的氣泡會溶解在液池內,最大的氣泡會上升到自由表面并發生破裂[12].Pumphrey等[13]提出了4種氣泡夾帶形式:不規則夾帶、規則夾帶、大氣泡夾帶和“Mesler夾帶”.其中,液滴撞擊液面后形成彈坑狀水花是發生規則夾帶和大氣泡夾帶的前提.

規則夾帶是指當液滴直徑和速度在一定范圍時,總能觀察到單個氣泡的夾帶現象.Pumphrey等[14]對規則夾帶的形成機制進行了研究,認為毛細波對夾帶的形成起著至關重要的作用.Oguz等[15]考察了We數和弗勞德數(Fr)對規則夾帶的影響,并以We=41.3Fr0.179和We=48.3Fr0.247作為臨界曲線,在這兩條曲線間的區域總能觀察到規則夾帶現象.

Zou等[16]和Wang等[17]對大氣泡夾帶現象進行了可視化實驗研究,發現扁長形液滴撞擊液面后能形成與液滴尺寸幾乎相當的氣泡.Deng等[18]研究了液體黏度和表面張力對氣泡夾帶的影響,發現夾帶氣泡尺寸隨著毛細數(Ca)的增大而減小.

Thoroddsen等[19]在實驗中發現低速液滴撞擊液池時,液滴與液池間存在半球狀的空氣層阻延兩者直接接觸.隨著時間變化,空氣層發生破裂并分裂為成百上千個小氣泡.Sigler和Mesler[20]較早地研究了這種現象,稱其為“Mesler夾帶”.Saylor和Bounds[21]通過研究發現硅油液滴比水滴更易觀察到“Mesler夾帶”現象.

通過上述文獻可以看到,學者們主要對液滴撞擊液膜的流動行為以及氣泡夾帶現象進行了較為系統且深入的研究,并對各現象的形成機制和不同物理參數的影響也進行了相應探討.

在潤滑系統工作過程中,油滴撞擊油膜時易發生氣泡夾帶現象.氣泡加速潤滑油氧化變質,使油膜出現破裂,造成機械部件的直接接觸并最終影響設備的服役壽命.此外,油滴往往是連續撞擊油膜,當第一個油滴撞擊油膜引起氣泡夾帶后,第二個油滴撞擊的是含有氣泡的油膜.因此,研究油滴與含氣泡油膜的碰撞行為可以為潤滑過程中氣泡影響研究提供參考,同時對于深化氣泡動力學規律的理解也具有重要意義.

當前較少見到公開發表的文獻中有關于液滴與含氣泡液膜碰撞行為的研究,在缺乏相關理論指導和實驗研究的前提下,數值模擬研究是揭示其運動規律的有效途徑.基于上述背景,本文采用CLSVOF方法模擬了油滴撞擊含氣泡油膜的運動過程,重點探討油膜層中氣泡在油滴撞擊作用下的變化規律.

2 數值模型計算方法

2.1 控制方程

油滴撞擊含氣泡油膜是典型的兩相流現象,氣液界面流動狀況十分復雜,且各物理參數間也存在較強的相互作用,求解這類問題的關鍵在于如何精確追蹤相界面.目前,流體體積法(volume of fl uid,VOF)和水平集(level set)方法廣泛應用于追蹤相界面的數值模擬中.

VOF方法通過定義流體體積函數β來捕捉相界面,當β=1時,目標區域僅存在液體;當β=0時,目標區域內僅存在氣體;當0＜β＜1時,目標區域為兩相界面.流體體積函數β的控制方程為

式中U為速度矢量,t為實際時間.

VOF方法在計算兩相流動時能保證質量守恒,但流體體積函數β不是連續函數,不能精確計算曲率和法向向量等幾何參數[22].

Level set方法通過光滑連續的水平集函數?(x,t)來追蹤相界面,能準確計算界面上的幾何參數.通過水平集函數?(x,t)可得法向向量n和界面曲率κ(?)為:

Level set方法的缺點在于求解兩相流動時容易造成質量不守恒,較難追蹤尖銳的界面變化.為了改善VOF和level set方法的缺陷,Sussman等[22]對這兩種方法進行了耦合并提出了CLSVOF方法.

近年來,CLSVOF方法廣泛應用于模擬液滴撞擊壁面現象[5,6,11,23]、液體中氣泡變化行為[24?26]等方面.在CLSVOF方法中,水平集函數?(x,t)用來精確求解界面上的幾何參數,流體體積函數β則用于計算流體體積分數.

通過CLSVOF方法可以得到氣-液界面處的密度ρ(?)和黏度μ(?)為:

式中 H(?)是Heaviside函數,下標g和l分別表示氣相和液相.

式中ε=1.5?r,?r為最小的網格尺寸.

表面張力的計算通過連續表面力(continuum surface force)模型[27]

來計算,式中F為表面張力,σ是表面張力系數.

兩相流動的質量方程和動量方程分別表示成

式中g為重力加速度,P為壓力.

2.2 計算方法可行性驗證

當前尚未見到有關液滴撞擊含氣泡液膜的實驗研究,為驗證CLSVOF方法是否能準確模擬油滴撞擊含氣泡油膜的變化過程,我們通過間接驗證的方法對液滴撞擊液膜過程中的流動行為和大氣泡夾帶現象進行試算.

數值計算過程中控制方程通過有限體積法進行離散,壓力速度耦合采用PISO方法,壓力求解使用PRESTO!方法,氣-液界面插值選擇幾何重構方法,動量求解采用二階迎風格式.圖1(a)為扁長形液滴撞擊液池后形成大氣泡夾帶現象的二維數值模擬結果與文獻[17]中實驗結果的定性對比.圖中扁長形液滴的等效直徑為5.65 mm,碰撞速度為0.953 m/s,We=75,Fr=16.由圖1(a)可以看出,計算結果與實驗觀察符合較好.此外,數值計算采用二維軸對稱模型,模擬結果能發現實驗中無法觀察的流場內部細節特征.

圖1(b)為液滴撞擊液膜后水花內徑Din隨實際時間t變化的模擬結果與文獻[28]中實驗結果的定量對比,其中液滴的雷諾數(Re)和We數分別為13676和667,無量綱液膜厚度δ可表示為

式中H為液膜厚度,D為液滴直徑.

圖1 計算結果與實驗結果對比 (a)定性對比;(b)定量對比Fig.1. Comparison of experiment and simulation:(a)Qualitative comparison;(b)quantitative comparison.

從圖1(b)中可以看到,撞擊前期模擬和實驗中水花內徑Din較為接近,均隨時間近似遞增;但在撞擊后期,模擬結果與實驗結果有所差異.分析認為其原因是模擬條件與實驗條件之間存在一定差異,此外,現有數字圖像處理技術也存在局限性.

綜上所述,CLSVOF方法既能準確捕捉液滴撞擊液膜過程中的氣泡形態變化,也能較精確計算出水花流動形態的特征參數變化.因此,本文采用CLSVOF方法模擬油滴撞擊含氣泡油膜這一動態過程.

2.3 參數選擇及數值模型建立

本文以潤滑油液滴作為研究對象,對油滴與含氣泡油膜的碰撞行為進行研究.在不同的潤滑方式中,油滴直徑大小差異明顯.例如,油-氣潤滑系統工作時,油滴顆粒直徑一般為10—100μm;噴油潤滑中的潤滑油直徑多為亞毫米級;油霧潤滑中的潤滑油顆粒直徑約為1—3μm.為方便研究,將模擬過程中的油滴直徑固定為50μm,其余相關參數如表1所列.

表1 計算時采用的相關參數Table 1.Parameters used in calculation.

圖2所示為初始時刻油滴撞擊含氣泡油膜的二維軸對稱計算域示意圖.為使模擬結果具有三維特征,對二維模型的對稱軸x添加軸旋轉對稱(axisymmetric swirl)條件.定義初始時刻油滴和油膜剛開始接觸,油滴中心和氣泡中心在一條直線上,氣泡位于油膜高度中心.油膜底部為不滑移壁面,環境壓力保持一個標準大氣壓.整個計算區域大小為400μm×150μm,并采用四邊形網格對計算區域進行離散.

圖2 計算域示意圖Fig.2.Schematic of computational domain.

2.4 網格無關性分析

為驗證網格的無關性,圖3給出了4種不同密度的計算網格下中心油膜高度h0隨實際時間t的變化情況,圖中氣泡直徑為10μm.為提高計算效率同時保證計算精度,對氣泡及周圍區域網格進行局部加密,圖3中網格尺寸為加密后的網格.

圖3 不同網格下的中心油膜高度對比Fig.3.Comparison of central fi lm height in diff erent grids.

由圖3可得,網格尺寸對中心油膜高度的變化規律有一定的影響.在30μs＜t＜45μs時,不同尺寸的網格所對應的中心油膜高度最小值所在時間不完全相同.隨著網格尺寸減小,網格數量增多,氣-液界面劃分將越精細.從油膜高度最小值所在時間來看,當網格尺寸達到0.0332μm2時,油膜高度最小值所在時間基本不發生變化.因此,選用尺寸為0.0332μm2的網格作為計算網格.

3 結果與分析

3.1 氣泡大小的影響

圖4為油滴撞擊含不同大小氣泡的油膜時氣泡的形態演化過程.圖4(a)、圖4(b)和圖4(c)中氣泡直徑分別為10,20,和25μm,油滴直徑固定為50μm,其余參數見表1.

由圖4(a)可以看出,當運動剛剛開始時,油膜表層受油滴撞擊作用有較高的豎直向下速度,油膜底層靠近壁面處速度較低,可見油膜在豎直方向上具有較大的速度梯度,具體表現為油膜內部豎直方向的黏性剪切力;而在油膜的水平方向上,油滴的部分動能轉換為油膜向兩側運動的動能,使油膜內部水平方向同樣存在著黏性剪切力作用.氣泡在兩種不同方向的黏性剪切力的共同作用下開始下沉并逐漸變為橢球形.當t=6μs時,氣泡的變形程度達到最大,氣泡上部較平,整體呈現半橢球形;在6μs之后,氣泡受黏性剪切力的影響減弱,表面張力開始主導氣泡變形;在表面張力的作用下氣泡逐漸變圓,并于13μs時形成橢球形;隨時間演化,當t=35μs時氣泡整體形狀近似于球形,且氣泡上表面位置接近于液面;當t=37μs時,氣泡發生破裂,位于氣泡上表面的油膜受表面張力作用形成膜液滴.本文將這種運動稱為“自由表面破裂”.

圖4 不同大小氣泡的演化過程 (a)10μm;(b)20μm;(c)25μmFig.4.Evolution process of diff erent sizes of bubbles:(a)10μm;(b)20μm;(c)25μm.

從圖4(b)可以看出,當氣泡直徑為20μm時,油滴撞擊油膜后氣泡變形過程與氣泡直徑為10μm時基本一致,但最終氣泡沒有發生破裂而是穩定存在于油膜內,本文將其稱為“穩定變形”.

從圖4(c)可以看出,當氣泡直徑為25μm時,氣泡變形過程出現明顯變化.當t=1μs時,氣泡受油滴撞擊發生變形,整體類似于碗狀;當t=6μs時,氣泡在黏性剪切力的作用下發生破裂,并在油膜中心形成一個橢球形小氣泡,小氣泡外圍分布著一些更小的氣環,原先氣泡的外側形成一個大氣環;在6μs到20μs間,油膜內的氣泡和氣環在表面張力的作用下逐漸變圓,油膜中心附近的氣泡及氣環位置幾乎保持不變,邊緣處的大氣環隨著油膜向外側移動;在20μs到100μs間,黏滯阻力和表面張力消耗了油膜向兩側傳播的動能,油膜高度開始回復,油膜中心的小氣泡和周圍氣環發生并聚,并形成體積稍大的氣泡,油膜內部外側的大氣環向油膜中心回縮;當t=200μs時,隨著液面升高,兩側的大氣環移動到油膜中心并形成大氣泡,油膜中心的小氣泡位于大氣泡下方.本文將這種運動過程稱為“油膜內部破裂”.

3.2 氣泡大小對氣泡變形特征參數的影響

由3.1節可見,氣泡大小對氣泡變形歷程及變形特征有顯著的影響.為了進一步研究氣泡大小對氣泡變形的影響規律,本節主要對不同直徑氣泡的變形特征進行分析.

為定量描述油滴撞擊含氣泡油膜壁面的運動特性,如圖5所示,定義氣泡的橫向長度為w,縱向長度為h,使用無量綱化參數形狀系數E來表征氣泡的形狀變化,并可表示為

形狀系數E越大,說明氣泡總體變形量越大.同樣將運動變化時間以無量綱時間T表示:

式中V0為油滴的碰撞速度,t為實際時間,d0為油滴直徑.

圖5 氣泡的橫向長度w和縱向長度hFig.5.Transverse length w and longitudinal length h of the bubble.

圖6 給出了氣泡直徑d在10—30μm條件下,油滴撞擊含氣泡油膜后氣泡形狀系數E隨無量綱時間T的變化情況(為了能清晰看出差別,圖中分成d 6 20μm和d>20μm兩種情況,這里不列出氣泡破裂后的數值變化情況).其中,油滴直徑d0為50μm,碰撞速度V0為20 m/s,其余參數見表1.

圖6 氣泡大小對形狀系數的影響 (a)d 6 20μm;(b)d>20μmFig.6.Infl uence of bubble size on shape parameter:(a)d 6 20μm;(b)d>20μm.

由圖6(a)可以看出,油滴在撞擊含氣泡油膜后,氣泡形狀變化過程分為兩個階段:顯著變形階段(T<10)和相對穩定階段(T>10).在顯著變形階段,隨著氣泡直徑增大,同一時間的氣泡形狀系數也相應增大,且氣泡達到形狀系數最大值的時間也增加,但不明顯.這是因為氣泡直徑越大,氣泡的表面積越大,越容易受黏性剪切力拖拽而產生變形.在相對穩定階段,各氣泡間的形狀系數差距較小.此外,氣泡直徑越大,氣泡破裂發生時間越晚,且破裂前的形狀系數呈現小范圍波動變化.

從圖6(b)中可以看出,當d>20μm時,氣泡破裂發生在顯著變形階段.在T<1.8時,隨著氣泡直徑的增大,同一時間的氣泡形狀系數也相應增大,但氣泡形狀系數最大值和氣泡直徑間無明的顯線性關系.此外,氣泡達到形狀系數最大值的時間和破裂發生時間均隨著氣泡直徑的增大而縮短.較為特殊的是,該階段的氣泡發生的是內部破裂并存在多個破裂位置,且各個位置破裂時間不完全相同.為便于分析,將氣泡最先發生破裂的時間點定為破裂發生時間.經過觀察發現,這一階段的氣泡破裂不是發生在氣泡形狀系數達到最大值時,而是發生在氣泡達到形狀系數最大值后的1μs內.為對這一現象進行分析,圖6(b)內插圖中給出了d=30μm,T=1.84和T=1.92時氣泡右半空間分布相圖.從內插圖可以看出,在T=1.84時氣泡形狀系數最大,而在T=1.92時氣泡靠近中心處的小范圍區域內出現收縮,使氣泡縱向距離h增大,引起氣泡形狀系數減小.隨時間變化,氣泡破裂發生在收縮區域內.

行為異常組在矛盾性和獨立性維度的得分明顯高于正常組（P＜0．05），而在親密度、情感表達、知識性、娛樂性、道德宗教觀、組織性及控制性維度的得分則明顯低于正常組（P＜0．05）。見表1。

綜合上述分析,氣泡大小將顯著影響氣泡的變形歷程.d=20μm是氣泡兩種破碎形式間的臨界點,此時氣泡能穩定地存在于油膜中(圖6(a)僅給出了T 6 40時氣泡的形狀系數);當d<20μm時,氣泡破碎發生在自由表面附近并處于相對穩定階段;當d>20μm時,氣泡破裂發生在油膜內部并處于顯著變形階段.

3.3 氣泡位置對氣泡變形特征參數的影響

3.2 節中不同大小的氣泡均位于油膜中心,但各氣泡頂端和油膜上表面的間距以及氣泡底端和壁面的間距不相等.可見氣泡大小的不同間接引起了氣泡位置的變化,因此本節研究氣泡位置對氣泡變形過程的影響.

圖7給出了d=10,20和30μm時,油滴撞擊含3種不同位置(氣泡頂端和油膜表面相切、氣泡位于油膜中心和氣泡底端和壁面相切)氣泡的油膜后氣泡形狀系數E隨無量綱時間T的變化關系.為方便說明,下文將這3個位置簡稱為頂端、中心和底端.此外,圖7中內插圖為氣泡位于油膜頂端時氣泡破裂相圖以及氣泡位于底端時氣泡穩定相圖,位于油膜中心的氣泡變化相圖可參考圖5.計算過程中的其余參數見表1.

圖7 氣泡位置對形狀系數的影響 (a)d=10μm;(b)d=20μm;(c)d=30μmFig.7.Infl uence of bubble position on shape parameter:(a)d=10μm;(b)d=20μm;(c)d=30μm.

從圖7(a)可以看出,當d=10μm時,位于油膜頂端和油膜中心的氣泡形狀系數隨無量綱時間的變化趨勢保持一致,氣泡均在自由表面發生破裂且破裂時間也較為接近.不同點在于T＜6時,位于油膜頂端氣泡的形狀系數始終大于位于油膜中心氣泡的形狀系數.

由圖7(b)可以看出,當d=20μm時,位于油膜頂端和油膜中心的氣泡形狀系數隨無量綱時間的變化趨勢出現顯著不同,其原因在于位于油膜頂端的氣泡在運動過程中出現油膜內部破裂現象,而位于油膜中心的氣泡在撞擊過程中保持較穩定的狀態.

從圖7(c)可以看出,當d=30μm時,位于油膜頂端和油膜中心的氣泡均在運動過程中發生油膜內部破裂現象,但油膜頂端的氣泡最大形狀系數遠大于油膜中心的氣泡最大形狀系數.

圖8 底層氣泡大小對形狀系數的影響Fig.8. Infl uence of bottom bubble size on shape parameter.

綜合圖7和圖8中4幅圖可以得出,氣泡位置對氣泡變形過程有一定的影響.當氣泡直徑相同時,氣泡越接近油膜頂端受到油滴撞擊的影響越強,越容易發生變形;位于油膜底端的氣泡受油滴撞擊影響較弱,同時壁面對氣泡也有一定吸引作用,因此位于油膜底部的氣泡整體變形趨勢較平緩,并最終依附在壁面上維持穩定的狀態.

3.4 氣泡破裂機制分析

氣泡破裂現象屬于氣泡動力學范疇.在進行氣泡動力學分析時常引入兩個無量綱參數,莫頓數(M o)和厄缶數(Eo):

式中g為重力加速度,μl為液體黏度,ρl為液體密度,ρg為氣體密度,σ為表面張力系數,d為氣泡直徑.

M o數反映氣泡外部液體黏性對氣泡變形的影響,Eo數表征氣泡浮力和表面張力的關系.Eo數越大,氣泡越易受浮力作用發生形變.本文中氣泡的M o=5.06,Eo=2.51×10?5—2.26×10?4,可見氣泡變形過程中受液體黏性作用較大,浮力對氣泡變形的影響較為微弱.

3.4.1 自由表面破裂

圖9給出了氣泡直徑d為10μm,位于油膜中心處的氣泡在破裂區域內的壓力云圖和速度矢量圖.其中,圖9(a)左側為36.0μs時氣泡分布相圖,圖9(a)、圖9(b)和圖9(c)為隨時間變化的壓力-速度分布圖.

從圖9(a)可以看出,在36.0μs時,氣泡頂端已經接近油膜表面,壓力最大值位于氣泡上部的薄膜區域,大小為53 kPa,壓力最小值位于氣-液界面的兩處位置,大小為?6 kPa.在圖片右側的速度分布中,油膜上方較遠處的空氣速度方向豎直向下,但在接近油膜層時速度方向發生變化.這是由于相互接觸的空氣和油膜存在相對運動,使油膜表面產生了毛細波[29].毛細波的存在使氣-液界面發生不穩定現象,這正是氣-液界面出現局部負壓的原因.此外,氣泡內表面速度方向主要受表面張力作用呈現回旋狀.

圖9 氣泡破裂區域壓力和速度分布圖(d=10μm) (a)36.0μs;(b)36.2μs;(c)36.4μsFig.9.Distribution of pressure and velocity in the region of bubble rupture(d=10μm):(a)36.0μs;(b)36.2μs;(c)36.4μs.

從圖9(b)可以看出,在36.2μs時,氣泡上表面的油膜形成了“頸部”區域,該區域內部壓力最小值為?13.31 kPa,兩側壓力為33.6 kPa.可見“頸部”區域承受較大壓力梯度作用有發生斷裂的趨勢.由速度分布可知,受毛細波的影響,氣泡上側油膜附近速度方向變化較為劇烈.氣-液相界面的曲率發生變化,造成相界面的表面張力方向發生改變并促進了“頸部”區域的形成.

從圖9(c)可以看出,“頸部”區域發生斷裂,氣泡內部壓力迅速降低,“頸部”區域內的油膜在表面張力的作用下形成膜液滴.從速度分布中可見,由于氣泡發生破裂,氣泡內的速度方向和大小都發生改變,氣泡內的氣體通過破裂口向外界逃逸,破裂口的速度最大值為11.72 m/s.

通過上述分析可知,表面張力是促使氣泡破裂的主要作用力.該階段的氣泡破裂主要是由氣-液界面的不穩定引起的,具體表現為毛細波傳播造成油膜頸部區域壓力大小和速度方向的改變.

3.4.2 油膜內部破裂

圖10給出了氣泡直徑d為25μm,位于油膜中心處的氣泡在破裂區域內的壓力云圖和速度矢量圖.從圖10(a)可以看到,在5.2μs時,氣泡中心處壓力為28.55 kPa,兩側壓力稍小,大小為25.12 kPa.速度分布顯示,氣泡上側和下側的油膜內部速度方向豎直向下,油膜兩側速度方向呈現傾斜向下,且氣泡上側速度大于下側速度.氣泡內部速度方向由中心向兩側,兩側速度較大,為10.2 m/s;氣泡中心處速度僅為2.02 m/s.可以看出,氣泡內部存在較大的速度梯度,該速度梯度是氣泡發生變形的主要原因.

由圖10(b)可以發現,在5.5μs時,隨著氣泡變得扁長,氣-液界面逐漸變得不穩定;氣泡邊界處出現負壓,大小為?34.36 kPa.壓力最大值位于氣泡內部且偏離中心,大小為224.85 kPa.圖中黑色虛線為壓力梯度較大的區域,該區域內氣-液界面附近的速度方向與其他區域內速度方向明顯不同.此外,在紅色實線內的速度局部放大圖中,氣泡較薄處的速度方向由水平方向轉為傾斜向下,大小為1.77 m/s;氣泡較厚處的速度方向仍保持水平方向,大小為11.36 m/s.

從圖10(c)可以看到,在5.6μs時,氣泡發生破裂,破裂位置發生在圖10(b)中壓力梯度較大區域.在速度局部放大圖中可以發現,破裂處左邊氣泡的速度方向呈現兩種模式:一種受外界油膜影響并與油膜速度方向保持一致;另一種是表面張力作用使氣泡有變圓趨勢,速度方向在氣泡邊緣處呈現回旋狀,而破裂處右邊氣泡內部的速度較大并保持水平方向.可以看出,氣泡破裂發生區域具有兩個特征:一是氣泡破裂前,該區域內存在較大的壓力梯度;二是該區域位于速度方向變化的過渡范圍內.

綜合上述分析,這一階段的氣泡在變形過程中主要受油膜內部黏性剪切力和表面張力的作用.其中黏性剪切力作用是造成該類型氣泡破裂的重要因素,具體表現為圖10(a)中氣泡內部較大的水平方向速度梯度,在速度梯度的作用下,氣泡變得扁長.扁長形氣泡的界面不穩定性使氣泡內部出現較大的壓力梯度,同時讓氣泡內部速度方向發生改變,最終破裂.

圖10 氣泡破碎區域壓力和速度分布圖(d=25μm) (a)5.2μs;(b)5.5μs;(c)5.6μsFig.10.Distribution of pressure and velocity in the region of bubble rupture(d=25μm):(a)5.2μs;(b)5.5μs;(c)5.6μs.

4 結 論

采用CLSVOF方法對油滴撞擊含氣泡油膜現象進行了數值模擬分析,研究了油膜層內部氣泡的運動演化過程,考察了氣泡大小及氣泡位置等因素對受撞擊作用的氣泡變形特征參數的影響規律,探索了不同破裂方式下的氣泡破裂機制,主要得到以下結論.

1)通過與實驗結果比較,驗證CLSVOF方法在求解液滴撞擊液膜過程中的氣泡形態變化和流動特征參數變化的可行性,并以此為依據,建立了油滴撞擊含氣泡油膜的數值模型.

2)氣泡大小是決定油滴撞擊油膜后氣泡呈現自由表面破裂、穩定變形以及氣泡在油膜內部破裂等不同運動形態的重要因素.

3)當10μm 6 d<20μm時,隨著氣泡直徑增大,氣泡在運動過程中的變形量也相應增大,并隨時間變化發生自由表面破裂,破裂發生時間隨氣泡直徑增大而延長;當20μm

4)當氣泡位置發生改變時,氣泡的變形過程也不完全相同.相同撞擊條件下,位于油膜頂端的氣泡比位于油膜中心處的氣泡更易產生變形;位于油膜底層的氣泡總體變形量最小,并最終依附在壁面上.

5)氣泡發生自由表面破裂和油膜內部破裂是氣-液界面的不穩定和表面張力共同作用的結果.油膜內部破裂時,除了受氣-液界面不穩定性和表面張力的作用外,黏性剪切力的影響也十分重要.

[1]Peduto D 2015 Ph.D.Dissertation(Karlsruhe:Karlsruhe Institute of Technology)

[2]Rioboo R,Bauthier C,Conti J,VouéM,de Coninck J 2003 Exp.Fluids 35 648

[3]Okawa T,Shiraishi T,Mori T 2006 Exp.Fluids 41 965

[4]Guo J H,Dai S Q,Dai Q 2010 Acta Phys.Sin.59 2601(in Chinese)[郭加宏,戴世強,代欽 2010物理學報 59 2601]

[5]Song Y C,Ning Z,Sun C H,Lyu M,Yan K,Fu J 2013 Chin.J.Theor.Appl.Mech.45 833(in Chinese)[宋云超,寧智,孫春華,呂明,閻凱,付娟 2013力學學報45 833]

[6]Liang G T,Guo Y L,Shen SQ 2013 Acta Phys.Sin.62 024705(in Chinese)[梁剛濤,郭亞麗,沈勝強2013物理學報62 024705]

[7]Bisighini A 2010 Ph.D.Dissertation(Bergamo:University of Bergamo)

[8]Rein M 1996 J.Fluid Mech.306 145

[9]Blanchette F,Bigioni T P 2009 J.Fluid Mech.620 333

[10]Chen X,Mandre S,Feng J J 2006 Phys.Fluids 18 051705

[11]Ray B,Biswas G,Sharma A 2010 J.Fluid Mech.655 72

[12]Thoraval M J,Li Y,Thoroddsen S T 2016 Phys.Rev.E 93 033128

[13]Pumphrey H C,Elmore P A 1990 J.Fluid Mech.220 539

[14]Pumphrey H C,Crum L A,Bj? rn? L 1989 J.Acoust.Soc.Am.85 1518

[15]Oguz H N,Prosperetti A 1990 J.Fluid Mech.219 143

[16]Zou J,Ji C,Yuan B G,Ren Y L,Ruan X D,Fu X 2012 Phys.Fluids 24 057101

[17]Wang A B,Kuan C C,Tsai P H 2013 Phys.Fluids 25 1518

[18]Deng Q,Anilkumar A V,Wang T G 2007 J.Fluid Mech.578 119

[19]Thoroddsen S T,Thoaval M T,Takehara K,Etoh T G 2012 J.Fluid Mech.708 469

[20]Sigler J,Mesler R 1990 J.Colloid Interface Sci.134 459

[21]Saylor J R,Bounds G D 2012 AIChE J.58 3841

[22]Sussman M,Puckett E G 2000 J.Comput.Phys.162 301

[23]Guo Y L,Wei L,Liang G T,Shen S Q 2014 Int.Commun.Heat Mass.53 26

[24]Wang Z,Li Y,Huang B,Gao D M 2016 J.Mech.Sci.Technol.30 2547

[25]Ohta M,Kikuchi D,Yoshida Y,Sussman M 2011 Int.J.Multiphase Flow 37 1059

[26]Fan W,Qi T,Sun Y W,Zhu P,Chen H 2016 Chem.Eng.Technol.39 1895

[27]Brackbill J U,Kothe D B,Zemach C 1992 J.Comput.Phys.100 335

[28]Cossali G E,Marengo M,Coghe A,Zhdanov S 2004 Exp.Fluids 36 888

[29]Feonychev A I 2007 J.Eng.Phys.Thermophys.80 961