球形物冲击下的颗粒缓冲性能研究(编辑修改稿)

球形物冲击下的颗粒缓冲性能研究(编辑修改稿)内容摘要：

、 球形物密度以及 冲 击速度等因素下进行离散元计算，分析各参数对颗粒系统缓冲性能的影响。 （ 2）设置试验方案，以高速摄像拾取图像并进行 Matlab 数值计算为思路， 在不同颗粒粒径、球形物大小 、 球形物密度以及 冲击速度 等因素下开展了试验研究。 （ 3）分析离散元模拟和试验结果，比较不同情况下盛颗粒容器底部受力以及球形物受力，探究球形物冲击 下颗粒物质的缓冲性能 影响因素，冲击力 变化趋势和冲击力峰值 异同。 东北大学毕业设计（论文） 第 1 章 绪论 4 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 5 第 2 章 冲击荷载下的接触模型理论分析 DEM 离散元 程序 流程 当颗粒在阻尼器中运动时，将所有颗粒都视为弹性小球以简化计算。 整个振动响应过程，被分隔成若干个极小的时间步 t ，在一个时间步内颗粒作匀速直线运动，如果获得第 i 个时间步 it 内颗粒的位置、速度、受力情况，便可由牛顿运动定理求得第 i+1 个时间步 1it 时颗粒的位置、速度、受力情况。 通过不断迭代，最终求得整个阻尼系统在一段时间内的响应。 程序的前处理部分主要工作是文件的读入和初始参数的计算。 本程序的读入文件 为。 在 中，定义了颗粒尺寸、颗粒个数、材料杨氏模量、泊松比、摩擦系数、计算时间步数等必要的程序运行参数。 读入文件后，程序的前处理部分会自动计算出程序循环计算部分的必要参数，如时间步时长 t 、颗粒之间的法向和切向刚度 ,NSKK、整个阻尼系统的背景空间大小等。 除此之外，给计算需要的数据数组开辟内存空间、分配颗粒的初始位置、给背景空间划分网格等步骤也都是在程序的前处理部分完成的。 程序的前处理部分为之后的循环计算部分做 好了铺垫。 循环计算部分主要包括三个大块，分别是：搜索模块、受力模块和运动模块。 搜索模块的主要作用是初步判断两颗粒之间还有颗粒与墙壁之间是否能够接触。 分为 GRID、 SEARCH、 SEARCH_WALL 三个主要子程序。 本程序的搜索模块采用的是基于背景网格的搜索模式，将固定的背景空间划分为一定数量的网格，网格大小为边长小于最小颗粒半径 sr 的 倍，每个网格内只能存在一个颗粒的球心。 如果两颗粒球心处于相近的网格内，便判断其可能接触。 受力模块主要包括 FORCE_SINGLE 和 FORCE_WALL 两个子程序，主要用于计算颗粒与颗粒之间还有颗粒与墙壁之间的相互作用力，计算方式会在 节详细说明。 运动模块主要包括 MOVE 子程序，用于计算颗粒及箱体在力的作用下所做的运动。 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 6 下图为 程序流程图示 ： 图 程序流程图 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 7 颗粒间接触受力 FORCE_SINGLE 子程序 用于计算颗粒间接触时，每个 颗粒 受到的合力、合 力矩。 基本思想为：系统中存在 NLIST 个“接触对”，计算每个“接触对”中的 A、 B 两球相互作用时所受到的力和力矩，然后通过对所有“接触 对”进行循环累加，得到每个小球所受到的合力、合力矩。 图 颗粒间接触受力示意图 两球心连线为法线方向 n ，过接触点且垂直法线的为切线方向  ； 两球相互作用力可分成三种：法向力，切向力和滚动摩阻 [1]。 速度求解均采用通过全局坐标系下三个坐标方向的分量进行，其小球球心速度、角速度以及它们的分量为（如图所示） A Ax Ay Az An AB Bx By Bz Bn Bv v v v v vv v v v v v         r r r r r rr r r r r r； A A x A y A zB B x B y B z          r r r rr r r r 程序采用如下模型，用弹簧、粘壶、耦合器和滑动器来模拟表示颗粒间三类接触作用 [10]，即法向、切向和滚动。 Av xyzAxvAyvAzvAnvAvAxAyAzA B n() 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 8 图 三类接触示意图 一、 法向力计算 1 o v e r la p D O T NnN n nF k c     其中：① 第一部分 nk overlap 为最简单的线性模型，在 overlap 很小的情况下成立，是Hertz 接触理论的一阶展开。 overlap 为两球重叠量，即两球半径之和（ DAB ）与两球心距离（ DS ）的差值。 ② 两球接触的法向有效刚度系数 πn RA RBkERA RB , RA 、 RB 为小球半径， E 为材料杨氏模量。 ③法向阻尼系数 2n n nc z k A M A SS ， AMASS 为两小球平均质量，   22ln lnn e p sz e p s ，eps 为恢复系数，由 input 文件给出。 ④ DOTN 为两球球心速度 Av ， Bv 的法向速度差 An Bnvv。 二、切向力计算 1n n ns s s s sF F f F k u       判 断 ：若 11nns s NFF ，则 11nns s NFF 若 11nns s NFF ，则 1n n ns s s s s s sF F f F k u c u        amp。 其中： （ 1） 切向刚度 snk k RNS ， RNS 为法向、切向刚度比，由 input 文件给出。 （ 2） su 为两球接触点切向相对滑移量，通过其分 量 ,sx sy szu u u   计算而得，以 sxu的计算为例进行说明： () () 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 9 DUTsxut   DUT 为 A、 B 球接触点切向方向速度差在 x 轴上的投影 Ax B xvv 程序中的公式： DUT=DUDOTN*XM+VRIATXVRIBTX DU： Av ， Bv 速度差在 x 轴的投影 Ax Bxvv DOTN： Av ， Bv 的法向速度差 An Bnvv DOTN*XM： Av Bv 法向速度差在 x 轴上的投影 Anx Bnxvv DUDOTN*XM： A、 B 球心平动引起的切向速度差在 x 轴上的投影，由总量差减去法向分量差而得， ()A x B x A x B x A n x B n xv v v v v v     VRIATX： A 球转动引起的切向速度在 x 轴的投影 VRIBTX： B 球转动引起的切向速度在 x 轴的投影 VRIATXVRIBTX：由转动引起的 A、 B 球切向速度差在 x 轴的投影。 （两球接触点的相对速度根据两球在接触点的绝对速度的差值求得，每个球在接触点的绝对速度由随着球心平动和绕着球心转动这两个速度矢量叠加而得。 球心（平动）速度对法向和切向速度均有贡献，而绕着球心转动引起的速度只对接触点的切线方向有贡献。 将  分解成三个分量，分别计算三个转动分量引起的速度在三个坐标轴上的投影，然后进行叠加。 x 不产生 x 方向的速度分量，以此类推）。 （ 3） 静滑动摩擦系数 s 在本子程序内部定义并给出数值。 （ 4） 切向阻尼系数 2s s nc z k A M A S S， snz rzs z， rzs 法向、切向阻尼比，由 input给出。 （ 5） su 也通过其分量计算而得，例如： sxu DUT。 三、 滚动摩阻计算： 滚动摩擦模型： m in ( , )r r r r nM k c r F   判断： 若 11nnr r NM r F ，则 11nnr r NM r F ，（此处 1nNF 不考虑粘滞力，与上同）。 其中： （ 1） 滚动摩阻的刚度系数 22rsR A R BkkR A R B  () () 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 10 （ 2） 静滚动摩阻系数   22 2 2224rD S R B R Ar A R F R BDS   ， DS 为两球球心距离，ARF 未知。 颗粒与边界间接触受力 Force_WALL 子程序 计算颗粒与边界接触时，边界受到的合力及颗粒受到的合力、合力距。 由于圆柱的特殊性，可把圆柱桶分成三个墙壁（ NWall 3 ），即下壁（ IW1 ）、圆 柱壁（ IW2 ）以及上壁（ IW3 ） . 接触判断：一个球最多可能和三个墙壁同时接触  ip = na_w ic , ic 为 接触对粒子编号。  iw = nb_w ic *, ic* 接触对墙编号。 粒子球心坐标  x(ip), y(ip), z(ip)， 接触点的坐标  x _ c o n t a c t , y _ c o n t a c t , z _ c o n t a c t. 图 颗粒与边界接触示意图 Fortran 程序中 上下墙壁 ， 如果 IW1 : x_contact = x(ip) y_contact = y(ip) z_contact = dx = x_contact x(ip)=0 dy = y_contact y(ip)=0 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 11 dz = z_contact z(ip) ds = sqrt(dx*dx + dy*dy + dz*dz) 如果 IW=3 ： x_contact = x(ip) y_contact = y(ip) z_contact = zPoint(1) （ zPoint(1)为上墙壁的 z 坐标） dx = x_contact x(ip)=0 dy = y_contact y(ip)=0 dz = z_contact z(ip) ds = sqrt(dx*dx + dy*dy + dz*dz) 如果 IW2 ： 图 圆柱壁受力示意图 R_Location = sqrt( x(ip)*x(ip) + y(ip)*y(ip) ) （粒子球心距原点距离） R_Contact = R_Cylinder （接触点距原点距离） x_contact = x(ip)*R_Contact/R_Location y_contact = y(ip)*R_Contact/R_Location z_Contact = z(ip) dx = x_contact x(ip) dy = y_contact y(ip) 东北大学毕业设计（论文） 第 2章 冲击载荷下的接触模型理论分析 12 dz = z_Contact z(ip)=0 ds = sqrt(dx*dx + dy*dy + dz*dz) 当 ds R_Particle 且 ds0 时，粒子和壁之间有接触，此时粒子与墙之间的重叠量Overlap = R_Particle – ds。 颗粒碰撞后运动状态 MOVE 子程序 用于计算颗粒碰撞之后的运动。 得到颗粒在碰撞之后的位置和运动参数（包括颗粒的速度 U(I)、 V(I)、 W(I)，坐标 X(I)、 Y(I)、 Z(I)，角速度 OM1(I)、 OM2(I)、OM3(I)，转角 THX(I)、 THY(I)、 THZ(I)， 将这些参数传递回 中的 mon公共区中，再提供给下一个时间步，重新计算颗粒碰撞产生的力和力矩，从而形成颗粒的运动和受力之间的循环模式。 在程序中，颗粒运动参数的计算方法采用如下中心差分格式 ， 用 Nt 时刻的加速度 /角加速度和12Nt时刻的速度 /角速度来更新12Nt时刻的速度 /角速度；用12Nt时刻的速度 /角速度和 Nt 时刻的位移 /转角更新 1Nt 时刻的位移 /转角。 计算单个颗粒的质量 ： 3432D S I Z EA M A S S D E N S   其中： AMASS 为 颗粒质量； DSIZE。