★ FEATURED ARTICLE 仿真做了一周峰值强度始终对不上实验室数据——那种感觉做过混凝土离散元模拟的人应该都懂。后来导师点了一句“你试件里面是圆颗粒还是碎石”我才反应过来用纯ball堆出来的模型本质上模拟的是“玻璃珠混凝土”而不是碎石混凝土。这篇文章把我在PFC5.0里做碎石混凝土单轴压缩模拟的完整思路整理出来重点围绕ball加clump颗粒建模、单轴压缩实验的伺服加载与边界控制以及最容易被忽略的声发射事件数代码的实现逻辑。适合正在做混凝土/岩石颗粒流模拟、被参数标定或者后处理统计卡住的研究生和工程师。1. 碎石混凝土建模的第一步为什么纯ball模型注定失败1.1 圆颗粒的“滚动假象”与骨料嵌锁缺失纯ball颗粒在PFC里是标准圆盘圆盘之间天生就会滚动。单轴压缩时纯ball试件的宏观表现往往是“颗粒旋转整体外鼓”而真实碎石混凝土破坏时粗骨料几乎是不转的裂纹主要在砂浆里扩展遇到骨料就绕行。这个差异直接决定了强度、破坏模式和声发射事件分布的失真。深一层讲这是因为纯ball模型缺少骨料嵌锁效应。碎石骨料是不规则多面体相邻骨料咬合在一起能抵抗滑移宏观上体现为更高的摩擦角和抗剪强度。PFC里要让这种嵌锁效应出现光调摩擦系数和粘结强度是不够的——你加再大的bond强度圆颗粒之间的接触面依然光滑试件在压缩过程中更像一团“胶水粘起来的珠子”。我最早用纯ball做的时候为了让峰值强度够到42 MPa只能把平行粘结的抗拉强度调到很高结果弹性段斜率跟着一起抬应力应变曲线整体“立”起来和实验完全对不上。这就是颗粒形状对宏观响应的系统性影响不是参数可以补偿的。1.2 clump的核心参数pebble数量、重叠量与面积占比clump是由若干个pebble圆形子颗粒刚性合并成的不规则块体内部pebble之间不定义接触整个clump作为一个刚体运动。用clump模拟碎石骨料能抓住“形状嵌锁”这个本质特征。建模时最需要控制的是三个参数pebble数量、重叠程度、骨料面积占比。pebble数量我建议单个碎石取5到10个左右。太少形状接近椭圆嵌锁效果弱太多计算量陡增而形状改善有限。我常用的做法是用6到8个半径不等的pebble拼一个凸多边形轮廓生成出来的块体既不规则又不至于太“扎手”。PFC5.0里可以用clump template在生成前先预览形状这一步值得做因为形状直接决定后面裂纹绕行路径是否真实。重叠量控制不好会出现两类问题pebble之间重叠太少clump轮廓凹坑太多骨料形状失真重叠太多clump等效密实度高和砂浆基体接触的周长偏大导致断裂行为整体偏“韧”。我的经验是相邻pebble重叠半径控制在1.3到1.5倍小pebble半径区间生成后先统计一下每个clump的表面积和等效粒径做一个分布直方图看看。面积占比这块正常碎石混凝土粗骨料体积率在40%左右二维模拟里可以近似用“骨料面积/试件总面积”等于这个比例。但要注意PFC里生成的clump内部是实心的所以投放时要预留砂浆厚度实际投放率可以比目标值低2到3个百分点否则后期平衡时挤压过大初始应力场就乱了。1.3 级配投放从筛分曲线到随机位置的完整流程级配投放流程比我第一次做时想象的要啰嗦核心是“先生成模板再随机放置再平衡”。第一步是根据筛分曲线确定粒径范围。比如5到10 mm、10到15 mm两档碎石先把每档的中值粒径换算成clump等效半径再按照筛分累计通过率分配每档的数量比例。这一步用EXCEL算完直接生成一个list没必要在PFC里再算一遍。第二步是随机位置生成。PFC5.0里可以用FISH写一个循环每次从模板库里随机抽一个clump给定随机的x、y坐标和旋转角放置前先检查与已生成clump的距离是否满足最小间距要求。最小间距一般取相邻颗粒半径之和的1.1倍太近后面平衡容易产生飞溅太远面积占比又达不到。第三步是给骨料和砂浆之间加“界面过渡区”。这是很多人忽略但非常关键的一步。真实混凝土的骨料-砂浆界面ITZ是薄弱区强度只有砂浆本体的60%到80%。在PFC里怎么模拟把clump上的pebble与周围ball之间的平行粘结参数设定为砂浆接触的0.6倍左右即可。这样压缩时裂纹会优先沿骨料界面扩展破坏模式一下子就贴近真实了。我这边用这个设置之后声发射事件的空间分布明显开始沿着骨料边界聚集这是纯ball模型完全给不出来的现象。2. 细观参数标定把宏观实验翻译成PFC能认的参数2.1 接触模型与粘结模型的选择逻辑混凝土在PFC里最常见的组合是“线性接触 平行粘结模型”linear parallel bond model。为什么不用点接触粘结contact bond因为点接触只能在接触点传力不能传递弯矩模拟出来的试件一旦出现局部破坏裂纹就像撕布一样快速贯通强度偏低、脆性偏强。而平行粘结在接触点附近覆盖一个有限尺寸的“胶层”既能传力又能传弯矩抗弯抗剪能力都更接近砂浆的胶结行为。这个模型里的参数分两组一组是接触本身的包括法向刚度、切向刚度、摩擦系数另一组是平行粘结的包括粘结有效模量、法向/切向强度、粘聚力和内摩擦角。设置时先分清楚哪些参数控制“弹道阶段”哪些控制“破坏阶段”标定才不会越调越乱。需要提醒的是clump内部的pebble之间不需要设接触和粘结因为clump是整体运动。而clump与clump之间、clump与砂浆ball之间需要设接触和粘结其中钢筋类的模拟还要考虑界面但碎石混凝土模拟这里不涉及。2.2 五步标定流程从弹模到破坏模式我一直用的标定顺序是从“弹性”到“强度”再到“破坏形态”一次只动一个参数第一步调平行粘结的有效模量。单轴压缩应力应变曲线的初始斜率主要由这个参数决定。PFC里的有效模量不是直接等于宏观弹模但它和宏观弹模是近似线性关系先试算一组值做一条EMOD-斜率标定线再找目标值。第二步调刚度比匹配泊松比。法向刚度与切向刚度的比值主要影响泊松比通常在1.0到4.0之间。混凝土的泊松比一般在0.2附近如果你算出来0.35说明刚度比太大或者加载端约束太强。第三步调平行粘结的抗拉强度和粘聚力。这两个决定峰值抗压强度。经验上抗压强度对粘聚力更敏感抗拉强度对接下来的劈裂破坏形态更敏感。先锁定粘聚力把峰值强度顶到差不多的高度再微调抗拉强度让应力峰后曲线不至于太脆或太韧。第四步调摩擦系数匹配残余强度。峰后曲线掉下来之后试件靠骨料和破碎颗粒之间的摩擦继续承担荷载。摩擦系数越大残余强度越高、峰后曲线越平。第五步检查破坏模式。把模拟的最终破裂图与实验室试件破坏照片对比看是斜向剪切带、竖向劈裂还是“X”型共轭破坏。不对就回到第三步小幅度调整粘聚力和抗拉强度的比值。2.3 参数标定中两个最容易翻车的变量翻车点一局部阻尼。PFC默认局部阻尼给到0.7的爱好者不少但那是为了快速稳定。单轴压缩时阻尼太高相当于给试件“泡在水里加载”应力波传播被抑制峰值强度和破坏形态都会变得乱糟糟。我建议标定时把局部阻尼固定在0.1到0.2之间并且在整个模拟过程中不要变。翻车点二刚度比的“方向性”。有些教程会把法向刚度设得远大于切向刚度来模拟材料不可压缩性但在混凝土单轴压缩里刚度比过大会导致压缩过程中侧向变形明显偏小泊松比偏低裂纹扩展以张拉为主最后变成竖直劈裂而不是剪切带。记住刚度比的标定必须和破坏形态一起看不能只看应力应变曲线的拟合程度。3. 单轴压缩模拟的实施细节伺服加载、数据采集与应力应变3.1 试件尺寸与加载速率的匹配原则试件尺寸方面二维模拟里我习惯用宽50 mm、高100 mm的标准比例保证高径比2:1这是和实验室试件一致的前提。clump粒径分布在5到10 mm时整个试件大概能放几十颗骨料ball数量控制在几万量级计算效率比较友好。加载速率的设置是许多初学者的盲区。PFC里墙的加载速度不能直接照抄实验室万能试验机的2 mm/min因为数值模拟没有“真实的加载时间”概念时间步长、阻尼、刚度会影响准静态效果。我做加载速率敏感性测试的经验先把墙速设为0.1 mm/s跑完一次再把墙速降到0.02 mm/s跑一次如果两次峰值强度差异不超过2%说明这个速率已经可以接受。如果差异在5%以上说明存在明显惯性效应必须继续降速。有一点很多人没意识到加载速率对声发射事件数的影响极大。高速加载时大量bond在极短时间窗内同时断裂事件数被“压缩”成一两个大事件低速加载时微破裂逐条出现事件分布更接近实验室结果。所以做声发射统计之前先做一个速率敏感性测试确定一个稳定窗口再跑正式模拟。3.2 伺服控制墙的逻辑与初始应力平衡单轴压缩实验的边界条件是底部固定、顶部加载、四周自由。PFC里最简单的做法是底部放一面固定墙顶部放一面加载墙加载方式有两种指定墙速度或者伺服控制力。要模拟准静态压缩我更推荐伺服控制因为固定速度加载在试件进入软化段的时候容易出现“局部垮塌式”卸载再加载的抖动。伺服控制的本质是每个计算步检测当前墙受到的应力与目标应力做差通过误差量调整墙的移动速度让应力加载过程平稳。增益系数的选择很关键增益太大墙会振荡太小加载半天爬不上去。PFC经典公式是def servo_control while_stepping local cur_stress wall.stress(wall_id) local err target_stress - cur_stress local wv gain * err wall.vel(wall_id) wv end增益gain怎么定有一个经验公式gain alpha * 试件宽度 / (接触数 * 法向刚度 * 时间步长)alpha取0.3到0.5。这个公式的物理含义是把“当前不平衡力需要多少速度才能在一个时间步内抵消掉”算出来。我通常在标定模型时先跑10000步看墙应力是否平稳逼近目标值不稳就减小alpha。另一个容易出错的地方是初始应力平衡。生成clump和ball之后颗粒之间难免有重叠接触力不均匀直接加载会让试件一开局就“内伤”。我的做法是生成后用较大阻尼0.7先跑一段平衡等最大不平衡力降到平均接触力的1%以下再把局部阻尼调回0.1到0.2开始正式伺服加载。这一步顺序颠倒了后面所有结果都会带“初始内应力”污染。3.3 应力应变的计算口径与历史曲线输出应力和应变的计算口径如果搞混整个模拟就白跑。应力不用测量圆去算直接用加载墙的反力除以试件截面积。二维PFC里默认厚度为1所以截面积就是墙的接触总力除以试件宽度。PFC里可以用wall.force.contact(wall_id)取总力再除以宽度。应变用墙的累积位移除以试件初始高度。我的历史数据记录习惯是每100步输出一次轴向应变、轴向应力、当前断裂总数、声发射事件累计数、最大不平衡力。PFC里用history命令可以记录但我更推荐在FISH里自定义一个记录函数把每个采样点的数据写入CSV文件因为后面做峰后曲线、声发射事件时序图的时候CSV处理起来比PFC内置history方便得多。需要特别注意的是加载终止条件。试件软化到残余段应力降到峰值的30%以下如果还在继续伺服加载墙可能反向拉动试件导致接触力分布异常。我一般在应力降到峰值20%时停止加载保存最终破裂图和接触力分布再做后续统计。4. 声发射事件数代码从bond断裂到AE事件的完整统计链路4.1 声发射在PFC中的物理对应平行粘结断裂与应变能释放PFC里不能直接模拟声波传播但声发射的物理本质是材料内部微破裂释放弹性波。在离散元框架里这个“微破裂”对应的就是平行粘结的断裂。每断一根bond局部储存的应变能瞬间释放宏观上对应一个AE源。所以“声发射事件数”代码的第一层就是统计bond断裂的数量。但这里有一个指标定义的分岔实验室里声发射仪一个“事件”往往对应若干微破裂的集合而不是每一根bond断裂都算一个独立事件。如果直接把bond断裂数当事件数数值会偏大而且不同加载阶段的比例关系也会失真。我在模拟里同时统计两个口径一是“累计断裂总数”二是“声发射事件数”。后者用一个时间窗口合并规则如果在一个时间窗口内比如100个计算步内连续检测到新增断裂且断裂位置之间的空间距离小于某一个阈值比如3倍平均ball半径就把这些断裂归为同一个声发射事件。这个口径和实验室AE仪的分辨率思路是一致的。4.2 事件数统计的核心代码逻辑与采样策略统计事件数的代码要解决两个问题怎么检测断裂发生以及怎么在不拖慢计算的情况下做检测。检测断裂最直观的方法是遍历所有接触检查平行粘结的断裂标志。PFC5.0里平行粘结断裂后接触属性pb_break会从0变成1。我的做法不是每步都遍历那样会占用大量额外计算时间而是每50步采样一次记录当前的断裂总数然后和后一个采样点的值做差得到这个时间窗内的新增断裂数。def ae_sampling local n_broken 0 local c contact loop while c # null if contact.prop(c, pb_break) 1 n_broken n_broken 1 endif c contact.next(c) endloop ae_new n_broken - ae_broken_prev ae_broken_prev n_broken end这段逻辑比较朴素不同版本里contact的遍历函数名可能有差异核心思路就是快照差分。采样间隔长一点效率高间隔短一点事件时间分辨率也高。我做声发射定位的时候用50步只做事件数统计用100步就够了。真正要小心的是“事件合并”的实现。我一般在FISH外部用Python处理先把每次采样的新增断裂列表连同坐标写进CSV然后按时间窗口和距离阈值聚类。这样把统计逻辑和计算过程解耦调试方便得多也方便做不同的窗口参数敏感性测试。4.3 空间定位、能量统计与震级-频度曲线如果要画声发射空间分布必须要记录每次断裂发生的位置。我建议在采样时顺便输出断裂的x、y坐标和对应计算步。PFC5.0里可以通过contact.pos(c)取接触坐标断裂之后接触还在只是粘结断开了所以坐标是拿得到的。能量统计是另一个层次。平行粘结断裂释放的能量可以粗略用断裂瞬间的接触力乘以相对位移来估算。严格一点可以读取断裂前一个计算步的接触力分量和位移增量做点积。这个做法计算量稍大但只要在采样间隔里插入一次“断裂前状态记录”大部分情况都能接受。有了事件数和能量之后可以做震级-频度统计把每个事件能量取对数作为“震级”M然后统计不同震级区间的事件频次N拟合logN与M的关系。斜率的绝对值就是b值。我在自己的模拟里发现临近峰值应力时b值明显下降也就是大事件比例增加这和实验室声发射监测到的“b值前兆”规律一致。能做到这一步模拟就从简单的强度复现进阶到破坏机理解释了。5. 调试现场十次模拟里八次会遇到的坑与排查清单5.1 计算速度与精度的平衡clump数量与时间步的取舍模拟碎石混凝土最大的痛点是计算量。clump虽然只有几十颗但每颗内部有6到8个pebble整个模型几万个颗粒接触数上万遍历一次就是一次不小的开销。如果再加上每50步一次的接触全遍历计算时间会明显上升。我的取舍原则是先跑小模型验证逻辑再放大到目标模型。小模型里clump数量减半、ball粒径放大把加载控制和声发射统计脚本全部跑通排掉逻辑错误正式模型只跑标定后的参数一遍过。这样总时间反而是最短的。时间步也要注意。PFC会自动计算临界时间步但当gjk接触或者clump形状比较尖锐时自动时间步会偏保守算得慢。不要为了提速去手动放大时间步一旦数值不稳定前一秒跑得飞快后一秒直接发散前面的墙、加载、断裂全废了。5.2 峰值强度和破裂模式不对时先按这个顺序排查现象一峰值强度偏低。不要急着加bond强度。先查是不是加载速率太快导致惯性卸载再查初始应力平衡是否充分。这两个问题不解决加多少强度都白费。确认没问题后再按“粘聚力-抗拉强度”的顺序调参。现象二破裂模式变成竖直劈裂没有剪切带。常见原因是刚度比太大或者界面过渡区强度太弱导致裂纹全部沿骨料边缘走。先降刚度比再调界面粘结强度一般能拉回剪切破坏模式。现象三峰后曲线掉得太快接近瞬断。多半是抗拉强度相对于粘聚力太低张拉裂纹过早贯通。把抗拉强度稍微抬一点通常能让峰后曲线变得平缓一些。现象四声发射事件数太少、全挤在一个时间点。查采样间隔是否太大导致细微信号没捕捉到再查加载速度是否过快把微破裂压缩到一个瞬间。我遇到过采样间隔200步时事件数看起来只有几个改成50步之后同一个模拟能数出几百个事件。5.3 一张可打印的检查清单我把平时调试用的一张检查清单放在下面每次跑正式模拟前过一遍能省掉大量的重复试算。检查项常见问题处理建议clump面积占比投放率不足或过高统计实际面积占比控制在35%到45%初始应力平衡最大不平衡力偏大平衡到平均接触力的1%以下再加载局部阻尼0.7过阻尼标定和正式模拟固定为0.1到0.2加载速率高速惯性效应做速率敏感性测试确认峰值差异小于2%伺服增益增益过大振荡用alpha 0.3到0.5的经验值试跑观察bond断裂统计采样间隔过大50步采样一次必要时同步输出坐标事件数口径未做时间窗口合并统一窗口和距离阈值记录口径说明输出文件只存应力应变缺断裂数至少包含轴向应变、应力、累计断裂数我自己最深的体会是PFC模拟最终能逼真多少不取决于软件本身而取决于建模人对材料破坏物理过程的理解程度。clump模拟碎石骨料、弱化界面过渡区、统计bond断裂并合并成声发射事件这些操作的本质都在回答同一个问题——“这个试件的破坏到底是从哪里开始、怎么扩展、什么时候失稳的”。把这个问题想清楚了参数标定和代码调试都只是时间问题。
阅读完成 · 觉得有帮助？
