PFC2D模拟层理岩体蠕变这事儿说难不难说简单也绝对不简单。早些年我做单轴蠕变模拟时最头疼的就是加载方式的选择——单级加载简单粗暴但容易“一刀切”分级加载贴近实际但控制逻辑繁琐剪切蠕变更是得单独搭建剪切盒模型。今天这篇就把这三种加载模式的实现思路和代码细节串起来讲希望能帮你少踩几个坑。1. 层理岩体模型的构建与参数选择1.1 层理结构怎么在PFC2D里“长”出来层理岩体和均质岩体最大的区别就是存在一组或者多组结构面。在PFC2D里模拟层理面常用的做法是在颗粒接触层面做文章——不是所有接触都用同一种模型而是在层理位置的两侧颗粒间赋予不同的接触本构也就是所谓的 smooth-joint 模型简称SJM。我个人的习惯是先生成完整岩石颗粒体再根据层理面的倾角和间距在指定位置“切割”接触。实际操作中这里有个细节特别容易忽略切割后两侧颗粒之间的SJM接触法向和切向刚度不能直接照搬颗粒接触的参数而是要单独标定。否则做出来的层理面要么过硬形同虚设要么过软整个试样沿着层理直接散掉。1.2 细观参数标定的三个关键点PFC2D里的细观参数和室内试验的宏观参数之间没有一一对应关系所以需要通过数值试验反复试算。做蠕变模拟前至少要标定以下几个方面完整岩石部分的接触模量、刚度比、黏结强度。层面SJM接触的法向刚度、切向刚度、内摩擦角、黏聚力。蠕变本构中涉及的黏性参数这取决于你用的是burger模型还是PSC模型。标定时的建议是先做单轴压缩再做巴西劈裂最后才是蠕变。很多人上来就标定蠕变参数结果发现弹性段都对不上那就是基础没打好。2. 蠕变模型的理论基础与代码实现路线2.1 蠕变不是“慢速加载”而是“恒定荷载下的时间效应”这是我见过最多的一个误区。很多新手误以为把加载速率放慢就是蠕变。实际上蠕变的核心特征是在应力不变的情况下变形随时间的持续增长。在PFC2D里要实现这一点需要给模型赋予时间相关的本构关系。PFC2D内置的Burger模型比较好地模拟连续介质蠕变但是对于离散元我们往往更关心微裂纹扩展、层理面的滑移错动所以实际中我更倾向于采用PSC模型并行黏结应力腐蚀模型因为它可以考虑亚临界裂纹扩展从而再现蠕变的三阶段特征。2.2 PSC模型的代码搭建思路PSC模型的工作原理说起来也不复杂当接触应力超过一定阈值后黏结直径会随时间逐渐减小这就是“应力腐蚀”。当黏结直径减小到零时微裂纹产生。核心FISH代码如下def stress_corrosion loop foreach cp contact.local ; 获取法向应力 sig_n contact.prop(cp,pn) if sig_n sig_threshold then ; 黏结半径衰减 rad_new radius * (1 - corrosion_rate * dt) contact.prop(cp,radius,rad_new) end_if end_loop end注意这里的腐蚀速率并不是一个定值它跟湿度、温度都有关系。如果只是做干燥环境下的短期蠕变腐蚀速率可以给一个很小的值比如1e-6/s量级。2.3 时间步与真实时间如何对应离散元的时间步是动态计算的通常非常小比如1e-8秒量级。但蠕变试验持续几小时甚至几天如果按真实时间算PFC2D根本跑不完。所以通常采用“时间加速法”因为PSC模型中的腐蚀速率是基于真实时间的我们需要在每次循环中把虚拟时间步放大。我的做法是在FISH中设置一个scale_factor用scale_factor乘以真实时间增量来加快腐蚀过程。需要注意的是这个值不能太大否则会出现数值失稳。3. 单级加载蠕变的实现与代码细节3.1 单级加载的基本流程所谓单级加载就是一次性把应力加载到目标值然后保持恒定观察蠕变变形。这在代码里实现起来相对简单主要分为两步首先加载阶段。我们可以通过伺服控制来给试样逐步施加荷载直到达到目标应力。其次蠕变阶段。关闭伺服机制保持应力恒定让模型在蠕变本构作用下继续变形。加载速率的控制很关键。如果加载太快会激发过大的动能力效应导致初始应变突变如果太慢又会浪费计算时间。我的经验是以伺服控制加载时将加载速度控制在每秒应变0.01%左右比较合适。3.2 单级加载的伺服控制代码伺服控制的本质是通过调整墙的移动速度使墙面的测量应力逐步逼近目标应力。PFC2D内置了伺服控制机制但我们也可以自己写FISH实现更精细的控制def servo_control w wall.find(1) st wall.stress(w) err target_stress - st if abs(err) tolerance then wall.vel(w) 0.0 status 1.0 ; 表示已达到目标应力 else g gain_coefficient * err wall.vel(w) g end_if end当status等于1.0的时候我们就把伺服关闭进入蠕变阶段。这个切换点是否准确直接决定了蠕变试验的初始条件是否准确。3.3 单级加载的蠕变阶段代码蠕变阶段不再需要控制应力我们只需要记录时间、应变等信息。为了防止模型漂移最好同时记录墙体位置和试样的平均应力确保应力没有出现明显回退。def creep_monitor ; 记录当前虚拟时间 time_now time.scale ; 计算应变量 strain_now (wall.pos_y - wall.pos_y0) / height0 history.add(strain, strain_now) history.add(time, time_now) end单级加载适合快速判断岩体在不同应力水平下的蠕变特征但缺点是无法获得完整的蠕变曲线族。这时候就需要分级加载上场了。4. 分级加载蠕变的实现与常见坑4.1 分级加载的核心路径切换分级加载的核心就是同一样品在多个应力水平下依次进行蠕变试验。每个应力水平的蠕变阶段结束后继续加载到下一个应力水平。实现上有两种常用路径逐级加载法加载、蠕变、再加载、再蠕变直接串联。阶梯加载法以恒定速率持续加载中间不加保持段通过应力历史重构蠕变。在PFC2D里我更推荐前者。逐级加载法更容易控制也方便观察每一级应力水平下的瞬时弹性应变和蠕变应变的分界。4.2 如何判断蠕变阶段是否结束分级加载最难判断的就是“上一级什么时候算完下一级什么时候开始”。标准做法是当应变速率小于某个设定阈值时认为当前蠕变阶段已经进入稳态或减速阶段可以进入下一级加载。这个阈值需要根据具体岩性调整。软岩可以放宽到1e-7/s硬岩则要更严格一些。设定的阈值过于严苛模拟时间会爆炸过于宽松各级变形数据就没有区分度。4.3 分级加载的FISH实现实现思路是通过一个状态变量state来控制当前处于加载还是蠕变阶段。当伺服控制完成且蠕变稳定后提高目标应力重复循环。def staged_creep target_stress target_stress stress_step while_stepping if state loading then servo_control else creep_monitor if strain_rate threshold then state loading end_if end_if end_loop end需要注意在切换到下一级加载时上一级蠕变存在的塑性变形不会消失。也就是说下一级加载的应力-应变曲线起点不是在原点而是在上一级终点的基础上继续。这其实是真实岩石试验中常见的“应变硬化”现象也是分级加载能获得完整蠕变曲线的优势所在。4.4 分级加载常见问题排查分级加载最容易出现的问题是蠕变阶段应力漂移尤其是加载墙上的应力无法完全恒定。这是因为PFC2D中颗粒间的接触数目有限单次蠕变过程中微裂纹扩展会导致局部应力重分布从而引起墙体应力的微小波动。解决手段有几个一是增加墙体接触颗粒的数量也就是加密模型二是在蠕变阶段不要完全关闭伺服而是开启一个弱伺服来调节三是记录蠕变阶段的平均应力而不是瞬时应力应力漂移的影响就会小很多。5. 剪切蠕变模拟的建模与实现5.1 剪切蠕变需要建什么模型剪切蠕变是岩体沿结构面发生剪切位移随时间增长的现象对层理岩体来说尤其重要。室内试验的设备是直剪仪数值模拟中也需要建立对应的双盒剪切模型。基本的建模思路是建立上下两个剪切盒中间为层理面位置。上盒施加恒定的法向应力下盒施加恒定的剪切力。监测剪切位移随时间的变化。在PFC2D中剪切盒通常用墙来模拟。上下盒的颗粒分别生成但接触属性需要由层理面SJM参数来控制。5.2 剪切蠕变中的法向与切向加载控制剪切蠕变模拟中法向应力和剪切应力都必须是恒定的。法向应力通过伺服控制上盒顶部墙体实现而剪切应力则通过给下盒施加恒定速度或恒定力实现。恒定剪切力的实现在PFC2D中比较微妙。可以直接给下盒墙体施加一个外力也可以每步检测墙体应力并修正墙体速度这个逻辑跟单轴蠕变中的伺服控制非常相似。这里给出一个恒剪力控制的FISH片段def shear_servo shear_stress wall.stress(w_shear) err shear_target - shear_stress gain 0.5 * err / sample_length wall.vel(w_shear) gain end注意剪切力的加载速率同样不能过快。过快会导致结构面发生瞬间剪切破坏那么后期蠕变阶段就没有观察价值了。5.3 剪切蠕变的典型结果是什么对于层理岩体剪切蠕变通常表现出三个阶段瞬时剪切变形阶段加载瞬间产生一定剪位移。蠕变阶段剪位移随时间增长速率逐渐衰减。加速破坏阶段如果剪切应力超过长期强度最终会进入加速蠕变直至剪断。在PFC2D中如果采用的是PSC模型可以捕捉到结构面上微裂纹的萌生—扩展—贯通全过程。这也是离散元相较于连续介质方法的最大优势——你能直接看到破坏是怎么沿着层理面演化的。5.4 剪切蠕变的后处理技巧剪切蠕变模拟的后处理重点是绘制剪切位移-时间曲线和剪切位移速率-时间曲线。速率曲线如果出现明显的弯折说明试样进入了加速蠕变阶段这个点对应的剪应力就可以推断为该应力水平下的长期强度。此外强烈建议输出结构面上颗粒的位移矢量。剪切蠕变的破坏模式往往是不均匀的——有的区域先滑移有的区域还锁住。单看位移-时间曲线很容易掩盖这种空间非均质性而位移矢量图能让机理更直观。6. 实操踩坑与经验参数速查6.1 六个最常见的模拟失败原因PFC2D蠕变模拟跑崩的情况我见过不少次了。归纳起来有六种情况最常见第一参数标定不闭环。做蠕变前没有充分标定瞬时力学行为弹性模量都偏差10%以上这时候蠕变结果没有意义。第二时间步设置过大导致接触力振荡。PFC2D稳定时间步是自动计算的但我们如果引入时间加速经常会破坏计算的稳定性。建议在加速因子超过100时逐步测试。第三模型尺寸效应。很多细观参数的标定结果依赖于模型尺寸。不同尺寸的模型蠕变结果可能差异巨大所以模拟时应统一试样尺寸与标定试验保持一致。第四墙的刚度设置不合理。墙刚度过大会产生高频振动过小则应力控制失效。第五蠕变阶段记录频率太低。如果记录步数太少蠕变曲线的早期减速阶段会看不清楚影响稳态蠕变速率拟合。第六忽略阻尼影响。PFC2D默认的局部阻尼会消耗额外的能量这个在弹性阶段问题不大但在蠕变阶段可能造成变形偏小。建议蠕变阶段把局部阻尼调到很低的水平。6.2 参数参考表下面是个人在模拟层理岩体蠕变时常用的参数范围供参考参数项取值范围备注颗粒半径0.5-1.5 mm不宜过小否则计算量爆炸颗粒刚度比1.5-3.0影响泊松比SJM法向刚度0.1-0.5倍颗粒接触刚度过大会失去层理面效应SJM切向刚度0.1-0.5倍法向刚度影响滑移响应腐蚀速率1e-6 ~ 1e-4 /s软岩取大值蠕变阶段阻尼0.05-0.15比标准值低时间加速因子10-500需要逐步校核6.3 参数校准怎么才能不返工标准做法永远是“先瞬时后蠕变”。先复现单轴压缩全过程确保破坏模式和峰值强度都大致吻合再加时间依赖性。如果连瞬时单轴都差很远那就不要指望蠕变能对上。更稳妥的做法是用多级单轴蠕变试验的应变-时间曲线一起校准PSC参数。目标不是拟合单个应力水平而是同一组参数同时复现多个应力水平下的蠕变曲线这种“跨工况校验”是判断参数可迁移性的有效手段。7. 写在最后的一点心得体会做了几年层理岩体蠕变模拟最大的感受是PFC2D这个工具看着自由度大实际上约束也很多。蠕变模拟尤其如此参数多、时间尺度跨越大、结果对细观参数异常敏感。但也正是这种敏感逼着你去真正理解岩体的变形破坏机制。如果真的打算用PFC2D在蠕变方向上做点东西建议从这些问题开始思考你的层理面是用SJM还是平行黏结替代蠕变本构是选用内置的Burger模型还是用PSC应力腐蚀模型时间加速因子是拍脑袋定的还是经过校验的这些问题的答案往往直接决定你模拟结果的可靠性。PFC2D是吃透机理的好工具但前提是我们得先用理性把它管住而不是让它用不可控的结果把我们教训一顿。