轉子-葉片系統非線性振動和動態穩定性分析

李炳強， 馬 輝,2， 曾 勁， 郭旭民， 崔 璨

(1. 東北大學 機械工程與自動化學院, 沈陽 110819； 2.東北大學 航空動力裝備振動及控制教育部重點實驗室，沈陽 110819)

旋轉系統包括軸、葉片和軸承組件等，應用于許多工程機械，如壓縮機、汽輪機、航空發動機等。在轉子-葉片模型中，軸承通常假定為剛性或帶有黏性阻尼的線性彈簧。隨著旋轉機械結構變得更薄，工作轉速更高，可能出現高振幅振動甚至失穩等危險現象，因此對轉子-葉片系統的穩定性和分岔進行研究是必要的。

考慮幾何非線性以及振動大變形，Khadem等[1]采用多尺度法研究簡支軸在主共振條件下的自由和受迫振動。此外，基于Hamilton原理和諧波平衡法，他們還研究了簡支旋轉軸的模態組合共振[2]。Shahgholi 等[3]考慮到軸的大變形和拉伸效應所引入的非線性，分析了旋轉軸在轉速波動影響下，主共振頻率附近的非線性振動和失穩現象。Chiu等[4]研究了軸扭轉和葉片彎曲振動對多盤轉子系統耦合非線性振動的影響。Wang等[5]建立了柔性轉子-葉片系統的時變非線性模型，討論了轉子的非線性振動對葉片響應的影響。Sanches等[6]研究了直升機螺旋槳的基礎共振現象，基于Floquet理論確定了該非線性系統分叉點的位置，描述了在不同條件下的失穩區域。Bab等[7]采用多尺度方法研究了柔性轉子和柔性/剛性葉片的耦合非線性共振，討論了圓盤偏心量和阻尼對系統失穩的影響。

碰摩力可以使轉子-葉片系統發生失穩。Parent等[8]研究了渦扇發動機風扇葉片-機匣的碰摩響應，預測了系統的失穩區域，討論了系統在失穩頻率區間的響應。高喆等[9]研究了碰摩轉子在隨機擾動下呈現出的非線性動力學特性。陶海亮等[10]建立了轉子-軸承系統的碰摩故障模型，分析了碰摩剛度以及靜子質量對轉子動力學響應的影響。劉昕等[11]采用Hertz接觸模型構建了葉片-機匣碰摩力模型，并運用該模型研究了轉子-軸承碰摩的非線性動力學響應。碰摩力存在多種形式，對于點碰摩或局部碰摩，碰摩力被認為是周期性脈沖力。Petrov[12-13]應用多諧波的非線性碰摩力研究轉子的響應，Sinha[14-15]提出若干碰摩脈沖力的數學表達式，如正弦半波脈沖，矩形脈沖等。Turner等[16-17]采用類半正弦波脈沖力研究葉片-機匣的碰摩。Kou等[18]應用正弦脈沖函數和連續正弦函數模擬葉片-機匣的法向碰摩力。Ma等[19]確定了周期脈沖碰摩力最大值的表達形式。

以上文獻中，有些學者將軸簡化為簡支梁，不考慮軸承支承剛度和阻尼的影響，另外一些學者雖然將軸承考慮為線彈性支承，但是沒有考慮支承剛度和陀螺效應對轉軸振型的影響。基于此，本文研究了主共振條件下轉子-葉片系統受周期性脈沖碰摩力沖擊的振動響應和穩定性，分析討論了軸承支承剛度、阻尼、陀螺效應、碰摩力、摩擦因數和圓盤偏心量等對系統響應幅值以及失穩區域的影響，并通過引入數值積分結果驗證本文攝動解的準確性。

1 轉子-葉片系統運動微分方程

圖1為轉子-葉片系統示意圖，其中X-Y-Z為慣性坐標系，X為未變形軸的中性軸方向，x-y-z為動態坐標系，附著于軸變形后的截面主軸上，xb-yb-zb為葉片的動態坐標系，xb沿著葉片的寬度方向，yb沿著葉片的厚度方向，zb沿著葉片的長度方向。軸左端的扭簧支承kφ=∞，也就意味著繞X軸的扭角為0，同時固定軸左端X方向的位移，kbz1,kbz2,kby1,kby2分別為軸承1和軸承2沿著Z方向和Y方向的支承剛度，cbz1,cbz2,cby1,cby2分別為軸承1和軸承2沿著Z方向和Y方向的阻尼，為了計算的簡化，令kbz1=kbz2=kby1=kby2=k,cbz1=cbz2=cby1=cby2=cbearing。旋轉軸的動能和勢能可表示為

圖1 轉子-葉片系統示意圖Fig.1 Schematic of rotor-blade systems

(1)

(2)

式中：l，m，I1,I2,I3,N11,α,D11,D22,D33分別為軸的長度、單位長度質量、單位長度極轉動慣量、單位長度關于y軸的轉動慣量、單位長度關于z軸的轉動慣量、軸向剛度、軸向應變、扭轉剛度、關于y軸的彎曲剛度、關于z軸的彎曲剛度;u,v,w,ω1,ω2和ω3為軸上任意點在慣性坐標系下沿著X,Y,Z三個方向的位移和角速度；“·”為對t求偏導；k1，k2，k3分別為截面的扭轉曲率、截面關于y軸的彎曲曲率和截面關于z軸的彎曲曲率。

盤假設為剛性，固定于軸上x=xd的位置，其勢能為0，動能可以表示為[20]

(3)

式中:mdisk,Ω,ey,ez,Idisk,Jdisk為盤的質量、旋轉角速度、圓盤關于Y方向的偏心距、圓盤關于Z方向的偏心距、直徑轉動慣量和極轉動慣量。

將葉片簡化為Euler-Bernoulli梁，對于第i個葉片，扭轉位移φi(t)表達為

φi(t)=Ωt+φ(xd,t)+?i,i=1, …,Nb

(4)

式中：Nb為葉片數;φ(x,t)為軸的扭轉變形;?i為第i個葉片的角位置

(5)

葉片的安裝角為γ,θ(x,t)和ψ(x,t)分別為軸繞Y軸和Z軸旋轉的兩個歐拉角，xb,Λi,yb,Rd和ζ分別為任意質點在橫截面內沿著葉片寬度方向的位置、葉片的彎曲振動幅值、葉片在橫截面內沿著厚度方向的位置、圓盤半徑、質點與葉根的距離，考慮到葉片動態坐標系向慣性坐標系的轉化，其上任意質點的位置向量為

(6)

則第i個葉片的動能為

(7)

(8)

對所有葉片，總體動能和勢能為

(9)

在葉尖處施加周期脈沖碰摩力，其表達式為

(10)

式中：Fn max,tc,tp,t0為最大法向碰摩力，接觸時間、周期以及開始碰摩時間。設定μ為摩擦因數，則第i個葉片的碰摩力所做的功為

WFni=μFniΛi(L,t)cosγ-μFni(Rd+L)φi(xd+t)+

Fni(-cosφi+μsinφi)vδ(x-xd)+
Fni(-sinφi-μcosφi)wδ(x-xd)

(11)

Hamilton原理的表達式為

(12)

其中,

在式(2)中，軸向應變的表達式為

(13)

式中：“′”為對x求偏導。

由圖1中支承方式可知，固定左端軸向位移，而右端可以沿著軸向移動，即α=0，可得

(14)

軸上任意點繞Y軸和Z軸旋轉的兩個歐拉角θ(x,t)和ψ(x,t)可以表達為

(15)

轉子-葉片系統的響應峰值和失穩區域均發生在彎曲共振頻率范圍內，而軸的扭轉頻率要遠高于彎曲頻率，因此，相較于彎曲慣性項和剛度項，扭轉慣性項的值很小，可以忽略。如果v=O(ε)且w=O(ε), 其中ε是無量綱的小尺度參數，則扭轉運動的非線性項數量級為O(ε4), 可得

(16)

將式(14)～式(16)代入式(1)～式(6)，可消掉u(x,t),φ(x,t)，θ(x,t)和ψ(x,t)，得到系統關于w(x,t),v(x,t)和Λi(ζ,t)的非線性運動微分方程。采用Coleman變換縮減系統自由度，引入變量ξ和η，表達式為

(17)

為了進一步對系統展開降維，采用如下復平面映射

(18)

通過以上變換，方程的數量由2+Nb降為2個，設定軸承，軸和葉片的阻尼系數分別為cbearing,c和cblade。進一步引入如下無量綱變量

(19)

值得注意的是，在攝動求解中，cbearing,c和cblade數量級為ε2，而Fni,ey和ez的數量級為ε3，其中ε是無量綱的小參數。綜合式(1)～式(19)得到轉子-葉片系統的無量綱運動微分方程為

(20)

(21)

式(20)為轉子運動微分方程，左側的項整理后可以分為5類，用5個中括號加以區分，分別代表轉子結構參數，圓盤的質量、轉動慣量以及偏心距的影響，葉片橫向振動的影響，軸承支承剛度及阻尼的影響，碰摩力的影響。式(21)為葉片的運動微分方程，左側的項整理后可以分為4類，用4個中括號加以區分，分別代表葉片結構參數，轉子旋轉效應的影響，碰摩力的影響以及轉子橫向位移的影響。可以看出，軸及葉片的運動是高度耦合的。

式(20)和式(21)都是偏微分方程，直接求解比較困難，可利用Galerkin方法分離變量。為研究主共振頻率范圍內轉子的響應和失穩狀態，需要求解其主共振狀態下的振型。

葉片可視作懸臂梁，其振型為

ψi(V*)=(sin(βiLV*)-sinh(βiLV*))-

β1L=1.875 1,β2L=4.964 1… 0≤V*≤1

(22)

將軸看做兩端彈性支撐的Timoshenko梁，如圖2所示，l1,l2為圓盤與兩端的距離，截面關于Y軸和Z軸的轉角用θy和θz表示，令θ=θy+jθz。盤片系統簡化為集中質量點，其質量為mD，轉動慣量為Jp，直徑轉動慣量為Jd，考慮到陀螺效應，軸的振型可以表達為[21]

φni(x)=c1isinαnix+c2icosαnix+
c3isinhβnix+c4icoshβnix

(23)

圖2 盤片軸的簡化模型Fig.2 Simplified model of shaft-disk-blade

式中：i=1, 2；xi∈[0,li]；c1i，c2i,c3i,c4i為待定常數;αni和βni需要滿足

(24)

其中,

(25)

式中：ρ,Is,κ,E分別為軸的密度、截面慣性矩、截面剪切系數和楊氏模量；ω為轉子-葉片系統的固有頻率，主共振條件下，可認為Ω=ω。定義Θi(xi)為θi的振型函數，則有

(26)

根據兩端彈簧支承的邊界條件，可得

(27)

在圓盤接界處的邊界條件可表達為

(28)

將式(26)～式(28)代入式(23),即可求得軸在彈支邊界下的臨界轉速及其各階主共振的振型。軸的各階整體振型為

(29)

某軸長度為0.595 2 m、截面面積為3.14×10-4m2、截面慣性矩為7.85×10-9m4、楊氏模量為2×1011Pa、密度為7 800 kg/m3,xd=0.407 7 m,mD=0.735 kg,Jp=6.25×10-4kg·m2,不同邊界條件下前兩階臨界轉速，如表1所示，支承剛度k=2×106N/m時的歸一化振型，如圖3所示，為了驗證解析方法的準確性，采用Nastran建模并進行正進動和反進動固有頻率對比。從表1中可知，當支承剛度非常大時，彈支軸與簡支軸的固有頻率比較接近，隨著支承剛度降低，固有頻率降低，在離心力作用下，兩端支承處的相對位移較大。在高階振動中，盤所在的一側振幅會相對較小。

用Galerkin方法分離變量，可設定p*(V*,t*)=ψi(V*)P(t*),z*(x*,t*)=φn(x*)Z(t*)，代入式(20)和式(21)，并從0～1積分，可消去變量x*和V*，從而轉化為常微分方程，采用龍格庫塔方法對若干固定的轉速值進行數值求解，從而驗證本文多尺度攝動解的準確性。

圖3 旋轉軸的彈支振型 (k=2×106 N/m) Fig. 3 Mode shapes of the rotating shaft (k=2×106 N/m)

邊界條件支撐剛度/(N·m-1)一階臨界轉速/(Rad·s-1)解析法FEM正進動反進動正進動反進動二階臨界轉速/(Rad·s-1)解析法FEM正進動反進動正進動反進動彈支2×106490.29-487.63491.91-488.521 835.10-1 796.221 836.28-1 795.45彈支2×108523.65-518.95523.34-518.482 265.93-2 190.572 261.72-2 189.84簡支526.83-521.77526.01-521.432 315.76-2 198.332 312.92-2 194.23

2 多尺度求解

采用多尺度方法進行求解，p*(V*,t) 和z*(x*,t) 可以展開為

p*(V*,t)=εp1(V*,T0,T2)+ε3p3(V*,T0,T2)+…

z*(x*,t)=εz1(x*,T0,T2)+ε3z3(x*,T0,T2)+…

(30)

式中：T0=t,T2=ε2t。為了表征共振頻率與正進動頻率的接近程度，引入無量綱調諧參數σ=O(1),使Ω*=βf+ε2σ，其中βf為無量綱正進動頻率。將式(30)代入式(20)和式(21)，使兩端ε和ε3項對應系數相等，通過ε項可確定解的表達形式，代入ε3項可得方程

(31)

(32)

式中：Hz,f(x*,T2),Hz,b(x*,T2) ,Hp,f(V*,T2)和Hp,b(V*,T2) 為久期項，N.S.T為非久期項；βb為無量綱反進動頻率；D0=?/?T0。設定φf(x*),φb(x*),ψf(V*),ψb(V*)分別為軸和葉片的正進動與反進動振型，根據之前的討論，可令ψf(V*)=ψb(V*)=ψi(V*),φf(x*)=φb(x*) =φn(x*), 如果方程存在非零解，則必有

(33)

其中，

式(33)左側Jacobian矩陣的特征值可以給出系統的穩定性信息：若所有特征值的實部都為非正值，系統穩定，反之，若存在任意一個實部大于0的特征值，系統不穩定。

對于其他碰摩形式，只要碰摩力在每個旋轉周期內對時間可積就可以應用本文的方法。非周期的碰摩力進行如下等效即可

(34)

式中：Fni(t)為作用于第i個葉片上的任意形式的碰摩力。

3 數值算例

圖4 轉子-葉片系統示意圖Fig.4 Schematic of rotor-blade system

各工況中均采用龍格-庫塔數值積分方法驗證多尺度攝動解的準確性。圖5為不同碰摩力下轉子-葉片系統的響應曲線，可以看出系統為硬式非線性，在無碰摩力或碰摩力較小的條件下，系統可以維持穩定。隨著碰摩力的增加，非線性增強，共振峰值提高，最大碰摩力的值為0.225時，響應發生分叉，系統存在兩個穩定解和一個不穩定的解。

圖5 不同法向碰摩力下的響應幅值Fig.5 Response amplitude under different normal rubbing force

圖6為不同支承剛度下系統的響應幅值，支承剛度較低時，非線性較弱。當-0.2<σ<0.533時，低支承剛度下，系統的響應幅值較大。但隨著支承剛度的增加，系統的跳躍頻率升高，不穩定區域和共振峰值也隨之增加。

圖6 不同支承剛度下的響應幅值Fig.6 Response amplitude under different stiffness

圖7為不同摩擦因數下系統的響應幅值，在AB區域，未達到二者跳躍頻率，提高摩擦因數，系統的響應只是輕微地增加，幾乎保持不變。然而，摩擦因數的增加會提高跳躍頻率，通過額外的BC區域，系統的共振峰值得到提高。

圖7 不同摩擦因數下的響應幅值Fig.7 Response amplitude under different friction coefficient

圖8為不同圓盤偏心量下的響應幅值，從圖中可以看出，圓盤偏心量的增加，與碰摩力增加的情形類似，系統的非線性增強，跳躍頻率以及共振峰值也隨之增加。

圖8 不同圓盤偏心量下的響應幅值Fig.8 Amplitude-frequency responses under different eccentricities

圖9為不同支承剛度和轉軸阻尼系數下，系統的穩態響應幅值。從圖9中可以看出，系統若要達到完全穩定，支承剛度越大，則需要更大的阻尼。相同阻尼下，在未達到跳躍頻率時，支承剛度越小，系統響應幅值越大。隨著調諧參數的增加，共振分支的響應幅值增加。另外值得注意的是，無論支承剛度的數值大小，過小的阻尼都會使系統發生失穩。

圖9 轉軸阻尼系數對穩態響應的影響Fig.9 Steady state responses under different damping

圖10 軸承阻尼系數對穩態響應的影響Fig.10 Steady state responses under different bearing’s damping

圖11 碰摩力對系統穩態響應的影響Fig.11 Steady state responses under different rubbing force amplitudes μ=0.1,k*=26 860)

圖12 不同支承剛度下碰摩力對穩態響應的影響Fig.12 Steady state responses under different support stiffness

圖13描述了不同支承剛度和轉軸阻尼條件下第一和第二分叉點的變化軌跡。隨著阻尼的增加，兩類分叉點的頻率都在降低，而第二分叉點的速度更快，最終兩條軌線交于一點。當實際阻尼大于該值時，將確保系統不發生分叉，達到完全穩定。從圖中可知，當阻尼較小時，兩個系統都會發生失穩，阻尼越小，失穩區間越大。當支承剛度較高時，需要更大的阻尼以使系統達到完全穩定。

圖13 第一和第二分叉點的軌跡Fig.13 Loci of the first and second bifurcatin points

4 結 論

本文基于多尺度攝動理論，研究了轉子-葉片系統在主共振狀態下的響應，得出的結論如下：

(1) 轉子-葉片系統為硬式非線性，隨著碰摩力和圓盤偏心量的增加，非線性增強，系統穩定性降低，第一和第二分叉點頻率升高，跳躍頻率和共振峰值升高。

(2) 隨著支承剛度降低，系統非線性減弱，分叉點頻率降低，共振峰值和跳躍頻率隨之降低。支承剛度較高時，需要更大的結構阻尼達到穩定狀態。在達到二者跳躍頻率之前，支承剛度較低時響應幅值較大。

(3) 隨著摩擦因數的增加，系統的跳躍頻率和共振峰值提高。

(4) 較小的結構阻尼會導致系統產生失穩，支承剛度較低時，軸承阻尼對系統振動幅值的影響較大。