龐超, 李猛, 高正紅

(西北工業(yè)大學 航空學院, 陜西 西安 710072)

旋翼作為直升機類飛行器的主要氣動部件,其氣動性能直接決定了整個飛行器的飛行性能。懸停和前飛是旋翼2個典型的工作狀態(tài),其中懸停性能(例如懸停效率,FoM)可以作為評價旋翼氣動性能的一個關鍵指標。隨著計算機算力的快速增長,直接從NS方程(空氣動力學中“第一性原理”)出發(fā)的計算流體力學(CFD)方法逐步被用在旋翼非定常流場的模擬中。朱正、盧叢玲等[1-2]利用運動嵌套網(wǎng)格方法,采用BL湍流模型及SA湍流模型對共軸剛性雙旋翼非定常懸停流場進行了模擬,分析了其流動干擾機理。張震宇等[3]利用多塊結構化滑移網(wǎng)格,采用k-ωSST兩方程湍流模型對剛性旋翼大前進比前飛狀態(tài)進行了模擬。美國航天航空學會(AIAA)自從2014年起每年召開一次旋翼懸停預測工作組會議(AIAA rotorcraft hover prediction workshop,HPW),相繼提供了3個旋翼標模[4-7]及相關的試驗數(shù)據(jù)作為統(tǒng)一的研究對象,包括S-76旋翼、PSP旋翼和HVAB旋翼。其中針對S-76旋翼[8-9]的討論主要集中在全湍流假設下,而HVAB旋翼模型目前缺乏正式的試驗數(shù)據(jù)作為支撐。PSP旋翼是當前HPW會議中唯一有槳葉表面試驗數(shù)據(jù)的計算標模,從2018年起國外研究人員開始對其進行初步的數(shù)值模擬研究[10]。在基于RANS方程的CFD方法中,轉(zhuǎn)捩模擬是非常重要的一部分,CFD2030愿景[11]以及HPW愿景中均將邊界層的轉(zhuǎn)捩模擬作為重要的研究內(nèi)容。對于固定翼類飛行器的二維、三維、高速、低速轉(zhuǎn)捩模擬已經(jīng)有相對成熟的方法和應用[12],然而這些轉(zhuǎn)捩模擬方法在旋翼模擬中的應用卻比較少見。其原因在于：一方面相比于固定翼模擬,旋翼模擬需要更大的計算代價,因為涉及非定常計算、尾跡捕捉等;另一方面轉(zhuǎn)捩模擬對計算網(wǎng)格尺度和網(wǎng)格量的要求要遠高于全湍流模擬。對于固定翼飛行器而言,在一些狀態(tài)下,其氣動、流場特性受到轉(zhuǎn)捩過程的影響，會有比較大的變化,所以對于這類飛行器在設計、評估過程中必須考慮轉(zhuǎn)捩的影響。而對于旋翼來說,人們需要明確是否也存在湍流模擬和轉(zhuǎn)捩模擬結果差異大這一現(xiàn)象,以確定在設計、評估過程中采用合適的模型和計算網(wǎng)格。

本文采用基于嵌套網(wǎng)格的自研RANS求解器,對HPW提供的PSP旋翼標模進行了懸停流場的全湍流和轉(zhuǎn)捩模擬,并從積分氣動力、槳葉表面壓力分布、槳葉載荷分布、槳葉表面流線、尾跡流場等方面詳細對比討論了全湍流模擬與轉(zhuǎn)捩模擬的差異,為設計人員在評估時提供依據(jù)。同時也進一步驗證了自研求解器在轉(zhuǎn)捩模擬和旋翼模擬方面的可靠性。

1 數(shù)值方法

1.1 RANS方程

本文自研求解器[13-17]基于動態(tài)結構化嵌套網(wǎng)格,對雷諾平均NS方程(RANS方程)進行數(shù)值求解,可壓縮RANS方程如(1)式所示。

(1)

式中:ρ為空氣密度;ui是速度分量;p為壓強;E為總能;H為總焓。

對于方程中無黏通量采用Roe迎風格式配合WENO插值實現(xiàn)5階精度,對于方程中的黏性通量采用2階中心差分格式計算。使用近似因式分解方法進行方程的時間推進,在非定常計算方面,采用隱式雙時間步長方法,并在子迭代中使用多重網(wǎng)格和當?shù)貢r間步長方法加速收斂。

物面采用無滑移邊界條件,遠場采用基于黎曼不變量的無反射邊界條件。

γeff=max(γ,γsep)

(2)

Dγ=ρΩγ(1-Gonset)

(3)

Reθt=F(Iinf)F(λθ)

(4)

1.3 嵌套網(wǎng)格方法

相較于其他方法而言,嵌套網(wǎng)格方法在網(wǎng)格處理方面有2個額外的步驟,即挖洞(hole-cutting)和插值(interpolation)。其中,挖洞是將侵入物面內(nèi)的網(wǎng)格單元和網(wǎng)格點進行標記,從而在計算中排除這些單元和點。插值是對于嵌套網(wǎng)格塊邊界處的網(wǎng)格點提供相應的插值模板單元和插值系數(shù),使得流場信息在相互嵌套的網(wǎng)格塊之間進行傳遞。本文挖洞基于Objective X-ray方法,插值采用Inverse Map與模板游走相結合的方法,具體實現(xiàn)方法參考文獻[13]。

2 求解器驗證

對于文中的轉(zhuǎn)捩模型采用標模平板算例(S & K平板[20]、T3A平板[21-22])和S809翼型[23]進行驗證,對計算得到的表面摩阻因數(shù)、壓力系數(shù)等結果與相應的試驗結果進行比對,說明求解器轉(zhuǎn)捩模擬的能力。

3個算例的來流湍流度及所采用的計算網(wǎng)格規(guī)模在表1中列出,其中J為流向,K為物面法向。S809翼型計算網(wǎng)格在圖1中給出。

表1 驗證算例計算狀態(tài)和網(wǎng)格規(guī)模

對平板標模算例,對比了CFD計算與試驗測得的表面摩阻因數(shù),而對S809翼型算例,對比了翼型表面壓力系數(shù)。對比結果如圖2～3所示,3個驗證算例的計算結果展示出了求解器在轉(zhuǎn)捩模擬方面良好的預測能力和精度。

由3個作用點處的截面直徑，換算出該截面處的慣性矩I的均值分別為：1.91e4，8.94e3，5.14e3mm4。彈性模量均值為24.41MPa，帶入式(7)分別得到以下3個關系式，即

圖1 S809翼型網(wǎng)格 圖2 平板標模算例表面圖3 S809翼型表面壓力系數(shù)對比結果 摩阻因數(shù)對比結果

3 計算模型與網(wǎng)格

3.1 計算模型

計算模型采用HPW工作會議提供的旋翼標模算例-PSP旋翼,其槳葉平面及截面翼型形狀如圖4所示。此旋翼共包含4片槳葉,半徑為1.689 m,參考弦長為0.138 m,槳葉實度σ為0.103 3,在0.95R處有30°的后掠角,槳葉轉(zhuǎn)速為1 150 r/min,對應的槳尖馬赫數(shù)為0.58,基于槳葉平均弦長和槳尖速度的雷諾數(shù)為1.87×106。

圖4 PSP旋翼槳葉平面形狀

該旋翼的風洞試驗于2016年在NASA-LaRC旋翼測試中心[24]完成,由于試驗論文中并未提及相應的風洞試驗湍流度,在本文轉(zhuǎn)捩計算中,湍流度設為0.7%。

3.2 計算網(wǎng)格

HPW工作組公開發(fā)布了PSP旋翼標模槳葉的幾何文件和表面網(wǎng)格文件,本文利用上述表面網(wǎng)格文件生成了相應計算網(wǎng)格,包含對應于4片槳葉的貼體網(wǎng)格和笛卡爾背景網(wǎng)格。背景網(wǎng)格在槳葉附近進行加密處理?？偩W(wǎng)格點數(shù)為5.1×107,其中貼體網(wǎng)格點數(shù)為4.4×107。貼體網(wǎng)格中槳葉表面法向首層網(wǎng)格距離為8.9×10-7m。背景網(wǎng)格最密的網(wǎng)格尺度為1/10倍的槳葉弦長。遠場網(wǎng)格距離旋轉(zhuǎn)中心約為20倍的槳葉半徑長度。圖5給出了槳葉貼體網(wǎng)格的3個截面和最密背景網(wǎng)格的1個截面,以及相應的挖洞結果。

圖5 PSP旋翼標模嵌套網(wǎng)格

4 計算結果與討論

4.1 旋翼積分力

在懸停狀態(tài)下,旋翼有2個最重要的氣動力系數(shù),即拉力系數(shù)CT和扭矩系數(shù)CQ。按照HPW工作組的統(tǒng)一定義,(5)式給出了這2個系數(shù)的計算公式。

(5)

在本文計算中,重點對比同拉力系數(shù)下的扭矩系數(shù)和相應的懸停效率。所有的計算結果均為旋翼最后一個旋轉(zhuǎn)周期內(nèi)的平均值。表2和圖6給出了懸停效率的計算結果與試驗結果的對比,圖中試驗值的誤差限取自文獻[24]中給出的最小值。該結果驗證了本文所采用的求解器對于旋翼懸停氣動力系數(shù)計算的可靠性,轉(zhuǎn)捩模擬結果和全湍模擬結果計算誤差均在5%之內(nèi)。

表2 CFD預測懸停效率及誤差

圖6 CFD預測旋翼氣動性能與試驗值對比

對于懸停效率這一評價旋翼氣動效率的重要參數(shù)而言,在低拉力狀態(tài)下,轉(zhuǎn)捩模擬得到明顯優(yōu)于全湍模擬的結果(約高10%),但是在大拉力狀態(tài)下,這一優(yōu)勢逐漸縮小(約高4%)。

4.2 截面壓力分布及槳葉載荷分布

(6)

式中:Ω為轉(zhuǎn)速；r為當?shù)匕霃剑籧為當?shù)叵议L。

4.2.1 截面壓力分布

沿槳葉展向從30%R處開始至80%R處結束,截取了6組壓力分布,對比結果如圖7～8所示。其中圖7為小拉力狀態(tài)下結果,圖8為大拉力狀態(tài)下結果。

圖7 小拉力狀態(tài)槳葉截面壓力分布對比 圖8 大拉力狀態(tài)槳葉截面壓力分布對比

相比而言,全湍流模擬結果得到了比較平滑的槳葉表面壓力分布,而由于計入了流動轉(zhuǎn)捩過程,轉(zhuǎn)捩模擬得到的表面壓力分布在流動轉(zhuǎn)捩區(qū)域附近出現(xiàn)了明顯波動。在其余區(qū)域,兩者壓力分布相似且重合度高。

通過槳葉截面上下表面壓力的波動可以判斷出槳葉轉(zhuǎn)捩發(fā)生的區(qū)域,經(jīng)過對比與試驗[24]中給出轉(zhuǎn)捩位置保持一致。進一步證明了本文所采用的求解器在旋翼懸停自然轉(zhuǎn)捩模擬上的可靠性。在小拉力狀態(tài)下槳葉表面出現(xiàn)較大范圍的層流流動,而隨著旋翼拉力的增大,該層流范圍逐步減小,這與4.1節(jié)中旋翼積分扭矩系數(shù)變化趨勢保持一致。

4.2.2 槳葉載荷分布

圖9給出了不同拉力狀態(tài)下,全湍流以及轉(zhuǎn)捩模型計算得到的槳葉載荷系數(shù)在槳葉展向上的分布,其中拉力系數(shù)分布呈現(xiàn)出合理的類橢圓形分布。

圖9 槳葉載荷分布對比

雖然截面壓力分布在轉(zhuǎn)捩區(qū)域有明顯的波動,但是對于沿槳葉展向的拉力系數(shù),除了槳葉根部小拉力區(qū)域,2種模擬方法并沒有表現(xiàn)出很大的不同。而對于沿槳葉展向扭矩系數(shù)來說,由于轉(zhuǎn)捩模擬中部分區(qū)域存在明顯的層流流動,在這部分區(qū)域轉(zhuǎn)捩模型預測得到的扭矩系數(shù)要低于全湍模型預測的值。在大拉力狀態(tài)下,由于表面層流范圍的縮小,全湍模擬和轉(zhuǎn)捩模擬得到的扭矩系數(shù)分布也趨于一致。

4.3 槳葉表面流線及旋翼尾跡流場形態(tài)

4.3.1 槳葉表面流線

圖10～11分別顯示出不同拉力狀態(tài)下槳葉上、下表面的表面流線,可看出全湍狀態(tài)下槳葉表面流場變化平緩,但是轉(zhuǎn)捩模擬得到了明顯不同的槳葉表面流線,其表面流線變化明顯,在槳葉內(nèi)側上表面發(fā)生了明顯的流動分離現(xiàn)象,這一現(xiàn)象在全湍模擬中并沒有出現(xiàn)。表面流線的明顯差異也進一步說明了4.2.1節(jié)中2種不同模擬方法得到的差異明顯的槳葉截面壓力分布。

圖10 小拉力狀態(tài)槳葉表面流線對比

圖11 大拉力狀態(tài)槳葉表面流線對比

4.3.2 尾跡流場

通過Q判據(jù)等值面,圖12給出了不同拉力狀態(tài)下旋翼槳盤下方的尾跡槳尖渦分布,利用渦量大小在對應的Q等值面上進行相應的染色。根據(jù)上述計算結果,轉(zhuǎn)捩過程引發(fā)的槳葉表面分離對槳尖渦的影響較小,但是對于槳盤下方小于槳盤半徑處的尾跡結構有比較明顯的影響。同時,隨著拉力的增大,槳尖渦的結構形式出現(xiàn)了一些變化,尾跡流場中開始出現(xiàn)圍繞著槳尖渦的二次渦結構,這一現(xiàn)象在全湍和轉(zhuǎn)捩模擬中均有所表現(xiàn)。

圖12 旋翼尾跡流場對比

5 結 論

本文利用自研RANS求解器,對AIAA HPW工作組提供的PSP旋翼標模,在懸停狀態(tài)下,分別進行了全湍流模擬和轉(zhuǎn)捩模擬。

計算結果首先表明本文CFD模擬得到的旋翼拉力-懸停效率曲線與試驗重合度高,計算誤差在5%之內(nèi),并且轉(zhuǎn)捩模擬中,CFD方法能正確模擬出槳葉表面流動的轉(zhuǎn)捩位置。由此,證明了本文求解器在旋翼懸停流場模擬方面的可靠性。

其次,通過全湍與轉(zhuǎn)捩模擬結果的對比發(fā)現(xiàn),轉(zhuǎn)捩過程對槳葉截面壓力分布、沿槳葉展向扭矩分布以及槳葉表面流線均有顯著影響,而對沿槳葉展向拉力分布影響較弱。轉(zhuǎn)捩過程給槳葉截面壓力分布帶來了明顯的波動以及槳葉表面流場的明顯分離。相比于全湍結果,轉(zhuǎn)捩模擬在層流區(qū)域得到了更低的扭矩系數(shù),從而帶來了懸停效率4%～10%的增量。

最后,對于旋翼槳盤下方槳尖渦,轉(zhuǎn)捩過程沒有明顯影響,但是由于轉(zhuǎn)捩過程帶來的流動分離,對于內(nèi)側槳葉下方尾跡區(qū)域,轉(zhuǎn)捩過程帶來了明顯的影響。