葉輪轉(zhuǎn)子碰摩的非線性動力學(xué)響應(yīng)

劉 昕， 張華彪， 孫小磊， 占傳林， 趙慶軍,3

(1．中國科學(xué)院 工程熱物理研究所，北京 100190；2．中國科學(xué)院大學(xué)，北京 100190；3．中國科學(xué)院 輕型動力重點(diǎn)實驗室，北京 100190)

?

葉輪轉(zhuǎn)子碰摩的非線性動力學(xué)響應(yīng)

劉 昕1,2， 張華彪1， 孫小磊1， 占傳林1， 趙慶軍1,3

(1．中國科學(xué)院 工程熱物理研究所，北京 100190；2．中國科學(xué)院大學(xué)，北京 100190；3．中國科學(xué)院 輕型動力重點(diǎn)實驗室，北京 100190)

葉輪轉(zhuǎn)子的碰摩是旋轉(zhuǎn)機(jī)械的常見故障，嚴(yán)重的碰摩可能導(dǎo)致轉(zhuǎn)子失穩(wěn)、葉片斷裂等后果。針對葉輪轉(zhuǎn)子-機(jī)匣的碰摩進(jìn)行研究，將葉片簡化為一端固支的懸臂梁，考慮其在旋轉(zhuǎn)狀態(tài)下的離心剛化效應(yīng)，機(jī)匣的局部變形采用Hertz模型表征，推導(dǎo)了葉片-機(jī)匣接觸力的表達(dá)式，分析了葉片長度、截面慣性矩、葉片截面積和轉(zhuǎn)速對接觸力的影響，結(jié)果表明當(dāng)葉片長度越小、截面慣性矩越大、轉(zhuǎn)速越高接觸力越大?；谒玫慕佑|關(guān)系，給出了葉輪轉(zhuǎn)子碰摩力模型，對葉輪轉(zhuǎn)子碰摩的穩(wěn)態(tài)響應(yīng)進(jìn)行數(shù)值模擬，結(jié)果表明系統(tǒng)在碰摩時的響應(yīng)頻率主要包括轉(zhuǎn)速激勵頻率、反向進(jìn)動頻率及其線性組合，葉片的擾動將導(dǎo)致主要響應(yīng)頻率附近出現(xiàn)邊頻，邊頻在不同轉(zhuǎn)速情況下存在多種不同形式。當(dāng)反向進(jìn)動頻率相對于轉(zhuǎn)速近似線性變化時，邊頻和轉(zhuǎn)速與葉片數(shù)目線性相關(guān)。

葉輪轉(zhuǎn)子；碰摩；柔性機(jī)匣；邊頻

碰摩是葉輪機(jī)械常見故障之一，碰摩的發(fā)生可能導(dǎo)致機(jī)匣變形、軸承磨損、葉片斷裂甚至轉(zhuǎn)子失穩(wěn)，最終造成嚴(yán)重的后果。

國內(nèi)外學(xué)者針對轉(zhuǎn)子系統(tǒng)的碰摩問題開展了廣泛的研究。ZHANG等[1]分析了全周碰摩的精確解和解的穩(wěn)定域。JIANG等[2]研究了同步全周/部分碰摩運(yùn)動的穩(wěn)定性。ZHANG等[3]采用了C-L方法對同步全周碰摩運(yùn)動的分岔響應(yīng)進(jìn)行了分析。CHU等[4]利用Fourier分析和Floquet理論對跨中單轉(zhuǎn)子的全局性分岔特性和穩(wěn)定性進(jìn)行研究，并通過數(shù)值模擬驗證了通往混沌的三種路徑。上述文獻(xiàn)在研究中使用了圓盤轉(zhuǎn)子模型，雖然便于分析，卻不能很好的描述葉片和機(jī)匣的碰摩現(xiàn)象。PADOVAN等[5]將葉片簡化為懸臂梁推導(dǎo)了葉片-剛性機(jī)匣碰摩力公式。PADOVAN等[6-8]針對不同情況的碰摩侵入量進(jìn)行分析。JIANG等[9]和MA等[10]分別考慮轉(zhuǎn)速和離心強(qiáng)化的影響，提出了不同的葉片機(jī)匣碰摩力模型，他們的研究更加側(cè)重于葉片本身的響應(yīng)。目前，針對工程中常見的葉輪轉(zhuǎn)子-機(jī)匣碰摩問題的研究并不多見。

本文將葉片簡化為懸臂歐拉梁，并考慮其在離心力作用下的剛度變化，利用Hertz模型描述葉片和機(jī)匣碰摩時的接觸關(guān)系，構(gòu)造了一種新的葉片-柔性機(jī)匣碰摩力的模型，并利用該模型對葉輪轉(zhuǎn)子-柔性機(jī)匣碰摩的非線性動力學(xué)響應(yīng)進(jìn)行研究。

1 葉輪轉(zhuǎn)子-機(jī)匣接觸力模型

將葉片假設(shè)為一端固定于輪盤的等截面懸臂梁，圖1中給出了葉片的受力情況。

圖1 葉片簡化模型Fig.1 Schematic diagram of rotating blade

在豎直方向存在離心力產(chǎn)生的分布載荷qω和法向接觸力FN，水平方向存在摩擦力FT，用庫倫模型描述，有

FT=μFN

(1)

式中：μ為摩擦因數(shù)。對分布載荷qω積分可得作用在距離葉根x位置處截面上的離心力，

Fω(x)=

(2)

式中：ρ,ω,A,Rr,lb分別為葉片密度，工作轉(zhuǎn)速，葉片橫截面積，輪盤半徑和葉片展向長度。根據(jù)歐拉梁方程可得葉片的撓度方程，有

(3)

把離心力表達(dá)式(2)代入，展開，有

y″+(c+bx+ax2)y=hx

(4)

式中:

(5)

式(4)是一個非齊次非線性的微分方程，直接求解比較困難，這里采用冪級數(shù)法進(jìn)行求解，設(shè)解為

(6)

考慮葉片葉根處撓度和轉(zhuǎn)角為0，有y(0)=0,y′(0)=0。將式(6)代入式(4)，忽略高次項，令其余項系數(shù)為0，可解得

(7)

式中:n為解的階次。下面對所得冪級數(shù)解進(jìn)行分析，如未特殊聲明，本文中的計算選取的參數(shù)由式(8)給出

ρ=8.1g/cm3,A=64mm2,Rr=73mm,μ=0.13

lb=45mm,I=210mm4,E=2.1×1011N/m2

(8)

圖2給出了不同階次的冪級數(shù)解與通過數(shù)值積分所得的撓度曲線的對比，其中橫坐標(biāo)為沿展向長度x，縱坐標(biāo)為撓度y，可以看到對于本文所研究的葉片，5階冪級數(shù)解和數(shù)值積分的誤差已經(jīng)小于1%，因此本文下面的計算取n=5。

圖2 葉片撓度曲線Fig.2 Deflection curve of the blade

由式(6)和(7)，利用近似關(guān)系計算葉尖的徑向變形量為

(9)

同時，考察轉(zhuǎn)子運(yùn)動過程中的轉(zhuǎn)子位移與機(jī)匣變形量之間的關(guān)系，如圖3。

圖3 葉片-機(jī)匣碰摩幾何位置關(guān)系Fig.3 Geometric model of contact between blades and casing

2δ(lb+Rr)cos(π-T))1/2

(10)

(11)

葉片和機(jī)匣的接觸采用Hertz模型，接觸力為

(12)

式中:α為與機(jī)匣和葉片材料相關(guān)的系數(shù)，反映了機(jī)匣的局部變形程度，其表達(dá)式為：

(13)

式中:E1,E2,η1,η2分別為機(jī)匣和葉片的彈性模量和泊松比，本文中分別都取2.1×1011N/m和0.2。

((lb+Rr)2+δ2-2δ(lb+Rr)cos(π-T))1/2-

Rc=α(FN)2/3+δb(FN)

(14)

式中:δb(FN)是關(guān)于FN的高次多項式。注意到上述表達(dá)式中左邊項實際上就是葉片和機(jī)匣變形量之和，記做δtotal，而右邊為一個關(guān)于FN代數(shù)多項式，使用葉片與機(jī)匣間的靜態(tài)間隙c對δtotal進(jìn)行無量綱化，有

(15)

(16)

2 葉片參數(shù)對接觸力的影響

根據(jù)上面求得的接觸力公式，給定不同葉片參數(shù)代入計算，可以得到葉片參數(shù)和轉(zhuǎn)速對接觸力的影響。

圖4～圖7分別給出了葉片截面積、葉輪轉(zhuǎn)速、葉片截面慣性矩和葉片長度對接觸力的影響?？梢钥吹讲煌~片截面積和不同轉(zhuǎn)速條件下相同變形量對應(yīng)的接觸力差異不大，而截面慣性矩和葉片長度對接觸力的影響較大。接觸力隨著葉片截面積的增大而略有增大。當(dāng)轉(zhuǎn)速較高時，離心強(qiáng)化效應(yīng)更強(qiáng)，葉片-機(jī)匣接觸的等效剛度更高，接觸力更大。對于一定的總體變形量，葉片長度越大，截面慣性矩越小，葉片就越容易彎曲變形，對應(yīng)的機(jī)匣的變形量就越小，而根據(jù)Hertz接觸模型，較小的變形量對應(yīng)較小的接觸力，因此較大的葉片長度或較小的彎曲剛度對應(yīng)較小的接觸力。

圖4 葉片截面積對接觸力的影響Fig.4 Contact force under different cross sections of blade

圖5 轉(zhuǎn)速對接觸力的影響Fig.5 Contact force under different angular speeds

圖6 葉片截面慣性矩對接觸力的影響Fig.6 Contact force under different moments of inertia

圖7 葉片長度對接觸力的影響Fig.7 Contact force under different lengths of blade

3 葉輪轉(zhuǎn)子運(yùn)動方程

考慮跨中葉輪轉(zhuǎn)子的碰摩，葉輪轉(zhuǎn)子-機(jī)匣系統(tǒng)的動力學(xué)方程可寫作

(17)

式中:轉(zhuǎn)子質(zhì)量為m，機(jī)匣質(zhì)量為M，Ce,Ke,eu分別為阻尼系數(shù)、支承剛度和不平衡量，角標(biāo)r和c分別代表轉(zhuǎn)子和機(jī)匣。Fx和Fy分別為碰摩力在x和y方向分量，當(dāng)無碰摩時為0，發(fā)生碰摩時對單個葉片有

(18)

Xc=x,Yc=y,τ=ωnt

(19)

(20)

為避免對全周葉片逐個檢測判斷是否發(fā)生碰摩，葉片/機(jī)匣碰摩事件的檢測采用PADOVA給出的方法[5]。圖8給出了發(fā)生碰摩時葉輪轉(zhuǎn)子和機(jī)匣的幾何關(guān)系。根據(jù)轉(zhuǎn)子進(jìn)動角φ與自轉(zhuǎn)角sτ，確定可能發(fā)生接觸的中間位置的葉片的角位置T，從而可以得到中間位置葉片與機(jī)匣的距離Dc以及相鄰兩葉片對應(yīng)的距離Dr和Dl，根據(jù)Dc、Dr和Dl即可判斷葉片與機(jī)匣是否發(fā)生碰摩。具體方法為：① 若Dc>0，說明沒有發(fā)生接觸;② 若Dc<0，說明中心位置的葉片與機(jī)匣發(fā)生了接觸;③ 再分別將兩側(cè)葉片與機(jī)匣的距離Dr和Dl視為Dc，重復(fù)步驟①，直到兩側(cè)葉片與機(jī)匣距離都大于0為止，從而得到所有與機(jī)匣接觸葉片的侵入量。按照前面所述的方法求出每個葉片的碰摩力FN，并按對應(yīng)角度位置求出總的矢量和，從而得到碰摩力Fx和Fy。將碰摩力代入方程，利用變步長數(shù)值積分方法即可對式(17)求解。

圖8 轉(zhuǎn)子機(jī)匣幾何位置關(guān)系Fig.8 Geometric model of rotor and casing

4 葉輪轉(zhuǎn)子-機(jī)匣碰摩響應(yīng)

在本文后面的計算中，如未特殊聲明，取阻尼系數(shù)Cre=2 000 Ns/m,Cce=6 000 Ns/m,，不平衡量eu=7×10-6m，葉尖與機(jī)匣的靜態(tài)間隙值為c=6×10-5m，支撐剛度Kre=2.5×108N/m,Kce=5×108N/m，轉(zhuǎn)子與機(jī)匣質(zhì)量取m=M=10 kg。

作為對照，本文計算了圓盤轉(zhuǎn)子與機(jī)匣的碰摩響應(yīng)，接觸力分別采用線性當(dāng)量剛度模型和Hertz模型進(jìn)行對比。由于本文中涉及到的葉片/機(jī)匣侵入量最大不超過間隙值c的15%，在此范圍內(nèi)對Hertz模型給出的侵入量-接觸力曲線進(jìn)行線性擬合，從而得到一個恒定剛度Kc近似代替Hertz模型。根據(jù)本文中的參數(shù)有Kc=1.5×109N/m，碰摩力表達(dá)式為

FN=Kcδtotal

(21)

式中:δtotal為葉片/機(jī)匣侵入量。圖9給出了兩種模型條件下接觸力的對比。

圖9 線性當(dāng)量剛度模型與Hertz模型的碰摩力關(guān)系Fig.9 Contact force under linear model and Hertz model

圖10和圖11分別給出了當(dāng)量剛度模型和Hertz模型條件下轉(zhuǎn)子發(fā)生碰摩時穩(wěn)態(tài)響應(yīng)的頻譜云圖，橫坐標(biāo)為無量綱頻率，s>0.95時轉(zhuǎn)子和機(jī)匣發(fā)生碰摩，位移響應(yīng)頻譜云圖上出現(xiàn)了非協(xié)調(diào)進(jìn)動的頻率成分。圖12給出了圓盤轉(zhuǎn)子碰摩的位移響應(yīng)頻譜和軸心軌跡。

圖10 當(dāng)量剛度模型圓盤轉(zhuǎn)子碰摩響應(yīng)頻譜云圖Fig.10 Spectrum contour of a rubbing 0-blade rotor with constant casing equivalent stiffness

圖11 Hertz模型圓盤轉(zhuǎn)子碰摩響應(yīng)頻譜云圖(無葉片)Fig.11 Spectrum contour of a rubbing 0-blade rotor with Hertz model

圖12 當(dāng)量剛度模型圓盤轉(zhuǎn)子的碰摩響應(yīng)(s=1.04)Fig.12 Rubbing response of the no-blade rotor with constant casing equivalent stiffness (s=1.04)

可以看到圖10中除了存在轉(zhuǎn)速激勵頻率fω外，還存在轉(zhuǎn)子反向進(jìn)動頻率fN及2fω-fN和2fN-fω等組合響應(yīng)頻率。將這些響應(yīng)記作fO，則fO可統(tǒng)一寫作如下形式

fO=n(fω-fN)+fN,

n=…,-2,-1,0,1,2,…

(22)

式中:n為整數(shù)。另外，在圖10中fO的兩側(cè)還存在一對幅值更小的等間隔分布的頻率fS。對比圖11中Hertz模型機(jī)匣的轉(zhuǎn)子系統(tǒng)的碰摩響應(yīng)，可以看出兩者的頻譜基本一致。

圖13給出了采用當(dāng)量剛度模型的葉輪轉(zhuǎn)子碰摩響應(yīng)的頻譜云圖，可以看到葉輪轉(zhuǎn)子的碰摩響應(yīng)頻譜和圓盤轉(zhuǎn)子的頻譜相比有較大區(qū)別。當(dāng)s<0.97時，頻譜和圓盤轉(zhuǎn)子類似只存在fO與fS，而當(dāng)s>0.97，fO與fS不再隨著轉(zhuǎn)速線性變化，且在這些頻率附近出現(xiàn)不同形式的邊頻fside。圖14給出了對應(yīng)的葉輪轉(zhuǎn)子碰摩的響應(yīng)頻譜和軸心軌跡。

圖13 當(dāng)量剛度模型葉輪轉(zhuǎn)子碰摩響應(yīng)頻譜云圖(葉片數(shù)目17)Fig.13 Spectrum contour of a rubbing 17-blade rotor with constant casing equivalent stiffness

圖14 當(dāng)量剛度模型葉輪轉(zhuǎn)子的碰摩響應(yīng)(葉片數(shù)目17，s=1.1)Fig.14 Rubbing response of the blade-rotor with constant casing equivalent stiffness (17 blades, s=1.1)

圖15 當(dāng)量剛度模型葉輪轉(zhuǎn)子碰摩響應(yīng)頻譜云圖局部，葉片數(shù)目17Fig.15 Details of the spectrum contour of a rubbing 17-blade rotor with constant casing equivalent stiffness

(23)

式中:n為整數(shù),Nblade為葉片數(shù)目。圖16給出了Hertz模型機(jī)匣的碰摩葉輪轉(zhuǎn)子的頻譜云圖，注意到不同形式的邊頻依然存在，說明邊頻并不依賴與特定的碰摩力模型，而是葉輪碰摩轉(zhuǎn)子的獨(dú)特性質(zhì)。

圖16 Hertz模型葉輪轉(zhuǎn)子碰摩響應(yīng)頻譜云圖，葉片數(shù)目17Fig.16 Spectrum contour of a rubbing 17-blade rotor with Hertz model

5 結(jié) 論

本文提出了一種新的葉輪轉(zhuǎn)子葉尖-機(jī)匣碰摩力模型，此模型考慮了高速旋轉(zhuǎn)葉片的離心強(qiáng)化作用，使用Hertz模型表征機(jī)匣局部變形情況下葉片和機(jī)匣的接觸。對該模型的分析表明，在葉片對機(jī)匣的侵入量相同的情況下，葉片展向長度越短、截面慣性矩越大、工作轉(zhuǎn)速越高，單個葉片上承受載荷越大。

基于葉輪轉(zhuǎn)子葉尖-機(jī)匣碰摩力模型，建立了葉輪轉(zhuǎn)子-機(jī)匣碰摩的運(yùn)動方程并進(jìn)行了數(shù)值模擬，計算結(jié)果表明在葉片與機(jī)匣發(fā)生碰摩時，位移響應(yīng)頻譜上存在反向進(jìn)動頻率和轉(zhuǎn)速激勵頻率的線性組合，對于葉輪轉(zhuǎn)子還在主要響應(yīng)頻率附近還存在較小的邊頻響應(yīng)。邊頻在不同轉(zhuǎn)速情況下存在多種不同形式。當(dāng)反向進(jìn)動頻率相對于轉(zhuǎn)速近似線性變化時，邊頻與主要響應(yīng)頻率之差相對于轉(zhuǎn)速線性變化，在位移響應(yīng)頻譜云圖上的邊頻的斜率和主要響應(yīng)頻率的斜率之差與葉片數(shù)目之比為有理數(shù)，體現(xiàn)了葉輪轉(zhuǎn)子碰摩的獨(dú)特性質(zhì)。

[1] ZHANG G F, XU W N, XU B, et al. Analytical study of nonlinear synchronous full annular rub motion of flexible rotor-stator system and its dynamic stability[J]. Nonlinear Dynamics, 2009, 57(4): 579-592.

[2] JIANG J, ULBRICH H. Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients[J]. Nonlinear Dynamics, 2001, 24(3): 269-283.[3] ZHANG H, CHEN Y. Bifurcation analysis on full annular rub of a nonlinear rotor system[J]. Science China Technological Sciences, 2011, 54(8): 1977-1985.

[4] CHU F, ZHANG Z. Bifurcation and chaos in a rub-impact Jeffcott rotor system[J]. Journal of Sound and Vibration, 1998, 210(1): 1-18.

[5] PADOVAN J, CHOY F. Nonlinear dynamics of rotor/blade/casing rub interactions[J]. Journal of Turbomachinery, 1987, 109(4): 527-534.

[6] PADOVA C, BARTON J, DUNN M G, et al. Experimental results from controlled blade tip/shroud rubs at engine speed[J]. Journal of Turbomachinery, 2007, 129(4): 713-723.

[7] PADOVA C, BARTON J, DUNN M G, et al. Development of an experimental capability to produce controlled blade tip/shroud rubs at engine speed[J]. Journal of Turbomachinery, 2005, 127(4): 726-735.

[8] PADOVA C, DUNN M, BARTON J, et al. Controlled fan blade tip/shroud rubs at engine conditions[C]// ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition. Vancouver,Canada, 2011: 951-961.

[9] JIANG J, AHRENS J, ULBRICH H, et al. A contact model of a rotating rubbing blade[C]// Proceedings of the 5th International Conference on Rotor Dynamics, echnische Universitat Darmstadt. Darmstadt, Germany, 1998: 478-489.

[10] MA H, TAI X, HAN Q, et al. A revised model for rubbing between rotating blade and elastic casing[J]. Journal of Sound and Vibration, 2015,337:301-320.

[11] JIANG J. The analytical solution and the existence condition of dry friction backward whirl in rotor-to-stator contact systems[J]. Journal of Vibration and Acoustics, 2007, 129(2): 260-264.

Nonlinear dynamics on rubbing of a blade-rotor system

LIU Xin1, 2, ZHANG Huabiao1, SUN Xiaolei1, ZHAN Chuanlin1, ZHAO Qingjun1,3

(1. Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing100190, China;2. University of Chinese Academy of Sciences, Beijing100190, China;3. Key Laboratory of Light-duty Gas-turbine, Chinese Academy of Sciences, Beijing100190, China)

As a type of common failure of the rotating machines, rubbing can lead to catastrophic consequences such as instability of the rotor or blade damages. Rubbing between blades and casing was investigated in this paper. The blades were simplified as cantilever beams with centrifugal stiffening effect, while the blade/casing contact was modeled by Hertz contact force and an expression of blade/casing contact force has been derived. Based on the derived expression, influence of various parameters on the contact force was analyzed, such as the rotating speed, cross section of blade, length of blade, and bending stiffness of blade. The results show that blades with shorter length and larger moment of inertia or under higher rotating speed will lead to greater contact force on the casing. Meanwhile, a blade-rotor/casing rubbing model was derived with the given expression of contact force and numeric simulations with this model were also conducted. The results show that the responses on the spectrum consist of synchronous frequency, reverse precession frequencies, and their linear combinations. It was also observed that the blades could lead to side frequencies near the main frequencies. Different kinds of side frequencies were observed under different angular speeds. When the reverse precession frequencies are nearly linear with the angular speed, the side frequencies are also linear with the angular speed as well as the number of blades.

blade-rotor; rubbing, flexible casing; side frequencies

國家自然科學(xué)基金(11302223)；國家科技支撐計劃資助項目(2012BAA11B01)

2015-06-11 修改稿收到日期：2015-09-24

劉昕 男，碩士生，1991年生

張華彪 男，博士，助理研究員，1984年生

TH113.1

A

10.13465/j.cnki.jvs.2016.20.028