王运涛，王光学，张玉伦

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

0 引言

随着大型专业前置处理软件、计算流体力学(CFD)技术、后置处理软件和计算机硬件技术的飞速发展，基于RANS方程(Reynolds Averaged Navier-Stokes)的CFD软件已经可以模拟高度复杂飞行器外形的绕流流场［1］。市场竞争、安全性和商业利润等多重的压力，使得CFD成为现代飞行器成功设计的重要因素，已经与风洞试验一道成为飞行器气动设计的最重要的研究手段［2］。我国中长期科学和技术发展规划纲要中，大飞机的研制已经被列为重大专项，在大飞机的气动设计过程中，CFD必将发挥重要的作用。

尽管飞行器设计者已经充分意识到了CFD应用的巨大潜力，但CFD软件尚没有像计算结构力学(CSM)软件那样被认为是一个成熟的工具，原因主要包括流动控制方程中包含各种假设、计算效率和计算精度需要进一步提高、复杂外形的计算结果强烈依赖于使用者的经验和网格质量等方面。因此，CFD软件的验证和确认问题(Verification＆Validation)依然是当前研究的热点［3］。为研究CFD模拟的精准度问题，国际上先后开展了许多专题研究，如 ETMA(Efficient Turbulent Models for Aeronautics)、ECARP(European Computational Aerodynamics Research Project)、AVTAC(Advanced Viscous Flow Simulation Tools for Complete Civil Transport Aircraft Design)等，其中比较具有代表性是AIAA的DPW(Drag Prediction Workshop)系列会议。

针对运输机构型，为研究CFD的阻力计算精度问题，明确CFD技术的现状和未来的努力方向，AIAA阻力计算工作小组在2001年6月召开了第一次工作会议(DPW I)，该次会议选择了DLR-F4翼身组合体作为标准算例，在本次会议上18家单位采用14中软件提供了计算结果［4］，2003年6月召开了第二次工作会议(DPW II)，该次会议选择了DLR-F6翼/身/挂/舱组合体作为算例，会议的重点是挂架/吊舱的安装阻力，在本次会议上共有22家研究机构提供了20种CFD软件的计算结果［5］，试验结果是90年代在法国ONERA S2MA 1.77m×1.75m跨声速风洞中完成的。2006年6月举办了第三次工作会议(DPW III)，此次会议提供了两类共四种构型，一类是翼身组合体基本构型(DLR-F6)及修型构型(DLR-F6_FX2B)［6］，一类是机翼的基本构型(DPW_W1)及简单优化构型(DPW_W2)。包含机翼构型的主要目的是鼓励学术机构参与研讨和开展网格收敛性的研究。本次会议的目的与上两次相同，但更强调外形的局部修改对总体气动特性的影响，组织形式也有所变化，最突出的不同在于没有事前提供相应的试验结果，只有不同软件的计算结果可供比较。本文采用的试验结果是2007年9月在NASA的国家跨音速设备(NTF)中完成的。

TRIP(TRIsonic Platform)是中国空气动力研究与发展中心自行研发的CFD软件，该软件采用结构网格技术和有限体积方法，通过数值求解三维任意坐标系下的雷诺平均NS方程，数值模拟亚跨超飞行器的气动特性和绕流流场。有关该软件的基本功能介绍请参考文献［7］。为进一步提高TRIP软件的计算效率，适应千万量级的网格规模的计算，满足精细化流场数值模拟需求，课题组成员又进一步开发了TRIP软件的并行版本，本文千万量级的计算结果均采用了并行计算技术。

在2007年的工作中，针对DPW II提供的翼身组合体构型和翼/身/架/舱复杂组合体构型，我们详细研究了网格密度、湍流模型对上述两种构型的总体气动特性影响，其中挂架和吊舱安装阻力的计算精度与相应的实验结果取得了较好的一致［8］。本文工作的主要目的是研究网格密度对运输机构型计算结果的影响，重点是CFD软件模拟局部修型后气动特性的改变量。利用DPW III提供的多块结构网格，本文详细研究了网格密度对两种机翼构型和两种翼身组合体构型气动特性的影响。本文采用波音公司的Edward N.Tinoco采用CFL3D软件得到计算结果和NFT(National Transonic Facility)的试验结果作对比分析。

1 机翼构型的计算网格和计算状态

本文采用的计算网格来自DPW III，该结构网格由ICEM软件生成，网格结构为多块对接网格(1-to-1)，网格分为粗网格(Coarse)、中等网格(Medium)和细网格(Fine)和极细网格(Very Fine)四种，DPW_W1构型网格的详细信息如表1所示，DPW_W2构型的网格序列与此相同，需要说明的是CFL3D软件采用的网格序列与本文不同。DPW_W2是DPW_W1的简单优化外形，二者的平面参数与厚度相同，对机翼的扭转进行了优化。两种典型运输机机翼的计算构型和表面网格分布(中等)见图1。

表1 标准网格统计列表(W1)Table1 GRID statistics for standard grids(W1)

图1 DPW机翼构型的表面网格(中等)Fig.1 Surface grids for wing configurations

本文的计算状态如下:

算例1:网格收敛性研究

·M=0.76，Re=5.0×106(基于平均气动弦长)

·α=0.5°

·全湍流计算

算例2:极曲线

·M=0.76，Re=5.0×106(基于平均气动弦长)

·α =-1°，0°，0.5°，1°，1.5°，2°，2.5°，3°

·中等网格，全湍流计算

其中平均气动弦长 C=0.1976m;参考面积 Sref=0.2903m2;压心位置，ΔX=0.154245m，ΔZ=0.00m(距翼根顶点)。

1.1 算例1计算结果与分析

本文采用TRIP软件，选择二阶精度的通量差分(FDS)类型的MUSCL差分格式和Menter的 SST两方程湍流模型［9］模拟了上述两种工况的流场，采用三重网格和并行技术加速收敛。对比计算结果采用了Tinoco采用CFL3D软件得到的结果(见DPW网站)。必须说明的是CFL3D采用的网格序列与本文并不相同(见表1)，CFL3D采用的网格要比本文采用的网格规模大，这会对计算结果产生一定的影响。

1.1.1 网格收敛性研究

图2给出了DPW_W1、DPW_W2两种机翼构型采用粗网格、中等网格、密网格和极密四套网格得到的网格收敛性研究结果，同时还给出了采用CFL3D的计算结果，气动横坐标为网格节点的(-2/3)次幂，纵坐标为阻力系数。由图2中可以看出，对于两种机翼构型，本文采用SST两方程模型均得到了网格收敛性结果。图3给出了DPW_W2构型压差阻力(CD_PR)和摩擦阻力(CD_SF)分量的网格收敛结果，可以看出随着网格规模的增加，压差阻力变化很大，而摩擦阻力变化较小，本文的计算结果与CFL3D相比较，同等网格规模下压差阻力很接近，本文的摩擦阻力略小一些。DPW_W1的结果与此类似。考虑到本文采用的网格序列与CFL3D采用的网格序列不同，本文算例2计算结果主要采用密网格(2，844，416)的计算结果与CFL3D采用中等网格(4，030，464)的计算结果相比较。

图2 DPWIII机翼构型的网格收敛性研究Fig.2 Grid convergence of DPW_W1/W2

图3 DPW_W2阻力收敛性研究Fig.3 Grid convergence of drag fraction for DPW_W2

1.1.2 压力分布的比较

图4、图5给出了DPW_W1和DPW_W2两种构型在55.1%、81.4%两个典型站位上的压力分布，其中的实线是CFL3D采用1435万网格得到的结果，点划线是TRIP软件采用284万网格得到的结果。由图中可以看出，本文采用百万量级网格得到的压力分布与CFL3D采用千万量级网格得到结果非常一致。显示压力分布对网格密度反映不敏感，或者说对于压力分布的计算，本文采用的密网格已经足够了。进一步比较图4和图5可以看到，在机翼中部和靠近翼梢站位上DPW_W2构型上表面的激波明显弱于DPW_W1构型，但比较0.5°攻角的气动特性比较可以看出(图2)，激波的减弱并没有使得该来流状态下总的阻力的减少，反而使得总阻力略有增加，根本原因在于压差阻力的略有增加，而摩阻变化不大。

图4 压力分布的比较(DPW_W1)Fig.4 Comparison of Cp distribution(DPW_W1)

图5 压力系数的比较(DPW_W2)Fig.5 Comparison of Cp distribution(DPW_W2)

1.2 算例2计算结果及分析

图6给出了DPW_W1/W2两种机翼构型的极曲线和俯仰力矩系数曲线，其中极曲线的横坐标采用了CD-Λ，其中Λ为展弦比，即总阻力减去诱导阻力，标识为CDP。TRIP软件采用的是284万网格点，CFL3D软件采用的是403万网格点，湍流模型均为SST两方程模型。从极曲线可以看出，CFL3D的结果在升力系数0.46以下，TRIP的结果在0.45以下，DPW_W2的阻力系数大于DPW_W1的阻力系数;CFL3D的结果在升力系数0.46～0.66之间，TRIP的结果在在升力系数0.45～0.66之间，DPW_W2的阻力系数小于DPW_W1的阻力系数;升力系数大于0.66以后，两种构型的阻力系数大致相当，采用两种软件得到的计算结果变化趋势与两种构型的阻力差量非常接近;两种软件得到的力矩特性有较大的差异，网格拓扑结构和计算方法细节处理的不同是导致差异的主要原因，但均反映出了以某一升力系数为分界线(TRIP软件在0.58，CFL3D软件在0.54)，两种机翼构型的俯仰力矩系数变化出现交叉。

图6 DPW-W1/W2构型气动特性比较Fig.6 Aerodynamic coefficients of DPW-W1/W2 configuration

2 翼身组合体构型的计算网格和计算状态

DPWIII提供了两种翼身组合体构型，一种是DLRF6构型，另一种是在DLR-F6_FX2B构型。其中DLRF6_FX2B构型是在F6构型的基础上，通过增加翼身接合部的整流部件得到的，目的是减少DLR-F6构型翼身接合部的分离区。DPW的网站上提供了两类多块对接网格，一类是由波音公司的ZEUS系统先进处理软件生成，另一类是由ANSYS公司的ICEM软件生成。为了与波音公司CFL3D软件的计算结果对比，我们采用了波音公司的多块对接结构网格(1-to-1)，网格分为粗网格(Coarse)、中等网格(Medium)和细网格(Fine)三种，DLR-F6构型网格的详细信息如表2所示，DLR-F6_FX2B构型的网格序列与此相同。图7给出了两种典型现代运输机翼身组合体的计算局部构型和表面网格分布的局部放大图。本文的计算状态如下:

表2 标准网格统计列表(DLR-F6)Table2 Grid statistics for standard grids(WB)

图7 两种翼身组合体构型的局部放大图Fig.7 Surface grids for wing-body configurations

算例1:网格收敛性研究

·M=0.75，Re=5.0×106(基于平均气动弦长)

·CL=0.500(±0.001)

·全湍流计算

算例2:极曲线

·M=0.75，Re=5.0×106(基于平均气动弦长)

·α =-3°，-2°，-1°，-0.5°，0°，0.5°，1°，1.5°

·中等网格，全湍流计算

其中平均气动弦长C=0.1412m;参考面积Sref=0.1453m2;压心位置，ΔX=0.5049 m，ΔZ=-0.05142m(据机头顶点)。需要注意的是，DPWII相同构型的雷诺数为300万，而DPWIII的雷诺数为500万。

2.1 算例1计算结果与分析

图8给出了DPW-F6、DPW-F6_FX2B两种翼身组合体构型采用粗、中和密三套网格得到的阻力系数、摩擦阻力系数的网格收敛性计算结果，同时还给出了采用CFL3D的计算结果。需要说明的是，CFL3D除了给出了如表2所示的三种网格的计算结果外，还给出了网格规模为1586万网格点数的计算结果。可以看到本文采用TRIP软件得到了网格收敛的计算结果。在升力系数0.5条件下，DPW-F6_FX2B的阻力系数略小于DPW-F6的计算结果，随着网格密度的增加，两种构型的阻力差量有变小的趋势。以F6构型中等网格的计算结果为例，TRIP软件得到的压差阻力比CFL3D大2个阻力单位，而摩擦阻力则小12个阻力单位，总体来说，TRIP软件的阻力计算结果比CFL3D小10个阻力单位左右。以F6构型为例，随着网格密度的增加，TRIP软件的阻力变化量在7个阻力单位左右，其中压差阻力变化量在10个阻力单位，摩擦阻力的变化量在3个阻力单位左右，这再一次说明了网格密度的变化主要影响压差阻力，而摩擦阻力对网格密度的变化不是很敏感。

图8 翼身组合体构型的网格收敛性曲线Fig.8 Grid convergence of wing-body configuration

在2008年的参考文献［10］中给出了NTF的试验结果。本文采用Richardson外推的方法，由中等网格与密网格的计算结果得到网格无关的气动特性。固定升力系数下，DLR-F6及其修型构型的Richardson外推值及相应的试验结果如表3所示。对于DLR-F6构型，试验得到的阻力系数为0.02752±0.00014，理想阻力系数为0.01915±0.00014，与相应的试验结果相比较，计算得到的阻力系数和理想阻力系数均小23个阻力单位;对于DLR-F6-FX2B构型，试验得到的阻力系数为0.02728±0.00019，理想阻力系数为 0.01851±0.00019，与相应的试验结果相比较，计算的阻力系数和理想阻力系数均小20个阻力单位。修型前后试验阻力差量为-0.00024±0.00033，计算得到的修型前后的阻力差量为0.00006，在试验精度范围之内。

图9 翼身组合体构型翼身接合部的表面流线Fig.9 Surface streamline of wing-body configuration

图9给出了 DPW-F6(左)、DPW-F6_FX2B(右)两种翼身组合体构型采用中等网格得到的翼身接合部的表面流线图。对DPW-F6修型的主要目的就是减少翼身接合部的分离区，可以看到本文的计算结果准确反映了修型前后翼身接合部分离区的变化。

2.2 算例2计算结果与分析

图10给出了TRIP软件和NFT试验给出的两种翼身组合体构型的纵向气动特性，其中极曲线的横坐标为CDP。可以看到，相同攻角下，计算得到的升力系数偏高，但计算与试验均反映出修型后相同攻角下的升力系数略有减少;相同升力系数下，计算得到的阻力值偏小，但计算与试验均反映出升力系数0.5以下，修型后阻力系数略有降低且计算与试验二者差量接近，升力系数大于0.5后，修型前后的阻力差量变小;相同升力系数下，计算得到的低头力矩偏大，但计算与试验均反映出修型后低头力矩增加且计算与试验二者差量基本一致。

图10 翼身组合体构型气动特性曲线Fig.10 Aerodynamic coefficients of wing-body configuration

表3 DLR-F6/FX2B构型计算结果外推值(M=0.75，C L=0.50)Table3 Case1 DLR-F6/FX2B data extrapolated to continuum(M=0.75，C L=0.50)

3 结论

采用TRIP软件，利用DPW III提供的多块对接网格数值模拟了DPW_W1/W2两种机翼构型和DLR-F6/F6_FX2B两种翼身组合体构型，本文主要研究了网格密度对总体气动特性和典型站位压力分布的影响，通过与CFL3D的计算结果和NTF相应的试验结果对比，得到以下一些基本结论:

(1)采用SST湍流模型，本文得到了网格收敛性的结果;网格密度的变化主要影响压差阻力，对摩擦阻力的影响较小。

(2)攻角α=0.5°时，DPW_W2构型机翼上表面的激波强度明显弱于DPW_W1构型;升力系数CL=0.5条件下，DPW-F6翼身组合体修修型后，翼身接合部的分离区明显减弱。

(3)对于机翼构型和翼身组合体构型，本文的计算结果反映了修型前后气动特性的变化量。相同升力系数下，本文计算得到的翼身组合体修型前后的升力系数、阻力系数和俯仰力矩系数变化量与试验结果相当。

［1］TINOCO E N，BOGUE D R.Progress toward CFD for full flight envelope［J］.Aeronautical Jounal，2005，109(1100):451-460.

［2］JOHNSON F T，TINOCO E N，JONG Y.Thirty years of development and application of CFD at Boeing commercial airplane SEATTLE［R］.AIAA 2003-3439.

［3］OBERKAMPF W L，TRUCANO T G.Verification and validation in computational fluid dynamics［J］.Progress in Aerospace Sciences，2002，38:209-272.

［4］LEVY D W，ZICKUHR T，VASSBERG J，et al.Summary of data from the first AIAA CFD drag prediction workshop［R］.AIAA 2002-0841.

［5］LAFLIN K R，KLAUSMEYER S M，ZICKUHR T.Summary of data from the second AIAA CFD drag prediction workshop［R］.AIAA 2004-0555.

［6］VASSBERG J C，SCLAFANIA J，MARK A D.A wing-body fairing design for the DLR-F6 model:a DPW-III case study［R］.AIAA 2005-4730.

［7］王运涛，王光学，洪俊武，等.TRIP 2.0_SOLVER的开发与应用［J］.空气动力学报，2007，25(2):163-168.(WANG Y T，WANG G X，HONG J W，et al.Development and application of TRIP 2.0_SOLVER［J］.ACTA Aerodynamica Sinica，2007，25(2):163-168.in Chinese)

［8］王运涛，王光学，张玉伦.TRIP 2.0软件的确认:DPW II复杂组合体的数值模拟［J］.航空学报，2008，29(1):34-40.(WANG Y T，WANG G X，ZHANG Y L.Validation of TRIP 2.0:numerical simulation of DPW II complex con figuration［J］.ACTA Aeronautica et Astronaautica SINICA，2008，29(1):34-40.)

［9］MENTER F R.Improved two-equation k-ω turbulence models for aerodynamic-flows［R］.NASA TM-103975，1992.

［10］VASSBERG J C，TINOCO E N，MANIM.Comparison of NTF exprrimental data with CFD predictions from the third AIAA CFD drag prediction workshop［R］.AIAA 2008-6918.