王鑫鑫， 葉天貴， 靳國永， 劉志剛

(哈爾濱工程大學 動力與能源工程學院，黑龍江 哈爾濱 150001)

葉片作為燃氣輪機的關鍵部件，其動力學性能直接影響著燃氣輪機運行的安全性和穩(wěn)定性。在實際運行中，由于葉片處于高速旋轉的工作狀態(tài)，不僅承受著強大的離心力載荷，而且需要抵抗高強度的氣流壓力，甚至渦輪葉片還會經受極大的熱載荷，工作環(huán)境十分苛刻。開展葉片動力學特性分析，掌握葉片的動力學性能具有重要理論價值和工程指導意義。

目前，工程中常用于求解旋轉葉片振動特性的方法有有限元法[1-3]、Galerkin法[4-5]、Ritz方法[6-7]等。牛宏偉等[8]通過有限元方法對航空發(fā)動機壓氣機轉子葉片進行模態(tài)分析，揭示了動頻曲線轉向趨勢和振型耦合規(guī)律。Sun等[9]建立了具有任意預置角的旋轉薄板動力學模型，結合Galerkin法進行求解，并對旋轉葉片的受迫響應進行了研究。李紅影等[10]用殼模型模擬凸肩葉片，運用Galerkin方法分析了凸肩葉片的分岔特性。Sinha等[11]從薄殼理論出發(fā)并結合Ritz方法研究了旋轉預扭葉片的振動特性。但是，旋轉葉片并不是一個獨立的個體，在工作過程中會受到氣動載荷的影響，尤其是軸流式旋轉機械，均勻靜葉柵尾流會使下一級動葉受到周期性的尾流激振力，是導致動葉片疲勞破壞的主要因素[12]。初世明[13]考慮尾流激振力的作用，建立了柔性帶冠旋轉梁模型并對其非線性動力學特征展開研究。王迪[14]建立了多載荷作用下的旋轉板模型，討論了旋轉葉片-機匣的碰摩問題。鄭麗娟[15]與羅驍?shù)萚16-17]采用數(shù)值方法分析了尾流激勵下的旋轉葉片氣動彈性行為。可以看出，許多學者針對各向同性材料葉片做了大量研究工作，已逐漸形成了較為系統(tǒng)的分析方法。

然而，隨著機械結構工程的迅速發(fā)展，尤其在航空航天領域，人們對葉片強度、剛度、耐熱、耐久性的要求越來越高，傳統(tǒng)材料已經不能滿足這些技術上要求。近年來，功能梯度材料作為一種兼具優(yōu)良力學性能和熱穩(wěn)定性能的新型材料，不僅在一定程度上彌補了傳統(tǒng)材料的不足之處，而且克服了多層復合材料層間應力集中、開膠以及斷裂的缺點，應用越來越廣泛。該材料在燃氣輪機葉片制造領域具有重要的發(fā)展前景和實用價值，吸引了大量學者們的關注。目前針對旋轉功能梯度葉片的研究主要集中在旋轉梁[18-21]與旋轉板[22-24]的自由振動，關于尾流激勵下具有任意預置角的旋轉功能梯度葉片動力學特性的研究較為少見，同時材料在空間上的連續(xù)梯度分布也對該類結構的動力學分析帶來了新的挑戰(zhàn)。

本文以多項正弦諧波疊加形式的尾流激振力作為氣動載荷，基于經典板理論建立了具有任意預置角的旋轉功能梯度葉片模型，通過罰函數(shù)將葉片的邊界條件轉化為邊界勢能的形式，以改進Fourier級數(shù)作為試函數(shù)表征葉片振動位移場。因此，本文利用Ritz方法對旋轉葉片動力學行為進行求解時邊界條件可以任意設置，且可根據(jù)實際工作選取合適的試函數(shù)以簡化運算，更有利于建立統(tǒng)一的求解方法。

1 模型與方法

1.1 模型的建立

如圖1所示，本文將具有任意預置角的旋轉功能梯度葉片簡化為一邊固支在剛性輪轂上的板模型。

圖1 旋轉功能梯度葉片模型Fig.1 The model of the functionally graded blade

其中，X-Y-Z為主體坐標系，x′-y′-z′為旋轉坐標系，x-y-z為局部坐標系，剛性輪轂的半徑為R，帶動葉片以Θ的角速度進行旋轉，葉片長l寬b厚h，在剛性輪轂上的預置角為θ。

假設葉片的材料組分按照Voigt模型沿厚度方向連續(xù)過渡，則楊氏模量和密度可表示為:

(1)

式中：下標c與m分別代表上層的陶瓷組分與下層的金屬組分，由于2種材料的泊松比μ相差無幾，因此將泊松比設為常數(shù)。

葉片振動位移場可表示為：

(2)

幾何方程為：

(3)

為簡便起見，簡化式(3)中的偏導數(shù)表達形式，例如εx將表示為εx=ux-zwxx。

物理方程為：

(4)

根據(jù)勢能的定義，葉片彎曲勢能Us表達式為：

(5)

罰函數(shù)方法引入的邊界彈簧勢能Up可表示為：

(6)

式中：k代表支撐彈簧剛度；K代表旋轉彈簧剛度；上標u、v、w代表約束形變位移的方向。

旋轉引起的離心力使得板結構內部產生特定的應力場，一般而言，薄板結構厚度方向尺寸遠小于其余2個方向尺寸。為了簡便起見，此處僅考慮x、y方向的離心應力分量，因此葉片旋轉所產生的附加離心力勢能Uc為[23]：

(7)

式中：χx與χy分別代表離心應力沿葉片x、y方向的分量：

圖2為葉片俯視圖與側視圖，其中α是在時間t內旋轉過的角度，A是葉片上一點。設A在旋轉坐標系內的位置為(x′,y′,z′)，在主體坐標系內的位置為(X,Y,Z)，在局部坐標系內的位置為(x,y,z)，通過A點可建立起3組坐標系之間的關系：

圖2 旋轉功能梯度葉片俯視圖與側視圖Fig.2 Top and side views of the functionally graded blade

(8)

同樣可以得到：

(9)

考慮到功能梯度葉片在旋轉過程中的位移變形U、V、W，對式(9)進行修正：

(10)

由式(8)、(10)可以得到主體坐標系X-Y-Z中旋轉梯度葉片上任意一點A的位置表達式：

(11)

由式(11)可得葉片上A點的速度表達式：

(12)

旋轉功能梯度葉片的動能可表示為：

(13)

將式(2)與式(12)代入式(13),可以得到動能T的具體表達式：

Θ2w2cos2θ-2zΘ2wywsinθcosθ-2zΘ2uwx+

(14)

對于軸流式旋轉機械來說，由于上游葉片尾緣和附面層的影響，引起導流葉柵出口氣流速度發(fā)生虧損，從而形成尾流，同時因為上游葉片與下游葉片間的相對轉動，使得上游尾流具有周期性激振效應[16]。如圖3所示。

圖3 尾流激振力示意Fig.3 The model of trail excitation

θ′為上游靜葉的預置角，假設該角度為定值30°；Fq為尾流激振力，可將其表示為多項正弦諧波疊加的形式[13]：

Fq=F0+F1sin(nΘt)+

F2sin(2nΘt)+…+Fksin(knΘt)

(15)

式中：Fq是以時間為自變量的周期函數(shù)；F0為一常值，與時間t無關；Fk為第k次諧波分量的幅值，在得知具體的Fq后，F(xiàn)k可以通過傅里葉級數(shù)展開得到[12]；Θ為葉片旋轉角速度；n為上級靜葉數(shù)目，決定尾流激振力的基頻。

由線性疊加原理可知，旋轉葉片在尾流激振力下總的響應等效于各激振力諧波單獨作用下響應的疊加。本文取第k次諧波進行分析，并只考慮垂直于葉片表面的尾流激振力。根據(jù)圖3中的幾何關系，可以得到垂直于葉片表面的第k次尾流激振力Fkτ的大小為：

(16)

則尾流激振力對葉片做功為：

(17)

1.2 位移試函數(shù)的選取

在本文方法中，構造適當?shù)奈灰圃嚭瘮?shù)是至關重要的。為克服傳統(tǒng)Fourier級數(shù)的弊端，本文引入改進Fourier級數(shù)來描述葉片的位移[25-26]，具體形式為：

式中：λlm=mπ/l，λbn=nπ/b；A、B、C是未知系數(shù)；M、N代表截斷級數(shù)。

在位移試函數(shù)中，雙重余弦項為位移主體項，ηi、ηj、ξi、ξj為輔助函數(shù)，其目的是為了保證容許函數(shù)在所有邊界點上連續(xù)可導，因此這些輔助函數(shù)的構造形式并不唯一。為了簡化后續(xù)的數(shù)學推導過程，構造輔助函數(shù)：

其中i=j=1, 2, 3, 4。

1.3 Ritz求解過程

旋轉功能梯度板的能量泛函可以表示為：

L=T-U+W0

(18)

式中：T為系統(tǒng)的總動能；U代表系統(tǒng)的總勢能；W0代表尾流激振力做功。

令能量泛函對未知系數(shù)求偏導為零，可以得到旋轉功能梯度葉片的振動微分方程：

(19)

式中：M為質量陣；C為阻尼陣；K為剛度陣；q為未知系數(shù)組成的列向量；F為廣義力向量。當外力功為零時，式(19)退化為求解自由振動的微分方程。

在處理式(19)的特征值問題時，可將正弦形式的激振力視為exp(ikΘt)形式簡諧力的虛部。由于尾流激振力是與葉片轉速Θ有關的，因此通過掃頻求解可以得到不同轉速下的多組未知系數(shù)q。為反映旋轉葉片整體的振動情況，引入均方根速度：

(20)

式中S為葉片的表面積。

計算結果表示為均方振速級Lv：

Lv=20lg(V/Vref)

(21)

其中Vref=10-9m/s。

本文采用復合楊氏模量的形式將結構阻尼包括在葉片模型中，即E^=E(z)(1+βi),其中β為結構阻尼系數(shù)，在本文中取β=0.01。

2 計算結果及分析

2.1 收斂性分析

由于本文采用基于改進Fourier級數(shù)的Ritz法分析葉片的振動特性，因此該方法的計算速度與精度和截斷級數(shù)M、N的選取有關。因而，通過收斂性分析來確定所需的截斷級數(shù)對后文的分析至關重要。本算例中：楊氏模量E=221 GPa，密度ρ=7 850 kg/m3，泊松比μ=0.3，葉片尺寸180 mm×85 mm×2 mm，邊界彈簧剛度、剛性輪轂半徑R、輪轂轉速Θ與預置角θ均設為0。表1給出了不同截斷級數(shù)下葉片的前8階固有頻率，此處沒有考慮剛體模態(tài)，且本文計算結果均保留5位有效數(shù)字。從表1中可以看出，隨著截斷級數(shù)的增加，各階固有頻率逐漸收斂。綜合考慮求解精度與計算速度，在后續(xù)的分析中令M=N=12。

表1 不同截斷級數(shù)下葉片前8階固有頻率Table 1 The first eight natural frequencies of the blade with respect to truncated numbers

此外，本文采用罰函數(shù)處理葉片的邊界條件，因此彈簧剛度的取值會直接影響求解結果。為準確模擬不同的邊界條件，定義剛度參數(shù)Γm，令k=K=10ΓmN/m，即假設支撐彈簧剛度k與旋轉彈簧剛度K同等變化，其他參數(shù)與表1算例相同。圖4展示了前3階固有頻率隨Γm的變化趨勢，從圖中可以直觀看出：在區(qū)域Ⅰ與區(qū)域Ⅱ內，各階固有頻率保持不變，分別對應自由邊界區(qū)間與固支邊界區(qū)間；在中間區(qū)域各階頻率會隨著Γm迅速增大，該區(qū)域對應彈性邊界區(qū)間。在后續(xù)的分析中，本文選擇Γm=0與Γm=12分別作為自由邊界值與固支邊界值。

圖4 邊界約束彈簧的收斂性Fig.4 Convergence of Γm for mechanical boundaries

2.2 正確性驗證

為驗證本文方法的正確性，表2給出了不同預置角下旋轉葉片的前5階固有頻率，同時列出了文獻[9]中的部分計算結果。在該算例中：楊氏模量E=221 GPa，密度ρ=7 850 kg/m3，泊松比μ=0.3，剛性輪轂半徑R=0.25 m，葉片尺寸180 mm×85 mm×2 mm，輪轂速度Θ=300π rad/s。

表2 不同預置角下旋轉葉片前5階固有頻率Table 2 The first five natural frequencies of rotating blade with respect to the presetting angle

通過表2數(shù)據(jù)的縱向對比可以發(fā)現(xiàn)，隨著預置角的增加，旋轉葉片的各階固有頻率均有增加的趨勢；通過橫向對比可以看出，當截斷級數(shù)取12時，無論在何種預置角下，本文方法計算的數(shù)據(jù)結果與文獻[9]的數(shù)據(jù)能很好的吻合，從而驗證了本方法計算旋轉葉片固有頻率的正確性。

2.3 動力學行為研究

假設葉片由上層的陶瓷氧化鋁與下層的TC4鈦合金組成，二者材料屬性如表3所示。選取尾流激振力中的第k次諧波進行分析，假設n=1，葉片的尺寸為180 mm×85 mm×2 mm，剛性輪轂半徑R=0.25 m。

表3 旋轉功能梯度葉片的材料屬性Table 3 Material properties of functionally graded blade

圖5給出了不同尾流諧波下激振力幅值對旋轉功能梯度葉片均方根速度響應的影響規(guī)律。其中，功能梯度指數(shù)P=0.5，預置角θ=60°。通過圖中3組曲線可以看出，葉片速度響應隨著尾流激振力幅值增大而增大，各階共振峰也逐漸陡峭；同時通過縱向對比發(fā)現(xiàn)尾流激振力諧波的諧次越高，所激起的共振峰越多，且各階共振峰對應的轉速會逐漸向左移動。

圖5 不同幅值激勵力下葉片的速度響應曲線Fig.5 Effect of amplitudes of excitation force on velocity response curves of functionally graded blade

為探究激振力諧次k對旋轉葉片振動特性的影響規(guī)律，作出不同功能梯度指數(shù)P下旋轉葉片的坎貝爾圖，如圖6所示。可以發(fā)現(xiàn)，葉片各階動頻會隨著轉速的增加而增大。轉速越高，葉片受到的離心力越大，使得葉片產生更大的附加拉伸變形，導致葉片抗彎剛度增加。同時可以看出，在該轉速范圍內，第1諧次尾流激振力頻率與葉片動頻沒有交點，因此不會激起葉片共振。從k=2開始，不同諧次的激振力頻率會與葉片動頻相交于不同的點，理論上這些交點所對應的轉速就是葉片的臨界轉速。尾流激振力的諧次越高，越容易與動頻產生交點，同時導致各階臨界轉速逐漸向左移動，圖6所反映的現(xiàn)象與圖5完美契合。同時，在實際工程中，由于旋轉葉片上游的靜葉數(shù)目n較多，導致尾流激振力的基頻較高，因此即便是低諧次的尾流也易與葉片動頻產生交點，導致葉片臨界轉速大多分布在低速區(qū)，因而在起動階段需要保證發(fā)動機很快的通過該危險區(qū)域。

圖6 功能梯度葉片坎貝爾圖Fig.6 Campbell diagrams of functionally graded blade

同時通過圖6的橫向對比可以發(fā)現(xiàn)，當其他參數(shù)相同時，葉片各階動頻與各階臨界轉速均會隨著功能梯度指數(shù)P的增加而降低。為直觀反映變化趨勢，圖7以尾流激振力第6次諧波分量為例，給出不同功能梯度指數(shù)P下旋轉葉片均方根速度隨轉速的變化曲線。當功能梯度指數(shù)為0時，葉片可認為只含有陶瓷氧化鋁，此時葉片的各階固有頻率最大。當功能梯度指數(shù)P增加時，TC4鈦合金在總材料中所占比例增大，葉片材料屬性的改變導致葉片各階固有頻率發(fā)生變化。因此在設計固定尺寸的轉子葉片時，可通過改變葉片材料分布情況來調整葉片固有頻率，使葉片的共振轉速不落在發(fā)動機常用工作轉速范圍內。

圖7 不同功能梯度指數(shù)下葉片的速度響應曲線Fig.7 Effect of different material parameters on velocity response curves of functionally graded blade

圖8給出了不同激振力諧次k下輪轂半徑R對旋轉葉片振動特性的影響規(guī)律。其中功能梯度指數(shù)P=0.5，預置角θ=60°。通過縱向對比可以發(fā)現(xiàn)，k越大，輪轂半徑R對低階共振峰的影響越小。同時可以直觀看出，在特定諧次的尾流激振力下，隨著輪轂半徑R的增大，旋轉葉片各階共振峰會向右移動，且頻率階次越高變化越明顯。這是因為當轉速相同時，旋轉半徑越大，葉片受到的離心力分量也越大，導致葉片的抗彎剛度顯著增加，從而使旋轉葉片各階共振轉速會向著較大的一方移動。同樣，當葉片展弦比(l/b)、厚長比(h/l)等參數(shù)改變時也會出現(xiàn)類似的現(xiàn)象，這里不再給出曲線圖。

圖8 不同輪轂半徑下葉片的速度響應曲線Fig.8 Effect of different radius of rotation on velocity response curves of functionally graded blade

圖9為不同預置角θ下旋轉功能梯度葉片均方根速度隨葉片轉速的變化曲線，其中，功能梯度指數(shù)P=0.5。

圖9 不同預置角下旋轉葉片的速度響應曲線Fig.9 Effect of different presetting angle on velocity response curves of functionally graded blade

可以發(fā)現(xiàn)，k越大，預置角θ對低階共振峰的影響越小。同時在特定諧次的尾流激振力下，隨著預置角θ的增大，各階共振峰向右偏移，且頻率階次越高變化越明顯。這是因為預置角會改變葉片旋轉時的姿態(tài)，葉片姿態(tài)不同，離心力產生的附加拉伸形變也稍有不同。

3 結論

1)本文引入罰函數(shù)法處理葉片的邊界條件，采用改進Fourier級數(shù)對旋轉葉片位移場進行擬合，推導了葉片的能量泛函并利用Ritz方法進行求解，數(shù)值結果算例表明本方法具有良好的計算精度，適用于旋轉葉片振動特性的預測和評估。

2)葉片轉速越高受到的離心力越大，導致葉片抗彎剛度增加，使得葉片各階固有頻率隨轉速的增加而增大。此外，輪轂半徑與葉片預置角也會改變葉片受到的離心力，對葉片的振動特性有不同程度的影響。

3)隨著尾流激振力幅值的增加，葉片速度響應也隨之增大，且尾流激振力的頻率越高，越容易與葉片動頻產生交點，各階臨界轉速也逐漸降低。

4)功能梯度指數(shù)會改變葉片的材料分布，從而導致葉片的固有頻率與臨界轉速發(fā)生變化，在實際工程中可通過改變葉片材料分布情況來調整葉片固有頻率，使葉片的共振轉速不落在發(fā)動機常用工作轉速范圍內。