企业官网建设流程全解析

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 参数参考表

下面是个人在模拟层理岩体蠕变时常用的参数范围，供参考：

6.3 参数校准怎么才能不返工

标准做法永远是“先瞬时后蠕变”。先复现单轴压缩全过程，确保破坏模式和峰值强度都大致吻合，再加时间依赖性。如果连瞬时单轴都差很远，那就不要指望蠕变能对上。

更稳妥的做法是，用多级单轴蠕变试验的应变-时间曲线一起校准PSC参数。目标不是拟合单个应力水平，而是同一组参数同时复现多个应力水平下的蠕变曲线，这种“跨工况校验”是判断参数可迁移性的有效手段。

7. 写在最后的一点心得体会

做了几年层理岩体蠕变模拟，最大的感受是：PFC2D这个工具，看着自由度大，实际上约束也很多。蠕变模拟尤其如此，参数多、时间尺度跨越大、结果对细观参数异常敏感。但也正是这种敏感，逼着你去真正理解岩体的变形破坏机制。

如果真的打算用PFC2D在蠕变方向上做点东西，建议从这些问题开始思考：你的层理面是用SJM还是平行黏结替代？蠕变本构是选用内置的Burger模型还是用PSC应力腐蚀模型？时间加速因子是拍脑袋定的还是经过校验的？这些问题的答案，往往直接决定你模拟结果的可靠性。

PFC2D是吃透机理的好工具，但前提是我们得先用理性把它管住，而不是让它用不可控的结果把我们教训一顿。