简介一套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之间。m1~2时强度分布非常分散试样会出现大量散布的微裂纹破坏形态偏碎裂m3~5比较接近多数硬岩的表现m8以上基本接近均质模型了。我自己的习惯是先用m4跑对比目标岩石的应力-应变曲线和破坏形态再微调。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的大致范围再做完整试样。m值强度离散程度典型破坏特征适用岩石类型1~2极高峰前大量微裂纹破坏面弥散整体呈渐进破坏软弱岩、煤岩、结构性岩体3~5中等峰前可见微裂纹萌生峰后呈剪切带或张拉-剪切复合破坏多数硬岩、花岗岩、砂岩6~10较低峰前裂纹很少脆断特征明显破坏面较单一致密均质岩石、高强度混凝土4.5 参数标定中的一处关键心得如果发现试样的泊松比怎么调都不对先别急着动kratio。多数情况下是边界条件的问题——上下加载板与试样端部之间摩擦设置不合理导致试样出现明显的鼓胀或端部约束效应。解决办法是在试样端部和加载板之间设置一组低摩擦墙接触同时保证试样侧向自由这样泊松比才标得准。我有一段时间标出来的试样弹性模量总偏高最后发现是试样端部约束过强导致侧向变形被限制等效刚度变大调整接触参数后一下就正常了。5. 一点实战体会这套流程我前后迭代过好几版最花时间的其实不是代码本身而是参数标定和结果验证。刚开始做声发射统计时我直接跑完整段solve再去数裂纹结果发现统计到的裂纹数总是等于总数完全看不出时间演化规律改成循环分段统计后才解决问题。后面做自动截图又遇到了视图范围不统一导致截图之间角度不一致的问题。如果正在做类似事情我的建议是先把骨架跑通用小模型比如5000个颗粒验证整个流程——非均质赋值、裂纹数阈值判断、云图数据导出——然后再上大模型正式计算。数据输出命名也尽量提前设计好比如stress_200.csv、stress_500.csv这样的格式后面批量处理会省很多事。这套工作流一旦跑顺整套单轴压缩仿真从加载到破裂分析再到后处理图表基本就是一键完成的事情了。本文还有配套的精品资源点击获取