简介：这份资源面向岩石力学与地下工程方向的科研人员和研究生，聚焦PFC2D5.0颗粒流离散元在热力耦合条件下的岩石热损伤模拟，重点解析空心圆盘岩样轴心受热损伤的案例代码。内容围绕材料属性定义、模型初始化、热载荷施加与结果分析等环节展开，涉及热传导、热膨胀及热应力与力学行为的耦合关系，并讨论热损伤区域的萌生与扩展规律，可服务于深部岩体工程与地热开发等场景的稳定性评估。资源包共10个文件，以txt说明、html网页和doc文档为主，分别承载案例解析、技术笔记与背景介绍，压缩包约23KB，体积轻便便于快速查阅。目前已有130人学习下载，适合希望借助案例代码理解热力耦合建模流程、补充离散元模拟思路的读者参考。
1. 从一次“温度一加，岩样就碎成渣”的模拟翻车说起
去年帮某实验室复核一组岩石热损伤的数值实验，对方用 PFC2D5.0 跑热力耦合，结果温度刚加载到 200℃，试样就整体崩解，颗粒飞得到处都是，跟真实花岗岩在热冲击下的渐进式开裂完全对不上。翻他们代码一看，热膨胀系数直接按钢的数值填，颗粒接触模型还是默认的线性刚度，热学模块和力学模块各跑各的，中间连个耦合时间步都没对齐。这套 PFC2D5.0 离散元热力耦合岩石热损伤案例代码，核心就是解决这类“热学力学两张皮”的问题——它把温度场演化、热膨胀应变、颗粒间接触力重分布串成一条链，用空心圆颗粒体系模拟岩石在温度梯度下的微裂纹萌生与扩展。适合做岩石力学、地热开采、核废料处置库围岩稳定性分析的从业者，尤其是那些已经在用 PFC2D 但热学模块一直跑不通的人。代码包里的空心圆方案不是随便选的，后面会讲清楚它跟实心圆在热传导和破裂模式上的本质差别。
2. 热力耦合在 PFC2D 里到底怎么串起来：从温度场到接触力链
2.1 热学模块与力学模块的耦合逻辑
PFC2D5.0 本身有热学选项，但很多人只把它当个“温度赋值器”用，这是最大的误解。热力耦合的本质是：温度变化 → 颗粒热膨胀（或收缩）→ 颗粒半径变化 → 接触重叠量变化 → 接触力变化 → 颗粒运动 → 孔隙率、热传导路径变化 → 温度场再分布。这个闭环里，时间步必须统一，热学时间步和力学时间步不能各走各的。常见做法是设置一个耦合步长dt_couple，每走一个力学步就更新一次温度场，或者反过来，但更稳的方式是用热学时间步去限制力学时间步，因为热扩散系数通常比力学波速慢好几个量级。
代码里用model configure thermal开启热学模块，然后通过thermal time-step设置热学时间步。关键参数是thermal conductivity（热导率）和specific heat（比热容），这两个值直接决定温度场扩散速度。岩石类材料常见取值：热导率 2.5～3.5 W/(m·K)，比热容 800～900 J/(kg·K)，热膨胀系数 8e-6～12e-6 /℃。这些值在代码里是以颗粒属性方式赋的，不是全局常数，因为不同矿物颗粒的热学参数差异很大。
; 开启热学模块并设置耦合时间步 model configure thermal thermal time-step 1e-4 ; 定义热学属性组：石英、长石、云母 thermal property conductivity 3.0 specific-heat 850.0 thermal property conductivity 2.2 specific-heat 780.0 thermal property conductivity 1.8 specific-heat 880.0 ; 热膨胀系数按颗粒分组赋值 fish define set_alpha loop foreach local p list if ball.group(p) == 'quartz' ball.extra(p,'alpha') = 1.0e-5 else ball.extra(p,'alpha') = 8.0e-6 endif endloop end逻辑说明：thermal time-step设得太大，温度场会震荡甚至发散；设得太小，计算时间成倍增加。我一般先用1e-4试跑 100 步，看温度场是否平滑，再决定是否放大。ball.extra是 PFC 里存自定义属性的常用手段，比全局变量安全，因为每个颗粒独立。热膨胀系数没有内置属性，必须自己用 FISH 存。
2.2 空心圆颗粒的热传导优势与接触模型选择
为什么这个案例用空心圆而不是实心圆？实心圆颗粒在热传导时，热量从表面传到中心需要的时间跟半径平方成正比，大颗粒内部温度滞后严重，导致热应力分布失真。空心圆把内部挖空，等效热传导路径变短，温度场响应更快，而且空心部分可以模拟岩石中的孔隙或弱胶结面，热裂纹更容易从这些位置萌生。代码里空心圆的内外径比一般取 0.6～0.8，太小了力学强度掉得厉害，太大了热传导优势不明显。
接触模型方面，热力耦合下不能用简单的线性接触，因为温度变化会改变接触刚度。常见做法是用linear模型加thermal膨胀，或者用parallel-bond模型，后者能模拟胶结破坏，更适合岩石。代码里用的是平行粘结模型，粘结半径系数设 0.8，粘结刚度比颗粒刚度略低，这样热膨胀时粘结先破，形成微裂纹。
; 空心圆颗粒生成与平行粘结设置 ball generate radius 0.5 0.8 ... id 1 2000 ... group 'quartz' ; 设置空心：用 FISH 修改颗粒密度和热学属性 fish define make_hollow loop foreach local b ball.list local r0 = ball.radius(b) ball.radius(b) = r0 * 0.8 ; 外径不变，内径通过密度折减模拟 ball.density(b) = ball.density(b) * 0.6 endloop end ; 平行粘结 contact model linearpbond contact method bond contact property pb_kn 1e8 pb_ks 1e8 pb_ten 5e6 pb_coh 1e7 pb_fa 30参数说明：pb_kn和pb_ks是粘结法向和切向刚度，岩石类一般取 1e8～1e9 N/m。pb_ten是抗拉强度，pb_coh是粘结力，这两个值决定热裂纹萌生的温度阈值。pb_fa是摩擦角，影响破裂后的残余强度。空心圆的密度折减系数 0.6 是经验值，对应内径约 0.63 倍外径，实际建模时最好用ball.radius直接改，但 PFC2D 不支持真正的空心几何，只能用密度和热学属性等效。
2.3 温度加载路径与边界条件设置
温度加载不能直接给所有颗粒赋一个高温，那样没有温度梯度，热应力出不来。代码里用的是边界温度斜坡加载：左侧边界固定 20℃，右侧边界以 5℃/min 的速率升到 300℃，中间颗粒通过热传导自然升温。这样温度梯度从右向左递减，热膨胀应变也递减，产生的热应力在中间区域最大，裂纹从右边界附近萌生然后向左扩展。
; 温度边界条件：右侧升温，左侧恒温 fish define apply_thermal_bc local t_current = thermal.time local t_right = 20.0 + 5.0 * t_current / 60.0 if t_right > 300.0 t_right = 300.0 endif loop foreach local b ball.list local x = ball.pos.x(b) if x > 0.9 * domain_max_x ball.temp(b) = t_right endif if x < 0.1 * domain_max_x ball.temp(b) = 20.0 endif endloop end ; 每 100 步调用一次 model solve time 0.01逻辑说明：thermal.time是热学模块的内部时钟，跟力学时间可能不同步，所以要用model solve time控制总时长。边界颗粒的温度直接赋值，内部颗粒靠热传导自动更新。注意ball.temp是内置属性，但只有在model configure thermal之后才有效。升温速率 5℃/min 是实验室常用值，太快了热冲击效应太强，太慢了计算时间太长。
3. 从零跑通案例：建模、加载、监测的完整操作链
3.1 颗粒生成与热学属性赋值
第一步是生成颗粒体系。代码里用的是ball distribute按高斯分布生成，半径范围 0.5～0.8 mm，孔隙率 0.15，对应中密砂岩。生成之后要压实，用model solve ratio-average 1e-5让体系平衡。压实之后才能赋热学属性，否则颗粒位置还在变，热学参数赋了也白赋。
; 颗粒生成与压实 model new model configure thermal ball distribute porosity 0.15 radius 0.5 0.8 ... box -10 10 -10 10 ball attribute density 2600 damp 0.7 contact model linearpbond contact method bond contact property pb_kn 1e8 pb_ks 1e8 ... pb_ten 5e6 pb_coh 1e7 pb_fa 30 model solve ratio-average 1e-5 ; 热学属性赋值 thermal property conductivity 3.0 specific-heat 850.0 fish define set_thermal_alpha loop foreach local b ball.list ball.extra(b,'alpha') = 1.0e-5 endloop end参数说明：porosity 0.15是目标孔隙率，实际生成后可能略有偏差，用ball porosity命令查看。damp 0.7是局部阻尼，压实阶段用大一点加快收敛，后续加载阶段要降到 0.1 以下，否则会抑制热裂纹扩展。ratio-average 1e-5是平衡判据，值越小越平衡，但计算时间越长，一般 1e-5 够用。
3.2 热力耦合加载与监测变量设置
加载阶段要同时监测力学和热学变量。力学监测：颗粒位移、接触力链、裂纹数量。热学监测：温度场分布、热流密度。代码里用history命令记录这些变量，每 10 步存一次，方便后续画曲线。
; 监测变量设置 history interval 10 history thermal time history mechanical ratio-average history ball temp average history contact force average ; 裂纹监测：平行粘结破坏数 fish define count_broken_bonds local count = 0 loop foreach local c contact.list if contact.model(c) == 'linearpbond' if contact.prop(c,'pb_state') == 3 count += 1 endif endif endloop count_broken_bonds = count end history fish count_broken_bonds ; 加载 model solve time 3600.0逻辑说明：history interval 10表示每 10 个力学步记录一次，太密了文件大，太疏了曲线不光滑。pb_state == 3表示粘结已经破坏，这是 PFC 内部状态码，不同版本可能不同，5.0 里 3 是破坏。model solve time 3600.0对应 1 小时物理时间，按 5℃/min 升温速率，温度从 20℃ 升到 320℃，但边界限制在 300℃，所以实际到 300℃ 后恒温。
3.3 结果提取与热损伤云图生成
跑完之后要提取结果。PFC2D 自带plot命令可以画云图，但热损伤云图需要自己算损伤因子。常见做法是用裂纹密度：单位面积内的破坏粘结数。代码里用 FISH 遍历所有破坏粘结，统计每个颗粒周围的裂纹数，然后映射到颜色。
; 热损伤云图：裂纹密度 fish define damage_factor loop foreach local b ball.list local crack_count = 0 loop foreach local c contact.list if contact.end1(c) == b if contact.prop(c,'pb_state') == 3 crack_count += 1 endif endif endloop ball.extra(b,'damage') = crack_count endloop end ; 输出到文件 fish define output_damage local file = file.open('damage.txt',1,1) loop foreach local b ball.list file.write(file,string(ball.pos.x(b)) + ' ' + ... string(ball.pos.y(b)) + ' ' + ... string(ball.extra(b,'damage'))) endloop file.close(file) end参数说明：contact.end1(c)返回接触的第一个颗粒，遍历所有接触就能统计每个颗粒的裂纹数。ball.extra(b,'damage')存损伤因子，后续可以用 Python 或 Tecplot 画云图。输出文件用空格分隔，方便导入。注意file.open的第二个参数 1 表示写模式，第三个参数 1 表示覆盖，别搞反了。
4. 避坑指南：热力耦合模拟里最容易翻车的五个地方
4.1 温度场震荡发散
现象：温度场跑着跑着出现负值或者剧烈震荡，计算直接中断。原因：热学时间步太大，超过了热扩散稳定性条件。热扩散稳定性要求dt < dx^2 / (2 * alpha)，其中dx是颗粒间距，alpha是热扩散系数。颗粒间距 0.5 mm，热扩散系数 1e-6 m²/s，算下来dt要小于 1.25e-4 s。解决：把thermal time-step降到 1e-5 以下，或者用model solve time控制总时长，别让热学步长自动放大。
4.2 热膨胀系数单位搞错
现象：温度升到 100℃，颗粒膨胀量肉眼可见，试样直接炸开。原因：热膨胀系数填成了 1e-3 级别，那是塑料的数值，岩石是 1e-5 级别。解决：填之前先查文献，花岗岩 8e-6，砂岩 1.2e-5，石灰岩 6e-6，别凭感觉。代码里用ball.extra存，赋值时检查数量级。
4.3 空心圆密度折减导致力学参数失真
现象：空心圆方案跑出来的弹性模量比实心圆低一半，跟实验对不上。原因：密度折减了，但接触刚度没跟着调，动力学响应变了。解决：密度折减后，接触刚度要按质量比重新标定。常见做法是pb_kn乘以密度折减系数的平方根，或者直接用ball.radius改小再重新生成接触。
4.4 边界温度直接赋值导致热冲击
现象：右侧边界颗粒温度瞬间从 20℃ 跳到 300℃，边界处裂纹密得像蜘蛛网，内部却没动静。原因：直接赋值没有升温过程，热冲击效应远大于真实情况。解决：用斜坡加载，每步增加 0.1℃，或者用thermal boundary命令设置对流换热系数，让边界温度渐变。
4.5 监测频率太高导致文件爆炸
现象：跑完发现 history 文件几十个 G，打都打不开。原因：history interval设成了 1，每步都记。解决：设成 10 或者 100，看总步数定。总步数 1e6 的话，interval 100 也有 1e4 个点，画曲线够了。另外history只记需要的变量，别把ball.pos全记下来。
5. 进阶技巧：用 FISH 做热损伤各向异性统计与参数反演
跑通基础案例之后，真正有价值的是从结果里挖出各向异性损伤信息。岩石热损伤不是均匀的，温度梯度方向跟裂纹扩展方向有固定关系。我一般用 FISH 写个各向异性统计函数，把裂纹按角度分到 12 个扇区，每个扇区统计裂纹密度，然后画玫瑰图。这样能看出热裂纹是沿温度梯度方向优势扩展，还是垂直于梯度方向。
; 各向异性损伤统计：12 扇区 fish define damage_anisotropy local sector_count(12) local sector_angle = 30.0 loop foreach local c contact.list if contact.prop(c,'pb_state') == 3 local dx = contact.pos.x(c) - 0.0 local dy = contact.pos.y(c) - 0.0 local angle = math.atan2(dy,dx) * 180.0 / math.pi if angle < 0.0 angle += 360.0 endif local idx = int(angle / sector_angle) + 1 if idx > 12 idx = 12 endif sector_count(idx) += 1 endif endloop local total = 0 loop local i = 1 to 12 total += sector_count(i) endloop loop local i = 1 to 12 local ratio = sector_count(i) / total io.out('Sector ' + string(i) + ': ' + string(ratio)) endloop end参数说明：math.atan2返回弧度，转成角度后分扇区。sector_angle是每个扇区的角度范围，12 个扇区就是 30°。io.out输出到控制台，也可以写文件。这个统计跑一次大概几秒钟，但能看出很多信息。比如温度梯度沿 x 方向时，如果扇区 1 和 7（0° 和 180°）的裂纹密度明显高，说明裂纹沿梯度方向扩展；如果扇区 4 和 10（90° 和 270°）高，说明垂直梯度方向扩展。
参数反演是另一个进阶方向。实验室做热损伤实验，能测到的是宏观弹性模量随温度的变化曲线，但不知道微观裂纹怎么长。可以用 PFC2D 做反演：先假设一组微观参数（粘结强度、热膨胀系数、热导率），跑模拟得到宏观模量-温度曲线，跟实验对比，用优化算法调微观参数。常见做法是网格搜索或者遗传算法，但 PFC2D 跑一次几分钟到几小时，优化效率很低。我一般先用少量样本建个代理模型，比如用多项式响应面，然后在代理模型上做优化，最后用 PFC2D 验证最优解。这样能把反演周期从几周缩到几天。
从那以后我每次跑热力耦合之前，都强制走一遍检查清单：热学时间步算过没有、热膨胀系数数量级对不对、边界温度是不是斜坡、监测频率会不会爆文件、空心圆密度折减后刚度调没调。这五条过一遍，基本不会翻车。希望帮到你。
本文还有配套的精品资源，点击获取
