非定常氣動力計算與顫振分析

蘇 澤 楊紅杰

（中航西安飛機工業集團股份有限公司西飛設計院，陜西 西安 710089）

在設計現代飛行器的過程中，必須要考慮氣動彈性的影響。傳統的氣動彈性問題主要指發散、舵面反效代表的靜氣動彈性問題與顫振代表的動氣動彈性問題。這類問題一般直接在頻域中求解氣動彈性方程，對亞音速飛行器來說，目前工程中普遍采用有限元方法直接進行模態分析，再利用偶極子格網法（DLM）求得頻域離散形式的非定常氣動力系數矩陣，從而實現氣動力與結構模態間的直接耦合。

當飛行器包括控制系統時，必須要著重考慮控制系統與彈性機體結構之間的耦合作用，也就是氣動伺服彈性力學（ASE）。傳統飛控系統通常采用SISO 控制方式，工程中仍采用頻域分析方法，基于經典控制理論進行穩定性分析。而現代飛行器控制回路之間的耦合作用明顯，其飛控系統模型也具有MIMO 的特征，由此產生了以狀態空間方法為基礎的現代控制理論。此時，通過有理函數擬合方法將頻域離散矩陣形式的非定常氣動力轉化為有理函數形式，使其能夠進行拉氏反變換，從而建立狀態方程，這是建立氣動伺服彈性狀態空間模型的關鍵[1]。

1 非定常氣動力

通常情況下，氣動彈性問題中的氣動力主要是指飛行器彈性模態、舵面剛體偏轉模態等廣義位移帶來的氣動下洗產生的非定常氣動力。

飛行器氣動彈性一般運動方程如公式（1）所示。

式中：q為飛行器彈性模態坐標向量；δ為控制面剛體偏轉坐標向量；ρ為大氣密度；V為飛行速度；Mqq和Mqδ分別為飛行器彈性模態和舵面剛體偏轉模態對應廣義質量矩陣；Cqq為廣義阻尼矩陣；Kqq為廣義剛度矩陣；Qqq和Qqδ分別為彈性模態和舵面剛體偏轉模態產生的廣義非定常氣動力系數矩陣。

對亞音速飛行器來說，工程中一般采用偶極子格網法來計算頻域離散矩陣形式的非定常氣動力。偶極子格網法是一種基于小擾動線化位勢流方程的面元法，當計算非定常氣動力時，需要合理地對氣動面進行網格劃分，將氣動面沿展向和弦向分成若干個兩側邊平行于來流方向的梯形面元，一般在面元網格的1/4 弦線上布置壓力偶極子，在網格的3/4 處布置控制點，通過物面邊界條件計算控制點的下洗速度[2]。

根據廣義氣動力的定義，其矩陣表達形式如公式（2）所示。

式中：R為廣義氣動力；NP為網格控制點處的模態矩陣，可以通過氣動控制點與結構節點之間的插值矩陣，由結構模態矩陣轉換得到；S為網格面積的加權矩陣，其對角項為各氣動網格的面積；?p為氣動面元網格的壓力分布。

根據非定常氣動力理論，根據網格控制點滿足的積分方程可以得到公式（3）。

式中：wq和wδ分別為彈性模態和舵面偏轉模態帶來的氣動網格控制點的下洗角；D-1為氣動力影響系數矩陣，也稱AIC矩陣。

當計算飛行器彈性模態與舵面剛體偏轉模態引起的非定常氣動力系數時，氣動網格控制點的下洗角與廣義坐標之間滿足公式（4）。

式中：b為參考半弦長；k為減縮頻率，k=ωb/V；Nq和Nδ分別為網格控制點處的彈性模態和舵面偏轉模態。

將公式（2）、公式（3）代入公式（4），得到廣義非定常氣動力的表達式，如公式（5）所示。

當給定氣動模型的網格劃分、飛行馬赫數以及減縮頻率時，可以確定飛行器彈性模態、舵面剛體偏轉模態等帶來的非定常氣動力。

2 有理函數擬合

當采用時域方法分析氣動伺服彈性問題和一些非線性問題時，要先對氣動彈性廣義運動方程進行拉氏變換，再將其轉換為狀態空間方程形式。而氣彈運動方程中采用偶極子格網法求得的離散矩陣形式的頻域氣動力是減縮頻率的復函數，無法直接進行拉氏反變換。此時，便需要引入有理函數擬合方法，在不損失計算精度的前提下，延拓頻域氣動力向拉氏域。

工程中最常用的有理函數近似方法包括ROGER 法、修正矩陣法（MMP）和最小狀態法（MS）等，這些方法均以最小二乘法為基礎。當采用這些方法擬合非定常氣動力建立狀態空間方程時，均需要引入氣動力擴充項作為狀態變量。不同方法引入的氣動力擴充項數量也不同，ROGER 法產生的氣動擴充項數量為結構模態數與氣動滯后根的乘積，MMP法對應的氣動力擴充項數量為廣義氣動力系數矩陣各列對應的氣動滯后根數量之和，而MS 法對應的氣動力擴充項數等于氣動滯后根的數量[4]。研究表明，當氣動力擴充項的數量相同時，MS 法的擬合精度最高，下面給出MS 法的主要擬合過程。

當給定飛行馬赫數時，頻域廣義氣動力系數矩陣可以為減縮頻率的復函數，統一用Qq（k）表示。采用MS 法進行有理函數擬合的表達式如公式（6）所示。

將公式（6）展開為頻域內實部和虛部的形式，如公式（7）所示。

當擬合精度較高時，選定的一系列減縮頻率下的氣動力擬合結果應與頻域計算結果近似相等。為了獲得更精確的有理函數擬合結果，需要假設在若干減縮頻率點處的擬合結果與頻域計算結果精確相等，即假設一定的約束條件：1）當約束減縮頻率為0 時，氣動力擬合結果與頻域氣動力相等。2）約束在減縮頻率k1處氣動力實部擬合結果與頻域氣動力實部相等。3）約束在減縮頻率k2處氣動力虛部擬合結果與頻域氣動力虛部相等。

由上述約束條件及其他約束條件可以將擬合公式中A0、A1和A2分別轉化為減縮頻率k、D和E矩陣的表達式，以進一步給出其他非約束減縮頻率點處的氣動力系數擬合公式。

為了進一步求得擬合公式的解，MS 法需要先給定R的矩陣元素，再由最小二乘法確定矩陣D和E。首先，給定矩陣E，按行擬合出矩陣D。其次，由現有的矩陣R和D，按列擬合求出矩陣E。最后，計算擬合的精度，如果擬合的精度不滿足要求，就重復前面的擬合過程，反復迭代計算D-E-D，直到得到滿意的擬合結果。一般情況下，迭代10次即可收斂。

當求出各系數矩陣后，令s=ik（s為拉普拉斯算子），將減縮頻率轉化為拉氏變量，氣動力擬合如公式（8）所示。

由此，利用有理函數擬合方法就可以將廣義非定常氣動力從坐標軸拓展到整個拉氏域。

3 算例

該文以一個典型翼面為例對上述非定常氣動力計算方法和有理函數擬合方法進行驗證。

3.1 計算模型

計算模型由機翼翼面與機翼根部連接組成，翼面采用板單元建模，根部連接為梁單元，并在梁單元根部固支約束。利用MSC.NASTRAN 軟件進行模態分析，得出前五階結構彈性模態見表1。

表1 模態分析結果

采用ZAERO 軟件建立機翼翼面的氣動網格，為了簡化步驟，沿弦向和展向均等分為5 段。為了驗證有理函數擬合方法的精確性，分別采用頻域法和狀態空間方法對上述模型進行開環顫振分析，顫振計算方法選取ZAERO 軟件中的g法，固定高度為海平面，馬赫數0.05，計算選取前五階彈性模態，不考慮阻尼，取減縮頻率為0.0～1.0。

當狀態空間法計算時，采用MS 法進行廣義氣動力系數矩陣的有理函數擬合，當約束減縮頻率為0 時，氣動力擬合結果與頻域氣動力相等，取擬合迭代次數為100 次，并根據經驗公式法選取10 個氣動滯后根，如公式（9）所示。

式中：Ri為第i個氣動滯后根取值；kmax為計算選取的最大減縮頻率；Nlag為氣動滯后根數量。

3.2 顫振結果對比

采用ZAERO 軟件中的g 法進行顫振分析，得到頻域法V-G、V-F 曲線，如圖1所示。2 種方法顫振頻率與顫振速度對比見表2。

由圖1 可知，計算得到一支顫振，其類型為典型彎扭耦合顫振，參與模態為機翼一階彎曲和一階扭轉模態，顫振頻率為8.112 Hz，顫振速度較高。在采用MS 法進行有理函數擬合后，顫振頻率及顫振速度與頻域計算結果偏差較小，說明氣動力擬合過程并沒有降低計算精度。

3.3 非定常氣動力擬合對比

頻域法計算不采用有理函數擬合，其得到的非定常氣動力系數為精確解，可為n×m維離散矩陣，每項矩陣元素均包括實部和虛部2 個部分（n為計算采用的彈性模態數量；m為彈性模態數量與舵面剛體偏轉模態數量的和）。該算例中不包括操縱面，機體彈性模態數量為5，因此頻域法計算得到的非定常氣動力系數矩陣為5×5 的復數矩陣。

第3.2 節得出顫振的主要參與模態為翼面的第一階模態、二階模態，不同減縮頻率的2 種方法得到的廣義氣動力系數矩陣元素Q11與Q22的對比如圖2所示，下標表示元素在矩陣中的位置。

圖2 矩陣元素隨減縮頻率變化趨勢對比

對比氣動力系數矩陣元素可以得出，在給定的減縮頻率下，采用MS 法進行非定常氣動力有理函數擬合，當選擇合適數量的氣動滯后根時，擬合得到的廣義氣動力系數矩陣元素與頻域直接計算值已基本一致，氣動力隨減縮頻率的變化規律也與頻域法一致，擬合過程整體偏差較小，擬合效果較好，結論與顫振對比結論一致。

4 結語

該文論述了亞音速飛行器非定常氣動力建模方法，給出了基于MS 法進行廣義氣動力有理函數擬合的計算方法，并根據典型翼面算例分析了頻域氣動力、擬合后的氣動力對顫振結果的影響。分析結果表明，當采用MS 法進行非定常氣動力擬合時，選擇合適數量的氣動滯后根后能夠得到較準確的氣動力計算結果，其擬合精度能夠滿足顫振分析需求。

實際上，針對線性系統的顫振分析、頻域ASE 分析、機動載荷分析以及離散陣風分析等傳統氣動彈性問題，工程中一般直接采用頻域計算得到的非定常氣動力；而針對時域ASE 分析及其他非線性動力學響應分析，例如非線性顫振、投放載荷分析等，難以直接使用頻域氣動力，此時，必須要進行有理函數擬合，通過狀態空間方法得到可靠結果。