适应锅炉调峰运行的水冷壁高温腐蚀预测模型

邓磊，袁茂博，杨家辉，岳洋，姜家豪，车得福

（西安交通大学能源与动力工程学院，陕西 西安 710049）

2022 年我国风能、水能、太阳能等可再生能源发电占总发电量的28.7%，该比例在“双碳”目标的驱动下仍会增长[1-2]。但由于风、光、水资源的时空分布不平衡，传统火电站在调配能源中的地位将显著增强。同时，大型电站煤粉锅炉普遍采用分级配风降低NOx排放[3]，这也导致主燃区形成了强还原性气氛，加剧了水冷壁高温腐蚀，直接影响到锅炉的安全运行[4]。Sun等[5]和Xiong等[6]均在实施空气分级改造后的锅炉中发现了腐蚀深度超过1.5mm 的水冷壁管。于英利等[7]通过调节双切圆锅炉同层二次风以缓解水冷壁高温腐蚀问题。王新宇等[8]以对冲燃烧锅炉为研究对象，模拟水冷壁近壁面的H2S 和CO 分布并提出高温腐蚀防治措施。然而前人研究主要为特定负荷下的炉膛烟气场模拟和腐蚀现象分析，难以直接应用于锅炉调峰运行水冷壁高温腐蚀程度的评估。Kung[9]通过实验得出了锅炉钢年腐蚀速率与壁温、H2S浓度以及合金中铬元素含量的拟合式，如式(1)所示。

式中，CR 为年腐蚀深度，mm；T为金属温度，K；CH2S为H2S 浓度，μL/L；ωCr为合金中铬元素的质量分数，%。

式(1)的年腐蚀速率是以实验腐蚀速率线性外推至全年得到的，但Xu等[10]的相关实验表明，H2S高温腐蚀速率与腐蚀时间是非线性。同时，水冷壁的年腐蚀深度是由不同工况下的腐蚀深度叠加组成，而各个工况的运行时长也存在差异，因此时间维度在锅炉调峰运行下水冷壁高温腐蚀深度的计算中尤为重要。对此，Yuan 等[11]提出了包含时间维度的H2S高温腐蚀数学模型，如式(2)所示。

式中，d-为腐蚀深度，μm；t为腐蚀时间，h；E为腐蚀活化能，kJ/mol；R为气体常数；系数a、b和指数n取决于金属性质。

式(2)将腐蚀深度与壁面温度、H2S 浓度以及腐蚀时间联系起来，根据水冷壁壁温分布、近壁面H2S 浓度分布和运行时间即可得到水冷壁高温腐蚀状态。在前期工作[12-13]的基础上，本文将水冷壁热流分布模拟计算与水动力特性计算相耦合，进一步结合管壁温计算，提出水冷壁壁温分布计算方法。同时，根据燃料硫释放以及含硫组分的相互转化特性，建立一种适用于还原性气氛的SOx生成模型，可计算炉膛内H2S 浓度分布。综合炉膛数值模拟、水冷壁壁温耦合计算以及高温腐蚀数学模型，构建出适应锅炉调峰运行的水冷壁高温腐蚀预测模型并基于Matlab GUI 平台开发对应预测软件。选取一台600MW 四角切圆燃煤锅炉为研究对象并建立炉膛数值模型，对100%锅炉最大连续蒸发量工况（BMCR）、75%热耗率验收工况（THA）、50%THA 以及35%BMCR 四种典型负荷开展数值模拟计算，获得壁面热流密度分布以及近壁面H2S 浓度分布；将数值模拟结果导入高温腐蚀预测模型可得到单一负荷下或多个负荷下的水冷壁高温腐蚀状态，进而为锅炉调峰运行水冷壁安全监测提供参考。

1 锅炉数值计算模型

1.1 物理模型与工况设置

本文选取一台超临界600MW 直流四角切圆燃煤锅炉为研究对象。由于实际锅炉结构复杂，如图1 所示对其几何形状和燃烧系统进行了合理简化。炉膛截面深度为17.696m、宽度为18.816m、锅炉总高为63.750m，水冷壁由倾角13.95°的螺旋管圈和垂直管圈组合构成。燃烧系统由煤粉喷嘴（PA）、 辅 助 风 喷 嘴（AUX）、 紧 凑 燃 尽 风（CCOFA）以及可分离燃尽风（SOFA）组成，辅助风喷嘴位于相邻的煤粉喷嘴之间，并由上下两只预置水平偏角的辅助风（CFS）喷嘴以及一只直吹风喷嘴组成。燃烧器中心线和炉壁的夹角为51°和48°，旨在炉膛内形成切圆燃烧，CFS 喷嘴和炉壁夹角为22°和19°。燃烧系统将下炉膛分为主燃区、还原区和燃尽区三个部分。炉膛上部从前到后依次布置有分隔屏、后屏、末级再热器和末级过热器，均简化为无厚度平面。100%BMCR、75%THA、50%THA 和35%BMCR 四个负荷的具体参数设置如表1所示。

表1 不同负荷下锅炉运行参数

图1 锅炉布置简图

1.2 数值模型选择

参照图1所示锅炉布置简图建立炉膛的三维模型，应用Fluent软件设置计算模型，对炉膛传热传质过程做稳态计算。选取Realizablek-ε湍流模型和DO辐射模型分别对流场和辐射传热进行计算[14]；采用灰色气体加权求和模型计算烟气的辐射吸收系数；采用基于欧拉-拉格朗日方法的随机轨道模型模拟煤粉颗粒的运动轨迹，粒径大小遵循Rosin-Rammler 分 布[15]，范 围 为1～100μm，平 均 粒 径61μm，煤样的元素分析和工业分析见表2。采用Multiple surface reaction 模型计算煤焦与O2/CO2/H2O三种气体的异相反应速率[16]。炉膛传热传质计算收敛后，加载UDF自定义SOx生成模型，进一步计算H2S、COS 和SO2的浓度场。SOx生成模型定义燃料硫向H2S、COS以及SO2转化的速率由式(3)～式(5)表示；含硫组分的相互转化为表3所示的10个总包反应[17]，最终三种气体的生成速率（rf）与消耗速率（rr）由式(6)～式(11)表示。

表2 煤样的工业分析与元素分析

表3 含硫组分相互转化反应方程[17]

式中，α和γ分别为挥发分和灰分中硫的质量分数，分别取值0.55 和0.1；rvol和rc分别为挥发分的挥发速率和焦炭的燃烧速率，kg/s；Mw,S为硫的摩尔质量，0.032kg/mol；Vcell为网格单元体积，m3；pc为焦炭氧化反应速率占焦炭总反应速率的比例，计算方法参见本文作者课题组[13]前期的工作。

1.3 数值模型验证

为选择合适的网格密度，对比了网格数量分别为1.34×106、1.75×106、2.07×106和2.48×106的四个模型沿炉高方向截面平均烟温，结果如图2所示。网格数量为207 万和248 万的温度曲线十分接近，继续增加网格数量对计算结果准确度的提高有限。综合准确性和经济性两方面因素，选用207 万网格数量的炉膛模型开展研究；选取100%BMCR 负荷下实炉测量数据与模拟结果进行对比，验证数值计算的可靠性。实测数据的获取方式为：水冷壁吸热量是根据省煤器出口和汽水分离器入口的工质焓差计算得到；炉膛出口氧浓度和后屏出口烟温是根据布置在烟道中的若干测点数据取平均得到；H2S浓度值是在水冷壁近壁面采用烟气分析仪测量得到。如表4所示，水冷壁吸热量、炉膛出口氧浓度和后屏出口烟温的相对偏差分别为7.50%、2.76%和4.79%。锅炉运行过程中壁面附近的H2S浓度会在一定范围内波动，因此数值模拟计算的H2S浓度与实际测量值之间存在一定的偏差。标高25m、29m 和35m 分别对应于A 层、C 层和F 层燃烧器高度位置，前墙中心线上H2S 浓度的相对偏差分别为7.76%、8.21% 和8.37%，可见各维度的计算偏差均在合理范围内。

表4 100% BMCR负荷下实炉测量结果与模拟结果对比

图2 网格无关性验证

2 水冷壁壁温耦合计算方法

螺旋管圈水冷壁为大容量电站锅炉下炉膛水冷壁的主要形式，传统的水动力回路划分是对炉膛四面墙进行分割，没有考虑螺旋管圈水冷壁的特殊布置形式，而工质在螺旋管圈水冷壁内盘旋上升的行程是连续的。本节提出了一种基于坐标变换的热量再分配方法，建立了水动力回路吸热量与工质在流动过程吸热量的映射关系。进一步耦合烟气侧与工质侧传热，结合管壁温计算，形成水冷壁壁温耦合计算方法。算法在Matlab平台上实现，具有集成度高、灵活性强的特点。

2.1 螺旋管圈水冷壁热量再分配方法

图3 为螺旋管圈水冷壁展开与重构的示意图。首先将水冷壁展开形成由前墙、右墙、后墙、左墙依次排列的平面，原三维螺旋管圈中某连续回路（ABCDEF）在展开平面上被分为ABCDE和EF的两段。为使回路重新拼合，对AE连线以上区域整体向右平移，直至E点与（E）点重合，移动距离为四墙宽度的总和。由于螺旋管圈普遍匝数在1.5 左右，因此E点落在左墙的高上，此时还需将（E）F连线以上区域向右平移两个四墙宽度总和。平移过程同样是壁面热流密度点坐标变换的过程。对变换后的平面进行网格划分，得到壁面重构后的热量分布。

图3 螺旋管圈展开与重构示意图

2.2 水动力计算方法

锅炉水动力特性计算的核心是获取锅炉各个管道中工质流量、压力和焓达到平衡时的各项参数。设第i回路压降为该回路中各段的重位压降（Δph）和摩擦压降（Δpf）的累加，如式(12)所示。

单相重位压降如式(13)所示。

式中，h为第i回路第j段垂直上升距离，m；ρ为第i回路第j段的工质密度，kg/m3。

汽水混合物重位压降如式(14)所示。

式中，ρl为第i回路第j段的饱和水密度，kg/m3；ρs为第i回路第j段的饱和蒸汽密度，kg/m3。

光管单相摩擦压降如式(15)所示[18]。

式中，L为第i回路第j段长度，m；D为管内径，m；G为单位面积质量流量，kg/(m2·s)；ν为第i回路第j段的比热容，m3/kg；f为湍流光管单相摩擦系数。

汽水混合物的摩擦压降如式(16)所示[19]。

式中，x为第i回路第j段含干度；νl为饱和水比热容，m3/kg；νs为饱和蒸汽比热容，m3/kg；ψ为摩擦阻力修正系数。

根据并联管组压降相等原则以及管路质量守恒，设并联管压降和总质量流量分别为ΔP和M，得式(17)压降方程和式(18)质量守恒方程。

压降方程和质量守恒方程可组成包含i+1 个方程的非线性方程组，采用离散牛顿法对方程组进行求解。由于摩擦压降和重位压降的计算方程中含有工质物性相关量，物性参数的确定也依赖于回路流量，因此方程组的解若不收敛，需用各回路流量更新方程组中的物性参数，直至计算收敛，获得最终的各管路流量和回路总压降（计算中水/蒸汽的物性参数从IAPWS-IF97 数据库[20]中获取）。

2.3 管壁温计算方法

亚临界圆管的单相传热系数采用Dittus-Boelter公式进行计算，如式(19)所示。

式中，λ为管的热导率，W/(m·K)；Re和Pr分别为雷诺数和普朗特数，实验验证范围为Re=104～1.2×105，Pr=0.7～120。

超临界光管的对流传热系数如式(20)所示[21]。

式中，αl为饱和水的对流传热系数，可采用式(19)计算得到，W/(m2·K)；μl为第i回路第j段的饱和水动力黏度，N/(m·s)；μs为第i回路第j段的饱和蒸汽动力黏度，N/(m·s)；p为工质压力，Pa；pcr为临界压力，22.115MPa；G为质量流量，kg/(m2·s)。

各回路各段的管内壁温度如式(22)[22]所示。

式中，Tf为流体温度，K；τ为管外径与内径的比值，De/D；Jn为向火侧内壁均流系数；q为壁面热流密度。

管外壁温度如式(23)[23]所示。

式中，δ为管壁厚度，m；λp为管壁热导率，W/(m·K)。

2.4 计算方法验证

研究选取的超临界600MW 四角切圆锅炉在螺旋管圈水冷壁出口高度沿周向均匀布置了70 个温度测点，为验证壁温耦合计算方法准确性，将螺旋管圈水冷壁划分为70 个回路。图4 对比了100%BMCR 负荷下螺旋管圈水冷壁出口的实测温度和计算温度。结果显示，100%BMCR 负荷下的计算温度与测点数据的最大相对偏差（即最大绝对偏差值与对应测点温度的比值）为2.7%。壁温耦合计算方法具有较高的计算精度。

图4 100%BMCR负荷下螺旋管圈水冷壁出口温度分布对比

3 水冷壁高温腐蚀预测模型

图5为水冷壁高温腐蚀预测模型的示意图。将锅炉各工况的运行参数作为边界条件输入炉膛数值模拟，结合SOx生成模型和壁温耦合计算方法，输出对应工况下的水冷壁壁温分布和近壁面H2S浓度分布。进一步将各工况的持续时间、壁温分布和H2S浓度分布代入高温腐蚀模型，累加即得到一段时间内的水冷壁高温腐蚀状态。本文选取的锅炉钢H2S 高温腐蚀数学模型如式(24)所示，各项系数及腐蚀活化能等参数由前期工作中[11]腐蚀实验数据的回归分析得到。

图5 水冷壁高温腐蚀预测模型示意图

式中，t为腐蚀时间，h；T为金属温度，K；CH2S为H2S浓度，μL/L。

根据图5的计算思路，基于Matlab GUI 平台开发了如图6 所示的水冷壁高温腐蚀预测软件，实现预测算法与用户的交互。软件界面主要包含三个模块：计算区域设定模块、水动力参数模块以及高温腐蚀预测模块。使用软件时首先确定计算区域与网格数量，如输入炉膛尺寸、高度区间、高度方向和宽度方向网格数量等；其次提供对应工况下的水冷壁内工质压力、工质的单位面积质量流量以及进入计算区域的工质焓；最后导入该工况下的数值模拟计算结果，设置工况运行时间。数值模拟数据中存储有壁面网格中心对应的三维坐标、热流密度和H2S 浓度，点击“计算腐蚀深度”按钮激活后台程序。后台程序会抓取计算区域设定模块内的炉膛边界和网格密度数据，对导入的数值模拟数据进行处理。同时抓取界面上的水动力计算参数，结合数值模拟结果，为不同工作压力下的水冷壁壁温分布耦合计算匹配对应函数，具体的函数及功能如表5 所示。根据壁温和H2S浓度分布获得特定工况下的水冷壁高温腐蚀状态。通过累加不同工况下的腐蚀深度，实现锅炉调峰运行下水冷壁高温腐蚀状态的预测。积累大量数据后，可利用机器学习算法建立锅炉运行参数（锅炉负荷、一次风率、燃烧器竖直摆角等）与水冷壁高温腐蚀程度的映射关系，从而实现水冷壁高温腐蚀状态的实时更新。

表5 壁温计算函数汇总

图6 水冷壁高温腐蚀预测GUI界面

4 模型计算结果

4.1 水冷壁高温腐蚀与壁温分布和H2S 浓度分布的关系

相关研究发现切圆锅炉下炉膛前、后墙高温腐蚀较为严重[24]，故选取100%BMCR 负荷主燃区水冷壁前墙为研究对象。图7 为该墙壁温、H2S 浓度以及腐蚀深度的分布云图。如图7(a)所示，壁温处于650～730K 之间，燃烧器区域壁温沿高度方向快速上升，且右侧壁温明显低于左侧；进入SOFA 燃尽区后壁温增速大幅降低甚至出现负增长。这是由于顺时针切圆燃烧方式下，射入炉膛的低温二次风对前墙右侧产生了冷却效果，主燃区燃烧温度高、辐射换热量大，壁温沿高度方向增速高；进入SOFA 区域后，辐射换热量降低，壁温增速缓慢。可见，热力-水动力耦合计算模型能够充分反映水冷壁的壁温分布，为高温腐蚀的准确预测提供基础。

图7 100%BMCR负荷水冷壁前墙壁温分布、H2S浓度分布以及年腐蚀深度分布

图7(b)显示高浓度H2S 主要位于燃烧器区域下部，随着炉膛高度的增加，H2S浓度波动下降，并在燃烧器层和SOFA 层之间区域出现局部的极大值。进入SOFA 区域后H2S 浓度迅速下降。通过数值模拟得到的在炉膛内H2S 浓度分布规律与孟繁兵等[25]的实炉测试规律吻合度高。图7(c)给出了前墙水冷壁高温腐蚀程度，云图显示底层燃烧器区域和SOFA 区域腐蚀程度较低，这两个区域分别对应于高H2S 浓度和高壁温。最大年腐蚀深度出现在燃烧器层和SOFA 层之间区域，为276μm。该区域特征是壁温和H2S 浓度都相对较高，当高壁温和高腐蚀气体浓度重叠时，水冷壁易发生严重的高温腐蚀。

4.2 不同负荷下水冷壁高温腐蚀特征

表1 中的四个负荷下水冷壁前墙的壁温分布、H2S 浓度分布以及高温腐蚀状态分别如图8～图10所示，年运行时间设为8760h。100%BMCR负荷下燃烧器与SOFA之间的区域腐蚀最为严重，最大年腐蚀深度为276μm。此外，燃烧器区域的高温腐蚀同样较为严重，75%THA负荷下的高温腐蚀较为严重区域与100%BMCR 负荷基本一致，但整体的腐蚀程度较低；50%THA 和35%BMCR 负荷下，高温腐蚀深度在F 层燃烧器高度位置迅速达到最大值，其中35% BMCR 负荷下最大年腐蚀深度达到256μm。各个负荷下前墙左侧的腐蚀程度均高于左侧，这与熊小鹤等[24]观测的现象一致。

图8 不同负荷下水冷壁前墙壁温分布

图9 不同负荷下水冷壁前墙H2S浓度分布

图10 不同负荷下水冷壁前墙高温腐蚀状态

结合图8 和图9 所示的壁温分布和H2S 浓度，高负荷下锅炉水冷壁的平均壁温较高，且燃烧器层和SOFA层之间区域的水冷壁缺少水平偏置二次风的保护，该区域H2S 浓度较高，使得燃烧器层与SOFA 层之间的区域成为高温腐蚀最为严重的区域。项岱军等[26]同样发现SOFA 风下方区域高温腐蚀比主燃区更严重。此外，为抑制氮氧化物的生成，高负荷下炉膛空气分级度高，一、二次风强度减弱，切圆燃烧偏斜程度大，导致燃烧器区域H2S浓度较高，高温腐蚀也较为严重。50%THA 和35%BMCR 负荷下工质处于亚临界状态，水冷壁内工质发生相变的高度位置容易出现壁温激升。如图8 所示，壁温在燃烧器区域内出现峰值，35%BMCR 负荷下壁温峰值位置更高，与壁面高浓度H2S 区域重合，这使得35%BMCR 负荷下燃烧器区域的上部腐蚀深度超过50%THA 工况。通过调整燃烧器摆角、主燃区空气过量系数等运行参数，错开高壁温区与高浓度H2S区，即可有效减缓低负荷下的水冷壁高温腐蚀程度。

4.3 调峰运行下水冷壁高温腐蚀状态

为模拟锅炉调峰运行，设计了锅炉在100%BMCR、75%THA、50%THA 和35%BMCR 四个负荷下各运行2190h 共计8760h 的情境。各工况的运行参数依照表1进行设置，使用水冷壁高温腐蚀预测GUI界面进行计算。该情境下的水冷壁高温腐蚀状态如图11 所示。其中，后墙的高温腐蚀最为显著，最大腐蚀深度为364μm。腐蚀严重区域位于F 层燃烧器高度位置以及燃烧器层和SOFA 层中间的位置。前墙、左墙和右墙的特征也与之类似，但腐蚀程度较轻。

图11 调峰运行下水冷壁高温腐蚀状态

结合图10 所示的不同负荷下水冷壁前墙高温腐蚀状态，多工况运行下水冷壁高温腐蚀状态表现为各个负荷运行时的高温腐蚀状态的时空叠加。通过本文的水冷壁高温腐蚀预测模型，能够掌握锅炉调峰运行下水冷壁高温腐蚀状态的空间分布和时间分布，为构建锅炉安全监测系统提供有利参考。同时，锅炉的运行参数能够影响水冷壁的壁温分布以及近壁面的H2S 浓度分布，进而影响水冷壁的高温腐蚀程度。运用高温腐蚀预测模型可对锅炉运行参数进行优化，为减缓水冷壁高温腐蚀提供支持。

5 结论

本文综合炉膛数值模拟、水冷壁壁温耦合计算以及包含时间维度的管壁高温腐蚀模型，提出一种适应锅炉调峰运行的水冷壁高温腐蚀预测模型。为预测验证模型可靠性，选取一台600MW 四角切圆燃煤锅炉为研究对象，计算了100%BMCR、75%THA、50%THA 以及35%BMCR 四种典型负荷下壁温分布和近壁面H2S浓度分布，进而预测出单一工况以及多工况叠加情境水冷壁的高温腐蚀状态，主要结论如下。

（1）根据螺旋管圈水冷壁的特征采用一种基于坐标变换的热量再分配方法，进一步耦合烟气侧与工质侧传热，结合管壁温计算，形成水冷壁壁温分布耦合计算方法。同时，建立了用于确定还原性气氛下的燃料硫释放以及含硫组分的相互转化过程的SOx生成模型。通过壁温耦合计算和SOx生成模型分别得到水冷壁的壁温分布和H2S浓度分布，为高温腐蚀的准确预测提供基础。

（2）不同负荷下水冷壁高温腐蚀特征存在区别，壁面腐蚀程度整体上随负荷降低而降低。100%BMCR 与75%THA 负荷下水冷壁前墙燃烧器层与SOFA层之间的区域腐蚀最为严重，最大年腐蚀 深 度 分 别 为276μm 和233μm；50%THA 与35%BMCR 负荷下高温腐蚀深度在燃烧器区域的上部迅速增加至最大值，分别为224μm和256μm。

（3）多工况运行水冷壁高温腐蚀状态表现为各工况腐蚀状态的时空叠加。在100%BMCR、75%THA、50%THA 以及35%BMCR 四个负荷下各运行2190h情境下，后墙高温腐蚀程度最严重，顶层燃烧器高度位置腐蚀深度达到364μm，其他三面墙的腐蚀特征与之类似。

（4）采用Matlab GUI 平台开发了水冷壁高温腐蚀预测软件，高度集成了炉膛数值模拟、水冷壁壁温耦合计算方法以及高温腐蚀数学模型，实现了通过锅炉运行参数和运行时间对多工况下水冷壁高温腐蚀状态的时空分布进行预测，为灵活调峰需求下锅炉安全监控系统的构建和开发提供了有力支持。