基于高階WENO格式的旋翼非定常渦流場數(shù)值模擬

印智昭， 招啟軍， 王博

南京航空航天大學 直升機旋翼動力學國家級重點實驗室， 南京　210016

?

基于高階WENO格式的旋翼非定常渦流場數(shù)值模擬

印智昭， 招啟軍*， 王博

南京航空航天大學 直升機旋翼動力學國家級重點實驗室， 南京210016

為提高直升機旋翼黏性渦流場計算流體力學(CFD)的計算精度及降低數(shù)值耗散，發(fā)展了一套基于五階加權(quán)本質(zhì)無振蕩(WENO)格式的旋翼非定常渦流場數(shù)值計算方法。采用運動嵌套網(wǎng)格方法生成圍繞前飛狀態(tài)旋翼的網(wǎng)格系統(tǒng)，以Navier-Stokes方程作為主控方程，湍流模型采用一方程Spalart-Allmaras模型。為了提高旋翼流場中對渦的形成、演化等發(fā)展過程的模擬精度，采用高間斷分辨率的Roe-WENO格式計算對流通量。對狀態(tài)變量的插值，選取適當?shù)募訖?quán)因子，通過多個重構(gòu)模板的凸組合來構(gòu)造五階WENO格式，在交界面附近獲得具有五階精度的狀態(tài)變量。為了提高計算效率，采用高效的隱式LU-SGS(Lower-Upper Symmetric Gauss-Seidel)推進方法和并行計算策略。最后，運用該方法分別對懸停狀態(tài)下C-T(Caradonna-Tung)旋翼和UH-60A旋翼以及前飛狀態(tài)下C-T旋翼和SA349/2旋翼的渦流場及氣動特性進行了數(shù)值模擬，并將數(shù)值結(jié)果與傳統(tǒng)二階精度格式的計算值進行了對比。結(jié)果表明：在計算時間增長10%～20%的前提下，五階WENO格式能夠捕捉到更精確的渦流場特性，反映了五階WENO格式在計算旋翼非定常渦流場時具有更高的計算精度與更低的數(shù)值耗散特性。

旋翼； 渦流； 加權(quán)本質(zhì)無振蕩格式； Navier-Stokes方程； 槳尖渦； 運動嵌套網(wǎng)格

直升機作為20世紀航空領(lǐng)域極具特色的創(chuàng)造之一，其填補了固定翼飛機在垂直飛行領(lǐng)域內(nèi)的空白。直升機所具有的垂直起降、空中懸停和近地飛行等特性，拓展了其在低空領(lǐng)域的適用范圍。但直升機的飛行環(huán)境復雜，尤其是旋翼始終工作在復雜多變的自身渦流環(huán)境中，對其氣動特性造成重要影響。因此為了能夠準確計算出旋翼的氣動性能，同時找到改進的氣動設(shè)計方法，關(guān)鍵在于是否能夠提高旋翼渦流場的計算精度以及準確模擬出槳尖渦的生成與演化過程[1]。因此，對上述問題的研究一直是直升機技術(shù)領(lǐng)域的一個重點和難點。

對于旋翼流場的研究方法，一般可分為試驗方法、理論方法以及旋翼計算流體力學(CFD)方法。相對于前兩種方法，旋翼CFD方法不僅能夠較精確地預估旋翼的氣動特性，且能夠降低旋翼設(shè)計的研究成本與周期，因此在過去的20年內(nèi)旋翼CFD方法得到了長足的發(fā)展。20世紀90年代，Srinivasan等[2]采用了Navier-Stokes方程計算了旋翼懸停流場，得到了較好的槳葉表面壓強分布，但是槳尖渦的模擬存在嚴重的數(shù)值耗散。Pomin和Wagner[3]采用了運動嵌套網(wǎng)格方法對7A模型旋翼懸停流場和C-T(Caradonna-Tung)模型旋翼前飛流場進行了數(shù)值模擬，取得了較好的氣動特性模擬結(jié)果。楊愛明和喬志德[4]采用了運動嵌套網(wǎng)格技術(shù)數(shù)值模擬了前飛非定常黏性流動，前飛數(shù)值結(jié)果與其他文獻計算結(jié)果吻合較好。Zhao等[5]使用混合方法計算了多種旋翼懸停與前飛狀態(tài)下的流場，獲得了較好的氣動特性與載荷分布，但并未對旋翼流場中渦的輸運特性進一步開展研究。從諸多學者的研究來看，傳統(tǒng)旋翼CFD方法日益完善，但是對于旋翼復雜的非定常渦流場，尤其是運動過程中的槳/渦干擾現(xiàn)象，數(shù)值模擬的精度受到低網(wǎng)格分辨率和低階數(shù)值格式的高耗散性制約。為了減弱傳統(tǒng)數(shù)值方法中出現(xiàn)的旋翼渦非物理失真現(xiàn)象，學者們開始將注意力轉(zhuǎn)移到提高數(shù)值方法的計算精度與數(shù)值耗散上。2005年，Hariharan[6]將五階和七階ENO(Essentially Non-oscillatory)格式用于模擬前飛狀態(tài)旋翼黏性渦流場中，計算顯示高階格式有助于提高槳尖渦的捕捉精度。美國國家航空航天局(NASA)對V-22傾轉(zhuǎn)旋翼飛行器的一份研究報告[7]中指出，與二階精度格式相比，五階空間精度格式可以有效提高渦的模擬精度，降低了二階方法中的耗散與色散誤差。

國內(nèi)方面，招啟軍[8]、徐廣[9]等將三階Roe-MUSCL格式應用于旋翼懸停流場計算中。Yu等[10]則針對可壓縮流場提出了基于加權(quán)本質(zhì)無振蕩(Weighted Essentially Non-oscillatory， WENO)格式的高精度的混合型激波捕捉格式，體現(xiàn)了WENO格式的低耗散性特性。李新亮[11]研究了WENO格式等高階保單調(diào)格式在高超聲速湍流數(shù)值模擬中的應用和最新進展。Xu和Weng[12]將WENO格式應用于C-T旋翼懸停流場數(shù)值模擬中，湍流模型采用了B-L模型，得到了WENO格式具有較低數(shù)值耗散的結(jié)論，但并未對其他具有復雜外形的旋翼以及不同前飛狀態(tài)進行拓展模擬。

針對高階格式能夠提高旋翼渦流場的模擬精度，降低數(shù)值耗散等特性，本文將五階WENO格式重點應用于旋翼非定常渦流場計算中，發(fā)展了一套結(jié)合運動嵌套網(wǎng)格方法，基于雷諾平均Navier-Stokes(RANS)方程求解旋翼非定常渦流場的計算方法。湍流模型采用了航空流動計算效果較好的一方程Spalart-Allmaras湍流模型，采用高效的隱式LU-SGS(Lower-Upper Symmetric Gauss-Seidel)雙時間推進方法。著重開展了具有先進氣動外形的UH-60A旋翼和具有非線性負扭轉(zhuǎn)的SA349/2旋翼非定常渦流場的數(shù)值模擬分析。結(jié)果表明五階WENO格式不僅能夠有效地提高旋翼渦流場的計算精度，降低了槳尖渦輸運過程中的數(shù)值耗散，而且其對于具有復雜外形的旋翼渦流場模擬也十分有效，表明了該方法的通用性和精確性，為今后研究復雜外形旋翼槳尖渦輸運特性以及氣動設(shè)計提供了較好的數(shù)值基礎(chǔ)。

1　旋翼運動嵌套網(wǎng)格方法

旋翼在旋轉(zhuǎn)過程中通常會伴隨著揮舞、扭轉(zhuǎn)等復雜運動。為了準確模擬旋翼運動，本文采用了運動嵌套網(wǎng)格方法來實現(xiàn)這一目標。生成三維旋翼貼體網(wǎng)格之前，需要選定旋翼展向的幾個特征剖面。根據(jù)Poisson方程生成旋翼特征截面處翼型的二維C型網(wǎng)格，通過Hilgenstock[13]源項控制網(wǎng)格的疏密和正交程度。特征剖面槳葉段之間的網(wǎng)格可以通過展向剖面網(wǎng)格間插值生成。在槳葉根部與尖部，采用繞翼型中弧線翻折方法生成圍繞旋翼的O-C-O型三維貼體網(wǎng)格。

嵌套網(wǎng)格系統(tǒng)需要解決兩個技術(shù)難點：洞邊界生成以及貢獻單元的搜尋。為解決上述問題，采用“擾動衍射”方法來快速確定槳葉在背景網(wǎng)格中的洞邊界，同時采用Inverse map方法[14]對邊界的貢獻單元進行快速的搜尋。當洞單元及其貢獻單元確定之后，嵌套網(wǎng)格系統(tǒng)之間可以采用三線性插值方法進行信息傳遞，保證傳遞信息的精度。圖1給出了本文所采用的嵌套網(wǎng)格系統(tǒng)和UH-60A旋翼網(wǎng)格生成示意圖。

圖1　運動嵌套網(wǎng)格以及UH-60A旋翼槳葉網(wǎng)格示意圖Fig.1　Schematic of moving-embedded grids and UH-60A rotor blade grids

2　旋翼非定常流場求解方法

2.1流場求解控制方程

為捕捉旋翼的非定常渦流場特性，包括激波位置、徑向流和槳尖渦等流動細節(jié)，保證精度，采用積分形式的RANS方程作為求解控制方程：

(1)

守恒變量、無黏通量和黏性通量的表達式分別為

定義坐標軸于固定坐標系下，絕對速度為V=[uvw]T，牽連速度為Vω=ω×r，ω為旋轉(zhuǎn)角速度，r為網(wǎng)格位置矢量。式中：ρ、p、E和H分別為流體單元密度、壓強、總內(nèi)能和總焓能；τ(·)和Θ(·)為黏性相關(guān)項。

2.2空間離散

采用有限體積法對控制方程進行離散，對于交界面上的對流通量采用Roe格式[15]計算無黏通量：

(2)

(3)

對于參與對比的傳統(tǒng)數(shù)值格式，采用了二階JST(Jameson)格式和二階Roe-MUSCL格式。其中Roe-MUSCL格式中限制器采用了Van Albada限制器。

2.3插值方法

為提高旋翼流場計算的精度，采用了五階WENO格式對Roe格式交界面處的左右狀態(tài)量進行重構(gòu)，本文中選擇原始變量進行重構(gòu)，如圖2 所示。1987年Harten等提出了本質(zhì)無振蕩(Essentially Non-oscillatory， ENO)格式[17]。在ENO格式的基礎(chǔ)之上，Liu等[18]于1994年提出了WENO格式。WENO格式的基本思想是通過多個模板的線性組合低階通量得到高階通量。通過構(gòu)造合適的加權(quán)系數(shù)以及選取適當?shù)墓饣蜃?，可以使得WENO格式的數(shù)值解在間斷附近保持基本無振蕩特性，而且能夠在光滑區(qū)域獲得一致高精度。

圖2　五階加權(quán)本質(zhì)無振蕩(WENO)格式插值示意圖Fig.2　Schematic of interpolation for the 5th weighted essentially non-oscillatory (WENO) scheme

三階ENO格式需要在所有模板中選擇一個最為光滑的模板；而五階WENO格式則是將所有模板進行了加權(quán)組合。對于三階ENO格式，有3個備選模板，每個模板包含3個格點信息。將它們組合在一起可以形成包含5個格點信息的插值區(qū)域。通過插值，可以獲得五階精度的WENO格式，插值表達式可以寫為

(4)

式中：ωk為相對應模板的權(quán)重，與每個模板的光滑性有關(guān)；fk為模板的狀態(tài)變量；qk為模板的插值多項式，一般采取Newton插值可得到。

將式(4)用于狀態(tài)變量的插值計算上，以i方向為例，對于i+1/2處的左狀態(tài)量，插值模板示意圖如圖3所示，q0、q1和q2分別為3個模板。

圖3　左狀態(tài)量插值模板示意圖Fig.3　Schematic of stencils for interpolation of left quantities

其加權(quán)重構(gòu)公式為

(5)

所使用的3個重構(gòu)模板的表達式為(Q為原始變量)

(6)

非線性加權(quán)系數(shù)按照Jiang和Shu[19]的計算方法，其表達式為

(7)

式中：權(quán)系數(shù)分別為C0=0.1，C1=0.6，C2=0.3；ε為避免分母為0而加入的小量，可取為10-6；ISk為模板的光滑度量系數(shù)，用于度量第k個模板的光滑性。采用文獻[18]計算方式，其表達式為

(8)

具體計算如下：

(9)

將式(6)、式(7)和式(9)代入式(5)中，就可計算得到具有五階精度的左狀態(tài)量。對于i+1/2處的右狀態(tài)量，計算方法與上文相同，根據(jù)對稱性便可計算得到。

2.4時間離散

為了提高流場的求解效率，使用基于雙時間推進方法的隱式LU-SGS[20]進行時間推進。對時間導數(shù)進行離散后得到

V(Tn+1-Tn)=Rn

(10)

式中：W、Fc和Fv分別為守恒變量、無黏通量和黏性通量；T為源項通量；τ為偽時間；t為物理時間；m為物理時間上的層數(shù)；n為偽時間上的層數(shù)

對式(10)采用LU-SGS格式進行求解，即

(L+D)D-1(D+U)ΔW=Rn-

(3Wm-4Wn+Wn-1)/2Δt

(11)

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

(12)

2.5湍流模型

(13)

式中：Fc,T、Fv,T和QT分別為黏性輸運方程的對流項、黏性項和源項。該方程求解方法與Navier-Stokes方程相同。

2.6并行策略

采用便于編程的OpenMP并行計算方法[22]，適用于單核多線程計算機，屬于多線程并行。在計算中，使用了基于四核CPU配置的計算機。對于C-T旋翼懸停流場，單核計算串行程序時，計算時長為9.5 h，而采用四核并行計算時只需耗時4.1 h，加速比達到2.32。采用OpenMP并行算法可以較為明顯地提高旋翼流場的計算效率。

3　算例及結(jié)果分析

為驗證本文所發(fā)展的數(shù)值方法的正確性，對懸停狀態(tài)下C-T旋翼和UH-60A旋翼、無升力前飛狀態(tài)下的C-T旋翼和前飛狀態(tài)下的SA349/2旋翼的非定常渦流場進行了數(shù)值模擬。計算狀態(tài)如表1所示(表中：SA349/2旋翼的變距與揮舞運動規(guī)律由式(14)給出；Matip為槳尖馬赫數(shù)；μ為前進比；θ為總距角；β為揮舞角)。

表1　計算狀態(tài)參數(shù)

3.1C-T旋翼懸停狀態(tài)流場模擬

懸停C-T旋翼模型作為經(jīng)典旋翼流場驗證模型，其包含兩片槳葉，槳葉外形為矩形。槳葉半徑為1.143 m，展弦比為6。該槳葉剖面翼型為NACA 0012翼型，無負扭轉(zhuǎn)。因計及黏性，槳葉物面第一層網(wǎng)格距離為5×10-6c(c為弦長)，使用的網(wǎng)格大小為225×49×57。對于背景網(wǎng)格，其大小為121×190×135。

圖4給出了五階Roe-WENO格式與JST格式和Roe-MUSCL格式計算的不同剖面槳葉壓強系數(shù)Cp對比。從圖中可以看出，3種方法的計算結(jié)果都與試驗值[23]吻合較好。在靠近槳尖位置處，五階Roe-WENO格式計算的激波位置要靠前一點，相比于Roe-MUSCL格式和JST格式更加接近試驗結(jié)果。

圖5給出了所計算得到旋翼槳尖渦渦核的軸向位置y/R和徑向位置r/R與試驗值[23]的對比(ψ為渦齡角)，選取渦量范圍為0.2～0.6。從徑向位置來看，兩種方法的計算結(jié)果均與試驗值吻合較好；但是在大渦齡角之后，JST格式與Roe-MUSCL格式表現(xiàn)出了較大的耗散。從軸向位置來看，在渦齡角達到180°之后，JST格式與Roe-MUSCL格式的計算結(jié)果開始不同程度地偏離試驗值，而Roe-WENO格式的計算結(jié)果與試驗值吻合較好，并且能模擬到更大渦齡角上的槳尖渦，說明了本文使用的五階WENO格式的數(shù)值耗散要低于Roe-MUSCL格式的數(shù)值耗散，Roe-MUSCL格式的數(shù)值耗散要低于JST格式。

圖4　不同剖面壓強系數(shù)分布計算值與試驗數(shù)據(jù)對比Fig.4　Comparison of calculated pressure coefficientdistributions and experimental data

圖5　槳尖渦渦核位置計算值與試驗數(shù)據(jù)對比Fig.5　Comparison of calculated position of blade-tip vortex core with experimental data

圖6　C-T旋翼等值渦量面對比Fig.6　Comparison of vorticity iso-surface of C-T rotor

圖7　WENO與JST格式計算兩個截面旋翼尾跡渦量對比Fig.7　Comparison of rotor wake vorticity contours at two cutting planes calculated by JST and WENO schemes

(B1、B2分別代表兩片槳葉)。從圖中可以看出90°方位截面下五階WENO格式可以捕捉到5個強度較為明顯的渦，而JST格式只能捕捉到3個明顯強度的渦，Roe-MUSCL格式捕捉到了4個明顯強度的渦；在0°方位截面下，Roe-WENO格式可以捕捉到5個渦，而JST格式在第2個渦之后就產(chǎn)生了明顯的耗散，出現(xiàn)了渦的擴散與消失；Roe-MUSCL格式的計算結(jié)果優(yōu)于JST格式，與Roe-WENO相比也存在明顯差距。

3.2UH-60A旋翼懸停狀態(tài)流場模擬

世界著名的戰(zhàn)術(shù)通用直升機UH-60A旋翼具有先進氣動外形布局特征。其由4片槳葉組成，槳葉展弦比為15.3。采用多段翼型優(yōu)化配置，主要包括SC1095和SC1094R8旋翼專用翼型等。槳葉具有非線性負扭轉(zhuǎn)，特別是槳尖位置處的“勺子型”扭轉(zhuǎn)分布，且采用了具有20°后掠角的新型槳尖。為了提高計算精度，在槳葉弦向進行了加密處理，旋翼網(wǎng)格大小為277×40×111，背景網(wǎng)格大小為121×108×88。

圖8給出了UH-60A旋翼不同展向剖面位置處壓強系數(shù)與試驗值[24]的對比。從圖中可以看出r/R=0.775位置處，Roe-WENO格式與JST格式的計算結(jié)果均與試驗值吻合較好；在r/R=0.865,0.920位置處，JST格式計算的上表面前緣壓力峰值略低于試驗值，Roe-WENO格式的計算結(jié)果與試驗結(jié)果吻合較好，表現(xiàn)出Roe-WENO格式的高精度特性。

圖9給出了UH-60A槳葉展向的升力系數(shù)CL與試驗值[24]對比。從圖中可以看出，在槳尖位置處，Roe-WENO格式的計算結(jié)果與試驗值比較吻合，而JST格式的計算結(jié)果偏大。圖10 給出了懸停效率FM隨著拉力系數(shù)CT的變化與試驗值[24]對比。從圖中可以看出，Roe-WENO格式的計算結(jié)果在大拉力的情況下與試驗值吻合更好。

圖8　不同剖面壓強系數(shù)分布計算值與試驗值對比Fig.8　Comparisons of calculated pressure coefficient distributions at different spanwise sections and experimental data

圖9　UH-60A旋翼槳葉展向升力系數(shù)分布計算對比Fig.9　Comparison of calculated spanwise lift coefficient distributions of UH-60A rotor

圖10　UH-60A旋翼懸停效率(FM)隨拉力系數(shù)變化Fig.10　Variation of figure merit (FM) curve as thrustcoefficient of UH-60A rotor

圖11　UH-60A旋翼渦量等值面對比Fig.11　Comparison of vorticity iso-surface of UH-60A rotor

3.3C-T旋翼無升力前飛狀態(tài)流場模擬

對于前飛試驗狀態(tài)下的C-T旋翼，其槳葉展弦比為7。計算使用的旋翼網(wǎng)格尺寸為265×49×96，背景網(wǎng)格尺寸為160×135×152。槳葉網(wǎng)格第一層高度為5×10-6c。計算物理時間步長取為1/120旋轉(zhuǎn)周期，即為3°方位角，CFL數(shù)取為5，計算圈數(shù)為3圈。

圖12給出了無升力前飛狀態(tài)下WENO格式與JST格式計算的r/R=0.893剖面壓強系數(shù)與試驗值[25]的對比。從圖中可以看出，在各個方位角下的計算值與試驗值均吻合較好。從方位角60°到方位角150°，旋翼上下表面均出現(xiàn)了激波；而且從對比中可以看出，WENO格式計算值比JST計算值更加貼近試驗值。計算結(jié)果表明了五階WENO格式可以有效地應用于旋翼前飛非定常渦流場數(shù)值模擬中，而且能夠表現(xiàn)出較好的數(shù)值精度。

圖12　槳葉不同方位角壓強系數(shù)計算值與試驗值對比 (r/R=0.893)Fig.12　Comparisons of calculated pressure coefficient distributions at different azimuthal angles with experimental data (r/R=0.893)

3.4SA349/2旋翼前飛狀態(tài)流場模擬

SA349/2旋翼[25]具有3片槳葉，槳葉形狀為矩形，展弦比為15。槳葉采用了OA209翼型，而且具有非常規(guī)的扭轉(zhuǎn)分布形式(非線性幾何負扭轉(zhuǎn))。計算使用的旋翼網(wǎng)格尺寸為237×45×61，背景網(wǎng)格尺寸為160×135×152。槳葉網(wǎng)格第一層高度為5×10-6c。計算物理時間步長取3°方位角，CFL數(shù)取為5，計算圈數(shù)為3圈。變距和揮舞運動規(guī)律為

(14)

圖13給出了SA349/2旋翼前飛狀態(tài)下剖面r/R=0.970處壓強系數(shù)計算值與試驗值對比[26]。從圖中可以看出，五階WENO格式與JST格式的計算結(jié)果均與試驗結(jié)果吻合較好，前者與試驗結(jié)果更接近一些，尤其是在激波位置和負壓峰值位置。

圖14(a)和圖14(b)分別從俯視角和仰視角下給出了渦量等值面圖對比。從圖中可以看出，渦量值取為0.015時，與JST格式計算結(jié)果相比，五階WENO格式能夠捕捉到更為完整的渦結(jié)構(gòu)和渦細節(jié)。

圖13　槳葉不同方位角壓強系數(shù)計算值與試驗值對比 (r/R=0.970)Fig.13　Comparisons of calculated pressure coefficient distributions at different azimuthal angles with experimental data (r/R=0.970)

為進一步顯示兩種方法在模擬旋翼渦量細節(jié)方面的差異，圖15則給出了前飛方向x上的若干個截面上的等渦量云圖，B1、B2和B3分別為3片槳葉。圖16則給出了前進方向上特征截面(x/R=1.0，1.6，2.6，3.4)上的渦量云圖。從這兩幅圖上可以看出，隨著槳尖渦的向后拖出，WENO格式計算的渦的強度減弱要明顯小于JST格式的計算結(jié)果；而且在x/R=1.0，1.6位置處，從WENO格式的計算結(jié)果可以看出不同槳葉拖出的槳尖渦會出現(xiàn)較為明顯的干擾現(xiàn)象(渦的融合與分離)，表明了五階WENO格式要比JST格式具有更小的數(shù)值耗散。

圖17則進一步給出了z/R=±0.87特征截面位置處的等渦量云圖分布，圖中的B1-1為第1片槳葉的第1圈。從圖17(a)中也可以清楚地看出，五階WENO格式可以計算出第1片槳葉與第2片槳葉的第2圈渦量，而JST格式則出現(xiàn)了較大的耗散。從圖17(b)中也可以看出該截面上WENO格式捕捉到了更多旋翼槳葉渦量的細節(jié)，這為旋翼槳/渦干擾特性的分析提供了良好的基礎(chǔ)。

圖14　前飛旋翼渦量等值面圖Fig.14　Comparison of vorticity iso-surface of rotor in forward flight

圖15　前飛時若干截面的旋翼等渦量云圖分布Fig.15　Vorticity contours of rotor at several cutting planes in forward flight

圖16　前飛情況下x/R=1.0, 1.6, 2.6, 3.4特征截面的旋翼等渦量云圖分布Fig.16　Vorticity contours of rotor at several cutting planes in forward flight (x/R=1.0, 1.6, 2.6, 3.4)

圖17　前飛時不同截面的旋翼等渦量云圖分布 (z/R=-0.87, 0.87)Fig.17　Comparison of vorticity contours of rotor at different cutting planes in forward flight (z/R=-0.87, 0.87)

圖18則給出了SA349/2旋翼槳葉特征剖面0.97R位置處在不同方位角ψ的剖面法向力系數(shù)(CN)與試驗值[26]對比。從圖中可以看出，JST和Roe-WENO格式均能夠較好地模擬出該狀態(tài)下氣動力幅值和相位相應的變化規(guī)律。

圖18　SA349/2旋翼槳葉r/R=0.97剖面法向力系數(shù)分布對比Fig.18　Comparison of normal forces coefficient distribution at the position of r/R=0.97 of SA349/2 rotor

3.5不同算例計算時間對比

本文所有算例均是在四核CPU，單核主頻3.4GHz計算機上進行并行計算。為比較WENO格式與JST格式的計算量效率，對相同計算條件下的計算時間進行了對比，具體結(jié)果如表2所示。從表中可以看出，五階WENO格式的計算量稍大于JST格式的計算量。在懸停與前飛狀態(tài)下，WENO格式的計算時間比JST格式的多10%～20%。從表中還可以看出，前飛狀態(tài)C-T與SA349/2網(wǎng)格總量相當，但計算時間相差較大。這是因為計算過程中旋翼流場部分占據(jù)大部分的計算資源與計算時間，而本文所使用的C-T旋翼網(wǎng)格量大小約為SA349/2旋翼網(wǎng)格量大小的2倍，因此時間也大約為2倍左右。

表2　不同計算算例時間對比

4　結(jié)　論

為了提高旋翼非定常渦流場的計算精度及降低數(shù)值耗散，本文發(fā)展了一套基于五階WENO格式、運動嵌套網(wǎng)格技術(shù)、Spalart-Allmaras湍流模型、隱式LU-SGS算法和并行策略的RANS方程求解方法。通過多種算例計算分析，結(jié)果表明：

1) 與試驗值相比，本文所發(fā)展的方法在旋翼渦流場的數(shù)值模擬中是有效的，能夠較為準確地計算旋翼非定常渦流場與氣動特性。

2) 與二階Jameson格式和Roe-MUSCL格式相比，五階WENO格式要表現(xiàn)出更高的精度與更低的數(shù)值耗散，能夠更好地保存渦結(jié)構(gòu)。

3) 在使用相同的計算網(wǎng)格和計算條件下，在獲得更高計算精度的同時，五階WENO格式的計算量稍大于JST格式的計算量，計算時間多10%～20%。

4) 五階WENO格式表現(xiàn)出了高階格式有助于渦流場中渦結(jié)構(gòu)細節(jié)捕捉，為后續(xù)的如旋翼槳/渦干擾、旋翼/機身干擾等復雜的渦流場特性提供了良好的數(shù)值基礎(chǔ)。

[1]LEISHMAN J G. Rotorcraft aeromechanics: Getting through the dip[J]. Journal of the American Helicopter Society, 2010, 55(1): 11001-1-11001-24.

[2]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.

[3]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.

[4]楊愛明, 喬志德. 用運動嵌套網(wǎng)格數(shù)值模擬旋翼前飛非定常流場[J]. 空氣動力學學報, 2000, 18(4): 427-433.

YANG A M, QIAO Z D. Numerical simulation of unsteady flow around a helicopter rotor in forward flight using moving overset grids[J]. Acta Aerodynamica Sinca, 2000, 18(4): 427-433 (in Chinese).

[5]ZHAO Q J, XU G H, ZHAO J G. New hybrid method for predicting the flowfields of helicopter rotors [J]. Journal of Aircraft, 2006, 43(2): 372-380.

[6]HARIHARAN N. High order accurate numerical convection of vortices across overset interfaces: AIAA-2005-1263[R]. Reston: AIAA, 2005.

[7]MAVRIPLIS D J, DARMOFAL D, KEYES D, et al. Petaflops opportunities for the NASA fundamental aeronautics program: AIAA-2007-4084[R]. Reston: AIAA, 2007.

[8]招啟軍, 徐國華. 使用高階逆風通量差分裂格式的懸停旋翼流場數(shù)值模擬[J]. 航空動力學報, 2005, 20(2): 186-191.

ZHAO Q J, XU G H. Numerical simulations for the flowfield of helicopter rotors in hovering based on high-order upwind flux-difference splitting scheme[J]. Journal of Aerospace Power, 2005, 20(2):186-191 (in Chinese).

[9]徐廣, 招啟軍, 王博, 等. 先進直升機旋翼懸停狀態(tài)氣動性能計算[J]. 航空學報, 2010, 31(9): 1723-1732.

XU G, ZHAO Q J, WANG B, et al. Prediction on aerodynamic performance of advanced helicopter rotor in hover[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(9): 1723-1732 (in Chinese).

[10]YU J, YAN C, JIANG Z H. A high resolution low dissipation hybrid scheme for compressible flows[J]. Chinese Journal of Aeronautics, 2011, 24(4): 417-424.

[11]李新亮. 高超聲速湍流直接數(shù)值模擬技術(shù)[J]. 航空學報, 2015, 36(1): 147-158.

LI X L. Direct numerical simulation techniques for hypersonic turbulent flows[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 147-158 (in Chinese).

[12]XU L, WENG P F. High order accurate and low dissipation method for unsteady compressible viscous flow computation on helicopter rotor in forward flight[J]. Journal of Computational Physics, 2014, 258(2): 470-488.

[13]HILGENSTOCK A. A fast method for the elliptic generation of three-dimensional grids with full boundary control[C]//Numerical Grid Generation in Computational Fluid Mechanics′ 88. Swansea, Wales: Pineridge Press, 1988: 137-146.

[14]MEAKIN R L. A new method for establishing intergrid communication among system of overset grids: AIAA-1991-1586[R]. Reston: AIAA, 1991.

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

[16]HARTEN A, HYMAN J M. Self adjusting grid methods for one-dimensional hyperbolic conservation laws[J]. Journal of Computational Physics, 1983, 50(2): 235-269.

[17]HARTEN A, ENGQUIST B, OSHER S, et al. Uniformly high order accurate essentially non-oscillatory schemes, Ⅲ[J]. Journal of Computational Physics, 1987, 71(2): 231-303.

[18]LIU X D, OSHER S, CHAN T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115(1): 200-212.

[19]JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1): 202-228.

[20]BLAZEK J. Computational fluid dynamics: principles and applications[M]. Netherlands: Elsevier, 2001: 204-208.

[21]SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston: AIAA, 1992.

[22]CHANDRA R. Parallel programming in OpenMP[M]. San Francisco: Morgan Kaufmann, 2001: 41-88.

[23]CARADONNA F X, TUNG C. Experimental and analytical studies of a model helicopter rotor in hover: NASA-TM-8122[R]. Cailifornia: Army Research and Technology Labs, 1980.

[24]LORBER P F. Aerodynamic results of a pressure-instrumented model rotor test at the DNW[J]. Journal of the American Helicopter Society, 1991, 36(4): 66-76.

[25]CARADONNA F X, LAUB G H, TUNG C. An experimental investigation of the parallel blade-vortex interaction: NASA-TM-86005[R]. California: Army Research and Technology Labs., 1984.

[26]HEFFERNAN R M, GAUBERT M. Structural and aerodynamic loads and performance measurements of an SA 349/2 helicopter with an advanced geometry rotor: NASA-TM-88370[R]. California: Army Research and Technology Labs., 1986.

印智昭男， 碩士研究生。主要研究方向： 直升機計算流體力學、 高精度CFD數(shù)值格式。

E-mail： yinzhizhao@nuaa.edu.cn

招啟軍男， 博士， 教授， 博士生導師。 主要研究方向：直升機計算流體力學、 直升機空氣動力學及流動控制、 氣動噪聲、 總體設(shè)計。

Tel: 025-84893753

E-mail： zhaoqijun@nuaa.edu.cn

王博男， 博士， 講師。主要研究方向： 直升機計算流體力學、 直升機空氣動力學、 新概念旋翼飛行器氣動設(shè)計等。

Tel: 025-84893753

E-mail： wangbo@nuaa.edu.cn

Numerical simulations for unsteady vortex flowfield of rotors based on high-order WENO scheme

YIN Zhizhao, ZHAO Qijun*, WANG Bo

National Key Laboratory of Science and Technology on Rotorcraft Aeromechanics, Nanjing University of Aeronautics and Astronautics, Nanjing210016, China

Aimed at improving the computational fluid dynamics (CFD) computation accuracy and reducing the numerical dissipation of the vortex flowfield of helicopter rotor, the highly-accurate CFD method based on fifth-order weighted essentially non-oscillatory (WENO) scheme has been established for numerical simulation on the unsteady flowfield of rotors. In this method, the moving-embedded grid technique has been employed in the numerical simulation, the Navier-Stokes equations have been taken as governing equations, and one equation Spalart-Allmaras turbulence model has been employed. In order to improve the computation accuracy of the rotor vortex formation and evolution, the discretization of convective fluxes and the time marching are conducted by Roe-WENO scheme and the highly efficient implicit LU-SGS (Lower-Upper Symmetric Gauss-Seidel) scheme respectively, and the parallel strategy is employed to accelerate the calculation as well. The computation of the left and right flow quantities is based on the fifth-order WENO scheme, which uses the appropriate weights to combine all the candidate stencils to obtain the fifth-order approximation. Finally aerodynamic characteristics and blade-tip vortex of the several rotors (C-T and UH-60A rotor in hover, C-T rotor in non-lifting forward flight and SA349/2 rotor in forward flight) have been simulated by the present method. Comparing numerical results with experimental data and calculated results using the JST scheme and Roe-MUSCL scheme, it is demonstrated that the present method with only 10%～20% increase of computational consumption has high accuracy and low numerical dissipation in simulating the aerodynamic characteristics and blade-tip vortex.

rotor; vortex flow; weighted essentially non-oscillatory scheme; Navier-Stokes equations; blade-tip vortex; moving-embedded grid

2015-08-18； Revised： 2015-09-28； Accepted： 2015-10-16； Published online： 2015-11-3009:22

s： National Natural Science Foundation of China (11272150); Priority Academic Program Development of Jiangsu Higher Education Institutions

. Tel.： 025-84893753E-mail: zhaoqijun@nuaa.edu.cn

2015-08-18； 退修日期： 2015-09-28； 錄用日期： 2015-10-16；

時間： 2015-11-3009:22

www.cnki.net/kcms/detail/11.1929.V.20151130.0922.002.html

國家自然科學基金 (11272150)； 江蘇高校優(yōu)勢學科建設(shè)工程基金

.Tel.： 025-84893753E-mail: zhaoqijun@nuaa.edu.cn

10.7527/S1000-6893.2015.0281

V211.3; V211.52

A

1000-6893(2016)08-2552-13

引用格式： 印智昭， 招啟軍， 王博． 基于高階WENO格式的旋翼非定常渦流場數(shù)值模擬[J]． 航空學報, 2016, 37(8): 2552-2564. YIN Z Z, ZHAO Q J, WANG B. Numerical simulations for unsteady vortex flowfield of rotors based on high-order WENO scheme[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2552-2564.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

URL： www.cnki.net/kcms/detail/11.1929.V.20151130.0922.002.html