分子动力学模拟研究混合碱土金属对无碱硼铝硅酸盐玻璃结构和性能的影响

王照猛，朱雪梅，陈淑勇，李从云，马立云

(1.中建材蚌埠玻璃工业设计研究院有限公司，浮法玻璃新技术国家重点实验室,蚌埠 233010 2.硅基材料安徽省实验室，蚌埠 233010)

0 引 言

无碱硼铝硅酸盐玻璃被广泛应用于电子玻璃、光致变色玻璃、玻璃纤维和放射性废物的固化等领域[1]。TFT-LCD基板玻璃属于无碱硼铝硅酸盐玻璃[2]。同时TFT-LCD液晶显示对超薄基板玻璃提出了高化学稳定性、高热稳定性、无碱无砷、低密度以及高弹性模量等性能要求[3]。因此，近年来人们对无碱硼铝硅酸盐玻璃的结构和性能进行了深入研究，并努力探索性能与结构之间的关系[4-5]。

然而，宏观性质是由微观结构决定的，为了更好地理解宏观性质，研究相应玻璃体系的组成、结构和性能之间的关系至关重要。分子动力学(Molecular Dynamics，MD)模拟是研究玻璃材料微观结构信息的一种有效方法，可以直接观察到原子结构，但很难通过实验方法获得[6]。Du等[7]采用分子动力学模拟研究了硼硅酸钠和铝硼硅酸盐玻璃的组成、表面结构和性能之间的关系，得到理论模拟数据与实验数据相吻合，同时其在研究钠-钙-锶-铝硼硅酸盐玻璃中也得到了类似的结果[8]。

碱土金属氧化物是玻璃形成外体中游离氧的主要提供者，起断键和电荷补偿作用，同时并对玻璃网络形成体的配位数有很大的影响。因此，碱土金属对无碱硼铝硅酸盐玻璃的物理化学性能的影响具有多面性。此外，不同碱土离子因其离子场强度不同而改变了铝硼硅酸盐玻璃的结构。Wu等[9]研究发现，具有高场强的碱土金属离子更容易促进非桥氧形成，同时，高电场强度的阳离子也会增加铝硅酸盐玻璃中五配位铝的形成。

因此，本研究的目的是通过MD模拟研究混合Ca/Mg对无碱硼铝硅酸盐玻璃的结构和性能的影响。具体计算了阳离子的局部环境和中程结构信息，并讨论了混合Ca/Mg对弹性模量和剪切模量的影响。

1 模拟过程

为了研究混合碱土金属对无碱硼铝硅酸盐玻璃结构和力学性能的影响，B2O3-Al2O3-SiO2的摩尔分数恒定，改变Ca/Mg比设计6组无碱硼铝硅酸盐玻璃样品。组分和密度如表1所示。

表1 玻璃的化学组成和密度Table 1 Compositions and densities of simulated glasses

采用LAMMPS软件进行分子动力学模拟。在MD模拟中，经验势对于生成真实的结构和描述玻璃性能至关重要。本研究使用的是部分电荷对势，结合了长程库伦相互作用和短程相互作用的Buckingham势。该方法在多组分玻璃的建模中具有良好的可移植性和计算效率。Buckingham势Φij(rij)[10]的计算公式如下：

(1)

式中：rij为原子i和j之间的距离；qi、qj为原子i和j的电荷；Aij、ρij和Cij为势参数[11](如表2)。

表2 Buckingham势参数Table 2 Parameters of the Buckingham potential

所有的模拟都使用了由3 210个原子组成的立方体模拟盒。在立方体模拟盒中，通过实验密度，如表1所示随机生成初始原子位置坐标。初始结构在0 K下进行能量最小化，然后通过正则过程(原子数不变、体积不变、温度不变，NVT系综)在300 K下弛豫60 ps。之后，通过NVT系综，系统在6 000 K熔化100 ps(目的是快速打破模型的初始结构，使原子充分混合均匀，使模型快速熔化。不考虑工业生产实际温度是因为在工业温度下很难打破模型的初始结构，模型熔化慢，增加模拟耗时)，在5 000 K下弛豫100 ps，然后以5 K/ps的冷却速率降温到300 K，冷却时长940 ps。在300 K下，玻璃在等温等压过程(原子数不变、压强不变、温度不变，NPT系综)下平衡20 ps，然后在NVT系综下平衡10 ps。最后，在微正则系综(原子数不变、体积不变、能量不变，NVE系综)下平衡10 ps，消除内部应力，获得玻璃的最终结构。

2 结果与讨论

2.1 径向分布函数

径向分布函数是描述液体和非晶体材料微观结构最常用的数学语言，是结构无序化程度的表征，反映的是从任意指定A原子为中心，半径为r到r+dr的球壳内发现B原子的数目为n(B)，公式如下：

(2)

式中：gαβ为径向分布函数，α、β代表两个原子；r为距中心原子半径;δr为球壳厚度;n(r)为球壳内粒子数;ρβ为理想晶体中B原子的密度；V为体积。

图1为O-Si、O-Al和O-B的偏径向分布函数，从图中可以得出，O-Si的键长约为1.63 Å(1 Å=0.1 nm)，O-Al键长约为1.75 Å，O-B键长约为1.39 Å，与文献值相符[12-14]，说明模拟选择的势函数是可靠的。混合碱土金属Ca/Mg对O-Si的键长变化几乎没有影响，对O-Al的键长有较小影响。与O-Si、O-Al偏径向分布函数不同的是，O-B偏径向分布函数有双峰，分别在1.39 Å和1.47 Å处，说明硼配位数存在三配位和四配位。1.39 Å处峰高于1.47 Å处峰，说明硼的配位数主要以三配位构型存在。从偏径向分布函数的峰宽以及尖锐程度可以得出，Al-O结构单元比Si-O和O-B结构单元松散、稳定性差。

图1 O-Si、O-Al和O-B的偏径向分布函数Fig.1 Partial radial distribution functions of O-Si, O-Al and O-B

2.2 键角分布函数(BAD)

图2为O-Si-O、O-Al-O和O-B-O的键角分布函数，O-Si-O键角分布在108.4°～108.6°范围内，说明在硼铝硅系玻璃中，[SiO4]四面体呈现为规则的四面体结构，且混合碱土金属的变化对[SiO4]四面体构型几乎没有影响。与O-Si-O键角分布相比，O-Al-O键角分布在106.4°～107.3°范围内，说明在硼铝硅系玻璃中[AlO4]四面体呈现不规则的四面体结构，且混合碱土金属的变化导致了[AlO4]四面体结构变形。与O-Si-O和O-Al-O键角分布函数相比，O-B-O键角分布函数分别在109.5°和120.0°附近呈现出单峰，说明硼配位数存在[BO3]平面结构和[BO4]四面体结构，与O-B偏径向分布函数的结果相一致。随着Ca/Mg比的降低，109.5°处峰高度降低同时120.0°处峰高度增加，说明[BO3]含量在增加，[BO4]含量在降低。O-Al-O的键角分布范围比O-Si-O分布范围宽，说明[AlO4]四面体结构的变形大，稳定性差，与偏径向分布函数结果一致。在铝硅酸盐玻璃中[AlO4]四面体结构周围存在碱金属和碱土金属离子来补偿电荷，导致[AlO4]四面体结构易变形，与模拟的结果相一致。

图2 O-Si-O、 O-Al-O和O-B-O的键角分布函数Fig.2 Bond angle distribution functions of O-Si-O, O-Al-O and O-B-O

2.3 配位数

图3(a)为O-Si配位数的百分含量图，从图中可以看出，O-Si的配位数主要为四配位，存在少量的五配位，没有发现三配位的存在。Stebbins[15]在1991年通过NMR发现在石英玻璃中存在Si的五配位结构。Ha和Garofalini[6]通过分子动力学方法研究表明，在硼铝硅酸盐玻璃中存在少量的Si的五配位结构。因此，在无碱硼铝硅酸盐玻璃中存在少量的Si的五配位结构。Horbach和Kob[16]在包含8 016个原子的SiO2玻璃系统中使用BKS势，发现在相对较高的2 750 K温度下只有99%的四配位硅。Hoang[17]在3 000个原子的SiO2玻璃系统中使用Morse势，观察到在低于700 K的温度下，多达2.7%的Si原子处于五配位。Lee等[18-19]发现在模拟SiO2玻璃系统中，缺陷浓度高达24.7%，这表明在模拟结构中不完全四配位的系统是一个不可避免的现象。Vollmayr等[20]通过对粒子A和B的任意二元混合径向分布分析表明，随着冷却速率的降低，系统的局部顺序变得更加明显，表明配位状态与冷却速率有关。

图3(b)为O-Al配位数的百分含量图，可以得出，O-Al的配位数以四配位主，存在少量的三配位和五配位。根据MD模拟结果表明，六倍配浓度很低，在统计误差范围内可以忽略不计，与以往MD模拟结果一致。随着Ca/Mg比的降低，五配位含量先降低后增加，在Ca/Mg=6/9时达到最小值。说明增加Mg2+的浓度对五配位铝的含量的影响不是呈线性关系，而是满足一个最值。Stebbins[21]在2005年通过NMR发现在硼铝硅酸盐玻璃中存在高配位的铝。Ha和Garofalini[6]通过分子动力学的方法研究表明，在硼铝硅酸盐玻璃中存在少量的高配位铝，得出模拟与实验的结构相一致。值得注意的是，由于三配位没有通过实验方法检测到，因此将计算得到的三配位作为模拟方法生成了铝结构缺陷。

图3(c)为O-B配位数的百分含量图，从图中可以得出，O-B的配位数以三配位主，存在一定量的四配位。随着Ca/Mg比的降低，三配位含量增加，四配位含量降低。这由于Mg2+比Ca2+电场强度大，给出游离氧的能力弱于Ca2+，导致玻璃组分中，随着Mg2+的增加，三配位含量增加，四配位含量降低。Stebbins[22]在2010年通过控制温度和退火速率研究硼铝硅酸盐玻璃中硼的配位，研究表明BO3/BO4比与玻璃中的网络修饰离子和样品中的B2O3含量有关。

图3 O-Si、O-Al和O-B的不同配位数含量Fig.3 Content of different CNs of O-Si,O-Al and O-B

2.4 氧种类

在玻璃体系中存在四种类型的氧：桥氧BO(氧原子与两个网络中间体相连)、非桥氧NBO(氧原子与一个网络中间体和外体相连)、自由氧FO(氧原子与两个网络外体相连)和TBO(氧原子与三个网络中间体相连)。如图4为硼铝硅酸盐玻璃体系中不同氧类型的百分含量图。随着Ca/Mg比值下降，BO呈下降的趋势，NBO呈增加的趋势。说明增加Mg2+的浓度对NBO的含量呈线性关系，由于Mg2+比Ca2+电场强度大，给出游离氧的能力弱于Ca2+，导致玻璃组分中，随着Mg2+的增加，NBO呈增加的趋势。与Wu等[9]研究发现具有高场强的碱土金属离子更容易促进非桥氧形成的结果相一致。改变Ca/Mg比会引起BO、NBO和TBO转变。模拟结果证实了FO原子的存在，其含量与BO、NBO和TBO相比非常低，可以忽略不计。在玻璃系统中，TBO已经在早期的核磁共振实验中被发现。例如，在无定形的Al2O3-2SiO2中发现了大量的TBO，占总的氧原子30%[23]。在硅酸盐玻璃中碱金属和碱土金属离子如Na+、Ca2+、Mg2+等会填充在网络周围，形成网络补偿体，破坏玻璃网络结构，增加非桥接氧数目。非桥氧进一步增加玻璃体系的热容和熵来改变离子的扩散和流动性，最终降低玻璃体系的粘度。因此，非氧数目是影响硅酸盐玻璃物理化学性质的重要因素。

图4 玻璃中氧类型及百分含量Fig.4 Oxygen speciation and content of glass

2.5 力学分析

在力学性质分析中，计算了硼铝硅酸盐玻璃的体积模量、剪切模量和杨氏模量。图5为不同组分样品的模量变化。可以得出，在低Mg成分中，随着Mg2+浓度的增加，杨氏模量和剪切模量减小。在高Mg2+成分中，随着Mg2+浓度的增加，杨氏模量和剪切模量增加。杨氏模量都大于65 GPa，满足玻璃基板在运输、操作和加工过程中的下陷要求。体模量无明显规律变化。

图5 玻璃的体模量、剪切模量和杨氏模量Fig.5 Bulk, shear and Young’s modulus of glass

3 结 论

(1)Ca/Mg比对O-Si和O-B的键长几乎没有影响；对O-Al的键长有较小影响。

(2)无碱硼铝硅酸盐玻璃体系中[SiO4]四面体呈现规则的四面体结构，[AlO4]四面体呈现不规则的四面体结构。

(3)O-Si和O-Al的配位数以四配位主，O-Al存在少量的高配位；O-B配位数以三配位主要，存在一定量的四配位。

(4)随着Ca/Mg比的降低，[BO3]含量增加，[BO4]含量降低；桥氧(BO)含量降低，非桥氧(NBO)含量增加。

(5)杨氏模量和剪切模量出现先降低后增加的趋势，在Ca/Mg=6/9时，杨氏模量达到最小值65 GPa。

综合本文模拟分析结果表明，目前的分子动力学方法可以产生合理的硼铝硅酸盐结构，其微观结构和力学性能与文献报道的相一致。因此，本文研究结果对实验设计无碱硼铝硅酸盐玻璃配方有一定的指导作用。