CRM-WBPN风洞试验模型数值模拟

2019-03-04孟德虹李伟王运涛孙岩

航空学报 2019年2期

孟德虹，李伟，王运涛,*，孙岩

1. 中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000 2. 中国空气动力研究与发展中心空气动力学国家重点实验室，绵阳 621000

AIAA阻力预测会议(Drag Prediction Workshop, DPW)已经连续成功举办了6届[1-6]，获得了世界范围内相关研究机构的广泛关注，已经成为CFD(Computational Fluid Dynamics)验证与确认研究领域最重要的国际合作之一。DPW系列会议的宗旨是评估各种CFD方法和软件在典型运输机构型气动特性预测方面的现状，明确CFD技术的发展方向。从2009年6月召开的第4届DPW(DPW IV)开始[7]， CRM(Common Research Model)构型[8]被DPW组委会选择为基准研究模型，并先后在多座低温跨声速风洞完成了多种组合构型的风洞试验，包括美国NASA Langley的NTF(National Transonic Facility) 风洞和NASA Ames的TWT(Transonic Wind Tunnel)风洞，欧洲的ETW(European Transonic Wind tunnel)风洞，试验结果包括了气动特性、表面压力分布及模型变形测量等[9-10]，以上风洞试验结果为CFD的确认工作提供了高质量的对比数据。

2012年6月召开的第5届DPW(DPW V)采用CRM翼身组合体(CRM-WB)构型作为基准研究模型，来自世界范围的22家研究机构共提供了57组网格收敛性计算结果和50多组抖振特性计算结果。通过与NASA NTF风洞和TWT试验结果的对比分析，上述计算结果与相应的风洞试验测力、测压试验结果均存在较大差异，尤其是升力特性与力矩特性差异明显。针对上述问题，Rivers等[11-12]采用非结构网格计算软件USM3D 6.0研究了CRM翼/身/平尾组合体(CRM-WBH)风洞试验模型的支撑机构和静气动弹性变形对数值模拟结果的影响。David[13]采用结构网格计算软件elsA研究了CRM翼身组合体(CRM-WB)风洞试验模型静气动弹性变形对数值模拟结果的影响，Keye等[14]采用非结构网格技术和流固耦合方法研究了静气动弹性变形对CRM-WB模型数值模拟结果的影响，王运涛等[15]采用结构网格技术和流固耦合方法研究了静气动弹性变形和模型支撑对CRM-WB模型数值模拟结果的影响，并采用高阶精度格式模拟了CRM-WB静气动弹性模型[16]，研究表明：计算模型中考虑风洞模型支撑装置和静气动弹性变形可以显著提高计算结果与试验结果的吻合程度，但上述研究模型均未包含挂架短舱。2016年召开的第6届DPW(DPW VI)将CRM翼/身/架/舱(CRM-WBPN)模型列为研究模型之一，主要评估CFD模拟复杂构型及挂架短舱阻力增量的能力，会议采用的CRM-WBPN模型仅考虑了机翼的静气动变形影响，没有考虑模型支撑对该构型的影响。

本文作者在前期工作的基础上[17]，采用多块对接结构网格技术和亚跨超CFD软件平台(TRIP3.0)[18-19]，基于DPW VI组委会提供的经过静气动弹性修正的CRM-WBPN模型，在网格收敛性研究的基础上，研究了模型支撑对挂架短舱阻力增量及总体气动特性影响。通过与NASA Langley NTF风洞试验数据以及CRM-WB模型数值模拟结果的比较，获得了一些有价值的研究结论。

1 CRM翼/身/架/舱风洞模型

CRM模型是典型的现代运输机构型，设计马赫数Ma=0.85，升力系数CL=0.50。该模型包括了翼身组合体、翼/身/平尾组合体和翼/身/挂架/吊舱组合体等不同构型。第6届DPW组委会选择了CRM-WB模型和CRM-WBPN模型作为研究模型，采用ETW风洞试验中测量的模型变形数据修正了CRM-WB模型与CRM-WBPN模型。为以下叙述方便，将包含静气动弹性影响的CRM翼/身/架/舱构型简称为CRM-WBPN，包含模型支撑和静气弹变形的CRM翼/身/架/舱构型简称为CRM-WBNPS。本文采用了与文献[14]中CRM-WBS相同的支撑处理方法，如图1所示，图中紫色部分为简化的支撑部分，绿色部分是为了计算中减小分离而增加的导锥。图2给出了迎角α=2.75°时，ETW风洞试验测量得到的机翼弯曲和扭转变形的展向分布，图中：η为无量纲展向站位，dy表示弯曲变形位移，Bengding LE和Bending TE分别为机翼前后缘变形量，dθ表示机翼剖面的扭转角，定义迎角变大扭转角为正。由图2可以看出，机翼弯曲变形沿翼展方向逐渐增大，而由于机翼后掠角的影响，机翼后缘的弯曲变形要大于前缘变形，因此向上的弯曲变形会引起机翼剖面产生负扭角。CRM构型计算外形的基本参数参考文献[8]，本文不再详细描述。

图1 CRM翼/身/架/舱计算模型Fig.1 Computational model for CRM wing/body/pylon/nacelle configuration

图2 CRM构型机翼弯曲和扭转变形沿展向分布(α=2.75°)Fig.2 Blending and torsional deformation distributionalong the wing span of CRM configuration (α=2.75°)

2 计算网格与计算方法

根据DPW VI组委会的网格生成指导原则，在文献[17]的基础上，采用商业软件进一步生成了不同规模的小、粗、中、细4套CRM-WBPN模型的多块对接结构网格，4套网格的详细信息见表1。表中：Nnode表示网格节点总数；nBL和λBL分别表示边界层网格数量和网格增长率；y+为第1层网格法向无量纲距离；Nblock为计算网格块的数量。

表1 CRM-WBPN模型网格参数Table 1 Grid parameters of CRM-WBPN model

采用与上述网格类似的网格拓扑及网格分布，根据CRM网站提供的模型支撑外形，进一步构造了CRM-WBPNS模型的计算网格(图3)，半模规模达到了157 673 088个节点。

数值模拟方法采用与文献[15，17]相同的方法，即采用亚跨超CFD软件平台求解雷诺平均Navier-Stokes(Reynolds-Averaged Navier-Stokes, RANS)方程，具体计算方法选择如下： RANS方程无黏项的离散采用二阶精度MUSCL(Monotonic Upwind Scheme for Conservation Laws)型ROE格式[20]，黏性项的离散采用二阶中心格式，湍流模型采用Menter’s SST(Shear Stress Transport)两方程模型[21]，离散方程组的求解采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel method)方法[22]，并采用多重网格技术和大规模并行技术加速收敛。数值模拟结果均采用全湍流方式计算，没有考虑流动转捩对计算结果的影响。

图3 CRM-WBPNS模型网格拓扑及表面网格(中等网格)Fig.3 Grid topology and surface grid of CRM-WBPNS model (medium grid)

3 网格收敛性

采用第2节的4套不同规模的计算网格，开展了CRM-WBPN模型固定升力系数下的网格收敛性研究。计算来流条件为：Ma=0.85，CL=0.5 ±0.000 1，Re=5.0×106。表2给出了采用小、粗、中、细4套不同密度的网格计算得到的CRM-WBPN模型的气动特性。可见迎角(α)、阻力系数(CD)、压差阻力系数(CDp)、摩擦阻力系数(CDf)和俯仰力矩系数(Cm)绝对值均随网格密度的增加而单调变化。从中等网格到密网格，迎角增加约0.01°、阻力系数减少0.3个阻力单位(1个阻力单位=10-4)、力矩系数减少0.000 1，这说明中等网格规模已基本消除网格依赖性，满足本文研究要求。

图4给出了来流条件Ma=0.85，CL=0.5，Re=5.0×106，不同网格密度下CRM-WB模型与CRM-WBPN模型的挂架短舱阻力增量ΔCD，其中CRM-WB模型的模拟结果来自文献[17]。可以看出，随着网格密度的增加，计算得到的阻力增量趋于定值。以下的研究中均采用中等网格。

表2 CRM-WBPN模型的气动特性 (CL=0.5±0.000 1)Table 2 Aerodynamic characteristics of CRM-WBPN model (CL=0.5±0.000 1)

图4 CRM模型挂架短舱阻力增量Fig.4 Drag increment of pylon/nacelle of CRM model

4 模型支撑对压力分布的影响

采用CRM-WBPN和CRM-WBPNS模型研究支撑装置对机翼和短舱压力分布的影响，来流条件为Ma=0.85，CL=0.5，Re=5.0×106。图5为2个模型表面压力分布等值线对比，从图中可以看出尾撑对机身、机翼上表面和短舱外侧压力分布影响较大，对机翼下表面和短舱内侧影响较小，模型尾支杆对模型后体气流有阻滞作用，使流速降低，压力增加，其影响区域前传至机翼，导致机翼上表面激波略有变化，从图6中机翼典型展向站位压力系数Cp的分布对比可以看出，上翼面的激波位置前移从翼根直到翼梢一直存在。

图5 CRM-WBPN构型表面压力分布等值线Fig.5 Surface pressure distribution contour of CRM-WBPN configuration

图6 CRM-WBPN构型典型展向站位压力系数分布Fig.6 Pressure coefficients distribution at different spa-nwise locations of CRM-WBPN configuration

图7为2个构型机翼剖面升力系数CLs展向分布与试验值比较。由图可以看出，模型支撑对机翼剖面升力系数略有影响。

图7 CRM-WBPN模型剖面升力系数展向分布Fig.7 Section lift coefficients distribution along the wing span of CRM-WBPN model

5 模型支撑对气动特性的影响

采用CRM-WB与CRM-WBPN、CRM-WBS与CRM-WBPNS 2组模型研究支撑装置对挂架短舱阻力增量的影响，来流条件与第4节相同。表3给出了在上述来流条件下，模型支撑对挂架短舱阻力增量的影响，由表可以看出，CL=0.5时，CRM-WBPN和CRM-WBPNS模型迎角增加，但2个构型的挂架短舱阻力增量基本接近，均在风洞试验的误差精度范围内。

采用CRM-WBPN和CRM-WBPNS模型研究不同来流迎角下，支撑装置对总体气动特性的影响。来流条件Ma=0.85，α=0°～4.0°，Re=5.0×106。需要说明的是，CRM-WBPN模型在不同迎角下的静气动弹性变形是不同的，本文的工作侧重于研究模型支撑对气动特性的影响，计算模型统一采用迎角2.75°测量得到的静气动弹性变形，这种做法与文献[10]相同。

图8给出了中等网格规模下，2种模型总体气动特性随迎角的变化，同时给出了NASA NTF风洞经过洞壁干扰修正的测力试验结果。从图8可以看出，模型支撑的对CRM-WBPN的影响与对CRM-WB的影响规律[14]相似，使得相同迎角下升力系数、阻力系数明显降低，俯仰力矩系数增加。采用CRM-WBPNS模型计算得到的升力系数和阻力系数的计算结果更加接近试验值；俯仰力矩系数计算结果与试验结果的吻合程度得到显著改善，但依然有较大差距。数值计算结果与风洞试验结果之间的差异需要从风洞试验数据的修正和数值计算方法2个方面进一步开展研究工作。

表3 CRM模型挂架短舱阻力增量Table 3 Pylon/nacelle drag increment of CRM models

图8 CRM-WBPN模型气动特性Fig.8 Aerodynamic characteristics of CRM-WBPN model

图9给出模型支撑对CRM-WBPN模型各部件气动特性的影响量。总体来看，模型支撑对机身和机翼气动特性影响较大，对短舱挂架的气动特性基本没有影响。从升力系数来看，模型支撑对机身和机翼的影响量基本相当，不同的是，支撑装置对机身影响量随迎角基本不变，对机翼影响量随迎角先增大后减小。从阻力系数和俯仰力矩系数来看，支撑装置对机身影响量比机翼的影响量要大，对机翼阻力系数影响量在2.5°迎角后变化较大。