霍知亮，郎瑞卿，閆澍旺

(1. 天津大學 建筑工程學院， 天津 300072; 2.天津市市政工程設(shè)計研究院，天津 300051; 3.天津市軟土特性與工程環(huán)境重點實驗室， 天津 300384)

在進行基礎(chǔ)或結(jié)構(gòu)物下沉貫入的數(shù)值分析時，常會遇到接觸處理困難、網(wǎng)格大變形造成的扭曲畸變等問題，引起程序收斂困難、計算失真。應(yīng)用何種數(shù)值分析方法，準確的模擬下沉過程結(jié)構(gòu)物和土的相互作用，已成為此領(lǐng)域的一個重要課題[1-2]。通常采用的拉格朗日分析方法，在模擬結(jié)構(gòu)物下沉貫入等大變形分析時，網(wǎng)格的扭曲常會導(dǎo)致計算的不收斂或者引起嚴重的誤差[3]。純粹的歐拉分析可能導(dǎo)致數(shù)值擴散或者邊界處的材料流動未知等問題，在此類分析中并不多見。ALE方法在一定程度上克服兩者的缺點。但進行下沉貫入分析時，常常為了計算的收斂性做出一些假設(shè)、以及不能很好的模擬土體“沖剪”破壞形態(tài)等問題，影響貫入物與地基土體的相互作用機理[1]。

耦合的歐拉-拉格朗日(Coupled Eulerian - Lagrangian，CEL)方法，采用了歐拉計算方法與拉格朗日計算方法各自的特點，通過固定的歐拉網(wǎng)格及可以在網(wǎng)格中運動的材料，可有效的處理網(wǎng)格畸變等大變形分析中存在的解法困難，在近年的分析研究中逐漸得到應(yīng)用。在海洋平臺樁靴下沉分析[4-7]、樁基礎(chǔ)沉貫分析和沉船問題[8]及其他海洋基礎(chǔ)承載力[9]中取得了一些成果。本文通過砂土中樁的靜壓下沉、吸力錨的自重下沉和海底管線自沉分析三個應(yīng)用實例，表明CEL方法在模擬下沉貫入問題時的適用性及其優(yōu)勢。其分析方法和成果可為其它類型基礎(chǔ)結(jié)構(gòu)的下沉分析提供有益參考。

1 耦合歐拉-拉格朗日(CEL)法

耦合歐拉-拉格朗日(CEL)法采用了歐拉計算方法與拉格朗日計算方法各自的特點。通常將基礎(chǔ)或結(jié)構(gòu)物模擬為拉格朗日網(wǎng)格，地基土體模擬為歐拉網(wǎng)格，通過固定的歐拉網(wǎng)格，材料可以在固定的網(wǎng)格中運動，從而可有效的處理網(wǎng)格畸變等大變形分析中存在的解法困難。并且，其相應(yīng)的接觸算法可以較好的模擬結(jié)構(gòu)物與土的相互作用，得到準確的分析結(jié)果。兩種有限元法中物體的變形原理如圖1所示。

圖1 不同有限元法中物體的變形

1.1 歐拉體積分數(shù)

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

圖2 通過歐拉體積分數(shù)確定材料的步驟

1.2 接觸屬性

在CEL有限元分析中，拉格朗日網(wǎng)格中的物質(zhì)材料與歐拉物質(zhì)材料兩者的接觸界面依據(jù)廣義接觸算法來實現(xiàn)，其接觸算法為：

FP=kP/dP

(1)

式中：FP為接觸界面上的接觸力；dP為相應(yīng)位置的罰位移；kP為剛度，數(shù)值的大小與接觸的材料特性相關(guān)。這種接觸算法可自動選定接觸屬性中的主面和從面，可較好的處理大變形問題中的非線性問題。

2 應(yīng)用示例

2.1 砂土中樁的靜壓下沉

示例1應(yīng)用CEL有限元法，對砂土地基中的靜壓樁進行貫入下沉分析。Sheng等[10]應(yīng)用ALE方法并引入自動荷載施加方案和平滑接觸離散技術(shù)，模擬砂土中樁的連續(xù)貫入問題。其中地基土采用摩爾-庫侖準則，相應(yīng)的物理力學力學指標參數(shù)：黏聚力c為1 kPa，砂土內(nèi)摩擦角φ為30°，剪脹角ψ為20°。砂土重度γ為20 kN/m3。在數(shù)值分析中彈性模量E取10 MPa，泊松比υ取0.3。樁的彈性模量E=100 GPa，泊松比υ=0.3。數(shù)值計算區(qū)域為消除邊界效應(yīng)對計算結(jié)果的影響，模型的長和高分別取為2.4 m和4.8 m，靜壓樁的直徑為0.8 m，樁長為3 m，根據(jù)Sheng等[10]計算模型，樁端角度為60°，摩擦系數(shù)μ取0.01，樁的貫入深度為2.5 m。

應(yīng)用CEL有限元法建立三維有限元計算模型，采用上述基本計算條件。計算得出的樁貫入阻力和貫入位移的關(guān)系，經(jīng)過歸一化處理，得到計算結(jié)果分別如圖3和圖4所示。從計算結(jié)果可以看出，其計算結(jié)果和Sheng等采用ALE有限元法的計算結(jié)果相近。

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

圖4 樁貫入時土體的變形

通常，應(yīng)用ALE有限元法進行基礎(chǔ)貫入分析時，為了避免在最初貫入下沉時網(wǎng)格變形過大，常把樁端模擬成錐形(例如Sheng等[10]、Henke[11])，并且為了模擬連續(xù)貫入，常使用zipper-type技術(shù)，即在樁端利用極細的光滑接觸剛性管，從而達到樁在貫入時與土平滑的接觸，更利于計算的收斂性(例如Mabsout等[12])。但做出這樣的簡化，在模擬樁下沉時會對樁-土相互作用和沉樁機理分析產(chǎn)生影響。應(yīng)用CEL有限元法無需對樁端截面進行修改，不需其他簡化即可模擬樁的連續(xù)貫入。

2.2 吸力錨自重下沉

示例2應(yīng)用CEL有限元法，對Luke[13]和Vásquez等[14]吸力錨自重下沉模型試驗和有限元分析進行模擬，并與其結(jié)果進行比較。在Vásquez等[14]對吸力錨自重下沉有限元分析中，通過開發(fā)一種摩擦接觸算法，并且應(yīng)用重劃分網(wǎng)格技術(shù)，模擬吸力錨自重下沉和負壓下沉過程，計算吸力錨下沉時孔隙水壓力變化以及下沉深度和土體抗力之間的關(guān)系。

在Luke[13]的模型試驗中，試驗槽高為1.829 m，直徑1.219 m。土體為高嶺土，其十字板抗剪強度su隨深度線性增加。吸力錨長L=0.914 m，直徑D=10.2 cm，長細比L/D=9，壁厚t=0.81 mm，裙壁與土的摩擦系數(shù)為0.16。試驗具體情況和詳細參數(shù)見Luke等[13，15]。應(yīng)用CEL有限元法建立三維有限元計算模型，模擬吸力錨自重下沉過程，計算得出的力和下沉深度的關(guān)系曲線如圖5所示。CEL有限元法計算結(jié)果在吸力錨下沉深度為0.3 m以內(nèi)時，與試驗值吻合較好，當下沉深度超過0.3 m時，CEL有限元法比試驗值略大，但計算結(jié)果與室內(nèi)試驗得出的力和下沉深度關(guān)系曲線趨勢基本相同。采用CEL有限元法無需開發(fā)摩擦接觸算法和重劃分網(wǎng)格技術(shù)，其應(yīng)用更為簡便。

需要指出，顯示動力計算方法與普通靜態(tài)應(yīng)力/位移分析相比，由于采用中心差分方法進行求解，在計算過程中對時間積分，其計算結(jié)果可能會產(chǎn)生波動或離散性(如圖5所示)。對于此種情況，分析者可根據(jù)計算結(jié)果和分析需要，采用例如Butterworth或者其它濾波函數(shù)[16]對曲線進行濾波處理。

圖5 吸力錨自重下沉時力和下沉深度曲線

2.3 海底管線自沉分析

示例3對海底管線自沉進行分析。海底管線作為傳輸石油、天然氣的介質(zhì)和通道，是海洋油氣開發(fā)的重要組成部分。在遠離海岸的深海區(qū)域，海底管線通常直接鋪設(shè)在海床上，無需特意溝埋等工程保護措施，通過自沉或配重自沉下沉至海床內(nèi)一定深度。通常，單孔的海底管線在鋪設(shè)過程中，其下沉深度約為0.1倍～0.5倍的海管直徑[17]。海管的自沉深度是分析海底管線屈曲等穩(wěn)定性的一個重要參數(shù)。

對于海管下沉的理論分析，國內(nèi)外許多學者依據(jù)承載力基本原理，結(jié)合管道基底特點進行了修正，提出了軟土地基海底管線下沉時豎向抗力的理論解[18-22]，其計算公式為：

(2)

(3)

式中：V為管線承受的豎向荷載，kN/m；su0為不同深度處地基土不排水抗剪強度，kPa；Nc為承載力系數(shù)；As為管線下沉時的截面積，m2；fb為浮力系數(shù)；D為管線的外徑，m。根據(jù)Merifield等[22]、Randolph等[23]的研究成果，因管線下沉時地基土體將發(fā)生隆起現(xiàn)象，建議浮力系數(shù)fb宜取1.5。管線下沉示意圖如圖6所示，其計算公式為：

圖6 管線下沉示意圖

(4)

(5)

Nc承載力系數(shù)與管線的直徑和下沉深度有關(guān)，根據(jù)現(xiàn)有的研究成果，其表達式為冪函數(shù)的形式[18,20,22]，如下所示：

(6)

式中：w為管線的下沉深度；a、b為擬合系數(shù)，其數(shù)值的大小和接觸面的粗糙程度、kD/sum的大小有關(guān)。系數(shù)a和b的取值見Aubeny等[20]。

采用CEL分析方法對飽和軟黏土地基中海管下沉進行模擬分析。取管線直徑D=0.8 m，模型區(qū)域為10D×8D，考慮安裝沉放和地基土的不排水特性，土體采用Tresca屈服準則，土體重度γ′取6.5 kN/m3，彈性模量和不排水抗剪強度的比值為500，即E/su=500，泊松比為0.49。地基土考慮了強度與埋深的線性關(guān)系，具體如圖6所示，其中地表土體sum取2.3 kPa，線性增加比例k取3.6 kPa/m。海管和地基土之間的接觸分別模擬為光滑和粗糙情況。

對于此類情況，Chatterjee等[24]采用考慮了大變形的RITSS方法進行模擬研究，此外，Tian等[25]采用二次開發(fā)程序，并結(jié)合大型商用軟件ABAQUS的“網(wǎng)格-網(wǎng)格的映射方法”，計算分析了海管自沉時土體節(jié)點的應(yīng)力和變形，并進行網(wǎng)格重劃分，處理海管下沉過程中的大變形的問題。采用CEL計算方法，并與式(2)的計算結(jié)果進行分析對比，結(jié)果如圖7所示。從計算結(jié)果可以看出，CEL計算結(jié)果與式(2)結(jié)果較為接近，并且Chatterjee等[24]和Tian等[25]的計算結(jié)果位于光滑界面和粗糙界面所得的模擬結(jié)果之間。通過分析可看出，CEL方法可有效模擬海底管線大變形問題，所得的計算結(jié)果對于此類問題有較好的適用性。

圖7 正常固結(jié)黏土的計算結(jié)果

3 結(jié) 論

(1) 通過歐拉-拉格朗日的接觸算法，可較好的處理大變形問題中的非線性問題，得到基礎(chǔ)或結(jié)構(gòu)物準確的應(yīng)力應(yīng)變狀態(tài)。

(2) 應(yīng)用歐拉-拉格朗日有限元法進行大變形分析時，無需對貫入物截面進行修改，不需其他簡化即可進行連續(xù)貫入分析。無需開發(fā)摩擦接觸算法和重劃分網(wǎng)格技術(shù)，其應(yīng)用更為簡便。