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

1 引 言

隨著工業(yè)化水平的提高,人們對(duì)于機(jī)械產(chǎn)品的安全性、可靠性、節(jié)能性以及服役壽命等提出了更高的要求.據(jù)統(tǒng)計(jì),在工業(yè)生產(chǎn)中約50%的設(shè)備損壞是由磨損引起的.潤(rùn)滑作為一種減小摩擦和降低磨損的有效手段,普遍應(yīng)用于工業(yè)生產(chǎn)中.在對(duì)旋轉(zhuǎn)工作的機(jī)械零部件和設(shè)備進(jìn)行潤(rùn)滑時(shí)將不可避免地涉及油滴撞擊油膜的現(xiàn)象.例如,在對(duì)航空發(fā)動(dòng)機(jī)的軸承進(jìn)行噴油潤(rùn)滑時(shí),為防止?jié)櫥托孤┮鸾够椭鸬痊F(xiàn)象,常采用空氣系統(tǒng)對(duì)軸承室進(jìn)行密封,使軸承室內(nèi)部處于復(fù)雜的油-氣混合狀態(tài),具體表現(xiàn)為油滴和油膜間的相互作用行為[1];在對(duì)齒輪、渦輪等減速器裝置進(jìn)行飛濺潤(rùn)滑時(shí),機(jī)箱中旋轉(zhuǎn)的零部件將潤(rùn)滑油霧化成油滴帶到摩擦副上形成自動(dòng)潤(rùn)滑,部分油滴將會(huì)與機(jī)箱中的油膜發(fā)生碰撞行為;在對(duì)高速電主軸進(jìn)行油-氣潤(rùn)滑時(shí),零部件表面受連續(xù)油滴撞擊并形成潤(rùn)滑油膜層,隨時(shí)間變化,后續(xù)油滴常與沉積的油膜層發(fā)生作用.

目前,在學(xué)術(shù)期刊上較常見(jiàn)的是非油類(lèi)液滴與液膜的碰撞行為研究,其研究?jī)?nèi)容主要分為兩個(gè)方面:一方面是關(guān)注液滴撞擊液膜的流動(dòng)行為;另一方面則集中考察氣泡夾帶現(xiàn)象.

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

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

當(dāng)液滴撞擊液膜形成彈坑狀水花時(shí),容易形成氣泡夾帶現(xiàn)象,這一現(xiàn)象已引起學(xué)者們的廣泛關(guān)注.在海洋和氣候科學(xué)中,氣泡夾帶現(xiàn)象有利于氣液界面氣體進(jìn)行分子交換;而在液體食品和藥劑加工過(guò)程中,氣泡夾帶現(xiàn)象則會(huì)導(dǎo)致泡沫產(chǎn)生,嚴(yán)重影響產(chǎn)品質(zhì)量和生產(chǎn)速率.由此可見(jiàn),確定夾帶發(fā)生的條件以及確定這種現(xiàn)象在數(shù)量上的發(fā)生概率是極為重要的.

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

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

Zou等[16]和Wang等[17]對(duì)大氣泡夾帶現(xiàn)象進(jìn)行了可視化實(shí)驗(yàn)研究,發(fā)現(xiàn)扁長(zhǎng)形液滴撞擊液面后能形成與液滴尺寸幾乎相當(dāng)?shù)臍馀?Deng等[18]研究了液體黏度和表面張力對(duì)氣泡夾帶的影響,發(fā)現(xiàn)夾帶氣泡尺寸隨著毛細(xì)數(shù)(Ca)的增大而減小.

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

通過(guò)上述文獻(xiàn)可以看到,學(xué)者們主要對(duì)液滴撞擊液膜的流動(dòng)行為以及氣泡夾帶現(xiàn)象進(jìn)行了較為系統(tǒng)且深入的研究,并對(duì)各現(xiàn)象的形成機(jī)制和不同物理參數(shù)的影響也進(jìn)行了相應(yīng)探討.

在潤(rùn)滑系統(tǒng)工作過(guò)程中,油滴撞擊油膜時(shí)易發(fā)生氣泡夾帶現(xiàn)象.氣泡加速潤(rùn)滑油氧化變質(zhì),使油膜出現(xiàn)破裂,造成機(jī)械部件的直接接觸并最終影響設(shè)備的服役壽命.此外,油滴往往是連續(xù)撞擊油膜,當(dāng)?shù)谝粋€(gè)油滴撞擊油膜引起氣泡夾帶后,第二個(gè)油滴撞擊的是含有氣泡的油膜.因此,研究油滴與含氣泡油膜的碰撞行為可以為潤(rùn)滑過(guò)程中氣泡影響研究提供參考,同時(shí)對(duì)于深化氣泡動(dòng)力學(xué)規(guī)律的理解也具有重要意義.

當(dāng)前較少見(jiàn)到公開(kāi)發(fā)表的文獻(xiàn)中有關(guān)于液滴與含氣泡液膜碰撞行為的研究,在缺乏相關(guān)理論指導(dǎo)和實(shí)驗(yàn)研究的前提下,數(shù)值模擬研究是揭示其運(yùn)動(dòng)規(guī)律的有效途徑.基于上述背景,本文采用CLSVOF方法模擬了油滴撞擊含氣泡油膜的運(yùn)動(dòng)過(guò)程,重點(diǎn)探討油膜層中氣泡在油滴撞擊作用下的變化規(guī)律.

2 數(shù)值模型計(jì)算方法

2.1 控制方程

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

VOF方法通過(guò)定義流體體積函數(shù)β來(lái)捕捉相界面,當(dāng)β=1時(shí),目標(biāo)區(qū)域僅存在液體;當(dāng)β=0時(shí),目標(biāo)區(qū)域內(nèi)僅存在氣體;當(dāng)0＜β＜1時(shí),目標(biāo)區(qū)域?yàn)閮上嘟缑?流體體積函數(shù)β的控制方程為

式中U為速度矢量,t為實(shí)際時(shí)間.

VOF方法在計(jì)算兩相流動(dòng)時(shí)能保證質(zhì)量守恒,但流體體積函數(shù)β不是連續(xù)函數(shù),不能精確計(jì)算曲率和法向向量等幾何參數(shù)[22].

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

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

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

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

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

式中ε=1.5?r,?r為最小的網(wǎng)格尺寸.

表面張力的計(jì)算通過(guò)連續(xù)表面力(continuum surface force)模型[27]

來(lái)計(jì)算,式中F為表面張力,σ是表面張力系數(shù).

兩相流動(dòng)的質(zhì)量方程和動(dòng)量方程分別表示成

式中g(shù)為重力加速度,P為壓力.

2.2 計(jì)算方法可行性驗(yàn)證

當(dāng)前尚未見(jiàn)到有關(guān)液滴撞擊含氣泡液膜的實(shí)驗(yàn)研究,為驗(yàn)證CLSVOF方法是否能準(zhǔn)確模擬油滴撞擊含氣泡油膜的變化過(guò)程,我們通過(guò)間接驗(yàn)證的方法對(duì)液滴撞擊液膜過(guò)程中的流動(dòng)行為和大氣泡夾帶現(xiàn)象進(jìn)行試算.

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

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

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

圖1 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比 (a)定性對(duì)比;(b)定量對(duì)比Fig.1. Comparison of experiment and simulation:(a)Qualitative comparison;(b)quantitative comparison.

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

綜上所述,CLSVOF方法既能準(zhǔn)確捕捉液滴撞擊液膜過(guò)程中的氣泡形態(tài)變化,也能較精確計(jì)算出水花流動(dòng)形態(tài)的特征參數(shù)變化.因此,本文采用CLSVOF方法模擬油滴撞擊含氣泡油膜這一動(dòng)態(tài)過(guò)程.

2.3 參數(shù)選擇及數(shù)值模型建立

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

表1 計(jì)算時(shí)采用的相關(guān)參數(shù)Table 1.Parameters used in calculation.

圖2所示為初始時(shí)刻油滴撞擊含氣泡油膜的二維軸對(duì)稱(chēng)計(jì)算域示意圖.為使模擬結(jié)果具有三維特征,對(duì)二維模型的對(duì)稱(chēng)軸x添加軸旋轉(zhuǎn)對(duì)稱(chēng)(axisymmetric swirl)條件.定義初始時(shí)刻油滴和油膜剛開(kāi)始接觸,油滴中心和氣泡中心在一條直線上,氣泡位于油膜高度中心.油膜底部為不滑移壁面,環(huán)境壓力保持一個(gè)標(biāo)準(zhǔn)大氣壓.整個(gè)計(jì)算區(qū)域大小為400μm×150μm,并采用四邊形網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行離散.

圖2 計(jì)算域示意圖Fig.2.Schematic of computational domain.

2.4 網(wǎng)格無(wú)關(guān)性分析

為驗(yàn)證網(wǎng)格的無(wú)關(guān)性,圖3給出了4種不同密度的計(jì)算網(wǎng)格下中心油膜高度h0隨實(shí)際時(shí)間t的變化情況,圖中氣泡直徑為10μm.為提高計(jì)算效率同時(shí)保證計(jì)算精度,對(duì)氣泡及周?chē)鷧^(qū)域網(wǎng)格進(jìn)行局部加密,圖3中網(wǎng)格尺寸為加密后的網(wǎng)格.

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

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

3 結(jié)果與分析

3.1 氣泡大小的影響

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

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

圖4 不同大小氣泡的演化過(guò)程 (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)可以看出,當(dāng)氣泡直徑為20μm時(shí),油滴撞擊油膜后氣泡變形過(guò)程與氣泡直徑為10μm時(shí)基本一致,但最終氣泡沒(méi)有發(fā)生破裂而是穩(wěn)定存在于油膜內(nèi),本文將其稱(chēng)為“穩(wěn)定變形”.

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

3.2 氣泡大小對(duì)氣泡變形特征參數(shù)的影響

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

為定量描述油滴撞擊含氣泡油膜壁面的運(yùn)動(dòng)特性,如圖5所示,定義氣泡的橫向長(zhǎng)度為w,縱向長(zhǎng)度為h,使用無(wú)量綱化參數(shù)形狀系數(shù)E來(lái)表征氣泡的形狀變化,并可表示為

形狀系數(shù)E越大,說(shuō)明氣泡總體變形量越大.同樣將運(yùn)動(dòng)變化時(shí)間以無(wú)量綱時(shí)間T表示:

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

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

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

圖6 氣泡大小對(duì)形狀系數(shù)的影響 (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)可以看出,油滴在撞擊含氣泡油膜后,氣泡形狀變化過(guò)程分為兩個(gè)階段:顯著變形階段(T<10)和相對(duì)穩(wěn)定階段(T>10).在顯著變形階段,隨著氣泡直徑增大,同一時(shí)間的氣泡形狀系數(shù)也相應(yīng)增大,且氣泡達(dá)到形狀系數(shù)最大值的時(shí)間也增加,但不明顯.這是因?yàn)闅馀葜睆皆酱?氣泡的表面積越大,越容易受黏性剪切力拖拽而產(chǎn)生變形.在相對(duì)穩(wěn)定階段,各氣泡間的形狀系數(shù)差距較小.此外,氣泡直徑越大,氣泡破裂發(fā)生時(shí)間越晚,且破裂前的形狀系數(shù)呈現(xiàn)小范圍波動(dòng)變化.

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

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

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

3.3 氣泡位置對(duì)氣泡變形特征參數(shù)的影響

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

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

圖7 氣泡位置對(duì)形狀系數(shù)的影響 (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āng)d=10μm時(shí),位于油膜頂端和油膜中心的氣泡形狀系數(shù)隨無(wú)量綱時(shí)間的變化趨勢(shì)保持一致,氣泡均在自由表面發(fā)生破裂且破裂時(shí)間也較為接近.不同點(diǎn)在于T＜6時(shí),位于油膜頂端氣泡的形狀系數(shù)始終大于位于油膜中心氣泡的形狀系數(shù).

由圖7(b)可以看出,當(dāng)d=20μm時(shí),位于油膜頂端和油膜中心的氣泡形狀系數(shù)隨無(wú)量綱時(shí)間的變化趨勢(shì)出現(xiàn)顯著不同,其原因在于位于油膜頂端的氣泡在運(yùn)動(dòng)過(guò)程中出現(xiàn)油膜內(nèi)部破裂現(xiàn)象,而位于油膜中心的氣泡在撞擊過(guò)程中保持較穩(wěn)定的狀態(tài).

從圖7(c)可以看出,當(dāng)d=30μm時(shí),位于油膜頂端和油膜中心的氣泡均在運(yùn)動(dòng)過(guò)程中發(fā)生油膜內(nèi)部破裂現(xiàn)象,但油膜頂端的氣泡最大形狀系數(shù)遠(yuǎn)大于油膜中心的氣泡最大形狀系數(shù).

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

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

3.4 氣泡破裂機(jī)制分析

氣泡破裂現(xiàn)象屬于氣泡動(dòng)力學(xué)范疇.在進(jìn)行氣泡動(dòng)力學(xué)分析時(shí)常引入兩個(gè)無(wú)量綱參數(shù),莫頓數(shù)(M o)和厄缶數(shù)(Eo):

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

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

3.4.1 自由表面破裂

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

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

圖9 氣泡破裂區(qū)域壓力和速度分布圖(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時(shí),氣泡上表面的油膜形成了“頸部”區(qū)域,該區(qū)域內(nèi)部壓力最小值為?13.31 kPa,兩側(cè)壓力為33.6 kPa.可見(jiàn)“頸部”區(qū)域承受較大壓力梯度作用有發(fā)生斷裂的趨勢(shì).由速度分布可知,受毛細(xì)波的影響,氣泡上側(cè)油膜附近速度方向變化較為劇烈.氣-液相界面的曲率發(fā)生變化,造成相界面的表面張力方向發(fā)生改變并促進(jìn)了“頸部”區(qū)域的形成.

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

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

3.4.2 油膜內(nèi)部破裂

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

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

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

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

圖10 氣泡破碎區(qū)域壓力和速度分布圖(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 結(jié) 論

采用CLSVOF方法對(duì)油滴撞擊含氣泡油膜現(xiàn)象進(jìn)行了數(shù)值模擬分析,研究了油膜層內(nèi)部氣泡的運(yùn)動(dòng)演化過(guò)程,考察了氣泡大小及氣泡位置等因素對(duì)受撞擊作用的氣泡變形特征參數(shù)的影響規(guī)律,探索了不同破裂方式下的氣泡破裂機(jī)制,主要得到以下結(jié)論.

1)通過(guò)與實(shí)驗(yàn)結(jié)果比較,驗(yàn)證CLSVOF方法在求解液滴撞擊液膜過(guò)程中的氣泡形態(tài)變化和流動(dòng)特征參數(shù)變化的可行性,并以此為依據(jù),建立了油滴撞擊含氣泡油膜的數(shù)值模型.

2)氣泡大小是決定油滴撞擊油膜后氣泡呈現(xiàn)自由表面破裂、穩(wěn)定變形以及氣泡在油膜內(nèi)部破裂等不同運(yùn)動(dòng)形態(tài)的重要因素.

3)當(dāng)10μm 6 d<20μm時(shí),隨著氣泡直徑增大,氣泡在運(yùn)動(dòng)過(guò)程中的變形量也相應(yīng)增大,并隨時(shí)間變化發(fā)生自由表面破裂,破裂發(fā)生時(shí)間隨氣泡直徑增大而延長(zhǎng);當(dāng)20μm

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

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

[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)[郭加宏,戴世強(qiáng),代欽 2010物理學(xué)報(bào) 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力學(xué)學(xué)報(bào)45 833]

[6]Liang G T,Guo Y L,Shen SQ 2013 Acta Phys.Sin.62 024705(in Chinese)[梁剛濤,郭亞麗,沈勝?gòu)?qiáng)2013物理學(xué)報(bào)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