多孔介质里的渗吸听着像是油藏工程或岩土方向的老课题但真要拿数值方法把它算明白尤其是把裂缝和基质之间的竞争关系刻画出来很多人会在界面追踪和网格重构上卡几个星期。我自己的经验是COMSOL的相场法Phase Field配合移动网格能把“裂缝遇上基质孔隙”这种复杂润湿过程处理得相当顺手。这篇文章就从建模思路、参数设置、求解调试到结果解读完整梳理一遍多孔介质渗吸的相场模拟流程适合正在做两相流、微流控、油气运移、甚至混凝土吸水模拟的同行参考。开裂多孔介质里的渗吸难点不在“流体进去了”而在“裂缝和基质谁先谁后、吸进去多少、界面怎么演化”。传统方法要么用水平集追踪界面要么用VOF重构相界面但在裂缝尖端、孔隙喉道这种几何突变位置容易出现界面破碎或质量不守恒。相场法把界面当成一个连续过渡层来处理用扩散界面代替锐利界面相当于给界面加了一个“厚度”换来的是拓扑变化随便发生裂缝分叉、孔隙闭合、液桥断裂都能自然过渡不需要手动处理网格撕裂。这套逻辑配合COMSOL的层流两相流相场接口确实适合做孔隙尺度到裂缝尺度的渗吸研究。1. 为什么选相场法做渗吸模拟1.1 多孔介质渗吸的物理本质先捋一下物理图景。渗吸是润湿相在毛细力驱动下自发替换非润湿相的过程比如水进入油饱和的裂缝性砂岩或者水银进入微孔陶瓷。对多孔介质来说毛细力的大小由拉普拉斯方程决定Δp 2σcosθ / r其中σ是界面张力θ是接触角r是有效毛细半径。孔隙越小毛细力越大所以水总是优先往细喉道里钻。但裂缝存在时情况就变了。裂缝的开度通常比孔隙大一到两个数量级毛细力反而比基质孔隙小按理说裂缝应该不易吸水。然而裂缝的渗透率极高一旦裂缝壁面亲水水会沿裂缝快速推进然后再从裂缝壁面侧向渗入基质。这就形成了两种渗吸模式顺向渗吸counter-current和逆向渗吸co-current。前者是润湿相沿大通道进入同时非润湿相从同一通道排出后者是润湿相和非润湿相走不同通道。相场模拟要能自然呈现出这两种模式的竞争而不是靠人为设置界面。1.2 相场法相比传统界面追踪方法的优势传统VOF方法在COMSOL里也很常用它通过求解相体积分数输运方程来重构界面优点是质量守恒严格缺点是对网格各向异性敏感界面曲率计算容易带噪声尤其在裂缝壁面锐角处尖角处的曲率半径可能跟网格尺寸同阶导致寄生毛细流。水平集方法用符号距离函数描述界面拓扑变化处理得好但水平集函数在演化过程中需要周期性重新初始化否则距离场会失真质量守恒也不如VOF。相场法用的是Cahn-Hilliard方程把界面能作为自由能泛函的一部分通过化学势驱动相场变量φ演化。界面不再是零厚度而是一个有厚度的扩散层这个厚度由毛细宽度ε控制。正因为界面有厚度相场法在处理裂缝分叉、液滴合并、薄膜破裂等问题时天然稳定。代价是需要在界面区域加密网格因为ε通常取最小几何特征尺寸的一半左右网格分辨率不够就会让界面看起来“糊”一片。对于多孔介质这种几何尺度跨度大的模型需要在裂缝壁面和孔隙喉道附近做边界层网格这一点COMSOL的网格模块可以直接处理。1.3 COMSOL中相场接口的适用场景COMSOL从5.x系列开始就把相场接口和层流两相流捆绑在“层流两相流相场”物理场中6.x版本更是强化了移动网格与相场耦合的稳定性。它适合模拟不可压缩、互不相溶的两相流体在多孔介质中的动态侵入过程可以输出每相的速度场、压力场、相场变量分布还可以用积分探针统计含水饱和度随时间的变化。这个接口求解的是Navier-Stokes方程、连续性方程、Cahn-Hilliard方程以及后处理所需的体积分数方程。体积分数Vf (1φ)/2φ从-1变化到1分别代表非润湿相和润湿相。在COMSOL里设置相场初始值、接触角、表面张力系数后模型会自动把毛细力以体积力的形式耦合进动量方程。对于渗吸过程还可以关闭惯性项只保留粘性力和毛细力这就是准静态的毛细主导流动计算稳定性更好。2. 模型搭建前必须想清楚的三件事2.1 几何模型的取舍裂缝与孔隙怎么抽象多孔介质的真实孔隙结构极其复杂全尺度直接建模不现实。我通常的做法是分层次简化基础版用二维截面把基质设为规则排列的圆孔或方形孔隙网络裂缝用一条或多条弯曲的窄缝表示。二维模型足以揭示裂缝-基质的渗吸竞争机制计算量小适合做参数扫描。如果要做定量对比再升级到二维轴对称或三维单裂缝模型。裂缝的几何抽象要注意两点。第一是裂缝开度不能小于界面厚度的2倍否则相场界面会与裂缝壁面发生数值重叠接触角边界条件失真。我踩过这个坑开度设为1微米界面厚度取0.5微米结果裂缝里根本没出现完整的相场过渡层水直接“瞬移”过去了。第二是裂缝表面不能做得太理想光滑至少要加一点锯齿或波纹因为真实裂缝壁面的粗糙度直接影响接触角滞后和渗吸速度。2.2 控制方程与无量纲参数的选择搭建模型之前先手写一遍控制方程这一步能省掉后面大量调试时间。不可压缩两相流的连续性方程∇·u0动量方程ρ(∂u/∂t u·∇u) -∇p ∇·μ(∇u∇u^T) Fst Fext其中Fst就是相场法引入的界面张力体积力表达式为G∇φG是化学势。Cahn-Hilliard方程写作∂φ/∂t u·∇φ ∇·(M∇G)其中M是迁移率G λ(-∇²φ (φ²-1)φ/ε²)λ是混合能密度。很多人以为M和λ需要自己输入其实COMSOL提供了由表面张力σ和界面厚度ε自动计算λ和M的选项。关键是选择合适的界面厚度ε和迁移率M。这里有个无量纲数需要特别留意Péclet数Pe UL/M其中U是特征速度L是特征长度。Pe太小扩散项占主导界面被严重抹平Pe太大对流项占主导数值不稳定。COMSOL默认的M大约是1×10^-10 m²/s量级但对微米级孔隙来说特征速度只有毫米每秒Pe会很高需要把M调大一到两个数量级来保持稳定。这也是“相场模拟不稳定”最常见的原因之一。2.3 边界条件与初始润湿状态的设置边界条件的选择决定了渗吸模式。入口边界设置为润湿相压力边界比如给定入口压力略高于出口或者直接设置为Dirichlet压力等于毛细压力。出口边界设为0压力或充分发展流出。裂缝壁面和孔隙壁面的润湿性通过接触角设置实现COMSOL的“壁”边界条件里有“润湿壁”选项可以指定本征接触角θ。更关键的是初始状态。通常把非润湿相充满整个多孔介质润湿相只从入口边界向内侵入。如果基质的初始含水饱和度不为零需要在初始值里用解析函数或变量表达式定义φ的分布。这里强烈建议用平滑的阶跃函数例如tanh来初始化界面避免阶梯状的突变导致Cahn-Hilliard方程初始迭代发散。另外接触角滞后这个真实效应在COMSOL中可以近似通过设置前进角和后退角来实现但要注意相场法本身用的是平衡接触角动态接触角需要用户自定义表达式。如果只是做趋势研究固定接触角足够没必要一上来就引入动态接触角模型。3. 在COMSOL里一步一步搭出相场渗吸模型3.1 创建组件与选择物理场接口打开COMSOL 6.x新建一个二维模型选择“流体流动 两相流 层流两相流相场”接口这会自动添加层流接口、相场接口和“多物理场耦合”里的两相流相场耦合。模型树里会同时出现“层流”和“相场”两个节点后续设置都在各自节点下完成。研究类型选择“瞬态”。虽然渗吸过程可能在几十秒内就结束但不要用稳态研究因为相场演化本质上是时间依赖的稳态求解器往往收敛不到物理正确的界面形态。设置物理场的时候建议把“相场”节点的“迁移率”和“界面厚度”改成手动指定并用参数全局变量统一管理这样后续做参数扫描时非常方便。3.2 几何建模与网格划分的实操细节以二维单裂缝模型为例基质区域是50μm × 100μm的矩形中间有一条从入口延伸到50μm处、开度3μm的倾斜裂缝。孔隙简化为周期性排列的圆形半径4μm孔隙比约0.35。这样的模型几何不算复杂但网格划分决定了成败。建议先划分“边界层网格”。在流体域的裂缝壁面、孔隙壁面和入口边界处添加边界层属性第一层厚度设为界面厚度的1/5约0.05μm层数5层。然后全局自由三角形网格最大单元尺寸设为2μm最小单元尺寸设为0.02μm。等网格划分完可以做一个快速检查在主菜单里选择“统计”看看界面区域内即相场初始过渡层所在位置的网格尺寸是否小于ε/2。如果大于ε/2界面会被网格强行“钉扎”表现为界面不前进或者呈锯齿状。这时需要局部加密用“尺寸”节点在入口裂缝处单独设置更小的最大单元尺寸。3.3 相场参数、润湿边界与求解器设置参数设置界面里重点看几个表面张力σ水-空气体系在20°C下约0.072 N/m油水体系约0.03~0.05 N/m根据自己的体系设置。接触角θ设为60°表示润湿相是水固体壁面中等亲水。入口压力设为500 Pa这个量级对应的是微米级别孔隙的毛细压力。相场节点下的“初始化”子节点默认会计算初始界面。这里要指定初始相场φ的分布我用一个全局定义阶跃函数step1(t)配合解析表达式。表达式写成phi_init -tanh((x-5)/epsilon)表示在x5μm处有一条竖直的油水界面界面厚度由epsilon控制。注意tanh函数的自变量需要除以epsilon否则初始过渡层厚度和Cahn-Hilliard界面厚度不匹配会导致初始阶段出现虚假的界面运动。求解器方面选择“瞬态”研究时间步设为“范围(0, 0.0001, 0.5)”单位秒。求解器配置里推荐使用PARDISO直接求解器并把“恒定(鲁棒)Newton”改为“自动(Newton)”非线性残差相对容差设为0.001。相场方程是4阶偏微分方程直接法比迭代法更稳虽然内存占用大一些但二维模型完全能承受。3.4 后处理怎么把渗吸过程可视化算完之后最能说明问题的是相场云图。在“结果 二维绘图组”里绘制φ的云图用冷色到暖色的渐变显示界面层就落在过渡色带里。如果想提取湿润前缘位置可以用“派生值 体平均值”计算含水饱和度即1减去探测器体积分数平均值或者用“积分”算子计算φ0的面积占总孔隙面积的比例。更高级一点的操作是创建一个“全局评价”探针记录渗吸前沿x位置随时间变化的曲线。做法是定义变量wf min(x*(Vf0.5))思路比较复杂建议直接用“三维绘图组”里的“动画”功能逐帧查看界面推进过程。我常用的是“结果 动画”里的“播放”把帧保存为png再合成gif方便放在论文或汇报里。4. 裂缝对渗吸的影响几个典型结果解读4.1 裂缝开度与毛细力竞争先做一组开度扫描1μm、3μm、5μm、10μm。结果非常直观开度较小的裂缝1~3μm本身毛细力与基质孔隙相当渗吸前沿比较均匀整体呈活塞式推进。开度5μm以上时裂缝内毛细力明显小于基质孔隙但裂缝渗透率高水先沿裂缝快速推进到末端然后再横向渗入基质形成“倒灌”现象。这时含水饱和度云图上会看到裂缝通道先变蓝基质区域随后才慢慢变蓝。这个现象的本质是毛管数与毛细管的竞争。裂缝开度增大局部毛细压力降低但流动阻力降低得更快导致渗吸模式从“基质控制”转为“裂缝控制”。如果你的模拟目标是揭示裂缝性油藏的自吸采收率那结果可以直接改写成“小开度裂缝促进均匀渗吸大开度裂缝诱导局部窜流”这类结论。4.2 裂缝网络与基质渗透率的耦合单裂缝做完之后自然要扩展到裂缝网络。我在模型里加入两条交叉裂缝和三条分支微裂缝基质的渗透率通过改变孔隙半径和密度来调节。结果呈现一个有趣的临界行为当基质渗透率很低半径2μm密度小时裂缝网络是主要的渗吸通道基质几乎只靠裂缝壁面的侧向吸渗速度慢、最终采收率低提高基质渗透率后基质和裂缝同时吸渗总吸入量明显增加。这提醒我们一个建模要点多孔介质渗吸模拟的孔隙几何不能随便乱画。孔隙半径和密度的组合直接决定了基质毛细压力和渗透率进而影响裂缝-基质耦合行为。最好预先用Kozeny-Carman公式估算一下渗透率再和裂缝开度做量级对比确保几何设置合理。4.3 接触角与裂缝表面粗糙度的影响把接触角从30°扫到120°最明显的趋势是初期渗吸速率单调下降但最终吸水量在60°~80°之间出现一个峰值。原因很简单接触角过小时润湿相在裂缝内推进太快前锋还没形成侧向渗吸就已经冲出出口接触角适中时前锋推进速度和侧向渗吸速度刚好匹配基质能得到更充分的时间吸收液体。这个“最优接触角窗口”如果不做参数扫描根本发现不了。表面粗糙度的影响更为微妙。在裂缝壁面上加三角形锯齿锯齿高度0.5μm、间距2μm相比光滑壁面初期渗吸速度降低约15%但后期渗吸更平稳。这是因为粗糙度增加了接触角滞后阻碍了前缘在裂缝壁面的快速滑移同时也提供了更多的毛细钉扎点使界面推进更像“爬行”而不像“冲刺”。5. 常见问题与排查技巧实录5.1 相场轮廓发散或界面模糊典型表现φ云图出现斑点状噪声或者界面厚度明显大于设定值。原因通常是初始界面宽度与Cahn-Hilliard界面厚度不匹配或者迁移率M设置过大导致扩散过快。解决办法是先检查初始表达式确保界面过渡宽度和epsilon在同一量级其次把M降低1~2个数量级如果界面又变得过窄并对流输运不稳定再适当增大界面厚度epsilon到最小网格尺寸的3倍以上。另一个特别隐蔽的原因是网格在界面附近太粗。相场法要求界面内有至少3~5个网格节点否则曲率计算失真局部界面能就会产生虚假的毛细力从而把界面“推”弯。这时候比起盲目细化整个区域不如只细化界面可能经过的通道比如裂缝中轴线和孔隙喉道周围。5.2 收敛困难与时间步长控制出现“求解器直到最大迭代次数仍未收敛”时第一反应不是调容差而是看时间步长。Cahn-Hilliard方程本身是刚性方程时间步长过大会导致界面位置振荡。把初始时间步长设为1×10^-5 s量级再用“自由步长”让求解器自动调整。我还习惯开启“基于先前步长的向后差分公式”这能显著提高稳定性。如果依然不收敛检查边界条件里有没有压力突变。入口压力一下子给到500 Pa初始时刻流体加速度很大很容易震荡。稳妥做法是用一个斜坡函数ramp(t)让压力在0.01s内从0平滑升到目标值表达式写成p_in 500*min(t/0.01, 1)。这一步能解决八成以上的启动发散问题。5.3 质量不守恒与参数校准渗吸模拟特别在意含水饱和度计算是否守恒。如果算出来的总水量随时间漂移优先检查迁移率M是否太大。Cahn-Hilliard方程虽然守恒但数值离散会引入微小误差M越大界面处的人工通量越大。可以把M减小并加密界面网格两者结合往往能把质量误差控制在0.1%以内。另一个很多时候被忽略的因素是体积分数变量的后处理方式。COMSOL里相场变量φ和体积分数Vf并不完全等价Vf (1φ)/2在数值上会有截断误差。做质量守恒统计时应该用“积分”算子直接对Vf积分而不是对φ积分再换算。这个坑我曾经在数据导出时踩过曲线看起来不守恒实际上是后处理算错了。5.4 网格依赖性问题的处理做参数扫描时如果发现不同网格密度下结论不一致比如开度3μm裂缝的渗吸时间在网格加密后变化超过20%就要考虑网格依赖性问题。相场法的一个优势是解对界面厚度ε的收敛性而不是简单对网格。正确做法是固定ε为某个物理合理值例如最小孔隙半径的1/4然后同时细化网格和ε看结果是否趋近于同一个值。实际操作中我习惯用“参数化扫描”功能扫网格最大尺寸2μm、1.5μm、1μm、0.8μm并同步调整ε。如果最终计算结果趋于稳定再取最粗的网格方案做后续批量计算因为精度和速度的平衡点往往在2μm左右。这也是一个很实用的论文复现技巧审稿人问到网格无关性直接贴出收敛曲线即可。6. 写在最后的几个建议做相场渗吸模拟这几年我踩过的坑远不止上面列的那些但最深的体会是不要迷信COMSOL的默认参数。默认的迁移率和界面厚度是为了通用稳定性设计的用在多孔介质这种多尺度几何里往往不够必须自己根据毛细尺度和特征速度调校。另一个建议是模型越简单越能说明问题。很多初学者上来就想建真实岩心的三维孔隙网络模型结果网格分分钟上千万算一个案例要一周根本没法做参数分析。不如先建一个二维理想化模型把机制摸透再逐步增加复杂度。后面的论文或实际项目如果需要更高精度再考虑用CT扫描数据重建真实几何。最后补充一个小技巧后处理时把入口压力和裂缝开度同时做成参数扫描用COMSOL“一维绘图组”的“全局”功能画出含水饱和度随时间的family曲线你会看到渗吸从“裂缝主导”到“基质主导”的清晰转变。这个图放在任何汇报里都是很有说服力的结果。我在实际项目里常常需要连夜改模型赶汇报最后发现大部分时间不是花在物理上而是花在调网格和求解器上。如果这篇文章能帮你少熬两个晚上那我觉得这次分享就值了。 
长河君 · 2026/10/4 6:21:52 · 阅读 1.2k COMSOL相场法模拟多孔介质渗吸全流程解析
多孔介质里的渗吸听着像是油藏工程或岩土方向的老课题但真要拿数值方法把它算明白尤其是把裂缝和基质之间的竞争关系刻画出来很多人会在界面追踪和网格重构上卡几个星期。我自己的经验是COMSOL的相场法Phase Field配合移动网格能把“裂缝遇上基质孔隙”这种复杂润湿过程处理得相当顺手。这篇文章就从建模思路、参数设置、求解调试到结果解读完整梳理一遍多孔介质渗吸的相场模拟流程适合正在做两相流、微流控、油气运移、甚至混凝土吸水模拟的同行参考。开裂多孔介质里的渗吸难点不在“流体进去了”而在“裂缝和基质谁先谁后、吸进去多少、界面怎么演化”。传统方法要么用水平集追踪界面要么用VOF重构相界面但在裂缝尖端、孔隙喉道这种几何突变位置容易出现界面破碎或质量不守恒。相场法把界面当成一个连续过渡层来处理用扩散界面代替锐利界面相当于给界面加了一个“厚度”换来的是拓扑变化随便发生裂缝分叉、孔隙闭合、液桥断裂都能自然过渡不需要手动处理网格撕裂。这套逻辑配合COMSOL的层流两相流相场接口确实适合做孔隙尺度到裂缝尺度的渗吸研究。1. 为什么选相场法做渗吸模拟1.1 多孔介质渗吸的物理本质先捋一下物理图景。渗吸是润湿相在毛细力驱动下自发替换非润湿相的过程比如水进入油饱和的裂缝性砂岩或者水银进入微孔陶瓷。对多孔介质来说毛细力的大小由拉普拉斯方程决定Δp 2σcosθ / r其中σ是界面张力θ是接触角r是有效毛细半径。孔隙越小毛细力越大所以水总是优先往细喉道里钻。但裂缝存在时情况就变了。裂缝的开度通常比孔隙大一到两个数量级毛细力反而比基质孔隙小按理说裂缝应该不易吸水。然而裂缝的渗透率极高一旦裂缝壁面亲水水会沿裂缝快速推进然后再从裂缝壁面侧向渗入基质。这就形成了两种渗吸模式顺向渗吸counter-current和逆向渗吸co-current。前者是润湿相沿大通道进入同时非润湿相从同一通道排出后者是润湿相和非润湿相走不同通道。相场模拟要能自然呈现出这两种模式的竞争而不是靠人为设置界面。1.2 相场法相比传统界面追踪方法的优势传统VOF方法在COMSOL里也很常用它通过求解相体积分数输运方程来重构界面优点是质量守恒严格缺点是对网格各向异性敏感界面曲率计算容易带噪声尤其在裂缝壁面锐角处尖角处的曲率半径可能跟网格尺寸同阶导致寄生毛细流。水平集方法用符号距离函数描述界面拓扑变化处理得好但水平集函数在演化过程中需要周期性重新初始化否则距离场会失真质量守恒也不如VOF。相场法用的是Cahn-Hilliard方程把界面能作为自由能泛函的一部分通过化学势驱动相场变量φ演化。界面不再是零厚度而是一个有厚度的扩散层这个厚度由毛细宽度ε控制。正因为界面有厚度相场法在处理裂缝分叉、液滴合并、薄膜破裂等问题时天然稳定。代价是需要在界面区域加密网格因为ε通常取最小几何特征尺寸的一半左右网格分辨率不够就会让界面看起来“糊”一片。对于多孔介质这种几何尺度跨度大的模型需要在裂缝壁面和孔隙喉道附近做边界层网格这一点COMSOL的网格模块可以直接处理。1.3 COMSOL中相场接口的适用场景COMSOL从5.x系列开始就把相场接口和层流两相流捆绑在“层流两相流相场”物理场中6.x版本更是强化了移动网格与相场耦合的稳定性。它适合模拟不可压缩、互不相溶的两相流体在多孔介质中的动态侵入过程可以输出每相的速度场、压力场、相场变量分布还可以用积分探针统计含水饱和度随时间的变化。这个接口求解的是Navier-Stokes方程、连续性方程、Cahn-Hilliard方程以及后处理所需的体积分数方程。体积分数Vf (1φ)/2φ从-1变化到1分别代表非润湿相和润湿相。在COMSOL里设置相场初始值、接触角、表面张力系数后模型会自动把毛细力以体积力的形式耦合进动量方程。对于渗吸过程还可以关闭惯性项只保留粘性力和毛细力这就是准静态的毛细主导流动计算稳定性更好。2. 模型搭建前必须想清楚的三件事2.1 几何模型的取舍裂缝与孔隙怎么抽象多孔介质的真实孔隙结构极其复杂全尺度直接建模不现实。我通常的做法是分层次简化基础版用二维截面把基质设为规则排列的圆孔或方形孔隙网络裂缝用一条或多条弯曲的窄缝表示。二维模型足以揭示裂缝-基质的渗吸竞争机制计算量小适合做参数扫描。如果要做定量对比再升级到二维轴对称或三维单裂缝模型。裂缝的几何抽象要注意两点。第一是裂缝开度不能小于界面厚度的2倍否则相场界面会与裂缝壁面发生数值重叠接触角边界条件失真。我踩过这个坑开度设为1微米界面厚度取0.5微米结果裂缝里根本没出现完整的相场过渡层水直接“瞬移”过去了。第二是裂缝表面不能做得太理想光滑至少要加一点锯齿或波纹因为真实裂缝壁面的粗糙度直接影响接触角滞后和渗吸速度。2.2 控制方程与无量纲参数的选择搭建模型之前先手写一遍控制方程这一步能省掉后面大量调试时间。不可压缩两相流的连续性方程∇·u0动量方程ρ(∂u/∂t u·∇u) -∇p ∇·μ(∇u∇u^T) Fst Fext其中Fst就是相场法引入的界面张力体积力表达式为G∇φG是化学势。Cahn-Hilliard方程写作∂φ/∂t u·∇φ ∇·(M∇G)其中M是迁移率G λ(-∇²φ (φ²-1)φ/ε²)λ是混合能密度。很多人以为M和λ需要自己输入其实COMSOL提供了由表面张力σ和界面厚度ε自动计算λ和M的选项。关键是选择合适的界面厚度ε和迁移率M。这里有个无量纲数需要特别留意Péclet数Pe UL/M其中U是特征速度L是特征长度。Pe太小扩散项占主导界面被严重抹平Pe太大对流项占主导数值不稳定。COMSOL默认的M大约是1×10^-10 m²/s量级但对微米级孔隙来说特征速度只有毫米每秒Pe会很高需要把M调大一到两个数量级来保持稳定。这也是“相场模拟不稳定”最常见的原因之一。2.3 边界条件与初始润湿状态的设置边界条件的选择决定了渗吸模式。入口边界设置为润湿相压力边界比如给定入口压力略高于出口或者直接设置为Dirichlet压力等于毛细压力。出口边界设为0压力或充分发展流出。裂缝壁面和孔隙壁面的润湿性通过接触角设置实现COMSOL的“壁”边界条件里有“润湿壁”选项可以指定本征接触角θ。更关键的是初始状态。通常把非润湿相充满整个多孔介质润湿相只从入口边界向内侵入。如果基质的初始含水饱和度不为零需要在初始值里用解析函数或变量表达式定义φ的分布。这里强烈建议用平滑的阶跃函数例如tanh来初始化界面避免阶梯状的突变导致Cahn-Hilliard方程初始迭代发散。另外接触角滞后这个真实效应在COMSOL中可以近似通过设置前进角和后退角来实现但要注意相场法本身用的是平衡接触角动态接触角需要用户自定义表达式。如果只是做趋势研究固定接触角足够没必要一上来就引入动态接触角模型。3. 在COMSOL里一步一步搭出相场渗吸模型3.1 创建组件与选择物理场接口打开COMSOL 6.x新建一个二维模型选择“流体流动 两相流 层流两相流相场”接口这会自动添加层流接口、相场接口和“多物理场耦合”里的两相流相场耦合。模型树里会同时出现“层流”和“相场”两个节点后续设置都在各自节点下完成。研究类型选择“瞬态”。虽然渗吸过程可能在几十秒内就结束但不要用稳态研究因为相场演化本质上是时间依赖的稳态求解器往往收敛不到物理正确的界面形态。设置物理场的时候建议把“相场”节点的“迁移率”和“界面厚度”改成手动指定并用参数全局变量统一管理这样后续做参数扫描时非常方便。3.2 几何建模与网格划分的实操细节以二维单裂缝模型为例基质区域是50μm × 100μm的矩形中间有一条从入口延伸到50μm处、开度3μm的倾斜裂缝。孔隙简化为周期性排列的圆形半径4μm孔隙比约0.35。这样的模型几何不算复杂但网格划分决定了成败。建议先划分“边界层网格”。在流体域的裂缝壁面、孔隙壁面和入口边界处添加边界层属性第一层厚度设为界面厚度的1/5约0.05μm层数5层。然后全局自由三角形网格最大单元尺寸设为2μm最小单元尺寸设为0.02μm。等网格划分完可以做一个快速检查在主菜单里选择“统计”看看界面区域内即相场初始过渡层所在位置的网格尺寸是否小于ε/2。如果大于ε/2界面会被网格强行“钉扎”表现为界面不前进或者呈锯齿状。这时需要局部加密用“尺寸”节点在入口裂缝处单独设置更小的最大单元尺寸。3.3 相场参数、润湿边界与求解器设置参数设置界面里重点看几个表面张力σ水-空气体系在20°C下约0.072 N/m油水体系约0.03~0.05 N/m根据自己的体系设置。接触角θ设为60°表示润湿相是水固体壁面中等亲水。入口压力设为500 Pa这个量级对应的是微米级别孔隙的毛细压力。相场节点下的“初始化”子节点默认会计算初始界面。这里要指定初始相场φ的分布我用一个全局定义阶跃函数step1(t)配合解析表达式。表达式写成phi_init -tanh((x-5)/epsilon)表示在x5μm处有一条竖直的油水界面界面厚度由epsilon控制。注意tanh函数的自变量需要除以epsilon否则初始过渡层厚度和Cahn-Hilliard界面厚度不匹配会导致初始阶段出现虚假的界面运动。求解器方面选择“瞬态”研究时间步设为“范围(0, 0.0001, 0.5)”单位秒。求解器配置里推荐使用PARDISO直接求解器并把“恒定(鲁棒)Newton”改为“自动(Newton)”非线性残差相对容差设为0.001。相场方程是4阶偏微分方程直接法比迭代法更稳虽然内存占用大一些但二维模型完全能承受。3.4 后处理怎么把渗吸过程可视化算完之后最能说明问题的是相场云图。在“结果 二维绘图组”里绘制φ的云图用冷色到暖色的渐变显示界面层就落在过渡色带里。如果想提取湿润前缘位置可以用“派生值 体平均值”计算含水饱和度即1减去探测器体积分数平均值或者用“积分”算子计算φ0的面积占总孔隙面积的比例。更高级一点的操作是创建一个“全局评价”探针记录渗吸前沿x位置随时间变化的曲线。做法是定义变量wf min(x*(Vf0.5))思路比较复杂建议直接用“三维绘图组”里的“动画”功能逐帧查看界面推进过程。我常用的是“结果 动画”里的“播放”把帧保存为png再合成gif方便放在论文或汇报里。4. 裂缝对渗吸的影响几个典型结果解读4.1 裂缝开度与毛细力竞争先做一组开度扫描1μm、3μm、5μm、10μm。结果非常直观开度较小的裂缝1~3μm本身毛细力与基质孔隙相当渗吸前沿比较均匀整体呈活塞式推进。开度5μm以上时裂缝内毛细力明显小于基质孔隙但裂缝渗透率高水先沿裂缝快速推进到末端然后再横向渗入基质形成“倒灌”现象。这时含水饱和度云图上会看到裂缝通道先变蓝基质区域随后才慢慢变蓝。这个现象的本质是毛管数与毛细管的竞争。裂缝开度增大局部毛细压力降低但流动阻力降低得更快导致渗吸模式从“基质控制”转为“裂缝控制”。如果你的模拟目标是揭示裂缝性油藏的自吸采收率那结果可以直接改写成“小开度裂缝促进均匀渗吸大开度裂缝诱导局部窜流”这类结论。4.2 裂缝网络与基质渗透率的耦合单裂缝做完之后自然要扩展到裂缝网络。我在模型里加入两条交叉裂缝和三条分支微裂缝基质的渗透率通过改变孔隙半径和密度来调节。结果呈现一个有趣的临界行为当基质渗透率很低半径2μm密度小时裂缝网络是主要的渗吸通道基质几乎只靠裂缝壁面的侧向吸渗速度慢、最终采收率低提高基质渗透率后基质和裂缝同时吸渗总吸入量明显增加。这提醒我们一个建模要点多孔介质渗吸模拟的孔隙几何不能随便乱画。孔隙半径和密度的组合直接决定了基质毛细压力和渗透率进而影响裂缝-基质耦合行为。最好预先用Kozeny-Carman公式估算一下渗透率再和裂缝开度做量级对比确保几何设置合理。4.3 接触角与裂缝表面粗糙度的影响把接触角从30°扫到120°最明显的趋势是初期渗吸速率单调下降但最终吸水量在60°~80°之间出现一个峰值。原因很简单接触角过小时润湿相在裂缝内推进太快前锋还没形成侧向渗吸就已经冲出出口接触角适中时前锋推进速度和侧向渗吸速度刚好匹配基质能得到更充分的时间吸收液体。这个“最优接触角窗口”如果不做参数扫描根本发现不了。表面粗糙度的影响更为微妙。在裂缝壁面上加三角形锯齿锯齿高度0.5μm、间距2μm相比光滑壁面初期渗吸速度降低约15%但后期渗吸更平稳。这是因为粗糙度增加了接触角滞后阻碍了前缘在裂缝壁面的快速滑移同时也提供了更多的毛细钉扎点使界面推进更像“爬行”而不像“冲刺”。5. 常见问题与排查技巧实录5.1 相场轮廓发散或界面模糊典型表现φ云图出现斑点状噪声或者界面厚度明显大于设定值。原因通常是初始界面宽度与Cahn-Hilliard界面厚度不匹配或者迁移率M设置过大导致扩散过快。解决办法是先检查初始表达式确保界面过渡宽度和epsilon在同一量级其次把M降低1~2个数量级如果界面又变得过窄并对流输运不稳定再适当增大界面厚度epsilon到最小网格尺寸的3倍以上。另一个特别隐蔽的原因是网格在界面附近太粗。相场法要求界面内有至少3~5个网格节点否则曲率计算失真局部界面能就会产生虚假的毛细力从而把界面“推”弯。这时候比起盲目细化整个区域不如只细化界面可能经过的通道比如裂缝中轴线和孔隙喉道周围。5.2 收敛困难与时间步长控制出现“求解器直到最大迭代次数仍未收敛”时第一反应不是调容差而是看时间步长。Cahn-Hilliard方程本身是刚性方程时间步长过大会导致界面位置振荡。把初始时间步长设为1×10^-5 s量级再用“自由步长”让求解器自动调整。我还习惯开启“基于先前步长的向后差分公式”这能显著提高稳定性。如果依然不收敛检查边界条件里有没有压力突变。入口压力一下子给到500 Pa初始时刻流体加速度很大很容易震荡。稳妥做法是用一个斜坡函数ramp(t)让压力在0.01s内从0平滑升到目标值表达式写成p_in 500*min(t/0.01, 1)。这一步能解决八成以上的启动发散问题。5.3 质量不守恒与参数校准渗吸模拟特别在意含水饱和度计算是否守恒。如果算出来的总水量随时间漂移优先检查迁移率M是否太大。Cahn-Hilliard方程虽然守恒但数值离散会引入微小误差M越大界面处的人工通量越大。可以把M减小并加密界面网格两者结合往往能把质量误差控制在0.1%以内。另一个很多时候被忽略的因素是体积分数变量的后处理方式。COMSOL里相场变量φ和体积分数Vf并不完全等价Vf (1φ)/2在数值上会有截断误差。做质量守恒统计时应该用“积分”算子直接对Vf积分而不是对φ积分再换算。这个坑我曾经在数据导出时踩过曲线看起来不守恒实际上是后处理算错了。5.4 网格依赖性问题的处理做参数扫描时如果发现不同网格密度下结论不一致比如开度3μm裂缝的渗吸时间在网格加密后变化超过20%就要考虑网格依赖性问题。相场法的一个优势是解对界面厚度ε的收敛性而不是简单对网格。正确做法是固定ε为某个物理合理值例如最小孔隙半径的1/4然后同时细化网格和ε看结果是否趋近于同一个值。实际操作中我习惯用“参数化扫描”功能扫网格最大尺寸2μm、1.5μm、1μm、0.8μm并同步调整ε。如果最终计算结果趋于稳定再取最粗的网格方案做后续批量计算因为精度和速度的平衡点往往在2μm左右。这也是一个很实用的论文复现技巧审稿人问到网格无关性直接贴出收敛曲线即可。6. 写在最后的几个建议做相场渗吸模拟这几年我踩过的坑远不止上面列的那些但最深的体会是不要迷信COMSOL的默认参数。默认的迁移率和界面厚度是为了通用稳定性设计的用在多孔介质这种多尺度几何里往往不够必须自己根据毛细尺度和特征速度调校。另一个建议是模型越简单越能说明问题。很多初学者上来就想建真实岩心的三维孔隙网络模型结果网格分分钟上千万算一个案例要一周根本没法做参数分析。不如先建一个二维理想化模型把机制摸透再逐步增加复杂度。后面的论文或实际项目如果需要更高精度再考虑用CT扫描数据重建真实几何。最后补充一个小技巧后处理时把入口压力和裂缝开度同时做成参数扫描用COMSOL“一维绘图组”的“全局”功能画出含水饱和度随时间的family曲线你会看到渗吸从“裂缝主导”到“基质主导”的清晰转变。这个图放在任何汇报里都是很有说服力的结果。我在实际项目里常常需要连夜改模型赶汇报最后发现大部分时间不是花在物理上而是花在调网格和求解器上。如果这篇文章能帮你少熬两个晚上那我觉得这次分享就值了。