迎角對撲翼運動氣動性能影響的高性能數(shù)值研究

祝銜琦，孫祥程，王 嫻，*

(1. 西安交通大學 航天航空學院，機械結(jié)構(gòu)強度與振動國家重點實驗室，西安 710049；2. 陜西省先進飛行器服役環(huán)境與控制重點實驗室，西安 710049)

0 引 言

經(jīng)過上億年的進化，鳥類、昆蟲、蝙蝠等生物均通過撲翼運動實現(xiàn)飛行。國內(nèi)外學者受到撲翼生物的啟發(fā)，研制了撲翼飛行器。撲翼飛行器主要在懸停、升力、能耗、控制等性能方面具有很大的優(yōu)勢。然而，目前撲翼飛行器的發(fā)展仍面臨無法提供更高升力，飛行效率低，尺寸無法達到撲翼昆蟲量級，無法做出更精確的飛行控制等問題[1]。因此，為了突破撲翼飛行器高升力、高效率、小尺寸等技術(shù)瓶頸，亟需開展更為細致的撲翼運動動力學研究，為設(shè)計更高效、高機動性的撲翼飛行器提供可靠的理論基礎(chǔ)。撲翼飛行器在飛行過程中，相對來流方向與翼弦間的角度為迎角。迎角是影響撲翼飛行器飛行時升阻力系數(shù)的重要因素。當迎角超過一定范圍，飛行器的飛行性能將會下降。因此，為了獲得撲翼飛行器飛行時最佳的飛行姿態(tài)，開展不同迎角對撲翼運動的影響機理研究尤為重要。

目前關(guān)于迎角對撲翼運動的影響機理研究主要以實驗為主：2007 年，寧琦等[2]對撲翼飛行器進行了風洞試驗，研究了迎角對其氣動特性的影響；2014 年，鮑鋒等[3]基于色流實驗與粒子圖像測速（particle image velocimetry, PIV）技術(shù)，開展了迎角對單自由度撲翼運動升力的影響機理研究。此外，對撲翼運動其他運動參數(shù)的研究中，1989 年，Ennos[4]等基于高速攝影技術(shù)，對撲翼昆蟲飛行時的運動學規(guī)律進行了探索；1999 年，Dickinson 等[5]通過模型實驗的方法確定了昆蟲飛行時的雷諾數(shù)范圍為40 與400 之間。然而，實驗方法普遍存在代價大、成本高等問題，且實驗?zāi)軌蛱峁┑臄?shù)據(jù)有限，難以精細地獲得能夠直觀、全面反映其空氣動力學機制的流場細節(jié)特征。因此，高效、準確地對撲翼運動進行數(shù)值模擬研究已成為近些年的熱點。國內(nèi)外學者主要通過求解N-S 方程的方法對撲翼運動開展數(shù)值模擬研究。1998 年，Liu等[6]通過求解三維非定常N-S 方程的方法模擬并分析了昆蟲飛行時的非定常流場結(jié)構(gòu)；2000～2004 年，Wang 等利用二維計算[7-9]、三維實驗[10]的方法研究了蜻蜓前飛時流場的變化與前后翼之間的影響，解釋了蜻蜓飛行高升力機制；2002～2015 年，孫茂等[11-14]通過求解N-S 方程的方法，對撲翼昆蟲運動時的穩(wěn)定性高、能耗低等特點進行了研究，并對氣動力機理進行了歸納；2008 年，Young 等[15]通過求解三維不可壓N-S 方程的方法模擬了昆蟲撲翼運動，研究了振幅、頻率等參數(shù)對撲翼運動性能的影響。

撲翼的運動方式復(fù)雜、外形復(fù)雜，可靠的模擬結(jié)果依賴于復(fù)雜邊界條件易實施的數(shù)值方法與高分辨率網(wǎng)格。近年來，高效準確的格子Boltzmann 方法（lattice Boltzmann method, LBM）[16]在計算流體力學（computational fluid dynamics, CFD）中得到普遍應(yīng)用，由于其算法簡單、邊界條件易實施、質(zhì)量動量嚴格守恒、適于并行等優(yōu)點，被應(yīng)用于各個領(lǐng)域，在撲翼運動的數(shù)值模擬中也得到快速發(fā)展。2014 年，Keisuke 等[17]采用浸沒邊界格子Boltzmann 方法針對類蜻蜓撲翼模型前后翼相位滯后角開展了數(shù)值模擬研究，結(jié)果表明，撲翼運動的移動方向取決于相位滯后角。2020年，彭連松等[18]采用基于LBM 的商業(yè)軟件Xflow，對串列雙翅撲翼和單對翅撲翼進行數(shù)值模擬，結(jié)果表明，雙翅懸停效率高于單翅拍動的懸停效率。此外，對撲翼運動進行數(shù)值模擬時，對運動復(fù)雜邊界的精確描述是十分重要的。動邊界識別方法總體上分為拉格朗日法和歐拉法，基于拉格朗日法的動邊界識別法，其重構(gòu)網(wǎng)格計算時間長，計算效率低。而Level-Set 方法[19]基于歐拉法，即使在規(guī)則分布的歐拉網(wǎng)格中也能精確捕捉到邊界點的位置。同時，由于網(wǎng)格的規(guī)則性，更適合并行計算，使得重構(gòu)網(wǎng)格的時間大大縮短，若在較高的網(wǎng)格分辨率下，可精準描述復(fù)雜物體邊界運動。如2005 年，李會雄等[20]基于Level-Set方法捕捉了氣-液兩相流運動相界面; 2006 年，陳凡紅等[21]基于Level-Set 方法追蹤了各種外形的物體在旋轉(zhuǎn)流場中的變形還原效果，研究表明，Level-Set 方法可以準確追蹤運動界面，且容易實現(xiàn)，具有較大通用性。另一方面，在硬件上，近些年來，圖形處理器（graphics processing unit，GPU）發(fā)展速度遠超中央處理器（central processing unit, CPU），GPU 開始逐漸應(yīng)用于非顯示的通用圖形處理器（general purpose graphics processing unit, GPGPU），出現(xiàn)了較多基于CPU/GPU 異構(gòu)系統(tǒng)的CFD 高性能計算在各類計算流體力學方法（如N-S、LBM 等）上的應(yīng)用研究，使得CFD 計算效率更高。其中GPU 的單指令多線程模式與LBM 算法完美匹配，兩者結(jié)合，較之CPU 可達百倍以上的加速比[22]。

基于以上種種，本文針對撲翼運動問題，采用LB 方法求解流場，Level-Set 方法追蹤復(fù)雜運動邊界，并結(jié)合GPU 并行加速技術(shù)，開發(fā)了撲翼運動三維高效數(shù)值模擬程序，對撲翼在不同迎角下，上下?lián)鋭訒r其氣動性能開展數(shù)值研究。通過精細捕捉其流體力學特征量—渦系結(jié)構(gòu)、壓力場等，分析其升阻力產(chǎn)生的機制，獲得其最佳運動迎角。同時，本文也為撲翼運動的研究提供了高效、可靠的數(shù)值方法。

1 計算模型與數(shù)值方法

1.1 計算模型

本文根據(jù)Liang[12]選用的撲翼模型進行建模，其幾何模型如圖1 所示。

圖1 撲翼幾何模型Fig. 1 Prototype of flapping-wing

翼面的展長L=4.90 cm ，平均弦長D=1.24 cm，厚度s=0.14 cm， θ為迎角。兩個撲翼分別繞各自的翼根作上下?lián)鋭舆\動，撲動方程為：

其中，α為 撲動角，αmax為 撲動角 振 幅，f為撲動頻率，φ為撲動初始相位角。雙翼同步撲動，撲動角振幅αmax=30°， 撲動頻率f=32 Hz ，撲動初始相位 φ=0°。撲動角速度用 ω表示。

雷諾數(shù)Re=UrefD/υ =100， 無量綱D=18，其中Uref=4αmaxfr，Uref為撲翼運動參考速度，r為撲翼的旋轉(zhuǎn)半徑。本文分別計算了 θ=0°～60°時不同迎角的流場情況。流體的計算區(qū)域為Lx×Ly×Lz，每個方向取1 6D，即 288×288×288，總計2 389 萬網(wǎng)格。撲翼的升力系數(shù)及阻力系數(shù)定義為：

其中，F(xiàn)z為 撲翼受到的沿z方向升力，F(xiàn)x為撲翼受到的沿x方向阻力，A為撲翼的面積。

1.2 數(shù)值方法

1.2.1 格子Boltzmann 方法

本文采用LBM-D3Q19 模型[16]（圖2）求解流場，其離散方程如下：

圖2 D3Q19 模型Fig. 2 D3Q19 model

對于流場，宏觀邊界條件為：進口u=U∞,v=w=0；出口以及其余邊界均采用充分發(fā)展邊界條件。分布函數(shù)f邊界條件為：所有邊界均采用非平衡態(tài)外推格式[23]，如下：

1.2.2 LBM 移動固壁邊界條件

1.2.2.1 曲面邊界非平衡態(tài)外推格式

對于移動的固體壁面，本文采用非平衡態(tài)外推格式[23]確定固體邊界。如圖3 所示，非平衡態(tài)外推格式將固體節(jié)點rw的分布函數(shù)fi(rw,t)分解為平衡態(tài)與非平衡態(tài)兩部分。

圖3 曲線邊界非平衡態(tài)外推格式Fig. 3 Schematic of extrapolation method for curved boundary conditions in lattice Boltzmann method

非平衡態(tài)部分中，q表示與壁面相交一個網(wǎng)格間距中處于流體區(qū)域的比例：

1.2.2.2 新生流體節(jié)點的重新填充

物體在流場中運動時，當前時刻是固體中的網(wǎng)格節(jié)點下一時刻可能會成為流體節(jié)點，此時需要重新填充新生流體節(jié)點[24]。如圖4 所示，藍色虛線為t時刻固體壁面，t+δt時刻，固體壁面運動至藍色實線處，p點由固體節(jié)點變?yōu)榱黧w節(jié)點。此時，新生流體節(jié)點p需要進行分布函數(shù)的計算。

圖4 新生流體節(jié)點的重新填充Fig. 4 Schematic of initialization of new fluid nodes

公式如下所示:

1.2.3 Level-Set 動邊界識別方法

Level-Set 算法是一種用隱函數(shù)描述邊界運動的方法[19]，利用等值面函數(shù) φ(x) ， 讓 φ以一定速度運動，零等值面即為固體邊界。在歐拉坐標系下， φ(x)可以用符號距離函數(shù)表示。符號距離函數(shù)表示當前節(jié)點與物體邊界之間的最短距離，它在邊界附近充分光滑，可以用來表示Level-Set 函數(shù)。符號距離函數(shù)d(x)的定義如下：

針對流/固動邊界的實時更新問題，采用Level-Set 動邊界識別方法對標準格子立方體網(wǎng)格的邊界進行實時更新。如圖5 所示，圖5（a）為物體運動前狀態(tài)，每個節(jié)點存儲 φ(x)。 將 φ(x)跟隨物體作平移、旋轉(zhuǎn)運動后，如圖5（b）所示，并將運動后的 φ(x)插值至歐拉背景網(wǎng)格中，更新每個節(jié)點的 φ(x)，物體隨即發(fā)生運動，從而識別出動邊界。

圖5 Level-Set 動邊界識別方法Fig. 5 Schematic of identification method of moving boundary of Level-Set

1.2.4 GPU 程序移植

GPU 并行模式是將數(shù)據(jù)分配于多條線程（thread），多條線程組成一個塊（block），GPU 中多個運算單元同時計算每個block 中的threads，而每條thread 中的數(shù)據(jù)則串行計算。

本文自主搭建了基于LBM-GPU 與Level-Set 動邊界識別方法的數(shù)值計算平臺，計算流程圖如圖6 所示。在流場計算的每一個迭代中，經(jīng)過分布函數(shù)的碰撞與遷移之后，對下一迭代步的邊界位置采用Level-Set 方法進行識別，并映射至LBM 網(wǎng)格中，對網(wǎng)格節(jié)點屬性，即固體內(nèi)部、固體壁面及流體，進行更新，并對新生流體節(jié)點的分布函數(shù)進行填充，完成一次循環(huán)。

圖6 計算流程圖Fig. 6 Flow chart of numerical method

此外，本計算采用一顆NVIDIA Tesla K20M GPU加速，以2 389 萬個網(wǎng)格計算工況為例，計算7 200 個LBM 步長，總耗時4.20 h；相同情況下使用一顆Intel Core i7-3770 CPU 計算共耗時298.2 h，加速比為71。

1.3 程序驗證

1.3.1 正確性驗證

本文利用圓柱沉降問題進行程序驗證，并與Li[25]的模擬結(jié)果進行對比。文獻采用應(yīng)力積分法求得物體受力，從而獲得物體運動的軌跡，實現(xiàn)追蹤動邊界的效果。如圖7 所示，圓柱在一個充滿流體的管道中自由下落，由于自身的重力以及流體的作用，圓柱會邊旋轉(zhuǎn)邊平移地下沉。管道寬度L=0.4 cm，圓柱直徑d=0.1 cm， 圓柱密度為 ρc=1.003 g/cm3。流體密度 ρ=1 g/cm3，運動黏度 υ =0.01 cm2/s，重力加速度G=980 cm/s2，圓柱在靠近左側(cè)壁面x=0.076 cm位置開始釋放。

圖7 圓柱沉降示意圖Fig. 7 Schematic of sedimentation of a circular cylinder

圖8 為圓柱下沉過程中圓柱下落軌跡、水平方向速度、豎直方向速度和角速度隨時間變化的曲線圖。可以看出，本文結(jié)果與文獻對比一致。此外，從圖8（a）可以看出，圓柱運動最終達到了以恒定速度up豎直下落的狀態(tài)。穩(wěn)態(tài)雷諾數(shù)定義為Re=dup/υ，計算得到該工況下的穩(wěn)態(tài)雷諾數(shù)為Re=1.03，文獻計算值為Re=1.03，兩者一致，本程序的正確性得以驗證。

圖8 數(shù)值計算與參考結(jié)果驗證[25]Fig. 8 Validation of numerical with reference case[25]

1.3.2 網(wǎng)格無關(guān)性驗證

本文使用不同三種網(wǎng)格分辨率模擬了圖1 所示的撲翼運動。計算參數(shù)θ =0°、Re=100、St=2Az f/U∞=0.15， 其中Az為拍動幅值。三個不同的網(wǎng)格分辨率分別為200×200×200 （Case A）、245×245×245 （Case B）、288×288×288 （Case C），三種網(wǎng)格分辨率對應(yīng)的撲翼所占的網(wǎng)格點數(shù)分別為4 513、5 528 和6 498。

圖9 為三種網(wǎng)格分辨率下的翼尖處渦量云圖。從圖中可以看出，渦量云圖基本一致。但隨著網(wǎng)格分辨率的增加，在撲翼尾跡處的旋渦分離趨勢更加明顯，渦的細節(jié)也更清晰。

圖9 翼尖渦量云圖Fig. 9 Vorticities contour at wing-tip

圖10 為不同網(wǎng)格分辨率下的升力系數(shù)與阻力系數(shù)隨時間變化曲線圖，其中t為無量綱時間：當前時間步數(shù)與一個周期時間步長的比值，撲翼上下?lián)鋭右淮螢橐粋€周期。可以看出，不同的網(wǎng)格分辨率獲得的結(jié)果基本一致。

圖10 網(wǎng)格無關(guān)結(jié)果圖Fig. 10 Grid independence results

本文計算6 個撲動周期，一個撲動周期對應(yīng)1 200個LBM 步長，共7 200 個LBM 時間步長，對于A、B、C 三種工況，計算時間分別為3.68 h、3.86 h、4.20 h。三種工況的計算均采用GPU 加速，因計算時間較短，本文在之后的計算中均采取288×288×288 的網(wǎng)格分辨率，以捕捉到更精細的流場細節(jié)。

2 結(jié)果與討論

2.1 迎角對升阻力系數(shù)的影響

如圖11 所示，以 θ=0°為例，在一個周期內(nèi)，撲翼在t=0 時 向下拍動，此時撲動角速度 ω最 大；t=T/4時轉(zhuǎn)為向上拍動，此時 ω=0；t=T/2時向上拍動，此時ω最 大；t=3T/4 時轉(zhuǎn)為向下拍動，此時ω =0。

圖11 一個周期內(nèi)撲翼運動示意圖Fig. 11 Schematic of flapping-wing movement in a period

圖12 為不同迎角下升力系數(shù)CL、阻力系數(shù)CD隨t變化的曲線。曲線取自撲翼運動的第5、6 個周期，此時流場已經(jīng)處于穩(wěn)定狀態(tài)，各個周期內(nèi)升阻力系數(shù)變化穩(wěn)定一致。如圖12（a）所示，在一個周期內(nèi)隨著時間推進，CL在撲翼下?lián)溥^程中達到最大值，在撲翼上撲過程中達到最小值。隨著 θ 的 增大，CL總體呈上升 趨 勢，直 到 θ=60°時 ，CL開 始 下 降。當 θ=40°與60°時 ，CL最小值接近于0。如圖12（b）所示，當θ=0°時，撲翼的迎風面積極小，CD接近于0，且CD變化較小。當迎角逐漸增大時，撲翼的迎風面積逐漸增大，CD呈上升趨勢，并在撲翼下?lián)溥^程中達到最大值，在撲翼上撲過程中達到最小值。

圖12 升阻力系數(shù)隨時間變化曲線Fig. 12 Variation curve of lift and drag coefficient with time

圖13 分別為平均升力系數(shù)CˉL、平均阻力系數(shù)CˉD以及升阻比CˉL/CˉD隨迎角變化曲線圖，此處，平均升阻力系數(shù)CˉL、CˉD為一個周期內(nèi)的算術(shù)平均值。從圖13（a）中可以看出，CˉL隨迎角的增大呈先增大后減小 的 趨 勢，當 θ ≥40°時 ，CˉL開 始 下 降。從 圖13（b）中可以看出，CˉD隨著迎角的增大逐漸增大。從圖13（c）中可以看出，CˉL/CˉD隨迎角增大呈先增大后減小的趨勢，當 θ=10°時 ，CˉL/CˉD最大，此時撲翼的氣動性能達到最佳。當 θ=60°時 ，CˉL/CˉD較 θ=0°時更小，該迎角下?lián)湟磉\動氣動性能反而下降。

圖13 平均升阻力系數(shù)隨迎角變化曲線Fig. 13 Variation curve of average lift and drag coefficient with angle of attack

下文通過壓力場和流場結(jié)構(gòu)對撲翼運動的氣動特性機理進行分析。

2.2 氣動特性機理分析

2.2.1 壓力場分析（圖1（c））受力面減小，反而使CL變小。

圖14 t =0上下翼面壓力云圖Fig. 14 Pressure contour of lower and upper surface of wing ont=0

圖15 t =T/2上下翼面壓力云圖Fig. 15 Pressure contour of lower and upper surface of wing ont=T/2

同樣地，分析圖16（b、d），從翼根到翼尖方向，Δp絕 對值也逐漸增大，但是隨著迎角的增大， Δp絕對值明顯比下?lián)潆A段小，定量說明了圖12（a）中高迎角時CL最小值接近于0 的情況。兩幅正壓曲線在翼尖處下降趨勢最明顯，說明在翼尖處壓力變化最大。

圖16 壓力系數(shù)和壓力差沿翼型弦長方向變化曲線圖（t =0 ，t =T/2）Fig. 16 Pressure coefficient and pressure difference distribution along chord length direction of airfoil on t =0 and t=T/2

2.2.2 流場及渦結(jié)構(gòu)分析

圖17 為不同迎角下?lián)湟磉\動渦系結(jié)構(gòu)隨時間的演變，圖中Ω表示LBM 格子單位下的渦量。從圖中可以觀察到前緣渦，撲翼拍動過程中，上翼面前緣開始產(chǎn)生前緣渦。前緣渦具有附著于撲翼表面不脫落的特性，附著于翼面的面積越大，提供的升力越大，而一旦當前緣渦從撲翼上脫落，升力系數(shù)將會大幅下降。可以看出，當t=0與t=T/4時，在翼尖處前緣渦最大，此時，翼尖處獲得的升力最大。此外，在同一迎角下，t=T/4 ～ 3T/4時，附著在翼面的前緣渦不斷從翼面前緣脫落，此時CL小于該迎角下整個周期的CˉL；t=3T/4 至 下一個周期的t=T/4，附著在翼面的前緣渦開始生成，此時CL大于該迎角下整個周期的CˉL，這解釋了圖12（a）所示的升力變化趨勢。

圖17 不同迎角下?lián)湟磉\動渦系結(jié)構(gòu)隨時間的演變Fig. 17 Variation of vortices structure of flapping-wing for different angle of attack with time

圖18 為不同迎角下，t=0時翼尖處的流線與壓力圖。從圖18 中可以看出，在拍動的過程中，撲翼下方壓強較大，上方壓強較小，空氣在翼尖處上下壓差更大，從撲翼下表面經(jīng)由翼尖繞至上表面，形成翼尖渦。隨著迎角的增大，翼尖渦逐漸由翼面處攀升至翼面上方，強度與影響范圍逐漸增大，從而導(dǎo)致渦致阻力的增加，使得CD逐漸增大。

圖18 t =0不同迎角流線與壓力圖Fig. 18 Streamline and pressure contour for different angle of attack on t=0

如圖17 所示，在t=0時，即撲翼處于下拍階段，翼尖渦從翼尖處分離并順著來流速度方向傳輸，導(dǎo)致t=0 ～T/2時 渦致阻力下降，此時撲翼CD下降；在t=T/2時，即撲翼處于上拍階段，翼尖渦又開始快速生成，使得t=T/2 ～T時渦致阻力上升，此時撲翼CD上升。

如圖17 所示，隨著迎角的增大，前緣渦附著于翼面的面積明顯增大，且與翼尖渦逐漸融合，并在撲翼尾跡處逐漸形成了封閉的渦環(huán)，不易耗散，且渦環(huán)與撲翼間的相互作用，使得渦環(huán)中的能量得以重復(fù)利用，使撲翼維持在一個高升力狀態(tài)。 θ=60°時，撲翼尾跡處的旋渦開始變得不穩(wěn)定并破碎形成一個個細小的渦，流場結(jié)構(gòu)紊亂，旋渦的能量被分解耗散，此時，撲翼的氣動性能下降。

3 結(jié) 論

本文結(jié)合LBM 與Level-set 動邊界識別方法，采用GPU 加速，建立了三維含有運動邊界流場的高效數(shù)值模擬方法并開發(fā)了相應(yīng)程序。基于此，開展了不同迎角下?lián)湟磉\動氣動特性機理分析。研究結(jié)果表明：

1） 隨著迎角的增加，撲翼的平均升力系數(shù)先增大后減小，平均阻力系數(shù)逐漸增大。升阻比在10°時達到最大。此時，撲翼能夠獲得最佳氣動性能。

2） θ =0°～40°范圍內(nèi)，隨著迎角的增加，撲翼上下翼面正負壓區(qū)面積增大，上下翼面壓差增大，撲翼獲得的升力增加；翼尖處壓差最大，給撲翼提供了主要的升力。

3） 隨著迎角的增大，前緣渦附著于翼面的面積明顯增大，撲翼后方脫落的渦旋也較難耗散，提高了撲翼的升力；同時，翼尖渦的強度與影響范圍逐漸增大，導(dǎo)致渦致阻力增加，阻力系數(shù)增加。

LBM 結(jié)合Level-Set 方法可精確捕捉具有復(fù)雜外形的運動邊界，同時，采用GPU 加速能夠獲得極高的加速比。本文的數(shù)值方法可為撲翼運動提供高效、可靠的研究手段。