简介：这份资源面向岩石力学、地下工程与地热开发方向的科研人员和研究生，提供PFC2D5.0颗粒流离散元在热力耦合作用下岩石热损伤模拟的案例代码与解析材料，适合已具备离散元基础、希望深入理解热传导、热膨胀与热应力耦合机制的进阶学习者。压缩包共10个文件，约23KB，以txt说明、html网页和doc文档为主，分别承载案例代码讲解、热力耦合机制梳理与空心圆盘岩样轴心受热损伤的分析思路，便于按主题检索对照。目前已有133人学习下载。读者可从中获取热损伤案例的建模流程、材料属性与边界条件设定方法、温度场与热应力分布的分析视角，以及热损伤萌生扩展的讨论，为地下工程支护设计与施工安全评估提供参考。

1. 从一组“温度升高、岩石开裂”的模拟需求说起

岩石在温度场作用下的损伤演化，是岩土与地质工程里绕不开的一类问题。核废料处置库围岩长期受热、地热开采中的储层岩石反复冷热交替、隧道火灾后衬砌背后岩体的力学响应，这些场景都需要回答同一个问题：温度到底怎么把岩石“搞坏”的。用物理实验做，周期长、成本高，而且高温高压下的内部裂纹扩展过程几乎没法直接观测。PFC2D5.0 的颗粒流离散元方法恰好补上这一环——它把岩石离散成大量颗粒，颗粒之间用接触模型绑定，温度变化引起颗粒膨胀或收缩，进而改变接触力，当接触力超过强度阈值时粘结断裂，裂纹自然萌生、扩展、贯通。这套“热力耦合作用下岩石热损伤”的案例代码，就是围绕这个逻辑搭出来的完整模拟流程，适合做岩石力学、地下工程、地质处置方向毕设或课设的同学直接参考，也适合需要快速验证热损伤趋势的从业者拿去做二次开发。

2. 热力耦合方案设计：从物理场到颗粒接触的映射关系

2.1 热力耦合在 PFC2D 里到底耦合了什么

PFC2D5.0 本身有热力学模块，但热力耦合不是把温度场算完再丢给力学模块那么简单。它的核心逻辑是：温度变化 → 颗粒半径变化 → 接触重叠量变化 → 接触力变化 → 粘结是否断裂。这个链条里有两个关键参数需要提前定清楚。

第一个是颗粒的线性热膨胀系数。PFC2D 里通过thermal.expansion命令给颗粒组赋予热膨胀系数，单位是 1/°C。这个值直接决定同样温升下颗粒半径变化多少。岩石类材料常见取值在 5e-6 到 12e-6 之间，石英含量高的取大值，长石和云母为主的取小值。案例代码里一般会按矿物组成给不同颗粒组分别赋值，而不是全场一个数。

第二个是热传导的时间步。PFC2D 的热求解器是显式差分，时间步受热扩散率控制。如果时间步取太大，温度场会震荡甚至发散；取太小，计算时间成倍增加。通常的做法是先用thermal timestep让系统自动计算一个建议值，再乘以 0.8 作为实际步长。这个 0.8 不是随便定的，是留出安全裕度，防止局部接触变化导致热参数突变时步长失稳。

提示：热力学模块必须在力学模块之前用model configure thermal激活，否则后面所有 thermal 开头的命令都会报错。这个顺序问题在初学阶段翻车率极高。

2.2 颗粒接触模型与热损伤判据的对应关系

案例代码里用的接触模型通常是平行粘结模型（parallel bond），因为它能同时传递力和力矩，适合模拟岩石这种有胶结强度的材料。热损伤的判据就挂在平行粘结的强度参数上：法向强度pb_ten和切向强度pb_coh。当温度引起的接触力超过这两个值中的任何一个，粘结断裂，微裂纹产生。

这里有一个容易混淆的点：平行粘结断裂后，颗粒之间并不是完全分离，而是退化成线性接触模型，仍然可以传递压力。所以热损伤模拟里看到的“裂纹”其实是粘结断裂的计数，不是颗粒飞散。代码里一般用fracture或者自定义 fish 函数统计断裂数量，再按时间步画曲线。

参数标定的思路是：先用单轴压缩模拟标定常温下的pb_ten和pb_coh，让宏观弹性模量和抗压强度对得上实验值；然后保持力学参数不变，只改热膨胀系数和温度边界，观察热裂纹萌生温度是否落在合理区间。如果热裂纹在 50°C 就大量出现，而实验数据显示 200°C 才开始，那就要回头检查热膨胀系数是不是给大了，或者颗粒刚度是不是偏软。

2.3 温度边界与力学边界的施加方式

温度边界一般有两种施加方式：定温边界和热流边界。定温边界是把模型边缘的一层颗粒温度固定，模拟围岩与热源的接触面；热流边界是给边界颗粒一个恒定的热流密度，模拟内部热源。案例代码里常见的是定温边界，因为核废料处置库的场景就是围岩一侧被加热。

力学边界则分两步：先让模型在自重或初始地应力下平衡，用model solve跑到不平衡力比小于 1e-5；然后再激活热力学模块，开始升温。如果顺序反了，先升温再加载，初始应力场会叠加温度应力，导致边界条件不干净，后面分不清哪些裂纹是热引起的、哪些是加载引起的。

代码里通常用一个fish函数记录每一步的断裂数、最大接触力和温度场平均值，输出到 history 里。这样后处理时可以直接画“温度-裂纹数”曲线，判断损伤阈值。

; 定义热力耦合主循环的 fish 函数框架 def thermal_cycle ; 记录当前温度场平均值 local temp_avg = thermal.temp.avg ; 统计平行粘结断裂数量 local frac_num = 0 loop foreach local ct contact.list if contact.model(ct) == 'parallelbond' then if contact.prop(ct,'pb_state') == 0 then frac_num += 1 endif endif endloop ; 写入 history 用于后处理 history(temp_avg, frac_num) end

这段代码的逻辑是每个循环步统计一次温度和断裂数。thermal.temp.avg是 PFC2D 内置的温度平均值函数，contact.prop读取接触的粘结状态，pb_state为 0 表示已断裂。history 记录后可以在 PFC2D 自带的绘图里直接出曲线，也可以导出到 Excel 做进一步处理。参数上唯一需要改的是contact.list的遍历范围，如果模型颗粒数超过 5 万，建议每 10 步统计一次，否则 fish 函数本身会拖慢计算速度。

3. 案例代码实现：从建模到热损伤统计的完整链路

3.1 颗粒生成与初始平衡的代码结构

案例代码的第一段是建模。PFC2D 生成颗粒有两种方式：ball distribute按孔隙率生成，或者ball generate按半径范围生成。热损伤模拟一般用ball distribute，因为要控制初始孔隙率在 0.1 到 0.15 之间，模拟致密岩石。

; 生成颗粒并设置初始孔隙率 model new model configure thermal ball distribute porosity 0.12 radius 0.3 0.5 box -0.05 0.05 -0.1 0.1 ; 设置颗粒密度和热膨胀系数 ball attribute density 2600 damp 0.7 ball attribute thermal.expansion 8e-6 ; 初始平衡 model cycle 1000 model solve ratio 1e-5

这段代码里，porosity 0.12是目标孔隙率，radius 0.3 0.5是颗粒半径范围（单位是米，对应实际岩石的细观尺度）。density 2600是花岗岩类岩石的常见密度。damp 0.7是局部阻尼，用于加速平衡，但注意阻尼只在初始平衡阶段用，热力耦合阶段要关掉或降到 0.1 以下，否则会吸收热膨胀产生的动能，裂纹数偏少。

初始平衡的判据是ratio 1e-5，即最大不平衡力与平均接触力之比小于这个值。如果跑了几万步还不收敛，通常是颗粒半径范围太宽或者孔隙率给得太低，颗粒重叠太大导致初始力就很高。解决办法是把半径范围收窄，或者把孔隙率提高到 0.15 再试。

3.2 热力学参数赋值与温度场初始化

初始平衡完成后，下一步是给热力学参数。PFC2D 的热力学参数包括比热容、热传导系数和热膨胀系数。比热容和热传导系数决定温度场扩散速度，热膨胀系数决定热应变大小。

; 热力学参数赋值 ball attribute thermal.specific_heat 800 ball attribute thermal.conductivity 2.5 ; 初始化温度场 thermal temp initialize 25 ; 设置定温边界：左右两侧颗粒温度固定为 200°C ball group 'left' range x -0.05 -0.045 ball group 'right' range x 0.045 0.05 thermal temp fix 200 range group 'left' thermal temp fix 200 range group 'right'

比热容 800 J/(kg·°C) 和热传导系数 2.5 W/(m·°C) 是花岗岩的典型值。温度初始化 25°C 是室温。定温边界设在左右两侧各 5mm 宽的颗粒组上，温度固定 200°C。这里要注意thermal temp fix和thermal temp initialize的区别：initialize 是给所有颗粒一个初始温度，fix 是把指定颗粒的温度在整个计算过程中锁定，不参与热传导更新。

如果边界温度设得太高，比如直接上 500°C，而热膨胀系数又是 8e-6，颗粒半径变化会非常大，接触力瞬间超过粘结强度，裂纹在第一步就全部断裂，曲线没有渐进过程。通常的做法是从 100°C 开始，分阶段升到目标温度，每阶段稳定后再升下一级。

3.3 热力耦合循环与裂纹统计的实现细节

热力耦合循环是整个代码的核心。PFC2D 里用model cycle或者model solve配合 fish 函数实现。关键是要在每一步同时更新温度场和力学场，而不是交替更新。

; 热力耦合主循环 model mechanical timestep auto model thermal timestep auto ; 分阶段升温 local target_temps = list(100, 150, 200, 250, 300) loop foreach local tt target_temps thermal temp fix tt range group 'left' thermal temp fix tt range group 'right' model solve time 10.0 thermal_cycle endloop

这段代码的逻辑是分五个温度阶段升温，每个阶段计算 10 秒的物理时间，然后调用前面定义的thermal_cycle统计一次。model mechanical timestep auto和model thermal timestep auto让 PFC2D 自动计算力学和热力学时间步，避免手动设置导致的不匹配。

这里有一个血泪经验：热力学时间步和力学时间步的量级差很多。热扩散的时间步通常是力学时间步的 100 到 1000 倍。如果强行用同一个时间步，要么力学计算太慢，要么热力学发散。自动时间步虽然省事，但每次model solve之前最好用model timestep查一下当前值，如果热力学时间步突然掉到 1e-8 以下，说明温度场梯度太大，需要减小升温幅度或者增加边界颗粒数量。

裂纹统计的另一个细节是区分“热裂纹”和“初始裂纹”。初始平衡阶段可能已经有少量粘结断裂，这些不算热损伤。代码里通常在升温前记录一次断裂数作为基准，后面统计时减去这个基准值。

4. 避坑与排查：热力耦合模拟里最容易翻车的五个点

4.1 温度场不收敛，温度值震荡或溢出

现象：运行几十步后温度场出现负值或者超过边界温度，history 曲线剧烈震荡。

原因：热力学时间步超过稳定性条件。PFC2D 的热求解器是显式差分，时间步必须小于颗粒热扩散特征时间的二分之一。如果颗粒半径很小或者热传导系数很大，临界时间步会非常小。

解决：用model thermal timestep auto让系统自动计算，然后手动乘以 0.8。如果仍然震荡，检查颗粒半径是不是小于 0.1mm，这种尺度下热传导系数要相应调低，否则时间步会小到计算无法推进。

4.2 裂纹数在升温第一步就爆表

现象：温度刚升到 100°C，断裂数就从 0 跳到几千，曲线垂直上升。

原因：热膨胀系数给大了，或者颗粒刚度太高。热应变等于热膨胀系数乘以温差，如果热膨胀系数是 1e-5，温差 75°C，应变就是 7.5e-4。对于弹性模量 50GPa 的岩石，热应力就是 37.5MPa，已经接近很多岩石的抗拉强度。如果颗粒刚度再偏大，接触力会更高。

解决：先把热膨胀系数降到 5e-6 试一次，如果裂纹数仍然在第一步就很高，检查平行粘结的法向强度是不是设得太低。通常热损伤模拟里pb_ten要比常温单轴压缩标定值高 20% 左右，因为热应力是三向的，约束比单轴强。

4.3 初始平衡阶段就跑不出收敛

现象：model solve ratio 1e-5跑了几十万步，不平衡力比还在 1e-3 附近晃。

原因：颗粒半径范围太宽，大颗粒和小颗粒之间重叠量差异大，接触力分布极不均匀。或者孔隙率给得太低，颗粒生成时重叠严重。

解决：把半径范围的上限和下限比值控制在 2 以内，比如 0.3 到 0.6。孔隙率不要低于 0.1，如果必须做低孔隙率，先用大半径生成再压缩，而不是直接生成小孔隙率。

4.4 热力学模块和力学模块的阻尼冲突

现象：升温后裂纹数偏少，岩石表现得像橡胶一样吸收能量。

原因：初始平衡阶段设置的局部阻尼没有关掉，热膨胀产生的动能被阻尼吸收，接触力达不到断裂阈值。

解决：在激活热力学模块之前，用ball attribute damp 0.1把阻尼降下来。如果模型仍然偏软，检查是不是用了damp local而不是damp combined，局部阻尼在热力耦合里会过度耗能。

4.5 后处理曲线和实际物理过程对不上

现象：温度-裂纹数曲线是光滑的，但实验数据显示热裂纹是阶梯式萌生的。

原因：升温速率太快，温度场还没均匀就进入下一阶段。或者 history 记录频率太低，漏掉了中间的裂纹萌生事件。

解决：每个温度阶段至少计算到温度场标准差小于 5°C 再进入下一阶段。history 记录频率至少每 100 步一次，如果模型颗粒数少，可以每 10 步一次。另外，裂纹统计要区分法向断裂和切向断裂，法向断裂对应拉裂纹，切向断裂对应剪裂纹，两者在热损伤里的萌生顺序不同。

5. 验收与答辩：怎么证明这套热力耦合模拟是可信的

验收这套模拟，核心不是看代码跑没跑通，而是看结果能不能复现已知的物理规律。我一般会从三个层面检查。

第一层是参数合理性。把热膨胀系数、热传导系数、比热容列成表，和文献里同类岩石的取值范围对比。如果热膨胀系数落在 5e-6 到 12e-6 之外，要么有特殊矿物组成，要么就是标定错了。这个表在答辩时直接放出来，评委一眼就能判断参数有没有离谱。

第二层是趋势验证。温度-裂纹数曲线应该呈现“缓慢萌生—加速扩展—趋于饱和”的三阶段特征。如果曲线是线性的，说明升温速率或者时间步有问题；如果曲线在某个温度点突然垂直上升，说明那个温度下发生了大规模粘结断裂，需要检查是不是边界温度施加方式太激进。通常的做法是把模拟曲线和实验得到的声发射累计曲线放在同一张图里，看拐点温度是否接近。

第三层是敏感性分析。选一个参数（比如热膨胀系数）上下浮动 20%，看裂纹萌生温度变化多少。如果变化超过 30°C，说明模型对这个参数过于敏感，结论的鲁棒性不够。答辩时评委很可能会问“如果热膨胀系数取 6e-6 结果会怎样”，提前做好这个分析就能直接回答。

答辩陈述的逻辑建议按“问题—方法—验证—结论”走，不要一上来就讲代码。先说明热损伤问题的工程背景，再讲为什么选 PFC2D 而不是有限元（有限元处理裂纹萌生需要预设路径，离散元不需要），然后展示参数标定和趋势验证的结果，最后给出热损伤阈值的估计。代码细节放在备用幻灯片里，评委问到再翻。

从那以后我每次做热力耦合模拟，都强制走一遍“参数表—趋势图—敏感性分析”这三步，哪怕时间再紧也不跳过。希望帮到你。

本文还有配套的精品资源，点击获取