1 引言 (Introduction)

湿法烟气脱硫(WFGD)技术是目前世界上最为有效可靠的脱硫工艺之一(Brogren et al., 1997)，其中，石灰石/石膏湿法脱硫技术应用最为广泛，截止到2011年底，占我国已投产脱硫机组份额的92%以上.2011年新颁布的《火电厂大气污染物排放标准》(GB132323—2011)对我国火电厂SO₂排放提出了更高要求，SO₂排放限值降为100 mg · Nm^-3，重点地区排放限值降为50 mg · Nm^-3.国内现行脱硫工艺面临新的挑战，燃煤锅炉原有脱硫设备大部分将面临升级改造(中国环境保护产业协会脱硫脱硝委员会，2013).因此，进一步提升脱硫效率，开发优化石灰石/石膏湿法烟气深度脱硫技术势在必行.

了解脱硫浆液内组分传质、溶解、反应等机理及深入优化脱硫温度等运行参数是进一步开发深度脱硫技术的关键.近年来，国内外学者针对脱硫过程中SO₂、石灰石等组分吸收溶解机制进行了一系列深入研究(Brogren et al., 1997；Kiil et al., 1998；Wang et al., 2008；Souza et al., 2010)，研究方法主要集中于实验研究、数值计算、工业试验等.由于目前测量手段或数值计算方法的局限性，研究对象的时空尺度主要侧重于微秒及微米以上，针对更小时空尺度的研究相对较少.

分子动力学(Molecular Dynamics，MD)是一种新兴的理论化学计算方法，能够从纳米尺度及纳秒尺度甚至更小的时空尺度深入研究气液/溶液中的分子微观结构、传质扩散机制、反应机理等(Garrett et al., 2006).针对N₂、O₂、CO₂、H₂O、NH₃等小分子气体在水中的溶解机理人们已经运用分子动力学进行了系统的研究(Somasundaram et al., 1999；cha et al., 2004；Paul et al., 2005；Dang et al., 2004).但针对SO₂溶解机制的研究相对较晚，部分先驱者初步研究了SO₂在纯水表面的反射、吸收行为(Matsumoto，1996)，以及水溶液内SO₂周围H₂O朝向、氢键等微观结构、SO₂溶解副产物对SO₂溶解的耦合影响(Tarbuck et al., 2006)等.Shamay等(2011)进一步研究了SO₂在气液界面的迁移途径及分布梯度，Moin等(2011)结合QM-MD方法研究了SO₂在水溶液中的水化结构及SO₂-H₂O作用位等，研究发现，SO₂表面形成了一层稳定的水化层.

前人主要以水溶液作为液相进行研究，针对含CaCO₃等吸收剂及其离子的浆液体系的研究还相对缺乏.从分子层面研究SO₂等组分在浆液中的传质机制及溶解吸收行为能够进一步为钙基湿法脱硫反应机理的深入研究及脱硫工艺优化提供有力的理论依据，具有重要的现实意义.基于此，本文将利用分子动力学计算研究浆液内组分扩散系数及SO₂在浆液中的溶解机制，并进一步分析脱硫运行温度对组分扩散系数及SO₂溶解吸收的影响作用.本文研究结果将为深入了解浆液内SO₂溶解微观结构及筛选脱硫运行温度提供有效理论依据.

2 计算模型及方法(Computational models and simulation method)

本文研究的SO₂-浆液体系主要针对液相内，不考虑气相及未溶解的固相的影响，故溶液元胞模型采用正方体的模拟盒子.由于钙基湿法脱硫吸收剂中的主要作用成分为CaCO₃，又考虑实际浆液环境中HCO^-₃、HSO^-₃等两性离子的复杂反应机制及力场适用性，为简化计算，浆液模型设计中只考虑Ca²⁺和CO^2-₃两种离子组分.石灰石-石膏湿法脱硫浆液中溶解盐的质量分数在1.2%左右，根据浓度换算结果，在模拟体系中加入864个H₂O分子、2个Ca²⁺和2个CO^2-₃来构建所需浆液模型.脱硫工艺稳定运行过程中，浆液单位时间的SO₂吸收量可以通过计算排出的固渣中的含硫物含量得到，以脱硫浆液内含固量约20%(王乃华等，2005)且其中90%为CaSO₄ · H₂O的条件计算，折合SO₂吸收量为约1.5 mol · L^-1.进一步考虑本文SO₂溶解特性的研究初衷，在“浆液”盒子中随机插入16个SO₂分子，构成SO₂吸收浓度近似为1 mol · L^-1的溶液体系.为了消除边界效应，盒子的3个方向都采用了周期性条件.初始盒子的尺寸为3.032 nm×3.032 nm×3.032 nm.本文中，H₂O分子的势能函数采用了著名的SPC模型(Berendsen et al., 1981).Ca²⁺、CO^2-₃、SO₂的分子/离子模型由Material Studio软件(Accelrys Inc.)构建优化，各原子的力场参数及电荷分布由COMPASS力场赋予，该力场已分别被证实适用于这3种粒子的溶液/液相特性计算(Yuan et al., 2013；Aguiar et al., 2013；Yang et al.，2000).

体系构建完成后，先对系统用Smart Minimizer法进行能量最小化处理；随后在等温等压系综(NPT)下将体系进行300~500 K、100 ps的退火弛豫以消除盒子内的局部应力并找到全局能量最低构象，退火弛豫循环5次，升温及降温过程中每次温差变化为10 K.然后对最优构象进行150 ps的NPT平衡弛豫，当系统势能波动在5%以内，计算体系达到稳定.最后对稳定构象进行600 ps的正则系综(NVT)模拟计算，每0.5 ps收集一次各原子的动力学计算全轨迹.为了消除计算运行的初始影响，取后500 ps的数据进行性质分析.NPT系综下运行的平衡压力为1 atm(1.01×10⁵ Pa)，NPT系综和NVT系综下的运行温度均为计算目标温度.

分子动力学运行计算中，选用的力场为COMPASS通用力场，NPT系综弛豫及NVT系综模拟中温度控制采用Andersen热浴法(Andersen，1980)，压力控制采用Bererdsen恒压法(Berendsen et al., 1984).计算非键作用中，范德华作用(van der Waals)采用Atom based方法(Karasawa et al., 1992)，截断半径为1.25 nm，静电作用采用Ewald方法(Allen et al., 1989).各分子的初始运动速度采用Maxwell-Boltzman分布随机取样，采用Velocity-verlet积分方法进行求解，时间步长为1fs.所有计算模拟过程均在Material Studio软件(Accelrys Inc.)上进行.

石灰石-石膏湿法脱硫工艺的运行温度在50 ℃左右，因此，本文设计的模拟温度为293.15~343.15 K，每间隔10 K计算1次.

3 结果与分析(Results and discussion)

图 1给出了293.15 K运行温度下的模型例图，其中，线性分子模型为H₂O分子，黄红球棍模型为SO₂分子，灰红球棍模型为CO^2-₃，绿色小球为Ca²⁺.经过不同温度的等温等压系综(NPT)平衡弛豫后，模拟体系密度达到稳定.在计算温度范围内随温度的升高293.15~343.15 K，体系密度缓慢下降(1.005~0.952 g · cm^-2)，与实际情况下水密度接近(CaCO₃溶质的影响可以忽略)，表明计算可靠(图 2).

3.1 溶液组分扩散系数

溶液中组分的扩散系数是表征体系内各物质传质输运特性的重要参数.湿法脱硫反应中，气液传质阻力主要来源于液相(钟毅，2008)，溶液组分扩散系数是影响整个脱硫反应工艺传质能力、脱硫反应速率及脱硫效率等的关键因素.分子动力学中计算粒子的扩散系数可以通过求解利用该粒子的均方位移(MSD，Mean Square Displacement)的Eintein关系式获得，即：

式中，D为扩散系数(m² · s^-1)，N表示该粒子的数目，r_i t 和r_i0 分别表示t时刻和初始时刻第i个粒子的位置向量(m)，<…>表示系综平均.由于计算MSD曲线时，已经对N个粒子进行了统计平均，当模拟时间t够长时，该粒子的扩散系数D可以通过求解MSD曲线对时间t的斜率获得.本研究中，组分扩散系数计算结果既是计算的重要内容，也是检测模型设计及力场选择可靠性的重要手段.

表 1中详细列出了293.15~343.15 K运行温度下H₂O、SO₂、Ca²⁺、CO^2-₃ 4种粒子的扩散系数计算结果.溶液中SO₂、Ca²⁺、CO^2-₃等粒子扩散系数一般由实验测得(丁红蕾，2010；Himmelblau，1964；Ebanks et al., 2010)，通过MD计算的文献较少.计算发现，H₂O分子扩散系数从293.15 K时的3.074×10^-9m² · s^-1逐渐上升到343.15 K时的6.079×10^-9 m² · s^-1.Mark等(2001)对SPC、TIP3P、SPC/E等多个H₂O分子模型及修正模型在298 K下的扩散系数进行了计算，结果分别为4.22×10^-9 m² · s^-1(SPC)、2.75×10^-9 m² · s^-1(SPC/E)、5.67×10^-9 m² · s^-1(TIP3P)，本文计算结果与该文献符合较好.溶液中SO₂扩散系数的计算结果与文献中参考值(20~40 ℃，1.62×10^-9 ~2.59×10^-9 m² · s^-1，Himmelblau，1964；25 ℃，1.83×10^-9 m² · s^-1，丁红蕾，2010)符合较好.而Ca²⁺和CO^2-₃的计算结果则与文献(Chen et al., 2013；Bruneval et al., 2007)的研究结果较为接近.

表 1中4种粒子的扩散系数均随温度的升高而升高，这是由于温度升高，粒子内能增加，动能增大，扩散能力增强.进一步而言，粒子之间碰撞几率加大，有效碰撞机会增多，从而加速脱硫反应.H₂O分子在溶液中是以氢键结合的分子团簇形式存在，温度升高使得H₂O分子在团簇内平衡位置振动增强，氢键结构被破坏，自由H₂O分子增多.这也会增大依赖于H₂O分子扩散的粒子的扩散能力，如H⁺等(Walbran et al., 2001)，进一步加大了脱硫反应速率.

溶液中组分扩散系数与温度的关系应符合扩散动力学的Arrhenius定律，即lnD=lnD_o-ΔE/(RT)，其中，D_o为指前因子(m · s^-1)，ΔE为扩散系数的表观活化能(kg · mol^-1)，R为气体常数，取值为8.314 J · mol^-1，T为温度(K).图 3中给出了H₂O、SO₂、Ca²⁺、CO^2-₃ 4种粒子的扩散系数的对数lnD与温度的倒数1/T的关系曲线.可以看到，4种粒子的lnD-l/T曲线线性较好.根据公式lnD=lnD_o-ΔE/(RT)计算得到了4种粒子的扩散活化能，其中，H₂O分子的扩散活化能为11.27 kJ · mol^-1，SO₂的扩散活化能为11.47 kJ · mol^-1，Ca²⁺和CO₃^2-的扩散活化能分别为16.74 kJ · mol^-1和21.32 kJ · mol^-1.

3.2 径向分布函数及水化分析

溶液内溶质分子或离子与溶剂分子产生相互结合，称为溶剂化现象，水溶液中这种现象又称为水化现象.溶质水化对溶质溶解、溶质扩散等都有重要影响，研究湿法脱硫中SO₂在浆液中水化、溶解的微观机制十分必要.分子动力学中分析溶液内粒子间微观结构的直接有效手段之一是分析径向分布函数(Radical Distribution Function，RDF)(Moin et al., 2011；Kerisit et al., 2005；Bruneval et al., 2007).通过径向分布函数曲线可以获得溶液体系内粒子间相互作用、粒子间有序性等，还可以进一步计算得到粒子间的配位能力，具体计算方法如下(Allen et al., 1989)：

式中，n_AB为原子B在原子A附近距离r的球形范围的配位数，ρ_B是原子B在该范围内的数密度，g_AB(r)即为原子B对原子A的径向分布函数.参照式(2)对溶质-水分子径向分布函数从溶质质心积分到第一峰谷处，即可得到溶质粒子周围的水分子的配位情况.

图 4给出了不同温度下293.15～343.15 K模拟体系内SO₂中S原子与周围的H₂O中O原子的径向分布函数(S(SO₂)-O_w RDF)，并由该径向分布函数进一步计算得到了SO₂在不同温度的浆液内的H₂O配位数变化情况，具体见图 4右上角附图.S(SO₂)-O_w RDF中只有一个明显单峰，第一峰位于约0.4225 nm处，第一峰谷位于0.5625 nm处，表明浆液内SO₂与H₂O之间只形成一个相对明显的水化层，这与Moin等(2011)的研究结果一致.不同温度下293.15～343.15 K第一峰与第一峰谷的位置基本没有变化，即温度对SO₂与H₂O结合的局部结构没有太大影响.峰值大小可以表征SO₂与周围H₂O之间作用的强弱程度，313.15 K与323.15 K下计算峰值最高，分别达到1.372与1.347，而其他4种温度下则峰值相对较低.进一步分析发现，SO₂周围的H₂O配位数在计算温度区间内形成明显单峰分布，313.15 K与323.15 K时H₂O配位数最大，分别达到20.30与19.89，即313.15~323.15 K之间SO₂与H₂O作用最强，最有利于SO₂溶解.

为了进一步研究浆液中Ca²⁺对SO₂的作用机制，图 5中给出了Ca²⁺与S(SO₂)的径向分布函数.RDF曲线中Ca²⁺与SO₂并没有形成稳定的峰型，只存在一个模糊单峰，表明模拟浆液中Ca²⁺与SO₂并没有形成稳定结构，Ca²⁺与SO₂之间的相互作用主要为范德华作用.313.15 K模拟温度下峰值最高，即323.15 K工况下Ca²⁺对SO₂作用最强，最利于吸收反应.这可能与本文之前分析发现的313.15～323.15 K温度下SO₂最易溶解有关.根据湿法脱硫工艺反应原理，SO₂溶解在浆液中并与H₂O反应形成H₂SO₃，H₂SO₃电离形成HSO^-₃、SO^2-₃等离子，Ca²⁺与这些离子相互结合然后进行化学反应.SO₂溶解的难易程度决定了其与Ca²⁺的作用强弱.

3.3 SO₂与浆液的相互作用能

通过引入相互作用能(Interaction Energy，IE)概念(Hanus et al., 2006；廖瑞金等，2012)进一步定量分析了温度对浆液中SO₂扩散、溶解能力的影响机制.依照相互作用能概念，本文定义如下计算式：

式中， E_total为体系总能量(kJ · mol^-1)；E_SO2为单个SO₂分子具有的能量(kJ · mol^-1)，通过单个SO₂在同一温度的真空中进行NVT模拟获得，n为SO₂的个数；E_slurry为纯浆液具有的能量(kJ · mol^-1)，通过扣除SO₂组分的同一浓度、分子数的浆液体系在相同条件的分子动力学模拟下获得；IE即为单个SO₂分子在浆液中受到的平均作用(kJ · mol^-1).进一步分析IE的动能部分，将E_total、E_SO2、E_slurry分别替换成体系总动能、单个SO₂分子动能、纯浆液体系动能，则可以获得浆液对单个SO₂动能的平均影响作用(IE_kinetic).其表征了单个SO₂分子在浆液内受到的捕集作用，IE_kinetic能量越小捕集作用越大，SO₂溶解能力越好.图 6中给出了本文计算的SO₂与所处浆液介质之间的温度关系曲线和IE_kinetic-温度关系曲线.可以看到，SO₂与浆液的相互作用能IE绝对值在计算温度区间293.15~343.15 K随温度的增大而增大.文献(Hanus et al., 2006)中指出，相互作用能为负值表明研究粒子与体系互相吸引结合，且其绝对值越大则结合作用越大.即随温度的升高，浆液与浆液内SO₂的结合能力增大.因为随温度升高，一方面浆液内粒子扩散系数增大，碰撞几率变大，有效碰撞机会增多；另一方面体系内活化SO₂分子增多，利于脱硫反应进行.而IE_kinetic-温度关系曲线则在293.15~343.15 K区间内呈一个明显的凹型分布，先随温度升高而减小，323.15 K左右之后又随温度的增大而增大，即浆液对SO₂分子的捕集作用在323.15 K左右达到最大.这与本文之前分析的293.15~343.15 K温度下SO₂水化作用最明显相吻合.

4 结论(Conclutions)

本文利用分子动力学方法研究了钙基湿法脱硫中浆液组分的扩散行为及SO₂在浆液中的溶解机制，并通过计算不同温度模型293.15~343.15 K，分析获得了温度对脱硫浆液组分扩散能力及SO₂溶解能力的影响规律.所得结论如下：

1)通过不同温度293.15~343.15 K下浆液模型的模拟计算，得到了浆液中H₂O、SO₂、Ca²⁺、CO^2-₃ 4种粒子在计算温度区间内的扩散系数，分别为3.074×10^-9~6.079×10^-9 m² · s^-1(H₂O)、1.649×10^-9~3.198×10^-9 m² · s^-1(SO₂)、 0.3157×10^-9~0.8863×10^-9 m² · s^-1(Ca²⁺)及0.2516×10^-9~0.9674×10^-9 m² · s^-1(CO^2-₃).4种粒子的扩散系数均随温度升高而增大.进一步通过Arrhenius公式计算得到了这4种粒子的扩散活化能，分别为11.27、11.47、16.74及21.32 kJ · mol^-1.

2)浆液内SO₂与H₂O之间只形成一个相对明显的水化层；SO₂周围的H₂O配位数在293.15~343.15 K计算温度下内呈单峰分布；313.15~323.15 K时SO₂水化作用最强，最利于SO₂溶解.

3)进一步的SO₂-浆液相互作用能分析表明，293.15~343.15 K温度内随温度升高，浆液与SO₂结合能力增大；而323.15 K左右浆液对SO₂捕集能力最强.