程鈺鋒，聶萬勝，胡永平

(1.裝備學(xué)院研究生院，北京懷柔 101416;2.裝備學(xué)院航天裝備系，北京懷柔 101416;3.河北滄州飛行訓(xùn)練基地，河北滄州 061036)

0 引言

現(xiàn)今油價(jià)大幅上漲，節(jié)油成為飛行器設(shè)計(jì)中必須考慮的一個(gè)環(huán)節(jié)［1］。螺旋槳發(fā)動(dòng)機(jī)因其在亞音速范圍內(nèi)的節(jié)油優(yōu)點(diǎn)引起了人們的普遍關(guān)注和高度重視，螺旋槳的研究工作也得到了進(jìn)一步的發(fā)展。螺旋槳?dú)鈩?dòng)性能的研究是螺旋槳?jiǎng)恿ο到y(tǒng)研制工作的關(guān)鍵，國內(nèi)外許多專家學(xué)者都采用計(jì)算流體力學(xué)(computational fluid dynamics，CFD)技術(shù)研究了螺旋槳的氣動(dòng)性能，取得了很多有意義的成果。CFD是近三十年來迅速發(fā)展的技術(shù)，但它已經(jīng)廣泛的滲入到社會(huì)和生活的各個(gè)方面，在研究流動(dòng)現(xiàn)象、解決流體工程問題等方面發(fā)揮了重要作用，尤其在航空航天領(lǐng)域已成為不可或缺的核心技術(shù)之一［2］。

螺旋槳和旋翼流場的復(fù)雜性使得螺旋槳、旋翼CFD技術(shù)總體落后于固定翼CFD技術(shù)［3］。起初，人們采用小擾動(dòng)勢方程法［4］研究螺旋槳的旋轉(zhuǎn)流場，但由于使用了小擾動(dòng)的假設(shè)，只對(duì)于薄翼等擾動(dòng)不太強(qiáng)的跨音速流動(dòng)才能給出較好的結(jié)果。后來，有人采用全勢方程法［5-7］研究螺旋槳的旋轉(zhuǎn)流場，全勢方程求解的復(fù)雜程度介于歐拉方程和小擾動(dòng)方程之間，在激度不太強(qiáng)的情況下，具有較好的模擬精度。為了更準(zhǔn)確地模擬旋翼流場中出現(xiàn)的激波和旋渦流動(dòng)，從20世紀(jì)80年代以來，許多研究以歐拉方程和N-S方程來求解旋翼流場［8］。我國的研究工作始于20世紀(jì)90年代。1993年，西北工業(yè)大學(xué)的楊國偉和何植岱［9］采用渦格升力面法，建立一種三維自由尾渦模型，構(gòu)造一種松弛迭代法，計(jì)算得到了收縮尾渦。1996年，西北工業(yè)大學(xué)的王立群、喬志德等人［10-11］，采用網(wǎng)格中心有限體積法和五步Ruuge-Kutta(5，2)時(shí)間推進(jìn)格式，通過求解三維歐拉方程，仿真得到了直升機(jī)旋翼懸停流場，在沒有任何尾跡模型的情況下，計(jì)算了兩旋翼在亞聲速和跨聲速時(shí)的壓力分布。1998年，上海交通大學(xué)的杜朝暉［12-13］為解決水平軸風(fēng)力渦輪設(shè)計(jì)與性能預(yù)估方法存在失速延遲現(xiàn)象的問題，以旋轉(zhuǎn)葉輪三維邊界層方程為理論基礎(chǔ)，在分析了風(fēng)力渦輪葉片表面產(chǎn)生失速延遲的機(jī)理后，建立了風(fēng)力渦輪設(shè)計(jì)與評(píng)估方法的三維失速延遲修正模型，提高了水平軸風(fēng)力渦輪的設(shè)計(jì)水平。20世紀(jì)初，國內(nèi)關(guān)于螺旋槳?dú)鈩?dòng)仿真的研究較少。2008年，中科院的聶營、王生等人［14-15］采用Gambit軟件對(duì)螺旋槳進(jìn)行幾何建模，再用Fluent軟件基于滑移網(wǎng)格，仿真研究了臨近空間螺旋槳的氣動(dòng)性能。2008年，西工大的許建華、宋文萍等人［16］采用雷諾平均N-S方程和嵌套網(wǎng)格技術(shù)對(duì)美國國家可再生能源實(shí)驗(yàn)室的CER實(shí)驗(yàn)型風(fēng)力機(jī)葉片進(jìn)行了仿真研究，湍流方程是B-L代數(shù)模型;計(jì)算中采用了有限體積空間離散法和改進(jìn)型隱式LU-SGS時(shí)間推進(jìn)格式，所用的嵌套網(wǎng)格技術(shù)有效的捕捉了螺旋槳的尾渦;在此基礎(chǔ)上，2009年［17］，他們采用在不同粗細(xì)網(wǎng)格上消除不同頻率誤差加速解的收斂的多重網(wǎng)格技術(shù)，仿真研究了螺旋槳側(cè)流粘性流場，其中多重網(wǎng)格采用非線性方程的全近似格式(FAS)。2010年，西工大的羅淞、楊永等人［18］采用多塊點(diǎn)對(duì)點(diǎn)網(wǎng)格生成技術(shù)，生成了螺旋槳的空間計(jì)算網(wǎng)格，分別研究了歐拉方程組和N-S方程組對(duì)不同進(jìn)距比下螺旋槳運(yùn)動(dòng)過程的仿真結(jié)果，與實(shí)驗(yàn)進(jìn)行對(duì)比后得出歐拉方程比較適合于螺旋槳的工程設(shè)計(jì)。

臨近空間低速飛行器大多采用螺旋槳作為其動(dòng)力裝置，螺旋槳?dú)鈩?dòng)性能的優(yōu)劣決定臨近空間低速飛行器的性能，因此，研究臨近空間環(huán)境下螺旋槳的氣動(dòng)性能是十分迫切的工作，但臨近空間環(huán)境特殊，很難開展實(shí)驗(yàn)研究［19］?；诖?，本文采用滑移網(wǎng)格模型，考慮RNG k-ε湍流模型，通過求解三維非定常N-S方程，在驗(yàn)證了數(shù)學(xué)模型可行性的基礎(chǔ)上，仿真研究了臨近空間螺旋槳非定常旋轉(zhuǎn)流場。仿真結(jié)果與文獻(xiàn)資料和理論分析一致，可以用于臨近空間螺旋槳的設(shè)計(jì)和研究。

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

1.1 控制方程

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

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

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

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

1.2 湍流模型

RNG k-ε湍流模型是基于k-ε標(biāo)準(zhǔn)兩方程的湍流模型，采用一種叫做重正規(guī)化群的數(shù)學(xué)方法對(duì)N-S方程進(jìn)行暫態(tài)推理得到的改進(jìn)型k-ε兩方程湍流模型。它由V.Yakhot和S.A.Orszag于1986年提出并逐步完善的［20-21］，其基本思想是認(rèn)為在流場中小渦是各項(xiàng)同性的，處于統(tǒng)計(jì)定常的和統(tǒng)計(jì)平衡的狀態(tài)。忽略浮力湍動(dòng)能的RNG k-ε湍流模型的輸運(yùn)方程如下:

其中:k是湍流動(dòng)能，ε是湍流耗散率;ui是速度分量，xi是坐標(biāo)分量;αk=αε=1.393 分別是 Prandtl數(shù)對(duì)k和ε的反饋?zhàn)饔孟禂?shù);ueff是有效粘性系數(shù)，Gk是由平均速度梯度引起的湍動(dòng)能;YM是由于可壓縮湍流脈動(dòng)膨脹對(duì)總的耗散率的影響;C1ε=1.42、C2ε=1.68是經(jīng)驗(yàn)常數(shù);Rε是湍流模型中數(shù)的解析項(xiàng)。

式中:Cμ=0.0845;η0=4.38;β =0.012;η =Sk/ε，S是漩渦大小。

由上可知，RNG k-ε湍流模型考慮了低雷諾數(shù)流動(dòng)粘性，改進(jìn)了標(biāo)準(zhǔn)k-ε模型的高雷諾數(shù)性質(zhì)，并且提供了Prandtl數(shù)的解析公式，考慮了湍流漩渦，因此更加適合于雷諾數(shù)不是很高和帶有強(qiáng)漩渦運(yùn)動(dòng)狀態(tài)的數(shù)值仿真。

1.3 計(jì)算方法

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

采用二階精度的有限體積AUSM(Advection Upstream Splitting Method)離散格式對(duì)粘性流體的控制方程和湍流方程進(jìn)行空間離散。AUSM格式是20世紀(jì)90年代Liou和Stefen提出并完善的高分辨率迎風(fēng)格式，融合了FVS穩(wěn)定性好的優(yōu)點(diǎn)和FDS高分辨率的優(yōu)點(diǎn)，具有良好的數(shù)值穩(wěn)定性和較高的間斷分辨率，其基本思想是認(rèn)為對(duì)流波的傳播與聲波的傳播是物理上不同的過程，前者與特征速度u為線性關(guān)系，后者與特征速度u+a和u-a有非線性關(guān)系，將無粘通量分解為對(duì)流通量和壓力通量。詳見文獻(xiàn)［22］。

圖1 耦合算法流程圖

1.4 滑移網(wǎng)格模型

滑移網(wǎng)格是在動(dòng)參考系模型和混合面法的基礎(chǔ)上發(fā)展起來的，常用于風(fēng)車、轉(zhuǎn)子、螺旋槳等運(yùn)動(dòng)的仿真研究。在滑動(dòng)網(wǎng)格模型計(jì)算中，流場中至少存在兩個(gè)網(wǎng)格區(qū)域，每一個(gè)區(qū)域都必須有一個(gè)網(wǎng)格界面與其他區(qū)域連接在一起。網(wǎng)格區(qū)域之間沿界面做相對(duì)運(yùn)動(dòng)。在選取網(wǎng)格界面時(shí)，必須保證界面兩側(cè)都是流體區(qū)域。

滑動(dòng)網(wǎng)格模型允許相鄰網(wǎng)格間發(fā)生相對(duì)運(yùn)動(dòng)，而且網(wǎng)格界面上的節(jié)點(diǎn)無需對(duì)齊，即網(wǎng)格交界面是非正則的。在使用滑動(dòng)網(wǎng)格模型時(shí)，計(jì)算網(wǎng)格界面上的通量需要考慮到相鄰網(wǎng)格間的相對(duì)運(yùn)動(dòng)，以及由運(yùn)動(dòng)形成的重疊區(qū)域的變化過程。

兩個(gè)網(wǎng)格界面相互重合部分形成的區(qū)域稱為內(nèi)部區(qū)域，即兩側(cè)均為流體的區(qū)域，而不重合的部分則稱為“壁面”區(qū)域(如果流場是周期性流場，則不重合的部分則稱為周期區(qū)域)。在實(shí)際計(jì)算過程中，每迭代一次就需要重新確定一次網(wǎng)格界面的重疊區(qū)域，流場變量穿過界面的通量是用內(nèi)部區(qū)域計(jì)算的，而不是用交界面上的網(wǎng)格計(jì)算。

下面，通過一個(gè)簡單的例子說明滑移網(wǎng)格是如何計(jì)算界面信息的。圖4是二維網(wǎng)格分界面示意圖，界面區(qū)域由面A-B、B-C、D-E和面E-F構(gòu)成。交界區(qū)域可以分為a-d、d-b、b-e等。處于兩個(gè)區(qū)域重合部分的面為d-b、b-e和e-c，構(gòu)成內(nèi)部區(qū)域，其他的面(a-d、c-f)則為成對(duì)的壁面區(qū)域。如果要計(jì)算穿過區(qū)域IV的流量，用面d-b和面b-e代替面D-E，并分別計(jì)算從I和III流入IV的流量。

圖2 二維網(wǎng)格分界面示意圖

1.5 模型驗(yàn)證

［15］，通過比較螺旋槳的靜拉力，驗(yàn)證上述數(shù)學(xué)模型的可行性。螺旋槳的幾何建模及網(wǎng)格劃分在gambit環(huán)境下完成，直徑為28 in，螺距為10 in，具體尺寸詳見文獻(xiàn)［15］。

表1是文獻(xiàn)［15］的實(shí)驗(yàn)和計(jì)算結(jié)果與本文計(jì)算結(jié)果的比較。由表1可見，文獻(xiàn)［15］的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的最大誤差為20%，本文計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的最大誤差為15.8%，本文計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)更為接近，說明本文所用計(jì)算模型是可行的。

表1 靜拉力的比較

2 螺旋槳模型及網(wǎng)格

圖3是本文所研究的螺旋槳局部及其網(wǎng)格示意圖，螺旋槳葉素選擇Eppler 387翼型，不同槳徑處的葉素弦長和安裝角都相同且分別為0.05m和20°，螺旋槳直徑D為0.8 m，兩個(gè)槳葉之間的輪轂是長0.06 m、直徑0.03 m的圓柱體。計(jì)算區(qū)域是一個(gè)長10 D、直徑8D的圓柱體。速度入口距螺旋槳4D，給定氣流速度及總溫，壓力出口距螺旋槳6D，給定總溫和總壓;遠(yuǎn)場距螺旋槳轉(zhuǎn)軸4D，給定氣流速度、總壓及總溫，如圖4所示。

由于滑移網(wǎng)格模型允許相鄰網(wǎng)格之間發(fā)生相對(duì)運(yùn)動(dòng)，而且網(wǎng)格界面上的點(diǎn)無需對(duì)齊，即網(wǎng)格是非正則的。利用這一特點(diǎn)，可以更好的分布網(wǎng)格的疏密度，既保證了計(jì)算流場所需要網(wǎng)格數(shù)又使網(wǎng)格總數(shù)減小，從而節(jié)約計(jì)算資源。對(duì)螺旋槳的仿真而言，螺旋槳近區(qū)域流場變化劇烈，所以應(yīng)該加密網(wǎng)格，遠(yuǎn)離螺旋槳的區(qū)域流場較為平緩，對(duì)網(wǎng)格數(shù)目要求不高，所以網(wǎng)格數(shù)目較少。

圖3 螺旋槳局部網(wǎng)格示意圖

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

圖4是整個(gè)計(jì)算區(qū)域的網(wǎng)格示意圖。整個(gè)計(jì)算區(qū)域共分為3個(gè)小區(qū)域。旋轉(zhuǎn)區(qū)域是中間包圍螺旋槳的一個(gè)小區(qū)域，該區(qū)域網(wǎng)格為非結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)約為60萬。旋轉(zhuǎn)區(qū)前后的兩個(gè)區(qū)域網(wǎng)格為結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)約為12萬。

對(duì)螺旋槳的仿真而言，螺旋槳近區(qū)域流場變化劇烈，所以應(yīng)該加密網(wǎng)格，遠(yuǎn)離螺旋槳的區(qū)域流場較為平緩，對(duì)網(wǎng)格數(shù)目要求不高，所以網(wǎng)格數(shù)目較少。基于此，本文將計(jì)算區(qū)域分為氣流入口區(qū)域、旋轉(zhuǎn)區(qū)域和氣流出口區(qū)域3個(gè)小的計(jì)算區(qū)域。螺旋槳前后的氣流入口和出口區(qū)域采用TTM［23］方法生成結(jié)構(gòu)網(wǎng)格，并在各自靠近螺旋槳的一端加密，網(wǎng)格總數(shù)約為12萬。旋轉(zhuǎn)區(qū)域采用TGrid網(wǎng)格劃分法生成非結(jié)構(gòu)網(wǎng)格，螺旋槳表面網(wǎng)格節(jié)點(diǎn)之間距離為1 mm，網(wǎng)格總數(shù)約為60萬。

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

本節(jié)在飛行高度為20 km，螺旋槳前進(jìn)速度分別為5 m/s和20 m/s的條件下，對(duì)不同轉(zhuǎn)速下螺旋槳的運(yùn)動(dòng)過程進(jìn)行仿真，分析螺旋槳拉力、效率隨轉(zhuǎn)速的變化規(guī)律及其原因，研究螺旋槳速度、渦流強(qiáng)度等參數(shù)在旋轉(zhuǎn)流場中的分布規(guī)律。

圖5是前進(jìn)速度為5 m/s和20 m/s時(shí)，螺旋槳的拉力和效率隨其轉(zhuǎn)速的變化規(guī)律比較圖。其中螺旋槳的效率η=TV0/2pinsM，T是螺旋槳總拉力，M是總扭矩，ns是轉(zhuǎn)速，V0是螺旋槳的前進(jìn)速度。由圖可見，當(dāng)螺旋槳的前進(jìn)速度不相同時(shí)，螺旋槳的氣動(dòng)性能隨轉(zhuǎn)速變化而變化的規(guī)律相同，但總的氣動(dòng)性能如拉力系數(shù)和效率還是有所差別，前進(jìn)速度為20 m/s時(shí)拉力和效率比前進(jìn)速度為5 m/s時(shí)稍大。由圖5(a)可見，在同一個(gè)前進(jìn)速度條件下，螺旋槳拉力隨轉(zhuǎn)速的減小而減小;當(dāng)轉(zhuǎn)速減小到一定的程度后，螺旋槳拉力變?yōu)樨?fù)值，這是因?yàn)槁菪龢墓ぷ鳡顟B(tài)已經(jīng)由前進(jìn)狀態(tài)轉(zhuǎn)為制動(dòng)狀態(tài)，此時(shí)螺旋槳不產(chǎn)生拉力，而產(chǎn)生與飛行方向相反的阻力。由圖5(b)可見，螺旋槳的氣動(dòng)效率隨轉(zhuǎn)速的增大先增大后減小，與文獻(xiàn)［19］的分析結(jié)果一致。

圖5 `拉力和效率隨進(jìn)距比的變化

圖6是前進(jìn)速度為5 m/s，轉(zhuǎn)速分別為60 rpm、180 rpm和1200 rpm時(shí)部分葉素周圍壓力系數(shù)分布比較圖，其中壓力系數(shù)Cp=2(p-p0)/ρV02，p 是流場壓力，p0是環(huán)境壓力。

圖6 前進(jìn)速度為5 m/s時(shí)葉素拉力系數(shù)分布比較圖

由圖6(a)可見，當(dāng)前進(jìn)速度為5 m/s，轉(zhuǎn)速為60 rpm時(shí)，葉素吸力面壓力系數(shù)大于升力面壓力系數(shù)，即槳葉迎風(fēng)面的壓力系數(shù)大于背風(fēng)面壓力系數(shù)，說明該工況下槳葉產(chǎn)生與飛行方向相反的阻力，螺旋槳處于制動(dòng)狀態(tài)。但壓力面和升力面的壓力系數(shù)差別不是很大，因此此時(shí)所產(chǎn)生的負(fù)拉力比較小，與圖5(a)所示一致。

由圖6(b)可見，當(dāng)轉(zhuǎn)速增大到180 rpm時(shí)，螺旋槳的工作狀態(tài)已經(jīng)由制動(dòng)狀態(tài)轉(zhuǎn)向了前進(jìn)狀態(tài)，葉素吸力面壓力系數(shù)小于升力面壓力系數(shù)，即槳葉背風(fēng)面壓力系數(shù)大于迎風(fēng)面壓力系數(shù)，螺旋槳產(chǎn)生與飛行方向一致的正拉力;在槳根部位葉素升力面和吸力面的壓力系數(shù)相差不大，隨著槳徑的延伸，在槳尖部位，葉素升力面壓力系數(shù)與吸力面壓力系數(shù)之間的差別越來越大，這是由于合成氣流的速度大小和攻角都隨著槳徑的延伸而逐漸增大的原因。

由圖6(c)可見，當(dāng)轉(zhuǎn)速繼續(xù)增大到1200 rpm時(shí)，槳葉迎風(fēng)面與背風(fēng)面的壓力系數(shù)之差越來越大，但槳根部位葉素升力面和吸力面壓力之差依然很小;隨著槳徑的延伸，葉素升力面和吸力面壓力之差逐漸增大，當(dāng)槳徑延伸到一定程度后，葉素升力面和吸力面壓力之差又逐漸變小，這可能是由于轉(zhuǎn)速增大，在槳尖部位存在誘導(dǎo)激波的原因。

比較圖6(a)、(b)、(c)可知，當(dāng)前進(jìn)速度不變，螺旋槳轉(zhuǎn)速增大時(shí)，葉素升力面和吸力面壓力之差越來越大，說明槳葉產(chǎn)生的拉力隨著轉(zhuǎn)速的增大逐漸增大。由葉素三角形理論［19］可知，當(dāng)螺旋槳前進(jìn)速度不變且轉(zhuǎn)速增大時(shí)，合成氣流相對(duì)于葉素的速度和攻角逐漸增大，因此槳葉產(chǎn)生的拉力逐漸增大，仿真結(jié)果與圖5(a)和理論分析一致。

圖7是螺旋槳前進(jìn)速度為20 m/s，轉(zhuǎn)速分別為300 rpm、900 rpm和1500 rpm時(shí)，部分葉素周圍流線分布比較圖。由圖7(a)可見，當(dāng)轉(zhuǎn)速為300 m/s時(shí)，從槳根到z=0.1 m處的槳徑區(qū)域，葉素處于負(fù)攻角狀態(tài)，說明有一部分葉素產(chǎn)生負(fù)拉力，這將減小螺旋槳的總拉力，從而降低螺旋槳?dú)鈩?dòng)性能，導(dǎo)致氣動(dòng)效率降低。由圖7(b)可見，當(dāng)轉(zhuǎn)速增大為900 rpm時(shí)，螺旋槳槳根部位流動(dòng)非常穩(wěn)定，但從z=0.2 m處，槳葉表面開始出現(xiàn)流動(dòng)分離現(xiàn)象，當(dāng)槳徑繼續(xù)延伸時(shí)，分離現(xiàn)象逐漸加劇。流動(dòng)分離的存在將會(huì)大大降低螺旋槳的氣動(dòng)性能，尤其是降低螺旋槳的氣動(dòng)效率。由圖7(c)可見，當(dāng)螺旋槳轉(zhuǎn)速繼續(xù)增大到1500 rpm時(shí)，在z=0.1 m的槳徑部位，葉素周圍就開始出現(xiàn)了流動(dòng)分離現(xiàn)象，而且分離渦隨著槳徑的延伸迅速增大，分離現(xiàn)象比轉(zhuǎn)速為900 rpm時(shí)嚴(yán)重很多，說明此時(shí)螺旋槳的氣動(dòng)性能低于轉(zhuǎn)速為900 rpm時(shí)。因此，螺旋槳效率隨轉(zhuǎn)速先提高后降低的原因是:當(dāng)轉(zhuǎn)速很小的時(shí)候，由于槳葉部分葉素處于負(fù)攻角狀態(tài)導(dǎo)致效率降低，當(dāng)轉(zhuǎn)速很大時(shí)，由于槳葉部分葉素處于大迎角工作狀態(tài)，葉素表面流動(dòng)分離嚴(yán)重，導(dǎo)致螺旋槳?dú)鈩?dòng)效率降低。

圖7 前進(jìn)速度20 m/s時(shí)葉素周圍流線分布比較圖

圖8是前進(jìn)速度為5 m/s，轉(zhuǎn)速分別為60 rpm、180 rpm和1200 rpm時(shí)，x=0的橫截面上，流場速度和渦流強(qiáng)度分布比較圖，其中速度單位為m/s，渦流強(qiáng)度的單位為1/s。由圖8(a)可見，當(dāng)轉(zhuǎn)速為60 rpm時(shí)，流場的速度很小，槳尖部位的氣流速度甚至小于槳根部位的速度，且槳根部位速度場的分布范圍也大于槳尖部位;渦流強(qiáng)度與速度分布一樣，槳根部位的渦流強(qiáng)度大于槳尖部位。由圖8(b)可知，當(dāng)轉(zhuǎn)速增大到180 rpm時(shí)，槳尖部位的速度大于槳根部位，槳尖部位速度場的分布范圍也大于槳根部位;渦流強(qiáng)度與速度分布一樣。由圖8(c)可見，當(dāng)轉(zhuǎn)速繼續(xù)增大到1200 rpm時(shí)，槳尖部位的速度進(jìn)一步增大，速度場的分布范圍也在增大，兩個(gè)槳葉速度場之間逐漸開始相互影響;渦流強(qiáng)度大小及分布與速度分布相同。比較圖8(a)、(b)、(c)可知，槳葉周圍速度及速度場的分布范圍，渦流強(qiáng)度及渦流強(qiáng)度分布范圍都隨轉(zhuǎn)速的增大而增大。

圖8 前進(jìn)速度為5 m/s時(shí)x=0處速度和渦量分布比較圖

圖9是前進(jìn)速度5 m/s、轉(zhuǎn)速1200 rpm時(shí)，流場速度分布沿氣流方向的演化圖，速度單位是m/s。螺旋槳中心位置在x=0的地方，因此x小于0的橫截面在螺旋槳前進(jìn)方向，x大于0的橫截面在螺旋槳的后面。由圖可見，在螺旋槳的前面，由螺旋槳槳葉轉(zhuǎn)速運(yùn)動(dòng)產(chǎn)生旋轉(zhuǎn)誘導(dǎo)速度的影響范圍較小，在x=-0.2的地方速度場的分布范圍就已經(jīng)比較均勻。當(dāng)x逐漸增大至靠近螺旋槳，即靠近x=0的橫截面時(shí)，速度場的分布變得越來越不均勻，由槳葉旋轉(zhuǎn)運(yùn)動(dòng)產(chǎn)生的旋轉(zhuǎn)流場對(duì)氣流的影響越來越大，并且逐漸變成由旋轉(zhuǎn)誘導(dǎo)速度為主。在x=0的橫截面上誘導(dǎo)速度的影響范圍最大，最大誘導(dǎo)速度在螺旋槳槳尖部位。當(dāng)x繼續(xù)增大時(shí)，旋轉(zhuǎn)誘導(dǎo)速度對(duì)氣流的影響又逐漸減小，至x=0.7 m處旋轉(zhuǎn)誘導(dǎo)速度的影響不明顯，可見在螺旋槳的后面，旋轉(zhuǎn)誘導(dǎo)速度的影響范圍遠(yuǎn)大于前面的影響范圍。

4 結(jié)論

本文基于滑移網(wǎng)格模型，考慮RNG k-ε湍流模型，通過求解三維非定常N-S方程，在驗(yàn)證了數(shù)學(xué)模型可行性的基礎(chǔ)上，仿真研究了螺旋槳非定常旋轉(zhuǎn)流場。所得流場清晰，速度、流線、渦流強(qiáng)度等流場參數(shù)分布正確，仿真結(jié)果與文獻(xiàn)資料和理論分析一致，可以用于臨近空間螺旋槳的設(shè)計(jì)和研究。主要結(jié)論有:

1)當(dāng)螺旋槳前進(jìn)速度不變，轉(zhuǎn)速增大時(shí)，合成氣流相對(duì)于葉素的速度大小和攻角逐漸增大，所以螺旋槳總拉力隨轉(zhuǎn)速的增大而增大。

2)當(dāng)螺旋槳前進(jìn)速度不變時(shí)，轉(zhuǎn)速很小則有一部分槳葉葉素處于負(fù)迎角狀態(tài)，導(dǎo)致效率降低;轉(zhuǎn)速很大則有一部分槳葉葉素處于大迎角狀態(tài)，葉素表面流動(dòng)分離嚴(yán)重，導(dǎo)致效率降低;所以螺旋槳的效率隨轉(zhuǎn)速的增大先提高后降低。

3)槳葉周圍速度大小及速度場的分布范圍、渦流強(qiáng)度大小及渦流強(qiáng)度分布范圍都隨轉(zhuǎn)速的增大而增大。

4)槳葉旋轉(zhuǎn)運(yùn)動(dòng)產(chǎn)生的誘導(dǎo)速度對(duì)螺旋槳前面流場的影響范圍較小，對(duì)螺旋槳后面流場的影響范圍較大。

參考文獻(xiàn):

［1］Arne S.Unsteady CFD Simulations of Contra-Rotating Propeller Propulsion Systems［C］//Hartford，44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference＆Exhibit.2008，AIAA 2008-5218.

［2］閻超，于劍，許晶磊，等.CFD模擬方法的發(fā)展成就與展望［J］.力學(xué)進(jìn)展，2011，41(5):562-589.

［3］徐國華，招啟軍.直升機(jī)旋翼計(jì)算流體力學(xué)的研究進(jìn)展［J］.南京航空航天大學(xué)學(xué)報(bào)，2003，35(3):338-344.

［4］Caradonna F X，Dessopper A，Tung C.Finite difference modeling of rotor flow including wake effects［J］.Journal of the American Helicopter Society，1984，29(2):86-33.

［5］Chang I C.Transonic flow analysis for rotors，Part 1:Three-dimensional quasi-steady full-potential calcula-tion［R］.NASA TP-2375，1984.

［6］Egolf T A，Sparks S P.A full potential rotor analysis with wake influence using all inner-outer domain techmque［C］.//In the 42th AHS Annual Forum，1986.

［7］Steinhoff J，Ramachandran K.Free wake analysis of compressible rotor ftows［J］.AIAA Journal，1990，88(3):426-431.

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

［9］楊國偉，何植岱.計(jì)及尾渦收縮的螺旋槳滑流計(jì)算［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1995，13(1):83-86.

［10］王立群，喬志德，鐘伯文.直升機(jī)旋翼懸停流場的歐拉方程計(jì)算［J］.空氣動(dòng)力學(xué)報(bào)，1998年，16(3):282-287.

［11］王立群，喬志德，張茹.直升機(jī)旋翼懸停流場的歐拉方程計(jì)算［J］.西北工業(yè)大學(xué)學(xué)報(bào)，1998年，16(4):594-598.

［12］杜朝暉.水平軸風(fēng)力渦輪設(shè)計(jì)與性能預(yù)估方法的三維失速延遲模型 -Ⅰ理論基礎(chǔ)［J］.太陽能學(xué)報(bào)，1999，20(4):392-397.

［13］杜朝暉.水平軸風(fēng)力渦輪設(shè)計(jì)與性能預(yù)估方法的三維失速延遲模型-Ⅱ模型建立及應(yīng)用［J］.太陽能學(xué)報(bào)，2000，21(1):1-5.

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

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

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

［17］許建華，宋文萍，韓忠華，楊旭東，基于CFD技術(shù)的螺旋槳?dú)鈩?dòng)特性研究［J］.航空動(dòng)力學(xué)報(bào)，2010年，第25卷第5期，1103-1109.

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

［19］程鈺鋒，聶萬勝.臨近空間螺旋槳?dú)鈩?dòng)性能分析［J］.直升機(jī)技術(shù)，2011，168:8-13.

［20］Yakhot V，Thangam S，Gatski T B，Orszag S A，Speziale C G.Development of Turbulence Models for Shear Flows by a Double Expansion Technique［R］.NASA Contractor Report 187611，ICASE Report No.91-65.

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

［22］Scandaliato A L.AUSM-based high-order solution for euler equations［C］//48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition 4-7 January 2010，Orlando，F(xiàn)lorida，AIAA 2010-719.

［23］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.