赵 华，周鹏展，张 亢，廖力达

(长沙理工大学能源与动力工程学院，湖南长沙 410114)

为了保障国家能源安全、保护生态环境、促进社会和经济的可持续发展，近年来中国风电产业迅猛发展，到2013年底，风电已经成为中国第三大能源，风电新增装机容量为16 088.7 MW，同比增长24.1%;累计装机容量为91 412.9 MW，同比增长21.4%，两项数据均居世界第一［1］。中国风资源丰富，主要集中在三北(东北、华北、西北)以及东南沿海及附近岛屿，而电力需求主要集中在东南沿海及中部人口密集地区，受电网“消纳”能力的限制，目前中国风电应用正在向沿海和中部风资源相对丰富、经济比较发达的地区扩展，但由于沿海和中部地区人口比较稠密，风力机在运行时产生的气动噪声会对风电场附近居民的生活造成影响［2］，而中国在2008年颁布的《风电场噪声限值及测量方法》标准中规定，居民生活区的噪声不得超过60 dB，因此风力机叶片的气动噪声问题开始引起研究人员的关注。

现阶段对风力机叶片气动噪声的分析主要有两类方法，一类是应用CFD软件对其气动噪声进行计算［3－4］，且费时费力［5－6］。另一类是基于大量实验数据，拟合得到气动噪声分析模型，其通用性和拓展性不佳［7］，两者都难以适应风力机叶片气动噪声随风速变化的动态分析需求［8－9］。

本文基于传统气动声学理论［10］，考虑风力机运行状态、来流风速和接收点位置等因素的影响，建立风力机叶片气动噪声的修正模型，运用时域分析方法对一种2 MW的风力机叶片气动噪声进行数值计算，并通过现场实测数据进行实验验证，为开发低噪声风力机叶片提供理论依据。

1 气动噪声形成机理

风力机叶片气动噪声的产生是由于空气在风力机叶片周围的高速流动引起的［11］，其流动如图1所示，其气动噪声的大小与风力机叶片的转速、入流速度、位置等因素有关，风力机叶片旋转过程中叶片厚度变化对叶片表面气流的影响和运动速度的变化，导致厚度噪声的产生;风力机叶片载荷的变化，导致载荷噪声的产生。

图1 风力机叶片气流流动示意图Fig.1 Schematic of the flow around the wind turbine blades

2 气动噪声时域分析方法

2.1 风力机叶片气动噪声修正模型

本文基于传统气动声学理论［12］，在关于面声源的时域解［13］的基础上，考虑风力机叶片的运行状态、来流风速以及接收点位置等因素的影响，将风力机叶片气动噪声分为载荷噪声［14－15］、厚度噪声［16－17］和扰动噪声3种进行分析，建立风力机叶片气动噪声的修正模型。

1)风力机叶片载荷噪声的计算公式

式中:l˙i为相对于时间 τ的载荷导数，下标i表示沿xi轴方向的分量为f=0上的当地速度;Mai=Vi/c0，M˙ai为相对于时间τ的Mai导数；=(xi－yi)/r为单位辐射矢径;ds为f=0面上的面积元素为接收点距离声源的距离;下角标ret表示方括号中所有对应于延迟时间的参量都取它们各自在其发射时刻τe的值;li为叶片表面对当地流体的法向作用力与xi轴方向上的分量;lr为l在接收点方向上的投影。

2)风力机叶片厚度噪声的计算公式

式中:V˙n=∂Vn/∂τ，Vn=－∂f/∂τ为叶片表面的法向速度。

3)风力机叶片扰动噪声的计算公式

根据风力机的运行状态，考虑来流风速对风力机叶片气动噪声计算结果的影响，设定扰动噪声谱级与风速的拟合关系为

式中:参数v为当地风速，单位m/s;参数p'为扰动噪声，单位dB;参数a，b为修正系数。

本文通过实际测量湖南风电场某2 MW风力机来流风速和气动噪声数据，根据最小方差准则，经过大量数据统计计算，确定取a=1.03，b=21.8，则风力机叶片扰动噪声的计算公式可写为

综上所述，本文基于传统气动声学理论建立的风力机叶片气动噪声修正模型为

2.2 风力机叶片气动载荷计算方法

风力机叶片上的气动力载荷主要分为摆振载荷Mb与挥舞载荷Mh［18－19］，基于动量叶素理论，把风力机叶片沿展向分成20个叶素，设第i个叶素所对应的摆振载荷和挥舞载荷分别为Mbi和Mhi，将各叶素上的对应载荷进行累加，得到整个风力机叶片的摆振载荷Mb与挥舞载荷Mh分别为

2.3 坐标变换

时域分析法预测风力机叶片气动噪声［13］的计算公式中，被积函数是延迟时间和空间位置的函数，需要用坐标转换确定被积函数中声源位置相对于固定坐标下接收点的位置。

固定坐标系O1A1B1C1固定在未扰动流体中，其原点O1位于风力机叶片的轮毂中心，A1轴与风力机塔筒平行，向下为正，C1与地面平行，逆风向为正，B1轴方向位于A1轴右侧，与A1轴和C1轴构成右手坐标系。

旋转坐标系O1A2B2C2的C2轴与固定坐标系中C1轴重合，A2轴和B2轴位于风力机叶片的旋转平面内，且B2轴为沿风力机叶片变桨轴线在旋转平面内的投影方向，θ为风力机叶片方位角。固定坐标系到旋转坐标系的转换关系为

式中:(i1，j1，k1)表示空间某点位于固定坐标系中的位置;(i2，j2，k2)表示空间某点位于旋转坐标系中的位置。

随叶坐标系O1A3B3C3的B3轴为沿风力机叶片变桨轴方向，与B2轴重合，ψ为桨距角，如图2所示。随叶坐标系到旋转坐标系的转换关系为

式中:(i3，j3，k3)表示空间某点位于随叶坐标系中的位置;(i2，j2，k2)表示空间某点位于旋转坐标系中的位置。

图2 旋转坐标系和随叶坐标系Fig.2 Rotating coordinate system and leaf coordinate system

2.4 延迟时间的求解

延迟时间方程为

延迟时间方程是用来确定叶片表面某点声源在离接收点不同距离处发射声波时刻的基本方程，延迟时间方程g=0表示一个半径为r=c0(t－τ)的球面。风力机叶片气动噪声的计算是在整个风力机叶片表面进行的，将风力机叶片的表面划分为若干个叶素，对于某一确定位置x和确定接收时刻t的接收点来说，风力机叶片表面每一叶素上的声源都对应一个不同的延迟时间τ，对面积进行积分，通过反复迭代计算得出延迟时间［4］。

将位于随叶坐标系中风力机叶片表面某点的坐标(a'，b'，c')转换为该点在固定坐标系下的轨迹y(τ)。设固定坐标系下接收点的位置矢量为x=(x1，x2，x3)，则其延迟时间方程为

3 数值仿真

根据风力机叶片气动噪声时域分析方法，对海拔高度在1 600～1 860 m之间的湖南风电场的某2 MW风力机叶片的气动噪声进行数值仿真［20－22］，计算所需的风力机基本参数如表1所示，其控制方式为变桨控制。

表1 风力机基本参数表Tab.1 Basic parameters ofwind turbine

运用风力机叶片气动噪声时域分析方法，根据风力机运行状态，计算得到风力机叶片气动载荷以及风力机叶片的速度分布，根据接收点的位置，运用坐标系变换关系，将固定坐标系下的风力机叶片速度分布转换为随叶坐标系下的速度分布，并根据扰动噪声计算公式，引入来流风速对气动噪声的影响，通过风力机叶片气动噪声的修正模型，计算得到风力机叶片气动噪声的理论值。风力机叶片气动噪声时域分析计算流程图如图3所示。

基于Matlab软件中的Simulink仿真模块，运用时域分析方法，建立如式(1)—式(4)所示的风力机叶片气动噪声修正模型，计算得到的风力机叶片气动噪声仿真图如图4所示。

图3 风力机叶片气动噪声时域分析流程图Fig.3 Flow chart of time domain analysis of aerodynamic noises from wind turbine blades

图4 风力机叶片气动噪声的时域计算仿真图Fig.4 Simulation diagram of the time domain calculation of aerodynamic noises from wind turbine blades

在对风力机叶片气动噪声进行计算时，将风力机叶片分为n个叶素，本文中取n=20，对各叶素相对于接收点距离10 m和各接收点时刻t=0 s所对应的各发射时刻τ，利用PPtime函数进行求解，在各叶素上声源的发射时刻确定后，由坐标系转换关系ZBZH函数求出风力机叶片各叶素在固定坐标系下的几何位置、运动状态，根据load函数计算得到风力机叶片气动载荷，确定风力机叶片载荷噪声和厚度噪声计算公式中的被积函数值，根据来流风速计算得到扰动噪声，利用建立的修正模型求出各叶素上各声源所对应的声场声压解。最后将其所有声源的声场解进行叠加，就构成了声压时间历程中的一个点，再考虑各个接收点的时刻t，求出以时间为横坐标，声压值为纵坐标的波形图。理论计算得到的风力机叶片气动噪声A加权时间历程图如图5所示。

图5为根据风力机叶片气动噪声的修正模型，运用时域分析方法得到90 s内风力机叶片气动噪声的时间历程图。风力机叶片在绕叶轮中心高速旋转时，在运行过程中使得气流发生剧烈的扰动，导致气动噪声的产生。风力机叶片在旋转平面的每个位置相对于距离塔筒中心10 m的正面接收点来说差别不大，正面接收点的声压值波动较小，其平均声压值约为57 dB(A);当风力机叶片旋转时，相对距离塔筒中心10 m的侧面接收点而言，其位置差别较大，这时侧面接收点的声压值波动也较大，但其平均声压值较小，大约为53 dB(A)。

图5 理论计算的风力机叶片气动噪声时间历程图Fig.5 Time sequence diagrams of theoretical calculations of wind turbine blades aerodynamic noise

4 实验结果分析

采用HT－8352噪声测量仪对该风力机叶片气动噪声进行测量，测量点的布置如图6所示，在距离风力机塔筒地面中心位置R为10 m的正面和侧面分别布置测试点MP1和MP2。测量时间间隔为1 s，并通过数据线导出到计算机，对数据进行后处理。

在风力机叶片气动噪声的测量过程中，通过风电场集控室获取90 s时间内风力机来流风速及对应风速下输出功率的具体数据，并绘制成曲线如图7所示。

图6 测量点布置示意图Fig.6 Schematic diagram ofmeasuring point

图7 风速与输出功率随时间的变化曲线图Fig.7 Variation diagrams of wind speed and output power with time

由于A加权声压级能较好地反映人对噪声的主观感觉，本文根据国家标准GB/T 2888—2008《风机和罗茨鼓风机噪声测量方法》给出的处理方法，绘制A加权声压曲线。在图7所示90 s内的风速下，对风力机叶片气动噪声进行测量，得到的风力机叶片气动噪声时间历程图如图8所示。

距离风力机塔筒地面中心10 m的正面接收点，在90 s内接收到的风力机叶片气动噪声的平均噪声水平为57 dB(A)，距离风力机塔筒地面中心10m的侧面接收点，在90 s内接收到的风力机叶片气动噪声的平均噪声水平为53 dB(A)，从图8可见，其正面接收点的气动噪声的平均值比侧面高，说明该风场的风力机气动噪声具有指向性。

理论计算结果与实验测试结果的比较如图9所示，实验测试是利用HT－8352噪声测量仪对风力机叶片气动噪声进行测量，并记录测量时间和接收点的位置，利用风电场集控室的监视结果获取测量时间内对应的风速和风力机输出功率，理论计算结果运用时域分析方法，根据建立风力机叶片气动噪声的修正模型，进行自编写程序计算得到实验测量时间对应风速下风力机叶片气动噪声值。从图9可见理论计算值与实验测试值的平均值比较吻合，说明基于传统气动声学理论，考虑风力机运行状态、来流风速以及接收点位置的影响，建立的风力机叶片气动噪声修正模型是正确的。但由于实际测试受外界环境的植被吸收、地面反射等因素的影响，以及测量设备、测量方法等原因使得实际测试结果与理论计算结果之间会存在一定的局部偏差。

图8 实测风力机叶片气动噪声时间历程图Fig.8 Time sequence diagrams ofmeasured wind turbine blades aerodynamic noise

图9 理论与实测风力机叶片气动噪声时间历程对比图Fig.9 Time sequence diagrams of theoretical calculations and measured wind turbine blades aerodynamic noise

5 结论

以湖南风电场某2 MW风力机叶片为研究对象，建立其气动噪声修正模型，运用时域分析方法对风力机叶片气动噪声进行数值计算，并通过与实验测试结果对比，得到了以下结论。

1)在传统气动声学理论基础上，考虑风力机运行状态、来流风速以及接收点位置等因素的影响，建立了风力机叶片气动噪声的修正模型。

2)根据最小方差准则，通过对湖南风电场某2 MW风力机的大量气动噪声实测数据的统计分析，确定其气动噪声修正模型中扰动噪声的计算公式为10 lg(v)=1.03*p'+21.8。

3)基于Matlab软件的Simulink模块，运用时域分析方法，通过自编程序计算得到的距离风力机塔筒地面中心10 m的正面接收点的气动噪声的平均声压值为57 dB(A)，其侧面接收点的气动噪声的平均声压值为53 dB(A)，说明该风力机叶片气动噪声具有指向性。

4)现场实验测试表明，测得的气动噪声平均值与计算结果吻合，证明本文建立的风力机叶片气动噪声修正模型是正确的。

/References:

［1］ 中国风能协会.2013 年中国风电装机容量统计［EB/OL］.http://www.cwea.org.cn/download/display_info.asp?id=55，2014－03－26.Chinese Wind Energy Association.China’sWind Power Installed Capacity in 2013［EB/OL］.http://www.cwea.org.cn/download/display_info.asp?id=55，2014－03－26.

［2］ 卢志才，米东，徐章遂，等.多尺度多元素形态滤波自适应降噪研究［J］.河北科技大学学报，2011，32(3):228－232.LU Zhicai，MIDong，XU Zhangsui，et al.Research in application ofmulti－scale andmulti－elementmorphological filtering in geomagnetic orientation［J］.Journal of Hebei University of Science and Technology，2011，32(3):228－232.

［3］ FARASSAT F.Derivation of Formulations Land 1A of Farassat［R］.Hampton:NASA Langley Research Center，2007.

［4］ 吉霞，吴亚锋，杨浩.管道泄漏噪声的数值仿真分析［J］.计算机仿真，2010，27(12):381－385.JIXia，WU Yafeng，YANG Hao.Numerical smiulation of pipeline leak noise analysis［J］.Computer Simulation，2010，27(12):381－385.

［5］ BAIY L，MA X Y，MING X.Lift enhancement of airfoil and tip flow control forwind turbine［J］.Applied Mathematics and Mechanics－English Edition，2011，32(7):825－836.

［6］ KIM H，LEE S，SON E，et al.Aerodynamic noise analysis of large horizontal axis wind turbines considering fluid－structure interaction［J］.Renewable Energy，2012，42:46－53.

［7］ GÖGMEN T，ZERDEM B.Airfoil optimization for noise emission problem and aerodynamic performance criterion on small scalewind turbines［J］.Energy，2012，46(1):62－71.

［8］ IANNIELLO S，MUSCARIR，MASCIO A D.Ship underwater noise assessmentby the acoustic analogy.Part I:Nonlinear analysis of amarine propeller in a uniform flow［J］.Journal of Marine Science and Technology，2013，18(4):547－570.

［9］ VAGLE S，WILLIAM G L，DAVID M F.An evaluation of the WOTAN technique of inferring oceanic winds from underwater ambient sound［J］.Journal of Atmospheric and Oceanic Technology，1990(4):576－595.

［10］ 张强.气动声学［M］.北京:国防工业出版社，2012.ZHANG Qiang.Aeroacoustics［M］.Beijing:National Defense Industry Press，2012.

［11］ YAO Huadong，DAVIDSON L，ERIKSSON L E，et al.Surface integral analogy approaches for predicting noise from 3D high－lift low－noise wings［J］.Acta Mechanica Sinica，2014，30(3):326－338.

［12］ 谢爱云，黄雅婷，冯涛，等.高频超声清洗中声场分布的仿真及实验研究［J］.河北科技大学学报，2014，35(4):324－329.XIE Aiyun，HUANG Yating，FENG Tao，et al.Simulation and experments research of acoustic field distribution in high frequency ultrasonic cleaning［J］.Journal of Hebei University of Science and Technology，2014，35(4):324－329.

［13］ AMITAIE，JEFFREY A.Nystuen in Precipitation:Advances in Measurement，Estimation and Prediction［M］.Editors:SilasMichaelides，2008.

［14］ CAO Tuanvu，WISLAND D T，LANDE T S，et al.Low phase－noise VCO utilizingmodified symmetric load and partial positive feedback for FDSM［J］.Analog Integrated Circuits and Signal Processing，2012，73(1):151－160.

［15］ VISHNEVSKIIV M，LYAKHOV A I.A wireless LAN under high load and in noise:Throughput evaluation［J］.Automation and Remote Control，2001，62(8):1298－1311.

［16］ RYU G B，KIM SS.Numerical analysis on thickness effect ofmagnetic thin filmson conduction noise absorption inmicrostrip line［J］.Metals and Materials International，2011，17(5):805－810.

［17］ MARBURG S.Developments in structural－acoustic optimization for passive noise control［J］.Archives of Computational Methods in Engineering，2002，9(4):291－307.

［18］ 曲佳佳.风力机叶片气动载荷的计算方法研究［D］.北京:中国科学院研究生院(工程热物理研究所)，2014.QU Jiajia.Research on the Calculation Methods ofWind Blade＇s Aerodynamic Load［D］.Beijing:Institute of Engineering Thermophysics，Chinese Academy of Sciences，2014.

［19］ 赵丹平，田德，韦丽珍，等.风力机叶片模型气动载荷的仿真与实验［J］.东北林业大学学报，2009，37(12):124－125.ZHAO Danping，TIAN De，WEILizhen，etal.Simulation and experimentof aerodynamic load forwind turbine blademodel［J］.Journal of Northeast Forestry University，2009，37(12):124－125.

［20］ 吴继飞，徐来武，范召林，等.开式空腔气动声学特性及其流动控制方法研究［J］.航空学报，2014，35:1－7.WU Jifei，XU Laiwu，FAN Zhaolin，et al.Investigation on aero－acoustic characteristics and flow controlmethod of open cavity flow［J］.Acta Aeronautica Et Astronautica Sinica，2014，35:1－7.

［21］ 张成军，阴妍，鲍久圣，等.多源信息融合故障诊断方法研究进展［J］.河北科技大学学报，2014，35(3):213－221.ZHANG Chengjun，YIN Yan，BAO Jiusheng，et al.Research progress in fault diagnosismethods based on multi－source information fusion［J］.Journal of Hebei University of Science and Technology，2014，35(3):213－221.

［22］ 王华君，何涛，王风涛，等.飞轮储能在风力发电系统中的应用及仿真［J］.河北工业大学学报，2014，43(4):16－21.WANG Huajun，HE Tao，WANG Fengtao，et al.Matlab/Simulink vesearch on stablity ofwind farms using flywheel energy storage system［J］.Journal of Hebei University of Technology，2014，43(4):16－21.