PFC单轴压缩仿真：非均质模型、声发射统计与自动截图全流程

简介：一套PFC单轴压缩模拟资源面向岩石力学、土木工程与材料科学方向的研究者与学习者，重点解决非均质模型下的单轴压缩仿真与声发射特征分析问题。资源内置完整PFC单轴压缩代码及配套说明，可在模拟过程中根据裂纹数量变化自动截图，并同步输出应力云图与位移云图数据，帮助用户追踪裂纹起裂、扩展与贯通的全过程，进而揭示材料在荷载作用下的破坏机制。资源共11个文件，以6个docx教学说明文档为主体，涵盖模型构建、参数设置与结果解读；辅以4个html可视化页面和1张jpg示意图，便于直观对照。整体压缩包约3.5MB，轻量易用；文档对声发射计数规则、裂纹扩展判定等关键细节作了解释，便于后续二次开发。目前已有249人学习参考，适合需要开展颗粒流离散元模拟或岩石力学数值试验的高年级本科生、研究生及工程技术人员快速上手。 搞PFC的同行应该都遇到过这个需求：想在单轴压缩仿真里把“非均质性”做进去，同时实时统计声发射事件，并且当裂纹数发展到指定数量时，自动截图、把应力云图和位移云图的数据一并导出来。这套东西听起来不复杂，但真做起来，涉及到模型标定、声发射事件定义、Fish编程实现、数据后处理好几个环节，每个环节都有不少坑。今天把这些经验完整拆一遍，代码逻辑和参数选型都会讲到，打算做岩石破裂机理研究、或者想通过数值模拟复现室内声发射试验的朋友，可以直接照着搭。

1. 任务拆解：非均质模型、声发射与按裂纹截图到底在解决什么问题

1.1 为什么非均质模型是这类模拟的刚需

很多刚开始做PFC单轴压缩的人，会直接用均质平行粘结模型跑一组试件。跑完之后发现破坏形态非常单调：试样基本上是“咔嚓”一下，沿某个对角面整齐断开，应力-应变曲线峰后几乎是垂直跌落。这种情况和真实岩石的破坏过程差异很大——真实岩石内部有矿物颗粒差异、微孔隙、胶结强弱不均，裂纹往往从某个薄弱的细观位置开始萌生，然后逐步扩展、汇聚，形成复杂的破裂网络，声发射事件也是从稀疏到密集再到峰值前后的活跃爆炸。

要做这些现象，必须引入非均质性。PFC里最常用的做法是基于Weibull分布，对颗粒间的平行粘结强度和刚度进行随机赋值。为什么选Weibull而不是正态分布？因为Weibull分布非常契合脆性材料的强度统计特征，它有一个形状参数m，可以直观控制强度分布的离散程度：m小，强度分布宽，试样里既有很弱的区域也有很强的区域；m大，强度分布窄，试样趋向均质。这个分布最初就是用来描述脆性材料强度的，用在岩石细观参数赋值上有理论依据。

1.2 声发射在PFC里的落地方式

声发射的物理本质是材料内部微破裂释放弹性波。放进离散元里看，最贴近的对应物就是“微裂纹的产生事件”。PFC里的平行粘结在达到强度极限时会发生脆断，这个断裂过程在微观上释放应变能，宏观上表现就是一次声发射事件。因此，我们在PFC里做声发射模拟，最常见的做法就是把裂纹增量当作声发射计数，把断裂释放的能量当作声发射能量。

这个思路虽然简化，但在工程领域是被广泛接受的。更高级的做法是采用矩张量（moment tensor）反演震源机制，利用裂纹点的位置和作用力方向，反推声发射事件的震级和机制类型。不过多数情况下，统计裂纹增量和能量已经足够分析破裂演化规律了。它天然带上了空间信息——每个裂纹都有坐标，所以声发射定位、时空分布图都能直接做出来。

1.3 “按裂纹数截图”这个需求的真实场景

为什么一定要“按裂纹数”截图而不是按时间固定截图？因为单轴压缩过程中，裂纹增长是不均匀的。弹性阶段基本不涨裂纹，接近峰值时裂纹快速增加，峰后彻底贯通。如果按时间步均匀截图，你会发现弹性阶段几百步内几乎没变化，而破坏阶段几步之内就全碎完了，均匀截图要么无效图太多，要么关键画面全错过。

所以更合理的做法是设定若干个裂纹数阈值，比如每增加200条裂纹就自动保存一次状态，既保证捕捉到从裂纹萌生到贯通的各个阶段，又能让输出文件数量控制在可处理的范围。同时把应力云图和位移云图的数据一起导出，这样后面画曲线、做后处理、对接ParaView都有据可查。

2. 模型构建与细观参数设计

2.1 试样生成和初始平衡

PFC单轴压缩模型的第一步，是把符合目标尺寸和孔隙率的颗粒试样生成出来。常规流程是：先建四面墙，在墙内按指定半径分布生成颗粒，或者先随机生成较大半径颗粒，然后通过半径膨胀法调整孔隙率。膨胀法的好处是颗粒分布均匀，不容易出现局部架空结构。生成后需要用solve aratio命令让试样达到力平衡，把颗粒间的重叠量和系统不平衡力降到很低水平。这一步没做好，后面加载一开始就会“爆炸”式调整，裂纹数根本没法看。

初始平衡完成后，把加载板用伺服控制，通过墙体速度伺服调节加载力，保证加载过程稳定。单轴压缩的加载速率需要反复试，核心原则是：加载过程不产生显著的惯性效应。一个粗略的判据是观察系统动能与应变能的比值，如果峰值破坏瞬间动能突然飙升到总能量的5%以上，说明加载过快，结果可能失真，需要把加载速度降下来。

2.2 Weibull非均质参数如何分配到细观单元

非均质赋值的核心代码思路，是用Weibull分布的逆变换给每一条平行粘结生成一个强度/刚度折减系数。Weibull分布的概率密度函数长这样：

f(x; m, x0) = (m / x0) * (x / x0)^(m-1) * exp(-(x / x0)^m)

其中m是形状参数，x0是尺度参数（与均值相关）。在PFC里实现时，不需要直接调用概率密度函数，而是对它做逆变换采样。对一条粘结，生成一个[0,1]均匀随机数u，然后令系数x = x0 * (-ln(1-u))^(1/m)，再把这个x归一化，作为该粘结参数与基准值的比值。

实操中，我习惯只对pb_ten（抗拉强度）和pb_coh（黏聚力）做非均质赋值，让它们共享同一个Weibull随机系数。这样做的原因是，这两个参数直接决定粘结何时破裂，对裂纹起裂位置影响最大。刚度参数（pb_mod）也可以做非均质，但离散度要控制得小一些，否则试样在加载前就会因为刚度差异过大产生初始应力集中。一个简单的代码框架是这样：

; 设置 Weibull 参数 def set_weibull global m = 4.0 global x0 = 1.0 end @set_weibull ; 遍历所有平行粘结 def assign_weibull loop foreach cp cp.list local u = math.random.uniform(0,1) local factor = x0 * math.pow(-math.ln(1-u), 1.0/m) ; factor 适当归一化，使均值保持在 1.0 附近 local norm_factor = factor / gamma_scale cp.prop("pb_ten") = base_ten * norm_factor cp.prop("pb_coh") = base_coh * norm_factor endloop end @assign_weibull

这里有个细节要特别注意：如果直接用原始Weibull随机数赋值，那这组数的均值不是基准值，而是x0乘以一个与m相关的系数，会导致试样整体强度偏低。所以赋值之前，最好先做一次统计，算出这组随机系数的均值，然后归一化处理。或者，事先用积分算出理论均值并与基准值换算。

从标定经验上看，m的取值一般是1到10之间。m=1~2时强度分布非常分散，试样会出现大量散布的微裂纹，破坏形态偏碎裂；m=3~5比较接近多数硬岩的表现；m=8以上基本接近均质模型了。我自己的习惯是先用m=4跑，对比目标岩石的应力-应变曲线和破坏形态，再微调。

2.3 接触模型选择与细观参数标定

接触模型方面，常规的岩石材料我做两套方案。一套是linear parallel bond，适合比较致密的硬岩，各向同性，标定效率高。另一套是flat joint模型，适合风化花岗岩、裂隙岩体这类需要考虑颗粒表面摩擦和转动阻力的岩体。Flat joint因为接触面积固定，破裂后还能保留部分承载能力，峰后行为更丰富，但参数更多，标定成本更高。

单轴压缩模型要标定的宏观指标通常包括单轴抗压强度、弹性模量、泊松比、破坏形态。标定顺序有讲究，建议按“刚度-强度-破坏形态”的顺序来：先用emod和kratio调弹性模量和泊松比，再调pb_ten和pb_coh匹配峰值强度，最后通过调整非均质系数m和加载条件让破坏形态靠近试验结果。不要上来就同时动四五个参数，不然参数之间的耦合会让你怀疑人生。

3. 声发射统计与自动截图数据输出的实现

3.1 用Fish实现声发射事件实时统计

PFC里声发射统计并不需要额外插件，核心思路就是追踪裂纹总数，每个时步或每隔固定步数计算增量。常用做法是在solve循环外加一个判断，或者用Fish callback在循环过程中周期性执行统计函数。

实际操作中，我是在主循环里这样组织的：

def ae_monitor global current_crack = c_num.text if current_crack > last_crack_total then local new_events = current_crack - last_crack_total local ae_count = ae_count + new_events ; 累计裂纹增量作为 AE 振铃计数 local ae_energy = ae_energy + new_events * avg_release_energy last_crack_total = current_crack endif end

这里的c_num.text是PFC里获取总裂纹数的Fish内建量。avg_release_energy可以简单地取单个平行粘结断裂时的应变能均值，如果有更高精度需求，可以在粘结断裂的回调函数中直接计算该条粘结断裂前后的应变能差。

需要注意的是，声发射时间序列要划分时间窗，不然你只知道总共发生了多少裂纹，不知道它们集中在什么时刻。我通常会把整个加载过程按固定的时步窗口划分，比如每500步统计一次这个窗口内的裂纹增量，保存成“时间-声发射计数”的两列数据，对应室内声发射试验里的振铃计数率曲线。如果再进一步，把每个裂纹的空间坐标按窗口输出，就能做声发射定位的时空演化图。

3.2 按裂纹数阈值自动截图的逻辑

自动截图的关键，是设置一个阈值数组，当总裂纹数跨过某个阈值时触发一次截图和数据导出。伪代码逻辑如下：

def auto_capture local thresholds = array.create(5) thresholds(1) = 200 thresholds(2) = 500 thresholds(3) = 1000 thresholds(4) = 1500 thresholds(5) = 2000 loop while current_crack < max_threshold cycle 200 local new_crack_total = c_num.text loop i from 1 to array.size(thresholds) if new_crack_total >= thresholds(i) then if not captured(i) then ; 设置视图并截图 io.out(string.build("截图：裂纹数达到 %1", thresholds(i))) screen_issue(string.build("crack_%1.png", thresholds(i))) ; 导出云图数据 export_stress_data(thresholds(i)) export_displacement_data(thresholds(i)) captured(i) = true endif endif endloop endloop end

截图前一定要记得先把图形窗口的显示范围设置好。PFC默认的视图可能只在局部，如果窗口范围没锁定，裂纹云图显示出来会忽大忽小，影响截图效果。建议在模型平衡后，锁定视图范围，同时把背景颜色、颗粒配色方案设置成统一标准，保证截图之间具有可比性。

实际操作中，我踩过一个坑：screen_issue截图时如果恰好碰到窗口刷新，可能输出空白或半截画面。后来我在截图前加了几个cycle来解决刷新时机的问题，实测稳定了很多。

3.3 应力云图与位移云图数据的完整导出流程

导出应力云图数据时，PFC本身不直接给你一张云图，而是要先计算每个颗粒的应力张量，然后把颗粒坐标和应力值写出来，再到后处理软件里画云图。经典做法是调用内置的ball stress计算功能，遍历所有ball，提取坐标和应力。代码大致长这样：

def export_stress_data(tag) local filename = string.build("stress_data_%1.csv", tag) local fp = file.open(filename, "write") loop foreach b in ball.list local x = b.pos.x local y = b.pos.y local sxx = b.stress.xx local syy = b.stress.yy local sxy = b.stress.xy file.write(fp, string.build("%1,%2,%3,%4,%5\n", x, y, sxx, syy, sxy)) endloop file.close(fp) end

位移云图更简单，把每个ball的x、y坐标和disp.x、disp.y写出来就行。一个比较推荐的输出格式是包含ball id、x、y、位移分量、应力分量和破坏状态的一张大表。别小看这个表，后面处理的时候方便得很，比如可以在ParaView里用Table to Points把颗粒变成点云，再通过Point Data的应力字段渲染云图，效果比PFC自带的截图更可控，也更适合论文出版要求。

导出频率也是要考虑的问题。每次导出都写所有ball的数据，模型颗粒数到5万以上的时候，文件会越来越大。建议只按事先设定的裂纹数阈值导出，而不是每个加载时步都导出，否则后处理时数据量大到内存不够。

3.4 后处理中的关键技巧

用ParaView或者Tecplot做云图时，PFC导出的原始点云数据是离散颗粒中心点，需要做插值才能形成连续云图。ParaView里最常用的就是Table To Points连接数据源，然后对点数据进行插值成面，再做Cut、Contour这些过滤器。需要留意的是，在插值前最好检查颗粒边缘区域的应力是否为“空心”，导致云图边界不完整。解决方法是把颗粒半径信息导出来，在插值时设置合适的半径范围参数，或者额外生成一层背景网格，把颗粒应力插值到网格节点上再显示。

一个现成的技巧：导出时把每个ball的半径也带上，后处理中直接按半径大小设置点的大小，能让模型显示更立体，比单纯用颜色表达应力更直观。

4. 常见问题与调参实录

4.1 裂纹数不增长或“突然爆裂”怎么办

裂纹数不增长，多半是加载速率过低或者粘结强度太高，导致应力一直低于起裂阈值。这时候可以先输出一个应力-应变曲线看看，如果应力一直在涨而裂纹不涨，说明还在弹性阶段，继续加载即可。如果应力也已经不动了，那就要检查是否早达到了平衡而伺服系统没继续加载。

另一个常见情况是裂纹“突然爆裂”——峰值前一条裂纹都没有，峰值后一瞬间冒出几千条裂纹。这说明模型过于均质，能量积累到一定程度后集中释放。解决思路是降低Weibull分布的m值，引入更多的强度薄弱点，让裂纹在峰前逐渐萌生；或者适当降低加载速率，让破坏过程更充分。

4.2 声发射统计和裂纹数对不上

这个问题基本都出在统计时机上。如果直接用solve命令一口气跑完，中间的裂纹增量都已经过去了，Fish里统计到的只是累计值。正确的做法是分段循环：cycle一定步数，然后统计增量，再继续cycle。我一般用200步作为一个统计窗口，既不会太频繁导致计数开销大，也不会漏掉突变段。

声发射能量统计还有一个容易忽略的坑：如果采用“裂纹数×平均能量”的简化方案，在峰后大破裂阶段会严重高估单次事件的能量。更合理的做法是在粘结断裂时实时计算释放能量，PFC里可以通过cb event监听平行粘结断裂事件，在那个回调里访问断裂时刻的力、位移信息。

4.3 云图数据输出异常或缺失

最常见的问题是导出的CSV文件出现了NaN值。原因通常是某些颗粒在计算过程中已经飞出边界或者已经删除，但导出的循环还在尝试访问它的位移或应力。解决办法是在遍历ball时先判断ball是否有效，或者获取最新颗粒列表再遍历。

位移云图数据异常还有一个特殊原因：如果在加载之前没有归零位移场，PFC会把初始平衡阶段的微小位移也计入。建议在伺服加载开始前，调用一次位移清零操作，确保后续输出反映的是加载引起的真实变形。

4.4 非均质系数m的标定经验表

我在实际操作中积累了几个常用m值的表现特征，可以给大家参考。标定时不需要一步到位，可以先跑小尺寸快速试探，确定m的大致范围，再做完整试样。

4.5 参数标定中的一处关键心得

如果发现试样的泊松比怎么调都不对，先别急着动kratio。多数情况下是边界条件的问题——上下加载板与试样端部之间摩擦设置不合理，导致试样出现明显的鼓胀或端部约束效应。解决办法是在试样端部和加载板之间设置一组低摩擦墙接触，同时保证试样侧向自由，这样泊松比才标得准。我有一段时间标出来的试样弹性模量总偏高，最后发现是试样端部约束过强，导致侧向变形被限制，等效刚度变大，调整接触参数后一下就正常了。

5. 一点实战体会

这套流程我前后迭代过好几版，最花时间的其实不是代码本身，而是参数标定和结果验证。刚开始做声发射统计时，我直接跑完整段solve再去数裂纹，结果发现统计到的裂纹数总是等于总数，完全看不出时间演化规律，改成循环分段统计后才解决问题。后面做自动截图，又遇到了视图范围不统一导致截图之间角度不一致的问题。

如果正在做类似事情，我的建议是先把骨架跑通，用小模型（比如5000个颗粒）验证整个流程——非均质赋值、裂纹数阈值判断、云图数据导出——然后再上大模型正式计算。数据输出命名也尽量提前设计好，比如stress_200.csv、stress_500.csv这样的格式，后面批量处理会省很多事。这套工作流一旦跑顺，整套单轴压缩仿真从加载到破裂分析再到后处理图表，基本就是一键完成的事情了。

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