林志遠(yuǎn)，熊傳祥，謝永寧，黃盛鋒

(1.福州大學(xué) 環(huán)境與資源學(xué)院， 福建 福州 350116；2.自然丘陵山地地質(zhì)災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室， 福建 福州 350116；3.福建省地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室， 福建 福州 350116)

廣義塑性模型由Pastor等[1]1990年系統(tǒng)提出，包含黏土和砂土兩部分，該模型因具有簡潔形式、模塊化特性以及易于數(shù)值實(shí)現(xiàn)的特點(diǎn)受到學(xué)者的青睞，但前人對(duì)廣義塑性模型的研究主要圍繞砂土模型展開[2-4]，對(duì)黏土部分的研究目前仍較為欠缺(包括理論和模型應(yīng)用)[5]。國內(nèi)關(guān)于黏土廣義塑性模型的研究仍存在一些亟待解決的問題，如：張衛(wèi)東等[6]在ADINA軟件中實(shí)現(xiàn)了二次開發(fā)，但該研究僅針對(duì)適用于正常固結(jié)黏土的理論部分開展數(shù)值試驗(yàn)驗(yàn)證，并未涉及在超固結(jié)狀態(tài)下的可行性論證，對(duì)理論模型與二次開發(fā)研究尚不完整；費(fèi)康等[7]在黏土廣義塑性模型框架中考慮溫度影響，使得改進(jìn)后的模型能夠反映黏土在溫度循環(huán)作用下的變形行為，但仍欠缺對(duì)改進(jìn)的模型在力學(xué)加載條件下的驗(yàn)證。故完善黏土廣義塑性模型在各種不同應(yīng)力歷史、不同應(yīng)力路徑中理論與應(yīng)用研究為本研究的目的之一。

目前，已有不少學(xué)者借助FLAC3D提供的二次開發(fā)平臺(tái)(UDM)開展二次開發(fā)研究，并獲得顯著的成果[8-15]。FLAC3D采用的顯式有限差分算法能準(zhǔn)確模擬材料的塑性流動(dòng)與破壞，在求解三維巖土問題、動(dòng)力問題以及高度非線性問題上有突出的優(yōu)勢，尤其適合作為黏土廣義塑性模型數(shù)值實(shí)現(xiàn)的載體，以應(yīng)用于實(shí)際工程分析。最重要的是，目前還未有基于FLAC3D軟件的黏土廣義塑性模型二次開發(fā)研究，其意味著對(duì)黏土廣義塑性模型的工程應(yīng)用價(jià)值的認(rèn)識(shí)尚存不足，對(duì)比砂土廣義塑性模型(常用于大壩變形分析[16-18])其應(yīng)用價(jià)值可能被低估，具有深入探究的意義。

鑒于黏土廣義塑性模型研究尚存在不足(包括可參考研究成果較少，缺乏超固結(jié)土部分理論與應(yīng)用的研究，工程應(yīng)用價(jià)值認(rèn)識(shí)不足)，以及黏土廣義塑性模型的FLAC3D二次開發(fā)研究還未被開展等問題，本文開展了關(guān)于黏土廣義塑性模型的FLAC3D二次開發(fā)研究。

1 黏土廣義塑性模型及其差分形式

1.1 黏土廣義塑性模型

黏土廣義塑性模型適用于求解各種不同應(yīng)力歷史(正常固結(jié)與超固結(jié))黏性土的力學(xué)特性，且得益于該模型無需明確給出屈服面和塑性勢能面函數(shù)形式的特點(diǎn)，極大程度降低了數(shù)值實(shí)現(xiàn)的難度。

黏土廣義塑性模型的增量形式為:

dσ=De∶dε-dλDengL/U

(1)

式中：dσ，dε分別為應(yīng)力、應(yīng)變?cè)隽肯蛄?；ngL/U為單位塑性流動(dòng)方向向量，“L/U”分別表示加載和卸載；De為彈性剛度矩陣。

De=K·(m?m)+2G·[I-1/3(m?m)]

(2)

向量m=(1,1,1,0,0,0)T，I為單位矩陣；?為克羅內(nèi)克積(Kronecker product)，如對(duì)于向量u和v，(u?v)ij=uivi。K、G分別為彈性體積和剪切模量，其與平均有效應(yīng)力p′的關(guān)系為:

(3)

dλ為塑性乘子(plastic multiplier)，為一標(biāo)量，其方程形式為:

(4)

式中:n為單位加-卸載方向向量；HL為塑性模量。假設(shè)采用相適應(yīng)的流動(dòng)準(zhǔn)則，則有n=ngL/U，在三維應(yīng)力空間中表達(dá)式為[19]:

(5)

式中：p′為平均有效應(yīng)力，q為剪應(yīng)力，θ為Lode角。其中p′,q與應(yīng)力張量之間的關(guān)系分別為:

(6)

(7)

假設(shè)應(yīng)力狀態(tài)與Lode角無關(guān)，將公式(7)代入式子(5)中，可得三維應(yīng)力空間中的塑性流動(dòng)方向和加-卸載方向向量的表達(dá)式，如下式所示:

(8)

上式中，d為剪脹性關(guān)系:

d=(1+α)(Mc-η)

(9)

式中：α為材料參數(shù)；Mc為臨界狀態(tài)應(yīng)力比；η為應(yīng)力比(=q/p′ )。

式(4)中塑性模量HL是一個(gè)標(biāo)量，其表達(dá)式為:

(10)

式中：γ為模型參數(shù)；ζ、ζmax分別為記錄應(yīng)力歷史的函數(shù)與其最大值。公式(10)中其余相關(guān)函數(shù)的表達(dá)式如下所示:

g(ξ)=βexp(-β1ξ)

β=β0(1-ζ/ζmax)

(11)

1.2 本構(gòu)模型的差分形式

由于FLAC3D采用的是顯式三維差分算法，為了更好地實(shí)現(xiàn)程序化，有必要推導(dǎo)黏土廣義塑性模型的差分形式。由式(1)可知總應(yīng)變?cè)隽颗c應(yīng)力增量的差分形式為:

Δσ=De∶Δε-ΔλDe∶ngL/U

(12)

分別用上標(biāo)“N”與“O”表示更新后的應(yīng)力狀態(tài)和更新前的應(yīng)力狀態(tài)，則更新后的應(yīng)力可表示為:

σN=σO+De∶Δε-ΔλDe∶ngL/U=

σI-ΔλDe∶ngL/U

(13)

σI為彈性嘗試獲得的應(yīng)力張量。

由式(13)可知，計(jì)算更新應(yīng)力σN還需解決塑性校正的問題。塑性校正主要涉及塑性流動(dòng)方向向量和塑性乘子的計(jì)算。由公式(8)可知，用于塑性校正的塑性流動(dòng)方向向量表達(dá)式應(yīng)為:

(14)

其次，塑性乘子的計(jì)算。由公式(4)可知塑性乘子的表達(dá)式，結(jié)合式(14)給出的加卸載方向與塑性流動(dòng)方向向量表達(dá)式可得用于塑性校正的塑性乘子為:

(15)

此時(shí)，式(13)轉(zhuǎn)變成如下式所示：

(16)

2 二次開發(fā)流程及其關(guān)鍵點(diǎn)

采用Visual Studio 2019集成開發(fā)環(huán)境(IDE)對(duì)FLAC3D軟件提供的二次開發(fā)模板中的源代碼(.cpp)和頭文件(.h)進(jìn)行修改和編寫，成功將黏土廣義塑性模型嵌入到軟件中。編譯成功的黏土廣義塑性模型會(huì)生成一個(gè)動(dòng)態(tài)鏈接庫(.DLL文件)供主程序通過命令流調(diào)用。調(diào)用用戶自定義模型時(shí)，F(xiàn)LAC3D需要?dú)v經(jīng)兩個(gè)步驟，首先需要通過zone config plugin命令調(diào)用插件模塊，然后通過命令zone cmodel load與完整的動(dòng)態(tài)鏈接庫存儲(chǔ)路徑來實(shí)現(xiàn)用戶自定義模型的調(diào)用。二者缺一不可。

黏土廣義塑性模型的二次開發(fā)流程如圖1所示，出于未考慮模型在循環(huán)加-卸載場景下的應(yīng)用，省略了應(yīng)力加-卸載判斷這一步驟，有必要在后續(xù)研究中進(jìn)行完善。實(shí)現(xiàn)圖1所示的開發(fā)流程，核心技術(shù)在于修改源文件的初始化函數(shù)initialize()和求解函數(shù)run()。其中初始化函數(shù)僅在初始化過程被調(diào)用一次，用于模型參數(shù)和部分中間變量的初始化，以及提示初始賦值過程出現(xiàn)的錯(cuò)誤；求解函數(shù)的功能則是根據(jù)主程序傳入的位移信息(或應(yīng)變?cè)隽啃畔?求解新應(yīng)力。

圖1 程序執(zhí)行流程圖

3 算例驗(yàn)證

3.1 單元試驗(yàn)

單元試驗(yàn)(One zone test)被廣泛用于驗(yàn)證模型二次開發(fā)的正確性[10，13]。本節(jié)通過對(duì)比單元試驗(yàn)結(jié)果、試驗(yàn)成果(Pastor等[1])和理論結(jié)果(Python編寫的計(jì)算機(jī)程序的計(jì)算結(jié)果)，成功證實(shí)了黏土廣義塑性模型二次開發(fā)的正確性。其中，Python程序的計(jì)算結(jié)果等效于有限差分程序中積分點(diǎn)上的計(jì)算結(jié)果，因此，Python程序的計(jì)算結(jié)果與單元試驗(yàn)的模擬結(jié)果應(yīng)基本沒有差別(約10-5量級(jí))。

單元數(shù)值試驗(yàn)針對(duì)Bangkok黏土和Weald黏土開展了固結(jié)排水與不排水三軸應(yīng)力路徑驗(yàn)證工作。兩種土樣的模型參數(shù)均取自文獻(xiàn)[1]，如表1所示。其中Bangkok黏土的前期固結(jié)壓力為414 kPa，超固結(jié)比為24的Weald黏土的前期固結(jié)壓力為120 psi(約827 kPa)。排水與不排水三軸數(shù)值試驗(yàn)的brick單元的邊界條件示意圖如圖2所示，頂部采用速度邊界條件(10-6/時(shí)步)模擬室內(nèi)試驗(yàn)應(yīng)變控制加載的情形，不排水三軸試驗(yàn)還應(yīng)調(diào)用流體計(jì)算模塊，用于流固耦合計(jì)算，水的體積模量取2.2 GPa。注意施加速度邊界條件時(shí)其量值不應(yīng)大于10-5/時(shí)步，否則將影響數(shù)值試驗(yàn)?zāi)M的精度。

表1 單元試驗(yàn)參數(shù)取值

正常固結(jié)Bangkok黏土的不排水和排水三軸試驗(yàn)的數(shù)值模擬結(jié)果、試驗(yàn)結(jié)果以及理論結(jié)果的對(duì)比如圖3—圖6所示。根據(jù)對(duì)比結(jié)果可以發(fā)現(xiàn)FLAC3D模擬結(jié)果與理論結(jié)果完全一致，與試驗(yàn)結(jié)果基本吻合，僅在有效應(yīng)力路徑(見圖4)與排水三軸試驗(yàn)的應(yīng)力比-體積應(yīng)變曲線(見圖6)的模擬結(jié)果上存在細(xì)微差別。圖7與圖8給出了超固結(jié)Weald黏土固結(jié)排水與不排水三軸應(yīng)力應(yīng)變行為的模擬結(jié)果。由結(jié)果可知，該模型較好地反映了超固結(jié)Weald黏土在排水條件下的軟化與剪脹特性，以及不排水條件下的硬化特性，在一定程度上補(bǔ)充了文獻(xiàn)[6]和文獻(xiàn)[7]未開展超固結(jié)土模型研究的不足。

圖2 單元試驗(yàn)示意圖

圖3 Bangkok黏土不排水三軸試驗(yàn)q-εs曲線

圖4 Bangkok黏土不排水三軸試驗(yàn)有效應(yīng)力路徑

圖5 Bangkok黏土排水三軸試驗(yàn)η-εs曲線

圖6 Bangkok黏土排水三軸試驗(yàn)η-εv曲線

3.2 路堤堆載數(shù)值模擬

本節(jié)通過路堤堆載數(shù)值模擬進(jìn)一步驗(yàn)證了黏土廣義塑性模型的工程適用性。

該數(shù)值試驗(yàn)?zāi)M了路堤填筑過程和填筑完成后的地基固結(jié)過程?？紤]到路堤填筑時(shí)長(約27 h)遠(yuǎn)小于地基固結(jié)的持續(xù)時(shí)長(約3 a)，且不涉及路堤填筑過程的研究，故簡化分層填筑過程，采用一個(gè)逐漸增大的均布荷載等效，其最大值為50 kPa，如圖9所示。由于路堤完成填筑的歷時(shí)較短，該過程可近似為不排水過程。

圖7 超固結(jié)Weald黏土q-ε1曲線數(shù)值模擬與試驗(yàn)結(jié)果對(duì)比

圖8 超固結(jié)Weald黏土U-ε1和εv-ε1曲線數(shù)值模擬與試驗(yàn)結(jié)果對(duì)比

圖9 路堤堆載計(jì)算模型示意圖

幾何模型采用30×1×20(x×y×z)的網(wǎng)格模擬20×1×10 m3的地基土體。地基土采用黏土廣義塑性模型(clayPZ)，設(shè)置地表排水。模型參數(shù)取值分別為：壓縮指數(shù)(λ)0.18，回彈指數(shù)(κ)0.05，臨界狀態(tài)應(yīng)力比(Mc)0.8，泊松比(ν)0.3，參考孔隙比(er)1.2，參數(shù)μ、β0、β1、γ、α取值分別為2、0.5、13、1.3、0.5。地基土干重度為20 kN/m3，前期固結(jié)壓力為147 kPa，側(cè)壓力系數(shù)K0為0.64。模型包括兩個(gè)求解過程，一是填筑不排水過程，采用局部不平衡力比值小于或等于10-4作為終止條件；二是地基固結(jié)排水過程，終止該過程的條件需滿足每次迭代的不平衡力比小于5×10-5，且總求解時(shí)長達(dá)到108s(約3 a)。

孔隙水壓力的演化情況如圖10至圖12所示。圖10展示了水平距離為0.5 m處填筑開始前與結(jié)束時(shí)模型預(yù)測的孔隙水壓力隨深度的分布情況。由于設(shè)置地表排水，地表處的孔隙水壓力仍為0 kPa，隨著排水條件影響的減弱，地表附近的孔隙水壓力逐漸增大，在地表以下0.5 m處達(dá)到峰值狀態(tài)(62 kPa)；地面以下0.5 m～1.5 m處脫離了上部荷載的直接影響范圍，超孔隙水壓力迅速減小；地表以下1.5 m～9.5 m的孔隙水壓力受到上部附加荷載的作用孔隙水無法及時(shí)排出，產(chǎn)生超孔隙水壓力，故其比未堆載前的孔隙水壓力大。圖11給出了路堤固結(jié)過程孔隙水壓力在不同深度的分布情況隨時(shí)間的演化?！? d”表示未進(jìn)行路堤填筑前孔隙水壓力沿深度方向的分布。據(jù)分析，填筑結(jié)束時(shí)的孔隙水壓力最大(見圖10)，而后隨著時(shí)間逐漸消散，恢復(fù)至未進(jìn)行填筑前的水平(見圖11)。由于固結(jié)時(shí)間為一年時(shí)的水壓分布已十分接近初始狀態(tài)，為提高插圖的可讀性未展示后續(xù)兩年的孔隙水壓力模擬結(jié)果。圖12給出了深度為1.5 m(z=8.5 m)、6.5 m、9.5 m處孔隙水壓力隨時(shí)間的演化情況，證實(shí)了上述分析，可以看到水壓在1 a時(shí)基本消散完全。

圖10 堆載過程孔隙水壓力的對(duì)比

圖11 固結(jié)過程不同深度孔隙水壓力隨時(shí)間的演化

圖12 路堤堆載全過程孔隙水壓力隨時(shí)間的演化

路堤堆載引起的地面變形(沉降或隆起)是路堤填筑過程值得關(guān)注的問題。如圖13所示，模擬結(jié)果表明堆載結(jié)束時(shí)，在填筑范圍內(nèi)的地面受到上部荷載作用發(fā)生沉降(最大值達(dá)0.24 m)，地基土被擠出，導(dǎo)致周圍地面發(fā)生隆起(峰值達(dá)到0.11 m)。其中堆載范圍內(nèi)地基土可近似為一維固結(jié)，因此其變形量可通過土體壓縮模量(Es)近似計(jì)算。根據(jù)經(jīng)驗(yàn)值可知彈性模量(E)通常為壓縮模量的2～5倍，故可根據(jù)土體的初始孔隙比、回彈指數(shù)、泊松比計(jì)算獲得地面附近1 m處的地基土壓縮模量介于67.2 kPa～168.0 kPa，所以壓縮變形范圍位于0.83 m～2.07 m之間，略大于數(shù)值模擬結(jié)果。其原因?yàn)槁返烫钪^程實(shí)為不排水過程，地基土并未完全固結(jié)，故其沉降量略小于估算值。而后隨著時(shí)間的推移，超孔隙水壓力逐漸消散，有效應(yīng)力增大，土體發(fā)生固結(jié)，分析區(qū)域內(nèi)(0～10 m水平距離)土體均發(fā)生下沉。符合預(yù)期分析的結(jié)果。

圖13 堆載結(jié)束時(shí)的地面沉降分布

4 結(jié) 論

(1) 本文在FLAC3D(7.0版本)提供的二次開發(fā)平臺(tái)上成功實(shí)現(xiàn)了黏土廣義塑性模型的數(shù)值嵌入。并利用單元試驗(yàn)?zāi)M了不同應(yīng)力歷史黏土的應(yīng)力應(yīng)變行為，成功證實(shí)了二次開發(fā)的可靠性與正確性。在一定程度上補(bǔ)充了前人在黏土廣義塑性模型中關(guān)于超固結(jié)土研究的不足(文獻(xiàn)[6-7])。

(2) 通過開展路堤堆載的數(shù)值模擬，驗(yàn)證了二次開發(fā)模型的工程適用性，表明了黏土廣義塑性模型的應(yīng)用價(jià)值，為以后研究飽和黏土相關(guān)問題時(shí)提供了新的工具。