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

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

Roe通量差分分裂格式[1]是以近似Riemann分解為基礎的Godunov類求解器，具有耗散小、接觸間斷分辨率高和激波捕獲性能強、可以較好模擬邊界層流動等優(yōu)點[2-7]，在亞跨超聲速流場模擬中得到廣泛應用。但對于特定問題Roe格式有時不滿足熵條件，會產(chǎn)生非真實的膨脹波，導致非物理解，例如在鈍體繞流中發(fā)現(xiàn)的“Carbuncle”現(xiàn)象[8-9]以及在研究運動激波經(jīng)過壓縮拐角問題時，發(fā)現(xiàn)的“扭曲馬赫桿”現(xiàn)象。Quirk指出，在強激波和聲速點附近，需要對原始的Roe格式進行熵修正。Harten通過對Jacobian矩陣中接近0的特征值，引入額外的耗散來修正奇異的特征值，提出了一種熵修正公式[10]。van Leer等給出了Harten熵修正的數(shù)學推導和引入耗散的建議值[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型等不同熵修正方法的能力進行了評估，指出這些熵修正方法沒有普適性，在聲速膨脹波附近無法完全消除非真實的膨脹波，并在Harten熵修正的基礎上，通過在熵修正區(qū)域增大熵修正的方法對熵修正進行了改進，從而在不影響其他計算區(qū)域的前提下，完全消除聲速膨脹波附近的非物理膨脹激波。2004年，Phongthanapanich和Dechaumphai[6]通過對存在數(shù)值激波穩(wěn)定性問題的研究，分析了van Leer、Sanders、Pandolfi和D’Ambrosio等提出的不同熵修正方法的能力和求解精度，并針對二維高速可壓縮流動存在強激波的問題，在van Leer等提出的Harten熵修正和Pandolfi、D’Ambrosio提出的多尺度耗散技術的基礎上，發(fā)展了一種更加適合于三角形非結構網(wǎng)格的熵修正方法，通過定常、非定常高速可壓縮流動對改進的熵修正方法進行了評估。國內(nèi)，2009年，周禹和閻超[2]從理論分析入手, 通過激波管問題、前臺階流動和運動激波的雙馬赫反射3個數(shù)值實驗，對針對Roe格式中較為流行的Muller型、Harten-Yee型和Harten-Hyman型等不同熵修正的性能做了深入研究，分析了不同熵修正方法的性能效果，并得出了Harten-Yee型熵修正方法對激波和非物理的膨脹激波均起修正作用，而Harten-Hyman型熵修正只對膨脹過程起修正作用，對激波情況無效的結論。

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

1 數(shù)值方法

1.1 控制方程

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

(1)

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

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

1.2 Roe通量差分分裂格式

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

|ARoe|IJ(WJ-WI)]

(2)

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

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

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

(3)

式中：

(4)

(5)

(6)

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

2 Harten熵修正及其改進

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

(7)

(8)

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

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

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

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

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

δ*′=kδ*

(9)

參數(shù)k的值根據(jù)網(wǎng)格位置不同取為0～1，當k=0時，表示此位置不使用熵修正，當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的表達式為

(10)

式中：d為所求面相鄰的兩個單元的體心之間的距離；s為該面的面積；Vn為界面的法向速度；V為總速度的標量值。對于平行于物面流向的界面來說，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翼身組合體跨聲速繞流流場進行計算，以驗證改進后的熵修正方法對計算結果的影響。 DLR-F4翼身組合體是一個典型的現(xiàn)代運輸類型飛機外形，具有大展弦比機翼和大型客機類型的機身。該外形在歐洲的ONERA-S2MA、NLR-HST和DRA-8 ft×8 ft(1 ft=0.304 8 m)等風洞中進行了全面系統(tǒng)的試驗研究，有詳實可靠的風洞試驗數(shù)據(jù)，是目前檢驗RANS流場解算器在復雜外形亞跨超聲速流動中模擬能力的一個很好的算例。在AIAA第1屆阻力會議上被作為標模采用各種解算器進行了計算[14]，計算與試驗取得了一致的結果。DLR-F4翼身組合體設計巡航馬赫數(shù)Ma為0.75，迎角α為0°，雷諾數(shù)Re為3.0×106，平均氣動弦長為0.141 2 m，機身長1.192 m，半展長為0.585 7 m，機翼參考面積為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)格。計算時采用半模，遠場邊界取約50倍的機翼平均氣動弦長，半模網(wǎng)格單元總數(shù)為2 164萬，其中三棱柱795萬, 四面體1 368萬。物面法向三棱柱層數(shù)為27層，第1層間距約為1.0×10-6m(y+≈1)，物面單元數(shù)為29.5萬。機翼后緣采用各向異性三角形網(wǎng)格，單元數(shù)為32個。

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

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

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

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

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

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

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

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

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

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

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

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

根據(jù)文獻[16]，第3屆阻力會議DLR-F6的阻力系數(shù)統(tǒng)計結果的可信值域為[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]，本文中、細兩套網(wǎng)格的計算結果均落在其中。

圖10是阻力以及將阻力分解為壓阻和摩阻的網(wǎng)格收斂圖，橫坐標取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)格量趨近于無窮大。由圖可見，本文計算結果當網(wǎng)格量趨于無窮時阻力收斂是單調(diào)的，而且收斂精度接近二階。壓阻隨網(wǎng)格加密減小，摩阻隨網(wǎng)格加密增大，由于壓阻減小得多，因而總阻力也隨網(wǎng)格加密而減小。

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

表3 DLR-F6計算結果Table 3 Numerical results of DLR-F6

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

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

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

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

5 結 論

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

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