赫連勃勃 張 權(quán) 周 錕

(武漢科技大學(xué)省部共建耐火材料與冶金國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430081)

鳥類和昆蟲具有極高的飛行技巧，例如前飛、懸停和滑翔。自然界中鳥類和昆蟲等飛行生物均采用撲翼飛行方式[1]。撲翼飛行既能像固定翼飛行器那樣快速前行，也能像旋翼飛行器那樣在空中懸停。鳥類和昆蟲這些超高的飛行技巧讓人們堅(jiān)信撲翼飛行器在環(huán)境勘察、軍情偵察等場(chǎng)景具有巨大的應(yīng)用潛力，與撲翼相關(guān)的研究受到了國(guó)內(nèi)外學(xué)術(shù)和工程應(yīng)用領(lǐng)域的高度關(guān)注[2]。

國(guó)內(nèi)外眾多學(xué)者對(duì)撲翼飛行進(jìn)行了大量的研究。國(guó)內(nèi)：葉軍科等[3]求解了加入擬壓縮項(xiàng)的不可壓納維–斯托克斯方程，研究了不同翼型厚度和雷諾數(shù)對(duì)翼型非定常氣動(dòng)特性的影響；韓偉等[4]研究了矩形波形、三角波形對(duì)俯仰運(yùn)動(dòng)能量獲取的影響，同時(shí)研究了斯特勞哈爾數(shù)、俯仰角、俯仰軸位置對(duì)能量獲取效率的影響；沈志偉等[5]應(yīng)用混合笛卡爾網(wǎng)格的方法，對(duì)單體翼型以及前后擺放的雙翼型的振幅影響進(jìn)行了數(shù)值模擬；劉葳興等[6]應(yīng)用浸入邊界法對(duì)三維柔性撲翼進(jìn)行了模擬，系統(tǒng)地研究了撲翼參數(shù)對(duì)三維柔性翼推進(jìn)性能的影響。國(guó)外：Pedro 等[7]對(duì)NACA 0012 水翼進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分并做了繞流數(shù)值模擬，研究了升沉頻率和俯仰頻率等參數(shù)對(duì)水翼氣動(dòng)性能的影響；Lee 等[8]對(duì)果蠅的三維流場(chǎng)進(jìn)行了數(shù)值模擬，研究了迎角對(duì)果蠅飛行速度的影響，并進(jìn)一步研究了尾流和翅膀的相互作用；Amiralaei 等[9]基于有限體積法通過求解非定常不可壓納維–斯托克斯方程對(duì)撲翼飛行進(jìn)行數(shù)值模擬，分析了翼型周圍的流場(chǎng)和渦流，研究了俯仰振幅、相位角等參數(shù)對(duì)撲翼飛行的影響。實(shí)驗(yàn)方面：王肇等[10]對(duì)純升沉運(yùn)動(dòng)進(jìn)行了水洞實(shí)驗(yàn)研究，實(shí)驗(yàn)結(jié)果展示了翼型升沉運(yùn)動(dòng)尾流特性和斯特勞哈爾數(shù)的關(guān)系，并和數(shù)值模擬結(jié)果做對(duì)比；邵立民等[11]通過風(fēng)洞實(shí)驗(yàn)，對(duì)微型撲翼飛行器的拍動(dòng)頻率、來流速度、翼型的彎曲度等進(jìn)行了研究；Anderson 等[12]在雷諾數(shù)1100 時(shí)對(duì)諧波振蕩箔片的推力進(jìn)行研究，并使用數(shù)字粒子圖像測(cè)速儀得到可視化數(shù)據(jù)；在雷諾數(shù)40 000 下獲得翼型力和功率數(shù)據(jù)。在三維模擬方面，朱霖霖等[13]模擬了撲翼鳥在黏性流體中自由飛行，研究了起飛俯仰角度以及尾巴和身體的夾角對(duì)整個(gè)飛行過程中俯仰角的影響；高赫等[14]建立了更加接近真實(shí)的三維鳥翅膀模型，編寫了三個(gè)自由度 Fluent 耦合控制函數(shù)，使用動(dòng)網(wǎng)格實(shí)時(shí)更新的研發(fā)方法，對(duì)翅膀的升力和阻力進(jìn)行計(jì)算，揭示了自由度對(duì)撲翼空氣動(dòng)力學(xué)性能的影響。

雖然國(guó)內(nèi)外學(xué)者對(duì)撲翼飛行的研究很多，但是對(duì)撲翼飛行推進(jìn)的研究較少，并且以往要獲得撲翼及周圍流場(chǎng)信息需采用貼體網(wǎng)格[15]，在復(fù)雜的流固耦合動(dòng)網(wǎng)格問題上，傳統(tǒng)的貼體網(wǎng)格不僅要建立復(fù)雜的網(wǎng)格，還要在每個(gè)時(shí)間步通過變化網(wǎng)格來適應(yīng)改變的邊界條件，不僅增加了計(jì)算難度[16]，還可能對(duì)模擬精確性造成影響。浸入邊界法（immersed boundary method,IBM）能夠很好地處理這一類具有復(fù)雜外形的動(dòng)網(wǎng)格問題。本文將撲翼飛行簡(jiǎn)化為二維運(yùn)動(dòng)與簡(jiǎn)單三維模型，從二維和三維兩個(gè)角度入手，通過求解納維–斯托克斯方程，應(yīng)用浸入邊界法數(shù)值模擬了低雷諾數(shù)下(Re=1100)撲翼運(yùn)動(dòng)的非定常流動(dòng)，系統(tǒng)地研究了純升沉運(yùn)動(dòng)、耦合升沉運(yùn)動(dòng)與俯仰運(yùn)動(dòng)時(shí)，升沉振幅、升沉頻率、俯仰頻率、俯仰角、相位角等參數(shù)對(duì)撲翼氣動(dòng)特性的影響，以及升沉翼型周圍渦流生成和脫落規(guī)律。三維方面對(duì)翅膀拍動(dòng)的時(shí)間非對(duì)稱性進(jìn)行了研究，系統(tǒng)地對(duì)翅膀快拍慢回特性進(jìn)行了模擬和分析。

1 數(shù)值方法

1.1 浸入邊界法

在浸入邊界法模擬中，對(duì)于流體的運(yùn)動(dòng)采用均勻的歐拉節(jié)點(diǎn)描述，而翼型邊界的運(yùn)動(dòng)是采用拉格朗日節(jié)點(diǎn)描述。浸入邊界法第一步是先構(gòu)造流場(chǎng)控制方程，假設(shè)整個(gè)計(jì)算區(qū)域?yàn)榕ｎD流體不可壓縮流動(dòng)，流體區(qū)域與固體區(qū)域全部用歐拉坐標(biāo)系下的納維–斯托克斯方程和質(zhì)量守恒方程描述

式中，u是流體的流速，p是流場(chǎng)的壓力，ρf是流體的密度，ν 是流體的運(yùn)動(dòng)黏度，?為哈密頓算子，fIBM為固體對(duì)流體的作用力，fIBM可以表示為

式中，B為插值函數(shù)，W為控制收斂的權(quán)重因子，F(xiàn)IBM為翼型表面拉格朗日離散點(diǎn)上的作用力，fIBM可由翼型表面拉格朗日離散點(diǎn)上的作用力FIBM經(jīng)過插值求得，而FIBM又可由翼型與流體在界面上的速度差與時(shí)間步長(zhǎng)的比率來確定，即

式中，xi，(i=1,2···n) 表示n個(gè)歐拉坐標(biāo)點(diǎn)坐標(biāo)，XI,(I=1,2···N)表示拉格朗日節(jié)點(diǎn)坐標(biāo)；UX為拉格朗日點(diǎn)上的速度，Δt為計(jì)算時(shí)間步長(zhǎng)，為插值所得拉格朗日網(wǎng)格點(diǎn)上的速度，使用插值函數(shù)求解的表達(dá)式為

式中，為歐拉網(wǎng)格點(diǎn)上的速度，h為歐拉網(wǎng)格的尺寸，插值函數(shù)δh為

式中，離散狄拉克函數(shù)ψ(r)的取值為

1.2 物理模型

本文的模擬從二維和三維兩個(gè)方面展開，撲翼鳥的三維拉格朗日離散模型如圖1(a)所示，模型由身體、左翼和右翼三部分組成。翅膀平均弦長(zhǎng)0.7，翼展2.5，采用無厚度剛性平板近似代替，忽略了翅膀厚度及彎度的影響。撲翼鳥身體由圓柱、圓錐及扇形尾部組成，身體長(zhǎng)度1.7。模型共由3115 個(gè)拉格朗日離散點(diǎn)組成，離散點(diǎn)間距Δx為0.04，翅膀可以繞x=0 軸和y=0 軸做拍動(dòng)和轉(zhuǎn)動(dòng)。因?yàn)楸疚牡哪M屬于低雷諾數(shù)的三維繞流問題，對(duì)于槽道大小的選擇，參考了周錕等[17]的先驗(yàn)誤差估計(jì)法，既保證了結(jié)果的準(zhǔn)確性，又不會(huì)造成計(jì)算資源的浪費(fèi)。槽道示意圖如圖1(b)，槽道的物理尺寸x×y×z為20×40×20，槽道的歐拉網(wǎng)格節(jié)點(diǎn)數(shù)為x/Δx，y/Δx及z/Δx。在槽道的頂部和底部設(shè)置速度邊界來模擬穩(wěn)定的空氣來流，來流速度U=1，流場(chǎng)y方向采用周期性邊界條件。對(duì)于二維模擬，這里將翅膀的拍動(dòng)和轉(zhuǎn)動(dòng)簡(jiǎn)化為升沉運(yùn)動(dòng)與俯仰運(yùn)動(dòng)，如圖2 所示，從右往左為不同時(shí)刻翼型的運(yùn)動(dòng)狀態(tài)，其中俯仰運(yùn)動(dòng)繞著靠近前緣1/4 弦長(zhǎng)處旋轉(zhuǎn)。選取NACA 0012 翼型作為二維模型，在Profili 軟件中獲取其坐標(biāo)信息，并在Matlab 中用其坐標(biāo)生成拉格朗日網(wǎng)格。二維翼型弦長(zhǎng)為1，拉格朗日離散點(diǎn)間距Δx根據(jù)精度要求調(diào)整。二維翼型升沉運(yùn)動(dòng)和俯仰運(yùn)動(dòng)方程為

圖1 三維模型和模擬區(qū)域示意圖Fig.1 Schematic diagram of 3D model and simulation area

圖2 升沉運(yùn)動(dòng)和俯仰運(yùn)動(dòng)Fig.2 Heave motion and pitch motion

式中，ha為升沉幅度，f為翼型的拍動(dòng)頻率，αa為俯仰運(yùn)動(dòng)的俯仰角，φα為升沉運(yùn)動(dòng)和俯仰運(yùn)動(dòng)的相位差。

1.3 氣動(dòng)參數(shù)

在撲翼運(yùn)動(dòng)模擬中，流場(chǎng)的雷諾數(shù)公式定義為

式中，U為流場(chǎng)的來流速度；ρf為流體的密度；μ為流體的黏性系數(shù)；d為特征長(zhǎng)度，取翅膀或翼型的弦長(zhǎng)。

在二維模擬中，翼型的斯特勞哈爾數(shù)、升力系數(shù)、推力系數(shù)分別定義為

式中，f為翼型的拍動(dòng)頻率，2ha為翼型后緣的總位移量[12]，U為來流速度，L為翼型所受的升力，D為翼型受到的阻力，在翼型的推進(jìn)研究中，從推力進(jìn)行分析更為方便，推力系數(shù)可以看作阻力系數(shù)的相反數(shù)。Ac為翼型或翅膀的投影面積。

把單位時(shí)間內(nèi)翼型克服流體所做的功定義為瞬時(shí)輸入功率，瞬時(shí)輸入功率、瞬時(shí)輸入功率系數(shù)及推進(jìn)效率計(jì)算公式分別為

式中，L(t)為翼型垂直方向的升力，(t) 為升沉運(yùn)動(dòng)的速度，M(t)為翼型升力和阻力在轉(zhuǎn)動(dòng)軸上產(chǎn)生的力矩，α˙(t) 為俯仰運(yùn)動(dòng)的轉(zhuǎn)動(dòng)角速度。CTm和CPm為一個(gè)穩(wěn)定周期內(nèi)的平均推力系數(shù)和平均輸入功率系數(shù)。

2 算例驗(yàn)證

在開始大量的計(jì)算前，需要先對(duì)IBM 程序進(jìn)行驗(yàn)證，并和現(xiàn)有文獻(xiàn)結(jié)果進(jìn)行對(duì)比，以確保程序的準(zhǔn)確性。這里和兩種不同方法的數(shù)值模擬結(jié)果進(jìn)行了對(duì)比。選取了文獻(xiàn)[7]使用的一個(gè)以單元為中心的基于壓力的有限體積流動(dòng)求解器的模擬結(jié)果以及文獻(xiàn)[18]的低雷諾數(shù)下有限差分法求解結(jié)果對(duì)本文的方法進(jìn)行驗(yàn)證。文獻(xiàn)[7]的模擬參數(shù)為：Re=1100，St=0.45，ha=1.0（以最大弦長(zhǎng)為無量綱長(zhǎng)度），f=0.225，αa=30°，φα=90°。文獻(xiàn)[18]的模擬參數(shù)為：Re=1100，St=0.3，ha=0.4，f=0.375，αa=20°，φα=90°。

在和文獻(xiàn)[7]的模擬結(jié)果對(duì)比時(shí)，也對(duì)網(wǎng)格大小的影響進(jìn)行了研究。在模擬時(shí)隨著拉格朗日網(wǎng)格點(diǎn)間距Δx的縮小，同一大小模型下拉格朗日點(diǎn)數(shù)目增加，更能準(zhǔn)確地描述模型的形狀，并且相對(duì)應(yīng)的歐拉網(wǎng)格數(shù)增加，模擬結(jié)果更加精確。但是這種增益不是線性的，即拉格朗日網(wǎng)格增加到一定大小時(shí)，再加密網(wǎng)格對(duì)模擬結(jié)果的影響變得非常小。這里選取了4 種Δx在同一邊界條件下進(jìn)行計(jì)算，圖3 為浸入邊界法計(jì)算的一個(gè)拍動(dòng)周期的瞬時(shí)升力系數(shù)和瞬時(shí)推力系數(shù)與文獻(xiàn)[7]的模擬結(jié)果對(duì)比圖。在瞬時(shí)升力系數(shù)圖中，Δx=0.05時(shí)結(jié)果略有偏差，Δx=0.012 5 和Δx=0.006 25時(shí)的瞬時(shí)升力系數(shù)基本重合并和文獻(xiàn)[7]的結(jié)果很好地吻合。對(duì)于瞬時(shí)推力系數(shù)，在Δx=0.05時(shí)結(jié)果誤差較大，Δx=0.012 5 和Δx=0.006 25的結(jié)果較接近并和文獻(xiàn)吻合。可以看出推力系數(shù)的精確性對(duì)網(wǎng)格大小有較高的依賴性，這與文獻(xiàn)[19]在固定翼型繞流驗(yàn)證中的結(jié)論一致。表1 為不同網(wǎng)格模擬結(jié)果的具體數(shù)據(jù)，這里重點(diǎn)關(guān)注了推力系數(shù)。在Δx=0.012 5 時(shí)，推力系數(shù)的相對(duì)誤差僅為0.021%，而計(jì)算時(shí)長(zhǎng)大約為Δx=0.006 25時(shí)的1/20，并且網(wǎng)格繼續(xù)加密到Δx=0.006 25時(shí)精確度提升不大。綜合考慮，后面的計(jì)算均取Δx=0.012 5。

表1 網(wǎng)格對(duì)浸入邊界法模擬結(jié)果的影響Table 1 Influence of mesh on IBM simulation results

圖3 一個(gè)周期的瞬時(shí)升力系數(shù)和推力系數(shù)與文獻(xiàn)[7]的結(jié)果對(duì)比Fig.3 Comparison of instantaneous lift coefficient and thrust coefficient for one cycle with Ref.[7]

圖4 為浸入邊界法計(jì)算的一個(gè)拍動(dòng)周期瞬時(shí)升力系數(shù)和瞬時(shí)推力系數(shù)與文獻(xiàn)[18]的模擬結(jié)果對(duì)比。文獻(xiàn)[18]的結(jié)果平均推力系數(shù)為0.35，浸入邊界法模擬的平均推力系數(shù)為0.34，相對(duì)誤差為3.09%。并且本文重點(diǎn)研究撲翼的推進(jìn)效率。從模擬結(jié)果看，在和兩種不同的數(shù)值方法對(duì)比中，浸入邊界法模擬結(jié)果和文獻(xiàn)結(jié)果變化趨勢(shì)相同，誤差較小，模擬結(jié)果可信度較高。

圖4 一個(gè)周期的瞬時(shí)升力系數(shù)和推力系數(shù)與文獻(xiàn)[18]的結(jié)果對(duì)比Fig.4 Comparison of instantaneous lift coefficient and thrust coefficient for one cycle with Ref.[18]

3 二維模擬分析和討論

3.1 升沉運(yùn)動(dòng)

首先對(duì)純升沉運(yùn)動(dòng)進(jìn)行了研究，這一部分模擬了Re=1100，0.2≤St≤0.6，0.1≤ha≤0.5 的情況，得到的推力系數(shù)等值線圖和推進(jìn)效率等值線如圖5 所示。

圖5 升沉運(yùn)動(dòng)平均推力系數(shù)等值線圖和推進(jìn)效率等值線圖Fig.5 Contour map of average thrust coefficient and propulsion efficiency of heave motion

在升沉運(yùn)動(dòng)中，St不變時(shí)，推力系數(shù)隨著振幅的增加而減小，而當(dāng)同一振幅時(shí)，推力系數(shù)隨著St的增大（f增大）而增大，推進(jìn)效率的峰值主要集中在0.3≤St≤0.4，0.2≤ha≤0.3 時(shí)。

特別地，對(duì)升沉運(yùn)動(dòng)中同一St時(shí)低頻率高振幅（情形（i））及高頻率低振幅（情形（ii））兩種特殊情況的壓力云圖和翼型表面壓力進(jìn)行了分析，選取了St=0.3，ha=0.4，f=0.375 和St=0.3，ha=0.1，f=1.5 兩種情況。情形（i）時(shí)一個(gè)周期內(nèi)的瞬時(shí)升力系數(shù)、推力系數(shù)如圖6 所示。

圖6 情形（i）時(shí)瞬時(shí)升力系數(shù)和瞬時(shí)推力系數(shù)Fig.6 Case (i) instantaneous lift coefficient and instantaneous thrust coefficient

圖7 為情形（i）時(shí)下?lián)溥^程的壓力云圖和翼型表面壓力分布圖。在低頻率高振幅時(shí)，渦旋中的壓力遠(yuǎn)低于翼型表面的壓力，由翼型表面壓力分布圖可以看出渦旋脫落時(shí)對(duì)翼型表面的壓力變化有較大的影響。當(dāng)翼型下表面的渦旋向后脫落時(shí)，上表面前緣也形成一個(gè)渦旋，并逐漸向后脫落。當(dāng)前一個(gè)渦剛剛脫落，新生成的渦旋到達(dá)翼型最大厚度時(shí)，升力系數(shù)會(huì)有一個(gè)突增，這與圖6瞬時(shí)升力系數(shù)在0.8T時(shí)的額外高峰相對(duì)應(yīng)。并且當(dāng)渦旋在翼型最大厚度之前時(shí)，對(duì)翼型的推力有較大的貢獻(xiàn)，對(duì)應(yīng)圖6 在0.7T時(shí)刻的瞬時(shí)推力系數(shù)。而當(dāng)渦旋移動(dòng)到翼型最大厚度之后時(shí)，由于這時(shí)渦旋對(duì)翼型的水平推進(jìn)作用力減小甚至到負(fù)值，從而導(dǎo)致推力開始減小。

圖7 情形（i）時(shí)的壓力云圖和翼型表面壓力分布圖Fig.7 Pressure nephogram and pressure distribution of airfoil surface in case (i)

圖8 為情形（ii）時(shí)下?lián)溥^程的壓力云圖和翼型表面壓力分布圖。在高頻率低振幅時(shí)，翼型表面的壓力低于渦旋中的壓力，并且渦旋的體積也比前一種情況要小。當(dāng)上表面的渦旋脫落時(shí)，下表面的渦旋并沒有隨之向后脫落，而是停留在翼型前緣，對(duì)推力有持續(xù)的貢獻(xiàn)。這種情況的推力峰值要高于前一種情況，這與文獻(xiàn)[18]在該方面的研究結(jié)果一致。

圖8 情形(ii)時(shí)的壓力云圖和翼型表面壓力分布圖Fig.8 Pressure nephogram and pressure distribution of airfoil surface in case (ii)

3.2 最大俯仰角對(duì)撲翼運(yùn)動(dòng)氣動(dòng)特性的影響

對(duì)于升沉運(yùn)動(dòng)耦合俯仰運(yùn)動(dòng)，簡(jiǎn)稱為撲翼運(yùn)動(dòng)。在本節(jié)研究最大俯仰角對(duì)推進(jìn)效率的影響，以便確定撲翼飛行的最佳俯仰角。模擬參數(shù)為：St=0.3，ha=0.5，f=0.3，升沉運(yùn)動(dòng)和俯仰運(yùn)動(dòng)的相位差φα=90°，最大俯仰角5°≤αa≤40°。

圖9 為平均推力系數(shù)和推進(jìn)效率隨最大俯仰角的變化。在研究范圍內(nèi)，推力系數(shù)的最大值在10°≤αa≤20°的范圍內(nèi)，推進(jìn)效率的最大值在20°≤αa≤30°的范圍內(nèi)，最大推進(jìn)效率為33.76%。在這種情況下，αa=20°是一個(gè)較好的選擇。在李寧宇等[20]和呂元博等[21]對(duì)翼型俯仰角的研究中，也均在20°時(shí)表現(xiàn)出較高的氣動(dòng)性能。

圖9 平均推力系數(shù)和推進(jìn)效率隨俯仰角變化Fig.9 Variation of average thrust coefficient and propulsion efficiency with pitch angle

3.3 相位差對(duì)撲翼運(yùn)動(dòng)氣動(dòng)特性的影響

改變升沉運(yùn)動(dòng)與俯仰運(yùn)動(dòng)的相位差會(huì)帶來不同的氣動(dòng)特性。這里對(duì)升沉運(yùn)動(dòng)和俯仰運(yùn)動(dòng)的相位差進(jìn)行研究，模擬參數(shù)為：St=0.3，ha=0.5，f=0.3，αa=30°，80°≤φα≤115°，相位角增幅為5°。

研究結(jié)果如圖10 所示，在研究范圍內(nèi)，平均推力系數(shù)隨著相位差增大而增大，推進(jìn)效率的最大值在80°≤φα≤90°之間，最大推進(jìn)效率為30.75%。

圖10 平均推力系數(shù)和推進(jìn)效率隨相位差變化Fig.10 Average thrust coefficient and propulsion efficiency vary with phase difference

3.4 斯特勞哈爾數(shù)對(duì)撲翼運(yùn)動(dòng)氣動(dòng)特性的影響

斯特勞哈爾數(shù)是推力產(chǎn)生機(jī)制中的一個(gè)關(guān)鍵參數(shù)，這里研究斯特勞哈爾數(shù)對(duì)撲翼運(yùn)動(dòng)的影響。0.3≤St≤0.6，0.1≤ha≤0.7，αa=30°，φα=90°，頻率在對(duì)應(yīng)的范圍內(nèi)變化。得出的平均推力系數(shù)和推進(jìn)效率如圖11 所示。推力系數(shù)隨著St的增加而增加，即在同一振幅時(shí)，推力系數(shù)隨著頻率的增大而增大。在同一St時(shí)，推力系數(shù)隨振幅的增大而減小。在同一St時(shí)，推進(jìn)效率表現(xiàn)出相反的變化趨勢(shì)，即振幅越大，推進(jìn)效率越大?？傮w上推進(jìn)效率大小依次為：St=0.4＞St=0.5＞St=0.6 時(shí)。而St=0.3 時(shí)在低振幅和高振幅下都表現(xiàn)出較高的推進(jìn)效率。最大推進(jìn)效率為39.16%。在研究范圍內(nèi)，對(duì)于St的選擇在0.3～0.4 為最佳，對(duì)于振幅的選擇既要保證較高的推力，又要兼顧推進(jìn)效率。

圖11 St 對(duì)撲翼運(yùn)動(dòng)氣動(dòng)性能的影響Fig.11 Effect of St on aerodynamic performance of flapping wing motion

4 三維模擬分析和討論

鳥和昆蟲通過翅膀周期性上下拍動(dòng)提供向上的升力。研究表明，昆蟲的翅膀上拍和下拍的過程中存在明顯的時(shí)間非對(duì)稱性。即在一個(gè)拍動(dòng)周期內(nèi)，翅膀下拍的速度大于上拍的速度，稱為快拍慢回特性。本節(jié)對(duì)昆蟲的這種特性進(jìn)行研究。定義k為下拍的時(shí)間占比，本節(jié)模擬了翅膀拍動(dòng)周期為1 s，拍動(dòng)角45°，雷諾數(shù)Re=1100，k分別為0.35，0.4，0.45 和0.5 的運(yùn)動(dòng)情形。翅膀拍動(dòng)速度分段方程[22]、瞬時(shí)輸入功率及升舉效率分別為

式中，k為下拍的時(shí)間占比，θ 為翅膀拍動(dòng)角，(t)為翅膀拍動(dòng)角速度，P3(t)為三維瞬時(shí)輸入功率，CLm為一個(gè)周期內(nèi)的平均升力系數(shù)，ηL為升舉效率。

k分別為0.35，0.4，0.45 和0.5 時(shí)一個(gè)周期內(nèi)的瞬時(shí)升力系數(shù)、推力系數(shù)和力矩系數(shù)如圖12所示。前半段為下拍過程，后半段為上拍過程。通過分析可知，隨著下拍速度的增快，下拍時(shí)瞬時(shí)升力系數(shù)的峰值也隨之增大。較快的下拍速度帶來大升力的同時(shí)，較慢的上拍速度并未產(chǎn)生很大的負(fù)升力。而推力系數(shù)在下拍過程隨著拍動(dòng)速度增大而增大，但在上拍過程中表現(xiàn)出相反的規(guī)律，并且在k較小時(shí)會(huì)在下拍上拍交替過程中產(chǎn)生較大的瞬時(shí)阻力，這對(duì)飛行會(huì)產(chǎn)生很大的負(fù)面影響。分析瞬時(shí)力矩系數(shù)，在k為0.35 時(shí)的下拍過程中，瞬時(shí)力矩系數(shù)明顯過大，這就會(huì)導(dǎo)致瞬時(shí)輸入功率過大，造成較高的能量損耗。

圖12 不同k 時(shí)的瞬時(shí)升力系數(shù)、推力系數(shù)和力矩系數(shù)Fig.12 Instantaneous lift coefficient,thrust coefficient and moment coefficient at different k

平均升力系數(shù)、推力系數(shù)和升舉效率如圖13所示，平均升力系數(shù)與瞬時(shí)升力系數(shù)的變化規(guī)律一致，都隨著k增大而減小。當(dāng)平均升力越大，代表一個(gè)拍動(dòng)周期內(nèi)提供的總升力更大。平均推力整體上隨著k減小而增大，但由于在k較小時(shí)會(huì)在下拍上拍交替過程中產(chǎn)生較大的瞬時(shí)阻力，抵消了一部分推力，這也導(dǎo)致k=0.35 時(shí)對(duì)平均推力的提升較小。升舉效率的變化規(guī)律相反，隨著k增大而增大，分析是由于在k較小時(shí)，翅膀下拍速度過快，下拍時(shí)的力矩較大，從而導(dǎo)致瞬時(shí)輸入功率過高，所需能量變高，因此升舉效率變低。

圖13 不同k 時(shí)的平均升力系數(shù)、推力系數(shù)和升舉效率Fig.13 Average lift coefficient,thrust coefficient and lifting efficiency at different k

翅膀上拍和下拍時(shí)，不斷有渦旋脫落，同時(shí)翅膀上下表面的壓力也在變化，這是升力的一個(gè)來源。這里對(duì)流場(chǎng)壓力云圖的變化進(jìn)行分析。圖14 為k=0.35 時(shí)翅膀下拍過程和上拍過程的展向壓力云圖。翅膀在下拍時(shí)，上表面壓力逐漸減小，同時(shí)下表面壓力逐漸增大。在0.18T時(shí)上下表面壓力差最大，這時(shí)的瞬時(shí)升力系數(shù)也達(dá)到最大。翅膀繼續(xù)向下拍動(dòng)，這時(shí)有渦旋從翅膀表面脫落，翅膀上表面壓力逐漸升高，升力也隨之降低。下拍過程結(jié)束后，翅膀上表面的低壓渦還未完全脫落。在上拍時(shí)，翅膀表面壓力情況與下拍相反，在水平位置上下表面壓力差最大，渦旋從翅膀下表面脫落，直到過程結(jié)束。在上拍過程中翅膀上下表面的壓力差明顯低于下拍過程。

圖14 k=0.35 時(shí)一個(gè)拍動(dòng)周期的翅膀展向壓力云圖Fig.14 Spanwise pressure contour for k=0.35

5 結(jié)論

本文應(yīng)用浸入邊界法從二維和三維兩個(gè)方面對(duì)撲翼運(yùn)動(dòng)進(jìn)行了一系列研究，在研究范圍內(nèi)得到如下結(jié)論。

（1）浸入邊界法能夠準(zhǔn)確地模擬撲動(dòng)翼型的升力系數(shù)，而推力系數(shù)的精確性對(duì)網(wǎng)格大小有較高的依賴性，當(dāng)Δx縮小到0.012 5 時(shí)滿足所需的精度要求。

（2）對(duì)于升沉運(yùn)動(dòng)，推進(jìn)效率的峰值主要集中在0.3≤St≤0.4，0.2≤ha≤0.3 時(shí)。在升沉運(yùn)動(dòng)中，當(dāng)渦旋在翼型最大厚度之前時(shí)，對(duì)翼型的推力有較大的貢獻(xiàn)。在低頻振動(dòng)時(shí)，渦旋脫落對(duì)翼型表面的壓力變化有較大的影響。在高頻振動(dòng)時(shí)，翼型表面產(chǎn)生的渦旋體積較小，并且渦旋停留在翼型前緣時(shí)間較長(zhǎng)，對(duì)推力有持續(xù)的貢獻(xiàn)。

（3）在對(duì)撲翼運(yùn)動(dòng)的研究中，推進(jìn)效率隨著俯仰角的增大先增大后減小，在25°達(dá)到峰值，最大推進(jìn)效率為33.76%。升沉運(yùn)動(dòng)與俯仰運(yùn)動(dòng)的相位角達(dá)到85°時(shí)，推進(jìn)效率最佳，最大推進(jìn)效率為30.75%。推力隨著St的增大而增大，推進(jìn)效率在0.3≤St≤0.4 時(shí)最高。

（4）在對(duì)三維撲翼拍動(dòng)時(shí)間非對(duì)稱性的研究中，翅膀的升力系數(shù)和推力系數(shù)隨著下拍速度的增大而增大，但是由于所需能量過高，導(dǎo)致升舉效率降低。在展向云圖分析中，翅膀上下表面壓力差隨著拍動(dòng)速度的增大而增大。