風、浪作用下海上風機單樁結構時域動力響應分析

王 銜，邱 松，陳 濤，顧祥林

（1.中交第三航務工程局有限公司，上海 200032；2.同濟大學工程結構服役性能演化與控制教育部重點實驗室，上海 200092）

0 引 言

風能作為一種清潔能源，是21世紀用來降低對化石燃料的依賴性的重要資源。我國具有豐富的近海風力資源，國家在《能源發展“十三五”規劃》［1］中指出，要積極開發海上風電，推動海上風電技術的進步和商業化運營，發展7～10 MW級風電機組。近年來，隨著裝機容量的不斷增大，風機葉輪的半徑和重量持續增加，對基礎支撐結構的安全性提出了更高的要求。海上風電支撐結構在使用壽命周期內需承受高達約109次風、浪荷載的作用［2］，因此保證其疲勞性能滿足要求尤為重要。

單樁支撐結構最早由荷蘭萊利（NL）公司于1994年提出，截至2020年底，采用單樁支撐結構形式的海上風電基礎占整個歐洲所有已裝機的海上風電基礎的81.2%［3-4］。這種適用于淺水海域且工藝相對成熟的支撐結構也是我國海上風電事業未來發展的首選。單樁與上部結構連接的主要方式是灌漿連接，其技術原理是通過向內外鋼管間的環形間隙內填充高性能灌漿料，連接直徑不同的過渡段與鋼管樁。由于結構具有細長性，單樁支撐結構中的灌漿連接段主要承受反復彎矩的作用。本文結合廣東某地的場地條件和實測風、浪數據，對某5 MW風機單樁支撐結構所受風、浪荷載進行模擬，并進行時域動力響應分析，從而提取出灌漿連接段的荷載邊界條件時程，為進一步分析海上風機灌漿連接段的疲勞性能奠定基礎［5］。

1 某5 MW風機單樁支撐結構動力響應模型的建立

某5 MW海上風機單樁支撐結構的設計參數見表1，結構示意圖及實際場地條件見圖1。

圖1 某5 MW風機單樁支撐結構及實際場地條件

表1 某5 MW海上風機單樁結構的設計參數

1.1 樁-土相互作用的等效非線性彈簧

在將海上風機單樁支撐結構的單樁打入海床時，需考慮樁體與土壤（簡稱“樁-土”）的相互作用，尤其是樁體側向位移下的土體對單樁的作用。挪威船級社（Det Norske Veritas，DNV）報告［6］指出，若單樁底部簡單地按固接考慮，最多可能使結構基頻產生20%的誤差。

對于樁-土的相互作用，文獻［7］提出采用非線性彈簧矩陣進行等效。圖2為樁-土相互作用等效彈簧示意圖，在樁頂施加水平力Fh和水平彎矩Mr之后，樁頂產生了水平方向的位移uh和轉角θr，施加的荷載與位移之間的關系可用非線性等效彈簧矩陣表示，即

圖2 樁-土相互作用等效彈簧示意圖

由于對稱性，本文采用ABAQUS軟件建立樁-土相互作用區域的有限元半模型，見圖3。土體邊界半徑為單樁直徑的10倍，模型總厚度在樁體長度的基礎上向下延伸4倍單樁直徑，由此降低非無限大土體邊界條件帶來的影響。考慮到模型中土體的楊氏模量遠小于單樁鋼材，變形以土體為主，鋼材的材料采用理想彈塑性材料即可。土體可根據其類別（砂土或黏土）的不同采用不同的材料模型。對于砂土，根據摩爾-庫倫破壞準則，在ABAQUS軟件中需輸入其楊氏模量、泊松比、摩擦角、剪脹角、黏結強度和對應的絕對應變等參數。對于黏土，根據Tresca破壞準則，在ABAQUS軟件中需輸入其楊氏模量、泊松比、屈服強度和對應的塑性應變等參數。

圖3 樁-土相互作用有限元半模型

土體和樁體都采用8節點實體減縮積分單元（C3D8R）模擬。土體的環向單元均勻分布，徑向單元中間密兩側疏。上下層土體的環向網格和徑向網格都應對齊，厚度方向的單元大小也應保持一致，但超過樁體長度的下部土體網格在厚度方向上可適當稀疏。單樁在厚度方向上采用四層單元，樁管與土體環向的網格節點應盡量對齊，單樁內外側土體環向的單元節點也應盡量對齊。不同層土體之間的接觸采用ABAQUS軟件中的“tie”連接完全耦合。砂土和黏土與單樁接觸面之間的正向接觸都采用“硬接觸”，并允許接觸面分離；切向接觸雖然都遵照庫倫摩擦準則，但砂土與樁體的摩擦因數為tan（2／3φ′），其中φ′為砂土的有效摩擦角，且不設置摩擦強度的上限值。同時，黏土與樁體的摩擦因數一般取為常數0.4，并設置摩擦強度的上限值為黏土的不排水剪切強度。

對于土體的邊界條件：底層土體底面采用固接約束；圓弧面上約束2個水平方向位移和沿豎直軸的轉角；所有土層和單樁的對稱面都施加對稱面約束。同時，在樁頂截面形心位置處設置參考點，將該參考點與樁頂截面耦合。先進行重力荷載加載，再分別將水平力或水平彎矩荷載加載到參考點上，由此即可得到非線性等效彈簧矩陣，如式（1）所示。

1.2 某5 MW風機單樁支撐結構動力響應模型

海上風機單樁支撐結構屬于細長結構，鋼管的徑厚比一般在1／100左右，因此很多學者都采用梁單元［8-9］或ANSYS有限元計算軟件中的Pipe單元［10-11］模擬風機結構。由于不關注葉片的振動，很多文獻［12-15］都采用塔頂集中質量的方式模擬風機葉片和機艙。本文采用ABAQUS軟件中的梁單元和頂部集中質量的方式對風機單樁支撐結構進行簡化。

此外，目前大多數學者在模擬風機單樁支撐結構時都未考慮對灌漿連接段進行模擬。然而，灌漿連接段處的局部剛度發生變化會對單樁支撐結構的動力特性產生影響。本文采用ILIOPOULOS等［16］提出的按比例分配法將灌漿連接段擬合成一種勻質材料的方法，根據3種材料的體積占比對密度、彈模和泊松比進行加權平均，采用均一化材料梁單元進行模擬。在考慮第1.1節所述樁-土相互作用等效彈簧約束之后，得到簡化的單樁支撐結構動力響應模型見圖4。對該模型和不考慮樁-土相互作用而直接固接的其他同類模型進行模態分析，得到標準振型和振型頻率分別見圖5和表2。通過比較可知，忽略樁-土的相互作用會對單樁結構的基頻產生約6.3%的誤差，會對結構的二階頻率產生約35%的誤差，由此可見，考慮樁-土相互作用是很重要的。

圖4 風機單樁支撐結構動力響應模型

圖5 風機單樁支撐結構標準化振型模態

表2 某5 MW海上風機單樁結構振型頻率

文獻［17］給出了海上風機結構中阻尼的來源，包括材料本身的阻尼、土體阻尼、水動力阻尼和風機運行的氣動阻尼。氣動阻尼作為對風機結構的影響最大的“阻尼作用”，其大小由風速和風機翼型等因素決定。考慮到該阻尼的復雜性，本文參照文獻［18］和文獻［19］將靜止的單樁支撐風機結構的阻尼比取為1%，參照文獻［17］、文獻［18］和文獻［20］將考慮氣動阻尼的單樁支撐風機結構的總阻尼比取為5%。

2 單樁支撐結構風、浪荷載模擬

處在海洋環境中的風機單樁支撐結構會受到腐蝕、閃電、船舶撞擊、地震、浮冰和洋流等因素的影響，但文獻［21］指出，單樁風機受到的荷載主要來自風對塔身的作用、風機的氣動作用和波浪對單樁的作用，因此本文只考慮這3種作用。同時，假定風荷載與浪荷載作用在同一方向上，且風對塔身的作用和波浪對基礎結構的作用在高度10 m以下的分段內均勻分布，積分之后集中作用在10 m分段的中點處，見圖6。

圖6 風機單樁支撐結構風、浪荷載簡化示意圖

2.1 廣東省某地測風數據及風速、波高聯合分布

下面給出廣東省某沿海區域的測風塔和水文氣象站測量得到的實際風速、海浪有效波高Hs和有效周期Ts等數據。測風塔風速測量時段為2016年1月1日00：00—2016年12月31日23：00；海洋站波浪測量時段為2016年3月1日—2017年2月28日。

由測風塔得到的年平均空氣密度為1.175 kg／m3；下文風廓線理論式中的粗糙區長度z0＝6.02×10－6m。風機輪轂高度90 m處平均風速的長期分布經過擬合，符合雙參數的Weibull分布，其概率密度函數為

式（2）中：β為尺度參數，其值為8.9 m／s；α為形狀參數，其值為2.67。

結合海洋站的觀測數據，實際風、浪的年觀測數據共對應385種工況，其中：風速小于風機切入風速（3 m／s）的工況共有37種，概率總和約為5.32%；風速超過3 m／s的工況共有348種，概率總和約為94.50%，在這些工況下需計算風機運行產生的氣動荷載。

2.2 風對塔身的荷載模擬

風場具有隨機性，經典理論［22］指出，任意一點的風速都可視為平均風速項與脈動風速項的疊加。文獻［6］認為，在10 min時間內，脈動風速服從高斯分布，且其紊流強度可認為是不變的。但是，長期來看，10 min時間內的平均風速服從雙參數的Weibull分布。由于海上風機單樁支撐結構具有細長性，可將脈動風場簡化為沿高度方向不同位置的單方向一維、多元平穩高斯過程。因此，風機塔身高度方向上任意一點水平方向的瞬時風速可表示為

式（3）中：vs（z）為高度z處水平方向的平均風速；vD（z，t）為高度z處t時刻水平方向的脈動風速。

對于水平方向的平均風速，文獻［22］給出了不同高度處風速的風廓線理論，各高度處的平均速度可表示為

式（4）中：k為馮·卡門常數，常取為0.4［23］；z0為粗糙區長度；u*為剪切速度。

對于脈動風速的模擬，DNV規范［24］建議采用Kaimal譜模擬脈動風速，其表達式［22］為

式（5）和式（6）中：f為Monin相似坐標參數；ω為風的圓頻率。

對于不同高度處脈動風速之間的互功率譜，可采用Davenport互譜［25］，其表達式為

式（7）和式（8）中：CZ＝10；CY＝16［22］。由于忽略了水平方向的維度，本文無需考慮（x1－x2）項。

對于塔身微段面積dA內塔身的風荷載dF（z，t），同樣由平均風荷載項dFs（z）和脈動項dFD（z，t）組成，即

式（10）中：CW為風荷載形狀系數；ρ為空氣密度；vs（z）為某點水平方向的平均風速。

對于塔身的風荷載的脈動項dFD（z，t），文獻［26］將其功率譜寫為

式（11）中：Svkvl為k點和l點的風功率譜矩陣；vsk和vsl分別為k點和l點的平均風速；?j（k）和?j（l）分別為第j階單位振型下k點和l點的振幅；n為需考慮的振型階數，對于海上風機單樁支撐結構，一般只需考慮第一階彎曲振型即可滿足精度要求。最后，結合風對塔身的荷載的功率譜矩陣SFkFl和文獻［27］提出的改進的多元平穩隨機過程譜表達法即可模擬得到風對塔身的荷載時程，見圖7。采用相同的方法也可模擬風速沿塔身的時程，對模擬得到的風速時程進行傅里葉變換，將得到的模擬譜與目標Kaimal風譜相對比（見圖8），由此即可證明模擬程序的正確性。

圖7 沿塔身不同高度處的風荷載模擬結果

圖8 Kaimal風譜模擬結果對比

2.3 風機運行氣動荷載——葉素動量理論

目前可用來計算風機運行氣動荷載的理論主要有葉素動量理論、渦旋理論、廣義制動盤理論和納維-斯托克斯解法等，其中葉素動量理論較為簡易，不依賴專業的分析軟件，應用最為廣泛。結合經典文獻［28］，對葉素動量理論進行相關編程，在此不再贅述。編程中使用的5 MW標準風機是美國國家可再生能源實驗室（National Renewable Energy Laboratory，NREL）在2009年的報告［19］中提出的，其詳細參數可參見該報告。

風機推力荷載時程模擬結果見圖9。為驗證程序的正確性，采用文獻［19］中給出的標準風機扭矩與風速的關系曲線同模擬結果相對比，結果見圖10。從圖10中可看出，模擬結果良好，驗證了葉素動量理論程序的正確性。

圖9 風機推力荷載時程模擬結果

圖10 風機扭矩與風速關系模擬結果

2.4 波浪荷載模擬

本文采用B-M海浪譜對海面波面高程η（x，t）進行模擬，該波浪譜由Bretschneider于1959年提出，后由日本學者Mitsuyasu（恒光易）［29］進行了改進，其表達式為

式（12）中：Hs和Ts分別為海浪的有效波高和有效周期。

采用一維平穩隨機過程的譜表達法可得到海面波面高程的時程η（x，t），采用線性波浪理論可建立η（x，t）與任意深度z處波浪內的水質點的速度水平分量vx和加速度水平分量ax的關系，即

式（13）和式（14）中：t為時間；k為波數；ω為圓頻率；g為重力加速度；d為水深。

波浪對單樁的作用力可采用Morison波浪力方程［30］表示，該方程是由Morison于1950年提出的半經驗半理論公式，認為波浪對柱體的作用主要是黏滯效應和附加質量效應，因此可將波浪力分為同速度的平方成正比的阻力項和同加速度成正比的慣性力項，單位長度圓柱體上的波浪力F為

式（15）中：vx和ax分別為水質點的水平速度分量和水平加速度分量；ρw為波浪的密度；Dp為圓柱體直徑；CD為拖曳力系數；CM為慣性力系數。德國勞氏船級社規范［31］附錄4.H給出了CD和CM的取值，當波浪的雷諾數小于2×105時，對于光滑的圓柱體，可取CD＝1.2，CM＝1.0。

塔身不同高度處的波浪荷載模擬結果見圖11。對模擬得到的海面波面高程時程進行傅里葉變換，得到的模擬譜與目標B-M海浪譜對比見圖12，由此可證明模擬程序的正確性。

圖11 塔身不同高度處的波浪荷載模擬結果

圖12 B-M波浪譜模擬結果對比

3 單樁支撐結構時域動力響應分析

對第2.2～2.4節所述風、浪荷載模擬程序進行匯總，并將所有工況信息寫入Excel文件中，采用MATLAB軟件自動讀取Excel文件的功能依次讀取所有工況下的平均風速、有效波高、有效周期和發生概率等參數，自動生成風、浪荷載時程匯總文件，見表3。本文在模擬風、浪荷載時，取時間間隔為0.2 s，模擬總時長為10 min。當風速低于3 m／s時（即2.1節中的37種風機不工作工況），匯總文件中沒有氣動荷載一列。

表3 模擬得到的風、浪荷載時程匯總

在不同工況下進行時域動力響應分析時只有輸入的荷載不同，而單樁支撐結構的動力響應模型的阻尼比只在風機不運行（37種工況）和風機運行（348種工況）2種狀態下不同。為在ABAQUS軟件中實現批量計算，運用Python語言編寫相應程序，自動讀取表3中的數據，并將其修改寫入“inp”文件中。對于各工況，按其風速判斷風機是否處于工作狀態，修改結構的阻尼系數。

將批量修改的“inp”文件批量提交給ABAQUS Command，等待逐個自動計算，在酷睿i5處理器、內存24G的硬件條件下，每種工況的平均計算耗時約為12 000 s。

然而，通過對提取的灌漿連接段荷載邊界條件進行整理，發現一些問題。本文采用第一步施加重力荷載，第二步施加風、浪隨機動力荷載的加載方式。當第二步開始時，風浪荷載相當于突加荷載，結構在該突加荷載下會產生振蕩的過程，結構的位移和內力出現突增現象，且單樁支撐結構的總阻尼比較小，該振蕩并不能很快衰減。該振蕩過程可通過風機運行和不運行2種工況下的灌漿連接段頂部軸力時程曲線（見圖13和圖14）。圖13給出了風機不運行工況下灌漿連接段頂部的軸力時程曲線，此時單樁風機結構的阻尼比為0.01。從圖13中可看出，振蕩過程約持續300 s，在300 s之后軸力趨于穩定。圖14給出了風機運行工況下灌漿連接段頂部的軸力時程曲線，此時單樁風機結構的阻尼比為0.05。從圖14中可看出，振蕩過程約持續120 s，在120 s之后軸力和水平彎矩趨于穩定。因此，在后續對灌漿連接段的疲勞性能進行評估時，采用的荷載邊界條件數據在風機不運行工況下取時程分析的后300 s范疇，在風機運行工況下取時程分析的后480 s范疇。

圖13 風機不運行工況下的灌漿連接段頂部軸力時程曲線

圖14 風機運行工況下的灌漿連接段頂部軸力時程曲線

此外，編寫相應的Python程序，實現對所有工況下的計算結果文件（灌漿連接段荷載邊界條件）的自動提取，為進一步對灌漿連接段進行疲勞分析奠定基礎。

4 結 語

本文結合廣東某地的場地條件和實測風、浪數據，綜合運用Kaimal與Davenport風譜、風機葉素動量理論、B-M波浪譜、線性波浪理論和Morrison波浪力理論，對某5 MW單樁風機支撐結構所受風、浪荷載進行了模擬，并對單樁結構進行了時域動力響應分析，從而提取出灌漿連接段的荷載邊界條件時程，為進一步對灌漿連接段的疲勞性能進行分析奠定基礎。單樁支撐結構動力響應模型綜合考慮了樁-土相互作用的等效彈簧和灌漿連接段部分剛度的勻質化。研究發現，忽略樁-土的相互作用會對單樁結構的基頻產生約6.3%的誤差。同時，考慮風機氣動阻尼的影響之后，風機在運行工況下的動力響應進入穩態階段的時長相比在停機工況下更短。