一種高安全性顫振邊界預測方法

鐘華壽，張偉偉，肖華，葉正寅

(1.西北工業大學 航空學院，陜西 西安710072;2.中國飛行試驗研究院 中航工業飛行仿真航空科技重點實驗室，陜西 西安710089)

0 引言

顫振是由于彈性結構在空氣動力、慣性力和彈性力的耦合作用下，產生的一種具有破壞性的氣動彈性不穩定狀態，整個飛行包線內都不能出現顫振。當前顫振研究主要有理論計算、風洞試驗和顫振試飛三種手段。作為在飛行包線內對這些結果的最終驗證手段，許多新型號的顫振飛行試驗必不可少。

顫振試飛需要花費巨大的人力和物力，同時具有一定的風險，提出一種高效的、低成本、低風險的顫振預測方法對減少顫振試飛的周期和費用，降低顫振試飛風險具有重大的意義。Dimitriadis等［1］將目前主要的顫振預測方法分為兩類:第一類主要是辨識結構運動方程，如NG法［2］，其采用系統辨識方法，通過兩個不同速度下的響應，獲得運動方程的系數，然后求解不同速度下已被辨識的系統獲得顫振臨界速度;第二類主要是基于穩定性判據的曲線外推，包括阻尼外推法［3］、顫振余度法［4］、包線函數法［5］以及 ARMA 方法［6］等，其中，阻尼外推法和顫振余度法在國內型號試飛中使用較多［7］。以上兩類方法都無法通過一個亞臨界速度點預測顫振臨界特性，第二類方法甚至需要多個點，由此給顫振試飛工作帶來了額外的工作量。當亞臨界響應動壓逐漸遠離臨界動壓時，這兩類方法的顫振預測精度都有所降低，如阻尼外推法要獲得較高的顫振臨界點預測精度，要求使試飛點盡量接近顫振臨界點，由此給顫振試飛工作帶來了一定的風險性。

本文研究之前，文獻［8-9］中曾采用相似的研究方法對具有沉浮和俯仰兩自由度機翼的顫振問題進行了分析和研究。本文在文獻［8］的基礎上進行改進和推廣，進一步對具有多模態的三維彈性機翼展開分析和研究，即利用三維彈性機翼的一個亞臨界響應，采用系統辨識方法獲得該馬赫數下的非定常氣動力模型，而該模型與高度(動壓)無關，接著耦合結構狀態方程和氣動力模型，建立閉環系統的氣動彈性穩定性分析方程，該方程中的動壓以參數形式出現，通過分析系統穩定性隨動壓的變化規律，求解系統的顫振臨界特性。

1 研究方法

1.1 獲得結構模態參數

通過地面模態試驗或有限元計算獲得模態振型、模態質量和頻率等相關參數。

1.2 測試模態位移響應

設Δt為采樣時間間隔，利用加速度與位移的二階中心差分關系，解算出模態位移響應為:

1.3 模態氣動力響應的求解

求解模態坐標下的結構運動方程，可獲得模態氣動力系數響應Fa。

式中:M為廣義質量矩陣;G為廣義阻尼矩陣，忽略;K為廣義剛度矩陣;q為來流動壓;Fb為模態空間的外激力，本文Fb=0。

1.4 非定常氣動力建模

若結構振動位移較小，流場動力學特性主要表現為線性動態特征。選用離散型輸入輸出差分模型進行亞臨界響應的氣動力建模:

式中:矩陣Ai和Bi為待辨識的矩陣，以模態位移作為輸入，模態氣動力系數作為輸出，采用最小二乘法進行估計。

將式(4)差分模型轉化為狀態空間模型，轉化過程和各參數含義可詳見文獻［10-11］:

1.5 獲得臨界顫振特性

其中:

耦合結構狀態方程式(5)和氣動狀態方程式(6)得到如下氣動彈性系統的穩定性分析狀態方程:

將氣動彈性穩定性分析轉化為求解狀態矩陣的特征值后，通過求解不同動壓下狀態矩陣的特征值，繪出系統阻尼和頻率隨動壓變化的V-g圖和V-ω圖，從而判讀出系統的顫振臨界特性。

2 算例與分析

本文采用數值方法［12］獲得有關數據，對所發展的顫振預測方法的有效性進行研究和分析。

以AGARD 445.6機翼為例，取前三階模態進行分析。通過改變同一馬赫數下的來流動壓(本文通過改變來流密度)，分別對馬赫數0.499，0.678，0.901，0.960，1.072進行臨界動壓的尋找，求解時間步長為4×10－3s，顫振速度和顫振頻率以無量綱化形式給出，機翼詳細參數及無量綱化過程詳見文獻［13］。圖1為顫振速度邊界的試驗值和計算值的對比。由圖可以看出，所使用的程序和網格計算的顫振速度與試驗結果基本吻合。

圖1 計算值和試驗值的對比Fig.1 Comparison between calculated and experimental values

在Ma=0.499，q∞=0.85q*時(q*表示對應馬赫數下的顫振臨界動壓，下同)，對每階模態給定一個初始擾動，模態位移響應如圖2所示，求解出的模態氣動力系數響應如圖3所示。

圖2 模態位移響應Fig.2 Modal displacement responses

圖3 模態氣動力系數響應Fig.3 Modal aerodynamic coefficient responses

經辨識后，耦合結構和氣動狀態方程，建立閉環系統的穩定性分析方程，通過特征分析獲得V-g圖和V-ω圖，如圖4所示。根據頻率靠近原則，可判別出第一階模態和第二階模態發生耦合，第一階模態首先失穩，圖中的無量綱顫振速度為0.439 1，頻率比為0.574 9。對前三階位移響應信號施加白噪聲，信噪比為20 dB，采用帶通濾波，濾波前后的響應對比如圖5所示。經濾波后分析得到V-g圖和V-ω圖，如圖6所示。其無因次顫振速度為0.440 7，頻率比為0.582 5，與之前的預測結果極為接近。

圖4 通過特征分析獲得的V-g圖和V-ω圖Fig.4 V-g and V-ω diagrams obtained by eigen-analysis

圖5 含有20 dB噪聲的前三階模態位移濾波前后響應對比Fig.5 Comparison of modal displacement responses of the 1st three order mode with 20 dB noise before and after filtering

圖6 經加噪濾波后分析得到的V-g圖和V-ω圖Fig.6 V-g and V-ω diagrams by using the signals with noise

重復上述方法和過程，計算Ma=0.678，0.901，0.960，1.072 對應的q∞(依次為0.875q*，0.857q*，0.85q*，0.82q*)的響應，加上 Ma=0.499，q∞=0.85q*的結果，5個馬赫數點的預測結果如圖7所示。由圖可以看出，跨聲速區出現明顯的凹坑，預測效果較為理想。

圖7 預測值和計算值的對比Fig.7 Comparison between predicted and calculated values

在Ma=0.678時，依次降低來流動壓，分別計算0.875q*，0.750q*，0.625q*，0.500q*下的亞臨界響應，以研究亞臨界點逐漸遠離顫振臨界點時對預測精度的影響。以0.500q*為例，分析得到的V-g圖和V-ω圖，如圖8所示。各動壓比下計算出的無因次顫振臨界速度和頻率比如圖9所示。從圖9可以看出，隨著亞臨界點逐漸遠離顫振臨界點，顫振臨界速度和頻率比的預測精度都有所下降，但下降的值較小。從圖9中還可以看出，對于本算例，在Ma=0.678時，在0.625q*～q*范圍內，所預測的顫振臨界速度和頻率比隨動壓比近似呈線性變化，并且其斜率相當小，即在距顫振臨界點一定范圍內，利用亞臨界響應點對顫振臨界點預測的結果精度較高。當來流動壓小于0.625q*時，預測的誤差增大得較快;但當來流動壓為0.500q*時，其預測的無因次顫振速度與計算值之間的誤差也僅為7.62%。

圖8 通過特征分析獲得的V-g圖和V-ω圖Fig.8 V-g and V-ω diagrams obtained by eigen-analysis

圖9 不同動壓比下響應預測的和Fig.9 Predicted values of and at different dynamicpressure ratios

3 結論

(1)該方法只需利用三維彈性機翼一個速度較低的亞臨界響應即可獲得對應馬赫數下的顫振臨界特性，適用于兩個及以上復雜模態的顫振分析，可極大降低顫振試飛或風洞試驗的成本和風險;

(2)當響應測試動壓與顫振臨界動壓的動壓比降低時，顫振臨界特性的預測精度有所下降，但即使動壓比在0.5時，對于本算例，所預測無因次顫振臨界速度的誤差仍可控制在8%以內;

(3)可利用兩個或多個不同動壓下的亞臨界響應預測和驗證顫振臨界特性;

(4)初步研究只是在仿真環境中進行，還需要在試驗環境中，噪聲強度較大的情況下施加適當的激勵后驗證其實用性。

［1］ Dimitriadis G，Cooper J E.Flutter prediction from flight flutter test data［J］.Journal of Aircraft，2001，38(2):355-367.

［2］ Nissim E，Gilyard G B.Method for experimental determination of flutter speed by parameter identification［R］.AIAA Paper 89-1324，1989.

［3］ Rhlin C L，Watson J J，Ricketts R H，et al.Evaluation of four subcritical response methods for on-line prediction of flutter onset in wind tunnel tests［J］.Journal of Aircraft，1983，20(10):835-840.

［4］ Zimmerman N H，Weissenburger J T.Prediction of flutter onset speed based on flight testing at subcritical speeds［J］.Journal of Aircraft，1964，1(4):190-202.

［5］ Cooper J E，Emmet P R，Wright J，et al.Envelope function:a tool for analyzing flutter data［J］.Journal of Aircraft，1993，30(5):785-790.

［6］ Matsuzaki Y，Ando Y.Estimation of flutter boundary from random responses due to turbulence at subcritical speeds［J］.Journal of Aircraft，1981，18(10):826-868.

［7］ 俱利鋒，梁坤鵬，梁海洲，等.顫振邊界預測技術在顫振試飛中的應用研究［J］.飛行力學，2010，28(5):79-83.

［8］ 張偉偉，于俊杰，全景閣，等.一種基于亞臨界響應的顫振邊界預測新方法［J］.航空工程進展，2012，3(4):390-396.

［9］ Jieun S，Taehyoun K，Seung J.Experimental determination of unsteady aerodynamic coefficients and flutter behavior of a rigid wing［J］.Journal of Fluid and Structures，2012(29):50-61.

［10］ Zhang W W，Ye Z Y.Reduced-order-model-based flutter analysis at high angle of attack［J］.Journal of Aircraft，2007，44(6):2086-2089.

［11］ Zhang W W，Ye Z Y.Effect of control surface on airfoil flutter in transonic flow ［J］.Aeta Astronautica，2010，66(7-8):999-1007.

［12］張偉偉.基于CFD技術的高效氣動彈性分析方法研究［D］.西安:西北工業大學，2006.

［13］盧學成，葉正寅，張陳安.基于ANSYS/CFX耦合的機翼顫振分析［J］.計算機仿真，2010，27(9):88-91.