張婷婷， 周健斌， 竇忠謙， 章俊杰

(上海飛機(jī)設(shè)計(jì)研究院，上海 201210)

飛機(jī)的跨音速顫振特性研究是民機(jī)顫振設(shè)計(jì)的重要部分。當(dāng)前，國(guó)內(nèi)民機(jī)的顫振設(shè)計(jì)工作主要是通過(guò)跨音速顫振模型風(fēng)洞試驗(yàn)來(lái)獲得飛機(jī)部件的跨音速壓縮性曲線[1-4]。在分析方面，由于跨音速區(qū)的氣流存在激波、流動(dòng)分離等問(wèn)題，精確的跨音速顫振計(jì)算(尤其在高馬赫數(shù)區(qū))存在很大的難度。

近年來(lái)，隨著計(jì)算流體力學(xué)技術(shù)的發(fā)展，基于歐拉方程/N-S方程以及結(jié)構(gòu)動(dòng)力學(xué)方程的時(shí)域顫振分析方法取得了一定進(jìn)展[5-8]，但是其計(jì)算量大、耗時(shí)長(zhǎng)。基于小擾動(dòng)線性諧振蕩假設(shè)的面元法，計(jì)算效率快、計(jì)算精度高，在工程顫振分析中發(fā)揮著重要作用，由于其不能考慮攻角、翼型的影響，在跨音速顫振計(jì)算時(shí)，通常需要先對(duì)氣動(dòng)力進(jìn)行一輪修正。近年來(lái)，國(guó)內(nèi)學(xué)者對(duì)面元法的氣動(dòng)力修正方法進(jìn)行了研究探索[9-11]，但大多用于理論分析，在民機(jī)顫振設(shè)計(jì)工程中少有應(yīng)用。

本文以MSC.NASTRAN軟件為平臺(tái)，提出一種適用于工程應(yīng)用的，基于多項(xiàng)式的氣動(dòng)力修正方法，以片條內(nèi)升力和力矩隨攻角變化斜率為修正目標(biāo)，用一組多項(xiàng)式方程模擬片條力矩分布，保證整個(gè)翼面的氣動(dòng)力大小和分布都與目標(biāo)相符，進(jìn)而使用修正后的氣動(dòng)力進(jìn)行跨音速區(qū)的顫振分析。該方法適宜于工程使用，快速有效，且已在某型民機(jī)的跨音速顫振模型風(fēng)洞試驗(yàn)中得到了應(yīng)用和驗(yàn)證，計(jì)算結(jié)果的跨音速壓縮性曲線趨勢(shì)與試驗(yàn)一致，凹坑馬赫數(shù)一致。

1 方法介紹

面元法無(wú)法考慮翼面厚度、非線性影響等因素，在跨音速區(qū)計(jì)算的氣動(dòng)力與實(shí)際有一定偏差。并且由于非定常氣動(dòng)力計(jì)算機(jī)理復(fù)雜，影響因素眾多，計(jì)算精度通常難以保證，直接修正非定常氣動(dòng)力數(shù)據(jù)存在較大難度。使用計(jì)算流體力學(xué)(computational fluid dynamics, CFD)可以求解非線性的流場(chǎng)，得到準(zhǔn)確的定常氣動(dòng)力數(shù)據(jù)，對(duì)面元法的定常氣動(dòng)力部分進(jìn)行修正，得到氣動(dòng)力修正系數(shù)矩陣Wkk，然后將Wkk應(yīng)用于所有頻率段中的非定常氣動(dòng)力修正，使用修正后的非定常氣動(dòng)力進(jìn)行顫振計(jì)算，可有效提高跨音速區(qū)顫振分析的準(zhǔn)確性。

MSC.NASTRAN軟件采用偶極子格網(wǎng)法(doublet-lattice method,DLM)計(jì)算氣動(dòng)力，本文使用CFD計(jì)算的定常氣動(dòng)力去修正DLM計(jì)算的定常氣動(dòng)力，得到Wkk，然后將Wkk應(yīng)用于所有頻率段中的非定常氣動(dòng)力修正，使用修正后的非定常氣動(dòng)力進(jìn)行顫振計(jì)算。修正氣動(dòng)力的跨音速顫振計(jì)算流程圖，如圖1所示。

圖1 氣動(dòng)力修正的跨音速顫振計(jì)算流程圖Fig.1 Flow chart of transonic flutter calculation based on aerodynamic force modification

1.1 氣動(dòng)力修正原則

為保證修正后的氣動(dòng)力與目標(biāo)盡可能相似，需要確保力的大小和分布均與目標(biāo)相同，即力和力矩與目標(biāo)相同。如果選擇以網(wǎng)格為單位進(jìn)行氣動(dòng)力修正，保證各網(wǎng)格氣動(dòng)力與目標(biāo)相同，由于DLM各網(wǎng)格的壓心位置是固定的(1/4弦長(zhǎng)處)，而CFD各網(wǎng)格的壓心位置都不同，網(wǎng)格的力矩?zé)o法滿足目標(biāo)要求。因此，將翼面沿展向劃分為若干片條，以片條為單位，通過(guò)協(xié)調(diào)各網(wǎng)格的氣動(dòng)力大小，來(lái)保證整個(gè)片條的力和力矩與目標(biāo)相同。

由于CFD計(jì)算的力和力矩的零升攻角不同，對(duì)應(yīng)于DLM的某一給定攻角，CFD的力和力矩處于不同的攻角狀態(tài),無(wú)法以某一攻角狀態(tài)的力和力矩為目標(biāo)來(lái)做修正。而在小攻角范圍內(nèi)，升力和力矩隨攻角變化近似為線性的，可以選取該線性變化的斜率作為修正目標(biāo)。同時(shí)，為確保力和力矩的斜率大于0，修正系數(shù)也必須大于0。

1.2 氣動(dòng)力修正對(duì)象

DLM的氣動(dòng)力和力矩隨攻角的變化曲線為過(guò)原點(diǎn)的直線，因此，1°攻角時(shí)的氣動(dòng)力和力矩即為升力和力矩隨攻角變化的斜率。同時(shí)，DLM網(wǎng)格的壓心位置是確定的，所以力矩分布是由氣動(dòng)力的分布決定的。因此，選取1°攻角狀態(tài)的各網(wǎng)格氣動(dòng)力作為修正對(duì)象。應(yīng)用MSC.NASTRAN軟件的靜氣動(dòng)彈性分析模塊，計(jì)算給定馬赫數(shù)下1°攻角狀態(tài)翼面上各網(wǎng)格的氣動(dòng)力L0。

1.3 氣動(dòng)力修正目標(biāo)

建立CFD模型，計(jì)算翼面在給定馬赫數(shù)和攻角狀態(tài)(選取小角度范圍-5°～5°)下的定常氣動(dòng)力，得到各片條站位下的壓力系數(shù)分布曲線，如圖2、圖3所示。

圖2 CFD計(jì)算模型示意圖Fig.2 The CFD model

圖3 單片條內(nèi)壓力曲線分布示意圖Fig.3 Pressure curve distribution in a single strip

對(duì)每個(gè)氣動(dòng)網(wǎng)格區(qū)域內(nèi)的CP進(jìn)行積分，得到對(duì)應(yīng)于每個(gè)氣動(dòng)網(wǎng)格的升力系數(shù)CLi、壓心位置Xi

(1)

(2)

再計(jì)算出各片條的升力L和對(duì)前緣點(diǎn)的力矩M

L=∑(CLi·Si)

(3)

(4)

M=L·(X-X0)

(5)

式中:Si為第i個(gè)計(jì)算網(wǎng)格的面積;X為片條壓心坐標(biāo);Xi為第i個(gè)網(wǎng)格的壓心坐標(biāo);X0為網(wǎng)格前緣點(diǎn)坐標(biāo)。

繪制片條升力和力矩隨攻角變化曲線，取小角度范圍內(nèi)的線性段斜率，作為氣動(dòng)力修正要滿足的目標(biāo)。

1.4 基于多項(xiàng)式的氣動(dòng)力修正方法

一個(gè)片條內(nèi)弦向通常劃分有多個(gè)網(wǎng)格，即修正系數(shù)為多個(gè)未知量，需同時(shí)滿足兩個(gè)等式方程和一組不等式方程，為多解問(wèn)題，求解困難。

本文提出一種基于多項(xiàng)式的氣動(dòng)力修正方法，用一組多項(xiàng)式方程來(lái)描述片條內(nèi)弦向網(wǎng)格的力矩分布，各片條的多項(xiàng)式方程根據(jù)片條力矩具體分布情況進(jìn)行調(diào)整，目的是將未知量降為3個(gè)，求解后可得到一組求解區(qū)間，在區(qū)間內(nèi)取值，都滿足力和力矩與修正目標(biāo)相同的等式方程。

為選出與目標(biāo)氣動(dòng)力分布符合性最好的最優(yōu)解，繪制出各片條內(nèi)升力和力矩的分布曲線，在求解區(qū)間中調(diào)整參數(shù)，選取升力和力矩分布曲線與目標(biāo)曲線趨勢(shì)相同，與目標(biāo)差值均方差最小的修正系數(shù)作為最優(yōu)解。修正后的片條內(nèi)弦向網(wǎng)格的升力斜率和力矩斜率分布示意圖如圖4、圖5所示。

M=ax2+bx+c

(6)

(7)

(8)

式中：a、b、c為多項(xiàng)式函數(shù)的常數(shù)項(xiàng)；x代表片條內(nèi)弦向網(wǎng)格位置。

圖4 片條內(nèi)弦向網(wǎng)格升力斜率分布示意圖Fig.4 Lift slope distribution of tangential grid

氣動(dòng)力修正系數(shù)矩陣Wkk的目的是使偶極子格網(wǎng)法計(jì)算的升力和力矩匹配CFD計(jì)算的升力和力矩，其元素Wij為翼面各網(wǎng)格的氣動(dòng)力修正系數(shù)。

圖5 力片條內(nèi)弦向網(wǎng)格力矩斜率分布示意圖Fig.5 Moment slope distribution of tangential grid

(9)

式中，Lij為確定各片條的a、b、c常數(shù)項(xiàng)后，得到的片條內(nèi)各網(wǎng)格修正后的升力。

1.5 氣動(dòng)力修正后的顫振分析

氣動(dòng)彈性一般運(yùn)動(dòng)方程寫為以下的形式

(10)

式中:q為廣義坐標(biāo)列陣;M為廣義質(zhì)量對(duì)角矩陣;K為廣義剛度對(duì)角矩陣;Q為廣義非定常氣動(dòng)力列陣。

將氣動(dòng)力修正系數(shù)矩陣Wkk用于非定常氣動(dòng)力的計(jì)算

(11)

式中:FP為網(wǎng)格氣作用點(diǎn)處的模態(tài)矩陣;S為面積加權(quán)陣，對(duì)角項(xiàng)為各氣動(dòng)網(wǎng)格的面積;D為非定常氣動(dòng)力影響系數(shù);FH為控制點(diǎn)的模態(tài)矩陣;k為減縮頻率;t為參考長(zhǎng)度。

使用修正后的廣義非定常氣動(dòng)力Q，用p-k法進(jìn)行顫振方程求解，得出顫振速度與顫振頻率。

2 風(fēng)洞試驗(yàn)驗(yàn)證

2.1 翼吊發(fā)動(dòng)機(jī)機(jī)翼構(gòu)型

某型民機(jī)機(jī)翼跨音速顫振風(fēng)洞試驗(yàn)?zāi)Ｐ蜑橐淼醢l(fā)動(dòng)機(jī)構(gòu)型的單機(jī)翼模型，主要部件包含主翼面、小翼、吊掛和發(fā)動(dòng)機(jī)。如圖6所示。試驗(yàn)風(fēng)洞截面為正方形，尺寸2.4 m×2.4 m，試驗(yàn)段長(zhǎng)7 m，是引射式、半回流、暫沖型跨聲速風(fēng)洞，風(fēng)洞可用馬赫數(shù)范圍為0.3～1.2。試驗(yàn)獲得了以下馬赫數(shù)的臨界顫振速度：0.60、0.70、0.75、0.80、0.82、0.85。

氣動(dòng)網(wǎng)格示意圖如圖7所示，使用本文中的分析方法對(duì)機(jī)翼翼面的氣動(dòng)力進(jìn)行修正，并進(jìn)行了顫振分析。將跨音速區(qū)的計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果進(jìn)行對(duì)比，結(jié)果如表1和圖8所示(顫振速度已經(jīng)過(guò)無(wú)量綱處理)。

圖6 機(jī)翼跨音速顫振模型風(fēng)洞試驗(yàn)照片F(xiàn)ig.6 The wing flutter model in the wind tunnel

圖7 機(jī)翼氣動(dòng)網(wǎng)格示意圖Fig.7 FEM model of the wing structure

從計(jì)算結(jié)果可以看出，在高馬赫數(shù)區(qū)，修正前的計(jì)算結(jié)果無(wú)法得到速度凹坑，通過(guò)氣動(dòng)力修正，計(jì)算結(jié)果較好地獲取了跨音速區(qū)顫振特性，速度隨馬赫數(shù)變化趨勢(shì)與試驗(yàn)曲線一致，且速度最低點(diǎn)馬赫數(shù)區(qū)域與試驗(yàn)結(jié)果一致，當(dāng)量顫振速度偏差在5%以內(nèi)，計(jì)算結(jié)果絕對(duì)值總體低于試驗(yàn)值，偏于保守。在低馬赫數(shù)區(qū)，氣動(dòng)力修正對(duì)顫振計(jì)算結(jié)果影響較小。

2.2 帶操縱面的尾翼構(gòu)型

某型民機(jī)平尾跨音速顫振風(fēng)洞試驗(yàn)?zāi)Ｐ蜑閹Р倏v面構(gòu)型的單平尾模型，主要部件包括水平安定面、升降舵、以及舵面連接機(jī)構(gòu)。如圖9所示。在風(fēng)洞試驗(yàn)前，應(yīng)用本文中的分析方法計(jì)算出顫振速度隨馬赫數(shù)變化曲線，以此為基礎(chǔ)，規(guī)劃吹風(fēng)試驗(yàn)點(diǎn),試驗(yàn)點(diǎn)規(guī)劃示意圖如圖10所示，相比傳統(tǒng)的逐步逼近亞臨界顫振點(diǎn)的吹風(fēng)方法，車次減少了近30%，試驗(yàn)成本和風(fēng)險(xiǎn)都顯著降低。

圖9 平尾跨音速顫振模型風(fēng)洞試驗(yàn)照片F(xiàn)ig.9 The horizontal tail flutter model in the wind tunnel

圖10 試驗(yàn)點(diǎn)規(guī)劃示意圖Fig.10 Test site planning

試驗(yàn)獲得了以下馬赫數(shù)的臨界顫振速度：0.60、0.70、0.75、0.78、0.80、0.82。顫振計(jì)算結(jié)果如表2所示(其中，顫振速度已經(jīng)過(guò)無(wú)量綱處理)。

表2 平尾模型顫振計(jì)算結(jié)果

從計(jì)算結(jié)果可以看出，在跨音速區(qū)，顫振速度隨馬赫數(shù)變化趨勢(shì)與試驗(yàn)一致，速度凹坑點(diǎn)馬赫數(shù)區(qū)域與試驗(yàn)一致，馬赫數(shù)在(0.60～0.75)區(qū)間內(nèi)，當(dāng)量顫振速度最大偏差小于3%，馬赫數(shù)在(0.78～0.82)區(qū)間內(nèi)，當(dāng)量顫振速度最大偏差在10%以內(nèi)，計(jì)算結(jié)果絕對(duì)值總體低于試驗(yàn)值，偏于保守。

3 結(jié) 論

本文提出了一種適用于工程的，基于多項(xiàng)式的氣動(dòng)力修正方法，準(zhǔn)確模擬了翼面的氣動(dòng)力分布，使用修正后的氣動(dòng)力進(jìn)行顫振計(jì)算，提高了跨音速區(qū)的顫振計(jì)算精度?？缫羲兕澱衲Ｐ惋L(fēng)洞試驗(yàn)表明：

(1) 該方法對(duì)翼吊發(fā)動(dòng)機(jī)構(gòu)型的機(jī)翼顫振型、帶操縱面的尾翼顫振型都有較高的精度，且分析結(jié)果偏保守。

(2) 計(jì)算獲得的跨音速壓縮性曲線和試驗(yàn)結(jié)果相比，曲線趨勢(shì)一致，凹坑位置一致。

本文主要針對(duì)大展弦比的翼面(如機(jī)翼、平尾)進(jìn)行了氣動(dòng)力修正，發(fā)動(dòng)機(jī)、機(jī)身的氣動(dòng)力修正方法及其對(duì)顫振分析結(jié)果的影響有待進(jìn)一步研究。