俯仰振蕩翼型推阻力轉(zhuǎn)變滯后機制數(shù)值研究

馬德川，邱 展，李高華，王福新

(上海交通大學(xué) 航空航天學(xué)院，上海 200240)

0 引 言

近年來，微型飛行器（Micro Aerial Vehicle, MAV）因其潛在的軍用和民用前景，成為航空領(lǐng)域研究的熱門[1-2]。MAV的雷諾數(shù)很小，通常在100～10 000之間。低雷諾數(shù)下流體的黏性效應(yīng)十分顯著，導(dǎo)致MAV的阻力大幅增加。然而，在相近的雷諾數(shù)下，自然界中的鳥類、昆蟲和魚類依賴其特有的撲翼飛行或游動方式實現(xiàn)自主推進[3]。為優(yōu)化MAV設(shè)計，人們圍繞撲翼的推力特性及其背后的流動機理開展了大量研究。

目前已開展的大部分工作基于大展弦比假設(shè)對繞撲翼流動進行了二維近似。Koochesfahani[4]利用激光多普勒測速（Laser Doppler Velocimetry, LDV）技術(shù)對NACA0012翼型作小幅度俯仰運動的尾跡流場進行了測量，發(fā)現(xiàn)縮減頻率較小時，尾跡形態(tài)為卡門渦街，其中順時針旋轉(zhuǎn)的渦處于逆時針旋轉(zhuǎn)的渦上方；隨著縮減頻率的增大，順時針渦與逆時針渦的縱向位置交換，前者位于后者下方，該尾跡形態(tài)稱為反卡門渦街。Jones等[5]及Von Ellenrieder和Pothos[6]對NACA0012翼型作浮沉運動的實驗研究表明，當(dāng)斯特勞哈爾數(shù)大于某一臨界值時，出現(xiàn)尾跡不對稱現(xiàn)象，即反卡門渦偏離水平線向上方或向下方傾斜。Godoy-Diana等[7]實驗研究了俯仰翼型尾跡不對稱性機制，建立了反卡門渦對稱性破壞準則。Andersen團隊[8-9]對俯仰翼型和浮沉翼型的實驗及數(shù)值研究發(fā)現(xiàn)，在改變振蕩頻率和幅度時，除觀察到卡門渦街（2S BvK）和反卡門渦街（2S RBvK）外，還出現(xiàn)2P型、2P+2S型、4P型、4P+2S型等更復(fù)雜的尾跡形態(tài)，其中“S”表示孤立渦，“P”表示由兩個旋轉(zhuǎn)方向相反的渦組成的渦對，“2”和“4”為每個振蕩周期內(nèi)從翼型后緣脫落到尾跡中的孤立渦或渦對的數(shù)目。

尾跡形態(tài)隨撲翼運動學(xué)、動力學(xué)的演變是撲翼流場的主要特征，其與推力產(chǎn)生的關(guān)系是研究的重點。Von Karman和Burgers[10]在對平板橫向振蕩的研究中發(fā)現(xiàn)，反卡門渦街的“噴流”現(xiàn)象能夠產(chǎn)生推力。Koochesfahani[4]的研究表明，當(dāng)俯仰翼型尾跡形態(tài)由卡門渦街轉(zhuǎn)變?yōu)榉纯ㄩT渦街時，對應(yīng)的尾跡流場由速度損失型變?yōu)樗俣仍黾有突颉皣娏鳌?。他進一步基于積分動量定理證實，反卡門渦街是一種推力產(chǎn)生型尾跡。Jones等[5]在對浮沉翼型的研究中也發(fā)現(xiàn)，隨著無量綱浮沉速度的增大，發(fā)生阻力型尾跡(卡門渦街)到推力型尾跡(反卡門渦街)的轉(zhuǎn)變。Lai和Platzer[11]在此基礎(chǔ)上研究了浮沉翼型尾跡的“噴流”特性，指出最大噴流速度是無量綱浮沉速度的線性函數(shù)。

盡管反卡門渦街的“噴流”理論揭示了撲翼產(chǎn)生推力的流動機制，但隨著研究的深入，人們發(fā)現(xiàn)反卡門渦街的出現(xiàn)與推力產(chǎn)生并不是完全對應(yīng)的。Godoy-Diana等[12]采用粒子圖像測速（Particle Image Velocimetry, PIV）技術(shù)觀察了俯仰翼型的尾跡形態(tài)（Re=1 173），首次發(fā)現(xiàn)隨著斯特勞哈爾數(shù)增大，卡門渦街-反卡門渦街的轉(zhuǎn)變先于阻力-推力轉(zhuǎn)變，原因是翼型撲動產(chǎn)生的推力不足以克服其受到的阻力，Deng等[13]基于計算流體力學(xué)（Computational Fluid Dynamics,CFD）方法對此作了類似的解釋（Re=1 700）；他們均未對阻力的來源進行深入分析。Das等[14]數(shù)值研究了雷諾數(shù)對俯仰翼型推力特性的影響，在10～2 000范圍內(nèi)均觀察到了這一滯后現(xiàn)象，但未就其原因進行討論。Bohl和Koochesfahani[15]對原有積分動量定理作了改進，并用于估計NACA0012俯仰翼型的推力大?。≧e=12 600），發(fā)現(xiàn)在推力為零時尾跡速度型分布已經(jīng)呈“噴流”狀，認為該額外的動量通量用以克服俯仰翼型下游的壓力損失。Ashraf等[16]進一步對NACA0015俯仰翼型的流場進行了PIV測量（Re=2 900），應(yīng)用改進的積分動量定理估計了推力的大小，結(jié)果顯示反卡門渦街的形成不能保證推力的產(chǎn)生，原因是其導(dǎo)致的尾跡流場動量增加被翼型運動引起的壓力損失抵消了。積分動量定理僅涉及位于翼型后緣下游某一特定位置的控制面上的速度分布，因此文獻[15-16]并未考慮控制體內(nèi)部的非定常尾跡流場對推力產(chǎn)生的影響。

目前對撲翼推力特性的研究取得了一些重要結(jié)論。Sarkar和Venkatraman[17]研究發(fā)現(xiàn)，俯仰翼型的時均推力隨縮減頻率的增大而增大，隨平均迎角的增大而減小。Mackowski和Williamson[18]的研究表明俯仰翼型阻力-推力轉(zhuǎn)變對應(yīng)的臨界縮減頻率隨雷諾數(shù)的增大而減小。Das等[14]進一步指出，俯仰翼型阻力-推力轉(zhuǎn)變對應(yīng)的臨界斯特勞哈爾數(shù)Sttr與雷諾數(shù)之

間存在簡單的比例關(guān)系Sttr～Re-0.37。Deng等[13,19-20]通過對二維基準流動施加展向周期性小擾動，結(jié)合Floquet理論分析了俯仰翼型尾跡流場的穩(wěn)定性，發(fā)現(xiàn)隨著振蕩頻率或幅度增大，主模態(tài)依次經(jīng)歷了由穩(wěn)定模態(tài)、臨界穩(wěn)定模態(tài)到不穩(wěn)定模態(tài)的演變過程，其中不穩(wěn)定模態(tài)反映了流場的固有三維線性不穩(wěn)定性，相應(yīng)的尾跡稱為三維尾跡，主模態(tài)為臨界穩(wěn)定模態(tài)的尾跡為二維尾跡-三維尾跡的轉(zhuǎn)變邊界；文獻[13]進一步開展了三維直接數(shù)值模擬，驗證了其與二維Floquet分析的一致性；文獻[20]討論了尾跡形態(tài)與推進效率的關(guān)系，指出俯仰翼型的推進效率峰值與二維尾跡-三維尾跡轉(zhuǎn)變邊界吻合較好。田偉等[21]的實驗研究發(fā)現(xiàn)，俯仰軸位置對尾跡渦演化及推力特性有較大的影響。戴龍珍和張星[22]對比研究了撲翼在連續(xù)式和間歇式俯仰運動下的能耗問題，指出間歇式俯仰運動在推進速度較高時能耗更低，連續(xù)式俯仰運動則在推進速度較低時更節(jié)省能耗。

盡管上述工作對撲翼的流場特征和推力特性的認識已相對深入，但目前還沒有相關(guān)文獻針對在低雷諾數(shù)下，隨著斯特勞哈爾數(shù)增大，俯仰振蕩翼型的推力產(chǎn)生滯后于反卡門渦街形成的現(xiàn)象進行系統(tǒng)性研究，尤其對其背后的流動機制的認識存在不足。為此，本文采用CFD方法對NACA0012翼型在1 000雷諾數(shù)下作簡諧俯仰運動的非定常流場進行了數(shù)值模擬，考察了阻力-推力轉(zhuǎn)變臨界點附近的尾跡形態(tài)特征和推力特性。采用翼型表面積分方法和基于有限控制體的氣動力估計方法分別研究翼面分布力特性和尾跡流場特性變化對推力產(chǎn)生的影響，以探究阻力-推力轉(zhuǎn)變滯后于卡門渦街-反卡門渦街轉(zhuǎn)變的物理機制。

1 數(shù)值計算和驗證

1.1 控制參數(shù)

考慮NACA0012翼型以正弦規(guī)律繞1/4弦長位置作純俯仰運動：

其中，α0為 俯仰幅度，f為俯仰頻率。利用雷諾數(shù)，無量綱幅度和斯特勞哈爾數(shù)表征繞撲翼流動，其中無量綱幅度和斯特勞哈爾數(shù)分別定義為：

其中，D為翼型厚度，A為在一個俯仰周期內(nèi)翼型后緣通過的垂直距離。

本文選取俯仰幅度（α0）為2°～10°，對應(yīng)的無量綱幅度（AD）為0.44～2.17；對每個α0，通過改變俯仰頻率使斯特勞哈爾數(shù)（St）在0.08～0.49范圍內(nèi)（對α0= 2°，St= 0.08～0.9）。由此形成包含81個算例的數(shù)據(jù)集，內(nèi)容涵蓋了無量綱幅度和斯特勞哈爾數(shù)在不同取值下的推力系數(shù)及其分量，以及相應(yīng)的流場信息。

1.2 求解設(shè)置

計算域為足夠大的圓形區(qū)域，直徑為100倍翼型弦長，圓心位于1/4弦長位置，如圖1所示。采用481（周向）×241（徑向）的O型結(jié)構(gòu)網(wǎng)格對計算域進行離散，第一層網(wǎng)格高度為0.001c，滿足y+＜ 1。對直徑為10倍弦長的內(nèi)部圓形區(qū)域加密（徑向），使其具有較高的流場分辨率；該圓形區(qū)域之外為較粗化的網(wǎng)格。兩塊網(wǎng)格的周向節(jié)點數(shù)一致，且在交界面上的節(jié)點一一對應(yīng)。

圖1 計算模型及網(wǎng)格Fig. 1 Computational model and the corresponding grid used

采用Fluent 17.2壓力基求解器數(shù)值求解二維不可壓N-S方程。遠場設(shè)置為速度入口邊界條件，來流速度為0.146 m/s，相應(yīng)的雷諾數(shù)為1 000。選擇該雷諾數(shù)的原因是，一方面Re=1 000處于MAV的雷諾數(shù)范圍內(nèi)[3]，另一方面已開展的大部分研究中雷諾數(shù)的選擇都集中在103量級[8-9,12-14,16,20]。翼型為無滑移壁面邊界，兩塊網(wǎng)格交界面為內(nèi)部網(wǎng)格面（Interior）。該雷諾數(shù)下流體黏性較強，選用層流模型進行模擬。翼型運動由用戶自定義功能（User Defined Function,UDF）文件指定，整個流體域固連在翼型上，隨翼型作相同的剛性運動。

壓力-速度耦合求解采用SIMPLEC算法，空間離散中對流項和黏性項分別采用二階迎風(fēng)格式和二階中心格式，時間項采用二階隱式雙步法，對每個俯仰周期分N個時間步進行時間推進，每個時間步內(nèi)進行25次內(nèi)迭代。為獲得周期性穩(wěn)定解，每個算例均至少計算60個俯仰周期。

1.3 相關(guān)驗證

采用四套疏密程度不同的網(wǎng)格驗證網(wǎng)格無關(guān)性，網(wǎng)格信息見表1。驗證算例中選取α0= 10°，StA=0.5。如圖2（a）所示，除粗化網(wǎng)格Grid 1外，兩套中等網(wǎng)格Grid 2、Grid 3與細化網(wǎng)格Grid 4的阻力系數(shù)和力矩系數(shù)的時間歷程曲線十分接近，Grid 3和Grid 4的周期平均阻力系數(shù)和力矩系數(shù)的相對誤差低于1%。為驗證時間步長無關(guān)性，選取4個依次加倍的N值，即時間步長依次減半。圖2（b）顯示，N為720和1 440時，計算結(jié)果差異較小。為保證較高的計算精度同時適當(dāng)減小計算量，后續(xù)所有計算中都采用細化網(wǎng)格，時間步長取為T/ 720。

表1 網(wǎng)格信息Table 1 Information of the grids

圖2 網(wǎng)格和時間步長無關(guān)性檢查Fig. 2 Grid and time step independence test results

為進一步驗證求解器的可靠性，將計算結(jié)果與文獻中的結(jié)果進行對比。Deng等[20]的CFD計算中Re=1 700， 選取其中AD=2S t=0.5的算例進行對比驗證，為保持一致性，將模型替換為NACA0015翼型，俯仰軸選在前緣位置。圖3顯示，阻力系數(shù)和力矩系數(shù)的時間歷程曲線與Deng等的結(jié)果十分吻合。

圖3 阻力/力矩系數(shù)與文獻結(jié)果[20]對比Fig. 3 Drag and momentum coefficients compared with the Ref. [20]

2 基于有限控制體的氣動力分解方法

為分析俯仰翼型誘導(dǎo)形成的非定常尾跡對推力產(chǎn)生的影響，采用基于有限控制體的氣動力估計方法建立局部尾跡與推力產(chǎn)生的聯(lián)系。對不可壓黏性流動，由動量定理，作用在翼型上的氣動力為：

其中，p、ρ、τ分別為壓力、流體密度和黏性應(yīng)力矢量，a=du/dt為流體加速度，u為 流體速度。Vf為包圍翼型的有限控制體，其內(nèi)邊界為翼面?B，外邊界為Σn為指向控制體外部的單位法向量。方程（4）是在翼，面直接積分求解氣動力的常用方法，被大部分CFD求解器所采用，根據(jù)該式可以分析壓力和黏性應(yīng)力分布對氣動力產(chǎn)生的貢獻。方程（5）通過包圍翼面的控制體間接求解翼型所受的氣動力。Wang等[23]在方程（5）基礎(chǔ)上利用積分變換推導(dǎo)得到：

其中，ω為渦量。上式中每一項都有明確的物理意義[23-24]：第（1）項為lamb矢量或渦力項；第（2）項和第（3）項分別為由翼型運動的非定常慣性效應(yīng)引起的局部流體加速和翼型占據(jù)的虛擬流體對氣動力的貢獻，在無黏無旋流動中二者之和稱為附加質(zhì)量力；第（4）至第（6）項依次為外控制面上的動量誘導(dǎo)力、壓力誘導(dǎo)力和黏性力。

通過實驗方法得到壓力場的空間分布并不方便，而對速度場的測量相對容易且測量技術(shù)較為成熟[25]。Wang等結(jié)合不可壓N-S方程進一步將壓力場轉(zhuǎn)化為速度場，同時選取特殊的矩形控制體以消除方程（6）中的壓力項。在CFD計算中容易得到壓力場，因此本文直接應(yīng)用方程（6）進行氣動力分解且不對控制體的選取作任何限制。為避免后處理中的數(shù)值插值和積分誤差，將方程（4）和方程（6）寫入UDF文件，基于原始網(wǎng)格在CFD求解計算的過程中同時完成對各氣動力分量的計算。為減小由渦量耗散引入的積分誤差和不確定性，控制體與圖1中計算網(wǎng)格的加密區(qū)一致，該范圍內(nèi)渦量耗散程度較小。

3 結(jié)果和討論

3.1 推力產(chǎn)生滯后于反卡門渦街出現(xiàn)

圖4給出了零推力點附近不同俯仰幅度下周期平均推力系數(shù)隨斯特勞哈爾數(shù)的變化，圖中曲線為對各離散點的二次函數(shù)擬合，零推力點由擬合函數(shù)近似給出，并經(jīng)計算驗證滿足推力系數(shù)在10-4量級。圖4（a）顯示，在各俯仰幅度下，周期平均推力系數(shù)隨St的增大而增大，在St減小至零時趨于某一定值，對應(yīng)靜止翼型所受到的阻力（圖中未標出）；發(fā)生阻力-推力轉(zhuǎn)變的臨界斯特勞哈爾數(shù)StT=0隨俯仰幅度的增大而減小。圖4（b）顯示不同俯仰幅度下推力系數(shù)相對于StA的 變化曲線較為接近，均在StA= 0.3附近出現(xiàn)零推力點?？紤]到StA=ADSt，因此阻力-推力轉(zhuǎn)變臨界斯特勞哈爾數(shù)StT=0與 無量綱幅度AD近似為反比例關(guān)系。

圖4 不同俯仰幅度下的周期平均推力系數(shù)Fig. 4 Cycle averaged thrust coefficients for different pitching amplitudes

圖5對比了俯仰幅度為10°時不同斯特勞哈爾數(shù)下的尾跡形態(tài)。St= 0.066時為典型的卡門渦街，其中順時針旋轉(zhuǎn)的渦位于逆時針旋轉(zhuǎn)的渦上方。隨St增大，順時針渦和逆時針渦沿縱向靠攏，在St=0.082時基本處于一條直線上。此后隨St增大順時針渦與逆時針渦的縱向位置交換，前者位于后者下方，呈反卡門渦街排列。隨St進一步增大，逐漸出現(xiàn)不對稱尾跡和正反渦配對現(xiàn)象（St= 0.185～0.247）。注意到在St= 0.103時反卡門渦街就已經(jīng)出現(xiàn)，但并未產(chǎn)生推力，直到St= 0.143 時才達到零推力點。因此推力產(chǎn)生明顯滯后于反卡門渦街的出現(xiàn)。

圖5 α0= 10°時尾跡形態(tài)隨 St 的變化Fig. 5 Variation of the wake pattern with St at α0= 10°

為較準確判斷卡門渦街-反卡門渦街轉(zhuǎn)變臨界點，根據(jù)渦量極值點確定了后緣下游0.25～3.25倍弦長范圍內(nèi)所有尾跡渦的渦核位置(xc,i,yc,i)，定義渦核縱向平均偏離量為該范圍內(nèi)所有渦核縱向位置絕對值的平均值：

圖6 渦核縱向平均偏離量隨 StA 的變化Fig. 6 Variation of the averaged longitudinal deviation of vortex cores withStA

圖7 不同俯仰幅度下的臨界尾跡Fig. 7 Neutral wakes for various pitching amplitudes

基于以上分析，圖8繪制了尾跡形態(tài)隨無量綱幅度AD和斯特勞哈爾數(shù)St的演變過程。由于StA=ADSt，因此當(dāng)StA取值一定時，AD與St的變化關(guān)系對應(yīng)AD-St參數(shù)空間中的一條反比例函數(shù)曲線。可以發(fā)現(xiàn)，除AD=0.44(α0= 2°）外，其他AD取值下卡門渦街（BvK）-反卡門渦街(RBvK)轉(zhuǎn)變臨界點即臨界尾跡點（Neutral Wake）大致分布在曲線StA=0.18上；類似地，曲線StA=0.37基本對應(yīng)反卡門渦街-不對稱尾跡（Asymmetic Wake）轉(zhuǎn)變臨界。曲線StA=0.18和StA=0.37可以將1 000雷諾數(shù)下NACA0012翼型俯仰運動誘導(dǎo)產(chǎn)生的尾跡形態(tài)劃分為三個區(qū)域：StA=0.18左下方為卡門渦街區(qū)域，StA=0.37右上方為不對稱尾跡區(qū)域，介于二者之間的帶狀區(qū)為反卡門渦街區(qū)域。圖中同時繪制了阻力-推力轉(zhuǎn)變臨界點（=0），大致分布在曲線StA=0.3上 。曲線StA=0.3處于反卡門渦街區(qū)域中間靠后位置，其左下方為阻力區(qū)，右上方為推力區(qū)。以上臨界線位置與Godoy-Diana等(Re=1 173)[12]及Andersen團隊(Re=1 320,2 640)[8-9]的實驗及數(shù)值結(jié)果接近，盡管選擇的翼型形狀不同。

圖8 尾跡形態(tài)在 AD-St 空間中的分布Fig. 8 Distribution of the wake pattern in the parameter space AD-St

圖8顯示，在每個俯仰幅度取值下（α0= 2°～10°或AD=0.44～2.17），卡門渦街-反卡門渦街轉(zhuǎn)變臨界斯特勞哈爾數(shù)St*和阻力-推力轉(zhuǎn)變臨界斯特勞哈爾數(shù)StT=0均不相等，且StT=0大于St*，從而推力產(chǎn)生滯后于反卡門渦街出現(xiàn)。定義無量綱數(shù)ΔSt*=StT=0-St*來表征二者之間的相對滯后程度，ΔSt*越大表明推力產(chǎn)生越滯后于反卡門渦街形成。圖9給出了ΔS t*隨俯仰幅度的變化，可以看出，ΔSt*隨α0增大而減小，從而推力產(chǎn)生與反卡門渦街出現(xiàn)的相對滯后程度隨俯仰幅度的增大而減小。ΔSt*在0.06～0.14之間，對應(yīng)俯仰振蕩頻率范圍為0.7～1.7 Hz。

圖9 ΔSt*與俯仰幅度的關(guān)系Fig. 9 Relationship between Δ St* and the pitching amplitude

3.2 黏性效應(yīng)對阻力-推力轉(zhuǎn)變臨界點的影響

不可壓黏性流動中氣動力的主要來源是翼型表面的壓力和黏性應(yīng)力分布，當(dāng)壓力分布不均產(chǎn)生的推力與剪切層的黏滯阻力相平衡時出現(xiàn)零推力點。為分析造成俯仰翼型推力產(chǎn)生滯后于反卡門渦街出現(xiàn)的主要因素，首先利用方程（4）考察翼型表面非定常壓力和黏性應(yīng)力分布對推力產(chǎn)生的貢獻。圖10給出了俯仰幅度為6°和10°時周期平均推力系數(shù)及其分量隨斯特勞哈爾數(shù)的變化。發(fā)現(xiàn)基于方程（4）計算的總推力系數(shù)隨St的變化曲線與CFD軟件直接給出的結(jié)果幾乎完全一致，表明該氣動力分解中引入的積分誤差很小。

圖10 不同俯仰幅度下的周期平均推力系數(shù)及其分量Fig. 10 Cycle averaged thrust coefficient and the corresponding components for different pitching amplitudes

圖10顯示，在兩個俯仰幅度下，翼面壓力分布產(chǎn)生的推力分量僅在St很小時為負值，為壓差阻力，當(dāng)St大于某一臨界值時均為正值，為壓力誘導(dǎo)推力，定義該臨界值為StTP=0。在整個St取值范圍內(nèi)，翼面黏性應(yīng)力分布產(chǎn)生的推力分量均在零推力線下方，為黏滯阻力，并隨St增大呈緩慢線性增大的趨勢。推力系數(shù)隨St的變化曲線與其壓力分量曲線接近，尤其在更大的俯仰幅度下（α0= 10°），表明該雷諾數(shù)下俯仰翼型的推力主要是由翼面附近的非定常壓力場產(chǎn)生的。

為更直觀地顯示，圖11將圖10中零推力點附近進行了放大，對黏滯阻力加了負號。圖中陰影區(qū)為反卡門渦街對應(yīng)的斯特勞哈爾數(shù)范圍，其左側(cè)為卡門渦街區(qū)域，右側(cè)為不對稱尾跡區(qū)域。俯仰幅度為6°時，推力的壓力分量在卡門渦街-反卡門渦街轉(zhuǎn)變臨界點之前就變?yōu)檎担珶o法克服黏性剪切層產(chǎn)生的阻力；之后在反卡門渦街區(qū)域中間靠左部分，壓力誘導(dǎo)推力隨St增大迅速增大，但仍不足以抵消剪切層的黏滯阻力，從而總的推力依舊為負值，表現(xiàn)為阻力；直到St增大到臨界值StT=0時，由翼面非定常壓力分布貢獻的推力恰好與剪切層的黏滯阻力相等，總的推力為零。類似地，α0= 10°時總推力及其分量的變化特征與α0= 6°一致，此時推力的壓力分量恰好在卡門渦街-反卡門渦街轉(zhuǎn)變臨界點為零，即StTP=0≈St*。

圖8中繪制了推力的壓力分量為零時對應(yīng)的臨界斯特勞哈爾數(shù)StTP=0隨無量綱幅度的變化，即圖中曲線。曲線左下方為壓差阻力區(qū)，右上方為壓力誘導(dǎo)推力區(qū)?？梢钥闯?，各俯仰幅度下，壓力分量臨界點均在卡門渦街-反卡門渦街轉(zhuǎn)變臨界點之前或附近出現(xiàn)，表明推力產(chǎn)生滯后于反卡門渦街形成不是由翼面的壓差阻力造成的。在曲線StA=0.18 和StA=0.3之間的反卡門渦街區(qū)域，壓力誘導(dǎo)推力無法克服剪切層的黏滯阻力，是造成零推力點后移的主要原因。從圖10和圖11中發(fā)現(xiàn)，對給定的St值，壓力誘導(dǎo)推力大小嚴重依賴于俯仰幅度，且在α0增大時大幅增大，而黏滯阻力的大小對α0的變化并不敏感，因此俯仰幅度更大時更容易產(chǎn)生足夠抵消黏滯阻力的壓力誘導(dǎo)推力，從而零推力點更靠近卡門渦街-反卡門渦街轉(zhuǎn)變臨界點。這就解釋了為什么推力產(chǎn)生相對于反卡門渦街出現(xiàn)的滯后程度隨俯仰幅度增大而減小。

圖11 對圖10中零推力點附近區(qū)域放大顯示Fig. 11 Zoomed-in view of the zero thrust region in figure 10

對給定的俯仰幅度，在推力為零和其壓力分量為零時可以分別唯一確定一個臨界斯特勞哈爾數(shù)StT=0和StTP=0，對應(yīng)一個數(shù)對（StT=0，StTP=0）；通過改變俯仰幅度就可以確定一系列數(shù)對。圖12反映了由這些數(shù)對確定的點在StT=0-StTP=0參數(shù)空間中的分布。發(fā)現(xiàn)這些點全部落在了StT=0-StTP=0空間中的一條直線上（實線），擬合的相關(guān)系數(shù)幾乎為1。在雷諾數(shù)為無窮大或忽略黏性效應(yīng)影響時，推力全部由翼面壓力分布產(chǎn)生，因此StT=0=StTP=0，對應(yīng)圖12中的虛線。這表明低雷諾數(shù)的黏性效應(yīng)將無黏時的零推力臨界斯特勞哈爾數(shù)StTP=0按照線性比例規(guī)律放大為StT=0，且當(dāng)俯仰幅度在2°～10°范圍內(nèi)時放大因子（系數(shù)）不依賴于俯仰幅度大小。對低雷諾數(shù)流動，翼型表面的壓力分布可以通過實驗測量，而對黏性應(yīng)力分布的測量比較困難。因此根據(jù)StT=0和StTP=0的線性依賴性可以直接由壓力分量臨界點確定總推力的臨界點，從而忽略對黏滯阻力的測量。StT=0和StTP=0之間的線性相關(guān)規(guī)律僅是針對NACA0012翼型在1 000雷諾數(shù)下作較小幅度（α0= 2°～10°）俯仰運動得到的，對不同雷諾數(shù)、不同翼型形狀及更大的俯仰幅度是否成立還需要進一步探究。

傳感器節(jié)點的總體結(jié)構(gòu)由圖4所示。主要包含感知模塊、處理器模塊、無線WIFI通信模塊、電源模塊。感知模塊由傳感器接口與AD模塊組成,負責(zé)測點數(shù)據(jù)的采集和轉(zhuǎn)換,處理器模塊與無線WIFI通信模塊通過串口連接,建立無線通信橋梁,與主控制器實現(xiàn)控制命令和數(shù)據(jù)的無線傳輸。

圖12 推力與其壓力分量的臨界斯特勞哈爾數(shù)的關(guān)系Fig. 12 Correlation between the critical Strouhal numbers of the thrust and its pressure component

3.3 非定常尾跡對阻力-推力轉(zhuǎn)變臨界點的影響

尾跡形態(tài)隨斯特勞哈爾數(shù)的演變是撲翼流場的主要特征，為分析其對推力產(chǎn)生的影響，基于方程（6）定量估計非定常尾跡對推力的貢獻。關(guān)于該方程中各分量的物理意義及控制體的選取已在第2節(jié)說明。圖13對比了俯仰幅度為10°時由CFD直接計算的各瞬時氣動力系數(shù)和基于有限控制體的尾跡流場積分結(jié)果。發(fā)現(xiàn)對不同的S t值，兩種方法得到的升力系數(shù)和力矩系數(shù)的時間歷程曲線幾乎完全一致，阻力系數(shù)除在峰值位置有微小差異外也吻合得很好。圖14對比了由CFD得到的周期平均推力系數(shù)和尾跡流場積分結(jié)果?？梢园l(fā)現(xiàn)，俯仰幅度為6°時，兩種方法得到的結(jié)果基本一致，尾跡積分結(jié)果略微超估了推力；α0= 10°時兩種方法的計算結(jié)果吻合得很好?；诳刂企w的積分結(jié)果對更大的俯仰幅度來說精確度更高；在俯仰幅度較小時，尾跡渦自身強度較弱，在向下游運動時渦量耗散得較快，從而導(dǎo)致了一定的積分誤差。

圖13 各氣動力系數(shù)的時間歷程曲線，α0=10oFig. 13 Time histories of the aerodynamic force coefficients atα0=10o

圖14 CFD得到的周期平均推力系數(shù)與尾跡積分結(jié)果對比Fig. 14 Comparison of the cycle averaged thrust coefficient computed by CFD with that by the wake integration

圖15對比了俯仰幅度為10°時不同斯特勞哈爾數(shù)下各推力系數(shù)分量的時間歷程曲線，圖中推力系數(shù)上標V表示渦推力，A和B分別為翼型運動的非定常慣性效應(yīng)引起的局部流體加速和翼型占據(jù)的虛擬流體對推力的貢獻，M、P、S分別代表動量誘導(dǎo)推力、壓力誘導(dǎo)推力和黏性剪切力。

圖15 各推力系數(shù)分量的時間歷程曲線，α0=10oFig. 15 Time histories of different components of the thrust coefficient atα0=10o

基于以上分析，非定常尾跡對推力的貢獻主要包括三部分，即渦推力、動量誘導(dǎo)推力和壓力誘導(dǎo)推力。圖16分別給出了相應(yīng)的時均流場，其中速度場為沿來流方向的速度分量；圖中分別對速度場和壓力場以來流速度和局部最大壓力（主流的壓力）為參考進行了無量綱處理。圖16（a）中藍色和紅色區(qū)域分別是順時針渦和逆時針渦的時間平均結(jié)果，與圖5中的瞬時渦量場對應(yīng)。圖16（b）中無量綱速度小于1的區(qū)域為速度損失區(qū)（藍色），大于1的區(qū)域為速度增加或“噴流”區(qū)（紅色）。結(jié)合圖5，發(fā)現(xiàn)卡門渦街導(dǎo)致尾跡流場速度損失(St=0.041,0.066)，反卡門渦街使尾跡流場速度增加或產(chǎn)生“噴流”現(xiàn)象(St= 0.103～0.164)，介于二者之間的臨界狀態(tài)(St=0.082)對尾跡速度場的影響很小。圖16（c）中無量綱壓力小于1的區(qū)域為壓力損失區(qū)（藍色），等于1的區(qū)域為主流區(qū)（紅色）。結(jié)合圖5，發(fā)現(xiàn)卡門渦街和反卡門渦街尾跡流場都出現(xiàn)壓力損失現(xiàn)象。

圖15顯示，St= 0.041～0.164范圍內(nèi)，在一個俯仰周期的大部分時間內(nèi)小于零，為渦致阻力；結(jié)合圖5和圖16（a），發(fā)現(xiàn)在零推力點StT=0=0.143（α0=10°）附近，卡門渦和反卡門渦自身都產(chǎn)生渦致阻力。St=0.082～0.164范圍內(nèi)，均在零推力線上方，為動量誘導(dǎo)推力，其幅度隨St增大逐漸增大并更加遠離零推力線。因此動量誘導(dǎo)推力隨St增大而增大，這可以通過圖16（b）中反卡門渦街的“噴流”強度隨St增大而增強解釋。St=0.041時，均為負值，是由卡門渦街尾跡速度損失產(chǎn)生的動量誘導(dǎo)阻力；St=0.066時，在一個周期內(nèi)分布在零推力線上下兩側(cè)，之后的周期平均結(jié)果顯示其對應(yīng)小量的動量誘導(dǎo)推力。總體上，隨斯特勞哈爾數(shù)的變化趨勢與圖5中尾跡形態(tài)變化和圖16（b）中尾跡速度損失或增加的變化特征一致。在一個俯仰周期的絕大部分時間內(nèi)小于零，為尾跡流場壓力損失引入的附加阻力；的幅度隨St的增大而增大，原因是St更大時，尾跡壓力損失程度更深，因此其導(dǎo)致的附加阻力更大，如圖16（c）。

綜上分析，在零推力點StT=0=0.143 （α0= 10°）附近，卡門渦和反卡門渦自身的渦致阻力及尾跡流場壓力損失引入的附加阻力是俯仰翼型阻力的主要來源，此外卡門渦街導(dǎo)致的尾跡速度損失也貢獻小部分阻力；反卡門渦街的“噴流效應(yīng)”產(chǎn)生的動量誘導(dǎo)推力是推力的主要來源。該結(jié)論是針對俯仰幅度為10°得到的，對其它俯仰幅度也成立。為達到零推力點，渦致阻力、壓力損失的附加阻力和動量誘導(dǎo)推力需達到平衡。

圖17顯示了α0= 6°～10°范圍內(nèi)各周期平均推力系數(shù)分量隨斯特勞哈爾數(shù)的變化，其中為基于尾跡積分得到的總推力系數(shù)；圖中的陰影區(qū)域為反卡門渦街對應(yīng)的St范圍；箭頭標注的數(shù)字為基于尾跡分析預(yù)測的零推力點（左側(cè)）與CFD給出的結(jié)果（右側(cè)）對比，二者存在小量差異，表明尾跡積分對推力臨界點的預(yù)測存在一定誤差。對不同俯仰幅度，在零推力點附近，外控制面上的黏性應(yīng)力分布貢獻的周期平均推力幾乎為零，翼型運動慣性效應(yīng)和翼型占據(jù)的虛擬流體對周期平均推力的貢獻也很小可以忽略不計和，這與前面的分析一致。

圖16 時均流場隨 St 的變化Fig. 16 Variation of the time averaged flow field withSt

圖17 各周期平均推力系數(shù)分量隨 St 的變化Fig. 17 Variation of different components of the cycle averaged thrust coefficient withSt

圖17（c）顯示，對α= 10°，St＜ 0.164時，均小于零并隨St的變化程度很小，為較穩(wěn)定的渦致阻力；在St很小時（0.041）為負值，為動量誘導(dǎo)阻力，之后在St增大時變?yōu)檎挡⒋蠓龃?，為較大的動量誘導(dǎo)推力。St＞ 0.164時，急劇增大并最終變?yōu)闇u推力；動量誘導(dǎo)推力大小基本保持不變。在圖中顯示的整個St范圍內(nèi)，均為負值，且隨St增大而減小，表明由尾跡壓力損失引入的附加阻力隨St的增大而增大。在卡門渦街-反卡門渦街轉(zhuǎn)變臨界點（St*=0.082），即圖中陰影區(qū)左側(cè)邊界位置，小量的動量誘導(dǎo)推力無法克服渦致阻力和尾跡壓力損失引入的附加阻力，總推力為負值，表現(xiàn)為阻力；之后進入反卡門渦街區(qū)域，動量誘導(dǎo)推力隨St增大逐漸增大，但仍不足以抵消渦致阻力和尾跡壓力損失的附加阻力，直到與和相抵消時達到零推力點。因此，反卡門渦街的“噴流效應(yīng)”產(chǎn)生的動量誘導(dǎo)推力無法克服反卡門渦自身的渦致阻力和尾跡流場壓力損失引入的附加阻力，是俯仰翼型推力產(chǎn)生滯后于反卡門渦街出現(xiàn)的主要原因。該結(jié)論對其它俯仰幅度也成立，如圖17（a）和圖17（b）。α0= 6°時尾跡流場壓力損失引入的附加阻力隨St增大先略有減小后急劇增大；α0= 8°時各周期平均推力系數(shù)分量隨St的變化規(guī)律與α0= 10°時基本一致。

4 結(jié) 論

本文針對在低雷諾數(shù)下(Re=1 000)，隨著斯特勞哈爾數(shù)增大，俯仰振蕩翼型推力產(chǎn)生滯后于反卡門渦街出現(xiàn)的現(xiàn)象及其背后的物理機制進行了分析，主要結(jié)論如下：

1） 俯仰幅度在2°～10°范圍內(nèi)時，隨St增大，卡門渦街-反卡門渦街轉(zhuǎn)變臨界點與阻力-推力轉(zhuǎn)變臨界點均不重合，且前者先于后者出現(xiàn)；二者之間的相對滯后程度隨俯仰幅度的增大而減小。

2） 俯仰翼型的推力主要由翼型表面壓力分布產(chǎn)生，黏性應(yīng)力分布對推力大小的影響有限。不同俯仰幅度下，翼面壓力分布積分得到的推力分量在卡門渦街-反卡門渦街轉(zhuǎn)變臨界點之前或附近就由負值變?yōu)檎?，因此零推力點后移不是由翼面的壓差阻力造成的。

3） 當(dāng)翼型振蕩參數(shù)進入對應(yīng)反卡門渦街尾跡形態(tài)區(qū)域時，翼面壓力分布產(chǎn)生的推力分量無法克服剪切層的黏滯阻力，是導(dǎo)致推力產(chǎn)生滯后于反卡門渦街出現(xiàn)的主要原因。

4） 基于控制體的氣動力分析表明，盡管反卡門渦街會產(chǎn)生“噴流效應(yīng)”，但在St較小時，其產(chǎn)生的動量誘導(dǎo)推力無法克服反卡門渦自身的渦致阻力和尾跡流場壓力損失引入的附加阻力，因此即使存在反卡門渦街也不能產(chǎn)生凈推力。

本文系統(tǒng)地討論了俯仰翼型推力產(chǎn)生與反卡門渦街形成的相對滯后特征，較深入地探究了該滯后現(xiàn)象的發(fā)生機制，但也存在其局限性，如暫未考慮來流速度及雷諾數(shù)的影響。此外，翼型形狀、俯仰軸位置的影響還需進一步探究。尾跡形態(tài)演變是撲翼流場的主要特征，深入認知反卡門渦街與推力產(chǎn)生的關(guān)系有利于揭示撲翼推力產(chǎn)生機理，在撲翼式MAV設(shè)計方面具有潛在的應(yīng)用價值。