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

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

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

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

綜上所述，學者們在結合能量方程求解油膜溫度時，由于受多因素的影響能量方程過于復雜，均會對能量方程進行簡化，而簡化后的能量方程未必能反映真實的油膜溫度?；诖?，本文作者通過數(shù)學推導，建立一種可代入溫度實驗數(shù)據或導入軟件計算溫度數(shù)據的黏溫效應一維計算仿真模型，以更好地反映真實的油膜溫度；考慮到主軸內部油膜的熱量及溫度受多因素復雜耦合作用，主軸真實運轉時油膜特性受主軸形變影響，建立Workbench平臺下雙向熱流固耦合模型，研究黏溫效應對主軸摩擦副油膜特性的影響；對黏溫效應影響下柴油潤滑主軸的油膜特性進行分析，為主軸設計提供依據。

1 MATLAB數(shù)值模擬

1.1 方程的量綱一化

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

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

(1)

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

對式(1)進行量綱一化。采用x的量綱一化形式為：x=rφ→φ=x/r，0≤φ≤2π，其中φ為由軸與軸承中心連線沿軸旋轉方向角度；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)

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

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

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

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

圖1 柴油黏溫特性曲線

(4)

其中計算黏度變化時采用基于柴油黏溫特性實驗數(shù)據[11]經由MATLAB擬合工具箱進行函數(shù)擬合后的函數(shù)關系式(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方程的差分形式

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

(6)

(7)

(8)

(9)

(10)

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

(11)

(12)

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

(13)

其中令：

E=A+B+C+D

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

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

(14)

超松弛迭代公式為

Pi,j=

(15)

其中取α=1.3。

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

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

2.1 模型與網格

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

圖2 凸輪軸仿真模型示意

圖3 油膜網格示意

2.2 邊界條件的設定

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

表1 油膜邊界條件

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

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

2.3 雙向熱流固耦合仿真設定

在Ansys workbench中將Transient structural模塊與Fluent模塊用System coupling模塊連接，分別在Transient structural與Fluent中導入固體網格與流體網格，并設置其邊界條件；Fluent中湍流模型選擇標準k-ε模型并勾選黏性發(fā)熱，導入自編譯UDF文件以模擬Fluent軟件中黏度隨溫度的變化情況，采用由實驗數(shù)據擬合后的函數(shù)關系式(5)。油膜初始偏心率為0.4，轉速分別設置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受黏溫影響的壓力計算結果如圖5所示，其中Fluent仿真結果平面展開圖基于Fluent仿真數(shù)據導出后由MATLAB程序重新繪制實現(xiàn)；另外，因MATLAB采用雷諾邊界，沒有負壓區(qū)，故對Fluent仿真結果負壓區(qū)進行了處理。

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

圖5 不同轉速下Fluent和MATLAB仿真結果對比

3.2 雙向熱流固耦合仿真結果驗證

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

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

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

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

4.1 忽略黏溫效應與考慮黏溫效應的雙向熱流固耦合計算結果及分析

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

圖8 不同轉速下油膜及主軸溫度分布

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

圖9 不同轉速下忽略與考慮黏溫效應油膜壓力分布

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

圖10 不同轉速下忽略與考慮黏溫效應的主軸形變

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

4.2 忽略與考慮黏溫效應時軸頸轉速對油膜特性的影響

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

圖11 忽略與考慮黏溫效應時最大油膜承載力隨轉速的變化

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

圖12 忽略與考慮黏溫效應時最小油膜厚度隨轉速的變化

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

5 結論

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

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

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