燃煤电厂煤粉燃烧过程中辐射传热计算方法研究进展与展望

刘鹏宇，李德波，刘彦丰，陈兆立，廖宏楷，冯永新

(1.华北电力大学 动力工程系，河北 保定 071003；2.南方电网电力科技股份有限公司，广东 广州 510080)

气候变化是我国乃至世界所面临的最严重挑战之一。目前煤炭资源仍是我国能源消费的主体，据统计，2020年煤炭占我国一次能源生产近7成，消费的一半以上[1]，而我国目前发电结构仍然以火力发电为主，80%的CO2排放来自电力和工业部门。为实现我国向世界承诺的碳达峰、碳中和目标，我国能源政策逐渐向风电、水电和核电等清洁能源发电方式倾斜。但清洁能源具有随机性、间歇性和不稳定性等特点[2-3]，造成清洁能源大规模并网难以消纳和电网调峰难题，因此，目前仍为我国主力发电形式的燃煤火力发电须承担电网调峰的重任[4-5]。

燃煤电厂在承担调峰任务时，须以较低的额定容量运行。低负荷时炉膛内火焰位置改变，充满度较差，容易发生偏烧现象，从而恶化水动力循环，增大水冷壁、过热器和再热器等局部过热爆管的可能性[6-8]；同时，低负荷运行会对炉内火焰中心位置造成影响，从而响锅炉折焰角和冷灰斗处积灰结渣特性[9]，干扰锅炉平稳安全运行。因此，现阶段需要进一步通过掌握燃煤机组锅炉内壁面热流值或进行炉内火焰监视来了解火焰燃烧的情况，动态掌握燃煤机组炉内火焰的燃烧情况和炉内各受热面温度分布。

开展炉内火焰监视需要电荷耦合器提供的高分辨率辐射图像，通过数值计算方法来分析得出炉膛内各受热面的温度分布情况，因此需要采用能够提供炉内高分辨率图像的辐射传热模型。典型煤粉炉中90%以上的热量以辐射传热的形式进行[10]，因此锅炉机组炉膛内辐射传热数值计算方法是目前研究的重要方向。

本文介绍燃煤电厂锅炉机组中辐射传热计算的基本方程及其假设，指出困扰辐射传递方程求解的难题，对基于求解受热面热流密度的辐射传热模型和基于空间方向的辐射传热模型进行综述和讨论，从获得精确的壁面热流和高方向分辨率辐射强度2个方面展开分析，为精确、动态掌握炉膛燃烧情况提供理论依据。

1 辐射传递方程及其假设

煤粉燃烧时炉膛内所发生的热交换形式包括传热的3种基本方式：热传导、热对流和热辐射。这3种基本传热方式协同耦合，通过燃烧放热和机械能做功，将蕴藏在煤粉中的化石能源转换为电能。相较于热传导和热对流2种方式，热辐射的机理明显不同。热辐射具有不需要通过任何介质传播、温度场可以不连续、有强烈方向性、能量按波长分布等特点，求解方程十分复杂，对其求解需要考量多方面影响因素并作不同类型的假设近似。众多辐射传热解法的假设近似可以概括性地可分为下述2类：

a)物性假设：燃煤锅炉中辐射传热与烟气的物性存在较强的耦合，因此忽略烟气的吸收、散射或发射都能大大简化辐射方程的求解过程。

b)数学假设：使用不同的数值求解方法，针对热辐射的某些特性(如温度场不连续和具有强烈方向性)，对其求解过程作各种假设，寻求在数学模型上的近似和简化。

(1)

式(1)为积分微分方程，利用传统方程式的解法很难获得其在工程问题中的解析解，尤其是式中等号右边最后一项含有积分形式的内散射，使其在较为理想的非散射性介质中也较难获得解析解，目前的计算方法大多是围绕该项开展研究而建立数学模型。

2 基于热流密度的辐射传热模型

2.1 区域法

区域法是基于净辐射法提出的，其基本原理为：将某封闭区域分为壁面和空间体积2种单元，通过计算定义为直接交换面积的各单元之间辐射直接交换，将温度和热流量作为未知量列出能量平衡式，求解得到净辐射热流。该方法通过求解单元热流密度的方法求解RTE，最先由H.C. Hottel和E. S. Cohen基于行列式方程的形式应用于无散射的灰性介质中，而后继续发展至有散射的非灰性介质中[13-14]。

2.1.1 基于一维传热分区法的区域法

为了精确研究燃煤电厂锅炉炉膛内沿高度方向的烟气温度和局部热负荷分布等物理量，并基于此对锅炉受热面的设计和运行进行优化及对水冷壁辐射受热面的工作作出合理的预测，阎维平等[13]、欧宗现等[14]均提出了如图1所示的沿炉膛高度一维热力计算方法，并应用于670 t/h、600 MW超临界锅炉机组，发现各种工况下的计算结果与机组的实际运行数据较吻合。

图1 炉膛一维分区[13-14]Fig.1 One dimensional partition of furnace[13-14]

商显耀等[15]基于炉内一维分区燃烧传热计算模型，采用区域法对二维小区进行传热建模，并作如下假设：①燃烧与传热独立进行，沿高度方向的烟气温度由一维传热分区模型得到，模型各区域中的温度由二维分区模型得到；②采用受热面为灰体；③忽略管壁积灰结渣带来的新增热阻。基于区域法建立的二维模型如图2所示。

图2 区域法分区[15]Fig.2 Schematic diagram of zone method partition[15]

基于上述假设推导建立能量平衡方程，并应用于1 000 MW超超临界锅炉，发现耦合辐射传热计算值的平均误差在100%工况下为4.69%。虽然计算误差已经较低，但仍高于通过选取管内换热系数后代入壁面热负荷的管内换热模型计算结果。

2.1.2 二维及三维区域法

二维及三维区域法的计算基于每个单元之间的直接交换面积求解，而直接交换面积的求解则需要求解介质单元吸收率和系统形状函数，还需引入有效辐射函数，使得区域法求解较为繁琐。因此在二维和三维计算中，简化直接交换面积及求解过程是区域法的核心要素。

马金凤等[16]将仅针对二维吸收性介质的广义对称性[17]推广至二维、三维灰性散射性介质中，广义对称性的分段示意图如图3所示，其中：s表示各表面微元；g表示各容积单元；r1、r2和w1—w4为表面微元间的距离。通过该种方法，可以有效地将图3(a)中12个区域的144个交换面积简化为20个，大幅降低了区域法的计算量。采用简化后的直接交换面积计算方法，针对某电厂自然循环锅炉中含弥散性介质烟气的炉膛进行数值计算。

图3 二维及三维广义对称性[16]Fig.3 Simplified schematic diagram of two-dimensional and three-dimensional generalized symmetry[16]

区域法具有严谨的计算机理，计算精度较高，通常被用于其他模型的精确性验证。但区域法仍具有较大的局限性：直接辐射面的计算带来高精度结果的同时引入了庞大的计算量，计算机编程实现时需要较大的内存空间储存交换区域的数据，炉内流动耦合燃烧的计算方式会显著增加区域法的耗时；无法处理非灰体的辐射传热问题，对散射性介质的计算误差较大，而燃煤锅炉炉膛中煤粉和燃烧后的烟气介质为散射性介质，若忽略其造成的影响将引入较大的误差。因此，进一步简化直接交换面积和拓展该方法在非灰体中的应用，以及将区域法与其他方法进行耦合求解，是未来重要的研究方向。

2.2 蒙特卡洛法

蒙特卡洛法并不特指某个具体的算法，任何一个采用合适的统计抽样方法求解数学问题的方法通常都称作蒙特卡洛算法[18]。蒙特卡洛法的提出最初是用来研究核武器潜在行为，而后被Howell引入辐射传热的求解中[19]。

蒙特卡洛法打破了将辐射能视为连续能量的思路，其基于统计抽样的方法耦合了介质吸收、介质发射、固体壁面发射、固体壁面反射和吸收的模拟。热辐射所发射的辐射能在与其他介质作用前，沿发射路径直线传播。根据该特性，蒙特卡洛法需从发射源开始追踪具有统计意义的抽样能束，直到该能束的能量被完全吸收。蒙特卡洛法是目前求解多维各项同性、异性RTE的有效手段之一，主要应用于求解壁面热流分布等不涉及高方向分辨率的辐射强度的研究[20]。

在蒙特卡洛法的计算时间研究方面，阮立明等[21]选取2个随机数生成程序，研究二维正方形介质腔内辐射传递的计算时间，发现在不考虑编译软件和计算及性能的情况下，介质的物性(如光学厚度和散射相函数等)、壁面的物性(如壁面发射和反射特性等)、能束数量等因素有密切关系。除此之外，阮立明等还研究了基于二维正方形各向异性散射、吸收性灰体介质，发现其温度的计算误差随能束数的增加明显减小。同类型的研究中，成珂等[22]在细致推导蒙特卡洛法的基础上，对处在无限大黑体环境中的等温辐射表面和黑体表面的辐射进行了研究，发现随着跟踪能束的增加，蒙特卡洛法的精度显著提升至小于0.3%。

在燃煤电厂锅炉炉膛内辐射传热的数值计算中，夏钧[23]等对330 MW旋流对冲燃烧锅炉进行计算，求得炉膛中心截面温度如图4(a)所示，炉膛中心温度可达1 600 K，炉内温度近似呈对称分布；炉膛中心纵向截面温度场计算值与实验值的比较如图4(b)所示，可以看出二者吻合较好，验证了蒙特卡洛法在燃煤电厂锅炉机组中的应用。

图4 炉膛中心截面温度分布及计算值与实验值比较[23]Fig.4 Temperature distribution of central section of furnace and comparison between calculated value and experimental value[23]

蒙特卡洛法作为一种数学统计的方法，难免存在误差，但随着电子信息技术的飞速发展，其随机抽样的能束已经可以达到几亿束，工程问题误差可忽略不计，是未来十分具有潜力的辐射传热计算方法。

3 基于辐射方向的辐射传热模型

3.1 基于离散坐标法的辐射传热模型

离散坐标法的提出最初是为了研究天体物理问题[24-25]，而后被用来求解一维辐射传递问题。

离散坐标法在求解辐射传递方程时，划分好网格节点后，根据权值的思想将整个4π立体角离散至既定的坐标方向上，并将其难于求解的积分形式的内散射相用数值积分近似表示，从而在每个离散坐标方向上都可以求解辐射传递方程，并得到最终的辐射强度和温度分布。

离散坐标在其角度划分方法上采取高斯积分的原则选择射线方向。若积分要达到N阶，则整个4π立体角将被划分为N(N+2)个区域(坐标)数，这种方法一般被称为SN法。因此，若要保证辐射方程求解的精度，需要保证较多的离散空间方向。

刘林华等[26]通过三维矩形炉膛计算来对比区域法和离散坐标法，验证了离散坐标法(S6法)的准确性，并根据离散坐标法的模拟结果分析辐射求解中精度的影响因素，发现辐射传热的数值计算精度受介质的吸收系数、散射系数、散射相函数、光谱特性和壁面积灰结渣情况的影响，且介质的散射系数和壁面积灰结渣层的厚度是主要影响因素。因此在燃煤电厂锅炉机组的辐射传热计算模型中，忽略介质的散射特性将会引入较大的误差。王成军等[27]采用离散坐标法模型对某旋流燃烧室开展数值模拟工作，得到燃烧室温度分布，通过与区域法的对比，验证了离散坐标法的精确度。李力等[28]编制离散坐标法的计算程序，对一个矩形空间的三维纯辐射传热进行数值计算，并与区域法计算结果进行对比，同样验证了离散坐标法的正确性。

RTE中较难求解的散射相函数还必须满足散射相能量守恒定律得出的归一化条件，即

(2)

若燃煤锅炉燃用会产生较多小颗粒的煤质或烟气中未燃尽碳含量较高时，散射相函数将会变得更加复杂，此时SN法很难满足归一性条件。魏小林等[25]在SN法求解RTE的基础上，提出了命名为QN法的新坐标离散方法。QN法采用球面的经纬度进行离散，N的物理意义不再是阶数，而是第一卦限内离散的球面数。经分析计算可知，QN法可以较好地满足复杂的散射相函数的归一化条件。采用QN法对充满灰粒、未燃尽碳等小颗粒的长方体煤粉炉膛进行计算，炉内中心温度的试验值和计算值比较如图5所示，可以看出二者总体吻合较好，在燃尽段计算值略低于试验值。

图5 基于QN法的炉内中心温度试验值和计算值比较[25]Fig.5 Comparison of experimental and calculated values of furnace center temperature based on QN method[25]

在对离散坐标法的改进方面，杨占春等[29]认为离散坐标法按离散坐标的方向而不是实际离散方向进行积分，会产生较大的误差，因此提出基于传播方向各面投影几何关系的控制体区域划分方法，并分段对RTE进行积分，该方法称为分段积分离散坐标法。由该方法的三维矩形炉膛计算结果与区域法计算结果的对比可知，该方法的计算值更接近区域法，在高温区吻合十分良好且能够提高运行效率。

离散坐标法将辐射强度的方向离散至坐标方向上，从而确保了在一定立体角内其方向不变且强度均匀，是离散坐标法的显著优势；用数值积分表示内散射相，较大程度上简化了RTE的求解。但离散坐标法中离散坐标的选取是人为干涉的，会在坐标选取和权值构造方面引入人为误差，假散射和低离散阶数带来的误差也是该方法的弊端所在，加之离散坐标法在数值模拟中无法处理炉膛内散射问题，因此对离散坐标法的改进是未来重点研究方向。

3.2 基于离散传递法的辐射传热模型

离散传递法用确定的一根特征辐射射线来代替该辐射射线所在立体角的全部辐射能，不必计算多重积分，只需将发射的辐射离散为尽可能多的能束，则计算精度可以达到多重积分精度。特征辐射射线从发射网格边界面穿过内部介质到达吸收边界面被吸收，其吸收、发射源均为边界网格面，并在发射、吸收边界面上达到辐射能平衡。

离散传递法求解的重点是确定特征辐射射线的路径与其在介质中传递的路径长度，因此需确定特征射线从发射源至吸收源的交点及交点之间的距离，并沿此射线积分辐射传递方程进行求解，其精度的控制依赖于发射点所离散的能束数，能束数越多计算精度越高。

顾明言等[30]基于三维圆柱腔体建立离散传递法的辐射传热求解模型，将圆柱腔体上发射的特征射线划分为端面发射和柱面发射2类，通过数值积分求解辐射源项，并与已知的精确解进行比较，发现其与精确解的误差不大于4%，验证了离散传递法的准确性。将该方法应用于煤粉燃烧喷口B的燃烧过程，得出如图6所示的预报温度与实验结果对比，验证了离散传递法在燃煤锅炉机组中的适用性。图6中：X为轴向距离(m)；d3为燃烧器二次风风壳半径，d3=0.056 m。

图6 预测温度分布与实验值比较[30]Fig.6 Comparison between predicted temperature distribution with experimental value[30]

于娟等[31]针对离散传递法提出一种降维后的离散传递法，基于如下假设可认为辐射传热只在横断面的二维区域内进行：①沿炉膛高度方向温度变化较小，即温度梯度较小；②减弱系数较大。降维后的离散传递法舍去了立体角和三维的辐射强度概念，采用平面角表征发射面法线与微元能束中心轴的夹角，采用单位时间、可见面积和平面角描述辐射强度。于娟等采用降维后的离散传递法，对圆筒形燃煤试验炉进行燃用3种褐煤的数值计算，计算值与测量值比较如图7所示(线状数据为计算值，点状数据为测量值)，可以看出降维后的离散传递法对具有二维特征的辐射传递问题具有较准确的计算结果。

图7 降维离散坐标法的计算值与测量值比较[31]Fig.7 Comparison of calculated and measured values of dimension reduction discrete coordinate method[31]

离散传递法用特征辐射射线的方法代替发射源表面沿立体角的所有辐射能，避免直接计算空间多重积分，简化了RTE的计算。但为确保一定的计算精度，需要取大量的能束数来无线逼近多重积分值，可能反而引入较大的计算量。为了减少大量能束数的运算量，可以适当增大网格密度，根据工程实际需要确定所需的能束数和网格密度。离散传递法在数值模拟中应用时同样存在无法处理炉膛内散射的问题，因此，基于离散传递法的散射改进模型也是未来数值模拟计算中的重要改进研究方向。

3.3 基于DRESOR法的辐射传热模型

黄志峰等[33]基于各向异性纯散射的一维板状灰介质对DRESOR法进行验证，通过变化光学厚度τL得到关于相函数的无量纲量Ψ，对比无量纲量Ψ在前向和后向散射性介质的辐射热流计算值发现，DRESOR法对各向异性散射计算具有很高的准确性。光学厚度为1.0右边界处的辐射强度Ψ值如图8所示，证明DRESOR法可以计算得到高方向分辨率的辐射强度。

图8 光学厚度为1.0右边界处的辐射强度[33]Fig.8 Radiative intensity at right boundary with τL=1.0[33]

程强等[34]利用此方法对一维灰性平行平板介质进行3种工况的求解。3种工况中，工况3的发射源是工况1和2的耦合，从3种工况的计算中得出了辐射强度的可加性原理，该原理为计算机编程实现DRESOR法提供了理论依据支撑。基于该算例，程强等进一步将边界处计算所得的辐射强度代入RTE边界条件表达式中进行对比验证，误差小于1%，说明DRESOR法对RTE的求解结果可以自动满足边界条件，大大提升了未来基于燃煤电厂受热面壁温、热流计算的便捷性。

黄志峰等[35]将DRESOR法从一维提升至三维空间内，并对各项同性和各项异性散射介质进行了推导计算，均具有很高的计算效率，而且可以改善离散坐标法中由于RTE离散方向与真实能束分布不一致所引起的散射效应。张险等[36]在圆柱坐标系下根据半径、圆柱底面投影与x轴方向夹角、圆柱所在高度z，推导了DRESOR数及其计算方法，并成功在各项同性、异性散射均匀折射率介质中进行了验算，验证了文献[37]提出的DRESOR法求解RTE的基本公式在圆柱坐标系下的应用，拓宽了DRESOR法的应用范围。

基于蒙特卡洛法提出的DRESOR法凭借DRESOR数仅与所研究对象结构和辐射参数相关的特性，避免了由于研究对象温度改变导致的不稳定性，通过高方向分辨率辐射强度的求解优势，可以较容易地得到空间内任何方向的辐射强度，是未来开展炉内火焰和受热面温度监测较好的理论方法。但目前该方法的研究仍然停留在三维圆柱炉内，在燃煤电厂锅炉机组中的应用较为鲜见，未来采用该方法对燃煤电厂锅炉机组开展辐射传热计算是重点的研究方向。

4 结论与展望

燃煤电厂锅炉机组中辐射传热计算模型可以依据空间方向求解划分为2类：一类为着重求解辐射强度或壁面热流的辐射传热模型，如本文列举的基于一维传热分区法的区域法模型和蒙特卡洛模型；另一类是考虑到辐射方向而对空间方向进行离散划分的辐射传热模型，如本文列举的离散坐标法、离散传递法、DRESOR法模型。上述模型大部分均在燃煤电厂锅炉机组中取得了较好的辐射传热计算结果，为锅炉的工程实际提供了一定的指导，但针对目前承担调峰重任的燃煤锅炉机组来说，还有如下问题需进一步研究：

a)我国燃煤电厂深度调峰时开展实时炉内火焰监测，需要提供炉内高分辨率图像以供电荷耦合器采用的辐射传热模型。因此，考虑辐射方向的离散坐标法、离散传递法及DRESOR法辐射传热模型在深度调峰燃煤锅炉中的应用是未来的研究热点；针对深度调峰燃煤锅炉的辐射强度求解或壁面热流辐射传热模型的改进也是重点研究方向。

b)离散坐标法、离散传递法目前在数值模拟中有一定的应用，但缺点是无法处理存在散射的问题，因此目前燃煤电厂数值模拟中普遍采用考虑散射修正且精度受模型几何复杂性影响的P-1辐射模型[38-40]。基于散射修正的离散坐标法、离散传递法和其他方法应通过用户自定义函数在数值模拟中进行耦合计算，以提高数值模拟的准确度，为燃煤电厂的实际运行提供指导。

c)区域法的缺点是在计算非灰体和散射性介质时会引入较大的误差且耗时较高。若结合区域法与燃煤锅炉机组辐射传热进行计算，需进一步简化直接交换面积的计算并开发区域法的改进算法，或将区域法与其他方法进行耦合计算，以减少计算所需耗时。

d)蒙特卡洛法的优点是在复杂介质和不同炉形应用时有较为优异表现，展示出良好的发展前景。缺点是基于能束追踪的求解方法在炉内温度、热流情况未知的情况下需要进行迭代逼近处理，为满足精度要求，其高达几亿束的能束追踪计算将大大增加计算耗时，对计算机处理器的性能要求较高。未来可以结合蒙特卡洛法和商业计算流体力学软件Fluent对炉膛内部温度进行数值模拟，在已知温度场的条件下开展辐射传热计算。

e)离散坐标法和离散传递法的优点是可以获取辐射能在空间方向上的计算结果，可利用这类方法开展炉内火焰和受热面温度的监测。但这类方法为了确保精度，离散的能束较多，造成迭代时间较长。因此，在开展实时监测方面，需开发基于这类方法的改进算法，并考虑离散坐标法的散射效应以及离散传递法的能束数与网格密度的影响，避免引入较大的误差；数值模拟方面，均需开发考虑基于散射影响的改进算法，以提高燃煤电厂锅炉机组数值模拟的精确性。

f)基于蒙特卡洛法的DRESOR算法具有与离散坐标法和离散传递法相同的获取任意空间方向上辐射强度的能力，其优点是不受研究对象温度的干扰，其较快的求解速度适用于开发炉内火焰、温度监测系统。但目前针对现役燃煤锅炉机组的研究较少，需在今后的研究中多开展相关研究，以配合燃煤锅炉机组深度调峰的安全运行需求。