CRM翼身組合體模型高階精度數值模擬

王運濤， 孫巖， 孟德虹,*， 王光學

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

CRM翼身組合體模型高階精度數值模擬

王運濤1， 孫巖2， 孟德虹1,*， 王光學1

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

基于五階空間離散精度的WCNS格式，開展了CRM翼身組合體模型的高階精度數值模擬，以評估WCNS格式對復雜外形的模擬能力以及典型運輸機巡航構型阻力預測的精度。首先依照DPW組委會提出的網格生成指導原則，利用ICEM軟件生成了粗、中、細、極細四套網格，網格規模從“粗網格”的2 578 687個網格點逐漸擴展到“極細網格”的65 464 511 個網格點。研究了設計升力系數下，網格規模對氣動特性、壓力分布和翼根后緣局部分離區的影響，采用“中等網格”開展了抖振特性的數值模擬研究。通過與二階精度的計算結果、DPW V統計結果和部分試驗結果的對比分析，高階精度數值模擬結果表明，阻力系數計算結果與DPW V統計平均結果吻合較好；網格密度對機翼上表面的激波位置和翼身結合部后緣局部分離區略有影響；迎角為4°時，升力系數下降的主要原因是機翼上表面激波誘導分離區和翼身結合部后緣局部分離區的增加。

RANS方程； WCNS格式； CRM模型； 流場模擬； 網格密度； 氣動特性

AIAA阻力預測會議DPW(Drag Prediction Workshop) 從2001年發起到現在，已經成功舉辦了5屆[1-5]并持續了十多年時間。DPW系列會議的宗旨是評估基于雷諾平均Navier-Stokes(RANS)方程的各種CFD(Computational Fluid Dynamics)方法在典型運輸機構型氣動特性預測尤其是阻力預測方面的現狀，明確CFD技術的發展方向，并逐步建立一個評估CFD可信度的國際交流平臺。通過提供標準研究模型、發布基準網格并公開試驗數據，DPW系列會議的影響日益擴大，獲得了世界范圍內相關研究機構的廣泛參與，積累了豐富的計算數據和試驗數據，已經成為CFD驗證與確認發展歷程中最重要的國際合作之一。

第5屆DPW(DPW V)會議于2012年6月在美國路易安娜州的新奧爾良市召開，這次會議采用了與DPW IV相同的CRM(Common Research Model)模型[6]，不同的是，DPW V的研究構型去掉了CRM模型的平尾，只包含了機身和機翼，簡稱為CRM-WB，計算狀態包括了網格收斂性研究和抖振特性研究兩個方面。來自世界各地的22家研究機構共提供了57組計算結果[5]。這些基于RANS方程的計算結果基本上采用了二階空間離散精度的計算方法，采用三階離散精度以上差分格式的數值模擬結果尚未見公開報道。

高階精度格式一直是CFD領域的研究熱點，但在復雜外形上的應用才剛剛起步[7]。鄧小剛和張涵信[8]提出的WCNS(Weighted Compact Nonlinear Scheme)具有五階空間離散精度，通過在幾何守恒律方面持續不斷的研究工作[9]，已經成功應用于典型運輸機構型的高階精度數值模擬[10-12]，并取得了良好效果。

本文采用五階空間離散精度的WCNS格式對CRM翼身組合體(CRM-WB)模型開展了高階精度數值模擬，依據DPW V約定的計算狀態，開展了固定升力系數下的網格收斂性研究和固定馬赫數下的抖振特性研究，通過與二階精度計算結果、DPW V的統計結果和相應試驗結果的對比分析，進一步確認了WCNS格式模擬典型運輸機構型的能力。

1 CRM-WB模型與高階精度計算方法

CRM模型由NASA的亞聲速固定翼(SFW)空氣動力技術研究小組和DPW組織委員會合作設計開發，主要目的是為CFD的驗證和確認工作提供基準外形。CRM模型是典型的現代運輸機構型，設計馬赫數為0.85，升力系數為0.50。該模型包括了翼身組合體、翼/身/平尾組合體和翼/身/平尾/掛架/吊艙組合體等不同構型，DPW V組委會選擇了翼身組合體模型(CRM-WB)做為共同研究模型。CRM-WB計算構型見圖1，計算外形的基本參數見表1，其中Sref為參考面積，Cref為平均氣動弦長，b為展長，λ為梢根比，AR為展弦比，xref、yref、zref為力矩參考點的坐標。

本文采用有限差分方法離散任意坐標系下的RANS方程組，控制方程的對流項離散采用五階精度的WCNS格式，黏性項的離散采用六階精度中心格式，邊界及近邊界條件采用單邊四階精度離散，以上方法的詳細介紹見文獻[10-12]；湍流模型采用Menter剪切應力輸運(SST)兩方程模型[13]，離散方程組的求解采用BLU-SGS方法[14-15]。

ParameterAmericanunitSI-unitSref594720in2383.690m2Cref275.80in7.00532mb2313.5in58.7629mxref1325.9in33.67786myref468.75in11.90625mzref177.95in4.51993mλ0.2750.275AR9.09.0

2 高階精度計算網格

為了降低網格差異對數值模擬結果的影響，DPW組織委員會給出了網格生成指導原則[4-5]，對網格規模、物面第一層網格高度、邊界層網格增長率等網格參數進行了約定，并提供了基準的結構和非結構網格。由于高階精度格式對計算網格質量要求更高，本文研究中并沒有直接采用DPW V組委會提供的基準網格，而是根據網格生成指導原則，采用ICEM軟件重新生成了不同規模的粗、中、細、極細四套多塊對接結構網格。四套網格的詳細信息參見表2，其中，Nnum表示總的網格節點數，y+為第一層網格法向無量綱距離，nBL和λBL分別表示邊界層內法向網格數量和網格增長率，Nblock表示計算網格塊的數量。

圖2給出了CRM-WB計算構型的網格拓撲和表面網格(中等)，空間網格整體采用H型網格，在機翼和機身周圍分別包了一層O型網格，用于模擬邊界層流動。

表2 CRM-WB模型網格參數Table 2 Grid parameters of CRM-WB model

3 CASE1 計算結果與討論

DPW V組委會的CASE1狀態主要是開展網格收斂性研究。采用粗、中、細、極細四套網格和高階精度計算方法，開展了固定升力系數下網格規模對氣動力系數、典型站位壓力系數和翼身結合部后緣局部分離區大小等3個方面的影響，計算采用全湍流模擬方式。來流條件為：Ma=0.85，Re=5.0×106，CL=0.500±0.001。

1) 氣動力系數

表3給出了固定升力系數下，采用粗、中、細、極細四套網格得到的CRM-WB組合體的來流迎角α、阻力系數CD、壓差阻力系數CDp、摩擦阻力系數CDf和俯仰力矩系數Cm，同時給出了采用二階精度MUSCL(Monotonic Upwind Scheme for Conservation Laws)格式[16]的計算結果、風洞試驗結果[4,17]以及采用RE(Richardson Extrapolation)方法[4]外推得到的網格無關性結果。

從基于粗、中、細、極細四套網格的計算結果來看，采用高階精度方法得到的固定升力系數下的迎角、阻力系數、俯仰力矩系數均隨網格規模的增加而單調變化；而采用二階精度方法得到的阻力系數隨網格規模則出現了非單調的變化。從基于粗、中、細三套網格的計算結果來看，采用高階精度方法得到的阻力系數隨網格規模增加而單調增加，而采用二階精度方法得到的阻力系數則隨網格規模增加而單調減少，俯仰力矩系數的變化規律一致。

采用外推方法得到的高階精度阻力系數(0.024 92)計算結果略高于NASA Langley國家跨聲速設備(NTF)的風洞試驗結果(0.024 89)和Ames Langley 11 ft風洞的試驗結果(0.024 14)，這與文獻[5]給出的56組統計分析結果是一致的(中值為0.024 96，標準方差為0.000 53)；采用外推方法得到的高階精度方法俯仰力矩系數(-0.112 53)計算結果低于兩座風洞的試驗結果。

表3 CRM-WB模型的氣動特性(CL=0.500±0.001)Table 3 Aerodynamic characteristics of CRM-WB model (CL=0.500±0.001)

Rivers和Hunter等[18-19]采用非結構網格技術研究了風洞試驗模型的支撐機構和模型靜氣動彈性變形對CRM翼/身/平尾組合體(CRM-WBH)構型數值模擬結果的影響， Sclafani等[20]采用重疊網格技術研究了固定轉捩和模型靜氣動彈性變形對CRM翼身組合體(CRM-WB)構型數值模擬結果的影響。以上研究表明計算模型中考慮固定轉捩影響和靜氣動彈性變化使得阻力系數降低0.000 49、俯仰力矩系數增加0.01。因此，計算模型中考慮固定轉捩位置、模型靜氣動彈性變化、模型支撐裝置等因素可以進一步提高數值模擬結果與試驗結果的吻合程度。

2) 表面壓力系數

采用粗、中、細、極細四套不同規模網格和本文的高階精度方法，圖3給出了計算得到的CRM-WB模型機翼3個典型展向位置上壓力系數Cp分布曲線，同時給出了NTF風洞相鄰升力系數的測壓數據[16]，其中橫坐標x為機翼流向無量綱距離，η為機翼展向無量綱距離。

從圖3可以看出，在靠近翼根(η=0.131)位置，網格規模的變化對壓力系數分布基本沒有影響，并與試驗結果吻合良好。在機翼中部(η=0.502)的位置，除了粗網格外，其他三套網格的計算結果相近，且與試驗結果吻合良好。在機翼梢部(η=0.950)位置，上翼面30%弦長以前，四套網格計算結果相近，且負壓明顯高于試驗結果；30%～60%弦長之間，粗網格的計算結果與其他三套網格差異明顯；60%弦長之后，四套網格的計算結果相近，且與試驗結果吻合良好。總之，網格規模對壓力分布的影響由翼根到翼梢逐漸增加；中、細、極細三套網格下的計算結果相近；除靠近翼梢的站位外，計算得到的壓力分布與試驗測壓數據吻合良好。靠近翼梢位置計算結果與試驗結果之間的差異主要是由于計算模型中沒有考慮風洞模型的靜氣動彈性變形導致的。

3) 局部分離區

表4給出了高階精度方法下，采用不同規模網格得到的翼身結合部后緣分離泡的尺寸，其中ΔBL和ΔFS分別為分離泡的展向長度與縱向寬度(見圖4(a))。圖4為采用不同規模網格得到的翼身結合部后緣物面流線的變化。從表4和圖4可以看出，采用中等網格、細網格、極細網格得到的分離泡展向長度、縱向寬度與分離泡的尺寸隨著網格規模的增加而單調減少，而采用粗網格得到的分離泡尺寸不遵循上述變化規律。在網格收斂性研究的意義下，這說明粗網格的網格規模沒有進入數值解的漸進收斂區域。

表4 不同網格下分離泡尺寸Table 4 Separation bubble dimension with different grids

4 CASE2 計算結果與討論

DPW V組委會的CASE2狀態主要是開展抖振特性分析，在固定馬赫數下，采用中等網格模擬氣動特性隨迎角的變化。來流條件為：Ma=0.85，Re=5.0×106，α=2.50°、2.75°、3.00°、3.25°、3.75°、4.00°。基于中等網格，采用高階精度計算方法和二階精度方法，圖5給出了計算得到的氣動特性隨迎角變化曲線，圖6給出了采用高階精度方法得到的3.75°和4.00°兩個迎角下CRM-WB構型上表面流線。

由圖5可以看出，采用高階精度方法和二階精度方法得到的氣動特性具有明顯的差異。從圖5(a)看出，采用高階精度方法得到的失速迎角在3.75° 左右，而采用二階精度方法得到的失速迎角則提前到3.50° 左右；相同迎角下，采用高階精度方法得到的升力系數普遍高于采用二階精度方法得到的升力系數。從圖5(b)看出，失速迎角以前，相同升力系數下，采用高階精度方法得到的阻力系數小于采用二階精度方法得到的阻力系數。從圖5(c)看出，相同升力系數下，采用高階精度方法得到的低頭力矩系數大于采用二階精度方法得到的低頭力矩系數。

從圖6可以看出，在計算較大迎角范圍內，機翼上表面的激波誘導了機翼外側后緣的大范圍分離區，激波誘導分離區和翼身結合部后緣局部分離區的增加是導致4.0° 迎角下升力系數下降的主要原因。

5 結 論

1) 在固定升力系數下，采用粗、中、細、極細四套網格和高階精度計算方法得到了具有網格收斂性的氣動力結果并與DPW V統計分析結果吻合良好。

2) 在固定升力系數下，網格規模對機翼外側上表面的激波位置和翼身結合部后緣局部分離區大小略有影響。除了靠近翼梢的站位外，高階精度方法得到的典型站位的壓力分布與試驗結果吻合良好。

3)Ma=0.85時，采用高階精度方法得到的氣動特性與二階精度方法差異明顯。迎角4° 時，激波誘導分離區和翼身結合部后緣局部分離區的增加是導致升力系數下降的主要原因。

致 謝

感謝張玉倫、洪俊武、張書俊、李偉、楊小川等同志在高階精度格式程序實現方面所做的研究工作，感謝中國航空研究院白文博士在數據分析方面提供的幫助。

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

[6] 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: AIAA, 2008.

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

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

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

[10] 王運濤， 孫巖， 王光學, 等. DLR-F6翼身組合體的高階精度數值模擬[J]. 航空學報， 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).

[11] WANG Y T, ZHANG Y L, LI S, et al. Calibration of aγ-Reθtransition model and its validation with high-order numerical method[J]. Chinese Journal of Aeronautics, 2015, 28(3): 704-711.

[12] 王運濤， 孟德虹， 孫巖, 等. DLR-F6/FX2B翼身組合體構型高階精度數值模擬[J]. 航空學報， 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).

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

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

[15] 王光學， 張玉倫， 王運濤， 等. BLU-SGS方法在WCNS高階精度格式上的數值分析 [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).

[16] VAN LEER B. Towards the ultimate conservation difference scheme II, monoticity and conservation combined in a second order scheme[J]. Journal of Computational Physics, 1974, 14: 361-370.

[17] RIVERS M B, DITTBERNER A. Experimental investigation of the NASA common research model (invited): AIAA-2010-4218[R]. Reston: AIAA, 2010.

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

[19] 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: AIAA, 2012.

[20] SCLAFANI A J, VASSBERG J C, WINKLER C, et al. DPW-5 analysis of the CRM in a wing-body configuration using structured and unstructured meshes: AIAA-2013-0048[R]. Reston: AIAA, 2013.

(責任編輯：鮑亞平， 李世秋)

＊Corresponding author. E-mail: mdh157@163.com

High-order numerical simulation of CRM wing-body model

WANG Yuntao1, SUN Yan2, MENG Dehong1,*, WANG Guangxue1

1.ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China2.StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

High-order numerical simulation on CRM wing-body model is presented with the fifth-order WCNS scheme to assess the ability of high-order WCNS scheme on complex configuration simulation and the precision in predicating cruise drag of transonic configuration. Four grids (coarse, medium, fine, and extra fine) are created with software ICEM according to the gridding guidelines provided by DPW organizing committee, and the grid sizes range from 2 578 687 cells for the “Coarse” level to 65 464 511 cells for the “Extra-fine” level. Computation and analysis on four grids are carried out to investigate the grid effect on aerodynamic characteristics, pressure distribution and the local separation bubble at the wing root trailing edge, and the “Medium” grid is used in the numerical simulation and study of buffet onset. Compared to second-order numerical results, the statistic results submitted by DPW V participants and some experimental data, the high-order numerical results show that the drag coefficient computational results agree well with statistic data from DPW V participants; the grid density has some influence on the location of the shock wave and the size of the local separation bubble at the wing root trailing edge; the enlargement of the size of the separation zone due to shock wave and the local separation bubble at the wing root trailing edge on the upper surface of the wing is the main reason of the lift lift curve break at 4° angle of attack.

RANS equations; WCNS scheme; CRM model； flow simulation; grid density; aerodynamic characteristics

2016-04-07； Revised：2016-05-26； Accepted：2016-06-06； Published online: 2016-06-15 15:40

URL：www.cnki.net/kcms/detail/11.1929.V.20160615.1540.002.html

National Key Research and Development Program (2016YFB0200700)

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0185

2016-04-07； 退修日期：2016-05-26； 錄用日期：2016-06-06； 網絡出版時間：2016-06-15 15:40

www.cnki.net/kcms/detail/11.1929.V.20160615.1540.002.html

國家重點研發計劃 (2016YFB0200700)

＊通訊作者.E-mail: mdh157@163.com

王運濤， 孫巖， 孟德虹， 等． CRM翼身組合體模型高階精度數值模擬[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.

V211.7

A

1000-6893(2017)03-120298-08