1. 放着现成的PFC 2D不用我为什么偏要写一个先说清楚背景我去年接到一个任务分析一组颗粒材料的单轴压缩力学响应说白了就是看这种材料在竖向压力下怎么变形、怎么破坏、峰值强度大概是多少。常规做法很直接——用PFC 2D建个模型设好颗粒半径、密度、刚度、摩擦系数生成试样然后让上墙往下压跑完导出应力应变曲线。但真上手之后我发现一个问题PFC 2D帮我算得确实很稳但“稳”得有点像黑箱。我改几个参数solve一下曲线出来了可我心里对每一步计算过程并没有建立起直觉。比如接触力到底是怎么从颗粒重叠量算出来的时间步为什么必须设成那个数量级墙体伺服加载的原理是什么如果只是用软件这些问题可以永远不回答也能出结果。但你想做参数规律分析、想给组里学生讲清楚原理、或者想改一些PFC默认行为的时候就非常被动。PFC的Fish语言我也试过写多了觉得不够顺手而且很多内部数据结构想dump出来分析绕来绕去挺麻烦。于是我做了一个比较“笨”的决定用Python从零手写一个简化版PFC 2D专门用来做单轴压缩试验模拟。不是要替代商业软件而是要搭一个自己能看懂每一步、能任意改、能输出任意中间量的DEM仿真框架。目标很小圆盘颗粒、线性接触、速度和位移积分、墙体加载、应力应变输出跑通单轴压缩得到的曲线趋势和PFC结果能对上就够了。这篇文章就记录一下整个实现过程。里面包含DEM离散元法核心原理的拆解、Python代码骨架、单轴压缩试验的建模细节、数据后处理以及我在实测中踩过的几个坑。适合的人大概是这三类第一类是真的想学离散元原理但不想被商业软件绑架的第二类是已经在用PFC但希望搞明白它内部逻辑以便更合理设参数的第三类是需要在Python里做快速原型验证不想每次试验都开重型商业软件的。先给你看结论我最终用两千多个颗粒、跑了几万步得到了一条像模像样的应力应变曲线峰值强度、峰后软化、破坏模式都能复现出来。代码量不大核心部分不到三百行。整个过程跑下来我对DEM的理解确实比以前深得多。2. 手搓DEM的第一道门槛接触力模型和时间积分你要是直接打开PFC手册看离散元理论大概率会被各种接触模型吓到。但单轴压缩试验最基础的需求其实只需要最简单的线弹性接触模型就够了。搞清楚两个核心方程整个程序就立起来了。2.1 颗粒与颗粒之间怎么算力线性接触模型DEM的基本假设是颗粒被认为是刚体但允许颗粒在接触点产生微小重叠重叠量用来度量变形进而计算接触力。这句话是整个方法的灵魂。你试想一下两个圆形的颗粒碰到一起如果完全刚性那它们之间力的传递是瞬时的没法计算弹性力学又太复杂得考虑应力分布。PFC的处理方式很工程化把接触当成一个“罚弹簧”重叠越多法向接触力越大。力的公式极其简单fn kn * overlap其中fn是法向接触力kn是法向刚度overlap是两颗粒半径之和减去两颗粒中心距离也就是重叠量。这个公式看起来简单到甚至有点“不真实”但它恰恰是DEM能以显式方式计算的核心因为力的计算只依赖当前时刻的几何状态没有任何迭代收敛过程。当然只有法向力还不够。现实中颗粒接触还有切向摩擦PFC的做法是引入一个切向弹簧记录接触点处的切向相对位移增量乘上切向刚度ks得到切向力。同时切向力受Coulomb摩擦极限约束fs_max mu * fn if fs fs_max: fs fs_max这个公式意味着颗粒之间的切向力不会无限增长一旦超过摩擦极限就会发生滑动。这是模拟颗粒材料摩擦特性的关键后面单轴压缩的剪切强度很大程度就靠它贡献。2.2 牛顿第二定律怎么积分时间步与Verlet积分有了接触力下一步就是把力变成加速度、速度、位移这个循环就是DEM的时间推进。对每个颗粒受到的合力F包括所有接触力、重力、墙体反力除以质量m得到加速度a。然后需要把加速度积分成速度速度积分成位移这就是运动方程F m * a数值积分方式我用的是Velocity-Verlet格式它分两步走先是更新半步步速度再更新位置最后用新的受力更新全步步速度。在代码里看起来更直观for p in particles: p.vx (p.fx / p.mass) * dt p.vy (p.fy / p.mass) * dt p.x p.vx * dt p.y p.vy * dt每次更新完位置后需要根据新的颗粒坐标重新计算接触力再做下一轮更新。这个循环往复就是整个DEM模拟的核心。2.3 时间步长怎么选瑞利波速与稳定条件手搓DEM最容易被忽略也最容易“炸”的就是时间步长。如果dt取得太大系统会明显不稳定颗粒会越飞越快最终“爆炸”取得太小计算量成倍增加明明几百步能跑完的偏要跑几千步。时间步长的理论上限跟颗粒的“接触刚度-质量”系统有关工程上最常用的估计公式是基于瑞利波速的dt_critical 0.2 * pi * r_min * sqrt(rho / kn)其中r_min是最小颗粒半径rho是颗粒密度kn是法向刚度。公式背后的物理含义是单个颗粒表面弹性波传播穿过颗粒需要的时间是这个接触系统的最小特征时间时间步如果超过这个尺度信息传递就会有明显误差。实际使用中安全系数通常取0.1~0.2也就是dt 0.1 * np.pi * r_min * np.sqrt(rho / kn)我实际测下来这个dt数值大约在1e-6秒量级具体取决于你选的参数。一个单轴压缩试验模拟到10%应变可能需要几万步但Python跑起来也就是几十秒到几分钟的事完全可接受。3. Python代码骨架颗粒类、接触搜索和主循环说到代码实现我觉得最大的挑战不是物理公式本身而是怎么把数据结构组织得既清晰又高效。PFC内部是怎么组织这些的我们看不到但自己写的时候自由度很大我采用了“类和numpy数组混用”的方式平衡可读性和计算速度。3.1 颗粒和墙体的数据结构每个颗粒本质上就是个圆盘在2D平面内有坐标、半径、速度、受力。用类封装起来非常清晰class Particle: def __init__(self, x, y, r, mass): self.x x self.y y self.r r self.mass mass self.vx 0.0 self.vy 0.0 self.fx 0.0 self.fy 0.0墙体更简单在单轴压缩里只需要上下两个水平墙每个墙体用一个y坐标加一个当前速度表示。上墙向下运动下墙固定。墙和颗粒之间的接触力计算方法与颗粒间接触类似只是把墙看作一个半径无穷大的平面重叠量等于颗粒半径减去颗粒中心到墙面的距离# 上墙处理 if wall_top - p.y p.r: overlap p.r - (wall_top - p.y) fn kn * overlap p.fy fn这样处理完所有墙体接触后每个颗粒上的合力就齐了可以更新速度和位置。3.2 邻居检测几百个颗粒可以暴力几千个最好用网格一开始我写了个最简单的实现遍历所有颗粒对逐一判断是否接触也就是O(N²)复杂度。当颗粒数只有四五百个的时候这个方案完全没问题几千步一眨眼就跑完了。但我的试样需要2400个颗粒左右O(N²)就变成大概290万次距离判断每步虽然也没到跑不动但明显感觉到慢了。于是升级成均匀空间网格Cell List方案把模拟区域划分成边长略大于最大颗粒直径的网格每个颗粒根据坐标落入某个网格只需要检查同一网格及相邻网格中的颗粒是否接触即可把每步的颗粒对判断数量降低一个到两个数量级。网格方案的代码写起来也很直接def build_cell_index(particles, cell_size): cell_dict {} for i, p in enumerate(particles): cx int(p.x / cell_size) cy int(p.y / cell_size) cell_dict.setdefault((cx, cy), []).append(i) return cell_dict然后在每个时间步里只对相邻9个网格做颗粒对接触检测。这个优化让2400颗粒的模拟从“能跑”变成“跑得飞快”。3.3 力计算与运动更新一个时间步里发生了什么整个模拟主循环核心逻辑可以总结成四步置零受力、计算接触力、更新速度位置、更新墙体位置。这四步对应代码如下for step in range(n_steps): # 1. 力清零 for p in particles: p.fx p.fy 0.0 # 2. 颗粒间接触力 cell_dict build_cell_index(particles, cell_size) for (cx, cy), ids in cell_dict.items(): for neighbors in get_adjacent_cells(cx, cy): for i_idx in ids: for j_idx in neighbors: if j_idx i_idx: continue compute_contact(particles[i_idx], particles[j_idx]) # 3. 墙体接触力 for p in particles: if p.y p.r wall_top_y: compute_wall_contact(p, wall_top_y, is_topTrue) if p.y - p.r wall_bottom_y: compute_wall_contact(p, wall_bottom_y, is_topFalse) # 4. 更新运动 for p in particles: p.vx (p.fx / p.mass) * dt p.vy (p.fy / p.mass) * dt p.x p.vx * dt p.y p.vy * dt # 5. 更新上墙位置伺服控制 wall_top_y wall_v_top * dt这个主循环的妙处在于它每一步的物理意义非常清楚先看当前几何状态下颗粒受力如何再根据这个力更新几何状态。几何变了下一轮的受力自然跟着变如此循环就模拟出了颗粒材料的力学响应。这里有个非常关键的细节力更新完后才能更新速度位置顺序不能反。如果先更新位置再算力相当于用旧的力去推新的位置会引入一个时间步的错位虽然差值不大但对于长期累计的模拟来说会造成能量误差。3.4 完整代码骨架上面这几个片段拼起来加上参数设置就是一个能跑的DEM模拟器了。这里给一份我最终整理出来的核心结构方便你照着搭import numpy as np class Particle: def __init__(self, x, y, r, mass): self.x x; self.y y; self.r r self.mass mass self.vx 0.0; self.vy 0.0 self.fx 0.0; self.fy 0.0 def compute_contact(p1, p2, kn): dx p2.x - p1.x dy p2.y - p1.y dist np.hypot(dx, dy) overlap p1.r p2.r - dist if overlap 0: nx dx / dist ny dy / dist fn kn * overlap p1.fx fn * nx p1.fy fn * ny p2.fx - fn * nx p2.fy - fn * ny def run_simulation(particles, ...): # 主循环 for step in range(n_steps): # 受力清零、计算接触力、更新运动... pass真正跑单轴压缩之前需要把颗粒试样生成好这个过程本身也有不少讲究下一节详细说。4. 搭一个能出结果的单轴压缩试验模型有了DEM内核接下来最需要花心思的就是“试验模型”本身怎么搭。单轴压缩试验的模拟流程和物理试验其实是完全对应的制备试样、放好上下压板、以一定速度压缩、记录位移和力。4.1 颗粒生成随机撒点与初始孔隙率控制生成颗粒试样最常见的是两种思路。一种是直接按目标孔隙率随机投放让颗粒随机落到矩形区域中遇到重叠就挪走反复尝试直到填满另一种是先生成一组较小半径的颗粒然后逐步放大半径直到达到目标紧密程度。我采用的是后者因为对于单轴压缩这种需要密实试样的试验来说更好控制。半径膨胀法的思路是这样先在矩形区域中随机投放半径很小比如目标半径的0.5倍的颗粒确保没有任何重叠然后每个时间步把所有颗粒半径略微增大一点增大的过程中颗粒间开始接触、产生接触力颗粒会慢慢被推开系统重新达到平衡。重复这个过程直到半径达到目标值试样就变得非常密实了。这个方法的实践里有个小技巧膨胀速率不能太快否则颗粒来不及调整位置会形成局部应力集中的“锁死”状态。我一般把每次膨胀量控制在一个时间步内颗粒最大位移的十分之一以下让系统有充足时间松弛。4.2 上墙怎么压应变率伺服与应力伺服单轴压缩的加载方式物理试验里通常是位移控制或力控制。模拟里最常用的是应变率控制让上墙以一个固定的速度向下移动相当于恒定应变率加载。上墙速度的计算target_strain_rate 0.1 # 1/s wall_top_v target_strain_rate * sample_height注意随着试样变形当前高度不断减小如果严格按照“恒定应变率”来算墙速也应该随高度递减。不过工程实践里只要压缩幅度不大峰值强度前一般小于5%应变直接取初始高度乘目标应变率作为常数墙速误差可以忽略。比应变率更接近真实试验的是应力伺服控制设定目标轴向应力动态调整墙速使墙反力逐渐逼近目标值。这在PFC里一般用伺服机制实现代码上的核心是一个比例控制器current_stress wall_force / sample_width error target_stress - current_stress wall_v_top probe_gain * error这个probe_gain的选取有讲究太大会导致墙速振荡太小则加载太慢。我用了一个简单估算把墙视作和接触弹簧串联的系统增益设为墙刚度对应临界阻尼的0.5倍左右实测下来收敛还行。4.3 边界条件光滑墙还是摩擦墙对结果影响巨大在真实单轴压缩试验里压板和试样端面之间的摩擦是一个极其重要的因素。如果压板和端面完全光滑试样被压缩时端面可以自由横向膨胀整个试样变形比较均匀破坏模式倾向于中间鼓出不发生剪切带时如果压板端面有摩擦端面处的颗粒不能自由横向滑动试样会形成比较明显的“X型”剪切带或上下锥形破坏区。模拟里对应的设置是墙上有没有切向摩擦。我的代码里默认墙体接触只计算法向力相当于理想光滑墙。要做摩擦墙就在墙体接触里同样加上切向弹簧和Coulomb摩擦极限。不要小看这个细节同等参数下摩擦墙得到的峰值强度可能比光滑墙高20%甚至更多。这也是很多新手用PFC做单轴压缩时曲线总对不上的重要原因——你的物理试验用的压板到底是什么摩擦情况模拟里必须匹配。在颗粒生成时还有一个细节试样在矩形区域内生成后两侧是自由边没有侧向约束所以试样会凭借颗粒间的接触力和摩擦力维持自身形状。这意味着初始阶段试样内部已经存在一定的初始接触应力场如果初始应力太高会直接影响后续加载曲线。所以颗粒生成完成之后一定要让系统先松弛几百步再开始正式加载否则曲线开头会有明显的不稳定跳变。5. 数据后处理从墙反力到应力应变曲线模拟跑完原始输出只是一堆颗粒坐标、速度、接触力数据离“应力应变曲线”还有一步转化。这一步看似简单其实也有不少讲究处理不好曲线数据就是一团噪点。5.1 应力怎么算墙反力除以截面积单轴压缩试验中的轴向应力标准定义是轴向力除以试样截面积。在2D模拟里“截面积”变成了“宽度乘以单位厚度”所以我用的是axial_stress wall_force_total / (sample_width * 1.0)这里sample_width是试样宽度单位厚度取1.0米2D模型默认厚度单位。wall_force_total是上墙受到的所有颗粒接触力的合力也就是所有颗粒施加到上墙的垂直力之和。每次时间步都记录这个值就得到了一条应力随时间变化的曲线。注意一个容易出错的地方墙反力需要取绝对值而且只用Y方向的合力。有些颗粒可能由于侧向挤出对墙体产生水平力但单轴压缩中我们只关心轴向竖向的分量水平分量应该忽略或仅用于分析。为了平滑曲线我会对原始应力信号做一个滑动平均窗口大概取每1%轴向应变对应的步数。不然由于颗粒接触的离散性曲线会像锯齿一样高频振荡很难读出准确的峰值强度。5.2 轴向应变怎么定义与输出轴向应变的计算相对直接axial_strain (initial_height - current_height) / initial_height其中initial_height是试样初始高度current_height是当前上墙高度减去下墙高度。注意这里用的是“工程应变”定义即高度变化量除以初始高度而不是对数应变。在应变小于10%的范围内两者差异不大但工程应变是颗粒材料试验的标准表达方式。输出的时间序列最好每隔几步存一次不需要每步都记否则数据量会比较大。我一般是每50步记录一次2400颗粒跑几万步后导出的应力应变曲线已经足够平滑。5.3 破坏模式怎么看位移矢量图和接触力网络应力应变曲线只是最终结果之一DEM模拟最大的优势是能直接看到颗粒层面的破坏过程。单轴压缩试验的破坏模式通常可以结合两个可视化来判断颗粒位移矢量场和接触力网络图。位移矢量场的画法是在最终状态图上以每个颗粒的初始位置为起点画一条指向当前位置的箭头。颗粒材料受压后会向两侧横向膨胀如果形成剪切破坏带就能看到某一斜向条带内的颗粒位移和周围明显不同。接触力网络的画法更简单在两个接触颗粒之间画一条线段线宽正比于法向接触力大小。峰值强度时接触力网络中会形成明显的“力链”力链贯通的方向和加载方向一致破坏后力链逐渐断裂、重新定向形成剪切局域化。这两张图用matplotlib都能画虽然3D效果比不了PFC自带的可视化但2D下已经足够说明问题了。我后面的参数敏感性分析很大程度上就是靠这两类图来判断破坏模式有没有变化。6. 我实测中踩过的坑与参数标定的经验这一节专门说坑。这些坑几乎都是我实际写代码、跑数据时踩过的很多参数试错花了好几天时间拿出来分享希望你能少走弯路。6.1 颗粒数太少波形锯齿严重太多跑不动颗粒数量这个参数直接决定了计算量和曲线的平滑程度。我一开始图快只生成了大概300个颗粒跑完曲线一看峰值附近的锯齿非常厉害几乎没有规律。后来逐渐增加颗粒数量到1000个左右时曲线趋势才清晰起来。我的建议是试样宽度方向至少保证有30个以上颗粒这样力链的统计行为才比较稳定。具体到我的试样宽50mm颗粒直径1~2mm大约1500~2500个颗粒比较合适。颗粒数再多的话Python端的速度就开始吃紧了。我曾经试过半径减半、颗粒数翻4倍到接近一万个单条曲线就要跑一晚上这对于参数扫描几乎是不可接受的。这个规模说明手搓方案是有边界的后面再细说。6.2 初始填充的接触振荡Damping的引入颗粒生成阶段最容易出现的现象是试样刚生成完毕各颗粒之间存在大量初始接触力系统像“弹簧阵”一样不停振荡如果能量不耗散可能跑几千步都稳定不下来。PFC里处理这个问题靠的是阻尼机制一般有局部阻尼和粘性阻尼两种。局部阻尼的实现方式和速度方向有关作用力方向总是与速度方向相反等效于持续耗散系统动能。代码实现for p in particles: damping_force -local_damping * abs(p.fx) * np.sign(p.vx) p.fx damping_force damping_force_y -local_damping * abs(p.fy) * np.sign(p.vy) p.fy damping_force_y局部阻尼的典型值在0.2~0.7之间。注意这里有一个重要问题阻尼会额外耗散能量导致宏观峰值强度偏低所以正式加载阶段我一般会把阻尼系数调得比试件生成阶段小一些或者保持一个较小的值比如0.1~0.2。到底取多少需要拿已知材料参数去标定而不是盲目照搬。6.3 参数标定宏观弹性模量不是颗粒刚度的简单缩放接触参数中颗粒法向刚度kn直接决定了宏观弹性模量E但两者并不是简单的线性关系。颗粒材料的宏观刚度不仅取决于单颗粒的接触刚度还取决于颗粒排列、配位数、颗粒尺寸分布等。实际标定的做法是先取一组kn值做单轴压缩得到宏观应力应变曲线的线性段斜率E然后调整kn直到模拟的E和目标材料吻合。我的经验里一个粗略的近似关系是E ≈ kn / (2~4) 取决于配位数单位保持一致但这绝对不是普适公式只是初值参考。真正要让曲线吻合还是得跑几次参数标定。别嫌麻烦这一步绕不过去颗粒材料的宏观响应是非常典型的“涌现”行为单颗粒参数和宏观参数之间没有简单映射。6.4 一个容易忽略的大坑摩擦系数对峰值强度影响巨大接触摩擦系数μ是另一个极其敏感的参数。在单轴压缩中摩擦系数从0.2增加到0.8峰值强度可能翻倍都不止。原因在于摩擦系数越大颗粒之间的剪切阻力越大试样抵抗侧向变形的能力越强整体表现出的强度自然更高。这引出一个必须重视的问题如果模拟结果和物理试验对不上不要急着调kn先看看μ是否合理。颗粒材料的摩擦系数一般通过直剪或三轴试验标定而不是随意选一个。我在最初做标定时因为偷懒用了手册里的默认值0.5和实际材料大约0.8的摩擦系数差得远导致模拟强度一直偏低花了整整一天排查才发现根源在μ上。6.5 时间步相关的稳定性问题最后再加一个跟数值稳定性有关的坑就是时间步和刚度不匹配引起的“颗粒爆炸”。如果kn取到1e7 N/m以上颗粒又比较小、质量比较低而dt还是按之前的安全系数来取就很容易出现个别颗粒速度突变、飞出试样的现象。解决办法是每个时间步开始之前重新计算当前的临界时间步动态取dt或者固定取一个足够保守的dt值。动态时间步虽然更精确但会引入每步的额外计算开销而且dt变化会让“每50步存储一次”的应变对应关系变得不统一。我最后选择了保守的固定dt把kn和粒径范围控制在合理区间内这样代码简单、数据输出也好处理。7. 手搓PFC和商业PFC 2D到底怎么选看到这里你大概会有一个很自然的疑问既然自己能写是不是就不需要商业PFC了我的看法是完全不是两者解决的问题不同。手搓Python DEM的最大优势是透明度和可定制性。你能看到每一步的力是哪里来的、为什么这个颗粒会往那个方向飞、阻尼是怎么影响稳定时间的你可以任意dump中间数据不需要在PFC里折腾Fish语言输出。对于教学、原理研究、参数规律探索、快速原型验证这个方案足够好用而且代码完全属于你自己想怎么改都行。但手搓方案的边界也很明显。我的模拟只用到了圆盘颗粒和线性接触模型而真实岩土工程中颗粒形状是不规则的PFC支持通过clump生成不规则形状颗粒岩石材料需要平行粘结模型来模拟胶结力因为颗粒之间不只是摩擦还有连接纤维加筋材料需要接触模型包含拉伸和弯曲行为。这些在PFC里都已经有成熟实现和验证自己从零手搓的话工作量是指数级上升的。另外PFC在大规模三维模拟、GPU加速、流体耦合、热力耦合等方面已经形成了完整的解决方案真要算实际工程问题直接用商业软件是更负责任也更现实的选择。我现在的实操路线是两者结合用Python写参数生成脚本和批量分析工具把PFC的输入文件、参数扫描、结果后处理串成一个自动化流程遇到需要深入理解接触机理、或者某个现象解释不了的时候用自己写的简化版DEM做针对性试验。这样既能借助PFC的成熟能力做工程计算又能保留一个“看得见内部”的学习工具。说实话自己手搓完这套代码之后再回去看PFC的manual很多以前只是“知道”的参数突然变得“懂”了。比如为什么PFC默认的局部阻尼是0.7、为什么建议最小颗粒数与试样尺寸比例、为什么时间步自动计算时要取安全系数——这些在设计手册里只是一句话但当你亲手炸过一次模型、看着颗粒飞满屏幕之后就永远不会忘了。如果你也想真正理解颗粒离散元试试花一个周末手写一个最简版本我保证收益比看任何教程都大。