球形颗粒随机排列过程的计算机模拟

D0I:10.13374/j.issn1001-053x.1999.0M.020 第21卷第4期 北京科技大学学报 Vol.21 No.4 1999年8月 Journal of University of Science and Technology Beijing Aug.1999 球形颗粒随机排列过程的计算机模拟 程远方果世驹 赖和怡 北京科技大学材料科学与工程学院，北京100083 摘要利用离散元方法模拟了均匀球形颗粒的随机排列过程.在模拟过程中考虑了重力、颗 粒间的接触力、摩擦力以及范德华力(VDWI):颗粒的运动包括平动和转动.研究表明，颗粒间 的作用力对球形颗粒随机排列体的结构影响很大.对于粒径为100μm球形颗粒，如果不考虑 颗粒之间的摩擦力和范德华力，排列体为随机密排和有序密排的混合体，其排列密度为0.8696： 考虑颗粒间的摩擦力和范德华力，排列密度降低为082！3.描述了排列体儿何结构随时间的演 化过程、特定颗粒在排列过程中的运动轨迹及颗粒配位数的分布规律.同时给出了不同条件下 二元系的排列结果。 关键词离散元：随机排列：排列密度：配位数 分类号TB44 由于球形颗粒随机排列的几何特征在许多 粒的随机排列过程并描述其几何特征， 技术领域都非常重要，特别是在粉末冶金及工 程陶瓷领域，压坯中颗粒的排列结构直接影响 1模型 着烧结后材料的性能.因此人们对此已进行了 离散元法是Cundall最先提出并用于分析 多年研究.60年代，通过镶嵌技术测量了随机 岩石力学问题的一种数值方法.Cundall和Hat 排列光滑钢球的位置坐标，研究了这种排列体 等人详细地介绍了此方法.这种方法允许离 的排列密度、配位数分布频率和径向分布函数 散体做有限的位移和转动，包括颗粒之间原接 等几何指标).公认有两种随机排列形式：随机 触的解离和新接触的形成.它是一种时间步长 密排，二维条件下排列密度为0.8225：随机松 增量算法，在每一步长内都重复应用牛顿第二 排，二维条件下排列密度为0.73.由于试验过程 定律和力一位移接触定律.牛顿第二定律确定颗 非常繁琐，并且测量精度和颗粒直径的微小变 粒在合力作用下的运动，而力-位移接触定律根 化都会使颗粒分布函数和排列体的几何形状表 据颗粒的位移确定颗粒接触点的接触力.本文 现出很大的不确定性.80年代，研究者都试图 假设颗粒为球形，用球心坐标和半径来表征.排 用计算机程序得到球形颗粒的随机排列体~， 列过程的物理机制通过受力分析来体现.分析 并提出了相应的（二维的或三维的）模拟均匀球 中考虑的力有：(1)重力：(2)颗粒之间及颗粒与 形颗粒随机排列的模型，如通过引进“中心场 盒壁之间的接触力；(3)颗粒之间及颗粒与盒壁 力”，不断填充“种晶”间隙可以形成随机密排结 之间的摩擦力：(4)邻近颗粒之间的分子间作用 构刃.通过简化沉降模型可以进行随机松排例. 力. 但是几乎所有的排列方法都忽略了颗粒排列的 本方法做如下假设： 物理机制，也不能给出颗粒体结构随时间的演 ●颗粒为刚性体： 化过程.实际上粉末颗粒表面是不光滑的，颗粒 ●接触面积与颗粒本身相比非常小，即为 之间存在着摩擦力和分子间作用力，这些力使 点接触： 得实际颗粒的随机排列结构与光滑小球的随机 ●接触力与颗粒重迭大小有关，由力-位移 排列结构存在明显的差异，本文试图考虑颗粒 接触定律确定， 排列过程的物理机制，用离散元法模拟实际颗 1.1颗粒上的作用力 1998-11-10收稿程运方男.34岁，副教授 对于任意颗粒a来说，作用其上的力有重 *国家“九五”攻关重点课题(N0.95-527-03-04-01) 力、接触力、摩擦力和范德华力.建立直角坐标

第 21 卷 第 4 期 1 9 9 9 年 8 月 北 京 科 技 大 学 学 报 J o u r n a l o f U n i v e r s i ty o f S e i e n e e a n d Te e h n o l o g y B e ij i n g V b l . Z I N 0 . 4 A u g . 1 9 9 9 球形颗粒随机排列过程 的计算机模拟 程远方 果世驹 赖和怡 北京科技大学材料科学 与工程学院 , 北京 10 0 0 8 3 摘 要 利用 离散 元方法 模拟 了均 匀球 形颗 粒 的随机 排列 过程 . 在模拟 过程 中考虑 了重力 、 颗 粒 间 的接 触 力 、 摩 擦力 以及 范德华 力 ( V D W D ; 颗粒 的运 动包 括平 动和 转动 . 研 究表 明 , 颗粒 间 的作用 力对 球形颗 粒随 机排列 体 的结构 影响 很大 . 对 于粒径 为 1 00 林m 球形 颗粒 , 如 果不考虑 颗粒之 间的摩擦 力和 范德华力 , 排 列体为 随机密排 和有序 密排的混合 体 , 其排 列密度为 .0 8 69 6 ; 考虑 颗粒 间的摩擦 力和范德 华力 : 排 列密度 降低 为 .0 82 1 3 . 描述 了排列 体儿 何结构 随时间 的演 化过 程 、 特定 颗粒在排 列过程 中的运 动轨迹 及颗 粒配位 数的 分布规 律 . 同 时给 出了不 同条件下 二元系 的排 列结果 . 关键词 离散 元 : 随机排 列 ; 排列 密度 ; 配位数 分类号 T B 4 4 由于 球形 颗粒随机排列的几何特征在许多 技术 领域都非常重要 , 特 别是在粉末冶 金 及 工 程 陶 瓷 领域 , 压坯 中颗粒 的排 列结构直接 影 响 着 烧结后 材料 的性 能 . 因 此人们对此 已进行 了 多年研 究 . 60 年代 , 通过镶 嵌技术测量 了 随机 排列 光滑 钢球 的位置 坐 标 , 研 究了这 种排列 体 的排 列密度 、 配位 数分布频率和 径 向分布 函数 等几何指标 ` 1,2] . 公 认 有两 种随机排列形 式 : 随机 密排 , 二 维条件下排列 密度为 .0 8 2 5 ; 随机松 排 , 二 维条件下 排列密度为 0 . 73 . 由于试验过程 非常繁琐 , 并且 测量 精度和 颗粒直 径的 微小变 化 都会使颗粒分布 函数和排列体的几 何形 状表 现 出 很大的不 确定性 . 80 年代 , 研究者 都试 图 用计算机程序得到球形颗 粒的随机排列体 { , 一 ” , 并提 出 了相应 的 ( 二 维的或三 维 的) 模拟均匀球 形颗粒 随机 排列 的模型 . 如通 过引进 “ 中心场 力 ” , 不断填充 “ 种 晶 ” 间隙可 以形 成随机密排结 构 ` 6,7] . 通过简化沉 降模型可 以进行 随机松 排 〔6] . 但是 几乎 所有 的排 列 方法都忽 略 了颗 粒 排列 的 物理机 制 , 也 不 能给 出颗粒体结构随时间的演 化过 程 . 实际 上粉末 颗粒表面 是不 光 滑的 , 颗粒 之 间 存在着摩 擦力和 分 子 间作用力 , 这 些力使 得实 际颗粒 的 随机排列 结构与光 滑小球 的随机 排 列结构存在 明显 的差 异 . 本文试 图考虑 颗粒 排列过程 的物理机 制 , 用 离散元 法 模拟 实际 颗 19 9 8 一 l 卜 1 0 收稿 程 运方 男 . 34 岁 , 副教授 * 国家 “ 九五 ” 攻 关重 点课 题 ( N o . 9 5 ~ 5 27 一 03 一 0 4 一 01 ) 粒 的 随机 排列 过程并描 述其 几 何特 征 . 1 模型 离散元法 是 C u n da l 最 先提出 并用于 分析 岩石力学问题 的一 种数值方法 LS] . c un da n 和 H a rt 等人 详细地 介绍 了此方法 〔9,1 01 . 这种 方法允许离 散体做有 限 的位移和 转动 , 包 括颗 粒之间原接 触 的解 离和 新接触 的形成 . 它 是一 种 时 间 步长 增量算法 , 在每一 步长 内都重 复应用 牛顿第二 定 律和 力一位移接触定 律 . 牛顿第二 定律确定颗 粒在合力作用 下 的运动 , 而 力 一 位移接触定 律根 据颗粒 的位移确 定 颗粒接触点 的接触力 . 本文 假设颗粒 为球形 , 用球心坐标和 半径来表 征 . 排 列过 程 的物理 机制 通过 受力分析来体现 . 分析 中考虑的力有 : ( l) 重 力 ; ( 2) 颗粒之 间及 颗粒与 盒壁之 间的接触力 ; (3 ) 颗 粒之 间及颗粒 与盒壁 之 间的摩擦力 ; ( 4) 邻近颗粒之 间的分子 间作用 力 . 本方法 做如下 假 设 : . 颗粒为刚性体 ; . 接触 面积 与颗粒本 身相 比非常小 , 即 为 点接 触 ; . 接触力 与颗粒重 迭大小有关 , 由力一 位移 接触 定律确定 . 1 . 1 颗粒上的作用 力 对 于任 意颗粒 口 来说 , 作用其 上 的力有重 力 、 接 触力 、 摩擦 力和 范 德华力 . 建立 直 角坐标 DOI: 10. 13374 /j . issn1001 -053x. 1999. 04. 020

-388· 北京科技大学学报 1999年第4期 系xoy,重力的方向为y轴的负向，则颗粒a的 「F,=m成 M,=I0 (9) 重力可表示为： G,=-mg (1) 式中I为颗粒a的转动惯量，运用中心差值法， 式中，m,为颗粒的质量，g为重力加速度，i为张 可得方程(9)的数值表达式.时间由t时刻增加 量表示法中的自由指标，此式i=2,若两颗粒近 到t+△t时刻，则t+△t时刻的速度为： 邻，则颗粒之间的范德华力可表示为： aa=-wm+EAr m R=器& (2) (10) 式中，A为Hamaker常数，对于钢来说等于 -rn+4y 21.1×10-J,R为颗粒半径，h为两颗粒表面之间 t+△：时刻的球心位置坐标为： 的距离，e,为两颗粒中心连线方向的单位矢量， x"a=xX0+*a2m+△t (11) 1.2接触定律 在模拟开始，颗粒之间没有接触而自由下 图1表示相互接触的两球a,b,其重叠量为 落，一定时间后，颗粒的速度达到最大值，颗粒 ,则法向接触力F由下式确定 之间开始出现接触.此后，在每一个时间步长开 F=KU (3) 始，根据颗粒和盒壁的位置确定颗粒之间的接 触关系，由力-位移接触定律确定每一接触的接 触力.然后根据作用在每个颗粒上的合力和合 力矩由牛顿运动定律确定每一时间步长后每一 颗粒的速度和位置. 2排列过程和几何特性描述 图1两球接触模型 式中，K为法向接触刚度，e,为接触面外法线的 设球形颗粒的直径为d,在宽和高均为20d 单位矢量.切向接触力表示为增量形式，有 的方形盒子内通过随机数产生器产生球心坐 F=F+△F (4 标，并把颗粒逐个投入盒中.若颗粒与已经存在 式中△为单位步长内的切向力增量，表达式 于盒中的颗粒重叠，就将它舍去.随机发生器一 如下： 共产生9×10个球心坐标，在盒内产生240个互 △F=-k△UF (5) 不重叠的颗粒，作为模拟的初始状态，即在模拟 式中△为单位步长内的位移增量，最大切向 初始时刻，所有颗粒均悬浮于盒中，如图2所示. 接触力即为颗粒间的滑动摩擦力，由下式表示： (F)as=Fs tan (6) 式中p为摩擦角.设排列体中任意颗粒a和n个 颗粒接触，则作用在颗粒a上的合力F,为： F=G.++ (7) 式中F,F分别为颗粒a第j个接触的法向接 触力和切向接触力.作用在颗粒a上的合力矩 M为切向接触力与接触半径的差积. M=∑Rt×Fg (8) ● m引 式中，，为颗粒a第j个接触的切向单位矢量， 图2颗粒体的初始状态 1.3运动定律 根据上述离散元原理，编制了计算机程序， 对于排列体中的任意颗粒a,其运动包括平 通过选择合适的时间步长，模拟颗粒的排列过 动和转动，颗粒的平动由颗粒中心的位置坐标 程.最初颗粒自由下落，经过一定时间后，颗粒 x,速度x,和加速度，描述.颗粒的转动由角速 间开始碰撞，产生接触，要保持运算的稳定性， 度8，和角加速度，表示.这样任意颗粒i的运动 时间步长必须很小，每一次计算都至少需要花 方程可表示为： 费几十个小时的CPU时间.为此我们改进了算

. 8 8 3 - 4 北 京 科 技 大 学 学 报 年 第 期 9 9 9 1 系 x o y , 重 力的方 向为 y 轴 的负 向 , 则颗粒 a 的 重 力可 表 示为 : G = 一 m 君 ( l ) 式 中 , m ,为颗粒 的质 量 , g 为重 力加速 度 , i 为张 量 表示法中的 自由指标 , 此式 =1 2 . 若两 颗粒近 邻 , 则颗粒 之 间 的范 德华力 { ’ 0] 可 表 示为 : _ . AR 君 = 若哥子丁 e ( 2) 12 h 之 “ ` 、 “ z 式 中 , A 为 H a m ak e r 常 数 , 对 于 钢 来 说 等于 ZI . l x 1o 一 z0J , R 为 颗粒半径 , h 为两 颗粒表 面之 间 的距 离 , e ,为两 颗粒 中心 连线方 向的单位矢量 . 1 . 2 接触定律 图 1 表 示相 互接触 的两 球 a , b , 其重 叠 量为 砂 , 则 法 向接 触力 尸 ” 由下 式确 定 尸 n = K ” 研 ( 3 ) } 只 一 im{ l城 “ I 叹 ( 9 ) 式 中 I 为颗粒 a 的转动 惯量 , 运用 中心 差 值法 , 可 得方程 (9 ) 的数值 表达式 . 时 间由 t 时刻增加 到 t + △t 时刻 , 则 t + △t 时刻的 速度为 : 、 一 、 一 十 臀 △, 。 虎 )一 、 一+)z 黔 , ( 10 ) t + △t 时刻 的球心位置 坐标为 : x,( `+ 酬 = 对 )t + 戈{ `+ 胡 , + △t ( 1 1) 在模拟 开 始 , 颗 粒之 间没 有接触而 自由下 落 , 一 定时间后 , 颗粒 的速 度达到最 大值 , 颗 粒 之间开 始出现接触 . 此后 , 在每一 个时 间步长 开 始 , 根据 颗粒和 盒 壁的位置 确 定颗粒之 间 的接 触关系 , 由力一位移接触定律确定 每一 接触 的接 触力 . 然 后 根据作 用在每个颗粒 上的合力和 合 力矩 由牛顿运动 定律确定每 一时间步长后 每一 颗 粒的 速度 和 位置 . 图 1 两球接触模型 式 中 , K , 为法 向接触 刚 度 , ie 为接触 面外法线 的 单位 矢量 . 切 向接触 力尸 表示 为增量形 式 , 有 月 , = 尸 ` + △尸 , ( 4 ) 式 中△尸 s 为单位步 长 内的切 向力增量 , 表 达式 如下 : 八尸 ` = 一 k s△鱿 (5 ) 式 中 △梦 为单位步 长 内的位 移增量 . 最大切 向 接触力 即 为颗 粒间 的滑动 摩擦 力 , 由下 式表示 : (尸 , ) m 。 = 尸 t a n 沪 ( 6 ) 式 中 沪为摩擦角 . 设排 列 体 中任 意颗 粒 a 和 n 个 颗粒接触 , 则作用 在颗 粒 a 上 的合力 兀 为 : 只 二 G ; + 艺月 十 艺卿 + 乏月 ’ ( 7) 了= I J二 l j = l 式中 衅 , 月 ’ 分别为 颗粒 a 第 j 个接触 的法 向接 触 力和切 向接触力 . 作用在颗粒 a 上 的合力矩 杯 为切 向接触 力与接触半径 的差 积 . 2 排列过程 和几何特性描述 设球形颗粒 的直 径为 d , 在 宽和 高均 为 20 d 的方形 盒 子 内通 过 随机 数产 生 器 产生 球 心 坐 标 , 并把颗粒逐个投入盒 中 . 若颗粒与 己 经 存在 于 盒 中的颗粒重叠 , 就将它舍去 . 随机发生 器一 共产 生 x9 l 少个球心 坐标 , 在盒 内产 生 2 40 个互 不 重叠 的颗粒 , 作为模拟的初始状态 , 即在模拟 初始时刻 , 所有颗粒均悬浮于 盒 中 , 如 图 2 所示 . 似 = 艺R心x 凡 , ( 8 ) 式 中 , 心为颗 粒 a 第 j 个接触 的切 向单位矢 量 . 1 . 3 运动定律 对于 排列体 中的任意颗粒 a , 其运 动包 括平 动和 转动 . 颗粒 的 平 动 由颗粒 中心 的位置 坐 标 x , 速度 戈 ` 和 加速度 父 ` 描述 . 颗粒 的转动 由 角速 度 以和 角 加速 度反表示 . 这样任 意颗粒 i 的运 动 方程 可 表示 为 : 图 2 颗粒体 的初始 状态 根据上述离散元原理 , 编制 了计算机程序 , 通过选择合 适的时 间步长 , 模拟颗粒 的排列 过 程 . 最初颗 粒 自由下 落 , 经过一 定 时间后 , 颗粒 间开 始碰撞 , 产生接触 . 要保持运算 的稳定性 , 时间步长 必 须 很 小 . 每一 次计算都 至 少需要 花 费几十个小 时的 C P U 时间 . 为此我 们改进 了 算

Vol.21 No.4 程远方等：球形颗粒随机排列过程的计算机模拟 ·389· 法：颗粒发生首次碰撞，令其速度为零，同时为 排和有序排列的混合体，本文未示出).图5是 减少达到稳定状态的时间，使用了局部非粘性 排列密度随时间变化的曲线 阻尼，来耗散颗粒体的能量.考虑阻尼后，方 程(10)改写为： -2="-a+(1士am4/ (12) rg=0+1±a國a 这样可以使时间步长大大增加，使模拟时 间缩短为2~4h.每隔一定迭代次数，将颗粒体 的相关数据存盘，然后通过计算排列密度、配位 数分布函数和平均配位数来表征颗粒体的几何 特征. 为消除盒子边界引起的误差，上述参数的 图31=0.0475s时的状态 计算是在盒子中心的一个小范围内进行的.这 一范围为一假想的小盒子，宽16d、高9d,底部 与原来盒子底部相距2d,左右边界与原盒左右 边界相距1d.上边界位于上部第一层颗粒以下. 排列密度为假想盒内颗粒所占的总面积与假想 盒子总面积之比.首先计算假想盒子内的颗粒 数，然后由此求出盒内颗粒的总面积.颗粒的配 位数是指与研究颗粒接触的相邻颗粒的个数. 假设当两颗粒中心之间的距离小于1.0005d时， 认为它们在接触.对上述假想盒内颗粒的接触 情况进行计算，可得到排列体配位数频度分布 图4t=0.147s时的状态 规律，平均配位数是指各颗粒配位数之和除以 颗粒总数 0.9 0.8 3应用实例 0.7 蒙 0.6 3.1单元系均匀球形颗粒的排列过程 0.5 模拟条件如下：如图2所示，在边长为2mm 0.4 0.000.050.100.150.20 的方盒内，随机放置240个互不重叠的球形颗 t/s 粒，颗粒直径为100μm,颗粒密度为7.86gcm3, 图5排列密度随时间变化曲线 颗粒间及颗粒与盒壁间的摩擦因数为0.364，初 上述模拟结果表明：对于无摩擦光滑球形 始排列密度均为0.4712.模拟开始，颗粒在重力 颗粒，稳定状态为随机密排和有序密排的混合 作用下向盒底下沉，排列密度迅速增加.当颗粒 体.但是对于实际颗粒，颗粒间存在摩擦力和范 首次发生碰撞时，将其速度设为零.此后致密化 德华力，稳态排列密度等于或小于随机密排密 速率使排列逐渐趋于稳定状态，图2~4表示排 度0.822，这与文献报道的实验结果是一致的. 列结构随时间的变化过程，模拟开始，初始排列 颗粒在降落过程中，不仅y坐标值在减小， 相对密度为0.4712：当模拟时间达到0.0475s 由于重排，x坐标值也在变化.图6给出了第10 时，排列密度为0.6894，此时的排列体是不均匀 个颗粒在模拟过程中的运动轨迹，可以看出，最 的（底部密，顶部松）：当时间为0.147s时，排列 初颗粒只有y向位移，为自由落体运动：当发生 达到稳定状态，排列密度为0.8236（利用本模拟 初始碰撞后，颗粒开始有x方向的运动，并且出 方法，若不考虑颗粒间的摩擦力和分子间作用 现往复运动，为颗粒重排过程， 力，其稳定状态的排列密度为0.8696，为随机密

V b l . Z I N o . 4 程 远方 等 : 球 形颗 粒随 机排 列过 程 的计算 机模 拟 . 3 8 9 - 法 : 颗 粒发 生 首次碰 撞 , 令其速度 为零 . 同 时为 减 少达到稳定状态 的 时 间 , 使用 了 局 部非粘 性 阻 尼 , 来耗散颗粒体 的 能 量 〔, 1] , 考虑 阻尼 后 , 方 程 ( 10 ) 改写 为 : 「一 * , 、 . ` , _ ` 消 、 , . 、 尸 t) 。 ! 戈产 `一 “ ` 厂 之 , = 交产 ` 一 创 ` , + ( l 士 a )一 △t l 一 ` 一 m 之 门 2 、 } 已“ 一 ` ’ 一 尽 ` 一引 十 (` 士 a 瑞六△, 这样 可 以使时 间步长 大大增加 , 使模 拟时 间缩短 为 2 一4 h . 每 隔一 定 迭代次数 , 将颗 粒体 的相关数据存盘 , 然后通过计算排列密度 、 配位 数分布 函数和 平均配位数来表征颗粒体 的几何 特 征 . 为消除盒子边 界引起 的误 差 , 上述 参数 的 计算是在盒 子中心 的一 个小范 围内进行 的 . 这 一范 围为一 假想的小盒子 , 宽 16 d 、 高 g d , 底部 与 原来盒子底部相距 Z d , 左右边界与原盒左右 边 界相 距 ld . 上边界 位于上部第一 层颗粒 以下 . 排列密度 为假想盒 内颗粒所 占的总面积 与假想 盒 子总面积之 比 . 首先计算假 想盒子 内的颗粒 数 , 然后 由此求 出盒 内颗粒的总面积 . 颗粒 的配 位数 是指 与研究颗粒接触 的相邻 颗粒 的个数 . 假设 当两颗粒中心 之间的距 离小于 1 . O0 5 d 时 , 认 为它们在接 触 . 对上 述假想盒 内颗粒 的接触 情 况进行计 算 , 可 得到排列体配位 数频度分布 规律 . 平均配位数是 指各颗粒配位 数之 和 除 以 颗粒 总数 . 排和 有序排列 的混 合体 , 本文 未示 出 ) . 图 5 是 排列 密度随 时 间 变化 的 曲线 . 瓢 口n 改〕〔 三 i 图 3 t = 0 . 0 4 7 5 5 时 的状 态 馨图 4 t = .0 1 4 7 5 时的状 态 3 3 . 1 应 单 用 元系 实 均 例 匀球形颗粒 的排列过程 侧不却赢一 厂一 模拟条件如下 : 如 图 2 所示 , 在边长 为 2 ~ 的方盒 内 , 随机放 置 2 40 个互 不 重叠 的球形颗 粒 , 颗粒直径为 10 0 林m , 颗粒密度为 7 . 86 创c m , , 颗 粒间 及颗粒与盒壁 间 的摩擦 因数为 .0 36 4 , 初 始排列密度均为 .0 4 7 1 2 . 模拟开 始 , 颗粒在重 力 作用 下 向盒底下沉 , 排 列密度迅速增 加 . 当颗粒 首 次发生 碰撞 时 , 将其速度设 为零 . 此后致密化 速 率使排列逐渐 趋于 稳定状态 . 图 2一 4 表示排 列结构随 时 间 的变化 过程 . 模拟 开 始 , 初始排列 相对 密度为 .0 4 7 1 2 ; 当 模拟 时 间达 到 .0 0 47 5 5 时 , 排列密度 为 0 . 6 8 9 4 , 此 时 的排列 体是不 均匀 的 ( 底部密 , 顶 部松 ) ; 当 时间为 0 . 14 7 5 时 , 排列 达 到稳定状态 , 排 列密 度 为 0 . 8 2 3 6 ( 利用本模拟 方法 , 若不 考虑颗粒 间 的 摩擦力和 分 子 间作用 力 , 其稳定状态的排列密 度为 0 . 8 6 9 6 , 为 随机密 0 . 0 0 0 . 0 5 0 . 10 0 . 1 5 0 , 2 0 t/ S 图 5 排列 密度 随时 间变化 曲线 上 述 模拟结果 表 明 : 对于 无摩擦 光滑球形 颗粒 , 稳定 状态 为随机密 排和 有序 密 排的 混 合 体 . 但是 对 于 实际颗粒 , 颗 粒 间存在摩擦力 和 范 德华 力 , 稳态排 列密度等于 或 小于 随机密排密 度 0 . 8 2 2 , 这 与文献 报道 的 实验 结果是 一 致 的 19] . 颗粒 在降落过程 中 , 不仅 y 坐 标值在减小 , 由 于 重 排 , x 坐标 值也 在变化 . 图 6 给出 了第 or 个颗粒在模拟过程中的 运动轨迹 , 可 以看 出 , 最 初颗粒 只 有y 向位移 , 为 自由落 体运动 ; 当发 生 初始碰 撞后 , 颗粒 开 始有 x 方 向的运动 , 并且 出 现往 复运动 , 为颗粒 重排过程

·390· 北京科技大学学报 1999年第4期 16 3.2二元系均匀球形颗粒的排列 14 8 仍然使用边长为2mm的方盒作为排列容 器，排列体由a,b两种颗粒混合而成，颗粒直径 12 O 均为100um,a种颗粒的密度为7.86gcm,b种 10 ao 颗粒的密度为8.92gcm',a,b两种颗粒的体积 分数分别为0.3和0.7，初始排列密度为0.4712. 颗粒间及颗粒与盒壁间的摩擦因数为0.364.考 10.1 10.210.3 10.410.5 虑如下4种条件：(1)颗粒间只有法向接触力，颗 x/102μm 粒的运动为平动：(2)颗粒间存在法向接触力、 图6颗粒运动轨迹 切向接触力和摩擦力，摩擦因数为0.364，颗粒 排列体中颗粒的配位数是变化的，颗粒的 的运动为平动加转动：(3)除条件(2)考虑的因 最大配位数为6，而最小配位数为2.图7中曲 素外，还考虑颗粒间的范德华力：(4)摩擦因数 线a是该模拟稳定状态的配位数频度分布规律， 由0.364增至0.727，其他条件同(3).模拟过程 曲线b是Sutherland的研究结果，两者具有很 同上，图8是这4种情况达到稳定状态的结构 好的一致性.但是该模拟的平均配位数为3.609， 体. 而Sutherland D N的结果为4.3l.这表明摩擦力 从图中可以看出：尽管这4种情况的初始 和范德华力使排列的平均配位数减少 状态一致，但排列达到稳定状态后排列体的相 0.45 对密度和第二相的分布规律不同.这一点对于 研究复合材料的烧结和第二相掺杂对烧结的影 0.35 响非常重要.因为在复合材料的制备过程中，第 尔 0.25 二相的连通性（定义为渗愈）对基体粉末的致密 化有明显影响.第二相连通性研究的重点是渗 0.15 愈临界值与临界体积分数的关系.文献[12]认为 0.05 渗愈临界值取决于两种颗粒的直径比，从上述 分析可以看出，颗粒的物理参数（如摩擦因数） 3 4 6 配位数 对渗愈也有明显影响，详细情况另文分析 图7配位数分布规律 (a) (b) (d) 图8条件14的稳定状态(a)条件1，(b)条件2，(c)条件3，(d)条件4 5结论 向接触力外，再考虑摩擦力和范德华力，排列密 度会降低，相对密度为0.8236.对于两元系颗粒 离散元是模拟球形颗粒随机排列演化过程 排列来说，颗粒的物理参数影响排列密度和第 的有效方法，利用计算机程序，得到了多种条件 二相分布规律， 下的均匀球形颗粒的随机排列结构.对于100μ 的颗粒来说，只考虑重力和法向接触力，得到 参考文献 的随机密排结构的相对密度为0.8696，它是随 I Mason G.Radial Distribution Functions from Small Pack- 机排列和有序排列的混合体.除考虑重力和法 ing of Spheres.Nature,1968,217:733

· 3 9 0 . 北 京 科 技 大 学 学 报 1 9 9 9 年 第4 期 9 000 l 6 l 4 妇,nQO2 日沈 0 一b X认 10 . 1 1 0 . 2 10 . 3 10 . 4 1 0 x / 1 0 2 件m 图 6 颗粒运动轨迹 排列体 中颗粒 的配位数 是变化 的 , 颗粒 的 最大配位数 为 6 , 而 最 小配位数 为 2 . 图 7 中曲 线 a 是该模拟稳定状态的配位数频度分 布规律 , 曲线 b 是 S u th e d an d 的研 究结 果 `引 , 两 者具 有很 好的一致性 . 但是该 模拟 的平均配位数为 3 . 6 09 , 而 S u t h e ir an d D N 的结 果为 .4 31 . 这 表 明摩擦力 和 范德 华力使排列 的平 均配位数 减少 . 0 . 4 5 1 — ~ — - 一一 l 2 5 15盯卜 / \ .3 2 二 元系均 匀球形颗粒的排列 仍然使用边 长 为 Z m m 的 方盒作为 排列 容 器 , 排列体 由 a , b 两 种颗粒 混合 而成 , 颗粒直 径 均为 10 0 林m , a 种颗粒 的密度 为 7 . 8 6 g/ e m , , b 种 颗粒 的密度 为 .8 92 留c m , , a , b 两 种颗粒 的体积 分数分别为 .0 3 和 0 . 7 , 初始排列密度为 .0 4 7 1 2 . 颗粒 间及 颗粒与盒壁间的摩擦 因数为 .0 36 4 . 考 虑如下 4 种条件 : (1 ) 颗粒间只 有法 向接触力 , 颗 粒 的运 动为平 动 ; (2) 颗粒 间存在 法 向接触 力 、 切 向接触 力和 摩擦力 , 摩擦因 数为 .0 3 64 , 颗粒 的运动 为平动加 转动 ; ③ 除条件 ( 2) 考虑 的因 素外 , 还考虑颗粒间 的范 德华力 ; (4) 摩擦因 数 由 .0 3 6 4 增 至 .0 72 7 , 其他条件 同 ( 3 ) . 模 拟过 程 同上 , 图 8 是这 4 种情 况达 到稳定状态 的结构 体 . 从 图中可 以看 出 : 尽管这 4 种情况 的初始 状态一 致 , 但排 列达到稳定 状态后 排列体的相 对 密度和 第 二相 的分布规 律不 同 . 这 一 点对于 研 究复 合材料 的烧结和 第二 相掺杂对烧结的影 响非 常 重要 . 因 为 在复合 材料的制 备 过程 中 , 第 二相 的连通性 ( 定义 为渗 愈 ) 对基 体粉末 的致密 化 有 明显 影 响 . 第二 相 连通 性研 究 的重 点是 渗 愈 临 界值与临界 体积 分 数 的关 系 . 文献「121 认 为 渗 愈临 界值取 决于 两 种颗粒 的直 径 比 , 从上 述 分析可 以看 出 , 颗粒 的物理参数 ( 如摩擦 因数 ) 对 渗愈 也有 明 显 影响 , 详细 情况另文 分 析 . 铃壕精尔 3 4 5 配 位数 图 7 配位数分布规律 (a) 磷蒸撇 T T I 卜3 1 卫以 冰勺气r r 丫 ] 书薰舞叭鬓铸书彝店烈刃蒸歹那 ( b ) ( d ) 鬓熬i 5 结 离 论 散元是 模拟 图 球 8 形 条 颗 毋 件 粒 双 l 一 随 4 的 机 袱 稳 排 似 定 列 状 演 黝 态 化 摊 ( a ) 过 条 】 程 掉 件 l , ( 又 b ) 刃 条 公口 件 1 1 2 1 尺 , 电 ( 飞 e )条件 3 , ( d )条件 瞻 的有效方法 , 利用计算机程序 , 得到 了多种条件 下 的均 匀球形 颗粒的随机排列 结构 . 对于 10 0 林 m 的颗粒来说 , 只考虑重力和 法 向接触力 , 得到 的 随机 密排 结构 的相 对密度为 .0 8 6 9 6 , 它 是 随 机 排列和 有序排 列的混合体 . 除考虑重力 和 法 向接触 力外 , 再 考虑摩擦力和 范德华力 , 排列密 度会降低 , 相对密度为 0 . 823 6 . 对于 两元 系颗粒 排列来说 , 颗粒 的物理 参数 影响排列密度和 第 二 相 分布规律 . 参 考 文 献 M a s o n G . R a di a l D i s t r i b u t i o n F u n e t i o n s 行o ln S m a ll P a e k - i n g o f S Ph e r e s . N at u r e , 19 6 8 , 2 17 : 7 3 3

Vol.21 No.4 程远方等：球形颗粒随机排列过程的计算机模拟 ·391· 2 Finney JL.Random Packings and the Structure of Simple 9 Cundall PA.Formulation of a Three-dimensional Distinct Liquids.Proc R Soc,1970,319A:479 Element Model?Part I.Int J Rock Mech,Min Sci Ge- 3 Bordia R K.A Theoretical Analysis of Random Packing omech Abstr,1988,25(3):107 Densities of Mono-sized Spheres in Two and Three Di- 10 Hart R,Cundall P A,Lemos J.Formulation of a Three-di- mensions.Script Metallurgica,1984,18:725 mensional Distinct Element Model?Part II.Int J Rock 4 Sutherland D N.Random Packing of Circles in Plane.J Mech,Min Sci Geomech Abstr,1988,25(3):117 Coll Interface Sci,1977,60:96 11 Lange F F.Effect of Interparticle Potentials on Particle 5 German R M.Particle Packing Characteristics.New Jersey: Packing for Ceramic Processing,Powder &Grains 93. Princeton,1989 Thorm T,ed.Balkema:Rotterdam,1998.187 6 Chen Jianhong.Computer Simulation of 2D Particle Pack- 12 Cundall P A.Distinct Element Models of Rock and Soil ing,Sintering Processes and Microstructure-Properties Re- Structure.In:Analytical and Computational Methods in lations:[PHD Thesis].Alfred:Alfred University,1990 Engineering Rock Mechanics,Brown E T,ed.London: 7 Zheng Jinmin.Computer Simulation of Particle Packing Allen Unwin,1987.129 and Sintering:[PhD Thesis].Alfred:Alfred University, 13 Bouvard D,Lange F F.Relation between Percolation and 1991 Particle Coordination in Binary Powder Mixtures.Acta 8 Cundall P A,Strack O D C.A Distinct Numerical Model metall Mater,1991,39(12):3083 for Granular Assemblies.Geotechnique,1979,29:47 Computer Simulation of Random Packing of Spherical Particles Cheng Yunafang,Guo Shijiun,Lai Heyi Material Science and Engineering School,UST Beijing,Beijing 100083,China ABSTRACT A computer simulation for evolution of random structures of spherical particles in two dimen- sions has been performed by the distinct element method.The forces among the particles considered were gravity,friction force and Van Der Waals Interaction(VDWI).The motion of particles involved translation and rotation.For 100 um particles,if friction and VDWI were not considered,the packing assembly was a mixture of random packing and order packing with packing density of 0.8696,if friction and VDWI were taken into account,the packing density was 0.821 3.To study the procedure of packing,packing densities with different time,particle trajectory and coordination number were examined.The packing structures of binary system with different conditions were also given. KEYWORDS particle random packing;distinct element method;coordination number;packing density

从, 1 . 2 1 N o . 4 程 远方 等 : 球 形颗粒 随机 排 列过 程的计 算机 模拟 · 3 9 1 · 2 Fi n n e y J L . R a n d o m P a e k i n g s an d ht e St ur e utr e o f S im p l e L i qu id s . P ocr R S o e , 19 7 0 , 3 19A : 4 7 9 3 B o dr i a R K . A T h e o er ti c a l nA a ly s i s o f R a n d o m P a e k i n g D e n s i t i e s o f M o n o 一 s i z e d SPh e r e s i n l切 0 an d T hr e e D i - m e n s i o n s . S e ir P t M e at ll u rg l e a , 19 84 , 1 8 : 7 2 5 4 S u ht e r lan d D N . R an d o m P a e ki n g o f C i r e l e s i n P lan e . J C o ll I n t e far e e S e i , 19 7 7 , 6 0 : 9 6 5 G e mr an R M . P art i e l e P a e k i n g C har a e t e r i s t i e s . N e w J e r s e y : P r i n e 以 o n , 1 9 89 6 C h e n J 1 a n h o n g . C o m P ut e r Si m u l at i o n o f 2D P art i e l e P a e k - i n g , Si n t e ir n g P r o e e s s e s an d M i e r o s trU e ut r e 一 P r o P e rt i e s R e - lat i o ns : [ P l刃D T h e s i s ] . A l fer d : A l fer d U n i v e r s iyt , 19 90 7 Z h e n g l i钊m i n . C om P ut e r S而u lat i o n o f p art i e l e p a e k i n g 叨d Si n t e ir gn : [ P hD hT e s i s ] . A l fre d : A l fre d U n i v e r s i ty, 19 9 1 8 C un d a ll P A , S tr a c k 0 D C . A D i st i n e t N u m e r i e a l M o d e l for G anr u lar A s s e m b li e s . G e Oet e hn i q u e , 19 7 9 , 2 9 : 4 7 9 C u n d a l l P A . F o rl l l u l a t i o n o f a T hr e e 一 di m e n s i o n a l D i s t i n e t E l e m e n t M o d e l? P 叭 1 . I n t J R o e k M e e h , M i n S e i & G e - o m e e h A b s tr, 1 9 8 8 , 2 5 ( 3) : 10 7 10 H叭 R , C un d a l l P A , L e m o s J . F o mr u l a t i o n o f a T hr e e 一 d i - m e n s i o n a l D i s ti n e t E l e m e in M o d e l? P a rt 11 . I n t J R o e k M e e h , M i n S e i & G e o m e e h A b s t几 1 9 8 8 , 2 5 (3) : 1 1 7 1 1 L an g e F F . E fe e t o f I n t e rp a rt i e l e P o t e n t i a l s o n P art i e l e P a e k i n g of r C e r am i e P r o e e s s i n g , P o w d e r & G r a i n s 9 3 . T h o mr T, e d . B a l k e m a : R o t e r d am , 1 9 9 8 . 1 8 7 1 2 C un d a l l P A . D i s t i n e t E l e m e n t M o d e l s o f R o e k an d s o il S trU e ut r e . I n : A n a ly ti e a l an d C o m P u at t i o n a l M e ht o d s i n E n g i n e e r i n g oR e k M e c h an i e s , B or wn E 工 e d . L o n d o n : A l l e n & U n w i n , 19 8 7 . 12 9 13 B o u v ar d D , L an g e F .F eR l a t i o n b e wt e e n P e r e o lat i o n a n d P art i e l e C o o r d i n at i o n i n B i n a勺 P o w d e r M ixt u r e s . A c at m e at ll M a t e r, 19 9 1 , 3 9( 1 2 ) : 30 8 3 C o m P u t e r S im ul at i o n o f R an d o m P a e k i n g o f S P h e r i e a l P art i e l e s hC e ng uY n af ng, G u o hS ij i u ,n La i H办 , i M a t e ir ia S e i e cn e an d nE g in e e r in g S e h o l , U S T B e ij in g , B e ij ing 10 0 0 8 3 , C h in a AB S T R A C T A e o mP ut e r s im u l at i o n of r e v o lut i o n o f r an d o m s trU e trU e s o f s Ph e ir c a l Part i c l e s i n tw o d im e n - s ion s h as b e e n Per of n n e d 勿 ht e d i s itn e t e l e m e nt m e ht o d . hT e of r e e s am o n g ht e Part i e l e s e o n s ide r e d w e r e g a v lty, 衍 e ti o n fo r c e an d 、 恤n D e r Wh a l s I ent ar e t i o n (V D机 ) . hT e m ot i o n o f Part i e l e s i n v o l v e d tr a n s lat i o n an d r o at i on . Fo r 10 0脚 Part i e l e s , if 示 e t i o n an d V D机 w e r e n o t C o n s id e r e d , ht e Pa c ik n g a s s e m b l y Wa s a m i x trU e o f r an d om Pa e k l n g an d or d e r Pa e ik n g w iht Pa e k i n g d e n s iyt o f 0 . 8 6 9 6 , if 创 e t i o n an d V D W I w e r e t a k e n int o a e e o u n t , ht e Pac ikn g d e n s iyt w a s 0 . 8 2 1 3 . oT s t u d y ht e Pr o e e d ur e o f Pa e k i n g , Pa e ik n g d e n s it i e s w iht d i fe r ent t im e , Part i e l e tr aj e c t o yr an d e o o dr i n at i o n n um b e r w e r e e x am i n e d . hT e Pa c k i n g s trU e trU e s o f b i n娜 s y s t e m W iht id fe r e n t e o n d it ion s w e r e a l s o g i v e n . K E Y W O R D S Part i e l e arn d o m p a e ik n g : d i s t i n e t e l e m e in m e ht o d : e o o r d i n a t i o n mnU b e r : p a c k i n g d e n s iyt