简介PFC单轴压缩离散元仿真代码资源面向岩石力学、土木工程与材料科学领域研究者聚焦非均质模型下单轴压缩的声发射分析与破坏机制模拟支持按裂纹数量变化实时截取图像并同步输出应力云图与位移云图数据便于直观追踪材料内部裂纹萌生、扩展直至贯通的全过程。压缩包内含11个文件以docx文档为核心辅以html说明页面与jpg成果图整体仅3.5MB文档系统梳理了代码使用方法、非均质模型构建思路、声发射特征提取方法以及应力/位移云图数据的解读路径结构清晰便于对照查阅。资源已有255人学习适合需要开展PFC二次开发、岩石破裂过程可视化分析或学位论文仿真建模的中高级用户可直接用于单轴压缩数值试验的快速落地与结果输出。1. PFC单轴压缩代码非均质模型、声发射与云图输出的联动逻辑很多人在PFC里做单轴压缩最烦的不是建模型而是三个数据对不上号声发射事件数、裂纹扩展图、应力峰值时刻总是三个时间轴。这套PFC单轴压缩代码把颗粒粘结强度按Weibull分布做了非均质赋参让破坏形态更接近真实岩石同时在主循环里统计裂纹增量一旦累计裂纹数达到预设阈值就自动截图顺手把应力云图数据和位移云图数据落成CSV。对刚上手PFC的研究生它能省掉一两周的排错时间对正在跑批量参数的人它提供了一套可以复用的输出框架。前提是把下面几章的伺服控制、阈值判断和接口差异看明白再往自己的模型里搬。2. 非均质模型构建Weibull强度分配与平行粘结参数映射2.1 默认均质模型为什么压不出真实破坏形态PFC里模拟岩石最常用的是线性平行粘结模型给每个接触一个固定拉伸强度和黏结力。均质模型的单轴压缩典型结果是线弹性段结束后突然形成一条贯穿剪切带应力跌落几乎是垂直的。室内岩石试验常见的峰前非线性、峰后软化、多分支裂隙在均质模型里基本看不到。破坏形态也稳定得不像真实材料每次都从几何上最薄弱的位置启动随机性只来自颗粒排列本身。非均质模型的思路是把强度参数做成空间随机场。最常见的做法是让每个平行粘结的拉伸强度、黏结力服从Weibull分布用均质度系数m控制分布宽窄。m越大强度越集中模型越脆m越小强度差异越大微裂纹在加载早期就开始弥散出现宏观曲线呈现渐进破坏。岩石类材料模拟里m取3到6比较常见具体取值取决于你要复现的岩性。有一个容易被忽略的坑强度随机之后粘结刚度最好不要随机。否则模型在初始平衡阶段就会出现局部应力集中产生伪裂纹。多数代码库的做法是只随机强度、刚度保持统一值这套代码也是这么处理的。2.2 代码给每个平行粘结批量赋Weibull强度下面这段基于PFC的Python接口遍历所有线性平行粘结接触用逆变换采样给每个接触赋不同的拉伸强度import itasca import numpy as np # 建立模型、生成颗粒和墙体默认线性平行粘结模型 itasca.command( model new model large-strain on domain extent -0.012 0.012 -0.012 0.024 ball distribute porosity 0.13 radius 0.3e-3 0.5e-3 box -0.01 0.01 0.0 0.02 wall create id 1 box -0.012 0.012 -0.012 0.024 contact model linearpbond property \ pb_ten 15e6 pb_coh 9e6 pb_emod 10e9 \ kn 1e8 ks 1e8 fric 0.5 ball attribute density 2650 cycle 1000 calm 10 solve steady-state ) # 非均质赋参 m 4.0 # 均质度系数越小越离散 sigma_0 15.0e6 # 特征强度单位Pa rng np.random.default_rng(42) # 固定随机种子便于复现 for c in itasca.contact.list(): cm c.get_contact_model() if cm.get_name() linearpbond: u rng.uniform(0.0, 1.0) # 均匀随机数 # Weibull逆变换采样sigma sigma_0 * (-ln(1-u))^(1/m) s sigma_0 * (-np.log(1.0 - u)) ** (1.0 / m) c.set_prop(pb_ten, s) # 拉伸强度 c.set_prop(pb_coh, s * 0.6) # 黏聚力按强度0.6倍给逻辑说明u来自0到1的均匀分布经过逆变换采样后得到的s服从Weibull分布。pb_ten是平行粘结的拉伸强度pb_coh是黏聚力两者比值取0.6是岩石类材料模拟里的常见区间实际可在0.5到0.7之间调整。粘结刚度pb_emod没有参与随机目的是避免弹性阶段出现非物理的空间刚度波动。参数说明m4.0时模型属于中等非均质破坏区能观察到明显的分支裂隙m降到2.0时破坏会非常弥散峰值强度也会明显降低m超过8.0后基本退化成均质模型。随机种子固定为42配合颗粒生成种子固定可以保证每次重建模型拿到同一套颗粒排列和强度分布方便对比不同加载方案。2.3 伺服加载与初始平衡的执行顺序非均质赋参的前后顺序直接影响模拟是否干净。常见做法是先生成颗粒、建立墙体做初始稳定让颗粒间重叠产生的接触力充分消散稳定之后再对接触赋Weibull强度最后才进入伺服加载阶段。顺序颠倒的话强度已经随机化初始平衡阶段可能就出现第一根裂纹。# 初始平衡不加外载只消除模型内部不平衡力 itasca.command( cycle 0 solve steady-state ) # 赋Weibull强度 # ... 上一节代码 ... # 伺服加载上下墙以恒定应变率相向运动 itasca.command( wall servo on wall servo stiffness 1e8 wall attribute velocity-z -0.05e-3 solve cycles 20000 )这里把墙的加载速度控制在每步0.05毫米量级对应应变率约每秒0.1%到0.5%这个量级下惯性效应可以忽略结果接近准静态。伺服开启后墙体速度由伺服算法动态调节但宏观压缩速率仍由目标应变率决定。若试样尺寸改大了速度要按比例降低否则动态效应会把应力曲线拉出毛刺。下面给出一组常用的参数参考表方便直接套用参数常见取值作用均质度系数m3~6控制强度分布离散程度特征强度sigma_010~30 MPa决定宏观峰值强度水平黏聚比pb_coh/pb_ten0.5~0.7影响剪切破坏占比孔隙率0.12~0.15影响弹性模量与破坏形态加载应变率0.1%/s~0.5%/s保证准静态条件这套参数没有绝对正确答案不同粒径分布、试样高度下都要微调。建议先用m8跑一个接近均质的对照组标定特征强度和宏观强度的关系再逐步调m匹配目标应力-应变曲线。3. 声发射模拟实现AE事件判据、裂纹统计与能量近似3.1 用裂纹增量近似AE事件的基本逻辑声发射在PFC里没有直接的物理实体常规做法是把平行粘结断裂当成微破裂信号。每出现一个新裂纹就认为是产生了一次AE事件。这个假设在岩石力学数值模拟里被广泛接受因为粘结断裂对应真实岩石中晶粒尺度的开裂。实际操作要区分两个概念累积裂纹数和AE事件数。累积裂纹数是单调递增的整数主循环里直接读裂纹总数即可。AE事件数一般按固定时间窗统计比如每2000步内新增了多少裂纹这个增量对应声发射计数率。把计数率画成柱状图叠加在应力-应变曲线上能直观看到峰值前后声发射活跃度的变化。判据设置的合理程度直接决定结果有没有说服力。判据太宽松每个微裂纹都算事件峰前曲线会噪点密集判据太严格事件数太少看不出破裂的渐进性。常见做法是引入能量阈值过滤只有粘结断裂释放的能量超过某个门槛才计入AE事件这样能模拟真实监测系统里只有信号幅度够大才触发的行为。3.2 裂纹统计与AE计数率的代码实现主循环里最省事的方式是用差值法每次保存当前裂纹总数下一次进入回调时做减法。下面这段示意代码展示的是标准的差值统计框架import itasca # 全局记录变量 step_window 2000 # 统计窗口每2000步汇总一次 prev_crack 0 ae_count 0 step_since_last 0 # 每次加载循环中调用 def crack_monitor(): global prev_crack, ae_count, step_since_last n itasca.crack.num() # 当前累计裂纹数 delta n - prev_crack # 本步增量 prev_crack n ae_count delta step_since_last 1 # 达到窗口记录一次AE计数率并清零 if step_since_last step_window: add_ae_record(step_since_last, ae_count) step_since_last 0 ae_count 0逻辑说明itasca.crack.num()返回当前模型中的总裂纹数不同PFC版本可能叫法不同但都会有类似的查询函数。prev_crack记录上一次的累计值两者相减就是本步新增裂纹。窗口内累加的ae_count就是这个窗口的AE计数率。参数说明step_window2000对应加载循环步数而不是真实时间。要换算成物理意义上的AE事件率需要结合每步的应变增量比如2000步对应0.01%应变那么AE计数率的横轴就是应变方便和应力曲线对齐。窗口太长会抹平细节太短又到处都是毛刺一般先跑一次试算看曲线平滑度再调整。3.3 能量近似的三种常见做法单纯数裂纹数只反映了AE事件的频率缺少能量维度。真实声发射监测里能量和幅值是判断破裂尺度的重要指标。PFC模拟里做能量近似一般有三种方式第一种是统计断裂粘结储存的弹性应变能。平行粘结在断裂前相当于一个弹性弹簧存储的应变能约等于0.5倍力与位移的乘积。每个粘结断裂时释放的这部分能量被认为是该AE事件的能量。PFC内置的能量追踪可以输出弹性能历史断裂时刻弹性能的下降量就是近似事件能量。第二种是统计接触力突变幅度。粘结断裂瞬间附近颗粒的接触力会重新分配法向力会出现明显阶跃。把阶跃量当作AE幅值好处是实现简单、不需要额外能量计算缺点是和真实AE能量的对应关系比较模糊。第三种是直接用裂纹密度近似能量即单位体积内新增裂纹数乘以平均断裂能。这种方法最粗但在做相对对比、比如不同m值下的AE活跃度排序时足够用。# 用PFC内置机械能追踪记录应变能变化 itasca.command( model energy mechanical history energy strain solve cycles 20000 )跑完后导出弹性能历史和裂纹增量曲线叠加如果弹性能下降的时刻与裂纹突然增多的时刻基本重合说明AE能量近似是自洽的。这段代码不需要额外插件所有版本基本都支持。4. 裂纹数触发截图与应力/位移云图数据联动输出4.1 触发判据处理好裂纹数跳跃问题标题里说根据裂纹数截图同时输出应力和位移云图数据核心是触发的时机控制。一个常见的错误写法是用if crack_num target判断相等但裂纹数在单步内可能从48直接跳到53永远等不到正好等于50的时刻。正确做法是判断已经越过阈值当前裂纹数大于等于目标值且这个目标还没被触发过。import itasca import csv # 预设触发阈值单位条 TARGETS [50, 100, 200, 500, 800] next_target 0 triggered [False] * len(TARGETS) def check_and_trigger(): global next_target n itasca.crack.num() if next_target len(TARGETS): return if n TARGETS[next_target]: tag TARGETS[next_target] export_state(tag) # 导出云图数据 itasca.command(fview screenshot crack_{tag}.png) next_target 1逻辑说明TARGETS数组定义了想观察的裂纹数节点next_target指向当前应当触发但还未触发的目标。因为裂纹数是单调递增的一旦n TARGETS[next_target]说明已经跨过这个节点马上执行导出与截图然后把目标指针后移一位。这样即使单步内增加了七八条裂纹也不会漏触发。参数说明阈值列表按实际破坏阶段设置一般取峰值前、峰值附近、峰后几个节点。例如试样全破坏大约产生600条裂纹阈值设50、100、200、400、600就能覆盖加载全过程。阈值太稀疏会错过微裂纹起裂阶段太密则输出文件过多通常取5到8个节点就够了。4.2 截图与云图数据的导出代码导出部分分两块应力云图数据我用测量圆阵列读取位移云图数据直接遍历颗粒坐标和位移。测量圆的布置密度决定应力云图的分辨率测量圆半径则影响应力平均的平滑程度。def export_state(tag): # ---------- 应力场数据测量圆网格 ---------- xs np.arange(-0.009, 0.010, 0.002) # 水平位置 ys np.arange(0.001, 0.019, 0.002) # 竖向位置 # 布置测量圆半径0.001m读取平均应力分量 itasca.command( measure delete all ) for x in xs: for y in ys: itasca.command( fmeasure create id {(int(x*10000)int(y*10000))} fcenter ({x},{y}) radius 0.001 ) with open(fstress_map_{tag}.csv, w, newline) as f: writer csv.writer(f) writer.writerow([x, y, sxx, syy, sxy]) for m in itasca.measure.list(): # 不同版本PFC读取应力分量的函数名有差异 # 这里按官方FISH函数measure.stress.xx/m.yy/m.xy读取 xx measure_stress(m, xx) yy measure_stress(m, yy) xy measure_stress(m, xy) writer.writerow([m.center_x(), m.center_y(), xx, yy, xy]) # ---------- 位移场数据直接遍历颗粒 ---------- with open(fdisp_map_{tag}.csv, w, newline) as f: writer csv.writer(f) writer.writerow([x, y, ux, uy]) for b in itasca.ball.list(): pos b.get_pos() disp b.get_disp() writer.writerow([pos[0], pos[1], disp[0], disp[1]])逻辑说明应力数据依靠测量圆阵列测量圆相当于在模型里埋了一圈应变片每个测量圆输出一个平均应力张量。测量圆半径0.001m、间距0.002m得到约9乘9的应力场网格足够绘制云图。位移数据直接读每个颗粒的位置和位移矢量粒度更细画出来是颗粒尺度的位移矢量场。参数说明measure delete all每次都清空旧测量圆避免上一轮数据残留在内存里。测量圆半径不能大于网格间距否则相邻测量圆数据高度重叠云图会显得模糊。位移数据文件里颗粒数量就是模型颗粒总数如果模型有5万颗粒每次触发都会生成一个几万行的CSV磁盘空间要提前算好。4.3 输出文件的作用与后续绘图方式拿到的CSV已经包含空间坐标和应力分量、位移分量可以直接满足两类需求应力云图数据用于观察应力集中的迁移规律位移云图数据用于分析破坏带上颗粒的错动方向和量级。常见做法是导入Origin或Python的matplotlib用tricontourf直接基于散点数据画云图不需要重新插值到网格。如果觉得测量圆网格分辨率不够可以加密到0.001m间距但要注意测量圆越小应力波动越大个别圆的应力值会因局部接触状态异常而明显偏离。绘图前先扫一遍数据去掉边界外的点和明显异常的应力值否则色标会被几个极端值拉坏。5. 避坑指南PFC单轴压缩模拟最常见的五个翻车点5.1 伺服控制不收敛试样像橡皮泥一样一直蠕变现象墙体一直缓慢移动但应力却涨不上去试样在不断变形曲线呈现塑性流动加载步数已经远超预期。原因伺服控制器的增益系数和墙体刚度匹配不当。墙面刚度设得太小墙体受到颗粒推力会产生明显位移伺服系统来不及补偿加载速度太大也会让伺服算法在收敛前就积累大量不可逆变形。解决先把加载速度降一个数量级再把伺服刚度提高到颗粒接触刚度的5到10倍。经验做法是让墙体刚度与颗粒平均接触刚度处于同一量级然后开启伺服并跑500步观察墙体速度是否平稳。如果墙体速度来回振荡说明增益过大如果速度持续缓慢漂移说明刚度不足。5.2 初始内应力过大加载前就出现第一根裂纹现象强度赋完之后还没开始加载裂纹数就不为0。声发射曲线从加载起点就有值破坏形态也乱。原因颗粒生成后颗粒间存在大量重叠产生的高接触力在初始平衡阶段没有完全消散。强度随机化之后这些残余接触力可能局部超过粘结强度直接拉断粘结。解决把初始平衡做彻底平衡判定标准用最大不平衡力和平均接触力的比值常见做法是要求比值低于1e-5再进入下一步。强度赋值也可以等到平衡完成后再做顺序不能颠倒。如果平衡后仍然有局部高应力先跑300到500个无阻尼步让应力重新分布然后再赋强度。5.3 Weibull特征强度和宏观峰值强度对不上现象特征强度设了20MPa压出来的峰值强度只有5MPa差距大得离谱。原因特征强度是微观粘结强度不是宏观强度。破坏总是从强度最弱的粘结处启动一旦裂纹连通整体失效时的宏观应力远低于微观平均值这是非均质模型的固有特性不是代码bug。解决先用m8接近均质的模型标定特征强度与宏观强度的换算关系通常宏观强度是特征强度的20%到40%。再逐步调低m观察宏观强度随非均质性增强而下降的幅度。最后根据目标宏观强度反推特征强度。这个过程至少要跑三组不同m的试算批量调参时这个标定环节最省时间。5.4 AE计数过密导致输出文件暴涨现象设置了AE阈值触发截图结果加载到中段时几乎每几百步就触发一次截图文件几十张CSV上百个磁盘很快爆掉。原因裂纹增量在峰值附近非常集中相邻两个阈值之间可能只隔几十步触发频率远高于预期。另一种情况是AE事件判据太宽松把太多微小裂纹都计入事件。解决把触发阈值间隔拉大比如从50的增量改成100或200的增量。同时给AE记录增加能量阈值过滤只有断裂释放能超过门槛的裂纹才计入事件计数。这样触发密度会回到合理范围。如果必须保留每一条裂纹的完整数据可以把CSV改成只输出测量圆汇总数据不要输出颗粒级位移文件。5.5 云图数据输出时刻和截图不一致现象截图显示的裂纹形态和CSV里的应力分布明显对不上一个像峰值前一个像峰后。原因主循环里触发逻辑执行时模型还在继续求解。截图命令需要渲染场景CSV导出需要遍历颗粒两者耗时较长等导出完模型又往前跑了一段实际输出的模型状态已经不是触发瞬间的状态。解决触发后先暂停求解把当前模型状态冻结再执行截图和数据导出。PFC里可以用cycle 0来停止加载循环或者在执行输出前把求解暂停开关打开。输出完成后再恢复求解。这样保证截图、应力CSV、位移CSV三条数据严格对应同一时刻。6. 进阶用法批量扫均质度参数把云图数据变成论文级图表单次模拟只能验证一套参数批量扫描才能把非均质的影响趋势说清楚。以Weibull均质度系数m为例我可以生成一组模型分别取m2、3、4、5、6固定特征强度、随机种子和其他全部参数跑完自动记录峰值强度、峰值应变、AE累计数。批量跑的时候把上一章的触发函数包进一个外层循环每个m值跑两三次取平均可以有效避免单个随机种子带来的偏差。for m_value in [2.0, 3.0, 4.0, 5.0, 6.0]: for seed in [11, 22, 33]: run_single_simulation(m_value, seed) peak_stress, peak_strain read_history_result() ae_total read_ae_total() save_summary(m_value, seed, peak_stress, peak_strain, ae_total)跑完之后把每次触发导出的应力云图CSV汇总起来画m2和m6的对比云图。一个常见现象是m越小峰值应力越低破坏区越分散AE事件数越多m越大峰值应力越高破坏越集中于一条贯通带。趋势吻合的话就能直接写进论文做参数敏感性分析。云图绘制方面由于测量圆导出的CSV是散点数据matplotlib里可以直接用tricontourf免插值绘图import matplotlib.pyplot as plt import pandas as pd df pd.read_csv(stress_map_200.csv) plt.figure(figsize(4, 8)) plt.tricontourf(df[x], df[y], df[sxx], levels20, cmapRdBu_r) plt.colorbar(labelsxx (Pa)) plt.axis(equal) plt.savefig(sxx_200.png, dpi300)画图前一定要做数据清洗。测量圆边缘可能超出试样边界读取到很小的接触数应力值会异常。把颗粒数少于某个阈值的测量圆剔除再检查应力值是否在合理区间异常点直接删掉。这步数据剪枝虽然听起来麻烦但能避免云图色标被极端值拉坏。我后来批量处理时通常先用pandas快速过滤一遍保留x、y在试样有效范围内的行再去掉三倍标准差以外的应力值。批量参数跑完后的整理工作同样有规律可循先把每个m值下的峰值强度、峰值应变、AE总数汇总成一个总表再单独挑出典型裂纹数对应的云图数据做细节分析。汇总表决定趋势判断细节图决定最终呈现质量两者分工明确。那以后我每次跑非均质单轴压缩模拟都会强制先跑一个近均质的对照组标定特征强度再进入参数扫描而且每个触发节点都检查一遍截图和CSV时刻是否严格对齐。这个习惯让我少丢掉好几批无效数据。希望帮到你。本文还有配套的精品资源点击获取
