臺風作用下大型風電結構破壞模式研究

章子華，周 易，諸葛萍

(寧波大學 建筑工程與環境學院，浙江 寧波 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.