基于精細(xì)積分法的軸承-拉桿轉(zhuǎn)子系統(tǒng)非線性動(dòng)力學(xué)分析

黑棣，鄭美茹

（陜西鐵路工程職業(yè)技術(shù)學(xué)院 機(jī)電工程系，陜西 渭南 714000）

軸承-轉(zhuǎn)子系統(tǒng)作為旋轉(zhuǎn)機(jī)械的核心部件，廣泛應(yīng)用于航天、電力、工程等領(lǐng)域。許多學(xué)者針對(duì)軸承-轉(zhuǎn)子系統(tǒng)的相關(guān)問題進(jìn)行了研究。文獻(xiàn)［1-2］研究了表面織構(gòu)對(duì)滑動(dòng)軸承性能的影響。文獻(xiàn)［3］求解了Reynolds方程，并且得到油膜壓力分布曲線，分析了各種參數(shù)對(duì)油膜力的影響。文獻(xiàn)［4］利用變分原理和分離變量法求得了油膜力的近似解析解。

以上研究都是為了能夠更準(zhǔn)確地進(jìn)行軸承-轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)分析。因?yàn)橄到y(tǒng)的穩(wěn)定性直接決定旋轉(zhuǎn)機(jī)械運(yùn)行的穩(wěn)定性，因此，軸承-轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)分析計(jì)算十分關(guān)鍵，最常用的軸承-轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)計(jì)算方法有Runge-Kutta法［5-6］、Wilson-θ法［7-8］、Newmark法［9-10］等，但均對(duì)計(jì)算步長(zhǎng)的依賴較大。

精細(xì)積分法對(duì)步長(zhǎng)的依賴較小，而且能保證較高的精度。文獻(xiàn)［11］基于2N類算法提出了精細(xì)積分法。文獻(xiàn)［12］將Wilson-θ法和精細(xì)積分法結(jié)合分析了結(jié)構(gòu)地震時(shí)程反應(yīng)。文獻(xiàn)［13］給出了單步4階精度的精細(xì)積分法，該方法精度高，計(jì)算量小。文獻(xiàn)［14］采用改進(jìn)的精細(xì)積分法計(jì)算了轉(zhuǎn)子系統(tǒng)的非線性動(dòng)力學(xué)行為。文獻(xiàn)［15］將數(shù)據(jù)庫方法與精細(xì)時(shí)程積分法相結(jié)合，研究了可傾瓦軸承-轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)響應(yīng)。

現(xiàn)基于改進(jìn)的精細(xì)積分法，將其與傳統(tǒng)的Wilson-θ法、Newmark法進(jìn)行比較，求解轉(zhuǎn)子系統(tǒng)的非線性動(dòng)力學(xué)響應(yīng)，研究步長(zhǎng)對(duì)轉(zhuǎn)子系統(tǒng)的非線性動(dòng)力學(xué)行為的影響。

1 轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)方程

拉桿轉(zhuǎn)子系統(tǒng)模型如圖1所示。圖中：O1和O2為2個(gè)圓盤的中心；mO1，mO2分別為兩圓盤的質(zhì)量；eO1，eO2分別為兩圓盤的偏心量；l為轉(zhuǎn)軸的總長(zhǎng)；a和b為兩轉(zhuǎn)軸的長(zhǎng)度；AO1＝a1，BO1＝b1，AO2＝a2，BO2＝b2；ma，mb分別為 AO，BO的集中質(zhì)量；fx和 fy分別為軸承處x，y方向上的油膜力。

圖1 軸承-拉桿轉(zhuǎn)子系統(tǒng)示意圖Fig.1 Sketch of bearing-rod fastening rotor system

式中：M為質(zhì)量矩陣；G為陀螺矩陣；K為剛度矩陣；F為非線性油膜力向量；FAx，F(xiàn)Ay-分別為軸承A處轉(zhuǎn)子x，y方向的無量綱非線性油膜力；FBx，F(xiàn)By-分別為軸承B處轉(zhuǎn)子x，y方向的無量綱非線性油膜力；Q為外激勵(lì)列向量；W為重力列向量；FR為圓盤間的非線性回復(fù)力；KR為量綱一的非線性回復(fù)剛度；Kb為拉桿量綱一的彎曲剛度。q為位移列向量；XA，YA，XO1，YO1，XO2，YO2，XB，YB分別為轉(zhuǎn)子在軸承A端、圓盤O1處、圓盤O2處、軸承B端 x和 y方向的位移；θx1，θy1，θx2，θy2分別為圓盤1，2繞x和y軸的擺角；Jd1，Jd2分別為圓盤1，2的極轉(zhuǎn)動(dòng)慣量；Jz1，Jz2分別為圓盤1，2的赤道轉(zhuǎn)動(dòng)慣量；k11，k22分別為AO軸段x，y方向的彎曲剛度；k′11，k′22分別為 OB軸段 x，y方向的彎曲剛度；k33，k44，k′33，k′44分別為由于圓盤 1，2繞 x和 y軸擺動(dòng)而引起的軸段剛度；k14，k23為AO軸段的交叉剛度；k′14，k′23為 OB軸段的交叉剛度，由于轉(zhuǎn)子的對(duì)稱性，有 k11＝k22＝k′11＝k′22，k33＝k44＝k′33＝k′44，k14＝k23，k′14＝k′［16］23；kR為非線性回復(fù)剛度；kb為拉桿的彎曲剛度。

2 改進(jìn)的精細(xì)積分法

2.1 算法介紹

運(yùn)用精細(xì)積分法求解轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)響應(yīng)，將（1）式化為

2.2 算法驗(yàn)證

為了驗(yàn)證該算法的有效性，針對(duì)以下算例，將其與Newmark法、Wilson-θ法進(jìn)行比較。

分別采用3種算法進(jìn)行計(jì)算，計(jì)算結(jié)果與精確解見表1。

由表1可以看出，文中方法的計(jì)算結(jié)果比Newmark法、Wilson-θ法更接近精確解。

表1 3種方法計(jì)算結(jié)果Tab.1 Calculation results of three kinds of methods

3 數(shù)值算例

軸承-對(duì)稱拉桿轉(zhuǎn)子系統(tǒng)參數(shù)如下。軸承為360°徑向滑動(dòng)軸承，軸承寬度和直徑之比D／d＝1，μ＝0.022 Pa·s，c＝0.000 49 m。轉(zhuǎn)軸參數(shù)：r＝0.14 m，l＝3.3m，a＝1.6m，b＝1.6m，a1＝1.4 m，b1＝1.8 m，a2＝1.8m，b2＝1.4m，mA＝mB＝768.458 7 kg。圓盤參數(shù)：h＝0.4 m，R＝0.3 m，eO1＝eO2＝0.01c，mO1＝mO2＝690.044 5 kg，Jz1＝Jz2＝31.052 kg·m2，Jd1＝Jd2＝15.526 kg·m2。2m＝mA+mB+mO1+mO2，m＝1 458.5 kg。轉(zhuǎn)軸剛度：k11＝k22＝k′11＝k′22＝9.992 5×107N／m，k14＝k23＝ -3.758 4×107N／m，k′14＝k′23＝3.758 4×106N／m，k33＝k44＝k′33＝k′44＝2.367 8×108N／m，kb＝kR＝2 k11。

以轉(zhuǎn)子的量綱一轉(zhuǎn)速為控制參數(shù)分析轉(zhuǎn)子的非線性動(dòng)力學(xué)行為，其計(jì)算步長(zhǎng)為 h＝2π／200。當(dāng)量綱一轉(zhuǎn)速的變化范圍為0.9～1.05時(shí)，轉(zhuǎn)子軸承處和圓盤處y方向位移隨轉(zhuǎn)速的運(yùn)動(dòng)分岔圖如圖2所示。

圖2 轉(zhuǎn)子軸承和圓盤處y方向位移隨轉(zhuǎn)速的運(yùn)動(dòng)分岔圖Fig.2 Bifurcation diagram of y versus at bearing station and disk stations

由圖2可以看出，轉(zhuǎn)速較低時(shí)，轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)為周期運(yùn)動(dòng)，隨著轉(zhuǎn)速的升高，轉(zhuǎn)子系統(tǒng)將發(fā)生分岔行為。

圖3 ＝0.92時(shí)轉(zhuǎn)子軸承處的軌跡圖、Poincaré映射圖、時(shí)間歷程、頻譜圖和圓盤1處的軌跡圖Fig.3 The orbit，Poincaré map，time series，spectrum of the rotor at bearing station and the orbit of disk 1 for＝0.92

圖4 ＝0.92時(shí)圓盤分別繞x，y軸的擺角Fig.4 The swing angle of disks around the x and y axis for＝0.92

圖5 ＝0.95時(shí)轉(zhuǎn)子軸承處的軌跡圖、Poincaré映射圖、時(shí)間歷程、頻譜圖和圓盤1處的軌跡圖Fig.5 The orbit，Poincaré map，time series，spectrum of the rotor at bearing station and the orbit of disk 1 for＝0.95

圖6 ＝0.95時(shí)圓盤分別繞x，y軸的擺角Fig.6 The swing angle of disks around the x and y axis for＝0.95

圖7 ＝1.04時(shí)轉(zhuǎn)子軸承處的軌跡圖、Poincaré映射圖、時(shí)間歷程和頻譜圖Fig.7 The orbit，Poincaré map，time series and spectrum of the rotor at bearing station for＝1.04

圖8 ＝1.04時(shí)圓盤繞x，y軸的擺角Fig.8 The swing angle of disks around the x and y axis for＝1.04

以上計(jì)算步長(zhǎng)取為h＝2π／200，為了研究計(jì)算步長(zhǎng)對(duì)系統(tǒng)非線性動(dòng)力學(xué)行為的影響，以下步長(zhǎng)分別取為 h＝2π／400，2π／600，轉(zhuǎn)子軸承處和圓盤處的分岔圖如圖9、圖10所示。對(duì)比圖2、圖9和圖10可知，轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)的變化趨勢(shì)基本相同，但分岔點(diǎn)有所不同，隨著步長(zhǎng)的減小，分岔點(diǎn)有少許增大，能夠考察的轉(zhuǎn)速范圍也增大，且轉(zhuǎn)子系統(tǒng)能夠表現(xiàn)出更豐富的非線性動(dòng)力學(xué)現(xiàn)象。由此可以看出，計(jì)算步長(zhǎng)的變化對(duì)于轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)具有一定影響。

圖9 h＝2π／400時(shí)轉(zhuǎn)子軸承和圓盤處y方向位移隨轉(zhuǎn)速的運(yùn)動(dòng)分岔圖Fig.9 Bifurcation diagram of y versus at bearing station and disk stations for h＝2π／400

圖10 h＝2π／600時(shí)轉(zhuǎn)子軸承和圓盤處y方向位移隨轉(zhuǎn)速的運(yùn)動(dòng)分岔圖Fig.10 Bifurcation diagram of y versus at bearing station and disk stations for h＝2π／600

圖12 ＝1.05時(shí)3種步長(zhǎng)下軸承處軌跡Fig.12 Comparison between the three steps at bearing station for＝1.05

根據(jù)計(jì)算結(jié)果可以看出，計(jì)算步長(zhǎng)對(duì)非周期運(yùn)動(dòng)影響較大，故從分岔點(diǎn)附近開始要縮小計(jì)算步長(zhǎng)。在周期運(yùn)動(dòng)階段，取不同的步長(zhǎng)轉(zhuǎn)子的運(yùn)動(dòng)軌跡有一定差別，但可以適當(dāng)加大計(jì)算步長(zhǎng)。

4 結(jié)束語

提出了一種改進(jìn)的精細(xì)積分法，將其與常用的Newmark法、Wilson-θ法進(jìn)行比較，發(fā)現(xiàn)該方法的計(jì)算結(jié)果比Newmark法和Wilson-θ法更接近精確解，驗(yàn)證了其可靠性。

運(yùn)用改進(jìn)的精細(xì)積分法計(jì)算拉桿轉(zhuǎn)子系統(tǒng)的非線性動(dòng)力學(xué)響應(yīng)，從計(jì)算結(jié)果看出，該拉桿轉(zhuǎn)子系統(tǒng)存在豐富的非線性運(yùn)動(dòng)現(xiàn)象。同時(shí)研究了拉桿轉(zhuǎn)子圓盤的擺動(dòng)，當(dāng)轉(zhuǎn)子系統(tǒng)處于準(zhǔn)周期運(yùn)動(dòng)的情況下，圓盤的擺角隨時(shí)間發(fā)散，最終導(dǎo)致轉(zhuǎn)子碰壁。而其他情況下，圓盤的擺角均隨時(shí)間呈周期性變化。

分析了計(jì)算步長(zhǎng)對(duì)轉(zhuǎn)子系統(tǒng)運(yùn)動(dòng)的影響，從計(jì)算結(jié)果可看出，計(jì)算步長(zhǎng)對(duì)于轉(zhuǎn)子系統(tǒng)非周期運(yùn)動(dòng)階段的影響較大，故在非周期階段應(yīng)當(dāng)采用較小的計(jì)算步長(zhǎng)。在周期運(yùn)動(dòng)階段，即使其他參數(shù)都相同，采用不同的計(jì)算步長(zhǎng)，轉(zhuǎn)子的運(yùn)動(dòng)軌跡仍有一定差別。計(jì)算步長(zhǎng)對(duì)轉(zhuǎn)子系統(tǒng)的影響有待進(jìn)一步深入研究。