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

梁志宏，李 倫，張玉光

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

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

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

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

1 非定常氣動理論

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

(1)

(2)

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

雷諾數可以表示為：

(3)

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

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

(4)

(5)

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

2 撲翼幾何模型

2.1 三維撲翼氣動分析模型

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

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

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

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

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

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

(6)

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

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

圖3 不同k時撲翼拍動規律曲線

2.2 流場域構建及網格劃分

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

圖4 流場域網格劃分

3 撲翼氣動力計算方法

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

(7)

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

(8)

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

(9)

圖6 撲翼氣動力計算模型

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

(10)

(11)

(12)

(13)

4 仿真結果和分析

4.1 撲翼氣動力特性分析

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

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

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

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

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

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

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

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

4.2 撲翼運動流場分析

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

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

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

5 結論

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