外电场作用下水化聚全氟磺酸钾膜中水分子的电渗迁移

朱素华 严六明 纪晓波 邵长乐 陆文聪

(1上海大学材料科学与工程学院,上海 200444; 2上海大学理学院化学系,上海 200444;3常熟理工学院物理与电子工程学院,江苏常熟 215500)

外电场作用下水化聚全氟磺酸钾膜中水分子的电渗迁移

朱素华1,2,3严六明2,*纪晓波2邵长乐2陆文聪2

(1上海大学材料科学与工程学院,上海 200444;2上海大学理学院化学系,上海 200444;3常熟理工学院物理与电子工程学院,江苏常熟 215500)

用分子动力学模拟方法研究外电场(简称电场)作用下水化聚全氟磺酸钾膜中水分子的电渗迁移运动,并分析探讨膜的结构与水分子的电渗迁移特性的关系.结果表明,无外加电场时水分子和K+的速度都服从麦克斯韦分布;施加适当电场时,水分子和K+在垂直电场方向上的速度分量仍服从麦克斯韦分布,但平行电场方向上的速度分量则服从峰值漂移的麦克斯韦分布.并且,峰值漂移速度可作为水分子和K+的平均迁移速度的近似值,从而计算得到水分子的电渗系数.结果还显示,K+第一配位层内平均含有约4.04个水分子,它们的平均迁移速度只有K+的57%.这部分水分子贡献的电渗迁移系数为总电渗迁移系数(2.97)的77%.

分子动力学模拟; 质子交换膜; 外电场; 电渗迁移; 速度分布函数

由于极佳的质子电导率、优越的化学和电化学稳定性、良好的机械性能,聚全氟磺酸(PFSA)是目前质子交换膜燃料电池(PEMFC)中最常用的膜材料.但PFSA作为PEMFC质子交换膜也有严重缺点[1]:一方面,水化PFSA膜在温度高于80℃时蒸发失水,严重影响膜的水化和质子的自由迁移;另一方面, PEMFC低温工作时,阴极区积聚了大量的反应水和从阳极电渗迁移过来的水,淹没电极催化剂,妨碍反应原料气体向电极催化剂表面扩散,影响电极反应的顺利进行.因此,研究水化PFSA膜中水分子的迁移,具有重要的实际意义.

PFSA质子交换膜由憎水的氟碳主链和亲水的醚氧磺酸侧链组成.充分水化的PFSA膜中,憎水的氟碳主链聚集在一起,形成憎水区域;亲水的磺酸基团和水分子聚集在一起形成亲水区域.亲水区域之间由直径约1.0 nm的通道相连,形成一个连续的亲水相[2].在亲水相,阳离子可以自由移动,质子更可以通过Grotthuss机制迁移[3],因此,PFSA膜的质子电导率特别高.同时,阳离子迁移时,也会携带一定量的中性水分子一起迁移,引起电渗迁移(electroosmotic drag).

电渗迁移只是水分子在质子交换膜中同时存在的多种迁移过程的一种,并且,分别测量这些迁移非常困难.同时,由于PFSA膜的结构复杂,要从分子水平上研究其中的阳离子对中性水分子的电渗拖动更是困难.现有的理论模型,如基于Nernst-Planck方程的统计力学模型[4]、基于Navier-Stokes方程的连续模型[5]等,只能从统计角度或宏观角度解释电渗迁移的某些方面的特征;而要在分子水平上建立完整的电渗迁移理论模型仍然是个挑战.另一方面,分子动力学模拟(MD)作为一种研究复杂分子体系各种结构特征和动力学性质的重要方法[6-8],也已被应用于研究电渗现象[9-10],计算各种条件下的电渗系数[11-13].如Freund[9]通过对MD模拟得到在施加151 V·μm-1的电场时系统中水分子和Cl-的位置、速度分别进行加和,得到水分子和Cl-的粒子数分布和密度分布,以及平均流速度分布.但是,在燃料电池等实际体系中,水分子的平均迁移速度通常不到1 m· s-1,而其平均热运动速率却可达642 m·s-1(25℃);因此,MD模拟方法通过对速度进行加和计算平均迁移速度往往难以收敛,有时即使收敛其可靠性也不高.正因如此,为了加快收敛速度,往往需要在MD模拟中施加100 V·μm-1或更高的电场,远远偏离实际情况[9-10].为了克服这些困难,Yan等[11]提出先从MD模拟结果统计得到外场作用下水分子和水合氢离子的速度分布函数,再计算平均迁移速度和电渗迁移系数(即每个阳离子迁移时所携带的平均水分子数)的方法.利用这种方法,即使外电场强度低至2 V·μm-1,也能得到很好的收敛效果.

目前,与实际应用关系最密切的是质子对水分子等中性分子的电渗拖动现象.K+的半径与水合氢离子的相当,而它们与水分子之间的相互作用特征却截然不同,如K+与水分子之间主要是静电相互作用,而水合氢离子与水分子之间虽也存在静电相互作用,但主要是氢键相互作用.因此,研究水合氢离子被K+取代后水化聚全氟磺酸钾膜(K-PFSA)的电渗现象,可以了解电渗现象与离子和水分子之间相互作用的关系,具有重要的理论意义.本文中,我们将应用MD方法研究水化K-PFSA膜的微观结构、磺酸基和K+周围水分子的分布;同时,也将采用Yan等提出的通过统计速度分布函数计算平均迁移速度的方法,得到水化K-PFSA膜中水分子的电渗系数,并探讨其电渗机理.

1 研究方法

模拟系统包含12个全氟磺酸寡聚体,48个K+, 384个水分子(图1).每个全氟磺酸寡聚体由1条64个碳原子组成的氟碳主链和4条醚氧磺酸侧链组成.侧链上的所有磺酸已经电离,生成磺酸基.初始模拟原胞按下列方法生成:首先在边长为10.0 nm的正方体盒子中逐一放入12个PFSA寡聚体离子,然后在盒子的空隙中依次随机插入384个水分子和48个K+.如此生成的体系结构松散,密度只有0.1 g·cm-3左右.在MD模拟过程中先以高压1.01×108Pa压缩,然后再在常压1.01×105Pa下模拟直至达到平衡密度.

模拟体系的总势能由成键相互作用和非键相互作用组成:

其中,前三项表示成键相互作用,分别是键的伸缩(stretching)能、键角弯曲(bending)能和二面角扭曲(torsion)能,kb、kθ和Vm为对应的力常数,b0、θ0是平衡键长和键角,m和δ是二面角旋转多重度和平衡相因子.后三项表示非键相互作用,分别是Lennard-Jones势、静电势和与外电场作用的静电能,其中εij为i、j原子间Lennard-Jones势能阱的深度,σij是i、j原子间Lennard-Jones势能为零时的距离,rij是i、j原子间距离,qi为i原子的电荷,ψ(ri)为ri处的外电场强度.本文中,PFSA分子采用OPLS全原子模型[14],并补充缺少的参数[15];水分子采用刚性TIP3P力场模型[16];K+采用最近被拟合并优化的参数[17]; Lennard-Jones势的交叉项参数由Lorentz-Berthelot混合规则得到.为了便于描述,文中用O2、O1和OW分别表示磺酸根—SO-3氧原子、醚氧原子、水分子的氧原子,以示区别(图1).

本文MD模拟使用DL_POLY 2.18程序[18],周期性边界条件,Ewald算法计算静电势,截断半径取0.95 nm,Verlet蛙跳法积分运动方程[19],时间步长1 fs,常温常压(NPT)系综(Nosé-Hoover算法[20-21],驰豫时间都为0.1 ps).模拟时,先在1 K和1.01×108Pa下模拟100万步,使松散的初始构型逐渐收缩至密度为1.83 g·cm-3左右;然后,温度升至300 K,压力降至1.01×105Pa,再模拟1000万步使系统达到平衡,密度为1.80 g·cm-3左右;最后,在x方向上分别施加0,1,2,4,5 V·μm-1的电场,分别模拟100万步.在最后的100万步模拟中,每10步记录一个构型,以统计分析体系的结构特征和迁移性质.

2 结果与讨论

为了用MD方法统计电场作用下中性水分子和K+的迁移速度,模拟时必须施加适当的电场.电场太小,迁移速度过小统计结果无法收敛;电场太大,对体系的扰动过大而影响模拟结果.为此,体系被分别施加0、1、2、4、5 V·μm-1的电场进行模拟.结果发现,未施加外电场时,体系形成明显的亲水相和憎水相,与文献报道[22]一致.憎水相主要由PFSA的氟碳主链组成,形成具有较大空穴和连通管道的微观结构[23].憎水相的空穴和连通管道中充满亲水的磺酸基团、水分子和K+,形成亲水相.由于亲水相的连通性,水分子和K+可以长距离迁移.模拟发现,施加2 V·μm-1或以下的电场时,系统的微观结构没有明显变化(比较图2(A)和2(B));但施加4 V·μm-1或以上的电场时,PFSA的氟碳主链被沿电场方向拉伸(即图2(C)中的箭头方向),系统中亲水性分子或离子与憎水性基团相互混杂,引起憎水相和亲水相相互分离的结构被破坏(图2(C)).为了获得最大的分子、离子迁移速度,而又不明显影响系统的微观结构,在后面的MD模拟中施加的电场强度不大于2 V·μm-1.

模拟结果发现,施加的电场强度不大于2 V· μm-1时,体系的径向分布函数(RDF)与电场的大小没有显著相关.图3是K+和H2O与其它原子间的RDF.结果显示只有在K+与H2O之间,以及K+和H2O与磺酸基的O2和S原子之间的RDF出现明显的结构,有关特征参数总结于表1.从中可以发现, K+与OW和O2之间有较强相互结合,K+-OW和K+-O2的RDF的第一峰分别出现在0.2675和0.2725 nm处,对应K+和OW(或O2)的Lennard-Jones碰撞参数之和.OW-OW的RDF的第一峰出现在0.2775 nm,谷的位置在0.4525 nm,较相应纯水的峰的位置(0.278 nm)和谷的位置(0.353 nm)[16]略大,说明亲水相中水的结构较纯水松散.O2-OW间的RDF第一峰和谷的位置分别在0.2875和0.4175 nm,排斥力较水分子间更强.结果也显示,H2O和K+远离构成PFSA骨架的C和F,除F-OW在0.3175 nm处出现一个很低的峰以外,其它没有峰值出现,也说明体系中存在憎水相和亲水相的分离.

表1 径向分布函数特征参数Table 1 Characteristic parameters of the radial distribution functions

我们把亲水相划分为三个区域:即磺酸基的第一配位层(I区),K+的第一配位层(II区),以及磺酸基和K+的第一配位层外(III区,也称为水的体相区).磺酸基和K+的第一配位半径取相应S-OW和K+-OW的RDF谷的位置,分别为0.5275和0.3525 nm.表2为I、II、III三个区域中水分子的平均数目及其在电场方向的平均迁移速度.没有施加电场时,所有的平均迁移速度在模拟的统计误差范围内近似为0.当施加1 V·μm-1的电场时,迁移速度很小,无法得到有效的统计结果.当施加2 V·μm-1的电场时, H2O和K+有明显的迁移,可以得到较可靠的迁移速度.在I区中,磺酸基略沿电场相反方向运动,速度-0.53 m·s-1,H2O沿电场方向迁移,速度为1.28 m·s-1.在II区中,K+和H2O的迁移速度分别为8.63和4.91 m·s-1,水分子的迁移速度约为K+的57%.在III区中,H2O的迁移速度为1.80 m·s-1,只比I区的略大.

表2 不同外电场(E)作用时水分子数(nW)和水分子、—SO-3、K+平均迁移速度(vW、vS、vK)Table 2 Number of water molecules(nW)and average transport velocities(vW,vS,vK)for water molecules,—SO-3,and K+,respectively,under various electric fields(E)

简单气体服从麦克斯韦速度分布,甚至连密度高达一般固体30倍的Lennard-Jones凝聚相Ar也服从麦克斯韦速度分布[24].复杂的水化PFSA体系,其中的小分子H2O、CH3OH或带电离子Na+、H3O+等也符合麦克斯韦速度分布[11-12,25].因此,我们用麦克斯韦方程拟合H2O和K+的速度分布函数:

其中,m和v分别为分子(离子)的质量和平动速度,k是玻尔兹曼常数,T是温度参数.施加电场时,在平行于电场方向的速度分布发生了峰值漂移,所以平行于电场方向的速度分量用峰值漂移的麦克斯韦方程拟合:

其中,v//是平行于电场方向的速度分量,v0是电场引起的峰值漂移速度即平均迁移速度,T//是相应的温度参数.在垂直于电场方向的速度分布没有发生峰值漂移,所以垂直于电场方向的速度分量仍用麦克斯韦方程拟合:

其中,v⊥是垂直于电场方向的速度分量,T⊥是相应的温度参数.

从图4可以看出,模拟得到的H2O和K+的速度符合麦克斯韦速率分布(式(2)),而平行于电场方向的速度分量符合峰值漂移的麦克斯韦分布(式(3)),垂直于电场方向的速度分量符合麦克斯韦速度分布(式(4)).并且,在没有电场或电场强度为1 V·μm-1时,拟合得到的各温度参数在误差范围内相互一致(表3).值得注意的是,当施加的电场强度加大到2 V·μm-1时,各温度系数明显高于没有电场或电场强度为1 V·μm-1的情形,这与电场的加热作用有关,与有关的实验观察一致[26].同时我们也可发现,无论那种情况,拟合得到的温度参数都比模拟时设定的平衡温度300 K略高,这与一般NPT系综MD模拟过程不完全遵守能量均分定理,即MD模拟时系统的能量不在各种不同的运动模式间平均分布有关.

利用H2O和K+在电场方向上的速度分量的分布函数,可以得到它们的平均迁移速度,并计算电渗系数[11],有关计算公式罗列如下:

表3 在x方向施加不同外电场时由公式(2-4)拟合得到的温度参数Table 3 Fitted temperature parameters for K+and H2O based on Eqs.(2-4)with various electric fields applied in the x-direction

其中,V是模拟原胞的体积,nW和nK分别是原胞中H2O和K+的数目,vW,0和vK,0分别是电场引起的H2O和K+的平均迁移速度,JW和JK分别为H2O和K+的流量密度,Kdrag是电渗系数.由式(3)和式(5-7)得到的不同电场下H2O和K+的平均迁移速度以及水的电渗系数列于表4.从中可以看出,电场强度为1 V· μm-1时,H2O的平均迁移速度与其统计误差相当((0.11±0.06)m·s-1),结果不可靠.当电场强度增大到2 V·μm-1时,H2O和 K+的平均迁移速度分别为(2.98±0.06)和(8.04±0.10)m·s-1,可得到可靠统计结果.由这些数据计算得到的水的电渗系数为2.97.

比较表2中直接平均法得到的K+的迁移速度,以及表4中峰值漂移法得到的K+的迁移速度,两者结果相互一致.根据I、II和III区中的H2O的数目(表2),如果II区(即K+的第一配位层)中的H2O全都随K+一起运动,则H2O的电渗系数(即每个K+平均可携带的水分子数)可以达到4.04-4.21(表2).但由于K+与H2O的相互作用较弱,如在外电场强度为2 V·μm-1时,K+的迁移速度为8.63 m·s-1,而K+第一配位层中的H2O(II区)的迁移速度只有4.91 m·s-1,其迁移速度约为K+的57%,跟不上K+(表2),根据K+第一配位层内平均有4.04个H2O(表2),可以得到K+第一配位层内的H2O贡献的电渗系数约为2.30,占总电渗系数的77%.另外,I区和III区中的水分子也略被拖动,发生电渗迁移,这两个区域H2O贡献的电渗系数约为0.67,占总电渗系数的23%.综合这些因素总电渗系数为2.97.而H3O+与H2O间的相互作用主要是较强的氢键相互作用,当外电场为2 V·μm-1时,每个H3O+可平均带动6.84个H2O分子一起迁移[11].

表4 根据峰值偏移麦克斯韦速度分布函数得到的不同外电场下的H2O和K+平均迁移速度(vW,0,vK,0)和水分子的电渗系数(Kdrag)Table 4 Average transport velocities of H2O and K+(vW,0,vK,0)evaluated from the peak shifted Maxwellian velocity distribution functions and electroosmotic drag coefficients(Kdrag)of water molecules under various electric fields

从动量传递过程分析:K+首先被电场加速,获得动量;然后,K+与周围(II区)水分子发生碰撞又把动量传递给水分子,使II区中的水分子较I区或III区中的水分子具有更高的迁移速度;这些具有较高迁移速度的水分子又与其它水分子(I区或III区)碰撞再把动量传递给这些水分子,引起更多的水分子的迁移.与H3O+或Na+不同的是,K+与水分子的相互作用力比较弱,K+第一配位层中水分子的迁移速度较K+本身的小许多,如在外电场强度为2 V·μm-1时水分子的迁移速度约为K+的57%,因此,K-PFSA膜中水分子的电渗系数比K+第一配位数小许多,也比PFSA膜和Na-PFSA膜中水分子的电渗系数小很多.比如,K+与H3O+的尺寸相当,但由于H3O+和周围水分子之间形成稳定的氢键网络,所以一个H3O+在电场中迁移时可以拖动6.25-6.84个水分子运动(不考虑质子的hopping运动)[11].再比如,K+和Na+与周围水分子的相互作用类型相同,但K+比Na+尺寸大,所以Na+对水分子的吸引力较K+强,因此,水分子在Na-PFSA膜中的电渗系数可达6.1-8.7[13].由此可知,电渗系数与阳离子的性质之间的关系十分密切,即使相同的体系,阳离子的性质稍有变化,水分子的电渗系数可以有很大的变化.

3 结 论

用MD方法模拟水化K-PFSA膜中H2O和K+在电场作用下的迁移,统计得到相应的速度分布函数;在此基础上,拟合得到H2O和K+平均迁移速度,并计算H2O的电渗系数.结果表明,模拟时施加4 V·μm-1或以上的电场时,K-PFSA膜的微观结构遭到破坏,结果没有代表性;施加1 V·μm-1的电场时,H2O和K+平均迁移速度太小,无法正确统计.施加2 V·μm-1的电场比较合适,既能较好地统计各种平均迁移速度,又没有破坏水化K-PFSA膜的微观结构,与文献报道一致.

模拟结果还表明,不同区域中H2O的平均迁移速度并不相同:其中,K+第一配位层内H2O迁移最快,体相中的H2O次之,磺酸根周围的H2O最慢.由于K+对H2O的相互作用较弱,即使K+第一配位层内的H2O也跟不上K+一起迁移运动,其速度只有K+的 57%左右.根据K+第一配位层内有 4.04个H2O,可以得到K+第一配位层内的H2O贡献的电渗系数约为2.30,约占总电渗系数的77%.而其它区域的H2O贡献的电渗系数约为0.67,约占总电渗系数的23%.通过与水化Na-PFSA中水的电渗系数比较可以发现,虽然K+和Na+的性质相似,但Na+对H2O的相互作用较K+的强,相应的电渗系数也较KPFSA体系大.由此进一步验证:阳离子对H2O的相互作用的强弱,决定了其对H2O电渗拖动作用的大小.

1 Paddison,S.J.Annu.Rev.Mater.Res.,2003,33(1):289

2 Gierke,T.D.;Munn,G.E.;Wilson,F.C.J.Polym.Sci.Polym. Phys.Ed.,1981,19(11):1687

3 Agmon,N.Chem.Phys.Lett.,1995,244(5-6):456

4 Schaetzel,P.;Nguyen,Q.T.;Riffault,B.J.Membr.Sci.,2004, 240(1-2):25

5 Nilson,R.H.;Griffiths,S.K.J.Chem.Phys.,2006,125(16): 164510

6 Zhao,X.;Huang,X.R.;Sun,C.C.Chem.J.Chin.Univ.,2008,29 (4):819 [赵 熹,黄旭日,孙家锺.高等学校化学学报,2008, 29(4):819]

7 Gao,T.H.;Liu,R.S.;Zhou,L.L.;Tian,Z.A.;Xie,Q.Acta Phys.-Chim.Sin.,2009,25(10):2093 [高廷红,刘让苏,周丽丽,田泽安,谢 泉.物理化学学报,2009,25(10):2093]

8 Tao,C.G.;Feng,H.J.;Zhou,J.;Lü,L.H.;Lu,X.H.Acta Phys.-Chim.Sin.,2009,25(7):1373 [陶长贵,冯海军,周 健,吕玲红,陆小华.物理化学学报,2009,25(7):1373]

9 Freund,J.B.J.Chem.Phys.,2002,116(5):2194

10 Qiao,R.;Aluru,N.R.Phys.Rev.Lett.,2004,92(19):198301

11 Yan,L.;Ji,X.;Lu,W.J.Phys.Chem.B,2008,112(18):5602

12 Ji,X.;Yan,L.;Zhu,S.;Zhang,L.;Lu,W.J.Phys.Chem.B,2008, 112(49):15616

13 Yan,L.;Shao,C.;Ji,X.J.Comput.Chem.,2009,30(9):1361

14 Jorgensen,W.L.;Maxwell,D.S.;Tirado-Rives,J.J.Am.Chem. Soc.,1996,118(45):11225

15 Yan,L.;Zhu,S.;Ji,X.;Lu,W.J.Phys.Chem.B,2007,111(23): 6357

16 Jorgensen,W.L.;Chandrasekhar,J.;Madura,J.D.;Impey,R.W.; Klein,M.L.J.Chem.Phys.,1983,79(2):926

17 Lenart,P.J.;Jusufi,A.;Panagiotopoulos,A.Z.J.Chem.Phys., 2007,126(4):044509

18 Smith,W.;Forester,T.R.J.Mol.Graph.,1996,14(3):136

19 Verlet,L.Phys.Rev.,1967,159(1):98

20 Nosé,S.Mol.Phys.,1984,52(2):255

21 Hoover,W.G.Phys.Rev.A,1985,31(3):1695

22 Blake,N.P.;Petersen,M.K.;Voth,G.A.;Metiu,H.J.Phys. Chem.B,2005,109(51):24244

23 Kreuer,K.D.J.Membr.Sci.,2001,185(1):29

24 Mohazzabi,P.;Helvey,S.L.;McCumber,J.Physica A,2002,316 (1-4):314

25 Arias,J.L.;Ruiz,M.A.;Gallardo,V.;Delgado,A.V.J.Control. Release,2008,125(1):50

26 Fang,C.;Wu,B.;Zhou,X.Electrophoresis,2004,25(2):375

Electroosmotic Drag of Water in Hydrated Potassium Perfluorosulfonated Polymer Membrane in External Electric Fields

ZHU Su-Hua1,2,3YAN Liu-Ming2,*JI Xiao-Bo2SHAO Chang-Le2LU Wen-Cong2
(1College of Material Science and Engineering,Shanghai University,Shanghai200444,P.R.China;2Department of Chemistry,College of Science,Shanghai University,Shanghai 200444,P.R.China;3College of Physics and Electronic Engineering,Changshu Institute of Technology,Changshu 215500,Jiangsu Province,P.R.China)

The electroosmotic drag and the corresponding mechanism of water molecules in hydrated potassium perfluorosulfonate electrolyte polymer membrane were studied using molecular dynamics simulations,and the relationship between the membrane structure and electroosmotic drag characteristics was analyzed.It is concluded that velocities of both H2O and K+obey the Maxwell velocity distribution function without external electric field applied.If an appropriate electric field is applied,the velocities of H2O and K+still obey the Maxwell velocity distribution in the direction perpendicular to the electric field,and obey the peak shifted Maxwell velocity distribution in the direction parallel to the electric field.The peak shifting velocities coincide with the average transport velocities of H2O and K+induced by the applied electric field,and could be applied to evaluate the electroosmotic drag coefficient of water.The results also show that the average number of water molecules in the first coordination shell of K+is 4.04,and the average transport velocity of these water molecules is about 57%of that of K+.The electroosmotic drag coefficient contributed by these water molecules is about 77%of total the electroosmotic drag coefficient(2.97).

Molecular dynamics simulation;Proton exchange membrane;External electric field; Electroosmotic drag; Velocity distribution function

O646

Received:April 29,2010;Revised:June 14,2010;Published on Web:July 23,2010.

*Corresponding author.Email:liuming.yan@shu.edu.cn;Tel:+86-21-66132405.

The project was supported by the National Natural Science Foundation of China(20873081),and the Nano project of Shanghai Municipal Science& Technology Commission,China(0952nm01300).

国家自然科学基金(20873081)和上海市科委纳米专项(0952nm01300)资助项目

ⒸEditorial office of Acta Physico-Chimica Sinica