Anand粘塑性模型的UMAT子程序及验证

Anand 粘塑性模型的UMAT 子程序及验证

高军

1.引言

电子封装及其组件在工艺或者服役过程中, 由于功率耗散和环境温度的周期变化, 会因为电子印制电路板、芯片和焊点的热膨胀失配，在合金钎焊焊点处产生交变的应力应变, 导致焊点的电、热或者机械失效。焊点的热循环失效(可靠性)是电子封装及组装技术中的关键问题之一, 受到了人们的普遍关注。焊点体积细小, 应力应变很复杂。为了准确模拟焊点在服役条件下的应力应变响应, 对可靠性进行评估, 必须建立合理有效的描述钎焊合金材料力学响应的本构方程。

SnPb 基焊锡钎料广泛应用于电子封装领域，作为电的连接和机械的连接。对于钎料的力学性能的试验和本构模型，许多学者都进行了研究。通常SnPb 基焊锡钎料具有很强的温度和加载速率的相关性，应该采用统一型粘塑性本构模型描述SnPb 钎料的变形行为。

在统一型粘塑性本构模型中，应用最广泛的是Anand 模型。具有形式简单，模型参数少等特点，在电子焊点的寿命预测中广泛应用。它采用与位错密度、固溶体强化以及晶粒尺寸效应等相关的单一内部变量S 描述材料内部状态对塑性流动的宏观阻抗，可以反映粘塑性材料与应变速度、温度相关的变形行为，以及应变率的历史效应、应变硬化和动态回复等特征。

目前，很多大型商用有限元软件，如ANSYS 、MARC 等都把Anand 本构模型嵌入到通用材料模型库中供用户使用，但是，ABAQUS 的通用材料模型库中缺少Anand 模型。因此，本报告目的在于通过ABAQUS 的用户子程序接口UMAT ，选择合适的算法，将Anand 粘塑性本构模型引入ABAQUS 中，以便后续的研究。

2.Anand 本构方程

统一型粘塑性Anand 本构模型有两个基本特征：(1) 在应力空间没有明确的屈服面, 故在变形过程中不需要加载/卸载准则, 塑性变形在所有非零应力条件下产生。(2) 采用单一内部变量描述材料内部状态对塑性流动的宏观阻抗。内部变量(或称变形阻抗) 用S 标记, 具有应力量纲。

粘塑性Anand 模型的流动方程采用双曲蠕变规律对材料的率相关性与温度相关性进行预测，如下式：

n

T R Q e A p p S S p p S S sign S S h S s s s

m S s s T R Q Ae p T T I p C ⎥⎥

⎥

⎥⎦

⎤⎢⎢⎢⎢⎣⎡∧=⋅⎪⎭

⎪⎬⎫⎪⎩⎪⎨⎧--=⋅

⎥

⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎣⎡-=---= εεεεαξεαεεσ:32*:32)*1(*10:23

1):23sinh()(23)]

([:

式中，σ为Cauchy 应力，s 为偏应力, C 为弹性张量，ε为总应变, α为热膨胀系数，

p ε

为非弹性应变速率，A 为常数，Q 为激活能，m 为应变敏感指数，ξ为应力乘子，R

为气体常数， T 为绝对温度，T 0为参考温度，h 0为形变硬化-软化常数，a 为与硬化-软化相关的应变敏感系数，S *为变量饱和值，

∧

S 为系数，n 为指数。

粘塑性Anand 本构方程中，共有9个材料参数:A, Q, ξ,m, n, h 0，∧

S ,a 以及初始形变阻抗

S 0。

为真实模拟钎焊材料内部损伤变化，引入损伤，演化率如下式：

)]0

([::)]0([T T I p C T T I p D ------=αεεαεε 加入损伤的Anand 模型方程如下：

)

5(:32*)

4(:32

)*1(*10)

3(:2

31):2

3

sinh()(23

)

2()]0

([::)]0

([)1()]

([:)1(n

T R Q e A p p S S p p

S S sign S S h S s s s m S s s T R Q Ae p T T I p

C T T I p

D T T I p C D ⎥⎥

⎥

⎥⎦

⎤⎢⎢⎢⎢⎣⎡∧=⋅⎪⎭

⎪⎬⎫⎪⎩⎪⎨⎧--=⋅⎥⎥

⎥⎥⎦

⎤⎢⎢

⎢⎢⎣

⎡

-=------=---⋅-= εεεεαξε

αεεαεεαεεσ

3.算法与计算流程

计算（2），（3），（4）式，主要有三种数值算法，向前显式Euler 方法，向后隐式Euler 方法，中点法。向前显式Euler 方法是条件稳定，具有一阶精度；向后隐式Euler 方法和中点法是绝对稳定，分别具有一阶和二阶精度，它们均为隐式方法，需要利用迭代解隐式方程。迭代主要有四种方法：普通迭代，牛顿法，弦位法，抛物线法。

本报告中主要采用数值绝对稳定的向后Euler 方法和中点法两种数值算法，迭代采用普通迭代和弦位法，进行试算比较方法的优劣。

3.1 向后隐式Euler 方法+普通迭代

向后隐式Euler 计算公式为：

)1

,1(10)0(+++=+=n y n x hf n y n y y x y 向后Euler 方法是隐式方法，计算1

+n y

时要解隐式方程，通常用到迭代法。例如，先用向

前显式Euler 方法的计算结果作为初值，再作迭代，计算格式为：

))(1,1()1(1),()

0(1k n y n x hf n y k n y n y n x hf n y n y

+++=+++=+ 普通迭代的格式为：)(1k x k x

ϕ=+，判别迭代过程收敛的条件为：ε≤+-++)

(1

)1(1k n x k n x

采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：

）

（）

（）

（）

（）（10)1())()1((:))()1((32)1(*9))()1((:))()1((32))1(*)1(1()1(*)1(10)()1(8)1(:)1(23)1(1

))1()1(:)1(23sinh())1((23)()1(7)]0

)1(()1()1([::)]0)1(()1()1([)()1(6)]0

)1(()1()1([:))

1(1()1(n n T

R Q e tA n p n p n p n p S n S n p n p n p n p n S n S sign n S n S h n S n S n s n s n s m n S n s n s n T R Q Ae t n p n p T n T I n p n C T n T I n p n t n D n D T n T I n p

n C n D

n ⎥⎥

⎥⎥⎦

⎤⎢⎢⎢⎢⎣⎡+∆-+-+∧=+-+-+⋅⎪⎭

⎪⎬⎫

⎪⎩⎪⎨⎧++-++-+=++++⋅

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++++-⋅∆+=+-+-+-+-+-+-+⋅∆+=+-+-+-++-=+ εεεεεεεεαξεεαεεαεεαεεσ

根据上述计算格式，UMA T 子程序的计算流程为：

（1） 读取由ABAQUS 传递给UMAT 子程序的)1(,)(,)()(,,)(,)(+∆n T n S n p

n D n n εεεσ， 和T ∆，作为计算初值；

（2）采用迭代法，联立方程（6）-（10）式，求解)1()1()1(,)1(++++n S n p

n D n 及，εσ； （3）更新应力及全部状态变量，更新Jacobian 矩阵。

其中，迭代法的计算流程具体如下：

（1）迭代循环开始，针对)1()1()1(,)1(++++n S n p

n D n 及，εσ赋计算初值； （2）将方程（6）-（9）式写成形如),(1k k k y x f y =+的迭代格式，由第k 步的

)1()1(,)1(+++n p n D

n εσ，及)1(+n S 计算第k+1步的)1()1(,)1(+++n p

n D n εσ，及)1(+n S

；

（3）分析比较第k 步与第k+1步的)1()1(,)1(+++n p

n D n εσ，及)1(+n S ，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，

迭代失败。

向后Euler 方法具有绝对数值稳定性，误差具有一阶精度。虽然是绝对稳定的，但是迭代步长仍要受到一定限制。