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

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

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

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

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

1 計算方法

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

1.1 渦格模型

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

圖1 渦格模型示意圖

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

(1)

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

對于定常問題:

vn=-V∞·n

(2)

對于非定常問題:

vn=(Vs-V∞)·n

(3)

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

1.2 載荷計算

渦格上作用的空氣動力與渦強度有關,它由兩部分組成,一部分是與渦強度時間導數(shù)相關的Ft,另一部分是與渦強度空間導數(shù)相關的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 增量法求結構靜變形

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

KTΔq=ΔP-ΔR

(7)

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

KT=K0+Kσ

(8)

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

1.4 模態(tài)分析

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

(9)

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

(10)

1.5 分析流程圖

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

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

2 非定常氣動力驗證

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

?=?msin(kt)

(11)

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

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

3 算例及結果分析

3.1 計算模型

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

圖4 結構有限元模型

表1 主要振型及頻率

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

3.2 線性顫振分析

來流密度ρ∞=2 kg/m3,通過調整來流速度,判斷結構響應是否等幅震蕩來判斷是否發(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°),說明結構非線性將給顫振帶來不利影響。圖8給出了機翼前緣線隨攻角的變化,由于結構非線性的影響,隨著攻角的增大,每1°攻角產生的位移增量逐漸減小。表3給出了振動模態(tài)的頻率隨攻角的變化,隨著攻角的增大,機翼彎曲頻率變化幅度不大,相比之下扭轉頻率減小的幅度很大,10°攻角時機翼一扭頻率只有0°攻角時機翼一扭頻率的61%,這也正是顫振速度隨攻角減小的主要原因。

3.4 結果比較

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

a) 結構非線性+氣動力非線性;

b) 結構非線性+氣動力線性;

c) 結構線性+氣動力線性

結果如圖9所示。可見,當只考慮結構非線性變形,氣動網格不隨平衡位置變化時,計算出來的顫振速度變化趨勢與考慮氣動網格變化的結果趨勢完全相反,原因解釋如下:

圖9 考慮非線性對計算結果的影響

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

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

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

4 結 論

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

1) 隨著攻角的增大，由于結構非線性的影響，每1°攻角產生的位移增量逐漸減小。

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

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

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

從研究結果來看, 傳統(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]∥第十一屆全國空氣彈性學術交流會，昆明：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]. 中國科學：技術科學, 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