NaX分子筛吸附天然气中酸性气的分子模拟

王 翀，鲁家荣，闫 昊，刘熠斌，陈小博

(1.中石化石油工程设计有限公司，山东 东营 257026；2.中国石油大学(华东)重质油国家重点实验室)

天然气的主要成分为CH4,一般还存在部分的CO2、H2S等酸性气。酸性气的存在不仅会在水蒸气的作用下腐蚀设备和管路,而且会影响天然气的品质及后续加工过程。因此,天然气中酸性气的脱除是天然气工业的重要过程。目前酸性气主要脱除技术为湿式化学脱除,具体做法是采用碱性液体作为吸收剂,使其与天然气中的H2S和CO2发生化学反应以达到脱除酸性气的目的[1]。尽管湿式化学法脱除酸性气的效率较高,但能耗高,且存在溶剂废液污染。近年来节能环保型酸性气脱除技术有了长足的发展,如出现了分子筛变压吸附技术、膜分离技术,其中分子筛变压吸附技术具有脱除率高、能耗低、无腐蚀、操作简单等优点。

研究天然气中不同组分在分子筛中的吸附扩散行为,对于优化分子筛设计、提升酸性气脱除能力具有重要的指导意义。分子模拟是通过理论计算模拟分子结构与行为的方法,可用于多孔材料内的吸附与扩散研究。Maurin等[2]以不同离子交换的X分子筛为研究对象,采用巨正则蒙特卡罗(GCMC)方法结合根据离子本征特性导出的力场模型进行了模型分子吸附行为的模拟,计算的等温线和微分等量热与试验结果具有良好的一致性。Demir等[3]采用GCMC方法对CH4、CO2体系在类沸石金属有机骨架(ZMOFs)中的吸附行为进行了研究,发现在二元吸附体系下CH4和CO2在ZMOFs中的吸附位点相同,ZMOFs对CO2的吸附优于对CH4的吸附。Sung Chunyi等[4]通过分子模拟方法研究了多价阳离子交换的Y分子筛对克劳斯反应尾气中微量H2S的选择性吸附,并通过密度泛函理论(DFT)方法计算了不同材料对H2S的吸附能,提出使用GaY分子筛进行H2S的吸附是有效的脱硫策略。丁雪等[5]采用GCMC方法研究了干气中不同组分在ZSM-5分子筛中的吸附行为,并计算了吸附过程的焓变、熵变等热力学数据。党宇等[6]则将GCMC方法与分子动力学方法结合,研究了压力对噻吩、吡咯、呋喃3种杂环化合物在HY分子筛中的吸附与扩散性能,获得了概率密度、径向分布函数、吸附能等。

作为一种具有广泛应用前景的X分子筛,NaX分子筛被广泛应用于天然气中酸性气的脱除,但是对天然气中不同组分在NaX分子筛中竞争吸附扩散行为的分子模拟,尤其是针对同时含有CO2、H2S的体系的研究则鲜有报道。本课题采用GCMC与分子动力学结合的方法,研究3个温度(273,293,313 K)下天然气中CH4,H2S,CO2在NaX分子筛中的吸附与扩散行为,为理解天然气吸附脱硫、脱碳的微观行为和吸附材料的理论设计提供指导。

1 计算模型与方法

1.1 分子筛吸附剂模型

X分子筛是一种硅铝比为1.0～1.5的多孔硅酸盐材料,其拓扑结构为FAU型,属于立方晶系,空间群为Fm-3m,晶胞参数a=b=c=2.474 nm,α=β=γ=90°,笼口孔径为0.74 nm[7]。13X分子筛的晶胞由β笼和六方柱笼构成,β笼与β笼之间通过六方柱笼连接,形成窗口直径为0.74 nm的十二元环和对应的超笼结构。

本研究采用Materials Studio 8.0软件平台,分子筛的空间群和晶胞参数如上文所述。为简化模型,在13X分子筛的建模中采用硅铝比为1的结构。由于分子筛骨架为负电荷,因此引入Na+进行电荷平衡,平衡后的晶胞组成为:Na96Al96Si96O384。模型建立后,根据文献[8-9]分别为各种原子进行电荷分配:Na(+0.70),Al(+1.75),Si(+2.35),O(-1.20)。结构优化后的NaX分子筛结构模型如图1所示。

图1 NaX分子筛的结构模型

采用Materials Studio 8.0软件中分子动力学模块对NaX分子筛结构进行优化,并获得最低能量态。为了验证晶体模型的合理性,使用Reflex模块模拟了该结构的X射线衍射图谱(图2)。将模拟结果与数据库[10]中相应结构的标准图谱(图2中FAU)对比,发现两者具有相同的特征峰,证明结构优化后的晶体类型为FAU结构。

图2 优化结构的XRD模拟结果与数据库中FAU结构特征峰的对比

NaX分子筛的结构由带有负电性的硅氧铝骨架结构和阳离子构成,其中阳离子的类型以及分布对于其吸附性能具有重要影响[11]。Frising等[12]详细总结了不同文献中合成的NaX分子筛的阳离子分布情况。由于各种文献对Na+位点的定义繁杂并且Na+在分子筛中的分布较为复杂,在总结各种文献后给出Na+一般分布情况:Na+在骨架中主要分布在如图3所示的3种位点,分别是Ⅰ(方钠石笼,另有Ⅰ′位点与该位点十分接近,统称为Ⅰ位点)、Ⅱ(方钠石笼和超笼之间的六边形窗口中心)、Ⅲ(超笼中靠近中心四边形窗口,另有超笼中靠近周边四边形窗口的位点Ⅲ′与之十分接近,文献中对Ⅲ和Ⅲ′的界定较为模糊,故统称为Ⅲ位点)。不同文献所报道的Na+分布情况如表1所示。将结构优化后的NaX分子式模型与表1的数据相对比,发现优化模型具有上述4种阳离子分布位点,位点Ⅰ,Ⅱ,Ⅲ的个数分别为24,30,31个,与文献中的结果近似,证明模型具有合理性。

表1 不同文献所报道的NaX分子筛骨架外阳离子分布情况

图3 NaX分子筛的阳离子分布情况

1.2 吸附质分子模型

选取天然气中的CH4,CO2,H2S作为吸附质分子,3种分子均采用Materials Studio软件建模并采用Dmol 3量子力学模块优化,获得其ESP电荷(见表2)。3种分子的分子尺寸、四极矩等数据如表3所示。由表3可知,3种分子的动力学直径均小于NaX分子筛的孔口直径,因此3种分子均能扩散进入分子筛孔道进行吸附。模拟中将调节不同分子的摩尔比,探究其在分子筛模型上的吸附行为。

表2 天然气不同组分分子中各原子的计算ESP电荷

表3 不同分子的性质

1.3 模拟参数设置

NaX分子筛对吸附质分子的吸附模拟采用Sorption模块,抽样方法选用Metropolis方法,精度选择Ultra-fine。客体分子和分子筛的相互作用包括静电作用和范德华作用,其中范德华作用采用Lennard-Jones势函数描述,其数学表达式为:

(1)

式中:ULJ表示客体分子和吸附剂之间相互作用能,eV;i和j表示不同原子;Rij表示原子间距,nm;Dij和(R0)ij为Lennard-Jones参数;qi和qj表示原子所带电荷,e。模拟中采用COMPASSⅡ力场,静电相互作用采用Ewald方法处理,非键相互作用的加和采用Atom based算法,范德华作用截断距离设置为1.25 nm,正好小于晶胞边长(2.502 8 nm)的一半。计算平衡步数为1×106步,生产步数为1×107步。

吸附质分子的Dmol 3结构优化参数:精度选择Fine,选用广义梯度近似泛函(GGA)中的PBE近似,DFT-D色散校正选用TS方法。核处理选用全部电子,基组选用DNF数值基组。

分子在分子筛中的扩散采用Forcite模块,根据牛顿力学的原理进行分子动力学计算。模拟中选用正则系综(NVT),采用Velocity-Verlet算法进行积分,模拟步长1 fs,模拟总时长200 ps,其中前50 ps用于平衡体系。扩散系数Ds采用Einstein方程计算:

(2)

式中:Ds为扩散系数,m2/s;N表示单位晶胞内的分子数;r(t)是吸附质分子在t时刻的质心坐标;r(0)是吸附质分子质心的初始位置坐标;t表示时间,ps;平方项是扩散分子均方根位移(MSD)的系综平均。

2 结果与讨论

2.1 吸附等温线和吸附热

为了分析3种组分在NaX上的吸附强度,采用GCMC方法模拟计算了工业天然气变压吸附脱除酸性气条件下,CH4,CO2,H2S在不同温度(273,293,313 K)下的吸附等温线数据(见图4中散点),并采用Freundlich吸附模型对3种分子在3个温度下的吸附等温线数据进行拟合,拟合曲线见图4中线条。

图4 纯组分CH4,CO2,H2S在NaX分子筛中的吸附等温线

由图4可知:3种物质的吸附等温线均为Ⅰ型等温线,表示3种分子在相应条件下在NaX分子筛上的吸附形式均为单分子层吸附;随着平衡压力的增加,3种组分的吸附量均表现出先增加后趋于饱和的规律;在压力较低(0～2 MPa)时,CH4的吸附量增加相对缓慢,且在7 MPa左右达到吸附饱和;相同压力范围内H2S和CO2的吸附量迅速增加,在2 MPa附近即达到饱和。

Freundlich吸附模型适用此模拟的压力范围,其吸附量计算式如式(3)所示。

(3)

式中:N为吸附量,mmol/g;Nm为单层饱和吸附量,mmol/g;K为吸附平衡常数,MPa-1,数值上等于吸附速率常数与脱附速率常数之比,表征分子筛对某物质的吸附强弱;1/n为指数,表示浓度对吸附强度的影响,其值一般小于1;P为模拟条件下的压力,MPa。不同温度下CH4,CO2,H2S分子在NaX分子筛中吸附的Freundlich模型拟合参数分别如表4～表6所示。

表4 不同温度下CH4分子在NaX分子筛中吸附的Freundlich模型拟合参数

表5 293 K下3种分子在NaX分子筛孔道中的自扩散系数

表5 不同温度下CO2分子在NaX分子筛中吸附的Freundlich模型拟合参数

表6 不同温度下H2S分子在NaX分子筛中吸附的Freundlich模型拟合参数

由表4～表6可以看出:Freundlich吸附模型可以很好地描述3种分子在NaX分子筛上的吸附等温线,拟合方程的决定系数R2均大于0.99。表4～表6中的Nm,K,1/n分别为根据Freundlich吸附模型拟合得到的数据,拟合的Nm结果与采用GCMC方法模拟得到的单层饱和吸附量十分接近,表明Freundlich吸附模型可以较为合理地描述3种分子在分子筛中的单层吸附情况。

采用GCMC方法模拟得到在温度273,293,313 K下,CH4的饱和吸附量分别为6.79,6.43,6.05 mmol/g,CO2的饱和吸附量分别为9.13,8.74,8.42 mmol/g,H2S的饱和吸附量分别为9.84,9.45,9.23 mmol/g,可见NaX分子筛对H2S和CO2的吸附量高于对CH4的吸附量。另外,3种分子的饱和吸附量均随着温度的降低而升高,表明高压和低温有利于CH4,CO2,H2S在NaX分子筛上的吸附。在拟合压力范围内,3种分子在NaX分子筛中吸附的平衡常数由大到小的顺序为H2S>CO2>CH4,说明NaX对3种分子的吸附强度由大到小的顺序为H2S>CO2>CH4,表明NaX分子筛对于天然气中的酸性气分子具有较强吸附能力。3种分子在该温度下的微分吸附热也具有上述规律。由于CH4分子的偶极矩和四极矩为零,CH4与NaX分子筛孔道的作用力主要为范德华力,而CO2和H2S分子除了上述作用外,还与分子筛骨架外Na+产生较强的静电相互作用。H2S由于具有较大的偶极矩和较强的四极矩以及极化率,其与骨架外阳离子静电相互作用的强度高于CO2,因此出现低压吸附平衡常数和饱和吸附量由大到小的顺序H2S>CO2>CH4。Maghsoudi等[21]通过试验的方法研究了天然气中CH4,CO2,H2S在CHA分子筛中吸附和扩散行为,结果表明CHA分子筛在CH4的存在下可以选择性地分离H2S和CO2,其测定的吸附等温线对饱和吸附量的结论与本研究结论一致。

在获得3种分子在NaX中的吸附等温线的基础上,采用GCMC模拟方法得到了3种分子等量吸附热随着压力变化的图线。等量吸附热ΔQ的计算参见文献[6]。

由于温度对等量吸附热的影响很小,故这里仅给出温度273 K下纯物质CH4,H2S,CO2的吸附热曲线(如图5所示)。由图5可知,3种分子的吸附热由大到小的顺序为H2S>CO2>CH4,表明分子筛对CO2和H2S的吸附作用较强,与吸附活性位作用后产生了较为明显的热效应,与吸附等温线的模拟结果相符合;另外,3种分子等量吸附热均随着吸附压力的增加而降低,最终达到稳定,这与很多小分子的吸附呈相似的规律[6]。

图5 273 K下纯组分在NaX分子筛中的吸附等量热曲线

2.2 吸附位分析

根据物质在不同温度下的吸附相互作用能量分布曲线可以得到吸附质分子在分子筛中的吸附位点分布,获得吸附质的最佳吸附位置,探究温度对吸附位的影响[13]。图6为CH4,H2S,CO2分子在273,293,313 K下与NaX分子筛的相互作用曲线。由图6可知,3个温度下3种分子的吸附峰的位置基本不发生改变,说明一定温度范围内温度对吸附位点的影响较小。CH4吸附能量分布曲线出峰位置在-32 kJ/mol和-20 kJ/mol,且在-20 kJ/mol处为强峰;CO2吸附能量分布曲线出峰位置在-42 kJ/mol和-20 kJ/mol两处,在-42 kJ/mol处为强峰;H2S吸附能量分布曲线只有1个出峰位置,为-45 kJ/mol。可见3种分子在分子筛内主要有3种吸附位点,其对应的吸附能量分别为-20,-32,-42 kJ/mol左右,分别对应NaX分子筛中孔道中的活性位点Ⅱ,Ⅰ,Ⅲ。CH4主要吸附在Ⅱ位,H2S和CO2主要吸附在Ⅲ位。

图6 不同分子与分子筛间的相互作用曲线

以温度293 K为例,根据模拟计算得到了3种分子在该温度下的吸附概率密度势能场分布,如图7所示(红色代表分子在该处出现的概率高,相互作用强;绿色代表出现的概率低,相互作用弱)。由图7可以看出,CH4分子主要分布在方钠石β笼中以及方钠石笼和超笼之间的交界;CO2分子和H2S分子则主要分布在超笼中,方钠石笼中的分布量较少。从吸附质与分子筛相互作用上分析,可以看出CH4分子与NaX骨架的相互作用力较小,且相互作用较强的位置位于β笼中;CO2和H2S与骨架的相互作用更强,且主要作用位点在超笼中的Na附近。CO2在NaX中的分布情况与作用位点与Cui Yongkang等[14]的研究结果相一致。

图7 293 K下不同分子的吸附密度势能场分布

2.3 混合组分的竞争吸附

为了探究实际混合天然气体系在分子筛中的吸附行为,首先构建了3种不同CO2含量的混合气体,其中H2S的摩尔分数固定为3%,CO2的摩尔分数分别为10%,15%,20%。293 K下的吸附等温线如图8所示。由图8可以看出,尽管混合气体中CH4的含量占主导,但是NaX分子筛对CH4的吸附量远低于对CO2和H2S的吸附量,体现出NaX分子筛具有较高的酸性气脱除选择性。随着CO2含量的增加,CO2的饱和吸附量逐渐增加并且超过H2S的饱和吸附量。产生这种现象的原因可能是两者共同竞争骨架的Ⅲ号吸附位,随着CO2分子数增加,其占据了更多的吸附位,并获得了更高的吸附选择性。随着CO2摩尔分数从10%增加到15%,CO2的饱和吸附量增加了18.2%;而当CO2摩尔分数从15%增加到20%后,其饱和吸附量仅增加了7.9%。说明随着CO2分子数的增加,吸附选择性从CO2浓度主导转变为由H2S、CO2分子与吸附剂的相互作用主导,因此CO2吸附量增加缓慢。

图8 293 K下CO2含量不同的混合气体中各组分分子的吸附等温线

将CO2的摩尔分数固定为3%,调节CH4和H2S的比例,控制H2S摩尔分数分别为3%,6%,9%,研究293 K下H2S含量对分子吸附行为的影响,获得的等温线如图9所示。由图9可以看出,分子筛对CH4分子的吸附量仍远小于对H2S分子和CO2分子的吸附量。随着H2S含量的增加,H2S和CO2的饱和吸附量均变化较小,且前者的饱和吸附量远高于后者,体现出NaX分子筛对H2S具有较高的吸附选择性。

图9 293 K下H2S含量不同的混合气体中各组分分子的吸附等温线

2.4 扩散过程分析

客体分子在吸附剂中的扩散行为是吸附过程中的重要性质。由于CH4,CO2,H2S分子的动力学直径均小于NaX分子筛的孔道直径,因此3种分子的扩散性质主要受到分子的偶极矩以及极性的影响。扩散系数是表征扩散能力和计算传质速率的重要参数。本研究采用分子动力学的方法计算了293 K下CH4,CO2,H2S在NaX分子筛中的扩散轨迹,并进一步得到了3种分子的均方根位移(MSD)随时间的变化曲线,如图10所示。

图10 293 K下3种分子在NaX分子筛孔道中扩散的均方根位移

根据式(2)Einstein方程计算CH4,CO2,H2S在分子筛中的扩散系数,结果如表5所示。由表5可知,3种分子在分子筛中的自扩散系数由大到小的顺序为CH4>H2S>CO2。由吸附性质计算可知,由于H2S和CO2与NaX分子筛表面具有较强的相互作用,因此较难在孔道中扩散。CH4与NaX分子筛表面仅有范德华相互作用,因此在孔道中的扩散阻力最小。扩散系数的结论与吸附性质的计算结果具有较好的符合度。

为了进一步分析3种分子与分子筛表面Na+的详细作用情况,对优化后的结构采用Forcite模块分析工具进行分析。径向分布函数分析法是分析体系微观结构和相互作用的常用方法。径向分布函数g(r)的数学定义是以某一分子或原子为中心,在半径为r到r+dr的球壳中找到其他粒子的概率[15],其表达式为:

(4)

式中:δr为球壳厚度,nm;r为客体原子与中心原子之间的距离,nm;n(r)为球壳中的客体粒子数。对293 K下CH4,CO2,H2S分子与分子筛中Na+的径向分布函数进行了分析,结果如图11所示。由图11可以看出:H2S与CO2中心原子与Na+的径向分布图在r=0.4 nm处均出现了较强的峰,说明两种分子的中心均与分子筛骨架外的Na+形成了相互作用,其中H2S的主峰的对应半径是最小的,这是因为H2S分子结构为折线形,S原子容易和孔道中的Na+接触并相互作用,空间位阻较小;CH4由于结构为正四面体,C原子空间位阻较大,不易与孔道中的活性位点相互作用,因此在作用范围内几乎不存在尖峰。除此之外,由3种分子的端位原子径向分布函数可以看出,H2S和CH4在r=1 nm的范围内几乎不存在峰,说明两者的端位原子与活性位点没有较强的相互作用。从CO2的径向分布函数可以看出,中心原子C和端位原子O都与骨架中的Na+产生相互作用,相对于C—Na的径向分布,O—Na的径向分布具有更近的主峰半径和更高的概率值,这说明相对于C原子,CO2中的O原子与分子筛中的活性位点有更强的相互作用,这也解释了CO2扩散系数最小的原因。

—C; —H

3 结 论

(1)在温度273,293,313 K下,CH4,CO2,H2S分子在NaX分子筛中的吸附量随着压力的增加而迅速增大,并趋于饱和吸附量。在一定温度范围内,温度的升高降低了3种分子的饱和吸附量。相同温度下吸附作用强度和对应饱和吸附量由大到小的顺序为H2S>CO2>CH4。

(2)在温度293 K下,混合气体中CO2含量提高,CO2的吸附量增加;随着混合气体中H2S含量提高,3种分子的饱和吸附量基本保持稳定,但H2S吸附量远高于CO2。NaX分子筛对混合气体中3种分子的吸附选择性由大到小的顺序为H2S>CO2>CH4,其表现出较强的选择性吸附脱硫性能。

(3)CH4,CO2,H2S分子在NaX分子筛内的扩散系数由大到小的顺序为CH4>H2S>CO2。径向分布函数分析结果表明,CO2的中心原子和端位原子均与孔道活性位产生较强相互作用,解释了CO2的强扩散阻力现象。