塑性擴(kuò)容梯度對(duì)孔隙介質(zhì)巖體中T-H-M耦合作用的二維有限元分析

張玉軍，張維慶

（1.中國(guó)科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430071；2.中鐵隧道勘測(cè)設(shè)計(jì)院有限公司，天津 300133）

1 引言

對(duì)于承受壓力作用的巖體，隨著荷載的增加，其體積變形要經(jīng)歷壓縮、擴(kuò)容（剪脹）、潰散等階段。關(guān)于巖體的擴(kuò)容，人們已進(jìn)行了相當(dāng)多的研究[1－10]。金濟(jì)山[11]研究了擴(kuò)容應(yīng)變和擴(kuò)容體應(yīng)變的變化，并給出了本構(gòu)方程，用該理論做出的數(shù)值模擬結(jié)果與巖石力學(xué)試驗(yàn)的結(jié)果比較吻合。侯文詩(shī)等[12]通過(guò)對(duì)8種巖樣的單軸、三軸壓縮試驗(yàn)結(jié)果分析，獲得低圍壓下的應(yīng)力–體積應(yīng)變曲線，揭示了巖石擴(kuò)容起始特性與峰值點(diǎn)特性的對(duì)應(yīng)關(guān)系。段艷燕等[13]認(rèn)為，巖石的剪脹非微觀結(jié)構(gòu)的產(chǎn)生、擴(kuò)展和匯集，其實(shí)質(zhì)為巖石內(nèi)部組構(gòu)特征發(fā)生了顯著的變化，是破裂塊體之間鑲嵌組合的一種結(jié)構(gòu)效應(yīng)，剪脹變形要比質(zhì)點(diǎn)連續(xù)移動(dòng)的彈塑性變形量大得多。閻巖等[14]通過(guò)滲流-流變耦合試驗(yàn)顯示，在低應(yīng)力作用下試件體積不斷縮小，滲透系數(shù)也隨之減小，應(yīng)力增大到某值后體積逐步擴(kuò)大而出現(xiàn)擴(kuò)容，滲透系數(shù)也隨之增大。張玉軍等[15]在所建立的熱–氣–應(yīng)力(T-H-M)耦合彈塑性模型中，嘗試將巖體的孔隙率和滲透性作為平均有效應(yīng)力和剪脹體積應(yīng)變的函數(shù)，以一個(gè)假定的地質(zhì)結(jié)構(gòu)為考察對(duì)象，對(duì)CO2注入過(guò)程中巖體內(nèi)的熱-氣-應(yīng)力耦合現(xiàn)象進(jìn)行了有限元計(jì)算分析。

然而，迄今為止對(duì)于巖體的擴(kuò)容及相伴效應(yīng)（如滲透系數(shù)增加）的研究還主要限于室內(nèi)試驗(yàn)和理論建模，在數(shù)值模擬方面的應(yīng)用尚相當(dāng)少，故相應(yīng)的現(xiàn)象和規(guī)律的揭示還很不充分。鑒于此，筆者以一個(gè)假定的位于飽和單一孔隙地層中的高放廢物地質(zhì)處置庫(kù)實(shí)驗(yàn)室尺度模型為算例，在相同的初始溫度和巖體應(yīng)力條件下，基于Mohr-Coulomb屈服準(zhǔn)則，針對(duì)考慮塑性擴(kuò)容梯度對(duì)巖體孔隙率及滲透系數(shù)影響不同的5種計(jì)算工況，進(jìn)行熱-水-應(yīng)力耦合有限元數(shù)值模擬，考察了處置庫(kù)近場(chǎng)的溫度、孔隙率和滲透系數(shù)、孔隙水壓力、地下水流速、應(yīng)力的變化、分布情況，得出了一些新的認(rèn)識(shí)。

2 孔隙介質(zhì)熱-水-應(yīng)力耦合方程

對(duì)于單一孔隙介質(zhì)，筆者建立了一種熱-水-應(yīng)力耦合模型(由相應(yīng)的雙重孔隙—裂隙介質(zhì)模型退化而來(lái))[16]，控制方程如下。

（1）應(yīng)力平衡方程

應(yīng)力平衡方程為式中：σ、ε 分別為應(yīng)力和應(yīng)變；D=C-1為彈性矩陣；mT=[1 1 0]為法向應(yīng)力單位列陣；Ks、βs、T 分別為孔隙介質(zhì)的體積模量、熱膨脹系數(shù)和溫度；sw、pw、Ds、C 分別為孔隙介質(zhì)的飽和度、水壓力、濕氣容量和柔度矩陣；t為時(shí)間。

（2）水連續(xù)性方程

根據(jù)質(zhì)量守恒原理，在dt 時(shí)段內(nèi)流入某一物體的水量應(yīng)等于其內(nèi)部?jī)?chǔ)水量的增加。設(shè)定水的滲流可以用達(dá)西定律來(lái)描述，則有

式中：k、krw分別為孔隙介質(zhì)的固有滲透率張量和比滲透率；ρw、μw、γw依次為水的密度、黏滯系數(shù)和重度；z為位置水頭；Dt為孔隙介質(zhì)的溫度梯度水分?jǐn)U散系數(shù)；A、B、F 均為常數(shù)項(xiàng)表達(dá)式。

（3）能量守恒方程

根據(jù)能量守恒原理，在dt 時(shí)段內(nèi)流入某一物體的熱量應(yīng)等于其內(nèi)能的增加，可得

式中：Cw、Cs分別為水及孔隙介質(zhì)的比熱；φ、ρs、λ 分別為孔隙介質(zhì)的孔隙率、密度和導(dǎo)熱系數(shù)矩陣；Va為孔隙水的平均表觀流速； ui、 uj為位移分量；δij為克羅內(nèi)克符號(hào)。

3 彈塑性分析

彈塑性計(jì)算時(shí)有效應(yīng)力增量可表示為

式中：σ′、ε 分別為效應(yīng)力列陣和應(yīng)變列陣；De為彈性矩陣；A為硬化參數(shù)；Q、F 分別為塑性勢(shì)及屈服函數(shù)。

當(dāng)Q=F 時(shí)，稱為相關(guān)聯(lián)的流動(dòng)法則；當(dāng)Q ≠F時(shí)，稱為不相關(guān)聯(lián)的流動(dòng)法則。

本程序中使用相關(guān)聯(lián)的流動(dòng)法則，并令A(yù)=0，取Mohr-Coulomb準(zhǔn)則為屈服函數(shù)[[17]，即

其中：

4 孔隙率-滲透系數(shù)演化律

圖1、2為根據(jù)試驗(yàn)得出的巖石應(yīng)力-應(yīng)變?nèi)^(guò)程曲線及體積應(yīng)變與滲透系數(shù)的關(guān)系。圖1中，A、B、C 依次對(duì)應(yīng)于為軸向、徑向及體積應(yīng)變）。

圖1 應(yīng)力-應(yīng)變?nèi)^(guò)程曲線[12]Fig.1 Full relationships between stress and strain[12]

圖2 體積應(yīng)變-滲透系數(shù)的關(guān)系[14]Fig.2 Relationships between volumetric strain and permeability coefficient[14]

需要指出，本文中的“塑性擴(kuò)容”是指巖體的“破裂”或“破碎”階段，此時(shí)實(shí)際巖體已是非連續(xù)介質(zhì)，其內(nèi)部組構(gòu)特征發(fā)生了顯著的變化，剪脹變形要比質(zhì)點(diǎn)連續(xù)移動(dòng)的彈塑性變形量大得多[13]，用塑性力學(xué)（連續(xù)介質(zhì)、剪脹角等概念）已不能很好地描述其特性，故引用和擴(kuò)展了馬士進(jìn)提出的巖體擴(kuò)容與最大主應(yīng)變的線性關(guān)系式，這是一種經(jīng)驗(yàn)處理。

馬士進(jìn)[18]提出，在單軸壓縮條件下（圍壓σ3=0）巖體擴(kuò)容與最大主應(yīng)變可以簡(jiǎn)化成線性關(guān)系：

式中：ΔV、V 分別為巖體的擴(kuò)容體積增量和初始體積；ε1為最大主應(yīng)變；η為擴(kuò)容梯度。

馬士進(jìn)根據(jù)Price的實(shí)驗(yàn)，認(rèn)為巖體擴(kuò)容梯度并非常數(shù)，而是和巖體組成、側(cè)壓系數(shù)（σ2、σ3不為0）等有關(guān)。但在計(jì)算中，為簡(jiǎn)化問(wèn)題，馬士進(jìn)取η為定值[18]。筆者在本文的數(shù)值模擬中，也取η為常數(shù)，但考察η 的變化對(duì)計(jì)算結(jié)果的影響。

當(dāng)假設(shè)巖體的固體顆粒為不可壓縮時(shí)，近似地有Δεvk=Δφ[19]，從而將Taron等[20]給出的孔隙率表達(dá)式修正為

式中：εv為彈塑性體積應(yīng)變；φ0為初始孔隙率。

由此可求出對(duì)應(yīng)的滲透系數(shù)[21]為

式中：k0為初始時(shí)刻的滲透系數(shù)。

將上述單一孔隙介質(zhì)熱-水-應(yīng)力耦合模型及孔隙率-滲透系數(shù)演化律引入筆者所研制的二維有限元程序中[22]，其計(jì)算步驟為

①開(kāi)始時(shí)，孔隙率、滲透系數(shù)和其他參數(shù)取初值。

②經(jīng) ti=Δt 時(shí)段的彈性或彈塑性T-H-M耦合分析后，得出相應(yīng)的溫度場(chǎng)、滲流場(chǎng)和應(yīng)力場(chǎng)。

③針對(duì)塑性巖體單元，計(jì)算擴(kuò)容體積應(yīng)變，求出新的孔隙率、滲透系數(shù)和其他參數(shù)。

5 算 例

如圖3所示，此為一實(shí)驗(yàn)室尺度的核廢料玻璃固化體埋置模型，其周圍的介質(zhì)是飽和的孔隙巖體。作為近似簡(jiǎn)化，認(rèn)為這是一個(gè)平面應(yīng)變問(wèn)題。取計(jì)算域的范圍為水平向4 m，垂直向8 m，有800個(gè)單元，861個(gè)節(jié)點(diǎn)。玻璃固化體為矩形，其尺寸為0.4 m×0.4 m、中心的坐標(biāo)x=y=0（y 軸向下為正）。從固化體邊緣向右的點(diǎn)號(hào)依次為432、433、434、435、436。

圖3 計(jì)算模型Fig.3 Computation model

對(duì)于邊界條件，計(jì)算域的頂面位移自由，其上作用有分布荷載53.4 MPa，左、右側(cè)面的水平方向位移約束，底面的垂直方向位移約束，所有邊界的溫度為40 ℃，各邊界的孔隙水壓力均為4.59 MPa，巖體的初始孔隙率為0.1。熱-水-應(yīng)力耦合的環(huán)境對(duì)巖體介質(zhì)要產(chǎn)生塑性擴(kuò)容作用，從而引起孔隙率和滲透系數(shù)的變化。有關(guān)的計(jì)算參數(shù)見(jiàn)表1。取玻璃固化體為彈性體，而巖體的黏聚力和內(nèi)摩擦角分別為c=12.0 MPa和φ=35°。初始狀態(tài)時(shí)，巖體的溫度均為40 ℃。核廢物以1 000 W的不變功率釋放熱量，時(shí)間經(jīng)歷了4年。

巖體介質(zhì)的水分特性曲線符合Van Genuchten模型[23]，即

式中：α=3.86×10-6m-1，β=1.41，γ=1 -1/β；ψ為水勢(shì)；sws=1.0，swr=0.19，其分別為最大飽和度和最小飽和度。

表1 主要計(jì)算參數(shù)Table 1 Main computation parameters

比滲透率與飽和度的關(guān)系為

參照文獻(xiàn)[18]中的參數(shù)敏感性分析，設(shè)定了5種工況：擴(kuò)容梯度η 依次為0.0、1.0、1.5、2.0、2.5。主要結(jié)果及分析如下:由于滲流和地應(yīng)力對(duì)核廢物釋熱的影響較小[24]，故5種工況條件下計(jì)算域中的溫度場(chǎng)基本相同。以工況1為例，432、433、434、435各點(diǎn)處的溫度隨時(shí)間的變化曲線如圖4所示。由圖可見(jiàn)，開(kāi)始約0.1年內(nèi)巖體的溫度快速上升，之后增加減緩，到計(jì)算終了時(shí)432、433、434、435各點(diǎn)的溫度依次為97.8、81.9、72.6 ℃和65.7 ℃。圖5為工況1在4年時(shí)計(jì)算域中的溫度等值線分布。

由圖6可見(jiàn)，各工況的計(jì)算域中應(yīng)力量值及分布產(chǎn)生了一定差別，如到4年時(shí)432、434點(diǎn)水平正應(yīng)力(MPa)/垂直正應(yīng)力(MPa)依次為，工況1：-10.334/-32.079、-14.484/-35.332；工況2：-10.386/32.100、-14.515/-35.362；工況3：-10.455/-32.084、-14.553/-35.376；工況4：-10.575/-32.032、-14.610/35.384；工況5：-10.710/-31.929、-14.677/-35.388。到4年時(shí)巖體中塑性區(qū)見(jiàn)圖7，其分布大致呈豎立的啞鈴形，而工況1～5的塑性區(qū)面積依次為12.28、13.28、13.56、14.20 m2和14.92 m2。由此可知，隨著巖體擴(kuò)容梯度η 的增加，巖體中的應(yīng)力量值及塑性區(qū)面積均上升。

圖4 溫度-時(shí)間曲線Fig.4 Curves of temperatures vs.time

圖5 4年時(shí)工況1計(jì)算域中溫度等值線（單位：℃）Fig.5 Temperature contours in calculation domain at 4 years(unit:℃)

圖8為到4年時(shí)巖體中孔隙率增量(φ4y-φ0)等值線分布。由圖可見(jiàn)，相比于不考慮擴(kuò)容的工況1，其余考慮擴(kuò)容的工況有很大的差別不同，特別是在彈、塑區(qū)的交界處有顯著的等值線集中現(xiàn)象，具有某種“剪切帶效應(yīng)”，并且考慮擴(kuò)容的各工況也依各自的擴(kuò)容梯度不同，其孔隙率增量等值線的量值和分布也有所變化。此時(shí)在固化體中心正下方y(tǒng)為0.3 m和1.1 m處的孔隙率對(duì)于工況1～5依次為：0.100 535、0.100 687，0.100 537、0.112 463 9，0.100 539、0.137 618，0.100 541、0.149 020，0.100 541、0.163 456，均比初始值0.1升高。但因y為0.3 m和1.1 m處分別位于彈性區(qū)和塑性區(qū)，只有在塑性區(qū)才產(chǎn)生體積擴(kuò)容并計(jì)入其對(duì)孔隙率的影響，故塑性區(qū)的孔隙率相比于其初始值有明顯增長(zhǎng)，并且隨所取擴(kuò)容梯度值的變大，各工況的孔隙率也逐漸上升。上述兩點(diǎn)處的孔隙率增量隨時(shí)間的變化曲線如圖9所示。

圖6 4年時(shí)計(jì)算域中正應(yīng)力等值線（單位：MPa）Fig.6 Normal stress contours in calculation domain at four years(unit:MPa)

圖7 4年時(shí)計(jì)算域中塑性區(qū)Fig.7 Plastic zones in calculation domain at four years

圖8 4年時(shí)巖體中孔隙率增量等值線（工況1，單位：10-4，工況2～5，單位：10-2）Fig.8 Contours of porosity increment in rock mass at four years（unit is 10-4 for case 1，and 10-2 for cases 2-5）

圖9 巖體中2個(gè)點(diǎn)的孔隙率增量與時(shí)間曲線(固化體中心y=0)Fig.9 Porosity increments vs.time at two nodes(y=0 at the center of vitrified waste)

由圖9可知，開(kāi)始約0.6年內(nèi)，因玻璃固化體的急劇釋放熱量和巖體上覆及自重荷載的施加，使得巖體孔隙率增量?jī)?nèi)快速上升，之后隨時(shí)間的推移孔隙率增量雖有所下降，但變化幅度較?。u趨平衡），并且對(duì)于不同的擴(kuò)容梯度值，孔隙率增量在彈性區(qū)變化不大，而在塑性區(qū)卻有明顯差別。

相對(duì)應(yīng)的滲透系數(shù)增量的等值線分布及隨時(shí)間的變化曲線分別見(jiàn)圖10、11。其與孔隙率增量的演化規(guī)律類似，在時(shí)間上表現(xiàn)了一定的非線性特征。上述兩點(diǎn)處的滲透系數(shù)對(duì)于工況1～5依次為：1.261 51×10-13、1.267 01×10-13m/s；1.261 61×10-13、2.538 00×10-13m/s；1.261 67×10-13、3.463 31×10-13m/s；1.261 75×10-13、4.533 29×10-13m/s；1.261 77×10-13、5.994 23×10-13m/s。由圖中也看到，固化體中心正下方y(tǒng)=1.1 m處進(jìn)入塑性區(qū)之后該點(diǎn)滲透系數(shù)也比彈性時(shí)有明顯的增加（體積擴(kuò)容的作用）。

圖10 4年時(shí)巖體中滲透系數(shù)增量等值線（工況1，單位：10-15 m/s，工況2～5的單位：10-13 m/s）Fig.10 Contours of permeability increment in rock mass at four years（unit is 10-15 m/s for case 1，and 10-13 m/s for cases 2-5）

圖11 巖體中2個(gè)點(diǎn)的滲透系數(shù)增量與時(shí)間曲線(固化體中心y=0)Fig.11 Permeability increments vs.time at two nodes(y=0 at the center of vitrified waste)

圖12為在固化體中心正下方y(tǒng)=0.3、1.1 m處的孔隙水壓力隨時(shí)間的變化曲線。

圖12 巖體中兩個(gè)點(diǎn)的孔隙水壓力-時(shí)間曲線Fig.12 Pore pressures vs.time at two nodes for two cases

由圖12可見(jiàn)，一開(kāi)始由于巖體荷載的突然施加，孔隙水壓力有瞬時(shí)的少許下降，之后在應(yīng)力場(chǎng)和溫度場(chǎng)的共同作用下，孔隙水壓力隨時(shí)間呈現(xiàn)一定幅度的上升。到4年時(shí)該兩點(diǎn)處的孔隙水壓力值對(duì)于工況1～5為：4.587 01、4.587 02，4.589 31、4.590 40，4.591 22、4.592 64，4.593 45、4.595 37，4.595 92、4.597 88 MPa。彈性區(qū)和塑性區(qū)的孔隙水壓力無(wú)明顯差別，其原因在于彈性區(qū)和塑性區(qū)是水力連通的，兩區(qū)的孔隙水壓力隨時(shí)間趨于平衡。同時(shí)也比較可知，隨所取擴(kuò)容梯度值的變大，各工況的孔隙水壓力也有所上升（僅管幅度很小）。計(jì)算終了時(shí)兩種工況中的孔隙水壓力等值線如圖13所示。

圖13 4年時(shí)巖體中的孔隙水壓力等值線(單位：kPa)Fig.13 Contours of pore pressures in rock mass at four years(unit:kPa)

由圖13可見(jiàn)，大致以過(guò)玻璃固化體中心的水平線為分界，在其上、下方形成了2個(gè)孔隙水壓力等值線族閉合環(huán)，與啞鈴形的塑性區(qū)對(duì)應(yīng)。并且隨所取擴(kuò)容梯度值的改變，各工況的孔隙水壓力等值線的量值和分布亦呈現(xiàn)一定的差別。

圖14為各工況在4年時(shí)計(jì)算域中的孔隙水流速矢量分布，工況1與其他工況的比例為500:1。由圖可見(jiàn)，其與塑性區(qū)的分布有一定的對(duì)應(yīng)關(guān)系，并且隨著所取擴(kuò)容梯度值的增加，孔隙水流速矢量也逐漸變大。在固化體中心正下方y(tǒng)=0.3、1.1 m處，到4年時(shí)該2點(diǎn)處的孔隙水流速矢量對(duì)于工況1～5依次為：6.81×10-12、2.15×10-12，6.18×10-10、1.22×10-10，1.04×10-9、1.08×10-11，1.54×10-9、9.48×10-10，2.04×10-9、6.41×10-9m/s 。

圖14 4年時(shí)計(jì)算域中的孔隙水流速矢量Fig.14 Flow vectors of pore water in calculation domain at four years for two cases

6 結(jié)語(yǔ)

筆者針對(duì)其所建立的雙重孔隙—裂隙介質(zhì)熱-水-應(yīng)力（T-H-M）耦合模型，經(jīng)退化處理，得出了相應(yīng)的單一孔隙介質(zhì)熱-水-應(yīng)力耦合模型，并引入巖體擴(kuò)容梯度，使用Mohr-Coulomb準(zhǔn)則，計(jì)入塑性擴(kuò)容對(duì)巖體孔隙率及滲透系數(shù)影響，以一個(gè)假定的實(shí)驗(yàn)室尺度且位于飽和孔隙巖體中的高放廢物地質(zhì)處置模型為例子，就不同的擴(kuò)容梯度值進(jìn)行了5種工況的熱-水-應(yīng)力耦合的二維有限元模擬，考察了巖體中的溫度、正應(yīng)力、塑性區(qū)、孔隙率及滲透系數(shù)、孔隙水壓力和地下水流速的變化、分布情況。計(jì)算結(jié)果表明，兩種工況中的溫度狀態(tài)基本相同，計(jì)算終了的4年時(shí)，巖體近場(chǎng)的溫度可達(dá)到40.0～98.0 ℃；塑性區(qū)的分布大致呈豎立的啞鈴形；相比于不考慮擴(kuò)容梯度的工況1，其余工況的正應(yīng)力、孔隙率及滲透系數(shù)、孔隙水壓力和地下水流速等的分布與塑性區(qū)的分布有明顯的對(duì)應(yīng)關(guān)系，呈現(xiàn)了某種“剪切帶效應(yīng)”；正應(yīng)力量值、塑性區(qū)面積、孔隙率及滲透系數(shù)、孔隙水壓力、地下水流速等均隨所取擴(kuò)容梯度值的變大而增加。

[1]COOK N G W.An experiment proving that dilatancy is a pervasive volumetric property of brittle rock loaded to failure[J].Rock Mechanics，1970，2:181－188.

[2]陳宗基，康文法.在巖石破壞和地震之前與時(shí)間有關(guān)的擴(kuò)容[J].巖石力學(xué)與工程學(xué)報(bào)，1983，2(1):11－21.TAN Tjiong-kie，KANG Wen-fa.Time dependent dilatancy prior to rock failure and earthquakes[J].Chinese Journal of Rock Mechanics and Engineering，1983，2(1):11－21.

[3]陳宗基，石澤全，于智海，等.用8 000 kN多功能三軸儀測(cè)量脆性巖石的擴(kuò)容、蠕變及松弛[J].巖石力學(xué)與工程學(xué)報(bào)，1989，8(2):97－118.TAN Tjiong-kie，SHI Zhe-quan，YU Zhi-hai，et al.Dilatancy creep and relaxation of brittle rocks measured with the 8 000 kN multipurpose triaxial apparatus[J].Chinese Journal of Rock Mechanics and Engineering，1989，8(2):97－118.

[4]陳宗基，康文法.巖石的封閉應(yīng)力、蠕變和擴(kuò)容及本構(gòu)方程[J].巖石力學(xué)與工程學(xué)報(bào)，1991，10(4):200－312.TAN Tjiong-kie，KANG Wen-fa.On the locked in stress，creep and dilatation of rocks，and the constitutive equations[J].Chinese Journal of Rock Mechanics and Engineering，1991，10(4):200－312.

[5]許東俊，李小春.應(yīng)力狀態(tài)和巖石擴(kuò)容特性[J].巖土力學(xué)，1992，13(2):11－21.XU Dong-jun，LI Xiao-chun.Stress path and dilatancy behavior of rocks[J].Rock and Soil Mechanics，1992，13(2):11－21.

[6]康紅普.巖石擴(kuò)容與巷道底臌[J].礦業(yè)快報(bào)，1992，11(3):15－17.KANE Hong-pu.Rock dilatancy and opening floor heave[J].Express Information of Mining Industry，1992，11(3):15－17.

[7]MAHMUTOGLU Y，VARDAR M.Effects of inelastic volume increase on fractured rock behaviour[J].Bull Eng.Geol.Env.，2003，62:117－121.

[8]YUAN SHINCHE，HARRISON J.An empirical dilatancy index for the dilatant deformation of rock[J].International Journal of Rock Mechanics and Mining Sciences，2004，41:679－686.

[9]李林峰，侯克鵬，余賢斌，等.巖石峰后擴(kuò)容特性及擴(kuò)容模型的探討[J].礦業(yè)快報(bào)，2007，26(9):15－17.LI Lin-feng，HOU Ke-peng，YU Xian-bin，et al.Study of rock dilatancy characteristic and dilatancy model after peak stress[J].Express Information of Mining Industry，2007，26(9):15－17.

[10]張長(zhǎng)科，李林峰.巖石峰值應(yīng)力后擴(kuò)容與圍壓的關(guān)系[J].有色金屬，2009，61(4):134－137.ZHANG Chang-ke，LI Lin-feng.Study on relation between rock post-peak dilatancy and confining pressure[J].Nonferrous Metals，2009，61(4):134－137.

[11]金濟(jì)山.巖石擴(kuò)容性質(zhì)及其本構(gòu)模型的研究[J].巖石力學(xué)與工程學(xué)報(bào)，1993，12(2):162－172.JIN Ji-shan.Study of the dilatancy behavior of rocks and its constitutive model[J].Chinese Journal of Rock Mechanics and Engineering，1993，12(2):162－172.

[12]侯文詩(shī)，李守定，李曉，等.巖石擴(kuò)容起始特性與峰值特性的比較[J].巖土工程學(xué)報(bào)，2010，32(5):657－660.HOU Wen-shi1，LI Shou-ding，LI Xiao，et al.Comparison between the dilatancy onset and the peak of different rocks[J].Chinese Journal of Geotechnical Engineering，2010，32(5):657－660.

[13]段艷燕，宋宏偉.巖石峰后剪脹效應(yīng)研究綜述[J].地下空間與工程學(xué)報(bào)，2007，3(6):1027－1030.DUAN Yan-yan，SONG Hong-wei.Review of dilatancy nature of the rock post-failure behaviour[J].Chinese Journal of Underground Space and Engineering，2007，3(6):1027－1030.

[14]閻巖，王恩志，王思敬，等.巖石滲流-流變耦合的試驗(yàn)研究[J].巖土力學(xué)，2010，31(7):2095－2103.YAN Yan，WANG En-zhi，WANG Si-jing，et al.Study of seepage-rheology coupling experiment of rocks[J].Rock and Soil Mechanics，2010，31(7):2095－2103.

[15]張玉軍，張維慶.孔隙率的應(yīng)力–剪脹依存性對(duì)CO2地層注入熱-氣-應(yīng)力耦合影響的數(shù)值分析[J].巖石力學(xué)與工程學(xué)報(bào)，2010，29(9):1774－1781.ZHANG Yu-jun，ZHANG Wei-qing.Numerical analysis of influence of stress-dilatancy dependence of porosity on thermo-gas-mechanical coupling in CO2injection in stratum[J].Chinese Journal of Rock Mechanics and Engineering，2010，29(9):1774－1781.

[16]張玉軍.遍有節(jié)理巖體的雙重孔隙-裂隙介質(zhì)熱-水-應(yīng)力耦合模型及有限元分析[J].巖石力學(xué)與工程學(xué)報(bào)，2009，28(5):947－955.ZHANG Yu-jun.Coupled thermo-hydro-mechanical model and finite elementanalyses of dual-porosity fractured medium for ubiquitous-joint rock mass[J].Chinese Journal of Rock Mechanics and Engineering，2009，28(5):947－955.

[17]OWEN D，HINTON E.Finite element in plasticity:Theory and practice[M].Swansea U.K.:Pineridge Press Limited，1980，

[18]馬士進(jìn).軟巖巷道圍巖擴(kuò)容軟化變形分析及模擬計(jì)算[D].遼寧:遼寧工程技術(shù)大學(xué)，2001.

[19]RUTQVIST J，B?RGESSON L，CHIJIMATSU M，et al.Thermo-hydro-mechanics of partially saturated geological media:governing equations and formulation of four finite element models[J].International Journal of Rock Mechanics and Mining Sciences，2001，38(1):105－127.

[20]TARON J，ELSWORTHD.Constraints on compaction rate and equilibrium in the pressure solution creep of quartz aggregates and fractures:Controls of aqueous concentration[J].Journal of Geophysical Research，2010，Vol，115:B07211.doi:10.1029/2009JB007118.

[21]LEE D，ELSWORTH D，YASUHARA H，et al.Experiment and modeling to evaluate the effects of proppant-pack diagenesis on fracture treatments[J].Journal of Petroleum Science and Engineering，2006，74:67－76.

[22]張玉軍.廢料地質(zhì)處置近場(chǎng)熱-水-應(yīng)力-遷移耦合二維有限元分析[J].巖土工程學(xué)報(bào)，2007，29(10):1553－1557.ZHANG Yu-jun.2D FEM analysis for coupled thermohydro-mechanical-migratory process in near field of geological disposal of nuclear waste[J].Chinese Jounal of Geotechnical Engineering，2007，29(10):1553－1557.

[23]CHIJIMATSU M，KURIKAMI H，ITO A，et al.Implication of T-H-M coupling on the near-field of a nuclear waste repository in a homogeneous rock mass[R].Tokyo:Hazama Corporation，2002：1－43.

[24]RUTQVIST J，CHIJIMATSU M，JING L，et al.A numerical study of T-H-M effects on the near-field safety of a hypothetical nuclear waste repository—BMT1 of the DECOVALEX III project.part 3.Effects of T-H-M coupling in sparsely fractured rocks[J].International Journal of Rock Mechanics and Mining Sciences，2005，42(5/6):745－755.