无论是在传统能源的高效清洁利用，还是在新能源的可再生转化方面，多相反应过程均扮演着重要的角色，其中超过80%的多相反应过程采用多孔催化颗粒以提高反应速率或者合成特定产品。多孔催化颗粒内外表面难以避免发生积碳^[1]，积碳会覆盖活性位点，引起孔道变窄甚至堵塞，造成催化剂活性的迅速下降甚至失活^[2]。为了实现工业上连续生产，积碳失活的催化剂需要及时进行氧化或气化等再生过程，以消除积碳并恢复催化活性^[3]。因此，催化剂消碳再生技术一直是研究者们关注的焦点问题。

随着计算机水平的不断发展，数值模拟已广泛应用于多孔颗粒失活及再生过程。Yang等^[4-5]采用MP-PIC方法，结合流体力学和积碳燃烧动力学模型，模拟了失活催化剂分布和水平挡板对工业紧凑型再生器性能的影响。研究表明，添加水平挡板和改善积碳催化剂分布，可以提高工业FCC再生器的性能。积碳消除过程不仅依赖于反应器宏观流动和传递特性，而且取决于颗粒内部孔隙结构^[6]。Kern等^[7]利用一维伪均质模型耦合再生反应动力学模型，分析了单催化剂颗粒内氧气和积碳负载量径向分布的瞬时变化，结果表明，在内扩散的作用下，氧气浓度和积碳负荷在催化颗粒内部径向呈现梯级分布。在固定床反应器内，孔扩散阻力对反应前沿宽度的影响比轴向扩散阻力更加显著。Zhang等^[8]考虑颗粒内外热质传递，建立了丙烷脱氢Cr₂O₃/Al₂O₃催化剂再生的非均相模型。结果表明，随着反应温度的升高和积碳含量的增加，孔隙扩散对反应的影响更加显著。Reshetnikov等^[9]基于反应器轴向和颗粒径向二维多尺度模型，研究了低氧浓度条件下绝热固定床反应器内CrF₃/MgF₂催化剂再生过程, 结果表明，在进口气体低于初始催化剂温度条件下，氧化再生反应只发生在催化剂颗粒外表面附近，易导致过热和催化剂损坏。

上述宏观模型往往忽略微观孔隙结构对反应和传质过程的影响，而积碳消除行为导致多孔催化颗粒内外表面孔隙结构的动态演变，传递-反应特性又直接依赖于真实复杂的孔隙结构。随着可控技术的发展，基于微观离散模型探索局部非均匀结构对反应和传递过程的影响，逐步成为研究热点^[10]。因此，本文基于格子玻尔兹曼方法，开展多孔催化颗粒消碳行为的孔隙尺度模拟，分析消碳行为引起的孔隙结构演变特性和微纳尺度流动特性，评估气体稀薄效应下的消碳行为。

1 数学模型 1.1 物理模型

本文开展简单孔隙结构和复杂随机孔隙结构的二维孔隙尺度模拟，如图 1(a)所示，简单孔隙结构的孔隙长L，高H，黑色代表积碳(固相)，顶部和底部边界的积碳厚度均为h。采用四参数随机构造方法构建的复杂孔隙结构如图 1(b)所示，该孔隙结构的结构特性与真实催化剂结构相似^[11]。

1.2 格子玻尔兹曼方程

采用格子玻尔兹曼(LB)方法描述气体稀薄效应下微纳孔隙传递过程，大孔存在条件下，催化颗粒孔隙内存在显著的黏性流动^[12]。为了保持高努森数条件下的数值稳定性，本文采用正则化多松弛(MRT)模型来描述流动过程，其表达式为

$ \begin{gathered} f_i\left(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t+\Delta t\right)=f_i^{\mathrm{eq}}(\boldsymbol{x}, t)+\bar{f}_i(\boldsymbol{x}, t)- \\ \boldsymbol{M}^{-1} \boldsymbol{S} \boldsymbol{M} \bar{f}_i(\boldsymbol{x}, t) \end{gathered} $

(1)

式中：$ f_i\left(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t+\Delta t\right)、f_i^{\mathrm{eq}}(\boldsymbol{x}, t)$和$\overline{f_i}(\boldsymbol{x}, t) $分别代表密度分布函数、密度平衡分布函数和正则化密度非平衡态分布函数；x为位置矢量，t为时间，M和S分别为转换矩阵和松弛矩阵^[13]。采用Bosanquet型有效黏度考虑努森层对黏度的影响^[14]。

密度分布函数$f_i^{\mathrm{eq}}(\boldsymbol{x}, t) $可表示为

$ f_i^{\mathrm{eq}}(\boldsymbol{x}, t)=w_i \boldsymbol{\rho}\left(1+\frac{\boldsymbol{u} \cdot \boldsymbol{c}_i}{c_{\mathrm{s}}^2}+\frac{\left(\boldsymbol{u} \cdot \boldsymbol{c}_i\right)^2}{2 c_{\mathrm{s}}^4}-\frac{\boldsymbol{u} \cdot \boldsymbol{u}}{2 c_{\mathrm{s}}^2}\right) $

(2)

式中：w_i、c_i和c_s分别代表格子权重、格子速度和格子声速，ρ和u为密度和速度。

正则化非平衡态分布函数f_i可表示为

$ \overline{f_i}=w_i\left[\frac{1}{c_i^2} H_{\mathrm{e}}^{(2)}\left(\boldsymbol{c}_i / c_{\mathrm{s}}\right) \sum\limits_{j=0}^8 f_j^{\mathrm{neq}} \boldsymbol{c}_j \boldsymbol{c}_j\right] $

(3)

式中，H_e⁽²⁾为二阶Hermite多项式. 非平衡态分布函数f_i^neq可通过$ f_i^{\mathrm{neq}}=f_i-f_i^{\mathrm{eq}}$求得。

对于孔隙尺度扩散和反应过程，采用被动标量模型来描述反应输运过程^[11]，对应的对流扩散方程如下：

$ \frac{\partial C}{\partial t}+(\boldsymbol{u} \cdot \nabla) C=\nabla \cdot(D \nabla C) $

(4)

式中，C和D分别为组分浓度和分子扩散系数。

基于LBM方法，其演化方程为

$ \begin{gathered} g_{C i}\left(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t+\Delta t\right)-g_{\mathrm{C} i}(\boldsymbol{x}, t)= \\ -\tau_{\mathrm{gs}}^{-1}\left[g_{C i}(\boldsymbol{x}, t)-g_{C i}^{\mathrm{eq}}(\boldsymbol{x}, t)\right] \end{gathered} $

(5)

式中：g_Ci为表示浓度的分布函数；g_Ci^eq为浓度的平衡态分布函数，且

$ g_{C_i}^{\mathrm{eq}}(\boldsymbol{x}, t)=w_i C\left(1+\frac{\boldsymbol{u} \cdot \boldsymbol{c}_i}{c_{\mathrm{s}}^2}+\frac{\left(\boldsymbol{u} \cdot \boldsymbol{c}_i\right)^2}{2 c_{\mathrm{s}}^4}-\frac{\boldsymbol{u} \cdot \boldsymbol{u}}{2 c_{\mathrm{s}}^2}\right) $

(6)

无量纲松弛时间τ_gs=3D+0.5。

宏观浓度C可通过式(7)计算：

$ C=\sum\limits_{i=0}^8 g_{C i} $

(7)

1.3 边界条件

对于流动过程以及进口和出口，采用Zou-He压力边界指定压力，具体表达形式可参考文献[15]。采用反弹-镜面反射组合边界条件描述气固交界面速度滑移现象^[16]：

$ f_i(\boldsymbol{x}, t)=(1-r) f_i^{\mathrm{sr}}(\boldsymbol{x}, t)+r f_i^{\mathrm{bb}}(\boldsymbol{x}, t) $

(8)

其中，$ f_i^{\mathrm{sr}}(\boldsymbol{x}, t)$和$f_i^{\mathrm{bb}}(\boldsymbol{x}, t) $分别表示镜面反射方向和反弹方向的密度分布函数，r表示反弹比^[17]。

对于扩散-反应过程，一般性浓度边界条件可表示为

$ b_1 \frac{\partial C}{\partial \boldsymbol{n}}+b_2 C=b_3 $

(9)

式中：$ \partial C / \partial \boldsymbol{n}$为边界处法线方向上的浓度梯度，而参数b₁、b₂和b₃取值决定了边界条件类型。相关边界条件的详细描述可参考文献[11]。

基于赵建平等^[18]的热重分析实验，积碳催化剂在O₂/CO₂气氛中的再生反应是一级反应。因此，这里将消碳反应表达式和反应速率简化为：

$ \text { Reactant }_{(\text {gas })}+\text { Coke }_{(\text {solid })} \rightarrow \text { Product }_{(\text {gas })} $

(10)

$ R=k_{\mathrm{r}} \cdot C_{\mathrm{R}} $

(11)

式中：R为反应速率，k_r为反应常数，C_R为浓度，反应表面的浓度可通过下式计算：

$ C_{\mathrm{w}}=\frac{C_{\mathrm{f}}+0.5 \Delta x \boldsymbol{n} \cdot \boldsymbol{c}_i b_3 / b_1}{1.0+0.5 \Delta x \boldsymbol{n} \cdot \boldsymbol{c}_i b_2 / b_1} $

(12)

式中：C_f为边界处气相格子的浓度，n为反应表面法向矢量。

1.4 固相更新方法

采用固相更新方法实现消碳反应导致的孔隙结构的改变。气相和固相均采用控制体积V_m标记。起初，固相和气相控制体积分别标记为1和0。当V_m≤0时，固相格子转变为气相格子。边界固相格子的控制体积随消碳反应变化可表示为^[11]

$ V_{\mathrm{m}}(t+\Delta t)=V_{\mathrm{m}}(t)-\overline{V_{\mathrm{m}}} a_{\mathrm{m}} R \Delta t $

(13)

式中a_m、R和V_m分别为比表面积、反应速率和摩尔体积。

1.5 模型验证

微尺度流动-扩散和反应模型在文献[11]中已得到验证。为了验证固相更新方法，模拟了二维简单孔隙的溶解过程，如图 2所示。当达姆科勒数(Dam)取0.000 1，积碳壁面发生均匀溶解时，该问题存在解析解。无量纲渗透率和孔隙率遵循指数为3的函数关系^[19]。比较渗透率-孔隙度理论关系与LBM溶解模拟值发现，模拟得到的关系式与理论值相符合，从而验证了固相更新方法的有效性。

通过网格独立性分析，选用250×100格子数目作为最终的计算网格。基于开源数值计算软件GNU Octave进行孔隙尺度模拟。

2 结果与讨论 2.1 简单孔隙的结构和传递特性变化

采用无量纲参数(初始努森数Kn₀)表征孔隙结构特性。图 3给出了基准工况(Kn₀=0.1)条件下消碳过程的孔隙结构瞬时变化。由图 3可知，入口附近积碳首先消除，随着时间的推移，入口消碳显著，沿轴向孔隙的变化幅度增大。

为了表征消碳过程孔隙结构演变对传递和反应特性的影响，图 4给出了基准工况下速度和反应物分布的瞬时变化。

由图 4(a)可知，随着消碳反应进行，整体流动阻力减小，导致整体流动速度升高。图 4(b)给出了入口处径向速度分布。由图 4(b)可以观察到努森层作用下的速度滑移现象。随着时间的推移，孔径增大导致滑移速度不断减小，这说明消碳过程削弱了气体稀薄效应。图 4(c)给出了生成物浓度的瞬时变化。由图 4(c)可知，随着反应的进行，生成物浓度沿轴向逐渐上升。但随着时间推移，生成物浓度上升的速率反而有所降低。这主要归因于孔隙的变宽削弱了气体稀薄效应，消碳反应速率也随之下降。

2.2 初始努森数的影响

用初始努森数Kn₀表征孔径大小，分析其对孔隙结构演变和主相气体输运机制的影响。孔隙结构随不同初始努森数的变化如图 5(a)所示。由图 5(a)可知，当初始努森数＜0.01时。初始努森数的增大导致入口处消碳反应更加显著，孔隙结构更加不均匀。这归因于气体高稀薄效应下，努森扩散对主相气体输运的贡献愈发重要，强化了局部传质速率。当初始努森数＞0.01时，孔隙局部变宽，削弱了气体稀薄效应的影响，抑制了上述气体高稀薄效应下消碳反应的入口效应。上述分析表明，“大孔”和“小孔”内均匀分布的积碳，更容易均匀地消除，而“中孔”内均匀分布的积碳，受到入口效应的影响更加显著。因此，对不同孔径的孔隙结构应采用对应的消碳方式，以提高消碳效率。引入本征达西数Da_in、表观达西数Da_app和渗透比R_Kn来评价主相气体输运特性变化。将表征黏性流的本征渗透率κ_in采用孔径D_p的平方无量纲化得到本征达西数，将表征气体主相总输运(即黏性流动和努森扩散之和)的表观渗透率κ_app采用孔径的平方无量纲化得到表观达西数，将努森扩散等价的渗透率κ_Kn和本征渗透率κ_in之比定义为渗透比，表达式分别如下：

$ \begin{aligned} & D a_{\text {in }}=\frac{\kappa_{\text {in }}}{D_{\mathrm{p}}^2}, \\ & D a_{\text {ap }}=\frac{\kappa_{\text {app }}}{D_{\mathrm{p}}^2}, \\ & D a_{\kappa n}=\frac{\kappa_{K n}}{\kappa_{\mathrm{in}}} \end{aligned} $

(14)

由图 5(b)可知，随着消碳反应进行，黏性流动在不同的初始努森数下均大大增强，本征达西数随孔隙率增大而不断升高。对应孔隙结构的变化，随着初始努森数的增大，本征达西数升高的速率先减小后增大。从图 5(c)可知，初始阶段的表观达西数随着初始努森数的增大而增大。这归因于较高的初始努森数意味着努森扩散对主相气体输运的贡献更加显著。随着消碳反应进行，黏性流动和努森扩散均得到增强，因此表观达西数的升高速率明显快于本征达西数。类似于本征达西数，表观达西数的升高速率随初始努森数的变化同样存在临界阈值，表明“中孔”内积碳消除引起的传质速率提高受到一定的限制。由图 5(d)可知，不同初始努森数条件下渗透比均随着孔隙率的增大而减小。随着初始努森数的增大，渗透比的下降速率也显著提高，表明气体高稀薄效应下消碳反应导致孔隙结构演变对黏性流动的促进作用比努森扩散更加显著。

2.3 达姆科勒数的影响

孔隙结构演变特性与反应特性密切相关，达姆科勒数对消碳过程中孔隙结构演变机制和气体输运机制的影响见图 6。图 6(a)为同一孔隙率下，孔隙结构随不同达姆科勒数的变化情况。由图 6(a)可知，当达姆科勒数＜0.010时，积碳消除过程中积碳分布较为均匀；当达姆科勒数进一步增大，尤其是达姆科勒数为1.000时，入口效应将十分显著。表明当反应物具有高反应速率或低扩散速率时，孔隙深处的积碳更加不容易消除。图 6(b)给出了不同达姆科勒数条件下，渗透比随着孔隙率的变化。由图 6(b)可知，随着达姆科勒数的增大，渗透比下降的速率也明显加快，并且更早达到稳定状态。表明高反应速率条件下，消碳反应导致孔隙结构演变对黏性流动的促进作用比努森扩散更加显著。

2.4 复杂孔隙的结构和传递特性

图 7给出了复杂孔隙结构特性的瞬时变化。由图 7可知，随着消碳反应进行，可以观察到孔隙的变宽和合并，并且入口附近积碳首先消除，这与单孔相似。但在积碳消除过程中，孔径始终呈钟形分布。随着孔隙变宽和合并，导致孔径频率分布的峰右移和降低，这表明孔径分布变得更加均匀。复杂孔隙结构的演变也引起了流动特性的变化。从图 8可知，孔隙变宽和合并，大大减小了复杂孔隙内的气体流动阻力，导致流速升高。速度的频率分布也随之发生变化，低速度的频率降低而高速度的频率升高。

3 结论

基于格子玻尔兹曼方法以及固相更新方法，构建了消碳反应的孔隙尺度模型，开展了孔隙尺度消碳行为的数值模拟，主要结论如下：

1) 简单孔隙内积碳均匀分布条件下，入口附近积碳首先消除，形成沿轴向随时间延长变化幅度明显的孔隙结构。消碳过程减弱了气体稀薄效应，降低了传质速率和反应速率。

2)“大孔”和“小孔”简单孔隙内均匀分布的积碳，更容易均匀地消除，而“中孔”内积碳，入口效应更加显著，消碳反应引起的传质速率升高也受到了抑制。

3) 气体高稀薄效应和高反应速率条件下，消碳反应导致简单孔隙结构演变对黏性流动的促进作用比努森扩散更加显著。

4) 对于复杂的随机孔隙结构，同样存在入口效应，并且积碳消除会导致孔径频率分布峰值的右移和下降。