翼型動態(tài)失速氣動力二次峰值數值模擬研究

井思夢，趙國慶，招啟軍

（南京航空航天大學直升機旋翼動力學國家級重點實驗室，南京 210016）

動態(tài)失速是一種非定常流動分離現(xiàn)象，通常發(fā)生在大拉力或機動飛行狀態(tài)的直升機旋翼上。動態(tài)失速一旦發(fā)生，會引起升力突降、旋翼反扭矩和槳葉振動突增等問題，從而限制直升機的最大飛行速度和機動性能。因此，動態(tài)失速機理一直是旋翼空氣動力學領域的研究熱點之一。

迄今，許多學者針對翼型開展了動態(tài)失速機理的研究［1?3］。通過試驗測量翼型的氣動載荷以及流場特性，可以直觀地了解翼型動態(tài)失速過程中氣動力的遲滯效應以及動態(tài)失速渦的發(fā)展。McCros?key 等［4］測量了NACA0012 翼型及其前緣修型翼型的動態(tài)失速特性，發(fā)現(xiàn)渦脫落現(xiàn)象是不同類型動態(tài)失速的主要共同特征。Wang 等［5］采用PIV 技術測量了OA209 和SC1095 翼型前緣渦的輸運速度，研究表明前緣渦的輸運速度主要受到翼型振蕩頻率的影響，當振蕩頻率增加時，前緣渦的輸運速度也會提升。Geissler 等［6］開展了OA312 翼型動態(tài)失速特性的試驗和數值模擬研究，發(fā)現(xiàn)渦量的發(fā)展、脫落和積累對翼型的動態(tài)失速特性有重要影響。

由于翼型動態(tài)失速的風洞試驗較為復雜，成本高，且受到測量設備、技術的限制，因此只能開展有限工況下的研究。隨著計算流體力學的發(fā)展，數值模擬逐漸成為翼型動態(tài)失速特性研究的主要途徑［7?10］。

為了研究翼型外形對動態(tài)失速特性的影響規(guī)律，一些學者開展了相應的研究。王清等［11］和鄒錦華等［12］開展了翼型前緣外形對動態(tài)失速特性的影響規(guī)律研究，結果表明：上翼面前緣外形對翼型的動態(tài)失速特性有顯著影響，并且合理的變形能夠抑制翼型的動態(tài)失速。王友進等［13］針對翼型厚度對動態(tài)失速特性的影響開展了數值模擬分析，研究發(fā)現(xiàn)：薄翼型的動態(tài)失速是由前緣分離引起的，而厚翼型的動態(tài)失速則是由后緣分離引起的，并在分離區(qū)域向翼型前緣擴展的過程中發(fā)生失速。Shar?ma 等［14］模擬分析了NACA 系列不同厚度對稱翼型在低雷諾數下的動態(tài)失速特性，結果表明：較薄的NACA0009 翼型的動態(tài)失速是由層流分離泡的破裂引起的，而最厚的NACA0018 翼型的動態(tài)失速則是由湍流邊界層的分離引起的。

此外，氣動參數以及翼型運動參數對翼型的動態(tài)失速特性也有顯著影響。Kim 等［15］和孔衛(wèi)紅等［16］開展了壓縮性效應對翼型動態(tài)失速特性的影響研究，分析表明：小馬赫數下壓縮性效應的影響較小，而在高馬赫數下，壓縮性效應對翼型的動態(tài)失速特性起主導作用。宋辰瑤等［17］、趙國慶等［18］和楊鶴森等［19］對翼型動態(tài)失速特性開展了參數影響研究，揭示了平均迎角、迎角振幅及縮減頻率等振蕩參數以及雷諾數和馬赫數等氣動參數對氣動力遲滯效應以及氣動力峰值的影響規(guī)律。

直升機旋翼翼型工作的縮減頻率一般在0.03～0.15 之間，趙國慶等［18］針對縮減頻率為0.05，0.10 和0.15 等3 個典型工況下的旋翼翼型動態(tài)失速特性開展了數值模擬研究。分析表明：縮減頻率表征了旋翼翼型動態(tài)失速特性的強弱程度，對旋翼翼型動態(tài)失速的遲滯效應、氣動力系數峰值和失速迎角具有重要影響。當縮減頻率為0.05 時，除了動態(tài)失速渦引起的氣動力峰值外，翼型的氣動力系數曲線還出現(xiàn)了二次峰值，然而文中并未開展對氣動力二次峰值的細致研究。Mcalister 等［20］針對翼型動態(tài)失速的試驗研究也表明：在某些動態(tài)失速狀態(tài)下，翼型的氣動力會出現(xiàn)明顯的二次峰值。Choudhry 等［21］對動態(tài)失速升力特征進行了系統(tǒng)的綜述研究，分析了低雷諾數下翼型以恒定速度抬頭至指定迎角過程中的動態(tài)失速現(xiàn)象和升力變化，發(fā)現(xiàn)了由二次渦結構引起的升力二次峰值現(xiàn)象。

目前，在翼型動態(tài)失速機理方面，相關學者已經做了大量的研究工作。然而，針對翼型動態(tài)失速氣動力二次峰值的研究工作相對較少。因此，本文建立了翼型非定常流場數值模擬方法，在此基礎上開展了翼型動態(tài)失速氣動力二次峰值的機理研究，并著重分析了翼型厚度、彎度及彎度位置等參數對氣動力二次峰值的影響，獲得了一些有意義的結論。

1 翼型非定常流場數值模擬方法

1.1 網格生成方法

通過求解泊松方程來生成圍繞翼型的C 型網格，將泊松方程離散化后可以得到

圖1 給出了圍繞SC1095 翼型的C 型網格。從圖1 可以看出，翼型網格具有良好的正交性和貼體性。

圖1 圍繞SC1095 翼型的C 型網格Fig.1 C?type grid around SC1095 airfoil

為了模擬翼型俯仰振蕩時的非定常氣動特性，采用運動嵌套網格技術實現(xiàn)翼型網格的運動以及翼型網格和背景網格之間的信息交換。采用Hole?Map 方法確定背景網格在翼型網格上的洞邊界，采用Inverse?Map 方法進行背景網格人工內邊界的貢獻單元搜索。圖2 給出了運動嵌套網格系統(tǒng)示意圖。

圖2 運動嵌套網格系統(tǒng)Fig.2 Moving?embedded grid system

1.2 流動控制方程及求解方法

采用積分形式的RANS 方程作為翼型流場的控制方程

式中：ρ表示密度；p表示壓強；E和H分別表示單位質量的總能和總焓；Vr和Vt分別表示相對運動速度和由網格運動引起的牽連速度，且Vr=VVt，其中V表示絕對運動速度；ix與iy分別為x與y方向的單位矢量；τij為黏性應力；Θi為描述黏性和熱傳導的作用項。

采用CLORNS 求解器［22］對流場控制方程進行求解?？臻g離散采用Roe 格式，網格面兩側的流場原始變量采用三階MUSCL 格式進行重構，時間推進采用LU?SGS 隱式格式。為了較好地模擬大迎角下翼型的氣流分離和失速特性，采用S?A 湍流模型模擬流場的黏性系數，實現(xiàn)控制方程的封閉。

1.3 方法驗證

采用所建立的CFD 方法對NACA0012 翼型在深度失速狀態(tài)下的氣動特性進行了模擬。迎角變化規(guī)律為α=14.91°+9.88°sin（2kt），來流馬赫數Ma為0.283，縮減頻率k=0.151，雷諾數Re=3.4×106，t為時間。在該算例中，計算域設置如圖3 所示，翼型表面采用無滑移和絕熱條件，背景網格外邊界采用無反射遠場邊界條件。在劃分翼型網格時，設置壁面第一層網格單元高度為4×10-6c，保證y+≈1。圖4 給出了3 套不同網格的氣動力系數計算值與試驗值［20］的對比，3 套網格的信息在表1中給出。從圖4 中可以看出：計算值與試驗值吻合較好，表明本文的數值模擬方法能夠準確地模擬翼型的動態(tài)失速特性；并且，在翼型上仰過程中，3 套網格的計算結果一致，在下俯過程中，中等網格和密網格的計算結果基本重合。因此，在后文的研究中，均采用了與此算例中等網格相同的網格量。

表1 NACA0012 翼型網格信息Table 1 Information about the three sets of grids

圖3 計算域示意圖Fig.3 Schematic diagram of computational domain

圖4 翼型氣動力系數計算值與試驗值[20]對比Fig.4 Comparison of aerodynamic force coefficients be?tween calculated results and test data[20]

2 動態(tài)失速氣動力二次峰值分析

為了探究翼型動態(tài)失速狀態(tài)下氣動力二次峰值的發(fā)生機理，對NACA0012 翼型在狀態(tài)Ma=0.293，α=14.82°+9.90°sin（2kt），k=0.051，Re=3.82×106下的動態(tài)失速特性進行了模擬分析。

圖5 給出了NACA0012 翼型的升力、阻力和力矩系數曲線，圖中A～F為6 個不同時刻。從圖5中可以看出，NACA0012 翼型在此狀態(tài)下發(fā)生了深度失速，并且失速后升力、阻力和力矩系數出現(xiàn)了二次峰值。

圖5 NACA0012 翼型氣動力系數Fig.5 Aerodynamic force coefficients of NACA0012 airfoil

圖6 給出了圖5 中6 個不同時刻（A～F）NA?CA0012 翼型的壓強系數分布云圖和流線圖。不同時刻對應的翼型迎角如表2 所示。不同時刻翼型表面的壓強系數分布如圖7 所示。

圖6 NACA0012 翼型壓強分布云圖和流線圖Fig.6 Pressure coefficient cloud and freestream lines of NACA0012 airfoil

圖7 NACA0012 翼型壓強系數分布Fig.7 Pressure coefficient distribution of NACA0012 airfoil

表2 不同時刻的翼型迎角Table 2 Angles of attack of airfoil at different time

從圖6 可以看出，在A時刻，翼型上表面已經形成了動態(tài)失速渦（Dynamic stall vortex，DSV），而且翼型上表面中部出現(xiàn)了由動態(tài)失速渦導致的局部負壓區(qū)域。動態(tài)失速渦在翼型前緣附近形成后，會沿著翼型上表面向后緣移動。一方面，當動態(tài)失速渦附著于翼型表面時，渦致升力的存在會使得翼型升力系數持續(xù)增加；另一方面，負壓中心也在向后緣移動，導致低頭力矩逐漸增加。接著，動態(tài)失速渦會從翼型表面脫落，渦致升力開始減小，同時，前緣負壓峰值減小，二者共同作用導致升力系數下降，阻力、低頭力矩進一步增加。在A時刻，動態(tài)失速渦即將從翼型表面脫落，此時翼型的升力系數達到峰值，阻力和低頭力矩系數接近峰值，如圖5 所示。

在B時刻，動態(tài)失速渦誘導形成了后緣渦（Trailing?edge vortex，TEV），導致后緣附近局部負壓峰值的出現(xiàn)，如圖7（a）所示。在主流和后緣渦的作用下，動態(tài)失速渦向上抬升，脫離翼型表面，使得由動態(tài)失速渦引起的局部負壓減小，進而導致升力系數的突降。從圖7（a）可以看出，與A時刻相比，B時刻翼型前緣負壓峰值減小，同時翼型中部動態(tài)失速渦引起的局部負壓減小，并且后緣渦引起了新的局部負壓，使得翼型的阻力系數和低頭力矩系數減小。

在C時刻，動態(tài)失速渦的強度進一步減弱，在后緣渦的作用下，動態(tài)失速渦又逐漸向翼型上表面靠攏。后緣渦上洗并膨脹，強度也逐漸減弱。此時，翼型的升、阻力系數和低頭力矩系數進一步減小。

在D時刻，后緣渦上洗，從翼型后緣脫落并向下游移動至尾跡區(qū)，由后緣渦引起的局部負壓消失。由于后緣渦的誘導作用，動態(tài)失速渦繼續(xù)向翼型上表面靠攏，并重新附著于翼型上表面。此時，動態(tài)失速渦較弱，在翼型表面并沒有形成明顯的局部負壓。但是，如圖7（b）所示，與C時刻相比，D時刻的前緣負壓峰值和翼型上表面負壓明顯增加，導致升力、阻力和力矩系數回升。

由圖7（b）可知，在E時刻，翼型前緣負壓峰值和上表面負壓較D時刻進一步增加。此時，翼型升、阻力系數和低頭力矩系數達到二次峰值。同時，從圖6（e）可以看出，由于動態(tài)失速渦的誘導作用，翼型后緣附近形成了新的較弱的后緣渦。

在F時刻，動態(tài)失速渦向上抬升，遠離翼型表面，分離區(qū)域擴大，氣流分離加劇，導致前緣負壓峰值和翼型上表面負壓的減小，從而導致升力、阻力和力矩系數的下降。

總體而言，動態(tài)失速渦的形成，導致第一個氣動力峰值的出現(xiàn)；動態(tài)失速渦誘導形成后緣渦，在后緣渦的作用下，動態(tài)失速渦向翼型表面靠攏，并重新附著于翼型上表面，導致翼型氣動力二次峰值的出現(xiàn)。

3 翼型外形對氣動力二次峰值的影響分析

為了探究翼型外形對氣動力二次峰值的影響，對NACA 系列不同厚度、彎度和彎度位置的翼型開展了氣動特性模擬分析。計算狀態(tài)同第2 節(jié)。

3.1 翼型厚度影響

為了探究翼型厚度對氣動力二次峰值的影響，對NACA0009、NACA0012 和NACA0015 這3 個不同厚度翼型的動態(tài)失速特性進行了模擬分析。翼型外形對比如圖8 所示。

圖8 不同厚度翼型外形對比Fig.8 Comparison of airfoils with different thicknesses

圖9 給出了NACA0009、NACA0012 和NA?CA0015 這3 個不同厚度翼型的氣動力系數對比。從圖中可以看出，翼型厚度對動態(tài)失速特性有顯著的影響。翼型厚度增加，動態(tài)失速迎角增加，升力、阻力和低頭力矩系數峰值均增加。在線性升力段，相同迎角下厚度較大的翼型產生的升力較小，且厚度對線性段升力線斜率的影響較小。此外，3 個不同厚度翼型在此動態(tài)失速狀態(tài)下均出現(xiàn)了氣動力二次峰值。本文稱兩個峰值之間的氣動力最小值為谷值，由圖9 可知，翼型厚度越大，二次峰值出現(xiàn)的越晚，且二次峰值相對谷值的增量（下文簡稱“二次峰值增量”）越大。

圖9 不同厚度翼型的氣動力系數對比Fig.9 Comparison of aerodynamic force coefficients of airfoils with different thicknesses

圖10 給出了3 個翼型的壓強系數分布云圖和流線圖。從圖10 中可以看出，在翼型上仰過程中迎角為18.00°的狀態(tài)下，NACA0009 翼型上表面已經形成了動態(tài)失速渦，且動態(tài)失速渦移動至后緣附近，誘導形成了后緣渦。此時，NACA0012 翼型前緣形成了動態(tài)失速渦，而NACA0015 翼型的流動仍附著在翼型表面，沒有發(fā)生明顯的氣流分離。當迎角增加至20.23°，NACA0012 翼型上表面的動態(tài)失速渦移動至后緣附近，誘導形成了后緣渦。此時，NACA0015 翼型后緣發(fā)生了氣流分離。當迎角繼續(xù)增加至21.78°，NACA0015 翼型上表面形成了動態(tài)失速渦，NACA0012 翼型的后緣渦脫落至尾跡區(qū)。

圖10 不同厚度翼型的壓強分布云圖和流線圖對比Fig.10 Comparison of pressure coefficient cloud and streamlines of airfoils with different thicknesses

NACA0009、NACA0012 和NACA0015 翼型上表面先后形成了不同強度的動態(tài)失速渦和后緣渦，動態(tài)失速渦和后緣渦相互作用，導致翼型氣動力出現(xiàn)二次峰值。翼型厚度增加，推遲了動態(tài)失速渦和后緣渦的形成，從而推遲了氣動力二次峰值的出現(xiàn)。此外，3 個翼型在此狀態(tài)下均發(fā)生了前緣失速，翼型厚度增加，形成的動態(tài)失速渦和后緣渦強度增加，導致氣動力二次峰值增量增加。

3.2 翼型彎度影響

為了探究翼型彎度對翼型動態(tài)失速狀態(tài)氣動力二次峰值的影響，對NACA0012、NACA2212 和NACA4212 這3 個不同彎度翼型的動態(tài)失速特性進行了模擬分析。翼型外形對比如圖11 所示。

圖11 不同彎度翼型外形對比Fig.11 Comparison of airfoils with different cambers

圖12 給出了NACA0012、NACA2212 和NA?CA4212 這3 個不同彎度翼型的氣動力系數對比。從圖12 中可以看出，彎度對翼型動態(tài)失速特性有明顯影響。迎角相同時，彎度較大的翼型產生的升力系數越大，低頭力矩系數也越大。隨著翼型彎度的增加，動態(tài)失速迎角、阻力和力矩發(fā)散迎角增加，氣動力系數峰值呈現(xiàn)先增加后減小的趨勢。NA?CA0012 和NACA2212 翼型的氣動力系數出現(xiàn)了明顯的二次峰值，而NACA4212 翼型的二次峰值則較小。而且，翼型彎度越大，翼型升力系數二次峰值出現(xiàn)的越晚，二次峰值增量呈現(xiàn)先增加后減小的趨勢，并且與NACA0012 和NACA2212 翼型相比，NACA4212 翼型的升、阻力和力矩系數二次峰值增量均最小。

圖12 不同彎度翼型的氣動力系數對比Fig.12 Comparison of aerodynamic force coefficients of airfoils with different cambers

圖13 給出了3 個翼型在不同時刻的壓強系數分布云圖和流線圖。從圖13 中可以看出，NA?CA0012 和NACA2212 翼型先后發(fā)生了前緣失速，即動態(tài)失速渦在翼型前緣附近形成。而NA?CA4212 翼型則先出現(xiàn)了后緣氣流分離，隨著迎角的增加，分離點向上游移動，并在翼型中后部形成了動態(tài)失速渦。由于動態(tài)失速渦的誘導作用，3 個翼型后緣附近形成了不同強度的后緣渦。當翼型發(fā)生前緣動態(tài)失速時，翼型彎度的增加會推遲動態(tài)失速渦和后緣渦的形成，從而推遲二次峰值的出現(xiàn)；且彎度增加會增強動態(tài)失速渦和后緣渦的強度，從而導致二次峰值增量增加。而當彎度增加到一定程度，翼型上表面后緣附近的曲率變化較大，導致翼型失速類型轉變?yōu)楹缶壥伲c前緣失速情況相比，后緣失速情況下形成的動態(tài)失速渦和后緣渦的強度較小，因此氣動力系數二次峰值增量較小。

圖13 不同彎度翼型的壓強系數分布云圖和流線圖Fig.13 Comparison of pressure coefficient cloud and streamlines of airfoils with different cambers

3.3 翼型彎度位置影響

為了探究翼型彎度位置對動態(tài)失速氣動力二次峰值的影響，對NACA2212、NACA2412 和NACA2612 這3 個不同彎度位置翼型的動態(tài)失速特性進行了模擬分析。翼型外形對比如圖14所示。

圖14 不同彎度位置翼型外形對比Fig.14 Comparison of airfoils with different positions of camber

圖15 給出了3 個不同彎度位置翼型的氣動力系數對比。在線性升力段相同迎角下，彎度位置靠后的翼型產生的升力系數較大，低頭力矩系數也較大。隨著彎度位置的增加，翼型的動態(tài)失速迎角、阻力和力矩發(fā)散迎角減小，氣動力系數峰值也有所減小。彎度位置越靠后，二次峰值出現(xiàn)得越早，且二次峰值增量越小。

圖15 不同彎度位置翼型的氣動力系數對比Fig.15 Comparison of aerodynamic force coefficients of airfoils with different positions of camber

圖16 給出了3 個翼型的壓強系數分布云圖和流線圖。從圖16 中可以看出，3 個翼型都發(fā)生了前緣失速，且彎度位置越靠后，動態(tài)失速發(fā)生的越早，后緣渦的形成也越早，導致二次峰值出現(xiàn)的越早。而且，彎度位置越靠后，動態(tài)失速渦和后緣渦的強度越小，從而導致翼型氣動力系數二次峰值增量越小。

圖16 不同彎度位置翼型的壓強系數分布云圖和流線圖Fig.16 Comparison of pressure coefficient cloud and streamlines of airfoils with different positions of camber

4 結論

基于URANS 方程，結合運動嵌套網格技術、LU?SGS 隱式格式、Roe?MUSCL 數值格式和S?A湍流模型，建立了翼型非定常氣動特性數值模擬方法，開展了翼型動態(tài)失速氣動力二次峰值機理研究和參數影響分析。綜合計算分析結果，得出結論如下：

（1）在動態(tài)失速過程中，翼型上表面會形成動態(tài)失速渦，導致第一個氣動力峰值的出現(xiàn)；動態(tài)失速渦誘導形成后緣渦，在后緣渦的作用下，使得動態(tài)失速渦向翼型表面靠攏，并重新附著于翼型上表面，從而導致翼型氣動力二次峰值的出現(xiàn)。

（2）翼型厚度增加、彎度位置前移以及一定范圍內彎度的增加，會推遲動態(tài)失速渦和后緣渦的形成，并增加動態(tài)失速渦和后緣渦強度，從而使得氣動力二次峰值推遲，二次峰值增量增加。而當翼型彎度較大時，翼型的失速方式從前緣失速變?yōu)楹缶壥佟Ｅc前緣失速情況相比，后緣失速情況下形成的動態(tài)失速渦和后緣渦的強度較小，因此氣動力系數二次峰值增量較小。