这几年在北方沿海、盐渍土地区跑项目，混凝土耐久性病害里最让人头疼的组合就是冻融循环和氯离子侵蚀同时出现。单独做冻融试验或单测氯离子扩散，结果往往偏乐观，现场却早早出现顺筋裂缝和表层剥落。原因在于这两个过程根本就不是简单叠加，而是挤在同一个孔系统里互相放大。要搞清里面的逻辑，最合适的工具就是数值模拟——把温度、水分、盐分、裂缝放在一起算，去看冰晶生长和氯离子迁移之间到底怎么“相爱相杀”。
这篇文章不是把理论方程抄一遍，而是从机理、建模思路、关键参数到可落地的实现过程，完整讲一遍混凝土冻融和氯盐耦合劣化的数值模拟方法。适合正在做混凝土耐久性仿真的研究生、想引入多物理场计算的结构工程师，以及看过一些文献但不知道如何把两类劣化模型真正接起来的同行。
1. 冻融和氯盐为什么总在同一个工地捣乱
1.1 冻融先动手：孔隙冰晶的“物理打击”
混凝土本质上是一种多孔材料，水化产物堆叠之后留下大量毛细孔和凝胶孔。表面看起来密实，内部其实像一块细密的硬质海绵。潮湿环境下，这些孔里吸饱了水；温度降到冰点以下，孔中的水开始结冰。水变成冰体积膨胀约9%，如果孔隙饱和度足够高，冰晶没有地方“伸展”，就只能向孔壁施压。
这里有个常被忽略的临界点：当孔隙水饱和度低于约91.7%时，结冰产生的体积膨胀恰好能被孔内剩余空间容纳，理论上不会产生破坏性压力。一旦超过这个值，混凝土内部就会真正感受到来自冰晶的挤压。这是所有冻融破坏计算的出发点——如果整体饱和度不够高，模型的力学响应就应该是温和的；只有局部或全局饱和度突破阈值，损伤才开始累积。
冻融对混凝土的打击并不是一次完成的。每次循环产生的压力都会让微裂缝扩展一点，反复冻融相当于无数个小型“内胀循环”在孔壁附近不断做疲劳损伤。回到工程现场，你能看到的现象就是构件表面一层一层地起皮、剥落。这种破坏不是从构件中心开始，而是从表层向内推进，因为表层最先吸水、最先受冷、也最先结冰。
1.2 氯离子趁虚而入：裂缝就是“高速公路”
氯离子在完好混凝土中的迁移速度其实非常慢。普通硅酸盐混凝土在20℃下的氯离子扩散系数在10的负12次方量级，也就是每秒钟只能移动几个纳米到几十个纳米。如果没有裂缝，氯离子要靠纯扩散穿透5cm保护层需要很多年，这也是混凝土本身具备耐久性的根本原因。
但冻融循环改变了这个格局。冻融产生的微裂缝在混凝土内部形成了一张连通网络，氯离子从“在致密浆体里慢慢蠕动”变成“在裂缝网络里快速流动”。更有实测数据表明，裂缝区域的氯离子表观扩散系数可以比完好区高出一个数量级甚至更多。这个影响没法用安全系数简单覆盖，因为裂缝是逐步发展的，扩散能力的提升也是动态的。
1.3 双人组合是“乘数效应”而不是简单相加
氯盐环境反过来又让冻融破坏变得更剧烈。盐溶液冰点比纯水低，理论上不太容易结冰，但氯盐的吸湿性很强，会让混凝土表面长期保持高湿度，饱和度居高不下。更麻烦的是，盐在孔隙中析出结晶时也会产生结晶压力，和冰晶压力形成“双重挤压”。
我曾经用一张关系图来理解这个耦合链条：
冻融循环产生微裂缝，微裂缝加速氯离子渗入，氯盐提高混凝土吸水能力和饱和度，高饱和度又让冻融更容易产生破坏。这个闭合循环意味着一旦启动，劣化速度就是指数级上升的。单因素模型很难复现这种“1+1远大于2”的现象，只有把物理场交叉耦合起来才能抓住本质。
2. 数值模拟的整体思路：先选尺度，再分物理场
2.1 建模尺度是第一个分岔口
做冻融和氯盐耦合模拟，第一件事不是打开软件，而是先想清楚你要在哪个尺度上回答什么问题。研究冰晶在某个纳米级孔隙里怎么形核、长大，要用分子动力学或者相场法；研究一根混凝土柱在20年服役期内保护层怎么失效，则要用连续介质尺度的多场耦合模型。
很多初学者把这两个尺度混在一起，试图用一个模型同时模拟“孔隙中冰晶的形貌演化”和“氯离子沿构件深度的浓度分布”。这种模型不是不能搭，但计算代价极高，而且参数通常多到无法标定。工程和寿命预测中更有价值的路线，是把冰晶生长机理“打包”成宏观本构关系——它不需要我们知道每一粒冰晶长什么形状，只需要知道在一定温度、一定孔隙率、一定含水量下会产生多大的结晶压。
从数值方法角度看，这和做地震波数值模拟、地下水流场模拟用的离散化思路是相通的。有限元、有限差分、有限体积这些方法处理温度场、浓度场、应力场都有成熟方案。我们真正要做的不是从头发明算法，而是把物理过程之间的耦合关系定义清楚。
2.2 把“黑箱”拆成五个相互交换数据的场
一个完整的劣化模型通常需要至少五个场协同工作：温度场、湿度场（反映孔隙水饱和度）、冰含量场、氯离子浓度场和力学损伤场。它们各自遵循不同的控制方程，但在混凝土内部不断交换信息。
温度场驱动冰晶生成？不准确，更严谨的说法是：温度场决定了孔溶液能否结冰以及结冰量；冰晶生成产生的压力进入力学场，带来应力与损伤；损伤变量反过来改变氯离子的输运速度；氯离子浓度又影响孔溶液的冰点和吸湿性能。谁都不能单独算完然后去推下一步，必须同步迭代。
下面的表格展示了各场之间的主要交互关系：
| 物理场 | 控制变量 | 主要驱动力 | 对其它场的影响 |
|---|---|---|---|
| 温度场 | 温度T | 环境温度、水化热 | 决定冰晶量、影响扩散系数 |
| 水分/湿度场 | 饱和度Sw | 干湿循环、盐的吸湿 | 饱和度影响冻融损伤是否发生 |
| 冰含量场 | 成冰率φ_ice | 温度、孔溶液冰点 | 结晶压、体积膨胀、损伤源 |
| 氯离子浓度场 | 浓度C | 扩散、对流、吸附 | 降低冰点、改变孔结构 |
| 力学损伤场 | 损伤变量d | 结晶压、膨胀应力 | 放大氯离子扩散系数 |
2.3 上来就做“全耦合”往往是灾难的开始
见过不少团队想把所有因素一把梭哈进模型，最后代码写了一万多行，参数数据库却根本填不齐。我的建议是分三步走：
第一步，先做单向耦合。比如只算“冻融损伤→氯离子扩散系数增大”，不考虑反作用。这样程序能跑通，结果也能和文献对照。第二步，加入“氯离子→冰点降低”的反馒头反馈，比较单向和双向的差距有多大。第三步，才考虑要不要加入吸水、盐结晶、钢筋锈蚀等附加过程。
强耦合数值模型的每个额外反馈环都可能带来收敛困难和标定难题。每加一个环，都要问自己一句：它能引起多大的变化？如果反馈可以忽略或者可以用系数近似，就不值得为它增加一个完整的偏微分方程。这个取舍能力，比会写方程更体现经验。
3. 冰晶生长的模拟要点：冰冻不是0°C一刀切
3.1 孔径越小，冰点越低
在绝大多数数值模型中，设置冰点边界的默认做法是把0℃当成冻结门槛。这在实际冻融模拟里会带来明显的偏差。混凝土孔隙中的水不是纯净水，而且受到孔壁表面张力和弯曲界面的约束，结冰温度会明显下降。描述这个现象常用Gibbs-Thomson关系做一个量级估算：
ΔT = 2·γ·T_m / (ΔH_f·ρ_i·r)
让我把这个公式拆开看。γ是冰水界面能，约0.04 J/m²；T_m是纯水冰点温度273K；ΔH_f是冰的熔化潜热约3.34e5 J/kg；ρ_i是冰的密度约917kg/m³；r是孔隙半径。代入后整理出一个很好用的近似关系：
ΔT ≈ 71.4 / r(nm) K
孔径100nm的毛细孔冰点下降约0.7K，换算下来大概在-0.7℃附近结冰；孔径10nm的凝胶孔冰点下降约7K，要到-7℃才会冻结。这解释了为什么温度刚降到-3℃到-5℃时，并不是所有孔隙水都会结冰——先冻的是大孔，小孔里的水还在亚稳态液态。
做混凝土冻融数值模拟时，正确做法是结合孔径分布曲线，给定一个温度值，反推出当前温度下能结冰的最大孔径，再统计该孔径以下或以上的孔体积占比，得到成冰量。这一步看起来只是在修正材料参数，实际效果却非常关键。忽略它会导致模拟中在-2℃时所有孔隙水瞬间结冰，损伤被严重高估。
3.2 氯离子对冰点的“反向作用”不能忽略
孔溶液中存在氯离子时，冰点还会进一步降低。这里用理想溶液的冰点下降公式能做初步估计：
ΔT_f = K_f · m · i
K_f是水的摩尔冰点下降常数，约1.86 K·kg/mol；m是溶质的质量摩尔浓度；i是范特霍夫因子，对NaCl约为2。举个例子，混凝土表面接触3%的NaCl溶液，相当于质量摩尔浓度约0.52mol/kg，理论冰点降低接近1.9℃，冰点约为-1.9℃。如果氯离子浓度上升到5%附近，冰点会降到-3℃以下。
这个效应在耦合模拟中有双重作用。一方面，高浓度氯盐会让孔溶液更难结冰，纯水冻融循环和盐水冻融循环的临界温度是截然不同的；另一方面，盐溶液的高吸湿性让混凝土表面的平衡饱和度提高，在宏观尺度上又是促进冻融破坏的因素。很多模型只做了冰点降低却忽略了吸湿增水，或者只考虑吸水却不改冰点，都是不完整的。
3.3 把冰晶生长“折算”成力学载荷
从力学角度看，模拟孔隙中每一粒冰晶的外形几乎不可行，工程上普遍采用等效体积应变法。核心思路是：假定孔溶液结冰后，固态冰的体积增量可以映射为一个各向同性的体积应变增量。
在有限元软体中实现时，最省事的方法是把它当成“温度应变”处理。先根据温度和孔径分布算出当前步的成冰率，再转换成等效热膨胀系数增量。软件只看到结构因为“某种场”产生了体积膨胀，从而计算出应力和开裂。这个方法虽然没有精细刻画冰晶的局部压力分布，但在保护层尺度的寿命预测中已经被验证是足够准确的。
损伤计算建议采用基于应变等效的损伤模型。每完成一次冻融，材料的弹性模量按损伤变量折减，残余强度和拉应力承载力同步下降。损伤变量从0到1变化，0代表完好，1代表完全丧失承载力。在冻融荷载下，损伤通常会率先发生在距表面10mm到30mm的区域内，然后逐层向内推进。
4. 氯离子迁移模拟：扩散系数只是一个开始
4.1 传输机制不止扩散一种
氯离子在混凝土里的迁移方式，严格来说有五种：浓度梯度引起的扩散、毛细吸附造成的对流、压力梯度驱动的渗透、电场作用下的电迁移、以及碳化或化学结合带来的固化。只写一个纯Fick第二定律，解析解确实漂亮，但和真实环境差得很远。
完整的氯离子输运方程至少应该包含扩散项、对流项和吸附源项：
∂C/∂t = ∂/∂x[D(C,T,d)·∂C/∂x] - v·∂C/∂x + S_bind
其中v是水分流动的达西速度，S_bind是结合氯离子的源项。在干湿交替区域，毛细吸附引起的对流往往是表层氯离子快速积累的重要途径。冻融循环让表层微裂缝增多，毛细吸附作用也会被放大。
4.2 扩散系数必须随损伤动态放大
这是整个耦合模型中最容易出彩、也最容易出错的地方。氯离子扩散系数D在模型中被当成常数，是绝大多数失败模拟的共同根源。实际条件下，D值受三个因素显著影响：
- 温度：温度升高，离子热运动增强，扩散加快，可用Arrhenius型公式修正；
- 水灰比和养护龄期：水化越充分，孔结构越致密，D越低；
- 损伤程度：冻融微裂缝和高应力损伤会成倍提升D。
在损伤场和浓度场耦合时，可以用指数形式描述损伤带来的扩散加速：
D(d) = D_0 · exp(a·d)
D0是完好混凝土的基准扩散系数，d是当前损伤变量，a是标定系数，通常在3到8之间。当损伤为0.3时，指数因子大约在2.5到10倍之间；当损伤发展到0.6以上，扩散能力可以放大到10倍量级。这正是冻融后氯离子快速渗透的数值体现。
4.3 自由氯离子和结合氯离子的区别
模型里算出的浓度若是总氯离子含量，会高估对钢筋的腐蚀威胁，因为实际引发锈蚀的是孔隙溶液中的自由氯离子。部分氯离子会与水泥水化产物发生物理吸附或化学结合，生成Friedel盐，这部分不再自由迁移。
计算时需要考虑结合等温线。常用Langmuir或Freundlich吸附模型描述结合氯离子量C_b与自由氯离子浓度C_f之间的关系：
C_b = α·C_f / (1 + β·C_f)
这种非线性关系会导致一个重要现象：在同等表面氯离子浓度下，表观扩散前沿比线性模型移动得更慢。若同时模拟冻融裂缝使自由氯离子通道增多，吸附的影响会相对减弱。
判断结果是否合理时，要留意浓度剖面的形状。若某处出现浓度大于表面浓度或者负浓度的“异常峰”，通常是从扩散方程源项设置出了问题，而不是真实物理现象。
5. 双向耦合的算法实现：让温度、损伤和浓度场真正“对话”
5.1 双向耦合的关键链条
两个过程一旦要在同一个模型里互相影响，至少要明确以下几条耦合链：
第一链：温度降低→成冰率上升→结晶压增加→应力损伤增加。第二链：损伤增加→扩散系数放大→氯离子浓度提升。第三链：氯离子浓度提升→孔溶液冰点下降→相同温度下成冰率降低。第四链：氯离子浓度提升→表面吸湿增强→饱和度上升→冻融更易发生。
这四条链同时存在时，冻融和氯离子侵蚀就形成了正负反馈交织的复杂系统。数值模拟的价值正是在于把这些反馈放在统一时间轴上迭代，看清哪一个过程在某个时间段占据主导地位。
5.2 一个实用的时间步算法框架
下面给出我在算例中常用的伪代码思路，实际工程中用Python+FEniCS或者商业软件二次开发都可以套用这个结构：
for step in range(total_steps): # 1. 求解温度场，更新每点温度值 T = solve_thermal(T_prev, bc_temp, dt) # 2. 结合氯离子浓度修正冰点，计算成冰率 T_freeze = freeze_point(T, C_free, pore_size) phi_ice = ice_fraction(T, T_freeze, pore_distribution) # 3. 用力学平衡求解应力，并由结晶压更新损伤场 sigma = solve_mechanics(phi_ice, d_prev) d = update_damage(d_prev, sigma) # 4. 更新损伤相关扩散系数 D = D0 * exp(a * d) # 5. 求解氯离子输运方程，包含吸附项和对流项 C_free = solve_transport(C_free_prev, D, sw, dt) # 6. 将自由氯离子浓度反馈给冰点计算，进入下一循环 T_prev = T d_prev = d C_free_prev = C_free如果使用COMSOL Multiphysics，可以在“系数型偏微分方程”接口中分别建立温度场和浓度场，再用“固体力学”接口计算应力场。损伤变量放在“域常微分方程”接口里，作为跨物理场共享变量。ABAQUS用户则适合写USDFLD子程序，在每个积分点更新损伤与扩散系数。选择哪个软件不重要，关键是时间步长内的耦合顺序要固定、清晰、非循环依赖。
5.3 数值稳定的“三件套”
跨物理场耦合模型最容易出现的问题就是不收敛或结果振荡。结合我自己的调试经验，下面三点只要做到，基本能避免80%的报错。
第一，时间步长必须取到足够小。氯离子扩散用显式格式时有稳定条件dt小于等于dx平方除以两倍D。网格取2mm、扩散系数取10的-11次方时，稳定时间步长约为200秒量级，太大就会发散。
第二，反馈源项要做平滑化处理，尤其是冰点切换和损伤突增。直接从0到1的阶跃变化极易造成迭代震荡，可以用双曲正切函数在冰点附近做一个宽度约0.5℃的过渡带。
第三，耦合变量建议采用上一时间步的旧值，也就是松耦合，不要在一个时间步内做太多次强迭代。对于寿命尺度模拟，物理上本来就没有必要让每条反馈链在每个增量步内完全同步收敛，松耦合效率更高也更稳。
6. 一个可复现的算例：100mm厚混凝土板的冻融+氯盐耦合模拟
6.1 工况与参数设置
这个算例模拟一维混凝土板在单面接触盐溶液并经历冻融循环时的劣化过程。模型宽度10mm或20mm均可，深度100mm，只分析厚度方向的传热与扩散。这样处理方便复核，也能快速对比不同参数的影响。
表里的参数都是常见取值，你可以在此基础上替换成自己的试验数据：
| 参数 | 取值 | 说明 |
|---|---|---|
| 板厚 | 100 mm | 模拟保护层与内部区段 |
| 水灰比 | 0.45 | 中等强度混凝土 |
| 初始孔隙率 | 12% | 假设均匀分布 |
| 表面氯离子浓度 | 0.8%（混凝土质量） | 相当于长期接触盐溶液 |
| 基准扩散系数D0 | 5e-12 m²/s | 20℃、未损伤状态 |
| 温度范围 | -20℃~20℃ | 单次冻融周期24小时 |
| 损伤放大系数a | 5 | 由冻融试验数据拟合 |
| 氯离子吸附系数 | Langmuir单参数简化 | 控制自由氯离子含量 |
6.2 建模分步清单
第一步，建立几何模型，100mm深度画成一根线段或矩形条，网格尺寸建议2mm。温度场在厚度方向是瞬态传热问题，如果只做一维模拟，横向网格可以很粗。
第二步，设置首尾温度边界和初始温度。温度边界按三角波变化，表面在24小时内从20℃降到-20℃再升回20℃，内部边界按绝热处理。对流换热系数在试验箱中很难标定，可直接采用“实测表面温度曲线”作为第一类边界条件，比估算热交换系数更可靠。
第三步，给材料添加孔径分布假设。如果没有压汞法实测数据，可用简化双峰分布代替：一部分孔径集中在20nm到200nm，另一部分集中在1μm到10μm。温度下降到某个值后成冰率先在大孔部分发展，随着温度降低逐步扩展到小孔。
第四步，在力学接口中把冰晶体积应变转化为“等效温度应变”。这里不是真的给结构降温，而是根据成冰率换算一个虚设的膨胀应变，再让有限元求出应力场和损伤场。
第五步，计算氯离子浓度场。注意把自由氯离子浓度作为状态变量，损伤变量通过扩散系数放大项实时耦合进去。表面浓度作为定值，内部末端设置为零通量边界。
模拟结束后分别输出第50、100、200个循环时刻的损伤分布和自由氯离子浓度剖面。观察损伤峰值的位置和移动规律。
6.3 结果怎么读
完好混凝土中的氯离子浓度剖面通常是单调递减的平滑曲线。引入冻融以后的典型变化是：近表面处浓度升高幅度比纯扩散模型更快，浓度曲线斜率明显改变，在损伤严重区域附近可能形成一个浓度平台甚至“局部上凸”。这代表裂缝网络已经实质性地改变了离子输运路径。
损伤方面，损伤峰值会从表面向内部逐层推进，而不是在所有位置均匀增长。同时要注意表面层剥离后的边界条件更新问题——如果模拟中混凝土表层已经剥离，理论上氯离子的“有效表面”也在内移，这时候需要考虑几何或表面边界更新。初版模型可以先不考虑剥蚀，用于比较相对趋势是完全足够的。
我一直建议先跑这个一维模型，再扩展到二维真实截面。一维模型计算量小，半小时内能出结果，能快速验证双向耦合代码是否正常。二维构件加钢筋、加箍筋以及非均匀水灰比之后，用同样的耦合框架逐步扩展即可。
7. 常见问题与排查经验速查表
7.1 从一次次不收敛中摸出来的对照表
| 现象 | 可能原因 | 处理办法 |
|---|---|---|
| 氯离子浓度出现负值 | 时间步过大或扩散系数突变 | 缩小时间步；把损伤放大因子改为指数平滑形式 |
| 冻融循环中损伤无增长 | 冰点未考虑孔径效应，所有孔一起结冰后压力被平均化 | 改用孔径分布计算成冰量，区分大孔小孔冻结顺序 |
| 损伤变量超过1 | 迭代过程中未限制上下限 | 在每个增量步末尾对损伤做截断处理，并强制不低于前一步 |
| 温度场计算发散 | 表面边界条件加载太快 | 初始几步采用斜坡加载，或在温度峰值附近做平滑 |
| 氯离子冰点影响不明显 | 冰点修正公式没有耦合到成冰率计算 | 检查求解顺序，确保当前步的浓度在计算冰点前已经更新 |
| 结果对网格太敏感 | 损伤局部化导致网格依赖性 | 引入特征长度正则化，或采用非局部损伤模型 |
| 模拟剖面比试验实测偏慢 | 未考虑干湿循环吸附对流 | 补充饱和度场和毛细吸附速度项 |
7.2 我的几条高频经验
如果只能给一条建议，那就是别急着堆代码。先用手算和文献数据校核单一物理场，比如不加耦合时模拟出来的氯离子剖面能不能和纯扩散解析解对上。解析解校不住，后面加了耦合只会更糊涂。
冻融和氯盐耦合模型本质是热-水-力-化学多场耦合问题，这类模型没有唯一正确的建模方式。如果目标是回答“构件寿命还剩多少年”，就没必要在孔隙尺度纠结冰晶形状；如果目标是写一篇机理研究论文，连续介质模型里的简化冰点公式又可能不够精致。先把目标讲清楚，再决定简化哪些物理过程。
还有一点容易被忽略：单位统一。这个模型牵涉温度、浓度、压力、损伤多个量纲，MPa和Pa混用、百分比和质量浓度混用都是现实发生过的问题。每次运行大算例之前，我习惯先挑一个单元做单步试算，手算确认新增量没有量级错误再铺开跑。这个习惯帮我省下的调试时间，至少以周计算。
