与涡轮喷气发动机相比，螺旋桨推进系统具有成本、能耗较低，以及工作效率较高的优势，尤其适用于低速轻型飞机及无人机。在概念设计阶段，螺旋桨的选型与设计至关重要，而基于雷诺平均纳维尔-斯托克斯（Reynolds averaged Navier-Stocks，RANS）方程的高可行度计算流体力学（computational fluid dyna-mics，CFD）模拟是目前常用的气动特性评估方法（如拼接网格、嵌套网格等求解技术）。Bryant等［1］采用Ansys Fluent中基于压力的RANS求解器进行了小型螺旋桨的推力性能评估，每个前进比计算都需要14~15 h（160线程，主频2.60 GHz）。Liu等［2］利用基于多重参考系（multiple reference frame，MRF）的RANS方法和低可信度的叶素动量理论（blade-element momentum theory，BEMT）对多旋翼小型无人机进行了气动性能预测，虽然RANS方法精度高于BEMT方法，但每个前进比计算均需消耗900个CPU小时（CPU hour，即3600秒内1GFLOP的机器所完成的浮点运算）左右的计算资源。可见，可信度（精度）和效率是相互矛盾的指标，基于RANS的高可信度螺旋桨气动特性计算方法效率普遍较低。
针对飞行器概念设计阶段的快速迭代需求，低可信度方法是较为合适的选择，但需要对误差特性给出定量评估。开源飞行器建模分析工具NASA OpenVSP包括涡格法、面元法等多种低可信度气动特性评估方法，在国内外研究机构和飞机设计部门均有应用。Zhao等［3］使用OpenVSP基于激励盘理论（actuator disk theory，ADT）的定常涡格法与Ansys Fluent软件基于k-ω SST湍流模型的定常RANS方法，分析了分布式螺旋桨滑流对机翼气动特性的影响，这两种方法得出的升力系数之间存在20%的差异，该值在可接受范围。Wu等［4］同样使用OpenVSP基于ADT的定常涡格法及Simcenter STAR-CCM+软件基于k-ω SST湍流模型的定常RANS方法，计算了机翼螺旋桨滑流，对比试验数据，涡格法升、阻力计算误差分别为3.6%和4.7%，而RANS计算误差分别为1.2%和5.9%。孙宗燕等［5］采用OpenVSP非定常涡格法和面元法对比分析了分布式螺旋桨-机翼模型与原型机翼的升阻特性，计算误差范围为5%~13%。尽管以上研究基本验证了OpenVSP低可信度计算方法的合理性，但是缺乏标模验证及对网格规模、计算步长、方法参数等影响的精细化分析。其次，对多重参考系准定常涡格法的计算特性研究尚属空白。
本文面向飞行器概念设计需求，研究了基于OpenVSP涡格法的螺旋桨气动分析方法。从兼顾计算效率和精度的角度出发，对网格收敛性和计算参数影响进行了基于标模试验数据的对比研究，并在公开文献中首次对基于多重参考系的准定常涡格法进行了验证，形成了若干螺旋桨快速气动分析应用指南，对低可信度气动评估方法工程应用具有一定参考价值。
1 NASA OpenVSP简介
2 OpenVSP非定常/准定常涡格法
2.1 涡格法基本原理
一旦定义了涡环和涡线，就可采用广义Biot-Savart定律［8］求得诱导速度。
经典势流理论未考虑可压缩性效应，为此VSPAERO进行了修正［9］， 长度rp 的平方为
式中：β=（1-M ∞ 2）1/2（其中M ∞<1）为Prandtl–Glauert修正因子。如果自由流平行于x轴，则x坐标随着马赫数的增加而延长，而y和z坐标保持不变，因此等效不可压缩速度势 可以定义为
则涡环对奇点P的诱导速度可表示为
式中：Г为涡环的强度；K为考虑可压缩性对计算影响的因子，如果M ∞<1，则K=2，否则K=1。rβ 由式（5） 定义。
将诱导速度 代入Neumann边界条件方程式［10］，此方程要求奇点处总速度在法向上的分量为零，从而确保流体不穿过升力面。
▽（Φ ∞+Φ'）· n =0
式中：▽为哈密尔顿算子；Φ ∞为自由流速度势；Φ'为扰动速度势； n 为升力面法向量。式（5） 变形后可得到线性方程组
A·Γ = b
式中： A 为系数矩阵，其元素aij= v ij ·n i 为涡环j在奇点i处诱导速度的法向分量（ v ij 为诱导速度， n i 为垂直于奇点i的法向量）； Γ 是由每个涡环强度组成的列向量。
b i =- V ∞ ·n i
2.2 非定常涡格法
VSPAERO“旋转叶片”模式（Advanced/Rotating Blades参数）采用非定常涡格法（unsteady vortex lattice method，UVLM）求解螺旋桨流动。控制方程为拉普拉斯方程。
▽2 Φ*=0
式中：▽2为拉普拉斯算子；Φ* 为流场中任意位置的速度势函数。
涡环在每个奇点的诱导速度、升力面速度和自由流速度的总和需要满足每个奇点的壁面不可穿透边界条件
Γ j=（ V s，i- V w，i- V ∞ ） ·n i
式中：等号左侧为涡环的诱导速度；nh 为涡环数量； V s，i 为奇点i处的表面涡诱导速度； V w，i 为尾涡对奇点i的诱导速度； V ∞为自由来流速度。在VSPAERO中，默认尾涡是不固定的，尾涡的位置用速度场迭代求解，因此尾涡也将表示流线。为了使每个时间步皆满足库塔条件，后缘的涡环切向脱落到尾涡，但其环量强度沿涡线保持恒定［11］。对所有表面涡和尾涡，尾迹的拉格朗日点 r 随着自由流速度和诱导速度移动，表达式为
V = V ∞ + V Blabe+ V Wake
式中： V 为自由流速度 V ∞ 、桨叶（表面涡）诱导速度 V Blade及尾涡诱导速度 V Wake的矢量和。在VSPAERO中，每个涡点位置和速度更新公式为
式中：tn 为当前时间步；tn +1为下一时间步，以此类推。
2.3 基于MRF的准定常涡格法
本文采用了一种基于MRF的准定常涡格法简称 MRF-VLM，该方法将 MRF 与 VLM 在 OpenVSP 平台中进行耦合实现，在目前已发表的文献或其他公开资料中尚未发现相关研究。通过与NASA OpenVSP研发团队交流，本文对该方法进行了应用与验证。MRF基本原理是将旋翼周围的流场包含在一个圆柱形旋转区域内，从而与飞机所在的流场区分开来，在该旋转流域内建立固连于旋翼的旋转坐标系，因此旋转区域与旋翼保持相对静止，进而将旋转区域内的流场转化为旋翼在某一相位角下的瞬态流场进行定常求解［12］。
与非定常方法不同的是，MRF方法只可采用命令行形式执行，计算工况等参数则通过图形用户界面进行设置［13］。
3 网格收敛性与参数影响分析
3.1 APC Electric螺旋桨标模
OpenVSP采用以典型部件为基本单位的参数化建模思想提高外形生成效率，各部件的参数定义方式符合飞机设计的专业习惯。本节以APC Electric 11×8螺旋桨标模为例，开展网格收敛性和迭代参数影响研究，其试验模型［14］和OpenVSP建模外形如图3所示。其中，“11”代表桨盘直径为0.279 4 m，“8”代表螺距为0.203 2 m，建模外形忽略了桨毂部分，截取r/R=0.2处作为翼根。设置螺旋桨转速为3 016 r/min、前进速度为9.915 m/s（对应马赫数0.029），雷诺数为6.741×104，螺旋桨表面积为4.63×10-3m²，平均气动弦长为1.657 m，其余状态量均为海平面大气参数。
3.2 网格收敛性
网格收敛性研究是为了找到OpenVSP求解某一类问题适当的网格规模，既能保证由网格引入的不确定性尽可能小，又能控制网格量避免过度的计算开销。OpenVSP涡格法计算时，翼面网格数量可通过周向节点数Num_W（即机翼或螺旋桨上下表面）和展向节点数Num_U进行调整。展向网格量为Num_U-1，对于涡格法，弦向的网格量为（Num_W-1）/2，两者相乘为表面网格量。但对于机翼，由于OpenVSP以X-Z平面进行对称建模，因此总的表面网格量实际为上述单侧网格量的2倍。为避免不整除的情况，OpenVSP规定机翼应满足Num_W=4n+1，螺旋桨应满足Num_W=8n+1（n为正整数）。
采用UVLM方法对各网格功率系数CP 和拉力系数CT 进行计算。
CP =
CT =
式中：P和T分别为计算出的功率和拉力值；ρ为密度；n为螺旋桨转速；D为浆盘直径。网格收敛性计算结果如图5所示。
从图5可以看出，经过前缘和桨尖加密处理的网格表现出更平稳的网格收敛特性，而无区域加密的网格则在中等网格规模时存在较大波动。综合考虑计算效率，选择总网格量为4 752（Num_W=97，Num_U=100）及加密的网格参数用于下文螺旋桨的计算分析。
3.3 迭代参数影响
采用VSPAERO进行UVLM和MRF-VLM计算时需设置迭代参数，其直接影响计算的稳定性、收敛性和效率。其中，UVLM方法需设置“螺旋桨旋转圈数”（Propeller/Num Revs参数），由于MRF-VLM方法是准定常方法，需设置“尾迹迭代步数”（Advanced/Num It参数）。采用不同迭代参数对APC Electric 11×8 螺旋桨不同前进比状态进行计算。图6、7分别为不同前进比下UVLM和MRF-VLM计算的CP 、CT 收敛曲线。由图6、7可知，最低前进比（J=0.171）状态下两种方法计算均收敛较慢，UVLM和MRF-VLM方法分别需要至少约6圈旋转和30步迭代结果才趋于稳定收敛。
UVLM基于势流理论而不考虑黏性效应，流动分离对计算影响较小，在高前进比下表现出了较好的数值稳定性和收敛性。同时，由于非定常性减弱，采用准定常假设的MRF法也能够快速达到收敛。综上，对于UVLM和MRF-VLM方法，分别采用7圈旋转和40步迭代对于各前进比计算状态是较为合适的设置。表1进一步给出了在上述迭代参数设置下不同数值方法的计算时间对比。本文使用Intel Xeon E5-2698v3 处理器，主频为2.30 GHz，16核32线程。可以发现，作为准定常方法的MRF-VLM计算效率更高。
表1 不同数值方法的计算时间对比（并行线程数均为4） |
| 数值方法 | 圈数/迭代步数 | 计算时间/min |
|---|---|---|
| UVLM | 7 | 20 |
| MRF-VLM | 40 | 4 |
3.4 计算精度和误差分析
采用UVLM和MRF-VLM方法及上述研究获得合适的网格参数和迭代参数对不同螺距的APC Electric 11×8和11×10螺旋桨进行计算，并与试验［14］拉力系数CT 、功率系数CP 和推进效率η。
η=J
从图9中还可以发现，在螺距为8、前进比为0.3~0.61的中等前进比范围内，计算结果误差较小，UVLM方法和MRF-VLM方法计算误差分别在6.1%和3.6%以内，可知在合适的网格和迭代参数下，计算精度满足概念设计阶段的要求。在较低和较高前进比下，误差有所增大。然而在螺距为10时，UVLM方法的误差普遍大于MRF-VLM方法，尤其在较低和较高前进比下误差波动较大。相比之下，MRF-VLM方法表现更加稳定，除最低前进比（J=0.175）外，其余误差控制在10.43%以内。结合前一小节对计算时间的分析，也基本证明了MRF-VLM方法具有较高的计算效率和精度。
4 MRF-VLM方法计算验证
4.1 NACA 5868-9螺旋桨标模
可得到螺距H分别为2.56、4.47、6.72 m。采用海平面大气参数，网格和迭代参数与上文APC螺旋桨的设置相同。
4.2 计算精度和误差分析
由于螺旋桨推力主要取决于升力，而功率主要取决于阻力，因此升阻力将影响螺旋桨的效率。当前进比较低，尤其是在高桨距角操作时，螺旋桨的叶片可能会进入失速状态，导致升力突然下降，阻力增加。结合上文对APC螺旋桨的计算可知，涡格法对于非线性失速阶段的预测并不准确。
5 结论
本文研究了OpenVSP非定常和准定常（多重参考系）涡格法两种计算方法在螺旋桨高效气动评估中的计算特性，通过标模算例分析了网格和迭代参数对收敛性的影响，给出了此类低可信度计算方法的计算精度和效率。主要结论如下：
1）对于单个螺旋桨计算，推荐单个叶片周向和展向布置100个左右节点，同时进行前缘和桨尖区域加密，迭代参数为7圈（对于UVLM方法）或40步（对于MRF-VLM方法），能够确保在不同工况下的计算稳定性和收敛性。
2）UVLM和MRF-VLM方法在低桨距角、中等前进比条件下的计算结果与实验数据较为吻合，但在高桨距角和较低、较高前进比下UVLM方法的误差较大（超过10%），而MRF-VLM方法表现更为稳定，且计算效率较高，符合概念设计阶段的精度和效率要求。
3）除计算模型与试验螺旋桨之间存在一定的几何差异外，涡格法基于势流理论的局限性使得UVLM和MRF-VLM方法无法准确模拟较强的流动非线性。针对这一问题，可提出适当的非线性修正方法，以提高涡格法在高桨距角、高低前进比下的预测精度。