海上風機大直徑單樁動態(tài)打樁過程數值仿真分析*

蔡舒鵬，張永康※，金 曄，薛 馳，張 笛，霍小劍，林 峰

（1.廣東工業(yè)大學機電工程學院，廣州 510006；2.中鐵建港航局集團有限公司，廣東珠海 519075；3.江蘇中天科技股份有限公司，江蘇南通 226000；4.武漢理工大學，武漢 430070；5.武漢船用機械有限公司，武漢 430080）

0 引言

隨著人類對能源資源需求的不斷增長以及陸地資源的逐步緊缺，人類活動不斷向海洋延伸，世界各國紛紛把海洋作為獲取能源資源和發(fā)展經濟的重要方向，通過海洋油氣開發(fā)、海上風電開發(fā)、港口碼頭建設、跨海大橋建設等發(fā)展向海經濟。我國海岸線長1.8萬km，風能資源豐富，據估計近?？砷_發(fā)的風能約7.5億kW，是陸地的2.96倍。“十四五”規(guī)劃中，我國將海上風電作為解決能源危機、降低環(huán)境污染、實現(xiàn)“雙碳目標”的國家戰(zhàn)略。

近年來，隨著國家“雙碳目標”戰(zhàn)略的不斷推進，國內海上風電場的開發(fā)建設速度也明顯加快，海上風電場的選址逐漸向離岸更遠、水深更大的方向發(fā)展[1-5]。在海上風電場建設中，基礎施工是項目建設中最重要的一個關鍵環(huán)節(jié)，單樁基礎作為樁基形式的重要種類，在海上風機基礎中占據相當比例。海上打樁作業(yè)一般由大型起重船或平臺作為載體，利用其上的起重機通過液壓錘吊打方式進行海上風電高樁承臺樁和導管架樁等施工[6]。隨著海上風電場離岸距離的不斷提升，單樁基礎中樁的尺度參數也在向著超長超大和超重方向發(fā)展，關于超大直徑單樁貫入土體過程中海水的涌入及土的力學響應的研究還相對較少[7-10]，因此研究超大直徑單樁動態(tài)打樁過程中樁-土-水相互之間的耦合效應有非常重要的意義。

國外學者對樁-土之間的相互作用進行了廣泛深入的研究，宋玉普等[11]考慮樁-土結構的耦合作用對導管架平臺進行了優(yōu)化設計，張兆德等[12]比較了ALE（Arbitrary Lagrangian Eulerian）和CEL（Coupled Eulerian Lagrangian）方法在計算巖土貫入問題時的優(yōu)缺點，發(fā)現(xiàn)ALE算法相比CEL算法更能夠保證土應力的穩(wěn)定性，并能夠保證材料-邊界-網格的一致性運動，但在預測樁端極限承載力時結果偏大。王志強等[13]基于LSDYNA通過ALE算法對靜壓單樁沉樁過程進行了數值仿真，并討論了樁頭錐角、樁土摩擦因數和土壤分層對實心樁沉樁擠土效應的影響。王娜娜等[14]基于LS-DYNA通過ALE算法研究了小尺寸混凝土樁打樁動力下沉過程中土體的應力、變形和擠土效 應。Hamann等[15]基 于ABAQUS通 過CEL算 法研究 了土壤在部分排水的情況下的沉樁過程，并研究了滲透率對沉樁和周圍土壤壓力的影響。Wang等[16]研究了超大直徑鋼管樁在動態(tài)打樁的軸向力作用下面臨的屈曲變形風險，發(fā)現(xiàn)控制橢圓度在30 mm以下對于減小鋼管樁屈曲具有重要意義。

大尺寸單樁的動態(tài)打樁過程如圖1所示。

圖1 大尺寸單樁的動態(tài)打樁過程

本文旨在通過多物質ALE流固耦合方法對海上風機大直徑單樁動態(tài)打樁過程進行數值仿真分析。以江蘇啟東H3海上風力發(fā)電場建設過程中的沉樁為工程背景，使用商用有限元軟件LS-DYNA進行建模，對動態(tài)打樁過程中土體的應力、應變、體積分數以及沉樁阻力的變化規(guī)律進行仿真分析。為海上風機大直徑單樁的動態(tài)打樁過程提供數值仿真建模參考和實用的分析方法。

1 多物質ALE流固耦合算法的原理

1976年，美國Lawrence Livermore國家實驗室J O Hallquist博士主持開發(fā)完成DYNA程序系列，最初是作為軍工上武器設計的分析工具。后來其版本和功能經過不斷地更新和完善，使得DYNA程序的應用范圍逐漸從軍事推廣到民用、商業(yè)等更多領域。在1988年，Hallquist創(chuàng)建了LSTC公司，將DYNA程序引入商業(yè)化發(fā)展軌道，并將之更名為LS-DYNA；1997年LSTC公司并與ANSYS公司合作，實現(xiàn)了LS-DYNA與ANSYS前后處理的連接，大大加強了LS-DYNA的前后處理能力和通用性，同時新開發(fā)了后處理程序LS-POST，極大地促進了LS-DYNA的發(fā)展。

LS-DYNA是非線性顯示動力學的鼻祖和先驅，其具有Lagrange、Euler和ALE算法，Lagrange算法的單元網格附著在材料上，隨著材料的流動而產生單元網格的變形。但是在結構變形過于巨大時，有可能使有限元網格造成嚴重畸變，引起數值計算的困難，甚至程序終止運算。ALE算法和Euler算法可以克服單元嚴重畸變引起的數值計算困難，并實現(xiàn)流體-固體耦合的動態(tài)分析。3種算法的主要特點如下。

ALE算法先執(zhí)行一個或幾個Lagrange時步計算，此時單元網格隨材料流動而產生變形，然后執(zhí)行ALE時步計算：（1）保持變形后的物體邊界條件，對內部單元進行重分網格，網格的拓撲關系保持不變，稱為Smooth Step；（2）將變形網格中的單元變量（密度、能量、應力張量等）和節(jié)點速度矢量輸運到重分后的新網格中，稱為Advection Step。

Euler算法則是材料在一個固定的網格中流動，在LS-DYNA中只要將有關實體單元標志Euler算法，并選擇輸運（Advection）算法。LS-DYNA還可將Euler網格與全Lagrange有限元網格方便地耦合，以處理流體與結構在各種復雜載荷條件下的相互作用問題。

為了說明問題，以一個2D的長方形變形為例來說明[17]。對于同一個物理過程，可以用不同的方式來描述：A Lagrangian，B Eulerian，C ALE，分別以上面3種方式來分別描述該物體的變形，分析其差別，如圖2所示[17]。

圖2 Lagrangian、Eulerian和ALE算法描述同一個物體變形時的區(qū)別

經過一個dt的時間變化后，比較3種描述的構形變化。

（1）A Lagrangian：對于拉格朗日描述，空間網格的節(jié)點與假想的材料點是一致的，也就是說，網格變形，材料也跟著網格變形，如圖2中A所示，所以對于大變形情況，網格可能發(fā)生嚴重畸變。

（2）B Eulerian：對于歐拉描述，兩層網格重疊在一起，一層空間網格固定在空間中不動，另一層附著在材料上隨材料在固定的空間網格中流動，并通過下面兩步來實現(xiàn)：首先，材料網格以一個拉格朗日步變形（像所描述的那樣）；然后拉格朗日單元的狀態(tài)變量被映射或輸送回到固定的空間網格中去。這樣網格總是不動和不變形的，相當于材料在網格中流動，如圖2中B所示，從而可以處理流體流動等大變形問題。

（3）C ALE：對于ALE（Arbitrary Lagrangian-Eulerian，任意拉格朗日-歐拉）描述，與歐拉描述一樣，有兩層網格重疊在一起，但空間網格可以在空間任意運動，其余與歐拉描述一樣，有物質的輸送在兩層網格中發(fā)生。如圖2中C所示，該方法可以整個物體有空間上的大位移，并且本身有大變形的非線性問題，如鳥撞飛機等問題。

海上風機基礎施工時的動態(tài)打樁過程會使樁周圍的土體產生大變形，若采用傳統(tǒng)的Lagrange算法，就會導致在模擬計算過程網格嚴重畸變、求解時間步長不斷減小，結果精度降低，甚至計算不收斂。因此，通過上述分析，可以采用ALE任意拉格朗日歐拉法將土體看成可流動的流體，而在空間上具有大位移的樁體為固體，土體網格的運動和樁體的運動分別獨立描述，從而避免了土體網格產生嚴重的畸變。

2 有限元ALE流固耦合模型的建立

2.1 材料本構模型

2.1.1 土壤材料模型

LS-DYNA中提供的可以用來模擬土壤的材料模型有20多種，其中最常用的是5號材料模型，對應關鍵字*MAT_SOIL_AND_FOAM。該模型能通過合理的試驗來確定其相關的材料參數，還能夠可靠地預測土壤的彈性性能、體積壓縮性能及其特殊的屈服特性，而且所需參數較少，在很多領域得到了廣泛的應用，因此本文也采用該材料本構模型。

該材料模型遵循屈服面無應變強化的修正Mohr-Coulomb塑性模型，其非線性D-P屈服函數φ是靜水壓力的二次函數，表示為：

J2為應力偏張量不變量，可由下式計算得到：

而應力偏張量σ′ij與材料的靜水壓力有關，表示為：

其中靜水壓力σm可以通過3個方向的主應力計算得到：

具體的材料參數如表1和表2所示。

表1 土壤材料的壓力-體積應變對應曲線值

表2 SOIL_AND_FOAM材料模型的參數及取值

2.1.2 樁體的材料模型

本文研究的樁體可以設置為線彈性材料或彈塑性材料，但因為主要研究對象是土體，而樁體的彈性模量比土體大得多，為了簡化模型，減少計算時間，樁采用線彈性剛體模型，對應LS-DYNA中的20號材料模型，對應關鍵字*MAT_020_RIGID，詳細的材料參數如表3所示。這里需要注意的是，根據實際工況，目前現(xiàn)場最大的打樁錘質量可達200 t，如果采用實心樁體建模，要控制其密度參數，使其質量達到約200 t，避免產生過大的慣性效應。

表3 樁體的材料模型參數

2.1.3 空物質材料模型

在樁體打入土體的過程中，由于占據了一部分土體的體積，土體上部靠近樁的部分會有隆起或下沉，引起“土拱效應”，為了讓土體可以向上方流動，需要在土體上方預先定義一定厚度的空物質層，該層包含的物質材料可以為水，也可以為空氣，其強度和剛度都為0，以此來容納上部土體的變形，因為本文研究的是海上的打樁過程，所以空物質材料選取水及對應的狀態(tài)方程?？瘴镔|材料在LS-DYNA中使用關鍵字*MAT_NULL定義，密度為1 000 kg/m3，狀態(tài)方程使用關鍵字*EOS_LINEAR_POLYNOMIAL定義，具體參數如表4所示。

表4 空物質材料水的狀態(tài)方程參數

2.2 邊界條件及ALE流固耦合定義

本文建立的有限元模型如圖3所示，其中樁體為圓柱實體模型，單元類型為SOLID164，而土和空氣單元類型也均為SOLID164。樁體的半徑為3 m，樁長度為25 m；通常情況下研究打樁過程地基土體應該為半無限大體，但在有限元建模中顯然不可按照實際情況建立，相關研究成果表明動態(tài)打樁對土體的水平影響范圍一般在10倍樁徑以內，因此土體的水平邊界取30 m，豎向取樁長度的1.2倍，即30 m；同時在土體上方定義了厚度為5 m的空物質層，意味著水深為5 m，同時讓土體可以向上方流動。以樁底端為坐標原點，樁受打擊方向為Z軸負方向建立有限元模型，為減少單元和降低計算量，樁和土體均采用1/4模型，在對稱邊界上分別設置對稱邊界條件，即在X面上限制X向的平動、Y和Z向的轉動，在Y面上限制Y向的平動、X和Z向的轉動；在土體外側使用關鍵字*NON-REFLECTING_BOUNDARY定義非反射邊界條件，這樣由樁產生的應力波在到達土體邊界時不會被反射，以達到半無限大體的效果。

在ALE流固耦合仿真中還需要通過關鍵字*ALE_MULTI_MATERIAL_GROUP定義多物質材料組分，這樣在計算中界面重構時這兩種組分可以互相混合。在本分析中，可流動的組分為土壤和水。在ALE流固耦合仿真中沒有摩擦和接觸部分的定義，但與之對應的是關鍵字*CONSTRAINED_LAGRANGE_IN_SOLID，其中SLAVE ID為拉格朗日描述的實體，即為樁，而MASTER ID則為ALE描述的實體，即為土壤和水，通過該關鍵字的定義，可以建立流體和結構之間相互作用的耦合機制。

建立好的模型單元總數為55 500個，節(jié)點總數為61 422個，樁與土壤接觸部分采用細化網格，其余部分采用均勻劃分的網格，如圖3的縮略圖所示。其中綠色實體為樁體，藍色實體為空物質層（水層），紅色部分為土壤。

圖3 1/4樁-土有限元模型及樁頭處的網格劃分示意圖

3 大直徑單樁動態(tài)打樁過程數值仿真

大直徑單樁動態(tài)打樁過程采用LS-DYNA進行顯式動力學分析，其中土壤的自重采用動態(tài)松弛方法施加，對土壤施加恒定的Z向重力加速度-9.8 m/s2，對應關鍵字為*LOAD_BODY_Z，并使用關鍵字*CONTROL_DYNAMIC_RELAXATION使重力載荷以動態(tài)松弛的方法施加；然后使用關鍵字*DEFINE_CURVE在樁頂施加三角波壓力載荷，最大錘擊力為2 000 kN，因為本文主要研究動態(tài)打樁過程中的樁-土相互作用，所以這里采用了簡化模型，沒有考慮液壓打樁錘的沖擊載荷作用在樁體時的能量和應力波的傳遞，而是將沖擊載荷直接作用在樁體上部，樁頂壓力載荷-時間曲線如圖4所示。

圖4 樁頂壓力載荷-時間曲線

3.1 土體應力場分析

大直徑單樁動態(tài)打樁過程中樁貫入不同深度時土體的徑向應力場（x向）分布如圖5所示?？梢钥吹皆跇抖烁浇耐馏w有明顯的應力集中區(qū)，應力分量的等值線具有類似于“應力泡”的層狀結構，這表明樁端土體具有類似球形的擴張形式。樁尖下的徑向應力為壓應力，其值隨著貫入深度的增加而增大，“應力泡”的影響范圍大致為5R，其中R為基樁半徑。在樁尖下面6R時，徑向應力迅速從高度壓縮的狀態(tài)減小到0。

圖5 單樁動態(tài)打樁過程中樁貫入不同深度時土體的徑向應力場分布（單位：Pa）

大直徑單樁動態(tài)打樁過程中樁貫入不同深度時土體的豎向應力場（z向）分布如圖6所示?？梢钥吹綐都庀碌呢Q向應力也為壓應力，其“應力泡”接近于橢圓形，相較于徑向“應力泡”，其在徑向的影響范圍要小，大致為2.5R，但其豎直方向的影響范圍更深，可達8R。

圖6 單樁動態(tài)打樁過程中樁貫入不同深度時土體的豎向應力場分布（單位：Pa）

大直徑單樁動態(tài)打樁過程中樁貫入不同深度時土體的切向應力場（xz向）分布如圖7所示?？梢钥吹綐抖松舷碌膽Ψ栂喾矗露藶槔瓚?，上端為壓應力，整體呈“X”形狀，并相較于樁尖附近，影響范圍大致為8R。

圖7 單樁動態(tài)打樁過程中樁貫入不同深度時土體的切向應力場分布（單位：Pa）

大直徑單樁動態(tài)打樁過程中樁貫入不同深度時土體的等效應力場分布如圖8所示?？梢钥吹降刃υ跇都庀虏咳〉米畲笾?，其“應力泡”形狀接近于心形，水平影響范圍為4R，豎直影響范圍最大可達8R。

圖8 單樁動態(tài)打樁過程中樁貫入不同深度時土體的等效應力場分布（單位：Pa）

橫向比較不同貫入深度下的各應力等值線的數值變化，發(fā)現(xiàn)隨著貫入深度的增加，樁尖處的壓應力也在逐漸增大，但各應力集中區(qū)的影響范圍和等值線的輪廓形狀變化不大，說明動態(tài)打樁過程對土壤的影響范圍是有限的。

3.2 土體應變場分析

大直徑單樁動態(tài)打樁過程中樁貫入不同深度時土體的等效塑性應變場分布如圖9所示?？梢钥吹綐抖烁浇牡刃苄詰兇笥跇秱鹊牡刃苄詰?，等效塑性應變在樁尖下部取得最大值，隨著樁端的貫入，樁周圍的土體逐漸受到壓縮，在達到屈服應力后開始產生塑性應變并形成塑性變形區(qū)，形成了類似于“水滴狀”的塑性變形區(qū)，隨著貫入深度的增加，塑性變形區(qū)也逐漸增大，證明有更大范圍的土體被擾動破壞，并發(fā)生了不同程度的塑性變形。

圖9 單樁動態(tài)打樁過程中樁貫入不同深度時土體的等效塑性應變場分布

3.3 土體體積分數分析

樁土在ALE流固耦合分析時，歐拉網格節(jié)點和與材料節(jié)點是分離的，因此網格節(jié)點的位置變化并不能反映材料節(jié)點的位置變化，所以不能通過觀察網格變形來得到物質的變形情況，但可以通過查看材料在網格中的體積分數的等值線來查看多物質材料的分界面情況，通過LS-PrePost后處理軟件得到的大直徑單樁動態(tài)打樁過程中樁貫入不同深度時土體的體積分數分布如圖10所示。圖中紅色部分代表該區(qū)域土體的體積分數為1，而藍色部分代表該區(qū)域水的體積分數為1。隨著樁體的不斷下沉，樁周的土體被不斷擠開，并且樁周附近開始有水貫入，在水和土壤界面處土壤呈“開口狀”，在距離樁身一定水平距離之外的土壤有微微隆起，這是受樁傳遞的應力沖擊波所擠壓的部分土體產生了位移。

圖10 單樁動態(tài)打樁過程中樁貫入不同深度時土體的體積分數分布

3.4 打樁阻力分析

樁-土之間的作用力可以通過二進制文件dbfsi輸出，在k文件中通過設置關鍵字*DATABASE_FSI進行激活，圖11所示為單樁動態(tài)打樁過程中樁-土之間的阻力變化?？梢婋S著每次錘擊的進行，樁貫入土體的深度增加，所受到的阻力也不斷增加，但每次擊打過后，由于應力波的反射作用和土體受壓縮出現(xiàn)的塑性屈服，樁會出現(xiàn)一定的反向位移，即回彈，從而導致沉樁阻力減小。如在0～1 s內樁上端壓力從最大值逐漸減小為零，樁體受到不斷減小的動態(tài)壓力和不斷增大的阻力的共同作用，因此會先向下移動一段距離再向上回彈一段距離，直到第二個階段1～2 s樁上端的壓力又開始逐漸增大，使樁重新開始向下移動。

圖11 單樁動態(tài)打樁過程中樁-土之間的阻力變化

需要注意的是，本文暫未考慮土的孔隙水壓力對沉樁阻力造成的影響。在實際工程操作中，打樁過程中樁周土體會積累很大的超靜孔隙水壓力，使樁周一定范圍內土體發(fā)生水力劈裂現(xiàn)象，導致樁周土體排水固結及強度恢復速度加快，停錘時間越長，土體強度恢復程度越大，造成后繼打樁困難甚至拒錘現(xiàn)象的發(fā)生。

4 結束語

本文使用商用有限元軟件LS-DYNA進行建模，通過多物質ALE流固耦合方法對海上風機大直徑單樁動態(tài)打樁過程進行了數值仿真分析。研究了大直徑單樁動態(tài)打樁過程中樁-土-水相互之間的耦合效應，得到了動態(tài)打樁過程中土體的應力、應變、體積分數以及沉樁阻力的變化規(guī)律，其結論如下。

（1）在樁端附近的土體有明顯的應力集中區(qū)，應力分量的等值線具有類似于“應力泡”的層狀結構，這表明樁端土體具有類似球形的擴張形式。發(fā)現(xiàn)隨著貫入深度的增加，樁端處的壓應力也在逐漸增大，但各應力集中區(qū)的影響范圍和等值線的輪廓形狀變化不大，說明動態(tài)打樁過程對土壤的影響范圍是有限的。

（2）樁端附近的等效塑性應變大于樁側的等效塑性應變，等效塑性應變在樁端下部取得最大值，隨著樁端的貫入，樁周圍的土體逐漸受到壓縮，在達到屈服應力后開始產生塑性應變并形成塑性變形區(qū)，形成了類似于“水滴狀”的塑性變形區(qū)，隨著貫入深度的增加，塑性變形區(qū)也逐漸增大，證明有更大范圍的土體被擾動破壞，并發(fā)生了不同程度的塑性變形。

（3）隨著樁體的不斷下沉，樁周的土體被不斷擠開，并且樁周附近開始有水貫入，在水和土壤界面處土壤呈“開口狀”，在距離樁身一定水平距離之外的土壤有微微隆起，這是受樁傳遞的應力沖擊波所擠壓的部分土體產生了位移。

（4）隨著每次錘擊的進行，樁貫入土體的深度增加，所受到的阻力也不斷增加。但每次擊打過后，由于應力波的反射作用和土體受壓縮出現(xiàn)的塑性屈服，樁體會先向下移動一段距離再向上回彈一段距離，從而導致沉樁阻力先增大后減小，直到下一錘擊又開始。