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

吳 琪，招啟軍，趙國慶，王 博

(南京航空航天大學直升機旋翼動力學國家級重點試驗室，江蘇南京 210016)

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

吳 琪，招啟軍*，趙國慶，王 博

(南京航空航天大學直升機旋翼動力學國家級重點試驗室，江蘇南京 210016)

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

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

0 引言

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

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

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

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

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

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

在剖面翼型網(wǎng)格的基礎(chǔ)上，對于槳葉段網(wǎng)格通過槳葉展向剖面間插值生成，而對于槳根、槳尖處考慮到C-H型網(wǎng)格易出現(xiàn)大變形、非光滑過渡等問題，則采用繞翼型中弧線翻折策略生成槳葉的C-O型結(jié)構(gòu)網(wǎng)格。為了避免槳葉網(wǎng)格在嵌套時超出背景區(qū)域邊界，靠近槳根位置應(yīng)采取相關(guān)塌縮措施。圖1則給出了針對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)格方法對旋翼懸停流場進行數(shù)值模擬，其中嵌套網(wǎng)格系統(tǒng)由圍繞旋翼槳葉的C-O型貼體結(jié)構(gòu)網(wǎng)格和圓柱背景網(wǎng)格兩部分組成?？紤]懸停流場具有軸對稱性質(zhì)，采用周期性邊界條件，只需取1/N(N為槳葉片數(shù))全場。為了解決嵌套網(wǎng)格技術(shù)中的相關(guān)難點，采用高效魯棒的“擾動衍射”方法(DDM)快速確定槳葉在背景網(wǎng)格中的洞邊界，同時結(jié)合Inverse Map方法對背景網(wǎng)格洞單元的貢獻單元進行快速搜尋。當采用DDM方法時，選取圍繞槳葉的包絡(luò)面因盡量避免靠近槳葉表面的非線性流動區(qū)域，從而保證三線性插值的精度。圖2和3分別給出了針對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 控制方程

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

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

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

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

3 數(shù)值計算方法

3.1 時間推進方法

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

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

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

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

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

3.2 空間離散

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

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

對于粘性通量采用二階中心格式進行計算。

3.3 初始條件與邊界條件

1)初始條件

流場守恒變量的初始值由來流給定，其初值為:

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

2)邊界條件

3.4 并行方法

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

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

4 算例計算與分析

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

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

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

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

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

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

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

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

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

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

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

4.1.2 NACA 0012翼型氣動力曲線計算

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

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

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

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

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

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

表3 不同旋翼懸停計算狀態(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。對于該旋翼，槳葉剖面翼型為NACA0012翼型，無負扭轉(zhuǎn)。算例的槳葉網(wǎng)格大小為225×49×57，第一層槳葉網(wǎng)格的壁面距離為5×10-6c(y+≈1)。對于背景網(wǎng)格，其網(wǎng)格大小為81×125×89。

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

圖9 C3狀態(tài)旋翼槳葉密度殘值隨時間收斂曲線對比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方程隱式計算的槳葉不同剖面的壓強系數(shù)分布的計算結(jié)果與試驗值［18］的對比，同時圖中也給出了采用B-L湍流模型的計算結(jié)果。相比于B-L湍流模型，耦合S-A湍流模型模擬的計算結(jié)果更接近于試驗值，特別的對于槳尖位置處負壓峰值的模擬具有相對更高的計算精度。

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

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

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

圖13 C3狀態(tài)采用粗網(wǎng)格和細網(wǎng)格計算的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)旋翼槳葉不同剖面壓強系數(shù)分布對比Fig.14 Pressure coefficient distributions on the rotor blade at state C4

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

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

圖16 C4狀態(tài)槳尖渦渦核位置的計算值與試驗值對比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ī)外形旋翼懸停流場的氣動特性以及流場細節(jié)特征。

4.2.2 UH-60A旋翼

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

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

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

圖18 C5狀態(tài)旋翼槳葉不同剖面壓強系數(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方程求解方法，用于求解不同旋翼翼型和懸停旋翼可壓粘性繞流，通過算例計算和分析，得到如下結(jié)論:

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

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

(3)根據(jù)渦核位置適當加密背景網(wǎng)格區(qū)域，能更好的捕捉懸停狀態(tài)旋翼槳尖渦的細節(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ù)值模擬懸停旋翼流場［J］.空氣動力學學報，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

國家自然科學基金(11272150)

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

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

吳琪，招啟軍，趙國慶，等.基于隱式算法的懸停旋翼粘性繞流高效CFD分析方法［J］.空氣動力學學報，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.