0 引言
随着固体推进技术的成熟,对于推进剂用氧化剂粉碎粒度以及称量精度要求的提高,以及超细粉碎技术的应用,使得固体推进剂用氧化剂颗粒粒径在2 mm左右[1]。而氧化剂本身不稳定存在安全的风险,并且粉碎后受潮变质导致推进剂性能下降[2],因此对粉碎后的固体推进剂用氧化剂进行自动称量装备设计尤为重要。
目前对于称量设备的研究主要可以分为容积式和重力式。而为了更精确的称量一定质量的氧化剂,本文中研究的主要是重力式称量。并且随着数值模拟技术的发展,可以通过仿真对颗粒的运动进行观测。研究人员在称量设备以及模拟物料在称量的过程上做了许多工作。毛艳[3]通过离散单元法分析了重力式自动称量设备的散粒流动;伍凌川等[4]对动态称重系统进行分析建立了数学模型并且进行了改进,弥补了传统称重系统低效低精度的问题;史慧芳等[5]利用Rocky软件对称料过程进行动态仿真,并通过Ansys Workbench对料斗进行了有限元分析;岳显等[6]使用真空输送以及振动加料称种的方式对物料进行自动定量称重;肖国先[7]使用离散元法对颗粒在料仓中的流动进行数值模拟分析流动规律以及仓壁应力情况;马利英[8]通过CFD-DEM耦合的方法分析了卸料斗中自由下落的微粒流场特性;ZhangYaxiong等[9] 和FengYong等[10]通过DEM方法研究了料仓内谷物颗粒的流动特性以及料仓内物料临界高度与料仓尺寸参数的关系。对于出料口质量流量的研究H.Navid等[11]使用激光线扫描仪测量谷物流动的质量流量。Jochen Mellmann等[12]利用弗劳德数刻画了质量流量并研究了颗粒特性和流动特性对其的影响。Huang Xingjian等[13-15]结合欧拉粒状物质连续体模型、遗传算法以及梯度下降法优化料仓几何参数,最大化地提高了质量流量并且改善了流动模式。
由于传统的出料口开闭形式为开关控制,即称料开始后出料口始终保持打开状态直至装料筒中质量达到传感器设置数值后关闭出料口。虽然称量过程也使用多级补料的方式,但受制于传感器精度以及颗粒在出料口闭合前始终保持流动状态,因此当出料口闭合时仍有一部分物料进入到装料桶中由此造成了称料的误差。
针对上述问题,本文中旨在提出一种采用正弦运动规律的出料口开闭控制形式,通过周期性的出料以更好地保障称料的精度和效率。由于研究的对象为粉碎后的固体推进剂用氧化剂颗粒,为观察颗粒的流动性并通过仿真手段探究出料工艺参数对质量流量以及称量效率的影响,从而解释其合理性,故选择了离散单元法(DEM)对模型进行仿真分析,研究出料口开闭对于质量流量的影响,并规划匹配不同的出料口开闭参数从而有效减少称量时间并降低称量误差。
1 数值模拟方法
1.1 DEM粒子和几何模型
本文中选取的固体推进剂用氧化剂为高氯酸铵(AP)。进入到称量过程的AP一般为经过粉碎处理后的产物,并且粉碎后颗粒为不规则形态,通过对颗粒粒径进行随机测量,获得的结果显示为AP颗粒的粒径分布符合正态分布,其正态分布的期望(μ)为10 μm,标准差(σ)为0.5 μm。但考虑到仿真的时间以及计算机的解算能力,将颗粒视为球形并进行缩放[16-18],放大后颗粒服从正态分布为N(1,0.25),其中μ=1 mm,σ2=0.25 mm2。
料仓的结构如图1所示,其实际的测量方法是通过料仓与楼板的压力传感器进行测量的。通过对料仓以及装料筒进行简化,获得的用于DEM仿真的模型如图2所示。在仿真过程中颗粒的物料参数、颗粒-颗粒间以及颗粒-料仓间的物理参数如表1所示。
表1 颗粒物理参数及其值
Table 1 Physical parameters and their values in simulations
参数数值AP颗粒密度ρp/(kg·m-3)1 950泊松比νp0.4剪切模量Gp/Pa1×108料仓密度ρs/(kg·m-3)7 800泊松比νs0.3剪切模量Gs/Pa7×1010粒子-粒子恢复系数epp0.35静摩擦因数μspp0.1滚动摩擦因数μrpp0.05颗粒-料仓恢复系数eps0.45静态摩擦因数μsps0.1滚动摩擦因数μrps0.01仿真参数时间步长Δt/s5×10-6
图1 称量设备结构示意图
Fig.1 Schematic diagram of the structure of the weighing equipment
图2 称量设备简化模型
Fig.2 Simplified model of weighing equipment
本文中研究的是微细AP,但为了便于仿真实验进行,对颗粒进行了一定的缩放,虽然这样的设置会降低仿真的精度,但实际生产过程中考虑到颗粒间的库仑力、粘性力的影响颗粒也并不是独立且离散的下落,而是以一些不均匀的团簇形式下落。可将细微颗粒组成的团簇视为整体,其中细微颗粒间受到的库仑力、粘性力视为颗粒的内部作用力,由此在仿真过程中忽略颗粒间的粘性。
1.2 离散元法(DEM)
对于粉体流动的研究方法可分为2类:一类是基于连续介质的数值方法;另一类是基于非连续介质的离散粒子近似法[19]。离散元素法是将研究对象视为离散的单元,通过接触模型与运动方程来获取实现对运动情况的预测,使用接触模型获取的力与位移的关系和依靠牛顿第二定律的运动方程来计算求解每一时刻各个颗粒间的接触力和运动参数[20-21]。
本文中使用的仿真分析方法为离散元法(DEM),DEM是一种拉格朗日方法,即粒子的运动可以由6个自由度的运动方程来确定,并考虑到作用在每个粒子上的各种力,由于DEM方法是在较短的时间步长下计算大量粒子间的摩擦碰撞因此需要很大的内存和计算能力。在DEM中对于颗粒间的碰撞可以简化为三维刚性颗粒模型,并且将颗粒间或颗粒与边界的接触表示为振动模型如图3所示。通过接触模型与运动方程来获取实现对运动情况的预测,使用接触模型获取的力与位移的关系和依靠牛顿第二定律的运动方程来计算求解每一时刻各个颗粒间的接触力和运动参数。
图3 接触模型及振动模型
Fig.3 Contact model and vibration model
颗粒接触间的法向振动运动方程式(1)为
(1)
颗粒接触间的切向振动运动可分为切向滑动与颗粒的滚动2个方程式(2)、式(3),分别为
(2)
(3)
由牛顿第二定律可以由位移计算颗粒受到的作用力,得到颗粒i的运动方程式(4)为
(4)
式中:m1,2为颗粒的等效质量;I1,2为颗粒的等效转动惯量;un、us分别为颗粒的法向和切向相对位移; θ为颗粒自身的旋转角度;Cn、Cs分别为接触模型中的法向和切向阻尼系数;kn、ks分别为接触模型中的法向和切向弹性系数;s为旋转半径; Fn、Fs分别为颗粒受到的法向和切向外力分量;M为颗粒受到的外力矩。
2 仿真与分析
2.1 仿真实验
仿真实验设置参数及仿真结果见表2。由于要探究出料口开闭频率对于质量流量的影响,因此在出料口处设置了以正弦运动为开闭形式的挡板。为简化仿真难度降低计算的时间,结合正弦运动规律,出料口在0.5个正弦运动周期完成一次开闭,故分别对单一出料口0.5个运动周期进行离散元仿真,并测量仿真结果从而分别计算单位时间内两出料口的质量流量。
表2 不同出料口频率对应的仿真结果
Table 2 Simulation results corresponding to different outlet frequencies
序号粗称出料口频率/Hz仿真称重值/kg粗称出料口质量流量/(kg·s-1)精称出料口频率/Hz仿真称重值/kg精称出料口质量流量/(kg·s-1)10.1003.2640.6530.5000.1000.10020.3001.2080.7251.0000.0590.11830.5000.7350.7351.5000.0460.13840.7000.5310.7432.0000.0380.15250.9000.420.7562.5000.0300.15061.1000.3520.7743.0000.0230.13871.6000.2450.7843.5000.0180.12682.1000.1910.8024.0000.0140.11292.6000.1560.8114.5000.0120.108103.1000.1240.7695.0000.0100.100113.6000.0970.6985.5000.0090.099124.1000.0780.6406.0000.0080.096134.6000.0640.5896.5000.0070.091145.1000.0530.5417.0000.0060.084155.6000.0450.5047.5000.0050.075166.1000.0380.4648.0000.0050.080
2.2 拟合分析
根据仿真实验测量的结果(见表2)显示随着频率的增长两出料口单次开闭落料的质量(见图4)都随之下降且下降趋势一致,通过对测量的散点进行多项式函数式(5)拟合。
图4 出料口单次开闭的质量曲线
Fig.4The weight curve of single opening and closing of the outlet
y=k0+k1x+k2x2+…+knxn
(5)
其中,k0、k1、…、kn为各项系数。结果显示在给定频率范围内,质量流量与出料口开闭频率的关系为一个四阶的多项式函数。对于两出料口的拟合曲线如图5所示。式(6)和式(7)分别为粗称出料口和精称出料口的频率与质量流量关系的多项式回归函数。
图5 质量流量相对于频率的拟合曲线
Fig.5 The fitting curve of the mass flow rate with respect to the frequency
y=0.000 841 3x4-0.005 054x3-0.029 23x2+
0.145 7x+0.657 3
(6)
y=-0.000 228 2x4+0.005 111x3-0.039 26x2+
0.107 1x+0.052 79
(7)
拟合结果显示,和方差(SSE)分别为2.9×10-3和1.3×10-4,确定系数(R-square)分别为0.983 4和0.986 0,SSE越趋近于0,R-square越趋近于1拟合效果越显著。因此通过数据的拟合回归得到了出料口质量流量与出料口开闭频率的函数关系,以供下一步目标优化进行数据的选择。
由图4和图5中拟合曲线所示,不同出料口质量流量均呈现先增后减的形式,这受到出料口开闭频率的影响。对于出料口开闭频率较低时,虽然单次周期内出料量大的但所耗时间也较长因此换算成质量流量则并不是较高水平;而对于出料口开闭频率较高的情况,由于频繁的开闭使得单次周期内出料量小,即使在很长一段时间内出料量也维持在较低水平,因此质量流量也随着出料口开闭频率的增大而逐渐下降。
3 多目标优化
三级称量优化结构如图6所示。根据上文的描述,为寻找到最佳的称量工艺参数方案,根据实际要求对称量精度和消耗时间这2个指标进行数学建模(见式(8),式(9))。并且希望获得称量精度的最大值和消耗时间的最小值。式(15)描述的是通过两级称量补料后的总质量,为了进一步刻画称量的精度,我们选择将最后1 kg通过两级称量补料完成,并比较称量质量与1 kg的适配度,将其定义为称量精度γ(式(10))。
图6 目标优化关系图
Fig.6 Objective optimization diagram
a1i·m1i+diag{bj1,bj2,…,bjn}·M2=Mi
(8)
a1i·t1i+diag{bj1,bj2,…,bjn}·T2=Ti
(9)
(10)
式中:a1i为粗称出料口选择第i (i=1,…,n)种单次开闭时间下的运动频次,并且将
记为粗称出料口单次开闭时间运动频次向量;diag{bj1,bj2,…,bjn}为对角阵,其对角线元素为与粗称出料口选择对应的精称出料口单次开闭时间下的运动频次,其中i=j;a1i为粗称出料口选择第i(i=1,…,n)种单次开闭时间下的运动频次,并且将记为粗称出料口单次开闭时间运动频次向量;m1i为粗称出料口选择第i (i=1,…,n)种单次开闭时间下的出料质量,并且将
记为粗称出料口单次开闭时间出料质量向量;M2为精称出料口单次开闭时间出料质量向量,
为粗称出料口选择第i(i=1,…,n)种单次开闭形式的时间,并且将
记为粗称出料口单次开闭时间向量;T2为精称出料口单次开闭时间向量,
分别为粗称出料口和精称出料口最小开闭时间的n×1阶列向量,其组成元素的计算方法(式(11));Mi为两级补料完成后的称量质量,是一个n×1阶列向量,将
为n×n阶的矩阵,其中任意元素表示为Mij。
(11)
本文中设置的称量质量达标条件为1 kg±0.005 kg,并且对于A、B中元素必须为整数。通过对粗称出料口质量流量进行选择,并获取单次开闭的落料质量,由于频率的调节精度为0.1 Hz,因此粗选质量流量大于0.735 kg/s的15组频率作为二级补料作业的出料口开闭频率,对于质量流量0.125 kg/s的15组频率作为三级补料作业的出料口开闭频率。
根据图7和图8数据分析结果显示,对于符合称量质量达标条件的向下取整方法满足精度要求的有35组,而对于向上取整方法满足精度要求的有33组。而其中称量质量可以达到千分之一的精度的分别有9组和5组,二级、三级补料的用时为1.51～4.33 s和1.52～2.75 s,数据参数如表3所示。因此对于工艺参数的选择粗称出料口频率为2.5 Hz或2.1 Hz精称出料口频率为1.6 Hz或3 Hz。
表3 目标优化结果参数
Table 3 Objective optimization result parameters
粗称出料口频率/Hz精称出料口频率/Hz称重结果/kg总用时/s粗称出料口频率/Hz精称出料口频率/Hz称重结果/kg总用时/s10.72.80.9993.9482.92.20.9991.6620.73.60.9994.3392.93.80.9991.7431.11.00.9993.41100.52.01.0012.7541.31.60.9991.78110.93.01.0012.2851.33.60.9991.85121.33.01.0011.8261.71.20.9991.59131.93.81.0012.2472.51.60.9991.51142.13.01.0011.52
图7 2种取整方法的精度
Fig.7 Accuracy of the two rounding methods
图8 2种取整方法的用时
Fig.8 Time for two rounding methods
由表3所示,对于粗称出料口与精称出料口处出料口开闭频率的选择多数为图5中所示质量流量较高水平所对应的频率,并且在以两出料口频率可调节基础上选择质量流量最大的情况,匹配获得的称量精度较高,并且完成称量作业所消耗的时间较短,并且根据图7和图8所示图线,在质量流量较高水平所对应的频率所对应的称量精度处于较高水平,称量作业所消耗的时间也较少。
4 结论
1) 使用DEM方法对固体推进剂氧化剂的称量过程进行仿真,研究了二级、三级补料过程出料口开闭形式对质量流量的影响,拟合确定了质量流量与出料口开闭频率的关系。
2) 获取了不同粗称出料口频率与精称出料口频率匹配后的称量精度与用时。结果显示在粗称出料口频率为2.5 Hz或2.1 Hz,精称出料口频率为1.6 Hz或3.0 Hz时所获得的精度达到99.9%,用时分别为1.51 s和1.52 s。
3) 通过数值模拟的方法对于称量工艺参数进行优化选择,减少了实验成本,为实际生产过程中工艺参数的选择提供新的思路。
[1] 邓国栋,刘宏英.超细高氯酸铵粉体制备研究[J].爆破器材,2009,38(1):5-7.
DENG Guodong,LIU Hongying.Study on preparation of the superfinepowder of ammonium perchlorate[J].Explosive Materials,2009,38(1):5-7.
[2] 万雪杰.机械研磨法制备球形超细高氯酸铵及其性能研究[D].南京:南京理工大学,2015.
WAN Xuejie.Study on the preparation and performances of spherical ultrafine ammonium perchlorate by mechanical grinding[D].Nanjing:Nanjing University of Science &Technology,2015.
[3] 毛艳.基于离散单元法的重力式自动称重方法研究[D].武汉:华中科技大学,2017.
MAO Yan.Study on automatic gravimetric weighing method based on DEM[D].Wuhan:Huazhong University of Science and Technology,2017.
[4] 伍凌川,李全俊,黄权.动态称量技术在发射药称重过程中的应用[J].兵器装备工程学报,2017,38(2):70-74.
WU Lingchuan,LI Quanjun,HUANG Quan.Application of dynamic weighing technology in propellant powder weighing[J].Journal of Sichuan Ordnance,2017,38(2):70-74.
[5] 史慧芳,郭进勇,胡翔,等.基于Rocky/Ansys Workbench的发射药自动称量包装过程仿真[J].兵工自动化,2021,40(8):56-60.
SHI Huifang,GUO Jinyong,HU Xiang,et al.Simulation of propellant automatic weighing and packing process based on rocky/ansys workbench[J].Ordnance Industry Automation,2021,40(8):56-60.
[6] 岳显,孔淼.粒状发射药自动定量称重系统[J].兵工自动化,2016,35(10):83-85,93.
YUE Xian,KONG Miao.Automatic quantitative weighing system of granular propellant[J].Ordnance Industry Automation,2016,35(10):83-85,93.
[7] 肖国先.料仓内散体流动的数值模拟研究[D].南京:南京工业大学,2004.
XIAO Guoxian.Numerical simulation study upon granular materials flow in silos[D].Nanjing:Nanjing Tech University,2004.
[8] 马利英.卸料斗几何参数对自由下落微粒流流场特性的影响研究[D].天津:天津商业大学,2015.
MA Liying.Study on effects of discharge hopper geometric parametersflow field characteristics of free-falling particle stream[D].Tianjin:Tianjin University of Commerce,2015.
[9] ZHANG Y X,JIA F G,ZENG Y,et al.DEM study in the critical height of flow mechanism transition in a conical silo[J].Powder Technology,2018,331(Volume).
[10] YONG F,ZIRAN Y.Discrete element method modeling of granular flow characteristics transition in mixed flow[J].Computational Particle Mechanics,2020(8).
[11] NAVID H,TAYLOR R K,YAZGI A,et al.Detecting grain flow rate usinga laser scanner[J].Transactions of the ASABE,2015,58(Volume).
[12] JOCHEN MELLMANN,THOMAS HOFFMANN,CHRISTIAN FüRLL.Mass flow during unloading of agricultural bulk materials from silos depending on particle form,flow properties and geometry of the discharge opening[J].Powder Technology,2014,253(Volume).
[13] HUANG X J,ZHENG Q J,LIU D D,et al.A design method of hopper shape optimization with improved mass flow pattern and reduced particle segregation[J].Chemical Engineering Science,2022,253(Volume).
[14] HUANG X J,ZHENG Q J,YU A B,et al.Optimised curved hoppers with maximum mass discharge rate-an experimental study[J].Powder Technology,2021,377(Volume).
[15] HUANG X J,ZHENG Q J,YU A B,et al.Shape optimization of conical hoppers to increase mass discharging rate[J].Powder Technology,2020,361(Volume).
[16] SHUICHI TANABE,SRIKANTH R.GOPIREDDY,et al.Influence of particle size and blender size on blending performance of bi-component granular mixing:A DEM and experimental study[J].European Journal of Pharmaceutical Sciences,2019,134(Volume).
[17] ALI HASSANPOUR,HONGSING TAN,ANDREW BAYLY,et al.Analysis of particle motion in a paddle mixer using Discrete Element Method (DEM)[J].Powder Technology,2010,206(Volume).
[18] AVIK SARKAR,CARL R.WASSGREN.Effect of particle size on flow and mixing in a bladed granular mixer[J].AIChE Journal,2015,61(Volume).
[19] 乔登攀,孙亚宁,王述红,等.散体移动体积连续性方程研究[J].中国矿业,2005(6):63-66.
QIAO Dengpan,SUN Yaning,WANG Shuhong,et al.Study on the volume continuum equation of moving granular materials[J].China Mining Magazine,2005(6):63-66.
[20] 胡国明.颗粒系统的离散元素法分析仿真:离散元素法[M].武汉:武汉理工大学出版社,2010.
HU Guoming.Analysis and simulation of granular system by discrete element method using EDEM[M].Wuhan:Wuhan University of Technology Press,2010.
[21] HORABIK J
ZEF,PARAFINIUK PIOTR,WI
CEK JOANNA,et al.DEM modelling of the influence of initial stress state on the discharge rate of spherical particles from a model silo[J].Powder Technology,2022,403(Volume).