張玉軍，楊朝帥,

(1. 中國科學(xué)院 武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢，430071)2. 中國中鐵隧道集團(tuán)有限公司 技術(shù)中心，河南 洛陽，471009)

當(dāng)?shù)責(zé)崽锘蚋叻派湫院藦U物處置庫位于地層深部的高溫、高壓及富水環(huán)境中時(shí)，對(duì)于孔隙?裂隙巖體而言，存在明顯的壓力溶解作用[1]，其主要使得巖體的孔隙率和滲透系數(shù)減少，從而對(duì)熱?水?應(yīng)力?化學(xué)耦合過程產(chǎn)生影響。對(duì)此，國外已經(jīng)進(jìn)行較多的研究[2?6]。其中，Yasuhara等[7]通過室內(nèi)試驗(yàn)，提出一種表現(xiàn)壓力溶解三過程的數(shù)學(xué)模型；Taron等[8]基于Yasuhara的工作，將壓力溶解、熱?水?應(yīng)力收縮與膨脹、礦物沉淀與溶解對(duì)裂隙開度張開與閉合的復(fù)雜影響進(jìn)行數(shù)學(xué)簡(jiǎn)化，并以此結(jié)合 TOUGHREACT和FLAC3D軟件分析天然裂隙巖體的滲透率變化機(jī)理；Taron等[9?10]還考慮擬靜水或有邊界流的情況，發(fā)展一種新的壓力溶解模型，并將其應(yīng)用于具有動(dòng)態(tài)滲透率的裂隙巖體工程地?zé)釒熘旭詈系牧W(xué)和化學(xué)過程的數(shù)值模擬。上述壓力溶解模型具有很高的學(xué)術(shù)價(jià)值和實(shí)用價(jià)值，但目前尚未見將其用于顆粒聚集體孔隙介質(zhì)的多場(chǎng)耦合數(shù)值分析中，而這種分析對(duì)于在諸如巖鹽、黏土類巖體中進(jìn)行高放廢物地質(zhì)處置等環(huán)境工程有著重要的指導(dǎo)價(jià)值。為此，本文作者建立壓力溶解模型，引入所研制的巖土介質(zhì)熱?水?應(yīng)力耦合有限元分析程序中(對(duì)其中的體積應(yīng)變求解式進(jìn)行一定修正)，給出顆粒聚集體孔隙率隨壓力溶解和滲透系數(shù)隨孔隙增縮的演化表達(dá)式，在計(jì)算時(shí)對(duì)孔隙率和及滲透系數(shù)實(shí)施動(dòng)態(tài)調(diào)整。以1個(gè)假定的位于飽和地層中的高放廢物地質(zhì)處置庫為算例，在相同的初始溫度和巖體應(yīng)力條件下，針對(duì)2種計(jì)算工況即邊界上水壓力相同(域內(nèi)無因邊界水頭差產(chǎn)生的穩(wěn)態(tài)流量，工況1)和左、右邊界上水壓力有較大差別(域內(nèi)有因邊界水頭差產(chǎn)生的穩(wěn)態(tài)流量，工況2)，進(jìn)行熱?水?應(yīng)力耦合有限元數(shù)值模擬，考察處置庫近場(chǎng)的溫度、顆粒界面水膜及孔隙中的溶質(zhì)濃度、遷移和沉淀質(zhì)量、孔隙率和滲透系數(shù)、孔隙水壓力、地下水流速、應(yīng)力的變化、分布情況。

1 混合質(zhì)量平衡

對(duì)于如圖1和圖2所示的顆粒聚集體化學(xué)壓縮模型，Taron等提出物質(zhì)在壓力作用下的溶解濃度平衡關(guān)系[10]：

圖1 孔隙介質(zhì)、裂隙及化學(xué)壓縮過程的概念化[10]Fig. 1 Conceptualization of chemical compaction process,porous media, and fractures[10]

圖2 壓縮及顆粒和生長(zhǎng)的膠結(jié)物半徑的關(guān)系示意[10]Fig. 2 Diagram of compaction and relationship between granular and growing cement radius[10]

其迭代求解式為

式中：ci是在顆粒接觸面水膜中溶解硅的即時(shí)濃度；cp是在顆粒周界(孔隙空間)處的即時(shí)硅濃度；k+是溶解速率常數(shù)；α為固體接觸部分所占據(jù)的顆粒邊界的比例(0≤α≤1)；為表征體顆粒交界處的反應(yīng)面積(REV)；為液態(tài)水中硅的溶解度；D為擴(kuò)散系數(shù)；Vi是顆粒接觸之內(nèi)表征體的體積(REV)；為表征體的孔隙反應(yīng)面積(REV)；Vp為表征體的孔隙體積(REV)；aSiO2為應(yīng)力作用下固體硅的化學(xué)活度；Qss為因邊界水頭差產(chǎn)生的穩(wěn)態(tài)流量；Δt為時(shí)間增量。

式(1)中各物理量可表示為：

式中：V0和V分別為表征體(可取單位量)的初始和即時(shí)體積；0φ和φ分別為初始和即時(shí)的孔隙率；Sag和Vg分別為平均的單個(gè)顆粒表面積和體積；ag為2個(gè)相互貫穿球形顆粒的接觸面積；Nc為每個(gè)顆粒的接觸數(shù)；Df為分子擴(kuò)散系數(shù)；ω為顆粒接觸面水膜的厚度；Δμ為化學(xué)勢(shì)梯度；R這氣體常數(shù)；T為絕對(duì)溫度。

并且有

式中：ri為初始時(shí)刻的平均顆粒半徑；h為累積的球形顆粒貫穿深度(初值為0)；σ′為有效應(yīng)力；Kb為顆粒材料的體積模量；為產(chǎn)生于顆粒接觸面處的遷移質(zhì)量(以溶解的數(shù)值為正)；Vm為固體的摩爾體積。

而

式中：σa和分別是真實(shí)應(yīng)力及達(dá)到平衡時(shí)的臨界應(yīng)力；Rc為接觸面積比；EA和TA分別為熔化熱和熔化溫度；βc為埋置常數(shù)，rgc為顆粒交界處沉淀膠結(jié)物外表面到顆粒中心連線的距離(如圖2所示)。

2 孔隙率?滲透系數(shù)演化律

Taron等給出的孔隙率表達(dá)式為

式中：εV為有限的體積應(yīng)變；R為作為源/匯項(xiàng)的隨時(shí)間求和的反應(yīng)體積。

并且

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

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

將上述顆粒聚集體的壓力溶解模型引入文獻(xiàn)[12]的孔隙介質(zhì)熱?水?應(yīng)力耦合有限元程序中。

3 算例

實(shí)驗(yàn)室尺度的核廢料玻璃固化體埋置模型如圖 3所示，其周圍的介質(zhì)是飽和的石英顆粒聚集孔隙巖體。作為近似簡(jiǎn)化，認(rèn)為這是一個(gè)平面應(yīng)變問題。取計(jì)算域?yàn)樗较? m，垂直向8 m，有800個(gè)單元，861個(gè)節(jié)點(diǎn)。從固化體邊緣向右的點(diǎn)號(hào)依次為432，433，434，435和436。

對(duì)于邊界條件，計(jì)算域的頂面位移自由，其上作用有分布荷載53.4 MPa；左、右側(cè)面的水平方向位移約束；底面的垂直方向位移約束；所有邊界的溫度為40 ℃；對(duì)于工況 1，各邊界的孔隙水壓力均為4.59 MPa；對(duì)于工況 2，左、右邊界的孔隙水壓力分別為4.59和2.295 MPa，而上、下邊界的孔隙水壓力在這二者之間線性過渡，巖體的初始孔隙率為0.3。熱?水?應(yīng)力耦合的環(huán)境對(duì)巖體顆粒要產(chǎn)生壓力溶解作用，從而引起孔隙率和滲透系數(shù)的變化。有關(guān)的計(jì)算參數(shù)如表 1和表 2所示(該表中主要數(shù)據(jù)參考文獻(xiàn)[1，10])。初始狀態(tài)時(shí)，巖體的溫度均為40 ℃。核廢物以1 kW的不變功率釋放熱量，時(shí)間經(jīng)歷4 a。

圖3 計(jì)算模型(工況1)Fig. 3 Computation model(case 1)

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

式中：α=3.86×10?6m?1，β=1.41，γ=1?1/β；ψ為水勢(shì)；sws=1.0，swr=0.19，其分別為最大飽和度和最小飽和度。比滲透率與飽和度的關(guān)系為

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

表2 壓力溶解計(jì)算參數(shù)Table 2 Parameters for pressure solution

取巖體介質(zhì)的溫度梯度水分?jǐn)U散系數(shù)為Dt=2.5×10?12m2/(s·℃)，針對(duì)前述的2種工況，其主要結(jié)果及分析如下。

由于滲流和地應(yīng)力對(duì)核廢物釋熱的影響較小[14]，故2種工況條件下計(jì)算域中的溫度場(chǎng)基本相同。以工況1為例，432，433，434和435各點(diǎn)處的溫度隨時(shí)間的變化曲線如圖4所示。從圖4可見：在開始的約0.1 a內(nèi)緩沖層的溫度快速上升，之后增加減緩，到計(jì)算終時(shí)432，433，434和435各點(diǎn)的溫度依次為97.8，81.9，72.6和 65.7 ℃。圖5所示為工況1在4a時(shí)計(jì)算域中的溫度等值線分布。

圖6所示為巖體中2個(gè)點(diǎn)的溶質(zhì)濃度?時(shí)間曲線(固化體中心x=0)。從圖6可見：在壓力溶解作用下，顆粒接觸面處和顆?？紫吨械娜苜|(zhì)濃度ci和cp在開始約0.2 a內(nèi)急速升、降，之后隨時(shí)間的推移仍在減小，但變化幅度不大。到計(jì)算終了時(shí)x為0.3 m和1.1 m處的溶質(zhì)濃度ci分另為14.65 mol/m3和15.03 mol/m3(工況 1)，15.19 mol/m3和 15.27 mol/m3(工況 2)；cp：13.03 mol/m3和 13.04 mol/m3(工況 1)，13.04 mol/m3和 13.05 mol/m3(工況 2)。

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

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

圖6 巖體中2個(gè)點(diǎn)的溶質(zhì)濃度?時(shí)間曲線(固化體中心x=0)Fig. 6 Solute concentrations versus time at two nodes (x=0 at center of vitrified waste)

圖7 巖體中2個(gè)點(diǎn)的遷移/沉淀量?時(shí)間曲線(固化體中心x=0)Fig. 7 Removal/precipitation masses versus time at two nodes (x=0 at center of vitrified waste)

圖7 所示為巖體中2個(gè)點(diǎn)的遷移/沉淀量?時(shí)間曲線(固化體中心x=0)。從圖7可見：在壓力溶解作用下，顆粒接觸面處的遷移質(zhì)量和顆?？紫吨械某恋碣|(zhì)量亦在開始約0.2 a內(nèi)急速升、降，之后隨時(shí)間的變化趨緩。正是因顆粒接觸面處的溶解和顆?？紫吨械某恋硎沟脦r體孔隙率減小。到計(jì)算終了時(shí)x為0.3 m和1.1 m處，分別為1.038 11 mol和1.279 06 mol(工況1)，1.379 43 mol和1.430 91 mol(工況 2)；分另為?1.038 13 mol和?1.279 08 mol(工況 1)，?1.379 45 mol和?1.430 94 mol(工況 2)。經(jīng)比較可知：，即顆粒接觸面處的溶解量基本變?yōu)轭w?？紫吨械某恋砹?。

圖8所示為源/匯項(xiàng)反應(yīng)體積R和顆粒貫穿深度h隨時(shí)間的變化曲線。從圖8可見：在壓力溶解作用下，R和h也是在開始約0.2 a內(nèi)急速上升，之后隨時(shí)間的推移R緩慢下降，而h還呈現(xiàn)明顯的增長(zhǎng)趨勢(shì)(孔隙率減小的原因之一)。到計(jì)算終了時(shí)x為0.3 m和1.5 m處，R分別為 8.330×10?5m3，8.227×10?5m3(工況 1)，8.341×10?5m3和 8.237×10?5m3(工況 2)；h分別為1.317×10?5m和1.502×10?5m(工況1)，1.512×10?5m和 1.615×10?5m(工況 2)。

圖8 巖體中2個(gè)點(diǎn)的反應(yīng)體積和顆粒貫穿深度?時(shí)間曲線(固化體中心x=0)Fig. 8 Reaction volumes and granular interpenetrations versus time at two nodes (x=0 at the center of vitrified waste)

圖9 2種工況在4 a時(shí)巖體中孔隙率等值線分布(單位：10?1)Fig. 9 Contours of porosity in rock mass at four years for two cases

圖 9所示為工況 1到 4 a時(shí)玻璃固化體周圍 2 m×2 m范圍內(nèi)巖體孔隙率等值線分布。此時(shí)在計(jì)算域水平對(duì)稱軸上距玻璃固化體中心0.3 m和1.1 m點(diǎn)處的孔隙率分別為 0.181和 0.154(工況 1)，0.145和0.133(工況 2)，分別約為初始值 0.3的 60%，51%和48%、44%。距離玻璃固化體越遠(yuǎn)，孔隙率減小越多，其原因?yàn)椋涸诓ＡЧ袒w附近溫度較高，由此產(chǎn)生的熱效應(yīng)使得真實(shí)應(yīng)力σa與臨界應(yīng)力之差較小，并且固體硅的化學(xué)活度與溫度T成反比，因而，該處孔隙率下降較少，反之亦然。上述2點(diǎn)處的孔隙率隨時(shí)間的變化曲線如圖10所示。從圖10可見：在開始的約0.2 a內(nèi)，因上覆及自重荷載的施加，壓力溶解使得巖體孔隙率內(nèi)急速下降，之后隨時(shí)間的推移孔隙率雖仍在衰減，但變化幅度較小(漸趨平衡)，且與溫度的增長(zhǎng)趨勢(shì)并不同步。將圖10與Taron等由解析及試驗(yàn)得出的巖體孔隙率?時(shí)間曲線[10](圖11)相比較，二者變化趨勢(shì)較一致。

圖10 巖體中2個(gè)點(diǎn)的孔隙率-時(shí)間曲線(固化體中心x=0)Fig. 10 Porosities versus time at two nodes(x=0 at center of vitrified waste)

圖11 Taron模型與試驗(yàn)數(shù)據(jù)的對(duì)比[10]Fig. 11 Comparison of Taron model with experimental data[10]

相對(duì)應(yīng)的滲透系數(shù)的等值線分布及隨時(shí)間的變化曲線見分別圖12和圖13，其與孔隙率的演化規(guī)律類似，在時(shí)間上表現(xiàn)了一定的非線性特征(由式(22)決定)。上述2點(diǎn)處的滲透系數(shù)分別為0.199×10?13m/s和 0.115×10?13m/s(工況 1)，0.093×10?13m/s 和0.070×10?13m/s(工況2)，分別約為初始值的16.0%，9.2%和7.5%，5.6%。

圖14所示為2種工況中432和434點(diǎn)處的孔隙水壓力隨時(shí)間的變化曲線。從圖14可見：開始時(shí)由于巖體荷載突然施加，孔隙水壓力有瞬時(shí)的少許下降，之后在壓力溶解(巖體孔隙率和滲透系數(shù)減小)、應(yīng)力場(chǎng)和溫度場(chǎng)的共同作用下，孔隙水壓力隨時(shí)間呈一定幅度上升。到4 a時(shí)該432和434點(diǎn)處的孔裂隙水壓力分別為4.580 1 MPa和4.580 5 MPa(比初值4.59 MPa略小)(工況1)，3.464 6 MPa和3.217 9 MPa(分別大于初值3.33 MPa和3.10 MPa)(工況2)。計(jì)算終了時(shí)2種工況中的孔隙水壓力等值線如圖15所示。

圖12 2種工況在4 a時(shí)巖體中滲透系數(shù)等值線分布(單位：10?14 m/s)Fig. 12 Contours of permeability in rock mass at four years for two cases

圖13 巖體中2個(gè)點(diǎn)的滲透系數(shù)?時(shí)間曲線(固化體中心x=0)Fig. 13 Permeabilities versus time at two nodes(x=0 at the center of vitrified waste)

圖16 所示為工況1，2在4 a時(shí)計(jì)算域中的孔隙水流速矢量分布，前者與后者的比例尺為50:1?？吹?種工況的流速矢量分布截然不同。以432，433，434和435點(diǎn)為例，孔隙水流速分別為：6.66×10?11，5.21×10?11，4.63×10?11，4.58×10?11m/s(工況 1)；5.43×10?9，5.07×10?9，4.85×10?9，4.72×10?9(工況 2)?？梢姽r2的流速比工況1約大2個(gè)數(shù)量級(jí)。

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

圖15 兩種工況下4 a時(shí)巖體中的孔隙水壓力等值線(單位：MPa)Fig. 15 Contours of pore pressures in rock mass at four years for two cases

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

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

圖17 所示為4 a時(shí)計(jì)算域中正應(yīng)力等值線。從圖17可見：2種工況的計(jì)算域中應(yīng)力量值及分布也差別很大，如到4 a時(shí)，432，433，434和435各點(diǎn)水平正應(yīng)力/垂直正應(yīng)力依次為：?9.9/?31.6 MPa，?12.5/?33.4 MPa，?14.2/?34.8 MPa 和?15.5/?35.9 MPa(工況 1)；?19.0/?41.3 MPa，?20.7/?42.9 MPa，?20.8/?43.6 MPa和?20.5/?43.7 MPa(工況 2)。其原因在于與工況 1相比，工況2中孔隙水的壓力較小但其流速較大，由此產(chǎn)生的動(dòng)、靜水壓力更顯著地?cái)_動(dòng)了巖體應(yīng)力。

4 結(jié)論

(1) 將 Taron等建立的顆粒聚集體的壓力溶解模型引入孔隙介質(zhì)熱?水?應(yīng)力耦合有限元程序中，以一個(gè)假定的實(shí)驗(yàn)室尺度且位于飽和顆粒聚集體孔隙巖體中的高放廢物地質(zhì)處置模型為例子，就考慮邊界上孔隙水壓力相等與否的兩種工況，通過熱?水?應(yīng)力耦合的二維有限元模擬，考察了巖體中的溫度、顆粒界面及孔隙中的溶質(zhì)濃度、遷移和沉淀量、孔隙率及滲透系數(shù)、孔隙水壓力、地下水流速和主應(yīng)力的變化、分布情況。

(2) 2種工況中的溫度狀態(tài)基本相同，計(jì)算終了4 a時(shí)，近場(chǎng)的溫度可達(dá)到 40.0～98.0 ℃；相比于工況 1的相同孔隙水壓力邊界條件，工況2因其計(jì)算域左、右邊界存在2.295 MPa 的孔隙水壓力差，巖體中有較大的流速(穩(wěn)態(tài)流量Qss)，促進(jìn)顆粒介質(zhì)的溶解、遷移和沉淀，使得孔隙率和滲透系數(shù)加快下降，從而對(duì)滲流場(chǎng)(孔隙水的壓力及流速)和應(yīng)力場(chǎng)產(chǎn)生顯著的影響。

[1]Yasuhara H, Elsworth D, Polak A. A mechanistic model for compaction of granular aggregates moderated by pressure solution[J]. J Geophys Res, 2003, Vol, 108, NO. B11, 2530.doi:10. 1029/2003JB002536.

[2]Tenthorey E, Cox S, Todd H. Evolution of strength recovery and permeability during fluid-rock reaction in experimental fault zones[J]. Earth and Plane Sci Lett, 2003, 206: 161?172.

[3]Moore E, Lockner A, Byerlee D. Reduction of permeability in granite at elevated temperatures[J]. Science 265, 1994:1558?1561.

[4]Lin W, Roberts J, Glassley W, et al. Fracture and matrix permeability at elevated temperatures[R]. Workshop on significant issures and available data, near-field/altered-zone coupled effects expert elicitation project. San Francisco, 1997,Novermber.

[5]Polak A, Elsworth D, Yasuhara H, et al. Permeability reduction of a natural fracture under net dissolution by hydrothermal fluids[J]. Geophys Res Lett, 2003, 30(20), 2020. doi:10.1029/2003GL017575.

[6]Durham W, Bourcier W, Burton E. Direct observation of reactive flow in a single fracture[J]. Water Resour Res, 2001, 37:1?12.

[7]Yasuhara H, Elsworth D. Evolution of permeability in a natural fracture: Significant role of pressure solution[J]. J Geophys Res,2004, 109, B03204. doi:10.1029/2003JB002663.

[8]Taron J, Elsworth D. Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs[J].Int J Rock Mech Min Sci, 2009. doi:1016/j.ijrmms.2009.01.007.

[9]Taron J, Elsworth D. Coupled mechanical and chemical processe s in engineered geothermal reservoirs with dynamic permeability[J]. International Journal of Rock Mechanics and Mining Sciences, 2010, 47: 1339?1348.

[10]Taron J, Elsworth D. 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.

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

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

[13]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[R]. DECOVALES III-Task3-Bench Mark Test 1(BMT1)-Subtask BMT1-B, 2002: 1?43.

[14]Rutqvist J, Chijimatsu M, Jing L, et al. A numerical study of THM effects on the near-field safety of a hypothetical nuclear waste repository—BMT1 of the DECOVALEX III project. Part 3: Effects of THM coupling in sparsely fractured rocks[J].International Journal of Rock Mechanics and Mining Sciences,2005, 42(5/6): 745?755.