霍知亮,郎瑞卿,閆澍旺
(1. 天津大學 建筑工程學院, 天津 300072; 2.天津市市政工程設計研究院,天津 300051; 3.天津市軟土特性與工程環境重點實驗室, 天津 300384)
在進行基礎或結構物下沉貫入的數值分析時,常會遇到接觸處理困難、網格大變形造成的扭曲畸變等問題,引起程序收斂困難、計算失真。應用何種數值分析方法,準確的模擬下沉過程結構物和土的相互作用,已成為此領域的一個重要課題[1-2]。通常采用的拉格朗日分析方法,在模擬結構物下沉貫入等大變形分析時,網格的扭曲常會導致計算的不收斂或者引起嚴重的誤差[3]。純粹的歐拉分析可能導致數值擴散或者邊界處的材料流動未知等問題,在此類分析中并不多見。ALE方法在一定程度上克服兩者的缺點。但進行下沉貫入分析時,常常為了計算的收斂性做出一些假設、以及不能很好的模擬土體“沖剪”破壞形態等問題,影響貫入物與地基土體的相互作用機理[1]。
耦合的歐拉-拉格朗日(Coupled Eulerian - Lagrangian,CEL)方法,采用了歐拉計算方法與拉格朗日計算方法各自的特點,通過固定的歐拉網格及可以在網格中運動的材料,可有效的處理網格畸變等大變形分析中存在的解法困難,在近年的分析研究中逐漸得到應用。在海洋平臺樁靴下沉分析[4-7]、樁基礎沉貫分析和沉船問題[8]及其他海洋基礎承載力[9]中取得了一些成果。本文通過砂土中樁的靜壓下沉、吸力錨的自重下沉和海底管線自沉分析三個應用實例,表明CEL方法在模擬下沉貫入問題時的適用性及其優勢。其分析方法和成果可為其它類型基礎結構的下沉分析提供有益參考。
耦合歐拉-拉格朗日(CEL)法采用了歐拉計算方法與拉格朗日計算方法各自的特點。通常將基礎或結構物模擬為拉格朗日網格,地基土體模擬為歐拉網格,通過固定的歐拉網格,材料可以在固定的網格中運動,從而可有效的處理網格畸變等大變形分析中存在的解法困難。并且,其相應的接觸算法可以較好的模擬結構物與土的相互作用,得到準確的分析結果。兩種有限元法中物體的變形原理如圖1所示。

圖1 不同有限元法中物體的變形
在CEL有限元方法中,通過各個單元體中的歐拉體積分數(Eulerian Volume Fraction,EVF)來確定歐拉體在網格中的運動軌跡。若所建立的參考體全部填充于所建立的歐拉單元,則對應的歐拉體積分數EVF為1.0;反之,若參考體完全未填充于歐拉單元,則對應的歐拉體積分數EVF為0.0。若所建立的參考體部分填充于所建立的歐拉單元,即EVF小于1.0時,此單元中未被填充的部位被無屬性的“空”占據。其基本原理如圖2所示。

圖2 通過歐拉體積分數確定材料的步驟
在CEL有限元分析中,拉格朗日網格中的物質材料與歐拉物質材料兩者的接觸界面依據廣義接觸算法來實現,其接觸算法為:
FP=kP/dP
(1)
式中:FP為接觸界面上的接觸力;dP為相應位置的罰位移;kP為剛度,數值的大小與接觸的材料特性相關。這種接觸算法可自動選定接觸屬性中的主面和從面,可較好的處理大變形問題中的非線性問題。
示例1應用CEL有限元法,對砂土地基中的靜壓樁進行貫入下沉分析。Sheng等[10]應用ALE方法并引入自動荷載施加方案和平滑接觸離散技術,模擬砂土中樁的連續貫入問題。其中地基土采用摩爾-庫侖準則,相應的物理力學力學指標參數:黏聚力c為1 kPa,砂土內摩擦角φ為30°,剪脹角ψ為20°。砂土重度γ為20 kN/m3。在數值分析中彈性模量E取10 MPa,泊松比υ取0.3。樁的彈性模量E=100 GPa,泊松比υ=0.3。數值計算區域為消除邊界效應對計算結果的影響,模型的長和高分別取為2.4 m和4.8 m,靜壓樁的直徑為0.8 m,樁長為3 m,根據Sheng等[10]計算模型,樁端角度為60°,摩擦系數μ取0.01,樁的貫入深度為2.5 m。
應用CEL有限元法建立三維有限元計算模型,采用上述基本計算條件。計算得出的樁貫入阻力和貫入位移的關系,經過歸一化處理,得到計算結果分別如圖3和圖4所示。從計算結果可以看出,其計算結果和Sheng等采用ALE有限元法的計算結果相近。

圖3 樁貫入時抗力和貫入深度曲線

圖4 樁貫入時土體的變形
通常,應用ALE有限元法進行基礎貫入分析時,為了避免在最初貫入下沉時網格變形過大,常把樁端模擬成錐形(例如Sheng等[10]、Henke[11]),并且為了模擬連續貫入,常使用zipper-type技術,即在樁端利用極細的光滑接觸剛性管,從而達到樁在貫入時與土平滑的接觸,更利于計算的收斂性(例如Mabsout等[12])。但做出這樣的簡化,在模擬樁下沉時會對樁-土相互作用和沉樁機理分析產生影響。應用CEL有限元法無需對樁端截面進行修改,不需其他簡化即可模擬樁的連續貫入。
示例2應用CEL有限元法,對Luke[13]和Vásquez等[14]吸力錨自重下沉模型試驗和有限元分析進行模擬,并與其結果進行比較。在Vásquez等[14]對吸力錨自重下沉有限元分析中,通過開發一種摩擦接觸算法,并且應用重劃分網格技術,模擬吸力錨自重下沉和負壓下沉過程,計算吸力錨下沉時孔隙水壓力變化以及下沉深度和土體抗力之間的關系。
在Luke[13]的模型試驗中,試驗槽高為1.829 m,直徑1.219 m。土體為高嶺土,其十字板抗剪強度su隨深度線性增加。吸力錨長L=0.914 m,直徑D=10.2 cm,長細比L/D=9,壁厚t=0.81 mm,裙壁與土的摩擦系數為0.16。試驗具體情況和詳細參數見Luke等[13,15]。應用CEL有限元法建立三維有限元計算模型,模擬吸力錨自重下沉過程,計算得出的力和下沉深度的關系曲線如圖5所示。CEL有限元法計算結果在吸力錨下沉深度為0.3 m以內時,與試驗值吻合較好,當下沉深度超過0.3 m時,CEL有限元法比試驗值略大,但計算結果與室內試驗得出的力和下沉深度關系曲線趨勢基本相同。采用CEL有限元法無需開發摩擦接觸算法和重劃分網格技術,其應用更為簡便。
需要指出,顯示動力計算方法與普通靜態應力/位移分析相比,由于采用中心差分方法進行求解,在計算過程中對時間積分,其計算結果可能會產生波動或離散性(如圖5所示)。對于此種情況,分析者可根據計算結果和分析需要,采用例如Butterworth或者其它濾波函數[16]對曲線進行濾波處理。

圖5 吸力錨自重下沉時力和下沉深度曲線
示例3對海底管線自沉進行分析。海底管線作為傳輸石油、天然氣的介質和通道,是海洋油氣開發的重要組成部分。在遠離海岸的深海區域,海底管線通常直接鋪設在海床上,無需特意溝埋等工程保護措施,通過自沉或配重自沉下沉至海床內一定深度。通常,單孔的海底管線在鋪設過程中,其下沉深度約為0.1倍~0.5倍的海管直徑[17]。海管的自沉深度是分析海底管線屈曲等穩定性的一個重要參數。
對于海管下沉的理論分析,國內外許多學者依據承載力基本原理,結合管道基底特點進行了修正,提出了軟土地基海底管線下沉時豎向抗力的理論解[18-22],其計算公式為:
(2)
(3)
式中:V為管線承受的豎向荷載,kN/m;su0為不同深度處地基土不排水抗剪強度,kPa;Nc為承載力系數;As為管線下沉時的截面積,m2;fb為浮力系數;D為管線的外徑,m。根據Merifield等[22]、Randolph等[23]的研究成果,因管線下沉時地基土體將發生隆起現象,建議浮力系數fb宜取1.5。管線下沉示意圖如圖6所示,其計算公式為:

圖6 管線下沉示意圖


(4)


(5)
Nc承載力系數與管線的直徑和下沉深度有關,根據現有的研究成果,其表達式為冪函數的形式[18,20,22],如下所示:
(6)
式中:w為管線的下沉深度;a、b為擬合系數,其數值的大小和接觸面的粗糙程度、kD/sum的大小有關。系數a和b的取值見Aubeny等[20]。
采用CEL分析方法對飽和軟黏土地基中海管下沉進行模擬分析。取管線直徑D=0.8 m,模型區域為10D×8D,考慮安裝沉放和地基土的不排水特性,土體采用Tresca屈服準則,土體重度γ′取6.5 kN/m3,彈性模量和不排水抗剪強度的比值為500,即E/su=500,泊松比為0.49。地基土考慮了強度與埋深的線性關系,具體如圖6所示,其中地表土體sum取2.3 kPa,線性增加比例k取3.6 kPa/m。海管和地基土之間的接觸分別模擬為光滑和粗糙情況。
對于此類情況,Chatterjee等[24]采用考慮了大變形的RITSS方法進行模擬研究,此外,Tian等[25]采用二次開發程序,并結合大型商用軟件ABAQUS的“網格-網格的映射方法”,計算分析了海管自沉時土體節點的應力和變形,并進行網格重劃分,處理海管下沉過程中的大變形的問題。采用CEL計算方法,并與式(2)的計算結果進行分析對比,結果如圖7所示。從計算結果可以看出,CEL計算結果與式(2)結果較為接近,并且Chatterjee等[24]和Tian等[25]的計算結果位于光滑界面和粗糙界面所得的模擬結果之間。通過分析可看出,CEL方法可有效模擬海底管線大變形問題,所得的計算結果對于此類問題有較好的適用性。

圖7 正常固結黏土的計算結果
(1) 通過歐拉-拉格朗日的接觸算法,可較好的處理大變形問題中的非線性問題,得到基礎或結構物準確的應力應變狀態。
(2) 應用歐拉-拉格朗日有限元法進行大變形分析時,無需對貫入物截面進行修改,不需其他簡化即可進行連續貫入分析。無需開發摩擦接觸算法和重劃分網格技術,其應用更為簡便。