12高精度旋翼非定常涡流场数值模拟方法史勇杰8

1、第二十八届（2012）全国直升机年会论文高精度旋翼非定常涡流场数值模拟方法史勇杰 魏 鹏 徐国华 招启军（南京航空航天大学直升机旋翼动力学重点实验室，南京，210016）摘 要：基于粘性涡方法和CFD方法，发展了一套适合于复杂旋翼涡流场的耦合拉格朗日-欧拉计算方法。方法中针对高雷诺数旋翼流场中桨尖涡的紧凑结构特点，引入粘性涡方法建立了高分辨率的尾迹求解模型，以此模拟尾迹的畸变运动及其影响；为捕捉激波等细节流场特征，在需要考虑细节流场的桨叶近体区域采用了CFD方程对其进行求解；而两者之间的信息交换分别采用集中涡源法和边界修正法完成。应用所建立的计算方法，对Caradonna-Tung模型旋翼悬停

2、状态下的旋翼流场进行了模拟，通过分析旋翼尾迹特征及耦合边界处流场特征，表明所建方法充分发挥了CFD方法和粘性涡方法各自的优势，在旋翼流场计算方面的具有独特的优势。最后通过对比桨叶表面压强系数分布，进一步验证了该方法在旋翼流场分析中的有效性。关键词：旋翼流场；旋翼尾迹；Euler 方程；离散涡方法；混合方法；直升机0 引言桨叶近体的复杂非线性流动和桨叶产生的强烈桨尖涡是旋翼非定常流场的重要特征，两者相互耦合、相互影响，对直升机的气动载荷、气动噪声以及动力学特性等有重要影响，关于它们的研究一直是直升机技术研究领域的一个重点。然而，同时准确预测桨叶附近的非线性、跨音速、可压缩流动和旋翼尾迹的畸变运动

3、具有很大难度，这给旋翼气动特性的深入研究提出了巨大挑战。迄今为止，对于旋翼非定常涡流场的研究，主要采用两类数值计算方法，即基于拉格朗日法的涡动力学方法1和基于欧拉法的计算流体力学（CFD）方法2。其中涡动力学方法开展较早，传统尾迹线涡法在经历了固定尾迹、预定尾迹以及自由尾迹的发展历程之后日臻成熟，目前基于桨尖涡卷起模型的自由尾迹方法已被广泛应用于各种理论分析及工程计算中。而近几年出现的高分辨率粘性涡方法成为继自由尾迹方法后的又一个重要发展趋势，Opoku，Anusonti，Zhao3等人先后开展了基于拉格朗日描述的旋翼粘性离散涡方法方面的研究，通过对不同旋翼流场的计算表明，所发展的方法可以更准

4、确的预测尾迹涡的畸变运动及扩散效应。然而，这些方法中均采用了线性气动力模型及不可压假设，因此对于桨叶近体非定常流动的模拟存在困难，无法满足现代复杂外形桨叶分析需求。上世纪80年代以来， CFD方法开始应用于旋翼气动问题的研究中，该方法在欧拉描述下通过数值方法求解Navier-Stokes方程，可以精确地捕捉桨叶表面气动力的细致变化，以及前行侧桨叶的强压缩和后行侧的反流、动态失速等复杂旋翼扰流特性。然而该方法在采用有限体积或有限差分法进行控制方程求解时，离散格式中固存的数值耗散会使得空间流场的涡出现非物理耗散，从而导致方法在无法准确捕捉旋翼尾迹。尽管很多学者已尝试采用高阶格式4、自适应网格5等技

5、术消除这一缺陷并取得了较好的计算结果，但在现有计算资源条件下，尾迹的捕捉能力仍显不足。考虑到动力学方法与CFD方法各自的优势，许多学者将两者相结合发展了混合计算方法1。早期的混合方法多基于线涡法和CFD方法，Berezin和Sankar等7基于CFD/全位势/自由尾迹建立了一套混合方法，而Bader等采用场速度方法直接将CFD与自由尾迹耦合在了一起，计算结果表明，这些方法均可有效的增加尾迹捕捉能力，从而提高计算精度，然而受自由尾迹方法本身缺陷限制，它们的通用性并不好。近两年Stock6，Zhao7等初步将高分辨率涡方法与CFD方法结合，发展了具有高精度的耦合方法，从已发表的文献看，耦合方法在尾

6、迹捕捉和旋翼气动性能预测等方面都要好于传统的CFD方法，而且计算量也有显著降低。目前国内尽管已有不少学者开展了自由尾迹/CFD混合方法研究8, 9，但在粘性涡/CFD耦合方法方面尚未有相关研究。为此，本文在先前研究工作的基础上，进一步将粘性涡方法与CFD方法相结合，建立了适合于复杂旋翼涡流场计算的耦合拉格朗日-欧拉计算方法。然后，采用Caradonna-Tung试验旋翼为算例对所建方法的有效性进行了分析及验证，结果表明该方法充分结合了两种方法各自的优势，无论是在旋翼尾迹捕捉方面还是在桨叶近体流场细节模拟方面都是非常有效的。1 计算方法图1给出了本文耦合拉格朗日-欧拉法的示意图，该机算方法将旋翼

7、气动问题分为尾迹预测和桨叶近体流场捕捉两部分：在离桨叶较远区域，尾迹的影响占据着主导地位，故采用基于拉格朗日描述的粘性涡方法进行尾迹预测；而围绕桨叶的旋翼近体流场则应用基于欧拉描述的CFD方法进行求解，以捕捉该流场区域的各种非线性特征以及近尾迹影响。图 1 耦合拉格朗日-欧拉法示意图1.1 旋翼尾迹求解模型直升机飞行时，旋翼桨叶拖出的尾涡发生卷起现象并形成桨尖涡，这些结构紧凑且强度很高的桨尖涡堆积在桨盘下（后）方，成为旋翼尾迹的主要组成部分。拉格朗日粘性涡方法将该尾涡系离散为一系列无序涡元，并在拉格朗日描述体系进行数值求解来预测尾迹的卷起、畸变及耗散等过程。不考虑桨叶附近流场，直升机流场问题可

8、认为不可压问题，离散涡元随时间的变化由运动学方程和涡量动力学方程进行求解 式中，下标为涡元索引号，表示位置，为环量值，为当地速度,为涡诱导速度，为时间，为运动粘性系数。式(1)中，涡量动力学方程的右边第一项表示流场的速度梯度引起的涡元的伸缩与弯曲，其中的自诱导速度及速度梯度项由Biot-Savart定理计算，并引入快速多极子加速算法(FMM, Fast Multipole Method)10对该计算过程进行加速。而第二项表示的是涡量的粘性扩散效应，本方法采用粒子强度交换法(PSE, Particle Strength Exchange)11对其进行求解。受桨尖涡诱导作用，旋翼流场存在着的强烈剪

9、切流导致长时间数值计算后离散涡元会出现严重的疏密失调，导致涡元间无法满足保证计算精度的重叠性要求，从而致使计算精度下降。为此，本文采用文献12的涡元重布方法来调整涡元的疏密度，该方法根据当地速度应变张量判断局部流场的拉伸状况，当涡元所在位置的局部拉伸达到特定值时，将涡元分裂成两个相同强度的涡元填补由于拉伸导致的“空白”。然而，随着新涡元从桨叶脱落以及旧涡元的分裂，计算过程中涡元数目会迅速增长，导致计算效率降低；为合理控制因涡元数目增长导致的效率问题，一方面可以将远离旋翼的涡元删除，另一方面则需要对距离很近且方向较为一致（涡矢量夹角小于特定值）的涡元进行合并处理。1.2 桨叶流场计算模型为捕捉桨

10、叶近体空间的非定常流动，采用以绝对物理量为变量的守恒形式的雷诺平均N-S方程对其进行数值求解。在惯性系(x,y,z)下，方程表示为 式中，守恒变量，为惯性系下三个方向上的速度，分别表示压强、密度和能量；为控制体单元体积，为面积，分别为对流通量喝粘性通量，为粘性项逻辑开关。本文旨在研究所建方法的耦合有效性，着重于两计算域在交互过程中的信息传递，考虑到粘性项对此没有直接影响，在实际计算过程中一方面出于方便对数值耗散进行研究的缘故，另一方面出于效率方面的原因，并未打开粘性项逻辑开关。计算过程中所有的物理量均以悬停桨尖速度、密度和桨叶弦长作为参考值而无量纲化。控制方程的空间离散采用二阶中心方法，而时间

11、推进上采用双时间方法。2 计算域间的信息交换对于本文所建方法，拉格朗日计算域与欧拉计算域之间的相互信息交换至关重要，合理的信息交换策略无论对计算精度还是对计算效率都有帮助，因此需要对其进行仔细处理。2.1 欧拉计算域至拉格朗日计算域CFD计算域至粘性涡计算域的信息传递，可采用两种方式：集中涡源法和分布涡方法。后者直接将CFD计算域计算的涡量分布传递给拉格朗日计算域，物理意义更为明确，但却容易将CFD计算域的数值耗散引入到拉格朗日计算域，同时也存在着计算工作量大的问题。为此，本文采用了集中涡源法，即先由CFD方法计算得到桨叶展向载荷分布，再根据库塔如科夫斯基（Kutta-Joukowski）环量

12、/升力定理来获得桨叶展向附着涡环量分布 其中是桨叶剖面升力系数，为附着涡环量值。桨叶附着涡随空间及时间的变化是桨叶脱落尾涡的主要原因，其中附着涡随空间变化脱落的尾涡为尾随涡，随时间变化脱落的桨尖涡为脱体涡。因此在求解得到桨叶附着涡后，拉格朗日计算域新生涡元的强度即可由下面公式确定 式中，表示新生涡环量值，为桨叶剖面相对气流速度，包括来流速度、旋翼旋转速度以及桨叶挥舞运动速度。2.2 拉格朗日计算域至欧拉计算域粘性涡计算域至CFD计算域的信息传递，常用的有几种方式，第一种方法是采用涡诱导速度直接修正桨叶当地迎角；第二种是采用广义网格速度法将影响施加在全体网格格点中，即将涡线产生的诱导速度等效成网

13、格畸变；第三种方式直接通过边界条件施加于CFD计算域外边界上。其中，第一种方式尽管效率很高，然而在迎角修正过程中会引入误差，精度并不高；第二种方法尾迹贡献施加在所有网格点上，桨尖涡对流场的强烈影响可以直接体现出来，但这也导致其效率很低；第三种方法较前两种方法更符合物理意义，在控制好网格质量的前提下可以得出满意结果（参考文献），而且仅需计算外边界处的诱导速度，效率较高，因此本文采用了第三种方式。计算过程中，所需边界条件除速度外，还需要密度和能量，本文假定CFD计算域大小恰当，桨叶网格外边界以外压缩性已经可以忽略，因此密度等于远场处密度，而能量边界条件为 3 方法验证及分析为验证本文方法对悬停流场

14、的计算能力，采用了有实验结果可供对比的“Caradonna & Tung”试验旋翼做为算例进行计算，为避免气弹作用对桨叶气动特性的影响，旋翼桨叶具有很高的设计挥舞刚度，因此被广泛用于旋翼流场计算的对比研究。该桨叶展弦比为6，弦长0.1905，翼型为NACA 0012，无扭转和尖削。试验进行了不同总距、不同桨尖马赫数的多个算例，本节中若无特别说明，选取总距8度、桨尖马赫数0.439的状态作为默认算例。3.1 黏性涡方法与CFD方法算例分析作为各自方法的验证，以及为本文耦合方法提供算例对比，首先分别采用黏性涡方法和CFD方法独立的对默认算例进行了计算。对于拉格朗日涡方法，影响其计算精度及计

15、算时长的主要参数有光滑参数（为重叠因子，为最小涡元尺寸），桨叶分段数，方位角步进步长，计算域控制参数（距桨毂中心大于该参数的涡元将被删除），以及涡元合并距离系数和涡元合并角度阀值。本文各参数取值为：、，计算结束时涡元数约为3万多。而对于CFD计算方法，网格数量是影响精度和效率的主要因素。在桨叶网格方面，为了方便对比，包括下一节耦合方法计算在内所有计算均采用了同一套网格，即213×37×89的CO型贴体网格，网格外边界距离桨叶物面约1倍弦长。另外，本节的全CFD求解模型采用如图2所示的嵌套网格系统，其背景网格采用关键部位加密处理的笛卡尔网格，为了对比研究背景网格对桨尖涡数值耗

16、散的影响，分别采用137×107×137和118×85×118两套背景网格进行了计算。图 2 用于全CFD计算的嵌套网格系统图3首先给出了黏性涡方法计算结束时的涡元分布图，其中涡元的体积表征环量值的大小，可以看出，受下方尾迹的诱导作用，旋翼桨尖涡在离开桨叶后迅速收缩。随着寿命角的增大，收缩运动逐渐减缓，而桨尖涡部分区域出现非规则畸变，这导致大约在两圈以后相互靠近的桨尖涡间开始出现严重干扰现象，桨尖涡间的这种干扰现象导致旋翼涡系发生破坏并最终形成紊流状态，文献13对此现象也进行了分析，且给出了本文一致的结论。(a) 黏性涡方法(b) CFD方法图3 流场空

17、间涡量等值面图旋翼桨尖涡的收缩畸变是悬停旋翼流场的重要特征，准确的捕捉该畸变对于旋翼气动分析具有重要意义。图3给出了本文耦合方法计算得到的桨尖涡径向位置及轴向位置随寿命角的变化趋势，由图可见，桨尖涡轴向位移的计算值在小寿命角时变化缓慢，当寿命角达到180度时，桨尖涡受另外一片桨叶桨尖涡的诱导影响，轴向位移明显加快，这与实际旋翼尾迹运动特性相符合，也与试验值吻合很好，同时可以看出径向收缩量也与试验值得到了很好吻合，这为进一步开展拉格朗日-欧拉法的研究提供了必要的前提保证。图4 计算与试验的桨尖涡径向位置及轴向位置对比图3.2 耦合拉格朗日-欧拉法的验证与分析图5为该算例状态下的涡量图，为显示直观

18、，CFD计算域的涡量用桨叶网格外表面的涡量云图表示，拉格朗日计算域的涡量采用涡量等值面显示。从图中可以看出，两个计算域的的涡量集中区域完全吻合，当拉格朗日计算域的桨尖涡量穿入CFD计算域时，涡量信息通过CFD计算域外边界的速度边界条件成功的完成了传递，涡量的这种信息交换使得本文方法可非常有效的改善CFD计算域的数值耗散问题，同时也有助于增强桨涡干扰特性的捕捉。图5 耦合拉格朗日-欧拉法计算得到的涡量图图6则给出了流线分布图，其中有压力云图的部分为CFD计算域的计算结果，白色背景部分为拉格朗日涡方法的计算结果，图中可以看出，除边界处有强烈的涡作用的位置外，两区域的流线平滑连接，具有很好的连续性，

19、表明拉格朗日计算域的流场信息通过边界处信息交换有效的传入到了CFD计算域。（a）r=0.4R(b) r=0.6R(c)r= 0.8R图6 桨叶剖面流线计算值最后给出桨叶不同剖面处表面压强系数分布曲线（图7），可以看出，计算得到的压强系数分布曲线与试验值在各个剖面处都吻合很好，表明本文方法在低桨尖马赫数下可以准确的计算桨叶细节气动特性。而对于高桨尖马赫数状态，桨叶近体流场的压缩效应更大，为验证本文方法的不可压边界条件假设对高桨尖马赫数算例状态的影响，又计算了桨尖马赫数为0.877的算例状态，图8给出了该状态下的桨叶表面的压强系数分布，从图中看出，本文建立的耦合方法可以捕捉到激波导致的压强突变，并

20、与试验值基本吻合，仅激波位置较试验值有所偏后，这可能是由于本文CFD计算域所采用了欧拉求解器，而在无激波位置，计算值与试验值吻合的则更好。(a) r/R=0.8(b) r/R=0.89(c) r/R=0.96图7 桨叶表面压强系数分布（Ma=0.439）(a) r/R=0.8(b) r/R=0.89(c) r/R=0.96图8 桨叶表面压强系数分布（Ma=0.877）4 结论本文基于粘性涡方法和CFD方法，发展了一套适合于复杂旋翼涡流场的耦合拉格朗日-欧拉计算方法。通过Caradonna-Tung模型旋翼悬停状态下的旋翼流场进行了模拟，通过分析旋翼尾迹特征及耦合边界处流场特征，表明所建方法充分

21、发挥了CFD方法和粘性涡方法各自的优势，在旋翼流场计算方面的具有独特的优势。参 考 文 献1Komerath N M, Smith M J, Tung C. A Review of Rotor Wake Physics and Modeling J. Journal of the American Helicopter Society, 2011, 56(2): 220061-2200619.2Strawn R C, Caradonna F X, Duque E P N. 30 years rotorcraft computational fluid dynamics research and

22、 development J. Journal of the American Helicopter Society, 2006, 51(1): 5-21.3Zhao J, He C. A Viscous Vortex Particle Model for Rotor Wake and Interference Analysis J. Journal of the American Helicopter Society, 2010, 55(1): 1-14.4Hariharan N, Sankar L. First-principles based high order methodologi

23、es for rotorcraft flowfield studies; proceedings of the Proceedings of 55th Annual Forum of AHS, F, 1999 C. American Helicopter Society, 1999.5Dietz M, Kessler M, Krämer E. Tip Vortex Conservation on a Helicopter Main Rotor Using Vortex-Adapted Chimera Grids J. AIAA Journal, 2007, 45(8): 2062-2

24、074.6Stock M J, Gharakhani A, Stone C P. Modeling Rotor Wakes with a Hybrid OVERFLOW-Vortex Method on a GPU Cluster; proceedings of the 28th AIAA Applied Aerodynamics Conference, Chicago, Illinois, F, 2010 C.7Zhao J, He C. A Hybrid Solver with Combined CFD and Viscous Vortex Particle Method; proceed

25、ings of the American Helicopter Society 67th Annual Forum, Virginia Beach, VA, F, 2011 C. American Helicopter Society.8Cao Y-H, Yu Z-Q, Su Y, et al. Combined Free Wake/CFD Methodology for Predicting Transonic Rotor Flow in Hover J. Chinese Journal of Aeronautics, 2002, 15(2): 65-71.9Shi Y, Zhao Q, F

26、an F, et al. A New Single-blade Based Hybrid CFD Method for Hovering and Forward-flight Rotor Computation J. Chinese Journal of Aeronautics, 2011, 24(2): 127-135.10Greengard L, Rokhlin V, 1997.11Mas-Gallic S. Contribution à l'analyse numérique des méthodes particulaires D, 1987.12

27、Mansfield J R, Knio O M, Meneveau C. Dynamic LES of Colliding Vortex Rings Using a 3D Vortex Method J. Journal of Computational Physics, 1999, 152(1): 305-345.13Leishman J G. Principles of helicopter aerodynamics M. 2nd ed. Cambridge ; New York: Cambridge University Press, 2006.High Precision Numeri

28、cal Method for Simulating Rotor Unsteady Votex Flow FieldSHI Yong-jie WEI Peng XU Guo-hua ZHAO Qi-jun(National Key Laboratory of Rotorcraft Aeromechanics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016)Abstract: A new coupled Eulerian-Lagrangian numerical method is developed for complex rotor vortex flowfiled simulation based on the viscous vortex method and the CFD method. Considering the compact structure of the tip vortex in the high Reynolds rotor fl