下沉過程耦合歐拉-拉格朗日有限元法分析

霍知亮，郎瑞卿，閆澍旺

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

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

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

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

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

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

1.1 歐拉體積分?jǐn)?shù)

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

圖2 通過歐拉體積分?jǐn)?shù)確定材料的步驟

1.2 接觸屬性

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

FP=kP/dP

(1)

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

2 應(yīng)用示例

2.1 砂土中樁的靜壓下沉

示例1應(yīng)用CEL有限元法，對砂土地基中的靜壓樁進(jìn)行貫入下沉分析。Sheng等[10]應(yīng)用ALE方法并引入自動荷載施加方案和平滑接觸離散技術(shù)，模擬砂土中樁的連續(xù)貫入問題。其中地基土采用摩爾-庫侖準(zhǔn)則，相應(yīng)的物理力學(xué)力學(xué)指標(biāo)參數(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àn)橄吔缧?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有限元法進(jìn)行基礎(chǔ)貫入分析時，為了避免在最初貫入下沉?xí)r網(wǎng)格變形過大，常把樁端模擬成錐形(例如Sheng等[10]、Henke[11])，并且為了模擬連續(xù)貫入，常使用zipper-type技術(shù)，即在樁端利用極細(xì)的光滑接觸剛性管，從而達(dá)到樁在貫入時與土平滑的接觸，更利于計算的收斂性(例如Mabsout等[12])。但做出這樣的簡化，在模擬樁下沉?xí)r會對樁-土相互作用和沉樁機(jī)理分析產(chǎn)生影響。應(yīng)用CEL有限元法無需對樁端截面進(jìn)行修改，不需其他簡化即可模擬樁的連續(xù)貫入。

2.2 吸力錨自重下沉

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

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

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

圖5 吸力錨自重下沉?xí)r力和下沉深度曲線

2.3 海底管線自沉分析

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

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

(2)

(3)

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

圖6 管線下沉示意圖

(4)

(5)

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

(6)

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

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

對于此類情況，Chatterjee等[24]采用考慮了大變形的RITSS方法進(jìn)行模擬研究，此外，Tian等[25]采用二次開發(fā)程序，并結(jié)合大型商用軟件ABAQUS的“網(wǎng)格-網(wǎng)格的映射方法”，計算分析了海管自沉?xí)r土體節(jié)點(diǎn)的應(yīng)力和變形，并進(jìn)行網(wǎng)格重劃分，處理海管下沉過程中的大變形的問題。采用CEL計算方法，并與式(2)的計算結(jié)果進(jìn)行分析對比，結(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)物準(zhǔn)確的應(yīng)力應(yīng)變狀態(tài)。

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