12高精度旋翼非定常渦流場數(shù)值模擬方法史勇杰8

1、第二十八屆（2012）全國直升機年會論文高精度旋翼非定常渦流場數(shù)值模擬方法史勇杰 魏 鵬 徐國華 招啟軍（南京航空航天大學直升機旋翼動力學重點實驗室，南京，210016）摘 要：基于粘性渦方法和CFD方法，發(fā)展了一套適合于復雜旋翼渦流場的耦合拉格朗日-歐拉計算方法。方法中針對高雷諾數(shù)旋翼流場中槳尖渦的緊湊結(jié)構(gòu)特點，引入粘性渦方法建立了高分辨率的尾跡求解模型，以此模擬尾跡的畸變運動及其影響；為捕捉激波等細節(jié)流場特征，在需要考慮細節(jié)流場的槳葉近體區(qū)域采用了CFD方程對其進行求解；而兩者之間的信息交換分別采用集中渦源法和邊界修正法完成。應(yīng)用所建立的計算方法，對Caradonna-Tung模型旋翼懸停

2、狀態(tài)下的旋翼流場進行了模擬，通過分析旋翼尾跡特征及耦合邊界處流場特征，表明所建方法充分發(fā)揮了CFD方法和粘性渦方法各自的優(yōu)勢，在旋翼流場計算方面的具有獨特的優(yōu)勢。最后通過對比槳葉表面壓強系數(shù)分布，進一步驗證了該方法在旋翼流場分析中的有效性。關(guān)鍵詞：旋翼流場；旋翼尾跡；Euler 方程；離散渦方法；混合方法；直升機0 引言槳葉近體的復雜非線性流動和槳葉產(chǎn)生的強烈槳尖渦是旋翼非定常流場的重要特征，兩者相互耦合、相互影響，對直升機的氣動載荷、氣動噪聲以及動力學特性等有重要影響，關(guān)于它們的研究一直是直升機技術(shù)研究領(lǐng)域的一個重點。然而，同時準確預測槳葉附近的非線性、跨音速、可壓縮流動和旋翼尾跡的畸變運動

3、具有很大難度，這給旋翼氣動特性的深入研究提出了巨大挑戰(zhàn)。迄今為止，對于旋翼非定常渦流場的研究，主要采用兩類數(shù)值計算方法，即基于拉格朗日法的渦動力學方法1和基于歐拉法的計算流體力學（CFD）方法2。其中渦動力學方法開展較早，傳統(tǒng)尾跡線渦法在經(jīng)歷了固定尾跡、預定尾跡以及自由尾跡的發(fā)展歷程之后日臻成熟，目前基于槳尖渦卷起模型的自由尾跡方法已被廣泛應(yīng)用于各種理論分析及工程計算中。而近幾年出現(xiàn)的高分辨率粘性渦方法成為繼自由尾跡方法后的又一個重要發(fā)展趨勢，Opoku，Anusonti，Zhao3等人先后開展了基于拉格朗日描述的旋翼粘性離散渦方法方面的研究，通過對不同旋翼流場的計算表明，所發(fā)展的方法可以更準

4、確的預測尾跡渦的畸變運動及擴散效應(yīng)。然而，這些方法中均采用了線性氣動力模型及不可壓假設(shè)，因此對于槳葉近體非定常流動的模擬存在困難，無法滿足現(xiàn)代復雜外形槳葉分析需求。上世紀80年代以來， CFD方法開始應(yīng)用于旋翼氣動問題的研究中，該方法在歐拉描述下通過數(shù)值方法求解Navier-Stokes方程，可以精確地捕捉槳葉表面氣動力的細致變化，以及前行側(cè)槳葉的強壓縮和后行側(cè)的反流、動態(tài)失速等復雜旋翼擾流特性。然而該方法在采用有限體積或有限差分法進行控制方程求解時，離散格式中固存的數(shù)值耗散會使得空間流場的渦出現(xiàn)非物理耗散，從而導致方法在無法準確捕捉旋翼尾跡。盡管很多學者已嘗試采用高階格式4、自適應(yīng)網(wǎng)格5等技

5、術(shù)消除這一缺陷并取得了較好的計算結(jié)果，但在現(xiàn)有計算資源條件下，尾跡的捕捉能力仍顯不足?？紤]到動力學方法與CFD方法各自的優(yōu)勢，許多學者將兩者相結(jié)合發(fā)展了混合計算方法1。早期的混合方法多基于線渦法和CFD方法，Berezin和Sankar等7基于CFD/全位勢/自由尾跡建立了一套混合方法，而Bader等采用場速度方法直接將CFD與自由尾跡耦合在了一起，計算結(jié)果表明，這些方法均可有效的增加尾跡捕捉能力，從而提高計算精度，然而受自由尾跡方法本身缺陷限制，它們的通用性并不好。近兩年Stock6，Zhao7等初步將高分辨率渦方法與CFD方法結(jié)合，發(fā)展了具有高精度的耦合方法，從已發(fā)表的文獻看，耦合方法在尾

6、跡捕捉和旋翼氣動性能預測等方面都要好于傳統(tǒng)的CFD方法，而且計算量也有顯著降低。目前國內(nèi)盡管已有不少學者開展了自由尾跡/CFD混合方法研究8, 9，但在粘性渦/CFD耦合方法方面尚未有相關(guān)研究。為此，本文在先前研究工作的基礎(chǔ)上，進一步將粘性渦方法與CFD方法相結(jié)合，建立了適合于復雜旋翼渦流場計算的耦合拉格朗日-歐拉計算方法。然后，采用Caradonna-Tung試驗旋翼為算例對所建方法的有效性進行了分析及驗證，結(jié)果表明該方法充分結(jié)合了兩種方法各自的優(yōu)勢，無論是在旋翼尾跡捕捉方面還是在槳葉近體流場細節(jié)模擬方面都是非常有效的。1 計算方法圖1給出了本文耦合拉格朗日-歐拉法的示意圖，該機算方法將旋翼

7、氣動問題分為尾跡預測和槳葉近體流場捕捉兩部分：在離槳葉較遠區(qū)域，尾跡的影響占據(jù)著主導地位，故采用基于拉格朗日描述的粘性渦方法進行尾跡預測；而圍繞槳葉的旋翼近體流場則應(yīng)用基于歐拉描述的CFD方法進行求解，以捕捉該流場區(qū)域的各種非線性特征以及近尾跡影響。圖 1 耦合拉格朗日-歐拉法示意圖1.1 旋翼尾跡求解模型直升機飛行時，旋翼槳葉拖出的尾渦發(fā)生卷起現(xiàn)象并形成槳尖渦，這些結(jié)構(gòu)緊湊且強度很高的槳尖渦堆積在槳盤下（后）方，成為旋翼尾跡的主要組成部分。拉格朗日粘性渦方法將該尾渦系離散為一系列無序渦元，并在拉格朗日描述體系進行數(shù)值求解來預測尾跡的卷起、畸變及耗散等過程。不考慮槳葉附近流場，直升機流場問題可

8、認為不可壓問題，離散渦元隨時間的變化由運動學方程和渦量動力學方程進行求解 式中，下標為渦元索引號，表示位置，為環(huán)量值，為當?shù)厮俣?為渦誘導速度，為時間，為運動粘性系數(shù)。式(1)中，渦量動力學方程的右邊第一項表示流場的速度梯度引起的渦元的伸縮與彎曲，其中的自誘導速度及速度梯度項由Biot-Savart定理計算，并引入快速多極子加速算法(FMM, Fast Multipole Method)10對該計算過程進行加速。而第二項表示的是渦量的粘性擴散效應(yīng)，本方法采用粒子強度交換法(PSE, Particle Strength Exchange)11對其進行求解。受槳尖渦誘導作用，旋翼流場存在著的強烈剪

9、切流導致長時間數(shù)值計算后離散渦元會出現(xiàn)嚴重的疏密失調(diào)，導致渦元間無法滿足保證計算精度的重疊性要求，從而致使計算精度下降。為此，本文采用文獻12的渦元重布方法來調(diào)整渦元的疏密度，該方法根據(jù)當?shù)厮俣葢?yīng)變張量判斷局部流場的拉伸狀況，當渦元所在位置的局部拉伸達到特定值時，將渦元分裂成兩個相同強度的渦元填補由于拉伸導致的“空白”。然而，隨著新渦元從槳葉脫落以及舊渦元的分裂，計算過程中渦元數(shù)目會迅速增長，導致計算效率降低；為合理控制因渦元數(shù)目增長導致的效率問題，一方面可以將遠離旋翼的渦元刪除，另一方面則需要對距離很近且方向較為一致（渦矢量夾角小于特定值）的渦元進行合并處理。1.2 槳葉流場計算模型為捕捉槳

10、葉近體空間的非定常流動，采用以絕對物理量為變量的守恒形式的雷諾平均N-S方程對其進行數(shù)值求解。在慣性系(x,y,z)下，方程表示為 式中，守恒變量，為慣性系下三個方向上的速度，分別表示壓強、密度和能量；為控制體單元體積，為面積，分別為對流通量喝粘性通量，為粘性項邏輯開關(guān)。本文旨在研究所建方法的耦合有效性，著重于兩計算域在交互過程中的信息傳遞，考慮到粘性項對此沒有直接影響，在實際計算過程中一方面出于方便對數(shù)值耗散進行研究的緣故，另一方面出于效率方面的原因，并未打開粘性項邏輯開關(guān)。計算過程中所有的物理量均以懸停槳尖速度、密度和槳葉弦長作為參考值而無量綱化?？刂品匠痰目臻g離散采用二階中心方法，而時間

11、推進上采用雙時間方法。2 計算域間的信息交換對于本文所建方法，拉格朗日計算域與歐拉計算域之間的相互信息交換至關(guān)重要，合理的信息交換策略無論對計算精度還是對計算效率都有幫助，因此需要對其進行仔細處理。2.1 歐拉計算域至拉格朗日計算域CFD計算域至粘性渦計算域的信息傳遞，可采用兩種方式：集中渦源法和分布渦方法。后者直接將CFD計算域計算的渦量分布傳遞給拉格朗日計算域，物理意義更為明確，但卻容易將CFD計算域的數(shù)值耗散引入到拉格朗日計算域，同時也存在著計算工作量大的問題。為此，本文采用了集中渦源法，即先由CFD方法計算得到槳葉展向載荷分布，再根據(jù)庫塔如科夫斯基（Kutta-Joukowski）環(huán)量

12、/升力定理來獲得槳葉展向附著渦環(huán)量分布 其中是槳葉剖面升力系數(shù)，為附著渦環(huán)量值。槳葉附著渦隨空間及時間的變化是槳葉脫落尾渦的主要原因，其中附著渦隨空間變化脫落的尾渦為尾隨渦，隨時間變化脫落的槳尖渦為脫體渦。因此在求解得到槳葉附著渦后，拉格朗日計算域新生渦元的強度即可由下面公式確定 式中，表示新生渦環(huán)量值，為槳葉剖面相對氣流速度，包括來流速度、旋翼旋轉(zhuǎn)速度以及槳葉揮舞運動速度。2.2 拉格朗日計算域至歐拉計算域粘性渦計算域至CFD計算域的信息傳遞，常用的有幾種方式，第一種方法是采用渦誘導速度直接修正槳葉當?shù)赜牵坏诙N是采用廣義網(wǎng)格速度法將影響施加在全體網(wǎng)格格點中，即將渦線產(chǎn)生的誘導速度等效成網(wǎng)

13、格畸變；第三種方式直接通過邊界條件施加于CFD計算域外邊界上。其中，第一種方式盡管效率很高，然而在迎角修正過程中會引入誤差，精度并不高；第二種方法尾跡貢獻施加在所有網(wǎng)格點上，槳尖渦對流場的強烈影響可以直接體現(xiàn)出來，但這也導致其效率很低；第三種方法較前兩種方法更符合物理意義，在控制好網(wǎng)格質(zhì)量的前提下可以得出滿意結(jié)果（參考文獻），而且僅需計算外邊界處的誘導速度，效率較高，因此本文采用了第三種方式。計算過程中，所需邊界條件除速度外，還需要密度和能量，本文假定CFD計算域大小恰當，槳葉網(wǎng)格外邊界以外壓縮性已經(jīng)可以忽略，因此密度等于遠場處密度，而能量邊界條件為 3 方法驗證及分析為驗證本文方法對懸停流場

14、的計算能力，采用了有實驗結(jié)果可供對比的“Caradonna & Tung”試驗旋翼做為算例進行計算，為避免氣彈作用對槳葉氣動特性的影響，旋翼槳葉具有很高的設(shè)計揮舞剛度，因此被廣泛用于旋翼流場計算的對比研究。該槳葉展弦比為6，弦長0.1905，翼型為NACA 0012，無扭轉(zhuǎn)和尖削。試驗進行了不同總距、不同槳尖馬赫數(shù)的多個算例，本節(jié)中若無特別說明，選取總距8度、槳尖馬赫數(shù)0.439的狀態(tài)作為默認算例。3.1 黏性渦方法與CFD方法算例分析作為各自方法的驗證，以及為本文耦合方法提供算例對比，首先分別采用黏性渦方法和CFD方法獨立的對默認算例進行了計算。對于拉格朗日渦方法，影響其計算精度及計

15、算時長的主要參數(shù)有光滑參數(shù)（為重疊因子，為最小渦元尺寸），槳葉分段數(shù)，方位角步進步長，計算域控制參數(shù)（距槳轂中心大于該參數(shù)的渦元將被刪除），以及渦元合并距離系數(shù)和渦元合并角度閥值。本文各參數(shù)取值為：、，計算結(jié)束時渦元數(shù)約為3萬多。而對于CFD計算方法，網(wǎng)格數(shù)量是影響精度和效率的主要因素。在槳葉網(wǎng)格方面，為了方便對比，包括下一節(jié)耦合方法計算在內(nèi)所有計算均采用了同一套網(wǎng)格，即213×37×89的CO型貼體網(wǎng)格，網(wǎng)格外邊界距離槳葉物面約1倍弦長。另外，本節(jié)的全CFD求解模型采用如圖2所示的嵌套網(wǎng)格系統(tǒng)，其背景網(wǎng)格采用關(guān)鍵部位加密處理的笛卡爾網(wǎng)格，為了對比研究背景網(wǎng)格對槳尖渦數(shù)值耗

16、散的影響，分別采用137×107×137和118×85×118兩套背景網(wǎng)格進行了計算。圖 2 用于全CFD計算的嵌套網(wǎng)格系統(tǒng)圖3首先給出了黏性渦方法計算結(jié)束時的渦元分布圖，其中渦元的體積表征環(huán)量值的大小，可以看出，受下方尾跡的誘導作用，旋翼槳尖渦在離開槳葉后迅速收縮。隨著壽命角的增大，收縮運動逐漸減緩，而槳尖渦部分區(qū)域出現(xiàn)非規(guī)則畸變，這導致大約在兩圈以后相互靠近的槳尖渦間開始出現(xiàn)嚴重干擾現(xiàn)象，槳尖渦間的這種干擾現(xiàn)象導致旋翼渦系發(fā)生破壞并最終形成紊流狀態(tài)，文獻13對此現(xiàn)象也進行了分析，且給出了本文一致的結(jié)論。(a) 黏性渦方法(b) CFD方法圖3 流場空

17、間渦量等值面圖旋翼槳尖渦的收縮畸變是懸停旋翼流場的重要特征，準確的捕捉該畸變對于旋翼氣動分析具有重要意義。圖3給出了本文耦合方法計算得到的槳尖渦徑向位置及軸向位置隨壽命角的變化趨勢，由圖可見，槳尖渦軸向位移的計算值在小壽命角時變化緩慢，當壽命角達到180度時，槳尖渦受另外一片槳葉槳尖渦的誘導影響，軸向位移明顯加快，這與實際旋翼尾跡運動特性相符合，也與試驗值吻合很好，同時可以看出徑向收縮量也與試驗值得到了很好吻合，這為進一步開展拉格朗日-歐拉法的研究提供了必要的前提保證。圖4 計算與試驗的槳尖渦徑向位置及軸向位置對比圖3.2 耦合拉格朗日-歐拉法的驗證與分析圖5為該算例狀態(tài)下的渦量圖，為顯示直觀

18、，CFD計算域的渦量用槳葉網(wǎng)格外表面的渦量云圖表示，拉格朗日計算域的渦量采用渦量等值面顯示。從圖中可以看出，兩個計算域的的渦量集中區(qū)域完全吻合，當拉格朗日計算域的槳尖渦量穿入CFD計算域時，渦量信息通過CFD計算域外邊界的速度邊界條件成功的完成了傳遞，渦量的這種信息交換使得本文方法可非常有效的改善CFD計算域的數(shù)值耗散問題，同時也有助于增強槳渦干擾特性的捕捉。圖5 耦合拉格朗日-歐拉法計算得到的渦量圖圖6則給出了流線分布圖，其中有壓力云圖的部分為CFD計算域的計算結(jié)果，白色背景部分為拉格朗日渦方法的計算結(jié)果，圖中可以看出，除邊界處有強烈的渦作用的位置外，兩區(qū)域的流線平滑連接，具有很好的連續(xù)性，

19、表明拉格朗日計算域的流場信息通過邊界處信息交換有效的傳入到了CFD計算域。（a）r=0.4R(b) r=0.6R(c)r= 0.8R圖6 槳葉剖面流線計算值最后給出槳葉不同剖面處表面壓強系數(shù)分布曲線（圖7），可以看出，計算得到的壓強系數(shù)分布曲線與試驗值在各個剖面處都吻合很好，表明本文方法在低槳尖馬赫數(shù)下可以準確的計算槳葉細節(jié)氣動特性。而對于高槳尖馬赫數(shù)狀態(tài)，槳葉近體流場的壓縮效應(yīng)更大，為驗證本文方法的不可壓邊界條件假設(shè)對高槳尖馬赫數(shù)算例狀態(tài)的影響，又計算了槳尖馬赫數(shù)為0.877的算例狀態(tài)，圖8給出了該狀態(tài)下的槳葉表面的壓強系數(shù)分布，從圖中看出，本文建立的耦合方法可以捕捉到激波導致的壓強突變，并

20、與試驗值基本吻合，僅激波位置較試驗值有所偏后，這可能是由于本文CFD計算域所采用了歐拉求解器，而在無激波位置，計算值與試驗值吻合的則更好。(a) r/R=0.8(b) r/R=0.89(c) r/R=0.96圖7 槳葉表面壓強系數(shù)分布（Ma=0.439）(a) r/R=0.8(b) r/R=0.89(c) r/R=0.96圖8 槳葉表面壓強系數(shù)分布（Ma=0.877）4 結(jié)論本文基于粘性渦方法和CFD方法，發(fā)展了一套適合于復雜旋翼渦流場的耦合拉格朗日-歐拉計算方法。通過Caradonna-Tung模型旋翼懸停狀態(tài)下的旋翼流場進行了模擬，通過分析旋翼尾跡特征及耦合邊界處流場特征，表明所建方法充分

21、發(fā)揮了CFD方法和粘性渦方法各自的優(yōu)勢，在旋翼流場計算方面的具有獨特的優(yōu)勢。參 考 文 獻1Komerath N M, Smith M J, Tung C. A Review of Rotor Wake Physics and Modeling J. Journal of the American Helicopter Society, 2011, 56(2): 220061-2200619.2Strawn R C, Caradonna F X, Duque E P N. 30 years rotorcraft computational fluid dynamics research and

22、 development J. Journal of the American Helicopter Society, 2006, 51(1): 5-21.3Zhao J, He C. A Viscous Vortex Particle Model for Rotor Wake and Interference Analysis J. Journal of the American Helicopter Society, 2010, 55(1): 1-14.4Hariharan N, Sankar L. First-principles based high order methodologi

23、es for rotorcraft flowfield studies; proceedings of the Proceedings of 55th Annual Forum of AHS, F, 1999 C. American Helicopter Society, 1999.5Dietz M, Kessler M, Krämer E. Tip Vortex Conservation on a Helicopter Main Rotor Using Vortex-Adapted Chimera Grids J. AIAA Journal, 2007, 45(8): 2062-2

24、074.6Stock M J, Gharakhani A, Stone C P. Modeling Rotor Wakes with a Hybrid OVERFLOW-Vortex Method on a GPU Cluster; proceedings of the 28th AIAA Applied Aerodynamics Conference, Chicago, Illinois, F, 2010 C.7Zhao J, He C. A Hybrid Solver with Combined CFD and Viscous Vortex Particle Method; proceed

25、ings of the American Helicopter Society 67th Annual Forum, Virginia Beach, VA, F, 2011 C. American Helicopter Society.8Cao Y-H, Yu Z-Q, Su Y, et al. Combined Free Wake/CFD Methodology for Predicting Transonic Rotor Flow in Hover J. Chinese Journal of Aeronautics, 2002, 15(2): 65-71.9Shi Y, Zhao Q, F

26、an F, et al. A New Single-blade Based Hybrid CFD Method for Hovering and Forward-flight Rotor Computation J. Chinese Journal of Aeronautics, 2011, 24(2): 127-135.10Greengard L, Rokhlin V, 1997.11Mas-Gallic S. Contribution à l'analyse numérique des méthodes particulaires D, 1987.12

27、Mansfield J R, Knio O M, Meneveau C. Dynamic LES of Colliding Vortex Rings Using a 3D Vortex Method J. Journal of Computational Physics, 1999, 152(1): 305-345.13Leishman J G. Principles of helicopter aerodynamics M. 2nd ed. Cambridge ; New York: Cambridge University Press, 2006.High Precision Numeri

28、cal Method for Simulating Rotor Unsteady Votex Flow FieldSHI Yong-jie WEI Peng XU Guo-hua ZHAO Qi-jun(National Key Laboratory of Rotorcraft Aeromechanics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016)Abstract: A new coupled Eulerian-Lagrangian numerical method is developed for complex rotor vortex flowfiled simulation based on the viscous vortex method and the CFD method. Considering the compact structure of the tip vortex in the high Reynolds rotor fl