非饱和土粒间毛细作用的微观不连续变形分析*

2021-07-19李同录乔志甜张常亮郭龙骁

工程地质学报 2021年3期

李强李同录乔志甜张常亮李萍沈伟郭龙骁

(①长安大学地质工程与测绘学院，西安 710054，中国)(②黄土高原水循环与地质环境教育部野外科学观测研究站，正宁 745399，中国)(③博洛尼亚大学地质和环境科学学院，博洛尼亚 40121，意大利)(④九州大学土木与结构工程学院，福冈 819-0395，日本)

0 引言

非饱和土是一个固-液-气三相体系，土-水相互作用是决定其力学性质的内在因素(Fredlund et al.，1993；陈鸿宾等， 2019; 高英等， 2019)。而土-水相互作用是一个复杂的物理化学过程，包括结合水的吸附作用以及毛细水的表面张力作用。其中：毛细作用是基质吸力的重要来源之一，相对饱和土，基质吸力增加了非饱和土颗粒间的有效应力，从而提高了土的抗剪强度，对非饱和土性质有重要影响(邢鲜丽等， 2015；蒋秀姿等， 2018；李威等， 2018；李同录等， 2019)。

由于毛细作用是一种微观现象，其本质需要从微观角度去探索。早期对非饱和土中毛细作用的微观研究基于理想的接触球模型(Fisher， 1926)。这种方法将土颗粒简化为等大或不等大的光滑球体，将其抽象为立方体或四面体的结构形式的土粒骨架，在土粒骨架间施加毛细水作用，分析土中的毛细水分布、表面张力、球间引力等。这种简化模型有助于定性地理解毛细水在非饱和土中的作用。随着数值计算能力的提高，为利用数值模拟方法来分析土粒间毛细作用提供了可能。土体在微观上是不连续的，因此目前主要采用离散元法(DEM)和非连续变形分析法(DDA)等离散介质模型来模拟其微观力学行为。根据颗粒形状，离散元(DEM)又可分为块体离散元和球体离散元2类。目前基于块体离散元对非饱和土的微观力学性质研究较少，大多基于球体离散元。对于球体离散元，颗粒的建模方式一般分为两种，第1种是使用单一的圆盘或球形，另一种是多个圆盘通过黏结力进行黏结形成类多边形(周伟等， 2009；李杨等， 2017)。这种方法能够有效地模拟颗粒的破裂过程，但形成的颗粒与实际颗粒的体积有一定差距，会导致计算时产生误差。在利用球体离散元模拟非饱和土中的毛细作用时，通常采用圆盘或球形颗粒构建土骨架，在相同大小颗粒接触处施加相同大小的应力以考虑毛细力的作用(Jiang et al.，2004； Shen et al.，2016)。这种算法的主要问题在于未考虑不同含水率下毛细水的分布状态，而毛细水的分布状态与土体力学性质有直接关系。DDA最初用于模拟块状结构岩体的变形破坏过程(邬爱清等， 1997；周少怀等， 2000；焦玉勇等， 2007)，由于土体中颗粒多呈多棱角的不规则形状，所以近年来开始用于分析土体的微观力学作用。其中张国新等(2000)利用规则的多边形构建微观砂土DDA模型，研究了砂土微观变形机理。郭培玺等(2008)通过将假想多边形随机填充的方式构建了粗粒料DDA模型，基于此模型研究了粗粒料的力学特性及颗粒运动规律。郭龙骁等(2017)将实际黄土颗粒的形状通过人工摆放的方式构建黄土微观模型，模拟了微观黄土单向固结实验。这些研究虽然对微观模型做了一些改进，并表明了DDA用于土粒微观力学作用研究的可行性。但都基于传统的二维DDA算法，尚未考虑毛细水作用。实际问题大多是三维，所以近几年开始发展三维DDA算法(Zhang et al.，2016； Meng et al.，2019)，但由于其复杂的接触关系造成自由度的大幅增加，使得三维DDA算法进展缓慢，目前还很不成熟。因此，本课题组近年来开始研究考虑非饱和土毛细水的二维DDA算法(Guo et al.，2019)。

本文在二维DDA算法的基础上，提出了一种用于非饱和微观结构力学分析的拓展DDA算法。该算法可算出不同饱和度时水分在土颗粒之间的分布形态及弯液面半径，再用Young-Laplace方程计算出粒间毛细负压，用毛细负压和表面张力表征毛细水对土颗粒的作用，这为进一步开展非饱和土各种力学行为的分析奠定了基础。基于该算法，构建了一个中粗颗粒微观模型，对不同含水率下土体中的毛细水分布及基质吸力进行计算，获得模型土的土水特征曲线(SWCC)。

1 非饱和土粒间毛细力的计算

当土处于非饱和状态时，土粒之间主要由水联结，在土中形成基质吸力。当土特别干时，如黄土含水率一般在7%以下，粒间联结以结合水为主(王铁行等， 2014)；当含水率增高时，结合水膜增厚，结合水联结减弱，土粒间以毛细水作用为主。天然状态下伏于地下的土都偏湿，以毛细水作用为主，因此暂先不考虑结合水的作用，只考虑土粒间的毛细水联结。

1.1 基本假设与算法流程

通常所说的毛细现象就是指液体在毛细管中的移动，表面张力和弯液面内外压差是造成这一现象的本质因素(高世桥等， 2010)。其中表面张力是水体在毛细管中移动的起因，内外压差是其动力来源。

在岩土工程领域，通常将土壤孔隙理想化为毛细管，并将弯液面的内外压差称为基质吸力，其与交界面几何形状的关系满足Young-Laplace方程，在二维情况下表现形式为：

(1)

式中：ua为孔隙气压；uw为孔隙水压；ua-uw为基质吸力；Ts为表面张力；θ为土水间的接触角；r液面曲率半径。

式(1)中，表面张力Ts与接触角θ在温度和材料等条件确定的情况下均为常量，所以基质吸力的大小主要取决于弯液面的曲率半径r。由于实际土体中毛细作用复杂，为了建立模型，对复杂条件需进行适当简化，本文的模型基于如下假设：

(1)土颗粒表面为理想表面，忽略增湿和减湿接触角的不同；不考虑矿物表面特征对接触角的影响，假定所有土颗粒的接触角都相等。

(2)不考虑渗流时间和渗流路径，认为土体系统内任意时刻毛细水分布均处于平衡状态。平衡状态的条件是：水压、气压处处相等；水和气在液面处凝结和蒸发平衡；所有弯液面曲率半径相等。

(3)土体中所有土颗粒为凸形颗粒。

基于上述假设，用Young-Laplace方程可计算出粒间毛细水的基质吸力。本算法是在原有DDA算法(Shi， 1992)的基础上进行拓展，拓展部分的算法如图1所示。

图1 拓展DDA算法中粒间毛细力计算流程

1.2 颗粒间毛细水量的计算

土粒间的毛细水分布在颗粒接触处和距离最近处。传统的DDA算法可以识别所有土粒间的接触点和接近点，对于任意接触点，其邻近有两个或多个土粒，首先判定两个土粒间毛细水的分布和面积，在此基础上可判定多个土粒间毛细水的分布和面积。

两个土粒间毛细水的分布是由曲率半径相等的两个弯液面和两个土粒的边界所限定。由于土粒位置固定，若给定弯液面的曲率半径r、确定出两个弯液面的圆心点，就可算出毛细水的分布区域。确定弯液面圆心可采用圆心轨迹交汇法(图2)，具体计算过程如下：

图2 双颗粒系统中毛细水分布示意图

(1)以选定的某一颗粒A为对象，假定颗粒A表面任意一点为液面的起始点SA1，过SA1作颗粒边的外法线，将外法线逆时针旋转接触角θ，在旋转后的线上确定与SA1的距离为r的一点OA1(图2a)。

(2)使线段SA1OA1长度与外法线夹角大小不变，沿颗粒A表面逆时针移动一周，端点OA1在颗粒A外围形成一条平行于颗粒表面直边的轨迹线。

(3)该轨迹线在角点处存在缺口，故颗粒角点处的轨迹需要单独处理。以角点为圆心，以r为半径，作一段圆弧，将缺口处的轨迹线连接起来(图2c)。按此方法处理后，颗粒A外围将形成一条完整的圆心轨迹线。

(4)当相邻两个颗粒之间的最近距离小于2r时，才可能形成毛细水弯液面。故在距颗粒A表面距离为2r的范围内寻找相邻颗粒B。在颗粒B上随机选取一点SB1，过SB1做颗粒表面线的外法线，将外法线顺时针旋转接触角θ，并确定与SB1距离为r的一点OB1，按如上相同的做法使SB1OB1顺时针绕着颗粒B表面一周，做出相邻颗粒B的圆心轨迹线。

(5)颗粒A和颗粒B的圆心轨迹线存在两个交点，分别以这两个交点为圆心，以A颗粒的边为起点，B颗粒的边为终点，逆时针画出两个圆弧。这两个当中有一个是凹液面，另一个是凸液面。凸液面不合理，应去除。

(6)做两个液面与A颗粒或B颗粒交点处的液面切线，设该切线在液体一方与固液交界线之间的夹角为β。当βi>90°时，液面为凸液面，将其删除；当βi<90°时，为凹液面，将其保留。

(7)将颗粒A和颗粒B的角色置换，重复步骤(1)～(6)则可确定另一侧的凹液面(图2b)。将两侧凹液面组合，则确定了毛细水分布区(图2d)。

图3 多颗粒中毛细水弯液面形式

在图3的3个土颗粒中，共有3个毛细区， 6个液面。当土中含水量较低时，各液面相互之间不会产生重叠(图3a)。当土中含水量逐渐增多时，弯液面C2、C3、C5会出现重叠区域，在未重叠区域内会出现一个封闭的气泡(图3b)。在实际土体中，这个气泡在二维平面上应当呈圆形，并非曲边三角形，但这并不影响其中含水量的计算以及对土性质的分析。

可以看出，弯液面曲率半径r的大小与土中的含水率Vw之间存在相关关系。若人为给定含水率Vw，需通过反复迭代计算出r的值，具体算法如下：

(1)先假定一个弯液面初始半径区间[ra，rb]，按上面所述的算法分别确定ra和rb所对应的所有土粒之间的弯液面，弯液面和所连接的土粒边界构成封闭的毛细水分布区。随即将毛细水分布区的像素指定为特定且与模型中其他任何要素不同的颜色属性。最后通过统计整个模型中毛细水像素所代表的面积，就可以计算毛细水的含量Va和Vb，这种方法优点在于计算简单，精度相对较高，并且不会重复计算毛细水的重叠部分。

(2)若计算出的Vb>Vw>Va，采用黄金分割算法逼近Vw。若Va>Vw或Vb

1.3 毛细作用力计算

在实际土体中，毛细水作用在土颗粒间的力包括凹透镜形水体内的毛细负压(ua-uw)作用在土颗粒上的力和两侧的表面张力Ts。其中毛细负压(ua-uw)是作用在土颗粒与水的接触面上的面力，表面张力Ts是作用在三相界线上的线力。据此，该算法在二维平面上进行了相应简化，将毛细负压作用在被毛细水浸湿的颗粒边界上，方向指向浸湿边界的外法线方向(图3)。计算时，将毛细负压(ua-uw)作用在土颗粒上的力简化为(ua-uw)×L的力，记为fc，其中，L为毛细水浸湿段长度，毛细负压(ua-uw)在r已确定的情况下，可由式(1)获得，作用点位于浸湿段中点，方向指向颗粒表面的外法线(图4)。

将表面张力Ts仅作用于固、液、气三相交界点处，在温度等条件确定的情况下，其值大小为一常量，方向沿弯液面的切线由两端指向中间(图4)。

图4 单颗粒上毛细作用力计算

在直角坐标系中，将所有力都分解到x轴和y轴方向上，再求和。同时还必须计算所有毛细力产生的力矩，以颗粒的形心为转动中心，计算其合力矩(图4)。将合力及合力矩添加到DDA的荷载向量中，对每个颗粒都做如此计算，整个土粒间的毛细作用就施加在计算模型中了。值得注意的是，封闭气泡中的吸力和表面张力未直接作用于土颗粒上，故不考虑气泡周围的力。

2 黄土微观结构模型的构建

现以黑方台L1黄土样为参考，建立其理想土骨架结构DDA模型，来研究黄土中的土-水相互作用。该土样的干密度为 1.32 g·cm-3，比重为2.71，孔隙比为1.06，液、塑限分别为27.5%和18.3%，图5为其粒度分布曲线。

图5 样品粒度分布曲线

考虑到计算效率，模型的颗粒数量不宜太多。根据图5， 20～30μm是颗粒组分相对较多的部分，本文只选取该范围内的粉粒建立理想黄土微观骨架模型。在放大100倍下显微镜下观察了大量的黄土中该范围内粉粒散粒的形状，图6a为其中部分颗粒，发现大部分粗颗粒多呈不规则多边形或扁平状。将镜下观察的颗粒进行矢量化，由于矢量化后的形状包含的角点数量较多，为保证计算效率，需简化。简化原则为当多个角点位于同一边或近似在同一边时，仅保留两个角点(乔志甜等， 2020)。简化后的多边形的边数尽量不大于8(图6b)。

图6 镜下粗颗粒形状及多边形简化后形状

从中随机选取350个大小为20～30μm不同形状的颗粒用于生成土的模型。然后将选取的颗粒随机丢落于内径为0.54mm×0.40mm的样盒模型中，土中颗粒会在重力的作用下达到相对稳定状态，形成的模型作为初始模型(图7)，相应的初始孔隙比e为0.39，小于实际黄土的孔隙比。在模型生成过程中未考虑空气阻力和静电引力的影响，还难以形成与实际土一致的孔隙结构，对此还有待进一步研究。

图7 微观黄土骨架DDA模型

由于黄土中20～30μm的粉粒大多为石英，密度可参考石英矿物的密度。此外黄土颗粒最初来自物源区的岩石，模型使用的弹性模量E可以参考岩石中(花岗岩)的相应参数。此外，本文只考虑毛细作用，故将黄土颗粒之间的黏聚力设为0。土样盒也为一组DDA单元，相对于黄土颗粒，样盒为刚体，故将样品盒的弹性模量设为黄土颗粒弹模的20倍。样盒与土粒之间的内摩擦角、黏聚力均设置为0。具体参数如表1。

表1 模型参数表

3 毛细作用下土水特征曲线模拟结果

利用本文拓展的DDA算法，给定一系列的含水率，可计算出该土样的饱和度及相应的毛细水弯液面半径，用Young-Laplace方程(1)可计算出不同半径相应的基质吸力，计算结果如表2所示，由此得到饱和度与基质吸力之间的关系，即SWCC曲线(图8)。图9给出了其中6种饱和度工况下的毛细水分布。

表2 不同饱和度下弯液面半径及吸力值

图8 SWCC曲线模拟结果与实测结果

图8中的实测数据来自黑方台L1原状黄土压力板仪实验结果。由图8中的DDA模拟结果可以看出，SWCC曲线存在着明显的边界效应区(0～8kPa)、过渡区(8～110kPa)、残余区(>110kPa)。在边界效应区，饱和度变化基本稳定，基质吸力在该范围内随饱和度的下降缓慢增加。该阶段土中水相分布状态如图9a，土体内部大部分孔隙均被水体充满，只有少量气泡存在于水中。当饱和度进一步下降，SWCC将处于过渡区，边界效应区和过渡区之间的界限值即进气值，模拟试验获得进气值约为8kPa，对应的饱和度为90.3%。在进气值附近，土体中的水相分布状态如图9b，此时土体中的气泡开始扩大并出现少量独立的弯液面，土体开始进入非饱和状态。图9c、图9d是完全处于过渡阶段时土体中水相分布状态，土体中气体开始占据大部分孔隙，独立的弯液面数量明显增多。此时，饱和度下降速度较快，基质吸力的增幅相对较大。随饱和度进一步降低，SWCC将处于残余区。过渡区与残余区同样存在一个界限吸力值，即残余值，该值大小约为110kPa，对应的饱和度6.31%。在该饱和度附近土体中水的分布状态如图9e，此时土体中的水仅分布于颗粒的接触位置及较小的孔隙中。图9f是完全处于残余阶段时土体中水的分布状态，土体中大孔隙均被气体占据，极少量的水分布在颗粒的接触位置，土体中水量较少，基质吸力较高。并且在该阶段，基质吸力随饱和度缓慢下降而迅速增加。

图9 不同饱和度下毛细水的分布状态

综上所述，模拟结果符合典型土水特征曲线(McQueen et al.，1974)的形态。饱和度相对较高的情况下，模拟结果和实测曲线较为吻合。原因在于：此时，毛细作用是基质吸力的主要来源。毛细水DDA算法以及模型虽对复杂条件进行了一些简化，但并未脱离实际情况。模型中使用的颗粒粒径分布以及其形状，都是参照实际黄土构建的。并且，构建的理想模型基本包含了实际土体中主要的孔隙结构。所以，在毛细作用起主要作用的区间内，能够较为准确地计算出基质吸力。但随着饱和度下降，与实测曲线有所偏差。这一差距主要是因为：实测的土水特征曲线中的吸力应包含毛细水的毛细作用和结合水的吸附作用两部分(Tuller et al.，2005)，一般在残余含水率之前，毛细作用起控制作用，之后吸附作用起控制作用。此外，真实黄土中存在一定量的黏粒，黏粒往往会吸附较多的结合水，这类吸附作用主要来源于范德华力、电场力和水合作用，而非毛细作用。并且在非饱和土中，接触角并非是个常数，会随外界条件的变化而变化。一般在饱和度降低的时候，接触角会变小(Yang et al.，2012；杨松等， 2015)，进而造成基质吸力变大。这些影响因素中，未考虑结合水的吸附作用是导致低含水率模拟结果与实际结果误差的主要原因，所以建立合理的结合水非连续变形算法是作者们将要进一步研究的问题。

总的来说，本文使用的模型仅考虑了20～30μm的粗粒，且未考虑吸附作用和接触角的变化，因此，随着饱和度下降，差距也随之明显。但从整体上来看，该算法能够较为合理地考虑毛细作用，并能计算不同饱和度条件下的毛细水分布状态。