基于用户自定义函数的直升机旋翼CFD模拟加速方法的设计和验证

喻延福+杨志刚+贾青+Ivan+Dobrev

摘要：为在直升机旋翼气动性能数值模拟时简化建模过程、缩减计算时间，利用用户自定义函数（User Defined Function，UDF）设计混合模型盘激励模型和线激励模型，并对简单旋翼的悬停工况进行模拟.与风洞试验结果的对比表明：所设计的混合模型在简化旋翼模型的同时，能有效地计算旋翼的气动特性，模拟旋翼悬停时的流场，具有正确性和可行性；盘激励模型作为定常计算模型能够较快地计算得到旋翼的气动性能，缺点是不能体现每个桨叶对流场的单独作用；所设计的线激励模型在计算时由于所用的诱导速度为平均值，所以计算结果中旋翼效率比实际值偏高；通过与粒子图像测速（Particle Image Velocimetry，PIV）测量结果对比发现，线激励模型能较好地模拟出桨尖涡的分布.

关键词：直升机； 旋翼； 气动特性； 风洞试验； 用户自定义函数； 粒子图像测速； 数值模拟

中图分类号： V211.52

文献标志码： B

Abstract：To simplify the modeling process and reduce computation time during the numerical simulation of aerodynamic performance of helicopter rotor， the actuator disk model and the actuator line model of hybrid model are developed using User Defined Function（UDF） and a simple rotor in hovering is simulated. The results are compared with the results of wind tunnel test. It is indicated that， while the rotor model is simplified， these hybrid models are effective and feasible to calculate the aerodynamic characteristic of helicopter rotor in hovering and simulate the flow field； the actuator disk model is a steady solution model which is able to calculate the aerodynamic characteristic of rotor in a short time， but it cannot separately present the effect of every blade on the flow field； as the fact that an average induced velocity is used in the actuator line model， the calculated result of rotor efficiency is a little higher than the actual value； by the comparison between the hybrid model results and the Particle Image Velocimetry（PIV） measurement results， it was found that the actuator line model could well predict the distribution of blade tip vortex.

Key words：helicopter； rotor； aerodynamic characteristic； wind tunnel test； user defined function； particle image velocimetry； numerical simulation

0 引 言

对风力机的风轮或直升机旋翼等旋转流体机械进行数值计算时，如使用传统的CFD方法，需将桨叶转子几何表面的边界条件设为壁面，在壁面生成边界层并在适当区域加密，这对模型的网格质量和数量有较高要求；如要研究桨尖涡和进场尾迹，还需进行非定常计算，不仅工作量大，且对计算资源的要求高.例如，在风力发电农场的气动研究中，上游风力机的桨叶转子产生的尾流会对下游风力机产生影响，为仿真研究这种现象，需要对多个风力机同时进行非定常数值计算，对计算资源有极高的要求，计算时间长.因此，有必要使用一种能加速计算同时可正确模拟桨叶与空气间相互气动作用、仿真尾流分布情况的仿真方法.

混合模型[1]是进行风力电厂CFD计算时应用较多的一种方法.该模型简化真实桨叶转子模型实现加速计算.在该方法中，真实的桨叶转子几何模型被外形较简单的特殊几何模型代替，模型的网格数量较少且无须生成边界层，能模拟真实的桨叶转子与空气的相互作用，对计算资源的要求较低.之所以称之为混合模型，是因为在数值计算时将动量理论和叶素理论与数值模拟时控制方程的迭代求解相联系：在每一步迭代结束时，通过理论方法计算气流流经桨叶转子时所受的力，再将此力以源项的形式代入到控制方程的下一步迭代中，从而模拟桨叶对气流的作用力.根据代替桨叶转子的几何模型不同，针对风力机数值模拟的混合模型主要发展出3种形式：盘激励模型 [2]、线激励模型[3-4]和面激励模型[5-9].在盘激励模型中，用一个类似桨盘的圆盘代替桨叶转子，对流经圆盘的气流施加单位体积力或单位面积力.因此力沿圆盘的方位角方向均匀分布而沿径向有变化，故其缺点是无法体现每个叶片的单独作用.在线激励模型中，每根桨叶被一条沿桨盘半径方向的线段或细杆代替，线段或细杆上所提供的力沿桨叶展向变化.与盘激励模型相比，此模型能体现每个桨叶对气流的单独作用以及桨尾和桨尖涡的作用，改善流场特征的体现.在面激励模型中，每个桨叶被其中心面所代替，在数值计算时通过边界条件在此中心面上提供压力差.与前2个模型相比，面激励模型模拟得到的尾流更接近真实情况，但在流速较大时，由于流场的分离作用，仿真结果不能与试验结果很好吻合.[10]

上述研究主要集中在对风力机的数值仿真方面，对于直升机旋翼这种通过旋翼旋转对气流加速、来流速度未知的桨叶转子情况并未涉及.本文将上述风力机风轮模拟混合模型的方法应用于直升机旋翼的数值模拟中，基于旋翼的气动理论，借助用户自定义函数（User Defined Function， UDF）设计一套针对旋翼气动模拟的盘激励模型和线激励模型，并通过对某小型旋翼进行悬停状态的数值计算和风洞试验，分析计算与试验结果，对所设计的模型予以验证.对面激励模型的流场分离现象和模拟桨叶振动时的局限性，本文不进行讨论.

1 混合模型设计

1.1 理论计算

在混合模型中，代替旋翼的简单几何模型通过对空气施加单位体积力或单位面积力模拟桨叶与空气的相互作用，UDF理论计算部分的主要目的是计算简单几何模型中每个微元上应施加的力.

以盘激励模型为例（见图1），用圆盘代替旋翼后，根据微分原理，桨叶上宽度为dr的叶素被圆盘上径向宽度为dr的同心圆环微元代替.半径为r处的叶素所施加的作用力被均匀分布在半径同为r的圆环微元上.针对图1中半径为r处的叶素，其二维翼型的受力分析见图2.

1.2 UDF编写

利用CFD求解商业软件FLUENT中的UDF功能实现混合模型中计算功能.在UDF中通过特定的函数得到每个网格单元的体积、流体速度矢量、三维坐标等变量的值，进而将其利用在自定义函数的计算中.模型中所用函数有x，y和z等3个方向的源项函数，计算旋翼的简单几何模型内单位密度力，每步迭代后计算旋翼升力、扭矩功率以及旋翼效率的函数.由于线激励模型的非定常计算量较大，故在UDF进行单节点试算成功之后，还需进行并行化改写，以满足多节点并行计算的需要.加入UDF后的控制方程的迭代求解步骤见图3.

利用叶素理论求解叶素的升阻力，需要在给定翼型和迎角时知道CL和CD.在UDF的编写过程中需提前录入所需翼型的升阻特性数据.为与试验保持一致，本算例选用NACA0015翼型，由于雷诺数较低，约为Re=0.1×106，故采用文献[11]中NACA0015在雷诺数较低情况下的试验数据，见图4和5.在升阻特性曲线上截取若干数据点，UDF中通过插值函数进行取值.线激励模型中计算源项时，对于诱导速度的获得，未用UDF获取网格节点处转轴方向的速度，原因是此处网格的局部气流速度不能代表宏观的旋翼上游诱导速度.本文利用动量理论计算旋翼的上游诱导速度的平均值，相当于将诱导流场看作一种流速均匀的流场.这种流场情况只有当桨叶有理想情况的扭转时才能实现，而当桨叶无扭转时，这是一种会使得旋翼的效率比实际值偏高的简化方法.[12]

2 数值计算

对具有2个无扭转桨叶的小型旋翼进行悬停状态的数值计算.旋翼半径R=290 mm，忽略桨毂部分的几何结构，桨叶在桨毂部分的无效区域径向长度为41 mm.针对盘激励模型和线激励模型建立相应的CFD模型，模型中网格区域分为2个部分，即旋翼区域和旋翼周边圆柱形计算域.在盘激励模型中，旋翼被厚度为0.01 mm的圆盘所代替，圆盘半径为R.在线激励模型中，旋翼区域为半径360 mm，厚度40 mm的圆盘形旋转区域，旋翼的2片桨叶被2个与桨叶等长（即L=R=290 mm）、截面半径r=5 mm的圆杆代替，被圆盘包裹，见图6.桨毂无效区半径均设为41 mm；计算域半径取3R=870 mm，上游长度取5R=1 450 mm，下游长度取15R=4 350 mm.在求解计算时借助网格滑移，使线激励模型的旋转区域网格以一定角速度绕z轴旋转，其余计算域网格静止.参考二维翼型压力分布，升力的产生集中在翼型前半部分，故圆杆截面的位置位于桨叶截面翼型的前端，见图7.由于盘激励模型结构简单，所以圆盘内和计算域均采用六面体结构网格，计算域靠近圆盘下游的区域适当加密，总网格数量约为165万个.对于线激励模型，在旋转区域的2个圆杆内部以及圆盘周围计算域采用结构网格，圆盘内其余区域采用四面体非结构网格，计算域下游桨尖涡区域网格适当加密，总网格数量约为442万个.

本算例马赫数小于0.3，求解不可压N-S方程.盘激励模型进行定常求解，湍流模型选择标准k-ε模型.计算域入口边界条件设为压力入口，出口为压力出口，其余壁面为无剪切应力壁面.圆盘内部网格区域设置源项UDF，圆盘表面面网格为内部面网格.线激励模型进行非定常求解，湍流模型选择k-ω SST模型.计算域出入口及壁面边界条件与盘激励模型相同，圆杆内部网格区域设置源项UDF，圆盘内的网格为动网格，以角速度ω绕旋翼转轴运动，周围计算域的网格静止，利用网格滑移，圆盘面网格边界条件设为交界面.动量、湍动能和耗散率的离散采用二阶迎风差分格式，压力-速度耦合使用SIMPLE算法求解.

通过改变桨距角θ，利用2种混合模型分别模拟θ分别为4°，6°，8°，10°和12°时旋翼在转速为匀速2 000 r/min时的情况.

3 风洞试验

3.1 试验台架布置

为验证所设计的混合模型的正确性，利用试验方法验证数值计算结果.试验在法国国立高等工程技术大学风洞实验室中的六分力天平上进行.试验风洞为半开口式回流风洞，试验段长度为2 m，喷口大小为1.35 m×1.8 m.旋翼总成及支架安装在试验段的六分力测力天平上，使旋翼转轴在试验段中间位置，方向水平，离地高度700 mm，见图8.旋翼模型尺寸与前文一致.通过更换桨毂上固定桨叶的铰链改变θ.旋翼安装在水平转轴上，转轴和电机之间通过圆柱轴扭矩传感器T20WN连接，通过该传感器测量扭矩和转速.通过在控制系统上改变电机的输入频率改变旋翼转速.由于旋翼尺寸较小且效率较低，旋翼产生的作用力不会在整个风洞中产生循环气流，故试验时风洞的风机不开启.开启旋翼电机前先将天平y和z方向受力调零.开启旋翼并达到稳定速度后，微调圆形托盘使得y方向受力为0，即保证旋翼转轴方向与z轴方向一致.试验时，在相同θ下测量旋翼转速不同时的升力、阻力矩和电机功率，然后更换铰链改变θ进行下一组测量.试验测量θ分别为2°，4°，6°和10°的数据，但2°时旋翼产生的升力较小，天平测量精度不够，故没有记录到最终结果.

3.2 测量

试验利用粒子图像测速（Particle Image Velocimetry，PIV）技术测量旋翼下游的流场分布.见图9.Litron Nano公司的YAG激光器安装在试验段的上方，可产生200 mJ的激光脉冲.激光透过下方的柱透镜可形成铅垂方向的平板状光路，调整角度使光路穿过旋翼转轴.在试验段的进风口导入由烟雾生成器生成的橄榄油微粒烟雾，微粒直径为1～10 μm.拍摄所用的CCD相机为Dantec FlowSense 4M，分辨率2048 dpi×2048 dpi，镜头为Nikkor AF-S 105 mm f/2.8G ED IF.为获得旋翼旋转的位置，在桨叶上粘贴光学位置传感器，连接到电脑的PIV控制系统，同步激光脉冲和相机快门，设定当旋翼转到某一角度时进行拍摄.根据流场速度，连续拍摄2张照片的时间间隔为20 μs.拍摄时同一工况拍摄200组照片.PIV测量现场图见图10，可以明显看到由于桨尖涡对烟雾微粒的离心作用形成的黑点.

4 结果与分析

4.1 旋翼气动性能

2种模型在不同θ下的气动性能的数值结果和真实旋翼的试验结果见图11和12.由图11a可知：旋翼相同θ相同时线激励模型计算得到的升力比盘激励模型的结果大，且θ越大二者相差越大.由图11b可知：在扭矩功率方面，两者结果基本重合，从数值上看在θ=12°时偏差为4%，θ<12°时偏差小于1%.究其原因，是由于线激励模型中诱导速度的计算采用平均值在整个桨盘平面上均匀分布的假设，从而使得旋翼效率比实际值偏高；而对于阻力矩的值，从图5中可看出迎角为0～13°时，翼型的CD值变化很小，所以由于诱导速度的偏差对基于CD计算的扭矩功率影响很小.与试验结果相比，混合模型得到的升力和功率在量级上和趋势上均与试验保持一致，数值上升力和扭矩功率试验值略微偏大，与数值结果的误差最大约为5%.分析原因，是由于实际试验中旋翼在转速2 000 r/min时总成会产生轻微共振，为避开此共振点，在试验时实际采集的结果是在旋翼转速约为2 200 r/min时的数据.

对比分析旋翼悬停状态3个无量纲因数的计算结果，见图12.盘激励模型与线激励模型相比：对于大小与升力相关的拉力因数CT和悬停气动效率FM，线激励模型的值均偏大，其原因与前面所述升力偏大的原因相同；对于与扭矩功率相关的扭矩因数，二者偏差均很小.与试验结果进行对比时，试验时实际转速偏大带来的影响在对结果的无量纲化时可被消除.因此可看出：当θ不同时，混合模型方法得到的无量纲因数在量级上和趋势上均与试验结果相一致，但数值大小上试验结果中的功率因数CP偏低较多，θ=10°时与数值结果偏差约25%.究其原因，可能是由于基于叶素理论进行计算的混合模型没有考虑实际旋翼旋转时桨叶振动以及桨尖涡导致的功率损失.

4.2 流场分析

θ相同且转速相同时，对2种混合模型仿真结果与风洞试验结果的流场进行对比分析.本文主要展示θ=10°时旋翼匀速转动使得流场稳定后旋翼周围的流场分布情况.

混合模型仿真结果在x=0平面的压力云图见图13，桨盘在z=0位置，转轴在y=0位置.

由图13可以看出：由于旋翼的作用，桨盘两侧形成压力差，从而对旋翼产生升力.盘激励模型由一个圆盘代替旋翼，压力沿方位角方向均匀分布，而线激励模型由圆杆代替桨叶，压力差集中在桨叶周围，故图13a中所示的压力差远小于图13b中圆杆两侧的压力差.桨盘中心由于桨毂无效区的存在，两侧无压力差.

虽然在盘激励模型对空气的作用力计算中可体现切向的分量，对流场沿旋转切向方向也有加速，但由于盘激励模型中旋翼对空气的作用力沿方位角方向上均匀分布，所以不能像线激励模型那样体现每个桨叶单独对空气的作用.

z=0平面的速度云图见图14，可更清晰地看出盘激励模型的这一缺陷.盘激励模型的优点是模型简单，数值计算时进行定常求解，故能在短时间内算出旋翼的气动性能参数.当需要短时间内计算旋翼的气动性能时，盘激励模型有一定优势.

观察桨尖涡，Q=1 200 s-2时线激励模型等值面涡量图见图15.由此可看出：线激励模型可模拟每根桨叶形成的桨尖涡，桨尖涡螺旋向下发展，螺旋半径减小，这与理论相一致.流场在x=0平面涡量不变量Q的计算结果云图见图16.从图16可看到桨尖涡的位置.试验进行PIV测量时，当旋翼转到同一桨叶处于铅垂向上位置时进行拍摄.用Dantec DynamicStudio 2.30对θ=10°，转速分别为1 100和2 200 r/min时的200组拍摄照片进行后处理，得到的平均速度场云图见图17.转轴与x轴重合，桨叶位于x=300 mm的位置.图17显示旋翼下游的平均速度场分布和桨尖涡的位置，图中从左向右数，每2个涡为一对，由旋翼旋转1周时2个桨尖分别产生.其桨尖涡的中心位置与仿真结果一致，但由于数值模拟的数值耗散以及网格滑移交界面的影响，仿真结果中第2对及之后的桨尖涡显示不明显.

5 结 论

借助UDF建立能够简化建模流程、缩减计算时间、针对旋翼CFD仿真的盘激励和线激励混合模型，对简单旋翼的悬停工作状态进行仿真计算，通过与风洞试验结果进行对比，主要结论如下.

1）所设计的混合模型在简化旋翼模型的同时能有效计算旋翼的气动性能并仿真旋翼的悬停状态，具有正确性和可行性，有进一步开发和应用价值.

2）作为定常计算模型，盘激励模型能够较快计算得到旋翼的气动性能，缺点是不能像线激励模型那样体现旋翼每个桨叶对流场的单独作用.

3）目前设计的线激励模型在计算时因所用的诱导速度为平均值，故计算结果中旋翼效率比实际值偏高.这一点需在今后的研究中进行改善.

4）通过与PIV测量结果进行对比表明：线激励模型能较好地模拟桨尖涡的分布情况，在研究桨尖涡对直升机工况的影响时有效.

参考文献：

[1] VERMEER L J， SRENSEN J N， CRESPO A. Wind turbine wake aerodynamics[J]. Progress in Aerospace Sciences， 2003， 39（6/7）： 467-510. DOI： 10.1016/S0376-0421（03）00078-2.

[2] MIKKELSEN R. Actuator disc methods applied to wind turbines[D]. Technical University of Denmark， 2003.

[3] SRENSEN J N， SHEN W Z. Numerical modeling of wind turbine wakes[J]. Journal of Fluids Engineering， 2002， 124（2）： 393-399. DOI： 10.1115/1.1471361.

[4] DOBREV I， MASSOUH F. Etude d'un modèle hybride pour représenter l'écoulement à travers un rotor éolien[C]//Procédure de 17ème Congrès Franais de Mécanique. Troyes， 2005.

[5] MASSOUH F， DOBREV I， RAPIN M. Numerical simulation of wind turbine performance using a hybrid model[C]//Proceeding of the 44th AIAA Aerospace Sciences Meeting and Exhibit. Nevada， 2006： 1-10. DOI： 10.2514/6.2006-782.

[6] DOBREV I， MASSOUH F， RAPIN M. Actuator surface hybrid model[J]. Journal of Physics： Conference Series， 2007（75）： 12019. DOI： 10.1088/1742-6596/75/1/012019.

[7] DOBREV I， MASSOUH F， MAALOUF B. Lifting surface method for prediction of rotor vortex wake[J]. 19ème Congrès Franais de Mécanique， 2009： 24-28.

[8] SHEN W Z， ZHANG J H， SRENSEN J N. The actuator surface model： a new Navier-Stokes based model for rotor computations[J]. Journal of Solar Energy Engineering， 2009， 131（1）：11002. DOI： 10.1115/1.3027502.

[9] WATTERS C S， MASSON C. Recent advances in modeling of wind turbine wake vortical structure using a differential actuator disk theory[J]. Journal of Physics： Conference Series， 2007（75）： 12037. DOI： 10.1088/1742-6596/75/1/012037.

[10] MEMON A A. Développement d'un modèle de surface active pour améliorer la représentation des charges aérodynamiques sur une pale éolienne[D]. Arts et Métiers ParisTech， 2012.

[11] SHELDAHL R E， KLIMAS P C. Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines： SAND-80-2114[R]. Sandia National Labs， Albuquerque， USA， 1981.

[12] LEISHMAN J G. Principles of helicopter aerodynamics with CD extra[M]. Cambridge， Eng： Cambridge University Press， 2006.

（编辑 武晓英）