表面电荷如何改变水-十烷-二氧化硅系统中的油接触角
水-油-岩体系中油润湿性对岩石表面吸附离子及化学剂诱导的表面电荷演化极其敏感。通过大规模分子动力学模拟，我们揭示了表面电荷对理想水-十烷-二氧化硅系统中油接触角的影响。结果表明油纳米液滴接触角高度依赖表面电荷。当表面电荷密度超过临界值0.992 e/nm²时，接触角达78.8°，呈现明显亲水状态。接触角变化可从油分子数密度分布得到确证。随表面电荷密度增加，油分子吸附减弱，纳米液滴与二氧化硅表面接触面积减小。此外，数密度分布、RDF分布和分子取向表明油分子在二氧化硅表面以平行取向层层吸附。然而，在较高表面电荷密度下，油分子靠近二氧化硅表面的层状结构逐渐变得模糊。 随着原油开采难度增加，提高采收率（EOR）愈发重要[1,2]。水驱因成本低、效率高成为主导手段。水驱是多尺度过程，涉及多孔储层中的指进现象或微流道中原油被注入水驱替的两相流过程。因此水驱受孔隙到储层尺度各因素控制。鉴于油藏大比表面积和低毛细管数，岩石表面油水润湿性对油水两相流影响巨大，对采收率至关重要[3,4,5]。 通常亲水岩面利于驱油，但岩面初始常因原油极性组分吸附呈亲油状态。因此将岩面润湿性从亲油转为亲水对EOR技术关键。水驱中常人工向注入水添加化学剂（表面活性剂、聚合物、离子、纳米颗粒）实现润湿性改变[6-11]。例如溶解表面活性剂可使储层岩石转为更亲水态；公认离子对形成和表面活性剂分子通过与岩面吸附原油组分相互作用是表面活性剂诱导润湿性改变的两大机制[12]。低盐水驱中，降低盐度致双电层膨胀和多组分离子交换是砂岩储层润湿性改变的两种竞争机制[13,14]。纳米颗粒水驱中，三相接触线区域纳米颗粒自层化和低维结构化对油相施加结构楔入压力，最终促使亲油岩面转为亲水态[15-18]。 无论何种水驱技术，润湿性改变的本质机制是注入化学剂诱导的水-油-岩三相相互作用变化。水-油-岩相互作用中静电作用影响深远，因化学剂吸附岩面改变表面电荷。例如二价阳离子促进酸性油组分与负电岩面结合，将岩面电势由负转正[19-22]。水化学本身也影响岩面电荷[23]。因此表面电荷无疑是影响岩面润湿性的最重要因素之一。Puah等人[24]表明表面电荷影响固体表面静态润湿性和润湿动力学。他们进一步证实静态接触角在零电荷点上下按Lippmann方式变化，动态接触角数据可用分子动力学理论描述。 因三相润湿性直接关联水-油-岩分子相互作用，分子层面研究揭示表面电荷对岩面润湿性影响十分必要[25,26]。本文构建二氧化硅表面模拟硅酸盐岩面，考虑表面电荷，用分子动力学（MD）模拟研究其三相润湿性。期望从分子数密度分布、岩面分子结构等方面揭示表面电荷影响及其深层机制。 我们研究模型体系中三相润湿性：油模型为十烷，岩石模型为二氧化硅。二氧化硅表面能良好代表砂岩储层[11,27,28]。为计算油接触角，油被包裹在水中形成圆柱液滴（见图1）。圆柱油滴可消除线张力对接触角计算的影响。构建二氧化硅衬底时，将一组α-石英三角晶胞分别沿[100]、[010]、[001]方向组装。然后将α-石英晶体从菱面体转为边长15×15×1 nm³的立方体。为获得真实表面结构，参考Puibasset等人[29]程序，晶体先沿(111)晶面切割，移除不完整四面体中的硅原子。然后用氢原子饱和仅与一个硅原子键合的非桥氧原子。即二氧化硅衬底顶底两面均饱和氢原子。因此我们通过向二氧化硅表面顶层氢原子添加部分电荷改变表面电荷，同时在底层氢原子添加等值相反电荷保持体系电中性。十烷和水分子结构见图1。 模拟体系示意图。油立方体被水包围置于二氧化硅衬底顶部（左），十烷和水分子原子视图（右） 如图1所示，水和油初始排列在二氧化硅衬底顶部，油立方体被三个不同大小的水立方体包围。油立方体尺寸5×15×5 nm³，位于衬底中央。左右水立方体尺寸与油立方体相同，顶部水立方体尺寸15×15×5 nm³。水和十烷分子数由立方体尺寸及300 K、1 bar下密度决定。水和油密度分别为1.000和0.725 g/cm³，对应水分子58,319个，十烷分子1,150个。分子排列在Material Studio软件实现，初始随机分布在特定盒子中。随后进行能量最小化计算优化分子几何。模拟达平衡态时，油盒演变为圆柱液滴，三个水盒融合形成包围油滴的连续相。 MD模拟在LAMMPS（大规模原子/分子大规模并行模拟器）平台进行。模拟在NVT系综下，温度300 K。总共运行500万步，时间步长1 fs。每10,000步记录原子坐标以观察油滴演变并最终计算接触角。体系约200万步达平衡态（图2所示）。因此采用最后100万步原子坐标平均分子数密度分布。x-和y-方向周期性边界条件，z-方向反射边界条件。模拟期间体相SiO₂原子固定，表面O和H原子在Nose-Hoover温控器控制下灵活。 圆柱油滴接触角计算流程 水-十烷-二氧化硅原子相互作用均用12-6 Lennard-Jones势耦合极性项建模，同时考虑范德华和库仑力[30-32]： 其中ε为能量参数，σ为长度参数，qᵢ和qⱼ为原子i和j电荷，C为静电常数，χ为介电常数。R_cut为范德华和短程静电力截断距离，设为10 Å。模拟中用粒子-粒子粒子-网格法处理长程静电力。式(1)势参数获取：水分子采用TIP3P模型[33]；十烷分子采用OPLS-AA力场[34]；二氧化硅衬底采用CVFF力场[35]。模拟体系非键合势参数列表见表1。跨原子势参数用Lorentz-Berthelot规则获得。单分子内原子间键合相互作用和四原子二面角相互作用也妥善考虑。键合相互作用用谐振子模型，二面角相互作用用OPLS势模型[36]。 二氧化硅表面润湿性以油接触角表征，基于纳米液滴中油分子密度分布获得。模拟早期油滴从立方体演变为半圆柱并达平衡；平衡过程通常耗时200万步，可从油滴质心z方向随时间变化确认。平衡态下质心z坐标保持恒定。因需平均平衡油水构型的油分子密度分布，为保险起见，仅收集最后100万步原子坐标获取密度分布。平均分子数密度时，将油/水区划分为大小2×2×2 ų的规则立方网格。按时间平均每个网格出现的原子数，获得油纳米液滴密度分布的3D等值面和x-z平面2D等值面。2D等值面通过沿y方向进一步平均3D等值面获得。利用2D等值面按以下规则识别油水界面：界面处网格油原子数密度为油滴体相网格密度（ρ_b）的一半。考虑界面厚度，密度在0.2ρ_b至0.8ρ_b范围内的网格定义为界面。识别离散界面点后，拟合为圆计算油滴接触角。拟合需得圆心(x₀,z₀)和半径r。接触角θ = arccos(d/r)，其中d为圆心与岩面基线z方向距离（d = z_b − z₀，z_b为基线z位置，z_b = 12.16 Å）。若距离为负值，圆心在岩面上方，接触角必超90°。图2显式描绘了接触角完整计算流程。该方法在笛卡尔坐标系下分析圆柱油滴即可得接触角。 因缺乏水包油滴在二氧化硅衬底上接触角有效数据，我们通过对比球形水滴在二氧化硅衬底上接触角验证模型。单独模拟获得4179分子水纳米液滴接触角。验证模拟中二氧化硅表面羟基固定，同其他代表性工作。测得接触角7.8°，与Pafong等人[37]报道7.0°吻合，相对偏差仅10.3%。该良好一致性表明模拟模型和方法可靠。 计算圆柱油滴接触角，获得接触角随表面电荷密度变化曲线（图3）。清晰可见接触角高度依赖表面电荷密度；从负电荷表面到正电荷表面，油滴接触角显著增大。这意味着岩面随表面电荷从负值增至正值，从亲油态转为亲水态。不同表面电荷密度下x-z平面平均密度等值线也插入图中，直观展示表面润湿性变化。 接触角随表面电荷密度变化。不同表面电荷密度下油滴等值线图亦插入图中 图3还显示：当表面电荷密度增至0.992 e/nm²时，岩面呈现极明显亲水态。低于此临界值，接触角随电荷密度变化不明显；超过临界值，接触角达78.8°，亲水态极显著。认为接触角随表面电荷变化源于水-二氧化硅和油-二氧化硅相互作用强度变化，因表面电荷影响水-二氧化硅相互作用和油-二氧化硅相互作用。特别注意：基于此理想体系的模拟结果不能直接与实验研究对比以揭示油-水-岩体系中盐度依赖润湿性。实验中水相阳离子可能将岩面电势由负转正，但绝非唯一影响。阳离子其他潜在影响包括多离子交换和双电层膨胀等[21,38-40]。因此改变表面电荷与改变水中阳离子浓度完全不同。接触角随表面电荷密度和阳离子浓度变化趋势可能不同。本模拟研究仅探讨二氧化硅表面部分电荷对理想油-水-岩体系三相润湿性影响。这些结果不能直接揭示低盐水驱机制并指导其EOR应用，但能提供表面电荷如何影响三相润湿性的基础理解。本研究仍具意义，因其揭示了表面电荷在水-十烷-二氧化硅体系三相润湿性中的深远作用。 为揭示电荷密度依赖接触角机制，分析油分子沿高度方向（z方向）密度分布。图4a给出沿高度方向每层宽0.2 nm的油分子平均密度。可见圆柱油滴底部（近衬底表面）油密度高且呈波动分布；这意味着油分子在衬底表面层状吸附。油滴顶部平均密度降至近零，因圆弧形顶部使平行层中油分子数随层接近顶部递减。同时发现，随表面电荷由负转正，吸附层数密度逐渐降低；表明油分子吸附强度随表面电荷增加而减弱。吸附减弱意味着衬底表面与油滴接触面积减小，从而油滴接触角增大。油滴高度随表面电荷变化也反映接触角增大。图4a显示油滴高度随表面电荷增加而增大；意味着高表面电荷下油滴更细长，接触角更大。正表面电荷下油分子吸附减弱可从水-十烷-二氧化硅体系能量确认。图4b显示时间平均体系能量随表面电荷增加而增大。正体系能量越高，油分子与二氧化硅吸引力越弱。插图显示中性二氧化硅表面体系能量随时间变化。 a 油分子沿z方向密度分布；b 不同表面电荷密度下时间平均体系能量 接触角随表面电荷变化也可从油分子沿x方向分子数密度分布反映。图5显示z方向高度1.9 nm处x-y平面内油分子沿x方向数密度分布呈中间平台状。平台宽度关联油滴x-y平面横截面积，高度关联油滴内分子数密度。不同表面电荷下平台宽高显著变化，对应不同圆柱油滴形状及接触角。低接触角平台较宽；高接触角平台较窄。注意表面电荷密度7.936 e/nm²时平台高度特别低，因近二氧化硅表面密度分布大幅波动（见图4a），某些情况下1.9 nm高度分子密度恰位于密度分布低谷。 z方向高度1.9 nm处x-y平面内油分子沿x方向数密度分布 我们还获得圆柱油滴内部原子负电荷密度分布（图6）。负电荷分布极不均匀。靠近二氧化硅表面区域负电荷密度极高，远离表面区域负电荷均匀分布。且靠近表面高密度区负电荷呈层状分布。负电荷层状分布直接关联油分子近表面层状结构，因负电荷直接关联油分子碳原子。此外，不同表面电荷密度下负电荷密度分布略有差异。随表面电荷密度增加，负电荷层状分布逐渐模糊。这关联纳米液滴与二氧化硅表面接触面积减小；高表面电荷密度下接触面积减小，溶剂化相互作用减弱致油分子层状结构不明显。 圆柱油滴内负电荷密度分布x-z平面。a −1.984 e/nm²。b −0.992 e/nm²。c −0.1984 e/nm²。d 0.0 e/nm²。e 0.1984 e/nm²。f 0.992 e/nm²。g 1.984 e/nm²。h 3.968 e/nm²。i 5.952 e/nm²。j 7.936 e/nm² 水-十烷-二氧化硅体系分子结构对理解三相润湿性也至关重要。我们从径向分布函数（RDF，g(r)）[41]和分子取向分析油水分子结构。图7给出中性二氧化硅表面上油水分子RDF分布和油分子取向。O-O和C-C原子对RDF分布截然不同（图7a），即C-C原子RDF波数更多、波峰值更高。结合RDF物理意义，易推断油分子在二氧化硅表面强吸附形成多层，水分子弱吸附仅现单层。该吸附现象也可见图7a插图，显示油水分子在中性二氧化硅表面平衡快照。图4a数密度分布确认油分子明显吸附层。通常液体分子在固体表面层状结构关联液-固相互作用。强液-固相互作用下液体分子沿法向表面层层分布；弱相互作用下层状结构不明显[42]。本研究中水-二氧化硅相对较弱，故油分子近固面层状结构比水分子更明显。 中性二氧化硅表面分子结构。a 油水分子RDF分布。b 不同z方向高度油分子取向 为进一步反映长链油分子在二氧化硅表面取向，我们获得不同z方向高度油分子与二氧化硅表面夹角（图7b）。角度θ定义为二氧化硅表面基线与油分子中心线夹角，中心线定义为C1-C6连线与C1-C10连线的中线。油分子碳原子序号从靠近二氧化硅表面一端开始。角度θ定义亦插入图7b。可见大多数油分子角度极小，尤其是靠近表面油分子。意味着油分子以平行表面取向吸附。图7a插图也可确认此现象。 利用MD模拟研究表面电荷密度对水-十烷-二氧化硅体系油接触角影响。结果显示：随表面电荷从负值增至正值，油纳米液滴接触角显著增大。表面电荷密度超临界值0.992 e/nm²时，接触角达78.8°，亲水态极显著。接触角变化可从油分子沿x-和z-方向密度分布确证。吸附层油分子数密度减小和油滴高度增加均意味着高表面电荷密度下油滴接触角增大。因油分子近二氧化硅表面层状结构，油滴内靠近二氧化硅表面负电荷亦层层分布。随表面电荷密度增加，负电荷（油分子）层状分布逐渐模糊，因油-二氧化硅相互作用减弱。油水分子RDF分布和分子取向分布也显示油分子以平行表面取向层层吸附。这些结果对理解水相离子及化学剂诱导润湿性改变，特别是理解表面电荷在水-十烷-二氧化硅体系三相润湿性中的深远作用具有重要意义。 一致价势场 提高采收率 大规模原子/分子大规模并行模拟器 分子动力学 液体模拟优化势全原子力场 径向分布函数摘要
背景
方法
模拟体系
MD模型
接触角计算
结果与讨论
接触角
密度分布
分子结构
结论
缩略语
纳米材料