摘" 要：該文基于多重網(wǎng)格法對(duì)齒輪的線接觸彈流潤(rùn)滑問(wèn)題進(jìn)行求解，探討不同參數(shù)對(duì)計(jì)算結(jié)果的影響，發(fā)現(xiàn)油膜壓力及油膜厚度對(duì)于載荷和彈性模量的變化較為敏感；在對(duì)齒輪等溫彈流潤(rùn)滑問(wèn)題求解的基礎(chǔ)上，考慮熱效應(yīng)、時(shí)變效應(yīng)及壁面粗糙度對(duì)潤(rùn)滑效果的影響。

關(guān)鍵詞：彈流潤(rùn)滑；直齒圓柱齒輪；熱效應(yīng)；粗糙度；時(shí)變效應(yīng)

中圖分類號(hào)：TH132.41" " " 文獻(xiàn)標(biāo)志碼：A" " " " " 文章編號(hào)：2095-2945（2023）13-0117-05

Abstract： In this paper， the problem of linear contact elastohydrodynamic lubrication of gears is solved based on multigrid method， and the influence of different parameters on the calculation results is discussed. It is found that oil film pressure and oil film thickness are sensitive to the change of load elastic modulus. On the basis of solving the problem of isothermal elastohydrodynamic lubrication of gears， the effects of thermal effect， time-varying effect and wall roughness on lubrication are considered.

Keywords： elastohydrodynamic lubrication; spur gear; thermal effect; roughness; time-varying effect

齒輪傳動(dòng)依靠輪齒之間的嚙合來(lái)傳遞速度和力矩，是傳動(dòng)系統(tǒng)中主要的承力部件。作為主要承力部件，齒輪失效的絕大多數(shù)情況都是潤(rùn)滑問(wèn)題導(dǎo)致。通過(guò)齒輪嚙合區(qū)域的潤(rùn)滑油油膜厚度及油膜壓力的分布情況，判斷齒輪的潤(rùn)滑情況進(jìn)而有效防止齒輪發(fā)生失效。由于齒輪潤(rùn)滑試驗(yàn)難度大，且精度難以保證，尤其是難以獲取齒輪嚙合區(qū)的相關(guān)數(shù)據(jù)，因此通過(guò)彈流的數(shù)值計(jì)算來(lái)得到齒輪嚙合區(qū)域的潤(rùn)滑情況是十分有必要的。

對(duì)于彈流潤(rùn)滑問(wèn)題主要是對(duì)基本方程進(jìn)行求解，國(guó)內(nèi)外研究學(xué)者主要圍繞著求解方法展開工作，主要的求解方法有直接迭代法、牛頓法、逆解法及多重網(wǎng)格法。

直接迭代法是根據(jù)初始?jí)毫?chǎng)計(jì)算得到潤(rùn)滑油油膜厚度、粘度、密度，然后求解雷諾方程得到新的壓力場(chǎng)分布情況，如此反復(fù)迭代直至得到滿足收斂精度的潤(rùn)滑油油膜壓力分布與油膜厚度分布。Petrusevich[1]于1951年首次借助直接迭代法得到了彈流潤(rùn)滑問(wèn)題的完整數(shù)值解。但是該方法收斂速度較慢，同時(shí)精度存在一定問(wèn)題。1962年Stephenson等[2]利用低松弛迭代算法對(duì)直接迭代法進(jìn)行了改進(jìn)，大大提高了直接迭代法的收斂速度。Hamrock等[3]采用直接迭代法對(duì)點(diǎn)接觸彈流潤(rùn)滑問(wèn)題進(jìn)行求解，并利用回歸得到了膜厚計(jì)算公式，為后續(xù)工程設(shè)計(jì)提供了參考。

牛頓法求解前需要對(duì)方程進(jìn)行離散，根據(jù)離散方法的不同，可以分為牛頓有限元法及牛頓有限差分法。1977年，小野京右[4]利用牛頓法對(duì)線接觸彈流問(wèn)題進(jìn)行求解，獲得了輕、中載荷下的數(shù)值解。Houpert等[5]在此基礎(chǔ)上對(duì)該算法進(jìn)行了改進(jìn)，大大提升了計(jì)算速度，同時(shí)使其適用于重載荷條件。OH[6]利用牛頓法對(duì)非穩(wěn)態(tài)變載荷工況下的點(diǎn)接觸彈流問(wèn)題進(jìn)行求解。Hughes等[7]對(duì)牛頓法進(jìn)行了改進(jìn)，大大提高了其在求解彈流問(wèn)題時(shí)的收斂速度與收斂的穩(wěn)定性。

逆解法首先根據(jù)初始?jí)毫Ψ植记蠼饽ず?，根?jù)膜厚對(duì)Reynold方程進(jìn)行求解獲得新的膜厚分布情況，再對(duì)壓力分布進(jìn)行修正，如此反復(fù)直至得到滿足收斂精度的潤(rùn)滑油油膜壓力分布與油膜厚度分布。Dowson等[8]利用逆解法對(duì)彈流潤(rùn)滑問(wèn)題進(jìn)行求解，結(jié)果發(fā)現(xiàn)該方法計(jì)算重載工況十分合適。同時(shí)基于計(jì)算結(jié)果通過(guò)推導(dǎo)得到了線接觸彈流問(wèn)題的油膜最小厚度的計(jì)算公式。Cheng等[9]使用逆解法對(duì)線接觸熱彈流潤(rùn)滑問(wèn)題進(jìn)行求解，結(jié)果發(fā)現(xiàn)該方法收斂速度較快且很穩(wěn)定。Hou等[10]對(duì)逆解法進(jìn)行了一定程度的改進(jìn)，使得其收斂速度得到了一定程度的提高。之后又對(duì)考慮熱效應(yīng)的點(diǎn)接觸熱彈流潤(rùn)滑問(wèn)題進(jìn)行了求解。

多重網(wǎng)格法是將同一計(jì)算域劃分為多套疏密程度不同的網(wǎng)格，然后將待求解的方程在每一層網(wǎng)格上均按照相同格式進(jìn)行離散，輪流在多套網(wǎng)格上進(jìn)行迭代求解，最終得到滿足進(jìn)度要求的結(jié)果。1986年，Lubrecht等[11]第一次利用多重網(wǎng)格法對(duì)點(diǎn)（線）接觸等溫彈流潤(rùn)滑問(wèn)題進(jìn)行求解，發(fā)現(xiàn)多重網(wǎng)格法無(wú)論是在收斂速度還是計(jì)算精度上都有較好的表現(xiàn)。后續(xù)，Lubrecht針對(duì)膜厚方程中的彈性變形項(xiàng)提出了多重網(wǎng)格積分法，該方法大大提高了計(jì)算效率。Venner[12]同時(shí)使用多重網(wǎng)格積分法與多重網(wǎng)格法對(duì)彈流潤(rùn)滑問(wèn)題進(jìn)行求解，計(jì)算速度大大加快的同時(shí)，在重載荷工況下表現(xiàn)也較為良好。Wang等[13]改進(jìn)了Venner的求解算法，大大提高了計(jì)算速度，并且對(duì)極重載荷下的彈流問(wèn)題進(jìn)行求解獲得了最大赫茲接觸壓力為9.2 GPa工況下的數(shù)值解。

多重網(wǎng)格法求解彈流問(wèn)題收斂快，而且較為穩(wěn)定，在極重載工況下也有較好的表現(xiàn)，被認(rèn)為是彈流計(jì)算中最為高效穩(wěn)定的方法，本文利用的方法就是基于多重網(wǎng)格法提出的。

1" 等溫線接觸彈流潤(rùn)滑分析

對(duì)于彈性流體動(dòng)力潤(rùn)滑問(wèn)題基本方程主要包括Reynolds方程、膜厚方程、載荷平衡方程、黏度-壓力方程和密度-壓力方程。

1.1" 基本方程

1.1.1" Reynolds方程

Reynolds方程用到了下列7個(gè)假設(shè)：①不考慮體積力對(duì)流體的作用；②固體表面處的流體速度為零；③認(rèn)為流體壓力在膜厚方向上不發(fā)生變化；④接觸表面的曲率半徑遠(yuǎn)大于膜厚；⑤流體為Newton流體；⑥流動(dòng)狀態(tài)為層流；⑦流體慣性力與黏性力相比可以忽略不計(jì)。

根據(jù)上述，假設(shè)則對(duì)于穩(wěn)態(tài)等溫線接觸條件下牛頓流體的雷諾方程為

（1）

式中：？籽為潤(rùn)滑油密度（kg/m3）；h為潤(rùn)滑油油膜厚度（m）；？滋為潤(rùn)滑油動(dòng)力黏度（Pa·s）；us為表面卷吸速度（m/s）。

邊界條件為

式中：xin，xout分別為計(jì)算域入口、出口。

1.1.2" 膜厚方程

對(duì)于不考慮壁面粗糙度的線接觸潤(rùn)滑問(wèn)題油膜厚度方程主要包括3項(xiàng)：中心膜厚、幾何膜厚、彈性變形

式中：δ（x）為接觸表面的彈性變形（m）。

式中：c為待定常數(shù)；E為彈性模量。

1.1.3" 載荷平衡方程

1.1.4" 黏度-壓力方程

式中：？滋0為初始工況下潤(rùn)滑油黏度（Pa·s）。

1.1.5" 密度-壓力方程

式中：？籽0為初始工況下潤(rùn)滑油密度（kg/m3）。

對(duì)上述方程進(jìn)行無(wú)量綱化及離散化，然后再基于多重網(wǎng)格法求解。

1.2" 線接觸彈流潤(rùn)滑特性分析

計(jì)算所采用的各項(xiàng)參數(shù)見(jiàn)表1。

通過(guò)多重網(wǎng)格法得到的線接觸彈流潤(rùn)滑問(wèn)題的壓力分布和油膜形狀如圖1所示。

由圖1可知，出口處間隙的形狀急劇擴(kuò)大導(dǎo)致出口區(qū)域的油膜壓力急劇下降進(jìn)而導(dǎo)致固體表面所產(chǎn)生的彈性形變量變小，引起了油膜厚度方向上的收縮，而油膜的收縮使得間隙中具有較大黏度的油液流動(dòng)受阻，所以產(chǎn)生了第二壓力峰。

不同載荷下潤(rùn)滑油油膜壓力分布情況與油膜厚度分布情況如圖2、圖3所示。隨著載荷的增加，第二壓力峰出現(xiàn)的位置不斷靠近出口區(qū)域；同時(shí)第二壓力峰的峰值變化量隨著載荷增加不斷減少；潤(rùn)滑油油膜厚度也逐漸減小。最小油膜厚度出現(xiàn)的位置不斷靠近出口區(qū)域，同時(shí)區(qū)域中的最小油膜厚度也逐漸減小。

不同卷吸速度下的潤(rùn)滑油油膜壓力分布情況與油膜厚度分布情況如圖4、圖5所示。卷吸速度對(duì)潤(rùn)滑油壓力分布情況的影響表現(xiàn)在：隨著卷吸速度的增大，第二壓力峰出現(xiàn)的位置不斷遠(yuǎn)離出口區(qū)域；同時(shí)第二壓力峰的峰值隨著卷吸速度的增大而不斷增大。卷吸速度對(duì)潤(rùn)滑油油膜厚度分布的影響表現(xiàn)在：隨著卷吸速度的增大，潤(rùn)滑油油膜厚度也逐漸增大，最小油膜厚度出現(xiàn)的位置不斷遠(yuǎn)離出口區(qū)域，同時(shí)區(qū)域中的最小油膜厚度也逐漸增大。

不同彈性模量線接觸彈流潤(rùn)滑問(wèn)題進(jìn)行計(jì)算，潤(rùn)滑油油膜壓力分布情況與油膜厚度分布情況如圖6、圖7所示?？梢园l(fā)現(xiàn)彈性模量對(duì)潤(rùn)滑油壓力分布情況的影響表現(xiàn)在：隨著彈性模量的增加，第二壓力峰出現(xiàn)的位置不斷遠(yuǎn)離出口區(qū)域；第二壓力峰的峰值隨著彈性模量的增加不斷增大。彈性模量對(duì)潤(rùn)滑油油膜厚度分布的影響表現(xiàn)在：隨著彈性模量的增加，接觸區(qū)域潤(rùn)滑油油膜厚度逐漸減小，但變化不明顯，最小油膜厚度出現(xiàn)的位置不斷靠近出口區(qū)域，同時(shí)區(qū)域中的最小油膜厚度也逐漸減小。

2" 考慮粗糙度和時(shí)變效應(yīng)熱彈流潤(rùn)滑

在齒輪嚙合過(guò)程中，由于齒輪輪齒表面發(fā)生相對(duì)滑動(dòng)引起大量的耗散熱，進(jìn)而引起潤(rùn)滑油及齒輪輪齒溫度的升高。潤(rùn)滑油的黏度會(huì)隨著溫度的升高而降低，導(dǎo)致油膜承載力大幅度下降。因此齒輪的彈流計(jì)算中有必要考慮熱效應(yīng)的影響。

2.1" 基本方程

在探討熱效應(yīng)對(duì)齒輪潤(rùn)滑情況的影響時(shí)，還需要引入能量方程，另外其余參數(shù)均要考慮溫度的影響。

2.1.1" 能量方程

。 （8）

式（8）中各項(xiàng)根據(jù)其物理意義從左到右依次代表了對(duì)流換熱項(xiàng)、導(dǎo)熱項(xiàng)、擠壓產(chǎn)熱項(xiàng)和粘性耗散熱項(xiàng)。

2.1.2" 黏度-溫度-壓力方程

式中：T0為環(huán)境溫度（K）。

2.1.3" 密度-溫度-壓力方程

。（10）

對(duì)式（8）、式（9）、式（10）進(jìn)行無(wú)量綱化、離散化，參照等溫彈流處理方法，在此不再贅述。采用逐列掃描法對(duì)潤(rùn)滑油的溫度分布進(jìn)行求解。

2.2" 熱彈流潤(rùn)滑特性分析

對(duì)齒輪熱彈流潤(rùn)滑進(jìn)行計(jì)算，可以得到齒輪節(jié)點(diǎn)處溫度的分布情況如圖8所示。

從圖8可以發(fā)現(xiàn)潤(rùn)滑油油膜溫度隨著油膜壓力的增大而快速升高。之后隨著油膜壓力的減小油膜溫度也逐漸降低，另外油膜溫度在膜厚方向上沿著上下表面中心近似于對(duì)稱分布，中心油膜溫度變化最為劇烈，越接近固體表面溫度越低，同時(shí)溫升也趨于平緩。

2.3" 考慮壁面粗糙度和時(shí)變效應(yīng)的熱彈流潤(rùn)滑

齒輪嚙合過(guò)程是一個(gè)時(shí)變過(guò)程，Reynolds方程中多了表示時(shí)變效應(yīng)的項(xiàng)。粗糙度對(duì)于熱彈流潤(rùn)滑問(wèn)題的影響主要是通過(guò)膜厚方程中引入表示表面粗糙的項(xiàng)。在考慮時(shí)變效應(yīng)的情況下，存在壁面粗糙度與光滑壁面的油膜特性如圖9所示。

從圖9看出，當(dāng)考慮壁面粗糙度時(shí)油膜壓力及油膜厚度分布情況會(huì)產(chǎn)生一定的波動(dòng)，但是整體上與光滑壁面的分布趨勢(shì)相同，彈流潤(rùn)滑的典型特征——第二壓力峰及頸縮現(xiàn)象均有出現(xiàn)。最小油膜厚度出現(xiàn)的位置基本相同，光滑壁面情況下的最小油膜厚度略大于粗糙壁面情況。

3" 結(jié)論

本文建立了直齒圓柱齒輪彈流潤(rùn)滑問(wèn)題計(jì)算模型，分析了熱效應(yīng)、壁面粗糙度及時(shí)變效應(yīng)對(duì)齒輪潤(rùn)滑效果的影響，獲得以下結(jié)論。

1）使用多重網(wǎng)格法對(duì)不同載荷下的線接觸彈流潤(rùn)滑問(wèn)題進(jìn)行求解，分析了不同參數(shù)對(duì)潤(rùn)滑效果的影響，油膜壓力及油膜厚度對(duì)于載荷的變化較為敏感。

2）熱彈流得到節(jié)點(diǎn)處潤(rùn)滑油溫升主要集中在靠近接觸中心區(qū)域處，而潤(rùn)滑油油膜壓力及油膜厚度是由入口處的卷吸作用所決定的，所以熱效應(yīng)對(duì)于潤(rùn)滑油油膜壓力及油膜厚度分布情況的影響并不明顯。

3）對(duì)考慮壁面粗糙度及時(shí)變效應(yīng)的齒輪熱彈流潤(rùn)滑問(wèn)題進(jìn)行求解，結(jié)果表明最大油膜壓力及最小油膜厚度沿著嚙合線的分布趨勢(shì)與光滑壁面情況下基本相同，壁面粗糙度的存在導(dǎo)致其在相應(yīng)位置出現(xiàn)了一定程度的波動(dòng)。

參考文獻(xiàn)：

[1] PETRUSEVICH A I. Fundamental conclusions from the contact-hydrodynamic theory of lubrication[J]. Izv. Akad. Nauk. SSSR （OTN）， 1951（2）：209.

[2] STEPHENSON R R， OSTERLE J F. A direct solution of the elasto-hydrodynamic lubrication problem[J]. ASLE TRANSACTIONS， 1962，5（2）：365-374.

[3] HAMROCK B J， DOWSON D. Elastohydrodynamic lubrication of elliptical contacts for materials of low elastic modulus I-Fully flooded conjunction[J]. ASME Journal of Lubrication Technology， 1978（100）：236-245.

[4] 小野京右.彈流的數(shù)值解析法[J].日本潤(rùn)滑協(xié)會(huì)第二十二回東京講習(xí)會(huì)所，1977，27（10）：726-732.

[5] HOUPERT L G， HAMROCK B J. Fast approach for calculating film thicknesses and pressures in elastohydrodynamically lubricated contacts at high loads[J]. Trans. ASME J. Tribol.， 1986（108）：411.

[6] OH K P. The numerical solution of dynamically loaded elastohydrodynamic contact as a nonlinear complementary problem[J]. Journal of lubrication technology， 1984，106（1）：88-95.

[7] HUGHES T G， ELCOATE C D， EVANS H P. A novel method for integrating first-and second-order differential equations in elastohydrodynamic lubrication for the solution of smooth isothermal， line contact problems[J]. International Journal for Numerical Methods in Engineering， 1999，44（8）：1099-1113.

[8] DOWSON D， HIGGINSON G R. A numerical solution to the elasto-hydrodynamic problem[J]. Journal of mechanical engineering science， 1959，1（1）：6-15.

[9] CHENG H S， STERNLICHT B. A numerical solution for the pressure， temperature and film thickness between two infinitely long， lubricated rolling and sliding cylinders， under heavy loads[J]. Trans. ASME.， J. Bas. Eng.， 1964（87）：695.

[10] HOU K P， ZHU D， WEN S Z. An inverse solution to the point contact EHL problem under heavy loads[J]. Journal of tribology， 1987，109（3）：432-436.

[11] LUBRECHT A A， TEN NAPEL W E， BOSMA R. Multigrid， an alternative method of solution for two-dimensional elastohydrodynamically lubricated point contact calculations[J]. Journal of tribology， 1987，109（3）：437-443.

[12] VENNER C H. Multilevel solutions of the EHL line and point contact problems [D]. Ph D thesis， Enschede： University of Twente， 1991.

[13] WANG J， QU S， YANG P. Simplified multigrid technique for the numerical solution to the steady-state and transient EHL line contacts and the arbitrary entrainment EHL point contacts[J]. Tribology International， 2001，34（3）：191-202.