颗粒-颗粒接触力的热力学模型∗

蒋亦民 刘佑

1)(中南大学物理与电子学院,长沙 410083)

2)(蒂宾根大学理论物理研究所,德国,蒂宾根 72076)

1 引 言

颗粒物质(诸如沙堆等)与传统晶体材料(例如半导体硅)的重要不同是多出一个介观(或细观)颗粒尺度.这使得从电子-原子核之间的二体库仑静电力出发,计算其物理力学性质的所谓“从头算(ab initio)”思想,虽然原则概念上是存在的,但实际操作上没有可能性.也就是说对硅晶体可以从头计算,但对颗粒物质最小尺度的分析计算只能从细观颗粒出发.另外为了计算数目庞大的颗粒系统,还不能像赫兹那样将颗粒处理成可变形的弹性固体[1],而是必须简化为只有六个运动学(kinematic)变量的、带转动的“质点粒子”.如果要求颗粒i的六个运动学量,即质心位矢和角位矢{ri(t),φi(t)}遵守牛顿力和力矩方程的话,则需要建立{ri,φi}与其他颗粒j之间的作用力fij关系模型来闭合牛顿方程组,然后对其数值求解,可得到该质点系统细观结构和力结构的时间演化(我们将省略与本文内容无关的边界条件).如果进一步选择合适的物理体积元做这些细观信息的统计平均(诸如粗粒化平均),还可得各种宏观平均量及其涨落的时空变化以及各种关联性质.这套称作软球离散元(DEM)模拟方法的一个优势是可以描述以长程关联为特征的弹性力链的呈现与消亡现象,又称固液转变或jamming转变.也就是说软球DEM适用于颗粒物质的所有类固液气态及其转变[2].

软球DEM面临的颗粒-颗粒接触作用力建模问题,在晶体物性的“从头算”领域里是不存在的.“从头算”认为库仑静电力能够解决所有凝聚态物理问题[3,4],但对颗粒物理而言接触力建模才是它的一个基础课题.由于这是个带耗散的不可逆非保守力,原则上需要在热力学概念基础上建模,具体地讲就是除了接触力外,还应该给出机械能功率与热功率的比例.或者给出计算弹性势能的表达式,因为接触时的总能量变化功率可以将作用力乘以运动学速度得到:fij(vi−vj)(这里vi≡dri/dt),动能的计算公式又是熟知的,故一旦明确了接触过程弹性势能的计算方式,问题也就彻底解决了.注意著名的赫兹法向接触力fn∼|ri−rj|3/2[1],是忽略耗散时的结果.有耗散时,特别是有塑性耗散时(这在切向运动时必须考虑),接触力是不能简单地表达为颗粒几何运动学变量的代数函数型这样的捆绑关系的.这个复杂性来自接触力中有弹性贡献,而弹性是需要用额外的独立变量来表征的概念[3,4]:即所谓的平移对称破缺后呈现的Goldstone量(弹应变或弹簧长度变化).很多模型都注意到了有塑性时弹簧长度变化与几何运动学相对量{ri−rj,φi−φj}的关系需要用微积分来描述,特别是切向.例如Luding模型在处理切向力时,就引入了一个切向弹簧长度变量,它与运动学量的联系是利用增量关系模型和一个历史累积变量来定义[5].有些人不用弹簧长度,直接建立切向力与运动学量的增量关系模型(参见文献[2])(工程领域常常将某些微积分方程写成增量关系形式,这个习惯很可能源自Truesdell[6]).长期以来,接触力建模研究一直处于这样的唯象力学阶段,并未去澄清热力学量,即热功率或弹性势能的定义.值得提到的是,很多(例如Hertz-Mindlin)三维接触力模型的可逆部分采用了虎克力乘以重叠量平方根的Boussinesq模型[7,8].这里虽然不涉及耗散,由于Boussinesq力不是保守的,我们仍无法知道弹性势能.

由于这些问题,虽然有少数DEM工作计算了弹性势能与动能比值[9,10],一个严格热力学的接触力模型仍值得建立.考虑到目前一个严格热力学的颗粒物质连续力学理论——颗粒固体流体动力学(GSH)已经存在[11−14],接触力的热力学模型完全可以用类似的方案得到(第二节).这时GSH中处理库仑屈服、塑性、滞迴等复杂力学现象的成功方法,可以直接继承到接触力上.例如前面提到的Boussinesq力没有势能的问题,完全可以连同库仑屈服一起,用一个类似GSH的势能模型来解决[15].由于各种能量耗散已经在细观接触的层面上明确定义,我们相信基于它的DEM计算能更好地描述长时间接触力链的热损耗和弹性势能及其演化(今后拟开展的工作).本文的第三节将针对二体力特有的碰撞恢复系数随碰撞初速度变化问题,做些简单的数值计算,以表明热力学方法能够计入这个实验现象.目前常用的Hertz-Mindlin和Luding模型的恢复系数都是材料常数,这个变化关系长期以来一直未能被描述(参见文献[2]).鉴于能量和耗散与接触力的关系比较复杂,我们将尽可能地给出详细解释和讨论,以方便读者了解热力学能量方法与以往的纯力学唯象方法的区别(第四节).

2 颗粒作用力模型

如图1,两个质量m1,m2,半径a1,a2的球形颗粒,实验室坐标下的质心位矢是r1,r2,相互作用力f.记相对位矢r=r1−r2,(法向)重叠量δ=d−r,其中r=|r|,d=a1+a2.两球位置的连线方向,即法向单位矢量n=r/r=(nx,ny).切向单位矢量t=(−ny,nx).记f的法向和切向投影为fn,t,有

为简单起见,本文假设图1的所有矢量都位于二维平面,并且忽略转动(即局限于考虑没有转动的二维运动,这些简化不影响将讨论的内容和概念).这时两球的运动将由牛顿方程

图1 两个运动颗粒的坐标和力示意图Fig.1.Coordinates and forces of two grains.

所支配.其中的v1,2=dr1,2/dt是速度,f1,2外力是两球受到的外力(考虑为已知,例如重力或边界墙力).另外记两球的相对速度为v=v1−v2,它的法向投影vn=vn,切向投影vt=vt.两球的切向相互位移是δt.用点表示时间导数,例如˙δ=dδ/dt.由(2)和(3)式可得相对速度的运动方程:

其中m=m1m2/(m1+m2)是折合质量.

显然一旦有了力fn,t=fn,t(r1,r2,v1,v2)的模型公式,上面的牛顿方程组就闭合了(外力是已知的).所谓的软球DEM模拟,就是在合适的边界条件下,对这样二体相互作用的大量颗粒系统进行数值求解.注意fn,t的建模不是一件容易的事情,因为除了保守的弹簧型作用力外,还需要考虑诸如库仑屈服打滑(slip)或塑性(plasticity)、黏力(dashpot)、滞迴(hysteresis)等一系列的不可逆耗散现象.我们面临的不是一个哈密顿保守力学问题,而是复杂的动-静摩擦力模型.另外根据赫兹的分析,法向力fn中的弹簧力部分肯定不是线性的.如引言里提到的,目前常用的DEM力模型的一个问题是能量的耗散情况不明确.例如含黏力Hertz-Mindlin模型:

其中的kn,t,ηn,t,µ是材料常数.文献中并没有明确给出相应的弹簧势能w弹的计算公式.这时如果期望分析计算力链的弹性势能就没有标准了.也许一个最好的做法是取

没有滑动时,它对δt的导数的确是ft公式(6)的弹性部分,但出现打滑就不好说了.另外它对δ的导数也不是fn公式(5)的弹性部分.总之对上述Hertz-Mindlin模型,我们觉得是不好定义弹簧势能.实际上Hertz-Mindlin模型可解释为用Boussinesq力描述法向和切向弹簧力,再加上黏力和库仑屈服.由于Boussinesq力不是保守力,因而无从知道相应的弹性势能.为改善以往颗粒-颗粒作用力模型的这个能量缺陷,下面将模仿GSH的二级不可逆耗散图像(见图2),提出一套基于能量守恒的建模方法.

一般而言如果要求弹性能w弹只是几何相对位置r(t)的代数函数,两球之间不会有滑动或塑性.也就是说含塑性变形的运动学(kinematic)变量不能直接用来描述弹性,需要额外引入表征接触处弹性变形或弹簧长度变化的位移矢量:

图2 颗粒-颗粒作用时的二级不可逆热力学示意图,图中的三个能量之和是守恒的Fig.2.Two-stage-irreversibility of grain-grain interaction.Sum of the three energies is conserved.

作为w弹的自变量(注意对本文考虑的二维运动,只讨论矢量的两个分量即可),即

弹 性 能w弹的 导 数πn=−∂w弹/∂un和πt=−∂w弹/∂ut将给出颗粒-颗粒相互作用力fn,t的弹性部分.可以不失一般地将它们写成:

这里的上标“热”和“松”,分别表示图2中左右两种耗散机制的力贡献.类似地可将弹性变形un,ut的运动方程写成

其中R松≥0是它的激发功率,I≥0是它向热能弛豫的衰减功率.松动能的激发功率是与πn,t的乘积,加上与vn,t的乘积:

相应地,有热功率

其中的I是来自缓冲区w松的贡献.最后,图2中机械能的定义是两个颗粒的动能与弹性能之和:

当没有外力时,能量守恒要求机械能w机械、松动能w松和热能之和保持常数,或它们的时间导数:

这可以直接计算证明.的确微分方程(16)有

代入(9)和(10)式,有

将(18)和(13),(15)式相加,即得(17)式.

与以往的模型对比,基于方程(7)—(17)理论构架下的接触相互作用除了给出弹性势能和热功率外,在处理弹性变形量与几何运动学量的关系上,以及机械能与热能的耗散关系上,引入了不那么僵硬的微分方程和能量缓冲区等措施.这将使得理论具有足够的灵活性来描述前面提到的各种复杂的接触力学现象(库仑屈服、打滑、塑性、黏力、滞迴等).

将取类似文献[15,16]的弹性势能模型

其中k,ξ,c是材料常数,c是描述湿黏力的常数[16].注意un和ut分别是法向和切向弹簧长度变化,量纲是米.弹性能量w弹量纲焦耳.弹簧系数k量纲为这里考虑了非线性的赫兹3/2幂率弹簧力.势能(19)给出的弹簧力是

这里要求法向弹簧的长度变化总是小于零:un＜0,但剪切弹簧的ut可正可负.后面将看到un＜0可以由热力学稳定来保障.

松动能的衰减可以写成弛豫时间模型的形式

(9),(10)和(11),(12)式中的耗散项,可以用标准的Onsager非平衡热力学来处理.仿照GSH的做法,有

是非对角(off-diagonal)迁移系数.如果把上面的矩阵方程写成非矩阵形式,有f松=ˆη松v,即类似地,有和和将它们代入前面的(7)—(17)式,得颗粒-颗粒相互作用力为:

弹簧变形量un,t的运动方程为:

松动能的运动方程为

另外热功率为

将(28)—(32)与牛顿方程(2),(3)联立求解,可计算两个颗粒运动轨迹.

以上接触理论的材料参数是:这里除弹性势能中的三个参数k,ξ,c必须是常数外,其他参数都是与耗散有关的迁移系数,并且可以是松动能w松、重叠量δ等变量的函数.考虑到接触力的复杂性,在这里保留了较多的迁移系数,以便理论有足够的灵活性来应对.

一个从GSH继承来的重要性质,是非线性弹性势能(19)有失稳现象,其对应的恰好是库仑屈服[15].势能(19)的力学稳定区域由不等式

给出,等号是稳定区域的边界方程.另外在稳定区域里总是有un＜0,故(19)—(21)不会出现虚数(因为实际运动总是发生在稳定区域内).利用(20),(21)式可将(35)式写成

失稳意味着一旦触及(35)或(36)式,接触力的弹性部分将快速松动.这个现象可以方便地通过让在屈服边界急剧增大来描述.

除失稳外,材料参数(34)式的进一步建模还需要顾及重叠量δ的影响,因为接触力、热功率、松动能等都是重叠时δ＞0才能被激发而出现的东西.一旦颗粒分离δ＜0,它们的激发源都会消失,同时也迅速衰减为零(假设没有湿黏力c=0).下面试用一组简单的材料参数值,以讨论碰撞恢复系数问题为例,对此做些具体说明.

3 碰撞恢复系数

作为一个简单例子,取下面的参数模型

和

其中Θ(x)是Heaviside阶梯函数,d=d1=d2和m=m1=m2分别是颗粒直径和质量(考虑两个相同颗粒).(37)式中假设了没有湿黏力,另外尽量略去了一些迁移系数.模型(38)意味着弹簧变形的激发源,即(30)和(31)式右边的第一项只在重叠时出现,分离后它们是零.黏滞系数(39)式也是这样.与GSH类似,(40)式假设切向弹簧变形的弛豫系数比例于(相当于GSH里的颗粒温度Tg).(41)式在满足稳定条件(35)时等于零,一旦不满足稳定条件,它将取适当的正值,起到尽快松弛切向弹簧变形、让系统迅速返回稳定区域的目的.注意上述模型的法向耗散与切向耗散有很大差异,后者比前者多考虑了(40)和(41)式两个系数.这显然有其合理性,因为切向耗散的确比法向的复杂,需要更多的迁移系数来描写.另外δ和un,t有长度量纲,为了让Heaviside函数的自变量是量纲1的,用颗粒粒径d对它们做了约化.另外发生颗粒重叠和屈服的数学描写总会涉及不连续性,需要用Heaviside函数来表达.这里是将它们反映在迁移系数(38),(39)和(41)式中,并且后者与弹性势能失稳判据(35)式有关.这个基于热力学概念的描述方法,与直接在力模型中唯象描写的传统做法是有差异的.

图3是上述参数值下数值计算的两个等质量和直径的、沿着x方向相对运动的颗粒做对心碰撞的情况.碰撞前后的速度分别是±v0和±vf(都沿着x方向),并且都是分离的,fn,t,un,t,w松都是零.颗粒重叠作用时,它们的激发、衰减、以及图2中三个能量之间的转换情况等,都能够明确计算(见图3).注意由于热功率R≥0,图3(c)的机械能(黑色曲线)总是随时间减小或不变.如果忽略碰撞过程细节,可以将这个接触作用模型看作以e=|vf/v0|为恢复系数的非弹性碰撞.

在分析长时间接触的力链演化问题时,除颗粒初始位置和速度外,还需要un,t,w松的初始值.通常可以取初始法向弹簧变形un等于初始重叠量δ.另外任何接触力模型(包括Hertz-Mindlin和Luding)都会出现Heaviside阶梯函数Θ(x). 数值分析时可用近似的解析函数来代替,例如取

本文的热力学接触力模型的恢复系数e=|vf/v0|不是材料常数,而是初速度v0和(斜碰撞(oblique impact)时的)碰撞长度的函数.对心碰撞时,本节考虑的模型的恢复系数随初速度增加而减小(图4i).这与实验定性符合[2],表明热力学建模方法有描述复杂观测细节的潜力.

图3 颗粒-颗粒对心碰撞过程的时间演化曲线 (a)位置坐标;(b)接触力(c)弹簧势能we,动能wk和机械能w;(d)松动能wrelaxation;是两个颗粒的初始能量(动能)Fig.3.Temporal evolution of normal collision of two grains:(a)Relative position;(b)contact force f=(c)potential of springs we,kinematic en-ergy wk,mechanic energy w;(d)relaxation energywrelaxation.The notation is initial(kine-matic)energy of two grains.

本节计算的目的仅限于用一个简单例子对热力学建模方法做些具体解释.实际颗粒的参数情况不一定是(37)—(41)式那样简单,参数值的标定也需要结合实验仔细处理.另外出于方便解释概念和方法,本文忽略了三维和颗粒转动运动φ(t).但在分析实际问题时,它们是应该计入的.

图4 对心碰撞恢复系数随初速度变化曲线Fig.4.Variation of normal restitution coefficient with collision velicity.

4 讨 论

为方便讨论本文热力学模型与其他接触力模型的概念差异,我们将弹性能方程(8)给出的力,和方程(28)—(31)写成下面的矢量形式:

1)弹性力π总是要求为保守力.也就是说有明确定义的弹性势能w弹(u).目前的DEM计算,虽然有些采用了保守的线性法向和线性切向弹簧,但有些是非保守的(例如Boussinesq力).也就是说对保守性的要求并不像热力学那样,是强制性的.非保守的弹簧模型会导致机械功率和热功率之间的分配比例不能清楚定义的问题.

2)弹簧应变u并不是纯几何的位矢r=∫vdt,而是由运动方程(43)描述的独立变量.这在目前的DEM计算中都注意到了(特别是切向分量),但u的运动方程与(43)式有很大不同.

3)接触力方程(44)中的黏力在当前DEM计算中大都有类似的对应(例如(5),(6)中的ηt,n项).但在对黏滞系数ˆη的处理方式上,没有图2所示的二级不可逆概念背景.

4)本文的接触力和几何运动学速度(f,v),与弹簧力和弹簧变形(π,u)之间的关系方程(43),(44)采用了典型的Onsager非平衡热力学形式.但这里除了熟知的机械能和热能外,还必须考虑一个中间缓冲能量w松去控制某些迁移过程和系数.这个称作“二级不可逆”的方案是颗粒物质热力学的基本特征.它不仅对宏观连续力学层面的GSH有效[11−14],细观层面的颗粒接触热力学也必须这样处理.它能把库仑屈服、弹性松弛、力与几何变形之间的滞迴和棘轮(ratcheting)等复杂力学现象,直接自然地与迁移过程联系起来,使得用热力学语言来理解和分析它们成为可能.从物理角度看这显然是合理的,因为颗粒物质的这些复杂力学现象都源自于耗散发热.

5)方程(43),(44)右边的第二项可分别理解为弹性弛豫和黏滞.右边的第一项涉及到以往未曾遇到的非对角迁移系数ˆG或ˆα=ˆ1−ˆG.如果忽略其他迁移系数,可将ˆα单独看作一个具有改变两端力和形变比值效果的“力学元件”.我们把它形象地称作“变速箱(gearbox)”,其能量和力学关系情况如图5所示.值得指出的是,如果变速比ˆG恒定(α=const.)并且u=(1−α)ε的话,可将图中的弹簧势能和力写成几何变形ε的函数:w弹=(1−α)2kε2/2和f=−(1−α)2kε,它们满足热力学关系f=−∂w弹/∂ε.这意味着如果不考虑α的动态效应,变速箱是不需要的,因为它的作用只是调整与其联接的弹簧的弹性系数,完全可以将其归入弹簧元件中.这也许是以往的各种力学元件都不曾提到变速箱的原因.但如果α是动态的,在考虑具有热力学背景的力学理论时,变速箱与弹簧就必须区别开来了.变速箱这个非对角迁移系数是我们近年来为理解颗粒物质宏观弹塑转变行为而提出的[17],它的统计物理和动理学意义还有待深入研究.唯象地看它的力学效果是:“在不增加热损失的前提下减弱弹性”．目前可以肯定的是,作为弹应变u的驱动源项的系数(见(43)式),它与弹性的出现与消亡密切相关,而且是动态的．本文为了简单,假设α随着由重叠量控制的接触弹性的出现与消亡,在0和1之间跳变,见模型(38)．当然实际情况可以复杂许多,例如还与松动能w松有关,在0和1之间变化等．

6)GSH热力学的一个重要特征是图2所示的二级不可逆耗散现象:除了左边的直接耗散为热的过程外,机械能还可以通过右边的缓冲能量w松进行耗散,同时迁移系数也可以受其控制(即ˆG,ˆλ,ˆη可以是w松的函数)．需要注意的是,库仑屈服现象虽然源自于弹性势能的静力失稳(见(35),(36)),我们还需要在某些迁移系数上采取相应的措施来保障任何动力学演化都不会进入该非稳区域,例如模型(41)．(因为任何观测到的实际系统都在热力学稳定区域内运动着,不会进入非稳区域．这与广为熟悉的范德瓦耳斯气体理论类似,它的热力学状态空间里也有一个实际气体不可能进入的非稳区域)．

图5 变速箱的能量与力学关系(注意α可以是动态的,但m,k假设为常数)Fig.5.Energetics and mechanics of gearbox.Notably α may temporally change,but m,k are assumed constant.

由于两个颗粒之间的很多相互作用行为是可以实验观测的,任何接触力模型原则上都应该用这些细观颗粒尺度的测量结果进行验证和标定．但以往的一些DEM工作发现,这个细观-宏观联系并不都是协调的,例如用线性弹簧计算的宏观行为,却比用更合理的赫兹弹簧要好(参见文献[12]的第六章)．也许这些矛盾和困难与所采用的接触力模型对热力学概念的描述不够全面有关．如果是这样的话本文的热力学建模方法将有助于问题的澄清．

目前颗粒物理研究方法可分为纯几何、DEM接触力、和GSH热力学三种．几何法致力于仅从颗粒的运动学变量(r,φ)出发,通过深入考察其几何时空结构与关联,以及几何构型熵(con fi guration entropy)或涨落耗散关系( fl uctuation-dissipation relation)等,来解释颗粒物质的各种物性(参见文献[18—21])．也就是说期望构建类似开普勒行星运动理论那样的颗粒物理．传统DEM方法的出发点是{r,φ,f},即比纯几何法增加了接触力f,并且期望在牛顿方程的基础上解释颗粒物质物性．由于接触力f不是可以用哈密顿描写的保守力,DEM一般都是直接数值求解大量颗粒的牛顿方程组．如果忽略f中的不可逆,仅保留保守的弹性力,也可以用传统的平衡态统计物理来分析这个哈密顿颗粒系统[22]．另外对f中的不可逆部分,还可以尝试用动理学[23]或mode-coupling[24]等非平衡统计技术来分析．但随着颗粒数密度增加,这些非平衡统计技术都面临困难,目前还无法满意地处理颗粒固体,特别是集体弹性行为的纳入和描述．因此数值求解牛顿方程仍是当前分析接触作用系统的最好办法(如引言里提到的,它能描述宏观弹性的呈现于消亡)．注意几何法和传统DEM法都可以不具体涉及热力学熵,可称作“非热(athermal)”类型理论．与之对比,热力学方法则强制性地要求无论是连续力学的宏观尺度,还是本文考虑的颗粒接触力细观尺度,热力学能量守恒和熵增加原理都应该得到明确的体现和遵守．当然三种方法之间有共同相似或互补的内容,但肯定还有不能互通的内容,因为牛顿力不可能取代热力学,也不可能被纯几何学取代．考虑到耗散发热是颗粒物质的基本性质,物理热力学应该是最为合理和可靠的方法．

值得提到的是,固体颗粒间的法向接触力一般没有塑性打滑等现象,模型(37)中取很可能是普遍合理的．但在Luding模型里法向滞迴是他的特色[5]．为方便今后与其对比,我们在Onsager矩阵(23)—(26)里保留了这几个迁移系数．迁移系数那些是无价值(irrelevant)可以去掉的,那些是有价值(relevant)需要保留的,及其相关的动力学意义等,目前还不能确定．这些问题需要结合DEM计算和实验情况来逐步澄清．模型(37)—(41)的参数数值仅作为参考,没有普适性．它们都是材料参数,应该依据实验来标定．另外将本文理论向三维推广时,只需要对文中的力、速度、位移等矢量补充一个切向分量即可．对球形颗粒,由于对称性两个切向的迁移系数是一样的,理论参数仍由(34)式列出．如果考虑(37)—(41)式的模型,这个三维理论给出的二维运动将退化为本文的公式．这时它给出的对心碰撞也就是本文的图3和图4．颗粒转动在一些接触力模型有所考虑,例如Luding模型[5]．但也有不少DEM工作省略了这个复杂性,例如文献[9]．本文出于简单和方便讨论热力学概念的原因也省略了颗粒转动．它的计入要复杂许多,建议参照Luding的做法．

常用接触力模型与本文建议的热力学方法有很多对应的和相似的地方,例如弹簧和黏滞阻尼．算例模型(37)—(41)式中的弹簧系数k(对应于刚度系数),也是常用模型中经常出现的常数．但热力学肯定还有自己的不同内容和改善的地方,例如将屈服和打滑放到了迁移系数中、给出了机械能和热功率的具体表达式、要求伴随任何迁移耗散过程的热功率都有“二次正定多项式”的形式(见(32),(33)式)等．注意“二次正定多项式”的要求来自热力学稳定性,它也会改善模型的数值稳定性．常用模型一般只有黏性的耗散热是平方正定的形式,其他复杂的不可逆现象往往没有这个性质,例如屈服．显然热力学内容的最大优势将体现在研究能量和耗散问题方面,这在最近已经开始受到重视(例如Yu小组的工作[25,26])．

本文建模方法与常用接触力模型的仔细对比,特别是具体改善和优越的地方,将在今后的定量对比工作中具体体现出来．这里可分为物理和工程两个方面来开展．物理方面,仅关心“二体力”这个基础层面即可,也就是说只是对比不同方法和模型的“两个颗粒动力学”情况即可,这样便于讨论和澄清与基本概念相关的优劣对比,同时避开大规模数值编程的麻烦．工程方面,主要考虑大量颗粒的DEM计算情况的对比,涉及大规模数值编程或商业软件．另外“二体力”和“DEM”两边都应该尽可能多地考虑与实验的对比．

在考虑计入颗粒转动因素时,首先是增加描述颗粒转动角的运动方程(即牛顿力矩方程),在机械能公式(16)中增加颗粒转动动能,以及可能的扭力势能等．然后是滚动力(rolling)和扭转力(torsion)的建模问题．由于热力学的影响主要体现在不可逆耗散部分,建模时涉及的纯几何运动学(kinematic)部分,以及滚动力和扭转力的可逆部分,基本可直接挪用Luding的相应公式[5]．但滚动力和扭转力的不可逆部分,需要用与本文类似的能量守恒和Onsager迁移系数等方法来处理．

5 结 论

虽然耗散发热是颗粒间接触运动时的重要现象,当前DEM采用的接触力模型仍属于唯象力学类型,对热力学要求的机械功率和热功率等的定义并不十分清楚．本文的主要结论是:基于热力学的接触力建模方法可以参照已有的GSH概念来实现．我们相信从唯象力学发展到热力学模型,不仅仅从物理概念上更加合理,对诸如恢复系数行为等的一些实验现象的描述,也会有所改善．当然在尝试将本文内容用于DEM数值计算之前,还需要先将其推广到三维和计入颗粒转动运动．这些将是今后需要陆续开展的工作．

感谢厚美瑛、孙其诚、程晓辉、刘晓星和张国华的讨论.

[1]Hertz H 1881J.Reine Angew.Math.92 156

[2]Thornton C 2015Granular Dynamics,Contact Mechanics and Particle System Simulations—A DEM Study(Particle Technology Series,Volume 24)(eBook DOI 10.1007/978-3-319-18711-2)(Switzerland:Springer International Publishing AG)

[3]Laughlin R B,Pines D 2000Proc.Natl.Acad.Sci.USA97 28

[4]Laughlin R B 2004A Different Universe(Changsha:Hunan Science and Technology Press)(in Chinese)[王文浩译2008不同的宇宙(长沙:湖南科学技术出版社)]

[5]Luding S 2008Granular Matter10 235

[6]Truesdell C 1972Rational Thermodynamics(Berlin:Springer-Verlag)

[7]Boussinesq J 1873C.R.Hebd.Seances Acad.Sci.77 1521

[8]Bonneau L,Andreotti B,Clément E 2007Phys.Rev.E75 016602

[9]Sun Q C,Jin F,Zhou G D 2013Granular Matter15 119

[10]Kumar N 2014Ph.D.Dissertation(Enschede:University of Twente.The Netherlands ISBN:978-90-365-3634-9,DOI:10.3990/1.9789036536349)

[11]Jiang Y M,Liu M 2009Granular Matter11 139

[12]Sun Q C,Hou M Y,Jin F 2011Physics and Mechanics of Granular Matter(Beijing:Science Press)(in Chinese)[孙其诚,厚美瑛,金峰 2011颗粒物质物理与力学(北京:科学出版社)]

[13]Jiang Y M,Liu M 2014Acta Mech.225 2363

[14]Jiang Y M,Liu M 2015Eur.Phys.J.E38 15

[15]Jiang Y M,Liu M 2003Phys.Rev.Lett.91 144301

[16]Jiang Y M,Liu M 2004Phys.Rev.Lett.93 148001

[17]Jiang Y M,Liu M 2007Phys.Rev.Lett.99 105501

[18]Torquato S,Stillinger F H 2010Rev.Mod.Phys.82 2633

[19]Edwards S F,Moun field C C 1996Physica A226 1

[20]Edwards S F,Moun field C C 1996Physica A226 12

[21]Edwards S F,Moun field C C 1996Physica A226 25

[22]Parisi G,Zamponi F 2010Rev.Mod.Phys.82 789

[23]Luding S 2009Nonlinearity22 R101

[24]Hayakawa H,Otsuki M 2008Prog.Theor.Phys.119 381

[25]Kuang S B,Zou R P,Pan R H,Yu A B 2012Ind.Eng.Chem.Res.51 14289

[26]Hou Q F,Kuang S B,Yu A B 2017Chem.Engineer.Sci.161 67