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