蒸发-降雨条件下膨胀土边坡裂隙演化模拟

靳福杰，王叶娇，徐永福，杨果林，孙德安

(1.上海大学力学与工程科学学院，上海，200444；2.上海交通大学土木工程系，上海，200240；3.中南大学土木工程学院，湖南长沙，410075)

膨胀土是一种吸水膨胀、失水收缩开裂的特殊土，其胀缩性及裂隙性等特殊性质经常诱发地质灾害和工程问题，造成巨大经济损失，甚至严重危害人民生命安全[1-2]。膨胀土的工程性质对外部气候环境变化较敏感，自然边坡及工程边坡长期暴露于自然环境中，其长期稳定性易受大气蒸发-降雨作用影响。气候变化引起的干湿循环过程不仅导致膨胀土的工程性质发生衰减，而且土体不均匀胀缩会加剧裂隙发育[3-4]。殷宗泽等[5-7]发现裂隙破坏边坡土体的完整性，显著降低土体抗剪强度，并为雨水入渗提供通道，因此，裂隙发育是影响边坡稳定性的根本原因[8]。随着人类社会发展，膨胀土地区的工程活动也愈来愈多，膨胀土边坡的长期稳定性问题日益突出，研究大气作用下膨胀土边坡裂隙演化过程对解决边坡稳定性有着重要的意义。

国内外学者系统研究了膨胀土裂隙发育问题。此前，研究膨胀土裂隙演化的试验方法主要采用图像采集装置、激光扫描仪和CT机等直接观测手段[9-11]。随着计算机技术发展，数值模拟成为研究裂隙演化的重要方法。众多学者基于有限元方法模拟裂隙的发育过程。沈珠江等[12]在非饱和土固结理论基础上采用轴对称模型模拟了黏性土干湿条件下竖向裂缝的开裂与闭合过程；VO 等[13]在干燥条件下模拟了考虑水力耦合的黏土裂隙发育机理。

但裂隙的形成过程本身是非连续的，采用连续假定的有限元方法并不能较好地模拟裂隙演化过程。近年来，适用于非连续体的离散元方法已逐渐被应用到土体收缩开裂中。PERON 等[14]基于EL YOUSSOUFI 等[15]提出的经验公式，模拟细粒土条干缩开裂过程，预测裂隙发生时间以及裂隙最终形态；司马军等[16]考虑了土体抗拉强度随含水率的变化，模拟薄层黏性土干缩开裂，定量分析了表面裂隙率，其模拟结果与室内试验结果较一致；GUO 等[17]模拟了不同厚度以及边界条件下黏土开裂过程，发现边界粗糙度越大，试样厚度越小时产生的裂隙越多。

目前采用离散元方法模拟单元试验尺寸黏性土的裂隙开展成果较多，但是，在现有模拟过程中，各单元均采用相同含水率变化或人为设定含水率梯度，这与实际情况不相符，另外，因气候引起的膨胀土边坡含水率变化并导致土体的裂隙演化过程的相关研究较少。本文采用离散元方法分析蒸发-降雨作用下膨胀土边坡的裂隙发育规律与机理，首先，通过室内膨胀土干缩试验和膨胀试验标定离散元单元胀缩系数，并考虑受含水率影响的膨胀土抗拉强度变化，结合有限元非饱和土渗流分析，实现气候作用下膨胀土边坡胀缩过程中裂隙演化的离散元模拟；其次，探讨离散元模拟膨胀土边坡各部位干燥过程裂隙发育深度、形态以及降雨过程裂隙的闭合；最后，与现有试验结果[10,18]进行比较，验证用离散元模拟蒸发-降雨条件下膨胀土边坡裂隙演化的可行性。

1 室内试验方法

室内试验所用的膨胀土取自河南省南阳市宛城区红泥湾镇，取土深度约6 m，其基本物性指标如表1所示，土样属于弱膨胀土，为低液限黏土(CL)。

表1 南阳膨胀土基本物理性质Table 1 Basic physical property of expansive soil in Nanyang

土体的土水特征曲线(soil-water characteristic curve，SWCC)与饱和渗透系数是计算膨胀土边坡非饱和渗流问题的重要参数[19]。根据膨胀土原状土样的干密度与含水率制备重塑土样，抽真空饱和后放置于室内环境中，将土体含水率自然风干至目标含水率；采用滤纸法与WP4C露点水势仪，测量该膨胀土的土水特征曲线。100 kPa 以下低吸力部分采用滤纸法测量，超过100 kPa 的吸力采用WP4C露点水势仪测定，测得的土水特征曲线如图1所示。膨胀土的饱和渗透系数采用变水头法进行测定，并基于VAN GENUCHTEN[20]提出的方法估算非饱和渗透系数，如图2所示。

图1 重塑膨胀土的土水特征曲线Fig.1 SWCC of remolded expensive soil

图2 非饱和膨胀土渗透系数Fig.2 Permeability coefficient of unsaturated expansive soil

膨胀土的膨胀-收缩曲线是离散元模拟膨胀土裂隙演化过程的重要关系曲线。采用静压法制备直径为50.0 mm、高为20.0 mm 的圆饼状重塑样，控制干密度为1.6 g/cm3，含水率为20%，然后抽真空饱和。将试样放置于恒温恒湿箱内(温度为25 ℃，相对湿度为30%)进行干燥自由收缩试验，在干燥过程中，每隔1 h称试样质量并测量土体体积。测得的膨胀土收缩特征曲线如图3(a)所示。

膨胀曲线由侧限条件下断续吸水膨胀试验测定，制备直径为61.8 mm、高为10.0 mm 的圆饼状重塑样，控制干密度为1.6 g/cm3，含水率为10%。采用针筒通过多孔板孔洞定量向土样加水，并盖上盖板以减少土体水分蒸发，6 h 内百分表读数变化在0.01 mm内时认为完成该阶段膨胀，试验测定膨胀土体应变与含水率关系如图3(b)所示。

图3 南阳膨胀土收缩和膨胀曲线Fig.3 Soil shrinkage and expansion curve of expansive soil in Nanyang

2 蒸发-降雨下非饱和渗流有限元模拟

2.1 计算模型与参数

为获得蒸发-降雨条件下膨胀土边坡各部含水率分布，使用有限元软件Geostudio 中SEEP 模块进行非饱和渗流数值计算，其采用的二维渗流偏微分控制方程[21]为

式中：H为总水头；kx和ky分别为x和y方向的渗透系数；Q为边界流量；mw为土水特征曲线的斜率；γw为水的容重；t为时间。

参考文献[10]中膨胀土边坡模型尺寸，本文中采用的计算边坡尺寸如图4所示，边坡的坡比设为1∶1。边坡土的初始含水率为20%，根据测定的土水特征曲线设置边坡体内孔隙水压力为-387.9 kPa。设置边坡边界bc，cd和de为不透水边界(总流量Q=0)，ef，fa和ab依照表2中参数设置为地表大气相互作用边界。查阅南阳地区气候数据，确定模拟工况为蒸发7 d 后降雨1 d，平均气温为27 ℃，相对湿度为75%，潜在蒸发量为7 mm/d，降雨强度为10 mm/d。

图4 边坡模型尺寸Fig.4 Slope model size

2.2 计算结果

利用有限元方法模拟边坡体内的非饱和渗流过程，计算得到蒸发-降雨工况下膨胀土边坡土体体积含水率分布，如图5所示。由图5(b)可见：在模拟蒸发过程1 d后，坡体表层土体的水分最先开始蒸发，距表层1.5 cm 深度内体积含水率变化超过10%，然而，18 cm深度以下的体积含水率未受蒸发影响。由图5(c)～(d)可见：随着蒸发时间延长，大气影响深度逐渐向坡体内部扩展，但依然仅有表层部分土体(2 cm深度范围内)的体积含水率出现明显下降，说明边坡体内的体积含水率分布受大气显著影响区域范围有限。蒸发结束后，降雨1 d的边坡体内部体积含水率分布如图5(e)所示。从图5(e)可见：降雨过程导致坡体表层的体积含水率迅速增加；随着深度增加，坡体内体积含水率增加量而逐渐减小。

图5 蒸发-降雨过程边坡体积含水率分布Fig.5 Volume water ratio distribution of slope during evaporation-rainfall process

坡顶水平位置x=1.4 m、不同深度处体积含水率随时间的变化量如图6所示。由图6可见：蒸发7 d后表层体积含水率减少16.3%；随着深度增加，体积含水率变化并不明显，深度30 cm处体积含水率变化仅为2.93%。这是由于表层土体直接与大气接触，土-气水分交换频繁，体积含水率变化较大。而深层非饱和土中水分向上迁移方式为扩散作用[22]，水分传输速率较小，导致在相同时间内，随着深度增加，深层非饱和土体积含水率变化幅度逐渐减小。在蒸发过程中，坡体体积含水率变化速率逐渐减小，表层土在蒸发前2 d内体积含水率降低14.7%，而在后5 d内仅降低1.41%。这是由于坡体初始体积含水率较高，土-气界面之间的蒸汽压梯度较大，蒸发速率较快，导致表层土的体积含水率短期内快速下降。随着蒸发进行，表层体积含水率逐渐减小，导致吸力迅速增加，而对应土体表面的蒸汽压迅速下降，土体蒸发速率变小。降雨6 h 后表层体积含水率增加13.9%，随着深度增加，体积含水率增加量逐渐减小，10 cm深度以下体积含水率几乎不变。这是由于非饱和土渗透系数较低，降雨后雨水入渗至深部土体需要一定时间。随着降雨时间增加，表层体积含水率逐渐趋近于饱和体积含水率，雨水入渗深度逐渐增大，深层土体体积含水率变化具有明显的滞后效应。

图6 坡顶x=1.4 m处不同深度体积含水率随时间变化情况Fig.6 Variation of volume water ratio with time at different depths of slope top x=1.4 m

3 膨胀土裂隙演化离散元模拟

3.1 离散元模拟胀缩变化

目前，殷宗泽等[5,9,16]探讨了膨胀土的裂隙演化机理。随着膨胀土含水率降低，土体吸力增大，体积收缩，内部产生拉应力，当拉应力大于土体的抗拉强度与自重侧压力时，便会产生裂隙。而膨胀土遇水时则会发生膨胀，受到水的浸润时导致裂隙收缩甚至闭合。

在离散元方法中，将每个颗粒单元看作是土颗粒和孔隙的集合体，如图7所示，通过控制颗粒单元的体积变化模拟膨胀土的胀缩现象。二维颗粒单元胀缩基于以下公式：

图7 离散元颗粒模拟膨胀土湿胀干缩Fig.7 Moisture expansion and drying shrinkage of expansive soil by discrete element particle simulate

式中：R和R0分别为颗粒半径和颗粒初始半径；α为胀缩系数，可按式(3)或式(4)计算：

式中：εv为初始含水率到目标含水率时体积应变；e和e0分别为目标含水率和初始含水率对应的孔隙比。离散元模拟过程中颗粒的胀缩系数α采用室内膨胀收缩试验得结果进行标定。

3.2 离散元接触模型

本文使用MatDEM软件建立离散元数值模拟，单元间接触模型采用线弹性颗粒接触模型[23]，颗粒与颗粒之间靠弹簧相互接触与作用，如图8所示。颗粒间的法向力Fn和法向变形Xn可以通过颗粒间的法向弹簧模拟。

图8 线弹性模型示意图Fig.8 Schematic diagrams of linear elastic model

式中：Fn为法向力；Kn为法向刚度；Xn为法向相对位移；Xb为断裂位移。当2个相邻颗粒间的相对位移未超过断裂位移时，颗粒受到拉力或者压力的弹簧力作用；当相对位移超过断裂位移后，颗粒间拉力消失，仅存在压力作用。

颗粒间的剪切力Fs用切向弹簧来模拟。

式中：Ks为切向刚度；Xs为切向位移。在切向方向上弹簧的破坏基于摩尔-库仑准则：

式中：Fsmax为最大剪切力；Fs0为颗粒间的抗剪力；µp为颗粒间的摩擦因数。当颗粒间切向力超过最大剪切力时，切向连接断裂，此时，颗粒间仅存在滑动摩擦力，即-µpFn。

3.3 模型建模及参数设置

MatDEM 建模方式是通过模拟颗粒单元自然沉积来构建初始的堆积模型[24-25]。基于边坡尺寸以及计算机计算能力，确定颗粒平均半径为2 mm，分布系数为0.15，单元半径在1.73～2.27 mm 之间服从正态分布。对所有颗粒单元施加随机方向的初速度，通过迭代计算使颗粒运动和相互碰撞至随机位置。随后对单元施加重力作用，使单元逐渐沉积，待重力沉积结束按照边坡尺寸切割堆积体，最终坡体如图9所示，颗粒数总量为88 105个。

图9 离散元边坡模型Fig.9 DEM model for slope

随后对堆积好的颗粒设置材料参数，该模型有5 个微观参数，包括法向刚度Kn、切向刚度Ks、断裂位移Xb、初始抗剪力Fs0和摩擦因数µp，由弹性模量E、泊松比v、抗压强度Cu、抗拉强度Tu和内摩擦因数µi等5个宏观力学参数计算得到[26]。模拟中采用的宏微观参数见表2。

表2 膨胀土宏观力学性质及对应离散单元力学参数Table 2 Macro mechanical parameters of expansive soil and mechanical parameters of discrete elements

为了实现离散元模拟蒸发-降雨作用下膨胀土边坡裂隙演化过程，需考虑含水率变化引起膨胀土胀缩以及抗拉强度变化。为颗粒单元设置含水率属性，并将有限元非饱和渗流计算结果赋予每个离散元颗粒。颗粒含水率变化引起颗粒胀缩变化。膨胀土抗拉强度随水分场变化按照室内直接拉伸试验得到的土体单轴抗拉强度与质量含水率确定。

式中：Tu为单轴抗拉强度，kPa；w为含水率，%。

3.4 裂隙演化模拟结果

蒸发-降雨过程中边坡含水率分布以及裂隙演化过程的离散元分析结果如图10所示。为更加清晰地展示裂隙的演化过程，将坡顶及坡面区域裂隙处局部放大，见图11与图12。由图10可见：在蒸发1 d后，在边坡含水率变化最明显的坡顶及坡面靠近坡肩的部位率先出现裂隙(图10(a))；随着蒸发进行，边坡坡顶、坡面和坡底均出现干缩裂隙(图10(b))，并随着大气影响区域加深而不断向深处扩展(图10(c))。坡顶及坡底裂隙沿竖向发展，坡面裂隙沿垂直于坡面方向发展，这是由于蒸发过程中坡体距表层相同深度处的含水率变化幅度、土体收缩变形相同，当收缩产生的拉应力大于土体抗拉强度与自重侧应力之和时，裂隙便会出现并沿该方向继续拓展。裂隙宽度沿发育方向呈现出上宽下窄的“V”形，其主要原因是上部土体失水较多且受到自重侧应力较小，导致裂隙宽度较大；随着深度增加，土体失水较少且受到较大自重侧应力，裂隙宽度逐渐变窄。蒸发结束后，坡面部位裂隙发展深度从低到高呈现递增的趋势，这是由于底部土体受到较大自重侧应力，对裂隙的发育起到了部分抑制作用。在7 d蒸发之后对1 d降雨过程进行模拟，可以看到边坡表层土的含水率增加，吸水膨胀，裂隙收缩甚至闭合，如图11(d)与图12(d)所示。

图10 蒸发-降雨过程边坡含水率分布以及裂隙演化Fig.10 Distribution of water ratio and crack evolution of slope during evaporation-rainfall process

图11 蒸发-降雨过程中坡顶区域含水率分布与裂隙演化情况Fig.11 Distribution of water ratio and crack evolution in slope top region during evaporation-rainfall process

蒸发-降雨过程中膨胀土边坡不同位置(坡顶、坡面与坡底)处的裂隙数量以及最大裂隙深度统计结果如图13和图14所示。由图13和14 可见：蒸发2 d后，坡面处裂隙数量为11条，最大裂隙深度为12.6 cm；坡顶处裂隙数量为6 条，最大裂隙深度为9.90 cm；坡底处裂隙数量为5 条，最大裂隙深度为5.36 cm；蒸发至7 d后，坡面处裂隙仅增加5 条，最大裂隙深度为20.6 cm；坡顶处没有新裂隙产生，最大裂隙深度为10.2 cm；坡底处裂隙仅增加1 条，最大裂隙深度为5.82 cm；降雨1 d 后，坡面裂隙数量减少至6 条，最大裂隙深度降低至10.1 cm；坡顶裂隙数量减少至1 条，最大裂隙深度降低至2.83 cm；坡底裂隙数量减少至3 条，最大裂隙深度降低至1.28 cm。

图13 蒸发-降雨过程中边坡各位置处裂隙数量Fig.13 Number of cracks at different positions of slope during evaporation-rainfall process

图14 蒸发-降雨过程中边坡各位置处最大裂隙深度Fig.14 The maximum crack depth at different positions of slope during evaporation-rainfall process

由此可见，在蒸发过程中，裂隙数量以及最大裂隙深度的变化规律与含水率变化速率规律相同。在蒸发早期，边坡表层含水率减少较快，裂隙发育速率也较快，边坡各部裂隙数量快速增加；但随着蒸发时间延长，表层含水率变化速率逐渐减小，裂隙发育速率以及裂隙数量增加趋势明显下降，最大裂隙深度与裂隙数量逐渐趋于稳定，这与张家铭等[10,14]的模型试验所观察到的竖向裂隙拓展方式相似。在蒸发过程中，坡面处裂隙深度大于坡顶与坡底部位裂隙深度，这是由于坡面存在临空面导致自重侧应力较小，在相同情况下，坡面处裂隙发育深度较大。在降雨过程中，裂隙数量以及最大裂隙深度逐渐减小，这是由于边坡表层土吸水膨胀，裂隙逐渐收缩甚至闭合，导致最大裂隙深度减小，裂隙数量下降。

4 结论

1)离散元法能够较好地模拟蒸发-降雨作用下非饱和膨胀土边坡裂隙演化过程，并可定量分析裂隙数量和裂隙深度。

2)在蒸发过程中，膨胀土边坡表层含水率变化明显，随着深度增加，含水率变化量迅速减小；随着蒸发时间延长，大气影响区域逐渐向内部扩展，但扩展速率逐渐减小；降雨后，边坡表层含水率显著增加，坡内沿深度方向含水率变化存在滞后。

3)膨胀土边坡干缩裂隙最先出现于含水率变化最明显的坡肩附近，随着蒸发进行，边坡各处均出现裂隙，其开展方向均垂直于表层；裂隙数量和最大裂隙深度的变化规律均与土体含水率变化规律相同；随着蒸发进行，裂隙向坡体内部扩展，宽度逐渐收窄；坡面处裂隙的发育深度由低到高呈现出递增的趋势，且高于坡顶与坡底处裂隙发育深度；降雨后边坡裂隙收缩甚至闭合，边坡各部裂隙数量及最大裂隙深度逐渐减小。