DPW4翼／身／平尾组合体的数值模拟

王运涛，张书俊，孟德虹

（中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000）

0 引 言

随着大型网格生成软件、大规模并行计算技术和后置处理软件的发展，CFD几乎可以模拟所有高度复杂飞行器外形的绕流流场［1］，CFD已经与风洞试验一道，成为飞行器气动设计工作者最重要的研究手段。飞行器设计者使用CFD工具面临的首要问题是软件的计算精度和效率问题。多重网格技术、低速预处理技术和并行计算技术的综合应用，极大地提高了CFD软件的计算效率，但CFD的验证与确认问题依然是当前计算空气动力学研究的热点［2］。

为研究CFD模拟的精准度问题，国际上先后组织了许多专题研讨会，其中比较具有代表性的如欧洲的计算空气动力学研究项目 ECARP［3］（European Computational Aerodynamics Research Project）和AIAA的DPW（Drag Prediction Workshop）系列会议。为了开展典型运输机构型的阻力预测精度研究，AIAA阻力计算的工作小组先后在2001年6月、2003年6月、2006年6月召开了三次阻力计算的工作会议，选择的研究构型包括DLR－F4翼身组合体构型、DLR－F6翼／身／架／舱复杂组合体构型、DLR－F6翼身组合体及其修型构型等。该系列会议通过规定统一的计算状态、统一的计算网格和精细的试验结果，对世界范围内各家研究机构提供的CFD软件计算结果进行评估，明确CFD的发展现状和下一步的努力方向。

AIAA的第四届阻力计算的工作会议DPW 4（AIAA 4thDrag Prediction Workshop）于2009年6月召开。本次会议的基本宗旨与前三次会议相同，在总结前三次会议经验的基础上，在组织方式上做了一些调整。第一是重新设计了研究模型，即共同研究模型 CRM（Common Research Model）［4］，采用了现代运输机的翼／身／平尾组合体构型，设计马赫数为0.85；第二是事先不提供任何风洞试验结果，即采用盲评的方式组织此次专题研讨活动。

本文基于TRIP3.0（TRIsonic Platform V3.0）软件和多块对接网格技术，采用粗、中、细三套结构网格，选择SST两方程湍流模型和MUSCL型的ROE格式，通过求解任意坐标系下的雷诺平均NS方程（RANS），数值模拟了DPW－IV翼／身／平尾组合体构型的气动特性。数值模拟内容包括：网格收敛性研究、下洗影响和阻力发散马赫数计算。在与National Transonic Facility（NFT）的试验结果对比基础上，本文详细研究了网格密度、来流马赫数和平尾偏角对总体气动特性和压力分布的影响。通过研究发现TRIP软件能够较好地模拟该构型的气动特性；网格密度对压差阻力的影响较大，对摩擦阻力的影响较小；网格密度对压力分布的影响不明显。

1 TRIP3.0软件的简介

TRIP是中国空气动力研究与发展中心自行研发的大型通用CFD软件，已经成功推广到国内19家气动研究单位，极大地促进了国内CFD应用水平的提高。TRIP软件从0.0版本、1.0版本、2.0版本，发展到现在的3.0版本。TRIP3.0版本在兼容2.0版本［5］所有功能的基础上，与本项研究相关的几个方面的改进包括：

·改进了多重网格收敛加速技术，克服了多重网格技术不能在并行计算中使用的缺点；

·对算法进行了优化，减少了冗余运算或重复运算，每一步的计算效率同比提高30%以上；

·提高了SA一方程湍流模型和SST两方程湍流模型的求解精度，增强了软件的鲁棒性；

·增加了双精度计算选项，进一步提高了数值模拟的精度。

2 计算网格

本文采用的计算网格由Boeing公司提供［6］，网格结构为多块对接网格（1－to－1），网格分为粗网格（Coarse）、中等网格（Medium）和细网格（Fine）三种。这三套网格的y＋分别约为1、1／3和4／9，壁面第一层网格分别不大于0.0375mm、0.025mm和0.0166mm。全机网格的详细信息如表1所示，网格拓扑和表面网格如图1所示。

计算采用的平均气动弦长MAC＝7.00532m，参考面积Sref＝383.6895552m2，力矩参考点：Xm＝33.67786m，Ym＝4.51993m，Zm＝－11.90625m。

图1 翼／身／尾组合体网格拓扑以及表面网格（中等）Fig.1 Surface grid for DPW4wing／body／tail（Medium）

3 计算结果和分析

本文的数值模拟采用TRIP3.0软件，选择二阶精度的通量差分（FDS）类型的MUSCL差分格式［7］，Menter’sk－ωSST 两方程湍流模型［8］，综合采用局部时间步长、三重网格加速收敛和大规模并行计算技术，数值模拟了两种工况的流场，其中 WBH0、WBHm2和 WBHp2的平尾迎角分别为0°、－2°和2°，WBHNo为CRM不带平尾的模型。

3.1 CASE 1计算结果与分析

3.1.1 网格收敛性研究（CASE1a）

网格收敛性研究的计算构型为 WBH0，平尾偏角为0°，数值模拟采用粗网格、中等网格、细网格，计算状态如下：

图2给出了WBH0组合体采用粗网格、中等网格和密网格的阻力系数、摩阻系数和压差阻力系数的计算结果，其中横坐标为网格点数的－2／3次幂。表2给出了CL＝0.500（±0.001）翼身组合体对应的迎角、阻力系数及其分量，可以看出，采用上述三套网格，在固定升力系数下，迎角、总阻力系数、压差阻力系数、摩擦阻力系数及俯仰力矩系数均是单调变化的。从粗网格到密网格，迎角从2.305°增加到2.32°，总阻力系数从0.02722减少到0.02682，压差阻力系数从0.01218增加到0.01232，摩擦阻力系数从0.01504减少到0.01450，俯仰力矩系数从－0.05196增加到－0.04357。图中同时给出了文献［9］的计算结果。文献［9］中采用了非结构网格，粗、中、细三种网格规模分别为410万、1170万和3410万，空间离散格式采用二阶中心格式，湍流模型采用SST两方程湍流模型。

根据Richardson外推公式，本文得到的固定升力系数下的迎角、阻力系数、力矩系数的网格无关性结果分别为：2.340、0.02701和－0.04025，与文献［10］的对比如表3所示。文献［10］提供了19家单位29个计算结果的算术平均值及标准方差，本文的计算结果均落在文献［10］的计算误差带内。文献［10］同时给出了试验结果的平均值及误差，与相应的试验结果的平均值相比较，阻力的计算误差接近0.0008，迎角相差0.687°，力矩出现了符号的差异。导致计算结果与风洞试验结果的差异的主要原因首先应该是风洞试验模型添加了转捩带，而数值模拟采用了全转捩的模拟方式；其次，是风洞试验中采用了垂直尾部支撑的方式，试验模型尾部的支撑形式对气动特性有明显的影响，对力矩特性的影响尤其显著。具体的原因需要通过更加精细的试验和数值模拟进行分析。

图2 DPW4翼／身／平尾组合体的网格收敛性研究Fig.2 Grid convergence study for DPW4wing／body／tail

表2 WBH0翼／身／平尾组合体的网格收敛性研究Table 2 Grid convergence study for WBH0wing／body／tail

表3 WBH0组合体外推结果与文献的比较Table 3 WBH0data extrapolated to continuum

图3给出了平尾迎角为0°的翼／身／尾组合体分别采用粗网格、中等网格和密网格，CL＝0.500±0.001时，平尾附近的压力分布及表面流线云图。可以看出，不同的网格密度下均存在一定的分离区，三套网格给出的流线形态基本一致，密网格得到的平尾后部的分离区要小一些。图4给出了其在13.06%、28.28%、39.71%、50.24%、60.28%、72.68%、84.56%、95%等8个典型站位上的压力分布比较。可以看出，在机翼典型站位上，粗网格、中等网格和密网格的压力分布几乎完全重合。

图3 压力分布及表面流线云图（CL＝0.500±0.001）Fig.3 DPW4wing／body／tail surface streamline

图4 典型站位压力分布（CL＝0.500）Fig.4 Cpdistributions for typical wing stations

3.1.2 配平影响（CASE1b）

计算状态如下：

·M＝0.85

·α＝0°，1.0°，1.5°，2°，2.5°，3.0°，4.0°

·平尾迎角：－2°，0°，2°，无平尾

·Re＝5.0×106（基于MAC＝7.00532m）

基于中等规模的网格，进一步研究比较了M＝0.85来流状态下，平尾迎角分别为－2°、0°、2°和无平尾翼身组合体的气动特性。根据平尾迎角为－2°、0°和2°时，Cm～CL曲线，可以插值得到不同升力系数下平尾的配平迎角，如图5所示。升力系数小于0.34时，配平迎角为正值，绝对值随着升力系数的增加而逐渐减小；升力系数大于0.34时，配平迎角为负值，绝对值随着升力系数的增加而增加，在升力系数0.6时，配平迎角绝对值达到最大；升力系数进一步增加，配平迎角绝对值逐渐减小。

图5 下洗研究Fig.5 Downwash study

根据配平迎角随升力系数的变化曲线和阻力系数随升力系数的变化曲线，进一步得到配平后阻力系数随升力系数的变化曲线，翼身组合体与配平后的阻力差量变化趋势见图6，升力系数小于0.42时，配平后引起的阻力增量基本不变，量值在20个阻力单位（1阻力单位＝0.0001）；升力系数大于0.42后，配平后引起的阻力增量随升力系数的增加而基本线性增加，在升力系数0.64左右达到最大，量值为55个阻力单位。

图6 不同升力系数下平尾配平引起的阻力增量Fig.6 Delta effects of drag and and corresponding lift values

3.2 CASE 2计算结果及分析

计算状态如下：

·M＝0.70，0.75，0.80，0.83，0.85，0.86，0.87

·CL＝0.4，0.45，0.5（±0.001）

·中等网格

·Re＝5.0×106（基于MAC＝7.00532m）

基于中等规模的网格，在不同升力系数下，WBH0的阻力特性随马赫数变化趋势如图7所示。可以看出，马赫数小于0.83，阻力系数变化平稳；马赫数继续增大，阻力特性则迅速恶化。

图7 阻力系数随来流马赫数变化曲线（ih＝0）Fig.7 Drag rise curves for various Mach numbers（ih＝0）

4 结 论

本文采用TRIP（TRIsonic Platform）软件和结构网格技术，数值模拟了DPW4（AIAA 4thDrag Prediction Workshop）翼／身／平尾组合体构型，详细研究了网格密度、马赫数、平尾迎角对总体气动特性和压力分布的影响，通过与相应的试验结果和计算结果的比较，得到以下一些基本结论：

（1）本文的计算结果与DPW4计算结果的平均值取得了较好的一致，风洞试验结果与计算结果在迎角、力矩特性方面的差异有待于进一步分析；

（2）网格密度对总阻力系数有一定的影响，相比较而言对压差阻力的影响大一些，对摩擦阻力系数影响较小；

（3）在合理的网格分布前提下，网格密度对压力分布的影响不明显，粗网格、中等网格和密网格的压力分布几乎完全重合。

致谢：感谢TRIP软件的其他开发人员，张玉伦、洪俊武和王光学等同志在多重网格技术和湍流模型方面的研究工作。

［1］TINOCO E N，BOGUE D R.Progress toward CFD for full flight envelope［J］.AeronauticalJournal，2005：451－460.

［2］WILLIAM L O，TIMOTHY G T.Verification and validation in computational fluid dynamics［J］.Progressin AerospaceSciences，2002，38：209－272.

［3］HAASE W，CHAPUT E，LESCHZINER M A.The European high lift project EUROLIFT Ⅱ－objectives，approach，and structure［R］.AIAA 2007－4296.

［4］VASSBERG J C，DEHAAN M A，RIVERS S M，et al.Development of a common research model for applied CFD validation studies［R］.AIAA 2008－6919.

［5］王运涛，王光学，张玉伦.采用TRIP2.0软件计算DLR＿F6构型的阻力［J］.空气动力学学报，2009，27（1）：108－113.

［6］ftp：／／cmb24.larc.nasa.gov／outgoing／DPW4／multiblock＿boeing

［7］Van LEER B.Towards the ultimate conservation differencescheme II，monoticity and conservation combined in a second order scheme［J］.J.Comp.Phya.，1974，14：361－370.

［8］MENTER F R.Improved two－equationk－ωturbulence modelsfor aerodynamic－flows［R］.NASA TM－103975，1992.

［9］ELIASSON P，PENG S H.Influence of turbulence modelling and grid resolution in computations of the DPW－4CRM configuration［R］.AIAA 2010－1416.

［10］VASSBERG J C，TINOCO E N，et al.Summary of the fourth AIAA CFD drag prediction workshop［R］.AIAA 2010－4547.