0.   引言

岩土工程、水利工程、道路桥梁工程等的附存环境或对象相当大一部分是由离散形式的颗粒材料组成,如砂土、碎石、粗粒土等,可统称为岩土颗粒材料。岩土颗粒材料的力学特性直接影响地（路）基、堤防、堆石坝等工程的设计、施工和运行。颗粒材料看似简单,但其物理力学性质复杂,存在压硬性、剪胀性、应力路径相关性等复杂力学特性^[1-2]。这些复杂的力学行为与其离散性、多尺度和能量耗散机制有关^[1]。对颗粒材料进行宏观唯象建模不足以反映其所有力学特性^[3-4],颗粒材料的微观动力学行为、接触网络、介观结构、细观组构等为揭示颗粒材料复杂的宏细观力学特性提供了新的途径^[5-7]。

细观尺度下颗粒材料的研究手段主要包括物理实验和数值模拟两大类。物理实验方面,常用的观测手段有光弹实验、扫描电镜技术、CT（computed tomography）扫描技术、DIC（digital image correlation）成像技术等^[4]。数值分析方法有离散单元法（discrete element method,DEM）^[5-6]、连续-离散耦合分析方法（combined finite and discrete method,FDEM）^[8-9]等。其中,颗粒离散单元法（DEM）发展迅速,计算简便,已广泛应用于岩土类颗粒材料的研究中,经过不断改进能够对颗粒形状^[10]、颗粒破碎^[11]、以及颗粒材料的临界状态^[12-13]、剪胀^[13]、细观接触^[14]和组构^{[6, 14]}等特性进行表征。

岩土颗粒材料是多体相互作用的耗散体系,颗粒间摩擦是耗散外力功的主要途径之一,是颗粒材料宏观剪切强度的主要来源之一,也是其变形和破坏特性的重要内部控制因素^[4]。颗粒间摩擦特性会影响颗粒间接触力的大小分布和颗粒重排列,进而引起颗粒体系宏观力学特性的变化。Dai等^[15]采用玻璃珠进行了一系列直剪试验,发现颗粒间摩擦对颗粒材料剪切行为的影响显著。Sandeep等^[16]指出颗粒间的摩擦系数与杨氏模量和表面粗糙度成反比。Senetakis等^[17]采用破碎的石灰石颗粒进行微观力学滑动实验并对颗粒间的摩擦系数进行量化。

细观数值模拟方面,Rothenburg等^[12]采用离散单元法研究了不同摩擦系数的试样孔隙比与配位数之间的关系；Antony等^[5]分析了摩擦特性对颗粒材料宏观抗剪强度的影响；Huang等^[18]评估了颗粒材料的临界状态对颗粒间摩擦的敏感性；Barreto等^[19]认为摩擦系数影响了强接触力链的稳定性。尽管关于颗粒间摩擦系数对颗粒材料宏细观力学特性的影响研究较多,但摩擦特性对颗粒材料在三维复杂加载路径下力学行为的影响,以及不同颗粒间摩擦系数下各类细观结构对宏观抗剪强度和体积变形的作用机制,还有待进一步的研究。

颗粒间接触形成的接触力网络是一个高度非均质、强非线性的的复杂网络。颗粒体系内大部分的力是通过相对稀疏的链状网络（称之为“力链”）传递^[20]。Radjai等^[21]、Thornton等^[22]将颗粒体系的接触网络划分为两个互补的子接触网络：强接触网络和弱接触网络。大于系统平均接触力的接触形成的强接触网络是传递力的主要通道；弱接触网络中的接触力小于系统平均接触力,为强接触网络的稳定提供支撑^[22]。分别对强弱接触网络的配位数、组构等细观参量进行分析,是研究颗粒材料强度和变形的有效途径^{[5, 23-24]}。

本文将基于岩土颗粒材料真三轴离散元数值试验,分析细观滑动摩擦与颗粒体系强度、变形之间的关系,研究不同摩擦系数下强、弱接触网络的配位数分布、接触力、组构张量及各向异性的特征,探讨强、弱接触网络对颗粒材料宏观应力-应变关系的贡献,揭示细观尺度摩擦特性对宏观力学行为的作用机制。

1.   DEM数值试验

1.1   数值试验及参数

本文利用开源软件LIGGGHTS进行离散元模拟^[25]。进行真三轴加载路径的DEM模拟之前,需要制备各向同性且具有一定围压的立方体数值试样。制样时,在400 mm×400 mm×400 mm的立方体空间中随机生成粒径服从Gaussian分布的圆球颗粒,共计31253个。为了避免制样产生各向异性,在试样的各个方向采用位移控制等速地压缩试样直至目标大小。随后,施加三向等压应力进行各向同性固结,直至达到预定的围压值0.5 MPa,得到试样的初始孔隙比为0.582。采用Hertz-Mindlin非线性接触模型^[26]描述颗粒间的接触力,颗粒间滑移服从库仑摩擦定律,即颗粒间发生滑动的条件为｜F_t｜≥μ｜F_n｜,μ为颗粒间的滑动摩擦系数。加载过程中的μ值分别取0.05,0.1,0.2,0.3,0.4,0.5。其它DEM模拟参数如表1所示。

表  1  数值模拟细观参数

Table  1.  Microscopic parameters of numerical simulation

参数	值
密度ρ/(kg·m-3)	2600
颗粒粒径d/mm	2.25~8.25
弹性模量E/GPa	65.0
泊松比$\nu$	0.12
恢复系数e	0.95
滑动摩擦系数μ	0.05,0.1,0.2,0.3,0.4,0.5

下载: 导出CSV 

| 显示表格

1.2   真三轴加载路径

采用等p（围压）等b（中主应力系数）加载路径对各向同性试样进行加载。加载时,为了消除边界效应,避免应变局部化,采用周期性边界^[27]。加载过程中,轴向加载速度为-0.01 mm/s。采用周期性边界的伺服机制,使试样的应力状态达到预定的值以满足真三轴应力路径^{[13, 27]},试样中主应力${\sigma }_{\text{2}}$和小主应力${\sigma }_{\text{3}}$的控制条件^[13]分别为

颗粒材料的广义剪应力q,平均静水压力p,体积应变${\epsilon }_{\text{v}}$,广义剪应变${\epsilon }_{\text{d}}$均可由应力应变张量主值导出：

2.   细观摩擦与宏观强度

2.1   加载过程中的应力应变

不同的摩擦系数下,剪应力q,体积应变${\epsilon }_{\text{v}}$随剪应变${\epsilon }_{\text{d}}$在三轴压缩路径下（b=0）的演化曲线如图1所示。图中体积应变以剪缩为正,剪胀为负。随着颗粒间摩擦系数μ的增大,峰值剪应力逐渐增大,且峰值后的应变软化程度愈加明显。体积响应表现为：随着μ值的增加,试样更快进入剪胀状态,且试样体积膨胀程度越明显。

图  1  应力应变曲线（b=0.0）

Figure  1.  Curves of stress-strain relationship（b=0.0）

下载: 全尺寸图片 幻灯片

将不同摩擦系数下等p等b真三轴数值试验得到的峰值应力状态点绘制在π平面,如图2所示,为大小不同的三角锥形。π平面内的应力演化反映出不同b值,μ值条件下,试样剪应力峰值的大小,可以看出,试样的剪应力峰值随着b值的增大而减小,随着μ值的增大而不断增大,且μ值较大时,试样对滑动摩擦系数的敏感度减弱,体现在剪应力的增加幅度逐渐减小。

图  2  不同摩擦系数下剪应力峰值状态

Figure  2.  Peak states of deviatoric stress of samples under different friction coefficients

下载: 全尺寸图片 幻灯片

2.2   三维强度准则

对于无黏结的颗粒材料,其抗剪强度${\tau }_{\text{f}}$=${\sigma }_{\text{n}}$tan(${\phi }_{\text{mob}}$),其中${\sigma }_{\text{n}}$是破坏面上的垂直压力,${\phi }_{\text{mob}}$=sin^-1$({\sigma }_1-{\sigma }_3)$/$({\sigma }_1+{\sigma }_3)$是宏观内摩擦角,是颗粒材料的抗剪强度指标,反映了颗粒材料的宏观摩擦特性。Skinner等^[28]指出峰值状态的宏观内摩擦角${\phi }_{\text{max}}$随着颗粒间摩擦系数μ的增大而增大。图3总结了多个学者得出的${\phi }_{\text{max}}$与粒间摩擦系数μ的关系^{[19, 28-32]}。结合本文所得不同b值下的峰值内摩擦角,可以发现随着摩擦系数μ的增大,峰值内摩擦角${\phi }_{\text{max}}$的增大变化幅度反而减小,与宏观抗剪强度的变化吻合（如图2所示）。当b=0.5时,峰值内摩擦角${\phi }_{\text{max}}$较大,且${\phi }_{\text{max}}$在三轴拉伸路径（b=1.0）下的值大于^{[23, 32]}三轴压缩路径（b=0.0）,与已有的DEM和物理试验研究结果基本一致。

图  3  峰值内摩擦角φ_max与滑动摩擦系数μ的关系

Figure  3.  Relationship between peak internal friction angle φ_max and inter-particle friction coefficient μ

下载: 全尺寸图片 幻灯片

岩土材料的强度理论中,最常采用的是Mohr-Coulomb^[33]和Drucker-Prager^[34]强度准则,但Mohr-Coulomb准则没有考虑中主应力对强度的影响,Drucker-Prager准则不能反映三轴拉伸和压缩应力状态下抗剪强度的不同。于是,国内外学者^[35-38]基于岩土材料的真三轴试验结果提出了更多的三维强度准则,其中常见的有Lade-Duncan准则^[35]和Mastsuoka- Nakai准则^[36]。为了探讨不同摩擦系数下各个强度准则的预测能力,以便于选取一个预测精准且形式简单的强度准则,定义以下指标反映强度准则的预测能力：

式中,${\sigma }_{\text{predicated}}^i$为第i个中主应力系数时的预测值,${\sigma }_{\text{measured}}^i$为相应数值试验得到的应力破坏值。

将不同摩擦系数下各强度准则的预测误差列于表2,并绘制摩擦系数为0.05,0.2和0.5时各强度准则以及数值模拟所得的破坏应力在π平面的破坏迹线^[37-38],如图4所示。图中TC和TE分别表示三轴压缩（triaxial compression）和三轴拉伸（triaxial extension）应力路径。Mohr-Coulomb准则的破坏迹线为不等边三角形,比较靠近静水压力轴,是所有强度准则中最保守的。Drucker-Prager准则在π平面的破坏迹线是一个圆,由于其过高的估计了中主应力对强度的贡献,导致其预测的强度与数值模拟结果相差较大,但随着摩擦系数的减小,其预测误差逐渐降低。Mastsuoka-Nakai准则预测的破坏点在b=0和b=1时与Mohr-Coulomb准则重合,在0<b<1区间内,Mastsuoka-Nakai准则破坏迹线随着摩擦系数的减小逐渐外凸。Lade-Duncan准则在μ=0.5时与数值试验结果最为接近,误差值只有0.95%。当μ≤0.2时,摩擦特性对不同b值下的强度影响较小,各强度准则在π平面的破坏迹线逐渐向Drucker-Prager准则靠拢,此时各强度准则的预测能力相差不大。当0.2<μ≤0.4时,Mastsuoka-Nakai准则预测能力较好。当0.4<μ≤0.5时,Lade-Duncan准则的预测精度较高。由于岩土颗粒材料属于摩擦性颗粒材料,其细观摩擦系数大约在0.3～0.5之间^{[18, 27]},上述关于强度准则的对比,进一步说明了Mastsuoka- Nakai准则和Lade-Duncan准则在一般岩土材料中的适用性。

表  2  三维强度准则预测误差表

Table  2.  Prediction errors of three-dimensional strength criteria

摩擦系数	误差/%
	Drucker-Prager	Mohr-Coulomb	Lade-Duncan	Mastsuoka-Nakai
0.05	6.32	5.18	3.33	2.57
0.10	10.57	6.64	7.86	6.81
0.20	13.60	7.00	8.57	5.46
0.30	15.16	8.05	5.94	1.88
0.40	15.75	10.03	3.91	1.23
0.50	16.33	12.05	0.95	5.07

下载: 导出CSV 

| 显示表格

图  4  破坏面在π平面上的迹线

Figure  4.  Trace of failure surface on π plane

下载: 全尺寸图片 幻灯片

3.   摩擦特性与细观结构演化

3.1   强弱接触网络内的配位数分布

按照Radjai等^[21]提出的方法,用平均接触力将颗粒体系的总接触网络划分为强接触网络和弱接触网络。本文以Г_all,Г_strong,Г_weak分别表示总、强、弱接触网络。图5表明不同加载路径下峰值状态强、弱接触网络中平均法向和切向接触力比值与摩擦系数μ的关系基本不受b值的影响。由接触力比值可知,强接触网络对总接触网络平均接触力的贡献较大。随着μ值的增加,强弱接触网络平均法向接触力的比值随之增大,而平均切向接触力的比值随之减小,这是因为μ值的增加使得强接触网络中平均法向接触力明显增大,而弱接触网络中平均法向接触力增加较小。

图  5  强弱接触网络之间平均接触力的比值

Figure  5.  Stress ratios of average contact stress between strong and weak contact networks

下载: 全尺寸图片 幻灯片

峰值状态下,强、弱和总接触网络的配位数分布情况如图6所示。图中的配位数是指单个颗粒的接触数目,与之对应的颗粒数目是指接触网络中具有相同配位数的颗粒个数。Bratberg等^[41]研究了二维状态下颗粒材料内部的接触连通性,结果表明配位数越少的颗粒所受约束少,各向异性程度高；相应配位数越多的颗粒越稳定,各向异性程度低。在总接触网络中,随着μ值的增大,配位数众数由5减小为4,且配位数大于5或小于4的分布曲线均发生左移,试样整体的平均配位数减小,各向异性程度增加。弱接触网络中颗粒配位数分布与总接触网络的分布规律相似,但配位数众数由4减小为3。

图  6  峰值状态配位数的分布

Figure  6.  Distribution of coordination numbers at peak stress

下载: 全尺寸图片 幻灯片

值得注意的是,强接触网络的配位数分布与总网络差异较大,μ值的改变对强接触网络配位数为1,2的颗粒数目影响甚微,这些颗粒是构成“力链”末端及“力链”主体的颗粒群^{[23, 41]},此现象间接印证了强接触网络中“力链”的存在,而且表明在峰值状态下,试样内形成“力链”的颗粒数目不随μ值的改变而波动,或许意味着此时“力链”结构已经达到最优的颗粒分配状态。同时,配位数大于3的颗粒数目减少,说明摩擦系数越大,试样内传递力的强接触网络越稀疏,接触各向异性越明显。因而试样宏观抗剪强度增加的原因可以认为是强接触网络法向接触力的明显增大。

对比不同加载路径下强、弱和总接触网络配位数的分布曲线,发现其基本表现出与加载路径无关的特性,与刘嘉英等^[6]的研究结果一致。在既定的加载路径下,总接触网络配位数的分布中弱接触网络占主导位置,弱接触网络中配位数高的颗粒数目较多,结合配位数与颗粒材料体积变形的关系,可以推断弱接触网络对试样整体的体积变形贡献较大。

3.2   组构张量及其偏值在π平面的分布

岩土颗粒材料的宏观力学特性与加载过程中细观组构的演化规律密切相关^[42]。文章采用Satake^[43]和Oda^[44]提出的组构张量表达方式：

式中,n是沿着颗粒间接触面法线方向的单位矢量,N_c为颗粒体系中的接触总数。组构张量偏张量${{\Phi }^{\prime }}_{ij}$的第二不变量可以表征颗粒体系接触方向的各向异性,用等效各向异性系数A计算如下：

强接触网络的组构张量定义为

式中,$n_i^s$是单位法向量在第s个强接触处的第i个分量。N_s是颗粒体系中的强接触总数。强接触网络的等效各向异性系数A^s如下表示：

类比于广义剪应力q和平均应力p的定义,${\Phi }_{\text{d}}$和${\Phi }_{\text{m}}$可以表示为

而强接触网络中${\Phi }_{\text{d}}^{\text{s}}$和${\Phi }_{\text{m}}^{\text{s}}$分别定义为

图7是峰值状态下总接触网络和强接触网络组构张量各向异性不变量A和A^s在π平面的分布。总接触网络的组构张量不变量A在π 平面内呈现出倒三角锥形,与应力张量分布（如图2）的形状相反。而强接触网络组构张量不变量A^s在π平面内的分布表现出与宏观应力响应较好的一致性,说明强接触网络对宏观应力的贡献更为显著。在剪切变形中,不同b值和μ值的峰值应力的变化主要是由强接触网络的组构各向异性主导。Thornton等^{[29, 45]}指出总接触网络的组构各向异性程度小于强接触网络（对比图7（a）,7（b））。这是由颗粒体系的传力方式造成的,强接触承担颗粒体系大部分的力,会选择有利的接触方向以平衡外荷载的施加,使得强接触的传力网络呈现出强烈的方向各向异性。

图  7  π平面峰值状态组构张量偏张量不变量

Figure  7.  Anisotropy invariants of fabric tensors on $\pi$ plane for peak stress states

下载: 全尺寸图片 幻灯片

应力比q/p与强接触网络组构比${\Phi }_{\text{d}}^{\text{s}}$/${\Phi }_{\text{m}}^{\text{s}}$的关系见图8。不同μ值下,q/p与${\Phi }_{\text{d}}^{\text{s}}$/${\Phi }_{\text{m}}^{\text{s}}$的比值始终约为1。Sazzad等^[24]指出不同加载路径下,应力比与强接触网络的组构比亦约为1。说明当只考虑强接触网络时,宏观应力比q/p与细观组构比${\Phi }_{\text{d}}^{\text{s}}$/${\Phi }_{\text{m}}^{\text{s}}$之间的关系始终呈线性唯一,该特性与粒间摩擦和加载路径均无关。同时也证明了颗粒体系的宏观力学特性与强接触网络组构各向异性的密切性。结合图2和图7可知,颗粒体系宏观抗剪强度随粒间摩擦的增大可能是由细观强接触网络组构各向异性的加强引起的。

图  8  不同μ值下q/p和${\Phi }_{\text{d}}^{\text{s}}$/${\Phi }_{\text{m}}^{\text{s}}$的关系（b=0.0）

Figure  8.  Relationship between q/p and ${\Phi }_{\text{d}}^{\text{s}}$/${\Phi }_{\text{m}}^{\text{s}}$ for different values of μ（b=0.0）

下载: 全尺寸图片 幻灯片

3.3   强弱接触网络各向异性的特征及其演化

Yimsiri等^[46]在定量分析颗粒体系的各向异性时,定义了两类不同机理的各向异性：几何各向异性和力学各向异性。颗粒材料的几何各向异性可以用接触法向向量和枝向量的分布表示,此外组构张量的表示形式还可以如下所示^[43-44]：

式中,Θ表示域内n的方向分布,E(Θ)是概率分布函数。一般采用E(Θ)的傅立叶展开式描述接触法向：

其中二阶各向异性张量$a_{ij}^{\text{c}}$是一个偏张量,表征接触方向的各向异性程度。表达方式为

式中,${{\Phi }^{\prime }}_{ij}$是${\Phi }_{ij}$的偏张量。

由于本文数值模拟采用的是理想圆球颗粒,枝向量几乎与接触法向同向,且颗粒尺寸分布相对较窄,故枝向量各向异性对总强度的贡献很小。力学各向异性包括法向接触力引起的法向接触力各向异性aⁿ和切向接触力引起的切向接触力各向异性a^t,其具体定义详见文献[47]。

试样初始状态是各向同性的,外荷载的施加使得试样内部结构急剧变化,因而,初期试样各向异性系数均随剪应变的增加而迅速增大,且摩擦系数μ越大,强接触网络各向异性系数增加速率越快,这是由于颗粒间摩擦越大,越能积极有效地调动强接触网络的颗粒参与外力的抵抗。从强、弱接触网络中a^c,aⁿ的演化曲线（如图9）可知,强接触网络中a^c,aⁿ值均大于弱接触网络,且均随着摩擦系数的增加而增大,说明强接触力沿加载方向的支撑作用明显,强接触网络的各向异性对总接触网络的各向异性贡献较大。

图  9  强弱接触网络中各向异性系数的演化

Figure  9.  Evolution of anisotropic coefficient in strong-weak contact system

下载: 全尺寸图片 幻灯片

4.   结论

（1）在颗粒体系真三轴加载条件下,随着摩擦系数的增大,宏观颗粒体系的剪应力峰值增大,进入剪胀状态更迅速,且后期剪胀程度更明显。当摩擦系数μ≤0.2时,Mohr-Coulomb、Drucker-Prager、Lade- Duncan和Matsuoka-Nakai准则的预测能力相差不大；当摩擦系数为0.2<μ≤0.5时,Lade-Duncan和Matsuoka-Nakai准则能较好地拟合不同加载路径下的数据点,更适用于一般岩土颗粒材料。

（2）粒间滑动摩擦对颗粒体系中强、弱接触网络的配位数分布影响差异明显。在三维加载条件下,各接触网络的配位数分布表现出与加载路径无关的特性。在剪应力峰值状态,弱接触网络的配位数分布受μ值的影响较大,但强接触网络中参与形成“力链”的颗粒数目基本保持不变。此外,强接触网络的组构张量偏不变量在π平面内的分布表现出与宏观应力响应较好的一致性,且颗粒体系宏观应力比q/p与细观强接触网络组构比${\Phi }_{\text{d}}^{\text{s}}$/${\Phi }_{\text{m}}^{\text{s}}$的关系呈线性唯一,与粒间摩擦和加载路径均无关,说明颗粒体系中的强接触网络在应力传递过程中发挥了重要且恒定的作用。

（3）滑动摩擦系数对强、弱接触网络的细观接触力及其组构各向异性影响显著。颗粒间摩擦系数的增加显著提高了强接触网络的法向接触力以及各接触网络的各向异性系数,其中强接触网络的法向接触力各向异性对整体力学各向异性贡献显著。