改進超固結狀態參量次加載面模型數值實現與應用

馬少坤，趙乃峰，2，潘柏羽，周 東，江 杰

(1. 廣西大學 土木建筑工程學院；工程防災與結構安全重點實驗室，南寧 530004；2. 同濟大學 地下建筑與工程系，上海 200092)

?

改進超固結狀態參量次加載面模型數值實現與應用

馬少坤1，趙乃峰1，2，潘柏羽1，周 東1，江 杰1

(1. 廣西大學 土木建筑工程學院；工程防災與結構安全重點實驗室，南寧 530004；2. 同濟大學 地下建筑與工程系，上海 200092)

為合理描述超固結土復雜的彈塑性力學行為，對現有Hashiguchi次加載面模型中的超固結狀態參量R進行修正，在硬化方程中，考慮塑性體應變與塑性剪應變的綜合作用，提出了修正超固結狀態參量的次加載面模型。同時，著重介紹了該模型的隱式積分算法及數值實現過程，編制了對應的接口子程序，實現了該模型在有限元軟件ABAQUS中的應用。通過不同工況和加載方式下的數值模擬驗證了程序的合理性，最后應用模型研究了Fujinomori 黏土的三軸壓縮力學特性并與UH模型的模擬結果、室內試驗研究進行對比。結果表明，子程序具有較高的計算精度和可靠性，模型能夠準確地模擬黏土的超固結特性。

次加載面模型；超固結土；ABAQUS；用戶子程序

在隧道、基坑等大型地下工程的開挖過程中，開挖面土體因所處位置不同而經受不同的加載、卸載等復雜應力路徑，應力狀態發生變化，由正常固結狀態轉變為超固結狀態。Mesri[1]、Tavenas等[2]及眾多巖土力學學者研究表明，自然沉積土多數為結構性土且處于超固結狀態。超固土體具有剪脹、應變軟化等特性，其應力應變關系與正常固結土存在較大差異。以劍橋模型[3-4]為代表的經典彈塑性理論認為：巖土介質只存在一個屈服面，土在卸載再加載過程中，應力應變關系為彈性，如圖1(a)所示。而實際上，正常固結土一旦卸載就處于超固結狀態，土在超固結狀態下的再加載過程仍具有塑性變形，如圖1(b)所示[5]。

為準確反映超固結土的以上特性，許多學者提出了不同的模型。如Krieg[6]、Dafatias等[7]提出的二面模型，Dafatias等[8]提出的邊界面模型，Hashiguchi等[9-12]提出的次加載面模型等。次加載面模型認為超固結狀態下無純彈性域[13]，塑性屈服由次加載面控制，該面始終過當前應力狀態點，隨加卸載擴大或縮小，以此來描述超固結狀態下的塑性應變。張鋒等[5]基于Nakai等[14]提出的土的密度的概念，建立了一個新的超固結重塑黏土的次加載面劍橋模型；Asaoka等[15-16]基于原始劍橋模型屈服準則，提出了可以考慮土的結構性和超固結性的上加載面模型，但模型中超固結狀態參數R和結構狀態參數R*的定義模糊，限制了模型的推廣與應用。

圖1 經典彈塑性理論土與實際土應力應變關系[5]Fig.1 Relationships of stress-strain in critical state theory and reality[5]

如上所述，現有次加載面模型理論復雜、形式多樣，在有限元軟件的本構模型庫中鮮有涉及，使該模型的應用與推廣受到極大限制。首先基于Hashiguchi[9-12]、Yamakawa等[17]的次加載面劍橋模型，改進超固結狀態變量R的演化規則，在硬化方程中考慮塑性體應變與塑性剪應變對超固結狀態參量R的共同作用，再利用隱式積分算法編制與改進模型對應的接口子程序，實現對有限元軟件ABAQUS的二次開發，建立適用于超固結土應力應變關系的數值模擬平臺；最后通過不同工況和加載方式的數值模擬，對Fujinomori 黏土數值模擬與室內試驗結果、UH模型模擬結果對比驗證程序的合理性，為該模型的實際工程應用奠定基礎。

1 次加載面修正劍橋模型

1.1 次加載面概念

次加載面理論認為超固結土存在兩個應力狀態面：正常固結屈服面和次加載面，如圖2所示。正常固結屈服面為修正劍橋模型的屈服面，大小由土在其固結歷史上所經受的最大應力水平來確定。次加載面[18-19]為土體卸荷至某一應力狀態時通過此應力狀態點與正常固結屈服面幾何相似的面，相似比為R，相似中心為p-q應力空間的原點。次加載面位于正常固結屈服面內，隨應力狀態的變化擴大或縮小。卸載時，當前應力狀態點遠離正常固結屈服面，R減小，次加載面縮小；加載時，當前應力狀態點向正常固結屈服面靠近，R增大，次加載面擴大，若當前應力狀態點處于正常固結屈服面上，則R=1，次加載面與正常固結面重合，土由超固結狀態轉變為正常固結狀態。

圖2 p-q空間中的次加載面示意圖Fig.2 Sketch of subloading surface in p-q space

1.2 屈服面方程與塑性勢函數

引入描述超固結程度的超固結比(OCR)概念，給出本文模型兩個應力狀態面相似比(R)的定義式：

(1)

式中：pc為次加載面大小；pnc為正常固結屈服面大小。

結合正常固結屈服面方程(修正劍橋模型屈服面方程[2])及式(1)可得次加載面在p-q應力空間的方程：

(2)

為準確描述土及軟巖的力學特征，近代土力學尤其是次加載面理論建議采用相關聯流動法則[3]，因此，式(2)即為本文模型的塑性勢函數表達式。

確定了塑性勢函數后，塑性應變率即可由以下流動法則給出：

(3)

1.3 協調方程

由式(2)等號兩邊取微分，得到滿足塑性一致條件的協調方程式(4)。

(4)

式中dpnc可用修正劍橋模型的硬化規律表達：

(5)

對式(4)中描述超固結程度的狀態變量R的微分(dR)而言，Hashiguchi[11]、Asaoka等[15-16]建議：

(6)

式中：mR為材料參數，其大小表示超固結狀態隨塑性體應變增加而消失的快慢。

式(6)中dR僅與塑性體應變相關，沒有考慮到塑性剪應變對超固結性的影響。考慮到超固結性的發展受塑性體應變和塑性剪應變綜合影響的情況，本文綜合Asaoka等[15-16]、Yamakawa等[17]的建議式，把塑性體應變增量與塑性剪應變增量對超固結性的影響做加權，得出本文模型dR的表達式

(7)

(7′)

式中：η為非負的材料參數，表示塑性剪應變增量對超固結性的發展貢獻比。

2 模型應力更新算法

圖3 完全隱式積分算法示意圖[20]Fig.3 Sketch of implicit integration algorithm[20]

2.1 初始變量計算

由狀態變量σN、RN及材料參數M可得：

(8)

式中：pc,N、pnc,N分別表示N增量步時次加載面和正常固結屈服面的大小。

2.2 彈性試算

(9)

模型彈性計算采用多孔介質非線性彈性，體積彈性模量和剪切模量分別為

(10)

式中：v為泊松比。

2.3 初始屈服判斷

首先由試探應力分量σtr計算試探平均應力ptr、試探剪切應力qtr，進而計算初始屈服函數ftr：

(11)

若屈服函數f小于某個容許誤差值ftol，則應力狀態處于彈性階段，進行下文2.6節的處理，否則進行塑性修正。本文設定ftol=1×10-5。

2.4 塑性修正

(1)更新第k步迭代時φ的值φk：

其中：

2.5 一致切線模量

2.6 變量更新與存儲

由下式進行應變、孔隙比等的更新，并進行狀態變量STATEV(NSTATV)的存儲：

(12)

3 程序驗證

為了驗證本文算法的可靠性與程序的精度，選取文獻[22]中所列土樣，采用一階八節點三維實體孔壓單元(C3D8P)，對試樣進行不同工況下不同加載方式的數值模擬，并進行結果比較。限于篇幅，本文僅列出正常固結(NC)和超固結(OC)工況下三軸排水壓縮(CD)、三軸不排水壓縮(CU)兩種加載方式的實驗結果。材料參數如表1所示。

表1 材料參數[22]

式(7)中的材料參數mR在文獻[5,15-16]、文獻[23]中的取值在[2, 10]之間，經驗證，當mR取8.0時，本文模型的模擬結果更合理；對式(7′)中與超固結狀態參數R的發展有關的參數η，根據本文模擬，取0.8較合適。

模擬過程為：在初始分析步中限定模型土樣底部X、Y兩個方向的位移，給模型土樣施加圍壓并在以后的分析步驟中圍壓保持不變；在荷載步驟中給模型土樣施加軸向位移15.2 mm(軸向應變的20%)。

圖5為本文算法與修正劍橋模型結果的應力路徑圖。如圖所示，固結過程中正常固結土的應力路徑從“濕側區域”達到臨界狀態線，超固結土的應力路徑從“干側區域”穿過臨界狀態線后又返回到臨界狀態線上，反映了土固結的一般規律。對比本文算法、修正劍橋模型的應力路徑(圖5(a))可知：在模擬正常固結土的固結過程時，本文算法與修正劍橋模型有著極強的一致性，說明本文算法能夠準確描述正常固結土固結過程的應力狀態變化；在模擬超固結土的固結過程時，本文算法的應力路徑更平滑，與實際中平滑過渡的應力路徑更接近。對比本文算法、Sheng等[22]算法的應力路徑(圖5(b))，可以明顯看出，在模擬正常固結土不排水固結過程時，本文算法模擬結果的應力路徑與實際平滑的應力路徑更接近。說明本文算法相對修正劍橋模型與Sheng等[22]算法，能更準確地描述實際正常固結土和超固結土固結過程中的應力狀態變化。

圖5 有效應力路徑圖Fig.5 Elevation Effective stress path diagram

圖6為關于孔隙比與有效應力關系的本文算法模擬結果(標記為UMAT)與文獻[22]結果(標記為Sheng)的對比圖。由圖可知，正常固結土(NC)在不排水固結(CU)過程中孔隙比保持不變，在排水固結(CD)過程中孔隙比隨固結壓力增大而減小，表現出剪縮性；超固結土(OC)在不排水固結(CU)過程中孔隙比保持不變，在排水固結(CD)過程中孔隙比隨固結壓力增大先減小后增大，表現出剪脹性，在一定程度上反應了超固結土的應變軟化特性。以上特性與臨界狀態理論完全一致，范慶來[24]也曾得出相似的結論，可見本文算法能夠較好地描述土的減縮、剪脹及軟化特征。

圖6 孔隙比與平均有效應力關系Fig.6 Relationship between porosity ratio and mean effective stress

圖7 超固結土偏應力軸向應變關系曲線Fig.7 Relationships between deviatoric stress and axial strain of overconsolidated soils

綜上所述，本文算法與修正劍橋模型的模擬結果有著良好的一致性，且能更準確地描述超固結土的彈塑性力學行為，說明本文算法合理，所編子程序正確。

4 本文模型應用

應用本文模型，對Fujinomori 黏土進行三軸壓縮試驗模擬，并與該土樣的室內試驗結果[25]、UH模型模擬結果[26]進行比較。土樣參數如表2所示，具體試驗見文獻[25]。

表2 Fujinomori clay主要材料參數[25]

表2中的材料參數Rcs為三軸壓縮過程破壞時的極限主應力比(σ1/σ3)cs，由文獻[25]中的公式換算為本文模型的參數M=1.364。

圖8 試驗結果與模擬結果對比圖Fig.8 Comparison between Simulation and experience

5 結論

建立了改進超固結狀態參量的次加載面模型，通過編制該模型對應的用戶子程序，實現了對有限元軟件的二次開發，建立了適用于超固結土的數值模擬平臺。隨后，應用模型分析了不同工況和加載條件下土的力學特性，并與UH模型的數值模擬結果及試驗結果進行了對比分析，得如下結論：

3)在描述超固結土的力學行為方面，本文模型相對UH模型能夠更準確地刻畫超固結土的應變軟化、峰值強度與殘余強度值以及應力應變關系等特征。

[1] Mesri G. Discussion of “New design procedure for stability of soft clays” [J]. Journal of the Geotechnical Engineering Division，ASCE，1975，1：409-412.

[2] Tavenas F，Leroueil S. Laboratory and in situ stress-strain-time behavior of soft clays：A state of the art [C]//International Symposium on Geotechnical Engineering of Soft Soils，Mexico City，1987，2：3-48.

[3] Roscoe K H，Schofield A N，Worth C P. On the yielding of soils [J]. Geotechnique，1958，8(1)：22-53.

[4] Roscoe K H，Schofield A N and Thuraira A. Yielding of clay in states wetter than critical [J]. Geotechnique，1963，13(3)：211-240.

[5] 張鋒，葉冠林. 計算土力學[M]. 北京：人民交通出版社， 2007.

[6] Krieg R D. A practical two surface plasticity theory [J].Journal of Applied Mechanics，ASME. 1975，42：641-646.

[7] Dafatias Y F，Popov E P. A model of nonlinearly hardening materials for complex loading [J]. Acta Mechanical，1975，23(3)：173-192.

[8] Dafatias Y F，Herrmamn L R. A bounding surface soil plasticity model [C]//International Symposium on Soils under Cyclic Trans. Swansea，Balkema，1980：335-345.

[9] Hashiguchi K，Ueno M. Elastoplastic constitutive laws of granular material [C]//Constitutive Equations of Soils，Tokyo，JSSMFE，1977：73-82.

[10] Hashiguchi K. Constitutive equations of elastoplastic materials with elastoplastic transition [J]. Journal of Applied Mechanics，1980，47(2)：266-272.

[11] Hashiguchi K. Subloading surface model in unconventional plasticity [J]. International Journal of Solids and Structures，1989，25(8)：917-945.

[12] Hashiguchi K，Chen Z P. Elastoplastic constitutive equation of soils with the subloading surface and rotational hardening [J]. International Journal for Numerical and Analytical Methods in Geomechanics，1998，22：197-227.

[13] 孔亮，鄭穎人，姚仰平. 基于廣義塑性力學的土體次加載面循環塑性模型(Ⅰ)：理論與模型[J]. 巖土力學，2003，24(2)：141-145.

Kong L，Zheng Y R，Yao Y P. Subloading surface cyclic plastic model for soil based on Generalizedplasticity(Ⅰ)：Theory and model [J]. Rock and Soil Mechanics，2003，24(2)：141-145.(in Chinese)

[14] Nakai T，Horii H. Plane strain compression test of cemented sand and measurement of strain localization [C]//Proc. Of 28th Annual Conference of Japan Geotechnical Society，1993：545-548.(in Japanese)

[15] Asaoka A，Nakano M，Noda T. Superloading yield surface concept for highly structured soil behavior [J]. Soils and Foundations，2000，40(2)：99-110.

[16] Asaoka A，Nakano M，Noda T，et al. Delayed compression/consolidation of naturally clay due to degradation of soil structure [J].Soils and Foundations，2000，40(3)：75-85.

[17] Yamakawa Y，Hashiguchi K，Ikeda K. Implicit stress-update algorithm for isotropic Cam-clay model based on the subloading surface concept at finite strains [J]. International Journal of Plasticity，2010，26(5)：634-658.

[18] Hashiguchi K. Plastic constitutive equations of granular material [C]//US-Japan Seminar on Continuum Mechanics and Statistical Approaches in the Mechanics of Granular Materials. Sendai，JSSMFE，1978：321-329.

[19] Asaoka A，Nakano M，Noda T. Soil-water coupled behavior of heavily overconsolidated clay near/at critical state [J]. Soils and Foundations，1997，37(1)：13-28.

[20] 黃雨，周子舟. 下負荷面劍橋模型在ABAQUS中的開發實現[J].巖土工程學報，2010，32(1)：115-119.

Huang Y，Zhou Z Z. Numerical implementation for subloading Cam-Clay model in ABAQUS [J].Chinese Journal of Geotechnical Engineering，2010，32(1)：115-119.(in Chinese)

[21] Borjar I，Lee S R. Cam-clay plasticity，Part 1：implicit integration of elasto-plastic constitutive relations [J]. Computers Methods in Applied Mechanics and Engineering，1990，78：49-72．

[22] Sheng D，Sloan S W，Yu H S. Aspects of finite element implementation of critical state models [J]. Computational Mechanics，2000,(6)：185-196.

[23] Zhu H H，Ye B，Cai Y C，et al. An elastoviscoplastic model for soft rock around tunnels considering overconsolidation and structure effects [J]. Computers and Geotechnics，2013，50：6-16.

[24] 范慶來，欒茂田，倪宏革. 循環荷載作用下軟基上大圓筒結構彈塑有效應力分析[J]. 水利學報，2008，39(7)：836-842.

Fan Q L，Luan M T，Ni H G. Elastoplastic effective stress analysis of soft soil foundation of large-diameter cylindrical structure subjected to cyclic loading [J]. Journal of Hydraulic Engineering，2008，39(7)：836-842.(in Chinese)

[25] Nakai T，Hinokio M. A simple elastoplastic model for normally and overconsolidated soils with unified materical parameters [J]. Soils and Foundations，2004，44(2)：53-70.

[26] 姚仰平，侯偉，羅汀. 土的統一硬化模型[J]. 巖石力學與工程學報，2009，28(10)：2135-2151.

Yao Y P，Hou W，Luo D. Unified hardening model for soils [J]. Chinese Journal of Rock Mechanics and Engineering，2009，28(10)：2135-2151.(in Chinese)

(編輯 胡 玲)

Numerical implementation and applications on subloading surface model with improved overconsolidation state variable

MaShao-kun1，ZhaoNai-feng1，2，PanBai-yu1，ZhouDong1，JiangJie1
(1. College of Civil Engineering and Architecture;Key Laboratory of Disaster Prevention and Structural Safety，Guangxi University，Nanning 530004，P.R.China;2. Department of Geotechnical Engineering，Tongji University，Shanghai 200092，P.R.China)

To describe the complex behavior of overconsolidated soil，an improved overconsolidation state variable subloading surface model was presented. The multiple impact of plastic volumetric and shear strain on overconsolidation state variable R was considered in hardening equation. The stress-update algorithm and the procedure of numerical implementation were introduced. The the model was applied by compiling the interface subroutine in finite element software ABAQUS. The rationality of the subroutine was verified by numerical simulations in different working conditions and load types. At last，Fujinomori Clay was simulated by the model in compression test，and the results were compared with those of UH model and experimentation data. The subroutine was precise and steady andthe model proposed can describe the properties of overconsolidation accurately.

subloading surface model; overconsolidation soils; ABAQUS; subroutine

10.11835/j.issn.1674-4764.2015.02.003

2014-09-15 基金項目：國家自然科學基金資助(51068002；41362016)；廣西巖土力學與工程重點實驗室資助課題(13-KF-02)；上海市青年科技啟明星計劃項目(13QB1404300)

馬少坤(1972-)，男，博士，教授，主要從事地下工程研究工作，(E-mail)mashaokun@sina.com 江 杰(通信作者)，男，博士，高工，(E-mail) jie_jiang001@126.com。

Foundation item：National Natural Science Foundation of China(No.51068002；No.41362016)；The Project Was Funded by Guangxi Key Laboratory of Geomechanics and Geotechnical Engineering(No.13-KF-02)；Sponsored by Shanghai Rising-Star Program(No.13QB1404300)

TU470

A

1674-4764(2015)02-0016-07

Received：2014-09-15

Author brief：Ma Shao kun(1972-), professor, main research interest: underground engineering，(E-mail)mashaokun@sina.com. Jiang Jie(corresponding author),PhD,seniorengineer，(E-mail) jie_jiang001@126.com