戴玉婷，嚴(yán) 慧，王林鵬

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

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

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

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

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

1 失速顫振計(jì)算方法

1.1 非線性氣動(dòng)力計(jì)算方法

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

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

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

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

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

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

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

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

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

式中， τr為回程的無(wú)量綱時(shí)間。

對(duì)應(yīng)的力矩系數(shù)：

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

1.2 失速顫振計(jì)算方法

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

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

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

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

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

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

式中：

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

1.3 數(shù)值計(jì)算及結(jié)果分析

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

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

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

1.4 失速顫振結(jié)果分析

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

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

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

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

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

2 失速顫振風(fēng)洞實(shí)驗(yàn)

2.1 試驗(yàn)設(shè)計(jì)

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

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

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

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

2.2 試驗(yàn)結(jié)果與仿真對(duì)比分析

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

圖4 翼段氣彈系統(tǒng)風(fēng)洞試驗(yàn)?zāi)Ｐ虵ig.4 Wind tunnel test model for wing airframe system

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

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

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

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

圖8 初始攻角16 度風(fēng)速10 m/s 試驗(yàn)和數(shù)值振蕩頻率曲線對(duì)比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 不同初始攻角下翼段俯仰角隨風(fēng)速變化的試驗(yàn)和計(jì)算分岔曲線對(duì)比Fig.9 Comparison of test and calculated bifurcation curves of wing pitch angle with wind speed under different initial angles of attack

3 結(jié)論

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

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

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