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

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

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

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

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

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

1 數(shù)值方法

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

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

(1)

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

(2)

(3)

(4)

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

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

2 Green-Gauss方法及其改進(jìn)

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

(5)

(6)

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

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

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

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

(7)

式中：diJ為頂點(diǎn)i相鄰的第J個(gè)控制單元體心到頂點(diǎn)i的距離；UiJ為頂點(diǎn)i相鄰的第J個(gè)控制單元的體心值。

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

(8)

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

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

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

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

3 算例驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表2 CRM六面體網(wǎng)格預(yù)測(cè)結(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ù)測(cè)結(jié)果Table 3 Prediction results of CRM using hybrid grids

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

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

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

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

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

5 結(jié) 論

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

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

[1] 閻超, 席柯, 袁武, 等. DPW系列會(huì)議評(píng)述與思考[J]. 力學(xué)進(jìn)展, 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] 王運(yùn)濤, 孟德虹, 孫巖, 等. DLR-F6/FX2B翼身組合體構(gòu)型高精度數(shù)值模擬[J]. 航空學(xué)報(bào), 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é)報(bào), 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ù)測(cè)[J]. 航空學(xué)報(bào), 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] 王運(yùn)濤, 李松, 孟德虹, 等. 梯形翼高升力構(gòu)型的數(shù)值模擬技術(shù)[J]. 航空學(xué)報(bào), 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é)報(bào), 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] 張耀冰. 運(yùn)輸機(jī)氣動(dòng)特性混合網(wǎng)格數(shù)值模擬研究[D]. 綿陽(yáng)：中國(guó)空氣動(dòng)力研究與發(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ù)測(cè)精度分析[J]. 航空學(xué)報(bào), 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).