石勝慶，張祥謀，程鈺鋒

(1.江西上饒海軍91551部隊，江西 上饒 334000;2.北京裝備學(xué)院研究生院，北京 101416)

0 引言

自1903年萊特兄弟發(fā)明飛機以來，螺旋槳動力系統(tǒng)便作為飛行器的主要動力裝置之一而被人們廣泛認(rèn)識和深入研究［1］。

近年來，國內(nèi)對于螺旋槳氣動方面的各種研究也逐漸增多。2008 年，中科院的聶營、王生等人［2，3］采用Gambit軟件對螺旋槳進(jìn)行幾何建模，再用Fluent軟件基于滑移網(wǎng)格，仿真研究了臨近空間螺旋槳的氣動性能。2008年，西工大的許建華、宋文萍等人［4］采用雷諾平均N-S方程和嵌套網(wǎng)格技術(shù)對美國國家可再生能源實驗室的CER實驗型風(fēng)力機葉片進(jìn)行了仿真研究，湍流方程是B-L代數(shù)模型;計算中采用了有限體積空間離散法和改進(jìn)型隱式LU-SGS時間推進(jìn)格式，所用的嵌套網(wǎng)格技術(shù)有效地捕捉了螺旋槳的尾渦;在此基礎(chǔ)上，2009年［5］，他們采用在不同粗細(xì)網(wǎng)格上消除不同頻率誤差加速解的收斂的多重網(wǎng)格技術(shù)，仿真研究了螺旋槳側(cè)流粘性流場，其中多重網(wǎng)格采用非線性方程的全近似格式(FAS)。2010年，西工大的羅淞、楊永等人［6］采用多塊點對點網(wǎng)格生成技術(shù)，生成了螺旋槳的空間計算網(wǎng)格，分別研究了歐拉方程組和N-S方程組對不同進(jìn)距比下螺旋槳運動過程的仿真結(jié)果，與實驗進(jìn)行對比后得出歐拉方程比較適合于螺旋槳的工程設(shè)計。

以螺旋槳動力系統(tǒng)作為動力裝置的飛行器在上升、下降或轉(zhuǎn)彎的過程中，合成氣流的實際速度方向與槳盤所在面并不垂直，這時螺旋槳處于有偏角的側(cè)流狀態(tài)，其運動特點會發(fā)生改變。本文基于滑移網(wǎng)格模型，通過求解三維非定常N-S方程，數(shù)值模擬了螺旋槳側(cè)飛狀態(tài)時的運動過程，仿真研究了側(cè)偏角對螺旋槳氣動性能的影響。

1 數(shù)學(xué)模型

對于N-S方程，連續(xù)方程、動力方程和能量方程的通用形式可以寫成如下形式。

其中:ρ是氣體密度，U是速度矢量，φ是通用變量，Γ是廣義擴散系數(shù)，S是廣義源項。對于連續(xù)方程、動力方程和能量方程，φ分別為1、ui和T;Γ分別為0、μ 和k/cp;S分別為0、-?p/?xi和ST。ui是速度分量，T是溫度，μ是粘性，k是流體的傳熱系數(shù)，cp是比熱容，ST是粘性耗散項，即流體的內(nèi)熱源及由于粘性作用流體機械能轉(zhuǎn)換為熱能的部分。

理想氣體狀態(tài)方程為:

式中R是氣體常數(shù)。

湍流模型采用RNGk-ε兩方程模型，它是基于k-ε標(biāo)準(zhǔn)兩方程的湍流模型，采用一種叫做重正規(guī)化群的數(shù)學(xué)方法對N-S方程進(jìn)行暫態(tài)推理得到的改進(jìn)型k-ε兩方程湍流模型，它是由V.Yakhot和 S.A.Orszag于1986 年提出并逐步完善的［7，8］，其基本思想是認(rèn)為，在流場中小渦是各向同性的，處于統(tǒng)計定常和統(tǒng)計平衡的狀態(tài)。RNGk-ε湍流模型考慮了低雷諾數(shù)流動粘性，改進(jìn)了標(biāo)準(zhǔn)k-ε模型的高雷諾數(shù)性質(zhì)，并且提供了Prandtl數(shù)的解析公式，考慮了湍流漩渦，因此更加適合于雷諾數(shù)不是很高和帶有強漩渦運動狀態(tài)的數(shù)值仿真。

采用耦合求解器，首先同時求解連續(xù)方程、動力方程和能量方程，然后求解湍流方程。在耦合算法中使用隱式格式，即通過求解方程組的形式求解流場變量，它是使用塊Gauss-Seidel法與AMG法(Algebraic Multi-Grid，代數(shù)多重網(wǎng)格法)聯(lián)合完成的。

采用二階精度的有限體積AUSM(Advection Upstream Splitting Method)離散格式對粘性流體的控制方程和湍流方程進(jìn)行空間離散。另外，本文采用了滑移網(wǎng)格模型來處理螺旋槳的運動過程，詳見文獻(xiàn)［9］。

采用滑移網(wǎng)格模型處理螺旋槳的旋轉(zhuǎn)運動過程。滑移網(wǎng)格是在動參考系模型和混合面法的基礎(chǔ)上發(fā)展起來的，常用于風(fēng)車、轉(zhuǎn)子、螺旋槳等運動的仿真研究。在滑動網(wǎng)格模型計算中，流場中至少存在兩個網(wǎng)格區(qū)域，每一個區(qū)域都必須有一個網(wǎng)格界面與其他區(qū)域連接在一起。網(wǎng)格區(qū)域之間沿界面做相對運動?；瑒泳W(wǎng)格模型允許相鄰網(wǎng)格間發(fā)生相對運動，而且網(wǎng)格界面上的節(jié)點無需對齊，即網(wǎng)格交界面是非正則的。在使用滑動網(wǎng)格模型時，計算網(wǎng)格界面上的通量需要考慮到相鄰網(wǎng)格間的相對運動以及由運動形成的重疊區(qū)域的變化過程。詳見文獻(xiàn)［9］。

2 螺旋槳模型及邊界條件

2.1 模型及網(wǎng)格

由葉素理論可知，葉素的氣動性能對螺旋槳整體氣動性能有很大影響。由于臨近空間環(huán)境大氣密度較低，葉素雷諾數(shù)較低，因此本文選擇低雷諾數(shù)翼型SD-8000PT作為螺旋槳的葉素［10］。螺旋槳的安裝角和弦長詳見文獻(xiàn)［10］。兩槳葉螺旋槳直徑D為1.5 m，兩個槳葉之間的槳轂簡化為圓柱面，其半徑和長度分別為0.01D和0.05D。

計算區(qū)域是一個長10D、直徑8D的圓柱體。速度入口距螺旋槳4D，給定氣流速度及總溫，壓力出口距螺旋槳6D，給定總溫和總壓;遠(yuǎn)場距螺旋槳轉(zhuǎn)軸4D，給定氣流速度、總壓及總溫。

由于滑移網(wǎng)格模型允許相鄰網(wǎng)格之間發(fā)生相對運動，而且網(wǎng)格界面上的點無需對齊，即網(wǎng)格是非正則的。利用這一特點，可以更好地分布網(wǎng)格的疏密度，既保證了計算流場所需要網(wǎng)格數(shù)，又使網(wǎng)格總數(shù)減少，從而節(jié)約計算資源。

對螺旋槳的仿真而言，螺旋槳近區(qū)域流場變化劇烈，所以應(yīng)該加密網(wǎng)格，遠(yuǎn)離螺旋槳的區(qū)域流場較為平緩，對網(wǎng)格數(shù)目要求不高，所以網(wǎng)格數(shù)目較少?；诖耍疚膶⒂嬎銋^(qū)域分為氣流入口區(qū)域、旋轉(zhuǎn)區(qū)域和氣流出口區(qū)域3個小的計算區(qū)域。氣流入口和出口區(qū)域采用TTM［11］方法生成結(jié)構(gòu)網(wǎng)格，并在各自靠近螺旋槳的一端加密網(wǎng)格。旋轉(zhuǎn)區(qū)域采用TGrid網(wǎng)格劃分法生成非結(jié)構(gòu)網(wǎng)格，螺旋槳表面網(wǎng)格節(jié)點之間相隔1 mm。

圖1 計算區(qū)域及網(wǎng)格示意圖

2.2 邊界條件

以螺旋槳動力系統(tǒng)作為動力裝置的飛行器在上升、下降或轉(zhuǎn)彎的過程中，合成氣流的實際速度方向與槳盤所在面并不垂直，這時螺旋槳處于有偏角的側(cè)流狀態(tài)，其運動特點會發(fā)生改變。本文主要研究螺旋槳處于此類側(cè)流工況下的氣動性能。

圖2是笛卡爾坐標(biāo)系下飛行器運動方向即螺旋槳速度方向與槳盤關(guān)系的示意圖。圖中，紅色加粗線段CO是飛行器的運動方向即飛行器和螺旋槳的速度方向，槳盤所在面是坐標(biāo)面y-z，BO是CO在x-z平面上的投影，O是坐標(biāo)系的原點。假設(shè)飛行器的運動速度矢量為V，螺旋槳速度方向與x-y平面的角度α=∠AOB，螺旋槳速度方向與x-z平面的角度β=∠BOC，則螺旋槳速度矢量在x、y、z三個方向上的分量分別如下:

圖2 笛卡爾坐標(biāo)系中螺旋槳速度矢量示意圖

本節(jié)主要研究在20km的臨近空間環(huán)境下保持螺旋槳速度|V|為5 m/s和轉(zhuǎn)速為300 rpm不變的條件下，偏角β=45°時，螺旋槳的氣動性能隨偏角α的變化規(guī)律。

3 仿真結(jié)果及分析

本文在保持偏角β為45°不變的情況下，數(shù)值模擬了不同偏角α的工況下螺旋槳的運動過程。表1是不同偏角條件下螺旋槳的氣動參數(shù)，其中Tx、Ty、Tz分別表示螺旋槳在x、y、z三個方向上拉力的平均值，Mx、My表示螺旋槳拉力在x、y方向上扭矩的平均值，ηx是螺旋槳在x方向上的效率，效率計算公式中的拉力取螺旋槳在x方向上拉力即軸向拉力的平均值，扭矩取螺旋槳在x軸上扭矩即軸向扭矩的平均值。

由表1可見，螺旋槳在x方向上的拉力分量Tx隨著偏角α的增大逐漸增大，由公式(3)可知，當(dāng)偏角β=45°時，螺旋槳在x方向上的速度分量vx隨著偏角α的增大逐漸減小，當(dāng)螺旋槳的轉(zhuǎn)速不變前進(jìn)速度減小時，螺旋槳產(chǎn)生的拉力增大［9］，所以Tx隨著偏角α的增大逐漸增大。螺旋槳在y方向上的拉力分量Ty有隨α增大而逐漸減小的趨勢。Ty是螺旋槳側(cè)飛時偏角帶來的阻力，Ty越大側(cè)向阻力越大，對螺旋槳的飛行越不利。螺旋槳在z方向上的拉力分量Tz隨α的變化規(guī)律不是很明顯。比較可見，Tx的約為Ty的2倍，約為Tz的100倍，可見，螺旋槳在z方向的拉力分量可能是力學(xué)噪聲，即可以認(rèn)為Tz為0。

由扭矩分布可知，槳葉對x軸的扭矩Mx隨α的增大逐漸增大，槳葉隨y軸的扭矩My隨α的增大有逐漸減小的趨勢，這可能是由拉力變化而引起的。由效率的分布可見，不同偏角工況中螺旋槳的效率變化不大，最大變化幅度小于0.0051。

圖3(a)是偏角β=45°、偏角α=0°時槳葉周圍的流線分布圖。由圖3(a)可見，當(dāng)槳徑大于0.3 m時槳葉迎風(fēng)面就會發(fā)生流動分離現(xiàn)象，且分離渦的尺度隨著槳徑的增大逐漸增大，雖然r=0.1 m葉素處于負(fù)攻角工況，但槳葉的背風(fēng)面并沒有發(fā)生流動分離現(xiàn)象。

表1 不同偏角下螺旋槳的氣動參數(shù)

圖3 槳葉周圍流線分布比較圖

圖3(b)是偏角 β=45°、偏角 α=30°時槳葉周圍的流線分布圖。由圖3(b)可見，在該工況下槳葉上發(fā)生流動分離的槳徑范圍是小于0.1r和大于0.3r的槳徑范圍。在小于0.1r的槳徑區(qū)域內(nèi)，葉素處于負(fù)攻角工況，流動分離發(fā)生在槳葉的背風(fēng)面，在大于0.3r的槳徑區(qū)域，流動分離發(fā)生在槳葉的迎風(fēng)面，與前述情況一樣，分離現(xiàn)象隨著槳徑的增大逐漸嚴(yán)重。

圖3(c)是偏角 β=45°、偏角 α=60°時槳葉周圍的流線分布圖。比較圖3(b)和圖3(c)可見，偏角α=60°時槳葉周圍的流動分離現(xiàn)象比偏角α=30°時槳葉周圍的流動分離現(xiàn)象更為嚴(yán)重，首先，r=0.1處葉素背風(fēng)面分離渦的尺度加大，其次，r=0.3處葉素迎風(fēng)面分離渦的尺度加大。但這兩個工況下槳葉上發(fā)生分離的區(qū)域是相同的，分離渦的尺度隨槳徑的變化規(guī)律也是相同的。

圖3(d)是偏角 β=45°、偏角 α=90°時槳葉周圍的流線分布圖。由圖3(d)可見，當(dāng)α=90°時，螺旋槳槳葉周圍的流場非常復(fù)雜，且與α為0°、30°和60°的工況有所不同。該工況下，螺旋槳槳葉背風(fēng)面沒有發(fā)生流動分離現(xiàn)象，但整個槳葉迎風(fēng)面都發(fā)生了強度很大的流動分離，且分離渦的尺度隨著槳徑的增大逐漸減小。

4 結(jié)論

本文基于滑移網(wǎng)格模型，數(shù)值研究了螺旋槳側(cè)飛狀態(tài)時側(cè)偏角對螺旋槳氣動性能的影響。主要結(jié)論有:由于螺旋槳軸向速度隨著偏角的增大逐漸減小，所以螺旋槳軸向拉力隨著偏角的增大逐漸增大。螺旋槳的側(cè)向阻力隨著偏角的增大逐漸減小。槳葉上的分離區(qū)域隨著偏角的增大逐漸增大。當(dāng)偏角為0°、30°、60°時，槳葉上分離渦的尺度隨著槳徑的增大逐漸增大。當(dāng)偏角為90°時，槳葉周圍的流場非常復(fù)雜，且分離渦的尺度隨著槳徑的增大逐漸減小，與其余工況有所不同。

［1］Roosenboom Eric W M ，Heider A，Schr?der A.Propeller Slipstream Development［C］.25th AIAA Applied Aerodynamics Conference，25-28 June 2007，Miami，F(xiàn)L，AIAA 2007-3810.

［2］聶營，王生，楊燕初.螺旋槳靜推力數(shù)值模擬與實驗對比分析［J］.計算機仿真，2009，26(3):103-107.

［3］聶營.平流層飛艇高效螺旋槳設(shè)計與試驗研究［D］.北京:中國科學(xué)院光電研究所碩士學(xué)位論文，2008.

［4］許建華，宋文萍，韓忠華.基于嵌套網(wǎng)格技術(shù)的風(fēng)力機葉片繞流數(shù)值模擬［J］.太陽能學(xué)報，2010，31(1):91-94.

［5］許建華，宋文萍，韓忠華，等.基于CFD技術(shù)的螺旋槳氣動特性研究［J］.航空動力學(xué)報，2010，25(5):1103-1109.

［6］羅淞，楊永，左歲寒.螺旋槳流場Euler/NS數(shù)值解對比分析［J］.航空計算技術(shù)，2011，41(1):74-77.

［7］Yakhot V，Thangam S，Gatski T B，et al.Development of Turbulence Models for Shear Flows by a Double Expansion Technique［R］.NASA Contractor Report 187611，ICASE Report No.91-65.

［8］Rubinstein R.Renormallzation Group Theory of Bolgiano Scaling in Boussinesq Turbulence［R］.NASA Technical Memorandum 106602，ICOMP-94-8;CMOTT-94-2.

［9］程鈺鋒，聶萬勝.基于滑移網(wǎng)格的臨近空間螺旋槳流場數(shù)值仿真［J］.直升機技術(shù)，2012(2).

［10］Colozza A J.APEX 3D Propeller Test Preliminary Design［R］.NACA/CR-2002-211866，2002.

［11］Thompson.J F，Thames F C，Martin C W.Automatic Numerical Generation of Body-Fitted Curvilinear Coordinate System for Field Containing Any Number of Arbitrary Two-Dimensional Bodies［J］.J.Comput.Phys.1974(15):299-399.