氣泡夾帶對壁面油膜流動特性的影響

童寶宏， 陳 武， 張國濤， 胡曉磊

(1.安徽工業大學 機械工程學院,馬鞍山 243032；2.特種重載機器人安徽省重點實驗室,馬鞍山 243032)

1 引 言

含氣泡液體流動普遍存在于自然界及工程領域中，如水輪機工作時造成的空化現象[1]、水面破碎產生的氣泡波浪[2]和進行廢水處理的氣泡[3]等。已有研究表明氣泡的存在具有一定的危害性，尤其在油-氣潤滑系統工作過程中。油-氣潤滑系統工作時，干燥壓縮空氣將微米級的油滴傳輸至潤滑點形成潤滑油膜。相關研究表明[4]，在高速氣流輸送潤滑油過程中，潤滑油與管壁之間容易形成微氣室，氣室內的氣體在壓力作用下破裂為許多氣泡，游離在潤滑油膜內的氣泡會降低油膜承載能力，破壞油膜的完整性，造成潤滑不良并最終損壞摩擦表面[5]。因此，考察含氣泡油膜沿壁面流動的動態演化過程，對于深入研究油膜的潤滑性能和查明含氣泡油膜流動力學機制具有重要的意義。

油膜的形成與流動過程本身是一個復雜的非定常行為[6]。廣大學者對沿傾斜壁面上的液膜流動特性及形態演變過程開展了一系列研究工作，通過實驗觀測和數值模擬考察了壁面傾斜角度、壁面微結構和雷諾數等對液膜波動特性的影響。研究發現，增大壁面傾斜角度和改變壁面形狀會降低液膜的穩定性[7-9]；減少入口流量會導致液膜出現溪流[10]和滴流[11]等現象，破壞液膜的完整性。氣流是油-氣潤滑過程中不可缺少的組成，已有研究采用數值模擬方法深入考察了氣流對液膜厚度及流動速度的影響。結果表明，逆向氣流趨向于降低和平坦化區域中的液體速度，同向氣流減薄液膜厚度[12-14]。

在油膜實際流動過程中，氣泡的形成、破裂及射流總是伴隨著與周圍流體強烈的耦合作用[15-18]。對氣泡靠近自由液面運動的研究發現，初始氣泡壓力和位置對自由液面的形態及破裂形式具有不同程度的影響[19,20]。氣泡破裂過程會改變周圍流體的流動形式，氣泡周圍流體的壓力分布也更為復雜。相關學者對液體中氣泡破裂過程產生的壓力變化開展研究，研究結果表明，氣泡初始位置影響了壁面壓力大小和壓力作用范圍[18]。

國內外學者已對液膜流動行為進行了系統且深入的研究，并對液膜的流動速度、形態、厚度及不同物理量對其影響進行相應考察，取得了豐富的研究成果，但涉及含氣泡油膜流動的動態行為研究鮮有報道。因此，為了深入探究氣泡夾帶影響下潤滑油膜的動態演化過程及內在流動機理，本文采用VOF方法對建立的含氣泡油膜沿傾斜壁面流動的模型進行三維數值仿真模擬，考察氣泡及其破裂行為對油膜的形態和流動特性的影響，進一步查明含氣泡潤滑油膜流動的內在機理。

2 數值方法及計算域

2.1 控制方程

傾斜壁面上含氣泡油膜流動是典型的氣液兩相流現象，自由液面的形狀一般不是規整的幾何面，且含氣泡油膜流動過程中，其相界面的位置隨著空間和時間不斷發生變化。而求解此類問題的關鍵在于如何將相界面進行離散化，以及如何精確追蹤相界面隨時間的變化。目前流體體積法VOF(volume of fluid)廣泛應用到追蹤自由相界面的數值模擬中。關于VOF模型計算的控制方程為

連續性方程

(1)

(2)

VOF方法框架中，在網格上定義了一個額外的場變量，即體積分數函數γ(x,y,z)，為了指示每一個計算單元內一種流體的體積分數,取值范圍在0～1。當函數值為1時，表示計算單元內僅有一種流體填充，γ為0時，則由另一種流體填充，介于0～1時，表示相界面必須穿過該單元。體積分數方程為

?γ/?t+·(υγ)=0

(3)

在目標域中任意點的流體密度和粘度為兩種流體體積分數的加權平均值，即

(4,5)

本文同時考慮到氣相剪切力和表面張力的影響，在數值計算求解過程中將氣相剪切力和表面張力合并成同一個動量源相添加到動量方程中，即

(6)

式中σ為表面張力系數，ki為界面曲率。

固液之間因表面張力而產生的力為壁面粘附力，CSF模型中[21]，對這種力的處理方法不是通過設置壁面處的邊界條件，而是直接將其并入到表面張力動量源相中，此時表面法向量為

(7)

2.2 計算方法可行性驗證

為檢驗VOF方法模擬含氣泡油膜沿壁面流動的可行性，根據文獻[22]的實驗方法和物性參數，建立無氣泡液膜沿壁面流動的三維模型。圖1為物理模型和邊界條件?？紤]到液膜沿傾斜壁面流動呈明顯的對稱性，本文計算時僅建立了1/2模型，對稱軸為y=0的xoz平面，并將其設置為對稱邊界條件，模型與水平面夾角固定在45°，采用無滑移壁面邊界條件，油液從進口流入，出口流出，液相入口設定為速度入口邊界條件。

圖1 幾何模型及邊界條件

本文選取文獻[22]的部分工況進行CFD計算，表1為根據文獻[22]的液膜參數而采用的數值參數。

表1 數值模擬中液膜物理參數

表1中液膜Re計算公式為

Re=4Q/(bυ)

(8)

式中Q為體積流量(m3/s)，b為板寬(m)，υ為運動粘度(m2/s)。

數值計算時，采用有限體積法對控制方程進行離散，壓力速度耦合的求解選用PISO方法，壓力求解采用PRESTO!方法，氣液界面插值選取幾何重構方法，動量求解選擇二階迎風格式。為了獲得更好的收斂性，使用10-5s～10-4s的可變時間步長和恒定全局庫朗數0.25進行瞬態模擬。收斂標準是絕對殘差，監控的變量是出入口液相質量流量和，殘差波動在不同的時間步之間沒有明顯變化且出入口液相質量流量之和小于10-4kg/s時，表示計算達到準穩態，此時計算模擬停止。

取數值模型中部分區域作為研究區域，將數值計算出的液膜平均厚度δs與文獻[22]實驗測量的平均膜厚δe對比。表2為對比結果。模擬結果和試驗之間的相對誤差為

δ=(δe-δs)/δe×100%

(9)

表2 平均膜厚

由表2可知，平均液膜厚度隨著Re數增加而增加，數值模擬結果與實驗結果的吻合性較好，這一規律與文獻[22-23]描述的液膜厚度與Re之間的變化規律相同。 由表2可知，雷諾數在107附近時，數值結果與實驗數據誤差最小，Re為15，24和202時，數值模擬結果與實驗值的相對誤差大于6%，其原因可能是模擬條件和試驗條件之間存在細微差別，此外，文獻[22]也指出，通過測量數據計算出的入口流量Qave與實際入口流量Q之間有著0.02 LPM～0.3 LPM的誤差，進而導致數值模擬時設置的入口流速存在一定的偏差?？傮w來看，VOF方法能較為準確地模擬無氣泡油膜沿壁面流動的動態過程，在此基礎上，本文將采用VOF方法進一步數值模擬含氣泡油膜沿傾斜壁面流動的動態行為。

2.3 參數選擇及模型確定

圖2 流動油膜添加氣泡

2.4 網格無關性驗證

為精確捕捉流動油膜的氣-液界面變化和提高計算效率，對氣泡附近區域網格進行無關性討論。網格劃分過程中，采用六面體結構網格劃分整個計算域，然后對氣泡及其周圍區域進行局部加密，網格最小尺寸分別為0.07 mm,0.06 mm,0.05 mm,0.04 mm和0.03 mm。定義氣泡剛添加到油膜內為初始時刻，即t=0 ms。氣泡貫穿整個油膜流動過程，其直徑變化在模擬過程中能夠直觀體現，此外氣泡直徑也能反映氣泡的氣-液界面變化過程，因此可以用不同網格尺寸下氣泡直徑來表明網格質量。

圖3給出了氣泡破裂前，不同網格尺寸下沿板寬(y)方向xoz平面上氣泡最大直徑d隨時間的變化情況。

圖3 不同網格尺寸下氣泡最大直徑變化對比(y =0 mm)

可以看出，網格尺寸不同，氣泡最大直徑的變化規律基本相似。網格尺寸對氣泡破裂時間有顯著影響，其中網格尺寸為0.07 mm和0.06 mm的最大直徑和破裂時間與網格尺寸0.04 mm和 0.03 mm 相比存在較大差距。而網格尺寸為0.04 mm和0.03 mm氣泡的最大直徑隨時間變化基本相同，且破裂時間也非常接近?？傮w來看，為保證數值計算結果具有參考意義的同時考慮計算效率，本文選擇最小網格尺寸為0.04 mm的網格對氣泡流動區域的油膜進行劃分。

3 結果與討論

數值研究過程中，將氣泡中心設置在沿油膜流動方向x=8 mm、板寬方向y=0 mm及垂直于壁面方向z=0.9 mm，此時氣泡中心與液膜表面距離為0.93 mm。

3.1 含氣泡油膜流動形態分析

圖4展示了含氣泡油膜流動過程中油膜形態的演化過程。初始時刻，氣泡與流體之間的相對速度為零。氣泡受粘性剪切力及自身浮力的作用，沿著壁面向下游運動，同時產生遠離壁面的上浮運動。油膜受氣泡上浮作用影響，在氣泡上游的油膜表面出現液面隆起現象。隨后氣泡與油膜表面接觸發生破裂，氣泡破裂后，在油膜中形成空腔，空腔在表面張力和浮力作用下逐漸減小，直至形成射流。這與文獻[16,25]描述的現象一致，也進一步證明了本文數值模擬的準確性。

圖4 含氣泡油膜形態演化

圖5展示了t=10 ms時有無氣泡夾帶油膜流動過程中氣液界面高度對比。氣泡的運動和形變誘發油膜的氣液界面發生不同幅度的波動，在氣泡上游區域油膜的氣液界面波動幅度明顯大于下游，另外，氣泡上游油膜氣液界面有外凸趨勢，界面高度增加了0.042 mm，出現波峰；氣泡下游油膜氣液界面有內凹趨勢，界面高度減小了0.012 mm，出現波谷。

圖5 油膜流動過程的氣液界面(t =10 ms)

為進一步分析氣泡存在對油膜形態帶來的影響規律，本文統計了氣泡破裂前氣泡上游油膜氣液界面峰值A隨時間的變化。如圖6所示，縱坐標(A-H1)代表不同時刻含氣泡油膜上游氣液界面峰值與初始油膜氣液界面的差值。隨著氣泡在油膜內部發生形變與遷移，氣泡受到周圍流體的粘性剪切力作用；與此同時，氣泡對周圍流體也具有相同的反作用力，使得氣泡上游油膜的氣液界面產生波動，油膜氣液界面峰值開始逐漸增加；在t=12.5 ms時，有無氣泡油膜氣液界面差值達到最大值0.048 mm；t>12.5 ms后，氣液界面差值沒有發生顯著變化，而是處于一定幅值范圍內的細微波動狀態，直至氣泡發生破裂，這主要是因為隆起的油膜受表面張力和重力作用需要保持原有的形態。

圖6 油膜氣液界面峰值隨時間變化

3.2 氣泡對油膜流動速度的影響

油膜流動表速度和油膜內速度是反映油膜形成速率和傳熱特性的兩個重要參數。圖7展示了有無氣泡夾帶時油膜表速度分布，其中區域I-III分別表示氣泡的上游、頂部和下游的油膜?？梢钥闯?，氣泡鄰域內的油膜流動表速度受氣泡影響產生擾動，在不同區域內存在速度梯度分布，這主要是因為隨著氣泡向下游運動，油膜與氣泡之間存在速度差，致使流動油膜受到氣泡表面粘性阻力作用，阻礙了氣泡周圍流體的運動。在I-III區域中，氣泡對區域III內的油膜流動速度阻礙程度最大，速度峰值減小了0.021 m/s；區域I內的油膜阻礙程度次之,速度峰值減小了0.018 m/s；區域Ⅱ內的油膜流動速度阻礙作用最小，尤其靠近氣泡中心位置，其影響程度最小。以上分析說明氣泡的存在能夠延遲油膜表面流動時間，不利于摩擦表面潤滑油膜快速形成。

圖7 含氣泡油膜表速度分布(y =10 mm，t =18.7 ms)

圖8展示了油膜內某點速度隨時間變化曲線，可以看出，當6.5 ms

圖8 含氣泡油膜內某點速度時變曲線(x =17.5 mm，y =0 mm，z =1.65 mm)

圖9展示了t=18.7 ms時，沿流動方向上、區域I-III有無氣泡夾帶時油膜流動各方向的速度分布情況?？梢钥闯?，區域I-III油膜減小的表速度主要體現在方向x和方向z的速度Vx和Vz，根據3.1節所述，油膜中夾帶的氣泡因其上浮運動誘發油膜表面出現隆起現象，受粘性剪切力作用致使油膜流動受阻，原本沿方向x和方向z的速度減小。氣泡阻礙的油膜因氣泡上浮運動而產生方向z的運動趨勢，導致氣泡周圍的油膜表面高于遠處油膜，流經氣泡附近的流體發生分流現象，引起沿板寬(y)方向的分速度，沿板寬(y)方向的速度Vy略有增加，這也進一步解釋了3.1節所提到的，氣泡上游油膜的氣液界面峰值在增加到最大值后沒有繼續增加的原因。

圖9 含氣泡油膜表面不同方向速度(y =10 mm，t =18.7 ms)

3.3 氣泡破裂階段對油膜壓力的影響

氣泡在壁面附近破裂時，會引起周圍流體產生復雜的物理現象，這將影響著油膜的潤滑性能。圖10 展示了氣泡破裂過程油膜壓力分布云圖，其中黑色曲線為氣液相界面曲線。氣泡發生破裂后，油膜中的空腔在表面張力和浮力作用下逐漸減小，直至形成射流。從圖10(a,b)可以看出，在氣泡破裂前期，氣泡內部壓力由原來的69.8 Pa減少到4.3 Pa，氣體迅速釋放，氣泡空腔體積快速減小，空腔外圍流體隨著空腔壁向中心運動，但外圍流體由于慣性不易迅速跟上，致使在其鄰域形成負壓區域[26]，負壓最小值為-48.5 Pa。隨后空腔受慣性效應減弱，空腔體積收縮速度放緩，如圖10(c)所示，空腔下方的負壓區域也進一步縮小，而負壓區域主要集中在空腔下游，上下游油膜壓力差值為38.5 Pa，這主要是由于空腔上下游非對稱性流動引起的，空腔上游流體運動方向與空腔收縮速度方向相同，空腔下游流體則相反。在氣泡破裂后期，如圖10(d)所示，t=32.3 ms時，空腔底部寬度基本消失，空腔類似倒置的圓錐狀，此時空腔上下游流體向空腔底部聚集，空腔底部流體形成局部高壓區域，該高壓區域也是空腔形成射流的主要原因。

圖10 氣泡破裂過程油膜壓力分布

為定量描述氣泡破裂階段對油膜壓力的影響情況，選取氣泡下方壁面壓力來間接分析油膜壓力變化規律。圖11展示了含氣泡油膜破裂過程中氣泡下方壁面壓力分布，其中壓力正負號代表壁面所受壓力方向(壓力值為負，代表壁面受z軸正方向的力)。t=30.4 ms時，氣泡為未發生破裂壁面壓力分布；t=30.7 ms時，空腔底部壁面形成一個負壓區；t=30.9 ms時，區域負壓值繼續增大，出現第一個壓力峰值，數值為-30.8 Pa；t=31.8 ms時，壁面壓力減??；t=31.9 ms時，壁面呈現正負壓同時存在的現象；t=32 ms時，與圖11(e)中壁面壓力相比，壁面處的壓力有所增加，壁面所受負壓力全部消失；t=32.3 ms時，壁面壓力出現第二個壓力峰值，數值為 60.3 Pa，這表明氣泡破裂過程可能對摩擦表面造成損害；t=32.7 ms時，隨著空腔射流演化，壁面壓力值逐漸降低。

圖11 氣泡破裂過程壁面壓力分布

圖12展示了氣泡破裂過程壁面上一點的壓力隨時間變化曲線。氣泡發生破裂后，氣泡空腔快速收縮引起壁面出現第一個壓力峰值，而后，壁面壓力隨收縮速度放緩而減??；第二個壓力峰值是由空腔射流沖擊引起的。兩個壓力峰值作用方向相反，壁面經歷正負壓力交替的變化過程，由此形成對壁面沖擊載荷。

圖12 氣泡破裂階段壁面壓力時變曲線(x =22.7 mm，y =0 mm，z =0 mm)

4 結 論

本文基于VOF數值方法，建立了含氣泡油膜沿傾斜壁面流動的三維兩相計算流體力學模型，考察了含氣泡油膜沿傾斜壁面的流動特性，并對比分析了油膜中夾帶的氣泡對其影響情況，得到如下結論。

(1) 氣泡夾帶行為引發氣泡上游油膜的氣液界面產生外凸趨勢，下游油膜產生內凹趨勢，隨著氣泡運動演化，氣液界面波動幅度(A-Hl)明顯增加。

(2) 油膜夾帶氣泡的形變和遷移誘發氣泡周圍微流場的速度擾動現象，導致氣液界面處產生非均勻速度梯度分布。氣泡抑制了鄰域內油膜沿流動方向和垂直壁面方向的流動速度，引發了沿板寬方向上的表速度及波動。

(3) 氣泡破裂時，油膜空腔鄰域內流體產生正負壓力交替的波動過程。在破裂初期，受慣性效應影響，空腔周圍的油膜形成負壓區域，隨空腔收縮速度放緩和射流形成，油膜逐漸由負壓轉向正壓?？涨簧淞鞯男纬蓪a生比氣泡破裂初期更強的壁面壓力沖擊，壁面承受一定的交變載荷作用。