PFC2D 接 OpenFOAM：官方走不通，我给它自建了一条 2D 耦合线
1. 先说结论
PFC 官方的 CFD 耦合模块只随 PFC3D 发布，PFC2D 没有。想让 2D 的颗粒跟 OpenFOAM 的流场双向耦合，中间有四件事要自己写：孔隙率场、二维拖曳力、力回传、单位厚度约定。
我全写了，纯 numpy，然后在一根 1.2 米的水槽里做了一场 30 秒的沙坡冲刷。
先看结果——30 秒被压缩成几十秒：床体从一座完整的梯形山脊，被水流一步步削平，颗粒向前摊开。下半屏是跟拍放大，能看到单颗的移动。
结果：全程 15000 多步稳定不发散，床没有被冲散，而是被慢慢削平；入口流速从 0 线性爬到 0.4 m/s，动画里能直接读出颗粒开始移动的那个流速——约 0.08 m/s 开始蠕变，0.21 m/s 明显削顶。整条链已经沉淀成可复现的技能。
2. 为什么 PFC2D 官方连不了
PFC 的 CFD 模块（cfdarray）同时干两件事：算孔隙率、算拖曳力，再通过 TCP socket 跟外部求解器交换数据。问题是它只存在于 3D。四条证据，都是实机测出来的：
安装目录里只有 module3dccfd*.dll，没有 2D 版；PFC2D 里敲 model config cfd，直接报Bad conversion；from itasca import cfdarray，ImportError；官方手册 ccfd 分册通篇只有 PFC3D。
但 PFC2D 留了两样能用的接口：p2pLinkServer（socket 通信层，跟 3D 同款）和 ballarray（批量读写颗粒、施加外力）。通讯和驱动力都现成，缺的是夹在中间的那层物理——socket 帮你传数据，但没人替你算孔隙率和拖曳。
路线其实有三条，选型一句话：要 3D 走官方路径（已有验证路，本篇不讲）；要 2D 但要求 OpenFOAM 特有的东西（真实密度、大变形流场、后续接湍流），走本篇这条自建线；要 2D 又想快、零编译，用 LBM 那条（D2Q9，已验证）。这篇讲第二条。
3. 方案A：中间层四件事，全部手写
上图左框就是要自己补的全部内容：3D 侧 cfdarray 帮你做的，2D 侧一个都没，只能手写。
3.1 第一件：孔隙率场
每个颗粒都是一个圆，落到哪几个矩形网格单元里，就按圆和矩形的重叠面积算：
ε = 1 − Σ(重叠面积) / 单元面积
这是解析几何，不用采样，误差 0.05% 以内。唯一要注意的是：单元里塞进两颗大球时 ε 会算成负数，送进求解器之前必须钳在 0.05 以上——负孔隙率会让源项立刻发散（实测 14 步内爆）。
3.2 第二件：二维拖曳力（全篇最关键）
这是最关键的一件。3D 的拖曳公式不能直接搬到 2D——二维有 Stokes 悖论（无 界域里定常 Stokes 解不存在）。所以低雷诺数段用 Oseen 二维圆柱解：
F = 4πμU / (0.5 − γ − ln(Re/8))，γ = 0.5772
中高雷诺数段才切回工程常用的 Di Felice 球公式（参考长度取 2r）。分段点放在 Re = 1.5。
为什么必须分段？看这张图：低雷诺数区，Di Felice 的阻力系数从两个方向同时爆表（Cd ∝ 1/√Re 跑到 600，孔隙率修正 ε^(−χ) 又放大近 400 倍），乘积让拖曳虚高一万倍，第一次数据交换就炸。Oseen 那条才是二维真正的低速极限，Di Felice 那条是工程近似——必须在 Re=1.5 处握手交班。
3.3 第三件：力回传
颗粒受的拖曳力，按单元聚合成体积力送回求解器：
f_cell = −Σ 拖曳 / (Δx·Δy·ρ)
浮力显式加：(πr²)·ρ·g。整个系统是单位厚度的 2D 约定——颗粒质量 m = ρπr²，一切"体积"用 πr²。这套约定要贯穿孔隙率、质量、拖曳、力回传四处，错一处，量纲就塌。
3.4 第四件：稳定性钳制（不写必发散）
显式双向耦合有个硬约束：一步之内流速的变化，不能超过驱动它的相对速度。这个约束两头都要压——颗粒侧压拖曳力，流体侧压源项。为什么流体侧也要压？因为流体看到的源项要除以局部孔隙率，在颗粒密集处被放大 6–20 倍，只压颗粒一侧，流体先炸（实测只压一侧 0.13 秒就发散，两侧都压稳稳跑完 30 秒）。再配合三路松弛、CFL 限步和发散守卫，才算完事。
4. 怎么跑：两步走，看一场 30 秒的沙坡冲刷
4.1 第一步：颗粒堆成沙坡
227 颗圆盘按三角晶格精确摆成一个中央梯形沙坡，粒径 3.0–4.5 mm，最大最小粒径比 1.5，颗粒摩擦系数 2.0，初始零重叠——特别强调零重叠，因为它直接决定后面读数可不可信。
然后按浮重沉降到位。什么叫浮重：颗粒泡在水里感受到的自重不是 ρgV，而是 (ρ−ρf)gV，所以沉降时的重力加速度取 g·(1−ρf/ρp) = 6.11 m/s²。这一步沉降到位非常重要：否则耦合一开始泥沙突然"变轻 38%"，整个床体会弹一下，前几秒的画面全是假动态——初版就栽在这里。后面所有数字，都是在修掉这个问题之后跑的。
多说一句为什么不直接 6 秒冲完：慢爬升要跑 30 秒（约 15000 步，实机一百来分钟），听起来很贵。但只有慢爬升能把"什么时候开始动"变成可读数——流速上限压到 0.4，阈值就落在中段，前后都有对照。
4.2 第二步：耦合，入口慢爬升
OpenFOAM 在 WSL2 里起一个 120×20 的单层网格，入口流速用一张表从 0 线性爬到 0.4 m/s，爬 30 秒。
刻意选慢爬升而不是一口气冲上去，是因为这场实验的目的不是"看它被冲垮"，而是读起动流速：什么时候床面开始动、什么时候明显侵蚀，每个阈值都对应一个入口流速，直接标在动画上。
这张图是全场最重要的数据：上半是颗粒平均速度（对数轴），下半是床质心和脊顶高度，横轴同时标了时间和入口流速。起动点读法：颗粒平均速度在 6 秒附近抬头（入口 0.08 m/s），16 秒附近加速（入口 0.21 m/s），脊顶同步开始掉。三个量同时变，才是真起动，单个速度抖一下不算。
5. 结果：床被慢慢削平，而不是冲散
三帧放在一起看：前期床形几乎不动，中期脊顶被削低，后期整床摊成薄层、开始有颗粒到下游出口。
关键数字：
最有意思的是起动之后床没散。颗粒没有一粒粒被卷走，而是表层不断被剥离、脊顶越削越低、颗粒向下游摊开——这是标准的渐进型侵蚀形态，不是突发失稳。fric=2.0 的床在这组流速下很稳，真正把颗粒送到出口需要更高的流速或更长的时间。
换句话说：起动 ≠ 冲走。0.08 m/s 让它动起来（蠕变），0.21 让它掉肉（削顶），但要把它整体搬走，0.4 爬 30 秒还不够——这三个阶段在动图里是肉眼可分的。
6. 一个需要交代的边界：这套流速不能直接当"真水"用
做冲刷的都想要一个"起动流速"的绝对值，比如"0.2 m/s 就冲跑了"。这里必须诚实：本算例用的流体运动黏度是 ν = 3e-3 m²/s，比水大了三个数量级——这是为了保证粗网格 + 颗粒尺度数值稳定，属于加厚黏性做法。
所以在"0.4 m/s 才能削顶"这类数字里，颗粒雷诺数全程都很低（拖曳全程走 Oseen 段），绝对值的意义要打折扣。它的价值在于量级、机制和趋势：慢爬升读阈值的方法、渐进而非突发的侵蚀形态、friction 主导的床体上限——这些结论不受黏度缩放影响。要贴近真水的绝对起动流速，下一步应该把黏度拆回原位并细化网格，这是明确的后续方向，不是糊弄过去的坑。
7. 沉淀
整条链已经封装成 pfc-openfoam-2d 技能：算例生成器、两个耦合端、出图脚本、以及这次踩过的全部坑（沉降门禁、拖曳分段、云图转置）都写进去了，改个床形就能复跑你自己的算例。
PFC2D 没有官方路，但自己铺一条轮子完全够用。改天把黏度拆回真水再跑一遍，看绝对阈值怎么动。
