基于浸入邊界法的撲翼飛行數(shù)值模擬1)

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

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

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

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

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

1 數(shù)值方法

1.1 浸入邊界法

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

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

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

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

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

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

1.2 物理模型

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

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

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

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

1.3 氣動參數(shù)

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

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

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

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

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

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

2 算例驗(yàn)證

在開始大量的計(jì)算前，需要先對IBM 程序進(jìn)行驗(yàn)證，并和現(xiàn)有文獻(xiàn)結(jié)果進(jìn)行對比，以確保程序的準(zhǔn)確性。這里和兩種不同方法的數(shù)值模擬結(jié)果進(jìn)行了對比。選取了文獻(xiàn)[7]使用的一個以單元為中心的基于壓力的有限體積流動求解器的模擬結(jié)果以及文獻(xiàn)[18]的低雷諾數(shù)下有限差分法求解結(jié)果對本文的方法進(jìn)行驗(yàn)證。文獻(xiàn)[7]的模擬參數(shù)為：Re=1100，St=0.45，ha=1.0（以最大弦長為無量綱長度），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é)果對比時，也對網(wǎng)格大小的影響進(jìn)行了研究。在模擬時隨著拉格朗日網(wǎng)格點(diǎn)間距Δx的縮小，同一大小模型下拉格朗日點(diǎn)數(shù)目增加，更能準(zhǔn)確地描述模型的形狀，并且相對應(yīng)的歐拉網(wǎng)格數(shù)增加，模擬結(jié)果更加精確。但是這種增益不是線性的，即拉格朗日網(wǎng)格增加到一定大小時，再加密網(wǎng)格對模擬結(jié)果的影響變得非常小。這里選取了4 種Δx在同一邊界條件下進(jìn)行計(jì)算，圖3 為浸入邊界法計(jì)算的一個拍動周期的瞬時升力系數(shù)和瞬時推力系數(shù)與文獻(xiàn)[7]的模擬結(jié)果對比圖。在瞬時升力系數(shù)圖中，Δx=0.05時結(jié)果略有偏差，Δx=0.012 5 和Δx=0.006 25時的瞬時升力系數(shù)基本重合并和文獻(xiàn)[7]的結(jié)果很好地吻合。對于瞬時推力系數(shù)，在Δx=0.05時結(jié)果誤差較大，Δx=0.012 5 和Δx=0.006 25的結(jié)果較接近并和文獻(xiàn)吻合?？梢钥闯鐾屏ο禂?shù)的精確性對網(wǎng)格大小有較高的依賴性，這與文獻(xiàn)[19]在固定翼型繞流驗(yàn)證中的結(jié)論一致。表1 為不同網(wǎng)格模擬結(jié)果的具體數(shù)據(jù)，這里重點(diǎn)關(guān)注了推力系數(shù)。在Δx=0.012 5 時，推力系數(shù)的相對誤差僅為0.021%，而計(jì)算時長大約為Δx=0.006 25時的1/20，并且網(wǎng)格繼續(xù)加密到Δx=0.006 25時精確度提升不大。綜合考慮，后面的計(jì)算均取Δx=0.012 5。

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

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

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

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

3 二維模擬分析和討論

3.1 升沉運(yùn)動

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

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

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

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

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

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

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

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

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

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

對于升沉運(yùn)動耦合俯仰運(yùn)動，簡稱為撲翼運(yùn)動。在本節(jié)研究最大俯仰角對推進(jìn)效率的影響，以便確定撲翼飛行的最佳俯仰角。模擬參數(shù)為：St=0.3，ha=0.5，f=0.3，升沉運(yùn)動和俯仰運(yùn)動的相位差φα=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°是一個較好的選擇。在李寧宇等[20]和呂元博等[21]對翼型俯仰角的研究中，也均在20°時表現(xiàn)出較高的氣動性能。

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

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

改變升沉運(yùn)動與俯仰運(yùn)動的相位差會帶來不同的氣動特性。這里對升沉運(yùn)動和俯仰運(yùn)動的相位差進(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ù)對撲翼運(yùn)動氣動特性的影響

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

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

4 三維模擬分析和討論

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

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

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

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

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

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

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

圖14 k=0.35 時一個拍動周期的翅膀展向壓力云圖Fig.14 Spanwise pressure contour for k=0.35

5 結(jié)論

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

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

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

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

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