水泥回转窑传热特性数值模拟

2019-07-08周武洲赵朋程

燕山大学学报 2019年3期

刘彬，周武洲，赵朋程，李瑞

(1.燕山大学信息科学与工程学院，河北秦皇岛 066004；2.燕山大学电气工程学院，河北秦皇岛 066004)

0 引言

随着近些年水泥行业的高速发展，熟料烧结过程中的高能耗和高排放的问题，引起越来越多的关注[1-2]。每生产一吨水泥，一般要消耗3 GJ能量，同时排放出比所生产水泥总量更多的温室气体[3]。研究回转窑工况对窑内流场、温度分布以及煤粉燃烧的影响，可以为水泥厂的节能减排、改进水泥生产工艺提供理论指导，同时也可以为燃烧器的优化提供理论依据。然而由于回转窑内传热过程复杂、测量手段不完善等原因[4]，无法准确得到回转窑内部的流场、温度分布信息，故而通过数值模拟的方法分析回转窑内的流场和温度场的分布情况具有很大研究价值。

针对回转窑内部环境的数值模拟，国内外已经做过很多研究，如Elattar等[5]建立回转窑的二维模型，研究了回转窑主要操作对不同燃料下回转窑内火焰行为的影响；Mujumdar等[6]提出将熟料看作一种伪流体，分别对回转窑内的空气和熟料进行三维建模，以伪均匀近似法来模拟固体-固态反应，建立了水泥回转窑内的综合CFD模型；Manju等[7]利用两相流方法对回转窑气动煤注入和燃烧过程进行了建模，研究了回转窑操作参数对注入煤粉颗粒的分布及其提供热负荷的影响；Kangwanpongpan等[8]通过优化灰度气体模型的加权和，并结合离散坐标法(Discrete ordinate，DO)模拟了富氧燃烧环境下的辐射热传递，气体火焰温度的预测效果得到了改进；Gaikwad[9]等建立了简化CFD模型，该模型采用DO模型与加权和灰度气体方法耦合的辐射热传递模型，以该模型研究了空气、氧蒸汽等燃烧环境对温度和NO浓度分布的影响；廖斌[10]等针对一种低NOX四通道燃烧器，利用FLUENT软件研究了影响回转窑内NOX产生的主要因素以及降低NOX排放的方法；马爱纯等[11]利用商业CFD软件，对装备四通道旋流燃烧器的回转窑进行了数值模拟，分析了不同工况下窑内流场分布、温度场分布、NOX生成和火焰的特点。Liedmann等[12]建立了固定斜度的3D回转窑模型，利用商业CFD软件研究了煤粉掺杂垃圾衍生燃料(refuse derived fuels，RDF)后，燃烧器操作参数对回转窑内温度场的影响。

虽然上述学者所建立的模型能较好地模拟窑内燃料燃烧的特性，然而由于未能综合考虑回转窑斜度的影响，所以窑内温度场与流场的数值模拟结果会和实际情况有较大的差别。

本文在以上研究的基础上，建立回转窑内湍流流动、辐射换热和煤粉燃烧的数学模型。为了研究回转窑斜度对窑内流场以及煤粉燃烧等的影响，保持燃烧器水平方向不变，对所建立的4个不同斜度的回转窑-燃烧器流场三维几何模型进行全结构化网格划分，采用数值模拟的方法，利用FLUENT软件研究斜度、旋流风速和轴流风速对窑内流场、窑内空气温度场与熟料表面空气温度场的影响。

1 数学模型

回转窑内的湍流流动、煤粉燃烧和辐射换热数值模拟，涉及湍流模型、煤粉燃烧模型以及辐射模型等一系列物理化学过程子模型。为了使仿真结果中回转窑内部的流场和温度场更符合实际情况，需要根据回转窑内的流场特点以及四通道旋流燃烧器的燃烧特点进行子模型的选择，下面对各子模型进行分析。

1.1 湍流模型

在文献[11]中发现，标准k-ε模型的计算结果与回转窑实验测试结果更吻合，因此气相湍流模型选择标准k-ε模型。标准k-ε模型是基于湍流动能k及其耗散率ε的模型传输方程的模型。k的模型传输方程是从精确方程得出的，而ε的模型传输方程是使用物理推理获得的，与其数学上的精确方程非常相似，标准k-ε模型的传输方程为

Gk+Gb-ρε-YM+SK,

(1)

(2)

(3)

式中，μt为湍流(或涡流)粘度，Gk表示平均速度梯度产生的湍流动能，Gb表示由浮力的影响产生的浮力动能，YM表示在可压缩湍流中的波动膨胀对总体耗散率的贡献，常量C1ε=1.44、C2ε=1.92、σk=1、σε=1.3。

1.2 煤粉燃烧模型

煤粉燃烧过程大致可以分为热解和煤焦燃尽两个阶段。对于热解反应，本文采用目前工程模拟上广泛使用的双竞争模型[13]，该模型将煤粉的热解看成两个热解反应，这两个热解反应同时进行但反应速率不同，表示如下：

(4)

(5)

式中，k1，k2为反应速率常数，服从Arrhenius公式，α1表示煤粉的低温析出性能，α2为煤粉的高温析出性能。在低温段k1>k2，而高温段k1

煤焦的燃尽过程比煤粉的热解反应要慢得多，因此煤粉的燃烧时间主要受煤焦燃尽的影响，现在主流的研究观点认为煤焦的燃烧速率同时受化学反应速率和氧气扩散速率影响[14]，煤焦燃烧速率为

(6)

扩散反应速率常数为

(7)

化学反应速率为

R=C2e-(E/RTP)，

(8)

式中，A为煤粉颗粒表面积，Pox为煤粉颗粒周围氧气分压，Do为扩散系数，R为化学反应速率常数。

1.3 DO辐射模型

DO辐射模型通过有限数量的离散立体角建立辐射传递方程(Radiation transfer equation，RTE)，每个离散立体角与全局笛卡尔系统中固定的向量方向相关联。因为DO模型是通过使用像素化实现的[15]，故DO模型最大的特点是空间方向的离散化和相点的处理。在DO模型中，任何位置处的4π角空间的象限被离散为Nθ×Nη个立体角ψ，称其为控制角，其中角度θ和η分别是极坐标和方位角。

DO模型将RTE视为场方程，因此RTE方程可以写成

▽·(I(r,s)s)+(α+σs)I(r,s)=

(9)

该模型可以使用灰度模型对非灰度辐射进行建模。光谱强度Iλ(r,s)[27]的RTE为

▽·(Iλ(r,s)s)+(αλ+σs)Iλ(r,s)=αλn2Ibλ+

(10)

式中，λ为波长，αλ是光谱吸收系数，Ibλ是由普朗克函数给出的黑体强度，并假设散射系数，散射相位函数和折射率n与波长无关。

2 几何模型与边界条件

2.1 回转窑-燃烧器几何模型的网格划分

为了模拟回转窑实际生产环境，本文回转窑几何模型采用Ф4.0 m×70 m尺寸，但考虑到回转窑内耐火砖的厚度，实际回转窑内径设置为3.6 m。

本仿真所采用的燃烧器为目前广泛使用的四通道旋流燃烧器，在按实际燃烧器和回转窑1∶1比例建立几何模型后，抽取其内部流体域模型，燃烧器内部流体域结构剖面图如图1所示。燃烧器4个风道由外向内依次为轴流分道、煤粉风道、旋流风道和中心风道。

熟料平面的位置可以由动态休止角β和物料中心角ω确定，如图2所示。物料中心角可以根据填充率φ计算：

2πφ=ω-sinω。

(11)

设置计算工况物料填充率为10%，回转窑转速为2 r/min，计算可得ω=93.16°。根据文献[16]，水泥物料的动态休止角β在25°～33°之间，取β=27°。根据上述计算，建立去除熟料的回转窑内部流体域模型。

图1 燃烧器内部流体域几何模型结构图
Fig.1 The structure of the geometric model of the fluid domain inside the burner

图2 回转窑剖面熟料区域示意图
Fig.2 A schematic diagram of the clinker area on the rotary kiln profile

实际生产中，回转窑斜度i一般控制在0.03～0.04之间，故本文建立的4个回转窑模型的斜度i分别取0、0.03、0.035、0.04。

本文采用工程上广泛使用的网格划分软件ICEM对燃烧器内部流场进行全结构化网格划分。与非结构化网格相比，结构化网格大大减小了网格数量，因此极大降低了仿真时的计算量，同时增加了计算结果的准确性并提高了计算的收敛性。此外，结构化网格的计算效率也更高，因此，本文网格划分采用全结构化网格划分方法。燃烧器网格如图3所示。

回转窑及燃烧器入口网格轮廓如图4所示。由于回转窑长度过长，故以斜度为0.03的回转窑为例，将燃烧器和该回转窑组成的整体网格的剖面图局部显示如图5所示。

图3 燃烧器网格
Fig.3 The mesh of burner

图4 回转窑及燃烧器入口网格轮廓
Fig.4 Rotary kiln and burner entrance mesh profile

图5 斜度为0.03的回转窑网格局部剖面图
Fig.5 Partial profile of rotary kiln mesh with gradient of 0.03

2.2 边界条件

本仿真参考实际回转窑生产工况进行仿真边界条件设置，设置的标准工况边界条件参数如表1所示。实际中煤粉颗粒直径大小为离散分布，假设煤粉颗粒直径服从Rosin-Rammler分布。为了使仿真更符合实际情况，将煤粉颗粒设置为9种不同直径的单元，按照工况煤流量总量计算各个煤粉单元所占煤流量大小，煤粉流量分布表格如图6所示，总煤流量设置为2.46 kg/s，回转窑标准工况斜度设置为0.03, 氧含量设置为0.2315。燃煤的元素分析和工业分析如表2所示。回转窑壁面设置为绝热无滑移边界。

表1 标准工况边界条件参数表Tab.1 Boundary condition parameter table under standard operating conditions

图6 不同煤粉颗粒直径流量分布
Fig.6 Flow distribution of different pulverized coal particle diameter

表2 燃煤的元素分析和工业元素分析Tab.2 Ultimate and proximate analyses of the coal burned

为了确保所划分的网格能够满足仿真计算的需求，需要进行网格无关性验证，以确保仿真结果不会因网格数量过少而影响计算精度。边界条件均按表1进行设置，划分3组不同密度的网格模型，对3组网格模型分别进行仿真计算，计算结果如表3所示。

根据表3，综合考虑计算精度和计算量的要求，选择网格数量为453 114的网格步长进行不同斜度模型的网格划分。

表3 网格无关性验证结果Tab.3 Simulation result of grid independent verification

3 仿真结果与分析

为了研究不同回转窑斜度以及不同回转窑工况对回转窑内流场、煤粉燃烧等的影响，本文在2.2节边界条件的基础上设置两组实验：斜度对流场和煤粉燃烧的影响实验，将回转窑运行参数按表2设置，用不同斜度的回转窑模型进行仿真模拟计算；燃烧器工况对流场和煤粉燃烧的影响实验，固定回转窑斜度为0.03，分别改变旋流风道风速和轴流风道风速，其他参数按表2设定，进行仿真模拟，将仿真结果与标准工况的进行对比分析。

基于第一部分建立的数学模型，利用FLUENT软件对每组实验进行数值模拟。先开启1.1节描述的标准k-ε模型，进行450次冷态无反应流仿真迭代，以使回转窑内流场趋于稳定，然后开启1.2节的煤燃烧模型来点火，进行1 350次热态反应流仿真迭代，再开启1.3节中的DO辐射模型，进行300次带辐射的反应流仿真，最后开启粒子辐射相互作用，进行300次带粒子辐射相互作用的反应流仿真。

3.1 斜度对流场和煤粉燃烧的影响

根据上节边界条件进行仿真实验，由于回转窑长度过长，这里仅显示窑头15 m范围温度场的仿真结果，窑头部分垂直方向轴线位置处的温度场仿真结果如图7所示。

图7 不同斜度下窑头部分垂直方向温度分布图
Fig.7 The temperature distribution in the vertical direction of the kiln portion at different slopes

从图7可以看出，相较于i=0时回转窑内的火焰，i≠0时回转窑内火焰存在黑火头，垂直方向上火焰形状不规则情况明显，然而在i≠0的3组实验结果中，i=0.035时的火焰形状明显比i=0.03及i=0.04的更细长且更规则。

图8给出了i=0时轴线方向YZ平面上的流场速度分布，由图8可以看出，中心风道附近存在回流区，该回流区是由旋流风的强旋流效应和旋流风与中心风之间的速度差共同引起的。而图8中的回流区方向明显偏下，这是由于熟料区的存在使回转窑内的流场在靠近熟料区附近时风速更高，进而吸引中心回流区向下偏斜，而斜度的存在会加剧这种偏斜，造成回流区的减弱甚至消失。

图9给出了i=0.035时的流场速度分布，在图9中，可以发现中心回流区消失，这正是回转窑内下方流场流速过高，与上方流场形成较大速差射流引起的。

中心回流区可以卷吸下游高温烟气，有助于煤粉的点火，中心回流区的减弱或消失使得煤粉的着火点后移，形成黑火头。黑火头的存在可以防止燃烧器温度过高，延长燃烧器寿命，但过长的黑火头会使煤粉燃烧效率下降、煤耗增加等不利于熟料烧结的因素产生。

图8i=0时轴线位置垂直方向流场速度分布
Fig.8 The flow velocity distribution at the axis position in the vertical direction wheni=0

图9i=0.035时轴线位置垂直方向流场速度分布
Fig.9 The flow velocity distribution at the axis position in the vertical direction wheni=0.035

图10给出了不同斜度下回转窑轴线上温度分布，由此可以看出，i=0.035时回转窑内的温度分布要明显高于i=0.03和i=0.04的。图11给出了不同斜度下熟料平面轴线上温度分布，由图11也可以看出，i=0.035时，熟料平面上温度高于1 600 K烧成带范围要明显大于其余两种斜度的。所以i=0.035更适合熟料的烧结。

图10 不同斜度下回转窑轴线上温度分布
Fig.10 Temperature distribution on the rotary kiln axis at different slopes

图11给出了不同斜度下熟料平面轴线上温度分布。

由图11可以发现i=0.03和i=0.04的温度曲线几乎重合，但图10中i=0.04的轴线温度曲线峰值温度明显滞后于i=0.03的，由此可以看出，提高回转窑斜度会使烧成带区域靠近窑头，进而使冷却区范围缩小。图11中回转窑出口位置熟料平面温度升高，这是由出口处高温烟气回流造成的。

图11 不同斜度下熟料平面轴线上温度分布
Fig.11 Temperature distribution on the plane axis of clinker at different slopes

3.2 燃烧器工况对流场和煤粉燃烧的影响

3.2.1旋流风速的影响

设置旋流风速分别为70 m/s、80 m/s、100 m/s和110 m/s，其余边界条件按标准工况进行设置，对比标准工况下的仿真结果。不同旋流风速下的温度分布对比结果如图12和图13所示。

图12 不同旋流风道风速下回转窑轴线温度分布图
Fig.12 Temperature distribution on rotary kiln axis at different swirling wind speeds

由图12可以看出，s=80 m/s时的回转窑轴线峰值温度最高，但是在图13中，s=80 m/s时的熟料温度分布整体并非最高的。造成这种现象的原因主要是：旋流风速的降低使入口切向速度降低，入口切向速度的降低有利于煤粉的聚集，这使得s=80 m/s时火焰有较高的温度，同时过低的旋流风速大大降低了窑内气相流场的旋流强度，进而减弱了旋流风对火焰的约束，降低了火焰的稳定性，使得火焰热量不能有效地用于熟料的烧结。图14为不同旋流风速下窑头部分水平方向温度分布图，图14(a)和(b)的火焰形状明显要比(c)和(d)的更不规则，也说明了减弱旋流风降低了火焰的稳定性。

图13 不同旋流风速下熟料平面轴线上温度分布
Fig.13 Temperature distribution on clinker plane at different swirling wind speeds

图14 不同旋流风速下窑头部分水平方向温度分布图
Fig.14 Temperature distribution in horizontal direction at different swirling wind speeds

在图12中，s=100 m/s时温度在初始阶段的升高速度明显要低于s=80 m/s和s=90 m/s的，而图14(e)中火焰的黑火头明显要比其余状态下的长，这是由于旋流风速的提高增加了入口流场的轴向速度和切向速度，切向速度的增加使煤粉可以与空气充分混合，同时增强对下游高温烟气的卷吸，但是由于斜度的存在使中心回流区消失，旋流风对下游高温烟气的卷吸能力被削弱，而入口流场轴向速度的增加结合旋流风卷吸高温烟气能力的减弱共同导致了着火点的后移。

在实际回转窑运行中，旋流风速的增加，会导致在燃烧器出口区域的煤粉颗粒扩散，这会使高旋流风速的火焰变粗，而低旋流风速的火焰更细长，对比图14 中各图可以发现，图14(c)中的火焰明显要更细长，这与实际运行情况相符，而图14(a)和(b)中的火焰由于旋流风速过低，火焰稳定性不足，造成火焰形状不规则。对比s=90 m/s和s=100 m/s在图12和图13中的温度曲线，可以发现二者十分相似，旋流风速的提高并没有明显提升回转窑内温度分布水平。结合以上分析，可以认为90 m/s是比较好的旋流风速选择。

3.2.2轴流风速的影响

分别设置回转窑轴流风速a=150 m/s、a=160 m/s、a=170 m/s和a=200 m/s，其余边界条件按标准工况进行设置，对比标准工况下a=180 m/s的仿真结果。轴流风速为150、160、170、180和200 m/s的温度分布对比结果如图15和图16所示。

由图15可以看出，随着轴流风速的提高，回转窑内火焰峰值温度明显下降，且峰值温度随轴流风速的提高而向窑尾偏移。轴流风速a由150 m/s增加到180 m/s，回转窑火焰峰值温度降低了近150 K。在实际回转窑运行过程中，过高的轴流风速会形成大速差射流，这种强烈的速差射流会在燃烧器轴线位置形成较大的负压，进而阻碍中心回流区的形成，使中心回流区变小,甚至消失。同时，提高轴流风速虽然可以卷吸更多的高温二次风，但是风速的增加也会使这些高温二次风与煤粉的混合位置后移，使煤粉的预热延迟。中心回流区的减小以及煤粉预热的延迟共同引起了回转窑内火焰峰值温度的后移。并且提高轴流风速，会加速煤粉的运动速度，这将分散煤粉的燃烧，使回转窑火焰温度降低。

图15 不同轴流风速下回转窑轴线温度分布图
Fig.15 Temperature distribution of rotary kiln under different axial flow rates

图16 不同轴流风速下熟料平面轴线上温度分布
Fig.16 Temperature distribution of clinker plane axis under different axial flow rates

图17为不同轴流风速下窑头温度轮廓图，对比图17(a)～(e)，可以看出，当轴流风速增加到200 m/s时，回转窑内的火焰严重变形，这是由于强烈的轴流风速破坏了旋流风对煤粉的卷吸，使煤粉严重扩散，难以聚集燃烧引起的。而对比图17(a)～(d)可以发现，随着轴流风速的增加，黑火头长度逐渐缩短，这也验证了提高轴流风会卷吸更多的高温二次风对煤粉进行预热的理论。