ABAQUS实战：cohesive单元与内聚力本构模型UMAT详解（附文件与教学视频）

在复合材料分层、胶接接头失效等界面问题模拟中，cohesive单元（内聚力单元）扮演着至关重要的角色。它能有效模拟界面在载荷作用下的逐渐脱粘过程，从初始的线弹性响应，到损伤萌生，直至最终完全失效。然而，ABAQUS内置的内聚力本构模型（Traction-Separation Law）虽然方便，但在面对一些特殊材料行为、复杂的损伤演化准则或需要与其它物理场耦合时，就显得力不从心了。这时，用户自定义材料子程序（UMAT）就成了解决问题的钥匙。

1. 为何需要自定义UMAT？内置模型不够用吗？

ABAQUS内置的牵引-分离法则通常是基于简单的最大应力或最大应变准则来初始化损伤，损伤演化则多采用基于能量或位移的线性/指数软化规律。对于大多数标准情况，这完全够用。但在以下场景，你就需要考虑自己动手写UMAT了：

• 材料行为特殊：你需要实现一种内置模型中没有的损伤起始准则（比如混合模式比判据更复杂）。
• 损伤演化复杂：损伤变量D的演化不遵循简单的线性或指数规律，可能和历程、速率或其他状态变量相关。
• 耦合分析需求：内聚力模型的参数需要随温度、湿度或其他场变量变化，实现完全耦合。
• 研究创新本构：你在开发或验证一种新的内聚力理论模型。

简单说，当内置的“黑箱”无法描述你材料界面的真实“脾气”时，UMAT让你拥有了自定义“脾气”规则的能力。

2. 内聚力模型UMAT的核心实现原理

内聚力模型的核心思想是界面应力随界面张开位移（分离）变化，并通过一个损伤变量D来表征刚度的退化。一个典型的线性软化本构关系可以表示为：

当损伤未起始时，处于线弹性阶段： $\sigma_i = K_0 \cdot \delta_i \quad (i = n, s, t)$ 其中，$\sigma$为牵引应力，$K_0$为初始刚度，$\delta$为分离位移。

损伤起始通常由一个混合模式判据决定，例如二次应力准则： $\left( \frac{\langle \sigma_n \rangle}{\sigma_n^0} \right)^2 + \left( \frac{\sigma_s}{\sigma_s^0} \right)^2 + \left( \frac{\sigma_t}{\sigma_t^0} \right)^2 = 1$ 这里$\langle \rangle$是麦考利括号，表示只考虑拉应力，$\sigma^0$是各方向的强度。

损伤开始后，引入损伤变量D (0 ≤ D ≤ 1)，应力计算变为： $\sigma_i = (1 - D) \cdot K_0 \cdot \delta_i$ D从0（无损）演化到1（完全失效）。D的演化往往基于一个等效位移$\delta_m$和断裂能$G_c$。例如，线性软化的演化公式为： $D = \frac{\delta_m^f (\delta_m - \delta_m^0)}{\delta_m (\delta_m^f - \delta_m^0)}$ 其中，$\delta_m^0$是损伤起始时的等效位移，$\delta_m^f$是完全失效时的等效位移，与断裂能相关。

在UMAT中，我们的任务就是根据传入的应变（在cohesive单元中对应分离位移）增量，依据上述理论，更新应力（STRESS）和材料雅可比矩阵（DDSDDE），并管理好状态变量（STATEV），比如损伤变量D和最大历史等效位移。

3. Fortran UMAT子程序代码详解（关键片段）

下面是一个简化的、基于最大应力准则和线性软化的cohesive单元UMAT框架，包含了关键部分的注释。

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS), 1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS), 2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3), 3 DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),JSTEP(4) C C MATERIAL PROPERTIES FROM PROPS ARRAY C PROPS(1) - K0, 初始刚度 C PROPS(2) - SIGMA_N0, 法向强度 C PROPS(3) - SIGMA_S0, 第一剪切强度 C PROPS(4) - SIGMA_T0, 第二剪切强度 C PROPS(5) - GIC, 模式I断裂能 C PROPS(6) - GIIC, 模式II断裂能 C STATEV(1) - DAMAGE, 损伤变量D C STATEV(2) - MAX_EQUIV_DISP, 历史最大等效位移 C REAL*8 K0, SIGMA_N0, SIGMA_S0, SIGMA_T0, GIC, GIIC REAL*8 DAMAGE, MAX_EQUIV_DISP, EQUIV_DISP, EQUIV_STRESS REAL*8 DELTA_N, DELTA_S, DELTA_T, SIGMA_N, SIGMA_S, SIGMA_T REAL*8 DELTA_N0, DELTA_S0, DELTA_T0, DELTA_M0, DELTA_MF REAL*8 DDSDDE_LOCAL(NTENS, NTENS) C C 初始化材料参数 K0 = PROPS(1) SIGMA_N0 = PROPS(2) SIGMA_S0 = PROPS(3) SIGMA_T0 = PROPS(4) GIC = PROPS(5) GIIC = PROPS(6) C C 获取上一步的状态变量 DAMAGE = STATEV(1) MAX_EQUIV_DISP = STATEV(2) C C 计算当前总分离位移（假设应变即位移，适用于零厚度单元） DELTA_N = STRAN(1) + DSTRAN(1) ! 法向 DELTA_S = STRAN(2) + DSTRAN(2) ! 第一剪切 DELTA_T = STRAN(3) + DSTRAN(3) ! 第二剪切 (NTENS=3 for cohesive) C C 计算当前等效位移和等效应力（基于二次应力准则的混合模式） C 注意：法向只考虑拉伸 EQUIV_STRESS = SQRT( (MAX(0.0, DELTA_N*K0)/SIGMA_N0)**2 + 1 (DELTA_S*K0/SIGMA_S0)**2 + 2 (DELTA_T*K0/SIGMA_T0)**2 ) EQUIV_DISP = SQRT( (MAX(0.0, DELTA_N))**2 + DELTA_S**2 + DELTA_T**2 ) C C *** 损伤起始判断 *** IF (EQUIV_STRESS .GE. 1.0 .AND. DAMAGE .LT. 1.0E-6) THEN DELTA_M0 = EQUIV_DISP ! 记录损伤起始位移 STATEV(3) = DELTA_M0 ! 假设STATEV(3)存储起始位移 END IF C C *** 损伤演化 *** IF (EQUIV_DISP .GT. MAX_EQUIV_DISP) THEN MAX_EQUIV_DISP = EQUIV_DISP END IF C IF (MAX_EQUIV_DISP .GT. STATEV(3)) THEN C 线性软化定律计算损伤变量 DELTA_M0 = STATEV(3) C 计算完全失效位移（基于断裂能，简化处理） DELTA_MF = 2.0 * (GIC + GIIC) / (K0 * DELTA_M0) IF (MAX_EQUIV_DISP .LT. DELTA_MF) THEN DAMAGE = (DELTA_MF * (MAX_EQUIV_DISP - DELTA_M0)) / 1 (MAX_EQUIV_DISP * (DELTA_MF - DELTA_M0)) ELSE DAMAGE = 1.0 END IF END IF C 确保损伤在0-1之间 DAMAGE = MIN(1.0, MAX(0.0, DAMAGE)) C C *** 更新应力 *** SIGMA_N = (1.0 - DAMAGE) * K0 * DELTA_N SIGMA_S = (1.0 - DAMAGE) * K0 * DELTA_S SIGMA_T = (1.0 - DAMAGE) * K0 * DELTA_T C 注意：法向为压力时，通常不考虑损伤（闭合） IF (DELTA_N .LT. 0.0) THEN SIGMA_N = K0 * DELTA_N END IF C STRESS(1) = SIGMA_N STRESS(2) = SIGMA_S STRESS(3) = SIGMA_T C C *** 计算材料雅可比矩阵 (DDSDDE) *** DO I=1, NTENS DO J=1, NTENS DDSDDE(I,J) = 0.0 END DO END DO C IF (DAMAGE .LT. 1.0E-6 .OR. DELTA_N .LT. 0.0) THEN C 无损或法向受压，弹性刚度 DDSDDE(1,1) = K0 DDSDDE(2,2) = K0 DDSDDE(3,3) = K0 ELSE C 损伤后，切线刚度矩阵（此处为连续线性化的简化形式） C 实际应根据损伤演化律求导得到精确切线，此处为示例。 DDSDDE(1,1) = (1.0 - DAMAGE) * K0 DDSDDE(2,2) = (1.0 - DAMAGE) * K0 DDSDDE(3,3) = (1.0 - DAMAGE) * K0 C 注意：对于复杂的演化律，需要考虑dD/dδ对非对角线项的影响 END IF C C *** 更新状态变量 *** STATEV(1) = DAMAGE STATEV(2) = MAX_EQUIV_DISP C RETURN END

4. ABAQUS中的建模流程与实例解析

有了UMAT，如何在ABAQUS中使用它呢？结合教学视频，关键步骤如下：

1. 单元类型选择：在Part或Mesh模块，为界面区域指定COH3D8（三维8节点）或COH2D4（二维4节点）等cohesive单元类型。
2. 材料属性定义：在Property模块，创建材料。在General->User Material中，输入我们在UMAT中定义的PROPS参数，例如初始刚度K0、各方向强度、断裂能等。将材料分配给cohesive单元截面。
3. 相互作用设置：如果cohesive单元是嵌入的或粘接的，确保其连接关系正确。对于零厚度cohesive单元，通常采用“Tie”约束或共享节点的方式与相邻实体单元连接。
4. 分析步与输出：在Step模块创建静力通用分析步。为了帮助收敛，可以打开自动稳定（Automatic Stabilization）或使用粘性正则化（Viscous Regularization）。在Field Output Request中，务必勾选SDV（状态变量）以输出我们定义的损伤变量D。
5. 提交作业：在Job模块创建作业时，在General选项卡下，将编译好的.for文件（或.obj文件）链接到User subroutine file。然后提交计算。

一个简单的.inp文件材料定义部分可能如下：

*Material, name=COHESIVE_UMA *User Material, constants=6 5.0E6, 50.0, 60.0, 60.0, 0.5, 1.0 ** *Depvar 2, **

其中，*User Material后的数字6表示有6个材料常数（对应PROPS数组），下面一行就是这些常数的值。*Depvar定义了状态变量的数量，这里是2个（损伤D和历史最大位移）。

5. 避坑指南与调试技巧

• 初始刚度K0取值过大：这是导致计算震荡、不收敛的最常见原因。K0需要足够大以保证界面在失效前的刚度，但又不能大到引起病态刚度矩阵。一个经验法则是取相邻实体材料等效模量的10到100倍。可以先取一个较小的值试算，再逐步增加。
• 损伤演化参数设置不当：断裂能Gc设置过小，会导致材料过于“脆”，失效过程剧烈，不易收敛。确保Gc值与实际材料测试结果相符。可以使用粘性正则化参数（在材料属性或分析步中设置）来软化失效过程，改善收敛。
• 监控损伤过程：在.inp文件中添加以下关键字，可以输出详细的单元状态信息到.msg文件，对于调试UMAT非常有帮助：

  *PRINT, ECHO=NO, MODEL=NO, HISTORY=NO, CONTACT=NO *DEBOG

  在*DEBOG后可以指定输出频率。计算后查看.msg文件，搜索你关心的单元号和变量。
• 状态变量初始化：确保在UMAT开始时正确读取STATEV，并在每一步结束时正确更新。不正确的状态变量传递是结果异常的主要原因之一。

6. 性能优化与网格敏感性

cohesive单元的模拟对网格尺寸较为敏感。因为断裂能Gc是在整个单元上耗散的。单元尺寸L_e会影响软化段的斜率。为了保证结果的可重复性，建议：

• 单元尺寸与最大损伤增量：在显式分析中，有一个经验关系：最大损伤增量 per increment ≈ (Gc / σ0) / (N * L_e)，其中N是单元厚度方向的积分点数。这有助于估算稳定时间增量。
• 进行网格敏感性分析：针对你的模型，用2-3种不同的界面区域网格尺寸进行计算，比较载荷-位移曲线和最终失效模式。当进一步细化网格对结果影响很小时，可以认为网格密度已足够。
• 一致化网格：cohesive单元最好与其粘接的实体单元在界面处具有协调的网格尺寸，避免因刚度突变导致应力奇异。

最后，一个开放的思考

我们目前实现的是一个准静态、率无关的内聚力模型。但在许多实际场景中，如冲击、疲劳载荷下，界面的失效行为与加载速率密切相关。如何将我们这个UMAT扩展至率相关损伤情况呢？

一个常见的思路是引入与应变率或分离速率相关的强度σ0(δ̇)和断裂能Gc(δ̇)。在UMAT中，我们需要通过TIME和DTIME变量估算当前的速率，然后根据速率修正损伤起始和演化参数。例如，可以定义一个类似于粘塑性框架的过应力函数，让损伤演化速率与当前“应力”和“率相关强度”的差值相关。这无疑会增加本构的复杂性和参数标定的难度，但对于模拟动态剥离、高速冲击等过程却是必要的进阶。这或许是你下一步可以深入探索的方向。

通过这次从理论到代码，再到建模落地的完整梳理，相信你对cohesive单元及其UMAT开发有了更扎实的理解。实践出真知，赶紧打开ABAQUS，用附带的文件和视频教程动手试一试吧！