毛文貴 韓 旭 劉桂萍

1.湖南大學(xué)，長沙，410082 2.湖南工程學(xué)院，湘潭，411101

0 引言

隨著機器向高速、大功率、大型化方向發(fā)展，滑動軸承越來越多地工作在高速重載的惡劣工況下，這對其性能設(shè)計提出了更加嚴格的要求。滑動軸承油膜動力特性系數(shù)直接影響著轉(zhuǎn)子—滑動軸承系統(tǒng)的臨界轉(zhuǎn)速、不平衡響應(yīng)、穩(wěn)定性等動力學(xué)行為。滑動軸承油膜動特性通?？捎?個系數(shù)來描述。早在1965年的時候，Lund［1］就采用4個剛度系數(shù)和4個阻尼系數(shù)建立了撓性轉(zhuǎn)子系統(tǒng)模型，揭示了油膜渦動失穩(wěn)的物理本質(zhì)。這些動力系數(shù)表達了油膜的動反力和軸頸的振動位移與速度之間的定量關(guān)系，精確計算滑動軸承的動力系數(shù)對轉(zhuǎn)子—滑動軸承系統(tǒng)動力學(xué)的研究和旋轉(zhuǎn)機械的設(shè)計起著非常重要的作用。

目前，國內(nèi)外學(xué)者大多通過有限差分法求解雷諾（Reynolds）方程來研究層流狀態(tài)下滑動軸承的動靜特性［2－3］。有限差分法通過動態(tài)雷諾方程對各項擾動參數(shù)求導(dǎo)，來得到各項擾動壓力的微分方程，通過求解這些擾動壓力的雷諾方程，得到各項擾動壓力分布，并積分得到滑動軸承的線性化剛度系數(shù)和阻尼系數(shù)。然而上述研究對滑動軸承動力特性系數(shù)計算的前提是，雷諾方程中的偏心率與偏位角要等于轉(zhuǎn)軸處于平衡位置時的偏心率與偏位角，而在實際求解過程中，轉(zhuǎn)子的載荷可以預(yù)先指定但轉(zhuǎn)軸平衡位置卻是未知的，且轉(zhuǎn)子從開始工作到穩(wěn)定在平衡位置這段時間內(nèi)，轉(zhuǎn)軸的偏心率和偏位角是時刻變化的。因此采用傳統(tǒng)靜態(tài)滑動軸承動力特性系數(shù)分析轉(zhuǎn)子穩(wěn)定在平衡位置之前這段時間內(nèi)的瞬態(tài)響應(yīng)是不準確的。李元生等［4］采用有限元法和有限差分法，通過在每個時刻用學(xué)科間迭代方法來實現(xiàn)轉(zhuǎn)子瞬態(tài)響應(yīng)與滑動軸承動力特性系數(shù)分析之間的耦合。有限差分法求解滑動軸承動力特性系數(shù)至少要4次調(diào)用雷諾方程，只有單元劃分足夠細，方可達到理想精度，故計算量很大，而小擾動速度的界限也不太好確定，因此不適于非線性分析。

高速重載滑動軸承中，流動油膜處于轉(zhuǎn)子動壁面和軸瓦靜壁面之間，流速呈線性分布，接近靜壁面區(qū)流速較小，為層流狀態(tài)，遠離靜壁面區(qū)（核心區(qū)）流速較大，為紊流狀態(tài)，油膜處于一種復(fù)合紊流狀態(tài)［5］。肖忠會等［6－7］針對層流狀態(tài)下非線性油膜力的求解提出了一種快速算法，該算法把二維變分不等方程降為一維，并將油膜力與動力特性系數(shù)聯(lián)合起來一并求解，可省去諸多中間環(huán)節(jié)，計算量大為減小，另外該算法便于求解復(fù)雜區(qū)域上的問題，對于事先未知的自由邊界或求解區(qū)域內(nèi)部不同介質(zhì)的交界面比較容易處理。本文將其拓展到紊流情況下的求解，從而利用有限變分法獲得高速下徑向滑動軸承層紊流下非線性油膜特性的一種快速計算方法，并通過短軸承理論公式和常用的有限差分法進行驗證。

本文針對高速轉(zhuǎn)子—滑動軸承系統(tǒng)從開始工作到穩(wěn)定在平衡位置這段時間的非線性油膜動力特性系數(shù)的求解問題，提出了高速下徑向滑動軸承非線性油膜動力特性系數(shù)動態(tài)分析方法。該方法考慮任意時刻轉(zhuǎn)子軸心瞬態(tài)響應(yīng)位移及速度與滑動軸承非線性油膜動力特性系數(shù)之間的耦合影響，采用有限變分法和Newmark－β法實現(xiàn)滑動軸承非線性油膜動力特性系數(shù)與轉(zhuǎn)子瞬態(tài)響應(yīng)之間的耦合，并以Jeffcott轉(zhuǎn)子瞬態(tài)響應(yīng)為例，對比研究了傳統(tǒng)靜態(tài)與本文動態(tài)非線性油膜動力特性系數(shù)分析方法。

1 基于有限變分法的油膜動力特性系數(shù)的計算

1.1 紊流瞬態(tài)下油膜變分方程

圖1 軸瓦及其平面展開

長徑比為

油膜流速為v時，雷諾數(shù)

紊流瞬態(tài)雷諾量綱一形式［8］如下：

其中，K1、K2為紊流潤滑系數(shù)，表達式為

油膜很薄時一般認為油膜厚度方向壓力不變，沿油膜厚度方向做功微小，所以流速的變化主要由速度流引起，層紊流區(qū)域用雷諾數(shù)來區(qū)分，臨界雷諾數(shù)為1000［8］。K1＝K2＝12時即可表征層流狀態(tài)下的油膜特性。

鄭鐵生等［9］指出，二維變分法求解Reynolds邊界條件下的油膜壓力分布等價于求解如下算子方程：

算子L的值域為油膜潤滑區(qū)域Ω：

壓力分布函數(shù)q為算子方程準確解的充要條件是：當(dāng)p＝q時，泛函J（q）取極小值。J（q）的表達式為

由圖1所示，軸瓦潤滑域為矩形，油膜破裂邊雖然微曲，但其上壓力及壓力導(dǎo)數(shù)皆為零，設(shè)其為直線，對總體油膜力的影響微乎其微，故壓力分布函數(shù)q可分離變量為

引入變分Euler方程：

從而壓力函數(shù)可以通過參數(shù)Kc［6］和函數(shù)s（θ）的相互迭代修正得到，壓力函數(shù)的求解轉(zhuǎn)化為函數(shù)s（θ）的求解，二維方程降為一維方程。壓力函數(shù)用Ritz法對s（θ）采用分段線性有限元插值函數(shù)向量r（θ）來處理，采用修正的追趕法［6］求解。油膜壓力的計算轉(zhuǎn)化為下式線性矩陣方程中向量a的求解：

1.2 系數(shù)矩陣分解

求解式（3）可獲得紊流瞬態(tài)下的油膜壓力特性，但其中的系數(shù)矩陣A涉及A1和A2的求解。求解時既有對角度θ的積分，還與油膜厚度有關(guān)，油膜厚度不僅與角度θ有關(guān)，也是軸心位移的函數(shù)。將有關(guān)角度θ的積分運算和有關(guān)軸心位移的運算分離開來，每次軸心位移變化時只需按疊加處理，不必要重新積分，可大大縮短計算時間。由于

則A1及A2轉(zhuǎn)化為

其中矩陣Gi及Bi分別為

1.3 油膜力的求解

作用于軸瓦的量綱一局部油膜力（在X、Y方向上）分別為

化為整體坐標系下的油膜力分別為

式中，φ為偏位角。

1.4 油膜剛度阻尼的計算

量綱一剛度矩陣K及量綱一阻尼矩陣C為

2 轉(zhuǎn)子—滑動軸承系統(tǒng)動態(tài)分析方法

任意時刻轉(zhuǎn)子在外激勵作用下的運動微分方程為

式中，u為轉(zhuǎn)子軸頸軸心位移；M為質(zhì)量矩陣；D為阻尼矩陣（含滑動軸承阻尼矩陣C）；B為剛度矩陣（含滑動軸承剛度矩陣K）；F為激振力。

本文采用Newmark—β法［10］分析轉(zhuǎn)子系統(tǒng)的瞬態(tài)響應(yīng)。

任意時刻，只要已知滑動軸承處軸頸中心的位移、速度，就可以通過有限變分法計算滑動軸承非線性油膜動力特性系數(shù)，同時只要知道滑動軸承的非線性油膜動力特性系數(shù)就可以通過Newmark—β法分析轉(zhuǎn)子系統(tǒng)的瞬態(tài)響應(yīng)（包括滑動軸承處軸頸中心的位移、速度）。通過以上的分析可知，轉(zhuǎn)子瞬態(tài)響應(yīng)分析與滑動軸承的非線性油膜動力特性系數(shù)分析是緊密耦合在一起的。所以本文將時間進行離散（時間步長為轉(zhuǎn)子瞬態(tài)計算的時間步長），在每個時刻對非線性油膜動力特性系數(shù)進行迭代分析。轉(zhuǎn)子—滑動軸承系統(tǒng)動態(tài)分析方法計算流程如圖2所示。

圖2 轉(zhuǎn)子—滑動軸承系統(tǒng)動態(tài)分析方法

t時刻收斂所得的非線性油膜動力特性系數(shù)矩陣K和C將作為t＋dt時刻的初值進入動態(tài)分析中；同理，t＋dt時刻非線性油膜剛度矩陣K′和阻尼矩陣C′將作為t＋2dt時刻的初值進入動態(tài)分析中。

非線性油膜動力特性系數(shù)的收斂條件為

式中，ε為迭代收斂需要滿足的精度。

3 計算實例

以某Jeffcott轉(zhuǎn)子系統(tǒng)為例進行動態(tài)分析，其轉(zhuǎn)子、滑動軸承參數(shù)如表1所示。

表1 Jeffcott轉(zhuǎn)子系統(tǒng)參數(shù)

3.1 非線性油膜動力特性系數(shù)的計算

表2 滑動軸承油膜力和動力特性系數(shù)計算結(jié)果對比

3.2 Jeffcott轉(zhuǎn)子系統(tǒng)動態(tài)分析

軸的彈性模量為200GPa，圓盤處的外阻尼系數(shù)為0.01kg／s。Newmark—β 法控制參數(shù)α ＝0.5，β＝0.25，動態(tài)分析過程的時間步長選為0.1s，總時長為12s。

圖3～圖5所示分別為采用本文動態(tài)滑動軸承動力特性系數(shù)方法和靜態(tài)滑動軸承動力特性系數(shù)獲得的12s時轉(zhuǎn)子軸承處軸心的時域響應(yīng)和軸心軌跡。與采用靜態(tài)滑動軸承動力特性系數(shù)求得的轉(zhuǎn)子軸心時域響應(yīng)和軸心軌跡相比，兩種方法得出的變化趨勢一致，但采用動態(tài)滑動軸承動力特性系數(shù)方法獲得的瞬態(tài)響應(yīng)幅值更小，故采用靜態(tài)滑動軸承動力特性系數(shù)設(shè)計的軸承更偏于安全，而動態(tài)分析方法精度更高。

圖3 12s時轉(zhuǎn)子軸承處軸心的x方向時域響應(yīng)

圖5 12s時轉(zhuǎn)子軸承處軸心軌跡

從圖3～圖5可知，轉(zhuǎn)子啟動后，在非線性油膜力的作用下，軸頸處軸心位移隨時間呈周期性擾動，但最終會穩(wěn)定在平衡位置。

3.3 滑動軸承動力特性系數(shù)時域響應(yīng)

圖6描述了滑動軸承動力特性系數(shù)隨時間的變化規(guī)律。從圖6可以看出，在開始時間段內(nèi)，剛度系數(shù)和阻尼系數(shù)變化很大；在啟動最初時間段內(nèi)，轉(zhuǎn)子還沒有達到平衡位置，軸心渦動引起了系數(shù)的較大變化；當(dāng)穩(wěn)定在平衡位置后，滑動軸承動力特性系數(shù)也就趨于常數(shù)。

4 結(jié)語

轉(zhuǎn)子—滑動軸承系統(tǒng)非線性動力學(xué)分析的一項核心內(nèi)容是油膜動態(tài)特性的計算。短軸承模型實用性差；有限元或有限差分法雖然精度較高，但對于非線性分析所需的大規(guī)模數(shù)值算法而言，計算量大，求解時間長。因此，建立一種快速而實用的方法極為重要。本文通過有限變分法和Newmark—β法實現(xiàn)了滑動軸承動力特性系數(shù)分析與轉(zhuǎn)子瞬態(tài)響應(yīng)之間的流固耦合，通過學(xué)科迭代方法獲得非線性油膜動態(tài)特性。與傳統(tǒng)的靜態(tài)分析方法相比，瞬態(tài)響應(yīng)幅值較小，精度更高。算例表明，在確保精度的情況下，有限變分法可使運算時間大幅度縮短。

圖6 滑動軸承動力特性系數(shù)時域響應(yīng)

［1］Lund J W.The Sstability of an Elastic Rotor in Journal Bearing with Flexible Damped Supports［J］.ASME Journal of Applied Mechanics，1965，32（4）：911－920.

［2］龐曉平，牛坤，魏思遠.用CFD方法研究具有通用膜厚方程徑向軸承的動壓特性［J］.中國機械工程，2012，23（21）：2521－2524.Pang Xiaoping，Niu Kun，Wei Siyuan.Research on Hydrodynamic Lubrication Performance of Journal Bearings Based on General Film Thickness Equation Using CFD Software［J］.China Mechanical Engineering，2012，23（21）：2521－2524.

［3］楊金福，楊昆，于迭仁.非線性動態(tài)油膜力及π油膜力分解的特性研究［J］.動力工程，2005，25（6）：770－779.Yang Jinfu，Yang Kun，Yu Daren.Characteristics of Nonlinear Dynamic Oil FiIm Forces and the Decomposition ofπOil FiIm Forces［J］.Chinese Journal of Power Engineering，2005，25（6）：770－779.

［4］李元生，敖良波，李磊，等.滑動軸承動力特性系數(shù)動態(tài)分析方法［J］.機械工程學(xué)報，2010，46（21）：48－53.Li Yuansheng，Ao Liangbo，Li Lei，et al.Dynamic Analysis Method of Dynamic Character Coefficient of Hydrodynamic Journal Bearing［J］.Chinese Journal of Mechanical Engineering，2010，46（21）：48－53.

［5］王小靜，蘇葒，張直明.基于復(fù)合型紊流潤滑理論的徑向滑動軸承紊流潤滑性能的研究［J］.機械工程學(xué)報，2003，39（3）：85－89.Wang Xiaojing，Su Hong，Zhang Zhiming.Study on the Performance of Journal Bearing in Turbulent Regime with the Combined Reynolds Stress Model［J］.Chinese Journal of Mechanical Engineering，2003，39（3）：85－89.

［6］肖忠會，劉大全，張文，等.發(fā)電機組軸承非線性油膜力的快速算法研究［J］.中國電機工程學(xué)報，2004，24（12）：144－148.Xiao Zhonghui，Liu Daquan，Zhang Wen，et al.Study on a Rapid Algorithm for Nonlinear Oil－film Forces of Journal Bearings in Turbo－generator Unit［J］.Proceedings of the CSEE，2004，24（12）：144－148.

［7］肖忠會，王麗萍，鄭鐵生.滑動軸承油膜力Jacobi矩陣的一種快速算法［J］.應(yīng)用力學(xué)學(xué)報，2005（1）：67－70.Xiao Zhonghui，Wang Liping，Zheng Tiesheng.A Rapid Algorithm for Jacobi Matrix Oil－film Forces of Journal Bearings［J］.Chinese Jouranl of Applied Mechanics，2005（1）：67－70.

［8］張直明.滑動軸承的流體動力潤滑理論［M］.北京：高等教育出版社，1986.

［9］鄭鐵生，張文，馬建敏，等.滑動軸承非線性油膜力的幾個理論問題及應(yīng)用［J］.水動力學(xué)研究與進展，2003，18（6）：761－768.Zheng Tiesheng，Zhang Wen，Ma Jianmin，et al.Some Theoretical Problems on the Nonlinear Oilfilm Forces of Sliding Bearing and Their Applications［J］.Joural of Hydrodynamica，2003，18（6）：761－768.

［10］徐斌.MATLAB有限元結(jié)構(gòu)動力學(xué)分析與工程應(yīng)用［M］.北京：清華大學(xué)出版社，2010.

［11］顧致平，劉永壽.非線性轉(zhuǎn)子系統(tǒng)中的傳遞矩陣技術(shù)［M］.北京：科學(xué)出版社，2010.