基于扩散路径映射的混凝土碳化细观模拟方法

李 越， 阮 欣， 尹志逸

（1.同济大学 土木工程学院，上海 200092；2.上海市城市建设设计研究总院（集团）有限公司，上海 200092）

近年来混凝土结构病害集中爆发，其中最为普遍的是混凝土碳化问题［1-2］。结构长期服役过程中，大气环境中的二氧化碳通过混凝土空隙进入结构内部，削弱钢筋表面碱性环境，诱发钢筋锈蚀膨胀和局部保护层开裂，进而导致结构病害和威胁服役安全［3-4］。实验室环境中通过材料加速碳化测试逐渐建立了考虑不同环境影响因素的经验公式，并在结构设计规范中不断更新和细化，但仍然难以解释目前大量结构局部随机退化的现象［5］。经典的复合材料中，由于骨料颗粒和水泥基体随机分布，碳化作用的空间分布也极为复杂［6］。传统的试验拟合公式方法，时间周期长、测试成本高，且难以反映碳化作用的随机性。基于碳化机理和材料特性的数值模拟方法在混凝土材料变异性的精细模拟方面优势突出［7］。

区别于传统混凝土宏观模拟中基于均质化假定对材料单元赋予统一的单元属性，混凝土细观建模方法将材料模拟为包含骨料、水泥浆体和界面层（interfacial transition zone， ITZ）的三相复合模型［8］。基于王宗敏［9］的随机骨料结构（random aggregate structure， RAS）模型，细观模型主要通过骨料随机生成与投放实现考虑不同粒径、形状、棱角的骨料颗粒空间随机分布建模［10］。基于模型单元所处位置的材料组分赋予对应物理化学属性，以实现混凝土材料内部响应和反应的场域模拟。采用与宏观模型一致的模拟原理和材料本构关系，混凝土细观模型在力学响应、温度传导、物质扩散等方面都能够实现精细化的仿真分析。混凝土碳化作用的细观模拟中，通过材料配比和水化参数确定材料基础参数，通过Fick扩散定律计算各关键物质的场域动态分布，基于化学反应动力学原理控制碳化反应速率参数。潘子超等［11］在细观层面建立混凝土碳化过程数值模拟方法，通过三维球形骨料模拟了混凝土内部非均匀碳化问题，研究表明骨料对于混凝土碳化影响显著。Ruan等［12］对比了混凝土结构碳化模拟中宏观与细观模拟差异，讨论了不同混凝土骨料参数影响下的碳化深度分布平均水平与变异程度。Jiang等［13］也在混凝土材料碳化试验的基础上模拟了不同骨料形状参数对碳化深度的影响，研究表明，随着骨料长细比的增大，碳化深度的变异性显著上升。近年来随着基础算法与建模能力的提升，混凝土材料精细化模拟在材料尺度模拟方面已经取得了长足的进步，但是精度提升带来的庞大模型矩阵求解也限制了该方法在实际工程尺度的应用。因此为了精准预测混凝土结构碳化问题的空间随机特性，开发高效精准的混凝土细观碳化快速计算方法十分必要。

在现有混凝土材料细观碳化模拟研究的基础上，本文建立了基于扩散路径进行映射的混凝土结构层面细观碳化作用的快速模拟方法，简称为扩散路 径 映 射（diffusion path dependent mapping， DPDM）方法，通过将结构细观模型扩散路径计算以及水泥浆体碳化作用求解相结合，实现了细观模型精度和宏观模型效率的优势互补。与传统基于偏微分方程（partial difference equation， PDE）求解的碳化作用模拟方法相比，该方法保持了高精度的局部碳化作用模拟，同时大幅提升了模拟效率并降低了结果存储的空间需求。

1 碳化作用细观映射原理

混凝土细观碳化过程可以分解为水泥浆体中矿物成分水化过程以及水化产物的碳化过程，其反应原理与现有宏观模拟方法相似，关键的差异在于细观尺度材料介质对于求解域Ω的影响。矿物成分水化的化学反应过程主要包括硅酸二钙（2CaO·SiO2， C2S）、硅酸三钙（3CaO·SiO2， C3S）、铝酸三钙（3CaO·Al2O3， C3A）以及铁铝酸四钙（4CaO·Al2O3·Fe2O3， C4AF）消耗孔隙水分形成水合硅酸钙（3CaO·2SiO2·3H2O， CSH）凝 胶 以 及 氢 氧 化 钙（Ca(OH)2，CH）晶体。这些反应主要发生于水泥浆体与ITZ介质中，与混凝土中骨料组分并不发生反应，其具体反应过程如下式所示：

式中：Ri，hydr为混凝土矿物成分i细观尺度的水化机理作用，考虑到骨料介质不参与水化反应作用，因此细观尺度的水化反应空间求解域为模型内部的水泥浆体范围Ωcm，是混凝土模型范围Ωcon的子集，在本研究中也将作为细观模型路径影响抽离的关键参数；fi，hydr为矿物成分i的水化反应速率。

在不考虑矿物成分在水泥浆体中微观不均匀分布的情况下，由式（1）可以发现，水化反应的进程主要受化学反应原料（水、矿物成分）的空间浓度以及水化速率控制。水化速率可以通过Papadakis等［14］提出的水化动力学模型计算，而混凝土中水分随时间空间分布变化，可以通过下式的Fick扩散模型进行计算：

式中：JH2O为混凝土中的水分传输通量，由于骨料致密可以视为不可渗透介质，因此水分扩散过程也局限于水泥浆体求解域Ωcm，将细观介质空间分布以求解域的形式进行下标表征，可以有助于后续对相同求解域进行合并处理，以降低水泥浆体空间部分复杂影响的重复计算；DH2O，cm为水泥浆体中的水分扩散系数；cH2O为水泥浆体中的含水量，∇为用于求解场域梯度的Laplace算子。

水分在环境表面的传输渗透影响材料孔隙相对湿度，涉及水泥浆体中水分的吸附与解吸附过程，具体计算方法可以参照潘子超等［12］相关研究。综上可以建立混凝土内部水化过程的系统矩阵如下式所示：

混凝土碳化过程主要为大气环境中的二氧化碳扩散溶解于孔隙溶液，水化产物CSH凝胶、CH晶体以及部分矿物成分C2S和C3S与孔隙溶液中的二氧化碳发生反应，生成碳化产物碳酸钙。其具体反应过程涉及二氧化碳气体溶解、氢氧化钙溶解，酸解中和反应以及碳酸钙生成，如下式所示：

式中：g、aq以及s分别表示化学成分的气态、溶解态以及固态状态。

根据化学反应表达式可以发现，碳化反应进程的主要控制因素包括作为反应原料的水化产物、二氧化碳、水分的空间分布，二氧化碳的渗透扩散，二氧化碳的气态与溶解态平衡，以及二氧化碳的反应速率。其中二氧化碳的扩散过程同样可以采用式（5）以Fick扩散定律为基础的方法进行模拟；气态与溶解态二氧化碳的转换可以通过Herry定律建立浓度平衡方程，而碳化反应进程控制的关键指标为反应速率fj，carb，该指标与反应原料浓度有关，同样可以依据Papadakis等［14］的相关研究进行模拟。

式中：JCO2为混凝土中的二氧化碳的传输通量，同样由于骨料致密不可渗透介质，因此二氧化碳扩散过程也局限于水泥浆体求解域Ωcm；DCO2，cm为水泥浆体中的二氧化碳扩散系数；cCO2为水泥浆体中的二氧化碳浓度。

对碳化过程Rcarb进行总结，如式（6）所示。式（6）中箭头左侧为模拟需要考虑的机理效应包括化学反应过程、反应原料、反应速率系数、二氧化碳渗透以及气态溶解态平衡方程；箭头右侧为模拟输出的相关结果，包括碳化反应消耗后的原料浓度、水分以及二氧化碳浓度的分布情况。

在混凝土碳化水化和碳化过程中另一项材料参数的场域变化也同样需要考虑水泥浆体的孔隙率指标。随着水化和碳化产物的生成，水泥浆体内部的孔隙率逐渐下降，并影响水分和二氧化碳的扩散系数，由于孔隙率变化同样局限于水泥浆体内部，因此孔隙率变化过程φ可以表述为

式中：φ为混凝土孔隙率变化，其变化过程如式（8）所示，由初始孔隙率φ0受水化反应影响φhydr以及碳化反应影响φcarb，得到最终孔隙率在水泥浆体范围中随时间变化。其中水化作用和碳化作用的具体影响数值可以通过现有研究进行模拟计算。

综上，对于混凝土细观碳化问题的模拟可以表述为水化反应、碳化反应以及材料孔隙率变化过程的方程联列。如式（9）所示。联列过程Rhydr⊕Rcarb⊕Rφ均发生于服役时间中水泥浆体的求解域Ωcm，t，对求解域进行合并，即可表述为统一范围内的多个反应耦合求解。对比混凝土碳化过程的传统宏观模拟方法，多反应同步的模拟过程在时间范围Ωt内的求解并不存在差异，因此对联列方程的求解域Ωcm，t进一步分解得到Ωt和Ωcm。可以发现，与现有宏观模拟方法相同，而[⋮]Ωcm则对求解的空间范围提出了复杂的要求，需要对混凝土模型内构建水泥浆体网络结构，并在其中开展化学反应过程求解。该过程对于建模精度提出了很高的要求，且消耗大量的计算成本。

从目前碳化反应模拟研究来看，细观尺度求解域对于反应机理作用并不构成影响，因此本研究对求解空间域Ωcm的影响进行分开考虑。将传统宏观模拟结果反映到细观水泥浆体网络中，关键问题在于反映骨料随机分布对于侵蚀扩散路径的影响，而上述结论表明，碳化过程主要涉及环境物质的Fick扩散渗透以及各类反应原料的消耗和产物的生成。在不考虑化学成分浓度堆积的影响下，Fick侵蚀扩散路径一般符合最短路径原则，因此可以对细观求解域[⋮]Ωcm进行路径求解，同时在宏观求解过程输出路径长度信息，最终基于相同路径长度实现混凝土细观碳化模型的时间空间域求解，路径求解用⊕path表示。

2 扩散路径映射方法

2.1 方法概述与求解步骤

按照前述方法对混凝土细观碳化模拟过程进行拆解，骨料颗粒空间随机分布形成的水泥浆体网络以及碳化过程的迂曲分布与反应扩散的机理无关，可以通过计算机路径求解算法直接计算路径网络；而细观碳化过程中物理化学反应机理与宏观方法本质是统一的，不考虑路径空间分布对水化碳化反应进行精细求解也是可行的。面向工程分析需求，式（10）中的细观碳化模型拆解公式并不能直接指导具体计算仿真，因此将模型计算转化为关键物质浓度随时间空间分布的求解。

如图1所示，细观路径求解过程通过路径求解算法得到水泥浆体内部节点到环境表面的最短路径长度，基于细观模型中材料组分划分，路径计算可以表述为式（11）中的函数形式LMeso(X)。其中输入变量X为表征求解模型中节点空间位置向量，取值范围为模型中水泥浆体求解域Ωcm；输出则为各节点到环境表面的传输路径长度。

针对水泥浆体中的水化反应、碳化反应以及孔隙率变化过程，通过图1右上方的代理模型进行求解。通过较大的网格划分和简化的条带状模型设计，实现对于各阶段关键物质浓度或者指标随深度和时间的分布，该过程可以表述为式（12）的函数形式cMacro，i(l，t)。其中输入变量为模型深度l和求解时间节点t，求解域Ωt为研究所需的时间范围，输出为关键物质i的状态指标。

图1 混凝土细观碳化DPDM方法模拟流程Fig.1 DPDM simulation procedure for mesoscopic concrete carbonation

基于两个模型中一致的路径长度进行结果映射，其函数意义可以表述为式（13）。将式（13）与式（10）～（12）联列，可以将求解过程转化为式（14）中的复合函数形式，细观模型路径长度函数值作为代理模型中的函数输入变量，最终得到模型中关键物质i浓度随时间空间分布结果cDPDM，i(X，t)。

该模型的求解主要包含3个步骤：①在结构截面内建立骨料颗粒空间分布，并映射到截面网格单元中，采用节点网络最短路径计算方法求解扩散路径；②基于水泥浆体材料属性建立代理模型，基于现有研究基础模拟水化、碳化的物理化学反应过程，获得关键物质在深度分布的时变结果；③根据具体分析需求，基于同样侵蚀扩散深度，将特定时间节点的关键物质浓度映射到结构截面内，即可得到关注时刻的结构全截面内碳化作用关键物质的场域分布特征。得益于模型求解中对于结构特征和复杂材料环境特征的分离处理，该模型对于真实混凝土结构和复杂环境条件均能够提供高效精准的解决方案。

2.2 细观扩散路径计算

混凝土结构细观碳化路径的求解主要包含以下流程：结构关键截面的选取，截面骨料颗粒建模，以及细观最短扩散路径计算。

研究表明，混凝土材料内骨料颗粒的几何形态和随机空间分布对于局部碳化效应影响显著，因此在截面细观建模中需要考虑实际工程材料采用的骨料配比参数、类型以及浇筑方案的影响。本研究采用基于真实桥梁结构混凝土材料截面图像进行改进的三参数骨料模型［10］，该模型通过随机参数抽样实现了颗粒粒径、长宽比以及棱角细节的精准还原。

截面内细观扩散路径的求解建立在模型网格划分的基础上，对模型按照所需模拟的精度需求进行自然网格划分，按照单元节点与细观模型几何位置的关系将节点分类为骨料节点、水泥浆体节点以及ITZ节点。在可以侵蚀退化的水泥浆体和ITZ节点中，采用计算机中常用的Dijkstra最短路径算法进行路径计算［15］。具体计算过程是将结构表面与环境直接接触的节点设置为路径起点（扩散距离为0），对于内部节点的扩散长度按照周边关联节点的扩散长度最小值进行更新，并逐步迭代直至所有节点扩散路径长度降至最低并稳定，即可得到水泥浆体中所有节点到结构表面的最短距离。

2.3 碳化作用代理模型求解

通过代理模型对服役过程中水泥浆体物理化学反应过程进行模拟，主要包括反应传输的传输扩散模拟，和反应原料消耗及产物生成的计算求解。该部分算法的具体实现形式是在单个时间步中，基于Fick扩散定律的偏微分方程对自由物质（水分、二氧化碳）分布进行求解，并对其他物质（矿物成分，水化产物、剩余水分、）的浓度和指标（饱和度、孔隙率）的状态进行更新。

区别于传统混凝土碳化数值模拟，本方法中水泥浆体模型仅作为碳化作用物理化学反应模拟的代理模型，与真实结构形式和尺寸并不直接相关，只需要满足结构退化深度方向上的模拟需求即可。模型可以通过建立狭长的代理模型大幅降低求解难度。一般模型宽度方向尺寸只要满足2倍以上于水泥浆体单元尺寸即可；模型长度尺寸需要略大于截面各点到环境表面的深度以满足后续映射的需求，考虑到骨料对于物质传输路径的迂曲影响，一般可以设置该代理模型长度为结构截面深度的1.2倍以上，其代理模型求解和结果提取如图2所示。

图2 水泥浆体代理模型碳化模拟结果示意图Fig.2 Illustration of carbonation simulation results of cement mortar agent model

2.4 碳化结果空间映射

扩散路径映射方法中最终的碳化作用全截面模拟，是基于同样的扩散路径长度和所关注的时间节点，将代理模型中的关键物质参数沿深度分布结果映射到结构整体截面中，即可得到细观碳化作用在结构全截面的分析结果。

（1）对于真实混凝土结构，由于混凝土细观截面和扩散路径求解均为一次性计算，且骨料颗粒建模和Dijkstra运算的时间和存储消耗均与模型尺寸正相关，因此混凝土结构碳化分析在面临众多截面分析时，该部分的计算成本并不会随模拟面积呈现指数增长。

（2）对于基础设施面临的复杂多变服役环境，不同时间尺度的物理化学反应则通过代理模型进行，因此能够考虑较为精细的环境气候变化影响，且不会造成模型求解和结果存储的负担。

2015年8月至2017年6月我院脊柱外科采用PVP治疗椎体压缩骨折患者105例,其中男性19例,女性86例,年龄45岁至97岁,平均69岁;所有患者术前均行MRI检查,提示新鲜骨折椎体132个,所有患者病椎术中均采用负压状态下同轴取材法采集标本,所采集标本进行病理检查,取材后常规行椎体成形术。

3 试验校验与模拟分析

为了校验DPDM方法的模拟准确性以及对比传统方法在计算精度和求解效率上的提升，本研究开展了两方面的对比分析：①将DPDM方法与混凝土试块碳化试验进行对比，通过对比碳化深度预测结果说明该方法的模拟准确性；②将DPDM方法与细观偏微分求解方法进行对比，通过对比碳化反应关键物质浓度分布以及计算成本说明该方法在求解性能方面的优势。

3.1 试验校验

混凝土试块碳化试验采用边长100 mm的标准立方体，使用工程C50混凝土材料配合比设计（表1），碎石骨料按照Fuller级配进行混合，粒径范围为0.5 ～ 25.0 mm，含量为200 kg·m-3。材料试块在浇筑成型后经2 d干燥，之后对两个侧面进行环氧封闭，放入碳化试验箱进行为期60 d的碳化作用试验，试块周围二氧化碳浓度为8.929 mol·m-3，环境温度和湿度相对稳定。

表1 混凝土配合比参数Tab.1 Mixing proportion parameters of concrete kg·m-3

在DPDM数值模拟中，采用同样的材料参数进行混凝土细观建模，具体多边形骨料的形状指标控制采用Ruan提出的多重抽样细观建模方法。考虑到ITZ薄层厚度为30 ～ 80 μm之间，模型采用100 μm的网格划分能够精细模拟ITZ界面层以及扩散路径。而对于混凝土各相组分在水化、碳化过程中物理化学属性通过砂浆试块碳化测试以及相关研究结论推算得到。代理模型求解过程中，水化过程求解以小时为时间步长，碳化过程求解按照天为时间步长。

试验碳化作用结束后，劈裂试块，对新鲜截面进行酚酞显色反应以及磨粉化学滴定。得到截面中关键化学成分CH的浓度随深度分布情况，对比DPDM方法映射得到的CH浓度分布，其中碳化结果空间分布情况如图3所示。模拟结果中环境物质从试块表面侵入，并围绕骨料颗粒向内扩散，碳化侵入平均水平与细观变异性均与真实试验结果较为接近。对模拟结果最终CH浓度随深度变化曲线进行提取，并与试块磨粉滴定结果进行比对，如图3c所示，浓度分布曲线与实测结果趋势接近，碳化影响临界深度也基本一致。浓度分布曲线在临近试块表面存在少量偏差，其中实测结果呈现局部浓度偏差并不符合梯度扩散的原理，这也说明了当下传统试验测试手段指导细观数值模拟的精度不足。

图3 DPDM模拟结果试验对比校验Fig.3 Experimental comparison and verification of simulation result of DPDM

3.2 数值模拟分析

为了校验DPDM方法相较于现有混凝土细观碳化过程模拟的优势，本研究采用与前文试验一致的模型尺寸和材料性能参数。其中细观模型同样采用100 μm的网格划分，针对为期2 d的干燥水化过程采用小时为时间步长，针对为期120 d的碳化过程采用天为时间步长。将传统细观模拟和DPDM模拟结果进行对比，可以得到不同水化和碳化阶段关键化学物质在截面内的浓度分布结果。如图4左侧两列所示，上部为水化过程中水泥浆体内部饱和度分布，从深到浅的颜色反应代表饱和度从100%至0，其干燥水化阶段混凝土内部水分从暴露表面流失，而呈现试块内部饱和度的场域变化；下部为碳化过程水泥浆体中CH浓度分布变化过程，深色代表水化产物中饱和CH溶液浓度而浅色表示碳化反应消耗殆尽的CH浓度，其碳化影响范围随着二氧化碳从暴露表面渗入混凝土孔隙，并逐渐向内扩散推进引发化学反应。两阶段中材料变化过程差异显著，水分扩散过程相对快速且均匀，二氧化碳扩散反应过程则较为缓慢且呈现显著的CH浓度突变锋面。DPDM方法对于水分扩散过程和碳化反应过程的场域演化的特性能够较为精准地模拟，对于具体的侵蚀范围和深度水平的预测结果也与广泛接受的细观模拟结果较为接近。从饱和度和CH浓度分布的云图来看，两种方法结果仅在大颗粒骨料之间狭窄扩散通道以及侵蚀锋面位置存在细微差异。

通过对同样位置的计算结果差异绝对值进行计算，对比结果之间的具体差异以及偏差主要发生的位置特征，具体计算方法如下式所示：

式中：e(X，t)为两个模型在局部位置X、时间t的结果偏差；cDPDM(X，t)为DPDM方法计算得到CH浓度随时间空间分布规律；cPDE(X，t)为偏微分方程求解方法得到模拟结果。

对比得到相对于传统细观方法的计算偏差分布，如图4右侧一列所示，颜色越深代表计算结果的偏差越大。干燥水化阶段，水泥浆体饱和度计算偏差小于10%，主要存在于大颗粒骨料对于水分扩散造成的阻碍堆积效应，以及颗粒间狭窄通道位置水分扩散的细微计算偏差。碳化作用阶段，DPDM方法模拟得到的CH浓度分布情况与传统细观模拟结果也基本一致，其在整体碳化深度以及骨料阻碍效应的模拟方面都较为精准，与现有方法的差异主要位于碳化封面的浓度突变相对平缓，浓度变化位置达到1 mm左右，这从图4最右列的偏差中也可以发现。导致该问题的原因在于代理模型网格稀疏，导致浓度深度曲线较为平缓，浓度突变结果也较为平缓。

图4 传统PDE方程求解结果与DPDM模拟结果对比分析Fig.4 Comparison and analysis of solution result of traditional PDE and simulation result of DPDM

基于单个细观模型可以获得混凝土材料内部非均匀碳化的局部差异，骨料颗粒导致的侵蚀路径迂曲导致了最终材料内部碳化影响深度的变异特征，这也是混凝土结构层面的局部劣化随机涌现的根源。依据上述细观模型的DPDM模拟方法，从单个模型中即可建立完整的碳化深度概率密度分布模式，这是传统简化的宏观均质模型所不具备的。

为比较DPDM方法对求解效率提升和存储空间缩减的实际效果，本研究基于同样的计算软件平台（Matlab R2020a）和硬件配置（Intel I7 7900， 32 GB），对方法计算效率和存储空间需求进行对比分析。对5个混凝土截面模型求解，所得参数如表2所示。由表2可知，PDE求解耗时约为70 min，而DPDM方法求解耗时普遍小于2 min，时间缩减97%左右。其中PDE计算耗时主要用于模型整体偏微分方程的逐步求解，而DPDM求解过程主要为模型扩散路径的单次求解，而代理模型计算和映射耗时仅为数秒。因此当模型规模增大，求解时间步增多，PDE计算耗时将快速增长，而DPDM方法的优势则更为显著，这也是该方法适用于真实混凝土结构长期服役模拟的关键原因。

表2 PDE与DPDM求解效率及存储空间需求对比Tab.2 Comparison of solution efficiency and storage requirement between PDE and DPDM

4 结论

（1）基于混凝土细观碳化过程的扩散路径依赖特性，通过单次扩散路径计算取代全模型迭代求解来适应混凝土结构截面分析需求，通过水泥浆体代理模型求解来考虑服役环境的精细变化，因此DPDM方法在大体积、长周期混凝土结构全断面细观碳化模拟中具备优势。

（2）通过对比材料碳化试验以及传统偏微分求解结果，DPDM方法能够较为精准地模拟混凝土材料水化效应和碳化过程，其中关键物质CH的浓度随空间和时间的变异性分布也能够精准模拟预测。在100 mm边长的正方形混凝土截面模型中，DPDM方法能够保持细观碳化模拟的精度，并将求解效率提升20倍的同时降低结果存储所需空间至传统方法的2%以内。

作者贡献声明：

李 越：试验设计，程序设计，数据分析，撰写初稿。

阮 欣：项目构思，课题指导，提出论文框架，论文修改与审定。

尹志逸：参与试验设计与试验结果分析。