力学学报  2018 , 50 (5): 1041-1050 https://doi.org/10.6052/0459-1879-18-061

固体力学

粒料固有各向异性的离散元模拟与细观分析1)

钱劲松2), 陈康为, 张磊

同济大学 道路与交通工程教育部重点实验室, 上海 201804

SIMULATION AND MICRO-MECHANICS ANALYSIS OF INHERENT ANISOTROPY OF GRANULAR BY DISTINCT ELEMENT METHOD1)

Qian Jinsong2), Chen Kangwei, Zhang Lei

Key Laboratory of Road and Traffic Engineering of the Ministry of Education, Tongji University, Shanghai 201804, China

中图分类号:  TU441

文献标识码:  A

通讯作者:  2) 钱劲松, 副教授, 主要研究方向: 交通岩土工程. E-mail: qianjs@tongji.edu.cn

收稿日期: 2018-03-12

网络出版日期:  2018-09-18

版权声明:  2018 力学学报期刊社 力学学报期刊社 所有

基金资助:  1) 中央高校基本科研业务费专项资金资助项目(22120170129).

展开

摘要

料在摊铺后形成的颗粒定向排列将导致其力学性质的固有各向异性. 依据粒料的实际不规则形状, 构建了可模拟粒间咬合嵌挤作用的三维离散元复杂形状颗粒; 生成了5 种不同沉积方向的各向异性试件和1种各向同性试件, 对比了各试件在三轴压缩试验中的宏观力学特性差异; 引入组构张量以量化各向异性程度, 利用玫瑰图表达接触法向与接触力的分布特征, 探究了粒料各向异性的细观发展规律. 结果表明: 颗粒长轴愈趋向水平排布, 峰值应力比愈大, 剪缩与剪胀程度愈明显; 相较于各向同性试件, 沉积角$\theta$为$0^\circ$时, 峰值应力比和最大体积压缩应变分别提高了12.6\%和18.8\%, 其原因在于加载过程中颗粒旋转和滑动百分比更小, 内部调整时间更短、更易被剪密; 固有各向异性对颗粒法向接触力分布的影响不大, 但显著影响接触法向分布特征; 剪切过程中, $\theta$为$90^\circ$时的接触法向各向异性系数先快速减小后逐渐增大, 而$\theta$为$0^\circ$到$60^\circ$时则呈现出增大后稍有回落或趋于稳定的趋势, 且$\theta$ 愈小的试件各向异性系数增大愈快.

关键词: 颗粒材料 ; 离散元 ; 三轴压缩 ; 固有各向异性 ; 组构

Abstract

The particles tend to be spatially arranged in directional orientation after the paving of granular materials, and thus leading to the inherent anisotropy of mechanical property. Based on the actual irregular shape of granular materials, three-dimensional complex shape particles were modelled utilizing distinct element method to simulate the interlocking between particles. Five numerical test specimens with different bedding angles and an isotropic specimen were established respectively, and the mechanical properties of various specimens were compared during the triaxial compression simulations. Besides, the fabric tensor was introduced to quantify the anisotropy, the rose diagram was drawn to exhibit the distribution characteristics of contact normal and contact force, and then the development of anisotropy was investigated. It is shown that, as the long axis of particles change toward the horizontal direction, the stress ratio and the shear dilatancy of specimen increase continuously. Compared with isotropic structure, the peak stress ratio and the maximal volume compression strain of anisotropic structure when the bedding angle $\theta=0^\circ$ is 12.6% and 18.8% larger respectively. This is because the rotation and contact sliding ratio of particles is smaller, the internal adjustment time is shorter, and specimen can be sheared more densely. The inherent anisotropy has little effect on the distribution characteristics of contact force, but significantly affects the distribution characteristics of contact normal. When $\theta$ is $90^\circ$, the contact normal anisotropy coefficient drops quickly and then gradually increases during the shear process. Otherwise, the coefficient shows a steady or slight drop trend after an increase, and the coefficient grows faster as the $\theta$ decreases.

Keywords: granular materials ; distinct element method ; triaxial compression ; inherent anisotropy ; fabric

0

PDF (15469KB) 元数据 多维度评价 相关文章 收藏文章

本文引用格式 导出 EndNote Ris Bibtex

钱劲松, 陈康为, 张磊. 粒料固有各向异性的离散元模拟与细观分析1)[J]. 力学学报, 2018, 50(5): 1041-1050 https://doi.org/10.6052/0459-1879-18-061

Qian Jinsong, Chen Kangwei, Zhang Lei. SIMULATION AND MICRO-MECHANICS ANALYSIS OF INHERENT ANISOTROPY OF GRANULAR BY DISTINCT ELEMENT METHOD1)[J]. Acta Mechanica Sinica, 2018, 50(5): 1041-1050 https://doi.org/10.6052/0459-1879-18-061

引 言

沥青路面粒料基层在施工摊铺过程中, 颗粒长轴往往会发生定向的空间排列, 形成接触分布的定向性, 进而显著影响其力学性质, 该现象称为粒料的固有各向异性[1]. 已有研究表明, 固有各向异性将对路面长期性能产生显著影响, Masad等[2]和Luo等[3]通过比较不同计算模型发现, 考虑固有各向异性的计算结果更接近路面实际损坏情况, 忽略固有各向异性, 将低估路面永久变形, 并高估路面疲劳寿命. 国内外对固有各向异性进行了一系列的研究, 常用理论推导、室内单元试验或者数值模拟的方法进行分析与解释.

在相关理论研究方面, Cao等[4]扩展了各向同性强度准则, 引入了各向异性强度变量, 量化了沉积面与最大剪应力比面之间的距离, 提出了各向异性粒状材料的强度准则. 姚仰平等[5]假定了各向异性土内部强度的分布情况, 认为在外部应力作用下, 各向异性土的破坏机制与各向同性土具有明显差异. 上述理论方法常基于连续介质的假设, 依赖基于现象学的本构模型, 忽略了材料的散粒体特征, 难以揭示材料细观变化引起的各向异性的演化.

室内单元试验方面, Guo[6]制作不同沉积角的试样来模拟材料的固有各向异性, 通过直剪试验, 发现摩擦角与沉积面剪切面的相对位置有关. Oboudi等[7]等将不同各向异性的试件进行三轴剪切, 发现即使改变围压大小, 颗粒水平排布的试件均具有更大的初始剪切刚度与更高的强度. 但是, 不同制样方法, 试样的初始状态以及测试方法都可能使室内试验的结果出现差异[8], 特别在考虑复杂应力路径时, 试验的可靠性会降低[9].

相比之下, 采用离散元方法能较好地克服上述研究地局限性[10-11]. 本质上, 粒料宏观各向异性的力学特征受控于颗粒的细观组构及其变化[12], 离散元数值模拟方法则能较好地反映粒料散体介质特性, 并可对颗粒的细观参数进行实时统计[13], 进而从细观角度揭示颗粒各向异性发展的规律[12-17].

同时, 粒料颗粒棱角分明, 形状复杂, 其形状特征对于粒料力学特性有很大的影响. 相关研究显示[18], 形状不规则的颗粒相比于球形颗粒会导致散体介质呈现出更强的各向异性. 邵磊等[19]采用圆球体颗粒进行了三轴模拟实验, 拟合了室内试验应力应变曲线, 但由于没考虑颗粒形状影响, 导致大应变剪胀偏大. Mahmood等[20]采用了椭圆形颗粒构建了具有固有各向异性的试件, 研究了试件在双轴平面应变试验下固有各向异性对应力应变以及剪切带形成的影响, 发现非圆颗粒可以较好地模拟天然粒料所表现出来的力学特性. Hosseininia[21]使用不规则多边形与椭圆形颗粒进行的双轴压缩试验发现, 颗粒棱角性将影响试样的剪切强度以及体积变化. 相对而言, 既有分析粒料固有各向异性的研究, 多集中于形状规则颗粒或者二维试验的模拟, 而较少构建复杂形状颗粒并开展三维试验.

因此, 本文采用离散元作为研究手段, 通过生成复杂形状颗粒, 制备不同初始长轴朝向的试件, 模拟定围压的三轴压缩试验, 以分析试件在试验过程中的宏细观特性, 探究颗粒方向性对粒料力学特性的影响, 并分析粒料固有各向异性的来源和机理.

1 离散元模型

1.1 颗粒级配

选用DGAB25混合料级配作为数值试样的粒径分布依据, 粒径分布如图1中曲线A所示. 由于骨架颗粒中悬浮着粒径过小的颗粒, 会导致模型收敛速度大幅度降低, 故参考国内外学者采用的删减细小颗粒的方法[22], 即删去2.36 mm以下的颗粒, 形成的粒径分布如图1中曲线B所示. 进一步, 为避免尺寸效应, 相关学者[23]建议颗粒最大粒径不大于容器的1/10 (本试验容器为边长15 cm的正方体), 故采用相似级配法将B曲线的粒径再缩小一倍, 保证级配曲率系数和不均匀系数不变, 得到粒径分布如图1中曲线C 所示. 最后, 为提高计算效率, 基于等体积原则, 将级配C中粒径范围为$1.18\sim4.75~\text{mm}$颗粒的体积均摊至各档集料中, 最终计算得到不同粒径范围颗粒数目如表1所示.

表 1   不同粒径颗粒数目

Table 1   Number of different particle sizes

Particle size/4.75 6.68.09.513.25
mm6.608.09.513.2515.75
particle number39251346878531204

新窗口打开

图 1   颗粒级配曲线

Fig.1   The particle grading curve

1.2 颗粒形状

本文根据粒料的实际形状特征, 通过三维建模, 构建了10种不同形状的颗粒, 并以此作为建模模板, 如图2 所示.

图 2   不同形状粒料三维模型

Fig.2   3D granular model with different shapes

将三维模型以stl格式导入到离散元软件中, 基于Taghavi[24]提出的气泡堆积算法, 利用相互重叠而不产生力的不同大小的球体来填充三维模型面, 生成复杂形状的团颗粒 (clump). 填充过程中, 采用距离控制球体间的重叠量, 采用比例控制最小球体和最大球体的比例. 距离越大比例越小, 生成的模型越接近真实形状, 但计算效率越低. 为保证计算精度及效率, 经试算比较后, 选取距离参数为100, 比例为0.3, 生成的模型如图3 所示. 相比于球形颗粒模型, 该模型更能体现粒料颗粒形状复杂、棱角多等特点, 从而考虑到粒间咬合、嵌挤等相互作用.

图 3   三维粒料模型的构建

Fig.3   The construction of 3D granular model

1.3 力学参数

本文选取线性接触模型作为颗粒体之间的接触法则. 同时, 参考国内外相似研究的参数取值情况[25-26], 选取相关力学参数, 如表2所示. 其中, $\kappa^*$表示法向和切向的刚度比, $E^*$ 为有效模量, 两者决定了单元间接触力与相对位移的关系; $\mu$为摩擦系数, 决定了单元间法向力和切向力的关系; 阻尼系数(damping constant) 用来控制能量的耗散.

表 2   数值模型参数

Table 2   Parameters of numerical model

StiffnessEffectiveFrictionDampingParticle
ratiomoduluscoefficientconstantdensity/
KE* /Pa(kg • m-3)
1.3331080.50.72600

新窗口打开

1.4 试样制备

采用重力沉积法制备试样, 其具体步骤如下: (1) 生成尺寸为 $15~\text{cm}\times 15~\text{cm}\times30~\text{cm}$ 的沉积空间. (2) 生成颗粒时, 每档级配的颗粒按前述建立的10种颗粒形状模型表征, 并均匀分配. (3) 生成固有各向异性试件时, 将长轴水平摆放的颗粒绕水平面旋转$\theta$ 角后锁定颗粒旋转, 使其在重力作用下沉积, 从而人为定义沉积的颗粒定向排列, 重新建立边长为15 cm大小的墙体, 如图4所示, 其中定义$\theta$为沉积角.生成各向同性试件时, 则随机分布颗粒初始朝向. (4) 解除颗粒旋转锁定, 通过伺服程序控制100 kPa围压进行各向同性压缩, 保证初始试样的密实性, 达到所需应力状态[27]. 试验所制备试件的基本特性如表3所示, 其中iso表示各向同性试件.需要指出的是, 散体颗粒在振动压实后的颗粒定向排列较为复杂, 长轴倾角在一定范围内呈现概率分布[28]. 本文采用人为固定沉积角的方法, 则是对倾角分布集中度的一种简化.

表 3   试件基本信息汇总表

Table 3   Basic information of specimens

0/(。)Void ratio0/(。)Void ratio
iso0.284450.276
00.280600.284
300.280900.281

新窗口打开

图 4   固有各向异性试件

Fig.4   Inherent anisotropic specimen

1.5 试件加载

在试样制备完成后, 保持围压恒定为100 kPa, 开始对试件进行三轴压缩试验. Andrade等[29]指出, 随着加载速度的增大, 粒料的力学性能会由于动力效应而发生显著改变. 选取5种加载速度, 对各向同性试件进行了敏感性分析. 如图5所示, 当加载速度越大时, 孔隙率在加载初期减小越多, 随着加载速度减小, 孔隙率变化逐渐趋于定值, 说明加载速度存在一个准静态与动态之间的阈值. 综合考虑试验的准确性及程序的运行效率, 保证整个试件在准静态模式下进行加载, 采用0.12 m/s作为本文三轴数值试验的加载速度.

图 5   加载速度敏感性分析

Fig.5   The analysis of the velocity sensitivity

2 宏观力学特征模拟结果与分析

2.1 应力比‒应变关系

在三向应力状态下, 试件平均应力$p$、剪应力$q$和应力比分别可以用表示为

$$p=\frac{(\sigma_1+\sigma_2+\sigma_3)}{3}(1)$$

$$q=\sqrt{\frac{1}{2}\Big\lbrack (\sigma_1-\sigma_2)^2+(\sigma_2-\sigma_3)^2+(\sigma_1-\sigma_3)^2\Big\rbrack}(2)$$

$$\eta=\frac{q}{p}(3)$$

式中, $\sigma_i$($i=1,2,3$)分别代表三个方向的主应力.

为比较相同平均应力作用下不同试件所能承受的剪应力, 对各工况应力比随轴向应变的变化关系进行分析. 同时, 为验证数值模拟结果的合理性, 选取Lam等[30]通过室内试验得到的结果进行对比. 该试验同样为对不同沉积角的各向异性散粒体试件, 在围压为100 kPa的条件下进行的三轴单向压缩. 绘制模拟试验与室内试验不同试件的应力比与轴向应变的关系如图6所示, 各工况的峰值应力比如图7所示.

图 6   应力比与轴向应变关系

Fig.6   Relationship between stress ratio and axial strain

图 7   各工况峰值应力比统计图

Fig.7   Peak stress ratio under different conditions

可知, 模拟试验及室内试验的结果在宏观力学规律上具有一致性. 随着加载的进行, 两种试验结果均呈现出应力比先增大后缓落, 最终平稳并趋于一致的趋势. 在加载后期, 由于残余强度逐渐趋于一致, 应力比‒应变曲线均观察到一定程度的交叉现象. 由图7可以看出, 随着$\theta$增大, 峰值应力比不断减小. 当$\theta$为$0^{\circ}$ 时, 试件峰值应力比最大, $\theta$为$90^{\circ}$时峰值应力比最小. 模拟试验及室内试验中, $\theta$为$0^{\circ}$的试件相比$\theta$为$90^{\circ}$试件, 峰值应力比分别提高了14.3%以及11.4%.

图7可知, 进一步对比各向同性与各向异性试件的差异, 粒料的各向异性对其力学性能有显著的影响. 各向同性试件应力比处于较低水平, 与$\theta$呈$60^{\circ}$的试件水平相当. $\theta$为$0^{\circ}$的试件峰值应力比相较于各向同性试件提高了12.6%.

2.2 体积应变‒轴向应变关系

对方形试件的真三轴试验, 可统计墙体位移变化来实现试件变形量的计算. 参考Guo等[31]的定义, 将试件的轴向应变$\varepsilon_1$和体积应变$\varepsilon_v$定义如下

\begin{equation*}\varepsilon_1=\text{ln}\dfrac{H_0 }{H}\tag*{(4)}\end{equation*}

\begin{equation*}\varepsilon_v=\text{ln}\dfrac{V_0}{V}\tag*{(5)}\end{equation*}

式中, $H_0$和$V_0$分别代表加载前试件的高度和体积, $H$和$V$分别代表即时所得的试件高度和体积值.

具有不同颗粒朝向的试件在加载过程中体积应变与轴向应变的关系如图8所示. 可以看出, 在试验过程中, 试件的体积经历了"剪缩‒剪胀"的变化. 并且, 试件被剪缩和剪胀的程度随着$\theta$的增大而减小, $\theta$为$0^{\circ}$时, 试件应变变化最为显著, 从而会与体积变化较小试件的应变曲线出现交叉. 该现象与Mahmood等[20]和Zhao等[32]的试验规律相同. 同时可以观察到, 各向同性试件剪缩剪胀程度接近$\theta$为$45^{\circ}$的试件, 处于中间状态. 相较而言, $\theta$为$0^{\circ}$时, 试件体积压缩应变峰值比各向同性试件高18.8%.

图 8   体积应变与轴向应变关系图

Fig.8   Relationship between volumetric strain and axial strain

3 细观力学特征模拟结果与分析

3.1 接触滑动百分比

粒料在变形过程中, 往往伴随着颗粒的大量接触与滑移. 参照文献[33]中的定义, 当满足式(6)时, 认为接触发生滑动

\begin{equation*}\left|f^{\text t}\right|/\Big(\mu f^{\text n}\Big)>0.999~9\tag*{(6)}\end{equation*}

式中, $f^{\text t}$和$f^{\text n}$分别是切向接触力和法向接触力, 而$\mu$为接触摩擦系数.

接触滑动百分比是反映粒料集合内部受到载荷作用时颗粒运动和内部结构稳定程度的重要指标, 其大小等于颗粒滑动接触对数与总接触对数的比值. 当接触发生滑动的时候, 颗粒重新排列, 从而引发永久变形. 随着加载的进行, 各试件接触滑动百分比变化如图9所示. 需要指出的是, 试件制备过程中已施加了100 kPa的围压, 因此开始进行压缩试验前, 试件内部已出现一定的初始滑动百分比. 并且, 受长轴方向与重力方向差异性的影响, 初始滑移百分比亦体现出一定的差异性; 沉积角越小的试件相对稳定, 表现出的滑动倾向越小.

图 9   试件接触滑动百分比统计图

Fig.9   Contact sliding percentage of specimens

图9可知, 随着轴向应变的增大, 滑动百分比呈现出先增大至峰值、后逐渐减小直至平稳的总体变化趋势; 且$\theta$愈小, 接触滑动百分比愈小, 峰值出现愈早. 这是由于颗粒长轴趋于水平时, 试件结构更为稳定, 加载过程中颗粒转动程度较小, 内部调整时间较短, 在外力作用下更容易被剪密. 而长轴竖直的试件在加载过程中, 易发生较大规模的转动, 颗粒相互嵌挤, 削弱了试件的剪缩剪胀, 从而表现出图8中呈现的体积应变变化规律.

进一步对比图7可观察到, 滑动百分比越小的试件, 因在外力作用下更易密实, 在宏观力学特征上表现出更高的峰值应力比.

3.2 接触法向各向异性及空间分布

Oda等[34]指出, 非球形颗粒的空间排列和相对孔隙比具有统计意义上的相关性, 并称颗粒的这种特性为"材料的组构". 通过组构张量 (fabric tensor) 可以表征颗粒集合体细观组构的宏观响应[35].

当粒料受力后宏观行为表现很大程度决定于粒料接触力和接触法向的空间分布, 而粒料接触力和接触法向的空间分布可以用组构张量进行定量描述[36], 例如, 接触法向可以用二阶张量$R_{ij}$ 表示

式中, $n_i$是接触点单位法向量的分量, $N$是总接触数目, $\varOmega$表征接触法向量$ n$相对于全局坐标系的方向, $E(\varOmega)$是三维空间内球体分布的概率密度方程, 可以用二阶傅里叶展式表示[37]

\begin{equation*}E(\varOmega)=\dfrac{1}{4\pi}\Big\lbrack 1+ {\alpha}^r_{ij}n_in_j\Big\rbrack\tag*{(8)}\end{equation*}

式中, $\alpha^r_{ii}$为表征组构各向异性的二阶各向异性张量, 反映了接触法向分布的密度, $\alpha^r_{ii}$的主值$\alpha^r_{11}$, $\alpha^r_{22}$, $\alpha^r_{33}$分别代表了3个主应力方向接触法向分布的聚集程度. $\alpha^r_{ii}$可以由组构张量的偏量部分($R^{'}_{ij}$)表示

\begin{equation*}\alpha^r_{ij}=\dfrac{15}{2}R'_{ij}\tag*{(9)}\end{equation*}

对于张量$R_{ij}$来说, 其偏量部分可以由下式计算得到

\begin{equation*}R^{'}_{ij}=R_{ij}-\dfrac{1}{3}\varTheta\delta_{ij}\tag*{(10)}\end{equation*}

式中, $\varTheta=\sigma_{ii}$表示张量的第一不变量, $\varTheta\delta_{ij}$/3为球形张量.

对于组构的各向异性程度可以用上述张量的第二不变量进行表征

\begin{equation*}\alpha_r=\sqrt{\dfrac{3}{2}\alpha^r_{ij}\alpha^r_{ij}}\tag*{(11)}\end{equation*}

式中, $\alpha_r$为接触法向各向异性系数, 其大小反应了接触法向分布的各向异性程度. $\alpha_r$ 随加载过程的变化曲线如图10所示.

图 10   接触法向各向异性变化图

Fig.10   Contact normal anisotropy of specimens

同时, 在离散元中, 每个接触都是以向量的形式表征的, 可通过提取的接触向量值, 计算该接触在空间中所处的位置. 根据对称原理, 取半球空间进行分析, 沿经度及纬度方向分别划分为12份, 每$15^{\circ}$一份, 统计落在交叉区域内的接触数目, 绘制接触分布三维玫瑰花图, 从而可以直观分析试件接触及接触力分布的变化特点. 统计试件初始、体积压缩最大点、强度最大点和最终状态的接触分布, 绘制出如表4所汇总的三维玫瑰花图. 其中, 玫瑰花图的三维坐标大小等于该方向的接触数目与接触平均值的比值.

表 4   接触分布三维玫瑰图汇总表

Table 4   Three-dimensional rose graph of contact distribution

新窗口打开

图10可知, 不同固有各向异性试件$\alpha_r$初始值大小基本相等, 表明不同长轴方向的试件初始接触分布方向单一, 各向异性程度相近. 仅由于制样过程中少数颗粒接触方向发生改变, 造成$\alpha_r$ 初始值略有差异. 从表4中的初始接触分布玫瑰图也可以看出, 各向异性试件的初始接触分布与长轴分布方向具有良好的互余性. 以$\theta$为$0^{\circ}$为例, 颗粒长轴方向水平, 此时初始接触分布多集中在竖直方向; 而$\theta$为$90^{\circ}$ 时, 颗粒长轴方向竖直, 初始接触分布则多集中于水平方向. 各向同性试件的$\alpha_r$初始值则为0, 初始接触分布玫瑰图呈现圆球状.

对于$\theta$为$0^{\circ}$到$60^{\circ}$的试件, 在剪切过程中, $\alpha_r$ 增大至峰值后略有回落或趋于稳定. 这主要是因为随着加载的进行, 颗粒朝水平方向旋转, 水平向接触逐渐减少, 竖向接触逐渐增加, 接触分布更为集中, 表现出更强的接触法向各向异性. 待试件发生强度破坏后, 竖向应力不再增大, 竖向的接触增加接近饱和, 后期接触法向各向异性略有所减小.

对于$\theta$为$90^{\circ}$的试件, 剪切过程中$\alpha_r$呈现出先减小后逐渐增大的趋势. 这是由于颗粒初始长轴主要沿垂直分布, 在加载初期, 颗粒大量朝着水平方向旋转, 导致接触法向各向异性程度减少, 玫瑰图也从扁平状逐渐转化成趋于各向同性分布的圆球状; 而随着竖向接触的不断增加, 各向异性程度随之回升, 玫瑰图又演化成椭球状.

并且, 由图10可明显看出, $\theta$愈小的试件, $\alpha_r$增长愈快. 这是由于随着$\theta$的减小, 颗粒长轴更趋向于水平, 剪切过程中的旋转和滑动相对减小, 竖向接触增加程度更为明显, 故各向异性愈强, 同时也表现出图9中更小的接触滑动百分比, 以及图7中更大的峰值应力比.

对于各向同性的试件, $\alpha_r$增加至峰值后保持稳定, 最终与$\theta$为$60^{\circ}$的试件接近. 这是由于初始接触法向分布是各向同性的, 但随着加载的进行, 竖向接触不断增加, 玫瑰图逐渐呈现出椭球状, 各向异性程度也不断增大, 最终随着竖向接触增大至饱和, 各向异性程度也逐渐趋于稳定.

3.3 法向接触力空间分布

选取不同初始颗粒倾角的试件, 遍历各试件在强度峰值时每一个接触的方向和接触力大小 (试件其他状态时接触力空间分布状态与强度峰值时相似, 故不赘述), 绘制该时刻试件法向接触力的空间分布三维玫瑰图, 如图11 所示. 其中, 坐标轴代表实际应力与平均应力的比值.

图 11   法向接触力三维玫瑰图

Fig.11   Three-dimensional rose graph of normal contact force

各试件在峰值强度时, 法向接触力分布的三维玫瑰图形状十分相似, 均呈现竖直向 ($z$方向) 法向接触力密集, 而水平向($x$, $y$方向)法向接触力分布稀松的现象. 在加载过程中, 法向接触力分布的优势方向也始终与加载方向平行. 可知, 颗粒的定向排列对法向接触力分布的影响有限, 主要通过影响接触法向分布致使粒料力学性能产生差异.

4 结论

通过制备具有复杂形状颗粒的固有各向异性试件, 进行三轴压缩试验的离散元模拟, 可以得到以下结论.

(1) 粒料的固有各向异性将影响其宏观力学特征. 随着颗粒长轴趋向水平排布, 结构将表现出更高的抗剪强度以及更明显的剪缩剪胀特征, 其原因在于内部结构调整时间较短, 接触滑动的比例较小, 呈现出的有序状态更容易被外力压密.

(2) 不同颗粒长轴朝向的试件, 接触法向分布具有明显的方向性, 且随着偏应力的改变而不断调整;而接触力分布并不随着颗粒长轴朝向变化而改变, 优势方向也始终与加载方向平行.

(3) 沉积角$\theta$对加载过程中接触法向各向异性系数的发展具有显著影响. $\theta$为$90^{\circ}$ 时呈现出先减小后逐渐增大的规律, 而$\theta$为$0^{\circ}$到$60^{\circ}$时呈现增加至峰值后略有回落或趋于稳定的规律.

The authors have declared that no competing interests exist.


参考文献

[1] Hosseininia ES.

Investigating the micromechanical evolutions within inherently anisotropic granular materials using discrete element method

. Granular Matter, 2012, 14(4): 483-503

[本文引用: 1]     

[2] Masad S, Little D, Masad E.

Analysis of flexible pavement response and performance using isotropic and anisotropic material properties

. Journal of Transportation Engineering, 2006, 132(4): 342-349

[本文引用: 1]     

[3] Luo X, Gu F, Zhang YQ, et al.

Mechanistic-Empirical Models for Better Consideration of Subgrade and Unbound Layers Influence on Pavement Performance

. Transportation Geotechnics, 2017, 13: 52-68

[本文引用: 1]     

[4] Cao W, Wang R, Zhang JM.

Formulation of anisotropic strength criteria for cohesionless granular materials

. International Journal of Geomechanics, 2017, 17(7): 04016151

[本文引用: 1]     

[5] 姚仰平, 祝恩阳.

横观各向同性土的简明破坏机制解释

. 岩土力学, 2014(2): 328-333

[本文引用: 1]     

(Yao Yangping, Zhu Enyang.

Concise interpretation of damage mechanism for cross-anisotropic soil

. Rock and Soil Mechanics, 2014(2): 328-333 (in Chinese))

[本文引用: 1]     

[6] Guo P.

Modified direct shear test for anisotropic strength of sand

. Journal of Geotechnical & Geoenvironmental Engineering, 2008, 134(9): 1311-1318

[本文引用: 1]     

[7] Oboudi M, Pietruszczak S, Razaqpur AG.

Description of inherent and induced anisotropy in granular media with particles of high sphericity

. International Journal of Geomechanics, 2016, 16(4): 04016006

[本文引用: 1]     

[8] Kawajiri S, Kawaguchi T, Yamasaki S, et al.

Strength Characteristics of Compaction Soil with Particular Reference to Soil Structure And Anisotropy

. International Journal of Geomate, 2017, 13(38): 178-185

[本文引用: 1]     

[9] Sazzad MM, Suzuki K, Modaressi-Farahmand-Razavi A.

Macro-micro responses of granular materials under different b values using DEM

. International Journal of Geomechanics, 2013, 12(3): 220-228

[本文引用: 1]     

[10] 章青, 顾鑫郁, 杨天.

冲击载荷作用下颗粒材料动态力学响应的近场动力学模拟

. 力学学报, 2016, 48(1): 56-63

[本文引用: 1]     

(Zhang Qing, Gu Xinyu, Yang Tian.

Peridynamics simulation for dynamic response of granular material under impact loading

. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 56-63 (in Chinese))

[本文引用: 1]     

[11] 费明龙,徐小蓉,孙其诚.

颗粒介质固‒流态转变的理论分析及实验研究

. 力学学报, 2016, 48(1): 48-55

[本文引用: 1]     

(Fei Minglong, Xu Xiaorong, Sun Qicheng.

Studies on the transition between soil-and fluid-like states of granular materials

. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 56-63 (in Chinese))

[本文引用: 1]     

[12] Jiang MJ, Li T, Shen ZF.

Fabric rates of elliptical particle assembly in monotonic and cyclic simple shear tests: A numerical study

. Granular Matter, 2016, 18(3): 1-14

[本文引用: 2]     

[13] Yang H, Xu WJ, Sun QC, et al.

Study on the meso-structure development in direct shear tests of a granular material

. Powder Technology, 2017, 314: 129-139

[本文引用: 1]     

[14] Dai BB, Yang J, Zhou CY, et al.

DEM investigation on the effect of sample preparation on the shear behavior of granular soil

. Particuology, 2016, 25: 111-121

[15] Lu XL, Zeng S, Qian JG, et al.

DEM analysis of the shear strength of cross-anisotropic sand with non-spherical particles

. Géotechnique Letters, 2017, 7(3): 230-236

[16] Jiang MJ, Sima J, Li LQ, et al.

Investigation of influence of particle characteristics on the non-coaxiality of anisotropic granular materials using DEM

. International Journal for Numerical and Analytical Methods in Geomechanics, 2017, 41(2): 198-222

[17] 张坤勇, 李威, 罗兴军.

基于PFC2D的砂土原生各向异性微观机理数值试验

. 岩土工程学报, 2017, 39(3): 518-524

[本文引用: 1]     

(Zhang Kunyong, Li Wei, Luo Xingjun, et al.

Numerical experiments of microscopic mechanism of inherent anisotropy for sand based on PFC2D

. Chinese Journal of Geotechnical Engineering, 2017, 39(3): 518-524 (in Chinese))

[本文引用: 1]     

[18] Cai Y, Yu HS, Li X, et al.

Noncoaxial behavior of sand under various stress paths

. Journal of Geotechnical & Geoenvironmental Engineering, 2013, 139(8): 1381-1395

[本文引用: 1]     

[19] 邵磊, 迟世春, 贾宇峰.

堆石料大三轴试验的细观模拟

. 岩土力学, 2009,30(S1): 239-243

[本文引用: 1]     

(Shao Lei, Chi Shichun, Jia Yufeng.

Meso-mechanical simulation of a large scale triaxial text of rockfill materials

. Rock and Soil Mechanics, 2009, 30(S1): 239-443 (in Chinese))

[本文引用: 1]     

[20] Mahmood Z, Iwashita K.

Influence of inherent anisotropy on mechanical behavior of granular materials based on DEM simulations

. International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34(8): 795-819

[本文引用: 2]     

[21] Hosseininia ES.

Discrete element modeling of inherently anisotropic granular assemblies with polygonal particles

. Particuology, 2012, 10(5): 542-552

[本文引用: 1]     

[22] Gu XQ, Huang MS, Qian JG.

DEM investigation on the evolution of microstructure in granular soils under shearing

. Granular Matter, 2014, 16(1): 91-106

[本文引用: 1]     

[23] Kim BS, Park SW, Kato S.

DEM simulation of collapse behaviours of unsaturated granular materials under general stress states

. Computers and Geotechnics, 2012, 42: 52-61

[本文引用: 1]     

[24] Taghavi R.

Automatic clump generation based on mid-surface

//11 Proc. 2nd International FLAC/DEM Symposium, Melbourne, 2011: 791-797

[本文引用: 1]     

[25] Goldenberg C, Goldhirsch I.

Friction enhances elasticity in granular solids

. Nature, 2005, 435(7039): 188-191

[本文引用: 1]     

[26] Minh NH, Cheng YP.

A DEM investigation of the effect of particle-size distribution on one-dimensional compression

. Geotechnique, 2013, 63(1): 44-53

[本文引用: 1]     

[27] 吴越, 杨仲轩, 徐长节.

初始组构各向异性对砂土力学特性及临界状态的影响

. 岩土力学, 2016, 37(9): 2569-2576

[本文引用: 1]     

(Wu Yue, Yang Zhongxuan, Xu Changjie.

Effects of initial fabric anisotropy on mechanical behavior and critical state of granular soil

. Rock and Soil Mechanics, 2016, 37(9): 2569-2576 (in Chinese))

[本文引用: 1]     

[28] 蔡氧, 张磊, 钱劲松.

细砂路基碾压密实的颗粒流模拟与试验验证

.同济大学学报(自然科学版), 2017, 45(4): 527-532

[本文引用: 1]     

(Cai Yang, Zhang Lei, Qian Jinsong.

Discrete element method simulation and experimental verification on roller compaction of fine sand

. Journal of Tongji University (Natural Science). 2017, 45(4): 527-532 (in Chinese))

[本文引用: 1]     

[29] Andrade JE, Chen Q, Le PH, et al.

On the rheology of dilative granular media: Bridging solid-and fluid-like behavior

. Journal of the Mechanics and Physics of Solids, 2012, 60(6): 1122-1136

[本文引用: 1]     

[30] Lam WK, Tatsuoka F.

Effects of initial anisotropic fabric and $\sigma_2$ on strength and deformation characteristics

. Soils and Foundations, 1988, 28(1): 89-106

[本文引用: 1]     

[31] Guo N, Zhao J.

The signature of shear-induced anisotropy in granular media

. Computers and Geotechnics, 2013, 47(1): 1-15

[本文引用: 1]     

[32] Zhao J, Guo N.

The interplay between anisotropy and strain localisation in granular soils: A multiscale insight

. Géotechnique, 2015, 65(8): 642-656

[本文引用: 1]     

[33] Gu X, Huang M, Qian J.

DEM investigation on the evolution of microstructure in granular soils under shearing

. Granular Matter, 2014, 16(1): 91-106

[本文引用: 1]     

[34] Oda M, Koishikawa I, Higuchi T.

Experimental study of anisotropic shear strength of sand by plane strain test

. Soils and Foundations, 1978, 18(1): 25-38

[本文引用: 1]     

[35] Wang R, Fu PC, Zhang JM, et al.

Evolution of various fabric tensors for granular media toward the critical state

. Journal of Engineering Mechanics, 2017, 143(10): 04017117

[本文引用: 1]     

[36] Sitharam T, Dinesh SV, Shimizu N.

Micromechanical modelling of monotonic drained and undrained shear behaviour of granular media using three-dimensional DEM

. International Journal for Vumerical and Analytical Methods in Geomechanics, 2002, 26(12): 1167-1189

[本文引用: 1]     

[37] Rothenburg L, Bathurst RJ.

Analytical study of induced anisotropy in idealized granular materials

. Geotechnique, 1989, 39(4): 601-614

[本文引用: 1]     

/