Abaqus热力耦合断裂仿真：UMAT/VUMAT子程序开发与工程实践

简介：Abaqus热力耦合断裂代码开发[可运行源码]是一份面向Abaqus二次开发学习者的可运行源码包，聚焦UMAT与UEL子程序协同实现相场-温度场耦合的断裂模拟，适合材料科学、结构工程领域研究者及具备一定Abaqus基础的中高级用户。压缩包共3个文件，包含源码工程（inscode）、项目说明页（html）及gitignore配置，整体仅5KB，轻量便携。目前已有118人学习，适合快速理解核心实现。代码包围绕相场变量phi建立力学场与热学场联系，覆盖弹性矩阵相场退化、热传导方程残差计算、子程序间数据交换等关键模块，并附调试中节点校验、相场阈值判断等典型问题的解决方案。通过运行源码可生成温度-相场云图，直观呈现裂纹扩展与温度场的动态交互，为后续二次开发提供可直接修改和扩展的基础框架。 近几年做高温结构失效分析的项目多了起来，手头这个热力耦合断裂的活，我折腾了差不多两个月才把整套Abaqus代码跑通。市面上能找到的断裂仿真案例，大多是纯力学加载下的裂纹扩展，一旦牵扯温度场就麻烦不少——材料属性随温度漂移、热应变和塑性应变耦合、高温下损伤演化规律完全不是常温那套。这段时间我把整个过程踩过的坑、试对的路子、代码组织的思路整理出来，给准备入坑热力耦合断裂仿真的人做个参考。

1. 为什么热力耦合断裂一定要自己写代码：内置能力的边界

先说个扎心的事实：Abaqus自带的断裂相关模型，在纯力学场景下确实够用，但只要温度一掺和进来，就有很多对不上的地方。

拿内置的Damage for Ductile Metals来说，它支持Johnson-Cook损伤准则，可以设置温度相关的失效应变，看起来好像能用于热力耦合。但问题在于，它把温度对材料的影响简化成了几组表格插值，损伤演化的过程本身并没有和温度历史耦合——也就是说，你没办法描述“材料先经历了一个温度循环，损伤累积速率因此改变”这种真实物理过程。再比如Brittle Cracking和Concrete Damaged Plasticity，前者基本不考虑温度，后者虽然能带上温度场，但也是以外挂材料参数的形式，想嵌入一个温度相关的断裂韧性演化方程，几乎不可能。

XFEM呢？Abaqus的XFEM在纯力学裂纹扩展里表现不错，但热力耦合场景下，裂纹面两侧的温度边界条件、接触传热、裂纹张开后的热阻变化，内置的交互定义很难处理干净。我做热冲击仿真时试过XFEM，裂纹一旦起裂，温度场在裂纹面附近就开始振荡，增量步缩得没法看。

再说说单元删除这个最常用的“伪断裂”做法。内置的单元删除（Element Deletion）配合韧性损伤，可以模拟裂纹扩展的宏观效果，但它的损伤演化律是固定的那几种形式。工程上经常需要自己定义损伤驱动力——比如同时考虑应力三轴度、温度历史和应变率，内置模型就没法直接做了。

所以结论很明显：要在Abaqus里做真正的热力耦合断裂分析，自己写UMAT或VUMAT是绕不开的。所谓“真正的”，指的是：

• 本构方程里同时出现温度变量和损伤变量，两者互相影响
• 热应变的增量处理和力学应变的增量处理在同一个积分点完成，不存在“先算热再算力”的割裂
• 损伤演化驱动力可以由用户完全自定义，而不是只能用内置的那几种公式
• 如果分析类型是完全热力耦合，还要保证每次迭代里温度自由度和位移自由度的残差同时收敛

这也是市面上很多所谓“热力耦合断裂源码”卖得贵的原因，它不是一个简单UMAT模板能搞定的，需要把这些耦合关系全部理顺。我的建议是：能自己写就自己写，抄来的代码如果不理解它的耦合逻辑，一改参数就全线崩盘。

2. 断裂表征与热力耦合的核心：你得先想清楚的三个物理问题

写代码之前，有三个物理问题必须想清楚。这三个问题直接决定了你的代码架构，后面改起来才不会伤筋动骨。

2.1 裂纹用什么方式表征

连续损伤力学（CDM）加单元删除是工程上最稳妥的方案。核心思路是引入一个损伤变量D，D从0增长到1，达到1时单元刚度趋于零，然后触发单元删除，裂纹就以“撕开”的方式向前扩展。

这里有个关键细节：损伤演化最好用能量型（Displacement-based）而不是纯应变型（Strain-based）。原因很简单，纯应变型损伤严重依赖单元尺寸，网格密一倍的单元会比疏的单元先“断”，完全不符合物理。能量型损伤引入一个特征长度h（通常取单元面积的平方根），把应力-应变软化关系换算成应力-位移关系，断裂能Gf被固定下来，网格敏感性就能压到工程可接受的范围。Abaqus里这个特征长度会自动算好通过CELENT传给子程序，自己写代码时直接用就行。

2.2 热力耦合分析怎么选

这里有个概念很容易搞混：顺序耦合和完全耦合。

顺序耦合（Sequential Coupling）的流程是先跑一个纯热分析，把每个节点每个时刻的温度历史存下来，再以体荷载或预定义场的方式加载到力学分析中算断裂。这个方案速度快、稳定，但适用条件是温度场不受变形和损伤的影响。焊后冷却残余应力、热处理过程中的开裂这类问题，温度场基本由传热决定，顺序耦合足够用。

完全耦合（Fully Coupled）则是在每个增量步里同时求解温度自由度和位移自由度。碰到这些情况就必须上完全耦合：高速变形中的塑性功生热、摩擦生热、裂纹张开后改变了传热路径从而反过来影响温度分布。热力耦合断裂的很多场景——比如热成形模具表面的热疲劳裂纹，裂纹张开后高温气体窜进去加热裂纹尖端——本质上就是完全耦合问题。

从代码角度说，完全耦合意味着UMAT里不仅要处理力学本构，还要通过*COUPLED TEMPERATURE-DISPLACEMENT分析步把温度自由度和位移自由度绑在一起迭代，数值上比纯力学麻烦得多。

2.3 损伤准则和材料失效域怎么划

这个要看你模拟的材料。

金属材料在高温高应变率下，Johnson-Cook损伤模型是主力（式子里同时包含应力三轴度、应变率和温度项，正好覆盖热力耦合需求）。它的积分形式是：

[ D=\sum \frac{\Delta \bar{\varepsilon}_p}{\bar{\varepsilon}_f(\sigma^*, \dot{\bar{\varepsilon}}_p, T)} ]

当累积D达到1时单元失效。这个式子最关键的地方在于分母上的失效应变(\bar{\varepsilon}_f)是应力三轴度、应变率、温度的函数，意味着不同温度下材料的“抗断裂能力”不同，正好是热力耦合断裂想要的。

准脆性材料比如陶瓷、混凝土、高温涂层，用CDPM或加一个温度相关的指数软化损伤演化会更合适。这类材料没有明显的塑性段，起裂就是脆断，损伤演化方程往往是应力线性/指数衰减型。

我建议在设计代码时，把损伤准则做成一个可插拔的子模块，而不是写死在主程序里。这样换材料只需换一个函数，工程上大量省时间。

3. 代码架构与运行细节：从UMAT到可运行源码

源码能不能“跑起来”，取决于你对Abaqus子程序的接口机制理解得有多清楚。下面把代码骨架、状态变量布局、编译连接这几个关键点一次讲透。

3.1 整体代码骨架

我用的是“主程序+本构模块+损伤模块+热应变模块”的四层结构。主程序负责和Abaqus的接口通信，顺序是：

1. 读取材料常数 PROPS 2. 根据温变增量 DTEMP 计算热应变增量 DSTRAN_T 3. 总应变增量 DSTRAN 扣除热应变增量，得到力学应变增量 4. 弹性预测应力，计算屈服函数 5. 若屈服，径向返回映射做塑性修正，更新等效塑性应变 6. 调用损伤模块，更新损伤变量 D 7. 按损伤比例折减应力，更新状态变量 STATEV 8. 组装雅可比矩阵 DDSDDE

这里有个非常容易出错的点：Abaqus传给UMAT的DSTRAN是总应变增量，里面包含了热应变，必须自己把它扣除后再做塑性计算。如果忘了扣，结果就是温度升高时材料“自己变形自己硬化”，应力场完全乱套。正确做法是严格按照热膨胀系数和温度增量算出热应变，从总应变里减掉，剩下的才是驱动力学响应的部分。对于VUMAT，DSTRAN同样处理，只不过接口名字略有区别。

3.2 状态变量布局

状态变量（SDV）是整个子程序传递内部变量的通道，布局设计不好后面后处理会一脸懵。我用的分配方案是：

这样设计的好处是，后处理时只需要看SDV(2)就能判断裂纹扩展路径，SDV(4)可以直观看到哪些单元进入了待删除状态。注意在inp文件的*DEPVAR里要声明对应数量的状态变量，并且指定删单元时要SAVE到.sta文件里，方便排查删除时刻。

3.3 隐式与显式：谁更适合热力耦合断裂

这一条我纠结了很久，最后结论是分场景。

Abaqus/Standard + UMAT 的问题是断裂软化段的负切线刚度很容易导致Newton迭代不收敛。尤其是完全热力耦合分析，温度场迭代和位移场迭代相互激荡，增量步会被削到微秒级甚至直接崩掉。如果你只做准静态热力耦合断裂（比如稳态温度场下的裂纹扩展），可以用UMAT，但建议加上黏性正则化：在损伤更新方程里加一个阻尼项，把软化段的突变拉平。阻尼系数取特征时间步长的(10^{-5})到(10^{-3})这个量级，太小没效果，太大结果会偏离物理。

如果你要模拟的是快速裂纹扩展、冲击载荷、或者塑性功生热显著的过程，直接上Abaqus/Explicit + VUMAT。显式方法不需要迭代收敛，损伤软化导致的单元删除处理起来非常从容，而且VUMAT不需要提供DDSDDE，省掉了一大块麻烦。我最终交付的源码包就是VUMAT为主版本，配合显式完全热力耦合分析步，通用性和鲁棒性都好得多。

3.4 编译连接和作业提交

Abaqus子程序的编译环境匹配是个隐形大坑。我用的是Abaqus 2021 + Intel oneAPI Fortran编译器（Classic版），Windows下打开Abaqus Command后直接执行：

abaqus make library=vtc_umat.f

然后提交作业时带上子程序：

abaqus job=thermal_fracture user=vtc_umat.f cpus=4

注意cpus数量不要超过许可核数，不然会报cpus exceeds available的错。提交前建议先单独编译一次库文件，确认没有语法错误再提交正式分析，否则日誌刷屏找错很痛苦。

3.5 单元选择

热力耦合单元有专门的命名规则：CPE4T（平面应变四节点）、C3D8T（三维八节点）、C3D8RT（三维八节点减缩积分）等等。我的实测经验是：如果涉及单元删除，尽量用完全积分单元CPE4T或C3D8T，减缩积分配合单元删除在某些变形模式下容易出现沙漏，顺着裂纹路径冒出一排锯齿状变形，看起来很影响判断。

4. 实测最容易翻车的三个场景及排查链路

代码写出来之后才是噩梦的开始。下面三个问题是我实测中最常翻车的，每个都能让计算直接白跑几小时。

4.1 单元删除后的温度场突变

这是热力耦合断裂特有的大坑。单元删除在力学上是“刚度消失”，但在热学上，删掉的单元如果还参与传热计算，矛盾就来了——被删除单元的热导率如果仍然按原值参与计算，热量就会穿过已经开裂的裂纹面继续传导，物理上完全错误。

但如果你直接把被删单元的导热系数设为零，又会带来另一个问题：热流路径瞬间中断，裂纹尖端附近温度梯度暴涨，温度场出现尖锐的局部热点，然后导致网格畸变进一步加剧。

我的处理办法是做一个“热虚拟单元”过渡：单元删除后的若干个增量步内，不直接把导热系数清零，而是按损伤变量D线性折减到原值的10%左右，让热流有个衰减过程。这样既保证了裂纹面热阻增大这个物理事实，又不会让温度场剧烈跳变。

4.2 损伤软化段的收敛崩溃

UMAT跑完全热力耦合时，单元进入软化段后Newton迭代特别容易发散。表面现象是日志里出现很多个等号连续缩步，最后报“TIME INCREMENT REQUIRED IS LESS THAN THE MINIMUM SPECIFIED”。

根因是DDSDDE矩阵在损伤软化段出现了负特征值，整体刚度矩阵不再正定。我的排查链路是：

1. 先用纯力学分析（恒定温度场）跑一遍，把收敛问题和温度场解耦，确认是不是本构软化导致的
2. 在UMAT里给损伤演化增加黏性正则化项，公式改为[ \dot{D}=\frac{1}{\eta}(D_{eq}-D) ]，η取(10^{-4})倍的特征时间步
3. 如果还崩，把损伤对DDSDDE的贡献项先关闭（只更新SDV，不折减应力），确认问题出在雅可比矩阵装配上，而不是物理模型本身
4. 最后才考虑是不是温度场和位移场的迭代耦合太强，需要放宽温度收敛容差

这套流程我屡试不爽，基本能定位90%的收敛问题。

4.3 网格依赖和“撕扯状”裂纹

如果你发现裂纹路径是沿着网格边界的锯齿状，而不是一条符合物理的平滑曲线，大概率是损伤演化选型错了。

我刚跑通代码时，裂纹永远顺着网格线走，试了很多组参数都改不动。后来查了文献才明白，问题出在用应变型损伤而不是能量型损伤。用能量型损伤，把断裂能Gf和特征长度CELENT结合起来算等效位移，裂纹路径会对网格方向的敏感度大幅降低，路径终于从“齿轮状”变成了“自然状”。

还有个小技巧：网格要尽量在预计的裂纹路径上加密，同时把单元长宽比控制在1比1到1比1.5之间。单元太扁会造成损伤局部化方向被网格长边带偏，最后的裂纹走向跟实验差得很远。

5. 验证案例：单边缺口拉伸试件的完整标定过程

理论说得再多，不如拿一个案例把流程走一遍。这里用一个高温合金单边缺口拉伸（SENT）模型做验证，目标是复现试验件在650°C下拉伸断裂的全过程。

5.1 几何与网格

试件长100mm，宽20mm，一侧开2mm深的缺口，宽度方向单侧。我用C3D8T单元，缺口附近网格尺寸0.5mm，其余区域过渡到2mm。单元总数约1.2万，这个规模在单机上跑完全没压力。模型两侧约束，底部固定，顶部给位移载荷，拉伸速率设为0.5mm/min（为了和准静态试验对应，用显式分析时注意质量缩放要控制好，动能与内能之比始终低于5%）。

5.2 材料参数

以典型热作模具钢为例，温度相关参数如下：

Johnson-Cook损伤参数D1到D5是根据多组不同温度、不同应力三轴度的圆棒拉伸试验标定出来的。标定过程要提醒一句：不要直接抄文献里的D值，不同冶炼批次的材料D参数差异很大，有条件就自己做试验。实在没条件，至少也要做一组650°C下的拉伸试验，反推D5（温度项）的量级，否则高温段的失效点会和试验差出一大截。

5.3 模拟与试验对比

我跑了三组分析：纯力学常温拉伸、650°C恒定温度场拉伸、650°C完全热力耦合拉伸（考虑塑性功生热和裂尖绝热温升）。

结果很有代表性。常温组裂尖应力三轴度大约在0.65左右，裂纹基本沿缺口正前方扩展，路径较直；650°C恒定温度场下，由于材料软化，峰值载荷相比常温下降了约35%；完全热力耦合下，因为裂尖附近塑性功生热导致局部温度比周围高30-50°C，失效位置略有偏移，裂纹起始角比纯温度场组偏了约7度。

载荷-位移曲线和试验对比，峰值载荷误差在8%以内，断裂位移误差约12%。考虑到损伤参数本身的不确定性，这个精度对于工程预判已经完全够用。

5.4 源码包怎么组织

交付的“可运行源码”我是这样组织的，方便复现：

src/ vtc_vumat.f # 显式VUMAT主程序 vtc_umat.f # 隐式UMAT版本 material_constants.inc # 材料常数公共块 jc_damage.f # Johnson-Cook损伤模块 thermal_strain.f # 热应变计算模块 inp/ sent_explicit.inp # SENT试件显式热力耦合输入文件 sent_standard.inp # SENT试件隐式版本输入文件 run/ run_windows.bat # Windows批处理 run_linux.sh # Linux脚本

运行脚本里写好了环境变量、编译命令和作业提交命令，解压后按顺序执行两个脚本就能复现案例结果。里面还带了一个后处理Python脚本，自动提取载荷-位移曲线和裂纹扩展路径，省得手动提数据。

一些实操心得

如果只让我送一条经验，就是调试热力耦合断裂代码时，一定要分步解耦。我先关掉损伤（D恒等于0），跑一遍，确认热力耦合力学结果和Abaqus内置的经典塑性模型吻合；然后打开损伤但固定温度场，验证断裂路径对不对；最后才把温度和损伤全部打开。这三个层次每过一个再往前推，出了问题能立刻锁定是本构的问题、断裂模块的问题、还是热力耦合交互的问题。很多同行一上来就跑完全热力耦合断裂，崩了之后无从下手，其实都是没做这个分步验证。

至于这个源码往后还能怎么扩展，我目前想的是把GTN空洞损伤模型也做进去，替代Johnson-Cook去处理高应力三轴度下的韧性断裂，同时在VUMAT的接口上增加位错密度相关的物理本构，让高温蠕变和疲劳裂纹扩展也能塞进同一套框架里。热力耦合断裂这个方向的门槛确实不低，但一旦把代码框架理顺了，后面都是在往框架里填新物理模型的事，工作量会小很多。

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