吳 琪，招啟軍，趙國(guó)慶，王 博

(南京航空航天大學(xué)直升機(jī)旋翼動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)試驗(yàn)室，江蘇南京 210016)

基于隱式算法的懸停旋翼粘性繞流高效CFD分析方法

吳 琪，招啟軍*，趙國(guó)慶，王 博

(南京航空航天大學(xué)直升機(jī)旋翼動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)試驗(yàn)室，江蘇南京 210016)

為顯著提高旋翼粘性繞流CFD模擬效率，建立了一套高效的基于隱式LU-SGS算法和OpenMP并行策略的旋翼懸停流場(chǎng)求解方法。首先，求解Poisson方程生成槳葉剖面翼型的貼體正交網(wǎng)格，并通過剖面網(wǎng)格插值、翻折方法生成槳葉C-O型貼體網(wǎng)格;在此基礎(chǔ)上，采用基于“擾動(dòng)衍射”挖洞方法與“Inverse map”相結(jié)合的洞邊界劃定與貢獻(xiàn)單元搜索方法，解決了嵌套網(wǎng)格技術(shù)中的相關(guān)瓶頸問題。然后，以非慣性系下耦合S-A湍流模型的RANS方程為流場(chǎng)主控方程，對(duì)流通量采用三階Roe-MUSCL格式進(jìn)行離散，時(shí)間推進(jìn)采用高效的隱式LU-SGS方法，同時(shí)采用基于數(shù)據(jù)共享的OpenMP并行策略加速流場(chǎng)求解。最后，運(yùn)用所建立的方法分別對(duì)不同旋翼翼型和懸停狀態(tài)“Caradonna-Tung”以及UH-60A旋翼的流場(chǎng)及氣動(dòng)特性進(jìn)行了計(jì)算，并給出了根據(jù)渦核位置加密網(wǎng)格來提高槳尖渦捕捉精度的方法，同時(shí)將計(jì)算結(jié)果與試驗(yàn)值進(jìn)行了對(duì)比，驗(yàn)證了該方法在旋翼CFD流場(chǎng)模擬中的高效性和高精度特征。

旋翼;氣動(dòng)特性;RANS方程;隱式算法;Spalart-Allamaras湍流模型;嵌套網(wǎng)格方法

0 引言

旋翼是直升機(jī)的關(guān)鍵部件，它提供了直升機(jī)飛行所需的升力、推進(jìn)力及操縱力，而懸停狀態(tài)作為直升機(jī)重要且特有的運(yùn)動(dòng)狀態(tài)，其旋翼流場(chǎng)具有根部為低速流動(dòng)、尖部為跨聲速流動(dòng)，并會(huì)出現(xiàn)強(qiáng)烈集中的螺旋槳尖渦，對(duì)旋翼氣動(dòng)特性有重要的影響，因而高效及高精度的模擬懸停狀態(tài)下旋翼流場(chǎng)的渦尾跡特征及氣動(dòng)特性一直是直升機(jī)技術(shù)領(lǐng)域的研究熱點(diǎn)與難點(diǎn)之一［1］，同時(shí)它也為直升機(jī)旋翼前飛流場(chǎng)的模擬提供了理論基礎(chǔ)。

相比于傳統(tǒng)的試驗(yàn)方法和理論方法，基于RANS方程的先進(jìn)CFD數(shù)值模擬方法具有周期較短、成本較低，并且能充分捕捉流場(chǎng)細(xì)節(jié)尤其是跨聲速流場(chǎng)的激波特性、渦的形成及輸運(yùn)特性等優(yōu)點(diǎn)。隨著計(jì)算機(jī)技術(shù)水平不斷發(fā)展，旋翼CFD方法取得了長(zhǎng)足的進(jìn)步。20世紀(jì)80年代，Agarwal等［2］通過求解無粘可壓Euler方程得到了懸停狀態(tài)下的旋翼流場(chǎng)。90年代以后，Srinivasan等［3］首次采用NS方程計(jì)算了模型旋翼懸停狀態(tài)下的流場(chǎng)，計(jì)算得到的槳葉表面壓強(qiáng)分布與試驗(yàn)值吻合的較好。但由于數(shù)值耗散過大，使得槳尖渦的模擬存在著失真現(xiàn)象。Pomin和 Wagner等［4］則采用嵌套網(wǎng)格技術(shù)和NS方程對(duì)ONERA7A旋翼的懸停流場(chǎng)及氣動(dòng)性能進(jìn)行了數(shù)值模擬，取得了良好效果。最近，Duraisamy等［5］發(fā)展了一種基于重疊網(wǎng)格和WENO格式的旋翼懸停粘性流場(chǎng)的數(shù)值模擬，計(jì)算結(jié)果揭示了高階格式對(duì)槳尖渦模擬的有效性。與國(guó)外相比，國(guó)內(nèi)學(xué)者在旋翼懸停流場(chǎng)的研究中也取得了很大的進(jìn)展。江雄等［6］在國(guó)內(nèi)較早采用雙時(shí)間方法和重疊網(wǎng)格技術(shù)求解NS方程獲得了懸停旋翼的流場(chǎng)特性。招啟軍等［7］采用嵌套網(wǎng)格技術(shù)和高效的貢獻(xiàn)單元搜索方法進(jìn)行了懸停旋翼繞流的數(shù)值模擬。楊愛明等［8］則采用多重網(wǎng)格方法對(duì)旋翼懸停流場(chǎng)進(jìn)行了數(shù)值模擬。之前這些研究揭示了運(yùn)用CFD方法模擬旋翼粘性繞流的可能性，但它們大多是以顯式方法為主，計(jì)算效率往往不高;同時(shí)對(duì)旋翼粘性繞流的模擬也往往局限于采用簡(jiǎn)單的代數(shù)湍流模型(B-L模型［9］)，難以捕捉復(fù)雜狀態(tài)下的流動(dòng)細(xì)節(jié)特征，特別是旋翼大拉力狀態(tài)下的失速特性及渦細(xì)節(jié)特征的高精度模擬。

鑒于此，為了高效及高精度的模擬旋翼粘性繞流的特征，本文發(fā)展了一套耦合嵌套網(wǎng)格技術(shù)、S-A湍流模型［10］、LU-SGS隱式算法［11］、OpenMP并行策略的RANS方程求解方法，其中對(duì)流通量采用高階Roe-MUSCL格式［12］進(jìn)行離散，特別的對(duì)于附加的渦運(yùn)動(dòng)粘性系數(shù)方程，本文著重考慮了方程源項(xiàng)對(duì)隱式離散方程的貢獻(xiàn)并進(jìn)行推導(dǎo)。通過對(duì)不同旋翼翼型和懸停狀態(tài)C-T以及UH-60A旋翼粘性繞流的數(shù)值模擬，結(jié)果表明了該方法相較于傳統(tǒng)顯式推進(jìn)方法能顯著地提高計(jì)算效率，同時(shí)相比于代數(shù)湍流模型，耦合S-A模型的RANS方程求解方法能更好的模擬旋翼翼型、旋翼三維流場(chǎng)細(xì)節(jié)特征及氣動(dòng)特性。

1 旋翼嵌套網(wǎng)格系統(tǒng)

1.1 旋翼C-O型網(wǎng)格生成方法

旋翼翼型網(wǎng)格是旋翼網(wǎng)格生成的基礎(chǔ)，考慮粘性對(duì)流場(chǎng)的影響，本文采用求解Poisson方程的網(wǎng)格生成方法生成旋翼翼型剖面網(wǎng)格，并通過Higenstock［13］源項(xiàng)修正策略調(diào)整網(wǎng)格的疏密及正交程度。為了能準(zhǔn)確的模擬附面層粘性流動(dòng)，第一層翼型網(wǎng)格的壁面距離為5×10-6c(y+≈1)，其中c為翼型弦長(zhǎng)。

在剖面翼型網(wǎng)格的基礎(chǔ)上，對(duì)于槳葉段網(wǎng)格通過槳葉展向剖面間插值生成，而對(duì)于槳根、槳尖處考慮到C-H型網(wǎng)格易出現(xiàn)大變形、非光滑過渡等問題，則采用繞翼型中弧線翻折策略生成槳葉的C-O型結(jié)構(gòu)網(wǎng)格。為了避免槳葉網(wǎng)格在嵌套時(shí)超出背景區(qū)域邊界，靠近槳根位置應(yīng)采取相關(guān)塌縮措施。圖1則給出了針對(duì)UH-60A后掠槳葉網(wǎng)格的插值和翻折過程示意圖。

圖1 UH-60A槳葉網(wǎng)格插值、翻折過程示意圖Fig.1 Schematic of the interpolation and folding process of UH-60A blade grid generation

1.2 嵌套網(wǎng)格生成方法

本文采用嵌套網(wǎng)格方法對(duì)旋翼懸停流場(chǎng)進(jìn)行數(shù)值模擬，其中嵌套網(wǎng)格系統(tǒng)由圍繞旋翼槳葉的C-O型貼體結(jié)構(gòu)網(wǎng)格和圓柱背景網(wǎng)格兩部分組成?？紤]懸停流場(chǎng)具有軸對(duì)稱性質(zhì)，采用周期性邊界條件，只需取1/N(N為槳葉片數(shù))全場(chǎng)。為了解決嵌套網(wǎng)格技術(shù)中的相關(guān)難點(diǎn)，采用高效魯棒的“擾動(dòng)衍射”方法(DDM)快速確定槳葉在背景網(wǎng)格中的洞邊界，同時(shí)結(jié)合Inverse Map方法對(duì)背景網(wǎng)格洞單元的貢獻(xiàn)單元進(jìn)行快速搜尋。當(dāng)采用DDM方法時(shí)，選取圍繞槳葉的包絡(luò)面因盡量避免靠近槳葉表面的非線性流動(dòng)區(qū)域，從而保證三線性插值的精度。圖2和3分別給出了針對(duì)UH-60A懸停旋翼的嵌套網(wǎng)格和背景網(wǎng)格洞邊界單元圖。

圖2 UH-60A槳葉嵌套網(wǎng)格圖Fig.2 Embedded grid system of UH-60A blade

圖3 UH-60A槳葉背景網(wǎng)格洞邊界單元圖Fig.3 Hole boundary cells in the background grid of UH-60A blade

2 控制方程

對(duì)于懸停旋翼，定義坐標(biāo)軸固連于旋轉(zhuǎn)槳葉上，采用準(zhǔn)定常流動(dòng)方法計(jì)算槳葉的流場(chǎng)，將守恒積分形式的可壓RANS方程與Spalart-Allmaras湍流模型進(jìn)行耦合，可得到如下的控制方程:

對(duì)于無粘通量、粘性通量、旋轉(zhuǎn)通量，其表達(dá)式分別為:

對(duì)于旋轉(zhuǎn)流動(dòng)，為了更好的限制S-A湍流模型可能引起的粘性過大的問題，本文在標(biāo)準(zhǔn)模型的基礎(chǔ)上采用了文獻(xiàn)［14］提出的對(duì)生成源項(xiàng)珘S的修正公式:

其中，Ωij代表旋轉(zhuǎn)張量，Sij代表變形張量，式中模型其它相關(guān)參數(shù)可參考文獻(xiàn)［10］的取值。

3 數(shù)值計(jì)算方法

3.1 時(shí)間推進(jìn)方法

為了提高流場(chǎng)求解的效率，本文的時(shí)間推進(jìn)方法采用隱式LU-SGS格式。此外，采用隱式格式也可以克服模型方程的剛性，提高計(jì)算的穩(wěn)定性。對(duì)時(shí)間導(dǎo)數(shù)進(jìn)行離散后得:

其中，ψn=Wn+1－Wn，代入上式，利用

對(duì)于懸停流場(chǎng)，由于計(jì)算的RANS方程耦合了渦運(yùn)動(dòng)粘性系數(shù)的輸運(yùn)方程，特別需要考慮該方程的源項(xiàng)。忽略源項(xiàng)Tsa中高階項(xiàng)對(duì)偏導(dǎo)數(shù)的影響，則可得

式(8)采用LU-SGS格式推進(jìn)求解，則

其中LU分解矩陣L、D、U的表達(dá)式分別為

3.2 空間離散

對(duì)于控制方程，空間離散采用有限體積法，對(duì)于交界面上的對(duì)流通量，采用ROE格式計(jì)算無粘通量:

網(wǎng)格面上的流動(dòng)變量采用三階MUSCL格式插值獲得:上式中，Δ－、Δ+分別為后差、前差算子。參數(shù)κ取為κ=1/3。為了避免在非線性流動(dòng)區(qū)域插值引起的數(shù)值振蕩，同時(shí)引入Venkatakrishnan限制器:

對(duì)于粘性通量采用二階中心格式進(jìn)行計(jì)算。

3.3 初始條件與邊界條件

1)初始條件

流場(chǎng)守恒變量的初始值由來流給定，其初值為:

其中，Re為雷諾數(shù)。

2)邊界條件

3.4 并行方法

為了進(jìn)一步提高流場(chǎng)計(jì)算效率，本文引入基于數(shù)據(jù)共享存儲(chǔ)體系結(jié)構(gòu)的OpenMP并行計(jì)算方法［15］。采用該方法時(shí)，只需在原有串行程序基礎(chǔ)上對(duì)每個(gè)計(jì)算循環(huán)采取互不相關(guān)單元分區(qū)同時(shí)求解，從而加速流場(chǎng)求解。

對(duì)于下文算例3，其槳葉網(wǎng)格大小為225×49× 57，背景網(wǎng)格大小為91×138×117。計(jì)算基于四核CPU配置(單核主頻3.4 GHz)的計(jì)算機(jī)，采用串行計(jì)算需耗時(shí)3.9 h，而采用并行計(jì)算只需耗時(shí)1.37 h，加速效率達(dá)到了65%，可以看出明顯加速了流場(chǎng)的求解。

4 算例計(jì)算與分析

4.1 二維旋翼翼型數(shù)值模擬

為了驗(yàn)證本文發(fā)展的耦合S-A模型的RANS方程隱式求解算法在旋翼翼型繞流的求解精度及計(jì)算效率，選擇了以下2個(gè)算例。其計(jì)算條件如表1所示，其中:Ma為馬赫數(shù)，α為迎角。

表1 二維旋翼翼型算例參數(shù)表Table 1 Parameters of calculation statesfor different rotor airfoils

4.1.1 RAE 2822翼型跨聲速粘性繞流模擬

對(duì)于跨聲速粘性繞流，算例采用C型結(jié)構(gòu)網(wǎng)格來計(jì)算翼型流場(chǎng)，網(wǎng)格大小為369×65，其中在翼型表面布置了302個(gè)網(wǎng)格單元;為了能較好的捕捉附面層內(nèi)流場(chǎng)細(xì)節(jié)，第一層網(wǎng)格的壁面距離取為5×10-6c (y+≈1)，遠(yuǎn)場(chǎng)邊界取大約20倍的弦長(zhǎng)。

圖4給出了該狀態(tài)下基于耦合S-A湍流模型隱式LU-SGS方法計(jì)算和顯式Runge-Kutta方法計(jì)算的密度殘值收斂時(shí)間曲線對(duì)比。當(dāng)采用隱式算法時(shí)，密度殘值收斂4.5個(gè)量級(jí)時(shí)需要的迭代步數(shù)約為顯式的1/4左右，而總的計(jì)算時(shí)間則約為顯式的3/10左右。表明該方法既增加了方程的收斂穩(wěn)定性，又提高了計(jì)算效率。

圖4 C1狀態(tài)顯/隱式密度殘值隨時(shí)間收斂曲線對(duì)比Fig.4 Comparisons of time convergence history of density residual for explicit and implicit scheme at state C1

圖5給出了RAE2822翼型基于耦合S-A湍流模型隱式求解的翼型表面壓強(qiáng)系數(shù)分布與試驗(yàn)值［16］的對(duì)比。作為比較，同時(shí)也計(jì)算了基于B-L湍流模型的壓強(qiáng)系數(shù)分布。從圖5可知，基于耦合S-A湍流模型的壓強(qiáng)系數(shù)分布相比于B-L模型的計(jì)算結(jié)果更加貼近于試驗(yàn)值，且能更好地捕捉翼型上表面的繞流與激波位置。

表2則分別給出了該狀態(tài)下基于耦合S-A和BL湍流模型計(jì)算的氣動(dòng)力系數(shù)與試驗(yàn)值的對(duì)比。從表2可知，與B-L模型相比，基于S-A湍流模型計(jì)算的翼型阻力系數(shù)CD和力矩系數(shù)Cm與試驗(yàn)值更加貼近些，對(duì)于翼型升力系數(shù)CL兩者相差不大。

圖5 C1狀態(tài)翼型表面壓強(qiáng)系數(shù)對(duì)比Fig.5 Comparisons of calculated pressure coefficient distributions with experimental data at state C1

表2 C1狀態(tài)翼型氣動(dòng)力系數(shù)計(jì)算值與試驗(yàn)值的對(duì)比Table 2 Comparisons of calculated aerodynamic force coefficient with experimental data of airfoil at state C1

4.1.2 NACA 0012翼型氣動(dòng)力曲線計(jì)算

在模擬翼型跨聲速粘性繞流的基礎(chǔ)上，算例2則對(duì)NACA 0012翼型在不同迎角下的升力及阻力系數(shù)曲線進(jìn)行了數(shù)值模擬。圖6和7分別給出了基于耦合S-A或B-L湍流模型隱式求解的翼型升力和阻力系數(shù)隨迎角的變化曲線與試驗(yàn)值［17］的對(duì)比。由圖6可知，基于耦合S-A湍流模型的計(jì)算要比采用B-L湍流模型能更好的模擬翼型升力系數(shù)的靜態(tài)失速特性;對(duì)于翼型阻力系數(shù)的模擬，圖7表明采用耦合S-A模型的計(jì)算結(jié)果更貼近試驗(yàn)值。圖8則對(duì)比了采用耦合S-A模型與B-L模型計(jì)算的翼型在18°迎角下的流線圖，可以看出采用S-A湍流模型能很好的捕捉旋翼翼型失速的分離特征，而此時(shí)采用B-L模型卻捕捉不到，這與圖6的結(jié)果是相對(duì)應(yīng)的。

圖6 C2狀態(tài)下翼型升力系數(shù)曲線對(duì)比Fig.6 Comparisons of airfoil lift coefficient curve at state C2

圖7 C2狀態(tài)下翼型阻力系數(shù)曲線對(duì)比Fig.7 Comparisons of airfoil drag coefficient curve at state C2

圖8 C2狀態(tài)18°迎角下不同模型計(jì)算的翼型流線圖Fig.8 Comparisons of calculated streamline of airfoil at 18 degree AOA and state C2 by different models

4.2 懸停狀態(tài)旋翼氣動(dòng)特性數(shù)值模擬

在旋翼翼型數(shù)值分析的基礎(chǔ)上，進(jìn)一步驗(yàn)證本文發(fā)展的耦合S-A湍流模型的RANS方程隱式求解算法在三維旋翼懸停流場(chǎng)的計(jì)算精度及效率。本文采用嵌套網(wǎng)格方法分別對(duì)“Caradonna-Tung”(C-T)模型旋翼和UH-60A旋翼的氣動(dòng)特性進(jìn)行了數(shù)值模擬，計(jì)算狀態(tài)如表3所示，其中:Matip代表槳尖馬赫數(shù)，θ代表總距角，Re代表雷諾數(shù)。

表3 不同旋翼懸停計(jì)算狀態(tài)參數(shù)表Table 3 Parameters of calculation states for different rotors in hovering flight

4.2.1 C-T旋翼

懸停C-T模型旋翼含有兩片槳葉，槳葉外形為矩形，槳葉半徑為1.143 m，槳葉的展弦比為6。對(duì)于該旋翼，槳葉剖面翼型為NACA0012翼型，無負(fù)扭轉(zhuǎn)。算例的槳葉網(wǎng)格大小為225×49×57，第一層槳葉網(wǎng)格的壁面距離為5×10-6c(y+≈1)。對(duì)于背景網(wǎng)格，其網(wǎng)格大小為81×125×89。

圖9給出了C3狀態(tài)下采用耦合S-A湍流模型和隱式LU-SGS算法計(jì)算與顯式Runge-Kutta算法計(jì)算的密度殘值時(shí)間收斂曲線對(duì)比。當(dāng)采用隱式算法密度殘值收斂到目標(biāo)標(biāo)準(zhǔn)只需迭代2 579步，而顯式則需迭代10 158步，且隱式的計(jì)算時(shí)間約為顯式的1/3左右，從而有效地提高了懸停流場(chǎng)的計(jì)算效率。圖10則給出了隱式算法求解的旋翼拉力系數(shù)(CT)及懸停效率(Figure of Merit，F(xiàn)M)的殘值收斂曲線，可見當(dāng)密度殘值收斂到目標(biāo)量級(jí)時(shí)，旋翼拉力系數(shù)與懸停效率也達(dá)到較好的收斂。

圖9 C3狀態(tài)旋翼槳葉密度殘值隨時(shí)間收斂曲線對(duì)比Fig.9 Comparisons of time convergence history of density residual of blade for explicit and implicit scheme at state C3

圖10 C3狀態(tài)旋翼拉力系數(shù)與懸停效率收斂曲線Fig.10 Comparisons of convergence history of rotor thrust coefficient and FM at state C3

圖11則給出了基于耦合S-A湍流模型的RANS方程隱式計(jì)算的槳葉不同剖面的壓強(qiáng)系數(shù)分布的計(jì)算結(jié)果與試驗(yàn)值［18］的對(duì)比，同時(shí)圖中也給出了采用B-L湍流模型的計(jì)算結(jié)果。相比于B-L湍流模型，耦合S-A湍流模型模擬的計(jì)算結(jié)果更接近于試驗(yàn)值，特別的對(duì)于槳尖位置處負(fù)壓峰值的模擬具有相對(duì)更高的計(jì)算精度。

為了能更好的捕捉旋翼的槳尖渦流動(dòng)，以粗網(wǎng)格計(jì)算得到的渦心軌跡位置為參考，在粗網(wǎng)格的基礎(chǔ)上在槳葉展向0.7R～1.05R以及槳葉下方位置0～0.6R內(nèi)的區(qū)域內(nèi)進(jìn)行了適當(dāng)?shù)木W(wǎng)格加密，加密后的細(xì)網(wǎng)格數(shù)目為91×138×117。圖12和13則給出了采用不同大小背景網(wǎng)格計(jì)算獲得的渦量等值面圖，以及渦齡角90度截面上的無量綱等渦量線圖，可以看到對(duì)于粗網(wǎng)格在90°截面上僅能捕捉到三個(gè)集中槳尖渦，而細(xì)網(wǎng)格則能捕捉到四個(gè)集中槳尖渦。圖中由上往下的槳尖渦的渦齡角相差180°，且隨著渦齡角的增大，漩渦的強(qiáng)度逐漸減弱。這表明根據(jù)渦心軌跡位置適當(dāng)加密背景網(wǎng)格能有效提高對(duì)懸停狀態(tài)渦尾跡的捕捉精度。

圖11 C3狀態(tài)旋翼槳葉不同剖面壓強(qiáng)系數(shù)分布對(duì)比Fig.11 Pressure coefficient distributions on the rotor blade at state C3

圖12 C3狀態(tài)采用粗網(wǎng)格與細(xì)網(wǎng)格計(jì)算的渦量等值面對(duì)比Fig.12 Comparisons of calculated vorticity iso-surface based on the coarse and fine grid at state C3

圖13 C3狀態(tài)采用粗網(wǎng)格和細(xì)網(wǎng)格計(jì)算的90°渦齡角截面上的無量綱等渦量線Fig.13 Comparisons of dimensionless vorticity line based on the coarse and fine grid at 90 wake age at state C3

圖14 C4狀態(tài)旋翼槳葉不同剖面壓強(qiáng)系數(shù)分布對(duì)比Fig.14 Pressure coefficient distributions on the rotor blade at state C4

圖14分別給出了C4跨聲速狀態(tài)下的槳葉不同剖面壓強(qiáng)系數(shù)分布的計(jì)算結(jié)果與試驗(yàn)值［18］的對(duì)比，并同時(shí)對(duì)比了S-A和B-L湍流模型的計(jì)算結(jié)果。相比于B-L湍流模型，耦合S-A湍流模型求解能更好的模擬旋翼在有激波狀態(tài)下的跨聲速流場(chǎng)特征，特別是激波位置捕捉更精確。圖15則給出了采用耦合S-A模型計(jì)算的槳葉表面流線圖以及壓強(qiáng)分布圖，由圖可見在當(dāng)前狀態(tài)下近槳尖位置槳葉上表面有激波出現(xiàn)，而對(duì)應(yīng)的區(qū)域也發(fā)生了激波誘導(dǎo)的分離以及分離后的再附現(xiàn)象。圖16則給出了采用本文方法計(jì)算細(xì)網(wǎng)格得到的旋翼槳尖渦渦核的軸向位置與徑向位置。對(duì)比試驗(yàn)數(shù)據(jù)，表明該方法能精確的捕捉槳尖渦的空間位置。圖17給出了模擬得到的C4狀態(tài)90°渦齡角下的旋翼渦量圖，圖中清晰的反映了懸停狀態(tài)旋翼尾跡的收縮特性及槳尖渦強(qiáng)度隨高度不斷衰減的特征，符合實(shí)際物理規(guī)律。

圖15 C4狀態(tài)槳葉表面流線與壓強(qiáng)分布Fig.15 Streamline and pressure contour of blade at state C4

圖16 C4狀態(tài)槳尖渦渦核位置的計(jì)算值與試驗(yàn)值對(duì)比Fig.16 Comparisons of the calculated position of blade-tip vortex core with experimental data at state C4

圖17 C4狀態(tài)90度渦齡角下的旋翼尾跡渦量圖Fig.17 Vorticity contour of rotor wake at 90 wake age and state C4

從C3和C4算例表明本文建立的基于耦合S-A湍流模型的RANS方程隱式求解方法能高效且較好的模擬常規(guī)外形旋翼懸停流場(chǎng)的氣動(dòng)特性以及流場(chǎng)細(xì)節(jié)特征。

4.2.2 UH-60A旋翼

UH-60A旋翼采用了先進(jìn)槳葉氣動(dòng)外形布局，它由四片槳葉組成，每片槳葉半徑為8.18 m，該槳葉的展弦比為15.51。槳葉翼型的配置方式分為三段，它的兩端采用SC1095翼型，中段采用SC1094R8翼型;且槳葉具有非線性負(fù)扭轉(zhuǎn)分布以及20°后掠槳尖，其具體參數(shù)見文獻(xiàn)［19］。

圖18分別給出了UH-60A旋翼在C5狀態(tài)下槳葉不同剖面壓強(qiáng)系數(shù)分布的計(jì)算結(jié)果與試驗(yàn)值［19］的對(duì)比，其中圖中同時(shí)給出了S-A和B-L湍流模型的計(jì)算結(jié)果。相比于B-L湍流模型，采用S-A湍流模型的流場(chǎng)計(jì)算精度更高，特別是能更好的模擬槳尖后掠位置處的壓強(qiáng)分布。圖19還給出了該狀態(tài)下槳葉展向的拉力系數(shù)分布，可知耦合S-A模型求解與試驗(yàn)值吻合的更好，特別是對(duì)于槳尖位置處的剖面拉力系數(shù)的模擬。

圖20給出了該狀態(tài)不同總距角下采用耦合S-A模型或B-L模型計(jì)算的旋翼懸停效率隨拉力系數(shù)的變化趨勢(shì)，耦合S-A模型的計(jì)算結(jié)果在整個(gè)試驗(yàn)范圍內(nèi)與試驗(yàn)值能更好的吻合。這表明本文方法不僅可以準(zhǔn)確模擬常規(guī)氣動(dòng)外形槳葉，也可以有效地對(duì)具有先進(jìn)氣動(dòng)外形的槳葉進(jìn)行流場(chǎng)及氣動(dòng)特性的數(shù)值模擬。

圖18 C5狀態(tài)旋翼槳葉不同剖面壓強(qiáng)系數(shù)分布Fig.18 Pressure coefficient distributions on the rotor blade at state C5

圖19 C5狀態(tài)旋翼槳葉不同剖面拉力系數(shù)分布曲線Fig.19 Prediction of the blade spanwise lift coefficient distribution at state C5

圖20 UH-60A旋翼懸停效率與拉力系數(shù)曲線Fig.20 Thrust coefficient and FM curve of UH-60A rotor

5 結(jié)論

為了提高旋翼粘性繞流的模擬精度和效率，本文發(fā)展了一套耦合嵌套網(wǎng)格技術(shù)、S-A湍流模型、LUSGS隱式算法、OpenMP并行策略的RANS方程求解方法，用于求解不同旋翼翼型和懸停旋翼可壓粘性繞流，通過算例計(jì)算和分析，得到如下結(jié)論:

(1)相比于基于B-L湍流模型的RANS方程求解，耦合S-A湍流模型的求解方法能更好的模擬旋翼翼型及三維旋翼懸停的流場(chǎng)與氣動(dòng)特性，特別能更好的模擬先進(jìn)氣動(dòng)外形旋翼槳尖位置處的流動(dòng)細(xì)節(jié)以及扭矩(阻力)。

(2)相比于傳統(tǒng)顯式Runge-Kutta算法推進(jìn)，采用LU-SGS隱式時(shí)間推進(jìn)和基于數(shù)據(jù)共享的OpenMP并行策略都可以有效地提高計(jì)算效率。

(3)根據(jù)渦核位置適當(dāng)加密背景網(wǎng)格區(qū)域，能更好的捕捉懸停狀態(tài)旋翼槳尖渦的細(xì)節(jié)特征。

［1］ Leishman J G.Rotorcraft aeromechanics:getting through the dip［J］.Journal of the American Helicopter Society，2010，55(1):1-24.

［2］ Agarwal R K，Deese J E.Euler calculations for flowfield of a helicopter rotor in hover［J］.Journal of Aircraft，1987，24(4):231-238.

［3］ Srinivasan G R，Baeder J D，Obayashi S，et al.Flowfield of a lifting rotor in hover-A Navier-Stokes simulation［J］.AIAA Journal，1992，30(10):2371-2378.

［4］ Pomin H，Wagner S.Navier-Stokes analysis of helicopter rotor aerodynamics in hover and forward flight［J］.Journal of aircraft，2002，39(5):813-821.

［5］ Duraisamy K，Baeder J D.High resolution wake capturing methodology for hovering rotors［J］.Journal of the American Helicopter Society，2007，52(2):110-122.

［6］ Jiang Xiong，Chen Zuobin，Zhang Yulun.Numerical simulation of a hovering rotor flowfield using a dual time method［J］.ACTA Aerodynamica Sinica，1998，16(3):288-296.

江雄，陳作斌，張玉倫.用雙時(shí)間法數(shù)值模擬懸停旋翼流場(chǎng)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1998，16(3):288-296.

［7］ Zhao Qijun，Xu Guohua，Zhao Jingen.Numerical simulations of the unsteady flowfield of helicopter rotors on moving embedded grids［J］.Aerospace science and technology，2005，9(2):117-124.

［8］ Yang Aiming，Yang Xiaoquan.Multigrid acceleration and chimera technique for viscous flow past a hovering rotor［J］.Journal of aircraft，2011，48(2):713-715.

［9］ Lomax H，Baldwin B S.Thin layer approximation and algebraic model for separated turbulent flows［R］.AIAA-78-257，1978.

［10］Spalartp R，Allmaras S R.A one-equation turbulence model for aerodynamic flows［R］.AIAA-92-0439，1992.

［11］Yoon S，Jameson A.An Multigrid LU-SSOR Scheme for Approximate Newton Iteration Applied to the Euler Equations ［M］.National Aeronautics and Space Administration，1986.

［12］ Roe P L.Approximate Riemann solvers，parameter vectors，and difference schemes［J］.Journal of Computational Physics，1981，43(2):357-372.

［13］Hilgenstock.A fast method for the elliptic generation of three-dimensional grids with full boundary control［C］//Proceedings of the Second International Conference，1988:137-146.

［14］Dacles-Mariani J，Zilliac G G，Chow J S，et al.Numerical/experimental study of a wingtip vortex in the near field［J］.AIAA Journal，1995，33(9):1561-1568.

［15］Chandra R.Parallel Programming in OpenMP［M］.Morgan Kaufmann，2001.

［16］Holst T L.Viscous transonic airfoil workshop compendium of results［J］.Journal of Aircraft，1988，25(12):1073-1087.

［17］Abbott I H，Von Doenhoff A E.Theory of Wing sections:Including a Summary of Airfoil Data［M］.Courier Corporation，2012.

［18］Caradonna F X，Tung C.Experimental and analytical studies of a model helicopter rotor in hover［J］.Vertica，1981，5(2):149-161.

［19］Abhishek A，Datta A，Chopra I.Prediction of UH-60A structural loads using multibody analysis and swashplate dynamics［J］.Journal of Aircraft，2009，46(2):474-490.

Highly-efficient CFD calculations on viscous flow of hovering rotor based on the implicit algorithm

Wu Qi，Zhao Qijun*，Zhao Guoqing，Wang Bo
(National Key Laboratory of Science and Technology on Rotorcraft Aeromechanics，Nanjing University of Aeronautics and Astronautics Nanjing 210016，China)

To improve the CFD compute efficiency on the viscous flow of helicopter rotor，a highlyefficient implicit LU-SGS algorithm and OpenMP parallel strategy have been employed for predicting the aerodynamic characteristics of hovering rotor.Firstly，the orthogonal and body-fitted grid around rotor blade is generated by using Poisson equations and folding approach;then the“disturbance diffraction method”(DDM)and“Inverse map”method are utilized in the identification of hole cells and searching the donor cells respectively，thus to solve the bottleneck problem of the embedded grid technique.On these bases，the RANS equations coupled S-A turbulence model in a rotating frame are used as the governing equations，and the three-order scheme of ROE-MUSCL is used for spatial discretization of convective fluxes，and the highly-efficient implicit scheme of LU-SGS is adopted for temporal discretization，and data sharing OpenMP parallel strategy is employed to accelerate the calculation as well.Finally，the flowfield and aerodynamic characteristics of several rotor airfoils and hovering rotors(C-T and UH-60A rotor)have been simulated by the presented method.At the same time，a refined grid strategy according to the position of vortex core has been conducted to improve the capturing accuracy of blade-tip vortex.By comparisons of numerical results with experiment data，high-efficiency and high-accuracy abilities of the present method on the simulation of the rotor flowfield is demonstrated.

rotor;aerodynamic characteristics;RANS equations;implicit algorithm;Spalart-Allmaras turbulence model;embedded grid method

V211.3

A

10.7638/kqdlxxb-2014.0035

0258-1825(2015)04-0454-10

2014-05-09;

2014-06-21

國(guó)家自然科學(xué)基金(11272150)

吳琪(1991-)，男，碩士研究生，主要從事旋翼計(jì)算流體力學(xué)與旋翼氣動(dòng)優(yōu)化設(shè)計(jì)等.E-mail:kimiwu91@gmail.com

招啟軍*(1977-)，男，教授，博士生導(dǎo)師，研究方向:直升機(jī)計(jì)算流體力學(xué)、直升機(jī)空氣動(dòng)力學(xué)及流動(dòng)控制、氣動(dòng)噪聲、總體設(shè)計(jì)等.E-mail:zhaoqijun@nuaa.edu.cn

吳琪，招啟軍，趙國(guó)慶，等.基于隱式算法的懸停旋翼粘性繞流高效CFD分析方法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33(4):454-463.

10.7638/kqdlxxb-2014.0035 Wu Q，Zhao Q J，Zhao G Q，et al.Highly-efficient CFD calculations on viscous flow of hovering rotor based on the implicit algorithm［J］.Acta Aerodynamica Sinica，2015，33(4):454-463.