岩质颗粒压缩破碎过程的离散元数值模拟研究

，，，，

(1.许昌学院 土木工程学院，河南 许昌 461000；2.河海大学 力学与材料学院，江苏 南京，210098)

天然岩质颗粒材料是土木、交通、水利等领域广泛应用的一种建筑材料，如碎石基础、铁路宕渣、土石坝等.随着我国基础建设的快速发展，工程规模越来越庞大，目前难以完全应用胶凝材料[1]，当颗粒在较大外力作用下会产生破碎现象，造成颗粒材料级配发生变化，影响到颗粒骨架结构、咬合程度、粒间摩擦等一系列物理性能，改变材料的变形和剪切强度[2]，关系到工程的安全使用.如巴西坎泼斯诺沃斯面板堆石坝，设计高度达202 m，由于未充分考虑坝体底部颗粒破碎的影响，导致坝体因变形过大而发生破坏的现象[3].

天然岩质材料是由多种矿物晶粒、胶结物和孔隙组成的复合材料，内部含有不同阶次随机分布的微观裂纹和孔隙，宏观尺度上又被各种地质构造面所切割，如断层、弱面等，造成岩质材料存在着宏-细-微观上的非连续性.因此，岩质材料在宏观上接近连续但又非完全连续，微观上近乎破断但又非完全破断的复杂材质，从而使岩质材料具有非线性、随机性、流变性等复杂的力学行为[4].

离散元数值模拟方法可以根据材料的属性，能提供详细的细观动态参数、力链、位置的变化等信息，为颗粒碎裂方式、破坏机理、级配变化等问题提供有效的研究手段[5].其中将离散的基本颗粒单元，通过黏结的方式组合成可断裂的组合单元，模拟颗粒破碎的过程，是颗粒破碎问题研究的重要方法，如胶凝颗粒模型(Bonded-particle Model，简称BPM).Thornton等[6-9]首先采用可破碎的BPM单元，对冲击破碎展开了一系列研究，发现BPM单元密实度对颗粒的断裂方式和破坏方向有直接影响.Mcdowell等[10,11]为了得到单颗粒破碎强度与数理统计Weibull分布规律相吻合，采用了删除BPM单元中主要颗粒的方法.并进行了砂样的一维压缩试验模拟分析，发现得到的屈服应力比试验中石英砂小[12]，随后大量学者采用该类破碎组合模型展开了一系列研究.但是，没有对颗粒组合的细观结构进行详细的探讨，该破碎组合模型从细观结构上吻合含孔隙率的多孔岩质材料，而大部分颗粒属于密实结构材料，因此模拟颗粒压缩破碎时会释放过多的孔隙率.Nguyen等[5]为了解决释放多余孔隙率的问题，采用多边形凝聚单元(BCM)模型对巴西圆盘进行了研究，但又产生新的问题，如无法检测破碎时的裂纹、单元划分不当造成内部应力集中等. 综上所述，为了全面的研究岩质颗粒的细观破碎机理及对宏观现象的响应等重要问题，BPM组合单元仍需要不断的细化和改进.

本文基于离散元方法对岩质颗粒压缩破碎过程进行数值模拟分析，数值模型的建立和计算采用离散元软件PFC2D，主要目的是展示应用离散元方法，建立颗粒压缩破碎试验仿真计算模型的方法，以及通过数值模拟得到仿真试验的数值结果，归纳颗粒破碎每个阶段的特征.

1 离散元模拟BPM单元的基本原理

颗粒流(PFC)[13]，为了提高计算效率，将每个颗粒简化成一个圆形颗粒.为了模拟宏观颗粒的破碎过程，可采用凝聚颗粒单元法(BPM)模型，构建密实状态结构的颗粒破碎单元，将岩质材料的矿物晶粒等效为细观单元，将连接方式等效为黏结模型.

PFC中有两种接触刚度模型：线性接触和非线性Hertz-Mindlin接触，同时提供了刚度本构模型编译接口，可以根据需要编译新的本构关系.PFC中有两种黏结模型：接触黏结和平行黏结.选用点接触的接触黏结模型模拟，细观单元间刚度模型选取线性接触.

线性接触模型是通过串联的方式定义两个相互接触球形颗粒，接触法向刚度由下式计算：

(1)

接触切向刚度Ks由下式计算：

(2)

对于线性接触模型，法向割线刚度模量kn等于法向切线刚度模量，即

(3)

式中：Fn和Un分别为两球体颗粒法向作用力和位移重叠量；Kn由公式(1)给出.

如果法向拉力大于法向黏结强度，黏结发生破坏，法向和切向接触力都为零.而当切向作用力大于切向黏结强度时，黏结也会发生破坏，但接触力不会发生改变.如果切向作用力没超过最大摩擦力，法向力是压力.在颗粒接触点处相对位移和法向、切向接触力分量之间的关系如图1所示：

图1 接触点上本构特性

2 离散元模拟的基本过程

2.1 细观参数的标定

对于离散单元法，细观参数的标定是复杂和繁琐的过程.主要原因在于细观参数的定义和尺度与宏观参数存在很大的差异性，对于大部分离散材料，很难用明确的数学关系式表达宏-细观参数之间的关系.如法向刚度kn、切向刚度ks等，影响材料属性的参数非常多，因此需要合理的选择标定方法.

对于确定模型细观物理参数值主要分为理论公式法和宏观参数反演法.理论公式法参数标定主要适应于规则的连续介质材料[14].而对于大部分材料，尤其离散颗粒材料，则需要采用反演法完成细观参数的标定.实际中应先确定模拟材料的力学行为，然后与现场试验结果进行对比，通过不断对参数反复进行校核，直至达到与真实材料力学效应逼近的效果，确定合理的参数值[15].

2.2 宏观颗粒模型的建立

对于离散元数值模拟过程中，计算时间是需要考虑的一个重要问题.针对密质结构材料，新的BPM单元的构成细观单元具有两方面的特点：粒径相同和最密实排列，该结构避免了颗粒破碎过程中BPM单元内部孔隙率的多余释放，并极大提高BPM单元的生成效率.由于几何关系的影响，采用在一个颗粒周围依次增加颗粒数目来生成预定的BPM单元，会积累越来越大的误差，并且生成形状各异的宏观颗粒非常困难，无法满足研究需要.为了解决以上问题，提出了一种新的BPM单元生成方法.步骤如下：

图2 宏观颗粒生成示意图

(1)首先，根据要生成BPM单元的尺寸自动设置足够多细观单元的列和行数.生成紧实颗粒组合(孔隙比e=0.102 6)如图2所示，几何关系为：

(4)

式中：e为孔隙比；Aa为孔隙面积；A为总面积；As为固体面积；xi为横坐标；yi为纵坐标；ri为细观单元半径.

可以通过坐标转换(公式(5))，生成具有一定倾角α的试样，满足内部结构变化要求，内部结构为倾角的宏观颗粒.

(5)

(2)根据预定设置的列和行，根据需要自动检索定位每一行单元的坐标，根据需要确定每一行的单元数目，依次检索每一行，完成模型试样.

2.3 宏观颗粒压缩破碎过程的数值模拟

图3 颗粒压缩示意图

单颗粒的压缩破碎强度分析是研究颗粒破碎问题的重要内容，由于颗粒属性的千差万别，同时又受加载方式、形状、尺寸、应力状态的影响，造成颗粒破碎强度值离散性非常大.岩质材料破碎主要为脆性破坏，即使同种材料同样粒径的情况下，测得的强度值有时也会超过平均强度值的倍，偏差远远大于弹塑性连续介质材料[16].

为了分析天然颗粒破碎的过程，采用上节方法生成27行由517个单元组成的内部密实，外部形状为多棱的宏观颗粒，加载前内部排列具有一定夹角，如图3所示.采用上节反演法确定单元的细观参数，最终得到细观参数，详见表1.得到力-位移曲线与文献[17]一致(图4所示)，破碎最大作用力为7.36 kN，和试验结果7.34 kN, 误差控制在0.3%以内.

表1 单颗粒压缩模型参数

图4 力-位移曲线对比图

离散元方法可以根据需要，通过Fish语言编译，并调用Fishcall程序，监测黏结断裂的位置、坐标、数目、类型等细观数据，分析颗粒破碎过程的变化规律.对于该类颗粒，压缩破碎的基本过程可归纳为四个阶段，如图5、6所示：

(1)预压阶段.如图5(a)所示，主要针对不规则形状的颗粒，由于平整度不够或者摆放位置的影响.在该阶段与加载扳的接触一般为点接触，接触点部位作用力相对较大且不均匀，一般不发生颗粒破碎或少量破碎.使摆放不平或接触面偏小的颗粒发生转动，主要特征是；作用力小，位移很大，而规则的颗粒或内部强度低的颗粒没有此阶段.

(2)稳定加载阶段.如图5(b)所示，颗粒与加载边界的加载方式转变为面加载，颗粒内部作用力分布均匀，并相对于颗粒中心对称.对于不规则颗粒，随着加载作用力的增大，会发生局部破碎现象.并且伴随少量的法向力链断裂，切向链断裂数目呈线性增加.对于规则形状的颗粒力-位移一般成线性关系.

(3) 颗粒破碎阶段.如图5(c)所示，该阶段从接触强度观察到法向力链断裂数目快速增加，裂纹上下扩展贯通，对于不同属性颗粒破碎方式有所不同.加载作用力达到最大值，颗粒断裂.破碎的细观原因是由于法向黏结断裂断裂数目的突变.

(4)破碎完成阶段.如图5(d)所示，颗粒分裂成多个碎块或散体颗粒，颗粒承载作用力急剧降低或变形加速变大.黏结断裂增加到一定数目达到稳定值，该单颗粒已完全破碎.

图5 加载过程破碎示意图

图6 接触点破碎数目-位移曲线

通过离散数值方法模拟颗粒破碎过程，通过反演法对细观参数进行标定，能够得到和试验相吻合的宏观响应；从细观角度能很好的反应颗粒破碎的过程，并发现颗粒的法向黏结断裂数目的突变是颗粒承载力失效的重要特征，而切向黏结断裂数目，在宏观颗粒完全失去承载力之前基本保持线性增长.

3 结论

应用离散元方法建立了模拟岩质颗粒的数值模型，实现了较好的仿真模拟.数值模拟可以监测压缩破碎过程任意阶段的作用力、黏结断裂位置、类型、数量等物理量的变化过程，能弥补室内试验限于技术原因无法得到某些需要的物理量的不足，提供强大的数据支持.

仿真模拟能得到与室内试验相一致的力-位移曲线及相应的压缩破碎阶段.

细观颗粒的法向黏结断裂数目的突变是宏观颗粒承载力失效的重要特征.