1. 武汉大学 土木建筑工程学院, 湖北 武汉 430072;
2. 西南交通大学 力学与工程学院, 四川 成都 610031
收稿日期：2021-01-04
基金项目：国家自然科学基金资助项目(11772237，11472196)。
作者简介：潘宇(1998-), 男, 安徽阜阳人, 武汉大学硕士研究生;
楚锡华(1977-), 男, 河南濮阳人, 武汉大学教授, 博士生导师。
摘要：采用离散单元法, 研究线弹性接触模型、Hertz接触模型、线性平行粘结模型对不同排列颗粒材料中波传播的影响, 分析了三种接触模型下波的波速变化和能量变化, 同时也研究了弹性模量对颗粒材料中波传播的影响.最后对颗粒材料样本采用Hertz接触模型, 探究了不同围压条件下颗粒材料中波的波速变化和能量衰减.研究结果表明: 不同接触模型下波传播的波速和能量衰减率不同; 当颗粒样本从规则排列到随机排列, 波传播的波速变小, 能量衰减率变大; 随着颗粒弹性模量的增大, 波速呈对数增大; 在波传播过程中能量快速耗散, 随着围压的增大, 波速和能量衰减率都增大.
关键词：颗粒材料离散元波传播接触模型围压
Influence of Contact Behavior on Stress Wave Velocity and Energy in Granular Materials
PAN Yu1, WANG Jiao1,2, CHU Xi-hua1
1. School of Civil Engineering, Wuhan University, Wuhan 430072, China;
2. School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu 610031, China
Corresponding author: CHU Xi-hua, E-mail: chuxh@whu.edu.cn.
Abstract: The discrete element method is used to study the energy attenuation and wave velocity change of granular system in the process of wave propagation with linear-elastic contact model, Hertz contact model and linear parallel bond contact model. Meanwhile, the influences of elastic modulus on wave propagation in granular materials are also studied. Finally, the Hertz contact model is used to study the wave velocity and the energy attenuation of granular materials under different confining pressures. The results show that the wave velocity and energy attenuation rate of wave propagation are different with different contact models; when the sample is arranged from regular to random, the wave velocity is smaller and the energy attenuation rate is larger; with the increase of elastic modulus of particles, the wave velocity increases logarithmically; in the process of wave propagation, the energy dissipates rapidly, and with the increase of confining pressure, the wave velocity and energy attenuation rate both increase.
Key words: granular materialdiscrete elementwave propagationcontact modelconfining pressure
尺度在1 μm~104 m范围的物质都可称为颗粒材料.在自然界、日常生活和生产中, 颗粒材料普遍存在, 且其具有许多没有被认识清楚的奇特性质和规律.随着计算机计算能力的快速提高和现代先进测试仪器的运用, 以及人们对自然现象认识的深入, 对颗粒材料的研究已经成为物理力学的前沿热点之一[1-2].对颗粒材料的研究由于受到实验条件等诸多因素限制, 往往很难进行或者很难得到正确的实验结果, 因此目前对颗粒物质的研究大量采用计算机模拟技术, 如离散单元法[3](discrete element method, DEM).
近年来, 颗粒材料中的波传播受到了研究人员的广泛关注[4-5].研究表明: 颗粒材料中波传播方向上的组构分布[6]对波速和波幅衰减都有重大影响[7-8]; 颗粒材料可能表现出反常色散[9], 即波速随激发频率的增加而增加; 不同接触模型的模拟结果与实验数据拟合的效果不同, 非线性接触模型在忽略阻尼时可获得与实验数据拟合较好的效果[10].在各种类型如单次冲击[11]、周期性[12]或随机冲击序列等的载荷作用下激发的颗粒材料中的波传播在一些领域(如无损检测、能量收集、确定生物材料的机械特性[13-14]等)的应用中具有重要的实际意义.同时, 当颗粒材料受到外载荷作用产生波传播时, 材料内部会发生强烈的挤压和碰撞, 使系统能量迅速有效地衰减, 从而起到良好的缓冲作用[15-16], 因此颗粒材料也是吸能减震的优选材料.
本文运用离散单元法, 主要研究不同接触模型下颗粒材料波传播过程中的波速变化和能量衰减.首先讨论线弹性接触模型、Hertz接触模型、线性平行粘结模型对不同排列颗粒材料中波传播的影响, 分析三种接触模型下的波速变化和能量变化; 同时研究颗粒弹性模量对颗粒材料中波传播的波速影响; 最后对颗粒材料样本采用Hertz接触模型, 探究了不同围压大小下颗粒材料中的波速变化和能量衰减.
1 离散元接触模型1.1 离散单元模型离散单元法将颗粒材料看作是各个满足运动方程的颗粒单元的集合.用离散单元法进行计算模拟时, 对每个颗粒单元建立模型, 计算出颗粒碰撞过程中受到的接触力和接触力矩, 通过追踪系统中每个颗粒单元来对整个物体进行分析.
接触模型是离散单元法的重要基础, 实质是准静态下颗粒的接触力学弹性分析结果, 对其分析计算直接决定了粒子所受的力和力矩的大小.
离散单元法的接触模型多种多样, 每种模型下接触力的计算方法也大多不同, 对不同的仿真对象, 必须建立不同的接触模型.下面简单介绍本研究中所用的接触模型.
1.2 线弹性接触模型图 1a是线弹性接触模型示意图.线弹性接触模型将两颗粒间的接触看作是刚度为k的线性弹簧与阻尼系数为β的阻尼器并联作用, 两颗粒之间的法向力Fn计算公式如下:
 | (1) |
 | 图 1 两种接触模型草图Fig.1 Sketch of two contact models (a)—线弹性接触模型；(b)—线性平行粘结接触模型. |
式中: kn是法向刚度; βn是法向阻尼系数; α是两颗粒的叠合量;
在线弹性接触模型中, 相互接触的两颗粒间的接触刚度与颗粒刚度间的关系可表示为
 | (2) |
1.3 线性平行粘结接触模型线性平行粘结接触模型在线弹性接触模型基础上增加平行粘结键, 如图 1b所示.与线弹性模型相比, 线性平行粘结模型增加了法向粘结拉力Fnc和切向粘结拉力Fsc, 其余参数与线弹性模型相同.当粘结断裂时, 线性平行粘结接触模型会变成如图 1a所示的线弹性接触模型.
平行粘结键可以设想为一组具有恒定法向刚度和剪切刚度的弹簧, 均匀分布在接触面内以接触点为中心的圆形截面上.这些弹簧与点接触弹簧平行(点接触弹簧通过线弹性模型在某个点上建立粒子刚度), 因此平行粘结键会导致颗粒材料内两个粘结颗粒间产生可以增大的力和力矩, 且应力在法向或剪切方向上满足临界力破坏标准.当粘结应力超过强度极限时, 粘结会断裂, 使界面不粘结, 接触行为变为线弹性.
1.4 Hertz接触模型Hertz模型是基于Mindlin和Deresiewicz理论[17]的非线性接触公式建立的近似计算模型.设两球形颗粒发生弹性接触, 颗粒间法向力Fn可由下式求出:
 | (3) |
设初始切向力为Fis, 在经过一个切向位移增量Δδs后(Δδs可由Vu-Quoc等[18]改进的切向力计算公式得出), 切向力记为Fsi+1, 则有
 | (4) |
 | (4a) |
2 数值样本为了比较不同接触模型在多种工况下对颗粒材料中波传播的影响, 本文采用三种不同排列方式的颗粒材料样品, 如图 2所示.其中, 等半径规则排列和等半径随机排列的颗粒半径为0.01 m, 不等半径随机排列的颗粒半径范围为0.007~0.013 m.
图 2(Fig. 2)
 | 图 2 不同排列颗粒材料样本抽样Fig.2 Sampling of granular materials with different arrangements (a)—等半径规则排列(R=0.01 m)；(b)—等半径随机排列(R=0.01 m)；(c)—不等半径随机排列(R=0.007~0.013 m). |
如图 3所示, 样本尺寸为2.4 m×1 m, 分别在其内部填充不同排列形式的颗粒材料, 并使整个颗粒材料受四周墙体边界约束.在样本中等间隔取8个点(从左到右依次设为P1, P2, …, P8)作为样本点.
图 3(Fig. 3)
 | 图 3 样本中的样本点和激发粒子Fig.3 Sample points and excited particles in the samples |
在加载之前, 系统处于静态平衡状态.通过对左边界部分粒子施加速度来激发应力波, 改变粒子间的接触模型, 并保证不同的接触模型下颗粒属性一致, 观察应力波的波速和动能变化.为了排除阻尼与摩擦的影响, 将粒子间的阻尼系数与摩擦系数都设为0, 颗粒材料密度为2 000 kg/m3, 弹性模量为50 GPa.
3 数值结果及分析3.1 波速变化波速v可通过测量脉冲传播一定距离l所需的时间t求出:
 | (5) |
图 4(Fig. 4)
 | 图 4 不同排列颗粒材料在不同接触模型下加载时波传播到样本点的波速Fig.4 Wave velocity of the wave propagating to the sample point under different contact models in different arrangement of granular materials (a)—规则排列；(b)—等半径随机排列；(c)—不等半径随机排列. |
由图 4分析可得, 波在颗粒材料中传播时, Hertz接触模型下的波速最小, 即波在Hertz接触下传播较慢.在颗粒规则排列时线性平行粘结模型下的波速大于线弹性接触模型; 但在颗粒随机排列时两者的波速大致相同, 这是因为随机排列的颗粒间距较大, 使线性平行粘结模型中的粘结键断裂, 此时其等效于线弹性模型.另外, 随机排列时各颗粒间的碰撞不规律, 碰撞后颗粒运动方向不一致, 导致波在传播方向外的各个方向弥散; 所以波在随机排列颗粒材料中传播时的速度低于规则排列时的波速.
3.2 能量衰减观察各样本点在加载过程中动能随时间的变化(见图 5), 以此比较波沿程的动能衰减.不同接触模型下各样本中各点动能峰值包络线如图 6所示.
图 5(Fig. 5)
 | 图 5 规则排列颗粒材料样本中各样本点的动能变化及峰值的包络线Fig.5 Kinetic energy change and peak envelope of each sample point in granular materials with regular arrangement (a)—线弹性接触；(b)—Hertz接触；(c)—线性平行粘结接触. |
图 6(Fig. 6)
 | 图 6 不同排列颗粒材料中波传播的动能变化Fig.6 Kinetic energy variation of wave propagation in different granular materials (a)—规则排列；(b)—等半径随机排列；(c)—不等半径随机排列. |
由图 5、图 6分析可得, 波在颗粒材料中传播时, 波的动能沿程呈指数衰减; 不同接触模型下波的动能衰减趋势和衰减幅度不同, 同时波的能量衰减也受颗粒材料排列方式的影响.由图 5可以看出, 在不同接触、不同排列颗粒材料中的波传播都在前两个点具有最大的动能衰减.为研究波的能量衰减特性, 选取图 5中各图前两个峰值进行比较, 计算能量衰减率:
 | (6) |
图 7(Fig. 7)
 | 图 7 不同排列颗粒材料加载时波的能量衰减率Fig.7 Energy attenuation rate of wave in different arrangements of granular materials |
由图 7可以看出, Hertz接触模型颗粒材料中的波传播在规则排列条件下的能量衰减率比其他两种接触模型高, 但其在随机排列条件下的能量衰减率比其他两种接触模型低; 线性平行粘结模型颗粒材料中的波传播在规则排列条件下能量衰减高于线弹性模型, 但在随机排列中两者基本一致, 这也印证了3.1节中随机排列下线性平行粘结模型由于间距过大导致粘结键断裂而等效于线弹性模型的结论; 随机排列下各接触模型中波传播的能量衰减高于规则排列, 这是因为随机排列中颗粒碰撞时的冲击角度比规则排列大, 导致冲击动能的能量耗散更多.
3.3 颗粒弹性模量的影响分别对图 2中三种排列方式的颗粒材料采用线弹性接触模型, 改变颗粒的弹性模量E(E∈(10 MPa~40 GPa)), 并保证三种排列方式颗粒材料中样本尺寸、加载方式、材料属性和接触模型都相同, 可得三种排列方式下波速随弹性模量E的变化曲线, 如图 8所示.
图 8(Fig. 8)
 | 图 8 不同排列方式颗粒材料中波速随颗粒弹性模量的变化Fig.8 Variation of wave velocity with elastic modulus in granular materials with different arrangements |
由图 8可知, 波在不同排列方式的颗粒材料中传播时, 波速都随着颗粒弹性模量E的增大而增加, 即v=exp(Aln E+B); 规则排列颗粒材料中的波速大于随机排列中的波速, 这与3.1节中的结论相同; 另外, 等半径随机排列颗粒材料中的波速大于不等半径随机排列中的波速, 这同样是颗粒碰撞的不规律性导致的, 不等半径随机排列中颗粒碰撞的不规律性更大, 波的弥散更快.
4 围压对颗粒材料中波传播的影响4.1 围压样本模型对颗粒材料采用规则排列, 颗粒之间的接触模型为Hertz接触模型.如图 9所示, 试样尺寸约为1 m×0.5 m, 其中颗粒半径为0.006 m, 整个颗粒系统受四周墙体边界约束, 在靠近左侧墙体处等间隔取6个点(从左到右依次设为P1, P2, P3, P4, P5, P6), 以测量平均波速和能量衰减.
图 9(Fig. 9)
 | 图 9 Hertz接触规则排列的颗粒材料样本模型及样本点和激发粒子Fig.9 Regularly arranged granular material samples with Hertz contact model, sample points and excited particles |
围压通过颗粒样本四周的墙体施加到颗粒体系中.首先将围压施加到墙体上, 围压分别取0.1, 0.5, 1, 3, 5 MPa, 并使施加到墙体的围压应力保持不变; 通过墙体内颗粒之间的相互作用, 使整个试件达到平衡状态.施加围压后, 在激发冲击载荷之前, 系统处于静平衡状态.通过对左边界部分粒子施加速度来激发波, 观察波传播过程中波速和动能变化.为了排除阻尼与摩擦的影响, 将粒子间的阻尼系数与摩擦系数都设为0, 颗粒材料密度为2 000 kg/m3, 弹性模量为10 GPa.
4.2 数值分析按照上述条件所得波速、动能变化和能量衰减率如图 10所示.
图 10(Fig. 10)
 | 图 10 不同围压下波传播的性能变化Fig.10 Performance change of wave propagation under different confining pressures |
由图 10a可以看出, 在采用Hertz接触模型规则排列颗粒材料样本中, 波速随着围压的增大而增加, 当围压在1~5 MPa范围内时, 每增大2 MPa围压, 波速增加约1%.由图 10b和图 10c可以看到, 不同围压加载条件下, 波在上述样本中传播时, 随着围压增大, 波沿程各点动能峰值包络线的下降幅度更大, 曲线更陡; 波传播的能量衰减率也会增大, 当围压在0.1~1 MPa范围内时, 每增大0.5 MPa围压, 能量衰减率增大约1%;当围压在1~5 MPa范围内时, 每增大2 MPa围压, 能量衰减率增大约2%.
5 结论1) 波在颗粒材料中传播时, Hertz接触下的波速最小, 线性平行粘结接触下的波速最大; 在随机排列颗粒材料中, 线性平行粘结接触由于间距过大导致粘结键断裂而等效于线弹性接触, 所以波在两者中的波速大致相同; 因为随机排列下各颗粒间碰撞具有不规律性, 碰撞后颗粒运动方向不一致, 导致波在传播方向外的各个方向弥散, 所以波在随机排列颗粒材料中传播时波速低于规则排列中的波速.
2) 波在颗粒材料中传播时动能呈指数衰减; 不同接触模型计算接触力的方式不同, 导致波的动能衰减趋势和衰减幅度不同; 同时排列方式影响粒子碰撞时的冲击角度, 由于随机排列中颗粒碰撞时的冲击角度比规则排列大, 导致冲击动能的能量耗散更多.
3) 改变颗粒的弹性模量E, 发现波速随着E的增大呈对数函数增大: v=exp(A ln E+B); 改变样本围压, 发现随着围压增大, 波速会增加, 能量衰减率也会增大.
参考文献
| [1] | 楚锡华. 颗粒材料的离散颗粒模型与离散-连续耦合模型及数值方法[D]. 大连: 大连理工大学, 2006. (Chu Xi-hua. Discrete particle model and discrete continuous coupling model and numerical method for granular materials[D]. Dalian: Dalian University of Technology, 2006. ) |
| [2] | 修晨曦, 楚锡华. 颗粒材料修正的微形态连续体模型研究[J]. 重庆大学学报, 2020, 43(6): 1-11. (Xiu Chen-xi, Chu Xi-hua. A modified micromorphology model based on micromechanics for granular materials[J]. Journal of Chongqing University, 2020, 43(6): 1-11.) |
| [3] | Cundall P A, Strack O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65. DOI:10.1680/geot.1979.29.1.47 |
| [4] | Wang J, Chu X H. Compressive wave propagation in highly ordered granular media based on DEM[J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2018, 19(5): 545-552. DOI:10.1515/ijnsns-2017-0213 |
| [5] | Wang J, Chu X H, Xiu C X, et al. Stress wave in monosized bead string with various water contents[J]. Advanced Powder Technology, 2020, 31(3): 993-1000. DOI:10.1016/j.apt.2019.12.027 |
| [6] | Sadd M H, Gao J, Shukla A. Numerical analysis of wave propagation through assemblies of elliptical particles[J]. Computers and Geotechnics, 1997, 20(3): 323-343. |
| [7] | Wang J, Chu X H. Impact energy distribution and wavefront shape in granular material assemblies[J/OL]. Granular Matter, 2019, 21: 23[2020-12-15]. https://link.springer.com/article/10.1007/s10035-019-0880-z. |
| [8] | 修晨曦, 楚锡华. 基于微形态模型的颗粒材料中波的频散现象研究[J]. 力学学报, 2018, 50(2): 315-328. (Xiu Chen-xi, Chu Xi-hua. Study on dispersion behavior and band gap in granular materials based on a micromorphic model[J]. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 315-328. DOI:10.7638/kqdlxxb-2016.0040) |
| [9] | Xiu C X, Chu X H, Wang J. Prediction on dispersion in elastoplastic unsaturated granular media[J]. Theoretical & Applied Mechanics Letters, 2020, 10(2): 74-78. |
| [10] | Sadd M H, Tai Q M, Shukla A. Contact law effects on wave propagation in particulate materials using distinct element modeling[J]. International Journal of Non-Linear Mechanics, 1993, 28(2): 251-265. DOI:10.1016/0020-7462(93)90061-O |
| [11] | 黄德财, 陈伟中, 杨安娜, 等. 孤立波在一维复合颗粒链中传播特性的模拟研究[J]. 物理学报, 2014, 63(15): 266-271. (Huang De-cai, Chen Wei-zhong, Yang An-na, et al. Simulation study on the propagation of solitary waves in a one-dimensional composite granular chain[J]. Acta Physica Sinica, 2014, 63(15): 266-271.) |
| [12] | Bi Z W, Huang W F, Sun Q C, et al. Experimental analysis of cyclic loading on a cohesionless granular system[J]. Granular Matter, 2016, 18(1): 1-9. DOI:10.1007/s10035-015-0597-6 |
| [13] | Daraio C, Nesterenko V F, Herbold E B, et al. Tunability of solitary wave properties in one-dimensional strongly nonlinear phononic crystals[J/OL]. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2006, 73(2)[2020-12-15]. https://journals.aps.org/pre/abstract/10.1103/PhysRevE.73.026610. |
| [14] | Wang J, Chu X H, Jiang Q H, et al. Energy transfer and influence of excitation frequency in granular materials from the perspective of Fourier transform[J]. Powder Technology, 2019, 356: 493-499. DOI:10.1016/j.powtec.2019.08.061 |
| [15] | 季顺迎, 樊利芳, 梁绍敏. 基于离散元方法的颗粒材料缓冲性能及影响因素分析[J]. 物理学报, 2016, 65(10): 168-180. (Ji Shun-ying, Fan Li-fang, Liang Shao-min. Buffer capacity of granular materials and its influencing factors based on discrete element method[J]. Acta Physica Sinica, 2016, 65(10): 168-180.) |
| [16] | Wang J, Chu X H, Zhang J B, et al. The effects of microstructure on wave velocity and wavefront in granular assemblies with binary-sized particles[J]. International Journal of Solids and Structures, 2018, 159(2): 156-162. |
| [17] | Mindlin R D, Deresiewicz H. Elastic spheres in contact under varying oblique forces[J]. Journal of Applied Mechanics, 1953, 20(3): 327-344. DOI:10.1115/1.4010702 |
| [18] | Vu-Quoc L, Zhang X. An accurate and efficient tangential force-displacement model for elastic frictional contact in particle-flow simulations[J]. Mechanics of Materials, 1999, 31(4): 761-762. |
