李國俊，白俊強，唐長紅，李宇飛，劉 南

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

失速顫振及其伴隨的分離流動等復(fù)雜流動現(xiàn)象對飛行器的飛行性能和飛行安全有一定程度的影響，但是研究人員對失速顫振的發(fā)生機理仍然缺乏系統(tǒng)的認識和深入的研究。通常在失速顫振問題研究中，與失速顫振相關(guān)的動態(tài)失速氣動力現(xiàn)象是風(fēng)洞實驗和數(shù)值模擬的關(guān)注重點。雖然近些年來針對失速顫振本身的研究逐漸增多，但是相比于廣泛存在的各種各樣的失速顫振現(xiàn)象，這些研究還是遠遠不夠的[1]。

早在第一次世界大戰(zhàn)之前，飛機機翼和尾翼的失速顫振現(xiàn)象就已經(jīng)引起研究人員的注意[2]。從20世紀(jì)40年代開始，國內(nèi)外研究工作者開展了一系列關(guān)于失速顫振以及與之密切相關(guān)的動態(tài)失速現(xiàn)象的研究[3-7]，獲得了許多有價值的結(jié)論。通過研究發(fā)現(xiàn)，失速顫振是一種由動態(tài)失速和機翼自身的彈性力、慣性力耦合產(chǎn)生的有限振幅的自激振蕩現(xiàn)象。在機翼振動過程中，伴隨著升力面上的大面積動態(tài)流動分離以及再附現(xiàn)象，此時升力面處于失速攻角附近，氣動力的非線性很強，是一種嚴重的氣動彈性不穩(wěn)定現(xiàn)象。與經(jīng)典顫振相比，能量從氣流中傳輸?shù)秸駝訖C翼上的機理，既不依賴于兩階模態(tài)的彈性耦合或空氣動力學(xué)耦合，也不依賴于位移與其氣動力響應(yīng)之間的相位滯后，而是與動態(tài)失速和分離流動密切相關(guān)。

早期針對失速顫振進行數(shù)值模擬大多偏向于計算精度驗證，隨著研究的深入，研究工作者的研究重點逐漸偏向于失速顫振現(xiàn)象及其誘發(fā)機理的研究。近些年來，Dimitriadis等[8-11]開展了一系列關(guān)于失速顫振的風(fēng)洞實驗，并采用fsiFoam進行了數(shù)值模擬，深入分析了NACA0012翼型的失速顫振現(xiàn)象以及分岔特性。而國內(nèi)早在2002年，金琰等采用非定常RANS方程結(jié)合低雷諾數(shù)兩方程q-ω湍流模型對NACA0015翼型的失速顫振現(xiàn)象進行了研究。目前在失速顫振現(xiàn)象研究方面，大量的研究結(jié)果表明基于非定常RANS方程的氣彈求解器足以捕捉與失速顫振相關(guān)的關(guān)鍵流動特性和物理現(xiàn)象。

在特定的流動情況下，當(dāng)分離漩渦的脫落頻率接近結(jié)構(gòu)自然頻率時，會發(fā)生頻率鎖定(frequency lock-in)現(xiàn)象[12]。在頻率鎖定范圍內(nèi)，結(jié)構(gòu)響應(yīng)的幅值會顯著增加從而引起結(jié)構(gòu)疲勞甚至失效。目前大多數(shù)針對鎖頻現(xiàn)象的研究主要集中在圓柱繞流中，只有少數(shù)研究涉及機翼在渦致振動時出現(xiàn)的鎖頻現(xiàn)象。Tang等[13]開展了一系列NACA0012翼型在大攻角和高雷諾數(shù)下有關(guān)鎖頻現(xiàn)象的風(fēng)洞實驗，研究結(jié)果表明NACA0012翼型在給定的強迫運動幅值和頻率下呈現(xiàn)出“V”形鎖頻區(qū)。Zhu等[14]采用數(shù)值模擬方法研究了一個薄翼型的鎖頻現(xiàn)象。Young等[15-17]針對NACA0012翼型在浮沉運動中發(fā)生的鎖頻現(xiàn)象進行了研究，目的在于實現(xiàn)翼型推進效率的最大化。上述工作大多是針對強迫運動中的鎖頻現(xiàn)象進行研究，鮮有涉及自由振動中的鎖頻現(xiàn)象。

本文基于非定常RANS方程和結(jié)構(gòu)運動方程，建立時域氣動彈性分析方法。首先驗證了本文采用的非定常氣動力求解器模擬動態(tài)失速和預(yù)測鎖頻區(qū)域的可靠性。其次，使用該方法對翼型剖面為NACA23012翼型的剛性矩形機翼的顫振邊界進行驗證。在此基礎(chǔ)上，改變初始攻角和自由來流速度，對該矩形機翼的失速顫振問題進行數(shù)值模擬，分析失速顫振的流動特性和誘發(fā)機制，并對失速顫振中出現(xiàn)的分岔和鎖頻現(xiàn)象進行了研究。

1 流場求解

1.1 非定常氣動力求解

本文采用課題組自研的CFD代碼-TeAM求解非定常氣動力，其控制方程是三維可壓縮非定常積分形式的N-S方程，其直角坐標(biāo)系的守恒形式積分方程為

(1)

1.2 動態(tài)失速氣動力響應(yīng)預(yù)測

準(zhǔn)確預(yù)測動態(tài)失速的氣動力響應(yīng)是研究失速顫振現(xiàn)象的前提。本文擬采用非定常RANS方程對NACA0012翼型的動態(tài)失速氣動力響應(yīng)進行數(shù)值模擬，并與實驗結(jié)果[18]進行對比，以驗證其預(yù)測動態(tài)失速現(xiàn)象的精度。

實驗給定翼型弦長為61 cm，馬赫數(shù)Ma=0.293，雷諾數(shù)Re=3 758 528。翼型強迫振動的攻角隨時間的變化方程為α(t)=α0+αmsin(ωt)，k=ωb/V∞，其中α0為翼型初始攻角，αm為翼型強迫振動振幅，ω為翼型強迫振動角頻率，參考長度b為翼型的半弦長，k為減縮頻率。翼型繞自身的1/4弦線做正弦運動，選定翼型的運動規(guī)律為：α0=15°，αm=5°，k=0.202 4。計算所用網(wǎng)格如圖1所示。

圖1 NACA0012翼型網(wǎng)格示意圖Fig.1 Mesh for NACA0012 airfoil

圖2、圖3分別為NACA0012翼型的非定常升力系數(shù)和力矩系數(shù)隨攻角變化曲線，圖中實線是本文結(jié)果，三角形為風(fēng)洞實驗結(jié)果。從圖中結(jié)果可知，本文求解得到的結(jié)果和風(fēng)洞實驗結(jié)果基本吻合，表明本文采用的非定常氣動力求解器可以捕捉動態(tài)失速的主要物理特征，適用于模擬失速顫振中的非定常氣動力。

圖2 NACA0012翼型升力系數(shù)隨攻角變化曲線Fig.2 Lift Coefficient versus AOA of NACA0012 airfoil

圖3 NACA0012翼型力矩系數(shù)隨攻角變化曲線Fig.3 Moment Coefficient versus AOA of NACA0012 airfoil

1.3 鎖頻區(qū)域預(yù)測

為了確保本文采用的求解器可以應(yīng)用于失速顫振中的鎖頻現(xiàn)象研究，對在不同強迫運動規(guī)律下NACA0012翼型的鎖頻區(qū)域進行預(yù)測，并與文獻結(jié)果進行對比。NACA0012翼型的弦長為0.255 3 m，展長為0.52 m，翼型繞自身的1/4弦線做正弦運動。測定鎖頻區(qū)域的計算狀態(tài)為：雷諾數(shù)Re=200 000，初始攻角α=40°。計算所用網(wǎng)格采用和1.2節(jié)中相同的網(wǎng)格分布，鎖頻區(qū)域?qū)Ρ冉Y(jié)果如圖4所示。

圖4 NACA0012翼型鎖頻區(qū)域Fig.4 The lock-in region of NACA0012 airfoil

圖4中的虛線代表Besem等采用諧波平衡求解器得到的鎖頻區(qū)域，圓點表示本文在給定強迫運動規(guī)律下發(fā)生鎖頻的結(jié)果，三角形點表示不發(fā)生鎖頻的結(jié)果。圖中的橫軸表示無量綱頻率fenforce/fshed，即強迫運動頻率與固定狀態(tài)下NACA0012翼型的漩渦脫落頻率之比，其中fshed=11.617 1 Hz；縱軸表示強迫運動的幅值。從圖中可以看出，本文預(yù)測得到的鎖頻區(qū)域與文獻結(jié)果一致，可以應(yīng)用于失速顫振中的鎖頻現(xiàn)象研究。

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

本文以具有浮沉h和俯仰α兩個自由度的典型的二元翼段作為顫振的研究對象，其無量綱形式的結(jié)構(gòu)運動方程為

(2)

(3)

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

(4)

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

3 算例與分析

3.1 顫振邊界預(yù)測

Bollay等[21]在哈佛大學(xué)的低速風(fēng)洞中對剖面翼型為NACA23012的矩形機翼的失速顫振特性進行了研究。該機翼模型材料為紅木，機翼展長為1 m，弦長為0.171 45 m，實驗所用結(jié)構(gòu)參數(shù)為：xα=0.228，ωh=80.3 rad/s，ωα=87.2 rad/s，μ=161.2，其中rα按照文獻給定為0.5[22]，彈性軸即旋轉(zhuǎn)中心位于弦長的35.5%處。

采用本文建立的時域氣動彈性分析方法對具有浮沉和俯仰兩自由度的NACA23012翼型的顫振邊界進行預(yù)測，以驗證本文建立的顫振數(shù)值模擬方法的可靠性。計算所用網(wǎng)格如圖5所示。

圖6和圖7分別展示了NACA23012翼型在給定計算狀態(tài)下的顫振速度邊界和頻率邊界。其中實線為本文預(yù)測得到的顫振邊界，虛線為實驗結(jié)果。數(shù)值模擬結(jié)果表明，本文預(yù)測得到的顫振速度邊界和實驗結(jié)果吻合較好。從圖7可知，當(dāng)α<15°時，顫振頻率隨著初始攻角的增大而增大；當(dāng)15°<α<17°時，顫振頻率突然下降，其結(jié)果接近翼型的結(jié)構(gòu)固有頻率；當(dāng)α>17°時，顫振頻率又突然增大。

3.2 失速顫振的流動特性和誘發(fā)機制研究

為了探究失速顫振的觸發(fā)機制，選取初始攻角α=15°、自由來流速度v=17 m/s時的響應(yīng)作為分析對象。首先對該狀態(tài)下的結(jié)構(gòu)運動響應(yīng)和流動特性進行分析，確認機翼此時發(fā)生失速顫振，翼型的浮沉運動響應(yīng)和俯仰運動響應(yīng)分別如圖8和圖9所示。從圖中可以看出，浮沉位移很小，振幅不超過10 mm，而俯仰位移較大，振幅在10°左右。結(jié)合圖10中翼型的靜態(tài)升力曲線可以得知，此時翼型的有效攻角在其靜態(tài)失速攻角附近上下變化。

圖5 NACA23012翼型網(wǎng)格示意圖Fig.5 Mesh for NACA23012 airfoil

圖6 NACA23012翼型顫振速度邊界Fig.6 Flutter speed boundary of NACA23012 airfoil

圖7 NACA23012翼型顫振頻率邊界Fig.7 Flutter frequency boundary of NACA23012 airfoil

圖8 初始攻角α=15°、自由來流速度v=17 m/s時浮沉運動響應(yīng)
Fig.8 Plunging response atα=15° andv=17 m/s

圖9 初始攻角α=15°、自由來流速度v=17 m/s時俯仰運動響應(yīng)
Fig.9 Pitching response atα=15° andv=17 m/s

圖10 NACA23012翼型靜態(tài)升力系數(shù)隨攻角變化曲線
Fig.10 Static lift coefficient versus AOA of NACA23012 airfoil

圖11和圖12展示了氣動力系數(shù)隨瞬時攻角變化曲線，結(jié)果表明，升力系數(shù)并未在靜態(tài)失速攻角處出現(xiàn)陡降，而是隨著俯仰位移的增大而增大，直至在正向最大攻角處(26°)出現(xiàn)跳躍，力矩系數(shù)則表現(xiàn)出類似的變化。因為翼型在浮沉自由度上的位移很小，對翼型整體做功貢獻很小，因此忽略浮沉自由度的做功，僅對一個周期內(nèi)的俯仰運動和力矩進行積分得到力矩做功為0.010 58 N·m，此時機翼的非線性氣動力和自身的彈性力、慣性力耦合，從氣流中吸收能量從而產(chǎn)生有限振幅的自激振蕩現(xiàn)象，即失速顫振。

圖11 升力系數(shù)隨瞬時攻角變化曲線Fig.11 Lift coefficient versus instantaneous AOA

圖12 升力系數(shù)隨瞬時攻角變化曲線Fig.12 Moment coefficient versus instantaneous AOA

圖13為一個周期內(nèi)不同時刻的瞬態(tài)流場，其中背景云圖為壓力系數(shù)云圖，t/T代表一個周期內(nèi)的無量綱時刻。圖14為不同時刻的上表面摩阻系數(shù)分布。結(jié)合圖13和圖14對不同時刻的流場和氣動力變化進行分析，可以將一個周期內(nèi)的流場現(xiàn)象分為三個大的階段：

(1)前緣漩渦產(chǎn)生階段(0

當(dāng)翼型正向(上仰)轉(zhuǎn)動時，即圖13(a)～(c)，翼型前緣處開始產(chǎn)生漩渦并逐漸擴大。從圖14(a)中可知，翼型上表面的摩阻系數(shù)分布在t/T=0和0.125時幾乎沒有小于零的區(qū)域，意味著此時翼型上表面幾乎不存在分離流動；當(dāng)t/T=0.25時，翼型前緣處的摩阻系數(shù)分布小于零，即出現(xiàn)分離流動，此時翼型前緣處的負壓增大。

(a)t/T=0

(b)t/T=0.125

(c)t/T=0.25

(d)t/T=0.375

(e)t/T=0.5

(f)t/T=0.625

(g)t/T=0.75

(h)t/T=0.875

(a)前緣漩渦產(chǎn)生階段

(b)漩渦向后移動階段

(c)漩渦逐漸脫落階段

(2)漩渦向后移動階段(0.25

當(dāng)翼型到達正向最大位移處開始負向(下俯)轉(zhuǎn)動時，即圖13(c)和(d)，結(jié)合圖14(b)可知，翼型上表面的摩阻系數(shù)分布在t/T=0.25、0.3、0.35和0.4四個不同時刻均出現(xiàn)小于零的區(qū)域，并且該區(qū)域向后緣移動，這意味著翼型上表面的漩渦開始向后緣移動。

(3)漩渦逐漸脫落階段(0.4

當(dāng)翼型接近初始平衡位置并繼續(xù)轉(zhuǎn)動時，即圖13(e)～(h)，結(jié)合圖14(c)可知，翼型上表面的摩阻系數(shù)分布隨著翼型的不斷運動，其小于零的區(qū)域逐漸減小，這意味著位于翼型上表面的漩渦開始逐漸脫離翼型表面。

圖15展示了一個周期內(nèi)翼型的俯仰力矩系數(shù)和壓力中心位置變化情況，其中Xcp為壓力中心沿翼型弦向的坐標(biāo)，橫向虛線為旋轉(zhuǎn)軸位置。

圖15 一個周期內(nèi)俯仰力矩系數(shù)和壓力中心變化曲線Fig.15 Moment coefficient and pressure center curve in a whole motion period

圖16 一個周期內(nèi)的瞬時能量傳遞系數(shù)Fig.16 The instantaneous energy transfer coefficient in a whole motion period

3.3 失速顫振中的鎖頻及分岔現(xiàn)象研究

失速顫振除了伴隨著分離流動外，由于其流動的特殊性和復(fù)雜性，同時還可能具有其他一些復(fù)雜現(xiàn)象。本小節(jié)將對該翼型在失速顫振中出現(xiàn)的鎖頻及分岔現(xiàn)象進行分析研究。

該翼型的分岔圖如圖17所示，分別展示了初始攻角為0°、5°、15°、16°、19°、20°時在不同自由來流速度下的位移幅值。其中振幅為0的點為收斂狀態(tài)，橫軸為自由來流速度，縱軸為某自由度下的位移幅值。從圖中可以得知，當(dāng)初始攻角為0°、5°時，浮沉和俯仰振幅在來流速度超過顫振臨界速度后急劇增大，出現(xiàn)超臨界Hopf分岔；當(dāng)初始攻角為15°、16°時，浮沉和俯仰振幅在來流速度超過顫振臨界速度后先發(fā)生跳躍，隨后出現(xiàn)小幅下降，接著又發(fā)生跳躍；而當(dāng)初始攻角為19°、20°時，浮沉和俯仰振幅在來流速度超過顫振臨界速度后持續(xù)緩慢增長。

在3.2節(jié)對翼型的顫振邊界進行了驗證和分析，從圖17可以得知，顫振頻率在初始攻角增大至15°時突然減小，隨后在17°時陡增。這種特殊現(xiàn)象與流場和結(jié)構(gòu)的耦合機制有關(guān)。圖18展示了不同初始攻角和自由來流速度下的顫振頻率，在自由來流速度超過顫振臨界速度后，當(dāng)初始攻角為0°和5°時，顫振頻率隨來流速度的增加而增大；當(dāng)初始攻角為15°和16°時，顫振頻率維持在11.3 Hz左右，接近結(jié)構(gòu)的固有頻率，與來流速度大小基本無關(guān)；當(dāng)初始攻角為19°和20°時，顫振頻率隨來流速度的增加而略有減小。為了進一步分析發(fā)生鎖頻時結(jié)構(gòu)振動頻率的變化，對初始攻角為15°時的俯仰力矩系數(shù)進行快速傅立葉變換得到功率譜密度(PSD)示意圖，如圖19所示。圖中結(jié)果表明，結(jié)構(gòu)振動頻率的主峰值不跟隨自由來流速度的變化而變化，而是鎖定在結(jié)構(gòu)固有頻率附近，即發(fā)生鎖頻。結(jié)合以上分析結(jié)果，翼型在初始攻角為15°、16°和17°時發(fā)生鎖頻，導(dǎo)致顫振頻率減小。

綜上所述，失速顫振中出現(xiàn)的分岔及鎖頻現(xiàn)象對其顫振特性有顯著影響，明顯有別于經(jīng)典顫振，因此在失速顫振計算研究中需要特別關(guān)注。

(a)浮沉位移幅值

(b)俯仰位移幅值

圖18 不同初始攻角和自由來流速度下的顫振頻率對比Fig.18 Comparison of flutter frequency at different initial AOAs and freestream velocities

圖19 功率譜密度(PSD)示意圖(α=15°)Fig.19 Power spectrum density(PSD) (α=15°)

4 結(jié) 論

本文基于非定常RANS方程和結(jié)構(gòu)運動方程，建立時域氣動彈性分析方法，對由分離流動誘發(fā)的失速顫振和鎖頻現(xiàn)象進行了研究。結(jié)果表明：

(1)本文采用的非定常RANS求解器可以準(zhǔn)確模擬動態(tài)失速氣動力響應(yīng)和鎖頻區(qū)域。

(2)采用本文建立的時域氣動彈性數(shù)值模擬方法對NACA23012翼型的顫振邊界進行預(yù)測，結(jié)果表明本文預(yù)測得到的顫振速度邊界和實驗結(jié)果吻合較好。

(3)通過對初始攻角α=15°、自由來流速度v=17 m/s時下的結(jié)構(gòu)運動響應(yīng)和流動特性進行分析，發(fā)現(xiàn)在失速顫振中前緣漩渦的產(chǎn)生和尾渦脫落是一種能量轉(zhuǎn)換和注入機制，用以維持翼型的等幅振蕩。

(4)不同的初始攻角和來流速度導(dǎo)致翼型的振動特性發(fā)生分岔，而翼型在初始攻角為15°、16°和17°時發(fā)生鎖頻現(xiàn)象導(dǎo)致其顫振頻率突然減小。