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

關(guān)鍵詞：次加載面模型；超固結(jié)土；ABAQUS；用戶子程序

中圖分類號(hào)：TU470 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：1674-4764（2015）02-0016-07

在隧道、基坑等大型地下工程的開挖過程中，開挖面土體因所處位置不同而經(jīng)受不同的加載、卸載等復(fù)雜應(yīng)力路徑，應(yīng)力狀態(tài)發(fā)生變化，由正常固結(jié)狀態(tài)轉(zhuǎn)變?yōu)槌探Y(jié)狀態(tài)。Mesri[1]、Tavenas等[2]及眾多巖土力學(xué)學(xué)者研究表明，自然沉積土多數(shù)為結(jié)構(gòu)性土且處于超固結(jié)狀態(tài)。超固土體具有剪脹、應(yīng)變軟化等特性，其應(yīng)力應(yīng)變關(guān)系與正常固結(jié)土存在較大差異。以劍橋模型[3-4]為代表的經(jīng)典彈塑性理論認(rèn)為：巖土介質(zhì)只存在一個(gè)屈服面，土在卸載再加載過程中，應(yīng)力應(yīng)變關(guān)系為彈性，如圖1（a）所示。而實(shí)際上，正常固結(jié)土一旦卸載就處于超固結(jié)狀態(tài)，土在超固結(jié)狀態(tài)下的再加載過程仍具有塑性變形，如圖1（b）所示[5]。

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

馬少坤，等：改進(jìn)超固結(jié)狀態(tài)參量次加載面模型數(shù)值實(shí)現(xiàn)與應(yīng)用

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

1 次加載面修正劍橋模型

1.1 次加載面概念

次加載面理論認(rèn)為超固結(jié)土存在兩個(gè)應(yīng)力狀態(tài)面：正常固結(jié)屈服面和次加載面，如圖2所示。正常固結(jié)屈服面為修正劍橋模型的屈服面，大小由土在其固結(jié)歷史上所經(jīng)受的最大應(yīng)力水平來確定。次加載面[18-19]為土體卸荷至某一應(yīng)力狀態(tài)時(shí)通過此應(yīng)力狀態(tài)點(diǎn)與正常固結(jié)屈服面幾何相似的面，相似比為R，相似中心為p-q應(yīng)力空間的原點(diǎn)。次加載面位于正常固結(jié)屈服面內(nèi)，隨應(yīng)力狀態(tài)的變化擴(kuò)大或縮小。卸載時(shí)，當(dāng)前應(yīng)力狀態(tài)點(diǎn)遠(yuǎn)離正常固結(jié)屈服面，R減小，次加載面縮??；加載時(shí)，當(dāng)前應(yīng)力狀態(tài)點(diǎn)向正常固結(jié)屈服面靠近，R增大，次加載面擴(kuò)大，若當(dāng)前應(yīng)力狀態(tài)點(diǎn)處于正常固結(jié)屈服面上，則R=1，次加載面與正常固結(jié)面重合，土由超固結(jié)狀態(tài)轉(zhuǎn)變?yōu)檎９探Y(jié)狀態(tài)。

1.2 屈服面方程與塑性勢函數(shù)

引入描述超固結(jié)程度的超固結(jié)比（OCR）概念，給出本文模型兩個(gè)應(yīng)力狀態(tài)面相似比（R）的定義式：

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

為準(zhǔn)確描述土及軟巖的力學(xué)特征，近代土力學(xué)尤其是次加載面理論建議采用相關(guān)聯(lián)流動(dòng)法則[3]，因此，式（2）即為本文模型的塑性勢函數(shù)表達(dá)式。

確定了塑性勢函數(shù)后，塑性應(yīng)變率即可由以下流動(dòng)法則給出：

1.3 協(xié)調(diào)方程

由式（2）等號(hào)兩邊取微分，得到滿足塑性一致條件的協(xié)調(diào)方程式（4）。

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

dR=-mRχMlnRdεpv（6）

式中：mR為材料參數(shù)，其大小表示超固結(jié)狀態(tài)隨塑性體應(yīng)變?cè)黾佣У目炻?/p>

式（6）中dR僅與塑性體應(yīng)變相關(guān)，沒有考慮到塑性剪應(yīng)變對(duì)超固結(jié)性的影響?？紤]到超固結(jié)性的發(fā)展受塑性體應(yīng)變和塑性剪應(yīng)變綜合影響的情況，本文綜合Asaoka等[15-16]、Yamakawa等[17]的建議式，把塑性體應(yīng)變?cè)隽颗c塑性剪應(yīng)變?cè)隽繉?duì)超固結(jié)性的影響做加權(quán)，得出本文模型dR的表達(dá)式

式中： η為非負(fù)的材料參數(shù)，表示塑性剪應(yīng)變?cè)隽繉?duì)超固結(jié)性的發(fā)展貢獻(xiàn)比。

2 模型應(yīng)力更新算法

本文采用計(jì)算精度和穩(wěn)定性都較高的隱式積分算法[20]，基本思路如圖3示σN點(diǎn)為已知的第N增量步的應(yīng)力狀態(tài)，把給定的總應(yīng)變?cè)隽喀う舙N+1全部視為彈性應(yīng)變?cè)隽?，?duì)應(yīng)的應(yīng)力增量為（dσ）tr，得試探應(yīng)力狀態(tài)σtrN+1。進(jìn)行屈服判斷：若f（σtrN+1）≤0，說明材料沒有屈服，上述試探應(yīng)力狀態(tài)即為實(shí)際應(yīng)力狀態(tài)；若f（σtrN+1）>0說明材料屈服，需進(jìn)行塑形修正。塑性修正過程中，根據(jù)最近點(diǎn)回映算法確定塑性應(yīng)變?cè)隽浚éう舙N+1）1，通過Hooke定律計(jì)算經(jīng)初次迭代后的應(yīng)力狀態(tài)σ1N+1，再進(jìn)行屈服判斷。重復(fù)以上步驟直至f（σkN+1）=0，此時(shí)的σkN+1即為第N+1增量步實(shí)際的應(yīng)力狀態(tài)，迭代結(jié)束。

2.1 初始變量計(jì)算

由狀態(tài)變量σN、RN及材料參數(shù)M可得：

面和正常固結(jié)屈服面的大小。

2.2 彈性試算

試探應(yīng)力分量σtrN+1由Hooke定律得出

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

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

2.3 初始屈服判斷

首先由試探應(yīng)力分量σtr計(jì)算試探平均應(yīng)力ptr、試探剪切應(yīng)力qtr，進(jìn)而計(jì)算初始屈服函數(shù)ftr：

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

若屈服函數(shù)f小于某個(gè)容許誤差值ftol，則應(yīng)力狀態(tài)處于彈性階段，進(jìn)行下文2.6節(jié)的處理，否則進(jìn)行塑性修正。本文設(shè)定ftol=1×10-5。

2.4 塑性修正

（1）更新第k步迭代時(shí)φ的值φk：

（4）R的更新。由于迭代步間的應(yīng)變?cè)隽亢蛻?yīng)力增量一般較小，故本文假定在迭代步間的R保持不變，為上一迭代步末的值。在完成pc、p、q的更新后，再利R

v2，作為下一迭代步的初始值進(jìn)行后續(xù)迭代。實(shí)際計(jì)算表明，以上方法不僅使結(jié)果準(zhǔn)確可靠，也增大了迭代效率。

2.5 一致切線模量

在完成每個(gè)積分點(diǎn)應(yīng)力應(yīng)變狀態(tài)的計(jì)算后，提供與材料彈塑性本構(gòu)關(guān)系及其積分算法一致的一致切線模量（Jacobian矩陣）[21]。

2.6 變量更新與存儲(chǔ)

由下式進(jìn)行應(yīng)變、孔隙比等的更新，并進(jìn)行狀態(tài)變量 STATEV（NSTATV）的存儲(chǔ)：

3 程序驗(yàn)證

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

式（7）中的材料參數(shù)mR在文獻(xiàn)[5，15-16]、文獻(xiàn)[23]中的取值在[2， 10]之間，經(jīng)驗(yàn)證，當(dāng)mR取8.0時(shí)，本文模型的模擬結(jié)果更合理；對(duì)式（7′）中與超固結(jié)狀態(tài)參數(shù)R的發(fā)展有關(guān)的參數(shù)η，根據(jù)本文模擬，取0.8較合適。

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

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

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

圖6 孔隙比與平均有效應(yīng)力關(guān)系

Fig.6 Relationship between porosity ratio and mean effective stress

圖7為本文算法與修正劍橋模型結(jié)果的偏應(yīng)力軸向應(yīng)變關(guān)系圖。如圖7（a）、（b）所示，修正劍橋模型（UMAT）模擬超固結(jié)土的應(yīng)力應(yīng)變關(guān)系曲線有突變點(diǎn)，而本文算法（TEST）模擬的應(yīng)力應(yīng)變關(guān)系曲線光滑連續(xù)，說明本文算法能夠更準(zhǔn)確地描述土實(shí)際固結(jié)過程中連續(xù)平滑的彈塑性應(yīng)力應(yīng)變關(guān)系。

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

4 本文模型應(yīng)用

應(yīng)用本文模型，對(duì)Fujinomori 黏土進(jìn)行三軸壓縮試驗(yàn)?zāi)M，并與該土樣的室內(nèi)試驗(yàn)結(jié)果[25]、UH模型模擬結(jié)果[26]進(jìn)行比較。土樣參數(shù)如表2所示，具體試驗(yàn)見文獻(xiàn)[25]

圖8為數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果對(duì)比圖。如圖所示，本文模型和UH模型[26]均能準(zhǔn)確刻畫超固結(jié)土應(yīng)力應(yīng)變關(guān)系的非線性、應(yīng)變軟化等一般特性。然而在反映超固結(jié)土的峰值強(qiáng)度、殘余強(qiáng)度以及剪切應(yīng)力比與軸向應(yīng)變關(guān)系的發(fā)展趨勢方面，本文模型顯然更具優(yōu)勢，與試驗(yàn)結(jié)果更吻合。由此看來，本文模型綜合考慮塑性剪應(yīng)變與塑性體應(yīng)變對(duì)超固結(jié)狀態(tài)的影響，能夠更準(zhǔn)確地反映超固結(jié)土實(shí)際狀態(tài)的應(yīng)力應(yīng)變關(guān)系。

5 結(jié)論

建立了改進(jìn)超固結(jié)狀態(tài)參量的次加載面模型，通過編制該模型對(duì)應(yīng)的用戶子程序，實(shí)現(xiàn)了對(duì)有限元軟件的二次開發(fā)，建立了適用于超固結(jié)土的數(shù)值模擬平臺(tái)。隨后，應(yīng)用模型分析了不同工況和加載條件下土的力學(xué)特性，并與UH模型的數(shù)值模擬結(jié)果及試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析，得如下結(jié)論：

1）本文模型在模擬正常固結(jié)土的三軸壓縮試驗(yàn)時(shí)，能準(zhǔn)確描述變形特征、孔隙比與有效應(yīng)力關(guān)系、應(yīng)力應(yīng)變關(guān)系等特性。

2）在描述超固結(jié)土的力學(xué)特性方面，本文模型相對(duì)修正劍橋模型能夠更為連續(xù)平滑地模擬實(shí)際超固結(jié)土的彈塑性應(yīng)力應(yīng)變關(guān)系。

3）在描述超固結(jié)土的力學(xué)行為方面，本文模型相對(duì)UH模型能夠更準(zhǔn)確地刻畫超固結(jié)土的應(yīng)變軟化、峰值強(qiáng)度與殘余強(qiáng)度值以及應(yīng)力應(yīng)變關(guān)系等特征。

參考文獻(xiàn)：

[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]張鋒，葉冠林. 計(jì)算土力學(xué)[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]孔亮，鄭穎人，姚仰平. 基于廣義塑性力學(xué)的土體次加載面循環(huán)塑性模型（Ⅰ）：理論與模型[J]. 巖土力學(xué)，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]黃雨，周子舟. 下負(fù)荷面劍橋模型在ABAQUS中的開發(fā)實(shí)現(xiàn)[J].巖土工程學(xué)報(bào)，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]范慶來，欒茂田，倪宏革. 循環(huán)荷載作用下軟基上大圓筒結(jié)構(gòu)彈塑有效應(yīng)力分析[J]. 水利學(xué)報(bào)，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]姚仰平，侯偉，羅汀. 土的統(tǒng)一硬化模型[J]. 巖石力學(xué)與工程學(xué)報(bào)，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）

（編輯 胡 玲）