吳 杰, 楊衛(wèi)東, 虞志浩

(南京航空航天大學(xué) 直升機(jī)旋翼動力學(xué)重點(diǎn)實(shí)驗室，南京 210016)

直升機(jī)旋翼動力學(xué)的挑戰(zhàn)之一是精確地預(yù)測槳葉結(jié)構(gòu)振動載荷。除了結(jié)構(gòu)動力學(xué)建模與氣動力建模精度對載荷有重要影響外，載荷計算方法也同樣是非常重要的因素。Bielawa[1]最早采用力積分法與模態(tài)疊加法兩種載荷計算方法對無鉸式旋翼開始了這方面的研究。文章指出力積分法雖然能得到更準(zhǔn)確的結(jié)果，并且能以更少的模態(tài)收斂，但其實(shí)現(xiàn)過程更為復(fù)雜。Thomas[2]利用CAMRAD對幾種直升機(jī)旋翼槳葉(BO105，CH-34及SA349/2)進(jìn)行了力積分法與曲率法的對比研究。他指出當(dāng)剖面之間的槳葉結(jié)構(gòu)特性相差不大時，曲率法能獲得較好的載荷；而力積分法則依賴積分步長以及槳葉分段數(shù)的選擇，并且彎矩預(yù)測比剪力預(yù)測需要更多的分段數(shù)和更小的積分步長。上世紀(jì)八十年代，反力法在著名的旋翼動力學(xué)綜合分析平臺2GCHAS中被引入。Lim等[3]利用2GCHAS以及CAMRAD/JA研究了UH-60A直升機(jī)結(jié)構(gòu)載荷數(shù)據(jù)。文章認(rèn)為曲率法的缺點(diǎn)在于高階導(dǎo)數(shù)引起的數(shù)值精度損失；同時力積分法很難處理近槳根處的載荷計算，因為槳根處剖面特性變化通常較為劇烈；而反力法采用了有限元方程中組集之前的單元矩陣計算節(jié)點(diǎn)載荷，預(yù)測精度較高。

本文著重比較上述三種載荷計算方法對于槳葉結(jié)構(gòu)載荷的影響。采用剛?cè)狁詈夏Ｐ蚚4]描述旋翼槳葉的動力學(xué)運(yùn)動關(guān)系。該模型將鉸接式旋翼鉸與軸承處的轉(zhuǎn)角抽象成三個獨(dú)立自由度，作為剛體轉(zhuǎn)角與槳葉彈性變形耦合，相較于采用小轉(zhuǎn)角假設(shè)的傳統(tǒng)有限元模型在瞬態(tài)響應(yīng)計算及碰撞研究中具有明顯的優(yōu)勢[5]。氣動力由準(zhǔn)定常氣動模型計算得到，并計入所有剛體轉(zhuǎn)角及彈性變形的影響。由于尾跡模型對于揮舞彎矩的預(yù)測至關(guān)重要[6]，本文采用自由尾跡入流模型計算槳盤入流。剖面氣動力系數(shù)從翼型數(shù)據(jù)表中根據(jù)有效迎角與來流速度插值得到。離散后的氣動力以廣義力的形式表達(dá)，運(yùn)動方程依據(jù)Hamilton原理推導(dǎo)建立。由于傳統(tǒng)積分方法對于瞬態(tài)響應(yīng)積分收斂較慢以及力積分法會引入數(shù)值累積誤差，本文發(fā)展了一種新的基于Richardson外插技術(shù)的自適應(yīng)變步長Legendre-Gauss數(shù)值積分算法，加速收斂速度，降低累積誤差。

1 運(yùn)動方程及求解方法

1.1 坐標(biāo)系

(1)

其中，u,v,w為槳葉在變形前坐標(biāo)系中度量的彈性變形，η,ζ為質(zhì)點(diǎn)位于剖面坐標(biāo)系中的坐標(biāo)。矢量在不同坐標(biāo)系中的坐標(biāo)可經(jīng)由轉(zhuǎn)角余弦矩陣相互轉(zhuǎn)換得到。比如變形前坐標(biāo)系中的矢量坐標(biāo)經(jīng)過轉(zhuǎn)換矩陣TDU轉(zhuǎn)到變形后坐標(biāo)系中，即：rD=TDUrU，其中TDU是關(guān)于Euler角的函數(shù)矩陣。

圖1 鉸與軸承的運(yùn)動

為了充分考慮鉸接式旋翼槳葉在鉸及變距軸承處的慣性運(yùn)動，引入三個剛體自由度，用以模擬鉸與軸承處的槳葉轉(zhuǎn)角，分別是剛體揮舞角(θ1+βp)、剛體擺振角(θ2)以及剛體變距角(θ3+θ0)，如圖1所示。βp指槳葉預(yù)錐角，θ0為操縱變距。θ1、θ2與θ3分別是因慣性力與外力綜合作用引起的揮舞鉸、擺振鉸以及變距軸承的剛體轉(zhuǎn)角。為方便列寫動能變分式，式(1)中的位置矢量通常需要經(jīng)余弦轉(zhuǎn)換矩陣轉(zhuǎn)換到慣性系中。槳葉動能變分項最終由剖面動能變分沿槳葉展向積分得到：

(2)

勢能變分項由中等變形梁理論[7]推導(dǎo)得到。槳葉抽象為細(xì)長柔性Bernoulli梁，變形過程中剖面與槳葉參考軸線始終保持垂直，不計橫向剪切。為了更精確地對彈性扭轉(zhuǎn)建模，本文在勢能變分項推導(dǎo)中未使用階次準(zhǔn)則。文獻(xiàn)[8]指出，保留彈性扭轉(zhuǎn)的高階項(三階或以上)對于保持力積分法與模態(tài)疊加法結(jié)果的一致性非常重要。

與傳統(tǒng)有限元建模方法不同，由于模型中剛?cè)狁詈系拇嬖?，方程的建立與槳葉剖面振動載荷的計算需要作出相應(yīng)的修改。氣動力虛功需要以對應(yīng)于剛體自由度與彈性自由度的廣義力形式給出，因為它對剛體自由度同樣做功。外力變分項可寫成：

(3)

其中，R為質(zhì)點(diǎn)在旋轉(zhuǎn)槳轂坐標(biāo)系中的坐標(biāo)；TRD是矢量從變形后坐標(biāo)系到旋轉(zhuǎn)槳轂坐標(biāo)系中的轉(zhuǎn)換矩陣，因此也是剛體轉(zhuǎn)角以及彈性變形的函數(shù)；φ為彈性扭轉(zhuǎn)角；Fa與Ma分別指變形后坐標(biāo)系中的氣動力與力矩。實(shí)際計算中，它們不僅包含采用準(zhǔn)定常氣動模型得到氣動環(huán)量力與力矩，還考慮了基于Greenberg薄板理論的非環(huán)量力與力矩分量。非環(huán)量部分是迎角變化率及彈性扭轉(zhuǎn)的函數(shù)，提供氣動阻尼，對氣彈響應(yīng)計算中的數(shù)值穩(wěn)定性至關(guān)重要。

1.2 動力學(xué)方程

在得到動能變分、勢能變分以及外力虛功后，槳葉非線性動力學(xué)微分方程組依據(jù)Hamilton原理建立起來。槳葉運(yùn)動由形函數(shù)離散成廣義節(jié)點(diǎn)自由度q的形式(由剛體自由度與Chopra 15自由度[9]組成)。動力學(xué)方程的最終形式可表示為：

(4)

(5)

為避免Jacobi陣的求逆操作，通過求解關(guān)于方程單步迭代差值的線性方程組推進(jìn)牛頓迭代算法[11]，減少了程序計算量。同時，修改的牛頓法只在程序開始時計算Jacobi陣并在隨后的迭代中當(dāng)作常量處理；但在本文中為保證收斂方向，當(dāng)收斂速度變慢時程序會適時更新Jacobi陣。

1.3 載荷計算方法

旋翼動力學(xué)中一般使用三種計算方法預(yù)測槳葉結(jié)構(gòu)載荷。力積分法[1,8,12]將剖面結(jié)構(gòu)載荷沿外段槳葉從計算參數(shù)點(diǎn)(即待求載荷的徑向位置)到槳尖積分得到參考點(diǎn)結(jié)構(gòu)載荷。因而這種方法會累積所有可能的由數(shù)值算法或剖面特性差異引起的誤差。根據(jù)Hamilton原理與牛頓定律的等價性，槳葉慣性載荷可直接由動能變分對應(yīng)的自由度變分分量得到[8]。即慣性載荷可表示為：

(6)

基于這種相對簡單的對應(yīng)關(guān)系，求解動能變分中的非線性項TNL子函數(shù)可在力積分法計算載荷時被再次調(diào)用，極大地降低了載荷程序?qū)崿F(xiàn)的復(fù)雜度。由于質(zhì)量陣也來自于動能變分，一種新的計算慣性載荷的方法可以推導(dǎo)出來：

(7)

其中，Q是由所有槳葉運(yùn)動自由度組成的列向量，在本文中寫成(θ1,θ2,θ3,u,v,w,φ,v′,w′)T的形式，二階導(dǎo)數(shù)為其對應(yīng)加速度項；M為槳葉運(yùn)動自由度對應(yīng)的離散之前的質(zhì)量陣；g為體現(xiàn)式(6)中關(guān)系的抽象函數(shù)。

力積分法計算中，外段剖面對計算參考點(diǎn)x0的彎矩包括兩部分：外段剖面力對計算參考點(diǎn)的力矩以及外段剖面內(nèi)力對于本地軸線形成的彎矩；其中剖面力又由慣性載荷與氣動外力組成。最終彎矩沿槳葉展向積分得到。以變形前坐標(biāo)系中的揮舞彎矩為例：

(8)

其中，Δw=w-w0, Δx=(x+u)-(x0+u0)，下標(biāo)0表示該變量是對應(yīng)于計算參考點(diǎn)處的運(yùn)動變形。因為傳統(tǒng)Legendre-Gauss數(shù)值積分方法在有瞬態(tài)氣動力作用時收斂速度較慢，本文將Richardson外插技術(shù)與之結(jié)合，發(fā)展了一種新的數(shù)值積分算法，以加快積分的收斂速度并盡可能地減少累積誤差。

曲率法[2,4]的思想可以簡潔地解釋為彎矩只與梁的彎曲程度有關(guān)。該方法需要對槳葉彈性變形作二次偏導(dǎo)數(shù)，因而基于傳統(tǒng)模態(tài)法求解響應(yīng)時須保留足夠多的模態(tài)才能達(dá)到合理的精度。本文采用隱式方法在位形空間中求解方程，因此不存在這種因模態(tài)截斷引起的問題。仍以揮舞彎矩為例：

(9)

其中，EIy為揮舞 剛度；EA為拉伸剛度；θt指剖面預(yù)扭角；eζ指剖面重心在剖面坐標(biāo)系中擺振方向的偏移。值得注意的是，與式(8)不同，式(9)中的彎矩是在變形后坐標(biāo)系中度量的；因此在與力積分法比較時，需要將二者的載荷矢量統(tǒng)一到同一坐標(biāo)系中。由于飛行實(shí)測載荷數(shù)據(jù)是在變形后坐標(biāo)系中得到的，實(shí)際計算時式(8)的結(jié)果同樣需要經(jīng)過TDU轉(zhuǎn)換到變形后坐標(biāo)系中，而式(9)則不必。

反力法[12]基于有限元原理，因此結(jié)果取決于動力學(xué)有限元建模的精度。在有限元組集過程中，節(jié)點(diǎn)內(nèi)力在相鄰單元之間是相互抵消的；但在最終方程求解之后，由組集之前的相應(yīng)單元矩陣形成的梁單元方程是非平衡系統(tǒng)，其余量即為相鄰單元作用于節(jié)點(diǎn)的結(jié)構(gòu)載荷。反力法不能有效預(yù)測有限元非節(jié)點(diǎn)位置處的結(jié)構(gòu)載荷。盡管如此，反力法仍是一種計算代價較小的載荷計算方法。節(jié)點(diǎn)反力可寫為：

(10)

其中，下標(biāo)EL指該變量是原始單元矩陣或向量，并非直接取自方程組集之后的總體矩陣或向量。揮舞彎矩是Freaction中對應(yīng)于w′自由度的分量，擺振彎矩即是對應(yīng)于v′自由度的分量。類似的，反力法的載荷結(jié)果是在變形前坐標(biāo)系中度量的，在與別的方法以及實(shí)驗結(jié)果比較時需要進(jìn)行轉(zhuǎn)換。

2 算例與討論

采用BO105模型槳葉[13]作為第一個算例對比在無氣動外力作用下三種載荷計算方法的結(jié)果。BO105直槳葉安裝在無鉸式槳轂上。給定槳尖揮舞方向初速度10 m/s，徑向位置21.6%R處的揮舞彎矩、擺振彎矩以及扭轉(zhuǎn)彎矩分別如圖2/3/4 所示?？梢钥闯觯词乖谶@種劇烈振動狀態(tài)下，三種方法的結(jié)果仍能吻合得很好。這說明了三種方法在對于無鉸式旋翼槳葉純結(jié)構(gòu)振動載荷具有相同的精度。

第二個算例采用SA349/2[15]直升機(jī)飛行實(shí)測載荷數(shù)據(jù)驗證三種載荷計算方法。SA349/2是裝有非線性減擺器的鉸接式旋翼，其槳葉具有先進(jìn)幾何外形(預(yù)扭及非對稱翼型等)及非線性剖面特性(剖面重心偏置等)。選取試驗中的飛行狀態(tài)3 (μ=0.198,CT/σ=0.062) ，比較研究三種載荷計算方法。

圖2 揮舞彎矩，剖面位置21.6%R, BO105

圖5 揮舞彎矩，剖面位置80%R, SA349/2

圖6 擺振彎矩，剖面位置46%R, SA349/2

徑向位置80%R處揮舞彎矩及46%R處擺振彎矩分別如圖5/6所示。三種方法都能預(yù)測到彎矩的大致趨勢。對于揮舞彎矩，力積分法在整個周期上都預(yù)測偏高；并且彎矩隨時間變化的曲線較平滑，相位細(xì)節(jié)有損失。這種預(yù)測偏高的現(xiàn)象同樣在擺振彎矩出現(xiàn)，如圖6所示。對于槳葉中段的擺振彎矩，曲率法與反力法都比力積分法獲得了更為合理的結(jié)果。在槳葉前進(jìn)邊，三種方法都預(yù)測過高。除了線性當(dāng)量阻尼模型對于非線性減擺器建模的不準(zhǔn)確以外，準(zhǔn)定常氣動模型可能不足以預(yù)測這種飛行狀態(tài)下槳葉前行邊氣動載荷的計算。因為中等前飛狀態(tài)下，槳渦干擾的存在降低了準(zhǔn)定常氣動模型的氣動力預(yù)測精度。這種情況需要采用非定常氣動模型[4]或計算流體動力學(xué)模型對氣動力進(jìn)行建模[8,15]。

3 結(jié) 論

本文建立了考慮剛體與彈性運(yùn)動耦合的旋翼剛?cè)狁詈蟿恿W(xué)模型。依據(jù)Hamilton原理與牛頓定律的等價性，發(fā)展了一種易于實(shí)現(xiàn)的力積分方法。為加快積分收斂速度及減少積分誤差，開發(fā)了一種新的基于Richardson外推法的Legendre-Gauss數(shù)值積分算法。著重比較直升機(jī)常用的三種載荷計算方法，以無鉸式BO105模型槳葉以及鉸接式SA349/2直升機(jī)為研究對象，分析對比了三種載荷方法的預(yù)測精度與適應(yīng)范圍?；谏鲜鲅芯浚疚目偨Y(jié)如下：

(1) 在純結(jié)構(gòu)系統(tǒng)給定初始仿真運(yùn)動條件下，由于不考慮氣動力，三種載荷計算方法都能勝任結(jié)構(gòu)振動載荷的預(yù)測，包括槳葉剖面揮舞彎矩、擺振彎矩及扭轉(zhuǎn)力矩。

(2) 曲率法與反力法易于實(shí)現(xiàn)。在有氣動外力作用條件下，二者較力積分法能獲得頻率更為豐富幅值更為合理的載荷結(jié)果。反力法無法預(yù)測有限元非節(jié)點(diǎn)處的結(jié)構(gòu)載荷。依靠相鄰節(jié)點(diǎn)之間的形函數(shù)插值得到的載荷是不準(zhǔn)確的，槳葉建模時需要人為添加有限元節(jié)點(diǎn)。

(3) 氣動模型對于旋翼振動載荷至關(guān)重要，氣動力的預(yù)測能力需要進(jìn)一步提高。從對比結(jié)果看，三種載荷計算方法都未能準(zhǔn)確預(yù)測槳盤前行邊的振動載荷，這與前行邊的槳葉剖面氣動力預(yù)測精確不夠有關(guān)。

參 考 文 獻(xiàn)

[1]Bialawa Richard L. Blade stress calculation-mode deflection vs. force integration[J]. Journal of the American Helicopter Society, 1978, 23(3),10-16.

[2]Maier T H. An examination of helicopter rotor load calculations[C]. The American Helicopter Society National Specialists’ Meeting on Rotorcraft Dynamics, Arlington, Texas, November 1989.

[3]Lim J W, Analytical investigation of UH-60A flight blade airloads and loads data[C]. The 51th Forum of the American Helicopter Society, Fort Worth, Texas, May 1995.

[4]WANG Hao-wen, GAO Zheng. Rotor vibratory load prediction based on generalized forces[J]. Chinese Journal of Aeronautics, 2004,17(1):28-33.

[5]韓東,高正，王浩文,等. 直升機(jī)槳葉剛?cè)狁詈咸匦约坝嬎惴椒ǚ治鯷J]. 航空動力學(xué)報, 2006,21(1):36-40.

HAN Dong, GAO Zheng, WANG Hao-wen. Analysis of the computational method and the helicopter blade characteristics with coupled rigid and elastic Motions [J]. Journal of Aerospace Power, 2006,21(1):36-40.

[6]Bousman W G, Maier T H, An investigation of helicopter rotor blade flap vibratory loads[C]. The 48th Annual Forum of the American Helicopter Society, Washington, D.C., June 1992.

[7]Hodges D H, Dowell E H. Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[R]. NASA TN D-7818, December,1974.

[8]Datta A, Fundamental Understanding, Prediction and validation of rotor vibratory loads in steady level flight[D]. PhD thesis in Aerospace Engineering, 2004.

[9]Chopra I, Sivaneri N T, Aeroelastic stability of rotor blades using finite element analysis[R]. NASA-CR-166389, August 1982.

[10]虞志浩, 楊衛(wèi)東, 張呈林. 基于Broyden法的旋翼多體系統(tǒng)氣動彈性分析[J]. 航空學(xué)報, 2012,33(12):2171-2182.

YU Zhi-hao, YANG Wei-dong, ZHANG Cheng-lin. Aeroelasticity analysis of rotor multibody system based on Broyden method[J]. Chinese Journal of Aeronautics, 2012,33(12):2171-2182.

[11]Ascher U M,Petzold L R. Computer methods for ordinary differential equations and differential-algebraic equations[M]. SIAM press, 1998.

[12]Floros M W, Elastically tailored composite rotor blades for stall alleviation and vibration reduction[D]. PhD thesis in Aerospace Engineering, December 2000.

[13]Bir G, Chopra I. University of maryland advanced rotorcraft code (UMARC) theory manual[M].Vol3.July 20, 1994.

[14]Heffernan R, Gaubert M, Structural and aerodynamic loads and performance measurements of an SA349/2 Helicopter with an Advanced Geometry Rotor[R]. NASA TM-88370, 1986.

[15]Abhishek A, Analysis, Validation, Prediction and fundamental understanding of rotor blade loads in an unsteady maneuver[D]. PhD thesis in Aerospace Engineering, 2010.