龙海斌,吴裕平
(中国直升机设计研究所,江西 景德镇 333001)
飞行器标准机身模型通常用于风洞测量系统校验、机身气动特性CFD计算验证与确认等。在直升机机身风洞试验和数值模拟计算发展过程中,相继发展了ROBIN、ROBIN mod7、HELIFUSE和NUAA等标准机身模型。ROBIN(Rotor Body Interaction,简称ROBIN)机身模型的三维外形可用超椭圆方程进行描述,主要被美国NASA等机构用于旋翼/机身干扰研究、直升机机身气动特性数值模拟验证与确认等工作。1979年美国针对ROBIN机身模型进行了无旋翼状态(光机身)风洞测量试验,得到了部分表面位置的压力系数[1]。机身的长度为旋翼半径R的2倍。压力系数的测量位置包括机身开始端截面、向凸台过渡段截面、凸台开始截面、凸台中段、凸台结束位置以及尾梁截面等典型位置,如图1所示。测量过程中包含-10°、-5°、0°、5°攻角和-15°、-5°、0°、5°侧滑角状态。由于有风洞试验结果可以对比分析,而且机身外形比较简单,全球很多科研人员对ROBIN机身模型的绕流进行了数值模拟。本文对近年来ROBIN机身模型的数值模拟技术进行了梳理与分析,包括网格划分之前的网格类型的选择,网格划分过程中的数量控制,求解过程中的计算方法选取,数值模拟结果分析等[2-32]。
图1 ROBIN机身外形与表面测压孔截面位置图[33]
在机身气动特性数值模拟之前,首先要选择机身计算域网格划分类型。目前常用的网格类型有结构网格、非结构网格和混合网格等。其中结构网格的质量比较高,划分速度也比较快,在部分情况下计算结果的准确度比较高。但是,结构网格通常只适用于外形比较简单的气动特性数值模拟。非结构网格对复杂外形的适应能力比较强,在对复杂外形流体域划分网格时,非结构网格的数量相对较少。
由于ROBIN机身外形相对比较简单,因此部分研究人员在CFD计算过程中采用结构化网格,而且结构网格划分得比较细密,如图2所示。
图2 结构网格示意图[5]
在统计的29例ROBIN机身模型数值模拟过程中,有20例采用非结构网格,占总数的69%,说明大部分研究人员还是倾向于应用非结构网格进行机身气动特性数值模拟。主要原因有两个:一是ROBIN机身数值模拟通常作为数值模拟方法的验证算例,后续该数值模拟方法还要用于更复杂的机身绕流、外挂物流场甚至旋翼/机身干扰流场等计算,为了保证数值模拟方法验证的有效性,需采用一致的网格划分方法。二是随着非结构网格的发展与应用,非结构网格的划分过程比较简单,耗费的时间也比较少,而且数值模拟结果的准确度和精度越来越高,因此非结构网格的应用越来越广泛。典型的非结构四面体网格如图3所示[6]。
图3 非结构四面体网格示意图
在机身气动特性数值模拟过程中,网格数量不仅影响数值模拟结果的准确度,而且与整个数值模拟耗费的时间密切相关。同时,网格数量与网格划分的时间也紧密相关。在型号研制过程中还要考虑整个项目周期的要求等,在部分情况下还要考虑计算过程中的内存占用,因此需要限制流体域划分的网格数量。
网格划分比较稀疏时,可能造成部分机身外形失真,保形能力比较差,导致数值模拟结果有误差,同时部分流场细节难以模拟出来。而网格划分比较细密时,有可能导致部分网格的质量比较差,比如在非四面体网格划分过程中产生“金字塔”形网格,计算过程中难以收敛。网格数量多时,不仅网格划分耗费时间多,而且计算过程中耗费的时间也长,占用内存等计算资源也比较多。
在ROBIN机身数值模拟过程中,划分的网格数量大多数在百万量级,最少的为40.9万四面体非结构网格,最多的达到1064万六面体结构网格。在公开发表的19个数值模拟过程中,网格或网格点数量在各数量级的分布如表1所示。从表中可以看出网格或网格点数量在100万~1000万之间的占57.9%。部分网格如图4和图5所示。
表1 网格或网格点的数量级分布
图4 46.9万非结构网格示意图[12]
图5 1064万结构网格示意图[13]
在-5°攻角时,采用不同的网格数量数值模拟得到的凸台结束位置(x/R=1.0008)的压力系数与风洞试验结果的对比如图6所示。
图6 网格数量变化时数值模拟结果对比图
从图中可以看出,在一定数量范围内,网格数量的变化对数值模拟结果的准确度影响比较小。
在ROBIN机身模型数值模拟过程中,主要采用求解Eluer方程和求解N-S方程两种方法。在早期,部分研究人员采用面元法[17,18]求解,这种方法通过布置流动的奇点来求解气动问题。面元法的优点是计算速度非常快,网格划分简单,但是数值模拟结果的精度比较低,现在已经很少使用。求解Euler方程的方法在数值模拟过程中没有考虑流体的粘性,因此该方法计算速度相对比较快,计算过程中占用的计算资源也比较少。因此,后续对更复杂的组合机身外形、加装外挂物和旋翼/机身干扰流动等进行数值模拟时更方便和快速。由于没有考虑空气的粘性,因此对部分流动复杂区域的模拟能力很弱。目前大多数人员采用求解N-S方程的方法来对机身模型气动特性进行数值模拟。在统计的31个算例中,有19个采用求解N-S方程的方法,占61.3%。目前求解N-S方程的方法主要有雷诺平均(RANS)方法、大涡模拟(LES)方法和直接数值模拟(DNS)方法等。其中雷诺平均(RANS)方法对网格数量的要求最低,消耗的计算资源也比较少,数值模拟所花费的时间也最少,因此应用最广泛。在运用雷诺平均(RANS)方法对机身模型气动特性进行数值模拟时,需要引入湍流模型来使得方程封闭。常用的湍流模型有B-L模型、S-A模型和k-ω模型等,其中S-A湍流模型应用最多,占到2/3以上。采用各湍流模型的数量和百分比如表2所示。采用不同湍流模型得到的数值模拟与风洞试验结果如图7所示。
表2 各湍流模型应用情况表
图7 不同湍流模型模拟结果对比图[20]
部分研究人员采用尺度自适应模拟(SAS)方法对ROBIN机身模型的气动特性进行了数值模拟[22]。该方法通过引入第二长度尺度,可以实现雷诺平均(RANS)和大涡模拟(LES)方法的混合数值模拟。由于该方法提出的时间比较短,目前应用比较少。还有个别研究人员采用涡量输运的方法对机身模型的气动特性进行数值模拟[23]。
目前检验气动特性数值模拟结果准确性的主要方法是将其与风洞试验结果、理论解析值或其他公开发表的数值模拟结果等进行对比分析。ROBIN机身模型风洞试验中给出了机身表面的压力系数,因此研究人员将自己计算得到的机身表面压力系数结果与风洞试验结果进行了对比分析。大多数数值模拟过程中的机身模型长度为2m;由于直升机在大速度前飞为低头姿态,因而计算过程中选取的攻角多为0°或-5°,侧滑角为0°;选取的截面位置通常为x/L=0.0517、0.3074、0.4669、0.6003、1.0008,1.1620等。从对比分析的情况来看,大多数数值模拟结果与风洞试验值比较接近,在靠近凸台和腹部支撑附近的压力系数计算值与风洞试验结果有一定的差别。这是由于在数值模拟过程中没有考虑腹部支撑机构,同时这些区域存在一定的流动分离,而目前数值模拟方法在模拟流动分离时存在一定的困难。
对机身模型进行数值模拟,除了得到机身表面压力系数,还能得到机身附近的速度、流线和涡量变化等流场细节。这些数据对分析机身周围流动的机理,直升机机身减阻增升设计等都非常有帮助。采用面元法和求解N-S方法得到机身表面流线图。从图8中可以看出两种方法得到的前机身和中机身流线图基本上一致,但是在尾梁后下端位置的流线有一定的差别。不同计算方法得到的涡量如图9所示,从图中可以看出两者的主要差别在凸台尾流区的涡量,说明目前数值方法对尾流区的复杂流动模拟还存在一定的误差。
图8 -10°攻角时不同方法得到的机身表面流线图[24]
图9 Q准则得到的瞬时涡结构等值图[22]
早期ROBIN机身模型风洞试验给出的是机身表面压力系数,而在实际工程应用中更关注的是数值模拟方法计算机身气动力和力矩结果的准确度和精度。南京航空航天大学唐韬等[33]在某开口回流式风洞中对ROBIN机身模型的力和力矩进行了测量。试验过程中的来流速度分别为9m/s、18m/s和27m/s。后续在ROBIN机身数值模拟技术的相关研究中可参考此次风洞试验的结果。
通过对ROBIN机身模型数值模拟过程中的网格类型选取、网格数量控制、求解方法选择和与试验对比分析等进行总结与分析,可得出如下结论:
1)机身模型数值模拟过程中有比较多的参数和变量,如网格类型、网格数量、湍流模型等。改变其中的一个参数或变量或许对提高数值模拟结果的准确度有一定的效果。但是如果需要大幅度提高数值模拟结果的准确度和精度,则需要综合考虑多方面因素的影响。
2)在机身模型数值模拟过程中,网格类型、网格数量和湍流模型等的选取不仅要考虑数值模拟结果的准确度等,同时还需要综合考虑可用的计算资源,计算时间的限制,计算方法的连续性等。在工程应用领域,对ROBIN机身模型划分100万左右的非结构网格可满足相应的需求。在数值模拟过程中网格类型、网格数量和湍流模型等可考虑配套使用。
3)未来需要开发多种类型的直升机标准机身模型,以支持不同类型直升机的发展。在机身模型风洞试验过程中,不仅测量机身表面压力系数,还需要测量机身的阻力、升力、俯仰力矩和偏航力矩等力和力矩系数,同时还可以考虑测量机身表面附近的速度分布、尾流区域的涡系等,以方便后续的数值模拟研究工作的开展。
4)在数值模拟研究中考虑制定相应的标准数值模拟方法。在直升机型号研制过程中,由于项目周期和可用计算资源的限制,在一定程度上可以说机身气动特性数值模拟的速度越快越好。这就需要在标准机身模型数值模拟研究中找出兼顾效率和准确度的标准方法。