28水化硅酸钙的分子动力学模拟

水化硅酸钙的分子动力学模拟王渊，张文生，叶家元(中国建筑材料科学研究总院绿色建筑材料国家重点实验室，北京 100024)摘要：以Hamid模型和Bonaccorsi模型为基础建立了初始结构，并采用分子动力学方法，模拟了不同钙硅比（Ca/Si=0.67、0.83、1.0、1.5）的水化硅酸钙的结构。

根据模拟结果，得到了原子间的距离、径向分布函数、配位数、均方位移等参数。

模拟结果表明：无定形态水化硅酸钙存在近程有序远程无序的结构特点；无定形水化硅酸钙的基本结构单元为硅氧四面体，且以Q2形式连接；钙硅比的变化影响了各原子的扩散系数；模拟得到的原子间距离、配位数等结构参数基本与实验值相符合。

关键词：水化硅酸钙；结构；分子动力学1 引言水化硅酸钙（CSH）是硅酸盐水泥的主要水化产物，是决定水泥石性能的关键组分，因此其组成、结构及性质自上世纪50年代Grudemo[1]和Taylor[2]的开创性工作以来一直是水泥科学研究中的重要内容。

各国学者都进行了深入研究，提出了一系列结构模型，如类托贝莫来石和类羟基硅钙石模型[3]、富钙和富硅模型[4]、固溶体模型[5]和纳米结构、中介结构假说[6]。

其中，类托贝莫来石和类羟基硅钙石模型认为，托贝莫来石和羟基硅钙石是无定形水花硅酸钙（CSH）的结构原型，在常温下由化学试剂合成的低钙硅比的C-S-H（I）结构类似于1.4nm托贝莫来石结构，而由纯C3S或纯β-C2S水化得到的高钙硅比的C-S-H （Ⅱ）结构类似于羟基硅钙石结构，只是由于桥[SiO4]4-四面体的缺失，而使得无限长的硅氧四面体链断裂及扭曲，形成了由3n-1个[SiO4]4-四面体构成的短链化合物[3,7,8]。

若n大于1，则[SiO4]4-四面体连接成链状，链中的[SiO4]4-四面体除两端的外其他的都以Q2形式链接。

因此，对托贝莫来石和羟基硅钙石结构的研究，有助于理解水化硅酸钙的物理化学性能、吸附交换机制等性能。

研究材料结构的实验方法通常有X-射线衍射、核磁共振、X-射线光电子能谱等，而对非晶材料而言，这些方法实施起来有一定的困难。

分子动力学(molecular dynamics,MD)模拟作为计算机模拟的一种基本方法，在研究液态和非晶结构方面起着重要的作用。

它根据粒子之间相互作用势，通过对系统运动方程组进行数值积分，得到体系的相轨道，并由此分析系统的各种性质，是联系物质微观信息和宏观性质的一种计算方法。

它可以不受实验条件限制，在任意温度下对无定形体系进行模拟，可分析原子运动轨迹得到径向分布函数(radial distribution function,RDF)、配位数(coordination number,CN)和均方位移等信息。

上世纪90年代中期，Faucon [9]就运用分子动力学方法成功地模拟了水化硅酸钙结构不稳定的原因、硅氧四面体链的断裂机理和阳离子（Al ）取代硅后对结构的影响。

2 模拟过程本文选择了Hamid [10]模型和Bonaccorsi 模型[11]分别作为托贝莫来石（Ca/Si ≤1.0）和羟基硅钙石（Ca/Si=1.5）的结构模型，其中，Hamid 模型的晶体学参数如下：(1)晶系：单斜；(2)晶胞参数：a m =0.669 nm ，b m =0.739nm ，c m =2.277 nm ，γ＝123.49°；(3)空间点群：P 21。

Bonaccorsi 模型的晶体学参数如下：(1)晶系：三斜；(2)晶胞参数：a m =1.0576nm ，b m =0.7265nm ，c m =1.0931nm ，α=101.30°，β=96.98°，γ=109.65°；(3)空间点群：P 1。

模拟尺寸设为4a m ×4b m ×2c m ，各原子对应的离子电荷采用经验电荷（formal charges ），如Si 的电荷取＋4。

模拟时相互作用势采用二体作用势和三体作用势，其中二体势采用BHM （Born-Huggins-Meyer ）势，其函数表达式如式（1），各参数见表1。

二体作用势截断半径取0.8 nm 。

86)exp()(ij ij ij ij ijijij ij r D r C r A r U +--=ρ （1）表1 BHM 作用势参数[12,13] Table1 Parameters for BHM potential [12,13]A ij ×10-11/Jρij /nm Si –Si 1.754 0.029 Si –O 2.516 0.029 O –O 0.372 0.029 Si –Ca 6.729 0.029 O –Ca 9.823 0.029 Ca –Ca 25.582 0.029 Si –H 0.0644 0.029 H –H 0.0318 0.035 O –H0.3720.029三体作用势是一种比van de Waals 力作用势更短程、更弱的作用势，其对总能量的贡献通常弱于1％。

本文所选三体作用势经验函数的表达式如式（2），其参数如表2所示，截断半径取0.345 nm 。

20)cos )(cos exp(),,(θθγγλθ--+-=ijk cikik ikcij ij ijijk ik ij r r r r r r U （2）表2 三体作用势的参数[12]Table 2 Parameters for three body potential [12]λ×10-13/J γ/nm r c /nm θ0 O-Si-O 17.757 0.28 0.30 109.5° Si-O-Si 0.280 0.20 0.26 109.5° H-O-H32.7100.130.16104.5°模拟采用三维周期边界条件(periodic boundary condition,PBC)，对长程相互作用采用Ewald 求和方法，运动方程积分采用“蛙跳”(leap frog)算法，积分步长为1 fs ，速度由Maxwell 分布给出。

为了得到无定形结构中的长程无序、短程有序及硅氧四面体链旋转、扭曲结构，先将体系在5000K 下的NVT 系综中预平衡，然后将体系的温度逐步降低到300K 。

体系在300K 重新达到平衡后，将其输出作为NVE 系综的输入，最终制得稳定的无定形水化硅酸钙，然后对模拟结果进行分析计算，得出水化硅酸钙在300 K 时的结构性质及动力学行为。

3 模拟结果与讨论3.1径向分布函数(RDF)分子动力学模拟无定形体系的结构信息可以用RDF 来分析，表示系统的区域密度与平均密度之比，计算公式如下ρπ24)(r dNr g =（3） 其中ρ为系统密度，dN 表示与中心的距离为r →r+dr 间的原子数目。

图1为总径向分布函数和Ca-O ，Si-O ，Si-Si 原子对的偏径向分布函数。

对比各原子对的偏径向分布函数可知，图1（a ）的第一峰为O-H ，第二峰为Si-O ，第三峰为Ca-O ，第四峰为Si-Si ，各自的位置分别为1.28 Å ，1.56 Å ，2.54 Å和3.10 Å，此即为对应的原子对间的原子距离，与实验值基本符合。

比较不同钙硅比的水化硅酸钙的总径向分布函数，可以看出Ca/Si=0.67和Ca/Si=1.5的比其他两个钙硅比的要多出了一个O-H 峰，这是因为Ca/Si=0.67和Ca/Si=1.5时，其初始结构中硅氧四面体上的O-H 键要明显比其他两个钙硅比的要多。

各偏径向分布函数曲线均存在较明显的第一近邻峰和明显宽化的第二近邻峰，且在较大半径r 范围内其值趋于1，这说明其结构处于近程有序、远程无序的状态，这是明显的无定形态特征。

比较各原子对在不同钙硅比下的偏径向分布函数，各曲线相差不大，这说明虽然输入的初始结构有差别，但是经过模拟达到平衡后，其初始结构的差异对最终稳定结构的影响减小了。

R D Fr/10-10mr/10-10m(a) 总径向分布函数(b) Ca-OR D Fr/10-10mr/10-10mR D F(c) Si-O (d) Si-Si图1 径向分布函数Fig. 1 Radial distribution function3.2 配位数配位数（coordination numbers, CN ）是反映体系结构特征的重要参数，能在一定程度上反映模拟对象的局域结构。

图2为Ca-O ，Si-O ，Si-Si 原子对的配位数分布曲线，曲线的第一个平台对应的纵坐标即为最近邻配位数。

从图中可以看出，Ca-O ，Si-O ，Si-Si 原子对的最近邻配位数分别为6、4和2。

这说明在水化硅酸钙中，Ca 与O 形成了钙氧八面体，Si与O形成了硅氧四面体，而Si与Q2形式2CNr/10-10mCNr/10-10m(a) Ca-O(b) Si-OCNr/10-10m(c) Si-Si图2 配位数Fig. 2 Coordination numbers distribution3.3 均方位移分析粒子运动轨迹随时间变化的统计平均——均方位移（Mean Squared Displacement, MSD），可以得到体系的扩散系数等动态信息。

其基本原理是通过用分子动力学记录一定温度下粒子的动态演变过程，得到粒子运动的瞬时坐标，求解微分运动方程，确定均方位移,可由式（4）得到∑=-=ninnortrntMSD12)()(1)(（4）式中，r n (t)为t时刻原子n的位置。

扩散系数可由式（5）得到：)(61MSD tD ∂∂=（5） 即扩散系数为均方位移斜率的1/6。

图3为不同钙硅比的水化硅酸钙中各种原子在300K 时的均方位移曲线。

从图中可以看出，每种原子的MSD 曲线均随着时间的延长而呈上升趋势，这说明在体系内部发生了原子扩散。

Ca/Si=0.83时各原子的MSD 曲线上升的幅度最大，Ca/Si=1.5时最小，这说明Ca/Si=0.83时各原子更容易产生位移而扩散。

图3（d ）是H 原子的MSD 曲线，从中可以看出，Ca/Si=0.83时其扩散系数远大于其他钙硅比。

原因可能是Ca/Si 由0.67增加到0.83时，硅氧链上的H 原子减少，使得体系内原子间相互作用力减小，各原子更易扩散。

而Ca/Si 继续增加，层间的Ca 增加，体系结构更紧密，原子间相互作用力增加，使得扩散变得困难。

M S DT/psM S DT/ps(a)Ca (b)SiM S DT/psM S DT/ps(c)O (d)H图3 均方位移曲线 Fig. 2 MSD curves4 结论以Hamid的托贝莫来石模型和Bonaccorsi的羟基硅钙石模型为基础，通过分子动力学方法计算得到了无定形水化硅酸钙的模拟结构，并分析了其径向分布函数、配位数和均方位移，得到了以下结论：（1）不同钙硅比的初始结构输入，模拟后均能得到无定形态的水化硅酸钙，且其各原子对的键长、配位数相等。