存在邊界流時壓力溶解對顆粒聚集巖體中THM耦合作用的有限元模擬

張玉軍，楊朝帥,

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

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

1 混合質量平衡

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

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

圖2 壓縮及顆粒和生長的膠結物半徑的關系示意[10]Fig. 2 Diagram of compaction and relationship between granular and growing cement radius[10]

其迭代求解式為

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

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

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

并且有

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

而

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

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

Taron等給出的孔隙率表達式為

式中：εV為有限的體積應變；R為作為源/匯項的隨時間求和的反應體積。

并且

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

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

將上述顆粒聚集體的壓力溶解模型引入文獻[12]的孔隙介質熱?水?應力耦合有限元程序中。

3 算例

實驗室尺度的核廢料玻璃固化體埋置模型如圖 3所示，其周圍的介質是飽和的石英顆粒聚集孔隙巖體。作為近似簡化，認為這是一個平面應變問題。取計算域為水平向4 m，垂直向8 m，有800個單元，861個節(jié)點。從固化體邊緣向右的點號依次為432，433，434，435和436。

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

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

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

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

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

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

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

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

圖6所示為巖體中2個點的溶質濃度?時間曲線(固化體中心x=0)。從圖6可見：在壓力溶解作用下，顆粒接觸面處和顆?？紫吨械娜苜|濃度ci和cp在開始約0.2 a內急速升、降，之后隨時間的推移仍在減小，但變化幅度不大。到計算終了時x為0.3 m和1.1 m處的溶質濃度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 溫度?時間曲線Fig. 4 Temperatures versus time at some nodes

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

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

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

圖7 所示為巖體中2個點的遷移/沉淀量?時間曲線(固化體中心x=0)。從圖7可見：在壓力溶解作用下，顆粒接觸面處的遷移質量和顆?？紫吨械某恋碣|量亦在開始約0.2 a內急速升、降，之后隨時間的變化趨緩。正是因顆粒接觸面處的溶解和顆粒孔隙中的沉淀使得巖體孔隙率減小。到計算終了時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所示為源/匯項反應體積R和顆粒貫穿深度h隨時間的變化曲線。從圖8可見：在壓力溶解作用下，R和h也是在開始約0.2 a內急速上升，之后隨時間的推移R緩慢下降，而h還呈現(xiàn)明顯的增長趨勢(孔隙率減小的原因之一)。到計算終了時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個點的反應體積和顆粒貫穿深度?時間曲線(固化體中心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時巖體中孔隙率等值線分布(單位：10?1)Fig. 9 Contours of porosity in rock mass at four years for two cases

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

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

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

相對應的滲透系數(shù)的等值線分布及隨時間的變化曲線見分別圖12和圖13，其與孔隙率的演化規(guī)律類似，在時間上表現(xiàn)了一定的非線性特征(由式(22)決定)。上述2點處的滲透系數(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點處的孔隙水壓力隨時間的變化曲線。從圖14可見：開始時由于巖體荷載突然施加，孔隙水壓力有瞬時的少許下降，之后在壓力溶解(巖體孔隙率和滲透系數(shù)減小)、應力場和溫度場的共同作用下，孔隙水壓力隨時間呈一定幅度上升。到4 a時該432和434點處的孔裂隙水壓力分別為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)。計算終了時2種工況中的孔隙水壓力等值線如圖15所示。

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

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

圖16 所示為工況1，2在4 a時計算域中的孔隙水流速矢量分布，前者與后者的比例尺為50:1。看到2種工況的流速矢量分布截然不同。以432，433，434和435點為例，孔隙水流速分別為：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個數(shù)量級。

圖14 巖體中2個點的孔隙水壓力?時間曲線Fig. 14 Pore pressures versus time at two nodes for two cases

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

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

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

圖17 所示為4 a時計算域中正應力等值線。從圖17可見：2種工況的計算域中應力量值及分布也差別很大，如到4 a時，432，433，434和435各點水平正應力/垂直正應力依次為：?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中孔隙水的壓力較小但其流速較大，由此產生的動、靜水壓力更顯著地擾動了巖體應力。

4 結論

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

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

[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]張玉軍. 廢料地質處置近場熱?水?應力?遷移耦合二維有限元分析[J]. 巖土工程學報, 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.