散体系统冲击破碎的动力学分析＊

姜世平，于海龙，芮筱亭，洪 俊，黎 超

（1.南京理工大学发射动力学研究所，江苏 南京210094；2.东南大学土木工程学院工程力学系，江苏 南京210096）

由大量没有机械原件联接的单个物体所组成的系统称为散体系统。散体系统广泛存在于自然界，如沙堆、矿石、岩土、碎煤、水泥和谷物等。散体系统是一种有别于固体和流体的物质形态，是由大量的固体单元组成的不连续系统。从微观角度来看，散体系统与固体一样，都是由固体颗粒组成，但是散体系统只能承受压力和一定的剪切力而一般不承受拉力，这些特性与流体相似。散体系统破碎动力学是研究散体系统运动、受力及破碎规律的力学学科，在地震、滑坡、泥石流、发射装药发射安全性等工程领域具有重要的理论与应用价值［1－2］。

散体系统破碎既是物体数量巨大的多体动力学问题，又是物体破碎的固体力学问题，因此对于散体系统破碎过程的数值模拟异常困难［2］。建立在传统连续介质力学基础上的有限元法、有限差分法等，适用于预测损伤和破坏的区域，却难直接计算材料及结构发生破坏的整个过程［3］。P.A.Cundall［4］提出了离散单元法，为散体系统破碎动力学提供了有力工具。离散单元法的基本思想是，将连续体分离为铰接而成的刚体系统，应用牛顿第二定律描述该多刚体系统的运动，求得连续体的整体运动形态。G.A.Addeta等［5］采用离散单元法，模拟了二维物体在静态和冲击载荷作用下的破碎过程；H.A.Carmona等［6］采用离散单元法，模拟了三维脆性球体单元在不同速度下的冲击破碎过程。

本文中，从离散元法的基本原理出发，将散体单元离散成弹簧－球单元系统，建立散体系统冲击破碎动力学模型，编制散体系统冲击破碎数值模拟程序，再现散体系统在冲击载荷下的破碎过程。

1 弹簧－球单元离散模型

离散单元法中的单元连接形式可分为接触型和连接型。接触型是散体系统特有的连接形式，由P.A.Cundall［4］提出，已发展成为计算边坡稳定性、散体系统运输等的重要工具。连接型模型中，单元间没有间隙且符合变形协调条件，用于处理连续介质力学问题。对于刚体运动、材料损伤和破坏等非线性力学问题，连接型模型具有明显的优越性。这是因为，离散单元的节点在单元的形心上，只需实行连接形式从连接型到接触型的转换，无需改变单元类型或重分网格，就可以实现连续体到非连续体的转变。

弹簧－球单元离散模型的基本思想是，将连续体离散成刚性球单元的集合，任意相邻的两个球体单元之间用一个弹簧组连接。如图1所示，弹簧组包括一个法向弹簧和两个切向弹簧，并用相应的准则判别弹簧的断裂，介质的损伤、破坏通过弹簧的变形、断裂体现，且小球单元是其发生破坏时的最小单位。

图1 连接球单元的弹簧组模型Fig.1The model of the springs linked with sphere elements

图2 球形颗粒的离散模型Fig.2The discrete model of spherical grains

如图2所示，将散体系统中的每个球体离散成大小相同的刚性球单元系统，任意相邻的球体单元之间由一个法向弹簧和两个切向弹簧连接。离散后的模型与原来的球形散体颗粒模型不可能完全一致，但随着离散单元半径的减小，这种缺陷对计算结果不会产生本质的影响。

2 离散单元的运动方程

采用离散单元法模拟散体系统冲击破碎过程。以牛顿第二定律为基础，计算由离散单元组成的系统在准静态和动态条件下的变形过程，尤其是大变形和大范围运动。对每个单元，计算单元的所有作用力，求出合力与合力矩，列出动力学方程，通过差分格式解出每个时步的速度和位移，并对时域积分，就可获得任意单元的速度和位置。

如图3所示，建立全局坐标系OXYZ，任取单元i，单元i可能同时和多个单元接触，任取接触单元j。取单元i的球心o作为局部坐标系的原点，以单元i的中心指向单元j的中心为x轴，过球心o取平行于OXY平面且垂直于x轴的直线为y轴，z轴由右手螺旋法则确定，由此建立局部坐标系oxyz。

图3 坐标系及单元接触情况Fig.3Coordinate system and contact status of elements

根据牛顿运动定律，单元i在坐标系OXYZ中的平动方程为：

由动量矩定理，单元i在坐标系oxyz中的转动方程为：

式中：mi为单元i的质量为单元i的质心线加速度，fij为单元i受到接触单元j的接触力，fi为单元i所受除接触力以外的外力，Ji为单元i的转动惯量，ωi为单元i的角速度，Μij为单元i与单元j接触所产生的接触力矩，Μi为单元i所受除接触力矩以外的外力矩。

3 弹性力计算模型

离散后的散体单元由多个等直径的小球组成，任意相邻的两个小球之间由一个弹簧组连接，当两个相邻的小球之间处于弹簧连接状态时，弹簧组变形产生的弹性力为：

式中：kn，ij、ks1，ij、ks2，ij分别为法向弹簧和两个切向弹簧的刚度系数，Δu1、Δu2、Δu3分别为法向弹簧和两个切向弹簧的变形量。对于弹簧组刚度系数的选取，文献［7 －8］中给予了详细的推导。

4 球体单元的接触模型

当单元i和j接触时，在法线和切线方向上，弹性效应等效于弹簧－阻尼器作用，如图4所示。法线、主切线及次切线方向的刚度系数分别为、和，阻尼系数分别为、和，具体计算公式见文献［9］。法向接触力可表示为［10］：

图4 接触模型Fig.4Contact model of spherical bodies

式中：Δδn＝ri＋rj－rij，为单元间法向相对变形量，ri、rj分别为单元i和j的半径，rij为球心之间的距离；指数α反映了接触的性质，对球体之间的接触，取α＝3/2；nij为法向单位矢量，方向由i指向j；为单元间接触点的相对速度矢量在法向的分量，可由下式求得：

式中：vi、vj分别为单元i和j的质心速度；ωi、ωj分别为单元i和j的角速度。

对于切向接触，由于在接触过程中有可能发生从静止到滑移或者由滑移到静止的过渡，所以接触模型应采用增量形式：

式中：μ为单元间滑动摩擦系数。

5 弹簧的断裂准则

采用Mohr －Coulomb破坏准则［8］，将材料的破坏形式简化为两种状态。在第一种状态下，单元间受压力、拉伸力及剪切力，当拉伸力达到拉伸破坏极限时被拉断，单元间不能承受拉伸力，转变为第二种状态，单元间只能承受压力和剪切力，即单元开始就处于分离状态或处在单元达到破坏极限后的状态。此状态即为单元间由连接关系转为接触关系。

弹簧的拉伸强度Fs、压缩强度Fc与黏着力F可以根据与应变率˙εc有关的破坏参数获得［11］：

式中：r为离散单元的半径，fs、fc、c分别为静载条件下散体单元的拉伸强度、压缩强度与黏着力强度。

弹簧断裂意味着产生裂纹，并且球体单元之间的关系由连接状态转化为接触状态，接触力可按式（4）、（6）进行计算，每个单元的运动由离散单元法描述。当连结某个离散球体单元的所有弹簧都断裂时，该球体从连续体中分离出来，散体单元发生了破碎。

图5 施加于活塞上的冲击载荷Fig.5The impact load acting on the piston

图6 散体系统的离散模型Fig.6The discrete model of granular system

6 散体系统冲击破碎数值模拟

应用上述方法，对圆筒内由脆性材料组成的散体系统在冲击载荷下的挤压破碎过程进行数值模拟，冲击载荷如图5所示。将圆筒内的球形散体单元离散成大量的小球，如图6所示。任意时刻每个单元的运动由牛顿－欧拉方程描述；对有弹簧连接的单元，弹性力计算采用式（3），当拉伸、压缩或剪切力达到弹簧最大承载极限时弹簧断裂。当球体单元之间接触时，法向接触力和切向接触力分别采用式（4）、（6）计算，离散单元与边界的接触力计算模型与离散单元之间的接触力计算模型一致；散体单元数量为424，散体单元半径为5.4mm，散体单元密度为3.6t/m3，离散单元数量为57×424，离散单元半径为1.08 mm，离散单元弹性模量为1.0GPa，泊松比为0.26，弹簧拉伸强度为1.2kN，弹簧压缩强度为3.6kN，黏着力为1.8kN，单元间摩擦系数为0.1，单元边界间摩擦系数为0，时间步长为10－7s，时间长度为3s。

图7为采用Open GL显示方法的计算结果。在1ms之前，由于活塞上的载荷小于散体系统的承载能力，散体系统处于稳定状态而没有发生破碎；在1～2ms之间，作用在活塞上的冲击载荷达到了最大值，2ms时散体系统发生了大量的破碎和变形；在2ms之后，冲击载荷迅速衰减并达到稳定状态，在外界载荷的继续作用下，散体系统中仍然有少量的散体单元发生破碎。

图7 散体系统的破碎过程Fig.7The fragmentation process of granular system

计算过程中，跟踪每一个散体单元在不同时刻的破碎情况。图8分别为圆筒底部、中部和上部的散体单元的破碎过程。

在散体系统冲击破碎的动力学数值模拟过程中，影响散体系统破碎过程的因素主要有：（1）离散单元半径的选取不能太大，因为离散单元是散体系统破碎的最小单位，若离散单元半径选取过大，则离散模型不足以完全刻画散体系统破碎过程中发生的微观破碎，导致冲击破碎过程中诸如粉末状的微小破碎颗粒没有被体现出来。一般来说，离散单元半径取得越小，计算结果越准确，但是半径选取过小，离散单元的数量变得非常庞大，计算时间也将呈指数形式增长，所以需综合考虑计算精度和计算时间之间的平衡。（2）散体系统不同的初始堆积构型对计算结果的影响非常明显。相同的计算条件下，不同的初始堆积构型导致散体系统破碎程度不同，这个现象经常出现。（3）计算参数能否准确地反映材料特性，对计算结果非常重要。

图8 散体单元的破碎过程Fig.8The fragmentation process of granular grain

7 结 论

针对球形粒子组成的散体系统，基于离散单元法，建立了散体系统冲击破碎的动力学数值模拟方法，有效地实现了对散体系统冲击破碎过程的计算。本方法可为其他散体系统接触、碰撞、破碎过程的分析提供参考。

［1］ 吴爱祥，孙业志，刘湘平.散体动力学理论及其应用［M］.北京：冶金工业出版社，2002.

［2］ 芮筱亭，贠来峰，王国平，等.弹药发射安全性导论［M］.北京：国防工业出版社，2009.

［3］ 江见鲸，陆新征，叶列平.混凝土结构有限元分析［M］.北京：清华大学出版社，2005：56 －64.

［4］ Cundall P A.A computer model for simulating progressive large scale movement in block rock system［C］∥Symposium ISRM.1971：129 －136.

［5］ Addeta G A，Kun F，Ramm E.On the application of a discrete mode to the fracture process of cohesive granular materials［J］.Granular Matter，2002，4（2）：77 －90.

［6］ Carmona H A，Wittel F K，Kun F，et al.Fragmentation process in impact of spheres［J］.Physical Review E，2008，77（5）：051302.

［7］ Shan Li，Cheng Ming，Liu Kai －xin，et al.New discrete element models for three －dimensional impact problems［J］.Chinese Physics Letters，2009，26（12）：120202.

［8］ 刘凯欣，高凌天.离散元法在求解三维冲击动力学问题中的应用［J］.固体力学学报，2004，25（2）：181 －185.Liu Kai －xin，Gao Ling －tian.The application of discrete element method in solving three －dimensional impact dynamics problems［J］.Acta Mechanica Solida Sinica，2004，25（2）：181 －185.

［9］ Mishra B K，Murty C V R.On the determination of contact parameters for realistic DEM simulation of ball mills［J］.Powder Technology，2001，115：290 －297.

［10］ Kruggel －Emden H，Wirtz S，Scherer V.Applicable contact model for the discrete element method：The single particle perspective［C］∥Pressure Vessels and Piping Conference（PVP2008）.Chicago，2008.

［11］ 目黑公郎，博野元彦.用粒状体模拟混凝土结构的破坏分析［J］.东京大学地质研究所汇报，1991，63（4）：409 －468.Mekuro K，Hakano T.Fracture analysis of concrete structures by granule simulation［J］.Research Report of Geology Institute of Tokyo University，1991，63（4）：409 －468.