付煒嘉, 李杰, 婁琪琳

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安　710072)

?

基于高階格式的旋轉(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬研究

付煒嘉, 李杰, 婁琪琳

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安710072)

摘要：基于結(jié)構(gòu)網(wǎng)格動(dòng)態(tài)面搭接技術(shù)、五階Roe-WENO格式求解非定常RANS方程，開發(fā)了用于旋轉(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬的計(jì)算分析程序，并針對(duì)典型旋轉(zhuǎn)類算例進(jìn)行了分析驗(yàn)證。對(duì)Caradonna-Tung旋翼的計(jì)算結(jié)果表明：計(jì)算得到的旋翼槳葉表面壓力分布與實(shí)驗(yàn)值吻合良好，計(jì)算程序可以用于旋翼懸停、簡(jiǎn)單前飛狀態(tài)的氣動(dòng)力數(shù)值模擬與分析；同時(shí)在計(jì)算網(wǎng)格分布相同情況下，五階Roe-WENO格式較三階Roe-MUSCL格式對(duì)于旋翼尾跡的分辨率更高，計(jì)算得到的旋翼槳尖渦尾跡更加清晰。對(duì)螺旋槳機(jī)翼干擾問題的研究結(jié)果表明，所開發(fā)的計(jì)算程序?qū)τ趧傂孕D(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬具有很好的適用性，對(duì)于尾跡的捕捉具有較高的精度，滿足工程實(shí)際應(yīng)用的要求。

關(guān)鍵詞：結(jié)構(gòu)網(wǎng)格動(dòng)態(tài)面搭接技術(shù);五階Roe-WENO格式;旋翼;螺旋槳;非定常流動(dòng)

當(dāng)前，計(jì)算流體力學(xué)(CFD)技術(shù)已經(jīng)成為空氣動(dòng)力學(xué)問題研究的主要手段，被廣泛地用于工程實(shí)際問題的計(jì)算分析。旋轉(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬一直是空氣動(dòng)力學(xué)研究的難點(diǎn)與熱點(diǎn)，其難點(diǎn)主要包括以下2個(gè)方面：①如何采用合理的網(wǎng)格技術(shù)實(shí)現(xiàn)旋轉(zhuǎn)運(yùn)動(dòng)；②如何提高尾跡模擬的精度。

旋轉(zhuǎn)類問題在工程中隨處可見，包括直升機(jī)旋翼、螺旋槳、風(fēng)力機(jī)葉片、航空發(fā)動(dòng)機(jī)葉片等問題。目前，研究者大多通過求解雷諾平均N-S方程對(duì)這類問題進(jìn)行數(shù)值模擬研究[1-4]。從國(guó)內(nèi)外有關(guān)旋轉(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬研究文獻(xiàn)來看，大致有2種研究?jī)A向：①比較側(cè)重于旋轉(zhuǎn)物體本身的氣動(dòng)力數(shù)值模擬分析，例如旋翼槳葉、螺旋槳槳葉和風(fēng)力機(jī)葉片的表面載荷分布及其相關(guān)力系數(shù)的計(jì)算。②比較側(cè)重于旋轉(zhuǎn)物體尾跡的高精度數(shù)值模擬，通過對(duì)數(shù)值方法的改進(jìn)完善來提高尾跡渦模擬的模擬精度，捕捉到更多的流動(dòng)細(xì)節(jié)。后者逐漸成為旋轉(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬發(fā)展趨勢(shì)，對(duì)于螺旋槳問題，螺旋槳的尾跡的模擬精度對(duì)于螺旋槳滑流與機(jī)翼或飛機(jī)其他部件的氣動(dòng)干擾研究有重要的意義。對(duì)于旋翼問題，旋翼槳尖渦的模擬精度對(duì)于深入研究旋翼繞流機(jī)理、槳渦干擾、旋翼機(jī)身干擾以及旋翼氣動(dòng)噪聲預(yù)測(cè)至關(guān)重要[5]。由此可見，有必要發(fā)展一套適用于旋轉(zhuǎn)問題的非定常高精度數(shù)值模擬研究的計(jì)算分析程序。

本文采用結(jié)構(gòu)網(wǎng)格動(dòng)態(tài)面搭接技術(shù)實(shí)現(xiàn)物體的旋轉(zhuǎn)，基于五階Roe-WENO格式和非定常RANS方程，開發(fā)了一套用于旋轉(zhuǎn)類問題非定常流動(dòng)數(shù)值模擬研究的計(jì)算分析程序，并針對(duì)旋翼和螺旋槳等典型算例進(jìn)行計(jì)算分析來驗(yàn)證該計(jì)算程序的可靠性。

1數(shù)值計(jì)算方法

1.1控制方程

本文對(duì)旋轉(zhuǎn)類問題流場(chǎng)的數(shù)值模擬是基于慣性笛卡爾坐標(biāo)系開展的，在慣性笛卡爾坐標(biāo)下，積分形式的非定常N-S方程組可寫作：

(1)

式中，Ω為控制體,S為控制體邊界,n為S微元的外法向單位矢量,其余各參數(shù)具體含義見參考文獻(xiàn)[6]。采用有限體積法離散求解上述方程,黏性項(xiàng)采用二階中心差分格式離散,流場(chǎng)分析基于全湍流假設(shè),采用k-ωSST兩方程湍流模型進(jìn)行湍流數(shù)值模擬,時(shí)間推進(jìn)采用雙時(shí)間迭代方法。

1.2空間離散格式

本文中所采用的數(shù)值計(jì)算方法對(duì)無黏通量項(xiàng)離散選用基于迎風(fēng)型通量差分分裂的Roe格式。針對(duì)Roe格式構(gòu)造了基于改進(jìn)的WENO方法對(duì)交界面兩邊狀態(tài)參量進(jìn)行重構(gòu)的五階Roe-WENO格式。

對(duì)于Roe格式,在三維貼體正交坐標(biāo)系下,穿過網(wǎng)格單元界面的無黏通量可表示為

(2)

從(2)式可以看出,在采用Roe格式對(duì)N-S方程空間項(xiàng)進(jìn)行離散求解時(shí),首先要求得網(wǎng)格單元體交界面兩側(cè)變量的值qL和qR。通常對(duì)網(wǎng)格交界面左右兩側(cè)的變量是通過插值得到的。利用基于平均的網(wǎng)格中心點(diǎn)處的值來插值得到網(wǎng)格交界面處的變量值。因此,對(duì)交界面處變量插值的精度就決定了整個(gè)空間離散格式的精度。

由MUSCL插值方法得到的三階迎風(fēng)差分格式為

(3)

本文中采用改進(jìn)的五階WENO格式[7]對(duì)交界面兩側(cè)狀態(tài)參量進(jìn)行重構(gòu),以一維標(biāo)量模型方程為例,五階精度的邊界外推值形式為

(4)

(5)

ωk為非線性加權(quán)因子,其定義如下

(6)

1.3動(dòng)態(tài)面搭接網(wǎng)格技術(shù)通量傳遞

如何保證流場(chǎng)信息在搭接面準(zhǔn)確傳遞,是動(dòng)態(tài)面搭接技術(shù)的關(guān)鍵。對(duì)于圖1給出的二維搭接面示意圖,(i,j)表示Zone1中網(wǎng)格的格心,(l,m)表示Zone2中網(wǎng)格的格心。2塊網(wǎng)格搭接面的網(wǎng)格點(diǎn)分別可表示為{xi+1/2,j;j=1,jmax}和{xl-1/2,m;m=1,mmax}。定義F為x方向上的通量,為了保證搭接面處的通量守恒,應(yīng)有

圖1　二維搭接面示意圖

(7)

(8)

(9)

(10)

可以得到(7)式搭接面處通量守恒的離散形式

(11)

三維與二維情況最大的不同處在于,對(duì)于二維情況,2塊網(wǎng)格的公共區(qū)域是一條線,而三維情況下公共區(qū)域是面。對(duì)于圖2給出的三維搭接面示意圖,Zone1的網(wǎng)格點(diǎn)用離散點(diǎn)集表示為

顯然Zone1的搭接面為ξ=常數(shù),可以寫為

Zone2的搭接面也為ξ=常數(shù),可以寫為

在搭接面處的通量守恒方程為如下形式

(12)

相應(yīng)的,(12)式的離散形式可以表示為

(13)

圖2　三維搭接面示意圖

2旋翼流場(chǎng)數(shù)值模擬

2.1Caradonna-Tung懸停流場(chǎng)數(shù)值模擬

Caradonna和Tung在1981年發(fā)表了一組旋翼壓力測(cè)量的實(shí)驗(yàn)結(jié)果,該旋翼的實(shí)驗(yàn)裝置及外形見參考文獻(xiàn)[11],該旋翼包括2片槳葉,每片槳葉都是以NACA0012翼型為剖面的無扭轉(zhuǎn)直機(jī)翼。槳盤半徑R為1.143 m,展弦比為6。

圖3　網(wǎng)格分區(qū)示意圖

采用ICEM CFD軟件對(duì)整個(gè)計(jì)算域進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分,整個(gè)流場(chǎng)分為旋轉(zhuǎn)區(qū)域及背景靜止區(qū)域2個(gè)部分,2個(gè)區(qū)域通過搭接面?zhèn)鬟f信息,旋翼物面附近采用“O-O”型拓?fù)浣Y(jié)構(gòu),其余地方采用“H-H”型拓?fù)浣Y(jié)構(gòu),整個(gè)計(jì)算域網(wǎng)格節(jié)點(diǎn)總數(shù)為1 500萬。圖3給出了旋翼計(jì)算網(wǎng)格分區(qū)示意圖,靜止區(qū)域和旋轉(zhuǎn)區(qū)域之間通過搭接面?zhèn)鬟f信息。選取的計(jì)算條件為:槳尖馬赫數(shù)Mtip=0.44,槳葉槳距角θc=8′,基于槳尖弦長(zhǎng)的雷諾數(shù)為Re=1.92×106。

圖4給出了旋翼槳葉不同位置截面上的壓力分布對(duì)比,虛線為三階Roe-MUSCL格式的計(jì)算結(jié)果,實(shí)線為五階Roe-WENO格式的計(jì)算結(jié)果。

圖4　槳葉表面壓力分布

2種格式計(jì)算所得壓力分布均與實(shí)驗(yàn)數(shù)據(jù)吻合得較好,僅槳葉上表面的壓力分布略有差異,五階Roe-WENO格式計(jì)算所得的前緣壓力峰值與實(shí)驗(yàn)值更加接近,尤其是在靠近槳尖的附近區(qū)域。圖5給出了渦量值為0.2時(shí)旋翼渦量等值面對(duì)比圖?？梢钥吹?在渦量值相同情況下,三階Roe-MUSCL格式計(jì)算捕捉到的旋翼槳尖渦僅發(fā)展到了尾跡角240°左右,而五階Roe-WENO格式計(jì)算捕捉到的旋翼槳尖渦發(fā)展到了尾跡角700°左右,從五階Roe-WENO格式計(jì)算結(jié)果還能看到槳尖渦在向下發(fā)展的同時(shí)具有一定的收縮性,其渦核位置隨著尾跡角增大向內(nèi)收縮。

圖5　渦量等值面對(duì)比圖

圖6給出了2種格式計(jì)算得到的旋翼空間截面上渦量等值線對(duì)比圖。計(jì)算采用的Caradonna-Tung旋翼模型由2片槳葉構(gòu)成,因此上下2個(gè)槳尖渦尾跡角相差180°?？梢钥吹?Roe-MUSCL格式捕捉到的旋翼槳尖渦沿尾跡角發(fā)展了500°左右后就逐漸被耗散消失,而五階Roe-WENO格式捕捉的旋翼槳尖渦沿尾跡角發(fā)展了900°后依然清晰可見。由于Roe-MUSCL格式的數(shù)值耗散較大,計(jì)算得到的槳尖渦外形有明顯失真,從尾跡角180°開始上下相鄰的槳尖渦粘連在一起,尾跡角360°后的槳尖渦形態(tài)基本難以辨認(rèn),同一尾跡角處的槳尖渦渦核位置與五階Roe-WENO格式有明顯差異。五階Roe-WENO格式計(jì)算捕捉到的槳尖渦結(jié)構(gòu)形態(tài)清楚,在同一尾跡角處的槳尖渦渦量明顯強(qiáng)于Roe-MUSCL格式。

圖6　空間截面渦量等值線對(duì)比圖

從以上對(duì)比可以看出,本文所采用的五階Roe-WENO格式數(shù)值耗散較低,對(duì)于旋翼槳尖渦的渦強(qiáng)度及發(fā)展過程預(yù)測(cè)要優(yōu)于三階Roe-MUSCL格式。

2.2Caradonna-Tung前飛流場(chǎng)數(shù)值模擬

本節(jié)主要檢驗(yàn)本文所開發(fā)的計(jì)算程序?qū)τ谛砬帮w狀態(tài)的數(shù)值模擬能力。計(jì)算采用的模型仍然為Caradonna-Tung旋翼,為了盡可能捕捉旋翼前飛過程中的槳尖渦,本文對(duì)尾跡區(qū)內(nèi)的網(wǎng)格進(jìn)行了加密,空間截面網(wǎng)格見圖7,整個(gè)計(jì)算域網(wǎng)格節(jié)點(diǎn)數(shù)約為1 300萬。

圖7　網(wǎng)格示意圖

空間離散采用五階Roe-WENO格式,選取的計(jì)算狀態(tài)為:槳尖馬赫數(shù)Mtip=0.44,槳葉槳距角θc=8°,基于槳尖弦長(zhǎng)的雷諾數(shù)為Re=3.66×106,前進(jìn)比μ=0.3。該狀態(tài)下旋翼做有升力前飛運(yùn)動(dòng),槳葉前飛過程中無揮舞及變距現(xiàn)象。圖8給出了旋翼槳葉展向r/R=0.89處不同槳葉方位角時(shí)的壓力分布圖,將本文計(jì)算結(jié)果與參考文獻(xiàn)[5]的計(jì)算結(jié)果進(jìn)行了對(duì)比,可以看到兩者吻合得較好,再一次證明了本文采用數(shù)值方法的可靠性。

圖8　壓力分布對(duì)比圖

圖9給出了旋翼分別在槳葉方位角90°時(shí)的渦量等值面圖,可以看到采用五階Roe-WENO格式可以清晰捕捉到旋翼前飛過程中的尾渦結(jié)構(gòu)。旋翼的槳尖渦在前飛過程中隨著來流向后發(fā)展,分別在槳葉前行及后退一側(cè)形成了一道集中渦束,在槳葉前行一側(cè),由于槳葉旋轉(zhuǎn)速度與來流速度的疊加,槳尖渦強(qiáng)度較強(qiáng),而在槳葉后退一側(cè),槳葉旋轉(zhuǎn)速度與來流速度方向相反,槳葉當(dāng)?shù)氐暮铣伤俣认鄬?duì)較小,導(dǎo)致向后發(fā)展的槳尖渦強(qiáng)度要弱。從圖中還能觀察到旋翼的槳-渦干擾現(xiàn)象,旋翼槳葉在前飛過程中都處于自身的尾渦干擾之中。

圖9　渦量等值面圖

3螺旋槳滑流機(jī)翼氣動(dòng)干擾問題

本節(jié)主要針對(duì)螺旋槳機(jī)翼組合構(gòu)型進(jìn)行非定常數(shù)值流動(dòng)分析,進(jìn)一步驗(yàn)證本文所采用的結(jié)構(gòu)網(wǎng)格動(dòng)態(tài)面搭接技術(shù)的可靠性。計(jì)算所采用的模型見圖10,該模型與德宇航(DLR)Stuermer等研究的螺旋槳機(jī)翼組合模型AGARD2類似[4],不同的是Stuermer等采用的螺旋槳包含4片槳葉,本文采用的螺旋槳包含6片槳葉。所采用的機(jī)翼為平直機(jī)翼,采用的翼型為NACA0012對(duì)稱翼型,機(jī)翼安裝角為0°,機(jī)翼弦長(zhǎng)為3 m,展長(zhǎng)為12 m。

圖10　螺旋槳機(jī)翼組合構(gòu)型網(wǎng)格拓?fù)浣Y(jié)構(gòu)圖

整個(gè)計(jì)算域網(wǎng)格節(jié)點(diǎn)總數(shù)約為900萬。計(jì)算的狀態(tài)為:來流迎角α=0°,來流馬赫數(shù)Ma=0.15,螺旋槳轉(zhuǎn)速為n=1 075 r/min,槳葉槳盤直徑D=4 m。圖11給出了螺旋槳機(jī)翼構(gòu)型在有無螺旋槳情況下升力沿機(jī)翼展向的分布圖,中心體兩側(cè)機(jī)翼的升力分布呈反對(duì)稱分布,右側(cè)由于螺旋槳槳葉向上旋轉(zhuǎn)使得當(dāng)?shù)赜窃龃髲亩a(chǎn)生正的升力,左側(cè)則反之。圖12給出了螺旋槳兩側(cè)對(duì)稱的2個(gè)截面展向站位的壓力分布對(duì)比。以上計(jì)算得到的結(jié)論與Stuermer等得到的結(jié)論一致,進(jìn)一步驗(yàn)證了本文計(jì)算分析程序的可靠性。

圖11　有無滑流機(jī)翼展向升力系數(shù)對(duì)比

圖12　有無螺旋槳滑流機(jī)翼表面壓力分布對(duì)比圖

4結(jié)論

針對(duì)旋轉(zhuǎn)類問題非定常流動(dòng)的數(shù)值模擬,本文基于結(jié)構(gòu)網(wǎng)格動(dòng)態(tài)面搭接技術(shù)、五階Roe-WENO格式和非定常RANS方法開發(fā)了一套計(jì)算分析程序,并采用該計(jì)算程序?qū)Φ湫托砗吐菪龢憷M(jìn)行了數(shù)值分析驗(yàn)證。分析結(jié)果表明:

1) 本文開發(fā)的計(jì)算程序具有較好的通用性,不僅可以用于旋翼懸停和簡(jiǎn)單前飛狀態(tài)的數(shù)值模擬,對(duì)于螺旋槳及其他剛性旋轉(zhuǎn)類問題也具有很好的適用性。

2) 所開發(fā)的計(jì)算程序具有較高的精度,在保證氣動(dòng)力分析可靠的同時(shí),大大提高了旋轉(zhuǎn)類問題的尾跡渦模擬精度,對(duì)于深入研究旋轉(zhuǎn)類問題的流動(dòng)機(jī)理具有重要的意義。

參考文獻(xiàn)：

[1]王適存, 徐國(guó)華. 直升機(jī)旋翼空氣動(dòng)力學(xué)的發(fā)展[J]. 南京航空航天大學(xué)學(xué)報(bào), 2001, 33(3):203-211

Wang Shicun, Xu Guohua. Progress of Helicopter Rotor Aerodymaics[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2001, 33(3): 203-211 (in Chinese)

[2]Srinivasan G R, McCroskey W J. Navier-Stokes Calculations of Hovering Rotor Flowfields[J]. Journal of Aircraft, 1988, 25(10): 865-874

[3]肖中云. 旋翼流場(chǎng)數(shù)值模擬方法研究[D]. 中國(guó)空氣動(dòng)力研究與發(fā)展中心, 2007

Xiao Zhongyun. Investigation of Computational Modeling Techniques for Rotor Flowfields[D]. China Aerodynamics Research and Development Center, 2007 (in Chinese)

[4]Stuermer A. Unsteady CFD Simulations of Propeller Installation Effects[R]. AIAA-2006-4969

[5]徐麗. 直升機(jī)旋翼非定常流動(dòng)的低耗散高精度數(shù)值方法研究[D]. 上海: 上海大學(xué), 2010

Xu Li. Numerical Research on High-Order Schemes with Low Dissipation for Unsteady Helicopter Rotor Flows[D]. Shanghai, Shanghai University, 2010 (in Chinese)

[6]周志宏.復(fù)雜流場(chǎng)NS方程及飛機(jī)結(jié)冰數(shù)值模擬[D]. 西安:西北工業(yè)大學(xué),2011

Zhou Zhihong. Navier-Stokes Analysis in Complex Configurations and Icing Simulation on Airplanes[D]. Xi′an, Northwestern Polytechnical University, 2011 (in Chinese)

[7]Borges R, Carmona M, Costa B, et al. An Improved Weighted Essentially Non-Oscillatory Scheme for Hyperbolic Conservation Laws[J]. Journal of Computational Physics, 2008, 337: 3191-3211

[8]Lerat A, Wu Z N. Stable Conservative Multi-Domain Treatments for Implicit Euler Solvers[J]. Journal of Computational Physics, 1996, 123: 45-64

[9]Thomas James L. A Patched-Grid Algorithm for Complex Configurations Directed towards the F-18 Aircraft[R]. AIAA-1989-

0121

[10] Rai M M. A Relaxation Approach to Patched-Grid Calculations with Euler Equations[R]. AIAA-1985-0295

[11] Caradonna F X, Tung C. Experimental and Analytical Studies of a Model Helicopter Rotor in Hover[J]. Vertica, 1981, 5(2):149-161

收稿日期：2015-10-22基金項(xiàng)目：國(guó)家自然科學(xué)基金(11172240)、國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(2015CB755800)及航空科學(xué)基金(2014ZA53002)資助

作者簡(jiǎn)介：付煒嘉(1985—)，西北工業(yè)大學(xué)博士研究生，主要從事飛行器設(shè)計(jì)空氣動(dòng)力學(xué)與計(jì)算流體力學(xué)的研究。

中圖分類號(hào)：V211.3

文獻(xiàn)標(biāo)志碼：A

文章編號(hào):1000-2758(2016)03-0431-06

Unsteady Numerical Simulation of Rotor′s Flowfields Based on High-Order Scheme

Fu Weijia, Li Jie, Lou Qilin

(College of Aeronautics, Northwestern Polytechnical University, Xi′an 710072, China)

Abstract:A computational analysis procedure, which based on structure dynamic patched grid, the fifth-order Roe-WENO scheme and unsteady RANS method, was developed for the numerical simulation of the unsteady flowfiled of the rotor. Typical rotor examples were analyzed to verify the reliability of the procedure. Calculation results of the Caradonna-Tung rotor show that the calculated pressure distribution of the rotor is in good agreement with the experimental values. The developed procedure can be used for the numerical analysis of the aerodynamic performance of the rotor in hover or forward flight. In the case of the same computational grid distribution, the fifth-order Roe-WENO scheme has a higher resolution for the simulation of the rotor′s wake than that of the third-order Roe-MUSCL scheme and the wake of the tip vortex can be captured more clearly. The results of the propeller wing interference show that the developed procedure has good applicability for the unsteady simulation of the rigid rotation type of problem. The procedure developed by this paper has high accuracy of wake vortex simulation, it can meet the requirements of engineering applications.

Keywords:Structure dynamic patched grid, fifth order Roe-WENO scheme, rotors, propellers, unsteady flow