不用色标：给 PFC 试样印一层网格，剪切带自己画出来
先说结论
位移云图在均匀变形段几乎一片同色，局部化带被整体 位移"淹没"。换成实验力学里的老办法——加载前给颗粒印上规则图案，加载后图案本身就是材料线：条纹被拉歪的地方就是剪切，方格被压扁的地方就是压缩，不需要色标。
我把它做成了 PFC 6.0 的 pfc-dye-pattern 技能：五种图案、一次加载算例、网格法反算局部应变，并用独立管线交叉验证。胶结岩石单轴、砂土双轴、三维单轴、砂土边坡都能直接用。
1. 一、位移云图的问题，不是颜色不好看
用 |u| 或 vx 给颗粒连续着色，是 PFC 后处理里最常见的做法。问题出在"连续色标"本身：
试样整体被压缩 4%，所有球的位移都叠加了这一个均匀分量，色标随之被拉满； 真正关心的剪切带只占试样的一小块，颜色差异被均匀分量压掉； 换一个色标范围，图看起来就变了，结论跟着"变色"。
实验力学早就遇到过这个问题，解法是网格法 / 条纹法：在试件表面印规则网格，变形后网格线的畸变直接就是应变场。它的好处是自参照——图案自己和自己比，不需要外部色标，也不受刚体 位移影响。
这套东西可以原样搬到离散元里。区别只在于：实验里印在表面，PFC 里印在颗粒标签上。
2. 二、让图案成为材料属性：三个机制
第一，标签在加载前就定死。
成样平衡之后、加胶结之前，遍历所有球写四个量：
之后不再改动。实测跨 model save/restore、加胶结、删侧墙、整个加载过程，标签漂移都是 0。
第二，初始坐标必须自己存。
PFC 没有 ball.disp 这种"相对初始位置"的 getter。把 x0、y0 写进 extra(1)/extra(2)，一次探针导出就同时拿到参考构型和当前构型，不用导两帧再逐球相减。
第三，图案是初始坐标的纯函数。
floor((x0·cosθ + y0·sinθ)/s) mod 2 这类定义只依赖初始坐标。所以一次加载算例跑完，想换图案、换间距、换角度，离线重算标签即可。文章里这五张图案图，没有多跑任何一次 PFC。
3. 三、五种图案，一次算例全出
图1 五种图案都用同一批颗粒的初始坐标重算，不用为每种图案重跑加载。
单元尺寸不是随便给的：先用目标尺寸，再取整到能整除试样尺寸，保证单元严格铺满、面积相等。否则边缘会留下半格，应变拟合在边界失真。
4. 四、胶结岩石单轴：材料线直接画出破坏
先看最经典的一个：50×100 mm、6447 个球的平行胶结试样，真单轴压缩（删掉侧墙）。
先看同一套棋盘着色在加载前后的变化：
图2 左为变形前，右为变形后。格子边界始终连续，试样整体被压矮；但破坏集中在哪里，连续色块还看不出来——这正是需要把网格线单独抽出来的原因。
图3 灰虚线是初始网格，橙/蓝实线是同一批颗粒在峰后（ε=0.759 %）的位置，位移放大 5 倍。
这张图能读出三件事：
试样的破坏不是均匀的。网格线在试样中部偏右形成一条斜向的错动带，带外的格子基本保持正交； 靠近上下压板的两端变形明显更小——那层球被摩擦咬住、又顶着弹簧墙，是典型的端部死区； 带的位置、宽度和倾角都能直接量，不需要先做位移阈值分割。
定量结果：带与加载轴夹角 148°，独立跑一遍 pfc-shear-band 是 149°，相差 1°，显著性 z=59.5；网格法轴向应变 0.781%，压板实际施加 0.759%，相对差 2.86%。这一例的快照在峰后 30%（峰值 51.29 MPa，快照时 15.39 MPa），所以看到的是已经破坏的试样。
5. 五、砂土双轴：剪胀与两个独立口径
无胶结砂土不能做单轴——没有侧向约束、没有胶结，试样一加载就散。所以砂土用双轴：围压 σ3=0.2 MPa，加载到轴向应变 4.5%。
图4 左为变形前，右为变形后。棋盘格整体被压矮、拉宽，格子边界仍然连续——说明标签跟着材料走，没有"换色"。
把方格边界单独抽出来，变形集中区就清楚了：
图5 中部明显向外鼓，靠近上下压板的两端横向变形小；网格线在试样中部被拉斜，方向与最终的剪切带一致。
图上看到格子歪了，能不能给出数？可以。把每个染色单元内颗粒的形心当作测点，对相邻单元做最小二乘，得到局部变形梯度 F：
一个 5 mm 单元里有几十个球在平均，噪声被压掉；而且与图一一对应：图上歪得最厉害的格子，就是 γmax 最大的格子。
图6 单元通道：20×10 个染色单元，每格一个 γmax 值。
图7 逐球通道：每颗球一个局部 γmax，分辨率高得多，噪声也更大。两张图给出的高应变区走向一致。
砂土这一例最有说服力的是两个完全独立的量测口径：
网格法测到的侧向应变 +6.712%，伺服墙独立测得 +6.735%，相对差 0.33%。前者来自颗粒标签，后者来自边界位置，互证了剪胀； 剪切带与加载轴夹角：单元通道 22°、逐球通道 24°、独立 pfc-shear-band24°，最大差 2°，z=28.2；轴向应变：网格法 4.479%，压板施加 4.500%，相对差 0.46%。
6. 六、边坡：滑面不用猜
同一条脚本链，把加载换成重力、把"胶结"换成"开挖"，就得到砂土边坡算例。染色在开挖之前就印在完整方块上；开挖只删球、不移动球，所以留下来的图案仍然精确铺满坡体。
图8 上为变形前，下为变形后。坡面轮廓已经改变，坡脚出现明显的堆积体。
把材料网格抽出来，滑面的位置就直接显形了：
图9 坡面 50°，重力驱动跑到碎屑静止。滑体整体沿坡面下滑，材料网格被拉成弧形，坡脚出现堆积——滑面位置不需要靠速度阈值去猜，图案本身已经把破坏区画出来了。
这一例最终有 64.3% 的颗粒越过移动阈值，最大位移 89.8 mm，相当于坡高的 1.5 倍，约 58 万步后按静止判据停止。
边坡加载阶段单独用了 damp_load=0.2：成样阶段的 0.7 阻尼是为准静态压实准备的，直接沿用到重力驱动会过度耗散动能，坡体只发生毫米级蠕变。把两个阶段的阻尼拆开，同一个模型才真正滑起来。
7. 七、定量口径与验证
这套方法不是"画得好看"，每一个通道都有可核对的口径：
网格单元通道：每个染色单元一个应变计，拟合局部变形梯度，适合报数； 逐球通道：复用 pfc-shear-band的位移梯度 + Radon 找带，分辨率高，适合找带的位置；侧向材料量规：轴向用首末两行单元形心，侧向用首末两列。不能用颗粒外轮廓——自由棱角的球会外翻，量到的是边缘剥落，不是材料应变； 独立交叉验证： pfc-shear-band用自己的入口、自己的读写列再跑一遍，抓"列选错、轴搞反、平面搞错"这类接线错误。
几个关键门控的实测：
三维单轴算例也在其中：棋盘图案沿 x/y/z 三个方向都正确翻色（违例 0/0/0）。它的试样偏矮胖，局部化很弱，带倾角一项按规范记"不适用"，而不是硬凑一个通过。
8. 八、怎么用
技能已经固化成一条链。生成算例、跑五个阶段、导出前后两帧、出图、算应变、交叉验证：
后处理：
最需要先想清楚的是三个参数：
pattern | dstripe，看二维畸变用 checker | |
cell | ||
test | ucs biaxial / slope | ucs；砂土必须 biaxial；边坡用 slope |
关于作者
超级大的 Lobby｜同济大学岩土硕士 · 10 年 PFC · 离散元技术讲师
做 PFC 入门与进阶课程，也在做 AI + PFC 的自动化分析技能库。这个号写数值模拟的算例复盘、可复用的工作流，以及建模过程中真正会遇到的问题。
诚实声明
① 砂土双轴这一例的快照停在轴向应变 4.5%，约为峰值的 87%，不是峰后状态——无胶结砂土没有明显脆性跌落，应变到点即停。
② 三维算例试样偏矮胖，局部化不显著，带倾角一项按规范记 n/a。
③ 网格法在试样两端有"死区"：贴着压板的球被摩擦咬住，变形小于内部，容差按 5% + 3d/L 计入这个物理效应。
9. 参考资料
技能目录： C:\Users\76338\.claude\skills\pfc-dye-pattern\逐球通道来源： pfc-shear-band（位移梯度拟合 + Radon 找带）图像后端：matplotlib（PFC 无头控制台不提供 plot）
