鴿狀撲翼飛行器氣動特性研究

劉旭博，于紀言

(南京理工大學 智能彈藥國防重點實驗室, 南京 210094)

撲翼飛行器是一種模擬鳥類飛行的仿生飛行器，在無人機發展中占有重要位置且應用需求巨大[1]，而相比于固定翼飛行器和旋翼飛行器，撲翼飛行器具有更高的飛行效率與機動性，因此撲翼飛行器也是近年來的研究熱點。

國外研究撲翼機起步較早，最早在1909年Knoller和Betz[2]發現二維翼型在撲動過程中可以同時產生升力和推力。1975年Weis-Fogh[3]解釋了昆蟲升力產生的拍合機制。1997年，Van Den Berg等[4]在對飛蛾飛行的研究中發現了“延時失速”機制，并在后來被Dickinso等[5]在果蠅模擬實驗中證實。近十年來，國外漸漸入手將理論轉化到實踐，較為典型的是在珠海航展上露面的smart bird，該撲翼機是由德國Festo公司在研究銀鷗的基礎上仿生制作的，是現代人類較為成功的仿生撲翼機。

國內對撲翼機的研究起步較晚，集中于近20年，但也取得了不少成就。曾銳與昂海松[6]對撲翼的撲動模型做了相關研究并指出了撲動頻率與撲動幅度對平均升力系數的影響。王姝歆[8]著重研究了仿生翅翼的中的尺度率問題。孫茂等[9]對蝴蝶與蜜蜂的飛行力學進行了相關的理論研究。另外，我國對撲翼的研究主要體現在實驗分析上，2006年至今，多位學者都對針對撲翼飛行器進行了相關的風洞實驗研究，并與不少關于撲翼機的理論研究相印證[10-13]。

盡管人類對撲翼機研究已久，然而現代撲翼機的飛行效率仍然遠遠小于鳥類的飛行效率，要想達到自然界鳥類的飛行效率，顯然仍然需要更多的研究積累。本文在前人對鳥類撲翼結構研究的基礎上，使用計算流體力學方法，采用三維鴿狀有限翼展撲翅模型研究翅翼拍動時的氣動特性，并對比不同攻角下的翅翼升阻力情況。

1 撲動函數與撲翼建模

1.1 撲動函數

在仿真之前需要先得到撲翼的撲動規律，撲翼的結構研究在近年來已趨于成熟，其中公認的易于實現、撲動過程對稱且振動較小的撲翼機構為雙曲柄連桿機構，如圖1：

圖1 雙曲柄連桿機構

這里為了簡化模型，設置上下撲動幅值在45°左右，且撲動對稱。由于撲動對稱，機構沒有急回特性，因此機構的撲動角度曲線可以近似成一個簡諧運動的函數曲線，另外根據文獻[8]，鳥類的撲動角度變化規律確實類似于一個三角函數的周期性變化。

撲翼角度函數的格式如下：

θ=bsin(2πft)

(1)

其中:θ為翅翼與水平面的夾角，b幅值，即翅翼撲動時翅翼能達到的極限位置，f為撲翅的頻率。t為撲動時間。

由于撲翼角度函數控制著翅翼的撲動規律，因此各參數對氣動特性有著明顯的影響。根據文獻6和文獻7，一個翼型的撲動幅值b越大，其周期內平均升力系數就越大，在0～60°范圍內基本呈線性相關。撲動頻率f對氣動特性的影響和翅翼翼型有關，基本規律是：升力系數在撲動頻率的某一點會達到峰值，在峰值前升力系數隨著撲動頻率的上升而增大，在峰值后升力系數隨著撲動頻率的上升而減小。

以鴿子為撲翼飛行器的仿生對象，通過觀察鴿子飛行的高速攝影錄像可知，鴿子的撲動頻率在6～10 Hz左右，具體設置函數時，采用8 Hz的撲動頻率。以θ表示翅翼與水平面的夾角，時間t為自變量，初始位置時翅翼水平，則撲翼角度函數如下：

(2)

求導得到其角速度函數：

(3)

1.2 撲翼建模

該撲翼機翅翼以鴿子為仿生對象，典型鴿子的飛行參數與翅翼特征參數數據如表1。

表1 典型鴿子參數

根據以上參數與具體鴿子的翅翼形狀(如圖2)，進行翅翼的建模：

圖2 鴿翅

根據鴿翅的生物學構造，鴿翅的前緣部分主要由肱骨支撐，量得肱骨直徑約為0.01 m。因此仿生翅的前緣半徑選擇為0.005 m。翅翼的后緣部分為羽毛尖端，厚度可以忽略不計，因此構建模型時后緣厚度為0。鴿翅的厚度最大處位于肱骨與副羽的相接部分，約在弦向20%的位置，厚度為0.016 m。翼展前緣長度約0.13 m，后緣長度為0.3 m。鴿翅的下表面展開時沒有明顯的弧度，為平面型下表面。鴿翅的翅翼模型如圖3，仿生翅模型的具體參數如表2所示。

圖3 翅翼模型

參數符號數值弦長c0.1 m展長b0.3 m前緣半徑r0.005 m最大厚度t0.016 m最大厚度位置p20%前緣長度q0.13 m后緣長度h0.3 m展弦比λ3

對翅翼建模完成后，構建外流域，外流域大小與流域流體的流速有關，由于是低速飛行，雷諾數較小，流場變化范圍不大，因此取外流域的長寬大約是翅翼的6倍即可。構建后的外流域如圖4：

圖4 計算域

2 計算前處理

1) 計算方法及網格劃分策略

Fluent在進行流場計算時，主要采用了有限體積法。有限體積法的基本思路是：將計算區域劃分為一系列不重復的控制體積，并使每個網格點周圍有一個控制體積；將待解的微分方程對每一個控制體積積分，便得出一組離散方程。

在CFD計算過程中，翅翼形狀相對復雜，且翅翼是不斷撲扇的，因此這里采用四面體網格，四面體網格有較好的復雜模型適應性，因此，本文采用四面體網格。在大尺度的變形下，計算迭代中會重構運動邊界附近的網格，從而適應部件新位置的流場運算。使用商業軟件ICEM CFD對：流動區域進行離散，最終生成的網格質量大于0.3，如圖5所示。

圖5 計算域網格

2) 計算模型

由于研究的是撲翼運動，氣體流速較慢，氣體密度基本無變化，采用常數值。流場以湍流為主，采用低雷諾數下的k-e方程進行計算。計算時主要使用的方程如下：

質量守恒方程：

(4)

式中：ρ為氣體密度，x,y,z為坐標，u,v,w為速度矢量，t為時間。

動量方程：

(5)

低雷諾數下的k-ε方程：

Gk+Gb-ρε-YM+Sk

(6)

(7)

其中ρ為流體密度，k為湍動能，ε為湍流耗散率，Gk是層流速度梯度而產生的湍流動能，Gb是由浮力而產生的湍流動能，YM是在可壓縮湍流中由于過渡的擴散而產生的波動。αk和αε是兩個方程的普朗特數，Sk和Sε是在計算中自定義的參數。另外其中的數值C1、C2、C3具體如下：

(8)

C2=1.9

C3=1.0

3) 邊界條件與動網格

為了模仿鳥類飛行時空氣流向，設置外流域正面及周邊四個面為速度入口邊界，流速為11 m/s，方向為x方向。設置翅翼后方的一個面為壓力出口流場，在有一個大氣壓的環境壓力情況下，壓力出口設置為0。翅翼部分默認為wall邊界。

動網格是本文研究中的一個重要手段，將1.1節中得到的角速度函數寫入UDF文件，而后在Fluent中編譯。由于網格變形較大，則同時采用動網格中的Remeshing和Smoothing的方法對網格進行更新。在最終選擇時間步長時應該考慮到網格特征尺寸與流速的影響，否則容易遇到更新網格失敗的情況，從而導致計算中斷。最終時間步長大小應該基本遵守以下公式：

(9)

其中:d為網格特征尺寸，v為流體流速。在以上公式計算的基礎上也可以適當調整從而節省時間與資源。在本文的計算過程中，時間步長取0.000 05 s。

3 鴿狀撲翅氣動特性的數值計算

3.1 翅翼拍動升、阻力變化特點

在Fluent中設置監測時，一般得到的是升阻力系數Cl、Cd。無量綱量升阻力系數就可以直接表征升阻力的大小。升阻力系數的定義如下：

(10)

(11)

其中:L為升力，D為阻力，ρ為氣體密度，v為氣體特征流速。

Fluent中設置氣體來流速度11 m/s，設置攻角為8°，可以監測到升阻力系數在周期內的變化情況如下：

由圖6、圖7可以看出，隨著翅翼的周期性撲動，升力與阻力也隨著時間出現周期性的變化，在翅膀下撲的過程中，當翅翼在水平位置時，翅翼的面積在豎直方向上的投影達到最大值，此時豎直方向上，即升力為最大值，當翅翼向上撲動復位時，翅翼與氣流的受力是向下的，因此產生了向下的力，則升力為負值。同理當翅翼的撲扇速度達到最大值時，翅翼與氣流的相互作用最強，此時在翅翼上的壓力也達到最大值，因此即使翅翼在豎直面上的投影面積不變，但由于作用于翅翼上的壓力的變化，阻力也因此發生相關的變化。繪制不同時刻下翅翼上下表面的壓力云圖，如圖8所示

從圖8可以看出，當翅翼撲翅速度在最大值附近時(t=0.005 s)，翅翼下翼面壓力遠遠大于上翼面，因此產生較大的升力。當翅翼接近極限位置時(t=0.030 s)，下翼面壓力仍然較上翼面大，但此時上翼面的壓力值與下翼面壓力值已經比較接近，從而產生的升力較小。當翅翼的撲動方向反轉時(t=0.045 s、0.060 s)，則上翼面壓力比下翼面大，但壓力差并不顯著，因此產生負的升力，但負升力相對產生的正升力而言并不大。為了研究撲翼的升力產生機理，對氣動渦的研究意義非凡。翅翼在撲扇過程中展向某平面渦，如圖9所示。

圖6 升力系數曲線

圖7 阻力系數曲線

圖8 不同時刻的翅翼壓力云圖

圖9 不同時刻x=10 cm處截面的速度矢量圖

由圖9可以觀察出在平行于YZ軸的某個截面處渦的發展情況，當撲翅在水平位置附近時(t=0.005 s)，由于翅翼上下壓差，造成了氣流從下翼面繞流到上翼面，從而產生了翼尖渦，當翅翼繼續下拍時(t=0.010 s、0.015 s)，上下翼面的壓差逐漸減小，渦強度也逐漸減弱，且渦從翼尖向翼根擴散，當翅翼達到極限位置附近時(t=0.025 s)，翼尖渦逐漸破裂消失。下面從三維空間具體研究翅翼周圍整個流場的渦的狀態：

由圖10可以看出，在翅翼剛開始撲扇時(t=0.005 s)，翅翼左側邊緣部分產生了前緣渦，翼尖產生一定強度的翼尖渦，并且在翅翼的邊緣前緣渦與翼尖渦匯合，新匯合的渦朝著翼根出發展并漸漸減弱破裂。隨著翅翼的繼續運動(t=0.010 s、0.015 s)，前緣渦逐漸減弱并朝著弦向發展。當翅翼撲扇到極限位置時，前緣渦強度已經非常小(t=0.025 s)，前緣渦與翼尖渦形成的新的匯合渦也逐漸脫離翼展并減弱。

圖10 不同時刻的翅翼三維渦核圖

3.2 不同攻角對升力、阻力的影響

鳥類在飛行過程中有往往有不同的飛行姿態，比較典型的姿態變化便是翅翼與來流方向的夾角的變化。本文借助Fluent，對不同攻角下的翅翼的升力與阻力做了數值計算，但在表征翅翼的升力阻力能力時，一般采用周期內的平均升阻力系數作為衡量標準，具體計算方法如下：

(12)

(13)

使用以上公式，對結果進行近似計算，最終得到-8°、0°、8°、16°、攻角時的翅翼周期內平均升阻力系數，如圖11和圖12：

由圖11、圖12可以看出，撲翼產生的升力隨著攻角的變大，不斷增大，攻角為負值時，升力也為負值。撲翼產生的阻力在0°附近時最小，之后隨著攻角的增大，阻力也增加。但不同的是，升力隨著攻角的增加基本上是線性的，而阻力隨著攻角的變化規律性并不強。

圖11 不同攻角下的平均升力系數

圖12 不同攻角下的平均阻力系數

以上的升阻力系數的變化與鳥類不同時刻的飛行姿態相印證，例如，當鳥類起飛時，身體總是與飛行方向有一個不小的角度差，這里是因為起飛時需要較大的升力，當鳥類在天空中快速飛行時，身體總是與飛行方向保持平行，則主要是為了獲得較小的阻力。許多鳥類在雨天快速接近地面捕捉昆蟲時，總是將身體下傾，此時升力為負，快速下降，將自身的勢能轉化為動能，通過調整翅翼角度快速掠過地面。

4 結論

1) 仿生撲翼飛行器飛行過程中的升阻力隨時間進行周期性變化，一個周期內影響升阻力變化的原因主要有撲動速度與翅翼所處位置在升阻力方向的面積投影。

2) 當翅翼在水平位置附近時，前緣渦強度較大，隨后朝著翼根處方向發展并擴散，當達到極限位置附近時，前緣渦與翼尖渦匯合，并逐漸脫離翼展。

3) 不同攻角下，在一定范圍內，翅翼的平均升力隨攻角的變化呈線性上升趨勢，但阻力的變化較為復雜，在攻角為0°時，阻力最小。另外得出的規律也與鳥類不同的飛行姿態相印證，解釋了鳥類起飛、掠飛時升阻力的具體情況。