煤炭开采机械化程度不断提高，粉尘量也明显增多，在煤炭掘进和开采过程中都会产生大量的煤尘、矽尘等混合性粉尘，会给作业人员的生命健康带来巨大的危害。因此对存在粉尘的长期作业场所，粉尘治理成了越来越重视的问题。

随着分子模拟技术的发展，经常使用分子动力学模拟研究分子之间的微观相互作用^[1-2]。高正阳等^[3]研究了4种不同煤阶煤分子表面水分子吸附的微观机理；XU等^[4]和HUANG等^[5]分别研制了针对建筑工地粉尘和岩石粉尘的新型环保抑尘剂，并通过试验得出新型抑尘剂具有突出的抑尘效果；CRAWFORD等^[6]研究了表面活性剂对不同煤表面特性的影响，发现添加表面活性剂可以明显降低煤的接触角，提高煤的润湿效果；ZHANG等^[7]通过分子动力学模拟研究了表面活性剂对褐煤吸附行为和润湿性的影响，发现表面活性剂可以改变水与褐煤之间的相互作用力；陈璐瑶^[8]研究了表面活性剂结构对SiO₂粉尘润湿性，通过分子动力学模拟，得出SiO₂表面基团性质是影响SiO₂润湿性的根本原因；孟筠青等^[9]探究不同种类表面活性剂及其复配对SiO₂润湿性的影响，同时结合沉降试验，找出SiO₂粉尘润湿效果最好的单体表面活性剂及复配溶液；邱玄等^[10]通过分子动力学模拟了 Pluronic F-68 在萤石和石英表面的吸附行为，结果显示 Pluronic F-68在石英上的吸附能力更强。

上述研究表明,分子模拟已成为研究润湿和吸附性能的重要方法^[11-12]。但研究都是针对单一性粉尘研究润湿性。为此，借助于Materials studio分子模拟软件构建模型，采用分子动力学模拟和量子化学的方法，以表面活性剂CTAB为例，研究褐煤煤尘、矽尘的润湿机理的差异性。

1.   模拟方法

1.1   模拟体系构建

采用Materials Studio 分子模拟软件，构建了水−表面活性剂−褐煤煤尘/矽尘模型。初始模型如图1所示。

图  1  初始模型

Figure  1.  Initial model

下载: 全尺寸图片 幻灯片

选用Wender经典褐煤模型，其煤分子结构如图1（a）所示；矽尘模型为SiO₂超晶胞（图1（b）），添加1.5 nm的真空层以消除周期性的影响；表面活性剂选用十六烷基三甲基溴化铵（CTAB）（图1（c））。利用MS中Amorphous Cell 模块构建各个盒子，盒子详细参数见表1。

表  1  盒子或超晶胞参数

Table  1.  Box or supercell parameters

名称	长/nm	宽/nm	高/nm	分子数量	α、β、γ/（°）	密度/ （g·cm−3）
水盒子	2.946 0	3.241 2	1.566 4	500	90	1.000
CTAB盒子	2.946 0	3.241 2	0.239 7	5	90	1.322
煤分子盒子	2.946 0	3.241 2	1.342 6	8	90	1.350
SiO2超晶胞	2.946 0	3.241 2	2.288 3	—	90	2.200
注：α、β、γ为三维空间中与x、y、z轴的夹角。

下载: 导出CSV 

| 显示表格

通过build layer构成水−表面活性剂−褐煤/矽尘（SiO₂）2个模型，并在体系上方都添加了0.5 nm真空层。模拟体系的构建如图2所示，模拟平衡后水的密度场分布图如图3所示。

图  2  模拟体系构建

Figure  2.  Simulation system construction

下载: 全尺寸图片 幻灯片

图  3  模拟平衡后水的密度场分布图

Figure  3.  Density field distribution of water after simulated equilibrium

下载: 全尺寸图片 幻灯片

1.2   计算模拟方法

1）分子动力学模拟。运用Forcite模块对褐煤分子、水分子、CTAB分子和SiO₂超晶胞进行结构优化。首先采用Geometry Optimization任务中通过Smart算法进行局部结构优化。然后，利用Amorphous Cell 模块构建500个水分子盒子、5个表面活性剂盒子、8个褐煤分子盒子，通过build layer构成模型，在模型上方添加0.5 nm真空层。然后，为了得到模型全局最小和稳定的能量，在Anneal任务中采用NVT系综（温度为298 K）进行退火动力学模拟。最后，在Dynamics任务中采用NVT系综进行1 ns的分子动力学模拟，其中力场采用COMPASSⅡ，电荷采用Forcefield assigned。

2）量子化学模拟。量子力学计算在DMol3模块中完成^[13]，在广义梯度近似(GGA)^[14]下，利用BLYP泛函模拟电子交换和相关相互作用^[15]，模拟精度为Fine，其中SCF收敛精度为0.000 001。选择所有由D极化函数（DNP）增强的电子双数值原子轨道作为基础集，并将轨道截止质量定为良好。smearing的值选择为0.005，并且在计算过程中不考虑基集叠加误差（BSSE）^[16]，因为在DMol3中实现的数值基集可最小化甚至消除基集叠加误差。

2.   模拟结果

2.1   相互作用能

在模拟体系中，水分子的形态和流动的差异性主要是受表面活性剂水溶液和褐煤分子或者SiO₂相互作用的结果。相互作用能可以数字化体现分子之间的相互作用的强度，如果相互作用能为负值，则表明吸附过程是自发进行的；若相互作用能为正值，则反之。相互作用能的绝对值越大时，分子之间的相互作用也越强。通过分子动力学模拟，计算体系相互作用能公式如下：

式中：$ {E_{{{\mathrm{inter}}} }} $为活性剂水溶液与煤或SiO₂的相互作用能；$ {E_{{\mathrm{total}}}} $为体系内的总能量；$ {E_{\mathrm{i}}} $为体系中褐煤分子或者SiO₂的能量；$ {E_{{\mathrm{w + s}}}} $为水和表面活性剂的能量。

根据式（1）计算出平衡后体系的相互作用能，结果见表2。由表2可以看出，褐煤煤尘体系（简称WD-DS）的相互作用能的绝对值比矽尘体系（简称SiO₂-DS）的大，说明相同条件下，褐煤分子和活性剂水溶液的相互作用强度更强，褐煤煤尘表面吸附水的能力更强。

表  2  水−CTAB−煤尘/矽尘的相互作用能

Table  2.  Interaction energies of water-CTAB-coal dust/silicon dust

体系	能量/（kcal·mol−1）
	总能量	水和CTAB	粉尘	相互作用能
水−CTAB−煤尘	−2 739.640 2	−4 053.856 8	2 107.936 5	−793.720 4
水−CTAB−矽尘	−24 706.914 0	−4 065.075 7	−20 451.980 0	−135.598 0
注：1 kcal=4.184 kJ。

下载: 导出CSV 

| 显示表格

2.2   径向分布函数（RDF）

径向分布函数（RDF）被定义为特定原子在r距离内的密度与容重的比率。换句话说，特定原子的密度随与参考分子距离的变化随体积密度的变化表示RDF。所以，它可以用来演示给定分子周围的密度分布，数学表达式如下：

式中：$ {\mathrm{d}}N $为与中心的距离为$ r $到$ r + {\mathrm{d}}r $之间的原子数目，个；$ \pi $为圆周率，取3.141 592 6；$ \rho $为该原子的平均数目密度。

根据分子动力学模拟的结果，分析得到的CTAB头基和H₂O的径向分布曲线如图4所示，CTAB尾基和H₂O的径向分布曲线如图5所示。

图  4  CTAB头基和H₂O的径向分布曲线

Figure  4.  Radial distribution curves of CTAB head groups and H₂O

下载: 全尺寸图片 幻灯片

图  5  CTAB尾基和H₂O的径向分布曲线

Figure  5.  Radial distribution curves of CTAB tail groups and H₂O

下载: 全尺寸图片 幻灯片

由图4可以看出，模拟平衡后，CTAB头基和水分子在WD-DS和SiO₂-DS中的0.619 nm和0.597 nm处出现峰值，且SiO₂-DS的峰值要比WD-DS的高，0.7 nm两曲线出现交点；WD-DS在0.885 nm处又出现第二峰，交点之后WD-DS的径向分布曲线一直在SiO₂-DS的上方，两者缓慢下降。根据峰的位置判断3处峰是因为CTAB头基和水之间的静电力作用而形成的，虽然SiO₂-DS出现的峰值要高，但WD-DS总体有2个峰值，并且交点以后WD-DS的径向分布曲线高于SiO₂-DS的，说明WD-DS中头基和水的相互作用更强。

由图5可以看出，0.7 nm之前SiO₂-DS的径向分布曲线略微高于WD-DS的曲线，0.7 nm之后相反。两曲线在1.4、1.6 nm处出现强峰，显然WD-DS的CTAB尾基和水分子的相互作用强于SiO₂-DS。

2.3   均方位移和自扩散系数

均方位移MSD用来描述时间变化下目标粒子空间位置与初始位置的偏差程度。通过计算水分子的MSD和扩散系数D，可以进一步阐明表面活性剂加入对水分子运动的影响^[17]。可以将MSD对时间的曲线拟合成直线$ y = ax + b $，然后利用以下公式计算扩散系数。

式中：$ N $为粒子数；$ {r_i}\left( t \right) $为$ i $粒子在$ t $时刻的位置；$ {r_i}\left( 0 \right) $为$ i $粒子的初始位置；$ a $为拟合之后直线的斜率。

水分子的MSD拟合直线如图6所示。根据图6均方位移MSD拟合直线的斜率计算可得自扩散系数D_WD=0.004 5 nm²/ps和$D_{{\mathrm{SiO}}_2 }$=0.004 2 nm²/ps，从D的取值来看，CTAB增强了褐煤煤尘和矽尘的润湿性，使得表面和水的相互作用增强，加速了水分子的运动，进而提高了表面的润湿性，但对褐煤的效果更好。

图  6  水分子的MSD拟合直线

Figure  6.  MSD fitted straight line for water molecules

下载: 全尺寸图片 幻灯片

2.4   数密度

界面区域是相变现象发生的主要场所，是其本质就是从一种相态密度快速变化至另一种相态密度的区域。根据10−90准则^[18-19]，界面厚度被定义为密度为其体积相密度的10%~90%的2个位置之间的距离。仅对z轴方向水分子的相对数密度进行分析，WD-DS中的水分子相对数密度如图7所示，SiO₂-DS中的水分子相对数密度如图8所示。

图  7  WD-DS中的水分子相对数密度

Figure  7.  Relative number density of water molecules in WD-DS

下载: 全尺寸图片 幻灯片

图  8  SiO₂-DS中的水分子相对数密度

Figure  8.  Relative number density of water molecules in SiO₂-DS

下载: 全尺寸图片 幻灯片

由图7和图8可知：WD-DS中水的最大相对数密度为6.49，并依据10−90准则得出最大相对数密度的10%和90%位点之间的距离约为1.71、3.30 nm，由此可以算出活性剂水溶液和煤的界面厚度约为1.59 nm；同理，算得活性剂水溶液和SiO₂的界面厚度约为0.971 nm。由此可得WD-DS的固液界面层厚度大于SiO₂-DS。

3.   静电势

静电势ESP是指从无穷远处移动单位正电荷至该点时所需做的功。分析分子表面静电势可预测其容易发生吸附的位点。各个分子或晶胞静电势图如图9所示。图中：红色区域表示正电荷，该区域的电子收缩表明它很难失去电子，且可以相对容易地接受电子形成氢键；蓝色区域表示负电荷，该区域较容易提供电子形成氢键；白色区域表示电位接近于0的区域，具有相对稳定结构^[20]。各个分子或晶胞的静电势对比表见表3。

图  9  各个分子或晶胞静电势图

Figure  9.  Diagrams of the electrostatic potential of individual molecules or crystals

下载: 全尺寸图片 幻灯片

表  3  各个分子或晶胞的静电势对比表

Table  3.  Comparison of the electrostatic potential of individual molecules or crystals

分子或晶胞种类	静电势/（Ha·e−1）
	正	负	尾链
H2O	0.091 56	−0.065 44	—
CTAB	0.203 40	—	0.021 89
WD	0.124 20	−0.070 96	—
SiO2	0.568 80	−0.020 57	—

下载: 导出CSV 

| 显示表格

CTAB头基的最大正静电势0.203 40 Ha/e，大于水分子的最大正静电势0.091 56 Ha/e，可以得出 CTAB的头基是亲水的。尾链的正静电势为0.021 89 Ha/e，小于水分子的最大负静电势0.065 44 Ha/e，因此，CTAB尾链是疏水的。

平衡后，CTAB头基一部分吸附于WD表面的负静电势位置，另一部分与水进行吸附。当头基和WD表面负静电势吸附时，头基会占据原先水分子吸附的界面位置，尾基存在疏水性，朝向水中。吸附前该位置与水分子的电势差为0.162 52 Ha/e，当头基吸附后，该位置与水的电势差变为0.092 85 Ha/e，WD的亲水性减弱。CTAB尾基一部分吸附于WD表面的负静电势位置时，吸附后电势差变为0.274 36 Ha/e, WD的亲水性增强。CTAB头基尾基吸附电势差变化是增大的，对WD分子亲水性是增强的，润湿性有所提升。同理，CTAB对SiO₂的润湿作用有所提升。相比之下SiO₂加入活性剂之后的电势差增幅更大，润湿效果应该更好。

通过对分子间相互作用能对比可知，WD与水分子之间的作用能更大，润湿性更好，可能主要是因为WD分子中存在大量官能团，比如醚、醇、羧酸、酚、酮等，官能团的存在一定程度上会增加氢键的形成，进而增强分子间的相互作用力，这可能是SiO₂和活性剂水溶液的相互作用较弱的原因，导致了SiO₂润湿性差。

4.   结　语

1）基于分子动力学模拟的相互作用能和均方位移定量分析表明：在同等质量的WD和SiO₂、相同表面活性剂CTAB水溶液条件下，WD和活性剂水溶液的相互作用能为−793.720 4 kcal/mol，SiO₂的为−135.598 0 kcal/mol，即CTAB对WD的润湿性较好；结合水分子自扩散系数D验证了CTAB对WD的润湿性较好，其中D_WD=0.004 5 nm²/ps 和$D_{{\mathrm{SiO}}_2 } $=0.004 2 nm²/ps，D值越大，水的扩散效果越好。

2）基于分子动力学模拟的径向分布函数和相对数密度的定性分析表明：相较于尾基，CTAB头基和水峰值高、作用强，矽尘和水作用相比于WD弱；由相对数密度所得的固液界面厚度可知，WD和水的界面厚度1.59 nm大于SiO₂和水的0.971 nm，所以WD中沿z轴方向润湿性效果好。

3）静电势分析结果表明：表面活性剂CTAB头基是亲水头基，尾基是疏水的。活性剂吸附使得表面发生电势差的变化是影响表面润湿性的因素，CTAB对SiO₂和WD 都有润湿效果，单单考虑电势差的变化程度，SiO₂的电势差变化程度要大，但是考虑到WD官能团的作用，WD和CTAB水溶液的相互作用能大，CTAB对WD润湿性更好。