斜下降旋翼桨涡干扰声源定位及表面压力

2020-12-29刘正江汪文涛林永峰曹亚雄

航空学报 2020年12期

刘正江，汪文涛，林永峰，曹亚雄

中国直升机设计研究所重点实验室，景德镇 333000

旋翼是直升机特有的主要升力面和飞行姿态控制舵面，同时也是直升机最大的噪声源[1]。尤其在诸如小角度下降、拉平着陆等飞行状态，因离地面近，其噪声对周边环境必然产生更强烈的扰民影响。在这些近地飞行状态，旋翼桨-涡干扰噪声占据直升机噪声能量主要部分[2-3]。

针对桨-涡干扰(Blade-Vortex Interaction,BVI)噪声形成机理、声源位置以及传播特性，国内外一直在持续不断地开展着深入的研究。美国、德国和法国多家机构先后在1990年和2001年联合发起HART I和II计划，以B0105为试验对象，在DNW(German-Dutch Wind)消声风洞系统地开展了高阶谐波控制(Higher Harmonic Control, HHC)、单片桨叶控制(Individual Blade Control, IBC)的降噪试验研究[4-5]。航班直升机公司先进技术研究所为评估新设计的AT-2模型旋翼的基本参数，于2000年2月～3月在DNW开展了第2阶段模型旋翼试验，研究对比了基本型、1型和2型桨叶多种飞行状态下HHC对桨涡干扰噪声控制和振动控制的影响[6]。2002年4月，欧盟11个合作者共同启动了HeliNOVI项目，相关试验在德国航空航天研究院的DNW声学风洞进行，主要研究目标是减小振动和降低噪声。试验测量得到了矩形桨尖标准桨叶的瞬时声压、桨叶表面压力、动力学以及其他性能数据，研究了桨叶表面压力和噪声之间的关系，同时获得了主/尾桨-涡干扰过程中的桨尖涡轨迹和桨叶的相互位置关系的定量信息，建立了一个载荷和噪声综合数据库[7]。国外开展的上述项目不仅研究了桨涡干扰噪声，还同时研究了相同试验状态下的桨叶载荷/表面压力、流场以及桨叶形变。

中国如中国空气动力研究与发展中心[8]、南京航空航天大学[9-11]以及中国直升机设计研究所[12-13]开展过旋翼噪声试验研究，但大多进行的是单独噪声试验，同步开展相同状态下的斜下降旋翼桨-涡干扰声源定位试验和桨叶表面压力相关性试验较少。

本文首先对桨-涡干扰状态的桨叶表面载荷变化进行数值计算，对旋翼噪声数据去噪方法、桨-涡干扰识别、分离以及声源定位方法进行了研究，在直升机风洞同步开展模型旋翼斜下降桨-涡干扰声源定位和桨叶表面压力测量试验，研究桨-涡干扰状态、桨叶表面压力载荷变化特性以及桨-涡干扰噪声源特性，从时间和空间维度上分析桨-涡干扰噪声和桨叶表面压力载荷之间的相关性。

1 桨涡干扰状态桨叶表面载荷数值计算

1.1 自由尾迹计算桨叶表面载荷方法

采用自由尾迹方法进行斜下降桨涡干扰状态的桨叶表面气动载荷计算[14]。自由尾迹方法将旋翼尾迹离散为涡元，而涡元可以随当地气流运动，能很好地描述旋翼运动尾迹。计算中将尾迹分为近尾迹和远尾迹分别处理，由近尾迹过渡至远尾迹的涡龄角为60°。近尾迹区域内桨尖涡处于卷起过程中，尾迹保持涡面的形式；远尾迹桨尖涡完全形成，以线涡单元的形式对桨尖涡进行处理，并采用涡核模型消除桨尖涡附近计算奇异性。计算时先将桨叶离散为有限的气动分析单元，再采用Biot-Savart公式计算桨叶在不同展向位置和不同方位角的诱导速度，进而实现不同剖面桨叶的非定常载荷计算。Biot-Savart公式为

(1)

式中：dv为微段诱导速度；Γ为微段涡量；dσ为桨尖涡涡线微段；r为微段到周围流场任意点的距离的矢量；r为微段到周围流场任意点的距离。

1.2 自由尾迹计算结果

风洞试验采用4 m直径量级测压模型桨叶，桨毂为四支臂铰接式，桨叶半径R为2 100 mm，桨叶弦长D为140 mm，桨尖形状为矩形，采用OA213翼型。计算时选取桨叶表面载荷一个典型观测剖面位置：桨叶径向0.85R处。同时，旋翼转速为749 r/min，前进比为0.191(对应风速为30 m/s)。数值计算的桨叶和状态与风洞试验保持一致，桨叶表面气动载荷计算状态如表1所示。桨叶0.85R剖面的气动载荷计算结果如图1所示。

表1 计算参数Table 1 Computation parameters

图1 计算桨涡干扰状态桨叶表面载荷Fig.1 Computed rotor surface load in blade-vortex interaction

从图1可以看出，桨叶表面载荷在320°方位角附近(图1中虚线圆)振荡剧烈，表明桨叶在此方位角(后行侧)附近发生了较强烈的桨-涡干扰现象。桨叶表面载荷在后行侧320°方位角附近会产生桨-涡干扰，这与国外采用相近似参数试验件及状态得出的结果基本一致[5]。

2 桨-涡干扰噪声数据处理方法

2.1 改进型整周期时域同步平均去噪方法

旋翼旋转产生的气动噪声具有典型的以转速为基频的周期性特征，因此数据处理时首先对采集到的旋翼噪声信号进行整周期平均以剔除非旋转噪声干扰成分。整周期平均的主要思路为：利用同步采集的0方位脉冲键相信号作为表面压力和桨-涡干扰噪声数据块的起始点和结束点，将连续采集到的数据分割成整周期信号，然后将截取的整周期信号进行同序点平均。该平均方法充分利用了信号的相关性，不仅可以有效抑制随机性干扰噪声并提高旋翼噪声的信噪比，同时能较好地保持旋翼噪声信号幅值和相位特征信息[15-17]。

设测量信号由信号成分S(t)(确定性分量)和噪声干扰n(t)组成，则进行N次测量后，在某特定时刻ti经同步平均后的信噪比为

(2)

式中：Sj和nj分别为第j次采集时候的信号成分和干扰噪声信号成分。

试验时测量采用定频采样，采样率为每秒32K个采样点。高采样时转速对应到单个周期后会有轻微波动，体现在每圈数据长度会稍有差别，如果不进行相应处理而是直接进行整周期同序点平均，会造成相位错乱并最终导致后续声源定位出现较大误判。针对这种情况，对算法进行了如下改进：先按照整周期时域同步平均的方法将采集的若干秒数据截取成多个整圈数据，再统计出每圈数据的长度，得到数据长度多数值，最后从中挑选出长度为多数值的若干圈数据(一般不少于50 r，否则需增加采集数据时间)进行整周期同序点的同步平均。

2.2 基于小波变换桨涡干扰信号识别和分离方法

桨-涡干扰状态下的桨叶表面压力和桨-涡干扰噪声信号具有典型的脉冲特性，采用常规的滤波方法无法保证时域波形和相位不失真，而失真必将导致声源定位计算结果严重失实。基于小波变换的小波包处理不仅具有良好的时频局部细化能力，而且不会造成原有信号特征丢失。小波包分解是将频带进行多层次划分，再通过多分辨率分析方法对没有细分的频段进行进一步分解，从而提高时频分辨率。而根据熵理论及意义，不同类型的旋翼噪声在各个频段的能量是不一样的，熵值随之也不同，熵值对于信号的时变比较敏感，由其定义的小波包特征熵可表征信号复杂度在时域的变化情况以及丰富的频域特征[18-20]，因此选用其对桨-涡干扰噪声进行识别和分离处理。

利用小波包进行桨-涡干扰信号识别的具体步骤如下：

第1步信号分解

在LabVIEW环境下利用小波包工具箱提供的函数对整周期平均去噪后的桨-涡干扰信号进行n层小波包分解，分别提取第m(m≤n)层从低频到高频的k个频带信号；

第2步信号重构

对第1步之后得到的若干频带的小波包系数xj,k进行重构，得到k个小波包重构信号Sm,k(t)，其中各个重构信号分别包含了信号从低频到高频的各个频带的信息，信息量既无冗余，也不疏漏。对于m层的小波包分解而言，它是按照从最低频到最高频等频率间隔将信号分解到k个频带上。这里假设k=8，m=3，n=10。

第3步计算各频带信号的能量

设S3,k(t)(k=0,1,…,7)对应能量E3,k(k=0,1,…,7)，则有：

(3)

式中：xj,k(j=1,2,…,7;k=0,1,…,n)表示重构信号S3,k(t)离散点的幅值。对各个节点能量做处理组成向量，设节点向量为

E=[E3,0,E3,1,E3,2,E3,3,E3,4,E3,5,E3,6,E3,7]

(4)

再将各个节点能量做处理组成新的向量：

C=[C3,0,C3,1,C3,2,C3,3,C3,4,C3,5,C3,6,C3,7]

(5)

第4步构造信号小波包特征熵向量

(6)

式中：H3,k表示经过3层小波包分解后在8个频带的小波包特征熵。分别求得8个小波包特征熵，进而以8个小波包特征熵为元素可以构成一个特征向量表示为

H=[H3,0,H3,1,H3,2,H3,3,H3,4,

H3,5,H3,6,H3,7]

(7)

2.3 波束成形声源定位方法及软件实现

常用的声源定位方法主要有声强法、近场声全息(Nearfield Acoustical Holography，NAH)和远场波束形成(BeamForming)等，其中波束形成声源定位方法因测点需求少、计算速度快、中高频分辨率好、测量距离远、覆盖面积大以及易于布置等特点获得广泛工程应用[21-22]。

波束形成算法通常分两步来计算得到声源位置：首先计算不同麦克风对之间的时延；然后根据计算得到的时延值以及麦克风阵列的几何形状，采用几何法或插值法等方法来确定声源具体的空间位置。

2.3.1 bartlett时延算法

选用bartlett法计算时延，bartlett法是一种基于平均周期图法的功率谱估计方法，基本原理是先将N点有限长度序列分段后求周期图并进行平均，计算平均数据段的功率谱，再将功率谱进行多项式拟合得到拟合后的多项式系数矩阵，最后计算矩阵特征值即获得最后所需的时延值[23]。具体步骤如下：

首先，记输入信号数据为Y(n)，取q段平均，每段数据长度为L，第i段数据记为

Y(i)(n)=Y(n+(i-1)L)

(8)

式中：i=1,2,…,q;n=1,2,…,L。

第i段的周期图为

(9)

则平均后的功率谱为

P(k)=10lgPxx(k)

(10)

式中：Pxx为自功率谱函数。

然后对其进行多项式拟合，即：

(11)

式中：q0为多项式的项数；αq为常数项。则功率谱特征可表示为A=[α0,α1,…,αq0]。

再利用式(12)计算得到矩阵T：

T=ATA

(12)

最后计算矩阵T的特征值即可得到时延di0。

2.3.2 基于球面插值法的声源定位方法

在计算得到阵列传声器时延di0的基础上，后续确定声源位置的方法大体上分两种，一种是二维空间搜索方式，另一种是几何计算方法。二维空间搜索方法虽然可以精确定位空间声源位置，但当所用阵列传声器较多时，其搜索非常消耗计算时间和资源，因此在工程应用上受计算条件限制一般只能用于事后处理。而为提高计算效率以保证试验时声源定位结果展示具有一定的实时性，对比分析后选用了几何方法中综合效能最优的球面插值法进行声源定位计算[24-26]。

图2为球面插值法的空间定义。根据矢量几何和三角形三边关系，有：

新泰地处鲁中腹地，全市有林地面积达到117.1万亩，森林覆盖率达到39.8%，是山东省森林资源比较丰富的县之一。近年来，新泰市以明晰所有权、稳定承包权、放活经营权为核心，吸引社会资本进山入林，培育林业新型经营主体，打造林地集体所有、家庭承包、多元经营新格局。新泰被评为国家园林城市、全国绿化模范市、全国经济林产业区域品牌建设试点单位，“5+5”造林和经营模式在全省推广。

(13)

(14)

式中：RS为声源到传感器0的距离；Ri为传感器i到传感器0的距离；ri为声传感器i的空间向量ri=[xi,yi,zi]；rS为声源的空间向量。

由于di0是由bartlett时延计算得到的，与真实值有一定的误差，为此引入误差函数：

(15)

将式(15)得到的N个方程写成误差矩阵为

ε=δ-2RSd-2SrS

(16)

式中：

R*RT-d*dT

(17)

d=[d10,d20,…,dN0]T

(18)

(19)

图2 球面插值法空间定义Fig.2 Space definition of spherical interpolation

其中：R为由Ri组成的N× 1矩阵；d成为由di0为组成的N× 1矩阵。

式(16)可表示成：

(20)

式中：

(21)

采用最小平方误差估计，即使误差的平方和最小，表示为

(22)

为了使eLS达到最小，需要

(23)

写成矩阵分块形式为

(24)

简化为

(25)

式中：

(26)

式中：I为单位矩阵。

定义投影矩阵：

(27)

则声源位置的估计值可以化简为

rS,SI=(STPdS)-1STPdb

(28)

2.3.3 桨涡干扰噪声分析软件实现

软件采用LabVIEW、高级信号处理工具包和声音振动工具包开发，支持的前端采集硬件包括德国Bustec公司的3424系列VXI总线动态采集卡和美国NI公司的4472系列PXI/PCI总线动态采集卡。系统运行后可以自动识别在线硬件并生成通道参数配置表，用户可通过该表对传感器系数、三维空间坐标等参数进行配置(见图3)。主界面(见图4)中可以选择时间历程、幅值谱、倍频程以及声源定位云图等多种显示方式，用户可修改角度范围、比例系数等声源定位计算参数，声源定位云图中给出最大声源坐标位置、最大声源频率和总声压级等信息。

图3 通道配置界面Fig.3 Channel configuration interface

图4 声源定位软件主界面Fig.4 Noise source location software main interface

3 桨叶表面压力传感器和声阵列布置

3.1 桨叶表面压力传感器布置

试验前，在4号桨叶上表面弦向0.1D，径向0.90R、0.85R和0.80R3个位置分别布置KULITE的LQ-062微型表面压力传感器(见图5和图6)。桨叶加工时按照压力传感器大小预留圆形凹槽，压力传感器采用CN(Cyanoacrylate)胶固化在圆形凹槽内，压力传感器线缆走线采用和桨叶表面粘贴应变花相同的走线工艺，压力传感器的恒温电阻以上下沟槽压紧方式固定在桨毂中心处的轻质铝盒中。

图5 表面压力测点布置Fig.5 Arrangement of surface pressure measuring points

图6 表面压力测点实物图Fig.6 Photograph of surface pressure measuring points

3.2 麦克风声阵列选择及布置

常用的麦克风阵列包括十字阵列、矩形阵列、同心圆阵列以及螺旋阵列等，相对来说螺旋阵列具有最优的空间分辨率(即定位精度)，因此选择螺旋阵列，阵列有效直径为3.0 m。

受到试验场地空间位置以及阵列直径限制，同时采用的是远场声源定位算法，综合考虑上述因素后决定将声阵列布置在风洞流场外桨盘平面正下方偏前位置处，参考相关资料[9]阵列放置在该位置只能进行斜下降状态后行侧的桨-涡干扰噪声声源定位。为尽量降低地面声反射影响，同时又便于试验时的阵列和传声的安装以及现场标定位置调整，专门设计了一个高度为1.8 m的可移动阵列安装支架(见图7)。该安装支架可以通过滑轮推动，支架上部声阵列安装平面可以绕轴旋转已调节声阵列的俯仰角。试验前根据设定的试验状态，取试验时斜下降状态主轴倾角(此时，相当于桨盘平面后倒，见图8)中间值5.5°，将声阵列前移，阵列平面前倒到5.5°，并使旋翼主轴和阵列平面中心垂直轴在一条直线上，从而保证阵列平面和试验时的桨盘平面尽量(不含周期变距)平行。试验前安装声阵列支架以及地面包裹、铺设吸声海绵。

图7 声阵列实物图Fig.7 Photograph of microphone array

图8 声阵列风洞布置示意图(左视)Fig.8 Schematic diagram of microphone array in wind tunnel (left view)

声阵列布置空间位置示意图见图8。5 500 mm为桨毂中心到驻室进气段最近点水平距离，实线声阵列摆放的位置为桨盘平面正下方，虚点表示现场标定时其位置朝进气段前移998 mm，即移至桨盘平面正下方偏前位置。

4 声阵列图像畸变修正及现场校准

4.1 图像畸变修正算法

由于3 m直径螺旋声阵列不是完全刚性平面，无法精确保证所有传声器都在一个平面内，即螺旋声阵列平面和桨盘平面实际上存在一个很小的夹角，致使声源定位图像有一定的畸变失真，因此需要对其进行声源定位图像畸变矫正。

假设P(u,v)和Pd(ud,vd)分别是失真和无失真的图像坐标，两者之间的坐标位置关系如下：

(29)

式中：(u0,v0)为坐标原点(声源定位云图中心点)；Ki为径向畸变系数；rp为图像点到圆心的距离。实际使用只需考虑前两个畸变项K1和K2。

4.2 声阵列现场校准及精度分析

试验前，将旋翼桨盘平面后倒5.5°，通过在0°、90°和270°桨尖位置附近的定频全向喇叭(1 000 Hz 单频声)对声阵列参数进行了校准和图像畸变修正，最终结果见图9。

图9 声阵列参数校正后声源定位结果云图Fig.9 Cloud image of sound source localization results after acoustic array parameter correction

从表2可以看出，经过校准和图像畸变修正后，声源定位最大方位差为2.883°、最大径向位置差为0.010R，满足试验对上述值误差要求。

表2 现场校准声源定位误差分析Table 2 Error analysis of noise source location calibration

5 典型试验数据分析

5.1 桨涡干扰状态BVI噪声数据分析

选取一组典型试验状态数据进行分析，试验状态为旋翼转速749 r/min、风速30 m/s、主轴倾角5.5°，取状态稳定后的10 s数据，按2.1和2.2节介绍的数据处理方法，得到平均后的一圈数据，再进行BVI噪声的小波识别和分离，处理前后的结果见图10。

对比图10(a)和图10(b)可以得到，采用2.1和2.2节介绍的方法能够清晰识别并提取出桨涡干扰噪声，同时能够较好地保留桨涡干扰噪声的幅值和相位波形特征。

图10 声阵列噪声整圈声压历程Fig.10 One rotating period noise data of microphone array

5.2 桨涡干扰状态桨叶表面压力数据分析

选用和5.1节相同试验状态数据以及数据处理方法，对旋翼转速749 r/min、风速30 m/s的不同主轴倾角试验状态(见表3)的桨叶表面压力数据进行平均及识别和提取处理。同时，已知粘贴微型表面压力传感器的4号桨叶方位角为266°，转速0脉冲传感器的方位角为356°，据此对提取出的表面压力数据进行方位角修正。图11给出典型试验状态下0.90R、0.85R和0.80R剖面的一整圈(0°～360°方位)的时间历程曲线，轨迹角为6.5°。

从图11以及表3可以看出：

1) 桨涡干扰引起的桨叶表面压力变化主要集中在320°方位和60°～70°方位附近，前者为后行侧桨涡干扰引起的桨叶表面压力变化，后者为前行侧桨涡干扰引起的变化。