表1(Table 1)
引用方式 |
 Download
PDF  |
2.中国科学院过程工程研究所多相复杂系统国家重点实验室, 北京 100190;
3.北京科技大学物理系, 北京 100083
2.State Key Laboratory of Complex systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China;
3.Department of Physics, University of Science and Technology Beijing, Beijing 100083, China
颗粒物质(granular matter 或 granular materials)由大量颗粒组成是多体相互作用的无序体系(Jaeger et al.1996,汪卫华 2013,Cheng & Ma 2011),具有独特力学性质(李家春 2014).颗粒物质广泛应用于工程建设、工业生产和自然界中. 比如,在水利水电工程中堆石坝由块石堆积而成, 工程中十分关心其静态与动态响应. 在火电领域,煤炭的储存和输运涉及颗粒流问题. 我国重点发展的球床模块式高温气冷堆(简称高温气冷堆), 几百万球形燃料颗粒紧密堆积, 在反应器中的径向和轴向速度分布是堆芯几何设计、反应堆物理以及核燃料循环计算的基础, 现阶段主要基于唯象分析和经验判断. 反应堆长年运行, 不会一成不变地遵守预设流速分布, 如果出现设备故障, 或者颗粒的性质变化乃至颗粒破裂, 都会改变球流特性, 造成过热或过燃耗, 运行工况将异常复杂,高温气冷堆领域对颗粒流的科学认识需求极为迫切. 还比如, 中国科学院“未来先进核裂变能--ADS 嬗变系统” 战略性先导科技专项中的重金属散裂靶设计, 人们创造性地提出了新型流态固体(比如毫米级钨) 颗粒靶概念并完成初步设计, 与电子束耦合的小型台架原理性实验获得成功. 该方案通过固体颗粒的流动实现了靶区外的冷却, 规避了液态铅铋合金靶放射产物毒害性高、温度-材料腐蚀效应严重以及固态靶热移除难等缺点, 物理上可承受几十兆瓦束流功率. 涉及到的核心技术是确定靶容器中的颗粒流动速度和颗粒间的传热等, 与之相关的科学问题就是确定颗粒流流变性质、颗粒传热系数等.
当研究问题的尺度远大于颗粒粒径时,颗粒体系通常处理成等效连续介质(Li & Li 2009). 长期以来,人们已经建立一批有实际应用价值的颗粒介质宏观力学性质的本构关系模型(亦即应力-应变(率)关系)也发展了比较齐全的测试设备和技术,并基于有限元、有限差分等数值方法,解决了大量工程问题.基于经验和实验数据而建立的本构关系反映的是颗粒体系在一定边界条件下的力学行为因此具有一定的适用范围,针对不同过程多种模型共存的状态将持续很长一段时间. 另外,人们在颗粒试验中发现由于试样的初始结构差异,每次试验结果都有变化特别是颗粒体系软化后的测量结果,如此看来,片面追求本构关系模型的精准性并不合适(黄文熙 1983,Kolymbas & Wu 2000,沈珠江 2000).
与经典弹塑性材料相比,颗粒介质具有显著的离散性,体系内部的结构不均匀性(structural heterogeneity),导致动力学的不均匀性(dynamical heterogeneity),如变形局部化(strain localization)、非局域流变(non-local rheology)等,不恰当的平均化处理(homogenization)掩盖了这些不均匀性,很难解释实际颗粒物质表现出的异常复杂的宏观性质,比如人们争论的核心是采用一个还是多个屈服面? 是否采纳Drucker 公设?用关联还是不关联的流动法则等等.在缺乏对颗粒材料内部物理过程、变形机制、失稳机理和能量耗散有较为清晰认识的情况下,基于经验建立本构关系模型的进一步发展将面临很大难度,争议也将不断上演(李广信 2006). 临渊羡鱼,不如退而结网开展颗粒介质内在物理过程的基础研究,才能建立科学合理的本构关系(孙其诚和王光谦 2009,孙其诚等 2011).
颗粒体系是多体相互作用体系,构成多尺度结构,具有多物理过程和复杂力学性质. 颗粒体系涉及3 个尺度,亦即颗粒内部分子、电子等的微观运动,以若干颗粒大小为特征长度的介观运动以及颗粒体系的宏观运动. 不同层次的运动具有不同物理机制,不同层次间的运动又有一定的关联. 颗粒表面摩擦和颗粒非弹性变形,机械能部分的转化为颗粒热能,因此颗粒体系介观尺度的运动与微观运动紧密关联; 导致与此同时,介观力链位形演化,引起宏观层次力学性质的改变,比如应力陡降、塑性变形和流变等. 介观层次的运动特征是颗粒物质与普通流体和固体的主要区别所在,后者主要涉微观和宏观两个层次的运动,可直接采用统计力学或动理学等微观理论,推导宏观热力学势和其他物理量,比如微观分子的无序运动在宏观层次上的表现形式就是熵或温度. 颗粒体系介观层次的运动主要由场力(如重力)、颗粒接触力、间隙流体作用力等控制. 颗粒间的相互作用能(或者自身重力势能) 远大于kBT (kB 是Boltzmann 常数,T 是热力学温度),微观层次的热运动不足以使得介观层次的颗粒发生运动,必须依靠振动和剪切等机械方式激发,颗粒间才会发生法向脱离或者切向滑动,导致颗粒体系从一个位形(configuration) 转变到其他位形. 我们认为,决定颗粒体系宏观性质的不仅是组成它的颗粒,更大程度上取决于这些颗粒相互作用形成的介观结构. 因此,开展多尺度分析,探究颗粒介质变形与破坏的机制,探测、表征复杂力学性质的介观结构起源具有重要的理论价值,是现阶段颗粒物质的力学本构关系研究的主要任务,并在设计、调控工程中的颗粒材料力学行为方面具有重要的现实意义.
目前,凝聚态物理学中的“应力/应变-驱动玻璃转变”的思想极大启发和推动了颗粒材料力学性质多尺度的力学研究,已经成为理解颗粒材料力学性质(强度和韧性)、塑性变形机制、剪切带起源和流动的主流思想之一.在结构表征、宏观特性、结构和性能的关系、基础理论、新技术探索方面报导了大量数据和现象,提出了大量概念和模型,同时也存在大量的问题有待澄清和研究.本文介绍了颗粒接触力模型、结构特征和演化规律等方面的研究进展,最后展望了可能取得突破的研究课题,也希望促进颗粒介质力学研究和工程应用的交叉与融合.
2 颗粒接触力模型颗粒离散元方法(discrete element method,DEM)是研究颗粒体系运动规律的一种数值方法,它建立在牛顿第二运动定律上,通过更新每个颗粒的速度和位置,进而确定性的演化整个颗粒系统(孙其诚和王光谦 2008).自20世纪70年代提出以来离散单元法迅速在散体材料领域得到广泛的应用,并形成了诸如PFC,EDEM,YADE,LAMMPS,Esys-particle,dp3D等多种商业和开源软件.相比于各种连续介质模型,离散单元法的优势是能够直接从颗粒尺度刻画颗粒系统的宏观动力学响应特性,这为认识颗粒材料介观尺度结构和洞察新现象等方面发挥了强大作用,为揭示颗粒材料复杂宏观力学的介观机理奠定了基础.
随着颗粒离散元方法在不同领域的广泛应,用颗粒间接触力的简化模型发展较快. 依据颗粒间作用处理方式的不同颗粒间接触力模型可分为硬球模型和软球模型.硬球模型假设体系中颗粒碰撞是两体、瞬时的发生接触后颗粒的速度、角速度等立即发生改变.软球模型最早由Cundall在1971年提出并由Cundall和Strack于1979年发展建立(Cundall,1971,Cundall & Strack,1979),假设相互接触的颗粒间可发生重叠计算颗粒间接触力和扭矩等. 迄今为止研究者建立了多种不同的软球接触力模型. 本节评述软球模型暂不涉及颗粒间长程作用及流固耦合作用力等模型.
2.1 法向接触力模型
颗粒接触力的分析源自经典的无粘着接触力学理论,亦即Hertz 接触理论.当粘着作用不可忽略时,Johnson-Kendall-Roberts(JKR)接触理论在Hertz接触理论基础上,考虑了颗粒接触表面内分子或原子间作用力,解释了卸载时颗粒在接触界面的颈缩现象. Derjagin-Muller-Toropov(DMT)接触理论考虑了接触面以外的颗粒表面间的范德华力,当颗粒间表面分离时作用力简化到Bradley 理论.这些经典接触理论的适用范围如图 1 所示,图中两个相互接触的球形颗粒有效半径为R,界面粘着能为 △γ,法向接触力为N. λ是Tabor 数,
|
| 图 1 Hertz接触理论、JKR接触理论、DMT接触理论和Maugis-Dugdale 理论的适用范围 |
严格来说,Hertz接触理论对于各向同性弹性球体静态接触是准确的.对于其他情况,比如多边形或多面体,缺乏相应的接触理论这就无法提出普适且有效的任意形状颗粒接触力公式. 在绝大部分情况下比如考虑颗粒形状、大小、表面摩擦材料参数等以及特定条件下,必须开展实验进行测量,建立相应的接触力计算公式.上面提及的所有接触理论都针对弹性接触,事实上颗粒在接触过程中发生或多或少的塑性变形. 在一次加载时,接触面的变形超过弹性极限而出现塑性变形,将产生残余应力.一般情况下,残余应力提高了下一次加载产生塑性变形的临界载荷,经过多次重复加载作用,塑性变形的临界载荷将达到稳定值.切向力不仅取决于法向力大小,还与切向力加载历史相关.由于形变和接触力之间是非线性的,法向力和切向力又耦合在一起,很难得到法向力和切向力的分析解,所以通常采用增量方法分别计算法向力和切向力.
法向接触力的简化模型大致分为四类: 连续势函数模型(continuous potential model)、线性黏弹模型(linear viscoelastic model)、非线性黏弹模型(non-linear viscoelastic model)和滞回模型(hysteretic model).
势函数模型: 该模型借鉴分子动力学(molecular dynamics,MD)势函数模型来计算颗粒间作用力,在物理类期刊中离散单元法有时也被称为分子动力学方法(Zhu et al.2007,Rognon et al.2008,Tejada & Jimenez 2014).Aoki和Akiyama(1995)首次将用于模拟原子或分子间作用的连续势函数引入到颗粒介质中,颗粒被处理为质点,颗粒间只存在法向接触力fn通过对势函数φ求导计算得到,同时引入与颗粒间相对速度呈正比的阻尼力来反映碰撞前后的动能损失
|
(1) |
式中,nc是发生接触的颗粒中心连线的单位矢量,η是阻尼vij是两颗粒间的相对速度.
人们发展了不同的连续势函数φ形式(Langston et al.1994,Baxter et al.1997, Liu et al.2008)共同特征是颗粒间接触势随颗粒距离rij的增大而衰减,
线性黏弹模型: 该模型是Cundall和Strack(1979)提出DEM方法时采用的接触力模型,也是目前应用较为广泛的模型,该模型中假设颗粒间的弹性接触力与颗粒间重叠量 δ n呈正比
|
(2) |
采用线性黏弹模型时,颗粒接触被处理为带有阻尼的谐振子,此时颗粒间弹性恢复系数(e)和颗粒接触时间(tc)具有解析解(Schafer et al.1996,Hoomans 1999)
|
(3) |
|
(4) |
其中,meff=mimj/(mi+mj).在实际应用中,为保证离散积分的准确性,通过式(4)预估出颗粒间最小接触时间tc极为重要,DEM模拟中的时间步长必须要远小于tc. Hoomans(1999)建议一次碰撞时间处理为20~40个时间步长,基本能保证碰撞后的颗粒速度与通过动量定律计算得到的数值近似相等.需要指出的是,采用式(2)所示的线性黏弹模型时,颗粒接触时间tc与颗粒间弹性恢复系数有关. 如果颗粒间发生过阻尼,在极端情况下颗粒接触时间tc会非常长,造成不真实的结果(Luding et al.1994).
与非线性黏弹模型相比,线性黏弹模型中没有考虑颗粒接触面对颗粒间接触力的影响,且弹簧刚度系数(stiffness)kn不能由颗粒材料的杨氏模量和泊松比直接计算得到.尽管线性黏弹模型偏离了Hertz理论,但是大量的数值模拟结果表明采用线性黏弹模型的DEM模拟仍能给出定性合理的宏观颗粒流动力学特性.这可能是因为颗粒材料是一个非线性、多自由度复杂系统,涉及更多的物理过程,单纯由于颗粒接触力学不精确造成的误差被减弱,显然仅仅依靠模拟现象与宏观现象的一致性,并不能必然得到颗粒接触作用力模型准确的结论.
非线性黏弹模型: 该模型一般是基于Hertz接触理论计算碰撞颗粒间的接触力,这类模型虽然较线性模型复杂,但可以更准确地描述颗粒间的接触作用,且模型中弹簧刚度系数可直接由颗粒的杨氏模量和泊松比计算得到
|
(5) |
其中,
不同非线性黏弹模型的一个主要区别在于阻尼力的计算,即指数a的选取. 阻尼力用于反映颗粒接触过程中的能量耗散,因此其与颗粒弹性恢复系数e密切相关. 然而对于实际的颗粒系统,颗粒弹性恢复系数e不仅取决于颗粒本身性质,还与颗粒间相对速度有关(Cooke & Bridgwater,1982,Hussainova et al.1999,Thornton & Yin 1991,Zhou 2013). 基于所关注的核心物理、力学问题的不同,不同学者针对式(5)中阻尼力的计算选择了不同的a值,包括0( Lee &Herrmann,1993),1/4(Tsuji et al.1992),1/2(Kuwabara& Kono 1987 ,Brilliantov et al.1996). 当a=0时,由式(5)得到的弹性恢复系数随相对速度的增大而增大,接触时间随碰撞速度的增大而缩短. 当a=1/4时,式(5)预测的碰撞恢复系数为常数,但颗粒间接触时间与碰撞速度相关.
滞回模型: 该模型主要用于模拟塑性颗粒体系. 在模型表达形式上其与线性/非线性黏弹模型的最大区别在于将颗粒间的碰撞作用分为不同的阶段,每个阶段采用不同的弹簧刚度系数,进而基于线性关系或非线性关系计算颗粒间接触力.典型的滞回模型包括以下3种形式.Walton和Braun(1986)将塑性颗粒间的接触作用分解为加载阶段和卸载阶段,并假设每个阶段力-重叠量满足线性关系,提出了以下法向接触力作用模型
|
(6) |
式中,kl,kul分别为加载、卸载过程中的弹簧刚度系数, δ n0为该碰撞接触所造成的颗粒塑性变形. 由于式(6)的线性形式采用该模型时颗粒塑性变形 δn0、弹性恢复系数e和颗粒接触时间tc具有解析解
|
(7) |
当kl和kul取常数时,该模型给出恒定的弹性恢复系数和颗粒接触时间,而颗粒塑性变形则还与碰撞前颗粒间相对速度v0有关.
Sadd和Tai(1993)提出了一种表达形式与Walton-Braun模型类似的滞回模型
|
(8) |
其中,p为经验参数,q则根据加载阶段的最大重叠量计算得出 q =(A δ n,max)2,其中A为可调参数. ψ由最大重叠量δn,max和塑性变形量 δ n0计算得到: ψ =( δ n- δ n0)/ δ n,max- δ n0).与Walton-Braun模型相比,Sadd-Tai模型所预测的颗粒弹性恢复系数与颗粒重叠量有关,且考虑了颗粒塑性变形对下次碰撞作用的影响,因此更接近于真实塑性颗粒之间的作用,但所需要的经验参数更多.
实际颗粒在发生塑性变形之前会经历一弹性变形过程,基于这一事实Thornton(1997)将塑性颗粒之间的碰撞作用分为三个阶段:加载过程的弹性变形阶段、加载过程中发生塑性屈服后的塑性变形阶段、卸载过程,进而提出以下模型
|
(9) |
在该模型中,颗粒间作用达到屈服应力Py后发生塑性作用. 屈服应力Py与材料性质有关,但在多数的DEM模拟中该参数一般以可调参数的形式出现.发生应力屈服时的颗粒间重叠量可由下式计算得到: 
目前,对滑动摩擦已有深刻的认识,而对转动摩擦以及转动－滑动摩擦转换机制的认识还很不足. 一般而言,转动摩擦的机理与滑动摩擦不同,当接触物体间发生相对转动或有相对转动趋势时,在接触面间产生阻止相对转动的转动摩擦力矩,对颗粒物质宏观力学行为有重要的影响,比如对颗粒体系的峰值强度起着控制作用,还会导致较大的体积膨胀量.硬质材料的圆柱体或者球体之间的滚动或者它们在平面上滚动时滚动摩擦系数测量非常困难,不同文献中的取值甚至有数量级的差别,相比之下滑动摩擦系数取值较为一致. 需要注意的是人们熟知的滑动摩擦定律只适用于很少数的情况,只有对于给定的材料和确定的工况条件(温度、适度、压力和滑动速度)摩擦系数才是确定数值,许多材料的静、动摩擦系数都与载荷和滑动速度有关. 因此,一定要谨慎地选用任何文献报道的摩擦系数.
相比于法向接触力,切向接触力的计算模型更为复杂.依据Mindlin-Deresiewicz理论(Mindlin & Deresiewicz,1953)切向力不仅与接触历史有关,还与法向力大小有关.切向接触力理论推导上分法向力恒定切向力变化和法向、切向力同时变化两种情况Vu-Quoc等 (2001)对此进行了详细的评述. Vu-Quoc和Zhang,(1999)及Di Renzo和Di Maio,(2004)依据Mindlin-Deresiewicz理论详细推导了DEM模拟中切向力的计算.然而在实际大多数应用中,由于Mindlin-Deresiewicz理论过于复杂,切向力更多时候是通过对Mindlin-Deresiewicz理论进行简化计算.与法向力类似,切向线性或非线性力-重叠量关系都有采用,不同简化模型的最主要区别在于切向弹簧刚度系数的选取和计算.
Walton和Braun(1986)基于Mindlin-Deresiewicz理论提出了以下切向弹簧刚度系数计算公式
|
(10) |
式中,kt(0)是刚发生接触时的切向弹簧刚度系数,μ 为摩擦系数,Ft*(t)初始时刻设置为0,当颗粒间相对运动方向发生改变后,设置为相对运动方向发生改变时刻的Ft(t); 在碰撞过程中如果法向力发生改变,则Ft*(t)同比例地变化. 基于式(10)计算得到的切向弹簧刚度系数随相对切向位移的增大而减小.
Thornton和Yin(1991)提出了一更复杂的切向弹簧刚度系数计算公式
|
(11) |
其中,a为接触面积,
同样是基于Mindlin-Deresiewicz理论Tsuji等 (1992)提出另一种切向弹簧刚度系数计算公式
|
(12) |
该模型仍反映了切向力与法向力直接相关这一特点,与Walton-Braun模型与Thornton-Yin模型相比,最主要区别是将颗粒间接触作用处理为单一过程. 在恒定法向力情况下,式(12)给出的切向弹簧刚度系数在加载和卸载阶段是相等的,因此理论上该模型只适用于完全弹性颗粒接触. 与Tsuji模型表达形式类似,同样被广泛采用的另一类切向弹簧刚度系数模型为Langston等 (1994)提出的模型:
|
(13) |
式中fn,e 表示法向接触力中的弹性项,
以上模型的共同特点是切向弹簧刚度系数是法向力或者法向重叠量的函数. 在实际应用中,另一类被广泛采用的方式是将切向弹簧刚度系数设置为常数,这也是Cundall和Strack(1979)提出DEM方法时采用的方法. 与线性黏弹法向接触力模型类似,尽管这种简化显著偏离了Mindlin-Deresiewicz理论,但大量的离散单元法数值模拟结果表明,采用恒定切向弹簧刚度系数仍能定性地复现实验中所观测到的颗粒材料各种复杂宏观现象,见Zhu等 (2008)关于离散单元法应用综述性文章中的相关文献. 值得注意的是,依据Hertz-Mindlin理论,对于给定的材料,在颗粒间发生碰撞的初始阶段切向与法向弹簧刚度系数的比值取决于材料的泊松比v(kt/kn = 2(1 -v)/(2-v))(Maw et al.1976). 常见材料的泊松比一般位于(0,0.5)区间因此理论上切向弹簧刚度系数应满足2/3≦kt/kn≦1.然而在很多离散单元法数值模拟中切向弹簧刚度系数的数值选取都带有一定的随意性,kt/kn 的取值从0.1到1.0都有. 譬如在Xu和Yu (1997)的模拟中,kt/kn =1;Hoomans(1999)则依据法向弹簧和切向弹簧振荡周期应该相等这一原则建议kt/kn = 2/7; 而Luding,(1998)则建议kt/kn 应小于2/7.切向弹簧刚度系数的选取会显著影响碰撞后颗粒速度,以及颗粒接触动力学(Luding,2008a,Kruggel et al.2008,Thornton et al.2011,2013). 因此当研究目标为颗粒尺度的力学响应特性时,切向弹簧刚度系数模型,及数值的选取必须慎重.
当切向力大于某一个临界值时,颗粒之间会发生摩擦滑移.因此对于所有的切向力模型,都有以下限制条件
|
(14) |
其中μ 为平动摩擦系数.
需要指出的是,尽管离散单元法存在多种不同的颗粒法向/切向接触力模型,然而不同模型的选取,以及法向和切向接触力模型的搭配上,目前并没有统一的标准. 定性上,采用不同接触力模型都能给出与实验结果相吻合的宏观模拟结果,见离散单元法应用方面综述性论文及相关文献(Zhu et al.2008,Lu et al.2015,Mishra,2003). 然而在颗粒尺度,采用不同接触模型所得到颗粒尺度动力学响应特性可能会存在差异,多位学者对此进行了详细的评述和模型比较(Vu-Quoc et al.2001;Vu-Quoc & Zhang,1999; Mishra,2003a,2003b; Kruggel-Emden et al.2007; Di Renzo & Di Maio 2004). 此外,几乎所有的离散单元法文献报道中,所采用的接触力模型都是基于两体接触力学理论. 当颗粒间发生多体接触且接触面积不可忽略,譬如非规则块体颗粒间的面接触,或者颗粒绑定在一起形成固体桥键时,即使是采用完整的Hertz理论和Mindlin-Deresiewicz理论,模拟结果仍可能与实际情况存在较大的定量差异(Jefferson et al,2002). 这主要是因为此时同一个颗粒不同接触面间存在应力耦合效应,为得到定量准确的模拟结果,模拟中这种应力耦合效应必须要加以考虑(Liu et al.2010,2011a,2011b).
2.3 转动阻力矩模型由于颗粒间切向接触力作用在颗粒表面而不是颗粒质心,因此颗粒间切向作用会导致颗粒发生旋转.早期的离散单元法模拟都假设颗粒具有完全旋转自由度,即颗粒在转动方向上不受限制,转动加速度完全取决于颗粒间切向作用力. 然而实际颗粒并不是理想光滑的球形颗粒,且颗粒接触是以面而非点的形式存在因此颗粒的旋转能力必然会受到限制.Zhou等 (1999)的研究结果表明,DEM模拟中如果不对颗粒的转动加以限制,则无法准确地预测颗粒材料静态堆积时的静态休止角. 另一方面如果完全限制颗粒转动则也同样得不到定性准确的结果.譬如Wang等 (2006)的研究结果表明: 采用DEM模拟脆性材料的破碎时,只有在考虑颗粒旋转自由度的情况下才能复现出实验中所观测到的翅型裂纹形式. 对于实际的颗粒材料,由于颗粒接触面上的微观滑移/摩擦作用、颗粒塑性变形、颗粒间粘附力、颗粒的非规则形状等因素(Ai et al.2011),颗粒在受到切向力作用时发生的往往是受限旋转,即颗粒受到一转动阻力矩Mr
|
(15) |
文献中转动阻力矩模型大致可分为两类,分别以澳大利亚新南威尔士大学余艾冰课题组提出的转动摩擦系数模型(Zhou et al.1999,Zhu & Yu,2003)和Iwashita和Oda,(1998)提出的转动弹簧模型为代表. 余艾冰课题组采用DEM模拟研究颗粒材料堆积现象时,为了反映颗粒的有限转动,提出了以下两种转动阻力矩模型
|
(16) |
其中μr 为转动摩擦系数. 模型A中转动阻力矩与法向接触力中的弹性项呈正比,且转动阻力矩始终与颗粒转动方向相反. 模型B与模型A的主要区别在于转动阻力矩不仅取决于法向接触力弹性项,还与颗粒间相对转动速度有关. 之后的工作中他们对模型A和模型B进行了修正,提出以下模型(Zhu & Yu,2003)
|
(17) |
式(17)实际上是考虑转动自由度上的摩擦效应后式(15)和(16)的组合,只是此时转动刚度系数μ'r取为常数. 需要指出的是式(17)中转动阻力矩是颗粒间相对转动速度而不是相对转动位移的函数,这也是该模型与Iwashita和Oda(1998)提出的转动弹簧模型的主要区别.
Iwashita和Oda(1998)在研究颗粒材料剪切液化现象时发现,只有在DEM中考虑颗粒有限转动才能合理复现实际剪切带中的大孔穴现象.他们借鉴切向力的处理方式,在转动方向上引入转动弹簧、阻尼器和转动摩擦器,进而提出以下模型
|
(18) |
式中,kr,ηr ,μr分别为转动弹簧刚度系数、转动阻尼系数、转动摩擦系数; θ r 为角位移. B为与颗粒接触面积相关的特征尺度(Iwashita & Oda 2000). 由于颗粒间转动摩擦系数测量比较困难,一般以可调参数出现,因此在一些离散单元法模拟中B直接取为颗粒半径(Tordesillas等,2008).
采用式(18)所示的转动弹簧模型时需要计算颗粒间的角位移θ r,目前文献中有两种方法.
一种是基于颗粒的旋转和颗粒间质心相对位置变化,将颗粒接触点的空间变化分解为切向位移和角位移(Iwashita & Oda 1998). Jiang等 (2005)详细推导了角位移的计算公式.
另一种方法是直接基于颗粒间相对角速度计算角位移,即
转动弹簧模型的另一个重要参数是转动弹簧刚度系数kr.通过假设转动阻力矩应该与切向力所导致的转动力矩处于同一量级,(Iwashita和Oda1998,2000)采用如下公式计算kr
|
(19) |
即转动弹簧刚度系数是切向弹簧系数和颗粒半径的函数. 对于不同粒径的颗粒接触对,Tordesillas等 (2008)建议式(19)中的颗粒平均粒径以小颗粒的半径代替.Jiang等 (2005)假设颗粒接触为系列弹簧的组合,基于理论分析认为转动弹簧刚度系数应与法向弹簧刚度系数有关
|
(20) |
其中ψ为与材料性质和颗粒形状有关的经验常数.基于Hertz理论,Bardet和Huang(1993)提出以下公式来计算转动弹簧刚度系数
|
(21) |
与式(19)和式(20)相比式(21)给出的转动弹簧刚度系数是法向接触力的函数.
与接触力模型类似,目前并没有统一的转动阻力矩模型及相关参数选择标准,各种转动阻力矩模型都能定性地复现实验所观测到的宏观现象,但所预测的颗粒尺度的动力学响应特性存在一定的区别(Ai et al.2011).
2.4 阻尼对于实际的颗粒材料,所输入的能量通过颗粒间的碰撞作用及摩擦滑移而耗散.DEM模拟中颗粒碰撞所耗散的能量通过阻尼器来实现,即法向和切向接触力以及转动阻力矩包含一弹性项和一阻尼项(contact damping),见前文中的相应公式. 阻尼项中的阻尼因子η一般基于以下公式计算
|
(22) |
式中i表示法向、切向或者转动方向. 对于转动方向的阻尼因子,式(22)中m表示转动惯量. α与弹性恢复系数有关,当采用线性黏弹模型时,其具有如下解析解
|
(23) |
由于单个颗粒的质量较小,因此一般情况下
当模拟准静态颗粒系统时,为确保系统处于拟平衡状态和减少迭代步数一般还需要引入局部阻尼(local damping). 在离散单元法商业软件PFC中采用以下方式实现局部阻尼
|
(24) |
式中,i表示方向,Fi 表示颗粒在i方向所受到的合力,Fd为局部阻尼力其大小与局部阻尼因子α有关,PFC中α的默认值为0.7.式(24)所示的局部阻尼模型的优点是只有处于非力平衡状态的颗粒才受到局部阻尼作用、局部阻尼因子α为常数、以及局部阻尼与频率无关.
另一类局部阻尼模型是将局部阻尼力设置为颗粒速度的函数(Luding,2008b & Strack,1979,Luding,2008b)
|
(25) |
式(25)与接触阻尼的表达形式类似,但采用的是颗粒速度而非相对颗粒速度. 与式(24)类似,式(25)中的局部阻尼因子β一般需要通过反复试算才能得到. β过小则增加达到平衡态所需的迭代次数,过大则会导致人为的过阻尼现象(Luding,2008b).Cundall,(1982) 应用伺服机理,通过考察能量变化对模拟中β进行自适应调整
|
(26) |
式中,Ed为局部阻尼所耗散的能量,
需要特别指出的是,采用局部阻尼的前提是颗粒系统处于惯性可忽略的准静态,对于处于快速流动的颗粒系统,或者目标参数与颗粒速度直接相关的颗粒系统,采用局部阻尼会导致不合理的结果.
2.5 时间步长
在确定接触力模型之后,离散单元法是把速度和加速度对时间步长进行离散积分,进而得到颗粒的位置和速度变化. 因此,时间步长△t的选取直接决定了计算效率和准确度. 在接触时间具有解析解的条件下,时间步长的选取一般基于以下公式: △t=tc/N,N表示处理一次碰撞需要的时间步数,一般取为20~40(Hoomans,1999); 在无解析解的情况下,一般通过以下公式近似颗粒接触时间: 
在实际应用中,受计算能力的限制往往需要将模拟时间步长放大以加快模拟进度.对于准静态颗粒系统,常见的策略是将颗粒质量放大(Liu et al.2012).对于动态颗粒系统,一般采用降低弹簧刚度系数的方式增大时间步长.对于给定的颗粒间相对速度,降低弹簧刚度系数意味着颗粒变"软"颗粒间重叠量变大. 不同学者对颗粒间重叠程度的上限有不同的看法,从要求重叠量小于颗粒直径的0.5%(Misra & Cheung,1999)、1%(Ye 2005)到小于颗粒粒径的10%(Walton,1993)都有.降低弹簧刚度系数带来的主要问题是,在计算与颗粒接触面积相关的物理量时,必须做相应修正.譬如在采用DEM模拟颗粒传热时,颗粒间接触面积必须要采用实际的杨氏模量进行校正以准确计算颗粒间传热量(Gan et al.2016). 此外,弹簧刚度系数的不同取值也有可能会直接影响宏观统计结果.譬如Chialvo等 (2012)的模拟结果表明,DEM模拟得到的颗粒材料宏观摩擦系数与所选择的弹簧刚度系数密切相关.因此,在采用较小弹簧刚度系数以增大时间步长时必须要慎重.
2.6 接触力测量技术在密集颗粒体系中,相互挤压接触的颗粒形成力链网络,其空间分布是不均匀的,强的接触力主要沿着链状路径(亦即力链)结构传递.人们认为力链网络是确定颗粒物质力学性质如稳定性、弹性和塑性等的重要因素. 近年来,人们已经从实验、数值模拟和理论方面对颗粒体系中接触力概率密度函数(P(F))进行了系统的研究,并取得了一些共识.颗粒体系力学性质的测量技术和数据采集技术取得了较大进步,开发了许多高精度设备. 通常在颗粒体系上有控制的施加干扰,研究颗粒体系所呈现的规律,进而出推算体系内部的接触力、应力分布.颗粒接触力的测量方法大致可分为两类:(1) 接触式检测方法,包括高精度电子传感称量法、显色灵敏压痕方法等,它们可检测颗粒体系中某一截面上的接触力分布情况,但是不可避免地对颗粒体系带来了或多或少的干扰,由于力链对局部力的变化反应极为敏感,检测引起的轻微改变足以使得力链结构发生很大变化;(2) 非接触式检测方法主要包括光弹方法、荧光共聚焦显微镜法和磁共振弹性成像法.它们的优势在于无干扰检测,不仅可以确定颗粒的几何位置,分析体系的结构特征,更重要的是可以"看到"力和力链,进而统计得到颗粒体系中接触力分布的全部信息.常用的颗粒接触力测量方法汇总见表 1 中.
近期,美国Duke大学 Behringer 课题组结合光弹技术,采用高速相机拍摄了颗粒受冲击过程中的力链变化,测量了力的传播速度演化速度(Clark et al.2015)如图 2 所示.
|
| 图 2 美国Duke大学开展的颗粒光弹实验采用高速相机观测到了力链结构的演化过程(Clark et al.2015) |
颗粒介质是多体相互作用体系,其结构特征是宏观均匀的,又具有若干颗粒大小的介观结构.结构表征是理解和认识颗粒介质物理问题的基础和依据,是当前制约颗粒介质力学和应用发展的瓶颈.随着研究的深入和结构分析实验手段的不断改进,人们发现颗粒体系的结构存在非均匀性,结构不均匀是认识动力学不均匀、力学性质转变等问题的基础.在颗粒介质的不同区域,差别很大几,乎每个颗粒周围的颗粒组成、体积分数、配位数和应力等与其他区域各不相同,这种局域差异影响、甚至决定颗粒介质的性质. 比如,在颗粒介质中出现一种局域的本征态,即某一位置颗粒的振动的振幅会随着距离很快衰减为零,即颗粒的振动态只存在于体系局域的范围内,这对颗粒介质的很多性能如输运特性等有重要影响.现有分析手段所描述的颗粒系统结构过于简单化,还不能建立结构和其宏观力学性能的对应关系,在基本理论和实验手段上的研究都非常困难.
人们将现有的最先进手段用于颗粒体系非均匀结构研究,并发现颗粒体系结构不均匀的确凿证据. 目前,径向分布函数和静态结构因子是物理学中表征结构信息最为常用的两个重要物理量,前者是在实空间,后者是在倒异空间,两者存在简单的傅里叶变换关系.径向分布函数和静态结构因子揭示了颗粒介质的结构特征.
3.1 径向分布函数径向分布函数g(r)(radial distribution function)仅考虑成对颗粒的相互作用,并假设各向同性且均匀.以任一颗粒为原点,颗粒介质中颗粒的分布仅与径向长度r 的大小有关g(r)定义为
|
(27) |
其物理意义是与原点颗粒相距r 处单位体积的颗粒数密度,从而获得颗粒介质的结构信息,△( rj-r)是狄拉克*函数. 式中ρ为系统的平均颗粒数密度,L为系统的尺寸,v为系统的体积. 由于颗粒介质的短程序,即有确定的最近邻及次近邻配位层. 如图 3 所示,颗粒介质g(r) 曲线有清晰的第一峰和第二峰并且第二峰,常出现分裂(比如,在
|
| 图 3 二维双分散颗粒系统在不同颗粒含量φ时的径向分布函数.插图是第一峰高度g1随φ的变化 |
静态结构因子定义为
|
(28) |
式中, k 是波矢量.由于采用了周期性边界条件,
|
| 图 4 不同的颗粒粒径分散度s, 周期边界条件下二维颗粒体系的静态结构因子S(k) 云图(张国华课题组提供) |
近年来,无序材料的力学和几何结构特性受到了较多关注. Silbert和Silbert(2009)对比分析了三维无摩擦和有摩擦颗粒体系的静态结构因子S(k),发现光滑颗粒体系在双对数坐标下 S(k) 曲线在低 k 区的线性行为,其斜率近似为1.Xu和Ching (2010) 研究了三维双分散光滑颗粒体系,发现静态结构因子强烈的依赖颗粒粒径比,以前发现的低 k部分的S(k)~k 仅仅是颗粒的粒径比接近1 时的特殊情况. 近年来,关于二维体系结构因子的研究也取得了很大进展.Meyer等 (2011) 发现在中间波矢区域二维聚合物溶液结构因子随 k呈幂律标度S(k)~kv,且随着链长度 N 的增大,幂律指数 v从 -2 变为 -3. Wen等 (2012)发现二维颗粒链的静态结构因子S(k)~(k-2),与密集的聚合物溶液结构相似.
张国华课题组研究了s=0,0.001,0.005,0.008,0.02,0.1,0.2,0.30.4和0.5时的二维颗粒体系的静态结构因子(冯旭等 2013)结果如图 5 所示,在高k区域尺寸分散度s< 0.1的s(k)曲线基本重合,并且在k' = 2kπ/L = 45,80,90附近各出现一个尖锐峰. 随着s增大,峰值趋于平缓而且在k'=80和100附近这两个峰逐渐合并成一个平滑的峰.在低k区域(k< 20),当s>0.1时,S(k)几乎不随k变化S(k)的值随着分散度的增大而增大; 当s ≤ 0.02时,二维体系的S(k)曲线在3<k´< 5区域遵从幂律标度,S(k)~k-4/3,如图 5 局部放大图所示.这一点符合二维线性聚合物链体系中在中间波矢范围的标度关系,S(k)~k-1/v,与文献(Peng et al.2010,Meyer et al.2011,Wen et al.2012)的结果类似.这暗示着二维单分散体系中存在颗粒链结构.图 5 插图还显示了三维单分散体系的静态结构因子在低k区域满足:S(k)~k, 与文献(Silbert & Silbert 2009)中的结论一致因此,可以判断,造成二维和三维体系S(k)在低k区域不同的原因可能是,在二维体系中更容易形成长程关联.
|
| 图 5 不同颗粒粒径分散度条件下的静态结构因子S(k)随k的变化.插图为二维和三维的单粒径颗粒体系的S(k)曲线低k部分数值拟合(冯旭等2013) |
先进的表征颗粒尺度的结构、性能及其结构和性能内在的相关性的实验技术,为直观表征颗粒介质局域有序结构提供有效手段,对颗粒介质的发展至关重要.实验技术上,现代结构分析手段主要依赖于同步辐射、中子散射等对颗粒介质结构的分析能力非常有限,很难同时实现高时间和空间分辨. 现有分析手段重构的三维颗粒结构只是基于一维信息,因此反演得到的颗粒介质结构变化不能被准确地探测到. 近几十年来不断有表征材料从颗粒尺度的结构和性能及其相关性的独特的实验技术被发展出来,比如,核磁共振技术(MRI)、μ -CT等可以测量颗粒介质结构特性.颗粒材料结构表征及其弛豫动力学研究方面,朱震刚和刘长松课题组采用力学谱技术(Chai et al.2014),研究了微颗粒振动能量耗散的振幅谱、深度谱和频率谱,得到了在微剪切振动下颗粒微观位形和力链结构转变的证据.孙其诚课题组利用光弹实验对二维颗粒物质中力链网络及应力传递进行了研究(Liu et al.2009,Sun et al.2015).张洁等利用光弹实验对2D颗粒剪切Jamming过程展开研究(Bi et al.2011).近期王宇杰课题组利用同步辐射X射线影像技术对颗粒材料体积弛豫动力学过程进行了研究(Xia et al.2015). 这些结构表征和性能研究的实验技术和计算技术,为系统地研究颗粒结构特征和力学性质提供了前所未有的实验条件、方法和研究基础.
颗粒介质的结构及其变化与模量密切关联, 弹性模量的变化可以敏感地反映颗粒介质结构的变化. 由于颗粒介质结构表征的困难, 因此, 弹性模量可以作为一个重要参数来研究和描述结构的变化和演化规律. 这是因为弹性模量包含着颗粒物质中颗粒间作用的最基本的信息, 颗粒介质中颗粒的各种运动包括振动、扩散、流变、弛豫都和颗粒之间的接触作用密切相关. 所以, 颗粒介质中颗粒的运动实际上是颗粒间接触性质的变化, 反映在宏观性质上就和弹性模量相关, 这是弹性模量与颗粒介质结构、特性和性能的关联性的物理原因. 颗粒物质的一个基本特征是其能够对外部机械扰动产生线性和非线性响应, 这使得声学测量成为研究颗粒体系弹性模量、内部结构和力学性质的基本手段(Zhang et al.2012,Merkel et al.2014,Zheng 2014)
4.1 声速测量根据经典弹性理论,各向同性弹性体中体积模量B和剪切模量G与纵波波速cL和横波波速cT的关系为
|
(29) |
其中,ρ是弹性介质的质量密度.Jiang和Liu(2009)发展了一个颗粒介质弹性能密度函数, 该函数是弹性应变张量的函数, 其对弹性应变分量的导数是对应的弹性应力分量, 在卸载时应变可自动恢复,-弹性应变关系是非线性的.弹性模量可以由颗粒介质弹性能密度函数推导得到,确定颗粒体系的弹性能表达式,是分析颗粒弹性的关键所在. Sun等结合GSH理论的弹性能模型,Jamming 点附近的弹性特性(Xu 2011,Liu & Nagel 2010),提出了新的弹性能模型(Sun et al.2015)
|
(30) |
其中,B0 是弹性系数,φc是Jamming点处的临界颗粒体积分数,φ-φc反映了偏离Jamming点的程度. 
|
(31) |
采用声速测试技术,可以确定式(30)中弹性参数,如图 6 所示. 对于玻璃珠组成的体系而言,B0 = 9.22×109 Pa,φc≈0.635,a=0.50,ξ=3.88,b=0.39.
|
| 图 6 不同压强下的波速. 实验测量结果(Domenico 1977,Makse et al.2004)与基于式(31)的计算结果对比 |
Jia等 (1999)和Lherminier等 (2014)发现玻璃珠样品的声速随压强的变化满足分段幂律规律,即压强较低(30~120 kPa)时满足:
图 7 是北京科技大学张国华课题组采用飞行时间法测得的玻璃珠样品的横波声速和纵波声速随压强的变化曲线(张攀等 2016). 可以看出,同一压强下,纵波声速(cL)明显大于横波声速(cT),且cL和cT随压强变化均呈现幂律标度,即cT ~p0.280 9cL ~p 0.381 7. 根据Hertz 接触理论颗粒体系中横波声速和纵波声速随外加压强p变化均呈现幂律规律:c~p1/6. 显然,实验测得cT和cL随压强的幂率标度指数与Hertz 接触理论预言存在差异.导致上述差异的因素可能有以下两个:(1) 理论模型要求颗粒是光滑的,但是实际颗粒的表面是粗糙的,这会使实际结果同Hertz 接触理论预言存在一定偏差;(2) 纵波对于应力场的各向异性更加敏感,而横波对于结构各向异性更加敏感(Khidas & Jia 2010),
这也会导致两随压强变化出现不同的幂率标度. 图 7的插图显示了纵波与横波声速比cL/cT随压强p的变化,以看出cL/cT大约为1.6近似是一个常数.Lherminier等 (2014)在研究声波通过受压的二维颗粒介质时也发现类似的现象.
|
| 图 7 声速和G/B的测量(a)横波声速和纵波声速随压强的变化实线为拟合曲线插图是声速纵横比随压强的变化(b)G/B随压强的变化(张攀等 2016) |
由式(29)可以得到
|
(32) |
图 7 显示了式(32)得到的玻璃珠样品的G/B随压强p的变化规律发现实验给出的G/B值在0.58 ~1.39之间且颗粒体系G/B随压强p的增大而减小,满足幂律关系: G/B~p-0.453 9,如图 7(b)中红色直线拟合所表示. 根据(O'Hern et al.2003,Tighe 2011,Basu et al.2014,Goodrich et al.2014)的结果,对于Jammed颗粒体系的G/B~(φ-φ j01/2~p1/3(其中,φj0对应J点的体积分数). 显然我们实验上发现的G/B随压强p的变化规律(G/B~p-0.453 9)与Hertz 接触理论不相符暗示连续弹性介质理论在低压强下可能不再成立.
最近Wang等 (2015)的研究表明在卸载过程中硬球胶体玻璃分别经历高压强的TL玻璃态(横模纵模都存在的状态)到低压强的L玻璃态(只存在纵模的状态)及从L玻璃态到Unjammed态(横模和纵模都消失的状态)的两次转变.并且卸载过程中硬球胶体玻璃的剪切模量与体积模量之比G/B随压强p变化呈分段幂率规律当p>pj (pj是相转变点J点)时,G/B~p1/3(与赫兹接触理论的结果相同); 当pg < p < pj(pg 表示玻璃化转化点)时,G/B~p-1/2(与赫兹接触理论预言的G/B~p1/3不同).
4.3 声衰减系数随压强的变化声波在颗粒介质中传播时总会伴随着不同程度的损耗,也就是随着传播距离增加而减弱的现象.弹性波在颗粒体系中的衰减一般有三部分组成:和颗粒介质密度有关的散射衰减、和声源方向有关的扩散衰减和介质粘滞摩擦有关的吸收衰减.在实验中,所测量的颗粒介质体系沿着声源方向的长度较短,因此一般不考虑扩散衰减. 为了获得体系的声衰减系数α,实验上同时测量初始激励的强度I0(正比于加速度振幅平方)及声波通过长度为x的颗粒样品后得到强度I(x),进而得到体系的声衰减系数(Vitelli 2010)
|
(33) |
由式(33)可知,声波衰减系数越大,声吸收越多,反之则越少.图 8 显示了不同压强下玻璃珠样品中纵波的声衰减系数α. 实验测得的衰减系数值在0.16~0.25之间,且随着颗粒体系压强p的增大幂率减小: α~p-0.187 9,如图 8 中的红色直线所示.由于声衰减系数的值是通过频谱法测定的,颗粒体系几何带来的影响在数据处理过程中已经被去除(Brilliantov et al.1996).
|
| 图 8 纵波声衰减系数随压强的变化(实线为幂率拟合曲线) |
玻璃珠样品的纵波衰减系数α随压强增大而减小. 人们一般认为,玻璃珠样品的纵波衰减系数α随压强的变化可能与系统内部结构随压强调整有关.即随着体系压强的增大,颗粒之间接触更加紧密,对应的平均接触数增加,颗粒体系内部结构更不容易发生调整,相应的声波耗散减小. 此外,刘长松课题组(Wang et al.2007,2008a)研究颗粒体系衰减耗散问题时,都发现声衰减与内部接触变化有关.
注意到,如果只考虑散射衰减,三维Hertz 接触球形体系的衰减系数随压强变化呈幂律关系: α~p-7/3(Vitelli 2010).实验得到的玻璃珠样品的纵波衰减系数随压强变化的幂率指数为-0.187 9,与-7/3相比差别特别大.造成实验得到的幂指数和文献(Vitelli 2010)预言幂指数差异很大的原因可能是声波在颗粒体系中的波长(λ= c/f≈0.2~0.5 m)远大于颗粒尺寸,体系中散射衰减非常小,颗粒介质本身的吸收衰减是主要机理.
4.4 二阶谐波振幅随压强p的变化由于系统非线性的影响,颗粒介质中声波除了一阶声波外,还出现二次甚至三阶的谐波. 非线性系数β是表征颗粒介质非线性强弱的重要参数. 通常情况下,非线性系数越大,系统非线性导致的波形畸变越严重. 根据非线性声学,体系的非线性系数(β)和二倍频振幅(μ2w)与基频振幅(μ1w)平方的比值有关(Brunet et al.2008)
|
(34) |
其中,μ2w代表接收二倍频振幅,μ1w代表接收基频振幅,c代表声速w代表频率. 为了研究非线性系数和压强的依赖关系,我们利用傅立叶变换法分析了玻璃珠样品的μ2w/μ1w2随压强的变化曲线,如图 9 所示.μ2w/μ1w2随压强的增大而幂率减小: 
|
| 图 9 二倍频振幅与基频振幅平方的比值 μ2w/μ1w2随压强的变化关系图(实线为拟合曲线) |
非线性现象还是研究弱力链规律的有效手段. 就现在而言, 各向异性的介质研究非线性声学主要以模拟为主. 这是非常狭隘的, 故开展新的理论和实验就非常的迫切和需要. 非线性效应对于弱力链的出现是非常敏感的. 一方面, 弱力链分数主要受到材料的应变的影响而变化; 但是另一方面, 这些接触在线性弹性力学中很难被阐明. 因此, 像颗粒填充物重排的早期阶段引起的弱应变波很难通过线性声学传播来评估. 相比之下, 非线性效应像高次谐波的产生表现出了对结构更高的敏感性. 非线性声学的研究固体结构受到越来越多的关注. 固体材料受力之后, 胡克定律几乎完美的描述了材料中的应力与应变(单位变形量) 之间成线性关系, 而非线性声学就是处理应力应变的偏差. 这个结果打破了弹性波的叠加效应, 使得非线性引起的波频谱分量水平增加可以作为这样"弱特征" 的存在的敏感指标, 因而可以用于颗粒材料弱力链网络的检测.
5 颗粒介质失稳机理的研究作为一种类似临界现象的复杂物理过程, 失稳研究是今后颗粒介质力学研究的重点之一, 因为颗粒介质体系的稳定性是其应用的基础. 颗粒介质的失稳研究重点主要涉及的是颗粒介质在应力场中的破坏行为. 颗粒介质体系失稳更为广泛地存在于自然界中, 如泥石流、碎屑流等灾害的起动, 以及塌方等, 因此颗粒介质失稳机制研究的突破有助于对常见地质灾害的认知和防治.
5.1 Jamming相图Jamming转变是指无序材料,包括颗粒介质、玻璃、泡沫和胶体悬浊液等,Unjammed态与Jammed态之间的非平衡转变,对于颗粒体系而言则对应类流体与类固体之间的转变,由Liu和Nagel首次于1998年首次提出,之后逐渐完善而成为刻画颗粒系统等无序体系状态转变研究和失稳研究的基本框架(Liu & Nagel 1998).
如图 10 所示,温度T、体积分数φ和剪切应力Σ为控制参数,与气体的温度T、体积V和压强P类似. 目前,物理学家分别研究了某一个(或两个)参数变化时,系统发生Jamming转变的规律. 比如,玻璃转变的正常相图位于(1/φ)-T平面,转变线在(1/φ)-T面上把体系分为Jammed态(如玻璃)和Unjammed态(如液体);而泡沫、乳胶或颗粒物质的相图位于(1/φ)-Σ平面以体积分数φ为变量的临界屈服应力Σ(φ)曲线将相图分为Jammed和Unjammed态,对于颗粒物质来说对应着颗粒固体(granular solid)和颗粒液体(granular fluid). 另外,在T-Σ面内,改变T或Σ时系统也会发生Jamming转变,目前还没有找到与之对应的实际系统.
|
| 图 10 颗粒物质等无序材料的Jamming 相图 |
从图 10 所示的Jamming相图可以看出,Jamming相转变特性依赖于达到相边界的路径,比如沿着φ轴,当各向同性的jammed颗粒体系φ值降低到某一点时,体积弹性模量和剪切模量为零,这个转变点称为J点(Jamming point),所对应的临界体积分数为φc,△φ=φ-φc定义为到J点的距离. 体积分数φ是影响颗粒介质jamming 转变最为主要的物理量,是颗粒总体积与体系体积之比,或者局部体积分数定义为包含单颗粒的Voronoi 多面体中,该颗粒体积与Voronoi 多面体体积之比. 人们发现,颗粒体系内局部体积分数符合Weibull分布,且颗粒体系越密集,局部体积分数越趋于均匀分布. 颗粒系统在J点表现出某些临界行为,在某种程度上接近于二级转变,比如体积弹性模量呈现△φ的幂律标度,并且存在随△φ发散的特征长度等.对于有限尺寸的体系,边界对φc的影响不可忽略,使得φc在一定范围内分布; 对于无限大尺寸的体系,φc是常数. 另外,光滑的刚性颗粒的形变恒为零,体系总处于J点亦即△φ=0.
无摩擦软球体系是研究Jamming转变的理想系统,其力学和几何学性质在J点附近随着△φ呈现奇特的幂律标度.对于由N个d维光滑弹性颗粒组成的处于静力平衡的颗粒体系,假设总接触力个数(总自由度)为NZ/2,力平衡方程数目(总约束数目)为Nd. 显然,当Z ≥ 2d时,接触力有解,Ziso 2d对应等静态(isostatic state)值,即此时体系中每个颗粒的平均配位数等于体系稳定时所需要的最小值.在J点,体系压强为零,颗粒没有形变jammed条件要求颗粒间距离等于它们的半径之和(刚好接触条件),体系有Nd个位置自由度和NZ/2个位置约束条件. 显然,当Z ≤ 2d时,体系有唯一解,Zc=2d对应J点的临界配位数.由上述分析可知,无摩擦软球体系在J点达到等静态,Zc=Ziso=2d. 在J点,Z具有明确数值而且在J点附近表现出奇特的幂律规律(Song et al.2008).考虑到颗粒体系的力学特性敏感地依赖于Z,可以预期弹性模量等在J点附近也表现出类似的标度行为.光滑颗粒组成的体系体积模量B和剪切模量G依赖于颗粒间作用势U的类型. 一般而言,有3种作用势:Hertizan作用势U~δ5/2,Harmonic作用势U~δ2,一般幂律作用势,U~δα.研究结果发现颗粒体系的G~(△φ)α-3/2且对于Hertizan作用势和Harmonic作用势(李广信 2006),
2011年,Ciamarra等 (2011)提出了有摩擦颗粒物质的Jamming相图,如图 11 所示,密度倒数、剪切应力和摩擦系数3个参数控制着有摩擦颗粒的Jamming转变.与Liu-Nagel的Jamming相图不同,摩擦的引入为相关控制参数,表征穿过不同Jamming转变线的结构变化，并且考虑在有限剪切应力下Jammed的颗粒体系表现出易碎(fragile)行为.在这个相图里,φJ1(σ,μ),φJ2(σ,μ)和φJ3(σ,μ)面包围不同的流动特性的区域.在这个相图里,假定Jamming表面沿着一个界限分明的Jamming转变线φJ(σ)汇聚于μ =0面,与在无摩擦的体系里显示的"flow and jam"现象结果一致.
|
| 图 11 有摩擦颗粒体系的Jamming相图(Ciamarra et al.2011) |
在外界扰动条件下,颗粒体系的动力学行为呈现明显的不均匀性,这取决于介观尺度上的结构的复杂性,比如自由颗粒(rattler)、微剪切带、漩涡、力链以及剪切带等,它们出现在不同的宏观力学阶段.自由颗粒是指与周围颗粒的接触点数目小于等静态配位数的颗粒,对于体系的力学稳定性贡献很小,小的扰动都会导致自由颗粒重排以及结构破坏(McNamara et al.2008);微剪切带在应变硬化阶段起着决定性作用(Kuhn 1999);漩涡持续时间非常短暂,可能在剪切带形成过程中起到很重要的作用;剪切带是颗粒体系破坏的前兆,吸引了很多学者对它进行研究(Utter & Behringer 2004). 动力学不均匀性的刻画及其根源的研究是重点之一.
颗粒材料等无序固体的另一个典型特征是能够对其外部扰动产生强烈的非仿射响应.当对晶体材料施加变形时,它趋向于仿射地遵循宏观应变.而对颗粒材料施加变形时,由于缺乏局部反转对称性,每个颗粒除了产生一个仿射应变之外,还产生一个附加的非仿射位移.最近,Zaccone和Terentjev(2014)的研究表明非晶体材料(例如玻璃合金)的体积模量表现得几乎与仿射变形的假设一致(即压缩响应几乎是仿射的),而剪切模量显著减少(比对应晶体低50%),暗示剪切变形导致无序体系中出现更多非仿射位移. 认为,原子间的体积关联导致了局部短程序,此局部短程序显著地减少了各向均匀压缩下的非仿射位移,因而导致了非晶体材料中体积模量的准仿射行为,这暗示了无序材料中剪切模式比压缩模式更软的事实可能与在振动态密度反常玻色峰中剪切模式占优势存在关联,泊松比可能在确定无序体系的易碎性中扮演重要角色.
宏观平均化的应变定义掩盖了颗粒尺度的运动:在大尺度上观察到变形是仿射的. 在观察局部化现象时,必须定义基于颗粒尺度的应变,亦即将应变定义在单颗粒和它周围一圈的颗粒上(Tordesillas et al.2008). 图 12显示了Delaunay三角剖分方法,将空间剖分成包围单颗粒的凸多边形,边界用S表示,面积用V表示,参考颗粒圆心到第一圈上的第c个颗粒圆心的距离矢量为lc.
|
| 图 12 变形前后的Delaunay三角剖分 |
定义第c个颗粒和第c+1个颗粒之间连接线的单位外法线向量 n为
|
(35) |
当体系受力变形后,保持之前的接触关系不变. 定义u为中间参考颗粒的位移向量,uc为参考颗粒第一圈上第c个颗粒的位移向量.则第c个颗粒相对于中间参考颗粒的位移向量 Pc表示为
|
(36) |
其中ω为参考颗粒的转角.单颗粒的局部应变定义为
|
(37) |
求和是对参考颗粒外面第一圈颗粒B逆时针方向进行的.基于颗粒的局部应变的平均值与宏观应变是完全吻合的.由于单颗粒的实际位移,与通过局部应变和颗粒接触方向的乘积求得的位移不等,因此,定义该差值为非仿射位移
|
(38) |
那么,参考颗粒的非仿射体应变为
|
(39) |
局部的仿射应变与宏观应变存在很好的对应关系,而局部的非仿射应变对应着变形的局部化现象,且体系的变形越不均匀,得到的局部非仿射应变越大,即非仿射应变率不仅可以刻画变形局部化,而且可以直接对应局部化变形的大小.
图 13 所示为二维颗粒体系在双轴压缩条件下,局部非仿射应变空间分布的演化过程,可以看出在初始阶段(如ε1 = 1.7%)非仿射应变大的颗粒分布比较分散,之后渐渐的发展成较为集中的局部化带状结构,这一阶段对应着土力学试验中微剪切带逐渐发展成为主剪切带的阶段.在轴应变达到ε1= 3.4% 时,体系中出现 "X"形剪切带并且在之后一直保持着.
|
| 图 13 二维颗粒体系在双轴压缩时不同轴向应变时的局部非仿射应变较大颗粒的空间分布 |
类似于液体,固体在切应力或切应变下可以流动.
晶体是通过其内部晶格缺陷来控制其重排流动,而在颗粒材料中颗粒的重排倾向于局部化,但是对于其中的结构缺陷非常难以识别.晶体流动由晶格缺陷演变引起,而颗粒体系流动通过局域化重排引起,最近提出的"软点"(soft spots)物理概念用来确定可能发生重排的颗粒.颗粒物质是无序系统,结构上与无序胶体以及液体类似,可以构造颗粒体系Hessian矩阵分析颗粒物质振动态密度.Hessian矩阵H是颗粒接触势能函数关于所有颗粒位置坐标变量的二阶偏导函数,在瞬态稳定位形中的取值按顺序排列而组成,其矩阵元
|
| 图 14 软点与颗粒重排图(深色点是软点插图箭头显示了颗粒非仿射位移)(Dong et al.2015) |
Cubuk等 (2015)和Dong等 (2015)成功地利用Manning方法确定在二维颗粒柱实验、三维Lennard-Jones玻璃以及Hertzian接触的颗粒系统中对重排敏感的流动缺陷颗粒，暗示Manning和Liu(2011)提出的模式分析技术是预言和识别颗粒材料等无序材料中对动力学不均匀性负责的介观结构的有利工具.此外Schoenholz等 (2014)跟踪不同温度下在一个二维剪切的热Lennard-Jones玻璃中单个软点随时间的演化，发现单个软点的寿命与结构弛豫时间尺度相当，且大部份软点能够持续经历许多重排，暗示软点是体系非仿射重排的很好的预言子，为基于单个软点动力学行为的颗粒材料等无定形固体的介观塑性理论铺平了道路.
局部塑性事件的空间结构分析是理解颗粒体系塑性行为的关键，不同的构成形式可能导致不同的塑性行为.最初的失效模式的空间结构分析为理解基本的塑性事件的形成提供了一种有效的方法.为了更直观地表达体系颗粒重排与重排前体系低频模式的相关性，图 15(a)显示了塑性发生前后颗粒归一化的非仿射位移矢量图和塑性发生前体系最低频模式的极化矢量图(图 15(b)).
|
| 图 15 (a)塑性事件发生后归一化的非仿射位移图,(b)塑性事件发生前最低频率模式的极化矢量图,(c)投影系数α2随模式数的变化(这些本征模式由发生颗粒重排之前的体系计算而得,颗粒重排矢量为颗粒重排前后的位移)(插图: 投影系数的积累通过计算低于模式数的贡献总和而得(Dong et al.2015) |
极化矢量模式和颗粒重排相关性可以被量化为本征矢量eωV(由2N个分矢量组成的一个完备集)与非仿射位移矢量uV(包含所有颗粒的非仿射位移分量)的归一化内积,这里N是体系颗粒的数目(不包括rattlers颗粒),投影系数定义为
|
(40) |
图 15(a)和图 15(b)显示了在剪切作用下二维颗粒体系极化矢量模式和颗粒重排相关性,可以看到低频模式(特别是最低几个模式)对非仿射位移的贡献明显较大.插图是投影系数的积累,
颗粒材料的一个典型特征是对外部扰动产生强烈的非仿射响应(亦即颗粒的无序运动).准静态颗粒介质中颗粒无序运动的持续时间极短(微秒级),产生的位移极小(微米级),常用的粒子图像测速(particle image velocimetry,PIV)或磁共振测速(magnetic resonance velocimetry,MRV)等都无法准确测定颗粒速度. 近年来,测量技术取得较大发展,Deen等 (2004)采用PIV技术测量了振动流化床的颗粒运动验证了数值模拟的结果.Chou等 (2009)采用PTV技术测量了滚筒中顶部流动层的颗粒速度分布.由于图像测速技术只能用于二维流化床,Wildman等 (2001)用正电子发射颗粒跟踪(positron emission particle tracking,PEPT)技术跟踪计算了三维振动流化床内颗粒的速度波动,但该方法存在统计颗粒样本不足的问题.Wang等 (2008)采用PTV技术跟踪测量了低速剪切室中的颗粒运动轨迹,发现颗粒的均方位移(mean square displacementMSD)与平均位移线性相关,其斜率即有效温度,但该方法存在统计颗粒样本不足的问题,但这些结果还远没有达到人们的预期. 近期,散斑能见度光谱法(speckle visibility spectroscopy,SVS)是研究颗粒材料微观动力学过程的新方法(Bandyopadhyay et al.2005),利用散斑图像具备时空遍历一致性通过计算散斑图像的空间自相关,计算出颗粒速度的脉动(Yang et al.2015,Li et al.2016),见图 16.目前常用测量颗粒运动方法如表 2 所示.
|
| 图 16 (a)散斑能见度光谱技术测量转筒中颗粒的脉动速度装置和(b)脉动速度二阶矩的空间分布(亦即颗粒温度)(Li et al.2016) |
颗粒介质力学的研究对象是多体相互作用的无序体系,很多基本问题都未解决,其基本范式还没有建立,系统理论框架缺失也严重影响了颗粒材料的应用和新发展.此类复杂体系的若干科学问题之间有紧密联系,必须对整个颗粒介质系统进行分析,即使这样较为粗糙的研究也是必要的.
颗粒介质广泛分布于工程建设和自然界中,其基础研究需要与其他学科最新的成果结合,建立新的概念和范式去认识颗粒力学的基本问题. 与此同时,颗粒体系研究的应用是以新技术和新工艺的发展、工程安全性评估以及地质灾害预测等等重大需求为背景的.对颗粒介质体系的科学研究需要交叉、融合以及不同学科参与对建立完整的颗粒介质体系理论,推动颗粒介质力学的发展极为重要.
目前,颗粒材料的力学研究和工程应用都进入一个关键的阶段,迫切需要梳理下一步研究和发展的思路,凝炼重要和关键的科学问题和技术难题.下面试图总结、归纳出几个核心问题,希望能起到参考作用与国内外同行共同加强颗粒介质研究.
(1) 到目前为止绝大多数离散单元法数值模拟仍处于定性或半定量阶段,造成这一现状的主要原因包括:(a)模型的简化.所采用的颗粒接触模型都是从Hertz-Mindlin理论简化而来,很少考虑多体接触时的应力耦合效应;(b)模拟参数难以确定.接触模型的参数仅与颗粒属性有关,然而由于某些模拟参数如转动摩擦系数和转动弹簧刚度系数等很难通过实验直接测量得到,且受限于计算能力,多数情况下难以直接采用颗粒的真实材料参数,使得参数选择带有一定的随意性;(c)颗粒形态多样性.实际颗粒表面都带有一定粗糙度,且形状一般都偏离标准球形而离散单元法中通常将颗粒处理为理想的光滑球体.即使某些离散单元法数值模拟采用粘结颗粒或者多边形颗粒来构造非规则形状颗粒,但也只能采用有限的几种形状,不可能完全复现实际颗粒系统中颗粒形状的多样性.这些假设和简化都会影响模拟结果的定量准确性. 因此,不应刻意追求离散单元法数值模拟结果与实验数据的准确性.发展精确的接触模型、开发高效准确的并行算法、注重与其他数值模拟方法如格子玻尔兹曼、连续介质模型等耦合可能是今后发展离散单元法需要努力的方向.
(2) 颗粒体系进入非弹性阶段后,屈服准则和流动法则等的机理认识非常肤浅,这主要是由于体系内部形态不均匀,并会随着塑性的发展出现应变局部化等现象,非连续介质的特点尤为显著.过度的平均化处理(homogenization)把颗粒处理成均质的、连续的介质掩盖了这些内部结构(包括几何结构和力结构)的各向异性.初始时刻平均性质大体相同的颗粒体系的峰值强度呈现显著的分散性,这是因为颗粒试样的介观力位形和结构位形在初始时刻就有差异作为局部非平衡因素,这些差异在介观结构演化过程中得以放大,导致宏观力学性质差异较大. 从介观结构入手,建立颗粒接触力-宏观应力、颗粒位移-宏观应变的关系掌握介观应力分布、变形机理及破坏规律,可以为工程中出现的问题给予更科学的解释,并给出相应的有效解决办法,是现阶段颗粒物质介观研究的主要课题.
(3) 无序固体结构分析的概念和方法影响了现有颗粒介质的结构研究,但是颗粒介质的结构复杂、动态、和跨尺度,对其描述需要全新的概念、方法和思路,从颗粒相互作用的角度来描述颗粒介质结构是颗粒介质研究的挑战. 因此,如何表征颗粒介质的无序结构特征是颗粒介质基础研究和工程应用的核心问题之一.如何确立颗粒结构时空演化基本特征对认识结构很重要,并且颗粒介质结构不均匀性和动力学不均匀的联系等问题都需要回答.
(4) 颗粒物质在介观尺度上的结构是不同而且复杂的,由于具体问题和研究者视野的差异,介观力学的实体和表征可能不同,可以是由相邻颗粒质心、颗粒间接触点连接而成、包含一个孔隙的最小封闭空间,能够在承受外载时保持静定的空间结构,也可以是其他介观单元.对于颗粒体系,介观结构及其运动具有或然性,引入统计方法至关重要,必须理解在非平衡演化过程中低层次的物理本质如何传递到更高的层次,与更高层次强耦合,以及如何用公式表示对更高层次的影响.我们认为它应该部分根植于实验,部分根植于对力学规律的深入分析.这就有必要开展颗粒物质的统计力学研究.有必要在完善上述介观结构研究基础上,主要分析介观结构某一特征物理量的空间分布在剪切带发生和发展中的演变,构建颗粒物质的介观力学,其主要任务就是定性和定量分析它们对宏观力学现象的影响,预言颗粒体系的应力-应变关系、流变关系以及理解局部化问题,亦即颗粒物质的介观力学注重对宏观现象及运动规律做出介观解释,得到剪切带形成的失稳判据,预测剪切带特征宽度,揭示剪切带的形成和演化机制等,从介观结构层次全面解释颗粒体系力学性质.
(5) 颗粒介质的失稳是一种类似临界现象的复杂物理过程,和颗粒介质体系的结构,动力学和热力学密切相关,具有不同的时-空尺度宏观失稳是从介观局域失稳聚集和发展而成,从局域失稳发展成为宏观失稳的过程是认识颗粒介质失稳的关键问题.目前,介观尺度局域化的变形在实验上还无法直接观察,由局域至宏观的失稳过程更不清楚,能否找到影响颗粒介质失稳的局域或整体结构特征一直是颗粒介质体系研究中的一个核心难题.该问题的解决将增进人们对地质灾害起动机理的认识发展更加有效的防治措施.
(6) 发生失稳流动的颗粒体系其黏性系数与剪切速率及力链长度的规律研究是动力学研究的重点课题.目前较为成功的是局域流变模型(local rheology),比如法国μ(I)流变. 对于刻画准静态颗粒流,则有问题. 此时颗粒之间密切接触,使得速度脉动和力链网络都存在大尺度的空间关联,应力紊动产生的弹性波在颗粒之间的传播长度,由快速流态下的颗粒粒径d,转变为了力链长度l因而颗粒介质内部流动区域产生的扰动,沿力链传递到远处区域,颗粒介质表现出明显的非局域流变(non-local rheology)特性影响了动量输运.这一特殊的流变现象可以通过在局域模型中叠加考虑包含非局域特征长度影响的"二级流变"(secondary rheology)来实现,引入特征长度构建的非局域流变模型,或通过对各流层中的颗粒紊动进行积分处理引入颗粒之间的长程作用.尽管一些模型对准静态颗粒流的刻画上取得了成功,可能抓住了问题的关键,但对"二级流变"的机理探讨还存在很多争议,需要从介观结构层面予以澄清.
(7) 发展声学测量技术,用线性声学手段测量体系中横波、纵波的声速,进而获得体系强力链骨架的剪切、压缩模量,利用非线性声学手段获得力空间弱力链的分布,通过体系声衰减谱测量获得体系散射平均自由程相关的介观尺度的特征长度,研究该特征长度随力学载荷的演化规律.
已有的颗粒介质力学研究已为进一步探索研究颗粒体系、为最终建立基本理论框架积累了大量丰富经验和基础.颗粒介质中很多基本力学问题都还悬而未解,这是因为它是无序多体相互作用体系,其基本范式和理论框架还没有建立,严重影响了颗粒介质的力学研究,工程应用的新发展. 从另一方面讲这也为后续研究提供了一个良好机遇.颗粒介质的基础研究需要和其他学科最新的成果结合,建立新的概念和范式从全新的角度、思路、理念去认识颗粒介质的基本问题,这样才能把颗粒介质研究从定性、唯象的研究中解脱出来. 在过去的10年间,我国颗粒介质的研究队伍不断壮大,基础研究逐渐接近国际主流.2016年暑期,我国组织了5次颗粒介质的国际会议,促进了我国颗粒力学的发展.
颗粒介质研究的生命力还取决于它与不同学科的结合与交叉.不同学科的相互影响和启发,可能导致颗粒介质研究快速发展.目前颗粒介质研究与其他学科的交叉远远不够,颗粒介质很多基本问题的解决需要物理、力学、数学和工程等多领域科学家和工程师共同合作.颗粒介质普遍应用于工程建设、工业生产和日常生活中除了纯粹的物理学探索外,颗粒介质的研究最好紧密结合实践应用在解决工程难题的同时,颗粒介质的研究也会有的放矢,更具生命力.我们熟悉的有,核能领域中相关颗粒流的精准控制技术,水利水电工程领域中堆石坝变形预测,地质灾害领域碎屑流起动机理和流动规律.我们相信还有更多的工程应用中的核心技术需要颗粒介质的基础研究.因此,发展适用于实际情况下的颗粒材料研究是当前工程领域的迫切需求,进一步的研究有助于推动我国相关领域的发展并对社会经济发展做出贡献.
| 冯旭, 张国华, 孙其诚. 2013. 颗粒尺寸分散度对颗粒体系力学和几何结构特性的影响. 物理学报, 62:184501. ( Feng X, Zhang G H, Sun Q C. 2013. Effects of size polydispersity on mechanical and geometrical properties of granular system. Acta Physica Sinica , 62:184501. ) |
| 黄文熙.1983. 土的工程性质. 北京: 水利电力出版社 . ( Huang W X.1983. Engineering Properties of Soil. Beijing: China Water & Power Press ). |
| 李广信. 2006. 土的清华弹塑模型及其发展. 岩土工程学报, 28:1–10. ( Li G X. 2006. Characteristics and development of Tsinghua elasto-plastic model for soil. Chinese Journal of Geotechnical Engineering , 28:1–10. ) |
| 李家春.2014. 中国学科发展战略: 流体动力学. 北京: 科学出版社 . ( Li J C.2014. Development Strategy of China Discipline: Fluid Dynamics. Beijing: Science Press ). |
| 沈珠江.2000. 理论土力学. 北京: 中国水利水电出版社 . ( Shen J Z.2000. Theoretical Soil Mechanics. Beijing: China Water & Power Press ). |
| 孙其诚, 厚美瑛, 金峰, 等.2011. 颗粒物质物理与力学. 北京: 科学出版社 . ( Sun Q C, Hou M Y, Jin F, et al.2011. Physics and Mechanics of Granular Materials. Beijing: Science Press ). |
| 孙其诚, 王光谦.2009. 颗粒介质力学导论. 北京: 科学出版社 . ( Sun Q C, Wang G Q.2009. Introduction to Granular Mechanics. Beijing: Science Press ). |
| 孙其诚, 王光谦. 2008. 颗粒流动力学及其离散模型评述. 力学进展, 38:87–100. ( Sun Q C, Wang G Q. 2008. Review on granular flow dynamics and its discrete element method. Advances In Mechanics , 38:87–100. ) |
| 汪卫华. 2013. 非晶态物质的本质和特性. 物理学进展, 33:177–351. ( Wang W H. 2013. The nature and properties of amorphous matter. Progress in Physics , 33:177–351. ) |
| 张攀, 赵雪丹, 张国华, 等. 2016. 垂直载荷下颗粒物质的声波探测和非线性响应. 物理学报, 65:024501. ( Zhang P, Zhao X D, Zhang G H, et al. 2016. Acoustic detection and nonlinear response of granular materials under vertical vibrations. Acta Physica Sinica , 65:024501. ) |
| Ai J, Chen J F, Rotter J R, Ooi J Y. 2011. Assessment of rolling resistance models in discrete element simulations. Powder Technology , 206:269–282. doi:10.1016/j.powtec.2010.09.030 |
| Aoki K M, Akiyama T. 1995. Simulation studies of pressure and density wave propagations in vertically vibrated beds of granules. Physical Review E , 52:3288–3291. doi:10.1103/PhysRevE.52.3288 |
| Basu A, Xu Y, Still T, Arratia P E, Zhang Z, Nordstrom K N, Rieser J M, Gollub J P, Durian D J, Yodh A G. 2014. Rheology of soft colloids across the onset of rigidity: scaling behavior, thermal, and non-thermal responses. Soft Matter , 10:3027–3035. doi:10.1039/c3sm52454j |
| Bandyopadhyay R, Gittings A S, Suh S S, Dixon P K, Durian D J. 2005. Speckle-visibility spectroscopy: a tool to study time-varying dynamics. Review of Scientific Instruments , 76:093110. doi:10.1063/1.2037987 |
| Bardet J P, Huang Q. Rotational stiffness of cylindrical particle contacts//Proceedings of the Second In-ternational Conference on Micromechanics of Granular Media, Birmingham, UK. Thornton C ed., A. A. Balkema: 39-44. |
| Baxter J, Tüzün U, Burnell J, Heyes D M. 1997. Granular dynamics simulation of two-dimensional heap formation. Physical Review E , 55:3546–3554. doi:10.1103/PhysRevE.55.3546 |
| Bi D, Zhang J, Chakraborty B, Behringer R P. 2011. Jamming by shear. Nature , 480:355–358. doi:10.1038/nature10667 |
| Brilliantov N V, Spahn F, Hertzsch J, Poeschel T. 1996. Model for collisions in granular gases. Physical Review E , 53:5382–5392. doi:10.1103/PhysRevE.53.5382 |
| Brunet T, Jia X, Johnson P A. 2008. Transitional nonlinear elastic behaviour in dense granular media. Geophysical Research Letters , 35:L19308. doi:10.1029/2008GL035264 |
| Chai L, Wu X, Liu C S. 2014. A universal scaling law of grain chain elasticity under pressure revealed by a simple force vibration method. Soft Matter , 10:6614–6618. doi:10.1039/C4SM00727A |
| Chen K, ManningML, Yunker P J, EllenbroekWG, Zhang Z, Liu A J, Yodh A G. 2011. Measurement of cor-relations between low-frequency vibrational modes and particle rearrangements in quasi-two-dimensional colloidal glasses. Physical Review Letters , 107:108301. doi:10.1103/PhysRevLett.107.108301 |
| Cheng Y Q, Ma E. 2011. Atomic-level structure and structure-property relationship in metallic glasses. Progress in Materials Science , 56:379–473. doi:10.1016/j.pmatsci.2010.12.002 |
| Chou H T, Lee C F. 2009. Cross-sectional and axial flow characteristics of dry granular material in rotating drums. Granular Matter , 11:13–32. doi:10.1007/s10035-008-0118-y |
| Ciamarra M P, Pastore R, Nicodemi M, Coniglio A. 2011. Jamming phase diagram for frictional particles. Physical Review E , 84:041308. doi:10.1103/PhysRevE.84.041308 |
| Chialvo S, Sun J, Sundaresan S. 2012. Bridging the rheology of granular flows in three regimes. Physical Review E , 85:0121305. doi:10.1103/PhysRevB.85.121305 |
| Clark A H, Petersen A J, Kondic L, Behringer R P. 2015. Nonlinear force propagation during granular impact. Physical Review Letters , 114:144502. doi:10.1103/PhysRevLett.114.144502 |
| Cooke M H, Bridgwater J. 1982. The simulation of a particle disperse. Powder Technology , 33:239–247. doi:10.1016/0032-5910(82)85062-6 |
| Cubuk E D, Schoenholz S S, Rieser J M, Malone B D, Rottler J, Durian D J, Kaxiras E, Liu A J. 2015. Identifying structural flow defects in disordered solids using machine-learning methods. Physical Review Letters , 114:108001. doi:10.1103/PhysRevLett.114.108001 |
| Cundall P A. 1971. A Computer model for simulating progressive large scale movements in blocky rock systems//Proceedings of the Symposium of the International Society of Rock Mechanics, Nancy, France. 1: 132-150. |
| Cundall P A, Strack O D L. 1979. A Discrete Numerical Model for Granular Assemblies. Geotechnique , 29:47–65. doi:10.1680/geot.1979.29.1.47 |
| Cundall P A. 1982. Adaptive Density-Scaling for Time-Explicit Calculations//Proceedings of the 4th Inter-national Conference on Numerical Methods in Geomechanics (Edmonton, 1982) , 1, Rotterdam: Balkema:23-26. |
| Deen N G W, Dijkhuizen W, Bokkers G A. 2004. Validation of the granular temperature prediction of the kinetic theory of granular flow by particle image velocimetry and discrete particle model//The Third International Symposium on Two-Phase flow Modelling and experimentation: 22-23. |
| Di Renzo A, Di Maio F P. 2004. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes. Chemical Engineering Science , 59:525–541. doi:10.1016/j.ces.2003.09.037 |
| Domenico S N. 1977. Elastic properties of unconsolidated porous sand reservoirs. Geophysics , 42:1339–1368. doi:10.1190/1.1440797 |
| Dong Y, Zhang G, Sun Q, Zhao X, Niu X. 2015. Analysis of low-frequency vibrational modes and particle rearrangements in marginally jammed amorphous solid under quasi-static shear. Chinese Physical Letters , 32:126201. doi:10.1088/0256-307X/32/12/126201 |
| Dzubiella J, Hoffmann G P, Löwen H. 2002. Lane formation in colloidal mixtures driven by an external field. Physical Review E , 65:021402. |
| Gan J, Zhou Z, Yu A B. 2016. Particle scale study of heat transfer in packed and fluidized beds of ellipsoidal particles. Chemical Engineering Science , 144:201–215. doi:10.1016/j.ces.2016.01.041 |
| Goodrich C P, Liu A J, Nagel S R. 2014. Solids between the mechanical extremes of order and disorder. Nature Physics , 10:578–581. doi:10.1038/nphys3006 |
| Helbing D. 2001. Traffic and related self-driven many-particle systems. Review of Modern Physics , 73:1067–1141. doi:10.1103/RevModPhys.73.1067 |
| Hoomans B P B. 1999. Granular dynamics of gas-solid two-phase flows. [PhD Thesis]. Nertherland:University of Twente. |
| Hussainova I, Kubarsepp J, Shcheglov I. 1999. Investigation of impact of solid particles against hardmetal and cement targets. Tribology International , 32:337–344. doi:10.1016/S0301-679X(99)00073-0 |
| Iwashita K, Oda M. 1998. Rolling resistance at contacts in simulation of shear band development by DEM. Journal of Engineering Mechanics , 124:285–292. doi:10.1061/(ASCE)0733-9399(1998)124:3(285) |
| Iwashita K, Oda M. 2000. Micro-deformation mechanism of shear banding process based on modified distinct element method. Powder Technology , 109:192–205. doi:10.1016/S0032-5910(99)00236-3 |
| Jaeger H M, Nagel S R, Behringer R P. 1996. Granular solids, liquids, and gases. Review of Modern Physics , 68:1259–1273. doi:10.1103/RevModPhys.68.1259 |
| Jefferson G, Haritos G K, McMeeking R M. 2002. The elastic response of a cohesive aggregate-a discrete element model with coupled particle interaction. Journal of the Mechanics and Physics of Solids , 50:2539–2575. doi:10.1016/S0022-5096(02)00051-0 |
| Jiang M J, Yu H S, Harris D. 2005. A novel discrete model for granular material incorporating rolling resistance. Computers and Geotechnics , 32:340–357. doi:10.1016/j.compgeo.2005.05.001 |
| Jia X, Caroli C, Velicky B. 1999. Ultrasound propagation in externally stressed granular media. Physical Review Letters , 82:1863–1866. doi:10.1103/PhysRevLett.82.1863 |
| Jiang Y M, Liu M. 2009. Granular solid hydrodynamics. Granular Matter , 11:139–156. doi:10.1007/s10035-009-0137-3 |
| Johnson K L.1985. Contact Mechanics. England: Cambridge University Press . |
| Khidas Y, Jia X. 2010. Anisotropic nonlinear elasticity in a spherical-bead pack: Influence of the fabric anisotropy. Physical Review E , 81:021303. doi:10.1103/PhysRevE.81.021303 |
| Kolymbas D, Wu W. 2000. Introduction to hypoplasticity. Modern Approaches to Plasticity , 1:213–223. |
| Kuhn M R. 1999. Structured deformation in granular materials. Mechanics of Materials , 31:407–429. doi:10.1016/S0167-6636(99)00010-1 |
| Kruggel-Emden H, Smisek E, Rickelt S, Wirtz S, Scherer V. 2007. Review and extension of normal force models for the discrete element method. Powder Technology , 171:157–173. doi:10.1016/j.powtec.2006.10.004 |
| Kruggel-Emden H, Wirtz S, Scherer V. 2008. A study on tangential force laws applicable to the discrete element method (DEM) for materials with viscoelastic or plastic behavior. Chemical Engineering Science , 63:1523–1541. doi:10.1016/j.ces.2007.11.025 |
| Kuwabara G, Kono K. 1987. Restitution coefficient in collision between two spheres. Japanese Journal of Applied Physics , 26:1230–1233. doi:10.1143/JJAP.26.1230 |
| Langston, P A, Tuzun U, Heyes D M. 1994. Continuous potential discrete particle simulations of stress and velocity-fields in hoppers|transition from fluid to granular flow. Chemical Engineering Science , 49:1259–1275. doi:10.1016/0009-2509(94)85095-X |
| Lee J, Herrmann H J. 1993. Angle of repose and angle of marginal stability|molecular-dynamics of granular particles. Journal of Physics A , 26:373–383. doi:10.1088/0305-4470/26/2/021 |
| Lherminier S, Planet R, Simon G, Vanel L, Ramos O. 2014. Revealing the structure of a granular medium through ballistic sound propagation. Physical Review Letters , 113:098001. doi:10.1103/PhysRevLett.113.098001 |
| Li R, Yang H, Zheng G, Zhang B F, Fei M L, Sun Q C. 2016. Double speckle-visibility spectroscopy for the dynamics of a passive layer in a rotating drum. Powder Technology , 295:167–174. doi:10.1016/j.powtec.2016.03.031 |
| Li X, Li X S. 2009. Micro-macro quantification of the internal structure of granular materials. Journal of Engineering Mechanics , 135:641–656. doi:10.1061/(ASCE)0733-9399(2009)135:7(641) |
| Liu A J, Nagel S R. 1998. Jamming is not just cool any more. Nature , 396:21–22. doi:10.1038/23819 |
| Liu A J, Nagel S R. 2010. The jamming transition and the marginally jammed solid. Annual Review of Condensed Matter Physics , 1:347–369. doi:10.1146/annurev-conmatphys-070909-104045 |
| Liu J, Sun Q, Jin F. 2009. Visualization of force networks in 2D dense granular materials. Frontiers of Architecture and Civil Engineering in China , 4:109–115. |
| Liu X, Ge W, Li J. 2008. Non-equilibrium phase transitions in suspensions of oppositely driven inertial particles. Powder Technology , 184:224–231. doi:10.1016/j.powtec.2007.11.045 |
| Liu X, Martin C L, Delette G, Bouvard D. 2010. Elasticity and strength of partially sintered ceramics. Journal of the Mechanics and Physics of Solids , 58:829–842. doi:10.1016/j.jmps.2010.04.007 |
| Liu X, Martin C L, Bouvard D, Di Iorio S, Laurencin J, Delette G. 2011a. Strength of highly porous ceramic electraodes. Journal of the American Ceramic Society , 94:3500–3508. doi:10.1111/j.1551-2916.2011.04669.x |
| Liu X, Martin C L, Delette G, Laurencin J, Bouvard D, Delahaye T. 2011b. Microstructure of porous composite electrodes generated by the discrete element method. Journal of Power Sources , 196:2046–2054. doi:10.1016/j.jpowsour.2010.09.033 |
| Liu X, Papon A, Muhlhaus H. 2012. A numerical study of structural evolution in shear band. Philosophical Magazine , 92:3501–3519. doi:10.1080/14786435.2012.715249 |
| Lu G, Third J R, Muller C. R. 2015. Discrete element models for non-spherical particle systems: From theoretical developments to applications. Chemical Engineering Science , 127:425–465. doi:10.1016/j.ces.2014.11.050 |
| Luding S, Clément E, Blumen A, Rajchenbach J, Duran J. 1994. Anomalous energy dissipation in molecular dynamics simulations of grains: The "detachment" effect. Physical Review E , 50:4113–4122. doi:10.1103/PhysRevE.50.4113 |
| Luding S. 1998. Collisions and contacts between two particles//Herrmann H J, Hovi J P, Luding S. (Eds.). Physics of Dry Granular Media, Kluwer Academic Publs, Dordrecht.: 285-304. |
| Luding S. 2008a. Introduction to discrete element methods. European Journal of Environmental and Civil Engineering , 12:785–826. doi:10.1080/19648189.2008.9693050 |
| Luding S. 2008b. Cohesive, frictional powders: contact models for tension. Granular matter , 10:235–246. doi:10.1007/s10035-008-0099-x |
| Makse H A, Gland N, Johnson D L, Schwartz L. 2004. Granular packings: nonlinear elasticity, sound propagation, and collective relaxation dynamics. Physical Review E , 70:061302. doi:10.1103/PhysRevE.70.061302 |
| Manning M, Liu A. 2011. Vibrational modes identify soft spots in a sheared disordered packing. Physical Review Letters , 107:108302. doi:10.1103/PhysRevLett.107.108302 |
| Martin C L. 2004. Elasticity, fracture and yielding of cold compacted metal powders. Journal of the Mechanics and Physics of solids , 52:1691–1717. doi:10.1016/j.jmps.2004.03.004 |
| Maw N, Barber J R, Fawcett J N. 1976. The oblique impact of elastic spheres. Wear , 38:101–114. doi:10.1016/0043-1648(76)90201-5 |
| McNamara S, Garcia-Rojo R, Herrmann H J. 2008. Microscopic origin of granular ratcheting. Physical Review E , 77:031304. doi:10.1103/PhysRevE.77.031304 |
| Merkel A, Tournat V, Gusev V. 2014. Directional asymmetry of the nonlinear wave phenomena in a three-dimensional granular phononic crystal under gravity. Physical Review E , 90:023206. doi:10.1103/PhysRevE.90.023206 |
| Meyer H, Schulmann N, Zabel J E, Wittmer J P. 2011. The structure factor of dense two-dimensional polymer solutions. Computer Physics Communications , 182:1949–1953. doi:10.1016/j.cpc.2010.12.003 |
| Mindlin R D, Deresiewicz H. 1953. Elastic spheres in contact under varying oblique forces. Journal of Applied Mechanics , 20:327–344. |
| Mishra B K. 2003a. A review of computer simulation of bumbling mills by the discrete element method:part I|contact mechanics. International Journal of Mineral Processing , 71:73–93. doi:10.1016/S0301-7516(03)00032-2 |
| Mishra B K. 2003b. A review of computer simulation of bumbling mills by the discrete element method:part II|practical applications. International Journal of Mineral Processing , 71:95–112. doi:10.1016/S0301-7516(03)00031-0 |
| Misra A, Cheung J. 1999. Particle motion and energy distribution in tumbling ball mills. Powder Technology , 105:222–227. doi:10.1016/S0032-5910(99)00141-2 |
| O'Hern C S, Silbert L E, Liu A J, Nagel S R. 2003. Jamming at zero temperature and zero applied stress:the epitome of disorder. Physical Review E , 68:011306. doi:10.1103/PhysRevE.68.011306 |
| Peng Y, Wang Z, Alsayed A M, Yodh A G, Han Y. 2010. Melting of colloidal crystal films. Physical Review Letters , 104:205703. doi:10.1103/PhysRevLett.104.205703 |
| Renaud G, Calle S. 2010. Defontaine M 2010. Dynamic acoustoelastic testing of weakly pre-loaded unconsolidated water-saturated glass beads. The Journal of the Acoustical Society of America , 128:3344–3354. doi:10.1121/1.3502461 |
| Rognon P G, Roux J N, Naaim M, Chevoir F. 2008. Dense flows of cohesive granular materials. Journal of Fluid Mechanics , 596:21–47. |
| Sadd M H, Tai Q M. 1993. A contact law effects on wave-propagation in particulate materials using distinct element modeling. International journal of Non-Linear Mechanics , 28:251–265. doi:10.1016/0020-7462(93)90061-O |
| Schafer J, Dippel S, Wolf D E. 1996. Force schemes in simulations of granular materials. Journal de Physique , 16:5–20. |
| Schoenholz S S, Liu A J, Riggleman R A, Rottler J. 2014. Understanding plastic deformation in thermal glasses from single-soft-spot dynamics. Physical Review X , 4:031014. |
| Silbert L E, Silbert M. 2009. Long-wavelength structural anomalies in jammed systems. Physical Review E , 80:041304. doi:10.1103/PhysRevE.80.041304 |
| Song C, Wang P, Makse H A. 2008. A phase diagram for jammed matter. Nature , 453:629–632. doi:10.1038/nature06981 |
| Sun Q, Jin F, Wang G, Song S, Zhang G. 2015. On granular elasticity. Scientific Reports , 5:9652. doi:10.1038/srep09652 |
| Tanguy A, Mantisi B, Tsamados M. 2010. Vibrational modes as a predictor for plasticity in a model glass. European Physical Letters , 90:16004. doi:10.1209/0295-5075/90/16004 |
| Tejada I G, Jimenez R. 2014. Impact of the timestep in some molecular dynamics simulations on compression of granular systems. European Physical Journal E , 37:15. doi:10.1140/epje/i2014-14015-4 |
| Thornton C, Yin K K. 1991. Impact of elastic spheres with and without adhesion. Powder Technology , 65:153–166. doi:10.1016/0032-5910(91)80178-L |
| Thornton C. 1997. Coefficient of restitution for collinear collisions of elastic perfectly plastic spheres. Journal of Applied Mechanics , 64:383–386. doi:10.1115/1.2787319 |
| Thornton C, Cummins S J, Cleary P W. 2011. An investigation of the comparative behaviour of alternative contact force models during elastic collisions. Powder Technology , 210:189–197. doi:10.1016/j.powtec.2011.01.013 |
| Thornton C, Cummins S J, Cleary P W. 2013. An investigation of the comparative behaviour of alternative contact force models during inelastic collisions. Powder technology , 233:30–46. doi:10.1016/j.powtec.2012.08.012 |
| Tighe B. 2011. Relaxations and rheology near Jamming. Physical Review Letters , 107:158303. doi:10.1103/PhysRevLett.107.158303 |
| Tordesillas A, Muthuswamy M, Walsh S D. 2008. Mesoscale measures of nona±ne deformation in dense granular assemblies. Journal of Engineering Mechanics , 134:1095–1113. doi:10.1061/(ASCE)0733-9399(2008)134:12(1095) |
| Tsuji Y, Tanaka T, Ishida T. 1992. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technology , 71:239–250. doi:10.1016/0032-5910(92)88030-L |
| Utter B, Behringer R P. 2004. Self-diffusion in dense granular shear flows. Physical Review E , 69:031308. doi:10.1103/PhysRevE.69.031308 |
| Vitelli V. 2010. Attenuation of shear sound waves in jammed solids. Soft Matter , 6:3007–3012. doi:10.1039/c000834f |
| Vu-Quoc L, Zhang X, Lesburg L. 2001. Normal and tangential force-displacement relations for frictional elasto-plastic contact of spheres. International Journal of Solids and Structures , 38:6455–6489. doi:10.1016/S0020-7683(01)00065-8 |
| Vu-Quoc L, Zhang X. 1999. An accurate and e±cient tangential force-displacement model for elastic fric-tional contact in particle-flow simulations. Mechanics of Materials , 31:235–269. doi:10.1016/S0167-6636(98)00064-7 |
| Wang X, Zheng W, Wang L, Xu N. 2015. Disordered solids without well-defined transverse phonons: the nature of hard-sphere glasses. Physical Review Letters , 114:035502. doi:10.1103/PhysRevLett.114.035502 |
| Wang P J, Li Y D, Xia J H, Liu C S. 2008a. Characterization of reflection intermittency in a composite granular chain. Physical Review E , 77:060301. |
| Wang P, Song C, Briscoe C. 2008b. Particle dynamics and effective temperature of jammed granular matter in a slowly sheared three-dimensional Couette cell. Physical Review E , 77:061309. doi:10.1103/PhysRevE.77.061309 |
| Wang P J, Xia J H, Li Y D, Liu C S. 2007. Crossover in the power-law behavior of confined energy in a composite granular chain. Physical Review E , 76:041305. doi:10.1103/PhysRevE.76.041305 |
| Walton O R, Braun R L. 1986. Viscosity, granular temperature and stress calculations for shearing assemblies of inelastic, frictional disks. Journal of Rheology , 30:949–980. doi:10.1122/1.549893 |
| Walton O R. 1993. Numerical simulation of cline chute flows of monodisperse, inelastic, frictional spheres. Mechanics of Materials , 16:239–247. doi:10.1016/0167-6636(93)90048-V |
| Wang Y, Steffen A, Latham S, Mora P. 2006. Implementation of particle-scale rotation in the 3-D lattice model. Pure and Applied Geophysics , 163:1769–1785. doi:10.1007/s00024-006-0096-0 |
| Warr S, Hansen J P. 1996. Relaxation of local density fluctuations in a fluidized granular medium. Euro-physics Letters , 36:589–594. doi:10.1209/epl/i1996-00273-1 |
| Wen P P, Zheng N, Li L S, Li H, Sun G, Shi Q F. 2012. Polymerlike statistical characterization of two-dimensional granular chains. Physical Review E , 85:031301. doi:10.1103/PhysRevE.85.031301 |
| Wildman R D, Huntley J M, Parker D J. 2001. Granular temperature profiles in three-dimensional vibroflu-idized granular beds. Physical Review E , 63:061311. doi:10.1103/PhysRevE.63.061311 |
| Xia C, Li J, Cao Y, Kou B, Xiao X, Fezzaa K, Xiao T, Wang Y. 2015. The structural origin of the hard-sphere glass transition in granular packing. Nature Communications , 6:1–9. |
| Xu B H, Yu A B. 1997. Numerical simulation of the gas-solid flow in a fluidised bed by combining discrete particle method with computational fluid dynamics. Chemical Engineering Science , 52:2785–2809. doi:10.1016/S0009-2509(97)00081-X |
| Xu N. 2011. Mechanical, vibrational, and dynamical properties of amorphous systems near jamming. Fron-tiers of Physics , 6:109–123. doi:10.1007/s11467-010-0102-y |
| Xu N, Vitelli V, Liu A J, Nagel S R. 2010. Anharmonic and quasi-localized vibrations in jammed solids-modes for mechanical failure. Europhysics Letters , 90:56001. doi:10.1209/0295-5075/90/56001 |
| Xu N, Ching E S C. 2010. Effects of particle-size ratio on jamming of binary mixtures at zero temperature. Soft Matter , 6:2944–2948. doi:10.1039/b926696h |
| Yang H, Li R, Kong P, Sun QC, Biggs MJ, Zivkovic V. 2015. Avalanche dynamics of granular materials under the slumping regime in a rotating drum as revealed by speckle visibility spectroscopy. Physical Review E , 91:042206. doi:10.1103/PhysRevE.91.042206 |
| Ye M. 2005. Multi-level modeling of dense gas-solid two-phase flows. [PhD Thesis]. Nertherland: University of Twente. |
| Zaccone A, Terentjev E M. 2014. Short-range correlations control the G/K and Poisson ratios of amorphous solids and metallic glasses. Journal of Applied Physics , 115:033510. doi:10.1063/1.4862403 |
| Zhang Q, Li Y, Hou M, Jiang Y, Liu M. 2012. Elastic waves in the presence of a granular shear band formed by direct shear. Physical Review E , 85:031306. doi:10.1103/PhysRevE.85.031306 |
| Zhang Z X, Xu N, Chen D T N, Yunker P, Alsayed A M, Aptowicz K B, Habdas P, Liu A J, Nagel S R, Yodh A G. 2009. Thermal vestige of the zero-temperature jamming transition. Nature , 459:230–233. doi:10.1038/nature07998 |
| Zheng H P. 2014. Properties of surface waves in granular media under gravity. Chinese Physics B , 23:054503. doi:10.1088/1674-1056/23/5/054503 |
| Zhou Y H. 2013. Modeling of softsphere normal collisions with characteristic of coe±cient of restitution dependent on impact velocity. Theoretical and Applied Mechanics Letters , 3:021003. doi:10.1063/2.1302103 |
| Zhou Y C, Wright B D, Yang R Y, Xu B H, Yu A B. 1999. Rolling friction in the dynamic simulation of sandpile formation. Physica A , 269:536–553. doi:10.1016/S0378-4371(99)00183-1 |
| Zhu H P, Yu A B. 2003. The effects of wall and rolling resistance on the couple stress of granular materials in vertical flow. Physica A , 325:347–360. doi:10.1016/S0378-4371(03)00143-2 |
| Zhu H P, Zhou Z Y, Yang R Y, Yu A B. 2007. Discrete particle simulation of particulate systems: theoretical developments. Chemical Engineering Science , 62:3378–3396. doi:10.1016/j.ces.2006.12.089 |
| Zhu H P, Zhou Z Y, Yang R Y, Yu A B. 2008. Discrete particle simulation of particulate systems: A review of major applications and findings. Chemical Engineering Science , 63:5728–5770. doi:10.1016/j.ces.2008.08.006 |
