颗粒材料侧向压力分布特性离散元数值模拟分析

来志强，赵连军，武彩萍，吴国英，张 磊

(1.黄河水利委员会黄河水利科学研究院 水利部黄河泥沙重点实验室，郑州 450003；2.中国科学院 地理科学与资源研究所，北京 100101；3.华北水利水电大学 水利学院，郑州 450045)

0 引 言

水库大坝工程中常见的颗粒材料有筑坝材料堆石体、土体、泥沙等。这些颗粒材料因自重或外荷载作用下会产生侧向压力，在高堆石坝、挡土墙、深埋隧道等工程建设时必须加以考虑。特别是有关土石侧向压力问题的研究已有200多年历史，目前工程上常用的侧向压力计算理论是朗肯和库伦土压力理论。两大经典土压力理论均假定侧向压力沿土石体深度呈线性分布。侧向压力线性分布的假定在国内外挡土墙设计规范[1]、基坑结构稳定分析[2,3]及深埋隧道[4]等工程领域广泛使用。然而土石等工程材料是非常复杂的颗粒材料, 传统的连续介质理论往往忽略了对其离散特性的考虑[5,6]。现有研究表明土石等颗粒材料堆积后将导致拱结构效应的存在，这使得颗粒材料的侧向压力呈现非线性分布[7-9]，这与目前侧向压力线性分布假定相悖。因此，有必要研究颗粒材料侧向压力分布的非线性及其与经典土压力理论中线性分布假定之间的差异，对于目前工程设计具有十分重要的意义。

近年来，离散单元法(Discrete Element Method, DEM)的发展和计算机性能的飞速提高，使得通过细观数值方法研究土石等颗粒材料的力学行为特性成为可能[10-14]。离散单元法模拟土石等颗粒材料时不需要定义宏观本构关系，其复杂的力学响应将由其内部颗粒间相互作用实现，同时不受连续介质理论中连续性和小变形假定的约束。Jiang等[12]采用DEM分析了不同墙体运动模式下刚性挡墙侧向压力分布规律；孟云伟等[15]研究了堆石挡墙位移、应力分布和破坏形式等力学行为；刘冬等[16]和万励等[17]分别探讨了无黏性土体在平动模式下的破坏形式以及土压力的分布特性；Nitka等[18]从细观层面上揭示了挡墙土体剪切带的破坏机制。以上数值研究表明采用离散单元法研究土石等离散颗粒材料的侧向压力分布特性具有独特优势。

本文采用离散单元法DEM对土、石颗粒材料堆积过程进行二维数值模拟，研究了颗粒粒径(土和石)、颗粒摩擦系数对颗粒材料的侧向压力分布特性的影响规律，分析了DEM侧向压力分布特性与经典土压力理论解的差异，基于分形理论探明了不同颗粒粒径和颗粒摩擦系数下侧向压力强度与作用位置分形分布特征，为进一步完善挡墙、基坑等工程设计理论提供参考依据。

1 侧向压力分布DEM数值模拟分析

1.1 DEM数值计算模型

采用基于DEM理论的离散元软件PFC2D模拟颗粒材料自然沉积过程：首先分别生成左、右、下三道刚性墙体，在挡墙内随机生成10 000个粒径为d的圆形颗粒，并在自重作用下堆积至平衡状态。颗粒最终的堆积高度为H0，堆积宽度为L0，如图1所示。考虑到计算效率，因此采用能够较合理反映土石力学性质的单一粒径分布的圆形颗粒。

图1 DEM数值计算模型Fig.1 DEM numerical simulation model

参照DEM模拟土、石的相关文献[15,16]，颗粒密度取为2 600 kg/m3，颗粒法向和切向刚度取为1.0×108N/m，墙体刚度与颗粒刚度保持一致，墙体摩擦系数选取为0以便于数值结果与朗肯土压力理论解作对比；本文分别选取石块、土体颗粒粒径d为0.40 m、0.01 m，在保证其他数值计算参数不变的条件下，通过改变颗粒粒径d和颗粒摩擦系数μ以研究两者对颗粒材料的侧向压力分布特性的影响规律。本文最终数值模拟参数选取如表1所示。

表1 DEM数值计算参数[15,16]Tab.1 DEM numerical simulation parameters

1.2 DEM侧向压力分布

为研究颗粒粒径d、颗粒摩擦系数μ对颗粒材料的侧向压力分布特性的影响规律，分别模拟d为0.40 m，0.01 m和μ为0.30、0.60、0.90条件下颗粒堆积过程，其他计算参数均按表1所取。当不同d、μ条件下颗粒堆积稳定后，采用PFC2D中内嵌FISH语言统计颗粒与左侧墙体接触力链的强度及其位置，得到DEM侧向压力分布如图2所示。可以看出，不同d、μ条件下DEM侧向压力沿深度均呈锯齿状分布，这与经典朗肯土压力理论中侧向压力沿深度呈线性分布的假定不一致。

图2 不同颗粒粒径d、颗粒摩擦系数μ下DEM侧向压力分布Fig.2 Lateral force distributions obtained from DEM with different particle size d and friction coefficient of particles μ

图3 DEM力链分布中的拱结构效应Fig.3 Arch structure effect in the force chains of DEM

颗粒在重力作用下相互接触形成力链，图3给出了d=0.01 m、μ= 0.30条件下颗粒堆积体力链分布图(右侧)和拱结构效应图(左侧)，图中力链的粗细代表了接触力强度的大小。强度较大的力链构成了颗粒集合体承载外力的拱结构，对颗粒系统的稳定起决定性作用。拱结构效应很好地解释了DEM侧向压力分布呈锯齿状的原因。当与墙体有接触的颗粒为拱结构组成颗粒时(如图3中颗粒A、D和E)，强度较大的力将通过这些颗粒传递给墙体；当颗粒没有位于拱结构时(如颗粒B和C)，强度较小甚至为0(如颗粒C)的力将传递给墙体，因此颗粒材料侧向压力在拱结构效应的影响下呈现锯齿状分布。

2 不同颗粒粒径和摩擦系数下侧向压力理论解和DEM数值解

2.1 朗肯土压力理论

经典朗肯土压力理论认为静止状态下墙后填土侧向压力分布沿深度呈线性分布，可表示为：

f=K0ρgH

(1)

式中：f为侧向应力分布强度；K0为静止土压力系数，H为堆积深度。

静止土压力系数K0可依据文献[19]定义如下:

K0=1-sinφ

(2)

式中：φ为颗粒集合体内摩擦角。

根据文献[20]可知，颗粒集合体内摩擦角φ约等于其休止角度。本文通过DEM分别模拟了μ为0.30、0.60、0.90条件下颗粒集合体在无侧限约束下的堆积过程，得到了相应的休止角度(见表2)。由于DEM采用圆形颗粒，没有考虑颗粒形状咬合作用，使得表2中数值计算的内摩擦角略小于真实物理试验测得的。静止土压力合力Fr可表达为：

(3)

根据朗肯土压力理论可知合力作用点与左侧墙底距离为Lr=H0/3。依据上式可得不同d、μ条件下朗肯土压力理论解，具体如表2所示。

表2 不同颗粒粒径d、颗粒摩擦系数μ下DEM和朗肯土压力理论侧向压力对比Tab.2 Comparison of lateral force distributions obtained from DEM and Rankine theory with different particle size d and friction coefficient of particles μ

2.2 理论解和DEM数值解对比

由表2可知，DEM侧向压力合力值小于朗肯土压力理论解，当颗粒粒径d较大时，两者的差异显著减小；DEM侧向压力合力作用点位置与理论解差异很小。随着颗粒摩擦系数μ的增大，颗粒之间由摩擦提供的咬合作用力增强，颗粒集合体内部成拱率增加，导致DEM侧向压力合力增大；由图2可知，DEM某些部分侧向压力值较周围值骤增(图中红色圆圈所注)，这是朗肯土压力理论无法计算得到的，工程设计时应考虑这些危险点的存在，例如采用DEM数值模拟相同工况以得到骤增侧向压力值及作用位置，在挡墙相应位置处进行锚固、灌浆等加固处理。

3 DEM侧向压力分布分形特性研究

分形理论目前被广泛地应用于岩土工程领域中以定量描述表面复杂无序而又内在具有规律的系统。传统的欧几里得几何学定义任意形状直线的维度为一维，平面的维度为二维。分形理论则采用非整数的形式来描述任意形状直线(折线)的分布状态。本文采用分形理论量化DEM侧向压力的分布特性，研究不同d、μ条件下DEM侧向压力强度、作用位置分布内在规律。

3.1 DEM侧向压力强度分布分形特性研究

根据分形理论的定义可知[21,22]，当采用尺度为r的尺子量取总长度为L的任意形状直线(折线)时有如下关系式：

L=Nr

(4)

式中：N为该直线长度所包含尺度r尺子的数目，如图4所示。

同时，r、L和该直线的分形维数D之间具有如下关系式[20,21]：

L=Cr(1-D)

(5)

式中：C为常数。

由式(4)及式(5)可得：

N=Cr-D

(6)

通过拟合logr与logN之间的线性关系得到斜率即分形维数D。

参照图4中量取任意形状直线(折线)的方法，分别统计出不同d、μ条件下DEM侧向压力强度分布(图2中的折线)所包含不同尺度r的数目N(见表3)，并对其分形维数D进行拟合，如图5所示。分形维数D越大，表明DEM侧向压力强度分布的波动性越大。由表3和图5分析结果可知，当d、μ发生改变时，DEM侧向压力分布强度分形维数D变化较小，分布在1.17～1.29范围。这表明DEM侧向压力分布的波动性受d、μ的影响很小，可能是颗粒材料的本质特性。

3.2 DEM侧向压力作用位置分布分形特性研究

颗粒材料的离散特性使得DEM侧向压力作用位置分布具有非连续性。采用分形理论研究不同d、μ条件下DEM侧向压力作用位置分布特性。图6给出了分析方法：图中黑色直线代表长度为20 m刚性墙体，蓝色短线代表刚性墙体上侧向压力作用位置；间距为r的红色短线将刚性墙体平均分成数目为20/r(20/r若为分数，则四舍五入取整)个子组，N为子组中包含作用点(蓝色短线)的子组数目。DEM侧向压力作用位置分布分形维数D值在0～1之间，0表示作用位置分布几乎集中在某一点处，1则表示作用位置分布为线性连续分布。D越大表明侧向压力作用位置分布越分散均匀。

图4 不同尺度r下任意形状直线(折线)分形分析方法示意图Fig.4 Schematic diagram of fractal analysis method of arbitrary shaped lines with different ruler r

表3 不同颗粒粒径d、颗粒摩擦系数μ下DEM侧向压力强度分布分形分析
Tab.3 Fractal analysis of lateral force intensities obtained from DEMwith different particle sizedand friction coefficient of particlesμ

颗粒粒径、摩擦系数N (r=0.1)N (r=0.2)N (r=0.3)N (r=0.4)分形维数 Dd= 0.01 m、 μ= 0.30139.1160.1230.0526.781.23d= 0.01 m、 μ= 0.60192.2692.9549.0232.531.29d= 0.01 m、 μ= 0.90220.5199.4060.1939.551.23颗粒粒径、摩擦系数N (r=2.0)N (r=8.0)N (r=16.0)N (r=32.0 )分形维数 Dd= 0.40 m、 μ= 0.30510.08107.6442.1019.561.19d= 0.40 m、 μ= 0.60362.9767.6024.5010.961.28d= 0.40 m、 μ= 0.90423.4985.3736.9620.131.17

图5 不同颗粒粒径d、颗粒摩擦系数μ下DEM侧向压力强度分布log r- log N线性拟合Fig.5 Linear fitting relationship between log r and log N for lateral force intensities obtained from DEM with different particle size d and friction coefficient of particles μ

图6 不同尺度r下侧向力作用位置分布分形分析方法示意图Fig.6 Schematic diagram of fractal analysis method of lateral force positions with different ruler r

参照上述方法，表4给出了不同d、μ条件下DEM侧向压力作用位置分布分形特性分析结果，通过拟合表4中logr与logN的线性关系(见图7)，根据斜率即可求出D。由表4和图7可知，DEM侧向压力作用位置分布具有良好的分形特性；当颗粒粒径d较大时，D较小，表明侧向压力作用位置分布较集中；当颗粒摩擦系数μ增大时，D减小，这表明侧向压力作用位置分布整体倾向于集中密集，刚性挡墙的某些部位将承载更多的侧向压力，这将易导致挡墙局部产生破坏，不利于其安全稳定。

表4 不同颗粒粒径d、颗粒摩擦系数μ下DEM侧向压力作用位置分布分形分析Tab.4 Fractal analysis of lateral force positions obtained from DEM with different particle size d and friction coefficient of particles μ

图7 不同颗粒粒径d和颗粒摩擦系数μ下DEM侧向压力作用位置分布log r- log N线性拟合Fig.7 Linear fitting relationship between log r and log N for lateral force positions obtained from DEM with different particle size d and friction coefficient of particles μ

4 结 论

(1)经典朗肯土压力理论假定土石等颗粒材料为连续介质，有悖于其离散特性的本质，计算得到侧向压力分布随堆积深度均匀增大且连续分布。由于颗粒材料堆积后内部呈现的拱结构效应，DEM模拟得到的侧向压力分布呈锯齿状且具有不连续性，部分侧向压力值较周围值骤增，在工程设计时应结合相同工况下DEM数值模拟分析得到的骤增侧向压力值及作用位置，在挡墙相应位置处进行锚固、灌浆等加固处理。

(2)DEM侧向压力合力值小于朗肯土压力理论解，当颗粒粒径增大时，颗粒侧向压力合力与朗肯土压力理论解差异显著减小。两者计算得到的合力作用点位置几乎没有差异。

(3)采用分形理论分析了不同颗粒粒径和颗粒摩擦系数条件下DEM侧向压力分布特性。侧向压力强度分布分形维数受颗粒粒径和颗粒摩擦系数影响不大，分布范围在1.17～1.29。侧向压力分布位置受颗粒粒径和颗粒摩擦系数影响较大；颗粒粒径和颗粒摩擦系数越大，侧向压力作用位置分形维数越小，作用力分布位置越密集，越不利于刚性挡墙的安全稳定。

□