搞混凝土耐久性的人大概都绕不开“碳化”这个话题。钢筋锈蚀的起点往往就是碳化前沿推进到了钢筋表面把原本稳定的钝化环境毁掉了。而COMSOL这类多物理场仿真软件恰好很适合把碳化这个涉及气体扩散、水分迁移、化学反应、甚至还带点应力损伤的复杂过程放进同一个模型里算。我自己断断续续试了两三个月踩了不少坑也把一些关键参数和求解思路捋顺了这篇就把整个建模探索过程梳理出来供打算入手的朋友参考。1. 混凝土碳化模型为什么非要多场耦合不可1.1 碳化过程的物理本质与经典化学路径碳化听起来是“二氧化碳把混凝土吃掉”但本质上是一个气-液-固三相反应大气中的CO₂先通过混凝土孔隙网络向内部扩散溶解到孔溶液里形成碳酸再跟水化产物氢氧化钙CH以及C-S-H凝胶等发生反应生成碳酸钙CaCO₃。用化学式表达就是CO₂ Ca(OH)₂ → CaCO₃ H₂O这个反应本身不复杂复杂的地方在于它发生在多孔介质里而且每一步的速率都受不同因素控制。CO₂进得来主要靠的是孔隙通道的开放程度和湿度状态——太干了CO₂扩散很顺畅但没有水反应没法进行太湿了孔隙被水堵得严严实实气体进不去碳化同样停滞。这就是为什么混凝土的碳化深度和表面湿度之间往往存在一个峰值关系而不是单纯的线性对应。搞清楚这一层就能理解为什么传统经验公式比如经典的开方时间律 x k√t虽然工程上够用但换一个环境条件就经常失准。因为它把所有影响碳化的因素都压进了一个综合系数k里一旦环境湿度、温度、混凝土配合比变了这个k就不再可靠。1.2 单场模型的局限性与耦合关系的建立只算CO₂扩散的单一场模型说白了就是解一个扩散方程把扩散系数取一个经验值然后看浓度分布。这种模型在实验室恒温恒湿条件下拟合得还不错但拿到室外或者实际工程里偏差就大了。核心原因在于各个物理过程之间是互相牵制的。湿度场直接决定孔隙水的饱和度而水饱和度最少通过三个路径反向影响碳化过程。第一影响CO₂的有效扩散系数水膜越厚气体扩散通道越窄第二影响碳化反应的速率水是反应物也是碳酸电离和钙离子溶出的介质第三影响Ca(OH)₂的溶解和CaCO₃的析出形态湿度过低时碳化产物会堵塞孔隙反而改变后续的扩散路径。更麻烦的是温度场也掺和了进来。温度一变气体扩散系数跟着变反应速率常数按Arrhenius公式指数变化。如果混凝土还承受荷载、有裂缝甚至微裂缝损伤那么扩散系数还要叠加一个应力/损伤因子。裂缝就像给CO₂开了一条高速公路让碳化前沿局部快速深推。这就引出了“多场耦合”的本质诉求与其分场地猜参数不如把CO₂浓度场、湿度场、温度场以及可能的应力场放在同一个框架里按时间步推进让每个场在每一时刻都把最新状态反馈给其他场。COMSOL恰恰是干这个的——它是少数普通工程师不用自己从头写有限元代码、就能比较轻松搭建多场耦合模型的平台。2. COMSOL中碳化模型的场结构与方程映射2.1 核心因变量与对应物理场接口的选择在COMSOL里搭碳化模型先要决定算哪几个场。我用的方案是“21”结构核心算CO₂浓度场c和气相/水相饱和度相关的湿度场θ温度场视问题需要决定是否打开——如果你做的是恒温实验对照完全可以关掉温度场如果模拟室外全年变化的外部环境建议打开因为温度对反应速率常数的影响确实不能忽略。CO₂浓度场我推荐用**“多孔介质中的稀物质传递”**接口Transport of Diluted Species in Porous Media。这个接口自带多孔介质有效扩散系数的定义框架支持把扩散系数做成孔隙率、饱和度的函数而且可以通过源项把碳化反应消耗CO₂的部分直接写进去。需要注意的细节是这个接口的浓度因变量单位一般默认是mol/m³而混凝土碳化文献里常用的单位是kg/m³或者体积分数做单位换算的时候务必认真我最初就在这上面吃了一次亏——所有的方程和边界条件都写了唯独单位没验计算结果碳化深度差了十几倍。湿度场有两种选择如果只关心水饱和度的宏观变化可以用**“水分在多孔介质中的传递”**接口求解一个Richards方程因变量取压力水头或者饱和度如果考虑温度伴随的湿分传递也有更复杂的“Moisture Transport”多物理场接口。实际建模时我建议采用以相对湿度RH为因变量的简化形式好处是碳化速率经验公式里很多都用RH做自变量可以直接引用不需要额外换算。2.2 化学反应源项与速率方程的表达式实现碳化反应的CO₂消耗项是模型里最容易飘、也最容易不收敛的部分。工程文献中常见的做法是把碳化反应看作对CO₂的一阶消耗速率表达式写成R k · [CO₂]这里的k不是常数而是跟CH含量、温度、湿度有关的综合反应速率。更严谨一点可以用Papadakis等人在上世纪九十年代提出的多步反应模型形式把CH的消耗也纳入控制方程写成两个载荷方程∂m_CH / ∂t -k_CH · m_CH · c ∂c / ∂t ∇·(D_eff · ∇c) - k_c · m_CH · c在COMSOL里第一个方程可以作为一个额外因变量添加到“系数型偏微分方程”接口里也可以直接用一个域的常微分方程ODE来表达——因为CH在每一点只随时间变化、不参与空间迁移。这样处理的好处非常明显你不需要引入整个化学模块计算量也小得多。反应速率常数k的温度依赖用Arrhenius公式写k k₀ · exp(-Ea / R / T)Ea取值参考石灰石烧制和水泥水化类似反应的经典文献碳化反应的活化能大致在7~10 kJ/mol这个数量级比水化反应低得多。温度影响没有想象中剧烈但能明显改变长期碳化深度的对比——温度每升高10℃反应速率大概提升差不多一倍左右。2.3 扩散系数与湿度损伤的耦合表达式多场耦合建模的核心体现就在有效扩散系数D_eff的取值上。纯粹的气相扩散系数和有效扩散系数差距很大后者要乘上孔隙率ε、曲折度因子τ再折减掉水膜阻塞的影响。一个常用的形式是D_eff D_CO2 · ε^α · (1 - S) ^ βS是水饱和度α和β都是配合比相关的经验参数。没有实测数据时α取1.5~2.0、β取2.0~3.0是文献里常见的区间算出来的扩散系数一般不会离谱。有实验配合比压实度数据的话最好根据压汞法孔径分布结果标定。另外碳化消耗了CH、生成了CaCO₃本身会改变孔隙结构——碳化前期生成的碳酸钙会细化孔径扩散系数下降后期一旦开裂损伤扩散系数又会反弹。要在COMSOL里实现这种动态变化可以把孔隙率本身写成一个随CH消耗量变化的变量再把它传入有效扩散系数表达式里。我实际试过用以下方式处理ε_eff ε₀ - Δε · (1 - m_CH / m_CH_initial)然后把ε_eff替换原来的ε₀放进D_eff表达式。这种做法比完全固定孔隙率的模型更接近实测但代价是数值稳定性变差需要配合更细的时间步长和更保守的求解器容差。3. 几何构建、网格策略与瞬态求解的一线实测记录3.1 从二维截面开始还是直接上三维混凝土碳化模型在几何上不建议一上来就建三维。绝大多数情况碳化问题沿试块暴露面法向占主导侧面边界效应只在边角处明显。我自己的建议是先建一维模型——没错就是一维几何一个从表面到内部的线段——把方程和耦合关系全部跑通再扩展为二维轴对称或者二维截面。一维和二维的COMSOL模型架构几乎一样改动只在几何和边界上。一维模型的好处是计算极快几十年的时间跨度几分钟就算完特别适合用来做参数扫描。等到所有参数都校准得差不多了再换成二维去看边角效应或者非均匀暴露面遇到网格问题也不至于手忙脚乱。3.2 网格划分碳化前沿的“追锋”策略混凝土碳化的浓度场有个特点碳化前沿CO₂浓度从高到低剧烈变化的区域非常薄可能只有几毫米而整个构件厚度可能是几十厘米甚至几米。如果全模型用均匀网格为了分辨这薄薄一层网格总量会爆炸。COMSOL里应对这种情况首选边界层网格——在暴露表面处生成多层紧密网格然后向内部逐渐过渡到较粗网格。实际测试下来前沿区域网格尺寸在0.5~1 mm、内部过渡到5~10 mm二维模型几十万自由度就能算得很流畅。如果遇到底面附近出现振荡我通常会先检查是不是网格不够密而不是急着调求解器。还有一个小技巧如果后期需要追踪碳化前沿位置可以在后处理里把“碳化深度”设成一个等值面变量——定义一个表达式比如Ca(OH)₂含量降到初始值的某个比例工程上常用初始CH含量的10%的位置即视为碳化前沿。配合一维变量提取 “Cut Line” 功能每一时间步都能直接出一张碳化深度-时间曲线。3.3 瞬态求解器配置与时间步控制碳化是一个跨度极长的时间过程室外构件按年算加速碳化实验也至少按周算。COMSOL默认的瞬态求解器用的是BDF向后差分公式自适应步长理论上能处理刚性问题但如果反应源项写得太“硬”k值偏大、时间常数极小BDF会不断把步长压缩到很小导致计算时间暴增。我实测下来比较好用的配置是求解器选择“使用自适应时间步进的BDF最大阶数2~3”不要用默认的5阶高阶BDF对非线性刚性问题稳定性反而差。容差从默认的0.01放宽到0.005或0.001的级别太紧会拖慢速度太松则碳化深度曲线会明显毛糙。输出时间步用“范围对数”或“范围步长”来控制我的习惯是列出几个特定时间节点比如1天、7天、28天、90天、180天、365天、2年、5年、10年让求解器在这几个节点之间自己加密插值。时间步长如果按小时算一年就是8760个小时但碳化速率其实是随时间的开方递减的——前期快、后期慢。用等间距时间步效率很低改成对数分布后同样的计算量能多覆盖好几倍的时间跨度。3.4 边界条件的物理标定与表面状态模拟边界条件设置直接决定模型算出来的是“实验场景”还是“实际场景”。室内加速碳化实验的边界很好办——暴露面固定CO₂浓度比如体积分数3%温度湿度恒定其他面设为通量对称或零通量。户外实际结构就不一样了表面CO₂浓度随季节和昼夜波动而且雨淋会显著改变表层湿度。如果要模拟干湿循环COMSOL的“事件”接口可以派上用场设定一个周期函数让表面湿度边界在某个区间内交替变化同时用“下喷淋”的富水表面条件去模拟降雨。我试过用简单的“方波”函数来近似雨季和旱季的切换效果比恒定湿度边界明显更接近真实结构的碳化分布——湿润期表层湿度升高CO₂扩散几乎停滞但内部湿度被“封住”到了旱季碳化又以更快的速度推进。这种往复切换会造成碳化深度呈锯齿状前进肉眼看着奇怪实际上这正是真实构件在干湿交替下的特征。如果还考虑应力场需要额外引入“固体力学”接口把碳化引起的体积膨胀或收缩作为初始应变加上去。不过坦白说应力-碳化完全耦合的模型收敛难度非常大我目前也只是做到了局部开裂试件的单向耦合应力改变扩散系数但扩散不影响应力场完整双向耦合的进展还在推进。4. 常见失败场景与参数标定的实战排查心得4.1 实验数据标定diffusion coefficient的三条路径模型建得再漂亮没有可靠参数就是空中楼阁。标定CO₂有效扩散系数的办法按优先级排序有三种第一种最可靠用实验室加速碳化深度数据反演。办法很朴素——把扩散系数D_eff设成全局参数用COMSOL的参数扫描功能跑不同D值取实测碳化深度与模拟碳化深度误差最小的那个。这个方法看似粗暴但极其有效。参数扫描在COMSOL 6.x里已经是内置功能配合批处理几乎不用额外写脚本。第二种用快速碳化后试件不同深度的碳化产物总量分布来标定这需要做热重分析或XRD定量。好处是不仅能标定扩散系数还能标定反应速率常数缺点是实验成本高。第三种纯经验估算适合没有实验数据的预研阶段。经典经验公式推荐的普通C30混凝土在自然环境下的表观扩散系数大约在1e-8~1e-7 m²/s量级加速碳化3% CO₂浓度时通常是自然环境的10~20倍左右。先用这个量级试算再根据现象修正总比靠猜强。4.2 模型不收敛与负浓度的七种对策用过COMSOL做瞬态非线性问题的人应该都见过那个红色报错框——解在某个时间步发散然后整个求解器直接挂掉。碳化模型最常见的失败模式我整理了一下基本逃不出以下七种一是反应源项过强导致时间步崩溃。对策是去掉源项先跑纯扩散模型确认扩散模块没问题再由小到大逐步加载源项。二是初始条件设置不合理比如初始湿度场设为100%饱和度但边界却暴露在40%湿度的空气中差值过大第一步迭代就崩了。这种时候用“平滑过渡”的初始表达式给湿度场一个自表及里的梯度渐变收敛会舒服很多。三是单位错误前面提过不再赘述。四是网格在反应前沿过疏浓度场和湿度场梯度无法被正确捕捉。五是求解器容差太紧我见过有人把相对容差改到1e-8来追求精度结果算了一个多小时还在第一步徘徊——真没必要。六是材料参数表达式里用了不连续函数比如if语句拐点处求导不连续极易引发振荡建议改用连续的平滑近似函数smoothstep或analytic函数。七是边界条件突变——干湿切换的方波函数如果跳变幅度太大建议加一段斜坡过渡时间。4.3 碳化深度提取与实验对照的数据匹配技巧模型跑出来了怎么跟实验对照也是有讲究的。实验室测量碳化深度用的是酚酞试剂喷雾后测变色边界而模型里“变色边界”对应的是孔隙溶液中pH骤降的位置。严格来说pH降到8.5~9.0左右的位置才是酚酞变色线而不是CH含量归零的位置。所以在后处理里提取碳化深度时我建议不要直接用CH含量0作为判据而是定义一个新的变量基于CH消耗程度和孔隙饱和度估算一个等效pH值取pH9作为等值面位置。这样算出碳化深度会和实验室酚酞实测值对得很齐。这个小技巧我自己摸索了好久才总结出来——如果直接拿CH含量阈值去提往往比实验值偏深。温度、湿度数据的时间序列也建议和气象站或者其他实测数据挂钩。室内加速实验还好环境稳定户外数据就一定要用真实的气象数据否则模型输出得再平滑也只是一条“好看的假曲线”。5. 模型扩展方向的个人实践体会碳化模型跑通之后可扩展的方向其实非常多。我目前尝试过且比较推荐的有三个方向一是跟氯离子侵蚀耦合做“碳化-氯盐”双重耐久性分析——碳化改变了孔隙液的pH直接影响氯离子结合能力两者在真实沿海环境下几乎是同时发生的二是跟钢筋电极过程耦合把碳化前沿推进到钢筋表面之后激活钢筋腐蚀的电化学反应描述这样就能直接输出锈蚀起始年限这是工程上最关心的话三是跟细观尺度的网格模型结合把宏观扩散计算的结果传递到细观骨料模型里从而研究界面过渡区对碳化路径的偏折作用。第三个方向计算量很大我还在优化阶段。但前两个方向在COMSOL里已经有比较成熟的模块接口尤其是腐蚀模块和变形接口的配合能让模型真正服务于寿命预测而不只是发一篇论文。最后说一个最实际的建议如果是从零开始不要去翻软件内置的案例库找混凝土碳化——里面确实没有现成的案例。建议从“多孔介质稀物质传递”的官方示例改起先把纯扩散跑通再逐步添加反应项和湿度耦合。每一步改动都建一个新分支保存否则改崩了连退回哪里都找不到。建好一个可复用的模板之后以后再换配合比、换环境条件基本上就是改参数的事了。
