李尚斌，羅 駿，江露生

(中國直升機設計研究所，直升機旋翼動力學重點實驗室，江西 景德鎮(zhèn) 333000)

0 引言

撲翼飛行器屬于仿生范疇，自然界中鳥類和昆蟲等飛行生物均采用撲翼飛行方式。

國內(nèi)外已有研究者對撲翼飛行進行過研究。國外，Jones等[1]針對撲翼尾跡，展開了試驗和計算相關研究；Smith等[2]采用面元法并結合有限元模型對柔性撲翼運動進行了研究；J.C.S.Lai等[3]針對不同的減縮頻率和振蕩位移展開試驗研究，詳細分析了振蕩翼型的尾跡結構及產(chǎn)生推力的原因。Angela等[4]對鴿子起飛和著陸動力學進行了研究。國內(nèi)，宋筆鋒等[5-7]開展了翼型厚度、翼型開孔對撲翼的氣動特性的影響以及微型撲翼的推進特性試驗等研究；昂海松等[8-9]開展了適合于撲翼的動態(tài)網(wǎng)格生成方法及仿生撲翼機構設計等研究。

本文將撲動簡化為俯仰、沉浮和平移三種運動，通過求解NS方程，數(shù)值模擬了二維翼型撲動非定常流動，分析了減縮頻率、平均俯仰角、俯仰振幅和沉浮振幅對翼型撲動氣動特性的影響規(guī)律以及流場渦分布特性。

1 計算模型

選NACA0012翼型為計算模型，將撲動簡化為俯仰、沉浮和前飛三種運動，運動規(guī)律用以下函數(shù)表示：

α=α0+α1×sin(2×π×t/T);

h=H×sin(2×π×t/T+φ);

v=V;

式中，t(s)為時間，T(s)為撲動周期；α0(°)為平均俯仰角，α1(°)為俯仰振幅，α(°)為t時間俯仰角；H(m)為沉浮振幅，h(m)為t時間沉浮距離；φ(rad)為相位角；v(m/s)為前飛速度；k為減縮頻率，ω(rad/s)為撲動角速度,C(m)為翼型弦長。

2 數(shù)值模擬

2.1 計算方法

網(wǎng)格為結構O型網(wǎng)格，第一層網(wǎng)格距物面的距離為1×10-5倍弦長。

以NS方程作為流動控制方程，并采用有限體積法對方程進行空間離散,差分格式為二階迎風格式，湍流模型采用SST模型。NS方程如下:

(1)

其中W是守恒變量，F(xiàn)(W)是無粘通量，G(W)為粘性通量，具體表達如下：

(2)

(3)

(4)

2.2 算例驗證

這里選取Lee等人[10]典型的低雷諾數(shù)試驗狀態(tài)來驗證該文的數(shù)值方法?？s減頻率k=0.1，雷諾數(shù)Re=1.35×105，來流馬赫數(shù)Ma∞=0.0385，平均攻角α0=10°，俯仰幅度α1=15°，無沉浮運動。

圖1-圖3分別為升力系數(shù)、阻力系數(shù)和力矩系數(shù)隨攻角在一個周期內(nèi)的變化圖。從圖中可以看出，計算值和試驗值變化趨勢一致，誤差在合理范圍內(nèi)，計算結果可信。

圖1 升力系數(shù)隨攻角變化圖

圖2 阻力系數(shù)隨攻角變化圖

圖3 力矩系數(shù)隨攻角變化圖

3 結果分析

結果分析中出現(xiàn)的公式如下：

(5)

式中，L為升力，垂直來流方向，向上為正，D為阻力，平行來流方向，負數(shù)表示向前推，M為力矩，力矩點在0.25弦線處；ρ為空氣密度，v為來流風速，S為面積，C為弦長；CL為升力系數(shù)，CD為阻力系數(shù)，CM為力矩系數(shù)。

結果分析中，俯仰和沉浮運動相位差為π/2，0～0.5t/T內(nèi)為上撲階段，0.5～1.0t/T內(nèi)為下?lián)潆A段，0t/T時刻處在最低位置，0.5t/T時刻處在最高位置，0.25t/T時刻俯仰角最大，0.75t/T時刻俯仰角最小。

3.1 減縮頻率影響分析

圖4-圖6為α0=10°，α1=15°，H/C=1狀態(tài)下不同減縮頻率的升力系數(shù)、阻力系數(shù)和力矩系數(shù)隨時間變化圖。k=0.2、0.25、0.30、0.35分別對應Re=82726、66207、55150、47253。

圖4 不同減縮頻率升力系數(shù)隨時間變化圖

圖5 不同減縮頻率阻力系數(shù)隨時間變化圖

從圖中可以看出，對于CL，k越小，平均升力越大；下?lián)鋾rCL比上撲的大；k值大時，上撲產(chǎn)生向下的升力，下?lián)洚a(chǎn)生向上的升力，k值小時，上撲和下?lián)渚a(chǎn)生向上的升力。對于CD，k越大，向前的平均推力越大；k值大時，上撲和下?lián)渚a(chǎn)生向前的推力，k值小時，上撲產(chǎn)生向后的阻力，下?lián)洚a(chǎn)生向前的推力。對于CM，0～0.25t/T和0.75～1.0t/T，k越小，CM越大;0.25～0.5t/T，k越大，CM越大;0.5～0.75t/T，不同k值CM差別不大。

圖6 不同減縮頻率力矩系數(shù)隨時間變化圖

3.2 平均俯仰角影響分析

圖7-圖9分別為Re=66207，α1=15°，H/C=1，k=0.25狀態(tài)下不同平均俯仰角的升力系數(shù)、阻力系數(shù)和力矩系數(shù)隨時間變化圖。

圖7 不同平均俯仰角升力系數(shù)隨時間變化圖

圖8 不同平均俯仰角阻力系數(shù)隨時間變化圖

圖9 不同平均俯仰角力矩系數(shù)隨時間變化圖

從圖中可以看出，對于CL，下?lián)鋾rCL比上撲的大；α0值大時，上撲和下?lián)渚a(chǎn)生向上的升力，α0值小時，上撲產(chǎn)生向下的升力，下?lián)洚a(chǎn)生向上的升力；當α0≤α1時，α0越大，平均CL越大，α0=20°大于α1時，CL反而開始減小。對于CD，α0越小，向前的平均推力越大，但CD的相對波動越大；α0值小時，上撲和下?lián)渚a(chǎn)生向前的推力，α0值大時，上撲產(chǎn)生向后的阻力，下?lián)洚a(chǎn)生向前的推力。對于CM，0～0.25t/T，α0越大，CM越大，0.25～0.65t/T，α0越小，CM越大，0.65t/T后規(guī)律不明顯。

3.3 俯仰振幅影響分析

圖10-圖12分別為Re=66207，α0=10°，k=0.25，H/C=1狀態(tài)下不同平均俯仰振幅的升力系數(shù)、阻力系數(shù)和力矩系數(shù)隨時間變化圖。

從圖中可以看出，對于CL，下?lián)鋾rCL比上撲的大；α1值小時，上撲產(chǎn)生向前的推力，下?lián)洚a(chǎn)生向后的阻力；α1值大時，上撲產(chǎn)生向后的阻力，下?lián)洚a(chǎn)生向前的推力；當α1=α0時，上撲和下?lián)渚墚a(chǎn)生向前的推力，但下?lián)浜箅A段產(chǎn)生向后的阻力。對于CM，類似k對CM的影響規(guī)律。

圖10 不同俯仰振幅升力系數(shù)隨時間變化圖

圖11 不同俯仰振幅阻力系數(shù)隨時間變化圖

圖12 不同俯仰振幅力矩系數(shù)隨時間變化圖

3.4 沉浮振幅影響分析

圖13-圖15分別為Re=66207，α0=10°，α1=15°，k=0.25狀態(tài)下不同沉浮振幅的升力系數(shù)、阻力系數(shù)和力矩系數(shù)隨時間變化圖。

圖13 不同沉浮振幅升力系數(shù)隨時間變化圖

從圖中可以看出，對于CL，H越小，平均升力越大；下?lián)鋾rCL比上撲的大；H值小時，上撲和下?lián)渚a(chǎn)生向上的升力，H值大時，上撲產(chǎn)生向下的升力，下?lián)洚a(chǎn)生向上的升力。對于CD，H越大，向前的平均推力越大；H值大時，上撲和下?lián)渚墚a(chǎn)生向前的推力，但下?lián)浜髸r間段也產(chǎn)生向后的阻力，H值小時，上撲產(chǎn)生向后的阻力，下?lián)淝皶r間段產(chǎn)生向前的推力，下?lián)浜髸r間段產(chǎn)生向后的阻力。對于CM，H大時，周期范圍內(nèi)CM值波動范圍大；H小時，上撲CM值變化相對平緩，下?lián)渥兓蟆?/p>

圖14 不同沉浮振幅阻力系數(shù)隨時間變化圖

圖15 不同沉浮振幅力矩系數(shù)隨時間變化圖

3.5 流場渦分布

圖16為Re=66207，α0=10°，α1=15°，H/C=1，k=0.25狀態(tài)下不同時刻渦線分布圖。t/T=0位于最低點，t/T=0.5位于最高點。從圖中可以看出，最低點時翼型前緣和后緣出現(xiàn)明顯的渦，上撲過程中，前緣和后緣脫出的渦逐漸耗散，但前緣拖出的渦耗散得更快；下?lián)溥^程中，前緣拖出新的渦，并逐漸向后移動，渦強也隨之增強。對比(a)和(h)圖可以看出，(a)圖中后緣渦是由前緣脫出的并移動到后緣，(h)圖中前緣脫出的渦強度明顯比后緣脫出的大。

圖16 不同時刻渦線變化圖

4 結論

通過數(shù)值模擬分析，得出了以下結論：

1)研究范圍內(nèi)，增大平均俯仰角、俯仰振幅或減小減縮頻率和沉浮振幅能有效提高升力；

2)研究范圍內(nèi)，減小平均俯仰角或增大縮減頻率、俯仰振幅、沉浮振幅能有效增大向前推力；

3)研究范圍內(nèi)，升力主要由下?lián)潆A段產(chǎn)生，α1≥α0時，可獲得升力及推力的綜合效果；

4)下?lián)溥^程中，前緣脫出渦并向后移動，強度逐漸變強，上撲過程中，渦逐漸耗散并繼續(xù)向后移動。