离散元法及其应用 教学提纲课件.ppt

离散元法及其应用-2014,一、引言,在自然界和工农业生产领域，大量存在着颗粒材料，如农产品、肥料、土壤、药品、煤炭和岩石等。据估计世界上50 % 的产品和75 % 的原材料都是颗粒材料。 在农业生产领域，耕地、开沟、播种、施肥、镇压、脱粒、分离、清选、粉碎、干燥、输送、仓储、分级、加工和包装等过程中，始终存在着颗粒材料与农机部件的接触作用和颗粒材料的流动过程。 在众多工业生产领域，如制药、食品、化工、冶金、采矿、能源、岩土工程等领域，也大量存在着颗粒材料与机械部件的接触作用和颗粒材料的流动过程。,一个好的农机部件设计应使: 颗粒材料按照预期的方式运动，如播种时种子流动， 减少流动过程中不必要的损伤，如播种时种子损伤， 节省动力消耗，如开沟和耕翻土壤时牵引动力消耗，等等，此时必须考虑机械部件与颗粒材料的接触作用及颗粒群体动力学问题 。 机械部件的优化需考虑颗粒动力学问题。,播种,开沟,颗粒材料的性质介于固体与流体之间，又称第四种物质形态，有着复杂的力学特性：非均匀尺寸偏析，如巴西果、反巴西果和三明治效应；粮仓效应；成拱现象；漏斗现象；自组织临界，等等。 自组织临界是Bak等1987年解释非线性复杂系统无序行为时提出的，即大的相互作用系统包含着众多短程相互作用的组元，系统自然地从随机状态演化到一种有序的临界状态，在该状态时小事件引起的连锁反应能够对系统中任何数目的组元产生影响，从而可能导致大规模事件的发生。,非均匀尺寸偏析,粮仓效应,成拱,固体,颗粒,颗粒材料通常指直径大于1m的颗粒组成。依据固体颗粒的浓度，可将颗粒材料分为密相颗粒材料、松散颗粒材料与稀薄颗粒材料。颗粒材料又可分为干颗粒材料(不含液体)与湿颗粒材料(含液体)。 颗粒材料流动可分为：准静态流动，流动的初始阶段，当颗粒承受的载荷超过颗粒间静摩擦力时，颗粒间仍保持接触但开始流动；快流，流动完全发展阶段的快速剪切流动；慢流，处于准静态流动和快流的中间阶段。,密相,松散,稀相,二十世纪70年代后，许多物理学家、力学家和应用数学家开始对颗粒运动的物理机制发生兴趣，建立了两类颗粒动力学理论：基于连续介质力学的理论，如颗粒动理论、摩擦塑性模型和光滑粒子法等；基于离散介质力学的理论，如硬颗粒模型、软颗粒模型和Monte Carlo方法等。 连续介质力学理论是把物质或其特性，假设成无论在时间还是在空间位置上，均是连续的或可用连续函数表示。因此物质可以无限分割而不失去其固有特性，不考虑粒子的特性，是描述物质整体及其特性的一种方法。,1. 基于连续介质力学的理论,颗粒动理论(kinetic theory) 研究发现快速颗粒流中单个颗粒的运动，与气体中的分子热运动非常相似。因此，借鉴非均匀的稠密气体分子运动理论，Ogawa定义了颗粒温度，Jenkins将气体的动理论扩展到颗粒材料，在考虑颗粒碰撞及摩擦所造成的能量损失的基础上，修正了Boltzmann方程，得到宏观的颗粒相输运方程，并导出动理论模型，由此可求得固体体积分数分布、颗粒速度分布和浓度分布等。 适合于稀薄颗粒的快流分析。,17世纪中期法国工程师Coulomb提出了土的抗剪强度和土压力滑动理论，其后被推广为散体极限破坏的Mohr-Coulomb准则，在此基础上发展成为土力学。 摩擦塑性模型，即是将Mohr-Coulomb准则应用于颗粒材料，当颗粒间载荷超过颗粒间的摩擦结合力，颗粒间开始滑移即屈服，但颗粒仍保持接触并相互摩擦。 人们已建立多种颗粒材料屈服条件，其中有双剪切模型、塑性势模型和双滑移自由转动模型等。 摩擦塑性模型主要应用于准静态颗粒流。,摩擦塑性模型,在小变形的情况下，可采用有限元法分析准静态的颗粒运动。此时，颗粒中心为单元节点。由接触建立节点间的联系，通过作用在节点上的力建立平衡方程。 有限元法适合模拟颗粒接触的拓扑结构不发生变化的静态颗粒系统，在动态和大变形的情况下，大量接触的丢失或产生，导致拓扑结构发生很大变化，这将需要耗费大量的时间重新生成单元，其缺点是：网格重构；网格变形较大还将产生计算不收敛；缺少合适的分析模型，如接触和变形模型、分离破裂模型等。,有限元方法,光滑粒子法思想是(SPH)：通过带质量的粒子离散计算域，粒子即代表颗粒，通过引入表征节点及其影响域内物理量间的关系核函数，来构造局部光滑的连续场，通过求解描述连续场对时间变化规律的常微分方程，来实现数值模拟。 光滑粒子法，在求解爆炸冲击及大变形问题等方面有不少应用，土壤切削等。,无网格方法(meshfree method),连续介质理论的基本控制方程是连续方程、动量方程和能量守恒方程。由于颗粒介质并不满足连续性的假定，并且由于连续介质模型没有考虑颗粒物性参数、粒径形状、大小及其分布等对颗粒流的影响，因此用连续介质模型分析颗粒流一般误差较大。 目前进行相关农机部件设计时，大都依靠经验和试验方法，既费时费力又得不到理想的效果。据估计仅由颗粒材料输送所造成的相关设备利用损失就达40%，远未达到优化设计和节省能源的要求。,随着计算机技术的发展，基于离散介质力学的理论，愈来愈引起人们的重视。 离散介质力学方法的思想源于较早的分子动力学，适用于模拟颗粒群体的接触或碰撞过程，它的出现补充了连续力学方法的不足。,2. 基于离散介质力学的理论,1985年Campbell提出硬颗粒模型，其思想是当颗粒表面承受的应力较低时，颗粒不产生显著的塑性变形，碰撞只在瞬间发生，在碰撞过程中颗粒本身不变形，并且只考虑两个颗粒的同时碰撞，而不计三个以上颗粒的同时碰撞，采用动量守恒或能量守恒计算碰撞后颗粒的速度和位置，广泛的应用于快速、低浓度颗粒流的模拟。,硬颗粒模型(hard/rigid sphere model),软颗粒模型又称为离散元法。 1971年Cundall提出适于岩石力学的离散元法(discrete /distinct element method，DEM)，1979年Cundall又提出适于土力学的离散元法，并推出二维圆盘程序BALL和三维圆球程序TRUBAL，后发展成商业软件PFC-2D/3D，形成较系统的模型与方法。,软颗粒模型(soft sphere model),块体模型,颗粒模型,离散元法的基本思想是，把散粒群体简化成具有一定形状和质量颗粒的集合，赋予接触颗粒间及颗粒与接触边界(机械部件)间某种接触力学模型和模型中的参数，以考虑颗粒之间及颗粒与边界间的接触作用和散粒体与边界的不同物理机械性质。,3. 离散元法的基本方法,离散元法的基本假设：单元是刚性的，即单元的几何形状不会因单元间的挤压力作用而改变；由于计算时步间隔取得足够小，单元的速度和加速度在一个时步内为常量，并且单元在一个时步内只能以很小的位移与其相邻单元作用，其作用力也只能传递到其邻接单元，而不能传递得更远；单元间的连接是靠相互接触实现的，圆形单元的接触为点接触等。,以i 和j 颗粒接触为例，设其法向叠合量为 ，由此产生的法向接触作用力 可如下计算（局部坐标胡克定律）,式中 为接触的法向刚度系数。 由于切向接触作用力与运动和加载历史有关，因此切向力通常采用增量形式计算，t 时刻的切向力 为,式中 为上一时步接触的切向作用力； 为接触的切向刚度系数； 为接触点的切向相对位移； 为计算时步。,（静摩擦力）,（斥力）,求i 颗粒质心作用的合力和合力矩为（全局坐标系下）,求i 颗粒的新位置有二种方法，静态松弛法和动态松弛法。,静态松弛法根据颗粒不平衡力达到再平衡时的力与位移关系建立平衡方程组，通过求解方程组得到颗粒的新位置，是一种隐式解法且需求解刚度矩阵。 动态松弛法采用牛顿第二定律求解颗粒的新位置如下。,动态松弛法需加入人工阻尼，使方程解收敛。因此，阻尼系数和时步选择对求解精度和收敛速度影响较大。,动态松弛法牛顿第二定律（球颗粒）：,遍历所有颗粒，然后进入下一时步； 显式解法，适合于求解非线性问题。,动态松弛法中心差分（数值积分）：,由于其离散的特点，在分析高度复杂的系统时，无论是颗粒还是边界均不需作大的简化；当赋予接触颗粒间不同的接触模型时，还可以分析颗粒结块、颗粒群聚合体的破碎过程、多相流动甚至可以包括化学反应和传热问题。 正是由于诸多优点，使得离散元法已成为研究颗粒群体动力学问题的通用方法，并在岩土工程及装备、采矿工程及装备、化工过程及装备、制药工程及装备、食品工程及装备和农业工程及装备等研究领域得到较多应用。,4. 离散元法的应用,图1 边坡稳定性分析,图2 块体拱的稳定性分析,图3 地下洞室稳定性分析,图5 料仓落料过程试验与仿真对比,(a) 试验,(b) 仿真,图10 材料压实过程模拟,图11 地基的夯实过程分析,图16 不同尺寸的筛分过程分析,实际球磨机,图21 混料过程的试验和模拟分析比较,(b) 模拟,(a) 试验,图23 滚筒式混合机的工作过程模拟,图27 苹果运输时损伤过程模拟,图40 铲装过程模拟,图42 圆锥体插入土壤过程模拟,图41 土壤和轮胎的相互作用分析FEM-DEM,图47 铲装过程模拟,图51 轨道与地基的作用过程分析,图50 机械部件与土壤的作用过程分析,图52 滑坡过程分析,图53 脆性材料剪切破坏过程模拟,材料的破坏过程，实质是力学模型从连续介质模型到离散介质模型的转变过程。,图54 复合材料受压破坏过程模拟,图64 导弹冲击脆性材料过程模拟,(a) 加强前,(b) 加强后,图68 颗粒的破碎过程模拟,图69 大钢球冲击过程模拟,图70 颗粒的破裂过程模拟,图71 Evolution of the coating process颗粒的包衣、磨损、腐蚀过程模拟,图72 脆性材料的切削过程模拟分析,图77 颗粒驱动水轮过程模拟分析,图78 冲击过程模拟分析,1.颗粒的离散元法建模方法 在采用离散元法分析颗粒材料的流动过程时，把颗粒简化成何种几何模型，如何给定几何模型中的参数建立其分析模型，将直接影响到分析精度和计算效率。,二、研究现状及存在问题,颗粒模型,块体模型,已报道的颗粒建模方法 超二次方程 统一建模方法自然界中70%以上的颗粒, 组合颗粒模型,统一计算方法,颗粒建模,非球颗粒,球颗粒超二次方程建模,规则形状大豆，超二次方程建模,非规则形状玉米、岩石，组合颗粒建模,玉米籽粒,大豆籽粒,规则形状颗粒组合充填建模,非规则形状颗粒充填组合建模方法,非规则形状颗粒,(b) 三点成球方法,冗余数据点去除充填过程优化质心和转动惯量计算,尺寸减小消除叠合,已报道的三种农作物颗粒模型,(a) 大豆籽粒模型,(b) 玉米籽粒模型,(c) 玉米籽粒模型,我们建立的颗粒模型（组合颗粒模型）,(a) 大豆种子模型,(c) 小麦种子模型,(b) 玉米种子模型,采用人机交互充填建模,基于点云的充填方法基于三角形网格面的充填方法,自动充填,(d) 稻米籽粒模型,籽粒采用人机交互充填或自动充填建模,我们建立的颗粒模型,颗粒材料的样本生成方法(1)随机生成法；(2)三角划分法。加密方法：(1)重力的方法；(2)加压和尺寸扩张方法。,分析计算区域,P,尺寸增大消除间隙,加压消除间隙,入料口,(1)当采用超二次方程建立颗粒模型时，如何高效、高精度地进行颗粒与颗粒之间、颗粒与边界之间的接触判定，还需要深入研究。 (2)当采用球组合方法建立颗粒模型时，如何建立一种非规则形状颗粒的组合方法，还需要深入研究。 (3)由于农业生产中的散粒物料种类繁多、形状复杂，而且即使是同一种散粒物料，其颗粒形状和尺寸参数差别也较大。因而在实际分析设计时，把散粒物料颗粒简化成何种几模型，如何给定几何模型中参数，是一个需要深入研究的问题，壤土、玉米植株、水稻等。,已报道的的边界建模方法(1)函数建模 简单边界的建模如一个圆柱面 国外著名商品软件PFC、 EDEM(2) 颗粒排列方法 简单边界的建模如一个平面(3) Kremmer等提出“有限壁”方法 复杂边界建模方法 国外著名商品软件EDEM,2. 边界的离散元法建模方法,边界,边界,“有限壁”方法是一种近似方法：求解颗粒与边界的接触点（接触力作用点）和接触叠合量（用于计算接触作用力）是近似的；在小三角形平面的相互连接处一阶导数不连续，致使所求颗粒与边界的法向和切向接触作用力具有突变性；在采用赫芝模型求解接触作用力时，曲面边界取无穷大的曲率半径，这些均与实际的边界情况差别较大。为了提高建模精度，往往还需要较多数量的三角形平面，由此增加了离散元法的计算时间。,基于CAD模型的边界建模方法,边界建模,非规则曲面不能用初等解析数函表示，离散成规则曲面,规则曲面平面、球面、 柱面、锥面等，可用初等解析函 数表示，当曲面上存在缺失部分时，提出实边 界和虚边界的方法,曲面上划分网格常用两种方法，一种是映射法另一种是直接法。采用了直接法中的推进波前法AFT（Advancing Front Technique）。,运动边界建模,运动形式,复杂运动一种联合收割机中分离筛等,简单运动平动、转动和平转动组合,平动、转动,复杂运动,播种单体,开沟,播种,铲车,履带车,平动 新位置 原位置+Vt,转动 新位置 原位置+ r t,平转动边界运动建模,复杂运动边界运动建模,运动形式,多体动力学方法,图解方法,分析方法,笛卡尔方法一般机械系统,拉格朗日方法 航天器系统,多体动力学方法及其软件的特点：自动建模、自动求解,为笛卡尔坐标阵的约束方程,为各刚体的质量矩阵,为约束方程的雅可比矩阵,笛卡尔位形坐标阵,为各刚体上受到的外力矩阵,为与约束方程对应的拉格朗日乘子阵,我们采用多刚体运动学中的笛卡尔方法，分别建立各种铰链约束方程，采用牛顿拉斐逊方法求解系统的约束方程、速度与加速度方程，得到每个时步、每个刚体的平动位移和速度、角位移和角速度等，实现了多刚体动力学（MBD）与三维离散元法（DEM）的耦合。,离散元法发展到今天，大部分研究都集中在颗粒的分析模型和接触力学模型等方面，对于边界建模的讨论还较少。 建立一种高效、精确和适应性广的离散元法边界建模方法，已成为离散元法实际应用急需解决的关键问题之一。,(1)如何建立一种非规则曲面边界的离散元法分析模型还需要深入研究，离散方法离散成三角形平面或球面或其它曲面；精度控制基于几何自适应；与CAD软件集成；自动建模等。 (2)如何建立复杂运动规律边界的离散元法分析模型还需要深入研究，非平动转动边界。 (3)如何建立弹性边界的离散元法分析模型还需要深入研究，弹性边界或橡胶材料边界等。,弹簧,在离散元法分析中，颗粒与颗粒之间、颗粒与边界之间的相互作用力，一般分为法向和切向两个方向，并分别采用不同的力学模型来计算。由于通常要分析计算的颗粒数量较多，采用较完备的力学理论来建立模型并由此计算作用力，不仅计算复杂，而且计算量较大，因此，通常都采用简化模型，常用的接触力学模型有以下几种。,3. 相互作用的力学模型及参数确定,线性粘弹性模型，法向接触作用力为（斥力） 式中Fn为接触两体间的法向作用力；kn为接触的法向刚度系数，n为接触两体的法向叠合量；cn为法向粘性阻尼系数，vn为两体接触处的法向相对速度。,3.1 线性粘弹性接触力学模型,Fn(N),（无粘干颗粒）,由于切向力的大小与加载历史有关，因而通常切向接触力的计算都采用增量形式，切向作用力为（静摩擦力） 式中Fs(t)为当前时步接触两体间切向作用力；Fsk(t)为当前时步接触两体间切向弹性力， ； 为上一时步接触两体间切向弹性力，t为计算时间步长；ks为接触的切向刚度系数，vs为接触处的切向相对速度； 为当前时步接触两体间切向阻尼力， ；cs为切向粘性阻尼系数。s=x, y，线性粘弹性模型适合于无粘干颗粒。,线性粘弹性模型虽然得到广泛应用，但实际上颗粒相互接触时，法向作用力都是非线性的，非线性粘弹性模型可由赫兹弹性接触理论得到为（斥力） 式中 ， ， ， R1、R2分别为接触两体接触处的曲率半径，E1、E2分别为接触两体的弹性模量， 、 分别为接触两体的泊松比。,3.2 非线性粘弹性力学模型,Fn(N),（无粘干颗粒）,切向接触力学模型为（静摩擦力）,s=x, y，非线性粘弹性模型适合于无粘干颗粒。,人们深入研究发现，在上述粘弹性模型中，当接触的法向叠合量n=0时，法向接触力Fn0，很显然这与实际情况相违背；而且实际接触的两颗粒间，粘性阻尼力也是非线性的。 为了克服这些缺点，人们又提出了非线性粘弹性模型的一般表达式为（斥力） 式中r可在12间取值；s可在01间取值。Mishra通过实验确定，对于钢球类颗粒，r=1.6、s=0.8比较合适；而在Tsuji等研究中，取r= 1.5、s=0.25。,(1)粘性阻尼力并非真正代表一种能量耗散机制，而只是使计算尽快稳定的方法，因而cn的确定具有较大随意性； (2)虽然cn可由e求得，但由于e不仅与材料特性有关，还与接触表面形状、接触碰撞速度等有关，这给正确地确定cn也带来难度； (3)当法向作用力较大时，在接触点处产生塑性变形，塑性变形将产生能量耗散，而此时e又较小，根据粘弹性模型计算的法向力已有较大误差。,3.3 弹塑性接触力学模型,为此，Mishra、Walton、 Vu-Quoc等又分别提出了几种弹塑性法向接触力学模型，其中Walton提出的半锁弹簧(又称双线性)模型为（斥力） 式中k1和k2分别为加载和卸载时的法向刚度系数；和0分别为接触两体的法向叠合量和残余法向叠合量，k1和k2还满足关系式 。,加载时,卸载时,切向粘弹性模型也存在着正确确定cs较难等缺点。为此，一些研究人员又根据Mindlin的弹性摩擦接触理论，并由不同简化得到几种不同的力学模型，其一般表达式为 式中切向刚度系数ks是一个变化值，需在每一时步不断重新计算，而ks的计算方法不同，力学模型则不同。Walton等提出的切向力学模型，其ks的计算式为 式中当接触的法向作用力由半锁弹簧模型计算时， ；为摩擦系数；Fn(t)为t时刻的法向作用力；Fs*为Fs曲线的转折点，图中切向位移S=vst。,Fs增大时,Fs减小时,上述Walton的切向力模型，只适用于Fn为常值时切向接触力的计算，但在多数情况下Fn是变化的，为此Vu-Quoc等又提出另一种切向力计算模型，其ks计算式为 式中当接触的法向作用力由半锁弹簧模型计算时， ；当接触的法向作用力由赫兹弹性接触的非线性粘弹性模型计算时， 。,Fs增大且|Fs|Fs*时,Fs增大且|Fs|Fs*时,Fs减小且|Fs|Fs*时,Fs减小且|Fs|Fs*时,Vemuri等还提出一种直接计算总切向力的模型，以代替上述增量模型，其计算公式为（静摩擦力） 式中Fs0为Fs的初值；S为某时刻的切向位移，S0为S的初值。,（无粘干颗粒）,3.4 湿颗粒(土壤)接触力学模型,当两颗粒中心距D R1+R2时，法向接触力由线性粘弹性模型计算为（斥力） 当两颗粒中心距D在 R1+R2 D (1+Cad)( R1+R2)间时，法向接触力为（吸力）,对上述模型进行了改进，当n0时,当 Ca n 0时,（吸力）,（吸力）,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为两颗粒间液体体积可由下式求得,切向接触力学模型为 当两颗粒间间距n0时，颗粒间既有液桥力，还有弹性力，此时总作用力是两者的叠加。,3.6 表面粘附力学模型,当颗粒尺寸比较小时，颗粒间可能产生粘聚力，此时可用粘聚力模型计算颗粒的法向作用力为（吸力）式中为粘连表面能； ； ； 。,3.7 气固耦合力学模型,在自然界和工农业生产中，气固两相的耦合作用非常普遍。目前的分析方法可分为两种，即基于欧拉方法的双流模型和基于拉格朗日方法的颗粒轨道模型。,气吹排种器,气吸排种器,清选机,气力输送,双流模型即计算流体力学方法（CFD），采用数值方法求解方程。,双流模型把气体和颗粒均假设为流体，分别求解气固两相的质量和动量守恒方程如下（气相方程）。,式中x、y、z为三维坐标；t为时间；u、v、w分别为 在x、y、z三维坐标轴的投影；p为静压；g为重力加速度； 、 、 ，且 、 和 分别为颗粒相给气体作用力 在x、y、z三维坐标轴的投影。,管道流动分析（稀相颗粒）,车辆行驶空气阻力分析,风机工作过程分析,CFD方法及其软件的应用,潜艇阻力和浮力分析,颗粒轨道模型又有两种方法，一种是不考虑颗粒间的接触碰撞力，适用于稀疏颗粒流。 另一种是考虑颗粒间的接触作用力，即DEM-CFD耦合方法，日本学者Tsuji于1993年首先提出。 当气体粘性为常数、气体不可压缩且将颗粒对气体的作用力作为源项时，三维气固耦合的气相控制方程为,两相间的耦合作用，当固相浓度较高时采用双向耦合，当固相浓度较低时，采用单相耦合即只考虑气体对颗粒的作用力。,颗粒运动方程（球颗粒）,颗粒与气体之间的相互作用力比较复杂，包括曳力、浮力、压力梯度力、虚假质量力、Magnus力、Basset力、Saffman力等。但对于密相颗粒、无粘干颗粒且颗粒的密度远大于气体密度时，可以只考虑曳力，单颗粒上的受气体作用力为（D i Felice公式）式中 为当前时步单颗粒i所受的气体作用力；d为颗粒i的直径； 为气体密度； 为颗粒i所在网格的气体速度，可通过求解气固耦合的气相控制方程求出； 为当前时步颗粒i运动速度，可由求解气固耦合的DEM计算求出； 为空隙率即网格内气体所占体积的百分比； 为单颗粒曳力系数（颗粒阻力系数）。,Ergun公式:,W en 和Yu公式:,DEM-CFD耦合方法的几个问题,流体相的流态,层流 求解N-S方程,湍流,直接数值模拟,湍流方程，如时均方程、大涡模拟等,直接求解N-S方程，网格尺寸和时步需取得非常小,对于密相颗粒流，也采用层流方法，求解N-S方程, 离散计算分析区域划分网格,块 网 格适应性强,同位网格,交错网格,离散格式,对流相离散格式：中心差分、一阶迎风、二阶迎风和Quick格式,非稳态相离散格式：显示积分方案、全隐式积分方案、半隐式,扩散项相离散格式：中心差分,给定参数求解气相控制方程方程,确定初始条件和边界条件,SIMPLE,Gauss-Seidel,TDMA,排种器CAD模型,排种器分析模型,排种器分析模型,非球颗粒孔隙率,求解气体给颗粒的作用力,非球颗粒阻力系数,非球颗粒垂直于流体速度的面积,我们采用结构网格和有限体积法，对流相中心差分、一阶迎风格式、二阶迎风和Quick格式等，Gauss-Seidel和TDMA方法求解离散方程，得到气场压力和速度，在此基础上实现了DEM-CFD的耦合。,FLUENT软件,自主研发软件（QUICK格式）,自主研发软件（二阶迎风格式）,直管道流,自主研发软件（一阶迎风格式）,FLUENT软件,自主研发软件（QUICK格式）,自主研发软件（二阶迎风格式）,阶梯管道流,自主研发软件（一阶迎风格式）,自主研发软件（QUICK格式）,自主研发软件（二阶迎风格式）,阶梯管道流,FLUENT软件,自主研发软件（一阶迎风格式）,3.8 连接力学模型,离散元法颗粒作用力计算的力学模型,接触作用力模型（斥力）,连接作用力模型（吸力）,接触和连接力学模型的耦合可用于粉碎、切断、脱粒等过程分析。,玉米籽粒接触,分析脱粒时，玉米籽粒与玉米籽粒的接触及作用力，玉米籽粒与脱粒机的接触及作用力，可以采用接触力学模型计算；玉米籽粒与玉米芯的接触及接触作用力，可以采用连接力学模型计算。,Y方向压缩力- 时力位移曲线,Y方向压缩力- 是刚度系数、断裂力、断裂位移,X方向压缩力 时刚度系数、断裂力、断裂位移,Z方向压缩力 时刚度系数、断裂力、断裂位移,先育335号玉米对比,先锋8号玉米对比,我们采用连接力学模型建立并分析了：,玉米脱粒过程三维仿真分析,麦穗脱粒过程三维仿真分析,库仑莫尔准则,当Fsk(t)jFnk(t)时，Fs(t) 由上式计算。 当Fsk(t)jFnk(t)时，Fs(t)=dFnk(t)。 此时Fs(t)符号取其修正前的切向力符号。j为接触两体间的静摩擦系数，d为接触两体间的动摩擦系数。,(1)对农业生产中的散粒物料进行实际分析时，选择哪个力学模型更好，或如何改进已有的力学模型，还需要深入研究。 (2)农业土壤壤土属于非饱和的多相体系，本构关系复杂，含水率、质地、有机质含量、粒度和坚实度等对其性质都有影响。壤土颗粒间的相互作用力学模型，还需要深入研究。,粮食干燥分析,颗粒材料和机械部件的物理机械性质不同，力学模型中参数，如刚度系数、弹性模量、泊松比、碰撞恢复系数、摩擦系数等也不同。 Klinker等通过测取颗粒作用力和变形的关系，并用最小二乘法拟合成直线，从而得到kn的值。 Negi等通过振动试验，测取了大豆的自振频率及在此频率下的振动响应，并通过计算得到大豆的kn和cn。 对于半锁弹簧模型中的刚度系数k1，则应取力和变形曲线接近屈服点处的斜率。,3.9 力学模型中参数的确定,Deshpande等通过试验测定了不同含水率(8.725%)时，大豆的尺寸分布、形状、球形率、容重、表面积、千粒重、密度等的变化。 LoCurto等研究了不同含水率的大豆种子，以不同速度碰撞塑料和金属板时，碰撞恢复系数e的变化规律。 Yang等还研究了不同散粒物料如大豆、小麦，当其形状、尺寸、含水率、碰撞速度、碰撞角度不同时的碰撞特性。 Zhang等还提出一种由简单力学试验，求解大豆的弹性模量、泊松比和屈服应力的方法。,Gupta等对葵花籽的物理机械性质，包括尺寸、球形率、密度、悬浮速度、摩擦系数等进行了研究。 Moya等对包括小麦、葵花籽在内的多种散粒农业物料的弹性模量、泊松比、内摩擦角、比重等物理机械性质进行了研究。 Molenda等还对小麦在不同条件下的摩擦系数进行了研究。,目前确定颗粒力学参数的方法主要有三种： 把颗粒当成一个整体，通过弹塑性理论分析和单颗粒试验得到颗粒的力学参数； 通过接触力学分析、简单试验和试凑方法得到颗粒的力学参数； 采用宏观力学试验，如三轴试验、双轴试验和直剪试验、坚实度试验等，得到颗粒群体的宏观力学参数，然后建立宏观参数与微观颗粒参数的关系，由此得到颗粒的力学参数（反求，多解）。,三轴试验,直剪试验,坚实度试验,(1)颗粒的物理机械性质影响因素多、分散较大，研究在离散元法分析时参数的选取方法，以使离散元法的分析结果与实际情况较吻合，还需要深入研究。 (2)土壤颗粒的力学参数很难直接测试，如何选取参数，还需要深入研究。 (3)玉米脱粒时，如何选取参数，还需要深入研究。,4. 邻居搜索与接触判断,在计算某颗粒与其它颗粒的接触作用力前，首先应判定与该颗粒接触的颗粒或边界，该过程称接触检查。 为了减小接触检查的计算量，一般又把它分成两步来完成，第一步称邻居搜索，即采用某种方法确定与该颗粒较为接近的颗粒或边界(称邻居元) ，第二步称接触判断，即采用某种方法判定该颗粒与邻居颗粒或边界是否真正接触。 在离散元法分析中，常采用三种邻居搜索方法，即邻居列表法(Neighbor List)、网格法(Lattice or Boxing)和边界盒法(Bounding Box)。,邻居列表法是Verlet 提出的一种方法，其基本思想以某一颗粒为中心，以一定长r(通常r3.1Rmax，Rmax为系统中最大颗粒的半径)为半径画一个圆或球，则在该圆或球内的颗粒均作为该颗粒的邻居元，检查该颗粒与其它颗粒或边界是否接触，只需检查该颗粒是否与邻居元接触。当N为系统内的全部颗粒数量时，该方法精确确定系统内全部颗粒接触的计算复杂度为O(N2)。但由于离散元法分析时，时步通常取得较小，使颗粒在一个时步内位移较小，因而一旦颗粒的邻居元链表被确定，则可用于几个时步的接触判断，因而该方法的平均计算复杂度为O(N)，比直接检查系统内每一个颗粒与系统其它颗粒接触的计算复杂度O(N2)要小。当系统颗粒数N较大时，该方法的计算复杂度接近O(N2)；当增大r时，虽然可以在一定程度上降低该方法的计算复杂度，但过大的r也将增加计算时间。,网格法是Alder和Wainwright 提出的一种方法，其基本思想是以一定尺寸L为边长，把分析区域划分成规则的网格，则每一个颗粒均在某一个网格内。对于某一个颗粒来说，其所在的网格及其邻接的8个(二维)或26个(三维)网格内的颗粒或边界均作为该颗粒的邻居元，检查该颗粒是否与其它颗粒或边界接触，只需检查该颗粒与邻居元是否接触即可。当系统颗粒直径相同，而且L选取适当时，该方法的计算复杂度为O(N)：当系统颗粒尺寸差异较大或非球形颗粒长短径比较大时，计算复杂度为O(N2)。,边界盒法是Cohen和Baraff同时提出的一种方法，其基本思想是把平行于坐标平面并与颗粒外接的盒(称边界盒)分别垂直投影到坐标轴x，y，z上，如果两颗粒相互接触，那么其边界盒在三个坐标轴上的投影必定重合。该方法当系统颗粒数N较少且颗粒密度较小时，计算复杂度为O(N)；当系统颗粒密度较大时，计算复杂度也为O(N2) 。一些研究人员还对上述方法进行了改进，还有把两种方法结合应用的报道（包围盒）。,对于颗粒间的接触判断方法，一般来说与颗粒的几何模型有关。对于圆形或球形几何模型颗粒，只须简单判断两颗粒的中心距与它们半径之和的差，如果中心距小于半径之和，则他们相互接触，且法向叠合量为半径和与中心距之差，否则未接触。,对于椭圆形几何模型颗粒，接触判断常用三种方法。一种是求两椭圆方程交点的方法，并把两交点中点作为接触点(见图)，但该方法有时很难求出精确解，当两椭圆叠合量非常小或两椭圆半轴共线时还可能无解，而且该方法不能直接推广到三维情况。 第二种称几何势方法，该方法是通过求解两接触椭圆分别进入对方最远点的两点中点作为交点。 第三种是公法线法，Lin等还证明一般情下几何势方法比公法线法精度更高、求解速度更快。在上述三种方法中，后两种方法也用于椭球的接触判断。,求两椭圆交点方法,求椭圆交点的几何势方法,公法线方法简图,(1)大长径比颗粒的邻居搜索方法。 (2)边界图元如何放入网格内以进行邻居搜索，还需要深入研究。 (3)当采用超二次方程建立颗粒模型时，如何高效、高精度地进行颗粒与颗粒之间、颗粒与边界之间的接触判定还需要深入研究。,一般的求解过程为：接触检查，以球颗粒为例，分别检查系统内每个颗粒与其它颗粒或边界的接触情况，求出接触点和接触叠合量n 。,5. 离散元法的求解过程和软件开发,接触判定,颗粒作用力的计算，即应用接触力学模型分别计算每个颗粒上作用的合力和合力矩。,求解颗粒运动速度和位移，即根据牛顿第二定律和每个颗粒上的作用的合力和合力矩，分别求解每个颗粒的平动速度和位移。,进入下一个时步，接触检查、求作用力和运动等。,三维颗粒转动求解欧拉动力学方程,在离散元法中，求解颗粒速度和位移的方法可分为显式解法和隐式解法两种。显式解法用于动力学问题的求解和动态松弛法的静力学求解，隐式解法用于静态松弛法的静力问题求解。 在采用显式方法时，求解颗粒速度和位移常用中心差分法，也可采用其它数值积分方法，如Runge-Kutta法、Predictor-Corrector法和Newmark法等。,在采用显式求解方法时，由于假定在每一时步内颗粒的作用力、加速度为常值，因而时步的确定对于计算的稳定性和精度非常重要。大部分研究人员均采用Cundall给出的公式计算临界时步为 式中mmin为系统中最小颗粒的质量；kmax为最大的接触刚度系数。 由上述求解过程可知，离散元法的最大缺点是计算量大、计算时间长，特别是分析数量较大的颗粒系统时，此时采用并行计算是很有必要的。,离散元法发展到今天，已有大量应用软件出现。如ITASCA咨询公司开发的著名UDEC和3DEC及PFC-2D和PFC-3D软件，DEM Solutions公司的EDEM，Thornton研究组开发的著名GRANULE软件，及其它研究组开发的软件等，某些软件还可以通过英特网下载可执行代码。国内也出现了GDEM、TRUDEC、SUPER-DEM等软件。,(1)软件集成方法，数据库、数据文件、命令流。 (2)离散元法的分布式并行算法、多核并行算法和GPU算法，OpenMP、MPI、CUDA、OpenCL。 (3)基于OpenGL的机械部件实体模型运动过程动态显示，数据文件格式STEP、IGES、STL、DXF等。,三、数字化设计软件开发及应用,我们从2003年开始离散元法方面的研究，采用离散元法（DEM）、计算流体力学方法（CFD）和多刚体动力学方法（MBD）及其耦合，研制出一种新型CAE软件，并实现了与CAD软件（PRO/E、UG）集成，从而开发出集设计与性能分析为一体的数字化设计软件AgriDEM。,图83 一种相关机械部件数字化设计新方法和集成分析设计软件体系结构,1. 软件的结构及功能,该软件由CAD设计子系统、工作部件离散元法分析模型建模子系统、离散元法计算子系统和性能分析仿真子系统四个部分组成。 CAD设计子系统是由CAD软件组成，设计者可以选取满足要求的设计，还可以对其进行修改，也可以重新进行设计，并将设计结果保存。该子系统还为建立工作部件的离散元法分析模型提供数据。,离散元分析法模型建模子系统，是在CAD环境下进行的二次开发，主要功能是通过人机交互读取机械部件的CAD模型（CAD设计图）中与颗粒接触的零部件表面（图元），并设定图元运动属性和材料特性等，建立“几何模型+运动模型+材料特性”的离散元法分析模型，并将模型存储到数据库中，为离散元法计算子系统和性能分析仿真子系统提供数据环境。,离散元计算子系统负责分析计算颗粒之间及颗粒与机械部件之间的接触情况、相互作用力和运动情况。用户首先在数据库中选取机械部件的分析模型，再根据颗粒的特性选择不同的力学模型和相关参数，然后进行离散元法计算，由此求出每个颗粒的运动速度和位移，并将这些信息以文件的形式记录下来。,性能分析仿真子系统，是在离散元法计算的同时或完成后，将计算所得数据以图形仿真的形式重现出来，并提供颗粒运动速度、颗粒上作用力、颗粒作用力场、颗粒运动速度场显示、颗粒流量统计、机械部件工作阻力分析等功能，用以分析机械部件的工作性能。,该软件可由农机部件的CAD模型自动生成其离散元法分析模型，还可自动生成多种几何模型的颗粒材料，可采用多种接触力学模型进行接触力的计算，并采用人机交互方式输入计算参数，因而具有广泛的适用性。,为检验所开发软件的可用性，采用该软件对精密排种器进行了设计和分析.,2. 软件的应用研究,图84 组合内窝孔精密排种器 的CAD模型(主视图),图86 改变排种轮直径时组合内窝孔精密排种器的工作过程仿真,(a) 直径为158 mm时,(b) 直径为218 mm时,基于DEM的仿真,图88 改变散粒物料几何模型时组合内窝孔精密排种器的工作过程仿真,(a) 圆颗粒,(b) 椭圆颗粒,(c) 多边形颗粒,(b)椭圆形颗粒模型的测试图,(a)圆角矩形颗粒模型的测试图,图89 改变散粒物料几何模型时组合内窝孔精密排种器的工作过程仿真,图90 改变排种轮角速度时组合内窝孔精密排种器的工作过程仿真,(a) = rad/s时,(b) = -3/ 2 rad/s时,图92 改变CAD模型和接触力学模型后一种推土板的工作过程仿真,图93 改变CAD模型和接触力学模型后一种推土铲的工作过程仿真,图94 改变CAD模型和接触力学模型后肥箱的排肥过程过程仿真,(a) 一种肥箱,(b) 另一种肥箱,图95 改变CAD模型后一种排 肥器过程过程仿真,图96 改变CAD模型后一种开沟器过程过程仿真,图98 组合内窝孔精密排种 器内种子运动试验,图99 组合运动边界的分析计算,图100 由CAD模型得到的播种机整机工作过程的二维仿真分析,图107 颗粒流量分析,图101单粒种子跟踪显示,图102单颗粒受力变化曲线,图103单颗粒速度变化曲线,图106机械部件工作阻力变化曲线,图104颗粒作用力场显示,图105 颗粒运动速度场显示,(a) 颗粒跌落碰撞过程实验,(b) 颗粒跌落碰