凝胶孔对水化硅酸钙(C-S-H)力学性能影响的分子动力学模拟

2023-12-01李宗利刘士达童涛涛肖帅鹏李云波

材料科学与工程学报 2023年5期

李宗利,刘士达,童涛涛,肖帅鹏,李云波

(1.西北农林科技大学水利与建筑工程学院, 陕西杨凌 712100; 2.西北农林科技大学旱区农业水土工程教育部重点实验室, 陕西杨凌 712100; 3.山东电力工程咨询院有限公司, 山东济南 250013;4.中水北方勘测设计研究有限责任公司, 天津 300222)

1 前言

混凝土具有多尺度结构特征,不同尺度性能对混凝土宏观性能影响不同。混凝土材料的力学性能主要决定于水泥水化产物,而水泥水化产物中最重要的成分是水化硅酸钙(C-S-H)凝胶体,其占硬化水泥浆体体积的60%～70%,对水泥基材料的行为起着决定性的作用[1]。C-S-H 是水泥基材料组成的基本单元,其特性是解释水泥基材料宏观性能的基础。

在水泥水化过程中,起初由水、空气和水泥熟料占据的部分空间会逐渐被水化产物所填满,未被水化产物填充的空间将逐渐被分割成大小不同的孔隙。水泥基材料的孔隙也具有多尺度特征:孔隙从毫米尺度到纳米尺度的大小跨度超过6个数量级。根据孔隙的大小,可以将水泥浆体中的孔隙分为凝胶孔、过渡孔、毛细孔和大孔四种类型。从分子水平上看,C-S-H 为无序的链状硅酸盐,长度为4～10 nm[2],其内部含有大量的凝胶孔,孔径为0.5～5 nm[2]。C-S-H 凝胶孔可以看作是C-S-H 纳米颗粒之间的孔隙,如图1所示。Jennings[3]在C-S-H 第一代模型CM-Ⅰ的基础上,建立了C-S-H 纳米结构的第二代模型CM-Ⅱ,如图1(a)所示。该模型可以很好的描述C-S-H 的孔隙和层间结构。Qomi等[4]在C-S-H 干燥收缩的研究中,对含有凝胶孔的C-S-H 模型的描述如图1(b)所示。模型中心为凝胶孔,周围是层状的C-S-H。图1(c)为含凝胶孔的C-S-H 分子模型[5]。

图1 C-S-H 纳米结构模型 (a)第二代C-S-H 模型CM-Ⅱ[3];(b)含凝胶孔的C-S-H 模型[4];(c)含凝胶孔的C-S-H 分子模型[5]Fig.1 Nanoscale model of hydrated calcium silicate (a) the C-S-H model CM-Ⅱ[3]; (b) C-S-H model with gel pores[4];(c) molecular model of C-S-H with gel pores[5]

近年来,学者们利用分子动力学对C-S-H 的力学性能进行了研究。Hou 等[6]通过单轴拉伸研究了Ca/Si比对C-S-H 力学性能的影响,研究表明随着Ca/Si比值的增加,C-S-H 的力学性能不断下降。为了了解C-S-H 凝胶体的抵抗载荷能力,Jin等[7]采用分子动力学模拟C-S-H 沿不同方向的单轴拉伸试验,结果表明,由于C-S-H 在不同方向上结构不同,拉伸过程中表现出各向异性。马彬等[8]和刘士达等[9]分别研究了水分子对C-S-H 结构和力学性能的影响,结果表明,随着水分子含量的增多,C-S-H 层间距离增大,体积模量、剪切模量和杨氏模量均减小。以上文献虽然对C-S-H 的力学性能展开了研究,但并没有考虑凝胶孔的影响。Bauchy等[10]通过手动删除原子的方法,建立了含不同长度椭圆微裂纹的C-S-H 模型,并采用分子动力学模拟了带裂纹的C-S-H 的单轴拉伸,计算了单轴拉伸过程中C-S-H 晶粒的表面能和断裂韧性。Hou等[5]建立了含直径0.5～5 nm 凝胶孔隙的C-S-H 模型,并模拟其单轴拉伸。结果表明,由于凝胶孔的存在,C-S-H 的杨氏模量和抗拉强度显著降低。以上研究虽然考虑了凝胶孔对C-S-H 力学性能的影响,但是都仅仅针对单轴拉伸状态下的力学性能。不同受力状态下C-S-H 力学特性不同,凝胶孔对不同受力状态下力学性能的影响程度也不同,需要进一步系统地展开研究。

本研究利用分子动力学,首先分析了凝胶孔对C-S-H 弹性力学参数的影响,然后模拟了含不同尺寸凝胶孔的C-S-H 在单轴拉伸、单轴压缩、剪切不同受力状态下的应力应变关系,对比分析了凝胶孔对不同受力状态下的影响程度,系统揭示纳观尺度下凝胶孔对C-S-H 力学性能影响规律。

2 模拟方法

2.1 力场的选择和模型的建立

力场是表示分子系统势能变化的函数,选择合理与否影响着仿真的精度[6]。水化硅酸钙力场(CSHFF)是基于粘土力场(ClayFF)[11]开发的,主要通过改进ClayFF 中Ca、Si和O 的参数,使C-S-H 的结构特征和高阶特性体现的更精确[12]。该力场已被证明可以正确预测水泥水合物的基本结构和物理特征,并且在C-S-H 力学性能的分子模拟中已被广泛认可和应用[13-15]。因此,本研究选用CSHFF 力场作为分子模拟力场。Ca、Si、O、H 原子间相互作用参数通过文献 [11,15]可获得。

基于Pellenq 等[17]所提供的方法建立了C-S-H模型。该方法基于Tobermorite模型删除部分SiO2所建立的“真实”模型能够很好地匹配实际C-S-H 的各种参数,因此在C-S-H 分子模拟中被广泛应用[13,18-19]。已有研究证明所建模型能够很好地预测C-S-H 的物理力学性能。建立过程如下:

(1)以Tobermorite 11 Å 晶体(图2(a))作为C-S-H 模型的初始结构,根据核磁共振(NMR)测试提供的Qn分布进行断链[20],即Q0=0%,Q1=66.01%,Q2=33.99%,平均链长(MCL=2(Q2/Q1+1)=3.03)与NMR 测试结果保持一致,并参照Allen等[21]通过实验所测得的C-S-H 的平均分子式,Ca/Si比为1.7。通过上述步骤得到的弛豫前C-S-H 模型(图2(b))在300 K、等温等压(NPT)系综下弛豫,得到弛豫后的C-S-H 模型。

图2 含凝胶孔的C-S-H 模型的建立过程 (a) 11 Å Tobermorite; (b) 弛豫前的C-S-H; (c)弛豫后的C-S-H; (d) 饱和C-S-H; (e) C-S-H 超胞; (f) 含凝胶孔缺陷的C-S-H 超胞Fig.2 Establishment process of C-S-H model with gel pores (a) Tobermorite 11 Å; (b) unrelaxed C-S-H; (c) relaxed C-S-H;(d) saturated C-S-H; (e) C-S-H supercell; (f) C-S-H supercell with gel pore defect

(2)采用蒙特卡洛方法对弛豫后的C-S-H 进行吸水,吸附过程在300 K、NPT 系综下完成,吸附步数设置为1×109步。当吸附过程在饱和蒸汽压下达到平衡状态时,得到如图2(d)所示的饱和C-S-H。最终建立的模型在分子式、密度、Qn物种分布和平均链长等方面与实验结果接近[20-22]。

(3)为了得到稳定的统计结果,获得可靠的实效模式[23],将上述建立的饱和C-S-H 模型扩展成1×3×4的超胞,如图2(e)所示。

构建凝胶孔的方法有两种。第一种方法是通过随机堆叠小的C-S-H 颗粒,在C-S-H 颗粒之间自然形成凝胶孔;第二种方法是在较大的C-S-H 超胞中删除一部分原子来建立凝胶孔[3,5,24]。本研究采用第二种方法构建凝胶孔。通过手动删除C-S-H 超胞y-z(坐标轴如图2所示)平面上中心原子的方法,建立直径分别为0、10、15、20、25和30 Å 的凝胶孔,进一步弛豫后得到含不同尺寸凝胶孔缺陷的C-S-H 超胞膜型(图2(f))。该模型的孔径分布在0.5～5 nm 的凝胶孔径范围内。在删除原子过程中,为了避免水分子产生影响,只删除C-S-H 中的钙、硅原子以及游离的氧原子,将水分子保留在凝胶孔中。

图2展示了含凝胶孔的C-S-H 模型的建立过程。在水泥基材料中,C-S-H 的Ca/Si比变化范围较大。Ca/Si比越大,C-S-H 中硅酸盐链会由长链变为短链,内部凝聚力降低,力学性能降低。这与相关研究结果一致[6,25]。

2.2 分子动力学模拟

受目前常用的分子模拟软件分析功能限制,C-S-H 弹性力学参数和不同受力状态的模拟分别采用不同软件进行。

C-S-H 的弹性力学参数和剪切模拟均采用Materials Studio软件进行模拟。模拟过程中通过对C-S-H 结构施加应变荷载,分析结构变形得到C-S-H的弹性常数。C-S-H 的体积模量、剪切模量根据弹性常数和柔性常数,采用Voigt-Reuss-Hill平均值法计算,杨氏模量根据柔性常数进行计算。模拟剪切时,模型的三个方向均采用周期性边界条件,模拟过程中对超胞在y-z平面上施加剪应变,应变速率为0.005/ps。在每次施加完成后,对超胞进行弛豫,弛豫时间为10 ps,并将充分弛豫后的模型作为下一次施加的初始模型。剪切和弛豫过程均采用NVT 系综,模拟时间步长设置为1.0 fs。

单轴拉伸和压缩采用LAMMPS(大规模原子/分子并行模拟器)软件模拟,该软件是一个经过良好测试和广泛使用的开源经典分子动力学软件[26],在模拟单轴拉伸和单轴压缩方面较为方便。对第2.1节所建立的含不同尺寸凝胶孔的C-S-H 超胞模型进行充分弛豫,作为单轴拉伸、压缩模拟的模型。拉压沿C-S-H的层状和链状方向进行,即y和z方向(见图2)。模型三个方向均采用周期性边界条件,模拟过程中对模型施加应变,应变速率保持0.08/ps不变。当沿一个方向进行拉伸或压缩时,其他两个方向的压力保持为零,变形自由,体现泊松效应。垂直于加载方向的无压力设置消除了对变形的人为约束,允许张力自由发展,不受任何约束[27]。整个系统的弛豫和模拟均在NPT系综下进行,模拟时间步长1.0 fs。

图3为C-S-H 在不同受力状态下的荷载示意图。通过拉伸、压缩、剪切得到不同受力状态下的应力-应变曲线,并通过梯形积分规则计算应力-应变曲线与坐标轴围成的面积,得到拉伸、压缩和剪切时的单位体积应变能密度,用应变能密度作为评价系统变形性能的间接指标。

图3 C-S-H 不同受力状态下的荷载示意图 (a)、(b)为y 和z 方向的单轴拉伸;(c)、(d)为y 和z 方向的单轴压缩;(e)为剪切Fig.3 Load diagrams of C-S-H under different stress states.uniaxial tension along (a) y and (b) z directions,uniaxial compression along (c) y and (d) z directions, (e) shear states

3 结果与讨论

3.1 模型和力场的检验

原子径向分布函数(RDF)可以解释为系统的区域密度与平均密度的比值,RDF的峰值所对应的位置表示两种原子之间的键长。对第2.1节所建立的不含凝胶孔的C-S-H 超胞模型进行充分弛豫之后,计算了层间钙与桥接氧Caw—O、层内钙与桥接氧Cas—O、硅原子与桥接氧Si—O 的RDF,如图4所示。结果显示,C-S-H 模型中硅原子与氧原子的键长为(1.61±0.1) Å,层内钙与氧原子的键长为(2.31±0.2) Å,层间钙与氧原子的键长为(2.37±0.22 ) Å。Richardson[28]的实验数据显示,Si—O 平均键长为1.639 Å,Ca—O 键长在2.2～2.9 Å 范围之内;X 光散射实验[29]表明,Si—O 平均键长为1.67 Å,Ca—O 平均键长为2.43 Å;Svenum 等[30]使用密度泛函理论计算得到C-S-H 中Si—O 平均键长为1.65 Å,Ca—O 键长为2.3～2.5 Å。由此可见,所建模型的Ca—O 和Si—O键长与实验数据和密度泛函理论数据一致,说明模型具有可靠性。

图4 C-S-H 中Si—O、Cas—O、Caw—O 径向分布函数Fig.4 Radial distribution function for Si—O, Cas—O and Caw—O in C-S-H

使用CSHFF力场初步计算了C-S-H 的弹性力学参数,包括体模量、剪切模量和杨氏模量。将计算值与之前实验和模拟结果进行比较,如表1所示。结果表明,所建模型的弹性力学参数均在合理范围内,说明该模型和力场具有一定的可靠性。

表1 C-S-H 的弹性力学参数与之前实验和模拟值的比较Table 1 Elastic mechanical parameters of C-S-H in comparison with the previous experiments and simulations

3.2 凝胶孔对C-S-H 弹性力学参数的影响

从图5 中可以看出,随着凝胶孔尺寸的增大,C-S-H 的体模量、剪切模量和杨氏模量均不同程度地下降,但对x轴的杨氏模量影响较其他轴要小。说明凝胶孔会降低C-S-H 的弹性力学参数,且凝胶孔尺寸越大,降低程度越大。

图5 含不同尺寸凝胶孔的C-S-H 的体积模量、剪切模量和杨氏模量Fig.5 Bulk modulus, shear modulus, and Young’s modulus of C-S-H with different size gel pores

3.3 凝胶孔对C-S-H 应力应变关系的影响

图6为含有不同尺寸凝胶孔的C-S-H 沿y、z方向拉伸时的应力-应变曲线。初始阶段应力随着应变近似线性增加,当达到峰值点后,应力开始下降,最终减小为零,说明C-S-H 已经被完全拉断。随着凝胶孔尺寸的增大,峰值应力不断减小,应力减小为零所需的应变也不断减小,断裂时最大应变值也不断减小,这说明凝胶孔的尺寸越大,对C-S-H 的力学性能和变形性能的影响程度越大,C-S-H 更易被拉断。图中还展示了凝胶孔直径为30 Å 的C-S-H 在拉伸时的结构变化。随着应变的增加,孔隙开始扩大变形,裂纹首先在凝胶孔的边缘出现,然后不断发展,直至整个C-S-H完全断裂。含其他尺寸凝胶孔的C-S-H 在拉伸时的结构变化与之类似。

图6 C-S-H 单轴拉伸应力-应变关系 (a) y 轴方向; (b) z 轴方向Fig.6 Stress-strain relations of C-S-H under uniaxial tensions (a) y direction; (b) z direction

如图7所示,首先应力在初始阶段几乎随应变线性增加,当达到峰值后,应力开始缓慢下降并最终趋于平稳,说明C-S-H 在压缩作用下并未完全破坏。同样,随着凝胶孔尺寸的增大,峰值应力不断减小。这表明,凝胶孔尺寸越大,C-S-H 的抗压强度越低。凝胶孔对压缩时不同方向影响程度也存在差异。比较图6和图7可以看出,凝胶孔对抗压强度的影响远远小于对抗拉强度的影响。从结构变化来看,在压缩过程中,C-S-H 在压缩方向逐渐致密,凝胶孔被逐渐压缩变形,甚至消失。与拉伸不同的是,压缩时C-S-H 最终不会完全破坏。

图7 C-S-H 单轴压缩应力-应变关系 (a) y 轴方向; (b) z 轴方向Fig.7 Stress-strain relations of C-S-H under uniaxial compression (a) y direction; (b) z direction

图8显示与拉伸和压缩时不同,当剪应力达到峰值后,并没有呈现下降趋势,而是在模拟范围内呈锯齿状波动并缓慢增加。峰值应力也随着凝胶孔尺寸的增大而减小,说明凝胶孔尺寸越大,对抗剪强度的影响程度越显著。从结构变化可以看出,C-S-H 在剪切过程中变形不断增大,凝胶孔也随着剪应变的增加不断变形。

图8 C-S-H 剪切应力-应变关系Fig.8 Stress-strain relations of C-S-H under shearing

以上通过对应力-应变曲线的简单分析可以看出,凝胶孔会降低C-S-H 的力学性能。同时,C-S-H 在不同的受力状态下表现出不同力学性能和变形有差异。为了进一步分析凝胶孔对不同受力状态下C-S-H 力学性能和变形性能的影响程度,表2列出了含不同尺寸凝胶孔的C-S-H 在三种受力状态下的应力应变特征参数。从峰值应力的减小程度来看,当凝胶孔直径从10 Å 增加到30 Å 时,在受拉状态下,y方向的降低程度从4.97%增加到23.09%,z方向的降低程度从4.07%增加到27.10%;在受压状态下,y方向的降低程度从1.56%增加到15.81%,z方向的降低程度从1.48%增加到16.20%;受剪时,降低程度从2.15%增加到18.10%。说明凝胶孔对抗拉强度的影响程度最大,对抗压和抗剪强度的影响程度较小且几乎相同;从方向上来看,不管是抗拉还是抗压强度,凝胶孔对z方向的影响程度都略大于y方向。应力-应变曲线和坐标轴所围成的面积表示应变能密度,反映了C-S-H破坏前吸收的能量大小和变形性能。从应变能密度的降低程度来看,当凝胶孔直径从10 Å 增加到30 Å 时,在受拉状态下,y方向的降低程度从11.12%增加到48.62%,z方向的降低程度从27.42% 增加到66.25%;在受压状态下,y方向的降低程度从0.04%增加到7.78%,z方向的降低程度从0.62%增加到7.90%;受剪时,降低程度从4.39%增加到22.30%。这说明凝胶孔对C-S-H 拉伸变形性能影响程度最大,其次是剪切,对压缩变形性能影响程度最小;从方向上来看,不管是拉伸还是压缩,凝胶孔对z方向的变形性能影响程度都大于y方向。