應力腐蝕和壓力溶解作用下雙重孔隙巖體中THMM耦合過程的有限元分析

張玉軍，楊朝帥

（1. 中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071；2. 中鐵隧道集團有限公司 技術中心，河南 洛陽 471009）

1 引 言

當具有流體通道的孔隙和裂隙的雙重巖體位于地下數(shù)千米的深度時，若在此處進行礦物和能源的回采、高放廢物及CO2的隔離等人工活動，則環(huán)境應力和溫度條件的改變可能影響這些通道的滲透特性。其綜合效應對裂隙而言分別稱之為應力腐蝕[1－4]和壓力溶解[5－8]，前者引起支撐糙度的破裂（或壓碎）、后者導致接觸糙度溶解，從而使得裂隙產生閉合和其滲透性的退化。

這里的應力腐蝕，是指作用在裂隙支撐糙度上的壓縮荷載引起局部拉應力時，將出現(xiàn)“次臨界”或“準靜態(tài)”的破裂,形成隨時間發(fā)展的漸進破壞，特別是有水存在時由化學反應促使次臨界裂紋的生長。而壓力溶解則包括了3個連續(xù)的過程：首先在承受應力的接觸糙度處礦物溶解；之后溶質沿著水膜擴散；最終在孔隙壁處溶質沉淀下來。

Dove[9]曾在較寬廣的溫度和 pH值范圍內進行室內試驗，以此嚴格地考察石英的溶解動力學，并得到了因化學溶解引起的I型裂紋擴展速率經(jīng)驗表達式。Yasuhara等[10]基于試驗數(shù)據(jù),考察一個含有單條天然裂隙的均密石英巖樣中的裂隙開度的演化，由此提出描述該現(xiàn)象的應力腐蝕和壓力溶解模型。Taron等[11]將壓力溶解、熱-水-應力收縮與膨脹、礦物沉淀與溶解對裂隙開度張開與閉合的復雜影響進行數(shù)學簡化，建立了一種新的熱-水-應力-化學雙重介質模型，并以此結合TOUGHREACT和FLAC3D分析了天然裂隙巖體的滲透率變化機制。筆者等也曾對其所建立的雙重孔隙—裂隙介質熱—水—應力耦合有限元程序作了改進，即引入Ki-Bok Min等[12]提出的簡化壓力溶解模型，對裂隙開度進行適時修正，建立裂隙的滲透系數(shù)隨壓力溶解的演化模式，并以一個假定的位于非飽和地層中的高放廢物地質處置庫為算例，針對3種裂隙開度的工況開展了數(shù)值模擬[13]。然而這其中并沒有考慮溶質濃度場。眾所周知，高放廢物地質處置中核素可能的泄漏與擴散是非常重要的研究內容[14]，故基于上述工作，完善已有的計算模型和有限元程序，進行熱-水-應力-遷移耦合分析勢在必行。

為此筆者首先在文獻[13]的控制方程中引入Yasuhara的應力腐蝕和壓力溶解模型，再加進溶質濃度場，即對于孔隙-裂隙雙重介質，認為應力場和溫度場是單一的，但存在著孔隙水壓力和裂隙水壓力、孔隙濃度和裂隙濃度，從而可模擬相應的熱-水-應力-遷移耦合過程。然后針對一個假定位于非飽和雙重孔隙-裂隙介質中的高放廢物地質處置庫，擬定2種工況：（1）裂隙開度隨應力腐蝕和壓力溶解而變化（基巖的孔隙率亦是應力的函數(shù)）；（2）裂隙開度和基巖的孔隙率均為常數(shù)，在一定的初始溫度、孔隙水壓力、巖體應力和核素釋放強度條件下，進行有限元分析，考察了處置庫近場的溫度、孔隙水壓力、水流速、飽和度、核素濃度和應力的分布與變化，得出了若干有意義的認識。

2 裂隙滲透特性的修正

2.1 應力腐蝕對裂隙開度的影響

如圖1所示，假定脆性材料的裂隙中粗糙面接觸是Hertzian接觸，則在這種接觸的內部或外部可能因拉應力σt產生圓形裂紋，此被稱為應力腐蝕。Dove[9]對石英材料定義其中I型裂紋的擴展速率為

式中：r為沿由σt產生的I型裂紋長軸方向的距離，假定裂紋的初始長度及r都是足夠小的；μ為材料的泊松比；σt為由σa引起的拉應力，其在裂隙的接觸邊緣達到最大值；σa為作用在裂隙接觸面積上的真實壓應力；為裂隙的名義面積（可取單位值）；Rc是裂隙的接觸面積，并且

Rc由下式求算：

式中：E、Er和E0依次為裂隙力學開度的平均值、殘余值和初始值；Rc0為裂隙接觸面積的初始值；a為經(jīng)驗常數(shù)。

從而裂隙力學開度因應力腐蝕的演化律為

圖1 微裂紋傳播引起的裂隙壓縮[10]Fig.1 Schematic of fracture compaction induced by miscrocrack propagation[10]

2.2 壓力溶解對裂隙開度的影響

Yasuhara等[10]定義了物質的溶解速率為

式中：d Mdiss/dt為從裂隙界面上進入溶液中的溶解質量速率；Vm為固體的摩爾體積；σc為臨界應力，當時體系達到平衡、反應停止；k+為固體的溶解速率常數(shù)；ρg為固體的密度；dc為粗糙面接觸的直徑。

則裂隙力學開度的壓力溶解閉合速率為

從而裂隙力學開度因壓力溶解的演化律亦為

2.3 裂隙滲透系數(shù)

設巖體中裂隙的間距為s，而單條裂隙開度為

故單條裂隙的水力開度為[15]

式中：JRC為裂隙粗糙度系數(shù)。

從而，等效后裂隙巖體的滲透系數(shù)為

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

2.4 應力對基巖滲透系數(shù)的影響

而對于基巖的孔隙率和滲透系數(shù)，當巖體中的應力發(fā)生改變時，可根據(jù)Davis等[16]提出經(jīng)驗公式進行如下修正

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

3 滲透-遷移方程引入

筆者對于參考文獻[17]中滲透-遷移方程進行了改進,其新意在于增加了孔隙巖體和裂隙介質之間因濃度差產生的溶質交換，其形式如下：

而擴散張量可表示為

式中：αiT為橫向彌散度；αiL為縱向彌散度；為表觀流速的絕對值；αim為分子擴散系數(shù)；τi為曲折率；δαβ為克羅內克符號。

筆者將Yasuhara的應力腐蝕和壓力溶解模型以及式（18）引入已有的雙重孔隙-裂隙巖體熱-水-應力耦合控制方程[13]中,其有限元算法可見文獻[18]。

4 算 例

如圖2所示，在地下深處埋入一個圓柱狀核廢料玻璃固化體，其周圍是非飽和的雙重孔隙-裂隙石英巖體。將該模型近似簡化為平面應變問 題。取計算域尺寸水平向為4 m，垂直向為8 m，有800個單元，861個節(jié)點。從固化體邊緣向右的點號依次為 432、433、434、435、436。其邊界條件為：頂面位移自由，其上作用有分布荷載σv= 26.7 MPa；左、右側面的水平方向位移約束；底面的垂直方向位移約束；所有邊界的孔隙水壓力為 -4.59 MPa、裂隙水壓力-0.46 MPa及溫度20 °C固定。巖體中發(fā)育有水平及垂直兩組裂隙，熱-水-應力耦合的環(huán)境對裂隙的開度要產生應力腐蝕和壓力溶解作用。有關的計算參數(shù)見表1～3（表3中除Rc0、和JRC外，其余參數(shù)取自文獻[10]）。初始狀態(tài)時，巖體的溫度均為20 °C。核廢物以1 000 W的不變功率釋放熱量，時間經(jīng)歷了4 a。

孔隙介質和裂隙介質的水分特性曲線符合 Van Genuchten模型，即

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

比滲透率與飽和度的關系為

取孔隙介質及裂隙介質的溫度梯度水分擴散系數(shù)相同，即 Dt= 2.5× 10-10m2/s ℃

貯存罐為孔隙巖體中的源項，某種放射性核素的釋放強度為 Qc=1.44× 10-10mol?kg/m3?s-1。設與核素的滲透遷移有關且在計算中不變的參數(shù)為：孔隙曲折率τ1、τ2分別為0.4、0.8，縱向彌散度α1L、α2L分別為1.0 m和2.0 m，橫向彌散度αiT=αiL/10，分子擴散系數(shù)α1m、α2m分別為 1.0×10-9m2/s和2.0×10-9m2/s，分配系數(shù)Kd1、Kd2分別為8.0 ml/g和 5.3 ml/g，干密度ρd1、ρd2分別為 2.3 g/m3和2.1 g/m3，參數(shù)為100.0 m-2，放射性核素的衰減常數(shù)χ= ln2/Thalf，其中Thalf為半衰期，取為1 000 a。

針對前述的2種裂隙開度變化工況，計算了巖體中的溫度、位移、孔隙水壓力、核素濃度和應力等的變化及分布情況，其主要結果及分析如下。

圖2 計算模型Fig.2 Computation model

工況1和工況2兩種條件下計算域中的溫度變化基本相同。以工況1為例，432、433、434、435各點處的溫度隨時間的變化曲線見圖3。由圖看到，在開始的約0.1 a內緩沖層的溫度快速上升，之后增加減緩，到計算終了時432、433、434、435各點的溫度依次為 77.8 °C、62.0 °C、52.5 °C 和 45.7 °C。

圖3 工況1溫度-時間曲線Fig.3 Temperatures versus time at some nodes for case 1

圖4、5分別為玻璃固化體右緣中點的水平和垂直裂隙開度因應力腐蝕、壓力溶解產生的閉合速率隨時間的變化。從中看到，兩種因素產生的閉合速率均先隨時間增加，到達峰值再衰減，之后逐漸趨于不變，并且應力腐蝕引起的閉合速率要比壓力溶解引起的閉合速率高 6個數(shù)量級，這與文獻[10]中呈現(xiàn)的規(guī)律類似；同時水平裂隙開度比垂直裂隙開度閉合速率大，這是由于巖體中垂直應力大于水平應力的緣故。

圖4 玻璃固化體右緣中點因應力腐蝕引起的-時間曲線Fig.4 caused by stress corrosion versus time at middle point of right margin of vitrified waste

圖5 玻璃固化體右緣中點因壓力溶解引起的-時間曲線Fig.5 caused by pressure solution versus time at middle point of right margin of vitrified waste

圖6、7分別表現(xiàn)了該中點處水平和垂直裂隙開度、粗糙面接觸率隨時間的演化，前者由初始開度下降而趨于殘余開度，后者由初始接觸率上升而趨于其名義接觸率（單位值），并且對應于水平裂隙的值的變化明顯要快。

圖6 玻璃固化體右緣中點裂隙開度-時間曲線Fig.6 Fracture apertures versus time at middle point of right margin of vitrified waste

圖7 玻璃固化體右緣中點粗糙面接觸率-時間曲線Fig.7 Contact-area ratios versus time at middle point of right margin of vitrified waste

從圖8看到，在該中點處水平裂紋上的應力強度因子比垂直裂紋上的應力強度因子大得多，且二者均隨時間減小。

圖8 玻璃固化體右緣中點應力強度因子-時間曲線Fig.8 Stress intensity factors versus time at middle point of right margin of vitrified waste

而圖9顯示了該中點處水平和垂直裂隙的臨界應力相同，且其在開始階段快速下降，之后逐漸趨于常數(shù)。產生上述現(xiàn)象的原因正是由溫度、應力與化學的綜合作用。

圖9 玻璃固化體右緣中點臨界應力-時間曲線Fig.9 Critical stresses versus time at middle point of right margin of vitrified waste

圖10、11為工況1、2中432、433、434、435各點處的孔（裂）隙水壓力隨時間的變化曲線。比較后看到，由于工況1考慮了應力腐蝕和壓力溶解對裂隙開度的減小作用及孔隙的滲透系數(shù)隨應力的變化，而工況2的裂隙開度及孔隙的滲透系數(shù)保持為常數(shù)，故工況1的負孔（裂）隙水壓力上升較高，特別是在應力腐蝕和壓力溶解效應最強烈的 432點，該處的負裂隙水壓力達到很大的值。在4 a時432點的負孔（裂）隙水壓力分別為工況1：-12.25、-7.95 MPa；工況 2：-6.03、-0.66 MPa。

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

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

圖12、13是工況1、2在4 a時貯存罐周圍2 m×2 m范圍內孔（裂）隙水壓力等值線的分布。從中看到，由于工況1受到應力腐蝕和壓力溶解的影響，與工況2相比，其中的裂隙水壓力在貯存罐周圍呈現(xiàn)明顯的增高區(qū)域。

圖14為兩種工況在4 a時計算域中的孔隙水流速和裂隙水流速矢量分布?？梢钥吹?，由于工況 1考慮了應力腐蝕和壓力溶解效應，與工況2相比，其中的裂隙水流速矢量有顯著的不同，特別是在貯存罐附近。以432點為例，孔隙水流速和裂隙水流速分別為工況1： 3.40 × 10-8、1 .52× 10-8m/s ；工況2： 2.32× 10-8， 2.77× 10-8m/s 。

圖15、16是2種工況中432、433、434、435各點處的孔隙核素濃度和裂隙核素濃度隨時間的變化曲線。比較后可知，由于工況1中考慮了應力腐蝕和壓力溶解對裂隙開度的減小作用、平均有效應力對孔隙率的壓縮功能，使得裂隙和孔隙的滲透系數(shù)降低，而工況2中的裂隙開度、孔隙率及其滲透系數(shù)均為常數(shù)，故工況1的裂隙和孔隙中核素聚集程度較高，到4 a時2種工況中432、433、434、435各點的核素濃度分別為孔隙巖體：20.18/6.82、15.42/3.60、12.06/2.55、9.36/1.76；裂隙介質：0.86/18.44、8.39/6.82、6.28/5.54、5.00/4.63（“/”的上、下分別是工況1、2的數(shù)據(jù)，其單位是 10-3mol/m3）。

圖12 工況1中4 a時2 m×2 m區(qū)域中的孔(裂)隙水壓力等值線(MPa)Fig.12 Contours of pore and fracture pressures in a 2 m×2 m area at 4 years for case 1 (MPa)

圖13 工況2中4 a時2 m×2 m區(qū)域中的孔(裂)隙水壓力等值線(MPa)Fig.13 Contours of pore and fracture pressures in a 2 m×2 m area at 4 years for case 2 (MPa)

圖14 4 a時計算域中的孔(裂)隙水流速矢量Fig.14 Flow vectors of pore and fracture water in calculation domain at 4 years

圖15 工況1中4 a 時巖體中核素濃度-時間曲線Fig.15 Nuclide concentrations versus time at some nodes for case 1

圖16 工況2中4 a時巖體中核素濃度-時間曲線Fig.16 Nuclide concentrations versus time at some nodes for case 2

圖17、18是2種工況4 a時貯存罐周圍2 m×2 m范圍內的核素濃度和裂隙核素濃度等值線分布。

圖17 工況1中4 a時2 m×2 m區(qū)域中核素濃度等值線(10-3 mol/m3)Fig.17 Contours of nuclide concentration in a 2 m×2 m area at 4 years for case 1 (10-3 mol/m3)

圖18 工況2中4 a時2 m×2 m區(qū)域中核素濃度等值線(10-3 mol/m3)Fig.18 Contours of nuclide concentration in a 2m×2m area at 4 years for case 2 (10-3 mol/m3)

由于不考慮負的孔（裂）隙水壓力對應力平衡的影響[19]，兩種工況的計算域中應力量值及分布差別很小，這里以工況1為例給出4 a時計算域中的正應力等值線，見圖19。從中看到，由于玻璃固化體的存在及放熱效應，使得計算域中的應力場不同于僅由巖體自重產生的應力場（后者的應力等值線為水平線）。到4 a時玻璃固化體右緣中點的水平應力和垂直應力分別為-0.124 MPa和-26.75 MPa。

圖19 4 a時工況1計算域中正應力等值線(MPa)Fig.19 Normal stress contours in calculation domain at 4 years for case 1 (MPa)

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

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

表3 應力腐蝕和壓力溶解的計算參數(shù)Table 3 Parameters for stress corrosion and pressure solution

5 結 論

針對已建立的雙重孔隙-裂隙介質熱-水-應力耦合控制方程，在應用了Yasuhara等提出的裂隙開度的應力腐蝕和壓力溶解算式的基礎上，引入溶質濃度場，發(fā)展為相應的熱-水-應力-遷移耦合模型。以一個假定的位于非飽和雙重孔隙-裂隙巖體中的高放廢物地質處置庫的核素泄漏為算例，就主要考慮裂隙開度是否隨應力腐蝕和壓力溶解而變化（同時基巖的孔隙率也作為是否是應力函數(shù)的處理）的2種工況，通過熱-水-應力-遷移耦合的二維有限元模擬，考察了巖體中的溫度、裂隙開度的閉合速率、閉合量、孔隙水壓力、地下水流速、核素濃度和正應力的變化、分布情況。計算結果表明：工況1、2中的溫度狀態(tài)差別不大,計算終了（4 a）時，近場的溫度可達到 30.0 °C～80.0 °C；應力腐蝕引起的閉合速率要高于壓力溶解引起的閉合速率6個數(shù)量級，且兩種因素產生的閉合速率隨時間先增加后減小，并趨于穩(wěn)定；隨著時間的推移裂隙開度由初始值下降并趨于殘余值，而粗糙面接觸率亦由初始值上升并趨于其名義值；作用在裂紋上的應力強度因子及臨界應力均隨時間減小，并趨于常數(shù)；當考慮應力腐蝕和壓力溶解時，與不計該作用時相比，近場的負裂隙水壓力上升很高，且其流速矢量場有顯著的不同；由于工況1中應力腐蝕、壓力溶解和平均有效應力使得裂隙和孔隙的滲透系數(shù)降低，而工況 2中的裂隙及孔隙的滲透系數(shù)均為常數(shù)，故前者的裂隙和孔隙中核素濃度較后者為高；但2種工況的巖體中的應力量值及分布基本相同。

[1]ATKINSON K. A fracture mechanics study of subcritical tensile cracking of quartz in wet environments[J]. Pure and Appl. Phys., 1979, 117: 1011－1024.

[2]DOVE M. The dissolution kinetics of quartz in sodium chloride solutions at 25° to 300°C[J]. Am. J. Sci., 1994,294: 665－712.

[3]CHESTER S, LENZ C, CHESTER M, et al. Mechanisms of compaction of quartz sand at diagenetic conditions[J].Earth Planet Sci. Lett., 2004, 220(3－4): 435－451.

[4]NARA Y, KANEKO K. Study of subcritical crack growth in andesite using the double torsion test[J]. Int. J. Rock.Mech. Min. Sci., 2005,42(4): 521－530.

[5]ROBIN F. Pressure solution at grain to grain contacts[J].Geochim. Cosmochim. Acta, 1978, 42: 1383－1389.

[6]DEWERS T, HAJASH A. Rate laws for water-assisted compaction and stress-induced water-rock interaction in sandstones[J]. J. Geophys. Res., 1995, 100: 13093－13112.

[7]REVIL A. Pervasive pressure-solution transfer: A poro-visco-plastic model[J]. Geophys. Res. Let., 1999,26: 255－258.

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

[9]DOVE M. Geochemical controls on the kinetics of quartz fracture at subcritical tensile stresses[J]. J. Geophys. Res.,1995, 100(B11): 22349–22359.

[10]YASUHARA H, ELSWORTH D. Compaction of a rock fracture moderated by competing roles of stress corrosion and pressure solution[J]. Pure. Appl. Geophys., 2008,165: 1289－1306.

[11]TARON J, ELSWORTH D. Thermal-hydrologicmechanical-chemical processes in the evolution of engineered geothermal reservoirs[J]. Int. J. Rock. Mech.Min. Sci., International Journal of Rock Mechanics and Mining Sciences, 2009, 01: j.ijrmms. doi: 07. 1016/200901.

[12]MIN K B, RUTQVIST J, ELSWORTH D. Chemical and mechanically mediated influences on the transport and mechanical characteristics of rock fractures[J].International Journal of Rock Mechanics and Mining Sciences, 2009, 02: j.ijrmms. doi: 10. 1016/200804002.

[13]張玉軍, 張維慶. 裂隙開度的壓力溶解對雙重孔隙介質熱-水-應力耦合影響的有限元分析[J]. 巖土力學,2010, 31(41): 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(41): 1269－1275.

[14]SONNENTHAL E, ITO A, SPYCHER N, et al.Approaches to modeling coupled thermal, hydrological,and chemical porcesses in the drift scale heater test at Yucca Mountain[J]. International Journal of Rock Mechanics and Mining Sciences, 2005. 42(5/6), 698－719.

[15]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.

[16]DAVIS J, DAVIS D. Stress-dependent permeability:Characterization and modeling[R]. [S.l.]: Society of Petroleum Engineers, 1999.

[17]NISHIGAKI M. Density dependent transport analysis saturated-unsaturated porous media—3 dimensional eulerian lagrangian method[R]. [S.l.]: Okayama University, 2001.

[18]張玉軍, 張維慶. 雙重孔隙介質 THMM 耦合模型及其有限元分析[J]. 力學學報, 2010, 42(4): 660－669.ZHANG Yu-jun, ZHANG Wei-qing. Coupled thermohydro-mechanical-migratory model and FEM analyses for dual-porosity medium[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(4): 660－669.

[19]CHIJIMATSU M, KURIKAMI H, ITO A, et al.Implication of THM coupling on the near-field of a nuclear waste repository in a homogeneous rock mass[M].[S.l.]: [s.n.], 2002.