考慮初始迎角影響的二維翼型跨聲速顫振邊界預(yù)測

高國柱, 白俊強, 李國俊, 劉南, 李宇飛

(1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2.中航工業(yè)空氣動力研究院, 遼寧 沈陽 110034)

在現(xiàn)代飛行器設(shè)計中，隨著巡航速度的逐步提高及結(jié)構(gòu)重量系數(shù)的降低，機翼顫振問題日益突出，極易造成結(jié)構(gòu)破壞而引發(fā)災(zāi)難性的后果，嚴(yán)重影響了飛行安全。所以，在飛行器設(shè)計過程中，嚴(yán)格要求在飛行包線內(nèi)不能出現(xiàn)顫振現(xiàn)象，且必須留有一定的顫振裕度，顫振邊界成為約束飛行包線的一個重要因素。然而，目前大量的民用和軍用運輸機都在跨聲速區(qū)域巡航飛行，跨聲速飛行時的非定常激波運動及附面層分離會導(dǎo)致顫振裕度遠遠低于其他速度區(qū)域，且經(jīng)常引發(fā)各種復(fù)雜的非線性顫振現(xiàn)象[1-3]，顯著地影響飛行安全和任務(wù)性能，對跨聲速顫振邊界的預(yù)測以及顫振特性分析成為近年的研究熱點。

傳統(tǒng)的飛行器的顫振分析一般在零度迎角下計算和試驗[4]，忽略了來流非零迎角對顫振邊界的影響。然而飛行器翼型和風(fēng)力機葉片經(jīng)常在非零迎角狀態(tài)下工作，平衡位置與零迎角不同，顫振邊界也會隨之改變。如果在飛行器的顫振設(shè)計校核階段不對迎角影響加以考慮，可能導(dǎo)致設(shè)計出來的飛行器存在顫振安全隱患，導(dǎo)致結(jié)構(gòu)破壞，嚴(yán)重時會造成機毀人亡。

20世紀(jì)90年代，基于Euler/RANS方程的顫振數(shù)值模擬方法逐漸成熟，研究人員針對不同的構(gòu)型進行了大量的數(shù)值模擬，但是針對流動過程中自由來流非零迎角影響的顫振問題則鮮有深入研究。首先有針對性、目的性和全面性地對此類問題展開研究的學(xué)者是NASA的Langley研究中心的Edwards等人[5]。早在1983年，Edwards基于跨聲速小擾動方程，采用時域方法研究了馬赫數(shù)和迎角對顫振邊界的影響。Edwards等人的計算結(jié)果表明在馬赫數(shù)0.775到0.800之間，增加1.5°迎角會導(dǎo)致顫振臨界速度降低50%。1999年，Dowell等人[6]研究了間隙非線性引起的極限環(huán)振蕩時，迎角逐漸增加會引起顫振臨界速度的改變以及由此引起的一系列非線性行為。2002年，Dowell等人[7]進一步研究了低速時迎角改變對顫振邊界的影響。2006年，Park等人[8]研究有間隙的三維全動控制舵面在亞聲速和跨聲速區(qū)內(nèi)，初始迎角變化對氣動彈性響應(yīng)特性的影響。2013年，Bichiou等人[9]采用非定常渦格法，研究了迎角對顫振臨界速度和極限環(huán)幅值的影響，結(jié)果表明，迎角增加會使得顫振臨界速度提高、極限環(huán)幅值減小。

國內(nèi)從20世紀(jì)90年代開始研究迎角對顫振臨界速度的影響。在1996年，楊永年采用非線性渦格法探究了大迎角顫振分析方法。在2004年，葉正寅[10]通過風(fēng)洞試驗和數(shù)值模擬研究了在低速情況下，隨著迎角的增加，顫振臨界速度先略有上升，在靜態(tài)失速點附近，顫振臨界速度大幅度下降，迎角進一步增加后，顫振臨界速度又上升。在2007年，葉正寅和張陳安的研究結(jié)果表明在較高雷諾數(shù)的全湍流條件下，NACA0012翼型在氣動力的線性段，顫振臨界速度隨迎角增加略有上升，但在接近失速迎角附近時，顫振臨界速度則隨迎角增加快速下降。在2008年，張偉偉等人[11]采用系統(tǒng)辨識方法計算削尖三角翼的顫振臨界速度，研究表明迎角增加使得顫振臨界速度減小。在2011年，張偉偉等人[12]基于當(dāng)?shù)鼗钊碚撗芯苛说湫统曀俣婷娴念澱裉匦裕Y(jié)果表明在舵面無間隙時，平均迎角為20°時的顫振速度比5°的顫振速度低8.3%；在舵面有間隙時，迎角越大，激發(fā)極限環(huán)運動所需的初始擾動越大。2012年，劉暢暢等人[13]通過低速風(fēng)洞試驗研究了迎角改變對普通翼型NACA0012和風(fēng)力機翼型NREL S809的顫振特性的影響。

雖然針對考慮迎角影響的顫振研究已經(jīng)取得了一定的成果，但是大部分研究工作仍然停留在考慮迎角影響的低速和超聲速顫振邊界預(yù)測，對于在跨聲速區(qū)域內(nèi)顫振邊界預(yù)測和影響機制缺乏深入的認(rèn)識。本文基于非定常RANS方程和結(jié)構(gòu)運動方程，建立時域氣動彈性分析方法。使用該方法對Isogai案例A模型[14-15]的零迎角和考慮非零迎角影響的跨聲速顫振邊界進行研究，并從激波和分離流動的角度分析迎角對顫振特性的影響機制。

1 流場求解

1.1 非定常氣動力求解

本文采用課題組自研的CFD代碼-TeAM求解非定常氣動力，其控制方程是三維可壓縮非定常積分形式的N-S方程：

?Ω?Q?t+?(F-Fv)?x+?(G-Gv)?y+?(H-Hv)?zdΩ=0

(1)

式中,Q=ρρuρvρweT為守恒向量,ρ,(u,v,w)和e分別為密度、直角坐標(biāo)系下的速度分量和單位質(zhì)量氣體的總能量。F,G,H是3個方向的無黏矢通量,Fv,Gv,Hv是3個方向的黏性矢通量。其中無黏項采用Roe格式進行離散,黏性項采用二階中心差分進行離散,時間推進采用雙時間步法進行迭代求解,為封閉方程引入k-ωSST湍流模型。為了進一步提高計算效率,采用多重網(wǎng)格加速收斂技術(shù)、并行計算技術(shù)。

1.2 跨聲速非定常氣動力的精度驗證

準(zhǔn)確預(yù)測非定常氣動力是氣動彈性分析的基礎(chǔ),首先以NACA64A010強迫俯仰運動[16]為例,研究網(wǎng)格對氣動力的影響。參照阻力預(yù)測會議發(fā)布的網(wǎng)格生成標(biāo)準(zhǔn)得到的3套網(wǎng)格:粗網(wǎng)格、中等網(wǎng)格和密網(wǎng)格,其中中等網(wǎng)格如圖1所示。算例以Davis發(fā)布的中等非線性算例(存在激波)的CT6為準(zhǔn)。翼型強迫振動的迎角隨時間的變化方程為α(t)=α0+αmsin(ωt),k=ωb/V∞,其中α0為翼型初始迎角,αm為翼型強迫振動振幅,ω為翼型強迫振動角頻率,參考長度b為翼型的半弦長,k為減縮頻率,xa/c為轉(zhuǎn)軸位置。計算狀態(tài)為:Ma∞=0.796,α0=-0.21°,αm=1.01°,k=0.202,xa/c=0.248。

參考文獻僅提供了CT6的非定常氣動力隨俯仰角的試驗結(jié)果。計算結(jié)果與試驗結(jié)果對比如圖2和圖3所示,升力系數(shù)隨迎角變化曲線的計算結(jié)果與試驗吻合較好。但是力矩曲線有一定的差別,從和試驗結(jié)果的對比可見,本文的計算結(jié)果與Bendiksen[1]吻合較好,關(guān)于和試驗結(jié)果的差別問題,Bendiksen[1]認(rèn)為通過目前的湍流模型和薄層N-S方程無法預(yù)測圖中非定常力矩的“膨脹部分”,當(dāng)然還有可能是風(fēng)洞試驗的不足。由CT6結(jié)果可見,中等網(wǎng)格對激波運動的模擬精度與密網(wǎng)格基本一致,而粗網(wǎng)格則有所欠缺。同時,該算例的計算結(jié)果驗證了本文發(fā)展的求解方法對跨聲速非定常氣動力的計算具有較高的精度。

圖1 NACA64A010翼型中等網(wǎng)格示意圖 圖2 NACA64A010翼型升力系數(shù)隨迎角變化曲線 圖3 NACA64A010翼型俯仰力矩系數(shù)隨迎角變化曲線

2 結(jié)構(gòu)運動方程

以具有浮沉h和俯仰α2個自由度的典型的二元翼段來作為顫振的研究對象,該二元機翼顫振系統(tǒng)如圖4所示。具有浮沉和俯仰兩自由度的二維翼型結(jié)構(gòu)運動方程為:

m+Sα+Dh+Khh=-L

(2)

Sα+Iα+Dα+Kαα=My

(3)

式中，m為機翼質(zhì)量,Sα為二維翼型對剛心的質(zhì)量靜矩,Dh為浮沉阻尼,Kh為翼型關(guān)于剛心的沉浮剛度,Iα為二維翼型對剛心的質(zhì)量慣性矩,Dα俯仰阻尼,Kα為翼型關(guān)于剛心的俯仰剛度,L為升力,My為俯仰力矩。

圖4 二維翼型沉浮/俯仰兩自由度顫振問題幾何和結(jié)構(gòu)參數(shù)定義[17]

針對上述二維兩自由度翼型結(jié)構(gòu)運動方程,進行無量綱化可得:

M+D+Kq=Q

(4)

式中

q=hb

αM=1xα

D=2ζhωh0

Q=-Lmb

Mymb2

無量綱顫振速度定義為:

Vf=U∞μωαb

(5)

式中，U∞為自由來流速度,μ=m/πρ∞b2為質(zhì)量比,ρ∞為自由來流密度。

為了便于時域求解,引入狀態(tài)變量x=q

=f(x,t)=A·x+B·F

(6)

式中，A=0I

-M-1K-M-1D,B=0

M-1,F=0

Q。

本文采用基于預(yù)估-校正技術(shù)的四階隱式Adams線性多步法[18]對方程(6)進行時域推進求解:

n+1=xn+Δt24(55fn-59fn-1+37fn-2-9fn-3)=xn+Δt24(55A·xn-59A·xn-1+

37A·xn-2-9A·xn-3)+Δt24(55B·Fn-59B·Fn-1+37B·Fn-2-9B·Fn-3)

xn+1=xn+Δt24[9f(n+1,t+Δt)+19fn-5fn-1+fn-2]=xn+Δt24(9A·n+1+19A·xn-

5A·xn-1+A·xn-2)+Δt24(55B·Fn-59B·Fn-1+37B·Fn-2-9B·Fn-3)

(7)

該方法既保證了方程的求解效率,又具有較好的魯棒性。

在有迎角的情況下,翼型除動態(tài)位移外,定常氣動力會產(chǎn)生靜變形位移。在計算過程中,必須考慮靜變形帶來的影響。計算靜變形時,只需將方程(4)的動態(tài)項去掉,即可得到靜變形位移量:

Kqs=Q

(8)

對于不同迎角的氣動彈性數(shù)值模擬,先采定常流場程序,計算定常狀態(tài)下的流場,等待流場基本穩(wěn)定后,啟動計算翼型靜變形程序。當(dāng)翼型靜變形計算收斂后,給翼型一個擾動,翼型就會在氣動力、慣性力和彈性力三者的共同作用下進入動態(tài)響應(yīng)過程。根據(jù)響應(yīng)的特征,判斷是否發(fā)生顫振,從而得到顫振臨界速度。

3 跨聲速顫振邊界驗證

跨聲速顫振凹坑[19-20]與空氣的可壓縮性和激波相位滯后效應(yīng)密切相關(guān),因此跨聲速顫振邊界預(yù)測的關(guān)鍵在于激波捕捉的準(zhǔn)確性。本文首先對Isogai案例A模型的跨聲速顫振邊界進行預(yù)測,以驗證本文顫振計算方法的可靠性。其中計算采用k-ωSST湍流模型,結(jié)構(gòu)參數(shù)采用文獻[21]中的參數(shù),雷諾數(shù)按照文獻[22]中對所有馬赫數(shù)給定為6×106。計算網(wǎng)格為上述非定常氣動力驗證時所用的中等網(wǎng)格,如圖1所示。

圖5和圖6分別為在零迎角條件下計算得到的顫振速度邊界和顫振頻率比邊界,與參考文獻[22-23]的結(jié)果吻合較好。顫振邊界表現(xiàn)為“雙臺階”形,在Ma=0.84和0.91處顫振速度出現(xiàn)2次躍升,與Timme等的計算結(jié)果吻合較好。可以把顫振邊界大致分為4個區(qū)域[24]:①亞聲速區(qū)域Ma<0.80;②跨聲速凹坑Ma=0.80～0.83;③第一次顫振速度躍升區(qū)域Ma=0.84～0.91;④鎖定區(qū)域Ma>0.91。

第一次顫振速度躍升(Ma>0.83)與翼型表面分離流動的出現(xiàn)和擴大密切相關(guān)。圖7展示了5個不同馬赫數(shù)下的表面摩阻分布對比,圖8展示了顫振邊界和準(zhǔn)定常升力線斜率隨馬赫數(shù)變化曲線。從圖中結(jié)果可以看出,當(dāng)0.83≤Ma≤0.87時,隨著馬赫數(shù)增大,翼型表面的分離區(qū)擴大,使得升力線斜率急劇減小,導(dǎo)致顫振邊界急劇增大。

圖5 顫振速度邊界(α0=0°) 圖6 顫振頻率比邊界(α0=0°)圖7 不同馬赫數(shù)下的表面摩阻分布對比(α0=0°)

能量概念在理解和解釋顫振問題中具有重要地位,Bendiksen[25]針對跨聲速顫振提出了“跨聲速穩(wěn)定性法則”(transonic stabilization laws),從氣動力做功的角度分析振動機理,研究了氣動力的幅值及相位和顫振穩(wěn)定性之間的關(guān)系,此處的相位指的是機翼俯仰力矩相對于俯仰位移的滯后相位角。本文采用類似的分析方法對Isogai案例A模型在零迎角和非零迎角條件下的跨聲速穩(wěn)定性進行分析,其中k為減縮頻率。

為了分析顫振速度隨迎角變化的原因,以Isogai翼型的簡諧強迫俯仰運動為研究對象,分析不同馬赫數(shù)下氣動力對翼型的做功情況。通過力矩系數(shù)的幅值和相位差能反映俯仰力矩對翼型的做功情況。由于俯仰軸位于翼型前緣前,必然在氣動中心之前,也就是說翼型抬頭產(chǎn)生的是低頭力矩,反之亦然,所以Isogai翼型的俯仰運動和縱向俯仰力矩在定常情況k=0時就具有180°的相位差。圖9、圖10和圖11分別為俯仰運動中力矩系數(shù)的幅值、相位以及做功情況隨馬赫數(shù)變化曲線。

圖8 顫振邊界和準(zhǔn)定常升力線斜率隨馬赫數(shù)變化曲線 圖9 俯仰運動中力矩系數(shù)的幅值隨馬赫數(shù)變化曲線(α0=0°) 圖10 俯仰運動中力矩系數(shù)的相位隨馬赫數(shù)變化曲線(α0=0°)

圖11 俯仰運動中氣動力做功隨馬赫數(shù)變化曲線(α0=0°)

當(dāng)來流馬赫數(shù)Ma<0.80時,流場無激波,但隨馬赫數(shù)的增加,流動壓縮性的影響使力矩幅值緩慢增大,相位差的變化不明顯(10°以內(nèi)),氣動力對翼型做功也逐漸增加,因此顫振速度隨馬赫數(shù)的變化趨勢由壓縮性主導(dǎo)。

當(dāng)Ma=0.80～0.84時,激波的往復(fù)運動導(dǎo)致力矩幅值明顯增大,且出現(xiàn)了顯著的相位滯后,氣動力對結(jié)構(gòu)輸送的能量顯著增加,顫振邊界出現(xiàn)“跨聲速凹坑”。

當(dāng)Ma=0.85～0.91時,激波誘導(dǎo)的附面層分離占據(jù)主導(dǎo)地位,激波后方分離區(qū)域不僅大幅降低了力矩幅值,還帶來了極大的相位滯后,導(dǎo)致氣動力對翼型的做功降低,使顫振速度大幅提高。

當(dāng)Ma≥0.91時,空間馬赫數(shù)云圖如圖12所示。

圖12 迎角為0°時馬赫數(shù)云圖

由圖可以看出,此時激波達到翼型后緣,流場進入“凍結(jié)區(qū)域”(freeze region),以致翼型運動對激波和壓力分布的影響很小,翼型上下表面非定常壓力分布的幅值和相位差都維持在較小的量級,力矩系數(shù)的幅值、相位差以及做功也非常小。從圖10中可以得知,在“凍結(jié)區(qū)域”的相位角均接近180°,這意味著此時氣動力基本不做功,系統(tǒng)處于臨界穩(wěn)定狀態(tài),翼型的顫振速度隨馬赫數(shù)變化很小。

4 考慮迎角影響的跨聲速顫振邊界預(yù)測

初始迎角對跨聲速激波位置、強度以及翼型表面分離流動有很大的影響,進而會影響跨聲速顫振邊界。本文對Isogai案例A模型進行考慮迎角影響的跨聲速顫振邊界的數(shù)值模擬。為了消除網(wǎng)格帶來的影響,采用與零迎角計算所用相同的氣動網(wǎng)格進行數(shù)值模擬。結(jié)構(gòu)參數(shù)以及氣動參數(shù)和零迎角計算保持一致。

圖13和圖14分別為零迎角和非零迎角條件下計算得到的顫振速度邊界對比圖和頻率比對比圖。數(shù)值模擬結(jié)果表明,考慮迎角影響計算得到的顫振邊界較零迎角計算得到的顫振邊界在Ma≤0.73時基本保持不變;當(dāng)0.73≤Ma≤0.76時,隨著初始迎角增加,顫振速度減小,最大可減小12.5%;當(dāng)0.77≤Ma≤0.86時,隨著初始迎角增加,顫振速度較零迎角時逐漸增加,在零迎角時跨聲速凹坑最低點的位置(Ma=0.83)顫振速度增幅最大,增加了124%,即初始迎角影響使得跨聲速凹坑程度削弱,且跨聲速凹坑最低點前移;當(dāng)0.87≤Ma≤0.90時,考慮初始迎角影響的顫振邊界較零迎角顫振邊界上移,第2個凹坑程度削弱;當(dāng)Ma≥0.91時,考慮初始迎角影響的顫振邊界較零迎角的偏大。

圖13 顫振速度邊界

圖14 顫振頻率比邊界

為了探究初始迎角影響引起跨聲速顫振邊界變化的原因,選取4個典型的馬赫數(shù)進行不同迎角發(fā)生顫振時的靜變形計算結(jié)果對比,分別為Ma=0.70,0.78,0.83,0.94。

當(dāng)Ma=0.70時,如圖15a)所示,顫振臨界速度隨初始迎角增加時先基本保持不變,當(dāng)初始迎角增加到4°時,顫振臨界速度開始減小。從圖16a)和圖17a)可以看出,當(dāng)初始迎角小于4°時,靜變形后激波強度變化不明顯,而且翼型表面無分離流動,此時由于激波和分離的共同作用,使得顫振臨界速度變化不明顯。當(dāng)初始迎角大于4°時,靜變形后激波強度迅速增加,翼型上表面出現(xiàn)小范圍的分離流動,此時激波強度變化占據(jù)主導(dǎo)地位,隨著迎角增加,顫振臨界速度逐漸減小。

當(dāng)Ma=0.78時,如圖15b)所示,顫振臨界速度隨初始迎角增加而先減小后增加。從圖16b)和圖17b)可以看出,當(dāng)初始迎角小于3°時,靜變形后翼型表面無分離,但是激波強度隨著迎角的變化而迅速增強,導(dǎo)致顫振臨界速度隨著迎角增加而逐漸減小。當(dāng)初始迎角大于4°時,靜變形后激波強度隨著迎角變化不明顯,但是分離區(qū)迅速擴大。此時分離流動的影響占據(jù)主導(dǎo)地位,使得顫振臨界速度增加。

當(dāng)Ma=0.83時,如圖15c)所示,處于零迎角時跨聲速凹坑內(nèi),顫振臨界速度隨著初始迎角增加而劇烈增加。從圖16c)和圖17c)可以看出,隨著初始迎角增加,靜變形后激波強度變化不大,但在翼型上表面出現(xiàn)分離流動,并且分離區(qū)迅速擴大。此時分離流動的影響占據(jù)主導(dǎo)地位,顫振臨界速度則隨著初始迎角逐漸增加。

當(dāng)Ma=0.94時,如圖15d)所示,顫振臨界速度隨著初始迎角增加而小幅增加。從圖16d)和圖17d)可以出,此時靜變形后激波強度變化不大,但在翼型上表面的后緣處出現(xiàn)分離流動。此時由于分離區(qū)的影響,顫振臨界速度則隨著迎角增加逐漸略有上升。

圖15 顫振速度邊界和顫振頻率比邊界隨迎角變化曲線

圖16 靜變形后表面摩阻分布隨初始迎角變化曲線

圖17 靜變形后壓力系數(shù)隨初始迎角變化曲線

5 結(jié) 論

本文基于RANS方程和結(jié)構(gòu)運動方程,建立時域氣動彈性分析方法,分別對零迎角和考慮非零迎角影響的Isogai案例A模型的跨聲速顫振邊界進行預(yù)測,數(shù)值模擬結(jié)果表明:

1)本文采用零迎角計算得到的顫振邊界和文獻的結(jié)果吻合很好,驗證了本文建立方法的可靠性。

2)在Ma≤0.73時迎角變化對顫振臨界速度基本沒有影響;當(dāng)0.73≤Ma≤0.76時,隨著初始迎角增加,顫振臨界速度減小;當(dāng)0.80≤Ma≤0.86時,考慮初始迎角影響,顫振速度較零迎角時增加,非零迎角影響使得跨聲速凹坑程度削弱,且跨聲速凹坑最低點前移;當(dāng)0.87≤Ma≤0.90時,考慮初始迎角影響的顫振邊界較零迎角顫振邊界上移,第2個凹坑程度削弱;當(dāng)Ma≥0.91時,非零迎角影響的顫振速度較零迎角的增加。

3)無論迎角是否為零,當(dāng)激波到達翼型后緣,進入凍結(jié)區(qū)域,氣動力對翼型幾乎不做功,系統(tǒng)處于臨界穩(wěn)定狀態(tài),相同迎角下的顫振速度隨馬赫數(shù)的增大而變化很小。

4)結(jié)果表明,在更高的馬赫數(shù)下,增加初始迎角,會增加系統(tǒng)的穩(wěn)定性。這也表明增加初始迎角可以提高顫振臨界速度,從而作為一種抑制顫振的手段。

綜上所述,激波和翼型表面的分離現(xiàn)象共同作用,影響了翼型的顫振特性,其中激波增強會降低系統(tǒng)的穩(wěn)定性,分離區(qū)域的擴大會增加系統(tǒng)的穩(wěn)定性。來流迎角對激波強度和翼型表面分離區(qū)域大小均有影響,使得翼型的氣動力非線性增強,顫振特性與零迎角條件下有明顯區(qū)別。

因此針對機翼進行顫振設(shè)計時,需要對來流迎角的影響加以考慮,以便得到精確的顫振邊界,確保飛機的飛行安全。

參考文獻：

[1] Bendiksen O O. Review of Unsteady Transonic Aerodynamics: Theory and Applications[J]. Progress in Aerospace Sciences, 2011, 47: 135-167

[2] Bendiksen O O. Transonic Flutter Characteristics of Advanced Fighter Wings[C]∥56th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Florida, 2015

[3] Cunningham A M. The Role of Non-Linear Aerodynamics in Fluid-Structure Interaction[C]∥29th AIAA Fluid Dynamics Conference, 1998

[4] 葉正寅, 趙令誠. 亞音速、時間域內(nèi)機翼的氣動彈性分析方法[J]. 西北工業(yè)大學(xué)學(xué)報, 1989, 7(4):456-463

Ye Zhengyin, Zhao Lingcheng. Aeroelastic Characteristics of Wings in Subsonic Flow[J]. Journal of Northwestern Polytechnical University, 1989, 7(4): 456-463 ( in Chinese)

[5] Edwards J W, Bennet R M, Whitlow W, et al. Time-Marching Transonic Flutter Solutions Including Angle-of-Attack Effects[J]. Journal of Aircraft, 1983,20(11): 899-906

[6] Kholodar D B, Dowell E H. Behavior of Airfoil with Control Surface Freeplay for Nonzero Angles of Attack[J]. AIAA Journal, 1999, 37(5): 651-653

[7] Dowell E, Tang D. Nonlinear Aeroelasticity and Unsteady Aerodynamics[C]∥40th AIAA Aerospace Sciences Meeting & Exhibit, Reno, NV, 2002

[8] Park Y, Yoo J, Lee I. Effects of Angle-of-Attack on the Aeroelastic Characteristics of a Wing with Freeplay[J]. Journal of Spacecraft and Rockets, 2006, 43(6): 1419-1422

[9] Bichiou Y, Nuhait A, Abdelkefi A, et al. Unsteady Aeroelastic Behaviors of Rigid Airfoils with Preset Angles of Attack[J]. Journal of Vibration and Control, 2016, 22(4): 1010-1022

[10] 葉正寅. 不同迎角下翼型的氣動彈性性質(zhì)研究[C]∥首屆全國航空航天領(lǐng)域中的力學(xué)問題學(xué)術(shù)研討會, 成都, 2004: 240-242

Ye Zhengyin. The Study of Aeroelastic Property of Airfoil at the Different Angles of Attack[C]∥The 1st Proseminar about the Mechanical Problems in the Fields of Aeronautics and Aerospace, Chengdu, 2004: 240-242 (in Chinese)

[11] Zhang W, Ye Z. Effects of Angles of Attack on Flutter Characteristics of a Cropped Delta Wing[C]∥49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Schaumburg, IL, 2008

[12] 張偉偉, 王忠波, 葉正寅, 等. 迎角對間隙舵面的非線性顫振特性影響研究[J]. 機械強度, 2011, 33(2):296-301

Zhang Weiwei, Wang Zhongbo, Ye Zhengyin, et al. Effect of Angle of Attack on Flutter Characteristics of a Control Surface with Free-Play Nonlinearity[J]. Journal of Mechanical Strength, 2011, 33(2):296-301 (in Chinese)

[13] 劉暢暢, 劉子強, 季辰. 不同迎角的翼型氣彈特性風(fēng)洞試驗研究[J]. 空氣動力學(xué)學(xué)報, 2012,30(2):271-276

Liu Changchang, Liu Ziqiang, Ji Chen. Aeroelastic Characteristics of Airfoils at Different Angles of Attack in Wind Tunnel Testing[J]. Acta Aerodynamica Sinica, 2012,30(2): 271-276 (in Chinese)

[14] Alonso J, Jameson A. Fully-Implicit Time-Marching Aeroelastic Solutions[R]. AIAA-1994-0056

[15] Zhang Zhichao, Yang Shuichi, Liu Feng. Prediction of Flutter anf LCO by an Euler Method on Non-Moving Cartesian Grids with Boundary-Layer Corrections[R]. AIAA-2005-0833

[16] Davis S S. NACA 64A010(NASA AMES model) Oscillatory Pitching[R]. AGARD Report No. 702

[17] Thomas J P, Dowell E H, Hall K C. Nonlinear Inviscid Aerodynamic Effects on Transonic Divergence, Flutter, and Limit-Cycle Oscillations[J]. AIAA Journal, 2002, 40(4): 638-646

[18] 葉正寅, 張偉偉, 史愛明,等. 流固耦合力學(xué)基礎(chǔ)及其應(yīng)用[M]. 哈爾濱: 哈爾濱工業(yè)大學(xué)出版社, 2010:171-173

Ye Zhengyin, Zhang Weiwei, Shi Ai′ming, et al. Fundamentals of Fluid-Structure Coupling and Its Application[M]. Harbin, Harbin Institute of Technology Press, 2010:171-173 (in Chinese)

[19] Koji Isogai. On the Transonic-Dip Mechanism of Flutter of a Sweptback Wing[J]. AIAA Journal, 1979, 17(7):793-795

[20] Koji Isogai. Transonic Dip Mechanism of Flutter of a Sweptback Wing: Part Ⅱ[J]. AIAA Journal, 1981, 19(9): 1240-1242

[21] Alonso J, Jameson A. Fully-Implicit Time-Marching Aeroelastic Solutions[R]. AIAA-1994-0056

[22] Zhang Zhichao, Yang Shuichi, Liu Feng. Prediction of Flutter and LCO by an Euler Method on Non-Moving Cartesian Grids with Boundary-Layer Corrections[R]. AIAA-2005-833

[23] Timme S, Badcock K J. Searching for Transonic Aeroelastic Instability Using an Aerodynamic Model Hierarchy[R]. Technical Report Deliverable D2.3, UK: University of Liverpool, 2009

[24] 劉南. 機翼跨聲速非線性顫振及高效分析方法研究[D]. 西安：西北工業(yè)大學(xué)，2016

Liu Nan. Investigation of Transonic Nonlinear Flutter and Efficient Analysis Approach[D]. Xi′an, Northwestern Polytechnical University, 2016 (in Chinese)

[25] Bendiksen O O. Transonic Stabilization Laws for Unsteady Aerodynamics and Flutter[R]. AIAA-2012-1718