基于多項式修正片條氣動力的跨音速顫振分析方法及其試驗驗證

2022-01-27 14:14:20張婷婷周健斌竇忠謙章俊杰

振動與沖擊 2022年1期

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

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

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

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

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

1 方法介紹

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

MSC.NASTRAN軟件采用偶極子格網法(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 氣動力修正原則

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

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

1.2 氣動力修正對象

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

1.3 氣動力修正目標

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

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

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

對每個氣動網格區域內的CP進行積分，得到對應于每個氣動網格的升力系數CLi、壓心位置Xi

(1)

(2)

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

L=∑(CLi·Si)

(3)

(4)

M=L·(X-X0)

(5)

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

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

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

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

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

為選出與目標氣動力分布符合性最好的最優解，繪制出各片條內升力和力矩的分布曲線，在求解區間中調整參數，選取升力和力矩分布曲線與目標曲線趨勢相同，與目標差值均方差最小的修正系數作為最優解。修正后的片條內弦向網格的升力斜率和力矩斜率分布示意圖如圖4、圖5所示。

M=ax2+bx+c

(6)

(7)

(8)

式中：a、b、c為多項式函數的常數項；x代表片條內弦向網格位置。

圖4 片條內弦向網格升力斜率分布示意圖Fig.4 Lift slope distribution of tangential grid

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

圖5 力片條內弦向網格力矩斜率分布示意圖Fig.5 Moment slope distribution of tangential grid

(9)

式中，Lij為確定各片條的a、b、c常數項后，得到的片條內各網格修正后的升力。

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

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

(10)

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

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

(11)

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

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

2 風洞試驗驗證

2.1 翼吊發動機機翼構型

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

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

圖6 機翼跨音速顫振模型風洞試驗照片Fig.6 The wing flutter model in the wind tunnel

圖7 機翼氣動網格示意圖Fig.7 FEM model of the wing structure

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

2.2 帶操縱面的尾翼構型

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