飛機(jī)機(jī)翼顫振特性測試仿真研究

于之靖 ，花 貞 ，王 爍 ，王 威 ，吳 軍 ，諸葛晶昌

（中國民航大學(xué)a.電子信息與自動化學(xué)院；b.航空工程學(xué)院，天津 300300）

機(jī)翼動態(tài)變形與顫振特性是飛機(jī)性能和安全的重點(diǎn)考慮因素，也是氣動彈性變形分析、飛行試驗(yàn)科目和顫振適航認(rèn)證的重點(diǎn)測試內(nèi)容之一[1-3]。因此需要獲得盡可能多的機(jī)翼變形和顫振測量數(shù)據(jù)，為飛機(jī)結(jié)構(gòu)設(shè)計和適航認(rèn)證測試及應(yīng)用提供足夠的試驗(yàn)分析數(shù)據(jù)，以確保飛機(jī)飛行安全。

機(jī)翼顫振是機(jī)翼動態(tài)變形的極端表現(xiàn)形式，它是機(jī)翼在氣流中受到氣動力、彈性力和慣性力的耦合作用而發(fā)生的一種動態(tài)氣動彈性現(xiàn)象，極端條件下將使飛機(jī)結(jié)構(gòu)在極短的時間內(nèi)遭到破壞，嚴(yán)重影響飛機(jī)飛行安全。從固定翼飛機(jī)問世開始，各種機(jī)翼變形和顫振問題一直備受關(guān)注[4-6]。

1 技術(shù)背景

隨著對機(jī)翼變形和顫振研究的逐漸深入，相關(guān)學(xué)者提出各種關(guān)于機(jī)翼適航顫振的理論。亞臨界試驗(yàn)技術(shù)的發(fā)展，使得各種適航測試技術(shù)得以迅速發(fā)展，完善了測試數(shù)據(jù)的來源。而顫振適航測試需要盡可能貼合真實(shí)的飛行環(huán)境來模擬測試機(jī)翼振動情況，針對性地測量激振點(diǎn)。顫振適航飛行試驗(yàn)可以準(zhǔn)確地確定適航顫振包線，但難以控制飛行環(huán)境，而且不能綜合分析顫振因素，從而激發(fā)了人們選用其他技術(shù)代替適航飛行試驗(yàn)。風(fēng)洞試驗(yàn)[7]基于相對運(yùn)動定理，創(chuàng)建地面條件模擬飛機(jī)在空中的飛行狀態(tài)，但考慮到模型和環(huán)境的局限性，試驗(yàn)結(jié)果容易失真。梁模擬機(jī)翼試驗(yàn)只能用于機(jī)翼顫振特性的簡單分析，精確度不高。有限元數(shù)值模擬可以計算機(jī)翼在理想或?qū)嶋H的環(huán)境下的各種狀態(tài)，是目前常用的機(jī)翼顫振特性分析手段。偶極子格網(wǎng)法（DLM）計算流場相對簡單，但不能模擬過于復(fù)雜的流場[8-9]。隨著計算流體動力學(xué)（CFD）技術(shù)的發(fā)展，流固耦合計算方法隨之得以深入研究[10-11]。

本文基于顫振適航測量的目的，采用流固耦合的計算方法分析飛機(jī)流場，進(jìn)而研究機(jī)翼顫振特性，可以有效提高機(jī)翼顫振適航測量激振點(diǎn)選取的針對性，減少測試工作量。

2 機(jī)翼振動特性理論分析

機(jī)翼的顫振形式通常表現(xiàn)為扭轉(zhuǎn)變形導(dǎo)致的結(jié)構(gòu)解體、應(yīng)力應(yīng)變過大而導(dǎo)致的結(jié)構(gòu)剛度不足，因此，機(jī)翼顫振的研究主要集中在其氣動彈性分析以及振動響應(yīng)分析。通常，機(jī)翼結(jié)構(gòu)的振動測試激振區(qū)域的選取以其應(yīng)力、應(yīng)變或者其動態(tài)變形的集中分布區(qū)域作為依據(jù)。機(jī)翼的顫振特性可用動力學(xué)分析來確定，為完整描述飛機(jī)在飛行過程中所遇到的流固耦合[12-14]問題，研究耦合流場壓力效應(yīng)作用下機(jī)翼的顫振特性規(guī)律，根據(jù)Ansys的求解理論，通過機(jī)翼靜態(tài)條件下的線性分析，將在耦合流場作用下求解得到的流體壓力邊界條件加載到機(jī)翼結(jié)構(gòu)的動力學(xué)平衡方程中，得到機(jī)翼耦合條件下的動力特性方程為

其中：M、C、K 分別為質(zhì)量、阻尼、剛度矩陣；Mf、R、Kf、Q分別為流體等效質(zhì)量、流固耦合、流體等效剛度以及非定常氣動力系數(shù)矩陣；F 結(jié)構(gòu)外載荷向量分別為節(jié)點(diǎn)位移、節(jié)點(diǎn)速度、節(jié)點(diǎn)加速度矢量；ρf為流體密度；p流體壓力矢量。

3 機(jī)翼顫振測試仿真研究方法

在上述理論分析的基礎(chǔ)上，對B737-800機(jī)翼進(jìn)行顫振特性分析。首先簡化機(jī)翼模型如下：機(jī)翼材料采用LY12合金，泊松比為0.33，彈性模量為7.1×1010N/m3，密度為 2 780 kg/m3。劃分網(wǎng)格，如圖 1所示。然后分別通過Block-Lanezos算法[15]和弱耦合算法分析靜態(tài)和耦合流場條件下的機(jī)翼模態(tài)，可得到對應(yīng)的耦合條件下機(jī)翼的振動響應(yīng)，并且通過對比可得耦合流場對機(jī)翼振動模態(tài)的影響；將耦合模態(tài)映射到其他動力特性分析中，可得到耦合流場條件下的振動特性規(guī)律，分析研究后，找到機(jī)翼顫振失穩(wěn)點(diǎn)。

圖1 B737-800有限元模型Fig.1 B737-800 finite element model

3.1 機(jī)翼耦合流場計算

借助Fluent流體分析軟件，采用流固耦合的方法計算機(jī)翼流場壓力。首先利用SolidWorks軟件，以構(gòu)建的B737-800有限元模型為基礎(chǔ)，生成三維模擬機(jī)翼流場；采用Gambit對生成的機(jī)翼周圍的流場進(jìn)行體網(wǎng)格劃分，再導(dǎo)入Fluent軟件中。在Fluent軟件中，選擇基于壓力的穩(wěn)態(tài)求解器，SA湍流模型，選擇理想氣體Ideal-Gas（滿足氣體狀態(tài)方程），并選擇Sutherland定律計算其粘性；設(shè)定速度進(jìn)口，自由出口為初始邊界條件，選擇耦合求解方法，指定對應(yīng)的控制方程，獲得的機(jī)翼表面的壓力分布與流場計算結(jié)果如圖2和圖3所示。

圖2 機(jī)翼壓力分布圖Fig.2 Wing pressure distribution

圖3 流場速度分布圖Fig.3 Flow field velocity distribution

3.2 機(jī)翼模態(tài)仿真研究

分別分析靜態(tài)條件和耦合條件下的模態(tài)響應(yīng)，觀察振型變化，可得到耦合流場對模態(tài)的影響。靜態(tài)條件下只考慮機(jī)翼自身的重力，采用Block-Lanezos算法分析該條件下機(jī)翼的固有振型和頻率。耦合條件下，修改Fluent中獲得的機(jī)翼表面壓力場的數(shù)據(jù)格式，將其通過Fluent的輸出接口加載到Ansys中，在Ansys中重新生成對應(yīng)的機(jī)翼表面壓力數(shù)據(jù)組，讀取該數(shù)據(jù)組，并將其加載到相應(yīng)的有限元模型表面，從而完成耦合力場的導(dǎo)入。首先分析加載耦合力場后的機(jī)翼靜態(tài)結(jié)構(gòu)力學(xué)，獲得相應(yīng)的預(yù)應(yīng)力效應(yīng)矩陣；然后進(jìn)入模態(tài)求解器，采用弱耦合算法，擴(kuò)展模態(tài)，分析機(jī)翼在結(jié)構(gòu)剛度矩陣改變后的模態(tài)振型。圖4為機(jī)翼在靜態(tài)條件和耦合流場條件下的模態(tài)對比（左側(cè)為靜態(tài)條件下，右側(cè)為耦合條件下）。

圖4 靜態(tài)條件和耦合條件下的前4階機(jī)翼模態(tài)對比Fig.4 Comparison of first 4 wing modals under static condition and coupled condition

觀察圖4機(jī)翼在不同條件下的模態(tài)，考慮耦合力場后，機(jī)翼振型明顯發(fā)生了變化：1階模態(tài)變化微弱，但機(jī)翼的相對位移值減少了；2階振型發(fā)生了明顯的變化，在耦合流場條件下，機(jī)翼整體彎扭模態(tài)被嚴(yán)重削弱了；3階模態(tài)靠近翼尖處的彎扭模態(tài)被嚴(yán)重削弱了；4階振型翼尖處的彎曲模態(tài)被嚴(yán)重削弱了。通過對比分析得知，耦合條件下機(jī)翼的模態(tài)振幅均有所下降，但考慮到機(jī)翼的模態(tài)頻率變化微弱，可以認(rèn)為這種變化是由機(jī)翼在空氣流場條件下的耦合作用引起的，如表1所示。

表1 兩種條件下模態(tài)頻率對照Tab.1 Modal frequency comparison under two conditions

本文只對這兩種條件下機(jī)翼的前4階模態(tài)進(jìn)行分析研究，可發(fā)現(xiàn)考慮耦合力場后，機(jī)翼的模態(tài)頻率略有提高，這說明在空氣流場作用下，機(jī)翼產(chǎn)生了剛度硬化現(xiàn)象。兩種條件下，機(jī)翼的變形位移都沿半翼展方向逐漸增大，在70%翼展位置急劇增大，且最大位移均集中在翼尖位置；彎扭振型隨頻率的增加變得復(fù)雜。

圖5為機(jī)翼耦合流場下的各階振型對應(yīng)的應(yīng)力、應(yīng)變響應(yīng)?？捎^測到，1階振動對應(yīng)的最大應(yīng)力、應(yīng)變集中在20%～70%翼展區(qū)，2階振動對應(yīng)的最大應(yīng)力、應(yīng)變集中在85%翼展區(qū)至翼尖位置；其余2階振動對應(yīng)對應(yīng)的最大應(yīng)力、應(yīng)變集中翼尖位置。

圖5 耦合條件下的應(yīng)力應(yīng)變云圖Fig.5 Stress&strain nephogram of modal analysis

3.3 機(jī)翼諧響應(yīng)仿真研究

以機(jī)翼耦合模態(tài)為基礎(chǔ)，選取頻率范圍0～100 Hz，將Fluent中求解得到的耦合力轉(zhuǎn)換格式，加載到機(jī)翼對應(yīng)表面上，采用模態(tài)疊加法求解諧波響應(yīng)特性曲線。分析頻譜特性曲線，尋找振動幅度最大時的頻率，作為全機(jī)翼顫振特性的重點(diǎn)分析頻率，由此可提高顫振測試的針對性。圖6為機(jī)翼X、Y、Z 3個方向的頻譜圖。觀察其響應(yīng)結(jié)果可知X、Y、Z 3個方向均是在頻率為3 Hz時的機(jī)翼振動幅度達(dá)到最大，因此重點(diǎn)觀察3 Hz時機(jī)翼的振型變化，如圖7所示，分析可知在該頻率下，機(jī)翼發(fā)生了嚴(yán)重的彎曲變形，并且最大變形位移位于翼尖處；最大應(yīng)力集中于23%～70%翼展區(qū)。

圖6 X、Y、Z方向頻譜圖Fig.6 X，Y，Z direction spectrum

圖7 后處理響應(yīng)云圖Fig.7 Postprocessing nephogram

3.4 瞬態(tài)動力學(xué)仿真研究

圖8 機(jī)翼表面壓力隨時間變化歷程Fig.8 Change of surface pressure with time

為了獲得機(jī)翼在耦合流場中的瞬態(tài)振動特性曲線以及機(jī)翼的受力分布情況，為顫振測試點(diǎn)的選取提供借鑒作用，通過Fluent模擬求解得到機(jī)翼表面壓力載荷隨時間變化的曲線，如圖8所示。將該時間歷程處理后導(dǎo)入到Ansys中，加載到相應(yīng)的機(jī)翼表面，進(jìn)行機(jī)翼振動瞬態(tài)動力學(xué)分析[16]，得到機(jī)翼表面壓力響應(yīng)，如圖9所示。分析可知，機(jī)翼在0.37 s時響應(yīng)最大，在該時刻對機(jī)翼模態(tài)進(jìn)行擴(kuò)展處理，觀察最大響應(yīng)時刻機(jī)翼的振動響應(yīng)情況。圖10為機(jī)翼在0.37s時刻的振型?？梢钥吹?，機(jī)翼在響應(yīng)最大時刻發(fā)生了嚴(yán)重的彎扭變形，并且最大位移變形和最大應(yīng)力均位于翼尖處。

圖9 等效壓力響應(yīng)Fig.9 Equivalent pressure response

圖10 0.37 s時刻機(jī)翼振型圖Fig.10 Airfoil model at 0.37 s

4 實(shí)例驗(yàn)證

4.1 靜氣動彈性試驗(yàn)

為了預(yù)測飛機(jī)的靜氣動特性，保證飛行安全，中國空氣動力研究與發(fā)展中心[17]對靜氣動彈性機(jī)翼模型進(jìn)行了高速風(fēng)洞試驗(yàn)研究。機(jī)翼模型與真實(shí)機(jī)翼幾何外形、剛度均相似，采用外雙梁內(nèi)三梁結(jié)構(gòu)，梁架與蒙皮由復(fù)合材料碳纖維加工而成。翼展長1.383 m，翼根弦長0.579 m。利用模型變形視頻測量（VMD）技術(shù)對其進(jìn)行變形測量，試驗(yàn)馬赫數(shù)為0.3～1.2。機(jī)翼模型的測量點(diǎn)分布在20%、30%、40%、50%、60%和70%半翼展位置，如圖11所示。可看出，本文的機(jī)翼顫振測試仿真結(jié)果和風(fēng)洞試驗(yàn)的測試點(diǎn)分布情況基本一致。

圖11 靜彈性測試機(jī)翼標(biāo)定點(diǎn)分布Fig.11 Target distribution of static aeroelastic testing model

4.2 動氣動彈性試驗(yàn)

為了測得大展弦比機(jī)翼的動氣動彈性變形和壓力場變化情況，西密歇根大學(xué)的Liu等[18]提出了集成VMD系統(tǒng)和快速響應(yīng)壓力敏感涂料（PSP）以同步評估模擬機(jī)翼模型變形的方法。測試系統(tǒng)及機(jī)翼測試點(diǎn)分布示意圖，如圖12所示。

圖12 風(fēng)洞試驗(yàn)的視頻系統(tǒng)及測試點(diǎn)分布示意圖Fig.12 Video system and test point distribution in wind tunnel test

沿機(jī)翼半翼展方向分布測試點(diǎn)，并測量各個位置的扭轉(zhuǎn)角，如圖13（a）所示，曲線表明扭轉(zhuǎn)角沿機(jī)翼半翼展方向逐漸變大，并且可以看出，在70%翼展處，扭轉(zhuǎn)角會急劇增大。在B737-800機(jī)翼翼尖位置施加不同的扭矩，模擬仿真得到機(jī)翼扭轉(zhuǎn)角沿半翼展方向的變化趨勢如圖13（b）所示。對比發(fā)現(xiàn)，仿真結(jié)果與試驗(yàn)結(jié)果的扭轉(zhuǎn)角趨向一致。

圖13 沿翼展方向的機(jī)翼扭轉(zhuǎn)角分布的測量Fig.13 Spanwise wing twist distribution measurements

由機(jī)翼的靜氣動彈性風(fēng)洞試驗(yàn)和動氣動彈性試驗(yàn)可知，本文的顫振特性分析結(jié)果與機(jī)翼的測試點(diǎn)的分布基本一致；并且仿真的扭轉(zhuǎn)角分布趨勢與動氣動彈性試驗(yàn)的扭轉(zhuǎn)角分布趨勢相符。

5 結(jié)語

本文針對機(jī)翼顫振適航測試的激振點(diǎn)區(qū)域選取的問題，以B737-800機(jī)翼有限元模型為研究對象，采用計算立體流場的方法獲得機(jī)翼表面壓力，將Ansys與Fluent聯(lián)合起來求解機(jī)翼流場耦合-預(yù)應(yīng)力耦合模態(tài)，分析機(jī)翼固有振型在耦合條件下發(fā)生的變化；再應(yīng)用模態(tài)法對機(jī)翼進(jìn)行諧波響應(yīng)分析和瞬態(tài)動力學(xué)分析，獲得機(jī)翼的頻率響應(yīng)和位移響應(yīng)特性，進(jìn)而得到機(jī)翼的振動規(guī)律，找出了機(jī)翼振動危險點(diǎn)，為機(jī)翼顫振適航測試提供了借鑒依據(jù)，提高了測試分析的針對性和可靠性，減少了顫振測試工作量，提高了工作效率。本文的顫振測試結(jié)果與NASA蘭利研究中心進(jìn)行的顫振變形測試分析中選用的激振點(diǎn)的分布區(qū)域一致，并且扭轉(zhuǎn)角沿半翼展方向的趨向也相符，驗(yàn)證了機(jī)翼顫振測試仿真結(jié)果的正確性。

[1]中國民用航空局.CCAR-25-R4，運(yùn)輸類飛機(jī)適航標(biāo)準(zhǔn)[S].中國民用航空局,2011.

[2]趙忠良,吳軍強(qiáng),李 浩,等.2.4 m跨聲速風(fēng)洞虛擬飛行試驗(yàn)技術(shù)研究[J].航空學(xué)報,2016(2):504-512.

[3]孫亞軍,梁 技,楊 飛,等.超臨界機(jī)翼跨音速顫振風(fēng)洞試驗(yàn)研究[J].振動與沖擊,2014(4):190-194.

[4]LI X L,XU L J,TAN C,et al.Real-Time Measurement of Aerodynamic Deformation of Wing by Laser Range Finder[C]//Instrumentation and Measurement Technology Conference(I2MTC),2010:1581-158.

[5]BARROWS, D A. Videogrammetric Model Deformation Measurement Technique for Wind Tunnel Applications(Invited)[C]//AIAA-2007-1163,45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada,Jan 8-11,2007.

[6]BA KUNOWICZ J,MEYER R.In-Flight Wing Deformation Measurements on a Glider[C]//RaeS Research Aircraft Opearations Conference.2014:1-12.

[7]NV NGUYEN,M TYAN,JW LEE,et al.Investigations on stability and control characteristics of a CS-VLA certified aircraft using wind tunnel test data[J].Proceeding of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering,2016,230(14):2728-2743.

[8]LIU J M,JIANG X H,QIN Y L,et al.Numerical simulation of bird impact on solid-element hollow blades by using fluid-solid coupling method[J].Journal of Aerospace Power,2010,25(10):2211-2216.

[9]宗 寧,楊廣珺,呂勝利.大展弦比機(jī)翼流固耦合迭代求解模擬.航空計算技術(shù)[J].2013,43(3):13-16.

[10]朱世權(quán),李海元,陳志華,等.彈性機(jī)翼靜氣動彈性數(shù)值研究[J].工程力學(xué),2017,34(S):326-332.

[11]汪亞真,劉 明,王成軍,等.可變結(jié)構(gòu)小型多旋翼無人飛行器的設(shè)計[J].機(jī)械設(shè)計與制造,2016(4):186-189.

[12]仲繼澤,徐自力.基于動網(wǎng)格降階算法的機(jī)翼顫振邊界預(yù)測[J].振動與沖擊,2017(4):185-191.

[13]聶雪媛,黃程德,楊國偉.基于CFD/CSD耦合的結(jié)構(gòu)幾何非線性靜氣動彈性數(shù)值方法研究[J].振動與沖擊,2016,35(8):48-53.

[14]楊秋明,朱永峰,劉 清.基于流-固耦合傳熱的熱氣防冰系統(tǒng)干空氣飛行蒙皮溫度場計算研究[J].空氣動力學(xué)報,2016,34(6):721-724.

[15]DING J M,WANG Z G,ZHONG C B.Model and harmonic response analysis of composite resin concrete[J].Advanced Materials Research,2011,391-392:349-353.

[16]秦可偉,馬貴春,席 園,等.基于ANSYS的機(jī)翼動力學(xué)分析[J].航空計算技術(shù),2014(2):106-109.

[17]楊賢文,余 立,呂彬彬,等.靜氣動彈性模型高速風(fēng)洞試驗(yàn)研究[J].空氣動力學(xué)學(xué)報,2015,33(5):667-672.

[18]LIU T,MONTEFORT J,GREGORY J,et al.Wing Deformation Measurements from Pressure Sensitive Paint Images Using Videogrammetry[C]//AIAA Fluid Dynamics Conference and Exhibit,2011:343-349.