基于PFC的土石混合体三维冲击碾压模拟，土采用ball，石采用rblock。 能观察到每次碾压...
基于PFC的土石混合体三维冲击碾压模拟，土采用ball，石采用rblock。 能观察到每次碾压后的颗粒运动。 可更改冲击碾压速度、质量、冲击轮形态，可监测应力、位移、孔隙率变化等。
最近在折腾PFC搞三维冲击碾压模拟的时候发现个有意思的活——用ball颗粒模拟土体，rblock模拟碎石块组成的混合体碾压。这玩意儿不仅能实时看到冲击轮碾过时颗粒乱窜的场面，还能改着花样调参数看效果，今儿就跟大伙唠唠怎么玩转这个。
先上段创建材料的核心代码：
rock_lib = rblock.library.create(num=20) for i in range(20): rblock.geometry.random(rock_lib[i], vertex=8, size_range=(0.1,0.3)) # 生成2000个ball土颗粒 ball_zone = zone.create_box((-5,-5,0),(5,5,10)) ballmat = mat.create(density=2650, kn=1e8) generator = ball.generate(zone=ball_zone, radius=(0.08,0.12), material=ballmat, count=2000)这里用rblock.geometry.random搞出20种不规则石块造型，vertex=8控制每个石块有8个顶点，size_range定尺寸范围。生成土颗粒时特意把接触刚度kn设到1e8，防止计算时颗粒过度穿透，这个值要根据实际材料刚度微调。
冲击轮参数化设计才是灵魂，咱们可以这么玩：
# 冲击轮参数组 roller_params = { 'speed': 2.0, # m/s 'mass': 1500, # kg 'diameter': 1.2, 'surface_ribs': 5 # 凸肋数量 } def create_roller(pos): roller = rblock.create(position=pos, material=mat.create(density=7850)) rblock.geometry.cylinder(roller, height=0.8, radius=roller_params['diameter']/2) # 添加表面凸肋 for i in range(roller_params['surface_ribs']): angle = 360*i/roller_params['surface_ribs'] rib = rblock.geometry.create_rect(length=0.1, height=0.15) rblock.attach(roller, rib, position=(0,0,0.4), rotation=(angle,0,0)) return roller这个createroller函数整了个带凸肋的圆柱形冲击轮，surfaceribs参数控制凸肋数量。实际模拟时发现凸肋数量超过7个后颗粒飞溅会明显加剧，这说明表面纹理对碾压效果影响比想象中还大。
基于PFC的土石混合体三维冲击碾压模拟，土采用ball，石采用rblock。 能观察到每次碾压后的颗粒运动。 可更改冲击碾压速度、质量、冲击轮形态，可监测应力、位移、孔隙率变化等。
碾压过程用循环控制挺方便：
for step in range(100): roller.set_velocity( (roller_params['speed'], 0, 0) ) model.cycle(500) # 每个碾压步计算500时步 record_displacement() # 记录颗粒位移 calculate_porosity() # 实时孔隙率计算 if step % 10 == 0: export_vtk(f'output/step_{step}.vtk') # 输出可视化文件每个碾压步后记录位移场和孔隙率变化，存成VTK文件方便用Paraview做三维动态展示。这里有个坑——cycle次数不能太少，否则冲击轮还没完全压实就进入下一步了，建议先做个小规模试算确定合理时步数。
说到孔隙率计算，自己写了个实用函数：
def calculate_porosity(): total_vol = zone.volume(calc_zone) # 计算区域总体积 ball_vol = sum(4/3*math.pi*r**3 for r in ball.radius()) rblock_vol = sum(rblock.volume(rb) for rb in rblock.list()) return (total_vol - ball_vol - rblock_vol) / total_vol这里要注意计算区域要避开边界影响区，通常取中间80%的区域。实测发现当冲击轮质量从1吨增加到2吨时，孔隙率降幅会从12%提升到18%左右，但继续增重效果就不明显了，存在明显阈值效应。
玩这个模型最大的乐趣是能看到各种"翻车现场"——比如当冲击速度超过3m/s时，边缘石块会被撞飞形成溅射现象；把冲击轮改成六边形会发现应力集中区域明显改变。最近在尝试用Python写个自动参数扫描脚本，批量跑不同参数组合然后提取应力云图，有搞类似项目的朋友可以交流下怎么优化计算效率。