帶外掛機翼顫振分析的傳遞函數方法

段靜波， 江 濤

(1. 西北工業大學 航空學院， 西安 710072；2. 軍械工程學院 無人機工程系， 石家莊 050003)

帶外掛機翼顫振分析的傳遞函數方法

段靜波1,2， 江 濤2

(1. 西北工業大學 航空學院， 西安 710072；2. 軍械工程學院 無人機工程系， 石家莊 050003)

將傳遞函數方法應用于帶外掛機翼的顫振分析。利用干凈機翼的彎扭振動微分方程和Therdorson非定常氣動力模型建立了機翼顫振方程。將外掛視為與機翼連接的具有質量及轉動慣量的剛體，且考慮機翼與外掛間的俯仰聯接剛度，通過變形協調和內力平衡等連接條件引入機翼顫振模型。運用傳遞函數方法，將顫振微分方程轉換為狀態空間形式，通過求解復特征值問題獲得了帶外掛機翼的顫振速度和顫振頻率。通過與已有文獻及有限元結果對比，驗證了該方法的正確性和有效性。最后，討論了外掛質量、轉動慣量、外掛位置、掛架俯仰聯接剛度等因素對帶外掛機翼顫振特性的影響。

帶外掛機翼；非定常氣動力；顫振；傳遞函數法

機翼是飛機氣動彈性分析的重要對象。通常機翼下會安裝發動機、導彈等外掛，這將會顯著改變機翼的顫振特性。某些外掛狀態甚至可能導致顫振臨界速度的急劇降低，嚴重影響飛機的安全性。因此，對帶外掛機翼顫振問題一直是飛機設計中受到廣泛重視的課題。國內外學者在此方面進行了諸多研究。20世紀90年代，楊智春等[1]就采用顫振分析重頻理論對帶外掛物二元機翼的顫振頻率及顫振邊界變化的一般規律進行了研究，并進行了實驗驗證。楊智春等[2-3]還建立實驗模型研究俯仰、偏航、側擺等聯接剛度對機翼/外掛系統顫振特性的影響。Yang[4]應用KBM方法的一次、二次近似理論研究了機翼外掛系統以及三角機翼兩種模型的極限環顫振。楊翊仁等[5]按照工程處理的思路，先用一階諧波平衡當量線化方法估計系統極限環顫振頻率，然后引用孤立外掛單自由度系統在簡諧迫振情況下的次諧分叉頻率條件來預估機翼帶外掛系統極限環顫振的次諧響應存在區域。葉煒梁[6]根據顫振運動方程，應用v-g參數法和非定常氣動力的偶極子格網法計算了CKI機翼翼尖帶外掛物時的顫振特性。近年來，Ozcan等[7]采用實驗和數值分析方法對帶外掛機翼的氣動特性進行了研究。Tang等[8]采用實驗的方法對帶外掛三角翼的顫振、極限環振蕩在低速風洞中進行了研究，并與理論結果進行了對比分析。Fazelzadeh等[9]采用擴展伽遼金法研究了飛機滾轉下帶外掛機翼的分叉和顫振問題。Librescu等[10]采用擴展伽遼金法研究了復合材料帶外掛機翼的氣動彈性穩定性及響應問題。Karpel等[11]提出一種新的模態耦合技術，可以高效地進行帶外掛機翼的顫振及氣動伺服彈性分析。邱志平等[12]、Tang等[13]分別研究了帶有間隙型非線性剛度的帶外掛機翼的顫振問題。Chen等[14]提出了一種基于諧波平衡法的增量法來研究帶外掛翼型的極限環顫振問題，并將該問題轉化為極值問題求解。許軍等[15]基于Hamilton原理推導帶外掛機翼的動力學方程，研究了大展弦比帶外掛機翼彎彎扭運動的顫振特性。王鋼林等[16]采用梁單元對雙梁式機翼和發動機吊艙進行了結構建模，研究了機翼在不同攻角時的定常氣動力及發動機推力聯合作用下固有振動和顫振特性的變化情況。

傳遞函數方法是一種基于控制理論的半解析計算方法。該方法求解過程簡潔和統一，邊界條件處理規范和方便，常用于動力學系統的穩定性、動態響應特性分析。Yang[17]最早提出了該方法，并應用于梁、桿等一維結構振動問題。周建平等[18-19]系統地將方法拓展應用于平面應力問題、薄板彎曲問題、厚壁圓柱殼和加筋柱殼等二維、三維的結構振動問題。之后，馮瑩等[20]將傳遞函數方法應用于光波導傳播特性分析。李恩奇[21]基于分布參數傳遞函數方法進行了被動約束層阻尼結構的動力學分析。趙雪川[22]將傳遞函數方法應用于非局部彈性和黏彈性結構的動力學分析。Shen等[23]等基于傳遞函數方法研究了碳納米管的力學問題。由于機翼顫振本質上也是系統動力學穩定性問題，因而，本文在前人研究的基礎上，將傳遞函數方法應用于帶外掛機翼的顫振分析。

1 顫振分析模型

1.1 機翼彎扭振動微分方程

圖1所示為一根長直機翼，其半展長為L，半弦長為b。取固支端與機翼剛軸的交點為原點建立坐標系，y軸沿機翼軸線從翼根指向翼梢，x軸沿機翼弦向由前緣指向后緣，與y軸正交，z軸與x，y軸構成右手系。在此坐標系下，機翼的彎扭振動微分方程可寫為

(1)

式中，h為機翼彎曲振動位移，α為機翼扭轉振動轉角，EI為機翼抗彎剛度，GJ為機翼抗扭剛度，m為機翼單位長度質量，Iα為單位長度機翼繞彈性軸的轉動慣量，xα為機翼彈性軸到機翼截面重心的距離，Lh為機翼單位長度的升力，Tα機翼單位長度的扭矩，y為機翼展向坐標，t為時間。

(a)機翼展向

(b)機翼剖面圖1 長直機翼示意圖Fig.1 A straight aircraft wing

1.2 非定常氣動力模型

在忽略機翼重力影響條件下，機翼顫振時的外力只有氣動力。本文采用片條理論進行非定常氣動力計算。根據Theodorson理論，單位展長的非定常升力與相應的俯仰力矩按下式計算[24]

(2)

由于減縮頻率k是圓頻率ω和空速v的函數，為了后續求解方便 ，將C(k)寫為C(ω,v)，其它符號含義同文獻[24]。將式(2)代入式(1)得到機翼顫振微分方程

(3)

1.3 機翼掛載的處理

在處理機翼掛載時，考慮外掛的質量慣量特性，以及與機翼間的俯仰聯接剛度，而且忽略外掛俯仰運動產生的鉛垂方向(圖1中z軸方向)的位移，忽略外掛氣動力對機翼顫振的影響。

對干凈機翼而言，將機翼可視為一根懸臂梁，其邊界變形協調條件和內力平衡條件為

(4)

當機翼存在掛載時，設沿機翼軸線方向上距翼根y0處存在一個外掛，外掛質量為m0，轉動慣量記I0，外掛距機翼彈性軸的距離為x0，β(t)為外掛相對翼根弦線的俯仰扭轉角，Kβ為掛架俯仰扭轉剛度，如圖2所示。在求解帶外掛機翼顫振問題時，以外掛為界，將機翼分為左右兩部分，中間通過變形協調條件和內力平衡條件聯系兩部分。為進行區分，外掛左側機翼所有物理量增加下標“1”，外掛右側機翼所有物理量增加下標“2”。

圖2 外掛處機翼剖面示意圖Fig.2 Wing section on the position of the external store

基于機翼彎扭振動模型，在外掛處機翼變形應滿足彎曲撓度、彎曲轉角及扭轉角三方面的變形協調條件，具體為

(5)

同時，在掛載處機翼內力應滿足彎矩、剪力及扭矩三方面的內力平衡條件，具體為

(6)

上述邊界條件中引入了一個新自由度β(t)，因而需要補充一個條件。對于外掛，其俯仰振動微分方程為

(7)

2 傳遞函數法求解

對機翼顫振微分方程式(3)進行Fourier變換，并整理可得

(8)

式中，A1(ω,v)、A2(ω,v)、B1(ω,v)、B2(ω,v)的具體表達式如下

(9)

根據傳遞函數方法，定義狀態變量向量如下

(10)

式中，T表示向量轉置。

從而，式(8)可寫成如下狀態空間方程的形式：

(11)

(i=1,2)。

由于式(8)為齊次微分方程組，因而式(11)

中g(ζ,ω)=0。

根據傳遞函數方法，帶外掛機翼邊界與外掛處的變形協調條件、內力平衡條件可以寫為矩陣形式

Mbη(y=0,ω)+Nbη(y=L,ω)+
R0η(y=y0,ω)=γ(ω)

(12)

式中，Mb為機翼根端邊界條件選擇矩陣，Nb為機翼梢端邊界條件選擇矩陣，由于機翼根端固支梢端自由，因而，Mb、Nb的表達式可寫為

R0為機翼外掛處連續條件選擇矩陣，聯立式(5)、式(6)、式(7)可得出其表達式

將Mb、Nb、R0的矩陣式代入式(12)，且根據機翼兩端邊界條件式(4)，可得到

(13)

根據傳遞函數理論，式(11)的解可寫為

(14)

式中，G(y,ζ,ω,v)為狀態空間方程的域內傳遞函數，H(y,ω,v)為狀態空間方程的邊界傳遞函數，其表達式分別為

(15)

式中，變量ζ∈(0,L)，為機翼展向的坐標。

因為g(ζ,ω)=0且γ(ω)=0，由式(14)可進一步得到

(16)

當機翼顫振時，彎曲振動位移h和扭轉振動轉角α的振幅為非零常數，即式(16)的η(y,ω)有非零解，則有：

(17)

令：

(18)

由于A為復矩陣，其行列式值等于零的必要條件為矩陣行列式值的實部與虛部均為零，即

(19)

矩陣A中有空速v和圓頻率ω兩個變量，而式(19)恰好有兩個方程，可以定解。求解上述式(19)，可能會得到滿足方程式(19)的解，即存在多組(v,ω)能滿足方程式(19)。根據機翼顫振時在某一空速時由穩定轉變為不穩定，因而，空速v最小的一組解(v,ω)應為機翼的顫振速度和相應的顫振頻率。

具體求解過程如下：

(j=0,1,2,3,...)；

步驟(3) 取空速v=v0，圓頻率ω依次取，ω0+j△ω(j=0,1,2,3,...)將空速v和圓頻率ω的取值及機翼各物理參數代入式(9)，得到系數A1(ω,v)、A2(ω,v)、B1(ω,v)、B2(ω,v)的值；

步驟(4) 利用系數A1(ω,v)、A2(ω,v)、B1(ω,v)、B2(ω,v)計算式(11)中矩陣F(ω,v)的值；

步驟(5) 利用機翼外掛參數計算矩陣R0；

步驟(6) 將矩陣F(ω,v)、Mb、Nb、R0的值代入式(18)，計算Re[detA]和Im[detA]的值；

步驟(7) 再依次取空速v=v0+j△v(j=1,2,3,...)，重復步驟(3)步驟(6)，計算Re[detA]和Im[detA]的值；

步驟(8) 確定滿足式(19)的空速v和圓頻率ω。

3 算例分析

3.1 正確性驗證

機翼不攜帶外掛實際上是機翼外掛的質量和轉動慣量退化零的情形。為了驗證本文方法和代碼的正確性，首先采用本文方法計算無外掛機翼的顫振特性。

文獻[24]給出一無外掛機翼的顫振特性，機翼物理參數如表1所示，該機翼顫振的計算結果為：機翼顫振速度為35.3 m/s，顫振頻率為24.0 Hz。文獻[24]算例中機翼顫振頻率單位為Hz，本文計算時只需將頻率單位換算一下即可。

表1 機翼的物理參數

圖3 和的等值線圖Fig.3 The contour of Re[det A]and Im[det A]

圖4 和的等值線圖Fig.4 The contour of Re[det A]and Im[det A]

其次，采用本文方法計算機翼攜帶外掛的情形。外掛主要物理參數如表2所示。為了方便方法驗證，設機翼重心到機翼彈性軸的距離xα=0，其它參數取表1中數據。在MSC.Patran中建立帶外掛機翼有限元模型，如圖5所示，采用殼單元，共劃分40個單元，采用集中質量單元模擬機翼外掛的質量及轉動慣量，并利用MSC.Nastran中的氣動彈性分析模塊MSC.FlightLoads進行氣彈分析。

表2 外掛的物理參數

圖5 帶外掛機翼有限元模型Fig.5 The FEM model of the wing with an external store

表3給出了有無外掛情形下機翼固有特性的對比情況。從表3中可以看出，機翼外掛對機翼的固有特性會產生了影響，導致機翼彎曲、扭轉固有頻率均有了降低。特別是，外掛的出現使機翼的一階扭轉頻率從38.03 Hz下降到了29.59 Hz。由于扭轉特性對機翼顫振影響較大，這必將會導致機翼顫振特性的改變。表4給出了本文方法和有限元方法在機翼有無外掛情形下顫振特性的對比情況。從表4中可以看出，無論是機翼有外掛還是無外掛，兩種方法得到的顫振速度與顫振頻率比較吻合。由于MSC.FlightLoads進行氣彈分析時，在結構方面進行了動力學降階，利用了機翼若干低階模態，模態的選用會影響結果的精度。在氣動力模型方面MSC.FlightLoads采用的渦格法，能考慮機翼的三維效應；本文方法在結構方面直接采用微分方程，避免了模態降階，但在氣動力模型方面采用的是片條理論，不能考慮機翼的三維效應。因而，兩種方法計算結果存在一定的差異。

表3 機翼有無外掛固有頻率對比

表4 機翼有無外掛顫振特性對比

3.2 外掛對機翼顫振特性的影響分析

從3.1節可知，帶外掛機翼的顫振特性會隨外掛特性變化而發生顯著變化，因此有必要研究各種外掛參數對機翼顫振品質的影響。本節主要研究外掛質量、轉動慣量、位置以及外掛與機翼間的聯接剛度等因素對機翼顫振特性的影響。分析中所取基本參數見表1和表2。

圖6給出了外掛質量變化對機翼顫振特性的影響。圖中橫坐標為機翼外掛質量m0與機翼單位長度質量m的比值，縱坐標為機翼顫振速度。圖7給出了外掛轉動慣量變化對機翼顫振特性的影響情況。圖中橫坐標為機翼外掛轉動慣量I0與機翼單位長度轉動慣量Iα的比值，縱坐標也為機翼顫振速度。圖6和圖7中均給出了外掛距機翼彈性軸弦向距離x0=-40%b、0、40%b三種下的顫振速度。從兩圖可以看出，采用傳遞函數法得到的外掛質量、轉動慣量對機翼顫振影響是符合常理的，這也表明傳遞函數方法應用于機翼顫振分析是可行有效的。

圖6 機翼外掛質量對顫振速度的影響Fig.6 The influence of the mass of the external store on the velocity of the wing flutter

圖7 機翼外掛轉動慣量對顫振速度的影響Fig.7 The influence of the moment inertia of the external store on the velocity of the wing flutter

由于本文帶外掛機翼顫振模型可以考慮機翼長度特性，因而本文給出了外掛在機翼展向上位置變化時對機翼顫振特性的影響，如圖8所示。圖上橫坐標為機翼外掛展向位置y與機翼半展長L的比值，縱坐標為機翼顫振速度。從圖中可以看出，隨著外掛從翼根移向翼梢，當外掛在機翼彈性軸前時，機翼顫振速度是減小的，而當外掛位于機翼彈性軸上或位于機翼彈性軸后時，機翼顫振速度則是增大的。但是，外掛位于翼根翼中段時，仍然是外掛前于機翼彈性軸布置，有利于機翼顫振穩定性提高。然而，當外掛位于翼梢段時，外掛對機翼顫振速度的影響不再是這樣了。從圖8中虛線圓圈標注處可以看到，在三個弦向位置中，外掛位于x0=0處，機翼的顫振速度最大。這也就是說，外掛位于翼梢時，機翼可能會在弦向某位置出現顫振速度極大值。這與文獻[6]中結論基本一致。因此，對于翼尖外掛，其弦向位置的確定就存在設計空間，需具體問題具體分析計算。

圖8 外掛在機翼展向位置對顫振速度的影響Fig.8 The influence of the position of the external store on the velocity of the wing flutter

外掛掛架剛度也是機翼外掛的重要參數。因此，本文研究了外掛掛架俯仰剛度對機翼顫振的影響，如圖9所示。圖中橫坐標為機翼外掛掛架俯仰剛度Kβ與機翼單位扭轉剛度GJ的比值，縱坐標為機翼顫振速度。當GJ不變，Kβ改變時，從圖中可以看出，當外掛掛架俯仰剛度Kβ較大時，對機翼顫振速度影響不大。但是，當外掛掛架俯仰剛度Kβ較小時，外掛掛架俯仰剛度對機翼顫振影響變得比較明顯，在進行機翼設計時需要加以關注。

圖9 外掛掛架俯仰剛度對顫振速度的影響Fig.9 The influence of pitch stiffness of the external-store pylons on the velocity of the wing flutter

4 結 論

(1)本文將傳遞函數方法應用于帶外掛機翼顫振穩定性問題，通過本文結果與相關文獻算例、MSC. FlightLoads結果對比驗證，表明了本文方法的正確性和有效性；

(2)采用傳遞函數法處理機翼顫振問題時，在模型建立上，可以考慮機翼的長度特性，更近物理實際，有助于更好地揭示帶外掛機翼的顫振機理；在計算過程中，該方法避免了動力學模型的降階工作，而且能得到解析解，有助于提高求解效率；

(3)本文方法可以進一步拓展應用于質量、轉動慣量、弦長等沿展向變化的更為普通的帶外掛機翼顫振穩定性及氣動彈性響應問題。

[1] 楊智春, 趙令誠. 帶外掛二元機翼顫振模態的轉變[J].航空學報，1992，13(9):A552-A554. YANG Zhichun, ZHAO Lingcheng. Transition of flutter mode of two-dimensional wing with external store[J]. Acta Aeronautica et Astronautica Sinica, 1992,13(9):A552-A554.

[2] 楊智春, 趙令誠. 聯接剛度對機翼/外掛系統顫振邊界特性的影響[J].應用力學學報，1993,10(2):1-6. YANG Zhichun, ZHAO Lingcheng. The effective elastic moduli of inhomogeneous media with randomly distributed inclusions[J]. Chinese Journal of Applied Mechanics, 1993,10(2):1-6.

[3] 楊智春,趙令誠,李伶,等. 偏航、側擺聯接剛度對帶外掛三角機翼顫振特性的影響[J].振動工程學報，1992,5(2): 168-172. YANG Zhichun, ZHAO Lingcheng, LI Ling,et al. Effects of pylon yaw and lateral stiffness on the flutter of a delta wing with external store[J]. Journal of Vibration Engineering, 1992,5(2): 168-172.

[4] YANG Y R. KBM method of analyzing limit cycle flutter of a wing with an external store and comparison with a wind-tunnel test[J]. Journal of Sound and Vibration, 1995,187(2): 271-280.

[5] 楊翊仁, 趙令誠. 二元機翼帶外掛系統極限環顫振次諧響應分析[J].航空學報, 1992, 13(7): 410-415. YANG Yiren, ZHAO Lingcheng. Subharmonic response of the limit cycle flutter of wing-store system[J]. Acta Aeronautica et Astronautica Sinica, 1992, 13(7): 410-415.

[6] 葉煒梁. 帶翼梢外掛的CKI機翼顫振分析[J].南京航空航天大學學報,1993,25(4): 554-560. YE Weiliang. Analysis for flutter of CK1 wing with tip external store[J]. Journal of Naijing University of Aeronautics & Astronautics, 1993,25(4)：554-560.

[7] OZCAN O, UNAL M F, ASLAN A R, et al. Aeroelastic chanracteristics of external store configurations at low speeds[J]. Journal of Aircraft, 2011, 32(1):161-171.

[8] TANG D, DOWELL E H. Flutter and limit-cycle-oscillations for a wing-store model with freeplay[J]. Journal of Aircraft, 2006, 43(2):487-503.

[9] FAZELZADEH S A, MARZOCCA P, RASHID E, et al. Effects of rolling maneuver on divergence and flutter of aircraft wing store[J]. Journal of Aircraft, 2010, 47(1):64-71.

[10] LIBRESCU L, SONG O. Dynamics of composite aircraft wings carrying external stores[J]. AIAA Journal, 2008,46(3): 568-578.

[11] KARPEL M, MOULIN B, ANGUITA L, et al. Flutter analysis of aircraft with external stores using modal coupling[J]. Journal of Aircraft, 2004, 47(4):892-901.

[12] 周秋萍,邱志平. 機翼帶外掛系統極限環顫振的區間分析[J]. 航空學報, 2010,31(3): 514-518. ZHOU Qiuping, QIU Zhiping. Interval analysis for limit cycle flutter of a wing with an external store[J]. Acta Aeronautica et Astronautica Sinica, 2010,31(3): 514-518.

[13] TANG D, DOWELL E H. Aeroelastic analysis and experiment for wing/store model with stiction nonlinearity[J]. Journal of Aircraft, 2011, 48(5):1512-1530.

[14] CHEN Y M, LIU J K, MENG G. An incremental method for limit cycle oscillations of an airfoil with an external store[J]. International Journal of Non-linear Mechanics, 2012,47:75-83.

[15] 許軍,馬曉平. 帶外掛機翼的極限環顫振分析[J]. 機械強度，2014,36(6) : 841-845. XU Jun, MA Xiaoping. Limit cycle flutter analysis of a wing with an extrnal store[J]. Journal of Mechanical Strength, 2014,36(6) : 841-845.

[16] 王鋼林，謝長川.考慮翼吊發動機推力的機翼顫振分析[J].航空科學技術,2014,25(6):22-27. WANG Ganglin, XIE Changchuan. Flutter analysis with the thrust effects of engine under wing[J]. Aeronautical Secience & Technology,2014,25(6):22-27.

[17] YANG B. Transfer function of constrained/combined one dimensional contimious dynamic systems[J]. ASME Joumal of Sound and Vibration,1992,156(3):425-443.

[18] 周建平,雷勇軍. 分布參數系統的傳遞函數方法[M]. 北京: 科學出版社, 2006.

[19] 雷勇軍. 結構分析的分布參數傳遞函數方法[D].長沙:國防科技大學,1998.

[20] 馮瑩, 李湘榮, 李海陽, 等. 擴散平面光波導的傳遞函數方法[J].光學學報,1999,19(1):50-56. FENG Ying, LI Xiangrong, LI Haiyang, et al. Tranfer function analysis of diffused planar optical waveguides[J]. Acta Optical Sinica, 1999,19(1):50-56.

[21] 李恩奇. 基于分布參數傳遞函數方法的被動約束層阻尼結構動力學分析[D].長沙:國防科學技術大學,2007

[22] 趙雪川. 非局部彈性和粘彈性結構元件的力學分析[D]. 長沙:國防科技大學,2006.

[23] SHEN Zhibin, LI Xianfang, SHENG Liping,et al. Transverse vibration of nanoyube-based micro-mass sensor via nonlocal Timoshenko beam theory[J]. Computational Material Science, 2012,53:340-346.

[24] 趙永輝. 氣動彈性力學與控制[M]. 北京: 科學出版社, 2006.

Transfer function method for the flutter of aircraft wings with an external store

DUAN Jingbo1,2, JIANG Tao2

(1. College of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China；2. Department of UAVEngineering, Ordnance Engineering College, Hebei 050003, China)

The transfer function method was applied to analyse the flutter of aircraft wings carrying an external store. The flutter differential equation of a clean wing was established by combining the bending-twist vibration equation of the wing and the Therdorson’s unsteady aerodynamics model. The external store hung below the wing by a pitch spring was regarded as a rigid body owning a certain mass and rotary inertia. Then, the influence of the external store on the wing flutter was introduced by considering the conditions of deformation harmony and internal force balance. Further, using the transfer function method, the control equations ware formulated in a state-space form by defining a state vector. Both the flutter velocity and flutter frequency were obtained by solving a complex eigenvalue problem. The results are in good agreement with the literature solutions and the finite element solutions, which indicates that the present method is accurate and efficient. Finally, the effects of the mass, rotary inertia, position and pitch stiffness of the pylons ware investigated.

aircraft wing with an external store; unsteady aerodynamics; flutter; transfer function method

中國博士后科學基金(2014M560803)

2015-11-02 修改稿收到日期： 2016-01-04

段靜波 男，博士后，講師，1982年4月生

V215.3

A

10.13465/j.cnki.jvs.2017.10.018