具体实施方式

下面结合附图和实施例对本发明进一步说明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的技术特征可以相互结合。需要指出的是，除非另有指明，本申请使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。本发明公开使用的“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。

本发明提供一种液态硫在碳酸盐岩纳米孔隙中微观吸附特征参数的确定方法，所述微观吸附特征参数包括吸附分子量、超额吸附量、吸附热、吸附相密度、径向分布函数、扩散系数，所述确定方法包括以下步骤：

S1：建立目标岩样的分子模型。

在一个具体的实施例中，以某深层海相高含硫碳酸盐岩的目标岩样为例，首先利用场发射扫描电子显微镜、原子力显微镜等技术，对目标岩样进行分析，得到该深层海相高含硫碳酸盐岩孔隙结构的实际特征，扫描电镜结果如图1所示；根据得到的实际特征建立符合目标岩样的碳酸盐岩矿物分子模型。需要说明的是，根据目标岩样实际特征建立分子模型为现有技术，具体步骤在此不再赘述。

S2：挑选所述分子模型的吸附实验面，将所述分子模型沿所述吸附实验面进行切割并扩展。

S3：将扩展后的吸附实验面进行组合调整，获得符合周期性镜像原则的微观吸附骨架模型。

在一个具体的实施例中，步骤S2中挑选分子模型的吸附实验面，以及步骤S3中对扩展后的吸附实验面径向组合调整均采用现有技术中的挑选及组合调整方法。

本发明的微观吸附骨架模型的建立采用了周期性镜像原则，模拟中使用的分子构型具有周期性边界条件：在模拟盒子中的粒子受力运动的过程中，粒子从模拟盒子的一侧离开后会再从对侧进入模拟盒子中，以保证该模型体系中粒子数目的恒定。周期性镜像的存在，一方面能够解除粒子通过模拟盒子边界时出现的边界效应的影响，另一方面可以使用有限数目的粒子来模拟无限数目的粒子运动，进而提高模拟结果的精度。

S4：根据所述微观吸附骨架模型的大小设置截断半径。

在一个具体的实施例中，所述截断半径为所述微观吸附骨架模型最小边长的二分之一。

S5：使用COMPASSⅡ力场计算所述微观吸附骨架模型的总势能。

在一个具体的实施例中，所述总势能为键伸缩能、键角弯曲能、二面角扭曲能、键角面外弯曲能，耦合运动形式的键伸缩能、键角弯曲能、二面角扭曲能、键角面外弯曲能，以及静电相互作用力和范德华作用力之和；可选地，所述总势能通过下式进行计算：

式中：E为总势能，kcal/mol；为键伸缩能之和；为键角弯曲能之和；为键扭转能之和；为键角面外弯曲能；为键伸缩-键伸缩耦合作用能；为键伸缩-键角弯曲耦合作用能；为键角弯曲-键角弯曲耦合作用能；为键伸缩-二面角扭转耦合作用能；为键角弯曲-二面角扭转耦合作用能；为键角弯曲-二面角扭转-键角弯曲耦合作用能；为非键作用能；

K_b1、K_b2和K_b3为键伸缩力常数；b为键长，nm；b₀为键长初始值，nm；K_a1、K_a2和K_a3为键角弯曲力常数；θ为键角，°；θ₀为初始键角，°；K_t1、K_t2和K_t3为二面角扭转势常数；φ为二面角，°；K_x为离开平面的振动势常数；x为离开平面的振动角度，°；x₀为初始偏平面角度，°；K_bb为键伸缩-键伸缩耦合项常数；K_ba为键伸缩-键弯曲耦合项常数；b'为交叉项的键长，nm；b₀'为交叉项初始键长，nm；K_aa为键弯曲-键弯曲耦合项常数；θ'为交叉项键角，°；θ₀'为交叉项初始键角，°；K_bt1、K_bt2和K_bt3为键伸缩-键扭曲耦合项常数；K_at1、K_at2和K_at3为键弯曲-键扭曲耦合项常数；k为键弯曲-键扭曲-键弯曲耦合项常数；ε_i、ε_j为原子或分子i、j的势阱，kj/mol；ε₀为介电常数，8.854×10^-12F/M；σ_i、σ_j为原子或分子i、j的碰撞直径，nm；r_ij为为电荷i与电荷j之间的距离，nm；q_i为电荷受体的电荷值，C；q_j为电荷给体的电荷值，C。

S6：根据需要模拟计算的地层压力及温度区间，设置模拟计算时的温度和压力间隔梯度值，确定模拟计算的步长。

S7：计算不同温度压力下液态硫单质分子的逸度。

逸度是化学热力学表示气体实际压强的物理量，后续模拟均用所述逸度代替压力。在一个具体的实施例中，所述逸度通过下式进行计算：

f^l＝p^l,vexp(V^l(p-p^l,v)/RT) (2)

其中：

p^l,v＝p_cexp(-9.741(1-T/T_c)+5.127(1-T/T_c)^1.5-0.205(1-T/T_c)³+9.502(1-T/T_c)⁷)T_c/T (3)

当温度小于432.6K时，其中：

V^l＝115.708+0.06735T (4)

当温度大于432.6K且小于等于577.96K时，其中：

V^l＝171.075-0.1442T+0.000193T² (5)

式中：f^l为逸度，Pa；p^l,v为饱和蒸汽压，Pa；V^l为体积，cm³/mol；p为压力，Pa；R为理想气体常数，8.314472×10⁶cm³·Pa·mol^-1·K^-1；T为温度，K；p_c为临界压力，10.4209MPa；T_c为临界温度，1115.0292K。

S8：进行吸附等温线的巨正则系综蒙特卡洛模拟，得到某一设定温度下压力区间内的吸附热、吸附分子量、吸附等温线、以及分子相对集中度。

蒙特卡洛方法(MC)是基于统计力学并结合系综来计算其平均值的一种随机模拟方法，用Metropolis抽样方法对产生的随机数进行抽样从而达到简化计算的目的，并获取物质的微观结构及其性质。Metropolis抽样方法主要是在相空间中构建Markov链从而能够实现重要性抽样，Markov链的构建须满足如下几个条件：第一，随机抽样所产生的点必须属于体系中各点所组成的集合中的一点；第二，随机抽样所产生的点只能与该点的前一个点相关，而与从何处开始抽样无关。Metropolis蒙特卡洛方法的原理如图2所示。

蒙特卡洛模拟在本发明中的主要实施步骤为：

(1)搭建矿物骨架模型，保证其结构合理即可；

(2)选择系综；用不同系综计算会导致不同的计算结果，目前常用的系综包括巨正则、Gibbs、微正则、等温、等压及恒应力—等温等几种系综；

(3)依据模型中的随机变量分布模式和系综规则，随机产生一个单质硫分子在孔隙结构中分布的分子构型；

(4)设定计算的实现方式，确定需要获得的物理量数据(如吸附量、等量吸附热和构型能量等)；

(5)对分子构型中的粒子进行无规则的改变，比如产生、消失、平移和旋转；

(6)通过对每一次新的构型的能量进行计算，对比其变化前后的能量变化：

若新的构型能量比原构型能量低，则接收新构型，并将这个新的构型重复进行下一次迭代计算；

若新构型的能量比原构型的能量高，则计算此构型的玻尔兹曼因子，并产生随机数：

若产生的随机数小于所计算的玻尔兹曼因子时，则接收当前新构型；

若随机数大于玻尔兹曼因子，则舍弃构型进行下一次重复计算；

(7)直到最后搜索出比所给能量条件低的所有构型时结束计算；

(8)对平衡模拟结果进行统计分析，需要注意的是这个结果一定要与初始条件无关，否则无效。

在本发明中，步骤(2)选择系综时选用巨正则系综。巨正则系综的蒙特卡洛方法包含了三种粒子运动形式，分别为粒子的插入、删除及移动。

(1)粒子插入

在模拟盒子中的任意位置随机插入粒子源中的一个粒子，该粒子被随机挑选，粒子插入盒子中被接受的概率是粒子源中的化学势的函数：

式中：V'为模拟盒子的体积，m³；Λ为粒子的德布罗意波长，m；k_B为玻尔兹曼常数；μ为粒子源的化学势，J/mol；U为体系的总能量，J；N为粒子序列，个。

(2)粒子删除

在插入粒子后的模拟盒子中，随机删除盒子中一个粒子，删除该粒子被接受的概率也为粒子源中的化学势的函数：

(3)粒子的移动、转动

随机选择模拟盒子中已插入的粒子，将移动到新的位置或转动该粒子，该粒子移动或转动被接受的概率为：

P(s→s′)＝min(1,exp{-β[U(s′)-U(s)]}) (13)

式中﹐s和s'为粒子移动、旋转前后在体系中所处的状态。

在一个具体的实施例中，在温度为403.15K、压力为45MPa的条件下，巨正则系综蒙特卡洛模拟结果如图3所示。

S9：根据所述步长更改设定温度，重复步骤S8，直至获得所有温度区间下各压力区间内的吸附热、吸附分子量、吸附等温线、以及分子相对集中度。

在一个具体的实施例中，目标岩样所有温度区间下各压力区间内的吸附热如图4所示。

S10：根据步骤S8和步骤S9获得的参数，计算各温度下的吸附相密度、超额吸附量。

在一个具体的实施例中，所述吸附相密度通过下式进行计算：

式中：ρ_b为吸附相密度，g/cm³；C为分子相对集中度，％；N_t为吸附分子量，moleculars/uc；M为液态硫单质的相对分子质量，1；V为模拟的晶胞体积，cm³；N_A为阿伏加德罗常数，无量纲；

所述超额吸附量通过下式进行计算：

式中：：n_ex为超额吸附量，mmol/g；V_b为吸附相的体积，cm³；M_a为吸附质分子的质量，g。

需要说明的是，步骤S8和步骤S9获得的吸附分子量为绝对吸附量，根据所述吸附分子量采用式(7)计算所述超额吸附量时，需对所述吸附分子量进行单位换算，具体可通过下式进行单位换算：

S11：根据各温度下的吸附等温线，获得各温度下达到饱和吸附状态时对应的饱和压力值。

S12：分别进行各饱和压力值和对应的温度值条件下的定压吸附蒙特卡洛模拟，并筛选出体系总能量最低的分子构型。

S13：根据体系能量最小化对所述分子构型进行体积优化，得到稳定吸附构型。

需要说明的是，步骤S11-S13中，获得饱和压力值、定压吸附蒙特卡洛模拟、筛选总能量最低分子构型、体系能量最小化体积优化均为现有技术，具体方法在此不再赘述。在一个具体的实施例中，所述稳定吸附构型如图5所示。

S14：对所述稳定吸附构型进行分子动力学模拟，获得径向分布函数和均方位移。

分子动力学(MD)模拟方法是一种用来计算经典多体体系的平衡以及传递性质的模拟方法，体系中粒子的运动满足牛顿运动定律。

分子动力学模拟方法的基本原理：根据体系中粒子的位置，求取体系的势能，根据体系的势能得到作用在每个粒子上的力，求解牛顿运动方程得出体系任意时刻的相关参数，再统计平均得到体系的热力学参数。MD模拟方法的原理如图6所示。

MD模拟在本发明中的主要实施步骤为：

(1)在矿物—硫最低能量构型的基础上，给定其初始条件：初始位置r(0)，v(0)；

(2)由硫分子的位置和势能函数求出系统的势能，据牛顿第二定律，计算硫分子在t时刻的所受到的作用力与加速度；

(3)计算t+Δt时刻硫分子的位置，并求取速度及系统势能；

(4)重复上述步骤，直到给定时间步长内所有的硫分子信息记录完毕，记录下各个时间下硫分子的信息，生成运动轨迹；

(5)基于硫分子的运动轨迹，可以对其径向分布函数、均方位移、自由能等进行计算，来分析单质硫分子在矿物表面的吸附性质。

求解粒子运动方程有Verlet方法、Leap frog方法、Beeman方法、Velocity-Verlet方法、Predictor-corrector方法等。Verlet算法是求解牛顿运动方程中最基本也是最简单的算法，该算法核心如下：

将体系中某时刻(如t+Δt时刻)各粒子的位置进行泰勒展开可得：

同样，对在t-Δt时刻粒子的位置进行泰勒展开可得：

两式相加可得：

忽略高阶项可以得到：

通过上式可知，在已知初始粒子的位置和下一个时刻粒子的位置及加速度情况下，可以预测粒子在新的时刻位置。两式相减后忽略高阶项，可以得到粒子的运动速度：

对于其他的求解牛顿运动方程的算法在预测粒子的位置和计算粒子的运动速度方面的精度或者执行效率更高。不同的算法都有着不同的优点，如Verlet方法与Leap frog方法相比较，Verlet方法能够直接计算体系的总能量，而Leap frog方法则不能，但Leapfrog方法在计算粒子的速度方面精度却是更高。实际使用时，可根据需要选择不同的算法。

径向分布函数(RDF)是一种描述体系结构粒子聚集有序性的普遍适用方法，也是一种研究体系热力学性质的重要手法。

RDF是指在距一中心原子的距离为r到r+dr区域范围之内找到另一目标粒子的概率和这另一目标粒子在这个区域范围内出现的概率之比，是表征液体微观结构特征的物理量。体系的径向分布函数等于整个体系中所有原子的径向分布函数代数平均值。定义径向分布函数g(r)为：

式中：dN表示距中心原子的距离为r到dr区域范围内找到的目标粒子的数目，ρ'为目标粒子的密度。

在一个具体的实施例中，所述径向分布函数如图7所示，所述均方位移如图8所示。

S15：绘制所述均方位移与时间的关系曲线，计算所述关系曲线的斜率，根据所述斜率计算得到扩散系数，所述扩散系数通过下式进行计算：

D＝a/6 (8)

式中：D为扩散系数，无量纲；a为斜率，无量纲。

在一个具体的实施例中，步骤S10中，还包括根据步骤S8和步骤S9获得的参数，计算吉布斯自由能，并根据所述吉布斯自由能计算吸附过程中的熵变的步骤；

所述吉布斯自由能通过下式进行计算：

式中：ΔG为吉布斯自由能变化量，KJ/mol；ΔH为吸附热，KJ/mol；T为吸附过程的绝对温度，K；p为吸附过程的压力，MPa。

所述熵变通过下式进行计算：

式中：ΔS为吸附过程的熵变，KJ·(mol·K)^-1。

在一个具体的实施例中，步骤S10和步骤S11之间，还包括以下步骤：

S11'：根据所述吉布斯自由能或所述吸附热判断吸附类型；根据所述吸附热判断吸附作用力类型；具体的：

当所述吉布斯自由能变化量在-20～0KJ/mol，或所述吸附热的绝对值小于40KJ/mol时，吸附类型为物理吸附；当所述吉布斯自由能变化量在-400～-80KJ/mol，或所述吸附热的绝对值在50～200KJ/mol时，吸附类型为化学吸附；

当所述吸附热在4～10KJ/mol时，吸附作用力中有范德华力；当所述吸附热在2～40KJ/mol时，吸附作用力中有氢键；当所述吸附热为40KJ/mol时，吸附作用力中有配位基交换；当所述吸附热在2～29KJ/mol时，吸附作用力中有偶极间作用力；当所述吸附热大于60KJ/mol时，吸附作用力为化学键。

上述吉布斯自由能、熵变、吸附类型、吸附作用力类型虽然不是液态硫在碳酸盐岩纳米孔隙中必要的微观吸附特征参数，但是其对分析液态硫的吸附机理具有很大的帮助。

本发明通过上述方法能够定量获得液态硫在碳酸盐岩矿物中的吸附特征参数：吸附分子量、超额吸附量、吸附热、吸附相密度、径向分布函数、扩散系数，与现有技术相比，具有显著的进步。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。