粗粒土应力腐蚀蠕变影响因素的离散元分析

王 瑶，魏进兵，温仁节，王 琛

(1.四川大学 水利水电学院， 四川 成都 610000；2.四川蜀禹水利水电工程设计有限公司， 四川 成都 610000)

粗粒土是一种工程中广泛应用的土体材料，大量的国内外试验结果和监测资料[1-4]均表明其具有明显的长期变形特性，粗粒土的蠕变性质研究一直是岩土工程领域的研究热点之一[5]。

土的蠕变性质研究一般包括蠕变机理、蠕变性质与规律、蠕变模型及其应用等方面，其中，蠕变机理研究有助于对蠕变性质与规律的认识，亦有助于建立更符合土体变形规律的蠕变模型。迄今，一般认为粗粒土的蠕变机理包括[6-10]：(1) 土颗粒棱角破碎；(2) 颗粒内部的微裂缝受应力作用随时间发展为颗粒劈裂；(3) 土颗粒材料本身的流变性质。前两种蠕变机理都认为粗粒土蠕变的主要细观机理为颗粒破碎，颗粒位置发生重分布进而导致了蠕变现象。

考虑颗粒破碎的粗粒土蠕变机理研究可基于宏观试验和细观数值试验展开。宏观试验现今主要为土样级配分析，其方法简便且可为细观数值试验提供细观参数标定标准；缺点是仅仅能获取某些特征点的颗粒破碎特性，无法获取整个蠕变过程中连续的颗粒破碎结果。细观试验(如图像技术)和数值试验(如离散单元法)可实现连续监测，从而获取试样蠕变全过程相关参数。离散单元法[11]等数值试验方法相对于图像技术等细观试验方法来说，经济迅捷，并可监测蠕变过程中的颗粒运动、组构和力链演变、破碎形式及破碎分布等。Tran等[12]在PFC2D中将堆石颗粒视为由刚性颗粒胶结成的团粒，并让粒间胶结键强度随时间增长而降低，成功模拟了粗粒土在荷载作用下随时间增长而产生裂缝直至发生颗粒破碎的现象，并由此来研究堆石坝的蠕变。Potyondy[13]基于亚临界裂缝扩展理论[14]认为胶结键半径受应力腐蚀而逐渐缩减，在PFC 4.0中编写了相应的应力腐蚀程序，成功模拟了岩石试样的蠕变问题。在PFC 5.0版本中，应力腐蚀程序的实现及其模拟粗粒土蠕变的适宜性尚需进一步研究，为进一步模拟粗粒土蠕变细观机理，也需探讨参数对蠕应变的影响并确定适合的细观模型参数。

本文基于应力腐蚀理论在PFC 5.0中编制了考虑颗粒破碎的粗粒土蠕变离散元程序，分析接触模型参数中颗粒摩擦系数、胶结键粘结强度与颗粒接触刚度及应力腐蚀模型参数中材料参数、激活应力对蠕变的影响规律及其敏感性；基于室内试验结果，验证应力腐蚀模型的应用，标定适用于红砂岩粗粒土的细观数值试验参数，以期为粗粒土蠕变细观机理的进一步研究奠定基础。

1 粗粒土应力腐蚀蠕变的数值模型

1.1 应力腐蚀理论与程序

应力腐蚀原理是由岩石断裂力学中的亚临界裂缝拓展理论类比而来，将土体视为一种含有微裂缝且不连续的材料，认为应力腐蚀是长期应力作用下导致土体裂缝发展的主要原因[14]。应力腐蚀机制中存在两个重要的界限值，即发生应力腐蚀与裂缝拓展突变的界限值，分别称为起裂韧度和断裂韧度。当应力达到起裂韧度时，微裂缝开始缓慢拓展；当应力接近或超过断裂韧度时，裂缝发展速度会产生突变，迅速发生颗粒破碎，宏观上表现为蠕变变形显著增长；当大量的颗粒破碎后，新的颗粒破碎发展速率趋于缓慢，宏观上表现为蠕变趋于稳定。

基于这一理论，可在PFC中通过胶结键拉应力的函数来模拟胶结键的胶结半径随时间缩减的过程，从而模拟岩土材料的裂缝随时间发展的过程[13]。胶结键直径缩减速度公式：

(1)

式中:dc为胶结键直径；σa为激活应力；σt为胶结键抗拉强度；σ为胶结键上的最大拉应力；β1和β2均为材料参数。作用于胶结键上的最大拉应力是判断是否发生应力腐蚀的基准，当σ≥σa时，才会开始发生应力腐蚀，当σ≥σt时，胶结键则会破坏进而发生颗粒破碎。

在进行蠕变的离散元分析时，蠕变时间一般有两种考虑方式：其一是用黏弹性接触模型(如：Burgers模型)模拟真实的物理时间。然而因离散元软件的物理时间步通常小于百万分之一秒，对几天乃至几年的蠕变模拟需借助云计算或超算才能实现；其二是通过考虑蠕变机理(如应力腐蚀模型)，用蠕变机理模型中的时间变量替代真实的物理时间来模拟蠕变，可大幅度节省计算时间，更具有可行性[13]，本文采用该方法进行分析。

在PFC 5.0中编写的应力腐蚀程序流程如下：(1) 对试样施加伺服偏应力，判定试样是否破坏，若试样破坏，则结束蠕变模拟，若试样未破坏，则执行应力腐蚀程序；(2) 判定是否达到设置的蠕变模拟时间，若当前蠕变时间已达到设置的模拟时间，则结束模拟，若未达到设置的模拟时长，则执行应力腐蚀函数；(3) 遍历每一个胶结键，判定是否发生应力腐蚀，若不发生应力腐蚀，则结束蠕变模拟，否则计算应力腐蚀速度，缩减胶结键半径，同时更新蠕变时间；(4) 再次施加应力边界条件，按上述相同步骤反复迭代计算，直至数值模拟满足设置的蠕变模拟时长、试样不发生应力腐蚀和试样破坏三个判定条件之一，即结束蠕变模拟。

1.2 试样建模

1.2.1 墙体生成

PFC中的墙体一般为刚性墙且无质量，墙体与墙体之间不存在相互作用，故不能运用牛顿第二定律或对墙体直接施加力或弯矩作用，只能通过对墙体赋予速度的方式给颗粒体系施加伺服力作用，墙体与颗粒间的相互作用则通过改变墙体刚度的方式进行调整。数值建模参考室内试验仪器尺寸，宽100 mm、高200 mm。如图1所示。试样包含四面墙体，上下两面墙分别模拟试样帽和试样底座，用于施加伺服偏应力；左右两面墙为柔性伺服边界，模拟橡皮膜对试样施加伺服围压。

图1 试样墙体示意图

1.2.2 颗粒生成

颗粒生成包括级配设计和Cluster团粒生成。

在颗粒级配设计中，为更好地分析蠕变过程中颗粒破碎与分维值的关系，同时避免孔隙结构的影响，根据分形理论[15]将土料设计为具有自相似分形结构的级配，有:

(2)

式中:M(d

图2 级配曲线

在颗粒生成时，每个粒组对应使用一种模板生成颗粒，为使颗粒破碎后体积守恒，模板颗粒设置为两两相切。忽略各粒组土颗粒形状的随机性，统一用一种模板颗粒代替。其中，20 mm～10 mm模板代表扁长颗粒，10 mm～5 mm模板代表尖锐棱角颗粒，5 mm～2 mm模板代表浑圆度较高颗粒，模板颗粒如图3所示。

图3 团粒模板

根据面积等效原理可得每个模板的等效半径与颗粒数目为：

(3)

式中:R为大颗粒的面积等效半径，n组成模板的小颗粒个数，ri第i个小颗粒的半径，通过半径最终计算等效粒径与各模板团粒个数见表1。

表1 团粒粒径和数目

建模时拟使用的Cluster团粒因其能破碎的特性可有效用于模拟颗粒破碎，但在PFC中无法直接生成。为解决这一问题先根据级配利用不可破碎的Clump模板生成试样，再用自编写的函数转换成可破碎的Cluster团粒，如图4所示。

生成试样后，由转换命令给每个接触两侧球体赋予额外变量，识别不同的接触，当接触两侧球体的额外变量相同时，该接触定义为颗粒内部接触，否则定义为Cluster之间的接触。Cluster内部接触采用平行粘结模型，模拟试样颗粒的胶结作用；Cluster之间采用线性接触模型，模拟粗粒土之间的接触关系。

图4 数值试样模型

1.3 基准参数

表2 模型参数值

2 粗粒土应力腐蚀蠕变影响因素分析

在粗粒土蠕变数值试验分析中，参数敏感性分析是极为重要的一环，通过分析不同参数对蠕变的敏感性，可定性确定对蠕变起主要影响的参数，排除低相关性参数，为今后进一步的蠕变机理研究提供一定的参考[19-21]。本文选取蠕应变作为目标函数值，在基准参数的基础上，通过改变单一参数，控制其它因素不变来分析变化参数的敏感性。由于各参数的量纲并不统一，且取值相差较大，故引入无量纲参数敏感度函数来度量参数敏感性。针对蠕变问题，参数敏感度函数可定义为蠕应变的相对变化值δε与影响参数相对变化值δx之比为：

(4)

(5)

(6)

式中:Δε为蠕应变的绝对变化值；Δx为参数的绝对变化值；S(x)为参数敏感度函数；ε为蠕变基准值；x为敏感参数基准值。

本文进行围压σ3=300 kPa，应力水平S=0.75条件下的参数敏感性分析，分析各参数对最终蠕应变εf的影响，计算敏感度函数值，其值越大参数越敏感。

2.1 接触模型参数对蠕变的影响

最终蠕应变εf-摩擦系数f关系如图5(a)所示。随着摩擦系数的增大，颗粒之间的相对运动愈加困难，所消耗的能量亦越多，蠕应变越小；但这种影响并不是无止境的，因影响粗粒土变形的主要因素还是颗粒破碎，摩擦系数仅能在很小范围内影响到蠕应变，所以蠕应变减小速率逐渐变缓。颗粒接触摩擦系数可定性的反映母岩颗粒的表面摩擦特征，颗粒表面愈粗糙的粗粒土蠕变越小。

图5(b)给出了最终蠕应变εf-胶结键抗拉强度σt关系，εf随σt的增大而减小。原因为胶结键抗拉强度的增大使得胶结键愈发难以破坏，同时也降低了颗粒应力腐蚀速度，这两种影响都会导致颗粒破碎量减少，进而使蠕应变降低。表明粗粒土的蠕应变随母岩抗拉强度的增大而减小。

颗粒刚度kn与颗粒密度ρ对最终蠕应变的影响较小，如图5(c)和图5(d)所示。颗粒刚度决定了颗粒的接触位移大小，故在一定程度上决定了颗粒间的相对位置，但颗粒刚度决定的弹性应变在伺服阶段已基本完成，而基于应力腐蚀的蠕变变形主要来自于颗粒裂缝发展导致的颗粒破碎，故颗粒刚度对蠕应变的影响很小；颗粒密度变化仅影响颗粒加速度，而伺服调节始终会使颗粒达到平衡位置，故颗粒密度对粗粒土蠕变的影响不大。

综上，粗粒土的母岩性质对其蠕变有决定性影响，粗粒土颗粒的表面摩擦系数越大蠕应变越小，抗拉强度越大蠕应变越小。

图5 接触腐蚀模型参数-蠕应变关系曲线

2.2 应力腐蚀模型参数对蠕变的影响

应力腐蚀模型参数包括激活应力σa和材料参数β1、β2。蠕应变随激活应力σa的增大而呈线性减小，如图6(a)所示。由粗粒土母岩性质决定的激活应力σa作为应力腐蚀的阈值条件，直接影响到参与应力腐蚀的胶结键数目，随着激活应力增大，颗粒破碎急剧减少，蠕应变随之减小。本文激活应力为0.3 MPa～0.9 MPa，该线性关系在此范围外或不成立。

蠕应变随着β1和β2的增大而增大，如图6(b)和图6(c)所示。由式(1)可知胶结键应力腐蚀速度与β1成正比，β2是应力腐蚀速度函数的指数参数，故当β1、β2增大时，应力腐蚀速度变快，使得相同时间内的颗粒破碎量显著增大，蠕应变随之增大。随着蠕应变的增大，粗粒土孔隙逐渐变小，颗粒位置调整受限，蠕应变增长放缓。

各参数的参数敏感度如表3所示，四个接触模型参数中f、ρ、kn的敏感度函数值均不大于0.05，仅σt大于0.2，总体表现为对蠕应变影响较不敏感；而应力腐蚀模型参数的敏感度函数值均大于0.2，尤其β2达到了0.3，综合可知应力腐蚀模型参数比接触模型参数更加敏感。

2.3 参数标定

在离散元数值试验中，为了确保计算模型最终能实际模拟土体的宏观物理力学性质，标定出一套适用于室内试验结果的细观参数至关重要，可为进一步分析粗粒土的蠕变机理奠定基础。

图6 应力腐蚀模型参数-蠕应变关系曲线

表3 参数敏感度

室内试验选取粗粒土为人工破碎弱风化红砂岩，比重由虹吸桶法测得Gs=2.68,制样干密度通过尝试法即在承膜桶内将分维值为1.5的试验土料击实至人工所能达到的最密实状态获得ρd=1.70 g/cm3，孔隙比e=0.57这一干密度较工程中偏小，其原因为试样剔除了1 mm以下土料，且分维值D=1.5偏小，致试样具有较大的孔隙结构。试样尺寸与级配与数模试样保持一致，采用SR-20型三轴蠕变仪进行三轴饱和固结排水试验，采用单级加载方式，不同围压下的破坏偏应力通过砝码间隔50 kPa加载，直至轴向应变达到10%或试样发生持续不稳定变形为止。本文以试验围压σ3=300 kPa，应力水平S=0.25、0.50、0.75和0.90为例进行分析。由图7室内试验蠕变结果可知：红砂岩无黏性粗粒土的蠕应变-时间曲线可分为衰减蠕变和匀速蠕变两个阶段，匀速蠕变阶段试样变形微小，在应力水平较低时迅速趋于稳定。相同围压下，蠕变随应力水平增大而增大，其原因为随着应力水平的增大，土颗粒接触应力相应增大后产生了更多的颗粒破碎，进而导致更大的位置调整，致使蠕变随之增大。这与数值模拟理论分析结果相类似。

图7 数值模拟与室内试验蠕变结果对比图

在接触模型和应力腐蚀模型参数敏感性分析的基础上，本文标定了适用于红砂岩的模型参数标定值，见表2。室内试验与数值模拟所得的蠕应变εc-时间t关系曲线(σ3=300 kPa)对比如图7所示，数值曲线的衰减蠕变时长比室内试验短，后期变形较室内试验偏小，但大体趋势一致。结果表明，基于应力腐蚀理论，采用离散元模拟粗粒土的蠕变是可行的，可基于该程序与参数开展粗粒土蠕变细观机理的进一步研究。

3 结 论

(1) 接触模型参数中，颗粒摩擦系数和胶结键黏结强度对蠕变影响较大，颗粒刚度和密度敏感性较差，粗粒土颗粒母岩表面摩擦系数和抗拉强度越大其蠕应变越小。

(2) 应力腐蚀模型中，激活应力对粗粒土的蠕变有显著影响，随着激活应力的增大，能发生应力腐蚀的胶结键减少，颗粒破碎量急剧减少，蠕应变减小；材料参数β1和β2均对粗粒土的蠕变有较大影响，蠕应变随着β1或β2的增大而增大。

(3) 基于应力腐蚀理论进行考虑颗粒破碎的粗粒土蠕变的数值试验是可行的，并基于粗粒土在不同应力水平下的蠕变室内试验结果标定了接触模型和应力腐蚀模型的细观参数，为粗粒土颗粒破碎蠕变机理的进一步研究奠定了基础。