风力发电机雷击次数的改进计算模型与方法

刘 健， 杨仲江， 卢慧慧

(南京信息工程大学 大气物理学院，江苏 南京 210044)

风力发电机雷击次数的改进计算模型与方法

刘 健， 杨仲江， 卢慧慧

(南京信息工程大学 大气物理学院，江苏 南京 210044)

随着风能在中国能源比例中的持续增长，对风力发电机可靠性的要求变得越来越高。因此，风力发电机的防雷保护已成为重要课题。准确的风电发电机雷击次数是确定合理防雷措施的前提，针对现有IEC风机年雷击次数计算公式存在的问题，利用蒙特卡洛方法计算不同引雷半径公式下的风机下行雷击次数; 同时借鉴Rizk提出的高层建筑物上行雷击次数计算模型，应用于风机上行雷击次数计算。仿真分析结果表明：风机下行雷击次数和上行雷击次数均随风机高度增加而增大，下行雷击占总雷击事件比例随风机高度增加而下降；IEC公式计算得到的风机年雷击次数偏大，改进后的年雷击次数计算方法较为符合实际观测情况，可为风场防雷分析提供可靠依据。

风力发电机;年雷击次数;等效引雷半径;上行雷

近些年来，中国风电产业得到迅猛发展。根据统计[1]，截止到2014年底，中国累计风电装机容量11 460.9万kW, 并网装机容量达到9 637万kW，年总发电量近1 900亿kW·h。为了尽可能地接收风能，风机一般安装在山顶或者非常开阔的地带，如草原、沿海滩涂、沙漠戈壁等。雷电是威胁风力发电场安全的主要问题之一，同时，风机的高度和叶片的长度还在不断增加，进一步增大了风机被雷电击中的概率。美国一项长达5年的观测数据指出[2]，508台风机每台平均8.4年遭受一次雷击，以风机寿命20年计算，其运行服役期间因雷击导致的叶片损伤将发生2～3次。中国风电场也发生过较多的雷电灾害，海南东方风电场的数据显示[3]，仅风机桨叶的雷击损坏率就达到5.56片/(百片·年)。

为了减少风电机组因雷电灾害造成的损失，有必要针对风机的雷击风险进行计算，从而为风电场前期的选址过程中考虑雷电灾害因素提供决策依据。目前，针对风力发电机遭受雷击次数的计算主要依据IEC 61400—24等国际标准[4]。但是IEC给出的计算方法对于传统的固定对象较为适用，而风机系统始终处于一个动态运行过程中，使得雷电击中风机的情况变得非常复杂，传统的计算方法不能给出较为准确的雷击次数。

笔者介绍IEC风机年雷击次数计算公式，分析公式计算存在的问题。考虑风机有效高度的变化，利用蒙特卡洛模拟计算不同引雷半径公式下风机下行雷击次数。同时，利用Rizk提出的高层建筑物雷击发生率计算模型[5]计算风机上行雷击次数，分析上行雷击事件发生特点。最后，将改进计算结果与实际观测数据进行比较，为风力发电机遭雷击的风险分析提供依据。

1 现有年雷击次数计算方法

风力发电机的防雷设计应该基于其被雷电击中的频率，这个频率与风机高度、雷电活动特征和地形环境有关。IEC[4]给出的风力发电机年雷击次数计算公式如下：

N=NG·Cd·π·Ra2·10-6。

(1)

式中NG为雷击大地密度，次/(km2·a)；Ra为风机的等效引雷半径，m；Cd为风机的位置因子。

IEC[4]给出的风机等效引雷半径定义为风机有效高度的3倍：

Ra=3hs。

(2)

式中hs为风机的有效高度，m。IEC[4]建议平原区的风机有效高度取轮毂高度加叶片半径。

Eriksson[6]通过对世界不同地区平原上高度范围在20～540 m的雷击事件的观察，总结了年雷击次数计算的经验公式：

N=NG·24·h2.05·10-6。

(3)

不考虑位置因子时，Eriksson经验公式和IEC公式是类似的，即

(4)

2 存在问题

2.1 引雷半径

雷击等效截收面积和引雷半径有关，引雷半径是指引雷的结构物存在一定吸引半径，如果雷电下行先导进入吸引半径范围内，结构物产生迎面先导以拦截下行先导，否则雷电先导击中地面。IEC标准中引雷半径取风机有效高度的3倍(式(2))，但是相关研究[7]认为等效引雷半径与结构物高度和雷电流有关，这点可以通过电气几何模型(EGM)加以推导验证。基于电气几何模型的引雷半径计算示意如图1所示，计算公式如下：

(5)

式中St为雷击塔筒击距，m；Sg为雷击大地击距，m。计算公式如下[7]：

St=10I0.65，

(6)

Sg=0.9rs。

(7)

将式(6)、式(7)代入式(5)得到引雷半径具体计算公式如下：

(8)

由式(8)可知，风机引雷半径Ra与风机高度h和雷电流幅值I直接相关。IEC标准中未考虑雷电流影响是不合理的。

图1 引雷半径计算

Borghetti[8]和Rizk[5]根据各自研究都提出了相应引雷半径计算公式：

Ra=0.028hI+3h0.6，

(9)

Ra=4.27I0.55h0.41。

(10)

30 kA雷电流情况下不同引雷半径随高度的变化如图2所示，可以看出，引雷半径随着风机高度的增加而增大，但IEC公式计算结果明显大于其他计算公式，高度越高，差异越明显，因为其未考虑雷电流因素影响。引雷半径计算结果偏大会造成年雷击次数的过估计，从而增加防雷保护成本。

图2 引雷半径随高度变化曲线(30 kA雷电流情况下)

2.2 风机有效高度

风机在转动过程中有效高度处于不断变化当中。假设在风机运转在过程中转速保持不变，叶片转动的角速度为ω，风机机舱高度为h0，叶片长度为l，θ为叶片旋转相位角，风机有效高度变化的函数为[9]

(11)

θ=αi+ωt。

(12)

t=0时，叶片的初相位分别为α1=π/2，α2=π/2+2π/3，α3=π/2+4π/3。风机有效高度在一个转动周期内的变化如图3所示，其在一个转动周期内变化3次，最小值为h0+l/2，最大值为h0+l。图3中风机机舱高度h0取80 m，叶片长度l取40 m。

图3 风机有效高度变化曲线

IEC公式在计算风机年雷击次数时，有效高度取定值而未考虑其随时间的变化。如图4(a)所示，当θ=π/2时，有效高度为120 m，根据IEC公式计算雷击等效截收面积Ad1=407 150 m2，当θ=2π/3时(图4(b))，有效高度为100 m，等效截收面积Ad2=282 743 m2，两者截收面积相差高达30%。因此，在计算风机的雷击次数时，必须考虑有效高度的变化，不能简单取定值。

图4 风机有效高度最高点及最低点

2.3 上行雷击

根据相关研究，建筑物越高，越容易触发上行闪电。同样由于风机的存在，在叶片尖端电场会产生严重畸变。畸变的电场强度足以产生一次上行先导，拦截雷电下行先导，使得叶片尖端成为雷击点。根据上行先导起始判据[5]计算感应电动势，来分析风机高度对上行雷击事件的影响。

(13)

式中Ui为风机叶片尖端电位，kV；R为风机高度，m；x0为是与上行先导速度成正比的参数，Rizk[5]取5 m；E∞是先导准稳态场强，取3 kV/m；Ei是正极性先导起始场强，取400 kV/m[10]。

高度为80 m的风机，叶片长度为40 m时，由于雷暴云形成下行先导的影响，在风机尖端产生的感应电动势可以高达1 535 kV，当叶片长度为60 m时，风机顶端产生的感应电动势为1 547 kV。分析表明：风机越高，在叶片尖端产生感应电动势越大，就越容易触发上行雷击。

此外，风机叶片的转动也影响了周围电场，更有利于上行先导的产生[11-12]。叶片的转动避免了叶片尖端空间电荷的积累。在叶片静止情况下，地面电场强度达到一定程度时产生起始电晕，进而产生上行先导。电晕产生的空间电荷抑制了叶片周围电场的增强。正极性空间电荷的移动速率比电子速率低得多，不能突破电晕层屏蔽。叶片转动时，如果叶片转动速率大于电晕产生的小粒子的迁移速率，正极性电荷能够从电晕层逃逸出来。叶片尖端的电荷量累积小于静止状态下，叶片尖端电场也越强，越有利于上行先导的产生。

IEC和Eriksson公式计算得到的风机年雷击次数包含了全部雷击事件，无法区分上行雷击次数和下行雷击次数，没法分析风机高度对上行雷击事件的影响。

3 计算方法改进

3.1 下行雷击次数

风机遭到雷击的过程是一次随机的过程，受到雷电流、风机高度、叶片旋转相位角、雷击点方位等因素影响。为了解决IEC公式中无法考虑雷电流因素和有效高度变化的影响，运用蒙特卡洛法[13]进行雷击过程的模拟。通过对多次模拟结果的统计，得到年雷击次数估计值。

考虑到可能影响雷击风机的随机过程, 将击中概率y作为若干随机变量的函数：

(14)

式中i为雷电流幅值；(x，y)为雷击点方位；θ为叶片旋转相位角。

对这些参数进行随机抽样，第k次抽样的结果记做：

(15)

如果雷电下行先导离风机距离小于风机引雷半径，则认为击中风机。若第k次抽样满足击中风机的条件，则记yk=1，否则yk=0。进行N次独立重复试验，当N充分大时，根据大数定理，击中风机的概率为

(16)

雷电流幅值i的随机抽样依据中国规程[14]推荐雷电流幅值分布概率公式:

(17)

式中i为雷电流幅值，kA；P为雷电流幅值大于i的概率。每次产生一个服从[0，1]均匀分布随机数，便可根据雷电流幅值分布概率公式随机获得一个雷电流幅值。

Na=NGPA。

(18)

以美国堪萨斯某风场风机为例[15]，塔筒高度90 m，叶片长度为35 m。风场所在地区平均雷击大地密度为2.98次/(km2·a)。叶片长度为35 m时，不同引雷半径公式下风机下行雷击次数随塔筒高度变化如图5所示,可以看出，不同引雷半径公式计算得到的风机年雷击次数随高度变化趋势一致。塔筒高度越高，风机越容易遭受雷击。

图5 风机下行雷击次数随塔筒高度变化

在塔筒高度为90 m时，Borghetti公式计算得到的下行雷击次数为0.69次/年，Rizk公式计算结果为0.80次/年，EGM推导公式计算结果为0.37次/年。根据堪萨斯风场4.5年的场地观测数据统计得到每台风机年平均下行雷击次数为0.92次/年[15]，与Rizk引雷半径公式计算结果最为接近。EGM推导公式计算结果和其他公式相差较大，因为高度较高时，EGM模型引雷半径计算明显偏小。

3.2 上行雷击次数

为了进一步分析上行雷击事件发生比例，应用Rizk提出的高层建筑物上行雷发生率计算模型[6]计算风机上行雷击次数：

Nu=25NG0.8·e(-(Egc-Ego)/Eg1)，

(19)

式中Eg0为起始场强；Eg1为形成先导场强，普通雷暴情况下均取2 kV/m。Egc为地面临界场强，与建筑高度有关：

(20)

考虑上行雷击后，总雷击次数的计算公式为

Nf=Nd+Nu

(21)

叶片长度为35 m时，风机上行雷击次数、总雷击次数随塔筒高度的变化如图6所示，其中，下行雷击次数计算采用了Rizk公式。同时与IEC年雷击次数公式进行比较。

在塔筒高度为90 m时，上行雷击次数为0.23次/年，IEC公式计算得到的总雷击次数为1.38次/年，风场实际观测数据为1.18次/年[15]，改进方法为1.03次/年，较为符合实际情况。由图6可以看出，随着高度的增加，上行雷击次数不断增加。在高度相对较低时，改进方法计算结果与IEC公式计算结果相差不大，随着高度的不断增加，差距随之增大。

图6 总雷击事件和上行雷击事件随高度的变化曲线

Eriksson通过观测给出了上行雷击事件发生比例的经验公式[6]：

Pu=52.8lnh-230

(22)

Eriksson经验公式和该文计算结果的比较如图7所示，可以看出，上行雷击事件比例随高度增加的趋势是一致的。Eriksson经验公式得到的90 m高塔筒，35 m长叶片上行雷击比例为24.9%，改进算法算得比例为22.3%，二者相差不大。

图7 上行雷击事件比例随高度的变化

4 结语

分析现有风机年雷击次数计算公式存在问题，探讨适合大型风力发电机的年雷击次数计算方法，得到如下结论：

1)IEC公式计算风机年雷击次数时引雷半径选取过大，无法考虑风机转动过程中有效高度的变化，无法给出具体上、下行雷击次数。

2)利用蒙特卡洛模拟考虑风机有效高度周期性变化特点，采用Rizk引雷半径公式计算风机下行雷击次数。利用Rizk高层建筑物上行雷发生率模型计算风机上行雷击次数，计算结果与实际观测数据较为接近。

3)风机下行雷击和上行雷击次数均随高度增加而增大，但下行雷击事件占总雷击事件比例随着高度增加而减小。

[1]郭子炘，李庆民，闫江燕，等．海上风电场雷击演化物理机制的研究综述[J]．电气工程学报，2015，10(5)：10-17． GUO Zi-xi,Li Qing-min,YAN Jiang-yan, et al. Summary of research on physical evolution mechanism of lightning discharge of offshore wind farms[J]. Journal of Electrical Engineering, 2015, 10(5): 10-17.

[2] Garolera A C, Madsen S F, Nissim M, et al. Lightning damage to wind turbine blades from wind farms in US[J]. IEEE Transactions on Power Delivery, 2014, 31(3): 1 043-1 049.

[3] 李谦，钟定珠，彭向阳．海南东方风电场II期风力发电设备故障及雷害分析[J]．风力发电，2010(3)：22-27． LI Qian, ZHONG Ding-zhu，PENG Xiang-yang. Analysis of fault and lightning damage of wind power equipment of Dongfang wind farm second stage in Hainan Province[J]. Wind Power, 2010(3): 22-27.

[4] IEC 61400-24-2010. Wind turbine generation system Part 24: Lightning protection[S].

[5]Farouk A M Rizk. Modeling of lightning incidence to tall structures part1: Theory[J]. IEEE Transactions on Power Delivery, 1994, 9(1): 162-171.

[6] Eriksson A J. The Incidence of lightning strikes to power lines[J]. IEEE Transactions on Power Delivery, 1987, 2(3): 859-870.

[7] Eriksson A J. An improved electrogeometric model for transmission line shielding analysis[J]. IEEE Transactions on Power Delivery, 1987, 2(3):871-885.

[8] Borghetti A, Nucci C A，Paolone M. Estimation of the statistical distributions of lightning current parameters at ground level from the data recorded by instrumented towers[J]. IEEE Transactions on Power Delivery, 2004, 19(3): 1 400-1 409.

[9] Dale Dolan, Charles Sao, Peter Lehn. Lightning exposure of wind turbines[C]. Conference on Electrical and Computer Engineering (CCECE), Canadian, 2006.

[10] Zhou H, Nelson Theethayi, Gerhard Diendorfer, et al. On estimation of the effective height of towers on mountaintops in lightning incidence studies[J]. Journal of Electrostatics, 2010, 68(5): 415-418.

[11] Long M, Garcia M B, Thottappillil R. On the interception of dart lightning leaders from wind turbine blades[C]. 2014 International Conference on Lightning Protection (ICLP), Shanghai，China, 2014.

[12] Montanyà J, Velde O, Williams E R. Lightning discharges produced by wind turbines[J]. Journal of Geophysical Research Atmospheres, 2014, 119(3): 1 455-1 462.

[13] 郝俊琦，李琳，王平．基于蒙特卡罗法的超高压输电线路绕击跳闸率的计算[J]．华北电力大学学报，2012，39(5)：35-41． HAO Jun-qi, LI Lin, WANG Ping. Calculation for shielding failure trip-out rate of ehv transmission line based on Monte Carlo method[J]. Journal of North China Electric Power University, 2012, 39(5): 35-41.

[14] GB 50064—2014.交流电气装置的过电压保护和绝缘配合设计规范[S].

[15] Anna C G. Lightning protection of flap system for wind turbine blades[D]. Copenhagen: Technical University of Denmark, 2014.

Improved calculation method of annual lightning flash number for wind turbine

LIU Jian， YANG Zhong-jiang， LU Hui-hui

(The School of Atmospheric Physics，Nanjing University of Information Science and Technology, Nanjing 210044, China)

With the increasement of wind power generation, reliability of wind turbine becomes more and more important. Lightning protection of wind turbine has become a hotspot recently. Accurate annual lightning flash number determination is on the premise of proper lightning protection measurements. The formula for the number of lightning flashes according to IEC was introduced and its problems were also analyzed. Combined with the rotation of wind turbine blades, the annual number of lightning flashes under different equivalent attractive radius formulas was calculated by the means of Monte Carlo simulation. Model proposed by Rizk for assessment of lightning stroke incidence in high building was applied to calculate the number of upward flashes for wind turbines. The simulation results show that the annual number of downward flashes and upward flashes both increase with the height of the tower, while the percentage of the downward flashes decreases. The calculation method proposed in this paper is in good agreement with actual situation compared with the IEC formula, which could provide reliable evidence for lightning protection in wind farms.

wind turbine; annual number of lightning flashes; equivalent attractive radius; upward flashes

2016-04-21

国家自然科学基金(41175003)；江苏高校优势学科建设工程项目(PADA)

杨仲江(1961—)，男，副教授，主要从事雷电灾害防护技术的研究；E-mail：liujiannuist@163.com

TM83

A

1673-9140(2016)04-0149-06