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

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

0 引言

顫振是由于彈性結(jié)構(gòu)在空氣動力、慣性力和彈性力的耦合作用下，產(chǎn)生的一種具有破壞性的氣動彈性不穩(wěn)定狀態(tài)，整個飛行包線內(nèi)都不能出現(xiàn)顫振。當(dāng)前顫振研究主要有理論計算、風(fēng)洞試驗和顫振試飛三種手段。作為在飛行包線內(nèi)對這些結(jié)果的最終驗證手段，許多新型號的顫振飛行試驗必不可少。

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

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

1 研究方法

1.1 獲得結(jié)構(gòu)模態(tài)參數(shù)

通過地面模態(tài)試驗或有限元計算獲得模態(tài)振型、模態(tài)質(zhì)量和頻率等相關(guān)參數(shù)。

1.2 測試模態(tài)位移響應(yīng)

設(shè)Δt為采樣時間間隔，利用加速度與位移的二階中心差分關(guān)系，解算出模態(tài)位移響應(yīng)為:

1.3 模態(tài)氣動力響應(yīng)的求解

求解模態(tài)坐標(biāo)下的結(jié)構(gòu)運動方程，可獲得模態(tài)氣動力系數(shù)響應(yīng)Fa。

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

1.4 非定常氣動力建模

若結(jié)構(gòu)振動位移較小，流場動力學(xué)特性主要表現(xiàn)為線性動態(tài)特征。選用離散型輸入輸出差分模型進(jìn)行亞臨界響應(yīng)的氣動力建模:

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

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

1.5 獲得臨界顫振特性

其中:

耦合結(jié)構(gòu)狀態(tài)方程式(5)和氣動狀態(tài)方程式(6)得到如下氣動彈性系統(tǒng)的穩(wěn)定性分析狀態(tài)方程:

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

2 算例與分析

本文采用數(shù)值方法［12］獲得有關(guān)數(shù)據(jù)，對所發(fā)展的顫振預(yù)測方法的有效性進(jìn)行研究和分析。

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

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

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

圖2 模態(tài)位移響應(yīng)Fig.2 Modal displacement responses

圖3 模態(tài)氣動力系數(shù)響應(yīng)Fig.3 Modal aerodynamic coefficient responses

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

(4)初步研究只是在仿真環(huán)境中進(jìn)行，還需要在試驗環(huán)境中，噪聲強度較大的情況下施加適當(dāng)?shù)募詈篁炞C其實用性。

［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］ 俱利鋒，梁坤鵬，梁海洲，等.顫振邊界預(yù)測技術(shù)在顫振試飛中的應(yīng)用研究［J］.飛行力學(xué)，2010，28(5):79-83.

［8］ 張偉偉，于俊杰，全景閣，等.一種基于亞臨界響應(yīng)的顫振邊界預(yù)測新方法［J］.航空工程進(jìn)展，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技術(shù)的高效氣動彈性分析方法研究［D］.西安:西北工業(yè)大學(xué)，2006.

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