钢铁冶金的能耗约占工业一次能源消费总量的17%,约90%的钢铁由高炉-转炉系统生产,铁前系统（高炉）能耗与污染物排放至少占钢铁生产能耗的70%^[1]。在推进“双碳”目标实现的进程中,中国炼铁工业面临着节能减排的重大挑战,传统高炉炼铁流程的碳利用和碳减排空间已几乎达到极限,因此,必须开展各项新的低碳炼铁工艺研究,以提高碳利用率,并降低碳排放^[2]。在此进程中,由传统高炉向低碳高炉转换的炉况过渡过程中开展炉内热化学行为演变规律研究,对炼铁技术绿色低碳发展和实现高炉冶炼的减污降碳目标具有重要的推动作用。

炼铁高炉内部炉料颗粒流、气体流和液体流等多相物质流之间在高温高压环境运动过程中会发生强烈的同向或逆向热质传输,同时持续发生相变转化,因此,炉内传输现象异常复杂,使得直接观察和测量高炉内部信息非常困难。随着计算机技术的不断发展,基于物理机理的数值模拟仿真已成为认知、解析和优化调控高炉运行的一种重要手段^[3]。由于高炉冶炼过程通常伴随着不同时空尺度的传输现象,因而需要运用不同尺度方法建立和解析高炉过程^[4-5]。相较于广泛用于描述多相流传输过程的双流体模型（TFM）^[6],计算流体力学和离散单元法（CFD-DEM）的耦合方法能够解析颗粒尺度状态,可以提供大量微观信息（如单个颗粒的受力、运动轨迹、传热和化学反应特征等）,以进一步认识颗粒基多相物质流传输过程的内在机理。基于该方法可以开发出更完善合理的高炉数学模型,用于分析解决高炉运行过程中出现的问题,并对其进行优化控制。因此,越来越多的研究者采用CFD-DEM耦合方法建立的基于物理机理的数学模型研究高炉炼铁过程^[7-8]。采用球形和非球形的DEM方法对4 000 m³高炉的布料过程进行数值模拟^[9],发现多面体炉料接触到溜槽后会比球形炉料损失更多的动量,导致落点更靠近炉心,而且多面体炉料的形貌分布、层间渗流、料层稳定性和孔隙率分布都不同于球形炉料。LIU等^[10]开发了一个基于DEM的高炉布料模型,研究了溜槽角和摩擦系数等关键变量对料层分布的影响,结果表明,通过优化这些变量可以获得更合理的料层结构。LI等^[11]联合TFM高炉炉身模型和DEM布料模型提出一种高炉上部料层分布的优化策略,该策略可在保持高炉运行稳定的同时提高冶炼效率。文献[12]通过建立氧气高炉操作条件下的三维数学模型,分析了氧气高炉下部的温度场、浓度场及煤粉的流动和燃烧特性,结果表明,氧气高炉操作条件下的回旋区温度显著升高,高温区面积扩大,CO₂含量提高,焦炭床内CO含量显著增加。

目前,关于布料对高炉炉况影响的研究多集中于稳定状态条件下料层结构分布对高炉运行过程参数的影响,以及炉顶布料和局部料层分布及形成机理。对于从颗粒尺度解析料批重对改变料层结构的作用,进而影响高炉内部热化学行为的动态演变规律,目前还未见报道,而炉况动态改变引起的炉内流动和传热传质现象的演化研究对由传统高炉操作向低碳冶炼转变具有重要的指导意义。本研究采用CFD-DEM耦合方法建立了颗粒尺度高炉数学模型,用于研究料批重转换对高炉内部传输过程的动态影响,为实际的高炉料批重调整提供了参考。

1   模型描述

采用基于CFD-DEM耦合方法建立并经验证的瞬态颗粒尺度高炉数学模型^[3,13-15],考虑流体（气相）和颗粒相（矿石和焦炭颗粒）流动、传热传质,以及主要的化学反应,简要介绍如下。

1.1   颗粒相控制方程

固体颗粒的平移和旋转运动由DEM方法^[12]描述,颗粒i与其相邻的颗粒和/或墙壁之间发生相互作用,发生动量和能量（导热传热）交换,同时与周围的流体进行热量传输（对流传热和辐射传热）和发生化学反应进行质量传输,具体控制方程如式（1）—式（4）所示：

式（1）—式（4）中：$m_i$ 为颗粒i的质量；v_i 为颗粒的平移速度；ω_i为颗粒的旋转速度；f_e,ij和f_d,ij分别表示颗粒之间的弹性力和黏性阻力；f_pf,i表示颗粒和流体之间的相互作用力；m_ig为颗粒i的重力；I_i为颗粒的转动惯量$(I_i=\mathrm{2}/\mathrm{5}{m}_i{R}_i^{\mathrm{2}})$；T_t,ij和T_r,ij分别代表颗粒j作用在颗粒i上力矩的2个分量,前者由切向力产生,使颗粒i旋转,后者为滚动摩擦力矩；$c_{\mathrm{p},i}$为颗粒i的比热容；${\dot{Q}}_{i,j}$和${\dot{Q}}_{i,\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}$分别为颗粒i和颗粒j及壁面导热传热量；${\dot{Q}}_{i,\mathrm{f}}$和${\dot{Q}}_{i,\mathrm{r}\mathrm{a}\mathrm{d}}$分别为颗粒i与气体之间的对流传热量和周围环境的辐射传热量；${\dot{Q}}_{i,\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}}$是为化学反应热；c_i,m为颗粒中组分m的浓度；s_i,m为颗粒中组分m的反应速率。

1.2   气体相控制方程

对于气体相传输过程采用CFD方法建模,通过采用平均化的Navier-Stokes方程进行求解^[7],各守恒方程如式（5）—式（8）所示：

式（5）—式（8）中：ρ_f为气体密度；ε_f为气体体积分数；u为气体速度矢量；S_mf为由于化学反应产生的气体组分源相；p为压力；F_fp为单位体积内气体颗粒间相互作用力（F_fp $={\displaystyle \sum _{i=\mathrm{1}}^{k_V}}({\mathit{f}}_{d,i}+{\mathit{f}}_{pg,i})/\mathrm{\Delta }V$）；τ为应力张量；g为重力加速度；T_f为气体温度；c_f为气体比热容；k_e为有效导热系数；$\dot{Q}$为单位体积内气体与颗粒、壁面及环境辐射的换热量（$\dot{Q}=({{\displaystyle \sum _{i=\mathrm{1}}^{k_{\mathrm{V}}}}\dot{Q}}_{\mathrm{f},i}+{\dot{Q}}_{\mathrm{f},\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}+{\dot{Q}}_{\mathrm{f},\mathrm{r}\mathrm{a}\mathrm{d}})/\mathrm{\Delta }V$）；${\dot{Q}}_{\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}}$为化学反应热（${\dot{Q}}_{\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}}=∆H_{\mathrm{298}}^{°}+{\int }_{\mathrm{298}}^T∆C_{p}dT$）；C_m为化学组分m的浓度分数；S_m为化学反应源相；Γ_m为组分m的扩散系数。

1.3   主要化学反应

铁矿石的还原反应及焦炭的气化反应基于单界面未反应收缩核模型（USCM）^[16]。对于铁矿石的还原反应,该界面考虑了3种还原过程：Fe₂O₃→Fe₃O₄→Fe_(1-x)O→Fe,还原反应速率的计算取决于铁矿石的热力学及动力学条件,考虑的主要化学反应和具体反应速率公式如表1所列^[17-21]。

表  1  化学反应和反应速率公式

Table  1.  Chemical reactions and reaction rate equations

化学反应	反应速率	参考文献
$\frac{1}{3}\mathrm{F}{\mathrm{e}}_2{\mathrm{O}}_3（\mathrm{s}）+\mathrm{C}\mathrm{O}（\mathrm{g}）=\frac{2}{3}\mathrm{F}\mathrm{e}（\mathrm{s}）+\mathrm{C}{\mathrm{O}}_2（\mathrm{s}）$	$R_1^{*}=\frac{\mathrm{\pi }{d}_0^2{\phi }_0^{-1}{N}_0\cdot 273p(y_{\mathrm{C}\mathrm{O}}-y_{\mathrm{C}\mathrm{O}}^{*})/(22.4T_{\mathrm{p}})}{0.5d_0[(1-f_{\mathrm{s}}{)}^{-1/3}-1]/D_{\mathrm{s}1}+[(1-f_{\mathrm{s}}{)}^{2/3}{k}_1(1+1/K_1{)]}^{-1}}$	[17]
$\frac{1}{3}\mathrm{F}{\mathrm{e}}_2{\mathrm{O}}_3（\mathrm{s}）+{\mathrm{H}}_2（\mathrm{g}）=\frac{2}{3}\mathrm{F}\mathrm{e}（\mathrm{s}）+{\mathrm{H}}_2\mathrm{O}（\mathrm{g}）$	$R_2^{*}=\frac{\mathrm{\pi }{d}_0^2{\phi }_0^{-1}{N}_0\cdot 273p(y_{{\mathrm{H}}_2}-y_{{\mathrm{H}}_2}^{*})/(22.4T_{\mathrm{p}})}{0.5d_0[(1-f_{\mathrm{s}}{)}^{-1/3}-1]/D_{\mathrm{s}2}+[(1-f_{\mathrm{s}}{)}^{2/3}{k}_2(1+1/K_2{)]}^{-1}}$	[17]
$\mathrm{C}（\mathrm{s}）+\mathrm{C}{\mathrm{O}}_2（\mathrm{g}）=2\mathrm{C}\mathrm{O}（\mathrm{g}）$	$R_3^{*}=\frac{\mathrm{\pi }{d}_{\mathrm{c}}^2{\phi }_{\mathrm{c}}^{-1}{N}_{\mathrm{c}}\cdot 273p{y}_{\mathrm{C}\mathrm{O}}/(22.4T_{\mathrm{p}})}{1/k_{\mathrm{f}3}+6/(d_{\mathrm{c}}{\rho }_{\mathrm{c}}{E}_{\mathrm{f}}{k}_3)}$	[18]
$\mathrm{C}（\mathrm{s}）+{\mathrm{H}}_2\mathrm{O}（\mathrm{g}）=\mathrm{C}\mathrm{O}（\mathrm{g}）+{\mathrm{H}}_2（\mathrm{g}）$	$R_4^{*}=\frac{\mathrm{\pi }{d}_{\mathrm{c}}^2{\phi }_{\mathrm{c}}^{-1}{N}_{\mathrm{c}}\cdot 273p{y}_{{\mathrm{H}}_2}/(22.4T_{\mathrm{p}})}{1/k_{f4}+6/(d_{\mathrm{c}}{\rho }_{\mathrm{c}}{E}_{\mathrm{f}}^{\mathrm{\text{'}}}{k}_4)}$	[18]
$\mathrm{F}\mathrm{e}\mathrm{O}（\mathrm{l}）+\mathrm{C}（\mathrm{s}）=\mathrm{F}\mathrm{e}（\mathrm{l}）+\mathrm{C}\mathrm{O}（\mathrm{g}）$	$R_5^{*}=k_5{C}_{\mathrm{F}\mathrm{e}\mathrm{O}}^2{A}_{\mathrm{c}}/V_{\mathrm{B}}$	[19]
$\mathrm{S}\mathrm{i}{\mathrm{O}}_2（\mathrm{l}）+2\mathrm{C}（\mathrm{s}）=\mathrm{S}\mathrm{i}（\mathrm{l}）+2\mathrm{C}\mathrm{O}（\mathrm{g}）$	$R_6^{*}=k_6{C}_{\mathrm{S}\mathrm{i}{\mathrm{O}}_2}{A}_{\mathrm{c}}/V_{\mathrm{B}}$	[20]、[21]

下载: 导出CSV 

| 显示表格

1.4   模拟条件

根据LKAB的8 m³工作容积实验高炉的设计参数构建了相应的三维薄片模型（图1（a））^[22]。前后薄片模型壁面采用周期性边界条件,保证生成更合理的颗粒床层堆积孔隙,厚度设置为4倍颗粒粒径大小。对于气体相计算,采用的结构化网格如图1（b）所示。通过质量守恒和热平衡计算给定初始时刻风口入口气体组分和气体温度。高炉内部铁矿石层和焦炭层交替分布,详细物理性质如表2所列。在平衡模型相对准确度和计算效率情况下,适当放大了实际高炉过程采用的矿石和焦炭颗粒粒径^[15]。依据实际LKAB高炉操作,本研究中的高炉料线设置在6.40 m处。根据指定的料批重,交替将铁矿石、焦炭由高炉顶部装填,模拟计算过程中实现矿焦分层结构及保证稳定的料线位置,如图1（c）所示。在模型模拟计算时,当铁矿石颗粒被完全还原且其温度高于熔化温度时（本文假设为1 400 ℃）将被排出。为减少达到高炉稳定状态的模拟计算时间,本研究采用了结合高炉过程特征的相似理论加速计算方法^[15],在一定程度上缓解了DEM方法高计算资源需求的压力,同时参考实验结果设置了固体颗粒床层的初始温度场（图1（d））。

图  1  模拟条件：（a）高炉几何；（b）CFD网格划分；（c）初始矿焦交替分布；（d）初始温度场

Figure  1.  Simulation conditions: （a） blast furnace geometry; （b） CFD grid division; （c） alternate distribution of initial ore coke; （d） initial temperature field

下载: 全尺寸图片 幻灯片

表  2  气体、铁矿石和焦炭颗粒的物理性质

Table  2.  Physical properties of gas, iron ore and coke particles

变量	数值
炉内颗粒数量	12 070
矿石/焦炭粒径	0.04 m
矿石/焦炭密度	3 500/1 081 kg/m3
矿石/焦炭批重	207/64 kg, 163/50 kg
颗粒-颗粒/壁面滑动摩擦	0.40/0.01
矿石/焦炭颗粒杨氏模量	1.0×107 kg/（m∙s2）
颗粒泊松比	0.30
恢复系数	0.80
入口速度	70 m/s
流体密度	0.14 kg/m3
流体分子黏度	1.80×10-5 Pa⋅s

下载: 导出CSV 

| 显示表格

参考LKAB实验高炉料批重设置^[7,22-23],本研究采用2种布料批重：（B1）矿批207 kg,焦批64 kg；（B2）矿批163 kg,焦批50 kg。主要考察布料批重2种动态转变过程（State Transition 1 （ST1）：B1→B2→B1和State Transition 2 （ST2）：B2→B1→B2）中高炉内部的热化学行为的演变规律。在模拟操作过程中,首先模拟计算给定料批重的高炉运行过程,待炉况达到稳定状态时改变其料批重设置,计算稳定后,再恢复到原来的料批重设置,直至再次达到稳定状态。

2   结果与讨论

2.1   模型验证

采用的CFD-DEM高炉模型是本研究组基于编程语言自主开发的,并已在先前的工作^[13-15]中验证,同时与WATAKABE等^[22]解剖实验高炉所获软熔带形状位置和实验测量的炉内温度分布进行对比验证（图2）。由图2可见,此模型预测结果与实验结果较一致,表明本研究建立的模型具有可用性。

图  2  模型验证:（a）实验高炉解剖的软熔带^[22];（b）模型预测的软熔带;（c）模型预测的炉料温度分布及探针布置;（d）模型预测和实验测量^[22]的温度分布对比

Figure  2.  Model validation: （a） softening zone dissected by experimental blast furnace^[22]; （b） soft melting zone predicted by the model; （c） temperature distribution of furnace burden predicted by the model and probe arrangement; （d） comparison of temperature distribution between model prediction and experimental measurement^[22]

下载: 全尺寸图片 幻灯片

2.2   温度场的动态变化规律

高炉内温度场分布对高炉顺行及冶炼效率有重要影响,炉料的温度场决定了软熔带的位置,进而影响煤气流的分布。首先,对高炉炉料温度场进行数值解析,高炉在稳定状态下料批重发生改变后炉内炉料温度场的变化情况如图3所示,在ST1工况条件下,即当矿批重从207 kg变为163 kg,焦批重从64 kg变为50 kg时,炉内温度降低,尤其是靠近炉壁的区域更明显,低温区向下延伸。这是由于料批重减小时,上部料层变薄,料层数量增加,当炉料下降到软熔带附近时,铁矿石层软化熔化,较薄的焦炭层形成了较小面积的焦窗,增大了高温煤气流通过的阻力,因而软熔带上部温度降低。当料批重又调整回到矿批重为207 kg,焦批重163 kg时,料层温度开始回升,炉温恢复到先前的状态（即高炉处于稳定状态时的温度场分布）。在ST2工况条件下,矿批重从163 kg变为207 kg,焦批重从50 kg变为64 kg时,料层变厚,料层数减少,焦炭窗面积增加,有利于高温煤气流的通过,因此高炉上部干区温度升高。

图  3  ST1和ST2工况条件下炉料温度场分布

Figure  3.  Solid temperature distribution under ST1 and ST2 conditions

下载: 全尺寸图片 幻灯片

图4定量描述了炉料均温的瞬态变化情况。由图4（a）可见,第1次达到稳定状态时炉内炉料平均温度在860 ℃ 左右,当料批重发生改变,炉内温度平衡被打破,炉料均温随着料批重减小而逐渐下降,最终稳定在810 ℃左右,下降约5.81%；当料批重调整回到矿批重为207 kg,焦批重163 kg时,炉料平均温度逐渐恢复到860 ℃左右。由图4（b）可见,在达到稳定状态后,增加料批重,炉料均温缓慢上升。在此过程初始阶段炉料均温出现了小的波动,可能是由于料批重的改变引起软熔带焦炭窗形貌发生变化,影响了煤气流分布,进而影响了炉内温度分布。

图  4  ST1（a）和ST2（b）工况条件下炉料均温瞬态变化

Figure  4.  Instantaneous variation of the average burden temperature in the furnace under ST1 （a） and ST2 （b） conditions

下载: 全尺寸图片 幻灯片

当料批重改变时,原高炉内部稳定状态被打破,需要一定时间才能达到新工况条件下的稳定状态,研究该时间的变化情况可为优化调控高炉提供指导。ST1和ST2工况下各转变过程中的时间如图5所示,其中T₁、T₂、T₃分别表示第1、第2、第3次动态过渡状态。结果表明,在稳定状态料批重发生改变后,重新达到稳定状态所需时间逐渐减少。在ST1工况条件下,第1次过渡时间为805 min,第2次过渡时间为539 min,减少约33.04%；第3次过渡时间为410 min,较第2次减少约23.93%。在ST2工况条件下,第1次过渡时间为968 min,第2次过渡时间为678 min,减少约29.96%；第3次过渡时间为435 min,较第2次减少约35.84%。可能的原因是研究的实验高炉在经历多次的动态平衡后其自动调节能力变强,使其从一种稳定状态转变为另一种稳定状态的进程加快。另一个重要发现是ST1工况条件下在每个过渡状态持续的时间都少于ST2工况条件下相对应的各阶段,根据前面温度场的分析,ST1工况条件下炉料开始均温高于ST2工况,这导致ST1工况条件下可以更快地达到稳定状态。以此推测,炉内炉料均温越高,会在料批重改变时更快地向稳定状态转变。

图  5  过渡状态持续的时间

Figure  5.  Duration of transition state

下载: 全尺寸图片 幻灯片

2.3   软熔带的状态变化规律

软熔带是高炉内部非常重要的一个物理区域,其位置和形状能够直接反映炉内温度场分布及矿石颗粒相变情况,其一般温度范围为1 200~1 400 ℃。料批重是影响软熔带形成及演化的重要参数,图6（a）呈现了料批重改变后软熔带的形状和位置的变化情况。在减少料批重后,软熔带位置下降（ST1工况）,底部位置变化不大,上部位置降低明显；增加料批重后,软熔带高度增加,底部位置变化不大,上部位置明显上升（ST2工况）。较小料批重的料层下降,形成了较小面积的焦窗,阻碍高温煤气向上流动,导致上部干区温度降低,软熔带下降,反之亦然。

图  6  料批重改变后,ST1和ST2工况条件下软熔带形状（a）和高度（b）的演变

Figure  6.  Evolution of the shape （a） and height （b） of the soft melting strip under ST1 and ST2 working conditions after changing the batch weight

下载: 全尺寸图片 幻灯片

软熔带高度随料批重改变的动态变化情况如图6（b）所示。结果表明,在ST1料批重减小工况条件下,软熔带的高度由1.95 m降到1.81 m,降幅约7.18%；在ST2料批重增加工况条件下,软熔带高度由1.81 m增加到1.96 m,增幅约8.29%,两者相近。总的变化趋势反映出在料批重第2次调整后,软熔带的形状和位置能够逐渐恢复到最初料批重对应的炉况水平,这说明在实际工业生产过程中,当软熔带的形状和位置没有向着预期方向发展时,在有效的恢复时间内,可重新设置最初的料批重,理论上是可以恢复到原炉况。

料批重对炉内压力的影响如图7（a）所示,当料批重减小时,由于焦窗面积较小,导致煤气流通透性变差,软熔带下部压力增大,反之亦然。为定量化描述炉内部压降的瞬态变化情况,取2.00~6.50 m之间的压降进行分析（图7（b））,结果表明,在ST1工况条件下,压降由稳定时的6 100 Pa逐渐增加,到达新的稳定状态时,压降在7 000 Pa左右；ST2工况条件下,炉内压力变化呈现相反的趋势。综上,在ST1和ST2工况下,不论料批重增加或减少,在批重恢复后,压力分布也将恢复至先前水平。

图  7  料批重改变后,ST1和ST2工况条件下炉内压力分布（a）和压降（b）变化

Figure  7.  Evolution of pressure distribution （a） and pressure drop （b） in the furnace under ST1 and ST2 working conditions after changing the batch weight

下载: 全尺寸图片 幻灯片

2.4   矿石还原度的演变规律

图8呈现了料批重改变对炉内铁矿石还原度的动态影响情况。在本研究工况条件下,较大的料批重布料会呈现较高的铁矿石还原度分布,靠近炉内中心位置的铁矿石还原度较高（图8（a））,这与文献[15]一致。为更准确地定量展现铁矿石还原度的动态变化情况,图8（b）描述了炉内铁矿石平均还原度的瞬态变化分布。在ST1工况条件下,铁矿石平均还原度较高,稳定在0.65,当料批重减小时,高温煤气流在软熔带内通透性变差,影响炉内干区铁矿石的还原,铁矿石平均还原度减至0.63,降幅约3.08%,这表明能量消耗将随之增加；在ST2工况条件下,料批重增加时,高温煤气流在软熔带内的流通性得到改善,干区铁矿石平均还原度从0.63增至0.65。

图  8  料批重改变后,ST1和ST2工况条件下炉内铁矿石还原度分布（a）和平均还原度（b）的瞬态变化

Figure  8.  Transient changes of iron ore reduction degree distribution （a） and average reduction degree （b） in the furnace under ST1 and ST2 working conditions after the change of batch weight

下载: 全尺寸图片 幻灯片

3   结论

基于CFD-DEM同步耦合方法,从颗粒尺度解析了布料批重动态改变对高炉内部温度场、软熔带以及铁矿石还原度分布的影响,重点分析了料批重改变后主要参数变化以及转变过程中的特征,为高炉的实际生产操作以及向低碳操作炉况转变提供参考。主要得出以下结论：

1）料批重改变对炉料温度、软熔带形状和铁矿石还原度分布均具有显著影响。当矿批重从207 kg降至163 kg,焦批重从64 kg降至50 kg时,炉料温度降低约5.81%,软熔带高度减小约7.18%,炉内铁矿石平均还原度降低约3.08%,反之亦然。这表明在研究工况条件下,不论增加或减少料批重,当料批重恢复时,高炉主要运行参数及特征可以恢复到先前水平。

2）高炉内部从一个稳定状态调整料批重过渡到另一个稳定状态,两个稳定状态过渡过程的时间逐渐缩短。从本研究操作工况可知,ST1工况条件下,第2次过渡时间较第1次减少约33.04%,第3次过渡时间较第2次减少约23.93%；ST2工况条件下第2次过渡时间较第1次减少约29.96%,第3次过渡时间较第2次减少约35.84%。此外,炉内炉料均温越高,在料批重改变时,向新的稳定状态转变越快。