王 娜, 葉 靚, 戚姝妮(中國航空工業(yè)空氣動力研究院, 遼寧 沈陽 110034)

0 引 言

傾轉(zhuǎn)旋翼飛行器由于特殊的構(gòu)型和工作條件，流場結(jié)構(gòu)和氣動力特征復雜。如懸停工作狀態(tài)時，旋翼尾跡受到機翼翼面阻擋，在機翼下方誘導出大范圍的分離，在機翼上表面中線附近匯集形成噴泉流動；前飛狀態(tài)旋翼脫出尾跡與機翼相互干擾，此時旋翼作用類似于螺旋槳；前飛和懸停的過渡狀態(tài)時，旋翼尾跡彎曲嚴重，非定常氣動現(xiàn)象顯著。

有關(guān)傾轉(zhuǎn)旋翼的計算研究工作已經(jīng)開展不少，計算方法包括自由尾跡分析[1-4]和求解Navier-Stokes方程計算[5-7]等。在采用CFD方法時，各種層次的模型也被廣泛使用。由于旋翼流場結(jié)構(gòu)復雜，因而相應計算資源要求也通常較高，選擇合適的計算方法，使得計算資源在可承受的范圍之內(nèi)，又能夠有效描述氣動現(xiàn)象，是一項很重要的課題。求解雷諾平均Navier-Stokes(Reynolds Averaged Navier-Stokes，RANS)方程是被廣泛采用的方法，也有一些研究者使用高計

算量的大渦模擬(Large Eddy Simulation，LES)[8]方法進行氣動干擾分析，但在工作中，他們把旋翼對流場的作用進行了簡化，即采用動量源模型模擬旋翼，使得整體計算的網(wǎng)格數(shù)量得到有效控制，這樣雖然更高精度地模擬了旋翼下洗流通過機翼時的干擾效應，但旋翼周期運動及槳葉外形影響的敏感性帶來的影響被削弱；由于計算量的考慮，一些研究者在旋轉(zhuǎn)坐標系下求解旋翼附近流場區(qū)域[9]，大大節(jié)約了計算資源，但由于忽略了非定常效應影響，且慣性系與非慣性系區(qū)域間數(shù)據(jù)交換通常會產(chǎn)生尾跡截斷，導致流場數(shù)值模擬的精確性下降；一些研究者僅考慮一側(cè)旋翼、短艙及半機身的模型，附加了對稱面進行計算[10]以削減計算量，這樣對于雙旋翼間相互作用的動態(tài)時間影響效果則被忽略。近年來，有關(guān)傾轉(zhuǎn)旋翼的數(shù)值模擬在國內(nèi)成為熱點，研究者們開展的典型研究工作如過渡態(tài)計算方法研究[11-12]和懸停狀態(tài)的氣動干擾分析等[13]。

本文主要對傾轉(zhuǎn)旋翼的懸停狀態(tài)進行數(shù)值模擬研究。與先前研究工作不同的是，本文對采用的數(shù)值計算方法、計算模型和網(wǎng)格系統(tǒng)進行了綜合考慮。首先計算模型真實模擬整個機翼和雙旋翼的所有槳葉外形，由于不采用對稱面、動量源條件并在地面坐標系下求解，所有氣動部件的運動及相互影響被更真實刻畫；其次計算網(wǎng)格方面采用了可自適應的直角背景網(wǎng)格，并預先剖分尾跡區(qū)域，基于流場變量的周期記錄統(tǒng)計，對高物理量梯度區(qū)域(特別是旋翼尾跡區(qū)域、噴泉流動區(qū)及機翼下方的流動分離區(qū)域)實現(xiàn)自動網(wǎng)格細分，因而網(wǎng)格數(shù)量得到有效控制；最后在控制方程求解方面，采用了DES方法進行計算，對比LES方法對網(wǎng)格數(shù)量及計算量的較高要求，DES方法并沒有比RANS方法增加過多額外的計算負擔，且其計算近壁面分離流動的能力更好，是一種當前可行的研究手段。

針對典型工作狀態(tài)，采用以上方案，對傾轉(zhuǎn)旋翼雙旋翼/機翼干擾模型流場和氣動力進行了數(shù)值計算，特別是對槳葉/機翼氣動干擾、機翼前后緣下方流動分離現(xiàn)象和機翼向下載荷變化等進行了定量分析。通過與相同狀態(tài)RANS計算結(jié)果的比較，分析了采用DES模型時，在模擬旋翼、機翼氣動載荷及機翼下方渦演化等流場方面的計算結(jié)果差異。

1 模型、計算狀態(tài)和網(wǎng)格

由于很難查閱到國外傾轉(zhuǎn)旋翼飛行器外形的精確數(shù)據(jù)，計算模型選用了雙旋翼/機翼的干擾類比模型。其中每個旋翼由三片槳葉組成，槳葉半徑R為5.8 m，有較大和不規(guī)則的扭轉(zhuǎn)分布，機翼選擇有上反和小前掠的平面矩形機翼模型，半翼展長L為6.2 m。

計算狀態(tài)參考典型傾轉(zhuǎn)旋翼工作情況，選擇旋翼槳尖馬赫數(shù)0.7。槳葉總距角7.5°，選擇此總距角的原因在于本文主要討論不同方法計算的旋翼/機翼間氣動干擾差異，此時，旋翼表面氣動分離情況不嚴重，不同方法計算的旋翼氣動力差異不會對旋翼/機翼間作用討論帶來更大的干擾影響。

嵌套網(wǎng)格系統(tǒng)由八塊網(wǎng)格組成，即分別圍繞六片槳葉的貼體網(wǎng)格塊，圍繞機翼的貼體網(wǎng)格塊及包圍機翼和槳葉網(wǎng)格塊的背景網(wǎng)格塊。貼體網(wǎng)格塊中由于物體外形并不復雜，均生成純六面體網(wǎng)格。背景網(wǎng)格采用了自適應直角網(wǎng)格，結(jié)合問題實際，進行了網(wǎng)格初始設(shè)定和自適應調(diào)整：一是初始網(wǎng)格設(shè)定，先給定一單元尺度較大的均勻背景網(wǎng)格；二是根據(jù)槳葉、機翼貼體網(wǎng)格塊上的物體表面空間位置(槳葉做一周期的旋轉(zhuǎn)運動)和網(wǎng)格尺度，尋找背景網(wǎng)格上對應網(wǎng)格區(qū)域進行細分和自動網(wǎng)格尺度調(diào)整；三是預定尾跡區(qū)域網(wǎng)格調(diào)整，根據(jù)經(jīng)驗給定旋翼尾跡大致范圍內(nèi)區(qū)域，按照槳葉弦長預估渦核尺度，細分可能的尾跡區(qū)域到合適尺度；四是進行流場初步計算，根據(jù)流場物理量反饋，劃分變量高梯度區(qū)域，生成最終計算網(wǎng)格。經(jīng)過以上步驟得到的計算網(wǎng)格如圖1所示。

定義向上為Y軸方向。旋翼槳盤向后方向為Z軸方向，X軸方向用右手定則確定。圖2給出了本文用于計算討論的槳葉方位角位置(俯視圖)。

圖1 傾轉(zhuǎn)旋翼嵌套網(wǎng)格示意圖Fig.1 Schematic of embedded grid for tilt rotor

圖2 傾轉(zhuǎn)旋翼槳葉方位角示意圖Fig.2 Schematic of blade azimuth angle for tilt rotor

2 數(shù)值方法

2.1 控制方程形式

在地面坐標系下，不計體力等產(chǎn)生的源項，采用格心格式的有限體積法，積分形式的雷諾平均Navier-Stokes方程可以寫為：

(1)

其中，Ω為控制體的體積，S為積分面面積，W為守恒變量，F(xiàn)c和Fv分別為對流和粘性通量項，ρ為密度，p為壓強，V為垂直網(wǎng)格交接面的對流速度，u、v、w為速度，nx、ny、nz為交接面單位外法矢，H為總焓，τ為應力項，T為溫度，K為傳熱系數(shù)。

渦粘性采用一方程SA模型[14]及其對應DES方法[15]進行，積分形式的方程可以寫為：

(2)

DES計算時，不采用最近壁面距離，而是采用混合尺度，特征長度計算采用下式

(3)

d為計算單元到所有壁面的最近距離。因計算采用格心格式，Δ為當前計算單元中心到其所有鄰居單元中心的距離最大值，CDES為可調(diào)節(jié)系數(shù)。由于采用了混合尺度，該方法綜合了大渦模擬方法和RANS方程的優(yōu)點，在近壁面區(qū)域，計算出的尺度與RANS方法相同，在遠離壁面區(qū)域，湍流封閉模式是亞格子模型。

需要注意的是，在旋翼流場計算中，隨槳葉運動，固定的背景網(wǎng)格上計算點到壁面的距離不斷變化，因而需要在每個物理時間，重新計算特征長度。

2.2 求解方法

無粘通量計算時，空間方向采用重構(gòu)方法[16]構(gòu)建二階ROE[17]格式計算交接面上值，粘性通量計算采用中心格式。主控和湍流方程計算采用隱式時間推進[18]，為提高非定常計算的效率，還耦合使用了雙時間方法[19]。物理時間步長為槳葉旋轉(zhuǎn)過方位角1/4°的時間，每個物理時間內(nèi)迭代20次，計算9個周期后，得到相對穩(wěn)定的計算解。

3 算例驗證

3.1 NACA0015翼型較大迎角狀態(tài)氣動力計算

為驗證本文DES方法的計算能力，進行了NACA0015翼型在較大迎角(17°)狀態(tài)下的氣動力計算。計算的翼型弦長是0.3048m，來流馬赫數(shù)是0.29。圖3是計算得到的翼型表面壓強與試驗值[20]的對比?？梢园l(fā)現(xiàn)對于較大迎角計算狀態(tài)，翼型上表面附近發(fā)生大范圍的分離流動時，相對RANS計算，采用DES方法可以更好地預測翼型上表面的壓強分布。

圖3 NACA0015翼型表面壓強系數(shù)比較Fig.3 Comparisons of Cp for NACA0015 airfoil

3.2 旋翼懸停氣動力計算

為驗證本文方法在旋翼氣動力計算方面的能力，選用了被研究者們廣泛采用的有試驗結(jié)果可供對比的“Caradonna & Tung旋翼”[21]為研究模型。該旋翼槳葉的翼型為NACA0012，展長R是1.143 m，展弦比為6，無扭轉(zhuǎn)尖削。為體現(xiàn)出DES方法與RANS方法的計算差異，計算狀態(tài)選擇為槳尖馬赫數(shù)0.794，槳距角為12°。在該狀態(tài)下，槳葉靠近槳尖的部分截面上出現(xiàn)激波，激波和附面層相互干擾，形成槳葉上表面的分離流動。圖4給出了采用DES和RANS方法計算的槳尖附近的槳葉上表面壓強分布。觀察DES結(jié)果，可以發(fā)現(xiàn)隨著槳葉旋轉(zhuǎn)(方位角變化)，槳葉上表面激波后方形成的氣流分離區(qū)域逐漸向槳葉尾緣推移(反映在槳尖附近的槳葉表面壓強隨方位角變化)，到300°方位角時幾乎消失。而采用RANS方法時，各個物理時刻計算出的表面壓強基本相同，表示其未能捕捉到激波后方時變的分離流動現(xiàn)象。

(a) ψ=0°

(b) ψ=60°

(c) ψ=120°

(d) ψ=180°

(e) ψ=240°

(f) ψ=300°

圖5為計算得到的槳尖附近截面(r=0.96R)表面壓強系數(shù)隨槳葉方位角變化結(jié)果?？梢园l(fā)現(xiàn)，在該截面，DES方法計算的槳葉翼型上表面壓強隨方位角變化明顯。采用RANS計算的結(jié)果幾乎隨方位角不變，為清晰起見，只給出其0°方位角時的結(jié)果。

圖5 槳葉截面段翼型表面壓強隨方位角變化(r=0.96R)Fig.5 Surface Cp variation with blade azimuth angle for blade section r=0.96R

圖6為計算得到的各個槳葉截面表面壓強系數(shù)與試驗值的對比(從旋翼旋轉(zhuǎn)第10個周期開始進行3周的氣動力平均結(jié)果)。在大部分截面上，計算和試驗結(jié)果都有較好的一致性。在非?？拷鼧獾慕孛?r=0.96R)，計算和試驗結(jié)果的激波位置存在一定差異。比較RANS和DES計算結(jié)果，發(fā)現(xiàn)其差別體現(xiàn)在靠近槳尖側(cè)的截面(r=0.89R，r=0.96R)，對于遠離槳尖的截面，兩者差異消失。究其原因在于，對應于該計算狀態(tài)，槳葉表面的分離僅發(fā)生在靠近槳尖的激波后方，其對遠離槳葉尖部的槳葉表面附近流動影響有限(r=0.68R，r=0.8R)。在靠近槳尖附近的區(qū)域，DES計算得到的激波位置略偏前緣，而激波后的槳葉上表面壓強高于對應的RANS結(jié)果。

圖6 槳葉截面段翼型表面壓強比較Fig.6 Comparisons of Cp for several blade sections

4 旋翼/機翼干擾狀態(tài)計算結(jié)果分析

4.1 流動現(xiàn)象比較分析

在傾轉(zhuǎn)旋翼的懸停狀態(tài)，流場中典型的流動狀態(tài)即噴泉流動、旋翼下洗及旋翼尾流和機翼干擾后的分離流動。其中標志性的噴泉流動是由于兩個旋翼的尾流沖擊機翼翼面后沿翼面向雙旋翼中線處匯集并卷起后形成。本文計算出的結(jié)果如圖7所示(左為DES結(jié)果，右為RANS結(jié)果)。由圖可以看出，采用DES方法計算出的流場無論是機翼上方中線附近的噴泉流動區(qū)域還是機翼下方的誘導流動，都體現(xiàn)出了左右不對稱的時變特征，而采用RANS方法計算出的流場則對稱性相對較好。

圖8中給出了槳葉方位角0°時，計算得到的旋翼尾流受機翼阻擋，在機翼下方不同翼展平面誘導出的分離流動情況(左為DES結(jié)果，右為RANS結(jié)果)。由圖可見，旋翼誘導的尾流沖擊機翼，受機翼阻擋，形成機翼下方的旋轉(zhuǎn)分離流動。采用RANS計算的結(jié)果一般是預測出靠近機翼前緣和尾緣下方的兩個較大范圍的渦，而采用DES方法則計算出較多的向下發(fā)展的連續(xù)小渦。DES結(jié)果給出了更細致的機翼下方流動分離情況的描述。在靠近機翼尖部位置(X/L=0.8)，兩種計算方法得到的流動形式相接近，這可能是兩個原因?qū)е碌?，一是該位置是旋翼槳根的下方，誘導的下洗速度相對較小，二是機翼尖部的三維橫向流動效應削弱了誘導分離的影響。

(a) DES (b) RANS

圖7噴泉流動示意圖

Fig.7Schematicoffountainflow

(a) X/L=0

(b) X/L=0.2

(c) X/L=0.4

(d) X/L=0.6

(e) X/L=0.8

圖9給出了采用DES方法計算得到的不同槳葉方位角時的特征空間截面(對應機翼位置X=0.4L)上的渦量，可以發(fā)現(xiàn)本文計算較為清晰地捕捉到了旋翼尾跡及機翼下方脫落渦的時間發(fā)展變化歷程。

(a)ψ=0° (b)ψ=30°

(c)ψ=60° (d)ψ=90°

圖9切平面渦量

Fig.9Vorticityincuttingplane

4.2 氣動力比較分析

表1給出了本文計算得到的旋翼(一個)和機翼整體氣動力結(jié)果。由對比得知，采用DES方法計算得到的旋翼拉力和扭矩都大于RANS結(jié)果，但旋翼拉力的差別極小。最大的差別在于機翼向下載荷結(jié)果，采用DES方法計算出的機翼向下載荷較小，比RANS結(jié)果小2.3%左右。導致此項差別的原因可能是采用不同方法時，計算得到的機翼下方誘導流動(特別是機翼前后緣的渦脫落)形式不同。由計算還可以得知，采用DES計算的機翼下洗載荷占整體旋翼拉力的8.34%左右，小于RANS方法對應結(jié)果。

表1 氣動力比較Table 1 Comparisons of aerodynamic

圖10給出了一片槳葉在一個旋翼旋轉(zhuǎn)周期內(nèi)沿截面段的拉力分布變化?？梢园l(fā)現(xiàn)在大部分槳葉方位角，兩種計算方法差異較小，兩者主要差異是槳葉240°到300°方位角之間(正好對應于槳葉通過機翼上方位置)，采用DES方法計算得到了更為劇烈的拉力變化。由圖還可以得知，在槳葉通過機翼上方時，拉力最小，在遠離機翼的方位角時，槳葉拉力得到恢復。

圖10 旋翼旋轉(zhuǎn)周內(nèi)槳葉拉力系數(shù)變化Fig.10 Blade CT in 1 rotor revolution

圖11給出了幾個特征方位角時的槳葉拉力系數(shù)比較，可以發(fā)現(xiàn)在240°到300°方位角之間，槳葉拉力下降主要來源于槳尖部分的氣動力損失。發(fā)生這種現(xiàn)象的原因可能是該時刻槳葉位于機翼上方，兩側(cè)旋翼槳葉靠近且均受到噴泉流動形成的向下氣流影響，槳葉靠近槳尖截面的有效迎角減小。

圖11 槳葉截面段拉力比較Fig.11 Comparisons of blade sectional CL

圖12給出機翼向下載荷的周期變化，由于每個旋翼都由三片槳葉組成，因而旋翼下洗作用在機翼上，形成了周期性的機翼向下載荷。DES方法計算出的機翼向下載荷相對較小，且出現(xiàn)峰值的相位也略靠前。

圖12 機翼向下載荷比較Fig.12 Comparisons of wing download

圖13給出了機翼上表面的壓強等值圖(左為DES結(jié)果，右為RANS結(jié)果)，可以看出在不同槳葉方位角時的機翼上表面高壓區(qū)域的變化，反映出旋翼對機翼的下洗作用影響。

(a) ψ=240°

(b) ψ=270°

(c) ψ=300°

(d) ψ=360°

5 結(jié) 論

在嵌套網(wǎng)格系統(tǒng)下，基于DES方法，進行了傾轉(zhuǎn)旋翼雙旋翼/機翼干擾模型的流場和氣動力計算，并與RANS結(jié)果進行了對比研究，結(jié)果表明：

1) 由于計算的槳葉總距角較小，DES方法得到的旋翼整體氣動力與RANS方法相差不大，但在機翼向下載荷計算方面，DES方法計算得到的向下載荷較??；

2) 在槳葉通過機翼上方時，拉力出現(xiàn)周期內(nèi)最小值，原因可能是向下的噴泉流動減小了槳尖部分截面的迎角。采用DES方法計算時，槳葉通過機翼上方過程的拉力變化更劇烈；

3) 采用DES方法計算的流場左右對稱性相對RANS方法要差些，非定常效應現(xiàn)象更加顯著。機翼下方靠近前緣和后緣的渦結(jié)構(gòu)形式描述也更為細致，表明了該方法在壁面附近存在較大分離流動時的模擬效果更好。

[1]Johnson W. Airloads and wake geometry calculations for an isolated tiltrotor model in a wing tunnel[C]//Presented at the 27th European rotorcraft forum, Moscow, Russia, September 11-14, 2001

[2]Sitataman J, Baeder J D. Analysis of quad tilt rotor blade aerodynamic loads using coupled CFD/Free wake analysis[R]. AIAA 2002-2813, 2002

[3]Li Chunhua, Zhang Jie, Xu Guohua. Computational analysis on tiltrotor aerodynamic characteristics for transitional flight[J]. ACTA Aerodynamica Sinica, 2009, 27(2): 173-205

[4]Yue Hailong, Xia Pinqi. A wake bending unsteady dynamic inflow model of tiltrotor in conversion flight of tiltrotor aircraft[J]. Sci China Ser E-Tech Sci, 2009, 39(12):1992-2000.(in Chinese)岳海龍, 夏品奇. 傾轉(zhuǎn)旋翼機在轉(zhuǎn)換飛行時的旋翼尾跡彎曲非定常動態(tài)入流模型[J]. 中國科學E輯:技術(shù)科學, 2009, 39(12): 1992-2000

[5]Sheng C H, Narramore J C. Computational simulation and analysis of Bell Boeing quadtiltrotor aero interaction[J]. Journal of the American Helicopter Society. 2009, 54: 042002

[6]Fejtek I, Roberts L.Navier-Stokes computation of wing/Rotor interaction for a tiltrotor in hover[J]. AIAA Journal, 1992, 30(11): 2595-2603

[7]Lee Y, Baeder J D. Vortex tracking in overset method for quad tilt rotor blade vortex interaction[R]. AIAA 2003-3531, 2003.

[8]Kjellgren P, Hassan A, Sivasubramanian J, et al. Download alleviation for the XV-15. Computations and experiments of flows around the wing[R]. AIAA 2002-6007, 2002

[9]Potsdam M A, Strawn R C. CFD simulation of tiltrotor configurations in hover[C]//Presented at the American Helicopter Society 58th Annual Forum, Montreal, Canada, June 11-13, 2002

[10]Kim C, Lee J Y. Numerical analysis of hovering tilt-Rotor UAV for minimum download and ground effect analysis[R]. AIAA 2007-1400, 2007

[11]Li Peng, Zhao Qijun, Zhu Qiuxian. CFD calculations on the unsteady aerodynamic characteristics of a tilt rotor in a conversion mode[J]. Chinese journal of aeronautics, 2015, 28(6): 1593-1605

[12]Zhang Ying, Ye Liang, Yang Shuo. Numerical study on flow fields and aerodynamics of tilt rotor aircraft in conversion mode based on embedded grid and actuator model[J]. Chinese journal of aeronautics, 2015, 28(1): 93-102

[13]Li Peng, Zhao Qijun, Calculations on the interaction flowfield and aerodynamic force of tiltrotor/wing in hover[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(2):361-371.(in Chinese)李鵬, 招啟軍. 懸停狀態(tài)傾轉(zhuǎn)旋翼/機翼干擾流場及氣動力的CFD計算[J]. 航空學報, 2014, 35(2): 361-371

[14]Spalart P R, Allmaras S R. A one-equation turbulencemodel for aerodynamic flows[R]. AIAA-92-439, 1992

[15]Spalart P R. Detached-eddy simulation[J]. Annual review of fluid mechanics. 2009, 41:181-202

[16]Frink N T. Recent progress toward a three-dimensional unstructuredNavier-Stokes flow solver[R]. AIAA-94-0061, 1994

[17]Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2):357-372

[18]Luo H, Baum J D. A fast, matrix-free implicit method for computing low mach number flows on unstructured grids[R]. AIAA-99-3315, 1999

[19]Jameson A. Time-dependent calculations using multigrid with applications to unsteady flows past airfoils and wings[R]. AIAA-91-1596, 1991

[20]Piziali R A. 2-D and 3-D oscillating wing aerodynamics for a range of angles of attack including stall[R]. NASA TM-4632, Washington DC,1994

[21]Caradonna F X, Tung C. Experimental and analytical studies of a model helicopter rotor in hover[R]. NASA TM-81232, Moffett Field, CA, 1981.