兩自由度微撲翼飛行器升阻力系數(shù)分析

李德寶， 舒寧， 王洪洋

（合肥工業(yè)大學(xué)機(jī)械與汽車工程學(xué)院，合肥230009）

0 引言

微型撲翼飛行器由于運(yùn)動(dòng)方式獨(dú)特、靈活性強(qiáng)等特點(diǎn)，越來(lái)越受到關(guān)注。在數(shù)值理論研究方面，研究人員［1］分析了二維機(jī)翼在一定頻率范圍內(nèi)推力產(chǎn)生的機(jī)理，研究發(fā)現(xiàn)在0.7 Hz附近推進(jìn)效果最佳；研究人員［2］分析研究了低雷諾數(shù)下，昆蟲(chóng)翅翼做非定常運(yùn)動(dòng)的氣動(dòng)性能，發(fā)現(xiàn)當(dāng)翅翼大攻角撲動(dòng)時(shí)，對(duì)大升力的產(chǎn)生十分有利；研究人員［3］通過(guò)數(shù)值模擬，對(duì)昆蟲(chóng)懸停時(shí)高升力機(jī)理有了更為深入的解釋。計(jì)算流體力學(xué)（簡(jiǎn)稱CFD）分析軟件Fluent是通過(guò)數(shù)值模擬計(jì)算，對(duì)研究流體運(yùn)動(dòng)、熱傳遞等現(xiàn)象具有良好的優(yōu)勢(shì)。研究人員［4］利用Fluent軟件分析比較了翅翼在三種不同拍動(dòng)模式下的升力和升力系數(shù)。研究人員［5］運(yùn)用Fluent軟件模擬分析了仿蜻蜓翅翼?yè)鋭?dòng)時(shí)的流場(chǎng)特性。研究結(jié)果表明，對(duì)比以前風(fēng)洞實(shí)驗(yàn)等研究方法，利用Fluent軟件使空氣動(dòng)力的計(jì)算變得便捷和準(zhǔn)確，且成本低、周期短，能分析模擬不同狀態(tài)并獲得完整數(shù)據(jù)。

利用Fluent軟件，根據(jù)所設(shè)計(jì)兩自由度撲翼結(jié)構(gòu)的撲動(dòng)運(yùn)動(dòng)規(guī)律，對(duì)兩種不同撲動(dòng)形式進(jìn)行了二維仿真模擬，分析研究?jī)煞N拍動(dòng)形式下流場(chǎng)特性，獲得翅翼在一個(gè)拍動(dòng)周期內(nèi)升阻力系數(shù)的變化情況，并進(jìn)行相互間分析比較。

1 基本控制方程

在對(duì)昆蟲(chóng)和鳥(niǎo)類的長(zhǎng)期觀察中，研究人員發(fā)現(xiàn)翅膀上下?lián)鋭?dòng)的過(guò)程中，不僅有上下平動(dòng)，而且同時(shí)伴隨著翅翼繞某一軸線的規(guī)律轉(zhuǎn)動(dòng)，可以認(rèn)為是低雷諾數(shù)下的湍流流動(dòng)。黏性不可壓縮流體（二維）的基本控制方程由如下兩個(gè)方程組成。

連續(xù)性方程為

式中：p 為壓力，Pa；v為流體的速度矢量，m/s；Re 為雷諾數(shù)，常規(guī)飛行器流場(chǎng)雷諾數(shù)Re一般都大于106，由于撲翼飛行器在低空飛行，尺寸和飛行速度都較小，故流場(chǎng)雷諾數(shù)Re要比常規(guī)飛行器Re小得多，一般都在104～105以下。撲翼飛行器雷諾數(shù)按照Ellington［6］的定義，可利用翅膀的平均弦長(zhǎng)D和平均翼尖速度V進(jìn)行計(jì)算：

式中：ψ為展翅角；λ為展弦比；f為頻率；l為翅長(zhǎng)；μ為運(yùn)動(dòng)黏度。

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

采用雙曲柄雙搖桿機(jī)構(gòu)，故對(duì)曲柄搖桿機(jī)構(gòu)進(jìn)行簡(jiǎn)要分析。如圖1所示，在機(jī)構(gòu)中，構(gòu)件L1為原動(dòng)件，構(gòu)件L3在構(gòu)件L2帶動(dòng)下做往復(fù)運(yùn)動(dòng)，翅翼固連在構(gòu)件L3上，從而實(shí)現(xiàn)翅翼的上下拍打運(yùn)動(dòng)。四桿尺寸參數(shù)如表1所示。

圖1 鉸鏈四桿機(jī)構(gòu)

根據(jù)機(jī)構(gòu)各桿所構(gòu)成的封閉矢量多邊形，可以得出以下矢量關(guān)系式：

求解可得角位移分量形式為：

其中 θ1、θ2、θ3分別為各構(gòu)件的方位角，對(duì)于給定任意一個(gè)曲柄位置θ1，均可以根據(jù)上式求出連桿和搖桿的方位角 θ2、θ3，相應(yīng)的角速度和角加速度可通過(guò)求一階和二階導(dǎo)數(shù)得到。

在微型電機(jī)驅(qū)動(dòng)下，通過(guò)減速機(jī)構(gòu)，構(gòu)件1以14 Hz的頻率勻速轉(zhuǎn)動(dòng)時(shí)，翅翼的上下?lián)鋭?dòng)角度及角速度隨時(shí)間變化規(guī)律如圖2所示，翅翼的最大撲動(dòng)角和最小撲動(dòng)角分別為 124°、50.37°，整個(gè)拍動(dòng)幅度約為 73°。

表1 四桿機(jī)構(gòu)參數(shù)mm

圖2 角度和角速度變化曲線

考慮到研究人員［7］所設(shè)計(jì)的轉(zhuǎn)動(dòng)機(jī)構(gòu)過(guò)于復(fù)雜，故在凸輪機(jī)構(gòu)的基礎(chǔ)上引入了球形鉸鏈，有效地簡(jiǎn)化了撲翼結(jié)構(gòu)，減輕了整機(jī)重量。如圖3所示，由微電機(jī)驅(qū)動(dòng)，經(jīng)過(guò)減速機(jī)構(gòu)帶動(dòng)曲柄轉(zhuǎn)動(dòng)，通過(guò)連桿傳遞實(shí)現(xiàn)翅翼上下往復(fù)運(yùn)動(dòng)，同時(shí)由于翅根處凸輪機(jī)構(gòu)的轉(zhuǎn)向作用，從而使得翅翼從單一沉浮運(yùn)動(dòng)演變?yōu)槌粮〖痈┭龅膬勺杂啥冗\(yùn)動(dòng)模式，圖4為凸輪扭轉(zhuǎn)角度變化規(guī)律，如圖所示，翅翼可以實(shí)現(xiàn)約為40°的扭轉(zhuǎn)。

圖3 三維實(shí)體模型

圖4 凸輪扭轉(zhuǎn)角度變化曲線

3 撲動(dòng)過(guò)程中的升阻力系數(shù)分析

3.1 網(wǎng)格劃分及動(dòng)網(wǎng)格設(shè)置

網(wǎng)格劃分之前，選取合理計(jì)算區(qū)域十分關(guān)鍵，合理的求解區(qū)域?qū)⒅苯佑绊懻麄€(gè)計(jì)算結(jié)果的準(zhǔn)確性以及計(jì)算時(shí)間。采用橢圓面為翅翼的弦向截面，從而建立計(jì)算模型。將翼型置于800 mm×400 mm的矩形區(qū)域內(nèi)，為提高網(wǎng)格適應(yīng)性、計(jì)算精度和節(jié)約計(jì)算時(shí)間，網(wǎng)格采用三角形非結(jié)構(gòu)化網(wǎng)格。利用Gambit生成網(wǎng)格，如圖5所示，共有30 396個(gè)節(jié)點(diǎn)數(shù)個(gè)，網(wǎng)格數(shù)為60 122個(gè)。

圖5 全局網(wǎng)格圖

利用CFD軟件Fluent中的用戶自定義函數(shù)（UDF），使用void DEFINE_CG_MOTION函數(shù)編譯動(dòng)網(wǎng)格UDF文件控制翅翼的上下拍動(dòng)的運(yùn)動(dòng)規(guī)律，從而建立所需的動(dòng)態(tài)網(wǎng)格模型。網(wǎng)格更新主要有三種動(dòng)網(wǎng)格方式：彈簧近似光滑模型、動(dòng)態(tài)分層模型和局部重構(gòu)模型。彈簧近似光滑模型計(jì)算精度高，但在計(jì)算區(qū)域變形較大時(shí)，會(huì)影響計(jì)算精度；動(dòng)態(tài)分層模型要求與運(yùn)動(dòng)邊界相鄰的網(wǎng)格必須為楔形或者六面體網(wǎng)格；局部重構(gòu)模型僅能用于四面體網(wǎng)格和三角形網(wǎng)格。根據(jù)翅翼模型實(shí)際運(yùn)動(dòng)情況綜合考慮，確定采用彈簧近似光滑模型和局部重構(gòu)模型相結(jié)合的方法進(jìn)行動(dòng)網(wǎng)格的生成。

3.2 沉浮運(yùn)動(dòng)的升阻力系數(shù)

在FLUENT中有很多湍流模型可供選擇，在撲動(dòng)過(guò)程中，翅翼的上下往復(fù)運(yùn)動(dòng)使得翅翼附近的空氣流動(dòng)復(fù)雜多變，因此擬采用RNG k-ε模型。針對(duì)撲翼模型在流場(chǎng)的撲動(dòng)飛行，設(shè)定無(wú)窮遠(yuǎn)處自由來(lái)流速度U∞=3 m/s，在控制體上利用simple算法［8］對(duì)其進(jìn)行數(shù)值離散。翅翼在流場(chǎng)中快速拍打，翼背（上表面）和翼面（下表面）產(chǎn)生壓力差，從而產(chǎn)生升阻力。升阻力系數(shù)計(jì)算公式為

式中：ρ為空氣密度；c為翼型的弦長(zhǎng)；U為周期內(nèi)翼型表面來(lái)流平均速度。

在Fluent中，可以直接對(duì)升阻力系數(shù)進(jìn)行監(jiān)測(cè)，圖6為Fluent軟件求解流程圖。

圖6 Fluent求解流程圖

設(shè)定拍動(dòng)頻率為14 Hz，時(shí)間步長(zhǎng)為0.001 s，計(jì)算300個(gè)時(shí)間步長(zhǎng)，迭代了6000次后，升阻力系數(shù)變化規(guī)律曲線如圖7所示。

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

從圖中可以看出在下拍過(guò)程中，升力系數(shù)急劇增加，約在0.25T時(shí)達(dá)到了峰值，此時(shí)在上表面形成的漩渦最大，上下表面的空氣流動(dòng)速度的差異使得上下表面形成了壓力差，上表面壓力小于下表面，便形成了升力。同時(shí)由于攻角存在，后緣渦強(qiáng)度大于前緣渦強(qiáng)度，表現(xiàn)為阻力系數(shù)的增加。上拍過(guò)程中升阻力系數(shù)變?yōu)樨?fù)值。在整個(gè)拍動(dòng)周期，凈升阻力系數(shù)為正值。下拍時(shí)流場(chǎng)速度及壓力分布如圖8～圖9所示。

圖8 下拍時(shí)流場(chǎng)速度分布

圖9 下拍時(shí)流場(chǎng)壓力分布

3.3 沉浮加俯仰運(yùn)動(dòng)的升阻力系數(shù)比較

在凸輪機(jī)構(gòu)的轉(zhuǎn)向作用下，可以實(shí)現(xiàn)翅翼約為40°的扭轉(zhuǎn)。翅翼在下拍過(guò)程中與來(lái)流方向呈10°夾角，到最低位置時(shí)經(jīng)過(guò)凸輪機(jī)構(gòu)的作用，使得翅翼在上拍過(guò)程與來(lái)流方向夾角為50°。可根據(jù)仿真數(shù)據(jù)將扭轉(zhuǎn)變化規(guī)律簡(jiǎn)化為如下方程：

運(yùn)用與上述相同的研究方法，將轉(zhuǎn)動(dòng)的運(yùn)動(dòng)模式加入翅翼的運(yùn)動(dòng)中，可以得到升阻力系數(shù)的變化曲線，如圖10所示。圖中，A為沉浮運(yùn)動(dòng)時(shí)的升阻力系數(shù)，B為沉浮加俯仰運(yùn)動(dòng)時(shí)的升阻力系數(shù)，通過(guò)曲線對(duì)比可知，沉浮加俯仰運(yùn)動(dòng)時(shí)，上拍過(guò)程中，升力系數(shù)明顯升高，同時(shí)阻力系數(shù)也有所降低，整個(gè)拍動(dòng)周期翅翼獲得的凈升力系數(shù)要大于單一的沉浮運(yùn)動(dòng)時(shí)的凈升力系數(shù)，凈阻力系數(shù)要小于單一沉浮運(yùn)動(dòng)時(shí)的凈阻力系數(shù)。上拍時(shí)流場(chǎng)速度與壓力分布如圖11～圖12所示。

圖10 升阻力系數(shù)對(duì)比圖

圖11 上拍時(shí)流場(chǎng)速度分布

圖12 上拍時(shí)壓力流場(chǎng)分布

4 結(jié)論

CFD軟件Fluent可有效模擬模型翼在流場(chǎng)中的運(yùn)動(dòng)情況，所提供的自定義函數(shù)（UDF）窗口把翅翼相對(duì)復(fù)雜的運(yùn)動(dòng)模式轉(zhuǎn)化成可編譯的外部C程序，有效解決了數(shù)值模擬中流場(chǎng)N-S方程所存在的“動(dòng)邊界問(wèn)題”，結(jié)果表明運(yùn)用Fluent軟件模擬分析流場(chǎng)相關(guān)特性是可靠的。

研究結(jié)果表明，翅翼在不同的運(yùn)動(dòng)模式下的升阻力系數(shù)存在著差異，由于凸輪轉(zhuǎn)向機(jī)構(gòu)的作用，兩自由度的機(jī)構(gòu)在撲動(dòng)中獲得了更大的凈升力系數(shù)，同時(shí)也降低了凈阻力系數(shù)，有效增大了升力，減小了阻力。

研究局限在機(jī)構(gòu)有無(wú)扭轉(zhuǎn)情況下的升阻力系數(shù)，且均采用二維簡(jiǎn)化模型，對(duì)于具體翼型產(chǎn)生的影響及更復(fù)雜的撲翼運(yùn)動(dòng)規(guī)律還需更深入的研究。

［1］ Wang J.Vortex shedding and frequency selection in flapping flight［J］.JFluid Meeh，2000，410：323-341.

［2］ 蘭世隆，孫茂.模型昆蟲(chóng)翼非定常運(yùn)動(dòng)時(shí)的空氣動(dòng)力特性［J］.力學(xué)學(xué)報(bào)，2001，33（2）：173-182.

［3］ Sun M，Hossein H.Force and flow structures of an airfoil performing some unsteady motions at small Reynolds number［J］.Aeta Aerodyamica Siniea，2000，18（Suppl）：96-102.

［4］ 朱洪波，謝進(jìn)，陳永.用Fluent分析仿生翅翼運(yùn)動(dòng)產(chǎn)生的升力和升力系數(shù)［J］.機(jī)械設(shè)計(jì)與制造，2011(5)：68-70.

［5］ 張孝松，基于蜻蜓翅翼的仿生微撲翼飛行器機(jī)翼的有限元分析［D］：哈爾濱：哈爾濱工業(yè)大學(xué)，2006.

［6］ E1lington C.The aerodynamics of hovering insect flight∶I.the quasi-steady analysis ［J］.Philosophical Transactions of the Royal Society of London B，1984,305∶1-15.

［7］ 朱寶，王姝歆.兩自由度撲翼機(jī)構(gòu)及其運(yùn)動(dòng)仿真研究［J］.中國(guó)制造業(yè)信息化，2009,38（21）：24-28.

［8］ 周建華.微型拍翅式飛行機(jī)器人翅運(yùn)動(dòng)及控制系統(tǒng)研究［D］.南京：東南大學(xué)，2005.