0.   引言

高放废物的安全处置是关系到核工业可持续发展和保护环境、保护公众健康的重大问题。目前国际上公认的高放废物处置方式为深地质处置，采用“多重屏障系统”进行处置库设计^[1]。缓冲材料作为废物罐和处置围岩之间重要的人工屏障，在放射性衰变热、地下水入侵和围岩应力等作用下，将经历复杂的热-水-力耦合过程。因此，评价其长期性能对高放废物处置库的安全至关重要^[2]。膨润土由于具有极低的渗透性、良好的膨胀性、持水性及长期稳定性，现已被国际上多个国家选作高放废物深地质处置缓冲材料，如MX-80、FEBEX、Kunigel V1、高庙子膨润土等。

为了揭示缓冲材料在高放废物处置库条件下的热-水-力耦合行为演化特征，开展工程尺度的模型试验是一种重要途径。然而，由于试验时间尺度的限制，开展相应的数值模拟研究是预测和评价缓冲材料长期耦合性能的必要手段。据此，国际上多个高放废物地质处置地下实验室采用工程尺度模型试验和数值模拟相结合的方式研究缓冲材料热-水-力多场耦合性能。针对MX-80膨润土，瑞典SKB先后分别采用不同数值软件/程序ABAQUS、CODE_BRIGTH、COMPASS，基于不同尺度的室内及原位试验，提出了模拟处置库条件下温度-水压-膨胀力行为演化数值模拟方法^[3-5]。并且在瑞士Mont Terri地下实验室开展了全尺寸原位试验（FE）和1∶2尺寸的加热试验（HE-E）研究，分别采用程序TOUGH 2、OpenGeoSys等模拟并验证了MX-80膨润土多场耦合性能^[6-7]。对于FEBEX膨润土，类似研究在西班牙Madrid和瑞士Grimsel分别开展了室内mock-up试验和现场原位试验，同样使用不同程序/软件确立了FEBEX膨润土多场耦合理论模型及相应的模拟方法^[8-10]。在加拿大AECL地下实验室，选择了膨润土与砂混合材料作为缓冲材料，通过工程尺度的现场试验（BCE和ITT）研究并模拟围岩与压实缓冲材料之间的长期演变关系^[11-13]。与此同时，在其它国家也开展了类似试验与数值模拟工作，例如：韩国的KENTEX试验、捷克Mock-Up-CZ试验、比利时的OPHELIE Mock-up试验等^[14-16]。国际上专门成立了DECOVALEX计划开展高放废物地质处置系统围岩中多个物理场的耦合过程研究和工程屏障的可行性评估，并提出了描述缓冲材料各种不同力学行为的本构模型，例如：BBM模型、BExM模型等^[17-19]。

目前中国已经筛选出内蒙古高庙子膨润土作为处置库缓冲材料首选矿床，并开展了大量关于渗透性、导热性、力学特性的试验研究^[20-22]。为了模拟缓冲材料在处置库条件下的行为演化特征，核工业北京地质研究院于2011年设计并建造了工程尺度的大型缓冲材料多场耦合性能模型试验系统（China-Mock-up）^[23]。Chen等^[24]、Zhao等^[25]提出了高庙子膨润土热-水-力多场耦合理论模型及数值模拟方法，初步验证了理论模型与数值模拟方法的可靠性。本文则是在前期研究成果的基础上，采用有限元数值模拟软件LAGAMINE，考虑试验过程中真实复杂边界加载条件和5种不同的材料属性，实现了近5 a试验数据的定量模拟，准确拟合了缓冲材料内部温度、相对湿度及总应力演化规律，建立并验证了缓冲材料多场耦合性能理论模型与数值模拟技术，为后续处置库条件下缓冲材料长期性能评价提供技术支撑。

1.   模型试验简介

根据中国高放废物处置概念模型，China-Mock-up试验被设计成圆柱体形状，主要组件包括外围的碳钢腔体（内径900 mm×高2200 mm，厚度25 mm）、膨润土破碎颗粒、压实膨润土块体、加热器（直径300 mm×高1600 mm）、数据采集系统、加水系统、温控系统、传感器等，如图 1所示，其中加热器又可分为加热芯和氧化镁导热材料。在加热器与膨润土块之间、膨润土块与碳钢台架之间分别添加了15，50 mm厚的破碎颗粒。压实膨润土多为扇形块体，干密度为1.71 g/cm³，初始含水率为8.7%，初始吸力为80 MPa。

图  1  China-Mock-up模型试验装置结构剖面图^[23]

Figure  1.  Sketch of China-Mock-up facility ^[23]

下载: 全尺寸图片 幻灯片

为了系统获得缓冲材料内部温度、湿度、应力、变形等数据，试验装置内部分别布设了温度传感器、湿度传感器、应力传感器、位移传感器、孔隙水压力传感器、电阻应变片等10种不同类型的传感器，布设过程中确保每一层及各层之间的传感器具有关联互补功能。该试验系统是中国目前尺寸最大的缓冲材料热-水-力-化学耦合试验台架，可用于模拟高放废物处置库环境下的缓冲材料长期耦合性能试验研究。

2.   热-水-力多场耦合数学模型

为了更好地模拟China-Mock-up试验过程中膨润土内部的热量运移、液态水传播、水蒸气扩散与冷凝、土体膨胀/收缩等复杂行为，Chen等^[24]提出了描述高庙子膨润土热-水-力多场耦合行为的理论模型。本文采用相同的理论模型，主要包括水分运移方程、热量运移方程及力学本构模型。

2.1   水分运移方程

膨润土内部水的存在形式分为两种：液态水和水蒸气。基于质量守恒定律，土体中水分运移规律的控制方程可以通过两者质量守恒方程叠加获得，即

式中：${\rho _\mathrm{w}}$和${\rho _\mathrm{v}}$为水与水蒸气的密度；n为孔隙度；t为时间；S_{r, w}和S_{r, g}为水、水蒸气的饱和度；f_w，f_g为水与水蒸气的速度；i_v为水蒸气弥散作用。为了确保数值模型具有较好的收敛性与稳定性，在模拟过程中假设试验台架内部气压保持不变，忽略气体扩散这一过程，故在多孔介质中只考虑液态水的运移过程，其满足广义达西定律：

式中：p_w为压力；k_int，k_{r, w}分别为饱和渗透率和相对渗透率；μ为动力黏滞系数。

2.2   热量运移方程

热量运移模型是基于能量守恒定律提出的，假设材料各组成部分处于热动平衡状态，所有成分可以由温度指标定义。在温度场控制方程中涵盖了固体、水、气体等各个介质的能量守恒，可表示为

式中：E_s，E_w，E_g分别为固体、液体及气体的内能；i_c为热传导作用；j_Es，j_Ew，j_Eg分别为固体、液体及气体中的热对流作用；Q为源汇项。

2.3   力学本构模型

本文力学本构模型采用是Alonso等^[21]提出的BBM模型，能够较好地反映饱和度对土体的力学特性显著影响。其LC屈服面在（p, q, s）三维空间中可以描述为

式中：p为平均净应力；q为偏应力；M为临界状态曲线倾斜度；p₀为前期固结压力；p_s为考虑吸力的存在对土体黏聚力的增强作用。在（p, s）二维空间存在LC屈服线和SI屈服线，分别可以表示为

式中：p₀^*为饱和时的固结压力；p_c为参考应力；s₀为前期固结时最大吸力；k为弹性压缩曲线斜率；$\lambda (0)$和$\lambda (s)$为饱和状态下吸力为s时的压缩参数，其函数关系为

式中：r为土体受吸力影响下的最大刚度参数；β为控制土体刚度系数随吸力增长速率的参数。其中模型硬化准则可表示为

式中：${\lambda _\text{s}}$和k_s分别为吸力引起的体积变形塑性刚度参数和弹性刚度参数；p_at为标准大气压力；$\varepsilon _{\text{v}}^{\text{p}}$为塑性体积应变。

3.   热-水-力多场耦合数值模型

3.1   模型结构与网格剖分

考虑到China-Mock-up模型试验装置的实际结构，采用2D轴对称模型，分别考虑加热芯、氧化镁导热材料、外围碳钢腔体、压实膨润土块和填充的破碎颗粒五种不同属性材料，结果如图 2所示。每一种材料属性都假设是均质各向同性的，忽略膨润土块体、碎块颗粒间的缝隙，其中膨润土碎块颗粒认为其饱和渗透率比压实的块体偏大，分别取值为2.5×10^-18，2.5×10^-20 m²，其它参数取值与膨润土块相同。本研究采用比利时列日大学开发的有限元数值软件LAGAMINE，其可实现不同类型网格剖分和局部加密，提供了一系列描述土体力学特性的本构方程（例如：莫尔库仑模型、BBM模型等），能够定量模拟多场耦合条件下土体内部热量运移、水分运移、气体扩散、土体变形等问题。模型采用的是矩形网格，网络单元总数为1482，结点总数为4244（图 2）。选取约5年的温度、相对湿度、总应力试验数据进行模型验证工作。

图  2  China-Mock-up试验数值模型与网格剖分结果

Figure  2.  Numerical model for China-Mock-up experiment and mesh grids

下载: 全尺寸图片 幻灯片

3.2   模型参数的确定

前期研究已开展了土水特征试验、非饱和渗透试验、压缩试验、导热性能试验等，为本研究提供了模型所需的多场耦合参数。首先对于非饱和土，吸力s与饱和度S_{r, w}的函数关系可表示为

式中：S_{r, u}为土体的最大饱和度，取1.0；S_{r, res}为在土体的残余饱和度，其值为0.1。根据土水特征试验数据，采用最小二乘法进行数值拟合，结果如图 3。式（8）中参数为a₁=7.0×10^-6 Pa^-1，a₂=0.9，a₃=70。

图  3  高庙子膨润土土水特征曲线

Figure  3.  Water-retention curve of Gaomiaozi bentonite

下载: 全尺寸图片 幻灯片

非饱和渗透试验采用瞬时截面法来测定膨润土的渗透系数，通过土样特定截面的流量以及该截面的水力梯度，得到土体在特定时间特定截面上的平均渗透系数。根据前期试验数据，相对渗透率k_{r, w}拟合公式为

拟合结果见图 4。

图  4  饱和度与相对渗透率之间的拟合曲线图

Figure  4.  Relative permeability of GMZ bentonite as a function of degree of saturation

下载: 全尺寸图片 幻灯片

基于热传导试验，高庙子膨润土饱和度S_{r, w}与导热系数k的关系为（图 5）

图  5  高庙子膨润土导热系数与饱和度的关系

Figure  5.  Thermal conductivity of GMZ bentonite as a function of degree of saturation

下载: 全尺寸图片 幻灯片

由于现阶段获得的模型参数相对比较有限，其余则是参考国外具有代表性膨润土的试验数据^[24]。

3.3   模型边界条件设置

在试验过程中实际温度边界分为加热芯和试验腔体外围两部分。加热芯温度是根据试验设计的温度逐步升高的，最终稳定保持在90℃。在试验开始前255 d，处于试验调式阶段，加热过程中加热芯的温度控制相对较差，使得试验初期实测数据的变化幅度相对较大。模型中对加热过程进行简化处理：0~135 d，温度保持在30℃；136~255 d，温度从30℃升高至90℃；256 d之后，温度维持在90℃。按照此方式采用阶梯函数方式模拟这一加热过程，并对加热芯处的结点进行温度赋值，如图 6所示。在试验腔体外围虽然添加了保温层，但是试验前1500 d保温效果相对较差受室温影响较为明显，呈季节性波动，之后由于试验环境的改变，温度保持相对较小的波动范围。基于此，数值模型刻画外围温度边界时主要是采用傅里叶函数与线性函数分段拟合试验数据进行边界设置，其中在700 d左右室温出现小幅度明显的波动，这一过程采用人为方式进行边界处理，如图 7所示。

图  6  试验加热过程及数值拟合结果

Figure  6.  Heating process and numerical results of heater strip

下载: 全尺寸图片 幻灯片

图  7  室温与数值拟合结果随时间的变化

Figure  7.  Variation in room temperature and numerical results with time

下载: 全尺寸图片 幻灯片

试验注水管固定于腔体内壁上，水是从试验台架底部注入，注水量随时间逐渐累积增大的，如图 8（a）所示。在实际的注水过程中，受重力作用影响，试验台架底部的水量相对于上部区域较多。为了更真实地模拟试验真实注水条件，在模型中注水边界分为两部分，见图 2。试验前期上部注水量相对较小，下部的注水量更大，随着下部膨润土逐渐趋于饱和，下部注水量随时间逐渐减少，上部则相应逐步增大。考虑试验过程中影响每天注水量的试验因素很多，实际每天注水量并不完全相同，数值模型采用平均注水速率的简化方式分别对上部和下部水量边界赋值描述整个的注水过程，如图 8（b）所示。

图  8  注水量及平均注水速率随时间的变化

Figure  8.  Water injection process in China-Mock-up test

下载: 全尺寸图片 幻灯片

4.   结果与讨论

4.1   温度场模拟结果

试验台架不同位置内部特征点温度随时间的变化如图 9所示。结果表明模拟值与试验数据的吻合程度较高，能够很好地反映试验台架内部温度随时间演变这一过程，这充分说明在本研究中各个材料设定的热学参数及边界条件是合理的。其次，图 9内部特征点的温度曲线可以说明，在试验台架内部的温度受室温波动非常明显，而在数值模型中，采用分段函数方式刻画外部温度边界条件，能够很好地反演出台架内部温度前期随室温波动及后期线性变化这一规律。同时在0~255 d期间，试验台架内部处于加热过程，各个特征点的温度是加热芯加热温度与室温综合作用的结果。例如：从135~255 d模型加热芯边界设置的温度是阶梯的方式逐渐增高，而外围边界室温的温度是逐渐降低的，模拟结果是这两方面综合作用的结果。对于图 9（a），（b）中的特征点距离加热芯相对较远，模拟结果主要受室温的影响呈现降低趋势；对于图 9（c）~（f）中的特征点距离加热芯相对较近，在这一时间段加热芯的温度从80℃升高到90℃，相应模拟温度也呈现升高趋势。之后加热芯影响的区域随时间逐步稳定，然后室温对模拟结果呈现主导作用，模拟温度呈现下降趋势，同时在174 d左右类似的现象也出现过。但是此期间内部特征点温度整体趋势是增加的，说明在这一过程中大趋势仍然受加热芯加热的影响。同时从700~720 d这一阶段，监测到的温度数据由于室温作用影响出现较小幅度的波动，该现象可以通过模型中边界条件的精细刻画得到准确地预测。在1500 d之后，室温波动幅度较小，本数值模型通过采用线性的方式也可较好地模拟这一规律。

图  9  试验台架膨润土内部不同区域特征点温度的模拟结果

Figure  9.  Comparison between predictive temperature and experimental results at different locations of China-Mock-up facility

下载: 全尺寸图片 幻灯片

图 9（e），（f）模拟结果还表明在试验台架中上部特征点反映出数值结果低于试验数据，其原因是加热器开始加热后，水分会随着温度升高发生蒸发现象，产生水蒸气，导致膨润土内部局部温度较高区域的气体压力升高，气体会携带着部分热量向上部和外部运移。由于试验台架装置采用的是封闭式设计，在台架外围还封装了一层保温材料，故试验台架上部区域的热量随着时间推移会逐渐聚集，导致上部温度逐渐升高。在数值模拟过程中，为了获取较好的数值收敛结果，整个模拟过程中假设膨润土内部不同位置的气体压力随时间保持不变，忽略了气体迁移过程中携带的热量，使得数值结果低于试验数据。同时，图 9（a），（b）模拟结果则说明在试验台架底部受这一现象影响不明显。

图 10为不同时刻试验台架内部温度场的分布云图，该结果可较好地反映不同材料导热性能差异。根据图 10可知，在加热芯处的模拟温度最高，其次氧化镁材料具有良好的导热性能，其温度与加热芯的温度差异不是很明显。相对而言，膨润土块体及破碎颗粒的导热系数相对较低，且靠近外围边界，受室温的影响，其内部温度场不同位置差异性较为明显。其次，模拟结果还表明在100 d时加热芯的温度保持在30℃，温度与室温接近，试验台架内部温度场整体差异不大。后续随着加热芯的温度升高，氧化镁导热材料内部温度快速升高，当达到300 d左右时，温度场已经趋近于较为稳定的状态，主要受室温变化发生相应的波动，例如：在300 d时，室温接近8.7℃，之后室温逐渐升高，到400 d时室温则达到21.5℃，通过温度场云图可较好地反映试验台架内部温度场随时间的变化特征。通过对比1600 d与1800 d温度场云图可知，在这一期间，室温呈近似线性波动，变化幅度较小，相应的温度场分布差异性不是很明显。

图  10  试验台架内部不同时刻温度场分布云图

Figure  10.  Temperature distribution at different time

下载: 全尺寸图片 幻灯片

4.2   渗流场模拟结果

图 11为试验台架内部不同区域特征点相对湿度随时间的模拟结果。通过数值结果与试验数据的对比可知，本研究中采用的模拟方法整体上能够反映膨润土内部水分在多场耦合条件下的运移规律。在试验开始运行阶段，尚未开始加热和注水，在这一阶段相对湿度随时间保持不变。随着加热器开始运行，靠近加热器的特征点，相对湿度发生明显降低，之后随着注水作用的影响，相对湿度相应地开始逐渐升高并趋于饱和，如图 11（b）~（d），（f）所示，在这一阶段，相对湿度曲线表现出先干燥后饱和的趋势。然而在图 11（a），（e）中上述现象并未出现，其原因主要是这两个特征点位置靠近注水边界，距离加热芯相对较远，受注水过程影响比较明显，相对湿度随时间逐渐升高趋势。图 11的模拟结果还可以反映出试验数据与数值结果在前期吻合程度较好，例如：对于图 11（b）~（d）中的前600 d而言，模拟结果能够较好反映试验过程中相对湿度随时间变化规律。但是随着时间推进，试验监测到的相对湿度增长的速率明显高于模拟值，其原因是膨润土砌块之间、砌块与传感器之间以及膨润土颗粒存在缝隙，在实际注水过程中，水分会优先沿着缝隙渗入到膨润土内部，扩散至传感器的位置。在数值模拟的过程中是假设膨润土块之间和颗粒之间不存在缝隙，是均质各向同性的整体，水分是按照不同材料的渗透性逐渐从外向内缓慢迁移，迁移速率远小于沿缝隙运移的速度，使得模拟结果小于试验监测数据。

图  11  试验台架膨润土内部不同区域特征点相对湿度的模拟结果

Figure  11.  Comparison between predictive relative humidity and experimental results at different locations of China-Mock-up facility

下载: 全尺寸图片 幻灯片

膨润土内部不同时刻水分迁移过程中渗流场模拟结果如图 12所示。在前100 d，由于暂时未注水，膨润土内部的水压与初始吸力保持相同，随着时间推移，加热器周围区域的水压受加热影响开始降低，在600 d最小水压可达-230 MPa左右，这一模拟结果与图 11（b）~（d）相对湿度随时间的变化规律是保持一致的。后续随着注水量累积增多，水分由外向内缓慢运移，加热产生的干燥区域范围逐步缩小，土体内部的水压逐渐升高。同时渗流场云图还表明经过近5年水分运移，试验台架内膨润土仍尚未完全达到饱和，最大水压约为-9.2 MPa，这充分说明了膨润土块体的渗透率极低，水分运移缓慢。

图  12  试验台架内部不同时刻水压分布云图（仅显示膨润土块体和破碎颗粒）

Figure  12.  Porewater pressure distribution at different times (Only the domains representing compacted bentonite blocks and crushed bentonite pellets are presented here)

下载: 全尺寸图片 幻灯片

4.3   应力场模拟结果

压实膨润土块体中应力的变化是膨润土遇水膨胀和加热器的热效应双重作用引起的，是一个相对复杂的演化过程。本研究假设了总应力在各个方向上是一样的，忽略了各向异性特征，图 13则为试验台架不同区域特征点总应力随时间的变化结果，其中图 13（a）试验数据测量的总应力是环向的，图 13（b）~（f）中测量的是径向的。从模拟结果可以看出，随着时间的推移，注水量的增多，总应力整体趋势是逐渐增大的，模拟结果整体上能够反映试验过程总应力的演化趋势。其次，对于膨润土中间部位，如图 13（d），（e）所示，受温度与水分运移速率的影响，这两个特征点的应力在试验开始增长缓慢，可能受加热的影响，且膨润土渗透系数极低，水分迁移较为缓慢，吸力逐渐减小，然后土体才开始发生明显的膨胀现象。但是对于试验过程中部分特征点的总应力出现两次先下降后升高（主要包括92~1176 d和1200~1500 d这两个试验段），其原因同样是由于传感器与压实膨润土块体安装过程中以及不同膨润土块体间存在缝隙，使得水分快速运移到这些特征点处，总应力升高，后续由于停止注水，水压降低，相应的总应力出现降低的情形，结果如图 13（a），（c）所示。对于数值模拟而言，由于现阶段很难将这些缝隙在模型中进行刻画，使得模拟的总应力试验开始阶段与实测数据的吻合程度相对较好，后续由于水分沿着缝隙快速迁移，模拟的总应力低于实测数据值，但模拟结果整体可以反映膨润土块体的膨胀演化趋势，该认识可为处置库条件下缓冲材料长期耦合性能评价提供参考依据。

图  13  试验台架膨润土内部不同区域特征点总应力模拟结果

Figure  13.  Comparison between predictive total stress and experimental results at different locations of China-Mock-up facility

下载: 全尺寸图片 幻灯片

5.   结论及建议

（1）综合考虑China-Mock-up试验系统内部5种不同材料属性，建立了工程尺度条件下高庙子膨润土热-水-力多场耦合数值模拟方法，能够准确模拟预测了近5 a试验台架内部不同特征点温度随时间演化规律。

（2）模拟结果表明经过近5 a运行试验台架内部膨润土尚未达到饱和状态，靠近加热器的特征点，相对湿度随时间表现出先干燥后饱和的趋势，远离加热器的特征点，相对湿度随时间逐渐增大，可较好地反映试验台架内部不同特征点相对湿度的整体行为演化趋势。

（3）模拟的总应力与实测试验数据间表现较好的一致性，整体上能够反映试验台架内部不同特征点的总应力随时间逐渐增大这一演化趋势，该认识可为处置库环境下缓冲材料多场耦合性能提供参考依据。

（4）针对膨润土块、颗粒之间缝隙，建议后续开展相关室内试验，分析土体在膨胀过程中缝隙处渗透性、变形规律，建立描述膨润土块间缝隙愈合机理的本构方程，提出相应的缝隙接触面数值模拟方法。其次，考虑到候选场址围岩和地下水的影响，建议后续在处置深度开展多场耦合条件下现场原位试验研究，为研究真实处置库环境下缓冲材料行为演化特征提供技术支撑。