楊智春，張 惠，谷迎松，宋 淼

(西北工業(yè)大學(xué) 航空學(xué)院結(jié)構(gòu)動(dòng)力學(xué)與控制研究所，西安 710072)

無(wú)人機(jī)機(jī)翼通常采用大展弦比構(gòu)型，具有重量輕、柔性大特點(diǎn)；但其在氣動(dòng)載荷作用下會(huì)產(chǎn)生較大結(jié)構(gòu)變形，致基于結(jié)構(gòu)小變形假設(shè)的常規(guī)線性氣動(dòng)彈性分析方法不再適用。大展弦比機(jī)翼氣動(dòng)彈性分析須考慮結(jié)構(gòu)大變形導(dǎo)致的幾何非線性效應(yīng)。Kim等[1]通過(guò)采用跨音速小擾動(dòng)理論與大變形梁理論耦合分析大展弦比機(jī)翼幾何非線性顫振特性，較好預(yù)測(cè)顫振邊界。Dowell[2]對(duì)陣風(fēng)響應(yīng)下大展弦比機(jī)翼氣動(dòng)彈性特性進(jìn)行理論分析與風(fēng)洞實(shí)驗(yàn)研究，認(rèn)為陣風(fēng)載荷較小時(shí)幾何非線性對(duì)垂向、扭轉(zhuǎn)響應(yīng)影響較小。Patil等[3]通過(guò)研究大展弦比機(jī)翼靜/動(dòng)氣動(dòng)彈性特性，認(rèn)為大變形幾何非線性對(duì)大展弦比機(jī)翼靜/動(dòng)氣動(dòng)彈性特性有重要影響。Choi等[4]采用MSC.Nastran與CFD軟件耦合分析翼型(不同根梢比及后掠角)對(duì)機(jī)翼氣動(dòng)彈性特性影響，認(rèn)為弦長(zhǎng)不變具有錐形外截面翼型的升力系數(shù)更高、氣動(dòng)彈性特性更優(yōu)。安效民等[5]用近似能量守恒幾何大變形下殼體結(jié)構(gòu)非線性動(dòng)態(tài)響應(yīng)算法結(jié)合求解雷諾平均N-S方程的CFD求解器發(fā)展的高效耦合格式，模擬流場(chǎng)非線性及結(jié)構(gòu)非線性氣動(dòng)彈性響應(yīng)。謝長(zhǎng)川等[6]用線性化方法計(jì)算大展弦比機(jī)翼飛機(jī)在平飛設(shè)計(jì)載荷及陣風(fēng)載荷作用下非線性靜變形，分析結(jié)構(gòu)幾何非線性對(duì)大展弦比機(jī)翼固有振動(dòng)特性及顫振特性影響。劉湘寧等[7]研究結(jié)構(gòu)非線性對(duì)復(fù)合材料大展弦比機(jī)翼顫振特性影響結(jié)果表明，用非線性模型計(jì)算所得顫振速度小于線性模型所得顫振速度。

由以上研究知，大展弦比機(jī)翼幾何非線性顫振由機(jī)翼在定常氣動(dòng)載荷作用下大變形所致。求解機(jī)翼非線性顫振特性前須先求解機(jī)翼在定常氣動(dòng)載荷作用下的靜變形。為此，本文提出考慮靜氣動(dòng)彈性效應(yīng)及大變形幾何非線性的大展弦比機(jī)翼氣動(dòng)彈性分析方法。即在給定飛行狀態(tài)及機(jī)翼初始構(gòu)型下求解機(jī)翼結(jié)構(gòu)的定常氣動(dòng)載荷；采用考慮靜氣動(dòng)彈性效應(yīng)迭代方法求解大展弦比機(jī)翼結(jié)構(gòu)非線性靜變形，在收斂的大變形構(gòu)型下獲得考慮機(jī)翼結(jié)構(gòu)幾何非線性效應(yīng)的剛度矩陣；再進(jìn)行機(jī)翼的準(zhǔn)線性模態(tài)及顫振特性分析。本文用MATALAB軟件平臺(tái)對(duì)MSC.Nastran軟件二次開發(fā)，實(shí)現(xiàn)計(jì)算流程，并對(duì)大展弦比機(jī)翼模型[3]進(jìn)行分析，計(jì)算該機(jī)翼非線性氣動(dòng)彈性特性，所得結(jié)果與文獻(xiàn)[3]一致，驗(yàn)證本文所建分析方法及平臺(tái)的正確性。

1 大展弦比機(jī)翼非線性氣動(dòng)彈性特性分析

1.1 考慮大變形幾何非線性結(jié)構(gòu)剛度矩陣

結(jié)構(gòu)幾何非線性對(duì)大展弦比機(jī)翼氣動(dòng)彈性特性影響主要體現(xiàn)于結(jié)構(gòu)剛度特性隨載荷工況改變及機(jī)翼非線性變形會(huì)影響機(jī)翼的氣動(dòng)力分布。機(jī)翼幾何非線性大變形中結(jié)構(gòu)應(yīng)變?nèi)詾樾×?因此用有限元法對(duì)大展弦比機(jī)翼結(jié)構(gòu)進(jìn)行幾何非線性大變形分析時(shí)材料應(yīng)力應(yīng)變關(guān)系仍為線性，其平衡方程由變形后幾何關(guān)系描述，但結(jié)構(gòu)剛度特性需用切線剛度矩陣描述[8],即

K=K0+Kσ+KL

(1)

式中：K0為小位移(彈性)剛度矩陣，與單元節(jié)點(diǎn)位移無(wú)關(guān)；KL為大位移剛度矩陣，與單元節(jié)點(diǎn)位移有關(guān)；Kσ為初應(yīng)力剛度矩陣。

1.2 機(jī)翼結(jié)構(gòu)準(zhǔn)線性模態(tài)分析

機(jī)翼在大變形平衡位置附近作微幅振動(dòng)時(shí)可認(rèn)為振動(dòng)仍為線性，其振動(dòng)方程為

(2)

式中：x為機(jī)翼結(jié)構(gòu)偏離大變形平衡位置位移列陣；M為機(jī)翼結(jié)構(gòu)質(zhì)量矩陣；K為機(jī)翼結(jié)構(gòu)在大變形平衡位置處切線剛度矩陣。

設(shè)機(jī)翼微振動(dòng)為簡(jiǎn)諧振動(dòng)：

x=Xeiωt

(3)

代入式(1)得：

KX=ω2MX

(4)

式中：ω2為特征值；X為特征向量。分別對(duì)應(yīng)機(jī)翼在大變形平衡位置處的準(zhǔn)線性模態(tài)頻率及模態(tài)振型。

1.3 機(jī)翼氣動(dòng)載荷

對(duì)大展弦比機(jī)翼分析時(shí)可采用亞音速偶極子格網(wǎng)法計(jì)算其氣動(dòng)載荷，在各升力面網(wǎng)格1/4弦線處布置一馬蹄渦，用加速度勢(shì)偶極子模擬機(jī)翼氣動(dòng)力非定常項(xiàng)。每一氣動(dòng)網(wǎng)格單元的非定常壓差為

(5)

1.4 非線性靜氣動(dòng)彈性變形分析

靜氣動(dòng)彈性變形分析基本方程[9]為

(6)

分析單獨(dú)機(jī)翼靜氣動(dòng)彈性變形時(shí)用根部固支邊界條件約束機(jī)翼結(jié)構(gòu)剛體自由度。氣動(dòng)面定常氣動(dòng)力為

(7)

本文通過(guò)Matlab編程由MSC.Nastran軟件的靜氣動(dòng)彈性計(jì)算結(jié)果中提取。

進(jìn)行機(jī)翼非線性靜變形分析時(shí)需將空氣動(dòng)力通過(guò)樣條插值矩陣變換為作用在結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)等效力完成。樣條矩陣Gkg為由結(jié)構(gòu)網(wǎng)格點(diǎn)位移ug到氣動(dòng)力網(wǎng)格點(diǎn)位移uk的變換，可表示為

uk=Gkgug

(8)

據(jù)虛功原理，空氣動(dòng)力Fk及作用在結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)的等效力Fg在相應(yīng)虛位移上做功相等，得：

(9)

式中：δuk，δug為虛位移。

將式(8)代入式(9)得：

(10)

據(jù)虛位移任意性，得：

(11)

1.5 顫振特性分析

在機(jī)翼結(jié)構(gòu)大變形平衡位置處準(zhǔn)模態(tài)分析基礎(chǔ)上可得機(jī)翼在模態(tài)坐標(biāo)系下的顫振運(yùn)動(dòng)方程為

(12)

方程(12)的解可表示為

(13)

(14)

在顫振臨界點(diǎn)，機(jī)翼小擾動(dòng)振動(dòng)為簡(jiǎn)諧運(yùn)動(dòng),即

(15)

由式(14)可得顫振行列式為

(16)

采用p-k法求解廣義特征值可得考慮機(jī)翼大變形幾何非線性的顫振臨界速度及顫振頻率。

2 考慮結(jié)構(gòu)幾何非線性氣動(dòng)彈性特性分析

為使分析過(guò)程自動(dòng)化便于工程應(yīng)用，本文在Matlab軟件平臺(tái)上對(duì)MSC.Nastran軟件進(jìn)行二次開發(fā)，實(shí)現(xiàn)大展弦比機(jī)翼大變形幾何非線性氣動(dòng)彈性分析方法計(jì)算流程。分析過(guò)程有三步，即考慮大變形幾何非線性氣動(dòng)彈性靜變形分析、大變形靜平衡位置處準(zhǔn)模態(tài)分析及大變形靜平衡位置處顫振分析。分析流程見圖1。計(jì)算過(guò)程為：① 在MSC.Patran軟件中建立大展弦比機(jī)翼結(jié)構(gòu)模型及氣動(dòng)力模型；設(shè)定機(jī)翼初始飛行狀態(tài)(來(lái)流速度、飛行高度、馬赫數(shù)等)，調(diào)用MSC.Nastran軟件進(jìn)行靜氣動(dòng)彈性變形分析與氣動(dòng)載荷分析。② 用DAMP語(yǔ)言編程，將氣動(dòng)面的定常氣動(dòng)力轉(zhuǎn)化為每個(gè)結(jié)構(gòu)節(jié)點(diǎn)的氣動(dòng)載荷，將氣動(dòng)載荷與慣性載荷輸出為各節(jié)點(diǎn)的外力載荷。③ 據(jù)所得節(jié)點(diǎn)外力載荷，用MSC.Nastran軟件對(duì)結(jié)構(gòu)進(jìn)行大變形幾何非線性分析，獲得節(jié)點(diǎn)靜變形。④ 判斷變形是否收斂，若此次非線性計(jì)算所得節(jié)點(diǎn)靜變形與前次計(jì)算結(jié)果相差達(dá)到預(yù)設(shè)精度，迭代收斂；否則，對(duì)變形機(jī)翼更新氣動(dòng)面計(jì)算網(wǎng)格，轉(zhuǎn)到步驟①繼續(xù)迭代。⑤ 變形收斂后提取非線性靜平衡位置處切線剛度矩陣、質(zhì)量矩陣，對(duì)機(jī)翼進(jìn)行靜平衡位置處準(zhǔn)模態(tài)分析。⑥據(jù)圖1顫振分析模塊計(jì)算獲得考慮大變形幾何非線性的機(jī)翼顫振速度。若顫振速度與靜氣動(dòng)彈性計(jì)算時(shí)來(lái)流速度不匹配，則更新靜氣動(dòng)彈性分析中速度進(jìn)行迭代計(jì)算，顫振速度與靜氣動(dòng)彈性變形計(jì)算時(shí)速度一致時(shí)，該速度即為給定初速度攻角的匹配顫振速度。

圖1 考慮幾何非線性效應(yīng)的機(jī)翼氣動(dòng)彈性特性分析流程

3 考慮結(jié)構(gòu)幾何非線性的氣動(dòng)彈性特性分析算例

3.1 機(jī)翼的線性顫振分析

為對(duì)比說(shuō)明大變形幾何非線性對(duì)大展弦比機(jī)翼臨界顫振特性影響、驗(yàn)證氣動(dòng)網(wǎng)格的準(zhǔn)確性，對(duì)大展弦比機(jī)翼模型[3]的線性模態(tài)特性及顫振特性進(jìn)行研究。模型數(shù)據(jù)見表1。機(jī)翼有限元模型用20個(gè)梁?jiǎn)卧?，機(jī)翼質(zhì)量用分布質(zhì)量。不考慮機(jī)翼大變形幾何非線性因素，直接用線性振動(dòng)分析所得機(jī)翼振動(dòng)頻率與文獻(xiàn)[3]結(jié)果比較見表2，兩者吻合較好。

表1 機(jī)翼模型數(shù)據(jù)[3]

表2 線性振動(dòng)分析模態(tài)頻率(rad/s)

機(jī)翼線性顫振特性見表3，與文獻(xiàn)[3]結(jié)果比較看出，顫振頻率基本相同，顫振速度相對(duì)誤差5.7%。分析認(rèn)為該誤差由本文計(jì)算非定常氣動(dòng)力時(shí)用亞音速偶極子網(wǎng)格法而文獻(xiàn)[3]用簡(jiǎn)化片條氣動(dòng)力理論所致。

表3 線性顫振特性比較

3.2 機(jī)翼氣動(dòng)彈性特性分析

利用文獻(xiàn)[3]算例模型,考慮大變形幾何非線性因素,采用本文方法進(jìn)行分析。機(jī)翼顫振速度、顫振頻率隨攻角變化曲線見圖2(a)、(b)。文獻(xiàn)[3]認(rèn)為根部攻角在0.61°附近時(shí)顫振速度有一個(gè)跳躍，跳躍點(diǎn)后機(jī)翼顫振速度隨攻角變大而平穩(wěn)降低。本文所得翼梢變形隨攻角變化規(guī)律見圖2(c)。分析結(jié)果與文獻(xiàn)[3]變化趨勢(shì)一致。數(shù)值上所存差異原因同前分析。

圖2 機(jī)翼根部攻角對(duì)非線性氣動(dòng)彈性特性影響

由算例結(jié)果表明，采用本文的基于MSC.Nastran軟件平臺(tái)、考慮結(jié)構(gòu)大變形幾何非線性效應(yīng)大展弦比機(jī)翼氣動(dòng)彈性特性分析方法求解大展弦比柔性機(jī)翼的幾何非線性顫振可行，亦為工程提供簡(jiǎn)便有效手段。

4 結(jié) 論

本文通過(guò)建立考慮大變形幾何非線性效應(yīng)的大展弦比機(jī)翼氣動(dòng)彈性特性分析方法，研究大變形幾何非線性效應(yīng)對(duì)機(jī)翼顫振特性影響，結(jié)論如下：

(1) 顫振特性隨根部攻角穩(wěn)定變化后，大展弦比柔性機(jī)翼大變形所致幾何非線性效應(yīng)會(huì)使機(jī)翼顫振臨界速度顯著降低；

(2) 翼梢變形隨攻角增大而增大，機(jī)翼扭轉(zhuǎn)頻率隨翼梢變形增大而降低，機(jī)翼幾何非線性顫振臨界速度隨攻角增大而降低。

[1] Kim K S, Lee I, Lee H K. Efficient numerical aeroelastic analysis of a high-aspect-ratio wing considering geometric nonlinearity[J]. Journal of Aircraft, 2010, 47(1):338-342.

[2] Dowell E H. Gust response for flexibly suspended high-aspect ratio wings [J]. AIAA, 2010, 48(10):2430-2444.

[3] Patil M J, Hodges D H. Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft [J]. Journal of Aircraft, 2001,38(1):88-94.

[4] Choi S, Park K,Kim J,et al. Aeroelastic characteristics of high-aspect-ratio wing according to wing shape[J]. AIAA, 2013,143:1-8.

[5] 安效民,徐敏. 一種幾何大變形下的非線性氣動(dòng)彈性求解方法[J].力學(xué)學(xué)報(bào),2011,43(1):98-104.

AN Xiao-min, XU Min. An improved geometrical nonlinear algorithm and its application for nonlinear aeroelasticity [J].Chinese Journal of Theoretical and Applied Mechanics, 2011,43(1):98-104.

[6] 謝長(zhǎng)川,楊超. 大展弦比飛機(jī)幾何非線性氣動(dòng)彈性穩(wěn)定性的線性化方法[J]. 中國(guó)科學(xué):技術(shù)科學(xué),2011,41(3):385-393.

XIE Chang-chuan, YANG Chao. A linearization method to analyze geometrical nonlinear aeroelastic stability of the high-aspect-ratio aircraft[J]. Science China,2011,41(3):385-393.

[7] 劉湘寧, 向錦武.大展弦比復(fù)合材料機(jī)翼的非線性顫振分析[J]. 航空學(xué)報(bào), 2006, 27(2):214-218.

LIU Xiang-ning, XIANG Jin-wu. Study of nonlinear flutter of the high-aspect-ratio composite wing[J].Chinese Journal of Aeronautics, 2006, 27(2):214-218.

[8] 呂和祥,蔣和洋. 非線性有限元[M]. 大連：化學(xué)工業(yè)出版社, 1991.

[9] Rodden W P, Johnson E H. MSC/NASTRAN version 68 aeroelastic analysis users guide[M]. New York: The MacNeal-Schwendler Corporation, 1994.

[10] Dowell E H,Traybar J,Hodges D H. An experimental-theoretical correlation study of non-linear bending and torsion deformations of a cantilever beam[J].Journal of Sound and Vibration, 1977, 50(4):533-544.