0.   引言

岩土工程中的许多工程问题都与非饱和土的水-力耦合特性密切相关,例如,非饱和边坡在降雨和地震作用下的垮塌、非饱和路基的不均匀沉降、非饱和膨胀土的开裂、垃圾填埋场中渗滤液的运移等。不同于饱和土（固、液两相体）,非饱和土（固、液、气三相体）颗粒间同时存在孔隙气和孔隙水,在水-气界面会形成月牙形的弯液面,这使得非饱和土颗粒之间会产生黏结效应,这也是非饱和土具有复杂水力-力学行为的根本原因。因此,建立一个充分考虑黏结效应的非饱和土本构模型,对于合理预测非饱和土的水力-力学行为以及分析与此相关的岩土工程问题具有重要意义。

在早期的研究中,非饱和土颗粒间的黏结效应被归因于基质吸力,其主要从两个方面影响着非饱和土的力学行为^[1]：①通过改变土孔隙中流体的平均压力进而改变土骨架应力^[2]；②通过弯液面在颗粒之间施加黏结力。因此,基质吸力在很长一段时间被作为描述非饱和土力学行为的基本变量,相关的非饱和土本构模型主要包括以下两类：①以Bishop^[3]有效应力作为单应力变量的本构模型；②以BBM模型^[4]为代表的采用净应力和基质吸力为双应力变量^[5]的本构模型。虽然基质吸力在非饱和土本构建模中具有重要地位,但在非饱和土气-液分界面弯液面的复杂作用下,颗粒间的黏结效应并不能简单利用基质吸力来解释。此后,很多学者认识到这一问题,开始更加全面地研究影响非饱和土黏结效应的因素,并针对这一概念进行了不同的定义。沈珠江等^[6]认为颗粒间的黏结效应归因于一切增加颗粒间抗滑阻力的因素,并定义其为广义吸力。汤连生等^[7]采用粒间吸力来定义非饱和土颗粒间的黏结效应,并将其分为湿吸力和结构吸力,认为当饱和度不同时,粒间吸力也会产生较大变化。Wheeler等^[8]发现,即使在相同的基质吸力作用下,非饱和土中弯液面的数量和大小不同（即饱和度不同）也会使得黏结效应发生显著改变。将饱和度对非饱和土颗粒间黏结效应的影响应用到土的本构关系中,建立了一个非饱和土弹塑性水-力耦合模型。此后,饱和度对非饱和土颗粒间黏结效应影响越来越受到学者们的关注,考虑饱和度影响的非饱和土弹塑性本构模型得到发展^[9-10]。为了更加直观反映非饱和土的黏结效应,Gallipoli等^[1]开创性地定义了一个全新的黏结变量ξ来作为本构变量,以此取代基质吸力,并建立了一个直接考虑黏结效应影响的非饱和土弹塑性本构模型Gallipoli-Gens-Sharma-Vaunat模型（GGSV模型）。该黏结变量同时考虑了颗粒间黏结力的大小以及弯液面数量的影响。近几年,Hu等^[11],Ma等^[12]也采用类似的思路建立了模型。

除了基质吸力和饱和度以外,孔隙结构对于非饱和土黏结效应的影响也不可忽视。土体内部不同的矿物成分以及制样方式均会对孔隙结构产生影响,通常认为压实黏土具有显著的双孔结构,大小孔隙的变化对非饱和土黏结效应贡献不同,Alonso等^[13]针对双孔结构土,建立了可区分大小孔隙不同贡献的有效饱和度表达式来反映微观孔隙结构的影响。蔡国庆等^[14]通过引入有效饱和度对Wheeler模型^[8]进行改进,建立了考虑微观孔隙结构影响的非饱和土水-力耦合弹塑性模型。陈正汉^[15]通过将细观结构数据与宏观反映数据相结合,建立了考虑细观结构影响的非饱和土本构模型。

综上所述,非饱和土颗粒间的黏结效应与基质吸力、饱和度、孔隙结构等均有着密切关系,现有的非饱和土本构模型大多以考虑基质吸力和饱和度的影响来间接反映黏结效应,且对孔隙结构的影响鲜有考虑。建立一个综合考虑基质吸力、饱和度与孔隙结构对黏结效应影响的非饱和土本构模型,对于更加准确合理预测非饱和土水-力耦合行为具有重要意义。此外,现有的大部分非饱和本构模型都是基于传统弹塑性理论对非饱和土的变形行为进行描述,近年来,边界面模型^[16]因其具有良好的预测能力和适应性得到了较为广泛的应用^[17-18],基于适合的映射准则可描述高度非线性的土体变形行为。对于水力-力学行为更为复杂的非饱和土,许多学者^[19-21]选择采用边界面模型进行预测并取得了良好效果。

基于GGSV模型,建立了一个考虑颗粒间黏结效应的非饱和土水-力耦合模型。采用全新的黏结变量和有效应力作为本构变量来同时考虑基质吸力、饱和度、孔隙结构对黏结效应的影响。结合考虑变形影响的土-水特征曲线模型,实现了对非饱和土水力滞后以及水力-力学耦合特性的描述。此外,引入修正剑桥模型的椭圆屈服面作为边界面,在边界面塑性理论框架内对土体塑性变形行为进行更加精准的预测。利用膨润土-高岭土混合土、重塑高岭土的试验结果,对模型进行了参数标定和模型验证,结果表明所建立模型能够合理预测非饱和土的水-力耦合特性。

1.   本构模型建立

1.1   本构变量

GGSV模型采用Bishop有效应力与黏结变量作为本构变量,有效应力中的参数$\chi$的取值为饱和度$S_{\text{r}}$,本构变量表达式为

式中,$s$为基质吸力,$S_{\text{r}}$为饱和度,$p_{\text{net}}$为净应力,$f(s)$为非饱和状态下单个弯液面对土颗粒施加的粒间应力与饱和状态下粒间应力的比值（土颗粒被假定为理想球体）,以此来反映单个弯液面对非饱和土颗粒间的黏结作用,粒间应力的大小可依据Fisher^[22]给出的解析解进行计算,Gallipoli等^[1]给出了$f(s)$与$s$的关系曲线,如图1所示,可依据图1获取不同基质吸力下$f(s)$的取值。$(1-S_{\text{r}})$代表非饱和土每单位体积固体中弯液面的数量^[1],二者的乘积即为黏结变量,反映非饱和土颗粒间的黏结效应,同时考虑了基质吸力和饱和度对黏结效应的影响。

图  1  $f(s)$与s关系曲线

Figure  1.  Curve of relation between$f(s)$ ands

下载: 全尺寸图片 幻灯片

GGSV模型定义的本构变量虽然能够较好地反映基质吸力与饱和度对黏结效应的影响,但并未考虑孔隙结构的影响。Alonso等^[13]将饱和度$S_{\text{r}}$分为宏观饱和度$S_{\text{r}}^{\text{M}}$和微观饱和度$S_{\text{r}}^{\text{m}}$以反映水对大、小孔隙的填充程度,并从非饱和土微观孔隙结构出发,引入微观饱和度定义了有效饱和度：

式中,$S_{\text{res}}$为残余饱和度,$S_{\text{r}}^{\text{m}}$为微观孔隙的饱和度。在数值上可认为等于残余饱和度。通过引入微观饱和度可以有效区分双孔结构土体中大、小孔隙的不同贡献,进而考虑孔隙结构的影响。

为了反映微观孔隙结构对非饱和土黏结效应的影响,引入有效饱和度替换饱和度来建立全新的黏结变量$\zeta$和有效应力$p^{″}$,以此做为模型的本构变量：

1.2   考虑黏结效应的屈服方程

Gallipoli等^[1]认为,在相同的有效应力作用下,非饱和状态下孔隙比与饱和状态下孔隙比的比值$e\text{/}{e}_{\text{s}}$是黏结变量的唯一函数,并利用多组试验数据^{[1, 23]}验证了这一观点,建立了如下关系式：

式中,$e$为非饱和状态下土体的孔隙比,$e_{\text{s}}$为饱和状态下土体的孔隙比,$\xi$为黏结变量,a,b为拟合参数。

本模型采用相同的思路,建立如下的关系式：

式中,$\zeta$为改进的黏结变量。

饱和状态下的孔隙比可依据$e_{\text{s}}$-$\ln p^{″}$空间内的饱和正常固结线求得：

式中,N为$p^{″}$=1 kPa对应的比体积,$\lambda$为$e_{\text{s}}$-$\ln p^{″}$空间内饱和正常压缩曲线的斜率。

依据式（7）可以得到非饱和状态下的孔隙比为

令$N(\zeta )=h(\zeta )N$,$\lambda (\zeta )=h(\zeta )\lambda$,可得

式（9）定义了在$e$-$\ln p^{″}$-$\zeta$空间的压缩状态面,如图2（a）所示,任意非饱和土的压缩曲线均在这一状态面上,其在$e$-$\ln p^{″}$空间的投影如图2（b）所示。考虑下列特定情况：在饱和状态下,由A点经历一个常孔隙比的增湿过程到达B′点,再由B′点经历一个弹性压缩过程到达B点,由图2（a）可知,从A点到B点的整个过程,只在A点到B′点产生弹性变形,由B′点到B点孔隙比不发生变化。从空间路径来看,从A点出发沿着曲线AB到达B点也产生同样的弹性变形。这样便形成了一个空间弹性面AB′B,其与压缩状态面的交线即为弹性区域的边界线,也就是屈服曲线。

图  2  压缩状态面及压缩曲线示意图

Figure  2.  Diagrams of compression state surface and compression curve

下载: 全尺寸图片 幻灯片

接下来推导考虑黏结效应屈服曲线的方程。首先,弹性变形增量可定义为

式中,$\kappa$为回弹曲线斜率。

依据式（12）,可计算从A点到B点的孔隙比变化为

式中,${p^{″}}_0$,${p^{″}}_{\zeta }$分别为饱和状态和非饱和状态下的前期固结应力。

依据式（10）可得A点到B点的孔隙比变化为

将式（13）、（14）联立可求得屈服方程为

由于屈服曲线的位置由饱和状态下的前期固结应力控制,故将饱和前期固结应力作为模型的硬化参数。

1.3   边界面

边界面模型允许在边界面内产生塑性变形,可根据设定适当的映射准则对土体的变形进行非线性描述,对于多种土样均表现出良好的适应性和预测能力。选取修正剑桥模型的椭圆屈服面作为边界面,将式（15）建立的屈服曲线作为$p^{″}$- $\zeta$平面边界线。黄茂松等^[19]也采用了类似的方式建立非饱和黄土的边界面模型,将巴塞罗那模型的LC屈服线作为$p_{\text{net}}$-s平面的边界线。

图3为边界面在$p^{″}$-q平面和$p^{″}$-q-$\zeta$空间的示意图。为简化计算,本模型将边界面中的纯弹区域省略为一点,与坐标原点重合,认为所有变形均为弹塑性,这样的简化对土体实际变形的影响是可以忽略的^[24]。边界面方程为

图  3  $p^{″}$-q-$\zeta$空间内的边界面

Figure  3.  Bounding surfaces in$p^{″}$-q-$\zeta$ space

下载: 全尺寸图片 幻灯片

式中,${\overline{p}}^{″}$,$\overline{q}$,${{\overline{p}}^{″}}_{\zeta }$分别为边界面上的改进有效应力、剪应力和屈服应力,$M$为临界状态线斜率。

初始边界面与$p^{″}$轴的交点为饱和前期固结应力${p^{″}}_0$,位于屈服曲线上,满足式（15）,在后续的发展过程中,边界面随着塑性变形的产生不断发生硬化扩张,与$p^{″}$轴的交点不断右移。

1.4   流动法则

本模型采用非相关联的流动法则,塑性势面方程采用与文献[4]一致的形式：

式中,$\alpha =\frac{\lambda M(M-9)(M-3)}{9(\lambda -\kappa )(6-M)}$。

塑性应变增量采用边界面塑性理论计算：

式中,$H_{\text{p}}$,${\overline{H}}_{\text{p}}$分别为当前应力点和边界面上映射点对应的塑性模量,$n_{\text{Q}}$,$n_{\text{F}}^{}$分别为塑性势面和边界面上的单位法向量,分别代表塑性流动方向和加载方向,${\sigma }^{″}$,$\overline{{\sigma }^{″}}$分别为当前应力状态变量和边界面上的应力状态变量,其中,${\sigma }^{″}={\text{(}{p}^{″},q\text{)}}^{\text{T}}$,$\overline{{\sigma }^{″}}={({\overline{p}}^{″},\overline{q})}^{\text{T}}$。

依据边界面方程,可推导塑性势面和边界面上的单位法向量为

1.5   硬化法则

由于边界面的大小由饱和前期固结应力来确定,因此硬化参数为饱和前期固结应力,其与塑性体变关系密切,直接影响后续边界面的演化规律,硬化法则可用如下方程描述：

边界面的一致性条件为

依据边界面的一致性条件结合流动法则可计算出边界面上映射点对应的塑性模量${\overline{H}}_{\text{p}}$,

式中,$A\text{=}2\overline{p}-{p^{″}}_{\zeta }$,$B=h(\zeta )\lambda -\kappa$,$C\text{=}N-\lambda \ln {p^{″}}_{\zeta }$。

1.6   映射准则

当前应力点的塑性模量与映射点的塑性模量满足映射法则,可依据映射法则计算当前应力点的塑性模量。Zienkiewicz等^[25]建立了一个形式简单、参数少且适用性较好的径向映射法则,本模型采用类似的方式建立了如下的径向映射法则：

式中,$\gamma$为映射指数。${\delta }_0$,$\delta$分别为映射中心到当前应力点与映射点的距离,如图4所示,本模型选取坐标原点作为映射中心。映射指数$\gamma$的表达式为

图  4  边界面映射法则

Figure  4.  Mapping rules of bounding surface

下载: 全尺寸图片 幻灯片

式中,${\gamma }_0$为初始值,$D$为比例常数（$D＞0$）,${\displaystyle \int \left|\text{d}{\epsilon }^{\text{p}}\right|}$为塑性应变增量的累积量。${\displaystyle \int \left|\text{d}{\epsilon }^{\text{p}}\right|}\text{=}0$时,$\gamma \text{=}{\gamma }_0$；${\displaystyle \int \left|\text{d}{\epsilon }^{\text{p}}\right|}\text{=}\infty$时,$\gamma \text{=0}$,$H_{\text{p}}={\overline{H}}_{\text{p}}$。

将依据式（23）计算出的当前应力点的塑性模量$H_{\text{p}}$代入式（18）即可得到塑性应变增量。

弹性应变增量可由胡克定律计算：

式中,$C^{\text{e}}$为弹性柔度矩阵,

其中,$K$,$G$分别为土的体积模量和剪切模量,

其中,$\nu$为泊松比。

1.7   持水特性

非饱和土的持水特性是指由基质吸力与应力变化引起的饱和度的变化,通常采用土-水特征曲线（SWCC）来建立基质吸力与饱和度之间的关系。土在变形过程中其孔隙和孔隙之间连通的部分会生尺寸变化,进而导致SWCC发生变化^[26],因此变形对SWCC的影响不可忽视。近年来,变形对SWCC的影响受到越来越多的关注,Gallipoli等^[26-27]、Hu等^[28]、蔡国庆等^[29]均建立了相关的土-水特征曲线模型。在van Genuchten模型^[30]基础上,Gallipoli等^[26]引入比体积以及与进气值相关的参数建立模型,本文基于该模型建立水力滞后土-水特征曲线模型,但采用有效饱和度表示为

式中,$m_1$,$m_2$,$m_3$,$m_4$为模型参数,其中,$m_1$控制进气值的大小。

有效饱和度的变化可由下式计算：

基于Gallipoli ^[27]的研究成果,模型中决定主吸湿和主脱湿曲线的参数主要为$m_1$,为了简化计算,认为主吸湿曲线和主脱湿曲线仅由参数$m_1$控制,其余参数相等,

式中,$m_{1\text{d}}$为控制主脱湿曲线进气值的参数,$m_{1\text{w}}$为控制主吸湿曲线进气值的参数。

Zhou等^[10]采用了一种类似边界扫描法则^[31]的方法确定扫描曲线,相较于Gallipoli^[27]与Hu等^[28]采用的线性扫描法则,该方法对扫描曲线的预测更加准确。本文采用此方法在$s-S_{\text{e}}$平面建立扫描曲线,其在湿化和干化路径下的表达式为

式中,S_w,S_d分别为与当前应力点同一有效饱和度下主吸湿曲线和主脱湿曲线上对应的基质吸力,η为控制扫描曲线形状的参数,可通过对扫描试验数据的拟合来确定。扫描曲线示意图如图5所示。

图  5  扫描曲线示意图

Figure  5.  Diagram of scanning curve

下载: 全尺寸图片 幻灯片

干化与湿化扫描线上有效饱和度增量形式为

将考虑变形的水力滞后土-水特征曲线模型与边界面模型结合,即形成完整的考虑颗粒间黏结效应的水-力耦合边界面模型。

2.   模型参数及确定方法

三维状态下的本构模型共有16个参数：$a$,$b$,$\lambda$,$\kappa$,$N$,${p^{″}}_0$,$M$,$G$,$m_{1\text{w}}$,$m_{1\text{d}}$,$m_2$,$m_3$,$m_4$,${\gamma }_0$,$D$,$\eta$。其中,a,b可以采用常基质吸力下的等向压缩试验数据对$\zeta$与$e\text{/}{e}_{\text{s}}$的关系进行拟合获得,$\lambda$,$\kappa$,$N$,${p^{″}}_0$可通过饱和各向同性压缩试验获得,$M$,$G$可通过排水剪切试验确定,$m_{1\text{w}}$,$m_2$,$m_3$,$m_4$可通过湿化路径下常基质吸力等向压缩试验中饱和度与孔隙比的关系拟合得到^[32],$m_{1\text{d}}$,$\eta$通过主脱湿试验和扫描试验进行确定。在常基质吸力下,参数$m_{1\text{w}}$,$m_{1\text{d}}$,由$m_1$取代,参数$m_1$,$m_2$,$m_3$,$m_4$由常基质吸力下等向压缩试验中有效饱和度与比体积的关系拟合得到,${\gamma }_0$,$D$为材料常数,可根据土体的物理性质确定。此外,初始孔隙比可依据式（12）求得^[1]：

式中,$e_0$为饱和状态下前期固结应力对应的孔隙比（$e_0=N-\lambda \ln {p^{″}}_0$）,${p^{″}}_{\text{i}}$为土体在初始应力状态下对应的改进有效应力。

3.   模型验证

采用两类非饱和土试验数据进行模型参数标定和模型验证：①第一类是Sharma^[33]对膨润土-高岭土混合土开展的等向压缩和干湿循环试验数据,②第二类是Sivakumar^[34]对重塑高岭土开展的排水剪切试验数据。

3.1   膨润土-高岭土混合土

Sharma^[33]对膨润土-高岭土混合土开展了常基质吸力的等向压缩试验和常净应力下的干湿循环试验。等向压缩试验选取编号为Test 9、Test 7的试验数据进行模型验证,干湿循环试验选取编号为Test 8和Test 5的试验数据进行模型验证。试验Test 9和Test 7的基质吸力分别为200,300 kPa,其对应的净应力加载范围分别为10 kPa→100 kPa,10 kPa→175 kPa。试验Test 8和Test 5的净应力分别为10,20 kPa,其对应的基质吸力的加载范围分别为：400 kPa→20 kPa→400 kPa,400 kPa→100 kPa→400 kPa→10 kPa。试验土样的干密度${\rho }_{\text{d}}=1.24{\text{ g/cm}}^{\text{3}}$,含水率$w=25\%$。模型参数见表1,其中参数$\lambda$,$\kappa$,$N$,${p^{″}}_0$采用Gallipoli等^[10]对相同土样进行模型预测时采用的参数,参数a,b依据试验Test 11的数据拟合$\zeta$与$e\text{/}{e}_{\text{s}}$的关系进行标定,标定结果如图6所示。等向压缩试验中土-水特征曲线参数$m_1$,$m_2$,$m_3$,$m_4$依据Test 11试验中饱和度与比体积的关系进行标定得到。干湿循环试验土-水特征曲线参数$m_{1\text{w}}$,$m_2$,$m_3$,$m_4$由湿化路径下试验Test 11、Test 9、Test 7中饱和度与比体积的关系进行标定,标定结果如图7所示,其与等向压缩试验中的土-水特征曲线参数一致。$m_{1\text{d}}$,$\eta$通过标定Test 16中的干化试验及扫描试验数据确定,${\gamma }_0$,$D$依据土体的物理性质确定。$f\text{(}s\text{)}$根据图1可以得到在基质吸力为200,300 kPa时的取值分别为1.15,1.18。

表  1  膨润土-高岭土混合土模型参数

Table  1.  Model parameters of bentonite-kaolin mixture

$a$	$b$	$\lambda$	$\kappa$	$N$	${p^{″}}_0\text{/kPa}$	$m_{1\text{d}}$	$m_{1\text{w}}$
0.349	1.366	0.144	0.04	1.759	17	0.0047	0.0195
$m_1$	$m_2$	$m_3$	$m_4$	$\eta$	${\gamma }_0$	$D$	$S_{\text{res}}$
0.0195	6.911	3.929	0.03636	0.45	5	10	0.05

下载: 导出CSV 

| 显示表格

图  6  参数a,b标定

Figure  6.  Calibration of parametersa,b

下载: 全尺寸图片 幻灯片

图  7  土水特征曲线参数标定

Figure  7.  Calibration of soil-water characteristic curve parameters

下载: 全尺寸图片 幻灯片

图8为等向压缩试验的模型预测与试验数据的对比,从图中可以看出模型可较好预测常基质吸力下等向压缩试验中孔隙比及有效饱和度的变化规律,孔隙比随净应力的增加逐渐降低,有效饱和度随净应力的增加逐渐增加。在加载初期塑性模量由映射法则计算,因此初始的塑性变形相对较小,孔隙比的变化较小,当应力点达到边界面后开始产生较大塑性变形,孔隙比显著增加,有效饱和度在加载初期增长较为缓慢,在应力点达到边界面后增长较快。

图  8  模型预测与试验数据对比

Figure  8.  Comparison between model predictions and experimental data

下载: 全尺寸图片 幻灯片

图9为干湿循环试验的模型预测与试验数据的对比,在湿化过程中,随着基质吸力的减小,孔隙比逐渐增加,有效饱和度逐渐增加,土体发生膨胀,颗粒间的弯液面逐渐被水填满,黏结效应降低；在干化过程中,随着基质吸力的增加,孔隙比逐渐减小,有效饱和度降低,土体压缩,颗粒间的弯液面逐渐增加,黏结效应增强。在干湿循环过程中土体变形和饱和度的变化存在明显的滞回效应,从图中的对比关系可以看出,本文模型可以很好描述试验中孔隙比和有效饱和度变化中的滞回效应,预测结果较好。

图  9  模型预测与试验数据对比

Figure  9.  Comparison between model predictions and experimental data

下载: 全尺寸图片 幻灯片

3.2   高岭土

Sivakumar^[34]对重塑高岭土开展了不同基质吸力下的排水剪切试验,本文选取编号为8C,9C,13C,17C,18C的5组试验数据进行验证,其中试验13C,8C,18C的初始围压为150 kPa,基质吸力分别为100,200,300 kPa；试验9C,17C的初始围压为100 kPa,基质吸力分别为200,300 kPa。所有试样均在50 kPa净应力下达到吸力平衡后固结到指定围压。剪切路径均为$\Delta q/\Delta p_{\text{net}}=3$。土样干密度${\rho }_{\text{d}}=1.20$ g/cm³,含水率$w=25\%$,模型参数见表2。其中,参数$\lambda$,$\kappa$,$N$,${p^{″}}_0$,$M$,$G$结合文献[34]中饱和等向压缩试验数据与Gallipoli^[35]对相同试验数据给出的参数确定,参数a,b采用基质吸力为100 kPa的等向压缩试验数据拟合$\zeta$与$e\text{/}{e}_{\text{s}}$的关系进行标定,标定结果如图10所示。土水特征曲线参数$m_1$,$m_2$,$m_3$,$m_4$通过标定基质吸力为100 kPa等向压缩试验饱和度与比体积的关系确定。参数${\gamma }_0$,$D$依据土体的物理性质确定。$f(s)$根据图1可以得到在基质吸力为100,200,300 kPa时的取值分别为1.1,1.15,1.18。

表  2  重塑高岭土模型参数

Table  2.  Model parameters of reconstituted kaolin

$a$	$b$	$\lambda$	$\kappa$	$N$	${p^{″}}_0\text{/kPa}$	$M$	$G\text{/kPa}$
0.359	1.286	0.128	0.02	1.64	18	0.73	10000
$m_1$	$m_2$	$m_3$	$m_4$	${\gamma }_0$	$D$	$S_{\text{res}}$
0.0378	8.033	3.747	0.03635	5	1	0.05

下载: 导出CSV 

| 显示表格

图  10  参数a,b标定

Figure  10.  Calibration of parametersa,b

下载: 全尺寸图片 幻灯片

图11为模型预测结果与试验数据的对比,可以看出本文模型可以有效预测非饱和土在三轴剪切过程中的力学行为,偏应力和有效饱和度随着轴向应变的增加而增加,比体积随着轴向应变的增加而降低,最终达到临界状态。随着基质吸力的增加,偏应力逐渐增加,比体积和有效饱和度逐渐减小。

图  11  模型预测与试验数据对比

Figure  11.  Comparison between model predictions and experimental data

下载: 全尺寸图片 幻灯片

4.   结论

建立了一个同时考虑基质吸力、饱和度、孔隙结构对颗粒黏结效应影响的非饱和土水-力耦合本构模型,得到4点结论。

（1）引入可考虑孔隙结构影响的有效饱和度作为状态变量,对GGSV模型中的黏结变量$\xi$和有效应力$p^{\prime }$进行改进,采用全新的黏结变量$\zeta$和改进有效应力$p^{″}$作为模型的本构变量,通过对膨润土-高岭土混合土、重塑高岭土的试验数据进行拟合,建立了黏结变量与$e\text{/}{e}_{\text{s}}$之间的关系式,进一步得到考虑黏结效应的屈服曲线。

（2）引入修正剑桥模型椭圆屈服面作为边界面,将考虑黏结效应的屈服曲线作为$p^{″}$-$\zeta$平面的边界线,在边界面塑性理论框架内对土体的变形进行更加精准的预测。

（3）建立了考虑变形影响的水力滞后土-水特征曲线模型,实现了对非饱和土水力-力学耦合特性的描述。

（4）采用Matlab计算软件编写计算程序,利用膨润土-高岭土混合土、重塑高岭土的试验数据进行了参数标定与模型验证,结果表明,模型可以合理预测非饱和土在各向同性压缩、干湿循环以及三轴剪切条件下的水-力耦合特性。