单一应力下义马原煤热解产物演变规律及动力学模拟*

邓奇根 刘朝思 姚萌萌 李 帅

(1.河南理工大学安全科学与工程学院,454003 河南焦作;2.河南省瓦斯地质与瓦斯治理重点实验室-省部共建国家重点实验室培育基地,454003 河南焦作;3.煤炭安全生产河南省协同创新中心,454003 河南焦作)

0 引 言

利用分子模拟技术表征煤体结构,有助于从微观分子尺度了解煤体结构演化过程[1-2],且该技术已被证实是研究煤体结构的有效手段之一[3]。自从FUCHS等[4]基于化学统计分析方法提出第一个煤分子模型,煤大分子模型构建被学者广泛应用,目前已构建的煤分子模型约达150种,不少学者依据分子模型开展了煤热解机理、C/H/O等元素迁移、升温速率及不同温度对热解产物的影响等相关研究。

煤热解是煤热化学转化的一个复杂过程,通过气相色谱法、热重法等研究高压升温过程中煤热解产物变化和失重变化的宏观方法[5]已被广泛运用。在微观研究方面,van DUIN et al[6]运用过程可视化的分子动力学模拟技术,极大缩短了实验时间;WANG et al[7]采用红外光谱和密度泛函理论方法对禹州褐煤建模,将热重分析结果和反应分子动力学模拟热解结果对比,验证了模型的合理性;冯炜等[1]通过对枣泉煤进行分子模型构建和热解模拟,发现热解过程中C15以下碎片种类较多,且通过追踪产物CO2,获得了3种不同的形成路径;徐芳[8]以霍林河褐煤为主体,研究了等温和非等温热解反应分子动力学模拟过程中热解产物变化及微观反应机理,结果表明等温热解时,2 200 K以下煤的一次分解反应占主导,2 200 K以上煤的二次分解反应占主导;非等温热解时,焦油的最佳产量受热解温度和升温速率的共同影响;QIU et al[9]采用ReaxFF MD方法模拟褐煤在O2/CO2气氛下燃烧的微观反应行为,并对原子运动轨迹、键断裂及形成进行追踪。

目前,有关褐煤热解的微观分子研究中,相对缺少对S元素的演化分析,同时煤热解是一个复杂过程,除了受温度和时间的影响,应力作用对煤裂解也有一定影响。王登科等[10]通过CT扫描系统可视化研究受载煤体三维裂隙结构域的演化规律,结果表明围压是影响受载煤体失稳破坏过程的重要因素;焦振华等[11]对煤样进行不同应变率和不同应力状态下的冲击压缩实验,结果表明应变率对煤样应力应变和破坏模式有显著影响;宋昱等[12]通过技术手段分析了次高温条件下应力和应变速率对煤中缩合芳环结构有序性的影响,揭示了煤缩合芳环结构有序化演化规律。但是目前对煤体的应力研究,主要是通过实验仪器分析煤体裂隙结构变化和应力-应变响应规律,缺少对煤体微观分子裂解和重组的可视化过程研究。因此,本研究基于常规褐煤热解流程,选取含硫褐煤--义马原煤,运用ReaxFF MD方法,在单一应力(外力)作用下分析煤热解产物演化规律,以期为后续不同应力大小和方向作用下的褐煤分子尺度热解反应原理及热解产物变化规律研究提供参考。

1 分子模型及优化

义马原煤属于含硫低阶褐煤,分子式为C178H146O32N2S2,单个分子的3D模型如图1所示[13]。

利用Materials Studio软件DMol3模块对煤分子进行红外光谱分析,为检验模型合理性,比较了红外光谱的实验结果和模拟结果[13],如图2所示。由图2可以看出,红外光谱实验结果和模拟结果的基本趋势大致相同,但个别峰的位置和强度有差异,这些差异多是官能团或者氢键伸缩振动引起的偏移[14];2 980 cm-1处和1 755 cm-1处的峰,分别是羟基O-H键和羧基C-O键的伸缩振动(见图3)。因此,可认为模型合理。

图3 官能团伸缩振动

对煤分子进行热解模拟之前,先进行大分子构建和预处理。利用Forcite模块,选用CompassII力场,对基础模型进行几何优化和退火分析,找出能量最稳定状态;采用Amorphous Cell模块对能量最稳态模型进行晶胞构建,密度设为0.5 g/cm3[15],晶胞包含5个稳态碳链;在Gulp模块上选用ReaxFF 6.0力场进行预处理;将预处理的模型导入LAMMPS软件[16],并通过扩大盒子边界,使模型分子都在盒子内;最终盒子边长为6.535 nm×5.268 nm×5.487 nm,采用周期性边界作为约束条件,得到模型优化及预处理结果(见图4)。

图4 模型预处理过程

2 ReaxFF MD方法模拟煤热解

对导入的LAMMPS预处理模型,选用含C,H,O,N,S元素的ReaxFF力场[17]进行最小化处理,使分子间能量降到最低,得到的模型压力稳定在109.78 GPa附近。研究表明,2 500 K是煤热解合适的模拟温度[15],应力作用会促进煤分子裂解,因此在5×109s-1应变率下,将模拟终温设为1 000 K,1 500 K,2 000 K,2 500 K和3 000 K,进行应力作用下煤热解。在ReaxFF力场下,选用周期性边界条件、NPT系综、0.1 fs时间步长,运行50 000步,使模型达到指定温度和势能平衡,结果如图5所示。

图5 温度和势能的优化平衡

在对晶胞模型进行动力学弛豫达到平衡状态后,对模型在5×109s-1额定应变率下进行升温热解模拟。选用ReaxFF力场、NPT系综、0.1 fs时间步长,运行100 ps,分别加热到设定终温,得到的不同终温热解时应力变化关系如图6所示。由图6可以看出,在不同终温下,不同方向上的应力是上下波动的(正为拉力,负为压力),且波动范围较大,通过标准基线细化,发现终温在2 000 K之前应力波动范围几乎都在-70 GPa～70 GPa之间,终温在2 000 K以后,压力波动范围逐渐变宽,在区间之下的波动较多(小于-70 GPa),但通过整体比较,在允许一定误差的情况下,可认为在5×109s-1额定应变率下,应力区间主要集中在-70 GPa～70 GPa之间,可以进一步说明温度与应力属于独立主体。模拟结果中三轴的应力大小几乎吻合,表明在不同方向的瞬时应力大小是一致的,同时在热解过程中拉力和压力的交替作用,最终促使晶胞扩展约50%。

图6 不同终温热解时应力变化

3 结果与讨论

3.1 总体产物演变

对上述LAMMPS模拟结果进行提取和可视化图截取,开展产物和过程分析,参照张秀霞等[18]对碳产物的划分原则,将碳产物分为3大类:C0(H2,H2S等和自由基(·H,·H3N等))、C1-4(烃类小分子(CH4等)和含碳自由基(·CH3等))及C5-40(焦油)。由于晶胞碳含量甚微,C40+含量难以检测,此处不对其进行分析。不同终温下3类产物的质量分数如图7所示。

图7 不同终温下产物的质量分数

由图8[20]所示的义马原煤热解实验产物变化可知,C0和C1-4的质量分数随温度升高而增加,C5-40的质量分数随着温度上升先增后减,在650 ℃达到最大,而其他产物的质量分数逐渐降低,与本模拟结果有较高吻合度。

图8 义马原煤热解实验产物变化

由图7可知,随着热解终温上升,C0的质量分数持续增加,C1-4的质量分数逐渐减少,究其原因是在2 000 K之前碳链上的较弱共价键和活泼官能团受热分解成小分子气体产物和自由基;2 000 K以后发生了缩聚反应,焦油产率得以上升,降低了C1-4产率,表明应力和温度正向促进了反应进程,导致C0含量持续增加;说明应力对反应进程有促进作用,但无法改变整体反应趋势。而C5-40的质量分数呈现先增后减的趋势,在2 000 K终温下达到最大,随后逐步降低,究其原因是发生了缩聚反应,与范各庄烟煤升温模拟结果[19]基本相同。

2 000 K和2 500 K下ReaxFF MD热解状态构型如图9所示。由图9可以看出,随着模拟时间增加,碳链发生裂解生成小分子自由基(·H等),随后自由基攻击碳链的大分子稳定结构,进一步促进了煤裂解,导致热解产物C0和C1-4含量增加,但在75 ps～100 ps之间,2 500 K下的焦油含量相对低于2 000 K下的焦油含量,说明在2 000 K～2 500 K发生了缩聚反应,进一步验证了图7中的产物分布,亦表明热解过程主要受控于热解终温[18]。

图9 2 000 K和2 500 K下ReaxFF MD热解状态构型

图10所示为不同终温下产物数量和分子种类。由图10可知,不同终温下,随着模拟时间的增加,分子种类在26 ps～33 ps处出现峰值。33 ps之后随着模拟时间增加,产物数量持续增加,分子种类逐渐降低,符合褐煤热解实验和模拟的一般规律[20-23]。

图10 不同终温下产物数量和分子种类

在0～26 ps之间,分子数减少,种类数增加,且不同终温初始模拟时间下分子数和种类数几乎相同,与煤热解规律不完全相符,此时可能是应力作用结果。由图6可知,应力在模拟过程中,会不断出现拉力和压力交替作用,发生拉伸盒子运动,促进分子裂解,在前期温度较低阶段,应力起主导作用,即反应前期单一应力作用效果优于温度作用效果,应力促进了整体的反应进程。同时,对不同终温下的产物数量进行分析,发现在1 500 K终温后分子数超过初始应力作用的结果,可推测出,在1 500 K之前温度对煤热解的影响小于应力对煤热解的影响,1 500 K之后,温度占主导地位,表明1 500 K是义马原煤微观热解模拟的起始温度。

3.2 小分子气体和自由基演变

在煤热解过程中,C0和C1-4的小分子气体及自由基是主要产物。碳链上的C-H键、C-O键、桥键等先断裂,释放出有机小分子和大量活性自由基(·H等)[15],自由基会进一步加速烷烃的脱氢反应,生成小分子气体H2,CH4等[24]。

图11所示为不同终温下活性自由基·H和小分子气体H2数量的变化规律。由图11可以看出,随着温度升高,·H和H2的数量都持续上升。H2的产生可分为3种途径[18,25],一是温度升高促进碳链断裂,大分子脱氢生成小分子,导致大量H2产生;二是·H自由基与H2O发生碰撞,生成H2和·OH自由基;三是大量游离·H自由基相互吸引生成H2。在1 500 K之前,应力作用大于温度作用,导致H2的最终产生量小于初始应力作用的结果;而·H的数量在温度升高时会出现波动现象,主要原因是·H是非常活泼的自由基,不仅会自身结合生成H2,还会与·CH3结合生成CH4,与·OH结合生成H2O,与S离子反应生成硫化氢等;后续缩聚反应促进了焦油含量上升,降低了·H含量,但由于应力的作用且S离子的消耗殆尽又导致·H含量升高,而由图11也可以看出应力作用释放的·H的数量要大于缩聚反应消耗的·H的数量。

图11 不同终温下·H和H2的数量变化

图12 不同终温下其他小分子和自由基的变化

义马原煤热解实验[20]的气体体积分数和脱硫率如图13所示。由图13可知,H2含量随温度升高而上升,CO含量较稳定,CH4含量逐渐下降,基本与模拟结果的趋势相符。脱硫率逐渐升高,但是未达到60%,与模拟结果存在差异。究其原因是义马原煤中硫元素主要转化为H2S,COS和SO2,但H2和·H对COS,SO2的生成有抑制作用,对·SH与·H结合生成H2S有促进作用[30],并且模拟过程附加了一个5×109s-1的应变率,促使模拟过程·H含量增多,从而导致H2S的生成率高于实验过程中硫的转化率。

图13 热解实验的气体体积分数和脱硫率

4 结 论

1) 不同模拟终温下,允许一定误差下,5×109s-1应变率导致的单一应力几乎保持在±70 GPa之间。说明此应力对义马原煤热解的作用受温度的影响较弱。

2) 反应前期单一应力作用效果优于温度作用效果,应力促进了整体的反应进程,但没有改变反应的最终趋势。

3) 此单一应力作用下,焦油裂解释放·H的数量大于缩聚反应消耗·H的数量,同时促进H2S的生成。