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

(1. 西安交通大學(xué) 航天航空學(xué)院，機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710049；2. 陜西省先進(jìn)飛行器服役環(huán)境與控制重點(diǎn)實(shí)驗(yàn)室，西安 710049)

0 引 言

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

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

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

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

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

1.1 計(jì)算模型

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

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

翼面的展長(zhǎng)L=4.90 cm ，平均弦長(zhǎng)D=1.24 cm，厚度s=0.14 cm， θ為迎角。兩個(gè)撲翼分別繞各自的翼根作上下?lián)鋭?dòng)運(yùn)動(dòng)，撲動(dòng)方程為：

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

雷諾數(shù)Re=UrefD/υ =100， 無(wú)量綱D=18，其中Uref=4αmaxfr，Uref為撲翼運(yùn)動(dòng)參考速度，r為撲翼的旋轉(zhuǎn)半徑。本文分別計(jì)算了 θ=0°～60°時(shí)不同迎角的流場(chǎng)情況。流體的計(jì)算區(qū)域?yàn)長(zhǎng)x×Ly×Lz，每個(gè)方向取1 6D，即 288×288×288，總計(jì)2 389 萬(wàn)網(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）求解流場(chǎng)，其離散方程如下：

圖2 D3Q19 模型Fig. 2 D3Q19 model

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

1.2.2 LBM 移動(dòng)固壁邊界條件

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

對(duì)于移動(dòng)的固體壁面，本文采用非平衡態(tài)外推格式[23]確定固體邊界。如圖3 所示，非平衡態(tài)外推格式將固體節(jié)點(diǎn)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表示與壁面相交一個(gè)網(wǎng)格間距中處于流體區(qū)域的比例：

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

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

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

公式如下所示:

1.2.3 Level-Set 動(dòng)邊界識(shí)別方法

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

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

圖5 Level-Set 動(dòng)邊界識(shí)別方法Fig. 5 Schematic of identification method of moving boundary of Level-Set

1.2.4 GPU 程序移植

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

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

圖6 計(jì)算流程圖Fig. 6 Flow chart of numerical method

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

1.3 程序驗(yàn)證

1.3.1 正確性驗(yàn)證

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

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

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

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

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

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

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

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

圖10 為不同網(wǎng)格分辨率下的升力系數(shù)與阻力系數(shù)隨時(shí)間變化曲線圖，其中t為無(wú)量綱時(shí)間：當(dāng)前時(shí)間步數(shù)與一個(gè)周期時(shí)間步長(zhǎng)的比值，撲翼上下?lián)鋭?dòng)一次為一個(gè)周期?？梢钥闯?，不同的網(wǎng)格分辨率獲得的結(jié)果基本一致。

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

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

2 結(jié)果與討論

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

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

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

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

圖12 升阻力系數(shù)隨時(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為一個(gè)周期內(nèi)的算術(shù)平均值。從圖13（a）中可以看出，CˉL隨迎角的增大呈先增大后減小 的 趨 勢(shì)，當(dāng) θ ≥40°時(shí) ，CˉL開(kāi) 始 下 降。從 圖13（b）中可以看出，CˉD隨著迎角的增大逐漸增大。從圖13（c）中可以看出，CˉL/CˉD隨迎角增大呈先增大后減小的趨勢(shì)，當(dāng) θ=10°時(shí) ，CˉL/CˉD最大，此時(shí)撲翼的氣動(dòng)性能達(dá)到最佳。當(dāng) θ=60°時(shí) ，CˉL/CˉD較 θ=0°時(shí)更小，該迎角下?lián)湟磉\(yùn)動(dòng)氣動(dòng)性能反而下降。

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

下文通過(guò)壓力場(chǎng)和流場(chǎng)結(jié)構(gòu)對(duì)撲翼運(yùn)動(dòng)的氣動(dòng)特性機(jī)理進(jìn)行分析。

2.2 氣動(dòng)特性機(jī)理分析

2.2.1 壓力場(chǎng)分析（圖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絕 對(duì)值也逐漸增大，但是隨著迎角的增大， Δp絕對(duì)值明顯比下?lián)潆A段小，定量說(shuō)明了圖12（a）中高迎角時(shí)CL最小值接近于0 的情況。兩幅正壓曲線在翼尖處下降趨勢(shì)最明顯，說(shuō)明在翼尖處壓力變化最大。

圖16 壓力系數(shù)和壓力差沿翼型弦長(zhǎng)方向變化曲線圖（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 流場(chǎng)及渦結(jié)構(gòu)分析

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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