王子珺，趙伯明

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

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

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

圖1 狀態(tài)參數(shù)與臨界狀態(tài)線Fig. 1 State parameter and critical state line

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

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

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

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

1.1 彈性性質(zhì)

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

1.2 剪脹性與相變轉(zhuǎn)換

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

1.3 塑性流動法則

1.4 加載塑性模量

2 砂土本構(gòu)模型的二次開發(fā)與驗證

2.1 用戶自定義材料子程序UMAT 與坐標轉(zhuǎn)換

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

2.2 三軸試驗?zāi)M與驗證

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

圖3 計算模型Fig. 3 Computational model

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

表1 統(tǒng)一廣義塑性模型參數(shù)Table 1 Parameters of unified generalized plasticity model

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

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

圖4 p′0=100 kPa 下Toyoura 砂三軸排水剪切試驗結(jié)果(圖例符號)與有限元計算結(jié)果(曲線)對比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 可知，根據(jù)建立的ABAQUS 有限元模型得到的曲線與試驗結(jié)果吻合度很高。對于砂土的排水剪切試驗，其力學(xué)特性因初始孔隙比和圍壓不同而有所變化。在相同的圍壓下，當初始孔隙比較大時(e=0.996/0.990)，在剪切過程中試樣一直表現(xiàn)為剪縮狀態(tài)，應(yīng)力-應(yīng)變曲線一直表現(xiàn)為硬化；當砂土的初始孔隙比較小時(e=0.831/0.810)，試樣表現(xiàn)為先剪縮后剪脹，且剪縮量較小而剪脹量較大，應(yīng)力-應(yīng)變曲線表現(xiàn)為初始硬化，達到峰值強度以后表現(xiàn)為少量的軟化；對于中密度砂(e=0.917/0.886)，應(yīng)力-應(yīng)變關(guān)系則介于兩者之間。盡管試樣的初始孔隙比不同，隨著應(yīng)變的不斷增加，應(yīng)力最終趨于穩(wěn)定值。試驗結(jié)果與預(yù)測結(jié)果對比顯示，模型能夠很好地預(yù)測砂土體積變化趨勢，反映出砂土的力學(xué)特性。

圖5 p′0=500 kPa 下Toyoura 砂三軸排水剪切試驗結(jié)果與有限元計算結(jié)果對比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 分別表示利用松散狀態(tài)(e=0.747～ 0.751)和密實狀態(tài)(e=0.496～0.519)下三軸排水剪切試驗結(jié)果(圖例符號表示)與模型結(jié)果(曲線表示)的對比。

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

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

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

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

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

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

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

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

3 結(jié)論

本文基于廣義塑性理論與臨界狀態(tài)概念，研究提出了一個統(tǒng)一三維砂土本構(gòu)模型，工作要點及結(jié)論如下：

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

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