艾 濤 朱珍德 田 源 朱 端

(1. 河海大學(xué) 巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室， 南京 210098; 2. 河海大學(xué) 江蘇省巖土工程技術(shù)工程研究所， 南京 210098)

0 引 言

工程巖體內(nèi)部蘊(yùn)含著大量的節(jié)理、裂隙，在外部荷載的擾動(dòng)下巖體內(nèi)裂隙將會(huì)逐漸擴(kuò)展及貫通，從而引起整個(gè)巖體的破壞．因此，在細(xì)觀層次對(duì)巖石進(jìn)行研究可以很好地揭示巖石在外部荷載作用下細(xì)觀損傷斷裂過(guò)程的破壞機(jī)理．計(jì)算機(jī)技術(shù)的發(fā)展和數(shù)值模擬方法的成熟，給復(fù)雜巖石結(jié)構(gòu)變形與破壞細(xì)觀分析提供了有效的方法，而它計(jì)算速度快、精度高等優(yōu)點(diǎn)也使得其在解決土木工程問(wèn)題上得到了廣泛的應(yīng)用．

T.X.Bai等[1]利用有限元軟件FRANC對(duì)裂隙巖體中影響裂紋擴(kuò)展規(guī)律的幾個(gè)重要因素如彈性模量、預(yù)制裂紋深度等進(jìn)行了數(shù)值研究．朱珍德等[2]利用掃描電鏡在巖石的細(xì)觀結(jié)構(gòu)參數(shù)量化方面進(jìn)行了大量的研究并編制了參數(shù)提取的程序．唐春安等開(kāi)發(fā)了巖石破壞過(guò)程分析系統(tǒng)RFPA(Rock Failure Process Analysis)[3]，在此系統(tǒng)基礎(chǔ)上進(jìn)行了巖石三軸破裂過(guò)程數(shù)值模擬[4]．魏玉寒等[5]基于擴(kuò)展有限元法研究了三維含裂隙巖體中裂隙傾角對(duì)于表面裂紋擴(kuò)展的影響．另外還有一些典型的數(shù)值模型，如格構(gòu)(Lattice)模型[6]、顆粒(Granular)模型[7]、梁一顆粒模型[8]等等，這些模型引入介質(zhì)的非均勻性，對(duì)研究非均勻性帶來(lái)的材料破裂過(guò)程的復(fù)雜性頗有幫助．現(xiàn)存的巖石受荷破壞理論中，傳統(tǒng)的宏觀經(jīng)典彈塑性理論力學(xué)模型忽略了巖石細(xì)觀結(jié)構(gòu)的非均勻因素，而細(xì)觀損傷統(tǒng)計(jì)模型雖然考慮了非均勻的影響，但依然不足以表達(dá)整個(gè)破壞過(guò)程的細(xì)觀損傷復(fù)雜性．要深入研究巖石的破壞機(jī)理需要將宏、細(xì)觀損傷破壞問(wèn)題結(jié)合起來(lái)，建立多尺度的力學(xué)模型．在非均勻性研究的基礎(chǔ)上，通過(guò)考慮宏、細(xì)觀結(jié)合的數(shù)值分析方法將有助于更深入地了解巖石受荷破壞機(jī)理．

本文考慮巖石介質(zhì)的非均勻性，結(jié)合細(xì)觀統(tǒng)計(jì)損傷模型的思想，采用Weibull分布來(lái)描述離散后巖石細(xì)觀單元的力學(xué)性質(zhì)參數(shù)，并應(yīng)用Monte-Carlo方法編制程序進(jìn)行單元力學(xué)參數(shù)的賦值．在此基礎(chǔ)上，通過(guò)ABAQUS有限元軟件的二次開(kāi)發(fā)實(shí)現(xiàn)對(duì)單軸壓縮荷載作用下巖石細(xì)觀單元的力學(xué)參數(shù)的影響的模擬研究．

1 有限元法數(shù)值模型

1.1 ABAQUS軟件平臺(tái)以及用戶子程序接口

ABAQUS作為當(dāng)前國(guó)際上最為通用的有限元計(jì)算分析軟件之一[9]，其非線性力學(xué)分析功能具有獨(dú)特的優(yōu)勢(shì)，從而在土木、水利、機(jī)械等領(lǐng)域得到了廣泛應(yīng)用．盡管如此，國(guó)內(nèi)巖土工程領(lǐng)域比較通用的一些本構(gòu)模型，現(xiàn)有的ABAQUS版本卻并不具備，這無(wú)疑限制了ABAQUS在國(guó)內(nèi)巖土工程領(lǐng)域的應(yīng)用與推廣．但是ABAQUS提供了二次開(kāi)發(fā)的接口，通過(guò)其強(qiáng)大的分析求解平臺(tái)[10]，可使困難的分析簡(jiǎn)單化，使復(fù)雜的過(guò)程層次化．

ABAQUS為了方便用戶針對(duì)自己研究領(lǐng)域開(kāi)發(fā)本構(gòu)模型，采用了FORTRAN語(yǔ)言接口，提供了若干用戶子程序(User Subroutines)和在編程時(shí)可以調(diào)用的實(shí)用工具(Utility Routines)[11-12]．常用的接口程序有：用于定義不同物理場(chǎng)的自定義場(chǎng)變量子程序USDFLD、用于定義材料力學(xué)行為的用戶材料子程序UMAT、用于定義不同單元模式的用戶單元子程序UEL等，對(duì)于材料、單元、算法等各方面均有涉及．對(duì)于細(xì)觀單元模型，即可利用用戶子程序USDFLD將表征單元損傷程度值定義為場(chǎng)變量，而場(chǎng)變量的值依賴于積分點(diǎn)上應(yīng)力、應(yīng)變狀態(tài)，在分析步的每一增量步中每一積分點(diǎn)都將調(diào)用USDFLD子程序進(jìn)行不斷的迭代更新．

1.2 單元力學(xué)性質(zhì)的賦值

在前人的研究基礎(chǔ)上，采用Weibull分布來(lái)描述離散后巖石細(xì)觀單元的力學(xué)性質(zhì)參數(shù)，其分布密度函數(shù)為：

(1)

式中，u為細(xì)觀單元體的力學(xué)性質(zhì)參數(shù)(強(qiáng)度、彈性模量等)；u0為與細(xì)觀單元力學(xué)性質(zhì)的平均值有關(guān)的參數(shù)，其數(shù)值并非該平均值；m為分布函數(shù)的形狀參數(shù)，其物理意義反映了巖石介質(zhì)的均勻性，即均質(zhì)度．

基于Monte-Carlo方法和統(tǒng)計(jì)描述相結(jié)合的方法，通過(guò)編制程序?qū)崿F(xiàn)對(duì)單元力學(xué)性質(zhì)(彈性模量和強(qiáng)度)的初始化賦值．以彈性模量的賦值進(jìn)行說(shuō)明，設(shè)所有單元的彈性模量的平均值為E0，則基于式(1)彈性模量weibull分布函數(shù)的積分為：

(2)

其反函數(shù)為：

(3)

1.3 數(shù)值模型的建立

數(shù)值試樣采用尺寸為50 mm×100 mm二維模型，如圖1所示，試樣加載方式為以剛性體控制的位移加載模式，在位移加載的分析步中，采用100個(gè)固定的增量步．?dāng)?shù)值模型采用CPS單元(平面應(yīng)力單元)，考慮到計(jì)算精度和計(jì)算資源的協(xié)調(diào)問(wèn)題，進(jìn)一步選用CPS4R，即細(xì)網(wǎng)格劃分的二維平面線性減縮積分單元進(jìn)行計(jì)算，其常用于大應(yīng)變分析模擬．

圖1 巖石單軸壓縮數(shù)值模型

建模過(guò)程中先運(yùn)用ABAQUS前處理軟件對(duì)巖石細(xì)觀損傷數(shù)值模型的幾何尺寸、單元?jiǎng)澐?、邊界條件等建立初始模型；然后根據(jù)數(shù)值模型文件(INP文件)的特點(diǎn)，通過(guò)編制程序，生成分別包括單元的信息、單元集合的信息、材料參數(shù)的信息的文本文件，并用參數(shù)嵌入INP文件中；接著利用Fortran語(yǔ)言編制USDFLD子程序和相關(guān)程序進(jìn)行各細(xì)觀單元力學(xué)性質(zhì)的賦值、修改，通過(guò)調(diào)用Utility Routines中GETVRM函數(shù)，就可以利用強(qiáng)度準(zhǔn)則進(jìn)行判斷，對(duì)單元的力學(xué)參數(shù)進(jìn)行折減以進(jìn)行下一增量步的計(jì)算；最后保存模型以待運(yùn)行分析．

2 數(shù)值模型的力學(xué)參數(shù)研究

數(shù)值模型中涉及到單元的大小以及細(xì)觀單元的力學(xué)參數(shù)包括：彈性模量和強(qiáng)度的均值E0、fc 0和均質(zhì)度m，壓拉比χ，殘余強(qiáng)度系數(shù)λ．對(duì)于細(xì)觀單元的力學(xué)參數(shù)，巖石破壞細(xì)觀試驗(yàn)研究提供了破壞的細(xì)觀機(jī)理，但由于目前巖石細(xì)觀結(jié)構(gòu)力學(xué)性質(zhì)的試驗(yàn)研究還很少，細(xì)觀單元的力學(xué)參數(shù)還很難通過(guò)細(xì)觀試驗(yàn)獲取，只有通過(guò)大量的數(shù)值模型研究來(lái)對(duì)細(xì)觀力學(xué)參數(shù)進(jìn)行標(biāo)定．如沒(méi)有特別說(shuō)明，數(shù)值模擬中材料泊松比μ和摩擦角θ分別采用0.25和30°．設(shè)置不同單元數(shù)目進(jìn)行數(shù)值破壞形態(tài)模擬，圖2為采用1 mm×1 mm的單元尺寸(5 000個(gè)單元)下的數(shù)值試樣破壞形態(tài)圖，從破壞形態(tài)上看，該單元?jiǎng)澐挚梢詽M足數(shù)值模擬要求．

圖2 5 000個(gè)單元數(shù)目下數(shù)值試樣的破壞形態(tài)

2.1 均質(zhì)度m對(duì)宏觀的影響

為得出m對(duì)宏觀應(yīng)力應(yīng)變曲線的影響，建立均質(zhì)度分別為1.5、3、6的模型，具體參數(shù)見(jiàn)表1，結(jié)果如圖3所示．由圖3可見(jiàn)，數(shù)值試樣的宏觀響應(yīng)對(duì)于細(xì)觀單元整體均質(zhì)度的變化很敏感，隨著m的增大，試樣的脆性明顯隨之增大，宏觀的彈性模量和強(qiáng)度也相應(yīng)提高．m越大，數(shù)值試樣越均勻，宏觀力學(xué)參數(shù)越接近細(xì)觀參數(shù)，隨著m的增大，曲線的峰值應(yīng)變也隨之增大．

表1 不同單元尺寸的數(shù)值試樣

圖3 不同m下的應(yīng)力應(yīng)變曲線

2.2 宏細(xì)觀彈性模量和強(qiáng)度的關(guān)系

為了得出一定均質(zhì)度下，數(shù)值試樣的細(xì)觀單元的E0，fc 0與宏觀彈性模量和宏觀強(qiáng)度的關(guān)系，建立的數(shù)值試樣以及宏觀彈性模量和強(qiáng)度的結(jié)果見(jiàn)表2．

表2 不同細(xì)觀彈性模量和強(qiáng)度下數(shù)值試樣及其結(jié)果

對(duì)m=1.5和3下的數(shù)值試驗(yàn)結(jié)果與細(xì)觀單元的強(qiáng)度和彈性模量的擬合結(jié)果如圖4～5所示，可見(jiàn)細(xì)觀單元強(qiáng)度、彈性模量和宏觀強(qiáng)度、彈性模量的擬合具有很好的線性關(guān)系，可以得到：

(4)

(5)

式(4)(5)中，fc 0和E0為Weibull分布下細(xì)觀單元強(qiáng)度和彈性模量的均值，fc和E為數(shù)值試樣宏觀的強(qiáng)度和彈性模量．

圖4 細(xì)觀單元強(qiáng)度與宏觀強(qiáng)度的擬合曲線

圖5 細(xì)觀單元彈性模量和宏觀彈性模量的擬合曲線

2.3 壓拉比以及殘余強(qiáng)度系數(shù)對(duì)宏觀的影響

1)壓拉比的影響

巖石細(xì)觀單元的的壓拉比χ反映了巖石的抗壓和抗拉能力的差異，不同壓拉比系數(shù)下數(shù)值試樣的應(yīng)力應(yīng)變曲線如圖6所示．所有試樣的應(yīng)力應(yīng)變曲線都經(jīng)歷了線性變形、非線性變形、不斷劣化和殘余變形4個(gè)階段．隨著χ的降低，數(shù)值試樣的峰值強(qiáng)度逐漸提高，所有曲線都顯示相同線性段，殘余強(qiáng)度在χ>5時(shí)，隨著χ的降低有所提高但不明顯，χ=3時(shí)殘余強(qiáng)度提高很明顯．對(duì)于相同的單元壓縮強(qiáng)度f(wàn)c 0=300 MPa，不同的壓拉比對(duì)應(yīng)不同的單元受拉強(qiáng)度．

圖6 不同壓拉比數(shù)值試樣應(yīng)力應(yīng)變曲線

對(duì)于本文中摩擦角為30°，當(dāng)χ=3時(shí)，如果單元發(fā)生受拉破壞，由于σ3≤ft=-fc 0/3?-3σ3≥fc 0?σ1-3σ3≥fc 0，可見(jiàn)，當(dāng)χ≤3時(shí)，單元如果滿足最大拉應(yīng)力準(zhǔn)則則必然滿足摩爾庫(kù)倫準(zhǔn)則，即單元在受剪破壞一定發(fā)生在受拉破壞之前．之所以隨著壓拉比系數(shù)的降低，數(shù)值試樣的宏觀強(qiáng)度逐漸增加，是因?yàn)閴豪仍叫。?xì)觀單元的剪切破壞越占主導(dǎo)，受拉破壞的單元數(shù)目越少．

2)殘余強(qiáng)度系數(shù)的影響

為研究細(xì)觀單元的殘余強(qiáng)度系數(shù)λ對(duì)宏觀應(yīng)力力學(xué)性質(zhì)的影響，分別建立殘余強(qiáng)度系數(shù)為0.1，0.3，0.5的數(shù)值試件．當(dāng)壓拉比為5時(shí)，3個(gè)數(shù)值試樣的應(yīng)力應(yīng)變曲線如圖7所示．

圖7 不同細(xì)觀殘余強(qiáng)度下的應(yīng)力應(yīng)變曲線

由圖可見(jiàn)，隨著細(xì)觀單元?dú)堄鄰?qiáng)度系數(shù)的提高，試樣的峰值強(qiáng)度有一定的增加，試樣的延性有一定的增強(qiáng)，數(shù)值試樣宏觀的破壞后殘余強(qiáng)度并沒(méi)有增加，幾乎保持不變，略有降低．這主要是因?yàn)闅堄鄰?qiáng)度系數(shù)只是對(duì)受剪破壞的單元起作用，受拉破壞的單元幾乎沒(méi)有承載力，所以壓拉比越小，單元強(qiáng)度殘余系數(shù)的影響越大．雖然在峰值應(yīng)力附近，承載力有一定的提高，但如果受剪破壞的單元受力狀態(tài)最終變?yōu)槔炱茐臅r(shí)，單元將不再具有承載力．因此，試樣最后宏觀殘余強(qiáng)度都很小，變化不大，甚至隨著λ的提高，試樣最終破壞后的宏觀殘余強(qiáng)度略有降低．

3 巖石單軸壓縮過(guò)程的數(shù)值模擬

通過(guò)對(duì)數(shù)值試驗(yàn)中細(xì)觀單元參數(shù)的研究可知，數(shù)值試樣細(xì)觀參數(shù)的選取取決于宏、細(xì)觀的試驗(yàn)．為研究巖石在單軸荷載下的破壞過(guò)程，結(jié)合武沂泉[13]的細(xì)觀壓縮SEM試驗(yàn)及其宏觀的應(yīng)力應(yīng)變曲線建立對(duì)應(yīng)的數(shù)值試樣．

武沂泉[13]利用采自四川錦屏二級(jí)水電站工程現(xiàn)場(chǎng)的大理巖作為試樣進(jìn)行單軸加載．使用SEM對(duì)試件采用分層式掃描拍攝以采集試件在各級(jí)荷載下的SEM圖像試件在各級(jí)荷載下SEM圖像，并通過(guò)圖像處理提取了裂紋信息，得到了細(xì)觀裂紋隨應(yīng)力加載的變化規(guī)律及應(yīng)力應(yīng)變曲線．由宏觀應(yīng)力應(yīng)變曲線得出的宏觀的峰值強(qiáng)度為77.46 MPa，峰值應(yīng)變?yōu)?.003 4，彈性模量為30 571 MPa，根據(jù)公式(4)或(5)可得到細(xì)觀單元的強(qiáng)度f(wàn)c 0=902.8 MPa和彈性模量E0=50 818.7 MPa．綜合考慮其它參數(shù)的影響并通過(guò)試算，最終單元的細(xì)觀參數(shù)情況確定見(jiàn)表3．

表3 壓縮數(shù)值試樣的參數(shù)表

不同應(yīng)力狀態(tài)下單軸壓縮試驗(yàn)的SEM圖片及數(shù)值試樣破壞情況的試樣破壞過(guò)程及形態(tài)如圖8所示．可見(jiàn)，當(dāng)加載到應(yīng)力20 MPa左右時(shí)，試件表面變化不是很明顯，由于局部初始損傷和細(xì)觀不均勻性的存在，只有少量的局部黑點(diǎn)處出現(xiàn)拉應(yīng)力集中，對(duì)應(yīng)SEM圖片中的白色發(fā)亮處，這些微裂紋方向均與軸向荷載方向一致或成小夾角，且在寬度上都很細(xì)．繼續(xù)加載達(dá)到60 MPa，此時(shí)由掃描圖片可以看出微裂紋數(shù)目和寬度都有所增加，裂紋擴(kuò)展方向依然與軸向荷載幾乎一致，0～60 MPa這一過(guò)程右面數(shù)值試樣的破裂情況處于微破裂均勻隨機(jī)發(fā)展階段，此階段細(xì)觀單元的破壞是隨機(jī)出現(xiàn)并比較均勻分布在試樣中．當(dāng)應(yīng)力超過(guò)60 MPa，數(shù)值試樣的破裂情況進(jìn)入微破裂非均勻發(fā)展階段，單元的破壞由隨機(jī)均勻的破壞轉(zhuǎn)為在某些單元附近發(fā)生貫通，并逐漸向集中的某幾個(gè)局部區(qū)域發(fā)展，破壞單元的數(shù)目持續(xù)增多，在對(duì)應(yīng)的左邊掃描圖像中裂紋的寬度急劇增大，裂紋進(jìn)入不穩(wěn)定擴(kuò)展階段(圖8(c))．應(yīng)力應(yīng)變曲線的峰值后下降段對(duì)應(yīng)數(shù)值試樣的局部及宏觀破裂發(fā)展階段，微破裂的局部化趨于顯著，單元的破壞數(shù)目迅速增多，破裂單元迅速貫通，當(dāng)達(dá)到峰值時(shí)，數(shù)值試樣完全破壞，形成破壞的剪切帶(圖8(d))．

圖8 不同應(yīng)力狀態(tài)下SEM圖片及數(shù)值試樣破壞情況

由以上巖石單軸壓縮的真實(shí)試驗(yàn)與數(shù)值模擬的結(jié)果對(duì)照可見(jiàn)，兩者在應(yīng)力應(yīng)變曲線以及破壞形態(tài)和機(jī)理上都存在一定的相關(guān)性，數(shù)值模擬的結(jié)果在一定程度上可以合理地反映巖石單軸荷載作用下應(yīng)力應(yīng)變曲線以及從細(xì)觀破壞直至宏觀破壞的整個(gè)過(guò)程．

單軸壓縮下數(shù)值試驗(yàn)得到的應(yīng)力應(yīng)變?nèi)€以及對(duì)應(yīng)的SEM試驗(yàn)試樣曲線如圖9所示，可見(jiàn)數(shù)值試樣和試驗(yàn)試樣的應(yīng)力應(yīng)變曲線大體上是很接近的，說(shuō)明數(shù)值試樣的建立是合理的，也驗(yàn)證了細(xì)觀參數(shù)對(duì)宏觀試樣的影響．從兩者的曲線對(duì)比來(lái)看，數(shù)值試樣的應(yīng)力應(yīng)變曲線反映了單軸壓縮荷載作用下巖石的非線性變形階段、弱化、以及殘余強(qiáng)度的特征．

圖9 數(shù)值試樣與試驗(yàn)試樣應(yīng)力應(yīng)變?nèi)€

4 結(jié) 論

本文基于巖石介質(zhì)的非均勻性理論以及細(xì)觀力學(xué)的基本思想，對(duì)細(xì)觀單元各主要參數(shù)的影響以及荷載作用下的破裂過(guò)程和機(jī)理進(jìn)行了分析．主要得出以下結(jié)論：

1)隨著單元均質(zhì)度的增大，試樣的脆性增加，數(shù)目越多，越能反映試樣的破壞形態(tài)．

2)細(xì)觀單元和宏觀試樣的彈性模量、強(qiáng)度之間存在很好的線性關(guān)系．

3)隨著細(xì)觀單元的壓拉比的降低、數(shù)值試樣的宏觀強(qiáng)度呈現(xiàn)增大趨勢(shì)，而細(xì)觀單元的殘余強(qiáng)度系數(shù)對(duì)數(shù)值試樣的宏觀的殘余強(qiáng)度影響甚微，隨著其提高，對(duì)應(yīng)力應(yīng)變曲線硬化段有一定的提高．

4)巖石單軸壓縮數(shù)值模型結(jié)果表明數(shù)值試樣和試驗(yàn)試樣的應(yīng)力應(yīng)變曲線和破壞情況大體上是很接近的，說(shuō)明數(shù)值試樣的建立是合理的．