几个月前我在跑一组碎石单轴压缩模拟用的还是PFC5.0里最顺手的clump。跑了三天出来一条特别光滑的应力-应变曲线当时还挺开心。结果对照组一放进去我傻眼了试验曲线在峰值之后有明显的跌落和应力平台碎屑越来越多级配一直在变而clump模型从头到尾一颗没碎整个试样像一块铁疙瘩被压变形的根本不是“碎石压碎”。那一刻我才意识到要用PFC5.0模拟碎石、矿渣混凝土这类材料真正的主角不是clump而是可以破碎的cluster底层则是无数个ball通过粘结形成的可破裂结构。这篇东西写给谁看主要给两类人一是做颗粒材料离散元模拟的研究生和工程师模型里需要体现颗粒破碎但不知道cluster怎么建、参数怎么给二是已经建过cluster但发现“不好用”的朋友——要么不碎要么一加载就炸要么碎片数量离谱。后面这些坑我都会讲到。文章里的命令和脚本都是PFC5.0的写法部分细节不同小版本有差异但建模逻辑和排错思路是通用的。1. 为什么非要用clusterclump和刚性颗粒模拟不了的事1.1 从一次失败的clump模拟说起PFC里clump的基本思想是让若干个ball重叠成一个刚性体内部接触不再参与计算只保留颗粒外轮廓和外界接触。它的优势是能快速生成任意不规则形状的“固体颗粒”计算量也比同等形状的ball集合小很多。问题是clump内部自由度为零任何荷载下都不可能发生颗粒内部破裂。这在模拟卵石、钢球、完整岩块这类高强不易碎材料时没问题。但碎石和矿渣混凝土的宏观力学行为恰恰和“内部破裂”强相关。以级配碎石基层为例振动压实过程中粗颗粒尖端会产生应力集中单颗粒可能在接触点附近被压碎或剪碎原来的一个粗颗粒变成几个小颗粒级配随之变细。矿渣混凝土在单轴压缩或劈裂试验中骨料与砂浆界面的裂纹会扩展到骨料内部产生穿晶断裂也就是骨料本身碎了。这些现象用clump完全看不到。我见过不少文章把clump当作“可破碎颗粒”来用其实是在移花接木。他们手动设定一些弱面让clump在特定位置分离——但这只是人为预置了断裂路径颗粒内部并没有物理意义上的渐进损伤和碎片生成。真正要模拟“这颗粒子在受力过程中自己决定在哪断、断成几块”PFC5.0的cluster是更合适的选择。1.2 clump、cluster和散体ball的定位差别三者定位我整理成一个表新手可以按这个表快速判断自己该用哪个对象内部结构是否可破碎内部接触参与计算适用场景单个ball无内部结构否无砂土、粉体、流动问题clump / rigid block多个ball刚性叠加内部自由度冻结否不参与卵石、砾石、完整岩块、边界设备cluster多个ball由bond连接可相对位移是参与碎石、矿渣混凝土、软岩、煤岩一个容易混淆的点是PFC5.0中ball是几何单元clump和cluster都是ball的集合体但内部力学完全不同。clump内部接触被冻结计算效率高cluster内部每一个pebble之间的接触都参与力学计算并且额外承担bond强度。换句话说cluster是把“一堆ball”用粘结剂粘起来外力够大时粘结先破坏然后颗粒组解体。1.3 什么时候必须上cluster我判断一个项目是否必须用cluster就看三个问题模型是否要输出级配演化也就是“破碎后各粒径质量占比”如果要必须上cluster因为只有bond断裂后的碎片才有天然的粒径分布。是否关注峰值后的残余强度和剪胀行为脆性颗粒材料的脆-延转化通常由内部破碎主导clump给出的残余平台往往是假的。是否关注局部应变集中和断裂网络比如碎裂带、剪切带怎么扩展cluster能让你直接看到bond断裂的空间分布这是试验里都很难直接观测的。如果你的模拟只是看整体应力应变趋势、不管颗粒破碎那用clump反而效率高。但如果“破碎”本身就是研究对象就不要绕路直接面对cluster的建模成本和参数标定问题值得。2. 破碎的内核cluster的细观结构和bond断裂逻辑2.1 一个cluster到底由什么组成cluster在PFC5.0里的本质是一组互相接触的pebble本质上就是ballpebble之间通过粘结模型连接。注意“粘结”和“接触”是两个概念接触是几何上的贴合只要有重叠就会产生排斥力粘结是额外的力学连接在接触点额外传递拉力和弯矩。我把cluster想成砖块和水泥的关系。pebble是砖块bond是水泥。没有bond的ball堆在一起就是一垛没砌的砖一碰就塌有了bond它们才能形成一块完整的“人造石”但它不像clump那样水泥终身坚固而是受力超限时水泥接缝会开裂整块石头碎裂成若干小块。这样的碎片传播过程正好对应真实材料里骨料破碎的细观机制。建模时需要在两个层面选择参数pebble层面密度、刚度kn/ks、摩擦系数fric这些决定了“砖块”本身的变形和摩擦力。bond层面法向刚度pb_kn、切向刚度pb_ks、抗拉强度pb_ten、黏聚力pb_coh、内摩擦角pb_fa这些决定了“水泥”的力学性能。模拟碎石和矿渣混凝土我default接触模型直接用linearpbond平行粘结。有些资料推荐contact bond但我个人不建议用于这类材料contact bond只在接触点传力不传力矩颗粒一旦有轻微转动就会很脆地断开碎片容易过碎平行粘结在接触处形成一个有面积的“胶接盘”可以同时传递力和弯矩更接近岩石、混凝土骨料之间胶结物的力学行为。2.2 平行粘结是怎么断的强度判据与破坏模式平行粘结的核心逻辑是在每一个bond接触面上法向力和切向力共同作用产生法向应力sigma和切向应力tau。当应力超过材料给bond设定的强度bond断裂这个接触从“粘结接触”退化成普通接触。用公式说PFC5.0的linearpbond会实时计算法向应力sigma_n F_n / A其中A是bond的截面积。切向应力tau F_s / A同时考虑扭转和弯矩产生的附加应力。判断条件分两种拉断当法向拉应力sigma_n超过抗拉强度pb_tenbond立刻断裂。剪断当切向应力tau超过剪切强度即pb_coh加上法向压力产生的摩擦项pb_fabond断裂。这里需要特别理解剪切判据tau_c pb_coh sigma_n * tan(pb_fa)。这意味着法向压力越大bond越不容易被剪断——材料内部的摩擦角pb_fa牵制了剪切破坏。这个特点很重要模拟碎石在高压三轴下的“压而不碎”和低围压下“脆性破碎”的分野很大程度上就是pb_fa在起作用。理解了断裂判据就能理解cluster的破坏模式为什么与实际实验对得上。单轴压缩时试样内部主应力是轴向的最危险的是那些垂直于加载方向的bond它们承受拉应力材料表现为劈裂破坏高围压时法向压力抑制了拉断剪切断裂变成了主导材料表现为更韧的破坏。这就是为什么cluster模型能复现很多试验中的破坏模式转变。2.3 为什么“碎片阈值”和pebble数量决定了模拟质量刚上手cluster最容易犯的错误是单颗cluster里的pebble数量太少比如只塞5到6个球。这种情况下bond断裂的模式极其有限cluster要么不碎要么断成两大块根本谈不上级配演化。pebble数量直接决定了一个cluster能够“碎成多少种可能的碎片组合”。数量太少破坏路径几乎没有选择空间数量太多内部接触和bond数量呈几何级数上涨计算成本飙升而且很多bond从一开始就处于高应力的“冗余”状态加载初期的破碎行为会被这些冗余bond干扰。从我自己的经验看模拟碎石或矿渣混凝土骨料单个cluster内pebble数量控制在20到40个比较合理。这个范围下既能保证足够的内部裂缝扩展路径计算量也可以接受。如果你要把几万个cluster放进一个试样那单个cluster的pebble数量还要再往下砍否则模型规模会失控。另外一个细节是pebble的粒径级配。cluster内部的pebble也不能全做成等径球最好有个1.5到2.5的半径比。等径球组成的cluster内部结构过于规整压碎时容易产生规则的贯穿裂纹不真实。适度引入粒径变化能让cluster内部结构更像真实骨料破碎模式会自然散乱许多。3. 从零搭一个能破碎的clusterPFC5.0命令流与FISH脚本3.1 建模流程总览PFC5.0里建cluster的完整流程我拆成五步生成pebble集合在目标几何区域内随机填充ball控制粒径范围和孔隙率。定义cluster模板告诉程序哪些ball属于同一个cluster以及该cluster在空间中的包络范围。实例化cluster在指定坐标位置生成一个cluster程序会把该位置区域内的ball组合成cluster。设置bond参数给cluster内部的pebble之间的接触赋予平行粘结属性。建立加载环境加上伺服墙或加载板进入求解循环。需要注意PFC5.0的cluster模板功能在不同小版本里位置有差异。我下面的脚本是常用版本写法如果你的版本提示不识别cluster template或cluster create去菜单里找一下rblock或cluster相关模块另一条路是用rblock创建模板再生成cluster逻辑是一样的。3.2 关键命令与脚本示例先看最核心的建模段; PFC5.0 可破碎cluster建模示例 model new ; ---- 1 设定全局参数 ---- def par_cluster global r_min 1.5e-3 global r_max 2.5e-3 global cx 0.0 global cy 0.0 global cz 0.0 end par_cluster ; ---- 2 在正方体空间生成pebble集合 ---- ball distribute box x [-0.04] [0.04] ... y [-0.04] [0.04] ... z [-0.04] [0.04] ... radius [r_min] [r_max] ... tries 20000 ; ---- 3 声明cluster模板 ---- cluster template create id 1 ... type box ... length_x 0.08 ... length_y 0.08 ... length_z 0.08 ; ---- 4 在指定位置生成可破碎cluster ---- cluster create x [cx] y [cy] z [cz] ... id 1 template-id 1这里说明一下参数含义r_min和r_max是pebble半径范围实际颗粒级配就是由这两个值控制的。tries 20000是生成球集合时的最大尝试次数如果初始密实度过高导致重叠球太多程序会在这里反复重试直到满足目标孔隙率。cluster template create定义了一个虚拟的边界形状稍后实例化时落在该边界内的ball会被组合成一个cluster。需要注意ball distribute生成的只是独立的ball它们之间是普通接触还没有bond。cluster的作用是把这些ball在逻辑上归组归组之后才有下一步的bond赋值。3.3 接触模型与bond参数的注入方式cluster生成后需要给cluster内部的pebble之间赋予平行粘结。最简单的做法是设置全局接触模型默认值让后续建立的接触都采用linearpbond; ---- 5 默认接触模型设为平行粘结 ---- contact cmat default model linearpbond ... property pb_kn 3e7 pb_ks 2.5e7 ... pb_ten 2e6 pb_coh 3e6 pb_fa 35 ; ---- 6 设置ball本身的接触属性 ---- contact cmat default property kn 2e7 ks 2e7 fric 0.6在PFC5.0中cmat default设置的是“之后再生成的接触”或“尚未显式设置接触模型的现有接触”的默认属性。所以顺序很重要先生成ball再设置cmat最后跑几个cycle让程序把默认属性写到现有接触上。有时候cluster内部的bond没有被正确继承这时需要主动遍历cluster内的pebble手动给它们之间的接触设置bond强度。下面是一段简单示例; ---- 7 手动检查并补齐cluster内bond属性 ---- fish define repair_bond local cp cluster.find(1) if cp ! null local count cluster.pebble.count(cp) loop local i while i count local p1 cluster.pebble(cp, i) if i count - 1 local p2 cluster.pebble(cp, i 1) local c contact.find(p1, p2) if c ! null contact.prop(c, pb_ten) 2e6 contact.prop(c, pb_coh) 3e6 contact.prop(c, pb_fa) 35 endif endif i 1 endloop endif end repair_bond这只是一个示意脚本真实建模时应遍历cluster内所有pebble对并给所有真实存在的接触赋bond属性。如果你用的是比较新的PFC5.0版本cluster模板本身就带有bond属性字段在定义模板时通过cluster template property一并设置后续生成的cluster直接继承会省事很多。加载环境的建立也提一下。模拟碎石或矿渣混凝土的单轴压缩常见做法是上下两个伺服墙wall generate box ... wall attribute normal_stiffness 1e8 ; 给上板一个向下速度 wall attribute yvelocity -0.005 range id 2加载速率需要刻意控制。如果速度太快惯性效应会把破碎模式带偏出现“试样中央炸裂”而不是“从加载端向内部传播裂纹”。我建议先把应变率控制在低于0.1/s具体计算公式是加载速度 / 试样高度。如果试样高0.1m那加载速度不应超过0.01m/s。4. 让模拟“像真材料”碎石和矿渣混凝土的参数标定流程4.1 标定目标与试验参照cluster模型的参数不能直接查材料手册它是细观参数只能通过宏观试验结果反演标定。标定的目标量一般是四个弹性模量E、泊松比ν、单轴抗压强度σc、残余强度或破坏模式。以级配碎石为例室内单轴压缩和常规三轴试验的典型结果是E在30到60 GPa量级σc在几MPa到几十MPa之间破坏模式表现为明显的劈裂和局部破碎。矿渣混凝土则不同材料内部有砂浆基体和骨料两个相如果只模拟骨料颗粒的破碎标定目标更多是“骨料在混凝土破裂后不立即粉化而是形成一定数量的粗碎片”。建模之前先把手头试验数据整理好。我习惯做一张这样的表物理量目标值对应PFC参数弹性模量E由试验曲线弹性段斜率得到kn、ks、pb_kn、pb_ks泊松比ν试验横向/纵向应变比值kn/ks 或 pb_kn/pb_ks峰值强度σc试验峰值应力pb_ten、pb_coh残余强度 / 剪胀角破坏后的应力水平fric、pb_fa4.2 标定顺序与经验公式标定顺序非常关键。我建议的流程是固定摩擦性质先把ball的fric设成0.5到0.7pb_fa设成30到35度附近这两个值后期尽量少动。调刚度匹配弹性段单轴压缩算到弹性阶段就停观察应力应变曲线的斜率。E偏大就同时降低kn、ks、pb_kn、pb_ksE偏小就同时增大。这一步可以先只调线性刚度不碰bond强度。调泊松比E匹配后看横向应变。如果ν偏大增大kn/ks的比例比如从1.0调到1.5如果ν偏小反过来减。调bond强度匹配峰值峰值强度对pb_ten和pb_coh最敏感。先保持pb_ten与pb_coh比例约1:1.5整体乘一个系数看峰值是否落在目标区间。调破坏模式如果试样表现为整体崩解说明bond强度分布太均匀需要给bond强度加一点随机波动如果加载早期就大量破裂说明初始应力集中太强需要检查初始平衡状态。一个粗略的经验关系峰值抗压强度大约正比于bond抗拉强度系数通常在3到8之间。也就是说如果试验峰值强度是10 MPapb_ten可以先试2 MPa量级然后微调pb_coh。这个线性关系能帮你把标定起点定准而不是靠瞎猜。4.3 一组典型参数参考下面给一组我标定过的碎石和矿渣混凝土参数范围供做初值参考。注意不同级配、不同含水率、不同加载方式都要重新标定下面数值并非万能参数级配碎石矿渣混凝土骨料说明ball密度 kg/m³26002400堆叠密实度会自然影响质量kn / ks (N/m)2e71e7弹性变形主导fric0.60.5控制残余强度pb_kn / pb_ks (Pa/m)3e7 / 3e71.5e7 / 1.5e7bond刚度pb_ten (Pa)2e68e5抗拉强度pb_coh (Pa)3e61.2e6黏聚力pb_fa (°)3530bond内摩擦角第一趟标定时bond强度宁低勿高。强度设置太高cluster完全不破碎算出来的结果是clump的翻版强度稍低至少能观察破碎行为再逐步往上调。4.4 cluster内部的级配与尺寸效应标定过程中我发现一个高频坑试样内的cluster数量太少或者单个cluster太大会导致宏观结果对“该cluster内部bble怎么排”异常敏感。换一组随机种子峰值强度能差20%以上。这不是数值错误而是代表性体积单元不够大。解决办法有两个一是增加试样内cluster的数量让统计平均覆盖单颗cluster的随机性二是在同一cluster内部让bond强度有适当空间波动比如按FISH随机乘一个0.8到1.2的系数模拟真实骨料内部微缺陷分布。后者对破碎模式的真实感提升非常明显。至于试样尺寸一个直径0.1m、高0.2m的标准圆柱试件颗粒数量至少到数千颗才谈得上稳定。如果完整试件上万颗甚至几十万颗颗粒跑不动建议先缩小到1/4或1/8标定参数确认力学响应稳定后再放大模型。5. 调试与排坑cluster不碎、炸裂、碎片异常怎么办5.1 不破碎先查力链还是先查bond参数cluster模型最常见的故障是“明明设了bond为什么跑完整个加载过程一颗都没碎”。第一优先要查的不是bond强度而是力链是否真正传到了cluster内部。很多时候模型里的试样在压缩过程中已经整体刚体运动了上下墙和cluster之间只有一两个接触点大力全部被墙消耗掉cluster内部的bond自然不碎。这种情况在初始堆积松散、孔隙率偏高的模型里特别常见。解决方法就是先把孔隙率降到合理区间比如0.30到0.40然后跑足够的cycle让试样先产生稳定的接触力链。排除了力链问题再检查bond强度。我之前常犯的错误是用宏观抗压强度去估算bond强度结果设高了10倍。记住4.2节里的换算经验bond抗拉强度大约是宏观峰值强度的1/5到1/8。如果你设的pb_ten已经低于这个量级还是不碎那把加载速度调快一点试试——这多半是时间步内变形量太小还没到bond断裂阈值而模型的阻尼又把应力波完全吸收了。5.2 一加载就炸初始均衡与阻尼的坑另一种常见故障是模型“开局即炸”。刚开跑几个cyclecluster就从内部崩开或者整个试样瞬间失稳。这时候问题多半不在bond参数而在初始状态没搞干净。PFC5.0里生成ball时随机填充会造成大量重叠接触产生巨大的初始不平衡力。这种不平衡力如果不消除就赋bond相当于给每个bond施加了额外预应力一旦开始加载这些初始应力会率先把bond推垮。所以在生成cluster之后、赋bond之前最好先跑几百到几千个cycle让系统把不平衡力降到足够小。标准做法是检查model unbalanced-force把它降到最大接触力的1%以下再继续。阻尼也值得单独说。PFC5.0默认的局部阻尼是0.7这个值对大部分静力问题是稳的但如果你的模型动态效应明显比如加载速度太大会出现应力波在试样内反复弹跳导致局部炸裂。把局部阻尼临时调到0.0用model solve跑静力或者降低加载速度都可以缓解。5.3 碎片数异常和性能问题碎片数异常有两个方向一个是“全碎成单颗粒”一个是“几乎不碎但是碎得很碎”。全碎成单颗粒多半是bond没有在cluster内部完全建立只是给少数几个接触赋了bond剩下的pebble本来就没粘住加载时自然全部脱落。检查方法是统计cluster内实际bond数量fish define count_bond local cp cluster.find(1) if cp ! null local n 0 local count cluster.pebble.count(cp) loop local i while i count local p1 cluster.pebble(cp, i) if i count - 1 local p2 cluster.pebble(cp, i 1) local c contact.find(p1, p2) if c ! null if contact.prop(c, pb_ten) 0.0 n n 1 endif endif endif i 1 endloop io.out(cluster 1 bond count: string(n)) endif end count_bond我见过一个案例cluster有200多个pebble实际bond数量却只有30结果一开始加载就全散了。补上bond之后破碎模式立即正常。“碎得很碎”通常说明pebble粒径比太大或bond强度太低。把半径比降到2以内、适当提高pb_ten和pb_coh会让碎片尺寸更集中在几个大块上。性能问题则是另一个烦恼。cluster模型比clump慢一个量级因为内部接触全参与计算。优化思路有三个减小pebble数量、用wall替代不关心的ball边界、设置合理的计算域。还有一个小技巧跑纯压缩标定时把模型尺寸缩小标定完参数再放大到最终模型能省下大量的调试时间。5.4 调试时的可视化检查调试cluster时我强烈建议在PFC5.0里把接触和bond的状态可视化出来plot create view_bond plot add contact type bond plot add contact type crack这样可以直观看到加载过程中哪些位置先裂、裂纹怎么扩展、碎片怎么分离。很多时候光看数值曲线根本看不出问题但stress-strain曲线一出现异常拐点配合裂纹云图马上就能定位原因。6. 数据不会骗人破碎率、配位数、能量曲线与碎片统计6.1 颗粒材料破碎的核心量化指标模型能跑出破碎之后下一个问题是怎么把“破碎了多少”量化出来。最常用的是Marsal破碎率。试验中通过筛分得到加载前后的级配曲线两个曲线之间围成的面积就是破碎率。PFC里实现这个思路有两个入口直接统计bond断裂数量用断裂bond数除以初始bond总数作为“内部损伤度”。这个方法简单但bond断裂并不完全等价于宏观破碎。真正意义上的Marsal破碎率需要把cluster的碎片按“等效粒径”归类做虚拟筛分。碎片等效粒径可以按碎片内pebble体积折算成球径然后统计加载前后级配曲线。第二种结果更容易和试验对比但实现复杂。如果只是做相对对比我个人认为bond断裂数加碎片数量组合已经够用。配位数coordination number是另一个直观指标表示每个ball平均有多少个接触CN 2 * N_c / N_b其中N_c是接触总数N_b是ball总数。配位数下降说明颗粒间接触变少骨架在“溶解”破碎刚开始时配位数往往先下降之后新生成的细小碎片重新接触配位数又有所回升。这个“先降后升”的形态是颗粒破碎的典型特征写论文或者做报告时可以重点展示。6.2 用FISH和history把数据取出来PFC5.0里提取这些指标不难。history命令可以记录想要的标量; 记录不平衡力 history id 1 model unbalanced-force ; 记录平均配位数需要FISH函数支持 history id 2 ball coordination-number ; 记录bond断裂数需要FISH函数维护计数器 history id 3 bond-break-count history interval 100bond断裂数在PFC5.0里可以用crack相关FISH函数统计如果你的版本没有内置也可以自己遍历接触检查bond是否已经失效每次失效就把全局计数器加一。注意history只能记录FISH函数返回的标量值所以必须先把统计逻辑写成FISH函数再被history引用。一个完整的“碎片统计”FISH示例如下fish define frag_stat local btot 0 local bbreak 0 ; 遍历所有bond接触统计总数和断裂数 loop foreach c contact.list if contact.type(c) 8 ; linearpbond类型编号因版本而异 btot btot 1 if contact.prop(c, pb_ten) 0.0 bbreak bbreak 1 endif endif endloop global frag_ratio bbreak / max(btot, 1) global frag_count bbreak end frag_stat注意PFC5.0里bond断裂后接触模型通常会退化为线性接触contact.type的编号会变化我这个例子只写逻辑具体类型编号需要对着你版本的FISH手册查。6.3 从曲线判断破坏阶段的经验处理这些数据时可以把整个过程分为三个阶段第一阶段加载初期配位数略有上升或稳定bond断裂数极少。这是试样压密和力链重组的阶段。第二阶段应力接近峰值时bond断裂数量开始陡增配位数明显下降。这是破碎活跃期对应试样内部裂纹大量扩展。第三阶段峰值后bond断裂增速放缓碎片开始重新排列和接触配位数可能略有回升应力曲线出现残余平台。这三个阶段的划分直接对应试验中的声发射特征。如果有条件把每千cycle的bond断裂数做成“事件率”曲线你会看到它和真实的声发射振铃计数曲线长得非常像这也能侧面验证你的模型抓到的是不是真实的破碎行为。最后分享一个笨办法。我在做碎石压实模拟时最怕的就是参数组合不唯一bond强度和摩擦角都能让峰值强度看起来差不多但残余阶段差很远。所以我的顺序永远是先固定摩擦性质再调bond强度把这个顺序固化到标定流程之后一次标定基本能控制在两到三轮试算内。如果哪次模型跑出来的破碎模式和试验对不上我第一件事是去检查初始孔隙率和初始平衡而不是急着调bond——因为内力链状态不对后面所有参数都白调。