仿鳥型撲翼飛行器非對稱運動氣動特性數(shù)值分析

梁志宏，李 倫，張玉光

(1.沈陽理工大學 裝備工程學院，沈陽 110159; 2.濰坊學院 信息與控制工程學院，山東 濰坊 261061;3.遼沈工業(yè)集團有限公司(國營724廠)，沈陽 110045)

自然界中的鳥類經(jīng)過億萬年的飛行進化，形成了出色的飛行能力，飛行效率、機動性都遠遠超過了目前所有的人造飛行器。鳥類飛行過程中翅膀的運動由上下?lián)鋭?、弦向扭轉(zhuǎn)、展向折彎和前后揮擺等幾種運動復合而成[1-3]。鳥類飛行時翅膀通過周期性的復合運動產(chǎn)生飛行所需的上升升力和前向推力，實現(xiàn)爬升、俯沖、盤旋、懸停等動作，從而具備高超的機動性能。而撲翼飛行器根據(jù)鳥類翅膀撲動而設(shè)計的一種仿生飛行器，具有體積小、質(zhì)量輕、能耗低、噪音低，機動性能優(yōu)越、隱蔽性強等特點。因此，小型的撲翼飛行器是未來無人飛行器發(fā)展的一個重要方向[4,6]。

目前，國內(nèi)外對撲翼飛行器的研究主要集中在2個方面：一方面是仿生撲翼飛行器的結(jié)構(gòu)設(shè)計與運動學分析，通過觀察鳥類翅膀的運動規(guī)律，建立撲翼運動學與動力學模型，設(shè)計高效驅(qū)動機構(gòu)，結(jié)合新型柔性致動材料，多種形式的仿生撲翼飛行器應(yīng)運而生[5]。如加利福尼亞工學院MicroBat、Aero Vironment公司研制的“蜂鳥”、哈佛大學研制的RoboBee X-Wing、荷蘭Delft大學研制的BatBot(B2)的仿蝙蝠撲翼飛行器等。另一方面對撲翼飛行機翼非定常氣動力特性開展研究，通過理論分析、數(shù)值模擬和實驗測試相結(jié)合的方式展開。Tuncer等[16]對撲翼拍動過程中的升阻力變化進行了分析。Delaurier等[7-12,14]采用改進的條帶理論建立了對稱撲翼在飛行過程中翅翼運動的空氣動力學模型。國內(nèi)的研究人員對撲翼的非定常理論也展開了深入研究，昂海松等采用非定常渦格方法對撲翼復合運動的氣動力特性展開研究。孫茂[13]通過Navier-Stokes equations(N-S)方程與渦動力學分析昆蟲翅翼的非定常運動氣動力特性。

目前對撲翼非定常運動的空氣動力特性的研究有很多，且形成了一定的理論基礎(chǔ)。但在氣動分析時，機翼大多為上下對稱撲動。而通過觀察鳥類翅膀的運動規(guī)律，飛行過程中翅膀下拍和上拍的時間占比不同，下拍時間較短速度較快，而上撲占整個撲動周期的時間較長速度較慢。本文針對撲翼非對稱拍動的運動特定，建立撲翼非對稱運動控制方程，構(gòu)建三維撲翼氣動分析模型，研究采用動態(tài)混合網(wǎng)格技術(shù)，并結(jié)合布西涅斯克假設(shè)，求解N-S方程，分析撲翼在運動非對稱時對撲翼氣動力特性、氣動效率以及機翼附近流場變化的影響。

1 非定常氣動理論

撲翼飛行速度較低，飛行過程中空氣的粘性力體現(xiàn)較為明顯，屬于低雷諾數(shù)范疇。撲翼非對稱拍動產(chǎn)生飛行所需的升力和前進推力，撲翼拍動過程中流場時刻發(fā)生變化，此時流場為非定常流動。本文研究撲翼非定常氣動力特性時采用有限體積法離散三維不可壓縮流動Navier-Stokes(N-S)方程，此時流場控制方程的數(shù)學表達式可用連續(xù)性方程與動量方程描述[15]：

(1)

(2)

方程(1)為連續(xù)性方程，(2)為動量方程，u、v、w分別為速度在x、y、z方向上的分量；p為流場壓力；ρ為流體密度；Re為雷諾數(shù)。

雷諾數(shù)可以表示為：

(3)

其中：U為一個運動周期撲翼拍動的平均速度；c為弦長；μt為流體運動黏度。

為使得流體控制方程組封閉，引入布西涅斯克假設(shè)(Boussinesq)：

(4)

(5)

其中：κ表示湍流動能；μt表示流體運動黏度系數(shù)；fυ為黏滯阻尼系數(shù)。

2 撲翼幾何模型

2.1 三維撲翼氣動分析模型

氣動分析時建立的三維撲翼氣動模型如圖1所示。整個模型分為3個部分：左右機翼和機身部分，分析時兩側(cè)機翼繞機身按照預定運動規(guī)律上下拍動。以機身質(zhì)心為坐標原點建立坐標系O-XYZ，X軸與機身軸線重合，由頭部指向尾部為正方向，Z軸垂直于X軸指向左側(cè)機翼，Y軸根據(jù)右手坐標系確定，U∞表示來流風速的方向。機翼拍動過程，機翼位于OXZ面上方，拍動幅值角度α為正，反之為負。機身與水平方向呈8°的攻角。機翼的翼型選用NASA翼型庫中的BE6356，如圖2(b)所示。分析時，由于兩側(cè)機翼的運動規(guī)律相同，為節(jié)約計算資源，數(shù)值計算時只分析單側(cè)機翼拍動時的氣動力特性。

計算采用的撲翼氣動模型尺寸參數(shù)如圖2(a)和表1所示。

圖1 三維撲翼氣動分析模型

圖2 BE6356翼型和撲翼氣動模型幾何參數(shù)示意圖

參數(shù)數(shù)值翼展R/mm1000翼面面積Swing/m20.2機身長度L/mm500弦長b/mm230機身直徑d/mm10

機翼通過拍動提供飛行所需的動力，拍動過程可以分為上拍和下拍過程中。撲翼拍動過程的數(shù)學形式如式(6)所示

(6)

式中：ym為拍動幅值；β0為拍動初始相位；k為拍動時間周期占比；T為拍動周期。

分析時通過調(diào)整k的值，調(diào)整機翼上拍和下拍在一個拍動周期內(nèi)的時間占比，從而實現(xiàn)撲翼快撲慢回動作，達到提供撲翼飛行過程中飛行升力。

圖3 不同k時撲翼拍動規(guī)律曲線

2.2 流場域構(gòu)建及網(wǎng)格劃分

為了避免邊界回流對撲翼拍動時氣動力的影響，構(gòu)建分析流場域時采用翼展長度的15～20倍構(gòu)建，采用四面體網(wǎng)格對整個流場域進行網(wǎng)格劃分，構(gòu)建的流場域及網(wǎng)格劃分如圖4所示，總的網(wǎng)格數(shù)量為 1 194 220，最小網(wǎng)格尺度0.04 mm。仿真分析時，為更好地獲取機翼的氣動力及機翼附近的流場變化情況，因此在處理機翼附近處的網(wǎng)格時，對機翼前后緣和機身頭部、尾部的網(wǎng)格進行了加密。對流場域網(wǎng)格及撲翼模型網(wǎng)格進行優(yōu)化，得到較高的網(wǎng)格質(zhì)量，以獲得更好地計算精度。選用網(wǎng)格的長寬比(Aspect Ratio)作為網(wǎng)格質(zhì)量的評價標準對構(gòu)建的流場域進行優(yōu)化，優(yōu)化后整體的網(wǎng)格質(zhì)量能夠達到0.53以上，如圖5所示。而對于采用動網(wǎng)格方法進行仿真時，網(wǎng)格質(zhì)量達到0.2以上就可以使用，越高的網(wǎng)格質(zhì)量數(shù)值計算時收斂精度越高，收斂速度越快。

圖4 流場域網(wǎng)格劃分

3 撲翼氣動力計算方法

氣動升力的計算，將翅翼分為若干微元，如圖6所示。假設(shè)撲翼飛行過程中，氣動力方向一直垂直于機翼翼面。撲動過程中，微元的瞬時速度為ωr,該微元產(chǎn)生的瞬時氣動升力為：

(7)

Γ=1/2ρair(ωr)2

(8)

其中：Γ為該微元在流體內(nèi)動壓；ρair為空氣密度;CL(α)為微元攻角下的氣動力系數(shù)；c(r)為撲翼展向位置處的弦長。通過對撲翼該微元進行展向積分，即可計算得到整個撲翼的瞬時氣動力為：

(9)

圖6 撲翼氣動力計算模型

同理可以得到撲翼撲動過程中的氣動阻力FD、推力FT、功率PW、推進效率η等參量為：

(10)

(11)

(12)

(13)

4 仿真結(jié)果和分析

4.1 撲翼氣動力特性分析

數(shù)值仿真分析時，為提高求解計算精度，采用采用雙時間步隱式迭代方法。計算采用壓力速度耦合求解，由于撲翼運動時會由于翅翼撲動會產(chǎn)生推力作用的湍流，因此本文數(shù)值計算過程采用湍流模型中的切應(yīng)力傳輸(SSTk-ω)模型。

撲翼的氣動力特性受多種因素的影響，拍動幅值、拍動頻率、來流速度等參數(shù)均能夠使得撲翼的氣動性能產(chǎn)生改變。而本文在分析撲翼的氣動力特性時主要考慮撲翼的非對稱拍動對氣動力產(chǎn)生的影響。因此，在數(shù)值仿真時，拍動幅值、拍動頻率、空氣來流速度、雷諾數(shù)均保持不變，分別設(shè)置為：ym=π/4，f=1/T=5 Hz，v=5 m/s，Re=2 500。

數(shù)值仿真殘差曲線收斂后，不同拍動時間占比k=0.2、0.3、0.4、0.5時的瞬時升力、瞬時阻力與瞬時力矩隨時間變化的曲線如圖7所示。

圖7 不同k值時撲翼瞬時氣動力特性

分析不同時間占比時撲翼的瞬時氣動力變化曲線得知，隨著機翼隨著下?lián)渌俣鹊脑隹?，機翼產(chǎn)生的升力越大。當k分別為0.2、0.3、0.4、0.5時，撲翼上下拍動產(chǎn)生的最大升力分別為28.15 N、14.67 N、11.64 N和10.07 N。下?lián)溥\動占整個撲動周期時間的20%，上撲占80%時，撲翼產(chǎn)生的升力較對稱拍動時產(chǎn)生的升力增加了18.08 N。由此可以看出，通過提高撲翼的下拍速度，可以有效提高撲翼飛行過程中產(chǎn)生的瞬時升力。而分析阻力變化曲線可知，撲翼拍動過程中的阻力隨著下?lián)渌俣鹊募涌熘饾u增大。當k=0.2、0.3、0.4、0.5時，撲翼拍動過程產(chǎn)生的瞬時阻力最大值分別達到了2.92 N、1.36 N、0.93 N和0.34 N。下拍速度的增大，使得撲翼拍動時產(chǎn)生的阻力也明顯增大。k=0.2時(下拍速度最大)的最大阻力較k=0.5(對稱拍動)時增大了2.58 N。因此，在設(shè)計撲翼上下拍動時間占比時，要充分考慮拍動時間占比所引起的升力與阻力的變化情況。

觀察分析圖6(c)中撲翼瞬時力矩變化曲線可知，當k=0.3和0.4時，撲翼產(chǎn)生的瞬時力矩基本相同，k=0.5時，撲翼產(chǎn)生的力矩與k=0.3和0.4時的力矩相差不大。而當機翼以較快的下拍速度時(k=0.2)，機翼產(chǎn)生的力矩明顯增大。由于機翼拍動時產(chǎn)生的力矩為繞機身軸線的力矩。過大的力矩會對撲翼飛行過程中的姿態(tài)穩(wěn)定產(chǎn)生影響。因此，圖6所示為撲翼拍動時產(chǎn)生是瞬時升力變化曲線，其代表的是撲翼在拍動過程中某時刻機翼瞬時氣動力特性，為更加全面地分析撲翼在整個拍動運動周期的升阻變化情況，圖8所示給出了撲翼不同k值時的平均升阻力及氣動效率的變化曲線。

圖8 不同k時撲翼平均氣動力及氣動效率變化曲線

由圖8所示不同k時撲翼的平均升力變化曲線可以得知，撲翼拍動產(chǎn)生的平均升力與平均阻力隨著k值的增大而逐漸減小。平均升力與平均阻力的變化與瞬時升力峰值的變化趨勢相一致。而平均力矩隨著k值的增大先減小后增大，當k=0.4時，撲翼產(chǎn)生的力矩最小。而分析不同k值時撲翼的氣動效率變化曲線可知，撲翼的氣動效率隨著k值的增大而逐漸增大。當k=0.5時，撲翼的氣動效率達到最大值，為29.35，這是由于當上下對稱撲動時，撲翼拍動過程中的阻力減小，使得撲翼在整個拍動過程中的氣動效率增大。k=0.2時撲翼產(chǎn)生的升力最大，但由于下?lián)渌俣冗^快，上拍速度較慢，導致上拍過程中的阻力增大很多，從而使得整體的氣動效率下降。因此在設(shè)計撲翼飛行器是，要確定合理的上拍與下拍的時間占比，從而獲得較大升力的同時保證具備較高的撲翼氣動效率。

4.2 撲翼運動流場分析

根據(jù)上述分析，通過調(diào)整上下拍動的周期時間占比k值能夠改變撲翼拍動過程中撲翼產(chǎn)生的升力及阻力特性，而在撲翼的非對拍動時，機翼附近的流場也時刻發(fā)生變化。如圖9所示為不同時刻撲翼上拍與下拍過程中上下翼面的壓力變化云圖。

分析圖9(a)所示的撲翼下拍過程翼面的壓力變化。隨著撲翼向下拍動，撲翼下方的壓力逐漸增大，當撲翼拍動到中間位置時刻時，此時撲翼下拍的速度達到最大值，上下翼面的壓差達到最大，氣動力出現(xiàn)峰值。隨著撲翼繼續(xù)向下拍動，拍動速度逐漸減小，上下翼面的壓差開始減小，當拍動到最低位置時刻時，上下翼面的壓差最小。從下拍過程中上下翼面的壓力變化中可以看出，下拍開始，齒根處開始出現(xiàn)前緣渦，隨著撲翼不斷向下拍動，前緣渦不斷向翼尖處發(fā)展并不斷變大，當拍動到最低位置時，渦旋從翼尖處脫落，從而產(chǎn)生向前的推力，從而推動撲翼向前飛行。圖9(b)所示為撲翼上拍過程的壓力變化與上拍過程中翼面的變化相反。上拍過程中上翼面的壓力較大，下翼面出現(xiàn)負壓力。上拍階段產(chǎn)生的為負升力，上拍時間過長會導致上拍過程中的阻力增加較多。因此在調(diào)整撲翼上下拍動的時間占比時，要平衡快撲慢回動作時的升阻力變化，從而達到較高的氣動效率。

圖9 一個拍動周期撲翼上下翼面壓力變化過程

5 結(jié)論

本文根據(jù)鳥類翅膀的快撲慢回的運動規(guī)律，建立不同時間周期占比的撲翼拍動運動模型，采用動態(tài)混合網(wǎng)格技術(shù)，數(shù)值求解N-S方程，同時結(jié)合布西涅斯克假設(shè)，分析不同拍動時間周期占比時撲翼的氣動特性的變化。通過數(shù)值仿真結(jié)果表明：撲翼下拍速度越快，撲翼拍動產(chǎn)生的升力越大，然然而下拍速度過快會使得拍動過程中產(chǎn)生的阻力也會隨著增大，導致?lián)湟碚w的氣動效率下降。因此，在設(shè)計撲翼上拍與下拍的時間周期占比時，要充分考慮拍動周期占比對撲翼升阻力的影響，設(shè)計合理的時間周期占比，保證撲翼具備較大升力的同時，使得產(chǎn)生的阻力較小，從而使得撲翼飛行過程中具備較佳的氣動效率。本文的數(shù)值仿真結(jié)果，能夠為仿鳥型撲翼飛行器驅(qū)動機構(gòu)設(shè)計提供理論參考。