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

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

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

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

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

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

不考慮溫度對(duì)潤(rùn)滑油黏度的影響，量綱一化的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 軸承潤(rùn)滑示意圖及單軸瓦展開(kāi)圖

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

(7)

(8)

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

[K]{a}=

(9)

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

圖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ù)項(xiàng)即為油膜剛度矩陣

(15)

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

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

(16)

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

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

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

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

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

Δf=frN-fr0

(17)

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

(18)

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

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

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

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

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

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

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

4 算例分析

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

圖5 軸承的供油和潤(rùn)滑方式

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

首先計(jì)算軸頸中心做如下軌跡的橢圓運(yùn)動(dòng)時(shí)的油膜力、油膜剛度、阻尼。

τ∈[0,2π]

(19)

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

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

圖7 量綱一油膜力

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

圖8 量綱一油膜剛度

圖9 量綱一油膜阻尼

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

(20)

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

(21)

令軸心速度為0，如圖10所示，為從式(21)所示的8個(gè)不同的初始位置出發(fā)，按照載荷逼近法，軸心靜平衡位置的Newton迭代過(guò)程?？芍瑥牟煌恢贸霭l(fā)，迭代最終均收斂于同一點(diǎn)(0.137 1,-0.057 5)，從而驗(yàn)證了文中所提出方法的有效性。

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

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

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

5 結(jié)論

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

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

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

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

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