随机散粒岩土材料自然堆积2D数值仿真分析

王海礼 石 崇 王盛年 张玉龙 陈凯华

（1.河海大学 岩土力学与堤坝工程教育部重点实验室，南京 210098；2.河海大学 岩土工程科学研究所，南京 210098）

散粒岩土材料是指由微尺度的固体粒子或粒子团构成的一种介于固相和液相之间的特殊物质组成系统，具有不连续、非均匀、各向异性、流动性等特点.其堆积稳定性问题广泛涉及堤坝、边坡、道路等工程，因而亦是岩土力学研究的基本问题之一［1－2］.

散粒材料在自重作用下降落堆积，因散粒几何特征（如形态、尺寸等）、散粒材料间的摩擦力、散粒材料与堆积垫层间的摩擦力、堆积方式等影响，将形成不同的堆积形态，表现出丰富的散粒体随机结构，所反映散粒自然堆积的休止角亦各不相同，故可从自然堆积休止角在不同因素影响下的变化规律出发研究散粒岩土材料的堆积特性.目前，针对散粒材料堆积方面的研究主要包括室内试验和理论研究两个方面.室内试验主要通过不同人造散粒材料的堆积研究自然休止角的变化规律［3］；理论研究则主要采用离散元法对比分析不同堆积形成方式［4］、孔隙率［5］、干湿［6］情况以及摩擦［7－9］性质下的散粒堆积休止角的变化规律.然而这些研究从材料种类出发或用圆形颗粒代表散粒进行分析，忽视了散粒堆积体内部结构的复杂性，且对于不同形状、级配等符合实际的非均匀散粒堆积分析亦缺乏系统研究.

本文从自然堆积过程的随机性出发，以由任意数量复杂颗粒组成、粒径服从对数正态分布的随机不规则多边形构建岩土散粒，实现了自然堆积过程的计算机数值仿真，并对不同参数特征下的自然堆积进行模拟，研究了不同散粒形态、表面粗糙程度及粒径尺寸等对自然堆积体休止角的影响规律.

1 散粒岩土材料仿真

1.1 散粒随机生成技术

尽管自然界中散粒岩土材料几何形状多样、棱角分明、凹凸有别，但若对其进行简化并采用随机仿真技术仍旧可以模拟该类材料几何形态的多样性和复杂性.目前，现有方法多是以正多边形、椭圆形或凸多边形来模拟散粒岩土块体的几何形态，即通过对选定的“基块体”进行随机缩放或延拓，从而表征粒状材料细观结构特征.但由于该类几何形态忽视了散粒间因凹凸不平引起的咬合作用，模拟仿真结果与实际有较大差异，造成试验结论不准确.

鉴于此，本文采用随机任意多边形模拟散粒岩土材料几何形态.如图1所示，选取长轴尺寸、轴比（长轴与短轴比）均服从正态分布的椭圆作为“基块体”，将其随机分成N个弧段，以弧端点连接形成的多边形代替椭圆，然后对该多边形的每个顶点做其与椭圆中心连线上或延长线上的随机移动，通过调整随机因子λ，则可生成满足要求的散粒块体形态.

图1 岩块随机形态

随机多边形最终顶点坐标可表示为

式中，Ni′（xi′，yi′）为椭圆上第i个顶点的坐标；Ni（xi，yi）为随机 N 多边形第i个顶点的坐标；O（x0，y0）为初始“基块体”中心坐标；a为长轴值；μ为长轴与短轴比；θi为长轴a到第i个顶点与中心坐标连线的逆时针夹角；α为岩块倾角；λ为随机因子.

1.2 自然堆积休止角

散粒材料在自重作用下堆积，其处于极限平衡状态的堆积自由面与水平地面间的最大夹角直接反应了散粒摩擦力或散粒群流动性的大小程度.粒径、摩擦系数越小，形状越趋近圆形，则休止角愈小，流动性亦越好.

散粒堆成形后，由于堆积顶点与投放中线不一定重合，故可通过堆积顶点至散粒堆积外轮廓的水平距离Hi和堆积顶点至地面的垂直距离V比的反余切所得的锐角均值表示自然堆积休止角β的大小，如图2所示，则休止角β为

图2 堆积轮廓线和休止角

2 颗粒离散元法基本原理

PFC2D（2DParticle Flow Code）即二维颗粒流程序［10］，主要用于颗粒材料力学形态分析.它将物体视为由代表性颗粒单元连接形成的组合结构，从而模拟散粒材料的结构问题、介质运动以及其相互间的作用，期望以这种局部的模拟结果来研究边值连续计算问题.该方法以离散单元法为理论基础，将材料力学响应问题从物理域映射到数学域，即采用显式差分算法，循环交替运用力－位移定律和牛顿运动定律更新接触颗粒间接触力和颗粒位置，构建颗粒间新接触，从而模拟散粒材料的运动及其相互作用过程，其计算原理如图3所示.

图3 PFC2D计算原理

2.1 接触力

设B1和B2是两个相接触于O点的均质散粒，其接触点法向矢量n的方向由B1指向B2，则单个时间步Δt内接触点在接触公共面法向和切向的位移增量可表示为

式中，VB2→B1为B2对于B1的相对速度，n为B1指向B2的法向矢量，Δδ为合位移增量，Δδn为接触公共面法向位移增量，Δδs为接触公共面切向位移增量.

根据力－位移定律，则可知散粒B1和B2间的法向和切向力增量分别为

式中，kn为散粒法向平均刚度，ks为散粒切向平均刚度，Sc为接触面积.

故作用在散粒上的总的法向和切向力矢量分别为

当接触过程中切向力大于最大静摩擦力时，散粒将沿着接触公共面产生滑移，因此需对切向力进行修正，即

2.2 运动方程

假设在单个时间步Δt内，散粒间的作用力不变，则根据牛顿运动定理可得到散粒的运动方程：

式中，i＝1，2，3分别表示沿x，y，z三个坐标方向的分量，Fi（t）为t时刻散粒受到的合力，gi为重力加速度在i方向的分量，α为全局阻尼系数，vi（t）为t时刻散粒平均速度在i方向的分量，m为散粒质量，Δvi为相对于前一时间步的速度，Ti（t）为由接触力引起的合力矩，ωi（t）为t时刻散粒的转动速度，I为散粒转动惯量.

由于

故迭代更新后的速度为

于是t时刻散粒的平均位移ui（t）和转动位移θi（t）为

3 散粒自然堆积数值仿真

3.1 散粒堆积模拟模型

为模拟散粒在自重作用下的堆积，需先构建一个散粒的松散堆积结构.该结构中散粒相互之间、散粒与边界之间可有一定距离，但不能存在重叠，以保证在开始堆积前散粒之间、散粒与边界之间不发生相互作用.本文采用一个开口矩形区域作为散粒堆积结构容器，通过在该区域内按散粒粒度级配随机生成相互不重叠、与边界无作用的不规则任意多边形作为散粒轮廓边界，然后将这些散粒轮廓边界与一个模型尺寸相同的隔层错位颗粒模型进行比对，进而寻找出构成散粒的颗粒组，则这些由颗粒组成的不规则形状的颗粒块体即构成了一个松散堆积结构，图4所示为本文所采用的散粒堆积结构模型示意图.

图4 散粒堆积结构示意图

3.2 散粒堆积垫层参数标定

采用不同粒径的人造颗粒进行自重堆积试验，每个试验重复进行3次，试验垫层采用玻璃板，其最终的试验结果见表1.从表中不难看出，同一种垫层条件下，散粒堆积体休止角随粒径的变化相差不大.因此，为了进行数值试验分析，可对垫层表面粗糙度进行标定，通过多次模拟试算，最终标定出玻璃垫层摩擦系数为0.65.

表1 不同粒径自然堆积休止角

3.3 散粒堆积数值仿真结果

采用相同级配下的4种散粒形态（包括圆形、椭圆形、正多边形和多边形）、4种散粒表面粗糙度（摩擦系数取0.2、0.5、0.8和1.2）和3种散粒粒径（粒径介于0.25～1.25cm、0.5～2.5cm 和1.0～5.0cm）进行自重堆积数值仿真分析，得到其相关参数取值如表2所示.

表2 散粒堆积相关参数取值

在进行散粒自重堆积时，为了避免散粒堆积过程中由于下落速度过快对已形成的构形产生冲击而导致自然休止角的减小，散粒堆积结构容器的去除采用低速逐步向上提升的方法，待所有颗粒离开容器后撤去容器，当所有散粒的最终速度均小于10－5m／s时结束计算.图5～7分别是不同形态、表面粗糙度和粒径下的散粒堆积体情况，其中不同表面粗糙程度和不同粒径下散粒的形态由圆形、椭圆形、多边形等随机自由组合，最终得到相应的堆积体自然休止角见表3，表中散粒的粗糙程度采用摩擦系数表征.

图7 不同粒径散粒自重堆积

表3 散粒堆积体自然休止角

4 散粒特征对自然堆积休止角的影响

从图5和表3可知，在相同粒度级配、相同粗糙程度条件下，不同形态散粒自然堆积休止角随着散粒形态复杂化呈逐步提高趋势，其最终构形边界则逐渐呈无规律，不规则变化.散粒形态由圆滑到凹凸不规则变化造成棱角逐渐增多，进而导致散粒间咬合作用力得到加强，散粒运动阻力提高，这是散粒构形外在表现复杂化的原因.而从图6和表3可知，在相同粒度级配、相同形态（圆形、椭圆形、多边形等随机自由组合）条件下，不同粗糙度散粒自重堆积体休止角亦随着散粒粗糙度的增加而增大，且随着散粒粗糙度的增加散粒流动性逐渐降低，伴随有“粒子拱”出现，散粒间空隙率增大.图7为相同级配、形态、粗糙程度条件下的不同粒径散粒自重堆积体休止角变化，结合表3统计数据可知，散粒粒径对堆积体休止角影响则并不明显.

此外，比较表3中不同形态和不同粗糙程度的堆积体休止角的变化可以发现，尽管形态的复杂化使得休止角提高了约25%，但散粒表面粗糙程度增加对休止角却提高近50%，其影响比形态更为明显；散粒粒径变化对堆积体休止角的影响则相对要小一些.因此，可以得出散粒表面粗糙度对自然堆积休止角影响最大，形态次之，粒径最小.

5 结 论

通过采用随机散粒生成技术构建散粒自然堆积数值模型进行仿真分析，其主要结论如下：

1）在相同粒度级配条件下，散粒形态越不规则，自然堆积休止角越大，同时随着散粒表面粗糙度的增加，休止角亦增大，散粒体流动性降低，孔隙率提高.

2）对比分析不同散粒特征下的自然堆积休止角变化表明：散粒表面粗糙度对自然堆积休止角影响最大，形态次之，粒径再次之.

3）通过随机生成技术模拟散粒自然堆积过程，可对松散岩土材料的流动、压实等动力学行为进行分析，为工程建设服务.

［1］ 姜世平，芮筱亭，洪 俊，等.散粒体系统动力学仿真［J］.岩土力学2011，32（8）：2529－2532.

［2］ 吴爱祥，孙业志，刘湘平.散体动力学理论及其应用［M］.北京：冶金工业出版社，2002.

［3］ 赵长冰，刘 军，云 斌.散粒岩土材料自然堆积特性的实验［J］.南京工业大学学报：自然科学版，2010，32（s1）：22－26.

［4］ 王建国，王 葳，赵晓彤，等.利用计算机模拟颗粒自然堆积过程的数学分析［J］.华北电力大学学报，1999，26（4）：19－22.

［5］ 王建国，赵晓彤，胡 斌，等.非均匀颗粒自然堆积过程的计算机仿真［J］.华北电力大学学报，2000，27（1）：33－36.

［6］ Barabasi A，Albert R，Schiffer P.The Physics of Sand Castles：Maximum Angle of Stability in Wet and Dry Granular Media［J］.Physica A，1999，266：366－371.

［7］ 刘 军，刘汉龙.用Monte Carlo方法模拟砂土的自然堆积过程［J］.岩土力学，2005（s1）：113－116.

［8］ Zhou Y C，Wright B D，Yang R Y，et al.Rolling Friction in the Dynamic Simulation of Sandpile Formation［J］.Physica A，1999，269：536－553.

［9］ Zhou Y C，Xu B H，Yu A B.An Experimental and Numerical Study of the Angle of Repose of Coarse Spheres［J］.Powder Technology，2002，125（1）：45－54.

［10］Itasca Consulting Group.Inc.PFC2D （particle flow code in 2dimensions）（Version 3.1）［R］.Minneapolis：Itasca Consulting Group.Inc.，2004.