分離流動(dòng)誘發(fā)的失速顫振和鎖頻現(xiàn)象研究

李國俊，白俊強(qiáng)，唐長紅，李宇飛，劉 南

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

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

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

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

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

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

1 流場求解

1.1 非定常氣動(dòng)力求解

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

(1)

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

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

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

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

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

圖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)象研究，對在不同強(qiáng)迫運(yùn)動(dòng)規(guī)律下NACA0012翼型的鎖頻區(qū)域進(jìn)行預(yù)測，并與文獻(xiàn)結(jié)果進(jìn)行對比。NACA0012翼型的弦長為0.255 3 m，展長為0.52 m，翼型繞自身的1/4弦線做正弦運(yùn)動(dòng)。測定鎖頻區(qū)域的計(jì)算狀態(tài)為：雷諾數(shù)Re=200 000，初始攻角α=40°。計(jì)算所用網(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ū)域，圓點(diǎn)表示本文在給定強(qiáng)迫運(yùn)動(dòng)規(guī)律下發(fā)生鎖頻的結(jié)果，三角形點(diǎn)表示不發(fā)生鎖頻的結(jié)果。圖中的橫軸表示無量綱頻率fenforce/fshed，即強(qiáng)迫運(yùn)動(dòng)頻率與固定狀態(tài)下NACA0012翼型的漩渦脫落頻率之比，其中fshed=11.617 1 Hz；縱軸表示強(qiáng)迫運(yùn)動(dòng)的幅值。從圖中可以看出，本文預(yù)測得到的鎖頻區(qū)域與文獻(xiàn)結(jié)果一致，可以應(yīng)用于失速顫振中的鎖頻現(xiàn)象研究。

2 結(jié)構(gòu)運(yùn)動(dòng)方程

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

(2)

(3)

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

(4)

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

3 算例與分析

3.1 顫振邊界預(yù)測

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

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

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

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

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

圖9 初始攻角α=15°、自由來流速度v=17 m/s時(shí)俯仰運(yùn)動(dòng)響應(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展示了氣動(dòng)力系數(shù)隨瞬時(shí)攻角變化曲線，結(jié)果表明，升力系數(shù)并未在靜態(tài)失速攻角處出現(xiàn)陡降，而是隨著俯仰位移的增大而增大，直至在正向最大攻角處(26°)出現(xiàn)跳躍，力矩系數(shù)則表現(xiàn)出類似的變化。因?yàn)橐硇驮诟〕磷杂啥壬系奈灰坪苄?，對翼型整體做功貢獻(xiàn)很小，因此忽略浮沉自由度的做功，僅對一個(gè)周期內(nèi)的俯仰運(yùn)動(dòng)和力矩進(jìn)行積分得到力矩做功為0.010 58 N·m，此時(shí)機(jī)翼的非線性氣動(dòng)力和自身的彈性力、慣性力耦合，從氣流中吸收能量從而產(chǎn)生有限振幅的自激振蕩現(xiàn)象，即失速顫振。

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

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

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

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

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

(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)漩渦向后移動(dòng)階段

(c)漩渦逐漸脫落階段

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

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

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

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

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

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

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

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

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

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

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

綜上所述，失速顫振中出現(xiàn)的分岔及鎖頻現(xiàn)象對其顫振特性有顯著影響，明顯有別于經(jīng)典顫振，因此在失速顫振計(jì)算研究中需要特別關(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)運(yùn)動(dòng)方程，建立時(shí)域氣動(dòng)彈性分析方法，對由分離流動(dòng)誘發(fā)的失速顫振和鎖頻現(xiàn)象進(jìn)行了研究。結(jié)果表明：

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

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

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

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