简介这份资源面向岩石力学、地下工程与地热开发方向的科研人员和研究生提供PFC2D5.0颗粒流离散元在热力耦合作用下岩石热损伤的案例代码与解析材料适合已具备离散元基础、希望深入理解热传导、热膨胀与热应力耦合机制的中高级学习者。压缩包共10个文件约23KB包含4个txt说明、4个html网页文档与2个doc分析文档分别用于代码注释、理论梳理与案例讲解便于按模块查阅。内容围绕空心圆盘岩样轴心受热损伤这一典型场景涵盖材料属性定义、模型初始化、热载荷施加及温度场与热应力响应分析等环节并延伸至热损伤萌生与扩展的空间分布讨论。已有130人学习下载可帮助读者快速掌握热力耦合建模思路、理解岩石热损伤演化规律并为支护设计与施工安全评估提供参考。1. 热力耦合下的岩石热损伤为什么值得用 PFC2D5.0 颗粒流重做一遍岩石在温度梯度下开裂这件事在实验室里不好直接看。你只能看到最终裂纹形态看不到颗粒之间什么时候开始脱开、热应力在哪一步超过胶结强度、微裂纹怎么从均匀分布突然聚集成宏观裂缝。PFC2D5.0 的颗粒流离散元把岩石打散成一颗颗圆盘胶结模型负责把它们粘住热力学模块负责让颗粒热胀冷缩两者一耦合热损伤的全过程就能在时间步里逐帧回放。这套思路适合做岩石力学、地热开采、核废料处置库围岩稳定性、高温岩体工程的研究人员和工程师。空心圆在这里不是装饰它是颗粒生成阶段控制孔隙率与粒径分布的关键几何参数也是热力耦合案例里最容易被忽略却直接影响初始应力场的一环。下面按“先立住理论、再跑通代码、最后避开坑”的顺序拆开讲。2. 热力耦合在 PFC2D5.0 里到底耦合了什么从热传导到胶结断裂2.1 温度场与力学场的双向传递路径PFC2D5.0 的热力耦合不是把两个物理场简单叠在一起。它的基本流程是热学模块先根据颗粒间的热阻网络求解温度场每个颗粒获得一个温度值然后热膨胀系数把温度变化转换成颗粒半径的增量半径变化挤压相邻颗粒通过接触模型产生接触力接触力再反馈到力学求解器更新颗粒位置和接触状态如果胶结键的应力超过强度键断裂产生微裂纹。整个过程在每个时间步内交替执行热学步和力学步可以设置不同的步长比例。这里的关键参数是热膨胀系数alpha、热阻thermal resistance和热学时间步thermal timestep。热膨胀系数直接决定同样温升下颗粒能胀多大岩石类材料一般在 1e-5 到 1e-6 每摄氏度量级。热阻控制热量在颗粒间传递的快慢数值越大温度梯度越明显。热学时间步必须小于临界值否则温度场会震荡发散。我一般先用model solve thermal ratio检查热学步是否稳定再开始正式加载。2.2 平行胶结模型为什么是热损伤的“裁判”热损伤的本质是胶结键在热应力下断裂。PFC2D5.0 里常用的平行胶结模型parallel bond在颗粒接触点附近形成一个有限尺寸的胶结块同时传递力和力矩。它有两个关键强度参数法向强度pb_ten和切向强度pb_coh以及刚度pb_kn、pb_ks。温度升高时颗粒膨胀导致胶结块受拉或受剪一旦某个方向的应力超过对应强度胶结键就断了。这里有个容易翻车的地方平行胶结模型的强度是正态分布的不是固定值。你设pb_ten 10e6只是均值实际每个键的强度在均值附近随机波动。热损伤案例里如果强度离散太大裂纹会随机出现在低强度键上掩盖热应力集中的真实位置。我一般把强度标准差控制在均值的 10% 以内先跑纯力学加载验证模型再加温度场。2.3 空心圆颗粒生成孔隙率和粒径分布怎么定空心圆在 PFC2D5.0 里通过ball generate命令的radius范围控制但“空心”不是指颗粒内部有空腔而是指颗粒生成时按圆环区域分布中心留空。常见做法是用range限定一个圆环区域只在这个区域内生成颗粒。这样做的好处是初始模型就带有可控的孔隙率热传导路径更接近真实岩石。粒径分布一般用distribute关键字指定均匀分布或高斯分布。热损伤案例里粒径太均匀会导致热应力分布过于理想化裂纹路径不自然粒径差异太大又会让接触数波动剧烈热学求解容易不收敛。我一般用 1:1.5 到 1:2 的粒径比最小粒径根据模型尺寸定保证整个模型有 2000 到 5000 个颗粒接触数在 8000 到 15000 之间。空心圆区域的半径取模型总半径的 0.6 到 0.8 倍留出外围颗粒作为约束边界。3. 从零跑通一个热力耦合岩石热损伤案例代码逐段拆解3.1 模型初始化与颗粒生成; 热力耦合岩石热损伤案例 - 初始化 model new model large-strain on model thermal on ; 定义模型域 domain extent -0.05 0.05 -0.05 0.05 ; 设置随机种子保证可复现 model random 10001 ; 生成空心圆区域颗粒 ; 外半径 0.04内半径 0.024粒径 0.0008-0.0016 ball generate radius 0.0008 0.0016 ... box -0.04 0.04 -0.04 0.04 ... range annulus 0.024 0.04 ; 设置颗粒密度和热学参数 ball attribute density 2600 ... damp 0.7 ... thermal-expansion 8e-6 ... thermal-conductivity 2.5 ... specific-heat 900这段代码先打开热学模块model thermal on这是热力耦合的前提。model random固定随机种子保证每次生成的颗粒位置一致否则你调参数时根本分不清是参数变了还是随机分布变了。ball generate的range annulus就是空心圆区域内半径 0.024、外半径 0.04颗粒只在这个环带里生成。颗粒属性里thermal-expansion是热膨胀系数thermal-conductivity是导热系数specific-heat是比热容这三个是热学求解的核心参数。密度和阻尼是力学参数阻尼 0.7 是局部阻尼用于耗散动能热损伤这种准静态过程一般用 0.7 到 0.8。3.2 接触模型与胶结参数设置; 设置接触模型为平行胶结 contact cmat default model linearpbond ... method deformability emod 30e9 kratio 2.0 ... property pb_ten 12e6 pb_coh 15e6 ... pb_kn 1e11 pb_ks 0.5e11 ... fric 0.5 ; 在颗粒间建立胶结 contact method bond gap 0.0002 contact method property pb_ten 12e6 pb_coh 15e6 ... pb_kn 1e11 pb_ks 0.5e11 ; 初始化接触 model cycle 1contact cmat default定义默认接触材料linearpbond是线性平行胶结模型。emod 30e9是颗粒弹性模量对应硬岩的 30 GPa 量级。kratio 2.0是法向刚度与切向刚度比岩石类材料一般在 1.5 到 3 之间。pb_ten和pb_coh是胶结键的抗拉和抗剪强度这里设 12 MPa 和 15 MPa属于中等强度岩石。fric 0.5是摩擦系数胶结断裂后颗粒间靠摩擦接触。contact method bond gap 0.0002表示颗粒间距小于 0.2 mm 时才建立胶结这个值要略大于最小粒径否则很多接触点不会生成胶结。3.3 热学边界与温度加载; 设置热学边界条件 ; 内圈固定温度 20 度外圈固定温度 200 度 ball attribute temperature 20 range annulus 0.024 0.026 ball attribute temperature 200 range annulus 0.038 0.04 ; 固定热学边界颗粒 ball fix thermal range annulus 0.024 0.026 ball fix thermal range annulus 0.038 0.04 ; 设置热学时间步 model thermal timestep fix 0.05 ; 记录温度监测点 history thermal temperature ... ball position 0.032 0.0这段是热损伤案例的核心。内圈颗粒温度固定 20 度外圈固定 200 度中间颗粒自由传热。ball fix thermal把边界颗粒的温度锁死不让它们参与热学更新否则边界温度会漂移。model thermal timestep fix 0.05是手动指定热学时间步这个值需要根据颗粒热扩散系数和粒径估算一般取粒径平方 / (热扩散系数 * 10)量级。如果设太大温度场会震荡设太小计算时间成倍增加。history thermal temperature在半径 0.032 处放一个监测点记录温度随时间变化用来判断温度场是否达到稳态。3.4 力学求解与热损伤演化监测; 设置力学求解参数 model mechanical timestep scale model solve ratio-average 1e-5 ; 记录裂纹数量 history mechanical cracks ; 分阶段加载先热传导再力学平衡 model solve thermal time 100.0 model solve mechanical ratio-average 1e-5 ; 输出裂纹分布 plot export bitmap filename thermal_damage.pngmodel mechanical timestep scale让力学时间步自动缩放保证求解稳定。model solve ratio-average 1e-5是力学收敛标准平均不平衡力与平均接触力之比小于 1e-5 时认为平衡。history mechanical cracks记录裂纹数量随时间的变化这是判断热损伤程度最直接的指标。先model solve thermal time 100.0让温度场传导 100 秒再model solve mechanical让力学场平衡。实际案例里这两个求解是交替进行的PFC2D5.0 会自动在热学步和力学步之间切换但手动分阶段更容易观察每个阶段的变化。4. 热力耦合案例里最容易翻车的五个地方4.1 温度场震荡发散现象、原因与解决现象温度监测曲线上下跳动不收敛甚至出现负温度。原因热学时间步超过临界值或者热阻参数设置不合理。热学求解的稳定性条件是时间步小于R * C量级其中 R 是热阻C 是热容。解决先用model thermal timestep fix设一个很小的值比如 0.001跑几步看温度是否单调变化然后逐步放大到刚好不震荡的值。另外检查thermal-conductivity是否设得太大导热系数越大热阻越小临界时间步也越小。4.2 胶结键全部断裂现象、原因与解决现象加载初期裂纹数量就飙升模型直接散架。原因胶结强度设得太低或者颗粒膨胀系数太大温度稍微升高就产生超过强度的热应力。解决先做纯力学单轴压缩模拟确认胶结强度能承受预期应力水平。热膨胀系数从 1e-6 开始试逐步增加到 8e-6。另外检查contact method bond gap是否太大导致部分接触点没有胶结应力集中到少数胶结键上。4.3 空心圆区域颗粒飞散现象、原因与解决现象空心圆内边界颗粒在热加载过程中向外飞散。原因内圈边界颗粒只固定了温度没有固定力学位移热膨胀时它们自由移动。解决如果内圈代表孔洞颗粒飞散是正常的如果内圈代表约束边界需要用ball fix velocity固定内圈颗粒的力学自由度。我一般在内圈加一层固定颗粒只让中间区域自由膨胀。4.4 温度场达到稳态太慢现象、原因与解决现象跑了很久温度监测点还是没到目标值。原因热学时间步太小或者导热系数太低。解决先估算热扩散时间L^2 / alpha其中 L 是模型特征尺寸alpha 是热扩散系数。如果模型半径 0.04 m热扩散系数 1e-6 m²/s扩散时间约 1600 秒。热学时间步 0.05 秒意味着需要 32000 步计算量很大。可以适当放大时间步或者用model thermal timestep scale让软件自动缩放。4.5 裂纹分布与预期不符现象、原因与解决现象裂纹随机分布没有出现预期的热应力集中区。原因颗粒粒径分布太宽或者胶结强度离散太大。解决把粒径比控制在 1:1.5 以内胶结强度标准差降到均值的 5% 到 10%。另外检查model random种子是否固定不同种子生成的颗粒排列不同裂纹路径也会不同。做参数研究时固定种子做统计研究时至少跑 5 个不同种子取平均。5. 让热损伤案例从“能跑”到“可信”三个进阶技巧5.1 用温度梯度驱动代替固定温度边界固定温度边界虽然简单但温度梯度是突变的和真实岩石加热过程有差距。更接近实际的做法是用热源颗粒让热量从内圈向外圈传导。具体做法是在内圈放一圈高温颗粒不固定温度只给初始温度然后让热学求解自然传导。这样温度梯度是连续的热应力分布也更真实。代价是达到稳态的时间更长需要更大的热学时间步或更多的计算步数。5.2 用history和plot做热损伤演化曲线光看最终裂纹图不够热损伤是一个过程。我一般同时记录四个量监测点温度、裂纹总数、胶结断裂率、模型总动能。温度曲线判断热学是否稳定裂纹曲线判断损伤阶段断裂率判断损伤程度动能曲线判断力学是否平衡。把这四个量画在同一张图上能清楚看到温度传到某个位置时裂纹开始萌生温度梯度稳定后裂纹停止扩展。这个曲线比单张裂纹图有说服力得多。5.3 参数敏感性分析的批量跑法热力耦合案例的参数很多热膨胀系数、导热系数、胶结强度、粒径分布、孔隙率。一个一个手动改太慢。我一般用 PFC2D5.0 的model save和model restore配合外部脚本批量跑。先存一个基准模型然后写一个循环脚本每次model restore后改一个参数跑完保存结果。关键是要在脚本里自动提取裂纹数量和温度监测值输出到文本文件。这样跑 20 组参数可能只需要半天手动改参数可能两天都跑不完。提示批量跑之前先确认单次运行时间如果单次超过 30 分钟先优化模型规模减少颗粒数或放大时间步否则批量跑会变成通宵任务。我自己的习惯是每次改参数前先存一个base_model.pfc所有批量脚本都从这个基准恢复。有一次忘了存基准改了一轮参数后想回到初始状态结果只能重新生成颗粒浪费了一个下午。希望帮到你。本文还有配套的精品资源点击获取