先说个真实的场景去年我做胶接接头失效分析模型本身不复杂两块铝合金板中间一层胶但一到损伤软化阶段Abaqus就是死活不收敛负特征值一条接一条增量步从0.01一路缩到1e-7。折腾了一周最后发现不是模型的问题而是我对Cohesive单元和它背后的内聚力本构模型理解错了。后来我换成自己写的UMAT把牵引-分离关系在代码里一步步捋清楚问题一下就通了。这篇内容准备从Cohesive单元是什么、内聚力本构模型怎么翻译成UMAT代码、以及怎么用一个最简单实例把它跑通这三个层面整体讲一遍。如果你正在做复合材料分层、胶接强度、裂纹扩展这类仿真或者论文里需要自定义界面损伤模型这篇应该能帮你省下不少弯路。1. 先搞清楚Cohesive单元到底在算什么东西1.1 内聚力模型是为了补上“断裂过程区”这一课传统断裂力学拿应力强度因子K和能量释放率G去描述裂纹在大部分金属结构里够用但遇到胶层、复合材料层间、混凝土这类材料时问题就来了真实裂纹尖端并不是一个数学上的奇点而是存在一个尺寸虽小但真实存在的损伤区域材料在这个区域内逐渐软化、失去承载力然后才形成新表面。这个区域叫断裂过程区FPZ。内聚力模型Cohesive Zone Model的思路很直接不去纠结裂纹尖端的应力场有多复杂而是在可能开裂的界面上定义一对“牵引力”和“分离位移”的关系。你把它想象成撕双面胶刚开始拉的时候力很大胶层慢慢变形到某个峰值后胶层开始脱粘拉力逐渐下降最后完全断开。这个“力-位移”过程就是一条牵引-分离曲线。Cohesive单元在有限元里就是专门用来承载这条曲线的。所以当你看到Abaqus里的Cohesive单元不要把它当成普通实体单元看。它不是一个有体积的“材料块”而是一个厚度方向上只有一对节点、专用来描述界面分离行为的“界面单元”。胶层、复合材料层间、岩石裂缝都适合用它来模拟。1.2 Abaqus里的Cohesive单元几何形态和坐标系Abaqus里常用的是COH2D4也就是二维4节点cohesive单元以及三维的COH3D8、COH3D6。从外观上看Cohesive单元有上下两个面中间有厚度方向。这个厚度在几何上可以很薄甚至可以为零厚度计算时真正生效的是“本构厚度”constitutive thickness。本构厚度这个概念很容易被忽略但它特别重要。Abaqus里Cohesive单元的应变不是普通意义下的应变更严格地说叫“名义应变”它等于界面分离量除以本构厚度ε δ / T0如果本构厚度T0取1那应变数值上就等于分离量写UMAT时会方便很多。如果取真实几何厚度0.1 mm那应变就要放大10倍。后面讲UMAT时你会看到这个T0没处理对整个软化段都会错位。另一个容易踩坑的是单元局部坐标系。Abaqus里Cohesive单元的默认局部1方向是厚度方向也就是界面的法向2、3方向是剪切方向。这意味着UMAT里的STRAN(1)是法向变形STRAN(2)和STRAN(3)是两个切向变形。如果你的网格扫掠方向不对法向和切向互换损伤计算就全乱了。还有一点Cohesive单元的名词里常听到“traction”和“separation”。Traction是单位面积上的力量纲是N/mm²也就是MPaSeparation是分离位移量纲是mm。常规应力-应变关系里应力应变相乘得到能量密度而牵引-分离曲线下的面积直接就是断裂能量纲是N/mm。2. 双线性内聚力本构三条曲线里的门道2.1 弹性段、损伤起始和线性软化最经典的内聚力本构是双线性模型Abaqus内置的Cohesive Behavior默认就是这类。曲线分三段上升段、软化段、失效段。上升段代表界面还没损伤牵引力随分离量线性增加斜率就是界面刚度K。到达损伤起始点后曲线掉头向下直到完全失效。我画不出图但可以用参数把这条曲线说得清楚。弹性段写成t K × δ其中K的单位是N/mm³或者叫MPa/mm。这里有一个常见的问号K和材料弹性模量E是什么关系如果你把胶层看成厚度为T0的一串弹簧那K就等于E除以T0K E / T0所以K不是一个独立变化的量它取决于你的本构厚度假设。ABAQUS内置cohesive behavior里用户直接输K值UMAT里也建议直接输K省得在E和T0之间来回换算。损伤起始点由最大牵引力决定。对法向来说当牵引力达到t0比如20 MPa对应的损伤起始分离量δ0 t0 / K软化段的终点是δf此时界面彻底失效。双线性模型假设软化段是直线断裂能就是三角形面积Gc 0.5 × t0 × δf所以给了断裂能Gc可以反算出完全失效的分离量δf 2 × Gc / t0举个具体例子设K1e5 N/mm³t020 MPaGc0.5 N/mm那么δ02e-4 mmδf0.05 mm。这个δf通常比δ0大两个数量级左右很常见。你别小看这三个参数它们决定了整个界面行为。K太高会让损伤起始前的弹性阶段太“硬”计算中容易出现收敛振荡K太低会让整体刚度偏弱影响未损伤阶段的变形。t0决定界面在什么载荷下开始损伤Gc决定损伤扩展需要的能量本质上决定裂纹是否稳定扩展。2.2 损伤变量、不可逆性和混合模态损伤一旦开始材料刚度就不能恢复了。双线性模型用损伤变量D来描述刚度退化D [δf × (δmax - δ0)] / [δmax × (δf - δ0)]其中δmax是历史上曾经达到的最大分离量。之所以要用δmax而不是当前δ是为了处理卸载和再加载。当你把位移往回拉时界面没有完全恢复它会沿着刚度退化后的路径卸载也就是刚度变成(1-D)K。如果不保存δmax卸载时D可能变小就会出现“损伤愈合”这种物理上不可能的事情。这也是自己写UMAT时最容易犯的错误之一。在单模态下损伤起始判断很简单对比当前分离量δ和δ0或者对比当前牵引力t和t0。Abaqus内置模型里常用的还有最大名义应力准则max(tn/tn0, ts/ts0, tt/tt0) ≥ 1以及二次名义应力准则sqrt((tn/tn0)² (ts/ts0)² (tt/tt0)²) ≥ 1实际胶接和复合材料层间开裂很少是单纯的法向或者剪切往往是混合模态。混合模态下损伤演化得额外定义断裂能与模态比的关系比如BK准则Gc Gn (Gs - Gn) × [Gs/(GnGs)]^η这里的Gn是法向断裂能Gs是剪切断裂能η是材料参数。Abaqus内置模型支持这些自己写UMAT时如果只想复现双线性建议先做单模态跑通了再扩展混合模态。不要一开始就想着把所有模态写全那样调试起来根本分不清问题出在本构还是出在代码。3. UMAT子程序把本构模型翻译给Abaqus听3.1 UMAT在Standard分析里到底被调用来干什么Abaqus/Standard在每一个积分点、每一次迭代里都会调用UMAT。你需要做的就是根据传入的总应变STRAN、应变增量DSTRAN、历史状态变量STATEV、材料参数PROPS更新应力数组STRESS并给出Jacobian矩阵DDSDDE。Jacobian在隐式分析里是牛顿迭代的切线刚度它的物理意义是当前应力对应变的偏导数。对cohesive界面来说就是当前切线刚度d(t)/d(δ)。弹性段给K软化段给负斜率完全失效后给一个很小的值而不是0否则刚度矩阵奇异收敛直接崩。在这里有个关键点对Cohesive单元STRAN存储的是名义应变不是分离位移。如前所述应变乘以本构厚度T0才能得到真正的分离位移δ。如果你在Abaqus的Section里把本构厚度设成1那就可以直接把STRAN当成δ用代码里不用再做换算非常推荐。如果非要用真实几何厚度那所有应力和Jacobian都要跟着缩放调试时极容易出错。写UMAT前还需要明确状态变量怎么分配。我习惯这样安排STATEV(1)损伤变量D初始0STATEV(2)历史最大分离量δmax初始0STATEV(3)损伤起始标志位0表示未起始1表示已起始为什么要有标志位因为损伤的发展方向是单向的damage_flag可以辅助判断是否已经进入软化段。当你后续扩展疲劳载荷、加卸载循环场景时状态变量的设计会直接影响程序可维护性。3.2 一个可运行的双线性内聚力UMAT核心框架我不建议直接贴一个几百行的完整代码因为接口冗长还容易把人绕晕。先看本构计算核心逻辑再把它套进标准UMAT模板里思路会清晰很多。以单模态法向、线性软化为例写成伪代码1. 读取材料参数: PROPS(1) Kn 界面法向刚度 PROPS(2) Tn0 法向损伤起始牵引力 PROPS(3) Gn 法向断裂能 2. 读取历史变量: D STATEV(1) DELTAMAX STATEV(2) 3. 计算当前分离量: DELN STRAN(1) * T0 ! T0为本构厚度建议取1 4. 计算损伤起始分离量和失效分离量: DELTA0 Tn0 / Kn DELTAF 2.0 * Gn / Tn0 5. 更新历史最大分离量: DELTAMAX MAX(DELTAMAX, DELN) 6. 判断状态并更新应力: IF (DELN DELTA0) THEN ! 弹性段 TN Kn * DELN DDD Kn ELSEIF (DELN DELTAF) THEN ! 软化段, 直接线性软化 TN Tn0 * (DELTAF - DELN) / (DELTAF - DELTA0) DDD -Tn0 / (DELTAF - DELTA0) ELSE ! 完全失效 TN 0.0 DDD 1.0e-6 * Kn END IF 7. 根据当前分离量计算损伤变量用于状态输出和卸载刚度: IF (DELN DELTA0 .AND. DELN DELTAF) THEN D DELTAF * (DELTAMAX - DELTA0) / 1 (DELTAMAX * (DELTAF - DELTA0)) ELSE D 1.0 END IF 8. STRESS(1) TN 9. DDSDDE(1,1) DDD 10. STATEV(1) D STATEV(2) DELTAMAX这段逻辑很直白核心就是第五步和第六步的顺序先更新历史最大分离量再算应力。如果顺序反过来软化段应力会算错。配合一个简化的Fortran片段会更有体感C T0设为适配Abaqus的cohesive单元本构厚度取1 T0 1.0D0 KN PROPS(1) TN0 PROPS(2) GN PROPS(3) C 当前法向分离量 DELN STRAN(1) * T0 C 特征分离量 DELTA0 TN0 / KN DELTAF 2.0D0 * GN / TN0 C 更新历史最大分离量 STATEV(2) MAX(STATEV(2), DELN) IF (DELN .LT. DELTA0) THEN STRESS(1) KN * DELN DDSDDE(1,1) KN ELSEIF (DELN .LT. DELTAF) THEN STRESS(1) TN0 * (DELTAF - DELN) / 1 (DELTAF - DELTA0) DDSDDE(1,1) -TN0 / (DELTAF - DELTA0) ELSE STRESS(1) 0.0D0 DDSDDE(1,1) 1.0D-6 * KN END IF C 损伤变量用于后处理和卸载刚度 IF (DELN .GT. DELTA0 .AND. DELN .LT. DELTAF) THEN STATEV(1) DELTAF * (STATEV(2) - DELTA0) / 1 (STATEV(2) * (DELTAF - DELTA0)) ELSEIF (DELN .GE. DELTAF) THEN STATEV(1) 1.0D0 ELSE STATEV(1) 0.0D0 END IF这个片段只写了法向单模态。要扩展剪切模态把PROPS(4)、PROPS(5)、PROPS(6)给Ks、Ts0、Gs然后对STRESS(2)和STRESS(3)做同样的处理。混合模态要再引入有效分离量和模态比代码复杂度会上一个台阶但基本框架不变。还有一点必须提醒UMAT里给的DDSDDE是d(应力)/d(应变)在Abaqus的cohesive单元里它对应的量纲是界面刚度乘以本构厚度这里容易绕。如果直接在本构厚度T01的情况下界面刚度K的数值上就等于dσ/dε所以上面代码里弹性段DDSDDEKN软化段为负斜率逻辑是对的。一旦T0不等于1DDSDDE必须乘T0否则应力对应变的切线会差一个数量级。我见过太多人查了几天也没找到刚度差在哪最后就是这个问题。3.3 内置Cohesive Behavior和UMAT怎么选Abaqus自带Cohesive Behavior不需要写代码使用起来非常简单在材料属性里定义K、损伤起始准则、损伤演化准则就行适用于绝大多数标准双线性模型。那为什么还要写UMAT因为内置模型给用户的自由度有限。比如你要做指数软化、梯形软化或者损伤起始应力不再是固定值而跟应力三轴度、温度、应变率有关或者界面在循环载荷下刚度退化方式很特别内置模型就覆盖不了。UMAT的所有逻辑都掌握在你自己手里理论上你想怎么写都行。代价是你要自己保证本构的物理合理性以及收敛性。我的建议是能复现问题用内置模型先跑通流程确定参数和边界没问题再替换成UMAT。不要一上来就写UMAT否则模型网格出问题、接触没设好、分析步设置不对你会误以为是UMAT写错了排查线索全被带偏。4. 单搭接剪切实例从单单元校准到整模型跑通4.1 先用一个单元把UMAT校准到理论曲线上写完后第一个测试一定不要直接上大模型先用一个单元验证。在Abaqus里建一个COH2D4单元单元节点可以取为(0,0), (1,0), (1,1), (0,1)也就是边长为1 mm的方形。给它分配UMAT材料本构厚度设成1。边界条件底边两个节点固定顶部两个节点施加竖直向上位移让cohesive单元受拉。这个模型简单到连收敛问题都很难发生。跑完后取顶部节点总反力除以横截面积1 mm²得到牵引力。再取节点的竖直位移得到分离量δ。把这两个量画成曲线和理论双线性曲线对比。理论曲线是t从0到20 MPa然后线性下降到0对应的δ范围从0.0002 mm到0.05 mm。这里有个操作细节如果你直接在step里给顶部节点一个大位移比如0.06 mmAbaqus会一上来就进入大变形软化阶段。建议分成几个分析步或者用一个幅值曲线从0逐步加到0.06这样能看到完整曲线。我自己调试UMAT时一定会做这一步并且把数值输出的每一行都跟理论值手算一次。弹性段第一个数据点如果t不等于K×δ检查T0软化段如果斜率不对检查DDSDDE失效段如果还有残余应力那正常数值上必须保留小刚度。4.2 单搭接剪切模型的建模步骤单单元验证通过后就可以搭一个完整的单搭接剪切模型。几何尺寸建议参考经典胶接试样上下两块铝板长75 mm宽25 mm厚1.5 mm搭接长度25 mm。中间胶层几何厚度0.1 mm但cohesive单元的本构厚度仍然设成1原因前面说过。材料参数这样配铝板线弹性E70 GPa泊松比0.33。分析重点是界面失效板的塑性可以暂时不加避免结果解读复杂化。胶层cohesive单元Kn1e5 N/mm³Ks1e5 N/mm³tn020 MPats020 MPaGn0.5 N/mmGs1.0 N/mm。这个参数组合模拟的是一种中等强度结构胶。网格划分是关键。Cohesive单元不能用普通的自由网格生成必须用扫掠Sweep方式沿厚度方向扫出。为什么因为cohesive单元要求上下两个面的节点一一对应而且厚度方向就是局部1方向。如果你用Tet网格去剖分cohesive单元厚度方向乱掉计算出来的法向和切向完全是错的。具体操作把胶层厚度方向划分成一层cohesive单元单元类型选COH2D4板用平面应变单元CPE4R或者平面应力CPE4都可以。单元尺寸控制在0.5 mm左右搭接区网格加密保证损伤过程区里有至少几个单元。边界条件建议上板左端固定下板右端施加水平向右位移0.5 mm强制搭接区发生整体剪切。这个设置比垂直拉伸更容易激发剪切主导的界面损伤。分析步用Static, General初始增量设小一点比如1e-5 mm的位移增量更准确地说设置初始增量步1e-5最大增量步0.01最小增量步1e-8。开启非对称求解器选项因为cohesive软化后刚度矩阵不再对称非对称求解器能明显改善收敛性。打开后记得在后处理里勾选SDV输出不然STATEV是空的你都不知道损伤发展到哪里。4.3 UMAT结果和内置模型互相验证跑完之后先在Visualization里看变形和状态变量。损伤变量D的云图应该是从搭接区两端开始逐渐向中间扩展最终形成一条贯穿的对角线损伤带这是单搭接剪切最典型的破坏形态。然后提取下板右端参考点的支反力RF和位移U画载荷-位移曲线。再把材料换成内置的Cohesive Behavior参数完全一样重新跑一遍把两条曲线叠在一起。理想情况下它们应该高度重合。如果曲线有差异优先检查UMAT里是否把本构厚度T0当成1但cohesive behavior内部也默认T0各不相同。内置模型里K直接定义但默认本构厚度是几何厚度需单独指定。这个“单位不统一”是两条曲线不一致最常见的原因。我在这个实例里踩过一个大坑内置模型算出来的初始刚度比UMAT模型低很多。查到最后发现内置cohesive behavior必须搭配截面属性里的constitutive thickness来设置而我没有给它指定导致Abaqus默认用了几何厚度0.1 mm来换算名义应变。UMAT里我写了T01两者刚度自然差了10倍。这个坑说明一个道理做对比之前先确认两个模型对“本构厚度”的约定完全一致。5. 常见问题与排查经验这些坑我都替你踩过5.1 收敛不了负特征值、增量步一降再降Cohesive单元做隐式分析最大的敌人就是收敛问题。典型症状是消息区不停出现负特征值警告增量步从0.1缩到1e-6然后直接中止。负特征值并不代表模型物理上真的不稳定很多时候是某个cohesive单元刚度退化成零导致局部刚度矩阵奇异。我处理这种问题优先级是这样的在UMAT中对完全失效单元不返回0刚度而是返回一个很小的残余刚度比如初始刚度的1e-6倍。这能显著减少数值奇异。在材料属性里加粘性正则化内置模型有visco选项UMAT里可以自己给损伤变量D做粘性更新。Abaqus内置的粘性系数一般取0.001就能压住振荡。调整增量策略初始增量调小使用自动增量打开非对称求解器。检查是不是网格畸变或者单元方向错乱。cohesive单元严重畸变时结果不如直接重画。还要注意一点如果结构里同时存在接触接触收敛和cohesive收敛会叠加在一起问题更复杂。调试阶段尽量先把接触去掉用绑定或共用节点保证网格连续等cohesive部分稳定了再加接触。5.2 损伤云图诡异、应力振荡有一种很令人抓狂的现象D已经显示接近1但应力场还在乱跳。大概率问题出在单元局部坐标系的朝向。所有cohesive单元默认1方向是厚度方向但如果网格扫掠时单元反转某些单元的1方向和其他单元反了一个受拉的界面变成受压损伤变量永远不会增长。排查方法是使用Abaqus的COORD选项输出单元坐标系或者在Model-Edit Attributes里查看单元方向。后处理里可以画出cohesive单元的S1法向应力、S2切向应力来直观判断看云图是否连续过渡。如果发现同一层cohesive单元里应力符号正负交替八九不离十是方向问题。另一个导致应力振荡的原因是软化段刚度过陡。比如你给了很高的断裂能但δ0和δf差距过小软化斜率非常陡单元刚度在几步内从K跳到接近零隐式迭代自然不稳定。这种情况要么微调Gc和t0让曲线更缓要么给损伤演化加粘性正则化。内聚力模型的好处是你可以通过调整Gc来控制软化斜率但注意Gc必须来自于真实的实验测量不能为了收敛而无限制放大。5.3 其他几个容易掉进去的坑我之前顺手整理过一个速查表直接分享出来问题现象可能原因解决办法SDV输出总是0分析步输出请求里没勾选SDV在Field Output中勾选SDV初始刚度比预期小10倍本构厚度T0设置不一致统一按T01处理并检查截面厚度损伤扩展方向不对cohesive单元厚度方向扫掠反了检查单元局部坐标1方向完全失效后单元刚度为0导致负特征值软化后刚度归零残余刚度设为初始刚度的1e-6倍收敛一步步缩小到1e-8还是不收敛模型里cohesive层被压溃而非拉伸检查边界条件是否造成局部压缩载荷-位移曲线震荡损伤起始强度过高、软化过陡降低t0或增大Gc或加粘性正则化还有一个容易被忽略的点单位制。前面所有参数我都是在mm、N、MPa、N/mm这套单位制下写的。如果你换成m、kg、s制断裂能Gc的单位是J/m²界面刚度K单位是N/m³数值会差得非常远。不少人把论文里的Gc直接抄进Abaqus结果差了三个数量级都不知道。5.4 网格细度和断裂过程区的关系内聚力模型对网格尺寸有一定敏感性但比纯应力-应变软化本构好太多。一个经验法则是断裂过程区的长度大致可以用这个公式估算L ≈ E × Gc / (t0)²其中E是界面附近的材料模量。如果过程区尺寸是0.2 mm你的网格尺寸就要明显小于这个值否则损伤带只在一个单元里发展结果表现得像脆断看不出渐进损伤。实际操作中我通常让过程区内至少有3到5个单元。这个估算不需要非常精确用来判断网格尺度的量级足够。如果模型太大、网格太细算不动怎么办优先考虑对称模型和子模型。做单搭接剪切时如果结构和载荷对称可以取半模型或者四分之一模型。当然要注意边界条件是否允许对称。内聚力模型本身计算量不大但软化阶段增量步很小模型小了收益非常明显。我个人在实际操作中的体会是UMAT的编写并不是最难的部分难的是理解cohesive单元的约定和调试时的耐心。每次看到负特征值先别急着改本构把单元坐标系、本构厚度、单位制这几样最基础的东西过一遍往往能找到罪魁祸首。我现在遇到界面开裂问题一定会先写一个单模态UMAT跑单单元确认曲线吻合理论再上完整模型这个习惯帮我省掉了至少一半的调试时间。如果你只是复现双线性模型内置cohesive behavior确实够用但后续要加温度、湿度、疲劳退化这类因素还是早点上手UMAT更省事。最后再分享一个小技巧在UMAT里多输出几个状态变量比如把当前牵引力、历史最大分离量、损伤起始标志分别存到不同的SDV里调试时对照云图就能一眼定位到某个积分点算到哪一步了这比只盯着D一个变量清晰得多。
Abaqus中Cohesive单元与UMAT实现双线性内聚力本构模型详解
先说个真实的场景去年我做胶接接头失效分析模型本身不复杂两块铝合金板中间一层胶但一到损伤软化阶段Abaqus就是死活不收敛负特征值一条接一条增量步从0.01一路缩到1e-7。折腾了一周最后发现不是模型的问题而是我对Cohesive单元和它背后的内聚力本构模型理解错了。后来我换成自己写的UMAT把牵引-分离关系在代码里一步步捋清楚问题一下就通了。这篇内容准备从Cohesive单元是什么、内聚力本构模型怎么翻译成UMAT代码、以及怎么用一个最简单实例把它跑通这三个层面整体讲一遍。如果你正在做复合材料分层、胶接强度、裂纹扩展这类仿真或者论文里需要自定义界面损伤模型这篇应该能帮你省下不少弯路。1. 先搞清楚Cohesive单元到底在算什么东西1.1 内聚力模型是为了补上“断裂过程区”这一课传统断裂力学拿应力强度因子K和能量释放率G去描述裂纹在大部分金属结构里够用但遇到胶层、复合材料层间、混凝土这类材料时问题就来了真实裂纹尖端并不是一个数学上的奇点而是存在一个尺寸虽小但真实存在的损伤区域材料在这个区域内逐渐软化、失去承载力然后才形成新表面。这个区域叫断裂过程区FPZ。内聚力模型Cohesive Zone Model的思路很直接不去纠结裂纹尖端的应力场有多复杂而是在可能开裂的界面上定义一对“牵引力”和“分离位移”的关系。你把它想象成撕双面胶刚开始拉的时候力很大胶层慢慢变形到某个峰值后胶层开始脱粘拉力逐渐下降最后完全断开。这个“力-位移”过程就是一条牵引-分离曲线。Cohesive单元在有限元里就是专门用来承载这条曲线的。所以当你看到Abaqus里的Cohesive单元不要把它当成普通实体单元看。它不是一个有体积的“材料块”而是一个厚度方向上只有一对节点、专用来描述界面分离行为的“界面单元”。胶层、复合材料层间、岩石裂缝都适合用它来模拟。1.2 Abaqus里的Cohesive单元几何形态和坐标系Abaqus里常用的是COH2D4也就是二维4节点cohesive单元以及三维的COH3D8、COH3D6。从外观上看Cohesive单元有上下两个面中间有厚度方向。这个厚度在几何上可以很薄甚至可以为零厚度计算时真正生效的是“本构厚度”constitutive thickness。本构厚度这个概念很容易被忽略但它特别重要。Abaqus里Cohesive单元的应变不是普通意义下的应变更严格地说叫“名义应变”它等于界面分离量除以本构厚度ε δ / T0如果本构厚度T0取1那应变数值上就等于分离量写UMAT时会方便很多。如果取真实几何厚度0.1 mm那应变就要放大10倍。后面讲UMAT时你会看到这个T0没处理对整个软化段都会错位。另一个容易踩坑的是单元局部坐标系。Abaqus里Cohesive单元的默认局部1方向是厚度方向也就是界面的法向2、3方向是剪切方向。这意味着UMAT里的STRAN(1)是法向变形STRAN(2)和STRAN(3)是两个切向变形。如果你的网格扫掠方向不对法向和切向互换损伤计算就全乱了。还有一点Cohesive单元的名词里常听到“traction”和“separation”。Traction是单位面积上的力量纲是N/mm²也就是MPaSeparation是分离位移量纲是mm。常规应力-应变关系里应力应变相乘得到能量密度而牵引-分离曲线下的面积直接就是断裂能量纲是N/mm。2. 双线性内聚力本构三条曲线里的门道2.1 弹性段、损伤起始和线性软化最经典的内聚力本构是双线性模型Abaqus内置的Cohesive Behavior默认就是这类。曲线分三段上升段、软化段、失效段。上升段代表界面还没损伤牵引力随分离量线性增加斜率就是界面刚度K。到达损伤起始点后曲线掉头向下直到完全失效。我画不出图但可以用参数把这条曲线说得清楚。弹性段写成t K × δ其中K的单位是N/mm³或者叫MPa/mm。这里有一个常见的问号K和材料弹性模量E是什么关系如果你把胶层看成厚度为T0的一串弹簧那K就等于E除以T0K E / T0所以K不是一个独立变化的量它取决于你的本构厚度假设。ABAQUS内置cohesive behavior里用户直接输K值UMAT里也建议直接输K省得在E和T0之间来回换算。损伤起始点由最大牵引力决定。对法向来说当牵引力达到t0比如20 MPa对应的损伤起始分离量δ0 t0 / K软化段的终点是δf此时界面彻底失效。双线性模型假设软化段是直线断裂能就是三角形面积Gc 0.5 × t0 × δf所以给了断裂能Gc可以反算出完全失效的分离量δf 2 × Gc / t0举个具体例子设K1e5 N/mm³t020 MPaGc0.5 N/mm那么δ02e-4 mmδf0.05 mm。这个δf通常比δ0大两个数量级左右很常见。你别小看这三个参数它们决定了整个界面行为。K太高会让损伤起始前的弹性阶段太“硬”计算中容易出现收敛振荡K太低会让整体刚度偏弱影响未损伤阶段的变形。t0决定界面在什么载荷下开始损伤Gc决定损伤扩展需要的能量本质上决定裂纹是否稳定扩展。2.2 损伤变量、不可逆性和混合模态损伤一旦开始材料刚度就不能恢复了。双线性模型用损伤变量D来描述刚度退化D [δf × (δmax - δ0)] / [δmax × (δf - δ0)]其中δmax是历史上曾经达到的最大分离量。之所以要用δmax而不是当前δ是为了处理卸载和再加载。当你把位移往回拉时界面没有完全恢复它会沿着刚度退化后的路径卸载也就是刚度变成(1-D)K。如果不保存δmax卸载时D可能变小就会出现“损伤愈合”这种物理上不可能的事情。这也是自己写UMAT时最容易犯的错误之一。在单模态下损伤起始判断很简单对比当前分离量δ和δ0或者对比当前牵引力t和t0。Abaqus内置模型里常用的还有最大名义应力准则max(tn/tn0, ts/ts0, tt/tt0) ≥ 1以及二次名义应力准则sqrt((tn/tn0)² (ts/ts0)² (tt/tt0)²) ≥ 1实际胶接和复合材料层间开裂很少是单纯的法向或者剪切往往是混合模态。混合模态下损伤演化得额外定义断裂能与模态比的关系比如BK准则Gc Gn (Gs - Gn) × [Gs/(GnGs)]^η这里的Gn是法向断裂能Gs是剪切断裂能η是材料参数。Abaqus内置模型支持这些自己写UMAT时如果只想复现双线性建议先做单模态跑通了再扩展混合模态。不要一开始就想着把所有模态写全那样调试起来根本分不清问题出在本构还是出在代码。3. UMAT子程序把本构模型翻译给Abaqus听3.1 UMAT在Standard分析里到底被调用来干什么Abaqus/Standard在每一个积分点、每一次迭代里都会调用UMAT。你需要做的就是根据传入的总应变STRAN、应变增量DSTRAN、历史状态变量STATEV、材料参数PROPS更新应力数组STRESS并给出Jacobian矩阵DDSDDE。Jacobian在隐式分析里是牛顿迭代的切线刚度它的物理意义是当前应力对应变的偏导数。对cohesive界面来说就是当前切线刚度d(t)/d(δ)。弹性段给K软化段给负斜率完全失效后给一个很小的值而不是0否则刚度矩阵奇异收敛直接崩。在这里有个关键点对Cohesive单元STRAN存储的是名义应变不是分离位移。如前所述应变乘以本构厚度T0才能得到真正的分离位移δ。如果你在Abaqus的Section里把本构厚度设成1那就可以直接把STRAN当成δ用代码里不用再做换算非常推荐。如果非要用真实几何厚度那所有应力和Jacobian都要跟着缩放调试时极容易出错。写UMAT前还需要明确状态变量怎么分配。我习惯这样安排STATEV(1)损伤变量D初始0STATEV(2)历史最大分离量δmax初始0STATEV(3)损伤起始标志位0表示未起始1表示已起始为什么要有标志位因为损伤的发展方向是单向的damage_flag可以辅助判断是否已经进入软化段。当你后续扩展疲劳载荷、加卸载循环场景时状态变量的设计会直接影响程序可维护性。3.2 一个可运行的双线性内聚力UMAT核心框架我不建议直接贴一个几百行的完整代码因为接口冗长还容易把人绕晕。先看本构计算核心逻辑再把它套进标准UMAT模板里思路会清晰很多。以单模态法向、线性软化为例写成伪代码1. 读取材料参数: PROPS(1) Kn 界面法向刚度 PROPS(2) Tn0 法向损伤起始牵引力 PROPS(3) Gn 法向断裂能 2. 读取历史变量: D STATEV(1) DELTAMAX STATEV(2) 3. 计算当前分离量: DELN STRAN(1) * T0 ! T0为本构厚度建议取1 4. 计算损伤起始分离量和失效分离量: DELTA0 Tn0 / Kn DELTAF 2.0 * Gn / Tn0 5. 更新历史最大分离量: DELTAMAX MAX(DELTAMAX, DELN) 6. 判断状态并更新应力: IF (DELN DELTA0) THEN ! 弹性段 TN Kn * DELN DDD Kn ELSEIF (DELN DELTAF) THEN ! 软化段, 直接线性软化 TN Tn0 * (DELTAF - DELN) / (DELTAF - DELTA0) DDD -Tn0 / (DELTAF - DELTA0) ELSE ! 完全失效 TN 0.0 DDD 1.0e-6 * Kn END IF 7. 根据当前分离量计算损伤变量用于状态输出和卸载刚度: IF (DELN DELTA0 .AND. DELN DELTAF) THEN D DELTAF * (DELTAMAX - DELTA0) / 1 (DELTAMAX * (DELTAF - DELTA0)) ELSE D 1.0 END IF 8. STRESS(1) TN 9. DDSDDE(1,1) DDD 10. STATEV(1) D STATEV(2) DELTAMAX这段逻辑很直白核心就是第五步和第六步的顺序先更新历史最大分离量再算应力。如果顺序反过来软化段应力会算错。配合一个简化的Fortran片段会更有体感C T0设为适配Abaqus的cohesive单元本构厚度取1 T0 1.0D0 KN PROPS(1) TN0 PROPS(2) GN PROPS(3) C 当前法向分离量 DELN STRAN(1) * T0 C 特征分离量 DELTA0 TN0 / KN DELTAF 2.0D0 * GN / TN0 C 更新历史最大分离量 STATEV(2) MAX(STATEV(2), DELN) IF (DELN .LT. DELTA0) THEN STRESS(1) KN * DELN DDSDDE(1,1) KN ELSEIF (DELN .LT. DELTAF) THEN STRESS(1) TN0 * (DELTAF - DELN) / 1 (DELTAF - DELTA0) DDSDDE(1,1) -TN0 / (DELTAF - DELTA0) ELSE STRESS(1) 0.0D0 DDSDDE(1,1) 1.0D-6 * KN END IF C 损伤变量用于后处理和卸载刚度 IF (DELN .GT. DELTA0 .AND. DELN .LT. DELTAF) THEN STATEV(1) DELTAF * (STATEV(2) - DELTA0) / 1 (STATEV(2) * (DELTAF - DELTA0)) ELSEIF (DELN .GE. DELTAF) THEN STATEV(1) 1.0D0 ELSE STATEV(1) 0.0D0 END IF这个片段只写了法向单模态。要扩展剪切模态把PROPS(4)、PROPS(5)、PROPS(6)给Ks、Ts0、Gs然后对STRESS(2)和STRESS(3)做同样的处理。混合模态要再引入有效分离量和模态比代码复杂度会上一个台阶但基本框架不变。还有一点必须提醒UMAT里给的DDSDDE是d(应力)/d(应变)在Abaqus的cohesive单元里它对应的量纲是界面刚度乘以本构厚度这里容易绕。如果直接在本构厚度T01的情况下界面刚度K的数值上就等于dσ/dε所以上面代码里弹性段DDSDDEKN软化段为负斜率逻辑是对的。一旦T0不等于1DDSDDE必须乘T0否则应力对应变的切线会差一个数量级。我见过太多人查了几天也没找到刚度差在哪最后就是这个问题。3.3 内置Cohesive Behavior和UMAT怎么选Abaqus自带Cohesive Behavior不需要写代码使用起来非常简单在材料属性里定义K、损伤起始准则、损伤演化准则就行适用于绝大多数标准双线性模型。那为什么还要写UMAT因为内置模型给用户的自由度有限。比如你要做指数软化、梯形软化或者损伤起始应力不再是固定值而跟应力三轴度、温度、应变率有关或者界面在循环载荷下刚度退化方式很特别内置模型就覆盖不了。UMAT的所有逻辑都掌握在你自己手里理论上你想怎么写都行。代价是你要自己保证本构的物理合理性以及收敛性。我的建议是能复现问题用内置模型先跑通流程确定参数和边界没问题再替换成UMAT。不要一上来就写UMAT否则模型网格出问题、接触没设好、分析步设置不对你会误以为是UMAT写错了排查线索全被带偏。4. 单搭接剪切实例从单单元校准到整模型跑通4.1 先用一个单元把UMAT校准到理论曲线上写完后第一个测试一定不要直接上大模型先用一个单元验证。在Abaqus里建一个COH2D4单元单元节点可以取为(0,0), (1,0), (1,1), (0,1)也就是边长为1 mm的方形。给它分配UMAT材料本构厚度设成1。边界条件底边两个节点固定顶部两个节点施加竖直向上位移让cohesive单元受拉。这个模型简单到连收敛问题都很难发生。跑完后取顶部节点总反力除以横截面积1 mm²得到牵引力。再取节点的竖直位移得到分离量δ。把这两个量画成曲线和理论双线性曲线对比。理论曲线是t从0到20 MPa然后线性下降到0对应的δ范围从0.0002 mm到0.05 mm。这里有个操作细节如果你直接在step里给顶部节点一个大位移比如0.06 mmAbaqus会一上来就进入大变形软化阶段。建议分成几个分析步或者用一个幅值曲线从0逐步加到0.06这样能看到完整曲线。我自己调试UMAT时一定会做这一步并且把数值输出的每一行都跟理论值手算一次。弹性段第一个数据点如果t不等于K×δ检查T0软化段如果斜率不对检查DDSDDE失效段如果还有残余应力那正常数值上必须保留小刚度。4.2 单搭接剪切模型的建模步骤单单元验证通过后就可以搭一个完整的单搭接剪切模型。几何尺寸建议参考经典胶接试样上下两块铝板长75 mm宽25 mm厚1.5 mm搭接长度25 mm。中间胶层几何厚度0.1 mm但cohesive单元的本构厚度仍然设成1原因前面说过。材料参数这样配铝板线弹性E70 GPa泊松比0.33。分析重点是界面失效板的塑性可以暂时不加避免结果解读复杂化。胶层cohesive单元Kn1e5 N/mm³Ks1e5 N/mm³tn020 MPats020 MPaGn0.5 N/mmGs1.0 N/mm。这个参数组合模拟的是一种中等强度结构胶。网格划分是关键。Cohesive单元不能用普通的自由网格生成必须用扫掠Sweep方式沿厚度方向扫出。为什么因为cohesive单元要求上下两个面的节点一一对应而且厚度方向就是局部1方向。如果你用Tet网格去剖分cohesive单元厚度方向乱掉计算出来的法向和切向完全是错的。具体操作把胶层厚度方向划分成一层cohesive单元单元类型选COH2D4板用平面应变单元CPE4R或者平面应力CPE4都可以。单元尺寸控制在0.5 mm左右搭接区网格加密保证损伤过程区里有至少几个单元。边界条件建议上板左端固定下板右端施加水平向右位移0.5 mm强制搭接区发生整体剪切。这个设置比垂直拉伸更容易激发剪切主导的界面损伤。分析步用Static, General初始增量设小一点比如1e-5 mm的位移增量更准确地说设置初始增量步1e-5最大增量步0.01最小增量步1e-8。开启非对称求解器选项因为cohesive软化后刚度矩阵不再对称非对称求解器能明显改善收敛性。打开后记得在后处理里勾选SDV输出不然STATEV是空的你都不知道损伤发展到哪里。4.3 UMAT结果和内置模型互相验证跑完之后先在Visualization里看变形和状态变量。损伤变量D的云图应该是从搭接区两端开始逐渐向中间扩展最终形成一条贯穿的对角线损伤带这是单搭接剪切最典型的破坏形态。然后提取下板右端参考点的支反力RF和位移U画载荷-位移曲线。再把材料换成内置的Cohesive Behavior参数完全一样重新跑一遍把两条曲线叠在一起。理想情况下它们应该高度重合。如果曲线有差异优先检查UMAT里是否把本构厚度T0当成1但cohesive behavior内部也默认T0各不相同。内置模型里K直接定义但默认本构厚度是几何厚度需单独指定。这个“单位不统一”是两条曲线不一致最常见的原因。我在这个实例里踩过一个大坑内置模型算出来的初始刚度比UMAT模型低很多。查到最后发现内置cohesive behavior必须搭配截面属性里的constitutive thickness来设置而我没有给它指定导致Abaqus默认用了几何厚度0.1 mm来换算名义应变。UMAT里我写了T01两者刚度自然差了10倍。这个坑说明一个道理做对比之前先确认两个模型对“本构厚度”的约定完全一致。5. 常见问题与排查经验这些坑我都替你踩过5.1 收敛不了负特征值、增量步一降再降Cohesive单元做隐式分析最大的敌人就是收敛问题。典型症状是消息区不停出现负特征值警告增量步从0.1缩到1e-6然后直接中止。负特征值并不代表模型物理上真的不稳定很多时候是某个cohesive单元刚度退化成零导致局部刚度矩阵奇异。我处理这种问题优先级是这样的在UMAT中对完全失效单元不返回0刚度而是返回一个很小的残余刚度比如初始刚度的1e-6倍。这能显著减少数值奇异。在材料属性里加粘性正则化内置模型有visco选项UMAT里可以自己给损伤变量D做粘性更新。Abaqus内置的粘性系数一般取0.001就能压住振荡。调整增量策略初始增量调小使用自动增量打开非对称求解器。检查是不是网格畸变或者单元方向错乱。cohesive单元严重畸变时结果不如直接重画。还要注意一点如果结构里同时存在接触接触收敛和cohesive收敛会叠加在一起问题更复杂。调试阶段尽量先把接触去掉用绑定或共用节点保证网格连续等cohesive部分稳定了再加接触。5.2 损伤云图诡异、应力振荡有一种很令人抓狂的现象D已经显示接近1但应力场还在乱跳。大概率问题出在单元局部坐标系的朝向。所有cohesive单元默认1方向是厚度方向但如果网格扫掠时单元反转某些单元的1方向和其他单元反了一个受拉的界面变成受压损伤变量永远不会增长。排查方法是使用Abaqus的COORD选项输出单元坐标系或者在Model-Edit Attributes里查看单元方向。后处理里可以画出cohesive单元的S1法向应力、S2切向应力来直观判断看云图是否连续过渡。如果发现同一层cohesive单元里应力符号正负交替八九不离十是方向问题。另一个导致应力振荡的原因是软化段刚度过陡。比如你给了很高的断裂能但δ0和δf差距过小软化斜率非常陡单元刚度在几步内从K跳到接近零隐式迭代自然不稳定。这种情况要么微调Gc和t0让曲线更缓要么给损伤演化加粘性正则化。内聚力模型的好处是你可以通过调整Gc来控制软化斜率但注意Gc必须来自于真实的实验测量不能为了收敛而无限制放大。5.3 其他几个容易掉进去的坑我之前顺手整理过一个速查表直接分享出来问题现象可能原因解决办法SDV输出总是0分析步输出请求里没勾选SDV在Field Output中勾选SDV初始刚度比预期小10倍本构厚度T0设置不一致统一按T01处理并检查截面厚度损伤扩展方向不对cohesive单元厚度方向扫掠反了检查单元局部坐标1方向完全失效后单元刚度为0导致负特征值软化后刚度归零残余刚度设为初始刚度的1e-6倍收敛一步步缩小到1e-8还是不收敛模型里cohesive层被压溃而非拉伸检查边界条件是否造成局部压缩载荷-位移曲线震荡损伤起始强度过高、软化过陡降低t0或增大Gc或加粘性正则化还有一个容易被忽略的点单位制。前面所有参数我都是在mm、N、MPa、N/mm这套单位制下写的。如果你换成m、kg、s制断裂能Gc的单位是J/m²界面刚度K单位是N/m³数值会差得非常远。不少人把论文里的Gc直接抄进Abaqus结果差了三个数量级都不知道。5.4 网格细度和断裂过程区的关系内聚力模型对网格尺寸有一定敏感性但比纯应力-应变软化本构好太多。一个经验法则是断裂过程区的长度大致可以用这个公式估算L ≈ E × Gc / (t0)²其中E是界面附近的材料模量。如果过程区尺寸是0.2 mm你的网格尺寸就要明显小于这个值否则损伤带只在一个单元里发展结果表现得像脆断看不出渐进损伤。实际操作中我通常让过程区内至少有3到5个单元。这个估算不需要非常精确用来判断网格尺度的量级足够。如果模型太大、网格太细算不动怎么办优先考虑对称模型和子模型。做单搭接剪切时如果结构和载荷对称可以取半模型或者四分之一模型。当然要注意边界条件是否允许对称。内聚力模型本身计算量不大但软化阶段增量步很小模型小了收益非常明显。我个人在实际操作中的体会是UMAT的编写并不是最难的部分难的是理解cohesive单元的约定和调试时的耐心。每次看到负特征值先别急着改本构把单元坐标系、本构厚度、单位制这几样最基础的东西过一遍往往能找到罪魁祸首。我现在遇到界面开裂问题一定会先写一个单模态UMAT跑单单元确认曲线吻合理论再上完整模型这个习惯帮我省掉了至少一半的调试时间。如果你只是复现双线性模型内置cohesive behavior确实够用但后续要加温度、湿度、疲劳退化这类因素还是早点上手UMAT更省事。最后再分享一个小技巧在UMAT里多输出几个状态变量比如把当前牵引力、历史最大分离量、损伤起始标志分别存到不同的SDV里调试时对照云图就能一眼定位到某个积分点算到哪一步了这比只盯着D一个变量清晰得多。
