基于离散元法的柔性片烟建模及仿真参数标定

金 浩，蒋明洋，徐 超，沈 凯，潘凡达，樊 虎，焦凯旋

（1.华东理工大学化学工程系，上海 200237；2.浙江中烟工业有限责任公司，杭州 310024）

中国作为全球烟草大国，烟草行业所产生的税收在国家财政收入中占据重要地位。在烟叶的打叶复烤过程中，松散润叶作为关键环节备受关注[1]。然而，当前观察烟叶在润叶筒内运动行为的直接方法存在困难，同时烟叶的润叶效果通常表现出不稳定性，这直接影响烟叶的柔韧性和口感[2]，从而制约了烟草行业生产过程的优化。此外，松散润叶过程中的烟叶呈现出片状柔性，而传统的离散元软件（如EDEM）常采用球体或多球体嵌合颗粒群作为片烟的颗粒模型[3]，存在颗粒模型设置复杂和仿真时间过长的问题。本研究首次在烟草仿真领域采用离散元（Rocky DEM）软件，创建壳体柔性片烟颗粒模型，并进行“虚拟试验”标定烟叶的离散元参数[4]。

近年来，离散元计算方法已广泛应用于散装物料的工程领域[5-6]。研究散装物料的堆积角可为颗粒运动行为的输送、混料、筛选等问题提供理论基础，因此国内外已经展开了许多相关研究。陶志影等[7]采用离散元法并通过正交试验设计对谷壳垫料颗粒进行参数标定，为翻抛刀具的设计和离散元参数设置作了理论指导；闫建伟等[8]将EDEM 与FLUENT 软件进行耦合仿真，对白萝卜排种器排种过程研究，通过Plackett-Burman 试验筛选出显著性颗粒物性参数，为白萝卜种子排种器设计和优化提供了理论依据；史瑞杰等[9]以胡麻茎秆为研究对象，在胡麻茎秆的收割仿真过程中研究了缺乏柔性模型和接触参数的问题，通过设计Plackett-Burman 试验和Central Composite 试验标定了离散元参数，准确表征了胡麻茎秆的运动特性；侯杰等[10]为标定叶鞘包裹的水稻茎秆的离散元参数，采用Plackett-Burman、Central-Composite 和Box-Behnken 试验方法获得了可靠性接触参数；朱新华等[11]以不同含水率畜禽粪肥为研究对象，研究了含水率对离散元仿真参数的影响，并通过Hertz Mindlin with JKR黏结模型直接预测畜禽粪肥的离散元参数。

本研究采用可灵活创建接近真实颗粒的Rocky DEM 2021R2软件作为数值模拟软件，并选择在松散润叶阶段进料口处具有不同尺寸的片烟作为研究对象。通过数字图像处理技术测得物理和仿真试验的片烟堆积角。通过Design Expert 软件设计Plackett-Burman、最陡爬坡和Box-Behnken 试验，筛选出影响片烟颗粒堆积角的显著性参数以及其拟合方程。通过堆积角的实际和仿真试验对比，以及对片烟堆积过程柔性的研究，验证了参数准确性。

1 材料与方法

1.1 材料

1.1.1 尺寸分布 试验中所用的材料均来自于“浙江中烟工业有限公司”松散润叶阶段进料口处片烟，此时片烟呈柔性片状。随机选择300片片烟进行尺寸筛分（图1），以最长边为测量标准进行归类，考虑到单片片烟存在折叠现象，每片片烟的平均厚度约为0.8~1 mm，不同片烟的平均尺寸和比例见表1。

表1 片烟大小分布Table 1 Tobacco strips size distribution

图1 真实片烟颗粒模型Figure 1 Real tobacco strips particle model

1.1.2 含水率和表面黏附力测定 采用CSY-L2片烟水分测定仪，先后5次测定柔性片烟物料，该仪器持续测量并即时显示样品丢失的水分含量（%）。根据5次测量结果统计，片烟含水率为12%~15%。在松散润叶设备内选取平衡相对湿度为75%~80%，平衡时间为24 h，平衡温度为30 ℃的稳定环境下片烟，采用TA.XT plus 质构仪（英国Stable Micro Systems 公司）测得烟叶表面黏附力为2.623~12.167 mN。

1.1.3 堆积密度测定 称取一定质量的片烟物料，并通过料斗自然落入到量筒中，通过质量除以量筒内的真实体积，重复5次试验，计算得到平均值为223 kg·m-3。

1.1.4 片烟堆积角的物理试验测定 堆积角是指大量物料堆积在水平面上形成锥体，且锥体料堆保持稳定状态的最大锥角，即物料堆积的边缘轮廓与水平面间的夹角。本研究采用固定漏斗法[12]进行片烟堆积角试验，通过图片处理软件ImageJ对片烟堆积角边缘轮廓进行图像读取，对图2分别删除背景、高斯模糊、二值化，提取轮廓边界处理。图3为堆积角边缘轮廓二值化，图4为Origin软件利用最小二乘法将提取出的片烟堆积轮廓进行线性拟合。重复试验10次并进行左右堆积轮廓统计，堆积角的平均值为35.83°。

图2 片烟堆积角边缘轮廓Figure 2 Contour of tobacco strips repose angle

图3 片烟堆积图二值化Figure 3 Image binarization of tobacco strips repose angle

图4 堆积角试验边缘拟合曲线Figure 4 Edge fitting curve of repose angle experiment

1.2 仿真原理及模型选取

Rocky DEM 离散元软件中包含3 类典型的接触力模型：法向力模型、切向力模型、黏附力模型。其中，选取“Hysteretic Linear spring”为法向接触模型，该模型计算颗粒能量耗散时间较短，能完美模拟颗粒流动[13]；选取“Linear Spring Coulomb Limit”为切向接触力模型，该模型更具普遍适用性，颗粒在摩擦开始前具有弹性，可在静摩擦值和动摩擦值下恢复库伦摩擦接触行为[14]。选取“Constant Adhesive Force”作为黏附力模型，该模型用于模拟不表现应力固结效应的黏附颗粒行为，如液桥力（该模型参考ESSS-Rocky 2021 R2 手册）。

（1）法向接触模型方程：

（2）切向接触模型方程：

式中：和分别为当前时间和先前时间的切向接触力；∆δt为时间步长内法向重叠变化量；μ为接触过程中滑动的摩擦系数。

（3）黏附力接触模型方程：

在两个不同质量的粒子接触的情况下，考虑较小的颗粒质量来计算重力，即：

式中：为法向接触黏附力；为法向接触重叠；m1和m2为接触颗粒的重量；g为颗粒的重力加速度;δadh为颗粒间或颗粒与润叶筒壁面间的黏附力距离；fadh为力分数。

1.3 片烟建模与形变分析

由表1可知，在选取片烟的样本中，中片片烟平均尺寸在6.5~12.5 cm间占比为75%，占大多数，因此该颗粒设置两个尺寸，分别为10.72 cm和8.04 cm；大片片烟设置尺寸为13.4 cm；小片片烟设置为5.3 cm。尺寸比例分布按照表1进行划分。

有关烟草类离散元仿真，着重以烟丝为研究对象，且颗粒模型大多为球形或球形黏结刚性颗粒，其与实际柔性片烟形态相差较远。为了更好地描述片烟的形状和柔性，本研究采用Rocky DEM 仿真软件，导入外部spaceclaim 创建的片状壳体颗粒模型文件stl.，设置壳体并对其进行网格划分，从而实现柔性颗粒的仿真计算（该模型创建参考ESSS-Rocky 2021 R2 手册）。模型由具有均匀厚度的二维三角网格组成，每个三角网格称为一个壳体颗粒元，每两个相邻壳体颗粒重叠面成为联结面（图5）。

图5 两个壳体元间的受力和力矩Figure 5 Forces and moments exerted by a joint between two shell elements

其中，为壳体颗粒元1具有的法向力和切向力，为壳体颗粒元1发生形变时的扭转和两种弯曲力矩，计算为：

式中：U、J、I1、I2分别为联结面的面积、二阶矩和壳体颗粒元1 形变分离时联结面的二阶矩，联结面尺寸为重合线长和壳体厚度组成的矩形；和分别为联结面在法向、切向、扭转和弯曲下所具有的弹性比值；E为壳体颗粒元的杨氏模量，一般被认为等于颗粒材料的杨氏模量；γ为联结面的特征长度，定义为重合线长度的倍数；v为颗粒材料的泊松比；θT、θB1和θB2分为壳体元形变扭转角度和两种弯曲角度。

每个壳体颗粒元为弹性体，形变仅发生在接触联结面。为了直观表示两个壳体元间的形变方式，以三角形表示壳体元，重合线为联结面，得到扭转形变和弯曲形变示意图(图6)，其中红色箭头表示壳体旋转轴。其中扭转形变发生在不同水平面，第二弯曲形变发生在同水平面（忽略壳体厚度）。因此，通过烟叶模型的每个壳体颗粒元间的相互作用实现了真实状态下烟叶的柔性形变特征。

图6 壳体形变示意图Figure 6 Schematic diagram of shell deformation

结合图1 颗粒形状，本研究将片烟设置为近似长方形，以长宽比为2∶1，建立具有粒径分布的柔性片烟模型。虽然颗粒网格数越多越好，考虑到计算机仿真时间，柔性片烟网格数划分为21（图7）。其对应粒径片烟的类型和比例如表2。

表2 仿真片烟粒径分布Table 2 Particle size distribution of simulated tobacco strips

图7 模拟片烟颗粒模型Figure 7 Simulatied tobacco strips particle model

1.4 标定参数范围确定

有关烟丝刚性颗粒的物料特性研究较多[15-16]，而关于柔性片烟的物料特性研究并不多。考虑到柔性片烟和烟丝本征参数接近，参考文献[14-17]，设定片烟的泊松比为0.25~0.4；杨氏模量为0.1~0.5 GPa；片烟间碰撞恢复系数0.2~0.45。此外，Rocky DEM 软件针对柔性片状颗粒无需设置滚动摩擦，仅考虑滑动摩擦（参考ESSSRocky 2021 R2 手册），片烟间的静摩擦系数为0.4~0.8[14,17]；滑动摩擦系数为0.2~0.7[14]（原则上滑动摩擦系数小于静摩擦系数）；片烟间黏附力系数结合式（5）计算为0.05~0.25；片烟与不锈钢间恢复系数为0.27；片烟与不锈钢间静摩擦系数为0.3~0.6；片烟与不锈钢间滑动摩擦系数为0.1~0.5。通过查阅不锈钢物性特征，泊松比为0.243；密度为7 920 kg·m-3；杨氏模量为460 GPa，重力加速度为9.81 m·s-2。软件中片烟产生方式为连续注入，片烟的质量生成速率为3.5 t·h-1，片烟总数量为1 398个，时间步长为0.05 s，模拟时间为10 s，计算机仿真时间为5 h。

2 结果与分析

2.1 Plackett Burman试验结果与分析

为确定影响片烟堆积角的显著性参数，本研究运用Design Export13软件进行Plackett-Burman试验。

Plackett-Burman 试验是基于响应参数与各因素参数间的关系，比较每个因素间两水平的差异并筛选出显著性影响因素。本试验以片烟堆积角为响应值，将表3 中X1~X8进行高低水平编码[20]，并增加3 个虚拟参数用于平衡模型误差，共计12组试验，试验设计与结果如表4。

表4 Plackett-Burman 试验设计和结果Table 4 Plackett-Burman test design and results

由Plackett-Burman 对堆积角进行显著性结果分析（表5），该模型p＜0.01,R2=0.993 8,表明模型很好地反映出影响因子显著性情况，说明该模型的拟合结果与实际情况相吻合，其中p值的大小表示试验参数对响应值堆积角的影响显著程度。分析可得：片烟-片烟间的静摩擦系数（X4）和片烟-片烟间的滑动摩擦系数（X5）影响极其显著（p＜0.01）；片烟-片烟间的黏附系数（X6）影响显著（0.01＜p＜0.05）。因此，只考虑这3 个显著性参数进行最陡爬坡试验，其余影响参数认为其显著性不明显。

表5 Plackett-Burman测试参数的显著性分析Table 5 Significance analysis of Plackett-Burman test parameters

2.2 最陡爬坡试验设计

根据Plackett-Burman 试验结果，对筛选出来的3个显著性参数片烟与片烟间的静摩擦系数（X4）、片烟与片烟间的滑动摩擦系数（X5）和片烟与片烟间的黏附系数（X6）进行最陡爬坡试验。选取仿真堆积角与物理试验堆积角间相对误差为评价标准，以快速逼近最优参数选取。取各参数范围中值（±0.05）为起始值，根据各组数据情况设置步长。按照最陡爬坡试验设计方案调整各参数，进行虚拟试验得到仿真堆积角，计算相对误差（表6）。由表6 可知，随着显著性参数增大的同时，相对误差呈先减小后增大的变化趋势，其中3 号试验参数组情况下的相对误差最小，为1.25%。因此可知最优参数取值区间在3号参数组附近。

表6 最陡爬坡试验设计及结果Table 6 Design and results of the steepest ascent search experiment

2.3 Box-Behnken试验结果与分析

Box-Behnken 试验设计属于响应曲面设计类型，将试验参数设计为位于试验空间边缘中点处的参数组合。根据表5 的Plackett-Burman 试验结果，不显著试验因素取范围中值分别为：片烟泊松比为0.32、片烟与片烟间恢复系数为0.32、片烟杨氏模量为0.25 GPa、片烟与钢板间静摩擦系数为0.5、片烟与钢板间滑动摩擦系数为0.3。根据“最陡爬坡试验”，以片烟与片烟间静摩擦系数（X4）、片烟与片烟滑间动摩擦系数（X5）、片烟与片烟间黏附力系数（X6）为试验因素并进行取值，堆积角为响应指标，用表7 进行三因素三水平试验设计，共计进行17组试验，试验设计与结果如表8。由表8可知，目标堆积角值35.83°也包含于试验结果内，这表明Box-Behnken试验高中低水平选取合适。

表7 Box-Behnken 试验参数Table 7 Box-Behnken test parameters

表8 Box-Behnken 试验设计及结果Table 8 Box-Behnken test design and results

对Box-Behnken试验结果建立响应参数堆积角θ与3个影响因素的三元回归模型，三元回归模型的表达式为：

该三元回归方程的p=0.001 7，决定系数R2=0.939 3，决定系数的校正值AdjR2=0.861 2，信噪比为11.317 7，说明该三元回归方程模型显著性较强；变异系数3.96%，失拟项p=0.071 1>0.05，方程拟合较好，可以对目标值进行合理预测。对该三元回归模型进行方差分析(表9)，结果表明：片烟与片烟间滑动摩擦系数（X5）对堆积角影响极显著（p＜0.01）；片烟与片烟间静摩擦系数（X4）系数和片烟与片烟间静摩擦系数的二次方（X24）对堆积角影响显著（0.01＜p＜0.05）；其余参数对堆积角影响不显著。

表9 Box-Behnken 试验三元回归模型方差分析Table 9 Variance analysis of Box-Behnken test ternary regression model

2.4 接触参数的标定结果与分析

运用Design export 软件进行模型优化，以柔性片烟堆积角实验值35.83°为目标值，在参数值范围内对式（11）三元二次回归模型寻优求解，以优化后的参数值进行片烟仿真堆积角试验，并通过与物理试验测定结果进行对比校正，得到堆积角轮廓较为相近的一组最优解为：片烟与片烟间静摩擦系数（X4）为0.67，片烟与片烟间滑动摩擦系数（X5）为0.31，片烟与片烟间黏附力系数（X6）为0.2，其余非显著性参数取平均值。采用上述参数作为Rocky DEM 软件的参数，进行三次仿真试验，得到片烟的堆积角分别为35.42°、35.97°、35.28°。物理试验的堆积角平均值为35.83°，仿真试验的堆积角平均值为35.55°，二者的相对误差为0.78%,从而验证了模拟试验的可靠性与真实性，对比图如图8。

图8 片烟堆积角物理试验与仿真试验对比Figure 8 Comparison between physical test and simulation test of repose angle of tobacco strips

3 讨论与结论

目前烟草的数值仿真以“柔性片烟”为研究对象的报道较少，而柔性片烟加工运动过程对“松散润叶阶段”至关重要。片烟作为大尺寸平面状的柔性颗粒，与刚性球体黏结烟丝模型相比[16]，几何结构、运动特征及形变方式差异较大，并存在模型复杂、变化多、计算量大且成本高的特点。松散润叶装置一定程度上可以很好对片烟进行松散润叶，但在优化设计过程中，缺乏相关片烟离散元模型参数，无法观察预测滚筒内颗粒运动，产生了设计进度慢、生产调试浪费、无从优化等一系列问题。

本研究以“浙江中烟工业有限公司”松散润叶阶段进料口处片烟为研究对象，采用离散元Rocky DEM 2021R2 仿真软件，创建柔性片状壳体颗粒模型，以堆积角为评价指标。通过对照仿真试验与物理试验的堆积角数值，建立了3 个显著性接触参数与堆积角之间的三元回归方程，模型求解得到片烟与片烟间静摩擦系数（X4）为0.67，片烟与片烟间滑动摩擦系数（X5）为0.31，片烟与片烟间黏附力系数（X6）为0.2，其余非显著性参数取平均值，仿真试验与物理试验的相对误差仅为0.78%，可为进一步计算片烟在松散润叶装置内的运动研究提供理论参考。

需要指出的是，柔性片状模型滚动特征不明显，多以滑动的方式运动输送，因此动摩擦主要考虑滑动摩擦（参考ESSS-Rocky 2021 R2 手册）。离散元法可准确高效地分析颗粒的运动情况，并具有可视化特点，但松散润叶阶段片烟含水率变化对片烟物理运动具有一定影响，本研究仅对松散润叶阶段进料口处含水率12%~15%的片烟进行离散元参数标定，后续将进一步研究不同含水率片烟离散元参数变化规律，以完善松散润叶阶段片烟的离散元仿真数据库。