气固反应可分为气固催化反应与气固非催化反应，前者主要是气相发生化学变化，固相主要起催化作用。而气固非催化反应，随着反应的进行，固体反应物逐渐消耗，并在固体反应物表面生成固体产物层，这一类反应简称为气固反应。气固反应在工业生产中占有很重要的地位，如铁的氧化，碳的燃烧与气化，氧化钙、氧化镁等与二氧化硫的反应，金属氧化物脱除H₂S的反应，金属氧化物的还原等。近年来，随着人们对清洁能源的高效利用以及温室气体CO₂减排的关注，研究者又将气固反应应用至CO₂捕集、CO₂转化利用、制氢、化工品合成以及纳米燃烧等方面，如利用固体吸附剂脱除CO₂^[1]、热化学储热^[2]、化学链燃烧^[3]、CO₂的还原^[4]、烷烃脱氢制烯烃^[5]、铁-蒸汽法制氢^[6]、煤或生物质的气化^[7]、天然气制合成气^[8]等。同时，还可以利用气固反应进行金属纳米颗粒燃烧^[9]、太阳能制氢^[10]等。气固反应在大气颗粒物及其环境效应^[11]、地球化学^[12]等领域也有着广泛应用。

虽然气固反应体系由于反应对象的不同而发生不同类型的反应，但所涉及的微观物理/化学步骤是相似的。固体反应物通常是多孔颗粒，气体分子可通过孔隙扩散到颗粒内部并发生反应，如图 1所示，气固反应步骤包括^[13-14]：1) 气体分子从气相主体扩散到颗粒外表面，即外扩散；2) 气体分子进一步通过颗粒内部的孔隙进行扩散，即内扩散；3) 气体分子到达局部固体反应物表面后发生化学反应；4) 反应物分子/离子通过产物层扩散后进一步反应；5) 固体产物逐渐生长并导致颗粒结构变化。上述气固反应步骤同时进行，并相互影响制约，是典型的多尺度、多物理/化学步骤耦合系统。化学反应发生在表界面上属于原子/电子尺度；而反应物分子/离子通过产物层的体相扩散发生在纳米或微米晶粒尺度上；气体分子的外扩散则发生在微米或毫米的颗粒尺度上。气固反应对象跨越原子、微米与毫米尺度，需采用不同的数学物理模型。同时气固反应是动态变化过程，固体产物在颗粒内各个位置处不断生长，固体产物与固体反应物的摩尔体积存在差别，故固体产物的生长将导致颗粒孔隙结构尺寸逐渐变化^[14]，进而影响气体内扩散阻力。因此，反应过程中固体产物的生长将改变反应控制步骤^[14]。例如，初始反应时，颗粒内部孔隙尺寸较大，气体传质阻力较小，此时反应界面处的表面化学反应为控制步骤；若该反应过程中生成的固体产物体积大于消耗的固体反应物体积，则颗粒内固体产物取代固体反应物后将导致颗粒内孔隙尺寸减小^[14]，甚至出现颗粒孔隙堵塞现象，最终气体内扩散取代表面化学反应成为控制步骤。跨越原子尺度到宏观颗粒尺度的气固反应动力学是能源环境、化工冶金、材料等领域的学科前沿，在碳中和场景下的能源清洁转化与利用、碳捕集转化与利用、热化学储能、污染物脱除等国家重大需求领域有重要应用价值。

气固反应动力学的研究最早可以追溯到1920年，当时Tammann^[15]在银、铅与镉氧化时，观察到氧化膜增厚与氧化时间平方根成正比。1923年，Pilling与Bedworth^[16]确认抛物线氧化速度定律具有普遍性，提出了氧化膜增长速度的经验方程。金属氧化中产物层生长的经典理论是由Wagner^[17]于1933年提出，以化学位梯度和电位梯度为驱动力，建立了高温氧化物产物层生长的经典理论模型，该理论不适用于产物层比较薄的情况。1949年英国的Mott和Cabrera^[18]借鉴量子力学隧道效应提出了低温氧化的薄膜生长理论，当产物层厚度超过100时，薄膜生长理论不再适用。Fromhold^[19]于1976年对Mott-Cabrera低温氧化理论进行了修正。一些学者也对产物层的扩散机理^[20-21]进行了相关研究，主要针对所谓的产物层扩散究竟是气体扩散，还是离子扩散等方面来进行研究。无论是Mott-Cabrera理论，还是Wagner理论，以及其他相关研究，均是以固体反应物表面上已经存在一定厚度的产物层为前提的，即通常所说的产物层扩散阶段，研究集中在扩散的机制与扩散过程的宏观描述上。关于气固反应的颗粒模型研究，在早期的研究中，人们常将固体反应物颗粒认为是致密无孔的^[22]，即不考虑其结构特征。大多数情况下，固体反应物颗粒是多孔物质，孔隙率、比表面及核孔径分布将显著影响反应进程，近年来，人们借助多相催化理论的发展，对孔扩散和反应以及多孔固体结构特性之间的耦合加以考虑，提出了相关的晶粒模型^[23]、孔隙模型^[24]和逾渗模型^[25]等。这些模型主要是研究气体在颗粒内扩散与多孔结构之间的耦合特性，采用表观的总包反应来处理化学过程，采用拟合实验数据的方法来获得反应动力学参数^[25]，而基本不涉及表界面上原子/电子尺度的化学反应路径与机理，更没有关注跨尺度关联及多物理/化学步骤耦合这一基本问题。

在气固催化反应研究领域，量子化学计算可以预测分析原子/电子尺度的化学反应路径与机理^[26]，表界面尺度的反应现象可以通过分子动力学方法进行计算分析^[27]。气固催化反应所研究的主要是材料的表界面化学/物理过程，材料的体相近似按惰性物质处理^[28]。而在如图 1所示的气固非催化反应体系中，不仅表界面的原子参与反应，而且体相的物质也通过固态扩散的方式与气体分子接触参与反应并生成固体产物，同时伴随着动态的结构演化，体系的特征尺度为10²~10⁴ nm；同时，上述反应过程与设备内的流动、传热、传质相互耦合，增加了系统复杂性。目前主要采用经验或半经验性的表观模型来拟合实验数据以建立不同物理量间的关联，表观方法无法刻画反应机理与反应本质^[14]；现有的量子化学计算和分子动力学模拟局限于微观尺度，难以应用于跨尺度相互作用的描述，构成固体颗粒的大量原子/分子所涌现的集体的气固反应行为无法从底层的量子化学计算和分子动力学理论中推导出来。现有理论对描述气固反应动力学的跨尺度现象与机理面临根本性困难，迫切需要新的理论方法。

气固反应跨越了“原子→晶粒→颗粒”多重尺度，新的理论方法必须考虑在不同尺度上的关联。近年来，将固体颗粒视为由大量晶粒构成的多孔系统，通过研究晶粒反应特性进而得到颗粒整体行为的方法备受关注^[14]，为刻画气固反应动力学的跨尺度关联和反应/扩散过程耦合提供了新思路。基于晶粒模型，针对气固反应中的跨尺度关联和反应/扩散耦合这一问题，本文从介观晶粒层次提出基于第一性原理的速率方程理论，构建“原子→晶粒→颗粒”跨尺度的表界面气固反应动力学与体相扩散耦合模型，揭示了固体产物在表面的离散岛状生长机制，设计高效可靠的计算和分析方法，并应用于低能耗CO₂捕集、热化学储热、污染物脱除等技术，有力推动了该领域的发展。

1 晶粒模型

气固反应中的固体颗粒通常是多孔的，将多孔颗粒视为由大量晶粒构成^[14]，如图 2所示，晶粒与晶粒之间存在一定的孔隙。这样就把宏观上的颗粒离散为相互独立的晶粒，通过研究晶粒反应特性进而得到颗粒的整体行为^[14]。

晶粒的初始半径r₀ (m)可通过测得到的固体颗粒比表面积S₀ (m²/m³)和孔隙率ε₀计算得到

$ r_{0}=\frac{3\left(1-\varepsilon_{0}\right)}{S_{0}}. $

(1)

以H₂与Fe₂O₃的反应为例，当气体分子H₂扩散到Fe₂O₃晶粒表面时发生吸附、表面反应以及产物H₂O的脱附等一系列表面反应步骤。表面反应过程中Fe₂O₃会被消耗而生成固体产物。当反应进行到一定时间后，固体产物将完全覆盖整个固体反应物Fe₂O₃的表面，形成产物层。固体产物层将固体反应物Fe₂O₃和气体分子H₂分隔开，阻碍了固体反应物Fe₂O₃活性表面与气体分子H₂的直接接触，此时反应进程进入到产物层扩散阶段。在产物层扩散阶段，反应物通过固体产物层扩散到反应界面处才能进一步发生化学反应。产物层扩散有两种模式，包括气体产物层扩散和离子产物层扩散。对于气体产物层扩散模式，气体分子通过产物层从固体产物/气体的交界面扩散到固体反应物/固体产物的交界面，然后与固体反应物发生化学反应，该反应界面在固体反应物/固体产物的交界面处；对于离子产物层扩散模式，由于产物层内部的氧离子浓度C_O* (mol/m³)远高于产物层外部的氧离子浓度C_{O*, shell}(mol/m³)，产物层内外存在氧离子浓度梯度，固体反应物的氧离子通过产物层从固体反应物/固体产物的交界面扩散到固体产物/气体的交界面，然后与气体分子发生化学反应，如图 3b所示，该反应界面在固体产物/气体的交界面。随着产物层扩散和反应界面处化学反应的进行，固体反应物和气体分子不断被消耗，同时固体产物层的厚度逐渐增加^[14]。

晶粒模型认为化学反应、产物层生长及产物层扩散均发生在晶粒表面。反应过程中晶粒的结构变化过程如图 3所示，r₀为晶粒的初始半径(m)，r₂为晶粒的未反应核半径(m)，r₁为整个晶粒的半径(m)。固体未反应核在反应过程中逐渐被消耗且尺寸逐渐收缩，即r₂在反应中逐渐减小。根据图 3的晶粒几何结构，晶粒的转化率X可表示为

$ X=1-\frac{r_{2}^{3}}{r_{0}^{3}}. $

(2)

则晶粒转化速率为

$ \frac{\mathrm{d} X}{\mathrm{~d} t}=-\frac{3 r_{2}^{2}}{r_{0}^{3}} \frac{\mathrm{d} r_{2}}{\mathrm{~d} t}. $

(3)

由于未反应核尺寸在反应过程中整体逐渐收缩，产物层覆盖在未反应核的外表面，故存在几何关系式

$ Z=\frac{r_{1}^{3}-r_{2}^{3}}{r_{0}^{3}-r_{2}^{3}}. $

(4)

式中Z是固体产物与固体反应物的摩尔体积之比。

2 表面反应机理的密度泛函(DFT)计算

晶粒表界面反应是整个气固反应过程的关键步骤，直接决定并影响气体反应物转化和气体产物选择性、固体反应物转化率、气固反应动力学等。随着量子力学和计算技术的发展，实验上无法获得或很难获得的微观反应机理与特性可以通过量子化学密度泛函理论(DFT)进行研究^[29]。DFT量子化学计算方法考虑原子核外电子间的相互作用^[30]，通过求解薛定谔方程，可以计算出固体表面基元反应的所有反应细节，即从表面基元反应的初态(initial state, IS)经过渡态(transition state, TS)搜索到达反应的末态(final state, FS)。DFT量化计算可以精确地预测物质结构和电子性质，能够揭示晶粒表界面反应的化学反应机理，被广泛应用于气固催化反应研究领域，并且已经取得了大量研究成果^[31-33]。以Fe₂O₃与H₂和CO的反应为例，使用VASP (vienna ab-initio simulation package)软件包^[34]对建立的Fe₂O₃表面结构模型探究其分别与H₂和CO反应所形成的可能路径，经过对不同反应路径的能量分析，进而获得最有可能发生的实际反应路径与机制。对于Fe₂O₃与H₂的反应，DFT计算发现H₂分子在表面发生分解反应，与O结合生成OH，其中位于O顶位的一个H从其O顶位迁移到另外一个O顶位H，并与其成键结合形成吸附态的H₂O分子，然后H₂O分子脱附进入气相，如图 4所示。对于Fe₂O₃与CO的反应，DFT计算发现CO分子首先吸附在表面，然后吸附态CO分子在表面与表面的活性O结合生成吸附态的CO₂分子，最后CO₂分子从固体表面脱附进入气相，如图 5所示。

在如图 4和图 5所示的反应路径基础上，可以写出Fe₂O₃与H₂和CO的表面反应方程式，如R1~R7所示(R表示Reaction，即化学反应方程式)，

$ \mathrm{H}_{2}(\mathrm{g}) {\underset{k_{-1}}{\overset{k_{1}}{\mathop{\rightleftharpoons }}}\,} \mathrm{H}_{2, \text { ad }}, $

(R1)

$ \mathrm{H}_{2, \text { ad}}+2 \mathrm{O}^{*} {\underset{k_{-2}}{\overset{k_{2}}{\mathop{\rightleftharpoons }}}\,}2 \mathrm{OH}, $

(R2)

$ \mathrm{OH}+\mathrm{OH} {\underset{k_{-3}}{\overset{k_{3}}{\mathop{\rightleftharpoons }}}\,} \mathrm{O}_{*}+\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}, $

(R3)

$ \mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}} {\underset{k_{-4}}{\overset{k_{4}}{\mathop{\rightleftharpoons }}}\,} \mathrm{H}_{2} \mathrm{O}(\mathrm{g}), $

(R4)

$ \mathrm{CO}(\mathrm{g}) {\underset{k_{-5}}{\overset{k_{5}}{\mathop{\rightleftharpoons }}}\,} \mathrm{CO}_{\mathrm{ad}}, $

(R5)

$ \mathrm{CO}_{\mathrm{ad}}+\mathrm{O}_{*} {\underset{k_{-6}}{\overset{k_{6}}{\mathop{\rightleftharpoons }}}\,} \mathrm{CO}_{2, \mathrm{ad}}, $

(R6)

$ \mathrm{CO}_{2, \mathrm{ad}} {\underset{k_{-7}}{\overset{k_{7}}{\mathop{\rightleftharpoons }}}\,} \mathrm{CO}_{2}(\mathrm{g}) . $

(R7)

DFT计算还可以得到每个基元反应正向反应和逆向反应的能垒(E_{A, fwd}和E_{A, bwd}，eV)。正向反应的能垒等于过渡态(TS)能量E_TS(eV)减去初态(IS)能量E_R(eV)，另外还需要校正零点能E_{ZP, fwd}(eV)，如式(5)所示。逆向反应的能垒等于过渡态(TS)能量E_TS减去末态(FS)能量E_P(eV)，另外还需要校正零点能E_{ZP, bwd}(eV)，如式(6)所示。

$ E_{\mathrm{A}, \text { fwd }}=E_{\mathrm{TS}}-E_{\mathrm{R}}+E_{\mathrm{ZP}, \text { fwd }}, $

(5)

$ E_{\mathrm{A}, \text { bwd }}=E_{\mathrm{TS}}-E_{\mathrm{P}}+E_{\mathrm{ZP}, \text { bwd }}. $

(6)

因为DFT计算的是0 K时的表面反应，且所建立的固体表面结构是理想晶体表面；而实际的固体表面有很多缺陷且表面的不同位置的组成和形貌会存在较大差别，所以式(5)和(6)计算得到的反应能垒通常与实际值存在误差，误差一般为0~0.2 eV^[35]。根据DFT计算得到的反应路径及反应能垒可以对反应机理、材料活性、电子结构等进行一些定性分析^[36]，但是，DFT计算还不能解决气固反应动力学的预测以及速率控制步骤的具体分析，也无法直接反映温度压力等对反应过程的影响^[37]。

3 反应动力学参数的过渡态理论(TST)计算

从DFT计算结果预测操作条件下的反应动力学以及速率控制步骤等实际信息需要计算每个基元反应的反应动力学速率常数，采用过渡态理论(TST)可以实现从DFT数据到速率常数的计算。TST是1935年由Eyring^[38]、Polany和Evans^[39]等提出，建立在量子力学和统计力学的基础上。在图 4和图 5中，基于DFT计算的势能面说明由反应物分子(IS)转变为产物分子的过程中，反应体系一定要经过一个能级较高的过渡态(TS)，过渡态所对应的组分称为活化络合物，然后活化络合物再分解为产物(FS)。活化络合物的存在时间很短，一经生成，就快速向产物分子转化。过渡态理论认为^[38-39]：1) 反应物(IS)与活化络合物(TS)处于化学平衡状态；2) 活化络合物通过能垒的速率等于其分解速率乘以浓度。化学平衡常数可以采用统计力学中的配分函数计算，然后根据平衡常数计算出活化络合物的浓度，由此建立反应速率常数的计算式。

对于Langmuir-Hinshelwood类型的表面反应(X^*+Y^*=Z^*)^[40]，正逆反应速率的计算式分别为：

$ k_{i}=\frac{Q_{\mathrm{vib}}^{\neq}}{Q_{\mathrm{vib}}^{X^{*}} Q_{\mathrm{vib}}^{Y^{*}}} \frac{k_{\mathrm{B}} T}{h} \exp \left(-\frac{E_{\mathrm{A}, \text { fwd }}}{k_{\mathrm{B}} T}\right), $

(7)

$ k_{-i}=\frac{Q_{\mathrm{vib}}^{\neq}}{Q_{\mathrm{vib}}^{Z^{*}}} \frac{k_{\mathrm{B}} T}{h} \exp \left(-\frac{E_{\mathrm{A}, \mathrm{bwd}}}{k_{\mathrm{B}} T}\right) . $

(8)

对于表面吸附(X_(gas)=X^*)，吸附速率和脱附速率的计算式分别为：

$ k_{i}=\frac{p_{X_{(\mathrm{gas})}} A_{\mathrm{st}}}{\sqrt{2 {\rm{ \mathsf{ π} }} m_{X_{(\mathrm{gas})}} k_{\mathrm{B}} T}} \text {, } $

(9)

$ k_{-i}=\frac{Q_{\text {vib, } X_{(\text {gas})}}^{\text {gas }} Q_{\text {rot, } X_{(\mathrm{gas})}}^{\text {gas}} Q_{\text {trans2D, } x_{(\text {gas })}}^{\text {gas }}}{Q_{\text {vib }, X}^{X^{*}}}\frac{k_{\mathrm{B}} T}{h} \exp \left(-\frac{E_{\mathrm{des}}}{k_{\mathrm{B}} T}\right). $

(10)

式中，Q_vib^≠、Q_vib^X*、Q_vib^Y*、Q_{vib, X}^gas_(gas)、Q_{rot, X}^gas_(gas)、Q_{trans2D, X}^gas_(gas)和Q_{vib, X}^X*是配分函数，具体含义以及计算方法参见文^[14]；kB为Boltzmann常数，h为Planck常数，T为温度(K)，E_des为脱附能垒(eV)。TST采用理论计算的方式，由分子的振动频率、转动惯量、质量、反应的能垒等基本参数，分别计算平动、转动、振动配分函数^[14]，代入式(7)~(10)就能得到每个基元反应的正逆过程的速率常数k_i和k_-i。

4 气固反应动力学速率方程理论

基于DFT计算数据，通过过渡态理论计算出反应速率常数后，可以利用质量作用定律计算每个基元反应的反应速率。质量作用定律是由Guldberg和Waage于1867年提出^[41]，表述为基元反应的反应速率与各反应物浓度的幂的乘积成正比。但固体表面上的反应与表面活性位以及表面组分所占活性位的份额有关，质量作用定律很难描述表面反应特性。Langmuir^[42]于1916年提出表面覆盖率θ概念并推导得到Langmuir等温吸附方程，这样就可以基于表面质量作用定律列出每个表面基元反应的反应速率表达式。对于Fe₂O₃与H₂和CO的反应，反应速率表达式如下：

$ \dot{r}_{1}=k_{1} p_{\mathrm{H}_{2}} \theta_{\mathrm{O}_{*}}^{n_{1}}-k_{-1} \theta_{\mathrm{H}_{2, \text { ad }}}^{n_{1}}, $

(11)

$ \dot{r}_{2}=k_{2} \theta_{\mathrm{H}_{2, \text { ad }}}^{n_{1}} \theta_{\mathrm{O}_{*}}^{n_{1}}-k_{-2} \theta_{\mathrm{OH}}^{2 n_{1}}, $

(12)

$ \dot{r}_{3}=k_{3} \theta_{\mathrm{OH}}^{2 n_{1}}-k_{-3} \theta_{\mathrm{O}_{*}}^{n_{1}} \theta_{\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}}^{n_{1}}, $

(13)

$ \dot{r}_{4}=k_{4} \theta_{\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}}^{n_{1}}-k_{-4} p_{\mathrm{H}_{2} \mathrm{O}} \theta_{\mathrm{O}_{*}}^{n_{1}}, $

(14)

$ \dot{r}_{5}=k_{5} p_{\mathrm{CO}} \theta_{\mathrm{O}_{*}}^{n_{1}}-k_{-5} \theta_{\mathrm{CO}_{\mathrm{ad}}}^{n_{1}}, $

(15)

$ \dot{r}_{6}=k_{6} \theta_{\mathrm{CO}_{\mathrm{ad}}}^{n_{1}}-k_{-6} \theta_{\mathrm{CO}_{2, \text { ad }}}^{n_{1}}, $

(16)

$ \dot{r}_{7}=k_{7} \theta_{\mathrm{CO}_{2, \text { ad }}}^{n_{1}}-k_{-7} p_{\mathrm{CO}_{2}} \theta_{\mathrm{O}_{*}}^{n_{1}} . $

(17)

式中，θ_{H2, ad}、θ_OH、θ_H2O、θ_CO、θ_CO2和θ_O*分别为表面组分H₂、OH、H₂O、CO、CO₂和O的表面覆盖率；n₁为反映晶粒表面存在非均匀结构特性的参数。

获得了每个表面基元反应的反应速率之后，就可以列出微观动力学方程来计算表面组分的生成与消耗^[43]。微观动力学方程是表面组分和气相组分的守恒方程，通过求解微观动力学方程，可以获得不同温度和压力条件下的表面上的每个组分的覆盖率θ_i及其随时间的变化情况，把θ_i代回式(11)—(17)就能得到每个反应的速率，进而可以计算得到单位时间内单个活性位点的转化特性，即转化频率^[44]，同时根据每个反应路径所占比重可以分析选择性以及速率控制步骤^[45]。微观动力学是连接微观理论计算和宏观实验现象的桥梁，是处理表面反应的有效理论工具，对于提升表面反应计算水平、深入分析表面反应过程具有重要的意义^[46]。

虽然微观动力学在催化领域获得了广泛应用，但其无法应用于气固反应动力学的计算^[46]。气固催化反应所研究的主要是材料的表界面化学/物理过程，材料的体相近似按惰性物质处理，虽然其与表面原子有相互作用，但体相物质不参与反应。而在如图 3所示的气固非催化反应体系中，不仅表界面的原子参与反应，而且体相的物质也通过固态扩散的方式与气体分子接触参与反应，并伴随着固体结构的动态演变。微观动力学只描述表面组分覆盖率θ_i随时间的变化情况，描述对象是二维固体表面^[46]。而图 3所示的气固反应不仅仅包含表面反应，还包含固体物质在晶粒内部的体相扩散，同时还伴随着结构变化，实际上是三维体系。所以，原有的研究只能描述二维表面的微观动力学理论，无法处理如图 3所示的三维气固反应对象。气固反应中的晶粒是由大量原子构成的三维结构，晶粒的特征尺度为10²~10⁴ nm。量子化学和分子动力学模拟可以细致地描述原子分子尺度的反应过程与特性，但构成晶粒的大量原子/分子所涌现的集体的气固反应行为却无法从底层的量子化学和分子动力学理论中推导出来，现有理论和方法还不能从原子分子尺度出发直接计算预测宏观的颗粒尺度的现象和动力学行为。因此，现有理论和方法对描述气固反应动力学的跨尺度现象与机理面临根本性困难，需要新的理论方法。

表面组分H₂、OH、H₂O、CO、CO₂的表面覆盖率随时间的变化速率方程为：

$ \frac{\partial \theta_{\mathrm{H}_{2, \mathrm{ad}}}}{\partial t}=\dot{r}_{1}-\dot{r}_{2}, $

(18)

$ \frac{\partial \theta_{\mathrm{OH}}}{\partial t}=2 \dot{r}_{2}-2 \dot{r}_{3}, $

(19)

$ \frac{\partial \theta_{\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}}}{\partial t}=\dot{r}_{3}-\dot{r}_{4}, $

(20)

$ \frac{\partial \theta_{\mathrm{CO}_{\mathrm{ad}}}}{\partial t}=\dot{r}_{5}-\dot{r}_{6}, $

(21)

$ \frac{\partial \theta_{\mathrm{CO}_{2, \mathrm{ad}}}}{\partial t}=\dot{r}_{6}-\dot{r}_{7}. $

(22)

对于产物层的离子扩散模式，由于产物层内部的氧离子浓度C_O*远高于产物层外部的氧离子浓度C_{O*, shell}, 存在浓度梯度，固体反应物的氧离子通过产物层从固体反应物/固体产物的交界面扩散到固体产物/气体的交界面，然后与气体分子发生化学反应，如图 3所示，所以在建立表面氧组分的守恒方程时，不仅需要考虑表面反应，还需要包含体相晶格氧的扩散，即

$ \frac{\partial \theta_{\mathrm{O}_{*}}}{\partial t}=-\dot{r}_{1}-\dot{r}_{2}+\dot{r}_{3}+\dot{r}_{4}-\dot{r}_{5}+\dot{r}_{7}+J_{\mathrm{O}_{*}} / \varLambda_{\mathrm{O}_{*}}. $

(23)

式(23)^[47]中，θ_O*是表面组分O的覆盖率，J_O*是氧离子通过产物层的扩散通量(mol/m²·s^-1)，Λ_O*是表面活性位密度(mol/m²)。式(23)是描述晶粒尺度的反应/扩散方程，实现了气固反应的跨尺度关联以及将反应与扩散耦合在一起。在该式中，固体反应物质向反应发生的表界面位置的扩散速率和表面基元反应速率之间的竞争决定了表面组分的生成消耗及表界面的结构变化。式(23)中扩散通量J_O*可以通过r₁处的氧离子浓度梯度进行计算，即

$ J_{\mathrm{O}_{*}}=-\left.D_{\mathrm{O}_{*}} \frac{\mathrm{d} C_{i}}{\mathrm{~d} r}\right|_{r=r_{1}} . $

(24)

氧离子通过产物层的扩散通过Fick扩散定律描述，

$ \frac{\partial C_{i}}{\partial t}=\frac{D_{\mathrm{O}_{*}}}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial C_{i}}{\partial r}\right) . $

(25)

扩散方程(25)的边界条件为：

$ r=r_{2} \text { 时, } C_{i}=C_{\mathrm{O}_{*}}, $

(26)

$ r=r_{1} \text { 时, } C_{i}=C_{\mathrm{O}_{*}, \text { shell }} \text {. } $

(27)

产物层外侧, 即晶粒表面处的氧离子浓度与表面组分的覆盖率之间存在如下关系：

$ \begin{gathered} C_{\mathrm{O}_{*}, \text { shell }}=\left(\theta_{\mathrm{H}_{2, \text { ad }}}+\theta_{\mathrm{OH}}+\theta_{\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}}+\right. \\ \left.\theta_{\mathrm{CO}_{\mathrm{ad}}}+\theta_{\mathrm{CO}_{2, \text { ad }}}+\theta_{\mathrm{O}_{*}}\right) C_{\mathrm{O}_{*}} . \end{gathered} $

(28)

此外，随着产物层扩散和反应界面处化学反应的进行，固体反应物和气体分子不断被消耗，同时固体产物层的厚度逐渐增加，即r₂在反应中逐渐减小，r₂随时间的变化满足^[48]

$ \frac{\mathrm{d} r_{2}}{\mathrm{~d} t}=-\frac{r_{1}^{2}}{r_{2}^{2}} \frac{\varLambda_{\mathrm{O}_{*}} M_{\mathrm{O}}}{\rho_{\mathrm{oc}} R_{\mathrm{O}}}\left(\dot{r}_{4}+\dot{r}_{7}\right) . $

(29)

式中，M_O、ρ_oc和R_O分别为氧的质量(kg)、晶粒密度(kg/m³)和载氧量(%)。

晶粒表面基元反应与固体物质体相扩散将导致晶粒结构发生改变，如图 3所示。即晶粒内部的反应物质以缩核方式收缩，未反应核外面存在一层固体产物，反应物以离子扩散形式向外传递，到达晶粒表面处与气体分子发生反应。式(29)就是描述表界面反应与体相扩散所导致的晶粒结构变化的方程，是晶粒尺度所特有的方程^[48]。晶粒的转化速率的计算式为

$ \frac{\mathrm{d} X}{\mathrm{~d} t}=\frac{S_{0}}{\left(1-\varepsilon_{0}\right)}[1+(Z-1) X]^{2 / 3} \frac{M_{\mathrm{O}} \varLambda_{\mathrm{O}_{*}}}{\rho_{\mathrm{oc}} R_{\mathrm{O}}}\left(\dot{r}_{4}+\dot{r}_{7}\right). $

(30)

表面组分的覆盖率之和应该为1，即

$ \theta_{\mathrm{H}_{2, \text { ad }}}+\theta_{\mathrm{OH}}+\theta_{\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}}+\theta_{\mathrm{CO}_{\mathrm{ad}}}+\theta_{\mathrm{CO}_{2, \text { ad }}}+\theta_{\mathrm{O}_{*}}+\theta_{\mathrm{V}}=1 . $

(31)

初始条件如下：

$ \begin{gathered} t=0, \\ \theta_{\mathrm{H}_{2, \text { ad }}}=\theta_{\mathrm{CO}_{\mathrm{ad}}}=0, \\ \theta_{\mathrm{OH}}=\theta_{\mathrm{H}_{2} \mathrm{O}_{\mathrm{ad}}}=\theta_{\mathrm{CO}_{2, \text { ad }}}=0, \\ \theta_{\mathrm{O}_{*}}=1, \theta_{\mathrm{V}}=0, X=0 . \end{gathered} $

(32)

式中，θ_V为氧空位覆盖率。

联立式(1)—(32)，可以计算出晶粒表面组分覆盖率θ、固体转化率X、晶粒结构r₁和r₂等参数随时间变化特性，因而获得气固反应动力学特性。反应速率采用式(11)—(17)来计算获得，其中的反应速率常数通过式(7)—(10)计算。过渡态理论中的分子振动频率、转动惯量、反应能垒等基本参数可以通过量子化学DFT计算获得。式(30)是属于晶粒尺度的理论，用于描述构成晶粒的大量原子/分子所涌现的集体的气固反应行为；式(30)建立了原子分子尺度理论与宏观颗粒尺度所呈现的气固反应现象之间的理论联系，实现了跨尺度关联与计算。

5 表面反应的速率控制步骤分析

表面反应机理涉及多个基元反应，其中多数基元反应是快反应，一般只有几个基元反应是慢反应，这几个慢反应则为速率控制步骤，决定整个表面反应的进程。对于一个复杂的气固反应系统，尤其是当表面基元反应数量较多时，有必要对反应机理与反应路径进行分析，找到速率控制步骤，从而为优化反应性能提供关键信息。

表面反应机理涉及多个基元反应，其中多数基元反应是快反应，一般只有几个基元反应是慢反应，这几个慢反应则为速率控制步骤，决定整个表面反应的进程。对于一个复杂的气固反应系统，尤其是当表面基元反应数量较多时，有必要对反应机理与反应路径进行分析，找到速率控制步骤，从而为优化反应性能提供关键信息。

表面反应的速率控制步骤分析方法有敏感性分析方法^[49]、基元反应的可逆性分析方法^[50]、速率控制度分析方法^[51]等。其中，速率控制度分析方法(degree of rate control, DRC)近年来获得了广泛应用^[51]，是一种比较有效的速率控制步骤分析方法，本文采用DRC分析气固反应中的速率控制步骤。在DRC分子中，第i个基元反应的速率控制度(x_{DRC, i})定义为^[51]

$ x_{\mathrm{DRC}, i}=\left.\frac{k_{i}}{r_{\text {overall }}}\left(\frac{\partial r_{\text {overall }}}{\partial k_{i}}\right)\right|_{k_{j \neq i}, K_{i}} \cdot $

(33)

其中，K_i是第i个基元反应的平衡常数，整体反应速率r_overall可以采用固体的转化速率^[47]，即

$ r_{\text {overall }}=\frac{\partial X}{\partial t}. $

(34)

参数x_{DRC, i}表征了速率常数k_i的相对变化对反应速率r_overall的影响程度；x_{DRC, i}值越大，说明第i个基元反应对r_overall的影响程度越大。

通过对式(33)的进一步推导可以得到第i个基元反应的速率控制度与其能垒间的关系，见式(35)^[51]。

$ x_{\mathrm{DRC}, i}=\left.\left(\frac{\partial \ln r_{\text {overall }}}{\partial \ln k_{i}}\right)\right|_{k_{j \neq i}, K_{i}}=\left.\left(\frac{\partial \ln r_{\text {overall }}}{\partial\left(\frac{-E_{\mathrm{A}, i}}{k_{\mathrm{B}} T}\right)}\right)\right|_{k_{j \neq i}, K_{i}} $

(35)

式(35)中的能垒E_{A, i}可以通过DFT计算得到，可以看出DRC分析方法直接建立了其与量子化学计算间的关系，可以从更微观尺度分析表面反应的控制步骤以及能垒对基元反应的定量影响。

对Fe₂O₃与H₂的反应的计算结果进行DRC分析，得到的结果如图 6所示，可以看出，R2反应的x_{DRC, 1}为1，而其他反应的x_DRC为0，说明在Fe₂O₃与H₂的表面反应过程中，R2反应即H₂在Fe₂O₃表面的分解反应是速率控制步骤。插入吸附态H₂在Fe₂O₃表面的分解属于Langmuir-Hinshelwood类型的表面反应。

6 体相扩散的影响

在如图 3所示的晶粒模型中，以Fe₂O₃与H₂的反应为例，对于产物层的离子扩散模式，氧离子通过产物层从固体反应物/固体产物的交界面扩散到固体产物/气体的交界面，然后与气体分子发生化学反应。氧离子向反应发生的位置传递的速率和氧离子与H₂的表面反应速率决定了整个反应过程、生成物及其表界面的结构，从而影响固体转化行为。随着产物层扩散和反应界面处化学反应的进行，固体反应物和气体分子不断被消耗，同时固体产物层的厚度逐渐增加。在式(24)和(25)中，产物层扩散用扩散系数D_O*来表征，固体产物层结构对反应传递的影响体现在与转化率X有关的函数项上, 见式(30)。

为了定量衡量反应与扩散之间的竞争关系，定义无量纲Thiele模数为^[48]

$ \phi^{2}=\frac{k_{2} r_{0} \varLambda_{\mathrm{O}_{*}}}{c_{0} D_{\mathrm{O}_{*}}}=\frac{k_{2}}{\frac{c_{0} D_{\mathrm{O}_{*}}}{r_{0} \varLambda_{\mathrm{O}_{*}}}}=\frac{k_{2}}{D^{*}} . $

(36)

不同Thiele模数下的固体转化率随时间变化的曲线见图 7，可以看出，当ϕ≤0.1时，固态扩散对气固反应的影响可以忽略；而当ϕ>0.1时，则需要在气固反应动力学的计算中考虑固态扩散，否则会引起较大误差。

7 简化模型

式(1)—(32)所建立的速率方程理论模型包含多个微分方程，这些非线性微分方程组需要数值求解，且算法中包含多步迭代求解过程。将上述需要迭代数值求解的速率方程应用于模拟反应器尺度对象的设计优化或数值模拟时，其数值求解过程将非常复杂，且计算量巨大。在实际应用情况下，研究者更倾向于采用具有简单解析解形式的模型。因此，有必要针对速率方程理论模型的数值求解算法进行改进，建立解析形式的快速计算简化模型，减少计算时间和成本。

例如，对于Fe₂O₃与H₂的反应，经过表面反应速率控制步骤分析，发现H₂在Fe₂O₃表面的分解反应是速率控制步骤，同时考虑到表面反应是非常快的过程，针对表面组分O，认为其处于准稳态状态，联立式(18)—(30)，并令(11)—(15)中的n₁=1时，得到如下简化模型：

$ \begin{gathered} \frac{\mathrm{d} X}{\mathrm{~d} t}=\frac{k_{1} \varLambda_{\mathrm{O}_{*}} {~S}_{0}}{\left(1-\varepsilon_{0}\right)} \frac{M_{\mathrm{O}}}{\rho_{\mathrm{oc}} R_{0}} \times \\ \frac{\kappa_{\mathrm{S}}(X)}{g_{\mathrm{D}}(X)+\frac{k_{1}}{D_{\mathrm{O}_{*}}} \frac{M_{\mathrm{O}}\varLambda_{\mathrm{O}_{*}}}{\rho_{\mathrm{oc}} R_{\mathrm{O}}} \frac{\left(1-\varepsilon_{0}\right)}{S_{0}} \cdot p_{\mathrm{D}}(X)}. \end{gathered} $

(37)

式(37)中的函数κ_S(X)、g_D(X)、p_D(X)分别表示为：

$ \kappa_{\mathrm{S}}(X)=(1-X)^{2 / 3}, $

(38)

$ g_{\mathrm{D}}(X)=\frac{(1-X)^{2 / 3}}{[1+(Z-1) X]^{2 / 3}}, $

(39)

$ \begin{gathered} p_{\mathrm{D}}(X)=\frac{3(1-X)^{1 / 3}}{[1+(Z-1) X]^{1 / 3}}\cdot \\ {\left[[1+(Z-1) X]^{1 / 3}-(1-X)^{1 / 3}\right]}. \end{gathered} $

(40)

采用简化模型式(37)计算Fe₂O₃与H₂的反应，并与详细模型的计算结果进行对比，如图 8所示。可以看出，简化模型计算结果与详细模型计算结果完全一致，说明在特定条件下(此例是n₁=1且R2是速率控制步骤)，简化模型完全可以准确计算气固反应动力学。

值得指出的是，简化模型即式(37)—(40)的形式与传统的缩核模型形式完全一致^[13-14]，说明式(1)—(32)所表示的速率方程是通用的理论模型，传统的缩核模型只是特定条件下速率方程理论的一个简化形式。简化模型只考虑了单步反应，无法描述比较复杂的表面反应机理与速控步骤；而详细的速率方程理论模型则能够细致地描述表面基元与固态扩散的耦合以及晶粒结构的演变，所以适用性更广。

8 固体产物岛状生长模式

在气固反应研究中，人们普遍采用连续产物层假设，即认为固体产物以连续形式进行生长，随着反应的进行，产物层厚度逐渐增大。然而，基于连续产物层的假设在实际应用中存在一些局限性。对于一些气固反应，如CaO碳酸化反应、CaO的硫化反应、Ca(OH)₂的脱硫反应等，研究发现转化率曲线存在由初始的快速反应阶段向之后的慢速反应阶段急剧转折的现象^[52-54]，如图 9所示：在反应初始阶段，反应速率非常快；当转化率到达一个临界值后，反应速率迅速变慢，转化率曲线呈现两阶段特征。此外，研究者发现温度对临界转化率有显著的影响，温度越高则相应的临界转化率越高^[52-54]。基于连续产物层假设的传统模型计算结果为渐变型特征，无法描述反应动力学从快速阶段向产物层扩散阶段发生急剧转折的实验现象，也无法解释动力学转折处临界转化率随温度升高而增加的物理化学机制，只能通过拟合实验数据的方式来描述反应动力学。

本文采用具有光滑平整表面的单晶进行气固反应实验，通过原子力显微镜观测反应过程中生成的固体产物形貌，研究对象包括MgO硫酸化反应^[55]、CaO碳酸化反应^[52]、CaCO₃硫化反应^[56]及Fe氧化反应^[57]。通过对不同对象的研究均发现，在反应初始阶段，固体产物并不是以连续形式覆盖在固体反应物表面，而是呈现分散的三维岛状结构，如图 10和11所示。随着反应的进行，固体反应物表面的产物岛数量逐渐增加，产物岛之间的空隙逐渐被填补。同时，由于大产物岛比表面积较低，表面能较小，因此反应生成的固体产物经由表面扩散后容易被大岛捕获，使得产物岛尺寸逐渐增大。最终产物岛将完全覆盖整个固体反应物表面，阻隔气体反应物和固体反应物的直接接触，反应进入产物层扩散阶段。加拿大渥太华大学的Arias等^[58]的研究中也指出反应进程发生转折的根本原因是分散的固体产物逐渐包裹固体反应物表面。Fang等^[55]的研究表明反应温度对产物岛的尺寸有明显的影响。在低温下，产物岛呈现尺寸小、数量多的特点；而高温下产物岛则呈现尺寸大、数量少的特点。气固反应动力学转折处的临界转化率与产物岛尺寸直接相关，产物岛尺寸随温度升高而增大，因此临界转化率随温度的升高而增大。

气固反应中的固体产物是以三维岛状形式在反应物表面生长的，固体产物在反应过程中是不稳定的。在反应初始阶段，成核与生长相互竞争，固体产物数量受到化学反应速率和扩散速率影响。当化学反应速率高于扩散速率时，固体产物数量增加、尺寸较小。而当扩散速率远大于化学反应速率时，固体产物数量减少、尺寸增大。化学反应速率和扩散速率的相互竞争影响着固体转化行为。基于离散三维产物岛观测结果，本研究团队^[58-59]建立了速率方程理论，从微观尺度详细描述产物岛生长过程。该理论综合考虑了表面化学反应、单体表面扩散及捕获与逃逸、产物岛成核与聚并等步骤来计算得到产物岛尺寸的分布。定义δ表示固体反应物表面裸露份额，(1-δ)则表示固体反应物表面被产物岛覆盖的份额。表面裸露份额δ变化速率的表达式为^[60]

$ \frac{\mathrm{d} \delta}{\mathrm{d} t}=-\frac{k_{\mathrm{s}}}{h_{\mathrm{c}}} Z \delta\left(C-C_{\mathrm{eq}}\right). $

(41)

$ \begin{gathered} \frac{\mathrm{d} X}{\mathrm{~d} t}=\frac{k_{1} \varLambda_{\mathrm{O}_{*}} S_{0}}{\left(1-\varepsilon_{0}\right)} \frac{M_{\mathrm{O}}}{\rho_{\mathrm{oc}} R_{\mathrm{O}}} \kappa_{\mathrm{S}}(X)\cdot \\ {\left[\delta+\frac{1-\delta}{g_{\mathrm{D}}(X)+\beta \cdot p_{\mathrm{D}}(X)}\right]}. \end{gathered} $

(42)

计算得到的转化率曲线呈现典型的动力学转折特征，如图 12所示，初始阶段固体反应物表面部分处于裸露状态，可直接与气体分子发生化学反应，此阶段为快速反应动力学阶段；随着反应的进行，产物岛逐渐完全覆盖固体反应物表面，即δ→0，之后反应进入产物层扩散阶段，其反应速率远慢于初始阶段反应速率。可以看出，在理论上采用离散岛状生长模式来描述固体产物生长，代替传统模型中连续产物层假设，可以合理刻画反应过程中固体微观结构变化，揭示反应动力学从快速阶段向产物层扩散阶段转变特征的本质机理。

9 速率方程理论与热质传递的耦合 9.1 颗粒内扩散

上述的速率方程理论描述了晶粒尺度的表界面反应与体相扩散行为，而实际的固体颗粒则是由很多晶粒组成，且晶粒之间存在孔隙。气体扩散进入颗粒内部到达晶粒表面时开始发生化学反应，并生成固体产物，在颗粒内部存在着气体浓度梯度。气体在多孔固体颗粒内部的扩散/反应过程需要求解单颗粒模型。单颗粒模型分为详细数值求解算法^[61]和快速求解算法^[62]，如图 13所示。在详细数值求解算法中，颗粒沿径向被划分为多个微元薄层，用于详细求解颗粒内气体浓度分布曲线，进而获取其他参数的分布值，如图 13a所示；在快速求解算法中，不考虑颗粒内部具体的气体浓度分布值，仅考虑颗粒表面处的气体浓度值，如图 13b所示，采用Thiele模数表征气体通过颗粒孔隙的扩散阻力对颗粒内气体浓度分布的影响。在颗粒表面处，反应气体的浓度为C_s，如图 13所示，则气体浓度驱动力为(C_s－C_eq)。根据颗粒整体的效率因子η的定义，可知颗粒整体平均的气体浓度驱动力可用η(C_s－C_eq)表示。因此，从颗粒平均角度出发，用于计算颗粒整体转化速率的气体浓度项可简化为^[63]

$ C-C_{\mathrm{eq}}=\eta\left(C_{\mathrm{s}}-C_{\mathrm{eq}}\right). $

(43)

快速求解算法中颗粒的效率因子具体计算式为

$ \eta=\frac{3}{\varphi^{2}}(\varphi \operatorname{coth} \varphi-1). $

(44)

其中：φ为Thiele模数，η为效率因子，C_s为颗粒表面处气体浓度，C_eq为气体平衡浓度，C为颗粒内部气体浓度。

在反应过程中，颗粒孔隙结构随固体产物的生长而逐渐变化。例如，对于固体产物体积大于反应物体积的情况，产物生长将导致颗粒内孔隙尺寸逐渐减小，使气体内扩散阻力增大。实际上颗粒内的孔隙尺寸变化是沿颗粒径向存在分布的，且颗粒外表面处的孔隙尺寸减小得最快，故最先出现孔隙堵塞现象。为了描述颗粒孔隙结构变化对气体内扩散行为的实时影响，在快速求解算法中采用颗粒表面处气体在孔内的有效扩散系数来计算Thiele模数^[63]。当颗粒表面发生孔隙堵塞现象时，颗粒表面处气体在孔内的有效扩散系数立即减小为0，计算得到的Thiele模数为无穷大，从而可描述颗粒反应进程停止的情况，如图 14所示。

9.2 颗粒外扩散

对于单颗粒固体反应物，颗粒周围气相主流的反应气体浓度为C_A，气体从气相主体向颗粒外表面扩散的过程中由于外扩散阻力的影响而存在气体浓度差，即颗粒表面气体浓度C_s小于C_A，如图 14所示。从气体外扩散角度，颗粒整体消耗气体的反应速率(mol/s)可表示为

$ F_{\mathrm{A}}=4 {\rm{ \mathsf{ π} }} R_{0}^{2} k_{\mathrm{g}}\left(C_{\mathrm{A}}-C_{\mathrm{s}}\right). $

(45)

从颗粒固体反应物本身来看，其消耗气体的反应速率可表示为

$ F_{\mathrm{A}}=\nu_{\mathrm{A}} \chi_{\mathrm{R}} \frac{\rho_{\mathrm{p}}}{M_{\mathrm{R}}} \frac{4}{3} {\rm{ \mathsf{ π} }} R_{0}^{3} \frac{\mathrm{d} X}{\mathrm{~d} t}. $

(46)

式(45)和(46)中，k_g为流体与颗粒间的传质系数，R₀为颗粒半径(m)，χ_R为颗粒中活性物质含量(%)，M_R为反应物分子量(kg/mol)，ρ_p为颗粒密度(kg/m³)，ν_A为化学计量系数。

联立方程(45)和(46)以及固体转化速率表达式就能推导得到包括外扩散影响的气固反应动力学速率方程，详细推导过程见文^[63]，此处不再赘述。

9.3 颗粒温度

在固体颗粒与气体发生气固反应过程中，伴随着由化学反应导致的热量释放或吸收，以及颗粒与周围介质的辐射传热和对流传热，使得颗粒温度发生改变。通过对固体颗粒列出能量守恒方程，

$ \begin{gathered} m_{\mathrm{p}} c_{\mathrm{p}} \frac{\mathrm{d} T_{\mathrm{p}}}{\mathrm{d} t}={\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2} h_{\mathrm{k}}\left(T_{\mathrm{g}}-T_{\mathrm{p}}\right)+\\ {\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2} \varepsilon_{\mathrm{p}} \sigma\left(T_{\mathrm{r}}^{4}-T_{\mathrm{p}}^{4}\right)-m_{\mathrm{p}, 0} \cdot \frac{\mathrm{d} X}{\mathrm{~d} t} \cdot \Delta H_{0}. \end{gathered} $

(47)

其中，m_p是颗粒质量(kg)，m_{p, 0}是颗粒初始质量(kg)，c_p是颗粒比热容(J/(kg·K^-1))，T_p是颗粒温度，d_p是颗粒直径，T_g是气体温度，T_r是辐射温度，ΔH是反应焓(J)。式(47)中，对流换热系数公式为

$ h_{\mathrm{k}}=\frac{\lambda_{\mathrm{g}}}{d_{\mathrm{p}, \mathrm{k}}}\left(2.0+0.6 {Re}_{\mathrm{g}, d_{\mathrm{p}, \mathrm{k}}}^{1 / 2} {Pr}_{\mathrm{g}}^{1 / 3}\right). $

(48)

方程(47)的计算需要数值求解算法，但数值算法的计算量较大，会耗费系统较大计算资源。可以采用近似解析方法求解方程(47)，即假设颗粒温度在一个很小的时间步长内变化较小，所以可以对方程(47)进行线性化处理，得到近似解析解，如下^[64]：

$ T_{\mathrm{p}}(t+\Delta t)=\alpha_{\mathrm{p}}+\left[T_{\mathrm{p}}(t)-\alpha_{\mathrm{p}}\right] \mathrm{e}^{-\gamma_{\mathrm{p}} \cdot \Delta t}, $

(49)

$ \alpha_{\mathrm{p}}=\frac{h_{\mathrm{k}} {\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2} T_{\mathrm{g}}+m_{\mathrm{p}, 0} \cdot \frac{\mathrm{d} X}{\mathrm{~d} t} \cdot \Delta H+{\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2} \varepsilon_{\mathrm{p}} \sigma T_{\mathrm{r}}^{4}}{h_{\mathrm{k}} {\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2}+h {\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2} \varepsilon_{\mathrm{p}} \sigma T_{\mathrm{p}}^{3}}, $

(50)

$ \gamma_{\mathrm{p}}=\frac{{\rm{ \mathsf{ π} }} d_{\mathrm{p}}^{2}\left(h_{\mathrm{k}}+\varepsilon_{\mathrm{p}} \sigma T_{\mathrm{p}}^{3}\right)}{m_{\mathrm{p}} c_{\mathrm{p}}}. $

(51)

10 速率方程理论的应用 10.1 氧化物的还原

氧化物的还原是一类重要的气固反应，是用还原剂将高价态金属还原至低价态或还原至金属。氧化物的还原广泛应用于钢铁冶炼、有色金属提取、化学链燃烧、制氢等过程。人们通常采用表观模型拟合实验数据的方法来获得反应动力学参数，采用总包反应来处理化学反应，没有考虑表面反应的微观机理。一些研究者采用DFT计算还原过程的反应路径、分析反应能垒，但缺乏将DFT计算结果与实验测量得到的宏观动力学关联，导致微观计算与宏观现象脱节。作者团队利用所建立的速率方程理论，从量子化学DFT计算获得Fe₂O₃与H₂表面反应路径(如图 4和图 5)，利用TST计算速率常数，然后采用速率方程理论计算还原反应动力学，计算中的模型输入参数来自于DFT计算数据，并将计算结果与实验数据进行对比，如图 15所示。可以看出，所建立的基于第一性原理的速率方程理论能够很好地预测Fe₂O₃与H₂反应动力学规律。

采用速率方程理论计算钙钛矿与CO的还原反应动力学，计算结果与实验数据的对比如图 16所示，可以看出模型计算结果与实验数据吻合较好。说明所建立的理论模型能够实现还原反应中不同物质与不同还原性气体间的反应动力学预测。

10.2 氧化反应动力学

氧化反应在工业生产、材料防护和日常生活中有着广泛的应用。高温下的氧化反应速率非常快，一般采用热重分析等仪器测量高温下的氧化反应动力学，但热重实验中的气体扩散对反应的影响难以消除，实验结果会低估氧化反应速率；基于气体组分信号测量的微型流化床在测量这类仅有气体消耗的快速反应时，存在气体信号延时和失真等困难，目前尚存在较大挑战；本研究团队提出并建立了新型的测量气固反应动力学的微型流化床热重分析方法，解决了流态化条件下固体颗粒质量信号的稳定、准确测量难题。基于微型流化床热重分析法，测试了钙钛矿的氧化反应动力学特性。利用量子化学DFT计算了钙钛矿与O₂表面反应路径并将DFT计算得到的能量值和振动频率等数据代入TST计算速率常数，然后采用速率方程理论计算钙钛矿的氧化反应动力学，计算结果与实验数据的对比如图 17所示，可以看出速率方程理论能够准确预测温度和O₂含量等对氧化反应动力学的影响特性。

10.3 CaO与CO₂的碳酸化反应

CaO与CO₂的碳酸化反应及CaCO₃的分解反应可以应用于烟气CO₂捕集、吸附增强式气化或制氢、热化学储热等^[63]。CaO与CO₂的碳酸化反应动力学是实现上述过程的关键。本研究团队采用微型流化床热重测试了不同条件下的CaO碳酸化反应动力学，利用速率方程理论计算了反应动力学特性^[63]，结果如图 18所示。该图表明，在同一温度下，改变CO₂分压并不影响快速阶段结束时的临界转化率。当气体体积分数从25%增加到75%时，快速阶段的反应速率随气体体积分数的增加而明显增加。同时研究了颗粒粒径的影响，测试了粒径在0.2~3 mm之间的CaO的碳酸化反应动力学，结果如图 18b所示。可以看出，当粒径由0.2~0.3 mm增加到2.5~3 mm时，快速阶段的碳酸化反应速率明显减慢。粒径较大时，气体外扩散和内扩散阻力均较大。从图 18可以看出，所建立的速率方程理论可以准确地描述CaO碳酸化反应动力学特性以及反应动力学由快速阶段向慢速阶段转折的现象^[63]。

10.4 Ca(OH)₂与SO₂和SO₃的反应

在燃煤电站中，随着烟气超低排放标准的实施，锅炉尾部需要加装SCR脱硝装置脱除烟气中的NO_x，SCR催化剂会导致烟气中SO₃含量增加，需要进一步脱除SO₃。一般采用向烟道喷射Ca(OH)₂粉末的方式脱除SO₃，但烟气中的SO₂也会与Ca(OH)₂发生反应，因此SO₂和SO₃与Ca(OH)₂的竞争反应动力学对于脱除工业烟气中SO₃的性能尤为重要。宋蔷课题组^[67]采用气固反应的速率方程理论研究了不同温度下Ca(OH)₂与SO₂和SO₃的竞争反应动力学特性，讨论了影响Ca(OH)₂转化率的因素，并与实验结果进行对比验证，结果如图 19所示。结果表明，所建立的速率方程理论可以准确地描述Ca(OH)₂与SO₂反应动力学特性以及反应动力学由快速阶段向慢速阶段转折的现象，也可以准确地预测出临界转化率随反应温度升高而增加的实验现象，如图 19a。理论也揭示出SO₃的产物层扩散速率快于SO₂的产物层扩散速率，如图 19b所示。

进一步分析Ca(OH)₂表界面反应与SO₃气体外扩散间的竞争关系，研究发现Ca(OH)₂粒径对SO₃的脱除效率有着重要影响。在喷入的Ca(OH)₂质量一定条件下，当Ca(OH)₂颗粒粒径大于20 μm时，SO₃脱除效率基本不受反应动力学控制，此时SO₃的外扩散是控制步骤^[68]。当颗粒粒径大于10 μm之后，进一步增加表面反应活性几乎不影响脱硫剂的净反应速率；而当颗粒粒径小于6 μm时，脱硫剂反应活性的改善可以较大幅度增加净反应速率。所以，在实际喷射脱硫剂时，应把其粒径控制在6 μm之内^[68]，这时的SO₃气体外扩散速率与气固表面反应速率几乎在同一范围内，可以有效发挥脱硫剂的高活性优势。

10.5 化学链燃烧流化床反应器的设计

化学链燃烧是通过载氧体在空气反应器和燃料反应器之间进行氧化还原循环，将空气中的氧传递给燃料，避免CO₂气被氮气稀释，实现了CO₂在燃料转化源头的自分离，具有较低的CO₂捕集能耗和成本。化学链燃烧的空气反应器通常设计为快速床、鼓泡床、湍动床、移动床等形式。利用跨尺度的速率方程理论计算分析表明，常规实验方法测试得到的载氧体氧化反应动力学偏慢，导致化学链空气反应器的设计不合理。本文的理论预测与实验结果表明：载氧体具有非常快的氧化速率、可以在5 s内实现完全氧化反应。基于载氧体的快速氧化反应动力学，设计了输送床形式的化学链燃烧空气反应器，如图 20所示，气速可以达到8~12 m/s，对应的固体循环量为29.5~44.5 kg/(m²·s)，解决了化学链双流化床间如何实现高循环通量的难题。

化学链燃烧系统的关键问题之一是空气反应器中载氧体的再循环比对氧化反应速率的影响。采用所建立的跨尺度的速率方程理论模型对化学链燃烧系统进行计算与分析，揭示了在空气反应器内经历了不同循环次数的载氧体颗粒的反应速率会随着循环比的增加而降低，如图 21所示。进而基于质量平衡、热量平衡、流态化模型、反应动力学和反应器模型，建立了化学链燃烧系统设计平台并编写了热力计算程序。运用该平台和程序设计了热输入为3~5 MWth的化学链燃烧热态示范装置，为化学链燃烧热态示范装置的建造和运行奠定了基础。

10.6 煤焦燃烧的CFD模拟

煤粉低氮燃烧中还原区的存在使得焦炭的整体燃烧过程被推迟，煤焦反应速率下降，焦炭面临燃尽的问题。因此，需要在焦炭燃烧动力学基础上加入反应速率下降的机理描述，以准确预测煤粉低氮燃烧导致的焦炭燃尽率下降与飞灰含碳量升高现象。除此之外，煤粉低氮燃烧中焦炭的粒径衍化特性格外复杂。不同特征区域(主燃区、还原区、燃尽区)内燃烧气氛差异较大，焦炭燃烧速率的控制步骤各不相同，粒径、密度随碳质消耗的衍化模式多次迁移。颗粒尺寸会显著影响煤粉在锅炉中的燃烧速率与停留时间，因此，需要在已有理论模型的基础上考虑粒径、密度随碳质消耗的衍化特性，建立可靠的机理描述。因此，本研究团队在焦炭燃烧本征动力学模型的基础上，通过引入灰层阻碍氧化剂扩散的数学描述考虑了燃烧过程中反应速率的逐渐下降^[71]；通过分析焦炭燃烧的一维详细模型，建立了粒径、密度衍化问题可靠的零维简化描述^[71]。通过编写UDF将焦炭燃烧改进模型嵌入Fluent，实现了新的焦炭燃烧模型与热解、流动、传热过程的耦合计算^[71]，如图 22所示。以沉降炉、一维炉的实验数据为基准，对焦炭燃烧改进模型进行了对比验证。结果表明，与传统模型相比，焦炭燃烧改进模型可以更为准确的计算燃烧全过程的燃尽率。所建立的理论模型可以合理描述不同燃烧条件下焦炭粒径分布随碳质消耗的衍化特性。可以准确计算低氮燃烧中主要气体组分浓度的沿程衍化特点^[71]。利用涵盖褐煤、烟煤、无烟煤等不同煤阶大量煤种的实验数据，积累了焦炭燃烧改进模型的参数数据库，可以为新煤种的参数选取提供经验指导。

11 总结

气固反应中的跨尺度关联和反应/扩散耦合是实现流动、传递和反应耦合的核心。针对这一问题，本研究团队在介观晶粒尺度提出了新的理论模型，揭示了气固表界面反应中固体产物形成与生长的现象与机理，设计高效可靠的计算和分析方法，并应用于能源清洁转化等技术领域。取得的主要结论如下：

1) 针对量子化学计算和分子动力学模拟局限于微观尺度而无法应用于跨尺度计算的困难，从介观晶粒尺度提出了基于第一性原理的表界面反应与体相扩散相互耦合的速率方程理论，从“原子→晶粒→颗粒”层面建立了微观原子与宏观颗粒之间的理论联系，解决了气固反应动力学的跨尺度理论关联与计算的难题。

2) 针对传统理论无法解释气固反应动力学呈现两阶段特征且临界转化率随温度升高而增加现象，率先提出基于固体单晶的直接观测气固表界面反应的分析方法，首次揭示固体产物的三维岛状生长机制及演化过程，进而提出描述固体产物岛状生长的新理论，成功地解释了CaO碳酸化、载氧体氧化还原、CaO硫化等一系列气固反应动力学现象。

3) 将跨尺度的速率方程理论与设备内的流动、传热、传质相互耦合，设计高效可靠的计算和分析方法，并应用于氧化/还原反应、低能耗CO₂捕集、热化学储热、污染物脱除以及大型电站锅炉燃烧等技术，有力推动了该领域的发展。

12 未来研究展望

1) 气固反应动力学理论的完善。

跨尺度理论计算的准确性较大程度上依赖于DFT计算的精度。DFT在气固反应领域的计算准确性还有待完善。目前的DFT计算还不能完整考虑固体产物的固体结构与形貌、还无法准确计算固体产物的形成与生长过程以及产物层扩散过程。固体产物的形成与生长及产物层扩散对气固反应进程和动力学有着较大影响，如何从原子/分子尺度准确合理描述固体产物的形成与生长及扩散是下一步需要解决的关键问题。

2) 气固反应动力学理论的实际应用。

气固反应动力学速率方程理论目前只是应用于一些反应机理比较简单的气固反应对象上，该理论方法完全可以应用于更加复杂的气固反应体系，如CO₂的还原、烷烃脱氢制烯烃、天然气制合成气、甲烷耦合反应等；也可以拓展应用在大气颗粒物及其环境效应和地球化学等领域的非均相反应的研究上。

3) 气固反应动力学与机器学习的结合。

目前碳中和场景下的能源转化利用、CO₂捕集转化利用、能量存储、新能源等研究的核心是新材料的寻找、筛选、合成、测试表征及大规模制备等。常规的材料筛选与制备是通过实验的方法来寻找合适的元素和优化的配比，工作量巨大且效率非常低。在后续研究中，可以将跨尺度的气固反应动力学速率方程理论与机器学习相结合，即使用从跨尺度理论计算获得的模型属性作为机器学习模型的输入，并利用跨尺度计算中使用的输入参数作为机器模型学习的输出，以“反向”的方式训练机器学习模型。这种将跨尺度计算与机器学习相结合的模式可以实现计算模型的快速优化，来预先筛选新材料，为开发各种新材料提供理论指导与筛选优化。