谐波平衡法在叶轮机非定常数值模拟中的应用

张恒铭, 唐弋棣, 卫海洋, 林彬彬

(1.中国民航飞行学院航空工程学院，广汉 618300；2.L.A.Turbine涡轮机械公司，瓦伦西亚 91355；3.北京动力机械研究所，北京 100074)

叶轮机中的流动现象具有强烈的非定常性，而引起这种非定常性的因素很多。其中，由叶排相对运动和流动参数周向不均匀性所引起的非定常排间干扰是研究的热点。排间干扰主要包含了周向静压不均匀所引起的势扰动和尾迹扰动。这些扰动对压气机的性能有显著的影响。中外的学者在非定常流动方面开展了大量的研究[1-5]。他们的研究表明排间干扰对叶轮机性能产生了不可忽略的影响，例如尾迹作用于下游流场时不但会与下游的边界层、尾缘涡和端壁二次流等发生相互作用，且自身还会在下游的输运过程中发生拉伸、变形，形成“耗散”和“恢复效应”。排间扰动还是发动机噪声的最主要来源之一，研究表明叶片排本身的湍流现象(上游尾迹、角涡和叶尖泄漏涡等)作用于相邻叶排时，会诱发湍流干涉噪声[6-8]。此外，排间扰动还会影响叶轮机叶片的气弹特性，诱发颤振，引起受迫响应问题[9-10]。因此，叶轮机的排间干扰效应十分重要，要进行相关的数值研究，必须进行非定常的数值仿真。

然而，对于叶轮机来说，常规的非定常计算通常采用时域的直接推进方法，往往需要进行多通道，甚至全通道的数值模拟[11-12]，计算量巨大，对计算资源的消耗往往比定常模拟高两个数量级以上。因此，时域方法虽然比较精确，可以捕捉扰动与主流之间的非线性作用。但时域方法的耗时问题，限制了其在工程中的应用。

另一类可用于叶轮机非定常求解的方法是频域方法，频域方法通常是指将非定常变量表达为周期性函数的形式，再通过推导和变换最终将方程组表达为频域形式的方程组，从而可以对非定常流场进行求解，其最大的优点在于求解速度快。目前，频域方法中，使用较为广泛的有He等[13]提出的非线性谐波方法(nonlinear harmonic method，NHM)，Hall等[14]提出的谐波平衡方法(harmonic balance method，HB)以及Gopinath等[15]提出的时间谱方法(time spectral method，TSM)。鉴于频域方法的求解效率较高，其已经被一些研究人员用来考察叶轮机中的非定常现象[16-17]。但是频域方法的精确性不如时域方法，精确度与采用谐波的数量关系较大。同时，采用过多的谐波数会带来计算资源的浪费。因此，谐波数量的合理选取是一个值得研究的问题，目前，这方面的研究较少。

将谐波平衡方法引入到叶轮机的非定常数值模拟中，并以一对转风扇为研究对象，通过对比常规的时域非定常结果，考察了谐波数量对谐波平衡法精度和可靠性的影响。

1 物理模型和计算方法

叶轮机中的非定常数值仿真方法可分为时域非定常和频域非定常两种类型。时域非定常算法计算精度相对较高，但计算量巨大，往往难以被工程应用所接受。而频域非定常算法虽然会损失部分精度，但其可以大大减少计算资源的消耗。频域非定常计算方法主要有时间线化方法、非线性谐波方法和谐波平衡法。其中，谐波平衡法能够捕捉主流和扰动之间的强烈非线性扰动，且前提假设较少，方程简单易于编程，进而得到广泛应用。

1.1 谐波平衡法

谐波平衡方法是一种用于求解周期性非定常流动的高效方法。该方法的核心思想是将变量分解成总体平均项和时间扰动项，然后再将扰动项以傅里叶级数进行近似。非定常流动变量可表示为

(1)

式(1)中：W为非定常守恒变量;W0为总体平均量;An和Bn为第n个谐波的傅里叶系数;ωn是第n个谐波的角频率。令W1，W2,…，W2N+1为2N+1个时刻的变量值，则有式(2)成立:

(2)

式(2)中:

(3)

W*=[W1W2…W2N+1]T

(4)

(5)

对于变量W*，其满足N-S方程，写成半离散形式为

(6)

式(6)中：V为网格单元体积；把式(2)代入式(6)有：

(7)

(8)

式(8)即为需要进行求解的谐波平衡控制方程，τ为虚拟时间。容易知道，相对于定常控制方程，其多了一个由谐波产生的源项VDW*，且其控制方程数目是定常的2N+1倍。因此，理论上讲，其计算量略高于定常的2N+1倍，相比时域非定常计算，其效率明显高得多。

1.2 双时间步方法

时域非定常计算，采用的是双时间步法，对于半离散形式的方程，在时间上应用二阶的向后差分有

(9)

为对式(9)进行求解，引入了虚拟时间步，将物理时间层的替换为虚拟时间层的,并添加虚拟时间推进导数项，得：

(10)

式(10)即为用于求解时域非定常的控制方程，虚拟时间步的迭代称为内迭代，当计算收敛时，Wn+1=Wm。

1.3 计算方法

在空间离散上，采用了中心差分的JST格式，为二阶精度，其计算效率较高，稳定性好，在工程中已经被广泛应用和验证。时间推进上采用了隐式方法，采用了一阶L-F分裂格式进行了对流通量的线化处理，黏性项线化进行了简化，保留了对角项。此外，对于谐波控制方程，还需对谐波源项进行线化处理，这里参考Sicot的处理方式[18]，完成这些线化以后就得到了求解的线性代数方程组，对整个线性代数方程组的求解采用了上下对称的逐次超松弛迭代方法(LUSSOR)。对于隐式离散的方程组，系数矩阵进行分裂后为

(L+D+U)ΔW*=-R*-VDW*-VDΔW*

(11)

式(11)中：L为下三角矩阵;U为上三角矩阵;D为对角阵。LUSSOR迭代进行求解时的每一步包含了Smax次的向前扫掠和向后扫掠。即:

(12)

式(12)中:s为最内层迭代;l为次内层迭代;m为外部迭代。其具体的算法过程如图1所示。

图1 谐波Lussor算法流程

对于双时间步法的非定常计算，方程的空间离散、数值方法跟谐波平衡法保持一致。而对应的离散线性方程组，也采用LUSSOR方法进行迭代。不同的是，双时间步方法中由于没有谐波源项，不需要对谐波源项进行特殊处理，因此，其算法流程中没有图1中的次内迭代部分(深色背景部分)，即l层的迭代。

2 研究对象及计算设置

研究对象为一对转风扇，该风扇的设计流量为230 kg/s，风扇进口处的外径为0.635 m，轮毂比为0.247。前转子转速为5 200 r/min，叶尖切线速度是345.8 m/s，后转子转速为-4 100 r/min，叶尖切线速度为-272.6 m/s。前转子的设计压比为1.53，后转子的设计压比为1.45，总设计压比为2.218 5。如图2所示为对转风扇的三维示意图。

图2 对转风扇三维图

时域计算采用了双时间步法，物理时间步取为后转子所受前转子扰动周期的1/20，其也是前转子所受扰动周期的1/30。而对于谐波平衡方法，按照采用谐波数目的不同共进行了7次数值模拟，其上下游扰动谐波数均保持相同，谐波数从1取到7，这样以分析随谐波数变化时，谐波方法精度的变化。

由于在周期性边界上使用了相位滞后条件，因此只需要使用一个叶片通道网格。而在固壁上应用了壁面函数来处理边界层的流动，其对壁面第一层网格大小的依赖性较低，则H形网格即可以满足数值计算的需要，对转风扇的三维网格如图3所示。

图3 对转风扇三维网格

网格在周向和径向的网格数分别为57、65，两排网格数总共约80×104。叶尖间隙部分采用了削尖处理，间隙在径向上的网格数为4。

谐波平衡法的计算中采用了并行策略，鉴于谐波方法的良好并行特性，在计算资源足够的情况下，其对计算时间的消耗随谐波数的增加几乎不变，完成对转风扇的非定常数值仿真耗时约为2.5 h。而同样采用并行计算的时域非定常数值模拟则消耗了约28 h。即时域非定常的耗时约为谐波平衡方法的11倍，这证明了谐波平衡方法的高效性。

3 计算结果和分析

排间干扰中非常重要的一种扰动是上游尾迹的扰动，熵等值线图可以用来显示尾迹的传播。图4和图5给出了20%和80%叶高的熵等值线对比。从图4、图5中可以看出，随着谐波数的增加，上下游转子通道交界面处熵等值线的连续性变得越来越好，对于20%截面来说，当谐波数为1时，尾迹在交界面处的连续性最差，后转子计算域中交界面附近处的尾迹形状几乎不可辨识，随着谐波数提高到2，交界面处的尾迹逐渐清晰可辨，但连续性还是稍差。谐波数到3时，其连续性也已经变得很好，之后，再随着谐波数增加，其连续性会继续变好，但整体的变化已经不大。而对于80%叶高处的熵图来说，当谐波数为1时，虽然交界面附近传播到下游的尾迹比较清楚，但是尾迹的宽度和强度都跟上游计算域有明显差别，跟20%截面的情况类似，当谐波数到3时，其连续性也得到了很大的改善，之后随着谐波数的增加，其连续性的提升也并不大。

此外，通过对比还可以发现，谐波计算和时域计算的尾迹在下游叶排通道中的形态基本保持一致。对于谐波数为1的情况，由于谐波数较小，对尾迹的形状刻画得不够细致，因此使得尾迹宽度模拟得较大，但对于谐波数3甚至更多的情况，其吻合得很好，这就证明了谐波方法的有效性，同时，这也跟Hall等[14]及Sicot等[19]的验证分析相吻合。

为了进一步分析谐波方法的准确性，图6、图7给出了前、后转子叶尖间隙和50%叶高尾迹附近处的静压变化对比。影响静压变化的扰动因素很多，势扰动、尾迹扰动都很大程度影响着其变化规律。

URANS表示时域非定常，HB N=1，2,…，7表示谐波数

图5 80%叶高瞬时熵等值线图

对于前转子来说，其仅受势扰动的影响，因此其变化规律相对来说比较单一。这也能从图6中反映出来。无论对于叶尖间隙处还是尾迹附近处的静压变化，其都是一个类正弦曲线。对于间隙内的监测点来说，其变化曲线的形状随着谐波数的提升并没有发生大的变化，只是其平均值逐渐上升，变得越来越接近时域的变化曲线，分析认为其主要原因在于随着谐波数的提高，第二排计算的精确程度越来越高，由此影响了第一排静压的平均值，但其随时间的变化规律基本没有发生改变。而对于前转子叶中尾迹附近的监测点来说，随着谐波数提高，其变化的规律有一些细微的改变，而平均值则没有太大的波动，而当谐波数达到3以上时，静压变化规律基本不再随谐波数的变化而改变。总体来说，前转子所受排间扰动的影响较为简单，用较少的谐波数就可以达到比较满意的结果，尤其对于叶中区域来说，2个谐波数已经基本可以满足工程应用要求。

后转子所受的扰动相对来说较为复杂，势扰动、熵扰动和涡扰动同时存在于后转子通道内，对比图6、图7也能发现其变化规律更为复杂。其次，后转子的静压变化幅值也比前转子大很多，尤其对于间隙部分的流动来说，几乎大一个数量级。

对于间隙内的静压变化来说，谐波数为1时的结果跟时域的结果依然相差较大。而当谐波数达到2时，计算结果准确度大大提升，而当谐波数达到4以后整个静压变化曲线形态已跟时域非常接近。总体来说，除了谐波数为1的情况，其静压变化相对于时域的误差并不大，而最大误差出现在0.4周期处，该处的变化是一个周期内变化最复杂的地方。这表明对于比较复杂的变化，谐波平衡方法更容易产生误差，但其变化的趋势还是跟时域结果保持了一致。

图6 前转子监测点的静压变化

图7 后转子监测点的静压变化

对于后转子尾迹附近的监测点，从图6(b)可以看出，其变化规律在所有监测点中是最为复杂的，它的变化曲线跟正弦形态差别很大。因此，谐波的计算误差会较大，对比时域的变化就可以发现这种情况，即使对于谐波数为7的情况，其仍有一定的误差，但这种误差相对于平均静压值来说并不算大(不到5%)；另一方面，谐波模拟的静压变化趋势基本能够与时域保持一致，且当谐波数达到2或3时，其变化规律就跟时域的静压变化规律达到了较好的相似。Gopinath等也发现了类似的现象[20]，并认为用谐波平衡方法来研究排间的干涉效应是可行的。

4 结论

将谐波平衡方法引入到了叶轮机械的非定常数值模拟中，通过谐波源项的隐式线化处理和对应隐式迭代方法的改进，实现了叶轮机的高效非定常数值模拟。

然后将谐波平衡方法应用到了某对转风扇的非定常数值模拟中，结果发现谐波平衡方法对计算时间的消耗仅为时域非定常计算的9%。进一步对比时域非定常的计算结果，分析了谐波平衡方法的有效性和可靠性。结果发现，谐波平衡方法能够较好捕捉排间的非定常扰动信息，能够清晰刻画势扰动和尾迹的传播过程。当谐波数达到3或者以上时，其计算精度就基本可以满足工程应用的需要。