垂直下降狀態下的旋翼三維流場數值模擬

曹 棟 曹義華

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

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

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

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

1 控制方程和湍流模型

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

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

2 計算方法

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

2.1 壓差分布計算

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

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

而懸停誘導速度vh可表示為

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

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

根據迎角、當地雷諾數和馬赫數，以及剖面的翼型型號，可以得到翼型的升力系數Cl和阻力系數Cd.本文通過FORTRAN程序計算了不同剖面迎角和下降速度下的升力系數和阻力系數，從而，槳葉剖面上的升力和阻力分別為

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

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

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

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

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

2.2 網格生成與邊界條件

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

圖3 計算域網格劃分

2.3 懸停算例驗證

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

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

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

3 算例計算與分析

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

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

3.1 小速度下降情況

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

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

3.2 渦環及湍流狀態

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

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

圖9 湍流狀態下的中截面處的速度矢量分布

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

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

圖10 給出了 Vc=5，7，9，11m/s時，槳盤下方0.660R處Vz沿槳葉徑向位置的分布，相比于小速度下降狀態，在渦環形成階段，z軸方向速度變化并不明顯，其絕對峰值并未因下降率的增加有所減少甚至有所增加，從－12.5 m/s變為－15 m/s(定義相對氣流方向向上為正，則下洗氣流方向為負)，非定常的氣流特征開始呈顯;當進入完全發展的渦環狀態，對應于Vc=9 m/s時，增加的下降率對Vz的影響趨于明顯，同時在此截面處沿槳葉徑向，速度方向發生改變，與懸停情況類似，在內側速度為正，氣流方向向上，外側速度為負，氣流方向向下.但上洗的范圍有所擴大(r/R≤0.55)，相反地，氣流方向導致空氣動力特性十分復雜，這也與其所對應的速度矢量圖相符合，靠近內側速度矢量向上，靠外側速度矢量向下.而當旋翼處于湍流狀態時(Vc=11 m/s)，Vz沿槳葉徑向變化趨于0且均為正，即在此處位置，向上的自由流已完全抵消了旋翼的下洗作用.而湍流狀態的速度矢量圖顯示，形成的環狀氣流已經發展到旋翼上方，槳盤下方的的流場較為穩定.

圖10 渦環及湍流時槳盤下方0.660R處的Vz分布

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

3.3 風車狀態

穩定的風車狀態(一般Vc≥|1.8vh|)對應于

圖12 文獻［7］利用PIV實驗得到的典型渦環狀態下單側槳葉的速度矢量分布

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

圖13 風車狀態下中截面處的速度矢量分布

圖14 風車狀態時槳盤下方0.660R處的Vz分布

3.4 垂直下降性能預測

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

圖15 旋翼載荷與誘導速度及垂直下降飛行的關系

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

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

圖17 模型旋翼垂直下降狀態的誘導速度分布

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

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

圖19 模型旋翼垂直下降狀態的需用功率分布

4 結論

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

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

References)

［1］普勞蒂.直升機性能及穩定性和操縱性［M］.高正，譯.北京:航空工業出版社，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］于子文，曹義華.前飛旋翼三維湍流場的數值模擬［J］.北京航空航天大學學報，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］.北京航空航天大學學報，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)