李永昌, 戴玉婷, 楊超

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

開裂式阻力方向舵(簡稱阻力舵)是目前飛翼式飛機(jī)的主要航向操縱裝置之一[1-3]。 飛翼式飛機(jī)在飛行時(shí),阻力舵會(huì)對(duì)局部流動(dòng)造成較大的擾動(dòng),與升力、俯仰力矩及滾轉(zhuǎn)力矩耦合較大且變化規(guī)律復(fù)雜。 大開裂角狀態(tài)下的流場(chǎng)類似于大迎角的流場(chǎng)繞流,可能會(huì)引起氣流分離和渦流等復(fù)雜的流場(chǎng)狀態(tài)[4-5],對(duì)飛機(jī)的氣動(dòng)特性造成很大的影響,可能會(huì)導(dǎo)致氣動(dòng)彈性問題。

針對(duì)這種特殊的操縱裝置,國內(nèi)外開展了不同方向的研究。 Mohamad 等[6]討論了阻力舵對(duì)飛翼式飛機(jī)偏航運(yùn)動(dòng)的影響,并給出了低速狀態(tài)下阻力舵的操縱效率隨開裂角的變化規(guī)律。 張子軍等[7]通過風(fēng)洞試驗(yàn)的方法探究了阻力舵對(duì)飛翼式飛機(jī)氣動(dòng)特性的影響,試驗(yàn)表明,在給定的迎角和來流馬赫數(shù)范圍內(nèi),阻力舵的舵效基本不變,并且阻力舵可以在對(duì)升力、側(cè)力、俯仰力矩影響較小的情況下產(chǎn)生較大的偏航力矩。 劉其琛等[8]對(duì)阻力舵動(dòng)態(tài)打開過程的非定常氣動(dòng)特性進(jìn)行了研究,得到了升阻力系數(shù)等數(shù)據(jù)在阻力舵打開過程中的變化趨勢(shì)。 韓鵬等[9]針對(duì)二維阻力方向舵模型的氣動(dòng)彈性問題開展了工作,為這種特殊的氣動(dòng)外形建立了新的顫振計(jì)算方法。

目前的流場(chǎng)分解手段主要有2 種:本征正交分解(proper orthogonal decomposition,POD)方法和動(dòng)力學(xué)模態(tài)分解(dynamic mode decomposition,DMD)方法。 20 世紀(jì)60 年代,Lumley[10]將POD方法應(yīng)用于流體力學(xué)。 POD 方法的核心思想是通過樣本數(shù)據(jù)得到一組恰當(dāng)?shù)恼换瘮?shù),使得樣本數(shù)據(jù)在該正交基上的投影分量依次迅速衰減,從而可以用較少的基展開獲得高維數(shù)據(jù)的近似描述,構(gòu)造低維動(dòng)力學(xué)模型。 DMD 方法是由Schmid 等[11-12]在Koopman 分析的基礎(chǔ)上提出的一種低維系統(tǒng)分解方法,基于流場(chǎng)的動(dòng)力學(xué)特性,對(duì)流場(chǎng)快照進(jìn)行分解,所提取的子結(jié)構(gòu)在時(shí)間、空間演化特性上相互正交,可以得到流動(dòng)在時(shí)間和空間上的主要特征。 DMD 方法經(jīng)過不斷改進(jìn)發(fā)展,近年來被廣泛地應(yīng)用于流動(dòng)現(xiàn)象的機(jī)理分析,如臺(tái)階流動(dòng)[13]、方腔流動(dòng)[14]等。 潘翀等[15]利用DMD 方法對(duì)復(fù)雜流動(dòng)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析,發(fā)現(xiàn)DMD 方法不僅能夠捕捉流場(chǎng)的主要特征,還可以捕捉其高頻諧波。 寇家慶和張偉偉[16]對(duì)現(xiàn)有的DMD 方法與POD 方法進(jìn)行了比較,說明了DMD方法在流場(chǎng)分析中的優(yōu)勢(shì),并討論了DMD 方法的研究現(xiàn)狀及未來發(fā)展方向。 Rowley 等[17]對(duì)射流流場(chǎng)采用DMD 分析,得到了流場(chǎng)中的剪切層渦和壁渦等結(jié)構(gòu)。 韓亞東和譚磊[18]應(yīng)用DMD 方法獲得空化流場(chǎng)的動(dòng)力學(xué)模態(tài),揭示了空化發(fā)展和脫落的結(jié)構(gòu)特征。 Menon 和Mittal[19]對(duì)俯仰運(yùn)動(dòng)的翼型繞流進(jìn)行DMD 分析,獲得了其主要的流場(chǎng)模態(tài),為流固耦合運(yùn)動(dòng)的DMD 分析提供了方法。

本文通過計(jì)算流體力學(xué)(computational fluid dynamics, CFD)方法得到了不同雷諾數(shù)及不同開裂角下的二維阻力舵的流場(chǎng)數(shù)據(jù),并基于DMD方法對(duì)各流場(chǎng)進(jìn)行模態(tài)分解,分析了各個(gè)模態(tài)的流動(dòng)特征及頻率變化。 同時(shí),通過流固耦合方法進(jìn)行了二維阻力舵模型俯仰運(yùn)動(dòng)的機(jī)理分析,得到了不同開裂角下流固耦合系統(tǒng)運(yùn)動(dòng)形式的變化規(guī)律。

1 計(jì)算對(duì)象及方法

本文所用計(jì)算模型如圖1 所示,是一個(gè)帶有阻力舵的NACA0015 翼型,弦長c約為0.3 m,阻力舵安裝在翼型的3/4 弦長位置,模型在1/2 弦長位置連接有扭轉(zhuǎn)彈簧。 流場(chǎng)來流速度為V,迎角為α,規(guī)定阻力舵的開裂角β為單片舵的偏轉(zhuǎn)角度。 流場(chǎng)網(wǎng)格如圖2 所示,其為AMI 網(wǎng)格,分為外部區(qū)域和內(nèi)部區(qū)域,在進(jìn)行流固耦合計(jì)算時(shí),外部網(wǎng)格不動(dòng),內(nèi)部網(wǎng)格沿著邊界進(jìn)行移動(dòng),從而實(shí)現(xiàn)模型的俯仰運(yùn)動(dòng)。 由于主要關(guān)注翼型附近和尾流區(qū)的流動(dòng)特征,主要提取內(nèi)部區(qū)域的計(jì)算數(shù)據(jù)進(jìn)行DMD 分析。 翼型附近靠近壁面的第1 層網(wǎng)格厚度y+約為1。 流場(chǎng)的入口條件為速度入口,出口條件設(shè)置為壓力出口,翼型表面為無滑移壁面,計(jì)算時(shí)間步長為0.000 1 s,計(jì)算時(shí)選擇k-ωSST 湍流模型。

圖1 阻力舵模型Fig.1 Drag rudder model

圖2 流場(chǎng)網(wǎng)格Fig.2 Flow field grid

為驗(yàn)證數(shù)值計(jì)算方法的有效性,圖3 給出了NACA0012 翼型在動(dòng)態(tài)失速狀態(tài)下隨迎角α變化的升力系數(shù)的仿真數(shù)據(jù)與試驗(yàn)數(shù)據(jù)的對(duì)比。 翼型迎角α隨時(shí)間t的變化規(guī)律為:α(t) =10°+15°sin (Ω t),其中,角頻率Ω約為18.67 rad/s。 可以看出,本文選用的轉(zhuǎn)捩SST 模型與文獻(xiàn)中的仿真數(shù)據(jù)[20]及試驗(yàn)數(shù)據(jù)[21]在20°迎角的范圍內(nèi)相差不大。而在20°迎角以上的范圍,2 組數(shù)值計(jì)算的結(jié)果與試驗(yàn)值均有不同程度的誤差。 可以認(rèn)為本文選用的SST 模型在迎角不大時(shí)的計(jì)算數(shù)據(jù)是可靠的。

圖3 不同迎角時(shí)升力系數(shù)的計(jì)算結(jié)果與試驗(yàn)值Fig.3 Numerical and experimental aerodynamic coefficients with different angles of attack

為了確保計(jì)算效率的同時(shí)保證計(jì)算的精度,對(duì)阻力舵的流場(chǎng)網(wǎng)格進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證。 選用了3 套流場(chǎng)網(wǎng)格,其外部網(wǎng)格基本一致,而內(nèi)部網(wǎng)格邊界層向外的網(wǎng)格增長率分別為1.05、1.1、1.15。 各自的網(wǎng)格數(shù)量分別為15 萬(細(xì)網(wǎng)格)、8 萬(中等網(wǎng)格)、3 萬(粗網(wǎng)格)。 計(jì)算選用轉(zhuǎn)捩SST 模型,來流速度為10 m/s,2°迎角,時(shí)間步長為1 ×10-4s。 3 套網(wǎng)格在相同流場(chǎng)條件下計(jì)算得到的翼型表面俯仰力矩隨時(shí)間變化曲線如圖4 所示,力矩中心為翼型0.5 倍弦長位置。 可以看出,中等密度的網(wǎng)格滿足計(jì)算精度。 為了保證計(jì)算精度的同時(shí)兼顧計(jì)算效率,因此,選擇中等密度、網(wǎng)格數(shù)為8 萬的網(wǎng)格進(jìn)行后續(xù)計(jì)算。

圖4 三套網(wǎng)格的俯仰力矩隨時(shí)間變化曲線Fig.4 Pitching moment of three grids versus time

2 DMD 方法

對(duì)于某個(gè)流場(chǎng)來說,一系列時(shí)間間隔為Δt的流場(chǎng)快照矩陣可以用前N+1 個(gè)時(shí)刻的流場(chǎng)信息向量xi表示為

DMD 方法的基本思想是用一個(gè)線性算子來近似復(fù)雜流動(dòng)的演化。 也就是尋找一個(gè)矩陣A,使得=,這樣矩陣A就代表了相鄰2 個(gè)時(shí)刻流場(chǎng)之間的線性關(guān)系。 但是在實(shí)際計(jì)算中,矩陣A是一個(gè)高維矩陣,因此,DMD 方法需要用一個(gè)低維矩陣來近似代替A。

式中:U、V為酉矩陣;Σ為奇異值對(duì)角矩陣;V*為V的共軛轉(zhuǎn)置矩陣。

從而矩陣A的近似矩陣可以表示為

該模態(tài)對(duì)應(yīng)的增長率gi和頻率ωi分別為

式中:Re 為復(fù)數(shù)的實(shí)部;Im 為復(fù)數(shù)的虛部。

3 DMD 結(jié)果分析

對(duì)5°、10°、20°開裂角狀態(tài)下的二維阻力舵模型進(jìn)行非定常氣動(dòng)力計(jì)算,在流場(chǎng)穩(wěn)定后提取網(wǎng)格節(jié)點(diǎn)處的壓強(qiáng)數(shù)據(jù)。 對(duì)采集的流場(chǎng)數(shù)據(jù)進(jìn)行DMD 分析。

圖5 為來流10 m/s、開裂角5°狀態(tài)下的原流場(chǎng)及前15 階DMD 模態(tài)重構(gòu)的流場(chǎng)壓強(qiáng)云圖,重構(gòu)流場(chǎng)與原流場(chǎng)基本一致,圖6 為DMD 分解得到的前4 階流場(chǎng)模態(tài),均為穩(wěn)定的周期模態(tài)。

圖5 流場(chǎng)壓強(qiáng)云圖Fig.5 Pressure contours of flow field

圖6 開裂角5°前4 階DMD 模態(tài)壓強(qiáng)云圖Fig.6 Pressure contours of the first four DMD modes at crack angle of 5°

各階模態(tài)的增長率與頻率數(shù)據(jù)如表1 所示。

表1 前5 階DMD 模態(tài)增長率與頻率Table 1 Frequency and magnification of the first five DMD modes

DMD 的1 階模態(tài)表示頻率為0 Hz 的平均流場(chǎng),2 階模態(tài)體現(xiàn)了阻力舵的上下后緣位置交替出現(xiàn)了脫落渦。 這2 個(gè)模態(tài)代表了流場(chǎng)的主要特征。 可以看出,3 ～5 階模態(tài)頻率大致為2 階模態(tài)的整數(shù)倍,這說明DMD 分解捕捉到了流場(chǎng)的主要結(jié)構(gòu),同時(shí)提取到了高階諧波模態(tài)。

圖7 展示了一個(gè)周期T內(nèi)開裂區(qū)的流場(chǎng)演化過程。 可以看到,在開裂區(qū)內(nèi)部共存在a、b、c這3 個(gè)渦,其中,a是一直存在于開裂區(qū)內(nèi)的駐渦,b、c分別為上下舵面后緣位置生成的脫落渦,它們分別與流場(chǎng)1、2 階模態(tài)對(duì)應(yīng)。 在0 ～T時(shí)段,渦b從上舵面后緣脫落進(jìn)入流場(chǎng),渦c由開裂區(qū)內(nèi)逐漸向下舵面后緣移動(dòng);時(shí)段,渦b離開后緣位置進(jìn)入流場(chǎng),渦c繼續(xù)向后緣移動(dòng),渦量逐漸增強(qiáng);時(shí)段,隨著渦c向后移動(dòng),上舵面逐漸生成了一個(gè)新的渦b;T～1T時(shí)段,渦b向后緣移動(dòng),同時(shí)新的渦c也在下舵面開始形成。

圖7 流場(chǎng)演化過程Fig.7 Flow field evolution

對(duì)不同來流速度下的流場(chǎng)數(shù)據(jù)進(jìn)行分析,結(jié)果發(fā)現(xiàn),當(dāng)來流速度在10 ～60 m/s 時(shí),前5 階流場(chǎng)模態(tài)云圖基本一致。 不同來流速度下的DMD模態(tài)頻率如圖8 所示。 可以看出,隨著來流速度增加,各階流場(chǎng)模態(tài)頻率基本線性增加,且高階模態(tài)頻率增加的更快。

圖8 不同來流速度下的DMD 模態(tài)頻率Fig.8 Frequencies of DMD modes versus velocity

在來流速度10 m/s 的狀態(tài)下,分別對(duì)5°、10°、20°開裂角時(shí)的DMD 流場(chǎng)模態(tài)分解結(jié)果進(jìn)行分析,頻率分解結(jié)果如圖9 所示。 開裂角的增大引起尾流區(qū)、回流區(qū)的范圍變大,但流場(chǎng)結(jié)構(gòu)沒有太大差別。

圖9 不同開裂角下的DMD 模態(tài)頻率Fig.9 Frequencies of DMD modes versus crack angles

DMD 分析的結(jié)果表明,在給定的速度范圍內(nèi),流場(chǎng)結(jié)構(gòu)主要表現(xiàn)為開裂區(qū)內(nèi)駐渦與后緣脫落形成的卡門渦的疊加。 由于卡門渦的頻率與速度成正比,與物體垂直于來流方向截面的特征尺度成反比,從圖8 中的數(shù)據(jù)可以看出,在固定的開裂角下,DMD 模態(tài)頻率與來流速度基本呈線性關(guān)系; 而圖9 則說明了在相同來流速度下,隨著開裂角的增大,翼型垂直來流方向的截面長度增加,因此流場(chǎng)模態(tài)頻率隨之下降;2 組結(jié)果反應(yīng)出流場(chǎng)頻率與流速和開裂角相關(guān),不同DMD 模態(tài)之間頻率的增長率不同。

4 流固耦合穩(wěn)定性分析

當(dāng)渦脫落的頻率降低至與結(jié)構(gòu)的固有頻率靠近時(shí),有可能會(huì)激發(fā)結(jié)構(gòu)的振動(dòng),引發(fā)流固耦合系統(tǒng)的失穩(wěn)。

使用Fluent 六自由度求解器進(jìn)行結(jié)構(gòu)求解計(jì)算,在模型的0.5 倍弦長位置安裝扭轉(zhuǎn)彈簧。 設(shè)置模型的慣性矩I為0.001 35 kg·m2,給定來流速度V為10 m/s,初始迎角為0.5°。 結(jié)構(gòu)固有頻率fn可由結(jié)構(gòu)慣性矩I與彈簧剛度K得出:

無阻尼系統(tǒng)俯仰運(yùn)動(dòng)方程可以表示為

式中:M為流場(chǎng)計(jì)算得到的氣動(dòng)力矩;¨α為俯仰位移。

4.1 5°開裂角阻力舵流固耦合響應(yīng)分析

圖10 為5°開裂角的模型在不同的結(jié)構(gòu)固有頻率下俯仰角及俯仰力矩隨時(shí)間變化的結(jié)果,其中f為頻域分析中的頻率。 從圖10 中的結(jié)果來看,結(jié)構(gòu)固有頻率為6.6 Hz 時(shí),結(jié)構(gòu)響應(yīng)的主要頻率約為67 Hz,α的運(yùn)動(dòng)幅值很小,約為0.02°。這與表1 中流場(chǎng)2 階DMD 模態(tài)的頻率結(jié)果基本一致,因此,可以判斷這是由于卡門渦的周期性脫落造成的渦致振動(dòng)。 而當(dāng)結(jié)構(gòu)的固有頻率為6.2 Hz時(shí),結(jié)構(gòu)的運(yùn)動(dòng)頻率約為2.5 Hz,α的運(yùn)動(dòng)幅值約為3°,渦致振動(dòng)的影響基本可以忽略,系統(tǒng)處于失穩(wěn)狀態(tài)。 此時(shí)結(jié)構(gòu)運(yùn)動(dòng)頻率與固有頻率不一致是由于翼型帶來的氣動(dòng)附加剛度導(dǎo)致的。

圖10 不同固有頻率下的響應(yīng)曲線(5°開裂角)Fig.10 Response curves with different fn values natural frequency(crack angle is 5°)

定義無量綱的折減速度U*[22]為

式中:V為來流速度;fn為結(jié)構(gòu)固有頻率;c為翼型弦長。

由此可以得到一系列計(jì)算狀態(tài)下系統(tǒng)的失穩(wěn)頻率fa隨U*的變化規(guī)律。

對(duì)多個(gè)狀態(tài)下的系統(tǒng)響應(yīng)曲線進(jìn)行分析,得到圖11 所示結(jié)果。

從圖11 中可以看出,模型的運(yùn)動(dòng)大致分為3個(gè)狀態(tài),當(dāng)U*較小時(shí),系統(tǒng)不會(huì)發(fā)生失穩(wěn),而是做小振幅的渦致振動(dòng),運(yùn)動(dòng)頻率與渦脫落頻率一致。 隨著U*增大至5.2 ～5.7 之間,結(jié)構(gòu)的運(yùn)動(dòng)形式表現(xiàn)為“淺失速”下的失速顫振,此時(shí)流固耦合系統(tǒng)在一定幅值下作俯仰振蕩。 由于氣動(dòng)剛度的影響,系統(tǒng)的運(yùn)動(dòng)頻率略低于結(jié)構(gòu)的固有頻率[23]。 當(dāng)U*繼續(xù)增大,系統(tǒng)進(jìn)入“深失速”顫振狀態(tài),此時(shí)的運(yùn)動(dòng)頻率約在6 Hz 左右。 這是由于系統(tǒng)無阻尼,同時(shí)結(jié)構(gòu)剛度比較小,運(yùn)動(dòng)幅值較大,此時(shí)的流場(chǎng)特性占主導(dǎo)。 由于來流速度不變,在運(yùn)動(dòng)過程中前緣渦以6 Hz 左右的頻率脫落、再附,結(jié)構(gòu)響應(yīng)頻率也在6 Hz 左右。

圖11 系統(tǒng)失穩(wěn)頻率隨折減速度變化(5°開裂角)Fig.11 Instability frequency variations with U*(crack angle is 5°)

4.2 20°開裂角阻力舵流固耦合響應(yīng)分析

當(dāng)開裂角為20°時(shí),運(yùn)動(dòng)規(guī)律如圖12 所示。該狀態(tài)下的運(yùn)動(dòng)與5°開裂角相似,結(jié)構(gòu)固有頻率為6.1 Hz 時(shí),運(yùn)動(dòng)形式為渦致振動(dòng),其頻率與流場(chǎng)DMD 的2 階模態(tài)一致,為30 Hz 左右,振動(dòng)幅值很小。 當(dāng)結(jié)構(gòu)固有頻率為5.3 Hz 時(shí),系統(tǒng)運(yùn)動(dòng)形式為流固耦合失穩(wěn),頻率大約為1 Hz,α的運(yùn)動(dòng)幅值在15°左右。

圖12 不同固有頻率下的響應(yīng)曲線(20°開裂角)Fig.12 Response curves with different natural frequency(crack angle is 20°)

圖13 展示了開裂角20°時(shí)fa隨U*的變化情況,其規(guī)律與5°開裂角基本一致。 但是系統(tǒng)發(fā)生失速顫振時(shí),折減速度的范圍在5.6 ～6.4 之間,高于5°開裂角的狀態(tài)。 這個(gè)現(xiàn)象說明,在阻力舵打開的過程中,系統(tǒng)的顫振邊界提高,也就意味著在大開裂角狀態(tài)下,阻力舵系統(tǒng)具有更好的穩(wěn)定性。

圖13 系統(tǒng)失穩(wěn)頻率隨折減速度變化(20°開裂角)Fig.13 Instability frequency variations with U*(crack angle is 20°)

5 結(jié) 論

1) 采用DMD 方法對(duì)二維阻力舵的流場(chǎng)特性進(jìn)行探究,結(jié)果表明二維阻力舵流場(chǎng)的主要特征為開裂區(qū)的駐渦與后緣脫落渦的疊加。 流場(chǎng)的DMD 模態(tài)頻率與來流速度成正相關(guān),與開裂角的大小成負(fù)相關(guān)。

2) 進(jìn)行阻力舵單自由度流固耦合運(yùn)動(dòng)仿真計(jì)算,結(jié)果表明,隨著折減速度的增加,系統(tǒng)運(yùn)動(dòng)的主要形式由渦致振動(dòng)發(fā)展為失速顫振,渦致振動(dòng)的幅值遠(yuǎn)遠(yuǎn)小于失速顫振。 顫振邊界隨著開裂角的增大而減小,這表明增大開裂角對(duì)提高系統(tǒng)的失穩(wěn)邊界有著積極作用。