PFC2D 接上自由面：颗粒入水，水下堆积自己长出来
1. 先说结论
上一篇文章讲了怎么给 PFC2D 自建一条 2D 耦合线（没有官方 CFD 模块，孔隙率、拖曳、力回传全部自己写，做了一场沙坡冲刷）。那条线有一个边界：流场没有自由面——算不了「颗粒从水上落入水里」这类问题，因为空气和水不分家。
这一篇就是把自由面接上去：流场升级成 VOF 两相（空气 + 水，官方 interFoam 求解器内核），然后让 80 颗颗粒分批入水，在池底逐层堆成一座锥形水下堆积。
看全过程（水面是蓝的，红点是颗粒）：
三个可量化的结果：
水下锥形堆真的长出来了：80 颗颗粒 71% 集中在投放带（中心 0.3 m 内），堆顶高出池底约 4.6 cm，剖面中心明显隆起——不是抛洒平铺； 水体积全程守恒：液相体积漂移 +5.8%（预算 10%），水没有凭空消失； 水面溅起后自己回落：溅起的不平整度峰值 ~4 mm，8 秒末收敛到 ~2 mm，液体最大流速从 0.83 m/s 衰减到 0.15 m/s，颗粒末态几乎静止（平均速度 0.001 m/s）。
2. 为什么入水堆积必须要有自由面
单相链路里没有空气-水界面：颗粒落下去，流场只有一个「水」，溅起来的水花没有重力拉它回来，动图里就是水花一直飘着、永不回落——这不是物理。
要让它像真的，流场必须是两相：上面空气、下面水，界面能变形、能溅起、能回落。做法是把 OpenFOAM 官方的 interFoam 求解器内核（MULES 界面压缩 + 静水压力 p_rgh）编译进耦合链，取代之前的手写单相动量方程。
两相动量方程用守恒形式，这是与单相最本质的区别：
守恒形式：水和空气的密度差 800 倍，动量方程必须对 ρ 加权，否则密度间断两侧的动量不守恒，界面会「放气」； p_rgh（去静压的动压）：把重力从压力里拆出来单独处理（ −g_h ∇ρ），静水平衡不依赖压力迭代的初值，空池静水面的流速能收敛到 1e-4 而不是振荡；界面压缩项： ∇·(μ∇U)之外还有 MULES 的通量限幅 + 界面压缩波，把几毫米厚的物理界面压进 2–3 个网格单元。
颗粒作用在流体上的是一个源项 f（PFC 每步算好拖曳力后按单元聚合发过去）。这个源项有两处关键细节：
按水密度归一，不按局部混合物性：f 已经是「颗粒对流体的加速度」，乘回水密度才还原成守恒方程要的每体积力； 不除以局部孔隙率：密堆单元里 n 会掉到 0.05，除下去源项直接放大 20 倍（实测 |U| 从 4.6 冲到 6.1 发散）。
图1｜PFC2D（Windows 侧 DEM）与 OpenFOAM（WSL2 里）通过 socket 每步交换力场、孔隙率、流速和水体积分数 α
3. 水怎么表示：α 场 + 界面压缩
两相流场用水体积分数 α 描述：α=1 是纯水，α=0 是纯空气，界面在 0→1 过渡的地方。界面由官方的 MULES 算法压缩保持锐利（几个单元宽），不会被数值扩散糊成一团。
α 的口径有一个容易搞错的点：它是「空隙中的水分数」，不是「单元里的水体积分数」。一个单元里如果塞满了颗粒，水的绝对体积会变小，但那是颗粒排开的孔隙空间，不是水丢了——α 只在「空隙」这个维度上跟踪水。这个口径是体积守恒能守住的关键（后面图里水体积全程漂移 <6%）。
颗粒侧怎么知道自己在水里还是在空气里？取颗粒中心所在单元的 α，按 α 加权拖曳和浮力：全湿（α=1）用水的密度和黏度，全干（α=0）用空气的（拖曳几乎为零），中间平滑过渡——这就是颗粒穿水面时「减速入水」的物理来源。
拖曳和浮力共用同一个润湿过渡函数（smoothstep，α 在 0.05–0.95 之间过渡），保证入水瞬间两个力同步、不出现「先被拖一下、再被浮一下」的台阶。颗粒-流体相对速度的雷诺数分段用两套拖曳模型：低速段用二维 Oseen 圆柱解（避开二维 Stokes 悖论），高速段切 Di Felice 工程关联式。
图2｜四个时刻：颗粒从水面以上分批落入，穿过界面减速，池底的堆逐层长高。虚线是初始水位
4. 颗粒怎么堆成锥：分批投放是关键
一开始用「一次性把 80 颗倒进去」（哪怕收窄成一条窄带）也堆不成锥——所有颗粒同时砸入，激起横向回流，把轻颗粒在水底卷平了：末态铺满 0.76 m 宽的池底，堆高只有 2.8 cm。
改成分批投放（每 0.24 秒投 4 颗，前一批先沉底、后一批才到），堆就开始一层层往上长：
图3｜颗粒集中在投放带的比例：分批后 68%→71% 集中在中心 0.3 m，一次性只有约一半
图4｜堆顶能长多高：分批后堆顶从 2.8 cm 长到 4.6 cm
图5｜末态剖面：中心 0.4–0.6 m 明显隆起成锥，边缘散落少量（持续投料的水面波动推开少量颗粒，是真实物理）
这个「分批投放」机制已经是耦合脚本的通用功能：一个 release.csv 调度表（每行一颗颗粒的释放时刻与坐标），耦合主循环按模拟时钟逐批 ball create。换投放节奏、颗粒数、池子尺寸只需改这张表/生成器参数。
技术细节：为什么批次要「从水面附近低处释放」而不是从高处砸进去。颗粒入水激起的水流强度正比于入水动能，批次任务从 0.27 m 高处落下 → 入水速度 ~2.3 m/s；改成水面以上 4 cm 释放 → 入水速度 ~0.9 m/s，动能小一个量级，溅起的水流就不再能把之前的堆冲平。投放节奏、释放高度、批次颗粒数三者一起决定了「堆的成长速度 vs 被冲散的风险」，是本案例最重要的可调参数组。
5. 把耦合本身做稳：三道措施
显式双向耦合（PFC 算力、OpenFOAM 算流，各自独立时间步推进）天然容易振荡发散，这个两相版本加了三个专门的稳定措施：
无过冲限幅（双向）：每步施加的拖曳力不得超过「让颗粒/流体在一步内追上相对速度」的物理上限——颗粒侧和流体侧各压一次（两处都压，只压一侧坏得更快）； 三路松弛：孔隙率、拖曳力源项、流速三个量在交换时都做 新 = 旧×(1−R) + 新×R的混合（R=0.4），滤掉一个耦合步长的锯齿；发散守卫：CFD 侧每步检查 |U| 与 NaN（ python3 -u实时输出），一旦超限立即停机——宁可少跑几步，不做无意义的长跑。
这三道是上一篇文章（单相冲刷）就验证过的组合，两相版沿用了并针对入水冲击放宽了流速阈值（入水瞬间的合法流速峰值比冲刷高一个量级）。
6. 水面会平稳吗：阻尼是真实存在的
溅起的水为什么最终会平静？三个机制叠加：
粘性扩散：水的层流黏性项把界面动能变成热； 界面压缩的数值耗散：MULES 在界面附近限幅通量，亚网格高频波纹每步被低通一次； 重力恢复力：静水压力场把倾斜的水面拉回水平。
实测衰减曲线：
图6｜水面溅起的不平整度：峰值 ~4 mm，8 秒末收敛到 ~2 mm（接近网格单元高度，已是界面锐度下限）
图7｜液体区最大流速：从 0.83 m/s 衰减到 0.15 m/s，水基本静下来
图8｜t=7.9 s：水面基本平稳（α=0.5 等值线），颗粒静止堆在池底，全程无发散
7. 适用范围
这套自由面能力是通用的，两点口径说清：
二维圆盘 + unresolved 粗网格：颗粒不解析，拖曳是工程关联式，定量精度看场景，定性现象（入水减速、锥形堆、水面回落）可靠； 无表面张力：界面是「锐利但无张力」的层流自由面，不做毛细波和微尺度液膜；对本问题（入水堆积）影响很小。 网格分辨率：1 cm 单元 vs 5–8 mm 颗粒（粒径约 0.6–0.8 个单元），已到 unresolved CFD-DEM 的常规量级；颗粒更细时需同步加密网格。
换场景只需改生成器参数：池子尺寸、水深、颗粒粒径/数量、投放批次节奏。整条链（C++ 求解器扩展、两相算例模板、耦合脚本、验证门）已沉淀为可复用技能。验证门覆盖：静水面不流（空池 30 步 |U|<1e-4）、水体积守恒（全程 <10% 漂移）、堆积形态（沉底、集中、静止）、全程无发散——四个门全过，才算一个算例真正跑通。
