水泥水化二维模拟的数值方法及其应用

袁 庆,郑建军,邢鸿燕

(1.浙江工业大学建筑工程学院,浙江 杭州 310014;2.中天建设集团有限公司,浙江 杭州 310014)

20世纪80年代以来,随着计算机技术的飞速发展,除了以往的试验手段外,计算机模拟已成为解决材料科学理论和工程问题的重要手段。通过模拟水泥水化来研究混凝土细观结构与宏观物理力学性能之间的关系已逐步为国内外学者所关注[1-4]。目前主要有2类有代表性的水泥基材料水化模型:数字图像基模型和连续基模型[2-3]。Wittmann等[5]最早提出应用计算机来模拟水泥或混凝土微观结构,建立数字混凝土模型。随后,Bentz等[6]提出基于数字图像处理的水泥水化模型,将每一个水泥颗粒描述成一些图像像素的集合体,通过一套元胞自动机规则操纵全部像素来模拟水泥的水化过程[3,7],可描述多尺度、多相、非球形水泥颗粒。其局限性在于它的解析度,因为每个像素所占的体积为1 μ m3,无法描述小于该尺度的特征量。Jennings等[1]于1986年首先将水泥颗粒作为基本单元来描述水泥浆体的微观结构[1],创立了连续描述法。van Breugel等[2,8]开发了HYMOSTRUC模拟系统,考虑了水泥矿物组成、水泥颗粒分布、矿物掺和料、水灰比、养护温度等技术参数对水化过程的影响。Stroeven等[9]开发了SPACE系统,以水泥颗粒的动态混合过程作为出发点,较理想地模拟了实际生产过程中水泥颗粒的堆积过程。这些水泥水化模拟都是在三维空间中实现的,三维模拟的优点是能够比较真实地重现水泥水化物理化学过程,缺点是对计算机内存和速度要求非常高,计算时间长,特别是随后的宏观物理力学性能分析极其困难,甚至目前还无法实现,因此,研究二维水泥水化尤为必要,尽管其模拟精度不如三维水泥水化的高,但计算时间短,对计算机内存要求低,为下一步宏观物理力学性能的预测提供了可能性。为此,本文在前人工作的基础上,在水泥水化二维数值模拟方面进行了初步探索。

1 水泥颗粒初始分布

在二维平面上分析水泥水化特性时,通常将水泥颗粒模拟成圆。根据体视学原理,关于水泥颗粒个数的累积分布函数PN(d)为[10]

式中:d为圆形水泥颗粒直径;Dm为最大水泥颗粒直径;D0为最小水泥颗粒直径;α1和 α2为参数,对大多数水泥而言,α1=0.038,α2=0.980[11]。

在模拟水泥颗粒初始分布时,取一边长为a的正方体单元,这样,当D0,Dm和水灰比w/c给定时,就可以采用文献[10]的方法生成各种尺寸的水泥颗粒。水泥颗粒分布的基本原则是任何2个水泥颗粒的中心距离大于或等于它们的半径之和。为了消除水泥颗粒分布的边界效应,在正方形单元中引入周期性边界条件,即如果某一水泥颗粒与正方形边界相交,将该水泥颗粒位于单元外面的那部分反射到对边的边界上。文中算例D0=1μ m,Dm=15μ m,a=200μ m,w/c=0.5,水泥颗粒初始分布如图 1所示,图中黑色圆圈表示水泥颗粒,白色部分表示水。

图1 水泥颗粒的初始分布

2 水泥水化模拟

水泥水化是一个极其复杂的物理化学过程,受水泥化学组成、水泥颗粒尺寸分布、水灰比、温度等诸多因素的影响。在水泥水化过程中,水泥表面各点的水化程度不一致,为方便起见,通常假设水泥表面各点的反应速率相同,这样在任一时间可用3个同心圆来表示水泥颗粒,如图 2所示(图中Rin为未水化水泥颗粒半径,Rout为水化产物的外半径,Rair为空气层的外半径)。需要指出的是定义空气层是为了确定水化过程中每个水泥颗粒的外水化产物表面与自由水直接相邻的部分,用以确定下一步水化程度及计算局部水的消耗。根据水化动力学原理,生成的水化产物面积ΔSg和所消耗的水面积 ΔSw与已经水化的水泥面积 ΔSc之间存在着如下关系[11]:

图2 以同心圆表示的水泥颗粒

式中:k0和k1为常数,一般k0介于1.7～2.2之间,k1介于1.25～1.4之间。

整个水泥水化过程由结晶成核与晶体生长、相边界反应和扩散 3种机理所控制,由于第1种机理复杂,而且仅在水化度在1%～2%之间起作用,模拟时可以忽略不计。水化产物的厚度为

试验研究表明,当 δ小于临界厚度δtr时,反应速率由相边界反应机理所控制;当 δ大于临界厚度δtr时,反应速率由扩散机理所控制。这样,未水化水泥颗粒半径的减小量可表示成[11]

式中:Δt为时间间隔;β为控制扩散过程系数;H为阶跃函数;KT为相边界反应机理控制阶段的反应速率,它与水泥石养护温度T之间存在如下关系[2]:

式中:K0为温度为20℃时的反应速率;E为激活能;R为气体常数。对于一般水泥,E/R可取5364K-1。

由于生成的水化产物体积大于已经水化的水泥体积,而所消耗水的体积也不断增加,因而随着水泥水化的进行,各水泥颗粒的水化产物层和空气层不断增大,导致相邻水泥颗粒之间相互干扰。将每个同心圆分为2个部分:不受相邻水泥颗粒影响的活动区域和受相邻水泥颗粒影响的非活动区域。为了量化水泥颗粒之间的这种相互干扰效应,引入修正参数 ω0,ω1和 ω2分别表示未水化水泥、水化产物和空气层表面活动区域长度与相应周长之比。很显然0≤ωi≤1(0≤i≤2),而且 ωi越小,水泥颗粒之间的干扰效应越大。对于给定的时间t,通过数值方法容易确定 ωi值。有了 ωi值,未水化水泥颗粒半径、水化产物外半径和空气层外半径的增量可通过式(6)～(8)计算:

这样,t+Δt时刻的未水化水泥颗粒半径、水化产物外半径和空气层外半径分别为

考虑图1所示的水泥颗粒初始分布,通过水泥水化模拟就可得到水化1d,14d和28d后的水泥石微观结构,如图3所示。在图3中,黑色部分表示未水化水泥颗粒,灰色部分表示水化产物,白色部分表示包含水和空气的毛细孔。从图3可以看出,随着水化时间的增加,未水化水泥和毛细孔部分不断减小,而水化产物部分不断增大。

3 试验验证分析

为了验证本文数值方法的有效性,选用文献[8]的试验结果与本文方法模拟结果进行比较。在该试验中,所用水泥的各矿物组分的体积分数如下:C3S为56.7%,C2S为17.2%,C3A为6.7%,C4AF为7.9%,水泥勃氏比表面积为312m2/kg,取水泥颗粒最大直径为15μ m,水泥颗粒最小直径根据所给定的水泥勃氏比表面积进行反算[12],结果为 2.12 μ m。整个试验在等温条件下进行,温度为 20℃,水灰比w/c分别为0.3,0.4和0.5,所测得的水化度 α与时间t的关系如图4所示。在模拟过程中,将每一时刻所确定的每个未水化水泥颗粒半径Rin代入式(13)计算水化度:

式中:Rin,0为初始时刻未水化水泥颗粒的半径;∑表示对所有水泥颗粒求和。

图4 α～t关系模拟结果与试验结果的比较

由图4可见,数值模拟结果与试验结果吻合良好,当w/c为0.3,0.4和0.5时,数值模拟结果与试验结果之间的平均相对误差分别为9.37%,6.38%和5.84%。因此,本文数值方法的有效性得到了试验结果的验证。下面讨论养护温度对水化度的影响 ,设D0=1 μ m,Dm=15μ m,w/c=0.5,养护温度分别为10℃,20℃和30℃,结果如图5所示。由图5可见,对于给定的水化时间,水化度随着养护温度的升高而增大,主要原因是反应速率随养护温度的升高而增大。当养护温度从10℃升高到30℃时,水化时间为1d,7d,14d和28d时的水化度分别增大126.58%,32.95%,19.15%和13.64%,因此,养护温度对水泥浆水化度的影响早期比晚期更为明显。

4 水泥石微观结构初步分析

应用计算机模拟生成水泥石微观结构后,可以通过进一步分析获得宏观物理力学性能预测所需要的一些微观结构参数。在下面的所有计算中,a=200μ m,w/c=0.5,T=20℃,D0=1μ m,Dm=15μ m 。

图5 养护温度对水化度的影响

首先考虑毛细孔、未水化水泥和水化产物的面积百分数Acap,Ac和Ag,这三相组分面积百分数是预测和评价水泥基材料宏观物理力学性能最基本的参数[13]。采用蒙特卡罗方法[14]进行计算,即先在模拟区域中产生N个随机点,再计算落入每一相的随机点数,则每一相的面积百分数就等于落入该相的随机点数与N之比,其计算结果如图6所示。从图6可以看出,毛细孔和未水化水泥的面积百分数随水化时间的增大而减小,而水化产物面积百分数则随着水化时间的增大而增大。图6还显示,当水化时间小于7d时,这三相组分面积百分数变化率较大,而当水化时间超过7d时,这三相组分面积百分数变化率相对比较稳定。就毛细孔而言,1d,3d,7d和15d的面积百分数比28d的面积百分数分别大170%,102%,54%和22%。

图6 各组分面积百分数与养护时间的关系

除了毛细孔面积百分数外,毛细孔两点概率函数S也是评价水泥石传递系数和弹性模量的重要参数,它可以衡量两点距离为r处的毛细孔之间的相关性[13],其计算结果如图7所示。从图7可以看出,所有曲线具有类似性,即毛细孔两点概率函数随着距离的增大而减小,当r大于 12μ m时,S趋于常数,表明毛细孔之间相关性非常小,而且该常数等于r=0时的S2,这与理论分析结果完全一致。另外,图7还显示,对于给定的r值,两点概率函数随着水化时间的增大而减小,表明毛细孔之间的相关性减小。

图7 毛细孔两点概率函数与两点距离的关系

最后分析单位面积毛细孔周长Ca,这也是一个重要的微观结构参数。计算时先给每个水泥颗粒的水化产物外半径一个增量P,再计算毛细孔面积减少量Scap,则Ca就等于Scap与P的商。在给定毛细孔面积百分数的情况下,Ca越大,水泥石中的毛细孔越曲折,离子和液体的渗透系数越小[13]。计算得出3d,7d,15d和28d的Ca值比1d的Ca值分别小6%,16%,30%和38%,Ca随着水化时间的增大不断减小,这表明毛细孔的曲折程度随水化时间的增大而降低,原因可能是毛细孔面积百分数随着水化时间的增大而不断减小。

5 结 论

a.利用水化动力学原理,研究了水泥水化二维模拟的数值方法,该数值方法的有效性得到了试验结果的验证。

b.分析了温度对水化度的影响,当养护温度从10℃升高到30℃时,水化时间为1d,7d,14d和 28d时的水化度分别增大126.58%,32.95%,19.15%和13.64%。

c.基于模拟所得的水泥石微观结构定量分析了毛细孔、未水化水泥和水化产物的面积百分数,毛细孔两点概率函数以及单位面积毛细孔周长随水化时间的变化规律。

:

[1]JENNINGS H M,JOHNSON S K.Simulation of microstructure developoment during the hydration of a cement compound[J].Journal of the American Ceramic Society,1986,69(11):790-795.

[2]van BREUGEL K.Numerical simulation of hydration and microstructural development in hardening cement-based materials:(Ⅰ)theory[J].Cement and Concrete Research,1995,25(2):319-331.

[3]BENTZ D P.Three-dimensional computer simulation of cement hydration and microstructure development[J].Journal of the American Ceramic Society,1997,80(1):3-21.

[4]BULLARD J W.A three-dimensional microstructural model of reactions and transport in aqueous mineral systems[J].Modelling and Simulation in Materials Science and Engineering,2007,15(7):711-738.

[5]WITTMANN F H,ROELFSTR A P E,SADOUKI H.Simulation and analysis of composite structures[J].Material Science Engineering,1984,68(2):239-248.

[6]BENTZ D P,GARBOCZI E J.Percolation of phases in a threedimensional cement paste microstructural model[J].Cement and Concrete Research,1991,21(2):325-344.

[7]BENTZ D P,COVENY P,GARBOCZI E J.Cellular automation simulations of cement hydration and microstructure development[J].Modeling andSimulation inMaterialsScience and Engineering,1994,2(4):783-808.

[8]van BREUGEL K.Numerical simulation of hydration and microstructural development in hardening cementbased materials:(Ⅱ)applications[J].Cement and Concrete Research,1995,25(3):522-530.

[9]STROEVEN M,STROEVEN P.SPACE system for simulationof aggregated matter application to cement hydration[J].Cement and Concrete Research,1999,29(8):1299-1304.

[10]ZHENG Jian-jun,LI Chun-qing,ZHOU Xin-zhu.Thickness of interfacial transition zone and cement content profiles around aggregates[J].Magazine of Concrete Research,2005,57(7):397-406.

[11]STROEVEN M.Discrete numerical modelling of composite materials[M].Delft:Melnema,1999.

[12]郑建军,庞宪委,熊芳芳.水泥浆体初凝时间预测的数值方法及其影响因素评价[J].硅酸盐学报,2008,36(7):927-933.

[13]TOR QUA TOS.Random heterogenerous materials:microstructure and macroscopic properties[M].New York:Springer-Verlag,2001.

[14]ZHENG Jian-jun,XIONG Fang-fang,WU Zhi-min,et al.A numerical algorithm for the ITZ area fraction in concrete with elliptical aggregate particles[J].Magazine of Concrete Research,2009,61(2):109-117.