考虑多尺度孔隙耦合效应的非饱和膨胀土本构模型

摘要

颗粒型膨润土材料具有颗粒间孔隙、集合体间孔隙和集合体内孔隙等多尺度孔隙结构特征，传统非饱和膨胀土双孔模型（BExM）难以准确描述其宏观力学行为。在现有BExM模型的理论框架基础上，从颗粒型膨润土材料的三重孔隙结构出发，通过双加载屈服面方程和宏微观耦合效应函数表征膨润土颗粒混合物的多尺度孔隙弹塑性变形行为及其相互耦合关系。通过对比试验数据与模型模拟结果发现，模型可较好地描述非饱和颗粒型膨润土材料的膨胀力特征、膨胀变形特征以及宏微观孔隙结构演化规律。

关键词

膨润土; 非饱和土; 本构模型; 多尺度孔隙结构

高放射性核废物（简称“高放废物”）具有放射性核素活度浓度高、释热量大、毒性大以及半衰期长等特点，如何安全、有效、永久性处置高放废物是当前各有核国家面临的世界性难题^［1］。目前，世界各国普遍认为最可行的处置方案是“深地质处置”，即将高放废物埋藏至距地表以下500~1 000 m的稳定地质体中，使之与人类生存环境永久隔离^［2-3］。高放废物深地质处置通常采用“多重屏障系统”的设计理念，即将盛有核废物的处置容器安置于深层开挖的围岩“洞穴”中，并在容器与围岩间填充缓冲/回填材料^［4-5］。作为处置库中最为重要的人工屏障，缓冲/回填材料承担着核素阻滞屏障（阻滞核素迁移）、水力学屏障（延缓地下水入渗和维护废物罐结构稳定）和热传导屏障（有效传递核素衰变热）等多种屏障作用，对处置库的长期安全至关重要。大量研究表明，压实膨润土材料具有高膨胀性、低渗透性、强核素吸附性以及良好的热传导能力，被认为是最合适的缓冲/回填材料^［6-7］。根据处置库工程屏障系统的设计，缓冲/回填材料通常以压实膨润土块体的形式填充于废物罐与围岩之间。然而，近些年的工程实践表明，膨润土块体在砌置过程中，将产生多种施工缝隙，包括块体与围岩间缝隙、块体与废物罐间缝隙及块体与块体间缝隙等，这些施工缝隙的存在将削弱膨润土块体的缓冲性能，进而影响处置库的长期安全^［8］。为此，相关学者提出采用颗粒型膨润土材料填堵施工缝隙，以提高屏障系统的整体性能。研究表明，由于膨润土颗粒具有制备简单、安装操作方便、可有效填充不规则缝隙等优点，颗粒型膨润土已成为处置库缓冲/回填材料的重要组成部分^［9-10］。

膨润土类屏障材料的力学性质是衡量和评价工程屏障缓冲性能的关键指标之一。为此，国内外学者对压实膨润土材料的力学特性开展了大量研究工作，如叶为民等^［11］、Villar等^［12］、Cui^［13］。然而，现有研究多关注于压实膨润土块体，而对颗粒型膨润土材料的研究报道较少。近年来，随着颗粒型膨润土材料在工程实践中的广泛应用，颗粒型膨润土材料的力学特性研究逐渐引起国内外学者的关注^［14-16］。研究表明，压实膨润土块体通常呈现双重孔隙结构特征，即集合体间孔隙（inter-aggregate pore）和集合体内孔隙（intra-aggregate pore），其宏观力学性质受双重孔隙结构及其耦合效应共同影响^［17-18］。为了描述具有双孔结构的非饱和膨胀性土体的体变特征，Alonso等^［17］从宏微观层面出发，通过耦合效应函数考虑微观结构胀缩作用对宏观孔隙结构的影响，建立了非饱和膨胀土本构模型（BExM）。研究发现，BExM模型可以较好地反映非饱和压实膨润土的宏微观变形特征^［19-20］。然而，同压实膨润土块体相比，颗粒型膨润土材料的力学性质更为复杂，这主要与其多尺度孔隙结构有关。Hoffmann等^［15］指出，颗粒型膨润土材料内部不仅包含集合体间孔隙和集合体内孔隙，还存在更大尺寸的颗粒间孔隙。因此，颗粒混合物可近似看作三重孔隙结构，即颗粒间孔隙、颗粒内集合体间孔隙和集合体内孔隙^［21-23］。而在水化或加载过程中，混合物的力学性质受多尺度孔隙演化及其耦合作用的共同控制。Zhang等^［24］研究指出，非饱和膨润土颗粒混合物的变形特征涉及颗粒胀缩变形对粒间孔隙、集合体胀缩变形对集合体间孔隙等宏微观耦合效应的共同影响。因此，如何综合考虑颗粒型膨润土材料的多尺度孔隙结构特征及其宏微观耦合效应关系，是准确表征土体宏观力学行为的关键。

本文基于非饱和颗粒型膨润土材料的多尺度孔隙特征，发展和建立了考虑多尺度孔隙耦合效应的非饱和膨胀土本构模型。不同于传统的非饱和土双孔模型，该模型将土体分为三重孔隙，提出采用两个加载坍塌屈服面方程来描述不同尺度孔隙结构的弹塑性变形特征。同时，模型通过建立宏观和微观耦合效应函数以表征膨润土颗粒胀缩变形对粒间孔隙、集合体胀缩变形对集合体间孔隙的影响。最后，利用该模型对高庙子（GMZ）膨润土颗粒混合物的力学性质进行模拟与分析，结果表明，本文模型能较好地反映非饱和颗粒型膨润土材料的膨胀力特征、膨胀变形特征以及孔隙结构演化规律。

1 模型理论推导

颗粒型膨润土材料可近似分为颗粒间孔隙（大孔）和颗粒内孔隙。其中，颗粒内孔隙又可细分为集合体间孔隙（中孔）和集合体内孔隙（小孔）。因此，混合物的变形特征由3种孔隙的变化共同控制。

膨润土颗粒混合物的总孔隙比e可写为

$$e=e_{\mathrm{p}\mathrm{e}\mathrm{l}}+e_{\mathrm{m}\mathrm{a}\mathrm{c}}$$

（1）

$$e_{\mathrm{p}\mathrm{e}\mathrm{l}}=e_{\mathrm{m}\mathrm{e}\mathrm{s}}+e_{\mathrm{m}\mathrm{i}\mathrm{c}}$$

（2）

式中：e_pel为膨润土颗粒孔隙比；e_mac为颗粒间孔隙对应孔隙比（大孔孔隙比）；e_mes为中孔孔隙比；e_mic为小孔孔隙比。相应地，混合物的总体变增量dε_v可表示为

$$\mathrm{d}{\epsilon }_{\mathrm{v}}=\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{p}\mathrm{e}\mathrm{l}}+\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{m}\mathrm{a}\mathrm{c}}=\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{m}\mathrm{i}\mathrm{c}}+\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{m}\mathrm{e}\mathrm{s}}+\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{m}\mathrm{a}\mathrm{c}}$$

（3）

式中：dε_v-pel为颗粒自身的体变；dε_v-mic、dε_v-mes和dε_v-mac分别对应小孔、中孔和大孔的体变。

在水‒力耦合条件下，土体的变形主要受吸力变化和应力变化影响，混合物的总体变增量dε_v也可写为

$$\mathrm{d}{\epsilon }_{\mathrm{v}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}$$

（4）

式中：dε_vs为吸力变化引起的总体变增量；dε_vp为应力变化引起的总体变增量。

在吸力变化或应力变化过程中，颗粒混合物均会产生弹塑性变形，以下分别描述吸力变化和应力变化对土体体变的影响。

1.1　吸力变化条件下土体的变形计算

吸湿或脱湿过程中，膨润土颗粒混合物的总体变增量dε_vs可分为弹性体变增量$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{e}}$和塑性体变增量$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{p}}$两部分，即

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{p}}$$

（5）

基于颗粒混合物的三重孔隙结构，$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{e}}$和$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{p}}$可分别表示为

$$\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{e}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}+\\ \mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}\end{array}$$

（6）

$$\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{p}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{p}}+\\ \mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}\end{array}$$

（7）

式中： $\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{e}}$和$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}$分别为膨润土颗粒的弹性体变增量和塑性体变增量；$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}$、$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{e}}\mathrm{和}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}$分别为小孔、中孔和大孔对应的弹性体变增量；$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{p}}$、$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}\mathrm{和}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}$分别为小孔、中孔和大孔对应的塑性体变增量。

参照Alonso等^［17］提出的假设，模型假设微观结构（小孔）为饱和弹性体且遵循有效应力原理，微观结构的变形受有效应力$\widehat{p}$控制，$\widehat{p}$=p+$\mathrm{\aleph }$s（ p为平均净应力，s为吸力，$\mathrm{\aleph }$为饱和度S_r的函数）。当S_r=100%时，$\mathrm{\aleph }$=1，即$\widehat{p}$=p+s。同时，模型采用双应力状态变量（p和 s）描述中孔和大孔的体变行为。参照BExM模型，各级孔隙的弹性体变增量可表示为

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}=\frac{\mathrm{d}{e}_{\mathrm{m}\mathrm{i}\mathrm{c}}}{\mathrm{1}+e}=\frac{\mathrm{d}\widehat{p}}{K_{\mathrm{m}\mathrm{i}\mathrm{c}}}$$ ，$$K_{\mathrm{m}\mathrm{i}\mathrm{c}}=\frac{\left(\mathrm{1}+e\right)\widehat{p}}{{\kappa }_{\mathrm{m}\mathrm{i}\mathrm{c}}}$$

（8）

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{e}}=\frac{\mathrm{d}{e}_{\mathrm{m}\mathrm{e}\mathrm{s}}}{\mathrm{1}+e}=\frac{\mathrm{d}s}{K_{\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}}$$ ，$$K_{\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}=\frac{\left(\mathrm{1}+e\right)s}{{\kappa }_{\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}}$$

（9）

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}=\frac{\mathrm{d}{e}_{\mathrm{m}\mathrm{a}\mathrm{c}}}{\mathrm{1}+e}=\frac{\mathrm{d}s}{K_{\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}}$$ ，$$K_{\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}=\frac{\left(\mathrm{1}+e\right)s}{{\kappa }_{\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}}$$

（10）

式（8）~（10）中：κ_mic为小孔的弹性体变系数；κ_s-mes和κ_s-mac分别为吸力变化影响下中孔和大孔的弹性体变系数。

根据假设，微观结构为饱和弹性体，因此 $\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{p}}$= 0。这意味着吸力变化引起的膨润土颗粒塑性体变$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}$均为中孔的塑性体变$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}$，即$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}$ = $\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}$。为了描述吸力变化过程中微观结构的变化对中孔结构的影响，模型借鉴Alonso等^{［17，19］}提出的计算方法，采用微观耦合效应函数f_Dm或f_Im表征中孔塑性体变增量与小孔体变增量间的关系，即

当吸力降低时，

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}=f_{\mathrm{D}\mathrm{m}}\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}=\frac{f_{\mathrm{D}\mathrm{m}}\mathrm{d}\widehat{p}}{K_{\mathrm{m}\mathrm{i}\mathrm{c}}}$$

（11）

当吸力增高时，

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}=f_{\mathrm{I}\mathrm{m}}\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}=\frac{f_{\mathrm{I}\mathrm{m}}\mathrm{d}\widehat{p}}{K_{\mathrm{m}\mathrm{i}\mathrm{c}}}$$

（12）

微观耦合效应函数f_Dm和f_Im与土体的应力状态有关^［16］，可表示为p/p_0m的函数（p_0m为膨润土颗粒的先期固结压力）。为避免微观耦合效应函数在p/p_0m较高或较低时出现不收敛问题，本文提出下列关系式：

$$f_{\mathrm{D}\mathrm{m}}=a_{\mathrm{D}\mathrm{m}}/(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(b_{\mathrm{D}\mathrm{m}}+c_{\mathrm{D}\mathrm{m}}(\frac{p}{p_{\mathrm{0}\mathrm{m}}})))$$

（13）

$$f_{\mathrm{I}\mathrm{m}}=a_{\mathrm{I}\mathrm{m}}/(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(b_{\mathrm{I}\mathrm{m}}+c_{\mathrm{I}\mathrm{m}}(\frac{p}{p_{\mathrm{0}\mathrm{m}}})))$$

（14）

式（13）、（14）中：$a_{\mathrm{D}\mathrm{m}}$、$b_{\mathrm{D}\mathrm{m}}$、$c_{\mathrm{D}\mathrm{m}}$、$a_{\mathrm{I}\mathrm{m}}$、$b_{\mathrm{I}\mathrm{m}}$和$c_{\mathrm{I}\mathrm{m}}$均为参数。

在宏观层面上，膨润土颗粒的变形将会对颗粒间孔隙产生影响，引起宏观孔隙发生塑性变形。参照上述方法，模型提出宏观耦合效应函数f_DM或f_IM来表征宏观孔隙塑性体变增量与颗粒体变增量间的关系，即

当吸力降低时，

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}=f_{\mathrm{D}\mathrm{M}}\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{p}\mathrm{e}\mathrm{l}}$$

（15）

当吸力增高时，

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}=f_{\mathrm{I}\mathrm{M}}\mathrm{d}{\epsilon }_{\mathrm{v}-\mathrm{p}\mathrm{e}\mathrm{l}}$$

（16）

式（15）、（16）中：f_DM、f_IM与混合物的应力状态有关，可表示为p/p_0M的函数（p_0M为膨润土颗粒混合物的先期固结压力），二者间的关系式为

$$f_{\mathrm{D}\mathrm{M}}=a_{\mathrm{D}\mathrm{M}}/(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(b_{\mathrm{D}\mathrm{M}}+c_{\mathrm{D}\mathrm{M}}(\frac{p}{p_{\mathrm{0}\mathrm{M}}})))$$

（17）

$$f_{\mathrm{I}\mathrm{M}}=a_{\mathrm{I}\mathrm{M}}/(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(b_{\mathrm{I}\mathrm{M}}+c_{\mathrm{I}\mathrm{M}}(\frac{p}{p_{\mathrm{0}\mathrm{M}}})))$$

（18）

式（17）、（18）中：$a_{\mathrm{D}\mathrm{M}}$、$b_{\mathrm{D}\mathrm{M}}$、$c_{\mathrm{D}\mathrm{M}}$、$a_{\mathrm{I}\mathrm{M}}$、$b_{\mathrm{I}\mathrm{M}}$和$c_{\mathrm{I}\mathrm{M}}$均为参数。

1.2　应力变化条件下土体的变形计算

在加载或卸载过程中，膨润土颗粒混合物的体变增量dε_vp可分为弹性体变增量$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{e}}$和塑性体变增量$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{p}}$两部分，即

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{p}}$$

（19）

应力变化引起的各级孔隙弹塑性体变增量可分别表示为

$$\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{e}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}+\\ \mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}\end{array}$$

（20）

$$\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{p}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}=\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{p}}+\\ \mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}\end{array}$$

（21）

式（20）、（21）中：$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{e}}$和$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}$ 分别为应力变化引起的膨润土颗粒弹性体变增量和塑性体变增量； $\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}$、$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{e}}\mathrm{和}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}$分别为应力变化引起的小孔、中孔和大孔的弹性体变增量； $\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{p}}$、$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{p}}\mathrm{和}\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{p}}$分别为应力变化引起的小孔、中孔和大孔的塑性体变增量。

小孔、中孔和大孔的弹性体变增量可写为

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{i}\mathrm{c}}^{\mathrm{e}}=\frac{\mathrm{d}{e}_{\mathrm{m}\mathrm{i}\mathrm{c}}}{\mathrm{1}+e}=\frac{\mathrm{d}\widehat{p}}{K_{\mathrm{m}\mathrm{i}\mathrm{c}}}$$ ， $$K_{\mathrm{m}\mathrm{i}\mathrm{c}}=\frac{\left(\mathrm{1}+e\right)\widehat{p}}{{\kappa }_{\mathrm{m}\mathrm{i}\mathrm{c}}}$$

（22）

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}^{\mathrm{e}}=\frac{\mathrm{d}{e}_{\mathrm{m}\mathrm{e}\mathrm{s}}}{\mathrm{1}+e}=\frac{\mathrm{d}p}{K_{\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}}$$ ， $$K_{\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}=\frac{\left(\mathrm{1}+e\right)p}{{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}}$$

（23）

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}^{\mathrm{e}}=\frac{\mathrm{d}{e}_{\mathrm{m}\mathrm{a}\mathrm{c}}}{\mathrm{1}+e}=\frac{\mathrm{d}p}{K_{\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}}$$ ， $$K_{\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}=\frac{\left(\mathrm{1}+e\right)p}{{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}}$$

（24）

当描述应力变化引起的土体塑性变形时，首先，需定义土体的屈服面。参照Alonso等^［25］的方法，模型采用加载坍塌屈服面（LC）方程来表征加载引起的塑性变形，即

$$p_{\mathrm{0}\mathrm{M}}=p_{\mathrm{c}\mathrm{M}}{(\frac{p_{\mathrm{0}\mathrm{M}}^{*}}{p_{\mathrm{c}\mathrm{M}}})}^{\frac{(\lambda \left(\mathrm{0}\right)-{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}})}{(\lambda \left(s\right)-\kappa )}}$$

（25）

$$\lambda (s)=\lambda \left(\mathrm{0}\right)[r+\left(\mathrm{1}-r\right)\mathrm{e}\mathrm{x}\mathrm{p}(-\beta s)]$$

（26）

式（25）、（26）中：p_0M为混合物的先期固结压力；p^*_0M为饱和状态下混合物的先期固结压力；p_cM为参考压力；λ（s）为混合物e‒ln p平面下压缩曲线的斜率； λ（0）为饱和状态下混合物e‒ln p压缩曲线的斜率；r和β分别为试验参数。根据上述屈服面方程，可得加载引起的塑性体变为

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{p}}=\frac{(\lambda \left(s\right)-{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}})\mathrm{d}{p}_{\mathrm{0}\mathrm{M}}}{(\mathrm{1}+e)p_{\mathrm{0}\mathrm{M}}}$$

（27）

为了表征加载过程中膨润土颗粒的塑性变形，参照非饱和土BExM模型^［25］，提出LC'屈服面方程，即

$$p_{\mathrm{0}\mathrm{m}}=p_{\mathrm{c}\mathrm{m}}{(\frac{p_{\mathrm{0}\mathrm{m}}^{*}}{p_{\mathrm{c}\mathrm{m}}})}^{\frac{(\lambda {\left(\mathrm{0}\right)}_{\mathrm{m}}-{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}})}{(\lambda {\left(s\right)}_{\mathrm{m}}-{\kappa }_{\mathrm{m}})}}$$

（28）

$$\lambda {(s)}_{\mathrm{m}}=\lambda {(\mathrm{0})}_{\mathrm{m}}[r_{\mathrm{m}}+\left(\mathrm{1}-r_{\mathrm{m}}\right)\mathrm{e}\mathrm{x}\mathrm{p}(-{\beta }_{\mathrm{m}}s)]$$

（29）

式（28）、（29）中：p_0m为膨润土颗粒的先期固结压力；p^*_0m为饱和状态下颗粒的先期固结压力；p_cm为参考压力；λ（s）_m为e_pel‒ln p平面下颗粒压缩曲线的斜率； λ（0）_m为饱和状态下膨润土颗粒e_pel‒ln p压缩曲线的斜率；r_m和β_m分别为试验参数。

同式（27）相类似，应力变化引起的膨润土颗粒塑性体变可表示为

$$\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}=\frac{(\lambda {(s)}_{\mathrm{m}}-{\kappa }_{p-\mathrm{m}\mathrm{e}\mathrm{s}})\mathrm{d}{p}_{\mathrm{0}\mathrm{m}}}{(\mathrm{1}+e)p_{\mathrm{0}\mathrm{m}}}$$

（30）

此外，吸力变化和应力变化均会导致LC和LC'屈服面发生硬化，其硬化准则可表示为

$$\frac{\mathrm{d}{p}_{\mathrm{0}\mathrm{M}}}{p_{\mathrm{0}\mathrm{M}}}=\frac{(\mathrm{1}+e)(\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}}^{\mathrm{p}})}{\lambda \left(s\right)-{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{a}\mathrm{c}}}$$

（31）

$$\frac{\mathrm{d}{p}_{\mathrm{0}\mathrm{m}}}{p_{\mathrm{0}\mathrm{m}}}=\frac{(\mathrm{1}+e_{\mathrm{p}\mathrm{e}\mathrm{l}})(\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{s}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{v}\mathrm{p}-\mathrm{p}\mathrm{e}\mathrm{l}}^{\mathrm{p}})}{\lambda {(s)}_{\mathrm{m}}-{\kappa }_{\mathrm{p}-\mathrm{m}\mathrm{e}\mathrm{s}}}$$

（32）

2 模型预测与验证

2.1　参数获取方法

基于上述理论推导，模型共涉及27个参数，见表1。

表1  模型参数

Tab. 1  Parameters of model

孔隙类型	体变系数	屈服方程参数	耦合效应参数
小孔	κmic
中孔	κp-mes、κs-mes	p*0m、pcm、λ（0）m、rm、βm	aIm、bIm、cIm、aDm、bDm、cDm
大孔	κp-mac、κs-mac	p*0M、pcM、λ（0）、r、β	aIM、bIM、cIM、aDM、bDM、cDM

2.1.1　参数κ_mic

通过开展自由膨胀条件下膨润土颗粒的压汞（MIP）试验，绘制不同吸力下e_mic ‒ ln s关系曲线，拟合求取参数κ_mic。

2.1.2　参数κ_p-mes、κ_s-mes、κ_p-mac和κ_s-mac

参数κ_p-mes表示应力变化对集合体间孔隙（中孔）的影响，为e_mes‒ln p曲线弹性段斜率。参数κ_s-mes表示吸力变化对集合体间孔隙（中孔）的影响，为 e_mes‒ln s曲线弹性段斜率。类似地，参数κ_p-mac和κ_s-mac分别表示应力变化和吸力变化对混合物颗粒间孔隙的影响，为e_mac‒ln p曲线和e_mac‒ln s曲线弹性段斜率。

2.1.3　屈服面方程参数

根据不同吸力下膨润土颗粒混合物的压缩曲线，获取先期固结压力p^*_0M和压缩指数λ（s）_M随吸力的变化关系，采用式（25）和式（26）拟合获取参数p^*_0M、p_cM、λ（0）、r和β。类似地，根据不同吸力下膨润土颗粒的压缩曲线，通过式（28）和式（29）拟合获取参数p^*_0m、p_cm、λ（0）_m、r_m和β_m。

2.1.4　耦合效应参数

耦合效应参数f_DM和f_IM（或f_Dm和f_Im）根据恒荷载下混合物（或膨润土颗粒）孔隙比‒吸力关系曲线计算获取。根据f_DM‒p/p_0M和f_IM‒p/p_0M（或f_Dm‒p/p_0m和 f_Im‒p/p_0m）关系曲线，采用式（13）、（14）（或式（17）、（18））拟合求取参数a_DM、b_DM、c_DM、a_IM、b_IM和c_IM（或a_Dm、b_Dm、c_Dm、a_Im、b_Im和c_Im）。

2.2　模型参数确定

本文以高庙子（GMZ）膨润土颗粒混合物为对象，通过现有试验数据获取相关参数，并对颗粒混合物的宏观力学性质进行模拟预测和比对分析。

膨润土颗粒采用块体破碎法制备，详细的颗粒制备方法见文献［9］。不同吸力条件下的GMZ膨润土颗粒的孔隙分布曲线如图1所示。根据Zhang等^［26］的大小孔划分方法，本文选取孔径d=200 nm作为集合体内孔隙（小孔）和集合体间孔隙（中孔）的划分界限。图1b绘制了小孔孔隙比e_mic随吸力s的变化曲线，根据试验结果，拟合求取小孔体变参数κ_mic为0.046，这一计算结果与Wang等^［27］通过干湿循环试验获取的κ_mic相接近。

图1  不同吸力条件下膨润土颗粒的孔隙尺寸分布曲线与微观孔隙演化规律

Fig. 1  Pore size distribution curves and microstructural evolutions of bentonite pellets at different suctions

根据Zhang等^［24］获取的不同吸力条件下膨润土颗粒混合物的压缩曲线，参数κ_p-mac取为 0.015。图2为混合物先期固结压力p^*_0M和压缩参数λ（s）_M随吸力的变化关系。通过式（25）和（26）拟合获取模型参数p^*_0M = 0.462、p_cM = 0.025、λ（0） = 0.165、r = 0.742和β = 0.252。此外，按照Alonso等^［19］的假设，参数κ_s-mac取为0。

图2  不同吸力条件下膨润土颗粒混合物的先期固结压力与压缩参数

Fig. 2  Pre-consolidation pressures and compression parameters of bentonite pellet mixtures at different suctions

本文采用不同吸力条件下高压实膨润土块体（初始状态与膨润土颗粒一致）的一维固结试验结果表征膨润土颗粒的压缩特征。图3为膨润土颗粒的先期固结压力p^*_0m和压缩参数λ（s）_m随吸力的变化关系。结合式（28）、（29），拟合获取模型参数p^*_0m = 0.435、p_cm = 0.183、λ（0）_m = 0.156、r_m = 0.267和β_m = 0.264。同时，根据卸载曲线可获取参数κ_p-mes为0.012。参考Alonso等^［17］，本文取κ_s-mes = 0.004。根据初始荷载0.1 MPa下颗粒混合物和膨润土颗粒的孔隙比‒吸力关系曲线，分别获取f_DM‒p/p_0M和f_Dm‒p/p_0m的变化关系，见图4。结合式（11）~（18），拟合获取宏、微观耦合效应参数a_DM = -1.281、b_DM = 2.302、c_DM = 10.500、a_Dm = 1.850、b_Dm =1.581和c_Dm = 18.021。

图3  不同吸力条件下膨润土颗粒的先期固结压力与压缩参数

Fig. 3  Pre-consolidation pressures and compression parameters of bentonite pellet mixtures at different suctions

图4  f_DM‒p/p_0M和f_Dm‒p/p_0m的变化关系

Fig. 4  f_DM‒p/p_0M versus f_Dm‒p/p_0m

2.3　模型模拟结果

2.3.1　膨胀力特征模拟

膨胀力模拟采用Wang等^［28］提出的反向压实法（back-compaction method），如图5所示。该模拟方法的计算流程为：① 设置初始荷载p_ini = 0.05 MPa，计算初始状态下各级孔隙的孔隙比和模型参数；② 给予较小的吸力步长Δs，控制试样在初始荷载下水化膨胀，计算各级孔隙的膨胀变形量；③ 在恒定吸力下对试样进行压缩，使其恢复至初始孔隙比，计算压缩所需荷载为Δp，此时试样的上覆荷载为p_ini+Δp；④ 再次给予较小的吸力步长Δs，使试样在p_ini+Δp荷载下水化膨胀，计算各级孔隙的膨胀变形量；⑤ 在恒定吸力下压缩试样至初始孔隙比，计算试样的上覆荷载值。重复上述步骤，直至吸力降至0.01 MPa，试样的上覆荷载即为土体的最终膨胀力。

图5  反向压实法计算原理图

Fig. 5  Schematic diagram of back-compaction calculation method

图6为不同吸力条件下膨润土颗粒混合物的膨胀力p_s试验结果与模型模拟结果。从图6可以看出，膨胀力的变化规律可近似分为3个阶段：在阶段Ⅰ（113.0 MPa→38.0 MPa），膨胀力随吸力的降低而显著增加；在阶段Ⅱ（38.0 MPa→4.2 MPa），膨胀力变化平缓，不再显著；阶段Ⅲ（s<4.2 MPa），膨胀力明显增加直至稳定。通过分析试验结果与模型模拟结果可以看出，该模型能较好地反映颗粒型膨润土材料在三阶段的演化规律。

图6  膨润土颗粒混合物的膨胀力模拟结果

Fig. 6  Simulation results of swelling pressure of bentonite pellet mixtures

2.3.2　膨胀变形特征模拟

图7绘制了不同上覆荷载条件下膨润土颗粒混合物的孔隙比‒吸力变化关系曲线。通过对比试验数据与模拟结果可以发现，本文模型能够较好地反映颗粒混合物的膨胀变形特征。同时，从小孔孔隙比e_mic、中孔孔隙比e_mes和大孔孔隙比e_mac的发展规律可以看出，e_mic和e_mes均随吸力的降低而增加，而e_mac则逐渐减小。这一现象表明，在恒定荷载条件下，膨润土颗粒湿化膨胀填充颗粒间孔隙，颗粒集合体内和集合体间孔隙体积增大，这与Hoffmann等^［15］所观察到的试验现象一致。

图7  不同上覆应力下膨润土颗粒混合物的膨胀变形模拟结果

Fig. 7  Simulation results of swelling deformation of bentonite pellet mixtures at different vertical stresses

2.3.3　宏微观孔隙演化特征

图8绘制了恒体积水化过程中膨润土颗粒混合物大孔孔隙比e_mac和颗粒孔隙比e_pel随吸力的变化曲线。通过试验数据与模拟结果对比发现，模型能够较好地反映恒体积条件下混合物孔隙结构的演化规律。模拟结果表明，膨润土颗粒的膨胀变形由集合体膨胀和集合体间孔隙膨胀承担。通过对比图6中膨胀力的发展曲线可发现，水化过程中膨润土颗粒孔隙比与宏观膨胀力的变化近似呈相反规律，即在高吸力（s > 38.0 MPa）和低吸力（s < 4.2 MPa）阶段，颗粒的填充作用不明显但膨胀力显著增加；而在4.2 MPa < s <38.0 MPa吸力范围，膨润土颗粒的膨胀填充效应明显但膨胀力变化平缓。究其原因，在恒体积条件下，混合物的水化膨胀主要表现在填充粒间孔隙和促进膨胀力增长两方面。不同水化程度下颗粒膨胀的宏观表现受颗粒强度和粒间孔隙填充程度的影响不同。在高吸力阶段，颗粒自身强度较大且集合体联结紧密，颗粒膨胀引起膨胀力显著增加；而随着吸力的降低，颗粒自身膨胀软化，混合物骨架结构坍塌，膨胀主要填充颗粒间孔隙，膨胀力变化不明显；但当吸力较低时，粒间孔隙填充基本完成，颗粒膨胀将再次导致膨胀力显著增长。

图8  恒体积水化过程中膨润土颗粒混合物孔隙结构演化规律

Fig. 8  Evolutions of pore structure inside bentonite pellet mixtures under constant-volume condition

3 结论

本文基于颗粒型膨润土材料的多尺度孔隙结构特征及其耦合变形机制，建立了考虑多尺度孔隙耦合效应的非饱和膨胀土本构模型。以膨润土颗粒混合物为研究对象，模型将颗粒混合物分为三重孔隙，采用双加载坍塌屈服面表征颗粒间孔隙和颗粒自身孔隙的弹塑性变形特征，通过引入宏、微观耦合效应函数来描述颗粒变形对颗粒间孔隙、集合体变形对颗粒自身孔隙的影响，进而推导了各向同性应力状态下本构方程表达式。通过对比试验数据与模拟结果可发现，该模型能够较好地描述非饱和颗粒型膨润土材料的膨胀力特征、膨胀变形特征以及孔隙演化规律。然而，由于现有试验成果有限，本文模型仅考虑各向同性应力状态，尚未将其扩展至三轴应力状态。同时，该模型将膨润土颗粒混合物看作连续介质材料，并未考虑高吸力段混合物的颗粒型结构，因此，如何考虑水化过程中混合物由颗粒型结构向连续型结构的过渡是下一步理论模型研究的重点。

倪雪倩：结果分析与模型构建。

张 升：研究思路与写作指导。参考文献：

［1］ 王驹. 高水平放射性废物地质处置： 关键科学问题和相关进展 ［J］. 科技导报， 2016， 34（15）： 51.

WANG Ju. Geological disposal of high level radioactive waste： key scientific issues and progress in China ［J］. Science & Technology Review， 2016， 34（15）： 51.

［2］ 潘自强， 钱七虎. 我国高放废物地质处置战略研究 ［J］. 中国核电， 2013， 6（2）： 98.

PAN Ziqiang， QIAN Qihu. The geological disposal of high-level radioactive waste strategy research in our country ［J］. China Nuclear Power， 2013， 6（2）： 98.

［3］ 李腾， 周小毛， 王玲钰， 等. 乏燃料处置的必要性及其处置库环境化学行为 ［J］. 核化学与放射化学， 2020， 42（6）： 513.

LI Teng， ZHOU Xiaomao， WANG Lingyu， et al. Spent fuel disposal and chemical behaviors at repository environment ［J］. Journal of Nuclear and Radiochemistry， 2020， 42（6）： 513.

［4］ 王驹， 凌辉， 陈伟明. 高放废物地质处置库安全特性研究 ［J］. 中国核电， 2017， 10（2）： 270.

WANG Ju， LING Hui， CHEN Weiming. Study on the safety functions of repository for geological disposal of high level radioactive waste ［J］. China Nuclear Power， 2017， 10（2）： 270.

［5］ CHEN L， LIU Y M， WANG J， et al. Investigation of the thermal-hydro-mechanical （THM） behavior of GMZ bentonite in the China-mock-up test ［J］. Engineering Geology， 2014， 172（7）： 57.

［6］ 叶为民， 钱丽鑫， 陈宝， 等. 侧限状态下高压实高庙子膨润土非饱和渗透性的试验研究 ［J］. 岩土工程学报， 2009， 31（1）： 105.

YE Weimin， QIAN Lixin， CHEN Bao， et al. Laboratory test on unsaturated hydraulic conductivity of densely compacted Gaomiaozi bentonite under confined conditions ［J］. Chinese Journal of Geotechnical Engineering， 2009， 31（1）： 105.

［7］ 刘月妙， 徐国庆， 刘淑芬， 等. 我国高放废物处置库缓冲/回填材料压实膨胀特性研究 ［J］. 铀矿地质， 2001， 17（1）： 44.

LIU Yuemiao， XU Guoqing， LIU Shufen， et al. Study on compactibility and swelling property of buffer/backfill material for HLW repository ［J］. Uranium Geology， 2001， 17（1）： 44.

［8］ 陈永贵， 贾灵艳， 叶为民， 等. 施工接缝对缓冲材料水‒力特性影响研究进展［J］. 岩土工程学报， 2017， 39（1）： 138.

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.

［9］ ZHANG Z， YE W M， LIU Z R， et al. Influences of PSD curve and vibration on the packing dry density of crushed bentonite pellet mixtures ［J］. Construction and Building Materials， 2018， 185： 246.

［10］ LIU Z R， YE W M， ZHANG Z， et al. Particle size ratio and distribution effects on packing behaviour of crushed GMZ bentonite pellets ［J］. Powder Technology， 2019， 351： 92.

［11］ 叶为民， SCHANZ T， 钱丽鑫， 等. 高压实高庙子膨润土 GMZ01 的膨胀力特征 ［J］. 岩石力学与工程学报， 2007， 26（S2）： 3861.

YE Weimin， SCHANZ T， QIAN Lixin， et al. Characteristics of swelling pressure of densely compacted Gaomiaozi bentonite GMZ01 ［J］. Chinese Journal of Rock Mechanics and Engineering， 2007， 26（S2）： 3861.

［12］ VILLAR M V， LLORET A. Influence of temperature on the hydro-mechanical behaviour of a compacted bentonite ［J］. Applied Clay Science， 2004， 26（1）： 337.

［13］ CUI Y J. On the hydro-mechanical behaviour of MX80 bentonite-based materials ［J］. Journal of Rock Mechanics and Geotechnical Engineering， 2017， 9（3）： 187.

［14］ IMBERT C， VILLAR M V. Hydro-mechanical response of a bentonite pellets/powder mixture upon infiltration ［J］. Applied Clay Science， 2006， 32（3）： 197.

［15］ HOFFMANN C， ALONSO E E， ROMERO E. Hydro-mechanical behaviour of bentonite pellet mixtures ［J］. Physics and Chemistry of the Earth， 2007， 32（8）： 832.

［16］ 张虎元， 王学文， 刘平， 等. 缓冲回填材料砌块接缝密封及愈合研究 ［J］. 岩石力学与工程学报， 2016， 35（S2）： 3605.

ZHANG Huyuan， WANG Xuewen， LIU Ping， et al. Sealing and healing of compacted bentonite block joints in HLW disposal ［J］. Chinese Journal of Rock Mechanics and Engineering， 2016， 35（S2）： 3605.

［17］ ALONSO E E， VAUNAT J， GENS A. Modelling the mechanical behaviour of expansive clays ［J］. Engineering Geology， 1999， 54（1）： 173.

［18］ SÁNCHEZ M， GENS A， GUIMARÃES L N， et al. A double structure generalized plasticity model for expansive materials ［J］. International Journal for Numerical and Analytical Methods in Geomechanics， 2005， 29（8）： 751.

［19］ ALONSO E E， ROMERO E， HOFFMANN C， et al. Expansive bentonite–sand mixtures in cyclic controlled-suction drying and wetting ［J］. Engineering Geology， 2005， 81（3）： 213.

［20］ LLORET A， VILLAR M V， SÁNCHEZ M， et al. Mechanical behaviour of heavily compacted bentonite under high suction changes ［J］. Géotechnique， 2003， 53（1）： 27.

［21］ MOLINERO-GUERRA A， AIMEDIEU P， BORNERT M， et al. Analysis of the structural changes of a pellet/powder bentonite mixture upon wetting by X-ray computed microtomography ［J］. Applied Clay Science， 2018， 165： 164.

［22］ LIU Z R， YE W M， CUI Y J， et al. Insights into the water retention behaviour of GMZ bentonite pellet mixture ［J］. Acta Geotechnica， 2021， 16： 3145.

［23］ ZHANG Z， YE W M， LIU Z R， et al. Mechanical behavior of GMZ bentonite pellet mixtures over a wide suction range ［J］. Engineering Geology， 2020， 264： 105383.

［24］ ZHANG Z， YE W M， WANG Q， et al. Investigation on healing behavior of unsaturated GMZ bentonite pellet mixture based on compressibility ［J］. Acta Geotechnica， 2022， 17（10）： 4461.

［25］ ALONSO E E， GENS A， JOSA A. Constitutive model for partially saturated soils ［J］. Géotechnique， 1990， 40（3）： 405.

［26］ ZHANG Z， CUI Y J， YANG J W， et al. Water retention and compression behavior of MX80 bentonite pellet ［J］. Acta Geotechnica， 2021， 17（6）： 2435.

［27］ WANG Y， YE W M， CHEN B， et al. Thermal-mechanical effects on volume-change behavior of compacted GMZ bentonite during cyclic wetting-drying processes ［J］. Environmental Earth Sciences， 2019， 78（17）： 1.

［28］ WANG Q， TANG A M， CUI Y J， et al. Investigation of the hydro-mechanical behaviour of compacted bentonite/sand mixture based on the BExM model ［J］. Computers and Geotechnics， 2013， 54： 46.

作者贡献声明

张 召：理论推导与论文撰写。