劉旭博，于紀(jì)言

(南京理工大學(xué) 智能彈藥國(guó)防重點(diǎn)實(shí)驗(yàn)室, 南京 210094)

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

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

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

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

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

1.1 撲動(dòng)函數(shù)

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

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

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

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

θ=bsin(2πft)

(1)

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

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

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

(2)

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

(3)

1.2 撲翼建模

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

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

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

圖2 鴿翅

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

圖3 翅翼模型

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

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

圖4 計(jì)算域

2 計(jì)算前處理

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

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

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

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

2) 計(jì)算模型

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

質(zhì)量守恒方程：

(4)

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

動(dòng)量方程：

(5)

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

Gk+Gb-ρε-YM+Sk

(6)

(7)

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

(8)

C2=1.9

C3=1.0

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

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

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

(9)

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

3 鴿狀撲翅氣動(dòng)特性的數(shù)值計(jì)算

3.1 翅翼拍動(dòng)升、阻力變化特點(diǎn)

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

(10)

(11)

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

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

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

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

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

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

圖8 不同時(shí)刻的翅翼壓力云圖

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

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

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

圖10 不同時(shí)刻的翅翼三維渦核圖

3.2 不同攻角對(duì)升力、阻力的影響

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

(12)

(13)

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

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

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

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

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

4 結(jié)論

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

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

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