摘 要: 针对球团矿生产过程生球堆积排列对球床渗透率、表面积等物理特性的影响, 根据由实际物理现象抽象出的“落滚法则”,考虑生球球体塑性变形和生球粒度分布, 建立了球团矿生球顺次累加堆积的改进数学模型. 对于不同的堆积高度和生球粒度分布进行了生球堆积模拟, 并对球床的空隙率、渗透率和单位体积表面积进行了精确计算.仿真结果表明该模型符合实际生产过程, 为球团矿生球干燥和烧结过程精确模型的建立奠定了坚实的基础.
关键词: 落滚法则；球体堆积；塑性变形；粒度分析
Improved model of iron-ore pellets packing based on dropping and rollingmethodPENG Jun, WANG Jian-hui, GU Shu-sheng, FANG Xiao-ke(College of Information Science and Engineering，Northeastern University，Shenyang 110819，China. Correspondent:PENG Jun，E-mail：impengjun@gmail.com)Abstract:In the iron-ore pellet induration process, the green ball packing influences the porosity, the permeability, and thesurface area of pellet bed. Take the plastic deformation and size distribution of green balls into consideration, an improvednumerical model for sphere packing process is presented by using dropping and rolling rules which are derived from a realisticphysical phenomenon. The simulation is taken for different packing heights and size distributions, and the correspondingaccurate porosity, permeability and surface area for normal distribution of sphere radii can be got. The simulation resultshows that the model is efficient and well reflects realistic packing process, and it also provide the foundation of accuratedrying and induration model.
Key words:dropping and rolling method；sphere packing；plastic deformation；size distribution
    0 引 言
    利用链篦机-回转窑球团法可将造球机造好的生球进行干燥预热、焙烧固结和冷却. 首先, 布球机将生球布在慢速运行的蓖板上, 利用环冷机三冷段余热及预热二段排除的热气流对生球进行鼓风机干燥及抽风干燥; 然后由二冷段和回转窑的热气进行预热氧化, 脱除吸附水或结晶水, 并达到足够的抗压强度后直接送入回转窑进行焙烧. 因为回转窑焙烧温度高且回转, 所以加热温度均匀, 不受矿石种类的限制, 并且可以得到质量稳定的球团[1].冶金领域的学者对球团矿烧结过程模型进行了大量的研究[2-6]. Voskamp等[2]首次提出了烧结过程模型; Young等[3]构建了链篦机-回转窑球团法的完整模型; Thurlby[4]将链篦机-回转窑过程分解为几个子过程, 并通过气体流速和气压变化分析了这几个子过程的联系. 相对于链篦机-回转窑球团法, 带式焙烧机球团法的研究文献较多. Barati[5]提出了完整的带式焙烧机过程模型, 包括干燥、烧结和冷却过程, 并考虑了球床的收缩率; Majumder[6]在现有的模型基础上进行了部分改进, 开发出带式焙烧机模拟器, 对烧结过程进行模拟并对一些关键参数进行了预测.
    在上述球团矿烧结过程模型中, 球床空隙率一般取为定值, 模型中对于球床的空隙率和单位体积表面积没有严格的计算. 实际上, 力学和流体动力学研究表明, 球体堆积排列方式等内部微观因素的变化对造粒系统的许多宏观的物理特性(如渗透率等)影响较大, 因而本文对生球堆积排列方式进行了微观研究,由实际生产数据得出生球粒度分布; 考虑到生球球团这种颗粒聚合体在堆积过程中发生的塑性变形, 在由实际物理现象中抽象出的“落滚法则”基础上, 得出球体顺次累加算法; 通过仿真实验模拟各个因素对球床的空隙率和渗透率的影响, 进而对精确的球团矿烧结过程模型的建立提出了理论和计算依据.
    1 球球球体体体堆堆堆积积积排排排列列列
    力学研究表明, 在重力场中等粒度球体的堆积排列只有4种方式: 立方体、正交、四方-楔形和菱面体,且以球体的球心互相连接所构成的立体图形来命名.球体堆积排列方式的参数主要有空隙率和配位数. 空隙率定义为空隙的体积占整体的比例, 其计算公式为𝜑=𝑉0𝑉=𝑉−𝑉1𝑉. (1)其中: 𝜑为空隙率, 𝑉0为空隙体积, 𝑉为整体体积,𝑉1为球体体积, 𝑑为直径. 以立方体堆积为例, 空隙体积为单位立方体减去8个顶角上的八分之一球体的体积, 即𝜑=(𝑑3−43π(𝑑2)3)/𝑑3≈0.4764.此外, 也常用堆积密度来描述堆积形态, 计算式为𝑉1/𝑉.配位数定义为与特定的球体相互接触的球体的个数, 如立方体堆积中, 每个球体的配位数为6.不难得出其他堆积方式的空隙率和配位数, 如表1所示.表1 堆积排列的空隙率和配位数堆积方式 配位数 空隙率/%立方体 6 47.64正交 8 39.54四方-楔形 10 30.19菱面体 12 25.95球团矿干燥烧结过程中, 造球机生产出的生球由布料机放置于蓖板上, 堆积成一定高度的多层的生球球床. 除了实验室测量外, 实际生产中无法对球床的堆积形态和空隙率进行测量和研究计算.
    2 基基基于于于顺顺顺次次次累累累加加加的的的球球球体体体堆堆堆积积积改改改进进进算算算法法法
    现有的堆积模拟算法大致可分3类: 1) 非顺次重排列, 球体无需按顺序排列, 尽可能地增加其堆积密度; 2) 动态模拟, 该方法以离散元法为基础, 将堆积过程视为动态过程, 并考虑了球体间的作用力; 3) 顺次累加[7], 在重力场设定尺寸的容器中依次落入球体.由于顺次累加最接近真实的物理现象, 本文在顺次累加算法思想的基础上, 以落滚法则作为判稳依据, 考虑生球粒度分布和塑性变形两大重要影响因素, 建立了球团矿生球堆积的精确模型.
    2.1 基基基于于于落落落滚滚滚法法法则则则的的的顺顺顺次次次累累累加加加落滚
    法则是球体落下堆积时的一种运动规则, 包含落下法则和滚动法则两部分. 在三维坐标系中,为原点, 该三维坐标系位于重力场中. 平面随机落下粒度符合一定分布的球体, 其球心坐标为(,,), 其中(,)坐标随机产生. 球体在落滚法则下进行堆积. 首先, 球体在重力作用下沿轴反向落下, 启动落下法则, 即球心坐标沿轴按计算步长减小, 直到与面接触(坐标等于球体半径), 球体稳定; 当落球与别的球体接触时, 则按照滚动法则达到稳定. 记录稳定状态下球心的三维坐标和粒度数据,循环顺次累加, 直到满足终止条件, 即球床高度达到设定值max, 程序终止, 输出所有的数据. 顺次累加算法流程如图1所示.YN!"#$%&( , ) x y'()*7- 89 1 ?7- 89 3 ?7- 89 2 ?7 89 xoy+ ?7 89 xoy+ ?7 89 xoy+ ?,- ./0 1,- ./0 2123456Zmax?NNYYNNNNYYYY!"#$图1 顺次累加算法流程对于落球没有直接接触平面的情况, 在此运用滚动法则进一步详细分析落球的状态. 在三维坐标俯视图(图2)中, 稳态球体用虚线表示, 坐标和半径分别为(1,1)和𝑟1; 待稳定的球体用实线表示,坐标和半径分别为(0,0)和𝑟0.30201000 10 20 30 xyxy( , ) x y1 1( , ) x y0 0r1r0o图2 球体堆积形态示意图如果两者平面圆心距小于半径和, 即(1−0)2+第4期 彭 俊等: 基于落滚法则的球团矿生球堆积改进模型 549(1−0)2<(𝑟1+𝑟0)2, 则两球处于接触状态. 只有与相邻3个稳定球体相接触, 落球才会稳定. 落球首先与球1接触后并不稳定, 沿着与垂直且包含落球与球1球心间连线的平面滚动, 若没有与面接触,则沿着如图3(a)所示方向滚动, 直到与球2接触.(a) 1 2 !"#" $%&'("!""2"1"3z’x’y’(b) 1 !"#" 、" * 2 $%!""2"1"3z’x’y’(c) !"#+"$%!""2"1"3zxyo图3 落球状态为便于计算, 以球1的圆心为原点, 建立一个与原始三维坐标平行的新坐标′′′, 不难得出落球的位置坐标(′0,′0,′0)为⎧⎨⎩′0=1+ (𝑟1+𝑟0)cos𝜃cos𝜑,′0=1+ (𝑟1+𝑟0) cos𝜃sin𝜑,′0=1+ (𝑟1+𝑟0)sin𝜃.(2)其中: 𝜃为球心连线与′轴的夹角, 𝜑为球心连线在′′平面上的投影与′轴的方位角.落球与两球接触后, 若没有与面接触, 则将继续沿着与包含3个球心的平面相垂直的平面的方向滚动. ′、1和2分别为落球、球1和球2的球心,′𝐵为球心滚动的弧形轨迹, ′𝐴为12的垂线, ′𝐴与𝐴𝐵的夹角为𝛼. 𝐵为落球的新的球心位置, 坐标可以由向量计算, 即⎧⎨⎩𝐵=𝐴+𝐴𝐵,𝐴=1+12⋅1′⋅ 12∣122∣,𝐴𝐵=𝐴′+𝐴′×12∣𝐴′×12∣∣𝐴′∣tan𝛼.(3)最终, 落球与3个稳定球体相接触, 达到稳定状态. 在累加过程中, 采用周期性边界方法消除球床边界的壁面效应, 即当 <0时, =+max; 当 > max时, =−max. 坐标同理可得.
    2.2 生生生球球球粒粒粒度度度分分分布布布
    在以往的球团矿干燥或者烧结模型中, 生球一般被假设为等粒度球体, 在此假设基础上计算出的球床空隙率或单位体积表面积与实际情况相差较大, 降低了后续模型的精确度.实际生产过程中, 造球机生产出的生球并不是单粒度的. 粒度既要满足高炉冶炼的需要, 同时也要考虑造球设备的条件和生产量. 一般高炉用的球团粒度中, 6∼16mm粒级的占80%∼95%[1]. 根据文献[8]提供的5.5小时获得的23个批次的数据, 筛分粒度从小到大分为10组. 原始数据中测量的是每个筛分粒度组中的质量百分比, 由粒度筛分组中生球的平均质量和造球机的产量, 可计算出23个批次、10个筛分组的生球个数.根据计算结果产生的数据样本, 对其概率分布进行分析, 从而得出造球机生产出的生球的粒度分布,虽然𝑃值小于0.005,但考虑到样本点中除两端的点偏离实际分布外, 大部分都符合分布. 在实际造球中,两端的样本点对应粒度过大或过小的生球都属于不合格产品, 会重新回收造球, 故两端可以视为离群值,不予考虑[9]. 因此在实际生产过程中, 生球的粒度符合近似正态分布, 本文采用的粒度数据符合正态分布𝑁(12.11,1.3).
    2.3 生生生球球球的的的塑塑塑性性性变变变形形形
    对于球团这种颗粒聚合物, 在堆积过程中重力和挤压力的作用下势必发生一定的塑性变形. 如图4所示, 塑性变形会导致球度变化, 球团间会有部分重叠进而导致空隙率的变化, 因此从实际情况出发, 将塑性变形因素引入模型中.o1o2oR1R2B2B1a34 56786412340xyz图4 互相接触的球团间的塑性变形图4中: 𝑎为接触区域半径, 𝐵1和𝐵2分别是球体1和2的塑性变形量, 𝑅1和𝑅2分别是接触球体的半径, 12是球心距. 球床由多层生球堆积而成, 为了避免堆积过程中上层生球的重力作用以及鼓风段的气流作用而导致生球破裂, 球床的高度一般控制550 控 制 与 决 策 第28 卷在18∼22cm范围内. 由弹塑性力学的知识可知, 生球塑性形变量与其所处的高度成反比.生球接触区域半径与其坐标的关系可表示为=−1𝑎2+2𝑎+3. (4)其中: 1和2为正实数; 3为球床高度上限, 取值范围为18 ∼22cm;破裂零界点处即球床高度为22 cm时, 最底层生球塑性形变量为1 ∼1.5mm[10],约为球体半径的0.1倍, 则不难计算接触区域半径约为球体半径的0.43倍.
    3 堆堆堆积积积过过过程程程模模模拟拟拟和和和结结结果果果分分分析析析
    模拟程序采用Fortran语言和NAG公司FortranBuilder 5.2编译器. 为选取适当的堆积尺寸, 首先进行单粒度刚性球体堆积模拟, 堆积高度采用球团烧结过程球床高度的下限, 即18 cm,长度和宽度选取5组数据, 计算其对应的空隙率如表2所示.表2 不同底面面积, 单粒度堆积下的空隙率底面面积/cm2252302352402502空隙率/% 42.57 42.38 42.36 42.37 42.36由表2数据可知, 随着底面面积的逐渐增大, 空隙率会变小, 直到面积达到35×35后, 空隙率稳定在42.36%,与文献[11]报道的堆积体积密度为0.58相近. 为避免随机性对计算结果的影响, 对应每组三维尺寸程序运行3次, 结果显示空隙率的误差在0.04%以内. 因而在随后的实际模拟中, 容器底部的尺寸分别设为max= 30cm,max= 30cm.
    3.1 不不不同同同球球球床床床高高高度度度下下下的的的生生生球球球塑塑塑性性性变变变形形形
    由式(4)提供的球床高度和生球塑性形变的关系, 以及2.2节统计结果, 模拟粒度均值为12.11 mm,标准差为1.3 mm的正态分布球体堆积, 并考虑生球的塑性变形, 破裂零界点处塑性变形接触区域半径取值为0.43𝑅. 图5分别为不同球床高度下, 堆积形态3D图和生球间的应力网络图. 随着堆积高度的增加,下层球体间的应力不断加大, 塑性变形加剧. 图中表现为3D图中下层球体与接触球间的重叠区域增大;而应力网络图中, 应力线更加密实. 由此可见, 本模拟算法有效地反映了不同堆积高度下生球塑性变形的情况.
    3.2 不不不同同同生生生球球球粒粒粒度度度范范范围围围下下下的的的堆堆堆积积积模模模拟拟拟
    根据文献[1]提供的生球参考粒度11 mm以及2.2节的统计结果, 模拟程序将其分别设为均值, 对不同的粒度范围分布和不同的球床高度, 计算其对应的空隙率和渗透率, 如表3所示.空隙率和渗透率的关系由经典Carman-Kozeny方程描述, 即10100020303020100(c) =18 cm z !"#$20x yz1010002030302010(d) =18 cm z !%&'($20x yz010100020303020100(a) =10 cm z !"#$zy x10100020303020100(b) =10 cm z !%&'($zy x10100020303020100(e) =20 cm z !"#$20x yz1010002030302010(f ) =20 cm z !%&'($20x yz0 图5 不同球床高度下的堆积图和应力网络图𝜅=𝑑2𝑠180⋅𝜑3(1−𝜑)2. (5)其中: 𝑑𝑠 =w𝑑3𝑓(𝑑)d(𝑑)w𝑑2𝑓(𝑑)d(𝑑)为表面平均直径, 𝜅为渗透率(mm2), 𝑓(𝑑)为球体粒度的分布函数. 由积分函数的定义, 可转换为𝑑𝑠 =∑𝑖𝑑3𝑖/∑𝑖∑𝑖𝑑2𝑖/∑𝑖. (6)空隙率和单位体积表面积的关系由下式表示:𝐴= 6(1−𝜑)/𝑑𝑠. (7)第4期 彭 俊等: 基于落滚法则的球团矿生球堆积改进模型 551表3 不同粒度范围下的空隙率和渗透率粒度均值 粒度范围 高度 空隙率 d𝑠 渗透率 单位体积表面积mm mm cm % mm mm2m−115 40.04 12.24 0.148 6 29.3911 6.0∼16 18 39.71 12.01 0.138 0 30.1220 39.55 11.95 0.134 3 30.3515 40.09 11.87 0.140 5 30.2811 7.0∼15 18 40.01 11.79 0.137 4 30.5320 39.96 11.67 0.133 9 30.8715 40.29 11.54 0.135 7 31.0511 8.0∼14 18 40.06 11.42 0.129 6 31.5020 39.99 11.35 0.127 1 31.7315 41.07 12.41 0.170 7 28.4912.1 9.5∼14.7 18 40.15 12.35 0.153 1 29.0720 40.28 12.28 0.153 5 29.18由表3可知: 相同的粒度范围内, 随着球床高度的增加, 生球的塑性变形加剧, 空隙率变小, 渗透率变小, 表面积增大; 而在相同的球床高度情况下, 粒度范围越宽, 空隙率越小, 渗透率也变小. 可见, 渗透率受粒度均值的影响更大一些. 表面平均粒度与粒度均值正相关, 与球床高度无关; 而粒度范围越宽, 表面平均粒度越大. 单位体积表面积更多地受表面平均粒度的影响. 因而在实际生球干燥过程中, 为了获得较好的球床渗透率, 造球机的生球粒度设定值在满足其他参数的情况下应稍高一些.
    4 结结结 论论论
    本文根据实际物理现象抽象出“落滚法则”,并考虑实际球团生产过程中生球的粒度符合近似正态分布, 以及堆积过程中球体产生的塑性变形, 对生球堆积过程进行了仿真模拟, 得出了球床精确的空隙率、渗透率和单位体积表面积, 这对后续的生球干燥以及烧结过程模型的建立奠定了坚实的基础.
    参考文献(References)[1] 张一敏. 球团矿生产技术[M].北京: 冶金工业出版社,2005: 69-70.(Zhang Y M. Technology of iron ore pelletsmanufacture[M]. Beijing: Metallurgical Industry Press,2005: 69-70.)[2] Voskamp J, Brasz J. Digital simulation of the steady statebehaviour of moving bed processes[J]. Measurement andControl, 1975, 8: 23-32.[3] Young R, Cross M, Gibson R. Mathematical model ofgrate-kiln-cooler process used for induration of iron orepellets[J]. Ironmaking Steelmaking, 1979, 6(1): 1-13.[4] Thurlby J. A dynamic mathematical model of thecomplete grate/kiln iron-ore pellet induration process[J].Metallurgical and Materials Transactions B, 1988, 19(1):103-112.[5] Barati M. Dynamic simulation of pellet induration processin straight-grate system[J]. Int J of Mineral Processing,2008, 89(1/2/3/4): 30-39.[6] Majumder S, Natekar P, Runkana V. Virtual indurator: Atool for simulation of induration of wet iron ore pellets ona moving grate[J]. Computers and Chemical Engineering,2009, 33(6): 1141-1152.[7] Visscher W M, Bolsterli M. Random packing of equaland unequal spheres in two and three dimensions[M]. NewMex: Los Alamos Scientific Lab, 1973: 98-99.[8] Thurley M J, Andersson T. An industrial 3D visionsystem for size measurement of iron ore green pelletsusing morphological image segmentation[J]. MineralsEngineering, 2008, 21(5): 405-415.[9] Navidi W C. Statistics for engineers and scientists[M]. NewYork: McGraw-Hill, 2004: 251-254.[10] Forsmo S P E, Apelqvist A J, Bj¨orkman B M T, et al.Binding mechanisms in wet iron ore green pellets witha bentonite binder[J]. Powder Technology, 2006, 169(3):147-158.[11] Tory E, Church B, Tam M, et al. Simulated randompacking of equal spheres[J]. The Canadian J of ChemicalEngineering, 1973, 51(4): 484-493.