王運濤，孟德虹，孫巖，李偉

1.中國空氣動力研究與發(fā)展中心 計算空氣動力學研究所， 綿陽 621000 2.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室， 綿陽 621000

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

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

在文獻[9]網(wǎng)格收斂性研究工作的基礎上，基于DPW Ⅵ組委會提供的CRM翼身組合體數(shù)值模擬和NASA NTF風洞的支撐裝置數(shù)值模擬，本文采用5階空間離散精度的WCNS和多塊對接結構網(wǎng)格技術，開展了CRM翼身組合體風洞試驗模型的高階精度數(shù)值模擬，主要目的是采用高階精度方法，評估風洞模型靜氣動彈性變形和支撐裝置對CRM翼身組合體構型氣動特性的影響。

1 CRM翼身組合體風洞模型

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

NASA NTF風洞試驗中，CRM翼身組合體試驗模型縮比為0.027，基本參數(shù)為：模型參考面積Sref=0.279 7 m2;平均氣動弦長c=1.891 m;展長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 Ⅵ組委會提供的包含不同來流迎角α下靜氣動彈性變形的CRM翼身組合體計算模型，同時給出了CRM翼身組合體“剛性”外形的計算模型(圖1中紅色部分)。在固定迎角下，機翼的靜氣動彈性變形沿機翼展向逐漸增加；隨著來流迎角的增加，機翼的靜氣動彈性變形逐漸增加；迎角α=3.0°時，翼梢處的彎曲變形為17.4 mm，扭轉變形達到-1.1°。為方便討論，將CRM翼身組合體“剛性”外形計算模型標識為CRM-WB，將包含靜氣動彈性變形的CRM翼身組合體外形計算模型標識為CRM-WB-A。

風洞模型采用安裝于機身后體的葉片尾撐方式固定于風洞迎角變換裝置(如圖2所示)。由于風洞試驗模型后部的迎角變換裝置對氣動特性的影響屬于二次影響[13]，本文的數(shù)值模擬沒有考慮模型支撐后部的迎角變換裝置；對模型葉片尾撐延伸段進行了局部修型處理以避免底部分離導致的計算收斂困難(圖2中紅色部分)。本文將包含支撐裝置和靜氣動彈性變形的CRM翼身組合體計算模型標識為CRM-WBS-A。

圖1 CRM-WB-A構型的計算模型Fig.1 Computational model of configuration ofCRM-WB-A

圖2 CRM-WBS-A構型的計算模型Fig.2 Computational model of configuration ofCRM-WBS-A

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

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

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

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

3 計算結果及討論

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

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

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

圖4 CRM翼身組合體構型氣動特性Fig.4 Aerodynamic characteristics of CRM wing-body configuration

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

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

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

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

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

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

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

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

4 結 論

基于5階空間離散精度的WCNS和SST兩方程湍流模型，開展了靜氣動變形和模型支撐裝置對CRM翼身組合體構型氣動特性的影響，通過與CRM-WB計算模型和NASA NTF風洞試驗結果的對比分析，基本結論如下：

1)α=3.0°時，模型的靜氣動彈性變形主要導致機翼上表面的激波位置前移及外側機翼激波位置以前的負壓下降。

2)α=3.0°時，模型支撐裝置主要導致機翼上表面激波位置前移(翼梢除外)，并使得翼梢位置上表面流動結構發(fā)生明顯變化。

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

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

致 謝

感謝張玉倫、洪俊武、張書俊和楊小川等同志在高階精度格式程序實現(xiàn)方面所作的研究工作，感謝中國航空研究院白文博士在數(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] 王運濤, 孫巖, 孟德虹, 等. CRM翼身組合體模型高階精度數(shù)值模擬[J]. 航空學報, 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] 孫巖, 鄧小剛, 王運濤, 等. RBF_TFI結構動網(wǎng)格技術在風洞靜氣動彈性修正中的應用[J]. 工程力學, 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] 王光學, 張玉倫, 王運濤, 等. BLU-SGS方法在WCNS高階精度格式上的數(shù)值分析[J]. 空氣動力學學報, 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).