基于非结构混合网格的CHN-T1标模气动特性预测

张耀冰, 唐 静, 陈江涛, 邓有奇

(中国空气动力研究与发展中心 计算空气动力研究所， 绵阳 621000)

0 引 言

近些年来CFD取得了很大的进步，包括网格生成技术、流场求解和高性能计算机等，评估三维真实运输机外形成为可能。CFD预测工业相关外形的气动性能取得了长足的进步[1-4],随之而来的CFD可信度、数值计算结果验证和确认以及数值结果的不确定度分析在国内外CFD业界受到广泛关注和高度重视。AIAA为此专门召开了六次阻力预测会议[5-7]和三次高升力预测会议[8]，邀请世界范围内的大学、研究所和工业部门共同评估当前CFD方法预测运输类飞行器气动力和力矩系数的发展水平，明确未来的发展方向。

国内也多次召开了相关的研讨会，中国空气动力研究与发展中心(CARDC)与中国航空工业空气动力研究院联合组织了DLR-F4翼身组合体、NLR7301两段翼型数值模拟研讨会和CT-1标模的大迎角数值模拟技术研讨会[9]，2010年和2013年分别召开了第一、二届航空CFD可信度开放式专题研讨会，极大地促进了国内CFD的发展[10-13]。

非结构网格因为其比较容易处理复杂外形而得到了广泛应用。对于复杂外形来说，生成非结构混合网格的时间比生成多块对接结构网格少很多。因此计算准备时间大大减少，而且只需要较少的人为干预。非结构另一个吸引人的方面是可以方便地使用基于流场的网格自适应技术[14]。超过一半的AIAA阻力预测会议参与者使用了非结构网格技术[6]。

为了评估国内航空CFD软件的技术状态，明确CFD方法与软件的下一步发展方向，推进国内CFD验证与确认工作稳步发展，为大飞机研制提供技术参考，CARDC组织召开了第一届CHN-T1标模CFD可信度研讨会(AeCW-1)，从基本气动力预测、气动弹性影响、网格技术等方面，对各种CFD求解器进行统一的确认研究，促进并发展CFD对复杂流动现象的计算模拟能力。为了评估自行研制的基于非结构混合网格的流场解算器程序MFlow[15-19]对飞行器力和力矩的预测能力，课题组参加了这次会议。本文是对课题组完成的工作，即CHN-T1标模的非结构混合网格生成和采用MFlow程序完成的计算结果的分析和总结。

1 计算模型及网格

计算外形为CARDC研制的用于风洞试验和CFD可信度验证的，具有窄机身超临界机翼特征的运输机标模CHN-T1。该模型包含机身、机翼、平尾、垂尾、短舱、翼吊、起落架整流包等部件，详细参数见文献[20]。在研讨会的标模计算时，不考虑短舱和翼吊，如图1所示。该模型的机身采用单通道窄体尺度，代表目前中短航程商用运输机的机身特征。机翼采用亚声速高气动效率的超临界翼型下单翼。设计巡航状态是Ma=0.78，CL=0.5。

图1 CHN-T1外形Fig.1 Configuration of CHN-T1

AeCW-1有两个必算算例：

Case1: 网格收敛性研究

Ma=0.78，CL=0.500(±0.001);

采用粗、中、细三套网格，网格量从约600万到5 000万。

Case2: 抖振特性研究

Ma=0.78;

α=-2°,-1°,0°,1°,2°,3°,3.5°,3.75°,4°,4.25°,4.5°;

采用中等规模网格;

Case2a：常规外形Config1;

Case2b：尾支撑装置影响外形Config2;

Case2c：尾支撑加静气动弹性影响外形Config3;

Config1、Config2和Config3的外形比较如图1所示。

所有计算状态的雷诺数均为Re=3.3×106(基于平均气动弦长)，参考温度为300 K。

气动系数的参考量见表1。

表1 气动系数参考量Table 1 Reference quantity of aerodynamic coefficient

AeCW-1的网格生成准则如下：

(1) 网格收敛性算例:

CHNT-1 Wing-Body-Tail外形要求3套不同层次的网格。

(2) 网格生成指南:

a. 物面法向网格尺度:

◆ 粗网格y+<1.0,

◆ 中等网格y+<2/3,

◆ 细网格y+<4/9;

b. 建议物面法向前两层网格尺度保持不变;

c. 网格收敛算例的网格量增长率为3倍;

d. 对于结构网格，在每个坐标方向网格增长率保持1.5倍;

e. 网格收敛算例的网格必须保持相同的网格分布，即保持相同的拉伸因子、相同的拓扑结构等等;

f. 附面层法向增长比<1.25;

g. 远场取100倍的平均气动弦长;

h. 对于中等网格：

◆ 机翼前后缘的弦向网格尺寸约为0.1%的当地弦长,

◆ 机翼翼根展向网格尺度约为0.1%的半展长,

◆ 机翼翼尖展向网格尺度约为0.1%的半展长,

◆ 机身头部和后体约为2%的平均气动弦长;

i. 对于粗网格和细网格，以上的值应该做相应的缩放;

j. 机翼后缘网格：

◆ 粗网格至少有8个单元,

◆中等网格至少有12个,

◆细网格至少有16个;

k.中等网格的网格量应满足工业应用阻力预测的要求;

l. 满足多重网格计算要求;

m. 对于结果网格，粗网格约为六百万个点，中等网格约为一千五百万个点，细网格约为五千万个点。对于基于节点型解算器的非结构网格，网格量要求与结构网格相同，网格空间尺度的要求为内部节点间的距离。对于基于格心型解算器的非结构网格，网格空间尺度的要求为格心或面心的距离，网格点的数量约为上面的要求的1/3。

计算网格由作者使用网格生成软件Pointwise生成，为三棱柱/四面体/金字塔组成的非结构混合网格，这套网格适合用于格心型流场求解器使用。图2为表面网格和空间网格示意图，在物面附近的附面层区域使用三棱柱网格，空间采用四面体，两者交界的区域使用金字塔作为过渡。

图2 CHN-T1网格Fig.2 CHN-T1 grid

网格收敛性研究需要的基础网格有3套，基准网格为中等网格，然后在此基础上分别进行细化和粗化，得到细网格和粗网格。图3是粗、中、细三套网格在机翼剖面η=0.17处空间网格的比较。由图可见，从物面到远场，网格过渡平缓，没有突变。从粗网格到细网格，网格整体加密。在附面层区域有足够多的三棱柱单元，附面层第一层网格的高度变小，法向增长率变小，而三棱柱区域的高度基本保持不变。在机翼的前后缘，网格加密，达到网格生成指南要求的当地弦长的0.1%的尺寸。

根据会议提供的网格生成指南，网格量的比例因子为3，粗网格的y+<1.0。网格的一维单元数比例因子为31/3≈1.44。远场取为100倍的平均气动弦长，约为20 m。如图4所示，在翼前后缘及翼根翼梢采用各向异性的三角形面网格单元，在满足模拟精度的同时尽量减少网格量，提高计算效率。

图3 η=0.17处的网格剖面Fig.3 Grids at η=0.17

图4 各向异性三角形网格
Fig.4 Anisotropic triangle grid

由于网格生成软件的技术限制，未能满足网格指南中(2)b条的要求，即附面层前两层法向取相同的尺寸。网格信息见表2。

2 数值方法

本文计算采用课题组自行发展的雷诺平均NS方程流场求解器MFlow。MFlow是基于非结构混合网格的亚跨超声速流场解算器，可以处理任意形状的网格单元，具有较强的灵活性。采用有限体积法对控制方程进行空间离散，未知变量存储于网格单元的体心。

表2 网格信息Table 2 Information of grids

在本文计算中，对流通量采用Roe通量差分裂方法进行离散，具有很高的间断和黏性分辨率，熵修正采用课题组自有的改进型Harten熵修正[15],修正参数取0.025。单元内使用线性重构使得空间离散具有二阶精度，使用Venkatakrishnan限制器来抑制间断附近的振荡，限制参数取50。变量梯度求解使用节点型Gauss方法[15,18]。黏性通量采用中心差分离散[15]。使用当地时间步长的一阶欧拉后差来达到定常状态。Jacobian通量使用一阶迎风格式得到，采用LU-SGS方法求解离散方程，CFL数取200。使用基于FMG(Full Multigrid)方法的4重“V”型多重网格方法加速收敛，粗网格的步数均为2000步。计算采用基于MPI的大规模并行计算技术加速收敛,CPU核数为128核。

计算假定流场为全湍流流场，湍流模型主要采用原始的SA一方程模型，同时采用QCR修正的SA模型计算了Case2状态作为对比。

如图5所示，在本文所有的计算中，收敛情况都比较好。密度的残差都下降了十个量级以上，气动力系数前六位有效数字保持不变。

图5 密度残差和力系数收敛历程Fig.5 Convergence history of density residual and force coefficient

3 网格收敛性研究

网格收敛性研究的主要目标是估算力和力矩的网格收敛解，也就是估计当网格量趋近于无穷大时力和力矩的值。

表3是本文三套网格的计算结果、CARDC FL-26风洞试验结果[21]和文献[22]中采用网格量为104亿的结构网格计算结果的比较。计算的阻力系数与试验结果比较接近，细网格的阻力系数与试验结果相差约10.2 counts(1 counts=0.0001)。俯仰力矩系数两者差异较大，根据Case2的计算结果分析，可能是尾支杆的影响造成的。迎角相比，细网格比试验结果小0.1568°。

“∞”表示采用中等网格和细网格线性插值得到的网格尺度无穷小，网格量无穷大的结果。与“∞”相比，细网格的阻力系数相差约-3.7 counts。“∞”与文献结果相比，总阻力相差-3.038 counts，压差阻力相差3.556 counts，摩擦阻力相差-0.517 counts，俯仰力矩相差-0.000 997 7，迎角相差0.114 37°。

图6是总阻力、压差阻力、摩擦阻力、俯仰力矩和迎角的网格收敛图。横坐标是N-2/3，N代表网格单元数。采用N-2/3是基于MFlow程序的数值方法为二阶精度，对于由粗到细的一套自相似网格，表示其网格尺寸的二次方，因此直线就表示空间网格收敛为二阶精度，N-2/3为0时表示网格量为无穷大。

由图可见，本文计算结果当网格量趋于无穷时力系数和迎角的收敛都是单调的，而且收敛精度接近二阶，证实流场解位于渐进解区域。压差阻力CDP随着网格加密减小，摩擦阻力CDV略有增大，因此总阻力减小。摩擦阻力随着网格变化较小，表明当附面层网格加密到一定程度(y+<1.0)后可以准确地捕捉边界层特性。

由于计算网格不同、解算器的计算参数等原因，文献[22]的气动特性计算结果与本文计算结果的网格收敛曲线不符，两者有一定的差距。

图7是网格对CHN-T1模型机翼的压力分布的影响，包括三个典型展向站位，分别为0.17、0.50和0.85，包含了从翼根到翼尖各个不同的区域。图8是平尾的三个展向站位，即0.28、0.50和0.95。

该状态(CL=0.5)下，机翼上的流动基本没有分离。随着网格加密，压力分布变化很小。网格密度对剖面压力分布的影响主要体现在吸力峰值和对激波的分辨能力上，在机翼的其他位置相差很小。在激波附近，网格越密，激波越陡。平尾的上翼面压力大于下翼面压力，从而产生负升力，对飞机起抬头作用。网格加密对平尾压力分布的影响很小。

表3 case1状态的CFD计算结果与试验结果Table 3 Results of CFD and experiment for case1

图6 网格收敛特性Fig.6 Grid-convergence properties

图7 不同网格机翼剖面压力分布的比较Fig.7 Comparison of Cp prediction with grid refinement of wing section

图8 不同网格平尾剖面压力分布的比较Fig.8 Comparison of Cp prediction on horizontal tail with grid refinement

图9是翼身结合处物面附近的分离情况比较。CFX表示表面摩擦系数在流向的分量，经常被用于分析近壁区的流动分离情况。图中蓝色区域表示CFX<0的位置，即存在分离的区域。随着网格的加密，模拟得到的翼身结合处分离气泡尺度越大，总的来说，分离泡的尺寸很小，流向长度不超过当地弦长的2.5%。

4 抖振特性研究

图10显示了外形变化对升力系数、阻力系数和俯仰力矩系数的影响。

升力和阻力系数与试验结果[21]符合得较好，Config1外形的俯仰力矩系数与试验相比有一个较大的平移量，约为0.04左右，Config2和Config3的俯仰力矩系数与试验比较接近。

考虑尾支撑(Config2)时，升力系数略有增大；阻力系数变化很小，在10 counts以内，其中压阻系数增大，摩阻系数减小；俯仰力矩系数有一个低头力矩增量，且随迎角增大增量值减小。

考虑静气动弹性变形(Config3)的情况下，机翼的实际迎角减小，造成升力系数减小，且随迎角增大，减小的量值增大，在3.5°以后，量值变化趋势改变；阻力系数减小，其中压阻系数减小明显，摩阻系数变化很小；俯仰力矩系数有一个很小的抬头力矩增量。

图11和图12是Config1和Config2外形的机翼和平尾压力分布剖面的比较，迎角为2°。由图可见，尾支撑对机翼的影响较小，对平尾的影响很明显。由于尾支撑的存在，平尾产生的负升力明显减小，从而飞机的总升力增大，由于平尾位于飞机质心之后，负升力的减小产生了一个低头俯仰力矩增量。

图13是Config2和Config3外形的机翼压力分布剖面比较，迎角为3°。由于静气动弹性变形，机翼上翼面的激波位置前移，吸力峰值变小，从而造成机翼产生的升力减小，俯仰力矩产生一个小的抬头增量。

图14和图15是Config1外形上下翼面不同迎角下的分离情况。在迎角3.5°时，机翼上翼面在激波后出现明显的沿展向分布的小分离气泡，在4°时，变为大分离，从而引起升力系数曲线发生拐折，升力线斜率变小。在迎角为-2°时，下翼面产生明显的大分离。

图9 翼身结合处的分离气泡比较
Fig.9 Comparison of separation bubble on wing body junction

图10 不同构型气动特性比较Fig.10 Comparison of aerodynamics characteristic among different configurations

图11 机翼压力分布比较
Fig.11 Comparison ofCpprediction of wing section

图12 平尾压力分布比较
Fig.12 Comparison ofCpprediction of horizontal tail

图13 机翼压力分布比较
Fig.13 Comparison ofCpprediction of wing section

图14 上翼面分离情况
Fig.14 Separation on wing upper surface

图16是Config1和Config2在迎角3°时平尾的压力分布云图。由于尾支撑的存在，平尾下翼面靠近机身的低压区范围明显变小，压力值增大；而上翼面的低压区范围变大，压力值减小。同时尾支杆对垂尾也有一定的影响，使垂尾的低压区范围变大，压力值减小。

图15 下翼面分离情况
Fig.15 Separation on wing lower surface

图16 尾支撑对平尾压力分布的影响Fig.16 Influence of horizontal with sting on Cp prediction

5 QCR修正的影响

表4是QCR修正SA模型和原始SA模型计算的Config1外形的升力系数、阻力系数和俯仰力矩系数的差量。在小迎角时，两者计算结果很接近，当迎角变大(大于3.5°)时，两者差异的差量明显增大。

表4 QCR修正差量Table 4 Increment of QCR on the prediction of aerodynamic characteristic

6 结 论

本文使用自行发展的非结构混合网格流场解算器MFlow对AeCW-1提供的CHN-T1客机标模进行了数值模拟，研究了网格收敛性、气动特性曲线、压力分布和表面流动分离情况，分析结果显示，MFlow解算器对CHN-T1标模具有较好的模拟能力。具体结论如下：

(1) 本文计算得到了近似线性的网格收敛特性。

(2) 网格加密对激波和分离气泡的模拟精度有比较明显的影响。

(3) 支撑系统对平尾压力分布和俯仰力矩系数预测影响较大，对机翼影响较小。

(4) 考虑静气动弹性之后，机翼的当地迎角变小，压力预测捕捉到了这一变化。

(5) QCR修正的SA模型与原始SA模型计算结果在小迎角时相差较小，大迎角下有明显差别。