PFC 2D 圆盘强度上不去？arrlinear 实测：φ +13°，c +161 kPa

PFC 2D 圆盘强度上不去？arrlinear 实测：φ +13°、c +161 kPa

做 2D 双轴标定的同学大概都撞过这堵墙：颗粒摩擦系数 fric 从 0.3 调到 0.7，宏观摩擦角只从 18.8° 爬到 23.6°。调到头也就 24° 左右，怎么也上不去。

原因不在参数，在几何。2D 圆盘能自由滚动，颗粒之间"打滑"的路径根本不经过摩擦——你给再大的 μ，它滚过去就行了。

PFC 6.0 内置的 arrlinear（Adhesive Rolling Resistance Linear，粘性滚动阻力线性）正好有两个新旋钮：rr_fric 抑制滚动，adh_f0/adh_d0 提供法向拉力。我把它们跑了一遍：28 个双轴算例，共享同一份成样，两个假设都站住了。

1. 先搞清楚这个模型在干什么

arrlinear 的力一共三部分：线性接触力 + 滚动抵抗矩 + 粘附力。

粘附力的形式很简单：

F_a = F0                    当 gs <= 0     = F0 * (1 - gs/D0)      当 0 < gs < D0     = 0                     当 gs >= D0  其中 gs = gap - rgap（表面分离距离） 

四点必须记住的行为，否则容易做错：

• 它永不破坏。没有 bond 概念，isBonded() 恒为 false。所以 crack_num、DFN、sigma_ci 全部不存在，破坏只能用剪胀、配位数、组构来刻画。
• 它不提高抗剪。库仑判据 crit = Fn * fric 和滚动截断 M* = μr * R̄ * Fn 都只用线性法向力，粘附力不在里面。粘附只能靠"维持接触不分离"间接影响强度。
• adh_d0 会改变接触集。它同时决定 CMAT 的 proximity（接触检测范围），改完必须 model clean 重新检测。
• adh_d0=0 且 adh_f0>0 是禁忌。F0/D0 会打进临界时步的等效刚度里，D0=0 直接数值爆炸。

2. 实验怎么做的

28 个算例，2D 双轴，试样 0.06×0.12 m、7707 个球。所有算例共用同一份成样（在共享成样器里只成一次样，复 制给每个算例），这样初态逐位相同，跑出来的差异只能归因于参数。

扫描分六组：验证对照、rr_fric 单因素、fric × rr_fric 交互、adh_f0 单因素、adh_d0 单因素、围压包线（100/200/300 kPa）。

基准点选的是退化配置（rr_fric=0, adh_f0=0, adh_d0=0），此时 arrlinear 严格等于 linear，可以和以前的研究直接对照。

先说验证结果，它比结论本身更让我踏实：

3. 结果一：滚动阻力确实解锁了摩擦角

 

图1：rr_fric 从 0 加到 0.8，σ₃=200 kPa 下的表观摩擦角从 22.1° 一路涨到 35.1°。

三点围压拟合出来的莫尔-库仑包线更清楚：

更关键的是交互效应。同样把 fric 从 0.3 调到 0.7：

rr_fric = 0.0  ->  φ 从 19.4° 到 24.0°   斜率 11.5 °/单位 rr_fric = 0.4  ->  φ 从 26.7° 到 37.2°   斜率 26.3 °/单位 

滚动阻力不是简单地"加一个偏移"，它是把 fric 的杠杆放大了 2.29 倍。颗粒摩擦越高，滚动阻力的收益越大。

不过它也有饱和点：rr_fric 超过 0.4 以后，0.4→0.8 只再涨 2°。建议工作区间取 0.1~0.4。

4. 结果二：粘附给出 161 kPa 等效粘聚力

 

图2：adh_f0 从 0 加到 100 N，峰值偏应力严格线性增长，R²=0.9998。

这条线的意义在于它是线性且有截距的——强度包线不再过原点，这就是宏观"粘聚力"的来源。

 

图3：三族参数的强度包线。粘附族的截距 0.66 MPa 是实打实的，换算成 c = 161 kPa。

5. 结果三：不看峰值，看曲线形状

峰值只是曲线上的一个点。28 条完整曲线还能回答另一个问题：这两个旋钮是怎么改变破坏过程的。

我把每条曲线拆成 18 个形状指标——峰前的屈服点、峰值的应变位置、峰后的软化速率与残余比、接触丢失率。挑两条最重要的说。

滚动阻力让试样更"韧"：

峰值应变涨了 69%，脆性降 30%，接触网络保留得更完整。机制是：没有滚动阻力时，颗粒靠自转就能让试样"垮下去"，变形集中在少数滚动的颗粒上；抑制滚动后颗粒必须互相爬升，变形分布更均匀。

粘附不是"变脆"，是"崩得更陡"：

 

图4：注意弹性段（0~4%）所有曲线完全重合——粘附不改变初始刚度，变化全在峰值之后。

F0 从 0 加到 100 N，软化模量从 −13.8 陡到 −92.5 MPa（6.7 倍），但残余强度同时从 0.456 涨到 0.759 MPa，接触丢失率反而从 22% 降到 5%。

所以它不是"让材料变脆"，而是换了一种失效方式：粘附接触被克服时集体失效（陡降），剩下的摩擦骨架仍有较高残余承载力。想要"高强度 + 明显脆性跌落"的胶结岩行为，这个组合正合适。

真正变脆的是 adh_d0：取到 2 倍平均粒径时，残余比掉到 0.10、脆性 0.90——峰值 4.6 MPa 崩到 0.47 MPa。不过那组配置的接触集已经失真（配位数 13，而 2D 密排的物理上限约 6），只能作机制方向的证据。

一个附带的结论：残余强度几乎只由摩擦骨架决定。28 例里峰值跨度 0.12~4.6 MPa，残余跨度只有 0.19~1.10 MPa。要标定残余，调 fric 和 rr_fric；要标定峰值，用 adh_f0。

6. 那它到底适合模拟什么材料

跑完这 28 例，再回头看模型的机制，适用边界就清楚了。

它的定位是"有短程吸引力的颗粒 集 合体"，不是胶结固体。 官方手册对它的描述是 cohesive granular material——关键词是 granular（颗粒的），不是 bonded（胶结的）。而且写得很直白：它的吸引力与 PFC 的胶结材料不同，“没有破坏的概念”。所以它给得出"表观粘聚力"，给不出"胶结断裂"。

四个参数恰好对应四类可测的物理机制：

• adh_f0 / adh_d0 → 范德华力（官方明说是其线性近似）、液桥力（毛细）、静电力
• rr_fric → 颗粒棱角度（棱角抑制滚动）
• fric → 表面摩擦

适合的：

不适合的：

• 岩石、混凝土、强胶结砂岩——它们需要可断裂的键。这里吸引力永不破坏，给不出"键断→强度骤降"，也没有裂纹与声发射输出
• MPa 级粘聚力的材料——本次实测 F0=100 N 对应 c≈161 kPa，再往上只能靠加大 adh_d0，而 D0 超过一个粒径后，proximity 造出的虚假接触会主导（实测配位数 13，而 2D 密排的物理上限才 6）
• 干砂——用 linear 就行；要显式考虑棱角则用 rrlinear

一句话判据：粘聚力来自毛细/范德华/静电、且不需要看裂纹的颗粒材料，用它；靠水泥或胶结剂把颗粒焊在一起、想看损伤演化的，换 linearpbond。

7. 三个意料之外的发现

跑之前我以为 rr_fric 和 adh_f0 都是"强度旋钮"，结果不是。

一、adh_f0 是纯强度旋钮，rr_fric 不是。

adh_f0 从 0 到 100 N，峰值强度翻倍，但 E 从 17.13 到 16.92 MPa（几乎不动）、ν 从 0.125 到 0.106。而 rr_fric 从 0 到 0.4，E 涨了 24%、ν 直接减半。

对标定来说这是好消息：想单独调强度，用 adh_f0，不会把刚度和泊松比带偏。

二、adh_d0 是个被低估的强旋钮。

它的作用是设定 proximity——哪些远处球对算"接触"。对照组 A1（D0=0）和 A2（D0=5e-4）其他参数完全一样，只是接触检测范围不同：

A1  D0 = 0       q_peak = 0.241 MPa   配位数 4.11 A2  D0 = 5e-4    q_peak = 1.060 MPa   配位数 6.13    <- 4.4 倍 

机制不复杂：proximity 把分离距离小于 D0 的球对登记为接触，试样被围压压密之后，这些接触就变成真实承载的接触，等效于提高了骨架密度。

这个发现有个直接的副作用：评估"粘附力本身"的贡献时，不能拿纯线性算例当基准。正确做法是用同一 D0 下 F0=0 的算例。按这个口径，F0=100 N 的净贡献是 +54%（0.866 → 1.337 MPa），而不是拿 A1 比出来的 5.5 倍。

顺带一个坑：D0 = 1e-3 时实测配位数到了 13，而 2D 圆盘密排的物理上限约 6——接触集已经失真，那组数据不能用于物理结论。adh_d0 建议不超过一个平均粒径。

三、lin_mode 的影响比想象中大。

lin_mode 是法向力更新模式（0=absolute，1=incremental）。两个算例参数完全相同、只有它不同：q = 0.330 vs 1.108 MPa，差了 3.4 倍。

跨研究比较时一定要先核对这个设置。

8. 几个工程上的坑

并行不等于更快。 我一开始以为 4 路并行能省 75% 时间，实测总吞吐只有单进程的 1.11 倍（438 vs 393 cycles/s）。机器是 10 核 16 线程，瓶颈不在核数，更像 PFC 多进程之间的资源共享。想快，减算例比加并行有效。

超时值要按最慢的算例定。 我给无头运行器传了 2 小时兜底，而常规算例要 2.5 小时——结果每个算例都在 70% 处被杀掉，白跑了一轮。更麻烦的是 adh_d0 大的算例时步小、需要的步数是常规的两倍，兜底得设到 12 小时。

别用 extract_macro.py 切这种链的 .his。history export 会在最前面插一列 step，导出的列数比历史脚本假设的多一列，结果它把孔隙率当成了 crack_num，输出的"裂纹数"其实一直是孔隙率。这种错误不会报错，只会静默给错数。

剪胀角的公式和取值区间都得对。 平面应变下是 sin ψ = -D/(2-D)，D = dεv/dε1，而且只能用峰值之后的区间算——弹性压缩段的斜率不是剪胀。

9. 复现

bash# 28 例矩阵定义在 run_arrlinear_sweep.py 的 SWEEP 里
python code/run_arrlinear_sweep.py --matrix      # 看矩阵
python code/run_arrlinear_sweep.py --run --jobs 4 # 跑（约 12 h）
python code/collect_sweep.py                      # 出 CSV + 图 + 报告

完整的关系映射表（敏感度矩阵 + 闭式律 + 坑清单）在 results/mapping_table.md，模型说明书在 references/arrlinear-model.md。

10. 参考资料

• PFC 6.0 内置模型手册：Adhesive Rolling Resistance Linear Model
• 内置模型源码 contactmodeladhesive.cpp（粘附分段函数在第 944-950 行）
• 官方验证算例：verifications/adhesive_rolling_resistance/repose_angle
• 本研究全部产物：本构/arrlinear-双轴-参数研究/（含 28 例原始数据与复现脚本）