萬 召1 王衛(wèi)國1 龍新華 孟 光

(1.中國航發(fā)商用航空發(fā)動機(jī)有限責(zé)任公司 上海 200241；2.上海交通大學(xué)機(jī)械與動力工程學(xué)院 上海 200240)

滑動軸承廣泛應(yīng)用于大型、重載轉(zhuǎn)子，如燃?xì)?蒸汽輪機(jī)轉(zhuǎn)子、水輪發(fā)電機(jī)組轉(zhuǎn)子等，滑動軸承油膜力及油膜剛度、阻尼特性對轉(zhuǎn)子-軸承系統(tǒng)的臨界轉(zhuǎn)速、動力響應(yīng)和穩(wěn)定性等都有著顯著的影響，求解滑動軸承的非線性油膜力和油膜特性系數(shù)是轉(zhuǎn)子-軸承動力學(xué)分析的核心內(nèi)容之一。

滑動軸承非線性油膜力的求解方法主要有2種：解析法和數(shù)值算法[1]。解析法是基于對油膜壓力分布的簡單假設(shè)而得到的，通常有長軸承模型、短軸承模型。解析法求解方便、計算速度快，但精度低。數(shù)值算法主要包括有限差分法[2]和有限元法[3]，二維有限差分和二維有限元法計算精度高，但效率低[4]。張偉忠等[5]對比分析了有限差分與長軸承、短軸承等模型,指出各種模型的誤差和適用范圍；楊金福等[6]分析了幾種解析模型的區(qū)別和聯(lián)系；ZHENG等[7-8]提出了Reynolds方程的一維有限元算法，基于該算法研究了求解油膜剛度、阻尼系數(shù)及軸心靜平衡位置的Newton迭代法，指出一維算法可以在求解油膜力的同時獲得油膜剛度、阻尼系數(shù)，而二維有限差分法則需要求解多次油膜力之后再做差分才能獲得取油膜動力特性系數(shù)；肖忠會[9]詳細(xì)研究了一維有限元算法的變分理論、算法實現(xiàn)等，對比了有限差分與一維有限元算法的計算效率、精度，表明兩者誤差小于1%，而一維算法的效率是有限差分法的數(shù)百倍；黃文虎等[10]在綜合對比分析長軸承、短軸承、有限差分、一維有限元等模型與算法之后，指出一維有限元算法兼具精度與速度優(yōu)勢。盡管一維算法有上述優(yōu)點，但該算法的變分理論基礎(chǔ)復(fù)雜、數(shù)值實現(xiàn)難度高，應(yīng)用有限；另外，Newton迭代法對初值十分敏感，求靜平衡位置時容易發(fā)散；此外，受傳統(tǒng)二維有限差分法計算效率限制，對外載荷連續(xù)變化時軸心靜平衡位置的變化趨勢的研究鮮有報道。

本文作者研究了Reynolds方程的一維有限元算法的降維方法，以及角向油膜壓力分布的一維離散方法。針對初值選擇不當(dāng)時，Newton迭代易不收斂的問題，提出了一種載荷逼近方法，能夠有效求解任意外載荷下的軸心靜平衡位置?；谠摲椒?，研究了某滑動軸承在軸頸中心做橢圓進(jìn)動時的油膜壓力分布、油膜剛度、阻尼特性，同時研究了在外載荷連續(xù)變化時，軸心靜平衡位置的變化趨勢。

1 Reynolds方程的一維有限元快速解法

不考慮溫度對潤滑油黏度的影響，量綱一化的Reynolds方程可寫作

(1)

油膜壓力分布寫成如下形式

p(θ,ζ)=g(ζ)·r(θ)

(2)

其中g(shù)(ζ)為軸向壓力分布函數(shù)，r(θ)為角向壓力分布函數(shù)，即角向壓力分布與軸向壓力分布在形式上是解耦的。進(jìn)一步假設(shè)軸向壓力分布滿足如下雙曲函數(shù)形式：

g(ζ)=cosh(kλ)-cosh(kζ)

(3)

其中k是與角向壓力分布r(θ)相關(guān)的迭代參數(shù)，通常迭代初值可取k=1。油膜力的積分形式為

(4)

進(jìn)一步可得到

(5)

其中參數(shù)c3滿足如下公式

(6)

圖1 軸承潤滑示意圖及單軸瓦展開圖

由此可知，只要求得角向油膜壓力分布函數(shù)r(θ)，即可求得油膜力。如圖2所示，將軸向油膜壓力分布作一維有限元離散，即角向壓力分布滿足

(7)

(8)

其中ai為角向坐標(biāo)為θi時的油膜壓力，Li為單元型函數(shù)，可采用線性Lagrange插值型函數(shù)，則油膜壓力分布可轉(zhuǎn)化為求解如下線性方程組

[K]{a}=

(9)

其中[K]為三對角帶狀矩陣，對其做LU分解，進(jìn)而通過線性迭代即可求得油膜壓力分布向量{a}，再通過積分即可求得油膜力。

圖2 角向油膜的一維有限元離散

2 求解軸心靜平衡位置的Newton迭代法

{fb(xs,ys,0,0)}+{fr}=0

(10)

(11)

圖3 軸頸受力與軸心靜平衡位置

f(u)={fb(x,y)}+{fr}=0

(12)

由Newton迭代法可知，迭代求解公式為

(14)

式(14)中的導(dǎo)數(shù)項即為油膜剛度矩陣

(15)

則迭代公式可進(jìn)一步寫成

uk+1=uk-(f′(uk))-1(fb(uk)+fr)

(16)

由以上可知，Newton迭代法求解外載荷下軸心靜平衡位置的步驟可歸納為

(1)給定迭代初值u0={x0,y0}，求解此時的油膜剛度矩陣f′(u0)和油膜力fb(u0)，代入迭代公式(16)，求解下一個軸心位置u1={x1,y1}；

3 Newton迭代的收斂性與求解軸心靜平衡位置的載荷逼近法

由計算方法可知，Newton迭代法是一種局部收斂的迭代方法，用Newton迭代法求解方程的根時，初值的選取對于迭代能否收斂具有決定性的影響，當(dāng)?shù)植扛姆较蜻M(jìn)行時，收斂速度很快(迭代是2階收斂的)；當(dāng)?shù)畴x局部解的時候，迭代發(fā)散，會陷入死循環(huán)。當(dāng)需要大量求解不同載荷下軸頸的靜平衡位置時，如何快速確定迭代初值，并保證迭代收斂，是一個必須考慮的問題。

為此，文中提出一種可以從任意軸心位置出發(fā)，尋找到任意載荷下軸心靜平衡位置的載荷逼近方法。如圖4所示，設(shè)要求解外載荷frN下軸心的靜平衡位置，任意選取某一軸心位置us0作為初始位置，計算軸心為us0時的油膜力fb0=fb(us0,0)，令初始外載荷fr0=-fb0，則初始載荷fr0到外載荷frN之間的矢量差Δf為

Δf=frN-fr0

(17)

為了從初始外載荷逼近目標(biāo)外載荷，可將矢量差等分為N份，則矢量步長為

(18)

這樣可以由初始載荷fr0通過N步來逼近目標(biāo)載荷frN，進(jìn)而求得任意外載荷frN下軸心的靜平衡位置。載荷逼近法的迭代流程總結(jié)如下：

(2)計算由初始載荷fr0到目標(biāo)載荷frN的距離矢量Δf=frN-fr0，載荷fbN可由初始載荷fr0分N步逼近，逼近矢量步長h=Δf/N；

(3)采用Newton迭代法，并取初值為us0，由Newton迭代公式(16)，求外載荷為fr1=fr0+h時的靜平衡位置us1；

(4)重復(fù)上述迭代過程，由外載荷fri=fr0+i×h時的靜平衡位置usi，按照公式(16)求載荷為fi+1=f0+(i+1)×h時的靜平衡位置us,i+1;

(5)循環(huán)迭代，最終即可求得載荷為frN=fr0+N×h時軸心的靜平衡位置usN，并可同步求得該靜平衡位置下的油膜剛度、阻尼系數(shù)。

可以看到，上述載荷逼近法求解靜平衡位置，事實上也是一種初值選取方法，只是其所選取的Newton迭代初值可以是任意軸心位置?？梢灶A(yù)見，只要選取適當(dāng)?shù)妮d荷逼近步長，總是可以從某一載荷下的靜平衡位置，逼近、并收斂到另外一個載荷下的靜平衡位置。根據(jù)載荷逼近法，理論上可以計算軸頸在任意載荷下的靜平衡位置，這同時意味著可以將任意軸心位置作為迭代起始位置。載荷逼近法也利用了油膜力的一維有限元算法的速度優(yōu)勢。

圖4 求解任意外載荷下軸心靜平衡位置的載荷逼近法

4 算例分析

如圖5所示為某試驗所采用的中心剖分式圓柱滑動軸承，軸瓦材料為耐磨的巴氏合金，軸承左右兩側(cè)各有一個進(jìn)油孔和儲油腔。軸承安裝時，剖分面處于水平位置，潤滑油經(jīng)過剖分面的進(jìn)油孔進(jìn)入儲油腔；具有一定偏心的軸頸在軸承內(nèi)轉(zhuǎn)動時，將儲油腔內(nèi)的潤滑油吸出，形成壓力油楔并建立動壓油膜實現(xiàn)潤滑；然后潤滑油通過軸承的前后端面泄出，并進(jìn)入回油油路。軸承相關(guān)參數(shù)如表1所示，其中滑油為32號透平油，動力黏度為0.027 6 Pa·s(40 ℃)，密度為850 kg/m3。

圖5 軸承的供油和潤滑方式

項目值項目值直徑40 mm寬度32 mm寬徑比0.8間隙比η0.2%進(jìn)油角(上)27°軸瓦包角126°

首先計算軸頸中心做如下軌跡的橢圓運動時的油膜力、油膜剛度、阻尼。

τ∈[0,2π]

(19)

圖6給出了軸心做上述橢圓運動時，上半瓦包角范圍內(nèi)的角向壓力分布，可以看到油膜壓力分布以及油膜破裂邊界隨軸心位置的變化情況。對油膜壓力進(jìn)行積分，可以得到上、下半瓦油膜力合力的變化曲線，如圖7所示，分析可知軸心做橢圓運動時，兩個方向上的油膜力不相等，是非對稱的；而兩個方向的油膜力fu、fv均滿足f(τ)=-f(τ+π)， 即沿軸承中心對稱的兩個位置處的油膜力大小相等、方向相反；軸心渦動一個周期，油膜力也波動一個周期。

圖6 角向油膜壓力分布(上半瓦)

圖7 量綱一油膜力

圖8、圖9進(jìn)一步給出了軸心作橢圓運動時的油膜剛度、阻尼變化曲線，可以看到剛度、阻尼系數(shù)在兩個正交方向上也是非對稱的，且是相互耦合的、均存在交叉項。從絕對值上看，交叉剛度kvu最大，且交叉剛度kuv≠kvu，直接剛度有kuu>kvv；油膜阻尼也存在交叉項，但兩個方向上的交叉阻尼是相等的cuv=cvu，這是油膜剛度與油膜阻尼之間的重要區(qū)別之一，交叉阻尼會使得油膜發(fā)生失穩(wěn)。分析各個特性系數(shù)發(fā)現(xiàn)，有k(τ)=k(τ+π)，c(τ)=c(τ+π)，即在沿軸承中心對稱的兩個位置處，各系數(shù)的大小、符號均分別相同，表明在軸心渦動一個周期時，剛度、阻尼系數(shù)變化兩個周期。

圖8 量綱一油膜剛度

圖9 量綱一油膜阻尼

進(jìn)一步，按照文中給出的載荷逼近與Newton迭代方法，計算軸頸受到如式(20)所示的隨時間變化的量綱一外載荷作用下，軸心靜平衡位置的變化情況。

(20)

初始載荷為fr(τ=0)={0,-0.05}T，方向向下，首先需要求解在初始載荷下軸心的靜平衡位置。為驗證文中所提出的載荷逼近方法，選取如下8個位置作為初始迭代的軸心位置

(21)

令軸心速度為0，如圖10所示，為從式(21)所示的8個不同的初始位置出發(fā)，按照載荷逼近法，軸心靜平衡位置的Newton迭代過程。可知，從不同位置出發(fā)，迭代最終均收斂于同一點(0.137 1,-0.057 5)，從而驗證了文中所提出方法的有效性。

圖10 不同起始位置時軸心靜平衡位置的迭代

進(jìn)一步，以初始載荷下的軸心靜平衡位置為起點，采用Newton迭代法，求解在式(21)所示載荷下，軸心靜平衡位置的變化趨勢，結(jié)果如圖11所示。可知：在軸頸轉(zhuǎn)速、旋轉(zhuǎn)方向(逆時針)不變，載荷方向向下也保持不變，當(dāng)軸頸所受外載荷較小時，軸頸中心靠近軸承中心，偏心率很小，油膜剛度、阻尼均較??；當(dāng)外載荷逐漸增大時，軸心靜平衡位置逐步向右(與軸頸旋轉(zhuǎn)方向有關(guān))下沉，遠(yuǎn)離外載荷方向；當(dāng)載荷超過一定值后，軸心位置又會向左下方運動，逐步靠近外載荷的方向，最后靠近正下方的位置，軸承偏心率逐步增大，剛度、阻尼系數(shù)均增大，這與工程中實際情況是一致的。

圖11 外載荷下軸心靜平衡位置變化(量綱為一)

5 結(jié)論

(1)提出一種求解任意外載荷下軸心靜平衡位置的載荷逼近法，對某中心剖分式圓柱滑動軸承軸心作橢圓運動時的油膜力、油膜剛度、阻尼進(jìn)行計算，結(jié)果表明：

(a)油膜力、油膜剛度、阻尼在兩個正交方向上是非對稱的。且油膜剛度、阻尼存在交叉項，油膜阻尼的交叉項cuv=cvu相等，油膜剛度的交叉項kuv≠kvu不等。

(b)軸心以軸承中心為圓心做周期性橢圓運動時，沿軸承中心對稱的兩點處的油膜剛度、阻尼的各個分量的大小、方向均相等，即在一個軸心進(jìn)動周期內(nèi)，剛度、阻尼各分量變化2個周期；油膜力各分量在中心對稱兩點處幅值分別相等、符號相反，即在軸心一個進(jìn)動周期內(nèi)，油膜力也變化1個周期。

(c)當(dāng)軸心所受外載荷方向不變、幅值持續(xù)增大時，軸心偏心率增大，軸心靜平衡位置的先向遠(yuǎn)離載荷方向運動，后向靠近載荷方向運動，最終軸心會靠近載荷方向并偏向某一側(cè)，但不會與載荷方向重疊，偏離方向受軸頸自轉(zhuǎn)方向的影響。

(2)文中所提出的載荷逼近法，能夠有效地從任意初始位置迭代求解得到在任意外載荷下軸心的靜平衡位置，較好地解決Newton迭代法求軸心靜平衡位置時迭代初值的選取問題。