基于COMSOL相场法的随机孔隙结构驱替模拟实践

做驱替模拟这件事，我前后踩过不少坑。最早用简单几何通道做，界面推进四平八稳，和真实的多孔介质完全是两回事；换上随机孔隙结构之后，水驱油的界面一下子“活”了——指进、绕流、死角里的残余油全出来了。这就是基于COMSOL相场法与随机孔隙结构的驱替模拟案例想讲清楚的核心：用相场法捕捉两相界面，用随机生成的孔隙结构逼近真实岩层和微流控芯片里的拓扑，然后把两者放在同一个求解框架里，看流体怎么穿过多孔骨架。这篇内容适合正在做孔隙尺度渗流模拟的学生、工程师，以及想把COMSOL相场接口用在实际案例里的人。我会把随机孔隙结构怎么生成、相场接口有哪些关键参数、完整水驱油算例怎么设置、结果如何判读，以及我踩过的坑都过一遍。

1. 项目思路：为什么要把随机孔隙结构与相场法放在一起

1.1 驱替模拟真正想回答的问题

驱替（displacement）在工程上最常见的就是“二相流体在多孔介质中相互排挤”。石油工业里水驱油是这一过程的经典场景，二氧化碳封存也涉及超临界相驱替地层水，燃料电池气体扩散层里的水管理本质也是气体驱水。研究这类问题，最关心的通常有三个点：第一个是驱替的稳定性，注入流体能否形成一个整齐的前沿，还是像手指一样窜出去；第二个是被驱替相的残留量，比如水驱过后还有多少油被圈闭在孔隙里；第三个是驱替效率，突破之后含水率如何快速上升，还能不能靠改变注入策略把残余油采出来。

传统实验手段能给出宏观的采收率、压差曲线，但看不到孔隙内部的局部界面运动。而如果你的目标就是搞清楚“喉道怎么截断油滴”“哪些死孔隙完全进不去水”，孔隙尺度的两相流模拟几乎是唯一选择。把孔隙结构设为规则排列的圆柱或平行通道，虽然好算，但没有非均质性就捕捉不到驱替的指进现象，模拟结果会偏乐观，参考价值有限。随机孔隙结构正是为了打破这种“完美排列”，它让每个喉道、每个孔径变化都成为影响界面移动的变量，模拟结果才更接近真实的多孔介质。

1.2 相场法的取舍：为什么不是水平集或VOF

在多孔介质这类复杂几何里模拟两相流动，界面追踪方法有好几种流派，各有短长。VOF成熟、守恒性好，但它靠几何重构确定界面，在三维孔隙拓扑和液滴分裂合并频繁的工况下，表面张力的计算容易出现伪流动，接触角处理也比较麻烦。水平集方法界面平滑，法向曲率计算方便，但它的致命弱点是质量不守恒，需要反复重新初始化距离函数，算到后期你会发现油相面积悄悄变少或者变多。

相场法走的是另一条路线。它不显式追踪界面，而是引入一个连续变量来描述相态，界面只是一个厚度很薄的梯度过渡层，界面运动完全由自由能驱动的Cahn-Hilliard方程控制。拓扑变化发生时不需要特殊逻辑，液滴断裂、合并、在喉道处被截断都是自然演化出来的。更重要的是它本身满足对流-扩散形式的守恒律，只要网格和迁移率参数控制得当，质量守恒比水平集靠谱得多。

代价也很明显：为了分辨厚度只有几个微米的界面层，网格必须局部细化；迁移率参数取不好还会引入人为的界面弥散。而且它方程数多、非线性强，计算规模比VOF要贵。但在两相不相容、表面张力主导、且存在大量拓扑变化的驱替问题里，相场法提供的是长期稳定、拓扑无约束的求解路径。这也是我在COMSOL里坚持用相场接口的原因。

1.3 随机孔隙结构的意义与课题定位

随机孔隙结构到底“随机”在哪里？并不是随便撒一堆圆就叫随机。它有明确的约束条件：孔隙率要达标，颗粒之间不能重叠，也不能出现造成布尔运算退化的相切点，颗粒尺寸分布要符合目标区间，而且整个区域不能被某个大颗粒切断流动通路。这样生成的结构，每一个都是真实多孔介质的抽象样本。它的价值在于：(1) 提供统计学意义上的非均质性；(2) 可以批量生成多组样本，看结果对孔隙结构的敏感性；(3) 复现实验芯片上的颗粒堆积结构。

从课题定位上，这个模拟算例属于“微观渗流力学 + 数值多相流”交叉方向。它既不追求真实岩心的全尺寸重现，也不停留在规则管束的理论验证，而是把随机结构、相场界面动力学和驱替效率三者定量地连起来。你拿它来做参数研究、孔隙结构评价，或者作为更复杂反应的流场基底，都是合适的。

2. 随机孔隙结构的生成：几何是模拟的第一步

2.1 随机障碍物生成的约束条件

随机孔的骨架我用了最简单的二维模型：矩形域内随机分布圆形固体颗粒，流体在颗粒间的间隙中流动。生成算法看起来是“丢骰子”，实际操作中约束条件比想象中多。核心约束有四条：

• 圆心不能落在边界附近，否则布尔差集会在边界上留下薄片状尖角，网格质量差，求解器容易报错；
• 任意两个圆的圆心距离必须大于半径之和再加上一个最小间隙，否则颗粒接触点会是奇异几何，后续布尔操作和网格剖分都容易退化；
• 固体面积分数必须收敛到目标值，孔隙率才能精确控制；
• 随机种子要可控，这样每跑一组模拟都能复现。

最小间隙这个参数不能设得太小。我一开始设置0.5 μm，颗粒间距和网格尺度同一量级，结果剖出来的网格在间隙处单元质量一塌糊涂，局部网格尺寸被迫缩到纳米级，计算代价飙升但物理没有变好。后来把最小间隙拉到1~2 μm，匹配1 μm量级的网格，整体效率提升非常明显。

2.2 坐标生成脚本的写法与孔隙率控制

下面给一段可以直接用的Python脚本思路。它按目标固体面积分数逐步生成圆，并在圆与圆之间保持间距：

import numpy as np Lx, Ly = 1e-4, 1e-4 # 域尺寸 100 μm × 100 μm r_min, r_max = 5e-6, 1.5e-5 # 颗粒半径范围 5~15 μm solid_target = 0.4 # 目标固体分数（孔隙率 = 1 - solid_target） gap = 2e-6 # 颗粒间最小间隙 centers, radii = [], [] solid_area = 0.0 max_iter = 5000 domain_area = Lx * Ly for _ in range(max_iter): if solid_area / domain_area >= solid_target: break r = np.random.uniform(r_min, r_max) # 为避免边界问题，圆心缩进一个半径加余量 margin = r + 0.5 * gap cx = np.random.uniform(margin, Lx - margin) cy = np.random.uniform(margin, Ly - margin) # 检查与已有圆是否重叠 overlap = False for (cx0, cy0), r0 in zip(centers, radii): d = np.hypot(cx - cx0, cy - cy0) if d < r + r0 + gap: overlap = True break if overlap: continue centers.append((cx, cy)) radii.append(r) solid_area += np.pi * r * r porosity = 1.0 - solid_area / domain_area print(f"生成的颗粒数: {len(centers)}") print(f"孔隙率: {porosity:.3f}")

孔隙率是逐步逼近的，不是精确等于目标。因为每个圆的大小在生成时才知道，最后总会略超出或低于目标值。解决方法是当最后一步超出量很小时接受它，或者只按随机尝试重新来过。多生成几组种子你会发现：同样的目标孔隙率，颗粒数、孔径分布、连通路径完全不同，这本身就是随机结构的价值。

2.3 几何导入与布尔差集处理

坐标文件生成之后，要把这些圆送进COMSOL几何模块。最省事的方式是把圆心和半径输出为txt文件，再用脚本或外部CAD接口批量生成圆。如果你只做单组算例，直接在COMSOL图形界面里逐个添加圆也行，但颗粒数一多就不现实。我的建议是输出一个DXF或者用App开发器写个循环批量创建圆，再做三步操作：

1. 把所有圆全部选中执行“并集（Union）”，合并为单个固体域；
2. 用矩形域减去固体域（差集）；
3. 对差集结果做一次“清理几何”和“移除小碎屑”，把布尔操作产生的碎边清掉。

这步最容易翻车的地方是布尔操作失败。圆与圆太接近甚至相交时，COMSOL的几何内核会报退化错误。所以生成阶段把“最小间隙”当作硬约束，比在几何阶段补救有用一百倍。几何弄干净之后，流体域是一个带有大量圆形空洞的连续区域，网格剖分才能顺利。

3. 相场法核心原理与COMSOL参数设定

3.1 Cahn-Hilliard方程的基本形式

相场法的核心是Cahn-Hilliard方程。在COMSOL的“两相流，相场”接口中，约化的相场变量φ通常是-1到1，其中1代表一相，-1代表另一相，界面处有一个连续的过渡。方程的基本形式是：

∂φ/∂t + u·∇φ = ∇·(γ·λ/ε²·∇ψ) ψ = -∇·(ε²·∇φ) + (φ²-1)·φ

第一个式子左边是时间项和对流项，右边是扩散通量；第二个式子ψ是化学势对应的辅助变量。ε是界面厚度参数，λ是混合能密度，γ是迁移率。表面张力σ通过λ与其他参数关联，COMSOL内部会按你给定的表面张力换算混合能密度。双势阱项(φ²-1)·φ让体系自由能只有两个极小值，即两个纯相状态，从而驱动界面自然演化。

理解这个方程不需要陷入太深的数学。你把φ想成一个“相指示剂”，它在油里等于1，在水里等于-1，界面处经历从1到-1的连续过渡。流动带着它走，但扩散项又让它倾向于形成平滑的过渡层。这个过渡层的厚度就是ε控制的。流场作用于φ，φ通过物性插值和表面张力反过来作用于流场，二者耦合在一起。

3.2 与流动方程的耦合项：表面张力与物性插值

流场本身用不可压缩Navier-Stokes方程描述，但它不是单相流。因为φ在空间各点不同，密度和粘度必须做插值，常见形式是：

ρ = ρ1 + (ρ2 - ρ1) * (1 + φ) / 2 μ = μ1 + (μ2 - μ1) * (1 + φ) / 2

要注意你选择的φ极性：如果φ=1是油、φ=-1是水，那么插值公式里ρ2就是油的密度，ρ1是水的密度。弄反了会让界面附近的物性偏差，驱替行为可能完全改变。

表面张力的实现方式不是直接给一个张力值，而是把它转换成作用在界面上的体积力。化学势作用于φ的梯度，形成一个分布在界面层内的等效体积力。这就是相场法的一个优势：你不需要在几何上“找到”界面再施加法向力，而是通过场变量自动计算，拓扑变化时也不会出现力的不连续。

3.3 界面厚度和迁移率怎么取

这两个参数是相场模拟里最考验经验的地方，也几乎是所有调试工作的焦点。

界面厚度ε不是真实物理厚度，而是数值参数。它必须大于当地网格尺寸，否则界面层内没有足够节点，数值上无法分辨；但也不能太大，太大会让界面变成一条数微米宽的过渡带，表面张力作用区域被展宽，毛细效应失真。我的经验是：ε取最大网格尺寸的2~3倍，同时保证界面处至少分辨4~6层网格。比如最大网格1 μm时，ε设2 μm，感觉界面“有半毫米宽”的问题就少很多。

迁移率γ控制界面在对流驱动下的松弛速度。γ太小时界面僵硬，相位变化不灵活，收敛慢；γ太大则会引入数值扩散，界面的厚度随时间越变越厚。COMSOL默认值往往能用，但跑多孔介质这种复杂几何时会发现需要调。一个相对可行的起点是让γ约等于入口流速乘以特征网格尺寸，再根据界面宽度是否随时间扩张进行微调。

4. 驱替案例实操：一个完整的水驱油算例

4.1 几何尺寸、物性与工况参数一览

为了让描述具体化，我给出一个自己跑过的算例参数。二维流动域取100 μm × 100 μm，随机圆颗粒半径在5~15 μm之间分布，固体面积分数约0.4，对应孔隙率约0.6，颗粒间最小间隙设为2 μm。

水相作为注入相，密度1000 kg/m³，粘度0.001 Pa·s；油相作为被驱替相，密度800 kg/m³，粘度0.005 Pa·s，油水界面张力0.02 N/m。粘度比M=5，油比水粘，这是典型的有利于粘性指进的条件。入口速度设为0.01 m/s，此时毛细管数Ca = μ_w·U / σ = 0.001×0.01/0.02 = 5×10⁻⁴，属于毛细数偏高的范围。

如果你要做低毛细数对比，把入口速度降到0.001 m/s甚至更低，Ca降到5×10⁻⁵数量级，你会发现驱替前沿形态完全不同。这个对比非常值得做，它直接对应实验里“快注”和“慢注”两种策略。

4.2 边界条件、初始条件与润湿设定

边界条件我按如下方式设置：左侧为速度入口，给定法向流入速度，相场变量的流入值为水；右侧为压力出口，压力设0；上下壁面为润湿壁，接触角按工况设定，例如亲水壁取60°，中性壁取90°。

初始条件要注意技巧。整个孔隙域初始充满油（φ=1），但如果入口直接是水，就需要在入口处强行把φ翻成-1，会在极短时间内产生剧烈界面重构。稳妥的做法是先在左侧留出5 μm宽的水相区域，让初始界面以一条平滑的垂直线存在，然后再开始注入。初始界面的过渡用tanh型分布表达更自然，避免初始数值震荡。

润湿壁的接触角设定和φ的极性直接绑定。同样的“接触角60”，如果你的φ定义是1为油，壁面条件解释出来的是水相接触角还是油相接触角，结果可能正好差“镜像”。跑示例前最好先在简单平行通道里验证一下接触角的极性方向，再做多孔结构，不然你以为是亲水的结构实际跑出来是疏水的。

4.3 网格剖分与时间步长控制

网格剖分策略直接影响相场界面质量。我的做法是：全局最大单元尺寸1 μm，界面区域局部细化到0.3~0.5 μm，颗粒表面附近加边界层网格，避免速度场在固体边界附近分辨率不足。相场对网格质量的要求比单相流高得多，因为界面厚度只有2 μm，如果网格在这个区域太粗，φ的过渡会撕裂。

时间步长方面，相场模拟必须满足界面运动的稳定性约束，特征时间尺度很小。经验上初始时间步从1×10⁻⁵ s开始，运行一段后如果稳定可以让求解器逐步加大。比较省心的方式是打开自动时间步进，让COMSOL根据非线性收敛情况动态调整。计算总时长取0.02~0.05 s，也就是从注水开始到水突破出口所需的时间量级。

4.4 求解器设置与运行策略

瞬态求解器的选择上，相位场与流场同时求解比分为两步顺序求解更稳定，但计算成本高。如果是调试阶段，建议先关闭相场与流场的耦合中部分非线性项，或者用更短的模拟时间跑通流程；正式计算时再用全耦合或分离式迭代。

我还习惯做一个“空跑”验证：先不加表面张力，看单相流场是否合理，再打开相场接口，逐步叠加物理过程。这样每一步出问题都能定位。多孔介质入口可能出现速度突增，局部雷诺数如果超过层流范围就要降低入口速度或调整尺寸。本例Re约1，层流成立。

5. 模拟结果解析：从云图到饱和度曲线

5.1 相场云图与速度场的空间对应

跑完后第一件事就是看相场云图和速度流线叠加图。你要关注的不是“界面好不好看”，而是界面位置与局部速度场的空间对应关系。

在随机孔隙结构里，水相会优先选择大孔隙通道突破，流线在小喉道处汇聚成高速射流，把油相挤离。界面在喉道处被拉长、变薄，一旦超过临界曲率，就会被夹断形成孤立油滴。这些油滴悬浮在水相中，随流场移动直到被下一个喉道卡住。云图上你能直接数出有多少个这种被圈闭的油滴，它们的面积总和就是残余油的主要来源。

把速度场和相场一起看，还能发现低流速死角。有些孔隙入口虽然连通，但流速极低，水相前沿绕过它继续往大孔走，留下大片油区不动。这种“旁路”现象在规则排列结构里几乎不会出现，这就是随机结构带来的额外信息。

5.2 饱和度随时间演化与突破点

定量分析的第一步是计算水相饱和度随时间的变化。在COMSOL后处理中用表面积分，表达式对相场变量做判别，例如if(phi<0,1,0)，可以得到当前时刻水相面积占孔隙面积的比例。

典型的水驱曲线有几个特征段：初始阶段饱和度迅速上升，对应自由水前沿推进阶段，界面还没有到达主要出口；中段上升变缓，因为水相开始在多个孔隙通道间分流，部分通道已经突破但还有大片油区未被波及；出口突破后曲线进入平台期，饱和度增长非常缓慢，后续大量注水只提高少量采收率。这条曲线里的“突破点”是驱替策略的重要指标：突破越早，意味着指进越严重，非均质性起的作用越大。

建议同时绘制出口边界处水相分数随时间的变化曲线。出口出现水的时刻就是突破时间，很好判断。

5.3 指进模式的判断与物理讨论

指进模式可以用一个简化框架判断。低Ca且M>1时，通常是毛细指进占主导，界面在喉道处延迟突破，大量油被毛细力截留在小孔隙中；高Ca且M>1时，粘性指进占主导，水相以树枝状通道快速延伸，到达出口时大量油仍然留在未接触的孔隙区域；如果M<1且Ca适中，前沿相对平整，驱替效率较高。

实际案例中，M始终是固定的油水粘度比，改变Ca就改变驱替模式。你可以做一个入口速度从0.001到0.05 m/s的参数扫描，把每组工况的界面形态、突破时间、最终洗油效率放到一张表里。通常会发现存在一个折中的毛细数窗口：太高则指进早突破，太低则毛细力锁死残余油。这个窗口位置对孔隙尺寸分布非常敏感，换一组随机种子结论就可能偏移。所以做参数扫描时，建议同时跑3~5组随机孔隙结构，取结果的平均值而不是单次结果。

6. 常见问题排查与避坑经验

6.1 质量守恒漂移与虚假界面

相场模拟最常被质疑的就是质量守恒。虽然Cahn-Hilliard方程本身守恒，但离散误差和迁移率不匹配会造成φ的积分值缓慢漂移。排查思路是先算域内φ对时间的积分，看是否单调变化。如果持续漂移，优先检查：

• 迁移率是否过大，过大会在数值层面人为输送相组分；
• 界面厚度是否相对网格过小，导致界面处梯度不足、通量计算偏差大；
• 时间步长是否过大，影响了守恒精度。

如果问题依旧，可以适当加密界面网格，并把迁移率下调一个量级重新测试。我踩过一个教训：为了追求界面锐利，把ε设成比网格还小，结果完全没有守恒可言，界面出现“渗漏”。界面厚度必须和网格尺寸绑定考虑，不能只追求视觉锐利。

6.2 网格诱发振荡和几何退化导致的不收敛

不收敛问题里，很大一部分根源在几何和网格。随机圆颗粒如果存在相切或极近间距，布尔差集后留下的尖角区域网格剖分会生成特别差的单元，速度场在这里震荡，相场界面频繁抖动。排查手段是把几何放大看，找到那些狭缝和尖角，用间隙参数把它们消灭在几何生成阶段。

另一种不收敛是物理性的。高粘性比工况下，界面在喉道处被剧烈拉伸，局部曲率变化极大，如果时间步过大就会立刻发散。处理方法是减小时间步，打开自动时间步进，或者对相场方程单独使用隐式格式来提高对角占优性。

6.3 其他高频问题速查表

这里特别说一下接触线问题。在无滑移壁条件下，三相接触线处的应力是奇异性的，数值上几乎无法收敛。COMSOL润湿壁设置里通常允许指定一个小滑移长度，让接触线附近的流体有一定的切向滑移。滑移长度一般取网格尺寸的几分之一，既缓解奇异性，又不至于改变整体宏观流动。有人跑出接触线钉在原地不动，多半是滑移长度设为零或者过小。

最后再分享一个我自己的体会：相场模拟的调参顺序很重要。先把几何和网格做到干干净净，再调流动，最后才动相场参数。几何上有尖角、网格上有坏单元，后面所有的物理调参都会是白费力气。我第一次跑随机孔隙水驱油，调试时间大部分花在莫名其妙的网格和几何退化上，真正调相场参数只花了很小一部分。如果你准备长期做这类模拟，强烈建议把随机结构生成脚本做成一个独立工具，既能批量出几何，又能输出孔隙率、粒径分布统计，这样后面做参数扫描、趋势分析都顺手很多。