简介基于PFC2D5.0离散元平台围绕热力耦合作用下岩石热损伤案例展开以空心圆盘岩样轴心受热为典型场景讲解材料属性定义、模型初始化、热载荷施加与热应力分析等关键环节适合岩石力学、采矿或地热工程方向的研究生与工程师参考。压缩包共10个文件包含txt文本说明、html技术页与doc文档分别覆盖背景介绍、案例实现思路、结果分析与技术博客式深入探讨整体仅23KB轻量便携。目前已有128人学习使用。资源重点呈现颗粒流离散元在热力耦合建模中的具体操作逻辑以及热损伤萌生、扩展的数值实现方法可帮助读者快速理解PFC2D热-力耦合模拟流程并迁移至不同温度边界或加载条件下的岩石损伤问题。 搞岩石热损伤数值模拟的同行应该对PFC2D5.0不陌生。这个版本最大的价值在于把颗粒流离散元和热力学计算放在同一个框架里跑能让我们直接从细观层面看温度怎么让岩石内部“裂开”。今天我就把自己做空心圆试件热力耦合模拟的一套流程和关键代码片段拿出来拆一拆重点讲清楚热参数怎么设、温度怎么加载、裂纹怎么监测哪些地方容易算飞掉都是我实际跑过的经验不是那种只有框架没有细节的东西。先说清楚这个模型的用途。空心圆试件在岩石力学里对应的是地下巷道、钻孔、或者地热井壁附近的围岩。外边界模拟远处原岩内孔壁模拟开挖临空面。当高温流体或者地热梯度作用在孔壁上岩石内部会产生温度梯度进而诱发热应力配合原有的地应力就可能在孔壁附近生成裂纹、剥落甚至整体失稳。PFC2D5.0里做这个本质上就是三件事生成一个带孔洞的颗粒集合体给它赋上热力学参数然后让温度场和应力场在颗粒尺度上互相作用。1. 为什么用离散元做热损伤连续方法算不到的地方1.1 岩石热损伤的细观机制温度对岩石的破坏不是简单的膨胀收缩。热传导进入岩石后不同矿物颗粒的膨胀系数不同颗粒之间会出现不协调变形即使对于均质岩石温度梯度本身也会在局部造成拉应力。当这些细观的拉应力或者剪应力超过颗粒间粘结强度微裂纹就开始萌生。随着温度持续作用和应力重分布微裂纹会沿着薄弱面扩展、合并最终形成宏观裂缝。传统有限元处理这种问题需要预设裂缝路径或者用损伤模型去“模糊”地描述开裂区域。但裂纹的方向、分叉、交汇这些在工程里恰恰非常重要。PFC的做法就直白得多——它把岩石当成一个一个的颗粒颗粒之间的粘结就是“断点”应力超限粘结直接断开裂纹就生成了不需要任何预设路径。1.2 离散元在热力耦合里的天然优势PFC2D5.0处理热问题的逻辑很有意思。它不是像有限元那样直接解温度场偏微分方程得到每个节点的温度而是把热传导也颗粒化了热流通过颗粒接触面传递每个颗粒储热颗粒温度变化后通过线膨胀系数换算成半径变化半径变了颗粒之间的重叠量就变了接触力跟着变这就实现了“热”到“力”的耦合。这个机制好处很明显。第一不需要额外划分网格热传导网络和力学接触网络是同一个东西第二裂纹生成是自发行为热应力导致粘结断开模型自动响应不会出现有限元里网格畸变、不收敛的问题。缺点也有微观参数标定麻烦计算量大这些后面细说。2. 模型搭建前的关键准备从几何到接触模型2.1 空心圆试件的几何与颗粒生成空心圆模型最好从简单的几何参数开始。我常用的尺寸是外径100 mm、内径20 mm这个比例能保证孔壁破坏有足够的围压作用范围。颗粒半径不宜过大否则模拟不出细观的裂纹扩展也不宜过小否则颗粒量太大、算不动。我一般取半径0.5~1.0 mm的均匀分布按照孔隙率0.10~0.15来生成算下来颗粒数大概在3000~4000个二维模型这个规模跑起来比较合适。如果你想让模型更贴近某类真实岩体可以要求生成后对颗粒按幂律或高斯分布重新赋值这样可以制造非均匀性。不过作为基础案例均匀级配更容易聚焦在热力耦合机制本身。2.2 接触模型选型为什么用平行粘结颗粒合模完成后下一步是赋予接触模型。模拟岩石我首选平行粘结模型Linear Parallel Bond。这个模型包含两个部分一部分是线弹性接触负责传递力和不产生粘结的颗粒间相互作用另一部分是平行粘结相当于在接触面上涂了一层“胶水”既能传力也能传力矩等效一个面积上的应力分布当法向或者切向应力超过设定强度这一层“胶水”就断了。为什么不用简单的contact-bond模型因为接触粘结只能传力不能传力矩模拟出的岩石破坏模式偏向脆性拉裂剪胀和压剪破坏的表现很差。平行粘结模型允许破坏后摩擦力继续存在这样裂纹两侧在压紧状态下还能承受剪应力更接近真实岩石的残余强度特性。2.3 热力学参数的物理对应关系在PFC2D5.0里热学属性分布在两个对象上颗粒属性包括密度、比热容、线膨胀系数接触属性包括热传导系数。需要特别提醒的是这里的“接触热传导系数”不是材料导热系数它表示两个颗粒单位温差下通过该接触的热流速率单位是W/°C需要和颗粒半径、材料导热系数之间做一个换算。给一个我常用的花岗岩参数组合作为起点参数数值备注颗粒密度 kg/m³2650花岗岩典型值颗粒半径 mm0.5~1.0均匀分布接触法向刚度 GPa8~12平面应变条件下标定接触切向刚度 GPa3~5约为法向的1/3平行粘结刚度 GPa8~12与接触刚度接近平行粘结法向强度 MPa35~50对应单轴抗压强度约100 MPa平行粘结切向强度 MPa25~35约为法向的0.7线膨胀系数 1/°C8e-6花岗岩典型值接触热传导系数 W/°C20~30由宏观导热系数换算比热容 J/(kg·°C)900花岗岩典型值这个表格里的强度参数不是一个固定值需要根据目标岩石的宏观力学参数反标定。我的经验是先定刚度再通过双轴压缩虚拟试验标定粘结强度最后测试抗拉强度。别指望一次性到位这部分工作往往要占整个项目三分之一的时间。3. 热力耦合代码骨架与核心片段解析3.1 先理顺热-力耦合的计算逻辑在展示代码之前先把顺序理顺。很多人上来就调用solve发现温度还没传导进去应力已经平衡了结果一团糟。正确的做法是先建立力学平衡的初始状态再开启热力学模块设置温度边界进行热力耦合循环计算。流程上是力学平衡 → 赋初始温度 → 激活热模块 → 加载温度边界 → 耦合循环。每一步都要确保上一步状态稳定了再进行下一步。3.2 空心圆试样生成与初始平衡代码PFC2D5.0的代码分命令流和FISH函数。我习惯混合使用主体流程用命令流后处理用FISH。下面这个片段生成空心圆试样并施加初始地应力平衡; ; 模型基本设置 ; new domain extent -0.08,0.08 set gravity 0.0 0.0 ; 外圆墙 wall id 1 arc center 0,0 radius 0.05 ... begin-angle 0 end-angle 360 ; ; 在圆环区域内生成颗粒 ; 用两个圆包围区域外圈边界 wall内圈边界 wall ; wall id 2 arc center 0,0 radius 0.01 ... begin-angle 0 end-angle 360 ball distribute porosity 0.12 radius 0.5e-3,1.0e-3 ... box -0.05,0.05 -0.05,0.05 ... tries 5000 ... range union [cone(0,0,0.01, 0,0,0.05)] ; 删除落入内孔区域的颗粒 ball delete range cylinder end1 0,0,0 end2 0,0,1 radius 0.01 ; ; 赋予颗粒属性 ; ball attribute density 2650 ball attribute kn 8e9 ks 3.5e9 ball property linear ... kn 8e9 ks 3.5e9 fric 0.6 ball property pb_kn 8e9 pb_ks 35e9 ... pb_ten 38e6 pb_coh 28e6 pb_fa 35 ; 接触模型设为平行粘结 contact cmat default model linearpbond ... ball kn 8e9 ball ks 3.5e9 ball fric 0.6 ... ball pb_kn 8e9 ball pb_ks 3.5e9 ... ball pb_ten 38e6 ball pb_coh 28e6 ball pb_fa 35 ; ; 初始力学平衡 ; cycle 2000 calm 100 solve arithmetic 0 1e6 10000这里有两个细节值得注意。第一ball distribute的range写法用union把外圈内的区域定义为生成范围然后再删除内孔颗粒比用严格圆环区域直接生成要稳不容易出现颗粒悬挂在边界上的情况。第二calm命令在初始平衡阶段特别有用它可以清除颗粒的初始动能让试样快速达到力平衡不然颗粒在重力域里会因为初始重叠乱蹦。3.3 温度场加载与热力耦合循环代码初始平衡稳定后进入热力耦合阶段。首先需要激活热配置设置材料热学参数然后指定孔壁温度或热流边界。PFC2D5.0里激活热力耦合的命令是model configure heat这个一定要放在初始力学平衡之后否则颗粒生成过程中的力学调整也会叠加热效应结果很难分开看。; ; 激活热学模块并设置热力学参数 ; model configure heat ; 颗粒热参数比热容、线膨胀系数、初始温度 ball attribute specific-heat 900.0 ball attribute expansion 8e-6 ball attribute temp 20.0 ; 接触热传导参数 contact property heat-conductance 25.0 ; 热时间步需要小于力学时间步的稳定条件 set thermal timestep 1e-5 set timestep 1e-7 ; ; 温度边界孔壁加热到200°C ; 使用FISH函数对半径小于内孔半径某阈值的颗粒设置温度 ; fish define set_heat local p ball.head loop while p # null local bx ball.pos.x(p) local by ball.pos.y(p) local br math.sqrt(bx*bx by*by) if br 0.011 then ball.temp(p) 200.0 endif p ball.next(p) endloop end set_heat ; ; 热力耦合求解分阶段循环 ; fish define thermal_step local n 0 while n 1000 command solve time 0.1 set_heat endcommand n n 1 endloop end thermal_step这段代码里set_heat函数是核心中的核心。它的作用是保证温度边界持续作用——每次力学计算结束后重新把孔壁附近颗粒温度强制设定为200°C否则热扩散进行后边界温度会被周围颗粒拉低相当于一个逐渐失效的边界条件这在长时间模拟里是一个非常隐蔽的坑。热时间步和力学时间步的匹配是另一个需要重点关注的细节。PFC2D5.0支持独立的thermal timestep和mechanical timestep。如果耦合循环里每一步都让热和力同步推进力学时间步要远小于热时间步否则温度变化引起的膨胀在单步里太剧烈会导致局部颗粒重叠量突变接触力震荡发散。我实际跑下来把热步取1e-5、力步取1e-7比较稳妥。如果发现温度扩散过快导致应力波震荡就把热步再缩小如果计算太慢可以把力学时间步适当放大但要密切监控不平衡力曲线不能为了速度牺牲稳定。4. 热损伤演化监测与结果分析4.1 裂纹自动监测与AE计数平行粘结模型的好处是每个接触的粘结断裂事件都会被记录为一条裂纹。PFC2D5.0内置了裂纹监测模块可以通过FISH函数统计裂纹总数、裂纹位置、断裂类型拉伸破坏或者剪切破坏。我习惯把裂纹数据实时导出等效成声发射监测。每个断裂事件的时刻就是一次AE事件累计裂纹数和时间的关系就是岩石热损伤的演化曲线。在热力耦合模拟里这个曲线能直观反映损伤进程初期平坦热能正在积累应力未超过强度没有裂纹快速上升局部应力超限裂纹集中萌生趋于平稳应力重新分布后达到新的平衡裂纹停止发展这个曲线形态对应实际试验里的声发射平静期、活跃期、衰减期非常像。4.2 损伤变量如何定义才算合理损伤变量的定义直接影响结果解释。学术论文里常见有几种方式基于裂纹数量D N / N_maxN是当前裂纹数N_max是最终裂纹数基于裂纹面积或长度把每个裂纹映射到连续介质单元用裂纹密度定义损伤基于弹性模量退化加载前后卸载刚度对比D 1 - E/E_0在PFC里最常用的是前两种。裂纹数量法实现最简单但忽略了裂纹位置聚集带来的局部损伤不均匀性裂纹密度法需要把圆环区域划分格子统计每个格子内的裂纹数能反映孔壁附近损伤集中的空间分布。我的建议是宏观上画裂纹数量-温度曲线空间上用裂纹密度云图。两者结合才能既看清损伤程度又看清损伤模式。4.3 从裂纹图到热损伤模式的判读热损伤的破坏模式在PFC里非常直观。空心圆孔壁周围如果以切向拉应力为主裂纹多为沿着孔壁方向扩展的径向裂纹如果围压较高、热应力以切向压应力为主破坏模式可能是片帮剥落裂纹平行于孔壁并逐渐向深部迁移。从应力角度解释更容易理解。孔壁是自由边界切向应力集中效应明显。加热孔壁后近壁区域膨胀受约束产生切向压应力远离孔壁区域温度没上来约束较弱。这种梯度应力场加上地应力叠加很容易在孔壁附近形成“受拉-受压”的交替区域裂纹模式也因此多样化。模拟结束后用PFC的可视化模块分别显示裂纹彩色云图、接触力分布图和颗粒位移矢量场。裂纹分布看模式接触力链看传力路径位移场看宏观变形趋势。三个图叠在一起才能把热损伤的机理讲清楚。5. 常见问题排查实录与参数标定心得5.1 热传导步长与计算发散前面提过热时间步和力学时间步不匹配的问题。如果遇到温度场还没扩散应力波已经开始震荡甚至颗粒飞散优先排查两步长设置。热时间步建议先取一个小值测试热传导稳定性比如1e-6跑几步看看孔壁温度是否按预期向远处扩散如果不发散再逐步放大。另一个常见问题是加热温度给得过高。颗粒在极小范围内半径突变相邻接触力瞬间超额导致局部力链崩溃。这种情况不是参数错了而是加载方式太激进。解决手段是分级升温比如每50°C作为一个加载阶段每阶段中间加力学平衡计算。5.2 微观参数标定的顺序问题颗粒流模拟最大的坑是以为微观参数等于宏观参数。实际上一套微观刚度、强度参数至少要通过三个虚拟标定试验才能锁定单轴压缩确定弹性模量和单轴抗压强度、直接拉伸或巴西劈裂确定抗拉强度、三轴压缩确定内摩擦角和残余强度。热学参数标定相对简单线性膨胀系数直接用矿物的实测值导热系数可以通过单颗粒链传热虚拟试验反推。注意一点接触热导率與颗粒半径有关同样的接触热导率颗粒越粗等效宏观导热系数越大。所以颗粒尺寸变了接触热导率要重新换算不能照搬。5.3 边界条件的隐藏坑热量损失与绝热控制热力耦合模拟里还有一个容易被忽略的问题模型外边界是绝热的还是有热流的默认情况下如果没有特别设置颗粒和墙体之间可能没有热传导热量会堆积在模型内部这相当于绝热边界。但现实中围岩深部是恒温的外边界应该设置成恒定温度或者热汇。我通常的做法是在模型外圈设置一圈温度恒定的边界颗粒模拟远场恒温条件。同时记录整个模型的能量平衡包括颗粒动能、应变能、摩擦耗能和热量流入量确保热量守恒合理如果热量在模型内异常累积边界条件一定有问题。另外热膨胀导致的颗粒重叠会在接触处积累非常大的法向力这种力在卸载后不会消失除非设置一个摩擦滑移或粘结断裂的释放通道。在热力耦合模拟中这种不可逆变形是真实的物理过程不需要额外处理但解释结果时要把这部分区分清楚有些“损伤”其实是塑性滑移不是开裂。我在实际项目中遇到的多数问题归根结底都是热步长、边界加载方式和颗粒尺寸三个变量没有协调好。如果你打算跑这类模拟我的建议是先把一套小尺寸模型跑通、跑稳再放大到目标尺寸不要一上来就全尺寸。还有一个小技巧每个阶段计算完成后导出颗粒坐标和裂纹文件这样即使后面算崩了也能从上一步继续而不是从头再来。空心圆热力耦合这套东西短期内不太可能被普通的连续方法完全替代因为裂纹自发萌生扩展这件事颗粒离散元做的事情确实更贴近物理本质。参数标定确实烦但一旦把流程跑顺了你会发现它对理解岩石在温度作用下的细观破坏机制帮助特别大。希望这份代码解析能帮你省掉一部分我当年踩过的弯路。本文还有配套的精品资源点击获取