童寶宏， 陳 武， 張國(guó)濤， 胡曉磊

(1.安徽工業(yè)大學(xué) 機(jī)械工程學(xué)院,馬鞍山 243032；2.特種重載機(jī)器人安徽省重點(diǎn)實(shí)驗(yàn)室,馬鞍山 243032)

1 引 言

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

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

在油膜實(shí)際流動(dòng)過程中，氣泡的形成、破裂及射流總是伴隨著與周圍流體強(qiáng)烈的耦合作用[15-18]。對(duì)氣泡靠近自由液面運(yùn)動(dòng)的研究發(fā)現(xiàn)，初始?xì)馀輭毫臀恢脤?duì)自由液面的形態(tài)及破裂形式具有不同程度的影響[19,20]。氣泡破裂過程會(huì)改變周圍流體的流動(dòng)形式，氣泡周圍流體的壓力分布也更為復(fù)雜。相關(guān)學(xué)者對(duì)液體中氣泡破裂過程產(chǎn)生的壓力變化開展研究，研究結(jié)果表明，氣泡初始位置影響了壁面壓力大小和壓力作用范圍[18]。

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

2 數(shù)值方法及計(jì)算域

2.1 控制方程

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

連續(xù)性方程

(1)

(2)

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

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

(3)

在目標(biāo)域中任意點(diǎn)的流體密度和粘度為兩種流體體積分?jǐn)?shù)的加權(quán)平均值，即

(4,5)

本文同時(shí)考慮到氣相剪切力和表面張力的影響，在數(shù)值計(jì)算求解過程中將氣相剪切力和表面張力合并成同一個(gè)動(dòng)量源相添加到動(dòng)量方程中，即

(6)

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

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

(7)

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

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

圖1 幾何模型及邊界條件

本文選取文獻(xiàn)[22]的部分工況進(jìn)行CFD計(jì)算，表1為根據(jù)文獻(xiàn)[22]的液膜參數(shù)而采用的數(shù)值參數(shù)。

表1 數(shù)值模擬中液膜物理參數(shù)

表1中液膜Re計(jì)算公式為

Re=4Q/(bυ)

(8)

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

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

取數(shù)值模型中部分區(qū)域作為研究區(qū)域，將數(shù)值計(jì)算出的液膜平均厚度δs與文獻(xiàn)[22]實(shí)驗(yàn)測(cè)量的平均膜厚δe對(duì)比。表2為對(duì)比結(jié)果。模擬結(jié)果和試驗(yàn)之間的相對(duì)誤差為

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

(9)

表2 平均膜厚

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

2.3 參數(shù)選擇及模型確定

圖2 流動(dòng)油膜添加氣泡

2.4 網(wǎng)格無關(guān)性驗(yàn)證

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

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

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

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

3 結(jié)果與討論

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

3.1 含氣泡油膜流動(dòng)形態(tài)分析

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

圖4 含氣泡油膜形態(tài)演化

圖5展示了t=10 ms時(shí)有無氣泡夾帶油膜流動(dòng)過程中氣液界面高度對(duì)比。氣泡的運(yùn)動(dòng)和形變誘發(fā)油膜的氣液界面發(fā)生不同幅度的波動(dòng)，在氣泡上游區(qū)域油膜的氣液界面波動(dòng)幅度明顯大于下游，另外，氣泡上游油膜氣液界面有外凸趨勢(shì)，界面高度增加了0.042 mm，出現(xiàn)波峰；氣泡下游油膜氣液界面有內(nèi)凹趨勢(shì)，界面高度減小了0.012 mm，出現(xiàn)波谷。

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

為進(jìn)一步分析氣泡存在對(duì)油膜形態(tài)帶來的影響規(guī)律，本文統(tǒng)計(jì)了氣泡破裂前氣泡上游油膜氣液界面峰值A(chǔ)隨時(shí)間的變化。如圖6所示，縱坐標(biāo)(A-H1)代表不同時(shí)刻含氣泡油膜上游氣液界面峰值與初始油膜氣液界面的差值。隨著氣泡在油膜內(nèi)部發(fā)生形變與遷移，氣泡受到周圍流體的粘性剪切力作用；與此同時(shí)，氣泡對(duì)周圍流體也具有相同的反作用力，使得氣泡上游油膜的氣液界面產(chǎn)生波動(dòng)，油膜氣液界面峰值開始逐漸增加；在t=12.5 ms時(shí)，有無氣泡油膜氣液界面差值達(dá)到最大值0.048 mm；t>12.5 ms后，氣液界面差值沒有發(fā)生顯著變化，而是處于一定幅值范圍內(nèi)的細(xì)微波動(dòng)狀態(tài)，直至氣泡發(fā)生破裂，這主要是因?yàn)槁∑鸬挠湍な鼙砻鎻埩椭亓ψ饔眯枰３衷械男螒B(tài)。

圖6 油膜氣液界面峰值隨時(shí)間變化

3.2 氣泡對(duì)油膜流動(dòng)速度的影響

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

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

圖8展示了油膜內(nèi)某點(diǎn)速度隨時(shí)間變化曲線，可以看出，當(dāng)6.5 ms

圖8 含氣泡油膜內(nèi)某點(diǎn)速度時(shí)變曲線(x =17.5 mm，y =0 mm，z =1.65 mm)

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

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

3.3 氣泡破裂階段對(duì)油膜壓力的影響

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

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

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

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

圖12展示了氣泡破裂過程壁面上一點(diǎn)的壓力隨時(shí)間變化曲線。氣泡發(fā)生破裂后，氣泡空腔快速收縮引起壁面出現(xiàn)第一個(gè)壓力峰值，而后，壁面壓力隨收縮速度放緩而減小；第二個(gè)壓力峰值是由空腔射流沖擊引起的。兩個(gè)壓力峰值作用方向相反，壁面經(jīng)歷正負(fù)壓力交替的變化過程，由此形成對(duì)壁面沖擊載荷。

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

4 結(jié) 論

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

(1) 氣泡夾帶行為引發(fā)氣泡上游油膜的氣液界面產(chǎn)生外凸趨勢(shì)，下游油膜產(chǎn)生內(nèi)凹趨勢(shì)，隨著氣泡運(yùn)動(dòng)演化，氣液界面波動(dòng)幅度(A-Hl)明顯增加。

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

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