趙國慶，招啟軍，王清

旋翼翼型非定常動態(tài)失速特性的CFD模擬及參數(shù)分析

趙國慶，招啟軍，王清

(南京航空航天大學(xué)直升機旋翼動力學(xué)國家級重點實驗室，江蘇南京210016)

構(gòu)建了一套基于運動嵌套網(wǎng)格技術(shù)和可壓縮RANS方程的旋翼翼型非定常流動特性模擬的高效、高精度的CFD方法。首先，發(fā)展了基于Poisson方程求解的圍繞翼型的粘性貼體正交網(wǎng)格生成方法，并提出了基于最小距離法(MDM)改進策略的運動嵌套網(wǎng)格生成方法，克服了彈簧法可能導(dǎo)致網(wǎng)格畸變的不足;其次，為準(zhǔn)確模擬由湍流分離和氣流再附引起的氣動力的遲滯效應(yīng)，基于RANS方程、雙時間方法和高階插值格式，建立了旋翼翼型非定常氣動特性分析的高精度數(shù)值方法，并采用能夠較好捕捉氣流分離現(xiàn)象的S-A湍流模型;再次，針對旋翼后行槳葉動態(tài)失速時槳葉剖面來流速度較低、迎角較大的特點，為解決低來流速度時L-B半經(jīng)驗?zāi)Ｐ驮谛硪硇头嵌ǔ討B(tài)失速計算中的局限性，并克服可壓縮方程對低速流場計算收斂困難和精度低的問題，建立了基于Pletcher-Chen低速預(yù)處理方法、FAS多重網(wǎng)格法和隱式LU-SGS方法相結(jié)合的高效數(shù)值方法。應(yīng)用發(fā)展的方法，分別針對NACA0012、SC1095旋翼翼型靜態(tài)和輕度、深度動態(tài)失速進行計算，精確捕捉了氣動力遲滯效應(yīng)以及翼型前緣脫體渦的產(chǎn)生、對流和脫落過程，驗證了本文方法的有效性;最后，著重針對NACA0012動態(tài)失速狀態(tài)，開展了振蕩參數(shù)對旋翼翼型非定常動態(tài)失速特性影響的分析，研究結(jié)果表明翼型迎角平均值、振幅及減縮頻率的變化均能引起遲滯效應(yīng)的改變并使得氣動力峰值發(fā)生有規(guī)律的前、后移現(xiàn)象等。

旋翼;翼型;動態(tài)失速;N-S方程;運動嵌套網(wǎng)格;參數(shù)分析

0 引言

旋翼工作在嚴(yán)重非對稱、非定常的渦流場中，旋翼槳葉的揮舞、周期變距以及畸變尾跡(誘導(dǎo))形成的非均勻入流，導(dǎo)致槳葉剖面在不同方位角處的迎角有很大差別。當(dāng)旋翼槳盤載荷很高時，后行槳葉工作在較大的迎角狀態(tài)，容易產(chǎn)生氣流分離進而出現(xiàn)復(fù)雜的非定常動態(tài)失速現(xiàn)象。旋翼動態(tài)失速現(xiàn)象雖然能夠增大升力峰值，但同時造成阻力和力矩的突增，且翼型氣動力中心不再穩(wěn)定，這對旋翼的振動特性有重要影響，從而嚴(yán)重制約著直升機氣動性能和飛行速度的提高［1］，因此其準(zhǔn)確預(yù)測有助于直升機飛行性能包線的劃定、振動水平的分析和設(shè)計等。與此同時，翼型作為旋翼的基本組成元素，其非定常動態(tài)失速特性反映了旋翼動態(tài)失速的基本特性。然而，旋翼及翼型的動態(tài)失速現(xiàn)象的機理仍在探索之中，對脫體渦的運動引起的氣動力遲滯效應(yīng)、逆壓梯度的產(chǎn)生及大小等規(guī)律都未能徹底掌握，關(guān)于旋翼翼型非定常動態(tài)失速特性的研究一直是直升機空氣動力學(xué)領(lǐng)域的研究焦點和難點。因此，進行旋翼翼型動態(tài)失速的研究具有重要的理論和實際應(yīng)用價值。

目前，旋翼翼型動態(tài)失速特性的研究主要采用數(shù)值方法和試驗方法。試驗研究主要是針對風(fēng)洞中振蕩的翼型進行氣動力、壓力的測量［2-3］，因而可以直觀地認(rèn)知旋翼翼型動態(tài)失速中氣動力的遲滯效應(yīng)。但是旋翼翼型動態(tài)失速的風(fēng)洞試驗復(fù)雜、代價大、周期長并受試驗測量技術(shù)的限制，并且很難捕捉動態(tài)失速過程中渦的流動細(xì)節(jié)、氣流分離及再附等現(xiàn)象，因而其只能進行原理性或特定工況的研究。數(shù)值計算方法可以克服試驗研究中存在的這些不足和難題，針對旋翼翼型動態(tài)失速可進行多種工況及流動細(xì)節(jié)的模擬，已成為旋翼翼型非定常動態(tài)失速特性研究的一條重要途徑。

旋翼翼型動態(tài)失速模擬的數(shù)值方法包括基于試驗數(shù)據(jù)的半經(jīng)驗?zāi)Ｐ秃陀嬎懔黧w力學(xué)(CFD)方法。國外學(xué)者已發(fā)展了一系列的半經(jīng)驗動態(tài)失速理論模型，如Leishman-Beddoes(L-B)［4］、ONERA［5］、Johnson［6］等模型。其中，最為著名的是L-B二維翼型動態(tài)失速模型，該模型是由幾個從翼型試驗數(shù)據(jù)中獲得經(jīng)驗參數(shù)，采用指數(shù)法建立的半經(jīng)驗?zāi)Ｐ?。L-B模型在達到一定精度的要求下，大大縮短了計算時間，在直升機旋翼氣動載荷、氣動彈性以及飛行力學(xué)等方面有著廣泛的應(yīng)用［7-9］。但L-B等半經(jīng)驗?zāi)Ｐ筒⒉皇菄?yán)格的預(yù)測方法［10］，僅對特定翼型適用，同時L-B模型對動態(tài)失速中氣流再附時氣動力的計算值與試驗值相比有很大偏差，并且針對如旋翼后行槳葉的低來流速度、大迎角狀態(tài)，其無法完全勝任翼型非定常氣動力的準(zhǔn)確模擬;此外，L-B模型只能針對氣動力進行有效計算，無法模擬翼型非定常流場中復(fù)雜的氣流分離及再附現(xiàn)象中的渦流動細(xì)節(jié)，從而很難更深入地探索動態(tài)失速的形成機理。

Leishman指出，只有通過基于Navier-Stokes方程的CFD方法才能完整地描述旋翼翼型的動態(tài)失速現(xiàn)象［11］。Jameson［12］較早地提出了非定常流動模擬的雙時間方法，并采用Eluer方程對翼型的動態(tài)振蕩過程進行數(shù)值模擬，有效地計算了翼型淺度失速時升力系數(shù)的遲滯效應(yīng)，顯示了CFD方法在模擬翼型非定常動態(tài)失速方面的優(yōu)越性。由此，國內(nèi)外學(xué)者相繼開展了翼型動態(tài)失速現(xiàn)象的CFD數(shù)值分析研究［13-17］?，F(xiàn)有CFD研究結(jié)果表明，動態(tài)失速的顯著特征是具有能量的翼型前緣干擾渦的運動，渦運動引起翼型壓力場、力和力矩等的瞬時變化，導(dǎo)致氣動力表現(xiàn)出明顯的非線性遲滯效應(yīng)。然而，這些研究大多針對固定翼飛機的跨音速、小迎角的工作狀態(tài)，對于旋翼后行槳葉翼型所處的低馬赫數(shù)、大迎角狀態(tài)的數(shù)值研究不多，并且針對翼型振蕩參數(shù)對失速特性的影響分析也很少涉及，而上述這些對旋翼后行槳葉動態(tài)失速特性的研究卻至關(guān)重要。

雖然通過CFD方法已經(jīng)能夠有效對旋翼翼型動態(tài)失速中氣動力特性進行預(yù)測，但旋翼翼型非定常動態(tài)失速特性分析的CFD方法的發(fā)展及應(yīng)用還遠(yuǎn)未成熟，對動態(tài)失速現(xiàn)象的模擬結(jié)果與試驗值仍存在較大差異［12］。這主要是因為:(1)旋翼翼型非定常渦流場中存在較大的渦量非物理耗散，主要來源于數(shù)值方法的截斷誤差;(2)目前針對翼型動態(tài)失速的研究多采用變形網(wǎng)格方法展開，雖然通過彈簧方法等網(wǎng)格變形方法能夠有效避免畸形網(wǎng)格的出現(xiàn)，但仍難保證網(wǎng)格正交性和光滑性，另外，變形網(wǎng)格技術(shù)需要滿足幾何守恒定律，這均影響計算精度。

針對這些問題，本文采用基于改進最小距離法和Hole Map方法生成繞旋翼翼型的運動嵌套網(wǎng)格;將Roe-MUSCL格式、隱式LU-SGS方法相結(jié)合，以基于S-A湍流模型的RANS方程模擬大范圍氣流分離現(xiàn)象以及翼型前緣脫體渦的流動過程;針對旋翼后行槳葉來流速度低、L-B模型在低速狀態(tài)的計算失真等問題，并擴大可壓縮N-S方程的應(yīng)用范圍，建立了FAS聚合式多重網(wǎng)格和Pletcher-Chen預(yù)處理方法結(jié)合的高效CFD方法。最后著重開展了旋翼翼型動態(tài)失速特性及深度失速狀態(tài)下振蕩參數(shù)影響分析，獲得了一些有意義的結(jié)論。

圖1 繞翼型的C型網(wǎng)格Fig.1C-type grid around airfoil

1 網(wǎng)格生成方法

1.1旋翼翼型貼體網(wǎng)格生成

旋翼翼型貼體網(wǎng)格是進行動態(tài)失速分析的基礎(chǔ)，網(wǎng)格的正交性、光滑性以及翼型表面附近的網(wǎng)格布置直接影響數(shù)值模擬的精度。為生成壁面網(wǎng)格尺度、正交性滿足湍流求解要求的光滑過渡的翼型網(wǎng)格，本文發(fā)展了基于Poisson方程求解的翼型網(wǎng)格生成方法。計算平面的Poisson控制方程為:

式中，α、β、γ為坐標(biāo)變換系數(shù)，ψp、ψq為控制源項，其分別控制翼型網(wǎng)格的正交性和法向網(wǎng)格密度。

圖1為生成的NACA0012翼型的整體C型網(wǎng)格示意圖，其中翼型表面布置239個網(wǎng)格點，翼型法向60層網(wǎng)格，并設(shè)置第一層網(wǎng)格距離壁面1×10-6c(c為翼型弦長)。由圖可以看出，網(wǎng)格正交性和貼體性良好，從而為N-S方程準(zhǔn)確捕捉粘性及動態(tài)失速過程中氣流分離提供基礎(chǔ)。

圖2 基于改進的MDM方法的貢獻單元搜索過程Fig.2Procedure of donor cell searching by improved MDM method

1.2 運動嵌套網(wǎng)格生成方法

目前國內(nèi)外針對翼型動態(tài)失速的模擬多采用變形網(wǎng)格方法，并通過幾何守恒律計入網(wǎng)格變形對流動通量的影響。但當(dāng)翼型運動幅度較大時，變形網(wǎng)格方法仍可能導(dǎo)致網(wǎng)格的畸變，從而影響計算精度。

基于此，本文采用運動嵌套網(wǎng)格方法進行旋翼翼型動態(tài)失速的非定常流場的數(shù)值模擬。首先，分別生成隨翼型運動的貼體網(wǎng)格和固定不動笛卡爾背景網(wǎng)格;然后，采用Hole Map方法確定背景網(wǎng)格在翼型網(wǎng)格上的洞邊界，在此基礎(chǔ)上，提出了最小距離法(Minimum Distance Method，MDM)的改進策略，并進行背景網(wǎng)格人工內(nèi)邊界的貢獻單元搜索，從而生成翼型與背景網(wǎng)格的運動嵌套網(wǎng)格。

MDM方法有一個不足，即搜索得到的與背景網(wǎng)格點P最近的翼型網(wǎng)格點可能是局部最近點。針對這一問題，文獻［18］在翼型上下表面各選取一個網(wǎng)格點作為搜索起點，但這同時增加了貢獻單元搜索的復(fù)雜性。本文提出了一種改進的MDM方法，采用該方法對背景網(wǎng)格上的點P進行貢獻單元的搜索過程如圖2所示，具體過程如下:

(1)初始化所有翼型網(wǎng)格點標(biāo)記為0，并選取翼型網(wǎng)格點S(II，JJ)為起始點，其中II和JJ分別為沿翼型弦向和法向的網(wǎng)格點編號;

(2)以S點為中心，判斷點S的四個相鄰網(wǎng)格點到P的距離與|SP|的大小關(guān)系。跳過標(biāo)記為1的點，若存在標(biāo)記為0的網(wǎng)格點距離P點最近且小于|SP|，則令該點為新的S點，并將判斷過的網(wǎng)格點標(biāo)記為1，否則，將S點記為O，轉(zhuǎn)到第(4)步;

(3)重復(fù)步驟(2)，直到找到O點;

(4)找出與O點關(guān)于尾跡奇異面對稱的網(wǎng)格點，若|OP|小于|P|，則進入第(5)步;反之則將點記為S，釋放O點，并返回第(2)步;

(5)通過矢量法判斷點P在哪一包含O點的網(wǎng)格單元中，若不在這些單元中，則拓展網(wǎng)格范圍判斷，直到找到點P的貢獻單元。

MDM方法搜索點是在一條曲線上進行的，類似于一維搜索，搜索次數(shù)為O［4×(IMAX+JMAX)］，可以有效地節(jié)省貢獻單元搜索的時間。

翼型網(wǎng)格邊界單元的貢獻單元由兩個一維搜索得到。嵌套網(wǎng)格之間的信息傳遞通過雙線性插值完成。為減小插值過程中引入的數(shù)值誤差，在挖洞時將背景網(wǎng)格洞邊界遠(yuǎn)離流場梯度較大的翼型表面。圖3給出了翼型網(wǎng)格運動、網(wǎng)格嵌套關(guān)系及貢獻單元搜索結(jié)果示意圖。

圖3 圍繞振蕩翼型的運動嵌套網(wǎng)格邊界單元及貢獻單元示意圖Fig.3Schematic diagram of boundary and donor cells of moving-embedded grids around oscillated airfoil

2 旋翼翼型非定常流場數(shù)值計算方法

2.1 流場控制方程求解

綜合考慮旋翼翼型動態(tài)失速的非定常流場的復(fù)雜性，在進行動態(tài)翼型氣動特性分析時采用雷諾平均N-S方程?；诿鎲卧獮閐S的控制體Ω，積分形式的控制方程為:

其中，W、Fv分別為守恒變量和粘性通量。

基于翼型網(wǎng)格的運動，修正無粘通量為:

式中，Vt=nx為控制體表面的逆變速度，nx、ny表示控制面單位外法矢分量，則對流通量為:

其中u、v分別為x和y方向速度。

本文采用格心有限體積法在空間方向上進行離散。無粘通量采用Roe-MUSCL格式進行計算:

其中，A-Roei+1/2為Roe格式算子，i+1/2為左右單元L、R的交接面，WL、WR分別由MUSCL格式插值獲得:

式中，Δ-、Δ+分別為向后、向前差分算子。κ為自由參數(shù)，本文取κ=1/3，即三階迎風(fēng)格式。

為模擬旋翼翼型流場的非定常特性，采用雙時間方法進行時間推進。時間導(dǎo)數(shù)采用向后二階差分進行離散，方程(2)可隱式離散為:

為準(zhǔn)確模擬旋翼翼型動態(tài)失速流場的非定常特性，要求在進行時間推進時取很小的時間步長，這會降低數(shù)值模擬的效率。針對這一問題，本文采用隱式LU-SGS格式進行時間推進，從而有效加大時間步長，提高計算效率。時間推進過程為:

其中，L、D、U分別為下三角、對角和上三角矩陣。

2.2 S-A湍流模型

本文N-S方程湍流粘性系數(shù)計算采用S-A湍流模型。相對于代數(shù)湍流模型(B-L模型)，S-A模型能夠更好地捕捉翼型振蕩引起的大范圍氣流分離現(xiàn)象以及翼型前緣脫體渦運動過程，因而更適合于旋翼翼型動態(tài)失速過程的分析。S-A模型的控制方程為準(zhǔn)運動渦粘性系數(shù)珓ν的輸運方程，其積分形式為:

式中，F(xiàn)c、Fν、QT分別為粘性系數(shù)輸運方程的對流通量、粘性通量和單元源項，表達式如下:

對流項和粘性通量分別采用逆風(fēng)、中心差分格式離散，時間推進方法與N-S方程一致。由珓ν最終求得運動渦粘性系數(shù)νt=珓νfv1，fv1為珓ν的函數(shù)。

2.3 多重網(wǎng)格方法

多重網(wǎng)格方法(Multigrid method)［19］能夠?qū)⒓?xì)網(wǎng)格上的低頻誤差轉(zhuǎn)化為粗網(wǎng)格上的高頻誤差，從而達到在各自網(wǎng)格上消除其高頻誤差的目的，其工作量小、計算效率高，尤其適合旋翼及翼型復(fù)雜非定常流場的加速計算。本文采用FAS格式兩重網(wǎng)格方法，首先，將細(xì)網(wǎng)格上求解得到流場守恒變量W+h和殘值傳遞到粗網(wǎng)格上:

將修正量δW2h=Wn2h+1－W(2h0)疊加到細(xì)網(wǎng)格守恒變量中，得到新的守恒變量值:

為簡化計算，本文在計算時只在網(wǎng)格數(shù)量大的翼型貼體網(wǎng)格上采用多重網(wǎng)格方法。

2.4 低速預(yù)處理方法

旋翼后行槳葉一般工作在低速狀態(tài)，本文采用可壓縮方程進行低來流速度狀態(tài)翼型非定常動態(tài)失速計算時，會遇到“剛性”問題。這是因為當(dāng)?shù)亓鲃铀俣群鸵羲傧嗖钶^大時，時間相關(guān)方程對流項特征矩陣的最大和最小特征值的比值也會很大，使計算收斂很慢。低速預(yù)處理方法能夠通過改變方程的特征值來解決這一問題。預(yù)處理后的N-S方程為:

其中，WT=(p，u，v，T)T，Γ為預(yù)處理矩陣，本文采用Pletcher-chen［20］預(yù)處理方法。預(yù)處理后的N-S方程進行離散得到的殘值為:

3 旋翼翼型非定常動態(tài)失速模擬

3.1 定常狀態(tài)翼型流場的計算

為驗證本文低速預(yù)處理、多重網(wǎng)格方法在流場模擬中的有效性，以NACA0012翼型為例，首先對可壓狀態(tài)Ma=0.55，α=8.34的流場進行了模擬，計算及試驗［21］結(jié)果見圖4。由圖可知，本文預(yù)處理方法能夠針對可壓狀態(tài)進行有效計算，并且結(jié)合多重網(wǎng)格方法能夠更有效地提高流場求解的效率。

然后，考慮到旋翼后行槳葉的低速狀態(tài)(來流馬赫數(shù)一般小于0.35)，針對低來流馬赫數(shù)Ma=0.1、α=3.59°的狀態(tài)進行了計算。圖5給出了收斂歷程與預(yù)處理前后的翼型壓強系數(shù)分布。由圖可以看出，采用預(yù)處理方法能夠消除可壓方程在低速狀態(tài)計算的壓力分布的振蕩，表明預(yù)處理方法能夠有效提高求解精度;同時低速預(yù)處理、多重網(wǎng)格方法也能明顯增強流場計算的收斂性。

圖4 可壓狀態(tài)翼型流場收斂歷程及壓強系數(shù)分布Fig.4Convergence history and distribution of pressure coefficient for airfoil at compressible state

圖5 低速狀態(tài)翼型流場收斂歷程與壓強系數(shù)分布Fig.5Convergence history and pressure coefficient distribution of airfoil under low-speed condition

為驗證本文方法對定常狀態(tài)失速特性的計算能力，對NACA0012翼型在Ma=0.3，Re=4×106狀態(tài)的靜態(tài)失速特性進行了數(shù)值模擬，圖6給出了翼型升力、力矩系數(shù)隨迎角的變化。由圖可知，本文計算獲得的翼型升力、力矩系數(shù)與試驗值［22］吻合較好，并且能夠有效捕捉翼型的失速迎角，表明本文方法在翼型定常狀態(tài)氣動特性模擬方面的有效性。

為驗證S-A模型對粘性計算的模擬精度，本文對RAE2822翼型的粘性流場進行了計算，計算狀態(tài)為Ma=0.75，α=3.19°，Re=6.2×106，網(wǎng)格尺寸為369× 65，法向第一層網(wǎng)格距離壁面1×10-6c，以滿足粘性計算的需要。圖7給出了S-A模型和B-L模型計算的翼型上表面摩擦力系數(shù)分布及0.9c處速度型與試驗值［21］的對比，從中可以看出，在捕捉分離區(qū)的能力上，S-A模型比B-L模型具有明顯優(yōu)勢，也表明本文計算方法的有效性。

圖6 升力、力矩系數(shù)隨迎角變化曲線Fig.6Lift and momentum coefficients vary with AOA

圖7RAE2822翼型上表面摩擦力系數(shù)分布及0.9c處速度型Fig.7Friction coefficient and velocity profile at 0.9c on the upper surface of RAE2822 airfoil

3.2 輕度失速算例分析

對于二維流動，翼型繞l/4弦長位置作正弦振蕩，迎角變化規(guī)律為:

其中，α0、αm分別為迎角的平均值和振幅，無量綱時間t*=tU∞/c，縮減頻率k=πfc/U∞。U∞為來流速度，c為翼型弦長。

為驗證本文方法的有效性，選取NACA0012翼型的非定常流場進行求解。計算狀態(tài)為典型前行槳葉工作狀態(tài):Ma=0.6，Re=4.8×106，其迎角變化規(guī)律為α=3.16°+4.59°sin(2×0.0811×t*)。

為準(zhǔn)確地模擬每個時間步翼型的非定常氣動力，將翼型一個振蕩周期分為288個時間步，子迭代為200步，并進行了5個周期的計算。圖8分別給出了采用B-L、S-A湍流模型計算收斂后的翼型升力、力矩系數(shù)的遲滯回線，由圖可以看出，S-A模型的升力、力矩系數(shù)的計算值較B-L模型計算值與試驗值［2］吻合更好，鑒于此，本文算例均采用基于S-A湍流模型的非定常N-S方程流場數(shù)值計算方法對振蕩翼型的氣動力變化進行有效模擬。

圖9給出了S-A模型計算的翼型在一個周期內(nèi)不同迎角的表面壓強系數(shù)分布，圖示箭頭向上表示翼型的上仰過程，箭頭向下表示翼型的下俯過程。由圖可以看出，本文計算值與試驗值吻合良好。

圖8 翼型升力、力矩系數(shù)遲滯回線計算值與試驗值的對比Fig.8Comparisons of lift and momentum coefficients of airfoil between calculated and experimental data

圖9 不同迎角時翼型表面壓強系數(shù)分布計算值與試驗值的對比Fig.9Comparisons of pressure coefficient distributions of airfoils at different AOA

3.3 深度失速算例分析

為探索翼型深度動態(tài)失速時流場的分布特征，本文針對旋翼后行槳葉典型工況:Ma=0.2、α0=15°、αm=10°、k=0.05的狀態(tài)進行了數(shù)值計算。圖10為計算得到的翼型升力遲滯回線與L-B模型計算值及試驗值的對比。從圖中可以看出，馬赫數(shù)低于0.3時，L-B模型計算得到的翼型升力系數(shù)與試驗值差距很大，而本文在升力系數(shù)峰值、遲滯回線面積以及氣流分離的再附的計算值與試驗值［21］更接近，表明本文建立的方法能夠較好地模擬較低來流速度下的翼型非定常動態(tài)失速特性。

圖11給出了翼型半個周期的振蕩過程中翼型表面渦量分布云圖。從圖中可以看出，在翼型上仰過程中，隨著翼型迎角的增大，翼型前緣渦逐漸加劇，并且有后緣渦產(chǎn)生。當(dāng)翼型迎角超過動態(tài)失速臨界值時，前緣產(chǎn)生脫體渦，并進一步導(dǎo)致翼型前緣的氣流分離。當(dāng)翼型從最大迎角位置下俯時，由于相對來流速度的增大，前緣脫體渦以及氣流分離現(xiàn)象逐漸加劇，當(dāng)前緣渦流過翼型后緣并進入尾跡區(qū)時，翼型上表面氣流完全分離，同時伴隨著升力的突然減小。當(dāng)翼型迎角恢復(fù)到足夠小時，氣流會在前緣再附。前緣渦的生成、沿翼型表面的對流直至脫落導(dǎo)致上仰和下俯過程翼型在迎角相同時的升力不對稱，從而引起翼型升力的遲滯效應(yīng)。

圖10 翼型升力系數(shù)遲滯回線Fig.10Lift coefficient hysteresis loop of airfoil

圖11 翼型上表面半個周期內(nèi)不同迎角時的渦量分布云圖Fig.11Vorticity contours on the upper surface of airfoil at different AOA in half a period

3.4SC1095旋翼翼型深度失速特性分析

開展了SC1095先進旋翼翼型深度動態(tài)失速時氣動特性的計算。該翼型為美國先進的通用直升機黑鷹旋翼采用的主要翼型。針對后行槳葉典型的工況Ma=0.302，α0=9.78°，αm=9.9°，k=0.099的進行了數(shù)值計算。圖12給出了計算得到的翼型升力遲滯回線的本文計算值、L-B模型計算值及試驗值的對比。從圖中可以看出，相對于L-B模型，本文計算值與試驗值［3］吻合更好，表明旋翼翼型非定常流場模擬方法能夠針對不同旋翼翼型在低馬赫數(shù)狀態(tài)的動態(tài)失速特性進行有效計算。

圖12SC1095翼型升力系數(shù)遲滯回線計算值與L-B模型及試驗值的對比Fig.12Comparisons of lift coefficient hysteresis loop of SC1095 airfoil among calculated，L-B model and experimental data

4 旋翼翼型運動參數(shù)影響分析

由旋翼翼型動態(tài)失速的運動方程(15)可知，翼型動態(tài)失速特性與俯仰運動的平均迎角、振幅及縮減頻率均有關(guān)。本文通過對NACA0012翼型在不同參數(shù)組合下的氣動力、力矩特性進行數(shù)值模擬，開展翼型動態(tài)失速特性的參數(shù)分析。

基準(zhǔn)狀態(tài)為Ma=0.3，α0=10°，αm=10°，k=0.1，進行參數(shù)分析時只改變所分析參數(shù)，其余參數(shù)與基準(zhǔn)狀態(tài)一致。

4.1 振幅的影響

針對旋翼后行槳葉的大迎角的動態(tài)失速特性，對三個不同振幅αm=5°、10°、15°分別進行了數(shù)值模擬。圖13給出了在不同振幅時翼型升力、力矩系數(shù)隨迎角的變化規(guī)律，并給出了定常狀態(tài)升力系數(shù)的計算值。從圖中可以看出，隨著振幅的增大，遲滯效應(yīng)增強。在翼型迎角增大時，不同振幅對應(yīng)的升力系數(shù)在相同迎角時計算結(jié)果基本相同，并且發(fā)生失速時的迎角相對于定常狀態(tài)均有所推遲。這是因為當(dāng)迎角大于靜態(tài)分離迎角時，分離點推遲及前緣脫體渦的影響使得升力不發(fā)生衰減，直到前緣渦從翼型后緣脫離。隨著前緣渦的脫離，由于翼型上表面氣流分離及渦誘導(dǎo)速度的損失，翼型升力系數(shù)會急劇衰減。翼型迎角減小時，同樣存在分離氣流再附著的延遲。隨著振幅的增大，氣流再附延遲更明顯，即升力遲滯回線的面積更大。與此同時，隨著翼型前緣渦的運動、脫落及氣流分離，壓力中心后移，低頭力矩有一個陡增，并且振幅的增大使得低頭力矩峰值在增大的同時發(fā)生后移;當(dāng)迎角減小時，壓力中心向翼型前緣移動，并且當(dāng)迎角減小到一定值時，會產(chǎn)生一個抬頭力矩峰值，隨著振幅的增大，遲滯效應(yīng)加劇，使得抬頭力矩的峰值發(fā)生前移。

圖13 不同迎角振幅時翼型氣動力系數(shù)遲滯回線Fig.13Aerodynamic coefficient hysteresis loop of airfoil at different amplitudes

4.2 縮減頻率的影響

直升機旋翼工作狀態(tài)的縮減頻率在0.05～0.15之間，因此，本文針對k=0.05、0.1、0.15三個狀態(tài)分別進行了數(shù)值模擬。圖14給出了在不同縮減頻率時翼型升力、力矩系數(shù)隨迎角的變化規(guī)律。由圖可以看出，隨著縮減頻率的增大，氣動力系數(shù)遲滯效應(yīng)更加強烈，遲滯回線的面積增大，升力系數(shù)峰值增大，并且失速發(fā)生位置后移。與此同時，低頭力矩的峰值隨著縮減頻率的增大而增大并且后移，抬頭力矩峰值隨縮減頻率的增大發(fā)生前移。

縮減頻率表征了旋翼翼型動態(tài)失速特性的強弱程度，頻率越大，動態(tài)特性越強烈，反之，動態(tài)特性越弱。隨著縮減頻率的增加，一方面，在翼型迎角增大過程中，分離點的滯后作用更加顯著，故動態(tài)情況下分離點的位置向后移動，即分離迎角增大;另一方面，在迎角減小的過程中，分離點再附著的延遲也更顯著，故升力系數(shù)的減小量也更大。

圖14 不同縮減頻率時翼型氣動力系數(shù)遲滯回線Fig.14Aerodynamic coefficient hysteresis loop of airfoil at different reduced frequencies

圖15 不同平均迎角時翼型氣動力系數(shù)遲滯回線Fig.15Aerodynamic coefficient hysteresis loop of airfoil at different average AOA

4.3 平均迎角的影響

針對翼型三個不同平均迎角α0=5°、10°、15°分別進行了數(shù)值模擬，圖15給出了在不同平均迎角時翼型升力、力矩系數(shù)隨迎角的變化規(guī)律。從圖中可見，隨著平均迎角的增大，升力系數(shù)峰值增大并且后移;此外，隨著平均迎角的增大，氣動力的遲滯效應(yīng)明顯增強，并且氣流分離后的再附現(xiàn)象隨著平均迎角的增大有逐漸消失的趨勢。

平均迎角同迎角振幅一起決定了翼型動態(tài)過程的幅值范圍，即動態(tài)特性的程度。最大迎角(平均迎角+正幅值)越大，動態(tài)特性越明顯。當(dāng)最大迎角略大于失速迎角時，翼型處于小分離狀態(tài)，且前緣渦的強度也不是很大，故升力系數(shù)的變化不是很劇烈，只形成一個小的遲滯環(huán)(如圖中平均迎角為5°的狀態(tài))，即輕度失速。但當(dāng)最大迎角遠(yuǎn)大于靜態(tài)失速迎角時，則翼型發(fā)生完全分離，前緣渦也得到了顯著的加強，故遲滯環(huán)的面積更大(如圖中平均迎角為10°以上的狀態(tài))，即深度失速。

5 結(jié)論

本文建立了基于運動嵌套網(wǎng)格技術(shù)的旋翼翼型動態(tài)失速數(shù)值模擬方法，并以NACA0012和SC1095旋翼翼型為研究對象，進行了翼型動態(tài)失速數(shù)值計算驗證，并與試驗值以及L-B模型計算結(jié)果進行了對比。在此基礎(chǔ)上，著重開展了NACA0012旋翼翼型動態(tài)失速特性的參數(shù)分析。綜合本文的計算分析結(jié)果，主要得出以下幾點結(jié)論:

(1)本文建立的CFD方法能夠精確模擬翼型定常、非定常狀態(tài)的氣動特性，適合旋翼后行槳葉低速、大迎角工作狀態(tài)的模擬分析;

(2)采用多重網(wǎng)格法和低速預(yù)處理法不僅能夠提高流場求解的效率，在低速狀態(tài)更能夠提高流場的求解精度;并且運用該方法在較低來流速度下能夠取得比L-B模型更精確的翼型升力和力矩特性;

(3)翼型在上仰過程中產(chǎn)生的前緣脫體渦和氣流分離會在翼型迎角減小時有所加劇，當(dāng)迎角減小到一定值時產(chǎn)生的脫體渦和氣流分離現(xiàn)象消失，氣流發(fā)生再附，并且再附時翼型迎角比分離時有所滯后，進而引起翼型氣動力的遲滯效應(yīng)。

(4)一方面，隨著翼型迎角振幅的增大，氣流分離和再附延遲現(xiàn)象更為明顯，因此氣動力遲滯效應(yīng)更加明顯，遲滯回線面積增大;另一方面，隨著振幅的增大，升力系數(shù)、低頭力矩的峰值也會增大并發(fā)生后移，抬頭力矩的峰值前移。

(5)隨著縮減頻率的增加，在翼型迎角增大過程中，分離點的滯后作用更加顯著，故動態(tài)情況下分離點的位置向后移動，即分離迎角增大。相同的，在迎角減小的過程中，分離點再附著的延遲也更顯著，故升力系數(shù)的減小量更大。

(6)最大迎角(平均迎角+正幅值)越大，動態(tài)特性越明顯。輕度失速時，翼型處于小分離狀態(tài)，前緣渦的強度不大，因此，升力系數(shù)的變化平緩，只形成一個小的遲滯環(huán);深度失速時，翼型發(fā)生完全分離，前緣渦強度很大，進而產(chǎn)生大面積遲滯環(huán)。

［1］Geissler W，Raffel M，Dietz G，et al.Helicopter aerodynamics with emphasis placed on dynamic stall springer［M］.Berlin Heidelberg，2007.

［2］Landon R H.Compendium of unsteady aerodynamic measurements［R］.AGARD-R-702，1982.

［3］McAlister K W，Pucci S L，McCroskey W J，et al.Experimental study of dynamic stall on advanced airfoil sections［R］.NASA Technical Memorandum 84245，1982，2:330-380.

［4］Leishman J G.Validation of approximate indicial aerodynamic functions for two-dimensional subsonic flow［J］.Journal of Aircraft，1988，25(10):914-922.

［5］McAlister K W，Lambert O，Pitot D.Application of the ONERA model of dynamic stall［R］.NASA TP-2399，1984.

［6］Dindar M，Kaynak U.Effect of turbulence modeling on dynamic stall of a NACA0012 airfoil［R］.AIAA 1992-0027.

［7］Leishman J G，Beddoes T S.A semi-empirical model for dynamic stall［J］.Journal of the American Helicopter Society，1989，34:3-17.

［8］Sheng M，Galbraith R A McD，Coton F N.A modified dynamic stall model for low mach numbers［C］//45thAIAA Aerospace Sciences Meeting and Exhibit，Reno，Nevada，2007:1-17.

［9］Song C Y，Xu G H.Computations of unsteady dynamic stall responses on rotor airfoils［J］.Acta Aerodynamica Sinica，2007，25(4): 461-467.(in Chinese)宋辰瑤，徐國華.旋翼翼型非定常動態(tài)失速響應(yīng)的計算［J］.空氣動力學(xué)學(xué)報，2007，25(4):461-467.

［10］邵明玉，楊茂，陳鳳明.旋翼翼型動態(tài)失速特性的數(shù)值仿真研究［J］.計算機仿真，2012，29(7):70-74.Shao M Y，Yang M，Chen F M.Numerical simulation study on dynamic stall characteristics of rotor airfoil［J］.Computer Simulation，2012，29(7):70-74.

［11］Leishman J G.Principle of helicopter aerodynamics［M］.2007，525-561.

［12］Jameson A.Time dependent calculation using multigrid with application to unsteady flows past airfoils and wings［R］.AIAA 1991-1596.

［13］Ko S，McCroskey W J.Computations of unsteady separating flows over an oscillating airfoil［J］.AIAA Journal，1997，35:1235-1238.

［14］Lee T，Gerontakos P.Investigation of flow over an oscillating airfoil［J］.J.Fluid Mech.，2004，512:313-341.

［15］Qian W Q，Leung R C K.Numerical simulation ofairfoil dynamic stall incorporating transition model［J］.Acta Aerodynamica Sinica，2008，26(1):50-55.(in Chinese)錢煒祺，Leung R C K.考慮轉(zhuǎn)捩影響的翼型動態(tài)失速數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，2008，26(1):50-55.

［16］Dumlupinar E，Murthy V R.Investigation of dynamic stall of airfoils and wings by CFD［R］.AIAA 2011-3511，2011.

［17］Wernert B E，Schreck P，Raffel S.Investigation of threedimensional dynamic stall using computational fluid dynamics［J］.AIAA Journal，2005，43:1023-1033.

［18］Yang W Q，Song B F，Song W P.Distance decreasing method forconfirming corresponding cells of overset grids and its application［J］.Acta Aeronautica et Astronautica Sinica，2009，30(2):205-212.(in Chinese)

楊文青，宋筆鋒，宋文萍.高效確定重疊網(wǎng)格對應(yīng)關(guān)系的距離縮減法及其應(yīng)用［J］.航空學(xué)報，2009，30(2):205-212.

［19］Mavriplis D J.Multigrid solution of the two dimensional Euler equations on unstructured triangular meshes［C］//AIAA 25th Aerospace Sciences Meeting，Reno，Nevada，1987:1-16.

［20］Pletcher R H，Chen K H.On solving the compressible Navier-Stokes equations for unsteady flows at very low Mach numbers［R］.AIAA-93-3368:765-775.

［21］Terry L H.Viscous transonic airfoil workshop compendium of results［R］.AIAA-87-1460，1987.

［22］Lee T，Gerontakos P.Investigation of flow over an oscillating airfoil［J］.Journal of fluid Mechanics，2004，512:313-341.

勘誤表

Zhao Guoqing，Zhao Qijun，Wang Qing
(National Key Laboratory of Science and Technology on Rotorcraft Aeromechanics，Nanjing University of Aeronautics and Astronautics，Nanjing210016，China)

A high-efficiency and high-precision CFD method for simulating the unsteady dynamic stall of rotor airfoil has been established based on moving-embedded grid and compressible RANS equations.Firstly，the generation method of viscous and orthogonal body-fitted grid around the rotor airfoil is developed by solving Poisson equations.Meanwhile，aiming at overcoming the shortcoming of spring simulation approach which may result in the distortion of grid，an improved Minimum Distance Method is proposed to generate the embedded grid around airfoil.Secondly，in order to simulate the hysteresis effect of aerodynamic forces caused by the turbulence separation and re-attachment of the flow，a high-precision method on the analysis of unsteady aerodynamic characteristics of rotor airfoil is developed by employing RANS equations and dual-time method.The S-A turbulence model is employed to capture the separation phenomenon of flow around airfoil.Thirdly，according to the conditions of low-speed inflow and high AOAs ofthe retreating blade，together with the limitation of L-B semi-empirical model on the calculation of unsteady dynamic stall of airfoil，a combination method of Pletcher-Chen preconditioning，F(xiàn)AS multigrid approaches and implicit LU-SGS scheme is established to overcome the problems of convergence difficulty and insufficient precision of compressible equations.The steady，mild and deep dynamic stall cases of NACA0012 and SC1095 rotor airfoils are calculated using this previously mentioned method，the hysteresis effect and the formation，convection，shedding of the vortical disturbance are well captured，the effectiveness of numerical simulation method on dynamic stall is verified.Finally，focus on the deep stall of NACA0012 airfoil，the influence analyses of parameters on the unsteady aerodynamic forces of rotor airfoil are carried out，and the results demonstrate that the exchanges of averaged AOA，amplitude and reduced frequency may cause a variational hysteresis effect and regularly changes of peak value of aerodynamic force.

rotor;airfoil;dynamic stall;N-S equations;moving-embedded grid;parameters analyses

序號年、卷、期、頁碼更正前更正后1《空氣動力學(xué)學(xué)報》2014年32卷第6期第740頁中文摘要第三行在此研究中，翼型做沉標(biāo)和俯仰運動在此研究中，翼型做沉浮和俯仰運動2《空氣動力學(xué)學(xué)報》2014年32卷第6期編輯部2014年紀(jì)要第五行……審稿和編輯出版工作，從2014第4期起連續(xù)刊載……審稿和編輯出版工作，從2014第5期起連續(xù)刊載3《空氣動力學(xué)學(xué)報》2014年32卷第6期編輯部2014年紀(jì)要第十六行目前已經(jīng)準(zhǔn)備就緒并將在2014年實施的項目有目前已經(jīng)準(zhǔn)備就緒并將在2015年實施的項目有

V211.52;V211.3

Adoi:10.7638/kqdlxxb-2013.0010

Simulations and parametric analyses on unsteady dynamic stall characteristics of rotor airfoil based on CFD method

0258-1825(2015)01-0072-10

2013-01-23;

2013-02-25

國家自然科學(xué)基金項目(11272150);江蘇高校優(yōu)勢學(xué)科建設(shè)工程資助項目

趙國慶(1984-)，男，山東五蓮人，博士研究生，研究方向:旋翼流動控制、旋翼計算流體力學(xué)、直升機空氣動力學(xué).E-mail: zgq198495@nuaa.edu.cn

招啟軍(1977-)，男，教授、博士生導(dǎo)師，研究方向:直升機計算流體力學(xué)、直升機空氣動力學(xué)及流動控制、氣動噪聲、總體設(shè)計等.E-mail:zhaoqijun@nuaa.edu.cn

趙國慶，招啟軍，王清.旋翼翼型非定常動態(tài)失速特性的CFD模擬及參數(shù)分析［J］.空氣動力學(xué)學(xué)報，2015，33(1):72-81.

10.7638/kqdlxxb-2013.0010.Zhao G Q，Zhao Q J，Wang Q.Simulations and parametric analyses on unsteady dynamic stall characteristics of rotor airfoil based on CFD method［J］.Acta Aerodynamica Sinica，2015，33(1):72-81.