劉旭博，于紀(jì)言

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

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

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

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

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

1 撲動函數(shù)與撲翼建模

1.1 撲動函數(shù)

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

圖1 雙曲柄連桿機構(gòu)

這里為了簡化模型，設(shè)置上下?lián)鋭臃翟?5°左右，且撲動對稱。由于撲動對稱，機構(gòu)沒有急回特性，因此機構(gòu)的撲動角度曲線可以近似成一個簡諧運動的函數(shù)曲線，另外根據(jù)文獻[8]，鳥類的撲動角度變化規(guī)律確實類似于一個三角函數(shù)的周期性變化。

撲翼角度函數(shù)的格式如下：

θ=bsin(2πft)

(1)

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

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

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

(2)

求導(dǎo)得到其角速度函數(shù)：

(3)

1.2 撲翼建模

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

表1 典型鴿子參數(shù)

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

圖2 鴿翅

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

圖3 翅翼模型

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

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

圖4 計算域

2 計算前處理

1) 計算方法及網(wǎng)格劃分策略

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

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

圖5 計算域網(wǎng)格

2) 計算模型

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

質(zhì)量守恒方程：

(4)

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

動量方程：

(5)

低雷諾數(shù)下的k-ε方程：

Gk+Gb-ρε-YM+Sk

(6)

(7)

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

(8)

C2=1.9

C3=1.0

3) 邊界條件與動網(wǎng)格

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

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

(9)

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

3 鴿狀撲翅氣動特性的數(shù)值計算

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

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

(10)

(11)

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

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

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

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

圖6 升力系數(shù)曲線

圖7 阻力系數(shù)曲線

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

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

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

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

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

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

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

(12)

(13)

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

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

圖11 不同攻角下的平均升力系數(shù)

圖12 不同攻角下的平均阻力系數(shù)

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

4 結(jié)論

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

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

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