考虑微观结构影响的混凝土界面过渡区裂隙渗流-溶蚀耦合模型

王志良，张 跃，申林方，李 泽

(昆明理工大学建筑工程学院，云南，昆明 650500)

水泥浆与集料间的界面过渡区(ITZ)是混凝土材料的最薄弱环节[1]，具有孔隙率高、强度低、Ca(OH)2富集等特点，易成为裂隙的产生和扩展区。对于水工建筑物、隧道及桥梁等实际工程而言，混凝土结构长期与水接触，其中的可溶性钙与水中的Ca2+在浓度梯度的作用下不断地溶解析出，导致混凝土结构孔隙率增大，物理力学特性降低。由于ITZ内的裂隙较为发育且Ca(OH)2含量相对较多，而Ca(OH)2的溶解性又远大于其它水化产物，故ITZ内固相钙的溶解速率较快[2−3]，溶蚀脱钙现象也较为严重。因此，有必要建立混凝土ITZ裂隙的渗流-溶蚀耦合模型，并研究其固相钙的溶解机制，为预测混凝土结构的服役性能提供理论依据。

针对水泥基材料的溶蚀问题，国内外学者开展了大量的研究工作，并取得了丰硕的成果。汤玉娟等[4]采用6 mol/L的氯化铵溶液对水泥砂浆内衬进行加速溶蚀试验，并分析了其微观孔隙结构的变化；方永浩等[5]采用室内试验研究了含裂缝水泥基材料的渗透溶蚀过程，以及溶蚀前后的微观结构变化；王晓梅等[6]研究了在水流作用下水泥基材料开裂表面的溶蚀特性，并建立了相应的Ca2+迁移和表面溶蚀模型；Segura等[7]基于Fick第二定律和未反应缩核模型对水泥砂浆的溶蚀脱钙过程进行了数值模拟；Wan等[8]考虑固液平衡曲线的变化，建立了水泥浆体在6mol/L硝酸铵溶液中的加速溶蚀模型，并采用有限差分法对该模型进行了数值求解。然而，目前关于水泥基材料溶蚀特性的研究大多处于宏观层面试验以及溶蚀前后的微观形貌观测，极少涉及溶蚀过程中微观结构的动态演变。另外，相比于普通的水泥基材料，ITZ的尺度较小，属于微米级别，所以大部分学者主要对其形成原理、微观结构特征及其对宏观力学性质的影响展开研究，如：李冬等[9]建立了考虑混凝土微观组分影响的宏观力学性能理论预测模型，并指出混凝土的单轴抗拉强度受界面过渡区界面特性的影响显著；Leemann等[10]基于微观图像分析，将靠近骨料表面并超过水泥浆体平均孔隙率15%的区域界定为ITZ；Kenny等[11−13]基于微观试验研究指出，ITZ的厚度分布在5µm～100µm；Barnes等[14]认为，随着龄期的增长，ITZ内Ca(OH)2的排布状态由平行于骨料发展为垂直于骨料表面。然而，ITZ的微观试验研究往往难以获取ITZ裂隙溶蚀过程的动态演化，故数值模拟方法成为解决ITZ裂隙渗流-溶蚀问题的重要途径。由于ITZ内部微观结构的复杂性，传统的宏观数值计算方法(有限单元、有限差分等)难以考虑ITZ微观形貌及复杂边界条件的影响，研究Ca(OH)2溶解过程与流体渗流间的耦合作用机制。而基于分子动理论发展起来的格子Boltzmann方法，是一种宏观离散、微观连续的介观模拟方法[15]，在许多传统数值方法难以胜任的领域得到了广泛应用，如多孔介质渗流、多物理化学场耦合作用等[15−18]。

鉴于此，本文在对混凝土ITZ进行微观结构重构的基础上，基于格子Boltzmann方法，采用双分布函数分别模拟流体速度场和溶质浓度场的演化过程，建立了ITZ裂隙渗流-溶蚀耦合模型。根据两个经典算例，分别验证了计算模型在处理溶质对流-扩散及反应-扩散问题方面的有效性。最后，考虑混凝土ITZ微观结构的影响，讨论了不同初始流速、Ca(OH)2含量及Ca(OH)2排布状态等因素对其裂隙渗流-溶蚀耦合作用机制的影响。

1 理论模型

1.1 裂隙渗流控制方程

假设ITZ裂隙中的流动处于层流状态，且流体为不可压缩的牛顿流体，则其流动过程满足：

动量守恒方程：

式中：ρ为流体密度；u为渗流流速；t为时间；P为压力；υ为流体的运动粘滞系数。

1.2 裂隙溶质传输方程

溶质在裂隙中传输主要由对流和扩散控制，则其传输过程满足[19]：

式中：C为溶质浓度；D为溶质扩散系数。

1.3 裂隙表面溶解机制

ITZ的Ca(OH)2含量非常高，且溶解度远大于其它水化产物，故其溶蚀过程主要由Ca(OH)2的溶解造成，其它水化产物的溶解对孔隙率的影响可忽略[20]。为此，本文将Ca(OH)2溶解作为ITZ裂隙溶蚀的主要原因，其溶解过程如下：

假定裂隙溶液中的Ca2+和OH−全部来自于Ca(OH)2的溶解，则OH−浓度C(OH−)为Ca2+浓度C(Ca2+)的2倍。根据一阶动力学反应模型[6,21]，Ca(OH)2的表面溶解过程可表示为：

式中：Cw为固-液界面处的Ca2+浓度；n为由壁面指向流体的法线方向，如图1所示；kd为溶解速率系数；Keq为化学平衡常数，可表示为：

图1 固-液界面溶解反应示意图Fig.1 Diagram of the dissolution reaction at the solid-liquid interface

式中，下标eq表示为离子饱和状态。

在裂隙表面处，Ca(OH)2的溶解反应使其固相体积减小，由式(5)可得到节点体积演化：

式中：Vw为固-液界面处浓度所对应的无量纲体积；A为反应表面面积；M为Ca(OH)2的摩尔体积。

假定初始时刻，Ca(OH)2固体节点的体积为1，则溶蚀过程中其体积不断减小。当固体节点体积减小为0时，将其更新为流体节点。

2 格子Boltzmann模型

2.1 速度场演化方程

采用二维九速(D2Q9)的单松弛格子Boltzmann方程描述流体流动，其演化方程为：

式中：fi(x,t)和fieq(x,t)分别为t时刻、x位置处沿i方向的速度分布函数和平衡态分布函数；∆t为时间步长；τf为无量纲松弛时间；e i为沿i方向的离散速度，其9个方向上的离散速度为：

采用Chapman-Enskog展开，可以由式(8)推导出宏观方程式(1)、式(2)，并得到宏观物理量与格子参数间的联系：

2.2 浓度场演化方程

Ca2+在裂隙中的对流-扩散过程与流速呈线性关系，只需较少的离散速度即可表征其运移规律[19]，故采用二维五速(D2Q5)模型模拟溶质的运移过程，其演化方程为：

同理，由式(12)可推导出到溶质的对流扩散方程式(3)，并得到宏观参数与格子参数间的关系如下：

2.3 边界处理

采用具有二阶精度的半反弹格式，来模拟裂隙渗流过程中流体与固体壁面间的相互作用，其速度分布函数可表示为[15]：

式中：xf为与壁面相邻的流体节点；i为流体节点指向壁面节点的方向；−i为i的反方向。

对于Ca(OH)2表面的溶解反应边界，采用Zhang等[22]提出的半反弹处理格式，其浓度分布函数可表示为：

本文采用压力驱动裂隙内的流体流动，其出、入口处的压力边界采用Guo等[23]提出的非平衡态外推格式，即将边界节点xb的分布函数分为平衡态和非平衡态两部分。平衡态部分fieq(xb,t)可由式(10)求得；而非平衡态部分则由与之相邻的流体节点xbf代替，故其分布函数可表示为[23]：

3 ITZ微观结构及裂隙生成

假定ITZ区域由可溶物(Ca(OH)2)、不溶物(其它水化产物)、孔隙构成。为了建立ITZ裂隙的微观结构模型，在构造区域(L×H)中间设置一条开度为a的裂隙，将整个区域分为ITZ固相区和裂隙区，其示意图如图2所示。本文采用Wang等[24]提出的四参数随机生长法构建其微观结构，实施步骤如下：

图2 ITZ裂隙示意图Fig.2 Schematic diagram of ITZ fracture

1)设计ITZ的构造范围，并划分裂隙区及固相区。在固相区内将Pc作为Ca(OH)2初始固体节点的分布概率，随机生成初始固相生长核，其分布概率Pc不得超过Ca(OH)2的体积分数。

2)以Ca(OH)2初始固相生长核作为中心，根据不同方向的生长概率Pi向相邻的8个节点随机生长，如图3所示。

图3 随机生长方向示意图Fig.3 Diagram of random growth direction

3)重复步骤2)的操作，直至Ca(OH)2的固体体积分数达到设计值。

4)在固相区的非Ca(OH)2范围内，以一定的概率随机生成不溶物的固相结构，直至体积孔隙率达到设定值。

5)在构造范围内生成完可溶物(Ca(OH)2)、不溶物(其它水化产物)后，剩余的节点即为孔隙。

假定ITZ的构造范围为L×H=200×50格子，在中心处设置一条裂隙，其开度为a=20个格子，固相区的孔隙率为ϕITZ=0.15，Ca(OH)2体积含量为35%。采用上述方法，分别构建了Ca(OH)2水平生长、均匀生长和竖向生长的ITZ的微观结构及相应裂隙。为了避免后续数值计算中，边界处固相对流体流动的影响[19]，在裂隙两端各增设了5个格子的孔隙，如图4所示。图中灰色物质为Ca(OH)2固体，黑色为不溶物，白色为孔隙。

图4 Ca(OH)2不同生长方式的ITZ裂隙Fig.4 ITZ fractureswith different growth patterns of Ca(OH)2

4 模型验证

4.1 半无限区域内溶质的对流-扩散问题

为了验证本计算模型处理溶质对流-扩散问题的准确性，在半无限区域内建立了相应的计算模型。假定初始时刻计算域内的溶质浓度为0，流速恒定为u。随后从左侧注入浓度为C0的溶质，则一维的溶质对流-扩散方程可表示为：

针对L×H=200×10的计算网格，采用D2Q5模型进行了溶质对流-扩散的数值计算。为了便于分析，计算参数均采用格子单位。假定入口处的溶质浓度为C0=1.0，渗流流速为u=0.002，溶质的扩散系数为D=0.288。图5为本文数值解与解析解的对比，从图中可以看出，本文的数值解与解析解具有较好的一致性，这充分说明本计算模型处理溶质对流-扩散问题的有效性。

图5 对流-扩散问题溶质浓度的数值解与解析解对比Fig.5 Comparison between numerical and analytic solutions of soluteconcentration for convection-diffusion problem

4.2 矩形区域内的反应-扩散问题

为了验证本计算模型处理溶质反应-扩散问题的有效性，对于稳态下矩形区域内的反应-扩散问题进行了数值计算。假定计算域内的浓度分布由Laplace方程来描述，即：

计算模型如图6所示，假定在上边界处(y=H)发生一阶化学反应；左边界(x=0)注入恒定浓度为C0的溶质；右边界(x=L)和下边界(y=0)均为零通量边界，其边界条件如下：

图6 稳态溶质反应-扩散问题示意图Fig.6 Diagram of solutereaction-diffusion problem in steady state

在数值计算中，本算例的计算参数均采用格子单位。计算网格为L×H=200×50，假定初始时刻溶质浓度处于平衡状态，C=Ceq=1.0，左侧入口处的浓度C0=10.0，溶质的扩散系数为D=0.288。针对化学反应速率分别为kr=0.01、kr=0.05及kr=0.1 这3种情况，基于本文数值计算方法，求解了稳态溶质反应-扩散问题，其浓度场的本文数值解与解析解的对比，如图7所示，图中实线为解析解，虚线为本文数值解。由图可知，本文数值解与解析解的吻合度非常高，两者间的最大相对误差仅为0.35%，这说明本文计算模型处理溶质反应-扩散问题的准确性。

图7 反应-扩散问题溶质浓度的数值解与解析解对比Fig.7 Comparison between numerical and analytic solution of soluteconcentration for reaction-diffusion problem

5 分析讨论

为了研究ITZ裂隙的渗流-溶蚀耦合作用机制，建立了如图2(示意图)及图4(微观结构图)所示的计算模型。将ITZ裂隙的计算域L×H=1000μm×250μm 划分200×50的网格。裂隙开度为a=100µm，模型中的计算参数取值[26]如表1所示。设在初始时刻裂隙内的Ca2+浓度处于饱和状态，即C=Ceq，然后在左侧入口处注入Ca2+浓度为0的软水，同时在压力驱动下裂隙发生渗流，并在固-液界面处发生Ca(OH)2的溶蚀现象。整个ITZ裂隙的初始孔隙率为ϕ=0.515，随着溶蚀的发生，其孔隙率逐渐增大。为了分析ITZ裂隙渗流-溶蚀间的耦合作用机制，采用渗透率表征其演化过程，相应的计算公式为：

表1 模型参数Table 1 Model parameters

式中：Q为裂隙面渗流量；dP为压力差。

在初始时刻ITZ裂隙的渗流为Poiseuille流，其平均流速为：

5.1 渗流流速影响

为了研究渗流流速对裂隙渗流-溶蚀耦合作用机制的影响，针对Ca(OH)2体积含量为35%且均匀生长的ITZ裂隙，当初始平均流速u¯0分别为4.167×10−6m/s、8.333×10−6m/s、1.250×10−5m/s时，进行了其渗流-溶蚀过程的数值计算。图8为裂隙的整体孔隙率与相对平均流速的时程演化曲线。由图可知，在初始阶段由于软水的突然注入导致Ca2+浓度显著下降，使得壁面处的溶蚀速度加快。随着时间的推移，裂隙内的Ca2+浓度逐渐趋于稳定，溶蚀速率逐渐减慢，其整体孔隙率和相对平均流速的变化趋势也逐渐变缓。同时，对于初始平均流速较大的裂隙，由于其Ca2+浓度的传输速率较快，溶蚀过程迅速趋于稳定，故受隙宽影响的相对平均流速也表现出相似的趋势。总体而言，渗流流速越快，裂隙的溶蚀速率越快，整体孔隙率和相对平均流速的变化趋势也相应的越快。此外，图9为溶蚀2.5×104s后裂隙中线处Ca2+浓度与渗流流速的分布图。由图可知，在裂隙中线位置处，Ca2+浓度随着渗流流速的增加而减小，且在入口附近浓度变化较快，变化趋势随着向裂隙内部的推进而逐渐变缓。同时由于受到Ca(OH)2溶蚀裂隙壁面微观结构变化的影响，渗流流速沿着裂隙长度方向并非均匀分布，而呈现出不规则的波动，同时流速的波动也对Ca2+浓度产生了一定的影响，在流速变化较大的位置处，Ca2+浓度的波动也较为剧烈。

图8 整体孔隙率与相对平均流速的时程曲线Fig.8 Timehistory curvesof overall porosity and relative average flow velocities

图9 t=2.50×104 s时裂隙中线处Ca2+浓度与渗流流速分布Fig.9 Distribution of Ca2+concentration and velocity at the central part of fracture at t=2.50×104 s

5.2 Ca(OH)2含量影响

为了分析Ca(OH)2含量对ITZ裂隙渗流-溶蚀特性的影响，在初始平均流速u¯0为4.167×10−6m/s时，针对Ca(OH)2均匀生长且体积含量ω分别为35%、45%及50%这3种工况下的渗流-溶蚀过程进行了数值计算。图10为不同Ca(OH)2含量情况下裂隙相对渗透率随时间的变化趋势，从图中可以看出，在溶蚀初期，由于注入软水使得裂隙内的Ca2+浓度显著下降，溶蚀速率加快，相对渗透率的增长速率也较快。随着时间的推移，裂隙相对渗透率的变化趋势逐渐趋于稳定，且K/K0(ω=50%)>K/K0(ω=45%)>K/K0(ω=35%)。图11为裂隙中截面(x=L/2)不同位置处的Ca2+浓度的时程曲线。由图可知，距离壁面不同位置处Ca2+浓度的时程变化具有相似的变化规律，仅仅在壁面附近的浓度稍大。在初始阶段，受注入软水的影响ω=35%的ITZ裂隙中Ca2+浓度下降幅度最大、ω=45%次之、ω=50%最小，主要是由于Ca(OH)2含量较高的ITZ裂隙，其Ca(OH)2溶蚀速率较快，致使裂隙内部的Ca2+浓度相对较高。随着裂隙渗流-溶蚀耦合作用的进行，3者的浓度变化趋势渐趋稳定，且稳定后的Ca2+浓度变化规律为C(Ca2+)(ω=50%)>C(Ca2+)(ω=45%)>C(Ca2+)(ω=35%)。总体而言，3种Ca(OH)2含量情况下，裂隙溶蚀演化特征的差异主要是由于Ca(OH)2含量越高，壁面处其与流体的接触面积越大，发生溶蚀的概率也相应增加，裂隙内Ca2+浓度分布也就越高，如果Ca2+得不到有效的运移，它又反过来降低Ca(OH)2的溶蚀速率。故裂隙的相对渗透率受控于Ca(OH)2含量与Ca2+浓度的综合作用，即，裂隙内的Ca2+能够得到及时运移，Ca(OH)2溶解加速，则ITZ内Ca(OH)2含量即使较高，也能够快速溶解；反之，则会抑制Ca(OH)2溶解。

图10 不同Ca(OH)2含量时相对渗透率的时程曲线Fig.10 Time history curvesof therelative permeability under different Ca(OH)2 contents

图11 裂隙中截面不同位置处的Ca2+浓度时程曲线Fig.11 Time history curvesof Ca2+ concentration at different positionsof thecentral section of fracture

5.3 Ca(OH)2排布状态影响

为了研究Ca(OH)2排布状态对ITZ裂隙渗流-溶蚀特性的影响，针对如图4所示的Ca(OH)23种排布状态(水平生长、均匀生长和竖向生长)，在初始平均流速为u¯0=8.333×10−6m/s时，进行了裂隙渗流-溶蚀过程的数值计算。图12为3种Ca(OH)2排布状态下ITZ裂隙整体孔隙率与相对渗透率的时程演化曲线。从图中可以看出：在溶蚀前期，Ca(OH)2水平生长整体孔隙率变化相对较小，也就是溶蚀过程相对较慢。这主要是由于渗流初期水平生长的Ca(OH)2与流体接触面积较大，致使其溶解释放出大量的Ca2+(如图13所示)，从而抑制Ca(OH)2后续的溶解反应。随着裂隙渗流过程的进行，Ca(OH)2水平生长的裂隙，其整体孔隙率发展较快，并逐渐超越另外2种排布状态。这主要是由于Ca(OH)2水平生长与裂隙渗流方向一致，溶解后其连通性较好，致使相对渗透率发展较快(如图12(b)所示)，从而更利于Ca2+的运输及溶蚀的发展。同时，裂隙的相对渗透率在溶蚀早期主要受Ca(OH)2溶解释放出的Ca2+影响，待稳定后的3种排布状态的规律为：水平生长最大，均匀生长次之，竖向生长最小。从渗流场的微观流线图(图14)来看，Ca(OH)2水平生长裂隙的内部流线迂回度最小，而竖向生长则由于固相不溶物的影响，其迂回度最大，从而影响了裂隙整体的渗流流速和溶蚀进程。故Ca(OH)2微观结构上的排布状态，是影响ITZ裂隙渗流-溶蚀特性的一个重要因素。

图12 不同Ca(OH)2排布状态时整体孔隙率与相对渗透率的时程曲线Fig.12 Time history curves of overall porosity and relative permeability under different arrangements of Ca(OH)2

图13 t=3.0×104 s时不同Ca(OH)2排布状态裂隙中线处的Ca2+浓度分布Fig.13 Ca2+ concentration distribution at the central part of fracture under different arrangementsof Ca(OH)2 at t=3.0×104 s

图14 溶蚀后ITZ裂隙局部渗流流线图Fig.14 Local streamline diagram of ITZ fracture after dissolution

6 结论

基于格子Boltzmann方法，建立了考虑微观结构影响ITZ裂隙渗流-溶蚀耦合演化过程的数值模型，并结合算例验证了本数值模型在处理对流-扩散及反应-扩散问题方面的准确性。在此基础上，讨论了ITZ裂隙不同初始渗流流速、Ca(OH)2含量及Ca(OH)2排布状态等因素，对其裂隙渗流-溶蚀耦合作用机制的影响，得到以下结论：

(1)ITZ裂隙的初始渗流流速越快，壁面的溶蚀速率越快，其整体孔隙率增大越快。而流体中的Ca2+浓度随着渗流流速的增加而减小，且在入口附近浓度变化较快，其变化趋势随着向裂隙内部的推进而逐渐变缓。

(2)Ca(OH)2含量越高，在壁面处其与流体的接触面积越大，溶蚀现象越易发生，导致裂隙内Ca2+浓度也相应的增加，这又反过来抑制Ca(OH)2的进一步溶蚀。故ITZ裂隙的渗流-溶蚀过程受控于Ca(OH)2含量与Ca2+浓度的综合作用。

(3)对于Ca(OH)2不同排布状态的ITZ裂隙，待溶蚀稳定后其相对渗透率为：水平生长最大，均匀生长次之，竖向生长最小。这主要是由于Ca(OH)2水平生长裂隙的内部流线迂回度最小，而竖向生长则由于固相不溶物的影响，则迂回度最大，从而影响了裂隙整体的渗流流速和溶蚀进程。