简介：面向岩石力学、土木工程与材料科学等领域研究者的PFC单轴压缩数值模拟资源，重点解决非均质模型构建、声发射信号捕捉与破裂过程可视化等问题，适合开展岩石破坏机理分析的中高级研究人员使用。资源包含运行代码、六份docx文档、四份html说明页与一张jpg示意图，docx用于查看原理推导与代码解读，html可快速浏览模拟结论，整体压缩包约三点五MB。代码内置根据裂纹数自动截图功能，可同步输出应力云图与位移云图数据，帮助用户直观追踪裂纹从萌生到扩展的演化路径，并理解非均质材料在单轴压缩下的力学响应。全部文件共十一个，结构组织清晰，方便按需查阅。资料目前已有二百五十三人学习，适合需要借助离散元方法深入分析声发射特征、绘制云图并提取应力位移数据的研究人员与研究生。
1. PFC单轴压缩代码：从非均质模型到声发射输出的完整闭环
PFC单轴压缩代码看起来简单，真正跑起来时却有一堆暗坑：均匀模型的宏观强度曲线偏脆，声发射事件分布完全不合理，等到想导出应力云图才发现坐标系都对不上。这套非均质模型的单轴压缩代码，把颗粒参数按Weibull分布离散赋值，用裂纹生成事件作为声发射源，并在加载过程中按裂纹数自动截取裂纹画面，同时输出应力和位移云图数据。适合用PFC做岩石破裂模拟、需要快速生成试验曲线和云图素材的从业者，也适合研究生课题里拿前期试算数据探路的场景。不需要从零堆代码，拿到后改参数就能跑通一条完整流程。
2. 非均质模型与声发射仿真：先搞清楚三个底层参数再动手
很多人在PFC里搭单轴模型，图省事直接给所有颗粒设同一个强度值，结果应力应变曲线是一条直线式脆断，和真实岩石的全应力应变曲线差了十万八千里。原因在于真实岩石内部矿物晶粒、微裂隙和胶结物强度并不均匀，宏观破裂是渐进式萌生、扩展、贯通的过程。所以非均质模型的第一个关键，是让颗粒层面的力学参数分布符合某种统计规律。
2.1 非均质颗粒参数：Weibull分布与pb粘接强度映射
我一般会建议用Weibull分布给平行粘结的强度赋初值。Weibull分布天然适合描述脆性材料内部缺陷强度分布，在PFC岩石模拟里几乎成了默认选项。形状参数m控制着离散程度：m越大，颗粒强度越集中，模型越接近均匀；m越小，内部强弱差异越明显，破裂越容易从弱颗粒位置起裂。做单轴压缩试算时，m取3～5是常见区间。
下面这段是PFC Python接口的脚本逻辑，用来给所有平行粘结赋Weibull型抗拉强度：
# PFC Python 脚本思路：按 Weibull 分布给平行粘结赋抗拉强度 import random import math def weibull_value(m, k): # m 为形状参数，k 为尺度参数，用于把均值标定到目标强度 u = random.random() while u == 0.0: u = random.random() return k * (-math.log(1.0 - u)) ** (1.0 / m) # 遍历模型中的所有接触（这里用 comp.find(1) 示意查询全部接触） for cb in comp.find(1): # 目标抗拉强度设为 6 MPa，形状参数取 3.0 target_strength = 6.0e6 pb_ten = target_strength * weibull_value(3.0, 1.0) cb.set_prop('pb_ten', pb_ten)逻辑说明：每个接触独立采样一次，这样粘结强度在空间上随机分布，不会出现区域性聚类偏差。set_prop('pb_ten', value)是把抗拉强度写进平行粘结的属性里，PFC的破裂判据会直接用这个值和其他力状态做比较。形状参数m是影响非均质程度的直接开关，尺度参数k在这里设为1.0，再乘上目标强度，实际含义是把分布的中位数/众数拉到目标量级。如果你的岩石单轴抗压强度是40 MPa，目标强度调成40e6 * weibull_value(m,1)也可以，但注意平行粘结强度与宏观强度的换算系数并不是1:1，要先用一套参数标定。
除了pb_ten，还要同步考虑粘结刚度比和法向/切向刚度。常见问题是只随机了强度，刚度保持均匀，这样模型内部应力分布并不自然。我一般会把pb_kn和pb_ks也按同样的Weibull样本做相关性缩放，让弱颗粒的刚度也低一点，破裂才会沿着弱面发展。
2.2 声发射事件记录：以裂纹生成事件作为离散AE源
PFC里没有直接的"声发射传感器"，但原理上可以拿到：声发射本质是微裂纹快速形成时释放的弹性波，而PFC中每一条平行粘结破裂，都会造成局部接触力重新分布，宏观上就是一次微破裂事件。所以最自然的做法是用"bond break事件"代表一个AE事件。
需要在循环里追踪bond断裂。常见做法是用Fish回调或者Python事件接口，也可以在每个加载步后对比当前接触状态列表。下面是一个简化的记录逻辑：
# 在每个 cycle 之后执行，记录新增的 bond break old_break_ids = set(current_break_ids) # 遍历所有平行粘结，检查 bond_state 是否从 unbroken 变为 broken for cb in comp.find(1): if cb.valid() and cb.get_prop('bond_state') == 'broken': cid = cb.id() if cid not in old_break_ids: # 记录该事件发生的 cycle、位置和破裂类型 x, y, z = cb.get_bond_center() ae_events.append({ 'cycle': cycle_counter, 'pos': (x, y, z), 'type': cb.get_bond_type() # 'tensile' 或 'shear' }) old_break_ids.add(cid)逻辑说明：这一步的核心是“算差集”。如果没有在加载前记录初始接触状态，后期tracking会漏掉第一波破裂。bond_state属性在PFC里会随着接触力/强度关系变化，判断它为broken是常见模式。记录的cycle值可以和应力应变history对齐，这样就能画出AE计数率随应变变化的曲线。
参数说明：AE事件的类型区分很重要。拉伸破裂和剪切破裂的比例可以用来判定试样破坏模式——单轴压缩下脆性岩石通常以拉伸破裂为主。如果你的代码里没有区分破裂类型，建议后续补上，因为声发射b值等统计量需要按破裂类型分开统计。
2.3 伺服压实与平衡：让初始模型真实可压
直接生成一堆颗粒并不等于一个岩样。颗粒间空隙大、初始接触少，直接加载会导致压密阶段曲线失真。所以正式加载前需要对试样做固结压实，通常用墙体伺服控制上下或四周边界，给试样一个预定的初始应力。
单轴压缩试验中，如果试样是先用模具围压固结，再拆掉侧向墙体，那么在固结阶段四周围墙要伺服到目标围压，然后移出侧墙。如果只有上下压盘，也要先伺服到一个很小的预压力，让上下压盘与试样完全接触。代码逻辑如下：
def servo_wall(wall_id, target_stress, gain=0.5): # 测量该墙面的平均应力 stress_now = measure_wall_stress(wall_id) stress_diff = target_stress - stress_now # 调整墙体速度使应力趋近目标值 wall_set_vel(wall_id, gain * stress_diff)逻辑说明：这本质是一个比例控制器。gain系数决定收敛速度，过大容易震荡，过小则压实效率低。目标应力设置为0.1 MPa到1 MPa量级即可，目的是让颗粒骨架自锁。伺服压实阶段要持续到模型不平衡力比小于一个阈值，比如1e-5，再进入卸载墙体阶段。我遇到过不少翻车案例，都是因为"平衡没跑到位"就直接加载，结果初始应力场不均匀，裂纹刚加载就在角落乱爆。所以该阶段务必用model.solve(1e-5)或者model.cycle配合unbal判断收敛。
注意：这里的平衡标准和伺服目标需要写进脚本开头，最好做成全局变量，方便换参数。后面章节的加载阶段还要复用这套墙面逻辑。
3. 单轴压缩加载与数据输出：命令流、云图导出和裂纹数截图
模型固结完成后，就进入了真正的单轴压缩阶段。这个阶段最容易出问题的不是加载本身，而是"想监控的数据没有按时保存"。PFC模拟跑一次可能几十分钟甚至几小时，退出后发现应力应变history忘了记录，裂纹图也没截，后悔药是没有的。所以这一章我会从加载方式、数据导出、触发截图三个维度展开，做到每个cycle都能留下东西。
3.1 加载方式选择：墙体伺服与速度加载的取舍
单轴压缩有两种常见加载方式：一种是对上下墙施加恒定速度，另一种是伺服控制轴向应力。岩石力学试验机多用位移控制，仿真里也推荐速度加载，因为速度加载得到的应力应变曲线更稳定，不容易出现应力控制下的急速失稳。
速度加载的代码逻辑很简单：
# 设定加载速率：底部固定，顶部向下移动 bottom_wall.set_vel(0, 0, 0) top_wall.set_vel(0, 0, -0.05) # 负号代表压缩方向 # 第一个 cycle 前记录试样高度作为初始高度 initial_height = top_wall.pos_z() - bottom_wall.pos_z() initial_length = top_wall.pos_z() - bottom_wall.pos_z()参数说明：这里的-0.05表示每秒压缩0.05单位长度，具体数值要根据模型尺寸和颗粒刚度缩放。一个经验法则是让加载速率小到试样内最大不平衡力比一直低于1e-4，也就是准静态条件。如果跑得太快，颗粒会因为惯性飞出去，应力曲线会剧烈抖动，甚至出现"先裂后压"的反常识现象。
有些情况下你想做恒定轴向应力加载，模拟蠕变或者应力路径，那就需要把顶部墙改成伺服控制，目标应力随时间变化。但单轴抗压强度试验一般不需要，直接用速度加载到试样破坏后继续压到残余强度段，再把墙体速度置零。
3.2 应力应变与位移云图数据提取
PFC里提取轴向应力不能只看墙体的力，因为试样截面在破坏时会改变，尤其侧向无围压时试样会膨胀。常见做法是用测量圆或者统计试样中间截面处的接触力来计算应力。简化逻辑：
# 计算轴向应力：统计所有与顶墙接触的颗粒，取接触力总和再除以截面积 contact_force_sum = 0.0 for pe in contact_list: contact_force_sum += pe.force_comp_z() area = math.pi * (radius ** 2) # 圆柱试样截面积，按初始几何计算 axial_stress = contact_force_sum / area # 应变 = 轴向位移 / 初始高度 axial_displacement = -top_wall.pos_z() - (-initial_height) axial_strain = axial_displacement / initial_height逻辑说明：这里截面积取初始面积是基于工程应变定义。如果你想要真实应力（瞬时截面积），需要额外测量侧向膨胀，但在PFC单轴尺度下，初始面积算出的应力曲线已经足够用于参数标定。
位移云图数据导出是另一个重点。PFC本身可以画矢量云图，但要输出可后处理的数据，最好把每个颗粒的ID、位置、位移分量写进文件。代码：
with open('disp_{0:06d}.csv'.format(step_number), 'w') as f: f.write('id,pos_x,pos_y,pos_z,disp_x,disp_y,disp_z\n') for bp in comp.find(1): pos = bp.get_pos() disp = bp.get_disp() f.write('{0},{1:.6f},{2:.6f},{3:.6f},{4:.6f},{5:.6f},{6:.6f}\n'.format( bp.id(), pos[0], pos[1], pos[2], disp[0], disp[1], disp[2]))注意PFC的get_disp()返回的是从某个参考位置累积的位移。如果模型在固结阶段已经发生过平移，那加载前的位移要清零或者减去初始值。我的习惯是压实后、正式加载前，循环一遍所有颗粒，把bp.set_disp((0,0,0))，这样后面导出的位移就是干净的单轴加载响应。
应力云图数据不能直接通过颗粒导出，需要插值。PFC里常用测量圆网格，把每个测量圆内的平均应力张量输出，再固定到空间网格上。如果你只需要后处理画云图，导出颗粒接触力分量然后自己插值也可以。但最省事的做法是直接用PFC内置的plot背景应力云图，再配合plot.store()保存成图片。数据文件是给二次绘图用的，两者都要留。
3.3 裂纹数截图：按加载步骤自动抓取
标题里写的"根据裂纹数截图"，其实就是当裂纹总数达到某个阈值时自动生成截图和数据文件，而不是每隔固定时间步截图。因为裂纹在峰值前后增长极快，固定时间步可能前期一张图没有，后期瞬间几十张。按裂纹数触发能保证从初始破裂到贯通，每个阶段都有代表性的画面。
实现逻辑可以挂在每个cycle的尾部回调：
last_crack_threshold = 0 crack_snap_interval = 50 # 每新增50条裂纹保存一次 def check_crack_snapshot(): global last_crack_threshold crack_count_now = get_crack_count() # 判断是否跨过新的阈值 if crack_count_now // crack_snap_interval > last_crack_threshold: last_crack_threshold = crack_count_now // crack_snap_interval # 保存当前模型快照 model.save('snap_{0:06d}.sav'.format(crack_count_now)) # 截取裂纹图并保存 plot.set('cracks') plot.store('crack_{0:06d}.png'.format(crack_count_now)) # 同时导出当前应力位移云图数据 export_disp_csv(crack_count_now) export_stress_cloud(crack_count_now)逻辑说明：crack_count_now // crack_snap_interval取整除，当跨过下一个50倍数时触发。这样第一次是50、第二次100，直到试样完全破坏。model.save保存的是可断点续跑的模型文件，万一后面发现某阶段参数需要调整，还能从中间状态重跑一部分。截图保存路径建议统一放在snapshot/子目录，避免和模型文件混在一起。
需要提醒一个细节：plot.store('xxx.png')执行的是后台渲染，不会弹出窗口干扰计算。如果是在图形界面跑，每次截图都会占用内存，建议每隔几张清理一次，或者干脆在headless模式下用set plot export命令来存。
4. PFC单轴压缩常见问题排查：五个翻车点与对应修复
我从实际模拟项目里踩过的坑里挑几条高频的，每条都按现象、原因、解决三个层面写。有些问题你光看error log根本看不出来，但结果曲线一眼就能判定错了。
4.1 裂纹与声发射计数问题：三条高频翻车记录
现象：加载了很久，应力从零升到峰值，但裂纹数一直是0。而应力曲线在峰后突然跌落到很低，但裂纹数还是0。
原因：平行粘结压根没开启，或者说虽然设置了pb_ten，但接触类型不是线性平行粘结（Linear Parallel Bond），而是线性接触或者单纯摩擦接触。在PFC里，只有平行粘结才有bond_state，而线性接触没有粘接力，不能存储拉压破坏状态。第二种可能是你计数的方式依赖于接触列表变化，但平行粘结断掉后接触可能仍然存在（变成无粘结接触），只是bond_state变了，如果只查接触数量变化自然永远是0。
解决：在模型生成阶段用model.contactmethod设置平行粘结属性，并确认pb_on为1。计数时不要等接触消失，而是轮询每个平行粘结的bond_state，从unbroken变成broken才算新增裂纹。代码里我建议维护一个全局字典，记录上一轮所有接触的bond state，两轮之间做对比，这样最稳。
还有一条低频但诡异的翻车：裂纹数超过1万之后，计数变得特别快，一秒内同时爆出几百条。原因是PFC的裂纹定义包括微裂纹分裂，一条bond break在极限速度下会连锁触发相邻bond失效。这本身符合物理，但如果你按固定阈值截图，会发现同一时刻保存的多张图几乎一样。解决方法是把触发条件从"裂纹总数阈值"改成"每个加载步内新增裂纹数阈值"，例如新增超过50条时才截图，能有效控制图片密度。
4.2 加载与云图输出问题：从速率到坐标的排查
现象：应力应变曲线前半段抖动剧烈，噪声比信号还大，曲线像锯齿一样。位移云图导出的数据在试样底部出现"负位移"或者边界处明显不连续。
原因：加载速率过快。PFC是显式求解，单步时间内颗粒位移不能太大，否则惯性力主导。特别是非均质模型里某些弱粘结颗粒会先发生局部破裂，释放的应变能转化为动能，造成整体应力波动。位移云图不连续则多半是位移参考点未清零，试样在固结阶段已经发生过整体沉降，直接沿用原始位置导致底部颗粒位移为负。
解决：把加载速率降到原先的1/10，然后看最大不平衡力比是否持续低于1e-4。云图问题更简单，加载前循环所有颗粒把disp归零，再保存一份初始位置表。导出时用current_pos - initial_pos作为位移。另外，应力云图如果有方向性错位，检查是否把局部坐标系的应力分量当作全局坐标导出，PFC测量圆返回的应力张量是全局坐标，不需要额外旋转。
5. 把模拟结果用到位：裂纹数触发截图与云图数据的自动化整理技巧
这一章我给一个我一直在用的整合脚本思路。它是前面所有代码段的组合，但做了一件关键优化：不是每个cycle都检查截图条件，而是每运行N步后检查一次，并且把检查代码放在一个独立的回调里，这样修改加载速率或试样尺寸时不需要改动核心逻辑。
# 全局配置区 crack_interval = 50 check_step = 10 # 每10个cycle检查一次 last_cutoff = 0 def crack_callback(): global last_cutoff if cycle_count % check_step != 0: return crack_now = get_crack_count() cutoff = crack_now // crack_interval if cutoff > last_cutoff: last_cutoff = cutoff snapshot_and_export(crack_now) # 在加载循环开头注册回调 model.callback(crack_callback, 'cycle_start')这段代码把截图策略和加载逻辑解耦。check_step设为10，意味着最多滞后10个cycle才响应，对宏观破坏过程来说足够及时。crack_interval是每50条裂纹出一次结果，你可以按试样总裂纹预计数量调整。比如你预期峰值时裂纹5000条，想在整个过程中截30张图，那crack_interval=5000//30≈166。
数据整理方面，我会在每次截图后把当时的应力、应变、裂纹总数同时写进一个summary.csv，这样后续画图时能直接关联：
def snapshot_and_export(crack_count): # 已有函数，补充一行写入 summary append_summary(cycle_count, axial_stress, axial_strain, crack_count)这么做的好处是，你可以把多次模拟的summary合并，直接对比不同非均质m值下的峰值强度和残余强度。云图数据我一般保存成CSV和VTK两个版本，CSV给自己python脚本做自定义配色，VTK给可视化工具做三维渲染。如果你要写论文，VTK导出的云图边界更平滑，但体积大；不需要渲染的话，CSV足够。
最后说一个我自己的教训：早期做模拟时，我图省事只保存了模型末状态，结果想回溯峰值前后的应力分布，只能重新跑一遍。从那以后我每次做单轴压缩模拟，都会强制走一遍非均质赋值、AE记录、裂纹数触发导出的完整流程，哪怕只是一个小算例也留着中间快照。因为参数标定阶段的试算，每一次结果都可能成为下一组参数的参照，没有中间数据就等于白跑。希望这套流程和避坑记录能帮你在PFC单轴压缩上少走几趟弯路。
本文还有配套的精品资源，点击获取
