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

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

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

中圖分類號：TU470 文獻標志碼：A 文章編號：1674-4764（2015）02-0016-07

在隧道、基坑等大型地下工程的開挖過程中，開挖面土體因所處位置不同而經受不同的加載、卸載等復雜應力路徑，應力狀態發生變化，由正常固結狀態轉變為超固結狀態。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*的定義模糊，限制了模型的推廣與應用。

馬少坤，等：改進超固結狀態參量次加載面模型數值實現與應用

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

1 次加載面修正劍橋模型

1.1 次加載面概念

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

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

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

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

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

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

1.3 協調方程

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

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

dR=-mRχMlnRdεpv（6）

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

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

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

2 模型應力更新算法

本文采用計算精度和穩定性都較高的隱式積分算法[20]，基本思路如圖3示σN點為已知的第N增量步的應力狀態，把給定的總應變增量ΔεpN+1全部視為彈性應變增量，對應的應力增量為（dσ）tr，得試探應力狀態σtrN+1。進行屈服判斷：若f（σtrN+1）≤0，說明材料沒有屈服，上述試探應力狀態即為實際應力狀態；若f（σtrN+1）>0說明材料屈服，需進行塑形修正。塑性修正過程中，根據最近點回映算法確定塑性應變增量（ΔεpN+1）1，通過Hooke定律計算經初次迭代后的應力狀態σ1N+1，再進行屈服判斷。重復以上步驟直至f（σkN+1）=0，此時的σkN+1即為第N+1增量步實際的應力狀態，迭代結束。

2.1 初始變量計算

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

面和正常固結屈服面的大小。

2.2 彈性試算

試探應力分量σtrN+1由Hooke定律得出

為彈性體積模量，μ為剪切模量。

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

2.3 初始屈服判斷

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

ftr=qtrN+12M2+ptrN+1ptrN+1-pc，N（11）

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

2.4 塑性修正

（1）更新第k步迭代時φ的值φk：

（4）R的更新。由于迭代步間的應變增量和應力增量一般較小，故本文假定在迭代步間的R保持不變，為上一迭代步末的值。在完成pc、p、q的更新后，再利R

v2，作為下一迭代步的初始值進行后續迭代。實際計算表明，以上方法不僅使結果準確可靠，也增大了迭代效率。

2.5 一致切線模量

在完成每個積分點應力應變狀態的計算后，提供與材料彈塑性本構關系及其積分算法一致的一致切線模量（Jacobian矩陣）[21]。

2.6 變量更新與存儲

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

3 程序驗證

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

式（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]算法，能更準確地描述實際正常固結土和超固結土固結過程中的應力狀態變化。

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

圖6 孔隙比與平均有效應力關系

Fig.6 Relationship between porosity ratio and mean effective stress

圖7為本文算法與修正劍橋模型結果的偏應力軸向應變關系圖。如圖7（a）、（b）所示，修正劍橋模型（UMAT）模擬超固結土的應力應變關系曲線有突變點，而本文算法（TEST）模擬的應力應變關系曲線光滑連續，說明本文算法能夠更準確地描述土實際固結過程中連續平滑的彈塑性應力應變關系。

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

4 本文模型應用

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

圖8為數值模擬結果與試驗結果對比圖。如圖所示，本文模型和UH模型[26]均能準確刻畫超固結土應力應變關系的非線性、應變軟化等一般特性。然而在反映超固結土的峰值強度、殘余強度以及剪切應力比與軸向應變關系的發展趨勢方面，本文模型顯然更具優勢，與試驗結果更吻合。由此看來，本文模型綜合考慮塑性剪應變與塑性體應變對超固結狀態的影響，能夠更準確地反映超固結土實際狀態的應力應變關系。

5 結論

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

1）本文模型在模擬正常固結土的三軸壓縮試驗時，能準確描述變形特征、孔隙比與有效應力關系、應力應變關系等特性。

2）在描述超固結土的力學特性方面，本文模型相對修正劍橋模型能夠更為連續平滑地模擬實際超固結土的彈塑性應力應變關系。

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）

（編輯 胡 玲）