楊云浩，陳鴻杰，王 偉

（1.中國水電顧問集團成都勘測設(shè)計研究院，成都 610072；2.河海大學(xué)a.巖土力學(xué)與堤壩工程教育部重點實驗室；b.巖土工程研究所，南京 210098）

彈塑性各向異性損傷模型的FLAC3D開發(fā)與數(shù)值驗證

楊云浩1，陳鴻杰2a，2b，王 偉2a，2b

（1.中國水電顧問集團成都勘測設(shè)計研究院，成都 610072；2.河海大學(xué)a.巖土力學(xué)與堤壩工程教育部重點實驗室；b.巖土工程研究所，南京 210098）

探討適用于硬巖的各向異性損傷模型在FLAC3D中的實現(xiàn)方法。首先，基于不可逆熱力學(xué)框架下的彈塑性損傷本構(gòu)基本理論，結(jié)合彈塑性力學(xué)塑性勢與非關(guān)聯(lián)塑性流動法則，導(dǎo)出了彈塑性各向異性損傷本構(gòu)模型的增量形式及損傷演化速率的計算公式，在此基礎(chǔ)上，采用復(fù)合M C準(zhǔn)則并考慮黏聚力隨損傷的逐步弱化，詳細給出了該損傷模型在FLAC3D中的實現(xiàn)流程并探討了相關(guān)細節(jié)處理方法，編寫了相應(yīng)的DLL文件以供調(diào)用。進而，通過圓柱試樣的壓縮數(shù)值試驗驗證其正確性與合理性，結(jié)果表明該損傷模型能夠較好地描述巖石材料的非線性力學(xué)行為；以一個概化的地下洞室分層開挖實例初步檢驗了模型的應(yīng)用效果，計算所得洞周損傷較大部位與實際工程中直接或間接觀測到的圍巖損傷情況較為一致，由此表明以拉應(yīng)變度量的損傷及其演化方程能夠合理地反映硬巖的損傷本質(zhì)。

不可逆熱力學(xué)；彈塑性；各向異性損傷；FLAC3D；本構(gòu)模型

1 研究背景

在巖體地下工程領(lǐng)域，施工期圍巖的穩(wěn)定性長期以來都是采用常規(guī)彈塑性本構(gòu)模型進行分析。對于硬巖地下洞室，開挖卸荷引起的原巖應(yīng)力重分布導(dǎo)致圍巖內(nèi)微裂紋大量萌生，微裂紋的逐步擴展使得圍巖的剛度和強度屬性發(fā)生劣化［1－2］，從而在宏觀的應(yīng)力應(yīng)變曲線上表現(xiàn)出非線性特征，而傳統(tǒng)的彈塑性模型通常不能反映這些特征，相應(yīng)的穩(wěn)定性分析成果的可靠性也就大打折扣。

近年來，不少學(xué)者引入損傷的概念來描述巖石受載過程中變形和強度屬性的劣化并建立了相應(yīng)的損傷本構(gòu)模型。文獻［3］研究了彈脆塑性損傷本構(gòu)模型，但未明確損傷的啟動與演化準(zhǔn)則。文獻［4－6］所提出的彈塑性損傷本構(gòu)模型將損傷變量取為標(biāo)量，以應(yīng)變能或等效塑性應(yīng)變作為損傷演化參量。然而，文獻［7－8］指出：張拉性裂紋的出現(xiàn)是低圍壓條件下硬巖損傷的主要形式，其本質(zhì)是拉裂紋萌生、擴展直至裂紋相互貫通的結(jié)果。因而，以應(yīng)變能或等效塑性應(yīng)變作為損傷演化參量比較適合于黏土巖和泥巖等軟巖，卻不能很好地反映硬巖損傷的本質(zhì)。文獻［1］提出的損傷本構(gòu)模型將體積應(yīng)變作為損傷參量，體現(xiàn)了硬巖損傷細觀上的裂紋萌生擴展導(dǎo)致宏觀上出現(xiàn)體脹這一典型特征，但該模型只能考慮各向同性損傷。文獻［9］通過大理巖三軸卸荷試驗發(fā)現(xiàn)試樣主破裂面周邊常有較多的近垂直于卸荷方向的微張裂隙，由此表明試樣軸向與徑向的損傷必然是不同的。因此，對于硬巖，有必要考慮損傷的各向異性。文獻［10－11］分別針對花崗巖和混凝土所建立的損傷模型考慮各向異性損傷，將二階損傷張量的損傷主值表達為拉應(yīng)變的冪函數(shù)，但該模型僅考慮彈性損傷且其損傷演化方程在理論上不夠嚴(yán)謹(jǐn)。文獻［12］所提出的彈塑性各向異性損傷模型將二階損傷張量與應(yīng)變張量的正錐直接關(guān)聯(lián)并假定二者共軸。該模型建立在不可逆熱力學(xué)基礎(chǔ)上，理論依據(jù)嚴(yán)謹(jǐn)，所采用的損傷準(zhǔn)則及損傷演化方程也能很好地反映硬巖損傷的物理機制，模型參數(shù)少且物理意義較明確，有良好的應(yīng)用前景。

借助大型商業(yè)軟件提供的二次開發(fā)接口，編制可供其調(diào)用的用戶自定義本構(gòu)模塊，是將最新理論研究成果快速應(yīng)用于工程實踐的有效途徑。目前，在硬巖各向異性損傷模型的FLAC3D二次開發(fā)方面所開展的工作還較少。鑒于此，本文在文獻［12］的理論公式基礎(chǔ)上，首先導(dǎo)出了彈塑性各向異性損傷本構(gòu)模型的增量形式及二階損傷張量增量的具體表達式，進而，采用復(fù)合M C準(zhǔn)則并考慮黏聚力隨損傷的逐步弱化，詳細給出了該本構(gòu)模型在FLAC3D下的實現(xiàn)流程并探討了相關(guān)細節(jié)處理方法。

2 彈塑性各向異性損傷本構(gòu)模型理論

2.1 基于Helmholtz自由能勢的彈塑性各向異性損傷本構(gòu)模型

考慮損傷的各向異性，Dragon等人采用了如下形式的Helmholtz自由能勢［12］，即

式中：λ，μ為拉梅彈性常數(shù)；εe為彈性應(yīng)變張量；D為二階損傷張量。符號“tr”表示對張量求跡；“：”表示張量的雙點乘；等式右邊后2項為損傷引起的對自由能的修正項，a1，a2為待定參數(shù)，用于描述損傷所造成的彈性屬性的劣化，其單位與λ，μ相同。

求ψe（εe，D）關(guān)于εe的偏導(dǎo)，可得式（2）所示彈塑性損傷本構(gòu)關(guān)系［12］

式中：C（D）為彈性損傷剛度陣，是一個四階張量，其各分量如式（4）所示；ε，εp分別為總應(yīng)變和塑性應(yīng)變張量；δij為Kronecker記號。

2.2 本構(gòu)模型的增量形式

為便于數(shù)值實現(xiàn)，對式（3）的第一式兩邊關(guān)于虛擬時標(biāo)求導(dǎo)，得到增量形式本構(gòu)方程為

通過引入塑性屈服準(zhǔn)則函數(shù)fp（σ，γp）和塑性勢函數(shù)g，可將式（5）表達為更為具體的形式，簡要推導(dǎo)如下。

首先，由塑性一致性條件：f·p（σ，γp）＝0，有

式中γp為廣義塑性剪應(yīng)變，作為塑性硬（軟）化參量，計算公式如下：

式中λp為塑性乘子，由式（9）可得

由式（10）可見，如果未發(fā)生損傷的演化（此時該式分子部分的第二項為零），則λp退化為常規(guī)彈塑性下的形式。

式中符號“”表示并矢運算。

與常規(guī)彈塑性本構(gòu)關(guān)系式對比可見，式（11）等號右邊第一部分即為常規(guī)的彈塑性應(yīng)力修正量，后兩部分則為損傷相關(guān)應(yīng)力修正量。同樣，若不存在損傷演化，該式退變?yōu)槌Ｒ?guī)的彈塑性本構(gòu)方程式。

2.3 損傷準(zhǔn)則函數(shù)與損傷演化方程

由上述可見，彈塑性各向異性損傷本構(gòu)模型的關(guān)鍵是損傷的度量，也即損傷準(zhǔn)則及損傷演化方程的建立。在熱力學(xué)框架內(nèi)，損傷演化法則是由和損傷張量相關(guān)的共軛損傷力空間內(nèi)的損傷耗散勢所決定的［13］。文獻［13］給出了損傷準(zhǔn)則的一般形式，文獻［14］進一步給出了一個較為具體的形式，但因過于復(fù)雜，數(shù)值實現(xiàn)較為困難。文獻［12］基于室內(nèi)試驗觀測，將損傷與拉應(yīng)變直接關(guān)聯(lián)，建立了如下?lián)p傷準(zhǔn)則函數(shù)，即

進一步假定損傷演化率張量D·與ε＋共軸，則可得如下?lián)p傷演化方程

損傷乘子λd的具體形式由損傷一致性條件求得，即

式（12）和式（13）所示損傷準(zhǔn)則與損傷演化方程形式簡單、易于計算，而且物理意義也較為明確，在一定程度上反映了硬巖受載過程中損傷發(fā)生與發(fā)展的內(nèi)在物理機制。鑒于此，本文采用該損傷準(zhǔn)則與損傷演化方程來描述硬巖受載過程中的各向異性損傷。

3 彈塑性各向異性損傷本構(gòu)的FLAC3D二次開發(fā)

3.1 屈服準(zhǔn)則的選取

考慮工程上的實用性，采用復(fù)合M C準(zhǔn)則及相應(yīng)的非關(guān)聯(lián)塑性流動法則（如式（15）所示），同時不考慮損傷對塑性屈服的影響。

式中：fs，ft分別為剪切、拉伸屈服函數(shù)；gs，gt分別為非關(guān)聯(lián)的剪切、拉伸屈服塑性勢函數(shù)；，ψ，σt分別為內(nèi)摩擦角、剪脹角和單軸抗拉強度；cini，cres分別為黏聚力的初始值和殘余值。需注意的是，為了考慮黏聚力隨損傷發(fā)展的逐步弱化，此處借鑒了文［1］的做法：將黏聚力c看作是廣義塑性剪應(yīng)變γp的函數(shù)，即，隨γp的增長，c呈指數(shù)函數(shù)形式從初值降至殘余值，參數(shù)b控制c的下降速度。

約定σ1≤σ2≤σ3（即拉為正、壓為負(fù)），則在（σ1，σ3）平面上，復(fù)合M C準(zhǔn)則屈服面曲線如圖1［15］所示。圖中：h＝0為拉、剪屈服區(qū)分界線；剪切屈服區(qū)定義為fs＜0且h＜0；拉伸屈服區(qū)定義為ft＞0且h＞0；彈性區(qū)為fs＞0且ft＜0。

圖1 復(fù)合M C準(zhǔn)則屈服面及屈服區(qū)定義［15］Fig.1 Composite M ohr Coulomb failure criterion and domains used in the definition of the flow rule［15］

3.2 彈塑性各向異性損傷本構(gòu)的數(shù)值實現(xiàn)流程

采用前述理論公式，借鑒彈塑性模型的增量式迭代計算格式，建立了如圖2所示的彈塑性各向異性損傷模型的數(shù)值實現(xiàn)流程。

3.3 模型關(guān)鍵變量的計算方法

由圖2可見，在彈塑性各向異性損傷本構(gòu)數(shù)值實現(xiàn)時，在每個計算時步都要更新總應(yīng)變和彈性應(yīng)變，塑性應(yīng)變在單元進入塑性屈服之后要更新，而損傷張量在有損傷的進一步演化時需進行更新。借鑒FLAC3D內(nèi)置應(yīng)變軟化模型中塑性硬化參量的計算方法，對總應(yīng)變、彈性應(yīng)變、塑性應(yīng)變及損傷的計算采用逐時步累加的方法求取。實際上，F(xiàn)LAC3D處理每個單元時是對其內(nèi)各個子單元逐一處理的，即主計算程序傳給本構(gòu)程序的應(yīng)變增量是子單元的應(yīng)變增量，本構(gòu)程序輸出的也是子單元的應(yīng)力張量。為求得計算時步內(nèi)各關(guān)鍵變量的增量，本文采取了體積平均的處理方法。

以時步內(nèi)總應(yīng)變增量為例，其計算公式為

式中：VZone為單元的體積；VSubZone（k）為單元內(nèi)第k個子單元的體積；ΔεSijubZone（k）為第k個子單元的應(yīng)變增量；Overlayes為單元的體積被使用的次數(shù)，取值為1或2；TotSubZones為單元內(nèi)包含的子單元總數(shù)。

圖2 彈塑性各向異性損傷本構(gòu)在FLAC3D內(nèi)的計算循環(huán)流程Fig.2 Com putation procedures of elastoplastic anisotropic damage constitutivemodel in FLAC3D

圖3 數(shù)值試驗獲得的不同圍壓下的應(yīng)力－應(yīng)變曲線Fig.3 Stress strain curves under different confining pressures obtained from numerical compression tests

圖4 不同圍壓下數(shù)值試樣內(nèi)某單元損傷演化過程Fig.4 Damage evolution under different confining pressures obtained from numerical compression tests

表1 彈塑性各向異性損傷本構(gòu)模型參數(shù)取值Table 1 Parameters used in elastoplastic anisotropic damage constitutivemodel

4 彈塑性各向異性損傷本構(gòu)模型的數(shù)值驗證

按前述數(shù)值實現(xiàn)流程以及相關(guān)考慮，編制了彈塑性各向異性損傷本構(gòu)的動態(tài)鏈接庫（DLL）文件以供FLAC3D調(diào)用。為檢驗數(shù)值實現(xiàn)的正確性，進行圓柱試樣壓縮的數(shù)值模擬實驗。試樣上下兩端施加方向相反的速度邊界條件模擬軸向加載，試樣徑向施加均勻應(yīng)力邊界條件以模擬圍壓施加。數(shù)值試驗采用的損傷相關(guān)力學(xué)參數(shù)來自文獻［12］（如表1所示），得到的不同圍壓下的應(yīng)力應(yīng)變曲線如圖3所示。

由應(yīng)力應(yīng)變曲線可見，在應(yīng)力達到峰值前有明顯的非線性段，這是由于損傷導(dǎo)致材料剛度逐漸弱化的結(jié)果，這一點與真實應(yīng)力應(yīng)變曲線峰前段特征比較吻合。同時也可看到，由于采用的復(fù)合M C屈服準(zhǔn)則考慮了黏聚力的劣化，所以模型亦能在一定程度上模擬巖石類材料的峰后應(yīng)變軟化力學(xué)行為。

圖4為不同圍壓條件下試樣內(nèi)某單元的損傷主值隨加載過程的變化（D2與D1的量值及演化過程基本相同，此處略），可見：在彈性損傷階段，損傷值增加平緩；塑性屈服后，損傷值呈加速增長趨勢。

分析可知，圓柱試樣壓縮試驗條件下，按FLAC3D的約定，應(yīng)變以拉為正、壓為負(fù)，則有ε1＝ε2＞0＞ε3，于是有ε＋1＝ε＋2，ε＋3＝0（ε＋的含義見2.3節(jié)所述），由于假定損傷張量D與拉應(yīng)變張量ε＋共軸，因此，理論上應(yīng)有D1＝D2，D3＝0。對比理論分析與數(shù)值模擬的結(jié)果，可見二者基本吻合，這與D3不嚴(yán)格為零的原因與數(shù)值計算的精度有關(guān)。進一步考察損傷的分布（如圖5所示），可見：較大的損傷出現(xiàn)在與試樣軸向呈一定夾角的“X”形的窄條帶上，與真實巖樣三軸壓縮試驗的最終破裂面形狀也較為吻合。

綜上所述，本文彈塑性各向異性損傷本構(gòu)模型的數(shù)值實現(xiàn)是正確的。

圖5 試樣內(nèi)損傷分布云圖Fig.5 Contours of damage obtained from numerical triaxial com pression test

5 彈塑性各向異性損傷模型初步應(yīng)用效果

建立一高跨比為2.5的城門洞型地下洞室分層開挖概化模型（為節(jié)省篇幅，模型圖略去）。施加側(cè)壓系數(shù)為3.0的初始地應(yīng)力，大主應(yīng)力方向垂直于洞室軸向，最大值為15.0 MPa。洞室分2層開挖，應(yīng)用彈塑性各向異性損傷模型模擬洞室開挖后的圍巖力學(xué)響應(yīng)，圍巖力學(xué)參數(shù)取表1所給值。

分層開挖過程中，洞周圍巖損傷演化情況如圖6所示（因損傷主值D2較小、D3近似為零，此處略去其圖）。本例地應(yīng)力場側(cè)壓系數(shù)及洞室高跨比均較大，因而邊墻部位拉應(yīng)變大于壓應(yīng)變，而頂拱部位則剛好相反，故而邊墻部位的損傷遠大于頂拱部位，且損傷程度隨開挖的進行而進一步增大。

圖6 分層開挖過程中的圍巖損傷演化Fig.6 Evolution of damages of surrounding rock during staged excavation

此外，損傷的發(fā)生與進一步演化必將弱化圍巖的剛度，從而導(dǎo)致圍巖變形較不考慮損傷影響時有所增加。為進行對比，在保持相同的變形與強度參數(shù)值的前提下，將損傷的啟動門檻值參數(shù)調(diào)高，使損傷始終不發(fā)生，進行不考慮損傷影響的開挖模擬計算。

圖7所示為特征點處的水平位移對比，可見：由于拱腳部位為壓應(yīng)力集中區(qū)，拉應(yīng)變量值很小，損傷也小，因此考慮損傷與不考慮損傷情況下的位移基本相當(dāng)；邊墻部位為卸荷松弛區(qū)，拉應(yīng)變量值大，損傷相應(yīng)也大，故而考慮損傷情況下的水平位移較不考慮損傷的情況有顯著增加。

圖7 關(guān)鍵點水平位移對比Fig.7 Horizontal disp lacements of four key points under damage and non damage conditions

6 結(jié) 語

（1）對于開挖于完整性較好的硬巖巖體內(nèi)的地下洞室，圍巖應(yīng)力與變形仿真分析結(jié)果的可靠性在很大程度上取決于所采用的本構(gòu)模型。常規(guī)的彈塑性模型因其難以很好地描述硬巖受載過程表現(xiàn)出的非線性力學(xué)行為，故而有必要引入彈（塑）性損傷本構(gòu)模型開展仿真分析。通過FLAC3D的二次開發(fā)接口編寫自定義損傷模型是開展考慮損傷效應(yīng)的圍巖穩(wěn)定分析的一條捷徑。

（2）彈塑性各向異性損傷本構(gòu)的應(yīng)力增量計算公式與常規(guī)彈塑性的計算式相比，增加了2個損傷修正項；在無損傷演化發(fā)生時，則在形式上可退化為常規(guī)彈塑性本構(gòu)模型。

（3）由于考慮了損傷的各向異性，需要在每個時步內(nèi)計算單元的應(yīng)變張量主值及其主向以便更新拉應(yīng)變張量，因而應(yīng)用該模型進行計算耗時較多，今后還需在模型的執(zhí)行效率方面進行改進，使其可以有效地應(yīng)用到大規(guī)模問題的仿真分析中。

［1］ SALARIM R，SAEB S，WILLAM K J，etal.A Coupled Elastoplastic Damage Model for Geomaterials［J］.Com puter Methods in Applied Mechanics and Engineering，2004，193（27／29）：2625－2643.

［2］ GOLSHANIA，YOSHIAKIO，MASANOBU O，et al.AMicromechanical Model for Brittle Failure of Rock and Its Relation to Crack Growth Observed in Triaxial Compres sion Tests of Granite［J］.Mechanics of Materials，2006，38（4）：287－303.

［3］ 劉洪永，程遠平，趙長春，等.采動煤巖體彈脆塑性損傷本構(gòu)模型及應(yīng)用［J］.巖石力學(xué)與工程學(xué)報，2010，

29（2）：358－365.（LIU Hong yong，CHENG Yuan ping，ZHAO Chang chun，et al.Constitutive Model for Elasto Brittle Plastic Damage of Coal Rock Mass Due to Mining and Its Application［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（2）：358－365.（in Chinese））

［4］ SHAO JF，JIA Y，KONDO D，et al.A Coupled Elasto plastic Damage Model for Semi brittle Materials and Ex tension to Unsaturated Conditions［J］.Mechanics of Ma terials，2006，38（3）：218－232.

［5］ 賈善坡，陳衛(wèi)忠，于洪丹，等.泥巖隧道施工過程中滲流場與應(yīng)力場全耦合損傷模型研究［J］.巖土力學(xué)，2009，30（1）：19－26.（JIA Shan po，CHENWei zhong，YU Hong dan，et al.Research on Seepage Stress Cou pling Damage Model of Boom Clay during Tunneling［J］.Rock and Soil Mechanics，2009，30（1）：19－26.（in Chinese））

［6］ 賈善坡，陳衛(wèi)忠，于洪丹，等.泥巖彈塑性損傷本構(gòu)模型及其參數(shù)辨識［J］.巖土力學(xué)，2009，30（12）：3607－3614.（JIA Shan po，CHEN Wei zhong，YU Hong dan，et al.Parameter Identification of New Elastoplastic Dam age Constitutive Model for Claystone［J］.Rock and Soil Mechanics，2009，30（12）：3607－3614.（in Chinese））［7］ MITAIM S，DETOURNAY E.Damage around a Cylindri cal Opening in a Brittle Rock Mass［J］.International Journal of Rock Mechanics and Mining Sciences，2004，41（8）：1447－1457.

［8］ CAIM，KAISER P K，TASAKA Y.Generalized Crack Initiation and Crack Damage Stress Thresholds of Brittle Rock Masses near Underground Excavations［J］.Interna tional Journal of Rock Mechanics and Mining Sciences，2004，41（5）：833－847.

［9］ 黃潤秋，黃 達.高地應(yīng)力條件下卸荷速率對錦屏大理巖力學(xué)特性影響規(guī)律試驗研究［J］.巖石力學(xué)與工程學(xué)報，2010，29（1）：21－33.（HUANG Run qiu，HUANG Da.Experimental Research on Affection Laws of Unloading Rates on Mechanical Properties of Jinping Mar ble under High Geostress［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（1）：21－33.（in Chinese））

［10］周維垣，剡公瑞，楊若瓊.巖體彈脆性損傷本構(gòu)模型及工程應(yīng)用［J］.巖土工程學(xué)報，1998，20（5）：54－57.

（ZHOU Wei yuan，YAN Gong rui，YANG Ruo qiong.Elasto brittle Damage Model for Rockmass Based on Field Tests in Laxiwa Arch Dam Site［J］.Chinese Journal of Geotechnical Engineering，1998，20（5）：54－57.（in Chinese））

［11］高政國，黃達海，趙國藩.碾壓混凝土的正交異性損傷本構(gòu)模型研究［J］.水利學(xué)報，2001，（5）：58－64.（GAO Zheng guo，HUANG Da hai，ZHAO Guo fan.An Orthotropic Damage Constitutive Model for RCC［J］.Journal of Hydraulic Engineering，2001，（5）：58－64.（in Chinese））

［12］CHIARELLIA S，SHAO J F，HOTEIT N.Modeling of Elastoplastic Damage Behavior of a Claystone［J］.Inter national Journal of Plasticity，2003，19（1）：23－45.

［13］JU JW.On Energy based Coupled Elastoplastic Damage Theories：Constitutive Modeling and Computational As pects［J］.International Journal of Solids and Structures，1989，25（7）：803－833.

［14］MURAKAMI S，KAMIYA K.Constitutive and Damage Evolution Equations of Elastic brittle Materials Based on Irreversible Thermodynamics［J］.International Journal of Mechanical Sciences，1997，39（4）：473－486.

［15］Itasca Consulting Group Inc.FLAC3DUser’s Manual：Theory and Background in FLAC3D［M］.Minneapolis：Itasca Consulting Group，2005.

（編輯：姜小蘭）

Implementation and Numerical Verification of Elastoplastic Anisotropic Damage Constitutive M odel in FLAC3D

YANG Yun hao1，CHEN Hong jie2，3，WANGWei2，3
（1.Hydro China Chengdu Engineering Corporation，Chengdu 610072，China；2.Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering，Hobai University，Nanjing 210098，China；3.Research Institute of Geotechnical Engineering，Hohai University，Nanjing 210098，China）

The implementation of an anisotropic damage constitutivemodel in FLAC3Dfor hard rock is presented in this paper.An incremental form of elastoplastic anisotropic damage constitutive equation is first deduced using the plastic mechanics theory and the elastoplastic damage constitutive theory within the framework of thermodynamics of irreversible process.Meanwhile the concrete form of evolution rate of damage tensor is determined by the consis tence condition of damage criterion.On the basis of these formulas，the detailed flowchart of numerical calculation by using the elastoplastic anisotropic damagemodel and M C yield criterion is established and some details about

the numerical calculation are also discussed.A DLL（Dynamic Link Library）file corresponding to the elastoplastic anisotropic damagemodel iswritten to be called by the FLAC3Dsoftware.The uniaxial and triaxial compression nu merical tests of cylinder sample are then conducted to test the validity of this user defined model.The results show that thismodel is capable of describing the nonlinear behavior of hard rock.The validity of elastoplastic anisotropic damagemodel in practice is further tested by applying it to a numerical simulation of staged excavation of a fictitious underground opening.Result show that the highly damaged zone around the opening determined by the simulation is in agreementwith the actual observation in project area.The result also show that using the tensile strain tensor to measure damage and establish damage criterion grasps the nature of hard rock damage.

thermodynamics of irreversible process；elastoplastic；anisotropic damage；FLAC3D；constitutivemodel

TU452

A

1001－5485（2013）12－0048－06

10.3969／j.issn.1001－5485.2013.12.009

2012－12－26；

2013－01－24

楊云浩（1977－），男，山西黎城人，工程師，博士，從事巖石力學(xué)與地下工程方面的研究，（電話）15928690461（電子信箱）haoyun yang2001＠163.com。