干过山区高填方、隧道弃渣场处理的人应该都有同感：土石混合体地基是块难啃的骨头。块石和土混在一起，级配差、强度不均匀，普通碾压设备根本压不密实，现场还容易出现测点合格、过段时间又回弹变形的怪事。冲击碾压这几年被大量用在这种场合，效果确实有，但它到底怎么在松散土石混合体里起作用——冲击波怎么传播、块石怎么被击碎、颗粒怎么重新排列——在现场靠挖探坑和弯沉测试根本看不透。
所以我把这套问题搬进了数值模型里。这篇文章想聊的，就是我怎么在PFC2D里建立一套“考虑颗粒破碎的cluster松散土石混合体地基冲击碾压二维模型”，以及在这个建模过程中踩过哪些坑、积累了什么经验。文章主要面向正在做离散元模拟的研究生、做地基处理设计的工程师，以及想在论文或项目里复现这类模型的同行。你不需要有太深的离散元基础，但至少知道PFC里颗粒、墙、接触这几个基本概念，读起来会顺畅很多。
1. 项目概述：为什么做这个模型
1.1 土石混合体地基处理的工程痛点
土石混合体广泛存在于山区机场高填方、山区高速公路路基、隧道弃渣场再利用这类场景里。它的最大特点是“同一面地基上，材料性质天差地别”：大块石可能几十厘米长，周边的土可能是粉质黏土或砂土，两者刚度相差几个数量级。这样的材料做地基，最怕的是不均匀沉降和局部剪切破坏，表面看着压平了，底下可能还有大孔隙没被挤密。
冲击碾压之所以对这种地基有效，是因为它不像普通振动压路机那样只在浅层施加循环荷载，而是通过非圆形的冲击轮在滚动中反复对地面产生低频大振幅冲击，单次冲击动能很大，能把能量传递到较深的土层，同时利用冲击产生的剪切波和挤压作用，把大块石“啃”碎并挤入土体孔隙中。我在建模前跟做现场的同事聊过，他们最关心的两个问题一个是压实后沉降量能不能达标，另一个是碎石会不会被过度破碎导致级配反而变差。这两个问题都落在颗粒尺度上，传统连续介质模型很难回答，必须上离散元。
1.2 颗粒破碎与冲击碾压的耦合关系
这个模型的第一个关键词其实不是cluster，而是“颗粒破碎”。冲击碾压过程中，松散土石混合体里的大块石承担了大量接触应力，在冲击荷载下可能发生两类破坏：一类是颗粒整体劈裂，比如块石受到两个方向挤压产生拉裂；另一类是颗粒棱角的磨损和剥落，逐渐变成小块。破碎一旦发生，颗粒的尺寸和形状改变，级配也随之改变，进而影响整个地基的压实特性和强度。
问题在于，常规的离散元模型如果直接把颗粒设为刚性的，就没法反映“块石破碎—级配改变—力学性能变化”这条因果链。模拟结果往往偏硬、偏保守，压实沉降量算出来比现场小不少，破碎耗能和颗粒重排效应也完全丢失。放到工程上，这会导致设计和施工参数选择失误，比如碾压遍数不够、冲击能选择偏小。所以模型必须把破碎这个机制显式地建进去，而且要在计算效率和破碎形似度之间找到平衡。
1.3 cluster方案为什么优于其余两类
目前离散元里模拟颗粒破碎的主流方案大概有三类：粘结颗粒模型BPM、颗粒替换法（Particle Replacement）、还有cluster（颗粒簇）。BPM是把每个大颗粒内部人为分割成若干子颗粒，子颗粒之间用平行粘结连接，受力超过粘结强度就断开，本质上是“连体婴儿”式的可破裂颗粒。颗粒替换法更粗暴：当颗粒受力超过阈值时，直接把它替换成一组预先设定好的小颗粒，破碎瞬间完成。
我最终选择cluster，原因是BPM对子颗粒排列的依赖性很强，容易产生不真实的碎片形状，而且子颗粒数量多，计算量上升明显；颗粒替换法虽然简单，但破碎前后质量不守恒，相邻颗粒的接触信息也会重置，冲击荷载作用下容易出现能量突变和应力波失真。cluster本质上是用一小簇相互粘结的子颗粒构成一个“超级颗粒”，破裂发生在粘结断裂的瞬间，子颗粒仍然参与后续接触计算，质量和动量都守恒，碎片形状也相对自然。它相当于BPM和颗粒替换法的折中方案，特别适合冲击碾压这种强动态、高频率荷载工况。
二维模型在这个问题上也很有优势：它能快速跑参数敏感性分析、大量试算不同碾压参数，计算量比三维少两个数量级。对工程预判和机理研究来说，二维模型把复杂问题压缩到了“一个剖面”，足够说明趋势和规律。
2. 模型构建关键：cluster颗粒体系怎么做
2.1 土石混合体细观结构与级配复现
在PFC2D里建土石混合体，第一步是确定“土”和“石”的比例。现场资料通常给的是质量含石率，比如40%、50%、70%几档，建模时要换算成面积含石率。换算需要知道土和石的密度，我的习惯是把二者密度都输入实测值，然后用一个简单脚本按面积比例随机投放。粒径分布则根据现场筛分或设计级配曲线来定，我这边常取块石粒径范围20～200mm，颗粒级配用Weibull分布或直接按紧凑级配生成，都能接受。
真正容易忽略的是“松散”两个字。数值模型默认会生成一个稳定平衡系统，颗粒之间接触紧密，这和真实松散地基的状态差得很远。我的做法是分两步：先生成一个目标孔隙率偏大的颗粒集合，让颗粒在自重下只形成少量接触，再在这个基础上降低颗粒间摩擦系数到0.4左右，人为模拟一种尚未充分咬合的松散结构。等冲击轮荷载上来，颗粒才开始互相接触、重新排列、逐渐密实，这样才符合“松散土石混合体被冲击碾压”的物理过程。
块石用cluster模拟，土颗粒怎么处理也需要想清楚。如果不关心土骨架内部的细观破坏，土体可以直接用小粒径圆盘颗粒表示，这部分颗粒不启用破碎，只贡献摩擦、弹性和阻尼。这样做的好处是大幅减少破碎判定的计算量，坏处是土的塑性压实行为不够细腻。在我的模型里，土颗粒粒径设为5～20mm，块石用cluster模拟，两种颗粒用不同的接触参数和组名区分，后续统计破碎率时只要按组提取就行。
2.2 cluster生成流程与粘结参数控制
cluster的生成逻辑并不复杂，核心是“先做模板，再替换”。我先按级配生成一批圆形块石，然后针对每个要模拟破碎的块石，在其内部填充若干子颗粒，子颗粒之间用平行粘结连接，最终形成一个整体。这个过程可以用PFC的clump template功能，也可以手动写循环脚本。子颗粒的直径取块石直径的1/4～1/3，数量不要太多，我试下来单个块石配6～10个子颗粒就足够，既能体现劈裂，也不至于让计算量爆炸。
子颗粒的排列方式直接影响破碎模式。我做过对比实验：同样是受力劈裂，六边形密排的子颗粒偏脆，容易碎成大片；随机排列偏韧，容易从薄弱处逐渐剥落。实际模型里我选择在块石内部采用“随机心+径向填充”的方式，兼顾劈裂和棱角磨损两种破碎形态，和室内单颗粒压碎试验的结果更接近，不至于一压就碎成渣。
cluster内部的粘结参数是整个模型的灵魂，差不多可以类比成“胶的强度”。我用了PFC里的平行粘结模型，因为平行粘结能传递力和弯矩，适合模拟岩石类材料的破裂。粘结刚度和颗粒刚度保持一致，粘结强度则和块石的抗压强度挂钩，通常在标定阶段反复调。提供一个参考值范围：法向粘结强度5～15MPa，切向粘结强度2～8MPa，粘结刚度比设为1～2，摩擦系数0.5～0.8。不过不同岩石差异很大，不要直接套，下面一节详细说标定流程。
2.3 颗粒破碎判据和参数标定
“颗粒什么时候算破碎”这个问题很容易被新手忽略。实际模型的输出是一堆粘结键断裂事件，不能把每个键断裂都算成颗粒破碎，否则破碎率会虚高。我的做法是统计每个cluster内部的粘结断裂数量，当一个cluster内断键比例超过50%（相当于主裂纹已经贯通）时，把这个cluster计入“已破碎颗粒”。这个阈值需要根据目标岩石的破裂形态调，拉裂式的取40%～50%，沿晶断裂明显的取60%以上。
参数标定的核心对照实验是单颗粒压缩。我从模型里随机取出不同粒径的单个cluster，放在两块平板之间加载，记录法向力和位移曲线，再跟室内岩石单颗粒压碎试验或者文献里的Weibull强度分布数据对比。如果模拟的破坏荷载偏低，调大粘结强度；如果曲线偏脆、峰值后瞬间失稳，调大粘结刚度比；如果碎片的碎片很多但峰值不高，说明子颗粒粒径太小或者排列太均匀，需要微调子颗粒尺寸分布。
我个人的经验是，标定阶段宁可多花两周，也不要直接跳到碾压工况。标定不好，后面所有结论都建立在一个错误的地基上，反复试算浪费的时间更多。而且标定数据最好覆盖至少三档粒径，这样才能把尺寸效应考虑进去，因为大块石内部缺陷多、表观强度低，这在数值模型里要靠不同粒径cluster的粘结强度差异来体现。
3. 冲击碾压工况设置与碾压参数选定
3.1 模型边界、阻尼和初始状态处理
二维模型的几何尺寸不能拍脑袋定。冲击轮作用范围是有限的，如果模型太窄，边界反射波会叠加到冲击区域，压实效果虚高；如果太宽，计算耗时翻倍。参考冲击碾压的影响深度一般在1～3m，我设计的地基尺寸是宽8m、高3m，冲击轮作用在模型顶部中央，两侧各留出至少3倍冲击宽度，这样能明显看到边界反射的影响逐渐衰减，主要统计区域则取中间2m范围，避开边界。
底部和两侧的墙体约束也要分清。底部墙全固定，模拟下卧硬层；两侧墙我用的是竖向自由、水平固定的约束方式，模拟半无限体，尽量减小侧向反射。阻尼参数直接影响冲击波的能量衰减速度，这是冲击碾压模拟里最容易被忽视的环节。我采用局部阻尼加粘滞阻尼的混合方案，局部阻尼系数在冲击加载时取0.1，低于常规静态计算常用的0.7，否则冲击波能量会被人工阻尼吃掉一大半，导致颗粒破碎不足。
初始状态处理上，地基生成后先运行重力平衡，等最大不平衡力比降到1e-5以下再开始碾压。这一步要特别注意颗粒间的初始重迭，重迭过大会产生虚假应力，一加载就爆掉。我通常先生成一个“松弛”的颗粒集合，再逐步缩放颗粒半径以达到目标孔隙率，整个过程用位移增量慢慢逼近，避免瞬间释放的弹性应变能。
3.2 冲击荷载怎么施加才贴近实际
冲击碾压机实际是绕轴旋转的非圆轮，每一转都会有一个冲击过程，轮子始终与地面接触，冲击是周期性的。在二维模型里，我把它简化为一个刚性矩形墙体（模拟冲击轮底部），这个墙体按照真实冲击轮的短半径做周期性运动，向下冲击后轻微反弹，再接着冲击。更粗暴的简化是一个给定速度向下运动到目标深度后反弹，但这会丢失冲击过程的波形形状，影响破碎分布。
我这里给出一个简化版的PFC2D加载控制逻辑供参考，实际写的时候要根据所用软件版本调整：
; 冲击轮加载伪代码示例，控制冲击墙的速度和位置 def impact_cycle if wall.pos.y(wall_imp) > target_y wall.vel.y(wall_imp) = -v_impact else wall.vel.y(wall_imp) = 0.0 endif end冲击速度和冲击能不是随便定的。以国内常用的三边形冲击碾压机为例，冲击轮质量约16t，设计冲击能在25～30kJ，冲击频率1.5～2Hz，碾压速度10～12km/h。模型里二维取单位厚度，所以荷载的换算要按单位宽度去做：比如模型厚度方向取1m，冲击能除以宽度后作为等效能量输入，再反推冲击轮墙体的等效质量或速度。我自己在调试中常用“三种冲击能工况”做对比，比如15kJ、25kJ、35kJ，观察破碎率和沉降量的差异，这样比较容易看出冲击能选型的趋势。
冲击次数方面，现场一般按5、10、15、20遍来控制碾压效果，模型里也要对应跑多次冲击。不同遍数之间不能简单重复同一个速度曲线，因为地基在逐遍压实，颗粒之间的接触状态变了，冲击响应必然不同，所以每遍冲击前要重新检查一次接触状态。这也是模型计算量大、需要耐心跑的原因所在。
3.3 破碎率和压实效果的统计评价
模型跑完，不能只搬一张“颗粒变多了”的图出来交差，得有量化指标。我习惯同时输出三套指标：第一套是压实指标，包括冲击区域顶面的沉降量、模型整体孔隙率变化、颗粒接触数的增加；第二套是破碎指标，包括按粒径分组的已破碎cluster数量比例、相对破碎率Br（用Hardin的相对破碎率概念，基于破碎前后粒级面积分布计算）、以及主裂纹的方位统计；第三套是能量指标，统计冲击输入能量里有多少被摩擦耗散、多少被粘结断裂耗散、多少转化为动能和弹性应变能。
我特别想提醒的是统计窗口的问题。很多新手直接把整个模型范围做统计，结果边界颗粒一动没动，把破碎率和压实效果都稀释了。正确做法是在冲击轮正下方划定一个统计窗口，比如宽2m、深1.5m的矩形区域，分别统计窗口内和窗口外的指标，再做一个沿深度方向的梯度分布。这样才能看出冲击碾压影响深度是不是跟现场经验对得上。
4. 常见问题排查：从跑不动到结果失真
4.1 计算效率低下的优化思路
第一个会碰到的问题是模型跑得太慢。松散土石混合体摩擦接触多、破碎事件多，一个2万个颗粒的二维模型运行20遍冲击，在普通工作站上可能要跑好几天。优化手段优先级是：先减少子颗粒数量，再增大接触检测步距，最后再考虑并行计算。子颗粒数量从10个降到6个，计算量能降低30%以上，而破碎模式不会发生本质改变。
另一个有效的优化是分阶段调整时间步。冲击加载阶段时间步要小，一般取自然时间步的0.5倍；冲击间隙的静置阶段可以用较大的时间步快速迭代，让系统重新平衡。这个“变速跑”的做法我在PFC里通过脚本控制，整体计算时间几乎能缩短一半。千万不要在冲击加载阶段为了赶速度强行放大时间步，冲击波传播需要捕捉，放大时间步会导致接触力失真，颗粒破碎模式完全乱掉。
4.2 破碎率异常与参数不收敛
第二类问题是破碎率要么太高要么太低。破碎率异常偏低，通常是粘结强度设得过大，或者阻尼太大吸收了冲击能量，可以按10%的幅度逐步下调粘结强度。破碎率异常偏高，尤其是大块石全部碎成粉末，多半是子颗粒粒径太小、粘结刚度过低导致内部应力集中，或者子颗粒排列太均匀形成规则薄弱面。这里我要特别强调，cluster内部的子颗粒半径一定要有随机分布，哪怕只是±10%的扰动，也能有效避免规则排列带来的“假破碎带”。
参数不收敛的情况也经常出现，我一般用“单参数变量法”排查：把颗粒摩擦系数、粘结强度、阻尼系数分别设为偏高和偏低两个值，跑一遍短工况，看哪个参数对结果影响最大。影响最大的参数优先精调，其他参数粗调即可。这个办法比无脑调参效率高太多，也能帮助理解模型的物理行为。
4.3 颗粒飞溅与能量不守恒的排查思路
冲击碾压模拟中最吓人的现象是颗粒“炸开”，模型算着算着就爆了。出现这种情况，十有八九是初始颗粒重叠太紧或者接触刚度突变导致局部应力集中，在冲击瞬间释放出来。我处理的方法是检查最大不平衡力云图，找出初始应力异常的颗粒群，重新生成局部区域，而不是整个模型重跑一遍。
能量不守恒是另一个隐蔽问题。每次冲击结束后，系统的总能量（动能+应变能+摩擦耗散+粘结断裂耗散）如果明显增加或减少，说明接触参数或阻尼设置有误。我每10步就输出一次能量和动量，画成曲线放在后台看趋势。如果冲击间隙阶段动能回不到接近零的状态，说明粘滞阻尼太小，系统还在残余振荡，这时候要让时间步继续跑下去，等稳定后再开始下一遍冲击。这个细节关系到每遍冲击初始状态的一致性，直接影响最终结果的可靠性。
5. 复盘与几个可以继续做的方向
模型基本跑通之后，我回过一次头来审视初始参数假设，发现最费时间的其实不是冲击工况本身，而是颗粒破碎参数的标定。cluster子颗粒数量、粘结强度、断裂比例阈值的组合千变万化，我最终稳定在“单块石6～8个子颗粒、平行粘结、断裂阈值50%”这一组配置，既保证破碎模块形似，又把计算量控制在实际可接受范围内。如果你也要做类似模型，我建议不要一上来就追求极致的颗粒细观形态，先把一个稳定的流程打通，再逐步细化。
这套二维模型的后续扩展方向很多。我目前在往两个方向继续做：一是把二维结果和现场压实度检测数据做对比，反过来修正接触参数，让模型具备一定的工程预测能力；二是考虑在块石cluster外围添加一层弱接触来模拟风化软壳层，这样颗粒破碎会先发生在软弱壳层而不是块石内部，可能更接近含风化碎石的实际情况。水分的影响目前还没办法用纯离散元直接考虑，如果要做得更深入，需要把离散元和孔隙水压力计算耦合起来，那也是后话了。至少从目前的应用效果看，这个模型已经能比较清楚地解释冲击碾压遍数、冲击能与块石破碎程度之间的对应关系，对设计方案的前期比选有实际参考价值。
