把寿命写进混合键合(一)：用 UMAT 把“损伤”和“寿命”写进材料
用 UMAT 把“损伤”和“寿命”写进材料（ΔW → D）
本篇目标：
用一个最小可行的疲劳损伤 UMAT 模型，把混合键合 Cu–Cu 区域从“只会出应力云图”，升级为“会记忆历史损伤的材料”，并能在 Abaqus 里直接画出 寿命云图（D 云图）
目录
为什么用 ΔW 驱动损伤？
最小损伤模型：两个状态变量就够
UMAT 里的数据流长什么样？
关键 UMAT 实现（Fortran 核心代码）
用 UVARM 把 Wacc / D 写进云图
数值稳定性与参数选取小贴士
如何对接混合键合 Cu–Cu pad 场景
小结与下篇预告
1. 为什么用 ΔW 驱动损伤？
在混合键合（尤其 Cu–Cu 界面附近），我们真正关心的不是：
“这一瞬间等效应力有多红？”
而是：
“在反复热循环 / 功率循环之后，这个区域累积了多少疲劳损伤，离失效还有多远？”
常见的疲劳驱动量有三类：
应力幅： Δσ应变幅： Δε（或者塑性应变幅Δε_p）应变能密度幅： ΔW
这里我们选 ΔW 做驱动量，有几个好处：
ΔW同时包含 应力 × 应变 的信息，对高应力 + 大变形的危险区更敏感。自然适合累积：
累积应变能密度：
W_acc(n+1) = W_acc(n) + max(ΔW, 0)在多轴、非比例加载、热-力耦合下，比只用单一的
Δε_p更好扩展。
所以在本篇，我们用一个非常“工程化”的设计：
用
ΔW → W_acc → D的顺序，把“累积损伤”写进 UMAT。
2. 最小损伤模型：两个状态变量就够
我们把 Cu–Cu 混合键合区域视作一种“等效疲劳损伤材料”，定义一个标量损伤变量 D：
2.1 应力与刚度退化
未损伤弹性关系：
σ_el = C : ε
其中：
σ_el：未损伤弹性应力；C：各向同性弹性刚度矩阵；ε：总应变。
考虑损伤后的应力与刚度：
σ = (1 - D) * σ_el
C_eff = (1 - D) * C
D = 0：完好，刚度不变；D → 1：完全失效，刚度趋近 0（数值上我们不会真的到 1，稍后说）。
2.2 状态变量设计（STATEV）
只用两个状态变量就能跑起来：
STATEV(1) = W_acc 累积应变能密度
STATEV(2) = D 损伤变量（0~1）
增量内：
根据应力和应变增量，计算能量增量 ΔW；累积到 W_acc；根据 W_acc更新损伤D；用 D去衰减刚度与应力。
2.3 ΔW 的计算方式
在某个时间增量内：
旧应力： σ_old这一增量算出的未损伤弹性应力： σ_el应变增量： Δε
采用一个常见的近似（能量平均）：
ΔW ≈ 0.5 * (σ_old + σ_el) : Δε
若出现负值（卸载或数值噪声），我们简单处理：
if ΔW < 0 → 取 0
避免把“卸载”当作“损伤愈合”。
2.4 W_acc 与 D 的关系
累积能量：
W_acc_new = W_acc_old + max(ΔW, 0)
简单但好用的损伤演化：
D = (W_acc / W_f)^α
W_f：临界应变能密度（等效“损伤满格”的能量）；α：形状参数（控制曲线形状）；
数值上需要做截断：
D_max = 0.99
如果 D > D_max → 取 D_max
如果 D < 0 → 取 0
这一套关系就是本篇的“最小疲劳损伤模型”。 核心只有三行：
ΔW、W_acc、D，但已经可以让材料“记住被虐了多少次”。
3. UMAT 里的数据流长什么样？
把上面的物理关系塞进 UMAT，可以用一个“九步法”来理解：
读取材料参数
PROPSPROPS(1) = EPROPS(2) = νPROPS(3) = W_fPROPS(4) = α读取旧状态变量
STATEVSTATEV(1) → W_acc_oldSTATEV(2) → D_old计算当前总应变
εEPS(i) = STRAN(i) + DSTRAN(i)构造各向同性弹性矩阵
C用 E、ν算λ和μ，填 6×6 矩阵。未损伤弹性应力
σ_el = C : ε计算能量增量
ΔW并累积到W_acc_new根据
W_acc_new更新D_new用
D_new衰减应力与刚度σ = (1 - D_new) * σ_el
C_eff = (1 - D_new) * C回写
STATEV(1:2)
这就是 UMAT 里面数据的完整闭环。
4. 关键 UMAT 实现（Fortran 核心代码）
下面是一段可以直接编译的 UMAT 核心示例（3D 小应变各向同性 + ΔW 损伤）。 如果后续要加弹塑 / 蠕变，可以在 SIG_EL 的计算部分替换为“弹塑积分”。
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,
& DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
& TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
& NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,
& DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
& NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
INCLUDE'aba_param.inc'
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),
& STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(*),DPRED(*),
& PROPS(NPROPS),COORDS(*),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
& DDSDDT(NTENS,NTENS),DRPLDE(NTENS)
INTEGER NDI,NSHR,NTENS,NSTATV,NPROPS,NOEL,NPT,LAYER,
& KSPT,KSTEP,KINC
DOUBLEPRECISION E,NU,Wf,ALPHA
DOUBLEPRECISION C(6,6),EPS(6),SIG_EL(6),STRESS_OLD(6)
DOUBLEPRECISION WACC_OLD,WACC_NEW,D_OLD,D_NEW,Dmax
DOUBLEPRECISION dW,SSE_loc
DOUBLEPRECISION ZERO,ONE,LAMBDA,MU
INTEGER I,J
ZERO = 0.0D0
ONE = 1.0D0
C 材料参数：E, nu, Wf, alpha
E = PROPS(1)
NU = PROPS(2)
Wf = PROPS(3)
ALPHA = PROPS(4)
C 读取旧状态变量
WACC_OLD = STATEV(1)
D_OLD = STATEV(2)
C 旧应力备份
DO I = 1,NTENS
STRESS_OLD(I) = STRESS(I)
ENDDO
C 当前总应变 = 旧应变 + 增量
DO I = 1,NTENS
EPS(I) = STRAN(I) + DSTRAN(I)
ENDDO
C 构造各向同性弹性矩阵 C
DO I = 1,6
DO J = 1,6
C(I,J) = ZERO
ENDDO
ENDDO
LAMBDA = E*NU / ((ONE+NU)*(ONE-2.D0*NU))
MU = E / (2.D0*(ONE+NU))
C(1,1) = LAMBDA + 2.D0*MU
C(2,2) = C(1,1)
C(3,3) = C(1,1)
C(1,2) = LAMBDA
C(1,3) = LAMBDA
C(2,3) = LAMBDA
C(2,1) = LAMBDA
C(3,1) = LAMBDA
C(3,2) = LAMBDA
C(4,4) = MU
C(5,5) = MU
C(6,6) = MU
C 未损伤弹性应力 SIG_EL = C : EPS
DO I = 1,NTENS
SIG_EL(I) = 0.0D0
DO J = 1,NTENS
SIG_EL(I) = SIG_EL(I) + C(I,J)*EPS(J)
ENDDO
ENDDO
C ΔW ≈ 0.5 * (σ_old + σ_el) : Δε
dW = 0.0D0
DO I = 1,NTENS
dW = dW + 0.5D0 * (STRESS_OLD(I) + SIG_EL(I))*DSTRAN(I)
ENDDO
IF (dW .LT. 0.0D0) dW = 0.0D0
WACC_NEW = WACC_OLD + dW
C 更新损伤 D
IF (Wf .GT. 0.0D0) THEN
D_NEW = (WACC_NEW / Wf)**ALPHA
ELSE
D_NEW = D_OLD
ENDIF
Dmax = 0.99D0
IF (D_NEW .GT. Dmax) D_NEW = Dmax
IF (D_NEW .LT. 0.0D0) D_NEW = 0.0D0
C 应力与切线刚度（忽略 dD/dε 项）
DO I = 1,NTENS
STRESS(I) = (ONE - D_NEW) * SIG_EL(I)
ENDDO
DO I = 1,NTENS
DO J = 1,NTENS
DDSDDE(I,J) = (ONE - D_NEW) * C(I,J)
ENDDO
ENDDO
C 回写状态变量
STATEV(1) = WACC_NEW
STATEV(2) = D_NEW
C 内能（可选）
SSE_loc = 0.0D0
DO I = 1,NTENS
SSE_loc = SSE_loc + 0.5D0 * EPS(I)*STRESS(I)
ENDDO
SSE = SSE_loc
SPD = 0.0D0
SCD = 0.0D0
RPL = 0.0D0
DO I = 1,NTENS
DRPLDE(I) = 0.0D0
DO J = 1,NTENS
DDSDDT(I,J) = 0.0D0
ENDDO
ENDDO
DRPLDT = 0.0D0
PNEWDT = ONE
RETURN
END
5. 用 UVARM 把 Wacc / D 写进云图
STATEV 本身在 Viewer 里看不到，所以我们通常再配一个 UVARM，把 W_acc 和 D 映射成用户输出变量：
5.1 UVARM 子程序
SUBROUTINE UVARM(UVAR,NUVARM,ARR,NRR,NOEL,NPT,LAYER,KSPT,
& KSTEP,KINC,TIME,DTIME,CMNAME,ORNAME,
& NFIELDV,FIELD,NSTATV,STATEV,NUMPROPS,PROPS,
& COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
& JMAC,JMATYP,MATLAYO,LACCFLA)
INCLUDE'aba_param.inc'
CHARACTER*80 CMNAME,ORNAME
DIMENSION UVAR(NUVARM),ARR(NRR),TIME(2),FIELD(NFIELDV),
& STATEV(NSTATV),PROPS(NUMPROPS),COORDS(*),
& DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
INTEGER JMAC(*),JMATYP(*)
C 约定：STATEV(1) = W_acc, STATEV(2) = D
IF (NUVARM .GE. 1) UVAR(1) = STATEV(1)
IF (NUVARM .GE. 2) UVAR(2) = STATEV(2)
RETURN
END
5.2 输入文件中的声明
材料卡中：
*Material, name=CU_CU_DMG
*User Material, constants=4
** E, NU, Wf, ALPHA
110000., 0.34, 5.0E5, 1.0
*Depvar
2
*User Output Variables, variables=2
输出请求中：
*Output, field, variable=PRESELECT
*Output, field
*Element Output, elset=PAD
S, UVARM
这样在 Viewer 里：
UVARM1= 累积应变能密度W_acc；UVARM2= 损伤变量D（可以视为“寿命消耗比例”）。
6. 数值稳定性与参数选取小贴士
6.1 如何估计 W_f 的量级？
简单工程做法：
先用纯弹性材料跑一个基准工况：例如一次热循环；
在高危区域记录
ΔW*的数量级；目标设计寿命为
N_design次循环，则：W_f ≈ N_design * ΔW*
这是一个粗估，后续可以用实验 / 文献数据微调。
6.2 α 的影响
α = 1：D 与W_acc线性关系，最简单；α > 1：前期损伤慢、后期突然加速，更接近很多疲劳现象；一般不建议 α < 1，初始就损伤太快。
建议：初版先用 α = 1，后续根据疲劳曲线斜率再调。
6.3 D_max 为何不取 1？
若 D = 1，刚度严格为 0，刚度矩阵会非常病态；取 D_max = 0.95 ~ 0.99，既保持“几乎失效”，又不至于数值崩掉；真正“失效判断”可以在后处理里用阈值：例如 D ≥ 0.8视作失效。
7. 如何对接混合键合 Cu–Cu pad 场景
在混合键合仿真里，这个 UMAT 的用法可以非常简单粗暴：
在模型中为 Cu–Cu 混合键合区域定义一个 elset，例如：
*Elset, elset=PAD
...只在
PAD区域使用CU_CU_DMG材料：*Solid Section, elset=PAD, material=CU_CU_DMG其余区域（Si、oxide、underfill 等）仍用常规弹性 / 弹塑材料。
施加热循环、功率循环、机械载荷之后，直接看：
可以当作“寿命消耗曲线”。 哪些 pad 的损伤更高； pad 中心 vs 边缘的差异； UVARM2（D）的空间分布：D 随 step / 循环数的演化：
基于这套 D 云图，我们可以做很多设计类分析：
调整 pad 尺寸 / pitch，看最大 D 是否下降； 调整堆叠结构或材料（例如加减中介层），观察 D 云图变化； 结合功率分布 / 热设计，找出“可靠性短板 pad”。
8. 小结与下篇预告
本篇，我们完成了几件关键的事：
选择
ΔW作为疲劳驱动量，用一个简单的关系：ΔW → W_acc → D设计了只用两个状态变量（
W_acc和D）的最小损伤模型；给出了一段可直接编译的 UMAT 核心代码；
用 UVARM 把
W_acc / D写进云图，为后续“寿命云图 + 统计分析”打好地基；讨论了
W_f、α、D_max的选择建议和数值稳定性问题。
到这里，我们已经可以在一个简化立方体 / pad 模型上测试这套 UMAT， 看看在热循环下，“材料会不会自己记账”。
下一篇我们要做什么？
在第三篇里，我们会：
搭一个最小的混合键合 Cu–Cu pad 结构模型（上、下 Si + 中间一层 Cu–Cu）； 给出完整可运行的 .inp（几何、材料、载荷全部配好）；挂上本篇的 UMAT + UVARM，跑一个简化热循环； 展示真实的 D 云图、W_acc 云图，以及如何用 Python 做简单统计。
如果把这段 UMAT 编译通过、在一个小 block 模型上跑起来，下一篇就可以直接平移到混合键合结构上，基本是“复 制粘贴级别”的工作量了。
