混凝土碳化这件事，做过耐久性的同行都知道：它不像强度那样拉个试块压一下就能出结果，而是要等几年甚至几十年才看得出端倪。所以我拿到一批加速碳化试验数据时，第一反应不是急着拟合平方根公式，而是想把背后的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.2；g(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还是一般型PDE
COMSOL的数学模块里，和PDE建模相关的接口主要有“系数型偏微分方程”和“一般型偏微分方程”两种。
系数型PDE的模板格式是：
e_a · ∂²u/∂t² + d_a · ∂u/∂t + ∇ · (-c∇u - αu + γ) + β · ∇u + a·u = f
只要把d_a设为1，c设为De，右边f写成反应源项，就能表示CO₂扩散方程。好处是扩散项、对流项都有现成系数位，定义起来快；坏处是当源项涉及另一个变量S时，f处写表达式会稍微绕一点，要在“因变量列表”里把第二个场一并建好。
一般型PDE的模板是：
e_a · ∂u/∂t + ∇ · Γ = f
形式更简洁，Γ是通量表达式，f是源项。碳化的第一个方程就可以直接写Γ = -De·∇C，f = -k·C·S。第二个方程没有扩散项，Γ = 0，f = -k·C·S/ν。
我的建议是：如果只是复现碳化深度这种简单问题，系数型PDE够用；如果以后还想加氯离子耦合、碳化-锈蚀联合模拟，直接用一般型PDE，自定义自由度大得多，不用中途改接口。
双物种系统在COMSOL里要添加两个因变量C和S。分别设置初始值：C初始0，S初始S₀（即可碳化物质含量，按200到400 mol/m³给）。这两步必须在“初始值”节点里同时写清楚，漏掉任何一个，求解器很容易把S算成负值。
3.2 几何、网格与边界条件的选择
碳化问题最常用的几何就是一根一维线，从混凝土表面指向内部钢筋位置。很多人觉得一维模型太简单，但工程保护层厚度就是10到50 mm的事，一维模型给出的碳化深度-时间曲线已经足够回答设计问题。二维、三维模型主要用在研究棱角处双向碳化、裂缝附近的集中侵蚀等场景，对新手来说没必要一上来就建三维，收敛慢还看不出额外信息。
几何长度取保护层厚度再留出一段余量，比如取0.1 m。边界条件设置如下：
- 暴露表面：C = C₀（Dirichlet边界），代表环境中CO₂浓度恒定；
- 内部远端：零通量，∇C · n = 0；这条边界相当于对称面或者无限厚墙体截断边界。
初始条件方面，碳化开始前混凝土内部没有CO₂，所以C(t=0) = 0；S(t=0) = 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₂浓度随时间的变化曲线。碳化前锋到达每个探针时，浓度曲线会出现一个明显的转变拐点。这个拐点时间如果和试验中预埋传感器的响应时间对不上，马上就能知道扩散系数的量级估计是否合理。这个方法比只盯着碳化深度曲线要灵敏得多，排查问题时帮了大忙。
