旋轉葉片剛柔耦合系統動力學分析

韓廣才，吳艷紅，王寅超，劉志強

(1.哈爾濱工程大學 航天與建筑工程學院，黑龍江 哈爾濱 150001;2.哈爾濱工程大學 機電工程學院，黑龍江 哈爾濱150001)

隨著汽輪機技術的發展，以及柔性多體系統動力學研究的深入，有關旋轉葉片縱向、橫向振動問題引起了眾多學者的關注.將葉片視為彈性梁的研究可見諸于許多文獻.Eringen等人［1］研究了無旋轉運動彈性梁的非線性振動問題，Woodall等人［2］只考慮了梁的橫向變形，推導出簡化了的彈性梁動力學運動微分方程.Anderson等人［3］推導了旋轉柔性梁的縱向和橫向的非線性微分方程，Hodges等人［4］對于Anderson的方程進行了修正.Ansari等人［5］分析了安裝在轉盤上的預扭非均質非對稱截面懸臂梁的非線性振動問題，表明在葉片穩態運動情況下，一個微小的擾動就會引起葉片的運動發生突變，且這種突變取決于擾動方向.2007年韓國Hong Seok Lim等人分析了集中質量沖擊作用下旋轉梁的振動問題，基于Kane方法推導了系統的動力學方程［6-7］.2007年臺灣學者研究了盤的柔性對軸盤葉系統耦合振動的影響，文中討論了軸的扭轉、盤的橫向振動及葉片的彎曲振動相互間的耦合關系［8］.

多年來準確預測經歷大范圍剛體運動和彈性變形的柔性葉片的動力學行為，始終是旋轉葉片動力學問題研究的主要課題，基于線性理論的傳統方法由于無法計及動力剛化效應，導致在許多實際應用中得到錯誤的結果.本文基于Hamilton原理建立了具有動力剛化效應的剛柔耦合旋轉葉片動力學微分方程，考慮了變截面、預扭角葉片模型，計及了截面離心慣性力、截面轉動和截面剪力等因素的影響，最后采用數值分析的方法通過算例分析了彈性體葉片端部振動位移響應.

1 動力學模型

如圖1所示，將葉片視為一均質旋轉鐵摩辛柯梁，安裝在半徑為e的繞空間固定軸旋轉的剛性圓盤上，考慮長度為L的雙邊楔形預扭葉片，設ρ為葉片單位體積的質量、E為楊氏模量，且寬度、厚度和截面扭轉角沿葉片長度方向線性變化.葉片左端O點固定于半徑為O1的剛性旋轉圓盤上，Oz為剛性圓盤中心，軸為葉片未變形時中性軸.圖1中，O1x1y1z1為慣性參考系，O1x'y'z'為固結于圓盤上的動參考系，繞慣性參考系O1y1軸以角速度Ω旋轉，Oxyz為動參考系，固結在未變形的葉片上，隨同葉片繞O1y1軸以角速度Ω旋轉，作用于圓盤旋轉平面內的力矩為τ，JD為剛性圓盤繞O1y1軸的轉動慣量.P為葉片未變形時中性軸上的一點，由于旋轉運動而發生了軸向變形及橫向變形，對于細長葉片橫向變形比軸向變形大得多，故只考慮橫向變形，以z表示葉片未變形時P點位置，葉片變形后P至P'點，P點在動參考系Oxyz內的位移分別用w(z)和v(z)表示，w(z)為xz平面內的位移，v(z)是yz平面內的位移.變形后P點在O1x'y'z'動參考系內的位置可用rP表示為

式中:i、j、k分別為動參考系O1x'y'z'沿坐標軸方向上的單位矢量，將式(1)在慣性參考系O1x1y1z1內求導，并注意到圓盤轉動的角速度矢量為Ω=θ·j(θ為剛性圓盤轉動角位移)，可得

則

圖1 旋轉柔性葉片剛－柔耦合動力學模型Fig.1 Dynamic model of disc and blade rigid-flexible coupling system

2 動力學方程的建立

本文采用Hamilton原理建立旋轉葉片系統的動力學方程，Hamilton原理的基本形式如下:

式中:T為系統的動能，U為系統的勢能，δW為外力所作的虛功.考慮P點在動參考系Oxyz內的位移w(z)和v(z)是由兩部分組成，即w=wb+ws，v= vb+vs，下腳標b表示由于彎曲引起的橫向位移，下腳標s表示由于剪切引起的橫向位移.旋轉葉片系統動能由兩部分組成，剛性圓盤的動能和柔性葉片的動能.在計算柔性葉片動能時，除考慮葉片微元隨質心平動動能外，還計及葉片微元繞其截面形心慣性主軸轉動動能，則系統動能可寫為

由彈性力學基本理論可知由彎曲產生的應變能Ub及剪切應變能Us分別為

考慮葉片變形后微元體的向徑為

則作用于z截面處軸向離心慣性力可寫為

作用于dz微元體上單位長度的離心慣性力x軸方向上的分量可寫為

則由于xz平面內的彎曲而產生的離心力勢能為

由于yz平面內的彎曲而產生的離心力勢能為

則旋轉葉片系統勢能為

分別計算動能T、勢能U和外力的功W對變形廣義坐標w、v和剛性圓盤轉動廣義坐標θ的變分，并代入式(2)中，得系統的動力學方程為

3 系統離散動力學方程

采用有限元法給出系統離散動力學方程，將葉片視為旋轉柔性梁，分成n個單元，第i個單元的長度設為li，如圖2所示，設單元的節點位移為

圖2 旋轉葉片單元節點位移Fig.2 Node displacements of a rotating flexible blade element

u1、u2、u3、u4為xz平面內彎曲引起的節點橫向位移及轉角位移，u5、u6為xz平面內剪切引起的節點橫向位移.u7、u8、u9、u10為yz平面內彎曲引起的節點橫向位移及轉角位移，u11、u12為yz平面內剪切引起的節點橫向位移，令

式中:U=Nu.

考慮雙側楔形帶預扭角變截面旋轉柔性葉片，葉片寬度、厚度和預扭角度沿葉片長度方向(軸方向)線性變化.如圖3所示，葉片根部和端部截面寬度分別為b1、b2，厚度分別為h1、h2，預扭角度分別為θ1、θ2，注意到由于單元截面及扭轉角沿z軸方向線性變化，故截面面積A及截面慣性矩Ixx、Ixy、Iyy是關于z的函數.將葉片單元離散，設單元局部坐標為，如圖4所示，則第i個單元中性軸上P點未變形時在Oxyz坐標系內的位置可表示為

圖3 帶有預扭角變截面旋轉葉片模型Fig.3 Sketch of a rotating blade with pre-twisted angle and variable cross-section

圖4 旋轉柔性葉片第i個單元局部坐標Fig.4 Local coordinate system of a rotating flexible blade element

令A=［1 0 0 0］，B=［0 1 0 0］，C=［0 0 1 0］，D=［0 0 0 1］.

葉片單元動能為

式中:

葉片單元勢能為

式中:

外力的虛功為

第i個單元的Lagrange函數為ψi=T－U并代入Lagrange方程得葉片第i個單元的動力學方程為

4 離散動力學方程的數值仿真

4.1 考慮預扭角變化葉片端部位移響應

葉片的參數取為L=0.6 m，b1=0.004 m，b2= 0.004 m，h1=0.15 m，h2=0.15 m，θ1=0°，ρ= 7 800 kg/m3，e=0.01 m，μ=0.3，E=2.07× 1011N/m2.

旋轉角運動的規律取為

取n=500 r/min，T=150 s.

1)分別取θ2=0°、15°、30°、45°時，葉片端部位移的響應見圖5.

由圖5(a)和(c)可知，預扭角的變化對于xz平面內橫向位移和轉角位移沒有影響.由圖5(b)和(d)可知當無預扭時，yz平面內不產生橫向位移和轉角位移，當有預扭角時，盡管葉片只在xz平面內轉動，也會在yz平面內產生橫向位移及轉角位移，且在15°～45°范圍內隨著預扭角增大，yz平面內橫向位移及轉角位移也隨之增大.

圖5 葉片端部橫向位移及轉角位移Fig.5 Flexual and slope displacement of blade tip

2)分別取θ2=45°、60°、75°、90°時，葉片端部位移的響應見圖6.

圖6 葉片端部橫向位移及轉角位移Fig.6 Flexual and slope displacement of blade tip

由圖6(a)和(c)可知，預扭角的變化對于xz平面內的橫向位移和轉角位移沒有影響.由圖6(b)和(d)可知，在45°～90°范圍內隨著預扭角的增大，yz平面內的橫向位移及轉角位移隨之減小.仿真結果表明，位移響應是穩定的并不隨時間發散.可見，預扭角的變化對yz平面內的橫向位移及轉角位移有較大影響.由于yz平面內的橫向振動會使葉片系統動靜部分發生碰擦，導致葉片通道堵塞，產生嚴重的摩擦聲，引發葉片斷裂事故.

4.2 考慮截面變化葉片端部位移響應

葉片的參數取為L=0.6 m，b1=0.004 m，h1= 0.15 m，θ1=0°，θ2=45°.旋轉角運動的規律取式(5)，取n=500 r/min，T=150 s.

考慮截面變化時葉片端部位移響應，見圖7.

圖7 葉片端部橫向位移及轉角位移Fig.7 Flexual and slope displacement of blade tip

由圖7(a)～(d)可以看出，葉片截面寬度和高度的變化對xz平面內的橫向位移和轉角位移影響較小，而對于yz平面內的橫向位移和轉角位移有較大影響.

5 結束語

基于Hamilton變分原理推導了帶有預扭角、變截面，考慮剪切變形、截面轉角及離心剛化效應的大范圍旋轉柔性葉片剛柔耦合系統非線性偏微分積分連續動力學方程組.基于有限元法建立了旋轉葉片系統非線性、時變及完全耦合的有限維離散動力學方程，基于數值仿真結果分析了截面參數對葉片橫向位移的影響，為變截面葉片的設計提供了理論依據，同時也為后續葉片位移響應的穩定性分析及幅頻特性研究提供了理論基礎及分析數據.

［1］ERIGEN A C.On the non-linear vibration of elastic bars［J］.Q Appled Mathematics，1952，9:361-369.

［2］WOODALL S R.On the large amplitude oscillations of a thin elastic beam［J］.Journal of Non-Linear Mechanics，1966，1:217-238.

［3］ANDERSON G L.On the extensional and flexual vibrations of rotating bars［J］.Journal of Non-Linear Mechanics，1975，10:223-236.

［4］HODGES D H.On the extensional vibrations of rotating bars［J］.Journal of Non-Linear Mechanics，1977，12:293-296.

［5］ANSARI K A.Non-linear vibrations of a rotating pretwisted blade［J］.Computers＆Structures，1975，5:101-108.

［6］LIM Hongseok.Impact analysis of a rotating beam due to particle mass collision［J］.Journal of Sound and vibration，2007，308:794-804.

［7］YOO Honghee.Model analysis and shape optimization of rotating cantilever beams［J］.Journal of Sound and vibration，2006，290:223-241.

［8］YANG Chiahao.The influence of discs flexibility on coupling vibration of shaft-dics-blades system［J］.Journal of Sound and Vibration，2007，301:1-17.

［9］PARK Junghun，PARK Hyunijong，JEONG Seokyong，et al.Linear vibration analysis of rotating wind-turbine blade［J］.Current Applied Physics，2010，10(2):S332-S334.