點接觸彈流問題快速直接算法研究*

夏伯乾 方保江 郭 紅 張紹林

(鄭州大學(xué)機械與動力工程學(xué)院 河南鄭州 450001)

相比線接觸彈流問題，點接觸彈流計算因為要在二維求解區(qū)域上求解Reynolds方程和彈性變形方程，問題更復(fù)雜，計算量也更大，因而研究和發(fā)展高效、使用參數(shù)范圍廣的點接觸EHL問題求解方法具有重要意義。

求解點接觸EHL問題，主要有逆解法[1-3]、直接迭代法[4-6]、Newton-Raphson法[7-9]和多重網(wǎng)格法[10-12]，這幾種方法各有優(yōu)點[13]，但都需要在二維求解域上進行求解，方程組維數(shù)很大，計算效率不高。

2010年，文獻[6]在復(fù)合直接迭代法的基礎(chǔ)上，引入求解大型方程組的分塊求解思想，將點接觸彈流的二維求解問題化為一系列一維問題求解，提出了按行分塊、逐行計算的計算方法，極大地降低了方程組的維數(shù)，使計算效率比復(fù)合直接迭代法提高至少10倍以上。

為進一步提高點接觸EHL問題的求解效率，本文作者在文獻[6]的基礎(chǔ)上，通過分析點接觸EHL求解時網(wǎng)格劃分的特點，提出了按列分塊、逐列求解的求解點接觸EHL問題的計算方法，使計算效率相比文獻[6]的按行分塊、逐行計算的求解方法進一步顯著提高。

1 基本方程

1.1 Reynolds方程

(1)

邊界條件：

1.2 油膜幾何方程

(2)

式中：Rx、Ry分別為接觸體在x、y方向上的曲率半徑；h0為待定常數(shù);δ(x,y)表達式為

(3)

其中，Ω為有壓力分布的整個求解域；E′為綜合彈性模量。

1.3 潤滑黏壓關(guān)系

η=η0exp{(lnη0+9.67)[(1+5.1×10-9p)z-1]}

(4)

1.4 潤滑劑密壓關(guān)系

(5)

1.5 載荷平衡方程

(6)

2 求解區(qū)域的網(wǎng)格劃分

點接觸EHL求解區(qū)域網(wǎng)格劃分如圖1所示。

圖1 求解區(qū)域及網(wǎng)格劃分Fig 1 Solution area and meshing

3 方程的量綱一化

采用文獻[14]中的量綱一化參數(shù)，量綱一化后的Reynolds方程和油膜厚度方程分別為

(7)

(8)

4 點接觸彈流問題的高效直接算法研究

4.1 復(fù)合直接迭代法

在求解域的內(nèi)節(jié)點采用三點中心差分格式將方程(7)進行離散，H項用式(8)代入，經(jīng)整理得到如下的差分方程組：

(9)

式(9)中各參數(shù)的含義詳見文獻[14]。

方程(9)可寫成如下格式：

(10)

求解方程組(10)的計算流程如圖2所示。

圖2 計算流程Fig 2 Calculation process

4.2 按行求解的點接觸EHL高效直接算法

文獻[6]將方程組(10)改寫成如下分塊矩陣形式：

(11)

其中：

并針對式(11)構(gòu)造了如下迭代格式

(12)

4.3 點接觸EHL求解域的離散特點

求解點接觸彈流潤滑問題時，需要將求解域離散，如圖1所示。

一般地，X方向(卷吸速度u的方向)需要劃分的節(jié)點數(shù)較多，而Y方向需要劃分的節(jié)點數(shù)較少，X方向與Y方向的節(jié)點數(shù)量相差很大。Y向一般只需劃分20～40個節(jié)點，最多50個節(jié)點就已足夠，X方向的節(jié)點數(shù)則至少要是Y方向節(jié)點數(shù)的2～3倍，甚至更多，而且載荷越大，為保證計算精度，X方向需要劃分的節(jié)點數(shù)越多。

4.4 不同迭代格式的計算量比較

復(fù)合直接迭代法同時在二維計算域上通過方程組(10)求解壓力，方程組(10)是一個m×n維的大型方程組，用Gauss法求解時的計算量約為(m×n)3/3[15]。文獻[6]按行求解的直接迭代法因為把方程組(10)化為m個n維的方程組逐個求解，使用Gauss法求解時的計算量約為mn3/3，僅是復(fù)合直接迭代法計算量的1/m2，所以按行求解的方法使計算效率大為提高。

注意到求解點接觸EHL時計算區(qū)域離散的特點：X方向節(jié)點數(shù)是Y方向節(jié)點數(shù)的至少2～3倍，甚至更多。如果將文獻[6]的按行分塊、逐行求解改為按列分塊、逐列求解，就可以把m×n維的方程組(10)化為n個更小型的m維的方程組，用Gauss法求解時的計算量將變成僅有nm3/3。相比按行求解，計算量僅為其m2/n2倍，按最小n=2m計算,計算量僅為逐行求解算法的1/4，因此可以進一步提高計算效率。

4.5 一種求解點接觸EHL的新方法——按列分塊、逐列求解

基于4.4節(jié)的分析，把方程組(10)重寫成以下分塊矩陣形式：

(13)

其中：

構(gòu)造如下新的迭代格式：

(14)

這樣方程組(10)就被降階成一組(n個)維數(shù)更小的(只有m維)的小型方程組，用直接法求解，速度更快，總的計算流程仍然如圖2所示。

5 新算法的正確性驗證

按照式(14)逐列求解格式，編寫了計算程序，并對多種工況條件進行了計算。

圖3—6示出了其中一種工況的計算結(jié)果。主要參數(shù)如下：量綱一載荷W=1.062×10-7,量綱一速度U=1.441×10-11,曲率半徑為Rx=Ry=25.4 mm，α=2 057.5 GPa，E′=146 GPa,η0=0.041 1 Pa·s。

圖3 三維壓力分布 圖4 三維膜厚分布 圖5 二維膜厚和壓力分布Fig 3 Three-dimensional pressure distribution Fig 4 Three-dimensional distribution of film thickness Fig 5 Two-dimensional pressure and film thickness distribution

圖6 膜厚等值線Fig 6 Contour of film thickness

從圖3—6中，可以看到點接觸彈流潤滑的4個典型特征：(1)油膜壓力二次峰現(xiàn)象；(2)接觸區(qū)域油膜厚度的平行區(qū)；(3)壓力二次峰處的頸縮現(xiàn)象；(4)膜厚等值線圖中的馬蹄形現(xiàn)象。這些特征驗證了文中計算方法的正確性。

6 新算法精確性與高效性驗證

為了檢驗新算法的計算精度與效率，對12種不同工況進行了求解，將計算結(jié)果與H-D公式進行了對比，并將計算耗時與文獻[6]的逐行計算的耗時進行了比較，結(jié)果見表1。這些算例的主要輸入?yún)?shù)是：α=2 057.5 GPa，E′=146 GPa,η0=0.041 1 Pa·s，Rx=Ry且取值范圍為15～25.4 mm。

表1 新算法計算精度與高效性驗證Table 1 The accuracy and efficiency of the new algorithm

從表1可以看到，文中算法計算得到的中心膜厚和最小膜厚與H-D公式的計算值相比較，相對誤差都在5%～10%，證明該算法具有令人滿意的計算精度；同時可以看到，在計算耗時上，在相同的計算條件下，文中算法比文獻[6]提出的逐行求解的快速直接算法在計算效率上又有了明顯提高，且網(wǎng)格越密，計算效率提高得越明顯。

7 算法的參數(shù)適用性檢驗

為檢驗文中算法的參數(shù)適用范圍，又對6種不同工況進行了計算，結(jié)果見表2。

表2 算法的參數(shù)適用范圍Table 2 The range of parameters applicable to the algorithm

過去一般認為，直接迭代法只適應(yīng)于輕、中載荷工況，在重載工況下難以獲得收斂解，但從表2可以看到，文中算法不僅在輕載工況、中載工況，而且在重載工況下都可以得到很好的收斂解，證明文中算法具有很寬的載荷參數(shù)適用范圍。

8 結(jié)論

在文獻[6]點接觸EHL快速直接算法的基礎(chǔ)上，通過分析點接觸EHL求解時求解區(qū)域離散的特點，提出求解點接觸EHL問題的按列分塊、逐列求解的新計算方法，并通過數(shù)值實例，證明所提出的新算法不僅具有更高的計算效率，而且具有更寬的參數(shù)適用范圍，是對文獻[6]所提出的點接觸EHL快速直接算法的進一步改進和發(fā)展。