王運(yùn)濤

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

由AIAA發(fā)起并主導(dǎo)的DPW(Drag Prediction Workshop)系列會議成功舉辦了6屆，吸引了世界范圍內(nèi)相關(guān)研究機(jī)構(gòu)的廣泛參與，已經(jīng)成為CFD驗(yàn)證和確認(rèn)領(lǐng)域最具代表性的國際會議之一。DPW系列會議的主要目的是評估現(xiàn)代CFD技術(shù)模擬運(yùn)輸機(jī)高速構(gòu)型阻力預(yù)測能力，提供一個公開、公正的研究平臺，為CFD技術(shù)下一步的發(fā)展提供意見和建議。

2001年6月召開的第1屆DPW會議(DPW I)選擇了DLR-F4翼身組合體構(gòu)型作為基準(zhǔn)模型，研究工況為固定升力系數(shù)的阻力預(yù)測(Case 1)與極曲線預(yù)測(Case 2)，采用14種軟件獲得了18組Case 1的計算結(jié)果[1]。

2003年6月召開的第2屆DPW會議(DPW Ⅱ)選擇了DLR-F6翼身組合體構(gòu)型和DLR-F6翼/身/架/艙兩種高速構(gòu)型作為基準(zhǔn)模型，必選研究工況包括固定升力系數(shù)下的網(wǎng)格收斂性研究(Case 1)和極曲線預(yù)測(Case 2)，可選工況包括固定升力系數(shù)下的轉(zhuǎn)捩位置影響(Case 3)及固定升力系數(shù)下的阻力發(fā)散特性(Case 4)。這次會議的另一個研究目的是評估CFD技術(shù)模擬掛架/短艙安裝阻力的能力，采用22種軟件獲得了29組Case 1的計算結(jié)果[2]。

由于DPW Ⅱ的數(shù)據(jù)提供者普遍碰到了翼身結(jié)合部流動局部分離導(dǎo)致的氣動特性收斂困難，2006年6月召開的第3屆DPW會議(DPW Ⅲ)選擇了DLR-F6和DLR-F6_FX2B兩種翼身組合體構(gòu)型作為基準(zhǔn)模型。其中DLR-F6_FX2B翼身組合體構(gòu)型對DLR-F6-WB構(gòu)型翼身結(jié)合部后緣進(jìn)行了修型處理,并將數(shù)值模擬的雷諾數(shù)提高到500萬，以避免翼身結(jié)合部的局部分離。與DPW Ⅰ和DPW Ⅱ事先公布試驗(yàn)結(jié)果不同，DPW Ⅲ采用了“blind test”方式開展研究工作，即數(shù)值模擬提供者事前并不知道試驗(yàn)結(jié)果。另外，為了鼓勵更多的研究人員參與，基準(zhǔn)模型還包括兩個單獨(dú)的機(jī)翼構(gòu)型(DPW-W1、DPW-W2)。研究工況包括固定升力系數(shù)下的網(wǎng)格收斂性研究(Case 1A、Case 2A)和極曲線預(yù)測(Case 1B、Case 2B)。此次會議另一目的是評估CFD技術(shù)模擬局部修型引起氣動特性差量的能力，采用18種軟件獲得了26組Case 1A的計算結(jié)果[3]。閻超等[4]對DPW Ⅰ～DPW Ⅳ會議的情況進(jìn)行了概述。

基于前3屆DPW會議的經(jīng)驗(yàn)和教訓(xùn)，以及CFD驗(yàn)證與確認(rèn)工作的迫切需求，從第4屆DPW會議(DPW IV)到第6屆DPW會議(DPW VI)，基準(zhǔn)模型采用了波音公司重新設(shè)計的CRM(Common Research Model)[5]，并在美國、歐洲的風(fēng)洞中開展了多項(xiàng)試驗(yàn)，積累了更加詳細(xì)、豐富的試驗(yàn)數(shù)據(jù)。DPW Ⅳ～DPW Ⅵ的研究內(nèi)容主要包括網(wǎng)格收斂性、掛架/短艙安裝阻力預(yù)測、機(jī)翼下洗流動影響、馬赫數(shù)/迎角影響、雷諾數(shù)影響、靜氣動彈性影響等。本文介紹了CRM及開展的風(fēng)洞試驗(yàn)，概述了DPW Ⅳ～Ⅵ會議的情況，從計算網(wǎng)格生成、計算方法與湍流模型、計算結(jié)果與試驗(yàn)結(jié)果的對比、優(yōu)化設(shè)計等方面總結(jié)CFD驗(yàn)證與確認(rèn)工作的進(jìn)展，并給出了進(jìn)一步開展CFD驗(yàn)證與確認(rèn)工作的思考和建議。

1 CRM及風(fēng)洞試驗(yàn)

CRM由波音公司主導(dǎo)設(shè)計，主要目的是為CFD的驗(yàn)證和確認(rèn)工作提供基本的、可公開的研究構(gòu)型和豐富完備的風(fēng)洞試驗(yàn)數(shù)據(jù)。該模型的設(shè)計綜合考慮了CFD驗(yàn)證和確認(rèn)研究工作的需求，工業(yè)部門、政府部門的建議，以及當(dāng)前的氣動設(shè)計和風(fēng)洞試驗(yàn)?zāi)芰ΑＴ撃Ｐ瓦x擇了機(jī)身/機(jī)翼/發(fā)房/平尾構(gòu)型(見圖1)，采用下單翼/翼吊雙發(fā)布局，是典型的寬體民機(jī)布局。CRM包括了多種構(gòu)型，主要有翼身組合體(CRM-WB)、翼/身/掛架/短艙構(gòu)型(CRM-WBPN)、翼身組合體加0°平尾(CRM-WBT0)、加2°平尾(CRM-WBT+2)、加-2°平尾(CRM-WBT-2)這5種構(gòu)型。該模型的設(shè)計馬赫數(shù)為0.85，設(shè)計升力系數(shù)為0.50。該模型的參考面積Sref=383.690 m2，平均氣動弦長c=7.005 32 m，展長b=58.762 9 m、梢根比λ=0.275，展弦比AR=9.0，1/4弦線后掠角Λc/4= 35.0°。

CRM的風(fēng)洞試驗(yàn)是在NASA FA/SFW項(xiàng)目的資助下完成的[6-7]，2010—2013年在NASA Langley的NTF(National Transonic Facility)風(fēng)洞完成了兩期試驗(yàn)(見圖2)，試驗(yàn)構(gòu)型包括CRM-WB、CRM-WBPN、CRM-WBT0、CRM-WBT+2、CRM-WBT-2。馬赫數(shù)Ma=0.7～0.87，雷諾數(shù)Re=5.0×106～30×106，迎角α=-3°～12°(Re=5.0×106)、α=-3°～6°(Re=19.8×106,30.0×106)，溫度T=-250～120 °F(1 °F=17.22 ℃) ， 試驗(yàn)結(jié)果包括氣動特性、表面壓力分布及模型變形測量數(shù)據(jù)等。

圖1 CRM構(gòu)型Fig.1 Common research model configuration

圖2 CRM在NTF風(fēng)洞中的安裝照片F(xiàn)ig.2 Photo of CRM in the National Transonic Facility (NTF) wind tunnel

2010年4月，在NASA Ames的11 ft (1 ft=304.8 mm) TWT(Transonic Wind Tunnel)風(fēng)洞中完成了一期對比試驗(yàn)，試驗(yàn)構(gòu)型同樣是5種，試驗(yàn)馬赫數(shù)Ma=0.7～0.87，雷諾數(shù)Re=5.0×106，迎角α=-3°～12°,溫度T=100 °F。試驗(yàn)結(jié)果除了氣動特性、表面壓力分布及模型變形測量數(shù)據(jù)外，還包括壓敏漆(Pressure Sensitive Paint，PSP)試驗(yàn)、摩阻測量和粒子圖像測速(Particle Image Velocimetry，PIV)試驗(yàn)。2014年2月，在歐盟項(xiàng)目的資助下(FP/2007-2013)，該模型在德國哥廷根的ETW(European Transonic Wind tunnel)又開展了一期對比試驗(yàn)[8]，試驗(yàn)構(gòu)型僅限于CRM-WBT0，試驗(yàn)馬赫數(shù)、雷諾數(shù)、迎角范圍與NTF相同，溫度范圍T=-249～83.93 °F，試驗(yàn)結(jié)果包括氣動特性、表面壓力分布及模型變形測量等。此外，2012年，日本宇航院(JAXA)也在JAXA 2 m×2 m跨聲速風(fēng)洞中完成了80%縮比的CRM風(fēng)洞試驗(yàn)[9]；近期，法國國家航天航空研究中心(ONERA)也在S1MA風(fēng)洞中完成了CRM的風(fēng)洞試驗(yàn)[10]。

NTF、TWT和ETW的風(fēng)洞試驗(yàn)均采用了相同的縮比(0.027)模型，唯有尾部支撐方式和風(fēng)洞試驗(yàn)數(shù)據(jù)修正方法略有不同，3座風(fēng)洞的試驗(yàn)結(jié)果均已經(jīng)在互聯(lián)網(wǎng)發(fā)布。特別需要指出的是，3座風(fēng)洞的試驗(yàn)結(jié)果均沒有進(jìn)行支撐干擾修正，這會對數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比產(chǎn)生重要影響。NTF與TWT的風(fēng)洞試驗(yàn)結(jié)果相比，相同試驗(yàn)條件下，升力系數(shù)與力矩系數(shù)的差別很小，阻力系數(shù)的差量在10 counts以內(nèi)(1 count=0.000 1)。ETW與NTF的風(fēng)洞試驗(yàn)結(jié)果相比，在Ma=0.85的條件下，NTF試驗(yàn)獲得的阻力系數(shù)比ETW試驗(yàn)獲得的阻力系數(shù)低16 counts；在Re=5.0×106的條件下，NTF試驗(yàn)獲得的低頭力矩系數(shù)小于ETW試驗(yàn)獲得的低頭力矩系數(shù)；兩個風(fēng)洞之間的測壓試驗(yàn)結(jié)果吻合良好；在數(shù)據(jù)誤差范圍內(nèi)，兩個風(fēng)洞之間的變形測量試驗(yàn)結(jié)果吻合。由以上對比可以得到如下結(jié)論：在Ma=0.85、Re=5.0×106的條件下，不同風(fēng)洞之間的升力系數(shù)、壓力分布、變形測量結(jié)果吻合良好，力矩系數(shù)有差異，阻力系數(shù)相差10 counts以上。

2 DPW Ⅳ～Ⅵ會議概述

DPW Ⅳ會議于2009年6月在美國德克薩斯州(Texas)的圣安東尼奧市(San Antonio)召開[11-12]。此次會議選擇了CRM-WBT構(gòu)型作為基準(zhǔn)研究模型，必選工況包括固定升力系數(shù)下的網(wǎng)格收斂性研究(Case 1A)、下洗影響研究(Case 1B)，可選工況包括阻力發(fā)散馬赫數(shù)模擬(Case 2)、雷諾數(shù)影響(Case 3)，采用“blind test”方式開展。來自世界各地的17個研究機(jī)構(gòu)采用16款軟件提供了28組Case 1A的數(shù)值模擬結(jié)果。其中，有17組結(jié)果采用了非結(jié)構(gòu)網(wǎng)格技術(shù)，11組結(jié)果采用了結(jié)構(gòu)網(wǎng)格技術(shù)；18組結(jié)果采用了Spalart-Allmaras (S-A) 一方程湍流模型及其修正形式，7組結(jié)果采用了Menter 剪切應(yīng)力輸運(yùn)(SST)兩方程湍流模型及其修正形式。主要結(jié)論如下：采用非結(jié)構(gòu)網(wǎng)格技術(shù)得到的數(shù)值模擬結(jié)果的散布度要大于采用結(jié)構(gòu)網(wǎng)格技術(shù)獲得的數(shù)值模擬結(jié)果；翼身結(jié)合部后緣的分離泡尺寸散布度較大；在機(jī)翼展向η=0.845 6站位，激波位置散布度較大，但隨著網(wǎng)格規(guī)模的增加散布度減少；大多數(shù)結(jié)果模擬CRM-WBT0、CRM-WBT+2兩種構(gòu)型的阻力與力矩特性相似，但模擬CRM-WBT-2構(gòu)型的阻力與力矩特性差別較大；在迎角為4°時，各組結(jié)果之間力矩特性差異明顯；與阻力系數(shù)和力矩系數(shù)的絕對量相比較，配平得到的阻力和水平尾翼偏角差別較?。桓鹘M結(jié)果之間阻力發(fā)散馬赫數(shù)預(yù)測能力與雷諾數(shù)影響預(yù)測能力相當(dāng)。

DPW Ⅴ會議于2012年6月在美國路易斯安那州(Louisiana)的新奧爾良市(New Orleans)召開[13-14]。此次會議選擇了CRM-WB構(gòu)型作為基準(zhǔn)研究模型，必選工況包括固定升力系數(shù)下的網(wǎng)格收斂性研究(Case 1)、固定馬赫數(shù)下的抖振特性研究(Case 2)，可選工況包括湍流模型的確認(rèn)(Case 3)。來自世界各地的22個研究機(jī)構(gòu)采用17款軟件提供了54組Case 1的數(shù)值模擬結(jié)果。其中，26組結(jié)果采用了非結(jié)構(gòu)網(wǎng)格技術(shù)，12組結(jié)果采用了結(jié)構(gòu)網(wǎng)格技術(shù)，16組結(jié)果采用了自行生成的網(wǎng)格；34組結(jié)果采用了Spalart-Allmaras一方程湍流模型及其修正形式，12組結(jié)果采用了Menter SST兩方程湍流模型及其修正形式。主要結(jié)論如下：對于Case 1，數(shù)據(jù)散布度要小于DPW IV的結(jié)果，但依然較大；阻力系數(shù)與試驗(yàn)結(jié)果的對比情況要好于迎角與俯仰力矩系數(shù)與試驗(yàn)的對比情況。對于Case 2，采用RST模型與k-ε-RT模型的結(jié)果不理想；迎角為2.5°時，計算結(jié)果之間的吻合度較好；迎角為4.0°時，計算結(jié)果之間的分散度較大；計算得到的俯仰力矩系數(shù)與試驗(yàn)結(jié)果相差較大；靜氣動彈性變形影響顯著；大迎角狀態(tài)下，流動由顯著的激波誘導(dǎo)分離所主導(dǎo)。

DPW Ⅵ會議于2016年6月在美國華盛頓哥倫比亞特區(qū)(Washington, D.C.)召開[15-16]。此次會議選擇了包括靜氣動彈性變形的CRM-WB構(gòu)型和CRM-WBPN構(gòu)型作為基準(zhǔn)研究模型，必選工況包括湍流模型的確認(rèn)(Case 1)、固定升力系數(shù)下掛架短艙的安裝阻力及網(wǎng)格收斂性研究(Case 2)、CRM-WB構(gòu)型的靜氣動彈性影響(Case 3)，可選工況包括CRM-WB構(gòu)型的網(wǎng)格自適應(yīng)技術(shù)(Case 4)、CRM-WB構(gòu)型的流固耦合模擬(Case 5)。來自世界各地的25個研究機(jī)構(gòu)采用25種軟件提供了48組Case 1的數(shù)值模擬結(jié)果。其中，32組結(jié)果采用了非結(jié)構(gòu)網(wǎng)格技術(shù),9組結(jié)果采用了結(jié)構(gòu)網(wǎng)格技術(shù),5組結(jié)果采用了笛卡兒網(wǎng)格；36組結(jié)果采用了Spalart-Allmaras一方程湍流模型及其修正形式，6組結(jié)果采用了Menter SST兩方程湍流模型及其修正形式。中國空氣動力研究與發(fā)展中心(CARDC)的王運(yùn)濤(TRIP(TRIsonic Platform)軟件，ID號T1)、陳江濤(Mflow軟件，ID號D1/D2)為此次會議提供了數(shù)值模擬結(jié)果[15]。主要結(jié)論如下：對于Case 2，掛架/短艙的安裝阻力在試驗(yàn)結(jié)果的誤差范圍內(nèi);除個別采用笛卡兒網(wǎng)格的結(jié)果外，網(wǎng)格類型、湍流模型和收斂水平對機(jī)翼和短艙的壓力分布基本沒有影響。對于Case 3，迎角為2.5°時，計算結(jié)果之間的吻合度較好；迎角為4.0°時，計算結(jié)果之間的分散度較大；數(shù)值模擬結(jié)果過度預(yù)測了外翼載荷，導(dǎo)致外翼截面低頭力矩過大、截面升力系數(shù)過大；計算模型中考慮靜氣動彈性變形大大提高了數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的吻合程度；大迎角狀態(tài)下，流動由顯著的激波誘導(dǎo)分離所主導(dǎo)。對于Case 4，模擬三維雷諾平均Navier-Stokes (RANS)方程的自適應(yīng)網(wǎng)格技術(shù)需要進(jìn)一步完善。對于Case 5，流固耦合數(shù)值模擬技術(shù)重要性凸顯，但距離工程實(shí)用尚有差距。

3 CRM計算網(wǎng)格生成

無論對于結(jié)構(gòu)網(wǎng)格還是非結(jié)構(gòu)網(wǎng)格，網(wǎng)格的質(zhì)量、規(guī)模及網(wǎng)格生成技術(shù)本身對于數(shù)值模擬結(jié)果的重要性不言而喻。2003年，Rumsey等[17]采用重疊網(wǎng)格技術(shù)，開展了典型運(yùn)輸機(jī)構(gòu)型抖振邊界附近計算結(jié)果影響因素分析，將局部網(wǎng)格細(xì)節(jié)處理作為一個重要的影響因素。2014年，Slotnick等[18]在總結(jié)過去10年CFD的發(fā)展概況中指出，網(wǎng)格生成和自適應(yīng)仍是CFD工作流程中的重大瓶頸問題之一。

從DPW Ⅰ會議開始，DPW組委會就開始嘗試針對所研究的基準(zhǔn)構(gòu)型開展網(wǎng)格生成規(guī)范研究[19-23]，并從DPW Ⅱ開始引入了網(wǎng)格收斂性研究，作為必選工況之一。經(jīng)歷了DPW Ⅱ和DPW Ⅲ兩次會議的不斷完善，到DPW Ⅳ及以后的會議逐步形成了典型運(yùn)輸機(jī)構(gòu)型的網(wǎng)格生成規(guī)范[24-27]。以CRM翼身組合體構(gòu)型為例(見圖3)，網(wǎng)格生成規(guī)范包括了黏性邊界層物面距離、黏性邊界層內(nèi)的網(wǎng)格增長率、機(jī)翼展向翼根和翼尖的網(wǎng)格分布、機(jī)翼后緣網(wǎng)格數(shù)目、遠(yuǎn)場邊界距離以及不同規(guī)模網(wǎng)格之間網(wǎng)格數(shù)量比例等。為了吸引更多的CFD工作者廣泛參與和最大限度地汲取網(wǎng)格生成經(jīng)驗(yàn)，DPW組委會一方面根據(jù)不斷完善的網(wǎng)格生成規(guī)范提供多種數(shù)據(jù)格式的基礎(chǔ)結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格，另一方面允許參與者根據(jù)自身的工程應(yīng)用經(jīng)驗(yàn)自行生成合適的網(wǎng)格。根據(jù)DPW歷次會議的經(jīng)驗(yàn)和教訓(xùn)不斷總結(jié)網(wǎng)格生成規(guī)范的目的，一是為網(wǎng)格收斂性研究提供系列網(wǎng)格，二是減少采用結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格等不同類型網(wǎng)格數(shù)值模擬結(jié)果之間的差異，三是逐步形成工程實(shí)用的網(wǎng)格生成規(guī)范。

圖3 CRM-WB構(gòu)型表面網(wǎng)格Fig.3 Surface grid of CRM-WB configuration

歷次DPW會議網(wǎng)格生成的一個顯著特點(diǎn)是網(wǎng)格規(guī)模的逐漸增加。以翼身組合體構(gòu)型的中等結(jié)構(gòu)網(wǎng)格為例，2003年DPW Ⅱ組委會提供的DLR-F6翼身組合體構(gòu)型的多塊對接結(jié)構(gòu)網(wǎng)格規(guī)模為390萬，2006年DPW Ⅲ組委會提供的DLR-F6-FX2B翼身組合體構(gòu)型的多塊對接結(jié)構(gòu)網(wǎng)格規(guī)模為810萬，2012年DPW Ⅴ組委會提供的CRM翼身組合體構(gòu)型的重疊結(jié)構(gòu)網(wǎng)格規(guī)模為570萬，2016年DPW Ⅵ組委會提供的CRM翼身組合體構(gòu)型的重疊結(jié)構(gòu)網(wǎng)格規(guī)模為2 470萬。CRM翼身組合體網(wǎng)格規(guī)模從2012—2016年增加了4.3倍，這一方面得益于高性能計算機(jī)的飛速發(fā)展，另一方面也反映了工程應(yīng)用領(lǐng)域?qū)W(wǎng)格規(guī)模無止境的需求。在這個方面， 2010年Sclafani等[24]采用重疊網(wǎng)格技術(shù)和24億的網(wǎng)格模擬了CRM-WBT0構(gòu)型，并得到了具有網(wǎng)格收斂性的數(shù)值模擬結(jié)果。2016年，中國空氣動力研究與發(fā)展中心的TRIP軟件開發(fā)小組完成了CRM-WBPN構(gòu)型139億網(wǎng)格的生成(見圖4)，并在CARDC的高性能計算機(jī)集群上采用6400 CPU核穩(wěn)定運(yùn)行7～24 h，完成12 000步迭代，平均殘差下降4.8個量級，獲得了收斂后的流場計算結(jié)果。

圖4 CRM-WBPN構(gòu)型表面網(wǎng)格(139億)Fig.4 Surface grid of CRM-WBPN configuration (13.9 billion)

總結(jié)DPW Ⅳ～Ⅵ會議網(wǎng)格方面的研究成果，主要結(jié)論如下：網(wǎng)格生成技術(shù)本身沒有質(zhì)的突破；網(wǎng)格規(guī)模達(dá)到一定程度后，結(jié)構(gòu)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格數(shù)值模擬結(jié)果相當(dāng)；如何構(gòu)造網(wǎng)格序列開展網(wǎng)格收斂性研究依然值得進(jìn)一步探索；自適應(yīng)網(wǎng)格技術(shù)是重要的研究方向之一。

4 計算方法與湍流模型

毫無疑問，基于RANS方程的數(shù)值模擬技術(shù)依然是現(xiàn)代飛行器氣動設(shè)計和評估的主要工具，DPW歷次會議上所提供的數(shù)值模擬結(jié)果基本上是采用RANS方程、有限體積方法和二階空間離散精度獲得的。2012年召開的第5屆DPW會議(DPW Ⅴ)，共有來自世界各地的22個研究團(tuán)隊(duì)提供了57組計算結(jié)果，全部采用了基于RANS方程的有限體積方法；2016年召開的第6屆DPW會議(DPW Ⅵ)，共有來自世界各地的25個研究團(tuán)隊(duì)提供了54組計算結(jié)果，其中23個團(tuán)隊(duì)采用了基于RANS方程的有限體積方法。近年來，基于五階空間離散精度的WCNS (Weighted Compact Nonlinear Scheme)[28]格式，通過在幾何守恒律方面持續(xù)不斷的研究工作[29]，在復(fù)雜構(gòu)型數(shù)值模擬方面取得了重要進(jìn)展[30-38]，采用WCNS格式模擬了DLR-F6、DLR-F6/FX2B、CRM-WBH0和CRM-WB等DPW系列會議的翼身組合體構(gòu)型。圖5給出了采用RANS方程、有限差分方法和WCNS格式模擬CRM-WBT0構(gòu)型的表面壓力系數(shù)Cp云圖(Ma=0.85，Re=5.0×106，升力系數(shù)CL=0.500)，所采用的多塊對接結(jié)構(gòu)網(wǎng)格規(guī)模為557萬。上述結(jié)果顯示了WCNS格式在復(fù)雜構(gòu)型方面的應(yīng)用潛力。

湍流模型是影響數(shù)值模擬結(jié)果精準(zhǔn)度的另一個重要因素。Rumsey等[17]開展了典型運(yùn)輸機(jī)構(gòu)型抖振邊界附近計算結(jié)果影響因素分析，研究的影響因素主要包括不同軟件、不同湍流模型和不同局部網(wǎng)格細(xì)節(jié)處理，其中湍流模型的影響最大。不同湍流模型對升力系數(shù)的影響量為3.8%，對阻力系數(shù)的影響量為2.9%，對俯仰力矩的影響量為23.6%。Slotnick等[18]在總結(jié)過去10年CFD的發(fā)展概況中指出，CFD軟件中復(fù)雜湍流模型的有效性和收斂性是工程應(yīng)用中的重大瓶頸問題之一。正是逐步意識到了湍流模型對數(shù)值模擬結(jié)果的極端重要性，從DPW Ⅴ會議開始，DPW組委會專門選擇了湍流模型的驗(yàn)證工況，用于驗(yàn)證各種CFD軟件中湍流模型的正確性及工程適用性，并從DPW V會議的可選工況(Case 3)上升到DPW Ⅵ會議的必選工況(Case 1)，具體工作可參考Levy[14]和Roy[39]等的文獻(xiàn)。

圖5 CRM-WBT0構(gòu)型高階精度模擬Fig.5 High-order precision simulation of CRM-WBT0 configuration

總結(jié)DPW歷次會議上所提供的數(shù)值模擬結(jié)果可以看出，Spalart-Allmaras一方程湍流模型[40]及其各種修正形式和Menter SST兩方程湍流模型[41]及其各種修正形式依然是目前工程應(yīng)用領(lǐng)域主要采用的湍流模型。DPW Ⅳ會議上，17個研究團(tuán)隊(duì)中提供的28組結(jié)果中，18組結(jié)果采用了S-A一方程湍流模型，7組結(jié)果采用了SST兩方程湍流模型；DPW Ⅴ會議上，22個研究團(tuán)隊(duì)中提供的57組結(jié)果中，38組結(jié)果采用了S-A一方程湍流模型，13組結(jié)果采用了SSTk-ε兩方程湍流模型；DPW Ⅵ會議上，25個研究團(tuán)隊(duì)中提供的54組結(jié)果中，36組結(jié)果采用了S-A一方程湍流模型，8組結(jié)果采用了SSTk-ε兩方程湍流模型。

圖6 CRM分離泡尺寸隨迎角的變化[14] Fig.6 CRM separation bubble size vs angle of attack[14]

2004年RAND公司的研究報告中指出[42]，在飛行器設(shè)計相關(guān)的8個復(fù)雜流動機(jī)理研究方面，黏性流動分離的起始與發(fā)展是CFD模擬尚未很好解決的6個問題之一。圖6給出了DPW Ⅴ會議上，Ma=0.85，Re=5.0×106，CL=0.500時，部分?jǐn)?shù)值模擬結(jié)果得到的CRM-WB構(gòu)型翼身結(jié)合部后緣局部分離區(qū)域的大小隨迎角的變化[14]。其中，不同顏色的點(diǎn)代表了不同迎角下的數(shù)值模擬結(jié)果，WBL(Wing Butt Line coordinator)代表機(jī)翼在Y軸方向的坐標(biāo)。匯總DPW Ⅳ～Ⅵ會議上CRM構(gòu)型翼身結(jié)合部后緣局部分離流動數(shù)值模擬結(jié)果可以看出，在這個局部區(qū)域的流動模擬是一個難點(diǎn)。不僅數(shù)值模擬結(jié)果之間差異明顯，而且數(shù)值模擬結(jié)果與相應(yīng)的試驗(yàn)結(jié)果之間同樣差異明顯(試驗(yàn)結(jié)果在此種工況下沒有明顯的分離)，更嚴(yán)重的情況是由于這個區(qū)域存在的虛假分離會導(dǎo)致計算結(jié)果很難收斂。

湍流模型的研究進(jìn)展主要體現(xiàn)在包含QCR[43](Quadratic Constitutive Relation)修正的S-A一方程湍流模型或SST兩方程湍流模型的應(yīng)用研究。2012年，Yamamoto等[44]采用引入QCR修正的S-A一方程和SST湍流模型模擬了CRM-WBH0構(gòu)型繞流流場，計算結(jié)果顯示，采用QCR修正有效減弱了翼身結(jié)合部的虛假分離。Sclafani等[24]采用引入QCR修正的S-A模型模擬了DPW V提供的CRM-WB構(gòu)型繞流流場，計算結(jié)果顯示，采用QCR修正有效提高了翼身結(jié)合部的模擬精度(見圖7[24],Ma=0.85，Re=5.0×106，α=4°)。在以上研究工作的推動下，2016年召開的DPW VI會議上，在36組采用S-A湍流模型的結(jié)果中，有15組采用了QCR修正。

總結(jié)DPW Ⅳ～Ⅵ會議湍流模型方面的研究成果，主要結(jié)論如下：S-A和SSTk-ε湍流模型是目前工程應(yīng)用的主導(dǎo)；在S-A和SSTk-ε湍流模型中引入QCR修正有效提高了翼身結(jié)合部的流動模擬精度；GoldbergRT、EARSM(Explicit Algebraic Reynolds Stress Model)、RSM (Reynolds Stress Model )等湍流模型值得進(jìn)一步開展研究。

圖7 QCR修正對表面流線的影響 [24] Fig.7 Effect of QCR correction on surface streamline [24]

5 計算結(jié)果與試驗(yàn)結(jié)果對比

將計算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果相比較，是驗(yàn)證數(shù)值模擬結(jié)果精準(zhǔn)度的重要手段，也是CFD確認(rèn)工作的主要內(nèi)容。必須指出的是應(yīng)用于CFD確認(rèn)工作的風(fēng)洞試驗(yàn)與一般型號的風(fēng)洞試驗(yàn)有很大不同，應(yīng)用于CFD確認(rèn)工作風(fēng)洞試驗(yàn)的顧客是CFD工作者，Oberkampf和Trucanob[45]提出了設(shè)計CFD確認(rèn)試驗(yàn)的7項(xiàng)指導(dǎo)原則。在不同風(fēng)洞中、不同時期開展的CRM風(fēng)洞試驗(yàn)結(jié)果為典型構(gòu)型運(yùn)輸機(jī)高速構(gòu)型的CFD確認(rèn)工作提供了豐富的風(fēng)洞試驗(yàn)數(shù)據(jù)，但在試驗(yàn)數(shù)據(jù)的分析和整理方面還需要進(jìn)一步開展工作。

Levy等[13-14]總結(jié)了CRM風(fēng)洞試驗(yàn)與數(shù)值模擬的不同點(diǎn)，主要包括：CRM風(fēng)洞試驗(yàn)存在洞壁干擾，而CFD模擬采用了自由來流條件；CRM風(fēng)洞試驗(yàn)存在支撐干擾，而CFD模擬不包含支撐裝置；風(fēng)洞試驗(yàn)采用固定轉(zhuǎn)捩方式，而CFD模擬一般采用全湍流模擬方式；風(fēng)洞試驗(yàn)結(jié)果包含了靜氣動彈性變形影響，而數(shù)值模擬采用了剛性1g外形；風(fēng)洞試驗(yàn)結(jié)果與數(shù)值模擬結(jié)果均存在不確定性和誤差；風(fēng)洞試驗(yàn)結(jié)果對已知的影響因素做了各種修正(支撐干擾除外)，而數(shù)值模擬結(jié)果沒有做修正。上述6個方面將對數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比分析產(chǎn)生重要影響。以下將重點(diǎn)從固定升力系數(shù)下的氣動特性、壓力分布對比，CRM模型掛架/短艙的安裝阻力對比，氣動特性隨迎角的變化，流固耦合數(shù)值模擬等5個方面歸納總結(jié)DPW Ⅳ～Ⅵ會議上數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比情況。

1) 氣動特性對比(CL≈0.500)

為鼓勵更多的CFD工作者參與會議，歷次DPW組委會均按預(yù)先約定的網(wǎng)格規(guī)范提供了至少3套不同規(guī)模的粗(Coarse)、中(Medium)、細(xì)(Fine)網(wǎng)格用以開展網(wǎng)格收斂性研究，其中，中等網(wǎng)格是目前工程應(yīng)用方面普遍采用的網(wǎng)格規(guī)模，以下各節(jié)數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的對比均以中等網(wǎng)格的統(tǒng)計分析結(jié)果為主。

DPW Ⅴ和DPW Ⅵ兩次會議上的CRM-WB構(gòu)型不同點(diǎn)在于，采用ETW風(fēng)洞試驗(yàn)數(shù)據(jù)，DPW Ⅵ會議提供的CRM-WB構(gòu)型包含了不同來流迎角下的靜氣動彈性變形。圖8給了CRM-WB構(gòu)型不同迎角下靜氣動彈性變形風(fēng)洞試驗(yàn)測量結(jié)果，其中，橫坐標(biāo)為機(jī)翼展向的不同站位，縱坐標(biāo)為靜氣動彈性導(dǎo)致的扭轉(zhuǎn)角，升力系數(shù)為0.5的數(shù)據(jù)是插值得到的；不同顏色的曲線代表了不同來流迎角下測量得到的靜氣動彈性變形(A2.5即迎角為2.5°)。基于風(fēng)洞試驗(yàn)中測量得到的風(fēng)洞模型靜氣動彈性變形數(shù)據(jù)，David[46]研究了CRM翼身組合體模型靜氣動彈性變形對數(shù)值模擬結(jié)果的影響；Keye等[47-48]采用流固耦合方法研究了靜氣動彈性變形對CRM翼身組合體模型數(shù)值模擬結(jié)果的影響；這些工作直接促成了DPW Ⅵ組委會在CRM-WB計算模型中考慮了靜氣動彈性變形的影響。

圖8 由風(fēng)洞測量數(shù)據(jù)獲得的靜氣動彈性扭轉(zhuǎn)分布[15] Fig.8 Distribution of static aeroelastic twist derived from wing tunnel measurement data[15]

表1給出了DPW Ⅴ和DPW Ⅵ會議上CRM-WB構(gòu)型的“core solutions”統(tǒng)計結(jié)果[49-50]，計算狀態(tài)為Ma=0.85、Re=5.0×106、CL=0.500±0.001。其中，“Median”代表各家數(shù)值模擬結(jié)果的中位數(shù)，“Deviation”代表統(tǒng)計結(jié)果的均方誤差，CL為升力系數(shù),CD為阻力系數(shù),CDp為壓差阻力系數(shù),CDf為摩擦阻力系數(shù),Cm為俯仰力矩系數(shù)。

從表1可以得到如下基本結(jié)論：① 網(wǎng)格規(guī)模的增加并沒有降低統(tǒng)計分析結(jié)果的均方誤差，反而普遍有所增加(DPW Ⅵ組委會提供的基礎(chǔ)網(wǎng)格規(guī)模是DPW Ⅴ組委會提供的基礎(chǔ)網(wǎng)格規(guī)模的4.3倍), 如阻力系數(shù)的均方誤差由DPW Ⅴ統(tǒng)計結(jié)果的3.2 counts增加到DPW Ⅵ統(tǒng)計結(jié)果的4.4 counts; ② 數(shù)值模擬得到的阻力系數(shù)散布度與試驗(yàn)結(jié)果之間差量相當(dāng)，來流迎角與俯仰力矩系數(shù)的散布度大于試驗(yàn)結(jié)果之間的差量。DPW Ⅵ會議上，阻力系數(shù)的均方誤差為4.4 counts，略小于風(fēng)洞試驗(yàn)結(jié)果之間7.5 counts的阻力差異；③ 來流迎角的均方誤差是試驗(yàn)結(jié)果之間差量的兩倍，俯仰的均方誤差比試驗(yàn)結(jié)果之間的差量大了一個量級; ④ 計算模型中考慮了靜氣動彈性變形后，氣動特性的絕對量與試驗(yàn)結(jié)果之間的差異明顯減少。

表1 CRM-WB構(gòu)型數(shù)值模擬結(jié)果的中位數(shù)及均方差和風(fēng)洞試驗(yàn)結(jié)果[49-50]

圖9 包含/不包含模型支撐裝置壓力系數(shù)差量云圖[52]Fig.9 ΔCp contours with/without support system of model[52]

由表1可以看出，計算模型中考慮了靜氣動彈性變形后，數(shù)值模擬得到的來流迎角、阻力系數(shù)和俯仰力矩系數(shù)與試驗(yàn)結(jié)果仍有較大差異，其主要原因是計算模型中沒有考慮支撐裝置與固定轉(zhuǎn)捩位置的影響。有關(guān)CRM支撐干擾的研究工作最早是由NASA Langley研究中心的Rivers等[51-52]完成的。文獻(xiàn)[51] 采用USM3D軟件研究了CRM-WBT0模型支撐裝置對氣動特性數(shù)值模擬結(jié)果的影響，文獻(xiàn)[52]在上述工作的基礎(chǔ)上，又進(jìn)一步考慮了CRM-WBT0模型靜氣動彈性變形對數(shù)值模擬結(jié)果的影響。Rivers等的數(shù)值模擬結(jié)果表明(見圖9，ΔCp為無量綱壓力系數(shù)差量)，在Ma=0.85、Re=5.0×106、α=2°條件下，在計算模型中同時考慮靜氣動彈性變形和模型支撐裝置，使得升力系數(shù)下降0.026 0、阻力系數(shù)下降3.7 counts、俯仰力矩系數(shù)增加0.035 9，數(shù)值模擬結(jié)果更接近風(fēng)洞試驗(yàn)結(jié)果。針對CRM-WBH0構(gòu)型，Zilliac等[53]采用OVERFLOW軟件研究了固定轉(zhuǎn)捩位置對阻力系數(shù)數(shù)值模擬結(jié)果的影響。研究表明：在Ma=0.85、Re=5.0×106、α=2°條件下，在風(fēng)洞試驗(yàn)機(jī)翼的10%固定轉(zhuǎn)捩位置前，采用“全湍流”模擬方式模擬得到的阻力系數(shù)大了6.4 counts，占總阻力系數(shù)的2.3%。

2) 壓力分布對比(CL≈0.500)

氣動特性數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比屬于宏觀量的對比，壓力分布數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比則屬于微觀量的對比。壓力分布計算結(jié)果與試驗(yàn)結(jié)果的對比可以進(jìn)一步揭示二者之間升力和力矩等宏觀量差異的原因。

圖10 來自DPW Ⅴ的壓力系數(shù)分布對比[13]Fig.10 Comparison of pressure coefficient distributions from DPW Ⅴ[13]

圖10給出了DPW Ⅴ會議上，CRM-WB構(gòu)型機(jī)翼η=0.727、0.950站位上，Ma=0.85、Re=5.0×106、CL=0.500條件下，采用中等網(wǎng)格獲得的壓力系數(shù)計算結(jié)果與相鄰迎角NTF測壓風(fēng)洞之間的比較[13]，其中，曲線代表采用不同網(wǎng)格形式獲得的計算結(jié)果，黑色點(diǎn)代表相近升力系數(shù)下的測壓試驗(yàn)結(jié)果。由DPW Ⅴ會議上壓力分布的統(tǒng)計結(jié)果可以看出，在機(jī)翼η=0.50站位以內(nèi)，數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果之間吻合良好；隨著機(jī)翼站位向翼梢方向移動，數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果之間吻合度逐漸變差，主要表現(xiàn)在：相對于風(fēng)洞試驗(yàn)結(jié)果，數(shù)值模擬得到的機(jī)翼上翼面激波位置普遍靠后；翼梢處，數(shù)值模擬得到的機(jī)翼前緣壓力分布普遍高于風(fēng)洞試驗(yàn)結(jié)果。導(dǎo)致上述差異的主要原因是DPW Ⅴ會議所提供CRM-WB構(gòu)型沒有包含風(fēng)洞試驗(yàn)中模型靜氣動彈性變形的影響。

圖11 來自DPW Ⅵ的壓力系數(shù)分布對比[14]Fig.11 Comparison of pressure coefficient distributions from DPW Ⅵ[14]

與DPW Ⅴ會議所提供CRM-WB構(gòu)型不同，DPW Ⅵ會議所提供的CRM-WB構(gòu)型包含了ETW中測量得到的靜氣動變形。圖11給出了DPW Ⅵ會議上，CRM-WB構(gòu)型機(jī)翼η=0.603、0.950站位上，Ma=0.85、Re=5.0×106、CL=0.500條件下，采用密網(wǎng)格獲得的壓力系數(shù)計算結(jié)果與相近升力系數(shù)Ames風(fēng)洞測壓結(jié)果之間的比較[14]，其中，不同顏色曲線代表采用不同湍流模型獲得的計算結(jié)果，藍(lán)色點(diǎn)代表相應(yīng)的試驗(yàn)結(jié)果。由DPW Ⅵ會議上壓力分布的統(tǒng)計結(jié)果可以看出，考慮了靜氣動彈性變形后，在機(jī)翼η=0.500站位以內(nèi)，數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果之間吻合良好；在機(jī)翼η=0.603站位上，數(shù)值模擬得到的機(jī)翼上翼面激波位置明顯比相應(yīng)的試驗(yàn)結(jié)果靠后；在機(jī)翼η=0.950站位上，數(shù)值模擬得到的壓力分布與風(fēng)洞測壓結(jié)果之間的吻合程度大大改善。導(dǎo)致η=0.603站位上數(shù)值模擬與試驗(yàn)之間激波位置差異的主要原因是DPW Ⅵ會議所提供CRM-WB構(gòu)型的數(shù)值模擬結(jié)果沒有包含風(fēng)洞模型支撐裝置的影響。

基于DPW Ⅵ會議所提供的CRM-WB構(gòu)型、CRM-WB風(fēng)洞試驗(yàn)?zāi)Ｐ偷挠邢拊Ｐ秃蚇ASA的CRM模型官方網(wǎng)站提供的NTF風(fēng)洞模型支撐數(shù)據(jù)，中國空氣動力研究與發(fā)展中心的TRIP軟件開發(fā)小組數(shù)值模擬了同時包含靜氣動彈性變形和模型支撐裝置的CRM-WB構(gòu)型氣動特性[54]。圖12給出了機(jī)翼η=0.727、0.950站位上，Ma=0.85、Re=5.0×106、α=2.75°條件下，CRM-WB構(gòu)型的CFD計算結(jié)果(CRM-WB_CFD)、包含支撐裝置的CRM-WB構(gòu)型的CFD計算結(jié)果(CRM-WBS_CFD)和包含支撐裝置的CRM-WB構(gòu)型的流固耦合計算結(jié)果(CRM-WBS_FSC)以及相鄰迎角下NTF測壓試驗(yàn)結(jié)果。其中橫坐標(biāo)x/c為流向無量綱距離。比較CRM-WB_CFD和CRM-WBS_CFD的結(jié)果可以看出，模型支撐裝置的影響主要使得機(jī)翼上翼面的激波位置前移。比較CRM-WBS_CFD和CRM-WB_FSC的結(jié)果可以看出，計算模型中考慮機(jī)翼的靜氣動彈性變形，不僅使得η=0.727站位上的激波位置進(jìn)一步前移，而且顯著降低了機(jī)翼上翼面激波位置前的壓力系數(shù)；η=0.950站位上，靜氣動彈性變形使得激波位置后移，并且顯著降低了機(jī)翼上翼面激波位置前的壓力系數(shù)?？傊C合考慮了機(jī)翼靜氣動彈性變形和模型支撐裝置后，數(shù)值模擬得到了壓力系數(shù)分布更加接近試驗(yàn)結(jié)果。

3) 掛架短艙的安裝阻力對比(CL=0.500)

評估現(xiàn)代CFD技術(shù)模擬計算構(gòu)型變化導(dǎo)致的氣動特性變化的能力，一直是DPW系列會議的研究內(nèi)容之一。DPW Ⅱ會議上評估了DLR-F6構(gòu)型掛架/短艙安裝阻力的模擬能力[2，55]，DPW Ⅲ會議上評估了DLR-F6構(gòu)型翼身結(jié)合部局部修型導(dǎo)致的氣動特性微小變化的模擬能力[3，37，56]。DPW Ⅵ會議上開展了CRM構(gòu)型掛架/短艙安裝阻力數(shù)值模擬能力的評估。

表2給出了文獻(xiàn)[50]統(tǒng)計得到的CRM-WBPN和CRM-WB兩種構(gòu)型氣動特性差量的中位數(shù)及均方差，同時給出了NTF和Ames風(fēng)洞試驗(yàn)結(jié)果，計算狀態(tài)為Ma=0.85、Re=5.0×106、CL=0.500。由表中可以看出，DPW會議上各數(shù)值模擬結(jié)果得到的來流迎角差量、阻力系數(shù)差量與相應(yīng)的試驗(yàn)結(jié)果吻合良好，數(shù)值模擬結(jié)果之間的散布度較小，阻力系數(shù)的差量主要來自摩擦阻力系數(shù)；俯仰力矩系數(shù)差量與試驗(yàn)結(jié)果有一定差距，但各種軟件數(shù)值模擬結(jié)果之間的一致性較好。

圖12 CRM-WB構(gòu)型的壓力系數(shù)分布對比Fig.12 Comparison of pressure coefficient distributions of CRM-WB configuration

表2 中等網(wǎng)格ΔCD計算結(jié)果的中位數(shù)及均方差和風(fēng)洞試驗(yàn)結(jié)果[50]

Table 2 Core medians and standard deviations of numerical solutions for ΔCD with medium grid and wind tunnel test results[50]

SourceΔα/(°)ΔCDΔCDpΔCDfΔCmDPWⅥMedian0.1640.002290.000560.001690.0056Deviation0.0070.000130.000130.000050.0007TestNTF-1970.1500.002300.0100Ames-2160.1510.002300.0100

與DPW Ⅱ和DPW Ⅲ會議上統(tǒng)計分析所得到的結(jié)論一致。目前的CFD軟件可以較好地模擬由于計算構(gòu)型的變化引起的氣動特性差量。

4) 氣動特性隨迎角的變化(Ma=0.85)

評估現(xiàn)代CFD技術(shù)模擬典型運(yùn)輸機(jī)構(gòu)型氣動特性隨迎角變化的能力，是CFD確認(rèn)工作的主要內(nèi)容，也始終是DPW系列會議的重要研究內(nèi)容之一。

DPW Ⅴ會議上采用了CRM-WB構(gòu)型作為基準(zhǔn)研究模型， 圖14給出了DPW Ⅴ會議上各家參與單位給出的升力系數(shù)、力矩系數(shù)與試驗(yàn)結(jié)果的比較[14]。其中， 不同顏色的曲線代表各家參與單位給出的數(shù)值模擬結(jié)果，空心圓圈標(biāo)記代表NTF、TWT風(fēng)洞的試驗(yàn)結(jié)果，實(shí)心圓圈標(biāo)記代表了“偽試驗(yàn)結(jié)果”(Pseudo Wind Tunnel Data，因?yàn)镽ivers等考慮了靜氣動彈性氣動特性數(shù)值模擬結(jié)果，根據(jù)在風(fēng)洞試驗(yàn)結(jié)果中扣除了氣動彈性變形影響后得到的結(jié)果)。本次會議上數(shù)值模擬得到的氣動特性與相應(yīng)風(fēng)洞試驗(yàn)結(jié)果的對比是非常令人失望的[13]，這種現(xiàn)象在DPW Ⅳ會議上同樣存在[51-52]。主要表現(xiàn)在：相同迎角下，數(shù)值模擬得到的升力系數(shù)普遍大于試驗(yàn)結(jié)果；相同升力系數(shù)下，數(shù)值模擬得到的低頭力矩系數(shù)普遍大于試驗(yàn)結(jié)果。主要受Rivers等[51-52]有關(guān)CRM-WBH0模型靜氣動彈性變形和模型支撐裝置對氣動特性影響的研究工作啟發(fā)， DPW Ⅴ組委會構(gòu)造了“偽試驗(yàn)結(jié)果”，以期在相近構(gòu)型的基礎(chǔ)上，開展數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比。

圖13 CRM-WBH0構(gòu)型的極曲線[12]Fig.13 Polars of CRM-WBH0 configuration[12]

圖14 CRM-WB構(gòu)型的升力系數(shù)和力矩系數(shù)曲線(DPW Ⅴ)[14]Fig.14 Lift and pitching moment coefficients curves of CRM-WB configuration (DPW Ⅴ)[14]

由圖14可以看出，部分?jǐn)?shù)值模擬結(jié)果在迎角3.0°～4.0°之間升力系數(shù)出現(xiàn)了過早的下降，這明顯與相應(yīng)的試驗(yàn)結(jié)果不符。其根本原因是，這些數(shù)值模擬結(jié)果在相應(yīng)的迎角狀態(tài)下，在翼身結(jié)合部和機(jī)翼后緣出現(xiàn)了大范圍的虛假分離。剔除這部分結(jié)果后(47組結(jié)果中剩余26組結(jié)果)，數(shù)值模擬結(jié)果與“偽試驗(yàn)結(jié)果”之間的比較見圖15[14]。其中，不同顏色曲線代表部分參與單位采用不同湍流模型得到的數(shù)值模擬結(jié)果，實(shí)心圓圈標(biāo)記代表了“偽試驗(yàn)結(jié)果”。由圖可以看出，計算結(jié)果之間的數(shù)據(jù)散布度隨著迎角增加而逐漸增加，在迎角為4.0°時，升力系數(shù)的數(shù)值模擬結(jié)果之間相差0.055，力矩系數(shù)的數(shù)值模擬結(jié)果之間相差0.043。造成這種現(xiàn)象的主要原因是不同的數(shù)值模擬結(jié)果在機(jī)翼上翼面激波位置和機(jī)翼后緣分離區(qū)的模擬上存在明顯的差異。

圖15 CRM-WB構(gòu)型的升力系數(shù)和力矩系數(shù)曲線(DPW Ⅴ，剔除不合理結(jié)果)[14]Fig.15 Lift and pitching moment coefficients curves of CRM-WB configuration (DPW Ⅴ, minus outliers)[14]

圖16[14]給出了數(shù)值模擬得到的極曲線與NTF、Ames試驗(yàn)結(jié)果的對比。其中，不同顏色曲線代表部分參與單位采用不同湍流模型得到的數(shù)值模擬結(jié)果，空心標(biāo)記代表NTF、TWT風(fēng)洞的試驗(yàn)結(jié)果。由圖看出，相同升力系數(shù)下，阻力系數(shù)數(shù)值模擬結(jié)果的散布度隨迎角的增加而增加；升力系數(shù)小于0.6時，阻力系數(shù)數(shù)值模擬結(jié)果的散布度與風(fēng)洞試驗(yàn)結(jié)果之間的差量相仿。

圖16 CRM-WB構(gòu)型的極曲線(DPW Ⅴ)[14]Fig.16 Polars of CRM-WB configuration (DPW Ⅴ) [14]

汲取DPW Ⅳ、DPW Ⅴ兩次會議的經(jīng)驗(yàn)和教訓(xùn)，采用ETW風(fēng)洞試驗(yàn)的變形測量數(shù)據(jù)，DPW Ⅵ組委會選擇了包含靜氣動彈性變形的CRM-WB構(gòu)型作為基準(zhǔn)研究模型。需要說明的是，風(fēng)洞試驗(yàn)中，在不同的迎角下，機(jī)翼所承受的氣動載荷不同，機(jī)翼的靜氣動彈性變形也不同(圖8)。因此，DPW Ⅵ會議上的Case 3本質(zhì)上包括了7種構(gòu)型。

圖17[15]給出了DPW Ⅵ會議上所有計算結(jié)果的升力系數(shù)和俯仰力矩系數(shù)曲線與風(fēng)洞試驗(yàn)結(jié)果的對比。其中，不同顏色曲線代表參與單位采用不同湍流模型得到的數(shù)值模擬結(jié)果，空心與實(shí)心的圓圈標(biāo)記分別代表NTF、TWT風(fēng)洞試驗(yàn)結(jié)果。由圖17可以看出，與DPW Ⅴ的相應(yīng)結(jié)果相比較(圖14)，計算模型中考慮了靜氣動彈性變形后，氣動特性數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的吻合程度顯著提高，但部分?jǐn)?shù)值模擬結(jié)果依然在迎角為3.0°～4.0°之間出現(xiàn)了升力系數(shù)過早的下降。升力系數(shù)過早失速問題也是DPW歷次會議重點(diǎn)研究的問題之一，討論的焦點(diǎn)主要集中在：較大迎角下，RANS方程是否適用？是否URANS方程或則DES方法更適合描述這種流動？剔除這部分結(jié)果后(40組結(jié)果中剩余21組結(jié)果)，數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果之間的比較見圖18[15]和圖19[15]。

圖17 CRM-WB構(gòu)型的升力系數(shù)和力矩系數(shù)曲線(DPW Ⅵ)[15]Fig.17 Lift and pitching moment coefficients curves of CRM-WB configuration(DPW Ⅵ) [15]

圖18 CRM-WB構(gòu)型的升力系數(shù)和力矩系數(shù)曲線(DPW Ⅵ，剔除不合理結(jié)果)[15]Fig.18 Lift and pitching moment coefficients curves of CRM-WB configuration (DPW Ⅵ, minus outliers)[15]

圖19 CRM-WB構(gòu)型的極曲線(DPW Ⅵ)[15]Fig.19 Polars of CRM-WB configuration (DPW Ⅵ)[15]

圖18給出了剔除“outliers”后，數(shù)值模擬得到的升力系數(shù)與俯仰力矩系數(shù)曲線與NTF和Ames試驗(yàn)結(jié)果之間的比較，圖中的各種標(biāo)識與圖17相同。由圖可以看出，與DPW V的統(tǒng)計結(jié)果相類似，數(shù)值模擬結(jié)果之間的數(shù)據(jù)散布度依然隨著迎角增加而逐漸增加，在迎角為4.0°時，升力系數(shù)的數(shù)值模擬結(jié)果之間相差0.063，力矩系數(shù)的數(shù)值模擬結(jié)果之間相差0.045。讓人沮喪的是，從DPW Ⅴ～DPW Ⅵ，雖然中等網(wǎng)格的規(guī)模增加4倍多，但4°迎角條件下，數(shù)值模擬結(jié)果之間的散布度不僅沒有減少，反而略有增加。圖19給出了數(shù)值模擬得到的極曲線與NTF、Ames試驗(yàn)結(jié)果的對比，圖中的各種標(biāo)識與圖17相同。由圖看出，與DPW Ⅴ的統(tǒng)計結(jié)果相類似，相同升力系數(shù)下，阻力系數(shù)數(shù)值模擬結(jié)果的散布度隨迎角的增加而增加；升力系數(shù)小于0.6時，阻力系數(shù)數(shù)值模擬結(jié)果的散布度與風(fēng)洞試驗(yàn)結(jié)果之間的差量基本相當(dāng)。

5) 流固耦合數(shù)值模擬(Ma=0.85)

靜氣動彈性變形對典型運(yùn)輸機(jī)構(gòu)型氣動特性的影響一直受到氣動工作者的關(guān)注[57-59]，2012年AIAA發(fā)起的第1屆氣動彈性預(yù)測研討會AePW (Aeroelastic Prediction Workshop)[60]，同樣包含了靜氣動彈性變形預(yù)測研究。目前，借助CFD手段研究靜氣動彈性變形對氣動特性數(shù)值模擬結(jié)果的影響一般采用兩種方式:① 通過風(fēng)洞試驗(yàn)測量得到模型在氣動載荷作用下的靜氣動彈性變形，通過幾何重構(gòu)方法構(gòu)建出模型變形后的外形，利用CFD方法計算變形前后外形的氣動力差異；② 采用風(fēng)洞試驗(yàn)?zāi)Ｐ偷臍鈩油庑魏拖鄳?yīng)的結(jié)構(gòu)有限元模型，通過氣動/結(jié)構(gòu)耦合數(shù)值模擬方法研究靜氣動彈性變形對氣動特性數(shù)值模擬結(jié)果的影響。

采用風(fēng)洞試驗(yàn)測量得到的靜氣動彈性變形，Rivers等[51-52]開展了CRM-WBT0模型支撐裝置和靜氣動彈性變形對數(shù)值模擬結(jié)果的影響，David[46]研究了CRM-WB模型靜氣動彈性變形對數(shù)值模擬結(jié)果的影響。Keye等[48]采用流固耦合方法研究了靜氣動彈性變形對CRM-WB模型數(shù)值模擬結(jié)果的影響。DPW Ⅵ會議上[15]的研究工況Cases 2～4本質(zhì)上是采用風(fēng)洞試驗(yàn)測量得到的靜氣動彈性變形重構(gòu)CRM-WB構(gòu)型的幾何數(shù)模，進(jìn)而研究靜氣動變形對氣動特性數(shù)值模擬結(jié)果的影響；研究工況Case 5是采用流固耦合方法研究靜氣動彈性變形對氣動特性數(shù)值模擬結(jié)果的影響。Keye和Mavriplis[16]總結(jié)了DPW Ⅵ會議上4組流固耦合數(shù)值模擬結(jié)果(圖20)。其中，J4、V5代表CFD++軟件提供的結(jié)果，L2代表Tau軟件提供的結(jié)果，T1代表TRIP軟件提供的結(jié)果。主要結(jié)論如下：4組流固耦合數(shù)值模擬結(jié)果之間差異很小，采用流固耦合方法得到的氣動特性與壓力分布結(jié)果合理，靜氣動彈性數(shù)值模擬方法遠(yuǎn)未成熟。

從圖18可以看出，計算模型中考慮了靜氣動彈性變形影響后，CRM-WB模型氣動特性數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果之間的吻合程度有了明顯的改善，但依然存在明顯差異。計算模型中沒有考慮支撐裝置影響是主要原因之一。利用DPW組委會提供的CRM-WB計算模型、結(jié)構(gòu)有限元模型(圖21)和支撐模型，采用結(jié)構(gòu)網(wǎng)格技術(shù)(圖22)，TRIP軟件開發(fā)小組[56]采用流固耦合計算方法開展了包括支撐裝置的CRM-WB構(gòu)型(CRM-WBS)數(shù)值模擬，詳細(xì)研究了風(fēng)洞模型支撐裝置和靜氣動彈性變形對CRM-WB構(gòu)型氣動特性數(shù)值模擬結(jié)果的影響。

圖20 CRM-WB構(gòu)型流固耦合計算結(jié)果[16]Fig.20 Simulation results of fluid/structure coupling of CRM-WB configuration[16]

圖21 CRM-WB風(fēng)洞試驗(yàn)有限元模型Fig.21 Finite element model for CRM-WB wind tunnel test

圖22 CRM-WBS構(gòu)型網(wǎng)格拓?fù)浜蛯ΨQ面網(wǎng)格 (局部)Fig.22 Grid topology of CRM-WBS configuration and grids at symmetric plane (local)

圖23 CRM 翼身組合體構(gòu)型的氣動特性Fig.23 Aerodynamic characteristics of CRM wing-body configuration

圖23給出了CRM-WB構(gòu)型CFD數(shù)值模擬結(jié)果(CRM-WB_CFD)、CRM-WBS構(gòu)型CFD數(shù)值模擬結(jié)果(CRM-WBS_CFD)和流固耦合(CRM-WBS_FSC)數(shù)值模擬結(jié)果，同時還給出了NTF測力試驗(yàn)結(jié)果。來流條件為：Ma=0.85，Re=5.0×106，α=0°～4.00°；流固耦合計算時，速壓q=61.062 kPa，載荷因子q/E=3.342×10-7。對比CRM-WBS構(gòu)型與CRM-WB構(gòu)型的CFD數(shù)值模擬結(jié)果可以看出，在α≤3.75°范圍內(nèi)，模型支撐裝置使得升力系數(shù)、阻力系數(shù)下降，低頭力矩減少。對比CRM-WBS構(gòu)型的CFD結(jié)果與CRM-WBS構(gòu)型的流固耦合數(shù)值模擬結(jié)果可以看出，在計算迎角范圍內(nèi)，靜氣動彈性變形使得升力系數(shù)、阻力系數(shù)進(jìn)一步下降，低頭力矩進(jìn)一步減少。采用流固耦合方法得到的CRM-WBS構(gòu)型升力系數(shù)和阻力系數(shù)的計算結(jié)果更加接近試驗(yàn)值；俯仰力矩系數(shù)計算結(jié)果與試驗(yàn)結(jié)果的吻合程度得到進(jìn)一步改善，但依然有一定的差距。另外值得關(guān)注的一點(diǎn)是，不包含支撐裝置的CRM-WB_CFD數(shù)值模擬結(jié)果失速迎角在3.75°，而包含支撐裝置的CRM-WBS_CFD與CRM-WBS_FSC的數(shù)值模擬結(jié)果直到迎角4.00°還沒有出現(xiàn)失速。數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果之間的差異需要從風(fēng)洞試驗(yàn)數(shù)據(jù)的修正和數(shù)值計算方法兩個方面進(jìn)一步開展研究工作。

6 CRM的優(yōu)化設(shè)計

由于DPW系列會議持續(xù)擴(kuò)大的影響力，與CRM相關(guān)的研究工作也陸續(xù)展開。其中，CRM的氣動優(yōu)化設(shè)計、氣動/結(jié)構(gòu)優(yōu)化設(shè)計是當(dāng)前非?；钴S的一個研究方向。

Lyu等[61]基于RANS方程求解器和離散伴隨方法算法，開展了CRM機(jī)翼在多種約束條件下的單點(diǎn)優(yōu)化設(shè)計、多點(diǎn)優(yōu)化設(shè)計工作(圖24[61])，相對于原始外形，單點(diǎn)氣動優(yōu)化設(shè)計使得阻力系數(shù)降低8.5%。Kenway和Martins等[62]采用RANS方程求解器和梯度優(yōu)化算法，完成了CRM-WBH構(gòu)型考慮抖振起始約束的多個單點(diǎn)和多點(diǎn)優(yōu)化算例(圖25[62])，獲得了較好的氣動優(yōu)化效果。陳頌等[63]針對CRM-WBH構(gòu)型開展了基于離散伴隨技術(shù)的氣動外形優(yōu)化設(shè)計方法研究，進(jìn)行了有/無力矩配平約束的優(yōu)化設(shè)計。劉峰博等[64]發(fā)展了基于離散伴隨方法的梯度計算模塊，對CRM翼身組合體構(gòu)型進(jìn)行了單點(diǎn)減阻優(yōu)化設(shè)計。

圖24 CRM機(jī)翼單點(diǎn)和多點(diǎn)氣動外形優(yōu)化[61]Fig.24 Single-point and multi-point aerodynamic shape optimization of CRM wing[61]

圖25 CRM-WBH構(gòu)型單點(diǎn)和多點(diǎn)氣動外形優(yōu)化[62] Fig.25 Single-point and multi-point aerodynamic shape optimization of CRM-WBH configuration[62]

Liem等[65]采用Kriging代理模型方法，實(shí)現(xiàn)了28種飛行條件、158個約束和448個氣動/結(jié)構(gòu)設(shè)計變量下，CRM-WBH構(gòu)型的氣動/結(jié)構(gòu)多點(diǎn)綜合優(yōu)化(圖26),以及平均燃油消耗降低8.8%的優(yōu)化設(shè)計目標(biāo)。Gaetan等[66]采用氣動/結(jié)構(gòu)耦合伴隨方法，分別以起飛重量最小和燃油消耗最少為優(yōu)化設(shè)計目標(biāo)，完成了CRM-WBH構(gòu)型476個設(shè)計變量多點(diǎn)優(yōu)化算例(圖27)，分別實(shí)現(xiàn)了起飛重量降低4.2%、燃油消耗降低6.6%和燃油消耗降低11.2%、起飛重量基本不變的設(shè)計目標(biāo)。

圖26 CRM-WBH構(gòu)型氣動/結(jié)構(gòu)分析(Ma=0.85)[65]Fig.26 Aerostructural analysis for CRM-WBH configuration at Ma=0.85[65]

圖27 CRM-WBH構(gòu)型氣動/結(jié)構(gòu)綜合優(yōu)化設(shè)計變量[66]Fig.27 Aerostructural optimization design variables for CRM-WBH configuration [66]

7 CFD驗(yàn)證與確認(rèn)工作的思考與建議

從2001—2016年，DPW會議已經(jīng)連續(xù)舉辦了6屆，其影響力已經(jīng)超出了典型運(yùn)輸機(jī)高速巡航構(gòu)型氣動特性的預(yù)測這一范圍。基于DPW系列會議的成功經(jīng)驗(yàn)，AIAA又發(fā)起了高升力預(yù)測研討會(High Lift Prediction Workshop， HiLiftPW)、氣動彈性預(yù)測研討會(AePW)，與CFD技術(shù)和CRM相關(guān)的其他研究工作也在逐漸展開。本文總結(jié)DPW系列會議的成功經(jīng)驗(yàn)，主要有以下思考和建議：

1) CFD的驗(yàn)證工作進(jìn)展緩慢，需要進(jìn)一步開展CFD驗(yàn)證方法的研究。CFD的驗(yàn)證工作本質(zhì)上是個數(shù)學(xué)問題，主要研究CFD計算模型(方法、軟件)、表達(dá)概念模型(如RANS方程)及其數(shù)值解準(zhǔn)確程度的測度過程。網(wǎng)格收斂性研究是DPW系列會議采用的主要驗(yàn)證手段，Richardson外插方法是獲得網(wǎng)格無關(guān)解的主要方法。這種方法存在的主要問題是：第一，如何構(gòu)造相容的網(wǎng)格序列開展網(wǎng)格收斂性研究？雖然，在典型運(yùn)輸機(jī)構(gòu)型網(wǎng)格生成規(guī)范方面，DPW組委會開展了大量的工作，但如何構(gòu)造相容的網(wǎng)格序列依然是一個懸而未決的問題。第二，Richardson外插方法是否有效？在RANS方程二階空間離散精度的前提下，網(wǎng)格規(guī)模的-2/3次冪僅僅是網(wǎng)格間距的定性度量，只能給出網(wǎng)格間距的一個趨勢；另外，即使按照DPW總結(jié)的網(wǎng)格生成規(guī)范生成不同規(guī)模的數(shù)值模擬網(wǎng)格，許多軟件并不能得到單調(diào)變化的氣動特性收斂曲線，這種情況下靠最密與次密的氣動結(jié)果外插得到的網(wǎng)格無關(guān)解的有效性就值得商榷。

2) 系統(tǒng)地開展CFD的確認(rèn)工作時，不僅所研究的飛行器構(gòu)型應(yīng)該從簡單的部件逐漸過渡到復(fù)雜組合體，研究的物理量也應(yīng)該從宏觀物理量(升力、阻力、力矩……)逐漸過渡到微觀物理量(邊界層速度型、局部表面分離流態(tài)、空間流動結(jié)構(gòu)、摩擦阻力分布、典型站位壓力分布……)的對比。數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果微觀物理量的對比分析，不僅可以進(jìn)一步解釋不同研究手段獲得宏觀量之間的差異(不同研究手段獲得的宏觀物理量的差異是必然存在的)，更可以為計算方法、湍流模型等數(shù)值模擬結(jié)果影響因素的改進(jìn)指明方向。這就要求進(jìn)一步開展精細(xì)的風(fēng)洞測試技術(shù)研究，尤其是非接觸測量技術(shù)研究，如壓敏漆技術(shù)、溫敏漆技術(shù)、粒子圖像測速技術(shù)、轉(zhuǎn)捩測量技術(shù)等，以便為CFD的確認(rèn)工作提供豐富的試驗(yàn)數(shù)據(jù)。

3) 將數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果對比分析是CFD確認(rèn)工作的主要手段，在進(jìn)行數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的對比前，必須先確認(rèn)：第一，數(shù)值模擬所采用的外形與相應(yīng)風(fēng)洞試驗(yàn)的模型是否一致；第二，數(shù)值模擬所采用的邊界條件究竟與風(fēng)洞試驗(yàn)的條件有哪些不同;第三，風(fēng)洞試驗(yàn)結(jié)果的同期或不同期的重復(fù)性精度是多少;第四，風(fēng)洞試驗(yàn)結(jié)果處理的基本方法及誤差。

致 謝

感謝張玉倫、洪俊武、王光學(xué)、張書俊、孟德虹、孫巖、李偉、楊小川等同志堅持不懈的努力，感謝李松同志收集了部分國內(nèi)研究資料，感謝國內(nèi)同行長期以來對TRIP軟件開發(fā)小組的堅定支持。

[1] LEVY D W, VASSBERG J C, WAHLS R A, et al. Summary of data from the First AIAA CFD Drag Prediction Workshop[J]. Journal of Aircraft, 2003, 40(5): 875-882.

[2] LAFLIN K R, VASSBERG J C, WAHLS R A, et al. Summary of data from the Second AIAA CFD Drag Prediction Workshop[J]. Journal of Aircraft, 2005, 42(5):1165-1178.

[3] VASSBERG J C, TINOCO E N, MANI M, et al. Abridged summary of the Third AIAA CFD Drag Prediction Workshop[J]. Journal of Aircraft, 2008, 45(3): 781-798.

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

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

[5] VASSBERG J C, DEHAAN M A, RIVERS S M, et al. Development of a common research model for applied CFD validation studies: AIAA-2008-6919[R]. Reston, VA: AIAA, 2008.

[6] RIVERS M B, DITTBEMER A. Experimental investigation of the NASA common research model: AIAA-2010-4218 [R]. Reston, VA: AIAA, 2010.

[7] RIVERS M B, DITTBEMER A. Experimental investigation of the NASA common research model in the NASA Langley transonic facility and NASA Ames 11-ft transonic wind tunnel: AIAA-2011-1126 [R]. Reston，VA: AIAA, 2011.

[8] RIVERS M B, RUDNIK R, QUEST J. Comparison of the NASA common research model European transonic wind tunnel test data to NASA test data: AIAA-2015-1093 [R]. Reston: AIAA，VA, 2015.

[9] UENO M, KOHZAI T, KOGA S, et al. 80% scaled NASA common research model wind tunnel test of JAXA at relatively low Reynolds number: AIAA-2013-0493[R]. Reston，VA: AIAA, 2013.

[10] CARTIERI A, HUE D, CHANZY Q, et al. Experimental investigations on the common research model at ONERA-S1MA—Comparison with DPW numerical results: AIAA-2017-0964 [R]. Reston，VA: AIAA, 2017.

[11] 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.

[12] VASSBERG J C, TINOCO E N, MANI M, et al. Summary of the Fourth AIAA Computational Fluid Dynamics Drag Prediction Workshop[J]. Journal of Aircraft, 2014, 51(4): 1070-1089.

[13] LEVY D W, LAFLIN K R, TINOCO E N, et al. Summary of data from the Fifth Computational Fluid Dynamics Drag Prediction Workshop: AIAA-2013-0046[R]. Reston，VA: AIAA, 2013.

[14] LEVY D W, LAFLIN K R, TINOCO E N, et al. Summary of data from the Fifth Computational Fluid Dynamics Drag Prediction Workshop [J]. Journal of Aircraft, 2014, 51(4): 1194-1213.

[15] TINOCO E N, BRODERSEN O, KEYE S, et al. Summary of data from the Sixth AIAA CFD Drag Prediction Workshop: CRM case 2 to 5: AIAA-2017-1208[R]. Reston，VA: AIAA, 2017.

[16] KEYE S, MAVRIPLIS D. Summary of data from the Sixth AIAA CFD Drag Prediction Workshop: Case 5 (coupled aero-structural simulation): AIAA-2017-1207[R]. Reston，VA: AIAA, 2017.

[17] RUMSEY C L, MORRISON J H, BIEDRON R T, et al. CFD variability for a civil transport aircraft near buffet-onset conditions: NASA/TM-2003-212149[R]. Washington， D.C.： NASA, 2003.

[18] 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.

[19] VASSBERG J, DEHAAN M, SCLAFANI T. Grid generation requirements for accurate drag predictions based on OVERFLOW calculations： AIAA-2003-4124[R]. Reston, VA: AIAA, 2003.

[20] TINOCO E N, WINKLER C, MANI M, et al. Structured and unstructured solvers for the 3rd CFD Drag Prediction Workshop：AIAA-2007-0255[R]. Reston, VA: AIAA, 2007.

[21] MAVRIPLIS D J. Results from the 3rd Drag Prediction Workshop using the NSU3D unstructured mesh solver：AIAA-2007-0256[R]. Reston, VA: AIAA, 2007.

[22] SCLAFANI A J, VASSBERG J C, HARRISON N A, et al. Drag predictions for the DLR-F6 wing/body and DPW wings using CFL3D and OVERFLOW on an overset mesh： AIAA-2007-0257[R]. Reston, VA: AIAA, 2007.

[23] BRODERSEN O, EISFELD B, RADDATZ J, et al. DLR results from the Third AIAA CFD Drag Prediction Workshop: AIAA-2007-0259[R]. Reston, VA: AIAA, 2007.

[24] SCLAFANI A J, VASSBERG J C, RUMSEY C, et al. Drag prediction for the NASA CRM wing/body/tail using CFL3D and OVERFLOW on an overset Mesh: AIAA-2010-4219 [R]. Reston, VA: AIAA, 2010.

[25] TEMMERMAN L, HIRSCH C. Simulations of the CRM configuration on unstructured hexahedral grids: Lessons learned from the DPW-4 Workshop: AIAA-2010-4670[R]. Reston, VA: AIAA, 2010.

[26] VASSBERG J C. A unified baseline grid about the common research model wing-body for the Fifth AIAA CFD Drag Prediction Workshop: AIAA-2011-3508[R]. Reston, VA: AIAA, 2011.

[27] MARTINEAU D, STOKES S, MUNDAY S, et al. Anisotropic hybrid mesh generation for industrial RANS applications: AIAA-2006-0534[R]. Reston, VA: AIAA, 2006.

[28] DENG X G, ZHANG H X. Developing high-order weighted compact nonlinear schemes[J]. Journal of Computational Physics, 2000, 165: 24-44.

[29] DENG X G, MIN R B, MAO M L, et al. Further studies on geometric conservation law and application to high-order finite difference scheme with stationary grid[J]. Journal of Computational Physics, 2013, 239: 90-111.

[30] 王運(yùn)濤， 孫巖， 王光學(xué), 等. 高階精度方法下的湍流生成項(xiàng)對跨聲速流動數(shù)值模擬的影響研究[J].空氣動力學(xué)學(xué)報，2015, 33(1): 25-30.

WANG Y T, SUN Y, WANG G X, et al. Numerical study of the effect of turbulent production terms on the simulation of transonic flows with high-order numerical method[J]. Acta Aerodynamica Sinica, 2015, 33(1): 25-30 (in Chinese).

[31] 王運(yùn)濤， 孫巖， 王光學(xué), 等. 湍流模型離散精度對數(shù)值模擬影響的計算分析[J]. 航空學(xué)報, 2015, 36(5): 1453-1459.

WANG Y T, SUN Y, WANG G X, et al. Numerical study of the effect of turbulent production terms on the simulation of transonic flows with high-order numerical method[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(5): 1453-1459 (in Chinese)

[32] 王運(yùn)濤， 孫巖， 李松, 等. 高階精度方法下的湍流生成項(xiàng)對低速流動數(shù)值模擬的影響研究[J]. 空氣動力學(xué)學(xué)報, 2015,33(3): 325-329.

WANG Y T, SUN Y, LI S, et al. Numerical analysis of the effect of turbulent production terms in low-speed numerical simulation[J]. Acta Aerodynamica Sinica, 2015,33(3): 325-329 (in Chinese)

[33] 王運(yùn)濤， 孫巖， 王光學(xué)， 等. DLR-F6翼身組合體的高階精度數(shù)值模擬[J]. 航空學(xué)報， 2015, 36(9): 2923-2929.

WANG Y T, SUN Y, WANG G X, et al. High-order numerical simulation of DLR-F6 wing-body configuration[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(9): 2923-2929 (in Chinese).

[34] 王光學(xué)， 張玉倫， 王運(yùn)濤, 等. BLU-SGS方法在WCNS高階精度格式上的數(shù)值分析[J]. 空氣動力學(xué)學(xué)報, 2015, 33(6):733-739.

WANG G X, ZHANG Y L, WANG Y T, et al. Numerical analysis of BLU_SGS method in WCNS high-order scheme[J]. Acta Aerodynamica Sinica, 2015, 33(6): 733-739 (in Chinese).

[35] 王運(yùn)濤， 孟德虹， 孫巖, 等. 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 numerical simulation of DLR-F6/FX2B wing-body configuration[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(2): 484-490 (in Chinese).

[36] 王運(yùn)濤， 孫巖， 孟德虹, 等. CRM翼/身/平尾組合體模型高階精度數(shù)值模擬[J]. 航空學(xué)報， 2016， 37(12)： 3692-3697.

WANG Y T, SUN Y, MENG D H, et al. High-order numerical simulation of CRM wing/body/horizontal tail model[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(12): 3692-3697 (in Chinese).

[37] 王運(yùn)濤， 孫巖， 孟德虹, 等. CRM翼身組合體模型高階精度數(shù)值模擬[J]. 航空學(xué)報， 2017, 38(3): 120298.

WANG Y T, SUN Y, MENG D H, et al. High-order numerical simulation of CRM wing-body model[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(3): 120298 (in Chinese).

[38] 王運(yùn)濤， 孟德虹， 孫巖, 等. CRM-WB風(fēng)洞模型高階精度數(shù)值模擬[J]. 航空學(xué)報， 2018, 39(4): 121642.

WANG Y T, MENG D H, SUN Y, et al. High-order numerical simulation of CRM-WB wind tunnel model[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(4): 121642 (in Chinese).

[39] ROY C J, RUMSEY C, TINOCO E N. Summary of data from the Sixth AIAA CFD Drag Prediction Workshop: Case 1 code verification[R]. Reston，VA: AIAA, 2017.

[40] SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston，VA: AIAA, 1992.

[41] MENTER F R. Two-equation eddy-viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[42] ANTON P S, JOHNSON D J, BLOCK M, et al. Wind tunnle and propulsion test facilities: Supporting analyses toan assessment of NASA’s capabilities to serve national needs: RAND_TR134[R]. California: RAND, 2004.

[43] SPALART P R. Strategies for turbulence modelling and simulations[J]. International Journal of Heat and Fluid Flow, 2000, 21: 252-263.

[44] YAMAMOTO K, TANAKA K, MURAYAMA M. Effect of a nonlinear constitutive relation for turbulence modeling on predicting flow separation at wing-body juncture of transonic commercial aircraft: AIAA-2012-2895[R]. Reston，VA: AIAA, 2012.

[45] OBERKAMPF W L, TRUCANOB T G. Verification and validation in computational fluid dynamics[J]. Progress in Aerospace Sciences, 2002, 38: 209-272.

[46] DAVID H U E. CFD investigation on the DPW-5 configuration with measured experimental wing twist using the elsA slover and the far-field approach：AIAA-2013-2508[R]. Reston，VA： AIAA, 2013.

[47] KEYE S, TOJITI V, EISFELD B, et al. Investigation of fluid-structure-coupling and turbulence model effects on the DLR results of the Fifth AIAA CFD Drag Prediction Workshop：AIAA-2013-2509[R]. Reston，VA：AIAA, 2013.

[48] KEYE S, BRODERSEN O, RIVERS M B, et al. Investigation of aeroelastic effects on the NASA common research model[J]. Journal of Aircraft, 2014, 51(4):1323-1330.

[49] MORRISON J H. Statistical analysis of the fifth drag prediction workshop computational fluid dynamics solutions[J]. Journal of Aircraft, 2014, 51(4):1214-1222.

[50] DERLAGA J M, MORRISONY J H. Statistical analysis of CFD solutions from the 6th AIAA CFD drag prediction workshop: AIAA-2017-1209[R]. Reston，VA: AIAA, 2017.

[51] RIVERS M B, HUNTER C A. Support system effects on the NASA common research model：AIAA-2012-0707[R]. Reston，VA： AIAA, 2012.

[52] RIVERS M B, HUNTER C A, CAMPBELL R L. Further investigation of the support system effects and wing twist on the NASA common research model：AIAA-2012-3209[R]. Reston，VA： AIAA, 2012.

[53] ZILLIAC G G, PULLIAM T H, RIVERS M B, et al. A comparison of the measured and computed skin friction distribution on the common research model：AIAA-2011-1129[R]. Reston，VA： AIAA, 2011.

[54] 王運(yùn)濤， 孫巖， 孟德虹， 等. 包含支撐裝置和機(jī)翼變形的CRM-WB模型氣動特性數(shù)值模擬[J]. 航空學(xué)報， 2017， 38(12)： 121202.

WANG Y T, SUN Y, MENG D H, et al. Numerical simulation of aerodynamic characteristics of CRM-WB model with support system and wing deformation[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(12): 121202 (in Chinese).

[55] 王運(yùn)濤, 王光學(xué), 張玉倫. 采用TRIP2.0軟件計算DLR-F6構(gòu)型的阻力[J]. 空氣動力學(xué)學(xué)報， 2009, 27(1): 108-111.

WANG Y T, WANG G X, ZHANG Y L. Drag prediction of DLR-F6 configuration with TRIP2.0 software[J]. Acta Aerodynamica Sinica, 2009, 27(1): 108-111 (in Chinese).

[56] 王運(yùn)濤， 王光學(xué)， 張玉倫. DPW Ⅲ機(jī)翼和翼身組合體構(gòu)型數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報， 2011, 29(3): 264-269.

WANG Y T, WANG G X, ZHANG Y L. Numerical simulation of DPW Ⅲ wing and wing-body configurations[J]. Acta Aerodynamica Sinica, 2011, 29(3): 264-269 (in Chinese).

[57] OWENS L R, WAHLS R A, RIVERS S M. Off-design Reynolds number effects for a supersonic transport[J]. Journal of Aircraft, 2005, 42(6): 1427-1441.

[58] KEYE S, RUDNIK R. Aero-elastic simulation of DLR’s F6 transport aircraft configuration and comparison to experimental data: AIAA-2009-0580[R]. Reston, VA：AIAA, 2009.

[59] MOUTON S, SANT Y L, LYONNET M. Prediction of the aerodynamic effect of model deformation during transonic wind tunnel tests[J]. International Journal of Engineering Systems Modelling and Simulation, 2013, 5(1-3): 44-56.

[60] HEEG J, CHWALOWSKI P, FLORANCE J P, et al. Overview of the aeroelastic predication workshop: AIAA-2013-0783[R]. Reston, VA： AIAA, 2013.

[61] LYU Z J, KENWAY G K W, MARTINS J R R A. RANS-based aerodynamic shape optimization investigations of the common research model wing: AIAA-2014-0567[R]. Reston, VA: AIAA, 2014.

[62] KENWAY G K W, MARTINS J R R A. Aerodynamic shape optimization of the CRM configuration including buffet-onset conditions: AIAA-2016-1294[R]. Reston, VA: AIAA, 2016.

[63] 陳頌, 白俊強(qiáng), 史亞云, 等. 民用客機(jī)機(jī)翼/機(jī)身/平尾構(gòu)型氣動外形優(yōu)化設(shè)計[J]. 航空學(xué)報, 2015, 36(10): 3195-3207.

CHEN S, BAI J Q, SHI Y Y, et al. Aerodynamic shape optimization design of civil jet wing-body-tail configuration[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3195-3207 (in Chinese).

[64] 劉峰博, 郝海兵, 李典, 等. 離散伴隨方法在氣動優(yōu)化設(shè)計中的應(yīng)用[J]. 航空計算技術(shù), 2017, 42(2): 33-40.

LIU F B, HAO H B, LI D, et al. Application of discrete adjoint method in aerodynamic shape optimization design[J]. Aeronautical Computing Technique, 2017, 42(2): 33-40 (in Chinese).

[65] LIEM R P, KENWAY G K W, MARTINS J R R A. Multi-point, multi-mission, high-fidelity aerostructural optimization of a long-range aircraft configuration: AIAA-2012-5706[R]. Reston, VA: AIAA, 2012.

[66] GAETAN K W, KENWAY G K W, MARTINS J R R A. Multipoint high-fidelity aerostructural optimization of a transport aircraft configuration[J]. Journal of Aircraft, 2014, 51(1): 144-160.