孙志杰, 史培阳, 范蕾
东北大学 多金属共生矿生态化冶金教育部重点实验室, 辽宁 沈阳 110819
收稿日期：2022-02-21
基金项目：多金属共生矿生态化冶金教育部重点实验室开放项目(NEMM2019002)。
作者简介：孙志杰(1996-), 女, 内蒙古赤峰人, 东北大学硕士研究生;
史培阳(1968-), 男, 辽宁铁岭人, 东北大学副教授。

摘要：如何有效控制溶液体系下石膏晶体的生长行为已成为脱硫行业亟待解决的问题之一, 而晶体生长过程的有效信息往往无法及时有效在线检测.利用分子动力学模拟技术, 计算出石膏晶体特定晶面的附着能、面心距离和面积等参数, 确定了其形态学上重要的生长晶面及其表面结构, 理论研究了氯化钠溶液与石膏晶体主要生长面之间的相互作用机制.并通过实验得到氯化钠对石膏晶体形貌的影响, 与分子动力学模拟结果进行对照分析.结果表明, 由于氯化钠溶液对石膏晶体不同晶面的作用变化, 可以显著改变石膏晶体形态, 模拟结果与实验一致, 分子动力学模拟可以为晶体的形态控制提供理论基础.
关键词：石膏晶体溶液体系晶体生长习性溶剂效应分子动力学
Molecular Dynamics Simulation on Effect of Sodium Chloride on Growth Habit of Gypsum Crystal
SUN Zhi-jie, SHI Pei-yang, FAN Lei
Key Laboratory for Ecological Metallurgy of Multimetallic Ores, Ministry of Education, Northeastern University, Shenyang 110819, China
Corresponding author: SHI Pei-yang, E-mail: shipy@smm.neu.edu.cn.

Abstract: How to effectively control the growth behavior of gypsum crystal in solution system has become one of the urgent problems to be solved in desulfurization industry, and the effective information of crystal growth process can't be detected online timely and effectively. In this work, molecular dynamics simulation (MS) was used to calculate adhesion energy, face-center distance and area of specific plane of gypsum crystal, to determine morphologically dominant crystal planes and surface structure, and to theoretically study the interaction mechanism between sodium chloride and the main growth planes. The effect of sodium chloride on the crystal morphology of gypsum was obtained through experiments, which was compared with the results of MS. Results showed that sodium chlorides have varied effects on each crystal planes of gypsum crystals, which further changing the crystal morphology of gypsum. The simulation results are consistent with the experiment. MS can provide a theoretical support for morphology control of gypsum crystals.
Key words: gypsum crystalsolution systemcrystals growth behaviorsolvent effectmolecular dynamics
钢铁工业作为中国的重要基础产业, 在国民经济中发挥重要作用^[1-2].国家统计局公布的数据显示, 2021年我国粗钢产量达到10.33亿t, 按中国钢铁协会公布的“吨钢SO₂排放量1.95 kg”测算, 2021年钢铁工业SO₂的排放量达201万t, 其中70%以上来自烧结工序^[3].因而, 烧结烟气脱硫已成为我国钢铁行业污染物控制的重中之重.石灰石-石膏湿法是技术相对成熟、应用最广的烟气脱硫工艺^[4], 但在湿法脱硫过程中可溶性盐的循环富集严重影响到石膏晶体的生长行为, 导致石膏晶体结构和生长形态发生显著变化, 不仅影响到系统脱硫效率和系统运行安全^[5], 还涉及到烧结烟气脱硫石膏的资源综合利用^[6-7].随着我国对环保监管力度的加大, 烟气脱硫和脱硫石膏产量也逐年提升, 但关于脱硫石膏晶体结构与应用之间的影响作用研究报道较少, 而关于可溶性盐对石膏晶体影响作用方面的文献报道更少^[8].深入探究盐溶液对石膏晶体生长行为的作用机制对其资源化综合利用具有重要意义.
在计算机技术和量子化学迅速发展的大环境下, 利用分子动力学模拟方法研究复杂多相溶液体系下无机晶体界面现象已然成为现实, 通过在微观尺度上针对多相界面现象的实时研究, 有助于获取无机晶体溶解和再结晶过程中的相关信息, 进而有利于实现对矿物溶解行为和再结晶过程中的界面调控^[9-10].石膏晶体作为自然界常见无机矿物, 其溶解和再结晶行为与环保和资源循环等一系列问题密切相关.因此, 研究石膏晶体在水溶液中赋存行为和演变机制对于其晶体结构和应用具有重要的实际意义.晶体的溶解和再结晶行为主要取决于两个因素, 即晶体内部结构和外部环境.Cui等^[11]利用分子动力学模拟, 在90 K和氩气保护下, 借助Lennard-Jones势计算了液体分子间、液体分子和固体分子间的相互作用关系, 得到杨氏方程的近似公式, 可定性分析纳米尺度粒子间作用系数.Chen等^[12]通过分子力学和分子动力学模拟, 研究了TNP与溶液分子间相互作用关系, 分析了表面分子和周期键链的排列方式, 并发现了在不同溶剂作用下TNP晶体形貌可发生明显变化, 其中对TNP结晶生长习性影响最大的晶面是{100}面.Chen等^[13]通过分子动力学模拟了Li₂CO₃晶体与水溶剂层之间的界面效应, 通过建立改进的AE模型, 研究了Li₂CO₃在水溶剂下晶体的生长过程, 计算出不同溶剂下附着能的变化规律, 因此, 溶剂效应是影响晶体在生长过程中形貌发生转变不可忽视的重要原因.
Material Studio^[14]作为在材料科学等领域被广泛应用的分子动力学模拟软件, 在探究外部环境对晶体生长的影响时有重要作用.本研究借助Material Studio模拟软件, 探究盐溶液条件下石膏晶体晶面的稳定性, 考察晶面结构演变规律, 通过计算晶面的表面能和建立双层结构模型, 阐明氯化钠溶剂分子与石膏晶体主要生长面间的作用关系, 并通过实验得到了氯化钠溶液对石膏晶体形貌的影响, 与分子动力学模拟预测结果一致, 可为石膏晶体结构调控提供理论基础.
1 分子动力学的模拟研究1.1 建立模型借助Findit数据库, 可找到石膏晶体晶胞参数, 其属于单斜晶系, 空间群为C12/C1, 晶胞参数如表 1所示, 根据石膏CaSO₄·2H₂O晶胞参数, 利用Material Studio软件建立晶胞结构, 如图 1所示.
表 1(Table 1)

表 1 石膏晶体晶格常数Table 1 The lattice constant of gypsum crystal 晶体名称 晶胞参数 空间群 晶体系统 a/nm b/nm c/nm α/(°) β/(°) γ/(°) CaSO4·2H2O 0.569 29 1.524 5 0.654 40 90 118.480 90 I121 I121

表 1 石膏晶体晶格常数 Table 1 The lattice constant of gypsum crystal

图 1(Fig. 1)

图 1 石膏晶胞结构图Fig.1 Gypsum crystals cell structure diagram

首先切割晶体主要表面, 创建宽度、厚度均为两个分子层的超晶胞, 并将其重建为三维周期性盒子, 对晶面结构进行优化；然后绘制H₂O分子、Na原子和Cl原子, 在此初级模型基础上, 使用Clean功能进行优化, 使其具备正确的原子间距、键角等基本属性, 再使用Forcite模块对Na原子和Cl原子模型进行电子及力场分配, 可以得到Na⁺和Cl^-；利用Amorphous Cell模块搭建氯化钠水溶液模型, 在Density模块中分别输入质量分数为2%, 4%, 6%, 8%, 10%氯化钠溶液在25 ℃条件下的密度, 根据氯化钠质量分数和密度来确定在Composition模块中H₂O, Cl^-和Na⁺的比值, 输入其个数, 确保溶剂模型与超晶胞的a, b值相近, 创建不同质量分数氯化钠溶剂盒子, 如图 2所示.
图 2(Fig. 2)

图 2 氯化钠溶液的初始构型Fig.2 The initial model of NaCl solution

利用Build Layers模块, 将溶剂层和晶面层分别作为双层结构模型的上、下层, 建立起层状结构单元, 得到模拟计算模型, 为了避免形成层状周期结构、消除层与层之间的相互影响作用, 需要在溶剂层上插入一较厚的真空层(3 nm), 如图 3所示.
图 3(Fig. 3)

图 3 石膏晶体-NaCl溶液双层结构模型Fig.3 The double layer structure model of gypsum crystal-NaCl solution

本研究选择Compass力场, 采用正则系综, 利用Ewald方法计算库仑力和范德华力, 每次模拟的总时间为1 000 ps, 计算步长1 ps, 温度设定为298 K, 具体参数如表 2所示.
表 2(Table 2)

表 2 参数设置Table 2 Parameter setting 力场 电荷指派 求和方法 迭代步数 Compass Forcefield assigned Group based 5 000

表 2 参数设置 Table 2 Parameter setting

1.2 溶剂效应对石膏晶体形态影响的分子动力学模拟石膏晶体稳定面参数如表 3所示, 可知[(020), (021), (11-1)和(110)]晶面为石膏晶体主要的四组晶面, 故本文选取此四组晶面作为主要研究对象, 来探究石膏晶体主要晶面在不同质量分数氯化钠水溶液下的生长变化规律.由表 3可知, (020)晶面、(11-1)晶面和(021)晶面多重度分别为2, 4, 4, 占晶体总面积的比例分别达到了87.37%, 6.69%和5.63%, 其中(020)晶面占比最大；通过Growth Morphology模块计算可得到多重度为2的(020)晶面附着能为-21.38 kJ·mol^-1, 同理可得到(11-1)晶面、(021)晶面和(110)晶面的附着能分别为-332.40, -279.50和-536.68 kJ·mol^-1.经过实践证明, 此模块的计算方法不仅可以考虑晶体表面能的影响, 其计算精度也较高.因(020)晶面面心距离最小, 附着能最低, 占晶体总面积比例最大, 故认为(020)晶面在自然环境条件下是石膏晶体生长最慢的晶面, 而(110)晶面占总面积比例最小, 晶面面心距离最大, 说明(110)晶面生长速度明显快于其他三个晶面, 从而决定了石膏晶体在自然环境条件下的形貌为平板型晶体结构.
表 3(Table 3)

表 3 石膏晶体稳定面参数Table 3 Gypsum crystal stability surface parameters 晶面米勒指数(hkl) 多重度 面心距离/nm 单个晶面占晶体总面积的比例/% Eatt/(kJ·mol-1) (020) 2 2.138 87.37 -21.38 (11-1) 4 27.950 6.69 -332.40 (021) 4 33.240 5.63 -279.50 (110) 4 53.668 0.31 -536.68

表 3 石膏晶体稳定面参数 Table 3 Gypsum crystal stability surface parameters

晶体表面能和晶体生长/溶解的动力学密切相关, 通过控制晶面的表面能便可调控动力学速度以及晶体的形貌.将界面现象与动力学的研究方法相结合, 为研究提供了新思路, 有望揭示对晶体界面性质以及动力学的新见解^[15].
石膏晶体表面能是确定晶体表面稳定性和吸附性的重要物理量，如表 4所示.晶面单位面积上自由能的增加称为晶体的表面能(J/m²):

(1)

表 4(Table 4)

表 4 石膏晶体主要晶面表面能Table 4 Gypsum surface energy of main crystal surface 晶面 Eslab/eV Ebulk/eV N a/nm b/nm Esurf/(J·m-2) (020) -7 935.661 4 -15 874.143 5 4 0.630 11 0.569 29 0.629 (021) -7 936.724 2 -15 874.143 5 4 0.630 11 0.813 66 0.108 (11-1) -7 935.951 8 -15 874.143 5 4 0.955 32 0.569 29 0.329 (110) -7 935.599 5 -15 874.143 5 4 0.654 40 0.813 66 0.442

表 4 石膏晶体主要晶面表面能 Table 4 Gypsum surface energy of main crystal surface

式中: E_slab为所截晶面弛豫后两个表面的总能量；A为优化后的表面积；N为表面模型中的原子个数；E_bulk为晶体内部原子数目相同时系统的能量值.
因为CaSO₄·2H₂O在晶面上存在构象和分子间距差异, 而氯化钠的加入可导致溶液与晶面间形成势垒差异, 引起晶面间的相对生长速率发生改变, 进而影响到石膏晶体的生长习性.设定计算步长为1 ps, 静电力和范德华力通过Ewald方法进行统计, 模拟计算时长为20 ps, 主要晶面参数如表 5所示.
表 5(Table 5)

表 5 晶面参数及分子数目Table 5 Crystal face parameters and molecular number 晶面 a/nm b/nm c/nm 分子总数 (020) 1.260 1.139 2.459 16 (11-1) 1.260 3.255 2.080 32 (021) 3.821 1.139 2.117 32 (110) 1.309 3.255 2.159 32

表 5 晶面参数及分子数目 Table 5 Crystal face parameters and molecular number

图 4为25 ℃, 2%氯化钠水溶液体系下石膏晶体主要晶面与溶剂层的双层结构模型.棍状模型为溶液, 球棍模型为晶体表面.优化此双层结构模型, 利用Forcite模块模拟计算得出总能量E_tot, 溶剂层能量E_sol, 溶剂与晶面的附着能ΔE按照式(2)计算^[16]：

(2)

图 4(Fig. 4)

图 4 石膏晶体-NaCl溶液双层结构模型Fig.4 The double-layer structure model of gypsum crystal-NaCl solution

1.3 结果与讨论利用分子动力学计算可以获得E_sol, E_tot，以及随着氯化钠溶液质量分数改变，晶面与溶剂层附着能ΔE发生的变化, 分别如图 5~7所示；随氯化钠溶剂质量分数的升高, 单个晶面占总晶面面积比例变化, 如图 8所示.观察图 5和6可知, E_sol, E_tot随氯化钠质量分数的升高, 均呈现增大的趋势, 在本实验条件下, 可将石膏晶体按照不同晶面与氯化钠溶液作用的附着能划分为三个区域.
图 5(Fig. 5)

图 5 溶剂层能量Esol变化值的绝对值Fig.5 The absolute value of solvent layer energy Esol change

图 6(Fig. 6)

图 6 总能量Etot变化值的绝对值Fig.6 The absolute value of total energy Etot change

图 7(Fig. 7)

图 7 溶剂与晶面附着能ΔE变化值的绝对值Fig.7 The absolute value of solvent and crystal plane attachment energy ΔE change

图 8(Fig. 8)

图 8 单个晶面占总晶面面积的比例Fig.8 Ratio of single crystal plane to total surface area

从图 7可知, 石膏各晶面附着能所占比例在氯化钠溶液质量分数较低时, 与自然环境下的接近, 说明氯化钠溶液质量分数越低对石膏晶体的形貌影响作用越小.当氯化钠溶液质量分数升高时, (020), (11-1)面占总面积的比例也逐渐增大, 而(021), (110)面占总面积的比例变化却截然相反, 导致石膏晶体的形貌向纵向发展.当氯化钠溶液质量分数>8%时, (021)晶面的附着能与其他三个晶面差距过大, 导致其溶解速率最快, 最终(021)晶面消失.
2 实验2.1 实验材料与方法实验设定氯化钠溶液质量分数为 < 1%, 2%, 6%, 10%, 采用硫酸钠和氯化钙完全混合作为CaSO₄·2H₂O及NaCl的来源, 质量分数不够的溶液添加氯化钠试剂补充.确定CaSO₄·2H₂O与氯化钠溶液的固液比为1∶10, 放置在烧杯中, 设定搅拌转速为100 r/min, 温度常温, 磁力搅拌12 h.对反应后的溶液进行固液分离, 对所得固体用蒸馏水洗净后, 在30 ℃的烘箱内进行烘干, 用扫描电镜观察不同质量分数氯化钠溶液体系下CaSO₄·2H₂O晶体形貌变化.
2.2 结果与讨论图 9为自然环境条件下石膏晶体形貌图, 呈片状平板型晶体结构.
图 9(Fig. 9)

图 9 自然环境下石膏晶体形貌Fig.9 Crystal morphology of gypsum in natural environment

图 10所示为不同质量分数氯化钠溶液体系下CaSO₄·2H₂O晶体形貌.从图中可以看出, 低质量分数氯化钠溶液体系下, 二水硫酸钙晶体形貌为片状平板形, 接近自然状态下石膏晶体形貌；随着质量分数的升高, 二水硫酸钙晶体某些晶面的生长受到抑制, 形貌由平板状变为柱状, 纵向发展, 与模拟预测结果一致.
图 10(Fig. 10)

图 10 不同质量分数氯化钠溶液条件下CaSO4·2H2O晶体形貌Fig.10 Crystal morphology of CaSO4·2H2O in different concentrations of sodium chloride solution (a)—未添加NaCl, 氯化钠质量分数 < 1%；(b)—氯化钠质量分数为2%；(c)—氯化钠质量分数为6%；(d)—氯化钠质量分数为10%.

3 结论1) 通过计算晶体稳定面参数, 可知在自然环境体系下石膏晶体主要是由(020), (021), (11-1)和(110)晶面组成, 根据各个晶面的面心距离、附着能以及占晶体总面积比例, 可知石膏晶体在自然环境下的形貌为片状平板型结构.
2) 随着氯化钠溶液质量分数的升高, 溶剂表面能、溶剂与晶面的结合能和附着能均呈现增加的趋势, 当氯化钠溶液质量分数较低时, 石膏晶体与自然环境下的晶体形貌相差不大；随着氯化钠质量分数的升高, (020), (11-1)面占总面积的比例增大, 而(021), (110)面占总面积的比例减小, 导致石膏晶体的形态纵向变大；当氯化钠溶液质量分数>8%时, 因(021)面的附着能远高于其他三个晶面, 最终导致石膏晶体中(021)晶面消失.
3) 通过实验可知, 低质量分数氯化钠溶液对石膏晶体形貌影响不大, 晶体形貌接近于自然环境条件下.随着氯化钠溶液质量分数的升高, 石膏晶体形貌由原来的平板状转变为柱状, 纵向发展, 与模拟得到的结果一致, 体现了各个晶面能量变化对石膏晶体形貌的影响.说明利用分子动力学模拟预测石膏晶体生长规律的方法是可行的.
参考文献

[1]	Liu F G, Cai M Y, Liu X L, et al. O3 oxidation combined with semi-dry method for simultaneous desulfurization and denitrification of sintering/pelletizing flue gas[J]. Journal of Environmental Sciences, 2021, 104(6): 253-263.
[2]	Tang L, Xue X D, Jia M, et al. Iron and steel industry emissions and contribution to the air quality in China[J]. Atmospheric Environment, 2020, 237(8): 2981-2994.
[3]	Bo X, Li Z, Qu J, et al. The spatial-temporal pattern of sintered flue gas emissions in iron and steel enterprises of China[J]. Journal of Cleaner Production, 2020, 266: 121667-121710. DOI:10.1016/j.jclepro.2020.121667
[4]	Zhu J, Ye S C, Bai J, et al. A concise algorithm for calculating absorption height in spray tower for wet limestone-gypsum flue gas desulfurization[J]. Fuel Processing Technology, 2015, 129: 15-23. DOI:10.1016/j.fuproc.2014.07.002
[5]	Melliti E, Touati K, Bruggen B, et al. Effect of Fe2+ ions on gypsum precipitation during bulk crystallization of reverse osmosis concentrates[J]. Chemosphere, 2021, 263: 127866-12877. DOI:10.1016/j.chemosphere.2020.127866
[6]	Koralegedara N H, Pinto P X, Dionysiou D D, et al. Recent advances in flue gas desulfurization gypsum processes and applications-a review[J]. Journal of Environmental Management, 2019, 251: 109572-109585. DOI:10.1016/j.jenvman.2019.109572
[7]	Cordoba P. Status of flue gas desulphurisation (FGD) systems from coal-fired power plants: overview of the physic-chemical control processes of wet limestone FGDs[J]. Fuel, 2015, 144(15): 274-286.
[8]	史培阳, 邓志银, 袁义义, 等. 利用脱硫石膏水热合成硫酸钙晶须[J]. 东北大学学报(自然科学版), 2010, 31(1): 76-79. (Shi Pei-yang, Deng Zhi-yin, Yuan Yi-yi, et al. Hydrothermal synthesis of calcium sulfate whiskers from desulfurized gypsum[J]. Journal of Northeastern University(Natural Science), 2010, 31(1): 76-79.)
[9]	Khorshidi Z N, Tan X, Liu Q, et al. Effect of inorganic salt contaminants on the dissolution of kaolinite basal surfaces in alkali media: a molecular dynamics study[J]. The Journal of Physical Chemistry C Nanomaterials & Interfaces, 2018, 122(9): 4937-4944.
[10]	Khorshidi Z N, Tan X, Liu Q, et al. Molecular dynamics study of the dissolution mechanism of kaolinite basal surfaces in alkali media[J]. Applied Clay Science, 2018, 152(2): 29-37.
[11]	Cui S W, Wei J A, Liu W W, et al. Approximate expression of Young's equation and molecular dynamics simulation for its applicability[J]. Chinese Physics B, 2019, 28(1): 527-531.
[12]	Chen L Z, She C C, Pan H X, et al. Habit prediction of 3, 4, 5-trinitro-1 H-pyrazole in four solvent mediums using a molecular dynamics simulation[J]. Journal of Crystal Growth, 2018, 7: 58-64.
[13]	Chen H, Duan S, Sun Y, et al. Molecular dynamics simulations of solvent effects on the crystal morphology of lithium carbonate[J]. RSC Advances, 2020, 10(10): 5604-5609.
[14]	Sharma S. Molecular dynamics simulation of nanocomposites using BIOVIA Materials studio, Lammps and Gromacs[M]. New Delhi: Elsevier, 2019.
[15]	Ren X T, Ye D Y, Ding N, et al. A molecular dynamics simulation of solvent effects on the crystal morphology of FOX-7[J]. Acta Armamentarii, 2015, 36(2): 272-278.
[16]	Yoreo J D, Gilbert P U, Sommerdijk N A, et al. Crystallization by particle attachment in synthetic, biogenic, and geologic environments[J]. Science, 2015, 349(6247): 498-508.