非定常旋涡-翼型相互作用问题数值模拟研究

赵文梅，赵 军，胡 偶

(中航工业直升机设计研究所，江西 景德镇 333001)



非定常旋涡-翼型相互作用问题数值模拟研究

赵文梅，赵 军，胡 偶

(中航工业直升机设计研究所，江西 景德镇 333001)

针对直升机旋涡-桨叶相互作用问题，以其二维简化模型旋涡-翼型相互作用问题为研究对象，发展了基于雷诺平均Navier-Stokes方程和自适应混合笛卡尔网格的流场数值模拟方法，通过给定合适的自适应判据，动态网格自适应技术能够准确捕捉旋涡位置和保持旋涡强度。通过两类旋涡-翼型相互作用问题的数值模拟研究，揭示了旋涡与翼型相互作用过程中产生的复杂流动现象和非定常气动力变化。数值模拟结果与相关文献中的试验数据和数值计算结果吻合，且对三维桨叶的气动载荷研究和振动分析具有指导意义。

混合笛卡尔网格；自适应；非定常；气动力；旋涡-翼型相互作用

0 引言

直升机旋翼在转动过程中桨叶的桨尖处会产生非常强烈的尾迹涡结构，且随着尾迹涡的发展会与后续桨叶发生相互作用，即直升机空气动力学中一个重要的研究课题：旋涡-桨叶相互作用问题(Blade-Vortex Interaction, BVI)。BVI问题对直升机的旋翼振动、气动噪声等都有显著的影响。目前，针对BVI问题，为简化问题模型和发展相关数值方法，一般将其简化为二维的旋涡-翼型相互作用问题(Airfoil-Vortex Interaction, AVI)。通过对AVI问题的研究，发展相应的数值模拟方法，揭示旋涡的发展过程以及旋涡与翼型相互作用过程中引起的复杂流动现象和气动力变化。

AVI问题涉及非常复杂的流动现象，例如旋涡的传播、激波的形成、二次涡的产生、旋涡碰撞、以及声波传播等等。诸多学者采用计算流体动力学(CFD)技术开展了AVI问题的数值模拟研究，主要有：Lee和Berchader[1]的高精度数值方法，Falissard和Lerat[2]提出的涡量保持格式(vorticity-preserving schemes)，Morvant[3]的涡量限制格式(vorticity-confinement schemes)，Lei Tang和Bader[4]基于三阶精度的TVD格式并结合网格重分布方法，Woo等人[5-6]的非结构网格自适应技术等。

传统的基于静态网格的CFD流场数值模拟技术，在处理AVI问题时，由于网格数目相对固定，因而不能准确地捕捉由于旋涡的运动引起的非定常流动现象，包括旋涡强度的保持、位置的捕捉以及碰撞过程中产生的复杂流动现象等。为此，本文基于雷诺平均Navier-Stokes方程，针对旋涡运动问题，发展了一套自适应混合笛卡尔网格(Adaptive Hybrid Cartesian Grid, AHCG)流场数值计算方法。

本文首先给出了AVI问题的几何模型；其次，建立AVI问题的数值模拟方法；最后，将发展的数值方法用于AVI问题的数值模拟研究。重点分析了旋涡的发展历程，旋涡与翼型相互作用过程中复杂流场现象，以及作用在翼型上的气动力变化。

1 AVI问题几何与数学模型

AVI问题的试验结果来源于Lee和Bershader[1]等人的研究，通过激波与NACA0018翼型衍射作用形成旋涡结构，然后与下游的NACA0012翼型碰撞。AVI问题的几何模型如图1所示。

图1 AVI问题几何模型

其中，(ΔX,ΔY)为旋涡的放置位置，Γ为旋涡的强度，V∞为来流速度，c为翼型弦长。

根据试验测量结果，Scully旋涡模型被用于定义数值模拟中的旋涡结构。在极坐标系(r,θ)下，Scully旋涡模型的切向速度Xθ的表达式为：

其中，r为网格点到涡核的距离，rc表示涡核半径。另外，还需要给出压强和密度关系，由于试验中旋涡是由激波衍射产生的，因此根据等焓假设和径向动量方程[3]，可得压力关于半径的解析表达式：

其中

对于(2)式，当r→∞时，p(r)→p∞。密度则可由等焓假设[3]计算得到。

2 AVI问题数值模拟方法

精确模拟AVI问题，需要重点关注的两大问题是：旋涡位置的捕捉和强度的保持，以及旋涡在与翼型相互作用过程中产生的复杂流动现象的准确模拟。因此，本文基于雷诺平均Navier-Stokes 方程(RANS)，发展了一套能够根据旋涡特征跟踪旋涡结构和动态自适应分布计算网格的自适应混合笛卡尔网格方法(AHCG)。该方法采用混合笛卡尔网格系统创建初始计算网格，在流场计算过程中，采用基于流场特征的网格自适应技术动态分布网格密度，从而达到流场精确数值模拟的目的。

2.1 混合笛卡尔网格方法

如图2所示，混合笛卡尔网格是由物面附近的贴体结构网格和填充流场其余区域的笛卡尔网格构成，两套网格之间信息传递是通过两套网格之间的网格交接面完成的。

与传统的混合网格方法类似，混合笛卡尔网格在笛卡尔网格与贴体网格之间需要进行数据传递。网格交接面的信息的传递直接影响计算的守恒性、数值模拟的精度以及计算效率。如图3所示，图中的实线表示笛卡尔网格的交接面，虚线表示贴体网格的交接面，当需要计算通过笛卡尔网格交接面上的网格边IJ的数值通量时，需要知道网格边IJ的左右状态的信息，而其左侧信息可以由笛卡尔网格单元I提供；右侧信息则从贴体网格中寻找一个“贡献单元J(donor cell)”提供。“贡献单元”满足需如下的几何条件：

|dIR|=min(|dIj|),且dIj与nI之间的夹角小于60°

(4)

其中，dIj表示网格单元j的中心相对于网格边IJ中心的位移向量，nI为网格边IJ的单位法向量。

图2 混合笛卡尔网格示意图

图3 交接面“贡献单元”的确定

2.2 基于流场特征的网格自适应技术

在流场数值模拟过程中可以根据流场特征，对流场计算网格进行动态的加密或粗化，达到流场的精确模拟，即基于流场特征的网格自适应技术。本文以速度散度和旋度为解自适应判据，其表达式如下：

然后就可以采用如下的判断标准，确定哪些网格单元需要加密，哪些网格单元需要粗化：

1)如果一个网格单元满足τcI>σc或者τdI>σd，那么这个网格就需要加密；

2)如果一个网格单元满足τcI<0.1·σc并且τdI<0.1·σd，那么这个网格就需要粗化。

采用如上的判据就可以实现计算网格随着流场特征的变化而动态的自适应，从而达到更为准确地捕捉流场信息的目的。

2.3 控制方程与数值方法

考虑可压缩雷诺平均Navier-Stokes方程(RANS)的积分守恒形式如下：

其中，W为守恒变量，Fc和Fv分别为对流通量和粘性通量[7]。在任意的有限控制体ΩI上对方程(7)进行空间离散可以得到如下的半离散格式：

上式中，m表示网格单元I和J之间的网格面，NF表示包围控制体ΩI的网格面的总数目，nm和ΔSm为网格面m的单位外法向矢量和面积。

本文采用有限体积法求解方程(7)，对流项使用具有低数值耗散特性的AUSM+格式[8]求解，并通过解的线性重构[7]获得二阶精度，同时采用Venkkatakrishnan限制器[9]抑制解的非物理振荡；粘性项采用二阶中心型格式离散；时间离散采用LU-SGS隐式迭代方法[10]，同时为提高非定常流动问题的时间模拟精度，采用了双时间步方法(dual-time stepping)[7]。湍流模型采用SSTk-ω湍流模型[11]。

3 AHCG方法数值验证

以二维非定常绕圆柱层流流动问题为例，分析采用AHCG方法和传统静态网格方法，计算由于非定常脱落涡引起的气动力变化方面的差距。基于来流速度和圆柱直径的雷诺数为200，来流马赫数为0.3。计算域为[-5.0,15.0]×[-5.0,5.0]，圆心坐标为(0.0,0.0)。

图4.a比较了采用本文AHCG方法(with AMR)和静态网格(without AMR)方法计算得到的阻力系数CD。从结果可以看出，AHCG方法得到的平均阻力系数为1.331，而静态网格方法得到的平均阻力系数为1.186。图4.b则比较了采用网格自适应技术(with AMR)和不采用网格自适应技术(without AMR)两种方法计算得到的升力系数CL，从图中可以看出，升力系数围绕零升力线振荡，采用网格自适应技术计算得到的Strouhal数为0.195，而不采用网格自适应技术计算得到的Strouhal数为0.166。表1将数值计算得到的阻力系数和Strouhal数与试验结果[12]进行了比较，从结果可以看出，动态网格自适应技术能够有效改善气动力参数的数值模拟结果，与试验结果更为吻合。

表1 圆柱绕流计算结果比较

图4 升阻力系数随时间的变化比较

4 AVI问题数值模拟结果与分析

4.1 低速Ma∞=0.5的AVI问题

基于AVI问题的几何与数学模型，考虑来流马赫数为Ma∞=0.5的流动状态，无量纲的涡核半径和旋涡强度分别为rc/c=0.018，Γ/(V∞c)=-0.283(顺时针方向)，V∞为自由来流速度，旋涡的位置参数为ΔX=-5.0,ΔY=0.0。计算网格采用本文第2节介绍的自适应混合笛卡尔网格。

图5显示了旋涡接近翼型以及与翼型发生碰撞的过程中，不同时刻的自适应计算网格和密度云图，并与试验观测结果[1]进行了对比。可以看出，当顺时针旋涡接近翼型前缘时，翼型前缘出现下洗速度，导致驻点上移。随后旋涡与翼型发生碰撞，旋涡沿翼型下翼面移动并在前缘附近诱导引起“袋状”超音速流动区，之后在翼型的下翼面形成一道弧形激波。当旋涡通过激波之后，会产生一个逆时针旋转的二次涡结构。随着这一对旋涡向下游移动，驻点位置开始下移，并伴有压缩波产生并随之向四周扩散。

图6给出了翼型上下表面不同位置(x=0.02,x=0.05,x=0.1)的压力系数随时间的变化过程，并与试验测量结果[1]和相关文献[6]中的数值模拟结果进行了对比。数值结果表明，当旋涡与翼型发生碰撞时，翼型表面的压力会发生阶跃式剧烈变化，且下翼面有压力脉动产生。这种翼型表面压力的剧烈变化反应到三维桨叶上，桨叶表面会产生压力分布的突变，而这种突变对桨叶的振动和强度都会有影响，因此在设计过程中必须加以考虑。

4.2 高速Ma∞=0.8的AVI问题

考虑跨声速流动状态，来流马赫数为0.8，基于翼型弦长的雷诺数为3.6×106，翼型攻角为0°，无量纲的涡核半径和旋涡强度分别为rc/c=0.05，Γ/(V∞c)=-0.2(顺时针方向)，其中c表示翼型的弦长，V∞为自由来流的速度，旋涡的位置参数为ΔX=-5.0,ΔY=-0.26。该算例没有试验数据，因此，将本文数值模拟结果与相关文献[5,6,13,14]的研究结果进行了对比。

图5 旋涡与翼型碰撞过程中不同时刻瞬态网格，以及密度云图与试验结果的比较

图7给出了翼型升力系数与旋涡涡核位置的之间的关系曲线，并与相关文献的研究结果进行了比较。数值结果显示，本文方法能够准确地模拟随着旋涡的运动而引起升力变化的过程。顺时针的旋涡在接近翼型的过程中产生下洗作用，从而导致初始升力值为负；而随着旋涡的运动，当旋涡靠近翼型前缘位置时，升力系数达到最小值，但随后升力系数迅速恢复；在旋涡沿着翼型下表面传播时，升力系数转变为正值。这种升力的变化，对三维桨叶的受力也具有显著的影响，因此，在桨叶设计过程中需要关注桨涡干扰问题。

图6 AVI问题翼型表面不同位置的压力系数随时间的变化

5 结论

本文针对直升机空气动力学中的旋涡-桨叶相互作用问题(BVI)，以其简化模型旋涡-翼型相互作用问题(AVI)为研究对象，发展了基于雷诺平均Navier-Stokes方程和自适应混合笛卡尔网格的流场数值模拟方法。通过给定合适的自适应判据，动态网格自适应技术能够实现网格的动态分布，进而准确捕捉旋涡位置并保持旋涡强度。

通过两类AVI问题的数值模拟研究，揭示了旋涡与翼型在相互作用过程中产生的复杂流动现象，包括诱导激波与二次涡的产生、非定常气动力的变化等。数值模拟结果表明，旋涡在与翼型相互作用的过程中，翼型所受的气动力产生明显的变化，而这一现象反映到三维桨叶上，即桨叶的气动力会产生变化。因此，二维AVI问题的数值模拟，对桨叶所受的载荷、振动和强度特性的研究都具有现实的指导意义。

图7 翼型升力系数与涡核位置的关系

[1] Lee S, Bershader D. Head-on Parallel Blade-Vortex Interaction [J]. AIAA Journal, 1994, 32(1): 16-22.

[2] Falissard F, Lerat A, Sides J. Computation of airfoil-vortex interaction using a vorticity-preserving scheme [J]. AIAA Journal, 2008, 46(7):1614-1623.

[3] Morvant R, Badcock K, Barakos G, et al. Airfoil-vortex interaction using the compressible vorticity-confinement method [J]. AIAA Journal, 2005, 43(1): 63-75.

[4] Tang L, Baeder J D. Adaptive Euler Simulations of Airfoil-Vortex Interaction [J]. International Journal for Numerical Methods in Fluids, 2007, 53(5): 777-792.

[5] Oh W S, Kim J S, Kwon O J. An Unstructured Dynamic Mesh Procedure for 2-D Unsteady Viscous Flow Simulations[R]. AIAA Paper 2002-0121, 2002.

[6] Oh W S, Kim J S, Kwon O J. Numerical Simulation of Two-Dimensional Blade-Vortex Interactions Using Unstructured Adaptive Meshes [J]. AIAA Journal, 2002, 40(3): 474-480.

[7] BLAZEK J. Computational fluid dynamics: principles and applications[M]. Amsterdam: ELSEVIER, 2001.

[8] LIOU M S. A sequel to AUSM: AUSM+ [J]. Journal of Computational Physics, 1996, 129(2):364-382.

[9] VENKATAKRISHNAN V. On the accuracy of limiters and convergence to Steady State Solutions[R]. AIAA Paper 93-0880, 1993.

[10] SHAROV D, NAKAHASHI K. Reordering of 3-D hybrid unstructured grids for vectorized LU-SGS Navier-Stokes calculations[R]. AIAA Paper 97-2102, 1997.

[11] Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[12] Bake W K. Dipole Sound from Cylinders[M]. Mechanics of Flow-Induced Sound and Vibration, 1st Edition, Vol.1 Academic Press, New York, 1996.

[13] Srinivasan G R, McCroskey W J, Beader J D. Aerodynamics of Two-Dimensional Blade-Vortex Interaction [J]. AIAA Journal, 1986, 24(10):1569-1576.

[14] Damodaran M, Caughey D A. Finite-Volume Calculation of Inviscid Transonic Airfoil-Vortex Interaction[J]. AIAA Journal, 1988, 26(11):1346-1353.

Numerical Simulation of Unsteady Airfoil-Vortex Interaction Problem

ZHAO Wenmei, ZHAO Jun, HU Ou

(Aviation Industry Corporation of China Helicopter Research and Development Institute, Jingdezhen 333001, China)

In the background of blade-vortex interaction problem, considering a 2D simplified model, airfoil-vortex interaction (AVI) problem, a novel numerical method for the AVI problem was developed, based on Reynolds average Navier-Stokes equations and adaptive hybrid Cartesian grid method. Given suitable adaptive criterion, the adaptive mesh refinement process could accurately capture the vortex position and keep the vortex intensity. On this basis, both unsteady flow phenomenon and aerodynamic force of two type AVI problems were investigated. Computed results showed good agreement with existing experimental data and computational results. This result has a guiding significance for researching on aerodynamic load and vibration analysis of the three dimensional blade.

hybrid Cartesian grid；adaptive, unsteady；aerodynamic force；AVI

2014-08-29

赵文梅(1986-)，女，甘肃敦煌人，硕士，助理工程师，主要研究方向：直升机旋翼动力学。

1673-1220(2015)01-006-07

V211.52

A