PFC 的平行键模型，为什么画不出真实岩石的强度包络？
PFC 的平行键模型，为什么画不出真实岩石的强度包络？
如果你用 PFC 的平行键模型（parallel bond，pb）做过岩石三轴，大概率见过这个现象：
围压从 5 MPa 加到 20 MPa，峰值强度几乎是一条直线。线性回归 R² = 0.998，看起来漂亮得没什么可挑的。
但真实岩石不是这样。
围压越高，强度增长得越慢；到了高应力段，曲线会明显变平，甚至接近水平。岩石力学里把这叫"强度包络的趋平"，背后是微裂纹贯通、颗粒破碎、剪胀被抑制，最终进入临界状态。
pb 模型画不出这条曲线。不是参数没调好，是它的接触本构里根本没有这个机制。
这篇文章从 pb 的接触力学讲起，说明它为什么只能给出直线，再讲我怎么用 AI 把一条"会拐弯"的摩擦律写进接触模型，最后用 10 个围压、20 个双轴算例把包络画出来。
1. 先说结论
真实包络是三段式：低围压近线性，中围压凹向下，高围压趋于水平。Westerly 花岗岩的切线摩擦角从约 53° 降到约 35°。 pb 的直线来自接触级的两条线性判据：拉破坏和剪破坏。剪破坏是 τ = pb_coh + tan(pb_fa)·σ_n，斜率pb_fa是常数，宏观包络就只能是直线。新模型 nlpb：只改pbShearStrength()一个函数，把常数摩擦角换成随接触法向应力衰减的有效摩擦角，增加两个参数。默认参数下与内置模型逐位一致。验证：10 个围压（5–100 MPa）× 2 个模型。pb 的包络切线斜率从 7.5 降到 3.8；nlpb 从 3.1 降到 0.07——80–100 MPa 段，强度几乎不再增长。 局限：当前参数下 nlpb 的低围压强度也降了 28%；而且幂律衰减没有下限，高围压一路趋到接近零，比真实岩石"趋平但保留约 30° 残余摩擦"更激进。下一版要加摩擦下限。
一句话：pb 的直线不是标定问题，是接触级本构的结构问题。把常数摩擦角换成随法向应力衰减的函数，宏观包络就会自己弯下来——改动只有一个函数。
2. 一、我给了三个指令
整个过程我只提了三个要求。现在回看，每一条都在收窄问题的边界。
第一条：调研真实岩石的强度包络线，制定 pb 的优化方案。
我没有直接让它写代码。先搞清楚"真实岩石到底长什么样"，再谈改模型——否则很容易做出一个数学上自洽、物理上没依据的东西。
第二条：不考虑下弯，只考虑斜率降低，做一个新模型。
这是最关键的一次范围裁剪。真实岩石在高围压下确实可能出现强度下降，但那是有条件的：高孔隙岩石的孔隙塌陷、高温熔融、约束拉伸、真三轴卸载路径。常规三轴压缩下，完整岩石的峰值强度几乎总是单调增长，只是增长越来越慢。
如果一开始就把"下弯"当成目标，模型会被迫加一个缺乏普适物理依据的帽盖，标定会变得非常困难。砍掉它，问题就只剩一件事：让包络的斜率随围压递减。
第三条：用 pb 和 nlpb 做至少 10 个围压的双轴实验，看强度包线的形状。
这条要求把"模型能跑"变成了"模型有用"。两三个围压点看不出曲率，10 个点才能把"斜率如何变化"画出来。
3. 二、真实岩石的包络：三段式，不是直线
先把目标形状说清楚。
图 1 真实岩石强度包络的三种描述：灰色虚线是线性 Mohr–Coulomb，蓝色虚线是 Hoek–Brown，红色实线在 Hoek–Brown 基础上叠加临界状态平台；空心圆是 Westerly 花岗岩的实测点。
最直接的一组实验数据来自 Westerly 花岗岩。Blake 等（2016）的常规三轴数据，我按 Mohr–Coulomb 换算成切线摩擦角（sinφ = k/(k+2)，k 为偏应力对围压的斜率）：
Lockner（1998）汇总到 200 MPa 的拟合给出同样的趋势：0–50 MPa 区间等效 φ 约 50°，100–200 MPa 区间降到约 36°。
另外两组数据把边界条件也说清楚了：
Carrara 大理岩：30 MPa 以下脆性；30–300 MPa 半脆性；300 MPa 以上完全塑性，强度基本与围压无关，包络水平。 高孔隙砂岩（Bleurswiller，孔隙率约 25%）：约 118 MPa 平均应力处进入压实帽，临界屈服应力反而随围压下降（20 MPa 时 54 MPa，100 MPa 时 17 MPa）。这是"帽盖闭合"，不是简单的峰值下弯。
所以真实岩石的包络可以概括成三段：
低围压：近直线、斜率大； 中围压：凹向下，斜率递减； 高围压：趋平，进入临界状态。
"趋平"是规律，"下弯"是例外。 这就是第二条指令为什么合理。
4. 三、pb 的本构：为什么它只能画直线
4.1 3.1 一个接触由两套弹簧并联
PFC 的平行键模型，描述的是两个颗粒接触面上的一小块"胶结梁"。2D 里可以把两个圆盘之间的接触画成这样：
图 2 平行键接触示意：法向力 F_n、切向力 F_s、弯矩 M，以及两条破坏准则。
一个 pb 接触由两部分并联组成：
线性部分：法向弹簧 kn、切向弹簧ks，加一个 Coulomb 摩擦系数fric。它负责未胶结、或胶结断裂之后的接触；平行键部分：在接触面上传递法向力、切向力和弯矩，有自己的刚度 pb_kn、pb_ks和强度pb_ten、pb_coh、pb_fa。
胶结完好时，两部分一起受力；平行键一旦破坏，它承担的力和弯矩清零，线性部分的摩擦继续承担残余强度。
4.2 3.2 破坏准则只有两条
平行键的破坏判据，在源码里也是两行。
拉破坏：接触面上的最大拉应力（法向应力 + 弯矩贡献）达到 pb_ten：
剪破坏：接触面上的最大剪应力（切向应力 + 扭矩贡献）达到 Mohr–Coulomb 强度：
其中 σ_n = -F_n/A_b，压缩为正。
这两条式子的结构是完全一样的：一个常数项（pb_ten 或 pb_coh）加一个线性项（斜率 tan(pb_fa)）。pb_fa 是接触属性，是一个常数，不与任何应力、应变或历史量耦合。
这就是问题的根源。
4.3 3.3 宏观直线是怎么来的
单个接触的强度是线性的，宏观包络为什么也是线性的？
图 3 左：不同接触的剪强度都是同一条斜率的直线；右：宏观包络继承同一条斜率。
可以把宏观强度理解成接触强度的统计平均。设接触法向应力与宏观围压近似成正比：
那么每个接触的剪强度是：
对接触网络做平均后，宏观偏应力对围压的斜率仍然由 tan(pb_fa)·α 决定。只要 pb_fa 是常数、α 随围压变化不大，宏观包络就只能是一条直线。
仓库里七种岩性的双轴标定也验证了这一点：宏观摩擦角与微观键摩擦角近似线性：
pb_fa 覆盖 22.7°–73.2°、七种岩性，这条关系都成立。宏观摩擦角几乎是微观摩擦角的线性映射——微观常数，宏观也常数。
那为什么 pb 的包络还是会有一点弯曲？因为 α 不是严格常数。围压升高后，接触数目、接触方向和力链结构都会变化，组构演化会带来一点曲率。但这是二阶效应，量级不足以产生临界状态的水平段。
4.4 3.4 三处堵死了非线性
把原因列清楚，改造的方向也就清楚了：
τ = c + μσ 线性 | ||
pb_fa | ||
真实岩石的曲率来自三件事：剪胀被抑制、颗粒破碎、孔隙塌陷。第 1、2 条可以直接在接触本构里改；第 3 条需要可破碎颗粒模型（clump/子球），计算量大得多。
所以这次只改第 1、2 条：让摩擦系数随接触法向应力衰减。这正好对应"只做趋平、不做下弯"的范围。
5. 四、nlpb：把常数摩擦角换成应力相关的
5.1 4.1 本构设计
新模型 nlpb（nonlinear-envelope parallel bond），从内置 linearpbond 派生。
核心改动只有一处——把常数摩擦角换成随接触法向应力衰减的有效摩擦角：
两个新参数：
n：幂指数，取值0 < n ≤ 1。n = 1时完全关闭修正；s_ref：参考法向应力，控制"从哪个应力水平开始明显衰减"。
图 4 左：有效摩擦系数随无量纲法向应力衰减；右：对应的接触剪强度包络。n = 1 是内置的线性模型，n < 1 开始趋平。
这个形式不是拍脑袋。Guo 和 Morgan（2004）的 DEM 剪切模拟给出颗粒摩擦的压力依赖：μ ∝ σ_n^(-b)，b 约 0.21–0.40，对应 n = 1-b 约 0.6–0.8。接触级的幂律衰减，正好产生宏观的凹向下包络。
三处设计取舍：
只改强度，不改刚度。时间步由刚度决定，所以模型的稳定性与内置 pb 完全一致； 驱动量是接触自身的法向应力，不是宏观围压。它是一条真正的局部接触本构，不依赖外部施加的边界条件； 默认 n = 1时跳过全部修正，与内置linearpbond逐位一致。这一点后来成了最重要的验收门槛。
5.2 4.2 为什么用接触应力，而不是宏观围压
一个自然的疑问是：为什么不直接用 σ_3 当驱动量？
因为接触模型看不到宏观围压。它只知道当前这个接触上有多少法向力。如果非要把 σ_3 传进来，就得做全局耦合——每一步把围压写到每个接触上。这样模型就变成了"边界条件的函数"，换一个加载路径就失效了。
用接触法向应力，模型是局部的、路径无关的：同一个接触，无论在单轴、双轴还是剪切盒里，本构都一样。非均匀应力场（比如剪切带内部）也能自然处理。
代价是：接触应力与宏观围压的比例 α 需要标定。这个比例不是常数，会随组构和破坏模式变化。所以 s_ref 不能靠公式估，只能靠试算——这一点在后面的结果里会看到。
5.3 4.3 实现链路
图 5 实现链路：机械改名 → 注入本构 → 编译 → 部署 → 验证。
具体做了五步：
派生源码：写一个转换脚本，把内置 linearpbond的约 40 个标识符机械改名（类名、宏、模型名字符串）。改名后做残留检查——漏一个就是静默 bug，尤其是propagateStateInformation里用模型名字符串做比较的那处；注入本构：写一个补丁脚本，在 pbShearStrength()里替换摩擦律。每处补丁都断言"锚点命中次数 = 1"，避免改错位置；编译：2D 用 Release2D配置，产出contactmodelmechanical2dnlpb006.dll。这是本仓库第一次实机验证 2D UDM 配置；部署：把 dll 复 制到 PFC 的 plugins/contactmodelmechanical2D；验证：加载 dll、读回自定义属性、单接触扫压、双轴多围压。
新属性要注册到九个地方：属性枚举、属性名列表、成员变量、构造函数、archive()、getProperty、setProperty、copy、setNonForcePropsFrom。漏掉任何一个，都会出现"能编译、能跑、但结果不对"的情况。
两个最阴的坑，都不报错。
第一次双轴跑完，nlpb 与内置模型逐位相同——连峰值应变都一样。查下来：新属性漏进了 archive()，设好的参数在 save/restore 之后被静默重置成默认值，整条链跑成了对照组。
补上之后又发现第二个：contact cmat apply 在已有接触上切换模型时，不会携带 UDM 新增的属性。必须在 apply 之后显式再设一遍。
这两个坑的唯一发现方式，是做一次 save/restore 往返，把属性读回来打印。
6. 五、10 个围压，20 个算例
验证用了同一批花岗岩参数、同一份成键试样，只改变围压和接触模型：
围压：5、10、15、20、30、40、50、60、80、100 MPa； 每个围压跑内置 linearpbond和nlpb各一个；nlpb 取 n = 0.7、s_ref = 300 MPa。
图 6 10 个围压下的峰值强度包络。真正重要的不是两条线的高低，而是高围压段的走向：80→100 MPa 区间，linearpbond 还涨了 76 MPa，nlpb 只涨了 1.4 MPa。
关键数据：
| 292.5 |
看最后两行：
pb：围压从 80 加到 100 MPa，偏应力从 519.5 升到 595.4 MPa，还涨了 76 MPa； nlpb：同样区间，从 291.1 只升到 292.5 MPa，只涨了 1.4 MPa。
切线斜率更直观：
| 0.07 |
pb 的斜率从 7.5 降到 3.8，有下降，但始终保留明显的压力敏感性；nlpb 从 3.1 一路降到 0.07，80–100 MPa 段的包络已经基本水平。
图 7 包络的切线斜率 d(σ₁−σ₃)/dσ₃。linearpbond 在 2.8–4.5 之间波动，始终有压力敏感性；nlpb 单调衰减，到 80–100 MPa 段只剩 0.07。
7. 六、结果分析：趋平实现了，但还差一步标定
7.1 6.1 pb 自身也有一点弯曲
pb 的切线斜率并不是常数：从 10–15 MPa 的 7.5 降到 60–80 MPa 的 2.8。这部分弯曲不是来自接触本构，而是组构演化——围压升高后接触数目、方向和力链结构都会变化，宏观强度自然会偏离严格的直线。
但它不够。到了 80–100 MPa，pb 的斜率仍有 3.76，等效摩擦角约 40°，而真实岩石在同一应力段通常已经降到 35° 左右甚至更低。
7.2 6.2 nlpb 的斜率下降来自接触级
把切线斜率换算成等效摩擦角：
5–10 MPa：k = 3.13 → φ ≈ 37.6° 60–80 MPa：k = 0.89 → φ ≈ 17.9° 80–100 MPa：k = 0.07 → φ ≈ 1.9°
包络从"有摩擦"走到了"几乎无压力敏感性"，这正是临界状态的形态。
图 8 10/40 MPa 两档的应力–应变曲线与峰值点：同一围压下 nlpb 更早达到峰值，峰后更脆。
7.3 6.3 和真实岩石比，趋平过头了
把 Westerly 花岗岩的数据放在一起看：
两个结论：
pb 的方向是对的，但高围压段降不下去——它始终保留 36–41° 的摩擦，无法表达临界状态； nlpb 的方向也对，但降过了头——真实岩石趋平后仍保留约 30–35° 的残余摩擦，而当前 nlpb 一路趋到接近零。
原因在公式本身：幂律 (1 + σ/s_ref)^(n-1) 当应力趋于无穷时衰减到 0，没有下限。真实岩石的"趋平"是趋于一个非零的临界状态摩擦，不是趋于零。
7.4 6.4 下一版：加一个摩擦下限
合理的改法是：
μ_floor 对应临界状态摩擦（约 tan 30°–tan 35°），s_ref 控制衰减的快慢。这样包络在低围压保持原斜率，高围压趋平到一个非零平台，形状才和实验数据对得上。
7.5 6.5 当前的另一个问题：低围压强度整体下移
nlpb 在 5 MPa 处的偏应力是 141.7 MPa，比 pb 的 199.3 MPa 低 28%。这说明 s_ref = 300 MPa 偏小，衰减在低围压就开始起作用了。
下一轮标定的目标很明确：低围压段与 pb 基本重合，高围压段才明显趋平。s_ref 需要往 1 GPa 量级调，同时加上 μ_floor。这只能靠试算确定——接触应力的放大倍数比宏观围压高得多，单靠公式估算不够准。
8. 七、AI 做本构：能做的是什么，不能做的是什么
这次工作里，AI 完成的部分很具体：
查文献，把真实岩石包络的三段式形态和定量数据整理出来； 读内置模型源码，定位到唯一的线性判据； 设计并实现 UDM，编译、部署、做单接触验证； 发现并修复两个"不报错但结果不对"的静默坑； 生成 20 个双轴算例、批量运行、提取峰值、画包络、算切线斜率。
这些事的共同点是：知识密集，但不创新。难点在于"知道框架怎么用、坑在哪、数据怎么读"，而这些恰好可以固化进技能，交给 AI 稳定执行。
但有三件事必须由人来做：
划定物理边界。"不做下弯"这一条，直接决定模型能不能标定； 判断什么算证据。“至少 10 个围压"这个要求，把"模型能跑"变成了"形状对不对”； 看出结果过火了。nlpb 把高围压摩擦压到接近零，是数据摆出来之后才看到的——AI 可以算出这个数，但"真实岩石不会到零"这个判断来自岩石力学常识。
这也是我现在用 AI 做本构的基本方式：人负责物理判断和验收标准，AI 负责查、写、跑、排错。
9. 八、加入「AI × PFC 本构」交流群
如果你也在用 PFC 做岩石/土的本构、参数标定和二次开发，欢迎加入我们的交流群。
AI × PFC 本构开发交流群
加 QQ 763388012
备注「vibePFC」
加 QQ 763388012，支付 150 元获取完整技能包
包含：工程生成器 / 编译器封装 / 验证器 / 单接触试验模板 / 11 个内置模型源码 / FISH 调试清单 / 踩坑记录
仅支持 PFC 6.0.30，其余版本尚未适配
含安装与使用说明，以及后续更新。
群内会持续分享 PFC 6.0 的实战经验：运行、探查、命令字典、裂纹后处理、边坡、隧道，以及本构二次开发的完整流程。本文的 nlpb 模型、单接触试验模板和 10 围压双轴算例也会作为案例在群里交流。
10. 参考资料
Blake et al. (2016). The effect of fracture density and stress state on the static and dynamic bulk mod uli of Westerly granite. JGR Solid Earth. Lockner (1998)；Byerlee (1967)：Westerly 花岗岩高压强度数据。 Fredrich, Evans & Wong (1989). Micromechanics of the brittle to plastic transition in Carrara marble. JGR. Rybacki et al. (2021). Semi-Brittle Deformation of Carrara Marble. JGR Solid Earth. Tembe, Baud & Wong (2008). Stress conditions for the propagation of discrete compaction bands in porous sandstone. JGR Solid Earth. Guo & Morgan (2004). Influence of normal stress and grain shape on granular friction. JGR Solid Earth. Itasca PFC 6.0 用户手册：Linear Parallel Bond Model / Contact Model Plug-ins。
