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

李戰東,龔昌全，王 巍，陶建國

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

1 引言

撲翼飛行器由于其獨特的仿生外形和靈活的機動性，在軍事和商業領域具有廣泛的應用前景。撲翼飛行器的氣動外形和運動方案是影響其氣動性能的主要因素[1]。與傳統飛行器相比，撲翼飛行器的運動主要由上下撲動、弦向扭轉、展向折疊和前后揮擺幾種運動耦合而成，產生飛行所需的氣動升力和推力，從而使其具備靈活的機動性能[2]。因此，研究運動模式對撲翼氣動性能的影響具有重要意義。

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

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

2 數值計算模型和方法

2.1 計算模型

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

表1 三維撲翼的幾何和運動學參數

圖1 NACA0011三維撲翼模型

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

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

(1)

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

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

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

(2)

式中：θ2m撲翼弦向的扭轉幅值。

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

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

(3)

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

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

2.2 數值計算方法

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

圖2 流域網格劃分及邊界條件設置

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

圖3 不同網格精度下的氣動力變化曲線

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

(4)

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

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

2.3 氣動力方法驗證

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

圖4 本文的氣動力方法與Yang實驗數據對比

從圖4中可以看出，本文采用的氣動力計算方法計算結果與實驗數據吻合效果較好，部分數據存在波動，這是因為Yang等人的撲翼模型數據是利用DIC技術采集得到，而本文的撲翼計算模型是根據Yang等人的撲翼幾何參數所建立，導致本文的撲翼模型在翼面后緣部分與實驗模型有微小差別，引起了小幅度的波動。總體來看，本文的仿真計算結果符合真實樣機的氣動力實驗結果，這說明本文的撲翼氣動力計算方法準確性較高。可以準確計算撲翼的氣動性能。

3 數值結果與分析

本文主要研究扭轉運動和水平運動對撲翼升力方向(Z)與推力方向(X)的氣動性能影響，研究在不同撲動頻率以及不同自由度下撲翼的氣動性能和流場結構。

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

基于圖1的撲翼模型，本文將撲翼撲動運動起始行程定義為由下撲行程起點向下撲動，撲翼的扭轉運動起始行程為向下扭轉，撲翼的水平運動起始行程定義為向后(-X)運動，撲動運動與扭轉運動的頻率一致，即上撲行程時，撲翼向上撲動，翼面向上扭轉，撲翼也同時向前運動，而下撲行程時，撲翼向下撲動，翼面向下扭轉，撲翼同時向后運動。撲翼水平運動的頻率調控參數k取值為1，對不同自由度的升力和推力系數進行了數值計算，得到撲翼在一個周期內不同自由度下的升力和推力系數變化曲線如圖5和圖6所示。

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

其中0T-0.5T為下撲行程，0.5T-1T為上撲行程。從圖5中可以看出，當加入扭轉運動后，兩自由度的升力峰值A2相比單自由度升力峰值A1提升了28%，而3自由度與兩自由度相比，兩者的升力變化曲線幾乎重合，這說明撲翼的扭轉運動可有效提升撲翼升力性能，而水平運動對升力的產生幾乎沒有影響。

在推力性能方面，如圖6所示，在單自由度撲動時，撲翼在整個撲動行程中推力系數都處于正值，撲翼的推力系數峰值B1都出現在下撲行程的水平位置，即圖5中的A1時刻；而3自由度下的撲翼在下撲行程時，推力系數為正，在上撲行程時出現較小負推力系數。推力峰值B3和兩自由度B2較單自由度B1而言提前了0.12T，即出現在撲翼下撲至水平位置之前。

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

相對于單自由度撲翼的下撲行程，加入弦向扭轉運動后，隨著撲翼向下扭轉的角度逐漸增加，在B2時刻(0.12T)推力系數增大到峰值，與該時刻單自由度的推力系數B1相比，撲翼的推力性能大幅度提升了4.69倍。在兩自由度撲翼的基礎上加入水平運動后，即3自由度下的推力峰值B3得到了更進一步的提升，與B1相比，具備3個自由度的撲翼推力系數提升了6.26倍，顯著提高了撲翼的推力性能。

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

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

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

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

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

在推力性能方面，從圖8中可以看出，10 Hz下最大推力系數為78，10 Hz和5 Hz的推力峰值都出現在上撲行程結束后開始下撲行程的起始位置附近。撲動頻率5 Hz與10 Hz相比，10 Hz下撲翼的推力提高了4.26倍，且在上撲行程起始位置附近的推力提升更大。

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

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

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

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

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

綜合來看，加入扭轉運動會使撲翼的升力小幅度增加，推力性能大幅度增加；而再加入水平運動即3自由度運動后，進一步提升了撲翼的推力性能。隨著撲動頻率的增加，撲翼3自由度運動模式的升力性能和推力性能得到了更大幅度的提升，這表明撲動頻率對撲翼的氣動性能具有正相關的顯著影響。此外，隨著撲翼的水平運動頻率增加，撲翼的升力略有提高，但同時也會增加負升力，而在推力性能方面，水平運動頻率的增大會小幅度減小撲翼下撲行程的氣動推力性能，但會顯著提高撲翼在上撲行程的氣動推力性能。

4 氣動機理分析

為了詳細闡述3自由度下撲翼的氣動機理，本文將不同自由度的數據計算結果導入CFD-POST軟件進行后處理，截取了一個完整撲動周期的計算結果，通過撲翼展向壓力場分布情況和渦結構變化過程對不同自由度的氣動特性進行分析。

4.1 撲翼壓力場分析

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

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

4.2 撲翼渦結構分析

為了進一步揭示不同自由度下撲翼氣動力的產生機理，本文還將A時刻升力峰值和B時刻推力峰值的渦結構進行了提取，得到了圖12撲翼弦向0.4R處截面渦量云圖。在A時刻，從圖12可以看出，當撲翼下撲到水平位置時，在撲翼上表面產生了強烈前緣渦，前緣渦具有附著于撲翼表面不脫落的特性，附著于翼面的前緣渦面積越大，撲翼的翼面壓差也就越大，這也解釋了圖11中該時刻(0.5T)的翼面壓差達到最大，由此產生升力峰值的原因。此外，當添加扭轉運動后，撲翼上表面的前緣渦明顯變得更強，渦形也變得更加飽滿，而加入水平運動后，A3相對于與A2的前緣渦強度和范圍變化不大。

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

當撲翼在上下撲動的轉換階段，會在撲翼的后緣產生一個尾緣渦隨著空氣的流動向后脫落，受到尾緣渦的誘導作用，在撲翼后緣附近將產生一個與來流相同方向的流動區域，此時漩渦中各處的誘導速度與來流方向相同，撲翼會受到推力。反之，當漩渦中各處的誘導速度與來流方向相反時，撲翼會受到阻力。

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

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

5 結論

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

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

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

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

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