PFC 2D双轴压缩中为何出现颗粒重叠过大？

Question

在PFC2D双轴压缩模拟中，颗粒重叠过大常导致数值失稳或不收敛。该问题多源于接触刚度参数设置不当，如法向刚度（kn）与切向刚度（ks）过高，使颗粒间相互作用力急剧增大，引发过度重叠；或因时间步长未满足显式积分稳定性条件，造成动力不平衡累积。此外，初始孔隙率过低或边界加载速率过快，也会加剧颗粒重排不充分，诱发局部密集重叠。如何合理选取刚度比（kn/ks）、控制应变率并优化初始生成密实度，成为抑制重叠、确保模拟可靠的关键技术难点。

Answer 1

1. 问题背景与现象描述

在PFC2D双轴压缩模拟中，颗粒重叠过大是导致数值失稳或不收敛的常见问题。当颗粒间发生显著重叠时，接触力迅速增大，系统动能剧烈波动，最终可能导致模拟中断或结果失真。该现象通常表现为：

• 模拟过程中出现“NaN”或“Inf”数值；
• 应力-应变曲线剧烈震荡；
• 颗粒速度异常升高；
• 模型整体能量不平衡率持续上升。

这些症状背后的核心因素包括接触刚度参数设置不合理、时间步长选择不当、初始孔隙率过低以及边界加载速率过快等。

2. 根本原因分析

从力学与数值计算角度出发，可将重叠问题归结为以下四类机制：

1. 接触刚度过高：法向刚度（kn）和切向刚度（ks）若取值过大，颗粒间微小位移即可产生巨大接触力，导致加速度剧增，引发非物理性穿透或反弹；
2. 刚度比（kn/ks）失衡：过高或过低的刚度比会影响剪切行为与体积变形的耦合关系，影响颗粒重排路径；
3. 显式积分稳定性不足：PFC采用显式中心差分法求解运动方程，其稳定时间步长需满足：
   Δt ≤ 0.1 × √(m/k_max)
   其中 m 为颗粒质量，k_max 为最大接触刚度。若 Δt 过大，则动力不平衡累积，诱发数值爆炸；
4. 初始密实度过高或加载速率过快：低孔隙率状态下的颗粒已处于紧密接触，缺乏足够空间进行协调重排；而高速加载使系统无法动态平衡，加剧局部应力集中。

3. 关键参数调控策略

参数	推荐范围	影响机制	调整建议
kn (法向刚度)	1e7 ~ 5e8 Pa	控制压缩响应	根据材料杨氏模量等效标定
ks (切向刚度)	0.1~0.7×kn	影响剪切强度	通过刚度比调节摩擦响应
kn/ks 比值	2~5	决定泊松比与剪胀性	避免小于1或大于10
时间步长 Δt	< 0.1×√(m/kn)	保证数值稳定	自动或手动设定安全系数
初始孔隙率	0.35~0.42	影响颗粒重排能力	使用雨落法生成均匀结构
应变率	1e-4 ~ 1e-6 /step	控制加载速率	低速加载以逼近准静态
阻尼系数	0.7~0.9	耗散动能	采用局部阻尼抑制振荡
颗粒密度	2500~2700 kg/m³	影响惯性效应	与真实材料匹配
接触模型	linear/coulomb	决定力传递方式	优先选用线性接触
边界速度	1e-6 ~ 1e-8 m/step	防止冲击加载	逐步增加速率测试敏感性

4. 数值优化流程图

```mermaid
graph TD
    A[开始模拟配置] --> B{检查初始孔隙率}
    B -- 孔隙率<0.35 --> C[重新生成样本
采用雨落法或膨胀法]
    B -- 孔隙率正常 --> D{设置kn与ks}
    D --> E[设定kn/ks ∈ [2,5]]
    E --> F[计算临界时间步长Δt_cr]
    F --> G[设定Δt ≤ 0.1×Δt_cr]
    G --> H{边界加载速率}
    H -- 应变率>1e-5/s --> I[降低加载速度]
    H -- 合理 --> J[启用局部阻尼]
    J --> K[运行预平衡阶段]
    K --> L{是否稳定?}
    L -- 是 --> M[进入双轴压缩阶段]
    L -- 否 --> N[调整刚度或阻尼]
    N --> K
```

5. 实践案例与调试技巧

某次砂土双轴压缩模拟中，初始设置 kn=1e9 Pa，ks=1e9 Pa（即 kn/ks=1），应变率为 1e-4/s，导致第2000步即出现颗粒穿透。经如下调整后成功收敛：

• 将 ks 调整为 0.3×kn，恢复合理刚度比；
• 降低边界速度至对应应变率 1e-6/s；
• 引入雨落生成方法，将初始孔隙率提升至0.40；
• 设置局部阻尼为0.8，增强能量耗散；
• 验证时间步长满足 Δt ≤ 5e-7 s。

调试过程中建议开启历史记录，监控总动能、势能、不平衡力比（unbal force ratio），确保其在加载前降至1e-4以下。