考慮油液流變時變特性和深度特征尺寸的靜壓轉臺流動分析

金玉峰, 段鸞芳, 吳偉偉

(1.江蘇省常州技師學院 智能制造學院, 江蘇 常州 213032; 2. 南京工業職業技術大學 電氣工程學院, 江蘇 南京 210023; 3.揚州大學 機械工程學院, 江蘇 揚州 225127)

引言

液壓油是數控機床加工過程中不可或缺的材料,主要起到潤滑與冷卻的作用[1-2]。當液壓油用于導軌、絲杠及軸承(如靜壓導軌、靜壓絲杠及靜壓軸承)中時,還需要具有一定的承載能力[3-7]。靜壓轉臺是數控機床中的常見部件,其以靜壓軸承為支承,液壓油作為工作介質,在壓力作用下形成壓力油膜,具有較好的承載剛性[8-9]。

靜壓轉臺的結構是影響剛度和承載能力等的主要因素。張瑞濤等[10]利用Fluent軟件對節流器孔徑及軸承間隙對剛度的影響進行了分析,基于分析結果設計了具有較高剛度的靜壓轉臺。仝志偉等[11]分析了支承布局對平面度誤差的影響,并對支承布局進行重構,可減小平面度誤差,提高轉臺的承載能力。王京等[12]發現雙環形油腔較傳統油腔具有更好的承載能力。王永柱等[13]對扇形轉臺中油腔形狀對承載能力的影響進行了分析對比,結果表明口字型油腔具有更好的承載能力。張慶鋒[14]分析研究了位移率和間隙對油膜承載能力的影響,根據結果優化了轉臺結構。申峰等[15]具體分析了封油邊形狀特征參數對轉臺承載能力的影響,利用分析結果優化封油邊參數可顯著改善承載能力。MICHALEC M等[16]利用CFD模型和MATLAB提出一種雙參數優化方法,對靜壓軸承的油腔形狀進行了優化,可有效降低能量損失。ZHANG Yanqin等[17]針對雙矩形腔的靜壓轉臺中載荷對油膜性能的影響展開了研究,研究結果表明,載荷大小與油膜的平均溫升之間存在非線性關系。YADAV S K等[18]利用有限元分析對比分析了圓形、橢圓形、方形及矩形油腔對油膜剛度系數的影響,結果表明油腔結構對油膜承載能力具有重要影響。

在靜壓轉臺的實際運行工作過程中,通過實驗發現液壓油的油液品質(油液的黏度)對機床零部件運作時的性能也具有一定影響[19-21]。液壓油在使用過程中,隨著時間的推移,油液的黏度會發生變化,需要了解該變化對承載能力的影響。本研究將以承載能力較好的雙環形油腔結構為研究對象,利用格子玻爾茲曼方法(LBM)分析討論封油邊厚度對油液黏度的時變特性的影響,為后續選擇合適的油液黏度和封油邊厚度提供理論依據。

1 流變測試

表1 兩種油品流變數據擬合結果

圖1 兩種油品實驗數據與擬合結果對比

根據表1結果,將其與流變實驗的數據相比較,結果顯示出較好的一致性,且表1中的相關系數無限接近于或等于1,因此擬合獲得的流變模型是可信的。由表1結果可知,1#和2#油品均呈現出非牛頓特性,具體為帶屈服特征的Herschel-Bulkley型流體,僅當剪切應力大于初始屈服應力時,流體才會發生剪切運動。

2 適應Herschel-Bulkley流體的LBM方法

2.1 虛外力項LBM方法

LBM的演化方程為:

(1)

為了解決由非牛頓特性引起的仿真離散及精度差的問題,引入“虛外力項”用以描述Herschel-Bulkley特征,對于一般含外力項的LBM演化方程為[22]:

feq(r,t)]+δtF′

(2)

式(2)中包含了LBM的兩個主要過程——粒子遷移和粒子碰撞,其中粒子遷移具體可描述為:

f(r+eiδt,t+δt)=f+(r,t)

(3)

粒子碰撞具體可描述為:

(4)

式中,r—— 位置矢量

t—— 某時刻

f—— 粒子分布函數

feq—— 平衡態分布函數

其具體表達式為:

(5)

式中,ρ—— 格子密度,通常取1

u—— 速度矢量

ωi—— 權重系數

對于LBM的D2Q9模型,其具體取值如下:

(6)

其中,ei為D2Q9模型的速度配置,具體取值為:

(7)

(8)

黏度一般可根據流變模型計算獲得,對于牛頓流體,黏度不隨剪切率的變化而變化,其為恒定常數,則松弛時間也為恒定值;而對于非牛頓流體,黏度隨剪切率的變化而變化,因此松弛時間為一變化值。對于液壓油所屬的Herschel-Bulkley流體,其流變模型為:

(9)

式中,τh—— Herschel-Bulkley流體的剪切應力

τh0—— 初始屈服應力

μh0—— 黏度系數

n—— 冪律指數

為了方便計算,前人針對Herschel-Bulkley流體提出了修正黏度-剪切率方程為:

(10)

(11)

其中,DII為應變率張量第二不變量,其可通過下式獲得:

(12)

式中,l—— 計算時的維數,此處l=2

Sαβ—— 應變率張量

其具體計算公式為:

(13)

對于式(4)中的虛外力項,具體計算方式如下:

(14)

(15)

結合式(14)和式(15),可知離散虛外力項具體可以描述為:

(16)

2.2 理論驗證

為了驗證上述提出的虛外力項方法的可靠性,利用Herschel-Bulkley流體在泊肅葉流中的理論解進行驗證,其理論解為[23]:

(17)

式中,H為泊肅葉流中兩平行平板之間的距離,yτ為一臨界點,表征了剪切應力是否達到屈服應力,具體值為yτ=τh0/(?p/?x),?p/?x為在x方向上的壓力梯度。具體參數設置如下:H=1,平板長度L=1,壓力梯度?p/?x=-4×10-3,黏度系數μh0=0.1,冪律指數n=0.9,初始屈服應力τh0=6×10-4。在數值模擬計算時,格子數取為200×200,此處僅為了驗證模型的可靠性,因此直接使用了無量綱數進行模擬與對比。經過仿真,對比結果如圖2所示,對比結果表明提出的虛外力項LBM具有較好的可靠性,可用于Herschel-Bulkley流體的仿真中。

圖2 泊肅葉流LBM數值解與理論解的對比

2.3 靜壓轉臺仿真步驟

將上述提出的虛外力項LBM用于液壓油在靜壓轉臺中的流動分析,其具體過程如下。

(1) 初始化參數。設定初始參數,主要包括格子數(由轉臺的尺寸決定)、入口壓力或入口速度、格子密度、初始粒子分布函數、初始黏度、初始松弛時間等。在此步驟中,需進行單位轉換,LBM是一種無量綱方法,此處以雷諾數為準則數,由實際尺寸和格子大小可確定長度比例,自定義一較小格子速度,由實際速度值和格子速度可確定速度比例,基于式(18),保證實際雷諾數與格子雷諾數相等,即可獲得運動黏度比例,以此類推可獲得其他比例,進而計算出其他格子物理量,完成實際物理量與仿真格子參數之間的轉換[24]:

(18)

其中,v為特征速度,d為特征長度,υ為運動黏度。

(2) 初始化平衡態分布函數。在初始階段可以自定義一個函數作為初始平衡態分布函數,后續根據式(5)實際計算迭代獲得新的平衡態分布函數。

(3) 執行碰撞步式(4)和遷移步式(5)。在執行碰撞步時,需根據式(13)～式(16)計算虛外力項,用以表征Herschel-Bulkley流體的非牛頓特性。

(4) 計算松弛時間。由式(13)計算應變率張量,結合式(11)、式(12)計算剪切率,由式(8)和式(10)計算下一次迭代所需的松弛時間。

(5) 邊界處理。不同于傳統有限元方法的邊界處理,LBM的邊界處理主要是對邊界粒子的狀態進行定義,目前比較成熟的方法有反彈格式、非平衡外推等,此處選用非平衡外推進行邊界處理。

(6) 判斷是否計算結束。設定一較小值為對照標準,以某一物理量為基準,計算上一循環與本次循環的結果差值的絕對值,當此值小于上述設定的對照標準時,即認為計算已滿足條件,跳出迭代計算,否則進行下一輪迭代循環計算。

(7) 計算宏觀物理量。參考式(19)進行物理量計算,并再次進行量綱轉換獲得宏觀物理量。

(19)

式中,ρ′ —— 仿真中油液密度

u—— 仿真中獲得的速度

3 不同靜壓轉臺結構中的流動分析

3.1 靜壓轉臺結構

基于雙環形靜壓轉臺較傳統靜壓轉臺具有更好的承載能力,采用如圖3所示的雙環形油腔結構的靜壓轉臺為研究對象。

圖3 雙環形靜壓轉臺結構

油液入口1的半徑r1設置為4 mm,油液入口2和3的半徑r2均設置為2 mm,內圓環的內圓半徑R2設置為20 mm,外圓半徑R3設置為24 mm,外圓環的內圓半徑R4設置為46 mm,外圓半徑R5設置為50 mm,油腔的總深度為H,封油邊深度為h,為了研究油液黏度時變特性下深度對油液在腔體內部流動的影響,其具體設置如表2所示。由于雙環形油腔是一典型的對稱結構,為提高運算效率,在以下的仿真中均取一半結構進行分析,仿真時徑向的格子數設置為1000,軸向的格子數由具體尺寸等比例設定。

表2 油腔深度的設置

3.2 1#油品流動分析

根據2.1節及2.3節中的描述,結合MATLAB 2017b軟件將1#油品的流變方程代入LBM方法中進行仿真分析,設定油液在3個入口的速度均為0.1 m/s,可以得到如圖4～圖7的結果,其中圖4為速度云圖,圖5為頂層油膜速度大小的分布情形,圖6為頂層油膜速度矢量的角度情形,圖7為頂層油膜的壓力分布結果。通過對比圖4a～圖4c可知,在保持油腔深度h2(h2=H-h)不變的情形下,不斷增加h,速度云圖的變化不明顯。根據圖4d～圖4f可以發現,在保持h不變的情形下,不斷增加h2,速度云圖差異較為明顯,在入油口1處,由1個主渦逐漸發展成2個主渦相連,在入油口2處,主渦強度逐漸增強。在環形臺階上方,速度增強的現象逐漸減弱。

圖4 1#油品在6種結構參數下的速度云圖

圖5 1#油品頂層油膜的速度大小分布

圖6 1#油品頂層油膜的速度矢量角度

圖7 1#油品頂層油膜的壓力分布

圖5a和圖5b分別描述了Case 1～Case 3和Case 4～Case 6情形下頂層油膜的速度大小分布情況,兩圖具有相似的變化趨勢,在兩個環形臺階上方的頂層油膜具有較大的速度,且入油口2輸入的油液增強了其左右兩側環形臺階上方頂層油膜的速度,造成同一臺階上方的速度大小有較明顯的不同。圖5b表明,相同的h下,較小的h2會導致環形臺階上方出現較大的速度,與圖4d～圖4f展現出的速度云圖相吻合。圖6中速度矢量與水平線的夾角展示了近乎一致的變化,該角度的變化與形成的渦流在頂層處的形狀相關,由圖4可看出在頂層油膜渦流具有相似形狀,與圖6量化的角度大小基本吻合。

圖7的壓力分布結果間接展示了轉臺的承載能力,6種情形的結果均表明,中心凹槽的壓力遠大于環形凹槽的壓力,間接說明了入油口2處油液的輸入形成了外圍渦流大大增加了中心凹槽的壓力,起到了增壓的作用,提升了轉臺的承載能力。此外,通過對比發現H值越小,頂層油膜所受的壓力越大,且呈現指數級增長的趨勢,中心凹槽的壓力與環形凹槽的壓力差值也越大。

3.3 2#油品流動分析

圖8～圖11為2#油品仿真分析的結果,圖8為速度云圖,其中圖8a～圖8c中,保證h2不變,逐漸增大封油邊h,速度云圖未有明顯差異。圖8d～圖8f中,保證h不變,h2較小時,中心凹槽內的渦流較明顯,隨著h2逐漸增大,中心凹槽內的渦流逐漸減弱,但均比圖8a～圖8c中的渦流明顯。

圖8 2#油品在6種結構參數下的速度云圖

圖9 2#油品頂層油膜的速度大小分布

圖10 2#油品頂層油膜的速度矢量角度

圖11 2#油品頂層油膜的壓力分布

圖9中展示了頂層油膜速度大小的分布情況,圖9a表明隨著h的逐漸增大,環形臺階上方的速度逐漸減小,圖9b表明隨著h2的逐漸增大,頂層油膜的速度在徑向分布上均有一定程度的減小,在環形臺階上方的速度減幅較大。圖10中所展示的速度矢量角度表明在6種情形下頂層油膜速度與水平線夾角的變化趨勢幾乎一致,與圖8所示的速度云圖基本吻合。

圖10中展示了頂層油膜受壓情況,由圖可知,入油口2中油液的輸入同樣增加了中心凹槽內的壓力,起到了增壓作用,且Case 4～Case 6的結果表明,H較小時,壓力的增幅較明顯。

通過綜合對比1#油品與2#油品的仿真結果可以發現,隨著時間的推移,油液中混入的雜質等改變了油液的流變特性,增加了油液的黏度,對比圖3與圖7、圖4與圖8可知,油液黏度的增加使油液的流速明顯減小,形成的主渦強度較弱。對比圖6與圖10的壓力分布可知,油膜的承壓能力也大大減弱。

4 結論

經過實際加工過程觀察及流變實驗分析發現,液壓油在經過一段時間的使用后,油液的流變特性發生了明顯變化,且剛使用與使用一段時間后的油液均呈現出典型的非牛頓特性,針對非牛頓油液流變隨時間變化的特點,提出了一種虛外力項LBM方法,對油液在雙環形轉臺的內部流動進行分析,通過仿真結果對比后,可以得出以下結論:

(1) 雙環形轉臺的封油邊深度h變化時,速度云圖及頂部油膜速度變化主要集中在環形臺階上方,但壓力分布差異較明顯,隨著封油邊深度h的增加,入油口2油液的輸入對中心凹槽的增壓作用減弱;

(2) 雙環形轉臺的油腔深度h2不斷增大時,兩種油品的頂層油膜速度沿徑向分布減小,且在環形臺階上方的速度減幅較大。兩種油品中心凹槽內的速度云圖有明顯不同,承載能力呈現“斷崖式”減小;

(3) 液壓油的流變特性隨時間的推移發生變化,黏度不斷增大,減弱了油液的流動,導致2#油品情形下,頂層油膜的速度減小,承載能力也大幅減弱。對于剛使用的1#油品,隨著油腔深度增加,中心凹槽的主渦區域不斷擴大,而2#油品的中心凹槽主渦區域不斷減小,此現象也是由油液黏度的變化導致。為保證轉臺的正常工作,應定期更換轉臺用液壓油。