2019, Vol. 37
2. 中国空气动力研究与发展中心 计算空气动力研究所, 绵阳 621000
2. Computational Aerodynamic Institute of China Aerodynamics Research and Development Center, Mianyang 621000, China
基于Navier-Stokes方程的数值模拟技术已经被广泛地用于多种飞行器复杂流场模拟，而对数值模拟技术进行验证确认以明确下一步发展方向的研究，一直是近年来CFD研究的热点[1]。国际上已经有很多组织举办了CFD验证确认的研讨会，如欧洲计算动力研讨会[2](European Computational Aerodynamics Research Project, ECARP)、美国的阻力预测会议[3-7](Drag Prediction Workshop, DPW)、高升力预测会议[8](High Lift Prediction Workshop, HiLiftPW)等。国内CFD工作者也进行了许多关于CFD验证确认的研究工作[9-14]。
为了进一步促进国内CFD验证和确认工作的稳步开展，评估CFD当前技术状态，探索CFD的发展方向。2018年4月，中国空气动力研究与发展中心(CARDC)计算空气动力研究所联合国内相关单位组织了第一届航空CFD可信度研讨会(1st Aeronautic CFD Credibility Workshop，AeCW-1)。
大会选择的运输机标模构型(CHiNa-Transport, CHN-T1)，是由CARDC自行设计的一种翼身平立尾(Wing/Body/Tails, WBT)单通道运输机巡航构型，其外形如图 1所示[15]。
|
图 1 CHN-T1翼身平立尾构型 Fig.1 CHN-T1 wing/body/tails configuration |
该构型于2016年9月在CARDC跨声速风洞(FL-26)完成了风洞试验[16]，该风洞是一座试验段尺寸为2.4m×2.4m的引射式、半回流、暂冲式跨声速增压风洞。CHN-T1模型的风洞试验照片如图 2所示[15]。另外，该模型还在欧洲跨声速风洞(European Transonic Windturnnel, ETW)进行了高雷诺数风洞试验，Re=15.0×106。
|
图 2 CHN-T1标模在2.4米风洞中试验照片 Fig.2 Photo of CHN-T1 standard model tested in 2.4 m wind tunnel |
大会提供的计算模型与风洞试验模型尺寸相同：展长为1.5482 m, 平均气动弦长为0.1937 m, 机身长度为1.5744 m。另外，组委会还提供了带支撑装置的计算模型(Wing/Body/Tails/Support, WBTS)以及通过风洞测量机翼的变形数据构造得到的静气弹计算模型(Wing/Body/Tails/Support/Deforming, WBTSD)。
另外，会议组委会确定了本次会议的研究内容：针对标模WBT构型进行网格收敛性研究(case1)；针对标模WBTS构型进行有/无支撑影响研究(case2a、case2b)；针对标模WBTSD构型进行静气动弹性研究(case2c)；针对标模WBTS构型进行雷诺数影响研究(case3)。
本文基于由CARDC自主研发的TRIP3.0软件平台，针对大会提出的计算状态，开展了CHN-T1模型的数值模拟，主要目的是评估支撑装置、机翼静气动弹性变形和雷诺数效应对CHN-T1模型气动特性数值模拟结果的影响。
1 计算方法与计算网格 1.1 计算方法本文的研究工作基于CARDC研发的亚跨超CFD软件平台(TRIP 3.0)开展。TRIP软件经过多年的应用发展[17-20]，已经成为一个非常成熟的高精度数值模拟软件平台。
本文的研究采用有限体积法和结构网格技术求解曲线坐标系下的雷诺平均Navier-Stokes方程(RANS)，控制方程的无粘项采用二阶MUSCL[22]差分格式进行离散，并利用Roe格式进行通量分裂；粘性项采用二阶中心格式进行离散；离散方程组的求解采用LU-SGS[23]方法。计算中采用SST两方程湍流模型进行全湍流模拟。另外，采用多重网格技术和大规模并行计算技术加速收敛。
1.2 计算网格本文研究中所使用的网格均为采用商业软件生成的多块结构网格，且提交给会议组委会作为大会的结构网格标准网格。进行网格生成时，单个网格块上三个维度各维度网格点数均为4N+1，以确保能够进行三重网格计算加速收敛。为确保网格能有效的捕捉流场细节，表面网格分布按照网格生成规范[24]进行：翼型前/后缘网格流向尺度约为当地弦长的0.1%，机翼翼根和翼尖约为半展长的0.1%，附面层内网格增长率以及流线的拉伸率应小于1.25，附面层第一层网格高度满足y+ < 1。空间网格拓扑采用H型，附面层网格拓扑采用O型，而在翼身结合部和平立尾与机身的结合部位附面层网格采用H型，以防止出现数值模拟流动分离区域偏小的现象[25]。
针对会议组委会确定的研究内容，分别生成了对应的网格：网格收敛性研究内容(case1)，首先对CHN-T1标模WBT构型生成了一套中等规模(Medium)多块对接网格作为基准网格，并以此为基础进行网格重构，从而分别获得粗网格(Coarse)和细网格(Fine)。(另外，本文还额外生成了一套百亿级的极细(Extremely fine)计算网格进行模拟，来进一步验证网格收敛性)。其网格拓扑及网格表面分布如图 3所示；有/无支撑影响研究内容(case2a、case2b)，选择WBTS构型在WBT构型中等网格基础上，生成了中等规模的带支撑网格，其网格拓扑及网格表面分布如图 4所示；静气动弹性影响研究内容(case2c)，选择WBTSD构型在case2b网格基础上，重新贴体生成了静气动弹性变形后的网格case2c，图 5给出了静气弹网格与未变形网格在机翼上的差异；雷诺数影响研究内容(case3)，在WBTS构型网格基础上，更改附面层内第一层距离与网格点数，得到高雷诺数下的网格case3。不同网格的具体参数见表 1。
|
图 3 CHN-T1构型case1网格拓扑及表面网格 Fig.3 Surface mesh and topology of CHN-T1 configuration for case1 |
|
图 4 CHN-T1构型case2b网格拓扑及表面网格 Fig.4 Surface mesh and topology of CHN-T1 configuration for case2b |
|
图 5 CHN-T1构型case2c表面网格 Fig.5 Surface mesh of CHN-T1 configuration for case2c |
| 表 1 CHN-T1构型不同网格参数 Table 1 Grid parameters for CHN-T1 configuration |
|
为了研究定升力系数下网格收敛性，采用标模WBT构型的粗、中、细3套网格以及百亿极细网格，为了便于表述分别记为case1-c、case1-m、case1-f，case1-e。四套网格的具体参数上一节已经给出。来流条件为：Ma=0.78，CL = 0.500(±0.001)，基于平均气动力弦长(MAC=0.1937 m)的雷诺数为Re=3.3×106。
图 6给出采用粗、中、细三套网格以及百亿极细网格获得的阻力系数收敛历程，其横坐标为迭代步数(Iter)。从图中可以看出，每套网格均可得到收敛的气动力结果。其中三套网格在5000步左右阻力系数已经接近收敛，收敛性良好(至10000步时阻力系数无变化)，且随着网格的加密，阻力系数单调减小；百亿极细网格在12000步左右阻力系数接近收敛，收敛性良好(至70000步阻力系数无变化)，且阻力系数收敛结果几乎与细网格的值重合。可以看出，阻力系数满足网格收敛性。
|
图 6 不同网格密度下的阻力系数收敛历程 Fig.6 Convergence history of drag coefficient |
表 2给出了不同网格密度下得到的迎角α、阻力系数CD、升力系数CL、压差阻力系数CD, p、摩擦阻力系数CD, f、理想阻力系数CD, idea=CD-CL2/(ARπ) 、和俯仰力矩系数Cm。表中同时给出了数值结果的插值解, 其值采用Richardson插值(Richardson Extrapolation, R.E.)方法[23, 25]基于中等网格和细网格的结果得到。图 7给出了不同网格密度下获得的阻力系数、定升迎角、俯仰力矩系数随网格密度的变化曲线，其中横坐标为网格单元数的-2/3次幂。
| 表 2 CHN-T1构型case1气动特性 Table 2 Aerodynamic characters of CHN-T1 configuration for case1 |
|
|
|
图 7 Case1气动力系数的网格收敛性 Fig.7 Grid convergence for case1 |
从图表中结果可以看出，对于粗中细三套网格，随着网格加密，阻力系数及其摩阻分量与压阻分量、定升迎角均单调变化(定升迎角增大，阻力系数减小，摩阻增大，压差阻力减小，力矩系数减小)，结果体现出良好的收敛性。因此，后续的研究内容均采用中等网格规模进行模拟。
百亿极细网格的结果除压差阻力分量外，其余结果均较插值结果略偏大，不过其阻力系数及定升迎角等结果均满足上述收敛性，只有力矩系数结果偏离了前面的单调性。由于百亿极细网格力矩系数的偏离幅度量值极小(较细网格力矩系数结果偏大0.00042)，而影响力矩系数变化的因素较多也十分复杂，而且此类超大规模网格的数值模拟也缺乏足够的数据对照，因此百亿极细网格力矩系数未保持单调变化的原因有待进一步研究。
3 CHN-T1标模支撑装置影响研究CHN-T1标模在进行风洞试验时，需在机身后体安装支撑装置以支持测量天平的数据测量。其外形如图 2所示。支撑装置的存在，会影响平尾附近的流动。直接利用标模WBT构型进行模拟与风洞试验的结果会存在一定程度的差异。需要进行标模WBTS构型的模拟以减小与风洞试验中的模型差异。
WBTS构型中的支撑装置与风洞试验中的支撑装置一致。为了研究支撑装置对模型气动特性的影响，选择网格case1-m与case2b进行模拟，来流条件为：Ma=0.78，α=-2°~4.5°，基于平均气动力弦长(MAC=0.1937 m)的雷诺数为Re=3.3×106。
图 8给出了CHN-T1标模有/无支撑构型气动特性随迎角变化曲线。从图中可以看出，有/无支撑构型的气动力特性差异不大，同迎角下升力系数有支撑构型较无支撑构型增大约0.015，而无支撑构型的升力与试验结果更接近；在零度迎角后随迎角增大，有支撑构型的阻力特性较无支撑构型增大幅度逐渐加大，在当前的计算范围内，阻力系数的变化约在2%，且有支撑构型的阻力与试验结果更接近。有支撑构型力矩系数较无支撑构型明显降低，同迎角下力矩系数大约减小0.05，且有支撑构型力矩系数更接近试验结果。
|
图 8 case2a和case2b气动力系数 Fig.8 Aerodynamic characters for case2a and case2b |
由于支撑安装在后机身平尾下方，直观感受可以知道支撑的存在对平尾流场具有很大的影响。为了更细致地研究支撑装置对全机的具体影响部件，本文进行了分部件积分。选择α=3°的结果进行比较分析。
首先给出了各个部件有/无支撑下的气动特性结果，以及各个部件的气动特性变化量在整体气动特性中的比例，如表 3。表中结果可以看出，平尾部件所占的比重最大：阻力增大4.03%，升力增大1.96%，俯仰力矩减小126.8%。即添加支撑以后，平尾的气动特性受到的影响最大。另外，从表中的数值也可以知道，添加支撑以后总体阻力增大1.91%，升力增大2.23%，俯仰力矩减小131.7%。即添加支撑对力矩特性的影响更大。
| 表 3 Case2a和case2b，α=3°各部件的气动特性 Table 3 Aerodynamic characters of each part at α=3° for case2a and case2b |
|
|
另外，图 9给出了有/无支撑条件下，平尾展向截面上/下表面压力系数分布曲线。从图中可以看出，添加支撑以后，平尾各截面上的压力曲线包围的面积均减小，即平尾上的向下的压力减小(也可认为是升力增大)。由于平尾距离飞机重心的力臂较大，平尾上的压力变化引起的力矩变化更明显，因此造成抬头力矩减小幅度较大。图 10给出了有/无支撑条件下，平尾上表面压力系数云图。从图中可以看出，添加支撑以后平尾上表面压力减小，平尾下表面压力增大，因此平尾上的升力增大。与图 9所描述的状态相同，平尾上产生一个较大的低头力矩，使得飞机的俯仰力矩减小。
|
图 9 case2a和case2b平尾上各截面表面压力系数分布曲线Cp分布 Fig.9 Cp distribution of the horizonal tail at each section for case2a and case2b |
|
图 10 case2a和case2b平尾上Cp云图 Fig.10 Cp contour of horizonal tail for case2a and case2b |
CHN-T1标模在进行风洞试验时，受气动压力的影响，随着迎角的增加，在机翼上会发生越来越大的弹性变形。机翼形变以后，其气动特性也会发生相应的变化。需要进行标模WBTSD构型的模拟，来研究静气动弹性变形的影响。
CHN-T1标模WBTSD构型中，不同迎角下机翼形变不同。为了研究静气动弹性变形的影响，选择网格case2b与case2c进行模拟，来流条件为：Ma=0.78，α=-2°~4.5°，基于平均气动力弦长(MAC=0.1937 m)的雷诺数为Re=3.3×106。
图 11给出了CHN-T1标模有/无静气弹变形气动特性随迎角变化曲线。从图中可以看出，机翼发生变形后，升阻力特性与力矩特性的变化幅度均不大。且随着迎角的增加，阻力系数与升力系数减小的幅度增大，力矩系数逐渐增大。
|
图 11 case2b和case2c气动力系数 Fig.11 Aerodynamic characters for case2b and case2c |
CHN-T1标模在进行风洞试验时，试验雷诺数(Re=3.3×106)远小于真实飞行雷诺数。为了研究雷诺数的影响，组委会还在ETW风洞进行了高雷诺数风洞试验(Re=15×106)。选择网格case2b与case3进行模拟，来流条件为：Ma=0.78，CL = 0.500(±0.001)，基于平均气动力弦长(MAC=0.1937 m)的雷诺数分别为Re=3.3×106、Re=15.0×106。
表 4给出了不同雷诺数下，定升力的气动特性结果与风洞试验结果，以及数值结果与试验结果的差量ΔS-E和不同雷诺数的数值结果差量ΔRe。从表中结果可以看出，不同雷诺数下，计算结果与试验结果均符合良好，且高雷诺数下的计算结果与试验结果的一致性更好；雷诺数增大以后，阻力系数及其摩阻分量和压阻分量均减小，定升迎角减小，力矩系数增大，其中摩阻分量相比压阻分量减少量更大。
| 表 4 case2b和case3，气动特性 Table 4 Aerodynamic characters for case2b and case3 |
|
|
图 12给出了不同雷诺数下，CL=0.500(±0.001)，机翼与平尾展向截面上/下表面压力系数分布曲线。从图中可以看出，雷诺数增大以后，机翼上表面压力减小，特别是激波后压力减小明显，机翼下表面压力略有增加，因此机翼上升力增大；对于平尾，上表面压力增大，下表面压力减小，因此平尾上产生的向下压力增大，即升力减小，抬头力矩增大。
|
图 12 Case2b和case3机翼与平尾上各截面表面压力系数分布曲线Cp分布 Fig.12 Cp distribution of the horizonal tail and wing at each section for case2b and case3 |
本文基于TRIP3.0软件平台，采用多块结构网格技术，通过求解曲线坐标系下的RANS方程，模拟了第一届航空CFD可信度研讨会提供的运输机标模构型。通过理论分析和试验结果对比，基本结论如下：
(1) 采用不同网格密度得到的结果，体现出良好的网格收敛性，百亿级的极细网格得到的结果依然符合网格收敛性；
(2) 支撑装置主要影响平尾上的气动特性，使得升力/阻力系数增大，俯仰力矩系数减小，其中对俯仰力矩系数影响更大；机翼的静弹性变形使得升力/阻力系数减小，俯仰力矩系数增大，然而对气动力的影响量均比较小；雷诺数增大使得阻力系数/俯仰力矩系数减小，主要影响气动力；
(3) 数值模拟结果与试验结果对比显示出良好的一致性，表明对于此类运输机标模构型，TRIP软件平台的模拟结果具有良好的适用性。
(4) 考虑了支撑装置、机翼静弹性变形后，数值模拟得到的结果与试验结果仍有差异，表明在对于运输机标模的模拟上，仍有值得深入研究的问题。
致谢: 本文的相关内容是在TRIP软件开发小组的共同努力下完成的，在此对课题组成员张书俊副研究员、杨小川工程师、孙岩助理研究员、岳皓研究实习员等同志表示衷心的感谢！
| [1] |
OBERKAMPF W L, TRUCANO T G. Verification and validation in computational fluid dynamics[J]. Progress in Aerospace Sciences, 2002, 38(3): 209-272. DOI:10.1016/S0376-0421(02)00005-2 |
| [2] |
HAASE W. ECARP-European computational aerodynamics research project: validation of CFD codes and assessment of turbulence models[R]. Friedrich Vieweg & Sohn Verlagsgesellschaft mbH, Stuttgart, Germany, 1997.
|
| [3] |
LEVY D W, ZICKUHR T, VASSBERG J C, et al. Summary of data from the first AIAA CFD drag prediction workshop: AIAA 2002-0841[R]. Reston: AIAA, 2002.
|
| [4] |
LAFLIN K R, KLAUSMEYER S M, ZICKUHR T, et al. Data summary from the second AIAA computational fluid dynamics drag prediction workshop[J]. Journal of Aircraft, 2005, 42(9): 1165-1178. |
| [5] |
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(6): 781-798. |
| [6] |
VASSBERG J C, TINOCO E N, MANI M, et al. Summary of the fourth AIAA CFD drag prediction workshop: AIAA 2010-4547[R]. Reston: AIAA, 2010.
|
| [7] |
LEVY D W, LAFLIN K R, TINOCO E N, et al. Summary of data from the fifth AIAA CFD drag prediction workshop: AIAA 2013-0046[R]. Reston: AIAA, 2013.
|
| [8] |
RUMSEY C L, LONG M, STUEVER R A. Summary of the first AIAA CFD high lift prediction workshop(invited): AIAA 2011-939[R]. Reston: AIAA, 2011.
|
| [9] |
阎超, 席柯, 袁武, 等. DPW系列会议述评与思考[J]. 力学进展, 2011, 41(6): 776-784. YAN C, XI K, YUAN W, et al. Review of the drag prediction workshop series[J]. Advances in Mechanics, 2011, 41(6): 776-784. (in Chinese) |
| [10] |
张耀冰, 邓有奇, 吴晓军, 等. DLR-F6翼身组合体数值计算[J]. 力学进展, 2011, 29(6): 776-784. ZHANG Y B, DENG Y Q, WU X J, et al. Drag prediction of DLR-F6 using MFlow unstructured mesh solver[J]. Acta Aerodynamica Snica, 2011, 29(2): 163-169. (in Chinese) |
| [11] |
王运涛, 孙岩, 孟德虹, 等. 包含支撑装置和机翼变形的CRM_WB构型气动特性数值模拟[J]. 航空学报, 2017, 38(10): 121202. WANG Y T, SUN Y, MENG D H, et al. Numerical simulation of aerodynamic characteristics of CRM-WB configuration with support system and wing deformation[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(10): 121202. (in Chinese) |
| [12] |
王运涛. DPW IV~DPW VI数值模拟技术综述[J]. 航空学报, 2018, 39(4): 021836. WANG Y T. An overview of DPW IV~DPWVI numerical simulation technology[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(4): 021836. (in Chinese) |
| [13] |
王运涛, 孟德虹, 孙岩, 等. CRM-WB风洞模型高阶精度数值模拟[J]. 航空学报, 2018, 39(4): 121642. WANG Y T, MENG D H, SUN Y, et al. High-order numerical simulation of CRM-WB wind tunnel model[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(4): 121642. (in Chinese) |
| [14] |
余永刚, 周铸, 黄江涛, 等. 单通道客机气动标模CHN-T1设计[J]. 空气动力学学报, 2018, 36(3): 505-513. YU Y G, ZHOU Z, HUANG J T, et al. Aerodynamic design of a standard model CHN-T1 for single-aisle passenger aircraft[J]. Acta Aerodynamica Sinica, 2018, 36(3): 505-513. DOI:10.7638/kqdlxxb-2018.0072 (in Chinese) |
| [15] |
洪俊武, 王运涛.关于组织"第一届航空CFD可信度研讨会——CHN-T1标模CFD可信度研讨会"的通知, 附件2: 候选模型[EB/OL]. http://www.cardc,cn/eReadNews.Asp.?Channelld=2&Classld=6&ld=119, 2018-05-10/2019-01-18.
|
| [16] |
李强, 刘大伟, 许新, 等. CHN-T1标模2.4米风洞气动特性试验研究[J]. 空气动力学学报, 2019, 37(2): 337-344. LI Q, LIU D W, XU X, et al. Experimental study of aerodynamic characteristic of CHN-T1 standard model in 2.4 m transonic wind tunne[J]. Acta Aerodynamica Sinica, 2019, 37(2): 337-344. DOI:10.7638/kqdlxxb-2018.0099 (in Chinese) |
| [17] |
王运涛, 洪俊武, 孟德虹. 湍流模型对梯形翼高升力构型的影响[J]. 空气动力学学报, 2013, 31(1): 52-55. WANG Y T, HONG J W, MENG D H. The influence of turbulent models to trap wing simulation[J]. Acta Aerodynamica Sinica, 2013, 31(1): 52-55. (in Chinese) |
| [18] |
王光学, 王运涛. TRIP2.0软件湍流模块的开发与确认[J]. 航空计算技术, 2007, 37(3): 20-23. WANG G X, WANG Y T. Development and validation of turbulence models of TRIP2.0[J]. Aeronautical Computing Technique, 2007, 37(3): 20-23. DOI:10.3969/j.issn.1671-654X.2007.03.006 (in Chinese) |
| [19] |
王运涛, 张玉伦, 洪俊武, 等. TRIP2.0_SOLVER的开发与应用[J]. 空气动力学学报, 2007, 25(2): 163-168. WANG Y T, ZHANG Y L, HONG J W, et al. Development and application of TRIP2.0_SOLVER[J]. Acta Aerodynamica Sinica, 2007, 25(2): 163-168. DOI:10.3969/j.issn.0258-1825.2007.02.005 (in Chinese) |
| [20] |
王运涛, 王光学, 张玉伦. 采用TRIP2.0软件计算DLR_F6构型的阻力[J]. 空气动力学学报, 2009, 27(1): 108-113. WANG Y T, WANG G X, ZHANG Y L. Drag prediction of DLR-F6 configuration with TRIP2.0 software[J]. Acta Aerodynamica Sinica, 2009, 27(1): 108-113. DOI:10.3969/j.issn.0258-1825.2009.01.020 (in Chinese) |
| [21] |
YOON S, Jameson A. Lower-upper symmetric Gauss-Sediel method for the Euler and Navier-Stokes equation[R]. AIAA Journal, 1988, 26(9): 1025-1026.
|
| [22] |
VAN LEER B. Towards the ultimate conservation difference scheme Ⅱ, monoticity and conservation combined in a second order scheme[J]. Journal of Computational Physics, 1974, 14: 361-370. DOI:10.1016/0021-9991(74)90019-9 |
| [23] |
VASSBERG J C, TINOCO E N, Mani M, et al. Summary of the Fourth AIAA CFD Drag prediction workshop: AIAA-2010-4547[R]. Reston: AIAA, 2010.
|
| [24] |
李伟, 孟德虹, 洪俊武, 等. 网格拓扑对DLR-F6构型数值模拟影响[J]. 航空学报, 2017, 38(2): 120177. LI W, Meng D H, Hong J W, et al. Numerical study of mesh topology on DLR-F6 configuration[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 120177. (in Chinese) |
| [25] |
TINOCO E N. Drag prediction with the Zeus/CFL3D system: AIAA-2004-0552[R]. Reston: AIAA, 2004.
|