高放废物处置库二维轴对称温度场模型构建及应用研究

周祥运; 孙德安; 许迅; 彭帆

doi:10.11779/CJGE20230253

高放废物处置库二维轴对称温度场模型构建及应用研究 English Version

周祥运^1,,
孙德安^2, ,,
许迅²,
彭帆³

1.
南京工程学院土木工程与智慧管理研究所, 江苏南京 211167
2.
上海大学力学与工程科学学院, 上海 200444
3.
衢州学院建筑工程学院, 浙江衢州 324000

基金项目:

国家自然科学基金项目 42077229

江苏省高等学校基础科学（自然科学）研究面上项目 23KJD560002

南京工程学院引进人才科研启动基金项目 YKJ202323

详细信息

作者简介:
周祥运(1992—)，男，博士，讲师，主要从事岩土工程的科研工作。E-mail: zhouxiangyun@njit.edu.cn

通讯作者:
孙德安, E-mail: sundean@shu.edu.cn

中图分类号: TU43
计量
- 文章访问数: 311
- HTML全文浏览量: 46
- PDF下载量: 81
出版历程
- 收稿日期: 2023-03-22
- 网络出版日期: 2023-09-21
- 刊出日期: 2024-10-31

Establishment and application of two-dimensional axisymmetric temperature field model for high-level radioactive waste disposal repository

1.
Institute of Civil Engineering and Intelligent Management, Nanjing Institute of Technology, Nanjing 211167, China
2.
School of Mechanics and Engineering Science, Shanghai University, Shanghai 200444, China
3.
College of Civil Engineering and Architecture, Quzhou University, Quzhou 324000, China

摘要

摘要: 高放废物处置库温度场长期演化特征是处置库缓冲层安全性能评估和设计、处置容器空间优化布局的重要依据之一。建立了处置单元二维轴对称双层热分析模型，以揭示处置库多重屏障系统近场温度演化规律。对传热控制方程应用拉普拉斯和有限傅里叶正弦变换，得到缓冲层和围岩层温度的拉普拉斯域解，借助Crump方法对拉普拉斯域解进行数值反演获得对应的时域解；通过与线热源解和数值解的对比验证了模型的正确性，分析了处置单元近场多屏障系统温度场的时空分布特征，探讨了不同参数对缓冲层峰值温度的影响；借助模型半解析解确定最小处置容器间距以及对原位试验结果进行预测。结果表明：①与线热源解相比，模型半解析解可以更准确地模拟缓冲层温度变化；②当处置隧道间距超过40 m时，其对缓冲层峰值温度影响较小；③在处置隧道间距取40 m时，最小处置容器间距为7.7 m；④模型半解析解可以较好地预测原位试验结果。
- 高放废物处置库 /
- 多屏障系统 /
- 热分析模型 /
- 峰值温度 /
- 处置容器间距
Abstract: The long-term evolution of temperature field in high-level radioactive waste disposal repository is one of the important bases for the safety evaluation and design of buffer layer and the spatial optimization of disposal container. A two-dimensional axisymmetric thermal analysis model with two layers of disposal unit is established to reveal the near-field temperature evolution of multi-barrier system. The Laplace and finite Fourier sine transforms are applied to the governing equations for heat transfer to obtain the temperature solutions to the buffer layer and the surrounding rock layer in the Laplace domain. The corresponding solutions in time domain are obtained by numerical inversion of the Laplace domain solutions with the help of the Crump method. The correctness of the model is verified by comparing with the linear heat source solution and the numerical solution. The temporal and spatial distribution characteristics of the temperature field in the near-field multi-barrier system of the disposal unit are analyzed, and the effects of different parameters on the peak temperature of the buffer layer are discussed. The semi-analytical solution of the model is used to determine the minimum disposal container spacing and predict the results of the in-situ tests. The results show that the semi-analytical solution of the model can simulate the temperature change of buffer layer more accurately by comparing with the linear heat source solution. When the tunnel spacing exceeds 40 m, it has small effects on the peak temperature of buffer layer. When the tunnel spacing is 40 m, the container spacing is 7.7 m. The semi-analytical solution of the model can well predict the results of the in-situ test.
- high-level radioactive waste disposal repository /
- multi-barrier system /
- thermal analysis model /
- peak temperature /
- disposal container spacing

HTML全文

0. 引言

高放射性核废物（简称高放废物）的安全处置是一个世界性的科学和技术难题，并已引起了广泛的社会关注^[1]。学者们几十年的研究探索表明，深层地质处置法是目前安全处置高放废物唯一技术可行的永久性工程解决方案。高放废物深层地质处置库一般采用的是“多重屏障系统”设计，即把废物储存在处置容器中，外面包裹缓冲材料，再向外为围岩。处置容器中核废物释放的衰变热会传递到缓冲层和围岩，造成近场屏障系统温度不断升高^[2-3]。过高的温度会影响地下水的流动性、缓冲层的膨胀性、渗透性和传热特性、以及放射性核素的吸附和扩散等^[4-5]。在处置库的深部环境中，缓冲层长期遭受复杂的热-水-力-气多场耦合作用。针对缓冲材料多场耦合特性，陈正汉等^[6-8]在仪器研发、试验研究、理论模型、轴对称问题的有限元表达式推导及软件设计、数值计算等方面做了大量工作。而在进行处置库处置容器间距设计时，处置库温度场分布特征是主要的依据之一，这对设计而言是偏于安全的。处置库的典型温度设计准则是缓冲区的峰值温度在处置库运营期间应保持在100 ℃以下^[9]。温度设计准则对处置容器间距设计是一种限制，因为它制约了在处置库的给定区域内处置容器的数量，即处置密度。处置容器间距和处置隧道间距应根据处置库热分析结果和温度设计准则来确定^[10]。因此，温度场的演变特征对处置库的安全性能评估和间距设计具有重要意义。

现有的处置库温度场解析方法都是基于Carslaw等^[11]提出的点热源作用下无限大岩体的温度增量模型。针对距离处置区一定距离的点（1倍~2倍隧道间距之间），Claesson等^[12]建立了基于矩形热源的解析解，分析了瑞典处置库的温度演化规律。在短距离内，基于点热源和矩形热源的解析表达式不能很好地描述圆柱形热源的温度分布^[13]。对此，Hökmark等^[14]由点热源解析式对处置容器的整个高度积分得到了线热源解析解，然后考虑到处置容器顶部、底部和中部热流量的差异，建立了复合线热源解析模型，并借助这些解析解分析了瑞典和芬兰两国处置库的温度场长期演化特征。刘东东等^[15]建立了处置容器分布的线热源解析模型，研究了不同处置隧道间距和处置容器间距组合条件下缓冲层峰值温度变化情况。然而，在以往的研究中，处置库围岩中任意一点的温度分布都没有统一的解析表达式。线热源模型通常把处置容器的半径设置为零，这与实际情况不相符，当处置容器长径比较小时，线热源解将会产生一定误差。此外，在建立处置单元热分析模型时，以往的解析模型很少考虑缓冲层与围岩之间传热特性的差异。

尽管数值方法可以快速获取处置库温度场分布特征，但在进行一些参数分析时，解析方法会更加方便，同时也可为数值方法提供验证手段。因此，首先建立了二维轴对称条件下处置单元的双层传热模型，对传热控制方程应用拉普拉斯和有限傅里叶正弦变换得到缓冲层和围岩层温度的拉普拉斯域解，通过与线热源解和数值解的比较，验证了本文解的正确性。然后，分析了不同几何参数和材料热参数对缓冲层峰值温度的影响。最后，借助模型半解析解确定了处置容器的最小间距，并对原位试验结果进行了预测。

1. 计算模型及基本假定

1.1 模型计算简图

以KBS-3V型（由瑞典开发的高放废物多屏障深地质竖直处置方案）处置库概念设计为例^{[14, 16]}，处置单元多重屏障系统传热计算模型如图 1所示。模型底部到顶部的距离为L，处置容器顶部和底部到地面的距离分别为z₁，z₂，处置容器高度为h，处置容器、膨润土块层和膨润土颗粒层的外半径分别为a，b，c。整个模型区域被限定在0 < z < L，0 < r < +∞的无限圆柱内。

图 1 处置单元多重屏障系统计算简图

Figure 1. Model for multi-barrier system in disposal unit

下载: 全尺寸图片幻灯片

在建立方程时作了如下假定：①由于缓冲材料和围岩的低渗透性，热量传递的方式主要为热传导，热辐射与热对流产生的热量忽略不计^[17]。②在进行处置容器间距设计时，要考虑到最差的工程状况。处置库运营初期，缓冲层和围岩处于非饱和或接近干燥状态，随着地下水缓慢入渗和温度的升高，两者导热系数会逐渐增大。因此，出于设计安全考虑，本文模型对应处置库刚刚封场时的状态，不考虑地下水和温度对缓冲材料和围岩热传导系数的影响^{[14, 16]}。③膨润土颗粒层主要用来填充膨润土块层和围岩层之间的施工缝隙，厚度大约在10~30 mm。因膨润土颗粒层厚度较薄，把膨润土块层和颗粒层当作1层处理^[18]。④Xu等^[19]通过COMSOL数值模拟发现，处置容器底部和顶部的膨润土块对缓冲层峰值温度影响很小。为了简化数学模型，忽略处置容器顶部和底部的膨润土。

1.2 传热控制方程

在直角坐标系下，缓冲层和围岩层的传热控制方程为^[20]

（ ） ，

$\begin{array}{c} \frac{{{\partial ^2}{T_{\text{B}}}(x,y,z,t)}}{{\partial {x^2}}} + \frac{{{\partial ^2}{T_{\text{B}}}(x,y,z,t)}}{{\partial {y^2}}} + \frac{{{\partial ^2}{T_{\text{B}}}(x,y,z,t)}}{{\partial {z^2}}}\\ =\frac{1}{{\chi }_{\text{B}}}\cdot \frac{\partial {T}_{\text{B}}(x,y,z,t)}{\partial t}（a < x < b,a < y < b,0 < z < L）\text{，} \end{array}$

(1)

（ ） 。

$\begin{array}{c} \frac{{{\partial ^2}{T_{\text{R}}}(x,y,z,t)}}{{\partial {x^2}}} + \frac{{{\partial ^2}{T_{\text{R}}}(x,y,z,t)}}{{\partial {y^2}}} + \frac{{{\partial ^2}{T_{\text{R}}}(x,y,z,t)}}{{\partial {z^2}}}\\ =\frac{1}{{\chi }_{\text{R}}}\cdot \frac{\partial {T}_{\text{R}}(x,y,z,t)}{\partial t}（x > b,y > b,0 < z < L）。 \end{array}$

(2)

由于模型是二维轴对称结构，为了简化计算，采用圆柱坐标表示，缓冲层和围岩层的传热控制方程为

（ ） ，

$\begin{array}{c} \frac{{{\partial ^2}{T_{\text{B}}}(r,z,t)}}{{\partial {r^2}}} + \frac{1}{r} \cdot \frac{{\partial {T_{\text{B}}}(r,z,t)}}{{\partial r}} + \frac{{{\partial ^2}{T_{\text{B}}}(r,z,t)}}{{\partial {z^2}}}\\ =\frac{1}{{\chi }_{\text{B}}}\cdot \frac{\partial {T}_{\text{B}}(r,z,t)}{\partial t}\text{ }（a < r < b,0 < z < L） \text{，} \end{array}$

(3)

（ ） 。

$\begin{array}{c} \frac{{{\partial ^2}{T_{\text{R}}}(r,z,t)}}{{\partial {r^2}}} + \frac{1}{r} \cdot \frac{{\partial {T_{\text{R}}}(r,z,t)}}{{\partial r}} + \frac{{{\partial ^2}{T_{\text{R}}}(r,z,t)}}{{\partial {z^2}}}\\ =\frac{1}{{\chi }_{\text{R}}}\cdot \frac{\partial {T}_{\text{R}}(r,z,t)}{\partial t}\text{ }（b < r,0 < z < L）。 \end{array}$

(4)

式中：下标B，R分别代表膨润土块层和围岩层；T为温度（℃）；χ为热扩散系数（m²/s），χ= $\lambda$ /(ρc)， $\lambda$ 为热传导系数（W/(K·m)）；ρ为密度（kg/m³）；c为比热容（J/(K·kg)），r为到处置容器轴线的半径距离（m）；z为到地表的垂直距离（m）；t为处置时间（a）。

1.3 初始及边界条件

处置库不同深度位置初始温度条件可以表示为

。

${T_{\text{B}}}(r,z,0) = {T_{\text{R}}}(r,z,0) = {t_0} + m \cdot z 。$

(5)

式中：t₀为地表的初始统一温度（℃）；m为处置库竖直方向的温度梯度（℃/m）。

在半径r=a处的热流量边界条件可以表示为

$\frac{{\partial {T_{\text{B}}}(a,z,t)}}{{\partial r}} = \frac{{ - \phi \cdot P(t)}}{{{\lambda _{\text{B}}} \cdot 2{\rm{ \mathsf{ π} }} a \cdot h}} \cdot [U(z - {z_1}) - U(z - {z_2})] 。$

(6)

式中：λ_B为缓冲材料的热传导系数(W/(K·m))；h为处置容器的高度(m)；z₁，z₂分别为处置容器顶部和底部的垂直坐标(m)；U(·)为单位阶越函数；P(t)为处置容器的衰变热功率（W）；ϕ为处置容器的侧表面面积与总表面积之比。

He等^[21]通过COMSOL数值模拟发现，处置容器上下表面的热流量对缓冲层峰值温度的影响较小。为了简化数学模型，只考虑侧面的热流量，每个处置容器初始衰变热功率修正值为ϕ·P(t)。处置容器衰变热功率P(t)的变化可借助ORIGEN软件模拟获得^[16]，其变化过程可通过指数叠加函数来拟合且精度满足工程需要。处置容器衰变热功率的指数叠加函数可以表示为

$P(t) = P(0) \cdot \sum\limits_{i = 1}^N {{a_i}} \exp ( - t/{t_i}) 。$

(7)

式中：P(0)为处置容器安装时的初始热功率；t为处置时间；t_i为时间常数，可以在20~20000 a之间任意选择，本文取20，50，200，500，2000，5000，20000；N为所选取的时间项数，N取值越大，拟合效果越好，一般N取5项，就可确保较好的精度，本文N取7项。a_i为与不同冷却时间相对应的系数，通过下式进行确定，可见满足a₁+a₂+a₃+a₄+a₅+a₆+a₇=1^[16]：

$\left\{ {\begin{array}{*{20}{c}} {P(0)} \\ {P({t_1})} \\ {P({t_2})} \\ \vdots \end{array}} \right\} = P(0) \cdot \left[ {\begin{array}{*{20}{c}} 1&1&1& \cdots \\ {{\text{e}^{ - \frac{{{t_1}}}{{{t_1}}}}}}&{{\text{e}^{ - \frac{{{t_1}}}{{{t_2}}}}}}&{{\text{e}^{ - \frac{{{t_1}}}{{{t_3}}}}}}& \cdots \\ {{\text{e}^{ - \frac{{{t_2}}}{{{t_1}}}}}}&{{\text{e}^{ - \frac{{{t_2}}}{{{t_2}}}}}}&{{\text{e}^{ - \frac{{{t_2}}}{{{t_3}}}}}}& \cdots \\ \cdots & \cdots & \cdots & \cdots \end{array}} \right] \cdot \left\{ {\begin{array}{*{20}{c}} {{a_1}} \\ {{a_2}} \\ {{a_3}} \\ \vdots \end{array}} \right\} 。$

(8)

相同温度和热流量作为缓冲层和围岩交界面的连续性条件，在半径r=b处的热边界条件可以表示为

${T_{\text{B}}}(b,z,t) = {T_{\text{R}}}(b,z,t) \text{，}$

(9)

${\lambda _{\text{B}}} \cdot \frac{{\partial {T_{\text{B}}}(b,z,t)}}{{\partial r}} = {\lambda _{\text{R}}} \cdot \frac{{\partial {T_{\text{R}}}(b,z,t)}}{{\partial r}} 。$

(10)

在半径无穷远处，温度可认为是常数且和相同深度处置库位置温度相同。在半径r→+∞时，热边界条件可以表示为

${T_{\text{R}}}( + \infty ,z,t) = {t_0} + m \cdot z 。$

(11)

在模型的上下边界处温度为常数，热边界条件可以表示为

${T_{\text{B}}}(r,0,t) = {T_{\text{R}}}(r,0,t) = {t_0} \text{，}$

(12)

${T_{\text{B}}}(r,L,t) = {T_{\text{R}}}(r,L,t) = {t_1} 。$

(13)

式中：t₁为地下深度L处岩体温度(℃)。

1.4 半解析解的推导

令

${T'_{\text{B}}}(r,z,t) = {T_{\text{B}}}(r,z,t) - ({t_0} + m \cdot z) \text{，}$

(14)

${T'_{\text{R}}}(r,z,t) = {T_{\text{R}}}(r,z,t) - ({t_0} + m \cdot z) 。$

(15)

对控制方程（3），（4）分别作关于变量t的拉普拉斯变换和关于变量z的有限傅里叶正弦变换可得

$\begin{array}{l} \frac{{{\text{d}^2}{{\tilde T'}_{\text{B}}}(r,w,s)}}{{\text{d}{r^2}}} + \frac{1}{r} \cdot \frac{{\text{d}{{\tilde T'}_{\text{B}}}(r,w,s)}}{{\text{d}r}} -\\ \left({w}^{2}+\frac{s}{{\chi }_{\text{B}}}\right)\cdot {\tilde{{T}^{\prime }}}_{\text{B}}(r,w,s)=0\text{ }（a < r < b） \text{，} \end{array}$

(16)

$\begin{array}{l} \frac{{{\text{d}^2}{{\tilde T'}_{\text{R}}}(r,w,s)}}{{\text{d}{r^2}}} + \frac{1}{r} \cdot \frac{{\text{d}{{\tilde T'}_{\text{R}}}(r,w,s)}}{{\text{d}r}} -\\ \left({w}^{2}+\frac{s}{{\chi }_{\text{R}}}\right)\cdot {\tilde{{T}^{\prime }}}_{\text{R}}(r,w,s)=0\text{ }（b < r）。 \end{array}$

(17)

式中：s为拉普拉斯变量，w为垂直坐标z的有限傅里叶正弦变换变量， ${\tilde T'_{\text{B}}}(r,w,s)$ ， ${\tilde T'_{\text{R}}}(r,w,s)$ 分别为变换后的缓冲层与围岩温度。

热边界条件经过拉普拉斯变换和有限傅里叶正弦变换后，可以表示为

${\tilde T'_{\text{R}}}( + \infty ,w,s) = 0 \text{，}$

(18)

$\frac{{\text{d}{{\tilde T'}_{\text{B}}}(a,w,s)}}{{\text{d}r}} = \frac{{ - Q(s)}}{{{\lambda _{\text{B}}} \cdot 2{\rm{ \mathsf{ π} }} a \cdot h}} \cdot \left[ {\frac{{\cos (w{z_1}) - \cos (w{z_2})}}{w}} \right] \text{，}$

(19)

${\tilde T'_{\text{B}}}(b,w,s) = {\tilde T'_{\text{R}}}(b,w,s) \text{，}$

(20)

${\lambda _{\text{B}}} \cdot \frac{{\text{d}{{\tilde T'}_{\text{B}}}(b,w,s)}}{{\text{d}r}} = {\lambda _{\text{R}}} \cdot \frac{{\text{d}{{\tilde T'}_{\text{R}}}(b,w,s)}}{{\text{d}r}} 。$

(21)

式中：

$Q(s) = P(0) \cdot \sum\limits_{i = 1}^7 {{a_i}} {t_i}/(s{t_i} + 1) 。$

(22)

式（16），（17）的通解可以表示为

$\left.\begin{array}{l}{\tilde{{T}^{\prime }}}_{\text{B}}(r,w,s)={C}_{1}\cdot {\text{I}}_{0}({q}_{\text{B}}r)+{C}_{2}\cdot {\text{K}}_{0}({q}_{\text{B}}r)\text{，}\\ {q}_{\text{B}}^{2}={w}^{2}+\frac{s}{{\chi }_{\text{B}}}\text{，}\end{array}\right\}$

(23)

$\left.\begin{array}{l}{\tilde{{T}^{\prime }}}_{\text{R}}(r,w,s)={C}_{3}\cdot {\text{I}}_{0}({q}_{\text{R}}r)+{C}_{4}\cdot {\text{K}}_{0}({q}_{\text{R}}r)\text{，}\\ {q}_{\text{R}}^{2}={w}^{2}+\frac{s}{{\chi }_{\text{R}}}。\end{array}\right\}$

(24)

式中：C₁，C₂，C₃，C₄为常数，可由边界条件确定。把式（23），（24）代入式（18）~（21），可得各系数为

$\left.\begin{array}{l}{C}_{1}=\mathit{\Delta } \cdot \frac{\Delta {f}_{1}}{\Delta f}\text{，}\\ {C}_{2}=\mathit{\Delta } \cdot \frac{\Delta {f}_{2}}{\Delta f}\text{，}\\ {C}_{3}=0\text{，}\\ {C}_{4}=\mathit{\Delta } \cdot \frac{\Delta {f}_{3}}{\Delta f}。\end{array}\right\}$

(25)

式中：

$\mathit{\Delta } = Q(s)(\cos (w{z_1}) - \cos (w{z_2}))/({\lambda _{\text{B}}}2{\rm{ \mathsf{ π} }} ahw) \text{；}$

(26)

$\Delta f = {\lambda _{\text{B}}}q_{\text{B}}^2{{\text{K}}_0}({q_{\text{R}}}b){n_1} - {\lambda _{\text{R}}}{q_{\text{B}}}{q_{\text{R}}}{{\text{K}}_1}({q_{\text{R}}}b){n_2} \text{；}$

(27)

$\Delta {f_1} = {\lambda _{\text{B}}}{q_{\text{B}}}{{\text{K}}_1}({q_{\text{B}}}b){{\text{K}}_0}({q_{\text{R}}}b) - {\lambda _{\text{R}}}{q_{{R}}}{{\text{K}}_0}({q_{\text{B}}}b){{\text{K}}_1}({q_{\text{R}}}b) \text{；}$

(28)

$\Delta {f_2} = {\lambda _{\text{B}}}{q_{\text{B}}}{{\text{K}}_0}({q_{\text{R}}}b){{\text{I}}_1}({q_{\text{B}}}b) + {\lambda _{\text{R}}}{q_{\text{R}}}{{\text{K}}_1}({q_{\text{R}}}b){{\text{I}}_0}({q_{\text{B}}}b) \text{；}$

(29)

$\Delta {f_3} = {\lambda _{\text{B}}}/b \text{；}$

(30)

${n_1} = {{\text{K}}_1}({q_{\text{B}}}b){{\text{I}}_1}({q_{\text{B}}}a) - {{\text{K}}_1}({q_{\text{B}}}a){{\text{I}}_1}({q_{\text{B}}}b) \text{；}$

(31)

${n_2} = {{\text{K}}_1}({q_{\text{B}}}a){{\text{I}}_0}({q_{\text{B}}}b) + {{\text{K}}_0}({q_{\text{B}}}b){{\text{I}}_1}({q_{\text{B}}}a) 。$

(32)

对式（23），（24）进行有限傅里叶正弦变换的逆变换，可得缓冲区和围岩区温度在拉普拉斯域内的解：

$\begin{array}{l} {\tilde T_{\text{B}}}(r,z,s) = \frac{2}{L} \cdot \sum\limits_{n = 1}^\infty {\left[ {\mathit{\Delta } \cdot \frac{{\Delta {f_1}}}{{\Delta f}} \cdot {{\text{I}}_0}({q_{\text{B}}}r) + \mathit{\Delta } \cdot \frac{{\Delta {f_2}}}{{\Delta f}} \cdot {{\text{K}}_0}({q_{\text{B}}}r)} \right]} \cdot \\ \sin (wz) + \frac{{{t_0} + 0.024z}}{s} \text{，} \end{array}$

(33)

$\begin{array}{l} {\tilde T_{\text{R}}}(r,z,s) = \frac{2}{L} \cdot \sum\limits_{n = 1}^\infty {\mathit{\Delta } \cdot \frac{{\Delta {f_3}}}{{\Delta f}} \cdot {{\text{K}}_0}({q_{\text{R}}}r)} \cdot \sin (wz) +\\ \frac{{{t_0} + 0.024z}}{s} 。 \end{array}$

(34)

式中： ${\tilde T'_{\text{B}}}(r,w,s)$ ， ${\tilde T_{\text{R}}}(r,z,s)$ 分别为缓冲区和围岩区在拉普拉斯域内的温度；r为半径； ${\lambda _{\text{B}}}$ ， ${\lambda _{\text{R}}}$ 分别为缓冲材料和围岩的导热系数；K₀(·)，I₀(·)分别为第一类和第二类零阶的修正贝塞尔函数；K₁(·)，I₁(·)分别为第一类和第二类一阶修正贝塞尔函数。

1.5 模型的验证

为了验证本文解的正确性，将本文理论解与数值模拟结果、线热源解析解结果进行对比分析。根据图 1所示的处置库概念设计，利用COMSOL软件建立的处置单元二维轴对称数值分析模型，如图 2所示。为了消除边界效应，处置容器距离模型上下边界的距离为500 m，径向围岩厚度设置为200 m。利用自由三角形网格将分析模型划分为3035个域单元和238个边界元。模型上、下边界初始温度设定为恒温，上边界温度设定为t₀ (℃)，下边界温度设定为t₁ (℃)，沿深度方向的温度梯度为m (℃/m)。热源以热流量的方式施加在缓冲层内边缘，热量在介质中的传热方式选择固体热传导。

图 2 处置单元多重屏障系统数值分析模型

Figure 2. Numerical model for multi-barrier system in disposal unit

下载: 全尺寸图片幻灯片

Ikonen^[16]和Hökmark等^[14]采用线热源解析解对芬兰和瑞典处置库温度场和处置容器间距进行研究，并通过CODE-BRIGHT数值解验证了线热源解析解在模拟围岩温度场方面的有效性。线热源解析表达式为^[14]

$\begin{array}{l} \Delta T(r,z,t) = \frac{1}{{{\rho _{\text{R}}} \cdot {c_{\text{R}}} \cdot 4{\rm{ \mathsf{ π} }} \cdot {\chi _{\text{R}}}}}\int\limits_0^t {\frac{{P(t')}}{{2h(t - t')}}} \cdot {\text{e}^{ - \frac{{{r^2}}}{{4{\chi _{\text{R}}}(t - t')}}}} \cdot \\ \frac{1}{2}\left[ {\text{erf}\left( {\frac{{h + {z_{\text{0}}}}}{{\sqrt {4{\chi _{\text{R}}}(t - t')} }}} \right) + \text{erf}\left( {\frac{{h - {z_{\text{0}}}}}{{\sqrt {4{\chi _{\text{R}}}(t - t')} }}} \right)} \right]\text{d}t' 。 \end{array}$

(35)

式中：ΔT(r, z, t)为围岩中任意一点的温度增量；P(·)为处置容器内部核废物的衰变热功率函数；h为处置容器的高度；r为距离处置容器轴线的半径距离；z₀为到处置容器中部的垂直距离；erf(·)为误差函数，t为处置时间；t' 为积分变量。

处置容器内部高放废物以瑞典基准燃料SVEA96为例^[13]，其燃烧值为38 MWd/kgU，冷却时间为40 a，初始热功率为1545 W。拉普拉斯数值反演的方法主要有Stehfest、Durbin、Crump方法等，其中Crump方法采用的ε算法可加快级数收敛速度。同时，与其他方法相比，该方法得到的数值解与精确解的误差更小^[22]。故本文采用Crump方法对本文半解析解进行数值反演，3种模型所需要的几何参数、材料热学参数以及核废物衰变参数取值如表 1所示^{[14, 16]}。图 3为3种模型解在深度500 m处岩壁位置（即缓冲层与围岩层交界面，r=0.825 m）和缓冲层内边缘位置（即缓冲层与处置容器交界面，r=0.525 m）的温度增量对比图。由图 3可以看出，在半径r=0.825 m处，本文半解析解、线热源解和COMSOL数值解计算得到的岩壁处峰值温度增量分别为30.45，31.12，29.63 ℃。本文半解析解与线热源解和COMSOL模拟计算结果吻合度较好，从而验证了本文半解析解的正确性。然而，在半径r=0.525 m处，线热源解明显低于本文半解析解和COMSOL数值解。这主要是由于线热源解忽略了膨润土缓冲层与围岩层热特性的差异所导致，可见线热源解并不能很好地模拟缓冲层温度的变化。

表 1 参数取值

Table 1. Values of parameters

参数	量值	参数	量值
λ_B/(W·(K·m)^-1)	0.68	z₂/m	502.625
λ_R/(W·(K·m)^-1)	2.50	t₀/℃	10
c_B/(J·(K·kg)^-1)	1014	t₁/℃	34
c_R/(J·(K·kg)^-1)	832	P(0)/W	1545
ρ_B/(kg·m^-3)	1 770	ϕ	0.91
ρ_R/(kg·m^-3)	2 630	a₁	0.049
a/m	0.525	a₂	0.696
b/m	0.825	a₃	-0.059
h/m	5.25	a₄	0.271
L/m	1000	a₅	0.027
z₀/m	0	a₆	-0.010
z₁/m	487.375	a₇	0.026

下载: 导出CSV

| 显示表格

图 3 本文半解析解与其他方法的对比

Figure 3. Comparison between semi-analytical method in this study and other methods

下载: 全尺寸图片幻灯片

2. 基于半解析解的处置单元近场温度分布

在KBS-3V型处置库概念设计中，一条主隧道一侧所有处置隧道的总和称为一个处置区^[23]，每个处置区的长度和宽度分别为1160，246 m。假设处置隧道间距为40 m，处置容器间距为6 m，则一个处置区可容纳29个隧道，每个隧道包含41个处置容器^[14]。一般来说，在近场温度分析中，只有位于有限半径内的处置容器才会对处置区中心处缓冲层的峰值温度产生一定影响。为了探究影响半径对缓冲层峰值温度的影响，考虑了5种不同的情况，即影响半径r_I为40，50，60，70，80 m。在40，50，60，70，80 m的影响半径范围内分别为35，39，51，61，75个处置容器。图 4分别给出了不同影响半径条件下缓冲层峰值温度随处置时间的变化情况。由图 4可以看出，影响半径对缓冲层峰值温度有显著影响。当影响半径为40，50，60，70，80 m时，峰值温度分别为出现在第6年的95.55℃，第9年的96.91℃，第10年的97.22℃，第10年的97.32℃以及第10年的97.37℃。当影响半径由60 m增加至70 m时，缓冲层峰值温度仅增加0.1℃。也就是说，当影响半径超过60 m时，其对峰值温度影响较小。为了确保计算的准确性，在下面温度场分析时影响半径取60 m。

图 4 不同影响半径下缓冲层温度变化

Figure 4. Temperature variation of buffer layer under different influence radii

下载: 全尺寸图片幻灯片

图 5给出了一个处置区中心位置处置单元多重屏障系统近场温度分布图。如图 5（a）所示，半径r为0.525，0.825 m处靠近处置容器的温度在处置前期急剧升高，之后由于核废物的热流量逐渐衰减，温度逐渐降低。处置容器表面峰值温度及其对应的时间取决于观测位置。半径r=0.825 m处的温度在处置前的20 a内仍呈上升趋势，这主要由于不同半径的处置容器核废物衰变热量热传导的延迟效应。如图 5（b）所示，在z=500 m时，处置容器中部温度最高、中部两侧的温度呈对称分布。同图 5（c）所示，处置后的第10年，缓冲层的温度梯度为79.86℃/m，围岩区温度梯度为7.58℃/m。缓冲层与围岩界面处存在较大的温度梯度。此外，相邻处置容器中心线两侧的温度分布几乎是对称的。

图 5 处置单元多重屏障系统近场温度分布图

Figure 5. Temperature distribution of multi-barrier system in disposal unit

下载: 全尺寸图片幻灯片

3. 参数分析

3.1 处置隧道间距

为了探讨处置隧道间距对缓冲层峰值温度的影响，这里考虑隧道间距20，30，40，50，60 m共5种情况。图 6给出了不同处置隧道间距（p_y）条件下缓冲层峰值温度（z=500 m）随处置时间的变化图。在处置隧道间距为20，30，40，50，60 m时，峰值温度分别为115.48，102.08，97.21，95.80，95.50℃。当处置隧道间距小于40 m，即20，30 m时，峰值温度随处置隧道间距的减小而显著上升。然而，当处置隧道间距大于40 m，即50，60 m时，其对峰值温度的影响很小。峰值温度对应的时间随着隧道间距的减小而向后移动。

图 6 处置隧道间距对缓冲层峰值温度的影响

Figure 6. Effects of tunnel spacing on peak temperature of buffer layer

下载: 全尺寸图片幻灯片

3.2 处置容器间距

图 7给出了在处置容器间距p_x为6，7，8，9，10 m时缓冲层峰值温度（z=500 m）随处置时间的变化图。由该图可以看出，不同处置容器间距对应的峰值温度分别为97.21，92.41，89.12，86.79，85.11℃。峰值温度曲线具有相似的变化规律，且峰值温度随处置容器间距的增大而减小，这主要是因为与缓冲材料（通常由压实膨润土制成）相比，具有更高导热性的围岩厚度增加。

图 7 处置容器间距对缓冲层峰值温度的影响

Figure 7. Effects of container spacing on peak temperature of buffer layer

下载: 全尺寸图片幻灯片

3.3 缓冲层厚度

缓冲层通常由压实的膨润土块组成，并需要一定的厚度以确保其优良的工程性能。图 8给出了不同缓冲层厚度d_b条件下缓冲层峰值温度（z=500 m）随处置时间的变化图。由图 8可以看出，峰值温度随着缓冲层厚度的增加而上升，缓冲层厚度0.1，0.2，0.3，0.4，0.5 m对应的峰值温度分别为86.29，92.08，97.21，101.70，105.85℃。缓冲层越厚，处置容器内的热流向外扩散越慢。因此，在保证缓冲层的膨胀性、渗透性和对核素的高吸附性等安全功能的基础上，应尽量减小缓冲层的厚度。

图 8 缓冲层厚度对缓冲层峰值温度的影响

Figure 8. Effects of thickness of buffer layer on its peak temperature

下载: 全尺寸图片幻灯片

3.4 缓冲材料热传导系数

压实高庙子膨润土因其渗透性低，对放射性核素的吸附能力强、力学性能优良，被中国选为缓冲/回填材料，但是与其他材料相比，其热传导系数相对较低。图 9给出了缓冲材料热传导系数λ_B为0.37，0.52，0.68，0.95，1.12 W/(K·m)时，缓冲层峰值温度随处置时间变化图。不同热传导系数对应的峰值温度分别为117.72，104.60，97.21，90.40，87.80℃。这主要是因为当缓冲材料热传导系数增大时，缓冲层的温度增量将会减小，这会导致峰值温度的下降。

图 9 缓冲层热传导系数对缓冲层峰值温度的影响

Figure 9. Effects of thermal conductivity of buffer layer on its peak temperature

下载: 全尺寸图片幻灯片

3.5 围岩热传导系数

图 10给出了不同围岩导热系数λ_R为2.1，2.3，2.5，2.7，2.9 W/(K·m)时缓冲层峰值温度（z=500 m）随处置时间的变化图。从图 10可以看出，围岩热传导系数对峰值温度有显著影响，峰值温度随围岩导热系数的增大而降低。5种情况下对应的峰值温度分别为104.74，100.67，97.21，94.22，91.60℃。围岩热传导系数的选择在温度场的计算中是个关键问题，可以通过室内和现场试验来综合确定。

图 10 围岩热传导系数对缓冲层峰值温度的影响

Figure 10. Effects of thermal conductivity of rock on peak temperature of buffer layer

下载: 全尺寸图片幻灯片

4. 温度场模型应用

4.1 确定处置容器最小间距

当缓冲层温度超过100℃时，缓冲材料的膨胀性能将受到恶化，此外，过高的温度将使地下水的盐分在处置容器表面发生富集，导致容器材料的腐蚀。因此，在处置库的运营期间，缓冲层峰值温度限制在100℃，但是考虑到10℃的温度安全裕度，允许的设计峰值温度设置为90℃^{[14, 16]}。因此，本章采用缓冲层峰值温度90℃的处置容器间距设计准则。图 11给出了不同围岩热传导系数条件下缓冲层峰值温度随处置容器间距的变化。考虑到围岩的非均质性和空间变异性，将围岩的导热系数设定为2.1，2.3，2.5，2.7，2.9 W/(K·m)。从图 11可以看出，当岩石导热系数为2.5 W/(K·m)时，合适的处置容器间距为7.7 m。岩石的比热容在计算中始终固定为常数832 J/(K·kg)。如果考虑导热系数和比热容的相关性，比热容的取值需要根据现场试验确定，半解析解的结果可以很容易进行修正。

图 11 缓冲层峰值温度随处置容器间距变化图

Figure 11. Variation of peak temperature of buffer layer with different container spacings

下载: 全尺寸图片幻灯片

图 11显示的结果是在处置隧道间距为40 m的条件计算得到的。在满足温度设计标准的前提下，应尽量减小处置隧道间距和处置容器间距，以实现处置容器布局优化，从而提高处置库存储容量。如果需要优化体积（或面积）而不是隧道总长度的布局，隧道间距应可能地短，通过解析模型也可以得到不同处置隧道间距条件下的类似结果。

4.2 预测原位试验结果

ATLAS (admissible thermal for argillanceous storage) Ⅲ原位加热试验是在比利时HADES地下实验室Boom黏土中进行的一项中等规模的原位加热试验^[24]。与之前进行的ATLAS Ⅰ和ATLAS Ⅱ原位加热测试相比，ATLAS Ⅲ加热测试增加了两个观测孔。ATLAS Ⅲ加热测试的目的是在更大的尺寸范围内获得更准确和更广泛的温度图像。

ATLAS Ⅲ加热测试开始于2007年4月，加热器功率逐渐增加，然后保持恒定值1400 W，并立即关闭，直到2009年11月测试结束。ATLAS Ⅲ试验装置由1个中央加热钻孔和4个观测钻孔组成。中央加热钻孔（AT89E）和3个观测钻孔（AT85E、AT93E、AT98E）都在同一水平面，而第4个观测孔（AT97E）略微向加热孔倾斜。AT97E钻孔向中央加热钻孔方向分别向下和向左倾斜10°。图 12显示了ATLAS Ⅲ加热测试的水平示意图^[24]。

图 12 ATLAS Ⅲ现场测试水平示意图^[24]

Figure 12. Horizontal view of ATLAS Ⅲ heating tests

下载: 全尺寸图片幻灯片

中央加热钻孔AT89E于1992年在试验隧洞89/90环处施工，井深19 m，直径230 mm。钻孔采用外径190 mm、内径160 mm的不锈钢管封闭。图 13为AT89E中央加热钻孔现场试验照片^[25]。

图 13 AT89E中央加热钻孔现场试验图^[25]

Figure 13. Field photo of central heating borehole AT89E

下载: 全尺寸图片幻灯片

AT85E和AT93E钻孔位于中央加热钻孔AT89E的两侧，两个观测钻孔分别位于中心钻孔左侧约1.5 m和右侧约1.3 m处。加热器位于中心钻孔AT89E中。加热器的长度为8 m，内径为160 mm，外径为190 mm。加热器的详细加热过程如表 2所示^[24]。

表 2 ATLAS Ⅲ试验加热和冷却信息

Table 2. Heating and cooling information of ATLAS Ⅲ heating tests

步骤	阶段	热功率/W	日期	天数序号	持续/d
1	加热	0→400	2007-04-02→2007-04-05	0→4	4
2	稳定	400	2007-04-06→2007-05-20	4→49	45
3	加热	400→900	2007-05-21→2007-05-25	49→54	4
4	稳定	900	2007-05-26→2007-07-29	54→120	66
5	加热	900→1400	2007-07-30→2007-08-03	120→125	5
6	稳定	1400	2007-08-04→2008-04-16	125→381	256
7	冷却	0	2008-04-17→2009-11-02	381→945	564

下载: 导出CSV

| 显示表格

在ATLAS Ⅲ原位加热试验中，热流量在钢管内壁施加，因此，在运用本章推导的半解析解时，把钢管和Boom黏土岩当成两种不同的传热介质。计算过程中所需的参数取值如表 3所示^[24]。图 14为AT85E和AT93E监测钻孔温度随时间变化的实测数据与不同解之间的对比图。从图 14中可以看出，在加热器功率加热稳定阶段，两个观测钻孔的温度都有所升高。在加热电源冷却阶段，观测钻孔温度迅速下降。同时，特别是在AT93E观测钻孔中，本文半解析解的结果可以准确预测现场测试的结果。半解析解的结果在AT85E观测钻孔的预测效果不如AT93E观测钻孔，但误差很小。这主要是由于在ATLAS Ⅲ原位加热试验中，Boom黏土岩是接近饱和的，地下水入渗造成黏土岩导热系数增大的影响可忽略不计，导致预测效果较好。

表 3 参数取值

Table 3. Values of parameters

参数	量值	参数	量值
λ_G/(W·(K·m)^-1)	16.3	ρ_Boom/(kg·m^-3)	2000
λ_Boom/(W·(K·m)^-1)	1.53	a/m	0.16
c_G/(J·(K·kg)^-1)	500	b/m	0.19
c_Boom/(J·(K·kg)^-1)	1434	t₀/℃	16
ρ_G/(kg·m^-3)	7850
注：表中下标G和Boom分别表示不锈钢钢管和Boom黏土。

下载: 导出CSV

| 显示表格

图 14 原位试验结果与不同解之间的对比

Figure 14. Comparison between in-situ test results and different solutions

下载: 全尺寸图片幻灯片

从图 14也可以看出，半解析解结果与Chen等^[24]的数值解进行了比较，半解析结果与数值模拟结果吻合较好，也验证了本文半解析解的正确性。与数值模型相比，解析模型也有自身优势。在进行参数敏感性研究时，特别是几何参数，例如，处置容器直径、缓冲层厚度、隧道或处置容器间距等，数值模型需要重新划分网格，涉及的计算工作量是繁重的。而解析法可以通过改变不同参数值而快速得到结果，不失是一种好的选择方案。同时，线热源解析解略低于本文半解析解，这主要是因为将圆柱形热源简化为线性热源，且忽略了不同介质间热扩散的差异。在ATLAS Ⅲ原位加热测试中，热功率在钢管内壁施加，且钢管厚度仅为15 mm。实测的温度数据都分布在Boom黏土岩，就无法评估线热源解在预测内层传热介质中温度变化的效果。在真实的处置库环境中，缓冲层的厚度大约在0.3 m左右，缓冲层与围岩之间热特性的差异不可忽视，如图 3（b）所示。因此，后续将借助原型处置库加热试验结果对不同模型解的预测能力进行评估。

5. 结论

建立了处置单元多屏障系统近场温度分布半解析解，研究了不同热参数和几何参数对缓冲层峰值温度的影响，并运用半解析解确定最小处置容器间距以及对原位试验结果进行预测。主要得到3点结论。

（1）当处置隧道间距小于40 m时，缓冲层峰值温度随处置隧道间距的减小而显著上升。然而，当处置隧道间距大于40 m时，其对峰值温度的影响很小。

（2）在一个处置板的中心位置，缓冲层峰值温度的影响半径为60 m。根据缓冲层温度设计准则，当隧道间距为40 m时，处置容器的最小间距为7.7 m。

（3）本文半解析解可以较好地预测原位加热试验结果。线性热源解略低于本文半解析解结果，同时半解析结果与数值模拟结果吻合较好，也验证了本文半解析解的正确性。

图 1 处置单元多重屏障系统计算简图

Figure 1. Model for multi-barrier system in disposal unit

下载: 全尺寸图片幻灯片

图 2 处置单元多重屏障系统数值分析模型

Figure 2. Numerical model for multi-barrier system in disposal unit

下载: 全尺寸图片幻灯片

图 3 本文半解析解与其他方法的对比

Figure 3. Comparison between semi-analytical method in this study and other methods

下载: 全尺寸图片幻灯片

图 4 不同影响半径下缓冲层温度变化

Figure 4. Temperature variation of buffer layer under different influence radii

下载: 全尺寸图片幻灯片

图 5 处置单元多重屏障系统近场温度分布图

Figure 5. Temperature distribution of multi-barrier system in disposal unit

下载: 全尺寸图片幻灯片

图 6 处置隧道间距对缓冲层峰值温度的影响

Figure 6. Effects of tunnel spacing on peak temperature of buffer layer

下载: 全尺寸图片幻灯片

图 7 处置容器间距对缓冲层峰值温度的影响

Figure 7. Effects of container spacing on peak temperature of buffer layer

下载: 全尺寸图片幻灯片

图 8 缓冲层厚度对缓冲层峰值温度的影响

Figure 8. Effects of thickness of buffer layer on its peak temperature

下载: 全尺寸图片幻灯片

图 9 缓冲层热传导系数对缓冲层峰值温度的影响

Figure 9. Effects of thermal conductivity of buffer layer on its peak temperature

下载: 全尺寸图片幻灯片

图 10 围岩热传导系数对缓冲层峰值温度的影响

Figure 10. Effects of thermal conductivity of rock on peak temperature of buffer layer

下载: 全尺寸图片幻灯片

图 11 缓冲层峰值温度随处置容器间距变化图

Figure 11. Variation of peak temperature of buffer layer with different container spacings

下载: 全尺寸图片幻灯片

图 12 ATLAS Ⅲ现场测试水平示意图^[24]

Figure 12. Horizontal view of ATLAS Ⅲ heating tests

下载: 全尺寸图片幻灯片

图 13 AT89E中央加热钻孔现场试验图^[25]

Figure 13. Field photo of central heating borehole AT89E

下载: 全尺寸图片幻灯片

图 14 原位试验结果与不同解之间的对比

Figure 14. Comparison between in-situ test results and different solutions

下载: 全尺寸图片幻灯片

表 1 参数取值

Table 1 Values of parameters

参数	量值	参数	量值
λ_B/(W·(K·m)^-1)	0.68	z₂/m	502.625
λ_R/(W·(K·m)^-1)	2.50	t₀/℃	10
c_B/(J·(K·kg)^-1)	1014	t₁/℃	34
c_R/(J·(K·kg)^-1)	832	P(0)/W	1545
ρ_B/(kg·m^-3)	1 770	ϕ	0.91
ρ_R/(kg·m^-3)	2 630	a₁	0.049
a/m	0.525	a₂	0.696
b/m	0.825	a₃	-0.059
h/m	5.25	a₄	0.271
L/m	1000	a₅	0.027
z₀/m	0	a₆	-0.010
z₁/m	487.375	a₇	0.026

下载: 导出CSV

表 2 ATLAS Ⅲ试验加热和冷却信息

Table 2 Heating and cooling information of ATLAS Ⅲ heating tests

步骤	阶段	热功率/W	日期	天数序号	持续/d
1	加热	0→400	2007-04-02→2007-04-05	0→4	4
2	稳定	400	2007-04-06→2007-05-20	4→49	45
3	加热	400→900	2007-05-21→2007-05-25	49→54	4
4	稳定	900	2007-05-26→2007-07-29	54→120	66
5	加热	900→1400	2007-07-30→2007-08-03	120→125	5
6	稳定	1400	2007-08-04→2008-04-16	125→381	256
7	冷却	0	2008-04-17→2009-11-02	381→945	564

下载: 导出CSV

表 3 参数取值

Table 3 Values of parameters

参数	量值	参数	量值
λ_G/(W·(K·m)^-1)	16.3	ρ_Boom/(kg·m^-3)	2000
λ_Boom/(W·(K·m)^-1)	1.53	a/m	0.16
c_G/(J·(K·kg)^-1)	500	b/m	0.19
c_Boom/(J·(K·kg)^-1)	1434	t₀/℃	16
ρ_G/(kg·m^-3)	7850
注：表中下标G和Boom分别表示不锈钢钢管和Boom黏土。

下载: 导出CSV

参考文献(25)

[1]	王驹, 陈伟明, 苏锐, 等. 高放废物地质处置及其若干关键科学问题[J]. 岩石力学与工程学报, 2006, 25(4): 801-812. doi: 10.3321/j.issn:1000-6915.2006.04.015 WANG Ju, CHEN Weiming, SU Rui, et al. Geological disposal of high-level radioactive waste and its key scientific issues[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(4): 801-812. (in Chinese) doi: 10.3321/j.issn:1000-6915.2006.04.015
[2]	叶为民, 王琼, 潘虹, 等. 高压实高庙子膨润土的热传导性能[J]. 岩土工程学报, 2010, 32(6): 821-826. http://cge.nhri.cn/cn/article/id/13418 YE Weimin, WANG Qiong, PAN Hong, et al. Thermal conductivity of compacted GMZ01 bentonite[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(6): 821-826. (in Chinese) http://cge.nhri.cn/cn/article/id/13418
[3]	陈航, 张虎元, 郭永强, 等. 混合型缓冲回填材料导热性能测试与预测研究[J]. 岩石力学与工程学报, 2014, 33(增刊2): 4312-4320. CHEN Hang, ZHANG Huyuan, GUO Yongqiang, et al. Measurement and prediction of thermal properties of bentonite-sand mixtures as buffer backfilling materials for high level radioactive waste[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(S2): 4312-4320. (in Chinese)
[4]	谈云志, 李辉, 王培荣, 等. 膨润土受热作用后的水-力性能研究[J]. 岩土力学, 2019, 40(2): 489-496. TAN Yunzhi, LI Hui, WANG Peirong, et al. Hydro-mechanical performances of bentonite respond to heat-treated history[J]. Rock and Soil Mechanics, 2019, 40(2): 489-496. (in Chinese)
[5]	ZHANG J R, SUN D A, YU H H, et al. Swelling of unsaturated GMZ07 bentonite at different temperatures[J]. Bulletin of Engineering Geology and the Environment, 2020, 79(2): 959-969. doi: 10.1007/s10064-019-01595-y
[6]	陈正汉, 秦冰. 缓冲/回填材料的热-水-力耦合特性及其应用[M]. 北京: 科学出版社, 2017. CHEN Zhenghan, QIN Bing. Thermal-Hydraulic-Mechanical Coupling Characteristics of Buffer/Backfill Materials and its Application[M]. Beijing: Science Press, 2017. (in Chinese)
[7]	陈正汉, 郭楠. 非饱和土与特殊土力学及工程应用研究的新进展[J]. 岩土力学, 2019, 40(1): 1-54. CHEN Zhenghan, GUO Nan. New developments of mechanics and application for unsaturated soils and special soils[J]. Rock and Soil Mechanics, 2019, 40(1): 1-54. (in Chinese)
[8]	陈正汉. 非饱和土与特殊土力学[M]. 北京: 中国建筑工业出版社, 2022. CHEN Zhenghan. Mechanics for Unsaturated and Special Soils[M]. Beijing: China Architecture & Building Press, 2022. (in Chinese)
[9]	CHO W J, KIM G Y. Reconsideration of thermal criteria for Korean spent fuel repository[J]. Annals of Nuclear Energy, 2016, 88: 73-82. doi: 10.1016/j.anucene.2015.09.012
[10]	刘月妙, 王驹, 蔡美峰, 等. 热-力耦合条件下高放废物处置室间距研究[J]. 铀矿地质, 2009, 25(6): 373-379. doi: 10.3969/j.issn.1000-0658.2009.06.009 LIU Yuemiao, WANG Ju, CAI Meifeng, et al. Study on disposal pit space for high-level radioactive waste in thermal-mechanical coupling conditiion[J]. Uranium Geology, 2009, 25(6): 373-379. (in Chinese) doi: 10.3969/j.issn.1000-0658.2009.06.009
[11]	CARSLAW H S, JAEGER J C. Conduction of Heat in Solids[M]. 2d ed. Oxford: Clarendon Press, 1959.
[12]	CLAESSON J, PROBERT T. Thermoelastic stress due to a rectangular heat source in a semi-infinite medium. Presentation of an analytical solution[J]. Engineering Geology, 1998, 49(3/4): 223-229.
[13]	HÖKMARK H, FÄLTH B. Thermal dimensioning of the deep repository[R]. Stockholm: Svensk Kärnbränslehantering AB, 2003.
[14]	HÖKMARK H, LÖNNQVIST M, KRISTENSSON O. Strategy for thermal dimensioning of the final repository for spent nuclear fuel[R]. Stockholm: Svensk Kärnbränslehantering AB, 2009.
[15]	刘东东, 项彦勇. 高放射核废处置库温度场的分布线热源解析模型[J]. 岩石力学与工程学报, 2019, 38(增刊1): 2816-2822. LIU Dongdong, XIANG Yanyong. A distributed line heat-source analytical model for the temperature field of a high level nuclear waste repository[J]. Chinese Journal of Rock Mechanics and Engineering, 2019, 38(S1): 2816-2822. (in Chinese)
[16]	IKONEN K. Thermal analysis of repository for spent EPR-type fuel[R]. Olkiluoto: Posiva Oy, 2005.
[17]	HARTLEY L, HOCH A, JACKSON P, et al. Groundwater flow and transport modeling during the temperate period for the SR-Can assessment. Forsmark area-version 1.2[R]. Stockholm: Svensk Kärnbränslehantering AB, 2006.
[18]	陈永贵, 贾灵艳, 叶为民, 等. 施工接缝对缓冲材料水-力特性影响研究进展[J]. 岩土工程学报, 2017, 39(1): 138-147. doi: 10.11779/CJGE201701012 CHEN Yonggui, JIA Lingyan, YE Weimin, et al. Advances in hydro-mechanical behaviors of buffer materials under effect of technological gaps[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(1): 138-147. (in Chinese) doi: 10.11779/CJGE201701012
[19]	XU X, HE L Q, SUN D A, et al. Validation of the fully-analytical solution with temperature superposition for the nuclear waste repository[J]. Nuclear Engineering and Design, 2023, 409: 112367. doi: 10.1016/j.nucengdes.2023.112367
[20]	陶文铨. 传热学[M]. 5版. 北京: 高等教育出版社, 2019. TAO Wenquan. Heat Transfer[M]. 5th ed. Beijing: Higher Education Press, 2019. (in Chinese)
[21]	HE L Q, ZHOU X Y, SUN D A. Fully analytical solution in time and space domains on temperature in multi-barrier nuclear waste repository[J]. Computers and Geotechnics, 2023, 154: 105164. doi: 10.1016/j.compgeo.2022.105164
[22]	CRUMP K S. Numerical inversion of Laplace transforms using a Fourier series approximation[J]. Journal of the ACM, 1976, 23(1): 89-96. doi: 10.1145/321921.321931
[23]	徐国庆. 国际高放废物处置研发工作在花岗岩地区的进展[J]. 世界核地质科学, 2016, 33(2): 119-124. doi: 10.3969/j.issn.1672-0636.2016.02.010 XU Guoqing. Abroad progress in R & D work on high-level radioactive waste disposal in granite areas[J]. World Nuclear Geoscience, 2016, 33(2): 119-124. (in Chinese) doi: 10.3969/j.issn.1672-0636.2016.02.010
[24]	CHEN G J, SILLEN X, VERSTRICHT J, et al. ATLAS Ⅲ in situ heating test in boom clay: field data, observation and interpretation[J]. Computers and Geotechnics, 2011, 38(5): 683-696. doi: 10.1016/j.compgeo.2011.04.001
[25]	DE BRUYN D, LABAT S. The second phase of ATLAS: the continuation of a running THM test in the HADES underground research facility at Mol[J]. Engineering Geology, 2002, 64(2/3): 309-316.