許 成,員?，|

(南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016)

高空長(zhǎng)航時(shí)無(wú)人機(jī)可晝夜長(zhǎng)時(shí)間巡航，執(zhí)行空中偵察、監(jiān)控、科研探測(cè)等任務(wù)，具有廣闊的應(yīng)用前景。該類(lèi)飛機(jī)一般采用大展弦比機(jī)翼結(jié)構(gòu)，在氣動(dòng)力作用下機(jī)翼變形明顯，幾何非線性特征顯著，非線性氣動(dòng)彈性分析技術(shù)成為研究該類(lèi)結(jié)構(gòu)氣動(dòng)彈性問(wèn)題的有效工具。

在幾何非線性氣動(dòng)彈性問(wèn)題的研究方面，Tang和Dowell[1]以線性化的方法對(duì)簡(jiǎn)化的梁模型進(jìn)行非線性顫振分析，結(jié)果與實(shí)驗(yàn)值吻合良好，表明該方法可以較好的預(yù)估顫振發(fā)生時(shí)的速度。謝長(zhǎng)川、楊超[2，3]等人借助動(dòng)力學(xué)線化的方式，對(duì)大展弦比機(jī)翼展開(kāi)氣動(dòng)彈性穩(wěn)定性分析，與線性分析結(jié)果對(duì)比表明，由于幾何非線性的影響，線性顫振計(jì)算方法不能準(zhǔn)確地預(yù)測(cè)臨界顫振速度。

除采取分步線化方法實(shí)現(xiàn)非線性氣動(dòng)彈性分析之外，也可采用非線性氣動(dòng)響應(yīng)方法進(jìn)行氣動(dòng)彈性響應(yīng)分析。冉玉國(guó)等[4]通過(guò)修改NASTRAN軟件求解序列，在非線性瞬態(tài)響應(yīng)求解模塊中引入氣動(dòng)力計(jì)算模塊，成功實(shí)現(xiàn)了非線性氣動(dòng)彈性響應(yīng)分析功能。崔鵬、韓景龍等[5]利用CFD/CSD方法，對(duì)切尖三角翼模型進(jìn)行非線性顫振研究，相對(duì)于實(shí)驗(yàn)結(jié)果，取得了良好的計(jì)算精度。

雖然CFD/CSD的方法是一種精確的方法，然而CFD方法計(jì)算量巨大，效率較低，本文嘗試采用擬合氣動(dòng)力與高精度動(dòng)響應(yīng)求解器結(jié)合的方法，建立起一套求解方法，以實(shí)現(xiàn)對(duì)幾何非線性氣動(dòng)彈性響應(yīng)問(wèn)題的高效求解。

1 氣動(dòng)彈性響應(yīng)分析

1.1 非線性有限元列式及動(dòng)響應(yīng)求解

本文構(gòu)造一種厚薄通用的四邊形板殼單元作為有限元建模的基本單元。該單元由板元、膜元兩部分組合而成。其中，板元部分通過(guò)合理假設(shè)剪切應(yīng)變場(chǎng)的方式避免出現(xiàn)剪切閉鎖現(xiàn)象[6]；膜元部分以Allman型二次膜位移插值模式為參考[7]，將面內(nèi)位移場(chǎng)寫(xiě)成雙線性協(xié)調(diào)函數(shù)與節(jié)點(diǎn)旋轉(zhuǎn)自由度確定的高次位移場(chǎng)之和的形式[8]，不僅提高了單元性能，而且增加了一法向轉(zhuǎn)動(dòng)自由度。

圖1 節(jié)點(diǎn)自由度的定義Fig.1 Definition of node DOF

根據(jù)Mindlin板理論及von Karman大撓度理論，將殼元對(duì)應(yīng)的Green應(yīng)變寫(xiě)成線性與非線性之和形式：

(1)

將上式拆散記為：

(2)

其中εm、εb和εs分別表示單元的膜應(yīng)變、彎扭應(yīng)變和剪切應(yīng)變。

根據(jù)式(1)，按照標(biāo)準(zhǔn)有限元推導(dǎo)過(guò)程，可得單元的幾何矩陣為：

B=B0+Bl

(3)

其中，線性部分B0=[BmBbBs]T，非線性部分Bl=AG=A[G1G2G3G4]。

在以上兩式中：

Bm=[Bm1Bm2Bm3Bm4]

(4)

Bb=[Bb1Bb2Bb3Bb4]

(5)

(6)

為了避免板殼單元在用于薄板分析時(shí)出現(xiàn)的剪切閉鎖現(xiàn)象[8]，此處重新假設(shè)剪切應(yīng)變場(chǎng)，經(jīng)過(guò)修正，板殼單元的剪應(yīng)變幾何矩陣表示為：

(7)

應(yīng)用時(shí)只需將B0中的Bs替換掉即可，其中Ns、Xs、Ys和Γ*在文獻(xiàn)[6,7]中已詳細(xì)給出。

動(dòng)響應(yīng)計(jì)算過(guò)程中，采用修正的UL方法進(jìn)行時(shí)間步推進(jìn)，即在每一時(shí)間步開(kāi)始時(shí)刻進(jìn)行參考構(gòu)型更新，由于時(shí)間步內(nèi)部迭代參考構(gòu)型不變，所以時(shí)間步內(nèi)部應(yīng)力可直接疊加，記S為Kirchhoff應(yīng)力：

S(k+1)=S(k)+△S

(8)

其中，D為組集后的彈性矩陣，△q為結(jié)構(gòu)節(jié)點(diǎn)位移增量。

每次時(shí)間推進(jìn)時(shí)，首先要將Kirchhoff應(yīng)力轉(zhuǎn)換為Cauchy應(yīng)力。

在Mindlin板假設(shè)與大變形情形下，板單元的切線剛陣由單元線性剛陣、初位移剛陣和初應(yīng)力剛陣組成，形式如下[8]：

(9)

將單元的切線剛度矩陣組集為整體切線剛度矩陣，代入系統(tǒng)的動(dòng)力學(xué)方程：

(10)

采用Newmark方法內(nèi)部嵌套Newton-Raphon迭代的方式，求解上述非線性動(dòng)力方程組[9]，將上式改寫(xiě)為等效靜力形式：

Keff△u=Peff

(11)

其中，Keff、Peff分別為有效切線剛陣、有效不平衡力，R為外載荷，F(xiàn)為節(jié)點(diǎn)等效內(nèi)力。

1.2 非定常氣動(dòng)力的時(shí)域耦合

對(duì)于高空長(zhǎng)航時(shí)無(wú)人飛行器來(lái)說(shuō)，由于其巡航速度一般處于亞音速范圍，機(jī)翼所受的氣動(dòng)力可視為線性氣動(dòng)力，并且可忽略非定常氣動(dòng)力的非平面效應(yīng)，從而可以選用偶極子格網(wǎng)法計(jì)算其頻域下的非定常氣動(dòng)力，再利用有理函數(shù)擬合轉(zhuǎn)換到時(shí)域表達(dá)。 作用在結(jié)構(gòu)有限元模型上的氣動(dòng)力表達(dá)式可記為：

(12)

(13)

(14)

其中：

yi(t+△t)=

i=1,2,3,4

(15)

至此，氣動(dòng)力表達(dá)式完全轉(zhuǎn)化到時(shí)域。

1.3 流固緊耦合實(shí)現(xiàn)

采用緊耦合方式在時(shí)域下對(duì)氣動(dòng)力、結(jié)構(gòu)動(dòng)力學(xué)方程分別求解，在每一時(shí)間步內(nèi)固結(jié)迭代，直到結(jié)構(gòu)變形與氣動(dòng)力變化到達(dá)平衡狀態(tài)，再解除時(shí)間步固結(jié)，向下一時(shí)間步推進(jìn)。時(shí)間步的推進(jìn)分為真實(shí)推進(jìn)(解除時(shí)間步固結(jié)后的推進(jìn))和假想推進(jìn)(固結(jié)的時(shí)間步內(nèi)部假想的推進(jìn))，具體的計(jì)算過(guò)程如下：

(1)初始化氣動(dòng)力求解系統(tǒng)，給定滯后項(xiàng)和實(shí)系數(shù)矩陣，給出初始y值；

(3)求解tn(m)(tn時(shí)刻的第m次配平迭代)時(shí)刻氣動(dòng)力向量Qn(m)(m=0,1,2,…)；

(4)將求得的氣動(dòng)力向量Qn(m)傳遞到結(jié)構(gòu)響應(yīng)系統(tǒng)，將結(jié)構(gòu)響應(yīng)求解推進(jìn)到tn(m+1)時(shí)刻；

簡(jiǎn)要計(jì)算流程：

圖2 程序流程簡(jiǎn)圖Fig.2 Flowchart of the program

2 算例

算例模型為一大展弦比機(jī)翼，選取NACA0012翼型，機(jī)翼的結(jié)構(gòu)有限元模型如圖3所示，機(jī)翼根部固支，半展長(zhǎng)為2.97 m，平均弦長(zhǎng)0.594 m，展弦比為10。在翼稍后緣結(jié)點(diǎn)作用一瞬態(tài)激勵(lì)，時(shí)間為0.01 s，大小為100 N，響應(yīng)分析時(shí)間步長(zhǎng)選取0.001 s。

圖3 機(jī)翼結(jié)構(gòu)的有限元模型Fig.3 The FEM model of wing

模態(tài)頻率/Hz一階7.16二階26.31三階39.761

分別以本文的氣彈響應(yīng)計(jì)算方法與傳統(tǒng)的Nastran顫振分析方法對(duì)上例進(jìn)行線性顫振分析。使用本文方法時(shí)，主體流程不變，在生成單元?jiǎng)偠葧r(shí)不考慮非線性項(xiàng)，即可對(duì)機(jī)翼模型做線性氣動(dòng)彈性響應(yīng)計(jì)算。

圖4為翼端節(jié)點(diǎn)在不同飛行速度時(shí)的線性響應(yīng)。對(duì)比線性氣彈響應(yīng)的結(jié)果與Nastran顫振分析的計(jì)算結(jié)果(圖5)可見(jiàn)，所得顫振速度幾乎相同，可見(jiàn)本文方法具有良好的工程適用性。

圖4 翼端節(jié)點(diǎn)在不同飛行速度時(shí)的線性響應(yīng)Fig.4 The liner responses of tip trailing point at certain velocities

圖5 Nastran顫振分析結(jié)果Fig.5 The result of Nastran flutter anlysis

仍然采用本文方法，在生成單元?jiǎng)傟嚂r(shí)考慮非線性項(xiàng)，算例不變，在翼端后緣結(jié)點(diǎn)作用瞬態(tài)激勵(lì)，時(shí)間為0.01 s，大小為100 N，翼端后緣節(jié)點(diǎn)位移響應(yīng)如圖6。

在來(lái)流速度不變(V=120 m/s)的情況下，給出不同大小的初始激振力f，翼端后緣點(diǎn)位移響應(yīng)結(jié)果如圖7。

由圖6可見(jiàn)，在給定來(lái)流速度的條件下，非線性氣動(dòng)彈性響應(yīng)幅值逐漸增大，在幾何非線性的作用下，響應(yīng)幅值逐漸穩(wěn)定，最終做等幅振蕩，系統(tǒng)出現(xiàn)極限環(huán)振蕩現(xiàn)象。來(lái)流速度越大，其穩(wěn)態(tài)振蕩幅值越大，且隨著來(lái)流速度的增加，系統(tǒng)進(jìn)入極限環(huán)響應(yīng)的時(shí)間越來(lái)越短。另外，在來(lái)流速度不變，初始擾動(dòng)不同的情況下，非線性氣動(dòng)彈性響應(yīng)幅值從不同的初始擾動(dòng)出發(fā)，最終會(huì)回到相同的等幅振蕩狀態(tài)，如圖7所示。可見(jiàn)，極限環(huán)響應(yīng)的穩(wěn)態(tài)幅值與系統(tǒng)的初始擾動(dòng)無(wú)關(guān)。

圖6 翼端節(jié)點(diǎn)在不同飛行速度時(shí)的非線性響應(yīng)Fig.6 The noliner responses of tip trailing point at certain velocities

3 結(jié)論

本文借助Nastran軟件求解并導(dǎo)出非定常氣動(dòng)力影響系數(shù)，通過(guò)Roger擬合，得到時(shí)域下的非定常氣動(dòng)力，并利用非線性有限元編程，采用緊耦合的方式將非定常氣動(dòng)力加載到結(jié)構(gòu)上，實(shí)現(xiàn)了大展弦比機(jī)翼的非線性氣動(dòng)彈性響應(yīng)計(jì)算。結(jié)果表明，本文提出的方法可有效預(yù)測(cè)機(jī)翼顫振的臨界速度?？紤]幾何非線性影響時(shí)，機(jī)翼的非線性氣動(dòng)彈性響應(yīng)系統(tǒng)會(huì)發(fā)生極限環(huán)振蕩現(xiàn)象，其振蕩幅值與來(lái)流速度有關(guān)，與初始狀態(tài)無(wú)關(guān)。

圖7 翼端節(jié)點(diǎn)在不同初始狀態(tài)下的非線性響應(yīng)Fig.7 The noliner responses of tip trailing point at certain initial state

參考文獻(xiàn)：

[1] Tang D, Dowell E H. Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings[J].AIAA journal, 2001,39(8):1 430-1 441.

[2] Xie C C, Leng J Z, Yang C. Geometrical nonlinear aeroelastic stability analysis of a composite high-aspect-ratio wing[J].Shock and Vibration, 2008,15(3):325-333.

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

[4] 冉玉國(guó),劉會(huì),韓景龍.大展弦比機(jī)翼的非線性氣動(dòng)彈性響應(yīng)分析[J].空氣動(dòng)力學(xué)學(xué)報(bào),2009,27(4):394-399.

[5] 崔鵬,韓景龍.基于 CFD/CSD 的非線性氣動(dòng)彈性分析方法[J].航空學(xué)報(bào),2010,31(3):480-486.

[6] 岑松,龍志飛.對(duì)轉(zhuǎn)角場(chǎng)和剪應(yīng)變場(chǎng)進(jìn)行合理插值的厚薄板通用四邊形單元[J].工程力學(xué),1999,16(4):1-15.

[7] Ibrahimbegovic A. A novel membrane finite element with an enhanced displacement interpolation[J].Finite elements in analysis and design, 1990,7(2):167-179.

[8] 朱菊芬,鄭罡.帶旋轉(zhuǎn)自由度C°類(lèi)任意四邊形板 (殼) 單元[J].計(jì)算力學(xué)學(xué)報(bào),2000,17(3):287-292,300.

[9] Bathe K J, Cimento A P. Some practical procedures for the solution of nonlinear finite element equations[J].Computer Methods in Applied Mechanics and Engineering, 1980,22(1):59-85.