基于Lattice-Boltzmann方法的水泥基材料溶液传输过程模拟

李东遥, 单钰涵, 韩宇栋, 丁小平, 侯东伟,3

(1.上海交通大学 土木工程系, 上海 200240; 2.中冶建筑研究总院有限公司, 北京 100088;3.上海市公共建筑和基础设施数字化运维重点实验室, 上海 200240)

混凝土结构通常处于气、水、土及其耦合的服役环境中.外界的侵蚀性介质,通过扩散、对流、渗透等输运机制到达材料内部,引发溶蚀、化学腐蚀及电化学作用等物理、化学过程,最终导致材料损伤、破坏,这是引发混凝土材料与结构耐久性问题的主要原因[1-2].在海洋环境中,氯离子侵入钢筋表面,与钢筋钝化膜结合并扩散,使钢筋脱钝[3].脱钝后的钢筋在氧气和水的作用下发生电化学腐蚀,进而引起钢筋锈蚀和结构胀裂.钢筋混凝土结构的寿命一般以氯离子到达钢筋表面并达到临界浓度为判据,因此氯离子随海水在混凝土中的传输过程决定了钢筋混凝土结构的耐久性能.海水在混凝土中的传输,通常采用水分传输系数或氯离子扩散系数来统一描述[4-6].

混凝土中,水蒸气和液态水的传输方式不同.水蒸气在浓度梯度的作用下向混凝土内部扩散,而液态水则在压力梯度下发生渗透[7].对于理想状态下单次干燥或湿润过程,可通过水分扩散方程予以描述[8-11];对于氯离子传输过程,传统方法是使用考虑其与孔隙表面的结合效应建立的改进扩散方程[12]描述.但水分迁移系数和氯离子扩散系数都只是一个唯象的综合性参数,并不能反映其微观作用机理;且现有研究还没有反映物理本质的、对溶液传输机制进行描述的方法.实质上,海水是一种“溶液”：从微观机理上,溶液中作为溶剂的水和作为溶质的各类盐离子虽然同时向混凝土内部迁移,但是其迁移规律不同,应被看作相互区别但又相互耦合的2个迁移过程[5].

为了考察溶液在水泥基材料中的传输过程,本文首先采用自生长模型建立了水泥基材料细观结构模型,然后基于格子-玻尔兹曼(lattice-Boltzmann)方法,计算模拟水泥基多孔介质中水分的多相传输和离子的迁移、扩散过程,初步考察以氯盐溶液为代表的“水分-离子”双重迁移机制,并分析材料参数的影响规律.

1 Lattice-Boltzmann方法简介

Lattice-Boltzmann方法(LBM)是近几十年发展的一种在微观尺度上求解流体力学方程的离散化方法.微观流体可看作空间微元中一系列的流体微粒,比分子尺度大,比宏观尺度小.在LBM中,计算区域被离散成标准化网格,时间被离散成时间步,而微粒的运动则采用Boltzmann方程或其离散形式描述.LBM模型通常包括格子(即离散速度集合)、平衡态分布函数及其演化方程3个部分.演化方程通常表示为：

fi(x+ciδt,t+δt)-fi(x,t)=Ωi(x,t)

(1)

式中：x为格点位置;ci为流体粒子离散速度集合,i为离散速度维数;δt为离散时间步长;t为当前时间步;fi为以速度ci运动的速度分布函数;Ωi为碰撞算子,表示微粒间的碰撞对速度分布函数的影响.

Ωi的表达式为：

(2)

LBM由分子动力学方法发展而来,其基本思想是基于粒子碰撞和Boltzmann统计来描述微观物理过程,因此不依赖于宏观连续性假设,可以模拟不规则形状介质中微观-细观尺度上的复杂流动问题.针对具体的物理场景,其无量纲参数依照约定选取,而物理场景中各个物理机制的重要性和主次分别,则根据物理过程和材料特性,通过构建相应的演化模型来实现.近年来,随着LBM理论的发展,学者们建立了多维度、多相流问题的计算方法[13-17].基于此,本文采用多种LBM模型,考察水泥基材料中的溶液传输问题.

2 水泥基材料微观结构建模

2.1 微观结构特征

水泥熟料为形状不规则,并根据用途和性能要求,具有一定的粒径分布的颗粒材料.水泥颗粒是C3S、C2S、C3A、C4AF等主要成分构成的固溶体,其中C3S占45%～60%.

在微观尺度上的水化硅酸钙(C-S-H)凝胶中,由于颗粒间水分通过扩散作用由外而内发生水化过程,在水化颗粒内外层生成的C-S-H凝胶中存在低密度和高密度组分,且在较大颗粒的内部存在未水化部分[18-21].而根据Jennings等[19]的理论,未水化的C3S平均密度为3.21g/cm3,水化后的C-S-H凝胶密度从1.75g/cm3至2.20g/cm3不等,此外Ca(OH)2密度为2.24g/cm3.综合来看,水化后的水泥基材料可看作密度由内而外逐渐减小的固相颗粒混合而成的多孔复合材料.因此,不能使用一般的多孔材料微观结构生成方法,而应该参考水泥水化过程,将水灰比和水化度作为输入量,建立具有一定密度梯度和孔隙率的水泥微观结构模型.

2.2 自生长法建模

本文提出了一种新的水泥基多孔介质细观模型生成方法,即基于Matlab软件,参考Wang等[22]提出的多孔介质模型生成步骤,改进其中单个颗粒的生成过程为水泥颗粒水化过程.全部步骤(见图1)如下：

(1)将粒径分布曲线离散为各个粒径对应的颗粒数量(图1(a)).

(2)以粒径为限制条件,在随机形状的多边形未水化内核外,生成水泥水化产物(图1(b)).

(3)水化产物的数量达到预设后,得到单个水泥颗粒水化后的模型(图1(c)).

(4)将单个颗粒的水化水泥颗粒模型随机投放至二维区域中(图1(d)、(e)).

图1 水泥微观结构模型生成步骤Fig.1 Cement microstructure model generation steps

下面使用数学描述图1(b)中的单个水泥颗粒的水化过程,即水泥颗粒模型自生长公式：

(3)

式中：di为生长出的第i个微元的密度;fρ没有统一的表达式形式,而是随着模型参数变化,用于控制各等级粒径的颗粒体现出的不同水化情况;rg为最终粒径;rg,min、rg,max为所有粒径的最大值和最小值;imax为待水化的所有微元的总数;KW/C与水灰比成正相关,这里直接取水灰比的值,KW/C=1.0的情况下不进行特殊处理,与Jennings等[19]的文献数据一致,在KW/C<1.0时,近似认为水泥颗粒边缘的最小密度随水灰比线性减小,KW/C=0时,水泥颗粒为未水化状态;Kα与水化度成正相关,这里直接取水化度的值,水化度越大,水化深度越深,在水化度为0时,水化深度为0.本模型中水化度统一取1.0,意为标准养护 28d 后的水泥水化度.

本方法旨在精细控制水泥颗粒的密度分布.式(3)引入三角函数项近似水泥密度的排布,可认为与Fick第二定律的双向扩散模拟结果近似.此外,本方法通过fρ对不同等级颗粒的密度分布进行单独设置.根据陈长久等[23]的研究成果,在28d标准养护后,水泥颗粒水化深度约为3.5μm,因此粒径在7μm以下的水泥颗粒将被完全水化,但由于反应物的浓度梯度仍然存在,这些颗粒还将持续进行生长和扩散.因此参考更大粒径的颗粒密度分布,减小其中央最大密度,并确保最终模型整体密度与水泥石密度仍然一致;而对于粒径在7μm 以上的颗粒,为统一模型最终整体密度与水泥石密度,需要增大大粒径颗粒的水化深度,直到模型整体密度满足要求.

自生长法生成的水泥基材料模型的密度(ρ)分布如图2所示.考虑到投放各等级粒径的颗粒数为正整数，估算模型边界为1420μm×1420μm。不同预设水灰比下，均可得到如图2所示的水泥基材料模型，从而进行扩散、渗流等问题的数值模拟。

图2 水泥基材料模型的密度分布Fig.2 Meso-structure model of cement

2.3 模型参数计算

2.3.1水灰比验证

设置式(3)中的控制参数KW/C,即可调整水泥细观模型的预设水灰比.根据Pommersheim等[24-25]的研究,对于单个水泥颗粒,在原始边界内外的产物体积之比约为1.0∶2.2.因此,对模型中的微元密度di进行统计分析,得到计算水灰比：

(4)

式中：dun为未水化水泥颗粒的密度;nun为模型中未水化部分的微元个数;nhy为模型中水化后的微元个数.

基于水泥颗粒生长参数,随机生成100个水泥基材料二维模型,得到的预设水灰比和计算水灰比及其误差,如表1所示.

表1 预设水灰比和计算水灰比及其误差

由表1可以看出,预设水灰比与计算水灰比的误差较小.因此在后文中使用0.40、0.45、0.50这3个预设水灰比,建立具有代表性的3个水泥模型,进行模拟计算.

2.3.2孔隙率统计

模型生成后,遍历模型中的孔隙,从而得到孔隙率φ,其表达式为：

(5)

式中：nd=0为孔隙微元的个数,即模型中密度为0的微元个数;ntotal为模型整体的微元个数.

使用Image J软件对模型中的毛细孔进行统计,并使用分水岭(watershed)算法对连通的孔隙进行分离,以孔隙的二维面积表征孔隙大小,并根据圆形面积公式反推其等效孔径.各模型的堆积孔隙率与最可几孔径如表2所示.

表2 各模型的堆积孔隙率与最可几孔径

由表2可以看出,当预设水灰比为0.45时,模型孔径显著增大;当水预设灰比为0.50时,模型整体孔隙率φ显著上升.其原因在于,预设水灰比从0.40上升到0.45,对于水泥模型的主要影响为水化度增加,颗粒增大,从而导致堆积孔相应增大;预设水灰比从0.45上升到0.50,水泥颗粒水化度已达到最大,此时水泥颗粒本身变化不大,但多余的水分形成孔隙,导致模型整体孔隙率上升.

3 水分、离子与溶液传输过程的LBM实现

3.1 水分传输过程的LBM模型

气、液两相的物质传输及相变过程,属于单组分多相流模型.对于复杂边界,应使用Shan等[26]的伪势模型进行多相模拟,其演化方程的形式为：

(6)

式中：τ为弛豫时间,与流体黏度有关;ωi为离散速度i的权系数;E为外力项,在边界处即作为边界条件.

平衡态分布函数的形式为：

(7)

式中：ρ为流体的宏观密度;u为流体的宏观速度;cs为格子声速.

本研究采用的离散速度格式为D2Q9模型,即二维模型中具有9个离散速度：

(8)

(9)

流体的宏观密度ρ和宏观速度u确定如下：

(10)

(11)

式中：F(x,t)为流体粒子所受的作用力,主要包括流-流作用力Fc和流-固作用力Fads,其表达式为：

(12)

(13)

式中：Gc、Gads分别代表流-流作用系数和流-固作用系数,决定了表面张力和固相的润湿性;s(x,t)为标示函数,在固体微元处取1,在流体微元处取0;Ψ(x,t)为序参量,表示x位置上t时刻流体的物质状态,出于计算稳定性的考虑,在本研究中取：

(14)

式中：ρ0为临界密度,即气相和液相的分界密度;Ψ0为常量.

流体格点上的压力值P可由非理想状态方程求得：

(15)

3.2 离子扩散过程的LBM模型

当研究钢筋混凝土结构在海工环境下的腐蚀过程时,通常采用氯离子扩散系数予以表征,这其实是假定了混凝土中的孔隙水处于饱和状态,且仅发生环境中的氯离子向内部扩散的过程.因此,仅考虑氯离子的扩散作用,忽略其他作用,LBM中的平衡态分布函数形式如下：

feq=ωiρ

(16)

在模型的半固体格点处需要采用部分反弹格式.经过比较,本研究决定采用Walsh等[27]改良的部分回弹方法.该方法在碰撞步执行前,重定向了需要回弹的速度分量,提高了并行计算的效率,且保证了在单个微元处的质量守恒,具有一定的优越性.其碰撞步演化方程为：

(17)

式中：0≤ns≤1是一个连续的变量,代表微元内传输介质的平均相对密度.

在模型中,水泥的细观结构仅包含毛细孔,而凝胶孔的存在体现为单个微元的密度减小.假设溶液可进入固相网格,其饱和浓度与该网格孔隙度成正比.因此将平衡态分布函数形式修正为：

feq(x,t)=(1-ns)ωiρ(x,t)

(18)

3.3 溶液传输过程的LBM模型

由于离子仅能通过孔隙液进行扩散,因此通过双分布函数的LBM模型进行计算,将水分传输与离子扩散同时进行.根据对流扩散方程(ADE),流场给定的条件下,扩散的平衡态分布函数可取为与浓度、对流速度均线性相关的形式：

(19)

4 结果与讨论

4.1 水分传输过程

使用0.40、0.45、0.50这3个水灰比的模型进行计算.模型左侧区域初始条件设置为液相密度,右侧初始条件为气相密度.为保证计算稳定,设定密度比为20∶1.模型左侧边界与右侧边界设置为恒定压力差以模拟渗透压的作用,均采用Zou-He压力边界条件[28].模型上下侧边界则采用循环边界条件.模型内部的非流体微元统一设置为密度大于1.80g/cm3的微元,避免产生过多的封闭孔隙.不同水灰比的水泥基材料中水分传输过程如图3所示.

由图3可见,在水分传输初期,即1001个时间步左右，液相(图中黑色部分)由于压力差以及毛细作用,迅速充满模型表面附近的孔隙.与此同时，可以看到在模型内部出现局部凝结的水滴.这是气相在蒸汽压作用下向小毛细孔中扩散，随后累积、附着在固相微元(图中灰色部分)表面形成的，随着传输进行到中后期，液相在压力差作用下,优先通过较大的孔隙通道,而后进入较小孔,并与包含凝结水的孔隙连通.模拟结果体现了液相与气相相互转化并共同传输的过程.纵向比较不同水灰比的模拟结果发现,在较高水灰比的模型中,水分传输过程发展较快.

图3 不同水灰比的水泥基材料中水分传输过程Fig.3 Moisture transmission process in cementitious materials with different water-cement ratios

4.2 氯离子扩散过程

仍然选取上述3个模型进行计算.模型的初始条件为：除左侧和右侧设置氯离子浓度C为1.0mol/L外,其他位置均为0mol/L,且无外力作用;模型孔隙中已完成饱水过程,以计算氯离子在多孔介质孔隙水中的扩散过程.模型的左侧及右侧均为恒定浓度,设置方法为借助通量(浓度×速度)守恒得到边界处未知的离散速度分布.模型上侧和下侧为循环边界条件.计算结果如图4所示.

由图4可见,在扩散初期,模型左侧边界附近聚集有一定量的氯离子;随着扩散的进行,边界处的氯离子通过孔隙通道持续向内部扩散,浓度前锋逐渐推进;扩散进行到后期,模型中大部分孔隙被氯离子溶液充满,此时的主导过程为氯离子从右侧溶出,与左侧的补充达到动态的平衡,模型整体的氯离子浓度分布不再发生变化,系统形成稳态扩散.纵向比较3个模型中的扩散过程发现,在扩散的初始阶段,水灰比大的模型扩散更快;随着时间的进行,介质中形成稳定的浓度梯度场,不同模型中的浓度分布趋于一致,均接近线性分布.

图4 不同水灰比模型的氯离子扩散过程Fig.4 Chloride ion diffusion process of different water-cement ratio models

4.3 溶液传输过程

仍选取上述3个模型进行模拟,右侧不设置浓度边界条件,以模拟氯离子与水分同步扩散.上下两侧设置为循环边界条件,其余边界条件相同.计算结果如图5所示.

由图5可见：在水分尚未传输到模型内部时,氯离子的传输主要在边界处进行,且迅速达到饱和状态,与水分传输几乎保持同步;随着时间的延长,水分在模型孔隙中已经以蒸汽的形式先一步到达较大的孔隙,而氯离子只能随液相水扩散,因此落后于水分迁移;参考图3可知,在约10001个时间步后,水分已经基本充满模型内部,而此时溶液中的氯离子仍然在扩散;80001个时间步时,水灰比为0.50的模型率先完成传输过程,达到稳态.

图5 不同水灰比模型的溶液传输过程Fig.5 Solution transfer process of different water-cement ratio models

4.4 扩散系数计算

由于扩散系数的计算需要稳定的扩散浓度梯度场,因此考察饱和介质中的离子浓度场变化过程.以浓度分布的变化率小于0.1%为稳态扩散的标志,得到稳态扩散所需的时间约为500000个时间步长.取第520001个时间步长时沿扩散方向的氯离子扩散速度和氯离子浓度.根据氯离子的有效扩散系数公式,可以计算该模型的氯离子扩散系数(Deff)：

Deff=JL/Δn

(20)

式中：J为稳态扩散中的氯离子通量;L为扩散距离;Δn为模型两侧的浓度差.

不同水灰比模型的氯离子扩散系数的计算结果如表3所示.

表3 不同水灰比模型的氯离子扩散系数

将文献收集到的水泥净浆的氯离子扩散系数[29-31],与本文的数值模拟结果进行对比,如图6所示.

由图6可见：尽管同样是对水泥净浆使用电通量法测定氯离子扩散系数,但不同文献的测定结果仍有一定的差别,这是由材料制备的离散性和电通量法本身的误差造成的.通过对比可以发现,LBM的模拟结果与试验测得的氯离子扩散系数在同一数量级,与特定试验的测量数据接近,同时随水灰比的变化规律近似;在水灰比为0.45和0.50的情况下,模拟结果呈现出与实测结果非常接近的趋势,但在水灰比为0.40时,与实测结果吻合稍差.原因在于,在水化度不变的前提下,水灰比低于一定的数值会导致毛细孔大幅减少,其作为模型渗透性的影响作用减弱.此时,水泥颗粒本身的致密性起到了决定作用.理想模型控制了水泥的水化度不变,在单个颗粒水化完成后再投放进模型,因此,低水灰比的理想模型保留了一定的渗透性.而在实际情况中,室内试验统一采用养护28d的试件.相对于高水灰比的试件,低水灰比的试件中,水泥颗粒水化程度较低,且在水化前就已经处于更为致密的状态,因此抗渗性能较好.

图6 不同文献中水泥净浆的氯离子扩散系数Fig.6 Chloride ion diffusion coefficient in different references

5 结论

(1)水泥基材料水分传输的模拟中,在蒸汽压作用下，液相和气相在材料内部孔隙中传输，并相互连通,达到饱和状态.水灰比大的模型中,水分传输过程较快.

(2)水泥基材料氯离子扩散的模拟中,由于模型两侧浓度差恒定,氯离子浓度前锋不断推进,直到形成稳定的浓度梯度场.水灰比大的模型中,这一过程发生得更早.

(3)在水分与离子共同迁移的溶液传输模拟中,水分的传输过程在时间和空间上都先于离子扩散过程，这是因为离子扩散作用需要在液相环境中才能进行.

(4)使用扩散系数描述氯离子的稳态扩散,在水灰比为0.45、0.50时,模型与试验结果相符,水灰比为0.40时模拟得到的扩散系数偏大,这是由于模型生成时预先统一规定了水化度,而没有考虑低水灰比时水化产物的水化度偏低所致.