剪切力在碳酸钙污垢形成机制中的分子动力学模拟研究

Molecular Dynamics Simulation of Shear Force in the Formation Mechanism of Calcium Carbonate Scaling

摘要: 本研究以CaCO₃为研究对象，运用分子动力学方法构建了碳酸钙水溶液与方解石晶面(104和1~10)的体系，模拟在不同温度(273~353 K)下剪切速率(0、10~32 s⁻¹)对该体系的影响，重点考察了剪切速率对体系的结合能(E_bind)和扩散系数(MSD)的影响。研究结果表明：温度和剪切速率对Ca²⁺和$C{O}_3^{2-}$ 在方解石两晶面上作用不同，表现为溶液中的Ca²⁺和$C{O}_3^{2-}$ 更倾向在方解石(104)晶面上有序生长，离子键在碳酸钙水溶液与方解石体系中起主要作用，而范德华力作用起辅助作用。

Abstract: This study focused on CaCO₃ as the research object. Using the molecular dynamics method, a system of calcium carbonate aqueous solution and the crystal faces (104 and 1~10) of calcite was constructed. The influence of shear rate (0, 10 to 32 s⁻¹) at different temperatures (273~353 K) on this system was simulated. The effects of shear rate on the binding energy (E_bind) and diffusion coefficient (MSD) of the system were particularly investigated. The research results indicate that temperature and shear rate have different effects on the interaction of Ca²⁺ and $C{O}_3^{2-}$ on the two crystal faces of calcite, Ca²⁺ and $C{O}_3^{2-}$ in the solution are more inclined to grow in an ordered manner on the (104) crystal face of calcite. Ionic bonds play a dominant role in the aqueous solution of calcium carbonate and the calcite system, while van der Waals forces play a supporting role.

1. 引言

循环水冷却器是工业生产中进行工艺介质冷却的重要设备。由于循环水冷却器使用的循环冷却水是含可溶性钙盐[1]的硬水，可溶性钙盐在长时间的运行过程中易生成微晶，微晶逐渐长大或者聚集后析出并逐渐沉积在金属表面形成水垢，水垢的生成很大程度上降低了循环水冷却器的热交换效率[2]，严重时可能导致非计划停车。研究表明水垢的生长特性主要取决于水中矿物质的种类和组成[3]、热量传递方式、流动情况和热力学条件等因素，同时也取决于构成传热表面的物质特性。近年来，污垢生长影响因素的研究为阻垢和抑垢提供了可靠的理论依据：如Al-Otaibi [4]等人根据污垢的不同，发现较高的表面温度将导致反可溶性盐的结垢更严重；Al-Roomi [5]等人发现升高温度，结垢量显著增加，诱导期显著缩短。Zhao [6]等人发现循环水温度升高，CaCO₃污垢的生长速度会加快；Misaghi M. [7]等人发现加速搅拌和升高温度会在更短的时间内得到更多的CaCO₃沉淀；Rahimi [8]等人发现温度较高时水垢的热阻曲线由线型转变为渐进型，晶体由文石变为向方解石；Madihi [9]等人研究了温度和不同阻垢剂对CaCO₃沉淀的影响，研究发现当温度从25℃升高到55℃时，结垢速率增加，而阻垢剂的存在使钙垢的形态发生了变化；Gao [10]等人发现增加管中流体的流速，水垢的沉积速率减小，剥蚀率变大；Raheem [11]等人对普通和功能化不锈钢基片进行了污垢沉淀试验(流速为10~30 mL/min)，发现流速和表面润湿性对非均相结晶过程有显著影响；Liu [12]等人考察了流速和Ca²⁺浓度对CaCO₃结晶的影响，发现提高流速可以促进亚稳态CaCO₃结晶的生成，且溶液中的晶体类型随Ca²⁺浓度的增加主要为球霰石；Yang [13]等人发现增强剪切速率可以促进成核，且随剪切速率的增加成核速率会达到一个最大值；Forsyth [14]等人发现搅拌和剪切流会促进二次成核，从而减小晶体的晶粒尺寸；Peng [15]等人发现剪切流始终抑制晶体–液体界面的生长。

诸多研究表明剪切力(对应流动速度)增加，污垢生长缓慢或停止生长(包括无机垢、有机垢、生物垢以及大部分晶体生长)，可见深刻认识剪切力对污垢生长的作用及其机理能为设计其它污垢抑制方法提供重要的理论启示。本研究以碳酸钙垢[16]为研究对象，采用分子动力学(MD)方法模拟剪切力对碳酸钙水溶液在方解石(104)晶面与方解石(1~10)晶面上结垢过程的影响。

2. 模型构建与模拟方法

2.1. 模型构建与优化

运用MS模拟软件进行分子动力学模拟，首先构建模型最底层建成方解石晶体的(104)晶面与方解石(1~10)晶面。据文献[17]可知方解石属于R3(-)c空间群，故本研究方解石(104)晶面超晶胞大小为：4.8558 × 2.9926 × 5.5652 nm³，共1440个原子(C: 288; O: 864; Ca: 288)，Ca²⁺和${\text{CO}}_3^{2-}$ 各4层；方解石(1~10)晶面超晶胞大小为：3.9900 × 5.1173 × 6.0161 nm³，共2400个原子(C: 480; O: 1440; Ca: 480)，Ca²⁺和${\text{CO}}_3^{2-}$ 各4层。其次构建包含1100个水分子、2个Ca²⁺离子和2个${\text{CO}}_3^{2-}$ 离子的溶液，其溶液浓度为0.1 mol/L，模型设置顶层30 Å的真空层，底层和中间层相连接的区域设置2 Å的真空层，保证体系能量合理，减少几何畸变[18]，连接三层模型即可得到初始模型，静态工况的初始模型如图1所示。

Figure 1. Calcium carbonate aqueous solution model for Calcite (104) and Calcite (1~10) crystal surface

图1. 方解石(104)和(1~10)晶面的碳酸钙水溶液模型

利用最陡下降法(Steepest descent)、共轭梯度法(Conjugate gradient) [19]和牛顿–拉弗森调整优化(Adjusted basis set Newton-Raphson)进行调优设计。

2.2. 分子动力学模拟

系统温度设定为273~353 K；通过Lees-Edwards边界条件使用等效剪切速率法将典型工业场景中湍流剪切速率转换为MD模拟中的剪切速率取0~32 s⁻¹，每隔2 s⁻¹取一次值，分为13组取值。利用COMPASS力场[20]、正则系综(NVT) [21]、Berendsen温度控制方法[22]，从这些轨迹文件中找到势能最低的轨迹文件[23]计算结合能，使用周期性边界条件，MD模拟细节参见表1。

Table 1. MD parameter settings

表1. MD参数设置

Simulation parameter	value	Simulation parameter	value
Ensemble	NVT	Themostat	Andersem
Time step	1.0 fs	Cutoff distance	1.25 nm
Forcefield	COMPASS	Shear direction	B(BC)

选取343 k、不同剪切速率下方解石(104)晶面与碳酸钙水溶液MD模拟为例，碳酸钙水溶液与方解石(104)晶面相互作用的动力学模拟如下：

Figure 2. MD simulation of 343 K, calcite (104) crystal surface with aqueous calcium carbonate solution (a) 12 s⁻¹; (b) 16 s⁻¹; (c) 20 s⁻¹; (d) 24 s⁻¹; (e) 28 s⁻¹; (f) 32 s⁻¹

图2. 343 K，方解石(104)晶面与碳酸钙水溶液的MD模拟：(a) 12 s⁻¹；(b) 16 s⁻¹；(c) 20 s⁻¹；(d) 24 s⁻¹；(e) 28 s⁻¹；(f) 32 s⁻¹

由图2可知，壁表面的CaCO₃以离子对形式存在，较小尺寸的离子反而有较高的扩散系数，Ca²⁺和${\text{CO}}_3^{2-}$ 会因热运动相遇并结合构成离子对。由于溶液浓度较低，不易缔合成较大尺寸的CaCO₃团簇，Ca²⁺和${\text{CO}}_3^{2-}$ 首先会形成离子对，不同的离子对相互结合促成较大CaCO₃团簇产生。因为存在分子间作用力，壁面附近的Ca²⁺和${\text{CO}}_3^{2-}$ 会优先吸附于壁面上，溶剂中的Ca²⁺和${\text{CO}}_3^{2-}$ 以做布朗运动为主。表面覆盖率(θ)定义为吸附离子数与表面最大吸附位点数之比，当剪切速率从12 s⁻¹增至32 s⁻¹时，表面覆盖率从0.32 ± 0.02显著下降至0.18 ± 0.01。

当离子向壁面附近扩散时，易于吸附在晶体表面，溶液中Ca²⁺与${\text{CO}}_3^{2-}$ 之间互相吸引，溶液中的离子会形成CaCO₃水合团簇结构。若增大浓度，Ca²⁺与${\text{CO}}_3^{2-}$ 形成水合团簇结构的概率增大、尺寸也更大。而团簇的尺寸越大，扩散系数越小，使其在换热器表面的吸附速率减小。

3. 结果与讨论

3.1. 系统平衡的判别

模拟中当系统的宏观物理量变化较小时，就说明系统达到了平衡[24]。本模拟体系的平衡由温度和能量平衡加以判别[25]，经计算体系弛豫后300 ps内温度和能量波动幅度均在±5%以内，符合平衡准则。此时方解石(104)晶面、方解石(1~10)晶面与碳酸钙水溶液相互作用的体系已达到平衡，并可以进行后续的模拟。

3.2. 水分子扩散

扩散系数是表征粒子扩散能力的重要参数，扩散系数越大，体系中的粒子发生定向或非定向位移的能力越强[26]。扩散系数D(t)用于描述扩散粒子的随机运动，通过分子动力学模拟分析原子轨迹，得到均方位移(MSD) [27]。均方位移是用于分析体系的扩散规律[28]，用体系内t时刻所有原子位移平方和的均值表示。

利用Einstein方程计算扩散系数，由于原子起始位置不停地移动，在每一瞬间时各原子处于不同位置，粒子的均方位移与扩散系数公式由式(1)和(2)表示[29]：

$\text{MSD}=\langle {\left[R_i\left(t\right)-R_i\left(0\right)\right]}^2\rangle$ (1)

$D=\frac{1}{6N_{\alpha }}\underset{x\to \infty }{lim}\frac{\text{d}}{\text{d}t}{\displaystyle \sum _{i=1}^{N_{\alpha }}\langle {\left[R_i\left(t\right)-R_i\left(0\right)\right]}^2\rangle }$ (2)

式中：$R_i\left(0\right)$ 表示粒子在0时刻的位置；$R_i\left(t\right)$ 表示粒子在t时刻的位置；$N_{\alpha }$ 是体系中粒子的数量。

由统计学可知，只要模拟粒子数量足够多，计算时间足够长，计算所得的平均值应与系统的任何一个瞬间相同，即系统的任何一个瞬间都可以作为计时零点。由Einstein方程得：

$\underset{x\to \infty }{lim}\frac{\text{d}}{\text{d}t}{\displaystyle \sum _{i=1}^{N_{\alpha }}\langle {\left[R_i\left(t\right)-R_i\left(0\right)\right]}^2\rangle }=6N_{\alpha }D$ (3)

式中D为粒子的扩散系数，假设计算时间足够长时，均方位移曲线的斜率即为$6N_{\alpha }D$ 。因此，扩散系数可由其斜率求出，通过模拟从而得到MSD随t的变化曲线。

图3为273~353 K时，各温度下水分子的扩散系数随剪切速率的变化曲线，为了便于对比，计算每组温度条件下剪切速率为0 s⁻¹时水分子的扩散系数。由图可见，在方解石(104)晶面上：水分子的扩散系数随着剪切速率的增大而增大。团簇状水分子因氢键结合，当温度升高时，水分子动能上升较快，平衡位置的振幅变大，导致部分氢键被破坏，使水分子数量增加，其扩散系数也随之增大。由于水分子在CaCO₃水溶液中的活性增强，Ca²⁺与${\text{CO}}_3^{2-}$ 相互接触的机率降低，Ca²⁺与${\text{CO}}_3^{2-}$ 很难形成微小的CaCO₃晶胚，此时增大剪切速率会阻碍碳酸钙晶体的生长。随剪切速率的增大，水分子的扩散系数呈线性关系。剪切速率上升会使溶液中的水分子更易扩散，活性增强，运动更激烈，更多束缚水脱离束缚以自由水的形式存在[29]，阻碍Ca²⁺与${\text{CO}}_3^{2-}$ 相互吸引形成团簇。对比剪切速率为0 s⁻¹时水分子的扩散系数，剪切速率增大使得水分子的扩散系数明显增大。

在方解石(1~10)的晶面上，不同温度条件下水分子扩散系数会随剪切速率的增大而增大。在32 s⁻¹时，不同温度下水分子自扩散系数随温度的升高而增大，在70℃时扩散系数达到最大值1.31 cm²/s；在0℃为最小值1.04 cm²/s时，CaCO₃的溶解度较低，在水分子数量足够多时，因水分子活性增强，Ca²⁺和${\text{CO}}_3^{2-}$ 接触反应的机率降低，因而很难生成微小的CaCO₃微晶，所以在32 s⁻¹时高温对CaCO₃晶体的生长具有一定阻碍作用。

体系中的最底层为方解石(1~10)晶面，因强极性的水分子存在，静电作用导致水分子的扩散系数偏小。

3.3. Ca²⁺与$C{O}_3^{2-}$ 的扩散性质

设定不同温度和剪切速率，用Einstein方程计算出MSD曲线，得到各曲线的斜率，得到Ca²⁺与${\text{CO}}_3^{2-}$ 在水溶液中的扩散系数，如图4和图5所示。

Figure 3. Self-diffusion coefficients of H₂O at different shear rates and temperatures: (a) Calcite (104) crystal plane; (b) Calcite (1~10) crystal plan

图3. 273~353 K温度下不同剪切速率H₂O的扩散系数：(a)方解石(104)晶面；(b)方解石(1~10)晶面

Figure 4. Diffusion coefficients for Ca²⁺、${\text{CO}}_3^{2-}$ : (a) 283 K; (b) 298 K; (c) 323 K; (d) 343 k

图4. 方解石(104)晶面上：Ca²⁺、${\text{CO}}_3^{2-}$ 的扩散系数：(a) 283 K；(b) 298 K；(c) 323 K；(d) 343 K

Figure 5. Diffusion coefficients for Ca²⁺、${\text{CO}}_3^{2-}$ : (a) 283 K; (b) 298 K; (c) 323 K; (d) 343 K

图5. 方解石(1~10)晶面上：Ca²⁺、${\text{CO}}_3^{2-}$ 的扩散系数：(a) 283 K；(b) 298 K；(c) 323 K；(d) 343 K

由图4可知，方解石(104)晶面上：Ca²⁺与${\text{CO}}_3^{2-}$ 扩散系数的平均值均相同，这是因为方解石(104)晶面不带电，Ca²⁺与${\text{CO}}_3^{2-}$ 在晶面附近做扩散运动，晶面对离子的吸引程度相同。

由图5可知，在方解石(1~10)晶面上，Ca²⁺与${\text{CO}}_3^{2-}$ 的扩散系数大于剪切速率为零时的扩散系数。剪切速率对Ca²⁺与${\text{CO}}_3^{2-}$ 在水中的扩散系数影响明显，最大值出现在20 s⁻¹处，最小值出现于10或32 s⁻¹处。由于方解石(1~10)晶面带正电，受引力作用的${\text{CO}}_3^{2-}$ 在晶面附近做扩散运动，而受斥力作用Ca²⁺的则远离晶面。

3.4. $C{a}^{2+}-C{O}_3^{2-}$ 与不同晶面的结合能

方解石(104)晶面与溶液中的结合程度${\text{Ca}}^{2+}{\text{-CO}}_3^{2-}$ 可以通过结合能来判断分析，结合能则是通过体系的相互作用来获取的，相互作用能$\Delta E$ 由式(4)表示[30]：

$\Delta E=E_{\text{tot}}-\left(E_{\text{surface}}+E_{\text{o}}\right)$ (4)

式中：E_tot是${\text{Ca}}^{2+}{\text{-CO}}_3^{2-}$ 与Ca(104)晶面的总能量，E_surface是方解石晶面的单点能，E_o是溶液中${\text{Ca}}^{2+}{\text{-CO}}_3^{2-}$ 的单点能，待系统平衡，去掉全部H₂O，目的是消除H₂O对计算$\Delta E$ 所产生的影响。定义相互作用能$\Delta E$ 为结合能E_bind的负值[31] [32]，即$E_{\text{bind}}=-\Delta E$ 。设置温度为273~353 K，计算不同剪切速率下${\text{Ca}}^{2+}{\text{-CO}}_3^{2-}$ 与方解石(104)晶面的结合能，结合能越大，${\text{Ca}}^{2+}{\text{-CO}}_3^{2-}$ 与方解石(104)晶面结合越牢固。

由图6可知，在方解石(104)晶面上：在273~353 K时体系的结合能随剪切速率的增大先逐渐升高，后缓慢减小。在本文所考察的温度下，结合能随剪切速率的变化均为非线性关系，在非零的中间速率下达到最大值，即在${\dot{\gamma }}_0$ (剪切速率)附近出现拐点。

Figure 6. Fluctuation curves of the system binding energy with shear rate: (a) 273 K; (b) 283 K; (c) 298 K; (d) 313 K; (e) 323 K; (f) 333 K; (g) 343 K; (h) 353 K

图6. 方解石(104)晶面系统结合能随剪切速率变化的波动曲线：(a) 273 K；(b) 283 K；(c) 298 K；(d) 313 K；(e) 323 K；(f) 333 K；(g) 343 K；(h) 353 K

Figure 7. Fluctuation curves of the system binding energy with shear rate: (a) 273 K; (b) 283 K; (c) 298 K; (d) 313 K; (e) 323 K; (f) 333 K; (g) 343 K; (h) 353 K

图7. 方解石(1~10)晶面系统结合能随剪切速率变化的波动曲线：(a) 273 K；(b) 283 K；(c) 298 K；(d) 313 K；(e) 323 K；(f) 333 K；(g) 343 K；(h) 353 K

由图7可知，在方解石(1~10)晶面上，结合能与剪切速率之间呈现非线性关系，剪切速率较低碳酸钙溶液系统结合能随剪切速率增大而增大，达到峰值后，系统结合能随剪切速率的增大而减小，这是由于平流增强与机械应变而导致能量势垒增加。剪切力对污垢生长的影响体现在分子运动和团簇生长的微观层面上。在低剪切速率下，系统结合能随剪切速率的增加而增大，这是因为离子的扩散速率增加；随着剪切速率的进一步增加，成核能垒也随之增加；当超过最大值后，随着剪切速率的进一步增加，系统结合能开始下降，晶面与${\text{CO}}_3^{2-}$ 之间的相互作用要比Ca²⁺、${\text{CO}}_3^{2-}$ 之间的相互作用强得多，这使得${\text{CO}}_3^{2-}$ 首先被吸收到晶体表面。

成核率存在一个峰值，低剪切速率下成核速率加快，达到峰值后，成核速率随着剪切速率的进一步增大而减小。对这种效应的机理性解释，即分子向原子核运输的平流增强与机械应变导致的能量势垒增加之间的竞争[18]。

3.5. 不同工况对两模型中$C{a}^{2+}-C{O}_3^{2-}$ 与方解石晶面结合能的影响。

将方解石(104)晶面与碳酸钙水溶液相互作用的体系规定为模型一，方解石(1~10)晶面与碳酸钙水溶液相互作用的体系规定为模型二，通过对比模型一与模型二，研究不同温度与不同剪切速率对两模型中污垢生长的影响。

图8为不同温度与不同剪切速率下模型一与模型二中离子与晶面的结合能。可以看出，两模型的结合能随剪切速率的变化趋势基本一致，即结合能随剪切速率先增加至最大值，然后又随剪切速率的增大而减小。在298 K，当剪切速率大于14 s⁻¹时，模型一中离子与晶面的结合能高于模型二(14 s⁻¹时：E_bind,I = 20.00 kcal/mol，E_bind,II = 20 kcal/mol)，污垢更容易在方解石(104)晶面上生长。同理，在283、333和343 K时，模型一中离子与晶面的结合能基本均高于模型二，因此，污垢更易在方解石(104)晶面上生长。

Figure 8. Binding energies of Ca²⁺, ${\text{CO}}_3^{2-}$ and calcite crystal surface in model I and model II. (a) 283 K; (b) 298 K; (c) 333 K; (d) 343 K

图8. 模型一与模型二中${\text{Ca}}^{2+}{\text{-CO}}_3^{2-}$ 与方解石晶面的结合能：(a) 283 K；(b) 298 K；(c) 333 K；(d) 343 K

4. 结论

本文模拟了不同温度、不同剪切速率下碳酸钙水溶液与方解石(104)晶面与方解石(1~10)晶面的相互作用。结果表明：

1) 在方解石(104)晶面上：温度和剪切速率的增大可促进水分子在溶液中的扩散，表现为水分子运动加剧，而高温和高剪切速率能抑制Ca²⁺和${\text{CO}}_3^{2-}$ 相互吸引形成团簇，不利于结垢的形成；当剪切速率较小时，剪切速率的增大在一定范围内能增加Ca²⁺、${\text{CO}}_3^{2-}$ 与方解石(104)晶面的牢固性，而超出这个范围，升高剪切速率反而会降低Ca²⁺、${\text{CO}}_3^{2-}$ 与方解石(104)晶面的结合程度。

2) 在方解石(1~10)晶面上：升高温度与增大剪切速率能促进水分子扩散，导致水的自由扩散系数增大，减小了Ca²⁺和${\text{CO}}_3^{2-}$ 相互接触反应的概率，使得碳酸钙污垢形成受到阻碍；Ca²⁺与${\text{CO}}_3^{2-}$ 的扩散系数受剪切速率的影响明显，均在20 s⁻¹处取得最大值，在此剪切速率时，离子活性最大，形成离子对的可能性最大；${\text{CO}}_3^{2-}$ 在带正电方解石(1~10)晶面受到强烈的静电作用扩散运动，而Ca²⁺在方解石(1~10)晶面则受静电作用远离晶面。

3) 当系统的剪切速率较低时，系统结合能会随剪切速率增大而增大，当达到峰值后，系统结合能随剪切速率的增大而减小。

基金项目

本课题得到了国家自然科学基金(项目编号：52176074)的资助。