XFlow格子玻尔兹曼方法模拟多孔介质两相流：从毛细管渗吸到参数标定实操

格子玻尔兹曼方法（LBM）这个名词，在油藏工程和渗流力学圈子里早就不是什么新鲜概念了，但真正能用工程级软件把它跑起来、而不是自己吭哧吭哧写两万行C++代码的人，确实不算多。我当年用XFlow做毛细管自发渗吸的模拟，说白了就是被逼出来的——甲方想看水怎么在微观孔隙里自己吸进去，传统VOF搞不定随机骨架里的接触线移动，而自己写LBM又面临生命周期和调试的双重压力。最后把宝押在XFlow的格子玻尔兹曼内核上，总算把双相流动、接触角滞后、渗吸深度与时间的关系这些硬骨头都啃下来了。

这篇文章就是冲着"能复现"三个字来的。如果你正在评估XFlow能否用于多孔介质两相流模拟，或者你已经打开软件但面对一堆物理模型参数不知道该从哪下手，那么下面这些从原理到实操、从参数标定到结果验证的内容，应该能帮你少烧掉一沓笔记本。

1. 为什么毛细管自发渗吸把传统CFD逼到了墙角

1.1 界面追踪的先天之痛

先说我为什么一开始没走常规路线。传统CFD里处理两相流，主流方案其实就是VOF（流体体积法）和Level Set（水平集法）。VOF的思路是把每个网格单元里"水的体积占比"当作一个标量场来追踪，这个思想很直观，实现也不复杂，所以OpenFOAM和Fluent里都有成熟求解器。但VOF有个绕不过去的坎：它必须对相界面做几何重构，也就是说，你每一时刻都要专门去算界面在哪个位置、什么形状，这需要额外引入界面压缩项和几何算法。当你的计算域是一块充斥着几百个随机喉道和孔腔的复杂多孔介质，这些重构计算就会变得异常繁琐，而且界面一拉伸就断裂。

Level Set的思路稍微优雅一点，它用一个距离函数来隐式表达界面，数值稳定性比VOF好，但它在相界面附近的体积守恒差，尤其是充满曲率极小的毛细管壁面附近，液面推进一点点都可能导致质量误差累积。

1.2 渗吸问题的硬约束：接触线

真正逼近传统CFD极限的，其实不是"怎么算界面"，而是"怎么算接触线"。毛细管自发渗吸的驱动力是弯月面处的Laplace压差，而这个压差的数值直接取决于动态接触角和毛细管半径。接触线在固体壁面的移动涉及微观分子尺度的物理化学作用，而宏观连续性方程对这部分无能为力。传统方法最常见的做法是人为指定一个接触角边界条件，但碰到真实岩心那样表面粗糙度极高的孔隙壁面，接触角实际上是时变的，这会让速度场和压力场在壁面附近出现非物理的奇异性。

1.3 LBM的介观哲学

格子玻尔兹曼方法能走通这条路，核心在于它的世界观完全不同。LBM不做界面追踪，也不解宏观N-S方程，它站在介观尺度上，用一堆离散速度分布函数的碰撞和迁移来重现流体行为。你设定一个"粒子群在格点上碰撞后趋向平衡态"的规则，然后让计算自动涌现出界面的曲率、表面张力、接触线滑移。

这种哲学对两相流简直是量身定制的：

• 界面不需要重构，它在密度场（或序参量场）里自然形成，自带梯度；
• 壁面处的相互作用力可以通过伪势模型非常自然地引入，接触角不再是"硬性边界条件"，而是"材料属性"；
• 多孔介质边界即使复杂到每个格点都是固体节点，LBM的反弹格式（bounce-back）也能天生处理。

我这么描述可能有点抽象，你试着把LBM想象成一场大型的交通仿真：每个格子是一个路口，分布函数是不同方向的车辆数，碰撞规则是红绿灯时序，迁移规则是车辆移向相邻路口。你不需要去追踪每一辆车的轨迹，只要规则设对了，早晚高峰的拥堵波自然而然就出现了。两相流模拟就是这个逻辑的高级版——不同的"车辆颜色"代表不同的流体相。

1.4 为什么项目最终落在XFlow上

既然LBM这么好，那自研代码是不是更可控？是的，学术圈大量论文都是自己写Lattice代码，而且D3Q19、D3Q27这些方案在开源社区也都有现成模板。但我的项目周期不允许我从头折腾并行计算和网格生成，尤其需要的是工程化后处理和模型参数的第一手解释。XFlow走的是无网格路线，它不需要传统意义上的"网格划分"，几何体直接以曲面形式导入，软件自动在计算域内填充格子，这对随机生成的多孔介质模型来说省掉了巨大的前处理工作量。

此外，XFlow的物理模型库中直接内置了多相流模型，能让你在图形界面里设置表面张力、壁面润湿性、接触角这些参数。它把学者们争论不休的微尺度参数封装成了工程参数，对做应用研究而不是做算法研究的人来说，这是效率上的降维打击。

2. XFlow里的两相模型：从Shan–Chen理论到可调参数

2.1 伪势模型的本质

XFlow的多相流模型其核心本质是Shan-Chen（SC）伪势模型，这是1993年提出来的经典框架。SC模型不去显式追踪界面，而是在每个格点上定义一个"伪势"函数，这个势函数依赖于局部的密度或序参量。邻近格点之间因为伪势差异会产生一个相互作用力，这个相互作用力在宏观尺度上等效为：同种流体之间趋向于抱团，不同种流体之间趋向于分离。界面张力就是这种微观力统计平均的宏观表现。

要注意，Shan-Chen模型里没有显式的"表面张力系数向量"，表面张力是从等温状态方程和相互作用强度参数$\psi(\rho)$中自发涌现出来的。这带来一个工程问题：你没办法直接填入"水的表面张力等于72.8 mN/m"，你得先做参数标定。这也是很多第一次用XFlow做两相流的人卡壳的地方。

2.2 XFlow中界面参数的三种设定逻辑

在XFlow的界面里，与表面张力相关的参数设置大致可以分成三类：

第一类是直接指定表面张力系数的（如果你选择连续的表面应力模型），这种情况比较少见，因为软件仍然需要在模拟过程把表面张力转化为体积力或格子力，内部做了复杂的换算。

第二类是接触角和壁面润湿性的设置。XFlow支持在几何壁面上赋予不同的润湿属性。这个润湿属性实质上是在壁面格点额外施加一个偏向某一相的伪势相互作用，导致接触角从90度向亲水或疏水方向偏移。如果你设置的壁面参数偏向水相，水就会更容易铺展；偏向气相，水就会收缩成球。

第三类是状态方程的参数。XFlow让用户可以修改流体的状态方程来改变密度比和压力密度关系。在LBM里，密度比不能设成真实水-气那样的夸张比例（真实水气密度比接近1000，SC模型稳定范围一般几百以内），所以工程上要做一个折中：让两相密度比尽量高一些，同时保证数值不发散。

2.3 元胞自动机式的直觉理解

这其实是理解LBM两相流最核心的思维转变：不要把界面张力当成一个"施加在已知界面上的外力"，而要把界面张力当成"微观规则运行后的统计结果"。就像你在站台上看人群，人多的区域人和人之间会自然形成有一定距离的排队间距，这个间距不是某个管理员画的黄线，而是所有人遵循同一套"别挤我"的微观规则涌现出来的宏观现象。Shan-Chen的界面也是一样，它没有一个追踪界面的网格，但界面却永远在那个"势函数梯度最大"的地方存在。

2.4 格子单位与物理单位的桥梁

LBM的一切计算都在格子单位中运转，而你关心的结果（渗吸速率、压力）是物理单位。这里的换算关系必须在建模前写清楚，否则结果一出来你都不知道自己在看什么。

XFlow软件本身会把大部分换算工作拦下来，它允许你在界面中直接定义流体的物理属性，然后在内部完成格子与物理单位的映射。但作为使用LBM的人，你至少要明确两件事：一是特征长度对应关系，即你所设定的1个格子长度实际代表多少微米，这决定了毛细管内径在计算域里占据的格点数；二是特征时间对应关系，这影响你能用多大的时间步长进行瞬态模拟。我的经验是，在XFlow中先用小模型把单位换算的逻辑彻底打通之后，再上完整的多孔介质模型，不然在结果阶段反复调整物理单位会让人崩溃。

3. 把"水的天性"标定进计算域：表面张力与接触角的实操校准

3.1 第一个案例：静止液滴测试

任何LBM两相模型，第一关永远是液滴测试。你不用设置复杂的边界，只需要在一个不大的计算域内放一个静止的球形液滴，然后观察它在表面张力作用下是否保持稳定形状，且是否满足Young-Laplace定律。XFlow有个好处，它可以直接用后处理功能输出压力场，你就能方便地读出液滴内外的压差。

具体操作上，我建议按如下步骤进行：

1. 建立一个小型三维计算域，例如每个方向100个格点，中央处放置一个球形水相液滴，其余区域为气相。
2. 液滴半径至少取25个格子，太小的半径会由于离散曲率误差导致压力场数值乱跳。
3. 设定两相密度比（例如设置成10到50），表面张力参数先行给一个默认值。
4. 运行直到液滴内外的密度场达到稳态，以总动能的减小为标准判断是否达到稳态。

然后，输出中心轴线上的压力分布，考察界面两侧的压力差$\Delta p$与理论预测公式$\Delta p=2\sigma / R$（对球形液滴）的匹配程度。如果压差偏大，说明表面张力参数偏高；压差偏小则反之。这个校准过程的目的不仅在于"设置一个正确的数值"，更关键的是找到格子单位下表面张力系数与物理表面张力之间的乘数因子。

3.2 接触角标定：从90度出发

接下来是接触角标定。XFlow允许你在壁面上设置不同的润湿性参数，我在实际操作中发现，最好在每个相似度的模拟系列开始前都自己做一个接触角参照测试。方法很简单：

1. 在计算域底部设置一个平面壁面。
2. 在壁面上放置一个半球形的液滴。
3. 运行状态达到平衡后，通过后处理输出密度等值面，量取液滴在壁面上接触角。

这个测试几乎是所有多孔介质毛细管模拟的"地基"，因为你后续所有关于亲水/疏水渗吸能力的结论，全部建立在壁面参数与接触角的对应关系上。如果这里标定不准，后面所有的渗吸深度、渗吸速率都会偏。

我在一次硅酸盐玻璃毛细管的模拟里，标定目标接触角是35度，第一次跑出来只有55度。反复排查才发现，壁面的润湿参数被我在模型复制时不小心覆盖回了默认值。这种错误极其典型——软件的几何树界面里"表面属性"和"求解器设置"经常被混用，新手尤其容易在看参数面板时漏掉"某个表面是否启用了物理属性"这一开关。

3.3 黏度的格子单位换算

格子玻尔兹曼方法里面，流体的运动黏度与松弛时间（relaxation time）直接相关，关系式为$\nu = c_s^2 (\tau - 0.5)\Delta t$。如果你在XFlow里直接输入物理单位下的黏度值，软件会负责换算，但你心里得有本账：

• 当$\tau$接近0.5时，数值黏度很小，模拟趋向于不稳定，界面可能波动；
• 当$\tau$过大时（例如大于1.0），数值耗散太大，界面运动会变得过于钝化。

我的建议是尽量把两个相的松弛时间控制在0.5到1.5之间，宁可把水相和气相的黏度比缩小一些，也不要让任何一个相的$\tau$值接近下限。我在模拟液态水-氮气体系时，真实黏度比大约为50，XFlow初始配方直接给了个默认黏度比达到了200，结果初始步长还没有跑到1024步，压力场就开始振荡。用XFlow做两相流模拟，你最终必须接受的现实是：它的范式面向"快速工程评估"，你必须在物理精度和计算稳定性之间做一个工程化的取舍。

3.4 密度比与计算效率的权衡

高密度比是LBM两相流模拟的永恒主题。对于水-气体系，密度比越接近真实值，数值稳定性越差；密度比压得过低，则"水"和"气"两相的密度差异不够显著，界面两侧惯性力对比失真。做自发渗吸模拟时，主要的驱动力是毛细张力而非重力，所以密度比对渗吸速率影响有限。基于这个前提，我认为把密度比设置在50到200之间是合情理的选择。密度比选用200时，水的格子密度约2.0，气相约0.01，收敛性和物理真实性达到了一个在工程上可以接受的平衡。

4. 搭一座自发渗吸的可视化实验台

4.1 单毛细管模型的几何设置

从最简单的算例开始——一根直的毛细管嵌入一块固体基体，毛细管一端与"液池"连接，另一端开放到气相。计算域我可以描述如下：

• 整体尺寸：长度300格点，截面100×100格点；
• 毛细管半径10到15个格点，位于计算域中央，沿z方向贯穿；
• 固体基体的材质设置为"亲水壁面"，接触角目标值由前面的液滴测试标定的壁面参数给定。

这里有个小技巧：毛细管入口和出口处不要做成突变截面，用一个斜坡或者锥形过渡，否则入口效应会干扰渗吸早期阶段的速度场，让初期渗吸深度$H(t)$曲线偏离理论的$t^{1/2}$规律。

4.2 边界条件怎么给才"自发"

自发渗吸的核心要求是：水相没有外部压力驱动，纯粹靠弯月面的Laplace压差把水拉进毛细管。实际操作中，常见做法是设为两端压力均为0（参考压力一致），出口侧允许开放流动，入口侧连接一个无限大液池。在XFlow中，这意味着要合理设置计算域的开放边界（open boundary）。

我不建议在这类瞬态渗吸模拟中把出口直接设成固定压力边界，因为初期弯月面形成时会产生压力脉冲，固定压力边界会把这种物理过程给反射回去。开放边界配合适当的缓冲区更加贴近实际的实验条件，而且缓冲区还能抑制压力波反射的数值伪影。

4.3 初始状态的设置技巧

初始化时，毛细管内预先填满气相，计算域入口侧的一段区域（也就是液池区域）填充水相。注意初始时刻不能把弯月面设置得太陡，否则界面处的密度场在第一步就会被强烈的伪势力冲散，产生所谓"数值爆炸"。稳妥的办法是给初始相分布一个平滑的过渡带，这个过渡带宽度大约为5-10个格子即可。

另一个容易忽视的地方是初始速度场。我见过太多人在LBM模拟里直接把速度初值设为0，这个逻辑在单相流里没问题，但在两相流里，初始时刻界面附近的密度梯度会立刻产生很大的局部力，若速度场完全没有缓冲，就会导致早期压力振荡。我的做法是让模拟前2000步采用较小的全局时间步长，或者干脆在最初的几百步内把表面张力从0线性增加到目标值，让系统"软启动"。XFlow如果支持时间相关的自定义场，这是一个非常实用的技巧。

4.4 多孔介质骨架构造的注意事项

做完单毛细管，就可以往真实的多孔介质模型迈进了。XFlow接受STL或STEP格式的几何文件，所以你可以用开源工具生成随机孔隙网络模型然后导出几何。我当时用的是竞争性生长算法生成的堆积球体骨架，再提取内部孔隙空间作为流体域。

入口和出口附近需要设置集流腔（manifold），否则多孔介质内的流动会受入口几何的强烈干扰。集流腔的尺寸至少是最大孔径的5倍。计算域的边界如果采用周期性边界条件，要注意XFlow在多相流下的周期性边界与单相流实现有细微差异——周期性边界下不相连通的两相压力自由浮动，此时验证体积守恒就格外重要。

5. 结果判读：渗吸深度、Washburn定律与易错陷阱

5.1 渗吸前沿的提取方法

模拟跑完之后，第一步是从后处理数据中提取渗吸前沿的位置。我的习惯是取密度等值面$(\rho_{liquid}+\rho_{gas})/2$与毛细管中心轴线相交的点作为瞬时弯月面位置。将不同时刻的前沿位置连成一条$H(t)$曲线。

这里要特别提醒：XFlow后处理的默认颜色映射未必能清楚显示两相界面，建议自定义密度值的user field。颜色映射范围一定要手动锁定，因为气液两相的密度差异可能高达两个数量级，自动映射会把低密度相的细节压到一片混沌之中。

5.2 和Washburn方程面对面

毛细管自发渗吸的经典理论是Lucas-Washburn方程：

$$ H(t) = \sqrt{\frac{\gamma , r \cos\theta}{2\mu} , t} $$

这个公式表达了渗吸前沿位置与时间的平方根关系，其中$\gamma$是表面张力，$r$是毛细管半径，$\theta$是接触角，$\mu$是液体黏度。

我把模拟结果和这个理论公式做了对比，发现早期阶段（初始几十个时间单位）模拟结果会出现与平方根关系偏离的情况，这主要是因为入口处的流动需要一定时间才能建立起稳定的弯月面形状。在毛细管半径小、初始界面曲率与平衡曲率差异大的情况下，这种偏离会更加明显。当你撰写研究报告时，一定要把这段过渡期单独标注出来，不要为了拟合理论曲线而隐藏它。

5.3 我的实测参数表

以直径20个格子、长度200个格点的单毛细管为例，一组能够稳定跑通并得到合理结果的重要参数如下：

对照组数据表明，对于亲水毛细管（接触角35°），模拟渗吸深度指数（即$\log H$对$\log t$的斜率）约为0.49，与理论值0.5高度吻合，这个结果在工程上可以说相当漂亮。

5.4 质量守恒是最容易被忽略的评判标准

做了这么多LBM模拟，我的最深体会是：很多人只盯着界面形状好不好看，却忽视了质量守恒。一个两相流模拟如果液体总质量随时间减少5%，那么你就不能用它来定量分析渗吸速率。在XFlow中，每次运行结束后我都建议输出两相的总质量随时间的变化曲线，确认其波动在百分之零点几的范围内。如果出现明显的质量漂移，通常原因包括：时间步长过大导致局部数值不稳定、壁面处浸润层厚度不足，或者开放边界压力反射导致非物理的泄漏。

在我自己跑的一组多孔介质随机骨架案例中，第一次质量漂移率高达8%，排查两小时后发现问题出在前处理的几何修复上：STL文件里有一些肉眼看不见的自交面，导致LBM在局部格点上的反弹格式出现错误，把流体"吞"掉了。用软件自带的几何清理工具重新修复后，质量守恒恢复到了0.5%以内，这一课值两小时。

5.5 从单毛细管到真实岩心的路还长

最后要坦诚一句，单毛细管与真实多孔介质之间，还存在一个很大的鸿沟。真实岩心的孔隙网络有着宽泛的孔径分布、顽固的角隅滞留效应、以及复杂的互连性，这些都是模拟难以完美复现的因素。但XFlow的价值恰恰在于它能让你快速淘汰那些明显错误的物性假设，从而把时间留在真正值得深入研究的物理问题上。

如果你手头正好要开始这类模拟，我给你的建议是不要一上来就追求和实验数据的完全一致。先用单毛细管模型跑通全流程，记录你设置的所有参数与对应的结果特征，然后加一个孔或者加一个弯，逐步逼近你要研究的真实结构。这种从简到繁的路子，比起直接堆一个巨大的随机结构模型，然后面对无数个不知从何调起的参数要省心得多。