2. 西北工业大学 航空学院, 陕西 西安 710072
旋翼桨尖涡的结构特征相对紧凑且在输运过程中伴随着黏性扩散及桨涡干扰等现象，想要对其进行准确的数值模拟比较困难。早在1956年Gray[1]就描述了旋翼尾迹的主要流场结构，主要包括旋翼桨尖涡和尾涡面2部分，其中对整个流场影响最大占主导地位的是旋翼桨尖涡，但是一直到今天旋翼桨尖涡的高精度数值模拟仍然是一个挑战，是旋翼CFD数值模拟中的关键技术之一。悬停状态的旋翼流场具有代表性，是研究旋翼其他状态流场的基础，针对悬停状态的旋翼桨尖涡开展先进数值模拟方法研究对于深入研究直升机旋翼绕流机理、桨涡干扰、旋翼/机身干扰以及旋翼气动噪声特性评估等有着重要的意义。
通过求解非定常雷诺平均N-S(RANS)方程的方法对直升机旋翼流场进行数值分析，能够比较准确地模拟旋翼尾涡的形成和输运，尾涡作为解的一部分得到，不需要附加尾涡模型。尽管该方法在旋翼桨叶气动力计算方面取得很好的效果，但是计算得到的旋翼桨尖涡形态还存在一定的偏差，一方面是传统的数值格式所隐含的数值耗散过大，导致旋翼桨尖涡在计算过程中被耗散减弱，从而引起计算对桨尖涡的描述存在非物理失真；另一方面是旋翼流场以桨尖涡为主要特征，RANS方法对于这种复杂湍流流动的模拟能力比较有限[2]。得益于近些年RANS/LES混合方法的快速发展，国外研究者已开始尝试采用RANS/LES混合方法对旋翼绕流流场进行数值模拟研究并且取得了一些突破[3-5]；国内关于旋翼桨尖涡RANS/LES混合方法数值模拟研究的公开文献还比较少。通过对国外旋翼RANS/LES混合方法的研究文献调研发现，大部分研究都是针对悬停状态的旋翼流场展开且采用的是DES方法，但是由于DES方法存在模型应力损失问题(modeled stress depletion)[6]，因此需要采用更为可靠的RANS/LES混合方法对旋翼流场进行数值模拟研究。
本文基于结构网格动态面搭接技术、五阶Roe-WENO格式和IDDES方法开发了一套用于悬停状态旋翼桨尖涡高精度数值模拟的计算分析程序，并分别采用RANS方法和IDDES(improvement of delayed detached-eddy simulation)方法针对悬停状态的旋翼的桨尖涡进行了数值模拟对比分析，并从旋翼桨尖涡的涡强度、涡核位置等方面对旋翼桨尖涡结构的细节特征进行了比较系统的分析，可为旋翼绕流机理研究和旋翼优化设计提供理论基础与技术参考。
1 数值分析方法 1.1 控制方程本文对悬停状态旋翼流场的数值分析是基于惯性笛卡尔坐标系开展的，在惯性笛卡尔坐标下，积分形式的非定常N-S方程组可写作：
|
(1) |
式中，Ω为控制体, S为控制体边界, n为S微元的外法向单位矢量。采用有限体积法离散求解上述方程, 黏性项采用二阶中心差分格式离散, 无黏项采用Roe格式离散, 流场分析基于全湍流假设, 采用k-wSST两方程湍流模型进行湍流数值模拟, 时间推进采用双时间迭代方法。
1.2 五阶Roe-WENO格式采用高阶格式(例如基于高阶WENO插值的格式)来有效降低计算方法所固有的数值耗散, 可以提高旋翼桨尖涡数值模拟的精度。对无黏通量项离散选用基于迎风型通量差分分裂的Roe格式。针对Roe格式构造了基于改进的WENO方法对交界面两边状态参量进行重构的Roe-WENO五阶格式。
1.3 结构网格动态面搭接技术参照国外有关旋翼流场RANS/LES混合方法的研究文献, 将围绕旋翼悬停流场的桨尖涡进行RANS/LES混合方法数值模拟研究。一方面旋翼悬停流场具有代表性, 是研究旋翼其他状态流场的基础。另一方面, 旋翼前飞状态下桨叶运动复杂, 包括桨叶的挥舞、变距等运动, 对这些运动处理不当会给数值方法的研究引入额外的干扰。
同时, 采用RANS/LES方法对流场进行数值模拟时, 在流动关注区域内要求计算网格最好为在3个维度方向上尺寸相同的各向同性立方体网格单元。进行计算网格设计时, 若采用多块点对接网格, 对流动关注区域计算网格的加密会延伸至整个流场计算空间从而引起计算网格量的极大增加, 造成计算资源的浪费。面搭接网格与点对接网格相比, 搭接面左右两侧网格节点数量和分布规律可以不相同, 即网格节点分布可以在搭接面处存在不连续, 这样一来可以大大减少不必要的计算量, 将计算网格尽可能集中在流动重点关注区域。
基于以上考虑, 本文采用结构网格动态面搭接技术来实现旋翼悬停流场的RANS/LES混合方法数值模拟。采用空间通量守恒插值法实现搭接面两侧的通量传递, 具体过程详见参考文献[7]。
1.4 IDDES方法IDDES方法作为RANS/LES混合方法的一种, 在流场中不同的区域分别采用不同的计算方法。在流场中漩涡运动所主导的区域采用DDES方法进行求解, 而在近物面边界层内则采用LES方法的壁面函数模型(WMLES)以较好的模拟边界层的速度型分布, 并且消除DDES方法在附面层内计算所引起的“对数层不符”(log-layer mismatch)问题。IDDES方法结合了DDES方法和WMLES方法的优点, 与其他DES类方法[8-9]相比对于小分离流动的模拟能力也较强。
与DDES方法相比, IDDES方法主要在2个方面进行了修改[10]。第一个是对网格尺度进行了重新的定义, 第二个是构造了新的RANS方法和LES方法的混合函数。下面对这2个方面分别进行介绍。
IDDES方法中对网格尺度Δ的定义不仅考虑了网格单元尺度的大小, 而且引入了物面距离的影响, 其定义如下:
|
(7) |
式中，dw为距物面的距离, hmax为网格单元在3个维度方向上的最大尺度, hwn为在物面法向方向上的当地网格尺度, Cw为常数取为0.15。
在IDDES方法中重新构造了新的湍流长度尺度lIDDES, 其定义如下:
|
(8) |
将上式中所定义的湍流长度尺度lIDDES替代k-w SST湍流模型中的湍流长度尺度就形成了IDDES混合模型。
(8) 式中混合函数fhyb包含了DDES分支和WMLES分支, 其构造形式如下:
|
(9) |
式中，fd是DDES方法中的延迟函数
|
式中，函数fstep只在WMLES模型被调用时起作用, 它使得在边界层内RANS方法能够快速向LES进行转换。其构造形式如下
|
(10) |
式中，参数α的定义为α=0.25-dw/hmax, fstep的构造能够使得在0.5275 < dw/hmax < 1.0区域内RANS方法快速向LES进行转换。
|
(11) |
式中，fhill的定义如下
|
在(11)式中函数famp的形式如下
|
(12) |
式中，ft, fl分别为:
|
k, ct为常数, 分别取为k=0.41, ct=1.87, =5.0。
2 算例验证串列双圆柱模型是BANC-Ⅰ[11-12]会议中的标准算例, 该模型的风洞实验数据详细, 国内外许多研究者都对该模型进行复杂湍流流动数值模拟以检验所采用的RANS/LES混合方法的可靠性。采用所开发的计算程序对串列双圆柱模型绕流流场进行了数值模拟分析, 检验所建立的基于k-w SST湍流模型的IDDES数值模拟方法及采用的结构网格面搭接技术的可靠性。
2.1 计算模型及网格在计算时2个圆柱之间的距离L选为圆柱直径D的3.7倍, 与风洞实验模型相同。在BANC-Ⅰ会议中定义串列双圆柱计算模型的标准展长为3倍圆柱直径[13]。串列双圆柱计算模型展长也取为圆柱直径的3倍。
对双圆柱流场进行数值模拟时, 流场计算域选取为:沿流向方向上总长度为60D, 上游边界到上游圆柱距离为20D, 下游边界距下游圆柱距离为40D; 沿法向方向上长度为20D; 展向长度为3D。图 1给出了串列双圆柱流场空间截面网格图, 在流动关注区域内, 网格单元在三个方向上的最大尺度约为0.02D。在近壁面区域, 物面法向方向上第一层网格到物面的距离为0.000 01D, 整个流场区域内计算网格节点约为2 000万。将网格搭接面设置在远离双圆柱的区域, 这样既能减少一部分的计算量, 同时又能避免搭接面穿过流动变化剧烈区域从而导致计算结果的不准确。
|
| 图 1 双圆柱空间截面网格图 |
计算条件参照实验状态取为:来流速度为44m/s, 基于圆柱直径的雷诺数为1.66×105。计算中展向2个端面为周期性边界条件, 圆柱表面为无滑移物面边界条件, 流场周围其余边界为远场边界。空间无黏项离散采用五阶Roe-WENO格式, 黏性项采用中心格式进行离散, 时间离散采用二阶隐式格式, 无量纲化时间步长取为0.01, 每个时间步长推进过程中取25步子迭代。
图 2给出了计算所得圆柱表面压力脉动系数分布与实验结果的对比, 计算所得上游圆柱表面压力脉动系数分布形态与峰值与实验结果吻合较好, 证明计算能够准确地捕捉到圆柱表面的流动分离现象。
|
| 图 2 圆柱表面压力脉动系数分布图 |
图 3给出了计算得的流场中空间截面展向涡量分布同实验结果的对比。实验得到的展向涡量分布为PIV测量得到的结果, 图中展示出了两圆柱之间既有大尺度的涡结构也有小尺度的涡结构。由上游圆柱所产生的剪切层十分薄, 同时剪切层在向下游发展过程中会上下飘动也增大了数值模拟的难度。从展向涡量分布图中可以看出, 计算结果揭示了上游圆柱尾迹区内由剪切层的发展过程并捕捉到了由剪切层所拖出的大尺度涡脱落现象。
|
| 图 3 瞬时流场展向涡量分布图 |
分别采用基于结构网格动态面技术的RANS方法和IDDES方法针对Caradonna-Tung旋翼[14]的桨尖涡进行了数值模拟对比分析研究, 该旋翼桨盘半径R=1.143 m, 展弦比为6, 在计算时对模型进行了适当简化, 忽略了桨毂的影响。为了能够捕捉到悬停状态旋翼流场的流动细节, 对旋翼下方尾迹区内的网格进行了精心设计, 在网格划分时确保关注区域内网格各向同性且尺度保持在桨尖弦长的5%左右。图 4给出了流场空间截面网格图(实际计算网格比图中密), 图中可以看到流场中静止区域及旋转区域的搭接情况, 整个流场网格总数量约为9 000万。
|
| 图 4 空间截面网格图 |
计算状态为:桨尖马赫数Ma=0.44, 桨叶总距角θ=8°, 基于桨尖弦长的雷诺数为Re=1.92×106。采用双时间方法推进, 每个旋转周期计算1 440步, 子迭代步数为25步, 共计算20个周期。为比较RANS方法与IDDES方法的计算结果, 2种方法采用的网格及空间离散格式(五阶Roe-WENO)均一致。
论文对2种计算方法所需的计算时间进行了统计, 在相同的计算机配置下(4台6核12线程的I7计算机并行), 采用RANS方法计算一个旋转周期需要约20 h, 而采用IDDES方法则需要约25 h。IDDES方法计算所需时间比RANS方法多25%。
图 5给出了旋翼悬停流场空间瞬态Q等值面图对比, 采用涡量值大小进行渲染。可以看出在相同的网格分布下, IDDES方法计算得到的旋翼悬停流场细节更为丰富一些, 并且所捕捉到的桨尖涡及尾涡面范围更广。图 6给出了某单片桨叶的实验图, 从图中可以看到旋翼桨尖涡附近存在许多小尺度的涡结构, 称为蠕虫涡结构, RANS方法计算结果中不能观察到这一细小涡结构。图 7给出了IDDES计算得到的旋翼Q等值面的剖视图, 从中可以清楚看到旋翼流场的结构特征, 与图 6给出的实验现象更为接近。这主要是因为IDDES方法在旋翼尾迹关注区域内的亚格子湍流黏性更小, 有利于小尺度涡结构的解析, 计算捕捉到了桨尖涡附近的细长型蠕虫涡结构, 这种小涡结构在旋翼桨尖涡向下发展到一定程度后开始出现并绕着桨尖涡分布在流场中。同时, IDDES方法计算得到的旋翼桨叶拖出的尾涡面结构也比RANS方法要丰富一些。细小的涡结构组成了尾涡面并随着旋翼的转动方向有序的排列。
|
| 图 5 Q等值面对比图 |
|
| 图 7 Q等值面剖视图(IDDES计算结果) |
图 8给出了RANS方法和IDDES方法计算得到的流场空间截面涡量分布图, 该图反映出桨尖涡向下发展的情况, 并且桨叶拖出的尾涡面较桨尖涡向下发展速度要快。
|
| 图 8 空间截面涡量对比图 |
从IDDES方法的计算结果可以观察到, 桨尖涡向下发展到一定程度后, 上下相邻的桨尖涡干扰作用变强, 涡核位置有了一定变化, 附近的蠕虫涡结构也相互混合到一起, 类似的现象在Chaderjian的研究中也能观察到, 称为“涡对”(vortex pairing)[3]现象, 且与参考文献[15]给出的某2片桨叶旋翼的实验烟流图的现象一致, 而采用非定常RANS方法得到的计算结果中很难捕捉到该现象。进一步表明IDDES方法对于旋翼桨尖涡的分辨率较高。
图 9给出了RANS方法及IDDES方法计算得到的桨尖涡不同尾迹角处涡核横向截面速度型分布对比图。由于IDDES方法计算得到的桨尖涡结构形态丰富并且涡量强度大于RANS方法, 因此从图中可以看到在同一尾迹角处, IDDES方法计算得到的Z方向速度峰值要大于RANS方法结果。而对比Z方向速度最大值与最小值位置之间的距离发现, RANS方法计算得到的桨尖涡涡核直径要大于IDDES方法。
|
| 图 9 不同尾迹角处涡核附近速度型对比图 |
图 10则给出了IDDES方法及RANS方法计算得到的旋翼桨尖涡涡核直径随尾迹角的变化图, IDDES方法计算得到的旋翼桨尖涡涡核直径远小于RANS方法, 与参考文献[3]中所说的桨尖涡涡核直径一般为旋翼桨叶桨尖弦长的10%左右更为接近(Caradonna-Tung旋翼的桨尖弦长为190 mm)。图 11则给出了RANS方法和IDDES方法计算得到的桨尖涡涡核位置与实验值的对比, 对涡核位置进行统计时进行了时均处理。从图中可以看到, RANS方法与IDDES方法计算得到的桨尖涡涡核位置均与实验值吻合得较好, 而在尾迹角增大到200°以后, IDDES方法的结果更好一些, 与实验更为接近。
|
| 图 10 桨尖涡涡核直径随尾迹角变化对比图 |
|
| 图 11 桨尖涡涡核位置随尾迹角变化对比图 |
为进一步验证本文所建立的旋翼桨尖涡IDDES方法的可靠性, 对某单片桨叶的旋翼悬停流场进行了数值模拟, 选取该桨叶进行计算分析是因为该桨叶有比较详细的涡核速度型分布实验结果, 该桨叶为单片直桨叶, 采用NACA2415翼型, 其具体参数及相关实验数据详见参考文献[16-17]。
采用文中开发的IDDES程序进行计算, 每个旋转周期计算1 440步, 内迭代步数取25步, 一共计算20个旋转周期。图 12给出了计算得到的该单片桨叶的Q等值面图, 采用涡量值大小进行渲染, 从图中可以看到计算得到的流场细节与上节的Caradonna-Tung旋翼同样丰富, 捕捉到了许多小尺度的蠕虫涡结构。图 13给出了不同尾迹角处涡核速度型的计算结果与实验值的对比图, 横坐标为径向位置, 采用旋翼桨叶弦长单位化, 纵坐标为径向不同位置处的旋转速度, 采用桨叶桨尖速度单位化。可以看到计算值与实验值吻合良好, 表明采用文中所建立的IDDES方法可以很好地计算出旋翼桨尖涡附近的速度型分布, 适用于旋翼桨尖涡流场的精细化数值模拟研究。
|
| 图 12 某单片桨叶Q等值面图 |
|
| 图 13 不同尾迹角处涡核附近速度型对比图 |
本文基于结构网格动态面搭接技术、五阶Roe-WENO格式，开发了一套用于悬停状态旋翼桨尖涡高精度数值模拟的IDDES方法计算程序，并采用该程序针对Caradonna-Tung旋翼和某单片桨叶旋翼的悬停流场进行了数值模拟研究，从桨尖涡的速度型分布、涡核直径、涡核位置等参数对悬停状态旋翼桨尖涡的流动进行了比较系统的研究，得到了以下结论：
1) IDDES方法与RANS方法计算得到的旋翼桨叶气动力基本一致，但IDDES方法计算得到的旋翼桨尖涡速度型分布、涡核直径、涡核位置与实际情况更为接近；
2) IDDES方法捕捉到了旋翼桨尖涡附近的细小蠕虫涡结构以及桨尖涡的"涡对"现象，与国外相关实验及研究结论一致，有利于深入研究旋翼绕流机理；
3) 与RANS方法相比，在相同的网格分布及空间离散格式下，IDDES方法的优势在于对旋翼空间流场的精细化数值模拟。
| [1] | GARY R B. An Aerodynamic Analysis of a Single Bladed Rotor in Hovering and Low Speed forward Flight as Determined from Smoke Studies of the Vorticity Distribution in the Wake[R]. Princeton University, Report, 1956: 356 |
| [2] | NATHAN Hariharan, ALAN Egolf, LAKSHMI Sankar. Simulation of Rotor in Hover: Current State and Challenges[R]. AIAA-2014-0041 |
| [3] | CHADERJIAN N M. Advances in Rotor Performance and Turbulent Wake Simulation Using Des and Adaptive Mesh Refinement[C]//7th International Conference on Computational Fluid Dynamics, Big Island, Hawaii, 2012 |
| [4] | CHADERJIAN N M, BUNING P. High Resolution Navier-Stokes Simulation of Rotor Wakes[C]//AHS 67th Annual Forum, Virginia Beach, 2011: 3-5 |
| [5] | CHADERJIAN N M, AHMAD J U. Detached Eddy Simulation of the UH-60 Rotor Wake Using Adaptive Mesh Refinement[C]//Proceedings of the American Helicopter Society 68th Annual Forum, 2012: 1-3 |
| [6] | BORGES R, CARMONA M, COSTA B, et al. An Improved Weighted Essentially Non-Oscillatory Scheme for Hyperbolic Conservation Laws[J]. Journal of Computational Physics, 2008, 337: 3191-3211. |
| [7] | RAI M M. A Relaxation Approach to Patched-Grid Calculations with Euler Equations[R]. AIAA-1985-0295 |
| [8] | TRAVIN A, SHUR M, STRELETS M. Physical and Numerical Upgrades in the Detached-Eddy Simulation of Complex Turbulent Flows[C]//12th Euromech Colloquium on LES and Complex Transitional and Turbulent Flow, Munich, Germany, 2000 |
| [9] | XIAO Z X, FU S. Studies of the Unsteady Supersonic Base Flows around Three after Bodies[J]. Acta Mechanica Sinica, 2009, 25(4): 471-479. DOI:10.1007/s10409-009-0248-4 |
| [10] | ANDREY K, TRAVIN A, MIKHAIL L S, et al. Improvement of Delayed Detached-Eddy Simulation for LES with Wall Modeling[C]//European Conference on Computational Fluid Dynamics, TU Delft, the Netherlands, 2006 |
| [11] | DAVID P. Summary of the Tandem Cylinder Solutions from the Benchmark Problems for Airframe Noise Computations-Ι Workshop[C]//49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition 2001, Orlando, Florida |
| [12] | JENKINS L N, KHORRAMI M R, CHOUDHARI M M, et al. Characterization of Unsteady Flow Structures around Tandem Cylinders for Component Interaction Studies in Airframe Noise[R]. AIAA-2005-2812 |
| [13] | SPALART P R, MEJIA R. Analysis of Experimental and Numerical Strudies of the Rudimentary Landing Gear to Validate Noise Predictions[R]. AIAA-2011-0355 |
| [14] | CARADONNA F X, TUNG C. Experimental and Analytical Studies of a Model Helicopter Rotor in Hover[J]. Vertica, 1981, 5(2): 149-161. |
| [15] | KOMERATH N M, SMITH M J. Rotorcraft Wake Modeling: Past, Present, and Future[C]//Proceedings of the European Rotorcraft Forum, Hamburg, Germany, 2009 |
| [16] | HAN Y O, LEISHMAN J G. Investigation of Helicopter Rotor-Blade-Tip-Vortex Allevation Using a Slotted Tip[J]. AIAA Journal, 2004, 42(3): 524-535. DOI:10.2514/1.3254 |
| [17] | HAN Y O, LEISHMAN J G. Hovering Performance of a Rotor with Slotted Blade Tips[C]//Proceedings of the 60th AHS Annual Forum, Baltimore, 2004 |
2. School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China
