简介这份资源面向岩石力学、地下工程与地热开发方向的研究生及科研人员聚焦PFC2D5.0颗粒流离散元在热力耦合条件下岩石热损伤的模拟实现。内容围绕空心圆盘岩样轴心受热这一典型场景涵盖材料属性定义、模型初始化、热载荷施加与结果分析等环节并涉及热传导、热膨胀及热应力与力学行为相互作用的处理思路适合具备一定离散元基础、希望深入理解热损伤萌生与扩展机制的读者参考。资源包共10个文件以txt说明、html网页与doc文档为主分别承载案例解析、技术笔记与背景介绍压缩包约23KB体积轻便便于快速查阅。目前已有133人学习下载。通过案例代码与配套解析读者可掌握温度场演变与热应力响应的记录方法理解热损伤区域的空间分布分析流程并借鉴不同温度梯度、加载速率下的模拟设置为自身课题建模与工程方案评估提供可复用的参考路径。1. 热力耦合岩石热损伤模拟从一份 PFC2D5.0 案例代码说起岩石在温度梯度下会开裂这件事在实验室里用声发射能测到但在数值模型里要复现出来靠的不是把热膨胀系数填进去就完事。我手上这份 PFC2D5.0 的案例代码核心就是解决这个问题用颗粒流离散元把岩石看成几千个可粘结的圆盘先跑热传导算温度场再把温度变化折算成颗粒的膨胀量最后让粘结键在拉剪组合应力下断裂。整个过程不需要连续介质假设裂纹是自发萌生的不是预设的。适合谁用做岩石力学、地热开采、核废料处置库围岩稳定性分析的研究生和工程师尤其是那些已经受够了有限元里“裂纹要预制”这套逻辑的人。这份代码给的是一个可跑通的基准案例不是黑匣子每个 fish 函数都能拆开看。2. 热力耦合在 PFC2D 里到底怎么实现原理与选型2.1 为什么是颗粒流而不是有限元热损伤的本质是细观尺度上矿物颗粒的热膨胀不匹配。有限元里要模拟这个得先假设裂纹在哪里、往哪扩展本质上是在验证一个已知答案。PFC2D 的做法完全不同它把岩石离散成几百到几千个圆形颗粒颗粒之间用平行粘结键连接。温度升高时每个颗粒按热膨胀系数增大半径但周围颗粒约束着它粘结键里就积累法向拉应力和切向剪应力。当应力超过粘结强度键断裂一条微裂纹就出现了。多条微裂纹连通宏观裂纹就形成了。这个路径不需要预设裂纹面所以能捕捉到热损伤最核心的特征——裂纹萌生的随机性和空间分布。选 PFC2D5.0 而不是 6.0 或 7.0原因很实际5.0 的 fish 语言虽然老但热力学模块的接口稳定网上能查到的热力耦合案例大多基于这个版本。6.0 之后命令体系改了很多老代码直接跑会报错。如果你手头只有 6.0也不是不能改但得做好花两天时间调命令格式的准备。2.2 热传导与力学求解的交替流程PFC2D 本身没有内置的热传导求解器这份代码的做法是自己写 fish 函数实现显式热传导。具体流程是先给每个颗粒赋初始温度然后进入热传导循环每个时间步里颗粒从相邻颗粒吸收或释放热量更新温度场。温度场稳定后再进入力学循环根据温度变化量计算颗粒半径膨胀量然后求解力学平衡。两个循环交替进行直到设定的热加载步数走完。这里有个关键参数热传导时间步和力学时间步的比例。如果热传导步太短温度场还没扩散开就进入力学计算结果就是局部热应力过大裂纹会集中在一个区域。常见做法是让热传导步数至少是力学步数的 10 倍具体看热扩散系数和颗粒尺寸。# 热传导与力学交替求解的 fish 伪代码框架 def thermal_mechanical_loop(): # 初始化温度场 for p in particles: p.temp initial_temp # 热加载循环 for thermal_step in range(total_thermal_steps): # 显式热传导颗粒间热交换 for contact in thermal_contacts: q heat_transfer_coeff * (contact.p1.temp - contact.p2.temp) contact.p1.temp - q * dt_thermal / (contact.p1.mass * specific_heat) contact.p2.temp q * dt_thermal / (contact.p2.mass * specific_heat) # 每 N 步热传导后做一次力学平衡 if thermal_step % mechanical_interval 0: # 温度变化引起的半径膨胀 for p in particles: delta_r p.radius * thermal_expansion_coeff * (p.temp - p.ref_temp) p.radius p.initial_radius delta_r # 力学求解 mechanical_solve() # 检查粘结键断裂 for bond in bonds: if bond.normal_stress bond.normal_strength or \ abs(bond.shear_stress) bond.shear_strength: bond.break()这段伪代码里heat_transfer_coeff控制热量传递速率thermal_expansion_coeff是热膨胀系数mechanical_interval决定多久做一次力学平衡。实际代码里这些参数都是从命令行读的改起来不用动核心逻辑。2.3 平行粘结模型的参数怎么定平行粘结模型是 PFC2D 里模拟岩石最常用的接触模型。它把两个颗粒之间的接触看成一根有厚度的梁能传递力和力矩。关键参数有五个法向刚度、切向刚度、法向强度、切向强度、粘结半径。前两个决定变形后三个决定什么时候断。定参数的方法一般是做单轴压缩数值试验反推参数。比如你要模拟花岗岩单轴抗压强度 150 MPa弹性模量 50 GPa那就先给一组初始参数跑单轴压缩看模拟出来的应力应变曲线和实验室数据差多少然后调刚度对弹性模量调强度对峰值应力。这个过程可能要迭代十几轮但一旦定下来同一组参数就能用于热力耦合模拟。注意热膨胀系数不要直接填岩石的宏观值。颗粒的热膨胀系数和粘结键的热膨胀系数可以不同常见做法是让颗粒的热膨胀系数略大于宏观值因为粘结键会约束颗粒膨胀实际表现出来的宏观膨胀量会小一些。3. 从零跑通案例文件结构、命令与参数调试3.1 代码包里的文件都是干什么的这份案例代码解压后大概有七八个文件核心是三个主程序文件、热传导函数库、力学求解函数库。主程序负责读参数、建模型、调循环热传导函数库里有热交换计算和温度更新力学求解函数库里有接触力计算和粘结断裂判断。另外还有两个辅助文件一个用来生成颗粒集合一个用来输出监测数据。常见做法是把所有 fish 函数放在一个文件里但这份代码拆开了好处是改热传导参数不用翻力学代码。如果你要改成自己的模型建议保持这个结构别图省事全塞一起。3.2 建模命令逐行拆解# PFC2D5.0 建模命令流简化版 new set random 10001 # 固定随机种子保证结果可复现 domain extent -0.05 0.05 -0.05 0.05 # 设置计算域单位米 gen ball distribute porosity 0.15 ... # 生成颗粒孔隙率 0.15 radius 0.0005 0.0015 # 颗粒半径范围 0.5-1.5 毫米 cmat default model linear ... # 默认接触模型 kn 1e8 ks 1e8 fric 0.5 # 法向刚度、切向刚度、摩擦系数 prop pb_kn 1e8 pb_ks 1e8 ... # 平行粘结刚度 pb_ten 1e7 pb_coh 1e7 # 粘结抗拉强度、抗剪强度 pb_rad 1.0 # 粘结半径系数 set thermal on # 开启热力学计算set random那行很关键。PFC2D 默认每次生成的颗粒排列都不一样如果不固定随机种子你这次跑出来的裂纹分布和下次跑出来的对不上根本没法做参数分析。domain extent设的是计算域大小一般比试样尺寸大 10% 左右给颗粒膨胀留空间。gen ball distribute里的孔隙率决定了颗粒的密集程度岩石类材料一般取 0.1 到 0.2。平行粘结参数里pb_ten和pb_coh是最容易搞混的。pb_ten是法向强度控制拉伸断裂pb_coh是切向强度控制剪切断裂。热损伤模拟里拉伸断裂通常先发生因为颗粒膨胀时粘结键主要受拉。所以调参数的时候先调pb_ten让拉伸裂纹在合理的温度下出现再调pb_coh控制后续的剪切破坏。3.3 热加载曲线的设置与监测# 热加载曲线设置fish 函数 def apply_thermal_loading(): # 线性升温从 20°C 升到 400°C用时 1000 秒 for step in range(thermal_steps): current_temp 20.0 (400.0 - 20.0) * step / thermal_steps # 只对边界颗粒施加热源 for p in boundary_particles: p.temp current_temp # 内部颗粒通过热传导自然升温 thermal_conduction_step() # 每 100 步记录一次数据 if step % 100 0: record_temperature_field() record_crack_count() record_energy()这段代码做的是线性升温。boundary_particles是边界上的颗粒只给它们赋温度内部颗粒靠热传导慢慢升温。这样做比直接给所有颗粒赋温度更接近实验室的加热条件。record_crack_count()统计的是断裂的粘结键数量这个数随时间的变化曲线就是热损伤演化曲线。record_energy()记录的是系统总能量包括应变能和摩擦能能量突变的地方往往对应裂纹大规模萌生。提示升温速率对结果影响很大。升温太快温度梯度大热应力集中裂纹会提前出现。常见做法是先跑一个慢速升温比如 1°C/s作为基准再跑快速升温做对比。代码里改thermal_steps就能调升温速率。3.4 结果输出与后处理PFC2D5.0 的输出格式比较原始一般是文本文件每行记录一个时间步的数据。后处理可以用 Python 读进来画图。裂纹分布图需要导出颗粒位置和粘结状态然后用 matplotlib 画散点图断裂的粘结用红色标出。# 用 Python 读取 PFC2D 输出并画裂纹分布 import numpy as np import matplotlib.pyplot as plt # 读取颗粒位置 data np.loadtxt(particle_positions.txt) x, y data[:, 0], data[:, 1] # 读取断裂粘结 bonds np.loadtxt(broken_bonds.txt) bx, by bonds[:, 0], bonds[:, 1] plt.figure(figsize(8, 8)) plt.scatter(x, y, s1, cgray, alpha0.5, label颗粒) plt.scatter(bx, by, s5, cred, markerx, label断裂粘结) plt.axis(equal) plt.legend() plt.title(热损伤裂纹分布) plt.savefig(thermal_crack_pattern.png, dpi300)这段代码里particle_positions.txt和broken_bonds.txt是 PFC2D 导出的。s1控制颗粒点的大小s5控制裂纹点的大小。画出来的图能直观看到裂纹是均匀分布还是集中在某个区域。如果裂纹全集中在边界说明边界热源太强得降低升温速率或者增大热传导系数。4. 避坑指南热力耦合模拟里最容易翻车的五个地方4.1 温度场震荡不收敛现象热传导循环跑了几百步温度场还在上下震荡没有稳定的趋势。原因热传导时间步太大超过了显式求解的稳定性条件。显式热传导的稳定时间步和颗粒尺寸、热扩散系数有关颗粒越小稳定时间步越小。解决把dt_thermal减小到原来的十分之一或者增大颗粒尺寸。如果不想改颗粒尺寸就得接受更长的计算时间。4.2 裂纹还没升温就出现了现象初始平衡阶段还没开始升温就有一批粘结键断了。原因初始颗粒排列太密重叠量太大初始接触力过高。PFC2D 生成颗粒时如果孔隙率设得太低颗粒之间会大量重叠初始应力就很大。解决把孔隙率从 0.15 调到 0.2或者用cycle命令让系统先平衡几千步再升温。常见做法是初始平衡阶段跑 5000 步以上确认最大接触力小于粘结强度的 10%。4.3 热膨胀系数填错导致裂纹数量不对现象升温到 400°C裂纹数量只有几十条明显偏少。原因热膨胀系数填的是宏观值比如花岗岩的 8e-6 /°C但颗粒尺度上这个值偏小。解决把颗粒的热膨胀系数调到 1.5 到 2 倍宏观值同时把粘结键的热膨胀系数设为 0。这样颗粒膨胀时粘结键不膨胀拉应力更集中裂纹更容易萌生。4.4 计算速度慢到无法接受现象跑 1000 个热加载步要十几个小时。原因颗粒数量太多或者力学求解的收敛精度设得太高。解决如果颗粒数量超过 5000考虑用 PFC2D 的并行计算功能或者把二维模型简化成更小的尺寸。力学求解的收敛精度一般设 1e-5 就够了设 1e-8 纯属浪费时间。4.5 重启后结果对不上现象中断后从保存点重启跑出来的裂纹分布和之前不一样。原因PFC2D5.0 的保存文件不包含随机数生成器的状态重启后随机种子变了。解决在保存之前用set random重新固定种子或者在 fish 里手动记录随机数状态。更稳妥的做法是每次跑完整流程别中断。5. 进阶技巧用热损伤曲线反推细观参数跑通案例只是第一步真正有用的是用模拟结果反推细观参数。具体做法是先跑一组不同pb_ten值的模拟记录每个值对应的热损伤起始温度裂纹数量开始明显增加的温度。然后画pb_ten对起始温度的关系曲线。如果实验室测得的起始温度是 250°C就在曲线上找对应的pb_ten值。这个过程叫参数标定比瞎试快得多。# 参数标定不同 pb_ten 对应的热损伤起始温度 pb_ten_values [5e6, 1e7, 2e7, 5e7] # 四组抗拉强度 onset_temps [] for pb_ten in pb_ten_values: set_pb_ten(pb_ten) run_thermal_simulation() onset_temp detect_crack_onset() # 检测裂纹起始温度 onset_temps.append(onset_temp) # 画标定曲线 plt.plot(pb_ten_values, onset_temps, o-) plt.xlabel(粘结抗拉强度 (Pa)) plt.ylabel(热损伤起始温度 (°C)) plt.xscale(log) plt.grid(True) plt.savefig(calibration_curve.png, dpi300)这段代码里detect_crack_onset()的逻辑是当裂纹数量超过总粘结数的 1% 时认为热损伤开始。这个阈值可以调但别设太低否则初始平衡阶段的噪声会被误判成裂纹。标定曲线画出来之后如果实验室数据落在曲线附近说明参数合理如果偏离太远就得检查热膨胀系数或者热传导系数是不是设错了。还有一个技巧用声发射模拟验证裂纹演化。PFC2D 里每个粘结断裂都可以记录断裂时刻和断裂位置把这些数据按时间排序就得到了模拟声发射事件序列。和实验室的声发射数据对比如果 b 值震级-频度关系的斜率接近说明裂纹演化模式相似。这个验证比单纯对比裂纹数量更有说服力。从那以后我每次跑热力耦合模拟都强制走一遍参数标定流程哪怕只是改了一个热膨胀系数。因为热损伤对参数太敏感了不标定根本不知道结果可不可信。希望帮到你。本文还有配套的精品资源点击获取