馬艷峰, 賀爾銘, 曾憲昂, 李俊杰, 唐長紅

(1.西北工業(yè)大學 航空學院, 陜西 西安 710072; 2.中航工業(yè)第一飛機設(shè)計研究院, 陜西 西安 710089)

大型運輸機、高空長航時無人機都具有展弦比大、柔性大、變形大等特點，飛機在飛行載荷作用下會產(chǎn)生彎曲和扭轉(zhuǎn)變形，一方面結(jié)構(gòu)平衡狀態(tài)相對于未變形結(jié)構(gòu)有一定差異，結(jié)構(gòu)存在預應力，各振型和頻率也有差別；另一方面，變形前后氣動面的差異也會導致非定常氣動力的不同。這兩方面綜合將使機翼的顫振特性發(fā)生變化[1]。

隨著數(shù)學計算工具的進步和對復雜系統(tǒng)非線性空氣動力學認識的加深，飛行器顫振問題得到了更深入研究。在20世紀90年代早期, Van Schoor等[2]基于完全線性理論對一個大柔性飛機進行了氣動彈性分析, 討論了結(jié)構(gòu)彈性對氣動彈性和飛行動力學穩(wěn)定性的影響。Patil等[3]研究了高空長航時飛行器的氣動穩(wěn)定性和飛行力學特性, 并在考慮幾何非線性的基礎(chǔ)上, 分析了大展弦比機翼的顫振特性。Mayuresh和Dowell[4]結(jié)合理論和實驗方法研究了大展弦比機翼的氣動彈性相應問題, 并由此掀起了大展弦比機翼非線性氣動彈性研究的熱潮[5-6]。國內(nèi)方面，楊智春等[7]研究了非線性大展弦比柔性機冀模型的靜氣動彈性特性, 并通過提取靜平衡位置的剩余剛度, 計算了顫振特性。謝長川、楊超[8]基于動力學小擾動假設(shè)，建立了具有大展弦比機翼柔性飛機的全機幾何非線性氣動彈性穩(wěn)定性分析的線性化方法和工程求解流程。崔鵬、韓景龍[9]考慮了幾何非線性，進行了切尖三角翼的非線性顫振和極限環(huán)振蕩的研究。

本文綜合考慮機翼結(jié)構(gòu)的幾何非線性和氣動力非線性，詳細地分析了大展弦比機翼的非線性顫振特性。以平板機翼為例，計算了機翼在不同攻角下的靜彈性變形、模態(tài)特性以及顫振特性。通過比較顫振計算線性解與非線性解的差別，說明對大展弦比機翼進行顫振分析時，必須同時考慮結(jié)構(gòu)幾何非線性和氣動力非線性。

1 計算方法

常規(guī)顫振計算采用偶極子網(wǎng)格法計算非定常氣動力，當結(jié)構(gòu)發(fā)生變形，偶極子網(wǎng)格無法跟隨結(jié)構(gòu)運動形成曲面，因此它不適用于大變形的非線性顫振分析。針對大變形計算氣動網(wǎng)格要跟隨結(jié)構(gòu)形成曲面的特點[10]，本文選用非定常渦格法來求解非定常氣動力，渦格的形狀較為靈活，沿展向2條端線也無須滿足平行于來流的限制，可形成任意曲面形狀[10-11]，另外一個主要差別是用渦格法計算非定常流動時采用時域法，比偶極子網(wǎng)格多出了從后緣到遠場的尾渦網(wǎng)格。

1.1 渦格模型

渦格模型示意圖如圖1所示[12]。將拓展到遠場的馬蹄渦拆分成若干個獨立渦格，而這些渦格一部分附著于物面上，一部分則形成尾渦，將尾渦網(wǎng)格布置到距翼面后緣30倍弦長為止。

圖1 渦格模型示意圖

設(shè)物面上的渦格共有N=Nc·Ns個(Nc為物面沿弦向的網(wǎng)格數(shù),Ns為沿展向的網(wǎng)格數(shù)),尾渦格共M=Nw·Ns個(Nw為尾渦沿弦向的網(wǎng)格數(shù)),則法向洗流速度和渦強度的關(guān)系可寫成如下的矩陣形式[11]:

(1)

式中：上標‘～’表示尾渦。

對于定常問題:

vn=-V∞·n

(2)

對于非定常問題:

vn=(Vs-V∞)·n

(3)

式中：Vs為物體在該控制點處的運動速度,對于非定常問題,單元法向矢量n隨物體的運動而不斷變化。

1.2 載荷計算

渦格上作用的空氣動力與渦強度有關(guān),它由兩部分組成,一部分是與渦強度時間導數(shù)相關(guān)的Ft,另一部分是與渦強度空間導數(shù)相關(guān)的Fs。

(4)

式中：ρ∞為來流密度,S為渦格面積。Ft作用于渦格的中心,方向沿單元的法向。

當渦格位于機翼前緣時:

i=(m-1)·Nc+1,m=1,2,…,Ns

(5)

當渦格不位于機翼前緣時:

i≠(m-1)·Nc+1,m=1,2,…,Ns

(6)

式中:b為渦格沿展向的跨度,FS作用于渦格前端渦線的中點,方向沿單元的法向,通過杠桿原理將Ft和Fs分配到角點上[11]。

1.3 增量法求結(jié)構(gòu)靜變形

本文只考慮翼面大變形產(chǎn)生的幾何非線性影響,不考慮材料非線性等因素,采用增量法計算,結(jié)構(gòu)的增量平衡方程如下:

KTΔq=ΔP-ΔR

(7)

式中：ΔP為載荷增量,ΔR為內(nèi)力增量,Δq為增量位移,KT為切線剛度矩陣:

KT=K0+Kσ

(8)

式中：K0為線彈性剛度矩陣,Kσ為幾何剛度矩陣(初應力陣)。

1.4 模態(tài)分析

在變形后的平衡位置求結(jié)構(gòu)響應還需求解結(jié)構(gòu)在該位置的動態(tài)特性,在平衡位置的結(jié)構(gòu)動力方程為:

(9)

通過模態(tài)法求解該方程,將(9)式表示為模態(tài)坐標下的減縮形式如下:

(10)

1.5 分析流程圖

本文使用改進的“松耦合”方式來聯(lián)合推進結(jié)構(gòu)和氣動方程。分析時尋求二者之間的匹配關(guān)系,顫振速壓即為產(chǎn)生靜彈變形的速壓,靜彈變形則為顫振分析的平衡狀態(tài)。反復進行靜彈和顫振計算,不斷迭代使二者關(guān)系達成一致,圖2給出了分析流程圖。

圖2 非線性顫振分析流程圖

2 非定常氣動力驗證

建立一個展弦比為24的長直機翼,網(wǎng)格數(shù)為20×40,給它施加一個強迫俯仰運動:

?=?msin(kt)

(11)

式中:?m為俯仰角的振幅,k為減縮頻率,t為無量綱時間,轉(zhuǎn)動的軸線為機翼沿展向的中心線,減縮頻率k=0.2時,升力和俯仰力矩響應歷程與經(jīng)典theodorson理論計算結(jié)果比較如圖3所示。二者的響應曲線吻合好,驗證了渦格法計算的正確性。

圖3 渦格法與theodorson氣動力的比較(k=0.2)

3 算例及結(jié)果分析

3.1 計算模型

計算選用的平板機翼根弦長為1 m,根梢比為1,展長為4 m,展弦比為4,前緣后掠角為8°。結(jié)構(gòu)模型如圖4所示,平板的厚度為0.03 mm,材料彈性模量E=71 GPa,密度為2 700 kg/m3,泊松比為0.33。主要前5階振動頻率見表1。

圖4 結(jié)構(gòu)有限元模型

表1 主要振型及頻率

氣動網(wǎng)格分為物面網(wǎng)格和尾渦網(wǎng)格兩部分,物面網(wǎng)格數(shù)量為20×30,與結(jié)構(gòu)網(wǎng)格數(shù)一致,為了便于數(shù)據(jù)交換,尾渦網(wǎng)格一致延伸到距后緣30倍弦長。

3.2 線性顫振分析

來流密度ρ∞=2 kg/m3,通過調(diào)整來流速度,判斷結(jié)構(gòu)響應是否等幅震蕩來判斷是否發(fā)生臨界顫振,得到Vf=171.47 m/s,ωf=6.115 Hz。在顫振速度下的廣義位移響應如圖5所示。由圖可見,顫振形態(tài)是典型的彎扭耦合型顫振,參與顫振的模態(tài)主要是機翼一彎、二彎和一扭。

圖5 廣義位移響應

3.3 非線性顫振分析

計算攻角取α=0°～10°,攻角增量的間隔取Δα=1°。當α1=0°時,取v1=171.47 m/s。按圖2所示的流程圖進行靜彈和顫振的匹配計算。顫振速度和頻率隨攻角的變化如圖6、圖7和表2所示。

表2 顫振速度和頻率隨攻角的變化

圖6 顫振速度隨攻角的變化

圖7 顫振頻率隨攻角的變化

圖8 不同攻角,顫振速壓下,機翼前緣的靜變形

表3 不同攻角下機翼平衡位置的模態(tài)頻率

如圖6可見,顫振速度呈現(xiàn)出隨攻角先略微增大,后逐漸減小的過程(拐點在α≈2°),說明結(jié)構(gòu)非線性將給顫振帶來不利影響。圖8給出了機翼前緣線隨攻角的變化,由于結(jié)構(gòu)非線性的影響,隨著攻角的增大,每1°攻角產(chǎn)生的位移增量逐漸減小。表3給出了振動模態(tài)的頻率隨攻角的變化,隨著攻角的增大,機翼彎曲頻率變化幅度不大,相比之下扭轉(zhuǎn)頻率減小的幅度很大,10°攻角時機翼一扭頻率只有0°攻角時機翼一扭頻率的61%,這也正是顫振速度隨攻角減小的主要原因。

3.4 結(jié)果比較

以下討論不考慮非線性或是部分考慮非線性對計算結(jié)果的影響。對以下3種情況的計算結(jié)果進行比較:

a) 結(jié)構(gòu)非線性+氣動力非線性;

b) 結(jié)構(gòu)非線性+氣動力線性;

c) 結(jié)構(gòu)線性+氣動力線性

結(jié)果如圖9所示?？梢?當只考慮結(jié)構(gòu)非線性變形,氣動網(wǎng)格不隨平衡位置變化時,計算出來的顫振速度變化趨勢與考慮氣動網(wǎng)格變化的結(jié)果趨勢完全相反,原因解釋如下:

圖9 考慮非線性對計算結(jié)果的影響

如圖10所示,非線性模態(tài)振型基本垂直于變形后的網(wǎng)格,若氣動網(wǎng)格隨結(jié)構(gòu)變化,則(2)式中vn=

圖10 氣動力線性與非線性差異示意圖

|Vs·n|≈|Vs|,而當氣動網(wǎng)格不變時,運動速度Vs和法向量n之間存在夾角,則vn=|Vs·n|<|Vs|,且隨著變形的增大,vn減小的幅度越大。由(1)式可知,vn減小會直接導致渦強度的減小,相應的非定常氣動力減小,從而使計算的顫振速度偏高。所以,在結(jié)構(gòu)非線性求解中考慮氣動面的變形是十分必要的,否則得出的結(jié)果不準確,甚至可能得出相反的結(jié)果。

4 結(jié) 論

本文針對大展弦比機翼進行了靜氣動彈性分析與顫振分析, 并通過考慮結(jié)構(gòu)變形和內(nèi)應力對系統(tǒng)剛度的影響, 形成了機翼幾何非線性氣動彈性分析的實施方法。計算結(jié)果表明：

1) 隨著攻角的增大，由于結(jié)構(gòu)非線性的影響，每1°攻角產(chǎn)生的位移增量逐漸減小。

2) 攻角的增大對機翼彎曲頻率影響較小，相比之下扭轉(zhuǎn)頻率減小的幅度很大。

3) 顫振速度隨攻角的增大呈現(xiàn)出先略微增大而后減小的趨勢，說明結(jié)構(gòu)非線性將給大展弦比機翼顫振帶來不利影響。

4) 氣動網(wǎng)格隨平衡位置變化與保持不變2種情形下計算結(jié)果表明考慮氣動網(wǎng)格變形的必要性。

從研究結(jié)果來看, 傳統(tǒng)的線性分析方法對于大展弦比機翼不適用,因此進一步研究大展弦比機翼的幾何非線性氣動彈性分析方法是很有必要的。

參考文獻：

[1] Dowell E. Some Recent Advances in Nonlinear Aeroelasticity: Fluid-Structure Interaction in the 21st Century[R]. AIAA-2010-3137

[2] Van Schoor M C, Von Flotow A H. Aeroelastic Characteristics of a Highly Flexible Aircraft[J]. Journal of Aircraft, 1990, 27: 901-908

[3] Patil M J. Nonlinear Aeroelastic Analysis, Flight Dynamics, and Control of a Complete Aircraft[D]. Atlanta:Georgia Institute of Technology, 1999

[4] Mayuresh J P, Dewey H H. Nonlinear Aeroelasticity and Flight Dynamics of High-Altitude Long-Endurance Aircraft[J]. Journal of Aircraft, 2001, 38: 88-94

[5] Dowell E, Edwards J, Strganac T. Nonlinear Aeroelasticity[J]. Journal of Aircraft, 2003, 40: 857-874

[6] Weihua Su, Carlos E, Cesnik S. Nonlinear Aeroelasticity of a Very Flexible Blended-Wing-Body Aircraft[J]. Journal of Aircraft,2010, 47: 1539-1553

[7] 楊智春，黨會學，李毅. 大展弦比機翼非線性氣動彈性特性的數(shù)值模擬研究[C]∥第十一屆全國空氣彈性學術(shù)交流會，昆明：2009

Yan Zhichun, Dang Huixue, Li Yi. Research of Nonlinear Aeroelastic Characteris on High-Aspect-Ratio wing[]∥CARS-2009, Kunming, 2009 (in Chinese)

[8] 謝長川，楊超. 大展弦比飛機幾何非線性氣動彈性穩(wěn)定性的線性化方法[J]. 中國科學：技術(shù)科學, 2011, 41(3): 385-393

Xie Changchuan, Yang Chao. Linearization Method of Nonlinear Aeroelastic Stability for Complete Aircraft with High-Aspect-Ratio Wings[J]. Sci China: Tech Sci, 2011, 41(3): 385-393 (in Chinese)

[9] 崔鵬,韓景龍. 基于CFD/CSD的非線性氣動彈性分析方法[J]. 航空學報, 2010, 31(3): 480-486

Cui Peng, Han Jinglong. Investigation of Nonlinear Aeroelastic Analysis Using CFD/CSD[J]. Chinese Journal of Aeronautics, 2010, 31(3): 480-486 (in Chinese)

[10] Ivan W, Chad Gibbsy S, Dowell E. Aeroelastic Analysis of a Folding Wing: Comparison of Simple and Higher Fidelity Models for a Wide Range of Fold Angles[R]. AIAA-2013-1635

[11] Bret K S, Philip S Beran. Analytical Sensitivity Analysis of an Unsteady Vortex Lattice Method for Flapping Wing Optimization[R]. AIAA-2009-2614

[12] Christopher C C, Tim F, Balakumar B. GPGPU Implementation and Benchmarking of the Unsteady Vortex Lattice Method[R]. AIAA-2013-0288, 2013