許 磊 張 翼 徐春龍 李斌茂 馬運(yùn)昌 唐詩(shī)澤 張 宇

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

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

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

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

1 MATLAB數(shù)值模擬

1.1 方程的量綱一化

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

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

(1)

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

對(duì)式(1)進(jìn)行量綱一化。采用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為量綱一油膜厚度。對(duì)于壓力p的量綱一化形式，先以一未定的p0值作為p的量綱一化數(shù)，則有：P=p/p0。

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

(2)

(3)

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

U=2πnr，由此可得參考?jí)毫0隨轉(zhuǎn)速升高而增加。

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

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

圖1 柴油黏溫特性曲線

(4)

其中計(jì)算黏度變化時(shí)采用基于柴油黏溫特性實(shí)驗(yàn)數(shù)據(jù)[11]經(jīng)由MATLAB擬合工具箱進(jìn)行函數(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è)定周向和軸向的步長(zhǎng)并為其劃分網(wǎng)格，將 Reynolds方程轉(zhuǎn)化為一組差分方程以此求得各節(jié)點(diǎn)的壓力值。i方向的步長(zhǎng)為Δφ，j方向的步長(zhǎng)為Δλ。按差分原理，將油膜厚度H沿周向i的一階差商等效為油膜厚度H在周向i的導(dǎo)數(shù)，如式(6)所示；將油膜壓力P和黏度的差分形式等效為2個(gè)方向上的偏導(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)

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

(14)

超松弛迭代公式為

Pi,j=

(15)

其中取α=1.3。

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

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é)點(diǎn)數(shù)496個(gè)，軸向節(jié)點(diǎn)數(shù)120個(gè)，網(wǎng)格總體質(zhì)量達(dá)到0.75以上。經(jīng)計(jì)算檢驗(yàn)，當(dāng)計(jì)算域網(wǎng)格總數(shù)為40萬左右時(shí)，油膜最大壓力、有效載荷計(jì)算值精度較好。

圖2 凸輪軸仿真模型示意

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

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

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

表1 油膜邊界條件

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

主軸材料為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è)主軸與周圍的空氣傳熱為自然對(duì)流，復(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í)驗(yàn)數(shù)據(jù)擬合后的函數(shù)關(guān)系式(5)。油膜初始偏心率為0.4，轉(zhuǎn)速分別設(shè)置1 500、2 000、2 500 r/min。

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

3 有限元模型的驗(yàn)證

3.1 Fluent仿真及MATLAB模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

如圖9所示分別為1 500、2 000、2 500 r/min轉(zhuǎn)速下忽略與考慮黏溫效應(yīng)的油膜壓力分布?？梢钥闯?，當(dāng)其他條件相同時(shí)，主軸油膜的壓力隨著主軸轉(zhuǎn)速的增加而提高，同轉(zhuǎn)速下受黏溫效應(yīng)影響壓力會(huì)下降。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)時(shí)，由于主軸油膜壓力下降，主軸形變會(huì)略微降低；3種轉(zhuǎn)速下考慮黏溫效應(yīng)的主軸形變較之于忽略黏溫效應(yīng)的主軸形變，下降幅度分別為0.10%、0.16%、0.26%。故考慮黏溫效應(yīng)時(shí)，主軸轉(zhuǎn)速變化對(duì)主軸形變的影響很小。

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

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

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

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

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

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

5 結(jié)論

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

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

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