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

秦愛芳, 賈 旭

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

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

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

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

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

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

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

1 THM 三相耦合機理

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

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

2 本構方程

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

2.1 熱本構方程

熱傳導假定服從Fourier 定律:

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

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

2.2 流體本構方程

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

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

以上計算采用的幾個關系式如下.

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

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

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

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

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

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

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

2.3 力學本構方程

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

式中,

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

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

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

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

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

屈服面2 為

式中, s0是土體受到的前期最大吸力. 當吸力從s0過渡到0 時, 土體也將從彈性狀態轉變到彈塑性狀態.

2.4 平衡方程

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

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

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

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

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

式中, σ 為應力狀態向量, B 為體力向量.

對于飽和多孔介質, 力學本構方程引入有效應力概念, 即

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

3 數值分析算例

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

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

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

表1 主要計算參數Table 1 Main computation parameters

圖1 HLWR 原位試驗模型Fig.1 In situ test model of HLWR

4 計算結果與分析

4.1 溫度的變化

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

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

圖3(a)和(b)為圖2 中a, b 兩點第21 天到第1 095 天的溫度變化曲線. 可以看出, 靠近花崗巖處的GMZ01 膨潤土初始孔隙率越小, 溫度就越高, 這是因為干密度越大導熱系數越大,接收的熱量就越快. 當干密度大于1.4 g/cm3時, 比熱容隨干密度的增加略有減小[8]. 這意味著初始孔隙率越小, 接收的熱量越多; 比熱容越小, 溫度越高. 靠近花崗巖處, 溫度變化相差雖然不大, 但是干密度越大, 熱擴散系數越大[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 膨潤土在熱源加熱下, 水分蒸發, 飽和度隨之降低. 之后由于巖石水分的不斷補充, 再次使飽和度上升. 孔隙率較低的GMZ01 膨潤土, 滲透系數較小, 水分更難滲透進來, 使得a 點處的飽和度出現了孔隙率越小飽和度越小的現象.

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

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

根據已有的對GMZ01 膨潤土吸力和滲透系數的研究可知: 當吸力大于65 MPa 時, 水的滲透系數隨吸力的增大而增大; 當吸力小于60 MPa 時, 滲透系數隨吸力的增大而減小; 當吸力小于20 MPa 時, 滲透系數隨吸力的增大而減小的幅度會更大[12]. 由圖4(c)和(d)可以看出:靠近熱源的GMZ01 膨潤土處于高吸力狀態下, 該點附近GMZ01 膨潤土的滲透系數隨吸力的增大而增大, 所以該點GMZ01 膨潤土的滲透系數隨著初始孔隙率的增大而減小; 靠近花崗巖處GMZ01 膨潤土的吸力大部分小于60 MPa, 所以該處的滲透系數隨著吸力的增加而減小,即滲透系數隨著初始孔隙率的增大而減小. 這也是造成圖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 應力的變化

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

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

4.4 數據驗證

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

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

5 結 論

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

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

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

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