垂直下降狀態(tài)下的旋翼三維流場數(shù)值模擬

曹 棟 曹義華

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

眾所周知，直升機(jī)流場的數(shù)值模擬最困難的問題在于旋翼的流場模擬.之前較流行的旋翼流場計算分析主要專注于懸停和前飛，對下降狀態(tài)的旋翼流場模擬較少.因為下降的旋翼，當(dāng)其處于渦環(huán)(vertex ring)狀態(tài)時［1］，由于來流速度朝上，遠(yuǎn)處下游尾跡在旋翼的上方，旋翼誘導(dǎo)速度與相對氣流方向相反，兩股反向氣流相遇形成紊亂的漩渦.同時旋翼的槳尖渦不能被吹離旋翼而聚集在一起，旋翼排出的部分氣流又被吸入，繞旋翼外緣打轉(zhuǎn)，流態(tài)復(fù)雜且很不穩(wěn)定.1982—2003年，3/4的直升機(jī)事故與進(jìn)入渦環(huán)狀態(tài)有關(guān)［2］.

對直升機(jī)旋翼流場的數(shù)值模擬，通常采用兩種方法.一種是致力于模擬旋轉(zhuǎn)槳葉附近的細(xì)致流動，包括捕捉槳葉附近的激波，分析槳尖渦的形成過程等.一般需要在每片槳葉周圍生成貼體網(wǎng)格，將大部分網(wǎng)格點集中于槳葉附近，尤其是槳葉的邊界層內(nèi)，這種方法可以準(zhǔn)確得到槳葉附近的流動特性.不過為了提高遠(yuǎn)場尾跡的計算精度也必須保證遠(yuǎn)場的網(wǎng)格密度，這將導(dǎo)致網(wǎng)格數(shù)目過大，計算代價很高.而采用一些尾跡模型來求解也會存在一些問題，因為尾跡結(jié)構(gòu)會隨著飛行條件和槳葉結(jié)構(gòu)的改變而變化.對于垂直下降飛行，特別是渦環(huán)，湍流狀態(tài)時由于氣流對沖，很難有合適的尾跡模型能很好地描述.

本文采用另一種方法，將旋翼簡化為一無限薄激勵盤，忽略旋翼槳葉附近的細(xì)致流動特征，采用動量葉素理論結(jié)合下降狀態(tài)時的誘導(dǎo)速度經(jīng)驗公式來計算槳葉的載荷分布，通過在N-S(Navier-Stokes)方程中加入此分布的載荷(槳盤壓力差)來代替旋翼對空氣的作用.由于用圍繞整個槳盤的網(wǎng)格取代圍繞槳葉的貼體網(wǎng)格，減小了網(wǎng)格生成的難度和網(wǎng)格數(shù)目，從而有效地節(jié)省了計算時間并能保證一定的計算精度.國外，文獻(xiàn)［3］等采用這種方法對旋翼的流場進(jìn)行了分析，而國內(nèi)文獻(xiàn)［4］等也使用這種思想進(jìn)行了旋翼三維流場的研究，但他們都僅限于懸停和前飛，本文通過對旋翼垂直下降進(jìn)行數(shù)值模擬，對旋翼下降時的流場和性能預(yù)測有一定的參考價值.

1 控制方程和湍流模型

由于垂直飛行時具有軸對稱特性，理論上不同方位角的流動特性相似，且本文采用的模型旋翼的槳尖速度較小且整個旋翼的流場速度也不大，因而將整個流場看作不可壓有粘流動，所以采用非定常、不可壓、湍流的N-S方程為控制方程.并結(jié)合使用S-A(Spalart-Allmaras)一方程的湍流模型，此湍流模型可以很好地適用于翼型、壁面邊界層等流動.控制方程的通用形式如下:

式中，Φ為廣義變量;ρ為空氣密度;V為速度矢量;Γ為廣義擴(kuò)散系數(shù);SΦ為廣義源項.對于特定意義的Φ，具有特定的Γ和SΦ.式(1)中的4項分別是非定常項、對流項、擴(kuò)散項和源項.由于不考慮熱傳導(dǎo)，SΦ僅存在動量方程中，表示單位體積內(nèi)旋翼賦予流場的氣動力，對于本文采用的無限薄作用盤而言，僅表示了單位面積的旋翼賦予流場的拉力(壓力)差.

2 計算方法

本文計算時利用商用軟件直接求解三維粘性不可壓N-S方程，微分方程的離散使用有限體積法，壓力速度耦合迭代選用著名SIMPLER(Semi-Implicit Method for Pressure-Linked Equations Revised)算法，對流項使用二階迎風(fēng)差分格式離散，離散得到的代數(shù)方程組用Gauss-Seidel迭代法求解.而旋翼槳盤的壓力分布則由FORTRAN編程得出，再通過多項式擬合成自定義函數(shù)(UDF)加載到商用軟件中.

2.1 壓差分布計算

從圖1中可以看出，為了求得槳葉翼型剖面上的升力和阻力，需要用到葉素相對速度U，可表示為

式中，UP，UT分別為槳葉剖面的切向速度和軸向速度;Ω，r分別為槳葉旋轉(zhuǎn)角速度和槳葉剖面徑向位置;Vc，vi分別為旋翼槳盤下降速度和與之對應(yīng)的誘導(dǎo)速度，考慮到在一定的下降速度下，經(jīng)典的動量理論不再成立，文獻(xiàn)［5］給出了一種近似的誘導(dǎo)速度分布，vi可表示為

而懸停誘導(dǎo)速度vh可表示為

式中，a，c，Nb，R 分別為升力線斜率、槳葉寬度、槳葉片數(shù)和旋翼半徑;θ為槳距角.

圖1 槳葉翼型剖面上的氣動力

根據(jù)迎角、當(dāng)?shù)乩字Z數(shù)和馬赫數(shù)，以及剖面的翼型型號，可以得到翼型的升力系數(shù)Cl和阻力系數(shù)Cd.本文通過FORTRAN程序計算了不同剖面迎角和下降速度下的升力系數(shù)和阻力系數(shù)，從而，槳葉剖面上的升力和阻力分別為

同時，來流角β和攻角α可分別表示為

則微段拉力dT和微段扭矩dQ可簡單地表示為

然后通過積分求解，可得到整個槳盤的拉力和扭矩，而槳盤上的壓強(qiáng)分布則可表達(dá)為

通過不同的槳盤徑向位置，可以得到不同的ΔP，本文是將槳盤沿徑向分為40個微段，這樣可使后面的曲線擬合獲得足夠精確的結(jié)果，圖2給出了壓差源項的計算流程.對于懸停狀態(tài)下的模型旋翼(槳盤半徑為0.914 m)，擬合后的公式為

圖2 動量源項槳盤壓力差的計算流程圖

2.2 網(wǎng)格生成與邊界條件

為了兼顧計算精度與速度，本文利用GAMBIT前處理軟件生成了六面體結(jié)構(gòu)化網(wǎng)格，計算域為包含旋翼的圓柱體，網(wǎng)格在旋翼附近分布較密，離開旋翼逐漸變稀.垂直飛行時，選取半徑為10R、高為30R的計算域，為了滿足垂直下降的流場分析，槳盤上部的網(wǎng)格計算域高度為下部的兩倍，網(wǎng)格數(shù)共有1 098 540個，見圖3.為了符合模型旋翼附近的不可壓流場狀況，將整個計算域外表面定義為速度進(jìn)口邊界，槳盤定義為風(fēng)扇邊界，其上將加載由前面計算得到的壓力分布源項，而槳轂則設(shè)定為壁面邊界條件.

圖3 計算域網(wǎng)格劃分

2.3 懸停算例驗證

為了能與現(xiàn)有的實驗結(jié)果對照，本文采用了S-A湍流模型計算了模型旋翼槳盤下方0.104R，0.215R，0.660R 及0.993R 位置處的動壓分布，旋翼參數(shù)見文獻(xiàn)［6］.圖 4、圖 5給出了 0.104R，0.660R處的計算結(jié)果與實驗結(jié)果［6］對比.從圖中可以看出，采用本文方法得到的懸停旋翼的下洗流場動壓分布與實驗值的變化情況基本吻合，為后文的垂直下降算例計算奠定了基礎(chǔ).

圖4 槳盤下方0.104R處的動壓分布

圖5 槳盤下方0.660R處的動壓分布

3 算例計算與分析

采用上文中提到的方法，本文模擬不同下降率的旋翼流場，從小下降率1 m/s到高下降率的19 m/s，并繪出了不同下降率下的流場流線特征(見圖6)，流場速度分布和旋翼性能預(yù)測圖(包括拉力和攻角變化)，并分析了產(chǎn)生變化的原因.

圖6 不同下降率下的旋翼流場流線圖

3.1 小速度下降情況

對應(yīng)于圖6中Vc=1 m/s和3 m/s的情況，可以看出由于下降率較小，向上的氣流很難抵消旋翼自身下洗氣流的影響，直到槳盤下方較遠(yuǎn)處下洗速度有所減弱后，兩股相對氣流對沖才卷起較大的渦，但由于其距槳盤較遠(yuǎn)，對旋翼的性能影響不是太大.圖7a給出了Vc=1 m/s時計算域中截面(y=0)處的速度矢量圖，明顯可以看出旋翼的下洗速度仍占主導(dǎo)地位，僅在離槳盤周邊很遠(yuǎn)處才有向上的速度分布.而圖7b則捕捉到了小速率下降時旋翼上方的速度鞍點(速度為0的點，位于正中央)，證明了旋翼在小速率下降的情況下旋翼上方由于流場不連續(xù)所導(dǎo)致的非定常特征.

圖8給出了懸停以及Vc=1 m/s和3 m/s時，槳盤下方0.660R處的垂向速度Vz分布，其中懸停時，槳盤中心區(qū)域(r/R≤0.3)的誘導(dǎo)速度為正(定義相對氣流方向向上為正)，即形成了一個較大的上洗區(qū)域，而隨著下降率的增加，此上洗區(qū)域消失.直到旋翼進(jìn)入渦環(huán)狀態(tài)，見圖9，又出現(xiàn)了部分上洗的情況.

3.2 渦環(huán)及湍流狀態(tài)

圖7 小下降率情況下的中截面(y=0)處的特征分布

圖8 小速度下降時槳盤下方0.660R處的Vz分布

圖9 湍流狀態(tài)下的中截面處的速度矢量分布

一般認(rèn)為當(dāng)Vc=0.7vh時，旋翼進(jìn)入渦環(huán)狀態(tài)，本模型旋翼計算得到的vh≈7 m/s，所以可以認(rèn)為圖6中的Vc=5 m/s已進(jìn)入渦環(huán)，由于下降率的增加，向上的自由流使槳尖渦螺旋線堆積在槳盤的下方，槳尖渦開始呈輻射狀從旋翼槳盤向外發(fā)展，形成了靠近槳盤的環(huán)狀情形.當(dāng)下降率繼續(xù)增加，對應(yīng)于圖6中Vc=7 m/s以及9 m/s的情況，槳尖渦環(huán)有所縮小并逐步靠近槳盤，形成了成熟的渦環(huán)狀態(tài).由于旋翼不斷旋轉(zhuǎn)，渦環(huán)強(qiáng)度不斷積累，甚至出現(xiàn)了左右不對稱的情況，所以整個流場的非定常特性十分明顯.

當(dāng)Vc=11 m/s時，形成的渦環(huán)上移至槳盤上方，旋翼進(jìn)入湍流狀態(tài)，在這種狀態(tài)下，氣流仍然會有高水平的擾動，但是，因為槳盤處的速度向上，所以穿過旋翼的環(huán)流少了很多.旋翼承受著由湍流造成的某些顛簸，但沒有渦環(huán)狀態(tài)時明顯，圖10中的垂向速度分布也證明了這點.圖9、圖11分別給出了Vc=11 m/s和9 m/s時計算域中截面(y=0)處的速度矢量圖，明顯可以看出旋翼向上的自由流速度已占主導(dǎo)地位.而圖11中的典型渦環(huán)狀態(tài)的速度矢量分布與文獻(xiàn)［7］采用PIV實驗得到的結(jié)果有明顯的相似(圖12)，對比圖9、圖11，前者的自由流強(qiáng)度更加明顯，且形成的速度湍流已移至槳盤上方.

圖10 給出了 Vc=5，7，9，11m/s時，槳盤下方0.660R處Vz沿槳葉徑向位置的分布，相比于小速度下降狀態(tài)，在渦環(huán)形成階段，z軸方向速度變化并不明顯，其絕對峰值并未因下降率的增加有所減少甚至有所增加，從－12.5 m/s變?yōu)椋?5 m/s(定義相對氣流方向向上為正，則下洗氣流方向為負(fù))，非定常的氣流特征開始呈顯;當(dāng)進(jìn)入完全發(fā)展的渦環(huán)狀態(tài)，對應(yīng)于Vc=9 m/s時，增加的下降率對Vz的影響趨于明顯，同時在此截面處沿槳葉徑向，速度方向發(fā)生改變，與懸停情況類似，在內(nèi)側(cè)速度為正，氣流方向向上，外側(cè)速度為負(fù)，氣流方向向下.但上洗的范圍有所擴(kuò)大(r/R≤0.55)，相反地，氣流方向?qū)е驴諝鈩恿μ匦允謴?fù)雜，這也與其所對應(yīng)的速度矢量圖相符合，靠近內(nèi)側(cè)速度矢量向上，靠外側(cè)速度矢量向下.而當(dāng)旋翼處于湍流狀態(tài)時(Vc=11 m/s)，Vz沿槳葉徑向變化趨于0且均為正，即在此處位置，向上的自由流已完全抵消了旋翼的下洗作用.而湍流狀態(tài)的速度矢量圖顯示，形成的環(huán)狀氣流已經(jīng)發(fā)展到旋翼上方，槳盤下方的的流場較為穩(wěn)定.

圖10 渦環(huán)及湍流時槳盤下方0.660R處的Vz分布

圖11 典型渦環(huán)狀態(tài)下的中截面(y=0)處的速度矢量分布

3.3 風(fēng)車狀態(tài)

穩(wěn)定的風(fēng)車狀態(tài)(一般Vc≥|1.8vh|)對應(yīng)于

圖12 文獻(xiàn)［7］利用PIV實驗得到的典型渦環(huán)狀態(tài)下單側(cè)槳葉的速度矢量分布

圖6中的Vc=15 m/s，此時氣流從槳盤下方向上穿過槳盤平面，流場又趨于穩(wěn)定，圖13給出了Vc=15 m/s時計算域中截面(y=0)處的速度矢量圖，由于自由來流速度較大，完全抵消了旋翼的下洗流速度，流場中相對氣流方向都已完全向上.而圖14則給出了Vc=13，15，17m/s時，槳盤下方0.660R處Vz沿槳葉徑向位置的分布，與渦環(huán)和紊流狀態(tài)相比，速度值皆為正(即速度方向向上)，且沿槳葉徑向分布區(qū)域平緩，流場變得穩(wěn)定起來.

圖13 風(fēng)車狀態(tài)下中截面處的速度矢量分布

圖14 風(fēng)車狀態(tài)時槳盤下方0.660R處的Vz分布

3.4 垂直下降性能預(yù)測

前面對流場的分析更主要是為了得到旋翼在下降情況下的性能特征，旋翼性能與流場誘導(dǎo)速度相互關(guān)系如圖15所示，可以看到兩者相互影響，存在一定的耦合.

圖15 旋翼載荷與誘導(dǎo)速度及垂直下降飛行的關(guān)系

本文采用兩種方法給出了下降飛行時固定槳距情況下的平均拉力，見圖16，可以看到兩種方法都在完全的渦環(huán)狀態(tài)(或湍流狀態(tài)初始)時出現(xiàn)最大的拉力損失，即 －1.6≤Vc/vh≤ －1.3，恰巧的是，本文計算和實驗［8－9］都顯示，在這個范圍里平均誘導(dǎo)速度也達(dá)到峰值，見圖17.因此不難理解入流的增加正是導(dǎo)致拉力損失的關(guān)鍵，這與文獻(xiàn)［10］的解釋類似.至于在 Vc≥15 m/s后，兩種方法的計算結(jié)果偏差明顯，可能的原因是在編制程序中并沒有涉及槳葉失速［11］的情況，因為按照文獻(xiàn)［5］中的初始誘導(dǎo)速度分布特性，在高下降率的情況下，槳葉剖面攻角早已超過了NACA0012翼型的失速攻角，如圖18所示，所以得到的結(jié)果偏大.

圖16 不同下降率下的槳盤平均拉力

圖17 模型旋翼垂直下降狀態(tài)的誘導(dǎo)速度分布

在得到槳盤平均拉力和誘導(dǎo)速度之后可以很容易計算出模型旋翼的功率，見圖19，與文獻(xiàn)［8］吻合較好.其中Ph為懸停時的功率值，可以看到雖然在渦環(huán)狀態(tài)的情況下，拉力有較大的損失，但由于誘導(dǎo)速度的增加，得到的需用功率值仍要大于懸停的需用功率，而在風(fēng)車狀態(tài)的情況下，特別是 －1.8≤Vc/vh≤ －1.6的范圍里，旋翼功率甚至小于懸停時的需用功率，旋翼向周圍的空氣吸收能量.

圖18 不同下降率下槳葉攻角的徑向分布

圖19 模型旋翼垂直下降狀態(tài)的需用功率分布

4 結(jié)論

本文通過引入壓力源項求解非定常、不可壓、湍流的N-S方程，來模擬直升機(jī)旋翼的垂直下降飛行情況.這種CFD(Computational Fluid Dynamics)分析方法，只需要知道旋翼槳葉的幾何特征和氣動特性參數(shù)，通過FORTRAN程序?qū)⒌玫降膲毫Ψ植家訳DF的形式代入軟件耦合求解就能得到旋翼的全部三維流動信息，比如，不同下降率情況下旋翼附近的流態(tài)特征和某一截面的速度分布.同時計算得到了不同下降率下旋翼誘導(dǎo)速度分布以及拉力和功率分布，計算結(jié)果與文獻(xiàn)吻合較好，并驗證了在渦環(huán)狀態(tài)下，旋翼拉力損失主要是由增加的誘導(dǎo)速度造成的.

此外，本研究還具有一定的實際意義:比如更改網(wǎng)格后可分析不同的旋翼參數(shù)(槳葉負(fù)扭度、槳葉半徑等)對下降飛行時旋翼性能的影響，包括對渦環(huán)發(fā)生邊界的準(zhǔn)確預(yù)測，這對新型旋翼設(shè)計以及如何避免進(jìn)入渦環(huán)狀態(tài)有一定的參考價值.目前正在進(jìn)行一些后續(xù)的工作，如采用準(zhǔn)確槳葉建模結(jié)合滑移網(wǎng)格取代激勵盤模型，可控的槳距變化取代固定槳距，以及機(jī)身對旋翼性能的影響等.

References)

［1］普勞蒂.直升機(jī)性能及穩(wěn)定性和操縱性［M］.高正，譯.北京:航空工業(yè)出版社，1990:77－79 Prouty R W.Helicopter performance，stability and control［M］.Translated by Gao Zheng.Beijing:Aeronautical Industry Press，1990:77－79(in Chinese)

［2］Varnes D J.Development of a helicopter vortex ring statewarning system through a moving map display computer［D］.Monterey，California:Department of Aeronautics and Astronautics，Naval Postgraduate School，1999

［3］Chaffin M S，Berry J D.Helicopter fuselage aerodynamics under a rotor by Navier-Stokes simulation［J］.Journal of the American Helicopter Society，1997，47(3):235 －243

［4］于子文，曹義華.前飛旋翼三維湍流場的數(shù)值模擬［J］.北京航空航天大學(xué)學(xué)報，2006，32(7):751－754 Yu Ziwen，Cao Yihua.Three dimensional turbulence numerical simulation of rotor in forward flight［J］.Journal of Beijing University of Aeronautics and Astronautics，2006，32(7):751 － 754(in Chinese)

［5］Leishman J G.Principle of helicopter aerodynamics［M］.Cambridge:Cambridge University Press，2006:81 －91

［6］McKee J W，Naeseth R L.Experimental investigation of the drag of flat plates and cylinders in the slipstream of a hovering rotor［R］.NACA TN 4239，1958

［7］Green R B，Gillies E A，Brown R E.The flow field around a rotor in axial descent［J］.Journal of Fluid Mechanics，2005(354):237－261

［8］Washizu K，Azuma A，Koo J，et al.Experiments on a model helicopter rotor operating in the vortex ring state［J］.Journal of Aircraft，1966，3(3):225 － 230

［9］Castles W，Jr，Gray R B.Empirical relation between induced velocity，thrust，and rate of descent of a helicopter rotor as determined by wind-tunnel tests on four model rotors［R］.NASA TN 2474，1951

［10］Azuma A，Obata A.Induced flow variation of the helicopter rotor operating in the vortex ring state［J］.Journal of Aircraft，1968，5(4):381－386

［11］蘇媛，曹棟，曹義華.計及大迎角失速的直升機(jī)旋翼的氣動力模擬［J］.北京航空航天大學(xué)學(xué)報，2010，36(2):168－171 Su Yuan，Cao Dong，Cao Yihua.Simulation of helicopter rotor aerodynamic force in conditions of high angle of attack and dynamic stall［J］.Journal of Beijing University of Aeronautics and Astronautics，2010，36(2):168 －171(in Chinese)