一種高安全性顫振邊界預(yù)測方法

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

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

0 引言

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

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

本文研究之前，文獻(xiàn)［8-9］中曾采用相似的研究方法對具有沉浮和俯仰兩自由度機(jī)翼的顫振問題進(jìn)行了分析和研究。本文在文獻(xiàn)［8］的基礎(chǔ)上進(jìn)行改進(jìn)和推廣，進(jìn)一步對具有多模態(tài)的三維彈性機(jī)翼展開分析和研究，即利用三維彈性機(jī)翼的一個亞臨界響應(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)試驗(yàn)或有限元計(jì)算獲得模態(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)運(yùn)動方程，可獲得模態(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)行估計(jì)。

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

圖1 計(jì)算值和試驗(yàn)值的對比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ù)上述方法和過程，計(jì)算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ù)點(diǎn)的預(yù)測結(jié)果如圖7所示。由圖可以看出，跨聲速區(qū)出現(xiàn)明顯的凹坑，預(yù)測效果較為理想。

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

在Ma=0.678時，依次降低來流動壓，分別計(jì)算0.875q*，0.750q*，0.625q*，0.500q*下的亞臨界響應(yīng)，以研究亞臨界點(diǎn)逐漸遠(yuǎn)離顫振臨界點(diǎn)時對預(yù)測精度的影響。以0.500q*為例，分析得到的V-g圖和V-ω圖，如圖8所示。各動壓比下計(jì)算出的無因次顫振臨界速度和頻率比如圖9所示。從圖9可以看出，隨著亞臨界點(diǎn)逐漸遠(yuǎn)離顫振臨界點(diǎn)，顫振臨界速度和頻率比的預(yù)測精度都有所下降，但下降的值較小。從圖9中還可以看出，對于本算例，在Ma=0.678時，在0.625q*～q*范圍內(nèi)，所預(yù)測的顫振臨界速度和頻率比隨動壓比近似呈線性變化，并且其斜率相當(dāng)小，即在距顫振臨界點(diǎn)一定范圍內(nèi)，利用亞臨界響應(yīng)點(diǎn)對顫振臨界點(diǎn)預(yù)測的結(jié)果精度較高。當(dāng)來流動壓小于0.625q*時，預(yù)測的誤差增大得較快;但當(dāng)來流動壓為0.500q*時，其預(yù)測的無因次顫振速度與計(jì)算值之間的誤差也僅為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)該方法只需利用三維彈性機(jī)翼一個速度較低的亞臨界響應(yīng)即可獲得對應(yīng)馬赫數(shù)下的顫振臨界特性，適用于兩個及以上復(fù)雜模態(tài)的顫振分析，可極大降低顫振試飛或風(fēng)洞試驗(yàn)的成本和風(fēng)險(xiǎn);

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

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

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

［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ī)翼顫振分析［J］.計(jì)算機(jī)仿真，2010，27(9):88-91.