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 。

8

6

)exp()(ij ij ij ij ij

ij

ij 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 –H

0.372

0.029

三体作用势是一种比van de Waals 力作用势更短程、更弱的作用势，其对总能量的贡

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

20)cos )(cos exp(),,(θθγ

γλθ--+-=ijk cik

ik ik

cij ij ij

ijk 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-H

32.710

0.13

0.16

104.5°

模拟采用三维周期边界条件(periodic boundary condition,PBC)，对长程相互作用采用Ewald 求和方法，运动方程积分采用“蛙跳”(leap frog)算法，积分步长为1 fs ，速度由Maxwell 分布给出。为了得到无定形结构中的长程无序、短程有序及硅氧四面体链旋转、扭曲结构，先将体系在5000K 下的NVT 系综中预平衡，然后将体系的温度逐步降低到300K 。体系在300K 重新达到平衡后，将其输出作为NVE 系综的输入，最终制得稳定的无定形水化硅酸钙，然后对模拟结果进行分析计算，得出水化硅酸钙在300 K 时的结构性质及动力学行为。

3 模拟结果与讨论

3.1径向分布函数(RDF)

分子动力学模拟无定形体系的结构信息可以用RDF 来分析，表示系统的区域密度与平均密度之比，计算公式如下

ρ

π2

4)(r dN

r 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 Å，此即为对应的原子对间的原子距离，与实验值基本符合。比较不同钙硅比的水化硅酸钙的总径向分布函数，