基于混合雷諾平均/高精度隱式大渦模擬方法的高升力體氣動噪聲模擬*

葛明明 王圣業 王光學2) 鄧小剛?

1) (國防科技大學空天科學學院,長沙 410073)

2) (中山大學物理學院,廣州 510275)

發展了基于七階精度混合型耗散緊致格式(HDCS)的混合雷諾平均(RANS)/高精度隱式大渦模擬(HILES)模型(HRILES),并結合Ffowcs Williams-Hawkings (FWH)聲比擬方法對30P30N 高升力體氣動噪聲問題進行了模擬.首先對雷諾數 Red=4.3×104 的單圓柱繞流算例開展驗證,并與傳統的延遲分離渦模擬(DDES) 模型進行對比.結果表明HRILES模型具有對亞臨界態尾跡區轉捩流動的模擬能力,平均阻力系數與阻力均方根值和實驗結果符合更好,結合FWH 聲比擬方法得到了合理的遠場聲壓級(SPL)的功率譜密度(PSD)分布.然后將其應用于30P30N 高升力體氣動噪聲算例模擬,結果表明HRILES模型準確預測縫翼凹腔剪切層各站位的平均速度、渦量和湍動能分布,壁面脈動壓力譜分布與實驗符合較好,近、遠場噪聲頻譜準確預測了縫翼低頻窄帶噪聲,并得到了合理的噪聲輻射指向性分布.

1 引 言

隨著航空運輸業的快速發展,航空器噪聲問題受到廣泛關注.歐盟早在2000年就在European Visions 2020中提出降低航空器噪聲50%的目標計劃[1].美國隨后也在AST (advanced subsonic transport)中提出安靜飛行計劃.伴隨著大涵道比發動機降噪技術的發展,機體噪聲的研究越來越受到重視,而在機體噪聲中高升力體部件貢獻了起降期間的大部分聲能.一般認為高升力體前緣縫翼的噪聲輻射是主要聲源[2],縫翼噪聲主要部分是寬頻噪聲,同時寬頻中會疊加多個窄帶峰值(NBP),此類高雷諾數氣動噪聲問題對數值模擬能力提出較高要求,需要發展兼顧計算精度與計算效率的計算模型與算法[3-5].

根據NASA 2030報告,大渦模擬方法將是航空工業下一代主流的湍流模擬方法[6].傳統的大渦模擬(LES)模型存在一些尚未被解決的問題,如亞格子項存在被截斷誤差掩蓋的可能性、空間濾波中的不確定度、在高雷諾數下的壁面流動問題中難以建立亞格子模型等[7].隱式大渦模擬作為一種簡單易用方法,被應用于各類問題計算,其中包括各向同性耗散湍流[8]、 自由剪切流[9]、 壁面流動[10],超聲速流動[11-13]、流動主動控制[14]等.高精度隱式大渦模擬(HILES)是基于HDCS-E8T7格式發展的一種隱式大渦模擬方法,利用格式的截段誤差代替顯式的亞格子模型,同時滿足幾何守恒律[15],已經成功地應用于各類復雜外形流動與噪聲預測[16-18],但對于高雷諾數壁面流動模擬,計算資源的需求仍然較大.

目前湍流模型實際應用仍然以非定常雷諾平均方程(URANS)、混合RANS/LES模型以及壁面建模的大渦模擬(WMLES)為主[19].混合RANS/LES類方法的發展通常分為分區與不分區兩種思路,分區方法人為將計算域設置為不同的區域采用不同的模型計算,流場信息通過交界面進行交換.不分區方法對全域統一求解,利用混合函數或者改變源項實現不同模型之間的切換.前者的關鍵在于處理交界面邊界條件,而后者在于交界面不確定性帶來的灰區問題.不分區方法中最著名的是Spalart等[20]發展的分離渦模擬(DES)模型,經典的DES模型引入網格尺度來修改Spalart-Allmaras模型(SA模型)破壞項中壁面距離,從而實現RANS向LES的轉換,在此基礎上改進發展出延遲分離渦模擬(DDES)、改進的延遲分離渦模擬(IDDES)等方法[21,22].Nichols[23]提出了構造混合模型的另一種思路,根據局部網格尺寸和湍流長度尺度從RANS平滑過渡到LES,稱之為多尺度模型.隨后將其應用于三維圓柱和武器填埋倉的流動模擬中,結果表明該混合方法在湍流尺度與網格尺度的比率接近2時,實現了較為合理的模擬結果[24].

本文利用Nichols[23]提出的多尺度模型構造思想,結合剪切應力輸運(SST)模型和HILES方法,希望實現對于高雷諾數流動的絕大部分區域使用HILES精確求解大尺度流動,僅邊界層內層利用SST 模型充當壁面模型.本文第2節主要介紹混合節點半節點型耗散緊致格式(HDCS-E8T7)與HRILES 方法,第3節模擬亞臨界態的圓柱流動以驗證HRILES方法,第4節采用HRILES和IDDES模型針對高升力體噪聲問題進行模擬,并將結果與實驗數據進行對比分析.

2 數值方法

2.1 高精度數值格式

本文所有的算例都采用七階混合型耗散緊致差分格式HDCS-E8T7,該格式采用顯式八階差分格式對對流項進行空間離散,通過在中心插值的基礎上增加奇數階導數來調高數值耗散,抑制數值色散,具體推導過程見文獻[25].對流項差分格式表示為

其中 α 為可調節耗散參數,根據色散關系保持格式優化取值為0.3085.黏性項采用八階中心差分格式加線性插格式進行計算.時間推進采用二階隱式雙時間步方法以節約計算資源.HDCS 格式已經應用于陣列圓柱、串聯圓柱-翼型、噴流噪聲[16-18]、三角翼[26]等典型算例中,為復雜外形氣動聲學問題模擬奠了定基礎.

2.2 HRILES模型

由于RANS方程與LES方程的相似性,可以通過過渡函數實現不同模擬方法的切換.Piomelli等[27]提出了基于SST的多尺度模型,對Red=8×106的圓柱繞流進行了模擬,認為多尺度模型相比于DES模型的優點主要在于能夠基于當地網格尺度平滑地從RANS方法過渡到LES方法.本文采用Morgan優化的過渡函數：

其中

(5)式中的 LT和 LG分別代表湍流長度尺度和網格長度尺度,原始文獻中定義如下：

本文計算中網格長度尺度參照Spalart等[21]在IDDES中的建議設置為

可以計算出湍流黏性系數 νt為

對于隱式大渦模擬,可認為亞格子黏性由數值格式隱式提供,因此湍流黏性系數可以簡單表示為νt=νtRANSfd.實際計算中,RANS模型只在邊界層內層啟動,大部分邊界層流動采用隱式大渦模擬,根據Wang等[26]的綜述,這類混合方法也可以看作為WMLES方法的一種.混合模型構造思路是想通過尺度判據將邊界層內層采用RANS 模型求解,它最顯著的優點是可以推廣到任意的RANS模型和SGS模型的組合.SST作為廣泛應用的工程湍流模型在高雷諾數復雜外形流動中表現良好,HILES方法具有高精度高效率的特點,HRILES模型可以結合二者的優勢.

2.3 噪聲計算方法

對于氣動聲學問題,遠場噪聲計算通常是采用聲比擬方法對近場計算得到的聲源面的時間序列數據進行積分,其中FWH方法最為常用.聲源面可以采用壁面或者空間任意可穿透面,區別在于后者包含了可能存在的四極子體聲源.由于本文算例屬于低亞音速流動,四極子聲源占比可以忽略不計,故采用固體壁面作為聲源面,采用Francescantonio提出的第二KFWH方程,積分求解遠場噪聲[28]：

其中

(11)式中的下標ret代表積分項對應的是推遲時間的數據,S代表積分聲源面.方程(10)等式右邊前兩項代表厚度噪聲(單極子聲源),后兩項代表載荷噪聲(偶極子聲源).對于二維構型的算例,由于網格量的限制,往往展向長度小于實驗長度,這會影響遠場噪聲結果,需要對聲學結果進行校正.本文采用緊致聲源校正常用的Kato公式[29]：

其中SPLexp為校正后聲壓級,SPLsim為數值積分得到的待校正聲壓級,Lsim,Lc和 Lexp分別代表數值模擬展向長度、流動展向相干長度、實驗展向長度.

3 單圓柱繞流

3.1 計算設置

采用HRILES模型對單圓柱算例對進行模擬,基于圓柱直徑d的雷諾數為 Red=43000 ,來流馬赫數0.21.該條件下流動處于亞臨界態,包括邊界層分離,尾跡區轉捩以及大尺度渦結構等復雜現象[30],對復雜流動結構的模擬精度直接影響阻力計算結果.計算采用O型網格進行,其中周向和法向分布為 181×181 ,展向均布31個點,壁面法向第一層網格無量綱高度 h/d=1.0×10-5.時間推進采用雙時間步法,真實時間步長 dt=1×10-6s.計算100 個對流周期得到近似擬序結構,又計算了200個周期用于統計平均與聲源面信息收集.

3.2 流動結果

表1中給出了流動結果的統計數據,并與經典SST-DDES模型的結果和Szepessy和Bearman[31]的實驗結果進行了對比.從阻力系數來看,相比SST-DDES方法的結果,HRILES方法得到的平均阻力系數 CD,ave和阻力均方根值CD,rms都更加接近實驗值.升力系數的斯特勞哈爾數Sr表征圓柱尾部的渦脫落頻率,本文計算得到的結果與Seo等[32]的LES 計算結果更加接近,采用DDES模型會高估壁面的渦脫落頻率.θsap代表分離起始位置角度,θ 計算以來流方向為起點,順時針為正,兩種模型計算得到的分離位置接近.

圖1展示了HRILES模型計算得到的壁面平均壓力系數 Cp分布結果,與實驗結果和SST-DDES模型的計算結果進行了對比.HRILES模型計算出的整個背風區壓力系數分布都要比SST-DDES模型結果更接近實驗值.從圖2給出的統計平均得到的流場流線分布來看,HRILES 模型與SSTDDES方法的結果差別較大,前者能捕捉到背風區尾跡的中等尺度的分離泡,而SST-DDES計算結果只得到了大尺度尾部回流區.

表1 單圓柱算例流動參數統計結果Table 1.Statistical results of aerodynamic coefficients for the single cylinder.

圖1 圓柱表面平均壓力系數分布Fig.1.Mean wall pressure coefficient distribution of the rod.

圖2 流線分布 (a) HRILES; (b) SST-DDESFig.2.Distribution of streamlines：(a) HRILES; (b) SSTDDES.

3.3 噪聲結果

Jacob等[33]在德國宇航中心(DLR)的聲學風洞進行了相似條件下的實驗研究,測量得到圓柱正上方,距離 r=185d 處觀測點的聲壓級功率譜密度分布.本文聲學結果計算采用FHW方程對壁面瞬態流動數據積分得到,并利用Kato公式修正展向長度差異的影響.每個非定常時間步采集壁面數據,共得到約3萬個采樣點,這些數據被分為7個窗,數據重疊率 50%,最終功率譜密度的分布通過平均7個窗的結果得到.圖3對比了兩種模型得到的對應實驗遠場觀測點的聲壓級功率譜分布,r代表觀測點至圓柱中心的距離.結果表明HRILES方法能夠準確捕捉尖頻噪聲,并且主頻幅值相較SST-DDES模型更接近實驗結果.

圖3 遠場 θ=90° ,r=180d 觀測點聲壓級功率譜密度Fig.3.Farfield acoustic result of the rod：PSD at (θ=90°,r=180d).

4 高升力體

4.1 計算設置

算例幾何外形與機體噪聲會議中給出的30P30N縫翼噪聲標準算例一致,氣動弦長C為18 in (1 in=2.54 cm).縫翼與襟翼的長度分別為0.15C與 0.3C.計算采用等效飛行條件,襟翼縫翼張開角度均為 30°,來流馬赫數0.17,基于弦長的雷諾數 ReC=1.7×106,攻角 5.5°.遠場邊界采用特征邊界條件以消除邊界聲波反射影響,展向邊界采用周期性邊界條件.

計算網格見圖4,采用日本宇航中心(JAXA)提供的多塊對接結構網格,計算域沿壁面法向方向延伸100個弦長.JAXA關于網格具體介紹參考文獻[34].本文為驗證方法,壁面法向網格雷諾數(y+)最大值約1.5.二維網格總量70000,展向長度設置為2 in,均布61 個網格點,空間網格總量為4200000.

圖4 30P30N計算網格Fig.4.Mesh of 30P30N airfoil.

4.2 流動結果

圖5展示了IDDES模型和HRILES模型計算得到的瞬態Q判據等值面(QC/U∞=5000),U∞為來流速度.Q等值面通過瞬態無量綱展向渦量 ωz進行著色,參考量為 U∞/C.結果清晰展示了縫翼凹腔渦流的尺度范圍和剪切層發展軌跡.縫翼下側凹腔內外壁面流動在過尖緣后混合,形成剪切流動,并迅速轉捩為湍流.湍流結構在空間充分發展后,最終在凹腔上壁面再附,再附點是數值模擬的重要指標.流動再附后部分流動向下游發展,部分沿上壁面回流進入凹腔形成渦流.IDDES和HRILES 模型模擬出的剪切層是相似的,但是HRILES模型能分辨出更精細的渦結構,特別是在縫翼凹腔內部.

圖6展示了x-y截面平均展向渦量 ωz,ave云圖,計算結果沿展向進行了平均.兩種模型計算結果分布相似,IDDES模型得到的剪切層渦量峰值比HIRLES 的更大.Pascioni[35]的實驗測量了凹腔剪切層七個不同站位的流動參數分布,在圖6(b)中分別對應 L1- L7.圖6(b)中 P1- P6代表高升力體表面非定常壓力脈動數據測點.圖7的x-y截面平均流線分布展示了高升力體流動中典型的流動特征,即縫翼和主翼凹腔中的兩個回流區.本文計算狀態條件下背風面均未出現流動分離現象.

圖5 QC/U∞=5000 等值面 (a) IDDES; (b) HRILESFig.5.The isosurfaces of the Q-criterion (QC/U∞=5000)：(a) IDDES; (b) HRILES.

圖6 平均展向渦量云圖 (a) IDDES; (b) HRILESFig.6.Contours of meanmean spanwise vorticity：(a) IDDES; (b) HRILES.

圖7 平均流線分布 (a) 縫翼; (b) 襟翼Fig.7.Distribution of streamlines：(a) Slat; (b) flap.

圖8給出了翼型表面平均壓力系數的分布結果與JAXA的實驗結果的對比[36],統計了0.05-0.10 s時間段的非定常數據,并沿展向進行了平均.HRILES模型結果大部分位置與JAXA的實驗結果符合得很好,縫翼上表面負壓略微偏低.IDDES模型的結果縫翼上表面明顯低估了背風區吸力效應.這與前人的模擬結果規律相似,Terracol等[34]采用DDES和分區DES(ZDES)模擬的結果也低估了背風面吸氣效應,他指出可能是由于縫翼上表面邊界層在當前雷諾數條件下是層流,使用全RANS 模型模擬引入的誤差; Gao等[37]采用通量重構(FR)方法進行的模擬認為是其計算幾何模型沒有考慮縫翼尾緣厚度造成的; Zhang等[38]的模擬結果也表明縫翼上表面壓力系數脈動的網格敏感性強于下表面.由于本文采用的也是JAXA的網格,幾何構型相同,因此得到與其類似的結果應該是由于IDDES 模型在背風面邊界層大部分區域也采用了RANS模擬.圖9給出了縫翼表面壓力系數脈動均方根的分布,并對比了實驗給出的 P2- P6五個點的數據.兩種模型的結果得到的再附點的位置(壓力脈動峰值處)接近,但HRILES模型得到的均方根值比IDDES模型更接近實驗值.

圖8 壁面壓力系數分布Fig.8.Distribution of wall pressure coefficient.

圖9 縫翼表面壓力系數脈動均方根分布Fig.9.RMS of the fluctuating pressure coefficient on the surface of the slat.

圖10 各個站位的平均速度分布Fig.10.Mean velocity magnitudes along the seven lines across.

圖11給出了縫翼凹腔剪切層對應站位的平均展向渦量 ωz,ave的結果.在剪切層起始位置,兩種模型的計算結果都要高于實驗數據,而隨著剪切層向下游充分發展,計算結果與實驗結果逐漸接近.近壁面剪切層的計算誤差可能是由于模擬剪切層轉捩現象滯后,剪切層未充分摻混.

圖12給出了縫翼凹腔剪切層對應站位的平均湍動能TKEave結果,整體上HRILES模型的結果明顯更加接近實驗數據.初始站位兩種模型結果都偏低,這意味著計算出的剪切層相比實驗存在轉捩推遲現象.不同的是HRLES模型在 L2站位后與實驗數據符合較好,而IDDES模型的湍動能始終偏低直至剪切層下游 L6站位,這也說明HRILES模型在空間剪切層發展的模擬上更具優勢.

圖13給出了壁面非定常壓力脈動功率譜密度分布.P1點位于主翼下表面,處于流場的穩態區域,此處的壓力譜可以近似表征噪聲結果.實驗數據中存在兩個窄帶噪聲峰值,對應的頻率分別是1330與1960 Hz,前人的研究認為這種低頻窄帶噪聲來源于凹腔內部的渦-聲耦合自激勵反饋機制.兩種模型的脈動壓力譜分布相似,但IDDES模型計算出的尖頻噪聲比實驗值偏高8 dB左右,而HRILES模型的結果與實驗數據非常接近.高頻寬帶噪聲則來源于縫翼后緣的渦脫落,其脫落頻率不僅取決于尾緣厚度還與邊界層位移厚度相關.P4點位于充分發展的湍流剪切層流動再附點附近,實驗結果分布表明該點為寬頻噪聲.兩種模型計算出的寬帶噪聲分布量級與實驗值符合較好,但是IDDES模型模擬出了顯著的窄帶峰值,HRILES模型基本保持了寬帶特性.在頻率高于10 kHz的部分,P1和P4點的計算結果存在幅值顯著衰減現象,這可能是由于空間網格密度不夠造成的,但由于遠場噪聲聲壓級統計要求只對256 Hz-10 kHz頻段進行寬帶濾波,因此對遠場噪聲計算結果影響不大.

圖11 各個站位的平均展向渦量分布Fig.11.Mean spanwise vorticity along the seven lines across.

圖12 各個站位的平均湍動能分布Fig.12.Mean turbulent kinetic energy along the seven lines across.

圖13 脈動壓力功率譜密度分布 (a) P1 ; (b)P4Fig.13.Frequency spectra of pressure fluctuations：(a) P1 ; (b) P4.

4.3 噪聲結果

圖14展示了HRILES模型計算得到的x-y截面瞬態壓力脈動云圖,采用來流動量q進行無量綱.圖中縫翼噪聲的強度遠大于其他區域,并有顯著的指向性特征,計算結果清楚地模擬出了組件之間存在的聲波反射、干涉現象.

圖14 瞬態脈動壓力云圖Fig.14.Contours of pressure fluctuation.

Pascioni和Cattafesta[39]測量了位于r=1 m,θ=287.5°處的噪聲譜,其中r以收起的縫翼前緣作為起點,θ 是以下游方向為起點逆時針為正.觀測點噪聲通過FW-H積分方法計算,聲源面選擇為壁面,聲源信息在非定常流動模擬時采樣,采樣頻率為10 kHz,共計8024個數據點.將采樣數據分為7個窗,數據重疊率為50%,采用Hamming窗對每個窗進行傅里葉變換,最終將每個窗的變換結果進行平均得到最終的聲壓級譜,計算聲壓級參考壓力值為 5×10-5Pa.由于實驗幾何模型展向長度為1 m,通過Kato公式對計算數據進行校正.從從圖15的結果看,觀測點獲得的遠場聲壓級譜也顯示出窄帶峰,這與 P1表面脈動壓力譜規律一致.計算結果在1.5 -10.0 kHz范圍內與實驗值相比符合較好,尤其是HRILES幅值和頻率都可以準確模擬,而IDDES的計算結果幅值偏大.

參考機體噪聲會議基準算例要求,圖16給出了位于 r=10C 處的聲壓級指向圖.不同于圖15實驗結果的展長1 m,r=10C 處噪聲計算要求展長統一為1 in,小于本文計算域展向長度(2 in),可以通過只統計z=0.5-1.5 in之間的壁面聲源數據實現,對頻段256 Hz 至10 kHz進行帶寬濾波得到最終聲壓級,聲壓參考值為 5×10-5Pa.IDDES 模型的結果表現出更強的偶極子聲源指向特性,第二第四象限的噪聲強度更大.Choudhari和Lockard[40]給出了 θ=270°處聲壓級的參考范圍為56.7-69.9 dB,平均值為61.5 dB.HRILES模擬得到的結果為61.5 dB,IDDES模擬得到的結果為66.6 dB,均在參考范圍之內.

圖17對比了分別采用縫翼、主翼、襟翼作為聲源面計算得到的遠場噪聲結果.不同組件計算結果表征的是各自壁面的非定常壓力脈動信息,其壓力脈動來源于流動非定常特性,但同時包括對其他組件產生的聲波的反射效應,可以作為研究聲源強度分布的參考.整體上縫翼噪聲在高升力體噪聲中占主導地位,其偶極子輻射指向特性也更顯著,其次為主翼噪聲,指向性不明顯.襟翼噪聲占比最小,也呈偶極子輻射特性,噪聲輻射在第一第三象限強度較大,聲源主要來自尾緣的渦脫落.

圖15 r=2.19C,θ=287.5° 觀測點聲壓級功率譜Fig.15.Power spectra density of sound pressure level at r=2.19C,θ=287.5°.

圖16 r=10C ,遠場聲壓級指向圖Fig.16.Directivity of SPL at r=10C.

圖17 各部件遠場(r=10C)聲壓級指向對比圖 (a) IDDES; (b) HRILESFig.17.Directivity of components' SPL at r=10C ：(a) IDDES; (b) HRILES.

5 結 論

本文發展了高精度混合RANS/HILES方法.對亞臨界態(Red=4.3×104)單圓柱繞流和30P30N高升力體基準算例進行了非定常模擬,結合聲比擬方法得到了氣動噪聲結果,并將HRILES 模型與DDES,IDDES模型進行了對比.

從單圓柱模擬的流動結果看,HRILES模型能夠模擬亞臨界態尾跡區層流轉捩現象,捕捉到背風面小尺度分離泡,因此對阻力系數的平均值和均方根值以及壁面壓力分布的預測結果比SST-DDES模型更接近實驗數據.遠場噪聲的功率譜密度分布的幅值與頻率都與實驗結果符合較好,證明了HRILES 模型具備預測復雜流動氣動噪聲問題的能力.

將HRILES模型應用于30P30N高升力體氣動噪聲問題,計算條件與機體噪聲研討會中規定的基準算例條件保持一致,并采用了JAXA提供的標準網格,同時采用IDDES模型模擬以對比模型差異的影響.從流動結果看,HRILES方法由于在分離區模型耗散更低,得到的縫翼凹腔中的渦結構更精細.對比剪切層各個站位的平均速度、展向渦量和湍動能分布,HRILES模型的結果與實驗數據更加接近,IDDES 模型的結果存在顯著的推遲轉捩的現象,導致剪切層上、中區域展向渦量偏大,湍動能偏低.二者計算出的壁面非定常壓力脈動的均方根值與實驗數據相符,并得到了正確的頻譜分布.從遠場噪聲結果看,噪聲頻譜分布特別是尖頻噪聲與實驗值符合較好,這說明數值模擬捕捉到了縫翼凹腔自激勵共振模態,但是IDDES模型的結果聲壓級幅值偏高.通過選取高升力體不同組件壁面作為聲源面,對比了各自對遠場噪聲分布的貢獻.結果都表明縫翼噪聲貢獻了大部分聲能,縫翼和襟翼噪聲都展現了顯著的偶極子特性,主翼的噪聲輻射指向特性不明顯.

綜上所述,本文通過算例應用與模型對比,確認了HRILES模型在典型氣動聲學問題中的適用性.HRILES作為廣義RANS/LES模型的一種,能夠在SST模型與HILES方法之間光滑過渡.對比HRILES和IDDES模型,二者在壁面區域都主要采用RANS模擬,且HRILES在大渦模擬的尺度判據計算上與IDDES 保持了一致,但后者在LES區域近似等效于Smagorinsky亞格子模型,這也是DES類方法的共性.HRILES方法具備了HILES方法對于含能尺度流動的高分辨率模擬能力,因此在混合層流動模擬結果優于傳統的IDDES方法.同時在邊界層采用SST模型降低隱式大渦模擬方法對壁面網格分布的要求,減輕計算資源需求,在中、高雷諾數流場模擬以及氣動降噪應用研究中具有優勢,未來將對模型在更高雷諾數的復雜外形應用中開展進一步研究.

感謝中山大學國家超級計算中心在計算資源方面提供的幫助.