台风作用下大型风电结构破坏模式研究

章子华，周 易，诸葛萍

(宁波大学 建筑工程与环境学院，浙江 宁波 315211)

作为再生清洁能源的风能越来越受到重视。我国风力发电产业起步较晚，但发展迅速，累计装机容量已达62.36 GW，跃居世界第一位[1]。

受季风影响我国东南沿海地区虽风能资源丰富，但易遭台风侵袭。台风风力较小时，风机运转时间长，能提高风电场发电量；而受台风正面袭击的风电场有可能遭受巨大损失。如2003年登陆广东的“杜鹃”台风，造成汕尾红海湾风电场13台风力发电机损坏，损失上千万元；2006年台风“桑美”登陆浙江苍南县直接穿过苍南风电场，致28台风力发电机组全部受损，其中5台倒塌，损失惨重。2010年台风“鲇鱼”正面登陆福建漳浦县六鳌镇，造成风电场三期Z13号风电机组倒塔、Z10号机组叶片折断[2]。

台风作用的风力发电结构破坏模式主要有：①由风致扭转颤振引起的叶片破坏,见图1(a)。在强湍流台风荷载作用下，柔长风机叶片发生强烈扭转振动，距扭心最远的叶片后缘扭转剪力达最大值并开裂，裂缝沿叶片横向向主梁扩展，达主梁与翼板交接处因刚度突变而改沿交界线纵向扩展[3]；②由极限风荷载引起的塔架失稳破坏及强度破坏，见图1(b)。强台风正面袭击风电场时，极限风速可达50 m/s甚至70 m/s以上。极端风荷载及由风载脉动分量引起的结构共振效应，即便采取停机顺桨措施，也难有效控制风电结构动力响应，不同塔筒段连接处、塔筒与基础连接处及下段塔筒极易发生失稳或强度破坏；③因桩基整体抗倾覆能力或单桩抗拔力不足引起的塔架倾覆，见图1(c)。在极限风速作用下风电结构倾覆力矩、迎风面及背风面最不利单桩竖向轴力达最大值。倾覆力矩超过桩基抗倾覆力矩则会发生塔架倒塌；单桩内力超过桩身正截面受拉承载力设计值则会发生单桩拉断破坏。

(a) 叶片破坏 (b) 塔筒失稳 (c)塔架倾覆

对台风所致风力发电结构破坏的有关研究已展开。汤炜梁等[4]总结计算风力机塔筒静态强度的理论方法，计算1.5 MW风力机在风速0～75 m/s变化过程中塔筒的静态挠度、弯矩、弯应力及底端连接螺栓的拉应力。为便于计算，该模型对叶轮、塔筒、基础进行较大简化，忽略动风荷载引起的放大效应，此对实际结构偏不利。文献[5]用Davenport脉动风速谱模拟某沿海风电场内的风场分布，建立风电结构-基础有限元模型，获得随机风荷载作用下风电结构动力响应时程。研究表明，结构低阶振型对风致振动响应有重要影响。徐旭等[6]基于台风基本特性，采用随高度变化的田浦台风风速谱及不随高度变化的石沅台风风速谱，用数值模拟方法仿真获得到与电视塔塔高竖向相关的43条脉动风速时程样本，为高耸结构抗台风研究提供重要借鉴。

本文用不随高度变化的台风脉动风功率谱，基于线性滤波法及竖向相关性简化表达式模拟某沿海风电场台风风场，建立风电结构-基础耦合有限元模型，计算台风荷载作用下风电结构动力响应，分析风电结构主体的可能破坏模式。

1 台风特性与风速条件

与高耸结构相关的台风特性包括台风平稳性与非平稳性、台风湍流脉动特性如湍流强度、湍流积分长度及阵风因子、功率谱密度函数及空间相关函数等[6]。本文所选台风风谱由平稳随机信号数据获得。

据国家标准《热带气旋等级》(GB/T 19201-2006)规定，热带气旋按中心附近地面最大风速划分为6个等级，其中底层中心附近最大平均风速32.7～41.4 m/s为台风，底层中心附近最大平均风速41.5～50.9 m/s为强台风，底层中心附近最大平均风速≥51.0 m/s为超强台风。本文研究的风电场位于江苏东台滩涂软土地区，中心位置东经120°54’，北纬32°47’。区内风资源分布东高西低，沿海边有一狭长风速急变带。受冬、夏季风影响，风力资源丰富，但也受台风威胁。本文取台风最大平均风速上、下限值v10=41.4 m/s，v10=32.7 m/s分别计算分析。

2 台风荷载模拟

2.1 脉动风功率谱

我国《建筑结构荷载规范》(GB5009-2012)采用沿高度不变的Davenport水平脉动风速谱，即

(1)

式中：Sv(n)为脉动风速功率谱；k=0.00215[7]为地面粗糙系数；v10为该地10 m高度处平均风速，m/s；n为脉动风频率；x为湍流积分尺度系数。

虽Davenport谱据大量不同地点、不同气候环境风速资料统计回归所得，具有一定代表性，但特定地区台风脉动风功率谱与普通脉动风功率谱有所区别。由于缺乏实测资料，本文暂用基于计算机拟合不随高度变化的台风脉动风功率谱经验公式[8],即

(2)

式中：各符号同式(1)。

2.2 脉动风时程模拟

统计分析实测风速样本可知，脉动风可用具有零均值的高斯平稳随机过程模拟，各态遍历性明显。脉动风场模拟方法中主要有谐波合成法(WAWS法)及线性滤波法(AR法)。Gerch等[9]提出将线性滤波技术用于生成时间序列等工程问题。AR法实质为由前面连续时刻随机量推导后面特定时刻随机量，具有时间相关性。模拟脉动风速时程公式为

(3)

式中：X，Y，Z为坐标向量矩阵；p为模型阶数；Δt为模拟风速时程时间步长，s；Ψk为自回归系数矩阵；N(t)为零均值且有给定协方差的正态分布随机过程。

2.3 空间相关函数

风电结构为典型高耸结构，空间中竖向相关性最显著，水平相关性不明显。考虑脉动风竖向相关性后所得空间点脉动风速时程更符合实际。常用相关函数有Davenport、日本规范AIJ建议及盐谷Shiotani。本文采用脉动风空间相关性简化表达式[10]为

(4)

式中：Lz值建议取60。

不同空间点脉动风风速功率谱密度间存在关系

(5)

式中：Sii(ω)，Sjj(ω)为点i，j脉动风自功率谱；Sij为点i，j互功率谱密度；ρij为点i，j相关性系数。

2.4 模拟风谱生成与验证

据上述理论编制程序生成风电场不同高度处台风脉动风速时程。限于篇幅，仅给出最大平均风速v10=32.7 m/s时风电场25 m、45 m及65 m处模拟台风脉动风速时程曲线，见图2～图4。

为验证模拟风谱的准确性，将模拟台风脉动风功率谱与目标功率谱(式(2))进行比较，见图5。由图5可见模拟功率谱与目标功率谱基本吻合，说明本文方法生成的模拟台风风谱在频域内的能量分布与实际台风风谱基本一致。

图2 25 m处模拟台风脉动风速时程

图5 台风脉动风功率谱密度(v10=32.7 m/s)

2.5 体系系数

风电结构主体为锥筒形构筑物，表面光滑，风荷载体型系数可按《建筑结构荷载规范》(GB50009- 2001)确定。据塔筒上、中、下三段平均直径与塔筒高度(H=63 m)比值，按线性内插法确定体型系数见图6[5]。

图6 塔筒体型系数μs

2.6 风荷载时程

据伯努利方程可得标准大气压、常温、干燥条件下风速风压关系式为

(6)

式中：w为风压值，kN/m2；γ为空气容重，kN/m3；g= 9.8 m/s2为重力加速度；v为风速，m/s。

塔筒任意高度处风荷载时程曲线计算式为

F(t)=μs(zi)A(zi)w(zi,t)

(7)

式中：μs(Zi)为Zi高度处体型系数；A(Zi)为Zi高度处迎风面积，m2；w(Zi,t)为Zi高度处风压，kN/m2。

台风天气风力机处于停机状态，作用于叶轮的水平轴向力[11]为

Fh=CDDρv2A

(8)

式中：CDD为阻力系数，取1.1；ρ为空气密度，kg/m3；v为风速，m/s；A为叶片迎风面积，m2。本风机叶片长37.5 m，平均宽度1.73 m，叶片平均迎风面积约65 m2。

综上所述，将考虑竖向相关性模拟风速时程曲线沿环向输入塔筒不同高度及叶轮形心处，即可据风速风压关系、体型系数及迎风面积等参数求出作用于塔筒、叶轮的风荷载时程曲线。

3 结构动力响应分析

3.1 有限元模型

风电场所用风力发电机轮毂高度为65 m。塔筒高度62.75 m，塔顶壁厚16 mm，塔底壁厚26 mm，采用空间壳单元shell63模拟。机舱(包括轮毂、叶片)总质量91.211 t，质心偏离塔筒轴线1.5 m，采用三维实体单元solid45模拟，自由度与塔筒顶部固接耦合。风机采用桩基础，混凝土承台以下共30根预制PHC桩，外径16.8 m，分布24根，内径4.1 m，分布6根，沿承台周向均匀布置，见图7。预制PHC桩采用三维弹簧单元combin14模拟。据试验确定单桩水平、竖向刚度分别为Eh=1.8×107N/m2，Ev=1.9×108N/m2。风电结构-基础耦合有限元模型见图8，该模型共14898个单元。材料参数见表1。其中机舱密度为据机舱总质量、体积计算的等效密度。塔筒材料为Q345-D钢，屈服强度345 MPa。

图7 桩基础布置

表1 材料参数

3.2 风电结构自振特性

风电结构自振特性见表2。计算结果表明，风电结构自振周期与台风脉动风变化周期(几秒到十几秒)接近，易产生共振效应。风电结构前4阶振型见图9。由图9看出，风电结构塔筒振动形式为前后弯曲、侧向弯曲及扭转。其中1阶振型为前后弯曲振动，2阶为侧向弯曲振动，3阶振型为2阶前后弯曲振动，4阶振型为2阶侧向弯曲振动。

表2 风电结构自振频率(单位：Hz)

图9 风电结构前4阶振型

3.3 结构风载动力响应

3.3.1 结构动力方程求解

结构动力平衡方程为

MΔu+CΔu+KΔu=ΔF

(9)

式中：M为结构质量矩阵；C为结构阻尼矩阵；K为结构总刚度矩阵；ΔF为增量荷载矩阵。

将每个时间步的台风荷载沿周向输入塔架不同高度及叶轮轮毂形心处，用修正的Newton-Raphson迭代法可得风电结构动力响应。需指出的是，GL规范与国内风电机组相关规范均未给出动力分析时阻尼比的建议值。马人乐等[12]通过环境脉动实测发现1.5 MW、轮毂高度65 m处风电结构1阶平动阻尼比约为1.75%。故本文计算时采用该实测值。

3.3.2 塔架变形及轮毂位移

对风电结构而言，塔架顶端水平位移对结构整体稳定性、安全性至关重要。最大平均风速分别为v10= 32.7 m/s及v10=41.4 m/s时，65 m轮毂处水平位移时程曲线见图10。由图10看出，v10=32.7 m/s时轮毂水平位移变化范围-0.203 6～1.221 1 m。v10=41.4 m/s时轮毂水平位移变化范围-0.327 4～2.027 2 m。按《高耸结构设计规范》(GB50135-2006)规定：在以风为主的荷载标准作用下，按非线性分析的高耸结构任一点水平位移不得大于该点离地高度的1/50，即该结构65 m高轮毂处的水平位移限值为1.3 m。因此，v10= 41.4 m/s时的风电结构水平位移已超出限值。

3.3.3 塔筒应力

在水平风荷载作用下风电结构塔筒变形呈弯曲型，竖向应力呈上小下大分布。虽塔筒壁厚按上薄下厚原则设计，但下段塔筒仍为屈曲危险区域。选塔筒底端迎风面单元及背风面单元为研究对象，v10=32.7 m/s时各单元应力时程曲线见图11、图12。由二图看出，t=75.75 s时塔筒迎风面、背风面底端出现最大Von-Mises等效应力，分别达947 MPa、913 MPa，远超塔筒钢材屈服强度345 MPa。因此，该风电结构塔筒已进入塑性变形阶段，损坏性极大。t=75.75 s、v=32.7 m/s时塔筒底部应力分布见图13、图14。由二图看出，塔筒屈服区域集中在塔筒底端约4m高度范围内，呈锥形分布。塔筒底端与基础预埋环连接处最不利，需重点校核塔筒与基础连接的螺栓强度。

图10 轮毂中心水平位移时程

图13 塔筒底部迎风面等效应力分布

图14 塔筒底部背风面等效应力分布

3.3.4 单桩轴力

在水平风荷载作用下风电结构桩基础通过弯矩分配承担不同程度倾覆弯矩。扣除自重荷载效应后，外圈桩迎风面处单桩竖向拉力达最大，背风面处单桩竖向压力达最大，上述2桩为最不利单桩。迎风面、背风面最不利单桩的轴力时程曲线见图15、图16。由二图看出，v10=32.7 m/s时迎风面最不利单桩竖向轴力变化范围-38～446 kN，背风面最不利单桩竖向轴力变化范围327～811 kN；v10=41.4 m/s时迎风面最不利单桩竖向轴力变化范围-313～486 kN，背风面最不利单桩竖向轴力变化范围286～1 085 kN。据设计资料试算，单桩竖向极限承载力标准值为1 250 kN，抗拉极限承载力标准值为470 kN，桩身正截面受拉承载力标准值为2 500 kN，最不利单桩轴力、内力均在安全范围内。但值得注意的是，迎风面最不利单桩竖向轴力出现拉压交替现象，加之风电结构在脉动风作用下振动强烈，对桩周土体扰动较大，可能导致桩基础极限承载力下降。

图15 迎风面最不利单桩轴力

图16 背风面最不利单桩轴力

4 结 论

(1) 本文采用不随高度变化的台风脉动风功率谱，基于线性滤波法及竖向相关性简化表达式模拟某沿海风电场台风风场模型。经验证，该模拟风功率谱密度与目标风谱较吻合。

(2) 在建立风电结构-基础耦合有限元模型及模拟风谱基础上计算风电结构在上、下限台风风速作用下的动力响应。对该风电结构而言，最可能发生的破坏模式为塔筒屈服。平均风速v10=41.4 m/s时虽塔顶水平位移超出规范限值，但单桩承载力安全余度较大，说明该结构布桩形式有较多优化空间；而部分单桩轴力发生的拉压交替变化会致桩基础承载力下降。

(3) 须重视塔筒进入塑性阶段后的屈服变形，以便准确掌握风电结构在极端风况下的倒塌规律。对不同地区的台风实测资料，需通过统计分析获得符合当地地形的台风风谱模型。

[1] 李俊峰.中国风电报告2012[M].北京：中国环境科学出版社, 2012.

[2] 王力雨, 许移庆.台风对风电场破坏及台风特性初探[J].风能, 2012, 5:74-79.

WANG Li-yu, XU Yi-qing.A preliminary study on typhoon damage to wind farm and the typhoon characteristics[J].Wind Energy, 2012, 5:74-79.

[3] 王景全,陈政清.试析海上风机在强台风下叶片受损风险与对策-考察红海湾风电场的启示[J].中国工程科学, 2010, 12(11): 32-34.

WANG Jing-quan, CHEN Zheng-qing.Analysis of risks and measures on the blade damage of offshore wind turbine during strong typhoons: enlightenment from red bay wind farm[J].Engineering Science, 2010, 12(11): 32-34.

[4] 汤炜梁,袁奇,韩中合.风力机塔筒抗台风设计[J].太阳能学报, 2008, 29(4): 422-427.

TANG Wei-liang, YUAN Qi, HAN Zhong-he.Withstanding typhoon design of wind turbine tower[J].Acta Energiae Solaris Sinica, 2008, 29(4): 422-427.

[5] 章子华,王振宇,刘国华.风电场脉动风模拟及风机塔架动力响应研究[J].太阳能学报,2011,32(7): 992-998.

ZHANG Zi-hua, WANG Zhen-yu, LIU Guo-hua.Simulation of fluctuating wind in wind farm and dynamic response of wind turbine tower[J].Acta Energiae Solaris Sinica, 2011, 32(7): 992-998.

[6] 徐旭,刘玉.基于台风风谱的电视塔风场数值模拟[J].特种结构, 2008, 25(2): 39-43.

XU Xu, LIU Yu.Wind field numerical simulation of television tower based on the typhoon spectrum[J].Special Structures, 2008, 25(2):39-43.

[7] 王修琼,崔剑峰.Davenport谱中系数K的计算公式及其工程应用[J].同济大学学报, 2002, 30(7): 849-852.

WANG Xiu-qiong, CUI Jian-feng.Formula of coefficientK in expression of Davenport sprectrum and its engineering application[J].Journal of Tongji University, 2002, 30(7): 849-852.

[8] 石沅, 陆威, 钟严.上海地区台风结构特征研究[A].第二届全国结构风效应学术会议论文集[C].1988:106.1-106.12.

[9] Gerch W, Yonemoto J.Synthesis of multivariate random vibration systems: a two-stage least squares ARMA model approach[J].Journal of Sound and Vibration, 1977, 52(4): 553-565.

[10] Shiotani M, Iwatani Y.Correlations of wind velocities in relation to the gusting loadings[C].Proceedings of the Third International Conference on Wind Effects on Buildings and Structures, Tokyo, 1971: 57-67.

[11] 李庆宜.小型风力机设计[M].北京: 机械工业出版社, 1986.

[12] 马人乐,马跃强,刘慧群,等.风电机组塔筒模态的环境脉动实测与数值模拟研究[J].振动与冲击, 2001, 30(5): 152-155.

MA Ren-le, MA Yue-qiang, LIU Hui-qun,et al.Ambient vibration test and numerical simulation for modes of modes of wind turbine towers[J].Journal of Vibration and Shock,2001,30(5):152-155.

[13] 汤炜梁,袁奇,高锐.风力机塔筒的三维有限元抗台风分析[A].中国动力工程学会透平专业委员会2007年学术研讨会论文集[C].2007.