基于仿真虚拟试验的堆积碎石土蠕变特性

陈 洋，董 辉，余美娟

(1.苏交科集团股份有限公司，南京 211112；2.湘潭大学，湖南 湘潭 411100)

0 引 言

在自然界中，常见的滑坡或斜坡其基体大部分是由于碎石土的堆积而构成[1-3]，因极端条件或实践工程产生滑坡灾害造成经济损害及人员伤亡严重[4]。在风化卸荷作用式、残坡积式、冲洪积式的土石二重介质，即堆积碎石土，在荷载、水、温度、湿度等外界环境下，内部土石体表现出明显的时效性(有流变特性)，例如蠕变破坏、疲劳、长期变形、应力松弛、应变率效应，及蠕变是关于碎石土流变特性的主要表征。碎石土蠕变特性与母岩质量、颗粒形状、基本物性(含水量，含石量及级配)、应力水平等紧密相关，其中在碎石土应力恒定时，其内部土体蠕变分成两阶段[5-7]：第一阶段变形在加载后立即形成，由于颗粒的迁移移位使孔隙率减小，称为瞬时变形；第二阶段持续很长时间，应力传递会导致变形，局部损坏或碎石破碎相继产生，碎石土体积减少(颗粒之间的相对位移、填充作用、结构调整)。目前，基本以试验为主(室内及现场)[8-11]来进行堆积碎石土流变特性的研究，从而得到流变参数并以此建立本构关系。由于碎石土的复杂性与试验周期较长，研究成果存在样本少、个体性强、工况单一等不足，难以满足区域变化、土石含量、块石成分、级配及风化程度等多因素影响的浅层堆积滑坡长期变形特征的理论研究和实践防控。

本文结合颗粒离散元细观方法，试图在解决碎石块体的随机生成、宏观与细观参数标定和蠕变加载机制的基础上，建立仿真试验平台，进行多因素变化下碎石土蠕变特征的分析和研究。

1 堆积碎石土的物理力学特性

针对堆积碎石土，试验取自自然堆积体碎石土斜坡，其位置位于湘潭市湘江东岸昭山地区。其斜坡性质为典型碎石土基体，基体以粉质粘土为主要构成，且内部含有砂质土体，其特征明显，基础母岩为泥质粉砂岩(图1)。

对采集到的昭山堆积碎石土进行取样，烘干、筛分试验后计算级配，查阅相关文献[12-15]，对不同地区的堆积碎石土级配进行分析和对比(图2和表1)，分析得出与其他地区相差极小，其堆积碎石土的研究具有代表性。

图2 碎石土的级配比对：昭山区、室内试验、其他地区

表1 不同粒径碎石含量的对比

根据室内固结不排水试验的试验方法，在制作三轴试样时，其几何特征参数为直径D=101mm，高度H=200mm，采用2～20的等量替代法处理含量为18.2%的超径料(>20mm)，图2为其级配曲线。当碎石土含石量为20%、40%、60%，含水率为7%、9%、11%、13%时，在不同围压下进行碎石土的室内试验(图3)。考虑到碎石含量为60%，含水率为11%的碎石土最接近昭山地区堆积碎石土的实际情况，故以此类堆积碎石土获得的数据对细观参数进行校准，以确定碎石土特性。

图3 试样模型的碎石土剪胀特征

2 蠕变试验虚拟模型的构建

2.1 模型的构建

采用细观的颗粒流(Partical Flow Code，PFC)方法，修改其虚拟双轴压缩试验的编程语言，对块石随机分布的模型进行构建和研究双轴压缩的伺服机制，使模型每循环10 000时步加载200kPa。虚拟试件模型的具体生成步骤分为：(1)调整土颗粒的参数；(2)基于碎石块的粒径变化级配得到信息球；(3)基于信息球随机生成任意碎石块；(4)基于土体的面积并填充土颗粒，模型产生各向同性应力，对材料细观参数赋值。生成的模型如图4所示。

图4 典型堆积碎石土模型的主体构建

为进一步研究堆积碎石土的内部力学特征，针对剪切过程用数值仿真的方法，基于室内试验及数值模拟构建了两个尺寸相异的数值模型。基于室内试验，几何尺寸参数为101mm×200mm，模型1与其一致；模型2的几何尺寸参数依据昭山区碎石土，以300mm×600mm的几何参数尺寸为级配的依据，对各种典型工况进行比较分析，可得到不同的堆积碎石土抗剪特性(见图2级配曲线)。

2.2 任意现状碎石土的随机生成

由于现场碎石的产状特征，其块体结构不规则，对此进行模拟时，采用PFC模型的clump编程模块让圆形颗粒依次连接为聚粒，其连接结构同样为不规则特征，以满足数值模拟的需要。模型典型模式碎石块的生成过程：(1)根据试验分析的典型碎石的级配和粒径、根据软件生成信息球；(2)根据运算程序在信息球随机确定n个(3

(1)

式中:r为模拟信息小球的半径；p为多边形信息小球填充体孔隙率；n为填充体内信息小球的数量；v为生成多边形的面积

图5 任意产状的随机性生成碎石块模型

根据仿真软件随机生成的任意产状模型碎石块与试验碎石形状的一致性是研究的关键。基于试验碎石块与仿真软件模型生成的任意产状颗粒随机进行比对:试验块石右边为产状颗粒簇，其形状变化与外力作用无关(表2)。

表2 模型的随机碎石块形状与室内试验块石的对比

3 堆积碎石土模型细观参数

为更好地体现试验与现场的碎石土不规则特征，在进行数值模拟其力学进程时，需要对其物理力学性质进行定义参数，其细观参数往往决定模型的正确性。采用PFC离散元软件中的平行粘结模型进行数值模拟，当进行校准过程时可选用与室内试验相同尺寸的模型1进行计算，可得到准确试验块石碎石土细观尺度校准参数。为达成以数值模拟模型运算的试验结果，可以对室内试验模型进行准确表征，对模型1的细观参数进行反复修改调节(见表3)，图6为虚拟试验应力-应变曲线与试验的对比(200kPa、400kPa)。经分析可以发现：数值模拟模型所得的结论规律能很好地表征室内试验，根据数值模拟内部模块特征，其模型试验界面参数采用刚性特征参数，在曲线的偏应力差达到峰值之前，主应力差以阶梯的形式上升；由于模型中颗粒的刚性更高，因此轴向应变小于室内测试中的应变[16]。

表3 碎石土模型的细观参数校准

图6 数值模拟与室内试验结果对比(不同围压)

为分析虚拟试验结果与室内试验的误差，绘出两者的莫尔破坏包络线(图7)，比对力学指标(表4)。

表4 数值仿真模拟试验和室内基本试验得到的c、φ值对比

图7 莫尔破坏包络线(虚拟试验与室内试验)

当进行虚拟试验时，其基本物性特征基本符合室内试验(粘聚力和内摩擦角等误差约5%)。受测试过程控制的影响，室内试验中获得的峰值主应力差与在100kPa围压水平下的虚拟测试结果有较大差异(约11%)。其数值模型基本正确，其计算结果误差合理，其结果(试验碎石土基本物性特征)由标定的细观参数可以基本反映。

由于模型2和模型1之间的大小差异很大，基于可靠性与准确性原理，对比分析模型1的细观尺度校准参数，发现其参数差异与模型2基本吻合。由于不同尺度对于模型数值仿真的结果影响较大，为防止几何参数对模拟结果产生较大影响，基于数值模拟的宏微观参数针对模型2进行参数影响分析。为了检验模拟模型校准细观参数的真实性与准确性，增加其应用合理性，通过进行同等摩擦系数(f=0.2)的室内基础试验以及数值模拟(模型2)，计算出其试验材料的c、φ值，并将所得的相应参数与室内测试试验和模型1数值模拟的结果进行比较(表5)。基于表中数据，数值仿真试验模型2得到的c、φ值其中对比来看标准对照误差较大，但数值仿真试验模型2的模型粒径参数和级配等级与原状土石体更加接近。根据尺寸效应使用模型2可以为后续的室内直接剪切试验提供理论依据。

表5 基本剪切试验参数比较(f=0.2)

4 堆积碎石土的蠕变特性分析

4.1 蠕变的试验环境

如本文第2节所介绍，蠕变虚拟试验中的碎石块并不会因外力而发生变形的颗粒簇组成，由此在分析其组合作用时，可只考虑粘滞性分量εp和弹性分量εe的组合作用。虚拟蠕变试验过程中为了增加模型的蠕变特性，同时考虑到现实中碎石土的摩擦系数大多在0.4～0.6之间，可以基于模型2(300mm×600mm)进行分析。当分阶段荷载(f=0.5)下围压为300kPa和400kPa时，各阶段的主应力差分别为200kPa、400kPa、600kPa和800kPa。该研究基于陈氏加载法(图8)来处理测试数据以获得蠕变曲线。

图8 测试数据的分析处理(陈氏加载法)

在测试中，伺服机制确保轴向应力(压盘)和围压(侧墙)恒定。压盘达到指定应力后，此时测试机制的压盘进入一个稳定值，即恒定的加载状态。为保证在下一级加载前模型达到稳定，虚拟仿真试验的稳定持续加载时间为1 000时步。在达到设定恒载时间后压盘加载进入下一个恒载过程,如图9所示。

图9 蠕变试验流程

4.2 轴向蠕变

图10显示了在300kPa和400kPa约束压力下堆积碎石土的分级载荷的轴向蠕变曲线，轴向蠕变特征规律与围压等级的大小成反比例关系。随着偏应力水平的增加，之间的最终轴向蠕变差更大，并且可以获得最终的轴向变形和偏应力，两者接近线性正相关。表6显示了堆积碎石土的平均轴向蠕变变形，在200kPa的偏应力下，在300kPa和400kPa的约束压力下，轴向蠕变的差异仅为0.183kPa，而当偏应力达到800kPa时轴向蠕变差为4.462kPa。

图10 堆积碎石土的轴向蠕变特性(300kPa和400kPa)

表6 碎石土的平均轴向蠕变值(偏应力不同，围压不同)

图10是根据较多使用的陈氏加载方法处理的，可以得出曲线基本趋势相同的,在300kPa和400kPa的围压下分别加载轴向应变-时间曲线(400kPa轴向应变略小于300kPa)和应力-轴向应变等时线(图11)。

图11 300kPa、400kPa等级围压应力-轴向应变等时曲线和轴向应变-时间曲线

分析应变-时间曲线有以下特征：(1)随着轴向荷载逐级施加，结果显示的蠕变曲线主要由两个状态阶段组成，第一阶段为试验衰减阶段，第二阶段为试验等速阶段。在试验时，需要加载某定值偏应力时，当加载时试验模型会立即发生瞬时弹性应变现象，这种现象具有一定的弹性特征，将时间持续进行时，瞬时弹性应变会蠕变应变，产生最终稳定的蠕变变形恒定值。(2)衰减蠕变阶段与偏应力两者之前具有促进关系。例如，在衰减蠕变阶段(约束压力为300kPa，偏应力为200kPa，时间步长为150)时，蠕变变化的速度显著减小，当减小到一定程度时转变为恒定蠕变状态；随着试验进行，当试验负载为偏应力为800kPa时，衰减蠕变状态阶段的总时间步为300时步。(3)偏应力小于400kPa时，瞬时轴向应变极小(土石之间连接的特殊作用)。当偏应力水平超过阈值(600kPa)时，碎石土将经历瞬时变形，瞬时轴向应变急剧增加。

分析应力-轴向应变等时曲线得：(1)围压和应力不同时的应力-应变等时曲线趋于非线性，偏应力400kPa为最大转折点。粘性变形使曲线更接近应变轴，并且时间步长与虚线的斜率呈正相关。(2)偏应力增加使等时曲线偏离原有直线趋势，表明堆积碎石土的蠕变非线性与应力水平成促进关系。

选用最常见的衰减曲线[17]对其变形过程进行拟合。其函数表达式为：

ε(t)=εi+εf(1-e-ct)

(2)

此函数为指数函数，其中εi=σ/E1为瞬时变形；εf=σ/E2为随时间发展的最终变形量；c则为初始相对变形率。

指数关系可以在数学上表示300kPa和400kPa轴向蠕变曲线，其拟合效果很明显。拟合参数和相对误差(表7)与堆积碎石土蠕变曲线与指数函数的函数关系一致，表明最大相对误差为0.33%。

表7 堆积碎石土的轴向应变函数拟合关系与其相对误差

4.3 体积蠕变

根据堆积碎石土体积蠕变曲线(图12)可知：当偏应力较小时，两条围压曲线基本不发生变化；当围压接近300kPa、偏应力800kPa时，主应力峰差为1 003.2kPa，体积应变小，曲线波动明显；相比轴向蠕变，300kPa、400kPa等级围压时其瞬时和最终体积变形两者基本无差别，而且最终轴向应变与偏应力之间存在线性正相关。

图12 堆积碎石土体积蠕变曲线(围压不同，加载分级)

与轴向蠕变等时曲线相比，采用陈氏加载法得到了体积蠕变曲线和应力-体积应变等时曲线(图13)，体积应变等时曲线表非线性明显，尤其是在第50步时呈台阶状，且随偏应力水平升高波动明显。

图13 300kPa、400kPa等级围压应力-体积应变等时曲线和体积蠕变曲线

堆积碎石土体积蠕变曲线的拟合函数关系推导及数值拟合均基于指数函数，且数值表现吻合度较好(拟合参数和相对误差见表8)。

表8 堆积碎石土的体积应变指数函数参数及相对误差

5 结论

(1)基于虚拟试验分析出轴向应变和体积应变曲线为蠕变衰减类型：基于指数函数构建经验本构关系式可精准表现碎石土的虚拟蠕变试验特征，轴向蠕变和体积蠕变的最大相对误差分别为0.33%和0.19%。

(2)随着堆积碎石土试验的进行，对于每级轴向负载时，轴向蠕变状态发展曲线与体积蠕变状态发展曲线均呈现出衰减蠕变阶段与等速蠕变阶段。衰减蠕变时间与偏应力成正比例增长变化，碎石土最终蠕变值与偏应力的促进关系与线性变化的促进关系一致。

(3)碎石土应力-应变等时曲线均呈折线型，转折点位于偏应力400kPa处；堆积碎石土的蠕变随时间的增长呈正相关，其蠕变特征与非线性的特征关系一致，且随偏应力水平的增加脱离直线发展。偏应力较高时体积蠕变曲线其趋势特征规律性，可以表征当时间增加时堆积碎石土内部结构的承载力下降趋势。

本文研究结论是以模型1的三轴试验标定其碎石土细观参数，其精确性仍需结合后续的大三轴室内试验进行提高。碎石仿真是以四种典型形状简化，尽管虚拟仿真体现了总体规律，但从三维角度真实模拟原状碎石土将是下一步的工作重点。