基于CFD仿真的撲動模式對撲翼氣動性能影響研究

李戰(zhàn)東,龔昌全，王 巍，陶建國

(1.沈陽航空航天大學(xué) 民用航空學(xué)院，沈陽 110136；2.哈爾濱工業(yè)大學(xué) 機電工程學(xué)院, 哈爾濱 150001)

1 引言

撲翼飛行器由于其獨特的仿生外形和靈活的機動性，在軍事和商業(yè)領(lǐng)域具有廣泛的應(yīng)用前景。撲翼飛行器的氣動外形和運動方案是影響其氣動性能的主要因素[1]。與傳統(tǒng)飛行器相比，撲翼飛行器的運動主要由上下?lián)鋭印⑾蚁蚺まD(zhuǎn)、展向折疊和前后揮擺幾種運動耦合而成，產(chǎn)生飛行所需的氣動升力和推力，從而使其具備靈活的機動性能[2]。因此，研究運動模式對撲翼氣動性能的影響具有重要意義。

撲翼飛行的運動模式對飛行過程中高升力和高推力的產(chǎn)生具有重大影響。近年來，國內(nèi)外的研究人員采用仿真計算和實驗，研究了撲翼飛行運動參數(shù)的變化規(guī)律，以探索撲翼的運動模式對其氣動性能的影響。對于剛性撲翼，研究人員通過將鳥類翅膀的上下?lián)鋭雍喕癁槎S翼型的沉浮運動,研究了運動參數(shù)對沉浮翼型產(chǎn)生推力大小的影響[3]，Yang等[4]將機翼置于二維層流中，通過改變扭轉(zhuǎn)和水平揮擺運動之間的相位差，發(fā)現(xiàn)當(dāng)扭轉(zhuǎn)角為30°，撲翼會產(chǎn)生高平均推力，而較小幅度的水平掃掠運動有利于提升撲翼的推進效率，且隨著揮擺運動幅度的增加，撲翼的升力也會隨之小幅度增加；Xijun Ke等[5]通過引入俯仰頻率與撲動頻率的頻率比值，研究了兩自由度三維撲翼俯仰角相對撲動相位偏移的可調(diào)規(guī)律，Khanh Nguyen等[6]采用計算流體動力學(xué)方法研究了機翼運動學(xué)對懸停氣動效率的影響，發(fā)現(xiàn)拍擊和甩打可以提高5%的升力；M.M De等[7]研究了非對稱1自由度運動撲翼模型在不同平面形狀下的氣動特性，發(fā)現(xiàn)橢圓形的翼面飛行效率更高；吳越等[8]采用了非定常渦格法(UVLM)計算撲翼氣動力，發(fā)現(xiàn)經(jīng)過優(yōu)化后的撲翼運動方案能夠使特定翼面發(fā)揮最佳的氣動性能；Sang-Hoon Yoon等[9]開發(fā)了一種三維撲翼流固耦合求解器，發(fā)現(xiàn)被動扭轉(zhuǎn)運動對弧面機翼氣動性能有顯著影響。謝鵬[10]基于ADAMS軟件，設(shè)計了一種折疊翼仿鳥撲動飛行器，并發(fā)現(xiàn)折疊翼在上撲階段受到更小的阻力，在撲動過程中能夠獲得更大的升力;Yang等[11]結(jié)合CFD數(shù)值計算方法，設(shè)計了一種具有完全自主操作的能力撲翼飛行器Dove，可實現(xiàn)半小時的巡航工作，并同時向4 km外的地面站發(fā)送實時穩(wěn)定的彩色視頻。撲翼的氣動性能與附著在機翼表面的渦結(jié)構(gòu)也密切相關(guān)，A.Boudis等[12]通過對二維兩自由度非正弦運動軌跡撲翼模型的渦脫落過程進行研究，發(fā)現(xiàn)相比于正弦撲動軌跡，非正弦軌跡每半個周期可以脫落2個反向渦，且增大了尾緣渦的強度，顯著提高了撲翼的推進效率。Chunlin Gong等[13]提出了一種基于浸沒邊界格子Boltzmann方法(IB-LBM)的大規(guī)模并行求解器，分析了三維撲翼前緣渦和尾流渦結(jié)構(gòu)的演化過程，發(fā)現(xiàn)隨著展弦比的增大，撲翼的推力系數(shù)先增大后減小。Deng等[14]搭建了六分力傳感器測力實驗分析了撲翼頻率和機翼柔性對氣動力的影響，然后使用粒子圖像測速技術(shù)揭示了撲翼后緣渦的脫落過程，并發(fā)現(xiàn)無翼肋約束的機翼比有翼肋約束的氣動性能更好。Yang等[15]通過數(shù)值模擬和實驗研究了單自由度撲翼樣機的氣動機理，發(fā)現(xiàn)撲翼在升力方向上的慣性力大小與氣動力的大小相同，而在推力方向上慣性力相對較小。

綜上所述，目前絕大部分的研究都是基于單自由度和兩自由度耦合運動模型展開，主要研究幾何和結(jié)構(gòu)參數(shù)設(shè)計對氣動性能的影響，并沒有完全揭示各個運動對氣動機理的影響。同時，大部分研究都是基于二維翼型開展，忽略了三維仿生撲翼流場強烈的非定常效應(yīng)、渦流動復(fù)雜的特點，沒有精確揭示運動模式對撲翼氣動性能的影響，而復(fù)雜渦系的精確捕捉是揭示撲翼飛行器產(chǎn)生高升力和大推力機理的重要基礎(chǔ)。本文建立了三維撲翼多自由度運動模型，采用數(shù)值模擬的方法將計算模型導(dǎo)入CFD求解器中，對撲翼的非定常N-S方程進行求解，獲得不同真實運動模式下翼的升力和推力數(shù)據(jù)，然后結(jié)合翼面壓力場分布情況和流場渦結(jié)構(gòu)的演變過程，對撲翼的升推力變化曲線、壓力場和渦結(jié)構(gòu)的結(jié)果進行分析討論，進一步揭示各個運動模式對翅翼氣動性能的影響。

2 數(shù)值計算模型和方法

2.1 計算模型

為了更真實的研究運動模式對撲翼氣動性能的影響，本文采用NACA00012翼型建立三維撲翼模型，翼面布局形式為反齊默爾曼機翼，研究其在單自由度和多自由度下的升力性能和推力性能問題。整個撲翼模型和運動坐標(biāo)系如圖1和表1所示。

表1 三維撲翼的幾何和運動學(xué)參數(shù)

圖1 NACA0011三維撲翼模型

翅膀的撲動方式和撲動規(guī)律對鳥類飛行過程中高升力和大推力的形成至關(guān)重要。典型的鳥類翅膀宏觀撲動包含撲動、扭轉(zhuǎn)、掃掠和折疊。本文將主要研究YZ方向的撲動、弦向XZ面的扭轉(zhuǎn)和水平XY面上的水平掃掠運動的耦合對撲翼氣動性能的影響。撲翼飛行器翅翼的上下?lián)鋭舆\動簡化為以X軸為旋轉(zhuǎn)軸作撲動的簡諧運動[4]，則翅翼的撲動規(guī)律方程為：

θ1(t)=θ1msin(2πft)

(1)

式中：θ1(t)為翅翼在垂直方向上的撲動角速度；f為撲翼運動的撲動頻率；θ1m為撲翼運動的撲動幅值。

對于撲翼的弦向扭轉(zhuǎn)運動，將其定義為以Y軸為旋轉(zhuǎn)中心的扭轉(zhuǎn)運動，簡化后的扭轉(zhuǎn)運動規(guī)律方程為：

θ2(t)=θ2 msin(2πft)

(2)

式中：θ2m撲翼弦向的扭轉(zhuǎn)幅值。

對于撲翼的水平運動，將其定義為撲翼沿X軸方向的水平運動，簡化后的水平運動規(guī)律方程為：

h(t)=hmcos(2kπft)

(3)

式中：h(t)為撲翼運動過程中，翅翼在水平方向上的位移；hm為撲翼運動的水平運動的幅值；k為撲翼水平運動的頻率調(diào)控參數(shù)。

本文將θ1(t)的撲動運動定義為單自由度運動；θ2(t)的扭轉(zhuǎn)運動和θ1(t)撲動的耦合運動定義為兩自由度運動；θ1(t)的撲動運動、θ2(t)的扭轉(zhuǎn)運動和h(t)水平運動共同控制的耦合運動定義為3自由度運動。

2.2 數(shù)值計算方法

基于圖1的三維撲翼模型，在開展撲翼運動數(shù)值計算模擬時，確定翅翼各個自由度的運動參數(shù)如下：翅翼的撲動幅值為120°，扭轉(zhuǎn)幅值為30°，扭轉(zhuǎn)中心為1/5c，撲翼運動頻率為f=5 Hz，水平運動的幅值為0.05c?？紤]到流動對氣動性能的影響，計算域的大小為13c×8c×8c，邊界條件的設(shè)置如圖2所示。網(wǎng)格劃分由ICEM CFD軟件完成，采用四面體非結(jié)構(gòu)網(wǎng)格來離散流體域，并對近機翼表面進行加密處理，通過計算網(wǎng)格總數(shù)量為118萬，最差網(wǎng)格質(zhì)量為0.35，滿足計算精度要求。

圖2 流域網(wǎng)格劃分及邊界條件設(shè)置

隨著網(wǎng)格精度的提高，研究計算結(jié)果的收斂性是至關(guān)重要的，因此，通過改變計算域網(wǎng)格的數(shù)量，對撲翼的計算結(jié)果進行了網(wǎng)格無關(guān)性驗證，為此還考慮了網(wǎng)格數(shù)量為57萬和271萬的計算結(jié)果，得到一個周期內(nèi)不同網(wǎng)格密度隨時間的升力變化曲線如圖3所示。從圖3中可以看出，采用網(wǎng)格數(shù)量為57萬和271萬所得的升力數(shù)據(jù)計算結(jié)果變化曲線與網(wǎng)格數(shù)量為118萬的計算得到了較好的吻合，這表明本文的撲翼數(shù)值計算結(jié)果基本不受網(wǎng)格密度變化的影響。

圖3 不同網(wǎng)格精度下的氣動力變化曲線

湍流模型采用Realizableε-epsilon模型，撲翼的流體控制方程采用有限體積法進行求解，其中控制方程中非定常項采用1階隱式格式進行離散，對流項采用2階迎風(fēng)格式進行離散，擴散項采用中心差分格式進行離散，速度和壓力的解耦采用SIMPLEC算法實現(xiàn)，然后結(jié)合網(wǎng)格光順和重構(gòu)的動網(wǎng)格技術(shù)，分別對不同自由度的撲翼運動模型進行仿真計算，得到撲翼的升力變化曲線和推力系數(shù)變化曲線。撲翼的推力系數(shù)定義為：

(4)

式中：FT為撲翼推力；ρ為撲翼周圍空氣密度；S為撲翼面積。

數(shù)值計算中，采用有限體積法，通過ANSYS FLUENT軟件對上述的控制方程進行數(shù)值仿真。撲翼模型的運動通過用戶自定義方程(UDF)控制實現(xiàn)。

2.3 氣動力方法驗證

為了驗證本文氣動力計算方法的準(zhǔn)確性，本文根據(jù)Yang[15]對撲翼樣機Dove的實驗數(shù)據(jù)，建立了幾何參數(shù)和運動學(xué)參數(shù)一致的撲翼模型，其中翼展R為300 mm，弦長c為100 mm，撲動角度為70°，撲動頻率為6 Hz；運用本文的氣動力計算方法和運動模型對撲翼的氣動性能進行仿真計算，得到對應(yīng)撲動行程的實驗數(shù)據(jù)對比結(jié)果如圖4所示。

圖4 本文的氣動力方法與Yang實驗數(shù)據(jù)對比

從圖4中可以看出，本文采用的氣動力計算方法計算結(jié)果與實驗數(shù)據(jù)吻合效果較好，部分?jǐn)?shù)據(jù)存在波動，這是因為Yang等人的撲翼模型數(shù)據(jù)是利用DIC技術(shù)采集得到，而本文的撲翼計算模型是根據(jù)Yang等人的撲翼幾何參數(shù)所建立，導(dǎo)致本文的撲翼模型在翼面后緣部分與實驗?zāi)Ｐ陀形⑿〔顒e，引起了小幅度的波動?？傮w來看，本文的仿真計算結(jié)果符合真實樣機的氣動力實驗結(jié)果，這說明本文的撲翼氣動力計算方法準(zhǔn)確性較高?？梢詼?zhǔn)確計算撲翼的氣動性能。

3 數(shù)值結(jié)果與分析

本文主要研究扭轉(zhuǎn)運動和水平運動對撲翼升力方向(Z)與推力方向(X)的氣動性能影響，研究在不同撲動頻率以及不同自由度下?lián)湟淼臍鈩有阅芎土鲌鼋Y(jié)構(gòu)。

3.1 不同自由度對氣動性能的影響

基于圖1的撲翼模型，本文將撲翼撲動運動起始行程定義為由下?lián)湫谐唐瘘c向下?lián)鋭樱瑩湟淼呐まD(zhuǎn)運動起始行程為向下扭轉(zhuǎn)，撲翼的水平運動起始行程定義為向后(-X)運動，撲動運動與扭轉(zhuǎn)運動的頻率一致，即上撲行程時，撲翼向上撲動，翼面向上扭轉(zhuǎn)，撲翼也同時向前運動，而下?lián)湫谐虝r，撲翼向下?lián)鋭?，翼面向下扭轉(zhuǎn)，撲翼同時向后運動。撲翼水平運動的頻率調(diào)控參數(shù)k取值為1，對不同自由度的升力和推力系數(shù)進行了數(shù)值計算，得到撲翼在一個周期內(nèi)不同自由度下的升力和推力系數(shù)變化曲線如圖5和圖6所示。

圖5 不同自由度的升力變化曲線

其中0T-0.5T為下?lián)湫谐蹋?.5T-1T為上撲行程。從圖5中可以看出，當(dāng)加入扭轉(zhuǎn)運動后，兩自由度的升力峰值A(chǔ)2相比單自由度升力峰值A(chǔ)1提升了28%，而3自由度與兩自由度相比，兩者的升力變化曲線幾乎重合，這說明撲翼的扭轉(zhuǎn)運動可有效提升撲翼升力性能，而水平運動對升力的產(chǎn)生幾乎沒有影響。

在推力性能方面，如圖6所示，在單自由度撲動時，撲翼在整個撲動行程中推力系數(shù)都處于正值，撲翼的推力系數(shù)峰值B1都出現(xiàn)在下?lián)湫谐痰乃轿恢茫磮D5中的A1時刻；而3自由度下的撲翼在下?lián)湫谐虝r，推力系數(shù)為正，在上撲行程時出現(xiàn)較小負(fù)推力系數(shù)。推力峰值B3和兩自由度B2較單自由度B1而言提前了0.12T，即出現(xiàn)在撲翼下?lián)渲了轿恢弥啊?/p>

圖6 撲動頻率為5 Hz下不同自由度的推力系數(shù)曲線

相對于單自由度撲翼的下?lián)湫谐蹋尤胂蚁蚺まD(zhuǎn)運動后，隨著撲翼向下扭轉(zhuǎn)的角度逐漸增加，在B2時刻(0.12T)推力系數(shù)增大到峰值，與該時刻單自由度的推力系數(shù)B1相比，撲翼的推力性能大幅度提升了4.69倍。在兩自由度撲翼的基礎(chǔ)上加入水平運動后，即3自由度下的推力峰值B3得到了更進一步的提升，與B1相比，具備3個自由度的撲翼推力系數(shù)提升了6.26倍，顯著提高了撲翼的推力性能。

3.2 不同撲動頻率對氣動性能的影響

在撲翼運動中，機翼的撲動頻率也是影響升力和推力性能大小的重要因素，較大頻率的運動不僅改變了撲翼附近的流場速度，同時還改變了撲翼表面的氣壓，從而改變氣流在整個翼展的分布情況。因此，本文還研究了撲翼在3自由度下不同撲動頻率的氣動特性。與上文的計算方法一致，本文計算了撲動頻率為10 Hz下的氣動性能，得到了不同撲動頻率下?lián)湟淼纳屯屏ο禂?shù)變化曲線如圖7和圖8所示。

圖7 3自由度下10 Hz和5 Hz的升力變化曲線

圖8 3自由度下10 Hz和5 Hz的推力系數(shù)變化曲線

從圖7中可以看出，撲動頻率為10 Hz時最大升力為0.224 N，5 Hz下最大升力為0.055 1 N，相比之下，在3自由度撲動模式下，撲動頻率10 Hz下升力比5 Hz提高了3倍。

在推力性能方面，從圖8中可以看出，10 Hz下最大推力系數(shù)為78，10 Hz和5 Hz的推力峰值都出現(xiàn)在上撲行程結(jié)束后開始下?lián)湫谐痰钠鹗嘉恢酶浇?。撲動頻率5 Hz與10 Hz相比，10 Hz下?lián)湟淼耐屏μ岣吡?.26倍，且在上撲行程起始位置附近的推力提升更大。

3.3 不同運動頻率對氣動性能的影響分析

為了進一步探索撲動模式對撲翼氣動性能的影響的力學(xué)，本文繼續(xù)研究了撲翼相對運動頻率對氣動性能的影響。在水平前飛工況下的3自由度撲翼基礎(chǔ)上，對水平揮擺運動引入了運動頻率調(diào)控參數(shù)k，分別計算了當(dāng)k=1、2和3時的氣動力參數(shù)，得到3自由度下不同運動頻率的升力和推力系數(shù)變化曲線如圖9和圖10所示。從圖9中可以看出，在下?lián)湫谐讨?，相對于k=1的常規(guī)3自由度運動而言，當(dāng)k=2時，撲翼的升力峰值增加了6%，隨著水平運動頻率增大，當(dāng)k=3時，撲翼的升力峰值減小了4%，而在上撲行程中，撲翼的負(fù)升力峰值隨著水平運動頻率的增大逐漸增大，當(dāng)k=2時，撲翼的負(fù)升力峰值增大了3.6%，當(dāng)k=3時，撲翼的負(fù)升力峰值增大了6%。

圖9 3自由度下不同水平運動頻率的升力變化曲線

圖10 3自由度下不同運動頻率的推力系數(shù)變化曲線

在推力性能方面，在下?lián)湫谐讨?，撲翼的推力系?shù)隨著k值的增加略有降低，而相對于k=1而言，k=2和k=3的推力峰值接近，推力峰值均降低了7%；在上撲行程中，隨著水平運動頻率的增大，撲翼的推力系數(shù)逐漸增加，當(dāng)k=2時，撲翼的推力系數(shù)峰值增加了17.8%，當(dāng)k=3時，撲翼的推力系數(shù)峰值增加了13%。

綜合來看，加入扭轉(zhuǎn)運動會使撲翼的升力小幅度增加，推力性能大幅度增加；而再加入水平運動即3自由度運動后，進一步提升了撲翼的推力性能。隨著撲動頻率的增加，撲翼3自由度運動模式的升力性能和推力性能得到了更大幅度的提升，這表明撲動頻率對撲翼的氣動性能具有正相關(guān)的顯著影響。此外，隨著撲翼的水平運動頻率增加，撲翼的升力略有提高，但同時也會增加負(fù)升力，而在推力性能方面，水平運動頻率的增大會小幅度減小撲翼下?lián)湫谐痰臍鈩油屏π阅?，但會顯著提高撲翼在上撲行程的氣動推力性能。

4 氣動機理分析

為了詳細(xì)闡述3自由度下?lián)湟淼臍鈩訖C理，本文將不同自由度的數(shù)據(jù)計算結(jié)果導(dǎo)入CFD-POST軟件進行后處理，截取了一個完整撲動周期的計算結(jié)果，通過撲翼展向壓力場分布情況和渦結(jié)構(gòu)變化過程對不同自由度的氣動特性進行分析。

4.1 撲翼壓力場分析

通過對數(shù)值計算結(jié)果的處理，得到撲翼YZ截面0.12c處不同拍動行程壓力云圖如圖11所示。

圖11 撲翼YZ截面不同拍動行程壓力云圖

4.2 撲翼渦結(jié)構(gòu)分析

為了進一步揭示不同自由度下?lián)湟須鈩恿Φ漠a(chǎn)生機理，本文還將A時刻升力峰值和B時刻推力峰值的渦結(jié)構(gòu)進行了提取，得到了圖12撲翼弦向0.4R處截面渦量云圖。在A時刻，從圖12可以看出，當(dāng)撲翼下?lián)涞剿轿恢脮r，在撲翼上表面產(chǎn)生了強烈前緣渦，前緣渦具有附著于撲翼表面不脫落的特性，附著于翼面的前緣渦面積越大，撲翼的翼面壓差也就越大，這也解釋了圖11中該時刻(0.5T)的翼面壓差達(dá)到最大，由此產(chǎn)生升力峰值的原因。此外，當(dāng)添加扭轉(zhuǎn)運動后，撲翼上表面的前緣渦明顯變得更強，渦形也變得更加飽滿，而加入水平運動后，A3相對于與A2的前緣渦強度和范圍變化不大。

圖12 撲翼弦向0.4R截面渦量云圖

當(dāng)撲翼在上下?lián)鋭拥霓D(zhuǎn)換階段，會在撲翼的后緣產(chǎn)生一個尾緣渦隨著空氣的流動向后脫落，受到尾緣渦的誘導(dǎo)作用，在撲翼后緣附近將產(chǎn)生一個與來流相同方向的流動區(qū)域，此時漩渦中各處的誘導(dǎo)速度與來流方向相同，撲翼會受到推力。反之，當(dāng)漩渦中各處的誘導(dǎo)速度與來流方向相反時，撲翼會受到阻力。

對于B時刻的渦量云圖，通過對比圖12(a)、(b)和(c)可知，在單自由度撲翼B1時刻，機翼后緣處行成了尾緣渦，準(zhǔn)備開始脫落行程，而在多自由度撲翼的B2和B3時刻，撲翼下表面的尾緣渦已完全脫落，且3自由度撲翼脫落的尾緣渦的形狀要比兩自由度更大，渦的強度也明顯強于前者。

綜上所述，通過比較升力和推力系數(shù)曲線，可以發(fā)現(xiàn)撲翼的扭轉(zhuǎn)運動會擴大前緣渦的形狀和增大前緣渦的強度，從而提高撲翼的升力性能。此外，扭轉(zhuǎn)運動還會加快撲翼下表面尾緣渦的脫落速度，使得多自由度撲翼提前達(dá)到推力峰值，且脫落的尾緣渦形狀更大，強度更大，這也是撲翼的推力性能得到大幅度提升的主要原因；而水平運動對前緣渦的形成幾乎沒有影響，但是撲翼加入水平運動使得撲翼脫落的尾緣渦強度獲得小幅度提升，形狀變大，致使3自由度撲翼的推力性能得到了進一步提升。

5 結(jié)論

針對扭轉(zhuǎn)和水平揮擺運動對三維撲動模型的氣動性能問題，基于FLUENT建立了三維撲翼氣動力仿真模型，計算了不同自由度和不同頻率下的升力、推力性能，并結(jié)合一個周期內(nèi)不同撲動行程的壓力云圖和渦結(jié)構(gòu)變化揭示了運動模式對三維撲翼氣動性能影響。研究表明：

1)本文建立的三維撲翼氣動力計算模型和計算方法與實驗測試結(jié)果吻合度較高，可準(zhǔn)確計算撲翼的氣動性能。

2)與單自由度撲翼相比，在迎風(fēng)水平前飛時，撲翼加入弦向扭轉(zhuǎn)運動(兩自由度)使撲翼的升力提升了25%，推力性能提升了4.69倍；而加入水平揮擺運動(3自由度)對撲翼的升力性能影響微小，但是使推力性能提升了6.26倍，且在不同的撲動頻率下，多自由度運動對撲翼氣動性能的提升效果更加明顯。

3)相對于常規(guī)3自由度撲翼運動而言，水平運動頻率調(diào)控參數(shù)k對撲翼氣動性能的影響比較明顯，隨著水平運動頻率的增加，撲翼的氣動升力先增大后減小，在下?lián)湫谐讨?，隨著k值的增大，撲翼的推力性能逐漸減小，而在上撲行程中，撲翼的推力性能先后增大后減小。當(dāng)k=2時，撲翼的升力峰值增加了6%，撲翼的推力系數(shù)峰值增加了17.8%，為取得較好的氣動性能，k=2的3自由度撲翼是一個比較好的選擇。

4)撲翼的弦向扭轉(zhuǎn)運動加速了撲翼下表面的前緣渦脫落進程，改變了前緣渦的脫落強度和形狀，增大了撲翼后緣與來流相同方向的流動區(qū)域強度，從而提高了推力性能；而水平揮擺運動會使兩自由度撲翼下表面脫落的尾緣渦強度得到進一步提高，脫落的尾緣渦形狀更大，從而使撲翼獲得更佳的氣動性能。