-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 时间 . 为此我 们改进 了 算