基于流固耦合的滑動(dòng)軸承非線性油膜動(dòng)特性研究

毛文貴 韓 旭 劉桂萍

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

0 引言

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

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

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

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

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

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

圖1 軸瓦及其平面展開

長徑比為

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

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

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

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

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

算子L的值域?yàn)橛湍櫥瑓^(qū)域Ω：

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

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

引入變分Euler方程：

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

1.2 系數(shù)矩陣分解

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

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

其中矩陣Gi及Bi分別為

1.3 油膜力的求解

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

化為整體坐標(biāo)系下的油膜力分別為

式中，φ為偏位角。

1.4 油膜剛度阻尼的計(jì)算

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

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

任意時(shí)刻轉(zhuǎn)子在外激勵(lì)作用下的運(yùn)動(dòng)微分方程為

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

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

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

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

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

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

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

3 計(jì)算實(shí)例

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

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

3.1 非線性油膜動(dòng)力特性系數(shù)的計(jì)算

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

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

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

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

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

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

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

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

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

4 結(jié)語

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

圖6 滑動(dòng)軸承動(dòng)力特性系數(shù)時(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］龐曉平，牛坤，魏思遠(yuǎn).用CFD方法研究具有通用膜厚方程徑向軸承的動(dòng)壓特性［J］.中國機(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］楊金福，楊昆，于迭仁.非線性動(dòng)態(tài)油膜力及π油膜力分解的特性研究［J］.動(dòng)力工程，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］李元生，敖良波，李磊，等.滑動(dòng)軸承動(dòng)力特性系數(shù)動(dòng)態(tài)分析方法［J］.機(jī)械工程學(xué)報(bào)，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ù)合型紊流潤滑理論的徑向滑動(dòng)軸承紊流潤滑性能的研究［J］.機(jī)械工程學(xué)報(bào)，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］肖忠會(huì)，劉大全，張文，等.發(fā)電機(jī)組軸承非線性油膜力的快速算法研究［J］.中國電機(jī)工程學(xué)報(bào)，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］肖忠會(huì)，王麗萍，鄭鐵生.滑動(dòng)軸承油膜力Jacobi矩陣的一種快速算法［J］.應(yīng)用力學(xué)學(xué)報(bào)，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］張直明.滑動(dòng)軸承的流體動(dòng)力潤滑理論［M］.北京：高等教育出版社，1986.

［9］鄭鐵生，張文，馬建敏，等.滑動(dòng)軸承非線性油膜力的幾個(gè)理論問題及應(yīng)用［J］.水動(dòng)力學(xué)研究與進(jìn)展，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)動(dòng)力學(xué)分析與工程應(yīng)用［M］.北京：清華大學(xué)出版社，2010.

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