初始孔隙率對高放射性廢物處置庫緩沖材料性狀的影響

秦愛芳, 賈 旭

(上海大學(xué)土木工程系, 上海200444)

近年來, 核能作為世界公認的清潔能源得到了充分的開發(fā)和使用. 與此同時, 安全、有效地處置高放射性核廢料也是世界各國面臨的一項重大課題. 目前, 在眾多方案中, 地質(zhì)處置[1]是一項被普遍接受的方案.

高放射性核廢料處置庫(high-level radioactive nuclear waste repsitory, HLWR)通常由廢物罐、玻璃固化體、緩沖材料、回填材料和圍巖組成. HLWR 關(guān)閉運行后, 通常需要長達千年甚至萬年以上處置期限. HLWR 封閉運行后高放射性廢料產(chǎn)生的熱量以及地下水壓產(chǎn)生的高靜水壓力, 會使廢物罐和圍巖之間的緩沖材料(主要為膨潤土)產(chǎn)生很復(fù)雜的熱-水-力(thermohydro-mechanical, THM)耦合作用. HLWR 中緩沖材料的選擇對整個處置庫是否能夠安全、可靠地運行至關(guān)重要.

目前, 許多歐美國家使用膨潤土作為HLWR 的緩沖材料, 利用膨潤土的低滲透性、膨脹性、充足的導(dǎo)熱性和高交換性[2]來防止核泄漏. 中國擬采用GMZ 膨潤土作為緩沖材料, 同時選擇北山花崗巖作為圍巖. GMZ 膨潤土除了具有高膨脹性和導(dǎo)熱性以外, 其高陽離子交換能力和高蒙脫石含量會進一步阻止核泄漏的發(fā)生[3].

Villar[4]和Komine[5]預(yù)測了不同孔隙率下飽和膨潤土的滲透系數(shù)變化曲線. Olivella等[6]和Gens 等[7]將高放射性處置庫周圍的FEBEX 膨潤土孔隙率分為兩部分: 宏觀孔隙部分和微觀孔隙部分, 并利用數(shù)值模擬分析了處置庫在側(cè)限條件下的飽和度、應(yīng)力、溫度變化情況.葉為民等[8]、劉月妙等[9]、張玉軍等[10]、孫德安等[11]也對GMZ 膨潤土的導(dǎo)熱性、力學(xué)特性、水力學(xué)特性、微觀孔隙特性等進行了研究. 宋杰等[12]探究了壓實度對非飽和土電阻率的影響.Ye 等[13-14]和Tang 等[15]研究了GMZ 膨潤土中, 溫度對吸力和液體滲透系數(shù)的影響.

膨潤土壓實度的不同, 會造成不同的初始孔隙率, 進而影響HLWR 關(guān)閉運行后吸力、飽和度、滲透系數(shù)、應(yīng)力、溫度等的變化. 若初始孔隙率太小, 可能會造成HLWR 關(guān)閉運行后內(nèi)部壓力過大, 側(cè)限條件下可能會對圍巖或開挖擾動區(qū)產(chǎn)生較大應(yīng)力和應(yīng)變, 致使圍巖開裂, 核素隨著地下水外泄等情況. 若初始孔隙率過大, 膨潤土無法起到緩沖和回填材料的作用, 可能會造成HLWR 關(guān)閉運行后膨潤土很難膨脹密實, 膨潤土和圍巖之間產(chǎn)生縫隙[16].

HLWR 作為各個國家的重大安全工程, 對人民生活和國家安全有著重要影響. 我國核廢料處置庫研究正處于初步階段, 大多針對緩沖材料GMZ 膨潤土的水-土特性和圍巖的水、力學(xué)特性展開, 處置庫現(xiàn)場試驗也在準(zhǔn)備階段. 因此, 考慮我國高廟子膨潤土特性, 建立一套合理可行的數(shù)值模擬, 用來預(yù)測或反映處置庫不同條件下的運行情況至關(guān)重要.

1 THM 三相耦合機理

當(dāng)HLWR 封閉之后, 膨潤土將接收由廢物罐不斷散發(fā)出的熱量, 使其溫度升高. 同時, 地下水在水壓的驅(qū)動下不斷向膨潤土內(nèi)滲透, 使膨潤土內(nèi)的水壓和飽和度產(chǎn)生變化. 膨潤土內(nèi)的溫度和飽和度發(fā)生改變后, 自身應(yīng)力也會產(chǎn)生復(fù)雜的變化. 在熱、水、力三相相互作用下, 膨潤土內(nèi)部就產(chǎn)生了一個非常復(fù)雜的THM 耦合現(xiàn)象. 例如, 距離廢物罐較近的膨潤土內(nèi)的水分受熱后將蒸發(fā)成水蒸氣, 該處的膨潤土由于失去水分使得土內(nèi)孔隙變小, 遠處的水分更難以滲透, 其飽和度及水壓力減小. 同時, 該處的膨潤土溫度較高, 又會產(chǎn)生較高的溫度應(yīng)力.

HLWR 中除了THM 耦合機理復(fù)雜外, THM 耦合作用的時間也較長. HLWR 內(nèi)部儲存的高放射性核廢料含有眾多放射性元素, 如銫、鈾、釙等都擁有數(shù)以萬年計的半衰期, 故HLWR中的THM 耦合會持續(xù)幾千年甚至幾萬年.

2 本構(gòu)方程

本工作中的熱、水本構(gòu)方程均引自文獻[7].

2.1 熱本構(gòu)方程

熱傳導(dǎo)假定服從Fourier 定律:

式中, ic為介質(zhì)傳導(dǎo)熱通量, ▽T 為介質(zhì)溫度的微分, λ 為介質(zhì)綜合熱傳導(dǎo)系數(shù)(與土熱傳導(dǎo)系數(shù)λs、水熱傳導(dǎo)系數(shù)λl、孔隙率以及飽和度有關(guān)), 其關(guān)系式為

其中sl為膨潤土的飽和度, λsat為膨潤土完全飽和時的熱傳導(dǎo)系數(shù), λdry為膨潤土干燥情況下的熱傳導(dǎo)系數(shù).

2.2 流體本構(gòu)方程

孔隙中液相遵循Darcy 流動定律, 即

式中: ql為液相的流量; Kl為液相滲透系數(shù)張量, 且Kl= KrlK/μl, 其中Krl為液相相對滲透系數(shù), μl為液相的動力黏滯系數(shù); g 為重力向量; K 為介質(zhì)固有滲透系數(shù); ▽pl為水壓力的微分; rl為液相的密度.

以上計算采用的幾個關(guān)系式如下.

(1) 固有滲透系數(shù)K 和孔隙率相關(guān), 采用Kozeny 模型, 即

式中, φ0為初始孔隙率, K0為φ0時的固有滲透系數(shù).

(2) 非飽和土的液相相對滲透系數(shù)Krl與其飽和度密切相關(guān), 即

式中: A 通常為1; n 為指數(shù), 通常取2～4; Se為有效飽和度. Se與實際飽和度Sl之間的關(guān)系如下:

式中, Slr為液體剩余飽和度, Sls為液相最大飽和度.

(3) 為了反映飽和度與基質(zhì)吸力之間的關(guān)系, 采用比較常用的Van Genuchten 模型, 即

式中, p=P0σ/σ0, 其中P0為參照溫度下的進氣值, σ0為參照溫度下的表面張力, σ 和p 為計算溫度下的表面張力和進氣值, m 為形狀參數(shù).

2.3 力學(xué)本構(gòu)方程

非飽和土的應(yīng)力應(yīng)變本構(gòu)模型非常復(fù)雜. 本工作采用巴塞羅那基本模型(Barcelona basic model, BBM). BBM 模型是在修正的劍橋模型基礎(chǔ)上, 由Alonso 等[17]進一步修改完成的非飽和土彈塑性本構(gòu)模型. BBM 模型中選用兩個屈服面, 考慮了吸力的影響, 并假定吸力增加會提高非飽和土的黏聚力但仍保持內(nèi)摩擦角為常數(shù). BBM 模型體積應(yīng)變的計算關(guān)系式為

式中,

其中κi0是初始零吸力時的平均應(yīng)力, κs0為初始吸力值, αi和αsp為試驗參數(shù), α0為彈性熱應(yīng)變的試驗參數(shù), Tref為初始溫度, dT 為溫度的變化量.

BBM 模型采用雙屈服面模型. 屈服面1 為

式中, p=(σ1+2σ3)/3-pa, q =σ1-σ3, M 為飽和土臨界狀態(tài)線的切線值,

式中：CM為測量成本；CZ為測量誤差成本；CMF為測量儀器成本；CMM為測量策略成本；CEA為棄真誤差成本；CEB為納偽誤差成本；CEA為單個棄真事件成本；cEB為單個納偽事件成本；N為待測區(qū)域關(guān)鍵部位的個數(shù)；P(∪TfA)為公差T時棄真事件發(fā)生的概率；P(∪TfB)為公差T時納偽事件發(fā)生的概率。

pc為一參照應(yīng)力, p*0為土飽和狀態(tài)下的前期固結(jié)壓力, p0為對應(yīng)于吸力s 情況下的前期固結(jié)壓力, λ(0)為飽和狀態(tài)下土ν-ln p 坐標(biāo)體系中的原始壓縮曲線斜率, λ(s)為吸力為s 情況下土在ν-ln p 坐標(biāo)體系中的原始壓縮曲線斜率, β 為控制λ(s)隨吸力增長速率的參數(shù), k 為因吸力導(dǎo)致張力變化的參數(shù).

屈服面2 為

式中, s0是土體受到的前期最大吸力. 當(dāng)吸力從s0過渡到0 時, 土體也將從彈性狀態(tài)轉(zhuǎn)變到彈塑性狀態(tài).

2.4 平衡方程

(1) 水的質(zhì)量平衡方程為式中, θwl和θwg分別為單位體積中液相中和氣相中水的質(zhì)量, φ 為孔隙率, Sl和Sg分別代表液相和氣相的飽和度, Jwl和Jwg分別代表水在液相和氣相中的總質(zhì)量通量, fw為外部水源供應(yīng)量.

(2) 空氣的質(zhì)量守恒方程為

式中, θal和θag分別表示單位體積液相和氣相中空氣的質(zhì)量, fa為單位體積介質(zhì)中外界提供空氣的質(zhì)量, Jal和Jag分別表示液相和氣相中空氣流動的質(zhì)量(相對于固定參照系).

(3) 在多項耦合條件下有多種物質(zhì)存在, 考慮到在非飽和土中有多種物質(zhì)轉(zhuǎn)移, 故采用如下能量平衡方程:

(4) 忽略慣性力的靜力平衡方程為

式中, σ 為應(yīng)力狀態(tài)向量, B 為體力向量.

對于飽和多孔介質(zhì), 力學(xué)本構(gòu)方程引入有效應(yīng)力概念, 即

式中, MT= (1,1,1,0,0,0). 但對于非飽和多孔介質(zhì)的本構(gòu)方程, 通常引入兩個獨立的應(yīng)力狀態(tài)量, 即凈法向應(yīng)力σ-pgM 以及基質(zhì)吸力s=pg-pl, 其中pg為氣壓力, pl為水壓力.

3 數(shù)值分析算例

本算例采用西班牙ENRESA 在瑞士核廢料地下試驗室進行的核廢料處置庫模型, 如圖1所示. 本試驗在地下約500 m 的巖層中開挖一個水平的直徑為2.4 m、長為18.8 m 的坑道作為處置庫; 沿水平軸向放置兩個半徑為0.45 m、長為4.54 m 的廢物罐; 廢物罐周圍用緩沖材料GMZ01 膨潤土填充. 軸對稱計算模型如圖1(c) 所示, 其中CL 表示中心線, CL 右側(cè)分別為熱源、GMZ01 膨潤土、花崗巖, 在GMZ01 膨潤土的左側(cè)為銅罐, 即熱源、邊界條件和初始條件如下: GMZ01 膨潤土的初始溫度為20?C, 初始各向應(yīng)力為0.5 MPa, 初始吸力為90 MPa,初始飽和度為0.46, 初始孔隙率分別取0.3, 0.4, 0.5, 大氣壓恒為0.1 MPa; 花崗巖的初始溫度為20?C, 初始各向應(yīng)力為28 MPa, 地下埋深為500 m, 初始水壓為0.5 MPa, 初始孔隙率為0.01, 邊界始終保持0.5 MPa 的水壓, 用來模擬地下水不斷向圍巖GMZ01 和膨潤土補充水分.GMZ01 膨潤土的基本參數(shù)如表1[18]所示. 本工作采用的基本控制方程為式(18) ～(21).

模型的熱、力學(xué)邊界條件如下: 左邊界在水平方向加約束, 上、下邊界在垂直方向加約束，右邊界在水平方向施加橫向水平壓力; 熱源從第一天開始, 不斷加熱, 第21 天時溫度達到100?C, 保持100?C 至第1 095 天.

為了監(jiān)測模型中溫度、飽和度、孔隙率、應(yīng)力、吸力等變化情況, 在徑向設(shè)定了2 個監(jiān)測點(見圖2) , 其中a 點距離熱源中心0.57 m, b 點距離熱源中心1.02 m(靠近花崗巖處). 該計算模型有180 個單元, 273 個節(jié)點. 圖2 顯示了計算結(jié)果輸出點和單元的位置.

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

圖1 HLWR 原位試驗?zāi)Ｐ虵ig.1 In situ test model of HLWR

4 計算結(jié)果與分析

4.1 溫度的變化

本工作研究了不同孔隙率(φ = 0.3, 0.4, 0.5)情況下溫度、飽和度、應(yīng)力等的變化. 根據(jù)表1 的參數(shù), φ=0.3 的GMZ01 膨潤土干密度為1.995 g/cm3, φ=0.4 的GMZ01 膨潤土干密度為1.710 g/cm3, φ=0.5 的GMZ01 膨潤土干密度為1.425 g/cm3.

圖2 計算結(jié)果輸出點和單元的位置(m)Fig.2 Position of point and unit of calculating output

圖3(a)和(b)為圖2 中a, b 兩點第21 天到第1 095 天的溫度變化曲線. 可以看出, 靠近花崗巖處的GMZ01 膨潤土初始孔隙率越小, 溫度就越高, 這是因為干密度越大導(dǎo)熱系數(shù)越大,接收的熱量就越快. 當(dāng)干密度大于1.4 g/cm3時, 比熱容隨干密度的增加略有減小[8]. 這意味著初始孔隙率越小, 接收的熱量越多; 比熱容越小, 溫度越高. 靠近花崗巖處, 溫度變化相差雖然不大, 但是干密度越大, 熱擴散系數(shù)越大[8], 所以在靠近花崗巖處的GMZ01 膨潤土初始孔隙率越小, 溫度也略高.

圖3 a, b 兩點GMZ01 膨潤土溫度隨時間變化曲線Fig.3 Variation curves of GMZ01 bentonite temperature with time at a and b points

4.2 飽和度和吸力的變化

圖4(a)為不同初始孔隙率膨潤土在a 點的飽和度隨時間變化曲線. 可以看出, a 點處GMZ01 膨潤土的飽和度先減小, 隨著時間的增長略有升高. 因為在靠近熱源處溫度較高,GMZ01 膨潤土在熱源加熱下, 水分蒸發(fā), 飽和度隨之降低. 之后由于巖石水分的不斷補充, 再次使飽和度上升. 孔隙率較低的GMZ01 膨潤土, 滲透系數(shù)較小, 水分更難滲透進來, 使得a 點處的飽和度出現(xiàn)了孔隙率越小飽和度越小的現(xiàn)象.

圖4(b)為不同初始孔隙率GMZ01 膨潤土在b 點的飽和度隨時間變化曲線. 可以看出, b點處GMZ01 膨潤土的飽和度總體隨時間呈上升趨勢. GMZ01 膨潤土中靠近花崗巖處的水分主要來自于兩個部分: 一部分來自于靠近熱源處的水分蒸發(fā)形成的水蒸氣的遷移; 另一部分來自于圍巖水分的滲透. 這兩個來源使得靠近花崗巖處的GMZ01 膨潤土的飽和度不斷上升. 由圖4(a)和(b)可以看出來, 不同的初始孔隙率在a, b 兩點表現(xiàn)出的規(guī)律是相反的. 在靠近熱源處, 初始孔隙率越大飽和度越大; 在靠近花崗巖處, 初始孔隙率越大飽和度越小. 由于熱源相同, 不同孔隙率的GMZ01 膨潤土的初始飽和度相同, 蒸發(fā)的水分是等量的. 在a 點, 孔隙率越小的GMZ01 膨潤土則飽和度降低得越多; 當(dāng)遷移至b 點, 孔隙率較小的GMZ01 膨潤土更加容易飽和.

圖4(c)和(d)為不同初始孔隙率GMZ01 膨潤土在a, b 兩點的吸力隨時間變化曲線. 可見:a 點GMZ01 膨潤土的吸力隨時間先增加后減小, 而b 點GMZ01 膨潤土的吸力是不斷減小的;靠近熱源處的吸力隨初始孔隙率的增大而減小, 而靠近花崗巖處的吸力隨初始孔隙率的增大而增大. 吸力和飽和度的關(guān)系是緊密相連的, 飽和度越大則吸力越小, 反之, 飽和度越小則吸力越大. 這也和飽和度的變化規(guī)律相互印證.

根據(jù)已有的對GMZ01 膨潤土吸力和滲透系數(shù)的研究可知: 當(dāng)吸力大于65 MPa 時, 水的滲透系數(shù)隨吸力的增大而增大; 當(dāng)吸力小于60 MPa 時, 滲透系數(shù)隨吸力的增大而減小; 當(dāng)吸力小于20 MPa 時, 滲透系數(shù)隨吸力的增大而減小的幅度會更大[12]. 由圖4(c)和(d)可以看出:靠近熱源的GMZ01 膨潤土處于高吸力狀態(tài)下, 該點附近GMZ01 膨潤土的滲透系數(shù)隨吸力的增大而增大, 所以該點GMZ01 膨潤土的滲透系數(shù)隨著初始孔隙率的增大而減小; 靠近花崗巖處GMZ01 膨潤土的吸力大部分小于60 MPa, 所以該處的滲透系數(shù)隨著吸力的增加而減小,即滲透系數(shù)隨著初始孔隙率的增大而減小. 這也是造成圖4(a)和(b)飽和度情況的原因之一.

圖4 a, b 兩點GMZ01 膨潤土飽和度和吸力隨時間變化曲線Fig.4 Variation curves of GMZ01 bentonite saturation and suction with time at a and b points

4.3 應(yīng)力的變化

圖5(a)和(b)為不同初始孔隙率的GMZ01 膨潤土在第21 天和第1 095 天的徑向應(yīng)力隨距離的變化關(guān)系. 可以看出: 第21 天時不同初始孔隙率的GMZ01 膨潤土的應(yīng)力變化不大; 靠近廢物罐和靠近花崗巖處的GMZ01 膨潤土應(yīng)力較大, 中間部分的應(yīng)力較小. GMZ01 膨潤土應(yīng)力的增加, 主要與兩個因素有關(guān): 一是由于飽和度的增加會導(dǎo)致GMZ01 膨潤土內(nèi)蒙脫石吸水,發(fā)生水化作用并迅速產(chǎn)生膨脹; 二是來自于熱源產(chǎn)生的熱應(yīng)力, 在靠近廢物罐處的GMZ01 膨潤土除了因水膨脹以外, 還受到來自于廢物罐的熱應(yīng)力. 在靠近花崗巖處的GMZ01 膨潤土受到的熱應(yīng)力雖然較小, 但由于飽和度較高會產(chǎn)生較大的膨脹應(yīng)力. 靠近廢物罐處的GMZ01 膨潤土初始孔隙率越大, 溫度越低, 熱應(yīng)力越小; 靠近花崗巖處的GMZ01 膨潤土初始孔隙率越大, 飽和度越低, 膨脹應(yīng)力越低, 所以GMZ01 膨潤土的應(yīng)力隨著初始孔隙率的增加而減小.

圖5 第21 天和第1 095 天時徑向應(yīng)力隨距離的變化關(guān)系Fig.5 Variation curves of stress with distance in the 21st and 1 095th days

4.4 數(shù)據(jù)驗證

圖6 為THM 耦合, 初始孔隙率為0.4 時, 距離熱源0.465, 0.8, 1.135 m 的c, d, e 處的飽和度計算結(jié)果以及與Gens 等[7]的結(jié)果對比. 可以看出, 二者的飽和度變化規(guī)律相同, 結(jié)果相似,從而驗證了本方法的合理性.

圖6 飽和度曲線對比結(jié)果Fig.6 Comparisions of the saturation curves

5 結(jié) 論

本工作通過有限元軟件Code-Bright 進行數(shù)值模擬, 研究了不同初始孔隙率下的溫度、飽和度、吸力、應(yīng)力的變化規(guī)律, 并得到以下結(jié)論.

(1) 就HLWR 內(nèi)GMZ01 膨潤土溫度場而言, 初始孔隙率或壓實度的不同會對靠近熱源處GMZ01 膨潤土的溫度產(chǎn)生明顯變化. 在靠近熱源處, GMZ01 膨潤土內(nèi)溫度隨初始孔隙率的減小而增大; 在靠近花崗巖處, 溫度隨初始孔隙率的變化不大.

(2) 就HLWR 內(nèi)GMZ01 膨潤土滲流場而言, 初始孔隙率的不同會導(dǎo)致明顯的滲流場變化. 在靠近熱源處, 初始孔隙率越大的GMZ01 膨潤土其飽和度越大, 吸力越小; 在靠近花崗巖處, 初始孔隙率越大的GMZ01 膨潤土其飽和度越小, 吸力越大.

(3) 就HLWR 內(nèi)GMZ01 膨潤土應(yīng)力而言, 初始孔隙率不同會造成應(yīng)力有較大不同. 第21天, 不同的初始孔隙率, 應(yīng)力變化較小; 第1 095 天, 初始孔隙率越大, 應(yīng)力越小.