肖國(guó)豐，馬思琪，劉 洋，張興瑞，范慶來(lái)

(魯東大學(xué) 海洋巖土工程研究所， 山東 煙臺(tái) 264025)

隨著科學(xué)技術(shù)的不斷進(jìn)步以及電子信息產(chǎn)業(yè)的發(fā)展，計(jì)算機(jī)技術(shù)也取得了前所未有的革新與突破，數(shù)值模擬也隨之得以快速發(fā)展、廣泛應(yīng)用。ABAQUS作為一套功能強(qiáng)大的工程模擬有限元軟件，被廣泛應(yīng)用于巖土工程的實(shí)踐與研究，為巖土工程領(lǐng)域線性以及非線性問(wèn)題的研究做出了突出貢獻(xiàn)。同時(shí)ABAQUS對(duì)巖土材料的模擬具有很強(qiáng)的適用性，經(jīng)過(guò)長(zhǎng)時(shí)間的版本升級(jí)其功能也愈加完善。在實(shí)際巖土工程中，土體經(jīng)過(guò)長(zhǎng)年累月的沉淀，已形成相對(duì)穩(wěn)定的地質(zhì)狀態(tài)。而使用有限元軟件對(duì)巖土工程問(wèn)題進(jìn)行求解的時(shí)候，為了使有限元軟件盡可能真實(shí)可靠的模擬實(shí)際工況，需在加載之前使其地應(yīng)力達(dá)到一種平衡狀態(tài)，即土體內(nèi)部應(yīng)力與施加的重力荷載相平衡。只有初始地應(yīng)力滿(mǎn)足平衡要求，后續(xù)工況模擬才能得到符合實(shí)際工況的解。

在巖土工程的數(shù)值模擬分析中，初始地應(yīng)力的平衡是必須予以重視的問(wèn)題。針對(duì)ABAQUS對(duì)初始地應(yīng)力平衡提出的多種方法，目前已有一些成果對(duì)其進(jìn)行分析探討。于新等[1]應(yīng)用自動(dòng)平衡法與ODB導(dǎo)入法分別進(jìn)行計(jì)算，認(rèn)為對(duì)簡(jiǎn)易模型來(lái)說(shuō)兩種方法差別不大；楊金尤等[2]則對(duì)帶帽剛性樁復(fù)合地基進(jìn)行了計(jì)算，其中定義關(guān)鍵字法受條件約束較大，而ODB導(dǎo)入法對(duì)復(fù)雜模型存在需反復(fù)導(dǎo)入的弊端；而對(duì)簡(jiǎn)易模型來(lái)說(shuō)，導(dǎo)入ODB方法同樣可以明顯改善地應(yīng)力平衡效果[3]；代汝林等[4]同時(shí)對(duì)不同復(fù)雜程度的模型分別采用多種方法進(jìn)行計(jì)算，對(duì)于地質(zhì)條件復(fù)雜的巖土體采用自動(dòng)平衡法和應(yīng)力提取法效果明顯優(yōu)于其他方法；郭亞然等[5]則提出一種粘彈性邊界模型來(lái)解決動(dòng)力響應(yīng)分析中地應(yīng)力平衡問(wèn)題。

以上研究多對(duì)簡(jiǎn)易模型進(jìn)行了分析探討，而對(duì)復(fù)雜尤其是存在復(fù)雜約束條件及相互作用的模型中，僅使用傳統(tǒng)的地應(yīng)力平衡方法很難滿(mǎn)足平衡要求。因此，本文提出了一種分步平衡法，將ABAQUS生死單元技術(shù)與自動(dòng)平衡法結(jié)合起來(lái)，來(lái)解決較復(fù)雜巖土性質(zhì)和接觸條件下的地應(yīng)力平衡問(wèn)題。

1 基于生死單元技術(shù)的分步平衡法的提出

進(jìn)行初始地應(yīng)力平衡的目的是忽略工程結(jié)構(gòu)物的施工貫入對(duì)地基土體帶來(lái)的影響，直接對(duì)已經(jīng)安裝就位的工程結(jié)構(gòu)物的受荷情況進(jìn)行分析，而此時(shí)就要求其工程結(jié)構(gòu)物在地基土中已經(jīng)處于一個(gè)只存在應(yīng)力而不存在位移的相對(duì)穩(wěn)定狀態(tài)[6]。在簡(jiǎn)單的巖土工程問(wèn)題中，其初始地應(yīng)力平衡有著多種方法，其中自動(dòng)地應(yīng)力平衡方法作為“一步平衡法”的代表，由于其操作簡(jiǎn)單得到了廣泛應(yīng)用。但如果實(shí)際模擬工況較為復(fù)雜，工程結(jié)構(gòu)物與土體接觸較多且存在用戶(hù)定義場(chǎng)時(shí)，僅使用自動(dòng)地應(yīng)力平衡方法有時(shí)很難滿(mǎn)足地應(yīng)力平衡要求，甚至?xí)龅绞諗坷щy的問(wèn)題[7-8]。

本文基于ABAQUS中的生死單元功能結(jié)合自動(dòng)平衡法對(duì)地基初始地應(yīng)力場(chǎng)進(jìn)行平衡。該方法主要分為兩步，因此本文稱(chēng)之為“分步平衡法”。在分別建立地基土體和結(jié)構(gòu)物部件后，第一步“殺死”結(jié)構(gòu)物單元，僅對(duì)地基土體進(jìn)行地應(yīng)力平衡，模擬工程結(jié)構(gòu)物未貫入狀態(tài)；第二步“殺死”結(jié)構(gòu)物占據(jù)區(qū)域內(nèi)的土體單元，同時(shí)激活裝配好的工程結(jié)構(gòu)物部件，并設(shè)置滿(mǎn)足實(shí)際工況的接觸條件，進(jìn)行結(jié)構(gòu)與地基相互作用后的地應(yīng)力平衡，模擬工程結(jié)構(gòu)物安裝就位后的穩(wěn)定狀態(tài)。分步平衡法相較于自動(dòng)平衡法來(lái)說(shuō)，能夠有效緩解收斂性困難，同時(shí)地應(yīng)力平衡與位移歸零精度更高。

2 分步平衡法的應(yīng)用

以吸力式桶型基礎(chǔ)為例來(lái)驗(yàn)證此方法的地應(yīng)力平衡效果?？紤]有限元模型的對(duì)稱(chēng)性及數(shù)值模擬的計(jì)算效率，取地基土體和桶體的一半進(jìn)行建模分析。吸力桶基礎(chǔ)外徑D0為3 m，壁厚ts為0.1 m，長(zhǎng)度L為6 m。為了避免邊界效應(yīng)對(duì)數(shù)值模擬結(jié)果產(chǎn)生影響，地基土半徑D取為7D0，深度H取3.33L。圖1為劃分的有限元網(wǎng)格。在地基底部邊界上，約束三個(gè)方向的自由度，在側(cè)面邊界上對(duì)水平方向兩個(gè)自由度進(jìn)行約束，而在對(duì)稱(chēng)面上則對(duì)法向自由度進(jìn)行約束。

圖1 有限元網(wǎng)格

土體采用符合Mohr-Coulomb屈服準(zhǔn)則的理想彈塑性本構(gòu)模型。目前基于ABAQUS有限元軟件采用Mohr-Coulomb本構(gòu)模型對(duì)砂土地基進(jìn)行數(shù)值模擬計(jì)算分析時(shí)，地基土的變形模量多采用恒定值[9-10]。而實(shí)際上，砂土的變形模量隨著圍壓水平而改變，并不是一個(gè)恒定值，這種壓硬性可以采用Janbu提出的冪函數(shù)公式來(lái)描述[11]，如公式(1)所示。

(1)

式中：Es是砂土地基的側(cè)限壓縮模量;大氣壓強(qiáng)σat取100 kN/m2;κ決定了基準(zhǔn)應(yīng)力狀態(tài)下的土體剛度;λ反映了土體剛度的應(yīng)力依賴(lài)性程度;σm為平均主應(yīng)力。

由于ABAQUS中的Mohr-Coulmb彈塑性模型默認(rèn)采用變形模量，所以根據(jù)以下土力學(xué)公式將側(cè)限壓縮模量Es轉(zhuǎn)換為變形模量E，如式(2)所示。

(2)

式中：E為變形模量；ν為土體泊松比。參考相關(guān)規(guī)范以及Achmus等的建議[11]，本次數(shù)值模擬采用κ=600，λ=0.55，具體模型參數(shù)見(jiàn)表1所示。通過(guò)Fortran語(yǔ)言編寫(xiě)子程序進(jìn)行二次開(kāi)發(fā)，實(shí)現(xiàn)土體模量的應(yīng)力相關(guān)性。

表1 材料參數(shù)

桶土界面接觸行為采用“摩擦接觸對(duì)”算法模擬，法向上采用硬接觸，切向上采用Coulomb定律描述其力學(xué)性質(zhì)，摩擦系數(shù)f=tan(0.6φ)，其中φ為砂土內(nèi)摩擦角。

對(duì)此，分別采用自動(dòng)平衡法和本文提出的分步平衡法來(lái)進(jìn)行地應(yīng)力平衡。使用地應(yīng)力自動(dòng)平衡方法，是在一個(gè)地應(yīng)力分析步完成地應(yīng)力平衡，吸力桶與地基土體直接在地應(yīng)力分析步中設(shè)置接觸條件及相互作用；而分步平衡法則需兩個(gè)地應(yīng)力分析步來(lái)分步完成初始地應(yīng)力平衡，吸力桶基礎(chǔ)與土體的接觸及相互作用只有在吸力桶結(jié)構(gòu)激活時(shí)才起作用，能夠更加真實(shí)地反映實(shí)際工況。

3 平衡效果對(duì)比

分別對(duì)砂土地基定模量和變模量?jī)煞N情況下的地應(yīng)力平衡效果進(jìn)行評(píng)價(jià)。評(píng)價(jià)時(shí)，主要通過(guò)豎向應(yīng)力大小和豎向位移歸零程度兩方面進(jìn)行衡量。

3.1 定模量模擬效果對(duì)比

對(duì)于定模量情況，砂土地基的變形模量E取值為60 MPa。分別采用自動(dòng)平衡法和本文提出的分步平衡法來(lái)進(jìn)行地應(yīng)力平衡，在分析步中施加體力加載。

圖2為豎向應(yīng)力S33分布的比較，圖中應(yīng)力單位為Pa。通過(guò)圖2可以看出，豎向應(yīng)力隨著土體深度的增加而增大。豎直方向的自重應(yīng)力S33要滿(mǎn)足S33=γ·H。在地基底面，豎向自重應(yīng)力的大小理論上約等于220 kN/m2，分步平衡法得到的結(jié)果比其約低18%，而自動(dòng)平衡法的誤差則達(dá)到了28%。從自重應(yīng)力分布來(lái)看，分步平衡法得到的結(jié)果更接近于理論值，隨后進(jìn)行的外力加載分析結(jié)果也會(huì)更接近于真實(shí)情況。

圖2 豎向自重應(yīng)力分布的對(duì)比

另外，在自動(dòng)平衡法中，為了能夠收斂，吸力桶的重度一般假設(shè)為土體的重度，這對(duì)于固定式海上風(fēng)機(jī)桶基工程的受力分析可能影響不大，但是對(duì)于桶形基礎(chǔ)用于承受拉拔荷載時(shí)會(huì)造成一定誤差。而對(duì)于分步平衡法，在第二步激活桶體單元時(shí)，可以施加真實(shí)的鋼材料的重度，從而更加真實(shí)地反映實(shí)際工況。

另外，從圖2還可以看出，分步平衡法所得到的豎向應(yīng)力分布更加均勻，而自動(dòng)平衡法得到的豎向應(yīng)力在吸力桶周?chē)a(chǎn)生了明顯的局部錯(cuò)位情況。

初始地應(yīng)力平衡效果的優(yōu)劣程度還取決于其豎向位移U3歸零是否能夠滿(mǎn)足要求。一般巖土工程領(lǐng)域?qū)Φ貞?yīng)力分析步的豎向位移U3歸零的要求是小于10-4m量級(jí)。圖3為兩種方法在定模量模擬情況下的豎向位移的對(duì)比，單位為m。

圖3 豎向位移的對(duì)比

從圖3中可以看出，無(wú)論是采用哪種方法對(duì)定模量進(jìn)行數(shù)值模擬，其豎向位移歸零均能滿(mǎn)足量級(jí)要求。但是就位移分布規(guī)律來(lái)看，自動(dòng)平衡法模擬結(jié)果中，桶體內(nèi)部土體產(chǎn)生了類(lèi)似條紋狀的彩帶，出現(xiàn)了不符合規(guī)律的位移變化，效果不是特別理想。而分步平衡法模擬出來(lái)的位移變化情況，桶基礎(chǔ)中土體因?yàn)橥暗呢炄攵l(fā)生沉降，桶蓋周?chē)鷷?huì)因?yàn)樨炄攵a(chǎn)生部分隆起，更加符合模型試驗(yàn)和原位測(cè)試觀察到的實(shí)際現(xiàn)象[12-13]。

3.2 變模量模擬效果對(duì)比

采用相同的模型尺寸及子程序?qū)煞N方法的地應(yīng)力平衡效果進(jìn)行對(duì)比分析。在ABAQUS應(yīng)用子程序時(shí)，土體變形模量是在給定的場(chǎng)變量FV1區(qū)間內(nèi)進(jìn)行線性插值賦予具體數(shù)值的。圖4為子程序在兩種方法中調(diào)用后的模擬情況，圖中模量單位為Pa。

圖4 土體模量分布的對(duì)比K0=1-sinφ=0.37

根據(jù)式(1)和式(2)可以估計(jì)地基土體底面的變形模量理論值。由于砂土地基一般滿(mǎn)足K0固結(jié)條件，其側(cè)壓力系數(shù)約等于：

(3)

因此地基底面的平均應(yīng)力可以計(jì)算為：

(4)

將式(4)代入式(1)可得地基底面的側(cè)限壓縮模量Es=68 607 kPa，進(jìn)而根據(jù)式(2)得變形模量E=50 965 kPa≈5.1×107Pa。兩種方法得到的地基底面變形模量值相差不大，比理論值約低6%。但對(duì)于地基頂面，此處的平均應(yīng)力σm=0，根據(jù)式(1)和式(2)其變形模量理論值應(yīng)該約為0。從圖4看出，分步平衡法準(zhǔn)確反映了地基淺層土體的低模量水平，而自動(dòng)平衡法得到結(jié)果約為3 961 kPa，錯(cuò)誤高估了地基淺層土體的變形模量，導(dǎo)致地基土體的變形模量分布規(guī)律與式(1)的冪函數(shù)關(guān)系相差較大，對(duì)后續(xù)桶體在拉拔荷載作用下的位移預(yù)測(cè)帶來(lái)很大影響。

圖5為變模量情況下豎向應(yīng)力S33分布的比較，圖中應(yīng)力單位為Pa。與定模量模擬情況類(lèi)似，從自重應(yīng)力分布來(lái)看，分步平衡法得到的結(jié)果更接近于理論值，隨后進(jìn)行的外力加載分析結(jié)果也會(huì)更接近于實(shí)際情況。

圖5 豎向應(yīng)力分布的對(duì)比

圖6為兩種方法在變模量模擬情況下的豎向位移的對(duì)比，單位為m。通過(guò)比較可以看到，采用自動(dòng)平衡法得到的豎向位移量級(jí)已經(jīng)達(dá)到10-4m，這對(duì)于固定式海上風(fēng)機(jī)桶基工程的位移分析影響不大，但對(duì)于桶形基礎(chǔ)承受拉拔荷載時(shí)[14]的位移會(huì)產(chǎn)生一定誤差。而且，自動(dòng)平衡法所得到的淺層土體豎向位移場(chǎng)產(chǎn)生了較大的波浪形趨勢(shì)，位移分布極不均勻，這與實(shí)際情況有較大的差異。

圖6 豎向位移的對(duì)比

采用分步平衡法的地應(yīng)力平衡效果明顯提高，其平衡后的豎向位移在10-5次方量級(jí)，比自動(dòng)平衡法結(jié)果低一個(gè)量級(jí)，因此更加滿(mǎn)足地應(yīng)力平衡的要求。與對(duì)應(yīng)的定模量情況模擬得到的豎向位移分布規(guī)律類(lèi)似，位移場(chǎng)僅在吸力桶基礎(chǔ)臨近區(qū)域受到影響，桶內(nèi)土體因貫入的影響產(chǎn)生沉降，桶蓋周?chē)馏w受到擠壓而發(fā)生少量隆起，與實(shí)際情況吻合度較高[12-13]。

通過(guò)上述分析可以看到，分析過(guò)程中存在結(jié)構(gòu)與土體相互作用或者巖土體材料性質(zhì)比較復(fù)雜時(shí)，僅采用自動(dòng)平衡法，其平衡效果往往不理想或不能滿(mǎn)足平衡要求，甚至有時(shí)會(huì)出現(xiàn)收斂困難而導(dǎo)致數(shù)值計(jì)算任務(wù)中斷。本文通過(guò)吸力桶工程實(shí)例的地應(yīng)力平衡效果的對(duì)比，可以看出分步平衡法可有效改善地應(yīng)力平衡效果，得到符合實(shí)際情況的數(shù)值模擬結(jié)果，對(duì)類(lèi)似數(shù)值模擬具有一定的參考價(jià)值。

4 結(jié) 論

在巖土工程數(shù)值計(jì)算中，對(duì)于復(fù)雜工況，要獲得一個(gè)只存在地應(yīng)力而不存在位移的平衡狀態(tài)，僅采用自動(dòng)平衡法等一步平衡法往往不理想，甚至有時(shí)會(huì)出現(xiàn)收斂困難。本文基于ABAQUS軟件中的生死單元功能，提出了一種分步平衡法，通過(guò)分步創(chuàng)建桶土之間的接觸以更好地模擬地基在桶土相互作用過(guò)程中的實(shí)際響應(yīng)。通過(guò)分析，得到如下結(jié)論：

(1) 對(duì)于定模量情況下的兩種方法對(duì)比，分步平衡法得到的豎向應(yīng)力場(chǎng)更接近于理論解，豎向位移量級(jí)盡管與自動(dòng)平衡法相差不大，但是得到的豎向位移分布規(guī)律更接近于已有的模型試驗(yàn)結(jié)果。

(2) 對(duì)于變模量情況下的兩種方法對(duì)比，分步平衡法不僅得到的豎向應(yīng)力場(chǎng)更接近于理論解，而且豎向位移量級(jí)也比自動(dòng)平衡法低1個(gè)量級(jí)。

(3) 對(duì)于壓硬性砂土地基中桶土相互作用問(wèn)題，分步平衡法能夠合理反映地基淺層土體的低模量特征，而自動(dòng)平衡法得到的結(jié)果則錯(cuò)誤高估了地基淺層土體的變形模量，導(dǎo)致地基土體的變形模量分布規(guī)律與Janbu公式相差較大。