基于流固耦合的滑動軸承非線性油膜動特性研究

毛文貴 韓 旭 劉桂萍

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

0 引言

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

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

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

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

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

1.1 紊流瞬態下油膜變分方程

圖1 軸瓦及其平面展開

長徑比為

油膜流速為v時，雷諾數

紊流瞬態雷諾量綱一形式［8］如下：

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

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

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

算子L的值域為油膜潤滑區域Ω：

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

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

引入變分Euler方程：

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

1.2 系數矩陣分解

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

則A1及A2轉化為

其中矩陣Gi及Bi分別為

1.3 油膜力的求解

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

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

式中，φ為偏位角。

1.4 油膜剛度阻尼的計算

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

2 轉子—滑動軸承系統動態分析方法

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

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

本文采用Newmark—β法［10］分析轉子系統的瞬態響應。

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

圖2 轉子—滑動軸承系統動態分析方法

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

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

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

3 計算實例

以某Jeffcott轉子系統為例進行動態分析，其轉子、滑動軸承參數如表1所示。

表1 Jeffcott轉子系統參數

3.1 非線性油膜動力特性系數的計算

表2 滑動軸承油膜力和動力特性系數計算結果對比

3.2 Jeffcott轉子系統動態分析

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

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

圖3 12s時轉子軸承處軸心的x方向時域響應

圖5 12s時轉子軸承處軸心軌跡

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

3.3 滑動軸承動力特性系數時域響應

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

4 結語

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

圖6 滑動軸承動力特性系數時域響應

［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］楊金福，楊昆，于迭仁.非線性動態油膜力及π油膜力分解的特性研究［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］李元生，敖良波，李磊，等.滑動軸承動力特性系數動態分析方法［J］.機械工程學報，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］王小靜，蘇葒，張直明.基于復合型紊流潤滑理論的徑向滑動軸承紊流潤滑性能的研究［J］.機械工程學報，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］肖忠會，劉大全，張文，等.發電機組軸承非線性油膜力的快速算法研究［J］.中國電機工程學報，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］.應用力學學報，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］鄭鐵生，張文，馬建敏，等.滑動軸承非線性油膜力的幾個理論問題及應用［J］.水動力學研究與進展，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有限元結構動力學分析與工程應用［M］.北京：清華大學出版社，2010.

［11］顧致平，劉永壽.非線性轉子系統中的傳遞矩陣技術［M］.北京：科學出版社，2010.