裂隙剛度隨應(yīng)力變化對(duì)雙重孔隙介質(zhì)熱-水-應(yīng)力耦合影響的有限元分析

張玉軍，徐 剛，楊朝帥

（中國科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點(diǎn)實(shí)驗(yàn)室，武漢 430071）

1 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

并且有

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

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

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

當(dāng)卸載（張開）時(shí)裂隙開度的計(jì)算式為

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

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

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

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

3 算例模型

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

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

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

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

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

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

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

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

4 結(jié)果分析

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

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

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

圖4是工況1在4 a時(shí)計(jì)算域中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時(shí)計(jì)算域中正應(yīng)力等值線見圖5、6。從圖看到，由于工況1中裂隙剛度是法向應(yīng)力的函數(shù)，與常裂隙剛度的工況2相比，兩者的巖體應(yīng)力場有較顯著的差別，如玻璃固化體右緣中點(diǎn)的水平正應(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時(shí)2組裂隙的剛度等值線（單位：GPa/m）Fig.4 Stiffness contours of two fracture sets at 4 years for case 1 (unit: GPa/m)

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

圖6 工況2，4 a時(shí)巖體中正應(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各點(diǎn)處的孔隙水壓力及裂隙水壓力隨時(shí)間的變化曲線。從圖看到，在開始的約0.03 a之內(nèi)，負(fù)裂隙水壓力迅速上升，并一度達(dá)到-3.22～-2.85 MPa，隨之由于來自地下水浸潤的水分補(bǔ)給大于溫度梯度引起的水分減少，使得孔隙和裂隙中的飽和度不斷增加，表現(xiàn)為負(fù)的孔（裂）隙水壓力絕對(duì)值較快下降，在約0.6 a之后孔（裂）隙水壓力僅有微小波動(dòng)，二者漸趨穩(wěn)定。裂隙水壓力一開始出現(xiàn)較大波動(dòng)的原因在于：初期時(shí)巖體應(yīng)力調(diào)整幅度大，致使較高的孔隙水壓力（絕對(duì)值）與較低的裂隙水壓力（絕對(duì)值）快速趨于平衡。因?yàn)榧訜釋?dǎo)致裂隙空間增大，而裂隙水體積變化相對(duì)較小引起的裂隙飽和度減少，進(jìn)而導(dǎo)致負(fù)壓絕對(duì)值增大。到計(jì)算終了的4 a時(shí)，上述 4個(gè)測(cè)點(diǎn)處的孔隙水壓力/裂隙水壓力依次為-0.785 MPa/-0.735 MPa、-0.773 MPa/-0.728 MPa、-0.756 MPa/-0.719 MPa和-0.739 MPa/-0.706 MPa，此時(shí)巖體中孔（裂）隙水壓力等值線分布見圖 8。對(duì)于工況 2，巖體中孔（裂）隙水壓力隨時(shí)間的變化、分布及量值與工況1相似，其絕對(duì)值略大一點(diǎn)，如計(jì)算終了時(shí)上述 4個(gè)點(diǎn)處的孔隙水壓力/裂隙水壓力依次為-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種工況條件下計(jì)算域中的孔隙水流速和裂隙水流速矢量分布也相似。以工況1為例，4 a時(shí)巖體中孔（裂）隙水流速矢量分布見圖 9。注意兩種流速矢量的比例為1:0.1。其中以432點(diǎn)為例，孔隙水流速和裂隙水流速分別為：工況1：3 .00× 10-9m/s ，2.41× 10-8m/s ；工況 2： 3.10× 10-9m/s ，2.55×10-8m/s 。

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

圖8 工況1，4 a時(shí)巖體中水壓力等值線 (單位：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時(shí)計(jì)算域中的水流速矢量Fig.9 Flow vectors of pore and fracture water in calculation domain for case 1 at 4 years

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

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

5 結(jié) 論

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

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

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

[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ī)律的實(shí)驗(yàn)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 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)變敏感的裂隙及裂隙巖體水力傳導(dǎo)特性研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 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]周維垣. 高等巖石力學(xué)[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]張玉軍, 張維慶. 裂隙開度的壓力溶解對(duì)雙重孔隙介質(zhì)熱-水-應(yīng)力耦合影響的有限元分析[J]. 巖土力學(xué),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.