简介：一套面向岩石力学、土木工程与材料科学研究者的PFC单轴压缩模拟资料包。资源围绕非均质模型构建与声发射分析展开，提供完整的PFC代码和配套文档，可模拟岩石在单轴压缩下的裂纹萌生、扩展与贯通全过程，并同步输出裂纹数截图、应力云图与位移云图数据，帮助研究者在细观尺度上直观理解岩石破坏机制。资料共11个文件、约3.5MB，以6个Word说明文档作为主体，涵盖非均质模型参数设定、声发射特征提取、应力与位移数据处理等内容；另有4个HTML可视化页面和1张示例图片，便于快速查看模拟效果和云图结果。目前已有253人学习下载，适合具备离散元或PFC基础、正在开展岩石压缩模拟研究的师生与工程师参考使用。整份资料将代码、说明与可视化结果融为一体，既可用于复现单轴压缩试验，也可作为构建非均质模型、分析声发射特性的实用范本。

1. PFC单轴压缩代码，最怕的不是跑不动，而是跑完了发现过程没留下

PFC单轴压缩代码，最怕的不是跑不动，而是跑完了发现过程没留下。很多人调好非均质模型、标定完参数，一加载就盯着屏幕等第一条裂纹、等峰后破坏，等到手忙脚乱截图的时候，应力云图没保存、位移数据没导出、声发射事件也没有按时间记录。这套脚本方案解决的就是这件事：在PFC里建一个颗粒级非均质岩样，加载过程中把声发射当作裂纹事件增量记录，同时按预设的裂纹数阈值自动触发截图，并同步输出应力和位移云图数据。适合两类人：一类是刚把PFC跑通、正被数据输出折腾的新手，另一类是批量跑参数扫描、需要过程图的熟手。

2. 非均质模型：让试样从“玻璃块”变成“岩样”

2.1 均质模型为什么压不出“真实岩样”的破坏曲线

很多人用PFC第一步建模的时候，所有平行键强度给同一个值，跑出来的单轴压缩曲线非常“干净”：线弹性段结束，应力到达峰值，然后一下子跌掉一半。这张图放论文里看着整齐，但和真实岩芯的实验曲线放在一起对比，问题立马上来——真实试样峰值前就有大量微破裂，应力曲线不是光滑直线，而是有小阶梯、有小起伏，峰值后也不是一次性失稳，而是分步跌落。原因不复杂：真实岩样是矿物颗粒的集合，颗粒大小不同、胶结强度不同、微裂隙位置随机，这些离散特征决定了破裂是从局部弱点开始，而不是全断面同时达到极限。

均质模型最直接的后果是声发射没东西可看。盯住裂纹计数器，前面几千步一个裂纹都没有，到峰值前几十分之一秒“唰”地一下冒出几百条，然后模型就过去了。这种事件序列既算不出b值，也没法做前兆分析，因为事件数量在前兆阶段太少，定位出来全是同一条带上的点。做单轴压缩的非均质模型，核心目的不是让曲线“更像”，而是让破裂过程有足够的前兆信息可以记录，让声发射事件在时间和空间上都有分布可查。

2.2 三种非均质实现方式，以及为什么优先选平行键

给模型加非均质，PFC里常见有三条路：

颗粒级非均质改变的是颗粒本身的刚度、密度或摩擦系数，对破裂模式的影响相对温和，更多用于模拟材料弹性的波动。区域级非均质需要把试样划分成多个区块，适合有明确结构面的情况，比如层理、矿物颗粒边界，但实现成本高，参数也难收敛。单轴压缩模拟完整岩样的破坏，我一般优先选接触级非均质，也就是对平行键强度做随机分布。原因有两个：一是接触强度直接决定了裂纹在哪里萌生，对声发射事件的空间分布影响最大；二是实现成本低，不需要重新划分区域，建模流程不变。

2.3 用Weibull分布撒强度：FISH赋值代码和参数边界

接触级非均质的标准做法是用Weibull分布给平行键强度做随机赋值。Weibull分布的形状参数m是这里的灵魂：m越大，强度分布越集中，模型越接近均质；m越小，低强度接触占比越高，试样越容易提前出现局部损伤。常用的逆变换采样公式是sigma = sigma0 * (-ln(1-u))^(1/m)，其中u是(0,1)均匀随机数，sigma0是特征强度。

; 生成Weibull随机数，m_weibull 和 m_sigma0 为全局参数 fish define weibull_sample local u = math.random.uniform weibull_val = m_sigma0 * (-math.ln(1.0 - u)) ^ (1.0 / m_weibull) end

这个函数每次调用产生一个满足当前m和sigma0分布的强度值。赋值到平行键时，我习惯先给每个ball的extra槽放一个随机强度，再遍历接触取两端ball的均值写入平行键抗拉强度。用平均值比用单端值更平滑，不容易出现相邻几个接触强度同时极低导致局部瞬间崩掉的情况。

; 将Weibull强度投影到每个平行键 fish define assign_hetero_pb local cl = contact.list local it = list.begin(cl) loop while it # null local cid = list.get(it) ; 只处理平行键接触：pb_ten 属性存在说明是平行键 if contact.prop(cid, 'pb_ten') # null then local b1 = contact.end1(cid) local b2 = contact.end2(cid) local s1 = ball.extra(b1, 1) local s2 = ball.extra(b2, 1) ; 取两端颗粒强度均值作为接触强度 contact.prop(cid, 'pb_ten') = 0.5 * (s1 + s2) endif it = list.next(it) endloop end

代码逻辑是先在ball层面把非均质“撒”下去，再投影到接触层面。这样每个ball只要生成时初始化extra槽一次，后续遍历接触时不需要反复生成随机数。调参数的时候，形状参数m默认我取5左右，模拟脆性岩石m取4到6；模拟强度离散程度大的软岩或者含初始缺陷的试样，m降到2到3，效果更明显。还要注意随机数种子，同样的参数、不同seed，破坏路径会不同。这不是bug，批量做参数扫描时正好用这个特性探测模型对非均质的敏感度。赋值要在伺服完成之后、正式加载之前进行，赋值后再跑一小段伺服重新平衡，否则伺服阶段经过大量接触力调整，非均质分布的意义会被洗掉。

3. 声发射模拟：没有传感器，就把裂纹当信号

3.1 PFC里的声发射不是波，是破裂序列

真实声发射试验里，传感器接收的是弹性波信号，波源就是材料内部裂纹快速扩展的那一瞬间。PFC是显式离散元，颗粒之间没有连续介质意义上的“波场”，它算的是颗粒运动，所以不能在PFC里直接装一个声发射传感器去“听”波形。但PFC里有一个真实试验拿不到的东西——每条裂纹产生的时间和位置。裂纹本身就是声发射事件的物理源头，把裂纹的时间、位置、类型按出现顺序记录出来，再做事件聚类，就能在统计意义上还原AE事件的时间序列和空间分布。

这里要明确边界：PFC模拟AE，得到的是事件的时空位置和相对频次，不是真实波形幅值。有人在论文里直接写“PFC模拟声发射能量”，严格说是拿裂纹数量或者接触断裂时的相对动能做代理，不是传感器标定后的能量。做分析时用裂纹数、事件数、事件密度这些相对指标没问题，拿去和真实声发射振幅、能量做绝对对比会站不住脚。

3.2 增量监听裂纹：用FISH把每个新生裂纹记进文件

记录裂纹最笨但最可靠的办法，是在每个计算循环中检查裂纹总数，和上一次记录的值比较，多出来的就是这段时间新增的裂纹。PFC里裂纹对象有全局递增id，所以监听逻辑可以写成：遍历裂纹列表，凡是id大于已记录最大值的就是新的，追加写入CSV。

; 增量监听裂纹并追加到csv fish define log_new_cracks local cl = crack.list local nc = list.size(cl) if nc <= crack_known then exit endif local f = file.open('crack_track.csv', 'a') local it = list.begin(cl) loop while it # null local cid = list.get(it) if crack.id(cid) > crack_known then local cx = crack.pos.x(cid) local cy = crack.pos.y(cid) local cz = crack.pos.z(cid) local ctype = crack.type(cid) file.write(f, string.build('%d,%d,%f,%f,%f,%d\n', step, crack.id(cid), cx, cy, cz, ctype)) endif it = list.next(it) endloop file.close(f) ; 更新已知裂纹id阈值 crack_known = nc end

这个函数的关键是每轮只追加新增记录，不做全量重复写。文件打开模式是'a'，即追加模式。最后把crack_known更新为当前裂纹总数，下一次再调用时新的裂纹id更大，不会写重。step字段是当前计算步，可以当作时间轴；如果关心真实时间，需要结合timestep换算。crack.type返回值区分张拉裂纹和剪切裂纹，具体数值不同版本可能不同，代码里先原样存下来，后处理再映射。

这个监听函数不要每个step都调用，效率太低。我的做法是每50到200个计算步调用一次。对单轴压缩来说，200步级别的时间分辨率足够，文件也不会膨胀得吓人。

3.3 把裂纹序列聚成AE事件，再用Python算b值

CSV里存的是一条条裂纹，但声发射分析关心的是“事件”——一次微破裂产生的那一组关联裂纹。把裂纹聚成事件需要一个简单可跑的聚类逻辑：时间上相邻、空间上相近的裂纹归为同一个事件。PFC里不需要做复杂聚类，一个前向链式判断就够用。

import pandas as pd import numpy as np df = pd.read_csv('crack_track.csv', names=['step','id','x','y','z','type']) df = df.sort_values('step').reset_index(drop=True) th_t = 100 # 时间窗：相邻裂纹计算步差阈值 th_d = 2.0 # 空间窗：事件半径取平均粒径的倍数 radius = 0.25 # 平均粒径，由PFC模型已知 event_no = 0 events = [] cx, cy, cz, ct = 0, 0, 0, 0 for idx, r in df.iterrows(): if idx == 0: event_no = 1 cx, cy, cz, ct = r.x, r.y, r.z, r.step else: dt = r.step - df.loc[idx - 1, 'step'] dd = np.sqrt((r.x - cx) ** 2 + (r.y - cy) ** 2 + (r.z - cz) ** 2) if dt <= th_t and dd <= th_d * radius: # 属于同一事件，更新事件中心 cx, cy, cz, ct = (cx + r.x) / 2, (cy + r.y) / 2, (cz + r.z) / 2, r.step else: # 新事件 event_no += 1 cx, cy, cz, ct = r.x, r.y, r.z, r.step events.append(event_no) df['event'] = events ae = df.groupby('event').agg( t=('step', 'min'), x=('x', 'mean'), y=('y', 'mean'), z=('z', 'mean'), n=('id', 'count') ) ae.to_csv('ae_events.csv')

这里的聚类是简化版，时间窗th_t取100步、空间半径th_d取平均粒径的2倍，对单轴压缩的破坏演化足够用。如果想更严谨，可以把参数喂给DBSCAN做密度聚类，但前向链式的好处是稳定、不会把时间上倒序的裂纹扯进来。聚类之后每个事件有出发时间、事件中心、事件内裂纹数，用裂纹数作为事件强度的代理，就可以做Gutenberg-Richter关系的近似回归。

要注意的是，这个b值是基于事件内裂纹数量的统计分布，不是真实震级。对比不同组的相对变化有参考价值，绝对不能把绝对值当成真实岩石的b值去和实验数据做定量对比。PFC的声发射分析定位在“相对趋势”上，拿它做前兆指标、做组间对比，比拿它做绝对标定靠谱得多。

4. 按裂纹数触发截图与云图数据导出：自动化的核心

4.1 主循环里放一个轮询函数，别用model solve一把梭

有人会问，用model solve直接跑到破坏，结束后再补一张图不行吗？单看最终云图当然可以，但只要做过一次参数扫描就会发现，这样出来的图只有“结果”没有“过程”，更没办法把某个裂纹数和某张云图对应起来。标题里说的“根据裂纹数截图”，本质是让求解过程在裂纹数跨过阈值时暂停一下，执行截图和导出动作，然后继续加载。实现方式就是自己控制主循环，每隔一段计算步检查裂纹数是否到达下一个目标值。

; 主循环：每200步检查一次阈值，触发输出 fish define run_ucs loop while list.size(crack.list) < crack_stop command model step 200 endcommand threshold_output() endloop end

loop while的条件是裂纹总数小于终止阈值crack_stop，也就是把加载过程结束条件绑定到裂纹数而不是应变。好处是脆性破坏一旦发生，裂纹数会迅速增长，循环能在破坏后不久自然退出，不需要事先知道峰值位置。model step 200放在command块里，意思是每轮跑200个计算步，跑完就交回给FISH做判断。阈值检查函数threshold_output放在每一次步进之后，驱动后续的截图和数据导出。

4.2 三种输出要同时触发：应力云图、位移云图、场数据CSV

threshold_output是整套脚本的核心动作函数，它做三件事：截图、导出云图数据、更新阈值。截图输出应力云图和位移云图各一张，数据导出则把每个ball的坐标、应力分量、位移分量写入CSV，保证图和数一一对应。

; 当裂纹数达到阈值时，截图并导出场数据 fish define threshold_output local nc = list.size(crack.list) if nc < out_target then exit endif ; 1) 截图：应力云图和位移云图 local f_stress = 'snap_stress_' + string(out_target) + '.png' local f_disp = 'snap_disp_' + string(out_target) + '.png' command plot output '@f_stress' width 1600 height 1200 plot output '@f_disp' width 1600 height 1200 endcommand ; 2) 导出ball应力位移场数据 local f = file.open('field_' + string(out_target) + '.csv', 'w') file.write(f, 'id,x,y,z,sxx,syy,szz,ux,uy,uz\n') local b = ball.first loop while b # null local id = ball.id(b) local x = ball.pos.x(b) local y = ball.pos.y(b) local z = ball.pos.z(b) local sxx = ball.stress.xx(b) local syy = ball.stress.yy(b) local szz = ball.stress.zz(b) local ux = ball.disp.x(b) local uy = ball.disp.y(b) local uz = ball.disp.z(b) file.write(f, string.build('%d,%f,%f,%f,%f,%f,%f,%f,%f,%f\n', id, x, y, z, sxx, syy, szz, ux, uy, uz)) b = ball.next(b) endloop file.close(f) ; 3) 更新阈值，进入下一档 out_target = out_target + threshold_gap end

这段代码里有三处要注意。文件名带裂纹数，每档输出都是独立文件，不会互相覆盖。ball.stress.xx这类属性接口是PFC 6.0之后的写法，如果你的版本老，改成对应的鱼眼函数读取颗粒平均应力，字段含义不变。CSV里导出的应力和位移就是云图画面的原始数据，后处理时直接在外部绘图工具里复现，不用再去PFC里重新截图。

threshold_gap的取值决定了输出密度。我一般设置成100条起步，之后每隔1000条裂纹输出一次。峰前阶段裂纹增长慢，100条一档能捕捉到起裂和扩展过程；峰后裂纹爆发式增长，1000条一档足够看到破坏面贯通。如果你关心声发射事件密集区，可以把前段间隔缩小到200条，后段间隔放大到2000条，避免文件刷爆。

4.3 固定视角与云图配置：不然前后截图对不上

截图输出最容易被忽略的是视角和配色范围。PFC里如果不显式固定camera位置，窗口会自动根据模型范围调整视角，前后两张图可能错位；云图颜色范围如果不去固定，峰前和峰后云图色带会随最大最小值自动变化，“红色代表高应力”这件事前后就不一致了。做过程对比图，这是最大的坑。

我一般在加载前把视角锁死，颜色范围也固定下来：

command plot set camera position 0 0 500 plot set camera direction 0 0 -1 end

视角设置按你的模型尺寸调整，关键是“设置一次之后不要改”。云图颜色范围以完整试样加载前的应力状态为基准，把color range的上下限固定成一个合理的应力区间，之后所有截图都用同一套色标。这样拼出来的序列图，颜色变化才是真实的应力演化，而不是每一张图各自按最大最小值自动归一化。

还有一个细节：云图里如果还叠加了裂纹，要把裂纹符号大小固定。PFC里裂纹符号会随窗口缩放变化，不同截图之间裂纹显得忽大忽小。统一缩放比例，或者干脆导数据之后在外部绘图里叠加裂纹位置，比在PFC里直接渲染更可控。

5. 非均质PFC单轴压缩最容易翻车的5个坑

5.1 裂纹数瞬间从几十涨到几千，曲线直接崩掉

现象：加载过程中裂纹数长期为零，到了某个时刻突然从几十条一口气跳到几千条，应力-应变曲线的峰后段几乎垂直跌落，完全看不到破裂渐进扩展的过程。

原因：加载速率太高或者阻尼设置不合理。PFC里加载速率过快时，动能来不及耗散，颗粒体系积蓄的弹性应变能在短时间内集中释放，裂纹呈雪崩式产生，这不是岩石的脆性破坏特征，而是数值体系的惯性效应。

解决：降低加载速率，比如把加载墙速度从每步1e-4量级降到每步1e-6量级，同时检查局部阻尼。PFC的局部阻尼默认值对单轴压缩来说偏保守的话，裂纹爆发会提前；阻尼太大则会抑制裂纹扩展。一般从默认值开始，用小模型先试出“峰前无裂纹、峰前几十条微裂纹”的合理区间，再放大到正式模型。

5.2 伺服结束后赋非均质，模型反而“散架”

现象：在伺服平衡阶段结束后，调用assign_hetero_pb给平行键赋强度，然后继续加载，发现模型还没到预期峰值就出现大量张拉裂纹，强度明显偏低。

原因：伺服阶段已经让颗粒之间建立了接触力和重叠量，此时再改平行键强度，相当于把受力状态下的“胶结”突然换成了另一个强度分布。本来处于临界状态的接触在赋值后直接失效，产生一批人为裂纹。

解决：把非均质赋值放到伺服之前，或者赋值后再做一段低速率伺服重新平衡。我习惯的顺序是：先生成颗粒和接触，再用均质参数伺服到平衡，然后暂停，执行assign_hetero_pb，最后以更低的墙速度跑很短一段伺服让体系适应新的接触强度分布，再进入正式加载。赋值后不建议用大步长或高速伺服，容易引入额外扰动。

5.3 云图颜色范围自动变化，前后截图不同色

现象：峰前截图里红色区域集中在试样中部，峰后截图红色区域好像到了端部，但对比数值发现中部应力明明没有降那么多——纯粹是每一张图的色带范围各自缩放导致的视觉错觉。

原因：PFC云图默认color range跟随当前模型的最大最小值更新。峰值前全场应力均匀，最大应力可能只有20MPa，云图显示很丰富；峰值后局部应力集中到50MPa，色标上限被拉高，之前的20MPa就变成了蓝色。前后图像无法做视觉对比。

解决：加载前固定云图色标范围。把color range上下限设为加载过程中预期的应力上下界，或者用伺服平衡完成后的最大应力作为上限，全过程中保持这个设置。位移云图也一样，提前估算端部位移极值，固定色标，输出的序列图才具有对比意义。

5.4 AE事件全部集中在加载板端部，试样中段反而很安静

现象：聚类后的AE事件位置图里，绝大多数事件分布在试样上下端部附近，中段破裂稀疏。这看起来像“试样两头在破坏”，和真实单轴压缩破坏模式对不上。

原因：端部效应。加载墙与颗粒之间的摩擦和刚度约束会导致端部颗粒受力状态异常，形成人为应力集中带，裂纹从那里优先起裂。

解决：加载墙用刚性墙时，把墙的摩擦系数调低，或者干脆用无摩擦墙；端部颗粒的刚度不做特殊非均质赋值，保持均匀。如果端部效应还是压不掉，后处理时把端部一定厚度范围内的AE事件排除掉，只统计试样中段的事件。这种截断处理要在论文里写清楚，否则审稿人会质疑数据可靠性。

5.5 截图时PFC卡死或图片文件损坏

现象：阈值触发截图后，PFC卡住几秒钟，或者生成的PNG文件打不开，文件大小为零。

原因：输出动作和求解循环没有错开。截图本身要消耗渲染资源，如果在高频率弹跳检查时频繁触发，或者文件名重复导致文件被占用，就容易出问题。有些版本下，屏幕窗口被遮挡或最小化时，截屏渲染也会异常。

解决：控制输出频率，阈值间隔不要小于500条裂纹，频率太低文件刷爆反而拖慢求解。文件名务必带裂纹数、每次都不同，避免覆盖正在写入的文件。PFC如果支持离屏渲染，就在配置里打开，不依赖系统当前窗口状态。截图触发后可以把文件打开和写入动作放到所有截图完成之后，避免穿插操作互相干扰。

6. 把截图和曲线串起来：一个能直接复用的云图复盘技巧

自动输出做好之后，手上会有一堆带裂纹数命名的PNG和一串CSV。怎么把这些东西变成能放进论文或者汇报里的“演化过程图”？我的做法是写一个Python脚本，把截图拼成GIF，同时把应力、裂纹数曲线和每个截图的位置对上。

import glob import pandas as pd import matplotlib.pyplot as plt from PIL import Image # 按文件名里的裂纹数排序 files = sorted(glob.glob('snap_stress_*.png'), key=lambda f: int(f.split('_')[-1].split('.')[0])) # 拼GIF，每帧400ms，可以直接拖进PPT frames = [Image.open(f) for f in files] frames[0].save('ucs_stress.gif', save_all=True, append_images=frames[1:], duration=400)

曲线对齐的原理很直接：文件名里的裂纹数就是横坐标。读每张图对应的CSV文件里的应力列，取轴向应力和裂纹数的关系画在一张图上，再用标注点标出每个截图的位置。这样一页图同时展示应力-应变曲线、裂纹增长曲线和对应破坏阶段的云图，比单独贴三张图直观得多。

这个复盘技巧对参数扫描特别管用。我后来做多组参数对比时，会把所有GIF放到同一画幅，逐帧对比同一裂纹数下不同参数模型的破坏形态。有一次某组参数跑出来声发射事件数比默认参数多了三倍，原本怀疑是参数问题，用这个逐帧对比方法检查发现是加载速率没调一致，端部先碎了。如果没有按裂纹数对齐截图，这个问题几乎不可能在后期检查数据时看出来。现在我每一个PFC脚本里，条件触发输出永远在solve之前写进去，边加载边落数据，不再有任何“跑完再补图”的侥幸操作。希望帮到你。

本文还有配套的精品资源，点击获取