首页/文章/ 详情

用 AI 技能做 LBM-DEM 耦合：从圆柱绕流验证到水力劈裂

lobby

50分钟前浏览0

PFC 6.0 里有一个 CFD 模块，但它只做数据交换，不含流场求解器。官方那条 Darcy 算例的渗透率，是从 Kozeny–Carman 经验式算出来的，不是解出来的。

想让颗粒在水里被冲走、被冲塌，或者被水压把缝撑开，就得自己接一个流体求解器进来——把格子 Boltzmann（LBM）和离散元（DEM）双向耦合起来。

这篇文章讲两件事：

方法：LBM 和 DEM 怎么接起来、怎么证明它是对的、接起来能做出什么；
载体：这套东西不是一篇教程，而是一个技能包——一段 AI 能直接调用、自己跑、自己验证的能力。

1. 耦的是什么：两套离散，一个回路

DEM 那边算颗粒：接触力、胶结、断裂。LBM 那边算流体：速度场、压力场。两边各自都不新鲜。

耦合的难点是它们必须互相知道对方：流体要知道颗粒在哪、多大，否则算出来的流场是空的；颗粒要知道流体施加了多大的力，否则颗粒永远不会被水推走。

图1 上半是这套耦合的物理图像。流体格子比颗粒细得多（Δx ≤ 0.025 d），所以每一个颗粒、每一条颗粒之间的间隙都被真正解析出来，而不是被当成一团孔隙率。流体绕过去，同时对颗粒表面施加一个力——这个力不是估出来的，是由动量交换直接算出来的。

图1 下半是闭合回路的数据通路：颗粒几何送进流体侧 → 流体推进、算出每颗粒的受力 → 力回传给颗粒 → 颗粒运动、可能断裂 → 下一轮。

颗粒位置决定流场，流场反过来决定颗粒怎么动——这就是"双向耦合"的含义。

1.1 1.1 先讲清楚方法的边界

两个硬约束决定了它能用在什么地方。

一是格子速度上限。 LBM 是显式时间推进的，格子速度必须远小于声速（否则不可压假设失效）；同时格子黏度有下限（否则 BGK 碰撞本身不稳定）。再叠上"格子必须细到能解析颗粒"这一条，物理流速的上限就定了：

颗粒直径	格子黏度	格子雷诺数	物理流速上限
1 mm	0.05（常用）	2	0.08 m/s
1 mm	0.005（BGK 极限）	20	0.80 m/s
2 mm	0.005	20	0.40 m/s

超过约 0.5 m/s，或者颗粒雷诺数 Re_d 超过约 100，这套方法就不适用了——那是压力基 NS 求解器的地盘。

二是网格量。 格子要细到 Δx ≤ 0.025 d。一个 30 mm 见方的二维试样、1.5 mm 颗粒，就是 800 × 800 = 64 万格点，这已经是纯 NumPy 能跑的上限；如果按真实的 0.1 mm 级颗粒算，是 576 万格点，单步 1.15 秒。

所以这套方法的地盘很清楚：几十毫米量级、流速在米每秒以下的颗粒–流体相互作用。出了这个范围，该换方法，而不是调参数。

2. 为什么它被做成一个「技能」

代码谁都能写。难的是让它可信，而且下次还能用。

这套 LBM-DEM 耦合没有停在"跑通了"，而是把代码、方法、边界、踩过的坑一起打包成了一个技能包（pfc-lbm-coupling）：

部分	内容
`SKILL.md`	契约：这套方法什么时候适用、硬约束、验证义务
`code/`	可运行的求解器与验证脚本：D2Q9、D3Q19、socket 服务端、三个 `validate_*`
`references/`	方法文档：socket 协议、流速上限推导、水力劈裂的方法学边界
`pitfalls.md`	45 条实机坑，每条带"错了会看到什么症状"

对做这件事的 AI 来说，技能包改变的是四件事。

一、不用凭记忆写命令。 PFC 6.0 的命令语法和 5.0 差得很远，而且写错往往不报错——比如 command 块里的局部变量会被原样送进解析器、生成一个名字奇怪的对象，日志里一声不响。所以"写每条命令前先查离线文档核验语法"被固定进了流程。

二、"验证"是硬性条款，不是提醒。 技能里写着：改了任何求解器代码，就跑对应的 validate_*。这条规则在本次是实打实生效的——下面验证那一节里，动量交换测力被单独验到了 0.000%，因为它就是整套耦合里所有拖曳力的来源。三维那条 validate_duct.py 还抓出过一个把渗透率压低 11.2% 的切片错误，而那个错误在流场云图上完全看不出来。

三、坑带着症状一起存。 45 条坑每条都写明"错了会看到什么"。举个真实的例子：出口不闭合质量收支时，阻力系数会线性漂移、永不收敛，但升力系数和回流长度全都正常。没有这句症状描述，这种坑事后根本无从查起；有症状，定位就从"到处怀疑"变成"对号入座"。

四、边界和失败也一起交接。 技能里明确写着哪些事没做到：裂隙内的流动没有被解析、裂纹压力通道还是分段的、冲刷算例里有一处已知的换算系数错误导致位移绝对值不可用。把这些写下来，比多做一个算例更有价值。

这个技能包不是为这篇文章新写的——圆柱绕流基准、水力劈裂耦合、以及中间所有的后处理配图，都是在这个框架下做出来的。下面按「先验证、再看应用、最后深入一个算例」的顺序讲。

3. 先把地基钉死：四级验证

自写求解器最大的风险，是「算出来像模像样，其实是错的」。

所以我们用了一条逐级加复杂度的验证链：能被解析解钉死的先钉死，再往上加物理。

3.1 ① 泊肃叶管流：求解器本体

压力驱动的槽道流动有解析抛物线解。

峰值误差 0.086%，而且网格加密后收敛阶 2.00——正是理论上的二阶精度。这一关说明：离散格式本身没问题。

3.2 ② 库埃特剪切 + 测力

上下壁面一静一动，流场是解析线性剖面。

这一关的真正目的不是看速度场（那个误差是 6.2e-14，机器精度），而是验证动量交换测力：

动量交换测出来的壁面剪应力，与解析值误差 0.000%。

这一步是整套耦合的地基。LBM–DEM 耦合里，颗粒受到的拖曳力全部来自动量交换——如果测力不对，后面所有的力都不可信。它精确到 0.000% 之后，DEM 侧才敢直接用这个力，不必再怀疑符号和量级。

3.3 ③ 方腔驱动流：加进回流

上壁面拖动，方腔内出现主涡和角涡。这一关有实验基准（Ghia 等，1982）。

内部速度剖面误差 < 0.5%。图上标注的 max err 0.106 是归一化速度的绝对偏差（顶盖速度取 1），它出现在紧挨顶盖的那一层——顶盖边界层没有被分够。加密即可改善，模型本身没有问题。

3.4 ④ 圆柱绕流：加进曲边界与涡脱落

这是最接近工程问题的一关，用的是 DFG 圆柱基准（Schäfer & Turek, 1996）：

通道 2.2 × 0.41 m，圆柱直径 D = 0.1 m，圆心在 (0.2, 0.2)；
入口抛物线速度，最大 0.3 m/s，平均 0.2 m/s；
圆柱故意偏心（下方 0.2 m，上方 0.21 m）——这个不对称是 Re = 100 时触发涡脱落的关键。

两个工况：Re = 20 定常（尾流是一对对称回流涡）和 Re = 100 周期（卡门涡街，升力做正弦振荡）。

分辨率取了 D/h = 20（圆柱直径占 20 个格子）。这个数字不是调出来的，是纯 NumPy 的算力逼出来的：D = 20 时每步 14–20 ms，D = 40 要 70 ms，单算例 35–60 分钟，不实用。

定量结果

量	本文模型	Schäfer–Turek (1996)	偏差
Re=20 阻力系数 C_D	5.793	5.58	+3.8%
Re=20 回流长度 L_a	0.0900 m	0.0847 m	+6.3%
Re=100 C_D 范围	3.310 – 3.610	3.10–3.16 / 3.22–3.24	约 +12%
Re=100 C_L 振幅	1.278	2.00	−36%
Re=100 Strouhal 数	0.3065	0.30	+2.2%

四个数字里有两个很准，两个偏得明显——这组偏差本身就说明了方法的性质。

准的是频率：Strouhal 数 St = fD/U 只差 2.2%。涡脱落的频率对分辨率不敏感，只要尾流不稳定，脱落的节奏就是对的。

偏的是曲边界的受力。原因在下面。

圆柱是「阶梯」的

半程反弹（half-way bounce-back）只认识与网格对齐的墙面。所以圆柱不是光滑圆，而是由一格格方块拼出来的阶梯圆（图5 里能看到锯齿）。

阶梯圆的周长比真圆长——周长长，阻力就大，所以 C_D 系统性偏高 3.8%～12%。这是一阶精度的典型代价。

而 C_L 振幅偏低 36%，是同一件事的另一面：阶梯边界把圆柱表面的涡量生成抹平了一部分，涡脱落的强度被低估。频率还在，力度不够。

这不是 bug，是明确的取舍：阶梯近似实现简单、与半程反弹天然兼容；要用插值反弹（Bouzidi）做到半精度，需要为每个边界节点维护到壁面的距离，是另一个量级的工程量。当前精度对「流动形态」和「频率」类问题够用；如果目标是颗粒受力精度或者自由沉降，就该换方法。

顺带一个有用的自检

Re = 20 是定常解，物理上应该上下对称。我们在中线上量了一下横向速度的残余：max |u_y| = 1.6e-3（格子单位）。相比主流 0.1 的量级，这个残余很小——说明流场确实收敛到了对称的定常解，没有残留的不对称扰动。

这种**「用物理对称性当自检」**的办法很便宜，也很有用：算例本身就告诉了你答案的一半。

4. 接起来能做出什么：三个实例

验证完之后，双向耦合真的用起来了。三个算例都在 LBM/ 下：

水下沙坡冲刷——100 颗颗粒堆成三角坡，逐级抬高入口流速，观察坡脚起动：

位移在入口流速 0.074 → 0.082（格子单位）之间由降转升，起动阈值确实存在。机理也清楚：驱动力正比于 u²（平方增长），而静摩擦有上限（常数），两条线的交点就是阈值。

不过这里要说清楚一个边界：这套方法的流速上限（第 1.1 节）恰好就压在这个阈值附近，而且这个算例早期有一处换算系数错误——时间步长 c_dt 被填成了正确值的 1/5，力因此被放大了 25 倍。所以这个算例的位移绝对值不可用，只有趋势和起动拐点可用。

边坡入水——坡体颗粒在水下失稳、沿坡面下滑：

水下颗粒坍塌——柱状颗粒床在水下坍塌、堆积：

这三张图里，背景色是流场速度，颗粒位置是 DEM 实时传过去的。流场不是事后叠上去的——颗粒每动一步，流场就跟着重算。

5. 水力劈裂：把一个算例做完整

前面讲的是怎么接、怎么验证、能做出什么。这一节把一个算例做到完整——从模型、假设、参数，一直到结果和它没做到的地方。

5.1 5.1 模型

试样是 30 × 30 mm 的方形胶结颗粒体，中心挖一个 φ10 mm 的圆孔作为注入孔。

颗粒直径 1.35–1.65 mm（4 档级配，Cu = 1.2），共 379 颗；
用平行粘结（parallel bond）连接；
孔内以恒定流量注水，水压升高到超过胶结强度时起裂。

孔为什么取 φ10 mm：它是 6.7 个粒径。最初想用 6 mm，但在 d = 1.5 mm 下那只有 4 个颗粒宽——注入边界会被颗粒尺度的不规则性主导，而不是被孔的形状主导。所以放大了。

5.2 5.2 五个假设，逐条说清楚

水力劈裂跟冲刷不是一回事。冲刷是「水推着颗粒走」，劈裂是「水把缝撑开」。这个差别决定了整套方法都要换。

假设一：过程是准静态的（蠕流）

依据：算一下两个时间尺度。

LBM 的声速 c = c_s Δx / Δt。本例 Δx = 37.5 μm、Δt = 7.03 μs，得 c ≈ 0.31 m/s。压力波扫过 30 mm 试样需要 约 1400 个 LBM 步。
而试样尺度的黏性扩散时间 L²/ν ≈ 900 秒，孔隙尺度 d²/ν ≈ 2.25 秒。

孔压的建立远慢于颗粒的力学响应。按真实时间同步推进既没有必要，也会让计算量失控。

做法：改成准静态增量——每注入一小份水量，流体求解器求出该时刻的稳态孔压场，把流体力施加到颗粒上，让 DEM 跑到平衡，再进入下一步。

代价：不模拟真实秒数。进度量改用「累计注入体积」，而不是物理时间。

假设二：裂隙内的压力 = 孔内压力（压力平衡假设）

依据：先看看裂隙里的流动能不能直接算。

裂隙开度约 0.1 mm = 2.7 个格子，而该处的压力梯度约 3×10⁷ Pa/m。按立方定律估流速：v ≈ 25 m/s，对应的格子速度约 4.7——超出 LBM 稳定上限约 100 倍。

在这个网格分辨率下，裂隙内的流动根本无法解析。

物理上，裂隙的导流能力远高于基体，压力沿裂隙接近均匀。所以改成：裂隙内部压力等于孔内压力，用裂纹位置处的「湿化区」在流体侧表达。

代价：这是近似，且有明确的适用范围——见 5.8。

假设三：基体几乎不渗漏

依据：在一个增量内，黏性扩散连一个孔隙都过不去（孔隙尺度扩散时间 2.25 s，而可模拟的物理时间约 0.5 s）。

验证：这一点被结果直接证实了——孔压随注入量严格线性上升（图12 上半）。如果水从基体渗漏流失，压力曲线会出现平台或回落。线性的压力曲线，就是「水全部用于撑开、没有流走」的证据。

假设四：二维

孔隙尺度与颗粒尺寸相当，二维圆盘堆积的渗透率比真实岩石高出约 7 个数量级（5.7 节详述它的后果）。

假设五：强度按同样比例缩尺

因为孔压天花板被压到了 kPa 级，胶结强度也必须跟着缩。这不是「凑数」，而是为了让「孔压超过强度才起裂」这个相对关系仍然成立。

5.3 5.3 做法：一次交换走一个循环

每一轮交换做四件事（左侧 DEM、右侧 LBM）：

DEM 侧收集所有颗粒几何（id、x、y、r），通过 socket 发给流体侧；
流体侧把颗粒体素化成固体掩码，重新确定注入源，然后推进 112 个 LBM 步——每步往源节点加一点密度；
流体侧用动量交换测出每颗粒受到的拖曳力，回传；
DEM 侧把力施加到颗粒上，松弛到平衡（ratio-average 1e-5），其间可能断键产生新裂纹；裂纹位置单独走文件通道传回流体侧。

然后进入下一轮。这次运行走了 120 轮，墙钟 54 分钟。

为什么裂纹要单独走文件通道？因为键断了不等于颗粒分开了——两个颗粒几何上仍然贴着，间隙约等于 0，流体侧的固体掩码照样把那个位置判为固体。纯几何的连通域分析看不见裂隙。

所以裂纹位置必须由 DEM 侧显式传出去，流体侧再据此湿化。

5.4 5.4 参数：三个表

DEM 侧（PFC）

参数	取值	说明
试样尺寸	30 × 30 mm	2D 方形
颗粒级配	1.35–1.65 mm，4 档，Cu = 1.2	379 颗
目标孔隙率	0.18	分层压缩成样，按级配精确计数
颗粒模量 E	2×10⁸ Pa
颗粒摩擦系数	0.5
密度	2600 kg/m³
胶结模量 pb_emod	2×10⁸ Pa	不变 ——刚度决定应力分布，不决定强度
胶结黏聚力 pb_coh	6 kPa	由 10 kPa 缩尺而来
拉压比 ten_coh	0.5	不变 ——它决定破坏模式，不能缩
胶结抗拉强度 pb_ten	3 kPa	= 6 kPa × 0.5
胶结摩擦角 pb_fa	40°
注入孔半径	5 mm	用临时墙环成孔，加胶结后拆除

一个容易被忽略的点：缩尺时只缩强度，不缩刚度、不缩拉压比。刚度管的是应力分布，拉压比管的是破坏模式（拉还是剪）——这两个一改，模型就不是原来那个物理问题了。

LBM 侧

参数	取值	说明
网格	800 × 800	Δx = 37.5 μm
格子黏度 ν_lb	0.005	BGK 稳定下限，它反过来决定了物理压力的天花板
物理黏度 ν	1.0×10⁻⁶ m²/s	水
时间步 Δt	7.03125 μs	= ν_lbΔx²/ν
注入流量 q	1.0×10⁻³ m²/s	每米厚度的体积流量
每轮 LBM 步数	112	热启动，不是冷启动求稳态
裂纹湿化半径 r_cr	25 格	模型参数，不是物理开度
交换轮数	120

换算系数（这三个必须对）

系数	公式	取值
长度	c_scale = 1/Δx	26666.67 格/m
时间	c_dt = ν_lbΔx²/ν	7.03125×10⁻⁶ s
力	c_fscale = ρΔx³/Δt²	1.0667 N/m（每格子单位力）
压力	p_scale = ρΔx²/Δt²	2.844×10⁴ Pa

力的换算是二维专用。2D 用的是 ρΔx³/Δt²（单位是「每米厚度的力」，与 PFC 二维的力的约定一致），3D 要多一个 Δx。这两个式子各自都对，只差一个维度——用 3D 的 Stokes 律去检验 2D 的式子，会得到一个「差一个 Δx」的假结论。

r_cr 是模型参数，不是几何量。 它表示「每条裂纹携带多宽的压力走廊」，直接决定压力场的形态。参考量级：25 格约等于一个粒径（d = 1.5 mm = 40 格）。

它必须全链同值——耦合运行用什么值，事后重建云图就必须用什么值。否则两张图根本不是同一个流场。

5.5 5.5 结果

孔压线性上升，破裂阶梯式发生

两条曲线的形态对比很有意思：压力几乎完全线性上升，而破裂是阶梯式的。

裂纹数停滞一段时间后突然增加几条，再停滞，再增加。每一次台阶对应一次局部破裂事件的级联，而不是连续开裂。

完整过程分三段：

阶段	注入交换	孔内超压	裂纹数
升压	0 – 16	0 → 4.1 kPa	0
起裂与停滞	16 – 30	4.1 → 7.7 kPa	1 → 6
阶梯式扩展	30 – 120	7.7 → 27.8 kPa	6 → 35

起裂点出现在第 16 次交换，孔压 4.1 kPa。

起裂之后有一段明显的停滞期：第 1 条裂纹出现后，十来个注入轮次（一直走到第 27 轮）都没有继续扩展，而孔压仍在稳定上升——直到第 30 轮、7.7 kPa 才等到第一次成批破裂（一轮之内从 2 条跳到 6 条）。这种「停滞—爆发」交替是脆性材料破裂的常见形态。

破坏与裂隙扩展过程

从六格演化里能看到几件事：

起裂位置不对称：第一条裂纹不是均匀地从孔周冒出来，而是先从孔壁某一侧开始。这是因为离散元模型里的孔壁不是光滑圆弧，而是几十颗颗粒拼成的粗糙边界，局部应力集中的位置由颗粒的随机排列决定；
裂纹沿径向向外辐射：符合孔内压力作用下孔壁受拉的预期；
最终扩展到接近试样边界：从孔壁的 5.6 mm 扩展到约 14.4 mm（试样半宽 15 mm）。

终态

孔周存在一定宽度的高压带，压力有沿裂隙向外延伸的趋势。

5.6 5.6 力都去哪了：一次必要的核算

胶结强度从 5 kPa 缩到 3 kPa 之后，有一个问题必须回答：孔压真的能把键拉断吗？

我们把最终几何（第 120 轮结束时的颗粒排布）连同同样的累计注入量重新装进流体求解器，单独跑到稳态，统计 LBM 实际回传给孔壁颗粒的力：

（两者的累计注入量相同，差别在收敛程度：耦合运行每轮只推进 112 步热启动，重建跑了 3000 步。所以这里的 30.98 kPa 比耦合终值 27.8 kPa 略高，更接近该几何下的稳态孔压。这不是矛盾，是两个收敛程度不同的读数。）

量	数值
孔内超压（重建稳态）	30.98 kPa
解析孔壁向外力 p·2πR	973.4 N/m
LBM 实际回传的径向合力	244.8 N/m
解析值 / 实测值	3.98
孔壁环单颗颗粒的径向力（取量值平均，共 81 颗）	6.47 N/m
单个键的抗拉能力 σ_t·d（3 kPa × 1.5 mm）	4.50 N/m
单颗径向力 / 键承载力	1.44

两条结论：

方向正确，量级只有解析值的 1/4。原因是孔壁粗糙（颗粒拼成的边界不是光滑圆弧），加上颗粒之间存在渗漏。换几何、换分辨率，这个诊断必须重跑——它不是一个可以外推的常数。
键确实断得掉：孔壁单颗颗粒受到的径向力平均 6.47 N/m，而一个键的抗拉能力是 4.50 N/m——比值 1.44，力超过了承载力。这与「孔压 4.1 kPa 就开始起裂」是自洽的。

这个诊断还有一个副产品：它解释了为什么缩尺的时候，缩的是强度而不是随便取个数。如果 pb_ten 仍按原计划的 5 kPa（承载力 7.5 N/m），比值就变成 0.86——键根本拉不断，模型只会一直升压。强度必须缩到「力能超过承载力」的位置，现象才出得来。

5.7 5.7 为什么所有数字都是 kPa

这是这条路线最重要的一条边界，必须说透。

链条是这样的：

二维圆盘堆积的渗透率由几何决定，k ≈ d²ε³/(180(1-ε)²)。代入 d = 1.5 mm、ε = 0.19，得 k ≈ 1.4×10⁻¹⁰ m²；而真实岩石约 10⁻¹⁷ m²——差了 7 个数量级。

在孔壁维持 MPa 级压力梯度，Darcy 定律要求流速 u = k∇ p/μ 达到数十 m/s——远远超过 LBM 的稳定上限。而 LBM 的黏度下限（τ > 0.5）决定了格子速度上限，这个上限反推回来，就是孔压天花板。本例实测约 30 kPa。

所以 pb_ten 必须跟着缩（5 kPa → 3 kPa），否则 2 MPa 的键用 kPa 级的孔压永远拉不断，模型只会一直升压。

现象是完整的：孔压积累 → 超过强度 → 起裂 → 裂隙扩展 → 压力场跟着破裂网络演化。只有绝对量级缩小了约 3 个数量级。
所以：报 kPa，不报 MPa。

5.8 5.8 还没做到的两件事

第一，裂纹湿化是「斑块」而不是「连续通道」。

现在用的是圆盘湿化（半径 r_cr），而裂隙是一条细线。圆盘的面积大部分落在固体颗粒里，被过滤掉之后只剩零散节点——压力场呈现分段斑块，不是连续走廊。

要连续，需要沿裂纹线段湿化：DEM 侧多传裂纹长度和走向，流体侧把「点 + 半径」改成线段扫描。这是下一步。

第二，裂隙内的流动没有被解析。

按 5.2 的假设二，裂隙压力直接取孔内压力。这在裂隙导流能力远高于基体时成立，但如果要研究裂隙内的流量分配或者滤失，就必须提高分辨率或引入专门的裂隙流模型。

另外要明确一点：这个流体求解器的定位是「准静态压力场求解器」，不是「模拟真实注水过程」。它回答的是「注入量—孔压—破裂」之间的关系，不是真实的时间尺度。

6. 这套技能怎么用

技能名 pfc-lbm-coupling，代码是纯 NumPy，无编译依赖、无 GPU 依赖。

bash# 任何求解器改动之后，自检必须重跑（几秒钟）
python validate_duct.py            # 3D，对 Boussinesq 解析解
python validate_poiseuille2d.py    # 2D，对泊肃叶解析解
python validate_couette2d.py       # 2D，★ 动量交换测力，必须报 0.000%

bash# 双向耦合：两个进程，先起服务端，再跑 PFC 驱动
python lbm_server.py --nx 400 --ny 140 --port 3334
python <pfc-run>/run_pfc.py example_two_way_slope.dat

技能里还有一条单向用法（不把力送回颗粒，只测量渗透率），代码是三维 D3Q19。它严格说不是耦合，但共用同一个求解器内核，写在这里备查：

bashpython <pfc-run>/run_pfc.py export_ball_geometry.dat --3d
python run_permeability.py --balls balls.csv --nx 128 --ny 64 --nz 64 \
       --nu 0.15 --g 1e-3 2e-3 --dump field.npz

验证义务不是可选项：

改了任何求解器代码，就跑对应的 validate_*。三维的 validate_duct.py 曾经抓到一个把渗透率压低 11.2% 的切片错误；
二维的 validate_couette2d.py 必须报 0.000%——它验的是动量交换测力，也就是双向耦合里所有拖曳力的来源。它一旦不是 0.000%，耦合的力就整体不可信；
单向那条用法必须做孔隙率交叉验证（体素化孔隙率 vs DEM 侧孔隙率），并且同时报重叠量——只有孔隙率数字没有鉴别力。

7. 小结

LBM-DEM 耦合解决的是"颗粒和水互相改变对方"这一类问题。地盘在几十毫米量级、流速米每秒以下；出了这个范围，该换方法而不是调参数。
验证要逐级做。解析解先钉死求解器，实验基准再加物理。动量交换测力那一关（0.000%）是整套耦合的地基。
圆柱绕流给出了一个诚实的精度画像：频率量（Strouhal）准到 2%，受力量的偏差来自阶梯边界的一阶精度，C_L 振幅低估 36%。知道哪一类量可信，比一个笼统的「误差 5%」有用得多。
水力劈裂的完整现象链已经能复现：孔压线性建立 → 超过强度起裂 → 裂隙阶梯式扩展 → 压力场同步演化。绝对量级缩尺约 3 个数量级，报 kPa。
把它做成技能，比跑通一个算例更有价值。代码、"什么时候不适用"的边界、45 条带症状的坑，一起沉淀下来——下一次做新算例，不用重新发现它们。
还有两件事没做完：裂纹压力通道的连续化，以及裂隙内流动的直接解析。

诚实声明

① 求解器是自写的（D2Q9 / D3Q19，纯 NumPy），不是 Itasca 官方实现，也没有对标任何商业 CFD 软件。它的正确性由第三节的四级验证背书，边界也在那里写明了。

② 圆柱算例用的是阶梯边界，曲边界只有一阶精度：C_D 系统性偏高 3.8%–12%，C_L 振幅低估 36%。频率量（Strouhal）准到 2%。知道哪一类量可信，比一个笼统的"误差 5%"有用。

③ 水力劈裂那一路的绝对量级是缩尺的。二维圆盘堆积的渗透率比真实岩石高约 7 个数量级，孔压天花板约 30 kPa——所以强度和结果都以 kPa 报，不以 MPa 报。复现的是现象链，不是真实岩石的绝对数值。

④ 裂隙内的流动没有被解析，用的是「裂隙压力 = 孔内压力」的压力平衡近似；裂纹压力通道目前是分段斑块，尚未连续化。

⑤ 冲刷算例里有一处已知的换算系数错误，那个算例的位移绝对值不可用，只有趋势与起动阈值拐点可用。这处错误连同修正都记在 LBM/README.md 里，没有藏起来。

⑥ 本文只做了二维。三维那条路由 PFC 的 CFD 模块 + 外部求解器完成（见另一个技能），与本文不是同一套东西。

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

8. 参考资料

格子 Boltzmann 方法：D2Q9 / D3Q19，BGK 碰撞 + 半程反弹（half-way bounce-back）
动量交换测力：Ladd (1994)；Aidun & Clausen (2010) 综述
方腔基准：Ghia, Ghia & Shin (1982), J. Comput. Phys.
圆柱基准：Schäfer & Turek (1996)，DFG 流动基准 2D-1 / 2D-2
渗透率参考：Kozeny–Carman 关系（也是 PFC 官方 Darcy 算例所用）
技能包：~/.claude/skills/pfc-lbm-coupling/（求解器、socket 协议、验证脚本、45 条实机坑）
求解器完整验证表：LBM/README.md；建设过程与决策记录：LBM/HANDOFF.md
水力劈裂算例留档：水力劈裂-lbm/（脚本链 + stage0–4 分阶段产物 + 后处理报告）

来源：超级大的lobby

断裂碰撞 UG python UM 离散元裂纹理论 PFC 材料

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2026-10-09

最近编辑：50分钟前

lobby

硕士 | 无擅长颗粒流PFC

获赞 963粉丝 5927文章 92课程 23

作者推荐

可试听

【PFC6.0】三维Cluster三轴模拟碎石离散元

¥39

还没有评论