温度进去了，力为什么没出来——PFC 6.0 里 PB/CB 模型热力耦合的双球实测

这篇文章记录一次针对 PFC 6.0 热力耦合的定向验证：三个待证论断、16 组固定双球算例、4 组自由双球对照，外加过程中挖出来的四个"不报错的坑"。其中第四个坑我从 C++ 源码里找到了根因，结论可能和你现在的认知不太一样。

1. 一、要验证的三条论断

热力耦合的链条是这样的：温度场 → 颗粒热膨胀（半径变大）→ 接触重叠改变 → 接触力改变 → 应力重分布 → 胶结破坏。

这条链上，胶结模型（PB 平行胶结 / CB 接触胶结）是最常被用在岩石类材料上的。所以问题很具体：

• 论断一：PB 模型的线性部分必须是全量模式，粒径变化才能转成接触力；
• 论断二：加温时颗粒之间必须有重叠，热才对胶结材料起作用；
• 论断三：或者，CB 模型在全量模式下也能有作用。

这三条都是"经验判断"，值得用最干净的算例钉死。而最干净的算例就是——两颗球。

2. 二、算例怎么搭：把"测热应力"变成"测接触力"

两颗半径 R = 0.5 m 的球，球心距 d = 2R + gs0，gs0 < 0 就是重叠量。关键设计是把两颗球的速度全部固定——完全约束。

参数（全部取整数值，方便手算复核）：

• 线性刚度 kₙ = kₛ = 1e9 N/m
• 平行键刚度 k̄ₙ·A = 1e9 N/m（A 为 2D 键截面积 2R̄，取 pb_rmul=1）
• 热膨胀系数 α = 3e-6 /°C，升温 ΔT = 100 °C
• 摩擦 0.5，键强度取 1e30（保证升温过程中键不断，先只看力学响应）

对应的闭式值：

这三个数就是后面的判据。算例脚本由 AI 自动生成（16 组参数矩阵一次跑完），全程零报错。

3. 三、结果：几何变了，力没变

图 1 是 16 组算例里 15 组（另一组没形成接触，见第六节）的线性力变化。红色虚线是理论值 3.0e5 N。

 

图 1 固定双球算例的线性部分法向力变化。横轴是力，红色虚线为闭式解 kₙ·2αRΔT = 3.0e5 N；深蓝为 PB、绿色为 CB、紫色为无胶结对照组。

只要球的膨胀系数不为零，每一组算例里球的半径都长了 150 μm、间隙都压紧了 300 μm、键的状态都还是"完好"（state=3）。 也就是说，几何层面温度确实起作用了。（唯一的例外是那个把球 α 设成 0 的对照例，它的半径当然不变，力自然也是 0。）

但力只有一半的算例有反应。

最直接的对照是前两组：

几何变化逐位相同，力一个 3.0e5、一个 0。 而且 3.0e5 和闭式解 kₙ·2αRΔT 的偏差是 0 位有效数字——不是"接近"，是完全相等。

换成 CB、换成没有胶结的纯 linear 模型，规律一模一样：全量模式 3.0e5，增量模式 0。这说明这条规律跟有没有胶结无关，是线性部分自己的事。

再看两组细节：

• "PB lin_mode 未设置"那一行，结果和 lin_mode=0 完全一致——所以 PFC 6.0 的默认值就是全量模式。但正因为默认是对的，一旦有人手动把它改成 1，症状会极其隐蔽：模型照跑、温度照升、半径照长、什么都不报错，只是力不出来了。
• "PB 全量 + thexp 已设"和"PB thexp 已设、球 α=0"这两行：前者 3.0e5，后者 0。这一对直接指向第五节的第四个坑。

4. 四、自由双球：更刺眼的一组

固定双球说明的是"力出不来"。把球放开、让它们自己找平衡，能看到更直观的版本——膨胀量本身也没了。

图 2 是四组自由双球算例的球心距变化。

 

图 2 自由双球（无外载）的球心距变化。理论自由热膨胀量 2αRΔT = 300 μm 用红色虚线标出。深蓝为 PB、绿色为 CB。

• CB + 全量：球心距恰好增加 300 μm，正好等于 2αRΔT，同时接触力回到 0。这是教科书式的自由热膨胀。
• CB + 增量：同一个模型，只把 lin_mode 从 0 改成 1，球心距变化 0。半径长了 150 μm，两颗球纹丝不动。
• PB + 全量：球心距增加 150 μm。比 CB 少一半，因为平行键把位移约束掉了一部分（平衡时线性受压 6.5e5 N、键受拉 −6.5e5 N，两者相消）。
• PB + 增量：同样是 0。

同一套参数、同一个模型，唯一的差别是 lin_mode。结果一个动、一个完全不动。

5. 五、四个坑（这部分比结论值钱）

5.1 坑 1：热步只在力学进程关闭时才会推进

这是最开始卡了我最久的地方。按直觉写：

bashmodel thermal on
model mechanical on
ball thermal attribute temperature-increment=100.0
model cycle 1

结果什么都不发生：半径不变、温度不变、力不变，而且没有任何报错。

原因是：热力与力学两个进程同时开着的时候，model cycle 只驱动力学进程，热步根本不走。正确写法是把力学关掉，让热步单独跑：

bashmodel thermal on
model mechanical off                      ; ← 关键
ball thermal attribute temperature 0.0
ball thermal attribute temperature-increment=100.0
model cycle 1                             ; 这一步才真正升温、改半径
model thermal off
model mechanical on
model cycle 1                             ; 再走一步力学，刷新线性力

model mechanical slave on 也能跑通。官方 3D 验证算例 datafiles3d/thermal/free_expansion 和 constrained_expansion 用的就是这套流程。

判据：跑完 model cycle 后，如果日志里只打印了力学时间步、没有热时间步和热时间，那热步就没走。

5.2 坑 2：球全固定，自动时间步算不出来

两颗球都 ball fix 之后，模型里没有任何自由体，PFC 直接报：

bash*** Unable to determine a valid timestep.

解决很简单：model mechanical timestep fix 1.0e-8。球本来就不动，时间步取多少都稳定。

5.3 坑 3：lin_mode 的语义和默认值

官方属性表写得很清楚，但很容易记反：

• lin_mode = 0：全量（绝对）更新 ← 默认值
• lin_mode = 1：增量更新

官方文档在属性表下面还有一句关键的注意事项：

When M_l = 0, the normal component of the linear force is automatically overridden during the next force-displacement calculation.

翻译过来就是：全量模式下，每一步都按当前几何（包括热膨胀之后的半径）重新计算法向力；增量模式只累加相对位移增量，而热膨胀不是相对位移增量，所以被完全忽略。

这就是"半径长了、力没出来"的全部原因。

5.4 坑 4：PB 的"键材料热膨胀"在 2D 球-球接触下恒等于 0

这一条最出乎意料。

PFC 文档在热力模块里专门写了一节：平行胶结模型的键材料本身也会热膨胀，公式是

bashΔF̄ₙ = −k̄ₙ·A·ᾱ·L̄·ΔT

其中 ᾱ 取两端颗粒膨胀系数的平均值，L̄ 是两球心距。按这个公式，上面那组算例里键力应该变化 3e5 N 左右。

实测：一点变化都没有。 我把球的 α 设成 3e-6、把热管（ThermalPipe）的 thexp 属性设成 3e-6（两种设置写法都试了），键力始终是 0。反过来把球的 α 设成 0、只留 thexp，结果也是 0——说明这条路径根本没被触发。

我去翻了安装目录里的接触模型 C++ 源码（ContactModelLinearPBond::thermalCoupling），原因在这里：

cppif (pb_state_ < 3) return;                       // 必须成键
int idx = ct->getModel()->getContactModel()->isProperty("thexp");
if (idx <= 0) return;                            // 热管必须有 thexp
double delUn = ts->length_ * thexp * ts->tempInc_;
if (delUn == 0.0) return;
double dthick = 0.0;
double Cmin1  = ms->end1Curvature_.x();          // 球：1/R ≠ 0
if (Cmin1 == 0.0) dthick = ms->end1Extent_.x();  // 只有零曲率的平面才给厚度
double br     = pb_rmul_ / max(Cmax1, Cmax2);
double pbArea = dthick <= 0.0 ? dPi*br*br : 2.0*br*dthick;
finc *= pb_kn_ * pbArea * delUn;                 // dthick=0 ⇒ pbArea=0 ⇒ 增量恒为 0

关键在 dthick：它只在端部曲率为零（也就是平面，比如 rblock 那种块体）的时候才被赋值。球-球接触的曲率是 1/R，不等于零，于是 dthick 保持 0，pbArea 算出来是 0，力增量恒等于 0。

（dthick 那个分支是给块体/rblock 留的，我这次没有在 rblock 上验证，先按源码这么记。）

6. 六、"必须有重叠"的精确版

论断二在实践中是对的，但精确表述要修正一下。

contact method bond 的默认成键区间是 gap ∈ (−∞, 0]，也就是只给重叠或相切的接触成键。实测三组：

• gs0 = +1e-3（有间隙）+ 默认成键：连接触对象都不存在，contact.list 是空的。没有接触，自然没有胶结，也谈不上热响应。
• gs0 = +1e-3 + contact method bond gap 1e-3（强行跨间隙成键）：接触和键都建起来了，但升温后线性力仍然是 0。因为间隙只从 +1e-3 压到 +7e-4，还没闭合；线性模型不能受拉，所以不出力。键那一路又是坑 4 的恒零。
• gs0 = 0（相切）：只要接触被检出，工作正常，全量模式同样给出 3.0e5 N。

顺便说一句相切：本次两颗球的算例里，gs0 = 0 恰好被检出了接触。但相切球对的接触检出本身不可靠（规则排列下实测约有一半漏检），工程上还是老老实实用受控重叠更稳。

7. 七、给工程使用的检查清单

1. 用 PB/CB 做热力耦合前，确认 lin_mode = 0。 本机默认就是 0，但别去改成 1——改了热就彻底不起作用，而且不报错。
2. 别指望 thexp。 2D 球-球 PB 的键材料热膨胀项算出来就是 0，热响应只能来自线性部分。
3. 升温用"热步（mechanical off）+ 力学步"交替。 一次性施加很大的 ΔT 会让半径突变、产生冲击，分步升温更稳。
4. 测热应力要测内部截面或内部球对。 自由边界最外层的球会过冲。
5. 建接触之前先把 contact cmat default 和 contact cmat thermal default 设好；specific-heat 必须在任何 model cycle 之前给，否则报 Specific Heat is 0。

8. 八、这次工作里 AI 做了什么

值得单独说一句：上面这套东西——16 组参数矩阵的脚本生成、无头运行、结果解析成表、出图、以及那四个坑的定位——都是 AI 自动完成的。人给的需求只有一句话。

四个坑里有三个是"不报错的错误"：热步不推进、增量模式忽略膨胀、thexp 恒零。这类问题靠"程序正常结束"是判不出来的，必须挂客观判据——这次用的是闭式解 kₙ·2αRΔT 和 2αRΔT。判据对得上，说明模型对；对不上，就得往下挖。

坑 4 最后是靠读安装目录里的 C++ 源码定案的。这也是我觉得 AI 做这类活儿的价值所在：它不会嫌翻源码麻烦。

参考资料

1. Itasca Consulting Group. PFC 6.0 Documentation — Thermal Option（Thermal Strains / ThermalPipe contact model）与 Linear Parallel Bond Model 属性表。
2. PFC 6.0 安装目录内置接触模型源码 ContactModelLinearPBond::thermalCoupling。
3. 本文全部算例脚本、16 组结果表与后处理代码，由 Claude Code 在 PFC 6.0 环境自动生成并验证。