第二代低阶面元法耦合自由尾迹计算旋翼桨叶气动载荷

王海鹏，王华明

（南京航空航天大学 直升机旋翼动力学重点实验室，江苏 南京 210016）

0 引 言

拉普拉斯方程是分析不可压、无粘、无旋流场的基本方程，用拉普拉斯方程计算飞行器的气动载荷，最常用的数值方法是面元法。面元法经过半个世纪的发展，由最初的求解流场中无升力体表面压力分布到现在的全机气动特性计算，其精度和稳定性已得到大幅度的提升［1］。20世纪80年代末，美国NASA推出了根据低阶面元法开发的计算程序“VSAERO”［2］，该程序通过引入内部速度势边界条件，在减少计算量的同时较大地提高了计算精度，受到世界航空界的广泛关注，被认为是第二代低阶面元法发展的起始。

在面元法发展的同时，旋翼自由尾迹计算方法也在不断进步，由最初的影响系数法［3］、时间步进松弛迭代法［4］、到现在的半隐式预估修松弛迭代法［5］，其计算稳定性与精度也在逐步提高。由于研究侧重点不同，国内外旋翼自由尾迹计算往往集中研究尾迹形状的求解，对于旋翼桨叶多采用简单的二阶升力线［5－6］或不考虑桨叶厚度的升力面模型［4］，这种处理办法不能完全模拟桨叶的三维几何特性，也不能计算桨叶上下表面的压力分布与速度分布。文献［7－9］采用低阶面元法耦合自由尾迹计算旋翼气动特性，但其尾迹计算均采用启动过程的时间步进法，稳定性较差，收敛速度慢。随着计算机技术的发展，雷诺时均N－S方程（RANS）也被运用于求解旋翼的气动特性，但由于旋翼气动环境的特殊性与计算方法本身的特点，计算仍需耗费大量的时间与资源。

本文采用第二代低阶面元法耦合自由尾迹预估修正法计算旋翼在轴流状态下的气动载荷。在计算中考虑桨叶的三维几何特性，采用矩形、三角形单位平面重构桨叶几何外形；桨叶尾迹采用全展尾迹，其中包括三圈自由尾迹（近尾迹）和四圈远尾迹；并利用桨叶翼型风洞试验数据对气动载荷进行粘性修正。通过和CFD方法进行对比分析，证明了本文所采用的计算方法的优点。

1 计算方法

1.1 第二代低阶面元法

面元法总控方程为拉普拉斯方程，其为椭圆型偏微分方程。对于三维拉普拉斯方程常利用格林定理，将其转化为升力体表面的二维问题［10］。

第一代低阶面元法通常只在升力体表面布置一种奇点，或在给定某种奇点分布下布置另一种未知强度奇点，通过物面不穿透边界条件求解未知奇点强度分布。上述计算方法往往在非控制点处存在法向速度“穿透”，对计算结果有较大影响。只有通过加密升力体便面的网格数目，才能获得较好的计算结果。

第二代低阶面元法采用矩形、三角形单位平面重构桨叶几何外形，桨叶面元上布置未知强度的面源与偶极子，桨叶后缘脱出偶极子尾涡（全展），桨叶和尾迹模型如图1所示。通过上述源项组合，空间中任意点P处速度势可表示为以下形式［10］。

其中s为物面坐标点矢量；R（s，p）为物面或尾迹点至计算点的距离；n为物面点处垂直于物面指向流域的法向量；Φ∞p为远方自由流在计算点P点处的速度势；S为桨叶表面；Sw为桨叶尾迹表面，δ（s）与μ（s）分别为布置于桨叶表面的源项与偶极子项。

图1 面元法的桨叶模型与尾迹模型Fig.1 Blade panels and its wake model

上述方程采用混合边界条件求解，具体即为桨叶外表面速度矢量不穿透（纽曼边界条件）和桨叶内表面速度势为常值（狄利克雷边界条件）。如下式所示：

新边界条件的引入使得桨叶表面的源强可直接显式求出而无需方程求解，极大地减少了计算量；内部速度势的引入减少面元间法向速度的“泄露”，降低了面网格分布对计算结果的敏感性，从而提高了计算精度。

尾迹强度使用后缘库塔条件确定，其等于后缘上下面元偶极子强度之差。在已知尾迹形状后，利用边界条件并在桨叶表面、尾迹处划分单位面元可将2式离散为线性方程组，求解此线性方程组可得未知源强与偶极子强度，进而反推出桨叶表面速度与桨叶气动载荷。

1.2 自由尾迹计算

利用尾迹不受力条件，得到尾迹总控方程［6］：

上式为一常微分方程，其中p为尾迹点坐标矢量，r为尾迹涡线，V（p）为在空间P点处的速度矢量合，其中包括桨叶面元，尾迹偶极子面（尾迹涡线）和远方来流的共同作用。本文采用半隐式预估修正法［6］求解尾迹节点位置。

在自由尾迹计算中，涡核模型的选取尤为重要，其直接决定尾迹计算的收敛。本文使用文献［6，11］中推荐使用的涡核模型，涡核半径表达式如下式所示。

式中rc为涡核半径；rc0为初始涡核半径，本文中选取0.1倍桨叶剖面弦长；a为 “Oseen”系数取1.225643；Rev为涡雷诺数；v为动粘性系数；t为尾迹历程时间；a1值在0.01至0.1之间选取［5］。

在自由尾迹计算中，评判两次迭代的尾迹均方根差值来判断尾迹形状是否收敛［5］。尾迹均方根差值计算如下式所示：

1.3 计算步骤与迭代方法

计算时首先划分桨叶面元，根据动量叶素理论估算平均诱导速度，生成初始尾迹，并划分尾迹面元；在此基础上根据第二代低阶面元法求解桨叶面元奇点强度，并求出该尾迹下桨叶气动载荷；接下来根据桨叶与尾迹奇点强度，运用半隐式预估修正法计算旋翼新尾迹；继而判断尾迹形状是否收敛，若满足收敛标准便结束计算，未满足则返回第二步继续计算。详细步骤如图2所示。

图2 计算流程框图Fig.2 Calculation flow chart

1.4 粘性修正

由于基于拉普拉斯方程的面元法无法考虑粘性力，在计算旋翼旋转力矩时，会造成较大误差。针对上述情况，通常采用“附面层修正”方法来考虑粘性作用影响，该方法引入有粘／无粘耦合迭代，用表面源模型模拟粘性带来的附面层位移厚度效应。该方法可较为准确地计算粘性作用，但增加了面元法的计算量与计算耗时。考虑到旋翼工作时，桨叶翼型剖面攻角较小，可忽略粘性力对翼型法向载荷的影响，本文即利用面元法所计算的法向气动载荷，根据相同雷诺数下桨叶剖面的二维翼型CFD计算数据或风洞试验数据，查表反推出对应法向力系数下的轴向力系数，用于修正粘性影响。此法也可快速估算桨叶不同半径处攻角。

2 计算结果对比

为验证上述计算方法的正确性，本文选用广泛引用的“Caradonna”旋翼作为计算算例，其翼型为NACA0012，展弦比为6，无扭转，无尖削，其他参数参见文献［12］。

在本算例计算中，对单片桨叶划分672个面元：其中桨叶沿展向划分10段；沿弦向上下表面划分48段；沿桨叶厚度方向划分4段。每片桨叶拖出3圈自由尾迹，4圈远尾迹，单周尾迹中等角度划分24个尾迹节点，共计1848个尾迹节点。

为便于对比，本文还采用CFD方法求解雷诺时均N－S方程计算该旋翼气动特性。在CFD计算中采用全结构化网格，网格数量约200万，采用剪切应力输运（SST）湍流模型。

图3为以上两种方法计算桨叶压力分布与试验值的对比。从对比结果可以看出，以上两种方法与试验结果有很好的一致性。

图4为以上两种方法计算旋翼在1250r／min时，对应5°、8°和12°总距角时旋翼拉力系数。通过对比可以发现，采用面元法耦合自由尾迹和CFD的计算结果与试验结果有较好的一致性。

图3 桨叶剖面压力分布对比Fig.3 The comparison chart of pressure distribution on blade sections

图4 旋翼拉力系数对比Fig.4 The comparison chart of rotor lift coefficient

图5为以上两种方法计算旋翼在1250r／min时，对应5°、8°和12°总距角时旋翼力矩系数。其中包括面元法修正前与修正后的计算结果。由于面元法只能考虑旋翼力矩中由桨叶诱导阻力造成的旋转力矩，无法考虑由桨叶摩擦造所成的旋转力矩，无摩擦修正的计算结果明显小于CFD计算结果。利用桨叶翼型风洞试验数据修正后，桨叶力矩计算结果与CFD计算结果有较好的一致性。

图5 旋翼力矩系数对比Fig.5 The comparison chart of rotor moment coefficient

图6为上述两种方法计算拉力系数与试验值的相对误差对比图。从图中可以看出，两种方法计算精度相当。

图7为上述两种计算方法耗时对比图，计算平台处理器为“Xeon QC E7330”。通过对比，可见本文方所述面元法耗时仅为CFD计算耗时的1／10左右，大大节省了计算时间。

图6 拉力系数相对误差对比图Fig.6 The comparison chart of lift coefficient relative error

图7 计算耗时对比图Fig.7 The comparison chart of calculation time cost

由于本文中尾迹模型为全展自由尾迹模型，其为从桨叶后缘拖出的涡面，而非一条集中涡线，直接比较尾迹形状有一定困难。图8给出了1250r／min时不同总距角的尾迹计算结果，通过对比可以发现随着总距角的增加，尾迹径向收缩更为剧烈，轴向移动速度也有所增加，这符合一般计算规律。

图8 旋翼尾迹形状（只显示3圈自由尾迹）Fig.8 The figure of rotor wake geometry

3 结论与展望

第二代低阶面元法耦合自由尾迹的计算方法能够模拟旋翼桨叶的三维几何特性，快速计算桨叶表面的压力分布和旋翼的气动载荷。

本文提出的的方法计算旋翼在轴流状态下的气动特性与CFD方法的计算精度相仿，但计算耗费的时间仅为CFD方法的十分之一。

通过引入桨叶的挥舞运动，本文提出的计算方法可拓展至直升机前飞时旋翼的气动特性计算。本文提出的计算方法和旋翼气动噪声分析方法相结合，可以预估旋翼的气动噪声。

［1］ERICKSON L L.Panel methods an introduction［R］.NASA TP 2995，1990.

［2］MASKEW B.Program VSAERO theory document［R］.NASA CR 4023，1987.

［3］BLISS D B.A new approach to the free wake problem for hovering rotors［A］.Annual Forum Proceedings of the American Helicopter Society［C］.1985，463－477.

［4］徐国华.应用自由尾迹分析的新型桨尖旋翼气动特性研究［D］.南京：南京航空航天大学，1990.

［5］BAGAI A.Contributions to the mathematical modeling of rotor flow－fields using a Pseudo－Implicit Free－Wake A－nalysis［D］.College Park：University of Maryland，1995.

［6］曾洪江，胡继忠.一种新的自由涡尾迹计算方法［J］.航空学报，2004，25（6）：546－550.

［7］尹坚平，胡章伟.旋翼表面非定常压力脉动计算的三维自由尾迹非定常面元法［J］.空气动力学学报，1997，15（2）：185－191.

［8］AHMED S R，VIDJAJA V T.Unsteady panel method calculation of pressure distribution on BO105model rotor blades［J］.Journal of the American Helicopter Society，1998，43（1）：47－56.

［9］仲唯贵.直升机旋翼尾涡与尾桨干扰噪声研究［D］.南京：南京航空航天大学，2006.

［10］KATZ J，PLOTKIN A.Low－speed aerodynamics，from wing theory to panel methods［M］. McGraw－Hill，1991.

［11］BHAGWAT M J，LEISHMAN J G.Generalized viscous vortex core models for application to free－vortex wake and aeroacoustic calculations［A］.Proceedings of the 58th Annual Forum of AHS International［C］.2002.

［12］CARADONNA F X，TUNG C.Experimental and analytical studies of a model helicopter rotor in hover［R］.NASA TM 81232，1981.