俯仰振蕩翼型推阻力轉變滯后機制數值研究

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

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

0 引 言

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

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

尾跡形態隨撲翼運動學、動力學的演變是撲翼流場的主要特征，其與推力產生的關系是研究的重點。Von Karman和Burgers[10]在對平板橫向振蕩的研究中發現，反卡門渦街的“噴流”現象能夠產生推力。Koochesfahani[4]的研究表明，當俯仰翼型尾跡形態由卡門渦街轉變為反卡門渦街時，對應的尾跡流場由速度損失型變為速度增加型或“噴流”。他進一步基于積分動量定理證實，反卡門渦街是一種推力產生型尾跡。Jones等[5]在對浮沉翼型的研究中也發現，隨著無量綱浮沉速度的增大，發生阻力型尾跡(卡門渦街)到推力型尾跡(反卡門渦街)的轉變。Lai和Platzer[11]在此基礎上研究了浮沉翼型尾跡的“噴流”特性，指出最大噴流速度是無量綱浮沉速度的線性函數。

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

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

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

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

1 數值計算和驗證

1.1 控制參數

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

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

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

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

1.2 求解設置

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

圖1 計算模型及網格Fig. 1 Computational model and the corresponding grid used

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

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

1.3 相關驗證

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

表1 網格信息Table 1 Information of the grids

圖2 網格和時間步長無關性檢查Fig. 2 Grid and time step independence test results

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

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

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

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

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

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

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

3 結果和討論

3.1 推力產生滯后于反卡門渦街出現

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

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

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

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

為較準確判斷卡門渦街-反卡門渦街轉變臨界點，根據渦量極值點確定了后緣下游0.25～3.25倍弦長范圍內所有尾跡渦的渦核位置(xc,i,yc,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繪制了尾跡形態隨無量綱幅度AD和斯特勞哈爾數St的演變過程。由于StA=ADSt，因此當StA取值一定時，AD與St的變化關系對應AD-St參數空間中的一條反比例函數曲線。可以發現，除AD=0.44(α0= 2°）外，其他AD取值下卡門渦街（BvK）-反卡門渦街(RBvK)轉變臨界點即臨界尾跡點（Neutral Wake）大致分布在曲線StA=0.18上；類似地，曲線StA=0.37基本對應反卡門渦街-不對稱尾跡（Asymmetic Wake）轉變臨界。曲線StA=0.18和StA=0.37可以將1 000雷諾數下NACA0012翼型俯仰運動誘導產生的尾跡形態劃分為三個區域：StA=0.18左下方為卡門渦街區域，StA=0.37右上方為不對稱尾跡區域，介于二者之間的帶狀區為反卡門渦街區域。圖中同時繪制了阻力-推力轉變臨界點（=0），大致分布在曲線StA=0.3上 。曲線StA=0.3處于反卡門渦街區域中間靠后位置，其左下方為阻力區，右上方為推力區。以上臨界線位置與Godoy-Diana等(Re=1 173)[12]及Andersen團隊(Re=1 320,2 640)[8-9]的實驗及數值結果接近，盡管選擇的翼型形狀不同。

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

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

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

3.2 黏性效應對阻力-推力轉變臨界點的影響

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

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

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

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

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

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

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

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

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

3.3 非定常尾跡對阻力-推力轉變臨界點的影響

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

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

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

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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

3） 當翼型振蕩參數進入對應反卡門渦街尾跡形態區域時，翼面壓力分布產生的推力分量無法克服剪切層的黏滯阻力，是導致推力產生滯后于反卡門渦街出現的主要原因。

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

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