基于离散元法的煤系土剪切带演化宏细观分析

张 榜，丰浩然，习明星，吴 灿，张 鸿

(1.南昌工程学院 土木与建筑工程学院，江西 南昌 330099；2.江西交通咨询有限公司，江西 南昌 330008)

0 引言

煤系土外观大多呈灰黑色，出露表面的土体呈浅灰色，层间黏结较弱，土体结构松散，开挖后容易风化崩解，水稳定性极差，大雨时水土流失，容易形成大的冲沟，煤系土的工程处治和利用已经成为当前亟待解决的难题。剪切带广泛存在于土质边坡的失稳破坏中，三轴压缩试验和剪切试验能明显观测到剪切带的形成，已有研究人员对剪切的形成及原理带做了大量试验和分析。从细观角度看，煤系土是非连续体，且破坏时会产生大变形，所以从细观角度分析试样破坏机理是非常有必要的。由Cundall教授提出的颗粒离散元方法[1]能较好地分析岩土体细观层面颗粒的运动变化，为从细观角度分析介质宏观力学现象提供了新思路。

毕忠伟[2]等利用PFC2D模拟双轴压缩试验，采用柔性伺服加载方式，分析了颗粒物质的宏观力学性能、剪切带的形成和发展。蒋明镜等[3-8]通过双轴压缩试验模拟了大量不同条件下剪切带的变化，比较不同黏结强度下颗粒体系的强度、变形、摩擦角、剪切带形成及位置差异，分析颗粒转动对于受力及剪切带的影响，模拟复杂土体深海能源土剪切带的形成。在实际的压缩试验或剪切试验中，土体在加载过程会产生颗粒的破碎，通过双轴压缩试验模拟颗粒破碎，发现颗粒易破碎的土体不易形成剪切带，而颗粒破碎对剪切带比较敏感主要发生在剪切带中[9-11]。研究表明力链能反映颗粒运动变化的宏观力学响应，对于模型受力性能有重要影响，稳定的力链承载能力较强且使得试样较为稳定，在力链较小的区域稳定性较差，易形成剪切带[12]。Tang等[13]模拟三轴压缩试验，分析了不同柔性墙体生成方式对试样受力的影响，发现伺服加载后模拟试验生成的剪切带，与实际三轴压缩试验生成的剪切带相符。

煤系土剪切带形成原因比较复杂，其发展过程受很多因素影响，由于物质的宏观力学现象是由组成介质颗粒的细观参数决定的，所以，要揭示煤系土的宏观变形破坏机理，就有必要从细观角度分析煤系土剪切带的形成和演化规律。目前现有对土体剪切带形成的研究文献多集中在二维双轴压缩试验模拟，存在一定的局限性，通过三维数值模拟煤系土剪切带演化的文献较少。因此，本研究采用离散单元方法，建立煤系土三轴试验三维数值计算模型，计算了煤系土在加载过程中试样的宏观力学响应，分析了颗粒旋转、体积分数及配位数的变化及力链的发展，从细观角度探讨了煤系土剪切带的形成和演化规律，为煤系土边坡工程的设计与防护提供了理论依据。

1 煤系土三轴试验数值模型建立

1.1 试验材料及颗粒级配

试验土样取自江西省宜春地区万载至宜春高速高速公路A2标段K12+110处边坡原状煤系土样，试样土体的天然物理力学指标为，土颗粒密度平均值为2.23 g/cm3，天然密度平均值为1.68 g/cm3，天然孔隙率为0.32，试件含水率W=18.2%，采用筛分法以及混合方法对土体试样进行颗粒级配测试，得到了煤系土不同粒径范围颗粒含量分布，如表1所示，土体颗粒主要集中在0.25～10 mm，占比85.41%，d10=0.129，d30=0.58，d60=1.27，不均匀系数Cu=9.84，曲率系数Cc=2.08。

颗粒级配是影响土体宏细观力学性质的重要指标，通过三维颗粒流程序(PFC3D)模拟煤系土的级配，生成的颗粒粒径分布如图1所示。从图中可以看出，数值模拟生成的级配曲线比较曲折，分布函数的曲线很光滑，通过实际试验的煤系土级配与按分布函数颗粒生成随机尺寸的颗粒、按级配模拟生成随机尺寸的颗粒进行对比，发现3条曲线总体上能较好吻合，说明数值模拟能有效生成合理的煤系土颗粒粒径，与实际的土样级配分布有很强的一致性，这表现出数值模拟现实土体的试验具有很高的可靠性，可以模拟土体在实际试验中的物理力学性质。

表1 煤系土室内颗粒级配试验结果

图1 按级配随机粒径生成的数值模拟级配曲线

1.2 三轴试验数值建模

利用PFC3D软件作为计算平台，自编程序建立三轴试验数值计算模型，试验模型墙体尺寸为高4.0 m，直径为2.0 m的三轴墙体模型。研究发现[14]，当试样较小时，边界摩擦对试验结果有很大影响，试样尺寸与颗粒平均粒径大于40时，对粗粒土的抗剪强度没有影响。因此，为了节约计算时间和提高分析效率，同时较好地分析煤系土剪切带的形成，将实际的土颗粒尺寸进行放大和集中处理。颗粒粒径在50～75 mm间等概率随机分布，生成颗粒7 558个，土颗粒密度平均值取2 230 kg/m3，试样孔隙率取值为0.32，生成模型如图2所示。

图2 三轴数值试验计算模型

根据煤系土本身的物理性质，对模型颗粒的属性进行近似赋值，颗粒的细观参数如表2所示。由于重力对试样加载试验几乎没有影响，故而未设置重力加速度。煤系土本身具有一定的黏聚力，所以该模型颗粒间设置了黏结强度。

表2 三轴试验的细观参数

2 计算结果分析

2.1 煤系土宏观力学特性分析

对墙体施加恒定围压0.8 MPa下进行伺服加载模拟计算，获取轴向压应力值，得到煤系土偏应力与轴应变曲线，如图3所示。从图中可以看出，在曲线A点以前，偏应力与轴应变曲线基本呈直线变化，试样强度呈线性增加，颗粒体系受力表现出线弹性。曲线AB段呈不规则变化，曲线斜率逐渐减小，颗粒体系强度增长逐渐减少，直到达到最大强度B处，此阶段颗粒体系表现为塑性变化；达到峰值强度之后，试样强度逐渐减小，曲线逐渐回落，偏应力下降较快，C点以后曲线逐渐趋于平缓回落。

图3 煤系土偏应力与轴应变变化曲线

随着对试样进行加载，侧向和轴向均会发生变化，体积的变化能反映试样物理力学特性，试验完成后，得到体积应变与轴应变曲线如图4所示。从曲线中可以发现，体积应变在加载初期是先下降，试样在侧向和轴向压缩下逐渐密实，A点前曲线呈直线变化，表现为线弹性，此阶段试样内部孔隙率逐渐减少，颗粒接触增加。加载到B点时，颗粒体系达到最大密实，此时体积应变εv=-1.14%到达最小值，轴应变εz=2.34%比偏应力达到峰值(2.8%)要早，AB段是非线性曲线，表现有塑性性质。B点之后，体积应变逐渐上升，颗粒体系表现软化性质，内部颗粒结构也逐渐发生改变，试样整体的密实度及承载能力逐渐下降。在初始等压加载时，试样整体密实，在达到最大密实度B点后，试样体积应变开始上升，试样出现了剪胀现象。

图4 煤系土体积应变与轴应变关系曲线

2.2 煤系土剪切带演化细观分析

煤系土的整体破坏必然从局部开始，随着局部破坏不断累积，最终达到承载的极限状态，从而引起整体破坏[15]，土颗粒作为土体的基本组成部分，了解颗粒体系局部变化特征，对于了解试样剪切带的形成及破坏机理很有必要。为了研究剪切带形成和演化特征，在剪切带内和剪切带外分别设置测量圆1和测量圆2，如图5所示。

图5 测量圆位置

2.2.1颗粒体积分数变化

如图6为煤系土体积分数与轴应变变化曲线，从图可以看出，剪切带内和剪切带外的两条体积分数变化曲线走势都是先上升之后再下降，体积分数在轴应变εz=2.8%左右达到最大值，说明在颗粒体系处于弹塑性阶段时，试样逐渐变得密实，峰值应力之后颗粒体系处于软化阶段，体积分数逐渐减小，试样逐渐变得疏松。

图6 体积分数与轴应变变化曲线

从图6可以看出，峰值应力(轴应变2.8%)之后剪切带内、外的两条体积分数曲线都呈下降趋势，在接近轴应变9.0%时出现交叉，剪切带内的体积分数开始比剪切带外的小。剪切带外的体积分数在轴应变为9.0%左右趋于平稳，但剪切带内的体积分数一直呈下降态势，而土体剪胀持续发生，这表明颗粒体系的剪胀现象与煤系土剪切带的生成和发展有很大关系。颗粒体系在初始受压时，颗粒之间发生相对移动，在未达到应力峰值前，试样整体较为稳定、强度持续增加。达到极限抗剪强度后，试样发生破坏，颗粒发生相互运动，颗粒间由于相互滚动及滑动，原本致密的结构逐渐疏松；剪切带外的颗粒在破坏之后结构重组，整体逐渐趋于平稳；剪切带内的颗粒随加载的进行，颗粒持续发生相互移动，孔隙逐渐增大，剪切带逐渐发展直至整体破坏。

2.2.2配位数变化

煤系土剪切带内外的配位数变化如图7所示。从图中可以看出，颗粒配位数在峰值应力之前呈逐渐增长的状态。试样处于峰值应力之后，剪切带内外的颗粒配位数均有减少，剪切带内的颗粒配位数减小幅度较剪切带外的减小幅度大。剪切带内配位数起伏较大，变化规律不明显。在极限破坏强度之前，颗粒体系内部变化基本一致，整体配位数增加，试样逐渐密实强度增加。达到抗剪强度之后，试样开始破坏，剪切带外颗粒接触下降不大且变化幅度较小，此区域颗粒相对比较稳定。剪切带内接触数剧烈下降，变化起伏较大，表明该区域颗粒间相互运动较为剧烈，颗粒结构不太稳定，容易破坏。

图7 配位数与轴应变关系曲线

处于峰值应力之前，颗粒间相互接触不断增加，试样的强度逐渐增加。试样处于软化阶段时，剪切带外会发生变化，但颗粒接触比较稳定；剪切带内颗粒接触变化波动较大，颗粒间发生剧烈不规则运动。这说明试样应变软化及残余强度形成与试样内部剪切带的形成密切相关。

2.2.3颗粒旋转分析

煤系土破坏从细观角度看是由于颗粒运动引起的，颗粒的移动量越大，试样的破坏也越明显。颗粒运动由滑动和转动组成，对与受约束较大的球形颗粒而言，颗粒的旋转相比于滑动更容易发生。随着外部加载，颗粒在外力作用下必然会发生运动，已有许多研究表明，颗粒旋转的分布对于试样破坏和剪切带的形成发展有着重要影响[16-17]。

颗粒的旋转分布可以通过检测欧拉角确定，如图8为煤系土试件在加载试验过程中颗粒旋转角度随轴应变分布图，颜色越深表明颗粒旋转角度越小，反之颗粒旋转角度越大。在初始加载时，颗粒开始相互运动，颗粒体系旋转分布不明显，但已有显现的趋势。在轴应变εz=1%时，试样内部受外力作用，整体运动变化不大，但局部颗粒已发生大的运动，局部剪切带逐渐发展。在εz=8%时，试样整体受力旋转，峰值应力后颗粒体系发生软化，颗粒整体运动重组，试样剪切带全部贯通。在εz=12.4%时，试样受力基本稳定，剪切带外的颗粒旋转趋于稳定，剪切带内的颗粒继续发生旋转。

图8 颗粒旋转欧拉角分布

从图8颗粒旋转分布可以看出，剪切带的形成和发展有一个过程，起初试样内部应力分布不均匀，应力较大的位置超过承载强度，造成颗粒局部运动而生产局部剪切带；随着加载的进行，试样破坏不断累积，局部剪切带不断扩散，直至最后达到临界强度后，试样失去平衡而破坏，局部剪切带贯通形成完整的剪切带(如图8虚线所示)。

2.2.4力链变化

力链能反映复杂颗粒运动的宏观力学性能，力链的粗细、位置变化、密集程度是颗粒体系在加载时受力变化直观的显现，为颗粒宏细观受力机理分析提供了有力依据。以往对于颗粒力链的分析多用二维模型分析，二维模型具有计算能力要求不高、颗粒运动及受力简单、力链易于观察等特点，在分析颗粒运动、受力变化、剪切带形成和模型破坏都取得了非常明显的效果[18]。但现实工程问题都是三维空间，很少能简化成二维平面问题，因此还需要在空间结构中进一步分析颗粒的物理力学特征。

煤系土三轴压缩试验数值模拟不同阶段力链演化如图9所示。颗粒在空间不同方向上接触受力，每个方向所受力大小均不一样，导致力链结构非常复杂。从图9中可以发现，在初始阶段，轴向和侧向受力处于等压状态，力链粗细在试样内部分布比较均匀，力链结构整体分布比较紊乱，有大量环状力链生成；随着加载的进行，侧向压力不变，轴向压力逐渐增加，力链的分布在轴向逐渐加粗加密，侧向力链还处于较细和较稀疏状态，但力链整体密度有了很大提升，试样整体也比较稳定，力链结构由环状逐渐向扁平状发展；在加载到轴应变εz=8.2%后，试样内部出现较为明显的贯通粗力链沿轴向分布，力链结构呈扁平状分布，侧向力链粗细基本没有大的变化，但力链密度不断提升。

图9 煤系土颗粒加载不同阶段力链演化

从图9中还可以看出，在荷载较大的地方出现强力链，在荷载较弱的地方出现弱力链，说明强力链沿轴向分布承受主要荷载，而弱力链沿侧向分布起着对强力链的支撑作用，让主要承载的力链不产生较大弯曲促使结构更好受力。强力链失去了弱力链支撑，颗粒结构将会不稳定；随着加载进行，试样承载能力逐渐提高，力链密度逐渐加大，力链结构经历了由环状向扁平状的转换，试样稳定性也得到了加强。总体力链分布两头较密较粗，而随着持续的加载，中间部分力链相对两头而言越发变得稀疏，说明在中间部分颗粒受到约束较少，此部分颗粒比较不稳定，更容易发生旋转，从而促使颗粒体系产生剪切带。

3 讨论

通过对煤系土三轴试验进行三维数值模拟，从细观角度分析剪切带的形成及发展，发现初始试样内部由于应力分布不均匀，颗粒局部开始产生剪切带；随着加载进行，局部剪切带不断发展，当达到抗剪强度后，局部剪切带贯通形成完整剪切带。结合试样宏观力学响应，发现煤系土软化现象、剪胀现象与剪切带的形成有直接关系。在实际岩土工程问题中，将工程问题简化为平面二维问题与实际存在差距，与毕忠伟[2]建立的二维模型相比，本研究建立三维空间模型相比二维平面模型而言更能反映实际复杂的颗粒受力，而且剪切带的形成和发展规律与二维模型高度一致。结合颗粒配位数、体积分数、颗粒旋转及力链的演变，从细观层面分析了试样剪切带的形成，使我们对煤系土的剪切破坏有了更深的认识，为处理复杂的实际工程问题提供了理论支撑。

本研究采用刚性墙体进行加载，无法观察到类似于实际三轴试验中的剪切破坏带，只能通过观察颗粒的旋转及力链演变来分析剪切带的形成和发展，这还有待进一步深入分析剪切带的破坏机理，今后可以从改刚性墙体为柔性墙体，改变不同颗粒形状、级配、接触模式及不同细观参数等方向进行综合分析。

4 结论

本研究采用离散单元方法，建立了煤系土三轴压缩试验数值模拟计算模型，对煤系土剪切带形成和发展进行仿真分析，计算了颗粒体系的宏观力学响应，从细观角度研究了煤系土在三轴试验数值模拟过程中颗粒的体积分数、配位数、颗粒旋转及力链的变化规律，得到如下结论：

(1)煤系土体剪切带的形成是剪胀现象产生的重要因素。试件峰值应力之后，剪切带内颗粒的体积分数一直处于下降状态且孔隙率大于剪切带外，这说明剪切带贡献很大部分体积的剪胀。

(2)颗粒旋转是煤系土剪切带形成的主要原因。通过颗粒相互旋转，剪切带内颗粒配位数降低，颗粒接触减少、孔隙率增加，使得试样强度下降，颗粒结构变得不稳定。此外，颗粒的几何形状将直接影响颗粒的旋转和颗粒之间的摩擦耗能，下一步将深入分析颗粒几何形状对土体宏观抗剪强度的影响。

(3)颗粒力链密度随着加载进行逐渐加大，力链结构经历了由环状向扁平状的转换，试样稳定性得到加强，在加载后期试件颗粒的力链分布为两头较密较粗，中间部分力链变得稀疏，颗粒受到约束较少，更容易发生旋转，从而促使颗粒体系产生剪切带。因此，力链能反映复杂颗粒运动的宏观力学性能，力链的粗细、位置变化、密集程度是颗粒体系在加载时受力变化直观的显现，为颗粒宏细观受力机理分析提供了有力依据。