劉文君 蔡 毓
(廣西大學計算機與電子信息學院,廣西 南寧 530004)
自然界中的鳥類能夠在復雜和高度動態的環境中飛行,具有穩定性、機動性和高效性,對環境干擾有很強的抵抗力,從而可以進行長距離遷徙。然而,由生物飛行靈感啟發研制的各種規模結構、傳感、穩定性控制的仿生撲翼飛行器還沒有自主地表現出這些特性[1]。仿鳥型撲翼飛行器和仿昆蟲型撲翼飛行器的飛行機理有較大差別。仿昆蟲型撲翼飛行器主要靠撲翼撲動產生升力和推力,由于其撲動頻率高,且機翼質量占比小,所以氣動、結構和飛行力學之間的耦合問題影響較小。對于仿鳥型撲翼飛行器而言,其具有低撲動頻率、高機翼質量占比和耦合關系性強等特點。
氣動效率是仿生撲翼飛行器研究中備受關注的領域,是該飛行器在飛行狀態下受到的升力、阻力、推力、力的方向、力的大小以及穩定性等影響的客觀反應,也是檢驗飛行器外形布局及氣動關系設計是否合理的重要依據[2]。目前,關于撲翼的氣動效率,在生物觀測、理論分析、數值模擬、飛行試驗4 個層面都在逐步開展研究。在數值模擬方面,Young J[3]采用可壓縮二維Navier-Stokes 求解器,對斜向振蕩的NACA0012 翼型的流場 進行了數值模擬,用數值粒子追蹤法對翼型的尾跡進行了可視化。Lai JCS[4]研究了Strouhal 數推進效率的峰值,得出僅用Strouhal 數不足以表征撲翼推進效率的結論。Trizila P[5]使用代理建模技術對前緣渦、尾流及射流進行了數值研究。肖天航等人[6]提出一種Delaunary 圖映射網格變形技術和非結構嵌套網格的數值計算方法。
雖然針對撲翼運動的研究已經初步揭示了產生升推力的飛行機制,然而定性的認識還不能完全滿足撲翼飛行器設計層面的需求。一個重要的俯仰運動學案例[7]在中沖程、中下沖程不同的迎角下,分析二維、三維CFD 模型流場非定常瞬態狀態,二者經驗系數的變化曲線基本一致。該案例表明二維模型為微型飛行器設計中的實際應用提供了重要的仿真依據。該文在進行二維數值模擬的同時,對結果進行了積分計算,得到整個機翼平面的氣動力和氣動效率。
在研制仿生撲翼飛行器時,可借助仿生學原理,參考鳥類的翼型構造、形狀[8],來研制機翼結構,如圖1 所示。飛行時產生升力克服其自身重量,產生推力克服飛行時的阻力。
該文在實驗中的飛行器機構簡圖,如圖2(a)所示,并將其機翼撲動基于葉素理論將機翼平面按展向分成若干翼帶,進行三維數學建模,如圖2(b)所示。

圖1 鳥翼的構造和不同階段翅膀的形狀變化
機構簡圖的運動角度位移是分析機構運動性能、優化機構的依據。機翼的運動分析,是根據原動件的運動規律,求解從動件的運動規律。機構運動如圖3 所示。

圖2 撲翼飛行器模型機構簡圖及機翼建模
在機構簡圖中的支鏈上引入閉環結構,利用閉環矢量方程建立連桿機構的約束方程。閉環方程為

式中:R1是曲柄,R2、R3、R4為機械連桿,R5為機構中的輔助桿。R1繞固定鉸鏈點O 轉動,即A 點是以O 為原點、R1為半徑的圓上運動的動點。α1是連桿機構的設定角度,α2是R1繞定點O 轉過的瞬時角,α3是翼桿AB 繞動的瞬時撲動角,α4是平板翼型中的機翼機構的位置角。沉浮運動參數示意圖如圖4 所示。沖程平面角為0。全局坐標系(X,Y,Z)固定在沖程平面的中心,Z 方向垂直于沖程平面,Y 方向垂直于機翼軸,X 方向與沖程平面平行。撲翼運動規律接近于簡諧運動,撲動飛行可以同時產生升力和推力。平板翼型沉浮運動位移被定義為

式中:hm是拍打振幅,f 是拍打頻率,t 是拍打時間,h(t)是t 時間內的位移。
該文的數值模擬基于ANSYS 軟件平臺進行。該平臺集成了N-S 方程流場解算器、動態網格生成技術、空氣動力學算法等多方面模塊,可用于理想和多種設置氣體/流體介質的定常/非定常數值模擬仿真計算。

圖3 撲翼機構運動簡圖
流場計算的基本過程是在空間上用有限體積法,將計算域離散成許多小的體積單元,在每個體積單元上對離散后的控制方程組進行求解。這也是流場計算方法的本質。采用分離式求解器(Segregated solver),逐一、順序地求解各變量。
假定初始壓力和速度,確定離散方程的系數及常數項;求解連續方程、動量方程、能量方程;求解湍流方程及其標量方程;判斷當前時間步上的計算是否收斂。

圖4 沉浮運動懸停參數示意圖
在非定常流場數值模擬的動態邊界運動中,動態網格生成是非常重要的組成部分。動網格計算中,網格的動態變化過程可以用3 種模型進行計算,即彈簧近似光滑模型(spring-based smoothing)、動態分層模型(dynamic layering)和局部重劃模型(local remeshing)。該文采用彈簧近似光滑法和局部重劃法。
彈簧近似法適用于該研究中的非結構化網格,網格拓撲不變,通過近似彈簧的壓縮或拉伸,實現網格和模擬流場計算區域的變化。網格邊界處位移變化后,彈簧系統會經過調整,使已破壞的原有平衡達到新的平衡,生成力的大小由胡克定律計算得到。在控制臺窗口掌握區域范圍、體積統計以及連通性信息。網格檢查沒有出現網格體積為負數的狀態。該次數值模擬實驗需對模型進行不同方向的簡諧運動,多次模擬之后得到相應良好的設定尺寸,根據預先設置的最大、最小網格規格,同時采用網格重劃法在網格更新計算區域進行網格重劃。在數值模擬中,還要注意更新區域的完整性。
數值模擬中的控制方程是具有恒定密度和黏度的非定常Navier-Stokes 方程[9]?;趧討B網格,采用二階迎風格式(second order upwind)有限體積法,考慮量值在物理上的值點分布,即曲線的曲率影響,具有二階精度截差,同時也會考慮氣流的流動方向,所以數值解在物理上是合理的。利用BLU-SGS 隱式算法進行時間推進,求解低速非定常雷諾平均N-S 方程。同一算例應用不同的湍流模型,得到的數值結果會存在很大不同,與此同時,采用湍流模型中試驗最成功的一種模型,Spalart-Allmaras 湍流模型[10],可以合理地處理邊界層的黏性影響區域。
文中分別進行了穩態和瞬時模擬計算,對于穩態問題,壓力修正法采用SIMPLEC 算法,對于瞬態問題,同時為了更好地滿足動量方程、連續性方程,采用PISO 算法。
該算例重點討論外部空氣動力學應用的無量綱系數,對氣動進行配置,計算典型的阻力系數?;谠撐难芯康牡退賳栴},周圍流場中的空氣密度保持不變,為不可壓縮流,并且不會出現馬赫數。氣動系數是關于雷諾數的函數CD=f(M,Re)。如圖5 所示,平板模型垂直于流動氣流,這種機構在任何常規配置中都會產生最大的阻力系數。
雷諾數計算公式為:

阻力系數計算公式為:

其中,S 是平板單位跨度的正面面積。文獻[7]中的阻力系數為2.0,該文數值模擬方法下,迭代收斂后的阻力系數為1.9986。數值模擬和理論方法對該算例進行比較,數據吻合,所以將利用該數值模擬方法進行進一步研究。
對不同參數的二維模型進行非定常流場數值模擬,得到已設置的氣流條件下的升力系數和阻力系數,與此同時,將每個翼帶的二維氣動力,通過積分計算形式,得到整個機翼平面氣動力。后續通過如下數據處理方法做進一步研究。
減縮頻率定義為

雷諾數表示遷移慣性力與黏性力的比值,定義為

經過無量綱化,拍打運動的升力系數和推力系數定義為

式中:CL、CT分別是升力系數和推力系數定義,L、T 分別是升力和推力,是來流密度,S 是平板跨度的正面面積。
依據文獻[11],平均推力系數和平均氣動力輸入系數分別定義為

同時平均氣動力系數定義為

依據文獻[12],升力效率和推進效率分別定義為

對于撲翼飛行,機翼產生顯著升力來抵消自身重力。氣動效率是與給定功率的升力和推力有關的運動效率,能保持在空中的能源供應。該文從上述2 個方面對撲翼氣動效率展開研究。參數空間包括撲打頻率f,風速,減縮頻率k,雷諾數Re。在二維撲翼無粘平板經典理論中,對于低振幅、低減縮頻率的情況,撲翼在周期內的平均推力不依賴初始攻角[13],并且對于三維流動的撲翼,只要運動不發生明顯的流動分離現象,同樣適用該理論[13]。該文數值研究中平均攻角α 設置為零。
低頻振蕩是導致流動失穩的一個主要因素[13],由飛行速度和撲打頻率共同作用。對于撲翼飛行,機翼產生顯著升力來抵消自身重力。在不同撲打頻率f 和風速下,總氣動力會隨撲打頻率的提高而增大。在飛行器起飛階段,提高撲打頻率,增加推力,升力會在達到一定撲打頻率后有增加的趨勢。結合如圖5 所示的氣動效率實驗結果,在撲打頻率為7 Hz、8 Hz 時,具有較高速度,達到最大升力效率。進入巡航階段,降低撲打頻率,達到最大推進效率。
該節實驗計算考慮有限翼展的拍打運動,針對k=0.2~1.5、Re=1.7×104~5.2×104的減縮頻率、低雷諾參數進行大量數值計算,得到定性以及定量的結果。
隨著減縮頻率k 的增加,實驗中雷諾數下的氣動力均增加。升力效率和推進效率分別如圖6 所示。計算出的氣動效率表明,隨著雷諾數Re 的增加,實驗中減縮頻率k 下的推進效率逐步下降至接近平緩,而升力效率的結果變化更陡峭,出現較高的峰值。
從續航時間考慮,應以較低的撲動頻率飛行,較高的撲動頻率對撲翼飛行沒有本質優化,卻會造成機翼架構的損壞,正負升力抵消。

圖5 氣動效率隨撲打頻率和風速變化的比較

圖6 氣動效率隨減縮頻率和雷諾數變化的比較
為了探討仿生撲翼飛行器產生升力和推力的飛行機理,該文在研究考慮撲翼氣動力和氣動效率數值計算方法的基礎上,利用數值模擬仿真,針對不同撲動頻率、風速、減縮頻率、雷諾數、撲動幅度等參數變量,研究揭示了該飛行器撲翼機構數值模擬計算的總氣動力和氣動效率的變化,從而研究分析主要飛行規律。為仿生撲翼飛行器的進一步研究提供了參考依據。