煤是一种非常重要的化石燃料。在中国，因为其储量相对丰富、价格相对低廉而被广泛使用^[1]。当前，煤贡献了中国电力供应的70%以上，已成为中国能源结构中不可或缺的部分，而且这种状况还将长期维持^[2]。工业设备中煤的燃烧通常可以分为煤的热解^[3]与煤焦燃烧两个阶段，因为煤焦燃烧是非均相反应，需要的时间比煤的热解长几个数量级，同时煤燃烧过程中绝大部分热量是在煤焦燃烧过程中释放的^[4-5]，所以煤焦燃烧在很大程度上决定了煤燃烧的动力学特性。

目前，对单颗粒煤焦燃烧的模拟常使用一些简化的方法，如单膜模型和双膜模型等^[6-7]。单膜模型只考虑碳与O₂生成CO₂的反应，假设该反应只发生在煤焦颗粒的表面，气相中没有火焰面，最高温度点就在颗粒的表面，燃烧产生的CO₂从颗粒表面向外扩散，而O₂则从无穷远处向颗粒表面扩散，整个过程中只有颗粒表面一层反应“膜”，因此称之为单膜模型。双膜模型假设碳与O₂反应只生成CO，该反应也发生在颗粒的表面，生成的CO会向外扩散，在离颗粒表面一定距离的地方与O₂发生燃烧反应，形成一个火焰面，CO与O₂的燃烧反应也仅仅局限在这个面上，因而整个过程中形成了火焰面与颗粒表面两层反应“膜”，称之为双膜模型。这两种方法因为相对简单而得到较广泛应用，但它们仅仅是近似计算方法，存在很多问题： 1) 煤焦燃烧是个复杂的过程，涉及多个化学反应，而单、双膜模型只考虑了其中的1到2个反应; 2) 煤焦颗粒是多孔介质^[8-10]，气固反应并不仅仅发生在煤焦颗粒外表面，颗粒孔隙内也会有反应发生，CO与O₂的燃烧反应并不仅仅在一个反应面上发生，而是发生在CO从颗粒表面向外扩散的整个过程中^[11]; 3) 颗粒周围气体有复杂的扩散与流动^[12]，而单、双膜模型忽略了流动，对扩散也是简化处理; 4) 单、双膜模型对整个燃烧过程还作了准稳态的假设，而实际过程则是一个非稳态的。

本文从多组分混合反应系统守恒方程^[13]出发，构建了描述单颗粒煤焦周围气体流场的方程组，考虑了颗粒周围各种气体的流动与扩散作用，以CO与O₂的反应作为方程中的源项。煤焦颗粒的燃烧属于非均相反应，它对周围流场的影响通过边界条件引入，同时流场的变化也会及时反馈到煤焦颗粒的燃烧上。煤焦颗粒也有自己独立的质量与能量守恒方程。煤焦燃烧速率通过He等^[4]构建的分形煤焦燃烧速率模型来描述，该模型考虑了煤焦颗粒孔隙内的扩散与反应的耦合作用。最后，实验测得了4种煤焦在沉降炉中燃烧的最终转化率，并用本文构建的方法对此进行预测，通过与实验结果的比较来验证方法的可靠性。

1 数值模拟方法的构建

煤焦燃烧过程中主要发生以下4个反应^[14]：

$\begin{align} & \text{C}\left( \text{s} \right)+1/2{{\text{O}}_{2}}\left( \text{g} \right)\to \text{CO}\left( \text{g} \right) \\ & \Delta {{H}_{1}}=-111.5\text{kJ}\cdot \text{mo}{{\text{l}}^{-1}}, \\ \end{align}$

(1)

$\begin{align} & \text{C}\left( \text{s} \right)+{{\text{O}}_{2}}\left( \text{g} \right)\to \text{C}{{\text{O}}_{2}}\left( \text{g} \right) \\ & \Delta {{H}_{2}}=-393.8\text{kJ}\cdot \text{mo}{{\text{l}}^{-1}}, \\ \end{align}$

(2)

$\begin{align} & \text{C}\left( \text{s} \right)+\text{C}{{\text{O}}_{2}}\left( \text{g} \right)\to 2\text{CO}\left( \text{g} \right) \\ & \Delta {{H}_{3}}=+170.8\text{kJ}\cdot \text{mo}{{\text{l}}^{-1}}, \\ \end{align}$

(3)

$\begin{align} & \text{CO}\left( \text{g} \right)+1/2{{\text{O}}_{2}}\left( \text{g} \right)\to \text{C}{{\text{O}}_{2}}\left( \text{g} \right) \\ & \Delta {{H}_{4}}=-282.2\text{kJ}\cdot \text{mo}{{\text{l}}^{-1}}. \\ \end{align}$

(4)

其中： 反应(1)-(3)为气固反应，而反应(4)为气相反应; ΔH为反应热，kJ·mol^-3。

1.1 描述颗粒周围流场的方程组

为了简化分析与计算过程，假设煤焦颗粒为各向同性的均匀球体。由于煤焦颗粒的尺寸很小，这样做不会引入太大的误差。显然，采用球坐标对颗粒周围流场进行描述更方便，其坐标变量分别为r、 Φ和θ。其中: r≥0，0≤Φ＜2π，0≤θ≤π。同时，这还是一个1维球对称问题，因而有：

${{u}_{\phi }}={{u}_{\theta }}=0,$

(5)

$\partial \left( \centerdot \right)/\partial \phi =\partial \left( \centerdot \right)/\partial \theta =0.$

(6)

其中： u_Φ和u_θ分别为流场中Φ方向和θ方向的速度，m·s^-1。

利用式(5)和(6)对Kuo等^[13]的球坐标下的多组分混合反应系统守恒方程进行简化，可以得到描述煤焦颗粒周围流场的微分方程组。

混合物的连续性方程为

$\frac{\partial \rho }{\partial t}+\frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( \rho {{r}^{2}}{{u}_{r}} \right)=0.$

(7)

其中：ρ为混合气体的密度，kg·m^-3; u_r为流体径向速度，m·s^-1; r为径向坐标，m; t为时间，s。

每种组分的连续性方程为

$\frac{\partial {{Y}_{i}}}{\partial t}+{{u}_{r}}\frac{\partial {{Y}_{i}}}{\partial r}-\frac{1}{\rho {{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}\rho D\frac{\partial {{Y}_{i}}}{\partial r} \right)=\frac{{{\omega }_{i}}}{\rho }.$

(8)

其中： Y_i为组分i的质量分数; D为扩散系数，m²·s^-1; ω_i为组分i在气相反应中的生成速率，kg·m^-3·s^-1。这里，煤焦颗粒周围气体组分为O₂、 CO₂、 CO和N₂，因此式(8)中i的取值为O₂、 CO₂和CO，而N₂的质量分数可以用式(9)计算，

${{Y}_{{{\text{N}}_{\text{2}}}}}=1-{{Y}_{{{\text{O}}_{\text{2}}}}}-{{Y}_{\text{C}{{\text{O}}_{\text{2}}}}}-{{Y}_{\text{CO}}}.$

(9)

因为煤焦燃烧过程中的气相反应只有式(4)表征的CO燃烧反应，所以式(8)中的ω_i便对应于此反应导致的各种气体的净产生速率。这里，引用Howard等^[15]实验测得的CO燃烧速率表达式，

$\begin{align} & {{r}_{4}}=1.3\times {{10}^{8}}\left( \frac{\rho {{Y}_{\text{CO}}}}{{{M}_{\text{CO}}}} \right){{\left( \frac{\rho {{Y}_{{{\text{O}}_{2}}}}}{{{M}_{{{\text{O}}_{2}}}}} \right)}^{0.5}}. \\ & {{c}_{{{\text{H}}_{2}}\text{O}}}^{0.5}\text{exp}\left( -\frac{125604}{RT} \right). \\ \end{align}$

(10)

其中： r₄为CO的燃烧速率，mol·m^-3·s^-1; M_i为组分i的摩尔质量，kg·mol^-1; T为气体温度，K; R为普适气体常数，8.314 J·mol^-1·K^-1; c_H2O为水蒸气的饱合摩尔浓度，在101.325 kPa、 308.15 K下为2.5 mol·m^-3[14]。根据CO燃烧反应的配比关系可以得到：

${{\omega }_{{{\text{O}}_{2}}}}=-{{r}_{4}}{{M}_{{{\text{O}}_{2}}}}/2,$

(11)

${{\omega }_{\text{C}{{\text{O}}_{2}}}}={{r}_{4}}{{M}_{\text{C}{{\text{O}}_{2}}}},$

(12)

${{\omega }_{\text{CO}}}=-{{r}_{4}}{{M}_{\text{CO}}}.$

(13)

颗粒周围气体的动量守恒方程为

$\begin{align} & \frac{\partial {{u}_{r}}}{\partial t}+{{u}_{r}}\frac{\partial {{u}_{r}}}{\partial r}=-\frac{1}{\rho }\frac{\partial \rho }{\partial r}+\frac{1}{\rho {{r}^{2}}}\frac{\partial }{\partial r}\left( 2\mu {{r}^{2}}\frac{\partial {{u}_{r}}}{\partial r} \right)- \\ & \frac{1}{\rho {{r}^{2}}}\frac{\partial }{\partial r}\left[ \frac{2\mu }{3}\frac{\partial }{\partial r}\left( {{r}^{2}}{{u}_{r}} \right) \right]-\frac{4\mu {{u}_{r}}}{\rho {{r}^{2}}}+\frac{4}{3}\frac{\mu }{\rho {{r}^{3}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{u}_{r}} \right). \\ \end{align}$

(14)

其中： p为气体压力，Pa; μ为气体的动力粘性系数，N·s·m^-2。

颗粒周围气体的能量守恒方程为

$\begin{align} & \frac{\partial \left( {{c}_{V}}T \right)}{\partial t}+{{u}_{r}}\frac{\partial \left( {{c}_{V}}T \right)}{\partial r}=\frac{1}{\rho {{r}^{2}}}\frac{\partial }{\partial r}\left( \lambda {{r}^{2}}\frac{\partial T}{\partial r} \right)- \\ & \frac{1}{\rho }\sum\limits_{i=1}^{N}{{{\omega }_{i}}\Delta h_{\text{f,}i}^{\text{o}}+}\frac{1}{\rho {{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}\rho T\sum\limits_{i=1}^{N}{{{c}_{p,i}}D\frac{\partial {{Y}_{i}}}{\partial r}} \right)+ \\ & \frac{2\mu }{\rho }{{\left( \frac{\partial {{u}_{r}}}{\partial r} \right)}^{2}}+\frac{4\mu }{\rho }{{\left( \frac{{{u}_{r}}}{r} \right)}^{2}}-\frac{2\mu }{3\rho }{{\left[ \frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}{{u}_{r}} \right) \right]}^{2}}. \\ \end{align}$

(15)

其中: T为气体的温度，K; c_V为气体的定容比热容，J·kg^-1·K^-1; λ为气体的导热系数，W·m^-1·K^-1; Δh_f，i^o为组分i的标准生成焓，J·kg^-1; c_p，i为组分i的定压比热容，J·kg^-1·K^-1。

为封闭方程组，还必须引入气体的状态方程，

$p=\rho RT\sum\limits_{i={{\text{O}}_{2}},\text{C}{{\text{O}}_{2}},\text{CO,}{{\text{N}}_{2}}}{{{Y}_{i}}/{{M}_{i}}}.$

(16)

这样，方程数与变量数便相等，方程组可解。

1.2 描述颗粒状态的方程组

煤焦颗粒与O₂的反应速率方程为

${{r}_{\text{C+}{{\text{O}}_{2}}}}={{k}_{\text{a}}}m{{S}_{\text{g}}}{{c}_{{{\text{O}}_{2}}}}.$

(17)

其中： r_C+O2为煤焦颗粒与O₂的反应速率，mol·s^-1; k_a为煤焦燃烧的表观反应系数; m为煤焦颗粒的质量，kg; S_g为煤焦颗粒的比表面积，m²·kg^-1; c_O2为煤焦颗粒表面O₂的摩尔浓度，mol·m^-3。

He等^[4]在其建立的煤焦燃烧模型中给出了k_a的表达式，

${{k}_{\text{a}}}=3{{k}_{\text{i}}}\left( \frac{\coth Q}{Q}-\frac{1}{{{Q}^{2}}} \right).$

其中：

$\begin{align} & {{k}_{i}}={{A}_{\text{i}}}\exp \left( -{{E}_{\text{i}}}/R{{T}_{\text{p}}} \right), \\ & Q=\frac{1}{2}{{\left[ \frac{{{k}_{i}}}{\varsigma \chi {{\left( {{T}_{\text{p}}}/0.032 \right)}^{1/2}}} \right]}^{1/2}}, \\ & \chi =\frac{2{{\theta }^{5/3}}}{S_{\text{g}}^{2}\rho _{\text{i}}^{2}d_{\text{p}}^{2}{{\left( 1-\theta \right)}^{2}}}\exp \left( -{{D}_{\text{f}}} \right), \\ & {{A}_{\text{i}}}=62.37,{{E}_{\text{i}}}=54000\text{J}\cdot \text{mo}{{\text{l}}^{-1}} \\ & \varsigma =1.975. \\ \end{align}$

(18)

其中： k_i为煤焦燃烧的本征反应系数; Q为中间变量; T_p为煤焦颗粒的温度，K; θ为孔隙率; D_f为煤焦孔隙的分形维数^[16-17]; d_p为颗粒直径，m; ρ_i为煤焦的真密度，kg·m^-3; χ为模数，由煤焦颗粒孔隙参数组合而成，表征了煤焦孔隙结构对气体扩散的阻力。该模型考虑了煤焦颗粒内部扩散与反应的耦合作用。

对于燃烧产物中CO/CO₂比例，引用Tognotti等^[18]的实验结果，

$\eta =50{{\left( {{p}_{{{\text{O}}_{2}}}}/101325 \right)}^{-0.21}}\exp \left( -3040/{{T}_{\text{p}}} \right).$

(19)

其中： η为燃烧产物中CO/CO₂比例; p_O2为煤焦颗粒外表面O₂的分压，Pa。

${{r}_{1}}={{r}_{\text{C+}{{\text{O}}_{2}}}}\cdot \eta /\left( \eta +1 \right),$

(20)

${{r}_{2}}={{r}_{\text{C+}{{\text{O}}_{2}}}}/\left( \eta +1 \right).$

(21)

其中： r₁和r₂分别为式(1)和(2)中化学反应的速率，mol·s^-1。

对于煤焦与CO₂的反应速率引用Annamalai和Ryan^[19]的实验结果，

$\begin{align} & {{r}_{3}}=0.034{{\left( {{p}_{\text{C}{{\text{O}}_{2}}}}/101325 \right)}^{0.34}}m{{S}_{\text{g}}}\cdot \\ & \exp \left( -\frac{109000}{R{{T}_{\text{p}}}} \right)/{{M}_{\text{C}}}. \\ \end{align}$

(22)

其中： r₃为式(3)中化学反应的速率，mol·s^-1; p_CO2为煤焦颗粒外表面CO₂的分压，Pa; M_C为碳的摩尔质量，kg·mol^-1。

煤焦颗粒的能量守恒方程为

$\begin{align} & \frac{m{{c}_{\text{char}}}}{4\pi r_{\text{p}}^{2}}\frac{\text{d}{{T}_{\text{p}}}}{\text{d}t}+\varepsilon \sigma \left( T_{\text{p}}^{4}-T_{\text{w}}^{4} \right)+h\left( {{T}_{\text{p}}}-{{T}_{\infty }} \right)= \\ & -\frac{{{r}_{1}}\Delta {{H}_{1}}+{{r}_{2}}\Delta {{H}_{2}}+{{r}_{3}}\Delta {{H}_{3}}}{4\pi r_{\text{p}}^{2}}. \\ \end{align}$

(23)

其中： c_char为煤焦颗粒的比热容，J·kg^-1·K^-1; r_p为煤焦颗粒的半径，m; T_w周围环境的温度， K; T_∞为煤焦颗粒周围流体的平均温度，K，取为(T_p+T_w)/2; h为煤焦颗粒的对流换热系数，W·m^-2·K^-1，根据大空间中静止煤焦颗粒的 Nusselt 数等于2可以求得h=λ/r_p; ΔH₁、 ΔH₂和ΔH₃分别为式(1)-(3)中化学反应的反应热，J·mol^-1。式(23)的物理意义是煤焦颗粒燃烧所产生的热量一部分被煤焦颗粒自身吸收，用于提高颗粒的温度，剩下的热量则分别以辐射和对流方式向周围传递。

煤焦颗粒的质量变化为

$\frac{\text{d}\alpha }{\text{d}t}=\frac{\left( {{r}_{1}}+{{r}_{2}}+{{r}_{3}} \right){{M}_{\text{C}}}}{{{m}^{0}}\left( 1-\text{as}{{\text{h}}^{0}} \right)}.$

(24)

其中： α为煤焦颗粒的转化率; m⁰为煤焦颗粒的初始质量，kg; ash⁰为煤焦颗粒的初始灰分质量分数。

1.3 边界条件

定义流场中最靠近煤焦颗粒的一层气体网格为边界网格，与该网格相关的各个物理量均有下标0。由于煤焦颗粒处于流场之中，为满足温度场的连续性，这里采用第一类热边界条件，即边界网格中气体的温度与煤焦颗粒的温度相等，

${{T}_{0}}={{T}_{\text{p}}}.$

(25)

因为有气固反应的发生，会有气体进出煤焦颗粒表面，所以对于流场边界网格的速度，则不能简单地设置为0。首先，根据气固反应的速率计算出进出煤焦颗粒的各种气体的通量(mol·s^-1)：

${{N}_{{{\text{O}}_{2}}}}=-{{r}_{1}}/2-{{r}_{2}},$

(26)

${{N}_{\text{C}{{\text{O}}_{2}}}}={{r}_{2}}-{{r}_{3}},$

(27)

${{N}_{\text{CO}}}={{r}_{1}}+2{{r}_{3}}.$

(28)

流出颗粒的气体其通量为正，流入颗粒的气体其通量为负。根据状态方程可以将各种气体的通量折算成体积流量，再除以煤焦颗粒表面积，便可以得到各种气体组分的流速。最后，将各种气体的流速以质量分数加权求和，便可求得流场边界网格中气体的流速(m·s^-1)，

${{u}_{r0}}=\frac{R{{T}_{\text{p}}}\left( \left| {{N}_{{{\text{O}}_{2}}}} \right|{{N}_{{{\text{O}}_{2}}}}{{M}_{{{\text{O}}_{2}}}}+\left| {{N}_{\text{C}{{\text{O}}_{2}}}} \right|{{N}_{\text{C}{{\text{O}}_{2}}}}{{M}_{\text{C}{{\text{O}}_{2}}}}+\left| {{N}_{\text{CO}}} \right|{{N}_{\text{CO}}}{{M}_{\text{CO}}} \right)}{4\pi r_{\text{p}}^{2}{{p}_{0}}\left( \left| {{N}_{{{\text{O}}_{2}}}} \right|{{M}_{{{\text{O}}_{2}}}}+\left| {{N}_{\text{C}{{\text{O}}_{2}}}} \right|{{M}_{\text{C}{{\text{O}}_{2}}}}+\left| {{N}_{\text{CO}}} \right|{{M}_{\text{CO}}} \right)}.$

(29)

其中，p₀为边界网格中气体的压力，Pa。

由于煤焦燃烧反应不断将固体碳转化为CO和CO₂而带入气相中，因而边界网格中气体密度会有一个增量，

$\Delta {{\rho }_{0}}=\left( {{r}_{1}}+{{r}_{2}}+{{r}_{3}} \right){{M}_{\text{C}}}\text{d}t/{{V}_{0}}.$

(30)

其中： dt为计算的时间步长，s; V₀为边界网格的体积，m³。

根据式(26)-(28)中给出的各种气体进出煤焦颗粒的通量，也能计算出边界网格中各种气体组分的质量分数变化：

$\Delta {{Y}_{{{\text{O}}_{2}},0}}={{N}_{{{\text{O}}_{2}}}}{{M}_{{{\text{O}}_{2}}}}\text{d}t/\left( {{\rho }_{0}}{{V}_{0}} \right),$

(31)

$\Delta {{Y}_{\text{C}{{\text{O}}_{2}},0}}={{N}_{\text{C}{{\text{O}}_{2}}}}{{M}_{\text{C}{{\text{O}}_{2}}}}\text{d}t/\left( {{\rho }_{0}}{{V}_{0}} \right),$

(32)

$\Delta {{Y}_{\text{CO,0}}}={{N}_{\text{CO}}}{{M}_{\text{CO}}}\text{d}t/\left( {{\rho }_{0}}{{V}_{0}} \right).$

(33)

其中，ρ₀为边界网格中气体的密度，kg·m^-3。

2 实验验证

将4种煤焦在沉降炉中进行燃烧实验，分别测得它们最终的转化率，并运用本文中构建的模拟方法对这些最终转化进行预测，通过与实验结果的对比来验证模拟方法的可靠性。

实验中所用沉降炉的结构如图 1所示，其炉管内径为50 mm，炉管高为750 mm，沉降炉注入煤焦颗粒的位置距离炉管进口200 mm。该沉降炉为电加热，取样点设在炉管的出口处。

实验样品煤焦所对应的原煤的工业分析与元素分析分别列于表 1和2中。将这些原煤在沉降炉中1 200 ℃高温N₂气氛下进行热解，以制得燃烧实验用的煤焦样品。

作煤焦燃烧实验时，沉降炉的壁温为1 000 ℃，炉内所通气体为空气，流量为0.5 N·m³·h^-1，流速为0.33 m·s^-1，煤焦颗粒的给粉速率为1 g·h^-1，煤焦颗粒在炉管内的停留时间为1.67 s。

沉降炉给粉速率足够小，以保证煤焦颗粒之间互不影响，这样本文构建的模拟方法可以近似适用。实验中，煤焦颗粒离开炉管后，使用热重基于灰分守恒法测量出煤焦的转化率。在模拟时，无需模拟煤焦颗粒的整个燃烧过程，只需计算到实际的炉内停留时间(1.67 s)即可，此时的煤焦转化可以与实验测量值进行对比。

模拟煤焦燃烧还需要知道煤焦样品的平均粒径与孔隙结构参数，它们分别由Mastersizer 2000激光粒径分析仪和Autoscan 33压汞仪测得，测量结果列于表 3中。可以看到，所有煤焦样品的平均粒径都很小，最大的也不到120 μm，这样煤焦颗粒在沉降炉内是跟随气流一起运动的，颗粒与气流之间的相对速度近似为零，因此本文所提出的方法可以用来近似模拟沉降炉实验。但是，当煤焦颗粒尺寸较大时，则这种方法的误差会明显增大，而本文不存在这个问题。

实验测得的煤焦最终转化率见表 4，同时表 4也给出了模拟预测的结果。可以看出，实测结果与预测结果吻合得较好，从而证明了本文所提出的模拟单颗粒煤焦燃烧的数值方法具有一定的可靠性。还可看出，不同煤焦的最终转化率相差很大，这与煤焦本身的属性(孔隙参数与粒径)及选用的燃烧速率模型有关。然而，这些煤焦预测结果都没有达到相对误差小于5%，如表 4所示。造成误差的因素有很多，如本文使用的1维假设、求解偏微分方程组时引入的数值误差、所用煤焦孔内燃烧模型存在缺陷以及实验中煤焦样品的不均匀性等。虽然如此，目前的误差对于工程应该来说还是可以接受的。在将来的工作中将会引入算法的敏感性分析，以进一步提高算法的精度。

本文还分别采用单膜模型与双膜模型对沉降炉的实验进行模拟，两种模型预测的煤焦最终转化率及相对误差列于表 5中。与表 4中的数据对比可以看出，本文构建的模型精度最高，其次是双膜模型，单膜模型的精度最低。这也说明，相比于单膜模型和双膜模型，本文构建的模型不仅在理论上具有优势，在实际预测效果上也是最好的。

3 对煤焦燃烧过程的模拟

本节将对一种煤焦在不同环境温度下的燃烧过程进行模拟，并给出燃烧过程中颗粒温度、煤焦转化率、总反应速率以及颗粒边界层中各种气体组分摩尔分数的变化规律，以加深对煤焦燃烧过程的了解。

该煤焦颗粒的平均粒径为30.0 μm，孔隙率为48.0%，比表面积为2.8 m²·g^-1，分形维数为 1.06。分别模拟该煤焦颗粒在空气气氛中外界环境温度分别为1 000、 1 200、 1 400和1 600 K条件下的燃烧过程。

图 2给出了燃烧过程中煤焦颗粒温度的变化规律。煤焦颗粒初始温度为室温300 K，由于环境温度较高，同时煤焦颗粒尺寸非常小，因而煤焦颗粒的升温速率很快。当煤焦颗粒温度足够高时，燃烧反应开始，这个过程会放出大量的热，导致煤焦颗粒温度高出环境温度。最后，由于反应速率减慢，煤焦颗粒通过对流与辐射的散热量大于燃烧的产热量，煤焦颗粒温度会缓慢下降。这个变化规律与He等^[20]的研究结果相符。

图 3为燃烧过种中煤焦颗粒的转化率曲线。可以看到，环境温度越高，煤焦颗粒燃烧越快，燃尽时间也就越短。

图 4为燃烧过程中煤焦颗粒总反应速率的变化规律。图 4中，r=r₁+r₂+r₃，不仅包括碳与O₂的反应，还包括碳与CO₂的反应。可以看到，燃烧过程中反应速率先增大后减小，有一个最大值出现。反应速率先上升是由煤焦颗粒的温度不断上升导致的，而后的下降则是由反应中后期煤焦的有效反应面积随着反应进行不断减小导致的，而此时，燃烧产热与颗粒散热达到平衡，颗粒温度变化较小。图 4还表明，煤焦颗粒的最大燃烧速率随着环境温度的升高而变大。

图 5-7分别为燃烧过程中煤焦颗粒边界层中O₂、 CO₂和CO摩尔分数的变化规律。可以看到，随着环境温度的升高，煤焦颗粒表面O₂的摩尔分数下降越多，而CO₂和CO的摩尔分数则上升越多。这是因为环境温度升高会使燃烧速率加快，从而增加了对O₂的需求量，这就使扩散与反应达到平衡时O₂的摩尔分数降得更低，进而作为燃烧产物的CO₂和CO的摩尔分数则会变得更高。到燃烧后期，反应速率减慢，煤焦颗粒边界层中各种组分的摩尔分数由于扩散的原因又有还原的趋势，即O₂的摩尔分数逐渐上升，而CO₂和CO的摩尔分数逐渐下降。

对比图 6与7还可以发现，煤焦颗粒表面CO的摩尔分数要比相同反应条件下CO₂的摩尔分数高，在所有温度条件下均是如此，这说明，CO是煤焦燃烧的主要产物，这也是双膜模型^[7]中假设的碳与O₂燃烧只生成CO的一个原因。随着环境温度的升高，CO与CO₂摩尔分数的差别会变小，这主要是由于高温使CO与O₂反应的速率加快，促进了产物中CO向CO₂的转化。

4 结论

本文从多组分混合反应系统守恒方程出发，构建了一种能较为准确地模拟单颗粒煤焦燃烧过程的数值方法。该方法考虑了煤焦燃烧过程中气固反应与气相反应，不仅能预测煤焦颗粒的温度、转化率、燃烧速率等与煤焦相关的量，还能较好地描述煤焦颗粒周围气体的流动、扩散及传热。

通过对4种煤焦在沉降炉中燃烧的最终转化率进行预测，并与实验测得的数据进行对比，证明了本文方法的可靠性，也说明该方法可以用于对一些流动形态较为简单的实际燃烧过程的模拟。

本文通过对煤焦颗粒燃烧过程的模拟，发现煤焦燃烧速率在整个燃烧过程中先增大后减小。随着环境温度的升高，煤焦燃烧速率加快，颗粒表面O₂的摩尔分数下降更多，而CO₂与CO的摩尔分数则升高更多。从产物的量来看，CO是煤焦燃烧的主要产物。