基于自激冰振的風力機海冰載荷分析

葉柯華， 李 春，2， 楊 陽， 孫 瑞， 張萬福，2

(1.上海理工大學 能源與動力工程學院，上海 200093；2.上海動力工程多相流與傳熱重點實驗室，上海 200093)

風力發電是最主要的新能源利用形式之一，正處于高速發展中。因為對風力發電持續的高需求，風力機作用風力發電載體，一方面技術尤其以大型化和離岸化得到迅猛發展[1-2]；另一方面裝機容量和安裝數目的高速增長，優質風能開發區日趨飽和，復雜嚴苛安裝環境成為不得不面臨的選擇。嚴寒海區已為風力機選擇之一。此處，兆瓦級風力機面臨氣動載荷、波浪載荷和海冰載荷等共同作用，導致整機處于嚴峻運行工況中，產生復雜的內應力和結構振動，承受劇烈的疲勞損耗中。

風力機為剛柔耦合多體結構，長期處于非定常交變載荷作用中，涉及多物相強耦合過程。Jonkmant等[3]采用葉素動量理論計算風力機湍流風況下氣動載荷，同時考慮波浪載荷、地震載荷和海冰載荷等作用，建立風力機多體動力學模型，對風力機整機動態響應進行仿真，并以此開發FAST系列程序。眾多研究主要針對風力機整機仿真和氣彈耦合問題，考慮葉片和塔架等大位移變形對氣動性能影響[4-5]。風力機多種載荷的研究，除考慮氣動載荷[6]和波浪載荷[7]外，集中在地震載荷[8]、風沙載荷[9]和海冰載荷等[10]。文獻 [11-12]對波浪載荷作用下海上風力機穩定性結構設計進行了系列研究。丹麥和加拿大研究機構對不同冰況對風力機作用情況進行實驗和模擬，得出抗冰錐角分別為55°和65°時，塔基剪切力最大為0.54 MN和0.63 MN[13-14]。

海冰載荷是嚴寒海區特有且常見的載荷，具有作用力大、持續時間長、振動特性復雜等特點。相比靜態海冰載荷作用，動態海冰載荷能導致海工平臺發生的冰激振動，具有更大的破壞力。1969年，我國渤海海域“海井一號”石油平臺因冰激振動發生傾覆[15]。1986年，俄羅斯庫頁島東北海域Molikpaq石油平臺受冰激振動影響，地基海床發生“砂土液化”，平臺嚴重損毀[16]。故嚴寒海區風力機受海冰載荷的影響不容忽視。

針對兆瓦級近海風力機在動態海冰載荷作用下冰激振動特性展開研究，本文以NREL 5 MW近海3葉片水平軸風力機為研究對象，支撐結構采用單樁式，采用Kaimal風速譜建立湍流風況作用氣動載荷條件，海冰載荷計算采用M??tt?nen 自激冰振模型。研究海冰載荷與氣動載荷聯合作用對風力機基礎結構的動力學響應的影響，為嚴寒海區風力機結構設計提供參考。

1 多體動力學模型

風力機為復雜的剛柔耦合多體結構，運行過程涉及氣彈和水彈等多物態強耦合，載荷主要包括：氣動載荷、波浪載荷、重力載荷和海冰載荷等。

風力機多體動力學方程為[17]

(1)

式中：M為風力機質量矩陣；C為風力機結構阻尼矩陣；K為風力機剛度矩陣；Faero為氣動力；Fhydro為波浪力；Fgravity為重力；Fice為冰-構作用力。

以凱恩(Kane)方法建立風力機多體動力學模型。凱恩(Kane)方法基于達朗貝爾原理，鑒于廣義坐標對非完整系的不獨立性，以廣義速率代替廣義坐標作為獨立變量描述系統運動。

凱恩動力學方程為

(2)

第Ni剛體，主動力施加在其質心Xi，FNi和MNi分別為作用在每個剛體Ni上的主動力與力矩，EωNi為剛體Ni角速度，EvNi為剛體Ni速度。

則風力機系統的廣義主動力為

(3)

廣義慣性力為

(4)

求解各參數時首先通過4階Adams-Beshforth預測-校正方法確定低階項的值，并以此構成方程的右邊項，然后采用高斯消元法求解系統自由度的加速度。計算得到的加速度值用于修正預測值并以此提高計算精度。迭代數次后采用4階Adams-Mounton預測-修正算法確定加速度值，并給出該時間步的最終解。由于該算法的非自發性，前4個時間步通過Runge-Kutta法求解。

2 研究對象和仿真方法

2.1 研究對象

選用NREL 5 MW 近海水平軸風力機為研究對象，采用單樁支撐，葉片數為3片。整機由：葉片、機艙、輪轂、塔架和樁柱等部分組成。機艙有安裝傳動系統、發電系統和控制系統等。葉片和塔架為柔性結構，采用超級單元方法建模。塔架為圓臺構型，底部直徑為5.6 m，頂部直徑為4.0 m。風力機結構如圖1所示。風力機參數如表1。

圖1 風力機結構Fig.1 The structure of wind turbine

海平面近似位于樁柱與塔架交界處，海冰在海平面上，作用于樁柱上部，塔架最下部。此處也是塔基剪切力計算位置。塔頂位移以大地坐標系為參考，取機艙和塔架連接處圓心。

塔架為懸臂梁，樁柱之下地基與海床剛性連接，限制各向自由度，地基位移忽略不計。

表1 風力機參數

2.2 仿真方法

仿真采用美國可再生能源實驗室(National Renewable Energy Lab，NREL)編譯的風力機整機動力學仿真程序FAST進行。湍流風場由FAST前處理程序根據Kaimal風速譜生成。以上程序皆為開源性質，采用Fortran語言編譯，通過動態鏈接庫實現數據調用和輸入輸出。風力機模型和控制方法完全參數化設計。動力學模型采用凱恩(Kane)方法。

海冰載荷和氣動載荷作用下在FAST環境下風力機動態響應仿真如圖2。海冰載荷和氣動載荷分別計算，相互獨立。

圖2 仿真方法Fig.2 The simulation method

3 海冰載荷

3.1 冰激振動

柔性海工平臺由海冰碰撞產生的交變載荷所導致結構發生一定幅值的穩態振動稱為冰激振動。相比靜態海冰載荷，冰激振動進一步加大結構疲勞損耗，造成內應力聚集，當振動頻率與結構固有頻率接近時，結構易發生失穩，造成毀滅性的破壞。

目前，國內外學者對冰激作用產生原因仍處于討論當中，觀點主要分為強迫振動理論和自激振動理論[18-19]。1966年， Peyton[20]依據室內冰池實驗和實地觀測，得出庫克灣海冰保持1 Hz破碎頻率，不受結構形式影響。1976年， Neil[21]進一步建立海冰載荷周期為冰排破碎長度和冰排前進速度的函數。1971年， Matlock等[22]采用彈簧-阻尼系統建立的柔性結構強迫冰振模型。但是實地觀測中，海工平臺發生穩態冰激振動時，海冰載荷并非必然存在周期性。故此， Blenkarn首次提出冰激振動是自激振動[23]。M??tt?nen[24]根據海冰強度隨加載速率變化規律進一步發展了自激振動理論。徐繼祖等[25]結合自激性和周期性提出了“冰力振子模型”。

強迫冰振理論主要依據實地觀測到的海冰載荷周期，再進行經驗公式推導[26]。

自激冰振理論則認為海冰載荷不具有周期性。低速冰況下海冰與結構接觸區域呈“漿狀”和“粉末狀”，近似為“黏性牛頓流體”。此碎冰層海冰強度與加載速率的非線性關系，導致了冰-構作用期間出現負阻尼[27]。而近海風力機多數處于低速冰況，因此，本文選擇自激冰振模型計算海冰載荷。

3.2 M??tt?nen 自激冰振模型

如圖3所示，海冰強度為加載速率的函數，隨加載速率增加歷經：延性區、過渡區和脆性區。在低加載速率區域，強度隨加載速率增加而增加，海冰主要為韌性破壞。當加載速率增加至一定程度，強度隨加載速度增加而降低，強度曲線出現“負坡段”，從而導致負冰阻尼。隨加載速率繼續增加，強度不再變化，此時海冰主要發生脆性破壞。

圖3 海冰強度隨加載速率變化Fig.3 Variation of the ice strength with loading rate

M??tt?nen 于1976年所提出模型是目前典型的自激冰振模型。

M??tt?nen 模型依據海冰強度與加載速率關系，建立冰-構作用力為冰-構相對速度的函數

(5)

冰阻尼系數φ為

(6)

根據第2節中風力機多體動力學方程，將冰-構作用力移至方程左邊，與結構阻尼共同構成系統阻尼

(7)

當進入“負坡段”時，且負冰阻尼絕對值大于結構阻尼絕對值將導致系統總阻尼為負，從而引發自激冰振發生。

3.3 冰-構作用力結果

自激冰振中海冰載荷同時由冰況和結構特性決定。本文為研究不同冰況對風力機動態響應影響，計算1.0 ～1.7 m，間隔0.1 m8種不同冰厚。各冰厚計算0.02 ～0.8 m/s，間隔0.02 m/s的40種不同冰速。

海冰對風力機造成的載荷最大值如圖4所示。

圖4 海冰載荷最大值Fig.4 The maxes of the ice loads

4 風場及氣動載荷

風場是氣動載荷計算前提，也直接影響氣動載荷計算準確性。風力機運行的湍流風場中風速具有時間非定常性和空間不均勻性。其中，風速脈動部分是平穩隨機過程，根據Kaimal風速譜建立。風速譜是功率譜，表征不同大小風速對脈動性貢獻大小。垂向風速變化依據風廓線，縱向風速變化采用空間相干模型[28]。

本文建立的風場大小為150 m×150 m，以輪轂高度(90 m)為參考建立平均風速11.4 m/s(額定風速)風場。圖5為輪轂處風速變化。

風作用于葉片上推動風輪旋轉，也產生阻力。氣動載荷計算采用葉素動量理論(Blade Element Momentum Theory, BEM)，根據葉素誘導速度，沿展向積分得出氣動載荷。

圖5 輪轂高度(90 m)風速時域變化Fig.5 The wind velocity on the height of hub(90 m)

仿真中氣動載荷均采用由Kaimal風場生成的11.4 m/s的湍流風況計算，組合1.0～1.7 m冰厚以及0.02 ～0.8 m/s冰速的320種不同冰況。

5 風力機結構動態響應分析

在FAST中分別計算NREL 5MW近海風力機在湍流風況和海冰作用下動態響應。仿真總時長為300 s，步長為0.01 s。

5.1 塔基載荷塔頂位移時域結果

塔基各向剪切力最大值變化如圖6所示。由圖6可知，塔基剪切力增加主要集中于中低冰速區域，增幅隨冰厚增加而減少。前后向剪切力在冰速到達0.1 ～0.2 m/s時，出現較大值；側向剪切力在冰速達到0.3～0.4 m/s時，出現較大值。在低冰速區，剪切力略有增加。隨冰速超過0.1 m/s，海冰強度逐漸進入“負坡段”，負阻尼使得風力機從海冰動能中汲取能量，導致塔基剪切力大幅增加。前后向剪切力在冰速進自激振動和出自激振動過程快速增加。側向剪切力隨冰速平緩增加，但出自激振動過程快速減小。進入中高冰速區后，塔基剪切力近似不變。此時，海冰強度位于脆性區，強度不再隨加載速率變化而變化。

圖6 塔基剪切力最大值Fig.6 The maxes of the shear forces on the tower base

標準差系數為標準差與均值之比無量綱系數，作為反映塔基各向剪切力波動劇烈程度的參考量。

塔基各向剪切力的標準差系數如圖7所示。由圖7可知，在冰速小于0.1 m/s和大于0.4 m/s時，塔基剪切力前后標準差系數小于側向，其余冰速前后向更大。塔基前后向剪切力標準差系數最大值約為側向8倍。塔基前后向剪切力均值約為620.7 kN遠大側向均值32.3 kN，故前后向波動幅度大于側向。在中低冰速區，標準差系數較大。前后向剪切力在冰速0.1 m/s后出現一個快速的增加，并達到最大值。側向剪切力在冰速0.1 m/s后出現一個平緩增加，逐漸達到最大值后再快速減少。前后向剪切力標準差系數隨冰厚增加而降低。

塔頂各向位移最大值如圖8所示。在中低冰速區塔頂位移明顯更大，當冰速超過0.4 m/s后，各向位移近似不變。在低冰速區時，塔頂前后向位移未出現小幅增加。塔頂前后向最大位移約為側向5.84倍。

圖7 塔基剪切力標準差系數Fig.7 Coefficient of variations of shear forces on the tower base

圖8塔頂位移最大值Fig.8 The maxes of the displacements on the tower top

5.2 塔頂位移頻域結果

對冰厚1.0 m下的冰速0.02 m/s、0.06 m/s、0.08 m/s、0.22 m/s、0.30 m/s和0.40 m/s，6種冰況造成的塔頂位移進行快速傅里葉變換(Fast Fourier Transform, FFT)，以研究自激冰振下的風力機結構動態響應頻率特性。

根據冰池試驗結果顯示當圓柱結構受過低和過高冰速作用，其頻域上均未能觀察到穩態振動，當冰速介于之間，在頻域上可以觀察到明顯頻率集中，視為出現穩態振動。

各冰速下塔頂前后向位移頻域結果如圖9所示。在塔架固有頻率0.31 Hz附近，塔頂前后位移幅值有所增加，且增幅隨冰速增加緩慢增加。冰速為0.40 m/s時幅值約為0.01 m，較冰速0.02 m/s增加約為0.002 m。當冰速0.22 m/s和0.30 m/s時，受自激冰振影響，塔頂前后位移頻率在1.1 Hz附近大幅增加，遠大于其他頻率下幅值，分別為0.168 m和0.319 m。此時，風力機發生可視為發生明顯的穩態振動。此時，內部應力變大，結構動態響應增加，更甚使用壽命將減少，事故發生概率激增。當冰速0.02 m/s和0.40 m/s時，塔頂前后位移在1.1 Hz下幅值增幅不明顯。

圖9 塔頂前后位移頻域結果Fig.9 The frequency domain results of the fore-aft displacement on the tower top

圖10 塔頂側向位移頻域結果Fig.10 The frequency domain results of the side-side displacement on the tower top

各冰速下塔頂側向位移頻域結果如圖10所示。側向位移幅值在1.1 Hz和固有頻率0.31 Hz附近明顯增加。塔架固有頻率處幅值隨冰速緩慢增加。冰速0.02 m/s時，幅值為0.008 1 m；當冰速0.40 m/s時，幅值為0.011 9 m。當冰速為0.22 m/s和0.30 m/s時，1.1 Hz處幅值大于其他冰速，分別為0.009 3 m和0.021 m。當冰速0.06 m/s和0.08 m/s時，1.1 Hz處幅值小于其他冰速。

6 結 論

隨著風力發電的發展，風力機不得不面臨越來越復雜的載荷環境。其中海冰載荷是安裝于高緯海域風力機將面臨載荷之一，具有作用力大、持續時間長、振動特性復雜的特點，會造成巨大的塔基剪切力，威脅塔架穩定，進一步影響的風力機安全運行。針對以上情況，采用M??tt?nen 自激冰振模型，以單樁作為支撐的NREL 5MW近海3葉片水平軸風力機為研究對象，研究不同冰況下海冰載荷對風力機影響，得出以下結論：

(1)受海冰強度“負坡段”作用，風力機在海冰載荷碰作用生自激振動，振動發生冰速主要集中在中低冰速(0.40 m/s以下)。

(2)塔基剪切力在低冰速(0.10 m/s以下)時有小幅增加之后隨冰速繼續增加逐漸進入自激振動，塔基剪切力大幅增加，增幅遠大于之前。當冰速超過0.40 m/s，塔基剪切力轉而減小至近似不變。

(3)自激振動時，塔基剪切力隨冰厚增加而減少。

(4)塔基剪切力前后向波動大于側向。

(5)受自激振動影響，塔頂位移顯著增加，前后向大于側向，前后向最大位移約為側向5.84倍。

(6)在海冰和湍流風共同作用下，塔頂位移在頻率上幅值主要集中在1.1 Hz和塔架固有頻率0.31 Hz附近。自激冰振(冰速0.22 m/s和0.30 m/s)造成塔頂前后向位移在1.1 Hz處的幅值急劇增加，數值遠大于其他頻率處幅值。此時，可視為塔架發生明顯的穩態的振動。

[ 1 ] 李春, 葉舟, 高偉, 等. 現代陸海風力機計算與仿真[M]. 上海: 上海科學技術出版社, 2012.

[ 2 ] 張豐豪, 何榕. 基于柔性多體動力學的風力發電機固有頻率研究[J]. 中國電機工程學報, 2014, 34(11): 1807-1814.

ZHANG Fenghao, HE Rong. Natural frequency analysis of wind turbine based on flexible mulit-body dynamics[J]. Proceedings of the CSEE, 2014, 34(11): 1807-1814.

[ 3 ] JONKMAN, JASON M, MARSHALL L, et al. FAST user’ guide: NREL/EL-500-38230[R]. National Renewable Energy Laboratory, 2005.

[ 4 ] GEBHARDT C G, ROCCIA B A. Non-linear aeroelasticity: an approach to compute the response of three-blade large-scale horizontal-axis wind turbines[J]. Rnewable Energy, 2014, 66(6): 495-514.

[ 5 ] ZHAO X, MAIBER P, WU J. A new multibody modelling methodology for wind turbine structures using a cardanic joint beam element[J]. Renewable Energy, 2007, 32(3): 532-546.

[ 6 ] MARTIN O L H. 風力機空氣動力學[M]. 肖勁松, 譯. 北京: 中國電力出版社, 2012.

[ 7 ] MICHAILIDES C, GAO Z, MOAN T. Experimental study of the functionality of a semisubmersible wind turbine combined with flap-type wave energy converters[J]. Renewable Energy, 2016, 93: 675-690.

[ 8 ] 楊陽, 李春, 繆維跑, 等. 湍流風場與地震激勵聯合作用下的風力機結構動力學響應[J]. 振動與沖擊, 2015, 34(21): 136-143.

YANG Yang, LI Chun, MIAO Weipao, et al. Response of structural dynamic charateristics of wind turbine operating in turbulent wind combined with Seismic Motion[J]. Journal of Vibration and Shock, 2015, 34(21): 136-143.

[ 9 ] 李新一, 祁文軍. 水平軸風力機葉片在風沙流場中的數值模擬[J]. 現代制造工程, 2014，1: 64-67.

LI Xinyi, QI Wenjun. Numerical simulation of horizontal wind turbine blad in the sand-flow field[J]. Modern Manufacturing Engineering, 2014,1: 64-67.

[10] 李靜, 陳健云, 陳小波, 等. 海上風力發電體系在風-冰作用下的動力反應[J]. 太陽能學報, 2013, 34(7)： 1234-1240.

LI Jing, CHEN Jianyun, CHEN Xiaobo, et al. Dynamic responses of analysis of wind turbine tower subjected to wind force and ice exction[J]. ActaEnergiae Solaris Sinica, 2013, 34(7)：1234-1240.

[11] 周國龍, 葉舟, 成欣, 等. 垂蕩板對傳統Spar平臺水動力特性的影響[J]. 水資源與水工程學報, 2015, 26(4): 143-148.

ZHOU Guolong, YE Zhou, CHEN Xin, et al. Influence of heave plate on hydrodynamic characteristics of traditional spar platform[J]. Journal of Water Resources & Water Engineering, 2015, 26(4): 143-148.

[12] 丁勤衛, 李春, 楊陽, 等. 極限海況下三種漂浮式風力機平臺的動態響應對比[J]. 水資源與水工程學報, 2015, 26(1): 159-165.

DING Qinwei, LI Chun, YANG Yang, et al. Comparison of dynamic response for three floating wind turbine platformsunder extreme sea situation[J]. Journal of Water Resources & Water Engineering, 2015, 26(1): 159-165.

[13] BARKER A, TIMCO G, GRAVESEN H, et al. Ice loading on Danish wind turbines Part 1: dynamic model tests[J]. Cold Regions Science & Technology, 2005, 41(1): 1-23.

[14] GRAVESEN H, S?RENSEN S L, V?LUND P, et al. Ice loading on Danish wind turbines Part 2: analyses of dynamic model test results[J]. Cold Regions Science & Technology, 2005, 41(1):25-47.

[15] 柳春圖, 段夢蘭. 海冰工程中的結構力學問題[J]. 機械強度, 1995, 17(3): 7-20.

LIU Chuntu, DUAN Menglan. Problem of structural mechanics in sea ice engineering[J]. Journal of Mechanical Strength, 1995, 17(3): 7-20.

[16] 黃焱. 冰激海洋平臺振動的動力模型試驗研究[D]. 天津: 天津大學, 2004.

[17] 劉延柱, 洪嘉振. 多剛體系統動力學[M]. 北京: 高等教育出版社, 1989.

[18] 徐繼祖. 海冰引起的結構振動[J]. 海洋工程, 1986, 4(2): 42-47.

XU Jizu. The structure vibration caused sea ice[J]. The Ocean Engineering, 1986, 4(2): 42-47.

[19] 張力. 導管架海洋平臺冰激振動控制的實驗研究[D]. 大連: 大連理工大學, 2008.

[20] PEYTON H R. Sea ice foreces: ice pressures against structures[C]∥ Ottawa:IPAS, 1966.

[21] NEILL C R. Dynamic ice forces on piers and piles. An assessment of design guidelines in the light of recent research[J]. Canadian Journal of Civil Engineering, 1976, 3(3):305-341.

[22] MARLOCK H, DAWKINS W P, PANAK J J. Analytical model for ice-structure interaction[J]. Journal of the Engineering Mechanics Division, 2014, 97(4):1083-1092.

[23] 王翎宇, 徐繼祖. M??tt?nen模型在多腿結構上應用[J]. 中國海上油氣(工程), 1992, 4(5): 43-49.

WANG Lingyu, XU Jizu. The M??tt?nen model application on multi-leg systems[J]. China Offshore Oil and Gas(Engineering), 1992, 4(5): 43-49.

[25] 徐繼祖, 王翎羽. 冰力振子模型及其數值解[J]. 天津大學學報, 1988(4): 28-34.

XU Jizu, WANG Lingyu. Ice force oscillator model and its numerical solutional[J]. Jouranl of Tianjin Unversity, 1988(4): 28-34.

[26] 韓雪松. 筒型基礎平臺冰激振動下地基的動力軟化研究[D].天津: 天津大學, 2004.

[27] 彭忠. 柱體剛度對冰激振動影響的實驗研究[D]. 天津: 天津大學, 2005.

[28] JONKMAN B J, KILCHER L. TurbSim user’guide[R].[S.l.]: National Renewable Energy Laboratory, 2012.