全氢竖炉内气流分布特性的数值模拟

余 珊， 孙 野， 陈 韧， 邵 磊， 邹宗树

（1. 辽宁科技学院 冶金工程学院， 辽宁 本溪 117004； 2. 辽宁省低品位非伴生铁矿优化应用重点实验室，辽宁 本溪 117004； 3. 东北大学 冶金学院， 沈阳 110819）

作为助力解决气候危机的重要手段，发展低碳经济、大幅降低碳排放已是全球共识，我国更是将“双碳”（即碳达峰、碳中和）大计提升到国家战略的高度.众所周知，钢铁行业是碳排放大户，因此大力推进钢铁行业低碳转型发展是保证“双碳”大计顺利实施和最终达成的关键.在此背景下，采用氢气替代碳素（焦炭、煤粉等）还原铁矿石的“氢冶金”成为现阶段人们普遍关注和热议的焦点，其中就包括100%“以氢代碳”的氢气竖炉工艺.

关于氢冶金，徐匡迪[1］早在2010 年就发起了加强相关基础理论研究和工程技术开发的前瞻性倡议.此外，考虑到当时我国的能源秉赋特点和钢铁行业状况，他还提出了以焦炉煤气作为原料量产制氢的建议.随着绿电（主要是太阳能、水能、风能等为代表的可再生能源以及核能发电）并网规模的不断扩大以及电解水制氢在技术和装备上的日臻成熟，近年瑞典率先提出和实施了涉及氢气竖炉和电炉的（近）零碳排放钢铁生产工艺研发项目，并将其命名为“突破性氢能炼铁技术”（HYBRIT）[2］.该项目已于2021 年6 月成功产出约百吨氢气直接还原铁（H2-DRI），展示了氢气竖炉炼铁在技术和装备上的可行性.

着眼于铁矿石的气体还原，若保持供气强度不变，仅将供气组成由传统“碳冶金”的CO＋H2改为全H2，带来的最大变化当属竖炉轴向温差增大，上部低温区扩展，还原反应的动力学条件恶化.这主要是因为采用H2还原铁氧化物在整体上是强吸热反应，与CO 刚好相反，如表1 所列.

表1 H2 还原铁氧化物的反应热[3］Table 1 Reaction enthalpy of iron oxide reduction with H2

因此，若想达到预定的生产要求，必须设法增加氢气竖炉的热收入，使加热固体炉料和进行还原反应的热需求得到满足后，仍能有效维持炉内温度.基于此，一些学者开展了相关理论研究，旨在明确供气温度、供气强度及球团入炉温度等决定系统热收入的因素对氢气竖炉生产效能的影响[4-7］.考虑到可直接使用电能供热，本文作者近期提出了在竖炉中上部增装微波设备的想法，设计了一款加热固体炉料的反应器，同时还进行了初步的模拟研究，结果与预期相符[8］.为进一步开展反应器尺度的研究，本文中采用数值模拟方法对氢气竖炉内的气流分布特性进行解析，具体考察竖炉底部中心辅助供气的影响，以期对氢气竖炉的设计和操作优化有所裨益.

1 模型构建

1.1 主要假设及简化

鉴于实际竖炉中气固两相间的热质传递过程极其复杂，为建立可用于解析炉内气流分布特性的数学模型，作出如下主要假设及简化：①假定竖炉连续稳定生产，模型中暂不考虑瞬态项的影响；②基于双流体假设构建模型，将气相作为流体处理，而将固相（颗粒群）视为与气相可相互穿透的拟流体；③模型中将气相（混合物）视为理想气体，其密度受温度和压力的影响，具体由气体状态方程计算，各物性参数通过理论或经验公式，由相关组分的物性参数计算得到；④由于竖炉内铁矿石的还原为气体的等摩尔反应，气相总流量在反应前后不发生改变，所以忽略化学反应对炉内气流整体分布的影响；⑤为节省计算资源，将竖炉内部空间简化为轴对称区域，模型仅考虑二维的计算域.

1.2 控制方程

为节省篇幅，以下仅给出模型所涉及的动量、热量、组分传输及湍流（k-ε 模型）相关控制方程的通式和主要符号的说明，各控制方程中更为细节的参数／系数及相间耦合关系的具体取值／表达式详见文献[9］和[10］.

式中：下标q 代表气相或固相；α 为相的体积分数；ρ 为相的密度，kg／m3；为速度矢量，m／s；ψ 为待求的独立变量.对于动量、热量、组分及湍流方程来说，分别对应速度、温度T、组分的质量分数Y，以及湍动能k 和耗散率ε；Γ 为广义扩散系数，根据独立变量ψ 的不同，可以取黏度系数、导热系数或扩散系数；S 为广义源项，可以为相间动量传输速率（采用矢量形式的Ergun 方程计算[11］）、热量传输速率（采用修正的Ranz-Marshall 方程计算[12］）或湍流相关速率参量.

1.3 计算方法及边界条件

本文中借助ANSYS-FLUENT 商业软件实现模型构建和计算，所采用的主要方法如下：①基于同位网格，采用SIMPLE 算法处理压力与速度的耦合；②各控制方程中的对流项和扩散项分别采用二阶迎风格式和中心差分（亦为二阶精度）格式离散；③为提高计算精度，各控制方程的收敛准则均为待求变量的残差且低于10-6.

根据竖炉自身的特点，计算域的顶部边界为固相（氧化球团）入口，但同时也是气相（尾气）出口，模型中固相视该处边界为速度入口，气相视其为压力出口，固相下行的速度根据球团的质量流量计算.计算域最底端边界为固相速度出口，气相质量流量设定为0.计算域中气相可通过侧壁围管和底部中心喷入，对应边界处气相的质量流量由预先给定的供气强度计算得到，但固相的质量流量设定为0.当通过相关边界进入计算域时，气固两相的温度均预先给定.另外，计算域中所有壁面均为无滑移绝热边界.

1.4 模型验证

鉴于目前还未出现关于氢气竖炉生产数据的公开报道，本文中将模型应用于八钢C-3000 煤气竖炉，通过对比竖炉围管气压的测量值和计算值进行模型精确度的验证.在验证性算例中，竖炉的顶压为300 kPa，煤气全部由竖炉围管喷入，入炉煤气的组成、温度及供气强度分别为65%CO ＋23%H2＋9%CO2＋3%H2O（体积分数），1 100 K，1 050 m3／t（标准态下氧化球团），氧化球团的质量流量为218 t／h，其他几何及操作参数详见文献[10］和[13］.

上述工况下C-3000 竖炉围管气压的测量值约为360 kPa，本文模型的计算值为357.2 kPa，两者相差很小，表明本文中的数学模型具有较高的精确度.

2 结果分析及讨论

2.1 计算方案

经上述验证后，将供气组分限定为H2，便可形成氢气竖炉模型.图1 示出了本文中考虑的二维轴对称计算域及相应网格的配置（网格数量和尺寸比）.

图1 本文考虑的氢气竖炉计算域（m）Fig.1 Computational domain of the H2 shaft furnace considered in the current paper（m）

从图1 中可以看到，供气可由侧壁围管和底部中心喷入炉内，由于存在很大的速度梯度，侧壁围管和底部中心处的网格均进行了局部加密，尺寸为其他未加密部位网格的1／100 左右.计算区域共包含5 000 个四边形网格，局部加密网格的尺寸约为5.4×10-5m2，其他未加密部位网格的尺寸约为5.3×10-3m2.需要说明的是，本文中所采用的网格配置已根据前期参数（采用围管平面的平均气压）的敏感度分析进行了优化.具体优化过程如下：在保持上述网格尺寸比（即1／100）不变的前提下，分别采用网格数量为3 000，4 000，5 000和6 000 的四种配置进行仿真模拟，计算所得围管平面的平均气压分别为315.4，321.1，323.4 和323.7 kPa.

在中心供气方面，相关的工程实践可参考八钢C-3000 的中心气体分配装置（CGD）.氧化球团常温入炉，等效直径、质量流量及床层空隙度分别为15 mm，100 t／h，0.4.气相入炉温度、总流量（体积）及炉顶压力分别为1 173 K，1 700 m3／t，300 kPa.本文中将竖炉中心供气比记作β，定义为通过底部中心供入气相的流量占总流量的比例.β作为被考察因素，变化范围为0 ～0.20.另外，考虑气相需增压才能被顺利喷入炉内，采用等熵公式计算压缩气相对应的（单位质量球团）能耗，具体表达式如下：

式中:n为气相的流量， mol／h； m 为固相质量流量，t／h；R 为理想气体常数，8.314 J／（mol·K）；P为目标气压，Pa；λ，T0，P0分别为等熵压缩效率、参考温度及参考气压，在以下算例中分别取0.7，298 K，100 kPa.

2.2 气流分布特性

图2 （a）（b）分别为无中心供气（β ＝0）条件下，氢气竖炉模型内的气相流线和压强分布情况.氢气由围管喷入后，在沿径向流动的同时会将初始动能迅速转化为轴向速度，从而逐渐转向炉顶出口，最终形成如图2 （a）所示的L 型流线，这与文献报道的COREX 预还原竖炉内（无中心供气）的气相流线类似[14］.由图2 （b）可知，除近围管的有限区域外，竖炉内气相压强沿径向整体分布较均匀. 轴向上围管平面的平均气压约为323.4 kPa，而竖炉顶压为300 kPa.也就是说，本文中所考虑的工况下围管平面以上的固相床层内压降较小，约为23.4 kPa.这主要是因为气体通过固相床层流动的惯性阻力与其密度成正比，而同等条件下（温度、压力、摩尔量）H2的密度仅为CO 的1／14.

图2 无中心供气（β ＝0）时氢气竖炉内的计算结果Fig.2 Calculated results ofthe H2 shaft furnace without central gas supply（β ＝0）

为获得更多流场信息，截取图2（a）中6 个轴向深度（即z ＝6.3，6.6，6.9，7.2，7.5，7.8 m）的气相（总）速度进行分析，结果如图3 所示.

从图3 中可看出，各轴向深度的气相速度虽然在临近围管的侧壁处较大，但都沿内径向衰减，在距竖炉中心约1.0 m 处基本达到稳定.需要注意的是，竖炉中心处（r ＝0）的气相速度沿轴向深度逐渐减小.据此可以推断，在竖炉底部炉芯部位存在一个气流迟滞区.为活跃该区域的气相流动，本研究中通过实施中心辅助供气，具体考察中心供气比对竖炉中心（r ＝0）气相速度（轴向）分布、炉内气相流线及压力分布的影响.

图3 无中心供气（β ＝0）氢气竖炉内不同轴向深度的气相速度径向分布Fig.3 Radial distributions of gas phase velocity at different axial depths in the H2 shaft furnace without central gas supply（β ＝0）

2.3 中心供气的影响

图4 为不同中心供气比下氢气竖炉内气相速度在轴向上（z ＝6.3 ～7.8 m）的分布图.由图可知，实施中心辅助供气可有效提高竖炉底部的气相速度.气相速度沿轴向深度先缓慢降低，后又迅速升高.当β 由0.05 提高至0.20 时，竖炉中心处气相速度的最小值从6.3 m／s 增大至9.3 m／s，增幅为47.6%.

图4 不同中心供气比下的氢气竖炉内气相速度轴向分布Fig.4 Axial distributions of gas phase velocity in the H2 shaft furnace at different proportions of central gas supply

图5 为不同中心供气比下氢气竖炉内的气相流线图.从图中可以看到，实施辅助供气可在竖炉中心形成一股集中气流.随着β 的增大，该股气流的初始动能也增大，其径向穿透距离呈现一定程度的延伸.

图5 不同中心供气比下的氢气竖炉内气相流线Fig.5 Gas phase streamlines in the H2 shaft furnace at different proportions of central gas supply

图6 为不同中心供气比下氢气竖炉内的气相压强分布图.由图可知，β 的增大对竖炉内围管平面以上气压分布的影响不明显，但是底部供气边界处的大压力梯度区却逐渐扩大，这就意味着压缩气体所需能耗在增加.

图6 不同供气比下的氢气竖炉内气相流线Fig.6 Gas phase pressure distributions in the H2 shaft furnace at different proportions of central gas supply

表2 列出了不同中心供气比下底部和围管入口的气相流量和从模型计算结果中提取的气相压强.基于表中数据，采用式（2）可计算出压缩气体能耗，结果如图7 所示.从图中可以看到，随着β的增大，底部供气能耗因气体流量的提高而逐渐增大，围管供气的能耗随之减小.但必须指出的是，压缩气体所需的总能耗随β 的增大先降低后升高，在β ＝0.10 时达到最低，约为345 MJ／t.

表2 不同中心供气比下底部和围管入口的气相流量和压强Table 2 Flow rates and pressures at the gas inlets under different proportions of central gas supply

图7 不同中心供气比下的压缩气体所需能耗Fig.7 Energy requirements for gas compression at different proportions of central gas supply

3 结 论

（1）实施中心辅助供气可有效提高氢气竖炉底部的气相速度.当中心供气比由0.05 提高至0.20 时，竖炉中心处气相速度的最小值从6.3 m／s增大至9.3 m／s，增幅为47.6%.

（2）在文中所考虑的变化范围内，中心供气比的提高对竖炉内围管平面以上气压分布的影响不明显.

（3）对竖炉生产而言，压缩气体所需的总能耗随中心供气比的增大先降低后升高，在中心供气比为0.10 时达到最低，约为345 MJ／t.