0 引言
以醚类为代表的含氧燃料可由生物质转化而来,具有可再生的优势,其分子中含有氧原子,相比烷烃含碳更少,更容易实现完全燃烧;产生的CO2更少,是一类清洁燃料[3]。
Fig. 1 Main paths of low temperature oxidation for DEE图1 乙醚低温氧化的主要路径 |
丙醚有正丙醚和异丙醚两种同分异构体,本文针对两者反应生成的代表性OOQOOH自由基,开展其最重要的反应路径——氢转移裂解生成氢过氧化物醛酮类物质（ketohydroperoxide, KHP）和OH的基元反应特征研究,获得相应反应的速率常数。
1 计算方法
采用量子化学方法对基元反应的反应物、过渡态及产物的分子结构进行优化,并计算其单点能。然后基于过渡态理论计算其速率常数。
1.1 量子化学方法
使用一种密度泛函[M062-X/6-311++g(d, p)]方法,对反应物、过渡态和产物的能量最低的构象进行结构优化,并利用CCSD(T)-F12a/jun-cc-pVTZ[21,22]方法计算电子能。同时,针对每个物种,使用5种密度泛函方法：M05-2X[23]、TPSSH、BHHLYP、WB97X-D[24]、B3LYP;并搭配选用6个不同基组：MG3S、def2-TZVPP、cc-pVTZ、aug-cc-pVTZ、jul-cc-pVTZ、jun-cc-pVTZ[25],即使用30种不同的计算方法,优化其分子结构并计算单点能。比较这些方法获得的单点能与基准方法获得的单点能差异的平均值,选择偏差最小的方法作为全局计算方法。各方法对比结果如表1所示,通过比较特定反应36个数据的平均误差,发现BHHLYP/cc-pVDZ与CCSD(T)-F12a/jun-cc-pVTZ的差值仅为1.55 kJ/mol,要明显低于其他计算方法的结果,因此选用BHHLYP/cc-pVDZ和基组作为最优的结构优化和频率计算方法。获得的单点能作为基准,用于筛选能够较好计算物种能量的密度泛函方法。
Table 1 Average error relative to CCSD(T)-F12a/jun-cc-pVTZ表1 相对于CCSD(T)-F12a/jun-cc-pVTZ的平均误差 |
| 方法/基组 | 相对能量/(kJ/mol) |
|---|---|
| CCSD(T)-F12a/jun-cc-pVTZ | 0.00 |
| BHHLYP-cc-pVDZ | 1.55 |
| WB97X-D-cc-pVDZ | 2.22 |
| BHHLYP-MG3S | 2.30 |
| BHHLYP-jun-cc pVTZ | 2.39 |
| BHHLYP-jul-cc-pVTZ | 2.39 |
| BHHLYP-aug-cc-pVTZ | 2.39 |
| BHHLYP-def2-TZVPP | 2.43 |
| M05-2x-cc-pVDZ | 2.64 |
| WB97X-D-MG3S | 2.97 |
| b3lyp-cc-pVDZ | 2.97 |
| WB97X-D-jun-cc-pVTZ | 3.06 |
| WB97X-D-jul-cc-pVTZ | 3.06 |
| WB97X-D-aug-cc-pVTZ | 3.06 |
| WB97X-D-def2-TZVPP | 3.14 |
| TPSSH-cc-pVDZ | 3.14 |
| M05-2x-MG3S | 3.39 |
| M05-2x-aug-cc-pVTZ | 3.52 |
| M05-2x-jun-cc-pVTZ | 3.52 |
| M05-2x-jul-cc-pVTZ | 3.52 |
| M05-2x-def2-TZVPP | 3.56 |
| b3lyp-MG3S | 3.81 |
| TPSSH-MG3S | 3.85 |
| b3lyp-jun-cc-pVTZ | 3.89 |
| b3lyp-jul-cc-pVTZ | 3.89 |
| b3lyp-aug-cc-pVTZ | 3.89 |
| TPSSH-jun-cc-pVTZ | 3.93 |
| TPSSH-jul-cc-pVTZ | 3.93 |
| TPSSH-aug-cc-pVTZ | 3.93 |
| b3lyp-def2-TZVPP | 3.93 |
分子空间构象的计算方法是：独立旋转每个除甲基外的C—C键二面角获得大量的初猜分子结构,然后利用选定的密度泛函方法,即BHHLYP/cc-pVDZ方法,对所有结构进行优化,删去收敛到相同构象的计算结果,即获得了优化计算后的所有空间异构的分子结构。同时,也获得了所有分子结构的单点能和振动频率。至此,完成①②③步。
1.2 热力学数据及速率常数计算
基于统计热力学相关理论,物种的热力学参数（熵、焓、定压比热容）可以通过分子的配分函数获得。分子的配分函数是分子的平动配分函数、振-转配分函数及电子能配分函数的乘积。考虑多结构效应的分子配分函数计算方法如式（1）所示：
${{Q}^{\text{MS-T}}}=\underset{j=1}{\overset{N}{\mathop{\sum }}}\,{{Q}_{\text{rot},j}}\exp \left( -\frac{{{U}_{j}}}{{{k}_{\text{B}}}T} \right)Q_{j}^{\text{QH}}\underset{\tau =1}{\overset{t}{\mathop{\prod }}}\,{{f}_{j,\tau }}$ (1)
其中：T为温度;${{k}_{\text{B}}}$为波尔兹曼常数;$\tau $为耦合扭转的势函数;${{Q}^{\text{MS-T}}}$为分子振-转配分函数;j为该分子的一个空间构象;${{Q}_{\text{rot},j}}$为该分子构象的转动配分函数;${{U}_{j}}$为该构象相对于最稳定构象的能量差值;$Q_{j}^{\text{QH}}$为该构象的振动简谐配分函数;${{f}_{j,\tau }}$为该结构振-转耦合效应的影响。使用该方法获得的热力学数据用MS-T表示;若不考虑振-转耦合效应,即${{f}_{j,\tau }}$取1,则记为MS-LH;若仅考虑分子能量最低的构象的配分函数,获得的结果记为SS-HO。
基元反应的速率常数基于过渡态理论计算,并利用多结构效应进行校正。速率常数${{k}^{\text{MS-T/TST}}}$由基于过渡态理论获得的速率常数${{k}^{\text{TST}}}$与多重构象校正因子$F_{\text{act}}^{\text{MS-T}}$乘积获得：
${{k}^{\text{MS-T/TST}}}=F_{\text{act}}^{\text{MS-T}}{{k}^{\text{TST}}}$(2)
式（3）给出过渡态理论计算速率常数的方法：
${{k}^{\text{TST}}}=\kappa \frac{{{k}_{\text{B}}}T}{h}\frac{{{Q}^{\ddagger }}}{{{Q}_{\text{R}}}}\exp \left[ -{{E}_{0}}/\left( {{k}_{\text{B}}}T \right) \right]$ (3)
式中：${{Q}^{\ddagger }}$和${{Q}_{\text{R}}}$分别对应最低能量构象的过渡态、反应物的配分函数;E0为过渡态和反应物两者之间的能量差;$\kappa $为隧道效应因子。
$F_{\text{act}}^{\text{MS-T}}$作为反应的多重构象校正因子,可由式（4）得到：
$F_{\text{act}}^{\text{MS-T}}=F_{\text{TS}}^{\text{MS-T}}/F_{\text{R}}^{\text{MS-T}}$ (4)
式中：$F_{\text{R}}^{\text{MS-T}}$、$F_{\text{TS}}^{\text{MS-T}}$分别为反应物和过渡态的多重构象影响因子,其计算方法如下：
$F_{X}^{\text{MS-T}}=\frac{Q_{\text{rovib}}^{\text{MS-T}}(X)}{Q_{\text{rovib}}^{\text{SS-QH}}(X)}$ (5)
其含义是考虑多重构象效应获得的配分函数与最低能量结构获得的配分函数的比值。
2 结果和讨论
2.1 分子构象
图2为所计算的三个基元反应。OOpROH16hs是异丙醚的低温氧化产物,而OONROOH4hs、OOEROOH5hs是正丙醚的脱羟基产物。
Fig. 2 Three OOQOOH species and their reaction products图2 三个OOQOOH物种及其反应产物 |
Fig. 3 The population of distinguishable isomers of different species within each relative potential energy range图3 物种在每个相对势能范围内空间异构体的数量布居 |
Fig. 4 Molecular structure of typical transition state图4 典型过渡态分子结构 |
从本质上来讲,分子不同构象对应基团绕C—C键转动一周产生的1个或数个极小值,也称为物种的局域周期。而由于C—C键之间的耦合作用,不同C—C键转动产生的极小值有相干性。因此,在如图5所示的反应物和过渡态的局域周期数是介于2 ~ 3之间的非整数。
Fig. 5 Local periodicity Mj,τ of species图5 物种的局域周期性Mj,τ |
2.2 热力学参数
根据分子的平动、振-转运动配分函数及其电子能配分函数获得其总配分函数,从而可以直接计算该分子的热力学参数,如熵、定压比热容、焓。图6给出了所涉及物种在298 ~ 1 200 K温度范围内的配分函数。图中可见,MS-LH和MS-T方法获得的配分函数相差较小,前者比后者高24% ~ 36%。而根据SS-HO方法获得的配分函数远低于MS-T方法的结果,并且随着温度升高,差异越大。基于不同方法下的配分函数,获得了各物种的热力学参数,如表2所示。在室温下,基于MS-T方法获得的焓值比MS-LH方法高0.63 kJ/mol,当温度升高到1 200 K后,前者比后者高0.35 kJ/mol,总体而言,差异变化较小。对比而言,SS-HO方法获得的焓值比MS-T方法获得的结果低9.2 ~ 18.4 kJ/mol,具有显著差别。熵和定压比热容的对比结果与焓值的结果一致。由此可见,对于长链物种,采用MS-T方法考虑振-转非谐效应是十分必要的。
Fig. 6 The partition functions obtained using different methods图6 利用不同方法计算获得的配分函数 |
Table 2 Standard entropy, heat capacity, and relative enthalpy calculated using different methods表2 使用不同方法计算得到的标准熵、定压比热容和相对焓 |
| 物种 | 温度/K | $S_{T}^{0}$/[J/(mol∙K)] | $C_{P}^{0}$/[J/(mol∙K)] | $H_{T}^{0}$/(kJ/mol) | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MS-LH | MS-T | SS-HO | MS-LH | MS-T | SS-HO | MS-LH | MS-T | SS-HO | ||||
| OONROOH4hs | 298 | 540.66 | 542.05 | 438.32 | 213.86 | 214.83 | 196.20 | 44.50 | 45.13 | 35.62 | ||
| 300 | 542.09 | 543.47 | 439.53 | 214.78 | 215.75 | 197.16 | 44.92 | 45.59 | 36.00 | |||
| 500 | 673.90 | 675.49 | 554.65 | 306.58 | 306.12 | 289.42 | 97.32 | 98.04 | 85.02 | |||
| 800 | 840.21 | 840.67 | 707.39 | 399.85 | 394.82 | 385.45 | 204.57 | 204.49 | 187.70 | |||
| 1 200 | 1016.95 | 1 014.23 | 873.66 | 469.59 | 458.66 | 458.53 | 379.88 | 376.57 | 358.11 | |||
| OOEROOH5hs | 298 | 538.82 | 539.95 | 434.46 | 225.33 | 226.09 | 196.32 | 44.96 | 45.67 | 35.04 | ||
| 300 | 540.33 | 541.46 | 435.72 | 226.09 | 226.84 | 197.33 | 45.38 | 46.09 | 35.41 | |||
| 500 | 675.20 | 676.33 | 552.59 | 308.76 | 307.67 | 291.35 | 98.92 | 99.58 | 84.68 | |||
| 800 | 841.80 | 841.47 | 707.69 | 399.55 | 394.03 | 388.29 | 206.29 | 206.03 | 188.12 | |||
| 1 200 | 1018.33 | 1 014.73 | 876.17 | 468.96 | 458.07 | 461.59 | 381.39 | 377.79 | 359.74 | |||
| OOpROOH16hs | 298 | 521.70 | 520.03 | 434.34 | 233.66 | 236.01 | 202.10 | 43.33 | 43.58 | 35.66 | ||
| 300 | 523.29 | 521.62 | 435.64 | 234.58 | 236.93 | 203.10 | 43.79 | 44.04 | 36.08 | |||
| 500 | 663.56 | 662.35 | 555.65 | 318.55 | 317.55 | 294.69 | 99.42 | 99.84 | 86.27 | |||
| 800 | 833.52 | 830.71 | 712.12 | 404.03 | 398.30 | 388.29 | 208.88 | 208.25 | 190.17 | |||
| 1 200 | 1011.21 | 1 005.27 | 880.69 | 470.59 | 460.59 | 459.75 | 385.11 | 381.26 | 361.29 | |||
2.3 反应速率常数
图7为利用BHHLYP/cc-pVDZ方法获得的三个基元反应的势能面。其中,用能量最低的构象代表每个物种,势能面的能垒计算考虑了零点能效应。图中可见,异丙醚相关的OOQOOH生成KHP + OH反应的能垒显著低于正丙醚相关物种的反应能垒。另一方面,计算结果也表明三个基元反应是显著的放热反应,也反映了该类反应对燃料低温氧化的促进作用。
Fig. 7 The potential energy surface of the lowest energy conformation at the BHHLYP/cc-pVDZ level图7 在BHHLYP/cc-pVDZ级别下最低能量构象的势能面 |
基元反应的速率常数受到反应物和过渡态结构配分函数的共同影响。表征其综合效应的参数$F_{\text{act}}^{\text{MS-T}}$随温度的变化结果如图8所示。图中可见,对于三个反应,该参数在300 K下取值约为0.1,而后迅速降低到600 K时的0.02,此后对温度不再敏感。其含义在于,相比于常规的仅考虑反应物和过渡态的最低能量构象并将振动视为谐振的计算方法,采用更为客观准确的多结构非谐耦合（MS-T）方法获得的结果低90% ~ 98%,显示了多结构非谐效应对基元反应速率常数的突出影响。因此,对于低温加氧反应及其他长链碳氢燃料的基元反应速率常数研究中,参与反应的各物种的空间构象及非谐效应不可忽略。
Fig. 8 Multistructural torsional anharmonicity factor of reactions图8 反应的多结构扭转非谐性因子 |
图9显示了在298 ~ 1 200 K温度范围内,三个反应$\text{OOEROOH5hs}\rightleftharpoons \text{HKPE}+\text{OH}$（简称RE5）、$\text{OONROOH4hs}\rightleftharpoons \text{KHPH}+\text{OH}$(简称RN4)、$\text{OOpROOH16hs}\rightleftharpoons \text{KHPp}+\text{OH}$（简称RP16）的速率常数。虚线是不考虑多重构象效应,仅利用最低能量结构计算获得的速率常数;实线是考虑多重构象效应后获得的结果。尽管经计算表明,本文采用的密度泛函方法获得的能垒与高精度耦合簇方法相比,偏差小于1.7 kJ/mol,为了谨慎评估计算不确定度,将各种影响因素折算成能垒偏差考虑,可认为能垒的计算不确定度为 ±4.2 kJ/mol。以速率常数最大的反应RP16为例,获得了该反应速率常数随温度变化的不确定度：在300 K下为3.8倍,而温度升高到1 200 K时降到1.5倍。因此,可认为本文中速率常数的计算结果不确定度不超过4倍。
Fig. 9 High pressure limit rate constants in this work (solid lines: MS-CVT/SCT method, denoted as MS; dash lines: SS-CVT/SCT method, as SS) compared with literature data (pQ16 from ref. [10] versus RP16; RN4 from ref. [15] versus RN4; R5-MS from ref. [30] versus RE5)图9 本文计算的高压极限速率常数（实线使用MS-CVT/SCT方法计算,简称MS;虚线使用SS-CVT/SCT方法计算,简称SS）与文献中相似反应的速率常数的比较（pQ16取自文献[10],与RP16比较;RN4取自文献[15],与RN4比较;R5-MS取自文献[30],与反应RE5比较） |
对比图9中虚线与实线发现,前者是后者的10 ~ 100倍,随着温度升高,二者的差异越大,表明考虑多结构效应后,基元反应速率常数变为原来的1% ~ 0.1%。这是由于考虑多重构象效应后,反应物和过渡态的配分函数增大,而反应物构象多,受到更大的影响,增幅更为显著。
另一方面,通过对比不同OOQOOH反应的速率常数可以发现三者的差异。反应RP16的速率常数最高,RN4的次之,RN5的最低,这主要取决于它们的能垒高低：反应RP16的能垒比反应RN4的能垒约低20.9 kJ/mol,而反应RN4的能垒比反应RN5的能垒低12.5 kJ/mol。另外随着温度升高至超过1 000 K,反应RN4的速率有反超反应RP16的趋势,这是由于反应RN4的配分函数对速率常数的正向影响比RP16的更大。
在图9中对比了本文计算的RN4反应（图中黑色实线与虚线）与SERINYEL等[15]使用的正丙醚相同反应的速率常数（图中黑色点线）。SERINYEL等借用了乙醚一级加氧产物ROO氢转移反应的速率常数,用于正丙醚二级加氧产物OOQOOH氢转移裂解反应。图中可见,该反应速率与本文中不考虑多结构效应获得的计算结果较为接近（黑色虚线）,而与考虑多结构效应后的结果有明显偏差（黑色实线）。这是由于,相比于乙醚,丙醚的碳链更长,碳基团之间的振-转耦合作用更强,空间异构构象的复杂度更高,对速率常数的影响更为突出。由此也说明,开展OOQOOH氢转移裂解反应速率常数的详细计算的必要性。同时,也与XING等[30]针对二乙醚自由基二级加氧反应产物的氢转移裂解反应速率常数（图9,红色方块）以及SAKAI等[10]针对二乙醚自由基一级加氧产物ROO氢转移反应（图9,红色三角）的速率常数进行了比较。由图可见,文献中二乙醚的结果与本文中相似反应的未考虑多结构效应的计算结果十分接近。由于丙醚的OOQOOH的多结构效应较为显著,对速率常数有明显影响。
Table 3 Rate fitting of Arrhenius three parameters under high pressure limit表3 高压极限下的速率拟合阿伦尼乌斯三参数 |
| 反应 | 计算方法 | A/[cm3/(mol∙s)] | n | Ea/(J/mol) |
|---|---|---|---|---|
| RN4 | SS-CVT/SCT | 9.63 × 10-11 | 6.62 | 47 912.29 |
| MS-CVT/SCT | 6.05 × 10-17 | 8.49 | 36 651.23 | |
| RE5 | SS-CVT/SCT | 2.13 × 10-9 | 6.25 | 71 028.09 |
| MS-CVT/SCT | 3.35 | 8.31 | 58 646.86 | |
| RP16 | SS-CVT/SCT | 3.64 × 10-8 | 5.55 | 33 918.53 |
| MS-CVT/SCT | 1.92 × 10-12 | 6.84 | 26 152.66 |
3 结论
利用量子化学方法构建了三种典型丙醚二级加氧反应产物OOQOOH发生氢转移反应生成氢过氧化物醛酮类物质和OH基的势能面,结果表明,异丙醚的关联产物的反应势能垒最低。
考虑了物种的多重构象影响,获得了反应物的配分函数,进一步计算了298 ~ 1 200 K范围内的热力学数据。结果表明,考虑多重构象效应时获得的焓值比不考虑时平均高8.3 kJ/mol。熵和定压比热容的差异也较为明显。
最后,基于过渡态理论,获得了三个基元反应的速率常数,并拟合成可用于构建详细反应机理的阿伦尼乌斯形式,获得了指前因子、温度指数和活化能三个参数。对比结果表明,考虑多结构效应后,速率常数降低10 ~ 100 倍。本文提供了计算结果小于文献中利用类比法获得的结果,表明考虑多结构效应对二级加氧产物进一步裂解反应的速率常数的显著影响,获得的速率常数更具可靠性。
