零厚度裂隙的两种建模方法，结果差了 10 倍
零厚度裂隙的两种建模方法，结果差了 10 倍
同样一条 45° 零厚度裂隙，PFC 里两种主流建模方法给出的强度折减差了将近 10 倍：一个掉 23%，一个只掉 2.5%。差别不在参数，在接触法向。
1. 问题的由来
用 PFC 做含裂隙岩石的双轴模拟，"零厚度裂隙"差不多是绕不开的东西。文献里常见两套做法：
方法一：用 fracture create 在试样里铺一条 DFN 裂隙，再 fracture contact-model 把与它相交的接触换成 smoothjoint 模型。不引入新几何，只改接触本构，所以裂隙是真正零厚度的。
方法二：成样阶段插一道 45° 的刚性 wall 把颗粒分到两侧，加完胶结后把墙删掉。删墙瞬间两侧颗粒恰好相切，界面就是一层无胶结的 ball-ball 接触（linear 模型）。
两套说法都自称"零厚度"，参数也都能配得差不多。问题是：它们算出来的东西一样吗？
我拿同一个试样（60 × 120 mm 平行胶结，围压 10 MPa）做了对比，另外跑了一个无裂隙的完整试样做基线。三个算例共用同一份成样——同种子、同撒球，7948 颗球的位置逐字节相同。
2. 结果一：宏观响应差了一个量级
完整试样的峰值偏应力 123.8 MPa。方法一 95.0 MPa，折减 23.3%。方法二 126.9 MPa，折减 2.5%——落在数值离散的范围内，几乎等于没影响。
方法二的曲线和完整试样基本重合，从初始刚度到峰后脆性断裂都对得上。
方法一则是另一副样子：0 到 0.25% 应变区间有一段零应力平台，之后斜率明显更缓（弹性模量 6.07 GPa，只有完整试样的 36%），峰值也推后到 1.23% 应变，峰后是渐进下降而不是脆断。
3. 结果二：破坏形态是两种东西
完整试样是标准的 X 型共轭剪切带，两条带对称发展。
方法一的微裂纹明显沿着 45° 裂隙面集中，拉成一条连续的宏观破裂面——蓝色虚线就是预制裂隙的位置，红黑色的裂纹基本贴着它走。
方法二的裂纹分布和完整试样几乎分不出来，还是 X 型共轭带，裂隙面上没有出现额外损伤。
定量一点看：把每条裂纹中心投影到垂直于裂隙面的方向，统计落在 ±2 mm 带内的比例。完整试样 3.9%（相当于随机水平），方法一 19.0%，方法二 7.4%。
所以方法二并不是完全没影响，只是影响比方法一小一个量级。
4. 结果三：差在接触的"力学法向"
这两套方法的单接触刚度是一样的量级（都是 1e7 N/m），界面接触数量也相当（105 个 vs 123 个）。既然这样，为什么宏观表现差 5 倍以上？
区别在接触的力学法向。
smoothjoint 模型的特点是把接触法向强制锁定成裂隙面方向。我实测了加载前那 123 个 smoothjoint 接触，它们的 dip 属性全部等于 135°，一个例外都没有——用 Oda 组构张量算各向异性度 A = 2.000，这是理论极限值，意味着界面在数学意义上是一个完美定向的面。
方法二的界面法向由删墙后两侧球的球心连线决定。虽然主方向也在 135° 附近，但有明显散布，A = 0.745。
这个差别在加载中会放大。smoothjoint 的法向不随变形改变，接触永远只能提供沿面法向的支撑，切向完全靠摩擦，所以界面可以持续滑移。方法二的接触法向跟着球心走，界面呈锯齿状互相咬合，宏观上滑不动。
5. 结果四：界面在加载中的演化完全不同
把每个接触沿其法向画成一条短线段、按接触模型着色，能看到更直观的东西。
上排是刚胶结完的状态。两种方法都只在 45° 方向上留下一条很细的异色接触线，几何上高度相似——这也是为什么光看初始模型，会觉得两者差不多。
下排是破坏后的状态，两边彻底分道扬镳：
方法一的弱面随破坏自动扩展。smoothjoint 接触从 123 个涨到 946 个，而且方向散开了（各向异性度掉到 0.539），沿裂隙面两侧铺成一条有宽度的破碎带。原因是裂纹追踪脚本把新生的每条微裂纹也写进 DFN，而 fracture contact-model ... activate 没有范围限制，对这些新接触同样生效。
方法二的界面始终固定，接触数还从 105 降到 77——颗粒失去接触就没了，界面只会缩小。
6. 结果五：把界面摩擦归零
到这里有个自然的怀疑：方法二既然和完整试样几乎一样，会不会那个界面根本没接上？
我查了。加载前跨裂隙面的接触 105 个，其中带完好平行键的是 0 个，全部是 linear。破坏后 77 个界面接触里，75 个（97%）有非零法向力，力的大小中位数 2.52e4 N，和同一试样里普通胶结接触的 1.88e4 N 是同一量级，甚至更大。还有 18 个已经压到了摩擦滑移极限。
界面一直都在，而且一直在承力。
那为什么强度不受影响？我做了个只改一个变量的对照实验：方法二的一切不变，只把界面接触的摩擦系数从 0.5 改成 0。
完整试样 123.8 MPa，界面 μ=0.5 时 126.9 MPa（−2.5%），界面 μ=0 时 93.4 MPa（−24.6%）。
把摩擦拿掉，方法二立刻掉到 24.6%，和方法一的 23.3% 几乎相等。
而两次方法二的弹性模量几乎没变（16.90 vs 16.77 GPa）。也就是说，摩擦作用于强度，不影响刚度——界面靠承压传递刚度、靠摩擦抵抗剪切，这个图景和数据完全对得上。
7. 怎么选
两个结论值得记一下。
第一，"零厚度"这个词本身不构成任何保证。两种方法都做出了真正零厚度的界面，宏观结果却能差 10 倍。判断一套建模方法是否合适，不能看它是不是"零厚度"，要看它给出的界面力学性质对不对。
第二，界面是否成为弱面，几乎完全由摩擦决定。方法二在 μ=0.5 时是安全的结构面，μ=0 时立刻变成控制性弱面。如果你的工况里节理摩擦角有变化（比如遇水软化、或有充填物），用方法二时这个参数是敏感项。
具体建议：
想复现真实无厚度节理（可沿面滑移、不承拉、随破坏扩展），用 DFN + smoothjoint。它的法向是可控参数，物理意义清楚，代价是弹性模量会受 sj_kn取值影响，需要标定。只想在模型里放一个"不胶结的面"、又不想动接触本构，用 wall + linear。但要知道它测得的影响会小一个量级，而且结果强烈依赖界面摩擦系数。 无论用哪种，都建议跑一个只改单一变量的对照。这次就是靠摩擦对照才把机理钉死的——不做对照，只看两条曲线几乎重合，很容易得出"界面不起作用"的错误结论。
方法备注：2D PFC 6.0，平行胶结试样 60 × 120 mm，围压 10 MPa，轴压应变控制加载。裂隙物理倾角 45°，贯穿试样。四算例共用同一份成样，球数均为 7948。宏观量（峰值、模量、泊松比、脆性度）用同一套脚本从 history 导出计算，口径一致。
完整报告（9 张图、含全部数据与自检闸门）随代码一起归档在工作目录 post/report.html。
