COMSOL多物理场仿真：多孔介质两相流与药剂扩散建模全解析

上个月帮一个水务设计院的朋友做仿真，课题名字起得挺长：“COMSOL仿真建模：自由多孔介质两相流物质传递与水池自重作用下药剂扩散”。说白了就是一件事：一个水池如果发生渗漏，药剂在周围土壤里怎么跑、扩散多远、浓度能到多少。这类问题在水处理、危废防渗、化工储罐区环评里非常常见，也是COMSOL多物理场耦合建模里非常典型的一个组合。

这个项目表面看是“多孔介质两相流+稀物质传递+固体力学”三个模块叠在一起，实际上里面有不少坑：物理场接口怎么选、水平集参数怎么给、自重带来的孔隙率变化怎么传到渗透率上、瞬态求解怎么保证不发散。这篇文章我把完整建模思路和实操步骤整理出来，包含可以直接套用的参数表和表达式，适合正在做多孔介质迁移、渗漏扩散、流固耦合相关课题的工程师和研究生参考。做这类仿真之前，先把下面三个问题想清楚，后面能少走一半弯路。

1. 项目全景：这个仿真到底在算什么

1.1 标题拆解：三个关键点

先别急着开COMSOL，把标题拆清楚。“自由多孔介质两相流物质传递”这个说法听起来绕，实际拆开就是三件事。

第一，“自由多孔介质”。这里指的并不是某个叫“自由”的材料，而是说多孔介质区域边界是开放的，流体可以从边界自由进出，或者模型里同时存在自由流动区域和多孔介质区域。在COMSOL里处理这种问题经常用Brinkman方程把自由流动和多孔介质流动统一起来，界面处不需要额外拼装，一个方程域内自动过渡。新手最容易犯的错，就是把自由水体区域和多孔介质区域当成两个接口分别建模，再手动耦合，结果界面条件怎么都对不上。

第二，“两相流”。多数实际问题就是水-空气两相流，比如药剂溶液下渗，孔隙里的空气要被顶出来，这就是两相驱替过程。物理上涉及毛细力、黏性力、重力三者竞争，数值上要追踪相界面。COMSOL里追踪界面有水平集、相场、VOF三种方式，后面我会详细对比。

第三，“物质传递”。药剂在水相里溶解、随水流对流、在浓度梯度下扩散，还可能被土壤颗粒吸附、发生降解反应。严格说是水相浓度场的问题，但因为土壤孔隙曲折，扩散系数要做修正，不能直接把自由液体里的扩散系数拿来用。

“水池自重作用下药剂扩散”是后半句，很多人一看到自重就懵了，以为要做什么复杂的流固耦合。实际上工程上最常见、最稳妥的做法是分步解耦：先用固体力学接口算水池自重引起的地基压缩和沉降，得到应变场，再用经验公式把应变映射为孔隙率和渗透率的变化，最后把这个更新后的渗透率分布喂给两相流和传质模型。只有在科研级、需要精细捕捉变形过程的时候，才需要用到移动网格做全耦合。这个帖子我会把两种路径都说清楚。

1.2 应用场景与模型定位

这类模型的实际工程背景非常广。最常见的是储罐、污水处理池、尾矿坝的渗漏评估：池体如果发生裂缝或防渗层失效，受污染的水往土壤里渗，污染羽怎么扩展，下游地下水会不会被污染。另一个典型场景是农田灌溉、盐碱地改良：灌溉水里的盐分或药剂在非饱和土壤中的运移问题，本质上就是两相流加传质。还有土壤修复工程，比如原位化学氧化，把氧化剂注射到地下，看它怎么在含水层里扩散覆盖。

这类问题的共性特点是：介质是非均质多孔介质，流动是两相甚至三相的，污染物有扩散、吸附、反应等多重过程，而且往往伴随着上部结构荷载引起的介质变形。能把这些因素统一到一个模型里评估的，COMSOL这种多物理场平台优势非常明显。做环评报告或方案设计时，这种仿真结果比简单的手算经验值强太多。

1.3 两步走的整体建模思路

我的模型组织方式分两步走，前后衔接，逻辑非常清晰。第一步建一个“基准模型”，不管水池自重，先模拟自由多孔介质里的水-气两相渗流和药剂对流扩散，验证物理过程正确、结果合理。第二步加“自重影响模块”，用固体力学算自重造成的地基压缩，把应变转成孔隙率和渗透率修正量，重新跑两相流和传质，对比两种情况下的扩散差异。

这样做有几个好处。第一，问题隔离，哪个环节出问题一眼就能看出来。如果基准模型都不收敛，就不用谈后面的自重耦合。第二，计算量可控。全耦合瞬态求解经常需要几百步甚至上千步，稍微有点网格质量问题就发散。第三，工程上可解释性强。分步计算依然能回答“自重会让扩散范围扩大多少”这个核心问题。

有人可能会问，为什么不用移动网格做一个流固耦合大模型？说实话我一开始也想这么干，但做完整个模型之后发现必要性不大。水池自重引起的土体应变通常是百分之几的量级，对孔隙率和渗透率的影响可以用经验公式精确表达，而且移动网格会让水平集方程的求解难度陡增，经常出现网格质量退化。所以除非水池本身是可变形结构、变形量很大，否则单向参数映射是性价比最高的方案。

2. 物理场接口选择与控制方程要点

2.1 自由多孔介质两相流：水平集接口与Brinkman方程

COMSOL里跟多孔介质两相流相关的接口主要就是两个方向：一个基于Darcy定律的两相流，一个基于Brinkman方程的两相流。两者的差别在于，Darcy定律忽略了惯性项，适用于流速很低、纯粹渗流的场景；Brinkman方程保留了黏性项，能同时描述自由流动区域和多孔介质区域的速度过渡，适合“自由多孔介质”这种边界开放或水流速度不能完全忽略的情况。

我选择的是带水平集的多孔介质两相流接口，原因在于水平集方法在界面捕捉方面相当稳定，对二维和三维模型都适用，计算资源消耗适中。相场方法物理上更精细，考虑接触角、界面张力更严格，但方程阶数高，更容易出现不收敛。VOF方法在大尺度自由表面流里非常成熟，但在多孔介质这种孔隙尺度微小、毛管力占主导的场景下并不占优。我把三者整理成了一张对比表，方便你按自己课题做决策。

水平集接口里有个界面厚度参数，直接决定界面锐度。这个参数很关键，取太小了网格跟不上，取太大了界面糊成一片。我实测下来，界面厚度取最小网格尺寸的1～2倍是最稳妥的，同时要在界面经过的区域做局部网格细化，否则算出来的驱替前缘会不真实。

2.2 稀物质传递：对流-扩散方程与有效扩散系数

药剂在水相中的迁移用稀物质传递接口来描述，控制方程是标准的对流-扩散方程，里面最关键的一项就是有效扩散系数。很多人直接把药剂在水中的分子扩散系数填进去，这其实是错的。多孔介质里的固体骨架挡着，扩散路径是弯弯曲曲的，有效扩散系数肯定小于自由液体里的扩散系数。

有效扩散系数的常用做法是引入弯曲因子和孔隙率修正，形式是有效扩散系数等于分子扩散系数乘以孔隙率再乘以弯曲因子。弯曲因子的经验公式有很多，比如Bruggeman型、Millington-Quirk型，它们的本质都是把孔隙率、饱和度折算成一个小于1的修正因子。实际工程中，如果手头有土柱实验数据，最好用实验值反推这个因子，比任何经验公式都靠谱。

这个传质方程里还要注意对流项的速度来源。在多孔介质两相流里，溶质只在水相里跑，所以对流速度应该是水相的孔隙流速，而不是两相的平均流速。有些初学者直接把达西速度放进去，遇到空气饱和度高的区域，对流项就会偏大，算出来的扩散前锋明显不对。COMSOL里可以通过饱和度变量把水相速度筛出来，这个细节后面实操环节再细说。

2.3 水池自重：固体力学与参数映射

水池自重这部分，物理上属于固体力学范畴，在COMSOL里用固体力学接口，材料按线弹性模型处理。水池底板和地基土壤都要建模，底部给固定约束，四周给法向约束或自由边界，整个结构承受重力体积力。如果要考虑水池装满水的工况，侧壁还要加静水压力载荷，这个很容易漏掉。

真正核心的部分是如何把自重带来的变形影响传给多孔介质。工程上我用的是“固体体积不变”近似：变形导致的体积应变会使得孔隙体积变化，从而改变孔隙率。设初始孔隙率为，体积应变为，则更新后的孔隙率可以表达为一减去固体体积分数再除以体积变化系数。压缩的时候体积应变为负数，孔隙率减小，符合直觉。

渗透率的变化用Kozeny-Carman经验公式处理。这个公式把渗透率和孔隙率关联起来，形式是修正渗透率等于初始渗透率乘以初始孔隙率相关项再乘以更新后孔隙率相关项。它的好处是只依赖孔隙率就能算渗透率变化，岩石、土壤、混凝土等多孔材料都能凑合用。要注意，这个公式是经验性的，绝对数值不一定准，但能表达“压缩导致渗透率下降”的趋势，工程评估够用。

2.4 多物理场耦合的两种路径：单向更新与移动网格

前面提到过，自重和多孔介质流动的耦合有两条路。第一条是单向参数更新路径：先算固体力学，提取应变场，映射为孔隙率和渗透率修正，再带着修正后的参数跑两相流和传质。第二条是全耦合路径：用移动网格或变形几何让流体域跟随固体变形，每一步更新网格和物理场。两条路适用场景不同，我做了个对比。

单向更新路径实现简单、求解稳定、计算速度快，适合变形量不大的工程评估。它的局限性是假设变形不影响几何形状，只影响材料参数，这在大变形、大位移问题里不成立。全耦合路径物理上更严格，能捕捉到变形对流动通道的直接改变，但计算代价高、收敛难度大，经常需要调时间步和阻尼因子。课题如果只是要回答“自重会不会让扩散距离增加”这种工程问题，单向更新路径已经足够了，没必要为难自己。

这里还有个小的技术细节，如果你要用变形几何做全耦合，网格变形过大会导致单元扭曲，计算直接崩掉。一个缓解办法是把网格控制得足够细，同时用超弹性平滑来分散变形，但这也只能撑住中等变形量。所以我的建议很明确：先做单向更新，结果存下来，如果评审或导师非要看动网格，再单独做一个简化版全耦合模型做对比。

3. 实操搭建：从几何到求解器的完整步骤

3.1 几何建模：水池与多孔介质域的快速搭建

几何模型我在二维平面里建，三维版本原理一样，就是计算量指数级上升。整个模型分成三个域：水池底板区域、水池上方水体区域（如果做满水工况）、周围和底部的地基多孔介质区域。用矩形和布尔运算就能搭出来，不需要导入任何复杂CAD几何。

这里我建议用工作平面来画图，在2D工作平面上把二维截面画好，再设置好各域的材料归属。工作平面的好处是后期调尺寸方便，改一个宽度参数，几何自动更新，不需要重新画。很多网友问“comsol工作平面的作用是什么”，其实就是相当于CAD里的草图平面，所有几何都在平面上定义，再通过拉伸、旋转变成三维，习惯之后改模型效率高得多。

水池和地基的界面要特别注意，如果两个域的材料属性不同，界面要单独选出来设置接触条件。多孔介质和固体水池之间，如果是“池底裂缝渗漏”的工况，还需要在界面处切出一个小口来代表裂缝，用边界条件控制泄漏速率。如果只是“池底材料本身渗漏”，那直接把水池底面指定为多孔介质透水边界就行，不需要建裂缝几何。

3.2 材料参数和变量定义

材料参数是整个模型的地基，我这里列一组可以从零开始跑的典型参数，你根据自己的实际情况替换数值。需要注意的是单位一定要统一，COMSOL虽然能识别部分单位，但建议所有输入都确认一遍。

参数建好后，自定义变量里要写几个关键表达式。孔隙率更新和渗透率更新这两个表达式是核心，我写成COMSOL变量形式：

text复制eps_p_updated = 1 - (1 - eps_p0) / (1 + solid.evol)
kappa_updated = kappa0 * (1 - eps_p0)^2 / eps_p0^3 * eps_p_updated^3 / (1 - eps_p_updated)^2

第一行是孔隙率修正，利用固体力学接口的体积应变变量来更新孔隙率；第二行是Kozeny-Carman渗透率修正。注意不同版本COMSOL的体积应变量名可能有差异，比如有的版本叫solid.evol，有的版本需要自己写表达式，稳妥的做法是在“固体力学>应变”节点里查看实际变量名。只要变量名对应上，表达式是通用的。

3.3 边界条件与初始条件设置

两相流模型的边界条件分几类。模型底部和两侧设成开放边界，允许流体自由进出，对应“自由多孔介质”的设定。顶部如果是地表，设大气压边界，压力为0。水池底部作为泄漏入口，水相饱和度设为1，压力设为池内水头产生的压力值。如果建模满水工况，水池侧壁还要加静水压力。

水平集接口的初始条件很关键。初始时多孔介质区域充满空气，水平集函数处处设为0；水池底部附近的介质区域初始设为水，水平集函数为1。这样重力作用下水开始往多孔介质里渗，空气被挤压排出，形成明显的两相驱替前缘。

传质模型的初始浓度设为0，水池内部区域初始浓度为初始浓度。泄漏入口处浓度固定为初始浓度，这就是药剂的连续扩散源。模型四周的传质边界可以是敞开边界或零通量边界，实际要按场地是不是封闭来选。需要注意传质方程里的速度项，应该选水相速度，而不是总速度，这是前面强调过的。

3.4 网格划分的细节与技巧

网格划分对两相流模型成败影响非常大。我的经验是要做两层网格处理。第一层是对整个多孔介质区域用自由三角形网格，尺寸设置为能够分辨饱和度梯度的尺度，通常用“较细化”预设就能满足。第二层是在可能的相界面前沿区域做局部加密，加密区的网格尺寸要小于水平集界面厚度的一半。

结构力学部分的水池可以用映射网格或四边形网格，这种规则区域用结构化网格能显著减少单元数量，提高计算速度。算自重变形时不需要特别细化，网格尺寸稍大也没关系，因为变形场是平滑的，网格相关性不强。如果后续要做移动网格，那整个区域的网格质量都要严格控制，最小雅可比不能出现负值。

我经常看到有人一上来就用“极细化”网格，然后跑了一天还没算完。其实多孔介质流动的解通常比较光滑，网格密度有一个收益递减点。先用默认网格跑通流程，确认物理结果合理，再逐步加密做网格无关性验证，这才是高效的做法。

3.5 求解器配置与瞬态时间步控制

求解设置这里有两个研究。第一个研究做稳态固体力学，计算水池自重下的变形和应变场。不需要过多设置，固定约束和重力载荷加好后，直接求解即可。第二个研究做瞬态两相流和传质，时间范围从0到3600秒，输出步长按需设置，比如每10秒输出一个结果。

时间步进的选择是最大难点。两相流水平集方程对时间步非常敏感，步长太大马上发散，步长太小计算量爆炸。我的做法是用BDF方法，初始步长设置很小，比如1e-3秒，然后让求解器自适应增大。同时启用了全耦合求解器，并设置阻尼因子为0.9，实测下来比分离式求解器稳定很多，虽然每步计算量稍大，但总计算时间反而短。

如果你还是遇到不收敛，那就启用“分离式求解器”，把两相流和传质分开解，先算流动，后用流动结果算传质，每个时间步内迭代几次。这样牺牲了一些精度，但稳定性大幅提升，工程结果是可用的。分离式求解器适合模型调试阶段，确认没大问题后再改回全耦合。

4. 常见问题与排查技巧实录

4.1 不收敛：从时间步到阻尼因子

模型不收敛是两相流仿真里最常见的坑，严重的时候几乎是必发一次。我总结的排查顺序是：先看网格有没有负雅可比，再看时间步是否太大，然后看水平集界面厚度和网格尺寸是否匹配，最后调整求解器阻尼因子。大多数情况都能在第四步之前解决。

有一回算一个带水池自重的全耦合模型，总是到第37步就发散。排查了很久发现是网格在变形区域被拉得太长，雅可比接近0。把该区域的网格加密，改用超弹性平滑，问题就解决了。所以一旦发散，先看网格质量，而不是盲目调小时间步，否则只会浪费时间。

如果是水平集参数引起的发散，症状通常是界面附近压力场出现锯齿状振荡。这时把界面厚度参数调大一点，或者强化局部网格加密，振荡基本就会消失。水平集界面厚度不要追求越小越好，要根据实际网格能力来定。

4.2 水平集界面模糊或振荡

界面模糊的典型表现是饱和度的过渡带特别宽，驱替前缘看不清楚。原因通常是网格太粗，分辨率不够。界面厚度参数设置的比局部网格尺寸大太多也会导致这个问题。经验是把界面厚度设为最小网格尺寸的1倍左右，同时保证加密区的网格尺寸均匀，不要突然从细网格跳到粗网格。

界面振荡则表现为饱和度在界面附近忽高忽低，出现过冲。这往往跟时间步长太大有关，导致界面位置在一个时间步内跨过多个网格。把时间步上限调小，或者启用自适应时间步进并设置更紧的误差容差，界面基本会稳定下来。工程上，界面有轻微数值振荡是可以接受的，只要不导致非物理的负饱和度就行。

4.3 solidworks导入警告与几何修复

很多人在做水池这类规则结构时会用SolidWorks建模再导入COMSOL，结果打开模型时弹出一堆警告。最常见的是“退化实体”、“间隙”、“无效面”之类的提示。这个问题本质上是中间格式转换时的数据丢失，SolidWorks的实体模型转成STEP格式时，面片、倒角、小孔都有可能出现拓扑错误。

处理方法是在COMSOL导入面板勾选“修复几何”，然后到“几何”节点里执行“删除实体”、“忽略小特征”、“形成联合体”这些操作。如果模型不大，更省事的办法是直接在COMSOL里重建几何，水池结构本身很简单，用参数化矩形拉伸就能建出来，不会引入任何外部格式问题。

还有一些人导入后发现“绘图为空”，这个不一定是几何丢了，很可能是软件版本不兼容导致导入不完全。先查看“几何”节点的日志信息，有没有几何实体被移除的提示，再决定是清空重新导入还是重建。

4.4 绘图为空与数据导出的正确姿势

“comsol提示绘图为空”这个问题在多个版本里都有人遇到。很多时候不是数据没算出来，而是绘图表达式的名称写错了，或者数据集没有选对。比如你在“二维绘图组”里选了“切面”，但模型是二维的，切面数据集自然为空。我的排查习惯是：先在“数据集”节点确认有可用的解数据，再检查“值”选项卡里表达式名称是否与变量名完全一致，注意大小写。

数据导出这个环节也容易出问题。结果后处理时想导出一条剖面上的浓度分布，经常有人把“计算”和“导出”混为一谈。应该用“导出”节点，先选择数据集、设定要导出的点集或线集，再选择表达式，输出CSV或TXT文件。如果导出数据点特别多，建议先插值到规则网格再导出，否则下游处理非常吃力。

另外，如果你要用Python或MATLAB做进一步分析，COMSOL支持Livelink功能，可以在外部脚本里直接调用模型和结果数据，比手动导出CSV再读入高效得多。我通常是用Livelink for MATLAB做参数扫描和批量后处理，省掉大量重复劳动。

4.5 移动网格变形过大导致失败

如果做了全耦合移动网格，变形过大会导致网格单元翻转，求解直接失败。典型症状是警告“网格的雅可比为负”或者“无法更新网格”。缓解手段有这么几种：第一，把变形量分摊到多个时间步，避免一步变形太大；第二，用映射网格的自适应平滑，让变形尽量均匀；第三，把变形几何和物理场求解分离，每个时间步先更新几何再重新计算，不要全耦合同时解。

但说实话，在工程评估场景我越来越不推荐做移动网格全耦合，原因前面也说了，单向参数更新的结果差异在1%以内，计算成本节省一个数量级。只有当你的课题本身就是研究变形对传输过程的影响机制，全耦合才值得做。

5. 模型结果的工程解读与扩展方向

5.1 怎么解读饱和度与浓度场

模型跑完以后，很多人盯着云图发愁，不知道怎么看。我的建议是先看饱和度场，再看浓度场，最后把两者叠起来分析。饱和度场能告诉你水的扩散范围和前缘位置，可以判断药剂溶液是否已经渗到某个边界。浓度场能告诉你药剂扩散的浓度分布和衰减趋势，重点看等浓度线的扩展距离。

工程评估时，单独看一个时间点的云图是没有意义的，要把若干时间点的结果导出来，画成“扩散前锋距离-时间”曲线，或者“下游监测点浓度-时间”曲线，这才是能给甲方或评审看的核心成果。我在项目里最常用的输出是特定深度处浓度随时间变化曲线，以及最终稳态或长时间后的浓度空间分布。

如果你发现浓度场扩散范围明显大于饱和度场范围，说明药剂主要靠分子扩散在迁移，对流作用弱。如果浓度场前缘和饱和度前缘基本重合，说明对流占主导。这个判断对工程方案选择非常重要，两种机制对应的防渗策略完全不一样。

5.2 自重影响到底有多大

对比有自重和无自重的两组结果，我发现自重通常会让渗透率下降，理论上应该导致扩散变慢、范围变小。这正是单向参数更新的意义所在：如果你忽略自重，结果偏保守，偏安全。在工程环评中偏保守有时是故意的，但如果偏保守太多，可能导致防渗方案成本不必要的提高。

实际算例中，自重引起的孔隙率变化如果是10%左右，渗透率可能下降20%到30%，扩散距离的差别可能在10%上下。这个幅度足以影响工程决策了，所以不能简单忽略。更重要的一点是，自重的影响并不是均匀的，靠近水池底板正下方的区域压缩最严重，远离水池的区域影响很小，所以渗透率在空间上变成非均匀分布，这对流动和扩散路径有显著影响。

如果要对自重做更精细的评估，可以考虑把水池拆成空池、满水、带水压三种工况分别计算。空池只有结构自重，满水工况还有水压力作用在地基上，地基实际受到的竖向应力更大，渗透率下降也更多。这些工况之间的差异，往往比模型参数的细微调整更影响最终结论。

5.3 后续扩展：吸附、反应、裂缝和AI辅助建模

这个模型后续可以扩展的方向非常多。最直接的是加吸附项，在多孔介质稀物质传递接口里增加Langmuir或Freundlich等温吸附模型，让固体骨架吸附一部分药剂。实际操作不难，就是在传质方程里增加一项存储项和反应项，但要注意吸附会让有效扩散延迟，浓度曲线会变缓，不能再用简单对流扩散模型来拟合。

加化学反应也很实用，比如药剂降解按一阶不可逆反应处理，只需加一个衰减项。这里要注意反应速率常数和药剂半衰期的换算关系，不要随便填一个数量级离谱的值。可以在材料参数里定义半衰期，然后用表达式计算反应速率，这样可读性更强。

还有一个很多人感兴趣的方向是裂缝建模。多孔介质里有一条裂缝，流体会优先沿裂缝流动，扩散前锋会变得极不规则。在COMSOL里可以用裂隙流接口或定义一个高渗透率薄层域来近似，裂缝的渗透率通常比基岩高出几个数量级，对结果的影响非常大。这个扩展很适合后续做一个专题模型。

最后提一个新趋势，现在COMSOL新版带了MCP服务，可以让我们通过外部大模型工具来辅助建模和脚本生成。比如用自然语言描述想改哪个参数，工具能帮你直接修改模型树、生成参数扫描脚本，省去大量手动操作。我试过用这个功能自动生成批量后处理脚本，确实节省了不少时间。不过它只是辅助工具，物理概念和建模思想还得自己把住关，不能完全交给工具代劳。

实操收尾：说点个人心得

整套模型做完，我最大的体会是“先分解，再耦合”这个思路比硬啃耦合方程要高效得多。物理上看起来复杂的多孔介质两相流、传质、固体变形交互问题，一旦拆成几个标准模块，每个模块都有成熟的物理场接口和求解策略，剩下的工作就是参数映射和数据流转。真正难的不是COMSOL操作，而是你对自己要算的物理过程有没有清醒的认识。

另一点经验是，耦合类仿真要想在公司或课题组里真正落地，必须把模型参数写得足够清晰、可替换。我习惯把所有关键参数都定义成全局变量，文件名、几何尺寸、材料参数一目了然，后续改一个数重新算一遍，几分钟就能出对比结果。这比每次都在界面里翻半天去找某个设置要效率高得多。

最后再分享一个小技巧：模型收敛性实在解决不了的时候，先算一个简化版，比如退化为单相流、把多孔介质当作均匀材料、忽略自重影响，把流程跑通。把简化版的结果和完整版结果对比，往往能精准定位是什么模块导致的问题。这个方法帮我解决过好几次“莫名其妙就是发散”的困境，建议你也试一试。