高放廢物處置庫緩沖層分區(qū)設(shè)計(jì)與導(dǎo)熱性能數(shù)值模擬

許 迅 明華軍，2 談云志

(1. 三峽大學(xué) 特殊土資源化利用宜昌重點(diǎn)實(shí)驗(yàn)室， 湖北 宜昌 443002； 2. 三峽大學(xué) 水利與環(huán)境學(xué)院， 湖北 宜昌 443002)

核能的開發(fā)利用產(chǎn)生了大量的高放射性廢物(簡稱“高放廢物”)．如何安全地處置高放廢物，已成為國家核能利用必須解決的問題[1]．深地質(zhì)處置方法被國際上公認(rèn)為永久處置高放廢物的最佳措施[2]，即將高放廢物固體化后，放置于距地表深約500～1 000 m的完整地質(zhì)體坑道中，在其間設(shè)置緩沖回填材料，形成人工屏障，阻隔放射性核素遷移，實(shí)現(xiàn)放射性核素與生物圈的永久隔離．緩沖層是填充在廢物罐和地質(zhì)體之間的重要人工屏障，是處置庫安全性和穩(wěn)定性的有效保障[3]，因此，對(duì)處置庫導(dǎo)熱功能的改進(jìn)，主要著重于對(duì)緩沖層的設(shè)計(jì)．

膨潤土具有低滲透、高膨脹和強(qiáng)吸附的特性，是高放廢物處置庫首選的緩沖材料[4]．膨潤土導(dǎo)熱性能相對(duì)較差，即使高壓實(shí)(ρd=1.80 g/cm3)狀態(tài)，也只有1.3 W/(m·K)，遠(yuǎn)低于常規(guī)地質(zhì)體熱傳導(dǎo)系數(shù)(2.43 W/(m·K))[5-6]，導(dǎo)致高放廢物所產(chǎn)生的衰變熱消散緩慢，引起熱量聚積．據(jù)Lee等[7]計(jì)算，處置庫封場后前10年，廢物罐表面溫度最高至107℃，超出處置庫安全運(yùn)行臨界值(t=100℃)[8]，導(dǎo)致水分蒸發(fā)，進(jìn)而引起膨潤土收縮，促使局部裂隙衍生，損傷緩沖層屏障功能[9]，因此，如何使衰變熱快速向周圍擴(kuò)散是高放廢物處置庫安全運(yùn)營的關(guān)鍵問題之一．

許多學(xué)者采用在膨潤土中摻入具有高導(dǎo)熱性能的材料(石英砂、巖屑、石墨等)的方法來提高膨潤土的導(dǎo)熱性能[4，10-12]．葉為民， 王瓊，潘虹，等[4]利用熱探針法研究了干密度、含水率對(duì)加砂高廟子膨潤土導(dǎo)熱系數(shù)的影響，發(fā)現(xiàn)在其他參數(shù)相同的條件下，導(dǎo)熱系數(shù)隨干密度、含水率及摻砂量的增加而增大．馬超[10]和謝敬禮，馬利科，高玉峰，等[11]均通過添加花崗巖巖屑改善高廟子膨潤土的導(dǎo)熱性能，結(jié)果表明添加巖屑提升導(dǎo)熱性能效果要劣于添加石英砂，但巖屑作為添加劑對(duì)處置庫水化學(xué)環(huán)境影響小，相容性較好．劉月妙，蔡美峰，王駒[12]研究了在GMZ01膨潤土中分別摻入石英砂、石墨后的導(dǎo)熱性能，發(fā)現(xiàn)石英砂和石墨均可以明顯提高GMZ01膨潤土的導(dǎo)熱性能．

石英砂等材料在提升膨潤土導(dǎo)熱性能的同時(shí)，也顯著地降低了其防滲能力．Cho等[13]進(jìn)行了不同摻砂率的Kuungju膨潤土-石英砂混合物滲透試驗(yàn)，發(fā)現(xiàn)摻砂率70%時(shí)混合物的滲透系數(shù)在10-10m/s數(shù)量級(jí)左右，不能滿足國際原子能機(jī)構(gòu)(IAEA)對(duì)緩沖材料滲透性能的要求(<10-12m/s)[14]．因此，擬將緩沖層分為兩層，即外部緩沖層采用純膨潤土材料，以保持防滲性能，而將摻入劑(如：石英砂、石墨)摻入內(nèi)部緩沖層的膨潤土中，以增加緩沖材料的導(dǎo)熱性能．本研究主要通過數(shù)值模擬的手段對(duì)處置庫工程屏障進(jìn)行熱力學(xué)分析，進(jìn)而評(píng)估緩沖層分區(qū)設(shè)計(jì)的可行性．

1 雙層緩沖層概念

為了有效提高緩沖材料導(dǎo)熱性能，同時(shí)保證其防滲性能也達(dá)到要求，提出了雙層緩沖層的概念．圖1為單層緩沖層和雙層緩沖層的概念模型．假設(shè)內(nèi)部緩沖區(qū)的大小是原始緩沖區(qū)的一半．在雙層緩沖層中，外部緩沖層為純膨潤土，內(nèi)部緩沖層為膨潤土-摻入劑的混合物．對(duì)于雙層緩沖層，由于外部緩沖層依然為純膨潤土，故保持了原有的防滲功能，可以有效減少地下水的滲入．

圖1 高放廢物處置庫的單層與雙層緩沖層概念模型(單位：mm)

2 數(shù)值模型

2.1 熱傳導(dǎo)方程

熱力學(xué)第一定律指出，微元體Ω內(nèi)動(dòng)能KΩ和內(nèi)能EΩ的變化率與機(jī)械功率Pext和微元體之間的熱量交換率Qexch相等，用數(shù)學(xué)語言可表示為:

(1)

在微元體Ω內(nèi)，假設(shè)單位質(zhì)量的內(nèi)能為E，則微元體的總內(nèi)能EΩ為：

(2)

其中：ρ為材料密度(kg/m3)；cp為材料比定壓熱容，表示單位質(zhì)量材料升高單位溫度所需要的熱量[J/(kg·K)]；T為材料溫度(K)．

由于計(jì)算不考慮外力影響，故Pext=0．

熱量的交換功率Qexch包括通過傳導(dǎo)、輻射、對(duì)流和熱源作用下的熱量交換，由于固相介質(zhì)中熱量的傳遞主要通過熱傳導(dǎo)進(jìn)行，所以在表面積為S的微元體Ω中交換功率可表示為：

(3)

其中：λ為材料熱傳導(dǎo)系數(shù)[W/(m·K)]．

結(jié)合式(1)～(3)，可得材料的熱傳導(dǎo)方程為：

(4)

由于土體材料具有三相性，且各相的熱力學(xué)特性均不同，由比定壓熱容定義，可得如下公式：

ρcp=ρdcp，s+ρwnSrcp，w+ρa(bǔ)n(1-Sr)cp，a

(5)

其中：ρd為土的干密度；ρw為水的密度；ρa(bǔ)為空氣的密度(kg/m3)；n為孔隙比；Sr為飽和度；cp，s、cp，w、cp，a分別為土體顆粒、水、空氣的比定壓熱容[J/(kg·K)]．具體物理參數(shù)見表1．

表1 熱力學(xué)參數(shù)

由于土體熱傳導(dǎo)系數(shù)與土顆粒、土中水、土中氣賦存狀態(tài)有關(guān)，故本文中純膨潤土[7]和摻砂膨潤土[15]的熱傳導(dǎo)系數(shù)可分別采用如下公式：

λ=4.588w+0.124

(6)

(7)

其中：w為初始含水率；λw和λa分別表示水和空氣的熱傳導(dǎo)系數(shù)；Φs、Φb、Φw、Φa分別為砂、膨潤土、水和空氣的體積分?jǐn)?shù)．

Q為體積熱源(W/m3)，其衰減函數(shù)可表示為[7]：

Q(t)=6.353×103(t+30)-0.758

(8)

其中：t為處置庫運(yùn)行時(shí)間(年)．

2.2 邊界及初始條件

圖2為緩沖區(qū)的俯視圖，根據(jù)該圖對(duì)模型水平方向的邊界條件進(jìn)行定義．

圖2 高放廢物處置庫的雙層緩沖層概念模型俯視圖

1)邊界條件

由上文已知，廢物罐釋熱為高放廢物的衰變熱，故模型內(nèi)邊界，即廢物罐表面(r=a)處熱通量隨時(shí)間衰減，相應(yīng)邊界條件方程見式(8)．

為避免外邊界對(duì)計(jì)算結(jié)果的影響，需將外邊界設(shè)置在足夠遠(yuǎn)處，定義20倍熱源直徑為外邊界，且為恒溫邊界，不受溫度場影響，表示為：

T(∞，0)=T0

(9)

其中：T0為初始溫度，由式(10)確定．

2)初始條件

模型中各處初始溫度均為所處地深的地溫，該處地表溫度假設(shè)為15℃，地層溫度梯度為3℃/100 m，即：

T0=15+0.03z

(10)

其中：z為埋深(m)．

3)連續(xù)性條件

假設(shè)內(nèi)部緩沖材料與外部緩沖材料在r=b處完全接觸，且外部緩沖材料與圍巖在r=c處完全接觸，則必須滿足連續(xù)性條件：

T1(b，t)=T2(b，t)

(11)

T2(c，t)=T3(c，t)

(12)

(13)

(14)

其中：T1、T2、T3分別為Buffer1、Buffer2、Rock的溫度；λ1、λ2、λ3分別為Buffer1、Buffer2、Rock的熱傳導(dǎo)系數(shù)．例：T1(b，t)為時(shí)間為t時(shí)在r=b處Buffer1的溫度．

2.3 幾何模型

由于本研究主要關(guān)注工程屏障的熱力學(xué)行為，而不是整個(gè)處置系統(tǒng)的溫度分布，因此，以單個(gè)處置庫為研究對(duì)象．根據(jù)瑞典KBS-3V概念模型建立單個(gè)處置庫系統(tǒng)模型，系統(tǒng)包括封裝高放廢物固化體的廢物罐、緩沖材料、回填材料和圍巖，將緩沖材料分為Buffer1、Buffer2兩部分．為了保證計(jì)算的可執(zhí)行性，考慮單個(gè)處置庫的空間對(duì)稱性，使用二維軸對(duì)稱模型代替三維模型進(jìn)行分析．

為了避免邊界條件對(duì)計(jì)算結(jié)果的影響，取計(jì)算模型尺寸為20 m×1 000 m，將坐標(biāo)原點(diǎn)(0，0)設(shè)置在廢物罐中心，如圖3所示．圖4是計(jì)算模型的網(wǎng)格劃分，為了精確分析緩沖、回填區(qū)的溫度分布情況，故靠近廢物罐的部分網(wǎng)格劃分更加密集．完整網(wǎng)格包含6 985個(gè)單元和7 264個(gè)節(jié)點(diǎn)．

圖3 雙層緩沖層處置庫幾何模型 圖4 圖3對(duì)應(yīng)幾何模型的網(wǎng)格劃分

3 模擬與討論

3.1 模型驗(yàn)證

Hodgkinson[16]針對(duì)隨時(shí)間指數(shù)衰減的線性熱源，建立實(shí)體圓柱型模型，給出了對(duì)應(yīng)的半解析解，模型中各點(diǎn)溫度隨時(shí)間變化的表達(dá)式如下：

(15)

其中：R，z分別為徑向坐標(biāo)和軸向坐標(biāo)；t為時(shí)間；b為圓柱體的半長；r為其半徑；Q為熱源衰變熱，見式(8)；q0為初始熱源衰變熱，由式(8)可知q0=6.353 kW/m3；erf(x)為誤差函數(shù)，表示為：

(16)

I0(x)為修正的貝塞爾函數(shù)[17]，表示為：

(17)

為了驗(yàn)證本文模型方法的適用性，使用COMSOL Multiphysics 5.3a構(gòu)建如圖5左側(cè)所示的二維軸對(duì)稱模型，根據(jù)表2參數(shù)進(jìn)行了熱力學(xué)數(shù)值模擬分析，將結(jié)果與解析解進(jìn)行了比較．

表2 驗(yàn)證模型相關(guān)參數(shù)

兩種方法計(jì)算的緩沖層在金屬罐中心的溫度變化對(duì)比如圖5所示．

圖5 Hodgkinson解析解問題[16]對(duì)應(yīng)數(shù)值模型與結(jié)果驗(yàn)證

解析解和數(shù)值分析計(jì)算得到的峰值溫度分別為71.77℃和72.64℃．峰值溫度的相對(duì)誤差為0.71%，且兩種方法得到達(dá)到峰值溫度均需要兩年時(shí)間．整個(gè)過程中的平均誤差為7.13%．由此可知，在該模型中，使用COMSOL Multiphysics 5.3a進(jìn)行熱力學(xué)分析得到的數(shù)值解與解析解較吻合，因此，可以使用COMSOL Multiphysics 5.3a采用上述物理方程對(duì)高放處置庫的熱力學(xué)行為進(jìn)行分析．

3.2 參數(shù)敏感性分析

(16)

參數(shù)敏感性函數(shù)可以比較系統(tǒng)特性對(duì)各參數(shù)的敏感程度，敏感性函數(shù)值越大，系統(tǒng)對(duì)該參數(shù)越敏感．每個(gè)基準(zhǔn)數(shù)據(jù)對(duì)應(yīng)6個(gè)分析數(shù)據(jù)，各參數(shù)的分析數(shù)據(jù)一般在參考數(shù)據(jù)的0.5～1.5倍范圍內(nèi)，個(gè)別參數(shù)在此范圍上有所縮放，以求涵蓋該參數(shù)的合理范圍．

熱力學(xué)計(jì)算中，采用常規(guī)單層緩沖材料進(jìn)行參數(shù)敏感性分析，兩個(gè)主要參數(shù)分別為熱傳導(dǎo)系數(shù)λ，比定壓熱容cp，故設(shè)定基準(zhǔn)狀態(tài)X*=(λ，cp)，系統(tǒng)特征F*=f(λ，cp)．由談云志，彭帆，錢芳紅，等[19]及陳航，張虎元,郭永強(qiáng)，等[20]的試驗(yàn)及預(yù)測，緩沖材料的熱傳導(dǎo)系數(shù)在0.4～1.5 W/(m·K)的范圍內(nèi)，比定壓熱容在500～1 700 J/(kg·K)的范圍內(nèi)，故取基準(zhǔn)數(shù)據(jù)λ=1.4 W/(m·K)、cp=1 000 J/(kg·K)，計(jì)算得到的處置庫峰值溫度作為基準(zhǔn)值，分析處置庫峰值溫度對(duì)上述兩個(gè)參數(shù)的敏感程度，結(jié)果見表3．

表3 敏感性參數(shù)計(jì)算表

由表3可計(jì)算熱傳導(dǎo)系數(shù)和比定壓熱容的平均敏感性函數(shù)值分別為24.83和0.11．熱傳導(dǎo)系數(shù)的敏感性函數(shù)值遠(yuǎn)大于比定壓熱容敏感性函數(shù)值．

由圖6的曲線斜率大小可以判斷，材料比熱容變化對(duì)峰值溫度幾乎沒有影響；熱傳導(dǎo)系數(shù)變化對(duì)峰值溫度的影響較大．

圖6 峰值溫度-參數(shù)變化率關(guān)系

但隨著導(dǎo)熱系數(shù)的增加，效率增益逐漸減小，這是因?yàn)榫彌_材料的導(dǎo)熱系數(shù)超過處置庫周圍的圍巖后，緩沖層和圍巖整體形成了熱的不良導(dǎo)體，廢物罐表面熱量雖然能夠快速通過緩沖層到達(dá)與圍巖的交界面，但是由于圍巖的導(dǎo)熱系數(shù)較小，熱量不能及時(shí)導(dǎo)出，導(dǎo)致溫度變化較?。虼耍趯?duì)比單層和雙層緩沖區(qū)對(duì)峰值溫度的影響時(shí)可忽略比熱容的影響，以簡化計(jì)算．

3.3 溫度分布及對(duì)比

根據(jù)第2節(jié)所述單層、雙層緩沖層數(shù)值模型，利用有限單元計(jì)算程序COMSOL Multiphysics 5.3a計(jì)算了100 a內(nèi)模型溫度演變．

圖7為填充單層緩沖材料的處置庫在運(yùn)行2 a后，處置庫坑道周圍溫度的分布情況．峰值溫度出現(xiàn)在廢物罐表面的中心位置，并通過廢物罐、緩沖區(qū)、回填區(qū)向圍巖擴(kuò)散．離廢物罐越近，溫度越高．且回填區(qū)的溫度明顯低于緩沖區(qū)．

圖7 單層緩沖材料處置庫運(yùn)行2年后溫度分布

圖8為填充單層緩沖材料的處置庫自開始運(yùn)行后，廢物罐表面中心位置(如圖中點(diǎn)A)溫度隨時(shí)間的演化情況．廢物罐表面A點(diǎn)溫度在處置庫開始運(yùn)行1.8 a后即達(dá)到峰值溫度110.65℃，然后由于衰變熱的降低，隨運(yùn)行時(shí)間而逐漸降低，處置庫運(yùn)行100 a后溫度降至57℃．

圖8 單層緩沖材料處置庫廢物罐表面A點(diǎn)溫度時(shí)間演化

圖9(a)、(b)分別為單層緩沖材料與雙層緩沖材料對(duì)應(yīng)處置庫緩沖區(qū)3個(gè)特征點(diǎn)(點(diǎn)A、B、C)的溫度隨時(shí)間的演化過程．由圖可知，兩種結(jié)構(gòu)處置庫均在運(yùn)行1.5 a左右在廢物罐表面中心位置(點(diǎn)A)達(dá)到峰值溫度．雙層緩沖材料處置庫峰值溫度為92.73℃，單層緩沖材料處置庫峰值溫度為110.65℃，可見，將處置庫緩沖層進(jìn)行分區(qū)設(shè)計(jì)可以明顯降低處置庫峰值溫度，溫度總體下降了16%；同時(shí)，隨著時(shí)間變化兩種結(jié)構(gòu)A點(diǎn)溫度趨于接近，處置庫運(yùn)行100 a后均趨于60℃．由此可見，雙層緩沖區(qū)結(jié)構(gòu)可以顯著降低處置庫溫度最不利時(shí)段的運(yùn)行環(huán)境溫度，對(duì)處置庫安全穩(wěn)定運(yùn)行提供了有效保障．

圖9 緩沖區(qū)監(jiān)測點(diǎn)溫度變化

4 結(jié) 論

本研究采用熱力學(xué)分析方法，使用COMSOL Multiphysics 5.3a有限元軟件，對(duì)瑞典KBS-3V處置庫概念模型進(jìn)行相應(yīng)分析，確定了熱力學(xué)敏感性參數(shù)，并通過對(duì)比單、雙層緩沖層處置庫設(shè)計(jì)方案對(duì)處置庫中峰值溫度的影響，得到了如下結(jié)果：

1)材料比熱容變化對(duì)峰值溫度幾乎沒有影響；熱傳導(dǎo)系數(shù)變化對(duì)峰值溫度的影響較大，但隨著導(dǎo)熱系數(shù)的增加，效率增益逐漸減小．

2)相對(duì)單層緩沖層結(jié)構(gòu)，采用雙層緩沖層設(shè)計(jì)可以將處置庫峰值溫度降低約16%，可以使處置庫在臨界溫度(100℃)以下安全運(yùn)行．