基于非線性氣動力的失速顫振計算與試驗研究

戴玉婷，嚴 慧，王林鵬

(北京航空航天大學航空科學與工程學院，北京 100191)

自1980 年以后，科學界對動態(tài)失速和失速顫振的仿真計算有很大興趣，各種計算模型不斷發(fā)展。在過去幾十年中，許多研究人員通過實驗，數(shù)值模擬等方式[1 ? 4]開展了研究。目前人們普遍認為沒有通用的動態(tài)失速模型，即使是先進的CFD 計算[5 ? 9]也不能提供一致的結果。Dimitriadis和Li[10]對NACA0012 翼型的俯仰和沉浮自由度做了一些試驗研究。證明沉浮運動對失速顫振的發(fā)生沒有影響，并得出失速顫振呈現(xiàn)出不同的分支的結論。Yabili 等[11]等利用自己制定的求解器對二元翼段低雷諾數(shù)下的動態(tài)失速和失速顫振進行了CFD 計算，并與Li 的試驗結果做了對比驗證。

對于工程應用而言，半經(jīng)驗動態(tài)失速模型在直升機和風力渦輪機空氣動力學中仍然普遍存在。在這些半經(jīng)驗動態(tài)失速模型中，Leishman-Beddoes 模型是最受歡迎的模型之一，主要適用于0.3 Ma～0.8 Ma。Sheng 等[12]提出了一種改進的動態(tài)失速模型，該模型適用于0.3 Ma 以下氣動力計算。由Galbraith 等[13]開發(fā)的一種新的失速發(fā)生標準已被用來取代基于Evans-Mort 相關性的原始標準。任勇生和劉廷瑞[14]研究了復材薄壁梁非線性失速顫振特性，得出結構阻尼在抑制顫振、增強結構穩(wěn)定性方面效果明顯。

失速顫振試驗相關研究也從動態(tài)失速[15]、試驗規(guī)律、影響因素[16]等方面有了新的進展。隨著計算機，激光和相關技術的發(fā)展，二維和三維粒子圖像測速(PIV)技術成功地用于測量翼型周圍的流場速度[17 ? 18]，分析非定常氣動力產(chǎn)生機理。Bhat 和Govardhan[19]使用PIV 和測壓元件測量了NACA 0012 俯仰振蕩翼型周圍流場和非定常氣動力。Uruba[20]通過振蕩模式分解(OPD)方法分析了俯仰振蕩的NACA 0012 翼型周圍的PIV速度場。?idlof 等[21]展示了一種自振蕩俯仰的翼型設計方法，該模型適用于中等雷諾數(shù)下流動引起的振動。

在Sheng 等[12]的修改模型基礎上，本文對標準Leishman-Beddoes 非線性非定常氣動力模型參數(shù)進行了修正，建立二維翼段氣動彈性模型進行失速顫振計算，并根據(jù)計算結果設計并完成了失速顫振風洞試驗，進一步驗證了修正的L-B 模型可以用來進行低速翼型的失速顫振工程分析與極限環(huán)振蕩評估。

1 失速顫振計算方法

1.1 非線性氣動力計算方法

原始L-B 模型適用的速度范圍在0.3 Ma～0.8 Ma之間，為了使該非定常氣動力模型適用于低速情況，本文在Beddoes 和Sheng[3, 11]研究基礎上，充分考慮翼型表面渦對氣動載荷的影響，在低馬赫數(shù)條件下對L-B 模型進行修正，使其適用于低速條件的非定常氣動力計算。

在上升沿失速過程中，渦的形成和分離靠近前緣位置，在低馬赫數(shù)時附面層在分離的時候仍有部分附著在翼型的上表面。因此，其附著效應將為上表面帶來一個附加升力，稱之為“超調”。超調升力是低馬赫數(shù)下特有的量，該部分的升力系數(shù)可由下式計算[22]：

式中：B1=1.0為與翼型有關的參量； f為弦向分離點的位置； f′′為初始模型中對于分離點 f的一階修正；Vx為由渦產(chǎn)生的法向力系數(shù)的形函數(shù)。

式中： τ為無量綱時間；Tv=7.2為沿弦向傳播的時間常數(shù)；Tvl=4.5為渦傳播距離的時間常數(shù)。

式中：B2=0.32為翼型相關參量； τv為渦傳遞的無量綱時間參數(shù)。

在下降途中，氣流重新附著到翼型表面時，會受到渦的影響使得附著延遲，進而實際氣動力要小于未做延遲處理的模型所得值。該差值為下降沿的超調值。αmin0=16.57°表示下降沿氣流再次附著階段的起始角度，Tr=7.02為氣流再次附著階段渦運動時間常數(shù)，當α<αmin0時氣流再次進入附著階段，翼型上表面渦運動停止，氣流由分離進入完全附著階段。在下降沿中將超調過程開始時對應的攻角設為αmin0，則超調過程與超調開始對應的攻角之間的關系見下式：

式中： Tr為重新附著過程中的時間延遲常數(shù)； q 為減縮頻率。

則回程過程中，在下降沿中法向力系數(shù)的超調值可寫成：

式中， Vxr為 回程中基于翼型的法向力系數(shù)的形函數(shù)，其計算式為：

式中， τr為回程的無量綱時間。

對應的力矩系數(shù)：

低速修正模型最終的法向力系數(shù)和力矩系數(shù)：

1.2 失速顫振計算方法

圖1 為二維翼段氣動彈性系統(tǒng)，該系統(tǒng)具有俯仰和沉浮自由度。綜合1.1 節(jié)非定常氣動力計算方法，則該系統(tǒng)的氣動彈性方程可以寫為：

式中： θ為翼型的俯仰角； h為翼型的垂向位移；m 為翼型單位長度質量； Kh和 Kθ分別為沉浮和俯仰彈簧剛度； Ch和 Cθ分別為沉浮和俯仰阻尼系數(shù)；Sθ為質量靜距； Iθ為翼型相對1/4 弦線處的極慣性矩； Qh和 Qθ分別為外部氣動力和外部氣動力矩；其表達式在1.1 節(jié)中已求得，具體如下：

圖1 二元翼段氣動彈性系統(tǒng)Fig.1 Two-dimensional airfoil aeroelastic system

式中：CL和CM分別為升力系數(shù)和氣動力矩系數(shù)；c 為弦長；l 為展長；xh為彈性軸距1/4 弦線的距離；ρ 為空氣密度；V 為來流速度。

L-B 模型中的攻角α 和俯仰減縮角速率q 可由下式進行計算：

將式(3)改寫為時域狀態(tài)空間方程形式，可得到：

式中：

當給定初始條件x ，利用顯式或者隱式積分可由求得其時域響應：在本文中，使用四階定步長龍格庫塔法求解該狀態(tài)空間方程。

1.3 數(shù)值計算及結果分析

為驗證低馬赫數(shù)下的修正模型，仍然選用NACA0012 翼型。并將計算結果與原始模型和試驗以及Sheng 的模型的結果[3]做了對比。驗證中來流速度為0.1 Ma。對比圖中給出了深度失速條件(攻角α=15?+10?sinωt)下修正L-B 模型的計算結果與試驗對比，對比結果如圖2 所示。

圖2 攻角為α=15?+10?sinωt 時翼型的法向力與力矩系數(shù)Fig.2 Normal force and moment coefficient of airfoil when α=15?+10?sinωt

從圖2 可以看出，在低馬赫數(shù)下，原始L-B模型計算數(shù)據(jù)較試驗數(shù)據(jù)有較大差異，主要體現(xiàn)在超調位置，原始模型和試驗數(shù)據(jù)重合性較差。而修正之后的L-B 模型在超調位置上與試驗數(shù)據(jù)吻合的較好。說明修正之后的L-B 模型有較好的精度，能夠用來計算翼段在低速情況下的氣動力特性。模型中俯仰力矩系數(shù)還受 B2、 Tvl等因素影響，所以力矩對的不是特別好。

1.4 失速顫振結果分析

本文選用參考文獻[10]中的試驗條件和試驗結果作為數(shù)值算例的參數(shù)，來驗證失速顫振計算結果的正確性。文獻展示了在低馬赫數(shù)下基于NACA0012 翼型的失速振蕩試驗。該試驗的重要參數(shù)見下表1。

表1 試驗重要參數(shù)Table1 Important parameters of the test

本文計算了二維翼段失速顫振的俯仰分岔圖和極限環(huán)振蕩曲線，并與試驗[10]的計算結果進行了對比，對比如圖3 所示。

圖3 氣動彈性系統(tǒng)的俯仰分岔圖和極限環(huán)振蕩曲線圖Fig.3 Pitch bifurcation diagram and limit cycle oscillation graph of aeroelastic system

從分岔圖中可以看出，當風速為12.8 m/s，翼段開始周期性俯仰振蕩，即開始分岔。當風速較小時，尤其當風速V 低于15 m/s 時，計算數(shù)據(jù)與試驗數(shù)據(jù)吻合較好，但是在風速V 大于15 m/s 時，俯仰角振蕩幅值明顯高于試驗值。且從試驗數(shù)據(jù)來看，隨著風速的增加俯仰振蕩幅值未有明顯變化。這與結構中沉浮剛度和俯仰剛度可能存在非線性或阻尼非線性等因素有關。

2 失速顫振風洞實驗

2.1 試驗設計

為了探究失速顫振現(xiàn)象的形成機理，根據(jù)風洞的參數(shù)，設計并制作了二維翼段失速顫振的風洞試驗模型，進行了多種狀態(tài)下的風洞試驗。失速顫振試驗在北京航空航天大學D4 開口風洞中完成，開口段尺寸為1.5 m。開口風速范圍小于70 m/s。

翼段翼型采用與仿真條件相同的NACA0012翼型，弦長為0.36 m，翼展為1 m。翼段翼肋之間相距0.125 m。另外，在翼弦向20%和70%處分別布置前梁和后梁。厚度和展長能保證最大瞬時風洞堵塞因子小于5%，滿足試驗要求。在前梁兩端處安裝加速度傳感器，轉軸處安裝角位移計。

翼段質量為2.363 kg，轉軸位于弦長的37.4%位置處，重心處于弦長的39.04%位置，繞轉軸的轉動慣量 Iθ為0.135 kg·m2。支撐系統(tǒng)為轉軸滑軌系統(tǒng)。圖4 為試驗裝置模型，其中，圖4(a)為試驗系統(tǒng)的原理圖，圖4(b)為試驗系統(tǒng)的實物圖。圖4(a)中：1 為八個線性彈簧；2 為四只初始攻角調節(jié)器；3 為兩組滑塊系統(tǒng)；4 為沉浮剛度調節(jié)的線性彈簧；5 為支撐系統(tǒng)；6 為兩組轉軸搖桿機構；7 為角位移傳感器；8 為翼段轉軸；9 為加速度傳感器；10 為翼段。最終翼段沉浮剛度 kh為4896.595 kg/m，俯仰剛度 kθ為135.28 N·m/rad。角位移傳感器的精度為0.1 度。

試驗目的是通過給定初始攻角，測量翼段發(fā)生顫振的速度，以及發(fā)生顫振后翼段的俯仰振幅。在試驗過程中，通過調節(jié)拉桿兩端連接器(圖4(a)中2)的長度，由角度傳感器測量結果，將初始攻角調節(jié)到要求的位置。然后給定一系列穩(wěn)定風速，測量該初始攻角和給定速度下的角加速度和角位移，判斷失速顫振特性。

2.2 試驗結果與仿真對比分析

試驗就四個初始攻角進行了測試，依次為7.5 度、10 度、13 度、16 度。試驗風速從6 m/s～20 m/s。圖5 為16 度初始攻角時，在翼段轉軸處的角度傳感器在不同風速下測得角位移動態(tài)曲線圖。從下列圖中可以看出在風速為6 m/s 時，曲線幾乎平直，幅值近乎為0。風速增加時，角度幅值信號隨風速的增加而增加，且呈等幅震蕩形式。將不同風速下的俯仰角振蕩曲線幅值提取出來，則可生成攻角為16 度時翼段俯仰角度幅值隨風速變化的分岔圖，分岔圖見圖6。圖6 為不同初始攻角下翼段俯仰角隨風速變化的分岔曲線圖，從圖中可以看出，在初始攻角為16 度時，翼段發(fā)生俯仰自由度等幅振蕩的風速為6.5 m/s，即分岔風速為6.5 m/s。在初始攻角一定的情況下，俯仰角的幅值隨風速的增加而增加。另外，當初始攻角增大時，俯仰角的分岔速度明顯降低。

圖4 翼段氣彈系統(tǒng)風洞試驗模型Fig.4 Wind tunnel test model for wing airframe system

根據(jù)第2 節(jié)二維氣動彈性計算模型與試驗模型參數(shù)，對試驗模型失速顫振特性進行了數(shù)值模擬。模型的阻尼大小按照經(jīng)驗取0.1 倍的剛度。圖7 為初始攻角為16 度時風洞試驗與計算模型在不同風速下的對比圖。其中：圖7(a)為風速8 m/s時俯仰角響應；圖7(b)為風速10 m/s 時的俯仰角響應。從圖中可以看出，仿真計算振幅和試風洞試驗振幅吻合度較好。在風速較小時，受風速調節(jié)控制系統(tǒng)、開口外界環(huán)境影響和支撐支座的影響，氣動彈性時域響應有一定的波動。圖8 為初始攻角16 度、風速10 m/s 試驗和仿真的振蕩頻率曲線對比，從圖中可以看出，振蕩頻率吻合較好。圖9 為不同初始攻角下，試驗測得的翼段俯仰角等幅振蕩的幅值隨風速變化的分岔曲線與數(shù)值計算的分岔曲線的對比圖。從圖中可以看出，除了在初始攻角為7.5 度時稍有偏差，數(shù)值計算的結果和試驗測得結果吻合較好。尤其是初始攻角較大時，吻合度較高。說明修正后的L-B 模型能夠較好應用于大攻角失速顫振分析。

圖5 初始攻角16 度不同風速下翼段失速顫振試驗曲線Fig.5 Wing stall flutter test curve of different wind speed at initial angles of 16 degrees

圖6 不同初始攻角下翼段俯仰角隨風速變化的分岔曲線Fig.6 Bifurcation curve of wing pitch angle with wind speed under different initial angles of attack

圖7 初始攻角16 度試驗和數(shù)值模擬計算曲線對比Fig.7 Test and numerical simulation calculation comparison curve at initial attack angle of 16 degrees

圖8 初始攻角16 度風速10 m/s 試驗和數(shù)值振蕩頻率曲線對比Fig.8 Comparison curve of test and numerical oscillation frequency with initial attack angle of 16 degrees and wind speed of 10 m/s

圖9 不同初始攻角下翼段俯仰角隨風速變化的試驗和計算分岔曲線對比Fig.9 Comparison of test and calculated bifurcation curves of wing pitch angle with wind speed under different initial angles of attack

3 結論

本文著眼于非定常氣動力的不確定性對氣動彈性帶來的影響，對標準的氣動力動態(tài)失速模型進行了低速修正，建立了用于考察失速顫振現(xiàn)象的二維氣動彈性模型，根據(jù)風洞試驗的結果，可得出以下結論:

(1) 修正的L-B 模型可以用于二元翼段或低速大展弦比飛機大迎角失速時的非定常氣動力分析。修正模型在低速狀態(tài)下和試驗數(shù)據(jù)吻合較好，尤其是在去程和回程的超調部位，明顯優(yōu)于原始模型。

(2) 設計了二元翼段大攻角失速顫振試驗，試驗表明初始攻角越大，失速顫振分岔速度越低，極限環(huán)振蕩幅值越大。風洞試驗結果與非線性失速顫振理論計算結果吻合較好。說明基于修正L-B失速氣動力模型具有工程上可以接受的精度。