直升机机身气动特性CFD计算研究进展

2021-09-27龙海斌刘正胜吴裕平

直升机技术 2021年3期

龙海斌，刘正胜，吴裕平

(1.中国直升机设计研究所总体气动研究室，江西景德镇 333001;2.32382部队，北京 100072))

0 引言

在直升机研制过程中，机身气动特性数据是飞行性能、品质和载荷等的设计输入，因此准确地获得机身气动特性数据非常重要。目前主要通过风洞试验和CFD计算两种方法来获得机身气动特性数据。其中风洞试验方法在航空航天等领域应用广泛，具有很成熟的数据修正和准确性验证方法，因此在工程应用领域有很高的可信度。但是风洞试验之前需要设计和制造风洞试验模型，风洞试验过程中又受到风洞的档期安排和测量系统的稳定性等的影响，因此成本比较高，获得机身气动特性数据的周期相对比较长。CFD计算是近几十年来发展起来的数值模拟方法。得益于计算机技术和数值模拟方法的进步，机身气动特性CFD计算的速度越来越快，可有效地加快直升机型号的研制进度。但是，目前尚未形成统一的机身气动特性标准CFD计算方法，CFD计算结果的验证与确认是公认的难题。针对近年来直升机机身气动特性CFD计算方面的进展情况，本文首先介绍直升机机身计算模型的常见构型和气动外形，之后分别从网格类型选择、网格数量控制、求解方法选取和计算结果分析等方面进行整理与分析，最后对机身气动特性CFD计算未来的发展进行展望。

1 机身模型简介

在直升机机身气动特性CFD计算中，计算模型的气动外形通常与风洞试验模型保持一致。用于理论研究和计算方法验证的标准机身模型通常只包含机身和尾梁部分，如ROBIN机身模型、NUAA机身模型等。型号研制过程中的机身模型包含机身、主桨毂、起落架、尾梁、平尾、垂尾和尾桨毂等。其中主桨毂通常为静止状态，部分机身模型的主桨毂为旋转状态。部分机身模型的尾桨毂上带有尾桨叶。武装型直升机机身模型还包含短翼或挂梁、外挂航炮、导弹和火箭弹等。部分直升机机身模型还带有天线、光电吊舱等外挂物。在部分机身气动特性CFD计算模型中还包含风洞支撑机构，计算域的边界为风洞试验段的内表面。

目前已经在飞行和试飞过程中的直升机中，机身与尾梁的过渡情形大致可分为两类：一是机身与尾梁平缓过渡，类似于ROBIN机身模型,如美国的UH-60“黑鹰”直升机、S-97共轴高速直升机、V-22倾转旋翼机以及武装直升机等；二是尾梁比较细长，或者机身尾部有舱门，如ROBIN mod7机身模型、小型无人机和运输直升机等。这些机身的腹部到尾梁的过渡段有比较大的流动分离，因此压差阻力比较大。机身模型如图1-图3。

图1 ROBIN标准机身模型

图2 ROBIN mod7标准机身模型(反装)[4]

图3 机身风洞试验模型[5]

本文中统计的机身类型包括无人直升机、单旋翼运输直升机、单旋翼武装直升机、共轴式直升机、倾转旋翼机等。按机身模型包含的部件来分，机身模型可分为简单的标准机身模型和型号光机身模型，带桨毂、起落架和平尾等部件的机身模型，部分机身模型的桨毂为旋转状态。

2 网格类型

对流体计算域进行网格划分是机身气动特性CFD计算的第一步，而在网格划分之前首先要选取网格类型。目前常用的网格类型有结构网格、非结构网格和上述两者组合的混合网格。结构网格在拓扑结构上相当于矩形区域内的均匀网格，每一层网格上的节点数都相等。针对外形简单的流体域，采用结构网格可以生成比较好的贴体网格，同时求解过程中的数值耗散比较小，因此外形相对比较简单的标准机身模型流场计算多用结构网格，或者进行理论研究需要对流场进行精细化计算时，也采用结构网格。工程应用中的直升机机身包含主桨毂、起落架、平尾和垂尾等部件，部分机型还包含短翼和外挂武器等，难以进行结构网格划分，因此基本上都采用非结构网格。非结构网格包含四面体非结构网格、笛卡尔网格等。部分情况下在机身表面附近划分边界层网格，以进一步提高对机身表面附近复杂流动的模拟能力。但是划分边界层网格容易导致网格质量降低，同时网格数量增长很多。多块网格综合了结构网格和非结构网格的优点，但是各块之间的数据传递存在一定的难度。在部分情况下，可以将非结构网格转化为多面体网格，或者直接划分得到多面体网格。对于相同的机身外形和计算域，多面体网格的数量比较少，可加快CFD计算速度。部分计算过程中需要考虑桨毂旋转，采用嵌套网格方法来计算旋转桨毂和机身周围的流场。在统计的52个划分网格的算例中，非结构网格、结构网格和混合网格的算例数量和所占百分比如表1所示。从表中可以看出，非结构网格占61.54%，说明非结构网格在机身气动特性CFD计算中应用比较广泛，其中采用非结构四面体的算例为22个。在网格划分过程中ICEM软件应用比较多。

表1 各类型网格数和所占比例

针对边界层网格划分，有13个算例划分了边界层网格，占总的算例数量的1/4，说明在直升机机身气动特性CFD计算过程中边界层网格划分相对比较少。这是由于在直升机机身气动特性计算过程中通常带有主桨毂、尾桨毂和起落架等部件，有比较多的拉杆等细小零件，表面形状比较复杂，生成边界层网格的难度比较大；同时，在生成边界层网格时，总的网格数量增长比较多，网格整体质量变差，后续求解计算时可能难以收敛。在划分边界层网格的算例中，第一层网格厚度多为0.05mm或0.005mm。结构网格、笛卡尔网格、非结构网格及其边界层网格分别如图4、图5和图6所示。

图4 结构网格示意图[6]

图5 笛卡尔体网格示意图[8]

图6 非结构及边界层网格示意图[15]

3 网格数量

根据CFD计算方法的基本原理，网格数量越多，则对流体域的划分越细，CFD计算模拟的流动与实际流动的情况越接近。因此在进行理论研究时，为了得到相关区域的流场细节、涡的变化情况等，通常将网格划分得比较细密，网格数量可能达到千万量级。在部分精细化计算时也不考虑网格数量的限制。但是在实际工程应用中，部分情况下机身气动特性计算完成之后，还需要进行旋翼/机身干扰等更复杂的流场计算，加上项目研制周期、可用计算资源等的限制，对机身计算域划分的网格相对比较粗糙。在部分情况下，划分的网格数量多，由于相邻区域的疏密控制等问题，在部分区域容易生成质量很差的网格，造成计算过程中收敛速度变慢甚至难以收敛。

在各类直升机机身气动特性计算过程中，划分的网格数量在几十万到几千万之间。在统计的38个公布的网格划分算例中，超过1000万的算例有12个，约占总数的31.58%，网格数量在100万量级的占50%。各量级网格数和所占比例如表2所示。通常来讲，计算模型的部件数量越多，则划分的网格数量就相对比较多。

表2 各量级网格数和所占比例

部分研究人员针对同一机身模型，划分不同的网格数量进行气动特性CFD计算，以研究网格数量变化对CFD计算结果的影响。文献[21]针对同一个计算模型，划分了4套不同数量的网格。根据CFD计算结果，在划分的网格数量范围内，机身的升力系数和阻力系数均随着网格数量的增加而不断减小，见图7。文献[4]针对ROBIN mod7标准机身模型，划分三套不同数量的结构网格，之后采用OVERFLOW求解器进行气动特性计算，计算过程中选择S-A湍流模型，得到的0°和-5°攻角时的阻力系数如表3所示。从表中的结果可以看出，随着网格数量增加，ROBIN mod7标准机身模型的阻力系数略有减小。

图7 网格数量对气动特性影响示意图[21]

表3 阻力网格数量相关性分析[4]

4 计算方法与湍流模型

在各类直升机机身气动特性CFD计算过程中，常用的计算方法是求解N-S方程方法。早期部分研究人员常用面元法和求解Euler方程的方法。其中面元法是通过在机身表面网格上布置流动的奇点来求解气动问题，计算结果的精度比较粗糙，目前应用比较少。求解Euler方程的方法没有考虑空气的粘性，因此求解速度相对比较快，对计算资源的需求也比较少。在工程应用中方便后续采用相同的方法对旋翼/机身干扰、旋翼/地面等复杂流场进行计算。但是由于没有考虑空气的粘性，因此不能计算空气的摩擦阻力，对比较复杂的流动区域也难以模拟。求解N-S方程的方法主要有雷诺平均(RANS)方法、大涡模拟(LES)方法、直接数值模拟(DNS)方法等。近年来还发展了部分混合方法，比如RANS和LES混合方法等。此外还有格子玻尔兹曼方法(LBM)等可用于机身气动特性CFD计算之中。受限于可利用的计算资源和求解周期限制等因素，目前雷诺平均(RANS)方法应用比较广泛，而大涡模拟(LES)、直接数值模拟(DNS)和格子玻尔兹曼方法(LBM)等在工程领域很少使用。在采用雷诺平均(RANS)方法求解N-S方程的过程中，需要增加方程来使得方程组封闭，即引入湍流模型。目前常用的湍流模型有零方程B-L模型、一方程S-A模型、两方程

和

模型等。其中一方程S-A和两方程

湍流模型应用比较多。在统计的50个采用雷诺平均(RANS)方法的算例中：采用S-A湍流模型的有22个，占44%；采用

湍流模型的有20个，其中采用SST

湍流模型的算例数量为13个，占26%。这是由于S-A湍流模型能给出比较准确的气动特性结果，同时在求解过程中只增加了一个方程，求解的速度相对比较快。

表4 主要的湍流模型应用对比

文献[22]分别采用面元法和CFD计算方法计算了ROBIN机身的绕流场，得到了机身表面压力系数，并与风洞试验结果进行了对比分析，部分结果如图8所示。其中CFD计算方法采用结构网格，求解器为CFL3D，求解过程中采用B-L湍流模型。结果表明CFD计算方法在分离流动和粘性流动中的模拟能力更强。文献[23]采用ICEM软件对ANSAT-M光机身模型进行了以六面体为核心的四面体网格划分，之中在求解N-S方程中分别采用Spalart-Allmaras(S-A), SST

(SST)和改型

(TR)湍流模型，计算结果如表5所示。从表中与风洞试验结果的对比情况来看，SST

湍流模型得到的阻力系数和升力系数与风洞试验结果比较接近。

图8 0°攻角时结果对比图[22]

表5 不同湍流模型计算结果与风洞试验结果对比

5 与试验结果对比

在直升机机身气动特性CFD计算完成之后，需要对计算结果的准确性和有效性等进行验证和确认。目前主要采用与风洞试验结果进行对比分析的方法进行验证与确认。在大部分情况下对比的是整个机身模型的力和力矩系数，部分重点关注位置的表面压力系数等。在风洞试验过程中可以采用PIV方法测量流场细节，同时在CFD计算结果也能比较方便地提取出流线和涡结构等流场细节信息。上述两者可进行比较与分析。文献[24]对某常规单旋翼直升机、某无人直升机和某共轴式直升机(图9中分别用C、W和G表示)光机身气动特性进行了CFD计算，并与风洞试验结果进行了对比分析。其中攻角和侧滑角的范围均为-16°至16°。文献[25]对某常规单旋翼直升机、某无人直升机和某共轴式直升机(图10中分别用C、W和G表示)的机身大攻角和大侧滑角气动特性进行了CFD计算，并与风洞试验结果进行了对比分析。文献[26]对某型武装直升机大侧滑角机身气动特性进行了CFD计算，并与风洞试验结果进行了对比分析(图11)。文献[35]采用Helios求解器对ROBIN mod7机身表面压力系数进行了CFD计算。计算过程中划分了不同的网格，其中：网格1为棱柱网格，数量为870万；网格2为四面体网格，网格数量为2400万；网格5为棱柱网格，网格数量为1550万。从目前CFD计算结果与风洞试验值的对比情况来看(图12)，当攻角或侧滑角变化时，机身气动特性CFD计算结果的变化趋势与风洞试验值基本一致，但在具体数值上有一定的误差，而且部分情况下的误差比较大。目前CFD计算方法可实现0°～360°范围内攻角或侧滑角变化时的机身气动特性计算。机身表面外形变化比较平缓的区域，CFD计算得到的机身表面压力系数与风洞试验值比较接近；而在主减速器整流罩尾部、机身与尾梁的过渡段等外形变化比较剧烈的区域，机身表面压力系数的CFD计算结果与风洞试验值相差比较大。