基于URANS/TPP 模型的风扇转静干涉单音噪声预测

2024-01-20蔡堉楠纪良李旦望夏烨倪臻

航空学报 2023年22期

蔡堉楠，纪良，李旦望，夏烨，倪臻

中国航发商用航空发动机有限责任公司先进技术研究部，上海 200241

航空发动机噪声是飞机最主要的噪声源，对其进行降噪设计是“安静型”飞机研制的关键。对于当代大涵道比构型发动机，风扇噪声占据主导［1］，其转子叶片的非定常周期性尾迹与出口导流叶片（Outlet Guide Vane，OGV）干涉会产生强烈的单音噪声。对于亚声速工况，风扇转静干涉噪声通常决定了噪声窄带频谱中声压级的最高值，即使是在超声速情况下，转静干涉噪声也依然是风扇离散噪声的重要组成部分［2］。随着对噪声抑制机理认知的深入和设计方法的进步，通过精细化叶型设计实现声源降噪成为进一步提升风扇低噪声性能的有效途径［1，3-4］，代表性的手段包括风扇转子低噪声设计［5-7］、弯掠静子设计［8-12］、仿生叶片设计［13-15］等。为了实现对设计的支撑，要求预测模型能够反映详细设计参数对声源强度的影响。

风扇噪声的产生伴随着流场和声场的耦合效应，给其声源的计算带来了极大挑战。近几十年来，该领域学者采用了计算流体力学（Computational Fluid Dynamics，CFD）方法和计算气动声学（Computational Aeroacoustics，CAA）方法对三维流场/声场进行真实建模［16-19］。由于由转静干涉产生的离散单音噪声扰动量级很高，而且只在特定的频率上出现，因此为采用非定常雷诺平均数值模拟（Unsteady Reynolds Average Numerical Simulation，URANS）模型的CFD 方法求解这类问题提供了可能。相比于直接模拟或者大涡模拟方法，URANS 模型计算量大幅降低，同时该模型能够刻画叶型设计参数对噪声的影响［20］。因此，URANS 模型成为目前工程中叶轮机单音噪声源预测的有效手段［21-24］。另外，风扇流动具有时-空周期性特征，非线性谐波（Non-Linear Harmonic，NLH）方法利用截断的傅里叶级数对流场进行逼近，可通过单通道计算实现非定常流场预测，具有比URANS 模型更高的计算效率［25-26］。但是当考虑周向非均布静子构型、进气攻角效应或者安装效应等存在非周期情况时，需要开展对所有风扇叶片数值建模的全环URANS 计算。

由于声波具有小尺度特征，其对数值方法的离散误差敏感，因此如何在计算量可接受的前提下保证声场预测精度，是进行风扇噪声预测的关键［27］。同时，风扇噪声在管道内声场以模态波形式传播。为了获取声源特性，需要从声源的脉动量中提取声模态信息［28-30］。

本文主要工作是针对大涵道比风扇构型，在对数值方法进行误差分析的基础上，基于URANS 模型进行风扇三维非定常流场模拟，结合三平面压力（Triple Plane Pressure）模态匹配方法（简称TPP方法）从流场中提取声源模态信息，并采用部件级声学试验数据验证了预测方法的可靠性。此外，通过本文数值模型，对比验证了弯掠静子叶型设计的降噪效果，分析了转静干涉噪声抑制基本规律。

1 预测模型

本文的风扇单音噪声预测模型以Navier-Stokes 方程和管道声模态传播方程为基础，结合URANS 模型和声源提取方法，实现了转静干涉单音噪声的预测。

1.1 基于URANS 模型的非定常流场计算

为了获得声源特性，需要通过CFD 方法计算获得风扇流场中的脉动压力场。描述非定常黏性流动的控制方程为可压缩的Navier-Stokes 方程组。在URANS 模型中，将Navier-Stokes 方程组进行时间积分平均，得到一组以时均物理量和脉动量乘积的时均值作为未知量的非封闭方程。本文采用k-ε模型，基于涡黏性假设，建立雷诺应力项与湍动能和湍流耗散率等物理量的关系，使方程组封闭。基于数值求解方法，能够获得流场压力脉动信息，从而为声源提取提供输入。

1.2 TPP 方法基本方程

对于声源信息的提取，Wilson［31］最早提出波分离（Wavesplitting）方法，通过在一个交界面上对压力和速度进行模态分解，实现对传播波和反射波的分离。基于该思路，Ovenden 和Rienstra［28］提出TPP 方法，实现了声模态的有效分解，同时可分辨传播波和反射波。近些年来，相关学者对TPP 方法做了进一步扩展，Wohlbrandt 等［30］提出了XTPP 模型，能够分离对流压力扰动对声源信号的影响；Guérin［32］提出改进TPP 模型，进一步考虑了旋流和径向流动的影响。

本文采用Ovenden 和Rienstra 的TPP 方法［28］，将发动机的进气道和外涵道简化成圆形或者环形管道，基于波动方程，某个频率的声波在管道内的传播可采用以下形式：

式中：z、θ、r分别表示轴向、周向和径向坐标；m和n分别为周向和径向模态阶数；为声模态幅值；为轴向波数；Ψ(κmnr)为径向分布函数；f为声波频率；c为声速；κmn为管道特征值；A和B为常系数；Jm和Ym分别为m阶的第1 类和第2 类贝塞尔函数。其中，上标“+”表示顺流传播，例如风扇外涵道出口的后传噪声；上标“-”表示逆流传播，例如风扇进气道的前传噪声。

常系数A和B的计算公式为

式中：上标“′”表示求导；rD表示管道半径。

TPP 方法通过设置3 个监测面，能够在非定常流场中提取周向、径向模态阶数及对应的幅值。对于某一周向模态m，式（1）可写为

式中：Pn为第n个声模态幅值；Ψn表示径向分布函数；Φn为轴向分布函数；i=1，2，3 表示TPP监测面序号。通过在式（6）等号左右两边分别乘以第j个模态的径向分布函数的共轭，并对径向取积分可得

式（7）可改写成以下形式：

定义代价函数：

其中：

为了使代价函数值最小，须满足以下条件：

式中：下标l表示目标声模态的序号。由于Pl为复数形式，因此有

根据式（12）和式（13）可得

根据式（14），得到最终形式：

式中：χln表示系数矩阵。式（16）共有N++N-个未知数。如图1 所示，P1、P2和P3分别为3 个监测平面，其中上标“+”和“-”分别表示沿轴向和逆轴向传播，且n≠0。此时，未知数个数与方程个数相同，通过求解式（16），可得到各阶模态对应的幅值，从而完成声模态提取。

图1 TPP 方法示意图Fig.1 Sketch of TPP method

对于（m，n）阶模态对应的声功率，可由式（17）得到。

式中：ρ0为流体密度；c0为流体声速；R和r分别为环形管道的外径和内径；k0和kmn分别为总波数和轴向波数。

2 声场计算误差分析

在管道中传播的声场相对于流场扰动值是小量，对数值离散产生的耗散误差极为敏感。而离散误差与空间和时间步长相关。过大的网格尺寸会导致计算失真，而网格过密会使计算量太大。因此，在保证计算可靠前提下尽量降低网格数量是本文模型进行工程级噪声预测的基础。

采用模态波在环形管道中传播的算例，讨论网格尺寸与时间步长对声场误差的影响。根据模态波传播特征，定义单位波长点数NPPW：

式中：Nz和Nm分别为单个波长轴向和周向投影距离所包含的网格点数。采用单频（3，1）模态声波入射，考虑Ma分别为0.5、0.3、0.1 和-0.3 时的工况，选取不同NPPW计算传播单位波长后的声压级误差ΔSPL，结果如图2 所示。可以看到，随着网格的加密，声压级误差逐渐减小，其与单位波长点数的拟合公式为

图2 网格尺寸与误差拟合曲线Fig.2 Fit line of mesh size and error

下面讨论时间步长的影响。定义无量纲单位周期时间步长数NT为关注单音声波的周期T与CFD 物理时间步长dt的比值，即

数值计算的耗散误差与NT存在相关性。图3为NT与声压级误差ΔSPL 的关系。可以看到，对于不同的单位周期时间步长数，在NT大于60 时，单位周期时间步长数对声压级误差的影响很小，并且继续增加单位周期时间步长数对精度提升有限，综合考虑计算经济性和准确性，在以下的计算中均采用NT=60。由于风扇单音噪声通常包含多个频率单音，越高频的成分要求的物理时间步长越小。因此在实际计算中，式（22）中T取关注的最高频率单音周期，以保证各频率单音分量的准确计算。

图3 时间步长与误差拟合曲线Fig.3 Fit line of time step and error

通过以上分析，验证了URANS 模型具备声场计算可行性，同时提出空间和时间步长选取方法，以指导风扇噪声计算。

3 风扇转静干涉单音噪声预测验证

采用国内设计的某型大涵道比风扇缩尺试验件构型，风扇的气动和声学试验在德国 Wildau 的AneCom 航空试验（AneCom Aero Test，ACAT）噪声试验台（如图4 所示）上完成。

图4 AneCom 航空试验噪声试验台Fig.4 AneCom AeroTest noise test facility

试验件的风扇转子叶片数为18，OGV 数为55。在进气道和外涵出口的管道外壁面分别布置1 圈100 个麦克风，采用壁面齐平安装方式并且周向非均布，用以进行声源周向模态分解。同时，在外涵出口布置旋转桶麦克风，旋转桶中周向等距设置有3 列麦克风，每列60 个且轴向等距分布，用以进行声源周向和径向模态分解。通过试验测得给定工况下的麦克风声压脉动信号和风扇轴的转速脉冲信号，结合转速脉冲信号对声压脉动信号开展重采样，并按叶片通过频率（Blade Passing Frequency，BPF）的整周期倍截断，以准确获取声压测点在BPF 位置的声压幅值和相位特性。在给定的频率下，管道内各阶模态在测点位置的声压特性是1 组线性无关基向量，基向量的线性组合即为测得的声压幅相特性。因此，利用最小二乘法求解测得的声压在给定声模态构成的线性空间下的坐标，即可得到各阶声模态在给定频率下的幅值，从而可获取声模态，这为噪声数值评估方法提供了试验数据对比验证。

风扇试验件直径约为0.45 m。管道声模态结果与风扇/增压级的转子、静子数相关，若改变叶片数进行约化计算，预测的模态阶数将与实际情况不符，因此本文采用全环计算方案。风扇转静干涉噪声主要由风扇转子与外涵出口导叶干涉产生，且为了尽量减少网格点数，对于内涵仅考虑第1 级静子导叶S0，风扇转子叶片数为18，OGV 数为55，S0 叶片数为67。风扇进口设置总温总压，外涵和内涵出口设置为静压出口。采用Numeca AUTOGRID 完成风扇部件网格划分，流场计算及后处理采用ANSYS CFX 完成。

计算网格设置如图5 所示，在风扇转子前进口和OGV 后出口通过拉伸网格的方法设置耗散段，以实现无反射边界。对于计算网格，根据式（21）可知，保证预测的最小目标波长的NPPW应大于25，同时每个径向波长内网格点数应大于15，最小目标声波周期的时间步长数NT为60。在计算网格前端和后端分别设置TPP 监测面，为前传和后传声源提取提供输入。采用URANS方法进行非定常流场计算，湍流模拟选用k-ε模型，在风扇转子和外涵静子表面y+接近于30，总网格量约为5×107，其中y+为无量纲数，表示网格节点到壁面的距离与流动特征尺度之比。

图5 计算网格设置Fig.5 Computing mesh settings

计算了适航规定的边线（OP1，103.7%）、飞越（OP2，93.9%）和进场（OP3，60.0%）转速工况。首先进行定常计算，将其结果作为非定常模拟的初场以保证计算稳定性。非定常计算总计时间步数约为4 000 步，在每个时间步内进行10 次内循环。采用了448 核CPU 并行运算，总计算时间为341 h。

图6 为103.7%转速下非定常计算迭代过程某一轴向位置不同径向高度处脉动压力随时间变化曲线。在迭代大于2 800 步时，脉动压力趋近于稳定周期解，表明非定常流在计算过程中已趋近于收敛。

图6 非定常压力随时间变化曲线Fig.6 Time history of fluctuating pressure

观察外涵出口处流量和压比值并与试验数据对比。如图7 所示，3 个适航工况流量误差分别为0.58%、-0.21% 和0.41%，压比误差为2.77%、0.12% 和-0.10%。通过试验对比分析，表明流场计算的正确性。

图7 流量和压比预测值与试验数据对比Fig.7 Comparison of prediction and experiment results of mass flow and pressure ratio

图8 为80%叶高处转子和静子间的静熵云图，图8 中采用了相同的云图标尺。可以观察到转子的尾迹进入外涵后与下游静子产生相互作用，它是非定常压力扰动的重要来源之一。转子尾迹被静子前缘分割，并与叶片表面边界层相互作用，在叶片表面产生高熵增区域，造成边界层的损失增加。在静子下游，转子尾迹与外涵通道中的主流继续发生掺混，这种影响一直发展到外涵出口。

图8 不同转速工况瞬时静熵云图Fig.8 Instantaneous contours of static entropy of different operation conditions

观察风扇静子后方外涵道内BPF 的瞬时脉动压力实部。从图9 中可以看到，由于转静干涉产生的脉动压力在周向和径向上呈明显周期分布，符合管道声模态波特征。需要指出的是，对于OP1 和OP2，其1 阶和2 阶BPF 均导通，但对于低转速工况OP3，1BPF 截止，因此观察了OP3的2BPF 和3BPF 下的脉动压力。

图9 OGV 出口脉动压力实部云图Fig.9 Contours of real part of fluctuating pressure at OGV outlet

下面对边线（OP1，103.7%）、飞越（OP2，93.9%）和进场（OP3，60.0%）3 个工况的后传噪声和进场（OP3，60.0%）工况前传噪声进行声模态提取。在对监测面时域压力脉动数据进行频域转换后采用TPP 方法进行声模态分解和声功率计算。最终，采用声模态阶数和幅值信息的形式，给出风扇转静干涉单音声源。

根据Tyler 和Sofrin［33］的转静干涉噪声理论，声波周向模态阶数m符合以下特征：

式中：j为谐波阶数；k为整数；F为转子叶片数；S为OGV 数或者S0 叶片数，分别代表转子与外涵OGV 或者第1 级内涵静子导叶的干涉噪声模态。

由于URANS/TPP 模型的限制，本文仅对BPF 下的主要占优模态进行讨论。图10 和图11为不同转速下风扇后传声模态与AneCom 试验数据的对比，干涉噪声周向模态阶数和误差如表1 所示。为了保证与试验数据一致，OP1 和OP2 工况声模态幅值采用声功率表示，OP3 采用声压级表示。可以观察到，由转静干涉产生的主要导通单音模态被准确捕捉到，声模态阶数试验结果一致。其中，在OP1、OP2 和OP3 后传工况，由转子与OGV 干涉产生的模态占主导；在OP3前传工况3BPF 下，除了转子/OGV 干涉产生的m=-19 模态，预测模型还捕捉了由转子/S0 干涉产生的m=-13 模态。计算出各工况下幅值绝对值平均误差为4.7 dB，验证了本文方法对占优单音噪声预测的有效性。

表1 各工况占优声模态阶数和预测误差Table 1 Order and prediction error of dominant acoustic modal under different operation conditions

图10 不同工况下后传声模态与试验数据对比Fig.10 Comparison of prediction data with experiment data under different operation conditions for rear acoustic mode

图11 OP3 前传声模态与试验数据对比Fig.11 Comparison of prediction data with experiment data at OP3 for front acoustic mode

4 低噪声弯掠静子设计的数值分析

风扇静子OGV 弯掠造型是抑制转静干涉噪声的有效手段。弯掠叶片主要增加了上洗速度沿径向的相位变化，从而使不同径向位置的非定常载荷相互抵消以降低声源强度。图12 为弯掠OGV 示意图，从叶根到叶尖，当叶片前缘随径向位置增大而向下游变化时，叶片后掠，掠角α取为正；反之，叶片前掠，掠角为负。当叶片随径向位置增大而沿周向相位角发生变化时，称之为弯叶片，弯角β与风扇转子旋转方向相同时为正。

图12 弯掠OGV 示意图Fig.12 Sketch of leaned and swept OGV

基于三维升力面理论对风扇弯掠静子对噪声的影响进行快速解析计算［11］，结果如图13 所示。考虑结构强度的限制，选取侧弯角5°和后掠角23°作为设计方案，进行三维数值分析。需要指出的是基准构型静子存在15°的后掠角，但无侧弯角。

图13 OGV 弯掠角对声功率影响分析Fig.13 Impact of leaned and swept angle of OGV on acoustical power

采用与第3 节中一样的URANS/TPP 耦合数值计算方法，计算了103.7%转速工况的风扇转静干涉噪声。首先观察弯掠设计对气动性能的影响。如表2 所示，采用弯掠设计后，外涵流量和压比变化分别为0.09%和0.26%，未对气动性能造成不良影响。

表2 弯掠静子设计对气动性能影响Table 2 Impact of leaned and swept design on aerodynamic performance

风扇干涉噪声的强度取决于叶片对上洗速度的响应而产生的表面非定常载荷分布。图14 为静子叶片表面非定常压力幅值分布，可以观察到1BPF 幅值大于2BPF。同时，吸力面幅值大于压力面。吸力面幅值热区集中在静子前缘和叶身靠近前缘处，压力面幅值主要集中在静子前缘。

图14 弯掠静子表面脉动压力幅值分布Fig.14 Fluctuating pressure amplitude contours on surfaces of leaned and swept OGV

提取弯掠设计和基准构型的静子叶片吸力面非定常力实部，从图15 可以观察到，非定常力在叶片的径向和轴向呈周期分布，在不同单音频率下分布呈现不同形式。根据Envia 和Nallasamy［8］的理论分析，声场强度由叶片表面压力差决定，可以通过对非定常力在叶片表面积分获得。径向相位变化越大，对叶片不同区域积分时的相位相消越明显，从而可降低管道声场强度。通过数值结果分析可以观察到，采用弯掠设计后非定常力在径向上的周期性相比于基准构型更加明显。例如在2BPF 下，在叶片中下段前缘，弯掠构型的力分布周期数显著增大。

图15 基准构型和弯掠构型静子表面脉动压力实部分布Fig.15 Real part of fluctuating pressure contours on surfaces of leaned/swept OGV and baseline OGV

为了验证降噪效果，对静子后传噪声进行TPP 声源提取。图16 给出声模态分布，可以观察到弯掠设计对主要单音模态阶数未产生影响，但是降低了声模态幅值。对1BPF 下m=18 以及2BPF 下m=-19 和m=36 主要模态有明显的抑制效果。通过功率计算，根据表3 预测结果，本文数值方法和三维升力面模型存在一定差异，主要原因可能在于建模方法上的区别，但是在1BPF 和2BPF 下，不同方法的结果均表明声功率幅值得到抑制。上述结果表明，通过三维数值分析，实现了对低噪声弯掠静子设计的降噪基本规律的分析和降噪效果验证。