为实现我国“碳达峰，碳中和”的目标任务，发展新型高效的能源技术是必不可少的，而固体氧化物燃料电池(SOFC)作为电化学能量转换系统，因其高效和环境友好的特性而受到越来越多的关注^[1-3]。SOFC在高温下(873.2~1 273.2 K)条件下工作，可以实现碳氢燃料的直接内部重整，因此具有高效率、燃料灵活性和长期稳定性等优势，能够在目前从以化石燃料为基础的能源结构向大规模开发可再生燃料的能源结构的过渡阶段发挥重要作用^{[4, 5]}。目前，利用CO₂将CH₄转化成合成气在作为SOFC的燃料被广泛研究，因为CO₂和CH₄都是当前单一碳源化学物质中即廉价又无毒的重要组成成分^[6]。然而，CH₄热分解以及一氧化碳的歧化反应会造成碳沉积在多孔阳极中造成镍等非贵金属催化剂严重失活^{[7, 8]}。陶瓷LSCM是一种混合导体，并且已经证明对CH₄重整具有高转化效率、优异的合成气以及较好的抗结焦性^[9]。然而，LSCM电极的性能仍然受电催化活性不足的限制，因此研究人员通过将金属纳米粒子掺杂到LSCM的主晶格中，然后在还原条件下，以金属纳米颗粒的形式在LSCM骨架的表面溶出，这些直接从LSCM骨架载体上原位生长的金属纳米粒子，由于纳米尺度锚定界面结构导致其表现出了更强的电催化性、高温稳定性以及抗结焦性^{[10, 11]}。镍对碳氢化合物具有较高的催化活性，但同时又易于形成焦碳。相反，铜对碳氢化合物的催化活性较低，但具有较高的抗结焦性能^{[12, 13]}。因此，利用镍和铜的各自催化特性，将镍与铜形成的合金纳米颗粒(Ni_1-x-Cu_x)并使其在LSCM骨架原位上生长，以释放锚定和受限界面结构。同时对不同镍/铜比掺杂的LSCM进行了性能研究，研究结果表明，生长的Ni_0.5-Cu_0.5纳米粒子深深地锚定在LSCM基底上，两相之间存在明显的异质结。并验证了在其对CH₄/CO₂的高效电化学重整过程中，运行300 h后仍具有较高的催化活性、稳定性以及抗结焦性^{[14, 15]}。康丁文等利用柠檬酸-硝酸盐法制备了LSCM阳极材料，以及化学共沉淀法制备出了其负载产物Cu/Ni-LSCM，并对其负载产物在固定床反应器下进行了CH₄干重整的催化实验。研究结果表明，Cu的引入能够降低高温下金属相的烧结，降低活性颗粒的晶粒尺寸。除此之外，Cu的引入能够降低催化反应后碳纳米管以及活性颗粒的直径，增强气体分子与活性颗粒的接触几率，降低沉积碳^[16]。

数值模拟技术作为对SOFC研究切实可行的有效方法以及可视化的性能被广大研究人员用于SOFC的内部反应的研究。本论文利用数值模拟的方法建立了LSCM以及Cu/Ni-LSCM中CH₄与CO₂干重整动力学模型，模型耦合了动量传递、质量传递、化学反应动力学、域微分方程以及气体在多孔介质中的传质模型。利用此模型对碳沉积对阳极的影响做了细致的研究。

1 模型描述与计算方法

在本研究中，具有物质的量之比为1∶1的CO₂/CH₄被供应到SOFC多孔阳极通道用于探究不同n(Cu)/n(Ni)对碳沉积的影响。为探究固体氧化物燃料电池阳极材料Cu/Ni-LSCM中不同n(Cu)/n(Ni)对碳沉积的影响，如图 1所示搭建了SOFC的阳极区域，本研究采用有限元软件COMSOL建立三维模型，模型尺寸如表 1所示。

有限元分析在数值模拟的计算过程中是一个必不可少的环节。而模型中网格的划分直接影响计算结果的精度。若网格划分较小即产生的网格数量较多，计算的结果也较为精确，但是同时计算的时间也会越长；若网格划分较大，可能会导致计算结果与现实结果不符。因此找到适中的网格划分对建模求解极为重要。对三维模型进行了网格划分如图 2所示。

1.1 模型描述

本研究建立了2个模型研究了CH₄干重整过程中CH₄在SOFC多孔阳极中热分解为氢气和固体碳的过程。在第1个模型中，通过反应工程中的化学反应接口模拟理想反应器中的等温过程，分析了碳沉积对催化活性的影响。在第2个模型中，通过浓物质传递以及Brinkman方程研究了碳沉积对孔隙率以及流体流动的影响，其中包含了时间和空间的依赖性。

1.2 计算方法 1.2.1 电极内的化学反应

在多孔的阳极支撑层内CH₄和CO₂发生干重整反应方程(1)，以及CH₄的热分解反应方程(2)，反应速率以及反应常数如表 2所示，反应方程如下：

$ \mathrm{CH}_{4}+\mathrm{CO}_{2}=2 \mathrm{CO}+2 \mathrm{H}_{2} $

(1)

$ \mathrm{CH}_{4}=\mathrm{C}+2 \mathrm{H}_{2} $

(2)

1.2.2 多组分传质模型

本研究SOFC的燃料是CH₄与CO₂在孔阳极内干重整，有多种气体在多孔阳极内发生反应、对流以及扩散等各种复杂的过程，因此，混合物质质量传递模型被用来计算SOFC中各种不带电荷气态物质的浓度。物种i在稳态下的摩尔守恒方程如式(3)和式(4)。

在气体通道内：

$ \nabla \cdot j_{i}+\rho(u \cdot \nabla) \omega_{i}=0 $

(3)

在多孔电极内：

$ \nabla \cdot j_{i}+\rho(u \cdot \nabla) \omega_{i}=R_{i} $

(4)

式(3)和式(4)中：ρ是气体混合物浓度；j_i代表物种i的摩尔扩散通量；R_i是作为电化学、化学过程的结果的物质i的摩尔变化率；u代表气体混合物的摩尔平均速度。

在气体通道以及多孔电极内物种i的扩散通量由Stefan-Maxwell计算所得^[20]。在气体以及燃料通道内分子的扩散或连续扩散起主要作用^[21]，而在多孔电极中，需要考虑分子扩散、分子与固体壁面的碰撞以及粘性流动，因此在多孔电极内采用Stefan-Maxwell以及努森扩散的耦合，被用来描述多孔电极的多物质传输过程：

$ j_{i}=-\left(\rho D_{i, \mathrm{~m}}^{\mathrm{eff}} \nabla \omega_{i}+\rho D_{i, \mathrm{~m}}^{\mathrm{eff}} \omega_{i} \frac{\nabla M_{\mathrm{n}}}{M_{\mathrm{n}}}-j_{\mathrm{c}, i}+D_{i}^{\mathrm{T}} \frac{\nabla T}{T}\right) $

(5)

$ D_{i}^{\mathrm{m}}=\frac{1-\omega_{i}}{\sum\nolimits_{k \neq i} \frac{X_{\mathrm{k}}}{D_{i \mathrm{k}}}} $

(6)

$ D_{i \mathrm{k}}=\frac{0.0143 T^{1.75}}{2 P\left(v_{i}^{1 / 3}+v_{\mathrm{k}}^{1 / 3}\right)}\left(\frac{1}{M_{i}}+\frac{1}{M_{k}}\right)^{1 / 2} $

(7)

$ j_{\mathrm{c}, i}=\rho \omega_{i} \sum\limits_{k} \frac{M_{i}}{M_{n}} D_{k}^{\mathrm{m}} \nabla X_{\mathrm{k}} $

(8)

式(5)~式(8)中：D_ik是物种i和k的二元扩散系数^[22]；ω_i表示物质i的摩尔通量；v_i代表气体分子i扩散体积；D_{i, m}^eff是考虑分子扩散和努森扩散以及多孔电极微观结构的校正后物种i和k的有效二元扩散系数。努森扩散系数D_i^k可以根据动力学理论^[23]计算，如式(9)所示:

$ D_{i}^{\mathrm{k}}=\frac{2 r_{\mathrm{p}}}{3} \sqrt{\frac{8 R T}{{\rm{ \mathsf{ π} }} M_{i}}} $

(9)

所以在气体通道内：

$ D_{i, \mathrm{~m}}^{\mathrm{eff}}=D_{i}^{\mathrm{m}}=\frac{1-\omega_{i}}{\sum\nolimits_{k \neq i} \frac{X_{\mathrm{k}}}{D_{i \mathrm{k}}}} $

(10)

在多孔电极内：

$ D_{i, \mathrm{~m}}^{\mathrm{eff}}=\frac{\varepsilon}{\tau}\left(\frac{1}{D_{i}^{\mathrm{m}}}+\frac{1}{D_{i}^{\mathrm{k}}}\right) $

(11)

式(9)和式(11)中：r_p代表多孔电极的平均孔径；ε/τ是孔隙率与弯曲度的比值，用来修正多孔微观电极结构。

动量守恒：计算流体力学模型用于模拟通道和多孔电极中的气体流动。连续性方程和经典的Navier-Stokes方程应用于燃料气体通道和多孔电极。

则在气体通道内的质量平衡方程为：

$ \nabla \cdot(\rho u)=0 $

(12)

气体通道内则采用连续的Navier-Stokes动量守恒：

$ \begin{gathered} \rho(u \cdot \nabla) u=\nabla \cdot[-P I+\mu(\nabla u+ \\ \left.\left.(\nabla u)^{T}\right)-\frac{2}{3} \mu(\nabla \cdot u) I\right]+F \end{gathered} $

(13)

由于电极的多孔结构以及电极之间的质量交换，连续性方程中包含了一个源项，以考虑从多孔电极中的质量转移。Brinkman方程，即带有达西项的修正牛顿流体方程，被用来描述多孔电极中的动量守恒。

多孔电极内的质量平衡方程如式(14)：

$ \nabla \cdot(\rho u)=Q_{m} $

(14)

动量输运的控制方程是Navier-Stokes方程。该方程用于同时计算多孔材料、燃料通道和空气通道中的动量输运。Navier-Stokes方程可以写成如下形式^[24]：

$ \begin{gathered} \frac{1}{\varepsilon_{\mathrm{p}}} \rho(u \cdot \nabla) u \frac{1}{\varepsilon_{\mathrm{p}}}=\nabla \cdot\left[-P I+\mu \frac{1}{\varepsilon_{\mathrm{p}}}(\nabla u+\right. \\ \left.\left.(\nabla u)^{T}\right)-\frac{2}{3} \mu \frac{1}{\varepsilon_{\mathrm{p}}}(\nabla \cdot u) I\right]- \\ \left(\mu \kappa^{-1}+\beta_{\mathrm{F}}|u|+\frac{Q_{\mathrm{m}}}{\varepsilon_{\mathrm{p}}^{2}}\right) u+\boldsymbol{F} \end{gathered} $

(15)

式(15)中：F为体积力矢量；I为单位矩阵；u是速度矢量；ρ是混合物密度；质量源Q_m占单位体积中所有气体物种质量的总生产或消耗量；κ和ε_p分别表示电极的渗透性和多孔性；μ表示气体混合物的动力黏度。因为电化学反应只发生在催化剂层中，所以除催化剂层外质量源Q_m为0。流体混合物的动力黏度定义为^[25]：

$ \mu=\sum\limits_{i} \frac{\omega_{i} \mu_{i}}{\sum\nolimits_{i}\left(\omega_{i} \sqrt{\frac{M_{j}}{M_{i}}}\right)} $

(16)

式(16)中：ω_i是物质i的质量分数；μ_i是其黏度. 可以使用物种M_i的相对分子质量。

1.2.3 碳沉积计算

当燃料电池运行时，对于多孔阳极，碳沉积会覆盖于活性催化颗粒镍的表面，降低了催化剂的比表面积，从而导致催化活性的丧失；而从燃料气体传质的方面来说，沉积的碳会堵塞固体氧化物燃料电池多孔阳极的微孔，造成多孔阳极的孔隙率降低，传质受阻。因而对SOFC的瞬态研究建模非常重要。其中，总体的碳沉积速率可以表示为^{[26, 27]}：

$ r_{\mathrm{C}}=R_{1}-R_{2} $

(17)

式(17)中：r_C表示总体的碳沉积速率。当R₁>R₂时，r_C=R₁-R₂；当R₁≤R₂时，r_C=0。

变化的沉积碳浓度：

$ \frac{\mathrm{d} c_{\mathrm{C}}}{\mathrm{d} t}=r_{\mathrm{C}} \frac{\mathrm{d} \rho_{\mathrm{C}}}{\mathrm{d} t}=r_{\mathrm{C}} M_{\mathrm{C}} $

(18)

式(18)中：c_C分别代表阳极端沉积碳浓度；ρ_C为沉积碳的密度；r_C为总体的碳沉积速率；M_C为沉积碳的摩尔质量。为了模型更好地收敛，我们设置初始阳极端沉积碳浓度为极小值0.000 01。

变化的孔隙率：

$ \frac{\mathrm{d} \varepsilon}{\mathrm{d} t}=\frac{\varepsilon r_{\mathrm{C}} M_{\mathrm{C}}}{\rho_{\mathrm{C}}} $

(19)

式(19)中：ε、r_C、M_C、ρ_C分别代表阳极端的孔隙率、碳沉积速率、沉积碳摩尔质量以及阳极端沉积碳密度。阳极端孔隙率初始值为0.4。

变化的催化活性：

$ \frac{\mathrm{d} a}{\mathrm{~d} t}=-k_{\mathrm{a}} r_{\mathrm{C}}{ }^{2} c_{\mathrm{C}} a $

(20)

式(20)中，

$ k_{\mathrm{a}}=k_{\mathrm{a} 0} \cdot \exp \left(\frac{1.356 \times 10^{5}}{R_{\mathrm{g}} T}-32.007\right) $

(21)

式(21)中：k_a0与a分别代表催化剂衰减常数与催化活性。催化活性初始值设为1。

SOFC的物性参数(孔隙率ε、密度ρ)以及边界条件中反应体系的环境压力P、环境温度T_ref、初始反应温度T、入口压降d_p以及气体黏度μ等参数如表 3所示。

2 结果与分析 2.1 模型验证

为了验证上述所建立数学模型的可靠性，将数值计算的结果与实验的结果进行比较，图 3给出了电流密度随电压变化的发电性能曲线的数值计算结果和文献[28]中给出的实验测试结果对比图。计算时所采取的电池操作条件与实验测试一致。

由图 3可以看出，模型中的电流密度随电压变化的曲线能够与实验测试结果可以较好地吻合，说明数学模型具有可靠性，可以用来预测和分析电池的性能。

2.2 不同n(Cu)/n(Ni)下催化活性随时间的变化

图 4显示了当温度为1 073.2 K，入口压降为6 Pa时，不同n(Cu)/n(Ni)下催化活性随时间的变化曲线图，可以看出，当n(Cu)/n(Ni)为1∶1运行10 h后，催化剂活性为0.989 7。而当n(Cu)/n(Ni)为5∶1运行10 h时，催化剂活性仅为0.528。研究结果表明，Cu的引入，极大的提高了电池阳极的抗积碳性能，当n(Cu)/n(Ni)为1∶1时，抗积碳性能最优，运行10 h后催化活性依然保持在较好的状态。Lu等^[29]分别对LSCM、LSCM-Cu、LSCM-Ni_0.25Cu_0.75、LSCM-Ni_0.5Cu_0.5、LSCM-Ni_0.75Cu_0.25及LSCM-Ni进行了CH₄干重整实验，实验表明在LSCM-Ni_0.5Cu_0.5作为SOFC的阳极材料情况下，使得CH₄干重整的效率最高，积碳速率最小。所以在n(Cu)/n(Ni)=1∶1时展开其他工况的研究。

2.3 n(Cu)/n(Ni)=1∶1时阳极内各气体摩尔浓度随时间分布

图 5显示了阳极燃料气体CH₄/CO₂发生干重整以及CH₄分解时，各气体以及沉积碳的瞬态摩尔浓度，显示了存在催化剂失活和恒定催化剂活性时各物质浓度。结果表明，从浓度随时间变化来看，CH₄/CO₂干重整以及CH₄分解在催化剂活性不变时的反应速率略大于催化剂活性失活的反应速率。

而且在恒定催化活性的条件下，当燃料电池运行到37 h时，多孔阳极内发生CO的歧化反应，致使CO的浓度略低于H₂的浓度，而此时生成的CO₂和CH₄继续发生干重整反应，从而使得平衡时恒定催化活性条件下CH₄摩尔浓度低于催化失活条件下的CH₄摩尔浓度，与此同时CO₂的浓度在恒定催化剂活性以及催化失活的条件下完全相同，也证实了CO歧化反应的发生。研究结果表明，在Cu/Ni-LSCM作为多孔阳极材料的较高催化活性的条件下，CH₄分解反应被有效抑制，CO歧化反应都是造成碳沉积的主要原因，而且CO歧化反应只在较高催化活性条件下发生。

2.4 积碳反应的影响

本模型通过孔隙率的变化来反应积碳反应的大小，图 6为在阳极/气体通道界面处x=0、5、10和15 mm处时，孔隙率随时间的分布情况。可以看出，在x=0 mm随着时间的进行，孔隙率下降最快。而在x=10 mm和x=15 mm时，随着时间的进行，x=15 mm处的孔隙率略低于x=10 mm处的孔隙率。研究表明，入口处的碳沉积最为严重。而且碳沉积速率与时间成正比的关系。谢静等^[29]研究结果都表明直接CH₄内重整SOFC积碳活性a在整个阳极范围内均大于1，在燃料入口处约达到了270。本研究结果与谢静等^[29]研究结果都表明了CH₄干重整反应易在入口处发生积碳，与此同时CH₄热分解是阳极积碳的主要原因。

图 7显示了10 000和20 000 s孔隙率沿电池长度方向上的变化，以及阳极/电解质界面上的孔隙率分布情况。

由图 7可以看出，孔隙率在燃料入口处降低到最小值，然后随气流方向逐渐升高。这也说明了CH₄热分解反应速率在入口处最大，随着CH₄干重整反应以及CH₄热分解反应的进行，CH₄浓度逐渐降低，生成的CO和H₂浓度逐渐升高。H₂浓度的升高同时也抑制了CH₄热分解反应的正向进行。

因此沿着气流方向，一方面，在燃料入口处附近的CH₄含量很高，CH₄热分解速率以及CH₄干重整反应速率较高，生成的碳以及H₂和CO的量也随之增加，生成的碳造成了孔隙率的减小。另一方面，随着反应的不断深入，CH₄在反应中逐渐被消耗，碳沉积的量也随之减少。10 000 s时，入口端处孔隙率减少了0.023；20 000 s时，入口端处孔隙率减少0.045。研究表明，随着碳沉积的增大导致孔隙率的减小，从而使得CH₄各反应速率随之减小。

由图 7(e)可以看出电池运行10 000 s以及20 000 s时，阳极/电解质界面以及阳极/气体通道界面处沿电池长度方向的孔隙率分布。可以看出，在阳极/气体通道的入口端为碳沉积最严重的区域，在x=15.6 mm至气体出口处，电池运行10 000 s以及20 000 s时，阳极/气体通道孔隙率分布近似相同。可以看出，此时CH₄已经从阳极/气体通道界面处扩散到阳极/电解质界面处，使其CH₄浓度均匀，造成碳沉积的量也大致相同。

3 结论

该模型研究了LSCM以及Cu/Ni-LSCM阳极材料的抗积碳性能、催化活性以及孔隙率随时间的变化情况。并在此模型的基础上探究了不同燃料入口比对碳沉积的影响。通过与实验的验证与对比得出结论。

1) Cu的引入，极大地提高了电池阳极的抗积碳性能，当n(Cu)/n(Ni)为1∶1时，抗积碳性能最优，运行10 h后催化剂活性最高。

2) 在n(Cu)/n(Ni)=1∶1的条件下，当电池催化性能随时间下降时，CH₄热分解是造成阳极积碳的主要原因。当电池阳极催化活性不变时，电池运行37 h后，CO会发生歧化反应，同时也会造成阳极发生积碳。

3) 碳沉积最为严重的部分为燃料入口处，沿着气流方向碳沉积逐渐减弱，因此加强阳极燃料入口处的抗积碳性能尤为重要。