混凝土碳化这件事做过耐久性的同行都知道它不像强度那样拉个试块压一下就能出结果而是要等几年甚至几十年才看得出端倪。所以我拿到一批加速碳化试验数据时第一反应不是急着拟合平方根公式而是想把背后的CO₂扩散和反应过程用数值模型完整跑一遍。选COMSOL的原因很简单它不需要自己从头写有限元程序又能通过修改PDE方程把碳化机理的细节塞进去。项目标题里提到的“PDE建模”其实就是把碳化过程翻译成一串带反应源的扩散方程让软件去解。这个过程我断断续续折腾了两周踩了不少坑也攒了不少经验。这篇文章把从方程设计到COMSOL实操的完整路线梳理出来希望能帮到正在做同类仿真的朋友。1. 碳化问题是什么为什么非要上COMSOL建模1.1 碳化过程的物理化学机制混凝土碳化听着复杂剥开看核心就三件事CO₂从表面往内部扩散溶解在孔隙水里形成碳酸然后和孔隙溶液中的氢氧化钙以及水化硅酸钙反应生成碳酸钙。CO₂扩散这一步由Fick定律描述。混凝土内部并不是密实的水化后残留的毛细孔、凝胶孔是CO₂进入的通道。扩散系数受孔隙率和含水饱和度影响极大孔隙水会堵住气路湿度越高扩散越难反过来湿度太低反应所需的水分不够反应又会受限。所以工程上常说“相对湿度50%到70%是碳化最凶的区间”既保证了CO₂扩散又有足够水分让反应发生。反应这一步生成的碳酸钙沉积在孔隙壁上会让孔隙率微微下降但反馈到扩散系数上的影响通常只有百分之几到十几除非是表层完全碳化形成致密壳。绝大多数简化模型会忽略孔隙率变化专注在CO₂浓度空间分布和反应前沿推进上。对钢筋来说最危险的还不是碳化本身而是碳化引起的pH下降。混凝土孔隙液正常pH在12.5以上钢筋表面有钝化膜碳化后pH能跌到9以下钝化膜失效只要水和氧气到位锈蚀就开始了。所以工程上关心的物理量其实非常集中碳化前沿走到哪里了什么时候穿过保护层到达钢筋表面。1.2 实验实测的真实成本让人想换数值手段实验室里最常见的做法是加速碳化试验把试件放进密闭箱通入浓度20%左右的高浓度CO₂恒定温度20℃、相对湿度70%28天或者更长时间后劈开试件喷酚酞指示剂测碳化深度。加速箱把自然界0.04%的CO₂浓度提高了几百倍用“时间压缩”的方式让28天等效于自然几年的碳化量。听起来效率很高但问题在于加速环境下扩散和反应的相对权重和自然暴露并不完全一致。高浓度CO₂下反应前锋对扩散的依赖变小碳化深度和时间呈更接近线性而非平方根的关系。直接拿加速箱数据外推50年寿命是有很大风险的。经验公式比如x k√t工程上很好用但k值完全依赖试验标定和配合比库换一种胶凝材料、换一个水灰比就要重新做一组试验。我想做的其实是机理层面的事情把碳化过程写成偏微分方程组用COMSOL求解器去推进时间这样能同时回答“碳化深度多少”和“过渡区有多宽”这类问题而且改变配合比、湿度、CO₂浓度只需要调参数不用重新做试验预判。1.3 为什么不用COMSOL现成的稀物质传递模块COMSOL自带的“稀物质传递”模块确实能加入反应源项边界条件、单位转换也都现成但真按碳化需求走下来会发现几个别扭的地方。第一碳化反应消耗的是CO₂同时也在消耗固体水化产物。如果只用稀物质传递描述CO₂浓度场固体反应程度的更新需要额外定义ODE标准模块里没有直接对应的变量。第二反应速率不是简单的线性消耗而是和“剩余可碳化物质含量”挂钩——你需要在PDE系统里同时求解CO₂浓度和固相反应程度两个场。第三后续如果想扩展温度场、湿度场耦合稀物质传递接口的自由度反而不如自己写的PDE直观。所以我最后选了“数学模块”里的PDE接口把碳化过程当作一个2×2的对流-扩散-反应系统来求解。这也是标题里叫“PDE建模”的原因COMSOL只是一个求解平台真正的功夫在方程怎么设、参数怎么取。2. PDE模型设计方程怎么选、怎么改2.1 三类碳化模型的思路差异碳化模型圈子里其实有三种典型做法。第一种叫瞬时反应前沿模型假设反应速度快到CO₂一进入未碳化区就被完全消耗碳化区和未碳化区之间存在一条明确界面。这种情况下碳化深度有解析解x √(2 · De · C₀ · t / a)x是碳化深度De是CO₂有效扩散系数C₀是表面CO₂浓度a是单位体积混凝土可碳化物质含量t是时间。这个公式工程界用了很多年优点是简单缺点也很明显它给不出碳化过渡区的浓度分布也没法反映反应速率变化的影响。第二种叫有限速率反应模型不假设反应瞬间完成而是给CO₂消耗速率一个显式的反应项常见形式是∂C/∂t ∇ · (De ∇C) - k · C · S其中S代表单位体积内剩余可碳化物质的量k是反应速率常数。这个模型能描述从完全碳化区域到未碳化区域之间还有一个部分反应的过渡带更接近真实情况。第三种是更复杂的多场耦合模型加入湿度场、温度场、孔隙率演化甚至钢筋腐蚀电流。这种模型COMSOL完全撑得住但对参数要求高很多参数没有实测值只能靠文献估计做起来容易变成“参数拟合大赛”。我做复现时选择的是第二种模型的扩展版——把固体反应程度S也作为一个独立变量求解形成双物种PDE系统。理由很实际收敛难度可控物理图像清楚参数数量也刚好能通过试验标定。2.2 双物种PDE系统的搭建我用的方程组是∂C/∂t - ∇ · (De ∇C) k · C · S 0∂S/∂t (k / ν) · C · S 0第一个方程描述CO₂气相浓度C的变化扩散项De ∇C让CO₂向内部输运源项-kCS表示反应消耗。第二个方程描述剩余可碳化物质S的变化S被反应逐步消耗ν是化学计量系数表示反应消耗1 mol mol的CO₂对应消耗多少mol可碳化物质。简化情况下ν取1。这里有个容易踩坑的点两个方程中的反应项符号必须严格相反否则总质量不守恒COMSOL算着算着就能看到碳化区和未碳化区之间出现奇怪的浓度跳变。De不是常数。我用的表达式是De D₀ · φ · g(RH) · ηD₀是CO₂在自由大气中的扩散系数约1.6×10⁻⁵ m²/sφ是孔隙率典型混凝土取0.1到0.2g(RH)是湿度影响因子很多文献用(1 - RH/100)^m的形式m在2到3之间η是孔隙曲折度相关修正项一般量级在0.01到0.1。乘下来De大约在10⁻⁹到10⁻⁸ m²/s量级和快速碳化试验反演的扩散系数能对上。表面浓度C₀也值得细算。大气CO₂体积浓度约400 ppm即0.04%20℃下用理想气体换算C₀ (1.013×10⁵ Pa × 400×10⁻⁶) / (8.314 × 293.15) ≈ 0.0166 mol/m³。不对这个数偏小因为这是分压直接换算实际孔隙中气相和液相溶解态CO₂还会重新平衡很多模型直接用0.03%到0.1%体积分率下的等效浓度参与计算。碳化箱加速试验里表面浓度可以到8到10 mol/m³量级和大气条件的差异一眼就能看出来。2.3 参数从哪来怎么标定参数标定是做数值模拟最花时间的一步我的顺序是先固定化学计量参数再用试验数据反演扩散和反应参数。可碳化物质含量a单位是mol/m³。水泥用量400 kg/m³的普通混凝土水化程度按0.7估算单位体积中能参与碳化的氢氧化钙和水化硅酸钙折算成CO₂消耗量大约在200到400 mol/m³之间。这个值决定了最终碳化深度的大基数直接影响平方根公式里的分母。De的获取有两条路。有条件的实验室直接做N₂/CO₂扩散试验测透气系数再换算有效扩散系数没条件就用快速碳化试验反演跑一组标准28天碳化数据用最小二乘法在模型里调De和k直到计算碳化深度和实测吻合。实际反演中De对深度的影响很大k对过渡区宽度的影响更大两者可以分开标定。k的绝对值不好定。它受湿度、温度影响强烈我参考文献常取1×10⁻⁵到1×10⁻⁶ m³/(mol·s)之间做量级估计再通过“碳化深度曲线形状”来细化。传感器多的话还可以用试件内部CO₂浓度时间来标定这个效果更准。实操提示碳化PDE的边界参数和初始条件一定要先做量纲检查。CO₂浓度用mol/m³时k的单位就得是m³/(mol·s)这两个单位对不上求解器不会立刻报错但结果会彻底飞掉。3. 在COMSOL里把PDE搭起来接口选择与实操细节3.1 系数型PDE还是一般型PDECOMSOL的数学模块里和PDE建模相关的接口主要有“系数型偏微分方程”和“一般型偏微分方程”两种。系数型PDE的模板格式是e_a · ∂²u/∂t² d_a · ∂u/∂t ∇ · (-c∇u - αu γ) β · ∇u a·u f只要把d_a设为1c设为De右边f写成反应源项就能表示CO₂扩散方程。好处是扩散项、对流项都有现成系数位定义起来快坏处是当源项涉及另一个变量S时f处写表达式会稍微绕一点要在“因变量列表”里把第二个场一并建好。一般型PDE的模板是e_a · ∂u/∂t ∇ · Γ f形式更简洁Γ是通量表达式f是源项。碳化的第一个方程就可以直接写Γ -De·∇Cf -k·C·S。第二个方程没有扩散项Γ 0f -k·C·S/ν。我的建议是如果只是复现碳化深度这种简单问题系数型PDE够用如果以后还想加氯离子耦合、碳化-锈蚀联合模拟直接用一般型PDE自定义自由度大得多不用中途改接口。双物种系统在COMSOL里要添加两个因变量C和S。分别设置初始值C初始0S初始S₀即可碳化物质含量按200到400 mol/m³给。这两步必须在“初始值”节点里同时写清楚漏掉任何一个求解器很容易把S算成负值。3.2 几何、网格与边界条件的选择碳化问题最常用的几何就是一根一维线从混凝土表面指向内部钢筋位置。很多人觉得一维模型太简单但工程保护层厚度就是10到50 mm的事一维模型给出的碳化深度-时间曲线已经足够回答设计问题。二维、三维模型主要用在研究棱角处双向碳化、裂缝附近的集中侵蚀等场景对新手来说没必要一上来就建三维收敛慢还看不出额外信息。几何长度取保护层厚度再留出一段余量比如取0.1 m。边界条件设置如下暴露表面C C₀Dirichlet边界代表环境中CO₂浓度恒定内部远端零通量∇C · n 0这条边界相当于对称面或者无限厚墙体截断边界。初始条件方面碳化开始前混凝土内部没有CO₂所以C(t0) 0S(t0) S₀。这个设定符合试验操作试件从标准养护室拿出来放到碳化箱之前的那一瞬间内部确实没有外来CO₂。网格划分是复现碳化模型最容易翻车的地方。碳化过程中CO₂浓度在扩散前沿附近形成一个陡峭下降的锋面如果网格太粗前锋会被数值抹平碳化深度结果偏大。我做过的经验是在表面和预期碳化深度区域做加密单元尺寸压到1到2 mm扩散前沿区域甚至可以到0.5 mm。COMSOL里可以开“自适应网格细化”但碳化锋面是随时间移动的自适应网格功夫不到位反而容易增加计算量不如直接在表层铺一层细密边界层网格。更稳妥的做法是先跑一个粗网格模型再跑一个加密一倍的网格模型看碳化深度曲线是否基本重合。如果结果差异超过5%就要继续加密。这类“网格无关性检查”虽然听着麻烦但能省掉后面一大堆对结果可信度的质疑。3.3 求解器配置与稳定性技巧COMSOL的PDE求解器默认配置对扩散-反应问题基本能跑但碳化问题的时间跨度太大碳化深度从毫米级推进到厘米级可持续几个月到几十年。这种多时间尺度问题需要对求解器做几个调整。时间步设置是从30天到50年步长不能给固定值。正确做法是让COMSOL用“自适应时间步”自动控制初始步长设置小一点比如1小时让反应初期快速扩散过程先稳定下来几分钟内完成的快速变化在自动步长控制下能捕捉到而不会让整体计算量爆炸。相对容差我习惯设到1×10⁻⁴比默认1×10⁻³稍微严格一点特别是当你关心过渡区内部浓度分布的时候。求解器类型选“隐式向后差分格式”即BDF对于刚性扩散-反应系统更稳。如果遇到收敛失败先检查两个方程的源项符号是否反了这是最常见的问题再检查S是否会变成负值——当S减小到0附近时反应项kCS理论上为0但数值迭代时可能因为过冲导致S略微为负负的S会令反应项反号把浓度推向正方向增长造成振荡。解决办法是在源项里加一个限制条件把S替换为max(S, 0)或者在方程里加上一个小的“保护项”。我在一般型PDE里直接写S if(S 0, S, 0)的辅助变量实测有效。关键经验最开始不要直接跑50年。先跑30天看碳化深度曲线是否和试验数据吻合再跑1年、10年逐步拉长时间。这样做的好处是能在每个时间尺度上发现参数设置偏差而不是最后等50年算完才发现结果完全不对浪费时间还难以定位问题。4. 模型复现中踩过的坑与排查实录4.1 负浓度和源项振荡第一个遇到的典型问题是CO₂浓度出现负值。现象很常规碳化前沿附近少数节点浓度在某个时间步跌到-0.01 mol/m³量级然后振荡。排查后发现原因有三层第一是S已经消耗到接近于0反应项kCS理论应该等于0但数值插值产生的过冲让S变成-1×10⁻¹⁰这类微小负值乘以C和k后产生了反向源项第二是隐式时间步在反应前锋移动较快时局部截断误差积累把浓度削成了负值第三是源项kCS中C和S都是变量非线性很强牛顿迭代在过渡区收敛半径不足。针对这三个原因逐个解决源项中S用max(S,0)保护保证S为负时不反应时间步上限做约束不让自适应时间步在反应前锋快速移动时跨得太大Newton迭代容差放宽松一点或者增加阻尼因子避免反复振荡。4.2 长时间求解突然卡死在某一步碳化模型在运行到后期特别是碳化深度穿越保护层到达远端边界附近时突然出现“求解器在时间步10.2年处无法收敛”的报错这是第二个高频问题。原因在于碳化前沿接近远端零通量边界时浓度场的梯度方向受到边界影响局部刚度矩阵的条件数变差时间步自动控制失效。解决方法是把几何长度从0.05 m加大到0.2 m让远端边界尽量远离碳化前锋这个做法等于用延长计算域的方式消除边界效应。计算量增加不大但收敛稳定性明显好转。还有一种情况是S接近完全消耗后第二个方程退化成∂S/∂t 0或更小量级此时方程组的雅可比矩阵接近奇异。处理方法是在第二个方程右边加入一个非常小的人工项比如-1×10⁻¹²·S让矩阵不至于奇异性拉满。这个做法听起来粗糙但数值稳定性提升非常明显。4.3 碳化深度怎么提取以及和文献对不上的问题算完模型后怎么把碳化深度提取出来又是一个麻烦问题。实验上是用酚酞试剂看变色边界pH低于9的区域不显色定为碳化区。数值模型里并没有直接的pH变量需要自行建立判据。我用过两种提取方式。第一种是设定S的剩余比例阈值S降到初始值S₀的5%以下的位置定义为碳化前沿。这个判据和试验酚酞边界比较一致因为酚酞变色的本质是氢氧根浓度骤降而氢氧根浓度和可碳化固相消耗高度相关。第二种是直接观察CO₂浓度等值面C浓度降到表面浓度1%以下的位置也可以作为碳化前沿指标。两种方法算出来的深度差异通常在1到2 mm量级对工程判断来说足够精准。实际做对标时经常发现计算碳化深度比加速试验实测深或浅50%以上。这时候先去查De和S₀是否偏离实际。我做参数校验时发现文献里给的De“有效扩散系数”往往是基于30年自然碳化数据回归的等效值和加速碳化28天环境湿度下的De并不是同一个值。湿度对De的影响是指数级的70%相对湿度下De比干燥状态能小一个数量级。在模型里必须明确给De加入湿度因子表达式而不是用单一常数。4.4 常见问题速查表问题现象可能原因解决方案浓度出现负值源项符号错误或S过冲为负源项里用max(S,0)检查两个方程源项符号相反求解器长时间不收敛几何边界过近、雅可比矩阵奇异加大计算域长度远离碳化前锋或对S方程加极小人工项碳化深度与试验差异大De湿度因子未考虑、S₀估算不准加入湿度影响项S₀根据胶凝材料用量和水化程度重新估算过渡区宽度异常α值偏小或网格太粗加密前锋附近网格校验α量级到1×10⁻⁵~1×10⁻⁶碳化深度随时间接近线性而非√t表面浓度过高反应速率控制降低表面浓度值或增大反应速率常数k让模型回到扩散控制网格加密后结果明显变化网格无关性未收敛对表层1~2 mm范围加密再加密一倍做对比以上表格里的问题我复现时几乎一条条都碰到了很多情况下不是方程错了而是参数的取值情境和假设条件没对上。数值模拟的一大好处是能快速定位坏处是如果参数体系不自洽报错或者结果不合理时会非常搞心态。5. 复现之后的经验沉淀我的体会是碳化PDE模型的复现本质上是把碳化机理从经验公式的“黑箱”变成可解释的“白箱”。COMSOL帮我们省了有限元求解的底层代码但模型的灵魂还是那组方程和参数。磨刀不误砍柴工花在参数标定和边界条件检查上的时间比直接跑模型的耗时更值得投入。最后分享一个操作习惯我在模型里布置了3个“探针”点分别在深度10 mm、25 mm、40 mm处记录CO₂浓度随时间的变化曲线。碳化前锋到达每个探针时浓度曲线会出现一个明显的转变拐点。这个拐点时间如果和试验中预埋传感器的响应时间对不上马上就能知道扩散系数的量级估计是否合理。这个方法比只盯着碳化深度曲线要灵敏得多排查问题时帮了大忙。