PFC 怎么加流体？单向、双向、接第三方，一次讲清

“PFC 能不能算渗流？”

这个问题我被问过很多次。它不好回答，因为答案不是一个"能"或"不能"，而是你在哪一层问——

• 
  
• 想知道一堆颗粒在静水里怎么沉？能，而且很准。
• 想知道水流怎么把沙床冲出一个坑？能，需要一个外部求解器。
• 想知道水怎么把岩石压裂开？抱歉，这条路在原理上就走不通。

这三件事的差别，全在"耦合机制"上。这篇文章就把机制从头拆一遍。

1. 先看全貌：一共四条路

一句话：PFC 自己不含流场求解器。它只负责"颗粒这一半"——算孔隙率、算拖曳力；流场那一半，要么你直接给定，要么交给外部软件去解。

PFC 里做流固耦合，翻来覆去就是四条路：

 

图 1 四条路的区别只在两个问题上：流场谁来解？信息是单行道还是双行道？

先别急着记，我们一层层往下拆。

2. 机制的核心：把流域切成一堆"大方格"

不管是哪条路，原理的底子都一样，叫粗网格法（unresolved CFD-DEM）。

它的做法可以用一句话概括：把整个流域切成一堆大方格，网格边长比颗粒大好几倍。

 

图 2 上：流域被切成大方格，颗粒被"平均"进去。下：挑一格放大，每一格只记三个数。

每个格子只记三个数：

分工是这样的：PFC 干两件事——把颗粒的位置翻译成每个格子的 ε，再把颗粒受到的拖曳力算出来。外部求解器干一件事——解流体力学方程，把 v 和 p 更新掉。两边隔一会儿交换一次数据。

2.1 为什么非得这么"糙"？

因为细算算不起。

要把颗粒周围的水算清楚，网格必须细到粒径的十分之一左右。一颗 1 mm 的颗粒就要占几千个格子；一捧沙就是几亿到几十亿个格点。

所以粗网格法做了一个交易：用经验公式换计算量。

这个交易的代价很具体——拖曳力不是"解"出来的，是"拟合"出来的：

 

图 3 上：公式拆成四块。中：拖曳系数随流态变化。下：颗粒越挤，同样的流速给出越大的力。

PFC 用的这套（Di Felice 关联式）本质上就两项：

拖曳力 = 单颗粒受力 × 密堆补偿

f_drag = f0 × ε(-χ)

• 
  
• f₀：一颗孤零零的球在无边界水里受到的力。它和"水与颗粒的相对速度"的平方成正比——跑得越快，力涨得越猛。
• ε^(−χ)：密堆补偿。颗粒越挤（孔隙率越小），周围颗粒对水流的"堵塞"越厉害，力会被放大。图 3 最下面那条曲线就在说这件事：孔隙率从 1.0 降到 0.4，同样的流速给出的力能放大十几倍。

这一点想清楚，很多争论就没了：粗网格法的精度天然受限于经验公式，它不可能比公式本身更准。所以用它的时候，先把"跟解析解对一遍"当纪律，不要直接上工程结论。

3. 单向耦合：水推颗粒，颗粒不挡水

搞懂了格子，再来看"单向"和"双向"就好懂了。

 

图 4 左边是单向：只有箭头从水指向颗粒，反方向那条线被划掉了。右边是双向：两个方向都在传。

单向耦合的意思很直白：水怎么流，早就定好了；颗粒只是被动挨推，它挡不挡水，没人管。

听起来很粗糙？但有一类问题它就够用，而且能对到解析解。PFC 官方自带一个这样的算例：

三颗半径 1、1.5、2 mm 的球，密度 2000 kg/m³，在密度 1000 kg/m³、黏度 1.5 Pa·s 的静水里自由下沉。流体速度全程人为设成 0——把流场当成一个"不动的水池"。

 

图 5 上：本机实跑的沉降曲线与公式手算的终速（同色虚线）。下：跟两种理论值比，差多少。

这张表值得多看两眼。为了拿到它，我们做了两件事：

1. 
   
2. 把概念反过来验一遍：照官方手册里那套 Di Felice 关联式，用 Python 独立实现（不调用 PFC），解出稳态终速；
3. 在本机把官方那个算例跑了一遍，导出沉降曲线。

结果是四位有效数字全同。换句话说，PFC 内部用的确实就是手册上写的那套公式，一个系数都没动。

而如果拿纯黏性的 Stokes 律去比，PFC 会"偏大" 0.8~3%。这不是误差失控，而是方向正确的偏差：Di Felice 的 Cd 公式含惯性项，Stokes 是纯黏性的 24/Re。球越大、Re 越高，惯性效应越明显，偏差就越大——趋势本身也说明公式用对了。

什么时候可以只做单向：颗粒体积占比很低（稀相）、颗粒根本挡不住水流的问题——单颗粒沉降、稀薄悬浮、低速搅拌。反过来，只要颗粒成堆、能明显改变流道，就得做双向。

4. 双向耦合：你推我，我也推你

双向耦合就是把那条被划掉的线接回来。

颗粒挡住水 → 水绕开 → 颗粒受力跟着变，这个圈得转起来。

 

图 6 一个耦合步里的六件事，转一圈。

每一步的数据交换是这样的：

外部求解器收到前三个之后要做的事，翻译成大白话就是：

在流体力学方程里加两项——一项是"这个格子里只有 ε 比例的空间能走水"，另一项是"这里插了一堆颗粒在拽水"。

写成公式，就是在标准 N-S 方程后面挂上 − f/n（f 是反作用力场，n 是孔隙率）。像 OpenFOAM 的 pyDemFoam 求解器，动量方程就是：

bashddt(U) + div(phi, U) - laplacian(nu, U) - f/n = -grad(p)

前两项是常规的时间项与对流项，第三项是黏性项，最后那个 f/n 就是颗粒插 进来的那一脚。

4.1 两个不能省的旋钮

双向耦合最容易踩的坑是发散。上游那份 OpenFOAM 求解器是一阶欧拉 + 中心差分，几百颗颗粒强耦合的时候，不松不弛几乎必然炸掉。所以有两个参数是设计的一部分，不是可选项：

• 
  
• 松弛因子 r = 0.5：孔隙率、拖曳力、速度三路都要松。新值和旧值各取一半，不让数据一步跳变。
• CFL 自适应时间步：水越快，步子越小。而且得留安全余量——因为用来估步长的是 PFC 手里那个滞后一步的速度场，会低估真实的 Courant 数。

5. 接第三方：一根 TCP 线，两种角色

到这里有个绕不开的问题：外部求解器怎么跟 PFC 说上话？

答案朴素得有点意外——一根 TCP socket。

 

图 7 PFC 当服务端先起，外部求解器当客户端后连。握手一次，之后每个耦合步交换六个场。

关键是角色不能搞反：PFC 侧调用 p2pLinkServer，它会一直阻塞到客户端连上来。所以顺序永远是 PFC 先起，求解器后连。

而客户端那一侧，其实任何 Python 环境都能当——pip install itasca 装上 p2pLinkClient 就行，不要求那台机器上装了 PFC。这也意味着：求解器跑在另一台电脑上、跑在 Linux 里、跑在 WSL 里，都行。

5.1 官方给了一个"傻客户端"

这是我觉得官方文档里最实用的一招。

它给了一个完整但什么都不算的客户端示例：每个耦合步老老实实把数据收下来，然后返回三组零（压力全 0、压力梯度全 0、速度全 0）。PFC 那边照常收下、照常继续跑。

为什么要这么干：因为裸 socket 是按序收发，没有消息边界——两边约好了先发三个、再收三个，全靠顺序对齐。谁多写一个、少收一个，就全线错位，而且不会报错，只会卡死。先用一个"返回零"的傻客户端把通道跑通，再往里面填真求解器，就把"通信对不对"和"物理对不对"这两件事拆开了。

5.2 这套机制，我们自己接了一次

光看图不够，我们照着官方文档把两侧都写了出来：

• 
  
• PFC 侧：program call 一个 Python 脚本，用 p2pLinkServer 起服务端；
• 外部侧：不装 PFC、不装任何 Itasca 的包，纯标准库 + numpy，按协议把客户端写出来。

20 个耦合步跑完，两侧零错误：

两个力看着不一样，其实是一回事：客户端收到的是体积力（已经除以单元体积和流体密度），2.2165e-03 ÷ 0.125 ÷ 1000 = 1.773e-05 ✓ 对得上。

还有一件事被顺带验证了：客户端每步都返回全零的压力和速度，跑完 20 步（0.02 s）后，球的速度是 −0.0980 m/s。而自由落体此时应该是 9.81 × 0.02 = 0.196 m/s —— 实实在在的一半。

一半不是巧合：颗粒密度 2000、流体密度 1000，cfd buoyancy on 给出的净加速度正好是 g × (1 − 1000/2000) = 4.905，乘 0.02 s 就是 0.0981 ✓ 浮力项确实生效了。

写客户端时踩到的一个坑：协议是「4 字节类型码 + 紧凑 payload」，两边都得按紧凑格式打包。而 Python 的 struct.pack("id", ...) 在默认模式下会为了对齐在 int 和 double 之间插入 4 字节填充——PFC 把填充当成了数据，读出来的浮点是 0。症状很有欺骗性：握手正常、前两个数据正常、第三个开始是垃圾。打包时写死小端紧凑格式（< 前缀）就对了。

接上真正的 OpenFOAM 之后，跑出来的是这样的东西——水流冲刷颗粒床（3D 双向耦合，547 颗碎石，2400 个流体单元）：

 

图 8 中跨剖面切片，颜色是流速。从 t=0.5 s 到 3.5 s，入口流速从 0.38 涨到 1.50 m/s，可以清楚看到脊体被冲蚀、颗粒在床面输运。

6. 2D 怎么办？自己写一个求解器

前面讲的这套，有一个硬限制：PFC 的 CFD 模块只有三维版本。

这不是"没配置好"，是物理上的缺席——模块的文件只有 3D 版（module3dccfd*.dll），PFC2D 里 from itasca import cfdarray 直接抛 ImportError，官方的 CFD 手册全篇 106 次提到 PFC3D、0 次提到 PFC2D。

所以二维要做流固耦合，路只剩一条：自己写求解器。

这条我们走过，用的是 FISH 自带的 socket，接一个自己写的 LBM（格子 Boltzmann） 求解器：

 

图 9 2D 水力劈裂：孔压场（蓝色→红色）与破裂时刻的拉裂纹网络（绿线）。

机制上和前面"接第三方"那节是一回事，只是"外部求解器"换成了自己写的：FISH 侧开 socket 客户端，Python 侧当服务端，每个交换周期把颗粒的位置、半径、速度发过去，求解器推完流场再把逐颗粒的力回传。

这里有一个反直觉的点，值得单独说：

交换频率是稳定性参数，不是性能参数。
两次交换之间，PFC 把流体力当成常数施加在颗粒上。实测：每 200 个 LBM 步才换一次力时，颗粒在恒定力下持续加速到 0.3 m/s，直接把流场击穿；改成每 20 步换一次就稳了。调这个数的目的不是"省时间"，是让颗粒的响应时间短于交换周期。

LBM 还有另一条更轻的用法——只做单向：让流体穿过一个固定的颗粒堆积，测出渗透率。

 

图 10 三维 LBM 单向测渗透率，并与 Kozeny-Carman 经验式对照。

为什么渗透率要用外部求解器测：PFC 自带的 CFD 模块不含流场求解器。官方那个 Darcy 算例里的渗透率，是从 Kozeny-Carman 经验式算出来的，不是流出来的。要一个真正从孔喉几何里"流"出来的渗透率，就得自己解流场。

7. 这张饼能画多大

回过头看，三条路线其实是一条连续谱，分档标准就一个：流体网格相对颗粒有多大。

 

图 11 从上到下，网格越来越粗，精度越来越依赖经验，但能算的规模越来越大。

而 resolved 这条高保真路，其实也有天花板。三条约束一联立就能算出来上限：

• 
  
• 格子速度要小于 0.1（不可压精度，马赫数约束）
• 格子黏度要大于 0.005（BGK 稳定）
• 网格要细到颗粒间隙的 1/2（密堆间隙约 0.05 倍粒径）

消掉时间步，得到的物理流速上限是每秒零点几米这个量级。

一条硬边界，必须说清楚
PFC 的 CFD 模块有三条明示假设：流体单元比颗粒大、单元内流体属性分段线性、单元不移动。
把它翻过来读——流体单元不移动、且必然大于颗粒，意味着它没法描述"水钻进一条比颗粒还窄的裂缝里把它撑开"。
所以：水力劈裂不能走这条路，这是原理问题，不是调参问题。

8. 不信就自己算一遍

前面讲的机制，能验的我们都验了。

验证一：单向耦合的终速（图 5）。把手册里那套 Di Felice 公式用 Python 重写一遍，和 PFC 实跑对照——四位有效数字全同。

验证二：密堆补偿这一项。这个更值得说，因为它是"经验公式"里最像玄学的部分。

做法很直接：一个测试球固定在流体单元正中央，流体速度恒定；然后往同一个单元里塞障碍球，把孔隙率从 0.996 一路压到 0.614，看测试球受到的力怎么变。

 

图 12 横轴是孔隙率，纵轴是拖曳力相对"单颗粒"的放大倍数。橙线是公式，蓝点是实测。

四个测点，七位有效数字全同。

验证三：通信回路（第五节）。PFC 侧用官方 p2pLinkServer，外部客户端是我们从协议写起的（不依赖任何 Itasca 的包）。20 个耦合步跑完，两侧数据逐位对得上，零错误——而且球的运动量正好等于"去掉浮力"后的净加速度，说明交换进来的场真的被用上了。

这个验证的意义：它说明 PFC 的"经验"不含糊——手册上写什么公式，代码里就执行什么公式，连孔隙率压到 0.6 这种密堆情形都不打折扣。
反过来也提醒一件事：这个力是从公式来的，不是从流场解出来的。公式适用范围之外的情形（超大颗粒、极端孔隙率、非球形颗粒），它不会自己变聪明。

9. 到底怎么选

诚实声明
① 机制与公式来自本机 PFC 6.0 官方文档（cfd_theory / cfd_third_party / one_way_cfd 等章节）的一手取证。
② 图 5 与图 12 是本机实跑的结果，不是照抄文档：图 5 跑的是 Itasca 官方单向耦合算例（数据在 oneway_speed.his，与官方文档给的终速一致）；图 12 是自建算例。两者的理论对照曲线都由 verify_drag.py 按手册公式独立手算，全程不调用 PFC。
③ 第五节的 socket 实测也是本机跑通的：PFC 侧脚本 + 自写的外部客户端，20 个耦合步零错误，脚本与日志在 CFD耦合验证/socket/。
④ 图 8、图 9、图 10 分别来自本机的 PFC3D×OpenFOAM 冲刷算例、LBM 水力劈裂算例、LBM 三维渗透率算例。
⑤ PFC 6.0 没有颗粒尺度的管网流（domain/pipe）模型，5.0 时代那套水力压裂脚本搬不过来——这一点我们单独取证过。
⑥ CFD 模块需要 license，model configure cfd 带星号（需额外授权）；model configure fluid（FLAC3D 那条）不带星号。

VibePFC 社区 · 一起把 AI 玩进 PFC
我成立了一个 VibePFC 社群，专门聊 AI 怎么落地到 PFC / 数值模拟——让大模型读懂模型、自己改脚本、自动排查"看不懂"的问题。收一个 50 元的门槛费。
想进群的同学，加 QQ 763388012，备注"vibePFC"即可。