多相催化微观动力学与催化剂理性设计

朱贻安，周兴贵，袁渭康

（华东理工大学化学工程联合国家重点实验室，上海 200237）

多相催化通过经济、高效以及环境友好的方式，将原材料转化为高附加值的化学品和燃料，在化学、食品、药、汽车以及石油等传统工业中具有众多应用，在当今世界经济发展中扮演着越来越重要的角色。至少有90%的化学过程中需要通过多相催化剂来完成[1]，因此，多相催化是化学工业的基础和未来。同时，一些新兴的领域正由多相催化衍生出来，包括燃料电池、绿色化学、纳米以及生物技术等。因此，作为多相催化的核心，催化剂的理性筛选和设计一直都是化学家以及化学工程师追求的目标。在权衡生产成本和经济效益以确定操作条件的前提下，选择具有较高活性、选择性、稳定性和低成本的催化材料是催化剂筛选的方向。在对催化剂进行理性筛选前，首先需要清晰地认识催化反应。而对于催化反应的认识过程，实际上是对其中包含的物理化学问题的理解过程。在多相催化反应过程中，由于影响反应最终结果的因素众多，涉及到催化剂的组成、结构、温度、反应气氛及其组成、气氛的分压和载体性质等，更重要的是这些因素之间还存在交互作用，例如温度与催化剂颗粒尺寸之间的关系，反应气氛与催化剂颗粒微观结构之间的关系、载体与活性组分界面处由相互作用而产生的特殊催化性质等，因此，研究者在很长一段时间内对于催化过程的理解和认识无法深入，无法抓住催化现象的本质。在没有定量甚至定性原则指导的前提下，过去催化剂的筛选完全依赖于经验和反复尝试（trial and error），因此效率不高。例如20世纪初，德国巴斯夫公司在最终确定合成氨催化剂的工业配方前，据说尝试了2 500多种催化体系，进行了6 500余次实验[2,3]，而今天，随着基于表面化学的微观动力学的建立和发展，这种情况正在逐步得到改观。

1 微观动力学

在化学反应工程中，宏观动力学方程的建立过程首先是依据有限的光谱实验数据构建一条反应路径，并通过尝试假定总包反应中的某一步反应（大多数情况下为非基元步）作为整个过程的速率控制步，而其他反应步骤处于准平衡状态下，通过求解方程组获得速率表达式。然后，将速率方程与动力学实验数据进行拟合验证，从中选择出能够描述反应动力学行为的速率方程。宏观动力学既可以是排除了质量传递、热量传递以及动量传递等物理影响因素，从而研究反应速率与反应物浓度、温度、催化剂和溶剂种类之间关系的本征动力学（即揭示化学反应自身的规律），亦可以是在传质、传热、催化剂失活以及反应器稳定性等因素影响下的表观动力学。其主要任务就是通过建立半经验的总包反应速率方程，指导反应器的设计和优化。然而，以上述方式获得的速率方程反过来却不能推演出整个反应的反应机理。这是因为不同的反应机理有可能获得相同的反应速率方程。因此，宏观动力学对于反应器设计而言是直接而有效的工具，但是由于它无法提供任何与反应机理相关的有价值信息，所以对于催化剂筛选没有帮助。

与之相反，微观动力学试图通过来自实验或者理论计算的表面热力学和动力学数据，在分子水平上描述化学反应包含的所有基元步。在描述过程中，对于整个过程的反应速率控制步和最丰表面物种不作任何假定，而是通过微观动力学分析得到反应中间产物的表面覆盖率，从而获得基元反应速率，最终确定反应的主要反应路径（dominant reaction pathway）和速率控制步（rate-determining step）。需要强调的是，这里提到的用于催化剂筛选的微观动力学是指包含了催化反应中涉及的基本表面化学的动力学，需要能够在原子水平上对催化过程进行描述。它起源于上个世纪80年代后期[4]，是由Dumesic[5-9]以及Froment等[10,11]分别独立创建并逐渐发展起来的。微观动力学包含微观动力学分析（microkinetic analysis）和催化反应综合（catalytic reaction synthesis）两个部分。微观动力学分析通过研究一个催化循环中发生于催化剂表面的一系列基元步相互之间，及其与表面之间的关系来考察催化反应[5]。催化反应综合是将来源于实验和理论的表面化学信息结合起来，阐述如何通过调变催化剂、催化反应循环和反应条件，提高目标产物的产率[12]。例如，在过渡金属催化剂的筛选过程中，前者用于研究单金属表面上特定反应的反应机理，后者用于完成合金催化剂的筛选工作。

1.1 微观动力学分析

在采用微观动力学分析一个给定催化过程的反应机理时，重要的一步是通过Arrhenius方程来估算基元反应的速率常数。而速率常数估算的前提是采用实验的或者计算的方法获取基元反应的指前因子（preexponential factor）和活化能。对于指前因子而言，通常采用碰撞理论（collision theory）或者过渡态理论（transition state theory）来获得其数量级。由于碰撞理论无法包含分子构型在估算中产生的影响，因此过渡态理论无论从合理性还是准确性的角度都更加适用。过渡态理论的基本假设是反应物与过渡态构型中的活性络合物之间存在准平衡关系，因此由反应物获得活性络合物的准平衡常数（K≠）可以定义为：

其中，ΔG°≠，ΔS°≠和ΔH°≠分别是由反应物生成活性络合物的标准自由能变、标准熵变和焓变，kB为Boltzmann常数，T为反应温度。在过渡态理论中，化学反应速率表示为活性络合物的浓度与频率因子（kBT/h）之间的乘积，其中h为Planck常数，因此化学反应速率常数（k）可以表示为：

其中，指前因子A为：

在过去相当长的一段时间内，由于无法合理地估算出上述标准自由能变、标准熵变和焓变，实验科学家往往是先依据化学反应速率和反应活化能获得某些化学过程的指前因子，然后将这些结果用以经验性地估算相似反应的指前因子。例如，分子吸附的指前因子为10～103Pa/s，而基于Langmuir-Hinshelwood机理发生的表面化学反应的指前因子为108～1013s-1。如果没有可以参考的经验性数据，经常采用的方法是在忽略熵变影响的前提下，采用频率因子对于不同温度下的指前因子进行估算。相对于指前因子，实验方法对于基元反应活化能的估算显得更加困难，主要遵循以下3步：首先，假设各基元步在气相中进行，估算其反应热。该估算主要依赖于实验测量获得的键能数据以及分子、自由基和离子的形成能；随后，通过引入气相物质的吸附热，将气相反应热转化为表面反应热。这一步中包含的表面成键和断键能来源于程序升温脱附和微量热实验等。事实上，这些键能是动力学模型中十分重要的参数；最后，依据基元反应活化能与其反应热之间的经验线性关系估算活化能。由于过去采用已有的实验数据进行反应速率常数的估算是唯一可行的手段，因此微观动力学分析被视为实验科学家的一种动力学研究工具。

从20世纪90年代开始，随着高性能计算机群的快速发展以及量子化学方法的普遍运用，采用团簇或者周期性模型获得与表面和吸附质之间相互作用相关的能量和几何结构已经成为了可能。一方面，超级计算机的硬件处理能力已经能够成功地处理复杂体系；另一方面，基于密度泛函理论（density functional theory，DFT）的第一性原理计算的计算精度足以描绘过渡金属及其合金对于特定化学反应的活性和选择性的变化趋势，已经成为描述多相催化过程的一个强有力工具。依据能量最小化原则进行构型优化以及采用NEB（nudged elastic band）等方法搜寻基元反应的过渡态，DFT计算不仅能够直观地描述吸附构型和过渡态结构，更加重要的是它能够提供能量和频率信息，结合统计热力学和过渡态理论，计算出基元反应的指前因子和活化能。然后，通过求解微观动力学模型获得各种反应中间体的表面覆盖率和各基元反应的反应速率，并最终获取反应的主要反应路径和相应的速率控制步，从而确立完整的反应机理，并为催化反应综合打下基础。同时，在采用微观动力学分析了解了反应机理的基础上，依据例如Langmuir-Hinshelwood或者Eley-Rideal机理，亦可写出总包反应的动力学方程，从而用于反应器的设计和优化。目前，通过DFT计算获得的反应热和活化能相较于实验结果，其误差为0.1～0.2 eV[13]。虽然这一误差导致定量重现实验中的反应速率还存在一定困难，但是这个计算精度足以预测不同反应体系中反应速率的变化，而且只要抓住了催化反应中关键表面化学的基元反应步，采用较为准确但可能更加耗时的计算方法或者结合一些已知的表面热力学和表面动力学实验数据对计算结果进行校正，就有可能半定量地获得复杂体系的动力学。

微观动力学分析不仅能够计算出反应速率，而且可以抽象出标识催化反应活性的描述符（descriptor）并确定其最优范围。对于一个基元反应来讲，在忽略熵贡献的前提下，其反应活化能（Eact）与反应热（ΔH）之间存在线性关系，即Brønsted-Evans-Polanyi（BEP）关系[14,15]：

其中，α是一个介于0和1之间的系数，E0是常数。换句话讲，如果忽略熵的影响，那么当基元反应的反应放热量增大或者吸热量减小时，其活化能会降低，反之亦然。因此，BEP关系构架起了动力学和热力学之间的桥梁，常常被用来解释实验中观察到的催化活性变化规律以及估算反应活化能。由于过去采用实验方法准确地测量基元反应的活化能和反应热比较困难，因此，缺少对于这一线性关系的直接和定量的证据支持。直到最近10年，一系列针对过渡金属表面上化学反应的DFT计算工作证明了BEP关系的正确性与可靠性[16-18]。

气固相催化过程的反应机理必然包含反应物的吸附以及产物的脱附。倘若催化表面与吸附质之间的相互作用过弱（解离吸附放热量低或者吸热量高），那么依据BEP关系，反应物的解离活化能将过高而降低总包反应速率；反之，倘若催化表面与吸附质之间的相互作用过强，则反应物的解离反应活化能会很低，然而此时脱附反应的活化能将过高从而降低脱附速率，最终造成表面活性位数量的减少而降低总包反应速率。因此对于一个特定的催化反应来讲，最优的催化剂应该与反应物、中间体和产物之间存在强度适中的相互作用，此即为Sabatier规则[19]。其定量表现为不同催化剂表面上的总包反应转化频率（turnover frequency，TOF）与速率控制步的反应热（抑或其它相关物理量）之间呈现火山形曲线（volcano curve）关系。通过火山形曲线，研究者能够进一步确定标识催化反应活性的描述符的最佳范围，从而完成微观动力学分析。

1.2 催化反应综合

依据微观动力学分析获得的信息，丹麦技术大学（现斯坦福大学）的Nørskov研究组在过渡金属催化剂筛选方面进行了大量研究工作[20-27]。研究组致力于采用基于标识化学反应速率的描述符进行催化剂筛选，该描述符通常是前面提到的控速步反应热或者简单吸附物种的吸附热。具体的策略是：首先，通过基于DFT计算的微观动力学分析或者实验研究获得不同单金属表面的化学反应转化频率以及速率控制步信息，而后尝试关联转化频率与描述符，并通过火山形曲线确定该描述符的最优取值范围；然后，通过计算一系列合金表面上相对应的描述符，在动力学分析确定的最优范围内筛选出具有高催化活性的合金催化剂；最后，运用DFT计算和微观动力学模拟对反应体系在该合金表面上的动力学进行系统计算，并通过合成相应的催化剂进行实验验证。后面两步为催化反应综合，是在微观动力学分析基础上筛选合金催化剂的过程。

如前所述，速率控制步的反应热通常作为标识反应速率的描述符。反应热的变化反映了吸附质与催化剂表面之间结合能的变化。而化学键的属性与强弱本质上与过渡金属表面不同的电子结构相关。对于后者，由Hammer等[20-22]提出的d带模型（d-band model）可知，在吸附质与过渡金属之间的电子耦合矩阵元（coupling matrix element）相近的前提下，吸附质分子吸附能的变化同过渡金属表面d电子带能量（d-band center）之间呈线性关系：

其中，Ed-hyb为过渡金属的d带与吸附质电子态之间的相互作用，C( fa,fd)与吸附质电子态和过渡金属d带的占据数相关，Vad为吸附质与过渡金属之间的电子耦合矩阵元，aε为吸附质电子态的能量，dε为过渡金属d带的能量，即通常所说的d带中心。在线性Muffin-Tin轨道（LMTO）理论中，Wigner-Seitz半径是影响d电子带能量的一个重要因素。例如，当原子的Wigner-Seitz半径变大时，来自相邻原子的嵌入电子密度降低，从而导致d带带宽减小。这就使得必须提高d带电子的能量从而保证d带的电子占据保持不变。随着d带电子能量的升高，由吸附质电子态和过渡金属d带之间杂化形成的反键态的能量将升高，使得至少部分电子态的能量将由低于费米能级变为高于费米能级，从而不能被电子所占据，这就加强了吸附质与过渡金属表面之间的相互作用。换句话讲，d带能量的迁移将直接调变催化剂的吸附能力，是影响宏观动力学的重要微观性质。因此，d带模型的贡献在于能够从电子结构的角度解释过渡金属表面吸附以及催化性质的变化规律，显著减少计算工作量，从而使得大量催化体系的研究成为可能。

对于合成氨反应的合金催化剂筛选是催化反应综合的一个应用实例[23-25]。研究者采用基于DFT的第一性原理计算和统计热力学计算，在过渡态理论框架下获得反应速率常数并求解微观动力学模型以获得不同单金属台阶表面上合成氨反应的TOF，并将其与N原子的吸附能进行关联。结果发现，两者之间呈现很好的火山形曲线关系，随着N原子吸附能（绝对值）的逐渐降低，金属对于合成氨反应的活性先升高后降低。其中，Ru和Os处于火山口，是最好的合成氨单金属催化剂。然而，这两种催化剂均为贵金属材料，工业应用成本过高，因此，目前工业生产中退而采用活性第三的Fe基催化剂。之所以存在上述火山形曲线关系，是由于N2分子的解离是合成氨反应的速率控制步骤，而且这个基元反应的过渡态在势能面上接近反应终态（N原子吸附态），因此过渡态构型的能量在不同金属表面上随着N原子吸附能的变化而变化，从而使得N2分子解离反应活化能与N原子的吸附能之间呈现线性关系。一方面，该线性关系表明N原子吸附能（绝对值）的升高会降低N2分子解离反应活化能，有利于提高反应速率；另一方面，解离反应活化能的持续降低会导致表面被过多的N原子所覆盖，表面空位的减少反而会降低反应速率。因此，N原子的吸附能存在一个使得合成氨反应速率最高的最优值。通过以上微观动力学分析，不仅获得了合成氨反应的机理，并且发现N原子吸附能可以作为标识催化剂对于合成氨反应活性的描述符，那么依据这个原则就可以通过催化反应综合来筛选合金催化剂。相对于Ru金属催化剂，N原子在Mo金属表面的吸附能过高，而在Co金属表面上的吸附能过低。将这两种元素的原子按照1∶1混合后，计算结果表明，Co-Mo合金拥有适中的N原子吸附能，恰好位于接近火山形曲线顶峰的位置。后续的实验结果证实，Co-Mo合金确实具有比Ru更好的合成氨反应催化性能，是替代Fe基催化剂的最佳候选。

另一个具有代表性的工作是对于甲烷化反应合金催化剂的筛选[26,27]。在该反应中，研究者首先通过DFT计算发现CO的解离反应在密堆积面上的活化能很高，而台阶和边缘可以显著地提高解离反应速率，是反应的活性中心。计算结果还显示CO解离吸附反应的过渡态构型中C-O键长被拉伸得较长，因此，与反应的终态构型（即C和O原子共吸附于表面）相似。在一部分金属表面上，CO分子的解离吸附控制总包反应速率，由于其反应活化能与反应热之间呈BEP线性关系，因此CO解离吸附热可以作为标识总包反应速率的描述符。当CO解离吸附热（绝对值）随着金属表面的改变逐渐增大时，CH4和H2O生成反应中初态由于含有C和O原子而变得更加稳定，因此反应活化能将逐渐升高，并会最终取代CO解离反应成为速率控制步。由于CH4和H2O生成反应的活化能与CO解离吸附热之间亦存在线性关系，因此，后者能够在所有金属表面上标识化学反应速率。通过将实验中获得的转化频率与CO解离吸附热关联，研究者获得了火山形曲线，并且发现最高反应速率对应的CO解离吸附热在0.06 eV左右。在此基础上，研究者首先筛选出117种能够在甲烷化反应条件下稳定存在的合金催化剂，然后计算出相应的CO解离吸附热，发现Ni3Fe对于CO的解离吸附热同位于火山形曲线顶端的Ru以及Co近乎相等。针对该体系的进一步的系统计算和实验都证明Ni3Fe具有非常好的甲烷化反应催化活性，同时价格十分低廉，已经被丹麦的Haldor Topsøe公司转化为实用催化剂[28]。

2 展 望

采用基于DFT计算和微观动力学模拟研究多相催化反应的反应机理，进而依据反应特性理性筛选过渡金属催化剂是一个切实可行的研究方向。由复杂的多相催化过程中剥离出标识催化活性的描述符并确定其最优范围，使得计算科学家有可能在比以往多几个数量级的催化材料中筛选催化剂，同时显著地减少了实验工作量。文中提到的两个成功案例已经充分展现了该方法的优越性。

然而即便如此，采用该方法进行催化剂筛选依然面临一些问题。首先，DFT计算的精度依赖于其交换-相关泛函（exchange-correlation functional）的精度。过去交换-相关泛函的发展主要聚焦于如何合理地描述分子的构型以及解离能，而对于决定化学反应动力学的反应能垒以及弱相互作用的描述（特别是van der Waals相互作用）却不甚精确。另外，由于DFT用于计算电子-电子排斥作用的交换-相关泛函（甚至DFT本身）是在单粒子近似的基础上发展起来的，因此，传统的电子结构计算方法不能很好地描述强电子相关（strongly correlated）体系，例如过渡金属氧化物或氮化物，这就局限了基于DFT计算的微观动力学模拟的应用范围。其次，采用多金属合金催化剂代替贵金属催化剂是一个发展方向，但是这种替换增加了反应系统的复杂性。如何在工业生产中保持催化剂活性，金属颗粒的热稳定性和催化稳定性是将要面临的一个课题。

尽管存在上述问题，但是应该相信随着高性能计算机、电子结构理论的进一步发展和催化剂材料制备水平的提高，在不久的将来，完全采用理论计算发现和优化催化剂将不再是遥不可及的目标。

[1] National Research Council. Catalysis looks to the future: panel on new directions in catalytic sciences and technology[M]. Washington D C: National Academy Press, 1992:1.

[2] Tamaru K. Catalytic ammonia synthesis: fundamentals and practice[C] // Jennings J R, ed. New York: Plenum Press, 1991:4.

[3] Mittasch A. Early studies of multicomponent catalysts[J]. Advances in Catalysis, 1950, 12: 81-104.

[4] Boudart M. From the century of the rate equation to the century of the rate constants: a revolution in catalytic kinetics and assisted catalyst design[J]. Catalysis Letters, 2000, 65(1/3): 1-3.

[5] Dumesic J A, Trevino A A, Milligan B A, et al. A kinetic modeling approach to the design of catalysts: formulation of a catalyst design advisory program[J]. Industrial & Engineering Chemistry Research, 1987, 26(7): 1399-1407.

[6] Dumesic J A, Rudd D F, Aparicio L M, et al. The microkinetics of heterogeneous catalysis[M]. Washington D C: American Chemical Society, 1993. 315

[7] Yaluris G, Rekoske J E, Aparicio L M, et al. Isobutane cracking over Y-zeolites II: catalytic cycles and reaction selectivity[J]. Journal of Catalysis, 1995, 153(1): 65-75.

[8] Cortright R D, Dumesic J A, Madon R J. Catalytic cycles for hydrocarbon cracking on zeolites[J]. Topics in Catalysis, 1997, 4(1/2):15-26.

[9] Cortright R D, Dumesic J A. Kinetics of heterogeneous catalytic reactions: analysis of reaction schemes[M]//Gates B, knoezinger H.Advances in Catalysis, vol. 46. USA: Academic Press, 2002: 161-264.

[10] Baltanas M A, Van Raemdonck K K, Froment G F, et al. Fundamental kinetic modeling of hydroisomerization and hydrocracking onnoble metal-loaded faujasites 1: rate parameters for hydroisomerization[J]. Industrial & Engineering Chemistry Research, 1989, 28(7):899-910.

[11] Froment G F. Kinetic modeling of hydrocarbon processing and the effect of catalyst deactivation by coke formation[J]. Catalysis Reviews-Science and Engineering, 2008, 50(1): 1-18.

[12] Rudd D F, Dumesic J A. Catalyst synthesis by analogy[J]. Catalysis Today, 1991, 10(2): 147-165.

[13] Hammer B, Hansen L B, Norskov J K. Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerh of functionals[J]. Physical Review B, 1999, 59(11): 7413-7421.

[14] Bronsted J N. Acid and basic catalysis[J]. Chemical Reviews, 1928, 5(3): 231-338.

[15] Evans M G, Polanyi M. Inertia and driving force of chemical reactions[J]. Transactions of the Faraday Society, 1938, 34: 11-24.

[16] Pallassana V, Neurock M. Electronic factors governing ethylene hydrogenation and dehydrogenation activity of pseudomorphic Pd(ML)/Re(0001), Pd(ML)/Ru(0001), Pd(111), and Pd(ML)/Au(111) surfaces[J]. Journal of Catalysis, 2000, 191(2): 301-317.

[17] Liu Z P, Hu P. General trends in CO dissociation on transition metal surfaces[J]. Journal of Chemical Physics, 2001, 114(19): 8244-8247.

[18] Nørskov J K, Bligaard T, Logadottir A, et al. Universality in heterogeneous catalysis[J]. Journal of Catalysis, 2002, 209(2): 275-278.

[19] Sabatier P. Hydrogénations et déshydrogénations par catalyse[J]. Berichte der Deutschen Chemischen Gesellschaft, 1911, 44(3):1984-2001.

[20] Hammer B, Norskov J K. Why gold is the noblest of all the metals[J]. Nature, 1995, 376(6537): 238-240.

[21] Ruban A, Hammer B, Stoltze P, et al. Surface electronic structure and reactivity of transition and noble metals[J]. Journal of Molecular Catalysis A: Chemical, 1997, 115(3): 421-429.

[22] Hammer B, Norskov J K. Theoretical surface science and catalysis: calculations and concepts[J]. Advances in Catalysis, 2000, 45:71-129.

[23] Logadottir A, Rod T H, Norskov J K, et al. The Bronsted-Evans-Polanyi relation and the volcano plot for ammonia synthesis over transition metal catalysts[J]. Journal of Catalysis, 2001, 197(2): 229-231.

[24] Jacobsen C J H, Dahl S, Clausen B S, et al. Catalyst design by interpolation in the periodic table: bimetallic ammonia synthesis catalysts[J]. Journal of the American Chemical Society, 2001, 123(34): 8404-8405.

[25] Honkala K, Hellman A, Remediakis I N, et al. Ammonia synthesis from first-principles calculations[J]. Science, 2005, 307(5709):555-558.

[26] Bligaard T, Norskov J K, Dahl S, et al. The Bronsted-Evans-Polanyi relation and the volcano curve in heterogeneous catalysis[J]. Journal of Catalysis, 2004, 224(1): 206-217.

[27] Andersson M P, Bligaard T, Kustov A, et al. Toward computational screening in heterogeneous catalysis: pareto-optimal methanation catalysts[J]. Journal of Catalysis, 2006, 239(2): 501-506.

[28] Norskov J K, Bligaard T, Rossmeisl J, et al. Towards the computational design of solid catalysts[J]. Nat Chem, 2009, 1(1): 37-46.