使用离散元法DEM进行转鼓内颗粒动力学研究

Investigation of particle dynamics in a rotary drum by means of experiments and numerical simulations using DEM

摘要

【研究对象】
无隔板转鼓，不同操作条件，不同形状颗粒。

【研究方法】

离散元法DEM（Discrete Element Mehod）：用来描述颗粒流模式。
聚块法（clumping method）：可以创建不同的颗粒形状。

【结论】

颗粒形状对颗粒运动和DEM颗粒流预测有重要影响。
转鼓前后壁距离越大，静置动角越小，活跃区域的颗粒速度越小。

1.介绍

【背景】
颗粒流在许多工业过程中起着重要作用。理解这种流运动对于设计和优化问题是必须的，旋转桶广泛用于：制粒、混合、发酵、干燥、包衣等工艺。

【描述颗粒流的方法】（啥意思？）

Euler-Euler方程：使用动力学颗粒流，模拟固态压力分布，两相被看作互相贯穿连续，其守恒方程在欧拉参考系上求解。
拉格朗日方程：使用离散元法，对所有粒子进行追踪，考虑粒子与粒子和粒子与边界的相互作用，对每个粒子的进行不平衡力的显式时间积分。

【相关研究】
由于计算成本较高，大多数DEM仿真使用球体，使用不规则球体的研究较少。

Shamsi和Mirghasemi使用块聚法，使用一组重叠球体，研究颗粒形状对颗粒行为的影响。

Suzzi等人和Ketterhagen也使用块聚法，创建不同的颗粒形状（椭圆形、双凸形、子弹形……）来评估他们在转鼓和涂层锅内的流动行为。

图片为DEM模型中使用重叠球体模拟米粒

这些研究表明，颗粒形状在颗粒流研究中很重要。

在一些使用DEM方法的研究中，没有考虑不同操作情况下模型的验证和对等比放大模型的验证。大多数数值研究局限于确定等比缩小模型下的最佳参数集，没有考虑操作情况和几何效应的影响

Just等人发现实验摩擦力的测量不能充分描述颗粒在转鼓中的运动，需要进行敏感性分析来调整或校准DEM模拟参数。

【研究目的】
通过实验和DEM数值模拟，评估颗粒形状和转鼓壁对转鼓内颗粒动力学的影响。

数值研究分为三步：

DEM模型标定。
在不同操作情况下对标定模型的验证。
对等比放大模型的验证。

2.材料和方法

2.1实验方法和转鼓装置

【转鼓装置】
3种转鼓几何构型，2种不同长度的转鼓装置。

使用D1进行数值模型标定和随后不同操作情况的验证。填充程度f和转速w使用中心组合设计技术确定。将实验测得的静置动角θ与相同条件下的模拟值进行比较。

内直径为0.39m的转鼓用来验证等比放大的校正模型，这种转鼓有0.02m（D2）和0.45m（D3）两种长度，来观察转鼓长度对颗粒动力学的影响。

【颗粒速度】
测量床面中弦（参考线）处的颗粒速度分布，来和模拟数据对比。

颗粒速度分为vpx和vpy，通过同一颗粒在不同帧的位置相减来计算，参考线之前的位置（x1，y1）和之后的位置（x2，y2），然后除以两帧对应的时间间隔（t2-t1）。

【其他参数】
转鼓由透明亚克力制作，无隔板，分批操作。用砂纸覆盖转鼓内壁防止滑移。
玻璃珠直径3.68mm，密度2460kg/m3；米粒的当量直径3.44mm，密度1465kg/m3。
玻璃珠和米粒的孔隙度分别为0.38和0.44。

2.2数值方法

2.2.1拉格朗日模型描述

【DEM方法】
使用PFC3D软件进行模拟，它是使用基于离散元法DEM的代码来追踪每个颗粒的轨迹，包括它们之间以及和刚性边界的碰撞。

DEM计算方法在粒子-粒子和粒子-鼓壁接触处使用力位移定律，来更新不平衡接触力；在颗粒上使用牛顿第二定律，来确定颗粒的运动，这种运动考虑了作用在它们身上的body force（不知道怎么翻译）。

颗粒的平动和转动可以用牛顿第二定律和运动学关系来计算。
力位移定律，一般采用线性粘弹性弹簧阻尼器模型（？？？），如图所示。（看不懂这张图）

这个模型考虑了一种震动力，它把接触处的刚度系数k和颗粒间的相对位移δ和阻尼力联系起来，而阻尼力与阻尼系数λ以及颗粒的相对速度uij成正比。

接触合力Fcij包括法向Fnij和切向Ftij分量，表达式如下，n和t分别代表法向和切向。

滑动摩擦用库伦准则表示，其考虑了摩擦系数uf。
阻尼力和滑动摩擦可以用来消耗动能。

许多离散元公式使用阻尼比β而不是阻尼系数k。

临界阻尼常数λc用刚度系数k和有效质量meff描述
meff = (mimj)/(mi + mj):

根据阻尼比β的不同，系统有不同的行为方式。

临界状态，β=1，收敛到0，不震荡
过阻尼，β＞1，指数衰减，不震荡
欠阻尼，β＜1，指数衰减，震荡

阻尼比β和恢复系数e的表达式如下

DEM参数的设置需要考虑系统特性，比如颗粒特性和系统运算消耗。
较大的刚度系数k会导致较少的碰撞次数，时间步长就越小。

Wachs和Hohner等人认为，刚度系数k主要用来控制颗粒间的最大重叠，只有恢复系数e（用阻尼比β和摩擦系数uf表示）对粒子动力学有重要影响。

Cleary认为，0.1%—0.5%的颗粒半径的平均重叠是必要的，它确保颗粒流的行为不受刚度系数k的影响。在三维情况下，典型的刚度系数量级在104-106N/m。（取决于颗粒大小和最大力）

使用PFC3DO软件实现弹簧缓冲器模型仿真时，不需要明确滚动摩擦系数，只需要指定摩擦系数，阻尼比和刚度系数。

通过中心复合设计选择摩擦系数和阻尼比，并用静置动角表示，再对静置动角进行数值分析

2.2.2仿真条件

前提条件：认为法向与切向的刚性系数都为k，法向与切向的阻尼系数都为λ。

【实验内容】

模型标定：由于摩擦系数uf和阻尼比β的变化导致静置动角的变化，对其进行灵敏度分析。
不同操作情况下的模型验证：使用前一步确定的合适的参数集，验证DEM模型。
等比放大模型的验证：验证其在等比放大模型中的准确性和可靠性。

【实验参数】
刚度系数k=2105N/m，时间步长110-6—1*10-7秒，忽略额外的静电力。

3.结果和讨论

3.1颗粒形状对颗粒流预测的影响

模拟实验是在米粒不断滚动的情况下进行的。

情况1、2、5的玻璃珠床面发生变形，无法测量静置动角

根据表2的情况1和5，床体发生了滑坡运动（slumping），即滚动状态（rolling regime）这种情况有两种不同的静置动角，不能反映真实的实验现象。

滑坡状态（slumping regime）：颗粒被转鼓作为刚体提起，到达高静置动角（upper angle of repose），然后发生雪崩（avalanche），直到到达低静置动角（lower angle of repose）

表2中的情况1的转鼓横平面，（a）米粒的低静置动角（b）米粒的高静置动角（c）玻璃珠床面变形

上图（c），从鼓壁到临近颗粒层的动量转移很低效，因为摩擦系数和阻尼比很低。这种情况比较容易发生床面变形，因为离鼓壁较远的颗粒没有“感觉”到鼓壁的运动，保持静止。

米粒的仿真实验使用较低的DEM参数，床面较平坦，这是因为发生了一个堆积机制，不规则的颗粒相互锁结，在上游形成一种稳定的结构，使颗粒共同运动。

因此颗粒形状对颗粒流动有较大影响。

对比表2的5和6，7和8。

发现静置动角受摩擦系数影响较大，受阻尼比影响较小。

增大摩擦系数或阻尼比→使粒子的平动和转动产生更大阻力→消耗更多动能→静置动角更高

由于不规则颗粒有更高的接触面积

所以不规则颗粒比球形颗粒静置动角更高

3.2转速和填充度对模型的影响

使用转筒D1进行试实验。

米粒的不同情况下的静置动角

观察上表可知，预测结果和实验值吻合较好。

对于不规则颗粒，未发现填充度和转速对标定参数的影响。

对于不规则颗粒，填充程度增大或旋转速度减小→静置动角减小

玻璃珠的不同情况下的静置动角

对比结果表明

对于规则颗粒，填充度和转速对模型预测有显著影响，因此在DEM模型中不能一般化。

填充度较低时，实验和模拟的差值较大

3.3等比放大模型的验证

这一节主要评估DEM模型的鲁棒性，之前的实验使用的是小转鼓D1，现在使用大转鼓D2和D3（两转鼓长度不同，D2：0.02m；D3：0.45m）。
本节使用米粒进行，因为其校准参数不受上一节的操作条件影响（填充度和转速），并使用之前确定的最佳参数集进行模拟。
观察颗粒的顺向效应（consequent effect）。

转鼓D2和D3的实验和模拟的静置动角对比：

填充度10%，转速3.6rpm

填充度10%，转速6.2rpm

填充度15%，转速3.6rpm

可以看出

（无论转速和填充度）转鼓长度越短→静置动角越大

无论长度，对于一定的填充度，转速越大→静置动角越大

转速不变时，填充度越大→静置动角略减小

这与在小转鼓D1中观察到的现象一致。

使用大转鼓D2进行试验，与预测值有4.0%的误差；D2与D3使用分别的仿真参数，平均误差为20.0%。

因此，为了减少计算时间，使用等比缩小的模型，会显著增加预测的错误。

不活跃区域的颗粒速度与转速可表示为线性关系。

vlin为线速度，w为转鼓转速，r为转鼓径向位置（可能不是半径）
通过这个线性关系，可以确定滚动状态（rolling regime）下的活跃区域层与不活跃区域层。

不同填充度与转速下，总体，x分量，y分量的颗粒速度分布：
填充度10%，转速3.6rpm

填充度10%，转速6.2rpm

填充度15%，转速3.6rpm

由图可以看出，滚动状态下存在两个不同的区域：床面附近较薄的活跃区域（这里的颗粒瀑布状滚落cascade down）和鼓壁附近较厚的不活跃区域（这里的颗粒被鼓壁带着一起向上运动）。

在所有的情况下，床表面的颗粒速度最大，颗粒速度随着径向位置的减小而减小，在速度剖面的拐点（或Mellmann指出的“涡点”）处减小为零，此后发生逆流。

在不活跃区域，颗粒速度与上面的线性速度计算公式重合，因此可以转鼓和颗粒之间的滑移可以忽略，该区域的粒子作为一个刚体一起运动。

对于同一长度的转鼓D2，不论操作条件如何，活跃区域的颗粒速度随着转速和填充度的增加而增加

然而，转鼓D2与转鼓D3，实验和模拟的在活跃区域的颗粒速度也出现了偏差。
最短的转鼓D2，**壁效应（wall effect）**更加明显。颗粒被转鼓壁带到顶点，获得了较高的势能，随后颗粒进入活跃层，向下滚动将势能转化为动能。

因此，转鼓越长，静置动角越小，活跃区域的颗粒速度越小
（其实不理解为什么转鼓长度也会对静置动角有影响，颗粒的运动不是基本在转鼓的一个切面吗，和长度有啥关系？）

活跃层的颗粒速度并不是抛物线，基本是线性的，在mid-chord（？）的位置。

4.结论

通过实验和DEM数值分析，评估了转鼓内颗粒流的特性。包括颗粒形状，转鼓的操作条件（填充度，转速），转鼓尺寸对颗粒速度分布和静置动角的影响。

静置动角受摩擦系数影响较大，受阻尼比影响较小，阻尼比的大小采用中心符合设计。
除了模型参数，颗粒形状对颗粒运动和DEM预测起重要作用。不规则颗粒，无论填充度和转速如何，校准模型都可以进行预测；球形颗粒，DEM模型参数需要标定特定的转速和填充度，不能一般化。
转鼓尺寸也有重要影响，转鼓越长，静置动角越小，活跃区域的颗粒速度越小。（不理解为何转鼓越长，颗粒的势能就越小，有啥关系？？）
使用DEM标定的等比缩小的模型预测的颗粒速度分布与等比放大的实验测得的不规则颗粒的速度分布有较好的一致性。这有助于缩小微尺度DEM模型与真实工业系统的差距。