破坏带藏在哪？换个画法它就自己跳出来了
先说结论：位移云图看不出来的破坏带，换一张"颗粒重排强度"的云图就能一眼看到。三个 PFC 算例——直剪、单轴压缩、巴西劈裂——三种完全不同的破坏形态（水平剪切带、X 形交叉带、竖向劈裂线）全都是自己跳出来的，不用调任何参数。破坏带上颗粒的重排强度，是其余颗粒的 6 万到 17 万倍。
1. 先说结论
这篇讲两种后处理方式，针对的都是已经破坏了的试样：
第一种，看颗粒怎么动。 叫 D²min（非仿射变形）。它把"整体挪动"从位移里扣掉，只留真正的错位重排。三个算例的破坏带全都自己跳出来了，还顺带白送三个量：剪应变、体应变、转动。 第二种，看力怎么传。 叫粗粒化连续场。它把颗粒之间的接触力"抹匀"成一张连续的应力云图，可以直接和有限元/FLAC 的应力场对比。
两种方式回答的是两个不同的问题——“哪里坏了"和"力是怎么走的”。下面分成两大块讲，每一块都先讲原理，再看三个算例的效果。
2. 方法一 · D²min：看颗粒怎么动
2.1 2.1 先看一个问题：位移云图看不出破坏带
做离散元模拟，最常用的云图就是给每个颗粒涂上"它移动了多远"。
这张图有个天生的毛病：它分不清"整体挪动"和"内部错位"。
想象一块橡皮，你把它整体平移 10 厘米——每个点的位移都是 10 厘米，位移云图上一片通红，但橡皮内部一点没坏。反过来，如果橡皮内部撕开了一道口子，只要整体没怎么动，位移云图反而看着挺平静。
直剪那个算例最能说明问题。剪切盒是上下两半：上半个盒子被推着往前走，下半个固定不动。所以位移云图上，上半部分整体很亮、下半部分很暗，中间一条水平分界线。
看着像破坏带，其实那条线只是盒子的接缝——它是边界条件画出来的，不是材料坏出来的。
2.2 2.2 原理：把"整体挪动"扣掉
D²min 这个名字里，D 是位移 displacement，min 是"最小差距"。做法用大白话讲是这样：
挑一个颗粒，把它和周围一圈邻居看成一个小团伙； 先记下这个小团伙一开始的站位——谁在谁的正东边、谁在斜后方、间距多大； 变形之后再看一遍，这伙人现在站哪儿了； 然后问一句：假设这块料是"规规矩矩地均匀变形"，这伙人现在"应该"站在哪儿？ 实际站的位置，和"应该站的位置"差多少，就是中间这个颗粒的 D²min。
差得少，说明这块地方只是整体挪了挪、或者均匀地被拉伸压缩了——这属于连续介质能解释的运动。
差得多，说明邻居们真的重新排队了——你被挤走了，或者旁边那个滚开了。这才是不可逆的损伤。
这个算法有个很漂亮的特性：整体平移、整体转动、均匀拉伸，这三件事不管多大，D²min 都精确等于零（不是"很小"，是数值上的零）。所以图上剩下的全是"真的错位"。
更划算的是，算 D²min 的过程会顺带把另外三个量也白送出来：每个颗粒被剪了多少（剪应变 γ）、胀了还是缩了（体应变 ε_vol）、转了多大角度（转动 ω）。所以下面每个算例都给两张图——一张重排强度主图，一张四联图。
2.3 2.3 效果 · 直剪：一条水平带，正好落在上下盒的接缝上
图上那条橙红色的带，位置正好就是上下剪切盒的分界。
有意思的是：这条带不是"接缝"本身。接缝只是给了个位置，真正在这条线上发生的是——颗粒被剪得剧烈重排、胶结大量断裂。离这条带越远，颗粒越"安分"，颜色迅速回落到深蓝。
把四个量并排看，破坏带的样子更立体了：
剪应变 γ（右上）：带内中位 0.102，被剪了 10% 左右。剪应变可以是弹性的、能恢复；而 D²min 一定不可恢复。两个量在同一位置都亮，说明这里既有弹性剪切、也有塑性错动。 体应变 ε_vol（左下）：带内中位 +0.104，而带外几乎是 0（−0.0006）。正号代表膨胀——颗粒在带内被剪得"松开"了，体积变大。这就是剪切带的剪胀。 转动 ω（右下）：带内颗粒的转动量是带外的 53 倍。颗粒不只是在滑动，还在滚。
2.4 2.4 效果 · 单轴压缩：X 形交叉的两条带
单轴压缩的试样是 50×100 毫米的竖长条，上下两块刚性压板对着压。破坏之后，图上出现的是一个 X 形：两条斜带从试样中部往两端散开。
这是岩石力学里非常经典的破坏模式——共轭剪切带。两条带的夹角接近 90 度，各自大致沿着最大剪应力的方向。数值模拟里它是自己长出来的，事先没有在哪儿画过线。
这张四联图值得多看两眼：四个量给出的带几乎完全重合。
剪应变带内中位 0.167，被剪了 17%，比直剪更狠； 体应变带内中位 +0.234，膨胀 23%，也是三个算例里最剧烈的； 转动的带内外比是 13 倍。
这四个量物理意义完全不同——一个算相对位移、一个算应变张量、一个算体积变化、一个算刚体旋转。它们指向同一条带，这种交叉验证比任何单独一张图都可靠。
顺带说，注意顶部和底部附近颜色偏深。这不是说那里没受力，而是那里靠近压板、侧向约束强，颗粒没多少空间重新排队。
2.5 2.5 效果 · 巴西劈裂：一条竖直的劈裂线
巴西劈裂是把一个圆盘竖着对径压，看它怎么裂。结果非常干净：一条几乎笔直的细线，从圆盘正上方一直贯到正下方，其余整个圆盘基本都是深蓝色。
这条线就是圆盘的劈裂面，而且它走的正是加载直径——上极点和下极点连起来的那条线。
巴西这张的亮点在体应变（左下）：劈裂面上是一条清楚的正条纹，缝两侧几乎为零。
体应变是正的，意味着材料在这里是膨胀的、体积变大的——通俗说，这条缝是"张开"出来的。
这一点把它和剪切带区分开了。直剪和单轴的带内，剪应变和体应变同时很高（0.102/+0.104、0.167/+0.234），那是典型的剪胀：一边剪、一边胀。而巴西的缝上，剪应变只有 0.034，体应变却很明显——它是被拉开的。
剪切带是"剪着剪着松了"，劈裂缝是"被拉开的"。 光看云图说不出这个区别，靠这几个颗粒级的量可以定量说出来。
2.6 2.6 方法一 · 三个算例的定量对照
破坏带 vs 背景（带内取最强 5% 的颗粒）
重排强度：直剪 17.4 万倍 / 单轴 6.3 万倍 / 巴西 12.2 万倍
转动：直剪 53 倍 / 单轴 13 倍 / 巴西 16 倍
体应变：带内 +0.104 / +0.234 / +0.040，背景全部 ≈ 0
三个算例都差着好几个数量级。这说明"破坏"和"没破坏"在数据上本来就是两拨人，只是位移云图把它们混在了一起。
三种破坏形态、同一个算法、一次都没调过参数——这一点比单个结果本身更有说服力。
3. 方法二 · 粗粒化连续场：看力怎么传
3.1 3.1 原理：把接触力"抹匀"成一张连续场
颗粒之间其实只有接触力，并没有"应力"这个概念。应力是连续介质的东西。要把离散的接触力变成一张连续的应力云图，中间必须做一次"抹匀"——这个操作叫粗粒化。
听起来简单，抹法却有两种，结果差很多。
粗略的抹法：把每个接触力当成作用在接触点上的一个力，再在空间上滑动平均。这有两个毛病：一是算出来的场到底等不等于"应力"说不清，二是动量守恒不逐点成立——换个平均半径，场就变了。
正确的抹法：把每个接触力沿着两颗球的连线铺开，再把这些"铺开的力"叠加起来。
为什么必须沿连线铺开？因为一个接触力本来就是两颗球之间传的力，它真正"作用"的范围就是这两颗球之间的那一段，不是一个点。把它当点力涂上去，等于谎报了力的作用范围。
这么做出来的场有两个实实在在的好处：
动量守恒逐点成立：场里任何一小块，进出它的力都是平衡的。这意味着它能直接和有限元、FLAC 算出来的应力场逐点相减对比。 场是光滑可微的：能求梯度、能看应力路径。
顺带说一句 PFC 自带的做法——测量圆（在试样里放若干个圆、算圆内的平均应力）。它有三个问题：只在"圆内应力均匀"的假设下才成立；给出来的是离散的点值、还得自己插值；圆半径取多大没有客观依据。粗粒化没有这三个问题。
3.2 3.2 效果 · 单轴压缩：应力场和重排场对上了
看图 7 的竖向应力（右上）：深蓝是压得最狠的地方，中部有两条斜向的深色带交叉——位置和图 4 里 D²min 那两条橙红带是对得上的。
同一个破坏面，两种完全独立的方法（方法一算相对位移、方法二算接触力）指向同一个位置。这是两种方法之间的一次交叉验证。
再看剪应力（左下）：也是同一套斜带，说明破坏带就是剪力集中通过的通道。
注意试样形状：四个面板里试样的左右两边已经鼓出来了，上下的角上还缺了几块。这说明破坏之后试样已经碎成块、整体鼓胀，不是原来那个规规矩矩的长方体了。位移云图上看不出这个，因为它被整体变形淹没了。
3.3 3.3 效果 · 巴西劈裂：它为什么从中间劈开
这张是巴西圆盘破坏后的大主应力场。红色是拉，蓝色是压。
关键就一句话：圆盘里绝大部分区域是在受拉的，只有上下两个受压点附近是蓝的。
岩石怕拉不怕压（抗拉强度通常只有抗压强度的十分之一量级）。圆盘被竖着对径压，看着是在"压"，实际上内部被挤出了一个很大的拉应力区——所以它就从中间垂直劈开了，而且劈的方向正好是拉力最大的那条直径。
方法一告诉我们"它在哪里裂了"，方法二告诉我们"它为什么在那里裂"。
竖向应力 σyy（右上）：沿加载轴形成一条高度压缩的竖带，越靠近上下两个受压点压得越狠（端部约 −40 MPa），越往两侧越接近零。这就是经典的巴西圆盘应力分布。 剪应力 σxy（左下）：在上下两个受压点附近出现红蓝成对的花瓣状分布——载荷就是从这两个点灌进去的。 水平应力 σxx（左上）：大面积偏暖，又是一个"内部受拉"的证据。
3.4 3.4 效果 · 直剪：一个反向的例子
方法一在直剪上干净利落，方法二在直剪上却读不出那条水平剪切带。
四个分量上都是一片规则的对角条纹，走向和实际那条水平带完全不搭。
原因在边界条件：这个剪切盒的侧壁是刚性墙、而且有摩擦，加上顶部伺服加载过程本身有过冲，竖向荷载有很大一部分是通过侧壁摩擦传下去的。结果是应力场被"盒子"主导了，成了边界条件的指纹，材料自己坏成什么样反而盖住看不见了。
（顺带说明：这组图里的应力数值到了 ±40 MPa 量级，也是同一个原因——不是材料这么强，是加载方式把力灌进去太猛了。）
这个反向例子本身很有用：当应力场出现和破坏模式对不上的规则条纹时，它其实是在提醒你——这组边界条件正在干扰结果。 这时候应力云图就不再是结果，而是算例质量的自检工具。
3.5 3.5 方法二 · 量级验算：对上了教科书公式
云图好看不算数，得能对上解析解。
巴西圆盘有一个闭式解——一百多年前 Timoshenko 就推出来了。这个解说：沿加载直径，水平方向的应力恒为拉应力，数值等于 2P/(πDt)；在圆心处，竖向压应力等于这个拉应力的 3 倍。
拿算出来的场去对：
圆心处的压应力：解析解 −18.82 MPa，粗粒化 −18.86 MPa，差 0.2% 圆心处的拉应力：解析解 +6.27 MPa，粗粒化 +5.87 MPa，差 6.5% 沿竖直直径的剖面形状也对：拉应力基本恒定，压应力从中部往两端急剧集中
这一条才是方法二真正的底气。 前面那些守恒检验只能证明"算法自洽"，而对上解析解证明的是"我的数字是对的"。
（拉应力那 6.5% 的偏差主要来自取样窗口——圆盘中心的应力梯度很大，取一小块区域平均本身就有误差，不是算法的问题。）
4. 两种方法，分工不同
| 哪里坏了 | 为什么坏、力怎么走 | |
两个配合着用：方法一定位破坏，方法二解释破坏。
5. 用法与边界
用法很简单：任意两个时刻的颗粒位置 + 两个时刻的接触力，就能出图。破坏前、破坏后、任意中间状态都能算，四个颗粒级的量一次出齐。
三个提醒：
第一，方法二的接触力必须取"总力"。 平行粘结的接触，内力是"线性部分 + 胶结部分 + 阻尼部分"三项之和，而 PFC 里 contact.force.shear 给的是幅值（永远非负），直接用它会有一半接触的切向反掉。必须把三个分量分开导出再相加，否则算出来的应力会系统性偏低。
第二，方法二要配边界条件一起读。 像直剪那种刚性侧壁的算例，应力场会被边界主导、读不出材料行为；四边自由的单轴和巴西就没这个问题。方法一没有这个顾虑——它只看颗粒怎么动、不看力怎么传，三个算例上都干净。
第三，两套算法目前都是二维的。 三维版本公式一样（矩阵从 2×2 变成 3×3），但还没实测。
适用范围：三个算例都是 PFC 2D 的平行粘结胶结试样，二维平面应变、单位厚度口径。结论里的各种"倍数"是同批模型内部不同区域的相对比较，不是绝对值。巴西圆盘的解析解对照只对二维圆盘成立。
6. 参考资料
Falk, M. L. & Langer, J. S. (1998). Dynamics of viscoplastic deformation in amorphous solids. Physical Review E, 57(6), 7192–7205. —— D²min 的原始定义（方法一） Goldhirsch, I. (2010). Stress, stress asymmetry and couple stress: from discrete particles to continuous fields. Granular Matter, 12, 239–252. —— 粗粒化应力场的严格推导（方法二） Timoshenko, S. P. & Goodier, J. N. Theory of Elasticity. —— 圆盘对径受压的经典解
