颗粒堆积体各向异性及宏细观力学特性的三维离散元模拟研究

莫品强，赵子露，胡裕琛，李国耀，林浩东

(中国矿业大学 a.深部岩土力学与地下工程国家重点实验室，b.力学与土木工程学院，徐州 221116)

颗粒物质广泛存在于自然界，其来源、形状、尺寸等特性对人类生活生产影响重大。然而，颗粒材料的运动学及力学问题尚未解决，被认为是当今世界最前沿的科学问题之一。颗粒系统的复杂性来源于颗粒间的几何组构、接触力、接触摩擦及丰富的集体运动现象；颗粒间的细观机理是揭示宏观力学现象的基础，而小尺度的颗粒又难以通过推理演绎预测大尺度的集体无序运动。在工程领域中，砂、石、煤等建筑或能源材料均属于颗粒材料，在我国经济建设中担任着重要角色。这些材料由于沉积作用而呈现各向异性特性，不同初始各向异性的颗粒体系在堆积过程中展现不同的力学特性，在工程中存在诸多安全隐患。例如，近年来一些大型储煤场煤堆滑坡事故频发，造成了严重的生命财产损失。因此探究初始各向异性对颗粒材料堆积体特性的影响有助于进一步揭示颗粒物质的宏细观机理，且对降低工程风险有着重要作用。

颗粒堆积体的休止角是评价堆积体稳定性的重要指标。学者们通过试验、数值模拟等方法[1-2]研究了堆积体休止角的影响因素[3]。颗粒堆积体的基底压力分布存在“中心压力下降”现象(pressure dip)，且堆积体的形成方式对该现象具有决定性作用[4]。此外，研究表明颗粒形状、沉积高度等因素也对压力下降现象有一定的影响[4-6]。随着试验技术和数值模拟的发展，学者们结合颗粒尺度机理分析研究了力链网络、颗粒接触方向、颗粒接触力分布等规律，认为力学各向异性是颗粒堆积体底部压力下降现象出现的重要原因；已有研究发现颗粒堆积体基底压力下降现象与方向相关的组构各向异性有关[7-8]。然而，目前初始各向异性程度对颗粒堆积体特性的影响尚缺乏系统研究。

本文将模拟不同围压固结状态与三轴排水剪切破坏状态下形成的具有不同初始组构的颗粒堆积体，通过分析其休止角与基底压力变化规律，结合颗粒系统细观力链分布及组构演化规律，建立宏细观联系，探究不同初始各向异性程度对颗粒堆积体宏观特性的影响。

1 三轴剪切试验

1.1 三轴压缩建模与参数设置

本文基于南京大学自主研发的岩土体高性能GPU矩阵算法离散元软件MatDEM[9]建立代表性体积单元的三轴数值试样，试样边界采用由颗粒组成的刚性墙，边界充当压力板，颗粒尺寸为0.14 mm；试样初始尺寸为0.016 m×0.016 m×0.016 m，随机生成颗粒，颗粒数量为5 289个。本研究中采用颗粒级配类似于Ottawa 20-30 sands[10]，颗粒最小直径0.5 mm，最大直径1.0 mm，平均粒径d50=0.725 mm，其颗粒级配曲线如图1所示，图2为三轴排水剪切模型。颗粒材料的微观力学参数如表1所示[10]。试样颗粒摩擦系数设为0.3，边界板摩擦系数设为0，颗粒的正向刚度和切向刚度设置为3.0×105Pa，由于线弹性模型在球形或非球形颗粒介质的微观和宏观力学行为方面与非线性Hertz-Mindlin模型相似[10]，线弹性模型可以用来研究颗粒土的微观力学行为，因此本模拟采用线弹性接触模型。

图1 试样的颗粒级配曲线Fig.1 Particle size distribution curve

图2 三轴排水剪切模型Fig.2 Triaxial drainage shear model

表1 颗粒材料的微观力学参数Table 1 Material properties of particles

试样的等压固结由0.01 MPa开始，围压σ0依次递增0.01 MPa，采用标准平衡直至稳定，分别得到0.1，0.5，1.0，1.5，2.0 MPa等压固结状态试样。采用应力张量σ表示固结和剪切过程中颗粒体系的宏观力学响应：

(1)

式中：V是试样的总体积，Fc是接触点C处的接触力矢量，lc是接触点C处两个接触颗粒中心的分支矢量。三轴试样平均应力p和偏应力q由应力张量σ定义如下：

(2)

(3)

式中：tr为张量迹运算，dev为偏张量运算。

为得到均匀致密的固结样本，同时加快模拟速度，将固结过程中颗粒间的摩擦系数设为0.参考文献[11]中的平衡标准，本文固结过程中平衡标准设置为：平均应力比αf≤0.025，偏应力比βf≤0.025，不平衡应力比γf≤0.001；其中，各参量的定义如下：

(4)

(5)

按照上述固结平衡后得到5个等压固结试样，固结后的状态参数如表2所示，根据试样颗粒间接触，配位数定义为与所考虑的颗粒相接触的颗粒数量，平均配位数〈Z〉是指颗粒组分中颗粒的平均接触密度，并记录颗粒体系中的平均接触力〈f〉.

表2 等压固结后的状态特性Table 2 Characteristics after isotropic consolidation

所有试样在固结后进行三轴排水剪切试验，整个剪切过程中摩擦系数设为0.3，采取位移控制方法，上板每次向下移动1.4×10-3mm，下板保持不变；而其他4个侧板单独移动，以保持恒定的围压0.1、0.5、1.0、1.5、2.0 MPa的应力控制伺服；同时控制轴向总位移为初始固结状态高度的50%.本模拟采用惯性数作为指标衡量准静态加载，惯性数为：

(6)

式中：εz是加载应变率，定义为v/H0，v为轴向加载速度，H0为剪切前的试样高度。

1.2 结果分析

跨尺度关联分析一直是颗粒研究方面重要的分析方法，是颗粒多尺度研究的关键，孙其诚等[7]提出颗粒物质力学研究框架(如图3所示)，指出颗粒性能参数决定力链形态及其稳定性，而力链网络的复杂动力响应决定了宏观土体的力学性质。因此从颗粒体系组构、力链网络的演化等细观尺度揭示宏观力学性能是必要且合理的。

图3 颗粒物质力学研究框架Fig.3 Research framework for mechanics of granular matter

1.2.1宏观应力-应变关系

三轴剪切试验的轴向应变εz和体积应变εv通过大变形假设计算，即：

(7)

式中：H和V分别是剪切过程中试样的高度和体积。

图4显示了所有试样偏应力比q/p和体积应变εv随轴向应变的变化。由图4(a)可以看出轴向应变达到20%以后偏应力比产生了一定差值，这是由于本模拟不同围压条件采用固定的速率进行剪切，虽满足准静态加载但惯性数有所不同，同时颗粒样本较少，模型样本存在一定不均匀性，如内部与边界处应力状态不一致等现象。当轴向应变达到30%时，偏应力比和体积应变都趋于稳定值，可以认为所有试样在轴向应变为40%时达到临界状态，因此本文取轴向应变为40%～50%这一阶段作为临界状态区间进行观察。

图4(a)给出了5个试样的应力路径图，分别标记了剪切前初始点、剪切过程峰值点(偏应力比到达峰值时，见图4(a)和临界状态点。其中，临界状态点取剪切轴向应变40%～50%的平均值，并通过直线拟合qc=Mcpc，获得了试样的临界内摩擦角φ=36.8°，与文献[12]中Ottawa干砂在中高应力水平下的内摩擦角φc=36.5°近似。无黏性颗粒材料孔隙比的临界状态线可以在e-(p/pa)x空间上近似线性绘制，其中pa是参考压力，取标准大气压1.01×105Pa，x是材料参数。LI et al[13]在对大量研究数据进行分析后，得出砂土临界状态线表达如下：

图4 三轴压缩试验结果Fig.4 Results of triaxial compression tests

(8)

式中：ec是临界状态时的孔隙比，eΓ是p=0时的孔隙比，λc是斜率，当x≈0.7时对拟合直线的影响较小[13]。图5(b)是e-(p/pa)0.7空间上的应力路径和拟合得到的临界状态线。

图5 应力路径与临界状态线Fig.5 Stress paths and the critical state line

1.2.2细观组构演化规律

颗粒间的接触被广泛认为是决定颗粒力学性质的重要因素之一[17]，颗粒之间的接触构成整个试样体系的接触网络。图6是2.0 MPa三轴压缩下试样的接触网络，由于对三维模型接触网络分析较困难，因此将接触网络投影至x-z平面进行分析。

图6 2.0 MPa三轴压缩接触网络Fig.6 Contact network for triaxial shear test under 2.0 MPa confining pressure

由于试样颗粒间接触过多，在图6中展示的数值是接触与平均接触力的比值，进行了归一化处理，图中颜色越暖、网络越粗代表接触力越大。可以看到三轴压缩初始状态下的接触力非常均匀的分布在整个试样，且接触力显示出各向同性的力学性质；峰值状态可以明显看到竖向接触力增大，颗粒间产生了竖向的各向异性，这是三轴压缩竖向位移作用的结果；临界状态时的接触网络出现了拱状接触，这是由于三轴压缩破坏时颗粒接触的重新分布，使接触力的方向发生了偏移，与附近接触偏移的颗粒形成了新的接触，从而形成拱状结构，也称为应力拱。其他围压下的接触网络与2.0 MPa相似不予显示。

颗粒集合体中两个接触颗粒间的相互作用可以通过接触法向以及垂直于接触平面的法向接触力和平行于接触平面的切向接触力来描述。离散单元法的一个重要优势就是能够详细地刻画出颗粒体系内部细观接触法向和接触力的分布情况，从细观角度剖析揭示宏观力学机制。

几何各向异性通常由接触法向和枝向量各向异性组成。本文采用ODA[15]提出的组构张量来量化接触法向：

(9)

(10)

(11)

(12)

力学各向异性主要由外力引起，取决于与接触平面方向相关的接触力，可分为法向接触力各向异性和切向接触力各向异性，分别定义如下[16]：

(13)

(14)

(15)

(16)

(17)

图7 接触法向分布Fig.7 Distribution of contact normals

类似地，图8是在不同状态下的法向接触力分布。有研究[17]表明法向接触力对偏应力的影响远远大于切向，因此着重关注法向接触力。在初始状态下，颗粒体系中法向接触力分布在各个方向是均匀的，法向接触力各向异性系数接近于0，并且随着围压的增大，法向接触力增大，而各向异性系数基本不变。法向接触力的各向异性程度随剪切增大，且临界状态的各向异性系数亦随围压增大而减小。

图8 法向接触力分布Fig.8 Distribution of normal contact forces

各向异性组构被认为是承担偏应力的主体，为颗粒体系提供抗剪强度，其与抗剪强度间存在的相关关系也被称为应力-力-组构(SFF)关系。GUO et al[16]针对球形颗粒体系提出了采用各向异性系数α表征偏应力比的表达式：

(18)

本文采用离散元宏细观数据验证了该SFF关系，如图9所示，可以看出各向异性和剪切强度之间存在很好的相关性。

图9 应力-力-组构关系Fig.9 Stress-force-fabric relationship

2 初始各向异性颗粒材料堆积试验

2.1 堆积试验建模与模拟方案

本文利用三轴排水剪切试验结果，探究不同初始各向异性程度条件对颗粒材料堆积体的休止角和基底压力分布规律的影响。采用边界去除法对上节中等压固结和剪切破坏的10个不同初始各向异性试样进行堆积试验。在LIU[14]的研究中提升速度介于慢速2.5 cm/s与快速7.0 cm/s间，发现提升速度越快，休止角越小，因此本模拟选取中等速度4.8 cm/s提升试样四周及上表面边界，让颗粒在重力作用下堆积。为加快模拟速度，将阻尼系数设置为0.01.颗粒在形成堆积体后，用不平衡应力比准则≤0.001来确定堆积的平衡状态[8]。

图10为边界去除试验堆积最终状态俯视、侧视及力链图。图10(c)中可以明显看到力链集中于堆积体内中心区域。颗粒与颗粒之间存在接触力，将颗粒之间接触力大于平均接触力的记为强接触，反之为弱接触。对力链强度归一化处理后，力链线条越粗颜色越暖代表力链越强，力链分布表明颗粒体系中弱力链个数远大于强力链，强弱力链共同构成了一个稳定的颗粒体系网络。

图10 试验最终状态及力链图Fig.10 Final state and force chain diagram

2.2 结果分析

2.2.1堆积体休止角与基底压力分布

本文针对每个堆积体的8个侧视方向进行分析，利用MATLAB图像处理技术获取单侧休止角，如图11所示，取其平均值作为该堆积体的平均休止角。

图11 MATLAB图像处理技术获取休止角Fig.11 Angle of repose is obtained using MATLAB

图12显示了在固结状态和剪切破坏状态下堆积体的休止角分布图，与文献[18]中圆形颗粒点沉积下得到的休止角12°～15°相似。文献[14,19]中探究了不同半径圆筒提升法对休止角的影响，圆筒直径越小，休止角越小，因此在考虑剪切前尺寸条件下，可以得出剪切破坏状态下的堆积体休止角均大于固结状态下的堆积体休止角，在图7、图8中堆积前剪切破坏状态初始各向异性系数大于固结状态，因此可以说初始各向异性程度越大，堆积体休止角越大。此外，固结状态下堆积体的休止角均随围压增大并略有升高，而剪切破坏状态下堆积体的休止角随着围压增大呈现不稳定的波动趋势。

图12 不同状态下堆积体的休止角分布Fig.12 Repose angle of sand pile under various states

为了分析基底压力分布，将基底分为1.4 mm宽的等距圆环，堆积体基底压力由每个区域内所含有的接触点z方向接触合力除以相应的圆环区域面积得到。图13显示了堆积体稳定后的基底竖向压力分布。

图13 基底压力分布Fig.13 Distribution of base pressure of sand pile

图13中数据经过归一化处理[6]，其中r为距离中心点水平距离，R为堆积体半径，ρwgHp为静水压力，ρw为水的密度，g为重力加速度，Hp为堆积体的高度。图13(a)结果表明固结状态下堆积体的应力峰值均出现在堆积体中心，与文献[8]中边界去除法得到的竖向压力分布趋势一致。这是由于等压固结试样的颗粒各向异性程度较低，堆积过程中顶点下方的颗粒产生的位移较小，靠近基底中心区域的力传递更接近垂直，因此产生了较强的竖向压力。而剪切破坏后的堆积体出现了压力下降现象(图12b)，与文献[8]点沉积法得到的竖向压力分布相似。这是由于剪切破坏后试样产生了具有一定强度的各向异性，且试样内部生成部分应力拱，导致中心顶点下方力的传递出现了一定倾角。模拟数据表明，初始围压固结后的堆积体对基低压力分布影响可以忽略不计，具有一定初始各向异性的颗粒材料采用边界去除法形成的堆积体会产生压力下降现象。

2.2.2力链、力链强度及力链方向

由图14(c)中的堆积体力链分布图中可以明显看到接触力力链，并且强接触力链少于弱接触力链，但强接触力链位于堆积体中心，承担着堆积体内部大部分力的传递。由于堆积体为三维模型，内部的力链难以从外部获取，因此需要将力链进行量化分析。根据OTTO et al[20]研究，采用力链形成的3个条件进行判别：1) 力链上的颗粒数要大于等于3个；2) 力链上颗粒之间的接触均为强接触；3) 力链上相邻接触之间的夹角要小于角度阈值θc，此阈值由模型内部平均配位数(〈Z〉)决定：

图14 力链角度识别Fig.14 Force chain Angle recognition

(19)

根据上述学者提出的3条力链判定条件，自编力链识别程序对堆积体内部接触进行识别，利用颗粒间强接触进行颗粒顺承接触搜索，将满足条件的接触记录下来形成完整的力链。

通过力链识别对颗粒堆积体稳定状态的力链分布进行分析，参考文献[21]中的力链概率公式进行拟合。力链长度分布及力链强度的统计规律如图14所示。

图15 堆积体的力链概率分布Fig.15 Probability distribution of force chains in sand pile

为了探究力链方向的分布规律，将力链方向投影在x-y平面进行角度分区，并获取每个区间力链数量及力链平均强度，得到的力链方向分布如图16所示。

图16 力链方向分布Fig.16 Distribution of force chain directions

图16中黑色虚线为每个区间的力链平均强度，采用归一化处理。红线为拟合曲线，采用ROTHEN-BURG et al[23]提出的公式：

(20)

其中，ωn为傅立叶系数，其值表征力链强度各向异性系数，T为最小正周期，θn为力链分布主方向。

由图16可知力链的方向分布具有一定的规律性，力链在竖直方向强度较高，力链分布主方向，θn均在90°左右。随着试样的初始围压升高，力链方向各向异性程度，ωn增大。这是因为围压大的情况下，颗粒间接触更加紧密，重力堆积时内部颗粒移动较小，其应力状态变化较小，而外部颗粒失去边界条件限制，产生较大位移，导致其堆积体内外部各向异性差异较大。同时发现固结试样的堆积体力链各向异性系数大于剪切试样的堆积体，由于应力拱分担了颗粒间力的传递，抵消了部分颗粒的重力，颗粒顺着应力拱结构分散到了堆积体外部，因此剪切破坏后形成的颗粒堆积体力链各向异性程度较小。

2.2.3堆积体的各向异性分析

为探究颗粒堆积体宏观表现与颗粒间细观组构各向异性之间的联系，从接触法向、枝向量、法向接触力及切向接触力4个方面对堆积体颗粒之间的各向异性进行分析。为探究颗粒堆积体宏观表现与颗粒间微观组构各向异性之间的联系，从接触法向、枝向量、法向接触力及切向接触力4个方面对堆积体颗粒之间的各向异性进行分析。图17为不同试样堆积体的各向异性系数，由于枝向量各向异性系数αd较小，均<0.01，不在图中显示。堆积体的综合各向异性系数α由式(18)确定，其值在稳定状态后随围压略有升高；堆积体的α≈0.7，高于等压固结后的α(≈0)，低于三轴剪切后的α(≈1.45).

图17 堆积体的各向异性系数Fig.17 Anisotropy of sand pile

图17可以看到堆积体的组构各向异性在固结状态和剪切状态存在明显差异，且剪切状态组构各向异性大于固结状态，文献[8]提到存在压力下降现象的砂堆中心区域出现了较大程度的组构各向异性，可以说组构各向异性导致其堆积体的基底中心压力下降。对比图17中各向异性系数可知，堆积体颗粒接触法向各向异性系数αc在三者中最小，并且随着围压的增大呈现出无规律波动，但其波动范围有限，因此接触法向各向异性不是堆积体宏观特征差异的主要原因。法向接触力各向异性系数αn占比最大，同时可以看到剪切试样的法向接触力各向异性程度大于固结试样的堆积体。剪切试样的法向接触力各向异性系数随围压增大而增大，说明堆积前的各向异性程度越大会导致堆积后法向接触力各向异性程度增大，这是导致固结试样堆积体宏观特征不同于剪切试样堆积体的主要原因。切向接触力各向异性系数αt大于接触法向，但固结试样与剪切试样的堆积体切向接触力各向异性大致相同，说明堆积前的各向异性不影响堆积后的切向接触力各向异性。

3 结论

本文基于三维离散元方法模拟了不同围压下的三轴排水剪切试验，利用固结、剪切后具有不同初始各向异性的试样采用边界去除法进行堆积试验，并对其宏观表现与细观组构进行分析，得到主要结论如下：

1) 本文研究所用的三轴排水试验模拟方法可以较好地反映不同围压下颗粒材料的应力应变行为、剪切变形及临界状态特性，且获得的临界状态摩擦角与真实材料试验结果一致；

2) 三轴试验的细观模拟结果揭示了试验的组构演化规律，并验证了应力-力-组构相关关系。颗粒材料具有的剪切强度与其发展的各向异性程度成正比，且力学各向异性在抗剪特性方面起主导作用；

3) 初始各向异性程度越大，堆积体的休止角越大，且初始各向异性程度大的试样内部的应力拱存在于形成的堆积体内，导致其出现了基底压力下降现象；

4) 堆积体力链分析证实了应力拱的存在，并且剪切试样的堆积体在应力拱作用下其力链强度更大，但其力链各向异性程度较小。颗粒堆积体的法向接触力各向异性是导致堆积体宏观特征差异的主要原因。

由于颗粒接触复杂程度高导致颗粒系统细微观研究较困难，但其细微观层面的机理分析对理解宏观力学行为至关重要。若想建立完整的颗粒体系研究，需要数值模拟与物理试验互为参照，进行颗粒体系多尺度研究。接下来研究工作将针对颗粒堆积动态过程的力链研究和力链稳定性分析，建立宏细观之间的理论联系，为对解决工程滑坡、沉降风险问题起到一定的促进作用。