许 迅 明华军，2 谈云志

(1. 三峡大学 特殊土资源化利用宜昌重点实验室， 湖北 宜昌 443002； 2. 三峡大学 水利与环境学院， 湖北 宜昌 443002)

核能的开发利用产生了大量的高放射性废物(简称“高放废物”)．如何安全地处置高放废物，已成为国家核能利用必须解决的问题[1]．深地质处置方法被国际上公认为永久处置高放废物的最佳措施[2]，即将高放废物固体化后，放置于距地表深约500～1 000 m的完整地质体坑道中，在其间设置缓冲回填材料，形成人工屏障，阻隔放射性核素迁移，实现放射性核素与生物圈的永久隔离．缓冲层是填充在废物罐和地质体之间的重要人工屏障，是处置库安全性和稳定性的有效保障[3]，因此，对处置库导热功能的改进，主要着重于对缓冲层的设计．

膨润土具有低渗透、高膨胀和强吸附的特性，是高放废物处置库首选的缓冲材料[4]．膨润土导热性能相对较差，即使高压实(ρd=1.80 g/cm3)状态，也只有1.3 W/(m·K)，远低于常规地质体热传导系数(2.43 W/(m·K))[5-6]，导致高放废物所产生的衰变热消散缓慢，引起热量聚积．据Lee等[7]计算，处置库封场后前10年，废物罐表面温度最高至107℃，超出处置库安全运行临界值(t=100℃)[8]，导致水分蒸发，进而引起膨润土收缩，促使局部裂隙衍生，损伤缓冲层屏障功能[9]，因此，如何使衰变热快速向周围扩散是高放废物处置库安全运营的关键问题之一．

许多学者采用在膨润土中掺入具有高导热性能的材料(石英砂、岩屑、石墨等)的方法来提高膨润土的导热性能[4，10-12]．叶为民， 王琼，潘虹，等[4]利用热探针法研究了干密度、含水率对加砂高庙子膨润土导热系数的影响，发现在其他参数相同的条件下，导热系数随干密度、含水率及掺砂量的增加而增大．马超[10]和谢敬礼，马利科，高玉峰，等[11]均通过添加花岗岩岩屑改善高庙子膨润土的导热性能，结果表明添加岩屑提升导热性能效果要劣于添加石英砂，但岩屑作为添加剂对处置库水化学环境影响小，相容性较好．刘月妙，蔡美峰，王驹[12]研究了在GMZ01膨润土中分别掺入石英砂、石墨后的导热性能，发现石英砂和石墨均可以明显提高GMZ01膨润土的导热性能．

石英砂等材料在提升膨润土导热性能的同时，也显著地降低了其防渗能力．Cho等[13]进行了不同掺砂率的Kuungju膨润土-石英砂混合物渗透试验，发现掺砂率70%时混合物的渗透系数在10-10m/s数量级左右，不能满足国际原子能机构(IAEA)对缓冲材料渗透性能的要求(<10-12m/s)[14]．因此，拟将缓冲层分为两层，即外部缓冲层采用纯膨润土材料，以保持防渗性能，而将掺入剂(如：石英砂、石墨)掺入内部缓冲层的膨润土中，以增加缓冲材料的导热性能．本研究主要通过数值模拟的手段对处置库工程屏障进行热力学分析，进而评估缓冲层分区设计的可行性．

1 双层缓冲层概念

为了有效提高缓冲材料导热性能，同时保证其防渗性能也达到要求，提出了双层缓冲层的概念．图1为单层缓冲层和双层缓冲层的概念模型．假设内部缓冲区的大小是原始缓冲区的一半．在双层缓冲层中，外部缓冲层为纯膨润土，内部缓冲层为膨润土-掺入剂的混合物．对于双层缓冲层，由于外部缓冲层依然为纯膨润土，故保持了原有的防渗功能，可以有效减少地下水的渗入．

图1 高放废物处置库的单层与双层缓冲层概念模型(单位：mm)

2 数值模型

2.1 热传导方程

热力学第一定律指出，微元体Ω内动能KΩ和内能EΩ的变化率与机械功率Pext和微元体之间的热量交换率Qexch相等，用数学语言可表示为:

(1)

在微元体Ω内，假设单位质量的内能为E，则微元体的总内能EΩ为：

(2)

其中：ρ为材料密度(kg/m3)；cp为材料比定压热容，表示单位质量材料升高单位温度所需要的热量[J/(kg·K)]；T为材料温度(K)．

由于计算不考虑外力影响，故Pext=0．

热量的交换功率Qexch包括通过传导、辐射、对流和热源作用下的热量交换，由于固相介质中热量的传递主要通过热传导进行，所以在表面积为S的微元体Ω中交换功率可表示为：

(3)

其中：λ为材料热传导系数[W/(m·K)]．

结合式(1)～(3)，可得材料的热传导方程为：

(4)

由于土体材料具有三相性，且各相的热力学特性均不同，由比定压热容定义，可得如下公式：

ρcp=ρdcp，s+ρwnSrcp，w+ρan(1-Sr)cp，a

(5)

其中：ρd为土的干密度；ρw为水的密度；ρa为空气的密度(kg/m3)；n为孔隙比；Sr为饱和度；cp，s、cp，w、cp，a分别为土体颗粒、水、空气的比定压热容[J/(kg·K)]．具体物理参数见表1．

表1 热力学参数

由于土体热传导系数与土颗粒、土中水、土中气赋存状态有关，故本文中纯膨润土[7]和掺砂膨润土[15]的热传导系数可分别采用如下公式：

λ=4.588w+0.124

(6)

(7)

其中：w为初始含水率；λw和λa分别表示水和空气的热传导系数；Φs、Φb、Φw、Φa分别为砂、膨润土、水和空气的体积分数．

Q为体积热源(W/m3)，其衰减函数可表示为[7]：

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

(8)

其中：t为处置库运行时间(年)．

2.2 边界及初始条件

图2为缓冲区的俯视图，根据该图对模型水平方向的边界条件进行定义．

图2 高放废物处置库的双层缓冲层概念模型俯视图

1)边界条件

由上文已知，废物罐释热为高放废物的衰变热，故模型内边界，即废物罐表面(r=a)处热通量随时间衰减，相应边界条件方程见式(8)．

为避免外边界对计算结果的影响，需将外边界设置在足够远处，定义20倍热源直径为外边界，且为恒温边界，不受温度场影响，表示为：

T(∞，0)=T0

(9)

其中：T0为初始温度，由式(10)确定．

2)初始条件

模型中各处初始温度均为所处地深的地温，该处地表温度假设为15℃，地层温度梯度为3℃/100 m，即：

T0=15+0.03z

(10)

其中：z为埋深(m)．

3)连续性条件

假设内部缓冲材料与外部缓冲材料在r=b处完全接触，且外部缓冲材料与围岩在r=c处完全接触，则必须满足连续性条件：

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的热传导系数．例：T1(b，t)为时间为t时在r=b处Buffer1的温度．

2.3 几何模型

由于本研究主要关注工程屏障的热力学行为，而不是整个处置系统的温度分布，因此，以单个处置库为研究对象．根据瑞典KBS-3V概念模型建立单个处置库系统模型，系统包括封装高放废物固化体的废物罐、缓冲材料、回填材料和围岩，将缓冲材料分为Buffer1、Buffer2两部分．为了保证计算的可执行性，考虑单个处置库的空间对称性，使用二维轴对称模型代替三维模型进行分析．

为了避免边界条件对计算结果的影响，取计算模型尺寸为20 m×1 000 m，将坐标原点(0，0)设置在废物罐中心，如图3所示．图4是计算模型的网格划分，为了精确分析缓冲、回填区的温度分布情况，故靠近废物罐的部分网格划分更加密集．完整网格包含6 985个单元和7 264个节点．

图3 双层缓冲层处置库几何模型 图4 图3对应几何模型的网格划分

3 模拟与讨论

3.1 模型验证

Hodgkinson[16]针对随时间指数衰减的线性热源，建立实体圆柱型模型，给出了对应的半解析解，模型中各点温度随时间变化的表达式如下：

(15)

其中：R，z分别为径向坐标和轴向坐标；t为时间；b为圆柱体的半长；r为其半径；Q为热源衰变热，见式(8)；q0为初始热源衰变热，由式(8)可知q0=6.353 kW/m3；erf(x)为误差函数，表示为：

(16)

I0(x)为修正的贝塞尔函数[17]，表示为：

(17)

为了验证本文模型方法的适用性，使用COMSOL Multiphysics 5.3a构建如图5左侧所示的二维轴对称模型，根据表2参数进行了热力学数值模拟分析，将结果与解析解进行了比较．

表2 验证模型相关参数

两种方法计算的缓冲层在金属罐中心的温度变化对比如图5所示．

图5 Hodgkinson解析解问题[16]对应数值模型与结果验证

解析解和数值分析计算得到的峰值温度分别为71.77℃和72.64℃．峰值温度的相对误差为0.71%，且两种方法得到达到峰值温度均需要两年时间．整个过程中的平均误差为7.13%．由此可知，在该模型中，使用COMSOL Multiphysics 5.3a进行热力学分析得到的数值解与解析解较吻合，因此，可以使用COMSOL Multiphysics 5.3a采用上述物理方程对高放处置库的热力学行为进行分析．

3.2 参数敏感性分析

(16)

参数敏感性函数可以比较系统特性对各参数的敏感程度，敏感性函数值越大，系统对该参数越敏感．每个基准数据对应6个分析数据，各参数的分析数据一般在参考数据的0.5～1.5倍范围内，个别参数在此范围上有所缩放，以求涵盖该参数的合理范围．

热力学计算中，采用常规单层缓冲材料进行参数敏感性分析，两个主要参数分别为热传导系数λ，比定压热容cp，故设定基准状态X*=(λ，cp)，系统特征F*=f(λ，cp)．由谈云志，彭帆，钱芳红，等[19]及陈航，张虎元,郭永强，等[20]的试验及预测，缓冲材料的热传导系数在0.4～1.5 W/(m·K)的范围内，比定压热容在500～1 700 J/(kg·K)的范围内，故取基准数据λ=1.4 W/(m·K)、cp=1 000 J/(kg·K)，计算得到的处置库峰值温度作为基准值，分析处置库峰值温度对上述两个参数的敏感程度，结果见表3．

表3 敏感性参数计算表

由表3可计算热传导系数和比定压热容的平均敏感性函数值分别为24.83和0.11．热传导系数的敏感性函数值远大于比定压热容敏感性函数值．

由图6的曲线斜率大小可以判断，材料比热容变化对峰值温度几乎没有影响；热传导系数变化对峰值温度的影响较大．

图6 峰值温度-参数变化率关系

但随着导热系数的增加，效率增益逐渐减小，这是因为缓冲材料的导热系数超过处置库周围的围岩后，缓冲层和围岩整体形成了热的不良导体，废物罐表面热量虽然能够快速通过缓冲层到达与围岩的交界面，但是由于围岩的导热系数较小，热量不能及时导出，导致温度变化较小．因此，在对比单层和双层缓冲区对峰值温度的影响时可忽略比热容的影响，以简化计算．

3.3 温度分布及对比

根据第2节所述单层、双层缓冲层数值模型，利用有限单元计算程序COMSOL Multiphysics 5.3a计算了100 a内模型温度演变．

图7为填充单层缓冲材料的处置库在运行2 a后，处置库坑道周围温度的分布情况．峰值温度出现在废物罐表面的中心位置，并通过废物罐、缓冲区、回填区向围岩扩散．离废物罐越近，温度越高．且回填区的温度明显低于缓冲区．

图7 单层缓冲材料处置库运行2年后温度分布

图8为填充单层缓冲材料的处置库自开始运行后，废物罐表面中心位置(如图中点A)温度随时间的演化情况．废物罐表面A点温度在处置库开始运行1.8 a后即达到峰值温度110.65℃，然后由于衰变热的降低，随运行时间而逐渐降低，处置库运行100 a后温度降至57℃．

图8 单层缓冲材料处置库废物罐表面A点温度时间演化

图9(a)、(b)分别为单层缓冲材料与双层缓冲材料对应处置库缓冲区3个特征点(点A、B、C)的温度随时间的演化过程．由图可知，两种结构处置库均在运行1.5 a左右在废物罐表面中心位置(点A)达到峰值温度．双层缓冲材料处置库峰值温度为92.73℃，单层缓冲材料处置库峰值温度为110.65℃，可见，将处置库缓冲层进行分区设计可以明显降低处置库峰值温度，温度总体下降了16%；同时，随着时间变化两种结构A点温度趋于接近，处置库运行100 a后均趋于60℃．由此可见，双层缓冲区结构可以显著降低处置库温度最不利时段的运行环境温度，对处置库安全稳定运行提供了有效保障．

图9 缓冲区监测点温度变化

4 结 论

本研究采用热力学分析方法，使用COMSOL Multiphysics 5.3a有限元软件，对瑞典KBS-3V处置库概念模型进行相应分析，确定了热力学敏感性参数，并通过对比单、双层缓冲层处置库设计方案对处置库中峰值温度的影响，得到了如下结果：

1)材料比热容变化对峰值温度几乎没有影响；热传导系数变化对峰值温度的影响较大，但随着导热系数的增加，效率增益逐渐减小．

2)相对单层缓冲层结构，采用双层缓冲层设计可以将处置库峰值温度降低约16%，可以使处置库在临界温度(100℃)以下安全运行．