生物质作为一种清洁的可再生能源，可以有效减少生产活动中的CO₂排放，从而应对目前较为严峻的全球变暖问题^[1]. 而生物质碱金属含量较高，容易引起积灰、结渣和腐蚀等问题，进而威胁设备安全性和可靠性^[2-3].

近几年，国内外学者就锅炉运行中关键的换热面积灰结渣问题展开了大量实验和模拟研究. 其中计算流体动力学（computational fluid dynamics, CFD）可以经济、有效地对换热面积灰倾向等方面进行预测和验证. 蔡永铁^[4]针对循环流化床尾部烟道的积灰问题进行了热态飞灰沉积实验和数值模拟工作，在稳态条件下研究了粒径、风速和探针管径对飞灰沉积的影响. 杨琦等^[5]对飞灰沉积数值模拟现状进行了分析和总结，对该领域主要研究成果进行了概括. Ma等^[6]在CFD模型中建立了相应的颗粒碰撞和粘附模型，用以模拟棕榈油渣和煤粉掺混燃烧过程的飞灰沉积，从而预测换热管壁面结渣特性. 在飞灰沉积过程中，换热管壁面外轮廓会随沉积过程进行而有所变化，进而会影响附近流场和颗粒的沉积. Tomeczek等^[7-8]采用Fluent中动网格技术模拟了对流区域换热管的动态沉积过程. García等^[9-10]针对余热锅炉换热面积灰进行了动网格沉积模拟，探究了不同管排列方式对积灰的影响.

本文针对煤粉和玉米秸秆（玉米秸秆掺比分别为0、5%和10%）掺烧形成的飞灰进行非稳态积灰模拟，并结合动网格技术，考虑灰渣表面的温度变化. 分别从模拟过程中的灰渣增长高度、灰渣表面温度、热流密度等方向将模拟结果和实验结果进行对比.

煤粉-玉米秸秆掺烧实验在300 kW中试实验台^[11]（见图1）中进行. 整体实验装置根据功能可以划分为若干个模块：送风系统、给料系统、点火系统、炉膛本体、烟气处理系统、温度检测平台、灰沉积探针系统、CCD相机在线监测系统和烟气采集系统. 炉膛高度为3 950 mm，内径为350 mm. 将1根冷却探针在炉温稳定时伸入炉膛开始结渣实验，另一侧由CCD相机实时捕捉结渣的图像. 热流参数等数据可以通过探针表面伸入的2根热电偶进行离线测量，其示意图如图2所示. 图2展示了结渣探针和其表面的结渣的相对位置，并对热电偶测控位置进行了说明.

Ma等^[12]针对生物质掺烧比例对积灰结渣的影响进行了炉膛内实验研究. 本文对该实验结果数据进行处理分析，包含灰渣表面温度计算和灰渣导热系数计算，将其用于后续模拟的改进和验证. 实验中使用的燃料为神华煤和玉米秸秆，其物性参数见表1. 实验过程中，给料速度保持45 kg/h恒定，该用量可用于近似计算模拟中的飞灰质量流率. 烟气来流方向上热流密度大小可通过该方向上布置的A孔和B孔2组热电偶分别测得管外壁和内壁的温度，然后按下式计算：

(1) $q = \frac{{\lambda \left( {{T_{\text{1}}} - {T_2}} \right)}}{{r{\text{ln }}({{{r}}_2}/{{{r}}_{\text{1}}})}}.$

式中：q为热流密度，λ=16 W/（m·K）为探针导热系数，T₁和T₂分别为内部热电偶测得的温度和外部测得温度，r₁=15 mm和r₂=19 mm分别为内部和外部热电偶至探针圆心的距离，r为探针半径（r=20 mm）.

另外，烟气来流方向上通过灰渣的热流可以通过沉积表面侧的辐射和对流换热或内部热传导进行理论计算，这里忽略了圆周方向的热传导. 根据傅里叶定律，径向上热流密度保持不变，可以列出如下公式：

(2) ${q_{\rm{o}}}=\overline {h}\left( {{{T}_{\rm{g}}}-{T_{\rm{s}}}} \right)+\varepsilon \sigma ({T_{\rm{g}}}^{4}-{T_{\rm{s}}}^{4}),$

(3) $\overline {h} = \overline {Nu} \times k/{D_{\text{p}}},$

(4) $\overline {Nu} {\rm{ = 0}}{\rm{.683}} \times R{e^{{0}{\rm{.466}}}} \times P{r^{{\rm{1/3}}}},$

(5) ${q_{\rm{i}}}={\overline \lambda _{{\rm{dep}}}}\frac{{{T_{\rm{s}}}-{T_{2}}}}{{\Delta {\delta _{{\rm{dep}}}}}},$

(6) ${q_{\rm{o}}}={q_{\rm{i}}}.$

式中：q_o和q_i分别为灰渣外侧和灰渣内部计算得到的热流密度， $\overline h$ 为烟气和结渣探针表面对流换热系数，T_g和T_s为烟气侧温度和灰渣表面温度，ε为灰渣黑度，σ为黑体的辐射常数（5.67×10⁻⁸ W/m²·K⁴）， $\overline {Nu}$ 为平均努谢尔特数，根据雷诺数和普朗特数计算，k为烟气导热系数，D_p为探针外径， ${\overline \lambda _{{\rm{dep}}}}$ 为灰渣有效导热系数， $\Delta {\delta _{{\rm{dep}}}}$ 为灰渣厚度.

可以发现，在以上参数中，除灰渣表面温度和有效导热系数外，其余参数均可根据实验获得，可以根据实验数据对这2个参数进行进一步分析，前者可通过式（2）进行迭代计算，后者可通过式（5）进行间接求解. 对玉米秸秆和煤粉不同掺烧比例（玉米秸秆掺烧比例R=0、5%和10%）实验工况中的灰渣表面温度和灰渣有效导热系数进行整理，图3即为处理得到的3个工况下灰渣表面温度随时间变化图. 可以看出，随着探针表面灰渣的逐渐形成，灰渣表面温度先升高然后基本保持不变. 3种工况下120 min时的灰渣表面温度分别为1 353、1 309和1 430 K，该结果将用于和后续数值模拟结果进行比较，以此判断数值模拟准确性. 图4为处理得到的3个工况下灰渣导热系数和灰渣最高厚度间的关系，由图4可以看出，不同玉米秸秆掺烧比例下，随着厚度的不同，灰渣导热系数的变化规律一致：均随灰渣厚度的增加而增加. 相较纯煤燃烧工况而言，R=5%工况下灰渣导热系数增加幅度最大，这是因为灰渣微观结构的不同. 值得注意的是，R=10%工况下灰渣导热系数较小，这是因为在该工况下，灰渣烧结层出现了较大的气孔. 灰渣导热系数的计算结果可通过UDF载入的方式，对数值模拟中换热过程进行校正，从而提高模拟精度. 需要指出的是，本研究未考虑灰渣厚度不均的情况，仅选择来流方向上顶点作为计算的依据. 考虑到实际灰渣具有分层结构（原生层、烧结层和渣层），本研究数值模拟中未考虑分层结构，仅考虑了“有效导热系数”. 灰渣的分层结构在数值模拟中呈现，可作为下一步工作重点.

本研究计算域的长度和宽度分别为1 200 mm和350 mm（如图5所示），烟气和注入飞灰从左至右进入计算域，横掠计算域中心放置的结渣探针（外径D为40 mm，内径为38 mm）. 该二维计算域对应于文献[12]实验中的一维炉炉膛结构（直径为350 mm）. 计算域整体的网格划分采用非结构化网格^[13]，为适当减少计算负担并且准确计算探针附近换热，在探针附近进行网格适当加密：其中环形探针固体域网格大小一致，为0.2 mm；流体域采用Gambit中Size Function进行局部加密，起始网格大小为0.2 mm，增长因子为1.27，最大网格大小为5 mm. 整体网格划分及探针附近网格划分效果图如图6所示，非结构化网格的好处是便于动网格技术的运用，本文模拟计算域共计49 636个计算单元. 在计算域中，左端面设置为速度入口，并作为释放飞灰颗粒的入口，设置烟气速度为5.6 m/s，烟气成分采用混合气体. 上、下端面设置为wall，温度为1 573 K. 右端面设置为出口outlet边界条件. 在壁面上的颗粒行为方面，炉膛壁面统一设置颗粒反弹，而探针外壁面可通过DEFINE_DPM_EROSION宏进行特殊化设置.

本研究只考虑生物质掺烧灰分对灰渣沉积的影响，因此将3种工况下的飞灰粒径分布视为一致：满足Rosin-Rammler分布，平均粒径为16 μm，粒径分布为1~62 μm，扩散系数为0.7. 入流烟气温度设置为1 573 K，采用混合气体，体积分数一定( $\varphi $(O₂)=5.0%， $\varphi $(O₂)=16.0%， $\varphi $(O₂)=0.2%， $\varphi $(N₂)=78.1%)，流速为5.6 m/s.

本文根据探针外壁面起始网格大小建立5个网格数不同的流域划分来实现对网格无关性的验证，选取的初始大小分别为1.00、0.75、0.40、0.20和0.10 mm，对应的网格数依次增加，分别为37 044、37 958、41 356、49 636和73 598. 在这5个网格中均执行一个时间步长的计算，记录下探针表面的飞灰沉积质量，结果如图7所示，其中，r_dep为沉积质量流率. 可以发现，随着网格细化，飞灰颗粒的探针表面充分换热，颗粒温度降低，形成的飞灰沉积减少，但是0.1 mm的网格过于细化，造成计算负担，因此选择0.2 mm作为探针表面网格大小.

本研究采用欧拉-拉格朗日气固两相流模型来研究烟气夹带飞灰的两相流问题. 采用欧拉模型描述连续相（高温烟气），采用拉格朗日随机轨道模型预测飞灰运动轨迹. 烟气入口流速为5.6 m/s，经计算，管内流动雷诺数为8 732，绕流雷诺数为998，均在湍流范围内. 采用可实现k-epsilon湍流模型来模拟连续相的流动，该模型可以同时处理剪切流、漩涡、分离等流动问题，较标准k-epsilon计算精度和准确度有明显提高.

飞灰颗粒由入口入射，入射速度为5.6 m/s，起始颗粒温度设置均为1 573 K. 颗粒在流体中的运动方程为

(7) $ \frac{{{\rm{d}}{{{u}}_{\rm{p}}}}}{{{\rm{d}}t}}=\frac{{{\rm{18}}\mu {{{C}}_{\rm{d}}}Re}}{{{\rho _{\rm{p}}}{{d}}_{\rm{p}}^{2}{\rm{24}}}} \times \left( {{{u }}-{{{u}}_{\rm{p}}}} \right)+{{ g}}\frac{{{\rho _{\rm{p}}}-\rho }}{{{\rho _{\rm{p}}}}}+{{ F}}. $

式中：u_p、 ${\rho _{\rm{p}}}$、d_p分别为飞灰颗粒的速度矢量、颗粒密度和粒径；u、ρ分别为烟气的速度矢量和密度；g为重力系数；C_d为曳力系数，将飞灰简化为球形颗粒，曳力系数的计算公式为

( ${C_{\rm{d}}}={{{a}}_1} + {{{a}}_2}/Re + {{{a}}_3}/R{e^2},$

其中，a₁、a₂和a₃为常数；F为Fluent中可设置额其余受力（比如热泳力、Saffman力等）.

除了曳力以外，本文还考虑了重力和热泳力对飞灰颗粒运动的影响. 尤其在近壁面，灰渣表面附近的烟气温度和灰渣表面温度会形成温度梯度，从而影响飞灰沉积，因此在模拟中需要考虑热泳力. Saffman力仅适用于颗粒为亚微米级的气固两相流^[14]，本文中颗粒最小粒径仅为1 μm，因此可以忽略Saffman升力的影响. 湍流脉动力对飞灰颗粒的影响通过Stochastic模型实现，即在Fluent中激活了Discrete Random Walk模型和Random Eddy Lifetime.

在沉积判定准则方面，采用基于熔融相比例的灰沉积模型，该模型被国内外许多学者使用^[15-16]. 采用软件Factsage5.2对不同温度（1 113～1 683 K）下3种燃烧形成飞灰中熔融相比例进行计算（如图8所示，η_p为熔融相质量占总质量的百分比），采用软件中的“平衡计算”模块，上使用SLAGA选择熔融相的成分，该相中包含的化合物有MgO，FeO，MnO，Na₂O，SiO₂，TiO₂，Ti₂O₃，CaO，Al₂O₃，K₂O. 输入的成分包括气体成分和灰成分，其中掺混工况下的飞灰成分可以根据神华煤和玉米秸秆两者的灰成分和灰分含量等进行加权计算. 可以看出1 300～1 550 K区间内，随着R增加，同温度下熔融相比例η_p有所增加. 计算得到的熔融相比例曲线可用于进行分段拟合，拟合得到的熔融相比例-温度函数可载入到UDF中，从而实现颗粒撞击壁面时的粘附判定，判定规则如下：

(8) ${\eta _{{\rm{stick}}}}={\eta _{\rm{p}}}\left( {{{{T}}_{\rm{p}}}} \right)+\left( {{\rm{1 - }}{\eta _{\rm{p}}}\left( {{{{T}}_{\rm{p}}}} \right)} \right){\eta _{\rm{s}}}\left( {{{{T}}_{\rm{s}}}} \right).$

式中： ${\eta _{{\rm{stick}}}}$ 为颗粒p的粘附概率； ${\eta _{\rm{p}}}\left( {{{T}_{\rm{p}}}} \right)$ 为颗粒温度为T_p时颗粒中熔融相比例； ${\eta _{\rm{s}}}\left( {{{{T}}_{\rm{s}}}}\right)$ 为撞击壁面温度T_s对应颗粒中熔融相比例.

采用Fluent中的动网格技术实现灰渣生长. 将粘附在每个面的颗粒质量进行累积加和，通过F_UDMI功能储存在每个面上. 网格节点位移可以通过毗邻的2个面计算得到的灰渣沉积量进行平均计算（如图9所示）：

(9) $\Delta {{{h}}_{{i}}}=\frac{{\rm{1}}}{{2}}\left( {\Delta {{{h}}_{{{j - 1}}}}+\Delta {{{h}}_{{j}}}} \right),$

(10) ${{A}}_{{i}}={{{A}}_{{{j - 1}}}}+{{{A}}_{{j}}}.$

式中： $\Delta {h_i}$ 为节点i的位移大小， $\Delta {h_{j - 1}}$ 和 $\Delta {h_j}$ 分别为面j−1和面j的存储增长高度，A_i为节点i的位移方向，A_j−1和A_j为面j−1和面j的法向矢量.

储存在UDM中的灰渣高度根据阶段累积飞灰质量、飞灰密度、孔隙率等进行计算：

(11) $\Delta {{{h}}_{{j}}}{\rm{ = 120}} \times \frac{{{m_{{\rm{dep}}}}}}{{{\rho _{\rm{p}}}\left( {{\rm{1 - }}{\varepsilon _{{\rm{dep}}}}} \right)}}.$

式中：m_dep为Fluent计算得到的面上累积的飞灰质量， ${\rho _{\rm{p}}}$ 为飞灰颗粒的密度，为2 500 kg/m³， ${\varepsilon _{{\rm{dep}}}}$ 为灰渣孔隙率：

(12) ${\varepsilon _{{\rm{dep}}}} = 1 - \left[ {\left( {1 - {\varepsilon _0}} \right) + \frac{{{R_{\rm{l}}}}}{{{R_{\rm{s}}}}}\left( {1 - {\varepsilon _0}} \right)} \right].$

式中： ${\varepsilon _0}$ 为原始孔隙率， ${R_{\rm{l}}}$ 和R_s分别为灰渣中液相和固相的比例，该比例的计算可通过化学平衡和液相密度等信息进行估计.

由于探针圆周方向上飞灰颗粒碰撞的随机性，探针各面间灰渣沉积量会产生差异；当差异过大时，会产生负体积致使计算中断. 本研究通过如下措施来避免负体积产生：1）“吸收”局部较大沉积量. 把探针圆周方向的面分成45组，每组内部沉积量平均，然后通过权重计算缩小各组间沉积量差距（根据式（13），其中下标m表示面组的编号，下标为k−2～k+2的值为相邻5个组的权重因子，主面组的权重最大，距离主面组越远的面组权重越小，5个权重因子分别为2、2、3、3和6）；2）灵活调整组间节点的位置. 即当组间节点对应夹角小于167°时，将该节点移动至所在三角形重心位置，坐标变换如式（14）所示.

(13) $\Delta {{h}}_{{m}}'=\frac{{{\varphi _{{{k - 2}}}}{{{h}}_{{{m - 2}}}}+{\varphi _{{{k - 1}}}}{{{h}}_{{{m - 1}}}}+{\varphi _{{k}}}{{{h}}_{{m}}}+{\varphi _{{{k + 1}}}}{{{h}}_{{{m + 1}}}}+{\varphi _{{{k + 2}}}}{{{h}}_{{{m + 2}}}}}}{{{\varphi _{{{k - 2}}}}+{\varphi _{{{k - 1}}}}+{\varphi _{{k}}}+{\varphi _{{{k + 1}}}}+{\varphi _{{{k + 2}}}}}},$

(14) $ {{x}}_{i}'=\displaystyle\frac{1}{{{3}}}\left( {{{{x}}_{{{i - 1}}}}+{{{x}}_{{i}}}+{{{x}}_{{{i + 1}}}}} \right),\;\; {{y}}_{i}'=\displaystyle\frac{1}{{{3}}}\left( {{{{y}}_{{{i - 1}}}}+{{{y}}_{{i}}}+{{{y}}_{{{i + 1}}}}} \right){{.}} $

如图10所示为Fluent中结渣模拟流程示意图. 数值模拟采用非稳态模式来模拟随时间动态变化的结渣过程. 在一个时间步长中，Fluent自带解算器先计算流场并求解颗粒运动方程，在此过程中，沉积在探针表面的飞灰质量存储在相应面上. 在经过一段时间周期后（200个时间步长，每个步长0.005 s），每个面上存储的飞灰质量值重置为0，并且执行动网格更新程序和表面温度的迭代计算更新. 通过该模拟方法，并行计算1个周期（1 s）需要1.5 h左右，因此持续100 min的沉积过程计算量过大. 为缩短模拟计算时间，假设每个时间阶段（2 min）内灰渣生长速率恒定，只需将每个周期的撞击量、沉积量等参数乘以120即能计算每2 min灰渣增长情况，与此同时，动网格更新和灰渣表面温度更新也是每2 min执行一次. 因此，本文中提及2种时间尺度，一种是模拟时间t（单位为s），另一种是实际沉积时间t′（单位为min）.

3种工况下飞灰颗粒在探针表面的撞击和沉积质量流率随时间变化如图11所示，其中，r_m为质量流率，r_imp为撞击质量流率，沉积效率的定义为同时间段内沉积质量和撞击质量的比值. 由于在数值模拟过程中，颗粒随机撞击，质量流率会随时间产生波动. 所有工况中撞击质量流率均在1.3×10⁻⁶~1.5×10⁻⁶ kg/（m·s⁻¹）证明三工况除灰分不同，其他条件均保持一致. 从曲线整体趋势来看，撞击质量流率随时间推进缓慢减小，而沉积质量流率先增加后保持稳定. 前者是由沉积物动网格引起的（无动网格工况中撞击质量流率基本保持不变，如图11中第三工况图中所示），后者是因为灰渣表面温度逐渐升高，粘附概率进一步提高. 在沉积效率方面对3个工况分别进行计算可以发现：R=0、5%和10%三种工况在t=50 s时的沉积效率分别为75%、80%和87%. 由此说明，随着玉米秸秆掺比的增加，沉积效率明显增加，该结果进一步说明了生物质掺烧引起的灰分差异会增加灰渣的形成倾向，而且影响较大.

3种工况下探针中心位置上灰渣增长高度随时间的模拟值和实验值如图12所示. 可以看出，在灰渣生长初期（0～20 min），模拟和实验值较为一致，但在100 min时灰渣厚度方面，三工况模拟值均比实验值大，这是因为本文在结渣机制中未考虑灰渣脱落的影响，目前对于数值模拟中灰渣脱落的研究尚未成熟. 在数值模拟中，飞灰质量流率恒定且考虑因素较实验而言比较单一，因此图12中灰渣增长的模拟值基本呈线性增长，纯煤燃烧工况下的灰渣生长率为0.055 mm/min，R=5%和R=10%工况下分别为0.062和0.069 mm/min，说明玉米秸秆掺烧形成的飞灰会加剧结渣的生成.

动网格执行下的灰渣形貌随时间的变化如图13所示. 在三工况中，灰渣主要生成在探针迎风面，这是由于大量飞灰颗粒撞击粘附所致，该现象与实验中拍摄到的图像保持一致^[12]. 通过使用动网格技术，模拟中充分考虑了灰渣形貌对流场的影响，具体影响在于流场变化引起的飞灰撞击质量流率的变化（详见3.1节）. 模拟中灰块形成位置和实验相同，均发生在迎风面^[12].

3种工况下形成的探针中心位置灰渣表面温度随时间变化如图14所示. 可以看出，随着沉积时间的增加，灰渣表面温度先迅速升高而后增速缓慢，该温度变化趋势和实验结果保持一致. 在数值模拟中，R=0、5%和10%工况在100 min时的灰渣表面温度分别为1 353、1 309和1 430 K，该结果和实验结果较为吻合. 通过探针中心位置灰渣的热流大小如图15所示. 可以看出，随着灰渣逐步形成，通过灰渣的热流密度逐渐减小，这是因为灰渣厚度的增加会引起热阻增加，从而导致换热特性有所下降，甚至导致传热恶化. 在0～100 min中，神华煤工况的热流密度降低率为52.13%，R=5%和10%工况下的热流密度降低率分别为46.96%和53.25%.

本文采用结合动网格技术的CFD模型对煤粉-玉米秸秆掺烧（玉米秸秆掺烧比例R为0、5%和10%）3个工况过程中的飞灰沉积进行了数值模拟，考虑了沉积过程中灰渣表面温度和灰渣导热系数的变化，并将数值模拟结果（灰渣增长高度、灰渣表面温度和热流密度等）和实验结果进行了比较. 模拟结果表明，飞灰撞击质量流率会随沉积进行而减小，这是因为沉积物形貌变化会影响探针周边的流场. 3种工况下100 min时的沉积效率分别为75%、80%和87%，这说明玉米秸秆掺烧的灰分更易于发生结渣. 在灰渣厚度方面，三者在0~100 min探针中心位置的灰渣增长率分别为0.055、0.062和0.069 mm/min，与实验数据吻合较好. 动网格执行下的灰渣形貌实时变化，会进一步影响飞灰沉积情况.

在换热特性方面，探针中心处灰渣表面温度随沉积时间的增加先迅速上升后缓慢上升，模拟得到的100 min时3种工况的灰渣表面温度分别为1 353、1 309和1 430 K；与之相反，热流密度逐渐减小，进而会导致传热特性下降，3种工况下0~100 min热流密度降低率分别为52.13%、46.96%和53.25%. 因此，在实际工业应用中需要考虑玉米秸秆掺烧引起的传热恶化问题.