有机膦系抑制方解石阶梯面生长的分子动力学研究

陈春钰，门丽娟，居殿春，郭 雷，姚富升，张 琪，俞鸿宇

（江苏科技大学张家港校区冶金与材料工程学院，江苏张家港215600）

现代工业采用的冷却水中随着水分的不断蒸发、浓缩，水中存在的Ca2+、Mg2+、CO32-、SO42-、SiO32-等无机离子浓度逐渐增大，会出现无机物结垢、金属腐蚀及由此衍生的微生物滋生等问题［1-3］，使得冷却水的传热效率大大降低，设备正常的热交换受到阻碍，严重威胁生产安全，影响企业的经济效益［4-6］。阻垢剂能够减缓和抑制水垢产生与生长，在循环冷却水处理中，含磷类聚合物阻垢剂是较为常用的阻垢剂，但极易水解［7-8］。若滥用含磷阻垢剂，会引起蓝藻爆发、赤潮等环境问题［9-10］。因此，人们开发出一些新型有机膦系阻垢剂来弥补这一缺陷。目前国内外阻垢剂的研发、筛选仍依赖于实验或经验外推，效率低且易造成资源、经济的浪费，根源在于阻垢剂的理论机理尚未弄清［11-12］。

Gill等［13］提出了“空间匹配理论”，认为阻垢剂与成垢晶体之间的相互作用主要是库仑力，二者功能基团之间的空间距离越接近则药剂的阻垢性能越好；刚性结构的阻垢剂分子对某种特定的晶垢起作用，而柔性结构的分子可对多种晶垢起作用。Shi等［14-15］运用该理论对阻垢剂的作用机理作了进一步研究，取得了一定的成果，并根据理论模型设计出了具备优良阻垢性能的冠醚类药剂。Marine等［16］采用MC模拟研究硅酸盐聚合早期的变量和扩散系数等，并取得了相关成果。Franca等［17］用羧酸基团取代膦酸基团，并比较膦酸基团的减少对阻垢的影响，也提出抑制能力取决于结构的不同及其他阳离子的存在，这些变化会改变阻垢剂的化学环境和配位强度。

虽然前人已经对评定阻垢性能及阻垢机理研究做了大量的工作，但目前阻垢机理研究并不完善，还存在大量问题［18-21］。本文采用周期性边界条件下的分子动力学（MD）方法，结合晶体生长PBC理论［22］，计算有机膦系与方解石（104）面4个不同阶梯面之间的相互作用，对有机膦酸阻垢剂的阻垢性能及机理进行研究，详细阐述该类药剂对CaCO3垢的阻垢分散机理，为新型阻垢剂开发应用提供合理的理论依据［23-27］。

1 计算方法

1.1 有机膦系阻垢剂

将方解石晶体的（104）面切割成4种不同的阶梯面，预先设定模拟条件后，运用分子动力学方法（MD）研究含N的亚甲基膦酸分子与方解石（104）4个不同阶梯面的相互作用，即有机膦酸分子在（104）面的4个阶梯面上的吸附情况，然后根据吸附情况从中选取5种亚甲基膦酸分子作为研究对象，分别为氨基亚甲基膦酸（AMP）、氨基二亚甲基膦酸（NDP）、氨基三亚甲基膦酸（ATMP），乙二胺四亚甲基膦酸（EDTMP）和己二胺四亚甲基膦酸（HDTMP）。由于研究的冷却水系统中的水为碱性水质，水温为70～77℃，故将有机膦酸分子中的膦酸基团去质子化为PO32-，5种有机膦酸分子的结构如图1所示。

图1 5种有机膦酸化合物的去质子化几何构型Fig.1 Deprotonation geometry of five organic phosphonic acid compounds

1.2 方解石晶体生长习性计算

采用Materials Studio 6.0版本中的Morphology模块对方解石晶体的生长习性进行模拟，选用COMPASS力场及NVT系统，温度设定为77℃［28］。无机盐的表面附着能是在晶体表面加上一层生长晶面所释放的能量，表面附着能越小，则意味着该晶面越稳定、生长速率越慢，晶体生长受到此晶面的控制就越大［29］。用Growth Morphology法假定晶体的表面附着能与晶面的生长速率成比例，预测或者根据前人研究输入一组合适的晶面（hkl），计算每个晶面的表面附着能。经模拟计算，得到的方解石主要生长面的表面附着能数据见表1。

表1 方解石主要生长面的表面附着能Table 1 Surface adhesion energy of the main growth surface of calcite kJ/mol

由表1看出，方解石（104）晶面的表面附着能Eatt（Total）最小，仅为-1 552.95 kJ/mol，所以方解石（104）晶面生长最为缓慢，方解石晶体的生长主要受到（104）面控制。方解石（104）生长面的侧视图见图2。

图2 方解石（104）面侧视图Fig.2 Side view of calcite（104）

由图2也可以看出，方解石晶体的（104）面是电中性面，每层晶面都是由一个Ca2+和一个CO32-交替排列而成，相对其他面最为稳定整齐［30］，说明（104）面是方解石最终能够稳定存在的生长面，因此其生长速度相较于其他晶面更为缓慢。

1.3 方解石（104）4种阶梯面的切割分配

本研究采用Materials Studio 6.0版本中的Discover模块进行分子动力学模拟，已知方解石晶体最主要的生长面是（104）面，并且Ca2+与CO32-交替排列［31-33］。根据晶体PBC理论，将方解石晶体（104）面切割成如图3所示的阶梯面。

图3 方解石（104）面的阶梯面侧视图Fig.3 Stepped lateral view of calcite（104）face

根据Ca-Ca面、CO3-Ca面、Ca-CO3面以及CO3-CO3面4个阶梯面所构建的模拟超晶胞真空盒尺寸为：3.228 nm×1.996 nm×3.370 nm，α=β=γ=90.0°。Ca-Ca面真空箱内一共有660个原子（Ca：148，C：128，O：384）；Ca-CO3面真空箱内一共有640个原子（Ca：128，C：128，O：384）；CO3-CO3面真空箱内一共 有620个 原 子（Ca：108，C：128，O：384）；CO3-Ca面 真 空 箱 也 存 在640个 原 子（Ca：128，C：128，O：384）。将模拟盒子的真空层厚度均设置为2.0 nm，把有机膦酸分子分别放置在这4个阶梯面上得到超分子模拟体系。在超分子模拟体系最优化的基础上，模拟200.0 ps后再优化1 000步，得到最终模拟结果。由上述模拟过程以EDTMP为例，给出了模拟计算前有机膦酸分子在真空盒内的初始位置图见（图4）。

图4 EDTMP与方解石（104）面4个阶梯面模拟前的侧视图Fig.4 Side view before simulating the four stepped faces of HDTMP with calcite（104）

2 结果与讨论

2.1 5种有机膦酸分子与4种阶梯面的结合能

阻垢剂与各阶梯面之间的相互作用能ΔE可以用下面的公式表示［34］：

式中：ΔE为阻垢剂与晶面的相互作用能，kJ/mol；Esystem为模拟后方解石阶梯面与有机膦系化合物组成的超分子体系的总能，kJ/mol；Einhibitor为模拟计算后有机膦系化合物单点能，kJ/mol；Esurface为模拟后晶面的能量，kJ/mol。

阻垢剂与晶面间的相互作用能主要为静电作用、范德华力作用和极少的氢键作用能。为了总体评价有机膦酸阻垢剂的阻垢能力，用修正结合能Ecor.来衡量有机膦与方解石晶体相互作用的强弱，将方解石（104）面的4个阶梯面综合一起考虑。

式中：Ecor.为 最终总的修正结合能，kJ/mol；Ei为某一种阻垢剂分子与方解石阶梯面i的结合能，kJ/mol；Ci为权重系数（0≤Ci≤1，∑Ci=1），即阻垢剂分子与阶梯面i的结合能占体系总结合能的比例。

先将整个系统进行优化，即将AMP、NDP、ATMP、EDTMP、HDTMP等有机膦酸分子与方解石（104）面的4个不同阶梯面的超分子体系最优化，在此基础上再进行分子动力学模拟计算200.0 ps，在得到最终模拟结果之前，还需要进一步进行优化1 000步。为了比较5种阻垢剂与4个阶梯面相互作用的强弱，计算列出相互作用体系的数据，各个有机膦酸分子与阶梯面的总结合能与修正结合能见表2。

表2 5种有机膦酸分子与4个阶梯面的结合能Table 2 Binding energy of five organic phosphonic acid molecules with four stepped surfaces kJ/mol

根据表2数据，比较各分子在4个阶梯面的结合能△E，发现AMP、NDP与4个阶梯面的结合能均小于1 087 kJ/mol，与其他3种分子相比较小。表明AMP与NDP的阻垢能力较弱，这是由于这2个分子中膦酸基团个数较少。由ATMP与4个阶梯面的结合能数据来看，推测其与阶梯面的吸附程度比AMP、NDP要大，其结合能由大到小依次为CO3-Ca、Ca-CO3、Ca-Ca、CO3-CO3，说明ATMP与CO3-Ca面结合得最为紧密。进一步比较EDTMP、HDTMP与4个阶梯面的结合能数据，相对于CO3-CO3面和CO3-Ca面的结合能而言，这2种分子与Ca-Ca面和Ca-CO3面的结合能都比较小并且比较接近，两个面的结合能相差不到125 kJ/mol；而与CO3-CO3面和CO3-Ca面的结合能较大，且都达到了2 926 kJ/mol以上，尤其是HTPMP与CO3-CO3面的吸附最强，结合能达到3 568.5 kJ/mol。

比较5种有机膦酸分子与（104）面总的修正结合能Ebind，cor.，AMP与NDP的修正结合能较小，均小于1 045 kJ/mol；EDTMP与HDTMP的 修 正 结 合 能较大，均大于2 090 kJ/mol。由此得到各分子吸附在（104）面上总Ebind，cor.由强到弱的顺序为HDTMP、EDTMP、ATMP、NDP、AMP。

这一结果也表明，阻垢效果随着阻垢剂中膦酸基团个数的增加而增强，阻垢剂中膦酸基团个数要达到3个及以上时阻垢效果才会比较明显；若阻垢剂中含有2个氮原子，阻垢效果随着2个氮原子间的碳链骨架的增长而增强，较长的碳链骨架会减小分子间的空间位置，使阻垢剂分子更容易形变，阻垢剂与晶面间的吸附也会更加紧密。

2.2 HDTMP与4种阶梯面的模拟作用

选择阻垢性能最优的HDTMP为例，图5给出了HDTMP分别与Ca-Ca面、Ca-CO3面、CO3-CO3面和CO3-Ca面上的相互作用结果图，进一步探讨有机膦酸分子在4种不同阶梯晶面上的吸附情况。

图5 HDTMP在方解石（104）4个阶梯面的吸附侧视图Fig.5 HDTMP adsorbs side view on four steps of calcite（104）

由晶体生长PBC理论可知，晶面的阶梯拐角位置与拐角点是晶体的活性生长点，也就是说晶体会优先从这2个位置开始生长［35-36］。图5中明显可见，HDTMP吸附在Ca-Ca面与Ca-CO3面的阶梯拐角位置；然而对于CO3-CO3面，HDTMP吸附在右边的晶体拐角点；对于CO3-Ca面，HDTMP吸附在左边的晶体拐角点。虽然吸附位置不同，但是都吸附在阶梯拐角位置或者终端拐角点，是晶体优先生长点，说明阻垢剂分子占据晶体生长位置，从而抑制晶体生长。

观察HDTMP在Ca-Ca面上的吸附，只有一个膦酸基团吸附靠近阶梯拐角位置，其他膦酸基团都因整个分子竖立而远离晶面，因此HDTMP与Ca-Ca面间的相互作用力必定较小，HDTMP在Ca-Ca面上的吸附结合能也说明这一点；对于HDTMP在Ca-CO3面的吸附情况，由于HDTMP含有6个亚甲基，碳链骨架明显地弯曲成拱桥形状，其中一端的膦酸基团依然吸附在阶梯拐角位置，另一端的膦酸基团也能与晶面吸附得较为紧密，再加上氢键作用，所以HDTMP的阻垢效果最强；HDTMP分别吸附在CO3-CO3面的右边终端拐角点，CO3-Ca面的左边终端拐角点，在Ca-O离子键与氢键的作用下，吸附得较为紧密，并且在如此强的相互作用力下HDTMP发生了不同程度的变形。结合HDTMP与阶梯面的结合能数据，4个阶梯面吸附由强到弱顺序为CO3-CO3、CO3-Ca、Ca-CO3、Ca-Ca，与模拟结果相一致。

3 结论

采用模拟软件（Materials Studio），结合晶体生长PBC理论，并采用分子动力学方法（MD）系统地研究了常用的5种有机膦阻垢剂对方解石（104）面4个不同阶梯面的阻垢机理，得到下述主要结论：1）模拟结果图以及结合能数据表明，方解石（104）阶梯面的正常生长过程遭到破坏，阻垢剂极大程度上阻碍了晶面的生长。阻垢剂抑制方解石晶面生长的能力由大到小顺序为HDTMP、EDTMP、ATMP、NDP、AMP；2）根据分子动力学的模拟结果与结合能数据，以阻垢剂分子吸附在Ca-CO3面为例，探讨阻垢剂分子与4种阶梯面之间的相互作用。以HDTMP在4种阶梯面上的吸附为例，进一步探讨不同阶梯面对阻垢效果的影响。结果表明，终端拐角位置与阶梯拐点位置均是阶梯面的活性生长位置，HDTMP均吸附在这2个位置，且稳定性越强吸附能力越弱。