考慮應(yīng)變軟化的T－bar承載力CEL有限元分析

張新全，于 龍，2

（1.大連理工大學(xué)水利工程學(xué)院，遼寧大連116024；

2.大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024）

考慮應(yīng)變軟化的T－bar承載力CEL有限元分析

張新全1，于 龍1，2

（1.大連理工大學(xué)水利工程學(xué)院，遼寧大連116024；

2.大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024）

通過歐拉-拉格朗日（Coupled Eulerian-Lagrangian，CEL）大變形有限元對(duì)不考慮應(yīng)變軟化的T形觸探儀（T－bar）承載力系數(shù)進(jìn)行計(jì)算，其計(jì)算結(jié)果與已有塑性理論解和RITSS（Remeshing and Interpolation with Small Strain，即在小變形有限元的基礎(chǔ)上通過網(wǎng)格重剖分和應(yīng)力插值技術(shù)實(shí)現(xiàn)大變形有限元分析）的計(jì)算結(jié)果均吻合較好，驗(yàn)證了CEL計(jì)算模型的可靠性。在ABAQUS有限元平臺(tái)的基礎(chǔ)上進(jìn)行二次開發(fā)，編寫應(yīng)變軟化子程序計(jì)算出考慮應(yīng)變軟化的T－bar承載力系數(shù)，其計(jì)算結(jié)果與RITSS吻合較好，驗(yàn)證了應(yīng)變軟化子程序的可靠性。對(duì)T－bar周圍土體的軟化程度和流動(dòng)機(jī)制分析發(fā)現(xiàn)，應(yīng)變軟化使土體抗剪強(qiáng)度減小的同時(shí)還改變了土體的流動(dòng)機(jī)制，這兩個(gè)因素的綜合作用導(dǎo)致T－bar的承載力系數(shù)減小，從而揭示了應(yīng)變軟化對(duì)T－bar承載力系數(shù)的影響機(jī)理。

耦合的歐拉-拉格朗日大變形有限元；T形觸探儀；ABAQUS；應(yīng)變軟化；流動(dòng)機(jī)制

靜力觸探在現(xiàn)場(chǎng)測(cè)試中明確海洋軟黏土地基強(qiáng)度特性方面獨(dú)具優(yōu)勢(shì)。經(jīng)過十幾年的發(fā)展，觸探儀已由傳統(tǒng)的錐形發(fā)展到觸探過程中使土體呈完全流動(dòng)機(jī)制的多種形式：T－bar、盤形和球形觸探儀等，如圖1所示。相對(duì)傳統(tǒng)錐形觸探儀，全流觸探儀存在以下優(yōu)勢(shì)：對(duì)黏土剛度敏感性低；能通過較少修正情況下準(zhǔn)確測(cè)出黏土的不排水抗剪強(qiáng)度；承載力系數(shù)具有穩(wěn)定的理論解和數(shù)值解。不考慮黏土應(yīng)變軟化的T－bar承載力系數(shù)NnosofteningT－bar見式（1）：

式中：qult為T－bar單位投影面積上受土的極限反力；su為黏土的不排水抗剪強(qiáng)度。最初是基于塑性理論解得到的［1］，該方法考慮圓柱探頭表面的粗糙程度算出的范圍為9.1（完全光滑）～11.9（完全粗糙）。一般取=10.5來分析T－bar測(cè)試結(jié)果［2］。

眾多學(xué)者對(duì)T－bar觸探試驗(yàn)的承載力系數(shù)進(jìn)行了研究。Randolph M F等［3］、Martin C M等［4］根據(jù)極限分析方法得到了T－bar承載力系數(shù)的塑性理論解。范慶來等［5］通過ABAQUS有限元分析了應(yīng)變軟化對(duì)T－bar承載力系數(shù)的影響，但由于T－bar位移行程不夠，未得到考慮應(yīng)變軟化的T－bar承載力系數(shù)。郭紹曾［6］用CEL方法模擬了T－bar在無軟化無摩擦無重度黏土中的貫入，得到的T－bar承載力系數(shù)比較離散。周琪等［7］采用FLAC3D研究了應(yīng)變軟化對(duì)水平埋置條形錨板承載力的影響。張輝等［8］采用PFC2D研究了砂土中埋設(shè)T－bar的豎向抗拔特性。Low W E等［9］分析3個(gè)陸上工程和7個(gè)海上工程勘察資料，發(fā)現(xiàn)液性指數(shù)和塑性指數(shù)對(duì)全流貫入儀貫入阻力有一定影響。

圖1 觸探儀

實(shí)際情況下T－bar貫入時(shí)，周圍土體要經(jīng)歷較大的變形從而發(fā)生應(yīng)變軟化，這對(duì)于高敏感度土體的T－bar承載力系數(shù)造成顯著的影響。若想利用T－bar觸探儀精準(zhǔn)地確定黏土的抗剪強(qiáng)度，必須對(duì)考慮應(yīng)變軟化的T－bar承載力系數(shù)進(jìn)行充分的研究。Zhou H等［10］在T－bar和球形觸探儀數(shù)值計(jì)算與理論分析的基礎(chǔ)上，提出了考慮應(yīng)變軟化效應(yīng)的土體抗剪強(qiáng)度計(jì)算模型，見式（2）：

式中：FV1代表應(yīng)變軟化效應(yīng)，其中各參數(shù)的物理意義為：su0是土體沒有擾動(dòng)的初始抗剪強(qiáng)度；δrem是完全擾動(dòng)后的土體抗剪強(qiáng)度與su0的比值（即為敏感度ST的倒數(shù)）；ξ是土體單元的累積絕對(duì)塑性剪應(yīng)變；ξ95是土體單元受到95%擾動(dòng)所需的ξ值，此值反應(yīng)土的塑性和軟化快慢，典型的取值為10～50［11］。Zhou H等［10］采用式（2）所示軟化模型利用大變形有限元方法（RITSS）對(duì)T－bar貫入過程進(jìn)行了數(shù)值分析，考慮了不同組合的軟化參數(shù)δrem和ξ95對(duì) T－bar承載力系數(shù)的影響，給出了與δrem和ξ95的關(guān)系見式（3）：

本文在ABAQUS程序CEL算法中通過用戶子程序?qū)崿F(xiàn)了公式（2）的應(yīng)變軟化模型，對(duì)T－bar貫入均質(zhì)黏土過程進(jìn)行了分析，細(xì)致討論了應(yīng)變軟化參數(shù)ξ95對(duì)土體軟化程度和流動(dòng)機(jī)制的影響，揭示了理論分析和大變形有限元擬合所得ξT－bar數(shù)值不同的原因和應(yīng)變軟化對(duì)T－bar承載力系數(shù)的影響機(jī)理。

1 CEL方法簡(jiǎn)介

在拉格朗日有限元中，連續(xù)體的變形被定義為物質(zhì)點(diǎn)坐標(biāo)與時(shí)間的函數(shù)，單元結(jié)點(diǎn)位置隨著物質(zhì)運(yùn)動(dòng)而發(fā)生變化，這種單元形式的優(yōu)點(diǎn)在于能夠準(zhǔn)確描述系統(tǒng)中不同介質(zhì)之間的接觸界面，但在大變形過程中將遇到網(wǎng)格嚴(yán)重扭曲而導(dǎo)致計(jì)算中斷的問題。而歐拉單元?jiǎng)t把連續(xù)體的運(yùn)動(dòng)表達(dá)為空間坐標(biāo)與時(shí)間的函數(shù)，在分析過程中，歐拉網(wǎng)格不發(fā)生變形，而物質(zhì)可以在歐拉網(wǎng)格內(nèi)運(yùn)動(dòng)，因此歐拉單元不會(huì)產(chǎn)生網(wǎng)格扭曲，可以有效解決大變形問題中網(wǎng)格畸形帶來的弊端。CEL有限元方法是在綜合拉格朗日單元和歐拉單元優(yōu)點(diǎn)的基礎(chǔ)上發(fā)展起來的。拉格朗日網(wǎng)格內(nèi)的物質(zhì)與歐拉網(wǎng)格內(nèi)物質(zhì)之間的界面通過廣義接觸算法來模擬，雖然不如動(dòng)態(tài)接觸方法嚴(yán)格，但是可以保證在困難情況下算法的收斂。CEL有限元采用基于中心差分準(zhǔn)則的顯式動(dòng)力學(xué)解法來求解非線性微分方程體系，該算法中穩(wěn)定時(shí)間增量步長(zhǎng)在計(jì)算過程中自動(dòng)獲取，因而其分析過程相對(duì)其他大變形有限元比較簡(jiǎn)單。

2 T－bar貫入CEL有限元模型

由于對(duì)稱性，黏土地基模型和T－bar模型只取一半，如圖2所示。T－bar和黏土地基模型均采用八節(jié)點(diǎn)方塊單元，黏土地基模型采用歐拉單元模擬，T－bar模型采用拉格朗日單元模擬，在土體表面線上方設(shè)置一層厚度為1D（D為T－bar直徑）的空歐拉單元以適應(yīng)T－bar貫入過程中土體表面的隆起。整個(gè)模型厚度方向只取一個(gè)網(wǎng)格，模擬平面應(yīng)變問題。地基模型為均質(zhì)黏土，采用摩爾庫倫屈服準(zhǔn)則，黏土的硬度su0/γD取0.5，其中γ為土體重度（硬度為控制T－bar臨界埋深率的指標(biāo)［13］，臨界埋深率即T－bar剛好達(dá)到全流破壞機(jī)制時(shí)對(duì)應(yīng)的埋深率，此時(shí)T－bar承載力達(dá)到極限值）。土體彈性模量和泊松比分別取為E=500su0和ν=0.495。T－bar被模擬為剛體，T－bar與土體之間為通用接觸，采用粘滯摩擦系數(shù)α模擬T－bar與土體之間的摩擦，即接觸面的切應(yīng)力τ=αsu0，在軟化計(jì)算中α取為ST的倒數(shù)，T－bar由黏土表面向下貫入4D。首先計(jì)算T－bar在無軟化黏土中的承載力系數(shù)，之后在ABAQUS平臺(tái)上進(jìn)行二次開發(fā)，通過FORTRAN語言將式（2）所示應(yīng)變軟化模型寫入用戶子程序，用應(yīng)變軟化效應(yīng)函數(shù)FV1控制土體的軟化程度，保持α=δrem=0.2不變，取不同的ξ95值計(jì)算考慮應(yīng)變軟化的T－bar承載力系數(shù)。

圖2 T－bar貫入CEL有限元計(jì)算模型

正常情況下T－bar的貫入過程是準(zhǔn)靜態(tài)的，在CEL中是用動(dòng)態(tài)顯示算法代替準(zhǔn)靜態(tài)加載，為了在節(jié)省計(jì)算時(shí)間的基礎(chǔ)上讓動(dòng)態(tài)顯示算法逼近靜態(tài)隱式進(jìn)行歐拉分析之前需要對(duì)拉格朗日體的加載速度和歐拉體的網(wǎng)格尺寸分別做速度收斂性分析和網(wǎng)格收斂性分析［14］。本文做了速度收斂性分析來研究拔出貫入速度v對(duì)計(jì)算結(jié)果的影響：先將最小網(wǎng)格尺寸設(shè)定為h=D/20，選取三個(gè)貫入速度為v=0.2 D/s、0.1D/s和v=0.05D/s進(jìn)行計(jì)算，得到T－bar的抗力位移曲線見圖3。計(jì)算結(jié)果表明，貫入速度越大，T－bar承載力系數(shù)稍微變大，v=0.1D/s和v= 0.05D/s對(duì)應(yīng)的荷載位移曲線相當(dāng)接近，本文接下來的計(jì)算都采取了v=0.1D/s的貫入速度，可以認(rèn)為是準(zhǔn)靜態(tài)加載。本文網(wǎng)格收斂性分析結(jié)果如圖4所示。將貫入速度設(shè)為v=0.1D/s，選取三個(gè)不同的最小網(wǎng)格尺寸h=D/10、D/20和D/40進(jìn)行計(jì)算。計(jì)算結(jié)果表明，h=D/20的網(wǎng)格尺寸足夠精細(xì)，接下來的計(jì)算最小網(wǎng)格尺寸都取h=D/20。

圖3 速度收斂性分析

圖4 網(wǎng)格收斂性分析

3 T－bar貫入抗力-位移計(jì)算結(jié)果

首先分析不考慮應(yīng)變軟化情況下T－bar貫入過程，其抗力-位移曲線如圖5所示與 RITSS方法［10］和塑性理論解［10］的計(jì)算結(jié)果比對(duì)見表1。α=0時(shí)，本文CEL得到的比下限塑性解小3.5%，比上限塑性解小4.0%。α=0.2時(shí)，本文CEL得到的比下限塑性解小3.0%，比上限塑性解小3.2%，比RITSS方法小3.5%。 α=1時(shí)，本文CEL得到的比上、下限塑性解小3.6%，比RITSS方法小4%。計(jì)算結(jié)果表明，不管是否考慮摩擦，本文CEL方法計(jì)算出的無軟化T－bar承載力系數(shù)與RITSS方法和上、下限塑性解吻合較好，說明了本文CEL模型的可靠性。

圖5 不考慮應(yīng)變軟化的T－bar抗力-位移曲線

表1 NnTo－sboafrtening不同方法的計(jì)算結(jié)果

其次，分析考慮應(yīng)變軟化情況下的T－bar貫入過程。在計(jì)算時(shí)，摩擦和應(yīng)變軟化參數(shù)的設(shè)定與Zhou H等［10］一致，即α=δrem=0.2不變，ξ95分別取10、15、25和50，共4組工況。本文CEL方法得到的考慮應(yīng)變軟化的T－bar抗力-位移曲線見圖6，得到的與Zhou H等［10］和公式（3）的結(jié)果（ξT－bar取3.7）對(duì)比見表2。ξ95=10時(shí)，本文CEL得到比RITSS方法計(jì)算結(jié)果大0.5%，比式（3）的結(jié)果大1.1%。ξ95=15時(shí)，本文CEL得到比RITSS方法計(jì)算結(jié)果小3.8%，比式（3）的結(jié)果小4.4%。ξ95=25時(shí)，本文CEL得到比RITSS方法計(jì)算結(jié)果小8.6%，比式（3）的結(jié)果小7.0%。 ξ95=50時(shí)，本文CEL得到比RITSS方小8.8%，比式（3）的結(jié)果小6.7%?？紤]應(yīng)變軟化的計(jì)算結(jié)果表明，本文CEL方法通過應(yīng)變軟化子程序計(jì)算出的與Zhou H等的RITSS方法計(jì)算結(jié)果吻合較好，充分驗(yàn)證了本文應(yīng)變軟化子程序的正確性，同樣與式（3）ξT－bar取3.7時(shí)的結(jié)果吻合較好，說明本文CEL計(jì)算結(jié)果擬合出的ξT－bar也為3.7。

4 土體的軟化程度和流動(dòng)機(jī)制分析

首先，分析ξ95對(duì)T－bar周圍土體軟化程度的影響。ξ95=10、ξ95=25和ξ95=50三種工況土體軟化程度見圖7，由圖7可知：圖7（a）對(duì)于ξ95=10、ξ95= 25和ξ95=50，T－bar上方均存在一定區(qū)域的土體已經(jīng)完全軟化（即su/su0=0.2=δrem所代表的區(qū)域），ξ95越小，完全軟化區(qū)越大，這是因?yàn)棣?5越小，土體在較小應(yīng)變情況下就發(fā)生了完全軟化；圖7（b）由T－bar的右側(cè)邊緣到土擾動(dòng)區(qū)邊緣，su/su0由0.2過渡到1，且ξ95越小，擾動(dòng)區(qū)的寬度越小，抗剪強(qiáng)度的分布越不均勻；圖7（c）對(duì)于不同ξ95取值情況下，T－bar正下方土的軟化程度差異不大，這是因?yàn)檫@部分土體累積剪應(yīng)變相對(duì)較小。

圖6 考慮應(yīng)變軟化的T－bar抗力-位移曲線（α=δrem=0.2）

表2 NsTo－fbte anring不同方法的計(jì)算結(jié)果

其次，分析ξ95對(duì)T－bar周圍土體流動(dòng)機(jī)制的影響。土體流動(dòng)機(jī)制矢量圖見圖8，由圖8可知：圖8（a）α=0.2無軟化對(duì)應(yīng)土體流動(dòng)機(jī)制為旋轉(zhuǎn)破壞模式，與Martin C M等［15］、Einav I等［12］在上限極限分析中所假設(shè)的位移破壞模式吻合。圖8（b）ξ95不同，對(duì)應(yīng)的流動(dòng)機(jī)制形狀和大小均發(fā)生了變化，ξ95越小，流動(dòng)機(jī)制范圍越小。這是因?yàn)棣?5越小，土體軟化越快而抗剪強(qiáng)度分布越不均勻，土體總是沿著最容易的方向流動(dòng)。Einav I等［12］中的理論分析是建立在相同流動(dòng)機(jī)制基礎(chǔ)上（假定考慮軟化流動(dòng)機(jī)制同不考慮軟化情況相同），僅變化剪切帶上的土體抗剪強(qiáng)度，得出ξT－bar=3.85；而本文和Zhou H等［10］的大變形計(jì)算結(jié)果都表明ξT－bar=3.7，理論分析和大變形有限元分析所得結(jié)果有所區(qū)別的原因就在于大變形有限元分析考慮了應(yīng)變軟化對(duì)流動(dòng)機(jī)制形狀和大小的實(shí)時(shí)影響。

圖7 土體軟化程度云圖

圖8 T－bar周圍土體流動(dòng)機(jī)制矢量圖

5 結(jié) 論

本文通過編寫用戶子程序在CEL算法中實(shí)現(xiàn)了考慮土體的應(yīng)變軟化特性，并對(duì)T－bar貫入過程進(jìn)行了模擬，揭示了應(yīng)變軟化對(duì)T－bar承載力系數(shù)的影響機(jī)理。主要結(jié)論為：

（1）通過對(duì)CEL中T－bar貫入模型的速度和網(wǎng)格分析，本文建議使用h=D/20的最小網(wǎng)格尺寸和v=0.1D/s貫入速度。

（2）無軟化計(jì)算結(jié)果與已有文獻(xiàn)塑性理論解和RITSS分析結(jié)果吻合較好?？紤]應(yīng)變軟化的計(jì)算結(jié)果與已有文獻(xiàn)RITSS方法計(jì)算結(jié)果吻合較好，驗(yàn)證了本文計(jì)算模型和應(yīng)變軟化子程序的正確性。

（3）考慮應(yīng)變軟化時(shí)，T－bar周圍土體抗剪強(qiáng)度不均勻，導(dǎo)致T－bar流動(dòng)機(jī)制形狀和大小發(fā)生改變。流動(dòng)機(jī)制的改變和土體強(qiáng)度的弱化綜合作用導(dǎo)致T－bar承載力系數(shù)減小。

［1］ Randolph M F，Martin C M，Hu Y.Limiting resistance of a spherical penetrometer in cohesive material［J］.Geotechnique，2000，50（5）：573－582.

［2］ Randolph M F，Hefer P A，Geise J M，et al.Improved Seabed Strenght Profiling Using T-Bar Penetrometer［C］// Offshore Site Investigation and Foundation Behaviour’New Frontiers：Proceedings of an International Conference.Society of Underwater Technology，1998.

［3］ Randolph M F，Houlsby G T.Limiting pressure on acircular pile loaded laterally in cohesive soil［J］.Géotechnique，1984，34（4）：613－623.

［4］ Martin C M，Randolph M F.Upper－bound analysis of lateral pile capacity in cohesive soi1［J］.Géotechnique，2006，56（2）：141－146.

［5］ 范慶來，欒茂田，劉占閣.軟土中T型觸探儀貫入阻力的數(shù)值模擬［J］.巖土力學(xué)，2009，30（9）：2850－2854.

［6］ 郭紹曾.靜力觸探測(cè)試技術(shù)在海洋工程中的應(yīng)用［J］.巖土工程學(xué)報(bào)，37（S1）：207－211.

［7］ 周 琪，于 龍.考慮黏土應(yīng)變軟化的錨板承載力數(shù)值分析［J］.水利與建筑工程學(xué)報(bào)，2014，12（4）：124－128.

［8］ 張 輝，于 龍，王 博，等.砂土中埋設(shè)管線豎向抗拔特性研究［J］.水利與建筑工程學(xué)報(bào)，2015，13（5）：156－160.

［9］ Low H E，RandolphM F，Lunne T，et al.Effect of soil characteristics on relative values of piezocone，T－bar and ball penetration resistances［J］.Géotechnique，2011，61（8）：651－664.

［10］ Zhou H，Randolph M F.Resistance of full－flow penetrometers in rate－dependentand strain－softeningclay［J］. Géotechnique，2009，59（2）：79－86.

［11］ Chung S F，Randolph M F.Penetration resistance in soft clay for different shaped penetrometers［C］//Proc.2nd Int. Conf.on Site Characterisation，Porto.2004，1：671－678.

［12］ Einav I，Randolph M F.Combining upper bound and strain path methods for evaluating penetration resistance［J］.International Journal for Numerical Methods in Engineering，2005，63（14）：1991－2016.

［13］ White D J，Gaudin C，Boylan N，et al.Interpretation of T－ bar penetrometer tests at shallow embedment and in very soft soils［J］.Canadian Geotechnical Journal，2010，47（2）：218－229.

［14］ Chen Z，Tho K K，Leung C F，et al.Influence of overburden pressure and soil rigidity on uplift behavior of square plate anchor in uniform clay［J］.Computers and Geotechnics，2013，52（7）：71－81.

［15］ Martin C M，RandolphM F.Upper bound analysis of lateral pile capacity in cohesive soil［J］.Geotechnique，2006，56（2）：141－145.

CEL Finite Element Analysis of T－bar Bearing Capacity with Strain Softening

ZHANG Xinquan1，YU long1，2
（1.School of Hydraulic Engineering，Dalian University of Technology，Dalian，Liaoning 116024，China；2.State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian，Liaoning 116024，China）

The T－bar bearing capacity coefficients without considering strain softening was calculated by CEL（Coupled Eulerian-Lagrangian），a large deformation finite element method.The results from CEL method are in good agreement with the results of existing plastic theoretical solutions and RITSS（remeshing and interpolation with small strain，in which the large deformation finite element analysis is realized by remeshing technique and stress interpolation technique on the basis of small deformation finite element），thus the reliability of the CEL calculation model is verified.Based on the ABAQUS finite element platform，we developed a new compile strain softening subroutine to calculate T－bar bearing capacity coefficient considering strain softening.The results by strain softening subroutine are in good agreement with the results of RITSS，so the reliability of the strain softening subroutine is verified.By analyzing the degree of softening and flow mechanism of soil around T－bar，it can be found that strain softening not only reduces the shear strength of soil，but also changes its flow mechanism.The combined effects of these two factors lead to the decrease of the bearing capacity coefficient of T－bar.So，the effect mechanism of strain softening on the bearing capacity coefficient of T－bar was revealed.

CEL large deformation finite element；T－bar；ABAQUS；strain softening；flow mechanism

TU447

A

1672—1144（2016）05—0050—05

10.3969/j.issn.1672-1144.2016.05.010

2016－05－16

2016－06－14

國(guó)家自然科學(xué)基金項(xiàng)目（51539008，51479027）；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金項(xiàng)目（DUT15LK36）

張新全（1991—），男，湖北荊州人，碩士研究生，研究方向?yàn)榇笞冃斡邢拊治觥－mail：zhangxinquan7919@163.com

于 龍（1979—），男，遼寧遼陽人，博士，副教授，主要從事海洋基礎(chǔ)承載力和大變形有限元分析方面的工作。E－mail：longyu@dlut.edu.cn