朱本瑞，孫 超，黃 焱

(1.天津大學 水利工程仿真與安全國家重點實驗室，天津 300072；2.路易斯安那州立大學 土木與環(huán)境工程學院,美國 70803)

隨著世界各國對能源需求的增加，海上風能作為一種清潔可再生資源備受關注。根據全球風能理事會(Global Wind Energy Council,GWEC)2018年數據統(tǒng)計，世界范圍內的風能發(fā)電量達到539 GW，并預測在未來5年內，以平均每年55 GW的裝機總量速度持續(xù)增長[1]。我國海上風電事業(yè)起步較晚，但發(fā)展迅速，截至2019年6月底，海上風電累計裝機容量已達403萬千瓦，預計2020年，實現并網裝機790萬千瓦[2]。目前，海上風電不斷向大型化發(fā)展，2019年，金風科技在江蘇大豐海上風電中，安裝了GW184-6.45 MW的海上風電機組，該機組搭載90 m葉片，是國內已投運機組中葉片直徑最大的海上風力發(fā)電機組，刷新了亞太地區(qū)投運機組的最大葉輪直徑記錄[3]。

不同于陸上風電結構，海上風電服役環(huán)境惡劣，除承受風、浪、流等海洋環(huán)境載荷，在冰區(qū)海域還受到海冰的威脅。相對于風浪流載荷，海冰是一種極其特殊的載荷，受其溫度、鹽分、冰厚、冰速以及與之作用結構的剛度、作用面積和形狀等參數的影響，具有復雜多變的失效模式和截然不同的作用機理。截至目前為止，學術界對冰與結構的相互作用機理仍未達成統(tǒng)一的認識[4]。文獻[5]研究表明，海冰載荷是冰區(qū)海上石油平臺設計的主要控制載荷，特別是其引起的冰激鎖頻振動現象(冰載頻率與結構基頻基本一致)，嚴重威脅結構的服役安全。位于我國遼東灣的JZ20-2平臺現場監(jiān)測表明，海冰與柔性直立樁結構相互作用易引發(fā)結構的強烈穩(wěn)態(tài)鎖頻振動[6]；JZ9-3海洋平臺的現場監(jiān)測數據亦觀察到了這一現象[7]；國外，美國阿拉斯加庫克灣的鉆井平臺[8]，加拿大大型沉箱式采油平臺[9]等都曾遭受到冰激鎖頻振動事件；而一系列的室內實驗亦證實了冰激鎖頻振動的存在[10-11]。相對海上石油平臺，海上風機結構，特別是被廣泛采用的單樁風機基礎，具有更柔的剛度，更易發(fā)生鎖頻振動現象[12]，從而致使風機塔筒發(fā)生劇烈的振動，大大降低風機的發(fā)電效率，甚至導致結構發(fā)生疲勞破壞，嚴重威脅結構的安全。黃焱等[13]率先對渤海單柱三樁式海上風機的冰激振動問題進行了研究，建立了基于綜合控制因子的冰振事件區(qū)劃及概率預判方法；葉柯華等[14]基于強迫振動原理，研究了風冰聯(lián)合作用下海上風機的振動響應。這些現存的文獻主要關注冰-結構相互作用的機理，而對發(fā)生鎖頻時，如何有效的避免或控制結構的振動研究較少，特別是對新興的海上風電結構。

在結構振動控制領域，目前常用的控制方法包括被動控制、半主動控制以及主動控制[15]。文獻[16]采用被動式調諧質量阻尼器(TMD)研究了冰區(qū)海洋石油平臺的振動控制問題；文獻[17]提出了半主動式冰錐結構用于進一步緩解錐體冰載引起的結構振動；許多學者利用不同阻尼器亦開展了海上風機結構的振動控制研究[18-22]，但分析載荷限于風、波和地震，而針對海冰載荷下的振動控制卻鮮有文獻報道。海冰運動主要取決于海流和海風，實現中，其運動方向往往與風載荷成一定夾角，在冰和風二者聯(lián)合作用下，風機塔筒在葉片旋轉平面內(面內方向)和平面外方向(面外方向)均會發(fā)生不同程度的振動，因此，需要對塔筒同時進行面內和面外振動控制。文獻[23]提出了一種三維擺式阻尼器，驗證了其相對于常規(guī)TMD具有明顯的優(yōu)勢，且更適用于在多種載荷非同向作用時的振動控制。為此，本文針對NREL 5 MW海上單樁風機，基于M??tt?nen-Blenkarn模型原理，采用ANSYS參數化設計語言，開發(fā)海上風機自激振動分析程序，分析NREL 5 MW風機在冰區(qū)海域的鎖頻振動問題，并將三維擺式阻尼器應用于該風機的冰激振動控制中，研究其對冰區(qū)海上風機振動響應的控制效果，以期為我國大型風機在冰區(qū)海域的投產與運行提供技術支持。

1 冰區(qū)風機動力平衡方程

位于冰區(qū)的海上風機，冬季海平面結冰時，主要遭受風、流和海冰載荷的聯(lián)合作用，其動力平衡方程可表示為

FT(t)+FI(t)+FC(t)

(1)

2 動冰載荷

2.1 自激振動原理

目前，針對海冰引起的結構動力問題，國內外學者提出了多種冰-結構相互作用模型。概括而言，可將其分為強迫振動和自激振動兩大類別[24]。其中，自激振動理論被廣泛應用于海洋工程結構冰激振動分析中[25]。該模型基于冰體連續(xù)破壞假設，充分考慮了海冰抗壓強度與應力速率的關系，如圖1所示。

圖1 海冰抗壓強度與應力速率關系Fig.1 Ice crushing stress vs.stress rate with different interaction regions

該曲線將冰與結構相互作用劃分為三個典型的區(qū)域：(1) 低冰速延性區(qū)，冰排發(fā)生延性擠壓破壞，結構響應為準靜態(tài)；(2) 中冰速延性-脆性轉換區(qū)，結構發(fā)生穩(wěn)態(tài)振動，響應最大；(3) 高冰速脆性區(qū)，冰排發(fā)生脆性擠壓破壞，結構響應為隨機振動。自激振動模型認為在延-脆轉換區(qū)內，由于海冰抗壓強度隨冰力加載速率的增大而降低，從而導致了負阻尼效應，當該負阻尼大于結構阻尼時，會導致結構振動趨于較大的穩(wěn)態(tài)振動，即發(fā)生鎖頻振動。

2.2 M??tt?nen-Blenkarn模型

基于實測數據，M??tt?nen-Blenkarn給出了海冰抗壓強度與應力速率的多項式表達式[25]，即

(2)

(3)

動冰載荷則可表示為海冰抗壓強度與作用面積的乘積，即

(4)

由于海冰載荷是風機結構響應速度的函數，因此式(1)為需要編程進行迭代計算。整個計算分析的流程圖，如圖2所示，其中氣動風載荷的生成詳見本文第3章。

圖2 程序求解流程圖Fig.2 The flow chart of programming

3 空氣動力載荷

不同于海洋平臺結構受力特點，海上風機依靠巨大的葉片將空氣動能轉化為機械能，因而，其風載荷所占的水平側向載荷比例較大，是其動力分析的重要載荷之一。

3.1 湍流風場模型

海上風機受到的風速可表示為平均風速和脈動風速之和。其中，任意高度z處的平均風速可采用對數或指數型式風剖面計算。采用對數風剖面型式時，則有

(5)

式中:v(z)為高度z處的平均風速，m/s；v0為參考高度H處的平均風速，m/s；H為風機輪轂處高度，m；z為距海平面的高度，m；z0為粗糙度長度，m，對于海平面，取0.03 m。

脈動風速則可采用IEC卡曼譜[26]進行描述，即

(6)

式中:Sv(f)為脈動風速的功率譜密度函數；f為頻率，Hz；LC為積分尺度參數；I為湍流強度；v為風速，m/s。

考慮風速的空間相干性，任意計算兩點i和j的互譜Sij可表示為

(7)

式中:Sii和Sjj分別為計算點i和j處的自功率譜；L為計算點i和j之間的距離，m；a0為相干性衰減系數；vH為輪轂處的平均風速，m/s；其余參數同上。根據IEC-1 1th推薦數值，本文取a0=12，LC=340.2 m。

根據式(5)～式(7)，基于TurbSim程序生成風機葉輪處的三維風場[27]，采用31×31的網格對整個風機掃掠面積的風場進行離散，然后利用MATLAB開發(fā)的程序得到作用于葉輪每段微元上的風速分布。

3.2 風機葉片載荷

將風機葉片沿展向離散為N個微元，假設每個微元之間互不影響，根據動量理論可推導作用于每段微元上的空氣推力和轉矩，采用葉素動量理論即可計算旋轉葉片上的空氣動力載荷[28]。

取任意第i個葉片微元分析，如圖3所示，相對來流風速vr可以表示軸向風速v(1-b)與切向風速Ωr(1+b′)的矢量和，即

(8)

(a) 葉片微元

式中:b和b′分別為軸向和切向速度誘導因子；r為微元距離輪轂旋轉中心的半徑，m；Ω為葉輪的轉速，rad/s。

于是，作用于每個葉素上的升力dFL和阻力dFD可以表示為

(9)

式中:c(r)為葉素弦長，m；CL和CD分別為升力系數和阻力系數；dr表示葉素長度，m；vr為相對于每個葉素的來流速度，m/s；ρ為空氣密度，kg/m3。

根據幾何關系，可以推導每個葉素上的軸向力和切向力為

dPN=dPLcos φ+dPDsin φ

dPT=dPLsin φ-dPDcos φ

(10)

式中:φ為相對來流風速vr與葉片旋轉平面的夾角，φ=a+θ，其中，a為風攻角；θ為翼型扭轉角。

由于軸向和切向速度誘導因子b和b′為非已知參數，采用式(8)～式(10)計算葉片風載荷時，需要通過迭代進行求解。為此，本文基于MATLAB開發(fā)了風機葉片軸向和切向風載荷時程的計算程序，并考慮了Prandtl葉尖損失修正和Grauert修正。圖4展示了風速為12 m/s，湍流強度為0.1時，計算得到的NREL 5 MW風機單個葉片上的風載荷時程。

圖4 葉片軸向風載荷與切向風載荷時程Fig.4 Normal and tangential aerodynamic loading applied on the blade

4 分析模型

4.1 5 MW風機參數

采用NREL 5 MW海上單樁風機作為分析對象，目標海域水深為20 m，樁基入泥深度36 m，塔筒總高90 m，過渡段距海平面高層為10 m，其基本尺度如圖5所示，其他參數詳見表1。將擺式阻尼器置于塔筒頂部，如圖5(a)所示。采用ANSYS軟件建立該風機振動響應分析模型，如圖5(b)所示，其中，風機葉片基于NREL 5 MW提供的剛度數據采用BEAM4單元進行等效模擬；機艙和輪轂等效為無密度剛性梁單元，其重量則采用MASS2進行模擬(分別為240 t和56.78 t)；塔筒采用BEAM188單元的變截面錐形梁模擬；單樁支撐結構采用PIPE59單元模擬；樁土相互作用則采用文獻[29]的方法，采用等效剛性梁和彈簧阻尼單元COMBIN14進行模擬，其中剛性梁等效長度取7.6 m，水平彈簧剛度為3.89×109N/m，轉動彈簧剛度為1.14×1011N·m/rad，轉動阻尼為9.34×108N·m·s-1/rad。

(a) 5 MW風機

(a) 面外位移

(a) vc=0.01 m/s

(a) 工況1 vc=0.06 m/s

(a) 面外方向

(a) 面外位移峰值

表1 5 MW海上風機基準模型參數Tab.1 Parameters of the NREL 5 MW baseline wind turbine

4.2 特征頻率

采用Block Lanczos法對該風機進行模態(tài)分析，求得其整體模型自振頻率，并與FAST軟件計算結果進行對比，如表2所示。

表2 模態(tài)對比結果Tab.2 Comparison between the finite element model and FAST

由表2可知，本文所建立模型的模態(tài)分析結果與FAST分析結果具有較好的一致性，但數值略小于FAST分析結果，其原因主要是邊界條件不同所致(文中模型考慮了土壤邊界約束，而FAST結果為固支約束條件下所得)。

4.3 擺式阻尼器參數

采用三維擺式阻尼器對該風機在海冰載荷下的振動響應進行控制，需要確定阻尼器的最優(yōu)參數，包括質量、最優(yōu)阻尼比以及最優(yōu)頻率比?；谖墨I[23]的研究成果，三維擺式阻尼器的最優(yōu)頻率比fo以及最優(yōu)阻尼比ζo的計算式為

fo=7.6μ2-2.5μ+1

ζo=-2.7μ2+μ+0.062

(11)

式中，μ為質量比，定義為阻尼器質量與風機結構總質量(不包含下部支撐結構和樁基部分)之比。

取質量比μ為2%，由式(11)求得fo=0.95，ζo=0.081；基于風機模態(tài)分析結果，求得擺式阻尼器的頻率fp為0.274 4 Hz，進而確定其擺長為3.27 m，阻尼為3 895.5 N·s/m，質量為13.9 t。在ANSYS有限元模型中，采用質量單元MASS21、LINK188以及COMBIN14單元分別模擬三維擺式阻尼器的質量、擺長和阻尼，即可建立帶阻尼器的分析模型，從而分析三維擺式阻尼器的減振效果。

4.4 分析工況

以我國渤海海域某風電場現場實測統(tǒng)計參數，對NREL 5 MW風機的鎖頻振動響應進行分析。該海域海冰抗壓強度為2.0 MPa，彎曲強度約為0.68 MPa，1年一遇平整冰厚為7 cm，50年一遇冰厚為32 cm?？紤]1年一遇海冰環(huán)境條件，對該冰厚下不同冰速以及不同環(huán)境載荷組合角度下的風機振動響應進行分析，以確定鎖頻發(fā)生的載荷工況。計算時，取風速為12 m/s，湍流強度為10%，風攻角保持垂直于風機葉片旋轉平面，即與x軸夾角為0°；海流剖面取線性剖面，表層流速取與冰速一致，底層流速為0 m/s，考慮4種不同的環(huán)境載荷組合方向，見表3。

表3 風冰流載荷組合方向Tab.3 Multi-direction combined of wind,ice and current load

5 結果分析

5.1 無阻尼器風機振動鎖頻分析

采用APDL開發(fā)基于M??tt?nen-Blenkarn模型的自激振動分析程序，取載荷步長為0.01 s，計算時長為300 s，針對表3中不同工況中等冰速進行搜索(冰速從0.01 m/s開始試算，每次增加0.01 m/s，直至冰載表現為脆性特征)，得到不同冰速下風機塔筒頂部的面內(x-z平面)和面外位移(y-z平面)響應，如圖6所示(限于篇幅，以工況2為例)。

由圖6可知，工況2下，即風攻角為0°，冰攻角為30°時，風機塔頂面外位移響應遠大于面內響應；當冰速為0.01 m/s，風機塔頂面內和面外位移振幅較小，表現為準靜態(tài)響應；隨著冰速的增加，當冰速為0.02 m/s、0.03 m/s、0.04 m/s，塔頂面外位移在50 s后表現為穩(wěn)態(tài)響應，塔頂面內位移則表現為典型的共振特性，即發(fā)生鎖頻振動，此時，風機塔頂振幅遠大于其他冰速下的響應，且隨著冰速的增加，其響應幅值逐漸增大，最大穩(wěn)態(tài)面外位移達0.93 m，最大面內位移達0.35 m，顯然，會影響風機的平穩(wěn)發(fā)電；當冰速相對較大時，即0.05 m/s和0.06 m/s時，隨著結構響應的衰減，海冰與風機結構相對速度增大，冰載荷降為較小值，塔頂響應最終趨于小幅振動，在穩(wěn)態(tài)響應階段受冰載荷的影響較小。

進一步去除瞬態(tài)分析初始效應的影響，取計算時間200 s后的風機塔筒頂部水平位移結果，對其進行傅里葉變換，分析其頻域特性，并與相應的冰載頻域特性進行對比，如圖7所示。

由圖7可知，當冰速為0.02 m/s、0.03 m/s、0.04 m/s時，風機塔筒振動主頻為0.28 Hz，略小于其結構一階頻率0.289 Hz(這與文獻[30]在實驗中發(fā)現的規(guī)律一致)，而此時，對應的冰力頻率亦在風機結構基本以及倍頻上出現峰值，即發(fā)生一階鎖頻振動；當冰速低于0.02 m/s時，冰載荷與塔筒位移的主頻率為0.27 Hz，小于其一階基頻；冰速為0.05 m/s，風機位移響應主頻為0.29 Hz，而冰載荷主頻率則演變?yōu)?.45 Hz；當冰速為0.06 m/s，風機振動主頻與風機葉片轉頻為0.2 Hz一致，即其振動響應主要受葉片運行及其風載荷控制，此時，冰載荷主頻比冰速為0.05 m/s時略有增大，變?yōu)?.48 Hz。綜合風機位移響應可知，工況2時，風機在冰速為0.02～0.04 m/s時，發(fā)生一階鎖頻振動，此時振動響應最大，而非鎖頻時，風機振動位移相對較小，運行平穩(wěn)，且在冰速大于0.05 m/s，冰載荷主頻跳躍為高頻，逐漸遠離風機基頻，向脆性區(qū)演變，冰載荷演化為風機振動的次要控制載荷。

表4給了其他工況下風機塔筒在計算時間200～300 s內的水平合位移響應最大幅值，以及發(fā)生鎖頻事件所對應的冰速。由表4可知，當海冰與風攻角相同時，發(fā)生一階鎖頻事件對應的冰速范圍最大，為0.02 m/s～0.06 m/s，且由于兩種載荷的組合疊加效應，風機位移響應振動最為劇烈，在冰速為0.06 m/s時，風機位移幅值達到1.424 m，嚴重威脅風機運行安全；隨著冰-風載荷夾角的增大，鎖頻事件對應的冰速范圍有所降低，風機位移響應幅值總體亦呈降低趨勢，當海冰方向為90°時，風機振動響應最小，發(fā)生鎖頻時的最大幅值位移降為0.76 m。

表4 風機塔筒水平位移振動幅值Tab.4 Horizontal displacement amplitude at tower top

5.2 三維擺式阻尼器振動控制評價

以4種環(huán)境組合工況中風機振動響應最大時的冰速為例，展示三維擺式阻尼器對風機冰激振動響應的控制效果，如圖8所示。

由圖8可知，工況1時，即冰-風同向時，三維擺式阻尼器對風機塔筒面外振動響應起到極佳的控制效果，穩(wěn)態(tài)響應的振動幅值減小可達90%；但在面內方向，僅在初始振動階段存在一定控制作用，當計算時間大于100 s后，部分時刻存在緩解振動的作用，部分時刻使風機振動略有增加，但應注意到，由于該工況風機面內載荷主要源自葉片的風載，因此，其面內位移振動很小，最大值不足0.2 m。工況4時，即當冰-風夾角為90°，三維擺式阻尼器在面內(冰載作用方向)起到良好的減振效果，在面外方向，則使得風機振動快速衰減至振動微小的穩(wěn)定響應；當冰-風夾角為30°和90°時，擺式阻尼器在起到雙向控制作用，在面內和面外均大大降低了風機塔筒的振動響應。

圖9展示了工況2下，冰速為0.04 m/s時有無3D-PTMD的風機塔頂面外和面內方向振動位移的頻率響應圖。由圖9可知，安裝3D-PTMD后，風機塔頂面外和面內的峰值頻率能量均被大大降低，3D-PTMD具有極佳的耗能效果；尤其在面內方向，3D-PTMD使得風機塔頂的面外峰值頻率由0.28 Hz調諧至0.2 Hz，減振動效果非常顯著。

計算表4中所有工況有、無3D-PTMD的風機塔筒響應峰值及標準差，進一步量化三維擺式阻尼器對風機振動控制的效果，如圖10所示。

由圖10可知，三維擺式阻尼器大大降低了風機塔筒頂部位移的峰值和標準差，面外方向最大分別降低為51.8%和93.7%，面內方向最大分別降低62.8%和84.2%，對風機塔筒的振動響應起到極佳的控制效果。由圖10(b)和10(d)可知，在面內振動風向，三維擺式阻尼器對部分計算工況存在負面影響，但這些工況主要為非鎖頻工況，此時的風機振動響應均相對較小(如工況2中冰速為0.06 m/s時，位移峰值和標準差分別放大6.4%和13.5%，而由圖6可知，該工況風機穩(wěn)態(tài)響應時的面外和面內位移振幅僅為0.141 m和0.096 m)，即三維擺式阻尼器的負面影響并不會干擾風機的正常運行。

6 結 論

基于自激振動原理，開發(fā)冰區(qū)風機振動分析程序，實現對NREL 5 MW海上單樁風機鎖頻振動的有效模擬，并對其振動響應進行控制，結果表明：

(1) 針對1年一遇冰情，綜合不同環(huán)境組合工況結果，該風機發(fā)生冰激鎖頻現象的冰速范圍為0.01～0.06 m/s；海冰作用方向對鎖頻事件的發(fā)生具有明顯的影響，當冰-風載荷同向時，對應的鎖頻冰速范圍最大，更易發(fā)生鎖頻事件，當冰-風非同向時，鎖頻所對應的冰速范圍變窄。

(2) 發(fā)生鎖頻時，隨著冰速的增加，風機塔筒振動響應逐漸增大，冰-風同向時，最大面水平振動幅值達1.424 m，即在常遇海冰作用下，一旦發(fā)生鎖頻振動，風機結構的運行安全性將受到嚴重威脅；冰-風非同向時，塔筒振幅降低，當冰-風為90°夾角時，振動響應最小。

(3) 三維擺式阻尼器可以對風機鎖頻振動響應起到極佳的控制效果，增加三維擺式阻尼器后，風機塔筒頂部面外和面內振動響應峰值和標準差均得到有效的抑制，在部分非鎖頻工況，三維擺式阻尼器存在負面效應，但影響甚微，綜合而言，在進行冰區(qū)海上大型風機安裝時，可以增設三維擺式阻尼器，以實現對風機振動響應的控制，提高風機運行的疲勞壽命。