开头

说实话，很多人一上来就跑PFC官方案例，跟着教程把平行粘结模型（parallel-bond model，也就是常说的pb 模型）的代码敲出来，也能跑出“看起来挺对”的应力应变曲线。但真正到自己的课题，需要根据实验数据去标定胶结参数、做破坏模拟的时候，问题就全来了：为什么峰值强度老是差一截？为什么试件破坏形态只有一条斜裂缝，完全没有实验里的粉碎区？为什么调个参数要跑二十几次试算，还是找不到方向？

这篇文章就是围绕这套流程来写的，我会把从胶结参数设置到破坏模拟的完整路径拆开讲清楚：pb 模型那几个关键参数背后到底对应什么物理含义、为什么标定顺序错了会浪费大量算力、单轴压缩和巴西劈裂的模拟应该怎么设置、破坏过程中的裂纹统计和声发射等价输出怎么做，以及我实际踩过的大小坑。

先说清楚一个前提：这里的 PFC 指的是 Itasca 的颗粒离散元软件Particle Flow Code，不是电源电路里那个功率因数校正 PFC。搜索时要小心，特别容易混。下面进入正题。

1. pb模型的细观结构：先搞懂它在模拟什么材料

1.1 平行粘结模型的物理构成

pb 模型在宏观上对应的是一类“有胶结”的脆性材料——岩石、混凝土、砂浆、陶瓷、部分冻土都可以用它来描述。它的核心思想很直观：两个颗粒接触的地方，除了常规的接触力学响应（挤压、摩擦、滑动）之外，还额外生成一个“粘结盘”或者“粘结梁”，这个胶结单元有一定的几何尺寸（由半径乘子pb_radius_multiplier控制）和力学属性（法向刚度、切向刚度、抗拉强度、粘聚力、内摩擦角）。

你可以把它想象成两块砖之间抹了一层水泥，水泥固化之后，砖块之间不仅能传压力，还能传拉力、传剪力。当外力在胶结体内部产生的应力超过胶结强度时，这个粘结点就发生断裂，宏观上对应材料内部微裂纹的萌生。PFC 中裂纹一旦生成，就是不可恢复的，这个特性和脆性材料损伤的不可逆性是对应的。

这是 pb 模型和线性接触模型linear的本质区别：线性接触模型不能承受拉力，只能通过颗粒间的互锁和摩擦提供抗剪强度，本质上模拟的是无粘结散体，比如砂土、堆石体。如果你拿线性接触模型去模拟岩石的单轴抗压，颗粒之间稍微一受拉就直接散架，根本不可能形成完整的试件。

1.2 pb模型与其它胶结类模型的选型边界

PFC 中能模拟胶结材料的模型不止 pb 一个，还有linearparallelbond（也就是常说的 flat-joint 出现之前的经典平行粘结）、flatjoint（平节理模型）、smoothjoint、jfr（破裂平节理）等等。我用下来，日常岩石类模拟需求中 90% 的情况选 pb 就够了，但要知道边界在哪。

pb 模型适合模拟的是：细观上颗粒内部为完整岩石块体、破坏基本沿颗粒边界或胶结面发生的材料。比如砂岩、大理岩、花岗岩在中等围压下的破坏，用 pb 模型可以得到比较合理的强度包络线和破坏模式。但如果是颗粒内部本身存在大量微裂隙的岩石，比如某些风化花岗岩、页岩，或者高围压条件下颗粒自身会发生破碎的情况，pb 模型就显得“太干净”了——因为 pb 模型默认颗粒本身是不可破碎的，破坏只发生在颗粒间的胶结点上。

flatjoint（平节理模型）相比 pb 多了“未胶结表面”的存在，颗粒之间的接触面积更大，更接近真实岩石中矿物颗粒的镶嵌结构，所以在模拟裂纹穿过颗粒本体、产生台阶状破坏面时有优势。但 flatjoint 参数更多，标定难度大不少。我自己的习惯是：先跑 pb 模型，如果结果的破坏形态和实验中切穿颗粒的现象明显不符，再升级到 flatjoint。这样分步走，排错更容易。

2. 胶结参数设置前的准备：试件生成与初始配比

2.1 确定粒径、孔隙率与颗粒数量

很多新手一开始就盯着pb_ten、pb_coh的数值往下调，这其实是个误区。在设置胶结参数之前，试件的几何形貌已经决定了很大一部分宏观力学响应。同样的胶结参数，放到孔隙率 0.4 的试件和孔隙率 0.15 的试件里，出来的单轴抗压强度可以差好几倍，因为孔隙率直接影响配位数（每个颗粒的平均接触数），而配位数决定了单位面积内有多少粘结点在承担载荷。

粒径的选取要参考你模拟的对象。一般来说，试件最小尺寸和最大颗粒粒径的比值建议保持在 20 以上，否则试件的破坏面会被颗粒尺寸主导，出现明显的尺寸效应。比如模拟直径 50 mm 的岩石试件，颗粒半径取 1.0~1.5 mm 是比较常见的做法，这样一个试件大概有几万颗粒。如果只是做参数规律性研究，颗粒可以大一点，到了正式标定时再细化。

孔隙率我一般控制在 0.15~0.25 之间。孔隙率太低会让试件过于致密，初始接触数过多，标定出来的胶结参数数值会偏小，而且在生成阶段颗粒反复“拥挤”会导致很大的初始不平衡力；孔隙率太高则试件太松散，泊松比容易偏大，破坏模式也会过于碎裂化。

2.2 试件生成的标准流程与浮粒消除

在 PFC 中生成一个带胶结的岩石试件，标准套路基本是四步走。

第一步，在墙体内按目标孔隙率随机生成颗粒，通过ball distribute或者半径膨胀法把颗粒布置到指定范围。这里要注意，刚生成的颗粒之间允许重叠，但重叠量不能太大，否则在后续循环中会积累巨大的排斥力。

第二步，让颗粒在没有胶结的状态下平衡，消除初始的颗粒重叠引起的弹性斥力。循环到最大不平衡力与平均接触力之比小于某个阈值（我常用 1e-5）就认为平衡完成。

第三步，删除所有“悬浮颗粒”。所谓悬浮颗粒，就是那些接触数为 0 或者只有 1 个接触的颗粒。这部分颗粒在施加胶结之后对承载贡献很小，但在裂纹统计时会制造大量虚假的“独立裂纹”，而且会拉低试件的整体弹性模量。手动检查时可以用 FISH 统计每种配位数颗粒的数量，然后把接触数小于 2 的颗粒删掉再重新平衡一次。

第四步，对所有接触施加平行粘结。执行contactmodel将接触模型替换为parallelbond，并通过property语句将胶结参数赋到每一个接触上。施加完胶结之后，还要再平衡一次，因为胶结的引入会在颗粒接触点上额外附加一个粘结内力，如果直接加载，可能出现初始裂纹。

我实际测过，如果不做浮粒消除，单轴峰值强度最多可偏低 5%~8%，而且破坏形态上会出现一些孤立的“飞散颗粒”，跟实验对不上。这一步虽然麻烦，但值得做。

3. 胶结参数的物理意义与标定顺序

3.1 刚度参数组：弹模与泊松比的决定因素

pb 模型的参数可以分成两组来理解：刚度组和强度组。

刚度组包括颗粒接触的法向刚度kn、切向刚度ks，以及平行粘结自身的法向刚度pb_kn、切向刚度pb_ks。新旧版本 PFC 的命名有差异，PFC 6.0 之后推荐用pb_kn、pb_ks这种形式，老版本是pb_stiffness等等。实际命令里，更常用的是给一个等效模量deformability选项，直接指定emod和kratio，由程序自动分配颗粒和粘结的刚度。我个人更推荐用emod和kratio的方式，因为这样避免手动换算刚度带来的不一致问题，而且和宏观力学参数的对应关系更直观。

弹性模量和泊松比主要由刚度组决定。数值试验中普遍的经验是：增大emod，宏观弹性模量近似线性增大；切向与法向刚度比kratio增大，泊松比增大。这个对应关系基本呈单调变化，所以在标定时可以把它们解耦处理——先用emod对准宏观弹模，再用kratio对准泊松比。这里有一个注意点：kratio对泊松比的影响范围有限，在 PFC2D 里一般只能把泊松比调到某个区间内，如果始终对不上实验值，就要回到孔隙率或者粒径分布上去调整，而不是死磕这个参数。

3.2 强度参数组：拉伸破坏与剪切破坏的竞争

强度组包括平行粘结的抗拉强度pb_ten、粘聚力pb_coh、内摩擦角pb_fa，以及半径乘子pb_radius_multiplier。

这里要重点讲一下它们各自的角色。当粘结受到法向拉应力时，抗拉强度负责“扛”；当粘结受到剪切应力时，粘聚力是主要的抗力来源，内摩擦角提供与正应力相关的附加抗力。所以这两个强度参数之间并不是独立起作用的，它们共同决定了材料在应力空间中的强度包络线。

我见过很多人标定时只调pb_ten，把单轴抗压强度强行调高。这其实是一种不合理的做法，因为单轴压缩下试件的破坏是复杂的剪张混合破坏，峰值强度是拉伸裂纹和剪切裂纹共同演化的结果。如果只拉高抗拉强度，虽然峰值上去了，但破坏模式可能从剪切破坏变成剧烈的脆性张拉破坏，应力应变曲线的峰后形态也会变得特别“脆”，跟实验对不上。

pb_radius_multiplier这个参数也容易被忽略。它控制粘结盘的半径与颗粒半径的比值，默认是 1.0。增大这个值意味着粘结的截面积变大，抗拉和抗剪的承载能力都会提高，同时也会影响裂纹的成核位置。在做微观参数标定时，我会固定它不动，一般在 1.0 到 1.5 之间取值。

3.3 参数标定的迭代逻辑：不要试图一次到位

参数标定最大的陷阱，就是细观参数与宏观响应不是一一对应的。你调pb_ten能影响单轴抗压强度，但同时也影响抗拉强度、影响脆延转变围压、影响裂纹类型比例；你调pb_fa对单轴峰值强度影响不大，但在围压条件下它的作用会非常明显。所以标定一定要有迭代顺序。

我的标准流程是：

1. 根据实验的单轴抗压弹性模量和泊松比，固定emod和kratio，先不去碰强度参数。
2. 用单轴抗压实验标定峰值强度：先设一个初值，比如pb_ten = 5 MPa，pb_coh = 20 MPa，pb_fa = 30°，跑完看峰值，按比例调整pb_coh；如果峰后脆性太强，再适当调高pb_ten与粘聚力的比值。
3. 用巴西劈裂实验（或者直接拉伸实验）标定抗拉强度：巴西劈裂试验的峰值强度对pb_ten最敏感，用这一步把pb_ten固定下来。
4. 用三轴压缩实验标定pb_fa和围压相关的强度包络线。
5. 回到单轴实验复核第一次的峰值强度，因为第 4 步中修改pb_fa后，低围压下的单轴强度也会受到轻微影响。这个过程通常要反复两三轮。

整个标定过程里，每次只动一个参数，记录它对所有宏观响应的影响，形成一张敏感度表。这比盲目随机试算高效得多。下表是我在某次砂岩标定中记录的模板，大家可以参考：

4. 从单轴压缩到巴西劈裂：室内试验的数值复现

4.1 单轴压缩模拟的设置与伺服机制

单轴压缩模拟看起来简单，但要做到和实验室可对比的精度，有几个细节很关键。

首先是加载方式。实验室的单轴压缩是刚性试验机加载，位移控制。数值模拟里最常用的方案是“墙伺服加载”：上下两个加载墙板以恒定速度向试件移动，同时通过伺服机制调节墙板的应力，使得墙板应力始终等于试件端部的目标应力。PFC 官方 fish 库里有现成的server伺服函数，但很多人直接把墙速开到很大，等曲线跑完就完事——这种做法会引入显著的动态效应。

加载速率到底取多大合适？一个常用的判断标准是惯性数I = ε·d / sqrt(P/ρ)，其中 ε 是应变率，d 是颗粒直径，P 是约束压力，ρ 是颗粒密度。惯性数要小于 1e-3 才能保证准静态条件。换算成直观的说法，就是颗粒半径在 1 mm 量级时，加载速度取 0.05~0.5 m/s 基本是安全的。如果你只是定性看破坏形态，加载速度稍大一点没关系；如果要做定量的强度对比，这个值必须控制在准静态范围内。

单轴模拟中建议输出的信息包括：轴压（用墙的接触力除以试件截面积）、轴向应变、侧向应变、裂纹总数、张拉裂纹数量、剪切裂纹数量、粘结断裂能消耗等。PFC 6.0 中可以通过history命令跟踪这些量。

4.2 巴西劈裂的加载方式与有效性问题

巴西劈裂实验在数值模拟中的实现方式比单轴要讲究。实验里试件是圆柱体，加载方向沿直径方向，通过两个相对的加载压条施加压缩载荷。数值模拟中做法类似，在试件两侧设置两个窄加载墙，以恒定速度向中心移动。

巴西劈裂对应的破坏机理是试件中心区域产生横向张拉应力，最终沿加载直径方向劈裂。这个破坏形态对pb_ten极其敏感，所以是标定抗拉强度的首选实验。但要注意，巴西劈裂模拟中如果颗粒粒径偏大，中心区域的张拉应力集中区被颗粒离散化得不够精细，会出现破坏不是从中心起裂，而是从加载点附近直接压碎的情况——这会让计算得到的“抗拉强度”虚高。我建议在巴西劈裂模拟中，加载点附近的颗粒尽量小一些，或者对加载区域做局部细化。

如果你有条件，也可以做直接拉伸模拟，把试件两端通过胶结连接到加载板上进行拉拔。这样得到的抗拉强度更“干净”，但操作稍微麻烦一些。两种方法的结果可以用来相互校验。

5. 破坏模拟的核心环节：裂纹演化、统计与声发射等价

5.1 裂纹分类与演化特征

pb 模型的破坏模拟，核心输出就是裂纹的“出生记录”。每断裂一个粘结，PFC 会记录一个裂纹事件，并根据断裂时刻的应力状态把它归类为张拉裂纹或剪切裂纹。判据非常直接：如果法向应力超过了抗拉强度，就是张拉裂纹；如果剪应力超过了剪切强度包络，就是剪切裂纹。

从岩石力学的角度，这两个裂纹类型比例直接决定了材料的破坏模式。单轴压缩下脆性岩石的破坏通常以剪切裂纹为主，但裂纹萌生阶段其实有大量的张拉裂纹在扩展，两者相互连接形成宏观破裂面。在倾角 60° 左右的单一剪切带破坏案例中，裂纹分布图上能看到一条明显的“裂纹密集带”，张拉裂纹散布在带两侧，剪切裂纹集中在带内。如果你跑出来的模型里张拉裂纹占比异常地高，比如超过了 70%，那就要警惕是不是胶结参数设置让试件过于脆性了。

在 PFC 中，追踪裂纹可以用crack tally命令查看裂纹总数，用crack map输出裂纹数据。如果要统计张拉和剪切的比例，需要在断裂事件中记录：crack_tension和crack_shear这两个计数变量在 PFC 6.0 里可以直接读取。

5.2 声发射（AE）的数值等价与能量追踪

实验岩石力学里会用声发射监测内部微破裂，数值模拟里没有声波，但有更直接的东西——粘结断裂事件本身就相当于一个 AE 事件。通过记录单位时间（或者单位加载增量）内的裂纹数量，就能得到和实验室 AE 振铃计数率类似的曲线。这个曲线在峰值前逐渐活跃、接近峰值时急剧上升，这个特征与室内 AE 试验是一致的。

更精细一点的做法是把裂纹事件按“震级”分类。PFC 中每次粘结断裂释放的弹性应变能可以计算出来，能量释放越大，对应的 AE 事件振幅越大。你可以用 FISH 在crack创建的 callback 事件中读取断裂时刻的力与位移，估算释放能量，然后按能量大小分级统计。这样你不仅能得到 AE 计数率，还能得到 b 值曲线——这可是岩石力学领域发论文的好素材。

能量追踪方面，PFC 内置了几种能量指标：边界能（边界力做功）、应变能、摩擦能、粘结能、动能。在破坏模拟中，判断系统是否稳定可以看动能与边界能之比。如果峰值附近动能突然占了很大比例，说明加载过快产生了动态破坏效应，结果不可信，需要降低加载速率。

5.3 破坏阶段划分：裂纹起裂、扩展与贯通

利用裂纹数量和能量数据，可以将破坏过程划分成几个特征阶段。这个划分对解释模拟结果非常有帮助。

第一阶段是裂纹起裂阶段。这个阶段应力应变曲线仍然基本线性，但已经有个别粘结达到强度极限发生断裂。在实验室中对应的是声发射事件的零星出现。第二阶段是裂纹稳定扩展阶段。裂纹数量开始线性增长，但应力还在上升，试件整体还具备承载能力。第三阶段是裂纹加速扩展阶段，对应峰值前后，裂纹增长速率急剧上升，宏观破裂面开始形成。第四阶段是峰后软化或破坏阶段，裂纹数量增长逐渐饱和，试件承载力下降。

在 PFC 中判断这些阶段的分界点，可以先看裂纹计数曲线，一阶导数发生突变的位置通常就是起裂点；再看应力应变曲线的峰值位置，两者结合就能定位扩展的各个阶段。另外也可以用体积应变拐点法来判断起裂应力与损伤应力。这套分析方法我已经在很多项目里用过了，比单纯看峰值强度有信息量得多。

6. PFC中pb模型实战：从胶结参数设置到破坏模拟的完整流程

6.1 一个完整的最小工作流示例

讲了这么多理论，下面给出一个可以直接上手的 PFC 3D 代码骨架。假设试件是圆柱体，直径 50 mm，高 100 mm，颗粒半径 0.8~1.2 mm，目标孔隙率 0.18。

第一步：定义模型尺寸

; 墙体生成 wall generate box [-0.05,0.05] [-0.05,0.05] [0,0.1] wall property kn 1e12 ks 1e12

第二步：生成颗粒并平衡

; 在墙体范围内生成颗粒 ball distribute radius 8e-4,1.2e-3 porosity 0.18 box [-0.05,0.05] [-0.05,0.05] [0,0.1] ; 设置颗粒密度和阻尼 ball property density 2650 damp 0.7 ; 消除初始重叠 cycle 2000 calm 50 solve equilibrium

第三步：删除浮粒

; 删除接触数小于等于1的颗粒（简化判断） ...

第四步：设置接触模型与胶结参数

contactmodel assign parallelbond contact property emod 5e9 kratio 2.0 pb_emod 5e9 pb_kratio 2.0 ... pb_ten 6e6 pb_coh 20e6 pb_fa 30 pb_radius_multiplier 1.0 solve equilibrium

第五步：实施单轴加载

; 删除两侧墙体，施加轴向加载 wall delete walls range id 1 wall generate ... wall attribute velocity-z 0.1

6.2 加载速率与准静态条件

加载速率 0.1 m/s 对于这个试件尺寸来说偏保守，不过稳定性优先。如果需要加快计算，可以逐渐提升到 0.2、0.3，但在加速时一定要观察动能曲线，如果出现突变就说明进入了动态加载区间。

另一个常用技巧是在加载初期施加一个小的阻尼，比如局部阻尼damp 0.7，这对抑制颗粒在峰值段的弹跳震荡很有帮助。峰值过后如果试件完全破坏，会出现大量高速飞散颗粒，这个时候求解器收敛性会变差，可以在程序中设置当裂纹数量达到某阈值时自动停机，避免无意义的循环。

这一步也可以配合history记录关键曲线，方便后续和实验曲线对比。

6.3 输出与后处理要点

PFC 自带的绘图功能比较基础，我习惯把数据导出到外部程序处理。在 PFC 中可以用 FISH 或者 Python 接口把应力应变数据、裂纹坐标、能量数据导出成文本文件，再用 Python 结合 matplotlib 绘制曲线和裂纹分布图。

要注意，裂纹数据默认记录的是裂纹的几何属性（位置、方向、类型、大小），绘制时用不同颜色区分张拉和剪切裂纹就可以了。一个实用的小技巧是，在导出的裂纹数据中加上“断裂顺序编号”，这样直接在图上按编号着色就能看到破坏是从哪个区域开始、向哪个方向扩展的。

7. 实操中踩过的坑与调试建议

7.1 参数与试件尺寸的耦合问题

第一个要说的坑，就是胶结参数的“尺寸依赖”。同一个pb_ten，放到颗粒半径 0.5 mm 的试件和 1.5 mm 的试件里，得到的宏观强度完全不同。这是因为粘结盘的面积和颗粒尺寸挂钩，颗粒越大，单个粘结的承载面积越大，整个试件的峰值强度也会更高。所以当你改变粒径分布时，原来的标定参数就全废了，需要重新标定。

这个问题的另一个侧面是：如果你在文献里找到一个标定好的参数组合，想直接用到自己的模型里，大概率不能直接用，必须搞清楚文献中试件的粒径和孔隙率。很多新手在这上面栽跟头，花了一个星期找参数，结果发现自己用的粒径是别人的三倍。

7.2 初始浮粒与初始裂纹导致的强度虚低

前面提到浮粒问题，这里再说一个更隐蔽的情况：生成试件后虽然做了平衡，但在施加胶结时，如果某些颗粒间距刚好在临界值附近，粘结生成后会在局部产生过大的内应力，直接造成“初始裂纹”。初始裂纹的存在会让试件的峰值强度偏低，而且会让应力应变曲线从一开始就出现非线性。

检查方法很简单：在施加胶结后立即统计裂纹数量，如果大于 0，就说明试件制备有问题。解决方法是减小颗粒重叠量，或者在施加粘结前略微扩大颗粒间距离。如果反复出现初始裂纹，可以先把弹性模量调低一档，生成成功后恢复，再重新平衡。

7.3 标定过程中的随机种子敏感性

PFC 中颗粒的位置由随机种子决定，不同的 seed 生成的试件几何形貌不同，即使完全相同的胶结参数和加载条件，宏观强度也会有几 percent 的波动。这个波动在标定时很烦人，因为你会分不清强度差异到底来自参数调整还是来自随机波动。

解决办法是每个参数组合至少跑三个不同的随机种子，取平均值作为该组参数的响应值。这样做虽然增加了计算量，但能大幅提高标定效率——否则你在微调参数时可能是在跟噪声搏斗。

7.4 初始裂隙或缺陷的引入方式

模拟含缺陷岩石时，需要在试件中预设裂纹或孔洞。PFC 中常用做法是提前在指定位置删除局部颗粒，形成空腔，或者直接删除部分粘结模拟初始裂缝。但有个坑：直接删除粘结核形成的“裂隙”在加载初期会因为应力集中产生大量寄生裂纹，这些是数值假象而非物理现象。更稳妥的做法是用smoothjoint模型替代目标面上的平行粘结，形成控制面，这样不仅可以控制裂隙的力学参数，还能避免局部应力突变。

8. 一些提高效率的实用技巧

最后分享几个我在实际项目中积累的小技巧。

第一个是“粗粒预标定”。在正式标定细颗粒试件之前，先用较粗的颗粒快速跑一遍参数范围，确定大致的强度趋势和参数敏感方向，再用细颗粒做最终精确标定。计算量可以减少一个量级。

第二个是“单轴压缩和巴西劈裂同步进行”。在标定程序中，同时准备好单轴和巴西劈裂两个模型的脚本，调整参数后并行提交，一次就能获得两组宏观响应。标定过程中时间最贵的不是计算本身，而是人工等待和判断。我这个流程把等待时间压缩了很多。

第三个是谨慎使用pb_radius_multiplier来微调强度。它在数值上确实能影响强度，但同时也改变了粘结的截面积和裂纹的形态，副作用比较大。如果需要微调强度，优先考虑调pb_ten和pb_coh的组合。

第四个是保存阶段的“快照”。每次标定出一个合理的阶段结果，就用save保存一份 p3sd 文件，同时记录当时的参数组合和对应的宏观响应。这样如果后面调整参数越调越差，可以快速回到之前的有效状态，而不是从头再来。

结尾

从胶结参数设置到破坏模拟，这个流程的每一步都有讲究。我自己刚开始接触 PFC 的时候也走过弯路，最大的体会是：数值模拟的“准”不是靠某一个参数调得好，而是靠试件制备、参数标定、加载条件、后处理分析各个环环相扣。如果你现在正卡在参数标定或者破坏形态对不上的问题上，建议不要急着堆算例，先把试件本身的平衡状态、浮粒、初始裂纹这些“底层卫生”检查一遍。底层干净了，上层参数的调整才会真正发挥作用。

最后再分享一个心得：PFC 模型不要太追求和实验结果“完全一致”。离散元本身的随机性和细观机制决定了它更适合做机制分析和趋势研判，而不是精确复现。把破坏模式、裂纹类型比例、强度包络线的趋势抓住，这个模型的价值就已经很高了。希望大家在跑模拟的时候少踩坑，多出活。