旋翼作为直升机类飞行器的主要气动部件, 其气动性能直接决定了整个飞行器的飞行性能。悬停和前飞是旋翼2个典型的工作状态, 其中悬停性能(例如悬停效率, FoM)可以作为评价旋翼气动性能的一个关键指标。随着计算机算力的快速增长, 直接从NS方程(空气动力学中“第一性原理”)出发的计算流体力学(CFD)方法逐步被用在旋翼非定常流场的模拟中。朱正、卢丛玲等[1-2]利用运动嵌套网格方法, 采用BL湍流模型及SA湍流模型对共轴刚性双旋翼非定常悬停流场进行了模拟, 分析了其流动干扰机理。张震宇等[3]利用多块结构化滑移网格, 采用k-ω SST两方程湍流模型对刚性旋翼大前进比前飞状态进行了模拟。美国航天航空学会(AIAA)自从2014年起每年召开一次旋翼悬停预测工作组会议(AIAA rotorcraft hover prediction workshop, HPW), 相继提供了3个旋翼标模[4-7]及相关的试验数据作为统一的研究对象, 包括S-76旋翼、PSP旋翼和HVAB旋翼。其中针对S-76旋翼[8-9]的讨论主要集中在全湍流假设下, 而HVAB旋翼模型目前缺乏正式的试验数据作为支撑。PSP旋翼是当前HPW会议中唯一有桨叶表面试验数据的计算标模, 从2018年起国外研究人员开始对其进行初步的数值模拟研究[10]。在基于RANS方程的CFD方法中, 转捩模拟是非常重要的一部分, CFD2030愿景[11]以及HPW愿景中均将边界层的转捩模拟作为重要的研究内容。对于固定翼类飞行器的二维、三维、高速、低速转捩模拟已经有相对成熟的方法和应用[12], 然而这些转捩模拟方法在旋翼模拟中的应用却比较少见。其原因在于：一方面相比于固定翼模拟, 旋翼模拟需要更大的计算代价, 因为涉及非定常计算、尾迹捕捉等; 另一方面转捩模拟对计算网格尺度和网格量的要求要远高于全湍流模拟。对于固定翼飞行器而言, 在一些状态下, 其气动、流场特性受到转捩过程的影响，会有比较大的变化, 所以对于这类飞行器在设计、评估过程中必须考虑转捩的影响。而对于旋翼来说, 人们需要明确是否也存在湍流模拟和转捩模拟结果差异大这一现象, 以确定在设计、评估过程中采用合适的模型和计算网格。
本文采用基于嵌套网格的自研RANS求解器, 对HPW提供的PSP旋翼标模进行了悬停流场的全湍流和转捩模拟, 并从积分气动力、桨叶表面压力分布、桨叶载荷分布、桨叶表面流线、尾迹流场等方面详细对比讨论了全湍流模拟与转捩模拟的差异, 为设计人员在评估时提供依据。同时也进一步验证了自研求解器在转捩模拟和旋翼模拟方面的可靠性。
1 数值方法 1.1 RANS方程本文自研求解器[13-17]基于动态结构化嵌套网格, 对雷诺平均NS方程(RANS方程)进行数值求解, 可压缩RANS方程如(1)式所示。
|
(1) |
式中: ρ为空气密度; ui是速度分量; p为压强; E为总能; H为总焓。
对于方程中无黏通量采用Roe迎风格式配合WENO插值实现5阶精度, 对于方程中的黏性通量采用2阶中心差分格式计算。使用近似因式分解方法进行方程的时间推进, 在非定常计算方面, 采用隐式双时间步长方法, 并在子迭代中使用多重网格和当地时间步长方法加速收敛。
物面采用无滑移边界条件, 远场采用基于黎曼不变量的无反射边界条件。
1.2 γ-Reθt-SA转捩模型γ-Reθt转捩模型最初是Langtry等人基于两方程k-ω SST湍流模型提出的, 之后Medida等[18-19]通过来流湍流度Iinf和速度梯度改变原始转捩模型中Reθt, Reθc和Flenght等参数计算方式, 将其与SA一方程湍流模型耦合, 提出γ-Reθt-SA转捩模型。其中主要控制方程的形式如(2)~(4)式所示。式中相关函数和参数的具体表达在文献[18-19]中有详细论述。
|
(2) |
|
(3) |
|
(4) |
相较于其他方法而言, 嵌套网格方法在网格处理方面有2个额外的步骤, 即挖洞(hole-cutting)和插值(interpolation)。其中, 挖洞是将侵入物面内的网格单元和网格点进行标记, 从而在计算中排除这些单元和点。插值是对于嵌套网格块边界处的网格点提供相应的插值模板单元和插值系数, 使得流场信息在相互嵌套的网格块之间进行传递。本文挖洞基于Objective X-ray方法, 插值采用Inverse Map与模板游走相结合的方法, 具体实现方法参考文献[13]。
2 求解器验证对于文中的转捩模型采用标模平板算例(S & K平板[20]、T3A平板[21-22])和S809翼型[23]进行验证, 对计算得到的表面摩阻因数、压力系数等结果与相应的试验结果进行比对, 说明求解器转捩模拟的能力。
3个算例的来流湍流度及所采用的计算网格规模在表 1中列出, 其中J为流向, K为物面法向。S809翼型计算网格在图 1中给出。
| 算例 | 计算状态 | 网格规模(J×K) |
| S & K | Iinf=0.03% | 201×81 |
| T3A | Iinf=3.3% | 201×81 |
| S809 | Iinf=0.02%, α=1°, Re=2×106 | 697×101 |
|
| 图 1 S809翼型网格 |
对平板标模算例, 对比了CFD计算与试验测得的表面摩阻因数, 而对S809翼型算例, 对比了翼型表面压力系数。对比结果如图 2~3所示, 3个验证算例的计算结果展示出了求解器在转捩模拟方面良好的预测能力和精度。
|
| 图 2 平板标模算例表面摩阻因数对比结果 |
|
| 图 3 S809翼型表面压力系数对比结果 |
计算模型采用HPW工作会议提供的旋翼标模算例-PSP旋翼, 其桨叶平面及截面翼型形状如图 4所示。此旋翼共包含4片桨叶, 半径为1.689 m, 参考弦长为0.138 m, 桨叶实度σ为0.103 3, 在0.95R处有30°的后掠角, 桨叶转速为1 150 r/min, 对应的桨尖马赫数为0.58, 基于桨叶平均弦长和桨尖速度的雷诺数为1.87×106。
|
| 图 4 PSP旋翼桨叶平面形状 |
该旋翼的风洞试验于2016年在NASA-LaRC旋翼测试中心[24]完成, 由于试验论文中并未提及相应的风洞试验湍流度, 在本文转捩计算中, 湍流度设为0.7%。
3.2 计算网格HPW工作组公开发布了PSP旋翼标模桨叶的几何文件和表面网格文件, 本文利用上述表面网格文件生成了相应计算网格, 包含对应于4片桨叶的贴体网格和笛卡尔背景网格。背景网格在桨叶附近进行加密处理。总网格点数为5.1×107, 其中贴体网格点数为4.4×107。贴体网格中桨叶表面法向首层网格距离为8.9×10-7 m。背景网格最密的网格尺度为1/10倍的桨叶弦长。远场网格距离旋转中心约为20倍的桨叶半径长度。图 5给出了桨叶贴体网格的3个截面和最密背景网格的1个截面, 以及相应的挖洞结果。
|
| 图 5 PSP旋翼标模嵌套网格 |
在本文具体计算中, 先将物理时间步长定为1°/步, 子迭代步数为2步, 共计算一个旋转周期, 得到合理的初始流场。并以所得计算结果为初始流场, 将物理时间步长缩短至0.5°/步, 子迭代步数增大至15步, 共计算6个旋转周期, 得到稳定的悬停流场。为了对比转捩过程带来的影响, 分别采用全湍流SA模型和γ-Reθt-SA转捩模型进行计算, 在等拉力系数的状态下进行2种模型计算结果的对比。在计算中通过调整旋翼总距实现上述等拉力系数的计算要求。
4.1 旋翼积分力在悬停状态下, 旋翼有2个最重要的气动力系数, 即拉力系数CT和扭矩系数CQ。按照HPW工作组的统一定义, (5)式给出了这2个系数的计算公式。
|
(5) |
式中: A为旋翼桨盘面积；R为桨叶半径; 最终悬停效率(FoM)为
在本文计算中, 重点对比同拉力系数下的扭矩系数和相应的悬停效率。所有的计算结果均为旋翼最后一个旋转周期内的平均值。表 2和图 6给出了悬停效率的计算结果与试验结果的对比, 图中试验值的误差限取自文献[24]中给出的最小值。该结果验证了本文所采用的求解器对于旋翼悬停气动力系数计算的可靠性, 转捩模拟结果和全湍模拟结果计算误差均在5%之内。
|
| 图 6 CFD预测旋翼气动性能与试验值对比 |
对于悬停效率这一评价旋翼气动效率的重要参数而言, 在低拉力状态下, 转捩模拟得到明显优于全湍模拟的结果(约高10%), 但是在大拉力状态下, 这一优势逐渐缩小(约高4%)。
4.2 截面压力分布及桨叶载荷分布桨叶载荷分布指桨叶拉力和扭矩沿着桨叶展向的分布, 根据HPW工作组建议, 旋翼桨叶截面压力系数Cp、拉力系数CdT和扭矩系数CdTor的定义在公式(6)中给出, 这几个系数均以当地旋转速度进行无量纲处理。
|
(6) |
式中: Ω为转速；r为当地半径；c为当地弦长。
4.2.1 截面压力分布沿桨叶展向从30%R处开始至80%R处结束, 截取了6组压力分布, 对比结果如图 7~8所示。其中图 7为小拉力状态下结果, 图 8为大拉力状态下结果。
|
| 图 7 小拉力状态桨叶截面压力分布对比 |
|
| 图 8 大拉力状态桨叶截面压力分布对比 |
相比而言, 全湍流模拟结果得到了比较平滑的桨叶表面压力分布, 而由于计入了流动转捩过程, 转捩模拟得到的表面压力分布在流动转捩区域附近出现了明显波动。在其余区域, 两者压力分布相似且重合度高。
通过桨叶截面上下表面压力的波动可以判断出桨叶转捩发生的区域, 经过对比与试验[24]中给出转捩位置保持一致。进一步证明了本文所采用的求解器在旋翼悬停自然转捩模拟上的可靠性。在小拉力状态下桨叶表面出现较大范围的层流流动, 而随着旋翼拉力的增大, 该层流范围逐步减小, 这与4.1节中旋翼积分扭矩系数变化趋势保持一致。
4.2.2 桨叶载荷分布图 9给出了不同拉力状态下, 全湍流以及转捩模型计算得到的桨叶载荷系数在桨叶展向上的分布, 其中拉力系数分布呈现出合理的类椭圆形分布。
|
| 图 9 桨叶载荷分布对比 |
虽然截面压力分布在转捩区域有明显的波动, 但是对于沿桨叶展向的拉力系数, 除了桨叶根部小拉力区域, 2种模拟方法并没有表现出很大的不同。而对于沿桨叶展向扭矩系数来说, 由于转捩模拟中部分区域存在明显的层流流动, 在这部分区域转捩模型预测得到的扭矩系数要低于全湍模型预测的值。在大拉力状态下, 由于表面层流范围的缩小, 全湍模拟和转捩模拟得到的扭矩系数分布也趋于一致。
4.3 桨叶表面流线及旋翼尾迹流场形态 4.3.1 桨叶表面流线图 10~11分别显示出不同拉力状态下桨叶上、下表面的表面流线, 可看出全湍状态下桨叶表面流场变化平缓, 但是转捩模拟得到了明显不同的桨叶表面流线, 其表面流线变化明显, 在桨叶内侧上表面发生了明显的流动分离现象, 这一现象在全湍模拟中并没有出现。表面流线的明显差异也进一步说明了4.2.1节中2种不同模拟方法得到的差异明显的桨叶截面压力分布。
|
| 图 10 小拉力状态桨叶表面流线对比 |
|
| 图 11 大拉力状态桨叶表面流线对比 |
通过Q判据等值面, 图 12给出了不同拉力状态下旋翼桨盘下方的尾迹桨尖涡分布, 利用涡量大小在对应的Q等值面上进行相应的染色。根据上述计算结果, 转捩过程引发的桨叶表面分离对桨尖涡的影响较小, 但是对于桨盘下方小于桨盘半径处的尾迹结构有比较明显的影响。同时, 随着拉力的增大, 桨尖涡的结构形式出现了一些变化, 尾迹流场中开始出现围绕着桨尖涡的二次涡结构, 这一现象在全湍和转捩模拟中均有所表现。
|
| 图 12 旋翼尾迹流场对比 |
本文利用自研RANS求解器, 对AIAA HPW工作组提供的PSP旋翼标模, 在悬停状态下, 分别进行了全湍流模拟和转捩模拟。
计算结果首先表明本文CFD模拟得到的旋翼拉力-悬停效率曲线与试验重合度高, 计算误差在5%之内, 并且转捩模拟中, CFD方法能正确模拟出桨叶表面流动的转捩位置。由此, 证明了本文求解器在旋翼悬停流场模拟方面的可靠性。
其次, 通过全湍与转捩模拟结果的对比发现, 转捩过程对桨叶截面压力分布、沿桨叶展向扭矩分布以及桨叶表面流线均有显著影响, 而对沿桨叶展向拉力分布影响较弱。转捩过程给桨叶截面压力分布带来了明显的波动以及桨叶表面流场的明显分离。相比于全湍结果, 转捩模拟在层流区域得到了更低的扭矩系数, 从而带来了悬停效率4%~10%的增量。
最后, 对于旋翼桨盘下方桨尖涡, 转捩过程没有明显影响, 但是由于转捩过程带来的流动分离, 对于内侧桨叶下方尾迹区域, 转捩过程带来了明显的影响。
| [1] |
朱正, 招启军, 李鹏. 悬停状态共轴刚性双旋翼非定常流动干扰机理优先出版[J]. 航空学报, 2016, 37(2): 568-578.
ZHU Zheng, ZHAO Qijun, LI Peng. Unsteady flow interaction mechanism of coaxial rigid rotors in hover[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(2): 568-578. (in Chinese) |
| [2] |
卢丛玲, 史勇杰, 徐国华, 等. 共轴刚性旋翼悬停状态气动干扰机理[J]. 南京航空航天大学学报, 2019, 51(2): 201-207.
LU Congling, SHI Yongjie, XU Guohua, et al. Research on aerodynamic interaction mechanism of rigid coaxial rotor in hover[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2019, 51(2): 201-207. (in Chinese) |
| [3] |
张震宇, 钱耀如, 王同光, 等. 基于非定常方程的刚性旋翼流场分析[J]. 南京航空航天大学学报, 2019, 51(2): 238-243.
ZHANG Zhenyu, QIAN Yaoru, WANG Tongguang, et al. Analysis on flow-field around helicopter rigid rotor blades based on unsteady RANS equations[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2019, 51(2): 238-243. (in Chinese) |
| [4] | HARIHARAN N S, EGOLF T A, SANKAR L N. Simulation of rotor in hover: current state, challenges and standardized evaluation[C]//52nd AIAA Aerospace Sciences Meeting, 2014 |
| [5] | HARIHARAN N S, NARDUCCI R P, REED E, et al. AIAA standardized hover simulation: hover performance prediction status and outstanding issues[C]//55th AIAA Aerospace Sciences Meeting, 2017 |
| [6] | PARWANI A, CODER J G. Effect of laminar-turbulent transition modeling on PSP rotor hover predictions[C]//56th AIAA Aerospace Sciences Meeting, 2018 |
| [7] | SINGH R, CORLE E, JAIN R, et al. Computation and quantification of uncertainty in predictions of HVAB rotor performance in hover[C]//AIAA Scitech 2019 Forum, 2019 |
| [8] | JAIN R. Sensitivity study of high-fidelity hover predictions on the sikorsky S-76 rotor[J]. Journal of Aircraft, 2017, 55(1): 1-11. |
| [9] | JAIN R. A comparison of CFD hover predictions for the sikorsky S-76 rotor[C]//54th AIAA Aerospace Sciences Meeting, 2016 |
| [10] | WEISS A, GARDNER A D, SCHWERMER T, et al. On the effect of rotational forces on rotor blade boundary-layer transition[J]. AIAA Journal, 2019, 57(1): 252-266. DOI:10.2514/1.J057036 |
| [11] | SLOTNICK J, KHODADOUST A, ALONSO J, et al. CFD vision 2030 study: a path to revolutionary computational aerosciences[R]. NASA/CR-2014-218178 |
| [12] |
杜一鸣. 涡黏性湍流模型修正与三维边界层转捩预测方法研究[D]. 西安: 西北工业大学, 2021 DU Yiming. Research on modification of RANS eddy-viscosity turbulence model and prediction method of three-dimensional boundary-layer transition[D]. Xi'an: Northwestern Polytechnical University, 2021 (in Chinese) |
| [13] | PANG Chao, GAO Zhenghong, YANG Hua, et al. An efficient grid assembling method in unsteady dynamic motion simulation using overset grid[J]. Aerospace Science and Technology, 2020(110): 106450. |
| [14] | PANG Chao, GAO Zhenghong. CFD simulation of unsteady interaction between rotor and fuselage for canard rotor wing in hovering using overset grids[C]//9th Asia Conference on Mechanical and Aerospace Engineering, 2018 |
| [15] | PANG Chao, GAO Zhenghong. Numerical investigation of aerodynamic performance for a canard rotor wing aircraft in late conversion mode[C]//Asia Pacific International Symposium on Aerospace Technology, Engineers, Australia, 2019 |
| [16] | PANG Chao, YANG Hua, GAO Zhenghong, et al. Enhanced adaptive mesh refinement method using advanced vortex identification sensors in wake flow[J]. Aerospace Science and Technology, 2021(115): 106796. |
| [17] |
杜一鸣, 庞超, 舒博文, 等. 基于嵌套网格的CHN-T1标模数值模拟[J]. 空气动力学学报, 2019, 37(2): 280-290.
DU Yiming, PANG Chao, SHU Bowen, et al. Numerical simulation of the standard model CHN-T1 based on overset grid[J]. Acta Aerodynamica Sinica, 2019, 37(2): 280-290. (in Chinese) |
| [18] | MEDIDA S, BAEDER J. Application of the correlation-based γ-Reθt transition model to the spalart-allmaras turbulence model[C]//20th AIAA Computational Fluid Dynamics Conference, 2011 |
| [19] | SHIVAJI M. Correlation-based transition modeling for external aerodynamic flows[D]. Maryland: University of Maryland College Park, 2014 |
| [20] | SCHUBAUER G B, KLEBANOFF P S. Contributions on the mechanics of boundary-layer transition[R]. NACA-TR-1289, 1955 |
| [21] | SAVILL A M. Some recent progress in the turbulence modelling of bypass transition//Near-Wall Turbulent Flows[M]. Netherlands: Elsevier, 1993: 829-848. |
| [22] | SAVILL A M. Evaluating turbulence model predictions of transition//Advances in Turbulence IV[M]. Dordrecht: Springer, 1993: 555-562. |
| [23] | SOMERS D M. Design and experimental results for the S809 airfoil[R]. NREL/SR-440-6918, 1997 |
| [24] | OVERMEYER A D, MARTIN P B. Measured boundary layer transition and rotor hover performance at model scale[C]//55th AIAA Aerospace Sciences Meeting, 2017 |
