从一句话到 110 m 工作面：AI 做 PFC 采空垮落模型的 4 版迭代

我只发了两句话。第一句给需求，第二句质疑结果。

帮我做一个完整的煤岩采空导致顶部岩层塌落的模拟，二维的。

最后拿到的东西是：140 m × 70 m 的 PFC2D 剖面、18 378 个颗粒、17 291 条裂纹、3.60 m 地表下沉（q≈0.60），以及一份 5.8 MB 的单文件报告。

中间 AI 推翻了 3 版模型。我只补了一句：

确定是有胶结的吧？感觉强度有点低。

先说结论：真正让模型垮掉的，不是把强度调到很低，而是把工作面从 64 m 放大到 110 m、覆岩从 37 m 加到 58 m。几何尺度够了，强度回到岩体尺度，基本顶仍然会破断。

1. 我发的两句话

1. “帮我做一个完整的煤岩采空导致顶部岩层塌落的模拟，二维的。” 这句话给了目标，没给参数。采高、埋深、工作面长度、覆岩层序、颗粒尺寸、边界条件、要出哪些指标，都得 AI 自己补。

2. “确定是有胶结的吧？感觉强度有点低。” 这是整个过程里最重要的一次人工干预。AI 当时交的模型能垮，但层内 pb_ten 只有 0.1–0.5 MPa，E 却给到 20 GPa。它去查了接触分配，确认层内确实全部胶结；然后承认强度是为了“让小模型垮”而压低的，决定重建更大尺度的模型。

2. AI 先做的不是写代码，是查家底

它先翻了本机已有的 PFC 技能：隧道开挖、边坡重力失稳、层理岩石、Voronoi 花岗岩、裂纹导出与监测、无头运行。

结论是：没有现成的“采空垮落”技能，但有四块可以复用——分层 CMAT、重力固结、裂纹追踪、结果探查。

于是它搭了一条新的脚本链：

0console.dat     初始化域与随机种子 1chengyang.dat   颗粒成样，顶部压板消除初始重叠 2gravity.dat     9 层分组 + CMAT 赋值 + 胶结 + 重力固结 3kaicai.dat      20 步开采，逐帧/逐裂纹输出 4extend.dat      采后继续松弛 50 000 步 postprocess.py   出图 + 单文件 HTML 报告 

 

图 1　初始剖面。9 层软硬互层，煤厚 6 m，覆岩 58 m，自由地表。

3. 前三个版本：为什么没成

第一版是 90 m × 45 m、4 m 采高、64 m 工作面、37 m 覆岩，层内直接给完整岩石强度（pb_ten 1.2–4 MPa）。

跑完只有 20 条裂纹，地表下沉 6 mm。

AI 先怀疑参数没赋上，写了个探针逐接触统计 CMAT：31 500 个球-球接触，层内强度分布和第一版的 7 个槽一一对上，只断了 21 个键。路由没问题。

真正的原因是物理尺度：二维圆盘颗粒在键断裂后会重新密排，碎胀被低估。4 m 采空区很快被冒落矸石填满，基本顶被托住，上覆岩层成拱，地表自然不下沉。

第二版把采高加到 6 m，第三版加到 8 m，同时把层间改成无胶结、层内强度降到 0.1–0.5 MPa。

8 m 版终于垮了：17 625 条裂纹、q=0.69、冒落带 31 m。

但参数组合变成了“高刚度 + 超低强度”：砂岩 E 20 GPa 配 pb_ten 0.5 MPa、泥岩 E 8 GPa 配 0.3 MPa，不像一块岩石。这就是我发第二句话的原因。

AI 的回应不是继续调低强度，而是把问题摊开：小模型 + 圆盘颗粒，必须靠超低强度才能补偿碎胀不足；如果要把强度放回自洽区间，就得把几何放大。

4. 第四版：把几何放大，让强度回到自洽

新模型 140 m × 70 m，6 m 煤，58 m 覆岩，110 m 工作面，两侧 15 m 煤柱，18 378 个颗粒。

层内 linearpbond，层间无胶结 linear。等效 UCS 约 10–30 MPa，E/强度比回到合理范围。

结果：

 

图 2　开采结束。覆岩破断、垮落，裂隙带贯通 58 m 覆岩。

 

图 3　地表下沉盆地。最大 3.60 m，q≈0.60，采空区上方形成宽缓下沉槽。

5. 模型具体怎么搭

1. 
   

2. 成样：ball distribute 在 140 m × 70 m 的盒子里生成 18 378 个颗粒，半径 0.30–0.45 m，目标孔隙率 0.16。顶部加一道压板，先把初始重叠压掉，避免颗粒把试样顶到盒外。

3. 分层与胶结：按颗粒 y 坐标分成 9 组，用 9 个 CMAT 槽逐层赋参数，判据是 range group 'Lx' matches 2（两端同层）。默认槽给层间无胶结 linear，摩擦 0.35；层内用 contact method bond gap 0.15 装平行键。装键前先清零接触力，让键在零应力状态下生成。

4. 重力固结：撤掉顶部压板，施加重力，跑 80 000 步到准平衡。把每个颗粒的坐标存进 ball.extra 作为位移零点，再重置位移。

5. 开采：每步把煤层组 L2 中 5.5 m 宽条带的颗粒标成 cut 并删除，跑 6 000 步让顶板响应。20 步采完 110 m。每步保存 .sav、写一行进度 CSV、存一帧颗粒位置。

6. 监测：fracture.monitor.p2fis 挂 bond_break 回调，每条断键记下拉/剪、位置、长度、开度。采完再追加 80 000 步，并从该状态续算 50 000 步；最后用独立探针导出逐裂纹 CSV 和键破损剖面。

层间为什么不胶结：沉积层理面近似无 cohesion。如果层间也胶结，9 层会组成一根厚梁，基本顶不容易破断，最后又得靠降低强度来凑。

6. 结果补图：应力、力链、配位数、裂纹方向

主图之外，AI 用探查技能从终态 3final_long.sav 又导了四组数据：接触 CSV、应力测量圆、配位数、裂纹走向。

 

图 4　主应力十字架。采空区上方平均应力降到 50–500 kPa（红/橙），两侧煤柱出现 1000–1800 kPa 的应力集中（蓝/绿），十字架在煤柱附近转为近竖直——典型的支承压力与应力拱。

 

图 5　强弱力链网络。深蓝是 fn 前 15% 的强链，集中在两侧煤柱和采空区下部；覆岩中上部以浅蓝弱链、灰色张拉链为主。承力骨架从煤柱向上收成拱形，拱内是垮落松散体。

 

图 6　完好键配位数场。intact 均值 2.59（全部接触为 3.77）。煤柱和完整岩层保持 4–6，采空区上方与地表附近降到 0–2；红色低配位区对应键断裂、接触丧失的损伤区。

 

图 7　裂纹走向玫瑰图（长度加权）。接近各向同性，仅南北向略强。说明覆岩是整体破碎，不是单一方向的剪切带；这也和“宽缓下沉盆地”的位移形态一致。

7. 动图：分层垮塌、位移与裂纹生发

 

图 8　推进过程。地表响应滞后于工作面，主要下沉发生在采后松弛阶段。

 

动图 1　颗粒分层垮塌。按 9 层岩性着色，能直接看到直接顶冒落、基本顶破断、上覆岩层弯曲下沉；橙色 区域是采空区范围。

 

动图 2　竖向位移演化。颜色是 0–4 m 的向下位移。随工作面推进，下沉槽从采空区向上、向地表发展；最后两帧是采后 130 000 步松弛。

 

动图 3　裂纹生发演化。颜色是出生时刻，紫早黄晚。早期裂纹从采空区顶板起裂，晚期向地表扩展；右侧累计曲线在 9–12 s 陡增，对应垮落雪崩。

8. 三个不显眼的坑

坑一：圆盘颗粒的碎胀被低估。 真实垮落矸石碎胀系数约 1.3–1.5，二维圆盘破断后能重新密排，碎胀接近 1。采高不够时，冒落带把采空区填满，覆岩不必下沉。这不是参数错，是模型尺度不够。

坑二：地表取粒阈值没有随模型高度缩放。 旧模板用 y > 40 m 取“地表”。换到 70 m 高的模型后，40–70 m 的深层颗粒也被算成地表，q 从 0.60 虚高到 0.92。修正为 y > 65 m 后，数字才对上。

坑三：冒落带判据不能取“最高的大于阈值条带”。 上部软弱层的局部高破损率会把冒落带顶界一路带到地表。改成“自煤层顶板向上连续、每 2 m 条带键破损率 > 25%”，31–40 m 的中砂岩硬层才成为屏障，冒落带落在 16 m。

另外几个 FISH 细节，都是跑长算例才会碰到的：

• 
  
• fracture.monitor.p2fis 不写 age 属性，直接读 fracture.prop(fr,'age') 会中止整个 dat；正解是在裂纹创建时写 fracture.extra(fnew,4)。
• FISH callback 随 .sav 保存。续算时重新注册会报 “already registered”，只要恢复计数全局量。
• FISH 的 else if 链只对应一个 endif，按普通 if 嵌套写会报 “No corresponding IF statement”。
• contact.prop(ct,'pb_state') 对无胶结的层间接触会报错，先用 contact.model(ct) 过滤。

9. 跑完不算完：独立探针复核

从 3final_long.sav 重新数了一遍：

• 
  
• 17 144 个颗粒，17 270 条拉裂纹 + 21 条剪裂纹 = 17 291 条，和 CSV 完全对上。
• 31 081 个层内键，断 8 868 个。
• 颗粒包围盒在模型内，没有颗粒飞出。
• CMAT 路由复查：same_layer 接触数 = linearpbond 接触数，层内没有漏胶结。

这里有一个必须说明的现象：每步只松弛 6 000 步，地表下沉滞后于工作面推进，主要下沉发生在采后 130 000 步。要采动过程中的动态下沉曲线，得把每步循环数加到 15 000–20 000，或者放慢推进。

10. 产物与复现

• 
  
• 报告：coal_mining_caving_large/postprocess/report.html（5.8 MB，单文件）
• 图件：初始剖面、终态覆岩、位移场、下沉盆地、推进演化、覆岩分区、动图
• 数据：mining_progress.csv、final_cracks.csv、bond_profile.csv、frame_*.csv
• 探查：3final_long_contacts.csv、3final_long_stress.csv、3final_long_ballfield.csv
• 脚本：gen_case.py、0console → 1chengyang → 2gravity → 3kaicai、4extend.dat、postprocess.py
• 本案例已收录进 VibePFC 社区。

复现命令：

bashcd coal_mining_caving_large
python ~/.claude/skills/pfc-run/code/run_pfc.py main.p2dat --timeout 7200
python ~/.claude/skills/pfc-run/code/run_pfc.py 4extend.dat --timeout 3600
python postprocess.py

11. 人和 AI 的分工

这次分工很清楚：

• 
  
• AI 负责查技能、写脚本、跑算例、写探针、出图、打包报告。
• 我负责给目标、看结果、提出“这个强度不像岩石”的判断、决定重建。

AI 能把链路跑通，但“这个数像不像岩石”这种判断，还是得人来看。

诚实声明：这是二维平面应变模型，参数未针对具体矿区标定。两带高度、下沉系数 q 的绝对值需要现场资料校准；圆盘颗粒对碎胀的低估、每步 6 000 步的快速推进，都会影响动态过程。本文给的是机理与趋势，不是某个矿的定量预测。

作者介绍

 

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

用到的技能：pfc-run、pfc-bedding、pfc-slope、pfc-crack-monitor、pfc-result-probe、pfc-contact-probe、pfc-contact-forcechain、pfc-contact-coord、pfc-measure-stress、pfc-crack-rose、pfc-crack-evolution。