基于LAMMPS模拟岩石表面润湿性

润湿性是指不相混的两相流体与岩石固相表面接触时，其中一相流体沿着岩石表面铺开的现象，该相称为润湿相。润湿性一般采用接触角法来确定，通常根据水在固体表面的角度θ来定义系统的润湿性，接触角为0°～75°为水润湿，75°～105°为中间润湿，105°～180°为油润湿。接触角是表征液滴对固体表面的润湿程度参数，宏观上可以通过接触角实验测量，但是，宏观与微观之间存在差异，如图1所示，传统的接触角在分子水平上变得不明确。为获得接触角分子模拟和实验结果间的定量比较与普适性解释，润湿通常以表面上的液滴形状为特征，极端情况是球形和完全铺展开，介于两者之间的所有其他形状都可以认为是球顶状。

图1 宏观和微观液滴示意图

分子模拟在研究润湿性微观作用机理方面发挥着越来越重要的作用，国内外研究者已通过分子模拟对不同材料表面润湿性进行了研究。以页岩为例，不同页岩储层的矿物组成存在较大差异，一般石英含量较大，因此，以石英为代表进行研究。选择α-石英单晶胞，参数来源于晶体结构数据库（COD），晶格参数为a=b=4.913Å，c=5.4052Å，α=β=90°，γ=120°，如图2所示。切取表面并裸 露出氧原子，扩胞作为石英矿物基底模型，储层条件下的石英通常为水湿，即羟基化修饰表面。再在基底构型顶部均添加一定厚度的真空层以减小相邻模拟单元中来自周期性构型间的相互作用。

图2 石英单晶胞结构

对于纳米水滴，首先构建单个水分子并进行几何优化，再建立密度为1g/cm3的含1000个水分子的立方体盒子，测量纳米水滴质心，并将其置于表面模型之上构建初始模型，如图3所示。

图3初始模型

采用LAMMPS进行分子动力学模拟时，将体系划分为多个同心圆环，圆环的轴心沿着z方向，坐标（x,y）为纳米水滴质心坐标，圆环在z方向的切块厚度dz为1Å，在xy平面内的最小半径为1Å，最大半径为盒子边长的一半，并划分为该数目的小块，运算过程中统计每个小圆环体内的水分子数密度，统计完成后计算出每个小体积元内水的密度，并以液滴底面的中心点为起始点由远及近地统计不同半径和不同高度的液滴密度，从而得到密度云图，具体原理图如下所示：

图4 同心圆环原理图[Huang C, Xu F, Sun Y. Effects of morphology, tension and vibration on wettability of graphene: A molecular dynamics study[J]. Computational Materials Science, 2017, 139: 216-224.]

图5 最终构型

模拟结束得到最终构型如图5所示，根据成品模拟阶段计算的二维密度，作出云图如图6所示，对密度云图外轮廓拟合得到气液界面轮廓线，做出切线得到切线斜率，通过斜率即可计算出润湿接触角。

图6 水分子密度云图