企业官网建设流程全解析

先说个真实的场景：去年我做胶接接头失效分析，模型本身不复杂，两块铝合金板中间一层胶，但一到损伤软化阶段，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²（也就是MPa）；Separation是分离位移，量纲是mm。常规应力-应变关系里，应力应变相乘得到能量密度，而牵引-分离曲线下的面积直接就是断裂能，量纲是N/mm。

2. 双线性内聚力本构：三条曲线里的门道

2.1 弹性段、损伤起始和线性软化

最经典的内聚力本构是双线性模型，Abaqus内置的Cohesive Behavior默认就是这类。曲线分三段：上升段、软化段、失效段。上升段代表界面还没损伤，牵引力随分离量线性增加，斜率就是界面刚度K。到达损伤起始点后，曲线掉头向下，直到完全失效。

我画不出图，但可以用参数把这条曲线说得清楚。弹性段写成：

t = K × δ

其中K的单位是N/mm³，或者叫MPa/mm。这里有一个常见的问号：K和材料弹性模量E是什么关系？如果你把胶层看成厚度为T0的一串弹簧，那K就等于E除以T0：

K = 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

举个具体例子：设K=1e5 N/mm³，t0=20 MPa，Gc=0.5 N/mm，那么δ0=2e-4 mm，δf=0.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/(Gn+Gs)]^η

这里的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，初始0
• STATEV(2)：历史最大分离量δmax，初始0
• STATEV(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单元里它对应的量纲是界面刚度乘以本构厚度？这里容易绕。如果直接在本构厚度T0=1的情况下，界面刚度K的数值上就等于dσ/dε，所以上面代码里弹性段DDSDDE=KN，软化段为负斜率，逻辑是对的。一旦T0不等于1，DDSDDE必须乘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 mm，Abaqus会一上来就进入大变形软化阶段。建议分成几个分析步，或者用一个幅值曲线从0逐步加到0.06，这样能看到完整曲线。

我自己调试UMAT时一定会做这一步，并且把数值输出的每一行都跟理论值手算一次。弹性段第一个数据点如果t不等于K×δ，检查T0；软化段如果斜率不对，检查DDSDDE；失效段如果还有残余应力，那正常，数值上必须保留小刚度。

4.2 单搭接剪切模型的建模步骤

单单元验证通过后，就可以搭一个完整的单搭接剪切模型。几何尺寸建议参考经典胶接试样：上下两块铝板，长75 mm，宽25 mm，厚1.5 mm，搭接长度25 mm。中间胶层几何厚度0.1 mm，但cohesive单元的本构厚度仍然设成1，原因前面说过。

材料参数这样配：

• 铝板：线弹性，E=70 GPa，泊松比0.33。分析重点是界面失效，板的塑性可以暂时不加，避免结果解读复杂化。
• 胶层cohesive单元：Kn=1e5 N/mm³，Ks=1e5 N/mm³，tn0=20 MPa，ts0=20 MPa，Gn=0.5 N/mm，Gs=1.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里我写了T0=1，两者刚度自然差了10倍。这个坑说明一个道理：做对比之前，先确认两个模型对“本构厚度”的约定完全一致。

5. 常见问题与排查经验：这些坑我都替你踩过

5.1 收敛不了：负特征值、增量步一降再降

Cohesive单元做隐式分析，最大的敌人就是收敛问题。典型症状是消息区不停出现负特征值警告，增量步从0.1缩到1e-6，然后直接中止。负特征值并不代表模型物理上真的不稳定，很多时候是某个cohesive单元刚度退化成零，导致局部刚度矩阵奇异。

我处理这种问题优先级是这样的：

1. 在UMAT中对完全失效单元不返回0刚度，而是返回一个很小的残余刚度，比如初始刚度的1e-6倍。这能显著减少数值奇异。
2. 在材料属性里加粘性正则化，内置模型有visco选项，UMAT里可以自己给损伤变量D做粘性更新。Abaqus内置的粘性系数一般取0.001就能压住振荡。
3. 调整增量策略：初始增量调小，使用自动增量，打开非对称求解器。
4. 检查是不是网格畸变或者单元方向错乱。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 其他几个容易掉进去的坑

我之前顺手整理过一个速查表，直接分享出来：

还有一个容易被忽略的点：单位制。前面所有参数我都是在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一个变量清晰得多。