应用多块对接结构网格方法的直升机涵道尾桨气动特性分析

倪同兵，招启军，赵国庆，高延达，徐国华

(南京航空航天大学直升机旋翼动力学重点实验室，江苏南京 210016)

0 引言

尾桨是一个十分重要的直升机气动部件，其流场和气动特性对直升机的性能、飞行品质、噪声特性等具有重要影响，关于它的分析一直是直升机技术研究的重点之一。尾桨分为常规尾桨与涵道尾桨两种。与常规尾桨相比，涵道尾桨具有安全性更高、气动性能更好、噪声更小等优点，因此，涵道尾桨目前在多种直升机型号中得到应用。然而，由于涵道的存在，涵道尾桨的滑流形式和气动特性都发生了很大的变化，并且涵道尾桨的气动特性与其飞行状态密切相关。因此，开展涵道尾桨气动特性的研究对于直升机设计和使用具有重要的价值。

目前，用于涵道尾桨气动特性研究的方法主要有试验方法［1－3］、动量理论与涡流理论相结合的理论分析方法［4－5］，新兴的计算流体力学方法也在涵道尾桨气动特性分析方面得到了初步应用［6－10］。试验方法是开展涵道尾桨气动特性研究的直接而可靠的方法，但试验方法成本较高，周期较长，从而限制了试验方法的应用。相比于试验方法，理论分析方法研究周期短、成本较低，因此基于动量理论与涡流理论相结合的理论分析方法得到了广泛应用，但是该方法难以准确地捕捉涵道尾桨周围流场的局部特性及结构参数变化的影响。而CFD方法能充分地捕捉流场特征尤其是跨音速流场的细节，这对于涵道尾桨气动特性的细致分析尤为重要。

CFD方法经过几十年的逐步发展已经有了长足的进步，应用越来越广。国外已经应用CFD方法对涵道尾桨进行了一些研究。Alpman等［3］利用Euler方程，结合作用盘理论研究了FANTAIL(科曼奇直升机尾桨)在悬停、侧飞及前飞等状态下的流场特点与性能，但未考虑粘性影响。Ruzicka等［9］利用结构嵌套网格，通过求解N-S方程研究悬停状态下FANTAIL的气动特性，较为细致地显示了流场信息。但他们的研究只针对于悬停状态，没有涉及更为复杂的直升机侧飞和前飞状态。国内于子文等［11］基于动量源方法对简化的涵道尾桨在悬停及侧飞时的气动特性进行了一些分析，未进行前飞状态的分析。为了进一步开展涵道尾桨的流场和气动特性分析，有必要发展一套更准确、更通用且能满足一定工程实用的CFD方法。然而，建立这一分析方法的困难主要在于网格的生成、粘性影响以及适应不同的飞行状态计算。

首先，网格的合理设计和高质量生成是CFD计算的前提和基础。由于涵道尾桨结构的复杂性以及飞行状态的多样性，流场网格的生成尤为复杂，对网格的生成要求更高。为了适应直升机涵道尾桨流场的特点，本文拟采用存储简单、计算效率高、易于构建高精度格式的结构网格，并按照涵道尾桨的结构进行分块，建立了一种多块对接网格方法，在保证网格生成质量和模拟复杂外形能力的基础上，减轻了网格的生成难度，从而较大地提高了计算效率。其次，为了准确模拟涵道尾桨的流场及气动特性，必须要考虑粘性影响，因此，在求解N-S方程时，引入能够很好地捕捉涵道壁附近出现旋涡、分离等现象的一方程S-A湍流模型。同时，分块结构网格采用求解Poisson方程的方法，以保证网格的贴体性及满足附面层计算的要求。最后，为适应直升机悬停、侧飞、前飞等不同飞行状态，必须考虑建立一种通用性好的CFD方法来模拟直升机涵道尾桨气动特性。因此，本文在整个流场网格分块及生成时充分考虑了这一要求，同时将动量源方法引入涵道尾桨气动特性分析中。采用该方法在保证尾桨下洗流场本质属性的前提下，舍去求解桨叶上下表面的流动细节。与此同时，由于用围绕整个桨盘的网格来取代围绕桨叶的贴体网格，减小了网格生成的难度和网格数目，从而有效地节省了计算时间，这对建立一套通用且高效的涵道尾桨气动特性CFD分析方法很有帮助。

应用上述方法，本文建立了一套直升机涵道尾桨流场和气动特性的CFD分析方法，并通过多种算例计算验证了该方法的有效性和准确性。在此基础上，进一步针对FANTAIL涵道尾桨的不同状态(悬停、前飞、侧飞等)下的气动特性进行分析研究，拟得出在上述不同工作状态下涵道尾桨的气动特性，为直升机尾桨气动设计提供参考。

1 对接网格系统

1.1 涵道尾桨建模和网格分块

本文建立了涵道尾桨的几何模型，并根据不同飞行状态下的流场特点及湍流模型对网格的要求，建立了一套涵道尾桨对接结构网格生成方法。由于使用多块网格方法，保证了整体网格生成质量，减小了网格生成难度。在块与块的交接面上采用点对点搭接传递信息。

整个涵道尾桨流场被分割为三个块类:块1为绕中心体的网格，块2为绕外涵道的网格，其他区域为块3。为减小计算量，规定块1与块2为粘性区，在其中求解湍流N-S方程，而块3为无粘区，只需求解Euler方程。块与块的边界为交接面，如图1所示。

图1 网格分区示意图(x=0剖面)Fig.1 Block layouts for the ducted tail rotor configuration(x=0 profile)

1.2 分块网格生成

因为需要计算侧飞状态，考虑到尾桨对流场信息的影响范围，并保证网格的通用性，本文选取了较大的计算域，桨盘上下方均为30R，径向也为30R，R为尾桨桨叶半径。生成的涵道尾桨整体网格由三部分组成，分别为绕中心体的网格、绕外涵道的网格及外部区域网格。针对流场划分的三个区域，由于涵道壁、中心体均为旋转体，本文首先对其截面采用求解椭圆方程(Poisson方程)方法生成贴体网格，然后将截面网格旋转分别生成围绕涵道壁、中心体以及块3的网格。

首先推导物理平面的Poisson方程得到计算平面的Poisson方程形式，如下式:

其中p、q为方程源项，在网格生成过程中用以控制壁面和外边界的网格疏密和正交程度。其中

然后采用中心差分格式离散式(1)，得到网格迭代公式，如下:

在此基础上运用Hilgenstock源项修正方法［11］，并且采用超松弛迭代方法迭代生成围绕中心体和外涵道的O型贴体网格，如图2所示。

采用本文方法生成的分块方法一方面由于很好地控制了网格的正交性，在块交接面处网格能够光滑过渡，从而有效降低由于网格畸变引起的计算误差;另一方面，由于网格生成过程中能很好地控制壁面网格疏密，在进行N-S方程流场模拟中能够有效捕捉附面层特性。

图3为涵道的分块网格，网格尺寸分别为147×50×119，266×50×119，236×50×119，第一层网格到物面的距离为5×10－5m。

图4给出的是桨盘平面网格。在桨尖和桨根处分别进行了加密。

图2 涵道尾桨网格分解示意图Fig.2 Schematic diagram of the grid decomposition around the ducted tail rotor

图3 涵道尾桨分块网格Fig.3 Multi-block grids around the ducted tail rotor

图4 桨盘平面网格Fig.4 Grids around the tail rotor disk plane

2 流场及气动特性计算方法

2.1 主控方程

建立直角坐标系(x，y，z)下的带动量源项的三维守恒型N-S方程，如下:

其中:

q=(u，v，w)T表示在直角坐标系下的速度，粘性相关项分别为xx=2μux－ (2/3)μ▽q，xy= μ(uy+vx)，φx=uxx+vxy+wxz+ κ∂T/∂x等，ρ表示密度。E 和 H分别为总能和总焓，V为控制体单元，S为单元面积，μ、κ、T分别为粘性系数、热传导系数和绝对温度，J为动量源项。

2.2 S-A湍流模型

湍流模型采用Spalart-Allmaras［12］一方程湍流模型计算。S-A模型是从经验和量纲分析出发，经过综合分析和利用二维混合层、层流以及平板附面层的试验结果而建立起来的一方程模型。相对于代数湍流模型，S-A模型对一定范围内的分离流动的模拟能力要比代数模型更好，计算结果比较理想，故选用S-A模型。

该方程的右端分别为生成项、扩散项、耗散项和转捩相关项，本文计算中忽略了转捩相关项，则计算为全湍流状态。

封闭函数及相关常数参见文献［12］。

2.3 动量源方法

本文采用了Rajagopalan发展的动量源方法［13］，将桨叶离散成沿展向分布的若干个微段，并假设微段处桨叶的弦长、翼型厚度、剖面面积等均不变化。桨叶旋转时，向微段中心通过的网格单元中添加动量源项。选取半径r处、沿展向长度为dr的桨叶微段。在计算域坐标系(x，y，z)下，桨叶剖面相对来流速度为:

式中，v’n为桨叶剖面法速度，v’θ为桨叶剖面的周向速度。

当地流场马赫数为M=v’/a，a为当地音速，此时诱导入流角为 β =arctan(v’n/V’θ)，已知桨叶根部安装角及桨叶负扭转，可以得到当地迎角α。由当地迎角及当地马赫数，可以确定桨叶微段处的升阻力系数CL、CD。设桨叶弦长为c，可以得到桨叶微段产生的升力及阻力为:

将dL和dD向n方向和θ方向分解，得到桨叶微段的受力为dF。根据牛顿第三定律，该桨叶微段对流场的反作用力为－dF。－dF是一个瞬时的力矢量，在本文中取各片桨叶旋转一周时对同一位置的平均作用。单位时间内，N片桨叶作用于该微面下方气流的力为:

式中，△P即动量源项，将其添加到桨盘下方的网格单元中，亦即式(1)右边的J项。△P在N-S方程的迭代过程中作为解的一部分存在，也在不断更新中。

2.4 通量求解

本文通过求解湍流N-S方程建立分析涵道尾桨气动特性的方法。该方法是在文献［14-15］的基础上发展而来。求解主控方程时，空间方向采用Jameson中心格式的有限体积法对控制方程进行空间离散，时间方向上采用五步Runge-Kutta迭代法进行求解，引入人工粘性以抑制解的振荡，保证数值计算过程稳定。物面处使用无滑移边界条件，远场采用无反射边界条件。采用当地时间步长，变系数的隐式残值光顺等技术加速收敛措施。

图5给出了进行涵道尾桨流场和气动特性计算的流程图。

具体步骤如下:

步骤1:根据涵道尾桨的几何参数、飞行状态下的流场特点及湍流模型对网格的要求，生成分块网格，并将各块网格进行对接。

步骤2:根据给定的飞行状态和飞行特点，设定初始条件和边界条件，进行涵道尾桨流场的CFD求解。

步骤3:计算按块循环，每块网格中按网格面循环。判断块的编号，若是块1或块2，则判断网格面是否为桨盘平面。若是桨盘平面，则计算网格面的动量源项，并加入处于桨盘下方的网格单元中;否则，直接计算网格单元的通量。

步骤4:判断块的编号，若是块1或块2，则进行N-S方程求解;若是块3，则进行Euler方程求解。

步骤5:根据求得的流场结果，判断涵道尾桨流场计算的精度。如果满足要求，则输出结果，结束计算;如果否，返回步骤3，继续迭代计算，直到满足精度要求。

图5 涵道尾桨流场计算流程图Fig.5 Calculation flowchart for ducted tail rotor

3 涵道尾桨气动特性分析

3.1 CFD方法验证

3.1.1 N-S方程求解程序有效性验证

本文通过计算三维ONERA M6机翼的绕流特性来验证求解湍流N-S方程的CFD方法有效性。

三维ONERA M6机翼的工作状态为马赫数0.8395，攻角3.06°，雷诺数1.172 ×107。图6 是计算得到的M6机翼在剖面2y/b=0.65处上下表面压强系数分布与试验值［17］的对比情况。从图6中可以看出，本文计算得到的激波位置、强度等均与试验值吻合良好，表明了本文求解方法的有效性。

3.1.2 涵道尾桨数值分析方法的验证

为验证本文方法对求解涵道尾桨流场和气动特性的有效性，计算了FANTAIL涵道尾桨［9］和TsAGI涵道尾桨模型［8］。

图6 M6机翼2y/b=0.65剖面压强系数与试验值对比Fig.6 Comparisons on calculated pressure coefficients with experimental data of M6 wing at 2y/b=0.65

图7 悬停状态时桨盘拉力随总距角的变化Fig.7 Variation of rotor thrust at different collective pitch in hover

图8 悬停状态时总拉力随总距角的变化Fig.8 Variation of total thrust at different collective pitch in hover

图7和图8分别为悬停状态下不同总距角时，本文计算获得的FANTAIL涵道尾桨桨盘拉力、总拉力与试验值的对比。对比看出，无论是大小还是变化趋势，计算值与试验值［9］都基本吻合，这表明本文建立的涵道尾桨气动特性CFD分析方法是有效的。

图9～图12分别为悬停状态、根部安装角为40°时，计算得到的TsAGI涵道尾桨模型桨盘处轴向诱导速度分布、桨叶展向载荷分布、桨盘处周向速度分布以及桨盘处诱导速度分布。贴体网格计算值和涡流方法计算值来自文献［18］。可以看出，计算结果与试验值或贴体网格计算值比较接近，从而进一步验证了本文建立的涵道尾桨气动特性CFD分析方法的有效性。

3.2 涵道尾桨气动特性的计算与分析

在以上验证本文CFD方法对涵道尾桨流场和气动分析有效性的基础上，进一步开展悬停、侧飞及前飞等状态下FANTAIL涵道尾桨气动特性的计算。

图13计算了不同飞行状态时桨盘拉力随总距角的变化。可以看出，随着总距角的增大，桨盘拉力均增大。

图9 TsAGI涵道尾桨轴向诱导速度分布Fig.9 Spanwise distributions of the induced downwash of TsAGI ducted tail rotor

图10 TsAGI涵道尾桨桨叶沿展向载荷分布Fig.10 Spanwise distributions of the blade sectional thrust of TsAGI ducted tail rotor

图11 TsAGI涵道尾桨的桨盘处周向速度分布Fig.11 Spanwise distributions of the circumferential induced velocity of TsAGI ducted tail rotor

图12 TsAGI涵道尾桨桨盘处诱导角分布Fig.12 Spanwise distribution of the induced angle-of-attack of TsAGI ducted tail rotor

图13 不同飞行状态时桨盘拉力随总距角的变化Fig.13 Variation of rotor thrust at different collective pitch at different flight states

图14计算了悬停状态根部安装角为40°时涵道尾桨拉力系数残值收敛曲线。可以看出，桨盘拉力系数只需较少的迭代步数就可以获得收敛结果，而涵道拉力系数振荡比较剧烈。计算得到的涵道壁产生拉力占总拉力的46%左右，与文献［7］吻合。

图15为计算的悬停状态，总距角为25°时y=0剖面的涵道尾桨合速度云图及流线。可以看出，由于强逆压梯度的作用，在涵道唇口处发生了气流分离，存在回流区。中心体下方，涵道上下方都有旋涡。图16为此时外涵道及中心体表面压强云图。在唇口处，气流速度较大，静压较小，这是涵道产生拉力的主要原因。

图14 悬停残值收敛曲线Fig.14 Residual convergence curve in hover

图17给出了前飞速度v=77.17m/s、总距角为20°时y=0剖面速度矢量图。由图17可看出，在涵道尾桨的抽吸作用下，气流发生了偏折。由于高速气流直接“撞击”到涵道壁面，涵道周围的流场比较复杂，出现了很多旋涡。图18则为涵道尾桨表面压强云图，自由来流沿X轴正向，处于气流运动方向下方的涵道内壁面上出现了很大的静压，这是产生涵道拉力的主要原因。

图15 悬停状态合速度矢量图Fig.15 The velocity vector contour map in hover

图16 悬停状态下涵道尾桨表面压强云图Fig.16 Pressure contour map in the surface of ducted tail rotor in hover

图17 前飞状态下速度矢量图Fig.17 Velocity vector map in forward flight

图19和图20分别给出了左侧飞(侧飞速度为23.15m/s，总距角为15°)时FANTAIL 涵道尾桨 y=0 剖面的流线和速度矢量，可以看出，由于侧飞速度较大，而总距角相对较小，气流刚通过涵道出口就会迅速向上偏折，在外涵道上下，中心体附近出现了很多的旋涡。

图21和22给出了右侧飞(侧飞速度为23.15m/s，总距角25°)时FANTAIL涵道尾桨y=0剖面的压强云图及速度矢量图。可以看出，由于气流具有比较大的来流速度，外涵道下方及侧面有很多回流区。涵道唇口处压强较小，由于气流直接“撞击”，涵道及中心体上方压强较大。

图18 前飞状态下涵道尾桨表面压强云图Fig.18 Pressure contour map in the surface of the ducted tail rotor in forward flight

图19 左侧飞时压强云图Fig.19 Pressure contour map in left sideward flight

图20 左侧飞时速度矢量图Fig.20 Velocity vector map in left sideward flight

图21 右侧飞时剖面压强云图Fig.21 Pressure contour map in right sideward flight

图22 右侧飞时速度矢量图Fig.22 Velocity vector map in right sideward flight

4 结论

本文建立了一个基于多块对接结构网格和动量源模型的直升机涵道尾桨流场及气动特性数值分析方法，通过多种算例进行流场及气动特性的影响分析，得出了结论如下:

(1)计算得到的FANTAIL涵道尾桨悬停状态下不同总距角时，桨盘拉力和总拉力同试验值吻合良好，速度分布、载荷分布等计算结果也与试验结果非常吻合，表明本文方法可有效地用于涵道尾桨的气动特性模拟。

(2)在悬停状态，随着总距角的增大，桨盘拉力和总拉力均增大。涵道唇口处气流加速与逆压梯度共同作用，是产生涵道拉力的主要原因。

(3)对于右侧飞状态，由于具有侧风速度，桨盘处相对入流速度增大，桨盘产生的拉力减小。与右侧飞相反，左侧飞时，桨盘处相对入流速度减小，产生的桨盘拉力增大。

(4)在前飞状态时，由于涵道尾桨的抽吸作用，流过涵道尾桨的气流会发生偏折。当高速气流“撞击”到涵道壁面时，涵道周围的流场会出现许多旋涡。

［1］KEYS C，SHEFFLER M，WEINER S，et al.LH wind tunnel testing:key to advanced aerodynamic design［A］.47th A-merican Helicopter Society Annual Forum ［C］，Phoenix，1991.

［2］AKTURKL A，SHAVALIKUL A，CAMCI C.PIV measurements and computational study of a 5-inch ducted fan for V/STOL UAV Applications［R］.AIAA-2009-332.

［3］李建波，高正，唐正飞，等.涵道风扇升力系统的升阻特性试验研究［J］.南京航空航天大学学报，2004，36(2):164-168.

［4］BOURTSEV B N，SELEMENEV S V.Fan-in-fin performance at hover computational method［A］.26th European Rotorcraft Forum［C］，2000.

［5］徐国华.直升机涵道尾桨与孤立尾桨的气动特性对比研究［J］. 空气动力学学报，1995，13(4):420-426.

［6］RAJAGOPALAN R G，KEYS C N.Detailed aerodynamic analysis of the RAH-66 FANTAILTM using CFD［J］.Journal of the America Helicopter Society，1997，42(4):310-320.

［7］ALPMAN E，LONG L N，KOTHMANN B D.Understanding ducted-rotor anti-torque and directional control characteristics part I:steady-state simulation ［J］.Journal of Aircraft，2004，41(5):1042-1053.

［8］LEE D H，KWON O J，JOO J.Aerodynamic performance analysis of a helicopter shrouded tail rotor using an unstructured mesh flow solver［A］.5th Asian Computational Fluid Dynamics Forum［C］，Busan，Korea，2003.

［9］RUZICKA G C，STRAWN R C，MEADOWCROFT E T，Discrete blade CFD analysis of ducted tail fan flow［R］.AIAA-2004-0048.

［10］于子文，曹义华.涵道尾桨的CFD模拟与验证［J］.航空动力学报，2006，21(1):19-24.

［11］张正科，庄逢甘，朱自强.两种椭圆型方程求源项方法在喷管内流场网格生成中的应用［J］.推进技术，1997，18(2):95-97.

［12］SPALART P R，ALLMARAS S R.An one-equation turbulence model for aerodynamic flows［R］.AIAA-92-0439.

［13］ZORI L A J，RAJAGOPALAN R G.Navier-Stokes calculations of rotor airframe interaction in forward flight［A］.The 48th Annual Forum of the American Helicopter Society［C］，Washington D C，1992.

［14］叶靓，招启军，徐国华.非结构嵌套网格的直升机旋翼/机身前飞流场数值模拟［J］.航空动力学报，2009，24(4):903-910.

［15］ZHAO Q J，XU G H，ZHAO J G.New hybrid method for predicting the flowfield of helicopter in hover and forward flight［J］.Journal of Aircraft，2006，43(3):372-380.

［16］HELLSTROM T，DABIDSON L，RIZZI A.Reynolds stress transport modeling of transonic flow around the RAE2822 airfoil［R］.AIAA-94-0309.

［17］SCHMICT，CHARPIN.Pressure distributions on the ONERA M6 wing at transonic mach numbers［R］.Experimental Data Base for Computer Program Assessment，Appendix B1 in AGARD-AR-138.1995.

［18］BOURTSEB B N，SELEMENEV S V.Fan-in-fin performance at hover computational method［A］.The 26th European Rotorcraft Forum［C］，2000.