螺旋桨非定常滑流的高效数值模拟研究

杨小川，王运涛，王光学，张玉伦

（中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000）

0 引 言

对于螺旋桨飞机而言，由于螺旋桨的存在，使得其在气动特性上与常规布局飞机相比，存在更为复杂的流动现象。其中，预测螺旋桨滑流对飞机机翼等部件气动性能的影响，一直是螺旋桨飞机研究的重点［1］。

关于滑流对全机气动特性影响的研究很多。常用的方法大多为简化的作用盘／激励盘模型［2－5］，特别是等效盘模型［6－7］，将螺旋桨桨叶视为无厚度的圆盘，经过盘面的气流被加速加转，能较好地考虑桨叶数目、形状、安装角和转速等多种因素对滑流的影响，得到了较好的计算结果。但在模拟螺旋桨滑流的非定常效应上，基于动量理论的简化模型不能满足计算要求。为了解决这个问题，国外Rooesnboom［8］等采用拼接网格技术模拟了螺旋桨转动，网格采用混合非结构网格生成，得到的计算结果与PIV试验对比较好；国内任晓峰［9］、张小莉［10］、许和勇［11］等均进行了非定常滑流的数值模拟，并进行了细致分析。

本文运用自主开发的大型“亚跨超CFD软件平台”（TRIP3.0）中的非定常计算模块［12］，并综合应用动态拼接网格技术、多重网格技术和大规模并行计算技术，其中方程采用非定常Euler方程，并运用双时间步推进方法求解控制方程，格式采用MUSCL＿Roe格式。

1 计算方法

1.1 控制方程

螺旋桨在旋转过程中，经过螺旋桨的气流是完全非定常的，因此，为了模拟螺旋桨滑流的非定常特性，必须采用非定常方法进行计算。假设oxyz是惯性笛卡儿坐标系，忽略彻体力，则非定常欧拉方程可表达为：

其中ρ、（u，v，w）、p、e 和h 分别代表密度、（x，y，z）方向的绝对速度分量、压力、内能和总焓。

1.2 时间离散

时间方向采用双时间步长推进［13］，非定常欧拉方程可表达为：

采用二阶后向差分，将式（1）在时间方向上展开：

这是一个关于Qn＋1的非线性方程。可通过引入一个伪时间步Δτ来求解Qn＋1。

在m→∞，Qn＋1＝Qm＋1，于是在伪时间上变成一个定常问题，而且在伪时间上没有时间精度要求，完全可以采用定常求解方法求解，各种加速收敛的措施也可以继续使用。只要伪时间上的定常问题收敛到一定程度，真实时间上的精度仍然是二阶的。

1.3 动态拼接网格

拼接网格又被称为面搭接网格，它的特点是在网格块的公共面上，网格线不需要连续对接，因此各块网格之间的拓扑关系在交接面上被切断，各块网格之间是基本独立的。在网格生成时，对利用拼接处理的网格块可以单独生成网格，不需要顾及相邻网格块的拓扑关系，从而使得网格生成的难度大大地降低。另外，由于切断了各块网格之间的网格线，因此对于局部进行加密的网格点就不会传递到不需要加密的空间或其它网格块中，从而节省了大量网格单元［14］。

动态拼接网格采用面积加权的线性插值方式，在计算过程中，螺旋桨在一个时间步长旋转3°，并进行一次网格拼接，每个周期旋转（120／桨叶数）次。该方法对亚、跨、超流场均能适用，同时程序的实现也比较简单。在网格生成中，既能很好的解决“剪刀缝”问题，又能将运动网格（特别是旋转的网格）和静止网格进行拼接处理，计算复杂的非定常流动问题。

1.4 多重网格技术

在NS／Euler方程求解过程中，通过误差分析可知，在迭代求解时，高频误差可以被快速地消除，而低频误差则慢得多，从而大大降低了流场的收敛速度。多重网格的思路实际上就是密网格经过网格粗化后，其低频误差会变成粗网格上的高频误差，利用一系列不断粗化的网格来迅速消除各层网格上的高频误差，即相当于消除了密网格上的低频误差，进而达到加速收敛的目的。

本文采用FAS方法［15－16］的三层网格加速技术，并与拼接网格技术相结合，大大提高了计算收敛速度。由于螺旋桨飞机的网格复杂，网格量较大，运用拼接网格技术在一定程度上减少了网格节点数，但仍不能达到快速收敛的目的。通过多重网格和拼接网格技术的综合运用，在模拟复杂外形上具有较强的实用性和高效性。

1.5 大规模并行计算技术

在计算空气动力学中，并行计算通常采用物理区域分割方法，实现粗粒度高效计算，在编程上采用单控制流多数据流（SPMD）模型，采用 MPI实现消息传递，编写具有多指令流多数据流（MIMD）风格的并行计算程序。

并行原理是：将整个流动区域分割成n个子区域，分配给n个CPU计算，把子区域的初始流场信息、几何信息（网格坐标、标识号）分别装载入各子区域对应的CPU内存中，在每一个CPU中启动计算进程，由主进程调度各CPU的计算。在每一次全场的扫描过程中，由各CPU完成子区域的计算并在边界完成数据交换（各CPU间的通信），由主进程收集全场数据，完成收敛准则判别，并按需要进行写盘等操作。

2 单独螺旋桨验证算例

2.1 计算条件

计算外形：单独螺旋桨，桨叶剖面属于Clark Y翼型，双叶，螺旋桨直径为0.76m，螺旋桨0.75r处安装角为16°或22°。网格采用拼接网格生成（如图1所示），网格量为430万。

计算状态：来流速度30m／s，转速分别为2420r／min和3200r／min，分别对应的安装角为22°和16°。

图1 螺旋桨计算网格Fig.1 Mesh of the propeller

2.2 与试验数据对比

本文采用非定常方法对文献［6］中单独螺旋桨算例进行了模拟，并与等效盘模型［17］计算结果进行对比分析，其中模型试验数据由南航内流研究中心提供。文中重点选取了两个状态进行验证，表1、表2给出了计算得到的结果与试验数据的对比情况。

表1 非定常方法与等效盘方法结果对比（状态一）Table 1 Comparision of unsteady and actuator model methods（case 1）

表2 非定常方法与等效盘方法状态二结果对比（状态二）Table 2 Comparision of unsteady and actuatormodel methods（case 2）

由表1和表2可以看出，采用非定常计算的扭矩值与试验值吻合很好，拉力较试验值有一定差距，等效盘方法得到的拉力值与试验值更接近，但扭矩值较试验偏大。

图2为等效盘模型的原理示意图，通过等效盘模型初步估算桨叶不同站位处的受力情况（如表3所示），可以看出桨叶上阻力对拉力和扭矩的贡献量，相对升力而言是小量，模拟的关键在于对桨叶升力的计算。

图2 等效盘模型原理示意图Fig.2 Schematic diagram of the actuator model method

表3 基于等效盘原理的不同桨叶站位受力分析预测情况（状态一）Table 3 Lift and drag coefficient on different sections based on actuator model methods

3 螺旋桨滑流对全机气动特性影响

3.1 计算模型

计算外形为某上单翼双发螺旋桨飞机，六叶螺旋桨，机翼参考面积约为70m2。其中螺旋桨为旋转部件，定义顺流方向看对螺旋桨逆时针旋转为正。

3.2 计算网格

由于机身和机翼等部件周围网格静止不动，而螺旋桨处网格需相对机身、机翼等高速旋转，故网格由旋转网格和静止网格两部分组成，包含螺旋桨的旋转网格，外形为一个有厚度的圆盘，圆盘的直径选取螺旋桨直径的1.3倍。在计算过程中，通过该区域网格的旋转来模拟螺旋桨的转动。圆盘以外的机翼、机身等区域为静止网格，静止网格与旋转网格采用拼接方式连接。

计算网格为结构化网格，半模网格量941万，其中旋转网格量240万。图3、图4给出了网格在拼接处的切面示意图，其中图3为静止网格切面，可以看到无网格的区域正好将螺旋桨包围，即为旋转网格部分（如图4所示），两个网格在交界面处进行拼接处理，图5给出表面网格示意图。

图3 静止网格切面示意图Fig.3 Mesh slice of the stationary subzone

图4 旋转网格切面示意图Fig.4 Mesh slice of the rotational subzone

图5 表面网格示意图Fig.5 Mesh of the surface

3.3 计算结果与分析

为了研究滑流对全机流场和气动特性影响，分别对不同迎角、不同螺旋桨转速进行计算分析。计算在银河集群中心进行，采用24核进行并行计算，单核网格量约为39万，计算周期约为2天。

图6给出了螺旋桨转速为600r／min，迎角为0°时，拉力系数随时间的变化情况。由图6可知，在整体上，拉力大小出现明显的周期性正弦变化。这主要是受到发动机短舱和机翼等部件的周期性影响，螺旋桨转动产生的拉力出现周期性变化。

3.3.1 滑流对全机流场影响

由于螺旋桨旋转过程中，桨叶是转动的，相位角也是随时变化，这种非定常的滑流对全机流场影响是动态变化的。

图6 拉力系数随时间变化Fig.6 The thrust coefficients change with time

图7、图8分别为上、下表面压力系数云图在有、无滑流情况下的对比图，其中螺旋桨转速为900r／min，迎角为0°。以上表面为例，在有滑流时，内侧（短舱与机身之间）机翼上表面出现明显的低压区，外侧（短舱与翼尖之间）机翼上表面压力分布出现明显变化。这说明整个螺旋桨后方都是明显的滑流区，且滑流对机翼压力分布产生较大影响。由于螺旋桨逆时针旋转，内侧机翼有效迎角增大，外侧机翼有效迎角较小，而离滑流区较远所受的影响则越小。

图9所示为机翼展向站位分布示意图。图10表示螺旋桨转速为600r／min，迎角为0°某瞬时状态下，各站位马赫数云图对比。其中图10（a）～图10（d）可观察到螺旋桨工作时桨叶拖出的涡系。

图7 上表面压力系数云图有、无滑流对比Fig.7 Comparison of the pressure coefficients on up－wing

图8 下表面压力系数云图有、无滑流对比Fig.8 Comparison of the lift coefficients on under－wing

图9 机翼展向站位分布Fig.9 Different sections on the wing

图10 有、无滑流马赫数云图对比Fig.10 Mach number contour of the different sections

图11 滑流对升力和俯仰力矩系数的影响Fig.11 Comparison of lift and pitch moment coefficient

3.3.2 不同迎角对全机气动特性影响

在螺旋桨转速为600r／min状态下，图11给出了不同迎角下滑流对升力和俯仰力矩（不含螺旋桨部件力）系数影响情况。可以看出，滑流使得飞机升力增大。仰角越大时，滑流对升力的影响越明显，在α＝20°时，升力系数增量达到0.36；在大迎角时，滑流对全机俯仰力矩影响最为明显。

图12给出了螺旋桨转速为600r／min，迎角为0°状态下，滑流对机翼展向站位压力系数影响情况。可以看出，滑流对机翼各站位压力影响明显，且越靠近螺旋桨后方，影响越明显。

图12 滑流对机翼展向站位压力系数的影响Fig.12 Comparison of pressure coefficient on different sections

3.3.3 不同转速对全机气动特性影响

图13和图14分别给出了迎角为0°，螺旋桨转速为400r／min～900r／min范围内全机升力和俯仰力矩系数的变化情况。可以看出，随着螺旋桨转速的增大，升力系数明显增加，俯仰力矩系数减小，低头力矩越大。

图15给出了非定常方法计算拉力系数随螺旋桨转速的变化曲线。随着转速的增加，螺旋桨拉力系数越大，且在400～900r／min范围内呈线性变化。

图13 不同转速下升力系数的变化情况Fig.13 Lift coefficient vs.rotational speeds

图14 不同转速下俯仰力矩系数的变化情况Fig.14 Pitch moment coefficient vs.rotational speeds

图15 拉力系数随转速变化曲线Fig.15 Thrust coefficient vs.rotational speeds

4 结束语

本文采用曲线坐标系下的非定常Euler方程，综合运用动态拼接结构网格技术和双时间步推进方法，模拟了某双发螺旋桨飞机在螺旋桨转动时的非定常效应。为了保证计算模拟的高效性，引入多重网格技术和大规模并行计算技术，达到了快速、实用地模拟螺旋桨转动时的非定常效应。

通过单独螺旋桨算例，初步验证了非定常方法的可行性；通过模拟不同迎角和螺旋桨转速下全机的气动特性，分析了滑流对流场的影响情况，计算结果表明：螺旋桨滑流对全机流场和气动特性影响明显，并使得升力系数增大，在工程上具有一定的应用前景。

由于本文控制方程采用非定常Euler方程，忽略了粘性项，对阻力和扭矩计算能力不足。因而，下一步的研究工作将采用雷诺平均N－S方程进行计算模拟。

［1］LIU P Q.Research on air propeller theory and its application［M］.Beijing：BUAA Publishing Company，2006：42－49.（in Chinese）刘沛清.空气螺旋桨理论及其应用［M］.北京：北京航空航天大学出版社，2006：42－49.

［2］MOENS F，GARDAREIN P.Numerical simulation of the propeller／wing interactions for transport aircraft［R］.AIAA－2001－2404，2001.

［3］DOUGLAS HUNSAKER，DERYL SNYDER.A lifting－line approach to estimating propeller／wing interactions［R］.AIAA－2006－3466，2006.

［4］STRASH D J，et al.Analysis of propeller－induced aerodynamic effects［R］.AIAA－98－2414，1998.

［5］ZUO S H，YANG Y.Numerical simulation of propeller／highlift system interaction［J］.Aeronautical Computing Technique，2007，37（1）：54－57.（in Chinese）左岁寒，杨永.螺旋桨滑流对带后缘襟翼机翼气动特性影响的数值分析［J］.航空计算技术，2007，37（1）：54－57.

［6］LI B，LIANG D W，HUANG G P.Propeller slipstream effects on aerodynamic performance of turbo－prop airplane based on equivalent actuator disk model［J］.ACTA Aeronautica Et Astronautica Sinica，2008，29（4）：845－852.（in Chinese）李博，梁德旺，黄国平.基于等效盘模型的滑流对涡桨飞机气动性能的影响［J］.航空学报，2008，29（4）：845－852.

［7］LI B，JIN J，LIANG D W.Research on numerical simulation of propeller airplane［A］.2007Academic Annual Meeting of Chinese Aviation Society［C］.Shenzhen：Chinese Aviation Society，The Chinese Academy of Engineering，Mechanical and Vehicle Engineering，2007：Special of Aerodynamics 60.（in Chinese）李博，金君，梁德旺.螺旋桨飞机的全机流场数值模拟研究［A］.中国航空学会2007年学术年会［C］.深圳：中国航空学会，中国工程院机械与运载工程学部，2007：气动专题60.

［8］ROOSENBOOM E W M，STÜRMER Arne，SCHRÖDER Andreas.Advanced experimental and numerical validation and analysis of propeller slipstream flows［J］.Journal of Aircraft，2010，47（1）：284－291.

［9］REN X F，YANG S P，DUAN Z Y，et al.The analysis of the slipstream on the wing about aerodynamic and characteristics by multiple reference frame model［A］.The 14th Computational Fluid Dynamics Conference［C］.2009：582－585.（in Chinese）任晓峰，杨士普，段卓毅，等.基于多参考坐标系的螺旋桨滑流对机翼气动特性影响分析［C］.第十四届全国计算流体力学会议论文集，2009：582－585.

［10］ZHANG X L，ZHANG Y F.Research on interaction of propeller and high－lift system［J］.Aeronautical Computing Technique，2011，41（4）：1－3.（in Chinese）张小莉，张一帆.螺旋桨滑流对增升装置气动特性影响研究［J］.航空计算技术，2011，41（4）：1－3.

［11］XU H Y，YE Z Y.Numerical simulation of unsteady propeller slipstream［J］.Journal of Aerospace Power，2011，26（1）：148－153.（in Chinese）许和勇，叶正寅.螺旋桨非定常滑流数值模拟［J］.航空动力学报，2011，26（1）：148－153.

［12］WANG Y T，WANG G X，ZHANG Y L.Drag prediction of DLR－F6 configuration with TRIP2.0software［J］.ACTA Aerodynamica Sinica，2009，27（1）：108－113.（in Chinese）王运涛，王光学，张玉伦.采用TRIP2.0软件计算DLR－F6构型的阻力［J］.空气动力学学报，2009，27（1）：108－113.

［13］MENG D H.Numerical simulation and investigation of aircraft flutter［D］.Mianyang：China Aerodynamic Research and Development Center，2008.（in Chinese）孟德虹.飞行器颤振问题数值模拟技术研究［D］.绵阳：中国空气动力研究与发展中心，2008.

［14］HONG J W，LIANG X P，WANG G X，et al.The research of application for cross－multigrid technology［J］.ACTA Aerodynamica Sinica，2007，25（1）：45－54.（in Chinese）洪俊武，梁孝平，王光学，等.多重拼接网格技术应用研究［J］.空气动力学学报，2007，25（1）：45－54.

［15］VATSA V N，WEDAN B W.Development of a multigrid code for 3－D Navier－Stokes equations and its application to a grid refinement study［J］.Computers and Fluids，1990，18：391－403.

［16］YANG X C，WANG Y T，WANG G X.Numerical simulation of propeller slipstream based on equivalent actuator disk model［J］.Physics of Gases Theory and Applications，2012，7（3）：40－47.（in Chinese）杨小川，王运涛，王光学.基于等效盘的螺旋桨飞机数值模拟［J］.气体物理－理论与应用，2012，7（3）：40－47.