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

(上海飛機設計研究院，上海 201210)

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

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

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

1 方法介紹

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

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

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

1.1 氣動力修正原則

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

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

1.2 氣動力修正對象

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

1.3 氣動力修正目標

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

圖2 CFD計算模型示意圖Fig.2 The CFD model

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

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

(1)

(2)

再計算出各片條的升力L和對前緣點的力矩M

L=∑(CLi·Si)

(3)

(4)

M=L·(X-X0)

(5)

式中:Si為第i個計算網(wǎng)格的面積;X為片條壓心坐標;Xi為第i個網(wǎng)格的壓心坐標;X0為網(wǎng)格前緣點坐標。

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

1.4 基于多項式的氣動力修正方法

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

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

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

M=ax2+bx+c

(6)

(7)

(8)

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

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

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

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

(9)

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

1.5 氣動力修正后的顫振分析

氣動彈性一般運動方程寫為以下的形式

(10)

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

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

(11)

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

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

2 風洞試驗驗證

2.1 翼吊發(fā)動機機翼構(gòu)型

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

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

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

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

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

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

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

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

圖10 試驗點規(guī)劃示意圖Fig.10 Test site planning

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

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

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

3 結(jié) 論

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

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

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

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