航空发动机和汽车发动机的动力都源自于燃料的燃烧.近年来, 因为能源安全、环境保护和气候变化等问题日益严重, 越来越多的研究者致力于提高发动机燃烧性能的研究.如今随着高性能计算技术, 复杂湍流流动数值模拟以及燃烧化学反应机理等方面的蓬勃进步, 计算机模拟技术正发展成为一个强大且高效的发动机辅助设计工具.其中, 研究燃烧化学反机理的目的是为了精确地预测和控制燃烧过程, 从而提高燃烧效率, 控制火焰稳定性以及优化排放.由于目前在燃烧反应模型中, 多采用简单的一步或多步总包反应机理, 所以其计算结果通常是有缺陷的, 因而须进行详细的、全面的燃烧化学反应机理研究.
对于在航空发动机中广泛使用的航空煤油, 关于其燃烧反应动力学的研究存在较大困难, 主要原因是它的化学成分非常复杂, 通常是由几百种碳氢燃料分子组成, 主要包括烷烃、环烷烃、烯烃和芳香族化合物等[1-2].因此通常的办法是将关于航空煤油的研究简化成由少数几种典型的碳氢燃料组成的替代模型的研究[3-6].此外, 对于目前备受研究者关注的生物燃料, 其中生物柴油凭借其优良的环保性, 可再生性和燃烧性能受到了研究者的青睐[7-11].更值得注意的是, 生物柴油良好的经济适用性使得它可以取代常规柴油或直接与常规柴油混合使用, 并且对柴油发动机的的特性及结构无额外要求[12-14], 也可以与航空煤油混合应用于航空发动机[15-16].
众多研究表明, 含10个碳原子以上的大分子燃料是组建航空煤油和生物柴油替代模型中最常用的组成部分, 为了更好地将航空煤油和生物柴油应用到发动机领域, 对其进行详细的燃烧化学动力学研究是十分必要的.目前相关的研究工作也在积极开展[17-20], 其目的是建立准确详细的关于航空煤油和生物柴油的反应机理[21-25], 用于以其为燃料的基于CFD的计算机辅助设计.
目前, 大分子燃料(航空煤油和生物柴油等)反应机理的发展遇到了瓶颈.因为大分子燃料普遍具有长碳链(一般含有12~20个碳原子)混合物属性、分子结构和物理化学性质的多样性, 所以之前大多数研究仅仅关注了含有较短碳链的典型分子.而相比于主要由碳氢燃料组成的航空煤油, 生物柴油的组分就显得更加复杂.因为一般采用植物油或动物脂肪等作为生物柴油的生产原料, 所以其主要成分由脂肪酸甲酯等含氧分子组成[26-27].
以生物柴油为例, 丁酸甲酯(MB, C3H7COOCH3)[28-32]和癸酸甲酯(MD, C9H19 COOCH3)[33-34]等具有代表性的替代燃料被用来模拟真实生物柴油的燃烧特性.其次, 建立一个关于替代燃料的详细反应机理通常需要为几百个化学组分指定热化学数据, 为上千个基元反应指定反应速率常数.但是, 应该注意的是, 对于如此大量的反应, 尤其是对那些在燃烧化学中十分重要而又难以通过实验来研究的反应和条件, 指定准确的与温度(压力)相关的反应速率常数是不现实的.
近年来, 理论化学动力学和电子结构理论的进步, 使得从头算化学动力学逐步发展成一个定量可靠的研究工具.对于一些相对较小的分子, 其反应速率常数的计算精度已经可以与实验相媲美[35-36].例如,沿反应坐标轴不同能垒的氢提取反应(RH+H→R+H2)的高压下速率常数能很好地用传统过渡态理论[37-38]定义
$ k(T) = \frac{{{k_{\rm{B}}}T}}{h}\frac{{{Q^† }}}{{{Q_{{\rm{RH}}}}{Q_{\rm{H}}}}}{e^{\frac{{ - V† }}{{{k_{\rm{B}}}T}}}} $ | (1) |
其中, V†表示能垒, Q†表示包括振动、转动及电子等因素的过渡态的配分函数, QRH, QH表示反应物单位体积的配分函数. T, kB和h分别表示温度, Boltzmann常数和Planck常数.能垒的高度和配分函数通过RH+R势能面上的电子结构计算得到.由式(1)可知, 理论上速率常数计算的不确定性主要依赖于能垒计算的准确性.例如, 当能垒(V†)被低估2 kcal/mol, 反应速率常数在典型的点火温度(500 K)下, 会被高估7.5倍, 在典型的燃烧温度(1 500 K)会被高估两倍[39].除此之外, 低温时的量子隧道效应[38-40]和在大分子体系中对扭转非简谐振动的处理[41-43]也会影响速率常数计算的不确定性.如果将存在较大不确定性的速率常数应用于反应机理, 会导致模型预测的燃烧参数(如层流火焰速度和点火延迟时间)大幅度偏离实验测量值.根据以上分析可知, 发展高精度的量化计算方法对于高精度反应动力学的研究是十分重要的.
如图 1所示, 在量子化学计算中, 方法与基组的配合决定了计算精度.我们可以把这种配合关系类比于流体力学计算中计算精度(更高阶)与网格质量(更密集)的关系.量子化学中的方法类比于计算精度, 基组类比于网格, 为了获得相对更精确的计算结果, 往往需要采用更高阶的计算精度和更密集的网格, 与此同时计算成本也会大幅增加.
![]() |
图 1 量子化学计算中关于方法配合基组的选择与计算精度关系的示意图 Fig.1 Illustration of the correlation between the method with basis set and calculated accuracy in quantum chemistry calculation |
目前,已有了一些高精度理论方法而且能行之有效地应用于研究一些燃烧化学关心的有机分子[40].比如单双取代包括非迭代三重激发的耦合族理论(the coupled cluster theory with single and double excitations and a quasi-perturbative treatment of connected triple excitations, CCSD(T))方法, 配合完全基组(complete basis set, CBS), CCSD(T)/CBS, 对于能垒和反应热的计算不确定性一般小于1.1 kcal/mol[44].单双取代包括三重激发的二次组态相互作用(the quadratic configuration interaction with singles, doubles and perturbative inclusion of triples,QCISD(T))方法配合完全基组(CBS), QCISD(T)/CBS,对于能垒和反应热的计算不确定性大约在1.0 kcal/mol左右[45], 如果考虑基于键可加性的校正格式, 计算不确定性可降低到0.6 kcal/mol[46], 但是由于计算量的问题, 这些方法都不适用于超过10个非氢原子的体系[40].因此, 发展可适用于大分子燃料体系的高精度量化计算方法不仅是必要的, 而且是急需的.
目前可适用于大分子体系的高精度量化计算方法有Zhang等[39]发展的一种两层的N层集成分子轨道和分子力学(our own N-layered integrated molecular orbital and molecular mechanics, ONIOM)方法, 可表示为ONIOM[QCISD(T)/CBS : DFT], Li等[47], Wang等[48]发展的一种基于能量的分块(generalized energy-based fragmentation, GEBF)方法和近期由Wu等[49-50]发展的基于ONIOM的基团加和法(cascaded group-additivity ONIOM, CGA-ONIOM)方法.
对于CGA-ONIOM方法, 目前还没有在实际的化学反应体系中进行可行性研究, 所以在此不做更多评价.而对于GEBF方法, Chi等采用该方法对生物柴油中的饱和甲酯体系进行了全面的研究, 并证明了该方法的可行性及计算精度[51-52].相比之下, Zhang等发展的两层ONIOM方法也已在饱和烷基酯体系和不饱和甲酯体系中得到了验证[39, 53-54], 其中包括饱和烷基酯与氢自由基的氢提取反应, 不饱和甲酯与氢自由基的氢提取反应, 氢加成反应, 异构化(分子内氢转移)反应和β位裂解反应.更主要的原因是相比于GEBF方法, ONIOM[QCISD(T)/CBS : DFT]方法在保证计算精度(对于处理类似反应体系, 这两种方法相对于QCISD(T)/CBS的计算误差基本都小于0.1 kcal/mol), 而且计算成本相比于QCISD(T)/CBS基本都可以节约90%的同时, 使用相对更为简单方便且可直接在Gaussian 09程序包中使用[55].
对于当前典型的大分子燃料(航空煤油和生物柴油等)来说, 由于大分子生物柴油的主要成分(脂肪酸甲酯等)相比于航空煤油中的主要成分(脂肪烃等), 其化学性质更为复杂, 主要原因在于生物柴油中不仅包括不饱和烃类的官能团, 如碳碳双键(C=C), 而且还包括多种含氧官能团, 如羟基(-OH), 羧基(-COOH), 羰基(C=O)等[56].所以, 能够处理生物柴油中含氧大分子燃料体系的高精度量化方法, 对于航空煤油中的碳氢大分子燃料也同样适用.因此, 本文主要对Zhang等[39]发展的ONIOM[QCISD(T)/CBS : DFT]方法关于处理生物柴油中的含氧大分子燃料(饱和和不饱和脂肪酸甲酯体系)进行了详细的介绍和评价.
1 ONIOM计算方法 1.1 结构优化计算反应势能面上固定点的几何结构和振动频率通过密度泛函理论Becke三参数泛函和Lee-Yang-Parr关联泛函(B3LYP), 配合6-311++G(d, p)基组计算得到[57-58].该方法由于计算准确性与经济性被广泛应用于燃烧化学动力学研究中.每个鞍点到相应的局部最小值还需使用内禀反应坐标理论对所有反应的过渡态进行验证.零点能校正由B3LYP/6-311++G(d, p)振动频率得到.
1.2 基准高精度单点能计算: QCISD(T)/CBS对于相对较小的分子体系, 其计算量(主要由非氢原子数量决定)也相对较低, 因此两种高精度QCISD(T)/CBS方法, 即[QCISD(T)/CBS]1和[QCISD(T)/CBS]2, 都可以用于计算, 其结果作为验证ONIOM方法正确与否和其精度的基准.第1种方法: [QCISD(T)/CBS]1或记为[QCISD(T)/CBS]TZ→QZ, 采用Dunning[59]的相关一致基组, cc-pVXZ, 其中X=D, T和Q, 这些基组分别简写为DZ(double-zeta), TZ(triple-zeta)和QZ(quadruple-zeta).由使用这些基组计算的结果可以外推得到完全基组(CBS)极限[60], 具体表示如下
E[QCISD(T)/CBS]1=E[QCISD(T)/CBS]TZ→QZ=E[QCISD(T)/QZ]+{E[QCISD(T)/QZ]-E[QCISD(T)/TZ]}×0.693 8
但是, 当研究的反应体系中非氢原子数超过6时, 因为这种方法的计算密集度过高, 所以计算量会大幅度增加, 此时可采用第2种方法[QCISD(T)/CBS]2[36, 39, 61]进行计算.该方法使用了2阶Møller-Ple-sset perturbation theory (MP2)方法配合cc-pVTZ和cc-pVQZ基组, 并根据该基组的计算结果外推得到完全基组(CBS)极限.一般情况下, 非氢原子数不超过9的反应体系, 都可以使用该方法处理, 具体表示如下
E[QCISD(T)/CBS]2=E[QCISD(T)/CBS]DZ→TZ+{[MP2/CBS]TZ→QZ-E[MP2/CBS]DZ→TZ}
其中
E[QCISD(T)/CBS]DZ→TZ=E[QCISD(T)/TZ]+{E[QCISD(T)/TZ]-E[QCISD(T)/DZ]}×0.462 9
E[MP2/CBS]TZ→QZ=E[MP2/QZ]+{E[MP2/QZ]-E[MP2/TZ]}×0.693 8
E[MP2/CBS]DZ→TZ=E[MP2/TZ]+{E[MP2/TZ]-E[MP2/DZ]}×0.462 9
[QCISD(T)/CBS]2方法避免了使用计算成本较高的QCISD(T)/QZ, 因而大大减小了计算量.该方法已在丁醇(C4H9O)[61]和烷基酯体系[36, 39]的研究中得到验证.采用QCISD(T)/CBS方法得到的计算结果可作为基准值用来验证ONIOM方法的计算准确性.
1.3 ONIOM高精度单点能计算目前,ONIOM[QCISD(T)/CBS : DFT]方法将反应体系的结构定义为两层, 针对不同的层级采用不同的理论方法进行处理. 图 2为采用ONIOM方法处理正壬烷(n-C9H20)与氢自由基的氢提取反应(C9H20+H→CH3(CH2)3CH ·(CH2)3CH3+H2)的示意图.对于高精度层, 即分子中的化学活性部分(chemically active portion,CAP)(图 2中红色圈中部分), 采用高精度QCISD(T)/CBS方法处理, 其余部分采用较低精度的DFT方法配合6-311++G(d, p)基组来处理.此外, ONIOM方法引进氢原子用于弥补由分层计算导致层级间出现的不饱和键.为了定量确定CAP的大小对于计算精度的影响, 引进了两个变量(N1, N2)对此问题进行了分析, 其中N表示主链上除发生氢提取反应的反应位置上的亚甲基(CH2)或甲基(CH3)外, 每一边碳或氧等重原子的数量, 而且如有连接这些原子的官能团也须纳入CAP之中.并且需要注意的是, 在处理过渡态时, CAP须考虑进攻的氢自由基以及邻近的亚甲基和甲基, 如有官能团(碳碳双键, 羟基和羰基等)也须考虑在内.以图 2中过渡态为例, 图中标识1, 9分别表示n-C9H20两端的甲基, 标识2到8分别表示n-C9H20中的亚甲基, 由图中圈内部分可知, 标识5为发生氢提取反应的位置, CAP(2, 2)包括了进攻的氢自由基和标识5以及其两侧的标识3,4和标识6,7的亚甲基, 它们共同组成了化学活性部分.关于选择合适CAP(N1, N2)尺寸的问题将在后文详细叙述.
![]() |
图 2 ONIOM/CAP(2, 2)方法关于C9H20+H→CH3(CH2)3CH ·(CH2)3CH3+H2氢提取反应的原理示意图 Fig.2 Illustration of ONIOM/CAP(2, 2) method for the hydrogen abstraction reaction by H-radical, C9H20+H→CH3(CH2)3CH ·(CH2)3CH3+H2 |
ONIOM方法近似求解系统的能量, 通过对较低精度系统下计算得到的系统能量, 采用CAP部分高精度能量与低精度能量的差值进行校正, 具体表示如下
EONIOM[high : low]=Elow(R)+Ehigh(CAP)-Elow(CAP)
其中, 低层低精度部分使用DFT方法处理, 高层高精度部分采用以DZ/TZ为基组的QCISD(T)/CBS方法处理.
ONIOM[QCISD(T)/CBS : DFT]能量由以下公式计算得到
E[QCISD(T)/CBS : DFT]=
EONIOM[QCISD(T)/CBS : DFT]DZ→TZ+
{EONIOM[MP2/CBS : DFT]TZ→QZ-
EONIOM[MP2/CBS : DFT]DZ→TZ}
EONIOM[QCISD(T)/CBS : DFT]DZ→TZ=
EONIOM[QCISD(T)/TZ : DFT]+
{EONIOM[QCISD(T)/QZ : DFT]-
EONIOM[QCISD(T)/DZ : DFT]}×0.462 9
EONIOM[MP2/CBS : DFT]TZ→QZ=
EONIOM[MP2/QZ : DFT]+
{EONIOM[MP2/QZ : DFT]-
EONIOM[MP2/TZ : DFT]}×0.693 8
EONIOM[MP2/CBS : DFT]DZ→TZ=
EONIOM[MP2/TZ : DFT]+
{EONIOM[MP2/TZ : DFT]-
EONIOM[MP2/DZ : DFT]}×0.462 9
能垒的零点能校正通过计算过渡态与反应物EONIOM[QCISD(T)/CBS : DFT]的单点能加上零点能的能量差得到, 反应热的零点能校正通过计算反应物和产物EONIOM[QCISD(T)/CBS : DFT]的单点能加上零点能的能量差得到, 具体表示如下
EB=EONIOM[QCISD(T)/CBS : DFT+ZPE]TS-
EONIOM[QCISD(T)/CBS : DFT]+ZPE]reactants
HR=EONIOM[QCISD(T)/CBS : DFT+ZPE]products-
EONIOM[QCISD(T)/CBS : DFT]+ZPE]reactants
所有相关计算可使用Gaussian 09程序包[55]进行.
2 ONIOM方法的验证与应用 2.1 两种QCISD(T)/CBS方法的对比与验证关于验证ONIOM方法的正确性与准确性, 首先需要得到可靠的基准值进行对比. Zhang等[39]使用这两种被广泛接受的高精度方法[QCISD(T)/CBS]1和[QCISD(T)/CBS]2对烃基脂分子的氢提取反应(CnH2n+1COOCmH2m+1+H, n=0~2, m=1~2, 所有反应详见表 1)进行了量化计算, 采用这两种方法计算得到的能量(能垒和反应热)差异基本都小于0.1 kcal/mol, 如图 3所示(图中R表示不同的反应, 详见表 1).由此, 我们可以判定这两种高精度方法的计算精度基本相同.同时, Zhang等[36]在关于丁酸甲酯(C5H10O2)氢提取反应从头算化学动力学的研究中, 通过与Liu等[35]采用CCSD(T)/CBS//B3LYP/6-311++G(d, p)的计算结果对比, 也证明了[QCISD(T)/CBS]2高精度方法的可靠性.由于计算成本的原因, [QCISD(T)/CBS]1高精度方法不适用于非氢原子数大于6的反应体系.因此, 可使用适用于相对更大体系的[QCISD(T)/CBS]2高精度方法, 建立更为全面的基准数据, 用于验证ONIOM[QCISD(T)/CBS : DFT]方法的可靠性.
![]() |
图 3 QCISD(T)/CBS1和QCISD(T)/CBS2方法计算CnH2n+1 COOCmH2m+1+H(n=0~2, m=1~2)的能量(能垒和反应热)的误差 Fig.3 Difference (in absolute value) of calculated results for the hydrogen abstraction reaction CnH2n+1COOCmH2m+1+H using QCISD(T)/CBS1 and QCISD(T)/CBS2 methods |
![]() |
下载CSV 表 1 用于对比[QCISD(T)/CBS]1和[QCISD(T)/CBS]2方法计算结果的反应体系 Tab.1 Reactions for comparison of calculated results using [QCISD(T)/CBS]1 and [QCISD(T)/CBS]2 methods |
对于ONIOM方法, 一个至关重要的问题就是在保证计算精度同时兼顾计算成本的前提下如何得到最优的CAP.
Zhang等[39]通过系统地测试己酸甲酯和正壬烷与氢自由基的氢提取反应中(C5H11COOCH3+H和n-C9H20+H)所有可能的(N1, N2)的组合, 总结了关于最优CAP选择的规律.这是目前采用高精度的[QCISD(T)/CBS]2方法可研究的最大体系, 并且考虑了C5H11COOCH3+H氢提取反应中全部可能的反应位置, 和n-C9H20+H氢提取反应中以正壬烷中心亚甲基作为进攻位的氢提取位置.
如图 4(a)所示, 分别对应以下5个反应:C9H20+H→CH3(CH2)3CH ·(CH2)3CH3+H2; C7H14O2+H→CH3(CH2)3CH ·COOCH3+H2(1), CH3(CH2)2CH ·CH2COOCH3+H2 (2), CH3CH2CH ·(CH2)2COOCH3+H2(3), CH3CH ·(CH2)3COOCH3+H2(4).采用ONIOM[QCISD(T)/CBS : DFT](图中标识位EB[ONIOM])方法和[QCISD(T)/CBS]2(图中标识为EB [QCISD(T)/CBS])方法计算能垒的差异是基于所选取CAP尺寸的大小, 即N1(N1=1~4)和N2(N2=1~4)的函数.在所有测试的CAP中, EB[ONIOM]与EB[QCISD(T)/CBS]的误差均小于0.8 kcal/mol, 值得注意的是, 相对较大的误差仅出现在CAP(1, 1), 即N1=1和N2=1.而在选择CAP(2, 2)时, 即N1和N2都为2, 两种方法对于能垒的差异基本都小于0.1 kcal/mol.所以对于类似的反应热, 两种方法计算得到的能量差异, 如图 4(b)所示.当采用CAP(2, 2)或更大的CAP尺寸时, HR[ONIOM]与HR[QCISD(T)/CBS]两种方法的计算误差基本小于0.1 kcal/mol.而相对较大的误差也与关于能垒的计算结果类似, 也发生在CAP(1, 1), 即N1=1和N2=1.通过对ONIOM方法采用不同CAP(N1, N2)尺寸与QCISD(T)/CBS方法的比较结果发现, CAP(2, 2)是研究反应时要求的最小需要尺寸, 而且CAP(2, 2)对于其他相类似的体系已经足够精确, 计算误差基本小于0.1 kcal/mol.
![]() |
图 4 5个目标反应能量(能垒和反应热)误差的变化.每个二维平面表示一个固定的N2平面(从左到右依次为N2=1~4)[39] Fig.4 Variation of the difference of the calculated energy barriers (EB and HR) for the five reactions. Each two-dimensional plane represents a constant N2 plane (N2=1~4 from left to right)[39] |
需要注意的是, 在处理含有官能团的反应体系, 尤其是包含酯类官能团(-COO-)时, 需要把官能团看作为一个整体来处理.以己酸甲酯与氢自由基的氢提取反应为例(C5H11COOCH3+H→CH3(CH2)2CH ·CH2COOCH3+H2), 见图 5, 图中标识Me1和Me2表示烷基链和酯基上的甲基, 标识1~4分别表示位于由近到远距离酯基位置的亚甲基.由图 5中圈内部分可知, 发生在标识2位置的氢提取反应的CAP(2, 2)需要把酯类官能团(-COO-)作为一个整体, 全部纳入CAP之中.所以对于此反应的CAP(2, 2)包括4个亚甲基(标识1~4)和酯类官能团(-COO-).
![]() |
图 5 ONIOM/CAP(2, 2)方法关于C5H11COOCH3+H→CH3(CH2)2CH ·CH2COOCH3+H2氢提取反应的原理示意图 Fig.5 Illustration of ONIOM/CAP(2, 2) method for the hydrogen abstraction reaction by H-radical, C5H11COOCH3+H→CH3(CH2)2CH ·CH2COOCH3+H2 |
在之前关于ONIOM方法介绍中, 研究对象是饱和烷基酯体系(己酸甲酯)和饱和烃(正壬烷)(C5H11COOCH3+H和C9H20+H的氢提取反应).而对于ONIOM方法在不饱和体系中的验证与应用, Zhang等也做了详细全面的研究[53], 相比饱和烷基酯体系的氢提取反应, 不饱和甲酯体系(CnH2n-1COOCH3)由于其分子结构相对复杂, 碳碳双键(C=C), 为了保持与之前ONIOM/CAP(2, 2)研究的一致性, 也为了考虑由于碳碳双键带来的共轭效应, 引进了一个新的描述, 即化学活性中心(chemically active center, CAC).
如图 6所示, 采用ONIOM/CAP(2, 2)方法处理油酸甲酯与氢原子的氢提取反应, 标识1~16分别表示位于由近到远距离酯基位置的亚甲基(或CH基团).标识Me1和Me2分别表示烷基链和酯基上的甲基.碳碳双键位于标识8与9, 两个CH基团之间.因为标识7和10位于相对于碳碳双键的β位, 所以位于σ轨道的电子将会与位于临近π轨道的电子发生相互作用, 从而产生一个扩展的分子轨道, 原因是碳碳双键产生的强共轭效应导致电子离域化增强.而且, 当氢提取反应发生在标识7或10的位置时, 受到强共轭效应的烯丙基也随之产生.所以, CAC也必须将碳碳双键纳入其中.对于图 6中, 发生在以标识7(亚甲基)作为氢提取位置的氢提取反应(CH3(CH2)7CH=CH(CH2)7COOCH3+H→CH3(CH2)7CH=CHCH ·(CH2)6COOCH3+H2), 亚甲基(标识7)和临近的碳碳双键(标识8和9)都需要纳入CAC.因此, CAC和其两侧的亚甲基(标识5, 6, 10和11)共同组成了CAP(2, 2).
![]() |
图 6 ONIOM/CAP(2, 2)方法关于CH3(CH2)7CH=CH(CH2)7 COOCH3+H→CH3(CH2)7CH=CHCH ·(CH2)6COOCH3+H2氢提取反应的原理示意图[53] Fig.6 Illustration of ONIOM/CAP(2, 2) method for the hydrogen abstraction reaction by H-radical, CH3(CH2)7CH=CH(CH2)7COOCH3+H→CH3(CH2)7CH=CHCH ·(CH2)6COOCH3+H2[53] |
为了进一步验证ONIOM/CAP(2, 2)方法对于处理不饱和体系中氢提取反应的可靠性, Zhang等[53]分别采用了ONIOM/CAP(2, 2)和[QCISD(T)CBS]2这两种方法计算并比较了CnH2n-1COOCH3+H(n=2~5)氢提取反应的能垒和反应热.通过比较发现, 使用这两种方法得到的能量结果, 误差基本都小于0.15 kcal/mol.
2.4 ONIOM方法处理氢加成反应在采用ONIOM/CAP(2, 2)方法处理氢加成反应中, CAC的选取方法与之前处理饱和烷基酯体系氢提取反应相同.
如图 7所示, 以油酸甲酯(CH3(CH2)7CH=CH(CH2)7COOCH3)的氢加成反应为例, 其中标识1~16分别表示位于由近到远距离酯基位置的亚甲基(或CH基团), 标识Me1和Me2分别表示烷基链和酯基上的甲基.碳碳双键位于标识8与9, 两个CH基团之间.对于图 7中, 发生在标识9(CH基团)的氢加成反应(CH3(CH2)7CH=CH(CH2)7COOCH3+H→CH3(CH2)8CH ·(CH2)7COOCH3), CAC(标识9, CH基团)和其两侧的亚甲基或CH基团(标识8, 7, 10和11)共同组成了CAP(2, 2).
![]() |
图 7 ONIOM/CAP(2, 2)方法关于CH3(CH2)7CH=CH(CH2)7 COOCH3+H→CH3(CH2)8CH ·(CH2)7COOCH3氢加成反应的原理示意图[53] Fig.7 Illustration of ONIOM/CAP(2, 2) method for the hydro- gen addition reaction of CH3(CH2)7CH=CH(CH2)7COOCH3+ H→CH3(CH2)8CH ·(CH2)7COOCH3[53] |
为了进一步验证ONIOM/CAP(2, 2)方法对于处理不饱和体系中氢提取反应的可靠性, Zhang等[53]采用了ONIOM/CAP(2, 2)和[QCISD(T)CBS]2这两种方法计算并比较了CnH2n-1COOCH3+H(n=2~5)氢加成反应的能垒和反应热的结果.通过比较发现, 使用这两种方法得到的能量结果, 误差基本都小于0.15 kcal/mol.
2.5 ONIOM方法处理异构化反应对于不饱和甲酯与氢自由基反应的体系中(CnH2n-1COOCH3+H)出现的异构化反应, 其反应特征是分子间氢的转移, 即从某一甲基或亚甲基转移到另一个CH基团.需要注意的是, 在采用ONIOM/CAP(2, 2)方法处理此类反应时, 所有位于氢原子转移之间的基团都需要当作化学活化部分考虑.所以, 对于此类反应, 存在两个CACs, 一个是失去氢原子的甲基或亚甲基基团, 另一个是得到氢原子的CH基团.因此, 由于两个CACs的存在, CAP(2, 2)也随之增大, 即在每个CAC两侧与之相邻的两个亚甲基(或甲基, CH基团)都需要纳入CAP之中.以油酸甲酯自由基(CH3(CH2)8CH ·(CH2)7COOCH3)的异构化反应为例, 如图 8所示, 其中标识1~16分别表示位于由近到远距离酯基的位置的亚甲基(或CH基团), 标识Me1和Me2分别表示烷基链和酯基上的甲基.
![]() |
图 8 ONIOM/CAP(2, 2)方法关于CH3(CH2)8CH · (CH2)7 COOCH3→CH3(CH2)13CH ·(CH2)2COOCH3异构化反应的原理示意图[53] Fig.8 Illustration of ONIOM/CAP(2, 2) method for the isomerization reaction of CH3(CH2)8CH ·(CH2)7COOCH3→ CH3(CH2)13CH ·(CH2)2COOCH3[53] |
对于图 8中所发生的异构化反应(CH3(CH2)8CH ·(CH2)7COOCH3→CH3(CH2)13CH ·(CH2)2COOCH3), 氢原子从标识3(亚甲基)转移到标识8(CH基团), 所以标识3与标识8都被认为是CAC, 因此CAP(2, 2)包括了第1个CAC(标识3)两侧的4个亚甲基(标识1, 2, 4和5)和第2个CAC(标识8)两侧的4个亚甲基(标识6, 7, 9和10), 因此CAP(2, 2)共包括了10个非氢原子, 这也是目前ONIOM方法可以处理的最大体系.
为了进一步验证ONIOM/CAP(2, 2)方法对于处理异构化反应的可靠性, Zhang等[53]采用ONIOM/CAP(2, 2)和[QCISD(T)CBS]2这两种方法计算并比较了CnH2nCOOCH3(n=2~5)异构化反应的能垒和反应热的结果.通过比较发现, 使用这两种方法得到的能量结果, 误差基本都小于0.10 kcal/mol.
2.6 ONIOM方法处理β位裂解反应在采用ONIOM/CAP(2, 2)方法处理β位裂解反应中, CAC需要包括位于自由基位的基团和与之相邻的α,β基团.以油酸甲酯自由基(CH3(CH2)15CH ·COOCH3)的β位裂解反应为例, 如图 9所示, 其中标识1~16分别表示位于由近到远距离酯基的位置的亚甲基(或CH基团), 标识Me1和Me2分别表示烷基链和酯基上的甲基.对于图 9中所发生的β位裂解反应(CH3(CH2)15CH ·COOCH3→CH3(CH2)13CH2 ·+H2C=CHCOOCH3), CH基团(标识1)和临近的亚甲基(标识2和3)都需要纳入CAC.因此, CAC和其两侧的亚甲基和酯类官能团(标识4, 5和-COO-)共同组成了CAP(2, 2).
![]() |
图 9 ONIOM/CAP(2, 2)方法关于CH3(CH2)15CH · COOCH3→CH3(CH2)13CH2 ·+H2C=CHCOOCH3, β位裂解反应的原理示意图[53] Fig.9 Illustration of ONIOM/CAP(2, 2) method for the β-scission reaction of CH3(CH2)15CH ·COOCH3→ CH3(CH2)13CH2 ·+H2C=CHCOOCH3[53] |
为了进一步验证ONIOM/CAP(2, 2)方法对于处理异构化反应的可靠性, Zhang等[53]采用了ONIOM/CAP(2, 2)和[QCISD(T)CBS]2这两种方法计算并比较了CnH2nCOOCH3(n=2~5)异构化反应的能垒和反应热的结果.通过比较发现, 使用这两种方法得到的能量结果, 误差基本都小于0.15 kcal/mol.
3 结论对于大分子燃料, 通常采用类比的方法来估测其反应速率常数以构建其燃烧反应模型, 但是无论采用何种方法都会为模型带来较大的计算误差.为了更好地将大分子燃料应用在发动机中, 对其进行高精度的化学反应动力学研究是十分重要的.因此, 本文较为详细地分析了现今在发展航空煤油和生物柴油燃烧化学反应动力学研究中遇到的主要困难, 即获得大分子的高精度反应能量.为此, 本文主要关注了目前可适用于大分子体系的高精度量化计算方法, 并对一种两层的ONIOM方法, ONIOM[QCISD(T)/CBS : DFT], 做了详细的阐述和使用说明的介绍.这些方法的提出, 不仅为研究大分子燃料体系的能量计算提供了准确与可行的计算方法, 而且为化学动力学计算提供了可靠的热化学数据, 对大分子燃料燃烧反应的研究具有重要意义.
致谢 本文的研究得到了国家自然科学基金“面向发动机的湍流燃烧基础研究”重大研究计划(91641105)和香港理工大学的研究基金(4-BCE8, G-YBXN)的资助.[1] |
Chung W M, Wang Q, Sezerman U, et al. Analysis of aviation turbine fuel composition by laser Raman spectros-copy[J]. Applied Spectroscopy, 1991, 45(9): 1527-1532. DOI:10.1366/0003702914335643 |
[2] |
Edwards T, Maurice L Q. Surrogate mixtures to represent complex aviation and rocket fuels[J]. Journal of Propulsion and Power, 2001, 17(2): 461-466. DOI:10.2514/2.5765 |
[3] |
Dahm K D, Virk P S, Bounaceur R, et al. Experimental and modelling investigation of the thermal decomposition of n-dodecane[J]. Journal of Analytical and Applied Pyrolysis, 2004, 71(2): 865-881. DOI:10.1016/j.jaap.2003.11.005 |
[4] |
Dagaut P, Cathonnet M. The ignition, oxidation, and combustion of kerosene:A review of experimental and kinetic modeling[J]. Progress in Energy and Combustion Science, 2006, 32(1): 48-92. DOI:10.1016/j.pecs.2005.10.003 |
[5] |
Dagaut P, El Bakali A, Ristori A. The combustion of kerosene:experimental results and kinetic modelling using 1-to 3-component surrogate model fuels[J]. Fuel, 2006, 85(7/8): 944-956. |
[6] |
Honnet S, Seshadri K, Niemann U, et al. A surrogate fuel for kerosene[J]. Proceedings of the Combustion Institute, 2009, 32(1): 485-492. |
[7] |
Marsh G. Biofuels:aviation alternative?[J]. Renewable Energy Focus, 2008, 9(4): 48-51. DOI:10.1016/S1471-0846(08)70138-0 |
[8] |
Basha S A, Gopal K R, Jebaraj S. A review on biodiesel production, combustion, emissions and performance[J]. Renewable and Sustainable Energy Reviews, 2009, 13(6/7): 1628-1634. |
[9] |
Giakoumis E G, Rakopoulos C D, Dimaratos A M, et al. Exhaust emissions of diesel engines operating under transient conditions with biodiesel fuel blends[J]. Progress in Energy and Combustion Science, 2012, 38(5): 691-715. DOI:10.1016/j.pecs.2012.05.002 |
[10] |
Lapuerta M, Armas O, Rodriguez-Fernandez J. Effect of biodiesel fuels on diesel engine emissions[J]. Progress in Energy and Combustion Science, 2008, 34(2): 198-223. DOI:10.1016/j.pecs.2007.07.001 |
[11] |
Sun J F, Caton J A, Jacobs T J. Oxides of nitrogen emissions from biodiesel-fuelled diesel engines[J]. Progress in Energy and Combustion Science, 2010, 36(6): 677-695. DOI:10.1016/j.pecs.2010.02.004 |
[12] |
Agarwal A K. Biofuels (alcohols and biodiesel) applications as fuels for internal combustion engines[J]. Pro-gress in Energy and Combustion Science, 2007, 33(3): 233-271. DOI:10.1016/j.pecs.2006.08.003 |
[13] |
Atadashi I M, Aroua M K, Aziz A A. High quality biodiesel and its diesel engine application:a review[J]. Renewable and Sustainable Energy Reviews, 2010, 14(7): 1999-2008. DOI:10.1016/j.rser.2010.03.020 |
[14] |
Knothe G. Biodiesel and renewable diesel:a compari-son[J]. Progress in Energy and Combustion Science, 2010, 36(3): 364-373. DOI:10.1016/j.pecs.2009.11.004 |
[15] |
Jenkins R W, Munro M, Nash S, et al. Potential renewable oxygenated biofuels for the aviation and road transport sectors[J]. Fuel, 2013, 103: 593-599. DOI:10.1016/j.fuel.2012.08.019 |
[16] |
Chuck C J, Donnelly J. The compatibility of potential bioderived fuels with Jet A-1 aviation kerosene[J]. Applied Energy, 2014, 118: 83-91. DOI:10.1016/j.apenergy.2013.12.019 |
[17] |
Starik A K, Titova N S, Torokhov S A. Kinetics of oxidation and combustion of complex hydrocarbon fuels:Aviation kerosene[J]. Combustion, Explosion, and Shock Waves, 2013, 49(4): 392-408. DOI:10.1134/S0010508213040023 |
[18] |
Coniglio L, Bennadji H, Glaude P A, et al. Combustion chemical kinetics of biodiesel and related compounds (methyl and ethyl esters):experiments and modeling-advances and future refinements[J]. Progress in Energy and Combustion Science, 2013, 39(4): 340-382. DOI:10.1016/j.pecs.2013.03.002 |
[19] |
Davis A C, Francisco J S. Ab initio study of key branching reactions in biodiesel and Fischer-Tropsch fuels[J]. Journal of the American Chemical Society, 2011, 133(47): 19110-19124. DOI:10.1021/ja205516s |
[20] |
Kohse-Höinghaus K, Oβwald P, Cool T A, et al. Biofuel combustion chemistry:from ethanol to biodiesel[J]. Angewandte Chemie International Edition, 2010, 49(21): 3572-3597. DOI:10.1002/anie.200905335 |
[21] |
Westbrook C K, Pitz W J, Herbinet O, et al. A comprehensive detailed chemical kinetic reaction mechanism for combustion of n-alkane hydrocarbons from n-octane to n-hexadecane[J]. Combustion and Flame, 2009, 156(1): 181-199. |
[22] |
Strelkova M I, Kirillov I A, Potapkin B V, et al. Detailed and reduced mechanisms of Jet A combustion at high temperatures[J]. Combustion Science and Technology, 2008, 180(10/11): 1788-1802. |
[23] |
Fisher E M, Pitz W J, Curran H J, et al. Detailed chemical kinetic mechanisms for combustion of oxygenated fuels[J]. Proceedings of the Combustion Institute, 2000, 28(2): 1579-1586. |
[24] |
Lai J Y W, Lin K C, Violi A. Biodiesel combustion:advances in chemical kinetic modeling[J]. Progress in Energy and Combustion Science, 2011, 37(1): 1-14. DOI:10.1016/j.pecs.2010.03.001 |
[25] |
Westbrook C K, Naik C V, Herbinet O, et al. Detailed chemical kinetic reaction mechanisms for soy and rapeseed biodiesel fuels[J]. Combustion and Flame, 2011, 158(4): 742-755. DOI:10.1016/j.combustflame.2010.10.020 |
[26] |
Hoekman S K, Broch A, Robbins C, et al. Review of biodiesel composition, properties, and specifications[J]. Renewable and Sustainable Energy Reviews, 2012, 16(1): 143-169. |
[27] |
Ma F R, Hanna M A. Biodiesel production:a review[J]. Bioresource Technology, 1999, 70(1): 1-15. DOI:10.1016/S0960-8524(99)00025-5 |
[28] |
Dooley S, Curran H J, Simmie J M. Autoignition measurements and a validated kinetic model for the biodiesel surrogate, methyl butanoate[J]. Combustion and Flame, 2008, 153(1/2): 2-32. |
[29] |
Gaïl S, Thomson M J, Sarathy S M, et al. A wide-ranging kinetic modeling study of methyl butanoate com-bustion[J]. Proceedings of the Combustion Institute, 2007, 31(1): 305-311. DOI:10.1016/j.proci.2006.08.051 |
[30] |
Hayes C J, Burgess Jr D R. Exploring the oxidative decompositions of methyl esters:Methyl butanoate and methyl pentanoate as model compounds for biodiesel[J]. Proceedings of the Combustion Institute, 2009, 32(1): 263-270. |
[31] |
Metcalfe W K, Dooley S, Curran H J, et al. Experimen-tal and modeling study of C5H10O2 ethyl and methyl esters±[J]. The Journal of Physical Chemistry A, 2007, 111(19): 4001-4014. DOI:10.1021/jp067582c |
[32] |
Yang B, Westbrook C K, Cool T A, et al. Fuel-specific influences on the composition of reaction intermediates in premixed flames of three C5H10O2 ester isomers[J]. Physical Chemistry Chemical Physics, 2011, 13(15): 6901-6913. DOI:10.1039/c0cp02065f |
[33] |
Diévart P, Won S H, Dooley S, et al. A kinetic model for methyl decanoate combustion[J]. Combustion and Flame, 2012, 159(5): 1793-1805. DOI:10.1016/j.combustflame.2012.01.002 |
[34] |
Herbinet O, Pitz W J, Westbrook C K. Detailed chemical kinetic oxidation mechanism for a biodiesel surrogate[J]. Combustion and Flame, 2008, 154(3): 507-528. DOI:10.1016/j.combustflame.2008.03.003 |
[35] |
Liu W, Sivaramakrishnan R, Davis M J, et al. Development of a reduced biodiesel surrogate model for compression ignition engine modeling[J]. Proceedings of the Combustion Institute, 2013, 34(1): 401-409. |
[36] |
Zhang L D, Chen Q X, Zhang P. A theoretical kinetics study of the reactions of methylbutanoate with hydrogen and hydroxyl radicals[J]. Proceedings of the Combustion Institute, 2015, 35(1): 481-489. |
[37] |
Truhlar D G, Garrett B C, Klippenstein S J. Current status of transition-state theory[J]. The Journal of Physical Chemistry, 1996, 100(31): 12771-12800. DOI:10.1021/jp953748q |
[38] |
Fernández-Ramos A, Miller J A, Klippenstein S J, et al. Modeling the kinetics of bimolecular reactions[J]. Chemical Reviews, 2006, 106(11): 4518-4584. DOI:10.1021/cr050205w |
[39] |
Zhang L D, Zhang P. Towards high-level theoretical studies of large biodiesel molecules:an ONIOM[J]. Physical Chemistry Chemical Physics, 2015, 17(1): 200-208. |
[40] |
Klippenstein S J, Pande V S, Truhlar D G. Chemical kinetics and mechanisms of complex systems:a perspective on recent theoretical advances[J]. Journal of the Ameri-can Chemical Society, 2014, 136(2): 528-546. DOI:10.1021/ja408723a |
[41] |
Yu T, Zheng J J, Truhlar D G. Multi-structural variational transition state theory. Kinetics of the 1, 4-hydrogen shift isomerization of the pentyl radical with torsional anharmonicity[J]. Chemical Science, 2011, 2(11): 2199-2213. DOI:10.1039/c1sc00225b |
[42] |
Zheng J J, Yu T, Papajak E, et al. Practical methods for including torsional anharmonicity in thermochemical calculations on complex molecules:The internal-coordinate multi-structural approximation[J]. Physical Chemistry Chemical Physics, 2011, 13(23): 10885-10907. DOI:10.1039/c0cp02644a |
[43] |
Zheng J J, Seal P, Truhlar D G. Role of conformational structures and torsional anharmonicity in controlling chemical reaction rates and relative yields:butanal+HO2 reactions[J]. Chemical Science, 2013, 4(1): 200-212. |
[44] |
Papajak E, Truhlar D G. What are the most efficient basis set strategies for correlated wave function calcula-tions of reaction energies and barrier heights?[J]. The Journal of Chemical Physics, 2012, 137(6): 064110. DOI:10.1063/1.4738980 |
[45] |
Zádor J, Taatjes C A, Fernandes R X. Kinetics of elementary reactions in low-temperature autoignition che-mistry[J]. Progress in Energy and Combustion Science, 2011, 37(4): 371-421. DOI:10.1016/j.pecs.2010.06.006 |
[46] |
Goldsmith C F, Magoon G R, Green W H. Database of small molecule thermochemistry for combustion[J]. The Journal of Physical Chemistry A, 2012, 116(36): 9033-9057. DOI:10.1021/jp303819e |
[47] |
Li W, Li S H, Jiang Y S. Generalized energy-based fragmentation approach for computing the ground-state ener-gies and properties of large molecules[J]. The Journal of Physical Chemistry A, 2007, 111(11): 2193-2199. DOI:10.1021/jp067721q |
[48] |
Wang K D, Li W, Li S H. Generalized energy-based fragmentation CCSD (T)-F12a method and application to the relative energies of water clusters (H2O)20[J]. Journal of Chemical Theory and Computation, 2014, 10(4): 1546-1553. DOI:10.1021/ct401060m |
[49] |
Wu J J, Ning H B, Ma L H, et al. Accurate prediction of bond dissociation energies of large n-alkanes using ONIOM-CCSD (T)/CBS methods[J]. Chemical Physics Letters, 2018, 699: 139-145. DOI:10.1016/j.cplett.2018.03.041 |
[50] |
Wu J J, Ning H B, Ma L H, et al. Cascaded group-additivity ONIOM:A new method to approach CCSD (T)/CBS energies of large aliphatic hydrocarbons[J]. Combustion and Flame, 2019, 201: 31-43. DOI:10.1016/j.combustflame.2018.12.012 |
[51] |
Chi Y W, You X Q, Zhang L D, et al. Utilization of generalized energy-based fragmentation method on the study of hydrogen abstraction reactions of large methyl esters[J]. Combustion and Flame, 2018, 190: 467-476. DOI:10.1016/j.combustflame.2017.12.021 |
[52] |
Chi Y W, You X Q. Kinetics of hydrogen abstraction reactions of methyl palmitate and octadecane by hydrogen atoms[J]. The Journal of Physical Chemistry A, 2019, 123(14): 3058-3067. DOI:10.1021/acs.jpca.8b08802 |
[53] |
Zhang L D, Meng Q H, Chi Y C, et al. Toward high-level theoretical studies of large biodiesel molecules:an ONIOM[J]. The Journal of Physical Chemistry A, 2018, 122(21): 4882-4893. DOI:10.1021/acs.jpca.8b02327 |
[54] |
Meng Q, Chi Y, Zhang L, et al. Towards high-level theoretical studies of large biodiesel molecules:an ONIOM/RRKM/Master-equation approach to the isomerization and dissociation kinetics of methyl decanoate radicals[J]. Physical Chemistry Chemical Physics, 2019, 21(9): 5232-5242. DOI:10.1039/C8CP05593A |
[55] |
Frisch M J, Trucks G W, Schlegel H B, et al. Gaussian 09[M]. Wallingford, CT: Gaussian, Inc., 2009.
|
[56] |
Pitz W J, Mueller C J. Recent progress in the develo-pment of diesel surrogate fuels[J]. Progress in Energy and Combustion Science, 2011, 37(3): 330-350. DOI:10.1016/j.pecs.2010.06.004 |
[57] |
Krishnan R, Binkley J S, Seeger R, et al. Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions[J]. The Journal of Chemical Physics, 1980, 72(1): 650-654. DOI:10.1063/1.438955 |
[58] |
Becke A D. Density-functional thermochemistry. III. The role of exact exchange[J]. The Journal of Chemical Physics, 1993, 98(7): 5648-5652. DOI:10.1063/1.464913 |
[59] |
Dunning Jr T H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen[J]. The Journal of Chemical Physics, 1989, 90(2): 1007-1023. DOI:10.1063/1.456153 |
[60] |
Martin J M, Uzan O. Basis set convergence in second-row compounds. The importance of core polarization functions[J]. Chemical Physics Letters, 1988, 282(1): 16-24. |
[61] |
Zhang P, Klippenstein S J, Law C K. Ab initio kinetics for the decomposition of hydroxybutyl and butoxy radicals of n-butanol[J]. The Journal of Physical Chemistry A, 2013, 117(9): 1890-1906. DOI:10.1021/jp400155z |