AI 从零做一个水合物双轴算例：从孔喉成核到 380 根拉伸断裂
先说结果：一个 20×40 mm 的二维水合物沉积物试样，从建模、跑算例到微观分析，全部由 AI 完成。峰值强度提高 38%（1.674 MPa，砂骨架 1.213 MPa）；加载中 998 根水合物–砂桥断了 380 根，全部是拉伸破坏、零剪切；破坏在空间上高度分散，没有形成剪切带。
这篇文章讲两件事：AI 是怎么把这个模型一步步搭起来的，以及结果里的每张图、每个数字到底说明了什么。
1. 这个案例在算什么
水合物在沉积物里有三种赋存形态：孔隙填充、颗粒包裹、骨架支撑。孔隙填充型的特征是——水合物以独立小颗粒充填在孔隙里，只在孔喉处和相邻砂粒连接，而砂粒骨架本身没有被胶结。
这个区别在离散元里很关键。很多水合物模型图省事，直接给砂–砂接触加胶结键，那是"胶结型"，力学上完全是另一回事。
所以这个算例要回答三个问题：
孔喉成核的几何构造，能不能真正复现孔隙填充形态？ 围压下，水合物桥对峰值强度和残余强度的贡献一样吗？ 这些桥到底是怎么断的、断在哪儿？
2. AI 是怎么搭起来的
2.1 水合物不是"键"，是真球
第一步就定了基调：不新增"胶结键"来代替水合物，而是真的生成独立的水合物 ball。接触模型按"谁和谁接触"分开：
砂–砂： linear，不胶结水合物–砂： linearpbond水合物–水合物： linearpbond球–墙： linear
2.2 孔喉成核：让每颗水合物都成为"桥"
第一版用 ball generate 在孔隙里随机放球，结果很难看：它只保证不重叠，水合物全悬在孔中央，和砂粒之间留着空隙。膨胀 0.5% 之后，100 个球只有 3 个形成了水合物–砂接触——剩下的都是孤立漂浮颗粒，不能传力。
正解是在砂粒对的孔喉处解析成核。对每一对邻近砂粒（不要求已经接触），求两个圆的交点：
delta 是 2 µm 的受控重叠，保证接触能被检测到。每个交点天然同时接触两颗砂，也就是天然保证 2 根桥。候选点洗牌后逐个放入，直到达到目标饱和度。
饱和度按二维面积定义：Sh = A_h / A_pore，目标 10%，实际放入 496 颗，回算 9.46%。
2.3 固结之后再成键
这是双轴和单轴最大的不同。顺序必须是：砂骨架先生成 → 伺服固结到 σ3 = 0.5 MPa → 在固结状态下放水合物、成键 → 保持围压加载。
如果反过来先成键再加围压，键会在固结过程中被压坏，后面测到的就不是"围压下形成的桥"。
3. AI 踩的坑（这部分比结果更有用）
整个过程里 AI 犯了至少七个错，每个都影响结果：
坑 1：2D 孔喉放不下 0.12 mm 的水合物球。 砂粒 r≈0.425 mm、孔隙率 0.16 时，孔喉内切圆半径只有 0.06–0.07 mm。球太大，几乎处处和第三颗砂重叠。最后把 rh 降到 0.09 mm 才填得进去。3D 文献里的 rh/r_sand ≈ 0.3 不能直接搬到 2D。
坑 2：ball generate 放出来的球没有接触。 见上面，100 个球只有 3 个接触。
坑 3：固结后成键，围压瞬间掉到 0。 这是最隐蔽的一个。在固结态给接触赋模型时如果执行 contact cmat apply，它会重新生成接触，把已经建立的围压接触力全部清零。正解是只对 h–s / h–h 接触用 contact model linearpbond range ... 逐类赋值，不碰砂–砂的 linear 接触。
坑 4：断键还原时多出 53 条"假接触"。 统计断键要对比加载前后的接触表。但 final 的接触表里混进了 53 条加载中新生成、从未成键的 h–s 接触（state = 0）。不剔除就会高估断裂数。
坑 5：配位数画错了。 第一版把全体 1671 颗球的配位数一起统计，均值 1.19——被 1175 颗砂粒稀释了。只统计水合物才是 2.01，而且初始 490/496 颗恰好 2 根键、6 颗 3 根。这个分布本身就是孔喉成核的几何签名。
坑 6：粒径画成了固定像素。 早期图件用 scatter 画球，点的大小和真实半径无关，只是示意。后来改成按 CSV 里的 r 画真实半径的圆——水合物在总览图里立刻变成小点，这才是真实比例。
坑 7：剪应变峰值落在试样角上。 网格化统计 γmax 时，贴近墙面的格点球数少、受边界约束，拟合出的高应变是端点效应。统计窗口内缩 3 mm 后，峰值才回到试样内部。
4. 结果：+38% 的峰值，和一条完整的断裂链
含合物试样的峰值是 1.674 MPa，砂骨架 1.213 MPa，提高 38%。但峰值出现的位置从 3.94% 应变提前到了 1.05%。
10% 应变时两者的残余强度分别是 0.990 和 0.923 MPa，只差 7%。
这条曲线的信息很直接：水合物提高的是峰值和脆性，残余强度还是砂骨架的摩擦咬合说了算。
5. 微观机制：为什么全部是拉伸
这是整个案例里最值得看的部分。
加载结束后，998 根水合物–砂键断了 380 根。用裂纹数据里的 ball1/ball2 和断键球对逐根比对，380 条裂纹对 380 根断键，一条不差。
而且这 380 根全部是拉伸破坏，剪切破坏为 0。
断键方向的各向异性 A = 0.38，主方向接近水平。原因不难想：竖向压缩让试样侧向膨胀（宽度从 20 mm 鼓到约 23 mm），水平排列的两颗砂粒彼此远离，夹在孔喉里的水合物桥就被直接拉断。平行键的拉伸强度和剪切强度都设成 5 MPa，但受力以法向张拉为主，所以剪切破坏一次都没发生。
配位数的变化同样直接：
初始平均 2.01 根，加载后剩 1.25 根。按颗粒算：21.6% 的水合物失去全部键，33.5% 至少断掉一半，45.0% 的键完整保留。
最后是"破坏是不是局部化"这个问题。视觉上，裂纹铺得很散；定量上：
裂纹最近邻指数（Clark–Evans）R = 0.93，完全随机是 1.0，只偏离 7%； 最大剪应变场的峰度高 2.06，正态分布是 3——比正态还平坦。
两条都指向同一个结论：这是弥散的张拉破坏，不是剪切带局部化。
6. 这套流程能复用什么
回头看，这个案例真正可复用的不是参数，是几个习惯：
每个数字都要能从 CSV 复算。 断键数、配位数、组构、位移场，全部从加载前后的两个快照算出来，脚本留在目录里。报告里写的 380 / 2.01 / 0.38，都能一次跑出来。
交叉验证比单看一张图可靠。 裂纹数和断键数一致、裂纹位置和断键位置互相印证——两个独立通道给出同一个数，才敢写进结论。
图要和数据分离。 先出 micro_summary.json 和几个中间 CSV，再画图。改配色、改布局时不动数据，改数据时不用重画。
诚实声明：这是二维圆盘模型，水合物粒径只有砂粒的 0.2 倍；只做了 Sh ≈ 9.5%、σ3 = 0.5 MPa 一个工况；键强度和刚度没有按真实水合物的力学性质标定，结果应看趋势而不是当定量预测。另外加载后期有 6 颗颗粒穿透墙体（位移达数十 mm），已在微观统计中剔除。
写在最后
这个算例最有意思的地方，不是那 38%，而是几何完全决定了力学：水合物球放在孔喉，就注定每颗有 2 根桥；桥的方向沿砂–砂连线，就注定竖向压缩时被侧向膨胀拉断；断得又散又均匀，就注定看不到剪切带。
AI 在这中间做的是把这条链子一环一环验出来——包括它自己一开始画错的粒径、算错的配位数。
VibePFC 社区 · 一起把 AI 玩进 PFC
我成立了一个 VibePFC 社群，专门聊 AI 怎么落地到 PFC / 数值模拟——让大模型读懂模型、自己改脚本、自动排查"看不懂"的问题。收一个 50 元的门槛费。
想进群的同学，加 QQ 763388012，备注"vibePFC"即可。
算例包获取 · 集赞 30
本算例支持集赞获取：
① 将本文分享到朋友圈；
② 集满 30 个赞；
③ 加 QQ 763388012（备注"水合物算例"），凭截图获取本算例完整包：建模脚本 + 加载脚本 + 后处理脚本 + 结果数据。
参考资料
Brugada J., Cheng Y.P., Soga K., Santamarina J.C. (2010). Discrete element modelling of geomechanical behaviour of methane hydrate soils with pore-filling hydrate distribution. Granular Matter, 12: 517–525. Jung J.W., Santamarina J.C., Soga K. (2012). Stress–strain response of hydrate-bearing sands: Numerical study using DEM. J. Geophys. Res., 117: B04202. Waite W.F. et al. (2009). Physical properties of hydrate-bearing sediments. Rev. Geophys., 47: RG4003. Itasca Consulting Group (2023). PFC 2D 6.0 Documentation.
