高壓供油泵凸輪軸潤滑特性的熱流固耦合分析*

許 磊 張 翼 徐春龍 李斌茂 馬運昌 唐詩澤 張 宇

(1.中北大學(xué)能源與動力工程學(xué)院 山西太原 030051；2.中國北方發(fā)動機研究所 天津 300400)

高壓供油泵作為內(nèi)燃機燃油供給系統(tǒng)的主要部件，在供油過程中，主軸、座圈等部件承受著大幅值的交變載荷，工作環(huán)境相對惡劣[1-2]。其中，柴油潤滑用于減少曲軸軸頸摩擦副之間的磨損。柴油在零件表面附著一層微米級的油膜，既能避免零件間的相互接觸，又能迅速帶走摩擦產(chǎn)生的熱量，其上油膜特性是影響零件壽命的一個重要方面[3-4]。宋蘭蘭[5]基于Reynolds提出的動壓潤滑理論，通過聯(lián)立 Reynolds方程、黏溫方程和能量方程分析了黏溫效應(yīng)影響下的滑動軸承動力特性及系統(tǒng)穩(wěn)定性。由于考慮軸承表面粗糙度等影響因素的油膜能量方程過于復(fù)雜，文獻[5]對能量方程進行了較大的簡化，未能很好地反映表面油膜溫度變化。PANG等[6]基于油膜厚度比理論，推導(dǎo)了潤滑狀態(tài)與扭矩載荷值的對應(yīng)關(guān)系。主軸發(fā)生形變后對真實油液分布也有較大影響，楊東鵬[7]運用Workbench建立了液體動靜壓軸承油膜和軸瓦的流固耦合分析模型，采用單向流固耦合分析法研究了油膜的相關(guān)特性對軸瓦的應(yīng)力變形情況所造成的影響。

在主軸工作時，主軸軸頸的旋轉(zhuǎn)作用使得油膜對軸頸產(chǎn)生摩擦熱，主軸內(nèi)部油膜的平均溫度會因此而升高[8]。當(dāng)軸徑內(nèi)部溫度升高后，黏溫效應(yīng)使得潤滑油黏度降低，導(dǎo)致油膜產(chǎn)生的熱量減少，同時由于潤滑油黏度的降低使得潤滑油流量增加，其端泄量增大，使得軸徑得到更好的冷卻，進而導(dǎo)致軸徑的溫度下降。故主軸內(nèi)部油膜的熱量及溫度受到黏性耗散、潤滑油黏溫特性、潤滑油通流特性的復(fù)雜耦合作用。柴油黏度較之于機油更小，受黏溫效應(yīng)影響，柴油黏度會更小，更易使摩擦副相對運動阻力增大，處于混合潤滑和干摩擦概率增大，從而導(dǎo)致主軸損傷。故針對主軸采用柴油潤滑時考慮黏溫效應(yīng)顯得極為重要。

綜上所述，學(xué)者們在結(jié)合能量方程求解油膜溫度時，由于受多因素的影響能量方程過于復(fù)雜，均會對能量方程進行簡化，而簡化后的能量方程未必能反映真實的油膜溫度?；诖耍疚淖髡咄ㄟ^數(shù)學(xué)推導(dǎo)，建立一種可代入溫度實驗數(shù)據(jù)或?qū)胲浖嬎銣囟葦?shù)據(jù)的黏溫效應(yīng)一維計算仿真模型，以更好地反映真實的油膜溫度；考慮到主軸內(nèi)部油膜的熱量及溫度受多因素復(fù)雜耦合作用，主軸真實運轉(zhuǎn)時油膜特性受主軸形變影響，建立Workbench平臺下雙向熱流固耦合模型，研究黏溫效應(yīng)對主軸摩擦副油膜特性的影響；對黏溫效應(yīng)影響下柴油潤滑主軸的油膜特性進行分析，為主軸設(shè)計提供依據(jù)。

1 MATLAB數(shù)值模擬

1.1 方程的量綱一化

計算流體動力學(xué)(CFD) 以計算機為載體，通過離散求解流動方程[9-10]。下文將根據(jù)離散求解流動方程的方式建立基于MATLAB的一維仿真模型。

常見的穩(wěn)定運轉(zhuǎn)的徑向滑動軸承，其一般使用的雷諾方程二維形式如下式:

(1)

式中：p為油膜壓力；h為油膜厚度；η為潤滑油的動力黏度；U為軸頸旋轉(zhuǎn)的線速度。

對式(1)進行量綱一化。采用x的量綱一化形式為：x=rφ→φ=x/r，0≤φ≤2π，其中φ為由軸與軸承中心連線沿軸旋轉(zhuǎn)方向角度；r為軸承半徑。采用z的量綱一化形式為：z=Lλ→λ=z/L，0≤λ≤1，其中L為軸承寬度；λ為量綱一軸承寬度。油膜厚度h的量綱一化形式為：H=h/c，其中H=1+εcosφ，ε=e/c，e為偏心量；c為軸承初始半徑間隙；ε為偏心率；H為量綱一油膜厚度。對于壓力p的量綱一化形式，先以一未定的p0值作為p的量綱一化數(shù)，則有：P=p/p0。

將以上量綱一化公式代入式(1)，并化簡得：

(2)

(3)

式中：φ為由軸與軸承中心連線沿軸旋轉(zhuǎn)方向角度；L為軸承寬度；n為轉(zhuǎn)速；η為潤滑油黏度。

U=2πnr，由此可得參考壓力p0隨轉(zhuǎn)速升高而增加。

1.2 黏度變化的量綱一化方程

當(dāng)考慮黏溫特性影響時，軸承間隙流場內(nèi)的黏度并非恒定，其黏度隨溫度變化曲線如圖1所示?？梢姡ざ入S流場溫度變化而發(fā)生變化。對黏度求偏導(dǎo)時，其計算值不為0，原量綱一化方程改進為以下形式：

圖1 柴油黏溫特性曲線

(4)

其中計算黏度變化時采用基于柴油黏溫特性實驗數(shù)據(jù)[11]經(jīng)由MATLAB擬合工具箱進行函數(shù)擬合后的函數(shù)關(guān)系式(5)，其中t為溫度。

η=-5.97×10-9×(t-273)3+1.58×10-6×

(t-273)2-0.000 155 2×(t-273)+0.006 557

(5)

1.3 考慮黏度變化的Reynolds方程的差分形式

采用有限差分法求解油膜的壓力分布，將軸承內(nèi)部的油膜沿周向展開，設(shè)定周向和軸向的步長并為其劃分網(wǎng)格，將 Reynolds方程轉(zhuǎn)化為一組差分方程以此求得各節(jié)點的壓力值。i方向的步長為Δφ，j方向的步長為Δλ。按差分原理，將油膜厚度H沿周向i的一階差商等效為油膜厚度H在周向i的導(dǎo)數(shù)，如式(6)所示；將油膜壓力P和黏度的差分形式等效為2個方向上的偏導(dǎo)數(shù)，如式(7)、(8)、(9)和(10)所示。

(6)

(7)

(8)

(9)

(10)

從而可以求得油膜壓力P沿周向和軸向的二階差商，并將其近似代替Reynolds方程中的二階偏導(dǎo)數(shù)可得式(11)和式(12)。

(11)

(12)

將式(11)、(12)代入量綱一化處理后的Reynolds方程中，便可得到Reynolds方程的差分形式，經(jīng)整理寫成如下形式：

(13)

其中令：

E=A+B+C+D

F=6ηi,jΔφ(Hi+1/2,j-Hi-1/2,j)

推得計算壓力分布的迭代公式如下：

(14)

超松弛迭代公式為

Pi,j=

(15)

其中取α=1.3。

基于以上差分方程可編寫考慮黏溫效應(yīng)影響的油膜壓力程序，建立基于MATLAB方法的仿真模型。該模型基于Reynolds方程的變形進行數(shù)學(xué)推導(dǎo)，實現(xiàn)可代入溫度實驗數(shù)據(jù)或?qū)胲浖嬎愕臏囟葦?shù)據(jù)的黏溫效應(yīng)一維計算仿真模型。與傳統(tǒng)方法簡化的能量方程相比，該方法得到的能量方程包含表面粗糙度等復(fù)雜因素對溫度的影響，能更好地反映實現(xiàn)油膜實際溫度，計算結(jié)果與三維仿真結(jié)果吻合更好。文中將Fluent軟件計算的溫度數(shù)據(jù)導(dǎo)入一維計算仿真模型，完成后續(xù)一維仿真計算。

2 雙向熱流固耦合仿真方法

2.1 模型與網(wǎng)格

轉(zhuǎn)子式高壓供油泵凸輪軸的結(jié)構(gòu)如圖2所示，由座圈、凸輪軸和油膜3部分組成。挺柱產(chǎn)生的交變載荷作用于主軸套上的平面區(qū)域。位于主軸和座圈間的油膜厚度為0.05 mm。為提高仿真準(zhǔn)確性，在0.05 mm厚度的油膜上劃分5層網(wǎng)格。如圖3所示，利用ICEM繪制流體域網(wǎng)格。油膜整體采用六面體網(wǎng)格，周向節(jié)點數(shù)496個，軸向節(jié)點數(shù)120個，網(wǎng)格總體質(zhì)量達到0.75以上。經(jīng)計算檢驗，當(dāng)計算域網(wǎng)格總數(shù)為40萬左右時，油膜最大壓力、有效載荷計算值精度較好。

圖2 凸輪軸仿真模型示意

圖3 油膜網(wǎng)格示意

2.2 邊界條件的設(shè)定

流體部分的邊界條件根據(jù)實際情況在表1中給出。固體部分根據(jù)凸輪升程公式[12]算出彈簧隨轉(zhuǎn)角的壓縮量，并結(jié)合座圈三面所承受柱塞壓力，進一步計算得到交變載荷隨時間變化的曲線如圖4所示。

表1 油膜邊界條件

圖4 座圈三面壓力隨時間的變化曲線

主軸材料為20CrMnTi[13]，泊松比為0.25，材料密度為7 800 kg/m3，比熱容為462 J/(kg·K)，彈性模量為207 GPa，熱膨脹系數(shù)為1.27×10-5℃-1，導(dǎo)熱系數(shù)為33.27 W/(m·K)。假設(shè)主軸與周圍的空氣傳熱為自然對流，復(fù)合傳熱系數(shù)為10 W/(m2·K)。環(huán)境溫度為303 K。

2.3 雙向熱流固耦合仿真設(shè)定

在Ansys workbench中將Transient structural模塊與Fluent模塊用System coupling模塊連接，分別在Transient structural與Fluent中導(dǎo)入固體網(wǎng)格與流體網(wǎng)格，并設(shè)置其邊界條件；Fluent中湍流模型選擇標(biāo)準(zhǔn)k-ε模型并勾選黏性發(fā)熱，導(dǎo)入自編譯UDF文件以模擬Fluent軟件中黏度隨溫度的變化情況，采用由實驗數(shù)據(jù)擬合后的函數(shù)關(guān)系式(5)。油膜初始偏心率為0.4，轉(zhuǎn)速分別設(shè)置1 500、2 000、2 500 r/min。

在Ansys Workbench中的Transient structural模塊輸入apdl命令以開啟多物理場耦合單元中的熱固耦合單元，利用System coupling模塊將其與Fluent中的流體域進行耦合，實現(xiàn)基于Ansys Workbench的雙向熱流固耦合。

3 有限元模型的驗證

3.1 Fluent仿真及MATLAB模型驗證

Fluent與MATLAB受黏溫影響的壓力計算結(jié)果如圖5所示，其中Fluent仿真結(jié)果平面展開圖基于Fluent仿真數(shù)據(jù)導(dǎo)出后由MATLAB程序重新繪制實現(xiàn)；另外，因MATLAB采用雷諾邊界，沒有負(fù)壓區(qū)，故對Fluent仿真結(jié)果負(fù)壓區(qū)進行了處理。

從圖5可以看出，F(xiàn)luent和MATLAB仿真計算的最大油膜壓力基本相同，誤差在2%以內(nèi)，壓力變化趨勢基本一致。Fluent仿真結(jié)果中，壓力峰值前的波動為油液進口處的壓力沖擊造成的；MATLAB仿真時，由于沒有類似實際油膜的壓力進口，使得計算結(jié)果存在誤差。Fluent計算結(jié)果表明，隨著轉(zhuǎn)速的增加，油膜最大壓力也會相應(yīng)增加。通過MATLAB模型驗證，F(xiàn)luent仿真結(jié)果具有可靠性，為下文的雙向熱流固耦合計算及驗證提供依據(jù)。

圖5 不同轉(zhuǎn)速下Fluent和MATLAB仿真結(jié)果對比

3.2 雙向熱流固耦合仿真結(jié)果驗證

根據(jù)上述雙向熱流固耦合計算邊界條件，分別在轉(zhuǎn)速1 500、2 000、2 500 r/min下基于MATLAB模型對雙向熱流固耦合模型計算結(jié)果進行驗證。熱流固耦合與MATLAB受黏溫影響的壓力計算結(jié)果如圖6、7所示。圖7中針對未考慮黏溫影響的熱流固計算結(jié)果驗證時，將MATLAB模型中黏溫關(guān)系式改為常數(shù)，由于熱流固耦合模型計算考慮了固體域形變對油膜的影響，因此造成了兩模型計算結(jié)果的差異，但固體域形變很小，故影響甚微，兩模型整體油膜壓力變化趨勢基本一致，驗證了熱流固耦合模型計算結(jié)果具有良好的可靠性。

圖6 考慮黏溫影響時不同轉(zhuǎn)速下Fluent和MATLAB仿真壓力對比

圖7 不考慮黏溫影響時不同轉(zhuǎn)速下Fluent和MATLAB仿真壓力對比

4 雙向熱流固耦合仿真及分析

4.1 忽略黏溫效應(yīng)與考慮黏溫效應(yīng)的雙向熱流固耦合計算結(jié)果及分析

如圖8所示分別為1 500、2 000、2 500 r/min轉(zhuǎn)速下油膜及主軸溫度分布。當(dāng)其他條件相同時[14]，油膜溫度從中間向兩端逐漸升高；主軸和油膜溫度隨主軸轉(zhuǎn)速的增加而升高，呈對稱分布；最高油膜溫度處于偏心位置的邊緣處。

圖8 不同轉(zhuǎn)速下油膜及主軸溫度分布

如圖9所示分別為1 500、2 000、2 500 r/min轉(zhuǎn)速下忽略與考慮黏溫效應(yīng)的油膜壓力分布?？梢钥闯觯?dāng)其他條件相同時，主軸油膜的壓力隨著主軸轉(zhuǎn)速的增加而提高，同轉(zhuǎn)速下受黏溫效應(yīng)影響壓力會下降。3種轉(zhuǎn)速下考慮黏溫效應(yīng)影響的油膜壓力較之于忽略黏溫效應(yīng)的油膜壓力，下降幅度分別為25.66%、27.56%、29.56%，隨主軸轉(zhuǎn)速增加，油膜壓力受黏溫效應(yīng)的影響愈發(fā)顯著。

圖9 不同轉(zhuǎn)速下忽略與考慮黏溫效應(yīng)油膜壓力分布

如圖10所示為3種轉(zhuǎn)速下忽略與考慮黏溫效應(yīng)的主軸形變情況，總形變?yōu)闊釕?yīng)變與油膜壓力共同作用下的結(jié)果。

圖10 不同轉(zhuǎn)速下忽略與考慮黏溫效應(yīng)的主軸形變

由圖10可見，主軸形變隨主軸轉(zhuǎn)速的增加略微升高；考慮黏溫效應(yīng)時，由于主軸油膜壓力下降，主軸形變會略微降低；3種轉(zhuǎn)速下考慮黏溫效應(yīng)的主軸形變較之于忽略黏溫效應(yīng)的主軸形變，下降幅度分別為0.10%、0.16%、0.26%。故考慮黏溫效應(yīng)時，主軸轉(zhuǎn)速變化對主軸形變的影響很小。

4.2 忽略與考慮黏溫效應(yīng)時軸頸轉(zhuǎn)速對油膜特性的影響

圖11所示為忽略與考慮黏溫效應(yīng)時不同轉(zhuǎn)速下的油膜特性分布?？梢园l(fā)現(xiàn)，當(dāng)軸頸轉(zhuǎn)速不斷增加時，主軸內(nèi)部油膜的承載力也會隨之增加，這是因為隨轉(zhuǎn)速提升，柱塞傳遞給座圈的三面壓力也會提高，且轉(zhuǎn)速越大，動壓效果越明顯，承載能力越好[15]。承載力隨主軸轉(zhuǎn)速增加的關(guān)系近似為一條直線，且黏度恒定的油膜承載力要遠(yuǎn)大于考慮黏溫效應(yīng)的油膜承載力。這是因為，受黏度影響，潤滑油不斷消耗軸頸旋轉(zhuǎn)所提供的機械能，不斷轉(zhuǎn)化為熱量使油膜溫度升高從而降低其黏度，導(dǎo)致潤滑油流速加快從而降低了油膜的承載力；且轉(zhuǎn)速越高，承載力降低越明顯。

圖11 忽略與考慮黏溫效應(yīng)時最大油膜承載力隨轉(zhuǎn)速的變化

圖12所示為忽略與考慮黏溫效應(yīng)時最小油膜厚度隨轉(zhuǎn)速的變化。相較于油膜的承載力，油膜最小厚度隨主軸轉(zhuǎn)速變化較小，隨主軸轉(zhuǎn)速增加油膜最小厚度略微增加。最小油膜厚度受黏溫影響略微減小，油膜由于受黏溫影響所產(chǎn)生的黏性剪切應(yīng)力減小，同油膜厚度下，黏度小的油膜所產(chǎn)生的油膜壓力也較小，從而最小油膜厚度減小。

圖12 忽略與考慮黏溫效應(yīng)時最小油膜厚度隨轉(zhuǎn)速的變化

圖11、12的結(jié)果表明，不考慮柴油黏溫特性的油膜承載力、最小油膜厚度均比考慮柴油黏溫效應(yīng)時高，而主軸在實際的工作過程中柴油的黏溫效應(yīng)是真實存在的，如果在主軸的設(shè)計環(huán)節(jié)忽略黏溫效應(yīng)，設(shè)計出的主軸軸徑承載力將達不到預(yù)期效果，主軸實際工作中摩擦副的潤滑形式處于混合潤滑與邊界潤滑的概率增大，干摩擦概率增大，從而導(dǎo)致主軸損傷。

5 結(jié)論

(1)提出一種可代入溫度實驗數(shù)據(jù)或?qū)胲浖嬎銣囟葦?shù)據(jù)的黏溫效應(yīng)MATLAB一維計算仿真模型。與傳統(tǒng)簡化后的能量方程推導(dǎo)方法相比，該方法包含表面粗糙度等復(fù)雜因素對溫度的影響，能更好地反映油膜實際溫度，使計算結(jié)果與三維仿真結(jié)果有更好的對應(yīng)?？紤]黏溫影響時采用MATLAB的油膜壓力仿真結(jié)果與采用Fluent的油膜壓力仿真結(jié)果變化趨勢一致，最大相對誤差2%，可作為考慮黏溫影響的Fluent及油膜-主軸雙向熱流固耦合仿真結(jié)果的一種有效驗證手段。

(2)考慮黏溫效應(yīng)影響時，隨著轉(zhuǎn)速增加，油膜壓力、溫度和承載力有顯著增加，主軸的形變量和最小油膜厚度變化較小。轉(zhuǎn)速增加導(dǎo)致的油膜溫度升高會使?jié)櫥宛ざ炔痪鶆颍M而導(dǎo)致主軸系統(tǒng)的失穩(wěn)。因此，在主軸設(shè)計時，需設(shè)置主軸合理的工作轉(zhuǎn)速以控制最高溫升在允許范圍內(nèi)，保證主軸的熱穩(wěn)定性。

(3)隨著主軸轉(zhuǎn)速增加，黏溫效應(yīng)對油膜壓力及承載力影響越趨明顯，忽略黏溫效應(yīng)的油膜壓力及承載力比考慮黏溫效應(yīng)影響的油膜壓力及承載力普遍偏高，且隨著轉(zhuǎn)速增加，差值將更加明顯。故高轉(zhuǎn)速下油膜特性分析時，需要考慮黏溫效應(yīng)的影響，否則分析結(jié)果會偏高。如果在主軸的設(shè)計環(huán)節(jié)忽略黏溫效應(yīng)，設(shè)計出的主軸承載力將達不到預(yù)期效果，主軸在實際工作中發(fā)生干摩擦概率將增大，主軸易因承載能力不足而損傷。