長期水平荷載對單樁式海上風機結構自振頻率的影響分析

孫毅龍， 許成順， 席仁強,2， 杜修力， 豆鵬飛,3

(1. 北京工業大學 城市與工程安全減災教育部重點實驗室，北京 100124； 2. 常州大學 機械工程學院， 江蘇 常州 213164； 3. 清華大學 水利水電工程系，北京 100084)

在海上風能的開發與利用中，風機結構整體自振頻率是海上風機設計中的主要難點之一。目前工程中均采用“軟-剛”模式進行海上風機結構和基礎的設計[1-2]，即海上風機結構系統自振頻率介于1P頻率(發電機轉動頻率)帶和3P頻率(葉片掃掠頻率)帶之間。但由于1P頻率帶和3P頻率帶之間范圍較小，因而對海上風機整體結構自振頻率的設計要求較高。

為準確地計算海上風機整體結構自振頻率，國內外學者開展了相關研究。Zaaijer[3]推導了簡化解析方法，研究了單樁式海上風機結構的振動特性及地基剛度對海上風機整體結構自振頻率的影響。Byrne[4]將地基與基礎、上部結構的相互作用簡化為一個旋轉自由度的彈簧，頂部簡化為集中質量，從而建立了海上風機結構整體自振頻率的簡化計算方法。Bhattacharya等[5-6]將基礎與地基間的相互作用簡化為水平剛度彈簧和搖擺剛度彈簧，上部結構等效為Euler-Bernoulli梁，從而建立結構自振頻率的求解方法，并與模型試驗、有限元計算進行了對比。黃茂松等[7]采用類似的樁土相互作用簡化方法，建立群樁基礎支撐的海上風機結構的頻率分析模型。上述學者建立的結構自振頻率計算方法形式簡單，但基礎剛度的確定較為復雜。因而Andersen等[8-10]將樁土相互作用簡化為溫克爾彈簧，建立了單樁式海上風機結構整體自振頻率計算方法。

然而海上風機結構所處的荷載環境特殊，在其運營周期內將長期承受風、波浪等循環荷載作用[11]，同時面臨風暴、地震等極端荷載的威脅，因而地基土物理力學特性會發生改變。Achmus等[12-14]基于室內土體循環三軸試驗，建立了考慮長期循環荷載效應的單樁數值模型，發現長期循環荷載作用會導致樁周土體參數發生變化。Nikitas等[15]采用動單剪儀，開展長期循環荷載下砂土模量、變形的變化規律研究，發現循環前期土體密實剪切模量增大，隨后逐漸趨于穩定，但該研究中未考慮樁基對土體參數變化的影響。Bayat等[16-17]基于建立二維飽和兩相介質計算模型，發現循環荷載作用會引起樁-土相互作用剛度的降低，但考慮的循環荷載加載次數較少。Cuéllar等[18-20]開展了地基土微觀機理變化研究，發現樁周土體反應在長期循環荷載作用前期表現為擠密，后期加載中樁周土體產生棘輪效應。在海上風機運營期間，樁周土體彈性模量等關鍵參數的較大改變，會對樁基承載特性、海上整體結構動力特性產生顯著的影響[21-22]。Lombardi等[23-24]開展了單樁式海上風機結構水平循環荷載模型試驗，驗證了風機結構自振頻率的偏移現象；最大偏移量達30%以上。Carswell等[25]基于p-y曲線模型，分析了黏土剛度變化對海上風機結構整體自振頻率的影響，發現黏土剛度的退化導致風機系統整體自振頻率向1P頻率偏移。

綜上所述，地基土物理力學特性的變化會導致風機結構系統自振頻率發生偏移現象，因而在海上風機結構自振頻率設計計算中，應考慮長期循環荷載對海上風機結構自振頻率的影響，但現有的研究卻鮮有報道。本文基于動力運動方程，考慮樁土相互作用，通過嵌入地基剛度衰減模型考慮長期循環荷載引起的地基剛度變化，建立單樁式海上風機結構自振頻率的簡化計算方法，探討循環荷載大小、加載次數對海上風機結構自振頻率的影響規律；為工程設計提供參考。

1 數值算法的建立

1.1 考慮樁土相互作用的自振頻率計算方法

單樁基礎是當前海上風電場中常用的基礎型式，因此本文計算方法的建立主要針對單樁基礎。由于風機塔筒為變截面，因此塔筒離散為若干段，樁土相互作用簡化為水平地基彈簧，風機、扇葉等簡化為集中質量，計算模型示意圖如圖1所示。

圖1 簡化模型示意圖Fig.1 Diagram of simplified model

根據結構動力學中，有阻尼體系的自由振動運動方程如下[26]

(1)

式中：M為質量陣；C為阻尼陣；K為剛度陣；u為計算方向的位移向量。

采用集中質量法，將海上風機上部結構化為多自由度體系，取樁基、結構等簡化為的節點質量mi構造質量矩陣M，如下

(2)

阻尼矩陣按照經典的Rayleigh阻尼構造為對角矩陣

(3)

因而在未考慮地基土阻尼時，阻尼矩陣如式(4)所示。其中ξ為鋼材料的阻尼系數，取為2%；本文主要研究海上風機結構系統的自振頻率，因而ω為結構的一階頻率。

C=2ξ?M

(4)

樁土相互作用簡化地基彈簧中地基阻尼取為[27]

(5)

式中：ξs為土體的阻尼比；ρs為土體的密度；Vs為土體的剪切波速；D為樁基直徑；a=ωD/Vs；ki為地基土剛度。

參考黃茂松等和楊春寶等的方法，構造剛度矩陣K為

(6)

(7)

(8)

(9)

樁土相互作用的水平地基剛度kxi選取p-y曲線的初始地基剛度，由于海上風機單樁基礎的樁徑較大，其尺寸效應顯著[28]，因而本文方法中的地基剛度參照作者在文獻[29]中建立的修正地基剛度，如式(10)所示。將kxi添加到剛度矩陣K地基部分的剪切剛度項中，形成考慮樁土相互作用的的整體剛度矩陣。

(10)

式中：nh在不同場地下的取值參照圖2；D0為相對樁徑，取1 m；zi為土體深度，m；z0為相對深度，取2.5 m；m，n為影響指數，取值參照文獻[30]。

圖2 地基模量常數nh與相對密實度的函數關系Fig.2 Constant of subgrade reaction nh versus relative density

將上述矩陣集成、組裝，代入式(1)中，得整體振動方程，先求解無阻尼下的自振頻率，再代入根據式(4)和式(5)計算有阻尼下的自振頻率。該計算方法與楊春寶等的自振頻率求解方法原理是一致的，不同的是地基剛度考慮了樁徑的增大對地基剛度的影響。

1.2 考慮砂土剛度衰減的自振頻率計算方法

海上風機單樁基礎在風、波浪等循環荷載作用下，樁周土體的模量等參數會發生變化，如圖3所示。為描述土體的該特性，Huurman[31]開展了土體三軸試驗，提出了砂土在長期循環荷載下的剛度衰減模型，砂土的割線模量與加載次數也為指數型函數關系，具體如下

圖3 循環荷載下砂土的割線剛度衰減示意圖Fig.3 Diagram of degradation of sand secant modulus under cyclic loading

(11)

(12)

式中：EsN為循環荷載加載N次下砂土的割線模量，kPa；Es1為循環荷載第一次下的土體割線模量(土體的初始彈性模量)；N為循環荷載加載次數；a，b為剛度衰減模型參數，其中，密實砂土，a為0.20，b為5.76，中密砂土，a為0.16，b為0.38[32]；X為砂土的循環應力比；σ1,cyc為土體循環周期的最大主應力，kPa；σ1,sf為土體靜力破壞時的最大主應力，kPa。

胡安峰等[33]通過引入循環折減系數λ(式(13))上述砂土模量的弱化效應嵌入到雙曲線p-y曲線中，通過與試驗結果、數值結果的對比，驗證了該循環弱化方法的正確性。

(13)

式中：kxiN為循環荷載作用N次后的水平地基剛度；t為循環荷載下的土體剛度弱化因子。

循環弱化因子t與土體的循環特性、應力狀態密切相關，其計算公式為

t=0.251a(X)b

(14)

(15)

(16)

(17)

式中：a，b的取值與式(11)是一致的；p1為靜力p-y曲線計算的樁周土反力；Kp為土體的被動土壓力系數；γ為土體有效重度；z為土體深度；D為樁徑。

根據式(13)～式(17)計算得到N次循環后的土體初始地基剛度kiniN，替換剛度矩陣K中的地基剛度部分，得考慮地基剛度弱化的整體剛度矩陣，代入式(1)，進行求解，得到單樁式海上風機結構的自振頻率f0和N次循環后的自振頻率fN，該方法計算流程圖，如圖4所示。

圖4 計算原理流程圖Fig.4 Flow chart of calculation principle

基于MATLAB數值計算軟件平臺，將上述理論公式，按照圖4的計算流程編寫計算程序，下面開展工程應用和對比研究，探討建立的算法的正確性。

2 工程應用與算法的驗證

本文計算方法以愛爾蘭Walney風電場單樁式3.6 MW海上風機為例，開展計算結果與現場監測結果、數值計算結果的對比研究，說明建立的簡化方法的有效性。

2.1 與工程監測結果的對比

該風機場地條件主要為密實砂土場地，砂土內摩擦角為40°，彈性模量為196 MPa，土體重度為21 kN·m-3[34]。該風機的具體參數如下：風機輪轂高度為83.5 m；本文研究對象為海上風機結構的自振頻率，因此將風機和扇葉等結構簡化為236 t的集中質量；單樁基礎樁徑為6 m，埋深為23.5 m；上述塔筒等結構密度為8 500 kg/m3，單樁密度為7 800 kg/m3，彈性模量均為210 GPa。具體尺寸參數如圖5所示。

圖5 Walney海上風機支撐結構尺寸Fig.5 Diagram of geometry for Walney offshore wind turbine support structures

同時根據規范推薦的模態分析法，基于FLAC3D有限差分平臺，建立單樁式海上風機結構整體數值計算模型，如圖6所示，地基大小為100 m×100 m×50 m，土體采用彈性本構模型模擬，本研究主要關注樁-土相互作用及場地土特性變化對結構頻率的影響，因而單樁、塔筒、輪轂等均采用彈性本構模擬，僅考慮了結構的阻尼比取2%，未考慮風機運行中的氣動阻尼等阻尼特性的影響。此外，依據圖2確定建立簡化方法的土體參數，保證與數值計算模型的土體參數是相互對應的。下面開展在考慮(SPI)與未考慮(No-SPI)樁-土相互作用時，本文建立的計算方法的計算結果與數值計算結果、現場監測結果的對比，如表1所示。

圖6 海上風機結構整體數值模型Fig.6 Overall numerical model for offshore wind turbine structure

表1 海上風機結構自振頻率計算結果的對比Tab.1 Comparison of natural frequency for offshore wind turbine structure 單位：Hz

由表1的結果對比可知，當考慮樁-土相互作用(SPI)時，本文建立的計算方法的計算結果與試驗結果、數值結果的誤差在1.5%以內，此外本文建立的方法計算結果較楊春寶等的方法更為接近測量值，同時當塔筒底部固定即未考慮樁-土相互作用(No-SPI)時，本文方法計算值與數值模擬計算結果差異在1%內。綜上所述，說明了本文建立的海上風機結構自振頻率計算方法的有效性。

2.2 考慮砂土剛度衰減的計算結果對比

下面采用數值模擬計算循環荷載下前述3.6 MW風機的自振頻率，并與建立的簡化計算方法進行對比，探討建立的簡化方法的合理性。

基于Flac3D內嵌的Fish語言和Achmus等的研究，進行二次開發，在Flac3D數值計算程序實現了砂土的剛度衰減模型，以反映循環荷載對場地土的物理力學特性的影響，其實現過程具體如下：

(18)

(3) 將Xc替換式(11)中的X，根據土體的初始彈性模量，可計算得到N次循環后的土體模量；基于Flac3D內置編程語言，將上述計算流程嵌入到數值模型中，可用于計算N次循環后大直徑單樁的受力特性。

采用該數值模型，與張紀蒙等[35]開展的模型試驗進行對比，樁基仍采用彈性本構模型，土體采用Mohr-Coulomb彈塑性本構模型，原型樁基樁徑為3 m，埋深30 m，土體相對密實度為53%，內摩擦角為33°，彈性模量取值采用式(19)計算。如圖7所示，數值計算結果與試驗結果基本吻合，說明該土體剛度衰減模型較好地反映了循環荷載下土體的循環弱化特性。

圖7 數值計算與試驗結果的對比Fig.7 Comparison of numerical calculation and results for tests

(19)

σ3=K0γz

(20)

式中：E0為土體的初始彈性模量，kPa；k′與λ由土體的土體密實程度決定，在此取為380和0.63；σ3為土體的小主應力，kPa；σat為參照應力取大氣壓力(105kPa)；z為土體深度，m；K0為土體的靜止土壓力系數。

將該土體剛度衰減模型應用于圖6中的數值計算模型中，采用自由振動衰減法計算自振頻率，即在塔頂施加一定的荷載，隨后將該荷載釋放，得塔頂結構的響應時程，然后將該時程進行快速傅里葉變換處理，得其自振頻率。

為保證兩種方法施加的荷載是一致的，根據式(21)對施加的荷載進行無量綱處理。在數值模型中計算極限水平荷載，然后根據式(21)計算荷載無量綱系數，上述施加的水平位移荷載的無量綱系數ξL為0.17。

(21)

式中：ξL為循環荷載大小的無量綱系數；Smax為循環位移荷載下樁基泥面處的循環水平位移幅值，m；Su為樁基泥面處的極限水平位移，m，大小為0.1D[36-37]。

基于編寫的MATLAB計算程序，計算該荷載下樁基不同深度處的樁基水平位移，根據樁基水平位移計算不同深度處的樁基土抗力p1，代入式(15)得到循環應力比X，從而計算N循環后的海上風機結構自振頻率。圖8是本文算法的計算結果與數值結果的對比，為了清晰地對比兩者的差異，橫坐標取加載次數的對數。結果表明，兩者的計算結果基本一致，均為隨著加載次數的增加，自振頻率減小，且兩者的數值差異在4%，說明了本文建立算法的可行性。

圖8 不同加載次數下風機整體自振頻率的對比Fig.8 Comparison of natural frequency under different cyclic number

3 參數影響分析

基于上述建立的海上風機結構自振頻率簡化計算方法，分析循環荷載大小、加載次數、樁徑對整體結構系統自振頻率的影響規律。

3.1 循環荷載大小的影響

樁基的累積變形發展，屬于樁基的疲勞極限狀態，所以式(21)中的荷載無量綱系數ξL最大約取0.30。采用本文建立的算法計算不同荷載幅值下的海上風機結構整體自振頻率，如圖9所示。由圖9可知，無論密實場地還是中密場地，隨著循環荷載幅值的增大，風機結構系統整體自振頻率線性減小，這是由于荷載的增大，樁附近土體的土抗力增大，使得特征循環應力增大，從而導致地基割線剛度衰減程度增大，即剛度矩陣中的地基剛度減小，因而導致風機結構系統自振周期延長。

圖9 不同加載次數下風機整體自振頻率的對比Fig.9 Comparison of natural frequency under different cyclic number

3.2 加載次數的影響

下面開展循環荷載加載次數對海上風機整體自振頻率的影響規律，計算了密實場地、中密場地兩類場地中，循環荷載的荷載無量綱系數為0.05～0.25時不同作用次數下風機結構系統的自振頻率，由于海上風機結構在其設計周期內承受108循環荷載的作用，因此本文計算中循環荷載次數取108，為了更好地分析循環荷載次數對風機結構整體自振頻率的影響規律，取前106的結果，如圖10所示。由圖10可知，循環荷載加載初期，隨著加載次數的增大，結構系統整體自振頻率變小；隨后當循環荷載次數N大于某一數值時，循環荷載作用次數對風機整體結構自振頻率影響較小。由式(11)可知，這是由于循環荷載作用次數的增多，樁周土體的割線模量減小，因而風機結構整體自振頻率變小。

圖10 不同加載次數下風機整體自振頻率Fig.10 Natural frequency of overall structures under different cyclic number

3.3 樁徑的影響

下面探討樁基樁徑對海上風機結構系統整體自振頻率的影響規律，如圖11所示。由圖11可知，隨著樁基直徑從2 m增大到10 m，中密和密實兩類場地下風機結構自振頻率分別從0.269 Hz增大到0.348 Hz、從0.286 Hz增大到0.353 Hz。結果表明，隨著樁基樁徑的增大，風機結構系統自振頻率對樁徑的敏感性減弱；在樁徑較小時，樁徑的變化對風機結構系統自振頻率影響顯著，但當樁基樁徑較大時(D>5 m)時，樁徑的變化對整體結構自振頻率影響較小，頻率變化在3%以內。因此在單樁式海上風機結構自振頻率設計時，應分析樁徑的影響，保證樁基直徑處于無影響范圍內。

圖11 不同樁徑下風機整體自振頻率Fig.11 Natural frequency of overall structure under different pile diameter

4 結 論

為研究長期循環荷載對海上風機結構自振頻率的影響，本文基于MATLAB數值分析計算平臺，采用結構的動力運動方程，在考慮長期循環荷載對地基剛度影響的基礎上，建立了單樁式海上風機結構的自振頻率的計算方法，通過與現場監測數據、數值計算結果對比，說明了建立的方法的可行性。最后探討了循環荷載大小、作用次數、樁徑對風機結構自振頻率的影響規律。得出如下結論：

(1) 在海上風機結構的疲勞荷載工況下，海上風機結構體系自振頻率隨著循環荷載幅值的增大而線性減小。

(2) 循環荷載作用初期，風機結構體系自振頻率隨著循環荷載次數的增大而顯著減小；此后其自振頻率減小趨勢變緩，最終基本穩定而不再隨加載次數的增加而變化。

(3) 單樁基礎樁徑較小時，對風機結構整體自振頻率影響顯著，存在一界值，當樁基直徑大于該數值時，樁基直徑的變化對風機結構整體自振頻率基本無影響。

(4) 考慮到長期循環荷載作用對風機結構自振頻率的降頻效應，為保證海上風機結構體系的長期運營安全，結構體系的設計自振頻率應稍偏移3P。

(5) 在缺乏風機結構體系自振頻率長期監測數據的情況下，本文提出的自振頻率簡化計算方法可以為單樁式海上風機結構的頻率設計提供一定的參考。

本文建立的算法考慮場地條件為均質中密砂土、密實砂土，對于其他密實程度的砂土場地、復雜多層場地、黏土場地下，海上風機結構自振頻率在長期循環荷載下的變化規律仍需進一步開展研究。