基于高性能樁單元法的海上風機單樁支撐結構水平承載力研究

張舒燁，劉思威，馬會環

(1.中山大學土木工程學院，廣東，珠海 519082；2.香港理工大學土木及環境工程學系，香港 999077)

1 研究背景

海上風電機由機艙葉輪、下部結構和基礎組成[1]。其支撐結構往往采用大直徑單樁基礎，約占已投入使用基礎形式的80%[2]。單樁支撐結構下部基礎插入土體深部以獲得足夠的水平抗力，上部塔筒根據彎矩沿其高度的變化規律常常設計為錐形。然而由于海上風電機尺寸巨大，給支撐結構的設計施工帶來了巨大的挑戰，從而極大地提高了支撐結構成本，其成本已占到總成本的20%～30%[3]。

數值分析是單樁支撐結構設計中的關鍵一環，單樁支撐結構的主要數值分析主要基于實體單元模型和基于離散彈簧法的梁單元模型[4]，二者的原理如圖1 所示，在實際工程設計中，單樁基礎常常被簡化為文克勒地基梁模型，使用連接在梁上的一組剛度不同的彈簧(p-y曲線)表示樁土相互作用，因此這種方法被稱為p-y曲線法，通過在單元節點上連接非線性彈簧來定義樁土相互作用的方法被稱為離散彈簧法，在選擇合適的p-y曲線以及劃分足夠多單元的前提下，通過該方法得到的樁基水平響應將與實體單元法具有較好的一致性[5]。

基于地基梁模型，許多學者提出了許多不同形式并適用于不同土質的p-y曲線。在砂土方面，MCCLELLAND 和FOCHT[6]首次提出了砂土的p-y曲線，REESE 等[7]提出的直曲線形式的柔性樁砂土p-y曲線模型被應用于早期的API 規范[8]中，此后，諸多學者提出了包括雙折線[9]、拋物線[10]、雙曲線[11]等形式在內的砂土p-y曲線。在粘土方面，MATLOCK[12]、REESE 和WELCH[13]分別提出了適用于硬粘土與軟粘土的p-y曲線。SULLIVAN等[14]則在此基礎上提出了一種統一化方法，以適應于任何類型的粘土。河海大學的王惠初等[15]也基于鄧肯模型提出了一種新的粘土p-y曲線統一法。

在諸如海上風電樁基等大直徑樁的工程設計中，學界發現。傳統p-y曲線逐漸顯現出了不適用性[16]，并對適用于大直徑樁的p-y曲線開展了研究。YAN 和BYRNE[17]、KIM 等[11]指出，API 規范中推薦的p-y曲線會顯著高估土體的初始剛度，使得設計偏于不安全。針對上述問題，一些研究者對傳統p-y曲線進行了修正，以使其適用于大直徑樁情況[18-20]，也有部分學者認為，基于靜力觸探(CPT)的p-y曲線能夠適用于大直徑樁計算，這是由于相關文獻[21 - 22]指出，CPT 中的的土層錐尖阻力qc主要取決于水平方向有效應力，認為可以通過土層錐尖阻力qc來反過來研究水平土抗力的特性，因此周全智等[23]在數值模擬中引入了基于CPT 的p-y曲線以考慮樁基大直徑的影響。

現有p-y曲線研究結果大大便利了海上風機單樁支撐結構的設計計算，但海洋土體往往具有分層、參數隨土體深度變化呈現高度非線性的特點，因此若采用傳統的離散彈簧法，需要布置足夠多的離散土彈簧以表征上述非線性，這增大了有限元模型需要劃分的單元數量。而同時，對于風機塔筒的錐形結構，常規設計中常考慮使用多段恒截面梁對變截面梁進行等效[24]，這也增大了單元劃分量，進一步降低了計算效率。

針對上述問題，LI 等[25]提出了一種考慮土-結構相互作用的歐拉-伯努利梁柱單元，該單元通過在能量方程和切線剛度矩陣中引入樁土作用項實現樁土相互作用關系的內置。在計算時，該單元可以不通過節點施加彈簧，所以無需劃分過多單元。并且對于錐形梁，該單元能夠考慮截面剛度隨梁長增加的變化，無需通過多個均勻截面單元進行等效。以上二者的結合將大大降低在同時考慮錐形梁和樁土作用情況下的單元劃分量，進而提高計算效率，因此該單元在后文中以高性能樁單元代稱。高性能樁單元與傳統單元的比較如圖2所示。

本工作基于高性能變截面梁柱單元方法，結合SURYASENTANA 和LEHANE[26]、GUO 等[27]提出的大直徑修正p-y曲線，研究了高性能變截面梁柱單元方法在海上風機單樁支撐結構計算中的適用性，并進一步研究了大直徑單樁基礎的水平承載性能，以期為海上風機支撐結構設計計算提供參考。

2 高性能樁單元分析方法

2.1 基本假設

高性能樁單元分析方法通過以下假定簡化單元推導過程及表達式：1) 忽略剪切變形與翹曲變形；2)荷載保守；3)樁身材料均勻且具有各向同性；4)樁單元符合歐拉-伯努利假設。

2.2 單元形函數

已知單元節點的位移如圖3 左側所示，選取單元中點作為坐標系原點，以單元軸向為x軸建立正交坐標系。圖中共有12 個參數，參數s代表平移位移、θ 代表轉角位移，下標第一項表示位移方向，第二項表示節點編號。對單元中的自由度按照圖4 進行壓縮，得到單元端部位移和節點力如下式：

圖3 單元節點位移和節點力分量Fig.3 Element node displacement and node force components

圖4 壓縮后的單元節點位移和節點力分量Fig.4 Element node displacement and node force components after compression

基于壓縮后的位移參數分別使用線性插值和埃爾米特插值方法構建單元形函數，使用節點位移表示沿錐形單元軸向任意一點的軸向、側向以及扭轉位移如下[28]：

2.3 總勢能

高性能樁單元法通過在勢能表達式中引入土體應變能實現樁土作用的內置計算，其總勢能表示如下：

式中：Π 為總勢能；UE為單元應變能；US為土體應變能；W為外力作功。

單元應變能可采用式(8)計算：

式中：E為樁體的彈性模量；G為樁體的剪切模量；A(x)為截面面積；Iy和Iz分別為關于y軸和z軸的截面慣性矩；J為極慣性矩。

土體應變能為單元不同深度處的樁側彈簧儲存勢能之和，在數值計算中為了得到土體應變能的顯示表達式，采用如下的高斯-勒朗德積分表達式計算：

式中：n為高斯積分點的個數；TFL和TFS分別為樁側水平抗力和豎向摩阻力的圓錐效應系數[29-30]；Wm為第m個高斯點的積分權重；ρm和xm分別為第m個高斯積分點的水平和豎向位移；κρ,m和κx,m分別為第m個高斯積分點的p-y曲線和t-z曲線在相應位移處的切線斜率。

外荷載作功可按式(10)計算：

式中：Fi和si分別為節點兩端12 個自由度上每一個自由度可能的外荷載和對應的位移。

最小勢能原理表明，結構單元內的全部勢能之和總傾向于達到最小，因此可以通過總勢能微商為0 的條件建立錐形單元的平衡方程：

2.4 切線剛度

通過對總勢能表達式進行二階變分可以得到有限元分析所需的切線剛度矩陣。高性能樁單元法中的切線剛度矩陣由線性剛度矩陣、幾何剛度矩陣和土體剛度矩陣組成，它們的矩陣形式和組合關系如下：

式中：kL為線性剛度矩陣；kG為幾何剛度矩陣；kS為土體剛度矩陣；E為樁身彈性模量；L為單元長度；αi、βi為錐形剛度因子；為幾何剛度因子，計算式見LIU 等[31]；n為高斯點個數、i為高斯點序號、α 為高斯點權重系數，計算式見LI 等[25]。

整合前述三個矩陣便能夠得到總體切線剛度矩陣如下：

2.5 求解方式

高性能樁單元法采用牛頓-拉夫遜法增量迭代算法求解單元剛度方程。并且對于樁受水平荷載時可能產生的大撓度，使用更新拉格朗日描述[32]進行剛度增量逼近，根據最后已知的單元狀態建立平衡條件，總變形是通過一系列旋轉變換而不是通過求和過程得到的，因此可以考慮任意大的構件變形。

綜上，高性能樁單元法能夠考慮樁土作用、變截面、大變形等因素，可以較好地滿足風機支撐結構的一體化計算需求，因此本文將在驗證上述算法正確性的基礎之上，對包括樁基和塔筒在內的海上風機單樁支撐結構開展計算分析。

3 算法驗證

3.1 單層土體有效性驗證

本小節選取均質粉砂層作為樁基埋置土層，大直徑鋼管樁作為樁基形式進行高性能樁單元的單層土體有效性驗證。土層和樁基的參數詳見表1 和表2。

表1 單層土體參數[23]Table 1 Single layer soil parameters[23]

表2 樁基參數Table 2 Pile parameters

樁土相互作用方面，本文使用基于靜力觸探方法(CPT)的p-y曲線以考慮海上風電單樁的大直徑效應，其中砂土p-y曲線選擇SURYASENTANA和LEHANE[26]基于有限元分析提出的剛性短樁p-y曲線；粘土則選擇GUO 等[27]基于實驗數據提出的粘土p-y曲線。p-y曲線分別如式(16)和式(17)所示：

式中：pcu為CPT 得到的土層錐尖阻力； σvD'為有效上覆土壓力。

1)0≤z≤3時:

2)z>3D時:

式中：p為單位長度土抗力；pcu為CPT 得到的土層錐尖阻力；y為樁體水平位移；pu為極限水平土抗力；D為樁基外徑；z為樁深。

根據上述公式計算得到土層不同深度砂土以及粘土的p-y曲線如圖5 所示。

綜上，對于臨床疑診為TIO，99mTc-HYNIC-TOC SPECT顯像結果陰性的患者，68Ga-DOTA-TATE PET/CT可作為補充檢查有效地檢出TIO致病腫瘤。TIO致病腫瘤68Ga-DOTA-TATE PET/CT圖像上均可見生長抑素受體高度表達。另外，同機CT具有一定的特點，骨組織腫瘤局部可見溶骨性骨質破壞或局灶性骨質密度增高，四肢長骨病變常呈偏心性生長，部分累及骨皮質；軟組織腫瘤多呈密度均勻等或低密度結節，具有以上影像學特點傾向于TIO致病腫瘤的診斷。

分別使用高性能樁單元和實體單元進行驗證計算。其中將高性能樁單元數目劃分為8、10、20、40、80 以研究單元數目對計算結果準確性的影響。同時在有限元分析軟件中建立基于實體單元的離散彈簧模型，模型共劃分3200 個單元，沿著軸向布置80 層彈簧以保證計算結果的準確性。基于上述模型對高性能樁單元的高效性和準確性進行驗證。

圖6 顯示了在50 000 kN 水平荷載，底部固支，不同單元數目下高性能樁單元的水平位移計算結果。可以看出隨著單元數目的增加，采用不同單元數目計算得到的水平位移曲線幾乎重合，這說明單層土體中，高性能樁單元計算結果幾乎不受單元數目影響，對于單元數目的要求是非常低的，具有高效性。而圖7 顯示了高性能樁單元數目為8 時與實體單元計算結果的對比曲線，兩者樁頂水平位移計算結果一致，這驗證了高性能樁單元在單層土體中計算的準確性。

圖6 不同單元數目下高性能樁單元水平位移曲線Fig.6 Horizontal displacement curve of efficient element under different number of elements

3.2 多層土體有效性驗證

分層土體中的p-y曲線相對于單層土體具有更高的非線性。在離散彈簧單元法中常通過增加彈簧和單元數量來考慮土體高度非線性的影響，從而增加了計算成本。而高性能樁單元法將土體抗力集成在單元內部，單元劃分不受土彈簧布置密度的影響，減少了單元劃分量，因此能夠提高非線性問題的計算效率。

研究選取淤泥層-黏土層-粉砂層作為樁基埋置土層，假設土層的每一層未均質，土層參數參見表3，樁基參數見表2。樁土相互作用基于2.1 節中的p-y曲線得到。

表3 多層土體參數[33]Table 3 Multi layer soil parameters[33]

高性能樁單元數目分別劃分為10、20、40、80；同時將實體單元分別置于10 層彈簧和80 層彈簧的條件下進行計算，以對比土壤高度非線性下不同單元計算結果受單元數目影響的差異。

水平荷載下的計算結果對比如圖8 所示，可以得知在高度非線性土體下，高性能樁單元計算結果仍幾乎不受單元數目影響，且與80 層彈簧下的實體單元計算結果符合良好。反觀實體單元模型，當彈簧布置層數為10 時，水平位移計算結果顯著偏大，需要劃分更多節點才能獲得滿足精度的計算結果，這從側面驗證了高性能樁單元于復雜土體中計算的優越性。

圖8 不同單元數目下兩類單元水平位移-荷載對比曲線Fig.8 Lateral displacement-load curve of two kinds of element under different element number

3.3 海上風機支撐結構算例驗證

本小節使用高性能樁單元對簡化的風機支撐結構進行計算。首先將風機支撐結構簡化為錐形圓筒和定截面圓筒的組合，分別用于模擬塔筒和單樁，塔筒和單樁均采用高性能單元進行模擬，因此不同組件間采用單元節點進行連接。 模型的幾何參數如表4 和圖9 所示，基礎埋置土層參數參照某風電場地質數據[23]選取，如表5 所示。采用高性能樁單元進行計算，將結構劃分為7 個單元，其中塔筒段使用5 個單元，基礎段使用2 個單元。同時設置基于離散彈簧法的實體單元模型進行對比驗證，實體單元模型如圖10 所示，包含61 161 個單元。

表4 風機支撐結構參數Table 4 Parameters of wind turbine support structure

表5 風機基礎埋置土層參數[23]Table 5 Parameters of buried soil layer for wind turbine foundation[23]

圖9 風機支撐結構參數Fig.9 Parameters of wind turbine support structure

圖10 實體單元有限元模型示意圖Fig.10 Schematic diagram of solid element finite element model

在海洋環境下，風機支撐結構所受的荷載通常包括葉輪氣動荷載、塔筒風荷載以及浪流荷載等，其中葉輪氣動力一般起控制作用[34]，本節基于常規的風機荷載值進行工況設計。葉輪氣動力通常包括水平作用力T與彎矩M，本研究基于莫繼華[35]使用葉素-動量理論計算得到的某近海2 MW 葉輪氣動力序列(如表6 所示)確定了荷載取值的大致范圍。設計了作用于塔筒頂端的500 kN水平荷載工況和500 kN 水平荷載、10 000 kN·m 彎矩組合工況。對模型底部采用兩種不同的約束狀態，分別為樁底固支和樁底無約束，以分別模擬端承樁與摩擦樁，共計4 種工況。可得四種工況下的樁身變位曲線如圖11(a)～圖11(d)所示。從圖中可知，對于風機單樁支撐結構計算，高性能樁單元的位移計算結果與實體單元趨近于一致，具有很高的精確性和適用性。并且對比兩種模型所使用的單元量，使用高性能樁單元對計算資源的消耗量遠小于實體單元，這充分體現出了高性能樁單元的高效性。

表6 基于葉素-動量理論得到的2 MW 葉輪氣動力值[35]Table 6 Aerodynamic value of 2 MW impeller based on blade element momentum theory[35]

圖11 4 種工況下樁體水平變形曲線Fig.11 Horizontal displacement curve of pile under four working conditions

4 大直徑單樁水平承載性能參數分析

圖12 不同樁長、樁徑樁水平位移曲線Fig.12 Horizontal displacement curves of piles with different pile lengths and diameters

根據《港口工程樁基規范》(JTS 167-4-2012)[36]，樁基的變形模式可分為：彈性長樁、中長樁和剛性樁，可以根據樁長(L)與樁基相對剛度特征值(T)的大小關系來判斷樁基的變形模式，相對剛度特征值通常采用式(18)計算，評判標準如表7 所示：

表7 樁基變形模式評判標準Table 7 Evaluation criteria for pile foundation deformation mode

式中：T為相對剛度特征值；Ep為樁材料的彈性模量；Ip為樁截面的慣性矩；m為樁側地基土抗力系數隨深度增長的系數；b0為樁的換算寬度。

為研究大直徑樁變形模式與相對剛度特征值的關系，研究將計算結果根據變形模式進行分類，并計算了各情況下相對剛度特征值與樁長的比值，計算結果如表8 所示，表8 中以深淺兩種顏色標記了不同情況下樁的變形模式，兩種變形模式的交界位置對應的L/T值可以作為判斷兩種變形模式的L/T臨界值。

表8 不同參數組合下樁基L/T 值Table 8 L/T value of pile under different parameter combinations

從表8 中可以得到大直徑樁的剛性樁-中長樁臨界L/T值為0.75～1.10，遠遠小于規范中規定的2.50，與前人研究中得到的臨界L/T值0.8～1.4[37]相較具有良好的一致性。

以樁頂水平位移為20 mm 作為樁基達到正常使用極限狀態的標志，計算不同樁長和樁徑條件下的樁基水平承載力如圖13 所示。從圖中可以看到，在給定樁長的條件下，水平承載力隨著樁徑的增大而近似呈線性增大，這說明直徑對樁基水平承載力的影響十分顯著，且不只體現在對剛度的增強效應上；在給定樁徑的條件下，水平承載力隨樁長的變化則呈現顯著的非線性。對于小直徑樁如3 m 樁，樁長增大時水平承載力增長趨勢趨于平緩，樁長達到25 m 后，水平承載力不再受樁長的影響，這說明對于大直徑樁，存在水平承載有效樁長問題，在水平荷載下，深部土體所起的水平抗力作用會隨樁長增加逐漸減小并趨近于零。本節中，求得3 m、4 m、5 m 樁基直徑對應的有效樁長分別約為25 m、35 m、45 m。

圖13 不同樁長-樁徑組合樁基水平承載力云圖Fig.13 Horizontal bearing capacity of pile with different lengths and diameters

5 結論與建議

本文提出了一種改進的歐拉-伯努利梁單元，以模擬變截面樁在外荷載下的內力與位移響應。本文通過對比傳統離散單元法與改進單元法在不同土層中的計算結果和計算資源消耗，驗證了改進單元在高度非線性土體中計算的高效性、準確性。在完成驗證的基礎上，將改進單元應用于海上風電單樁支撐結構的水平承載力研究，得到以下結論：

(1) 在進行海上風機單樁支撐結構計算中，高性能單元法相較于傳統梁單元方法具有更高的效率、相當的精確度以及更高的一體化建模程度。計算設置中，只需設置單元本身的參數，而不必對模型進行特殊處理，例如施加節點彈簧、變截面段等效等。基于合適的p-y曲線，高性能樁單元法能進行快速、準確的風機支撐結構一體化計算。

(2) 當前工程界常使用樁長-相對剛度特征值比值判斷樁體在水平荷載作用下的變形特性，其中劃分剛性短樁-中長樁的規范推薦臨界值為2.5，但該值對于大直徑樁適用性較差。本文和前人研究中得到該值的范圍分別為0.75～1.1 和0.8～1.4，二者具有較好的一致性。

(3) 大直徑單樁在水平荷載作用下，隨著樁長的增加，其水平承載力增加速率會逐漸減小直至趨近于0，存在有效樁長。實際工程中應當考慮上述效應，確定水平承載有效樁長，以幫助進行合理的單樁尺寸設計。