, 陈康为, 张磊同济大学 道路与交通工程教育部重点实验室, 上海 201804
SIMULATION AND MICRO-MECHANICS ANALYSIS OF INHERENT ANISOTROPY OF GRANULAR BY DISTINCT ELEMENT METHOD1)
QianJinsong, ChenKangwei, ZhangLei中图分类号:TU441
文献标识码:A
通讯作者:
收稿日期:2018-03-12
网络出版日期:2018-09-18
版权声明:2018力学学报期刊社力学学报期刊社 所有
基金资助:
展开
摘要
关键词：
Abstract
Keywords：
-->0
PDF (15469KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNoteRisBibtex收藏本文-->
引 言
沥青路面粒料基层在施工摊铺过程中, 颗粒长轴往往会发生定向的空间排列, 形成接触分布的定向性, 进而显著影响其力学性质, 该现象称为粒料的固有各向异性[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所示.Table 1
表 1
表 1不同粒径颗粒数目
Table 1Number of different particle sizes
| Particle size/ | 4.75 ? | 6.6? | 8.0? | 9.5? | 13.25 ? |
|---|---|---|---|---|---|
| mm | 6.60 | 8.0 | 9.5 | 13.25 | 15.75 |
| particle number | 3925 | 1346 | 878 | 531 | 204 |
新窗口打开
显示原图|下载原图ZIP|生成PPT图 1颗粒级配曲线
-->Fig.1The particle grading curve
-->
1.2 颗粒形状
本文根据粒料的实际形状特征, 通过三维建模, 构建了10种不同形状的颗粒, 并以此作为建模模板, 如图2 所示.
显示原图|下载原图ZIP|生成PPT图 2不同形状粒料三维模型
-->Fig.23D granular model with different shapes
-->
将三维模型以stl格式导入到离散元软件中, 基于Taghavi[24]提出的气泡堆积算法, 利用相互重叠而不产生力的不同大小的球体来填充三维模型面, 生成复杂形状的团颗粒 (clump). 填充过程中, 采用距离控制球体间的重叠量, 采用比例控制最小球体和最大球体的比例. 距离越大比例越小, 生成的模型越接近真实形状, 但计算效率越低. 为保证计算精度及效率, 经试算比较后, 选取距离参数为100, 比例为0.3, 生成的模型如图3 所示. 相比于球形颗粒模型, 该模型更能体现粒料颗粒形状复杂、棱角多等特点, 从而考虑到粒间咬合、嵌挤等相互作用.
显示原图|下载原图ZIP|生成PPT图 3三维粒料模型的构建
-->Fig.3The construction of 3D granular model
-->
1.3 力学参数
本文选取线性接触模型作为颗粒体之间的接触法则. 同时, 参考国内外相似研究的参数取值情况[25-26], 选取相关力学参数, 如表2所示. 其中, $\kappa^*$表示法向和切向的刚度比, $E^*$ 为有效模量, 两者决定了单元间接触力与相对位移的关系; $\mu$为摩擦系数, 决定了单元间法向力和切向力的关系; 阻尼系数(damping constant) 用来控制能量的耗散.Table 2
表 2
表 2数值模型参数
Table 2Parameters of numerical model
| Stiffness | Effective | Friction | Damping | Particle |
|---|---|---|---|---|
| ratio | modulus | coefficient | constant | density/ |
| K | E* /Pa | (kg • m-3) | ||
| 1.333 | 108 | 0.5 | 0.7 | 2600 |
新窗口打开
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]. 本文采用人为固定沉积角的方法, 则是对倾角分布集中度的一种简化.Table 3
表 3
表 3试件基本信息汇总表
Table 3Basic information of specimens
| 0/(。） | Void ratio | 0/(。） | Void ratio |
|---|---|---|---|
| iso | 0.284 | 45 | 0.276 |
| 0 | 0.280 | 60 | 0.284 |
| 30 | 0.280 | 90 | 0.281 |
新窗口打开
显示原图|下载原图ZIP|生成PPT图 4固有各向异性试件
-->Fig.4Inherent anisotropic specimen
-->
1.5 试件加载
在试样制备完成后, 保持围压恒定为100 kPa, 开始对试件进行三轴压缩试验. Andrade等[29]指出, 随着加载速度的增大, 粒料的力学性能会由于动力效应而发生显著改变. 选取5种加载速度, 对各向同性试件进行了敏感性分析. 如图5所示, 当加载速度越大时, 孔隙率在加载初期减小越多, 随着加载速度减小, 孔隙率变化逐渐趋于定值, 说明加载速度存在一个准静态与动态之间的阈值. 综合考虑试验的准确性及程序的运行效率, 保证整个试件在准静态模式下进行加载, 采用0.12 m/s作为本文三轴数值试验的加载速度.
显示原图|下载原图ZIP|生成PPT图 5加载速度敏感性分析
-->Fig.5The 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所示.
显示原图|下载原图ZIP|生成PPT图 6应力比与轴向应变关系
-->Fig.6Relationship between stress ratio and axial strain
-->
显示原图|下载原图ZIP|生成PPT图 7各工况峰值应力比统计图
-->Fig.7Peak 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％.
显示原图|下载原图ZIP|生成PPT图 8体积应变与轴向应变关系图
-->Fig.8Relationship 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的围压, 因此开始进行压缩试验前, 试件内部已出现一定的初始滑动百分比. 并且, 受长轴方向与重力方向差异性的影响, 初始滑移百分比亦体现出一定的差异性; 沉积角越小的试件相对稳定, 表现出的滑动倾向越小.
显示原图|下载原图ZIP|生成PPT图 9试件接触滑动百分比统计图
-->Fig.9Contact 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所示.
显示原图|下载原图ZIP|生成PPT图 10接触法向各向异性变化图
-->Fig.10Contact normal anisotropy of specimens
-->
同时, 在离散元中, 每个接触都是以向量的形式表征的, 可通过提取的接触向量值, 计算该接触在空间中所处的位置. 根据对称原理, 取半球空间进行分析, 沿经度及纬度方向分别划分为12份, 每$15^{\circ}$一份, 统计落在交叉区域内的接触数目, 绘制接触分布三维玫瑰花图, 从而可以直观分析试件接触及接触力分布的变化特点. 统计试件初始、体积压缩最大点、强度最大点和最终状态的接触分布, 绘制出如表4所汇总的三维玫瑰花图. 其中, 玫瑰花图的三维坐标大小等于该方向的接触数目与接触平均值的比值.
Table 4
表 4
表 4接触分布三维玫瑰图汇总表
Table 4Three-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 所示. 其中, 坐标轴代表实际应力与平均应力的比值.
显示原图|下载原图ZIP|生成PPT图 11法向接触力三维玫瑰图
-->Fig.11Three-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] | . |
| [2] | . |
| [3] | . |
| [4] | . |
| [5] | . . |
| [6] | . |
| [7] | . |
| [8] | . |
| [9] | . |
| [10] | . . |
| [11] | . . |
| [12] | . |
| [13] | . |
| [14] | . |
| [15] | . |
| [16] | . |
| [17] | . . |
| [18] | . |
| [19] | . . |
| [20] | . |
| [21] | . |
| [22] | . |
| [23] | . |
| [24] | // |
| [25] | . |
| [26] | . |
| [27] | . . |
| [28] | . . |
| [29] | . |
| [30] | . |
| [31] | . |
| [32] | . |
| [33] | . |
| [34] | . |
| [35] | . |
| [36] | . |
| [37] | . |
