黏土廣義塑性模型在FLAC3D中的數(shù)值實現(xiàn)與應用

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

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

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

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

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

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

1.1 黏土廣義塑性模型

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

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

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

(1)

式中：dσ，dε分別為應力、應變增量向量；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為單位矩陣；?為克羅內克積(Kronecker product)，如對于向量u和v，(u?v)ij=uivi。K、G分別為彈性體積和剪切模量，其與平均有效應力p′的關系為:

(3)

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

(4)

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

(5)

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

(6)

(7)

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

(8)

上式中，d為剪脹性關系:

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

(9)

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

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

(10)

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

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

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

(11)

1.2 本構模型的差分形式

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

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

(12)

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

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

σI-ΔλDe∶ngL/U

(13)

σI為彈性嘗試獲得的應力張量。

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

(14)

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

(15)

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

(16)

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

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

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

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

3 算例驗證

3.1 單元試驗

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

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

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

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

圖2 單元試驗示意圖

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

圖4 Bangkok黏土不排水三軸試驗有效應力路徑

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖13 堆載結束時的地面沉降分布

4 結 論

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

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