核電循環泵行星齒輪熱彈流潤滑研究

董慶兵，汪漢明，馮成程，魏靜，何東，楊海師

(1.重慶大學機械傳動國家重點實驗室，重慶 400044；2.重慶大學機械與運載工程學院， 重慶 400044；3.重慶齒輪箱有限責任公司，重慶 402263 )

核電循環水泵是沿海核電站重要設備，其作用是向汽輪機的凝汽器和輔助冷卻水系統提供冷卻水，降低核反應堆運行溫度，避免高溫熔毀核反應裝備，其可靠性直接影響到核電機組的安全運行。核電齒輪箱的良好潤滑性能是核電循環泵可靠運行的重要保障，而潤滑不良會給核電機組帶來巨大的經濟損傷和安全隱患。據不完全統計，因潤滑不良導致的設備事故約占25%，與潤滑因素相關的設備故障高達40%[1]。因此，研究提升核電齒輪箱的潤滑性能具有重要意義。然而，如何減小齒面潤滑不良造成的接觸疲勞，仍是當前提升核電齒輪箱服役性能關注的熱點問題。

齒面潤滑是集熱-固-流等多個物理場的典型彈流問題，國內外學者對齒輪的熱彈流潤滑問題已展開了廣泛研究。早期漸開線齒輪嚙合中的彈性流體動力潤滑(EHL)研究模型大多集中在直齒輪，且假設接觸表面光滑的情況[2-4]。隨著齒輪精度和平穩性要求的逐漸提高，斜齒輪的應用更加廣泛，且粗糙表面對齒輪彈流潤滑的影響已不容忽視。斜齒輪在嚙合過程中嚙合線長度時刻變化，摩擦力沿著接觸線方向也不統一，相較于直齒輪的嚙合特性和潤滑性能更復雜多變。莫云輝和巢桐[5]提出有限寬線接觸正反圓錐滾子模型作為斜齒圓柱齒輪的等效模型。楊萍和楊沛然[6]研究了斜齒圓柱齒輪瞬時的準穩態熱彈流潤滑行為，分析了端部修形對潤滑參數的影響。CLARKE等[7]根據斜齒輪嚙合瞬態彈流潤滑數值結果，分析了輪廓誤差對斜齒輪潤滑性能的影響。YANG等[8]考慮真實三維粗糙表面，分析了斜齒輪嚙合過程中的混合EHL特性、摩擦力和閃溫。上述研究考慮了齒輪修形、粗糙表面等特性對其外嚙合斜齒輪彈流潤滑結果的影響，但綜合因素下內/外嚙合的彈流潤滑行為有待深入研究。

本文作者以核電循環泵人字齒行星傳動系統為研究對象，忽略平行斜齒輪相互的軸向作用，考慮齒廓修形造成的嚙合誤差，以及齒面磨合期形貌的塑性變形等因素，利用移動平均濾波方法對未經磨合的初始形貌進行光滑處理，建立了典型工況下核電循環泵行星傳動系統斜齒輪熱彈流潤滑模型，研究了內/外嚙合齒輪副接觸區域內油膜厚度、油膜壓力、摩擦應力、閃溫等潤滑參數的變化規律，分析了齒面形貌和齒廓修形對其潤滑行為的影響，為核電循環泵齒輪箱的運行可靠性和接觸性能優化提供了技術手段和理論基礎。

1 行星齒輪三維結構及齒面接觸分析

1.1 行星輪系結構參數

核電齒輪箱通常采用立式NGW結構、功率四分流人字齒行星輪系一級傳動，如圖1所示，由于行星輪系的人字齒為對稱結構，因此對單側斜齒輪進行彈流潤滑計算。

圖1 核電循環泵齒輪箱三維結構Fig.1 Three dimensional structure of gear box of nuclear power circulating pump

齒輪的主要參數如表1所示，其中齒寬表示人字齒兩側齒寬及中間軸長度之和。

由于裝配誤差、制造誤差及齒輪受載后彈性變形等因素的影響，齒輪嚙合過程中會出現偏載和嚙合沖擊，工程應用中主要采用齒廓修形來減小嚙合沖擊和應力集中。

文中針對核電循環泵齒輪采取齒頂線性修形，其端面示意圖如圖2所示。圖中：Ca表示最大修形量，La表示修形長度。線性修形指的是在齒輪嚙合線上，任意修形位置x與修形量Δ之間呈線性關系，具體公式如式(1)所示?？紤]實際加工難度、經濟成本等因素，一般只對太陽輪或行星輪進行修形。文中按照表2所示參數對太陽輪和行星輪進行齒頂線性修形，在Fortran程序中可更改這些修形參數，從而求解不同修形量下彈流潤滑數值結果。

圖2 齒頂修形端面示意Fig.2 Schematic of tooth tip modification end face

表2 行星齒輪齒頂修形參數Tab.2 Planetary gear tooth tip modification parameters

Δ=Ca(x/La)

(1)

1.2 行星齒輪內/外嚙合齒面接觸分析

斜齒輪齒面接觸線的長度在嚙合過程中是時刻變化的，考慮齒輪重合度的影響，實際斜齒輪在嚙合時同時存在多條時變接觸線。文中模型在齒輪嚙合3/4齒交替的齒面上每條嚙合線都具有時變性。文中采用切片法，將斜齒輪沿齒寬方向劃分成若干份斜齒輪薄片，對單個齒輪薄片的剛度進行累加，從而得到斜齒輪副的嚙合剛度。實際傳動中輪齒嚙合剛度受到多種因素影響，文中考慮斜齒輪接觸變形、齒廓偏差、輪體變形以及軸向力的影響，得到輪齒誤差影響下的多齒綜合嚙合剛度表達式[9]如下：

(2)

基于上述綜合嚙合剛度可推導出齒輪薄片j的載荷分配系數表達式為

(3)

斜齒輪嚙合副的瞬時嚙合狀態熱彈流潤滑問題可以近似等效為圓錐滾子的準穩態熱彈流潤滑問題，如圖3所示。內嚙合副可等效為兩同向的圓錐滾子內切接觸問題，即內嚙合線上某點與齒圈基圓相切的曲率半徑對應半徑R1，與行星輪基圓相切的曲率半徑對應半徑R2。外嚙合副可等效為兩反向的圓錐滾子外切接觸問題，即外嚙合線上某點與太陽輪基圓相切的曲率半徑對應半徑R3，與行星輪基圓相切的曲率半徑對應半徑R4。

由圖3可得行星輪與內齒圈的曲率半徑R1、R2表達式為

圖3 行星輪內/外嚙合等效示意Fig.3 Schematic of internal and external meshing pairs of planetary gear

(4)

同理，太陽輪與行星輪外嚙合的曲率半徑R3、R4由以下公式得到：

(5)

式中：a1和a2表示內/外嚙合副中心距；rb1、rb2、rb3分別表示行星輪、齒圈、太陽輪的基圓半徑。

以內嚙合為例，可以確定y方向的等效曲率半徑Ry為

(6)

求得曲率半徑R1、R2后，可以得到嚙合點大、小齒輪線速度u1、u2為

u1=ω1R1，u2=ω2R2

(7)

卷吸速度ue、滑動速度us和滑滾比ξ的定義如下：

ue=(u1+u2)/2，us=u1-u2，ξ=us/ue

(8)

式中：下標“1、2”表示行星輪和齒圈。

2 彈流潤滑理論

2.1 潤滑模型基本方程

文中采用彈流潤滑有限長線接觸模型，每一嚙合瞬時等效為兩反向/同向圓錐滾子的接觸，熱彈流潤滑控制方程如下：

(9)

式中：p表示流體壓力；h表示油膜厚度；ux和uy表示沿x和y方向的入口速度；ρ、η分別表示潤滑劑的密度和黏度。

膜厚方程：

h(x，y)=h0+hs(x，y)+v(x，y)

(10)

式中：h0表示兩接觸體的剛性位移；hs表示接觸表面的原始幾何間隙；v表示齒面彈性變形。

hs(x，y)=ha+hb+hc

(11)

(12)

hb計算同理。

式中：ha、hb為齒輪a、b的接觸表面的原始幾何間隙；hc為齒頂修形引起的幾何間隙；Ra、Rb為齒輪a、b的曲率半徑。

由壓力引起的表面彈性變形的解為

(13)

載荷平衡方程為

(14)

熱彈流潤滑考慮了溫度對密度、黏度等參數的影響，在溫升較大的極端工況條件下，建模時需在每次求解溫度后更新這些參數，從而準確描述相應工況下的潤滑行為。但針對溫升較小的工況條件可忽略溫度對這些參數的影響。PU等[10-11]研究了橢圓點接觸任意卷吸方向混合潤滑條件下的摩擦行為和閃溫分布，并實驗驗證了摩擦因數計算的準確性；隨后又采用同樣的處理方法，進一步研究了供油不充分的乏油潤滑條件下的溫升。

文中以核電循環泵人字齒行星傳動系統為研究對象，設計合理的齒輪嚙合時溫升通常較小(≤100 K)，由此造成的溫度和黏度變化可忽略不計。盡管忽略了溫升的影響，文中模型仍可相對準確獲得齒輪的潤滑行為。文中密壓關系采用Dowson-Higginson公式：

(15)

黏度遵循Roelands黏度-壓力關系公式：

η(x，y)=η0exp[(lnη0+9.67)(1+p/(1.96×108))z-1]

(16)

指數z由以下公式得到：

(17)

式中：ρ0、η0分別為潤滑劑環境密度和環境黏度；α為黏壓系數。

2.2 摩擦力和閃溫模型

滑動摩擦和溫度上升主要產生在赫茲接觸區域，潤滑膜剪切應力是溫度的函數，溫度上升是由表面滑動摩擦引起的，混合潤滑時還需考慮固-固直接接觸產生的摩擦力。采用考慮溫升效應的非牛頓流體模型計算潤滑油膜中的剪應力：

(18)

式中，極限剪應力τL和極限剪切彈性模量G∞是壓力和溫度的函數，可以通過經驗估計。對于典型的礦物油，可使用以下經驗公式：

(19)

τL(p，T)=0.25G∞

(20)

假設潤滑油的剪切應變率可以按下式近似地計算出：

(21)

于是有：

(22)

通過迭代方法求解計算網格中每個節點的上述非線性方程，可以求解油膜中剪應力τ的分布。由于潤滑劑的特性對溫度敏感，而溫度直接受到摩擦產生熱量的影響，因此表面溫度和摩擦是相互影響的。文中根據半無限固體上的運動熱源理論，忽略沿表面垂直于熱源速度方向的熱流，以及潤滑劑的對流換熱作用?；谏鲜黾僭O，可通過求解以下方程式計算齒輪表面閃溫[12]。

(23)

(24)

式中：q是潤滑區域的剪切或界面摩擦產生的熱量，它可以被估計為q?τ|u2-u1|；Tb1和Tb2為兩接觸表面初始溫度；ρ1和ρ2為密度；C1和C2為比熱容；kf為熱傳導系數；其中剪應力τ來自公式(20)。

齒輪表面閃溫計算參數如表3所示。

表3 潤滑劑和齒輪的材料性能

2.3 數值方法

董慶兵團隊[12-13]采用量綱一化方法建立了彈流潤滑模型，研究了瞬態工況下潤滑油膜的動態變化規律，為核電循環泵斜齒輪傳動系統的潤滑仿真提供了理論依據。文中研究對象為斜齒輪，在嚙合過程中接觸線長度具有時變性，齒輪和坐標系建立方法的差異性會具有不同的計算區域，因此采取具體邊界數值，將齒輪每個嚙合瞬時沿齒廓切線方向設為x軸，接觸線沿齒寬方向設為y軸，接觸線中心點為坐標原點，則計算區域如下：

{(x，y)|xs≤x≤xe，ys≤y≤ye}

(25)

為保證不同嚙合瞬時的接觸區域都被包含在計算范圍內，外嚙合計算區域邊界為(xs=-0.8 mm，xe=0.5 mm；ys=-0.5b/cosβb，ye=0.5b/cosβb)，內嚙合計算區域邊界為(xs=-2 mm，xe=1 mm；ys=-0.5b/cosβb，ye=0.5b/cosβb)，b表示齒寬；βb為基圓螺旋角。將計算域沿x軸方向劃分為256個均勻網格，沿y軸方向劃分為512個均勻網格。運用二階中心差分格式計算雷諾方程的左端項，運用二階后差分格式計算公式右端楔形項，利用線松弛迭代法和逐行掃描法求解油膜壓力，齒面的彈性變形采用快速傅里葉變換算法求解。壓力計算的迭代因子ωp取0.01，壓力收斂精度取1×10-5。

3 潤滑特性計算結果及分析

3.1 接觸線各節點運動參數求解結果

根據人字齒行星傳動系統的幾何參數和計算公式，可以求得單側齒(即斜齒輪)嚙合過程中每一瞬時嚙合線中點的卷吸速度ue、滑動速度us和滑滾比ζ等時變參數的變化規律，如圖4所示。內/外嚙合副一個齒面從進入嚙合到嚙合結束的過程中，嚙合線中點處滑滾比先線性增加，然后保持不變，隨后繼續線性增加，滑動速度與卷吸速度呈現同樣的變化趨勢，并且內/外嚙合副接觸線中心點滑滾比為0的時間段不一致。

圖4 嚙合線中心節點運動參數Fig.4 Motion parameters of meshing line center node： (a) for internal mesh；(b) for external mesh

圖5所示為內/外嚙合線各點在一個周期內的滑滾比云圖，其縱坐標表示嚙合線方向上，嚙合線始端到末端所占用網格點數。從圖中可以發現嚙合線上各點的滑滾比都是在不斷變化的，以滑滾比為0的等高線作為分界點，分別選取內/外嚙合A、B、C、D 4處不同時刻的計算結果進行對比分析。其中B位置表示內/外嚙合節點附近，此時嚙合線中點處滑滾比為0。

圖5 嚙合副接觸域滑滾比分布云圖Fig.5 Distribution nephogram of slip-roll ratio in contact region of meshing pair：(a) for internal mesh；(b) for external mesh

3.2 穩態彈流潤滑數值求解結果

圖6所示為內嚙合B、D、C 3個位置的油膜壓力和油膜厚度在嚙合線中點x、y方向上的變化曲線。從圖6(a)中可以看出，沿x方向上，嚙合前期B位置處，行星輪與齒圈接觸區域前期油膜壓力較大，油膜壓力整體呈現先增加后減小的趨勢。由于二次壓力峰的存在，在出口區域的油膜形狀收縮，油膜厚度會局部減小，如圖6(b)所示。沿y方向上，內嚙合嚙合線長度先增加后減小，油膜壓力與油膜厚度變化呈相反趨勢的變化規律，且嚙合線兩端存在應力集中，如圖6(c)(d)所示。這是因為內嚙合的等效曲率半徑逐漸增加，導致沿著嚙合線方向的油膜壓力和油膜厚度會呈現向一端的傾斜趨勢。

圖6 內嚙合彈流潤滑結果Fig.6 EHL results of internal engagement：(a)oil film pressure in x direction；(b)oil film thickness in x direction； (c)oil film pressure in y direction；(d)oil film thickness in y direction

圖7所示為外嚙合A、B、C 3個位置的油膜壓力和油膜厚度在嚙合線中點x、y方向上的變化曲線。從圖7(a)中可以看出，沿x方向上，外嚙合前期嚙合線中點處油膜壓力較小，隨著齒輪轉動載荷增大導致中心壓力增加，二次壓力峰向出口方向偏移。油膜厚度在接觸區域中心處無明顯變化，但在出口區域二次壓力峰仍會導致油膜厚度局部減小，如圖7(b)所示，這是由于等效曲率半徑、載荷、嚙合線長度變化的綜合影響所導致。沿y方向上，由于外嚙合的等效曲率半徑先增加后減小，沿y方向的油膜壓力和油膜厚度則會呈現對稱分布，嚙合線兩端仍存在應力集中，如圖7(c)(d)所示。

圖7 外嚙合彈流潤滑結果Fig.7 EHL results of external engagement：(a)oil film pressure in x direction；(b)oil film thickness in x direction； (c)oil film pressure in y direction；(d)oil film thickness in y direction

3.3 齒輪表面形貌對彈流潤滑性能的影響

在齒輪的實際加工過程中，齒輪表面的粗糙度數值與彈流潤滑分析中的膜厚在同一數量級上，因此齒面形貌不能被忽略。文中建立的彈流潤滑計算模型同樣適用于粗糙表面，在公式(11)中加入實驗測量得到的幾何間隙，從而得到考慮齒輪表面形貌后的幾何間隙，再對原始齒輪表面的粗糙度經過磨合后所形成的形貌數據進行分析處理。原始粗糙度數據采用重慶大學機械傳動國家重點實驗室的高集成多功能摩擦磨損試驗機(MFT-5000)，使用白光干涉儀實測6級加工精度齒表面得到。

文中考慮齒輪嚙合過程中的磨合效應，采用移動平均濾波方法使表面變得平滑[14]，采用2D過濾方法，在x、y方向上對粗糙度數據采取固定點數進行平均濾波；測量表面粗糙度時，繪制成Abbot曲線[15]，用于表示磨合前后的粗糙表面變化情況[16]。濾波后的表面粗糙形貌最大值降低了30%左右，采用移動平均濾波方法過濾后的表面形貌可以代表齒輪表面磨合后的形貌進行彈流潤滑計算。

以嚙合線中心為坐標原點，中心油膜厚度為計算域的坐標原點，計算得到中心油膜厚度、摩擦因數、太陽輪和行星輪嚙合線的中心溫升以及齒面最大溫升在一個周期內的變化曲線，如圖8所示。對于考慮表面形貌的彈流潤滑分析，中心膜厚在相應位置產生波動的同時，比光滑表面有所減小，如圖8(a)所示，從而使得潤滑效果變差。摩擦因數是整個計算區域所有節點剪切應力的積分和與載荷的比值，在嚙合線中點滑滾比為0的位置，摩擦因數趨近于0；而在嚙合中期波峰與其對應的波谷積分求和的值與光滑表面積分求和數值相差很小，因此摩擦因數波動量較小，如圖8(b)所示。

圖8(c)(d)所示為太陽輪和行星輪的中心溫升以及齒面最大溫升。在考慮粗糙度時，油膜厚度減小，剪切應力會增加，但中間位置的摩擦因數變化較小，因此中心溫升變化與摩擦因數分布特點相似，各個瞬時的中心溫升在剪切應力為0的位置依然為0，在其他嚙合位置略微升高。而摩擦因數的整體增加，導致剪切應力增大產生更多熱量，油膜溫度快速升高，齒面最大溫升會明顯增大。因此，在彈流潤滑分析中考慮齒面形貌是十分必要的。

圖8 考慮表面形貌影響的彈流潤滑數值解Fig.8 Numerical solution of EHL considering the influence of surface topography：(a) central oil film thickness；(b) friction coefficient；(c) temperature rise of sun gear surface；(d)temperature rise of planetary gear surface

3.4 齒廓修形對彈流潤滑性能的影響

由于未修形齒輪的嚙合線兩端存在應力集中，齒廓修形可有效改善其嚙合狀態，因此有必要分析齒廓修形對其潤滑性能的影響。在一個嚙合周期的不同位置，參與嚙合的修形端是不斷變化的，因此齒頂修形對彈流潤滑求解的影響也是不同的。按照表2所示參數對太陽輪和行星輪進行齒頂線性修形，計算內/外嚙合在B位置修形前后油膜壓力在y方向的結果，如圖9所示。若嚙合過程中修形端未參與嚙合，則嚙合線長度減小，載荷會有所增加，導致油膜壓力增加，如內嚙合B位置，見圖9(a)；若修形端兩側都參與嚙合，油膜壓力兩側都會降低，如外嚙合B位置，見圖9(b)。

圖10所示為外嚙合副在B位置嚙合線計算域的剪切應力以及齒面溫升。圖10(a)所示為y=40 mm處的齒面溫升在x方向的分布曲線，T1表示太陽輪表面溫升。同方向的剪切應力和齒面溫升分布特點相似，修形后的溫升和剪切力在x方向都會稍微增加，齒面溫升出口大于0。這是因為處于邊界區域時太陽輪和行星輪接近環境溫度，中間過程通過剪切運動產生熱量，導致齒面溫度升高；而由于二次壓力峰的存在，出口油膜剪切力比入口大，產生更多的摩擦熱，溫度降低緩慢，所以出口溫度高于入口溫度，但x方向的剪切力入口和出口邊界都為0。圖10(b)所示為y方向剪切應力變化曲線，根據修形結果圖9(b)所示，嚙合線兩端應力集中減小，則剪切應力也隨之減小，齒面溫升也會有所下降。因此，在彈流潤滑分析中考慮齒輪修形是十分必要的。

圖9 修形對嚙合位置B油膜壓力的影響Fig.9 Influence of profile modification on oil film pressure at position B：(a)for internal mesh；(b)for external mesh

圖10 齒頂修形對剪切應力和齒面溫升的影響Fig.10 Influence of tooth tip modification on shear stress：(a) xdirection tooth surface temperature rise；(b) y direction shear stress

4 結論

針對行星輪內/外嚙合齒輪副，考慮非勻布的載荷分配特性解析得到內/外斜齒輪嚙合過程中接觸區域的運動參數變化規律，將行星輪齒輪副嚙合狀態幾何等效為圓錐滾子的接觸問題，建立典型工況下斜齒輪線接觸熱彈流潤滑計算模型，分析了計算域中各節點的油膜厚度、油膜壓力、剪切應力和溫升分布情況，研究了輪齒表面形貌和齒廓修形對彈流潤滑性能的影響。主要結論如下：

(1) 內/外嚙合副在嚙合過程中的等效曲率半徑、卷吸速度和滑滾比等參數的變化有所不同。在嚙合前期和后期，內/外嚙合的嚙合線長度較短，油膜壓力較小，接觸線上不存在滑滾比為0的點。在嚙合中期，油膜壓力較大，嚙合線兩端存在應力集中，膜厚較小，內/外嚙合接觸線上存在滑滾比為0的點，剪切應力在嚙合線上呈現正負交替變化。

(2) 考慮齒面粗糙形貌時，中心油膜厚度比光滑解減小，粗糙度的存在弱化了油膜承載能力。在輪齒嚙合前期和后期，摩擦因數增大比較明顯，中心摩擦因數波動較??；考慮粗糙度的齒表面最大溫升有明顯的增大趨勢，這與接觸壓力的增加導致剪切應力增大有關。因此，在熱彈流潤滑模型的精密計算中考慮粗糙度的影響是十分必要的。

(3) 對太陽輪和行星輪進行齒廓修形，研究發現油膜壓力和油膜厚度減小與參與的嚙合修形端數量有關；修形后剪切應力與齒面溫升變化規律相似，但兩者受到摩擦因數、修形量和載荷分配等多種因素的綜合影響；合理的修形可以有效降低嚙合端的應力集中，減小油膜壓力、油膜厚度、剪切應力和齒面溫升，改善嚙合線終端的潤滑狀態。