王子珺，趙伯明

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

相關(guān)研究表明砂土的變形性質(zhì)依賴于砂土密度以及剪切過程中所受的圍壓大小[5-7]，具體表現(xiàn)為松砂具有剪縮特性而中密砂則具有剪脹特性。為統(tǒng)一反映砂土密度和所受圍壓對(duì)砂土變形特性的影響，Bean 和Jefferies[5]首先提出通過一個(gè)狀態(tài)參數(shù)來反映這兩者的耦合作用。狀態(tài)參數(shù)ψ =e-ec，即當(dāng)前有效平均正應(yīng)力下，實(shí)際孔隙比e與臨界狀態(tài)孔隙比ec之差。臨界狀態(tài)線定義為剪切過程中砂土的體應(yīng)變?cè)隽繛?，而偏應(yīng)變?cè)隽繛闊o窮大的一個(gè)狀態(tài)。對(duì)于砂土的本構(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]。當(dāng)ψ ＞0時(shí)，當(dāng)前實(shí)際孔隙比大于相同壓力下的臨界孔隙比，材料處于松散狀態(tài)，如圖1 中A 區(qū)域所示；當(dāng) ψ ＜0時(shí)，當(dāng)前實(shí)際孔隙比小于相同壓力下的臨界孔隙比，材料處于密實(shí)狀態(tài)，如圖1 中B 區(qū)域所示。 ψ值越大，材料越松散， ψ值越小，材料越密實(shí)。圖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è)砂土的本構(gòu)模型，通過一組參數(shù)統(tǒng)一考慮不同初始密度與不同圍壓條件下砂土的應(yīng)力-應(yīng)變關(guān)系，但是該模型僅限于二維p′-q平面問題的求解。因此，本文在該方法的基礎(chǔ)上，將模型推廣應(yīng)用于三維p′-q-θ空間問題的求解，使之可以考慮應(yīng)力洛德角 θ的影響，從而反映三維應(yīng)力狀態(tài)下土的變形和強(qiáng)度特性。為了能夠利用有限元軟件前后處理方便、計(jì)算穩(wěn)定性與精度高等優(yōu)點(diǎn)，通過ABAQUS 提供的用戶自定義材料子程序UMAT，基于Fortran 語言編制了上述改進(jìn)的三維砂土模型的接口程序，實(shí)現(xiàn)該彈塑性本構(gòu)關(guān)系，進(jìn)而分別通過3 種砂，即Toyoura 砂、Fuji River 砂和Tokachi 砂的三軸排水剪切試驗(yàn)結(jié)果對(duì)模型的有效性進(jìn)行對(duì)比驗(yàn)證。研究成果可進(jìn)一步應(yīng)用于相關(guān)巖土工程的數(shù)值分析，為巖土工程分析提供更加快捷實(shí)用的解決方案。

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

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

1.1 彈性性質(zhì)

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

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

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

1.3 塑性流動(dòng)法則

1.4 加載塑性模量

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

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

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

2.2 三軸試驗(yàn)?zāi)M與驗(yàn)證

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

圖3 計(jì)算模型Fig. 3 Computational model

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

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

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

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

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

圖5 p′0=500 kPa 下Toyoura 砂三軸排水剪切試驗(yàn)結(jié)果與有限元計(jì)算結(jié)果對(duì)比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

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

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

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

圖7 密實(shí)狀態(tài)下Fuji River 砂三軸排水剪切試驗(yàn)結(jié)果(圖例符號(hào))與有限元計(jì)算結(jié)果(曲線)對(duì)比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)變曲線與試驗(yàn)結(jié)果吻合程度很高，表明有限元計(jì)算結(jié)果可以反映圍壓和砂土初始密度對(duì)應(yīng)力-應(yīng)變的影響，模型能夠準(zhǔn)確描述砂土體積的變化趨勢(shì)。

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

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

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

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

3 結(jié)論

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

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

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