HiLiftPW-3高升力构型数值模拟

洪俊武，王运涛，李伟,2, *，杨小川

1. 中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000 2. 国防科技大学 空天科学学院，长沙 410073

目前，尽管计算流体力学(Computational Fluid Dynamics,CFD)方法和计算机技术都获得了长足进步与发展，但基于雷诺平均Navier-Stokes(Reynolds Averaged Navior-Stokes,RANS)方程模拟复杂高升力构型的可信度依然有待提高，采用主流CFD方法获得高升力构型的最大升力系数及失速迎角，并准确预测和模拟分离流动的产生和发展，依然是CFD面临的巨大挑战[1-2]。为了对高升力构型的流动机理进行研究并深入展开CFD软件对此类构型的验证与确认工作，许多CFD可信度专题会议也以此作为研究主题[3-7]，其中，最具代表性的会议就是AIAA高升力构型性能预测会议。

2017年6月，第3届AIAA高升力构型(HiLiftPW-3)预测会议在美国丹佛召开。HiLiftPW-3会议选择的2组标模分别是NASA(National Aeronautics and Space Administration)提供的HL-CRM(HiLift-Common Research Model)标模和JAXA(Japan Aerospace eXploration Agency)提供的JSM(Jaxa hilift configuration Standard Model)标模。会议组委会确定的研究工况分别是：针对HL-CRM标模，开展网格收敛性研究(Case1)；针对JSM标模，开展气动力特性对比研究(Case2)；针对DSMA661二维翼型，开展湍流模型的标定(Case3)。

Case1的主要内容是对HL-CRM标模进行典型状态的网格收敛性研究，组委会没有提供计算状态的风洞试验结果，目的是对各CFD软件开展验证工作。Case2的主要内容是对JSM标模有/无挂架短舱的2种构型进行模拟，并将计算结果与风洞试验数据进行对比，目的是考察挂架短舱对气动特性的影响，确认CFD结果的模拟精度。Case3的主要内容是研究DSMA661翼型压力分布、速度型分布和尾迹区流动，主要目的是验证所采用湍流模型的有效性。本文重点讨论Case1和Case2的数值模拟结果。

采用TRIP3.0软件，王运涛等[8-10]研究了网格规模、湍流模型、转捩模型等多种因素对TrapWing全展长襟翼高升力构型数值模拟结果的影响，洪俊武等[11]采用拼接结构网格技术模拟了TrapWing半展长襟翼高升力构型的气动特性和压力分布特性。采用五阶空间离散精度的WCNS(Weighted Compact Nonlinear Scheme)格式，王运涛[12-13]和李松[14]等，开展了典型运输机高升力构型数值模拟技术研究。此外，国内一些CFD工作者[15-20]采用结构网格技术、非结构网格技术对高升力构型开展了一系列数值模拟技术研究工作。

本文作者作为HiLiftPW-3的参研团队之一(参会序号为015)，利用自主开发的CFD软件平台TRIP3.0，根据HiLiftPW-3两种高升力构型标模的外形和流场特征，采用拼接结构网格技术和有限体积方法，通过求解三维曲线坐标系下的RANS方程，对会议要求的Case1和Case2工况进行了数值模拟。研究结果表明，本文方法对典型运输机三段翼布局的低速问题模拟具有良好的适用性，在失速迎角之前，数值模拟得到的气动力系数和压力分布均与试验结果高度吻合，同时也取得了合理的网格收敛性结果。

1 计算方法

本项研究采用有限体积法和结构网格技术求解曲线坐标系下的RANS方程。离散方程组的求解采用LU-SGS[21]方法；无黏项离散采用MUSCL格式[22]，并利用Roe格式进行通量分裂；黏性项采用中心差分格式进行计算。计算中，湍流模型采用目前应用较为广泛的SA(Spalart-Almaras)一方程模型[23]，该模型计算量相对较小，而且具有良好的鲁棒性。

采用了多重网格技术、低速预处理技术、并行技术加速收敛。并行计算采用域分解法进行数据划分，并使用基于消息传递的并行编程模型来设计开发并行计算模块，迭代过程中采用MPI消息传递库，通过显式地发送和接收消息来完成各处理器之间的数据交换，具有较好的通用性和可移植性，可以获得较高的并行效率。实际计算中针对不同算例，在国产飞腾CPU上采用了16～640个CPU进行并行运算。

采用了拼接网格技术(Surface to Surface)来处理前缘缝翼和后缘襟翼偏转带来的剪刀缝，针对此类拓扑结构，拼接网格技术可以有效提高网格质量并大大降低网格规模，同时能在拼接面上保证通量守恒，本文的拼接网格方法在此前的大量应用中也已经得到了充分验证。

2 HL-CRM标模网格收敛性

HL-CRM是NASA提供的高升力构型标模，该构型是翼身组合体布局，机翼采用三段构型，多段翼之间无固定连接装置。HL-CRM构型的前缘缝翼与后缘襟翼的偏角分别为30°和37°，前缘缝翼的展向长度基本相当于机翼展长，从翼根延伸到翼尖；后缘襟翼的展向长度大约是机翼展长的2/3，从翼根延伸到展向外侧，该构型为典型的着陆构型。

2.1 计算网格

本文在研究中首先采用商业软件生成了一套基准的多块拼接网格,并在这套基准的中等规模(Medium)网格基础上，对网格进行重构，获得粗化的Tiny、Coarse网格和细化的Fine网格。这4套不同密度网格的单元规模分别为432万，1 458万、3 455万和11 661万，各套网格的具体参数见表1，中等网格表面网格及截面网格见图1。

图1 HL-CRM标横网格Fig.1 Grid for HL-CRM model

表1 HL-CRM标模网格参数Table 1 Grid parameters of HL-CRM model

2.2 网格收敛性

本节采用前述的4套计算网格开展了网格收敛性研究，来流条件为：马赫数Ma=0.2，迎角α=8°、16°，基于平均气动力弦长(MAC=7 005.32 mm)的雷诺数为Re=3.26×106。

图2给出了极粗、粗、中、细4套不同网格密度下升力系数CL和平均残差RESave的迭代收敛历程，其中横坐标为迭代步数(Iteration)，所有网格均采用三重网格进行计算。可以看到，4套网格的计算结果均获得了良好的收敛性，极粗网格(Tiny)收敛速度最快，细网格(Fine)收敛最慢；极粗网格约3 000步残差降到最低，粗网格需要5 000 步，中等网格需要接近10 000步，而细网格需要20 000步左右。由于网格越密，近壁面区的网格尺度也越小，导致附面层网格单元的当地时间步长也越小，这就使得流场收敛速度大大降低。

图2 不同网格密度下的升力系数和平均残差收敛历程Fig.2 Convergence history of lift coefficient and average residual under different grid densities

图3给出了采用不同密度网格得到的HL-CRM气动特性随网格单元数量的变化曲线，图中横坐标中N是网格单元数目。可以看到，不同迎角下升力系数(CL)、阻力系数(CD)和俯仰力矩系数(Cm)随网格密度的增加都是单调变化的，计算结果具有良好的网格收敛性。另外，从图3中可以看到，α=16°时粗网格的升力系数和力矩系数都和密网格结果相差较大，而α=8°时粗、密网格的升力系数和力矩系数差量则小得多。这是由于小迎角下，机翼周围流场变化相对较小，此时的粗网格基本能捕捉到该状态下空间流场的主要特征，所以粗网格的计算结果和密网格结果差距不大；而大迎角下，机翼周围流场变化较为剧烈，一些区域甚至出现了小的分离气泡，此时的粗网格由于尺度过大已经无法再现真实流场的主要特征，因此和密网格结果差距明显。

图3 Case1气动力系数的网格收敛性曲线Fig.3 Grid convergence curves of aerodynamic coefficients for Case1

3 JSM标模有/无挂架短舱对比

JSM是JAXA提供的高升力构型标模，该构型也是翼身组合体布局，机翼采用三段翼构型，前缘缝翼通过8个小连接块与主翼固定，后缘襟翼通过3个大的滑轨与主翼固定。另外，该构型分为无发动机挂架短舱的Case2a和有发动机挂架短舱的Case2c，本次会议主要研究该构型在有/无短舱情况下气动力特性与试验的对比以及有/无短舱情况下二者的气动特性差异。该构型的前缘缝翼与后缘襟翼的偏角分别是为25°和30°，前缘缝翼的展向长度基本相当于机翼展长，从翼根延伸到翼尖；后缘襟翼的展向长度大约是机翼展长的3/4，从翼根延伸到展向外侧，该构型也是典型的着陆构型。

3.1 计算网格

Case2计算网格的拓扑结构与Case1类似，Case2a(无短舱)和Case2c(有短舱)2套网格的具体参数见表2。图4给出了Case2a构型在风洞中的安装照片和计算网格，图5给出了Case2c构型在风洞中的安装照片和计算网格。由于Case2包含滑轨和多个连接部件，构型复杂程度比Case1高得多，因此网格量也大得多。

表2 JSM标模网格参数Table 2 Grid parameters of JSM model

图4 Case2a的风洞安装照片和计算网格Fig.4 Installation picture in wind tunnel and computational grid for Case2a

图5 Case2c的风洞安装照片和计算网格Fig.5 Installation picture in wind tunnel and computational grid for Case2c

3.2 Case2a无挂架短舱构型计算

Case2的来流条件为：Ma=0.172，α=4.36°,10.47°,14.54°,18.58°,20.59°,21.57°，基于平均气动力弦长(MAC=529.2 mm)的雷诺数为Re=1.93×106。图6给出了JSM无挂架短舱(Nacelle/Pylon off)构型气动特性随迎角变化曲线。

从升力曲线来看，计算和试验的失速迎角十分吻合，都在20.59°附近；在19.59°之前，二者几乎重合；失速迎角附近，计算得到的升力系数略高。作者认为，这个差异是试验模型的弹性变形造成的。由于试验模型在大迎角时翼尖上翘明显，造成升力系数偏低；而计算模型采用的是刚性外形，这使得计算值相对试验结果偏高。

从阻力曲线来看，计算结果比试验值明显偏高，一方面本文采用了SA模型进行全湍流模拟，而SA的数值模拟结果通常会导致阻力略微偏大；另一方面这也是本届高升力会议出现的一个普遍现象，图7是本届高升力会议组织方发布的Case2两组升阻比曲线统计结果，可以看到，除个别明显异常的结果外，30余家科研机构提供的40余组升阻比曲线全部位于试验曲线的右侧，阻力普遍偏大，所有数值结果均存在不同程度的平移。与这一现象不同的是，第1届高升力构型会议中主流软件的阻力计算结果均与试验值高度吻合[7]；第2届会议中，绝大多数软件获得的阻力计算结果也明显超过试验值，而且迎角越大差异越明显[24]。虽然高升力构型的研究重点并不在于阻力的精确模拟，但是从数据对比来看，这个问题也是CFD工作者需要密切关注的重点问题之一。

从力矩曲线来看，计算结果也与试验值吻合良好。产生这一现象的主要原因是该构型没有平尾，所以全机力矩的主要来源变成了机翼，这也意味着只要计算得到的机翼压力分布与试验吻合的话，升力和力矩系数都应该是吻合的。下文的压力分布对比也证明了这一点。

图6 Case2a气动力系数计算与试验结果的比较Fig.6 Comparison of aerodynamic coefficients of calculation and test results for Case2a

图8为组委会给出的Case2a所有参会结果与试验结果的气动特性对比，图中黑色为试验结果，红色为参会的全部数值模拟结果，绿色曲线为本文结果。可以看到，本文结果具有良好的可信度。

图9给出了α=10.47°时，典型机翼站位的压力系数(Cp)分布计算值与试验结果对比，4个展向站位分别是η=16%(A-A)，41%(D-D)，56%(E-E)和89%(H-H)。可以看到，该迎角下，不同站位前缘缝翼、主翼和后缘襟翼上压力分布的试验和计算值高度吻合，除个别试验跳点以外，几乎重合在一起。

图7 Case2升力随阻力变化曲线统计结果与试验结果的比较Fig.7 Comparison of lift-drag curves between statistical and test results for Case2

图8 Case2a所有参会结果与试验的比较Fig.8 Comparison of results for Case2a including all simulations and test

图9 Case2a α=10.47°典型站位Cp分布对比Fig.9 Comparison of distribution of Cp at typical station α=10.47° for Case2a

图10给出了α=18.58°典型站位的Cp分布计算值与试验结果对比，可以看到，该迎角下位于翼根和机翼中部3个不同站位的试验和计算值高度吻合，几乎重合；而接近翼尖的89%站位试验值与计算值差异明显。产生这一现象的原因前面已经阐述了，这是由于风洞试验使用的是弹性模型，数值计算模拟的是刚性模型，试验的弹性模型在大迎角下翼尖变形量较大，因此Cp差异较为明显；而小迎角下的所有站位和大迎角下的非翼尖站位变形量都较小，Cp分布都吻合良好。

图10 Case2a α=18.58°典型站位Cp分布对比Fig.10 Comparison of distribution of Cp at typical station α=18.58° for Case2a

3.3 Case2c有挂架短舱构型计算

图11给出了Case2c有挂架短舱(Nacelle/Pylon on)构型气动特性随迎角变化曲线。从升力曲线来看，计算和试验的失速迎角较为吻合，试验结果在20.09°附近，计算结果略为提前，失速迎角为19.59°。在19.59°之前，二者几乎重合。从阻力曲线来看，计算结果也比试验值偏高，与Case2a规律一致；从力矩曲线来看，在失速迎角之前，计算结果也与试验值吻合良好。图12为组委会给出的Case2c所有参会计算结果与试验结果的气动特性对比，图中绿色曲线为本文结果。可以看到，本文结果具有良好的可信度。

图13给出了α=10.47°典型机翼站位的Cp分布计算值与试验结果对比。可以看到，该迎角下不同站位的试验和计算值高度吻合，几乎重合在一起。图14给出了α=18.58°典型机翼站位的Cp分布计算值与试验结果对比，可以看到，该迎角下位于翼根和机翼中部3个不同站位的试验和计算值也高度吻合，而接近翼尖的站位试验值与计算值差异明显，产生这一现象的原因与Case2a完全相同，此处不再赘述。

图11 Case2c气动力系数计算与试验结果的比较Fig.11 Comparison of aerodynamic coefficients of calculation and test results for Case2c

图12 Case2c所有参会结果与试验结果的比较Fig.12 Comparison of results Case2c including all simulations and test

图13 Case2c α=10.47°典型站位Cp分布对比Fig.13 Comparison of distribution of Cp at typical station α=10.47° for Case2c

图14 Case2c α=18.58°典型站位Cp分布对比Fig.14 Comparison of distribution of Cp at typical station α=18.58° for Case2c

3.4 有/无挂架短舱构型的气动力差量对比

图15给出了在失速迎角之前，有/无挂架短舱2种构型气动力差异的试验和计算结果对比。为方便比较，图中给出了组委会提供的所有参会结果，绿色曲线为本文结果。

由于2种构型的气动力差量是小量，所以图中本文的升力差量曲线与试验值虽然不完全吻合，但是2条曲线的平均差量只有0.02左右，二者实际上是很接近的。可以看到，小迎角情况下，试验和计算得到的各气动力系数差量都相当吻合，大迎角时差距有所增大，但也在合理范围之内。这也进一步说明了，本文的CFD方法对于此类高升力构型的模拟是非常适用的，无论是气动力系数还是不同构型之间气动力系数的差量，计算和试验结果都具有良好的一致性。

图15 Case2两种构型气动力差量所有参会结果与试验结果的比较Fig.15 Comparison of aerodynamic dispersion between two configurations for Case2 including all simulations and test results

4 结 论

本文采用拼接结构网格技术，通过求解曲线坐标系下的RANS方程，数值模拟了第3届高升力构型会议提供的2组标模。通过理论分析和与试验结果的对比，结论如下：

1) 本次对比研究中，风洞模型翼尖在大迎角下的变形是导致试验结果和计算结果差异的主要来源之一。

2) 对此类三段翼典型布局的高升力构型，数值方法在失速迎角之前，可以获得可信度相当高的流场结果。对比结果表明，无论是压力分布还是气动力系数，试验和计算都取得了良好的一致性。

3) 本文的数值实践表明，TRIP软件平台采用的拼接结构网格技术和计算方法对此类外形和低速流场具有良好的适用性；本文的计算结果在与全球40余组数值结果的各项对比中均名列前茅，在各大科研机构和CFD软件的同台竞技中表现优秀，计算结果与试验结果合理地吻合，可以为大型飞机增升装置的气动设计提供可靠的技术支持。