简介：这套PFC单轴压缩代码采用离散元方法，针对岩石、土体等散体材料构建非均质模型，能够模拟单轴压缩下的裂纹萌生与扩展过程，并同步输出声发射信号、裂纹数截图以及应力/位移云图数据。资源面向岩土工程、岩石力学、土木工程方向的研究者与学生，可用于室内试验对照、破坏机理分析或DEM教学案例。压缩包共有11个文件，其中6份docx文档系统讲解模型搭建、参数标定及结果分析，4个html文件提供可视化交互页面，1张jpg图片直观展示模拟效果，包体仅3.5MB，轻量便于浏览。目前已有255人学习查看，适合快速上手。读者可获得完整的PFC非均质建模思路、声发射统计方法、裂纹演化截图生成流程，以及应力位移云图的数据输出与后处理技巧，对开展颗粒流数值模拟研究具有实用价值。
1. PFC单轴压缩：非均质模型、声发射与裂纹触发的自动截图方案
做过PFC单轴压缩的人都有这种体验：试样快破坏时，一只眼睛盯屏幕上的裂纹扩展，另一只手疯狂按截图键，结果要么截早了裂纹还没贯通，要么截晚了试样碎得不成形。这个项目就是解决这个痛点的——一份完整的PFC单轴压缩代码包，内置非均质模型生成、声发射监测、按裂纹数自动截图，以及应力云图和位移云图数据同步导出。你只需要设置好材料参数和阈值，就能在裂纹数达到设定值时自动保存当前模型状态，同时把对应时刻的应力场、位移场数据写进文件，方便后处理画云图或做定量分析。对于用PFC做岩石、混凝土细观模拟的研究生和工程师，这是可以直接改参数复用的基础模板。
2. 模型从零搭建：非均质试样是怎么“生成”的
2.1 颗粒集合的生成与初始压实
这个项目在PFC 5.0/6.0环境里跑，我习惯用PFC2D建模。试样尺寸是标准的50mm × 100mm，对应PFC默认单位下的0.05m × 0.10m。生成颗粒集合有两个路径：一是用ball distribute命令一步到位，二是在domain内循环创建颗粒再用循环删除悬浮颗粒。项目中常用的是后者，因为方便在生成后立即控制颗粒数量级和孔隙率。
new set random 10001 model domain extent -0.05 0.05 -0.10 0.10 set dt timestep fix 1e-7 ; 生成颗粒前先定义墙（上下加载板） wall create box -0.05 0.05 -0.10 0.10 ; 在墙体包围区域生成颗粒 ball distribute radius 0.0005 0.0008 porosity 0.12 range ... ; 用两阶段压实：先大阻尼跑平衡，去掉悬浮颗粒 solve equilibrium上面代码里，set random 10001是固定随机种子，保证每次建模颗粒位置可复现。porosity 0.12是目标孔隙率，颗粒半径0.5mm到0.8mm对应细观尺度模型。我一般会先用ball distribute快速生成，再用cmat给颗粒和墙定义接触模型。注意ball distribute生成时可能有些颗粒悬浮在空中，必须跑一段equilibrium让它们在重力下落稳，否则后续加载时会突然塌陷——这是很多新手刚跑PFC就翻车的第一个点。
悬浮颗粒去除的做法，常见的是循环检测每个球周围的接触数：
def remove_floating loop foreach b ball.list if ball.contactcount(b) < 2 then ball.delete(b) endif endloop end remove_floating这个函数的核心是ball.contactcount(b)，它返回颗粒接触数量。如果颗粒只有0个或1个接触，就视为悬浮颗粒直接删除。删除后再跑一次equilibrium，让系统重归稳定。这样处理完，颗粒集合才真正具备受压传递的骨架结构。
2.2 非均质性赋值：用 Weibull 分布模拟矿物相
岩石/混凝土内部不是均匀的，因为矿物晶粒、微裂隙、胶结差异等原因，强度在空间上天然分布不均。所以这个项目没有简单给所有颗粒的平行粘结赋一个定值，而是引入Weibull分布来模拟非均质性——这是细观数值模拟里最常用的做法。
Weibull分布的概率密度由形状参数alpha（也叫m）和尺度参数beta决定。alpha越大，强度越集中在均值附近，材料越“均匀”；alpha越小，强度分散，材料越“脆”且不均匀。在PFC里，我们给每个球的平行粘结强度取一个Weibull随机数：
def weibull_random(alpha, beta) local u = math.random.uniform local v = -math.ln(1.0 - u) local s = beta * math.pow(v, 1.0/alpha) return s end def assign_heterogeneity local pb_ten_mean = 8e6 ; 平行粘结抗拉强度基准值（Pa） local alpha = 4.0 ; Weibull形状参数，越大越均匀 local beta = pb_ten_mean / math.gamma(1.0 + 1.0/alpha) local upper = 1.5 * pb_ten_mean local lower = 0.3 * pb_ten_mean loop foreach b ball.list local s = weibull_random(alpha, beta) ; 剪枝：把超出物理范围的值截断 if s > upper then s = upper endif if s < lower then s = lower endif ball.setprop(b, "pb_ten", s) ball.setprop(b, "pb_coh", s) ; 平行粘结粘聚力同上 endloop end assign_heterogeneity说明一下参数：pb_ten_mean是抗拉强度的基准值，比如8MPa，对应中砂岩级别；alpha=4是我做室内单轴测试标定时常用的值，能让强度大致在0.3倍到1.5倍均值之间波动。math.gamma在PFC6的Fish里可以直接调用，如果版本不支持，可以写一个阶乘近似函数替代。
拿到随机强度后，我用upper和lower做了剪枝处理，这就是关键词“剪枝”的落点。为什么要剪枝？Weibull在生成极端值时有长尾效应，偶尔会产生一个几乎是均值5倍的高强度颗粒。这种颗粒在细观模型中会形成一个“超级钉子户”，裂纹绕开它走，导致模拟结果出现不真实的局部硬化。截断到1.5倍均值后，既保留了非均质性，又去掉了统计上不该出现的异常值。剪枝的另一个用途：当你想模拟“含弱层”时，可以对特定区域的颗粒强度整体乘以一个折减系数，比如0.5，形成明显的弱相带。
2.3 材料参数标定：从微观到宏观不是线性换算
这里要给刚接触PFC的读者提个醒：pb_ten、pb_coh这些微观参数并不能直接用实验室的单轴抗压强度来赋值，因为颗粒集合体的宏观强度和微观粘结强度之间存在一个放大系数。常见做法是先跑一个不含裂纹的完整试样做单轴压缩实验，不断调整pb_ten和pb_coh，直到模拟出来的宏观抗压强度、弹性模量与实验曲线接近。
项目代码里一般会把微观参数集中写在文件顶部的“材料参数区”：
; 材料参数 set g_pb_ten = 8e6 ; 平行粘结抗拉强度（剪枝前基准） set g_pb_coh = 30e6 ; 平行粘结粘聚力 set g_pb_fa = 25.0 ; 平行粘结摩擦角 set g_kn = 1e8 ; 法向刚度 set g_ks = 1e8 ; 切向刚度 set g_fric = 0.5 ; 摩擦系数这些值写在set全局变量里，比硬编码在函数里好用得多。后期做参数敏感性分析时，只需要改这里，再跑一遍就可以。记得我第一次跑通时，直接把实验室的单轴压缩强度给了pb_ten，结果模型一压就裂得稀碎，后来才意识到微观参数和宏观强度隔着两个量级。
3. 加载与声发射：单轴压缩过程中都在记录什么
3.1 伺服控制：不伺服，试样直接“炸”给你看
单轴压缩最关键的是加载方式。如果不做伺服控制，墙就会以恒定速度一直往下压，试样在极短时间内积累大量弹性能，破坏瞬间颗粒飞溅，应力波传播失真——模拟就废了。项目里的做法是边加载边测应力，动态调整墙的速度，让加载速率维持在一个准静态水平。
def servo_wall local stress = wstress / (w * thick) ; wstress是墙反力 local sstress = target_stress - stress ; 当前应力与目标差 local vscale = 0.3 wall.vel.x(wp, -sstress * vscale) end这段代码的关键在于wall.vel.x的更新频率。PFC主循环每时步都会调用servo_wall，根据当前轴向应力与目标应力的差值调整墙速度。差值大，墙速度快；差值小，墙速度慢。vscale是速率比例系数，我一般取0.3或0.5，太小收敛慢，太大容易震荡。
完整的加载流程是：先给试样施加一个极小的初始应力（比如1MPa）让试样端面和墙贴紧，再进入伺服加载阶段。加载速度通常用轴向应变速率控制，比如1e-5 /s，也就是每秒应变0.001%。我们设置墙速度时，用strain_rate * sample_height / 2计算上墙速度，下墙固定。注意单轴压缩时试样侧面是自由边界，没有围压，所以侧墙要么删掉，要么设置成只允许法向接触但不给约束力。
3.2 应力-应变曲线：history命令的灵活用法
PFC里记录曲线最简单的是history命令，项目里会当作一个独立模块来处理：
history add id 1 wall.getforce(wp_top) history add id 2 wall.getforce(wp_bottom) history add id 3 measure.getstress(sensor_mid, 1)history add id后面的数字是曲线编号，方便绘图时调用。第3条measure.getstress是从试样中间位置的测量圆读轴向应力，比墙反力换算更稳定。因为墙反力在破坏瞬间会剧烈波动，而测量圆相当于一个“滤波器”，读出来的应力更平滑。
应变记录就不单独用history了，直接通过墙的位置增量换算。常见做法是在一个Fish函数里每N步记录一次：
def monitor local strain = (wall.pos.y(wp_top) - init_height) / init_height local stress = measure.getstress(sensor_mid, 2) io.write("strain,t,stress", strain, time.total, stress) end注意应变用的是上墙位移除以试样原始高度，应力是从测量圆读取的轴向应力。如果后面想算弹性模量，就取应力-应变曲线线性段的斜率；想算泊松比，可以加一个横向测量圆读侧向应变。
3.3 声发射事件：把每个裂纹当成一次“微震”
声发射（Acoustic Emission，AE）在颗粒流模拟里的定义很简单：每一个平行粘结的断裂，就是一次微破裂事件，对应室内试验的一次声发射信号。PFC里没有直接的AE对象，项目通过监听crack事件来实现：
def crack_callback if this.total < crack.count then local n = crack.count - this.pre_count this.pre_count = crack.count ; 记录本次裂纹增量 loop local i = 1 to n local c = crack.list(i) local xpos = crack.pos.x(c) local ypos = crack.pos.y(c) local typ = crack.type(c) ; 1=拉伸,2=剪切 io.write("AE.csv", time.total, xpos, ypos, typ) endloop endif end这里用了一个Fish事件回调，每时步检查crack.count有没有增加。注意crack.list(i)在PFC5.0里的写法是crack.list(i)，在PFC6.0里则是crack.list.index。我写代码时会用一个兼容函数封装，避免换版本时到处改。
记录AE数据时，除了坐标、时间、类型，我还会加一个累计裂纹数crack.count，方便事后按时刻对齐。这个AE.csv文件可以直接导入到Excel或Python里分析，比如画声发射时空演化散点图，或者统计剪切裂纹和拉伸裂纹的占比变化。
4. 裂纹数触发截图与云图数据导出
4.1 设定阈值：按裂纹总数而不是按加载位移
项目最核心的独特点是“根据裂纹数截图”。我想实现的是：当累积裂纹数达到预设的阈值（比如200、500、1000）时，自动保存当前视图的截图，同时导出应力云图和位移云图数据。为什么要按裂纹数而不是按时间？因为不同模型的破坏过程差别很大，按固定时间截的图可能全都集中在破坏前后，或全都在弹性阶段，抓不到重点。按裂纹数可以均匀覆盖“微裂纹萌生-扩展-贯通”全过程。
阈值在代码里用数组维护：
def setup_crack_threshold this.crack_thresholds = list(200, 500, 1000, 2000) this.crack_reached = array(4, 0) ; 记录每个阈值是否已触发 end每次裂纹数增长时检查一遍，如果当前crack.count已经超过了某个阈值，且该阈值还没触发过，就执行截图和数据导出。
4.2 截图命令与Plot切换
PFC的截图核心是plot save命令。我们需要先准备好一个“位移云图”视图和一个“应力云图”视图，然后动态切换并保存：
def snapshot_by_crack local total = crack.count loop local i = 1 to array.getsize(this.crack_thresholds) local t = array.get(this.crack_thresholds, i) if total >= t if array.get(this.crack_reached, i) = 0 array.set(this.crack_reached, i, 1) ; 保存默认模型截图（含裂纹） command plot save current "crack_" + string(t) + ".png" size 800 600 plot name "displacement" select plot save current "disp_" + string(t) + ".png" size 800 600 plot name "stress" select plot save current "stress_" + string(t) + ".png" size 800 600 endcommand export_field_data(t) endif endif endloop end注意plot name "displacement" select这行，我在建模时给不同视图起了不同的名字，比如默认视图显示颗粒和裂纹，另建两个视图分别显示位移云图和应力云图。plot save current会把当前选中的视图保存为PNG文件，文件名为crack_200.png、disp_500.png这种。
有一点必须踩过才知道：plot save current保存的是当前屏幕显示的缩放状态，如果你之前在屏幕上放大过某个区域，保存出来的图就不是全试样。我一般会先plot restore恢复默认视角再截图，或者用plot camera显式设置坐标范围。
4.3 应力云图与位移云图数据导出
截图是给肉眼看，但发论文或做定量分析时还需要原始数据。项目里的export_field_data函数会在每次截图的同时，把测量圆的应力分量和每个颗粒的位置/位移写入CSV文件：
def export_field_data(step) local fn = "field_data_" + string(step) + ".csv" io.open(fn, "w") io.write(fn, "x,y,stress_xx,stress_yy,stress_xy,disp_x,disp_y\n") loop foreach m measure.list local px = measure.pos.x(m) local py = measure.pos.y(m) local sx = measure.getstress(m, 1) ; σxx local sy = measure.getstress(m, 2) ; σyy local sxy = measure.getstress(m, 3) ; σxy io.write(fn, px, py, sx, sy, sxy, 0, 0) ; 位移另从ball读 endloop io.close(fn) end这里有个细节：测量圆只能读应力，位移我是每个球单独读的。更常见的做法是布置一个规则网格（比如5mm间隔），在每个网格位置创建一个测量圆，同时用measure.rel_por来判断该位置是否有接触。位移云图的实际数据可以从ball.disp.x和ball.disp.y取出，在导出时用近邻插值映射到网格点上。
输出的CSV文件第一行是表头，包含坐标、应力分量和位移分量。这样一份文件就同时覆盖了应力云图和位移云图的原始字段。你后处理时可以直接用Python的pandas读进来，用matplotlib.contourf画云图；也可以导入到Tecplot里做更复杂的显示。
为了画图方便，我习惯把测量圆网格也存成一个grid_pos.csv，这样就不会在每次导出时重复计算网格坐标。这个文件在模型生成时就写好了，后续所有时刻都套用同一套网格，保证云图之间可对比。
4.4 自动截图的触发频率与性能开销
最后提醒一个性能问题：如果每个时步都检查阈值列表，模型跑到几万时步时会有额外开销。项目里是每10步检查一次，避免频繁调用crack.count。修改时直接改call fishcallback的步长即可。我实际跑一次从零破坏到完全碎裂，裂纹数到2000，大约需要30到50分钟，其中截图和数据导出的耗时可忽略不计。
5. 避坑/常见问题：参数标定、边界和截图时机里的“血泪教训”
5.1 试样一加载就瞬间碎成粉
现象：伺服加载刚启动，模型颗粒四散飞溅，应力曲线直接归零。
原因：墙的初始速度给得太快，或者伺服系数vscale过大，导致加载波在试样内反射，产生动态扰动。另一个常见原因是颗粒刚度太低，弹性波传播速度慢，墙已经压紧了但试样内部还没反应过来。
解决：把墙速度初值设为0，让伺服函数从前几步的零速度逐步调高；vscale从0.1开始试，观察轴向应力是否平滑增长；同时把时间步dt固定在一个较小的值（如1e-7），避免大时间步导致颗粒穿透墙壁。
5.2 裂纹数从头到尾都是0
现象：模型跑了几十分钟，应力-应变曲线一直线性上升，但crack.count始终不变。
原因：平行粘结强度设置过高，试样在设定的最大位移内根本没达到破坏条件。常见于用pb_ten直接赋了宏观强度的新人。
解决：先单独做一个“单颗粒拉伸试验”——把两个球用平行粘结连起来，匀速拉开，测出实际断裂拉力，再反算微观强度。严格做法是先跑一版不记录裂纹的模型，观察破坏应力；如果峰值应力远高于实验值，就把pb_ten往下降一半再跑。
5.3 截图保存出来的云图是空白的
现象：plot save生成的PNG文件打开后是一张全白/全空图，或者只有几个孤立的点。
原因：PFC的云图显示依赖于测量圆数据，如果你在plot save之前没有激活对应的数据源，云图就不会被渲染。另一个可能是当前视图的显示范围不对，把试样区域的云图截成了空白区域。
解决：在建模时先创建好两个云图视图，并在视图中显式绑定measure和ball disp数据源。保存前调用plot refresh或者等几个时步再保存，确保渲染缓冲更新。我一般会先手动执行一遍plot save看输出文件尺寸是否正常，再挂到自动触发的流程里。
5.4 声发射记录少了一大半
现象：后处理发现AE数据只记录了大约一半的裂纹增量，尤其是快速破坏阶段丢失严重。
原因：检查函数从crack.list(i)读取时，PFC在单时步内可能新增多个裂纹，但我没考虑同一时步内裂纹列表的动态变化。简单用循环取1 to n时，列表索引会越界或跳过。
解决：改用“前后差值”方法——每时步只记录新增数量，而不是遍历当前的crack.list。或者在事件回调里使用crack.list的迭代器。另一个坑是crack.list(i)的索引在高版本PFC中从0开始，需要先查版本文档再写。
5.5 Weibull 随机种子不固定导致结果不可复现
现象：同一套参数，每次跑出来的应力-应变曲线都略有差异，峰值应力波动明显。
原因：math.random.uniform使用的是全局随机数流，每次启动模型时如果不重置，随机序列不同，非均质分布也随之变化。
解决：在模型最开始执行set random 10001固定种子，并且保证生成颗粒、赋强度的时候没有其他函数额外消耗随机数（比如相机随机扰动之类）。如果要统计规律，就固定同一个种子多跑几次，改种子时记录下随机种子编号，方便复现。
6. 二次开发：从单轴到三轴、批量化和 Python 接口
这套代码的扩展性集中体现在伺服函数和阈值触发函数上。最简单的改法就是把单轴改成三轴压缩：在试样两侧增加左右墙，并为左右墙分别写一套servo_wall_lr，目标应力设为围压值（比如2MPa、5MPa）。这样上下墙负责轴向加载，左右墙负责保持恒定围压，一套代码就变成常规三轴试验模拟。
批量做参数敏感性分析时，我一般会改用PFC的Python API。PFC6.0的Python接口可以直接加载模型文件，循环修改Weibull参数和墙速度：
from itasca import itascacmd itascacmd("set random 10001") for alpha in [3.0, 4.0, 5.0]: for seed in range(3): itascacmd("python import alloy") # 修改全局参数，重建模型并运行 itascacmd(f"set g_alpha {alpha}") itascacmd("call 'build_model.p3dat'") # 结束后提取峰值应力写入汇总表这段Python脚本的关键是每次循环都重置随机种子、重建模型，保证不同参数之间的对比不受试样随机性的干扰。跑完后把每次的峰值应力、峰值应变、最终裂纹模式汇总成一个DataFrame，直接出敏感性分析图。
关于震级，你可以给声发射事件加入震级参数：用裂纹的剪切位移或法向位移增量作为震级标度，在crack_callback里把对应的位移差也写进AE.csv。这样能更真实地模拟室内试验的AE事件大小分布，也比单纯统计裂纹数更有说服力。
最后说一个我自己的习惯：每次在自动截图触发之前，我都会先手动跑一小段模型，用plot save验证一遍截图路径和文件名没问题，再谈批量跑。曾经有一次因为工作目录环境变量没配对，程序跑了一晚上，几千张PNG全都写进了临时目录，第二天发现文件全部丢失——从那以后，我每次跑新模型都强制走一遍“先截两张图再正式模拟”的流程，再也不把自动截图当作黑匣子直接信任了。希望帮到你。
本文还有配套的精品资源，点击获取
