给岩石加软硬互层层理：一个 PFC 技能，跑通单轴、巴西劈裂、直剪

给岩石加软硬互层层理：一个 PFC 技能，跑通单轴、巴西劈裂、直剪

先说结论：层理不要切几何——切开会掉颗粒、破坏接触平衡。正确做法是分层分组 + 三槽 CMAT：把颗粒按几何打上软/硬标签，再用三个 CMAT 槽让"层内软、层内硬、层间界面"各自拿到不同的胶结参数。这套东西与实验类型完全解耦，所以同一套代码能跑单轴、巴西劈裂、直剪三种试样。

1. 一、为什么不切几何

做层理岩石的数值模拟，第一反应往往是"把试样切开、层与层之间留缝"。我在 PFC 里试过，问题很实在：切割会丢颗粒，切口附近的接触平衡被破坏，重新平衡后试样轮廓已经不是你想要的形状了。

换个思路——试样始终是一整块颗粒，层理只体现在两件事上：

• 
  
• 颗粒属于哪一层（分组，纯几何判定）；
• 层与层之间的接触拿什么胶结参数（参数分配）。

具体分四步：按试样几何算出层理面（存成点 + 法向量）→ 在每个层理面位置插一条零厚度墙、平衡后两侧颗粒被推开 → 按几何公式把颗粒分组成 soft / hard、然后删墙（两侧回贴，层间接触天然存在）→ 用三槽 CMAT 施加胶结。

层理面用点加法向量描述，好处是任意倾角、任意层厚都一套公式。法向取 n = (sin α, cos α)，球心在法向上的投影 s = n·x，层号就是 1 + #{ j : s > s_j }。厚度用分数归一化到试样在法向上的投影跨度——这样"厚度"天然就是垂直层面量的真厚度，与倾角无关。

2. 二、三槽 CMAT：让"层间"自动落位

这是整个方法最巧的一步。

我原本想让层内软、层内硬、层间界面各拿一套参数，直觉写法是给 cmat default 加个范围条件。但 PFC 6.0 里这条路走不通：contact cmat default 根本没有 range 位，能带范围的只有 contact cmat add（而且范围必须放在最后），关键字还是 matches（复数）不是 match。

真正的解法是利用槽优先级——PFC 会先匹配可选槽 1、2……只有当所有可选槽都不匹配时才落到 default 槽。于是：

前两个槽只匹配"同组"的接触对，所以软–硬接触两个槽都不中，自动落到 default 槽——一行 range 都不用写。

这里有个必须注意的细节：承压板、剪切盒这类 ball-facet 接触要单独设、而且保持 model linear。因为 contact method bond 对不支持胶结的接触模型是静默忽略的——只要 ball-facet 槽不是 linearpbond，边界就永远不会和试样成键。要是设错了，边界会参与受力，峰值强度直接虚高。

3. 三、两个一等参数

技能对外只有两个主要参数：

• 
  
• --dip：层理倾角，0 表示水平层理；
• --thickness：真厚度（米），层数由它自动派生——nl = round(span / thickness)。

span 是试样在层理法向上的投影跨度。矩形试样是 W·|sin α| + H·|cos α|，而圆盘试样恒为 2R（各个方向投影都一样，所以巴西劈裂的层厚天生与倾角无关）。

除了这两个，还留了 --n-layers（直接给层数）和 --fractions（不等厚，比如 0.3,0.2,0.2,0.3）。

4. 四、一个技能，三种试样

层理那两步在三种实验里是同一件事，都插在"成样之后、胶结之前"。所以 2bedding 和三槽 3jiajiaojie 三个实验完全共用，--test 只切换成样、边界、加载这三步。

三种试样的几何与同一条 30° 层理：

 

图1　同一条 30°、10 层的层理，画在单轴棱柱、巴西圆盘、直剪方盒三种试样上

圆盘的裁剪比矩形还简单：直线 n·x = s 上离原点最近的点就是 s·n（因为 |n|=1），切向 t = (-ny, nx)，半弦长 L = √(Rout² − s²)，两端点直接写成 s·n ± L·t。闭式解，不用迭代，也没有除法保护的问题。

5. 五、单轴：21 个算例扫出各向异性

用 7 个倾角 × 3 个层厚跑了 21 个算例，全部共用同一份成样——初始状态逐位相同，差异只来自层理。下面这张矩阵图每行一个倾角、每列一个层厚，蓝的是软层、橙的是硬层：

 

图2　单轴扫描全部 21 个算例的分组云图（行 = 倾角，列 = 层厚）。层理在所有算例里都严格交替，并随倾角正确旋转

UCS 随倾角的变化是一个清楚的 U 形：15° 附近最强，60° 最弱，90° 又回升：

 

图3　UCS 随层理倾角的变化（每条线 = 一个层厚）

拿最薄的那组（5 mm，软层占比在全部 7 个角度都是 49.7%~50.1%，倾角是唯一变量）来读：

最强的 90° 与最弱的 60° 相差 2.2 倍（各向异性比 0.454）。

5.1 两个容易读错的量

别直接比界面接触数：它随倾角单调增（5 层时从 454 涨到 896），看着像"越陡越破碎"。其实是竖直层理面在试样里的长度恰好是水平面的 2 倍。除以长度后恒定在约 2240 /m，没有任何角度趋势。

软硬比也不是常数：从一端交替软硬时，层数为奇数则软层会多占一层——0° 时 5 层的算例软层占 60.1%，10 层只有 50.2%。所以"层厚变薄"的扫描会同时改变软硬配比。

5.2 层厚和配比，到底谁说了算

我一开始按层厚分组求平均，算出"层厚影响不到 5%"，结论是配比为主、层厚次要。这个结论是错的——各层厚组覆盖的倾角并不相同（5 mm 组 7 个角度全有，20 mm 组只有 4 个），求平均时把倾角效应混了进来。

改成同角度、同配比的配对比较（只保留三个层厚软层占比跨度小于 3 个百分点的角度）后：

修正后的结论：两个因素都重要，但机制不同。软层占比决定"有多少弱材料"，是单调的；层厚决定"界面约束能发挥多强"，效果依角度而变、并不单调——30°/45° 是薄层更强，60° 反过来是薄层更弱，75° 则是中间层最弱。
所以 要读层厚效应，必须先把倾角和配比都固定住，否则量到的是三者的混合。

6. 六、巴西劈裂

同一个层理步骤，换成圆盘试样加对径加载。10 条层理面把圆盘切成 5 软 5 硬，两端被切出尖角：

 

图4　巴西圆盘的分组云图：10 条层理面把圆盘切成 5 软 5 硬（层理 30°）

 

图5　巴西劈裂的抗拉强度曲线（层理 30° / 10 层）

曲线是典型的巴西响应：线性上升到峰值后脆性跌落。峰值 σt = 2.47 MPa @ 0.279% 应变，累计裂纹 594 条。

和无层理的同一算例比（σt 2.21 MPa @ 0.437%、裂纹 370 条）：层理让峰值提前了约 36%、裂纹多了 60%——界面提供了额外的破坏路径，方向是对的。

7. 七、直剪

直剪的剪切盒是永久边界（不像单轴的成样盒会被删掉），但层理墙用 101 以上的独立编号、只删自己建的墙，两者不冲突：

 

图6　剪切盒的分组云图：层理斜穿上下盒，跨过 y=0 的剪切面（层理 30°）

 

图7　直剪的剪应力-位移曲线（层理 30° / 10 层）

这里要如实说一个问题：直剪算出的剪应力绝对量级偏高（正应力只有 100 kPa，剪应力却到几十 MPa）。我查了无层理的原版算例，它同样如此——起点就是 1.98 MPa、跑完 9.58 MPa 时曲线还在上升。所以这是这个剪切盒算例本身的既有问题，与层理无关。本轮验证的是"层理能正确装进直剪"，剪应力的定量可靠性需要单独再查。

8. 八、踩过的三个坑

① 别软化 ball-facet 刚度：插层理墙时会想把边界墙调软来吸收冲击。但层理墙和成样盒墙共用同一个 CMAT 槽，软化前者必然一起软化后者——实测 310 个颗粒被挤到试样轮廓之外，整个试样向外膨胀，云图上看着就是"外侧松散"。正确做法是不软化，用激进的 calm 清速吸收插入瞬态。

② proximity 成键后必须归零：否则加载中张开的裂缝一旦落回阈值内就会重新生成接触（"裂缝愈合"），接触列表无 界增长、峰后响应失真。归零只改未来的检测距离、不删已有接触，而且归零后别再 model clean。

③ 密度要按实验取：单轴基线用 2.7e3，而巴西和直剪的基线都用了 1000 倍质量缩放（2.7e6）。时间步正比于 √m，取错密度会让时间步小 31.6 倍，同一个算例要多跑 31 倍的循环——症状是"加载怎么跑都不完"。实测误用 2.7e3 时巴西的时间步是 7.8e-8，正确时是 2.47e-6。

9. 九、怎么用

单个算例：

bashpython gen_bedding_case.py --test shear --case-dir "工况/层理直剪" \
    --dip 30 --n-layers 10
python run_pfc.py "工况/层理直剪/main.p2dat" --timeout 21600

"倾角 × 层厚"矩阵加汇总出图：

bashpython run_bedding_sweep.py --root "工况/sweep" \
    --dips 0,15,30,45,60,75,90 --thicknesses 0.02,0.01,0.005
python collect_sweep.py --root "工况/sweep"

批量时所有算例共用同一份成样（先跑 main_base.p2dat 得到 1chengyang.sav 再复 制到各个算例），成样只付一次代价、且各算例初始状态逐位相同。串行执行、失败自动重试。

VibePFC 社区 · 一起把 AI 玩进 PFC
我成立了一个 VibePFC 社群，专门聊 AI 怎么落地到 PFC / 数值模拟——让大模型读懂模型、自己改脚本、自动排查"看不懂"的问题。收一个 50 元的门槛费。
想进群的同学，加 QQ 763388012，备注"vibePFC"即可。

参考资料

1. 
   
2. Itasca PFC 6.0 文档：contact cmat 命令族与 using_cmat/cmat3.p2dat 三槽范例
3. 本文数据来自 21 个单轴 + 1 个巴西 + 1 个直剪算例的全部日志与曲线，报告与高清图见配套仓库