多段翼型高精度数值模拟技术研究

王运涛，孟德虹，邓小刚

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

0 引 言

随着计算机硬件技术和CFD技术本身的发展，采用基于雷诺平均NS方程（RANS）的数值模拟软件已经可以模拟真实飞行器的复杂外形及全机的复杂流场，包括二维高升力翼型和三维带增升装置的全机构型［1］。基于RANS方程预测巡航构型气动特性变化趋势和模拟全湍流附着流动的能力已经逐步得到飞行器设计工程师的认可，但高升力构型的数值模拟可信度水平依然很低。高升力构型的主要模拟难点在于：边界层转捩、激波／边界层干扰、粘性尾迹区干扰、尾迹与边界层掺混和分离流动等等。

为了提高CFD软件的数值模拟精度，研究高升力构型的流动机理，空气动力学的实验工作者和CFD工作者付出了巨大的努力［2－4］，高升力构型的数值模拟也是许多CFD可信度专题会议的主题。Rumsey等人［5］对多位国外研究者的工作进行了综述，但这些研究工作基本上全部采用二阶空间离散精度的差分格式。三阶以上空间离散精度的差分格式受到几何守恒率、对接边界的高精度处理等技术的影响，在复杂外形的应用上一直处于探索阶段，目前高精度格式的应用主要集中于简单构型的复杂流动机理研究，如气动噪声机理、超声速底部流动等。

邓小刚团队从20世纪90年代致力于高精度计算方法研究，先后构造了3至9阶的线性耗散紧致格式（Dissipative Compact Schemes，DCS）和一类高阶紧致非线性格式（Compact Nonlinear Schemes，CNS），通过引入加权思想，又构造了一系列加权紧致非线性格式 WCNS（Weighted Compact Nonlinear Scheme），包括隐式和显式两种类型，并通过典型算例验证了WCNS格式的耗散特性和色散特性。近期在几何守恒率、对接边界高精度处理等复杂外形应用技术方面取得重要进展［6－7］，为高精度格式的工程应用奠定了良好的基础。

本文基于空间五阶精度的 WCNS格式，综合利用几何守恒率和对接边界的高精度处理技术，通过求解任意坐标系下的雷诺平均的N－S方程，采用两种多段翼型的标准算例，考核了WCNS格式的模拟低速复杂构型的能力。详细研究了SA一方程湍流模型、SST两方程湍流模型、转捩位置对该翼型压力分布和典型站位速度型的影响。研究结果表明，高精度WCNS格式具有模拟低速复杂构型流场的能力，采用全湍流模拟方式可以较好地模拟多段翼型的压力分布，但对边界层速度型和阻力的模拟精度较差；模拟实验的转捩位置可以显著提高附面层内速度型的模拟精度和阻力的预测精度。

1 高精度计算方法

控制方程为任意坐标系的雷诺平均的NS方程：

其中Q为守恒变量，E、F、G为对流通量，Ev、Fv、Gv为粘性通量。对流项离散格式采用5阶精度加权紧致非线性格式（WCNS），粘性项离散格式为4阶精度中心格式，对流项离散及边界与近边界的具体格式为：

（a）内点格式

对流项离散采用原始变量型的 WCNS－E－5格式，设网格间距为h，以ξ方向为例

（b）边界格式和近边界格式

近期的研究发现，几何守恒率与对接边界的处理是制约高精度格式应用于复杂外形的关键因素，这方面的最新研究进展可参考文献［6］和［7］，这里不再赘述。

表1 NLR7301两段翼型的气动特性（GA P2．6%）Table 1 Aerodynamic character of NLR7301 two element airfoil

2 NLR7301两段翼型的高精度数值模拟

NLR－7301两段翼型的风洞实验是在70年代末期在NLR3m×2m低速风洞中完成的，实验结果包括了总体气动特性、压力分布、典型站位的速度型分布等多种数据［8］。该翼型的襟翼偏角为20°主翼／襟翼重叠区均为5．3%c；实验构型包含了两种不同的缝隙宽度，一种为2．6%c，另一种为1．3%c。该翼型是典型的起飞构型，主翼的后缘经过仔细的修型，全流场没有明显的分离。本文采用了缝隙宽度为2．6%c的构型，计算构型、典型测量站位见图1。

图1 NLR7301构型及速度型站位Fig．1 NLR7301profile and velocity locations

本文采用多块对接网格技术模拟该两段翼型，网格单元数为20万，物面第一层距离7．5E－6，远场150c。，局部网格分布见图2。

图2 NLR7301的计算网格（局部放大）Fig．2 NLR7301computational grid（local）

·湍流模型的比较

采用WCNS格式、Spalart－Allmaras一方程湍流模型和Menter's k－ωSST两方程模型，全湍流方式数值模拟了该两段翼型的绕流流场。来流马赫数0．185，攻角为6°，雷诺数2．51×106。主翼与襟翼的压力分布与实验结果的比较见图3，气动特性的数值模拟结果见表1。由图中看出，采用两种湍流模型得到压力分布基本一致，且与实验结果吻合良好。但从总体气动特性上来看，采用两种湍流模型得到的阻力系数明显比实验值偏高，其中，采用SA模型阻力系数大43个阻力单位（1个阻力单位＝0．0001），采用SST模型阻力系数大19个阻力单位。说明采用高精度格式和全湍流方式，可以较好地模拟压力分布，但阻力系数明显偏大。

图3 攻角为6°时压力分布比较图Fig．3 Surface pressure coefficient（fully turbulent）

图4 给出了采用两种湍流模型得到的典型站位的速度型与实验的比较。由图中看出，两种湍流模型得到速度型模拟结果基本一致，但在边界层厚度、主翼尾迹区的模拟精度与实验结果差别较大，这是导致计算得到的阻力系数与实验结果差别较大的主要原因。

图4 湍流模型对典型站位速度型的影响Fig．4 Effect of turbulence model to velocity profiles

·转捩位置的影响

NLR－7301的边界层测量的实验结果给出了典型迎角下，机翼上翼面、下翼面和襟翼上表面的转捩位置［8］，攻角6°的转捩位置见表2。需要说明的是，表2中给出的主翼上表面的转捩位置实际上是主翼前缘的层流分离到湍流再附的起始点和终止点，数值模拟的过程中将主翼上表面层流分离的起始点作为转捩模拟的起始点。本文采用SA一方程模型，直接模拟实验的转捩位置，重点研究了转捩位置对阻力特性、边界层速度型的影响。

表2 NLR－7301两段翼型转捩位置（GA P2．6%）Table 2 Transition location of NLR 7301two element airfoil

图5给出了是否模拟转捩位置对边界层内速度型的影响，包括了主翼上表面、下表面及襟翼上表面的典型站位。与采用全湍流方式相比较，模拟了实验的转捩位置后，主翼上下表面的边界层厚度，尤其是主翼下表面的速度型改善尤为明显，襟翼上表面的边界层厚度和主翼的尾迹区的模拟结果也有明显的改善。模拟转捩位置后的速度型分布与实验结果取得了较好的一致。边界层内速度型模拟精度的提高大大改善了阻力的预测精度，与实验结果几乎相同（见表1）。

图5 转捩位置对典型站位速度型的影响Fig．5 Effect of transition location to velocity profiles

3 30P－30N三段翼型的高精度数值模拟

麦道航空公司的30P－30N三段增升构型是被CFD工作者广泛采用的多段构型之一，该翼型的前缘缝翼和后缘襟翼的偏角均为30°，前缘缝翼的缝道宽度为2．95%，外伸量为－2．5%；后缘襟翼缝道宽度1．27%，外伸量为0．25%，是典型的着陆构型。计算构型、速度型测量位置如图6所示。

本文的数值模拟采用了多块对接结构网格技术，“C”型网格为主导，计算区域的远场边界取为120倍弦长。壁面的第一排网格达到了1．0×10－5弦长，网格在各个剪切层附近均进行了适当的加密，以保证附面层内和剪切层的数值模拟精度，网格规模达到了13．8万。图7给出了计算网格的局部放大图。

图6 30P－30N构型及速度型站位Fig．6 30P－30Nprofile and velocity locations

图7 30P－30N的计算网格（局部放大）Fig．7 30P－30Ncomputational grid（local）

·湍流模型的比较

采用 WCNS格式、Spalart－Allmaras一方程湍流模型和Menter's k－ωSST两方程模型，全湍流方式数值模拟了该三段翼型的绕流流场。本文模拟的状态为：马赫数M＝0．2，攻角α＝19°，雷诺数Re＝9×106。该攻角已接近失速攻角，这使得该状态的流动模拟更具挑战性。计算得到的压力分布与实验结果的比较见图8，气动特性的数值模拟结果见表3。由图中看出，采用两种湍流模型得到压力分布基本一致，且与实验结果吻合良好。但从气动特性上来看，采用两种湍流模型得到的阻力系数相差67个阻力单位，由于没有相应的阻力实验结果，这里没办法进行定量上的比较。

图8 攻角为19°时压力分布比较图Fig．8 Surface pressure coefficient（fully turbulent）

表3 30P－30N三段翼型的气动特性Table 3 Aerodynamic character of 30P－30Nthree element airfoil

图9给出了采用两种湍流模型得到的典型站位的速度型与实验的比较。由图中看出，两种湍流模型得到速度型模拟结果基本一致，但在缝翼尾迹区、边界层厚度、主翼尾迹区的模拟精度与实验结果差别较大。

图9 湍流模型对速度型的影响（全湍流）Fig．9 Effect of turbulence model to velocity profiles

·转捩位置的影响

文献［5］中给出了相对于收缩构型弦长c的转捩区域（见表4）。其中实验测量值中“n／a”表示在缝翼下表面型面尖点以后、主翼下表面型面尖点以后和襟翼下表面均没有观察到确定的转捩位置。为了研究转捩位置对速度型模拟精度的影响，本文选用SA一方程模型数值模拟了实验的转捩位置，对应于实验值没有明确转捩位置的点，本文分别选择了缝翼下表面尖点前、主翼下表面尖点前一小段距离为转捩区域，而在襟翼下表面则为层流区域。

表4 计算与实验的转捩位置Table 4 Transition location for calculation and test

图10给出了采用SA一方程模型是否模拟实验转捩位置的典型站位速度型的比较。可以看到模拟了实验的转捩位置后，典型站位的速度型得到了明显的改善。首先是主翼上边界层的厚度和与实验结果更加接近，其次是缝翼和主翼的尾迹区速度值和尾迹区的宽度均有明显改善，但与实验结果相比较还存在一定差距。文献［5］中指出，推迟缝翼上表面的转捩位置可以进一步提高缝翼尾迹区的数值模拟精度。从表3中可以看出，模拟了实验的转捩位置后，阻力系数降低了78个阻力单位。

图10 典型站位速度型的比较（转捩位置A）Fig．10 Velocity profile on typical station（Transition location A）

4 结 论

本文采用五阶精度的WCNS格式和结构对接网格技术，综合应用几何守恒率和对接边界高精度处理技术，通过求解任意坐标系下的RANS方程，数值模拟了NLR7301二段翼型和30P－30N三段翼型的复杂流场，主要研究了SA和SST两种湍流模型、转捩位置对数值模拟结果的影响。通过与相应的实验结果相比较，得到以下一些基本结论：（1）采用几何守恒率和对接边界高精度处理技术，WCNS格式可以模拟多段翼型的复杂流场。（2）采用全湍流模拟方式，两种湍流模型均可以较好地模拟多段翼型的压力分布，但速度型的计算结果与实验结果差别明显。（3）模拟实验的转捩位置，对气动力系数，尤其是阻力系数影响显著，计算结果更接近实验值；可以明显提高边界层、主翼尾迹区的数值模拟精度，但缝翼尾迹区的计算精度仍然与实验结果存在明显差距，需要进一步研究。

致谢：本文的研究工作是基于WCNS高精度数值模拟软件完成的。在此，对软件的开发者表示感谢。

［1］ TINOCO E N，BOGUE D R．Progress toward CFD for full flight envelope［J］．Aeronautical Jounal，2005，109：451－460．

［2］ PAUL L JOHNSON，KENNETH M JONES，MICHAEL D MADSON．Experimental investigation of a simplified 3Dhigh lift configuration in support of CFD validation［R］．AIAA 2000－4217．

［3］ STUART E ROGERS，KARLIN ROTH，STEVEN M NASH．CFD validation of high－lift flows with significant wind－tunnel effects［R］．AIAA 2000－4218．

［4］ 朱自强，陈迎春，吴宗成，等．高升力系统外形的数值模拟计算［J］．航空学报，2005，26（3）：257－262．

［5］ CHRISTOPHER L．RUMSEY，SUSAN X．YING．Prediction of high lift：review of present CFD capability［J］．Progress in Aerospace Sciences，2002，38：145－180．

［6］ XIAOGANG DENG，MEILIANG MAO，GUOHUA TU et al．Extending weighted compact nonlinear schemes to complex grids with character－based interface condition［J］．AIAA Journal，2010，48（12）：2840－2851．

［7］ XIAOGANG DENG，MEILIANG MAO，GUOHUA TU．Geometric conservation law and applications to high－order finite difference schemes with stationary grids［J］．Journal of Computational Physics，2010，230（4）：1100－1115．

［8］ VAN DEN BERG B．Boundary layer measurements on a two－dimensional wing with flap［R］．NLR TR 79009U，1979．