基于致动线/BPM模型结合的风力机气动噪声研究

2019-08-27杨建刚

噪声与振动控制 2019年4期

程志，杨建刚

（东南大学火电机组振动国家工程研究中心，南京 210096）

进入21世纪以来，全球风力机的发电功率在不断増加，风资源开发已从沿海向内陆展开[1]。风资源开发距离居民生活区越来越近。同时，为了提高风力机的效率，风力机也一直在向大型化发展，其气动噪声也不断增加，不可避免影响到居民的生活。

Nii和Yoshinor等[2]通过研究发现风电场在规划和选址时，气动噪声是必须要考虑的影响因子之一。荷兰、德国和丹麦等国家对风力机产生的噪声水平是否影响周边居民和动物生活进行过测量，发现风力机气动噪声对候鸟迁徙有很大影响[3]。同时音量和频率等声音参数的变化将会直接影响居民的身心健康。因此，风力机气动噪声也成为近来的研究热点。

气动声学的研究主要分为直接模拟[4-6]和混合模拟[7-9]两类。直接模拟方法虽然精度高，但其计算量巨大。混合模拟方法将流场与声场的计算分开，更适合于复杂的工业模型，其中各种软件中流场的计算方法已经相对成熟，但声场模拟仍然相对滞后。除了经验方法外，目前被广泛使用的声场计算方法是Lighthill声理论[7]以及其延伸理论，可对翼型噪声[8]、风力机近远噪声开展研究[9-12]。但该方法是基于 Lighthill、Lighthill-Curle 和 Ffowcs Williams-Hawkings方程，经过一定简化后推导得到的解析解[13]，会忽略掉流场中一些压强波动，造成误差[14]，具有局限性。相比之下，利用原始声学方程[15-16]进行微分计算更准确且适用面更广泛，并且能够展现出声波向外传播的过程[17]。目前国内的声场研究仍然广泛基于常用的积分法开展。本文声场计算基于声学扰动方程微分法展开。

在风力发电机的流场仿真计算中，致动线模型因为不用对风力机实体模型进行解析，免去了在叶片表面边界层进行高精度解析的要求，对网格要求较低、计算速度快[18]，在风场远场以及尾流计算中有很大优势并被广泛使用[19]。但致动线模型无法解析叶片边界层，边界层的翼型自噪声声源无法通过致动线方法获得。BPM模型不包括尾流波动带来的噪声源，但可以提供在致动线模型中缺失的翼型自噪声声源，同时基于声学波动方程计算声场具有精度高、能够展现声波传播过程等优点。因此本文将致动线方法和BPM模型进行结合来计算声学扰动方程的声源项。

本文以OpenFOAM为平台，模拟点声源扩散，并通过与声学理论的对比，验证声学扰动方程以及代码编写的正确性；然后，基于致动线模型计算并分析某型号风力机的流场，同时基于BPM模型计算同一型号风力机的声压情况；最后，将根据致动线模型计算出的流场与根据BPM模型计算得出的声源结合，作为完整声源项放入声学波动方程来计算声场，对声场的计算结果进行分析。

1 数值方法

1.1 不可压缩混合RANS/APE声场计算

本文声场计算基于混合方法来开展，在该方法中，详细的流体参数被分解成不可压缩的流动变量及可压缩的波动变量

式中：u′、p′和ρ′为瞬时不可压缩速度、瞬时不可压缩压强和瞬时不可压缩密度，U、ρ0和P分别为流场计算中得到的不可压缩速度、不可压缩密度和不可压缩压强。其中不可压缩流动的参数包括了湍流流场的信息，后面的可压缩波动参数代表声场情况，本文风力机流场致动线计算基于RANS湍流模型。

当获得不可压缩流场信息后，便可以基于声学扰动方程进行声场的计算。Ewert在论文中提及了多种不同形式的声波扰动方程[20]，本文使用的声波扰动方程在某一维方向上表达如下

式中：i、j、k分别表示三维向量方向，APE方程左边的波动参数包括了非稳态流动中声场传播的信息，方程的右边涵盖了声源项，这里是来自流场的压强变化为中间变量，c为声速，δ为张量转换因ij子为黏性力。关于声学波动信息的传播方程是从原始的可压缩流动方程推导而来，需要从可压缩流动方程中将不可压缩部分除去[21]。由于空气中的压强扰动引起耳膜振动并经听觉神经传至大脑被确定为声音[13，22]，得到流体波动参数后便可以得到声音信息。APE声学扰动方程已经在低马赫数流动问题、湍流问题以及点声源扩散问题中被验证正确[23-24]。

1.2 致动线模型

致动线模型是一种计算风力发电机流场的常用方法。风力机的叶片旋转时会施加气动力给周围的空气，基于BEM叶素模型将叶片和周围流体相互作用的体积力提取出来，然后施加到致动线模型的三维纳维-斯托克斯控制方程中，用来替代叶片表面的负载。

式中：f为施加的体积力。

通过计算BEM叶素模型以及翼型参数，可以得到施加的体积力f[19]。BEM叶素模型将叶片沿着翼展划分成连续的小段，然后将每一段叶素视为一个二维翼型。每段叶素的攻角由该段所在的速度三角形决定，每段叶素的升力和阻力由该段所在升阻力系数决定。该模型运用广泛，数学原理在此不做阐述。

由于致动线模型无法解析叶片边界层，边界层的翼型自噪声声源无法通过致动线方法获得。所以引入BPM模型来考虑翼型自噪声声源的影响。

1.3 制动线模型-BPM模型结合

对于翼型自噪声，BPM模型是最常用的半经验模型，由Brooks、Pope和Marcolini在实验测试的基础上提出[25]。翼型自噪声又分为：湍流边界层尾缘噪声、分离流噪声、层流边界层涡脱噪声、叶尖涡噪声和钝尾缘涡脱噪声。每种噪声类型都有自己的半经验公式，在每段叶素上，通过半经验公式提取得到10 Hz到20000 Hz中三分之一声压频段的5种噪声来源。来自于致动线模型的噪声源仍然保留在式（5）中。这里通过高斯正则化对来自于BPM半经验模型的噪声源进行处理，同时为了避免数据计算产生奇异性，不能用声源点声压替代，研究发现距离各段叶素一个弦长位置范围内的声压已经保持不变，本文用距离各段叶素十分之一弦长位置的声压来替代该叶素所在位置的声压，然后转换为波动压强

然后式（5）被改写为

式中：m为风力机叶片的数目，n为每段叶片被划分的叶素数量，l为10 Hz到20000 Hz之间三分之一倍频程频率的数量，根据三分之一倍频程的划分规则，文中l为31。d(=|x-se|)仍为正则化过程中某坐标为(x,y,z)的网格中心点到对应致动线中心点的距离，ϵ为BPM声源波动压强的高斯光滑系数。

2 算例与结果分析

2.1 点声源扩散

先对基于声学扰动方程的声场计算进行验证分析，将计算结果与空气中点声源的耗散理论进行对比验证。计算域模拟边长为200 m的的正方形扩散域，在正方形扩散域中心设置点声源。时间格式采用Crank-Nicholson格式，其它如梯度格式、散度格式及拉普拉斯项格式等采用Gauss linear 2阶格式。网格设置为4001×5×5，可以看出，Y方向和Z方向网格十分稀疏，在这两个方向传播的声参数由于无法被网格识别会直接耗散。对声音在X方向上的传播情况进行分析，采集点在Y及Z方向上与声源点的距离为0，在X方向上与声源点的距离分别为2、4、6、8、10、20、30、40、50，单位是米，如图1所示。

图1 点声源示意图

将点声源频率设置为80 Hz，幅值一定，当计算时间达0.5 s时，声场传播到计算域边界。t=0.4 s时的声压传播截面图如图2所示，可以看出Y及Z方向的声波被粗糙网格快速耗散。点声源向四周传播时波动压强的数值越来越小。

图280 Hz点声源扩散t=0.4 s时压强波动截面图

根据点声源在干燥空气中传播的斯托克斯声衰减定律[22]，声压的衰减满足公式

式中：A0可以为某点的声压，A()d为与该点距离为d的另外一点的声压，α与传播介质、声源形式以及传播方式等有关，与声源频率的平方成正比，满足α∝ω2。

将点声源频率分别设置为40 Hz、80 Hz、200 Hz及300 Hz，幅值相同，为一定值，其他模拟参数相同。表1是不同频率点声源传播模拟中相邻采集点间衰减指数α的变化情况。

图3（a）是不同频率点声源传播计算结果及其拟合曲线。结合表1及图3（a）可以看出，低频点声源拟合情况优秀，高频点声源由于传播衰减太快，远距离时网格难以捕捉声参数，模拟得出的高频衰减指数在远距离时有误差。在另外的网格无关性验证中发现，在工程噪声问题分析中，由于高频声音衰减很快，进行远场分析时可以不考虑，可以根据情况调整网格精度。

表1 不同频率声源相邻采集点间衰减指数α的变化情况

图3（b）展示了不同频率噪声衰减指数与频率之间的关系，说明在某一网格精度下，噪声衰减模拟满足α∝ω2。可以看出，用文中基于APE编写的OpenFOAM声场求解器模拟声源传播时，符合衰减指数α与频率平方成正比的斯托克斯声衰减定律，基于声学扰动方程的声场计算结果正确。

2.2 NREL5MWRef风力机流场

验证了声场计算方法的正确性后，这部分以及下节将研究分析NREL5MWRef型号风力机的流场和声场。NREL5MWRef型号风力机以及计算风场域的参数如表2所示。风场计算域如图4所示。

表2 NREL5MWRef风力机参数

图3 不同频率点声源传播中采集点波动压强

图4 风场计算域示意图

由于致动线模型和BPM模型的植入都不需要对叶片边界层局部解析，并且下文基于声学扰动方程的声场计算需要在远场保持相同的网格精度，结合研究经验，平衡计算精度和时间后，沿整个计算域3个方向采用等距的200×200×60网格布置。流场计算中时间步长δt为0.05 s，时间微分采用Crank-Nicholson格式，其它离散如梯度格式、散度格式及拉普拉斯项格式等采用Gauss linear2阶格式。在时间步t=80 s时，计算域内的非稳态流场已经完全被解析。此时流场的速度云图以及湍动能云图如图5所示。

速度云图分为3个扇区，对应3个叶片，反映了风力机叶片在旋转过程中的速度变化规律。由图5（a）可以看出，风力机叶片旋转区域具有后面有明显的速度下降、随着距离越来越远速度逐渐恢复的趋势；本文设定风力机来流为均匀来流，故如图5（b）所示，整体湍流强度低。但在风力机后面的位置仍然有明显的湍流存在，并随着距离越来越远，湍流强度减轻。湍流的云图趋势和速度的云图趋势保持一致，并与实际情况相吻合。

图6为风力机尾流的涡量等值面，能够清晰反映尾迹流场的特点。由于NREL5MWRef型号风力机尺寸大，叶尖速比高，叶尖涡会迅速破碎耗散，且由于此时泄涡频率很大，尾涡的螺旋状结构相对小风力机不太明显。

图5 风力机流场

图6 风力机涡量图

该型号风力机暂时没有对应的实验数据，但致动线模型是一种很成熟的风力机流场计算模型。其正确性在以往大量研究中都有验证，本课题组关于致动线模型流场也有大量的基础工作。

2.3 NREL5MWRef风力机声场

在前文得到的风力机流场的基础上，加上从BPM模型中得到的声源，进行声场计算。时间微分采用Crank-Nicholson格式，其它离散如梯度格式、散度格式及拉普拉斯项格式等采用Gauss linear 2阶格式。受计算量与计算时间的限制，时间步长δt取0.002 s，这对高频噪声的解析有一定影响，不过实际情况中高频噪声消减很快，在远场分析时可以放在次要位置。该声场计算以完全解析流场为起始点，采用并行计算，调用超算共75个2 G的CPU，计算时间3天。

人耳能够识别的极限声源频率范围为10 Hz到20000 Hz，故在从BPM模型中抽离声源时，选择10Hz到20000 Hz范围内的三分之一倍频程。图7展示了从BPM模型中抽离出的20 Hz、200 Hz、1000 Hz和20000 Hz频率声源，由于叶片的偏转，某一垂直截面上3个叶片的压强云图不一样。

由图7中可以看出，不同频率情况下声源的量级有很大差距，中间频率部分的声压级水平更高；另外，叶片上声源重心位置不同，这体现出叶片位置、声源频率以及发声机理之间的关系。虽然不同频率下声源位置分布不一，但整体上呈现频率越高、声源重心越靠近叶尖位置的现象，这是因为叶尖位置与空气的相对速度大，涡脱落频率高，流场变化相对叶根更为剧烈，声源频率高。

图7 不同频率的BPM声源声压云图

图8是波动压强在不同时刻的云图，体现了声压逐渐向外传播的过程，最后再覆盖到全域。声音在空气中扩散时，一方面随着传播球面的表面积越来越大，单位面积声压强度越来越小，同时由于空气介质的吸收作用，携带声音信息的粒子在介质中震荡能量被消耗，高频噪声震荡越快，消耗速度越快，所以高频噪声消散很快，从图8中可以看出随着声压的扩散，强度逐渐减轻。同时在靠近地面附近，由于声压的反射可以看到干涉情况出现，各截图的底部为地面。

图9为不同观察点（-200，0，0）及（-100，0，0）处波动压强的频谱图，点（-200，0，0）在风力机正上方。由于风力机的转速是12.1 r/min，可以看到，各位置声压级的频谱在该频率附近能量很高，交织在一起，并且由于低频衰减较慢，低频声压级没有什么差别。随着频率增大，可以看到点（-200，0，0）仍然呈现一定的波峰，这些波峰有相当一部分来自于从BPM半经验模型中添加的噪声源，但较远处声压级的高频部分迅速衰减，能量变小，这与声源传播的规律相一致。很多调查显示，风场附近的居民对高频噪声的感知较弱，但风力机低频噪声产生的连续不断的“嗡嗡嗡”声音会使他们极为难受，尤其在夜间使人难以入眠。

图8 不同时刻声压波动p′

图9 声压级频谱图

图10为声场被完全解析后，在风力机轴心处3个垂直法向截面上SPL声压级的云图，可以看出沿着下风向的声压级强度比上风向的声压级强度更大，风力机周围的云图不完全对称，向下风向轻微偏移。

图10 风力机周围SPL声压级云图

风力机周围的声压级可以达到100 dB，但在向外扩散过程中衰减，风力机垂直于风向两侧的声压级云图较弱，这与基于各种半经验模型[26]画出的声压级指向性图特征相似。由于网格精度的限制，远场的真实声压级会更强。致动线模型没有考虑塔架的影响，在后面的研究中可以基于致动线相似的原理，将塔架简化成固定的体积力加入控制方程，模拟塔架对声场以及流场的作用，以提高精度。

图11给出了风力机周围半径100 m和半径240 m处的SPL指向性分析。

图11 不同半径处声压指向性图

可以看出，由于采集点离风力机很近，噪声还未开始进行对称式扩散衰减，半径100 m处的声压指向性图分布相对均匀，这一点与图10（c）一致；而到了半径240 m处，声场在扩散的过程中逐渐变得不均匀，但呈现对称分布的特征，且在风力机两侧声压级收缩，在下风向位置声压级强度突出。240 m处声压级强度整体小于100 m处的声压级强度，但在局部由于地面反射作用声压强度更强，地面反射的作用在远场会被弱化。

由于该风力机体型很大，难以开展相关的流场以及声场实验，缺乏实验数据对照。但本文对该风力机进行仿真计算时，除了成熟的制动线模型外，采用的声场计算方法在点声源扩散模拟中已经得到验证。

在以后的风场设计、城市小风力机等设计时，需要综合考虑距离、方位的影响，结合噪声类型、频率对人的听觉干扰，根据噪声传播特征布置合适的隔音板，建立绿色、清洁、舒适的生活工作环境。