KBS－3V高放廢物處置概念庫溫度場的解析求解方法簡介及其啟示

潘小青， 李玉曉， 魏望和

(江西理工大學理學院，江西 贛州341000)

KBS－3V高放廢物處置概念庫溫度場的解析求解方法簡介及其啟示

潘小青， 李玉曉， 魏望和

(江西理工大學理學院，江西 贛州341000)

文章簡要介紹Probert和Claesson針對瑞典KBS-3V高放廢物處置概念庫開展的溫度場解析求解方法研究，包括全域解和局域解.全域解將處置庫區域內規則排列在網格單元上的廢物罐熱源簡化為無限大平面熱源得到，局域解由無限多平行巷道排列的線熱源、無限大平面上網格點熱源，及沿廢物罐高度方向上的有限長線熱源產生.最終所有附加場疊加在處置庫原始溫度場得到總溫度場.應用解析法計算處置庫溫度場與用數值法計算溫度場相比，得到的表達式簡潔，物理意義明確，計算時間更快，適用的處置時間更長.最后，提出了開展中國高放廢物處置概念庫溫度場解析求解的研究思路.

高放廢物處置庫;溫度場;解析解;全域解;局域解

0 引 言

高放廢物的最終處置目前國際上通用的做法是深部地質處置，即將高放廢物埋藏于地下深部地質體中，通過構建工程屏障與地質屏障一起實現高放廢物與人類和生態環境的長久隔離[1].高放廢物的放熱導致處置庫內部溫度的升高，升高的幅度與廢物的放熱量、處置庫的結構形式、工程材料與圍巖的熱學性質有關，也與處置平巷的間距、平巷或鉆孔中廢物放置的密度有關，因此，處置庫溫度場分析是結構設計中的重要環節[2].求解溫度場的方法主要有解析法和數值法兩種，因熱源放熱模式和處置庫結構的復雜性，數值模擬一直是核廢物處置庫設計中的主要方法.但數值法只能計算處置后1000年內的情況[3]，而處置庫的運營和需要的評價時間通常要1萬年甚至百萬年.為此，Probert和Claesson對KBS-3V概念庫溫度場計算開展了解析方法研究，得出了簡化的解析公式[4]，公式適用處置庫關閉1萬年以內，這遠遠超出了數值計算的時間限制.解析方法可對處置庫內各種不同的物理過程和相互作用得到明確的物理圖像，并能獲得更簡潔的數學表達式以方便計算，同時能獲得對各種相互作用之間關系的更好理解[5].因此無論在時間和空間上，解析法都是數值法有益的補充.本文簡要介紹Probert和Claesson對KBS-3V概念庫的解析求解研究方法和研究結果及H?kmark和F?lth對簡化模型的修正，據此提出開展我國高放廢物處置概念庫溫度場解析法求解研究的思路.

1 高放廢物處置KBS－3V概念庫的數學模型

1.1 KBS－3V處置庫模型

瑞典核燃料與核廢物管理公司(SKB)提出的處置庫概念模型KBS-3V模式如圖1所示[6]，廢物體為退役的乏燃料，乏燃料由內襯鑄鐵的銅廢物罐裝載，暫存30～50年后再進行地下處置，處置庫為長2L寬2B的矩形區域.廢物罐埋放于地下H深度處，廢物罐高度Hc，直徑d，放于寬度D的平行巷道下方深度為h、直徑為d′的豎直孔洞中[7]，處置巷道間距為D′.在罐與巖石之間填充一定厚度的膨潤土為緩沖材料，在處置期廢物罐與膨潤土間存在間隙，膨潤土與巖石之間也有一定間隙.

圖1 KBS－3處置庫示意圖

1.2 KBS－3V處置庫熱傳導方程

以處置庫中心廢物罐為坐標原點建立直角坐標系，處于地下H深度處的廢物罐中心z=0，x方向為垂直巷道方向，-L＜x＜L.y方向沿巷道方向，-B＜y＜B.設處置庫原始地溫為T0(z)，T(x，y，z，t)為由釋熱造成的溫度增量，則熱傳導方程:

式(1)中λ為導熱系數(假定地下導熱系數均勻)，pc為體積熱容，a為熱擴散系數，a=λ/pc，Q(x，y，z，t)為廢物罐熱功率，廢物罐排列于矩形區域分割成的小網格中心.

2 Probert和Claesson對KBS－3V概念庫溫度場的解析方法及其修正

2.1 疊加原理

KBS-3V處置庫的溫度場看成原溫度場與由廢物罐放熱而產生的溫度場疊加而成，而處置庫溫度場是由一系列復雜的熱源產生的，由于方程是線性的可用疊加法分析.若設Q1和Q2分別表示總熱源中任兩個成分，分別引起的溫度增量為T1和T2，則由Q=Q1+Q2產生的溫度增量為T1+T2.

在時間上由于廢物罐的放射性衰變持續釋熱，根據Duhamel’s原理，設在t=0時所有熱源瞬時放熱一個單位得到的溫度解為Tinst(x，y，z，t)(坐標原點在中間廢物罐中心位置)，則由隨時間變化的熱源Q0(t)產生的溫度為:

2.2 全域解

為簡化問題，將處置庫的溫度增量看成空間上兩個部分的疊加:一部分是總體的處置庫規模，而另一部分為局部的處置單元規模，所以總解分為全域解Tgl和局域的處置單元規模的穩態解Tloc.全域解由所有廢物罐凈熱源釋放而產生，局域解則是處置庫內部局域引起的平衡熱源產生的雙周期變化的溫度場.

考慮廢物罐區域與處置庫大小相比很小，又將所有熱源等效為整個處置庫區域的均勻平面熱源q0(t)=Q0(t)/DD′產生的溫度場(如圖2)，另外由于地表不受熱源影響即邊界條件(2)，需在z=2H處疊加一個同樣熱功率的鏡像負熱源-q0(t)，由此產生的全域溫度場Tgl(x，y，z，t):

式(4)中erf為誤差函數.

圖2 等效均勻平面熱源

2.3 局域解

局域解由3對凈釋熱平衡的熱源產生，如圖3所示[7].

圖3 局域解平衡熱源

A1:無限長無限多線源，線源單位長度熱功率為q(t)=Q0(t)/D，線源間距D′.

A2:無窮大平面負熱源，單位面積熱功率-q0= -Q0(t)/DD′.

B1:無限長點源，點源熱功率為Q0(t)，點源間距D.

B2:無限長負線源，線源單位長度熱功率為-q(t)=-Q0(t)/D.

C1:高度為Hc的有限長線源，單位長度熱功率為qc=Q0(t)/Hc.

C2:負點源-Q0(t).

其中A1的中心線與B2、B1的中心點與C2兩對平衡熱源相互抵消，剩下的3個熱源組為:

C1的解:

無限長沿巷道的單一線源產生的溫度解:

沿巷道的無限長點源產生的溫度解:

式(7)中，K0為修正貝塞爾函數.A1減去中心線和B1減去中心點剩下的源產生的附加溫度與A2和C1疊加得到總的局域解表達式為:

2.4 總溫度場

根據疊加原理，KBS-3V處置庫溫度解的表達式為:

在KBS-3V中Q0(t)隨時間變化很慢，最短的時間常數為46年，達到穩態的時間只需幾年，所以可以當做恒定熱源得到穩態解，再用實際的Q0(t)得到緩變解.

上述簡化模型存在以下缺陷:其一假定處置庫空間導熱系數和比熱均為常數，實際在廢物罐的周邊不同材料具有不同的熱特性;其二報告設定熱源的半徑為零，實際上廢物罐有限定的半徑，而且熱量的釋放也不均勻;其三設邊界無限遠離熱源，實際邊界上有廢物罐釋熱;另外廢物罐和緩沖回填材料之間具有空氣間隙，緩沖回填材料和主巖之間具有水充填的間隙，這些在簡化熱源模型中都被忽略.

2.5 H?kmark和F?lth的修正

Harald H?kmark和Billy F?lth對上面的熱源簡化模型做了修正[8]，修正模型將廢物罐表面溫度看成三部分的疊加，即巖石原始溫度T0+硐壁溫度增量Trω+廢物罐與巖石之間的熱對流和熱傳輸特性決定的溫度偏差ΔT.在求解處置硐壁溫度增量Trω時采用了復合線源模型代替前面的C1，如圖4所示，即將中心廢物罐附近三個熱源看成標準高度的線熱源與一有效高度的負熱源的疊加，負熱源的有效高度通過原型處置庫實驗數據得到.

圖4 中心罐局域熱源修正模型

考慮廢物罐熱輸出模式由熱輻射和熱傳導兩種方式，另外考慮各種間隙(如廢物罐和緩沖材料之間、緩沖材料和圍巖之間等)的存在，溫度偏差ΔT為:

式(11)中R1、R2分別為廢物罐和處置硐室的半徑，λb(eff)為罐與巖石之間的有效導熱系數，q(t)為罐表面單位面積的熱功率.

3 Probert和Claesson解析方法的驗證及應用

3.1 解析結果與數值結果的對比

Claesson和Probert分別將數據代入解析法得到的溫度場結果并用MATLAB程序計算任一點的局域溫度，計算用時僅需3 s，同時討論了不同處置時間、不同方向的剖面的溫度分布圖及全域場從1000年到10000年的溫度演化規律[7].H?kmark和F?lth應用修正的解析表達式計算了KBS-3V豎直處置概念庫中廢物罐表面的最高溫度，并用解析法與數值法相結合計算了KBS-3H水平處置概念庫中廢物罐表面的最高溫度[8]，用解析法得到的廢物罐表面最高溫度與用有限差分計算程序FLAC2D得到的結果僅相差0.2°C，而用Claesson和Probert最初的線源模型計算結果則高出修正模型的溫度2.3°C.H?kmark和Claesson應用溫度場解析解計算了廢物罐表面中心附近的溫度[9]，與用三維離散元法(3DEC)[10]計算的結果在不同位置處基本一致，處置5000年后解析法與數值法計算的溫度僅相差0.5°C.

3.2 解析方法在處置庫設計中的應用

KBS-3V處置庫溫度場的解析結果給處置庫設計提供了思路.從方程中可知改變巷道間距和廢物罐間距對溫度的影響，比如保持單位面積處置罐數目不變的條件下，減小總巷道長度33%將造成處置庫最高溫度上升10°C，原因是增大了巖石的熱阻.而為保持最高溫度不超過上限，需要比基礎情況增加250%的巷道總長度.另外，如果將一個處置庫分解為分別位于地下500 m和600 m的兩個層疊的分處置庫，計算表明在處置20～300年間，溫度將比基礎情況增加10°C[9].解析法以巖石的理想熱傳輸性質為基礎能快速有效地得到處置庫巖石在全部時間和空間上的熱演化規律，可以作不同的處置庫布局和巖石熱性質假設得到數以千計的結果，可以演示處置時間分布不同時中心和邊緣緩沖材料的峰值溫度差別，因此在處置庫戰略設計布局中提供巷道間距和處置硐室間距的合理范圍. KBS-3V正是以解析解為基礎確定巷道間距為40 m，以解析解提供的處置室間距范圍通過進一步數值計算確定處置室間距[11].

4 中國開展處置概念庫研究的思路

4.1 我國處置概念庫的可能結構及與KBS－3V的比較

中國高放廢物(HLW)地質處置概念庫的處置對象是玻璃固化體、廢物罐的材料為低碳鋼，確定使用膨潤土作為處置庫的回填材料，圍巖為花崗巖，處置庫初步設計為豎井一坑道型[12].與KBS-3V相比，中國處置概念庫的結構和圍巖及緩沖材料選取與瑞典都相同，所不同的是廢物罐內是經過后處理的玻璃固化體，經過30年暫存后熱功率約559 W，為

KBS-3V乏燃料廢物罐熱功率的1/3.概念庫模型假定處置庫深度為-500 m，處置區域分為4個，每個區域有300對處置巷道，每對巷道處置80個廢物罐.處置庫中廢物罐分別有豎直放置和水平放置兩種方式.水平處置方式中廢物罐外有5 mm空氣間隙，空氣間隙外為500 mm膨潤土塊，膨潤土塊與主巖間有55 mm組合間隙.豎直處置方式巷道中線間距為9.5～10 m，沿巷道底部中線等間距開挖圓柱體處置室，廢物罐置于處置室下部中央，其上、下部和左右兩側放置緩沖材料處置巷道內放置回填材料，處置室中線間距由最高溫度限制進行數值計算[13-14].

4.2 中國處置概念庫溫度場解析法求解的途徑和難點

目前趙宏剛等根據現有的玻璃固化體剩余衰變熱、工程材料的熱物理參數、巖石的熱特性及各種間隙的熱性質等相關參數對我國高放廢物處置庫進行了溫度場數值計算[13-15]，劉文崗等在處理單個罐時用了簡單的軸對稱模型得到單個處置硐室圍巖的溫度解析表達式[16]，在計算全域溫度場時仍采用數值方法.中國HLW概念庫溫度場的解析求解需要對現有處置庫設計模型中各部分的幾何參數、熱源性質、工程材料和圍巖的熱參數等有可驗證的實驗數據，在作簡化時才有據可依，這是問題的難點.在KBS-3V溫度場解析求解過程中隨著最新數據不斷地修正公式，因此對我國處置庫為得到熱方程解析解，可能途徑如下:

1)疊加原理應用.考慮到熱傳導方程是線性方程，處置庫溫度場可看成不同形式熱源產生的溫度場疊加而成，具體疊加形式可采用H?kmark和F?lth的形式.

2)處置硐室內部的溫度偏差可參考已有的結果[16]，但需要用建立廢物罐與巖壁之間的有效熱導率公式.

3)巖壁的溫度增量求解要考慮全部處置庫空間二維周期性熱源引起的邊界條件特點.

4)時間上考慮熱源釋熱模式的時間演化規律和熱源熱功率的大小對達到穩態時間的影響等.

[1]羅嗣海，錢七虎，王駒.高放廢物地質處置庫的特點及其結構型式[J].地質科技情報，2007，26(5):83-90.

[2]羅嗣海，錢七虎，賴敏慧，等.高放廢物處置庫的工程設計[J].巖石力學與工程學報，2007，26(增2):3904-3911.

[3]Lars Ageskog，Patrik Jansson.Heat propagation in and around the deep repository-Thermal calculations applied to three hypothetical sites:Aberg，Beberg and Ceberg[R].SKB，TR-99-02，1999.

[4]Jhan Claesson，Thomas Probert.Temperature field due to timedependent heat sources in a large rectangular grid-Derivation of analytical solution[R].SKB，TR96-12，1996.

[5]Jhan Claesson，Thomas Probert.Thermoelastic stress due to rectangular heat source in a semi-infinite medium--Derivation of analytical solution[R].SKB，TR96-13，1996.

[6]SR 97-Waste，repositary design and sites[R].SKB，TR-99-08，1999.

[7]Thomas Probert，Jhan Claesson.Temperature field due to timedependent heat sources in a large rectangular grid-Application for the KBS-3 repository[R].SKB，TR97-27，1997.

[8]Harald H?kmark，Billy F?lth.Thermal dimensioning of the deep repository[R].SKB，TR-03-09，2003.

[9]Harald H?kmark，Johan Claesson.Use of an analytical solution for calculating temperatures in repository host rock[J].Engineering Geology，2005，81(3):353-364.

[10]Harald H?kmark.Acceptance of emplacement hole positions-stage 1 thermomechanical study[J].Engineering Geology，1998，49(2): 215-222.

[11]Harald H?kmark，Margareta L?nnqvist，Ola Kristensson.Strategy for thermal dimensioning of the final repository for spent nuclear fuel[R].SKB，R-09-04，2009.

[12]Ju Wang.High-level radioactive waste disposal in China:update 2010[J].Journal of Rock Mechanics and Geotechnical Engineering，2010，2(1):1-11.

[13]趙宏剛，王駒，蘇銳.HLW水平處置硐室間距分析[J].世界核地質科學，2014，31(增1):305-312.

[14]劉月妙，王駒，蔡美峰等.熱力藕合條件下高放廢物處置室間距研究[J].鈾礦地質，2009，25(6):378-379.

[15]趙宏剛，王駒，劉月妙，等.HLW豎直處置熱分析[J].世界核地質科學，2013，30(1):44-51.

[16]劉文崗，王駒，周宏偉，等.高放廢物處置庫花崗巖熱-力藕合模擬研究[J].巖石力學與工程學報，2009，28(增1):2875-2883.

Introduction and thoughts on the approach of analytical solution to temperature field in KBS－3V HLW concept repository

PAN Xiaoqing，LI Yuxiao，WEI Wanghe
(Faculty of Science，Jiangxi University of Science and Technology，Ganzhou 341000，China)

Method study of the analytical solution to the temperature field in KBS-3V HLW repository in Sweden by Probert and Claesson is introduced.The solution consists of global temperature field and the local one.The grid points within the repository are simplified as the heat source distribution to a rectangular plane heat source for global solution.Local solution consists of a line heat source，an infinite rectangular grid of point sources and a balancing plane sink.The total temperature is the sum of the global temperature，the local temperature and the undisturbed rock temperature.The analytical approach gains simple formulas and physical understanding.It is valid to longer repository time and shorter calculating time compared with numerical approach.The approach of analytical solution to temperature field in HLW concept repository is proposed in the end.

HLW repository;temperature field;analytical solution;global solution;local solution

2095-3046(2015)01-0104-05

10.13265/j.cnki.jxlgdxxb.2015.01.018

TL942

A

2014-09-05

江西省教育廳科技項目(GJJ11467)

潘小青(1966- )，女，教授，主要從事計算物理、力學等方面的研究，E-mail:panxqecit@163.com.