面向混合網(wǎng)格高精度阻力預(yù)測的熵修正方法

張培紅，張耀冰，周桂宇，陳江濤，鄧有奇

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

Roe通量差分分裂格式[1]是以近似Riemann分解為基礎(chǔ)的Godunov類求解器，具有耗散小、接觸間斷分辨率高和激波捕獲性能強(qiáng)、可以較好模擬邊界層流動等優(yōu)點[2-7]，在亞跨超聲速流場模擬中得到廣泛應(yīng)用。但對于特定問題Roe格式有時不滿足熵條件，會產(chǎn)生非真實的膨脹波，導(dǎo)致非物理解，例如在鈍體繞流中發(fā)現(xiàn)的“Carbuncle”現(xiàn)象[8-9]以及在研究運(yùn)動激波經(jīng)過壓縮拐角問題時，發(fā)現(xiàn)的“扭曲馬赫桿”現(xiàn)象。Quirk指出，在強(qiáng)激波和聲速點附近，需要對原始的Roe格式進(jìn)行熵修正。Harten通過對Jacobian矩陣中接近0的特征值，引入額外的耗散來修正奇異的特征值，提出了一種熵修正公式[10]。van Leer等給出了Harten熵修正的數(shù)學(xué)推導(dǎo)和引入耗散的建議值[11]。隨后，CFD工作者針對不同的計算問題，發(fā)展了形式多樣的Roe格式熵修正方法，目前較為流行的有Muller型、Harten-Yee型和Harten-Hyman型3類[2-4]。

雖然熵修正對Roe格式非常重要，但國內(nèi)外專門針對熵修正研究的文章并不多見，缺乏系統(tǒng)的研究。2001年，Kermani和Plett[4]采用無黏Buergers方程和激波管兩個算例，對Harten-Hyman型和Hoffmann-Chiang型等不同熵修正方法的能力進(jìn)行了評估，指出這些熵修正方法沒有普適性，在聲速膨脹波附近無法完全消除非真實的膨脹波，并在Harten熵修正的基礎(chǔ)上，通過在熵修正區(qū)域增大熵修正的方法對熵修正進(jìn)行了改進(jìn)，從而在不影響其他計算區(qū)域的前提下，完全消除聲速膨脹波附近的非物理膨脹激波。2004年，Phongthanapanich和Dechaumphai[6]通過對存在數(shù)值激波穩(wěn)定性問題的研究，分析了van Leer、Sanders、Pandolfi和D’Ambrosio等提出的不同熵修正方法的能力和求解精度，并針對二維高速可壓縮流動存在強(qiáng)激波的問題，在van Leer等提出的Harten熵修正和Pandolfi、D’Ambrosio提出的多尺度耗散技術(shù)的基礎(chǔ)上，發(fā)展了一種更加適合于三角形非結(jié)構(gòu)網(wǎng)格的熵修正方法，通過定常、非定常高速可壓縮流動對改進(jìn)的熵修正方法進(jìn)行了評估。國內(nèi)，2009年，周禹和閻超[2]從理論分析入手, 通過激波管問題、前臺階流動和運(yùn)動激波的雙馬赫反射3個數(shù)值實驗，對針對Roe格式中較為流行的Muller型、Harten-Yee型和Harten-Hyman型等不同熵修正的性能做了深入研究，分析了不同熵修正方法的性能效果，并得出了Harten-Yee型熵修正方法對激波和非物理的膨脹激波均起修正作用，而Harten-Hyman型熵修正只對膨脹過程起修正作用，對激波情況無效的結(jié)論。

關(guān)于Roe格式熵修正方法研究的文章，要么是對已有熵修正的性能做評估，要么是通過改進(jìn)聲速膨脹波附近的熵修正以消除非物理膨脹激波，而有關(guān)通過改進(jìn)熵修正方法提高附面層內(nèi)阻力預(yù)測精度，尤其是提高三維非結(jié)構(gòu)網(wǎng)格阻力預(yù)測精度的文章還未見到。文獻(xiàn)[12]對傳統(tǒng)Green-Gauss梯度求解方法進(jìn)行了改進(jìn)，在前期研究的基礎(chǔ)上，本文基于非結(jié)構(gòu)混合網(wǎng)格Roe格式熵修正的特點，通過對傳統(tǒng)的Harten-Yee熵修正進(jìn)行改進(jìn)，既保留了使用熵修正帶來的程序魯棒性等優(yōu)點，同時把熵修正對阻力預(yù)測精度的影響降到最低，提出了一種可提高非結(jié)構(gòu)混合網(wǎng)格阻力預(yù)測精度的Roe格式熵修正方法。通過DLR-F4翼身組合體改進(jìn)方法前和改進(jìn)方法后的計算結(jié)果比較，驗證了方法的有效性。采用改進(jìn)后方法對第3屆AIAA阻力預(yù)測會議[13-20]的計算模型——DLR-F6 翼身組合體進(jìn)行了詳細(xì)的模擬，分析了網(wǎng)格收斂性和雷諾數(shù)的影響，進(jìn)一步驗證了改進(jìn)方法的可靠性和阻力預(yù)測能力。

1 數(shù)值方法

1.1 控制方程

本文采用中國空氣動力研究與發(fā)展中心(CARDC)自主開發(fā)的基于格心的并行非結(jié)構(gòu)混合網(wǎng)格雷諾平均Navier-Stokes(RANS)方程流場解算器MFlow。它可以使用任意形狀的網(wǎng)格單元，具有較大的靈活性。采用有限體積法對空間進(jìn)行離散，未知變量位于網(wǎng)格單元的體心。守恒形式的非定常可壓縮Navier-Stokes方程可寫為

(1)

式中：Ω為控制體體積；S為控制面面積；?Ω為控制體封閉面面積；W表示守恒變量；Fc表示無黏通量；Fv表示黏性通量。

本文主控方程對流項采用二階迎風(fēng)Roe通量差分分裂格式進(jìn)行離散，黏性項采用中心差分格式離散，并采用多重網(wǎng)格技術(shù)進(jìn)行收斂加速。湍流模型采用S-A(Spalart-Allmaras)一方程湍流模型，湍流控制方程空間離散采用一階迎風(fēng)格式。主控方程和湍流方程時間迭代均采用LU-SGS(Lower Upper-Symmetric Gauss-Seidel)方法。在進(jìn)行改進(jìn)熵修正方法前后對計算結(jié)果的影響分析時，梯度求解方法為節(jié)點型Green-Gauss法[12]，熵修正的參數(shù)δ*取0.3。

1.2 Roe通量差分分裂格式

Roe通量差分分裂格式以近似Riemann分解為基礎(chǔ)，具有高分辨定態(tài)激波的能力，能在1個網(wǎng)格之內(nèi)捕獲激波[5]。另外較其他一些通量分裂格式(如van Leer格式)，該格式的耗散較小，對于邊界層流動模擬是較好的選擇。Roe格式在控制體單元面上的通量表達(dá)式為

|ARoe|IJ(WJ-WI)]

(2)

式中：|ARoe|為Roe平均矩陣。

|ARoe|和左右狀態(tài)差的乘積的求法為

|ARoe|IJ(WJ-WI)=|ΔF1|+|ΔF2,3,4|+|ΔF5|

(3)

式中：

(4)

(5)

(6)

當(dāng)特征值很小時，Roe格式會違反熵條件，產(chǎn)生非物理解，如膨脹激波、Carbuncle現(xiàn)象等。為避免出現(xiàn)非物理解，需要對Roe平均矩陣的特征值進(jìn)行熵修正[21-25]。

2 Harten熵修正及其改進(jìn)

熵修正的方法很多，根據(jù)相關(guān)研究，對于存在激波的亞跨超聲速流場進(jìn)行數(shù)值模擬，比較好的是Harten-Yee熵修正方法[2],即

(7)

(8)

式中：δ*是一個小值，一般取0～0.4。

熵修正對Roe格式很重要，但是使用熵修正會增大耗散，特別是在附面層內(nèi)對物面法向速度的求解精度影響較大，從而導(dǎo)致有/無熵修正對阻力的求解精度影響很大[26-28]。

通常的Harten-Yee熵修正方法，由于翼面法向的速度是小量，如圖 1的A處所示，導(dǎo)致平行于流向的界面處出現(xiàn)很小的特征值，Harten-Yee

熵修正可能會顯著提高這一特征值，以避免出現(xiàn)非物理解。但是這樣就加大了該界面上物理量的耗散，導(dǎo)致不能準(zhǔn)確地再現(xiàn)附面層的速度剖面和壁面附近的切應(yīng)力，使摩擦阻力的預(yù)測出現(xiàn)偏差。針對Harten熵修正存在的這一缺陷，對Harten熵修正做出修正，思路是在基本不改變其他區(qū)域熵修正的同時，在平行于翼面流向的界面上不做熵修正，不改變此處的真實耗散，這樣既保留了使用熵修正帶來的程序魯棒性等優(yōu)點，同時把熵修正對阻力預(yù)測精度的影響降到了最低。

采用非結(jié)構(gòu)網(wǎng)格進(jìn)行黏性流場計算時，在物面附近附面層內(nèi)采用的是扁長型的三棱柱網(wǎng)格，見圖 2。對于每個網(wǎng)格面定義一個參數(shù)k，將每個面的值乘以式(8)中Harten-Yee熵修正的參數(shù)δ*，得到新的δ*′值為

δ*′=kδ*

(9)

參數(shù)k的值根據(jù)網(wǎng)格位置不同取為0～1，當(dāng)k=0時，表示此位置不使用熵修正，當(dāng)k=1時，表示此位置采用傳統(tǒng)的Harten-Yee熵修正，不影響熵修正的效果。希望在物面附近類似圖 1中A位置的網(wǎng)格區(qū)域不使用熵修正或者熵修正盡量小；但對物面附近類似圖 1中B位置的網(wǎng)格區(qū)域，雖然同樣有扁長型網(wǎng)格，因為這個位置可能會是駐點，或者激波所在區(qū)域，不希望在這些網(wǎng)格的任何界面上減小Harten-Yee熵修正，需要將這些區(qū)域剔除，因此使用界面上的法向速度與界面上速度的比值來剔除這些界面。為此，定義的參數(shù)k的表達(dá)式為

(10)

式中：d為所求面相鄰的兩個單元的體心之間的距離；s為該面的面積；Vn為界面的法向速度；V為總速度的標(biāo)量值。對于平行于物面流向的界面來說，d2?s，因此k≈0；對于垂直于翼面的界面，d2?s，因此k≈1；對于附面層外的網(wǎng)格，一般來說k≈1，不影響熵修正的效果；對于激波區(qū)域或駐點附近的網(wǎng)格，式(10)中第1項雖然很小，但是由于Vn≈V，第2項的值接近1，因此k≈1。

3 算例驗證

本文通過對DLR-F4翼身組合體跨聲速繞流流場進(jìn)行計算，以驗證改進(jìn)后的熵修正方法對計算結(jié)果的影響。 DLR-F4翼身組合體是一個典型的現(xiàn)代運(yùn)輸類型飛機(jī)外形，具有大展弦比機(jī)翼和大型客機(jī)類型的機(jī)身。該外形在歐洲的ONERA-S2MA、NLR-HST和DRA-8 ft×8 ft(1 ft=0.304 8 m)等風(fēng)洞中進(jìn)行了全面系統(tǒng)的試驗研究，有詳實可靠的風(fēng)洞試驗數(shù)據(jù)，是目前檢驗RANS流場解算器在復(fù)雜外形亞跨超聲速流動中模擬能力的一個很好的算例。在AIAA第1屆阻力會議上被作為標(biāo)模采用各種解算器進(jìn)行了計算[14]，計算與試驗取得了一致的結(jié)果。DLR-F4翼身組合體設(shè)計巡航馬赫數(shù)Ma為0.75，迎角α為0°，雷諾數(shù)Re為3.0×106，平均氣動弦長為0.141 2 m，機(jī)身長1.192 m，半展長為0.585 7 m，機(jī)翼參考面積為0.145 4 m2，力矩參考點為(0.504 9 m，0 m，0 m)。

圖 3是DLR-F4翼身組合體的計算網(wǎng)格示意圖，圖 3(a)為表面網(wǎng)格和對稱面網(wǎng)格，圖 3(b)為頭部表面網(wǎng)格和附面層內(nèi)三棱柱網(wǎng)格。計算時采用半模，遠(yuǎn)場邊界取約50倍的機(jī)翼平均氣動弦長，半模網(wǎng)格單元總數(shù)為2 164萬，其中三棱柱795萬, 四面體1 368萬。物面法向三棱柱層數(shù)為27層，第1層間距約為1.0×10-6m(y+≈1)，物面單元數(shù)為29.5萬。機(jī)翼后緣采用各向異性三角形網(wǎng)格，單元數(shù)為32個。

圖 4給出了無熵修正、Harten熵修正和改進(jìn)后的Harten熵修正3種方法對主控方程變量ρ的殘差收斂曲線影響。可以看出，無熵修正方法的ρ的殘差只下降了不到4個量級，而使用了熵修正方法后下降了7個量級，在6 000步后還將繼續(xù)下降，并且Harten熵修正和改進(jìn)后的Harten熵修正殘差收斂曲線基本一致，收斂性都很好。

表1列出了無熵修正(No modified)、Harten熵修正(Original)和改進(jìn)后的Harten熵修正(Modified)3種方法計算升力系數(shù)CL、阻力系數(shù)CD、俯仰力矩系數(shù)Cm的值，以及知名CFD軟件CFD++、CFL3D、NSU3D、USM3Dns和TAU、NLR、ONERA、DRA等風(fēng)洞得到的結(jié)果。為了更直觀表示出方法改進(jìn)前后氣動力和力矩的差別，圖 5給出了表1中結(jié)果的柱狀圖。可以看出，本文3種不同熵修正方法計算得到的結(jié)果與知名CFD軟件計算得到的結(jié)果基本一致，與風(fēng)洞試驗結(jié)果也吻合較好。

Table1Influenceofentropycorrectiononaerodynamiccoefficients

MethodCLCDCmOriginal0.544 00.031 15-0.162 1Modified0.540 30.030 79-0.161 1No0.540 30.030 79-0.161 1CFD++0.544 00.034 3-0.166 0CFL3D0.536 00.031 4-0.160 0NSU3D0.553 00.030 9-0.161 7USM3Dns0.539 00.029 3-0.157 1TAU0.548 00.030 6-0.165 9NLR0.481 20.027 9-0.131 2ONERA0.476 00.028 0-0.127 0DRA0.476 40.027 2-0.137 9

總的來說，兩種熵修正方法對氣動力系數(shù)和力矩系數(shù)影響較小，相較而言，改進(jìn)后的熵修正對氣動力和力矩影響更小，主控方程殘差收斂特性更好，下降量級更多，計算得到的氣動力系數(shù)和力矩系數(shù)與試驗結(jié)果更接近。改進(jìn)后的Harten熵修正方法計算結(jié)果與原始Harten熵修正方法計算結(jié)果相比，升力減小約0.7%，阻力減小3.6個阻力單位，低頭俯仰力矩減小約0.6%。改進(jìn)后的Harten熵修正方法計算結(jié)果與無熵修正計算結(jié)果相比，在取4位有效數(shù)字情況下氣動力和力矩系數(shù)值完全相同，升力系數(shù)只有在第5位有效數(shù)字，阻力系數(shù)在第6位有效數(shù)字才不相同，說明本文的修改將Harten熵修正對氣動力的影響降到最低限度，基本上可以忽略不計。

通過以上分析可以看出，改進(jìn)后的熵修正，軟件魯棒性和原始的Harten熵修正方法基本一樣，但計算精度提高，計算結(jié)果和無熵修正時的計算結(jié)果基本一致，達(dá)到了增加軟件魯棒性，同時提高軟件阻力預(yù)測準(zhǔn)度的目的。

圖 6是站位y/b=0.512處壓力系數(shù)Cp剖面的比較。從整體上講，熵修正對壓力剖面的影響很小，如圖 6(a)所示；圖 6(b)是激波位置的放大圖，原始熵修正使激波位置后移，而修改后的熵修正位置與無熵修正完全重合；圖 6(c)是上翼面后緣的局部放大圖，原始熵修正使這里的負(fù)壓增加。綜上所述，原始熵修正使激波后移，上翼面負(fù)壓增大，從而使得升力增大，阻力也增大。從這幾個放大圖還可以看出，本文改進(jìn)的熵修正壓力剖面與不使用熵修正的壓力剖面幾乎完全重合，因此對總體的氣動力影響很小。

第2節(jié)講到，改進(jìn)熵修正是希望在類似圖 1中A處的翼型中部區(qū)域附面層內(nèi)盡量減小或者不加熵修正，而在類似圖 1中B處的翼型前緣區(qū)域附面層內(nèi)熵修正盡量與原始的Harten熵修正保持一致。圖 7是機(jī)翼剖面弦向中段速度型的比較。熵修正對速度型影響比較小，但是從放大圖可見，原始的熵修正不僅影響速度矢量的大小，而且影響方向，而本文改進(jìn)的熵修正，速度型與沒使用熵修正的結(jié)果基本重合在一起。

圖 8是前緣附近的速度型，由圖可見，改進(jìn)后的熵修正與原始的熵修正基本重合，而與無熵修正的結(jié)果不同。因此，本文對熵修正的改進(jìn)達(dá)到了所希望的效果。

4 DLR-F6數(shù)值模擬結(jié)果及分析

DLR-F6翼身組合體[13,15]是第2、3屆阻力預(yù)測會議規(guī)定的計算外形，有詳細(xì)的風(fēng)洞試驗結(jié)果和各種知名CFD軟件的計算數(shù)據(jù)，被廣泛用于研究分析跨聲速流動現(xiàn)象，可以用來驗證各種CFD代碼的正確性和可靠性。DLR-F6翼身組合體構(gòu)型翼展為1.17 m，參考面積為0.145 4 m2，平均氣動弦長為0.141 2 m，設(shè)計巡航狀態(tài)為Ma=0.75,CL=0.5,Re=5×106(基于平均氣動弦長)。

為開展網(wǎng)格收斂性研究，采用第3屆阻力會議提供的粗、中、細(xì)三套網(wǎng)格進(jìn)行了計算，網(wǎng)格信息見表2，其中中等網(wǎng)格作為計算的基準(zhǔn)網(wǎng)格。網(wǎng)格單元包括三棱柱、金字塔和四面體，物面附近三棱柱的層數(shù)是變化的，最大為31層，在生成的三棱柱質(zhì)量不好的區(qū)域，例如翼身結(jié)合處，三棱柱的層數(shù)減少，用金字塔過渡到四面體，圖 9為DLR-F6翼身組合體構(gòu)型表面網(wǎng)格和對稱面網(wǎng)格。

表2 DLR-F6網(wǎng)格信息

表3給出了本文采用粗、中、細(xì)三套網(wǎng)格計算得到氣動力結(jié)果,其中:α為迎角，CL為升力系數(shù)，CD為阻力系數(shù)，CDp為壓差阻力系數(shù)，CDf為摩擦阻力系數(shù)，Cm為俯仰力矩系數(shù)。

根據(jù)文獻(xiàn)[16]，第3屆阻力會議DLR-F6的阻力系數(shù)統(tǒng)計結(jié)果的可信值域為[0.026 3，0.029 2]，其中壓差阻力系數(shù)為[0.013 8，0.016 8]，摩擦阻力系數(shù)為[0.011 3，0.012 8]，迎角為[-0.15°，0.36°]，俯仰力矩系數(shù)為[-0.158，-0.130]，本文中、細(xì)兩套網(wǎng)格的計算結(jié)果均落在其中。

圖10是阻力以及將阻力分解為壓阻和摩阻的網(wǎng)格收斂圖，橫坐標(biāo)取N-2/3，將阻力作為N-2/3的函數(shù)，N表示網(wǎng)格單元的數(shù)量，N-1/3表示網(wǎng)格的一維空間尺度，采用N-2/3是基于代碼的數(shù)值離散方法為二階精度，因此圖中若為直線就表示空間網(wǎng)格收斂為二階精度，(N-1/3)2為0時表示網(wǎng)格量趨近于無窮大。由圖可見，本文計算結(jié)果當(dāng)網(wǎng)格量趨于無窮時阻力收斂是單調(diào)的，而且收斂精度接近二階。壓阻隨網(wǎng)格加密減小，摩阻隨網(wǎng)格加密增大，由于壓阻減小得多，因而總阻力也隨網(wǎng)格加密而減小。

通過對以上3套網(wǎng)格的計算結(jié)果使用Richardson外插法[29]，可以估算出當(dāng)網(wǎng)格量趨于無窮時的阻力值。使用中等網(wǎng)格和細(xì)網(wǎng)格的插值結(jié)果不一定等于粗網(wǎng)格和中等網(wǎng)格的插值結(jié)果，如果插值結(jié)果相等，則圖 10中線條為直線。如圖 10(a)中箭頭所指，DLR-F6的插值結(jié)果為262.9 counts，與細(xì)網(wǎng)格阻力值差-9.6 counts。

表3 DLR-F6計算結(jié)果Table 3 Numerical results of DLR-F6

圖 11是本文預(yù)測結(jié)果與阻力會議中各知名CFD軟件預(yù)測結(jié)果以及試驗結(jié)果的對比圖，柱狀圖中給出的各軟件預(yù)測結(jié)果是進(jìn)行網(wǎng)格收斂性分析后插值得到的。從圖中可以看出，本文的預(yù)測結(jié)果與各知名CFD軟件的預(yù)測結(jié)果相當(dāng)，比試驗值略小，但處于合理范圍，進(jìn)一步驗證了本文方法的可靠性和正確性。

為了進(jìn)一步考核方法的可靠性，研究和分析雷諾數(shù)對翼身組合體氣動特性的影響，計算了Ma=0.75,Re=3×106, 5×106時DLR-F6的氣動特性和流場。圖 12給出了雷諾數(shù)對升力的影響曲線，可以看出，由于增加雷諾數(shù)減小了邊界層的厚度，從而使機(jī)翼的有效彎度增加，隨雷諾數(shù)增大，翼身組合體的升力系數(shù)增大，升力曲線斜率略增大。

圖 13給出了雷諾數(shù)對阻力的影響曲線，雷諾數(shù)的變化對壓差阻力和摩擦阻力均有影響，隨雷諾數(shù)增大，壓差阻力和摩擦阻力均減小，但相對來說，雷諾數(shù)對摩擦阻力的影響更明顯。圖 14給出了雷諾數(shù)對機(jī)翼剖面壓力分布的影響。隨雷諾數(shù)增大，機(jī)翼上表面前緣的壓力峰值減小，激波位置前移，y/b=0.150剖面的機(jī)翼上表面后緣壓力分布不同，主要是由于雷諾數(shù)變化引起分離氣泡位置和大小不同導(dǎo)致的。

圖 15給出了雷諾數(shù)對分離的影響。Re=3.0×106的分離氣泡相比Re=5.0×106要稍大一些，也說明由于增加雷諾數(shù)，使得機(jī)翼和機(jī)身的附面層變薄，從而延緩了流動分離。雷諾數(shù)的變化范圍不到兩倍，本文方法仍然計算出了雷諾數(shù)的影響，進(jìn)一步說明了本文計算方法的高精度和高可靠性。

5 結(jié) 論

1) 本文提出了一種可提高非結(jié)構(gòu)混合網(wǎng)格黏性計算精度的Harten-Yee熵修正改進(jìn)方法，提高了數(shù)值方法的魯棒性和阻力預(yù)測精度。

2) 通過對DLR-F6的數(shù)值模擬，詳細(xì)分析了網(wǎng)格收斂性和雷諾數(shù)的影響，進(jìn)一步考核和驗證了軟件的計算精度和可靠性。