翼型撲動氣動特性數值模擬研究

李尚斌，羅 駿，江露生

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

0 引言

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

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

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

1 計算模型

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

α=α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 數值模擬

2.1 計算方法

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

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

(1)

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

(2)

(3)

(4)

2.2 算例驗證

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

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

圖1 升力系數隨攻角變化圖

圖2 阻力系數隨攻角變化圖

圖3 力矩系數隨攻角變化圖

3 結果分析

結果分析中出現的公式如下：

(5)

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

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

3.1 減縮頻率影響分析

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

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

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

從圖中可以看出，對于CL，k越小，平均升力越大；下撲時CL比上撲的大；k值大時，上撲產生向下的升力，下撲產生向上的升力，k值小時，上撲和下撲均產生向上的升力。對于CD，k越大，向前的平均推力越大；k值大時，上撲和下撲均產生向前的推力，k值小時，上撲產生向后的阻力，下撲產生向前的推力。對于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 不同減縮頻率力矩系數隨時間變化圖

3.2 平均俯仰角影響分析

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

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

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

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

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

3.3 俯仰振幅影響分析

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

從圖中可以看出，對于CL，下撲時CL比上撲的大；α1值小時，上撲產生向前的推力，下撲產生向后的阻力；α1值大時，上撲產生向后的阻力，下撲產生向前的推力；當α1=α0時，上撲和下撲均能產生向前的推力，但下撲后階段產生向后的阻力。對于CM，類似k對CM的影響規律。

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

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

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

3.4 沉浮振幅影響分析

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

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

從圖中可以看出，對于CL，H越小，平均升力越大；下撲時CL比上撲的大；H值小時，上撲和下撲均產生向上的升力，H值大時，上撲產生向下的升力，下撲產生向上的升力。對于CD，H越大，向前的平均推力越大；H值大時，上撲和下撲均能產生向前的推力，但下撲后時間段也產生向后的阻力，H值小時，上撲產生向后的阻力，下撲前時間段產生向前的推力，下撲后時間段產生向后的阻力。對于CM，H大時，周期范圍內CM值波動范圍大；H小時，上撲CM值變化相對平緩，下撲變化更大。

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

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

3.5 流場渦分布

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

圖16 不同時刻渦線變化圖

4 結論

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

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

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

3)研究范圍內，升力主要由下撲階段產生，α1≥α0時，可獲得升力及推力的綜合效果；

4)下撲過程中，前緣脫出渦并向后移動，強度逐漸變強，上撲過程中，渦逐漸耗散并繼續向后移動。