1 数值计算方法
1.1 旋翼气动模型
流动控制方程采用任意拉格朗日欧拉方法描述的可压缩N⁃S方程，允许网格的任意运动和变形，在绝对坐标系下的积分形式如下
式中：
空间格式采用Ro
e[13] 通量差分分裂方法，单元面上的原始变量用三阶MUSCL插值得到，并加以连续可微限制器进行限制。湍流模型为两方程k⁃ω SST模型，为了加速流场收敛，计算采用了三重W循环的多重网格方法。应用双时间法模拟干扰流场的非定常流动现象，同时，在伪时间方向上使用隐式LU⁃SGS格式[14] 进行时间推进以提高流场的计算效率。在复杂外形上生成结构网格，常见的网格分区处理技术包括多块对接网格和重叠网格，前者的各子块之间对接面相互重合，对接单元一一对应；后者的网格块之间包括相互重叠的部分，通过插值进行流场信息传递。本文采用多块重叠网格方法模拟上、下旋翼桨叶的运动。该网格系统由两部分组成：第一部分为围绕旋翼生成的桨叶网格，第二部分为静止的背景网格，计算中生成的网格系统将在后文给出。
1.2 旋翼配平方法
共轴旋翼悬停时要求双旋翼扭矩平衡，需要进行总距配平。共轴刚性旋翼单副旋翼的桨距随方位角的变化可表示为
式中：
操纵量为
目标量为
式中：CT为总拉力系数；CQ为总扭矩系数；CL为滚转力矩系数；CM为俯仰力矩系数；上标u，l分别表示上旋翼、下旋翼。实际计算中，升力偏置量可表示为滚转力矩与升力的比值，即
操纵量和目标量间的关系可用非线性方程组表示
将式（6）进行泰勒级数展开至一阶并化简，可得目标量和操纵量的迭代关系式为
式中：
本文发展了基于差量法的共轴刚性旋翼高效配平方法。简单模型主要包括动量叶素理论（Blade element momentum theory, BEMT）和用来求解诱导速度的入流模型。本文用到的主要有均匀入流、Drees入流等入流模型。
基于差量法的配平计算由简单气动力模型完成，CFD求解器只对气动力进行修正，不直接参与配平迭代。图1给出了本文配平方法的计算流程，具体步骤如下：
（1）采用简单气动模型，根据目标量进行配平获得初始操纵量。
（2）将上述操纵量代入CFD求解器计算共轴旋翼旋转一圈，得到旋翼气动力。
（3）根据CFD计算的旋翼气动力对简单模型计算的气动力进行修正。
（4）给定新的简单气动模型目标量
（5）重复步骤（2）至（4），直至收敛。
2 结果分析
2.1 Harrington⁃2共轴旋翼计算
Harrington⁃2共轴旋
翼[15] 由两片刚性矩形桨叶组成，展弦比8.33，采用NACA 00XX翼型，翼型厚度沿展向线性变化，在0.2R位置厚度为27.5%，在R位置相对厚度15%。桨尖马赫数0.352，桨尖雷诺数3.5×106 。本文计算中采用结构重叠网格技术，网格单元总量约6 236万。该网格由两部分组成：第一部分为围绕旋翼生成的桨叶网格，第二部分为静止的背景网格，如图2所示。考虑到在桨叶旋转运动过程中，桨叶尖部、根部分别存在桨尖涡和桨根涡，对桨叶气动特性影响显著。这种三维效应的存在使得桨叶尖部和根部的网格质量至关重要。为了保证网格在桨叶端部能保持良好的正交性及物面间距，该桨叶子网格采用了C⁃O型拓扑结构，单个桨叶的网格量为309万，四片桨叶网格量为1 236万，第一层网格到桨叶表面的距离为0.005 mm。图3给出了计算拉力系数CT与扭矩系数CQ变化曲线与试验结果的对比。图中可见，两者吻合较好，表明建立方法的有效性，可用于共轴刚性旋翼气动干扰的分析。
图3 Harrington⁃2共轴旋翼气动性能计算与试验值对比
Fig.3 Comparison of calculated results with experimental data for Harrington⁃2 coaxial rotor
图4给出了的Harrington⁃2共轴双旋翼空间涡量分布。可见，本文计算能够分辨的桨尖涡尾迹角约为650°，在当前的网格密度下，模拟得到的桨尖涡耗散较小。桨叶在拖出桨尖涡同时拖出了桨根涡和尾迹涡层，其中尾迹涡层达到了2层以上。图3也反映了桨尖涡强度随高度不断衰减以及向内发展的特征，符合实际物理规律。这些充分说明本文采用的数值计算方法对悬停状态下的涡尾迹特征具有较高的捕捉精度，可以有效模拟直升机旋翼飞行状态下的复杂流场及其细节特征。
2.2 2 m直径共轴刚性旋翼计算
该2 m直径共轴模型旋翼模型的桨叶扭转角为-12°线性扭转。桨叶片数为4片，桨叶基本弦长0.07 m，旋翼直径为2 m。计算网格采用重叠网格方法。在桨叶周围生成贴体O型网格，网格维数为61×105×305（法向×展向×周向），每片桨叶的网格单元总数为195万，如图5所示。为了更好地模拟旋翼尾迹，背景网格在旋翼平面及旋翼下方进行了加密，同时在展向方向上对应桨尖和桨根位置也进行了加密，这些加密区域的网格密度为桨叶弦长的0.06倍。
表1，2分别给出了上旋翼和下旋翼拉力系数CT和扭矩Mk的计算结果与试验的对比值，值得一提的是为了方便对比，扭矩给出的是量纲值。可以看出，在上下旋翼保持相同总矩的情况下，旋翼拉力要小于上旋翼，这意味着下旋翼的气动效率明显要低于上旋翼，上、下旋翼的拉力占总拉力比率基本保持在58%和42%左右。从表中也可看出，随着总矩角θ0的增大，拉力系数和扭矩都在增大。
表1 共轴模型旋翼的上旋翼悬停计算结果
Tab.1 Computational results of the upper rotor for coaxial rotor model in hover
试验值CT 计算值 误差/
%
试验值
|Mk|/ (N∙m)
计算值/ (N∙m) 误差/
%
0.007 42 0.007 79 4.9 -53.84 -56.055 4.1 0.011 00 0.010 70 3.0 -72.34 -74.360 2.8 0.014 70 0.014 70 0.3 -104.52 -108.70 4.0 表2 共轴模型旋翼的下旋翼悬停计算结果
Tab.2 Computational results of the lower rotor for coaxial rotor model in hover
试验值CT 计算值 误差/
%
试验值
|Mk|/ (N∙m)
计算值/
(N∙m)
误差/
%
0.005 09 0.005 30 4.0 43.18 42.35 2.0 0.008 20 0.008 56 4.4 67.93 69.65 2.5 0.010 90 0.011 30 4.3 94.12 97.70 3.8 图6给出的是3个工况下旋翼的物面压力云图，图中桨叶位于正x轴方向时周向角为0°，上旋翼旋转方向为逆时针方向，下旋翼旋转方向为逆时针方向。可见，拉力系数越大的工况在桨尖处的负压值愈大且范围更宽，说明拉力的提升主要来自于桨尖侧迎角的增大。
图7计算了模型旋翼悬停状态，尾迹空间Q旋涡判定准则的等值面图。图中显示的是无量纲 Q值大于0. 5的涡量分布。由于双旋翼有一定的轴向间距，从上旋翼桨尖处脱落的集中涡呈螺旋状向下运动，在短暂的时间后会与下旋翼桨叶发生碰撞，从而出现“上/下旋翼桨⁃涡干扰”现象。这类干扰属于垂直桨⁃涡干扰，即干扰过程中尾迹涡轴线与桨叶剖面垂直，另外随着总矩角的增大，尾迹涡的强度和范围都在增加，且尾迹涡的收缩向内发展的速度也在加快。
表3给出了共轴模型旋翼前飞算例设置，其中桨尖马赫数为0.578，前进比为0.2。α为旋翼轴倾角，前倾为正。表中所列为试验中进行前飞垂向力和阻力配平状态下的操纵量，计算中采用与试验相同操纵量。此时双旋翼采用相同的总距
表3 共轴模型旋翼前飞算例设置
Tab.3 Settings of coaxial rotor model cases in forward flight
算例 α/(°) Case 1 4 6.54 -4.22 0.17 1.22 Case 2 4 6.67 -4.34 0.09 1.33 Case 3 -4 7.69 -3.81 1.15 2.14 Case 4 -4 8.93 -3.88 1.16 2.30 表4,5分别给出了CT和Mk的计算结果与试验的对比值，扭矩仍然给出的是量纲值。可见，旋翼拉力计算值与实验值误差在5%以内。
表 4 共轴模型旋翼的上旋翼前飞计算结果
Tab.4 Computational results of the upper rotor of coaxial rotor model in forward flight
试验值CT 计算值 误差/
%
试验值
|Mk|/ (N∙m)
计算值/ (N∙m) 误差/
%
0.012 1 0.011 8 2.1 -66.38 -69.50 4.7 0.012 1 0.011 9 1.7 -69.25 -72.37 4.5 0.011 9 0.011 6 3.3 -92.51 -96.30 4.1 0.012 1 0.011 9 1.0 -94.03 -97.14 3.3 表 5 共轴模型旋翼的下旋翼前飞计算结果
Tab.5 Computational results of the lower rotor of coaxial rotor model in forward flight
试验值
CT
计算值 误差/
%
试验值
|Mk|/ (N∙m)
计算值/
(N∙m)
误差/
%
0.012 3 0.011 8 4.6 42.86 44.92 4.8 0.012 4 0.011 9 4.0 44.61 46.39 4.0 0.012 5 0.011 3 4.6 65.60 68.43 4.3 0.012 5 0.012 4 1.5 63.54 65.44 3.0 图8是计算的前飞状态共轴旋翼尾迹空间Q旋涡判定准则的等值面图，显示的是无量纲Q值大于0.5的涡量分布。从图中可见，涡量比较集中的区域包括3个方面：一是从桨叶桨尖拖出的集中涡，二是桨叶后面的尾迹层，三是物面边界层区域。由桨尖拖出的集中涡和桨叶尾迹贴近后继桨叶的表面经过，对后继桨叶流场造成较大的干扰，比如图中在周向角180°位置的桨叶上就可以清楚地看到有集中涡经过。尾涡捕捉对网格密度的要求较高，空间网格又不可能做到处处加密，在尾涡进入网格稀疏的区域时，涡的耗散加剧，从图中可以看到桨尖涡在周向上最多发展了180°左右。在旋翼左右两侧，由桨尖涡卷起形成了集中涡，呈不对称分布。另外，虽然4个工况下的拉力系数相差不大，但总距越大，对应的上下旋翼扭矩值越大，则涡量强度越大，脱出的尾迹愈明显。
3 结论
本文先对Harrriton⁃2共轴双旋翼的悬停状态及气动特性进行了计算，得到了不同总距下拉力系数、扭矩系数的气动数据并与试验进行了对比，证明了程序的可靠性。然后以某2 m共轴双旋翼模型为研究对象，采用CFD方法分别对旋翼悬停和前飞状态进行了数值模拟，获得了旋翼悬停性能和典型工况下的前飞气动特性，得到以下结论：
（1）基于多重嵌套网格和CFD技术的共轴双旋翼气动干扰特性分析方法，在典型状态，共轴刚性旋翼拉力系数的计算结果与试验值误差在3%以内，扭矩系数的误差基本在5%以内。所建立方法能够有效计算共轴刚性旋翼典型干扰状态下的气动性能。
（2）悬停状态从上旋翼桨尖处脱落的集中涡呈螺旋状向下运动，会与下旋翼桨叶发生碰撞，从而出现“上/下旋翼桨⁃涡干扰”现象。这类干扰属于垂直桨⁃涡干扰，干扰过程中尾迹涡轴线近似与桨叶剖面垂直。
（3）所建立的数值计算方法对旋翼涡尾迹特征具有较高的捕捉精度，尚可有效模拟共轴刚性旋翼悬停和前飞下的复杂流场及其细节特征。
参考文献
- 1
BAGAI A. Aerodynamic design of the X2TM technology demonstrator main rotor blade[C]// Proceeding of the 64th Annual Forum of American Helicopter Society. Montreal, Canada:AHS, 2008.
- 2
COLEMAN C P. A survey of theoretical and experimental coaxial rotor aerodynamic research[C]//Proceedings of the 19th European Rotorcraft Forum. Italy:
[s.n.], 1993.
- 3
LEISHMAN J G, ANANTHAN S. Aerodynamic optimization of a coaxial proprotor[C]//Proceeding of the 62nd Annual Forum of the American Helicopter Society. [S.l.]:AHS, 2006.
- 4
周国仪, 胡继忠, 曹义华. 共轴式直升机双旋翼载荷计算模型研究[J]. 航空动力学报, 2003, 18(3): 343⁃347.
ZHOU Guoyi，HU Jizhong，CAO Yihua. Mathematical model for twin rotor loads of a coaxial helicopter[J]. Journal of Aerospace Power, 2003, 18(3): 343⁃347.
- 5
LEISHMAN J G, BHAGWAT M J, BAGAI A. Free⁃vortex filament methods for the analysis of helicopter rotor wakes[J]. Journal of Aircraft, 2002, 39(5): 759⁃775.
- 6
GRIFFITHS D A, LEISHMAN J G. A study of dual⁃rotor interference and ground effect using a free⁃vortex wake model[C]//Proceedings of the 58th Annual Forum and Technology Display of the American Helicopter Society International. [S.l.]: AHS, 2002.
- 7
黄水林, 徐国华, 李春华. 基于自由尾迹方法的共轴式双旋翼流场分析[J]. 南京航空航天大学学报, 2008, 40(6): 721⁃726 .
HUANG Shuilin, XU Guohua, LI Chunhua. Flow field analysis of coaxial twin rotors based on free wake[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2008, 40(6): 721⁃726 .
- 8
童自力, 孙茂. 共轴式双旋翼气动力特性的计算研究[J]. 航空学报, 1999,20(4): 348⁃350.
TONG Zili, SUN Mao. Navier⁃Stokes analysis of the aerodynamic properties of coaxial rotors[J]. Acta Aero⁃nautica et Astronautica Sinica, 1999,20(4): 348⁃350.
- 9
LAKSHMINARAYAN V K, DURAISAMY K, BAEDER J D. Computational investigation of coaxial rotor aerodynamics in hover[C]//Proceeding of the 63rd Annual Forum of the American Helicopter Society. Virginia Beach: AHS, 2007.
- 10
LAKSHMINARAYAN V K, BADER J D. Computational investigation of small scale coaxial rotor aerodynamics in hover[R]. AIAA 2009⁃1069, 2009.
- 11
许和勇, 叶正寅. 悬停共轴双旋翼干扰流动数值模拟[J]．航空动力学报，2011, 26(2): 453⁃457.
XU Heyong, YE Zhengyin. Numerical simulation of interaction unsteady flows around coaxial rotors in hover[J]. Journal of Aerospace Power, 2011,26(2): 453⁃457.
- 12
叶靓, 徐国华. 共轴式双旋翼悬停流场和气动力的CFD计算[J]．空气动力学学报，2012, 30(4): 437⁃442.
YE Liang, XU Guohua. Calculation on flow field and aerodynamic force of coaxial rotors in hover with CFD method [J]．Acta Aerondynamica Sinica,2012, 30(4): 437⁃442.
- 13
ROE P L. Approximate riemann solvers, parameter vectors, and difference schemes[J]. Journal of Computational Physics, 1981, 43(2): 357⁃372.
- 14
LUO H, BAUM J D. A fast, matrix⁃free implicit method for computing low Mach number flows on unstructured grids[R]. AIAA 99⁃3315, 1999.
- 15
HARRINGTON R D. Full⁃scale⁃tunnel investigation of the static⁃thrust performance of a coaxial helicopter rotor[R]. NACA TN⁃2318, 1951.
- 1
摘要
建立了一个适用于共轴刚性旋翼气动特性分析的数值模拟方法。该方法采用任意拉格朗日欧拉方法(Arbitrary Lagrange Euler ,ALE)描述的可压缩Navier⁃Stokes（N⁃S）方程求解流场，采用低数值耗散的Roe格式进行空间离散；使用多重嵌套网格方法以模拟双旋翼的运动。针对共轴刚性旋翼配平，引入“差量修正”策略解决了传统配平中雅克比矩阵计算复杂的问题。首先，对Harrington⁃2共轴双旋翼的悬停气动性能进行了计算，然后，对某2 m直径共轴双旋翼的悬停及前飞状态进行了计算，并与试验值进行了对比。结果表明：在典型状态下拉力系数的计算结果与试验值误差在3%以内，扭矩系数的计算结果与试验值误差基本在5%以内；所采用的数值计算方法对旋翼涡尾迹特征具有较高的捕捉精度，可以有效模拟共轴刚性旋翼悬停和小速度前飞下的复杂流场及其细节特征。
Abstract
A numerical method is developed to predict the aerodynamic characteristics of rigid coaxial rotor. The method uses the compressible Navier⁃Stokes (N⁃S) equation described by arbitrary Lagrange Euler (ALE) to solve the flow field. A low numerical dissipative Roe format for spatial dispersion is adopted. Multiple overset grid method is used to simulate the motion of the coaxial rotor. A delta trim method is implemented to simplify the calculation of Jacobin matrix for rigid coaxial rotor trim. First, the aerodynamic performance of the Harrington⁃2 coaxial rotor in hover is calculated. Then, the hover and forward flight states of a 2 m diameter coaxial rotor are calculated and compared with the experimental values. The results show that, the error of the calculation result and the test value is less than 3% under the typical state, and the error of the torque coefficient is less than 5%. The numerical method has a high capture precision for the vortex wake characteristics of the rigid coaxial rotor, and can effectively simulate the detailed characteristics of the complex flow fluid in hover and the low⁃speed forward flight states.
共轴刚性旋翼复合式直升机是未来高速直升机发展的主要构型之一，该直升机最大特点是采用了“前行桨叶概念旋翼”，也称共轴刚性旋
虽然目前关于共轴刚性旋翼气动问题的研究还比较少，但是在常规共轴旋翼气动特性方面，已开展了诸多研究。Colema
鉴于此，本文开展了适用于共轴刚性旋翼气动特性分析的数值模拟方法研究。主要研究工作包括：建立了适合共轴刚性旋翼非定常干扰流场模拟的多重嵌套网格系统和CFD分析方法；针对共轴刚性旋翼的操纵特征，发展了一种基于差量法的高效配平策略；对共轴刚性旋翼的悬停和前飞状态进行计算，并与试验值进行了对比，验证了所建立方法的有效性。
