張玉軍，徐 剛，楊朝帥

（中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071）

1 引 言

在地層深部，由于地應(yīng)力、溫度及水頭均較高，使得巖體的物理力學性質(zhì)與淺埋時相比有顯著的不同[1]。對巖體裂隙而言，其變形剛度、開度及滲透性都要受到熱、水、應(yīng)力、化學等廣義力源場的強烈作用，從而與這些因素呈現(xiàn)出明顯的非線性關(guān)系。對此國內(nèi)外學者已進行了大量研究。如Rutqvist等[2－3]提出了立方-塊體概念模型，將裂隙的開度表示為法向應(yīng)力的指數(shù)函數(shù)，用其對裂隙的滲透系數(shù)進行修正；Min等[4]基于室內(nèi)試驗結(jié)果，認為巖體裂隙滲透率的改變起因于多種因素，從而提出了加載（閉合）和卸載（張開）時單條裂隙開度變化的數(shù)學模型；Olivella等[5]在裂隙孔隙率與巖體體積應(yīng)變（其取決于平均應(yīng)力和溫度）之間建立起關(guān)系，以此給出了裂隙的滲透系數(shù)-孔隙率經(jīng)驗公式；Yasuhara等[6－8]考慮應(yīng)力腐蝕、壓力溶解和自由面溶解、沉淀對裂隙擴張與閉合的作用，建立了相應(yīng)的裂隙開度隨時間演化的表達式；劉才華等[9]通過對人工充填砂裂隙實施剪切試驗，提出了剪應(yīng)力和法向有效應(yīng)力耦合作用下裂隙巖石的滲流模型；陳益峰等[10]通過將巖體單裂隙視為非關(guān)聯(lián)理想彈塑性體，導出單裂隙在壓剪荷載作用下，其機械開度和水力傳導度的解析模型。而在描述裂隙的力學性質(zhì)、應(yīng)力與位移的關(guān)系方面，現(xiàn)在有著名的常剛度的 Goodman模型[11]、變剛度的連續(xù)屈服模型[12]和Barton-Bandis模型[13]。

在國內(nèi)外已有的針對增強地熱系統(tǒng)（EGS）、CO2地質(zhì)封存、高放核廢物地質(zhì)處置等所涉及的熱-水-應(yīng)力（THM）耦合分析中，一般計入了上述巖體裂隙的開度及滲透性隨廣義力源的演化[3，4，6]，如筆者對其所建立的遍有節(jié)理巖體雙重孔隙-裂隙介質(zhì)熱-水-應(yīng)力耦合模型作了改進，即利用 Rutqvist等的立方-塊體概念模型、Taron的單條裂隙開度變化模型[14]以及 Davis等的基巖孔隙率和滲透系數(shù)演化模型[15]，根據(jù)即時的溫度和應(yīng)力，對裂隙開度進行熱-力張閉、壓力（化學）溶解的調(diào)整，從而不斷地對孔隙率和滲透系數(shù)實施修正，并研制出相應(yīng)的二維有限元程序[16]，以一個假定的高放廢物地質(zhì)處置庫為對象進行了數(shù)值模擬。但這些計算中均是將裂隙的剛度作為不變量。前已講過，在巖體埋深較大的條件下，裂隙的剛度應(yīng)隨應(yīng)力呈非線性變化，這不僅對地應(yīng)力場有影響，而且因力、水互動也將對滲流場產(chǎn)生作用。因此，在相應(yīng)的THM耦合分析中考慮裂隙的剛度隨應(yīng)力等因素的變化應(yīng)更接近實際情況。

在本文中，筆者基于文獻[16]的工作，引入并修正了變剛度的連續(xù)屈服節(jié)理模型，以一個假定的高放廢物地質(zhì)處置庫為算例，在相同的初始溫度、孔隙水壓力和巖體應(yīng)力條件下，設(shè)定裂隙的開度及滲透系數(shù)、基巖的孔隙率和滲透系數(shù)分別遵從Taron、Rutqvist和Davis等提出模型規(guī)律變化，針對以下2種工況：①裂隙剛度是法向應(yīng)力的指數(shù)函數(shù)（限定變化范圍）；②裂隙剛度是常數(shù)，進行熱-水-應(yīng)力耦合平面有限元計算，考察了處置庫近場的溫度、孔隙水壓力、裂隙水壓力、水流速度、裂隙剛度和正應(yīng)力的分布與變化，得出了若干有意義的認識。

2 裂隙剛度和開度及滲透系數(shù)

筆者所建立的遍有節(jié)理巖體雙重孔隙-裂隙介質(zhì)熱-水-應(yīng)力耦合模型及其有限元方法可見文獻[16]?，F(xiàn)在此基礎(chǔ)上對裂隙剛度作如下修正。

Cundall等[12]在1984年曾提出了一種連續(xù)屈服的節(jié)理模型（continuously yielding joint model），其變形特性與荷載的關(guān)系可表示為

式中：kn、ks分別為法向剛度和切向剛度；σn為法向應(yīng)力；an、as、en和es為模型參數(shù)。

筆者在引入該模型的同時，為了體現(xiàn)0應(yīng)力狀態(tài)時的節(jié)理具有初始剛度、以及拉應(yīng)力對剛度值的減弱作用，將式（1）修改為

式中：kn0、ks0分別為初始法向剛度和切向剛度；α分別取“+1”及“-1”，對應(yīng)于壓應(yīng)力和拉應(yīng)力。

為了避免應(yīng)力絕對值較高時產(chǎn)生節(jié)理剛度過大、過小甚至出現(xiàn)負值從而使計算不穩(wěn)定的問題，筆者亦對剛度的變化限定了最大量和最小量[13]，即kn,max、ks,max和kn,min、ks,min。

再引入裂隙的開度與其粗糙度系數(shù)的關(guān)系式[17]

式中：JRC為裂隙粗糙度系數(shù)；B和e分別為巖體內(nèi)裂隙的力學開度和水力開度。

從而，得等效裂隙介質(zhì)（均化后）的滲透系數(shù)為

式中：g為重力加速度，g =9.81 m/s2；ν為運動黏度（對于20 ℃的純水，ν= 1.0×10-6m2/s ）；s為裂隙間距。

對于二維問題，在xy平面內(nèi)巖體裂隙的孔隙率及滲透系數(shù)可表示為[2]

分別為巖體裂隙孔隙率及滲透系數(shù)的修正因子。

并且有

式中：Bx、By分別為初始狀態(tài)時沿坐標軸方向的裂隙開度；Bxt、Byt分別為t時刻的相應(yīng)裂隙開度。

而加載（閉合）時單條裂隙開度變化的計算式為[14]

式中：Bmc為在給定溫度T、有效應(yīng)力σ′時的裂隙開度；分別為極端力學應(yīng)力、極端化學應(yīng)力時的殘余裂隙開度；為最大可能的裂隙開度改變量（或初始與殘余的裂隙開度之差）；ω、β、χ為經(jīng)驗參數(shù)，依次表征裂隙開度的應(yīng)力、化學和溫度依存性。

當卸載（張開）時裂隙開度的計算式為

式中：Bmc(σ′max)為最大荷載時的裂隙開度，根據(jù)卸載前的最大應(yīng)力由式（9）確定；Rm為裂隙開度的恢復率，是通過熱-力卸載恢復的力學（非化學）閉合量。

將式（9）、（10）中的σ′分別用t時刻的σ′x、σ′y取代可得到相應(yīng)的Bxt和Byt。

而基巖的孔隙率和滲透系數(shù)的計算式為[15]

式中：φ0、k0分別為0應(yīng)力狀態(tài)的基巖孔隙率和滲透系數(shù)；φr為高應(yīng)力狀態(tài)下的殘余基巖孔隙率；σ′m為平均有效應(yīng)力；a、c分別為由試驗確定的參數(shù)。

3 算例模型

如圖1所示，有一個圓柱狀核廢料玻璃固化體被埋入地下一定深處，其周圍的巖體為非飽和的雙重孔隙-裂隙介質(zhì)。作為近似簡化，可認為這是一個平面應(yīng)變問題。取計算域尺寸水平向為4 m，垂直向為8 m，有800個單元，861個節(jié)點。從固化體邊緣向右的點號依次為432、433、434、435、436，彼此間隔0.2 m。

圖1 計算模型Fig.1 Computation model

對于邊界條件，計算域的頂面位移自由，其上作用有分布荷載 σv= 26.7 MPa ；左、右側(cè)面的水平方向位移約束；底面的垂直方向位移約束；所有邊界的 孔 隙 水 壓 力 Pw1=-4 .59 MPa、 裂 隙 水 壓 力Pw2=- 0.46 MPa及溫度T=20 °C固定。巖體中發(fā)育有水平及垂直兩組裂隙。有關(guān)的計算參數(shù)見表1、2。初始狀態(tài)時，巖體的溫度為20 °C。核廢物以1 kW的不變功率釋放熱量，時間經(jīng)歷了4 a。

孔隙介質(zhì)和裂隙介質(zhì)的初始飽和度分別為0.44和0.02，其水分特性曲線符合Van Genuchten模型，即

式中：對于孔隙介質(zhì)：α = 3.86× 10-6m-1，n=1.41；對于裂隙介質(zhì) ： α= 5.26× 10-4m-1，n=2.55；m=1- 1/n；ψ為水勢；sws為最大飽和度，其值為1.0；swr為最小飽和度，其值對孔隙介質(zhì)和裂隙介質(zhì)各取0.19和0.01。

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

取孔隙介質(zhì)及裂隙介質(zhì)的溫度梯度水分擴散系數(shù)相同，其為

針對巖體裂隙剛度變化2種工況，計算了巖體中的溫度、孔隙水壓力、裂隙水壓力、飽和度、水流速度和正應(yīng)力的分布與變化情況。

4 結(jié)果分析

2種工況條件下計算域中的溫度變化有一定差別。以工況 1為例，圖 1中 432、433、434、435各點處的溫度隨時間的變化曲線見圖 2?？吹皆陂_始的約1 a內(nèi)緩沖層的溫度快速上升，之后增加減緩,到計算終了時 432、433、434、435各點的溫度依次為80.2、64.3、54.7、47.6 °C，而工況2的對應(yīng)點溫度為 85.6、69.4、59.5、52.0 °C，產(chǎn)生差別的原因在于：在工況1中，因為裂隙的剛度隨法向應(yīng)力而變化（壓應(yīng)力時剛度增加），使得巖體的整體變形較小，如計算域頂部中點的垂直向下位移位為0.25 m（工況2中達到0.72 m），則由于變形-熱耦合效應(yīng)導致巖體中溫度較低。但2種工況中溫度分布相似，圖3是工況1在4 a時計算域中的溫度等值線。

圖2 溫度-時間曲線Fig.2 Temperatures versus time at some nodes

圖3 4 a時的溫度等值線 (單位：°C)Fig.3 Temperatures contours at 4 years (unit: °C)

圖4是工況1在4 a時計算域中2組裂隙的剛度等值線，其變化范圍：水平裂隙的kn=2.94～2.95 GPa/m、ks=1.47～1.48 GPa/m，垂直裂隙的kn=1.65～1.80 GPa/m、ks=0.82～0.90 GPa/m。而工況2的2組裂隙剛度在全域內(nèi)為常數(shù)，即kn、ks分別為1.0、0.5 GPa/m。2種工況在4 a時計算域中正應(yīng)力等值線見圖5、6。從圖看到，由于工況1中裂隙剛度是法向應(yīng)力的函數(shù)，與常裂隙剛度的工況2相比，兩者的巖體應(yīng)力場有較顯著的差別，如玻璃固化體右緣中點的水平正應(yīng)力和垂直正應(yīng)力之比σH/ σV分別為，工況1為-0.17 MPa/-26.33 MPa，工況2為-0.09 MPa/ -26.66 MPa?？煽闯龉r1中裂隙剛度“場”與正應(yīng)力分布有明顯的相似性。

圖4 工況1在4 a時2組裂隙的剛度等值線（單位：GPa/m）Fig.4 Stiffness contours of two fracture sets at 4 years for case 1 (unit: GPa/m)

圖5 工況1, 4 a時巖體中正應(yīng)力等值線 (單位：MPa)Fig.5 Normal stress contours in rock mass at 4 years for case 1 (unit: MPa)

圖6 工況2，4 a時巖體中正應(yīng)力等值線 (單位：MPa)Fig.6 Normal stress contours in rock mass at 4 years for case 2 (unit: MPa)

圖7是工況1中432、433、434、435各點處的孔隙水壓力及裂隙水壓力隨時間的變化曲線。從圖看到，在開始的約0.03 a之內(nèi)，負裂隙水壓力迅速上升，并一度達到-3.22～-2.85 MPa，隨之由于來自地下水浸潤的水分補給大于溫度梯度引起的水分減少，使得孔隙和裂隙中的飽和度不斷增加，表現(xiàn)為負的孔（裂）隙水壓力絕對值較快下降，在約0.6 a之后孔（裂）隙水壓力僅有微小波動，二者漸趨穩(wěn)定。裂隙水壓力一開始出現(xiàn)較大波動的原因在于：初期時巖體應(yīng)力調(diào)整幅度大，致使較高的孔隙水壓力（絕對值）與較低的裂隙水壓力（絕對值）快速趨于平衡。因為加熱導致裂隙空間增大，而裂隙水體積變化相對較小引起的裂隙飽和度減少，進而導致負壓絕對值增大。到計算終了的4 a時，上述 4個測點處的孔隙水壓力/裂隙水壓力依次為-0.785 MPa/-0.735 MPa、-0.773 MPa/-0.728 MPa、-0.756 MPa/-0.719 MPa和-0.739 MPa/-0.706 MPa，此時巖體中孔（裂）隙水壓力等值線分布見圖 8。對于工況 2，巖體中孔（裂）隙水壓力隨時間的變化、分布及量值與工況1相似，其絕對值略大一點，如計算終了時上述 4個點處的孔隙水壓力/裂隙水壓力依次為-0.796 MPa /-0.743 MPa、-0.784 MPa/-0.737 MPa、-0.766 MPa /-0.727 MPa和-0.748 MPa/-0.714 MPa。

2種工況條件下計算域中的孔隙水流速和裂隙水流速矢量分布也相似。以工況1為例，4 a時巖體中孔（裂）隙水流速矢量分布見圖 9。注意兩種流速矢量的比例為1:0.1。其中以432點為例，孔隙水流速和裂隙水流速分別為：工況1：3 .00× 10-9m/s ，2.41× 10-8m/s ；工況 2： 3.10× 10-9m/s ，2.55×10-8m/s 。

圖7 工況1孔(裂)隙水壓力-時間曲線Fig.7 Pore and fracture water pressures versus time at some nodes for case 1

圖8 工況1，4 a時巖體中水壓力等值線 (單位：MPa)Fig.8 Contours of pore and fracture water pressures in rock mass at 4 years for case 1 (unit: MPa)

圖9 工況1在4 a時計算域中的水流速矢量Fig.9 Flow vectors of pore and fracture water in calculation domain for case 1 at 4 years

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

表2 裂隙組的計算參數(shù)Table 2 Parameters for fracture sets used in calculation

5 結(jié) 論

（1）對于所建立的遍有節(jié)理巖體雙重孔隙-裂隙介質(zhì)熱-水-應(yīng)力耦合模型，在可考慮應(yīng)力拉壓和壓力（化學）溶解對裂隙開度的綜合影響的基礎(chǔ)上，改進了裂隙剛度的計算方法（其為法向應(yīng)力的指數(shù)函數(shù)），從而實現(xiàn)了裂隙剛度、孔隙率和滲透系數(shù)的適時修正，使得相應(yīng)的二維有限元程序更能體現(xiàn)THM耦合現(xiàn)象的非線性特點。

（2）以一個假定的高放廢物地質(zhì)處置庫為算例，進行了巖體裂隙剛度變化的2種工況的數(shù)值分析，其結(jié)果顯示：與裂隙剛度不變的情況相比，裂隙剛度隨法向應(yīng)力變化時計算域中溫度場量值較低，巖體應(yīng)力的大小也有一定的不同，其分布與裂隙剛度“場”有明顯的相似性；2種工況的孔（裂）隙水滲流場形態(tài)相似，但裂隙剛度隨法向應(yīng)力變化時負孔（裂）隙水壓力的絕對值略小。這表明裂隙剛度改變會對THM耦合過程及其場量產(chǎn)生一定的影響。

（3）在筆者改進的裂隙剛度隨法向應(yīng)力變化的指數(shù)函數(shù)中，也可考慮拉應(yīng)力對裂隙剛度的減弱效應(yīng)，然而在本文計算的特定條件下，巖體中全是壓應(yīng)力，故其對裂隙剛度只有增大作用，這對 THM耦合的影響從總體上看是不太大的。筆者推測，巖體中若有范圍大、量值高的拉應(yīng)力存在，對 THM耦合的影響將可能更明顯。對此應(yīng)進一步的研究。

[1]U. S. National Committee on Tunneling Technology and U.S. National Committee for Rock Mechanics. Advances in technology for the construction of deep-underground facilities[J]. Tunnelling and Underground Space Technology, 1988, 3(1): 25－44.

[2]JONNY RUTQVIST, CHIN-FU TSANG. Analysis of thermal-hydrologic-mechanical behavior near an emplacement drift at Yucca mountain[J]. Journal of Contaminant Hydrology, 2003, 62－63: 637－652.

[3]RUTQVIST J, BARR D, DATTA R, et al. Coupled thermal-hydrological-mechanical analyses of the Yucca mountain drift scale test—comparison of field measurements to predictions of four different numerical models[J]. International Journal of Rock Mechanics and Mining Sciences, 2005, 42(5/6): 680－697.

[4]MIN K B, RUTQVIST J, ELSWORTH D. Chemically and mechanically mediated influences on the transport and mechanical characteristics of rock fractures[J].International Journal of Rock Mechanics and Mining Sciences, 2009, 46(1): 80－89.

[5]OLIVELLA S, GENS A. Double structure THM analyses of a heating test in a fractured tuff incorporating intrinsic permeability variations[J]. International Journal of Rock Mechanics and Mining Sciences, 2005, 42(5/6):667－679.

[6]YASUHARA H, ELSWORTH D. Compaction of a rock fracture moderated by competing roles of stress corrosion and pressure solution[J]. Pure Applied Geophysics, 2008,165(7): 1289－1306.

[7]YASUHARA H, ELSWORTH D, POLAK A. Evolution of permeability in a natural fracture: Significant role of pressure solution[J]. Journal of Geophysical Research,2004, 109(B3): B03204. doi:10.1029/2003JB002663.

[8]YASUHARA H, ELSWORTH D, POLAK A, et al.Spontaneous switching between permeability enhancement and degradation in fractures in carbonate:Lumped parameter representation of mechanically and chemically-mediated dissolution[J]. Transport in Porous Media, 2006, 65: 385－409.

[9]劉才華, 陳從新, 付少蘭. 二維應(yīng)力作用下巖石單裂隙滲流規(guī)律的實驗研究[J]. 巖石力學與工程學報, 2002,21(8): 1194－1198.LIU Cai-hua, CHEN Cong-xin, FU Shao-lan. Testing study of seepage characteristic of a single rock fracture under two-dimentional stresses[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(8): 1194－1198.

[10]陳益峰, 周創(chuàng)兵, 盛永清. 應(yīng)變敏感的裂隙及裂隙巖體水力傳導特性研究[J]. 巖石力學與工程學報, 2006,25(12): 2441－2452.CHEN Yi-feng, ZHOU Chuang-bing, SHENG Yong-qing.Strain-dependent hydraulic conductivity for single rock fracture and fractured rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(12): 2441－2452.

[11]周維垣. 高等巖石力學[M]. 北京: 水利電力出版社,1990: 276－282.

[12]CUNDALL P A, HART R D. Analysis of block test No. 1 inelastic rock mass behavior: Phase 2—A characterization of joint behavior (final report)[R]. Minneapolis: Itasca Consulting Group, 1984.

[13]BANDIS S C, BARTON N R, CHRISTIANSON M.Application of a new numerical model of joint behaviour to rock mechanics problems[C]//Fundamentals of Rock Joints (Proceedings of the International Symposium on Fundamentals of Rock Joints. Lulea, Sweden: Centek Publishers, 1985: 345－356.

[14]TARON J, ELSWORTH D. Thermal-hydrologicmechanical-chemical processes in the evolution of engineered geothermal reservoirs[J]. International Journal of Rock Mechanics and Mining Sciences, 2009,46(5): 855－864.

[15]DAVIS J P, DAVIS D K. Stress-dependent permeability:Characterization and modeling[R]. [S. l.]: Society of Petroleum Engineers. 1999, SPE Paper no. 56813.

[16]張玉軍, 張維慶. 裂隙開度的壓力溶解對雙重孔隙介質(zhì)熱-水-應(yīng)力耦合影響的有限元分析[J]. 巖土力學,2010, 31(4): 1269－1275.ZHANG Yu-jun, ZHANG Wei-qing. Finite element analysis of influence of pressure solution of fracture aperture on T-H-M coupling in dual-porosity medium[J].Rock and Soil Mechanics, 2010, 31(4): 1269－1275.

[17]OLSSON R, BARTON N. An improved model for hydromechanical coupling during shearing of rock joints[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(3): 317－329.