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

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

阻力的高精度預(yù)測對于現(xiàn)代飛行器的設(shè)計，特別是民用飛機和運輸類飛機的設(shè)計極其重要。一直以來，國外就十分重視CFD計算方法的阻力預(yù)測能力。為評估CFD阻力預(yù)測能力的發(fā)展現(xiàn)狀，AIAA應(yīng)用空氣動力學(xué)技術(shù)委員會舉辦了一系列阻力預(yù)測研討會[1-5]，截止目前，AIAA阻力預(yù)測研討會已經(jīng)舉辦了6屆，為世界范圍內(nèi)的大學(xué)、研究所和工業(yè)部門共同評估當前CFD方法預(yù)測運輸類飛行器跨聲速流場氣動力和力矩系數(shù)的發(fā)展水平，提供了一個公平的平臺，從而促進了CFD應(yīng)用水平的提高。

隨著計算機技術(shù)的不斷發(fā)展，CFD在湍流模型、網(wǎng)格技術(shù)、數(shù)值算法、可視化、并行計算等方面取得飛速發(fā)展，已成為型號設(shè)計部門的常規(guī)手段之一[6-9]。特別是非結(jié)構(gòu)網(wǎng)格生成技術(shù)的迅速發(fā)展，減少了網(wǎng)格生成人工工作量，大大提高了網(wǎng)格對復(fù)雜構(gòu)型的適應(yīng)能力和靈活性，推動了CFD解決復(fù)雜工程問題的能力[10-13]。但是，由于非結(jié)構(gòu)網(wǎng)格在黏性區(qū)域很難使用類似結(jié)構(gòu)網(wǎng)格的大長寬比網(wǎng)格，導(dǎo)致阻力預(yù)測精度降低，嚴重阻礙了非結(jié)構(gòu)網(wǎng)格的應(yīng)用和發(fā)展[14-16]。因此，如何解決非結(jié)構(gòu)網(wǎng)格阻力預(yù)測精度低與現(xiàn)代飛行器設(shè)計對阻力預(yù)測精度要求高的矛盾，是當前CFD領(lǐng)域亟需解決的關(guān)鍵技術(shù)和難題。

CFD計算過程中在網(wǎng)格單元面上重構(gòu)標量的值，以及計算擴散項和速度導(dǎo)數(shù)時都需要用到梯度，因此梯度的求解對于非結(jié)構(gòu)混合網(wǎng)格計算的精度和魯棒性起非常重要的作用[17-22]。通常的梯度求解方法有Least-Squares方法和Green-Gauss方法，由于Least-Squares方法對于曲面大壓縮比網(wǎng)格得到的結(jié)果很差，甚至?xí)葘嶋H值小一個量級，因此，Green-Gauss方法是目前應(yīng)用最廣的一種梯度求解方法[18]。Green-Gauss方法的關(guān)鍵是求解面心的值，其最大優(yōu)點是與計算通量的過程相似，不需要增加新的數(shù)據(jù)結(jié)構(gòu)，可以提高計算效率，但由于只使用了與當前單元共面的單元的信息，周圍其他單元的信息并沒有用到，在進行混合網(wǎng)格的梯度計算中存在缺陷，在不同類型單元交界處，得到的梯度變得極度不準確。為了克服Green-Gauss方法的缺點，提高計算準度，人們嘗試對Green-Gauss方法進行了多種改進[18-22]。一種改進是采用距離或者體積進行加權(quán)，這種方法雖然在一定程度上可以提高計算的準度，但當面心不在體心的連線上時，計算精度會很不理想[20]。另一種改進是先按照Green-Gauss方法求出梯度的初值，然后采用迭代法進行迭代，得到最終梯度的值，這種方法可以較好地改進梯度計算的準確性，但是由于迭代運算，可能會導(dǎo)致產(chǎn)生嚴重的振蕩[21]。Blazek進一步對Green-Gauss方法進行了改進[22]，他將與控制體有公共點的相鄰控制體的體心相連，形成新的控制單元，并將形成的控制單元作為新的控制體。這種改進方法的缺點是需要增加一個新的基于單元的數(shù)據(jù)結(jié)構(gòu)，用來提供體心和新控制體面之間的聯(lián)系，因此，這種方法不再是網(wǎng)格透明的(與單元信息無關(guān))，并且不能使用高效的循環(huán)[22]。

本文基于非結(jié)構(gòu)混合網(wǎng)格梯度求解的特點，在傳統(tǒng)的Green-Gauss梯度求解方法基礎(chǔ)上，通過增加計算所用的譜，采用更多周圍控制單元的信息，提出了一種適用于非結(jié)構(gòu)混合網(wǎng)格黏性計算的節(jié)點型Green-Gauss梯度求解方法。通過DLR-F4翼身組合體改進方法前后的計算結(jié)果比較，驗證了方法的有效性。采用改進后的梯度求解方法對第5屆AIAA阻力預(yù)測研討會的通用研究模型(Common Research Model，CRM)開展了詳細的計算研究，進一步驗證了改進方法的可靠性和阻力的預(yù)測能力。該方法克服了Green-Gauss方法和以往改進的Green-Gauss方法的缺點，提高了阻力的預(yù)測精度和程序的魯棒性。

1 數(shù)值方法

本文計算采用了課題組自主開發(fā)的MFlow解算器。MFlow是基于非結(jié)構(gòu)混合網(wǎng)格的用于求解亞跨超聲速流場的專業(yè)CFD軟件，軟件編寫采用模塊化設(shè)計思想。非結(jié)構(gòu)混合網(wǎng)格一般由四面體、六面體、三棱柱或者金字塔組成，本文在生成混合網(wǎng)格時，采用一種改進的推進層方法在物面黏性作用區(qū)生成大伸展比三棱柱形和金字塔形網(wǎng)格；在其他流動區(qū)域仍采用陣面推進方法生成四面體網(wǎng)格。這樣既保持了非結(jié)構(gòu)網(wǎng)格易于控制網(wǎng)格單元的大小、形狀及網(wǎng)格點的位置，具有靈活性大，對復(fù)雜外形的適應(yīng)能力強，并且可以方便地作自適應(yīng)計算和大規(guī)模分布式并行計算等優(yōu)點，又克服了非結(jié)構(gòu)網(wǎng)格進行黏性計算時計算效率和計算精度低的缺陷。

采用有限體積法對空間進行離散，未知變量位于網(wǎng)格單元的體心?？刂品匠虨槭睾阈问降姆嵌ǔ？蓧嚎sNavier-Stokes方程，即

(1)

式中：Ω為控制體的體積；?Ω為控制體封閉面的面積；W為守恒變量；Fc為無黏通量；Fv為黏性通量。定義分別為

(2)

(3)

(4)

式中：(x,y,z)為笛卡兒坐標系下的坐標分量；(u,v,w)為笛卡兒坐標系下的速度分量；(nx,ny,nz)為微元面的單位法向矢量在笛卡兒坐標系下的分量；τij為切應(yīng)力；ρ、p、E、H和V分別為密度、壓力、總能量、總焓和速度。

本文采用三重“V”循環(huán)多重網(wǎng)格方法進行收斂加速，主控方程對流項采用二階迎風Roe通量差分格式進行離散，黏性項采用中心差分格式離散。假定流場為完全湍流流場，湍流模型采用S-A一方程湍流模型[23]，湍流控制方程空間離散采用一階迎風格式。主控方程和湍流方程時間迭代均采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法。在使用原始Green-Gauss方法和節(jié)點型Green-Gauss方法對計算結(jié)果影響分析時，采用Roe格式，Harten-Yee熵修正。

2 Green-Gauss方法及其改進

對于本文采用的格心型網(wǎng)格，Green-Gauss方法的公式為

(5)

(6)

式中：NF為離散的控制體Ω的面數(shù)；SIJ為第J個面方向向外的面積矢量；UIJ為第J個面的面心值；UI為控制體單元I的體心值；UJ為控制體單元J的體心值。

Green-Gauss方法的關(guān)鍵是求解面心值，式(6)中使用面左右單元的體心值來求面心值，其最大優(yōu)點是與計算通量的過程相似，不需要增加新的數(shù)據(jù)結(jié)構(gòu)，可以提高計算效率，但由于只使用了與當前單元共面的單元信息，周圍其他單元的信息并沒有用到，如圖1所示，在進行混合網(wǎng)格的梯度計算中準度偏低。

為提高非結(jié)構(gòu)網(wǎng)格中梯度的求解準度，針對傳統(tǒng)Green-Gauss方法和Blazek改進方法求解梯度時存在的缺陷，本文提出了一種稱為“節(jié)點型Green-Gauss方法”的梯度求解方法。其基本思想是：通過增加計算所用的譜，采用更多周圍控制單元的信息來提高梯度求解的準度，計算過程中面心值用該面所有頂點的算術(shù)平均得到，而頂點值又利用它周圍所有單元的體心值通過距離加權(quán)平均得到，這樣就充分利用了周圍控制體的所有已知信息，如圖2所示。

假設(shè)控制體單元I,共有Ni個頂點，對于頂點i,有NI個相鄰的控制體單元。采用改進后的方法求解控制體單元I的梯度時，首先求解控制體單元的頂點的值，對于頂點i,其值UIi由其相鄰的控制單元的體心值通過距離加權(quán)得到，即

(7)

式中：diJ為頂點i相鄰的第J個控制單元體心到頂點i的距離；UiJ為頂點i相鄰的第J個控制單元的體心值。

控制體單元I的第J個面的面心值由式(7)中求得的頂點值算術(shù)平均得到,即

(8)

圖1 Green-Gauss方法Fig.1 Green-Gauss method

圖2 節(jié)點型Green-Gauss方法Fig.2 Node-based Green-Gauss method

式中：Nj為控制體單元I的第J個面的頂點個數(shù)。

最后把式(8)得到的面心值代入式(5)得到控制體單元I的梯度值。

3 算例驗證

為驗證改進后的梯度求解方法對計算結(jié)果的影響，本文對DLR-F4翼身組合體跨聲速繞流流場進行了計算。DLR-F4[24-25]是一個典型的現(xiàn)代運輸類飛機外形翼身組合體，具有大展弦比機翼和大型客機類型的機身，設(shè)計巡航馬赫數(shù)為0.75，機身長1.192m,半展長0.5857m,機翼參考面積為0.1454m2,平均氣動弦長為0.1412m,詳細參數(shù)見文獻[24]，在歐洲三大風洞ONERA-S2MA、NLR-HST和DRA-8ft×8ft (1ft=0.3048m)中進行了全面系統(tǒng)的試驗研究，是目前檢驗雷諾-平均Navier-Stokes(RANS)流場解算器在復(fù)雜外形亞跨超聲速流動中的模擬能力的一個很好的算例，作為算例，在AIAA第1屆阻力研討會上采用各種解算器進行了計算，并取得了一致的結(jié)果[2]。計算狀態(tài)：馬赫數(shù)Ma=0.75,迎角α=0°,雷諾數(shù)Re=3.0×106,其中雷諾數(shù)取基于機翼的平均氣動弦長0.1412m。

圖3為DLR-F4翼身組合體的計算網(wǎng)格示意圖，圖3(a)為表面網(wǎng)格和對稱面網(wǎng)格，圖3(b)為邊界層網(wǎng)格。計算時采用半模進行計算，半模網(wǎng)格單元總數(shù)為2164萬，其中四面體為1368萬，三棱柱為795萬。物面單元數(shù)為29.5萬，物面法向三棱柱層數(shù)為27層，第1層間距約為1.0×10-6m(y+≈1),前4層單元采用相同的間距，從第5層開始，變化率為1.2。機翼后緣采用各向異性三角形網(wǎng)格，單元數(shù)為32個。遠場邊界取約50倍的機翼平均氣動弦長。

圖3 DLR-F4翼身組合體網(wǎng)格Fig.3 Grid of DLR-F4 wing-body configuration

表1給出了梯度求解方法改進前后阻力計算結(jié)果與其他軟件和風洞試驗結(jié)果比較?？梢钥闯?，所有軟件計算結(jié)果都比試驗值偏大，改進梯度求解方法后，阻力減小了23.4個阻力單位，與試驗值更接近。

圖5為方法改進前后不同站位y/b剖面壓力系數(shù)Cp分布比較，c為弦長，由圖可見，改進后的節(jié)點型方法計算得到的前緣吸力峰比較光滑，并且峰值大，而原始Green-Gauss方法計算得到的結(jié)果有跳躍；兩種方法壓力剖面的區(qū)別主要在前緣吸力峰和激波的位置，以及波后的激波邊界層干擾區(qū)。由于這些位置都是流場變化劇烈的地方，梯度值較大，變化劇烈，對梯度的求解準度要求較高，準確求解這些位置的變量梯度非常重要。

圖4 梯度求解方法對收斂曲線的影響Fig.4 Influence of gradient calculation method on solution convergence

圖6為方法改進前后得到的翼根后緣區(qū)域分離氣泡的比較。很明顯，兩種方法計算得到的分離區(qū)域的大小和輪廓形狀基本相同，改進后的方法計算得到的翼上的大渦略偏前一點，細節(jié)上刻畫更明顯，中間區(qū)域出現(xiàn)了兩個小的分離氣泡，說明改進后的方法計算精度更高。

表1 不同梯度求解方法對阻力系數(shù)的影響

圖5 梯度求解方法對壓力分布的影響Fig.5 Influence of gradient calculation method on pressure distribution

通過以上收斂曲線、阻力系數(shù)、壓力分布和流場細節(jié)的對比分析可以看出，改進后的節(jié)點型方法可以更好地提高梯度的求解準度，更好地模擬梯度變化劇烈區(qū)域的壓力分布特征以及流場細節(jié)特征，使計算過程收斂更好，阻力預(yù)測精度更高。

圖6 梯度求解方法對分離氣泡的影響Fig.6 Influence of gradient calculation method on separation bubble

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

NASA的通用研究模型(CRM)是第5屆AIAA阻力預(yù)測研討會的計算外形，由機身和機翼組成，其中機身是典型的商業(yè)運輸機機身，機翼是跨聲速超臨界機翼, CRM表面網(wǎng)格見圖7。該外形在NASA Langley National Transonic Facility和NASA Ames 11-ft風洞都進行了詳盡的試驗，可以為各種CFD代碼提供驗證計算的標準數(shù)據(jù)。會議主辦方針對CRM提供了包括多塊結(jié)構(gòu)網(wǎng)格、重疊網(wǎng)格和非結(jié)構(gòu)網(wǎng)格在內(nèi)的標準網(wǎng)格系列，分為六面體網(wǎng)格和三棱柱/四面體混合網(wǎng)格兩種形式。六面體網(wǎng)格共有5套，網(wǎng)格單元數(shù)從638 976～40 894 464，分別用L1T、L2C、L3M、L4F、L5X表示。混合網(wǎng)格是在六面體網(wǎng)格的基礎(chǔ)上生成的，在壁面附近將一個六面體剖分成兩個三棱柱，在其他區(qū)域?qū)⒁粋€六面體剖分成6個四面體，編號與剖分前的六面體網(wǎng)格相同，分別為L1T、L2C、L3M、L4F。

第5屆AIAA阻力預(yù)測研討會中，來自9個國家的22家單位針對該模型完成了阻力預(yù)測并提交了計算結(jié)果，對固定升力巡航狀態(tài)的阻力預(yù)測能力進行了檢驗。計算狀態(tài)：Ma=0.85,CL=0.500 0,Re=5.0×106。

首先對CRM進行了網(wǎng)格收斂性研究，主要目標是估算氣動力的網(wǎng)格收斂解，也就是估計當網(wǎng)格量趨近于無窮大時氣動力的值。表2和表3分別給出了不同六面體網(wǎng)格和混合網(wǎng)格計算得到的，在設(shè)計點CL=0.5時的迎角、升力、阻力等，并給出了NASA Langley National Transonic Facility的試驗結(jié)果。從表3可以看出，計算得到的定升力迎角與試驗值均有一定的差異。在前幾屆的阻力預(yù)測會議中也發(fā)現(xiàn)了這種情況。根據(jù)Shahyar和Neal[25]的分析， 差異產(chǎn)生的原因包括風洞洞壁干擾和阻塞影響、支架影響、氣動彈性效應(yīng)，機翼后緣鈍度以及固定轉(zhuǎn)捩的影響，因此這里比較定升力時的阻力更有意義。

圖7 CRM表面網(wǎng)格(L1T)Fig.7 Surface grid of CRM (L1T)

表2 CRM六面體網(wǎng)格預(yù)測結(jié)果Table 2 Prediction results of CRM using hex grids

參數(shù)L1TL2CL3ML4FL5Xα/(°)2.2632.2102.1892.1742.165CL0.50000.50060.50030.50010.4996CD/counts263.0254.6252.2250.8250.1CDp/counts152.3141.4138.0135.8134.9CDf/counts110.7113.3114.3114.9115.1

表3 CRM混合網(wǎng)格預(yù)測結(jié)果Table 3 Prediction results of CRM using hybrid grids

從表2可以看出，使用六面體網(wǎng)格時，L4F和L5X得到的阻力跟試驗預(yù)測只有不到2 counts(1 counts=阻力系數(shù)0.000 1)的差別，達到了很高的阻力預(yù)測精度。由于混合網(wǎng)格是由同樣編號的六面體網(wǎng)格經(jīng)過剖分得到的，因此比較同樣編號的六面體網(wǎng)格和混合網(wǎng)格，可以得到不同網(wǎng)格形式對阻力預(yù)測的影響。在本文的研究中發(fā)現(xiàn)，同樣編號的混合網(wǎng)格預(yù)測的摩擦阻力系數(shù)CDf和六面體網(wǎng)格基本沒有差別，但是預(yù)測的壓差阻力系數(shù)CDp要比六面體網(wǎng)格稍大。

圖8為采用六面體網(wǎng)格計算得到的阻力、壓差阻力和摩擦阻力的網(wǎng)格收斂圖，橫坐標取N-2/3,N表示網(wǎng)格單元數(shù)，采用N-2/3是基于數(shù)值方法為二階精度，對于由粗到細的一套自相似網(wǎng)格，表示其網(wǎng)格尺寸的二次方，因此直線就表示空間網(wǎng)格收斂為二階精度，N-2/3為0時表示網(wǎng)格量為無窮大。由圖可見，當網(wǎng)格量趨于無窮時阻力收斂是單調(diào)的，而且收斂精度接近二階。壓阻隨網(wǎng)格加密減小，而摩阻隨網(wǎng)格加密增大，由于壓阻減小得多，因而總阻力也隨網(wǎng)格加密而減小。

圖9給出了第5屆阻力會議組委會統(tǒng)計的阻力預(yù)測值柱狀圖[26]，包括OVERFLOW、FLUENT、CFD++、EDGE、FUN3D、HIFUN、NSU3D、TAU等知名CFD軟件的計算結(jié)果，以及NTF風洞和Ames風洞的試驗結(jié)果。柱狀圖中給出的各軟件的計算結(jié)果是采用Richardson外插方法[27]進行網(wǎng)格收斂性分析后，插值得到的。第4列和第5列分別是MFlow軟件使用六面體和混合網(wǎng)格預(yù)測的阻力，從圖可以發(fā)現(xiàn)，MFlow的結(jié)果位于兩個風洞試驗結(jié)果的區(qū)間內(nèi)，而且與計算的平均值很接近，特別是六面體網(wǎng)格結(jié)果。這也證明了程序在翼身組合體外形阻力預(yù)測方面有著可靠的計算精度。

圖8 阻力系數(shù)的網(wǎng)格收斂性Fig.8 Grid convergence characteristics of drag coefficients

圖9 AIAA第5屆阻力預(yù)測研討會CRM阻力預(yù)測結(jié)果比較[26]Fig.9 Results comparison of drag prediction with CRM in the 5th AIAA Drag Prediction Workshop [26]

5 結(jié) 論

1) 通過改進傳統(tǒng)Green-Gauss梯度求解方法，發(fā)展了節(jié)點型Green-Gauss方法，克服了現(xiàn)有梯度求解方法的缺陷，更好地模擬了梯度變化劇烈區(qū)域的壓力分布特征，提高了數(shù)值方法的魯棒性和阻力預(yù)測精度。

2) 通過對NASA通用研究模型的數(shù)值模擬，以及網(wǎng)格收斂性和數(shù)據(jù)對比分析，進一步考核和驗證了軟件的計算精度和可靠性。

[1] 閻超, 席柯, 袁武, 等. DPW系列會議評述與思考[J]. 力學(xué)進展, 2011, 41(6): 776-784.

YAN C, XI K, YUAN W, et al. Review of the drag prediction workshop series[J]. Advances in Mechanics, 2011, 41(6): 776-784 (in Chinese).

[2] LEVY D W, ZICKUHR T, VASSBERG J C, et al. Summary of data from the first AIAA CFD drag prediction workshop: AIAA-2002-0841[R]. Reston，VA: AIAA, 2002.

[3] LAFLIN K R, KLAUSMEYER S M, ZICKUHR T, et al. Summary of the second AIAA CFD drag prediction workshop[J]. Journal of Aircraft, 2005, 42(5): 1165-1178.

[4] VASSBERG J C, TINOCO E N, MANI M, et al. Summary of the third AIAA CFD drag prediction workshop[J]. Journal of Aircraft, 2008, 45(3): 781-798.

[5] VASSBERG J C, TINOCO E N, MANI M, et al. Summary of the fourth AIAA CFD drag prediction workshop: AIAA-2010-4547[R]. Reston，VA: AIAA, 2010.

[6] SLOTNICK J, KHODADOUST A, ALONSO J, et al. CFD vision 2030 study: A path to revolutionary computational aerosciences: NASA/CR-2014-218178[R]. Washington, D. C.: NASA， 2014.

[7] BARTH T J. An overview of combined uncertainty and a posteriori error bound estimates for CFD calculations: AIAA-2016-1062[R]. Reston，VA: AIAA, 2016.

[8] 王運濤, 孟德虹, 孫巖, 等. DLR-F6/FX2B翼身組合體構(gòu)型高精度數(shù)值模擬[J]. 航空學(xué)報, 2016, 37(2): 484-490.

WANG Y T, MENG D H, SUN Y, et al. High-order accuracy numerical simulation of DLR-F6/FX2B wing-body configuration[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(2): 484-490 (in Chinese).

[9] JOHNSON F T, TINOCO E N, YU N J. Thirty years of development and application of CFD at Boeing Commercial Airplanes, Seattle[J]. Computer and Fluid, 2005, 34(10): 1115-1151.

[10] 張健, 鄧有奇, 李彬, 等. 一種適用于三維混合網(wǎng)格的GMRES加速收斂新方法[J]. 航空學(xué)報, 2016, 37(11): 3226-3235.

ZHANG J, DENG Y Q, LI B, et al. A new method to accelerate GMRES’s convergence applying to three-dimensional hybrid grid[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3226-3235 (in Chinese).

[11] 徐嘉, 劉秋洪, 蔡晉生, 等. 基于隱式嵌套重疊網(wǎng)格技術(shù)的阻力預(yù)測[J]. 航空學(xué)報, 2013, 34(2): 208-217.

XU J, LIU Q H, CAI J S, et al. Drag prediction based on overset grids with implicit hole cutting technique[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(2): 208-217 (in Chinese).

[12] 王運濤, 李松, 孟德虹, 等. 梯形翼高升力構(gòu)型的數(shù)值模擬技術(shù)[J]. 航空學(xué)報, 2014, 35(12): 3213-3221.

WANG Y T, LI S, MENG D H, et al. Numerical simulation technology of high lift trapezoidal wing configuration[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(12): 3213-3221 (in Chinese).

[13] 康忠良, 閻超. 適用于混合網(wǎng)格的約束最小二乘重構(gòu)方法[J]. 航空學(xué)報, 2012, 33(9): 1598-1605.

KANG Z L, YAN C. Constrained least-squares reconstruction method for mixed grids[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(9): 1598-1605 (in Chinese).

[14] LEVY D W, CHAN M S. Comparison of viscous grid layer growth rate of unstructured grids on CFD drag prediction workshop results: AIAA-2010-4671[R]. Reston，VA: AIAA, 2010.

[15] MAVRIPLIS D J. Viscous flow analysis using a parallel unstructured multigrid solver[J]. AIAA Journal, 2000, 38(11): 2067-2076.

[16] ALLMARAS S R, BUSSOLETTI J E. Algorithm issues and challenges associated with the development of robust CFD codes[R]. New York: Springer, 2009: 1-19.

[17] TINOCO E N, BOGUE D R, KAO T J, et al. Progress toward CFD for full flight envelope[J]. The Aeronautical Journal, 2005, 109(1100): 451-460.

[18] EMRE S, CHRISTOPH B. Gradient calculation methods on arbitrary polyhedral unstructured meshes for cell-centered CFD solvers: AIAA-2014-1440[R]. Reston,VA: AIAA, 2014.

[19] SHIMA E, KITAMURA K, HAGA T. Green-Gauss/weighted-least-squares hybrid gradient reconstruction for arbitrary polyhedral unstructured grids[J]. AIAA Journal, 2013, 51(11): 2740-2747.

[20] ANDREW W C, ANDREW J D. Towards accurate flow predictions using unstructured meshes:AIAA-2009-3650[R]. Reston, VA: AIAA, 2009.

[21] 張耀冰. 運輸機氣動特性混合網(wǎng)格數(shù)值模擬研究[D]. 綿陽：中國空氣動力研究與發(fā)展中心， 2014.

ZHANG Y B. Numerical simulation of transport aircraft’s aerodynamic characteristics using mixed grid[D]. Mianyang：China Aerodynamics Research and Development Center， 2014 (in Chinese).

[22] BLAZEK J. Computational fluid dynamics: Principles and applications[M]. Oxford: Elsevier Science Ltd.， 2001：129-179.

[23] SAXENA S K, MANOJ T N. Implementation and testing of Spalart Allmaras model in a multi-block code: AIAA-2002-0835[R]. Reston,VA: AIAA, 2002.

[24] REDEKER G. DLR-F4 wing body configuration, in Chapter B of a selection of experimental test cases for the validation of CFD codes: AGARD-1994-AR-303[R]. Paris: AGARD, 1994.

[25] SHAHYAR Z P, NEAL T F. Assessment of the unstructured grid software TetrUSS for drag prediction of the DLR-F4 configuration: AIAA-2002-0839[R]. Reston,VA: AIAA, 2002.

[26] TINOCO E, LEVY D. DWP 5 summary of participant da-ta[C]∥5th CFD Drag Prediction Workshop. Reston, VA: AIAA, 2012.

[27] 李釗, 陳海昕, 張宇飛. 基于廣義Richardson外插方法的阻力預(yù)測精度分析[J]. 航空學(xué)報, 2015, 36(7): 2105-2114.

LI Z, CHEN H X, ZHANG Y F. Accuracy analysis of drag prediction based on generalized Richardson extrapolation[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(7): 2105-2114 (in Chinese).