搞过有限元界面仿真的人应该都绕不开 Cohesive 单元。Abaqus 自带的牵引-分离Traction-Separation模型做常规粘接失效分析确实方便但如果你想把自定义软化曲线、温度依赖、率相关效应塞进去内置模型很快就会变成一堵墙。这篇博客就是用一个小例子把一个带 UMAT 的内聚力本构从原理到代码、从建模到调试全过程拆开。适合刚接触 UMAT、想理解 Cohesive 单元背后力学逻辑的工程师和学生也适合那些想摆脱“内置模型黑盒”的人。我去年做胶接接头强度预测时就一度卡在内置 Cohesive 模型的损伤演化太理想化上。后来花了几天时间把 UMAT 子程序的写作、编译和调参流程整理清楚回头看才发现整个链条里真正的难点不在 Fortran 语法而在你能不能把“牵引力、分离位移、断裂能、损伤状态”这套逻辑串起来。下面我就按自己的实战路径把这个主题讲透。1. 为什么需要自定义 Cohesive 本构 UMAT1.1 内置牵引-分离模型够用但不永远是够用的Abaqus 内置的 Cohesive 行为本质上是一个“牵引力-分离位移”的本构关系。施加在外载荷下胶层或界面的上下两个面分开一点点距离内部会产生对应的粘聚力抵抗分离。模型里常见的双线性、指数型、多项式型等已经覆盖了大部分线性软化问题的需求。所以很多人会问内置模型都在了还自己写 UMAT 干什么我的回答是内置模型照顾的是“标准化”场景。比如损伤起始准则通常写成最大应力准则或二次应力准则损伤演化用能量或位移控制。但实际工程里界面本构的损伤起始和演化往往具备非对称拉压行为、混合模式耦合更复杂、温度场影响明显、甚至还要考虑固化度变化。这时候如果只在 Abaqus 界面里点鼠标你只能接受软件预设的数学形式无法进行深度定制。UMAT 的意义就是把“材料行为”这个模块真正交到你手里让你自己定义应力更新和切线刚度矩阵。1.2 什么场景下必须写 UMAT结合我自己接触过的项目下面几类需求基本得上 UMAT需要把实验测得的力-位移曲线直接作为软化段输入而不是强行用双线性或指数型曲线拟合。界面损伤受到温度和湿度联合影响损伤阈值和断裂能都是场变量函数。需要加入速率相关项模拟不同加载速率下的粘接失效行为。要研究非均质界面比如胶层中间夹杂颗粒或局部改性区域。模型要集成到优化流程里频繁修改本构参数用 UMAT 可以完全绕开 GUI 的参数入口限制。写 UMAT 当然有门槛它要求你具备 Fortran 代码阅读能力理解 Abaqus/Standard 的增量求解流程还要掌握切线刚度矩阵的推导。但一旦你把第一个简单模型跑通后面换准则、换演化方程都是在一个框架内增加分支而已收益会非常明显。2. 内聚力本构模型的核心原理2.1 Cohesive 单元的“应变”到底意味着什么先把概念理清。Cohesive 单元和普通实体单元最大的区别在于它模拟的不是一个体积材料的变形而是两个面之间的分离行为。单元上下表面的相对位移通常用法向分离量 δn 和两个切向分离量 δs、δt 来描述。Abaqus 在 UMAT 中传递给子程序的“应变”并不是真正的几何应变而是经过本构厚度归一化的名义应变。也就是说εn δn / T0εs δs / T0εt δt / T0。这个 T0 是你设定的“本构厚度”可以等于单元的几何厚度也可以单独指定。使用默认设置时Abaqus 会取 Cohesive 单元的初始几何厚度作为本构厚度。如果单元初始厚度很小即使同样的分离位移名义应变也会变得很大。这一点是很多新手在给 UMAT 配参数时分不清单位、导致结果异常的根本原因。在我自己的 UMAT 里我倾向于先把 T0 通过 PROPS 数组传入然后直接将 STRAN 乘以 T0还原成真实的分离位移 δ再进行损伤判断和应力计算。这样思维上更贴近牵引-分离定律也便于直接和复合材料层间断裂实验数据对照。2.2 双线性内聚力模型的数学表达最经典的内聚力模型是双线性模型。它的曲线分两段第一段是线弹性阶段界面刚度保持恒定牵引力随分离位移线性增加第二段是软化阶段牵引力线性下降直到完全失效。曲线下的面积就是单位面积上的断裂能 Gc。如果用符号描述法向初始刚度为 Kn法向峰值强度为 tn0完全失效对应的法向分离量为 δn_f。那么第一段应力为 σn Kn·δn达到峰值后损伤变量 D 从 0 增加到 1应力变为 σn (1 - D)·Kn·δn。损伤演化率与断裂能相关常见形式是D [δmax·(δf - δ0)] / [δ·(δf - δ0_等效)]实际上双线性模型在混合模式下通常使用有效位移 δm 来统一损伤变量。有效位移的定义是δm sqrt(⟨δn⟩² δs² δt²)这里的麦考利括号表示法向压缩时不计入损伤驱动量。也就是说纯压缩不会引起混合模式损伤这和工程直觉一致压紧界面不会惩罚粘接强度。2.3 初始刚度、强度与断裂能的标定写内聚力 UMAT 时核心参数就三类界面刚度、峰值强度和断裂能。界面刚度一般取 10⁵~10⁶ N/mm³ 量级。取值不宜过大过大容易造成数值病态迭代矩阵条件数恶化取值过小则界面在加载初期就出现“虚假分离”整体结构刚度偏低。峰值强度对应损伤起始点最好通过单搭接剪切或 DCB 试验获得。断裂能则决定软化段的长度和曲线下降斜率由实验力-位移曲线面积计算。混合模式下的断裂能很多学者用幂指数准则或 BK 准则描述。我这里用一个简单例子只考虑法向和剪切两个方向所以按照二次应力准则判断损伤起始采用 BK 准则计算混合模式断裂能即可代码中需要你传入 G_I_c 和 G_II_c以及 BK 指数 η一般 η 在 1~2 之间。3. 从零走一遍 UMAT 实例3.1 UMAT 接口与基本框架Abaqus 的 UMAT 接口是固定的无论你做什么材料子程序入口都一样。关键参数是STRAN当前应变数组DSTRAN应变增量数组STRESS应力数组进入子程序时是旧应力出来时要更新为新增量步后的应力DDSDDE一致切线刚度矩阵即 ∂Δσ / ∂ΔεSTATEV状态变量数组用于保存历史信息比如最大有效位移和损伤变量Fortran 代码骨架大致如下SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, 1 RPL, DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, 2 TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED, CMNAME, 3 NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, COORDS, 4 DROT, PNEWDT, CELENT, DFGRD0, DFGRD1) INCLUDE ABA_PARAM.INC DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS,NTENS) DIMENSION STRAN(NTENS), DSTRAN(NTENS) DIMENSION PROPS(NPROPS), COORDS(3) CHARACTER*80 CMNAME C 这里开始写本构逻辑 ... RETURN END对于 Cohesive 单元NTENS 通常等于 3分别对应于一个法向和两个剪切方向。当然如果你用的是平面或轴对称 Cohesive 单元分量顺序要参考 Abaqus 单元手册不能凭借实体单元的经验想当然。3.2 双线性本构的 Fortran 实现我在这里给出一段简化但可运行逻辑的双线性内聚力 UMAT 核心代码片段基于混合模式有效位移判断损伤。PROPS 参数表如下参数含义PROPS(1)本构厚度 T0PROPS(2)法向初始刚度 KnPROPS(3)剪切初始刚度 KsPROPS(4)法向峰值强度 tn0PROPS(5)剪切峰值强度 ts0PROPS(6)法向断裂能 G_I_cPROPS(7)剪切断裂能 G_II_cPROPS(8)BK 准则混合指数 eta代码核心部分C T0 PROPS(1) C KN PROPS(2) C KS PROPS(3) T0 PROPS(1) KN PROPS(2) KS PROPS(3) C Abaqus传入的是名义应变换算成分离位移 DN STRAN(1) * T0 DS STRAN(2) * T0 DT STRAN(3) * T0 C 计算有效位移法向压缩不贡献损伤 DN_POS MAX(DN, 0.0D0) DM SQRT(DN_POS*DN_POS DS*DS DT*DT) C 读取历史变量 DM_MAX MAX(STATEV(1), DM) STATEV(1) DM_MAX D STATEV(2) IF (DM .GT. 0.0D0) THEN IF (D .EQ. 0.0D0) THEN C 判断损伤起始二次应力准则 IF ((DN_POS/TN0)**2 (DS/TS0)**2 (DT/TS0)**2 .GE. 1.0D0) 1 THEN D 1.0D0 END IF END IF END IF C 计算软化段参数 C DM0为损伤起始位移DMF为完全失效位移都需要根据断裂能计算 DM0 0.0D0 DMF 0.0D0 C 这里省略具体公式实际代码中需要调用BK准则计算混合模式断裂能 C 如果已经进入损伤更新D IF (D .GT. 0.0D0) THEN DM_EFF DM_MAX GR GE / (DMF - DM0) IF (DM_EFF .LT. DMF) THEN D DM0*(DM_EFF - DM0) / (DM_EFF*(DMF - DM0)) ELSE D 1.0D0 END IF END IF C 更新应力 STRESS(1) (1.0D0 - D) * KN * DN STRESS(2) (1.0D0 - D) * KS * DS STRESS(3) (1.0D0 - D) * KS * DT C 法向受压时损伤不折减法向刚度 IF (DN .LT. 0.0D0) THEN STRESS(1) KN * DN END IF C 更新切线刚度DDSDDE此处代码略需要求导得到一致切线矩阵 ... RETURN END注意上面这段代码为了展示核心逻辑省略了初始损伤判断时需要确定 DM0、DMF 的细节。真实可用的代码中DM0 由当前混合比下的界面强度反算DMF 由当前混合比下的断裂能反算。每个增量步都要先求当前应力比的混合模式断裂能再更新 D。最终 DDSDDE 矩阵也需要根据 D 对 δ 的导数推导。建议先在纸上推导清楚再写到 Fortran 里不要直接盲写。3.3 单搭接剪切试件的建模实操我这个实例选的模型是单搭接剪切试件这是粘接结构中最常见的验证算例。两块 25mm×25mm 的铝板厚度 2mm搭接长度 25mm中间胶层厚度 0.2mm。一侧自由端固定另一侧施加位移控制载荷模拟剪切破坏。几何建模不用复杂操作直接在 Abaqus/CAE 里拉伸三个 Part上板、下板、中间胶层。胶层的 0.2mm 厚度可以单独建一个薄层 Part。材料方面铝板用弹性模量 70GPa、泊松比 0.33 的线弹性材料胶层部分不赋予传统材料而是给 Cohesive 单元截面。网格划分注意三点胶层厚度方向只划分一层使用 COH3D8 单元上下铝板用 C3D8R 实体单元避免沙漏问题胶层与上下板交界面的节点要网格对齐。最稳妥的做法是在同一个 Part 里把板子和胶层分别切分出区域再用 Tie 约束绑定或者在划分网格时共享节点。在定义界面单元时截面类型要选择 Cohesive Section并把响应类型设为 Traction Separation。由于我们使用 UMAT材料行为不再用 Abaqus 自带的内聚行为而是在材料模型里选择 User Material。分别设置材料常数的数目和数值Kn、Ks 取 10⁶ N/mm³tn0 取 30MPats0 取 25MPaG_I_c 取 0.5 N/mmG_II_c 取 1.0 N/mmη 取 1.5。分析步用通用静态分析Static, General。打开几何非线性求解器默认 Newton-Raphson。加载位移总量建议 2mm初始增量步 0.005最小增量步 1e-6最大增量步 0.1。这样能保证软化阶段不至于跳步太多。3.4 结果怎么看、怎么判计算完成后先看胶层单元的状态变量 SDV即 STATEV 中保存的损伤变量 D分布。损伤会从搭接端部开始萌生逐步向中间扩展。再看载荷-位移曲线典型形态是初始线性段、达到峰值后逐渐软化、最终载荷降到接近零。如果你发现曲线峰值对应的位移远低于预期大概率是界面刚度 Kn、Ks 取值偏小导致结构在加载初期就产生了软变形。如果峰值强度对应位移正常但下降段很陡、甚至跳变可能是断裂能设置偏低或增量步不够细。对比实验曲线时重点关注峰值强度、曲线下降斜率、以及失效位移三者之间的匹配关系不要只盯着某个节点应力看。4. 常见问题与排查技巧实录4.1 UMAT 编译不过去怎么办每一位写 UMAT 的人都会碰到编译报错最常见的三种情况。第一种是环境变量或 Fortran 编译器版本不匹配。Abaqus 对 Fortran 版本有严格规定有的版本只能搭配 Intel oneAPI 老版本。我建议先在命令行单独编译一个最简 UMAT 验证环境再回头编译正式子程序。第二种是 Fortran 语法细节。例如 ABA_PARAM.INC 必须放在变量声明之后、所有可执行语句之前数组维度不能用变量动态声明长行续行符必须在第 6 列使用数字或字符。Abaqus 的预编译逻辑和标准 Fortran 略有差异很多新手会栽在旧式固定格式代码的缩进上。第三种是子程序内部变量类型不匹配。UMAT 中的 PROPS 数组默认是单精度如果你的材料常数超过单精度可表示范围需要显式转换成双精度。更简单的做法是在子程序开头把 PROPS 值赋给局部双精度变量。我的调试手段是在子程序入口加几个 WRITE 语句把 STRAN、DSTRAN、TIME 等关键变量输出到外部文件然后用一个小模型跑第一个增量步检查有没有 NaN 或异常负刚度。4.2 软化段不收敛的典型解法内聚力模型最大的数值难点是损伤软化阶段切线刚度从正值变成负值容易导致整体切线矩阵非正定Newton 迭代不收敛。我自己试过位移加载比力加载收敛性好得多所以强烈建议模型里用位移控制。第二个有效手段是粘性正则化。Abaqus 内置 Cohesive 模型有 viscosity 参数UMAT 里需要自己实现。思路很简单在损伤变量更新时引入粘性使其按指数松弛逼近无粘性解也就是D_new (D_unregularized - D_old) / (1 μ/Δt) D_old其中 μ 是粘性系数通常取 1e-4~1e-3 量级。粘性系数不宜过大否则峰值强度和断裂能会被明显钝化结果失真。如果仍然不收敛把最小增量步再调小一个量级同时检查是不是有个别 Cohesive 单元过度畸变。Cohesive 单元允许的变形范围有限一旦完全失效单元刚度接近零Abaqus 可能报告负特征值。这时可以在失效单元上通过删除单元或设置失效后的最小刚度来规避但不能影响整体应力分布。4.3 网格敏感性与初始刚度选择Cohesive 模型本身对网格尺寸有敏感性原因是软化带宽受单元大小影响。如果你用内置模型Abaqus 会按单元特征长度正则化断裂能让最终耗散能量保持稳定。UMAT 的好处是你可以自己控制这个正则化过程但代价是实现复杂度增加。初始刚度选取不当是另一个高频问题。界面刚度太小整体模型刚度和实验对不上太大则DDSDDE矩阵接近病态收敛困难。我通常采用一个经验值Kn α·E_interface / T0其中 α 在 10~100 之间。E_interface 可取胶层材料的弹性模量T0 为胶层实际厚度。这样初始刚度和真实厚度关联起来可以有效避免几何层厚度变化带来的不一致。调试时多做一个对比很有帮助先用 Abaqus 内置双线性模型跑一遍同样的几何和边界条件再切换到 UMAT 版本。两者的载荷-位移曲线应该在合理容差内重合。如果差距很大那说明 UMAT 参数表或者单位换算有问题这时候优先检查 T0 的取值。最后再分享一个实用小技巧写 UMAT 时第一轮先不做损伤把应力更新写成纯线弹性验证 DDSDDE 和应力更新是否一致。然后加上损伤起始判断再逐步引入损伤演化。这样每一步出错都能定位到具体逻辑块而不是在完整代码里大海捞针。我自己就是通过这种“三步走”策略把最初到处是 bug 的内聚力 UMAT 调成了一个稳定可复用的版本。这套逻辑后续还能轻松扩展出温度依赖子程序、率相关子程序只要理解了 Cohesive 单元和牵引-分离本构的底层关系写其他自定义材料也不是难事。