砂土統一本構模型研究及其三維數值實現

王子珺，趙伯明

(1. 北京交通大學城市地下工程教育部重點實驗室，北京 100044；2. 北京交通大學土木建筑工程學院，北京 100044)

相關研究表明砂土的變形性質依賴于砂土密度以及剪切過程中所受的圍壓大小[5-7]，具體表現為松砂具有剪縮特性而中密砂則具有剪脹特性。為統一反映砂土密度和所受圍壓對砂土變形特性的影響，Bean 和Jefferies[5]首先提出通過一個狀態參數來反映這兩者的耦合作用。狀態參數ψ =e-ec，即當前有效平均正應力下，實際孔隙比e與臨界狀態孔隙比ec之差。臨界狀態線定義為剪切過程中砂土的體應變增量為0，而偏應變增量為無窮大的一個狀態。對于砂土的本構模型，相關研究[8]表明其e-p′平面的臨界狀態線更適合于非線性形式描述，即：

式中：eΓ和λc分別為臨界狀態線的截距和斜率；pˉ=p′/pa為有效應力與大氣壓之比； ξ為常數，上述模型也得到相關學者的采納與應用[3,9-11]。當ψ ＞0時，當前實際孔隙比大于相同壓力下的臨界孔隙比，材料處于松散狀態，如圖1 中A 區域所示；當 ψ ＜0時，當前實際孔隙比小于相同壓力下的臨界孔隙比，材料處于密實狀態，如圖1 中B 區域所示。 ψ值越大，材料越松散， ψ值越小，材料越密實。圖2 為根據式(2)所建立的幾種不同砂 土(Toyoura 砂、Tung-Chung 砂、Fuji River 砂)的臨界狀態線，其中， ξ建議取為0.7[8,12]。

圖1 狀態參數與臨界狀態線Fig. 1 State parameter and critical state line

圖2 幾種不同砂土的臨界狀態線Fig. 2 Critical state lines of several different sandy soils

基于廣義塑性理論與臨界狀態概念，Ling 和Yang[4]建立了一個砂土的本構模型，通過一組參數統一考慮不同初始密度與不同圍壓條件下砂土的應力-應變關系，但是該模型僅限于二維p′-q平面問題的求解。因此，本文在該方法的基礎上，將模型推廣應用于三維p′-q-θ空間問題的求解，使之可以考慮應力洛德角 θ的影響，從而反映三維應力狀態下土的變形和強度特性。為了能夠利用有限元軟件前后處理方便、計算穩定性與精度高等優點，通過ABAQUS 提供的用戶自定義材料子程序UMAT，基于Fortran 語言編制了上述改進的三維砂土模型的接口程序，實現該彈塑性本構關系，進而分別通過3 種砂，即Toyoura 砂、Fuji River 砂和Tokachi 砂的三軸排水剪切試驗結果對模型的有效性進行對比驗證。研究成果可進一步應用于相關巖土工程的數值分析，為巖土工程分析提供更加快捷實用的解決方案。

1 基于廣義塑性理論和臨界狀態的砂土彈塑性模型

為了將狀態參數引入廣義塑性理論，對建立模型所需要的主要塑性參數，包括：剪脹性、塑性流動方向、加載方向和塑性模量等進行相應修正。

1.1 彈性性質

式中，G0和K0分別為彈性剪切模量與彈性體積模量。

1.2 剪脹性與相變轉換

顆粒材料的剪脹性是砂土體積與形狀變化之間的耦合現象[16]，應力剪脹方程描述了塑性體應變和塑性剪應變之間的關系：

1.3 塑性流動法則

1.4 加載塑性模量

2 砂土本構模型的二次開發與驗證

2.1 用戶自定義材料子程序UMAT 與坐標轉換

針對非線性問題，通過ABAQUS 軟件提供的用戶自定義材料子程序UMAT(Use-defined-Material Mechanical Behaviour)與 ABAQUS/Standard 主求解程序的接口對接實現與ABAQUS 的數據交換，完成有限元數值計算。UMAT 子程序通過ABAQUS 主程序傳入應變增量以及當前已知狀態的應力、應變及其他與求解過程有關的變量。基于本構方程，通過迭代計算獲得真實應力以及彈性、塑性應變，從而更新應力增量和狀態變量，通過DDSDDE 數組提供材料本構模型的雅克比(Jacobian)矩陣 ?Δσ/?Δε，即應力增量對應變增量的變化率[19-20]，供主程序進行Newton-Raphson非線性迭代。

2.2 三軸試驗模擬與驗證

基于ABAQUS 提供的UMAT 子程序接口，使用Fortran 語言編制了上述改進的三維砂土本構模型，通過有限元模擬三軸試驗并與3 種砂(即Toyoura 砂[21]、Fuji River 砂[22-23]以及Tokachi 砂[24])的試驗結果進行對比，以驗證模型的有效性。模擬的有限元模型如圖3 所示，計算網格采用C3D20RP型。由于本研究旨在模擬驗證常規三軸的應力-應變關系，因此，無需建立圓柱形試樣模型[25]。

圖3 計算模型Fig. 3 Computational model

提出的模型共需要12 個參數，分別為彈性參數(G0、K0)、臨界狀態線參數(eΓ、λc)、剪脹性參數( α、Mg、mg)、塑性流動參數(Mf、mf)、加載模量參數(m0、HL0、mb)。對于Toyoura 砂和Fuji River 砂，參數G0、K0、eΓ、λc、Mg、HL0、α由文獻[4]獲取；對于Tokachi 砂，上述參數由文獻[24]獲取，其余各項參數可分別通過第2 節中相關公式求得。此外，單位大氣壓力pa取為100 kPa。因此，上述UMAT 中的模型參數數組PROPS 共包括12 個分量，Toyoura 砂、Fuji River砂以及Tokachi砂的參數值見表1。

表1 統一廣義塑性模型參數Table 1 Parameters of unified generalized plasticity model

由于用于存儲與解有關的狀態變量的數組需要不斷更新，將每一步的結果作為下一步的初值，因此，還要設置相應的狀態變量數組STATEV。在初始分析步(Initial)后定義2 個分析步，Geostatic(地應力)和Load(加載)。其中，Geostatic 分析步賦予初始圍壓，Load 分析步模擬施加軸向荷載過程。分別約束底面(1-2-3-4)以及側面(1-4-8-5)、(1-5-6-2)3 個面上的法向平移自由度以及另外兩個方向的旋轉自由度，對其它3 個面施加初始圍壓。然后在模型頂面按照位移控制施加豎向荷載，直至模型軸向應變達到ε=25%。采用ABAQUS/Standard中自動劃分時間增量步長的方法，可以為UMAT輸送較為合理的 Δε值。最后提交任務，讀取包含UMAT 子程序的Fortran 文件，進行后處理。

根據調整的模型參數結合有限元計算結果對砂土彈塑性本構模型二次開發的準確性和有效性進行驗證。對于Toyoura 砂，圖4 和圖5 分別表示利用100 kPa 和500 kPa 下不同孔隙比(e=0.831～0.996；e=0.810～0.990)的三軸排水剪切試驗結果(圖例符號表示)與模擬計算結果(曲線表示)的對比。

圖4 p′0=100 kPa 下Toyoura 砂三軸排水剪切試驗結果(圖例符號)與有限元計算結果(曲線)對比Fig. 4 Comparison between experimental results (symbols)and finite element analysis results (curves) on drained triaxial responses of Toyoura sand with initial confining pressure p′0=100 kPa

由圖4 和圖5 可知，根據建立的ABAQUS 有限元模型得到的曲線與試驗結果吻合度很高。對于砂土的排水剪切試驗，其力學特性因初始孔隙比和圍壓不同而有所變化。在相同的圍壓下，當初始孔隙比較大時(e=0.996/0.990)，在剪切過程中試樣一直表現為剪縮狀態，應力-應變曲線一直表現為硬化；當砂土的初始孔隙比較小時(e=0.831/0.810)，試樣表現為先剪縮后剪脹，且剪縮量較小而剪脹量較大，應力-應變曲線表現為初始硬化，達到峰值強度以后表現為少量的軟化；對于中密度砂(e=0.917/0.886)，應力-應變關系則介于兩者之間。盡管試樣的初始孔隙比不同，隨著應變的不斷增加，應力最終趨于穩定值。試驗結果與預測結果對比顯示，模型能夠很好地預測砂土體積變化趨勢，反映出砂土的力學特性。

圖5 p′0=500 kPa 下Toyoura 砂三軸排水剪切試驗結果與有限元計算結果對比Fig. 5 Comparison between experimental results (symbols)and finite element analysis results (curves) on drained triaxial responses of Toyoura sand with initial confining pressure p′0=500 kPa

對于Fuji River 砂，圖6 和圖7 分別表示利用松散狀態(e=0.747～ 0.751)和密實狀態(e=0.496～0.519)下三軸排水剪切試驗結果(圖例符號表示)與模型結果(曲線表示)的對比。

對比圖6 和圖7 可知，對于松砂，當初始圍壓較大時(p=294 kPa)，試樣一直表現為剪縮；當初始圍壓較小時(p=98 kPa)，試樣表現為先剪縮后剪脹；當初始圍壓p=196 kPa 時，試樣變形介于二者之間。對于密砂，在考察的三種圍壓條件下均表現出先剪縮后剪脹性質，初始圍壓越小的試樣剪脹量越大；隨著圍壓的增大，砂土的強度和剛度都明顯提高，密砂的剪切特性有像松砂轉變的趨勢。

圖6 松散狀態下Fuji River 砂三軸排水剪切試驗結果(圖例符號)與有限元計算結果(曲線)對比Fig. 6 Comparison between experimental results (symbols)and finite element analysis results (curves) on drained triaxial responses of loose Fuji river sand

圖7 密實狀態下Fuji River 砂三軸排水剪切試驗結果(圖例符號)與有限元計算結果(曲線)對比Fig. 7 Comparison between experimental results (symbols)and finite element analysis results (curves) on drained triaxial responses of dense Fuji river sand

同樣，根據建立的ABAQUS 有限元模型得到的偏應力與軸向應變曲線以及體應變與軸向應變曲線與試驗結果吻合程度很高，表明有限元計算結果可以反映圍壓和砂土初始密度對應力-應變的影響，模型能夠準確描述砂土體積的變化趨勢。

對于Tokachi 砂，圖8 和圖9 分別表示利用70 kPa 和30 kPa 下不同孔隙比(e=1.01 和e=0.832)三軸排水剪切試驗結果(圖例符號表示)與模型結果(曲線表示)的對比。

圖8 中密Tokachi 砂(e=1.01)三軸排水剪切試驗結果(圖例符號)與有限元計算結果(曲線)對比Fig. 8 Comparison between experimental results (symbols)and finite element analysis results (curves) on drained triaxial responses of medium dense Tokachi sand (e=1.01)

圖9 密實Tokachi 砂(e=0.832)三軸排水剪切試驗結果(圖例符號)與有限元計算結果(曲線)對比Fig. 9 Comparison between experimental (symbols) and finite element analysis results (curves) on drained triaxial responses of dense Tokachi sand (e=0.832)

同樣，模擬結果能夠很好地反映試驗結果。當初始孔隙比較大時(e=1.01)，對于所受初始圍壓較大的試樣，在剪切過程中一直表現為剪縮，而對于所受初始圍壓較小的試樣則表現為先剪縮后剪脹；當初始孔隙比較小時(e=0.832)，兩種圍壓下的試樣均表現出一定程度的先剪縮后剪脹現象。

上述結果表明有限元計算結果可以反映圍壓和砂土初始密度對應力-應變曲線的影響。由于砂土應力-應變之間不成單值函數關系，因此，反映土體的數學模型一般形式復雜難于準確，而本文基于現有大型有限元軟件ABAQUS 平臺進行二次開發，可以有效利用其強大的非線性有限元分析功能及優秀的前后處理界面，使開發的砂土統一本構模型能夠進一步應用于相關巖土工程的數值分析，為所涉及的復雜工程問題提供更加快捷實用的解決方案。

3 結論

本文基于廣義塑性理論與臨界狀態概念，研究提出了一個統一三維砂土本構模型，工作要點及結論如下：

(1)該模型可通過一組參數實現砂土由壓縮至剪切過程中狀態參量的統一表述，且這些參數均能夠通過常規土工試驗獲得。基于ABAQUS 提供的用戶自定義材料子程序UMAT 接口實現了該彈塑性本構關系，同時可以結合軟件自動步長搜索功能，具有很好的穩定性和較高的計算精度。

(2)分別利用Toyoura 砂、Fuji River 砂與Tokachi 砂的剪切試驗與模型模擬結果進行對比驗證，結果表明有限元計算模型可以有效反映不同砂土類型在不同圍壓和砂土初始密度條件下的應力-應變曲線關系，能夠準確描述密砂的剪脹特性與應變軟化特性以及松砂的剪縮特性與應變硬化特性，從而更加真實地反映三維應力狀態下土的變形和強度特性。研究成果可進一步應用于相關巖土工程的數值分析，為巖土工程分析提供更加快捷實用的解決方案。