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

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

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

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

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

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

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

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

0 引言

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

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

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

1 旋翼嵌套網格系統

1.1 旋翼C-O型網格生成方法

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

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

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

1.2 嵌套網格生成方法

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

圖2 UH-60A槳葉嵌套網格圖Fig.2 Embedded grid system of UH-60A blade

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

2 控制方程

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

對于無粘通量、粘性通量、旋轉通量，其表達式分別為:

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

其中，Ωij代表旋轉張量，Sij代表變形張量，式中模型其它相關參數可參考文獻［10］的取值。

3 數值計算方法

3.1 時間推進方法

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

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

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

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

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

3.2 空間離散

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

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

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

3.3 初始條件與邊界條件

1)初始條件

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

其中，Re為雷諾數。

2)邊界條件

3.4 并行方法

為了進一步提高流場計算效率，本文引入基于數據共享存儲體系結構的OpenMP并行計算方法［15］。采用該方法時，只需在原有串行程序基礎上對每個計算循環采取互不相關單元分區同時求解，從而加速流場求解。

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

4 算例計算與分析

4.1 二維旋翼翼型數值模擬

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

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

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

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

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

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

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

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

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

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

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

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

圖6 C2狀態下翼型升力系數曲線對比Fig.6 Comparisons of airfoil lift coefficient curve at state C2

圖7 C2狀態下翼型阻力系數曲線對比Fig.7 Comparisons of airfoil drag coefficient curve at state C2

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

4.2 懸停狀態旋翼氣動特性數值模擬

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

表3 不同旋翼懸停計算狀態參數表Table 3 Parameters of calculation states for different rotors in hovering flight

4.2.1 C-T旋翼

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

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

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

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

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

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

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

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

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

圖14 C4狀態旋翼槳葉不同剖面壓強系數分布對比Fig.14 Pressure coefficient distributions on the rotor blade at state C4

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

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

圖16 C4狀態槳尖渦渦核位置的計算值與試驗值對比Fig.16 Comparisons of the calculated position of blade-tip vortex core with experimental data at state C4

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

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

4.2.2 UH-60A旋翼

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

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

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

圖18 C5狀態旋翼槳葉不同剖面壓強系數分布Fig.18 Pressure coefficient distributions on the rotor blade at state C5

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

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

5 結論

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

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

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

(3)根據渦核位置適當加密背景網格區域，能更好的捕捉懸停狀態旋翼槳尖渦的細節特征。

［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.

江雄，陳作斌，張玉倫.用雙時間法數值模擬懸停旋翼流場［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-)，男，碩士研究生，主要從事旋翼計算流體力學與旋翼氣動優化設計等.E-mail:kimiwu91@gmail.com

招啟軍*(1977-)，男，教授，博士生導師，研究方向:直升機計算流體力學、直升機空氣動力學及流動控制、氣動噪聲、總體設計等.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.