English

• 油页岩是一种内嵌固体有机质的沉积岩^[1]，可以用来直接燃烧或者通过干馏得到页岩油^[2]。油页岩的主要有机成分是干酪根，是一种有机高分子聚合物，在热解过程中发生复杂的物理化学变化^[3-6]。单一实验技术研究手段很难改善现有油页岩综合开发利用技术，分子模拟在油页岩干酪根分子结构研究中处于日益重要地位^{[7, 8]}。Ru等^[9]构建了桦甸油页岩干酪根化学式为C₂₄₃H₄₀₇O₃N₂₅S₂的桦甸干酪根二维分子模型，利用量子化学和分子模拟优化该模型得到干酪根三维分子模型。王擎等^[10]表征了中国五种油页岩干酪根分子结构，Guan等^[11]进而以桦甸油页岩为例，构造出一套油页岩干酪根碳骨架同分异构体，总结出同分异构体中相对较为稳定的油页岩干酪根三维分子结构模型。Faulon等^[12]获取干酪根相关分子结构参数，构建了化学式为C₈₇₃H₈₉₉O₃₅₈N₁₃的干酪根三维分子模型。Orendt等^[13]利用量子化学和分子力学计算方法，构建了美国Green River油页岩干酪根的三维分子模型。Collell等^{[14, 15]}采用分子模拟技术，构建出化学式为C₂₄₂C₂₁₉O₁₃N₅S₂的油页岩干酪根三维分子模型，并通过分子动力学退火优化了分子结构，模拟了干酪根分子真实的存在条件，进而利用分子动力学和蒙特卡罗法模拟计算出干酪根分子结构模型的孔隙结构。Ungerer等^[16]针对三种不同成熟度的干酪根，采用分子动力学模拟计算，构建了六个干酪根三维分子结构模型。Gao等^[17]对油页岩干酪根的分子模型和结构变化进行深入的研究，得出随着干酪根热解过程的进行脂肪结构裂解明显并且在裂解残渣中形成新的芳香结构，芳香碳代替脂肪族碳并且成为残留物中的主要碳。本研究以抚顺油页岩干酪根结构模型为基础，通过分子动力学优化得到三维结构模型，计算维持其结构稳定性的能量特性，进一步采用量子力学方法计算干酪根键级、键长和静电势，分析干酪根热解前期反应，并从干酪根大分子中抽离出小分子进行碳骨架热解过程的分析。

  1.   研究方法

  1.1   分子动力学模拟参数设定

  本研究在先前研究^[18]的基础上，进一步通过分子模拟研究中国抚顺油页岩干酪根三维结构的热解特性，并预测干酪根热解过程中化学键断裂顺序。分子动力学计算使用MS模拟软件中Forcite模块，为了计算干酪根分子中氢键能，采用Dreiding力场对干酪根模型进行模拟计算。分子动力学模拟计算采用经典牛顿力学原理，通过不断改变分子中原子在空间几何中的位置，找到分子能量最低构型的过程。

  分子动力学-结构优化参数设定如下：Setup界面：Task选择Geometry Optimization，Quality选择Fine；More界面：Algorithm选择Smart，Energy、Force和Displacement是Fine精度下的默认值，为了使得迭代到最优结果，最大迭代次数Max.iterations选择50000；Energy界面：力场选择Dreiding，电荷采用当前电荷分配方式，Electrostatic和Van der Waals都选择Atom based。

  通过以上结构优化得到的并不一定是其能量最低结构构型，而是当前环境的极小能量构型，为了跳出模型构建的局部最优陷阱，获得全局能量最低构型，需要对模型进行模拟退火计算。分子动力学-模拟退火参数设定如下^[19]：Setup界面：Task选择Anneal，Quality选择Fine；More界面：Annealing cycles选择10，Initial temperature选择300 K，Mid-cycle temperature选择600 K，Heating ramps per cycle选择5，Dynamics steps per ramp选择500；Dynamic界面：Ensemble选择NVT，在微正则系综条件下，系统的粒子数N、体积V、温度T都保持恒定。在温度保持一致的情况下，通过改变分子中原子的位置打破当前能量最小的壁垒，找到全局能量最小的原子位置参数，即全局能量最优构型；Energy界面：力场选择Dreiding，电荷采用当前电荷分配方式，Electrostatic和Van der Waals都选择Atom based。

  对模拟退火后干酪根构型进一步结构优化，并进行各个干酪根模型的能量计算。能量计算参数设定如下：Setup界面：Task选择Energy任务，Quality选择Fine精度；Energy界面：力场选择Dreiding，电荷采用当前电荷分配方式，Electrostatic和van der Waals都选择Atom based。

  1.2   量子力学模拟参数设定

  对干酪根结构模型进行分子动力学结构优化及模拟退火后，为了研究深层次细微之处干酪根的化学稳定性，量子力学计算使用MS模拟软件中DMol₃模块。为了计算干酪根键长、键级电子密度等性质，在对干酪根三维结构进行量子力学结构优化的同时进行了电子密度、静电势、布局分析等性质的计算。

  量子力学模拟计算参数设定如下^[14]：Setup界面：Task选择Geometry Optimization，Quality选择Medium，对于干酪根大分子有机分子处理交换相关能方法采用GGA下的PBE方法，为了计算键长、键级性质，不勾选Use symmetry选项，同时为了使得大分子容易收敛，降低模拟计算步长Max.step size为0.2，增大最大迭代次数Max.iterations为1000；Electronic界面：自洽循环次数Max.SCF cycles设置为200，新旧电子混合Density为0.1，同时勾选Use preconditioner和Use smeaning选项，其他参数选择默认值；Propeties界面：勾选Electron density、Electrostatics和Population analysis选项。

  2.   干酪根最优化模型的能量

  抚顺油页岩干酪根模型经过分子动力学结构优化和模拟退火后，原子之间键长改变，键角、二面角等发生扭转，芳香结构层片发生侧移，二维结构模型转变为更加接近实际的三维结构模型(图 1)。经过结构优化和模拟退火后，干酪根大分子结构模型三维立体效果明显，模型中分子结构相互折叠挤压，原子之间的化学键发生明显的弯曲，键角、二面角变化明显，尤其是脂肪结构中的化学键。同时，芳香结构由于本身的稳定性，环本身基本不发生明显形变，芳香环之间电负性的影响，使得芳香结构排列相对整齐，因为芳环之间的共轭大π键对周围原子的强作用力原因。

  图 1

  

  图 1.  抚顺干酪根三维结构及原子编号

  Figure 1.  Three-dimensional structure and atomic number of Fushun kerogen

  下载: 全尺寸图片 幻灯片

  Dreiding力场总能量包括价电子能和非成键能，价电子能由键伸缩能、键角能、扭转能和反转能组成，非成键能由范德华能、库伦能和氢键能组成^[20]。抚顺油页岩干酪根模型结构优化和模拟退火后，总能量降低，键伸缩能和范德华力下降尤为明显(表 1)。由此可知，对模型稳定性作用较大的是键能和范德华能。

  表 1

  

  表 1  抚顺结构模型结构优化前后的势能变化

  Table 1.  Changes of potential energy before and after optimization of Fushun structural model

  下载: 导出CSV

  Structure	Total energy/ (kcal·mol-1)	Valence electron energy/(kcal·mol-1)					Non-bond energy/(kcal·mol-1)
  		bond energy	angle energy	torsion energy	inversion energy		van der Waals energy	electrostatic energy	hydrogen bond energy
  Initial structure	13441.84	2802.15	197.48	130.84	7.33		10286.74	17.31	0.00
  Optimized structure	652.29	89.184	219.50	51.43	0.82		352.87	-56.27	-5.25

  结构优化后干酪根的非成键能居总能量的主要地位，并且主要以范德华能为主，说明维持干酪根分子结构稳定性的主要是范德华能；本身氢键较少，结构优化后有少量氢键能有所降低，说明优化使得干酪根结构更加紧凑，羟基之间距离减少氢键能降低。干酪根优化后原子之间相对距离变短，键伸缩能降低明显；优化后的结构模型更加紧凑使得键角能略微增加，同时扭转能和反转能降低。

  3.   干酪根结构中的化学键

  3.1   分子结构RDF分析

  径向分布函数(RDF)表示从任意指定的“中心”原子出发，在半径为R处，厚度为δr的球壳内的找到其他粒子的几率，用g(r)表示^[21]。

  图 2为抚顺油页岩干酪根原子径向分布图像。碳碳间第一个峰在0.132-0.138 nm，且烯烃类碳碳双键之间键长为0.135 nm左右，该峰主要是碳碳双键之间的分布。第二个峰是0.138-0.149 nm，苯环结构或含杂原子环结构中碳原子之间键长为0.14 nm左右，该峰主要是芳香环结构碳碳键之间分布。第三个峰是0.149-0.155 nm，脂肪结构中无支链碳碳单键键长为0.15 nm左右，所以该峰主要是脂肪直链中碳原子之间的分布。第四个峰是0.155-0.16 nm，含有支链位置的碳碳单键键长为0.155 nm左右，所以该峰主要是支链处碳碳单键的分布。由峰面积大小可知，干酪根中碳骨架主要是脂肪结构，其次是芳香结构，最后含有少量的碳碳双键结构。第五个峰超出了碳原子之间成键范围，由于碳原子在干酪根中主要以极性分子存在，所以碳原子间可能产生分子间作用力范德华力。

  图 2

  

  图 2.  抚顺原子径向分布函数图像

  Figure 2.  Atomic radial distribution function image of Fushun kerogen

  (a): radial carbon-carbon distribution function; (b): radial distribution function between hydrocarbons; (c): carbon and oxygen radial distribution function; (d): hydrogen-hydrogen radial distribution function

  下载: 全尺寸图片 幻灯片

  碳氢间距离主要集中在0.11 nm左右，碳氢键键长为0.11 nm左右，碳氢之间主要是以共价键的形式存在；原子间距在0.22 nm左右的第二个峰表示碳氢间非成键分布，极性碳与极性氢之间相互吸引或排斥，可能产生范德华力作用。

  羰基、羧基的碳氧双键键长主要是0.12 nm左右，碳氧间RDF中第一个峰表示有少量碳氧以碳氧双键形式存在。由于芳香结构中大π键的作用对链接它的原子具有强烈的吸引力，芳香碳原子与氧原子之间化学键键长主要是0.14 nm左右，并且醚键中碳氧原子之间化学键键长在0.15 nm左右。因此，第二个、第三个峰依次是芳香碳与氧原子、脂肪碳与氧原子之间形成共价键成键分布规律。

  氢原子之间距离普遍大于0.17 nm，是非成键范围。由于氢原子在干酪根分子中一般都呈正价，氢原子之间产生一定的排斥力，可能也是范德华力组成的一部分。

  3.2   化学键的键长

  化学键键长指成键原子原子核间距离。一般来说，相同类型原子间键长越长，键能越低^[22]，在化学反应中更容易先断裂。

  由图 3可知，抚顺油页岩干酪根化学键键长在0.10 nm左右有个小峰，主要是氧氢键和极少的氮氢键；键长更多地集中在0.11 nm左右，是碳原子与氢原子成键范围，键长较短更为稳定；键长0.135 nm左右的第三个峰主要是脂肪链中的碳碳双键和吡咯中碳氮化学键，键长较短，结构相对碳碳单键稳定；0.14 nm左右第四个峰是除吡咯以外的芳香环内碳碳键长及芳香结构与相邻氧原子之间化学键键长；0.150 nm左右主要是脂肪结构中碳氧单键、芳香结构与脂肪链之间化学键和与烯烃结构相连接的碳碳单键，它们较一般碳碳化学键短的原因，在于氧的强非金属性和芳香结构的特性；0.155 nm左右为第六个主要峰，主要是脂肪碳之间的化学键和醚键，其中, 含有支链的碳碳化学键较一般直链碳碳化学键键长要长一些，较直链碳碳单键更为不稳定。除此之外，一些小峰如0.185、0.192 nm等处主要是碳硫单键，更为不稳定，预测在热解过程中首先会断裂生成相应的自由基。

  图 3

  

  图 3.  抚顺干酪根化学键键长分布

  Figure 3.  Chemical bond length distribution of Fushun kerogen

  下载: 全尺寸图片 幻灯片

  3.3   化学键的键级

  对分子动力学优化后的抚顺油页岩干酪根进行量子力学优化，计算其化学键键级，统计原子间化学键键级，并将除氢原子外其他原子化学键键级从小到大排序(表 2)。其中，S原子在干酪根脂肪结构中是以一个砜一个亚砜的形式存在，两个结构中S与碳原子之间化学键键级分别是0.6185、0.6319和0.7638、0.7132，键级最小；S与氧原子间双键键级都在1.1以上，相对两侧碳硫键更为稳定，热解过程中更容易断开的是S与两侧C原子间化学键，这是热解过程中较早释放含硫化合物的可能原因^[23]；含S芳香结构噻吩由于芳香结构共轭大π键的特殊性，S与周围碳原子间键级在1.05左右，热解前期化学键较稳定。

  表 2

  

  表 2  抚顺干酪根化学键键级

  Table 2.  Chemical bond level of Fushun kerogen

  下载: 导出CSV

  Chemical bond	Bond order	Chemical bond	Bond order	Chemical bond	Bond order
  C118-S119	0.6185	C265-O266	0.9054	C150-C151	0.9140
  S119-C120	0.6319	C74-C200	0.9058	C85-C172	0.9146
  S177-C178	0.7123	C97-C196	0.9062	C169-C260	0.9146
  C176-S177	0.7638	C75-C76	0.9063	C255-C256	0.9156
  C263-O264	0.8726	C223-O270	0.9064	C209-C210	0.9159
  C144-C145	0.8743	C82-C183	0.9066	C126-C127	0.9161
  C267-O268	0.8774	C253-C254	0.9070	C87-C153	0.9163
  C116-C117	0.8876	C175-C176	0.9091	C105-C242	0.9164
  C164-C165	0.8923	C214-C215	0.9093	C131-C189	0.9166
  C140-C260	0.8963	C113-C114	0.9107	C86-C87	0.9170
  C82-C137	0.9027	C159-C160	0.9120	C144-C211	0.9172
  C82-C161	0.9027	C144-C261	0.9127	C83-C188	0.9180
  C131-C132	0.9030	C197-C198	0.9130	C148-C149	0.9186
  C105-C106	0.9053	C164-C263	0.9137	C196-C250	0.9187

  直接连接在芳香结构上的碳氧单键键级都在1.0左右，较为稳定，原因可能在于芳环的π-共轭效应。但是，氧原子与另一侧的氢原子间化学键稳定性就会遭到弱化，键级为0.8左右。作为脂肪链终端的氧与碳之间化学键键级一般都在0.9以下，氢氧之间键级为0.75左右，因此, 在热解过程中可能更容易断裂氢氧键生成H·自由基，使得氧未配对电子向碳氧中间靠拢，相邻C-O化学键更为稳定。同时，间位C-C键由于电子密度降低化学键稳定性减弱，预测容易断裂。此时，碳原子两侧两个半键重新组合与原来碳氧键一起生成碳氧双键，进一步氧化生成CO，因为C-O键同样不稳定，亦可能断裂生成·OH自由基。在羧基结构中心碳原子与脂肪碳间化学键键级为0.85左右，羧基结构中的氢与其相连的氧原子间化学键键级为0.769，在热解过程此处氢氧化学键和结构另一侧碳碳单键可能会同时断裂，两个半键重新组合生成CO₂。由以上氧原子存在状态及键级可以推测，热解初期氧原子可能一般以CO、CO₂和H₂O气体形式挥发，与先前实验研究结果相一致^[24]。

  作为季碳、叔碳的中心碳原子与周围碳原子的化学键键级一般在0.9左右，较脂肪直链中碳碳单键键级低。含有杂原子等官能团在干酪根热解初期断裂后，随着热解过程的进行，碳骨架中含有支链的C-C键相对直链更不稳定较早断裂。对于含有碳碳双键的结构，形成双键的碳原子之间化学键更为稳定，但是双键结构间位上的碳碳键却遭受弱化，键级一般也在0.9左右，可能由于双键强烈的电负性吸引电子导致间位C-C键共用电子稀薄引起的。除芳香硫中碳硫键、吡咯中碳氮键键级在1.0左右。

  3.4   化学键的静电势级

  通过基于密度泛函理论的量子力学计算抚顺油页岩干酪根结构模型以及电子密度和电势，将电势赋予电子密度图上，形成结构电势分布图，具体见图 4。由图 4可知，电势较低的结构主要集中在含氧或含氮官能团上，包括羧基、芳醚、脂醚、羟基、亚砜、砜和含氮芳香环。氧和氮的原子半径较小非金属性强，更容易吸引电子表现出负极性。碳碳双键电势同样较低，在碳原子杂化理论中，形成单键的碳原子电子轨道都参与了杂化，是sp³杂化轨道，而含有双键的烯烃碳原子是sp²杂化，sp²杂化碳原子比sp³杂化的电负性要大，由于电子云的偏移电荷在分子结构分布上不均匀更容易发生极化，使得化学键两侧一端带正电，一端带负电，产生较多的静电力。如表 1所示，抚顺油页岩干酪根中库伦能占非键能的13%左右。

  图 4

  

  图 4.  抚顺干酪根电势

  Figure 4.  Fushun kerogen potential

  下载: 全尺寸图片 幻灯片

  电势低的地方容易发生亲核反应，电势高的地方容易发生亲电反应。电势较高区域集中在氢原子上，并主要是与氧原子或氮原子直接相连的氢原子，还包括亚砜、砜相邻碳原子上的氢原子，这些氢原子在热解过程中断裂以后生成的·H自由基更容易与电势较低的部位发生亲电反应，如对羟基部位的氧原子发起攻击，使得C-O键断裂并与氢结合生成H₂O，对羰基结构C-C键进行攻击并与碳重新组合生成甲醛并进一步氧化得CO等；还包括在热解后期共轭二烯烃中供电子基团和单烯烃中吸电子基团发生双烯合成，生成环状化合物。

  4.   干酪根碳骨架结构稳定性

  油页岩干酪根碳骨架结构由于其复杂性，直接进行整体研究缺乏系统性与完整性。通过小分子模型来分析大分子化合物，被认为是一种研究大分子化合物化学性质的有效方法^[22]。因此，本节将碳骨架中主要的碳结构及主要官能团抽离出来详细分类研究。化学键稳定性受到周围化学环境的强烈影响，同一种化学键在不同化学环境中有不同的稳定性。本节研究了碳骨架中主要的C-C键和C-O键在不同官能团和不同连接方式下的化学性质，键级越大越稳定。因此，以键级为标准分析化学键的稳定性是合理可靠的。

  4.1   碳骨架中C-C键稳定性

  图 5为甲基取代基个数与位置对C-C键稳定性影响。

  图 5

  

  图 5.  甲基取代基个数与位置对C-C键稳定性的影响

  Figure 5.  Influence of number and position of methyl substituents on stability of C-C bonds

  下载: 全尺寸图片 幻灯片

  由图 5可知，通过改变脂肪链中的支链个数与位置得到相邻C-C化学键键级变化，与支链相连接的C-C化学键键级为0.888-0.975，平均键级为0.931。由图 5((a)-(c))可知，有一个甲基时相邻C-C键键级为0.975，有两个、三个时键级分别是0.953、0.937，即碳直链上增加一个支链邻位上C-C键，键级平均下降0.019。但是，间位上的C-C键键级分别是0.974、0.966和0.970，相对没有支链时键级0.975变化不大，说明支链个数对间位和更远处上的C-C键只有微小影响。

  图 5(d)相对图 5(c)中C-C键左侧增加一个支链，键级为0.911，较图 5(c)中0.937下降0.026。图 5(e)相对图 5(c)是在研究C-C键的间位上增加一个支链，键级为0.932，相对0.937只下降0.005，验证了上述支链对间位上C-C间有微小影响的结论。图 5(f)相对图 5(d)中C-C键左侧又增加一个支链，键级由0.911下降到0.888，降幅为0.023，相对同侧增加一个支链键级下降0.026，下降幅度降低。相同条件下，两侧增加支链相对一侧单独增加来说更能降低化学键键级，原因在于化学键两侧碳结构电子分布效应所致的空间位阻，同时在相同化学环境中C-C键处于终端时比中间具有更大的键级，较处在中间时稳定。

  图 6为不同官能团对C-C键稳定性的影响。

  图 6

  

  图 6.  不同官能团对C-C键稳定性的影响

  Figure 6.  Effect of different functional groups on stability of C-C bonds

  下载: 全尺寸图片 幻灯片

  由图 6(a)可知，当脂肪直链中含有C=C键时会对周围化学键产生强烈影响。根据图 6(a)与图 5(a)相同位置化学键相比可知，C=C键邻位上的化学键键级由0.975增加到了0.986，增长幅度为0.011；间位上的C-C键键级由0.974降低到0.955，降低幅度为0.019。C=C双键中包括一个π键和一个与C-C单键一样的σ键，在碳原子杂化理论中，C-C单键是sp³杂化轨道，而C=C双键是sp²杂化轨道。所以，C=C双键两侧中两个碳原子都会存在剩余的一个未参与杂化的p轨道，他们的p轨道对称轴都垂直于σ轨道对称轴所在平面形成π键，由于π电子的离域效应，双键π轨道与邻位原子的p轨道在侧面相互交盖，构成p，π-共轭体系，分子内能降低更加稳定。

  由图 6((b)、(c))可知，脂肪链中碳原子被N原子取代后，邻位和间位的C-C键基本没有变化，可能由于C、N在元素周期表中相邻位置，并且C-N键跟C-C键都属于σ键，具有相似化学性质。如图 6(d)-6(g)所示，含硫官能团对相邻C-C键影响较大，相邻C-C键键级为0.915-0.952，其中，含硫官能团对相邻C-C键削弱最明显。除此之外，作者发现这几类含硫官能团不仅对邻位C-C化学键有削弱影响，对整条直链中间化学键都有明显的削弱作用，且随着距离含硫官能团的变远削弱幅度减小。对于含氧官能团对相邻C-C键同样有削弱影响，其中羰基最明显，邻位C-C键键级为0.918(图 6(i))；两类含氧官能团作用稍小，邻位C-C键键级分别为0.947和0.951(图 6(h)、6(j))；脂肪醚、芳醚对C-C键基本没有影响(图 6(k)、6(l))。

  图 7为芳香结构对C-C键稳定性影响。

  图 7

  

  图 7.  芳香结构对C-C键稳定性的影响

  Figure 7.  Effect of aromatic structure on stability of C-C bonds

  下载: 全尺寸图片 幻灯片

  由图 7可知，脂肪链上连接的芳香结构对C-C键稳定性同样具有强烈的影响。六元环苯和吡啶对脂肪链与芳香环间直接相连C-C化学键具有削弱的作用，除图 7(d)结构中键级为0.960外，其他键级都在0.94左右(图 7((a)-(d)))。芳香结构对与六元环间接相连的C-C键影响不大，但是对更远一个单位距离的C-C键稳定性有削弱作用，但不如与环直接相连的键明显。对于五元环吡咯，脂肪链与吡咯的连接方式并没有对直接相连C-C键稳定性产生影响，两种连接方式键级均在0.96左右，与六元环情况相一致的是距离环第三位置的C-C键稳定性同样被削弱。五元环噻吩同样对直接相连C-C键稳定性具有削弱作用，与其他环结构不同的是，噻吩对相连的整条脂肪链C-C键都有一定的削弱作用，与硫醚结构对C-C键稳定性影响相类似，可能是硫本身的原子半径较大，与其他原子成键时电子云的推动作用有关系(图 7((g)、(h))。

  4.2   碳骨架中C-O键稳定性

  在干酪根碳骨架中，C-O键主要与C-C双键和芳香结构为主。图 8为碳碳双键对C-O键稳定性影响。由图 8(a)可知，C-O键在脂肪直链中的键级平均为0.850。由图 8(b)可知，当O与C=C键中一侧碳原子成键时，键级为0.890，比在脂肪直链中键级0.850高0.040，C=C键对相邻C-O键稳定性具有增强作用，因为C=C双键π轨道与邻位原子的p轨道在侧面相互交盖，构成p，π-共轭体系，分子内能降低更加稳定。由图 8((b)-(d))可知，C=C键对其间位上C-O键稳定性是具有削弱作用，从距离自身三个以上化学键距离的键几乎没有影响。

  图 8

  

  图 8.  碳碳双键对C-O键稳定性的影响

  Figure 8.  Effect of carbon-carbon double bond on stability of C-O bond

  下载: 全尺寸图片 幻灯片

  图 9为芳香结构对C-O键稳定性影响。由图 9可知，当氧原子与芳香结构直接或者间接相连时，氧原子与邻近芳香结构一侧碳原子间形成化学键键级较远侧C-O键中高，与远离芳香结构一侧的碳原子形成的C-O键键级低，稳定性差。由图 9中(a)分别与(b)、(d)与(e)、(g)与(h)、(g)与相比较，发现芳香结构对直接相连的氧原子形成的C-O键比间接相连时稳定性影响更大，当氧原子与芳香结构距离更远时，如图 9((c)、(f)、(i)、(l))，芳香结构对C-O键基本没有影响。

  图 9

  

  图 9.  芳香结构对C-O键稳定性的影响

  Figure 9.  Effect of aromatic structure on the stability of C-O bonds

  下载: 全尺寸图片 幻灯片

  4.3   碳骨架中自由基对邻近化学键稳定性

  由上述化学键稳定性分析可知，所谓较弱的化学键并不是某一类型固定的化学键，还受周围化学环境的影响，其中，包括各类官能团种类和不同连接方式。除此之外，在油页岩干酪根热解过程中会产生一定的自由基，碳骨架上自由基种类和连接方式同样影响周围化学键的稳定性。因此，进一步研究热解过程中产生的不同类自由基在不同连接方式下对碳骨架中C-C键稳定性的影响。

  图 10为自由基对C-C键稳定性影响。由图 10(a)可知，伯碳自由基邻位的C-C键键级为1.024，较图 5(a)相同环境中C-C键键级0.975提高了0.049，而其间位上键级为0.926，相对图 5(a)中相同位置化学键键级下降了0.049，说明伯碳自由基能够增强邻位C-C键稳定性，而降低间位C-C键稳定性。如图 10((b)、(c))所示，仲碳叔碳自由基对邻位C-C都有加强稳定性的作用，邻位化学键键级平均达到1.0，但是对间位C-C键影响不大。支链对相邻C-C键具有弱化作用(图 5((b)、(c))，C-C化学键键级提高是受到自由基及支链共同影响(图 10((b)、(c)，因此, 实际自由基对此处C-C键稳定性强化作用更为强大。因为碳自由基具有一个空p轨道，与碳自由基相连的烷基碳原子上的H原子中σ键轨道可以与此空轨道在一定程度上相互交盖而引起离域效应，使得σ电子离域并扩散到空p轨道上，形成σ，p-共轭，并且碳自由基处支链数越多越稳定。

  图 10

  

  图 10.  自由基对C-C键稳定性的影响

  Figure 10.  Effect of free radicals on stability of C-C bonds

  下载: 全尺寸图片 幻灯片

  由图 10((d)、(e)、(f))所示，氧自由基邻位的C-O键与间位的C-C键键级随相连C原子支链个数增多而下降，说明氧原子越处于支链多的位置本身C-O键越不稳定，同时导致支链与主链相连C-C键亦不稳定，C-C键不稳定，一方面，由于氧自由基的影响；另一方面，支链越多，本身支链处C-C键不稳定(图 5)。

  由图 10((g)-(l))可知，S和N自由基与氧自由基具有极其类似的影响作用，都是杂原子与邻位碳化学键及其间位C-C键的稳定性随着支链个数增多而下降。

  4.4   碳骨架热解中后期化学键断裂顺序预测

  通过上述分析，键级较低的化学键并不是固定不变的。各类化学键类型中，同类化学键的稳定性均有一定程度的波动，用某一种化学键的类型来直接预测干酪根中化学键的断裂优先级是不准确的。因为化学键稳定性还受周边的化学环境的影响，包括不同官能团和自由基。因此，通过构造与干酪根大分子化学环境相同的小分子来预测干酪根大分子的稳定性是可行的。

  干酪根碳骨架的热解首先是化学键的断裂。为了进一步探求干酪根反应活性位点预测热解反应过程，计算干酪根化学键断裂的难易程度是研究干酪根热解过程的关键。通过计算不同官能团和自由基对C-C、C-O键稳定性的影响，作者得出，脂肪链中支链只对邻位C-C键稳定性有削弱影响；C=C键对邻位C-C键稳定性具有增强效果，对间位有削弱作用；芳香结构对直接相连接C-C具有削弱作用，对含有氧原子的脂肪链中近侧C-O键具有增强效果，远侧具有削弱作用；干酪根中自由基对相邻的化学键同样有很大的影响，处于脂肪直链末端的碳自由基对邻位C-C键稳定性有增强效果，对间位的C-C键具有严重削弱作用，对于支链处碳原子自由基则会加强邻位的C-C键，使其相对稳定，而杂原子自由基增强杂原子与碳间化学键稳定性，削弱邻位C-C键稳定性。

  干酪根碳骨架中稳定性较差的化学键优先发生断裂。首先是含有芳香结构的并且远离芳香结构一侧的脂肪链中的C-O键，因为C-O键相对C-C键稳定性差，在芳香结构的影响后稳定性更差；其次为脂醚键的断裂，尤其是在碳碳双键间位上或者有支链的邻位上，因为双键或支链对C-O键稳定性的削弱作用；然后是脂肪链中的支链处和有碳自由基的终端碳碳键，由于热解过程中产生一定的自由基，自由基同样影响相邻化学键的稳定性。

  通过以上分析，终端有碳自由基结构容易发生连锁反应，而支链处某一个C-C键断裂后自由基加强了原来其他支链的稳定性，不会继续断裂，因此, 热解中碳链支化度不会明显下降反而由于直链的连锁反应导致亚甲基碳变少，使得支化度提高，与实验情况一致^[24]。然后是脂肪链与芳香结构连接的化学键和双键间位上化学键的断裂，将芳香结构与脂肪结构分开。随着热解温度的提高，含有自由基的脂肪支链处化学键继续断裂生成直链，支化度降低。最后为稳定性相对较高的化学键的断裂和一些自由基的重组，主要为含有杂原子的芳香结构的开环反应和带有自由基的烯烃结构的芳构化反应。在整个过程中，自由基含量先上升后下降，与实验结果一致^[25]。

  5.   结论

  经过结构优化和模拟退火后，油页岩干酪根大分子结构模型三维立体效果明显，模型中分子结构相互折叠挤压，化学键发生明显的扭转和弯曲，键角、二面角变化明显，尤其是脂肪结构。因为芳环中共轭大π键的强作用力导致其结构的稳定性，环基本不发生明显形变，又由于芳香环之间电负性的影响，使得芳香结构排列相对整齐。

  油页岩干酪根结构模型优化后，结构总能量明显降低。优化后干酪根的非成键能居总能量的主要地位，且以范德华能为主，说明维持干酪根结构稳定性的主要是范德华能。干酪根优化后更加紧凑，原子间相对距离变短，键伸缩能降低明显，优化过程中键角改变使键角能略微增加，但扭转能和反转能降低。

  根据油页岩干酪根化学结构特性预测，干酪根分子热解初期首先是亚砜、砜官能团与脂肪碳链间化学键的断裂生成SO₂等挥发性气体，其次是位于脂肪碳链终端的碳氧键的断裂生成·OH自由基与芳香结构相邻氧原子上断裂的·H自由基组合生成H₂O。

  随着热解温度的升高，羰基结构中碳原子与两侧碳原子间化学键同时断裂，生成CO或者羰基中心碳两侧半键与·H自由基结合生成CH₂O中间体，进而氧化成CO。在羧基结构中，首先断裂的是氧氢键，其次是碳碳键，断裂的两个半键重组生成CO₂，与实验初期干酪根热解会生成少量CO₂、CO、SO₂和H₂O相匹配。

  随着热解过程的进行，一些不稳定氧、硫官能团逐渐脱离，剩下的主框架结构的热分解是由氧桥的解离引起的。首先是含有芳香结构的并且远离芳香结构一侧的脂肪链中的C-O键的断裂；其次为脂醚键的断裂；然后是脂肪链中的支链处和有碳自由基的终端碳碳键的断裂；接着是脂肪链与芳香结构连接的化学键和双键间位上化学键的断裂；然后含有自由基的脂肪支链处和直链化学键继续断裂；最后为稳定性相对较高的化学键的断裂和一些自由基的重组。

• 1. [1]

     战金辉, 赖登国, 许光文. 油页岩:固体石油[J]. 科学世界, 2016(12): 68-73. ZHAN Jin-hui, LAI Deng-guo, XU Guang-wen. Oil shale:Solid petroleum[J]. Sci World, 2016, (12):  68-73.

  2. [2]

     朱煜凯, 张宇哲. 国外油页岩资源的利用分析[J]. 化工设计通讯, 2016,42,(1): 152-152. ZHU Yu-kai, ZHANG Yu-zhe. Analysis of utilization of foreign oil shale resources[J]. Chem Des Commun, 2016, 42(1):  152-152.

  3. [3]

     王忠鑫, 赵丹丹. 油页岩资源开发潜力评价研究与实践[J]. 露天采矿技术, 2016,31,(5): 83-87. WANG Zhong-xin, ZHAO Dan-dan. Research and practice of oil shale resource development potential evaluation[J]. Ope Min Technol, 2016, 31(5):  83-87.

  4. [4]

     王擎, 许祥成, 迟铭书, 张宏喜, 崔达, 柏静儒. 干酪根组成结构及其热解生油特性的红外光谱研究[J]. 燃料化学学报, 2015,43,(10): 1158-1166. doi: 10.3969/j.issn.0253-2409.2015.10.002WANG Qing, XU Xiang-cheng, CHI Ming-shu, ZHANG Hong-xi, CUI Da, BAI Jing-ru. FT-IR study on composition of oil shale kerogen and its pyrolysis oil generation characteristics[J]. J Fuel Chem Technol, 2015, 43(10):  1158-1166. doi: 10.3969/j.issn.0253-2409.2015.10.002

  5. [5]

     王越, 高燕, 白向飞, 武琳琳. 桦甸油页岩有机岩相特征及其富集特性[J]. 燃料化学学报, 2016,44,(3): 321-327. WANG Yu, GAO Yan, BAI Xiang-fei, WU Lin-lin. Characteristics and enrichment characteristics of organic lithofacies of huadian oil shale[J]. J Fuel Chem Technol, 2016, 44(3):  321-327.

  6. [6]

     KUMAR R, BANSAL V, BADHE R M, MADHIRA I S S, SUGUMARAN V, AHMED S, CHRISTOPHER J, PATEL M B, BASU B. Characterization of Indian origin oil shale using advanced analytical techniques[J]. Fuel, 2013, 113:  610-616. doi: 10.1016/j.fuel.2013.05.055

  7. [7]

     LⅡV S, KAASIK M. Trace metals in mosses in the estonian oil shale processing region[J]. J Atmos Chem, 2004, 49(1):  563-578.

  8. [8]

     ZHOU B, SHI L, LIU Q, LIU Z. Examination of structural models and bonding characteristics of coals[J]. Fuel, 2016, 184:  799-807. doi: 10.1016/j.fuel.2016.07.081

  9. [9]

     RU X, CHENG Z Q, SONG L H, WANG H Y, LI J F. Experimental and computational studies on the average molecular structure of Chinese Huadian oil shale kerogen[J]. J Mol Struct, 2012, 1030(4):  10-18.

  10. [10]

      王擎, 黄宗越, 迟铭书, 石聚欣, 王智超, 隋义. 油页岩干酪根化学结构特性分析[J]. 化工学报, 2015,66,(5): 1861-1866. WANG Qing, HUANG Zong-yue, CHI Ming-shu, SHI Ju-yi, WANG Zhi-chao, SUI Yi. Analysis of chemical structure of oil shale kerogen[J]. CIESC J, 2015, 66(5):  1861-1866.

  11. [11]

      GUAN X H, LIU Y, WANG D, WANG Q, CHI M S, LIU S, LIU C G. Three-dimensional structure of huadian oil shale kerogen model:An experimental and theoretical study[J]. Energy Fuels, 2015, 29(7):  4122-4136. doi: 10.1021/ef502759q

  12. [12]

      FAULON J L, VANDENBROUCKE M, DRAPPIER J M, BEHAR F, ROMERO M. 3D chemical model for geological macromolecules[J]. Org Geochem, 1990, 16(4):  981-993.

  13. [13]

      ORENDT A M, PIMIENTA S O, BADU S R, SOLUM M S, PUGMIRE R J, FACELLI J C. Three-dimensional structure of the siskin green river oil shale kerogen model:A comparison between calculated and observed properties[J]. Energy Fuels, 2013, 27(2):  702-710. doi: 10.1021/ef3017046

  14. [14]

      COLLELL J, UNGERER P, GALLIERO G, YIANNOURAKOU M, MONTEL F, PUJOL M. Molecular simulation of bulk organic matter in type Ⅱ shales in the middle of the oil formation window[J]. Energy Fuels, 2014, 28(12):  7457-7466. doi: 10.1021/ef5021632

  15. [15]

      COLLELL J, GALLIERO G, VERMOREL R, UNGERER P, YIANNOURAKOU M, MONTEL F, PUJOL M. Transport of multicomponent hydrocarbon mixtures in shales organic matter by molecular simulations[J]. J Phys Chem C, 2015, 119(39):  22587-22595. doi: 10.1021/acs.jpcc.5b07242

  16. [16]

      UNGERER P, COLLELL J, YIANNOURAKOU M. Molecular modeling of the volumetric and thermodynamic properties of kerogen:In fluence of organic type and maturity[J]. Energy Fuels, 2015, 29(1):  91-105. doi: 10.1021/ef502154k

  17. [17]

      GAO Y, ZOU Y R, LIANG T, PENG P. Jump in the structure of Type Ⅰ kerogen revealed from pyrolysis and 13C DP MAS NMR[J]. Org Geochem, 2017, 112:  105-118. doi: 10.1016/j.orggeochem.2017.07.004

  18. [18]

      王擎, 程枫, 潘朔. 油页岩干酪根化学键浓度与能量密度研究[J]. 燃料化学学报, 2017,45,(10): 1209-1218. doi: 10.3969/j.issn.0253-2409.2017.10.008WANG Qing, CHENG Feng, PAN Shuo. Study on chemical bond concentration and energy density of oil shale kerogen[J]. J Fuel Chem Techol, 2017, 45(10):  1209-1218. doi: 10.3969/j.issn.0253-2409.2017.10.008

  19. [19]

      ZHANG Z, JAMILI A. Modeling the Kerogen 3D Molecular Structure[C]//SPE/CSUR Unconventional Resources Conference. 2015.

  20. [20]

      马延平. 柳林3#煤的超分子构建及分子模拟[D]. 太原: 太原理工大学, 2012.MA Yan-ping. The supramolecular construction and molecular simulation of Liulin 3# coal[D]. Taiyuan: Taiyuan University of Technology, 2012.

  21. [21]

      张世良, 戚力, 高伟, 冯士东, 刘日平. 分子模拟中常用的结构分析与表征方法综述[J]. 燕山大学学报, 2015,39,(3): 213-220. ZHANG Shi-liang, WEI Li, GAO Wei, FENG Shi-dong, LIU Ri-ping. A Summary of structural analysis and characterization methods commonly used in molecular simulation[J]. J Yanshan Univ, 2015, 39(3):  213-220.

  22. [22]

      GUAN X H, WANG D, WANG Q, CHI M S, LIU C G. Estimation of various chemical bond dissociation enthalpies of large-sized kerogen molecules using DFT methods[J]. Mol Phys, 2016, 114(11):  1705-1755. doi: 10.1080/00268976.2016.1143983

  23. [23]

      WANG Q, LIU Q, WANG Z C, LIU H P, BAI J R, YE J B. Characterization of organic nitrogen and sulfur in the oil shale kerogens[J]. Fuel Process Technol, 2017, 160:  170-177. doi: 10.1016/j.fuproc.2017.02.031

  24. [24]

      王擎, 张岩, 迟铭书. 干酪根热解过程中的热解特性分析[J]. 石油学报(石油加工), 2017,33,(3): 507-514. WANG Qing, ZHANG Yan, CHI Ming-shu. Analysis of pyrolysis characteristics in pyrolysis of kerogen[J]. Acta Pet Sin(Pet Process Sect), 2017, 33(3):  507-514.

  25. [25]

      李慧莉, 邱楠生, 金之钧, 何治亮, 朱映康. 利用干酪根自由基浓度反演碳酸盐岩地层热历史[J]. 石油与天然气地质, 2005,26,(3): 337-343. doi: 10.11743/ogg20050312LI Hui-li, QIU Nan-sheng, JIN Zhi-jun, HE Zhi-liang, ZHU Ying-kang. Inversion of thermal history of carbonate formation using kerogen radical concentration[J]. Oil Gas Geol, 2005, 26(3):  337-343. doi: 10.11743/ogg20050312

• 

  图 1  抚顺干酪根三维结构及原子编号

  Figure 1  Three-dimensional structure and atomic number of Fushun kerogen

  下载: 全尺寸图片 幻灯片
  

  图 2  抚顺原子径向分布函数图像

  Figure 2  Atomic radial distribution function image of Fushun kerogen

  (a): radial carbon-carbon distribution function; (b): radial distribution function between hydrocarbons; (c): carbon and oxygen radial distribution function; (d): hydrogen-hydrogen radial distribution function

  下载: 全尺寸图片 幻灯片
  

  图 3  抚顺干酪根化学键键长分布

  Figure 3  Chemical bond length distribution of Fushun kerogen

  下载: 全尺寸图片 幻灯片
  

  图 4  抚顺干酪根电势

  Figure 4  Fushun kerogen potential

  下载: 全尺寸图片 幻灯片
  

  图 5  甲基取代基个数与位置对C-C键稳定性的影响

  Figure 5  Influence of number and position of methyl substituents on stability of C-C bonds

  下载: 全尺寸图片 幻灯片
  

  图 6  不同官能团对C-C键稳定性的影响

  Figure 6  Effect of different functional groups on stability of C-C bonds

  下载: 全尺寸图片 幻灯片
  

  图 7  芳香结构对C-C键稳定性的影响

  Figure 7  Effect of aromatic structure on stability of C-C bonds

  下载: 全尺寸图片 幻灯片
  

  图 8  碳碳双键对C-O键稳定性的影响

  Figure 8  Effect of carbon-carbon double bond on stability of C-O bond

  下载: 全尺寸图片 幻灯片
  

  图 9  芳香结构对C-O键稳定性的影响

  Figure 9  Effect of aromatic structure on the stability of C-O bonds

  下载: 全尺寸图片 幻灯片
  

  图 10  自由基对C-C键稳定性的影响

  Figure 10  Effect of free radicals on stability of C-C bonds

  下载: 全尺寸图片 幻灯片

  表 1  抚顺结构模型结构优化前后的势能变化

  Table 1.  Changes of potential energy before and after optimization of Fushun structural model

  Structure	Total energy/ (kcal·mol-1)	Valence electron energy/(kcal·mol-1)					Non-bond energy/(kcal·mol-1)
  		bond energy	angle energy	torsion energy	inversion energy		van der Waals energy	electrostatic energy	hydrogen bond energy
  Initial structure	13441.84	2802.15	197.48	130.84	7.33		10286.74	17.31	0.00
  Optimized structure	652.29	89.184	219.50	51.43	0.82		352.87	-56.27	-5.25

  下载: 导出CSV

  表 2  抚顺干酪根化学键键级

  Table 2.  Chemical bond level of Fushun kerogen

  Chemical bond	Bond order	Chemical bond	Bond order	Chemical bond	Bond order
  C118-S119	0.6185	C265-O266	0.9054	C150-C151	0.9140
  S119-C120	0.6319	C74-C200	0.9058	C85-C172	0.9146
  S177-C178	0.7123	C97-C196	0.9062	C169-C260	0.9146
  C176-S177	0.7638	C75-C76	0.9063	C255-C256	0.9156
  C263-O264	0.8726	C223-O270	0.9064	C209-C210	0.9159
  C144-C145	0.8743	C82-C183	0.9066	C126-C127	0.9161
  C267-O268	0.8774	C253-C254	0.9070	C87-C153	0.9163
  C116-C117	0.8876	C175-C176	0.9091	C105-C242	0.9164
  C164-C165	0.8923	C214-C215	0.9093	C131-C189	0.9166
  C140-C260	0.8963	C113-C114	0.9107	C86-C87	0.9170
  C82-C137	0.9027	C159-C160	0.9120	C144-C211	0.9172
  C82-C161	0.9027	C144-C261	0.9127	C83-C188	0.9180
  C131-C132	0.9030	C197-C198	0.9130	C148-C149	0.9186
  C105-C106	0.9053	C164-C263	0.9137	C196-C250	0.9187

  下载: 导出CSV
•