林志遠，熊傳祥，謝永寧，黃盛鋒

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

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

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

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

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

1.1 黏土廣義塑性模型

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

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

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

(1)

式中：dσ，dε分別為應(yīng)力、應(yīng)變增量向量；ngL/U為單位塑性流動方向向量，“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)，如對于向量u和v，(u?v)ij=uivi。K、G分別為彈性體積和剪切模量，其與平均有效應(yīng)力p′的關(guān)系為:

(3)

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

(4)

式中:n為單位加-卸載方向向量；HL為塑性模量。假設(shè)采用相適應(yīng)的流動準則，則有n=ngL/U，在三維應(yīng)力空間中表達式為[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)力空間中的塑性流動方向和加-卸載方向向量的表達式，如下式所示:

(8)

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

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

(9)

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

式(4)中塑性模量HL是一個標量，其表達式為:

(10)

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

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

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

(11)

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

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

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

(12)

分別用上標“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)可知，計算更新應(yīng)力σN還需解決塑性校正的問題。塑性校正主要涉及塑性流動方向向量和塑性乘子的計算。由公式(8)可知，用于塑性校正的塑性流動方向向量表達式應(yīng)為:

(14)

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

(15)

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

(16)

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

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

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

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

3 算例驗證

3.1 單元試驗

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

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

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

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

圖2 單元試驗示意圖

圖3 Bangkok黏土不排水三軸試驗q-εs曲線

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

圖5 Bangkok黏土排水三軸試驗η-εs曲線

圖6 Bangkok黏土排水三軸試驗η-εv曲線

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

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

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

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

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

圖9 路堤堆載計算模型示意圖

幾何模型采用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。模型包括兩個求解過程，一是填筑不排水過程，采用局部不平衡力比值小于或等于10-4作為終止條件；二是地基固結(jié)排水過程，終止該過程的條件需滿足每次迭代的不平衡力比小于5×10-5，且總求解時長達到108s(約3 a)。

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

圖10 堆載過程孔隙水壓力的對比

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

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

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

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

4 結(jié) 論

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

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