三方程线性弹性-阻尼DEM模型及碰撞参数确定

　　1　引言

　　自然界中存在着众多颗粒状物质,在外力与内部应力的共同作用下会产生类似于流体的运动状态。虽然颗粒运动与流体运动有些类似,但由于颗粒尺度远大于分子尺度,造成颗粒流动中会出现很强的不连续性,因此很难用传统的连续介质理论进行研究[1]。随着计算机技术的发展,近年来计算颗粒力学逐渐发展为研究颗粒体系的重要研究手段,其中软球模型DEM(又称离散单元法,Discrete el-ement method)因在处理高浓颗粒流中具有较高的效率而被广泛地应用到装填[2]、堆积[3]、混合[4]、粉碎[5]、仓储[6]、气固两相流[7]及流态化[8]等领域。自从Cundall和Strack提出DEM方法以来[9],很多研究者又发展了众多改进模型,但在这些模型中,大多数还很不完善,方法的准确性不高,比如在绝大多数现存的模型中,滚动摩擦力矩的影响没有被考虑在内,这与物理事实不符,因此无法准确地描述包含滚动的颗粒体系,比如在堆积过程的模拟中,不考虑滚动摩擦就无法得到稳定的沙堆。直到最近,滚动摩擦才被引入到DEM模型中,其重要性也得到充分的证实[10-12]。

　　根据弹性-阻尼的性质,DEM模型可分为线性模型和非线性模型,其中线性模型是最简单也是应用最广泛的模型。本文首先建立了完整的三方程线性弹性-阻尼模型,并将其应用到颗粒材料的三维模拟中,通过对几个典型颗粒系统的模拟来验证模型及算法的准确性。

　　2　离散单元模型

　　2.1　线性弹性-阻尼模型

　　颗粒在运动过程中主要受两种力作用,即自身重力和颗粒间的法向及切向碰撞接触力。根据牛顿第二定律,每个颗粒的平动运动方程为

　　在上述力和力矩的作用下,颗粒发生移动和滚动。式中m和I分别表示颗粒的质量和转动惯量,ni表示与颗粒i接触的颗粒总数,g为重力加速度。

　　当两个球形颗粒发生碰撞时,首先在接触点处发生弹性变形,颗粒在前进方向受到阻力,阻力的大小与法向变形位移、颗粒硬度成正比,达到最大位移变形时,颗粒停止运动,沿原来运动的方向反弹,碰撞后颗粒的动能会产生一定的损失,损失的大小与颗粒的弹性阻尼系数及颗粒间的相对速度有关。当两个颗粒发生偏心碰撞时,相撞点处的接触力可分解为法向分力和切向分力,其中切向力的极值受到颗粒表面摩擦系数与法向力乘积的限制,当切向分力大于该极值时,两颗粒在接触表面将发生滑动。同时,颗粒间还存在着滚动摩擦力矩,该力矩与滑动摩擦力类似,受到滚动摩擦系数的限制。上述模型考虑了法向接触力、切向接触力以及滚动摩擦力矩,每种作用力和力矩都可以简化为一个弹簧、一个阻尼以及一个滑动器,可用三组弹性-阻尼方程式来描述,称之为三方程线性弹性-阻尼模型,其中各种力和力矩的数学描述见表1。