(优选)离散元法及其应用课件

1、（优选）离散元法及其应用第1页，共207页。一、引言 在自然界和工农业生产领域，大量存在着颗粒材料，如农产品、肥料、土壤、药品、煤炭和岩石等。据估计世界上50 % 的产品和75 % 的原材料都是颗粒材料。 在农业生产领域，耕地、开沟、播种、施肥、镇压、脱粒、分离、清选、粉碎、干燥、输送、仓储、分级、加工和包装等过程中，始终存在着颗粒材料与农机部件的接触作用和颗粒材料的流动过程。 在众多工业生产领域，如制药、食品、化工、冶金、采矿、能源、岩土工程等领域，也大量存在着颗粒材料与机械部件的接触作用和颗粒材料的流动过程。第2页，共207页。 一个好的农机部件设计应使: 颗粒材料按照预期的方式运动，如播

2、种时种子流动， 减少流动过程中不必要的损伤，如播种时种子损伤， 节省动力消耗，如开沟和耕翻土壤时牵引动力消耗，等等，此时必须考虑机械部件与颗粒材料的接触作用及颗粒群体动力学问题 。 机械部件的优化需考虑颗粒动力学问题。VV播种开沟第3页，共207页。 颗粒材料的性质介于固体与流体之间，又称第四种物质形态，有着复杂的力学特性：非均匀尺寸偏析，如巴西果、反巴西果和三明治效应；粮仓效应；成拱现象；漏斗现象；自组织临界，等等。 自组织临界是Bak等1987年解释非线性复杂系统无序行为时提出的，即大的相互作用系统包含着众多短程相互作用的组元，系统自然地从随机状态演化到一种有序的临界状态，在该状态时小事件

3、引起的连锁反应能够对系统中任何数目的组元产生影响，从而可能导致大规模事件的发生。非均匀尺寸偏析粮仓效应成拱固体颗粒第4页，共207页。 颗粒材料通常指直径大于1m的颗粒组成。依据固体颗粒的浓度，可将颗粒材料分为密相颗粒材料、松散颗粒材料与稀薄颗粒材料。颗粒材料又可分为干颗粒材料(不含液体)与湿颗粒材料(含液体)。 颗粒材料流动可分为：准静态流动，流动的初始阶段，当颗粒承受的载荷超过颗粒间静摩擦力时，颗粒间仍保持接触但开始流动；快流，流动完全发展阶段的快速剪切流动；慢流，处于准静态流动和快流的中间阶段。密相松散稀相第5页，共207页。 二十世纪70年代后，许多物理学家、力学家和应用数学家开始对颗

4、粒运动的物理机制发生兴趣，建立了两类颗粒动力学理论：基于连续介质力学的理论，如颗粒动理论、摩擦塑性模型和光滑粒子法等；基于离散介质力学的理论，如硬颗粒模型、软颗粒模型和Monte Carlo方法等。 连续介质力学理论是把物质或其特性，假设成无论在时间还是在空间位置上，均是连续的或可用连续函数表示。因此物质可以无限分割而不失去其固有特性，不考虑粒子的特性，是描述物质整体及其特性的一种方法。 连续介质离散介质第6页，共207页。1. 基于连续介质力学的理论颗粒动理论(kinetic theory) 研究发现快速颗粒流中单个颗粒的运动，与气体中的分子热运动非常相似。因此，借鉴非均匀的稠密气体分子运动

5、理论，Ogawa定义了颗粒温度，Jenkins将气体的动理论扩展到颗粒材料，在考虑颗粒碰撞及摩擦所造成的能量损失的基础上，修正了Boltzmann方程，得到宏观的颗粒相输运方程，并导出动理论模型，由此可求得固体体积分数分布、颗粒速度分布和浓度分布等。 适合于稀薄颗粒的快流分析。第7页，共207页。 17世纪中期法国工程师Coulomb提出了土的抗剪强度和土压力滑动理论，其后被推广为散体极限破坏的Mohr-Coulomb准则，在此基础上发展成为土力学。 摩擦塑性模型，即是将Mohr-Coulomb准则应用于颗粒材料，当颗粒间载荷超过颗粒间的摩擦结合力，颗粒间开始滑移即屈服，但颗粒仍保持接触并相互

6、摩擦。 人们已建立多种颗粒材料屈服条件，其中有双剪切模型、塑性势模型和双滑移自由转动模型等。 摩擦塑性模型主要应用于准静态颗粒流。摩擦塑性模型F第8页，共207页。 在小变形的情况下，可采用有限元法分析准静态的颗粒运动。此时，颗粒中心为单元节点。由接触建立节点间的联系，通过作用在节点上的力建立平衡方程。 有限元法适合模拟颗粒接触的拓扑结构不发生变化的静态颗粒系统，在动态和大变形的情况下，大量接触的丢失或产生，导致拓扑结构发生很大变化，这将需要耗费大量的时间重新生成单元，其缺点是：网格重构；网格变形较大还将产生计算不收敛；缺少合适的分析模型，如接触和变形模型、分离破裂模型等。有限元方法第9页，共

7、207页。 光滑粒子法思想是(SPH)：通过带质量的粒子离散计算域，粒子即代表颗粒，通过引入表征节点及其影响域内物理量间的关系核函数，来构造局部光滑的连续场，通过求解描述连续场对时间变化规律的常微分方程，来实现数值模拟。 光滑粒子法，在求解爆炸冲击及大变形问题等方面有不少应用，土壤切削等。无网格方法(meshfree method)第10页，共207页。 连续介质理论的基本控制方程是连续方程、动量方程和能量守恒方程。由于颗粒介质并不满足连续性的假定，并且由于连续介质模型没有考虑颗粒物性参数、粒径形状、大小及其分布等对颗粒流的影响，因此用连续介质模型分析颗粒流一般误差较大。 目前进行相关农机部件

8、设计时，大都依靠经验和试验方法，既费时费力又得不到理想的效果。据估计仅由颗粒材料输送所造成的相关设备利用损失就达40%，远未达到优化设计和节省能源的要求。 连续介质离散介质第11页，共207页。 随着计算机技术的发展，基于离散介质力学的理论，愈来愈引起人们的重视。 离散介质力学方法的思想源于较早的分子动力学，适用于模拟颗粒群体的接触或碰撞过程，它的出现补充了连续力学方法的不足。2. 基于离散介质力学的理论第12页，共207页。 1985年Campbell提出硬颗粒模型，其思想是当颗粒表面承受的应力较低时，颗粒不产生显著的塑性变形，碰撞只在瞬间发生，在碰撞过程中颗粒本身不变形，并且只考虑两个颗粒

9、的同时碰撞，而不计三个以上颗粒的同时碰撞，采用动量守恒或能量守恒计算碰撞后颗粒的速度和位置，广泛的应用于快速、低浓度颗粒流的模拟。硬颗粒模型(hard/rigid sphere model)第13页，共207页。 软颗粒模型又称为离散元法。 1971年Cundall提出适于岩石力学的离散元法(discrete /distinct element method，DEM)，1979年Cundall又提出适于土力学的离散元法，并推出二维圆盘程序BALL和三维圆球程序TRUBAL，后发展成商业软件PFC-2D/3D，形成较系统的模型与方法。软颗粒模型(soft sphere model)块体模型颗粒模

10、型第14页，共207页。 离散元法的基本思想是，把散粒群体简化成具有一定形状和质量颗粒的集合，赋予接触颗粒间及颗粒与接触边界(机械部件)间某种接触力学模型和模型中的参数，以考虑颗粒之间及颗粒与边界间的接触作用和散粒体与边界的不同物理机械性质。3. 离散元法的基本方法第15页，共207页。 离散元法的基本假设：单元是刚性的，即单元的几何形状不会因单元间的挤压力作用而改变；由于计算时步间隔取得足够小，单元的速度和加速度在一个时步内为常量，并且单元在一个时步内只能以很小的位移与其相邻单元作用，其作用力也只能传递到其邻接单元，而不能传递得更远；单元间的连接是靠相互接触实现的，圆形单元的接触为点接触等。

11、 R1R2第16页，共207页。ij 以i 和j 颗粒接触为例，设其法向叠合量为 ，由此产生的法向接触作用力 可如下计算（局部坐标胡克定律）式中 为接触的法向刚度系数。 由于切向接触作用力与运动和加载历史有关，因此切向力通常采用增量形式计算，t 时刻的切向力 为式中 为上一时步接触的切向作用力； 为接触的切向刚度系数； 为接触点的切向相对位移； 为计算时步。R1R2ij（静摩擦力）n R1R2ijFn（斥力）第17页，共207页。ij 求i 颗粒质心作用的合力和合力矩为（全局坐标系下） 求i 颗粒的新位置有二种方法，静态松弛法和动态松弛法。第18页，共207页。ijij 静态松弛法根据颗粒不平

12、衡力达到再平衡时的力与位移关系建立平衡方程组，通过求解方程组得到颗粒的新位置，是一种隐式解法且需求解刚度矩阵。 动态松弛法采用牛顿第二定律求解颗粒的新位置如下。第19页，共207页。 动态松弛法需加入人工阻尼，使方程解收敛。因此，阻尼系数和时步选择对求解精度和收敛速度影响较大。动态松弛法牛顿第二定律（球颗粒）：第20页，共207页。 遍历所有颗粒，然后进入下一时步； 显式解法，适合于求解非线性问题。动态松弛法中心差分（数值积分）：第21页，共207页。 由于其离散的特点，在分析高度复杂的系统时，无论是颗粒还是边界均不需作大的简化；当赋予接触颗粒间不同的接触模型时，还可以分析颗粒结块、颗粒群聚合

13、体的破碎过程、多相流动甚至可以包括化学反应和传热问题。 正是由于诸多优点，使得离散元法已成为研究颗粒群体动力学问题的通用方法，并在岩土工程及装备、采矿工程及装备、化工过程及装备、制药工程及装备、食品工程及装备和农业工程及装备等研究领域得到较多应用。第22页，共207页。4. 离散元法的应用图1 边坡稳定性分析图2 块体拱的稳定性分析图3 地下洞室稳定性分析无 粘干颗粒接触作用力为斥力FnFn第23页，共207页。图4 颗粒沉降和输送过程试验与仿真对比(a)试验(b)仿真第24页，共207页。图5 料仓落料过程试验与仿真对比(a) 试验(b) 仿真图6 料仓落料过程试验与仿真对比(a) 试验(b

14、) 仿真第25页，共207页。图7 料仓落料过程模拟图8 颗粒材料在平面的堆积过程模拟图9 颗粒材料由料仓的出流过程模拟(a) 圆颗粒(b) 长扁颗粒(c) 圆柱颗粒第26页，共207页。图10 材料压实过程模拟图12 颗粒冲击过程模拟与试验对比(a) 冲击前(b) 冲击过程模拟(c) 冲击过程试验图11 地基的夯实过程分析第27页，共207页。图13 香蕉筛筛分过程分析图14 香蕉筛筛分过程分析图15 糙米的筛分过程模拟 (a) 模拟(b) 试验第28页，共207页。模拟试验图17 滚筒工作过程试验与模拟比较图16 不同尺寸的筛分过程分析第29页，共207页。图18 球磨机工作过程试验与模拟

15、比较(b) 模拟(a) 试验实际球磨机第30页，共207页。图19 滚筒式混料机的工作过程分析(a)(b)(c)(d)第31页，共207页。图20 一种混料机工作过程的试验和模拟分析比较(b) 模拟(a) 试验第32页，共207页。图21 混料过程的试验和模拟分析比较(b) 模拟(a) 试验图22 混合机的工作过程模拟(a) 试验(b) 模拟第33页，共207页。图23 滚筒式混合机的工作过程模拟第34页，共207页。图24 滚筒式混合机的工作过程模拟(a) (b)(c)(d) (e) (f)第35页，共207页。图27 苹果运输时损伤过程模拟图26 颗粒材料螺旋混合过程模拟图25 散粒物料的

16、混合过程分析第36页，共207页。图28 颗粒材料螺旋输送过程模拟(a)(b)(c)图29 颗粒材料带式输送过程模拟图30大豆在斜槽中流动过程分析第37页，共207页。(a) 自然堆积(b) 墙角堆积图31 颗粒材料在平面的堆积过程模拟图32 颗粒材料在容器内的堆积过程分析(a) 圆颗粒(b) 多边形颗粒第38页，共207页。图33 颗粒在容器内的堆积过程模拟图34 颗粒材料在容器内的堆积图35 颗粒材料在容器内的堆积第39页，共207页。图36 颗粒在容器内的堆积过程模拟图37 肥料撒施过程模拟第40页，共207页。图38 推土铲的工作过程模拟(a)(b)(c)(d)(e)第41页，共207

17、页。图39 铲装工作过程的试验与仿真对比(a)试验(b)仿真第42页，共207页。图40 铲装过程模拟接触作用力为斥力无 粘干颗粒FnFn第43页，共207页。图42 圆锥体插入土壤过程模拟图41 土壤和轮胎的相互作用分析FEM-DEM接触作用力为斥力Fn接触作用力为吸力FnFn湿颗粒粘颗粒+FnFnFn第44页，共207页。图43 金属杆插入土壤时的阻力和变形模拟和试验比较图44 摆式铲的切土过程试验与模拟比较(a) 模拟(b) 试验(b) 模拟(a) 试验第45页，共207页。图46 车轮的通过性分析模拟图45 一种带刺轮滚动过程模拟(a) 平地面(b) 坡路第46页，共207页。图47

18、铲装过程模拟第47页，共207页。图48 轮胎与土壤的作用过程模拟分析图49 犁耕地过程模拟分析第48页，共207页。图51 轨道与地基的作用过程分析图50 机械部件与土壤的作用过程分析第49页，共207页。图52 滑坡过程分析第50页，共207页。图53 脆性材料剪切破坏过程模拟 材料的破坏过程，实质是力学模型从连续介质模型到离散介质模型的转变过程。 斥力接触模型FnFnFnFn+吸力连接模型颗 粒聚合体第51页，共207页。图54 复合材料受压破坏过程模拟图55 地震时高坝肩的破坏过程模拟(a)(b)(c)第52页，共207页。图57 烟囱定向爆破过程模拟图56 地震时桥墩倒塌破坏过程分析

19、(a)(b)(c)(a)(b)(c)第53页，共207页。图58 混凝土结构框架的破坏过程模拟图59 海冰撞击钻井平台过程模拟图60 破冰船的破冰过程模拟第54页，共207页。图61 脆性材料剪切破坏过程试验与仿真比较(b) 仿真(a) 试验第55页，共207页。图62 钢弹冲击脆性材料过程模拟图63 导弹冲击脆性材料过程模拟(a) 速度300米/秒(b) 速度50米/秒(a) 速度300米/秒(b) 速度50米/秒第56页，共207页。图64 导弹冲击脆性材料过程模拟(a) 加强前(b) 加强后第57页，共207页。图65 玻璃冲击过程模拟(a) 单层(b) 夹层第58页，共207页。图66

20、 颗粒结块及冲击图67 结块冲击破裂过程模拟图68 颗粒的破碎过程模拟第59页，共207页。图69 大钢球冲击过程模拟图70 颗粒的破裂过程模拟第60页，共207页。图71 Evolution of the coating process颗粒的包衣、磨损、腐蚀过程模拟第61页，共207页。图72 脆性材料的切削过程模拟分析图73 磨损过程模拟分析Periodical boundariesUpper wallLower wallAbrasive particles图74 岩石磨削过程模拟分析第62页，共207页。图75 岩石切削过程模拟分析图76 颗粒冲击叶轮过程模拟分析图77 颗粒驱动水轮过程

21、模拟分析第63页，共207页。图78 冲击过程模拟分析第64页，共207页。图79 自然界风沙流动过程模拟图80 颗粒材料水中沉降过程模拟接触模型 +考虑流体力作用力CFD模型图81 气力清选过程模拟(a) 清选机结构和原理(b) 清选过程模拟第65页，共207页。图82 气吹式排种器工作过程模拟图83 气力输送过程模拟接触模型 +考虑流体力作用力CFD模型第66页，共207页。图84 颗粒的燃烧过程模拟图85 颗粒的干燥过程模拟（化学反应）接触模型 +化学反应过程化学反应模型干燥过程传热传质模型其他外力如电磁力等模型第67页，共207页。 1.颗粒的离散元法建模方法 在采用离散元法分析颗粒材

22、料的流动过程时，把颗粒简化成何种几何模型，如何给定几何模型中的参数建立其分析模型，将直接影响到分析精度和计算效率。二、研究现状及存在问题 颗粒模型块体模型第68页，共207页。已报道的颗粒建模方法 超二次方程 统一建模方法自然界中70%以上的颗粒n1,n2=1.5n1,n2=0.5n1,n2=220202211第69页，共207页。 组合颗粒模型统一计算方法第70页，共207页。颗粒建模非球颗粒球颗粒超二次方程建模规则形状大豆，超二次方程建模非规则形状玉米、岩石，组合颗粒建模玉米籽粒大豆籽粒 规则形状颗粒组合充填建模第71页，共207页。非规则形状颗粒充填组合建模方法非规则形状颗粒(b) 三点

23、成球方法(a) 一点相切方法(c) 离散元方法(d) 最大覆盖方法(e) 松弛方法冗余数据点去除充填过程优化质心和转动惯量计算尺寸减小消除叠合第72页，共207页。 已报道的三种农作物颗粒模型(a) 大豆籽粒模型(b) 玉米籽粒模型(c) 玉米籽粒模型 我们建立的颗粒模型（组合颗粒模型）(a) 大豆种子模型(c) 小麦种子模型(b) 玉米种子模型采用人机交互充填建模基于点云的充填方法基于三角形网格面的充填方法自动充填(d) 稻米籽粒模型第73页，共207页。(e) 小麦果穗模型(d) 玉米果穗模型(f) 小麦植株模型 籽粒采用人机交互充填或自动充填建模 我们建立的颗粒模型+=+=小麦果穗小麦茎

24、秆小麦叶片穗轴节片籽粒果柄基于颗粒聚合体方法质点弹簧模型果穗采用建模第74页，共207页。颗粒材料的样本生成方法(1)随机生成法；(2)三角划分法。加密方法：(1)重力的方法；(2)加压和尺寸扩张方法。分析计算区域P尺寸增大消除间隙加压消除间隙入料口第75页，共207页。 (1)当采用超二次方程建立颗粒模型时，如何高效、高精度地进行颗粒与颗粒之间、颗粒与边界之间的接触判定，还需要深入研究。 (2)当采用球组合方法建立颗粒模型时，如何建立一种非规则形状颗粒的组合方法，还需要深入研究。 (3)由于农业生产中的散粒物料种类繁多、形状复杂，而且即使是同一种散粒物料，其颗粒形状和尺寸参数差别也较大。因而

25、在实际分析设计时，把散粒物料颗粒简化成何种几模型，如何给定几何模型中参数，是一个需要深入研究的问题，壤土、玉米植株、水稻等。n第76页，共207页。已报道的的边界建模方法(1)函数建模 简单边界的建模如一个圆柱面 国外著名商品软件PFC、 EDEM(2) 颗粒排列方法 简单边界的建模如一个平面(3) Kremmer等提出“有限壁”方法 复杂边界建模方法 国外著名商品软件EDEM2. 边界的离散元法建模方法 边界边界第77页，共207页。 “有限壁”方法是一种近似方法：求解颗粒与边界的接触点（接触力作用点）和接触叠合量（用于计算接触作用力）是近似的；在小三角形平面的相互连接处一阶导数不连续，致使

26、所求颗粒与边界的法向和切向接触作用力具有突变性；在采用赫芝模型求解接触作用力时，曲面边界取无穷大的曲率半径，这些均与实际的边界情况差别较大。为了提高建模精度，往往还需要较多数量的三角形平面，由此增加了离散元法的计算时间。FnFn第78页，共207页。基于CAD模型的边界建模方法 边界建模非规则曲面不能用初等解析数函表示，离散成规则曲面 规则曲面平面、球面、 柱面、锥面等，可用初等解析函 数表示，当曲面上存在缺失部分时，提出实边 界和虚边界的方法虚边界实边界球颗粒 曲面上划分网格常用两种方法，一种是映射法另一种是直接法。采用了直接法中的推进波前法AFT（Advancing Front Techn

27、ique）。物理域 参数域 波前法平面边界球面柱面椭球面第79页，共207页。运动边界建模VV分离筛筛面运动形式复杂运动一种联合收割机中分离筛等简单运动平动、转动和平转动组合平动、转动复杂运动播种单体开沟播种铲车履带车第80页，共207页。nij边界n平动 新位置 原位置+Vt转动 新位置 原位置+ r t平转动边界运动建模第81页，共207页。复杂运动边界运动建模分离筛筛面运动形式多体动力学方法图解方法分析方法笛卡尔方法一般机械系统 拉格朗日方法 航天器系统多体动力学方法及其软件的特点：自动建模、自动求解为笛卡尔坐标阵的约束方程 为各刚体的质量矩阵 为约束方程的雅可比矩阵 笛卡尔位形坐标阵

28、为各刚体上受到的外力矩阵为与约束方程对应的拉格朗日乘子阵 第82页，共207页。振动筛筛面质心位移的变化曲线 （a）X方向（b）Y方向分离筛筛面 我们采用多刚体运动学中的笛卡尔方法，分别建立各种铰链约束方程，采用牛顿拉斐逊方法求解系统的约束方程、速度与加速度方程，得到每个时步、每个刚体的平动位移和速度、角位移和角速度等，实现了多刚体动力学（MBD）与三维离散元法（DEM）的耦合。 第83页，共207页。振动筛筛面质心加速度的变化曲线 （a）X方向（b）Y方向分离筛筛面振动筛筛面质心速度的变化曲线 （a）X方向（b）Y方向第84页，共207页。基于DEM和MBD耦合的筛分过程仿真分析(a) 二维

29、分析(b) 三维分析第85页，共207页。 离散元法发展到今天，大部分研究都集中在颗粒的分析模型和接触力学模型等方面，对于边界建模的讨论还较少。 建立一种高效、精确和适应性广的离散元法边界建模方法，已成为离散元法实际应用急需解决的关键问题之一。第86页，共207页。 (1)如何建立一种非规则曲面边界的离散元法分析模型还需要深入研究，离散方法离散成三角形平面或球面或其它曲面；精度控制基于几何自适应；与CAD软件集成；自动建模等。 (2)如何建立复杂运动规律边界的离散元法分析模型还需要深入研究，非平动转动边界。 (3)如何建立弹性边界的离散元法分析模型还需要深入研究，弹性边界或橡胶材料边界等。输送

30、带弹簧第87页，共207页。 在离散元法分析中，颗粒与颗粒之间、颗粒与边界之间的相互作用力，一般分为法向和切向两个方向，并分别采用不同的力学模型来计算。由于通常要分析计算的颗粒数量较多，采用较完备的力学理论来建立模型并由此计算作用力，不仅计算复杂，而且计算量较大，因此，通常都采用简化模型，常用的接触力学模型有以下几种。 3. 相互作用的力学模型及参数确定第88页，共207页。 线性粘弹性模型，法向接触作用力为（斥力） 式中Fn为接触两体间的法向作用力；kn为接触的法向刚度系数，n为接触两体的法向叠合量；cn为法向粘性阻尼系数，vn为两体接触处的法向相对速度。3.1 线性粘弹性接触力学模型 Fn

31、(N)KNKYKZCNCZCYYXZ单元1单元2n R1R2ij（无粘干颗粒）第89页，共207页。 由于切向力的大小与加载历史有关，因而通常切向接触力的计算都采用增量形式，切向作用力为（静摩擦力） 式中Fs(t)为当前时步接触两体间切向作用力；Fsk(t)为当前时步接触两体间切向弹性力， ； 为上一时步接触两体间切向弹性力，t为计算时间步长；ks为接触的切向刚度系数，vs为接触处的切向相对速度； 为当前时步接触两体间切向阻尼力， ；cs为切向粘性阻尼系数。s=x, y，线性粘弹性模型适合于无粘干颗粒。R1R2ij第90页，共207页。 线性粘弹性模型虽然得到广泛应用，但实际上颗粒相互接触时，

32、法向作用力都是非线性的，非线性粘弹性模型可由赫兹弹性接触理论得到为（斥力） 式中 ， ， ， R1、R2分别为接触两体接触处的曲率半径，E1、E2分别为接触两体的弹性模量， 、 分别为接触两体的泊松比。3.2 非线性粘弹性力学模型 Fn(N)n R1R2ij（无粘干颗粒）第91页，共207页。 切向接触力学模型为（静摩擦力）s=x, y，非线性粘弹性模型适合于无粘干颗粒。R1R2ij第92页，共207页。 人们深入研究发现，在上述粘弹性模型中，当接触的法向叠合量n=0时，法向接触力Fn0，很显然这与实际情况相违背；而且实际接触的两颗粒间，粘性阻尼力也是非线性的。 为了克服这些缺点，人们又提出了

33、非线性粘弹性模型的一般表达式为（斥力） 式中r可在12间取值；s可在01间取值。Mishra通过实验确定，对于钢球类颗粒，r=1.6、s=0.8比较合适；而在Tsuji等研究中，取r= 1.5、s=0.25。 n R1R2ij第93页，共207页。 (1)粘性阻尼力并非真正代表一种能量耗散机制，而只是使计算尽快稳定的方法，因而cn的确定具有较大随意性； (2)虽然cn可由e求得，但由于e不仅与材料特性有关，还与接触表面形状、接触碰撞速度等有关，这给正确地确定cn也带来难度； (3)当法向作用力较大时，在接触点处产生塑性变形，塑性变形将产生能量耗散，而此时e又较小，根据粘弹性模型计算的法向力已有

34、较大误差。3.3 弹塑性接触力学模型 第94页，共207页。 为此，Mishra、Walton、 Vu-Quoc等又分别提出了几种弹塑性法向接触力学模型，其中Walton提出的半锁弹簧(又称双线性)模型为（斥力） 式中k1和k2分别为加载和卸载时的法向刚度系数；和0分别为接触两体的法向叠合量和残余法向叠合量，k1和k2还满足关系式 。加载时卸载时n R1R2ij第95页，共207页。 切向粘弹性模型也存在着正确确定cs较难等缺点。为此，一些研究人员又根据Mindlin的弹性摩擦接触理论，并由不同简化得到几种不同的力学模型，其一般表达式为 式中切向刚度系数ks是一个变化值，需在每一时步不断重新计

35、算，而ks的计算方法不同，力学模型则不同。Walton等提出的切向力学模型，其ks的计算式为 式中当接触的法向作用力由半锁弹簧模型计算时， ；为摩擦系数；Fn(t)为t时刻的法向作用力；Fs*为Fs曲线的转折点，图中切向位移S=vst。 Fs增大时Fs减小时第96页，共207页。 上述Walton的切向力模型，只适用于Fn为常值时切向接触力的计算，但在多数情况下Fn是变化的，为此Vu-Quoc等又提出另一种切向力计算模型，其ks计算式为 式中当接触的法向作用力由半锁弹簧模型计算时， ；当接触的法向作用力由赫兹弹性接触的非线性粘弹性模型计算时， 。 Fs增大且|Fs|Fs*时Fs增大且|Fs|F

36、s*时Fs减小且|Fs|Fs*时Fs减小且|Fs|Fs*时第97页，共207页。 Vemuri等还提出一种直接计算总切向力的模型，以代替上述增量模型，其计算公式为（静摩擦力） 式中Fs0为Fs的初值；S为某时刻的切向位移，S0为S的初值。（无粘干颗粒）第98页，共207页。 3.4 湿颗粒(土壤)接触力学模型 当两颗粒中心距D R1+R2时，法向接触力由线性粘弹性模型计算为（斥力） 当两颗粒中心距D在 R1+R2 D (1+Cad)( R1+R2)间时，法向接触力为（吸力）n R1R2ij-n R1R2ij第99页，共207页。对上述模型进行了改进，当n0时 当 Ca n 0时nn（吸力）（吸

37、力）V第100页，共207页。3.5 湿颗粒液桥接触力学模型 当两颗粒间间距n (=s1+s2)在0 n(1+0.5 )V1/3时，法向接触力为（吸力）式中 ， ；Ri为颗粒半径(m)；S为颗粒间距(m)；1 为液桥半径(m)；2为液桥颈部半径(m)；i为接触半角 (rad)；i为半填充角(rad)；液体粘度Pa s；为流体表面张力(N/m)； vn为法向相对速度；V为两颗粒间液体体积可由下式求得 n第101页，共207页。 切向接触力学模型为 当两颗粒间间距n0时，颗粒间既有液桥力，还有弹性力，此时总作用力是两者的叠加。 第102页，共207页。3.6 表面粘附力学模型 当颗粒尺寸比较小时，

38、颗粒间可能产生粘聚力，此时可用粘聚力模型计算颗粒的法向作用力为（吸力）式中为粘连表面能； ； ； 。第103页，共207页。3.7 气固耦合力学模型 在自然界和工农业生产中，气固两相的耦合作用非常普遍。目前的分析方法可分为两种，即基于欧拉方法的双流模型和基于拉格朗日方法的颗粒轨道模型。气吹排种器气吸排种器清选机气力输送第104页，共207页。 双流模型即计算流体力学方法（CFD），采用数值方法求解方程。 双流模型把气体和颗粒均假设为流体，分别求解气固两相的质量和动量守恒方程如下（气相方程）。式中x、y、z为三维坐标；t为时间；u、v、w分别为 在x、y、z三维坐标轴的投影；p为静压；g为重力加

39、速度； 、 、 ，且 、 和 分别为颗粒相给气体作用力 在x、y、z三维坐标轴的投影。第105页，共207页。管道流动分析（稀相颗粒）车辆行驶空气阻力分析风机工作过程分析CFD方法及其软件的应用潜艇阻力和浮力分析第106页，共207页。 颗粒轨道模型又有两种方法，一种是不考虑颗粒间的接触碰撞力，适用于稀疏颗粒流。 另一种是考虑颗粒间的接触作用力，即DEM-CFD耦合方法，日本学者Tsuji于1993年首先提出。 当气体粘性为常数、气体不可压缩且将颗粒对气体的作用力作为源项时，三维气固耦合的气相控制方程为 第107页，共207页。 两相间的耦合作用，当固相浓度较高时采用双向耦合，当固相浓度较低时

40、，采用单相耦合即只考虑气体对颗粒的作用力。颗粒运动方程（球颗粒） 第108页，共207页。 颗粒与气体之间的相互作用力比较复杂，包括曳力、浮力、压力梯度力、虚假质量力、Magnus力、Basset力、Saffman力等。但对于密相颗粒、无粘干颗粒且颗粒的密度远大于气体密度时，可以只考虑曳力，单颗粒上的受气体作用力为（D i Felice公式）式中 为当前时步单颗粒i所受的气体作用力；d为颗粒i的直径； 为气体密度； 为颗粒i所在网格的气体速度，可通过求解气固耦合的气相控制方程求出； 为当前时步颗粒i运动速度，可由求解气固耦合的DEM计算求出； 为空隙率即网格内气体所占体积的百分比； 为单颗粒曳

41、力系数（颗粒阻力系数）。第109页，共207页。Ergun公式:W en 和Yu公式:第110页，共207页。DEM-CFD耦合方法的几个问题流体相的流态层流 求解N-S方程湍流直接数值模拟湍流方程，如时均方程、大涡模拟等直接求解N-S方程，网格尺寸和时步需取得非常小对于密相颗粒流，也采用层流方法，求解N-S方程第111页，共207页。 离散计算分析区域划分网格块 网 格适应性强非 结 构 网 格非规整、复杂边界适应性强贴体网格结构网格 规 整适应性差均匀网格同位网格交错网格第112页，共207页。方程离散方法有限差分方法 只是用于结构网格、均匀网格有限体积方法 适用于结构网格、非结构网格离散

42、格式对流相离散格式：中心差分、一阶迎风、二阶迎风和Quick格式非稳态相离散格式：显示积分方案、全隐式积分方案、半隐式扩散项相离散格式：中心差分第113页，共207页。给定参数求解气相控制方程方程确定初始条件和边界条件SIMPLEGauss-SeidelTDMA排种器CAD模型排种器分析模型排种器分析模型第114页，共207页。非球颗粒孔隙率求解气体给颗粒的作用力非球颗粒阻力系数非球颗粒垂直于流体速度的面积第115页，共207页。 我们采用结构网格和有限体积法，对流相中心差分、一阶迎风格式、二阶迎风和Quick格式等，Gauss-Seidel和TDMA方法求解离散方程，得到气场压力和速度，在此

43、基础上实现了DEM-CFD的耦合。第116页，共207页。FLUENT软件 自主研发软件（QUICK格式） 自主研发软件（二阶迎风格式） 直管道流自主研发软件（一阶迎风格式） 第117页，共207页。FLUENT软件 自主研发软件（QUICK格式） 自主研发软件（二阶迎风格式） 阶梯管道流自主研发软件（一阶迎风格式） 第118页，共207页。 自主研发软件（QUICK格式） 自主研发软件（二阶迎风格式） 阶梯管道流FLUENT软件 自主研发软件（一阶迎风格式） 第119页，共207页。第120页，共207页。气力输送仿真分析气吹排种器工作过程仿真分析第121页，共207页。3.8 连接力学模型

44、离散元法颗粒作用力计算的力学模型接触作用力模型（斥力）连接作用力模型（吸力）FnFnFn第122页，共207页。 接触和连接力学模型的耦合可用于粉碎、切断、脱粒等过程分析。玉米籽粒接触脱粒机与玉米籽粒接触OO玉米籽粒与玉米芯接触 分析脱粒时，玉米籽粒与玉米籽粒的接触及作用力，玉米籽粒与脱粒机的接触及作用力，可以采用接触力学模型计算；玉米籽粒与玉米芯的接触及接触作用力，可以采用连接力学模型计算。第123页，共207页。OZYXY方向压缩力- 时力位移曲线Y方向压缩力- 是刚度系数、断裂力、断裂位移第124页，共207页。OX方向压缩力 时刚度系数、断裂力、断裂位移Z方向压缩力 时刚度系数、断裂力

45、、断裂位移ZYX第125页，共207页。OZYX-FZFY=Fmax拉-FY=Fmax压FZFX-FX第126页，共207页。OZYX0246815304560夹角?/试验值计算值先育335号玉米对比 先锋8号玉米对比 F 0246815304560夹角?/试验值计算值第127页，共207页。(e) 小麦穗模型(d) 玉米穗模型(f) 小麦植株模型 我们采用连接力学模型建立并分析了：玉米脱粒过程三维仿真分析麦穗脱粒过程三维仿真分析第128页，共207页。库仑莫尔准则 当Fsk(t)jFnk(t)时，Fs(t) 由上式计算。 当Fsk(t)jFnk(t)时，Fs(t)=dFnk(t)。 此时Fs

46、(t)符号取其修正前的切向力符号。j为接触两体间的静摩擦系数，d为接触两体间的动摩擦系数。R1R2ij第129页，共207页。 (1)对农业生产中的散粒物料进行实际分析时，选择哪个力学模型更好，或如何改进已有的力学模型，还需要深入研究。 (2)农业土壤壤土属于非饱和的多相体系，本构关系复杂，含水率、质地、有机质含量、粒度和坚实度等对其性质都有影响。壤土颗粒间的相互作用力学模型，还需要深入研究。粮食干燥分析第130页，共207页。 颗粒材料和机械部件的物理机械性质不同，力学模型中参数，如刚度系数、弹性模量、泊松比、碰撞恢复系数、摩擦系数等也不同。 Klinker等通过测取颗粒作用力和变形的关系，

47、并用最小二乘法拟合成直线，从而得到kn的值。 Negi等通过振动试验，测取了大豆的自振频率及在此频率下的振动响应，并通过计算得到大豆的kn和cn。 对于半锁弹簧模型中的刚度系数k1，则应取力和变形曲线接近屈服点处的斜率。 3.9 力学模型中参数的确定 准静态压缩第131页，共207页。 Deshpande等通过试验测定了不同含水率(8.725%)时，大豆的尺寸分布、形状、球形率、容重、表面积、千粒重、密度等的变化。 LoCurto等研究了不同含水率的大豆种子，以不同速度碰撞塑料和金属板时，碰撞恢复系数e的变化规律。 Yang等还研究了不同散粒物料如大豆、小麦，当其形状、尺寸、含水率、碰撞速度、

48、碰撞角度不同时的碰撞特性。 Zhang等还提出一种由简单力学试验，求解大豆的弹性模量、泊松比和屈服应力的方法。 第132页，共207页。 Gupta等对葵花籽的物理机械性质，包括尺寸、球形率、密度、悬浮速度、摩擦系数等进行了研究。 Moya等对包括小麦、葵花籽在内的多种散粒农业物料的弹性模量、泊松比、内摩擦角、比重等物理机械性质进行了研究。 Molenda等还对小麦在不同条件下的摩擦系数进行了研究。第133页，共207页。 目前确定颗粒力学参数的方法主要有三种： 把颗粒当成一个整体，通过弹塑性理论分析和单颗粒试验得到颗粒的力学参数； 通过接触力学分析、简单试验和试凑方法得到颗粒的力学参数； 采

49、用宏观力学试验，如三轴试验、双轴试验和直剪试验、坚实度试验等，得到颗粒群体的宏观力学参数，然后建立宏观参数与微观颗粒参数的关系，由此得到颗粒的力学参数（反求，多解）。 准静态压缩vvv三轴试验直剪试验坚实度试验第134页，共207页。 (1)颗粒的物理机械性质影响因素多、分散较大，研究在离散元法分析时参数的选取方法，以使离散元法的分析结果与实际情况较吻合，还需要深入研究。 (2)土壤颗粒的力学参数很难直接测试，如何选取参数，还需要深入研究。 (3)玉米脱粒时，如何选取参数，还需要深入研究。 第135页，共207页。4. 邻居搜索与接触判断 在计算某颗粒与其它颗粒的接触作用力前，首先应判定与该颗

50、粒接触的颗粒或边界，该过程称接触检查。 为了减小接触检查的计算量，一般又把它分成两步来完成，第一步称邻居搜索，即采用某种方法确定与该颗粒较为接近的颗粒或边界(称邻居元) ，第二步称接触判断，即采用某种方法判定该颗粒与邻居颗粒或边界是否真正接触。 在离散元法分析中，常采用三种邻居搜索方法，即邻居列表法(Neighbor List)、网格法(Lattice or Boxing)和边界盒法(Bounding Box)。 第136页，共207页。 邻居列表法是Verlet 提出的一种方法，其基本思想以某一颗粒为中心，以一定长r(通常r3.1Rmax，Rmax为系统中最大颗粒的半径)为半径画一个圆或球，

51、则在该圆或球内的颗粒均作为该颗粒的邻居元，检查该颗粒与其它颗粒或边界是否接触，只需检查该颗粒是否与邻居元接触。当N为系统内的全部颗粒数量时，该方法精确确定系统内全部颗粒接触的计算复杂度为O(N2)。但由于离散元法分析时，时步通常取得较小，使颗粒在一个时步内位移较小，因而一旦颗粒的邻居元链表被确定，则可用于几个时步的接触判断，因而该方法的平均计算复杂度为O(N)，比直接检查系统内每一个颗粒与系统其它颗粒接触的计算复杂度O(N2)要小。当系统颗粒数N较大时，该方法的计算复杂度接近O(N2)；当增大r时，虽然可以在一定程度上降低该方法的计算复杂度，但过大的r也将增加计算时间。 第137页，共207页

52、。 网格法是Alder和Wainwright 提出的一种方法，其基本思想是以一定尺寸L为边长，把分析区域划分成规则的网格，则每一个颗粒均在某一个网格内。对于某一个颗粒来说，其所在的网格及其邻接的8个(二维)或26个(三维)网格内的颗粒或边界均作为该颗粒的邻居元，检查该颗粒是否与其它颗粒或边界接触，只需检查该颗粒与邻居元是否接触即可。当系统颗粒直径相同，而且L选取适当时，该方法的计算复杂度为O(N)：当系统颗粒尺寸差异较大或非球形颗粒长短径比较大时，计算复杂度为O(N2)。 圆弧段边界直线段边界圆形颗粒邻居元范围第138页，共207页。 边界盒法是Cohen和Baraff同时提出的一种方法，其基

53、本思想是把平行于坐标平面并与颗粒外接的盒(称边界盒)分别垂直投影到坐标轴x，y，z上，如果两颗粒相互接触，那么其边界盒在三个坐标轴上的投影必定重合。该方法当系统颗粒数N较少且颗粒密度较小时，计算复杂度为O(N)；当系统颗粒密度较大时，计算复杂度也为O(N2) 。一些研究人员还对上述方法进行了改进，还有把两种方法结合应用的报道（包围盒）。 第139页，共207页。 对于颗粒间的接触判断方法，一般来说与颗粒的几何模型有关。对于圆形或球形几何模型颗粒，只须简单判断两颗粒的中心距与它们半径之和的差，如果中心距小于半径之和，则他们相互接触，且法向叠合量为半径和与中心距之差，否则未接触。 第140页，共2

54、07页。 对于椭圆形几何模型颗粒，接触判断常用三种方法。一种是求两椭圆方程交点的方法，并把两交点中点作为接触点(见图)，但该方法有时很难求出精确解，当两椭圆叠合量非常小或两椭圆半轴共线时还可能无解，而且该方法不能直接推广到三维情况。 第二种称几何势方法，该方法是通过求解两接触椭圆分别进入对方最远点的两点中点作为交点。 第三种是公法线法，Lin等还证明一般情下几何势方法比公法线法精度更高、求解速度更快。在上述三种方法中，后两种方法也用于椭球的接触判断。 求两椭圆交点方法 求椭圆交点的几何势方法 公法线方法简图第141页，共207页。 (1)大长径比颗粒的邻居搜索方法。 (2)边界图元如何放入网格

55、内以进行邻居搜索，还需要深入研究。 (3)当采用超二次方程建立颗粒模型时，如何高效、高精度地进行颗粒与颗粒之间、颗粒与边界之间的接触判定还需要深入研究。 第142页，共207页。 一般的求解过程为：接触检查，以球颗粒为例，分别检查系统内每个颗粒与其它颗粒或边界的接触情况，求出接触点和接触叠合量n 。 5. 离散元法的求解过程和软件开发 圆弧段边界直线段边界圆形颗粒邻居元范围邻居搜索接触判定分离筛筛面求叠合量n边界n第143页，共207页。 颗粒作用力的计算，即应用接触力学模型分别计算每个颗粒上作用的合力和合力矩。第144页，共207页。 求解颗粒运动速度和位移，即根据牛顿第二定律和每个颗粒上的

56、作用的合力和合力矩，分别求解每个颗粒的平动速度和位移。 第145页，共207页。进入下一个时步，接触检查、求作用力和运动等。三维颗粒转动求解欧拉动力学方程第146页，共207页。 在离散元法中，求解颗粒速度和位移的方法可分为显式解法和隐式解法两种。显式解法用于动力学问题的求解和动态松弛法的静力学求解，隐式解法用于静态松弛法的静力问题求解。 在采用显式方法时，求解颗粒速度和位移常用中心差分法，也可采用其它数值积分方法，如Runge-Kutta法、Predictor-Corrector法和Newmark法等。第147页，共207页。 在采用显式求解方法时，由于假定在每一时步内颗粒的作用力、加速度为

57、常值，因而时步的确定对于计算的稳定性和精度非常重要。大部分研究人员均采用Cundall给出的公式计算临界时步为 式中mmin为系统中最小颗粒的质量；kmax为最大的接触刚度系数。 由上述求解过程可知，离散元法的最大缺点是计算量大、计算时间长，特别是分析数量较大的颗粒系统时，此时采用并行计算是很有必要的。 第148页，共207页。 离散元法发展到今天，已有大量应用软件出现。如ITASCA咨询公司开发的著名UDEC和3DEC及PFC-2D和PFC-3D软件，DEM Solutions公司的EDEM，Thornton研究组开发的著名GRANULE软件，及其它研究组开发的软件等，某些软件还可以通过英特

58、网下载可执行代码。国内也出现了GDEM、TRUDEC、SUPER-DEM等软件。第149页，共207页。 (1)软件集成方法，数据库、数据文件、命令流。 (2)离散元法的分布式并行算法、多核并行算法和GPU算法，OpenMP、MPI、CUDA、OpenCL。 (3)基于OpenGL的机械部件实体模型运动过程动态显示，数据文件格式STEP、IGES、STL、DXF等。 第150页，共207页。三、数字化设计软件开发及应用 我们从2003年开始离散元法方面的研究，采用离散元法（DEM）、计算流体力学方法（CFD）和多刚体动力学方法（MBD）及其耦合，研制出一种新型CAE软件，并实现了与CAD软件（

59、PRO/E、UG）集成，从而开发出集设计与性能分析为一体的数字化设计软件AgriDEM。 第151页，共207页。图83 一种相关机械部件数字化设计新方法和集成分析设计软件体系结构第152页，共207页。1. 软件的结构及功能 该软件由CAD设计子系统、工作部件离散元法分析模型建模子系统、离散元法计算子系统和性能分析仿真子系统四个部分组成。 CAD设计子系统是由CAD软件组成，设计者可以选取满足要求的设计，还可以对其进行修改，也可以重新进行设计，并将设计结果保存。该子系统还为建立工作部件的离散元法分析模型提供数据。 第153页，共207页。 离散元分析法模型建模子系统，是在CAD环境下进行的二

60、次开发，主要功能是通过人机交互读取机械部件的CAD模型（CAD设计图）中与颗粒接触的零部件表面（图元），并设定图元运动属性和材料特性等，建立“几何模型+运动模型+材料特性”的离散元法分析模型，并将模型存储到数据库中，为离散元法计算子系统和性能分析仿真子系统提供数据环境。 yxz第154页，共207页。 离散元计算子系统负责分析计算颗粒之间及颗粒与机械部件之间的接触情况、相互作用力和运动情况。用户首先在数据库中选取机械部件的分析模型，再根据颗粒的特性选择不同的力学模型和相关参数，然后进行离散元法计算，由此求出每个颗粒的运动速度和位移，并将这些信息以文件的形式记录下来。第155页，共207页。 性