CRM-WB風(fēng)洞模型高階精度數(shù)值模擬

王運(yùn)濤，孟德虹，孫巖，李偉

1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力學(xué)研究所， 綿陽(yáng) 621000 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 綿陽(yáng) 621000

AIAA阻力預(yù)測(cè)研討會(huì)[1-4](Drag Prediction Workshop， DPW)吸引了世界范圍內(nèi)空氣動(dòng)力學(xué)研究工作者的廣泛參與，已成為CFD驗(yàn)證與確認(rèn)研究領(lǐng)域最重要的國(guó)際合作之一。2012年6月召開的第5屆DPW(DPW Ⅴ)采用了CRM (Common Research Model)[5]翼身組合體構(gòu)型作為基準(zhǔn)研究模型，并在NASA Langley的NTF(National Transonic Facility)風(fēng)洞、NASA Ames的TWT(Transonic Wind Tunnel)中開展了相應(yīng)的風(fēng)洞試驗(yàn)[6-7]。DPW Ⅴ上，基于雷諾平均Navier-Stokes(RANS)方程和2階空間離散精度的計(jì)算方法，來(lái)自世界各地的22家研究機(jī)構(gòu)共提供了57組CRM翼身組合體構(gòu)型計(jì)算結(jié)果[8]，對(duì)這些計(jì)算結(jié)果的統(tǒng)計(jì)分析表明，設(shè)計(jì)升力系數(shù)下，來(lái)流迎角的計(jì)算結(jié)果普遍低于相應(yīng)的試驗(yàn)結(jié)果，低頭力矩的計(jì)算結(jié)果普遍大于相應(yīng)的試驗(yàn)結(jié)果；相同迎角下，升力系數(shù)、低頭力矩的計(jì)算結(jié)果普遍大于相應(yīng)的試驗(yàn)結(jié)果。本文作者在前期研究工作[9]中采用5階空間離散精度的WCNS(Weighted Compact Nonlinear Scheme)開展了CRM翼身組合體構(gòu)型的數(shù)值模擬，計(jì)算結(jié)果與試驗(yàn)結(jié)果之間同樣存在較大的差異。

Levy等在文獻(xiàn)[8]中總結(jié)了導(dǎo)致CRM翼身組合體構(gòu)型氣動(dòng)特性風(fēng)洞試驗(yàn)與計(jì)算結(jié)果之間差異的主要因素包括：數(shù)值計(jì)算模型中沒有包括試驗(yàn)?zāi)Ｐ驮跉鈩?dòng)載荷作用下的靜氣動(dòng)彈性變形，數(shù)值計(jì)算模型沒有包括風(fēng)洞模型尾部支撐裝置，數(shù)值模擬中沒有考慮風(fēng)洞試驗(yàn)中固定轉(zhuǎn)捩位置等。基于風(fēng)洞試驗(yàn)中測(cè)量得到的風(fēng)洞模型靜氣動(dòng)彈性變形數(shù)據(jù)，David[10]采用二階精度計(jì)算方法和結(jié)構(gòu)網(wǎng)格技術(shù)研究了CRM翼身組合體構(gòu)型靜氣動(dòng)彈性變形對(duì)數(shù)值模擬結(jié)果的影響；Keye等[11]采用2階精度格式、非結(jié)構(gòu)網(wǎng)格技術(shù)和流固耦合方法研究了靜氣動(dòng)彈性變形對(duì)CRM翼身組合體構(gòu)型數(shù)值模擬結(jié)果的影響；但上述研究工作所采用的計(jì)算模型中并沒有包含模型支撐裝置。2016年6月召開的第6屆DPW(DPW Ⅵ)選擇了包括靜氣動(dòng)彈性變形的CRM翼身組合體構(gòu)型作為基準(zhǔn)研究模型，以期提高數(shù)值模擬結(jié)果與風(fēng)洞試驗(yàn)結(jié)果的吻合程度，來(lái)自世界各地的25家研究機(jī)構(gòu)共提供了54組計(jì)算結(jié)果[12]。DPW Ⅵ數(shù)值模擬結(jié)果的統(tǒng)計(jì)分析表明，考慮了CRM翼身組合體構(gòu)型的靜氣動(dòng)變形后，氣動(dòng)特性計(jì)算結(jié)果與試驗(yàn)結(jié)果的吻合程度有顯著改善，但升力系數(shù)、低頭力矩系數(shù)的計(jì)算結(jié)果普遍高于試驗(yàn)結(jié)果的現(xiàn)象依然存在，計(jì)算方法基本采用2階空間離散精度的差分格式依然沒有包含模型支撐裝置的影響。

在文獻(xiàn)[9]網(wǎng)格收斂性研究工作的基礎(chǔ)上，基于DPW Ⅵ組委會(huì)提供的CRM翼身組合體數(shù)值模擬和NASA NTF風(fēng)洞的支撐裝置數(shù)值模擬，本文采用5階空間離散精度的WCNS和多塊對(duì)接結(jié)構(gòu)網(wǎng)格技術(shù)，開展了CRM翼身組合體風(fēng)洞試驗(yàn)?zāi)Ｐ偷母唠A精度數(shù)值模擬，主要目的是采用高階精度方法，評(píng)估風(fēng)洞模型靜氣動(dòng)彈性變形和支撐裝置對(duì)CRM翼身組合體構(gòu)型氣動(dòng)特性的影響。

1 CRM翼身組合體風(fēng)洞模型

CRM構(gòu)型是由NASA和DPW組織委員會(huì)聯(lián)合設(shè)計(jì)開發(fā)的寬體運(yùn)輸機(jī)構(gòu)型，主要目的是為CFD的驗(yàn)證和確認(rèn)工作提供基準(zhǔn)外形，設(shè)計(jì)馬赫數(shù)為0.85、升力系數(shù)為0.50。CRM系列構(gòu)型包括了翼身組合體、翼/身/平尾組合體和翼/身/平尾/掛架/吊艙組合體等不同構(gòu)型，DPW Ⅴ組委會(huì)和DPW Ⅵ組委會(huì)均選擇了CRM翼身組合體構(gòu)型做為基準(zhǔn)研究模型，不同的是DPW Ⅵ組委會(huì)選擇的CRM翼身組合體構(gòu)型包含了不同迎角下風(fēng)洞試驗(yàn)測(cè)量得到的靜氣動(dòng)彈性變形。

NASA NTF風(fēng)洞試驗(yàn)中，CRM翼身組合體試驗(yàn)?zāi)Ｐ涂s比為0.027，基本參數(shù)為：模型參考面積Sref=0.279 7 m2;平均氣動(dòng)弦長(zhǎng)c=1.891 m;展長(zhǎng)b=1.587 m;梢根比λ=0.275;展弦比AR=9.0;1/4弦線后掠角Λc/4= 35.0°;馬赫數(shù)Ma=0.85;雷諾數(shù)Re=5.0×106。圖1給出了DPW Ⅵ組委會(huì)提供的包含不同來(lái)流迎角α下靜氣動(dòng)彈性變形的CRM翼身組合體計(jì)算模型，同時(shí)給出了CRM翼身組合體“剛性”外形的計(jì)算模型(圖1中紅色部分)。在固定迎角下，機(jī)翼的靜氣動(dòng)彈性變形沿機(jī)翼展向逐漸增加；隨著來(lái)流迎角的增加，機(jī)翼的靜氣動(dòng)彈性變形逐漸增加；迎角α=3.0°時(shí)，翼梢處的彎曲變形為17.4 mm，扭轉(zhuǎn)變形達(dá)到-1.1°。為方便討論，將CRM翼身組合體“剛性”外形計(jì)算模型標(biāo)識(shí)為CRM-WB，將包含靜氣動(dòng)彈性變形的CRM翼身組合體外形計(jì)算模型標(biāo)識(shí)為CRM-WB-A。

風(fēng)洞模型采用安裝于機(jī)身后體的葉片尾撐方式固定于風(fēng)洞迎角變換裝置(如圖2所示)。由于風(fēng)洞試驗(yàn)?zāi)Ｐ秃蟛康挠亲儞Q裝置對(duì)氣動(dòng)特性的影響屬于二次影響[13]，本文的數(shù)值模擬沒有考慮模型支撐后部的迎角變換裝置；對(duì)模型葉片尾撐延伸段進(jìn)行了局部修型處理以避免底部分離導(dǎo)致的計(jì)算收斂困難(圖2中紅色部分)。本文將包含支撐裝置和靜氣動(dòng)彈性變形的CRM翼身組合體計(jì)算模型標(biāo)識(shí)為CRM-WBS-A。

圖1 CRM-WB-A構(gòu)型的計(jì)算模型Fig.1 Computational model of configuration ofCRM-WB-A

圖2 CRM-WBS-A構(gòu)型的計(jì)算模型Fig.2 Computational model of configuration ofCRM-WBS-A

2 計(jì)算網(wǎng)格與高階精度計(jì)算方法

根據(jù)DPW組織委員會(huì)給出的網(wǎng)格生成指導(dǎo)原則，文獻(xiàn)[9]采用粗、中、細(xì)和極細(xì)4套多塊對(duì)接結(jié)構(gòu)開展了CRM-WB模型的網(wǎng)格收斂性研究，并獲得了具有網(wǎng)格收斂性的氣動(dòng)特性高階精度計(jì)算結(jié)果。采用文獻(xiàn)[9]中CRM-WB構(gòu)型的中等網(wǎng)格作為基準(zhǔn)網(wǎng)格開展本文的研究工作。CRM-WB構(gòu)型半模網(wǎng)格規(guī)模為8 440 223，物面第1層法向無(wú)量綱距離為0.94。利用DPW Ⅵ組委會(huì)提供的不同迎角下的靜氣動(dòng)彈性變形(如圖1所示)，采用與基準(zhǔn)網(wǎng)格相同的網(wǎng)格拓?fù)浼熬W(wǎng)格分布，通過(guò)基于徑向基[14]與超限插值[15]的復(fù)合型動(dòng)態(tài)網(wǎng)格變形方法RBF-TFI[16]，構(gòu)造了不同來(lái)流迎角下的CRM-WB-A構(gòu)型的計(jì)算網(wǎng)格。在此基礎(chǔ)上，進(jìn)一步構(gòu)造了不同來(lái)流迎角下CRM-WBS-A構(gòu)型對(duì)稱面計(jì)算網(wǎng)格如圖3所示，半模網(wǎng)格規(guī)模達(dá)到了16 133 343。

本文采用有限差分方法離散任意坐標(biāo)系下的RANS方程組，控制方程的對(duì)流項(xiàng)離散采用5階精度的WCNS，黏性項(xiàng)的離散采用6階精度中心格式，邊界及近邊界條件采用單邊4階精度離散，以上方法的詳細(xì)介紹見文獻(xiàn)[17-18]；湍流模型采用Menter 剪切應(yīng)力輸運(yùn)(SST)兩方程湍流模型[19]，數(shù)值模擬采用“全湍流”方式；離散方程組的求解采用BLU-SGS方法[20-21]。

圖3 CRM-WBS-A構(gòu)型對(duì)稱面計(jì)算網(wǎng)格Fig.3 Computation grid on symmetric plan of CRM-WBS-A configuration

3 計(jì)算結(jié)果及討論

采用CRM-WB、CRM-WB-A和CRM-WBS-A共3種計(jì)算模型和高階精度計(jì)算方法，開展設(shè)計(jì)馬赫數(shù)下靜氣動(dòng)彈性變形和模型支撐裝置對(duì)數(shù)值模擬結(jié)果的影響研究，對(duì)比數(shù)據(jù)為NASA 2.5 m×2.5 m NTF風(fēng)洞試驗(yàn)測(cè)力和測(cè)壓數(shù)據(jù)[4]。數(shù)值模擬來(lái)流條件為：馬赫數(shù)Ma=0.85;基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Re=5.0×106;迎角α=2.5°～4.0°，迎角間隔Δα=0.25°。

3.1 氣動(dòng)力系數(shù)隨迎角的變化

圖4給出了采用高階精度計(jì)算方法得到的CRM翼身組合體3種計(jì)算模型的升力系數(shù)CL、阻力系數(shù)CD和俯仰力矩系數(shù)Cm隨來(lái)流迎角的變化，同時(shí)給出了NASA NTF風(fēng)洞的測(cè)力試驗(yàn)結(jié)果。

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

在計(jì)算迎角范圍內(nèi)，CRM-WB構(gòu)型與CRM-WB-A構(gòu)型氣動(dòng)特性的計(jì)算結(jié)果相比較，靜氣動(dòng)彈性變形的影響使得相同迎角下的升力系數(shù)、阻力系數(shù)和低頭力矩系數(shù)-Cm均下降，數(shù)值模擬得到的2種構(gòu)型的失速迎角均在3.75°左右。CRM-WB-A構(gòu)型與CRM-WBS-A構(gòu)型氣動(dòng)特性的計(jì)算結(jié)果相比較，進(jìn)一步考慮了模型尾部支撐后，相同迎角下的升力系數(shù)、阻力系數(shù)和低頭力矩系數(shù)進(jìn)一步下降，失速迎角后移。3種構(gòu)型的氣動(dòng)特性計(jì)算結(jié)果相比較，失速迎角以前靜氣動(dòng)變形對(duì)升力系數(shù)和俯仰力矩系數(shù)的影響量略高于模型支撐的影響量，而靜氣動(dòng)彈性變形和模型支撐對(duì)阻力系數(shù)的影響量基本相當(dāng)；CRM-WBS-A構(gòu)型的氣動(dòng)特性計(jì)算結(jié)果最接近風(fēng)洞試驗(yàn)測(cè)力結(jié)果。計(jì)算結(jié)果和試驗(yàn)結(jié)果仍存在偏差的原因主要是2個(gè)方面：① 風(fēng)洞試驗(yàn)應(yīng)給出試驗(yàn)結(jié)果的誤差帶并對(duì)風(fēng)洞試驗(yàn)結(jié)果的各種修正方式作進(jìn)一步的研究。② 風(fēng)洞試驗(yàn)在機(jī)翼、機(jī)頭粘貼了轉(zhuǎn)捩帶，而本文的數(shù)值模擬采用的是全湍流方式。

3.2 機(jī)翼表面壓力系數(shù)分布 (α=3.0°)

圖5給出了迎角α=3.0°時(shí)，采用高階精度計(jì)算方法得到的CRM翼身組合體構(gòu)型3個(gè)典型展向位置機(jī)翼剖面的壓力系數(shù)Cp分布曲線，同時(shí)給出了NTF風(fēng)洞相鄰迎角的測(cè)壓數(shù)據(jù)[4]，圖5中橫坐標(biāo)x為機(jī)翼流向無(wú)量綱距離，η為機(jī)翼展向無(wú)量綱距離。

從圖5可以看出，在內(nèi)側(cè)機(jī)翼(η=0.131、0.397)的位置，靜氣動(dòng)彈性變形使得機(jī)翼上翼面激波位置略微前移，對(duì)其他位置的壓力分布基本沒有影響；而模型支撐裝置使得機(jī)翼上翼面的激波位置明顯前移。在η=0.727站位，靜氣動(dòng)彈性變形使得機(jī)翼上翼面激波位置前移、激波位置以前的負(fù)壓明顯下降， 機(jī)翼下翼面50%弦長(zhǎng)前的壓力略有增加；模型支撐裝置使得上翼面的激波位置顯著前移。在機(jī)翼梢部(η=0.950)的位置，靜氣動(dòng)彈性變形對(duì)壓力分布的影響與η=0.727站位類似；而模型支撐裝置使得上翼面流動(dòng)結(jié)構(gòu)發(fā)生顯著變化，機(jī)翼上翼面出現(xiàn)了明顯的雙激波結(jié)構(gòu)。總之，由于靜氣動(dòng)變形導(dǎo)致機(jī)翼負(fù)扭轉(zhuǎn)角由翼根到翼梢逐漸增加，靜氣動(dòng)彈性變形對(duì)機(jī)翼表面壓力系數(shù)分布的影響也由翼根到翼梢逐漸增加，主要使得機(jī)翼上表面激波位置前移和外側(cè)激波位置以前的負(fù)壓值下降；支撐裝置引起機(jī)翼上表面激波位置前移(翼梢處除外)，并在機(jī)翼翼梢處導(dǎo)致機(jī)翼上表面激波結(jié)構(gòu)發(fā)生變化。

圖5 CRM翼身組合體構(gòu)型典型站位壓力系數(shù)分布曲線 Fig.5 Curves of pressure coefficients distribution at different spanwise locations of CRM wing-body configuration

3.3 機(jī)翼上表面流線 (α=4.0°)

從圖4(a)中看出，迎角α=4.0°時(shí)，與CRM-WBS-A構(gòu)型的計(jì)算結(jié)果不同，采用CRM-WB構(gòu)型與CRM-WB-A構(gòu)型得到升力系數(shù)均出現(xiàn)明顯下降。圖6給出了迎角α=4.0°時(shí)，采用高階精度計(jì)算方法得到的CRM翼身組合體3種構(gòu)型的上表面流線，機(jī)翼用壓力分布著色。由圖6可見，迎角α=4.0°時(shí)，CRM-WB構(gòu)型與CRM-WB-A構(gòu)型升力系數(shù)下降是由于翼身結(jié)合部后緣分離區(qū)突然增加導(dǎo)致的，而CRM-WBS-A構(gòu)型在翼身結(jié)合部后緣則沒有明顯的分離區(qū)。

圖6 CRM翼身組合體構(gòu)型機(jī)翼上表面流線(α=4.0°)Fig.6 Streamlines on wing upper surface of CRM wing-body configuration (α=4.0°)

4 結(jié) 論

基于5階空間離散精度的WCNS和SST兩方程湍流模型，開展了靜氣動(dòng)變形和模型支撐裝置對(duì)CRM翼身組合體構(gòu)型氣動(dòng)特性的影響，通過(guò)與CRM-WB計(jì)算模型和NASA NTF風(fēng)洞試驗(yàn)結(jié)果的對(duì)比分析，基本結(jié)論如下：

1)α=3.0°時(shí)，模型的靜氣動(dòng)彈性變形主要導(dǎo)致機(jī)翼上表面的激波位置前移及外側(cè)機(jī)翼激波位置以前的負(fù)壓下降。

2)α=3.0°時(shí)，模型支撐裝置主要導(dǎo)致機(jī)翼上表面激波位置前移(翼梢除外)，并使得翼梢位置上表面流動(dòng)結(jié)構(gòu)發(fā)生明顯變化。

3)α=4.0°時(shí)，計(jì)算模型中沒有包含支撐裝置是導(dǎo)致翼身結(jié)合部后緣出現(xiàn)局部分離進(jìn)而導(dǎo)致升力系數(shù)下降的主要原因。

4) 計(jì)算模型中同時(shí)包含靜氣動(dòng)彈性變形和模型支撐裝置顯著改善了數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的吻合程度。

致 謝

感謝張玉倫、洪俊武、張書俊和楊小川等同志在高階精度格式程序?qū)崿F(xiàn)方面所作的研究工作，感謝中國(guó)航空研究院白文博士在數(shù)據(jù)分析方面提供的幫助。

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

[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, DITTBERNER A. Experimental investigation of the NASA common research model (invited): 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 (invited): AIAA-2011-1126[R]. Reston, VA: AIAA, 2011.

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

[9] 王運(yùn)濤, 孫巖, 孟德虹, 等. CRM翼身組合體模型高階精度數(shù)值模擬[J]. 航空學(xué)報(bào), 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).

[10] DAVID H. 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.

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

[12] TINOCO E N, BRODERSEN O P, 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.

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

[14] RENDALL T C S, ALLEN C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2009, 228(17): 6231-6249.

[15] SONI B K. Grid generation for internal flow configuration[J]. Computers & Mathematics with Applications, 1992, 24(5-6): 191-201.

[16] 孫巖, 鄧小剛, 王運(yùn)濤, 等. RBF_TFI結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)在風(fēng)洞靜氣動(dòng)彈性修正中的應(yīng)用[J]. 工程力學(xué), 2014, 31(10): 228-233.

SUN Y, DENG X G, WANG Y T, et al. Application of structural dynamic grid method based on RBF_TFI on wind tunnel static aero-elastic modification[J]. Engineering Mechanics, 2014, 31(10): 228-233 (in Chinese).

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

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

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

[20] CHEN R F, WANG Z J. Fast, Block lower-upper symmetric Gauss-Seidel scheme for arbitrary grids[J]. AIAA Journal, 2000, 38(12): 2238-2245.

[21] 王光學(xué), 張玉倫, 王運(yùn)濤, 等. BLU-SGS方法在WCNS高階精度格式上的數(shù)值分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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).