引用本文

卢绍文, 蔚润琴, 崔玉洁. 磨矿破碎过程粒度分布的分布式参数蒙特卡洛动力学模拟及加速方法. 自动化学报, 2019, 45(9): 1655-1665. doi: 10.16383/j.aas.c180020

LU Shao-Wen, YU Run-Qin, CUI Yu-Jie. A Distributed Parameter Kinetic Monte Carlo Simulation Algorithm of Grinding Process and Its Acceleration. ACTA AUTOMATICA SINICA, 2019, 45(9): 1655-1665. doi: 10.16383/j.aas.c180020

http://www.aas.net.cn/cn/article/doi/10.16383/j.aas.c180020

关键词

磨矿破碎过程，分布式参数模型，粒度分布，蒙特卡洛动力学模拟，模拟加速 

摘要

本文针对磨矿破碎过程，提出一种分布式参数蒙特卡洛动力学方法的粒度分布预测模型和模拟算法.该算法采用了分段思想，将磨机沿着轴向分为若干个虚拟的子磨机；根据破裂、前向和后向移动三类微观事件定义了倾向函数和系统状态矩阵，并设计了分布式算法的调度策略.此外，针对蒙特卡洛动力学算法效率低的问题，提出了基于τ-leap的磨矿过程分布式参数蒙特卡洛模拟加速算法.为了解决分布式参数更新过程中状态不一致的问题，创新性地提出了一种基于缓冲区的同步方法.通过对仿真案例的分析表明，本文提出的分布式参数蒙特卡洛动力学算法具有较高的精度，提出的基于τ-leap的加速算法能够显著提高计算效率，同时保持较好的精度.

文章导读

选矿过程中, 颗粒破碎是在外力作用下使物料由大变小的过程, 使得原料能够满足后续的选别作业对物料的粒度要求, 便于物料中有用矿物与脉石解离.颗粒破碎研磨是高耗能工业过程, 运行优化的主要目标是在满足工艺要求的情况下降低能耗, 并使产品粒度稳定[1-3].粒度分布是破裂过程的关键工艺指标, 它是不同粒度颗粒在总物料中占比.粒度分布过大不能满足工艺要求, 但是如果物料磨得过细会导致能耗过高.粒度分布可以通过激光粒度仪在线检测, 但是这类仪表成本高, 而且受到生产环境干扰影响, 需经常校准.目前通常采用离线筛析, 其主要问题是滞后大, 不能满足在线运行优化的要求.由于破碎过程的主要扰动来自于原料的粒度、硬度的变化, 利用破碎过程模拟技术来预测在当前的生产条件下最终的粒度分布, 可通过运行操作优化来补偿干扰[4-6].此外, 经过验证的模拟模型还可以用于检验工艺设计, 对于优化磨矿生产过程控制和设计都具有重要应用价值[7-10].

破碎过程粒度分布的模拟主要有三类方法.第一类是稳态模拟方法[11], 该方法假设生产的边界条件不变, 将粒度分布的变化近似为线性过程, 主要用于初期的工艺设计, 不能用于在线指导操作.由于磨矿过程机理极为复杂, 难以建立精确的数学模型, 一般采用唯像学建模法[11]. 20世纪40年代, 利用统计原理研究磨机内物料粒度分布, 提出了选择函数和破裂函数两个概念[12].其中选择函数描述颗粒发生破裂的概率, 破裂函数刻画颗粒破碎后得到的子颗粒的粒度分布. Austin等经过实验给出了选择函数和破裂函数的经验公式.文献[13]提出矩阵动力学模型与总体平衡模型, 为磨矿的研究开辟了一个新的领域.文献[14]提出磨矿的动态矩阵理想模型, 进一步推动了磨矿数学模型的发展.基于能级关系的Bond模型是比较经典的球磨机模型, 但这种模型得不到各粒级的质量分数.上述这些模型主要是稳态模型, 无法刻画颗粒破碎过程中粒度分布的演变.

第二类方法为动态模拟方法.磨矿过程的动态模型需要建立粒度分布随时间变化的过程模型, 输入物料的性质、操作条件, 模型输出每一时刻粒度分布的模拟值作为输出.该方法主要采用总量平衡原理建立粒度分布的主方程模型[15].通过磨矿粒度分布的群体平衡动态模型, 记录磨机中不同位置各个粒径的颗粒每时每刻的状态, 从而描述出磨矿过程中物料在时间和空间上的粒度分布.但该模型一般是积分微分方程形式, 多数情况下难以获得解析解[16].

第三类是动力学模拟方法, 它是一种随机模拟方法[17-18], 它用随机过程来刻画颗粒的破碎过程, 能够得到颗粒破碎过程的演化过程.蒙特卡洛动力学模拟方法(以下简称MC方法)是常用的随机模拟方法, 它通过模拟每一次破碎事件的发生, 记录研磨过程中粒度分布的动态变化情况, 较真实模拟了物理过程.文献[17-18]提出了基于MC方法的颗粒破碎过程模拟算法.但是在颗粒破碎过程中, 大粒径颗粒会破碎为多个小粒级的颗粒, 因此随着仿真时间的推进, 需要模拟的颗粒总数以指数速度增加, 导致计算量快速增大, 效率降低.针对这一问题, Khalili等和Smith等提出了一种定数量MC (Constant-number Monte Carlo, CNMC)方法[19-20], 该算法在处理颗粒的凝并和破碎事件时, 通过随机从系统中增加或移除颗粒, 目的是使系统总的颗粒样本数保持不变, 从而提高了原有MC方法的效率.赵海波等提出了一种多重MC (Multi Monte Carlo, MMC)方法, 通过引入加权的"虚拟颗粒"的概念, 将实际颗粒由一颗或者几颗属性相近的虚拟颗粒代表并赋予一个数目转换权值, 该值就是该虚拟颗粒所代表的当地实际颗粒的个数, 模拟时通过调整数目转换权值来保证在发生凝并和破裂事件之后虚拟颗粒总数目仍然保持不变.卢绍文在CNMC基础上[21]提出一种带精度反馈和速度控制的动力学模拟算法(ED-KMC-AC), 能够在满足一个给定的精度阈值约束下, 最大地提高模拟速度.这些方法能够加快MC模拟的速度, 其本质上是通过人为减少系统内颗粒数目实现的, 这将降低模拟的准确程度. Gillepsie在文献[22]中提出了τ-leap近似加速算法, 他假设磨矿过程中多个颗粒破碎信息可以合并, 不需要每一次事件都更新系统, 通过一个调节因子来控制系统更新的粒度, 从而从总体上调节模拟的速度和精度.

上述模型都属于集总参数模型, 即假设磨机是一个质量分布均匀的质点, 物料在磨机内处处均匀.但实际工业生产中, 工业磨机的尺寸日趋大型化, 物料在磨机内的驻留时间更长, 物料性质的空间分布难以忽略, 集总参数模拟方法已经不能准确刻画物料在磨机腔体内的驻留和变化[23].针对这个问题, Kis等[23]提出了分布参数的磨矿过程模型.该模型不仅研究磨机径向方向的破裂, 同时也考虑了磨机轴向方向物料的移动, 考虑了对流前移和扩散后移两种颗粒消失事件.该方法将磨机分为若干段虚拟的子磨机, 从而得到每一段子磨机在连续磨矿过程中物料的变化情况, 最后合并得到产品的粒度分布.但是基于该方法的模型需要采用数值近似方法进行求解, 由于模型方程非常复杂, 不能保证解的存在.此外, 该方法只能解出粒度分布随时间变化的近似解, 无法描述分布参数下磨矿模拟过程中颗粒的破裂和移动的随机特性.针对这些问题, 本文提出一种基于分布式参数破碎过程模拟的粒度分布终点预报方法.我们采用分布参数模拟方法, 将磨机在轴向分段为若干串级连接的虚拟的子磨机, 由于每个子磨机的轴长足够小, 可以假设认为每段子磨机内物料为充分混合状态, 从而将连续磨矿过程看作多个串级连接的批次子磨过程, 但这种方法面临计算代价过大的问题.为了提高模拟效率, 在批次磨矿过程的蒙特卡洛动力学模拟方法基础上, 给出一种模拟加速方法.

图 1  离散径向分布参数磨机模型示意图

图 2  分段子磨机缓冲区示意图

图 3  分布参数MC算法(图中简称MC)和τ-leap加速算法计算得到的平均粒径随时间变化趋势

磨矿是选矿过程中最重要的工艺过程之一, 对于选矿过程的运行品质至关重要.而物料的粒度分布是反映磨矿过程中颗粒破碎情况的直接工艺指标.本文以磨矿过程为背景, 介绍了两种颗粒研磨过程的动力学模拟算法.首先提出一种分布式参数的连续磨蒙特卡洛动力学模拟算法, 进而针对该算法效率较低的问题, 提出了基于τ-leap的模拟加速算法.本文的主要贡献是:

1) 提出分布式参数蒙特卡洛动力学模拟算法.将磨机分为多个串级连接的子磨机.对于每段子磨机, 径向上划分为若干粒级, 假设子磨机内物料均匀混合而且颗粒在子磨机之间移动是在破裂之后发生的.那么子磨机内的破裂过程的模拟采用批次磨的时间驱动蒙特卡洛模拟方法实现.根据分段思想, 相邻的子磨机间有物料的相互传递.因此, 子磨机内颗粒数目的增加是由大颗粒破裂或者子磨机间颗粒的流入而导致的.从微观角度来观察, 磨机内颗粒可能发生三种事件:向前移动事件, 是由对流和物料本身的循环流动引起; 向后移动事件, 是由扩散作用引起; 破裂事件, 是由磨机研磨作用而导致颗粒破裂至比其粒级小的粒级内.将上述三种事件统称为颗粒消失事件.以此为基础定义系统状态矩阵、倾向函数和基于蒙特卡洛的模拟计算方法.

2) 基于τ-leap的磨矿过程分布式参数蒙特卡洛模拟加速算法.分布式参数蒙特卡洛动力学模拟算法能够精确刻画每一个子磨机的动态变化过程, 进而得到整个磨机中不同时间不同位置各个时间点的粒度分布.但由于该种方法需要记录每一次事件的发生, 计算量极大, 只能应用于小体系的模拟.针对这个问题, 我们提出基于τ-leap的模拟加速方法.将模拟时间划分为数个时间片作为步进时间, 同时假设在每个步进时间内有多个反应事件同时发生.即不必每次事件都更新系统状态, 而是根据步进时间片内多个事件的累积效应对系统状态进行更新, 从而提高模拟速度.为了解决状态不一致的问题, 我们创新性地提出了一种基于缓冲区的解决方法.考虑到移动事件和破裂事件直接的独立性, 对每个子磨机系统开辟两个状态同步缓冲区:前向移动缓冲区和后向移动缓冲区.在计算过程中, 将前向移动和后向移动的颗粒分别暂存在缓冲区, 等到达下一个同步时间点后, 再将缓冲区的颗粒与相邻子磨机进行交换.本质上, 这相当于将前/后向的移动事件和破裂事件分别进行处理.相应地, 在进行破裂事件的处理中, 将子磨机和其附属的缓冲区视为一个独立的批次磨过程, 并更新状态转移矩阵和倾向函数.

对本文工作加以总结如下:

1) 本文提出径向离散化的模型框架来解决磨机连续磨过程的分布式参数蒙特卡洛动力学模拟问题, 从仿真结果来看, 算法能够精确地刻画颗粒的随机破碎和径向方向的流动; 但是蒙特卡洛动力学模拟的计算速度远大于差分数值算法.

2) 针对分布式参数蒙特卡洛动力学算法的速度问题, 提出基于τ-leap加速算法来解决这一问题.通过仿真实验以及算法精度与速度的比较分析可知, 基于τ-leap的方法在本质上是一种对蒙特卡洛动力学算法的精度和速度的平衡机制, 通过调节参数ϵ可以实现速度和精度的取舍.

3) 本文所提出的方法, 可以用于球磨机或棒磨机的分布参数动力学模拟.

作者简介

蔚润琴

东北大学信息科学与工程学院硕士生, 2012年获太原理工大学自动化专业学士学位.主要研究方向为蒙特卡洛仿真.E-mail:18842300547@163.com

崔玉洁  

东北大学秦皇岛分校控制工程学院副教授.主要研究方向为机电一体化, 机器人技术.E-mail:cyjhjn@163.com

卢绍文  

东北大学流程工业综合自动化国家重点实验室教授.于2006年获得伦敦大学皇后玛丽学院电子工程学博士学位.主要研究方向为多尺度随机建模方法, 模拟软件设计和数据可视化方法.本文通信作者.E-mail:lusw@mail.neu.edu.cn